Dynamical Friction and the Evolution of Supermassive Black Hole Binaries: The Final Hundred-parsec Problem

The general reduced two-body problem where Newtonian gravity is the only force acting on the two bodies in the system is described by the equation

$$\begin{eqnarray}&&\ddot{{\boldsymbol{r}}}=-\displaystyle \frac{{GM}}{{r}^{3}}{\boldsymbol{r}},\end{eqnarray} \tag{ 2 }$$

where G is the gravitational constant, M is the total mass of the system, and r is the relative position between the two bodies.

Any dynamical interaction between the two bodies introduces an extra force to the binary which acts as a perturbation to the classic two-body problem. Under the effect of a perturbing force ${\boldsymbol{F}}({\boldsymbol{r}},\dot{{\boldsymbol{r}}})$, the equation of motion for the perturbed two-body problem can be written as

$$\begin{eqnarray}&&\ddot{{\boldsymbol{r}}}=-\displaystyle \frac{{GM}}{{r}^{3}}{\boldsymbol{r}}+{\boldsymbol{F}}({\boldsymbol{r}},\dot{{\boldsymbol{r}}}),\end{eqnarray} \tag{ 3 }$$

where the perturbing force ${\boldsymbol{F}}({\boldsymbol{r}},\dot{{\boldsymbol{r}}})$ depends in principle upon both the relative position ${\boldsymbol{r}}$ and velocity ${\boldsymbol{v}}=\dot{{\boldsymbol{r}}}$.

Equation (3) can be solved assuming that at each instant of time, the true orbit can be approximated by an instantaneous ellipse that is changing over time through its now time-dependent orbital elements C_i(t). Here ${C}_{i}=(a,e,i,{\rm{\Omega }},\omega ,f)$ are namely the semimajor axis, eccentricity, inclination, longitude of the ascending node, argument of periapsis, and true anomaly f. We also define $n={({GM}/{a}^{3})}^{1/2}$ as the mean motion. At each moment of time, these orbital elements describe the orbit the body would follow if perturbations were to cease instantaneously. We refer to these elements as osculating orbital elements. The time-evolution equations for the osculating orbital elements decomposed in the reference system ${K}_{R}(\hat{{\boldsymbol{r}}},\hat{{\boldsymbol{\tau }}},\hat{{\boldsymbol{n}}})$, where the unit vector $\hat{{\boldsymbol{r}}}$ is along the relative position vector between the two bodies in the binary, become

$$\begin{eqnarray}&&\displaystyle \frac{{da}}{{dt}}=\displaystyle \frac{2}{n\sqrt{1-{e}^{2}}}[{F}_{r}e\sin f+{F}_{\tau }(1+e\cos f)],\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{de}}{{dt}}=\displaystyle \frac{\sqrt{1-{e}^{2}}}{{na}}\left[{F}_{r}\sin f+{F}_{\tau }\left(\cos f+\displaystyle \frac{e+\cos f}{1+e\cos f}\right)\right],\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{di}}{{dt}}={F}_{n}\displaystyle \frac{\cos (f+\omega )\sqrt{1-{e}^{2}}}{{na}(1+e\cos f)},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{\rm{\Omega }}}{{dt}}=\displaystyle \frac{1}{{{na}}^{2}\sqrt{1-{e}^{2}}\sin i}{F}_{n}\displaystyle \frac{a(1-{e}^{2})\sin (f+\omega )}{1+e\cos f},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}\displaystyle \frac{d\omega }{{dt}} & = & \displaystyle \frac{\sqrt{1-{e}^{2}}}{{nae}}\left(-{F}_{r}\cos f+{F}_{\tau }\sin f\displaystyle \frac{2+e\cos f}{1+e\cos f}\right)\\ & & -\displaystyle \frac{\cos i}{{na}\sin i}{F}_{n}\displaystyle \frac{\sqrt{1-{e}^{2}}\sin (f+\omega )}{1+e\cos f},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{df}}{{dt}}=\displaystyle \frac{n{(1+e\cos f)}^{2}}{{(1-{e}^{2})}^{3/2}}-\displaystyle \frac{d\omega }{{dt}}-\cos i\displaystyle \frac{d{\rm{\Omega }}}{{dt}}.\end{eqnarray} \tag{ 9 }$$

Equations (4) and (5) indicate that for a perturbing force always vertical to the orbital plane, i.e., ${F}_{r}={F}_{\tau }=0$, the semimajor axis and the eccentricity do not change while Equations (6) and (7) show that the inclination i and the longitude of the ascending node Ω evolve only for a non-zero vertical to the orbital plane component of the perturbing force, i.e., ${F}_{n}\ne 0$. On the contrary, from Equation (8) we see that the periapsis is precessing for any non-zero perturbing force, i.e., ${\boldsymbol{F}}\ne 0$.

In the following section, we apply the formalism described above to study the effect of dynamical friction on the orbital evolution of a test mass moving inside a cluster of stars around a central SMBH. In this case, dynamical friction acts as a perturbing force on the evolution of the inspiraling object. We begin with a brief introduction to dynamical friction as described initially by Chandrasekhar (1943) and further studied by Antonini & Merritt (2012) in the case of a test mass moving inside a cluster of stars around a more massive central SMBH.

In what follows we consider the evolution of a binary comprising a massive object moving near a SMBH of considerably larger mass and which sits at the center of a galaxy. Gravitational interaction with stars results in the loss of energy and angular momentum by the massive object. We assume that the galaxy is spherically symmetric and isotropic, and we describe it using a power-law stellar density profile of the form $\rho (r)={\rho }_{0}{(r/{r}_{0})}^{-\gamma }$, where r₀ is a characteristic radius and γ is the slope of the density profile. In what follows, we set ${r}_{0}={r}_{\mathrm{infl}}$, with ${r}_{\mathrm{infl}}$ the radius containing a mass in stars twice the mass of the central SMBH (${M}_{\bullet }$), i.e., the SMBH influences radius. Unless otherwise specified, in what follows we use units such that ${M}_{\bullet }={r}_{\mathrm{infl}}=G=1$.

Assuming that the gravitational potential Φ is dominated by the central SMBH and neglecting the effect of the surrounding stars, we can write ${\rm{\Phi }}\approx -{{GM}}_{\bullet }/r$. Eddington's formula then uniquely leads to the following distribution function of the field-star velocities (Merritt 2013, p. 544):

$$\begin{eqnarray}&&f({\upsilon }_{\star })=\displaystyle \frac{{\rm{\Gamma }}(\gamma +1)}{{\rm{\Gamma }}\left(\gamma -\tfrac{1}{2}\right)}\ \displaystyle \frac{1}{{2}^{\gamma }{\pi }^{3/2}{v}_{{\rm{c}}}^{2\gamma }}\ {(2{v}_{{\rm{c}}}^{2}-{\upsilon }_{\star }^{2})}^{\gamma -3/2},\end{eqnarray} \tag{ 10 }$$

where ${\upsilon }_{\star }$ is the star velocity, ${\upsilon }_{{\rm{c}}}=\sqrt{{{GM}}_{\bullet }/r}$ is the circular velocity, and the normalization corresponds to unit total number.

The general formula for the dynamical friction force also including the contribution from stars moving faster than the massive body is

$$\begin{eqnarray}{{\boldsymbol{F}}}_{\mathrm{df}} & \approx & -4\pi {G}^{2}m\rho (r)\displaystyle \frac{{\boldsymbol{\upsilon }}}{{\upsilon }^{3}}\times \left\{\mathrm{ln}\,{\rm{\Lambda }}{\displaystyle \int }_{0}^{\upsilon }d{\upsilon }_{\star }4\pi f({\upsilon }_{\star }){\upsilon }_{\star }^{2}\right.\\ & & \left.+{\displaystyle \int }_{\upsilon }^{{\upsilon }_{{esc}}}d{\upsilon }_{\star }4\pi f({\upsilon }_{\star }){\upsilon }_{\star }^{2}\left[\mathrm{ln}\left(\displaystyle \frac{{\upsilon }_{\star }+\upsilon }{{\upsilon }_{\star }-\upsilon }\right)-2\displaystyle \frac{\upsilon }{{\upsilon }_{\star }}\right]\right\},\end{eqnarray} \tag{ 11 }$$

where ${\boldsymbol{\upsilon }}$ is the velocity of the massive body and m its mass. The quantity $\mathrm{ln}\,{\rm{\Lambda }}$ is the Coulomb logarithm defined as

$$\begin{eqnarray}&&\mathrm{ln}{\rm{\Lambda }}=\mathrm{ln}\left(\displaystyle \frac{{b}_{\max }}{{b}_{\min }}\right)\approx \mathrm{ln}\left(\displaystyle \frac{{b}_{\max }{v}_{{\rm{c}}}^{2}}{{Gm}}\right),\end{eqnarray} \tag{ 12 }$$

where ${b}_{\max }$ and ${b}_{\min }$ are the maximum and minimum impact parameters, respectively.

We can rewrite the dynamical friction force in a more compact form as

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{\mathrm{df}}=\kappa (r)[\mathrm{ln}\,{\rm{\Lambda }}\ \ \alpha (v)+\beta (v)+\delta (v)]\displaystyle \frac{{\boldsymbol{v}}}{{v}^{3}},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{\mathrm{df}}=\epsilon (r,v)\displaystyle \frac{{\boldsymbol{v}}}{{v}^{3}},\end{eqnarray} \tag{ 14 }$$

where we defined

$$\begin{eqnarray}&&\kappa (r)=-4\pi {G}^{2}\rho (r)m,\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\epsilon (r,v)=\kappa (r)[\mathrm{ln}\,{\rm{\Lambda }}\ \alpha (v)+\beta (v)+\delta (v)],\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\alpha (v)=4\pi {\int }_{0}^{v}f({\upsilon }_{\star }){\upsilon }_{\star }^{2}d{\upsilon }_{\star },\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\beta (v)=4\pi {\int }_{v}^{{v}_{\mathrm{esc}}}f({\upsilon }_{\star }){\upsilon }_{\star }^{2}\left[\mathrm{ln}\left(\displaystyle \frac{{\upsilon }_{\star }+\upsilon }{{\upsilon }_{\star }-\upsilon }\right)\right]d{\upsilon }_{\star },\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\delta (v)=4\pi v{\int }_{v}^{{v}_{\mathrm{esc}}}f({\upsilon }_{\star })(-2{\upsilon }_{\star })d{\upsilon }_{\star }.\end{eqnarray} \tag{ 19 }$$

The distribution function (10) can be rewritten as

$$\begin{eqnarray}&&f({\upsilon }_{\star })=\displaystyle \frac{{\rm{\Gamma }}(\gamma +1)}{{\rm{\Gamma }}\left(\gamma -\tfrac{1}{2}\right)}\ \displaystyle \frac{1}{{2}^{\gamma }{\pi }^{3/2}{v}_{{\rm{c}}}^{3}}{[2-{x}^{2}]}^{b}\end{eqnarray} \tag{ 20 }$$

where we defined $x={\upsilon }_{\star }/{v}_{{\rm{c}}}$ and $b=\gamma -3/2$. Integrals (17) and (19) have an analytic form while integral (18) demands numerical manipulation. Using Equation (20), we can rewrite the above integrals in dimensionless form as

$$\begin{eqnarray}\alpha (\xi ) & = & \displaystyle \frac{{\rm{\Gamma }}(\gamma +1)}{{\rm{\Gamma }}\left(\gamma -\tfrac{1}{2}\right)}\ \displaystyle \frac{4}{3}{\pi }^{-1/2}{2}^{b-\gamma }{\xi }^{3}\\ & & \times \,{}_{2}{F}_{1}[3/2,-b,5/2,{\xi }^{2}/2],\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&\beta (\xi )=\displaystyle \frac{{\rm{\Gamma }}(\gamma +1)}{{\rm{\Gamma }}\left(\gamma -\tfrac{1}{2}\right)}4{\pi }^{-1/2}{2}^{-\gamma },\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&\times {\int }_{\xi }^{1.4}{x}^{2}{(2-{x}^{2})}^{b}\mathrm{ln}\left(\displaystyle \frac{x+\xi }{x-\xi }\right){dx},\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&\delta (\xi )=\displaystyle \frac{{\rm{\Gamma }}(\gamma +1)}{{\rm{\Gamma }}\left(\gamma -\tfrac{1}{2}\right)}8{\pi }^{-1/2}\displaystyle \frac{{2}^{-\gamma -1}}{b+1}\xi ,\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&\times [{0.04}^{b+1}-{(2-{\xi }^{2})}^{b+1}]\end{eqnarray} \tag{ 25 }$$

where $\xi =v/{v}_{{\rm{c}}}$, and we made use of that fact that ${v}_{\mathrm{esc}}/{v}_{{\rm{c}}}=\sqrt{2}$.

The orbital velocity and dynamical friction force (14) decomposed in the reference system ${K}_{R}(\hat{{\boldsymbol{r}}},\hat{{\boldsymbol{\tau }}},\hat{{\boldsymbol{n}}})$ mentioned above can be written as

$$\begin{eqnarray}&&{\boldsymbol{\upsilon }}={v}_{r}\hat{{\boldsymbol{r}}}+{v}_{\tau }\hat{{\boldsymbol{\tau }}}+{v}_{n}\hat{{\boldsymbol{z}}},\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{\mathrm{df}}={F}_{\mathrm{df},r}\hat{{\boldsymbol{r}}}+{F}_{\mathrm{df},\tau }\hat{{\boldsymbol{\tau }}}+{F}_{\mathrm{df},n}\hat{{\boldsymbol{z}}}.\end{eqnarray} \tag{ 27 }$$

Substituting Equation (26) into Equation (14), we have ${F}_{\mathrm{df},r}=\epsilon (r,v)\tfrac{{v}_{r}}{{v}^{3}}$, ${F}_{\mathrm{df},\tau }=\epsilon (r,v)\tfrac{{v}_{\tau }}{{v}^{3}}$ and ${F}_{\mathrm{df},n}={v}_{n}=0$.

Using the dynamical friction components derived above and Equations (4)–(9), the osculating orbital element time-evolution equations of an inspiraling massive body due to dynamical friction are

$$\begin{eqnarray}&&\displaystyle \frac{{da}}{{dt}}=\displaystyle \frac{2\epsilon (r,v)}{{n}^{3}{a}^{2}}\displaystyle \frac{{(1-{e}^{2})}^{1/2}}{{(1+{e}^{2}+2e\cos f)}^{1/2}},\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{de}}{{dt}}=\displaystyle \frac{2\epsilon (r,v)}{{n}^{3}{a}^{3}}{(1-{e}^{2})}^{3/2}\displaystyle \frac{e+\cos f}{{(1+{e}^{2}+2e\cos f)}^{3/2}},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{di}}{{dt}}={F}_{\mathrm{df},n}\displaystyle \frac{\cos (f+\omega )\sqrt{1-{e}^{2}}}{{na}(1+e\cos f)}=0,\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{\rm{\Omega }}}{{dt}}=\displaystyle \frac{1}{{{na}}^{2}\sqrt{1-{e}^{2}}\sin i}{F}_{\mathrm{df},n}\displaystyle \frac{\cos (f+\omega )\sqrt{1-{e}^{2}}}{{na}(1+e\cos f)}=0,\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\omega }{{dt}}=\displaystyle \frac{2\epsilon (r,v)}{{n}^{3}{a}^{3}}\displaystyle \frac{{(1-{e}^{2})}^{3/2}}{{(1+{e}^{2}+2e\cos f)}^{3/2}}\displaystyle \frac{\sin f}{e},\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{df}}{{dt}}=\displaystyle \frac{n{(1+e\cos f)}^{2}}{{(1-{e}^{2})}^{3/2}}-\displaystyle \frac{d\omega }{{dt}}.\end{eqnarray} \tag{ 33 }$$

Equations (28)–(33) verify the aforementioned comment that in the absence of a vertical to the orbital plane component of the perturbing force, which is true in the case of dynamical friction (${F}_{\mathrm{df},n}=0$), the inclination i and the longitude of the ascending node Ω remain constant in the absence of other perturbing forces. Although according to Equation (32) dynamical friction can in principle induce a precession to the orbit, this precession has a negligible effect on the evolution of e and a.

We focus here on how dynamical friction affects the binary eccentricity. We begin by investigating qualitatively the expected eccentricity evolution using a simplified physical picture of the problem. This picture focuses on the eccentricity changes near periapsis and apoapsis. This analysis is useful to understand the link between the expected eccentricity evolution of the system and the physical origin of the dynamical friction force. The time evolution of the eccentricity vector ${\boldsymbol{e}}$ induced by a perturbing force ${\boldsymbol{F}}$ (in our case ${\boldsymbol{F}}={{\boldsymbol{F}}}_{\mathrm{df}}$) is given by

$$\begin{eqnarray}&&\dot{{\boldsymbol{e}}}=\displaystyle \frac{1}{{GM}}({{\boldsymbol{F}}}_{\mathrm{df}}\times {\boldsymbol{h}}+{\boldsymbol{\upsilon }}\times \dot{{\boldsymbol{h}}}),\end{eqnarray} \tag{ 34 }$$

where ${\boldsymbol{h}}={\boldsymbol{r}}\times {\boldsymbol{\upsilon }}$ is the angular momentum per unit reduced mass $\mu ={M}_{\bullet }m/({M}_{\bullet }+m)$, and the dot indicates a derivative with respect to time. In the absence of other perturbing forces in the binary, the angular momentum changes are only due to dynamical friction.

The dynamical friction force given by Equation (11) is a decelerating drag force (i.e., ${{\boldsymbol{F}}}_{\mathrm{df}}=-g(\upsilon ,r){\boldsymbol{\upsilon }}$, where $g(\upsilon ,r)$ is a function of the massive body position and velocity) always acting in the direction opposite to the direction of motion of the body. In addition, the time evolution of the angular momentum vector is given by $\dot{{\boldsymbol{h}}}={\boldsymbol{r}}\times {{\boldsymbol{F}}}_{\mathrm{df}}\Rightarrow \dot{{\boldsymbol{h}}}=-g(\upsilon ,r){\boldsymbol{r}}\times {\boldsymbol{\upsilon }}$. This leads to $\dot{{\boldsymbol{e}}}=-2g(\upsilon ,r)[{\upsilon }^{2}{\boldsymbol{r}}-({\boldsymbol{\upsilon }}\cdot {\boldsymbol{r}}){\boldsymbol{\upsilon }}]$. The eccentricity vector ${\boldsymbol{e}}$ is defined as always pointing in the direction of periapsis. This indicates that the eccentricity rate induced at periapsis is ${\dot{e}}_{{\rm{p}}}=\dot{{\boldsymbol{e}}}\cdot \hat{{\boldsymbol{r}}}\lt 0$, tending to decrease the eccentricity and circularize the orbit, while at apoapsis ${\dot{e}}_{{\rm{a}}}=\dot{{\boldsymbol{e}}}\cdot \hat{{\boldsymbol{r}}}\gt 0$, tending to increase the eccentricity (note that ${\boldsymbol{\upsilon }}\cdot {\boldsymbol{r}}=0$ at both periapsis and apoapsis). In addition, Equation (34) shows that for a massive body in an elliptical orbit, the eccentricity decrease rate due to dynamical friction is maximized at periapsis (f = 0) while the eccentricity increase rate is maximized at apoapsis ($f=\pi$). This implies that the eccentricity evolution depends on the relative time the massive body spends near periapsis and apoapsis along the orbit. Using this simplified picture, we calculate below the expected eccentricity evolution of the binary orbit.

Equation (29) gives the instantaneous change of the eccentricity

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{de}}{{dt}}\right)}_{{\rm{p}}}=\displaystyle \frac{2{(1-{e}^{2})}^{3/2}}{{n}^{3}{a}^{3}}\displaystyle \frac{{\epsilon }_{{\rm{p}}}}{{(1+e)}^{2}},\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{de}}{{dt}}\right)}_{{\rm{a}}}=-\displaystyle \frac{2{(1-{e}^{2})}^{3/2}}{{n}^{3}{a}^{3}}\displaystyle \frac{{\epsilon }_{{\rm{a}}}}{{(1-e)}^{2}},\end{eqnarray} \tag{ 36 }$$

with ${\epsilon }_{{\rm{p}}}$ and ${\epsilon }_{{\rm{a}}}$ the value of $\epsilon (r,v)$ calculated at periapsis and apoapsis, respectively. The time spent by the massive body near periapsis $({\rm{\Delta }}{\tau }_{{\rm{p}}})$ or apoapsis $({\rm{\Delta }}{\tau }_{{\rm{a}}})$ is proportional to

$$\begin{eqnarray}&&{\rm{\Delta }}{\tau }_{{\rm{p}}}\propto \displaystyle \frac{1}{{\upsilon }_{{\rm{p}}}}\sim \sqrt{\displaystyle \frac{1-e}{1+e}},\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{\tau }_{{\rm{a}}}\propto \displaystyle \frac{1}{{\upsilon }_{{\rm{a}}}}\sim \sqrt{\displaystyle \frac{1+e}{1-e}}.\end{eqnarray} \tag{ 38 }$$

The expected increase or decrease of the eccentricity can then be determined by the ratio

$$\begin{eqnarray}&&\zeta \equiv \displaystyle \frac{{\rm{\Delta }}{e}_{{\rm{p}}}}{{\rm{\Delta }}{e}_{{\rm{a}}}}=\displaystyle \frac{{({de}/{dt})}_{{\rm{p}}}\ {\rm{\Delta }}{\tau }_{{\rm{p}}}}{{({de}/{dt})}_{{\rm{a}}}\ {\rm{\Delta }}{\tau }_{{\rm{a}}}},\end{eqnarray} \tag{ 39 }$$

where ${\rm{\Delta }}{e}_{{\rm{p}}}$ and ${\rm{\Delta }}{e}_{{\rm{a}}}$ are the induced eccentricity changes near periapsis and apoapsis, respectively. If $| \zeta | =1$ the eccentricity remains constant, if $| \zeta | \gt 1$ the contribution near periapsis dominates and the eccentricity decreases, while if $| \zeta | \lt 1$ the apoapsis contribution dominates and the eccentricity increases.

We note that Equations (35) and (36) reduce to analytic expressions if we take into account only stars moving slower than the massive body. Under this consideration, we have $\beta (\xi )=\gamma (\xi )=0$ and Equations (16) and (17) lead to

$$\begin{eqnarray}&&{\epsilon }_{{\rm{p}}}\propto {(1-e)}^{-\gamma }{(1+e)}^{3/2}\ {}_{2}{F}_{1}^{p},\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&{\epsilon }_{{\rm{a}}}\propto {(1+e)}^{-\gamma }{(1-e)}^{3/2}\ {}_{2}{F}_{1}^{a}.\end{eqnarray} \tag{ 41 }$$

Substituting Equations (40) and (41) into Equations (35) and (36), we have

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{de}}{{dt}}\right)}_{{\rm{p}}}\propto -\displaystyle \frac{{(1-e)}^{-\gamma }}{{(1+e)}^{1/2}}\ \ \ {}_{2}{F}_{1}^{p},\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{de}}{{dt}}\right)}_{{\rm{a}}}\propto \displaystyle \frac{{(1+e)}^{-\gamma }}{{(1-e)}^{1/2}}\ \ \ {}_{2}{F}_{1}^{a},\end{eqnarray} \tag{ 43 }$$

where the hypergeometric function ${}_{2}{F}_{1}$ is always positive for all $0\lt e\lt 1$. The different signs in Equations (42) and (43) confirm the fact that dynamical friction tends to circularize the orbit at periapsis (${\dot{e}}_{{\rm{p}}}\lt 0$) while at apoapsis tends to increase it (${\dot{e}}_{{\rm{a}}}\gt 0$). Combining Equations (37), (38), (42), and (43), the ratio (39) is simplified to

$$\begin{eqnarray}&&| \zeta | ={\left(\displaystyle \frac{1-e}{1+e}\right)}^{3/2-\gamma }\displaystyle \frac{{}_{2}{F}_{1}^{p}}{{}_{2}{F}_{1}^{a}}.\end{eqnarray} \tag{ 44 }$$

As expected, for $\gamma =3/2$, we find $| \zeta | =1$ and the eccentricity remains constant; for $\gamma \gt 3/2$, the first term on the right-hand side of Equation (44) is $\gt 1$ and ${}_{2}{F}_{1}^{p}\gt {}_{2}{F}_{1}^{a}$. In this case, Equation (44) leads to $| \zeta | \gt 1$. On the other hand, for $\gamma \lt 3/2$ we have ${}_{2}{F}_{1}^{p}\lt {}_{2}{F}_{1}^{a}$ and $| \zeta | \lt 1$.

Figure 1 demonstrates the expected eccentricity evolution. When the contribution of fast stars is included, the critical value of γ below which the eccentricity increases is no longer 3/2 but $\approx 2$. This critical value of γ is found to depend slightly on the initial eccentricity while a smaller $\mathrm{ln}{\rm{\Lambda }}$ increases the parameter space for which the binary eccentricity increases. The fact that the critical value of γ in this case is greater than 3/2 reflects the fact that the relative contribution of the fast stars is larger near apoapsis where the massive body is moving slower than the local circular velocity. This results in an enhanced drag force at apoapsis and higher eccentricities.

Zoom In Zoom Out Reset image size

The orbital element time-evolution Equations (28)–(33) of a perturbed binary are in principle phase dependent and undergo physical oscillations with the orbital period. For a perturbation with a long timescale compared to the orbital period, these oscillations can be smoothed out by adopting orbit-averaging techniques. The orbit-averaged orbital element time-evolution equations then describe the secular evolution of the system.

In order to orbit-average a quantity along the orbit, we need to know how the true anomaly is changing over time. This is described by Equation (9), which shows that apart from the unperturbed Keplerian evolution described by the first term on the right-hand side of this equation, the true anomaly can also evolve due to possible precessions and specifically the periapsis precession $\dot{\omega }$ and the longitude of the ascending node precession $\dot{{\rm{\Omega }}}$. Given that due to dynamical friction we have $\dot{\omega }\ll 1$ and that $\dot{{\rm{\Omega }}}=0$, we can compute the secular evolution of the orbital elements neglecting the second term in Equation (9) and use

$$\begin{eqnarray}&&{df}=n\displaystyle \frac{{(1+e\cos f)}^{2}}{{(1-{e}^{2})}^{3/2}}{dt}\end{eqnarray} \tag{ 45 }$$

when integrating Equations (28) and (29) over the orbit. In order to incorporate higher-order effects in the time-evolution equations, we have to proceed one step further and include all the terms in Equation (9) when applying the orbit-averaging process. In this paper, we derive first-order secular time-evolution equations where the general orbit-averaging rule for the phase-dependent quantity $(...)$ is defined by the integral

$$\begin{eqnarray}&&\langle (...)\rangle =\displaystyle \frac{{(1-{e}^{2})}^{3/2}}{2\pi }{\int }_{0}^{2\pi }\ (...)\ \displaystyle \frac{{df}}{{(1+e\cos f)}^{2}}.\end{eqnarray} \tag{ 46 }$$

Under these considerations, the secular time-evolution equations for the semimajor axis and the eccentricity of the massive body can be written as

$$\begin{eqnarray}&&\langle \displaystyle \frac{{da}}{{dt}}\rangle =\displaystyle \frac{2{(1-{e}^{2})}^{2}}{\pi {n}^{3}{a}^{2}}{\int }_{0}^{2\pi }\ \displaystyle \frac{{(1+e\cos f)}^{-2}\epsilon (r,v)}{{(1+{e}^{2}+2e\cos f)}^{1/2}}{df},\end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray}\langle \displaystyle \frac{{de}}{{dt}}\rangle & = & \displaystyle \frac{2{(1-{e}^{2})}^{3}}{\pi {n}^{3}{a}^{3}}\\ & & \times {\displaystyle \int }_{0}^{2\pi }\ \displaystyle \frac{(e+\cos f)\epsilon (r,v)}{{(1+{e}^{2}+2e\cos f)}^{3/2}{(1+e\cos f)}^{2}}{df}.\end{eqnarray} \tag{ 48 }$$

Making use of Equation (47), we plot in Figure 2 the secular semimajor axis evolution of an inspiraling body with mass ratio $q=m/{M}_{\bullet }={10}^{-3}$. We see that for $\gamma \lt 1$ the dynamical friction timescale becomes much longer, although the contribution to dynamical fiction from the fast stars always increases the orbital decay rate. We discuss the implications of the long dynamical friction timescale for SMBH binaries in early-type galaxies in Sections 4 and 5 below.

Zoom In Zoom Out Reset image size

The secular evolution of the binary eccentricity is described by Equation (48). Using this equation we plot in Figure 3 the average fractional change in the quantity $R=1-{e}^{2}$ in one orbital decay time, i.e., $| a/\dot{a}|$. We plot this change as a function of the initial eccentricity.

Zoom In Zoom Out Reset image size

In the upper panel of Figure 3, we see that adding the contribution of fast stars makes the eccentricity higher compared to the case where only the slow stars were taken into account. The results shown in Figure 3 clearly imply that the change in eccentricity can be of order unity if $\gamma \lt 1$.

In the lower panel of Figure 3, we compare the results of Equation (48) with the eccentricity change predicted assuming that the stellar velocities follow a Maxwellian distribution. As expected, the use of a Maxwellian velocity distribution is inadequate to describe the evolution of the binary. More specifically, adopting a Maxwellian distribution always leads to a shorter timescale for the eccentricity evolution compared to what we obtain by using the distribution function (10).

We generate equilibrium N-body models of stars around an SMBH by adopting the truncated mass model

$$\begin{eqnarray}&&\rho (r)={\rho }_{0}{\left(\displaystyle \frac{r}{{r}_{\mathrm{infl}}}\right)}^{-\gamma }\exp (-r{/r}_{\mathrm{tr}}).\end{eqnarray} \tag{ 52 }$$

Monte Carlo initial positions and velocities were assigned to the N-body particles by numerically generating the distribution function corresponding to the isotropic equilibrium model of Equation (52), while at the same time taking into account the gravitational potential due to the central SMBH. A massive particle was placed in these models with mass $m\gg {m}_{\star }$, with corresponding mass ratio $q\ll 1$ at an initial galactocentric distance $r\lesssim {r}_{\mathrm{tr}}$.

The initial conditions were evolved forward in time using the direct-summation code ϕGRAPE (Harfst et al. 2007). This code uses a fourth-order Hermite integrator with a predictor-corrector scheme and hierarchical time steps. The performance and accuracy of the code depend both on the time-step parameter η and on the smoothing length . We set $\eta =0.01$ and $\epsilon =5\times {10}^{-4}$. With these choices, energy conservation was typically of order $0.1 \%$ over the entire length of the integration. The calculations were carried out in serial mode using graphics processing units combined with the sapporo library (Gaburov et al. 2009; Bédorf et al. 2015).

Figure 4 shows the evolution of the semimajor axis and eccentricity of the massive binary in N-body models with $\gamma =(0.6,1.5,2)$, two values of binary mass ratio $q\,=(2\times {10}^{-4},5\times {10}^{-5})$, and two different values for the mass of the field particles ${m}_{\star }=(5\times {10}^{-6},1\times {10}^{-5})$. In all the models shown in Figure 4, we set ${r}_{\mathrm{tr}}=0.2$ and place the secondary massive body at $r={r}_{\mathrm{tr}}$ with a velocity half the circular value. With this choice of parameters, the initial values of semimajor axis and eccentricity are ${a}_{0}\approx 0.1$ and ${e}_{0}\approx 0.7$.

Zoom In Zoom Out Reset image size

The resulting eccentricity and semimajor axis evolution in the simulations shown in Figure 4 are in good agreement with the theoretical predictions. In all cases, the semimajor axis evolves more rapidly than the eccentricity does. Indeed, we find no systematic change in the binary eccentricity over the simulated timescale for any value of γ. Although this behavior seems largely consistent with the predicted evolution, we also observe random-like variations of the orbital eccentricity. Such fluctuations in e are due to the hard scattering off surrounding stars which causes the angular momentum of the massive body to random walk with an amplitude that decreases with increasing mass ratio $m/{m}_{\star }$ (e.g., Matsubayashi et al. 2007). We note that in a real galaxy the mass ratio between the secondary massive body and field stars could be much larger than in our simulations, so such an effect would be essentially absent.

Figures 5 and 6 show simulations where the orbit of the massive body was followed until a had shrunk by a factor of $\approx 100$ below its initial value. In these additional simulations, we increased the mass of the secondary body and field particles, and set ${r}_{\mathrm{tr}}=1$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the simulations shown in Figure 5, ${a}_{0}\approx 1$ and ${e}_{0}\approx 0.3$. We see from the left panels that the orbital evolution of the binary at $a\gtrsim {a}_{{\rm{h}}}$ is essentially independent of m_⋆, or equivalently of the number of particles, N, used to represent the galaxy. At later times, $t\gtrsim 200$ or at $a\lesssim {a}_{{\rm{h}}}$, the binary hardening rate becomes significantly longer and shows a clear N-dependence, in the sense of slower hardening for larger N. In this phase, the hardening of the binary requires a repopulation of the depleted orbits through collisional loss-cone repopulation, which is an N-dependent process (e.g., Yu 2002).

The right panels of Figure 5 show the evolution of the massive binary in a set of integrations with the same value of N but for different values of q. From these plots we see that the dynamical friction timescale of a massive binary above ${a}_{\mathrm{crit}}$ scales approximately linearly with the mass of the secondary body (see also Figure 4). This is also expected given that the dynamical friction acceleration decreases proportionally with m, if we neglect the logarithmic dependence of the frictional drag on m through $\mathrm{ln}\,{\rm{\Lambda }}$. As before, the massive binary orbit is observed to stall at $\approx {a}_{{\rm{h}}}$.

Figures 5 and 6 show that the evolution of the binary eccentricity can be divided in two distinct phases. At separations $a\gt {a}_{\mathrm{crit}}$, the eccentricity of the binary changes steadily with time: it increases for $\gamma =(0.6,\,1)$ and decreases for $\gamma =2$. In the second phase, at $a\lt {a}_{\mathrm{crit}}$, the eccentricity established at early times tends to persist, remaining approximately constant as the orbit keeps shrinking. Because ${a}_{\mathrm{crit}}$ increases with m (see Equation (51)), lighter inspiraling bodies have higher eccentricities by the time their orbit has shrunk to $a\approx {a}_{{\rm{h}}}$. This fact appears evident in the top-right panel of Figure 5 where the general trend is clearly toward higher final eccentricities for smaller q.

As can be seen from Figures 4–6, the evolution of e and a at $a\gt {a}_{\mathrm{crit}}$ in the N-body simulations show good agreement with the theoretical predictions. For the cases we considered, the results of the N-body simulations confirm the expected qualitative result that for $\gamma \leqslant 1$, the eccentricity will increase during the inspiral while for $\gamma \geqslant 2$, the eccentricity will decrease. We conclude that the results of the N-body simulations support the correctness of the analytical treatment developed in Section 2 and consequently of Chandrasekhar's physical picture of dynamical friction.

Although the agreement with the theoretical prediction appears fairly good, a difference is also apparent: both a and e evolve more slowly in the N-body simulations than predicted. This results in a small displacement towards the right of the corresponding analytical curves shown in the figures. The discrepancy is caused by the fact that in our analytical treatment the contribution of the stars to the gravitational potential is ignored. Consequently, the massive body moves slower relative to the N-body, which leads to an artificially stronger frictional drag due to the ∝v⁻² dependence that appears in Equation (11). This also explains why for smaller initial separations the agreement appears to improve—for smaller a₀ the contribution of the stars to the gravitational potential can be more safely neglected given that the potential is more strongly dominated by the primary SMBH.

We have shown that after the binary semimajor axis has decreased to ${a}_{\mathrm{crit}}$, the eccentricity remains approximately constant in time. Therefore, the value of e at ${a}_{\mathrm{crit}}$ is also approximately the eccentricity the binary will have at the time it has decayed to ${a}_{{\rm{h}}}$. Below ${a}_{{\rm{h}}}$, the evolution becomes more uncertain. Scattering experiments typically suggest a slow but steady growth of eccentricity (Mikkola & Valtonen 1992).

The predicted eccentricity of the massive binary at ${a}_{\mathrm{crit}}$ is given in Figure 7 as a function of the binary mass ratio, different initial values of the orbital eccentricity, and for $\gamma =(0.6,\,1)$. In this plot we compare the theoretically predicted results with the results from the N-body simulations, which are shown as solid triangles and which were obtained as the average value of e at radii $a\lt {a}_{\mathrm{crit}}$. The agreement with the analytical predictions is again good. We see that for low mass-ratio binaries and a moderate initial eccentricity, ${e}_{0}\gtrsim 0.3$, the binary will reach $e\gtrsim 0.9$ by the time it has decayed to ${a}_{\mathrm{crit}}$. The results of this analysis indicate that due to dynamical friction, the eccentricity of a comparatively low-mass test body moving in the center of a massive galaxy will be high during the time it spends inside the sphere of influence of the central SMBH.

Zoom In Zoom Out Reset image size

The eccentricity evolution itself is especially important in the case of SMBH binaries since the energy losses due to gravitational-wave emissions depend strongly on e. How much the binary must shrink by stellar-dynamical processes before gravitational-wave emission takes over is very sensitive to the eccentricity of the binary. In what follows we apply the formalism developed in Section 2 and confirmed with N-body simulations to describe the evolution of SMBH binaries in early-type galaxies.

Modeling the main galaxy as a singular isothermal sphere, the timescale for a satellite SMBH of mass m to decay towards the center is (Binney & Tremaine 1987)

$$\begin{eqnarray}&&{T}_{\star }^{\mathrm{bare}}=17\displaystyle \frac{6.6}{\mathrm{ln}\,{\rm{\Lambda }}}\ {\left(\displaystyle \frac{{R}_{{\rm{e}}}}{10\mathrm{kpc}}\right)}^{2}\left(\displaystyle \frac{\sigma }{300\,\mathrm{km}\ {{\rm{s}}}^{-1}}\right)\left(\displaystyle \frac{{10}^{8}{M}_{\odot }}{m}\right)\,\mathrm{Gyr},\end{eqnarray} \tag{ 53 }$$

where R_e is the effective radius of the main galaxy.

Some previous work made use of Equation (53) to describe the formation and evolution of SMBH binaries during galaxy mergers (e.g., McWilliams et al. 2014). However, Equation (53) can lead to a significant overestimate of the real dynamical friction timescale and to erroneous conclusions about the survival timescale of SMBH pairs at the scale of R_e. Equation (53) neglects the fact that the satellite SMBH is brought in during the course of a galaxy merger and it will retain some of its host galaxy's stars until late in the merger process. Considering the stellar mass bound to the inspiraling SMBH leads to a significantly shorter timescale for the formation of a bound pair.

By assuming a strict proportionality between the mass of the stellar bulge of the satellite galaxy, M_s, and the mass of its central SMBH ${M}_{{\rm{s}}}={10}^{3}m$ (Merritt & Ferrarese 2001), and replacing m with M_s in Equation (53) gives

$$\begin{eqnarray}&&{T}_{\star ,1}^{\mathrm{gx}}=0.06\displaystyle \frac{2}{\mathrm{ln}\,{\rm{\Lambda }}^{\prime} }\ {\left(\displaystyle \frac{{R}_{{\rm{e}}}}{10\mathrm{kpc}}\right)}^{2}\left(\displaystyle \frac{\sigma }{300\,\mathrm{km}\ {{\rm{s}}}^{-1}}\right)\left(\displaystyle \frac{{10}^{8}{M}_{\odot }}{m}\right)\,\mathrm{Gyr},\end{eqnarray} \tag{ 54 }$$

where the argument of the Coulomb logarithm is taken to be ${\rm{\Lambda }}^{\prime} ={2}^{3/2}\sigma /{\sigma }_{{\rm{s}}},$ with ${\sigma }_{{\rm{s}}}$ the stellar velocity dispersion of the secondary galaxy.

Equation (54) does not consider that the satellite galaxy can be stripped of its stars due to the strong tidal field of the primary galaxy. Following Binney & Tremaine (1987), we assume that the satellite galaxy can also be modeled as a singular isothermal sphere, so that its mass is related to its velocity dispersion through the relation

$$\begin{eqnarray}&&{M}_{{\rm{s}}}(r)\approx \displaystyle \frac{1}{2G}{\alpha }^{2}{\sigma }_{{\rm{s}}}^{2}{r}_{t}.\end{eqnarray} \tag{ 55 }$$

Setting the Hill radius ${r}_{t}={({{GM}}_{{\rm{m}}}{r}^{2}/4{\sigma }^{2})}^{1/3}$ with M_m the mass of the main galaxy and $\alpha =2$ as appropriate for a sharp truncation, the dynamical friction timescale becomes

$$\begin{eqnarray}{T}_{\star ,2}^{\mathrm{gx}} & = & 0.15\displaystyle \frac{2}{\mathrm{ln}\,{\rm{\Lambda }}^{\prime} }\ \left(\displaystyle \frac{{R}_{{\rm{e}}}}{10\,\mathrm{kpc}}\right){\left(\displaystyle \frac{\sigma }{300\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{2}\\ & & \times {\left(\displaystyle \frac{100\mathrm{km}{{\rm{s}}}^{-1}}{{\sigma }_{{\rm{s}}}}\right)}^{3}\,\mathrm{Gyr}.\end{eqnarray} \tag{ 56 }$$

A good approximation of the timescale for a secondary SMBH to decay from R_e to the influence radius of the primary SMBH is

$$\begin{eqnarray}&&{T}_{\star }=\max ({T}_{\star ,1}^{\mathrm{gx}},{T}_{\star ,2}^{\mathrm{gx}}).\end{eqnarray} \tag{ 57 }$$

As the secondary SMBH enters the sphere of influence of the primary, the SMBHs are bound to each other, and the formulae given above, which are only strictly valid for a Maxwellian distribution of velocities and a self-gravitating cluster, can no longer apply. Nevertheless, most works in the past (e.g., Kelley et al. 2016) have neglected such a complication and applied Equation (53) until a_h.

Although such a simplification is reasonable for major mergers, we find that for $q\lesssim 0.1$ it necessarily leads to an erroneous evaluation of the binary evolution timescale as well as the evolution of its orbital eccentricity.

Using Equation (28), it can be shown that the characteristic dynamical friction timescale to decay from the primary SMBH influence radius, ${r}_{\mathrm{infl}}$, to a much shorter separation, $r=\chi {r}_{\mathrm{infl}}$, is

$$\begin{eqnarray}{T}_{\bullet }^{\mathrm{bare}} & = & 1.5\times {10}^{7}\displaystyle \frac{{[\mathrm{ln}{\rm{\Lambda }}\alpha +\beta +\delta ]}^{-1}}{(3/2-\gamma )(3-\gamma )}({\chi }^{\gamma -3/2}-1)\\ & & \times {\left(\displaystyle \frac{{M}_{\bullet }}{3\times {10}^{9}{M}_{\odot }}\right)}^{1/2}{\left(\displaystyle \frac{m}{{10}^{8}{M}_{\odot }}\right)}^{-1}{\left(\displaystyle \frac{{r}_{\mathrm{infl}}}{300\mathrm{pc}}\right)}^{3/2}\,\mathrm{years}.\end{eqnarray} \tag{ 58 }$$

The coefficients α, β, and δ can be easily computed from Equations (17)–(19), assuming a circular orbit, i.e., setting $\xi =1$ in these equations as justified by the fact that the orbital decay timescale is not significantly affected by the binary orbital eccentricity.

Equation (58) does not consider that part of the host galaxy can remain bound to the secondary SMBH even inside ${r}_{\mathrm{infl}}$. Assuming as before that the satellite galaxy can be modeled as a singular isothermal sphere, then the timescale to decay from ${r}_{\mathrm{infl}}$ to a smaller radius $r=\chi {r}_{\mathrm{infl}}$ can be obtained by replacing m with M_s(r), where now ${r}_{t}\approx {({M}_{{\rm{s}}}/2{M}_{\bullet })}^{1/3}r$. In this case, Equation (28) leads to

$$\begin{eqnarray}{T}_{\bullet }^{\mathrm{gx}} & = & 1.2\times {10}^{7}\displaystyle \frac{{[\mathrm{ln}{\rm{\Lambda }}^{\prime} \alpha +\beta +\delta ]}^{-1}}{{(3-\gamma )}^{2}}({\chi }^{\gamma -3}-1)\\ & & \times \left(\displaystyle \frac{{M}_{\bullet }}{3\times {10}^{9}{M}_{\odot }}\right){\left(\displaystyle \frac{100\mathrm{km}{{\rm{s}}}^{-1}}{{\sigma }_{{\rm{s}}}}\right)}^{3}\,\mathrm{years}.\end{eqnarray} \tag{ 59 }$$

The decay timescale from the influence radius of the primary SMBH is then

$$\begin{eqnarray}&&{T}_{\bullet }=\min ({T}_{\bullet }^{\mathrm{bare}},{T}_{\bullet }^{\mathrm{gx}}).\end{eqnarray} \tag{ 60 }$$

Note that Equation (60) is strictly valid only until the secondary hole reaches ${a}_{\mathrm{crit}}$ (see Equation (51)). Below this radius, the central cusp starts to be significantly modified by the secondary SMBH (Baumgardt et al. 2006) and time-dependent galaxy models are required in order to describe the details of the evolution of the binary orbit (Antonini & Merritt 2012). Here we ignore such complications and set $\chi ={a}_{\mathrm{crit}}/{r}_{\mathrm{infl}}$. Our choice is clearly conservative.

In order to quantify the error one would make by employing the standard Chandrasekhar formula, we plot in Figure 8 the dynamical friction timescale from Equation (58) with $\beta =\delta =0$ divided by ${T}_{\bullet }$. For $\gamma \lesssim 1$ the standard Chandrasekhar theory can lead to significant deviations from our more accurate formulation. Neglecting the fast-moving stars' contribution leads to an overestimate of the dynamical friction timescale that can be longer than ${T}_{\bullet }$ by about one order of magnitude for $\gamma \approx 0.5$. In Figure 8 we also compare our estimate to that obtained by assuming that the velocity distribution of the field stars follows a Maxwellian distribution. This latter approximation underpredicts the decay timescale by about one order of magnitude for $\gamma \approx 0.5$. For $\mathrm{ln}\,{\rm{\Lambda }}\gtrsim 6$ and/or $\gamma \gtrsim 1$ both approximations give a good estimate of the decay timescale, as the contribution of the fast stars to the drag force is smaller in this case.

Zoom In Zoom Out Reset image size

The phase of binary evolution determined by dynamical friction comes to an end when the binary's semimajor axis reaches $\sim {a}_{{\rm{h}}}$. Gravitational encounters will continue supplying stars to the binary at a rate that depends on the host galaxy morphology.

Dry mergers between luminous galaxies result in triaxial remnants, which leads to an efficient hardening of the binary (Khan et al. 2011; Vasiliev et al. 2014). After a major merger, a galaxy is likely to retain some degree of triaxiality or flattening during its subsequent evolution. The timescale to decay from ${a}_{{\rm{h}}}$ to coalescence is (Vasiliev et al. 2015)

$$\begin{eqnarray}&&{T}_{{\rm{h}}}\approx 1.2\times {10}^{9}{\left(\displaystyle \frac{{r}_{\mathrm{infl}}}{300\mathrm{pc}}\right)}^{\displaystyle \frac{10+4\psi }{5+\psi }},{\left(\displaystyle \frac{{M}_{\bullet }+m}{3\times {10}^{9}{M}_{\odot }}\right)}^{\displaystyle \frac{-5-3\psi }{5+\psi }},\end{eqnarray} \tag{ 61 }$$

$$\begin{eqnarray}&&{\phi }^{\displaystyle \frac{-4}{5+\psi }}{\left(\displaystyle \frac{4q}{{(1+q)}^{2}}\right)}^{\displaystyle \frac{3\psi -1}{5+\psi }}p(e)\,\mathrm{years},\end{eqnarray} \tag{ 62 }$$

with

$$\begin{eqnarray}p(e) & = & (1-{e}^{2})[k+(1-k){(1-{e}^{2})}^{4}],\\ k & = & 0.6+0.1\mathrm{log}([{M}_{\bullet }+m]/3\times {10}^{9}M\odot ).\end{eqnarray} \tag{ 63 }$$

The parameters ϕ and ψ parameterize the evolving hardening rate with values estimated from Monte Carlo simulations. In what follows we adopt $\phi =0.4$ and $\psi =0.3$, which are the values derived by Vasiliev et al. (2015) for triaxial galaxies.

We describe how much the binary must shrink by stellar-dynamical processes before gravitational-wave emission takes over using the ratio (Vasiliev et al. 2015):

$$\begin{eqnarray}\displaystyle \frac{{a}_{{\rm{h}}}}{{a}_{\mathrm{GW}}} & \approx & 55\times {\left(\displaystyle \frac{{r}_{\mathrm{infl}}}{30\mathrm{pc}}\right)}^{\displaystyle \frac{5}{10}}{\left(\displaystyle \frac{{M}_{\bullet }+m}{{10}^{8}{M}_{\odot }}\right)}^{-\displaystyle \frac{5}{10}}\\ & & \times f{(e)}^{\displaystyle \frac{1}{5}}{\left(\displaystyle \frac{4q}{{(1+q)}^{2}}\right)}^{\displaystyle \frac{4}{5}},\end{eqnarray} \tag{ 64 }$$

with

$$\begin{eqnarray}&&f(e)=\displaystyle \frac{{(1-e)}^{-\tfrac{7}{2}}}{1+\tfrac{73}{24}{e}^{2}+\tfrac{37}{96}{e}^{4}}\end{eqnarray} \tag{ 65 }$$

and ${a}_{\mathrm{GW}}$ the separation at which gravitational-wave radiation takes over. From Equation (64) one finds that ${a}_{{\rm{h}}}\lt {a}_{\mathrm{GW}}$ for $q\approx {10}^{-3}$ and moderate eccentricities. Binaries with mass ratio lower than this value never enter the hardening phase. The quoted approximation formula therefore breaks down in this case because these systems are never in a properly stellar-dynamical hardening regime. These binaries transit directly from the dynamical friction phase to the phase where gravitational-wave radiation leads to their rapid coalescence. Accordingly, in what follows we set ${T}_{{\rm{h}}}=0$ for $q\lesssim {10}^{-3}$.

Equation (61) implies that stellar-dynamical interactions are able to drive the binary to coalescence on a timescale $\leqslant 1$ Gyr in any triaxial galaxy, and that the coalescence timescale weakly depends on the mass of the binary and its mass ratio. Coalescence times also depend quite significantly on the binary eccentricity falling in the range from a few gigayears for almost circular binaries to $\approx {10}^{8}\,\mathrm{years}$ for very eccentric ones.

We note that Equation (61) was derived for major mergers, $q\gtrsim 0.1$. However, additional simulations performed using the Monte Carlo code RAGA (Vasiliev 2015) showed that Equation (61) works remarkably well for all binaries with $q\gtrsim {10}^{-2}$. At $q\approx {10}^{-3}$, Equation (61) starts to break down, but even for $q={10}^{-3}$ it overestimates by only a factor of $\approx 3$ the binary's coalescence timescale (E. Vasiliev 2017, private communication). Hence, unless otherwise specified, in what follows we adopt Equation (61) in the full range of mass ratios ${10}^{-3}\leqslant q\leqslant 1$. The fact that T_h is only an order of magnitude estimate for the lowest mass ratio binaries we considered does not affect our calculations below, given that the total lifetime of these binaries is largely dominated by their dynamical friction timescale.

Finally, we note that the assumption of triaxiality in the field-star distribution made in Equation (61) is obviously inconsistent with the spherical density model used to compute the dynamical friction timescales above. However, galaxy mergers produce remnants with deviations from isotropy that are small both in velocity and configuration space (Vasiliev et al. 2015). Hence, we might expect our estimates of the dynamical friction timescales to give a reasonable approximation.

Here we estimate the timescales associated with the three stages of binary evolution defined above for real galaxies and consider the possibility of "stalled" mergers in these systems, where a smaller satellite SMBH resides in the outer regions of the main host galaxy. We focus on massive early-type galaxies because (i) galaxy mergers are thought to play a crucial role in the late growth of these systems and (ii) these galaxies are often observed to have extended density profile cores within the sphere of influence of the SMBH, which implies a long dynamical friction timescale. For these reasons, luminous early-type galaxies are likely hosts of stalled satellite SMBHs. On the other hand, our results imply that in low-mass ellipticals, SMBH binaries have a short timescale, i.e., shorter inspiral times, possibly in the decreasing or near constant eccentricity regime.

We start by considering the case of a widely studied massive galaxy such as M87, for which ${M}_{\bullet }\approx 6.15\times {10}^{9}{M}_{\odot }$, $\sigma \approx 330\,\mathrm{km}\ {{\rm{s}}}^{-1}$ (Kormendy & Ho 2013), ${R}_{{\rm{e}}}\approx 9\,\mathrm{kpc}$ (Lauer et al. 1995), and ${r}_{\mathrm{infl}}\approx 300\,\mathrm{pc}$ (Marconi & Hunt 2003). Given these structural parameters, the calculation of the decay timescales can be simply performed by assuming that the secondary galaxy velocity dispersion is related to the mass of its central SMBH through the ${M}_{\bullet }\mbox{--}\sigma$ relation defined below in Equation (69) setting z = 0 (Ferrarese & Merritt 2000; Gebhardt et al. 2000)

An additional parameter that needs to be defined is γ, the slope of the deprojected spatial density profile of stars within ${r}_{\mathrm{infl}}$. We adopt the two representative values $\gamma =0.6$ and 1, as motivated by the facts that the M87 density profile is observed to be nearly flat inside ${r}_{\mathrm{infl}}$ and values $\gamma \leqslant 0.5$ are excluded by our assumption of kinematical isotropy.

The left panel of Figure 9 gives the orbital decay timescale as a function of the satellite SMBH mass inside M87 at the different stages of the evolution of the massive binary. The dynamical friction timescale to decay from the effective radius to the SMBH influence radius, ${T}_{\star }$, is shorter than $\approx {10}^{9}\,\mathrm{years}$ for $m\gtrsim 5\times {10}^{6},$ implying that SMBH binaries are unlikely to stall near R_e. The hardening timescale to reach coalescence from the hard binary separation is quite short, $\lesssim {10}^{9}\,\mathrm{years}$, even for the moderate eccentricity we adopted, e = 0.3. The dynamical friction timescale to decay from ${r}_{\mathrm{infl}}$ to ${a}_{{\rm{h}}}$ can instead be extremely long. For $m\lesssim {10}^{8}{M}_{\odot }$, ${T}_{\bullet }$ becomes longer than either T_⋆ or ${T}_{{\rm{h}}}$, and it is longer than $10\,\mathrm{Gyr}$ for $m\lesssim {10}^{6}{M}_{\odot }$ ($m\lesssim {10}^{7}{M}_{\odot }$ ) when $\gamma =1$ ($\gamma =0.6$).

Zoom In Zoom Out Reset image size

In the right panel of Figure 9, we show the timescale of orbital decay ${T}_{\star }$ for the sample of 31 core-Sérsic early-type galaxies in Dullo & Graham (2015; purple points). These systems are bright elliptical and lenticular galaxies with extended density profile cores. Dullo & Graham (2015) give the measured structural parameters of these galaxies, including R_e and σ, which we used to compute ${T}_{\star }$. This timescale is plotted as a function of the primary SMBH mass for $q={10}^{-3}$. Larger symbols are systems for which a direct SMBH measurement is available in the literature (Kormendy & Ho 2013). The resulting ${T}_{\star }$ is $\lesssim {10}^{9}\,\mathrm{years}$. The hardening timescale ${T}_{{\rm{h}}}$ (dotted–dashed line) also appears to be quite short $\lesssim {10}^{9}\,\mathrm{years}$, and, as expected, weakly dependent on the primary SMBH mass. Finally, the hatched red region gives ${T}_{\bullet }$, which was computed through Equation (60) by setting (e.g., Merritt et al. 2009)

$$\begin{eqnarray}&&{r}_{\mathrm{infl}}\approx 35{\left(\displaystyle \frac{{M}_{\bullet }}{{10}^{8}{M}_{\odot }}\right)}^{0.56}\,\mathrm{pc},\end{eqnarray} \tag{ 66 }$$

with ${M}_{\bullet }$ given by Equation (69) at z = 0.

The analysis shown in the right panel of Figure 9 confirms and generalizes our statement that the dynamical friction timescale inside the influence radius of an SMBH can become of the order of the Hubble time in luminous spheroids. We considered that ${T}_{\bullet }$ can be significantly larger than either ${T}_{\star }$ or ${T}_{{\rm{h}}}$, and can become longer than the Hubble time even for relatively high mass-ratio binaries not just in the core of M87 but for most galaxies in the sample.

In Figure 10, we investigate further the dependence of the lifetime of an SMBH binary on its total mass and mass ratio. The total lifetime of an SMBH binary is given by ${t}_{{\rm{L}}}={T}_{\star }+{T}_{\bullet }+{T}_{{\rm{h}}}$. We plot the lifetime of a massive black hole binary in the three different stages of the binary evolution as a function of the binary total mass and mass ratio. For the effective radius of the host galaxy, we use ${R}_{{\rm{e}}}=1.35{({M}_{\bullet }/{10}^{8}{M}_{\odot })}^{0.73}\,\mathrm{kpc}$, where the mass dependence and normalization were obtained from the observed mass–radius relation of galaxies at redshift zero (Figure 7 in Forbes et al. 2008 for their sample of elliptical galaxies), and after expressing the galaxy mass as ${M}_{{\rm{m}}}={10}^{3}{M}_{\bullet }$, assuming the latter relation holds at all redshifts.

Zoom In Zoom Out Reset image size

Based on Figure 10, the amount of time the secondary SMBH spends to decay from the effective radius of the galaxy R_e to the influence radius of the central black hole is short, typically less than ∼3 Gyr. Adding the time the binary spends in the dynamical friction phase until ${a}_{\mathrm{crit}}$ leads to significantly longer binary lifetimes. Figure 10 shows that including the dynamical friction timescale for high total-mass and low mass-ratio binaries expands the parameter space for long-lived binaries with lifetimes greater than ∼3 Gyr, while there is a considerable amount of these binaries that have lifetimes greater than ∼10 Gyr. This implies that binaries with high total mass and low mass ratio are expected to live long in the evolutionary stage between the SMBH influence radius and hardening radius of the binary. Based on Figure 7, we also expect that these binaries will have high eccentricities. The final phase in the evolution of an SMBH binary before coalescence is characterized by the hardening timescale that we add in the binary total lifetime in the right panel of Figure 10. The hardening timescale contributes mostly to the lifetime of high mass-ratio binaries, always being shorter than $\approx 1\,\mathrm{Gyr}$.

The observational evidence of small-orbit SMBH binaries is still scarce. Electromagnetic tracers of post-merger galaxy cores are hard to identify. This makes the study of the post-merger dynamics of binary SMBH systems difficult; pairs of SMBHs are usually observed during the early stages in their dynamical evolution when still at ∼kpc separations. Although observing SMBH binaries at parsec scales is challenging since they cannot be spatially resolved in optical and X-ray, more work in detecting SMBH binaries at subkiloparsec scales ($\lesssim 100\ \mathrm{pc}$) is needed since discovering such systems and obtaining a census of their population properties would serve as a test of galaxy evolution models and would provide valuable constraints for stellar and gas dynamical models for the decay of the binary orbit.

If accretion is triggered along the course of the merger, direct evidence of SMBH binaries comes from the presence of active galactic nuclei (AGNs). Accretion by SMBHs can give rise to a rich phenomenology that goes from dual to binary and off-set AGNs, in the radio-loud or quite mode, according to their dynamics, habitat, and merger type. Specifically, if only one SMBH is accreting, it will be visible as an off-center AGN accreting from a small disk. If both binary components of the SMBH binary accrete, then a dual or double AGN might be observed.

In this work, we are interested in mergers that occur in gas-poor environments where collision processes are at play. Core galaxies are promising systems in which to search for SMBH binaries since they have experienced a number of minor mergers during their lifetime and are minimally affected by extinction. For example, a number of displaced SMBHs have been observed as off-set AGNs in host elliptical galaxies (e.g., Lena et al. 2014). A characteristic example is the galaxy M87 (NGC 4486) which has a measured 7.7 ± 0.3 pc projected displacement of the photocenter relative to the AGN (Batcheldor et al. 2010; Janowiecki et al. 2010). As a fraction of the rather large core radius ${r}_{{\rm{c}}}$, the weighted mean displacement is only $\approx 0.01{r}_{{\rm{c}}}$. The expected observed displacement ${\rm{\Delta }}r$ of the primary SMBH is defined as ${\rm{\Delta }}r=q\ R/(1+q)$ for $q={10}^{-4}$ and $q={10}^{-3}$, with R the orbital separation between the two SMBHs. Based on Figure 13, the secondary SMBHs are expected to stall at a distance ${a}_{\mathrm{stall}}\approx 0.1{r}_{\mathrm{infl}}$. At this radius we have ${\rm{\Delta }}r\approx q\ {a}_{\mathrm{stall}}/(1+q)$, which gives ${\rm{\Delta }}r/{r}_{\mathrm{infl}}\sim {10}^{-4}$ $(q={10}^{-3})$ and ${\rm{\Delta }}r/{r}_{\mathrm{infl}}\sim {10}^{-5}$ $(q={10}^{-4})$. These values are small and not consistent with the observed off-center displacements mentioned in Lena et al. (2014) which are of order $\sim {10}^{-1}\mbox{--}{10}^{-2}$ for the massive elliptical galaxies shown in Figure 13 following the sample of Lena et al. (2014). However, if we assume that the secondary, rather than the primary, SMBH is accreting, then the predicted offsets, in this case the galactocentric distance of the satellite SMBH, appear to be consistent at least with those observed in NGC 4278 and NGC 5846.

Although difficult to spatially resolve in the optical and X-ray, giant ellipticals often host strong radio sources. For example, a pair of AGNs within the elliptical host galaxy 0402 + 379 and with a projected separation of 7.2 pc has recently been discovered (Rodriguez et al. 2006, 2009; Burke-Spolaor 2011). This is the closest binary SMBH yet discovered within a core elliptical and may be the tip of the iceberg of SMBH binaries with parsec-scale separations. Binary AGNs within host core ellipticals at subkiloparsec separations like the one in Rodriguez et al. (2006) could be explained as stalled holes in a slowly evolving orbit inside low-density cores.

Early observations have suggested the presence of double nuclei in core elliptical galaxies with the second nucleus at off-center subkiloparsec separations. A characteristic example is the double nucleus of the core elliptical galaxy NGC 5419 with the second nucleus at an off-center separation ∼70 pc (Mazzalay et al. 2016). The projected separation between the two nuclei is smaller than the estimated SMBH influence radius (${r}_{\mathrm{infl}}\approx 100$ pc) and much smaller than the core radius (${r}_{{\rm{c}}}\approx 500$ pc). The same appears to be true in the case of NGC 4696, which has a secondary nucleus at ∼30 pc from the center (Laine et al. 2003), which is slightly smaller than the estimated SMBH influence radius (${r}_{\mathrm{infl}}\approx 40$ pc) and much smaller than the core radius (${r}_{{\rm{c}}}\approx 250$ pc). These features in NGC 4696 have been interpreted as evidence of a recent minor merger with a gas-rich galaxy (Sparks et al. 1989; Farage et al. 2010). The case of NGC 6876, which has a double nucleus at ∼30 pc from the center (Lauer et al. 2002; Dullo & Graham 2012), is similar. This is slightly smaller than the estimated SMBH influence radius (${r}_{\mathrm{infl}}\approx 40$ pc) and much smaller than the galaxy inner core radius (${r}_{{\rm{c}}}\approx 119$ pc).

More recent observations suggest the presence of multiple nuclei located (in projection) within the core elliptical galaxy A2261-BCG (Bonfini & Graham 2016). The estimated mass of the core and the central SMBH in A2261-BCG is ${M}_{\mathrm{core}}\sim {M}_{\bullet }\sim 2\times {10}^{10}$. This implies that if an SMBH is present at the center of this galaxy, then such nuclei are residing well inside its influence radius.

The double/multiple nuclei considered above have galactocentric separations that are below ${r}_{\mathrm{infl}}$ but well above the separation at which the binary can be considered a "hard" binary. We conclude that the equilibrium and stability of such double/multiple nuclei could be understood in terms of the long dynamical friction timescale within the influence radius of the host galaxy SMBH predicted by our models.

We conclude by noting that the stalled satellites predicted by our models are likely to be on highly eccentric orbits (see Section 2). The highly eccentric orbit of the satellite can also have interesting observational consequences. As the satellite is disrupted, it will form an eccentric disk-like structure around the nucleus of the primary galaxy. A characteristic example could be the nucleus of the low-luminosity elliptical galaxy NGC 4486B (companion to M87). NGC 4486B is unusual for a galaxy of its luminosity in having a well-resolved core. Two brightness peaks are observed in the core separated by $\approx 12$ pc (Lauer et al. 1996). Neither peak is coincident with the galaxy photocenter, which falls between them. This double-peak structure has been interpreted as evidence for an eccentric disk of stars where the peaks would be the ansae of the disk orbiting an SMBH. The disk might be related to the disruption of a star cluster on an eccentric orbit by the tidal field of the primary SMBH; its eccentric nature would naturally result by the dynamical friction process in the flat density core. For example, taking ${M}_{\bullet }=2\times {10}^{8}{M}_{\odot }$ (Lauer et al. 1996), an inspiraling star cluster with total mass $m={10}^{6}{M}_{\odot }$ will reach the center in $\approx 0.2\mbox{--}0.5\,\mathrm{Gyr}$ starting from ${r}_{\mathrm{infl}}$. If we assume a reasonable value for the cluster core radius of $\approx 2\,\mathrm{pc}$ then the cluster will be disrupted by the SMBH tidal field when it reaches a separation of $\approx 10\,\mathrm{pc}$ from the center. This distance appears to be consistent with the observed offsets of the two nuclear brightness peaks in NGC 4486B.