Greatly enhanced merger rates of compact-object binaries in non-spherical nuclear star clusters

TL;DR

This paper demonstrates that non-spherical nuclear star clusters significantly enhance the merger rates of compact-object binaries through a cluster-driven Lidov-Kozai mechanism, providing a new channel for gravitational wave sources.

Contribution

It introduces a novel mechanism where non-spherical clusters increase eccentricity excitation, boosting merger rates of compact binaries compared to spherical cluster models.

Findings

Eccentricity excitation is most efficient when precession timescales are comparable.

Cluster perturbations increase extreme eccentricity fractions by 10-100 times.

Estimated merger rates are up to 15 Gpc^{-3}yr^{-1} for BH-BH binaries.

Abstract

The Milky Way and a significant fraction of galaxies are observed to host a central Massive Black Hole (MBH) embedded in a non-spherical nuclear star cluster. We study the secular orbital evolution of compact-object binaries in these environments and characterize the excitation of extremely large eccentricities that can lead to mergers by gravitational radiation. We find that the eccentricity excitation occurs most efficiently when the nodal precession timescale of the binary's orbit around the MBH due to the non-spherical cluster becomes comparable (within a factor of ~10) to the timescale on which the binary is torqued by the MBH due to the Lidov-Kozai (LK) mechanism. We show that in this regime the perturbations due to the cluster increase the fraction of systems that reach extreme eccentricities ($1-e\sim10^{-4}-10^{-6}$) by a factor of ~10-100 compared to the idealized case of a spherical cluster, increasing the merger rates of compact objects by a similar factor. We identify two main channels that lead to this extreme eccentricity excitation: (i) chaotic diffusion of the eccentricities due to resonance overlap; (ii) cluster-driven variations of the mutual inclinations between the binary orbit and its center-of-mass orbit around the MBH, which can intensify the LK oscillations. We estimate that our mechanism can produce black hole-black hole and black hole-neutron star binary merger rates of up to $\approx 15$ $\rm{Gpc}^{-3}yr^{-1}$ and $\approx 0.4$ $\rm{Gpc}^{-3}yr^{-1}$, respectively. Thus, we propose the cluster-enhanced Lidov-Kozai mechanism as a new channel for the merger of compact-object binaries, competing with scenarios that invoke isolated binary evolution or dynamical formation in globular clusters.

Tables

Table 1. Table 1 Results of the population synthesis models of massive binary stars.

	${ϵ}_z$	$f_2(\%)$	$f_2(\%)$	$f_m(\%)$	$f_m(\%)$
		$\alpha =1.7$	$\alpha =2.3$
BH-BH	0.001	4.5	2.5	0.17	0.17
…	0.1	4.5	2.5	1.5	5.8
…	0.3	4.5	2.5	2.6	12
…	1	4.5	2.5	3.9	17
BH-NS	0.001	1.3	0.87	0	0
…	0.1	1.3	0.87	0.072	0.51
…	0.3	1.3	0.87	0.26	1.2
…	1	1.3	0.87	0.29	1.5
NS-NS	0.001	0.22	0.28	0	0
…	0.1	0.22	0.28	0.045	0.045
…	0.3	0.22	0.28	0.18	0.27
…	1	0.22	0.28	0.41	0.50

Equations

$${ϵ}_{\mathrm{oct}}=\left(\frac{m_1-m_2}{m_1+m_2}\right)\left(\frac{a_{\mathrm{in}}}{a_{\mathrm{out}}}\right)\frac{e_{\mathrm{out}}}{1-e_{\mathrm{out}}^2}$$

$$\Phi (\mathbf{r})$$

$${\rho }_s(\mathbf{r})={\rho }_{s,0}{\left(\frac{r}{r_0}\right)}^{-\gamma }\left[1+\widetilde{ϵ}-{\widetilde{ϵ}}_z\frac{(\mathbf{r}\cdot {\hat{\mathbf{n}}}_z{)}^2\cdot }{r^2}-{\widetilde{ϵ}}_y\frac{(\mathbf{r}\cdot {\hat{\mathbf{n}}}_y{)}^2}{r^2}\right],$$

$${\widetilde{ϵ}}_i={ϵ}_i\frac{\gamma (5-\gamma )}{(2-\gamma )(3-\gamma )},$$

$$\widetilde{ϵ}=\frac{2}{(2-\gamma )(3-\gamma )}({ϵ}_y+{ϵ}_z).$$

$${\left(\frac{\tilde{a}}{\tilde{c}}\right)}^2={\left(\frac{\tilde{b}}{\tilde{c}}\right)}^2=\frac{1+\widetilde{ϵ}-{\widetilde{ϵ}}_z}{1+\widetilde{ϵ}}=\frac{1-{ϵ}_z}{1+{ϵ}_z},$$

$$\mathcal{E}=1-{\left(\frac{1-|{ϵ}_z|}{1+|{ϵ}_z|}\right)}^{1/2}.$$

$${\mathbf{r}}_{\mathrm{in}}$$

$${\mathbf{r}}_{\mathrm{out}}$$

$${\tau }_{\text{mbox}LK}$$

$${\rho }_s(r)=0.8\times 10^5{M}_{\odot }\mathrm{p}{\mathrm{c}}^{-3}{\left(\frac{r}{1\mathrm{pc}}\right)}^{-1.3},$$

$${\tau }_{\mathrm{sc},\mathrm{z}}$$

$${ϵ}_{\mathrm{sc},\mathrm{z}}$$

$${ϵ}_{\mathrm{oct},\mathrm{sc}}$$

$${\tau }_{2b}$$

$${\tau }_{\mathrm{coll}}$$

$$\Theta =\frac{4\mathrm{G}(m_1+m_2+m)}{(\sqrt{3}\sigma {)}^2{a}_{\mathrm{in}}}=0.08\left(\frac{m_1+m_2+m}{20M_{\odot }}\right)$$

$${\left(\frac{\sigma }{100\mathrm{km}{\mathrm{s}}^{-1}}\right)}^{-2}{\left(\frac{a_{\mathrm{in}}}{10\mathrm{AU}}\right)}^{-1}.$$

$$N_{\text{mbox}LK,coll}$$

$$\times$$

$${\tau }_{\mathrm{evap}}$$

$$N_{\text{mbox}LK,evap}$$

$${\tau }_{\mathrm{GW}}$$

$$1-e_{\mathrm{merger}}\lesssim 3\times 10^{-5}\times \left(\frac{a_{\mathrm{out}}(1-e_{\mathrm{out}}^2{)}^{1/2}}{0.1\text{mbox}pc}\right).$$

$${ϵ}_{\mathrm{GR}}\equiv \frac{{\tau }_{\text{mbox}LK}}{{\tau }_{\mathrm{GR}}}$$

$$\frac{d{\mathbf{e}}_{\mathrm{out}}}{dt}$$

$$\frac{d{\hat{\mathbf{j}}}_{\mathrm{out}}}{dt}$$

$$\frac{{\bar{\Phi }}_{\mathrm{rot}}}{{\varphi }_0}$$

$${\varphi }_0=\frac{3G{M}_{\mathrm{BH}}{a}_{\mathrm{in}}^2}{4a_{\mathrm{out}}^3(1-e_{\mathrm{out}}^2{)}^{3/2}}.$$

Full text

1. Introduction

Most nearby galaxies contain a compact stellar cluster residing at their kinematical and photometric center (e.g., Neumayer et al. 2011; Turner et al. 2012; den Brok et al. 2014; Georgiev & Böker 2014; Georgiev et al. 2016). These nuclear clusters (NCs) have masses in the range $\approx 10^{5}-10^{8}~{}M_{\odot}$ and half mass radii of only a few parsecs, making them the most massive and densest stellar clusters observed in the local Universe (Forbes et al. 2008; Misgeld & Hilker 2011). Most NCs are believed to host a central massive black hole (MBH); this is the case for the Milky Way, the location of a $\sim 3\times 10^{7}M_{\odot}$ NC (Schödel et al. 2014, 2017) and a $\approx 4\times 10^{6}M_{\odot}$ MBH (Ghez et al. 2008; Gillessen et al. 2009). Examples of galaxies hosting both a compact NC and an AGN are known across a wide range of galaxy masses (Seth et al. 2008; Neumayer & Walcher 2012), suggesting that NCs and MBHs often co-exist in galaxies.

The dynamics of stars near a MBH in a NC is characterized by high velocity dispersion, extreme stellar densities and short relaxation time (Merritt 2013). On even shorter timescales stars move around the MBH on nearly Keplerian orbits. Binaries that reside in this environment can therefore interact with the MBH for many orbits and perturbations resulting from their interaction with the MBH can accumulate for a long time before the binary properties are significantly altered by encounters with other field stars. In particular, binaries with highly inclined orbits with respect to their orbit around the MBH can be strongly affected by secular processes, namely by Lidov-Kozai (LK) oscillations, which periodically change their eccentricities and inclinations (Lidov 1962; Kozai 1962). At the highest eccentricities during a LK oscillation the components of the binary can merge (Antonini et al. 2010; Antonini & Perets 2012). This process can lead to a number of potentially important phenomena, including the merger of stellar binaries which can result in the formation of massive stars and/or G2-like objects (Antonini et al. 2011; Prodan et al. 2015; Stephan et al. 2016; Liu et al. 2017); and the formation of X-ray binaries. Particularly exciting is the possibility that many, or maybe even most, of the merging compact object binaries that are soon to be detected by Advanced LIGO (aLIGO) (Abbott et al. 2016a, b) are formed near MBHs through the LK mechanism (Antonini & Perets 2012; VanLandingham et al. 2016).

In the standard theory of hierarchical triple systems, for an eccentric outer orbit, the inner binary can reach high eccentricities and undergo chaotic evolution of its orientation for a large range of initial mutual inclinations of inner to outer orbit (e.g., Ford et al. 2000; Lithwick & Naoz 2011; Katz et al. 2011; Li et al. 2014). However, in the environment of a MBH, extremely high eccentricities are difficult to achieve for most binaries (Antonini & Perets 2012). One reason for this is that the octupole-order terms are effectively vanishing in this case (see also Section 3.1.1 below). Following Naoz et al. (2013), we quantify the importance of the octupole-level coefficients using the quantity

[TABLE]

with $m_{1}$ and $m_{2}$ the mass of the binary components, $a_{\rm in}$ the binary semi-major axis, and $e_{\rm out}$ and $a_{\rm out}$ the eccentricity and semi-major axis of the orbit of the binary around the MBH, respectively. Even for large values of $e_{\rm out}\sim 0.9$, and for $a_{1}\sim 1$ AU the octupole-order terms are effective (i.e., $\epsilon_{\rm oct}\gtrsim 10^{-3}$) only within a short distance from the galactic center ($\lesssim 0.01\rm pc$), which for a Milky Way like galaxy is inhabited only by a tiny fraction of the overall population of stars bound to the MBH (see Figure 14 in Feldmeier-Krause et al. 2015). 111This conclusion is in contradiction with pevious work from Stephan et al. (2016) who argued for a key contribution of the octupole-order terms. This in spite of their Figure 1 (top panel) which shows that the octupole dynamics can only become relevant inside $\approx 0.0\rm 1pc$, in agreement with our conclusion. Moreover, Stephan et al. (2016) and most previous work (e.g., Prodan et al. 2015) did not take into account the mass precession of the outer orbit which further suppresses the contribution of the octupole-order terms (See section 3.1.1).

A direct consequence of the negligible octupole contribution is that the magnitude of the inner orbit’s angular momentum has a lower limit equal to its constant component, $j_{z}$, along the total angular momentum. Thus, for isotropic initial conditions, the probability of achieving the high eccentricities required for a gravitational-wave mediated merger, $1-e_{\rm in}\lesssim 10^{-5}$ (see Eq. [25]), is only $\lesssim 0.003$.222Only binaries with inclinations relative to their orbit around the MBH in the range $[89^{\circ}.8,90^{\circ}.2]$ can reach $1-e_{\rm in}\lesssim 10^{-5}$.

From our discussion and consistent with the population synthesis study by Antonini & Perets (2012), it follows that in a spherical star cluster, where $j_{z}=\rm const$, the contribution from the standard LK mechanism to the number of merging binaries in galactic nuclei is negligible. However, there is ample evidence that a substantial fraction of NCs, including the one in the Milky Way, have a non-spherical shape (e.g., Georgiev & Böker 2014; Chatzopoulos et al. 2015; Feldmeier-Krause et al. 2017). One important feature of axisymmetric (and triaxial) systems is that the angular momenta (orientation and/or magnitude) of stellar orbits are not conserved. In particular, the orientation of the stellar orbits around the MBH can change gradually with time, implying that $j_{z}\neq\rm const$. Similarly, close to the MBH angular momenta can also evolve due to resonant relaxation (Rauch & Tremaine 1996). Consequently, the available fraction of phase space for achieving high values of $e_{\rm in}$ can potentially be much larger than that obtained in the standard LK theory. Thus, the binary merger rate in a real NC can be well above that obtained for an idealized model of a spherical star cluster.

In this paper, we consider the evolution of binaries residing in a non-spherical NC hosting a central MBH and argue that the potential from the cluster causes precession of the binary’s orbit in timescales that are comparable to the LK timescale of typical binaries (Eqs. [10]-[12]). As a result, $j_{z}$ can change dramatically and even lead to chaotic behavior so that extreme eccentricities can be attained for a wide range of orbital configurations. We quantify the how this mechanism enhances the merger rates of binaries comprising neutron stars (NSs) and black hole (BHs) in galactic centers.

The paper is organized as follows. In §2 we discuss our basic set up with the equations of motion fully expressed in the Appendix. In §3 we describe the relevant timescales in the problem. Our results are presented in §§5 and 4, where the former section presents the general dynamical behavior and simplified model, while the latter section shows the statistical results from a population of binaries changing different parameters. We summarize the results of these two sections in §5.6. In §6 we discuss the implications of our results for the merger rates of compact-object binaries. We summarize our main results in §7.

2. Basic equations

We model the evolution of a stellar binary with component masses of $m_{1}$ and $m_{2}$ with semi-major axes $a_{\rm in}$ orbiting a MBH of mass $M_{\rm BH}$ in a Keplerian orbit with semi-major axis $a_{\rm out}$. For the sake of brevity, in some expressions we drop the sub-index “in” when referring to the inner orbit. Thus, $e\equiv e_{\rm in}$, $a\equiv a_{\rm in}$, $\omega\equiv\omega_{\rm in}$, and so on.

We express the potential per unit mass from the MBH of mass $M_{\rm BH}$ and stellar cluster potential that governs the dynamics of a binary with masses $m_{1},m_{2}\ll M_{\rm BH}$ as

[TABLE]

The second term in the this equation represents the potential from a triaxial cluster. From Poisson’s equation we get the following density distribution (see e.g, Chandrasekhar 1969; Merritt 2013):

[TABLE]

where

[TABLE]

with $i=y,z$ and

[TABLE]

For $r\sim r_{0}$ the density distribution corresponds corresponds roughly to an ellipsoidal mass distribution and it is useful to relate the quantities above with axes ratios or ellipticities. In particular, for an axisymmetric potential ($|\epsilon_{y}|\ll|\epsilon_{z}|$) and $\gamma=1$, approximations widely used in this paper, we can express the axis ratios of the mass distribution as

[TABLE]

where $\tilde{a}$, $\tilde{b}$, and $\tilde{c}$ correspond to the axis length scales of the ellipsoidal mass distribution along $\hat{\bf n}_{x}$, $\hat{\bf n}_{y}$, and $\hat{\bf n}_{z}$ respectively. Thus, for $\epsilon_{z}>0$ the density distribution represents an oblate spheroid ($\tilde{a}=\tilde{b}<\tilde{c}$), while $\epsilon_{z}<0$ it represents a prolate spheroid ($\tilde{a}=\tilde{b}>\tilde{c}$). Moreover, we can define the ellipticity333We define the ellipticity as $1-q$ with $q$ being the ratio between the minor axis to the major axis of the ellipsoid. Also, the eccentricity of the ellipsoid can be computed as $\left[2|\epsilon_{z}|/(1+|\epsilon_{z}|)\right]^{1/2}$. of the density distribution for either oblate or prolate distributions as

[TABLE]

Thus, $\epsilon_{z}\simeq 0.05-0.5$ is equivalent to a distribution with $\mathcal{E}\simeq 0.05-0.4$, while the limiting case of a fully flattened distribution with $\mathcal{E}\to 1$ corresponds to $|\epsilon_{z}|\to 1$. Throughout this paper, we will take the coefficients $\epsilon_{y}$ and $\epsilon_{z}$ as free parameters and without loss of generality we will assume that they are positive444The sign of $\epsilon_{z}$ only changes the direction of the nodal precession and has no effect for the statistical distributions of orbital elements of the binaries.. We remark that formally the mass distribution in Equation (3) does not correspond to ellipsoidal mass distribution (unless $\gamma=0$) and the ellipticity is only defined locally for $r\sim r_{0}$ as a reference.

Within the sphere of influence of the MBH the orbits are nearly Keplerian and we can treat the force of the stellar cluster as a small perturbation that gradually torques the orbit, while not changing its energy–the so-called secular approximation. As depicted in Figure 1, we use the following coordinate system:

[TABLE]

with the MBH at the origin (${\bf r}_{\rm BH}=0$). The Keplerian orbits are fully specified by the specific angular momentum vectors ${\bf j}_{i}=\sqrt{1-e_{i}^{2}}\hat{\bf j}_{i}$ and the eccentricity vectors ${\bf e}_{i}=e_{i}\hat{\bf e}_{i}$ with $i=\{{\rm in},{\rm out}\}$, while $\hat{\bf q}_{i}=\hat{\bf j}_{i}\times\hat{\bf e}_{i}$ completes the orthogonal triad. The equations of motion are fully described in the Appendix.

3. Relevant timescales in the problem

The environment of a NC is a complicated, multi-timescale, system and there are many ingredients coming into play when considering the evolution of a binary. In what follows, we quantify the most relevant (typically the shortest) timescales in the problem and compare this to the LK timescale (e.g., Holman et al. 1997):

[TABLE]

where $P_{\rm in}$ is the period of the inner binary orbit.

We start in this section by discussing some of the main dynamical processes that can affect the evolution of the inner orbit of the binaries as well as their orbit around the MBH.

For our galactic center we shall assume that $M_{\rm BH}=4\times 10^{6}M_{\odot}$ and that the spherical component of density profile (setting $\tilde{\epsilon}=\tilde{\epsilon}_{z}=\tilde{\epsilon}_{y}=0$ in Eq. [3]) of the nuclear star cluster is given by

[TABLE]

consistent with recent observational constraints (Fritz et al. 2016; Schödel et al. 2017).

3.1. Evolution of the binary orbit around the MBH

3.1.1 Cluster-induced nodal and apsidal precession

A key new ingredient in our paper is the effect that the cluster gravitational potential has on the binary orbit around the MBH. There are two main dynamical effects that we describe next.

Nodal precession and possible changes

of $j_{z}$

Provided that the cluster is not fully spherically symmetric ($\epsilon_{z}>0$ in Equation 2), then ${\bf j}_{\rm out}$ will precess around ${\bf\hat{n}}_{z}$. From Equations (28) and (A15) we find that the characteristic timescale for this nodal precession is

[TABLE]

and for a generic density profile we have that $\tau_{\rm sc,z}\propto(a_{\rm out})^{\gamma-3/2}$. Since $\tau_{\mbox{\tiny LK}}\propto a_{\rm out}^{3}$, we generally have that $\tau_{\mbox{\tiny LK}}\ll\tau_{\rm sc,z}$ for small $a_{\rm out}$, while the opposite regime is reached for large $a_{\rm out}$. Thus, there is region in which $\tau_{\mbox{\tiny LK}}\sim\tau_{\rm sc,z}$.

Assuming the density distribution in Equation (11) we can express the nodal precession rate ($=1/\tau_{\rm sc,z}$) in units of the LK timescale as

[TABLE]

As we will show in §4, this ratio of timescales dramatically changes the behavior of the inner binary and for $\epsilon_{\rm sc,z}\gtrsim 0.1$ significant changes of $j_{z}=\sqrt{1-e_{\rm in}^{2}}\cos i_{\rm in,out}=\sqrt{1-e_{\rm in}^{2}}~{}\hat{\bf j}_{\rm in}\cdot\hat{\bf j}_{\rm out}$ occur, leading to new interesting phenomena. Geometrically the changes in $j_{z}$ are due to changes in $i_{\rm in,out}$ from the fast precession of ${\bf j}_{\rm out}$ around ${\bf\hat{n}}_{z}$.

Finally, we note that in triaxial clusters ($\epsilon_{y}>0$ in Equation 2) the orbit will also precess around ${\bf\hat{n}}_{y}$ with characteristic timescale $\tau_{\rm sc,y}\equiv\tau_{\rm sc}/\epsilon_{y}$.

Apsidal precession and suppression

of the octupole-level dynamics

Another important effect from the cluster is the orbit apsidal precession induced by the dominant spherical mass distribution. This important effect was neglected in most previous work (Prodan et al. 2015; Stephan et al. 2016). For the purposes of our model, the main dynamical effect is to quench the non-axisymmetric contribution (octupole-level or higher order terms) from the outer orbit. In other words, if the outer obit precesses fast compared to $\tau_{\mbox{\tiny LK}}$, its time-averaged potential will be effectively axisymmetric, even for $e_{\rm out}>0$.

To quantify this effect, we define the ratio between the octupole timescale, $\sim\tau_{\mbox{\tiny LK}}/\epsilon_{\rm oct}$ with $\epsilon_{\rm oct}$ from Equation (1), and the precession due to the spherical cluster using from Equation (A15) as

[TABLE]

where we assumed a cluster density profile given by Equation (11). For $\epsilon_{\rm oct,sc}\gg 1$ we expect that the octupole-level perturbations do not contribute significantly to the excitation of large eccentricities of the binary.

We note that regardless of the cluster-driven apsidal precession, the octupole-level perturbations are expected to be only relevant if $\epsilon_{\rm oct}\gtrsim 10^{-3}$. This condition can be written as $a_{\rm in}/[10\mbox{AU}]\gtrsim a_{\rm out}/[0.02\mbox{pc}]$ for $m_{1}\gg m_{2}$ and $e_{\rm out}\sim 0.5$. In order to compare which condition dominates at suppressing the octupole, we write the condition $\epsilon_{\rm oct,sc}\gtrsim 1$ as $a_{\rm in}/[10\mbox{AU}]\gtrsim(a_{\rm out}/[0.02\mbox{pc}])^{8.4/5}$. Thus, given that the latter condition scales more steeply with $a_{\mathrm{out}}$ we expect that the cluster is the dominant source of suppression for $a_{\mathrm{out}}\gg 0.02$ pc, while the condition from triple dynamics dominates for $a_{\mathrm{out}}\ll 0.02$ pc. In either case, we expect that the octupole is highly suppressed for $a_{\mathrm{out}}\gg 0.02$ pc.

We conclude that the octupole-level potential from the outer orbit is expected to have a negligible contribution at exciting extreme eccentricities at separations of $a_{\rm out}\sim 0.1-1$ pc. It can, however, be an important ingredient in the dynamics inside $a_{\rm out}\lesssim 0.01$ pc.

3.1.2 Two-body and resonant relaxation

In the dense stellar environment of a NC, stellar orbits can be significantly affected by their frequent interactions with field stars, which cause their orbital energy and angular momentum to random walk. The timescale associated with this process, referred to as two-body relaxation, is (Spitzer 1987):

[TABLE]

where $\Lambda^{\prime}\sim M_{\rm BH}/m$ and $\Lambda\sim(m_{1}+m_{2})/m$ are the Coulomb logarithms with $m$ the typical mass of the stars in the cluster. Note that $\Lambda^{\prime}\gg\Lambda$, since encounters on all scales contribute to the relaxation process. This timescale is clearly much longer than the other relevant timescales in the problem (e.g., $\tau_{\rm coll};\ \tau_{\rm LK}$). A massive binary will also segregate to the center via dynamical friction. The associated timescale is $\approx\tau_{\rm 2b}m/(m_{1}+m_{2})$; an expression roughly valid for a steep density profile cusp (Dosopoulou & Antonini 2016). Thus, dynamical friction acts on a timescale which is typically much longer than $\tau_{\rm coll}$, while it can become of order $\tau_{\rm evap}$ for relatively massive binaries ($\gtrsim 10\ M_{\odot}$), with small orbital semi-major axis $a_{\rm in}\lesssim 1\rm AU$, and that reside at large galactocentric distances ($\sigma\lesssim 100\rm km\ s^{-1}$). Such binaries are, however, unlikely to be affected by LK-like oscillations given the strong quenching caused by relativistic precession (see §3.2.2).

Based on our discussion, we conclude that collisional (gravitational-encounter-driven) relaxation processes are likely to become important only after the binaries have already evaporated from the NC. However, very near a MBH orbital angular momenta can evolve as a consequence of resonant relaxation (Rauch & Tremaine 1996; Hopman & Alexander 2006). Due to the finite number of stars in the vicinity of a MBH and because their Keplerian orbits maintain their orientation fixed in space for a long time, torques acting on a test star from all other stars do not cancel exactly. The resulting $\sqrt{N}$ deviation from spherical symmetry induces a net residual torque on the test star orbit.

The frequency of the orbital plane precession associated with resonant relaxation is: $\tau_{\mathrm{RR}}\approx\tau_{\mathrm{sc}}\sqrt{N}$ with $N$ approximately the number of stars within a sphere of radius $a_{\mathrm{out}}$. Comparison with the timescale associated with the nodal precession in a triaxial potential (Eq. [12]) shows that $\tau_{\mathrm{RR}}/\tau_{\mathrm{sc,z}}\approx\epsilon_{z}\sqrt{N}$. From Equation (11) setting the typical stellar mass in the nuclear cluster equal to $0.5M_{\odot}$ and $\gamma=1$, we have $N\approx 10^{6}\left(a_{\mathrm{out}}/\mathrm{1pc}\right)^{2}$. It follows that precession due to the aspherical cluster will dominate over resonant relaxation at radii $a_{\mathrm{out}}\gtrsim 10^{-3}\mathrm{pc}/\epsilon_{z}$. If we set $\epsilon_{z}=0.1$, then resonant relaxation can be ignored at $a_{\mathrm{out}}\gtrsim 0.01\mathrm{pc}$.

The evolution of binaries due to resonant relaxation couple with the LK process is a complex phenomenon and will be investigated in future work. Here we limit to notice that there is clearly a connection between the evolution of orbits in non-spherical clusters and due to resonant-relaxation. As also noted in Merritt et al. (2011) and Merritt & Vasiliev (2011), during the coherent phase of resonant relaxation, the gravitational potential from the star cluster can be represented in terms of a multipole expansion; if the lowest-order non-spherical terms in that expansion coincide with the potential generated by a triaxial cluster, then the behavior of orbits in the coherent regime would be the same as that derived below for orbits in a triaxial nucleus. We therefore expect the reorientation of orbital planes due to resonant relaxation to result into an evolution similar to the one described below in the case of a non-spherical star cluster.

3.2. Evolution of the binary inner orbit

3.2.1 Binary evaporation

Due to the large velocity dispersion near a MBH, most binaries in a nuclear cluster are very soft, meaning that $|E_{b}|/(m\sigma^{2})\ll 1$, where $m$ is the typical stellar mass, $E_{b}=-G(m_{1}+m_{2})^{2}/(2a_{\rm in})$, and $\sigma$ is the 1D stellar velocity dispersion. Encounters with field stars will make binaries softer and disrupt them (Heggie 1975).

The timescale on which a binary experiences an encounter with a field star is $t_{\rm coll}=(\nu\Sigma\sqrt{3}\sigma)^{-1}$, where $\nu$ is the number density of stars, and $\Sigma=\pi r_{\rm coll}^{2}\left[1+2G(m_{1}+m_{2}+m)/(\sqrt{3}\sigma)^{2}r_{\rm coll}\right]$ is the interaction cross-section for periapsis distances $\leq r_{\rm coll}$. The timescale for a binary to experience a close encounter with a field star is then (e.g., Binney & Tremaine 2008):

[TABLE]

with

[TABLE]

In the expressions above, we set $r_{\rm coll}=a_{\rm in}/2$. In more distant encounters the effect of the field particle on one of the binary components almost cancels the effect on the other component (e.g., Heggie 1975). While more distant encounters can still (secularly) change the eccentricity and orientation of the binary, we reserve a careful treatment of this potentially important effect to future work.

Most relevant for our work is the ratio between the collision timescale and the Lidov-Kozai timescale:

[TABLE]

as long as $N_{\mbox{\tiny LK}}\gtrsim 1$ the effect of encounters only becomes important after many LK cycles.

Strictly speaking, since we neglect the effect of encounters with field stars in our calculations, on timescales that are longer than $\tau_{\rm coll}$ our treatment might not longer be valid. However, the timescale over which impulsive encounters will fully disrupt a binary is longer than $\tau_{\rm coll}$. After only multiple encounters ($N_{\mbox{\tiny LK}}>1$), the orbital elements of the binary might change significantly and eventually lead to the evaporation of the binary. The evaporation time can be estimated as (Binney & Tremaine 2008):

[TABLE]

The ratio between the evaporation timescale and the Lidov-Kozai timescale for the galactic center ($m\nu=\rho_{s}$ from Eq. [11]) and assuming \sigma=280\rm km\ s^{-1}$$\sqrt{0.1\mbox{pc}/a_{\mathrm{out}}} (e.g., Kocsis & Tremaine 2011) is

[TABLE]

which is an estimate of the number of LK cycles that the binary can undergo before it fully “evaporates” due to encounters with stars in the NC.

We remark that the effect from encounters in the binary orbit is a stochastic process and the statistical fluctuations governing it are not taken into account in our work (see Bradnick et al. 2017 for a recent work that incorporates these stochastic processes). Instead we study how the results behave as a function of the number of cycles (see §5.3).

3.2.2 Relativistic perturbations

Two compact-objects reaching sufficiently close to each other lose energy by emitting gravitational wave (GW) radiation.

The merger timescale mediated by gravitational radiation is (Peters 1964):

[TABLE]

which at very large eccentricities ($1-e^{2}\simeq 2[1-e]$), relevant for our work, can be written as

[TABLE]

In order to have a rapid merger such that the binary orbit shrinks significantly in one LK cycle we require the change in semi-major axis due to GW loses, $\tau_{\rm GW}$, to be shorter than the change in the periapsis distance $a_{\rm in}(1-e_{\rm in})$ due to MBH, which at higher eccentricities is $\simeq\sqrt{1-e_{\rm in}^{2}}\tau_{\mbox{\tiny LK}}$. For BH binaries with $m_{1}+m_{2}=10M_{\odot}$ initially at $a_{\rm in}=10$ AU, this condition for rapid merger reads:

[TABLE]

It is well-known that LK oscillations can be quenched by additional sources of apsidal precession of the inner binary orbit (e.g., Blaes et al. 2002). This can prevent the binary from reaching the high eccentricities required for a rapid merger to occur. We parametrize the importance of relativistic precession by the ratio between the LK timescale and the GR precession timescale as (e.g., Fabrycky & Tremaine 2007)

[TABLE]

As we show below in Section §5.5, the eccentricities do no reach $1-e\lesssim 10^{-3}$ ($1-e\lesssim 10^{-5}$) for $\epsilon_{\rm GR}>0.1$ ($>0.01$).

3.3. Available phase-space for

eccentricity excitation.

In Figure 2 we show the regions in $a_{\rm in}-a_{\rm out}$ space where our model can potentially excite large eccentricities for two values of $m_{1}+m_{2}$.

Provided that the binary can undergo several LK cycles before evaporating ($N_{\mbox{\tiny LK,evap}}\gg 1$; Eq. [21]), the main limiting timescale to achieve extreme eccentricities is the relativistic precession. For reference, we show the shaded region that corresponds to $\epsilon_{\rm GR}<0.1$ (Equation [26]). As we show later in §5.5 we require that $\epsilon_{\rm GR}\lesssim 0.1$ ($\lesssim 0.01$) to reach $1-e_{\rm max}\lesssim 10^{-3}$ ($1-e_{\rm max}\lesssim 10^{-5}$).

The regions where the effect from the cluster can potentially enhance the LK mechanism are indicated in diagonal black lines. As we will show, the cluster mainly enhances the excitation of large eccentricities when the timescale at which it forces the nodal precession of the outer orbit becomes comparable to the LK timescale. We parametrize this with the ratio of timescales $\epsilon_{\rm sc,z}\equiv\epsilon_{z}\tau_{\mbox{\tiny LK}}/\tau_{\rm sc}$ defined in Equation (13). The black diagonal lines indicate $\epsilon_{\rm sc,z}=\{0.1,1,10\}$ and should serve only as reference since they assume a cluster with a particular density distribution and shape ($\epsilon_{z}=0.1$ or ellipticity of $\sim 0.1$).

From this Figure we conclude that our mechanism is most likely to play a significant role for binaries at separations of $a_{\rm out}\gtrsim 0.05$ pc from the MBH and with semi-major axes $a_{\rm in}\gtrsim 1$ AU. However, we note that as $a_{\rm out}$ increases, the number of LK cycles allowed before the binary evaporates due to encounters from passing stars $N_{\mbox{\tiny LK, evap}}$ in (Eq. [21], red lines) decreases. Since the behavior at $\epsilon_{\rm sc,z}\sim 1$ is often chaotic, the number of cycles limits the phase-space exploration to reach $e\to 1$ and, therefore, the merger rates.

Finally we show the region where our secular equations of motion in Appendix are expected to become invalid when the orbit loses enough angular momentum ($1-e\lesssim 10^{-4}$) due to the break down of the secular approximation (region above green line from Equation 18 in Antonini et al. 2014). When this occurs the orbital elements of the inner binary can change in dynamical timescales due to the torque from the outer orbit, while the semi-major axes remains fixed. This behavior is not captured by our secular code, but we expect that it does not invalidate our conclusion that the merger rates are high when $\epsilon_{z}\sim 1$ because the angular momentum changes in dynamical timescales generally lead to extremely large values of the eccentricities, as previously noted by several authors (e.g., Antonini & Perets 2012; Katz & Dong 2012; Antognini et al. 2014; Antonini et al. 2014).

4. Dynamical behavior

We describe the dynamical behavior that arises from including the effect from the cluster potential in the dynamics of binaries. We compare our results with those from evolving the binaries in the absence of a cluster (or a having a fully spherical cluster), in which case the dynamics reduces to the well-known Lidov-Kozai mechanism at quadrupole-order (expansion up to $a_{\rm in}^{2}/a_{\rm out}^{3}$).

We focus on the regime in which the main effect of the cluster is to force the outer orbit to undergo nodal precession (see Figure 1 for reference), but not change the magnitude of its angular momentum substantially. We remark that the model in the Appendix (Eqs. [A10]-[A13]) is not subject to this simplification and it can handle strong variations of the angular momentum of the outer orbit due to the triaxial potential.

4.1. Simplified two-degree of freedom model

We consider the regime in which the outer orbit has enough angular momentum and/or the cluster is not sufficiently flattened (or triaxial), so that the magnitude of the outer orbit angular momentum $j_{\rm out}=(1-e_{\rm out}^{2})^{1/2}$ does not change substantially. In practice, for mathematical simplicity, we consider the following limiting regime:

• •

nearly circular outer orbit: $e_{\rm out}\ll 1$;

• •

weakly axisymmetric cluster potential: $\epsilon_{y}\ll\epsilon_{z}\ll 1$ in Equation (2).

In this limit, the Equation (A12) becomes

[TABLE]

where $\tau_{\rm sc}$ is the timescale of variations of the outer orbit due to the spherical component of the stellar cluster given in Equation (A15). This last equation implies that ${\bf e}_{\rm out}$ only undergoes apsidal precession and $e_{\rm out}$ remains small. In turn, this small $e_{\rm out}$ implies that Equation (A13) becomes:

[TABLE]

meaning that $\hat{\bf j}_{\rm out}$ is simply undergoing nodal precession around $\hat{\bf n}_{y}$ and $\hat{\bf n}_{z}$ with rates proportional to $\epsilon_{y}$ and $\epsilon_{z}$, respectively. Simplifying things further, for an axisymmetric cluster ($\epsilon_{y}\ll\epsilon_{z}$) we have

[TABLE]

meaning that $\hat{\bf j}_{\rm out}$ simply precesses around $\hat{\bf n}_{z}$ at a constant rate $-\epsilon_{z}/\tau_{\rm sc}\left(\hat{\bf n}_{z}\cdot\hat{\bf j}_{\rm out}\right)$. Conveniently, we can move to the rotating coordinate system that co-precesses with the outer orbit and write the potential in Equation (A9) due to the MBH and the star cluster with normal vector $\hat{\bf n}_{z}$ as (e.g., Tremaine & Yavetz 2014)

[TABLE]

where

[TABLE]

We note that the dimensionless term multiplying $\left(\hat{\bf n}_{z}\cdot\hat{\bf j}_{\rm out}\right)\left({\bf{j}}\cdot\hat{\bf n}_{z}\right)$ is simply $\epsilon_{\rm sc,z}\equiv\epsilon_{z}\tau_{\mbox{\tiny LK}}/\tau_{\rm sc}$ defined in Equation (13). This parameter will be used throughout this paper and, as we will show, it critically determines the dynamical behavior of the inner binary.

In order to make the connection with a “perturbed” Lidov-Kozai mechanism more familiar to the reader, we write the potential in Equation (LABEL:eq:phi) in terms of orbital elements. By doing this, we arrive at the following dimensionless two-degree of freedom Hamiltonian:

[TABLE]

where we assume that $m_{1},m_{2}\ll M_{\rm BH}$, so we can define our coordinate system in the rotating frame as $\hat{\bf x}\equiv\hat{\bf e}_{\rm out}$ and $\hat{\bf z}\equiv\hat{\bf j}_{\rm out}$. For the sake of brevity, the variables of the inner orbit have no sub-index and all the angles are defined relative to the outer orbit (e.g., $i\equiv i_{\rm in,out}$). The angles $i_{\rm out,z}^{\prime}$ and $\Omega_{\rm out,z}^{\prime}$ are the inclination and the longitude of the ascending node of the cluster symmetry axis $\hat{\bf n}_{z}$ relative to the outer orbit in the rotating frame, both of which are constant in this frame. Note that this Hamiltonian can be re-written in terms of the dimensionless555The actual Delaunay’s variables have the momenta $j$ and $j_{z}$ multiplied by $m_{1}m_{2}/(m_{1}+m_{2})\sqrt{G(m_{1}+m_{2})a}$. canonical coordinate-momentum pairs: $\left\{\omega,j=\sqrt{1-e^{2}}\right\}$ and $\left\{\Omega,j_{z}=\sqrt{1-e^{2}}\cos i\right\}$.

This Hamiltonian shows explicitly that the momentum coordinate $j_{z}={\bf j}\cdot\hat{\bf j}_{\rm out}=\sqrt{1-e^{2}}\cos i$ can change due to the cluster potential and it varies as:

[TABLE]

where the unit time is given in Lidov-Kozai timescales, $\tau=t/\tau_{\mbox{\tiny LK}}$. This variation of $j_{z}$ is a crucial new ingredient in the dynamics of the stellar binary and as we will show in the next sections it can lead to the excitation of extremely high eccentricities for a wide range of initial conditions.

We note that this effect of forcing $j_{z}$ is similar the effect from high-order terms in the triple interaction potential (octupole and higher-order, see Naoz 2016), but unlike the three-body interaction, this forcing is not limited to happen in timescales much longer than $\tau_{\mbox{\tiny LK}}$666The octupole timescale is longer than $\tau_{\mbox{\tiny LK}}$ by a factor $\sim a_{\rm out}(1-e_{\rm out}^{2})/(a_{\rm in}e_{\rm out})$, while the timescale from higher-order terms are even longer than the octupole timescale. (i.e., $\epsilon_{\rm sc,z}$ is not necessarily a small parameter). In other words, for reasonable parameters of the triple system and the cluster potential, $\epsilon_{\rm sc,z}$ can take a wide range of values, including values much larger than unity.

Numerical example

In Figure 3 we show an example of the orbital evolution of a binary integrating the full model (Eqs. [A10]-[A13]) for a case with $\epsilon_{\rm sc,z}=0.3$ (black lines).

From panel a we observe that the eccentricities reach values of $1-e\sim 10^{-5}$ compared to $e_{\rm max}=\sqrt{1-5/3\cos^{2}(72^{\circ})}=0.917$ for the case with $\epsilon_{\rm sc,z}=0$ (red line). The binary flips its orbit relative to the outer orbit from prograde to retrograde and back during the large eccentricity periods (panel b). This behavior resembles the orbit flips in the three-body interactions with octupole-level perturbations (e.g., Katz et al. 2011; Lithwick & Naoz 2011) and are due to the forcing of $j_{\rm z}$, which periodically changes sign (panel c). Given the parameters of this example (initial $i_{\rm out,z}=70^{\circ}$ and $\epsilon_{\rm sc,z}=0.3$), from Equation (33) we have that $dj_{z}/d\tau\sim 0.1\sin i\sin\Omega$, meaning that $j_{\rm z}$ changes by $\lesssim 0.1$ in every LK cycle, consistent with small-amplitude oscillations observed in panel c. The argument of periapsis librates initially at around $90^{\circ}$ (panel d) and quasi-periodically switches to libration around $-90^{\circ}$ and back.

In conclusion, this first example shows that the effect of having nodal precession of the outer orbit can be understood (at least in this specific regime) as having an external perturbation in the three-body interaction that forces $j_{\rm z}$.

4.2. Dynamical behavior for

different values of $\epsilon_{\rm sc,z}$

In this section, we discuss the dynamical behavior for different values of $\epsilon_{\rm sc,z}$ and use Figures 4 and 5 to illustrate the general trends.

In Figure 4 we show a set of characteristic examples of the evolution of the binary for small, intermediate, and large values of $\epsilon_{\rm sc,z}$ (left to right panels, respectively). We integrate the Equations of motion from the full model (Eqs. [A10]-[A13]) with different initial conditions (see the corresponding labels).

In Figure 5 we show the surfaces of section of the two-degree of freedom Hamiltonian $\tilde{\mathcal{H}}_{\rm rot}$ in Equation (32) in $\omega-(1-e^{2})^{1/2}$ space to show how the increase in the perturbation strength parametrized by $\epsilon_{\rm sc,z}$ leads to the overlap between the LK resonance and a new resonance that appears at lower $j$. This resonance overlap, and possibly others, gives rise to the widespread chaos observed as $\epsilon_{\rm sc,z}$ approaches unity. We remark that a detailed study of the chaotic behavior in our model is beyond the scope of our work as $\tilde{\mathcal{H}}_{\rm rot}$ depends on two parameters ($\epsilon_{\rm sc,z}$ and $i_{\rm out}$) and sections for various energy manifolds should be considered. Thus, this figure should only be considered as a representative example.

To create the surfaces of section, we evolve the binaries using the equivalent potential in the rotating frame from Equation (LABEL:eq:phi) for which the equations of motion can be written as

[TABLE]

and vary $\epsilon_{\rm sc,z}$, while fixing $\tilde{\mathcal{H}}_{\rm rot}=0.1$, $i_{\rm out,z}^{\prime}=45^{\circ}$, and $\Omega_{\rm out,z}^{\prime}=0$.

4.2.1 Limit of small $\epsilon_{\rm sc,z}$:

slow nodal precession

For $\epsilon_{\rm sc,z}\ll 1$, the Hamiltonian in Equation (32) roughly coincides with the integrable Lidov-Kozai Hamiltonian. Thus, we expect that the binaries starting from low-eccentricity orbits and mutual inclinations $i_{\rm in,out,0}=\cos^{-1}\hat{\bf{j}}_{0}\cdot\hat{\bf j}_{\rm out}$ reach a maximum eccentricity roughly given by

[TABLE]

In the left panels of Figure 4 we show the evolution for a binary with $\epsilon_{\rm sc,z}=0.03$ and find that indeed the maximum eccentricity coincides with this expression (panel d). As expected, the evolution of the mutual inclinations $i_{\rm in,out}$ follow the regular LK oscillations, while the evolution of the inclination relative to $|\hat{\bf n}_{z}|$, $i_{\rm in,z}$, follows from the nodal precession of the binary around $\hat{\bf j}_{\rm out}$ and the LK inclination oscillations.

An exception to the evolution depicted above occurs when $|j_{\rm z}|$ is small, $|j_{\rm z}|\lesssim\epsilon_{\rm sc,z}$. From Equation (33), we estimate that for small $|j_{\rm z}|$ ($|\sin i_{0}|\sim 1$) the changes in $j_{\rm z}$ over one LK cycle are

[TABLE]

and, therefore, the cluster potential can drive the systems with $|j_{\rm z}|\lesssim\epsilon_{\rm sc,z}$ to cross $j_{\rm z}=0$. By crossing $j_{\rm z}=0$, the eccentricities can greatly overcome $e_{\rm max}$ in Equation (36). We will come back to this point below, in the population synthesis study of §5.2 (see yellow line in Figure 8).

In turn, the surface of section in Figure 5 with $\epsilon_{\rm sc,z}=0.01$ (panel a) shows that the trajectories are similar to those from the well-known integrable Lidov-Kozai Hamiltonian with the librating paths for low $j$ at around the fixed point at $|\omega|=90^{\circ}$ and circulating trajectories outside this region. The small perturbation due to the cluster potential introduces two main changes. First, the trajectories around the separatrix become chaotic (gray dots) similar to the ocupole-level pertubations in the LK mechanism seen by Holman et al. (1997) with similar perturbation strengths with $\epsilon_{\rm oct}\sim 0.008$ in Equation (1). Second, a new libration region appears at $j\lesssim 0.1$ with a fixed point at $|\omega|=90^{\circ}$ and $j\sim 0.03$. From the Equation (34), we observe that a fixed point ($d{\bf e}/d\tau=0$) involving the term due to the cluster $\epsilon_{\rm sc,z}\ll 1$ have to occur for $j\sim\epsilon_{\rm sc,z}$, explaining why the resonance appears at such small values. As we discuss next, this new resonance plays a major role in the evolution of the system as it expands with increasing $\epsilon_{\rm sc,z}$ and leads to global chaos.

Some intuition regarding the chaos around the separatrix might be gained by looking at the dynamically equivalent and time-dependent one-degree-of-freedom Hamiltonian in Equation (38), where $\hat{\bf j}_{\rm out}(t)$ changes steadily due to the nodal precession (Eq. [39]). Here, the variation of $\hat{\bf j}_{\rm out}(t)$ changes the location of the separatrix in time, while repeated separatrix crossings can act to drive stochasticity (Lichtenberg & Lieberman 1983). Thus, chaos might be expected to occur only for the regions in phase-space that are close enough to the separatrix and such distance depends on the strength of the perturbation.

In panel b, we show the surface of section for $\epsilon_{\rm sc,z}=0.05$ and observe that the phase-space structure is similar to the case with $\epsilon_{\rm sc,z}=0.01$ and the two libration islands remain separated in action-space. The main effect of increasing $\epsilon_{\rm sc,z}$ is to increase the area of the low$-j$ separatrix, whose respective fixed point moves from $j\sim 0.05$ for $\epsilon_{\rm sc,z}=0.01$ to $j\sim 0.2$ for $\epsilon_{\rm sc,z}=0.05$. Also, the area covered by the chaotic region around the high$-j$, or Lidov-Kozai, sepatrarix increases by a factor of a few and reaches down to $j\sim 0.4$. As a result, the distance between the chaotic region and the separatrix of the low$-j$ resonance decreases from $\Delta j\sim 0.4$ for $\epsilon_{\rm sc,z}=0.01$ to $\Delta j\sim 0.1$ for $\epsilon_{\rm sc,z}=0.05$. As we increase $\epsilon_{\rm sc,z}$ this distance approaches zero and the resonance overlaps leading to widespread chaos as observed in panels c and d. We discuss this overlapped regime next.

Finally, we remark that the chaos in the small $\epsilon_{\rm sc,z}$ regime occurs only around the separatrix and does not promote the excitation of large eccentricities. The excitation to extreme eccentricities happen for orbits with $|j_{\rm z}|\lesssim\epsilon_{\rm sc,z}$, in a similar fashion as the octupole-level perturbations of the LK mechanism that drives $j_{\rm z}$ towards 0.

4.2.2 Intermediate values of $\epsilon_{\rm sc,z}$

For intermediate values of $\epsilon_{\rm sc,z}$ the behavior is more complicated than the limiting cases of small and large $\epsilon_{\rm sc,z}$. In particular, the dynamical evolution of the binary can be governed by chaotic behavior in a large fraction of the available phase-space.

In the middle panels of Figure 4 we show the evolution for a binary with $\epsilon_{\rm sc,z}=1$ and observe that the inclinations behave in a chaotic-looking fashion (panel b). Most remarkably, starting from a low mutual inclination, $i_{\rm in,out}=15^{\circ}$, the binary orbit becomes retrograde after $\sim 5$ cycles. In turn, the evolution in $e-\omega$ space (panel e) shows rapid intermittence between librating and circulating trajectories, typical of chaotic systems, and the eccentricity reaches $1-e\sim 10^{-4}$ after $\sim 5$ cycles.

In Figure 5 we show the surfaces of section for $\epsilon_{\rm sc,z}=0.2$ (panel c) and $\epsilon_{\rm sc,z}=1$ (panel d) in which wide chaotic regions can be observed.

For the case with $\epsilon_{\rm sc,z}=0.2$ we observe that there exists three libration islands separated in action space. Presumably, the two rounder libration regions with corresponding fixed points at $j\sim 0.75$ and $j\sim 0.4$ are the same libration regions observed in the lower $\epsilon_{\rm sc,z}$ cases (panels a and b), but have significantly shrunk in area because these resonances overlap and they are well-embedded in a chaotic sea. The banana-shaped libration island at higher $j$ is a new resonance and, admittedly, there can be other very narrow resonances that are not captured by our set of initial conditions. Most interestingly, we observe that outside these libration islands, which correspond to a small area of the available phase-space indicated by the red lines, the motion is chaotic. Thus, most initial conditions can diffuse towards low angular momentum: $j\to 0$ or, equivalently, $e\to 1$.

As we increase $\epsilon_{\rm sc,z}$ from 0.2 to 1 (panel d), we observe that chaos is widespread and it reaches all the available phase-space. Thus, all initial conditions reach $e\to 1$ given enough time for the stochastic diffusion to populate the low$-j$ regions of phase-space. We come back to this issue of the diffusion timescale in section §5.3, where we show in Figure 9 that the fraction of low$-j$ orbits depends smoothly on the evolution time.

Finally, we remark that widespread chaos is not a condition to reach $1-e\ll 1$ and quasi-periodic orbits can reach extreme eccentricities simply by crossing $j_{\rm z}=0$, similar to the behavior described in the previous section for small $\epsilon_{\rm sc,z}$. The main difference with the small $\epsilon_{\rm sc,z}$ case is that such crossing can occur for initially larger volume of the phase-space since a rough necessary condition is to have $|j_{\rm z}|\lesssim\epsilon_{\rm sc,z}$.

4.2.3 Limit of large $\epsilon_{\rm sc,z}$:

fast nodal precession

In the limit of large $\epsilon_{\rm sc,z}$ the binary outer orbit precesses around $\hat{\bf n}_{z}$ on a timescale that is much shorter than $\tau_{\mbox{\tiny LK}}$. We can express the potential in Equation (LABEL:eq:phi) in the reference frame of the cluster principal axes

[TABLE]

with

[TABLE]

and average this potential over the fast varying angle $\Omega_{\rm out,z}(t)\propto t$ (i.e., keeping e and ${\bf{j}}$ constant in one precession cycle) as777An equivalent expression is easily found using the Hamiltonian written in terms of orbital elements (see, e.g., Equation 20 in Naoz et al. 2013) after averaging over the fast varying longitude of the ascending nodes.

[TABLE]

Remarkably, this potential coincides with the original one in Equation (38) but multiplied by the constant term $(1-3/2\sin^{2}i_{\rm out,z})$ and replacing $\hat{\bf j}_{\rm out}$ by $\hat{{\bf n}}_{z}$. This potential is axisymmetric with respect to $\hat{\bf n}_{z}$ and, therefore, $\ell_{\rm z}\equiv\sqrt{1-e^{2}}\cos i_{\rm in,z}$ with $i_{\rm in,z}=\cos^{-1}\hat{\bf{j}}\cdot\hat{\bf n}_{z}$ is a constant of motion. Thus, from this potential we can determine the maximum eccentricity, which is reached at $\omega=\pi/2$ and $3\pi/2$ (i.e., $[\hat{\bf e}\cdot\hat{\bf n}_{z}]^{2}=\sin^{2}i$), starting from $e=0$ and inclination $i_{\rm z,0}$ as

[TABLE]

meaning that we recover the classic expression of the maximum eccentricity from the LK mechanism but with the inclination of the inner orbit relative to $\hat{\bf n}_{z}$, instead of $\hat{\bf j}_{\rm out}$.

The precession amplitude, which is determined by $\sin i_{\rm out,z}$, only changes the magnitude of the potential and, therefore, these new LK eccentricity oscillations occur in a timescale

[TABLE]

where $\tau_{\mbox{\tiny LK}}$ is defined in Equation (10). We have checked numerically that this expression for the eccentricity oscillations scales well with different values of $i_{\rm out,z}$. In particular, we checked that the timescale of these oscillations indeed diverges as $\sin^{2}i_{\rm out,z}\to 2/3$.

In the panels c and f of Figure 4 we show an example of the evolution for $\epsilon_{\rm sc,z}=50$ and find that, as explained above, $i_{\rm in,z}$ (red curve) shows regular LK cycles, while the inclination $i_{\rm in,out}$ (black line) oscillates rapidly around the $i_{\rm in,z}$. The maximum eccentricity is given roughly by $e_{\rm max}$ in Equation (41) and we checked that this expression is exact for much larger $\epsilon_{\rm sc,z}$ (not shown). Also, we observe that the period of the oscillation (red curve in panel c) is $\sim 10\tau_{\mbox{\tiny LK}}$, consistent with the longer period due to the extra multiplicative term of $|1-3/2\sin^{2}i_{\rm out,z}|^{-1}\simeq 8.35$ in Equation (42). As expected from the potential in Equation (40) one recovers the standard LK dynamics, but now the orbital elements have to be defined relative to the cluster orientation (see, e.g., the circulating trajectory in $e-\omega_{\rm in,z}$ from panel f).

In summary, for large $\epsilon_{\rm sc,z}$ the cluster induces fast nodal precession so the time-averaged potential of the outer binary acting in the inner one is axi-symmetric around $\hat{\bf n}_{z}$ and, mathematically equivalent to that of the LK mechanism. Thus, this regime can only excite extreme eccentricities for binaries with $i_{\rm in,z}\sim 90^{\circ}$.

4.3. Effect from triaxiality

We briefly discussed the main effects from including a triaxial component to stellar cluster, meaning that $0<\epsilon_{\rm y}<\epsilon_{\rm z}$.

We expect that in the regime of small $\epsilon_{\rm sc,z}$, adding an extra source of slow nodal precession around $\hat{\bf n}_{y}$ will not change the behavior relative to the axisymmetric case and the maximum eccentricity will be given by the classic LK expression. On the contrary, as we will show, for large $\epsilon_{\rm sc,z}$ adding a slower precession around $\hat{\bf n}_{y}$ on top of the dominant precession around $\hat{\bf n}_{z}$ can have important effects. We ignore the complicated regime in which $\epsilon_{\rm sc,y}\lesssim\epsilon_{\rm sc,z}\sim 1$ because there are too many parameters that might change the phase-space structure ($\epsilon_{\rm sc,y}$, $\epsilon_{\rm sc,z}$ and orientations $\Omega_{\rm out,z}$, $i_{\rm out,z}$) and leave it for a future study.

In Figure 6 we show an example in the large $\epsilon_{\rm sc,z}$ regime, similar to that in panels c and f of Figure 4, which we discussed in §4.2.3. The only difference is that we set $\epsilon_{y}=10$, so $\epsilon_{y}/\epsilon_{z}=0.2$. This means that the time evolution of $\hat{\bf j}_{\rm out}$ is described by the following equations (Eq. (A13) with $1-e_{\rm out}\ll 1$):

[TABLE]

which describes the precession of $\hat{\bf j}_{\rm out}$ around both $\hat{\bf n}_{y}$ and $\hat{\bf n}_{z}$. We observe that adding this slower precessional frequency around $\hat{\bf n}_{y}$ creates a extra modulation in the inclinations (panel a) compared to the regular motion in panel c of Figure 4, and the orbit flips every $\sim 4$ oscillations, roughly coincident with the ratio of precession timescales $\epsilon_{\rm z}/\epsilon_{\rm y}=5$. This longer modulation is better observed in panel c where we show the evolution of $\ell_{z}$, which is constant for $\epsilon_{y}=0$, and observe that in this example it can cross $\ell_{z}=0$. During these crossings, the eccentricity reaches extreme values (panel b) and reaches up to $1-e\sim 10^{-7}$ after a time of $\sim 70\tau_{\mbox{\tiny LK}}$. Recall that in this particular example, the oscillations take extra time compared to other choices because of the dilation term $|1-3/2\sin^{2}i_{\rm out,z}|^{-1}\simeq 8.35$ in Equation (42). Other choices of $i_{\rm out}$ can make the timescales of the oscillations comparable to $\tau_{\mbox{\tiny LK}}$.

In order to have an intuitive understanding of the dynamics we give an heuristic argument based on the previous section and the assumption that the level of triaxiality is small, $\epsilon_{\rm y}/\epsilon_{\rm z}\ll 1$. Thus, we can treat this new precession source as a perturbation to the simple axisymmetric potential in Equation (40). In this regime, we expect that averaging the LK potential over the motion of $\hat{\bf j}_{\rm out}(t)$, as we did in §4.2.3, creates a strong axisymmetric potential in the direction of $\hat{\bf n}_{z}$ with a small ($\mathcal{O}[\epsilon_{\rm y}/\epsilon_{\rm z}]$) component in the $\hat{\bf n}_{y}$ direction. If so, this means that the resulting potential should have a terms $\propto\epsilon_{\rm y}/\epsilon_{\rm z}({\bf{j}}\cdot\hat{\bf n}_{y})({\bf{j}}\cdot\hat{\bf n}_{z})$ plus other terms involving ${\bf{e}}$. Thus, $\ell_{z}={\bf{j}}\cdot\hat{\bf n}_{z}$ is no longer fixed and would change as $\dot{\ell}_{z}\propto\epsilon_{\rm y}/\epsilon_{\rm z}\sin i_{\rm out,z}\ell_{z}$, which would explain the longer period oscillations of $\ell_{z}$ observed in panel c.

In summary, adding the triaxial component to the model affects most dramatically the regime of fast nodal precession ($\epsilon_{\rm sc,z}\gg 1$) because $\ell_{z}$ can change slowly in time and cross [math], allowing the eccentricity to reach extremely high values.

5. Population synthesis

We explore the effect from the cluster tidal field for a large population of binaries and how it depends on various parameters. We focus on the statistical distributions of maximum eccentricities ($e_{\rm max}$) because in the secular approximation its evolution mainly determines whether a merger is possible or not.

With the exception of our results in §5.5 we ignore the effect from relativistic precession in our calculations.

5.1. Fiducial case with $\epsilon_{\rm sc,z}=1$

In Figure 7, we show the maximum eccentricities for randomly oriented binaries in a cluster with $\epsilon_{\rm sc,z}=1$ ($\epsilon_{\rm sc}=10$ and $\epsilon_{\rm z}=0.1$ in Eq. [13]). We set $e_{\rm in}=0.1$ initially so the excitation of eccentricities is clearly visible in the plots, while at the same time it avoids extremely small $e_{\rm in}$ that lie near the separatrix, leading to long-timescale LK cycles and possibly to enhanced chaotic behavior.

In panel a, we show how the maximum eccentricities depend on $i_{\rm in,out}$ (angle between $\hat{\bf j}_{\rm in}$ and $\hat{\bf j}_{\rm out}$) and compare these with the expression $e_{\rm max}=\sqrt{1-5/3\cos^{2}i_{\rm in,out,0}}$ applicable in the absence of the cluster potential ($\epsilon_{\rm sc,z}=0$). Probably not surprisingly, we observe that most binaries reach $1-e_{\rm max}\ll 1-\sqrt{1-5/3\cos^{2}i_{\rm in,out,0}}$ (most black dots are well below the solid red line). The distribution of $1-e_{\rm max}$ shown in panel b depicts more clearly this behavior and we observe that the case with $\epsilon_{\rm sc,z}=1$ has a much shallower cumulative distribution towards smaller $1-e_{\rm max}$ compared to the runs with $\epsilon_{\rm sc,z}=0$. In particular, $\sim 50\%$ ($\sim 30\%$) of the binaries in the runs with $\epsilon_{\rm sc,z}=1$ reach $1-e_{\rm max}\lesssim 10^{-3}$ ($1-e_{\rm max}\lesssim 10^{-4}$) compared to $\sim 3\%$ ($\sim 1\%$) in the absence of the cluster potential ($\epsilon_{\rm sc,z}=0$).

In panel c we show the initial inclinations for two samples of our fiducial run: binaries that reach extreme eccentricities $1-e_{\rm max}<10^{-4}$ (green dots) and systems with $1-e_{\rm max}>10^{-2}$ (black dots). We observe that extreme eccentricities are almost always reached when the initial $\cos i_{\rm in,out}\sim 0$, regardless of the value of $i_{\rm out,z}$ (inclination of the outer orbit relative to $\hat{\bf n}_{z}$). This behavior is expected because these systems start with very small $j_{z}$ and from Equation (33) we expect that the forcing due to the cluster almost always allows $j_{z}$ to cross 0 at which point the eccentricity can be highly excited. On the contrary, if the initial $|j_{z}|$ is closer to unity ($|\cos i_{\rm in,out}|\sim 1$), then there is a more limited set of initial conditions that allow for extreme eccentricity excitation and this reflects on the limited range of $\cos i_{\rm out,z}$ that reach $1-e<10^{-4}$, as observed in the figure.

5.2. Effect from the cluster

field strength: dependence on $\epsilon_{\rm sc,z}$

In Figure 8 we explore the effect of varying $\epsilon_{\rm sc,z}$. We do so by varying $\epsilon_{\rm sc}$ and fixing $\epsilon_{\rm y}=0$ and $\epsilon_{\rm z}=0.1$ (see Equation 13). Recall that $\epsilon_{\rm sc,z}=\epsilon_{\rm z}\epsilon_{\rm sc}$ and that small values of $\epsilon_{\rm z}$ ensure that the outer orbit remains nearly circular and only undergoes nodal precession around $\hat{\bf n}_{z}$.

Figure 8 gives the distributions of $e_{\rm max}$ for different values of $\epsilon_{\rm sc,z}$. We observe that our fiducial set-up with $\epsilon_{\rm sc,z}=1$ is the most efficient at exciting extreme eccentricities followed by $\epsilon_{\rm sc,z}=10$ and $\epsilon_{\rm sc,z}=0.1$. This behavior might be expected because when $\epsilon_{\rm sc,z}\sim 1$ the binary is subject to global chaos for a wide range of initial conditions, while for smaller values of $\epsilon_{\rm sc,z}$ the chaotic regions shrinks significantly and might be bounded to a minimum $j$ (see Figure 5). We expect that the chaotic region also shrinks for larger $\epsilon_{\rm sc,z}$ so there is a smooth transition towards regularity that we find for large $\epsilon_{\rm sc,z}$. Regardless of these differences, for all the cases above ($\epsilon_{\rm sc,z}=\{0.1,1,10\}$) the excitation of extreme eccentricities, $1-e_{\rm max}\lesssim 10^{-4}$, occurs in $\sim 10-50\%$ of the binaries, compared to only $\sim 1\%$ in the absence of the cluster potential ($\epsilon_{\rm sc,z}=0$).

On the extreme ends of $\epsilon_{\rm sc,z}$ ($0.01$ or $100$ in Figure 8) the behavior is more easily explained. First, for very small values of $\epsilon_{\rm sc,z}$, only systems with very small values of $j_{\rm z}$ can be appreciably affected by the cluster field. In particular, from Equation (37) we get that the change in $j_{\rm z}$ over one LK cycle is $\sim\epsilon_{\rm sc,z}$, meaning that the initial conditions starting from $|j_{\rm z,0}|\lesssim\epsilon_{\rm sc,z}$ can go through $j_{\rm z}=0$, unless a possible libration of $\Omega$ continues to coherently force $j_{\rm z}$. Since reaching $j_{\rm z}=0$ is a necessary condition to reach extreme eccentricities and the initial distribution of $j_{\rm z}$ is uniform (isotropic inclinations), we expect that roughly a fraction $\sim\epsilon_{\rm sc,z}$ will be strongly affected by the cluster tides. Consistent with these expectations, we observe from Figure 8 that in the $\epsilon_{\rm sc,z}=0.1$ and $\epsilon_{\rm sc,z}=0.01$ cases, $1-e_{\rm max}$ departs from the $\epsilon_{\rm sc,z}=0$ line (no cluster case) and flattens to to the left to reach much higher eccentricities in $\sim 10\%$ and $\sim 1\%$ of the binaries (CDF of $\sim 0.1$ and $\sim 0.01$, respectively).

Second, as discussed in §4.2.3, for very large values of $\epsilon_{\rm sc,z}$ the maximum eccentricity is given by $e_{\rm max}=\sqrt{1-5/3\cos^{2}i_{\rm in,z,0}}$, where $i_{\rm in,z,0}=\cos^{-1}\hat{\bf{j}}\cdot\hat{\bf n}_{z}$ initially (Equation [41]). Since $\hat{\bf{j}}$ is initially chosen to be isotropically distributed, we expect that the distributions of $e_{\rm max}$ for $\epsilon_{\rm sc,z}\gg 1$ and $\epsilon_{\rm sc,z}=0$ coincide. From Figure 8 we observe that the curves $\epsilon_{\rm sc,z}=100$ and $\epsilon_{\rm sc,z}=0$ coincide for $1-e_{\rm max}\gtrsim 10^{-2}$, while for $1-e_{\rm max}\lesssim 10^{-2}$, $\epsilon_{\rm sc,z}=100$ is above the runs with $\epsilon_{\rm sc,z}=0$. This departure at very large eccentricities is due to small eccentricity oscillations taking place in the short timescales $\tau_{\tiny\rm KL}$ when $\epsilon_{\rm sc,z}$ is large, but finite so the approximation to get Equation (41) is not fully fulfilled (see the wiggles around $e_{\rm max}$ in panel f of Figure 4). These oscillations disappear as we increase $\epsilon_{\rm sc,z}$ to values $\gg 100$.

5.3. Effect from maximum evolution time: number of LK

cycles $N_{\mbox{\tiny LK}}^{\rm max}$

In the cluster environment the binaries can be subject to stellar evolution or evaporation before a merger can happen. Thus, it is important to understand how the maximum evolution time affects the merger fractions.

We give the maximum integration time in units of LK cyles as $N_{\mbox{\tiny LK}}^{\rm max}=t_{\rm max}/\tau_{\mbox{\tiny LK}}$, which can be easily converted to absolute timescales using Equation (10) or compared directly to the collision time $N_{\mbox{\tiny LK, coll}}$ and evaporation time $N_{\mbox{\tiny LK, evap}}$ in Equations (LABEL:eq:Ncoll) and (21), respectively.

In Figure 9 we show the distribution of $e_{\rm max}$ for different values of $N_{\mbox{\tiny LK}}^{\rm max}$ from 5 to 300 and observe that the fraction of systems reaching $1-e_{\rm max}\sim 10^{-4}-10^{-5}$ increases gradually with $N_{\mbox{\tiny LK}}^{\rm max}$. This gradual increase is expected from the chaotic diffusion of the eccentricities that survey the allowed parameter space eventually reaching $e\to 1$ (see panel d in Figure 5). Evidently, $e_{\rm max}$ can only increase with $N_{\mbox{\tiny LK}}^{\rm max}$, but we would expect that regular or quasi-periodic orbits saturate this growth after a number of cycles. We do not see such saturation.

This figure shows that $N_{\mbox{\tiny LK}}^{\rm max}$ critically determines the fraction of systems that can merge. Therefore, the binary evaporation or collision times are expected to play a major role at setting the merger rates in real astrophysical systems. On the other hand, we observe that even after $5-25$ cycles (blue and yellow dashed lines) the effect of the cluster becomes significant compared to the no cluster case (red line): after 25 cycles the number of systems reaching $1-e_{\rm max}\sim 10^{-4}$ is boosted by a factor of $\sim 20$.

5.4. Fully triaxial case: $0<\epsilon_{y}\lesssim\epsilon_{z}$

We briefly explore the effect from having a triaxial potential by having a finite value of $\epsilon_{\rm y}$.

In Figure 10 we show the distribution of maximum eccentricities for a set of integrations with $\epsilon_{\rm sc,z}=\{1,10,100\}$ in the axisymmetric potential (solid lines with $\epsilon_{\rm sc,y}=0$) and triaxial potentials with $\epsilon_{\rm sc,z}/\epsilon_{\rm sc,y}=\epsilon_{\rm z}/\epsilon_{\rm y}=5$ (dashed lines). We observe that adding the triaxial component to the axisymmetric potential increases the fraction of systems reaching higher $e_{\rm max}$ for all $\epsilon_{\rm sc,z}$ (all dashed lines are above of their corresponding color-matching solid lines). This behavior might be expected since adding precession around another axis can increase the phase-space volume of the chaotic regions for $\epsilon_{\rm sc,z}=1$ and drive the angular momentum $\ell_{z}={\bf j}\cdot\hat{\bf n}_{z}$ for the large $\epsilon_{\rm sc,z}$ cases (see panel c in Figure 6).

From this figure we also observe that the effect from the triaxial component at exciting large eccentricities is much more dramatic for larger $\epsilon_{\rm sc,z}$. In particular, for $\epsilon_{\rm sc,z}=1$ the triaxial integrations produce a distribution of $e_{\rm max}$ that is only slightly above its corresponding axisymmetric case. On the contrary, for $\epsilon_{\rm sc,z}=100$ we observe that the fraction of systems reaching $1-e_{\rm max}\lesssim 10^{-4}$ ($1-e_{\rm max}\lesssim 10^{-5}$) is $\simeq 25\%$ ($\simeq 10\%$) for the triaxial case compared to only $\simeq 3\%$ ($\simeq 1\%$) in the axisymmetric case. This dramatic effect of the triaxial component for the large $\epsilon_{\rm sc,z}$ was discussed in §4.3 and it is due to the extra slower precession around $\hat{\bf n}_{y}$, which can force $\ell_{z}={\bf j}\cdot\hat{\bf n}_{z}$ to cross [math] time at which the eccentricity $e\to 1$ (see Figure 6).

5.5. Effect from general relativistic precession

It is well known that external sources of apsidal precession such as general relativistic precession can limit the maximum eccentricity growth of the LK oscillations (e.g., Blaes et al. 2002; Liu et al. 2015). We quantify the effect from GR precession using the dimensionless ratio of timescales $\epsilon_{\rm GR}=\tau_{\mbox{\tiny LK}}/\tau_{\rm GR}$ in Equation (26).

In the absence of the cluster perturbations we can determine the maximum eccentricity reached in a LK cycle starting from low eccentricity and $i_{\rm in,out}=90^{\circ}$ as (e.g., Eq. 35 of Fabrycky & Tremaine 2007):

[TABLE]

We incorporate the first post-Newtonian corrections in our model by adding the following term to the Equation (A10):

[TABLE]

In Figure 11 we show the distribution of maximum eccentricities for different values of $\epsilon_{\rm GR}$, where $\epsilon_{\rm GR}=0$ corresponds to our fiducial simulation. As expected, the larger $\epsilon_{\rm GR}$ (faster GR precession) the smaller $e_{\rm max}$ from the ensemble of binaries. For comparison we show the corresponding values of $e_{\rm max,GR}$ from Equation (44) for the LK cycles with no cluster. We observe that the distributions show a clear drop at a maximum eccentricity and that this maximum roughly coincides with $e_{\rm max,GR}$. For instance, if we set $\epsilon_{\rm GR}=0.01$ (0.1) then $1-e_{\rm max,GR}\simeq 2.2\times 10^{-5}$ ($\simeq 2\times 10^{-3}$), while the maximum value in the ensemble of binaries is $1-e_{\rm max}\simeq 10^{-5}$ ($\simeq 0.8\times 10^{-3}$). Thus, from these examples it seems that $1-e_{\rm max,GR}$ is a factor of $\sim 2$ larger than the numerical integrations including the cluster potential.

In conclusion, relativistic precession quenches the eccentricity growth up to a maximum value that is only slightly higher ($1-e_{\rm max}$ smaller by a factor of $\sim 2$) than that derived for the LK mechanism in the absence of the cluster perturbations ($e_{\rm max,GR}$ in Equation [44]). In order for BH (stellar) binaries to merge we require that $\epsilon_{\rm GR}\lesssim 0.01$ ($\lesssim 0.1$) so $1-e_{\rm max}\lesssim 10^{-5}$ ($1-e_{\rm max}\lesssim 10^{-3}$).

5.6. Summary of the dynamical results

Our main result from the dynamical evolution of binaries is that when the nodal precession timescale of the binary’s center of mass due to stellar cluster becomes comparable (within a factor of $\sim 10$) to the torquing timescale of the binary orbit by to the MBH (Lidov-Kozai timescale), the binary often ($\sim 10-50\%$ of the cases) reaches extreme eccentricities ($1-e\lesssim 10^{-4}$). In the absence of the non-spherical cluster potential, this rate decreases by at least an order of magnitude. Thus, we refer to this mechanism as the cluster-enhanced Lidov-Kozai mechanism.

We study the dynamical behavior for different values of $\epsilon_{\rm sc,z}$ (Eq. [13]; ratio between the LK timescales and the nodal precession due to the cluster) and identify the following regimes that result in extreme eccentricity excitation:

• •

for small $\epsilon_{\rm sc,z}$ (slow nodal precession), the dynamics resembles the octupole-level perturbations in the standard LK mechanism as the main effect of the cluster is to change the, otherwise conserved, angular momentum in the $z-$direction $j_{\rm z}$, which can drive the binary to cross $j_{\rm z}=0$ and lead to $e\to 1$ (see Equation [33]).

• •

for intermediate values of $\epsilon_{\rm sc,z}$, the eccentricities and inclinations often undergo widespread chaotic evolution as the result of overlapping resonances and the eccentricities can diffuse towards $e\to 1$ (see surfaces of section in Figure 5);

• •

for large $\epsilon_{\rm sc,z}$, the binary orbit around the MBH undergoes fast nodal precession around $\hat{\bf n}_{\rm z}$ creating an effective axisymmetric potential around this axis. This potential coincides with that from the triple interactions (LK potential), but defining the orbital elements relative to $\hat{\bf n}_{\rm z}$ instead of $\hat{\bf j}_{\rm out}$ and changing its amplitude (i.e, the timescale) by a constant factor $1-3/2\sin^{2}i_{\rm out,z}$. Thus, the binary undergoes LK cycles with a maximum eccentricity given by the angle between ${\bf j}$ and $\hat{\bf n}_{\rm z}$ (Eq. [41]). This channel is expected to contribute most significantly to the extreme eccentricity excitation for fully triaxial clusters for which $\ell_{\rm z}={\bf j}\cdot\hat{\bf n}_{\rm z}$ can cross zero due to an extra slower precession source around $\hat{\bf n}_{\rm y}$ (see Figure 6).

Other findings include:

• •

the dynamics for an axisymmetric cluster with small $\epsilon_{z}$ can be well-described by a simple two-degree-of-freedom Hamiltonian (Equation [32]). The surfaces of section of this Hamiltonian show that the LK resonance overlaps with other resonances arising from the cluster-driven nodal precession as $\epsilon_{\rm sc,z}$ approaches near unity values (see Figure 5).

• •

the distribution of maximum eccentricities depends critically on the maximum evolution timescale, possibly linked to a chaotic diffusion of the eccentricities (see Figure 9);

• •

the eccentricity growth is quenched by short-range forces similar to the LK mechanism. In particular, when GR precession is included the minimum periapsis distance $a_{\rm in}(1-e_{\rm max})$ in our model is typically a factor of $\sim 2$ below that from the standard LK mechanism ($1-e_{\rm max,GR}$ from Equation [44]).

We note that, although different in nature, the behavior described above resembles the dynamics of four-body systems, where a fourth body (the nuclear cluster in our case) can change $j_{z}$ in a triple system (Pejcha et al. 2013; Hamers et al. 2015). Recently, Hamers & Lai (2017) have studied in more detail the dynamics of four-body systems and have independently arrived to similar conclusions than ours regarding the general dynamical behavior of triple systems affected by the nodal precession of a tertiary.

Our dynamical results are expected to have important effects on the mergers rates of stellar binaries in general. In what follows, we quantify the relevance our our mechanism for the compact-object (neutron stars and black holes) binaries and differ the study of other combinations of binary systems for a future study.

6. Merger of compact-object binaries and

implications for gravitational wave detections

In this section we compute the number of compact-object binary mergers that are formed during a star formation episode similar to the one that gave birth to the young massive stars observed in the Galactic center. We consider BH-BH, NS-NS, and BH-NS binaries. We use the stellar population synthesis code BinaryStarEvolution (BSE; Hurley et al. 2002) to evolve a population of massive binaries stars and compute the properties of compact-object binaries formed from this population. The compact-object binaries are then evolved in a model of the Galactic center and the dependence of the fraction of merging binaries on the degree of flattening of the NC is quantified. The results of the simulations are then used as a set of baselines for making predictions about the binary merger rates for ground-based laser interferometers such as Advanced LIGO and VIRGO.

We assume in what follows that the binary population in galactic nuclei follows the population in the field. Although little is known about the binary properties in the Galactic center there are observations suggesting that our assumption is reasonable. Long-term spectroscopic and photometric surveys of the most luminous and massive stars in the vicinity of the massive black hole reveal the existence of three OB/WR binaries in the inner 0.2 pc of the Galactic center. Using radial velocity change upper limits, Pfuhl et al. (2014) find that the spectroscopic binary fraction in the Galactic center is $\sim 0.3$, close to that observed in nearby young star clusters (Sana et al. 2012). These authors also find that the fraction of eclipsing binaries in the central parsec, $\sim 3\%$, is consistent with the fraction of such binaries in local OB star clusters (Lefèvre et al. 2009). Overall, observations suggest that the Galactic center binary fraction seems to be similar to the binary fraction in nearby young clusters.

In all our models, we assumed Solar metallicity and sampled the mass of the most massive star in the binary from the initial mass function (IMF) $dN/dm\propto m^{-\alpha}$, with $\alpha=1.7$ as observed in the galactic center (Lu et al. 2013), or with the canonical value for field stars of $\alpha=2.3$ (Kroupa 2002). We adopted a mass ratio $m_{2}/m_{1}$ distribution that is uniform between 0 and 1. This choice is consistent with the observed mass ratios of massive binary stars (Duchêne & Kraus 2013; Sana et al. 2012; Kobulnicky et al. 2014). Stellar masses were sampled in the range $10$ to $100\ M_{\odot}$. The upper limit on the mass comes from the fact that the stellar evolution tracks used in BSE are not valid above $100\ M_{\odot}$.

The distribution of orbital periods was $dN\propto(\log P_{\rm in}/\rm days)^{-0.55}$ in the range $(0.15-4)$, while the orbital eccentricities were sampled from a thermal distribution, $dN\propto e_{\rm in}$. These distributions are both consistent with observations of Galactic massive stars in nearby young clusters (Duquennoy & Mayor 1991; Sana et al. 2012).

We evolved the stellar binaries until two compact-objects were formed. When a compact-object was formed we assigned a natal kick from a Maxwellian distribution with dispersion $v_{\rm K}=250\rm km\ s^{-1}$, in agreement with observational constraints (Hobbs et al. 2005). If the remnant was a BH, the kick velocity was normalized to the mass of the remnant, assuming linear momentum conservation. Stars with zero-age main-sequence mass above $40M_{\odot}$ do not receive any natal kick as in Fryer & Kalogera (2001).

The distribution of the binary orbits around the MBH we adopted is also motivated by observations. Spectroscopic studies of the Galactic center reveal a population of $\approx 200$ very massive early-type stars, including Wolf-Rayet stars and O and B type main-sequence stars, giants, and supergiants (e.g., Bartko et al. 2010; Yelda et al. 2014). These massive blue giant stars will probably end their lives as an equal number of BHs and NSs. The majority of these massive stars are found at distances of $0.04$ to $0.4\rm pc$ from the MBH, and follow the surface density distribution $\Sigma\propto R^{-1.5}$ (Bartko et al. 2010). About $\approx 20\%$ of them reside in a stellar disk exhibiting a clockwise motion pattern on the sky (Yelda et al. 2014). We assume that the distribution of the compact-object binary orbits follows the distribution of their stellar progenitors. This is likely a good approximation for systems that remain bound after their components have received birth kicks (Bortolas et al. 2017). Accordingly, we sample the semi-major axis of the binary external orbits in the range $0.04\leq a_{\rm out}\leq 0.4\ \rm pc$ and from the space density distribution $\rho(r)\propto r^{-2.5}$.

To simplify the calculation, we suppose that the orbits of the young stars randomize in orientation and eccentricity on a timescale shorter than the relevant evolutionary timescales. Several mechanisms might achieve this: e.g., resonant-relaxation, torques from another disk. Moreover, as also mentioned above, the majority of the young stars do not reside inside the disk, but are more isotropically distributed around the center (Yelda et al. 2014). Accordingly, we assume an isotropic distribution of inclinations, and a thermal distribution of eccentricities, i.e., $dN\propto e_{\rm out}$. Our results are somewhat insensitive to this latter choice and we checked that by setting $e_{\rm out}=0$ we get similar results for the merger fraction. The longitude of the ascending node, and the argument of periapsis of both inner and outer orbit were also randomly distributed between [math] and $2\pi$.

We set $M_{\rm BH}=4\times 10^{6}\ M_{\odot}$, and assume that the field stars follow the density distribution in Equation (11), consistent with recent observational constraints (Fritz et al. 2016; Schödel et al. 2017). The initial conditions were evolved until a maximum timescale $\tau_{\rm max}$, which we set equal to the evaporation timescale $\tau_{\rm evap}$ in Equation (20). As discussed in §3.2.1, the binary might significantly change its orbital elements before evaporating, so we also provide the results for the much more conservative case of setting $\tau_{\rm max}$ equal to the collision time $\tau_{\rm coll}$ in Equation (17). We stop the integration if the binary periapsis distance to the MBH satisfied $a_{\rm out}(1-e_{\rm out})<r_{\rm bt}$, with $r_{\rm bt}=a_{\rm in}\left[M_{\rm BH}/(m_{1}+m_{2})\right]^{1/3}$ the binary tidal break-up radius. During the integration, we checked whether $\tau_{\rm LK}\sqrt{1-e_{\rm in}^{2}}\leq\tau_{\rm GW}$ (see Eq. [25]). If this condition is satisfied the gravitational merger will take place within one LK cycle and the binary was considered to have merged. We included the effect from relativistic precession using Equation (45). We considered four models characterized by different values of $\epsilon_{z}=(0.001,0.01,0.1,1)$. In all models we set $\epsilon_{y}=0$.

In total we evolved $10^{6}$ massive binaries. The fraction, $f_{2}$, of these binaries that produced stable BH-BH, NS-NS, and BH-NS binaries are given in Table 1. Only a small fraction of massive binary stars produce a bound pair of compact objects; most binaries are either disrupted due to natal kicks or merge during their main-sequence evolution.

Table 1 summarizes the results of our calculations, showing the fraction of formed BH binaries that merge, $f_{\rm m}$, as a function of the parameter $\epsilon_{z}$ that measures flattening of the cluster potential (see Equation [7] for its relation to cluster ellipticity $\mathcal{E}$ for $r\sim r_{0}$). When computing $f_{\rm m}$ we did not include binaries that underwent a phase of common-envelope evolution, since these binaries have $\epsilon_{\rm GR}>1$ and the LK process is strongly quenched. Our rate estimates below are therefore only for mergers due to the cluster-enhanced LK mechanism, and do not include those binaries that would merge as a result o stellar evolution processes.

The results shown in Table 1 exemplify the importance of the mechanism described in this paper, by showing a strong dependence of the merging fraction on the NC morphology. For BH-BH binaries we show that nearly spherical models ($\epsilon_{z}=0.001$ or ellipticity $\mathcal{E}\simeq 0.001$) the merger fraction is only $\approx 10^{-3}$. We expect therefore a small merger rate of BH binaries in spherical NCs, in agreement with previous work (Antonini & Perets 2012). In contrast, for $\epsilon_{z}=0.1$ ($\mathcal{E}\simeq 0.1$) the merging fraction increases up to $\approx 6\%$, and for $\epsilon_{z}=1$ (fully flattened, $\mathcal{E}=1$) to $\approx 17\%$. For BH-NS and NS-NS binary mergers we observe that no mergers occur for the spherical clusters, while the fraction reach $f_{\rm m}\sim 1\%$ for flattened clusters.

We remark that we start evolving the orbits of the compact-object binaries under the gravitational influence from the MBH and the cluster once they form. One potential concern from this simplification is that binaries might merge during the main sequence before they become compact objects (e.g., Antonini et al. 2017). We have checked, however, that most BH-BH mergers, $77\%$ and $75\%$ for $\epsilon_{z}=0.1$ and 1 respectively, occurred after a time that was longer than the main sequence lifetime of both binary components, meaning that these might have avoided an unwanted merger prior to the formation of the two BHs. Having relatively late mergers (after many LK timescales and continuous in time) is expected from the chaotic diffusion of the eccentricities (see Figure 9).

We conclude that allowing for a more realistic model in which the shape of the NC is not perfectly spherical will significantly increase the merger rates of compact-object binaries, which for BH-BH binaries it does so by factor of $\sim 30-100$ for $\epsilon_{z}\gtrsim 0.1$ (ellipticities $\gtrsim 0.1$).

6.1. Merger rate estimates

The total mass of the young stars currently observed in the inner $\approx 1\rm pc$ region is a few times $10^{4}\ M_{\odot}$ (Lu et al. 2013). The starburst that created this population of massive stars may have been just the most recent episode of an ongoing process of continuous star formation. Observations suggest in fact that a large fraction of the Milky Way NC was built up gradually in time through in-situ star formation episodes similar to the one from which the observed massive stars were formed (Figer et al. 2004; Pfuhl et al. 2011). Young stellar populations are ubiquitous in NCs that are found in external galaxies (e.g., Carson et al. 2015), suggesting continuous/episodic star formation as a generic mode for NC formation and growth (Rossa et al. 2006; Antonini et al. 2015a, b).

Observations of the central regions of the Milky Way show that this region is dominated by a dense NC consisting of a population of old stars with a total mass of $\sim 3\times 10^{7}M_{\odot}$ (Schödel et al. 2014, 2017; Feldmeier-Krause et al. 2017) with most young stars residing inside $\sim 0.5$ pc (Schödel et al. 2014, 2017; Feldmeier-Krause et al. 2017). We assume here that a fraction of this mass was formed via continuous in-situ star formation throughout the age of the Galaxy. Under this assumption we compute the rate at which compact objects must form, $\Gamma_{\rm co}$, such that the total mass formed in stars and stellar remnants after 10 Gyr is $10^{7}M_{\odot}$, approximately the total stellar mass within the sphere of influence of SgrA*. The value of $\Gamma_{\rm co}$ depends on the assumed IMF and on the star formation history of the NC. Table 1 in Löckmann et al. (2010) gives the composition of an evolved stellar population for different star formation histories and IMFs; we use this table to compute $\Gamma_{\rm co}$.

For a top-heavy IMF similar to the one adopted here and a constant star formation, Löckmann et al. (2010) find that $10^{5}$ BHs and NSs are formed for every $10^{6}M_{\odot}$ in stars and stellar remnants. If we require that $10^{7}M_{\odot}$ in stars and stellar remnants are formed after 10 Gyr, this gives $\Gamma_{\rm co}\sim 10^{6}/10\,\rm{Gyr}=10^{-4}\rm yr^{-1}$. In other words, the predicted total number of compact objects formed in the NC after $10$ Gyr in this case is $\sim 10^{6}$. For a canonical IMF and a constant star formation, Löckmann et al. (2010) find that $1.7\times 10^{4}$ BHs and NSs are formed for every $10^{6}M_{\odot}$ in stars and stellar remnants, which gives $\Gamma_{\rm co}\sim 2\times 10^{-5}\rm yr^{-1}$. Thus, in this latter case the predicted total number of compact objects formed in the NC after 10 Gyr is $\sim 2\times 10^{5}$. The number of compact objects predicted to form by a top-heavy IMF is therefore approximately 10 times larger than for a canonical IMF, given an equal total mass formed after 10 Gyr.

Given the rate of formation of compact objects in the central cluster, we can compute the merger rate of compact-object binaries as:

[TABLE]

where we have assumed that all massive stars form in binaries, similarly to what found from observation of nearby young clusters (Sana et al. 2012).

Next, we estimate the merger rate for compact-object binaries. In order to do so we further assume that all NCs are characterized by some degree of spherical asymmetry (ellipticities $\gtrsim 0.1$ or axis ratios $q\lesssim 0.9$). This assumption appears to be reasonable. The best-fitting Schwarzschild model of the Galactic center in Feldmeier-Krause et al. (2017) has an axis ratio $q\sim 0.3$ in the nuclear region, increasing up to $\sim 0.8$ at $\sim 1-2$ pc. Similarly, Fritz et al. (2016) measure a mean axis ratio of $q\sim 0.8$. Schödel et al. (2014) find $q\sim 0.7$ from Spitzer/IRAC photometry. NCs in external galaxies are also often observed to be strongly flattened, despite of projection effects that can make the clusters look more spherical. The nuclei of the nearby galaxies M31, M32, and M33 all have substantial ellipticities $\gtrsim 0.2$ (Lauer et al. 1998), while the majority of the NCs in the galaxy samples of Georgiev & Böker (2014) and of Carson et al. (2015) have measured ellipticities. Typical values are of order $\approx 0.2-0.3$ but values up to $\approx 0.9$ are also found.

Finally, we stress that the analysis presented here makes several simplifying assumptions which will be relaxed and detailed in future work. Here, we simply note that the merger rates we derived should be viewed as an order-of-magnitude estimate indicating that the cluster-enhanced LK mechanism is an important new merger channel for compact-object binaries, possibly competing with traditional scenarios which invoke either isolated binary evolution, or dynamical formation in globular clusters.

6.2. BH-BH binaries

BH-BH are the most common form of compact-object binaries produced in our models.

In panel a of Figure 12, we show the semi-major axes $a_{\rm in}$ and $a_{\rm out}$ of the BH-BH binaries that merge in our simulations and also show the relevant timescale ratios from panel b of Figure 2. We note that mergers happen preferentially for $\epsilon_{\rm sc,z}\sim 1$, which skews the distribution of $a_{\rm out}$ to wider separations: compare initial distribution (black dashed line) on panel b with those that merge (gray line). We also observe that the merged systems are expected to undergo changes in their orbital parameters in secular timescales only (green line in panel a) and the dynamics is, therefore, properly captures by our equations of motion.

Taking the merging fraction, $f_{\rm m}$, from Table 1, we obtain from Equation (46) a merger rate in the range $\Gamma_{\rm m}\approx(0.03,\ 0.8)\rm\ Myr^{-1}$, where the lower limit corresponds to $\epsilon_{z}=0.1$, and a canonical IMF; the upper limit is for $\epsilon_{z}=1$, and a top-heavy IMF. This merger rate can be compared to that corresponding to standard binary formation models. These latter models predict that the evolution of isolated binary stars in the Galactic field results in a merger rate in the range $0-30\rm\ Myr^{-1}$, assuming Solar metallicity (Dominik et al. 2012). We conclude that our new mechanism gives a local merger rate for BH-BH binaries that is comparable or even larger to that predicted by other standard scenarios (e.g., Abbott et al. 2016).

In order to estimate a rate per volume we use a number density of galaxies of $\approx 0.02\rm\ Mpc^{-3}$ (e.g., Conselice et al. 2005; Kopparapu et al. 2008). This gives a merger rate in the range

[TABLE]

Realistic values are likely to be somewhere in the middle of this range. For instance, by taking $\epsilon_{z}=0.1$ we have $\Gamma_{\rm BH-BH}\approx 5\rm\ Gpc^{-3}yr^{-1}$ for a top-heavy initial mass function. If we take $\epsilon_{z}=0.3$, we have $\Gamma_{\rm BH-BH}\approx 1\rm\ Gpc^{-3}yr^{-1}$ for a canonical initial mass function.

Remarkably, the merger rate we derived is similar to the predicted local merger rate of BHs formed in globular clusters $\Gamma=(2,\ 20)\rm\ Gpc^{-3}yr^{-1}$ (Rodriguez et al. 2016). This result can appear surprising: according to our models, the cluster-enhanced LK mechanism will realistically only produce about one thousand BH mergers over 10 Gyr. Over the same timescale one globular cluster alone will produce $\approx 100$ BH mergers, resulting in a total of $\approx 2\times 10^{4}$ mergers produced by $200$ globular clusters in the Galaxy (C. Rodriguez; private communication). Part of the reason why the globular cluster rate is relatively low is because a large fraction of the mergers happen early on in the evolution at high redshift, i.e., the merger rate decreases considerably with time (see Figure 12 in Rodriguez et al. 2016). Our more realistic assumptions result in $\lesssim 1/10\rm th$ of the BH mergers produced in globular clusters, however, these are produced at a $constant$ rate giving a relatively higher local merger rate.

Finally, we note that Antonini & Rasio (2016) derived a merger rate in nuclear star clusters without central massive black holes of $\sim 1\rm\ Gpc^{-3}yr^{-1}$, where the BH binaries form dynamically as in GCs. Our estimates show that even when a MBH is present in the NC, BH-BH binaries can still form and merge at a high rate.

6.3. BH-NS mergers

About $1\%$ of the massive binaries we evolved produced a BH-NS binary. Taking the merging fraction, $f_{\rm m}$, from Table 1, we obtain from Equation (46) a Galactic merger rate in the range $\Gamma_{\rm m}\approx(10^{-3},\ 2\times 10^{-2})\rm\ Myr^{-1}$, where, as before, the lower limit corresponds to $\epsilon_{z}=0.1$, and a canonical IMF, while the upper limit is for $\epsilon_{z}=1$, and a top-heavy IMF. For comparison, the standard binary formation channel produces a Galactic merger rate of BH-NS binaries in the range $0-10\rm\ Myr^{-1}$, assuming Solar metallicity (Dominik et al. 2012). Dynamical formation in globular clusters is unlikely to result in an important number of BH-NS binaries, as NSs do not efficiently segregate to the center, and most of them are ejected at birth from the parent clusters due to their high natal kicks.

The corresponding merger rate for BH-NS binaries per unit volume due to the cluster-enhanced LK mechanism is in the range

[TABLE]

6.4. NS-NS mergers

Only about $0.3\%$ of massive binaries in our models form a NS-NS binary. The merging fractions given in Table 1, lead to a Galactic merger rate in the range $\Gamma_{\rm m}\approx(3\times 10^{-5},\ 10^{-3})\rm\ Myr^{-1}$, where, as before, the lower limit corresponds to $\epsilon_{z}=0.1$, and a canonical IMF, and the upper limit is for $\epsilon_{z}=1$, and a top-heavy IMF. The resulting merger rate per unit volume for NS-NS binaries predicted by our models is $\Gamma_{\rm NS-NS}\lesssim 0.02\rm\ Gpc^{-3}yr^{-1}$, and it is therefore negligible compared with the predicted merger rates from the evolution of isolated binaries in the Galactic field (e.g., Dominik et al. 2012).

We remark that in our calculation we have not included mergers favored by stellar evolution processes that can occur prior to the formation of the two NSs. Evolution driven by mass transfer (e.g., Dosopoulou & Kalogera 2016a, b), followed by a common envelope phase (e.g., Kalogera & Webbink 1998) can lead to the formation of short period NS-NS binaries which can merge in less than one Hubble time (e.g., Kalogera et al. 2001; de Mink & Belczynski 2015). For NS-NS binaries, we expect these latter type of mergers to be more frequent than those due to the cluster-enhanced LK mechanism.

7. Conclusions

Galactic centers are often observed to host a central massive black hole embedded by a non-spherical nuclear star cluster. These clusters are composed by a wide range of stellar populations, including very massive stars.

We have studied the orbital evolution of stellar binaries in these environments, focusing on the excitation of large eccentricities that can lead to mergers of compact-object binaries by gravitational radiation.

We show that the excitation of extremely large eccentricities depends critically on the ratio between the Lidov-Kozai (LK) timescale due to the MBH and the timescale of the nodal precession of the binary around the MBH due to an axisymmetric (or triaxial) star cluster, which we define as $\epsilon_{\rm sc,z}$ (Eq. [13]).

For $\epsilon_{\rm sc,z}\sim 0.1-10$, reached by a wide range of binaries in the galactic center, we find that $\sim 10-50\%$ of randomly-oriented binaries reach $1-e\sim 10^{-4}-10^{-6}$, compared to only $\sim 0.1-1\%$ in the idealized case of an spherical star cluster ($\epsilon_{\rm sc,z}=0$). Thus, we identify the cluster-enhanced Lidov-Kozai mechanism as a efficient new channel to drive binary mergers in galactic nuclei.

We study the dynamics for different values of $\epsilon_{\rm sc,z}$ and identify two main channels for the extreme eccentricity excitation: (i) chaotic diffusion of the eccentricities due to the overlap of the LK resonance and resonances arising from the cluster-induced nodal precession; (ii) cluster-driven variations of the mutual inclinations between the binary orbit and its center of mass orbit around the MBH, which drive the otherwise conserved LK constant of motion $j_{z}$ to zero.

We have carried out a population synthesis study to estimate the merger rate of compact-object binaries in NCs hosting a central MBH. Under the observationally motivated assumption that most nuclear star clusters have finite ellipticities ($\gtrsim 0.1$), we derived the following merger rates: (i) BH-BH binaries in the range $(0.6,\ 15)\rm\ Gpc^{-3}yr^{-1}$; (ii) BH-NS binaries in the range $(2\times 10^{-2},\ 0.4)\rm\ Gpc^{-3}yr^{-1}$. These merger rates are comparable to those from either isolated binary evolution or dynamical formation in globular clusters. The wide range of predicted rates is explained mainly by the uncertainties in the star formation history of the galactic nuclei.

Our results suggest a strong correlation between the merger rates of black hole binaries in galactic nuclei and the morphology of their central clusters.

We thank the referee Ben Bar-Or for providing with a careful and constructive report that has helped us to improve the presentation of our paper. We thank Adrian Hamers for sharing his related work on quadruples and commenting on an early version our paper and Scott Tremaine for providing useful feedback on the manuscript. We are grateful to Chris Matzner, Dan Tamayo, Diego Muñoz, Fred Rasio, Juna Kollmeier, Katerina Chatziioannou, Norm Murray, Renu Malhotra, Yanqin Wu, and Aaron Zimmerman for stimulating discussions and comments. CP acknowledges support from the Gruber Foundation. FA acknowledges support from a CIERA postdoctoral fellowship at Northwestern University.

Appendix A Equations of motion for a binary embedded in a nuclear star cluster with

a central massive black hole

We derive the orbit-averaged equations of a motion from the potential per unit mass in Equation (2):

[TABLE]

We start by averaging over the outer orbit the terms that appear on the triaxial potential. First, as shown in Ivanov et al. (2005) and Merritt & Vasiliev (2011) the averaging of the spherical component results in

[TABLE]

where ${}_{2}F_{1}$ is the ordinary hypergeometric function and the constants terms can be approximated as (e.g., Merritt & Vasiliev 2011):

[TABLE]

for $0\leq\gamma<2$.

Second, the non-spherical in the potential involves the averaging of the following terms:

[TABLE]

where $\phi$ is the azimuthal angle in the ${\bf\hat{e}}_{\rm out}-{\bf\hat{q}}_{\rm out}$ plane indicating the position of the bodies along their Keplerian orbits. In principle, we can average the above expression over one period of this phase angle $\phi$ arbitrary value of $\gamma$ (see e.g., Sridhar & Touma 1999; Tremaine & Yavetz 2014 for the averaging method), but the mathematical expressions are complicated and not very illuminating. Instead, we provide the explicit expression for $\gamma=0$ (uniform ellipsoid) and $\gamma=1$ in which cases the algebra is greatly simplified and they capture the dynamical behavior for the general case. We note that the previous study by Merritt & Vasiliev (2011), based on the previous work by Sridhar & Touma (1999) and Sambhus & Sridhar (2000) using Delaunay variables, modeled the cluster potential as the sum of an spherical cusp for an arbitrary $\gamma$ and a triaxial uniform ellipsoid ($\gamma=0$), which is somewhat more general at the expense of having to define an extra density profile for the ellipsoid. Our approach and that in Merritt & Vasiliev (2011) should give similar results.

For $\gamma=0$ the averaging of Equation (A4) results in:

[TABLE]

implying that the averaged potential due to the cluster becomes

[TABLE]

For $\gamma=1$ the averaging of Equation (A4) results in:

[TABLE]

implying that the averaged potential due to the cluster becomes

[TABLE]

In what follows, we present the results for $\gamma=1$ only, while those for $\gamma=0$ follow trivially. Following a similar procedure for the other term in the potential we can then write the averaged potential over the inner and outer orbits:

[TABLE]

where $\mu_{\rm in}=m_{1}m_{2}/(m_{1}+m_{2})$ and $\mu_{\rm out}=(m_{1}+m_{2})M_{\rm BH}/(m_{1}+m_{2}+M_{\rm BH})$ are the reduced masses of the inner and outer orbits, respectively. In order to simplify the dynamics we have ignored the effect from the cluster on the inner orbit, which is a reasonable approximation only well-inside the sphere of influence of the MBH. In particular, for $\gamma=0$ the ratio between the amplitude of potential due to the cluster and that from the MBH acting in the inner orbit is $M_{s}(a_{\mathrm{out}})/(3M_{\rm BH})\ll 1$, where $M_{s}(a_{\mathrm{out}})$ is the enclosed mass of stars in the cluster in a sphere of radius $a_{\mathrm{out}}$. Incorporating the potential of the cluster acting in the inner orbit is straightforward and it simply follows from exchanging the outer orbit by the inner one in the Equations (A6) and (A8).

The first term in the potential corresponds to the well-known Lidov-Kozai interaction potential only up to quadrupole level. We have avoided writing the octupole term in the potential to keep the equations simpler, but they are included in our code and can be found elsewhere in the vectorial formalism (e.g., Boué & Fabrycky 2014; Liu et al. 2015; Petrovich 2015; Petrovich & Muñoz 2017). We note, however, that in the regime that we explore in our paper with $a_{\rm in}\sim 1-10$ AU and $a_{\rm out}\sim 0.1-1$ pc the octupole term has a negligible effect on the dynamics because (i) $\epsilon_{\rm oct}$ in Equation (1) is generally too small ($10^{-5}-10^{-3}$), while mostly values $\gtrsim 0.01$ are likely to contribute significantly to the dynamics (Li et al. 2014); (ii) the spherical cluster induces fast apsidal precession so the octupole contribution, which depends on the apsidal orientation of the outer orbit ($\propto{\bf e}_{\rm in}\cdot{\bf e}_{\rm out}$), averages out to zero (see Eq. [14] for reference).

Finally, the equations of motion can then be conveniently written in a dimensionless form as (e.g., Milankovitch 1939; Tremaine et al. 2009):

[TABLE]

where the unit time is given in Lidov-Kozai timescales, $\tau=t/\tau_{\mbox{\tiny LK}}$, with

[TABLE]

and we have defined $\epsilon_{\rm sc}\equiv\tau_{\mbox{\tiny LK}}/\tau_{\rm sc}$ with

[TABLE]

where $\rho_{s}(a_{\rm out})\equiv\rho_{s}(|{\bf r}|=a_{\rm out})$ with $\epsilon_{z}=\epsilon_{y}=0$ or the density for spherical component. Thus, the ratio of timescales $\epsilon_{\rm sc}$ (sc stands for stellar cluster), assuming $\epsilon_{z},\epsilon_{y}\ll 1$ and $m_{1}+m_{2}\ll M_{\rm BH}$, becomes

[TABLE]

where we have defined the enclosed mass with $\epsilon_{z}=\epsilon_{y}=0$ as:

[TABLE]

For reference, we can take the density profile in Equation (11) with $\gamma=1.3$ for the Milky way’s galactic center and set $M_{\rm BH}=4\times 10^{6}M_{\odot}$ to get

[TABLE]

Using a general density profile changes the scaling to $\epsilon_{\rm sc}\propto a_{\rm out}^{9/2-\gamma}$ and introduces a multiplicative factor of order unity that depends on $\gamma$.