Secular dynamics of binaries in stellar clusters I: general formulation and dependence on cluster potential

TL;DR

This paper develops a general theoretical framework for understanding the long-term orbital evolution of binary systems in axisymmetric stellar clusters, accounting for various potential shapes and orbital configurations.

Contribution

It derives a universal, quadrupole-order secular Hamiltonian for binaries in axisymmetric potentials, extending previous models and unifying different dynamical effects.

Findings

Reproduces known secular effects like Lidov-Kozai and Galactic tides.

Introduces parameters A and Gamma that encapsulate cluster potential effects.

Provides a basis for analyzing binary evolution in diverse stellar environments.

Abstract

Orbital evolution of binary systems in dense stellar clusters is important in a variety of contexts: origin of blue stragglers, progenitors of compact object mergers, millisecond pulsars, and so on. Here we consider the general problem of secular evolution of the orbital elements of a binary system driven by the smooth tidal field of an axisymmetric stellar cluster (globular, nuclear, etc.) in which the binary orbits. We derive a secular Hamiltonian (averaged over both the inner Keplerian orbit of the binary and its outer orbit within the cluster) valid to quadrupole order for an arbitrary cluster potential and explore its characteristics. This doubly-averaged 'tidal' Hamiltonian depends on just two parameters, which fully absorb the information about the background cluster potential and the binary's orbit within it: a dimensional parameter $A$ setting the secular timescale, and a dimensionless parameter $\Gamma$ which determines the phase portrait of the binary's inner orbital evolution. We examine the dependence of $A$ and $\Gamma$ on cluster potential (both spherical and axisymmetric) and on the binary orbit within the cluster. Our theory reproduces known secular results - such as Lidov-Kozai evolution and the effect of the Galactic tide on Oort Cloud comets - in appropriate limits, but is more general. It provides a universal framework for understanding dynamical evolution of various types of binaries driven by the smooth tidal field of any axisymmetric potential. In a companion paper (Hamilton & Rafikov 2019b) we provide a detailed exploration of the resulting orbital dynamics.

Figures

Tables

Table 1. Table 1: Details of the outer orbits 𝐑 g subscript 𝐑 g \mathbf{R}_{\mathrm{g}} used for numerical verification of the ‘axisymmetric annulus’ approximation in § 6.2 , § 7.1 and Figure 7 . An orbit in a spherical potential (with mass M 𝑀 M and scale radius b 𝑏 b ) is uniquely specified by its peri/apocentre distances r p / a subscript 𝑟 p a r_{\mathrm{p/a}} , or equivalently by its generalised semi-major axis a g subscript 𝑎 g a_{\mathrm{g}} and eccentricity e g subscript 𝑒 g e_{\mathrm{g}} . We also provide the orbit’s azimuthal period T ϕ subscript 𝑇 italic-ϕ T_{\phi} around the cluster, its analytical A ∗ superscript 𝐴 A^{*} and Γ Γ \Gamma values, and the corresponding numerically computed values A num ∗ , Γ num subscript superscript 𝐴 num subscript Γ num A^{*}_{\mathrm{num}},\Gamma_{\mathrm{num}} , all to 3 significant figures.

Orbit of ${𝐑}_g$	Potential	$(r_p/b,r_a/b)$	$a_g/b$	$e_g$	$T_{\varphi }\sqrt{\frac{GM}{b^3}}$	$A^{\ast }$	$A_{\mathrm{num}}^{\ast }$	$\Gamma$	${\Gamma }_{\mathrm{num}}$
(I)	${\Phi }_{\mathrm{iso}}$	$(5.29,11.2)$	$8.2$	$0.36$	$171$	$0.00106$	$0.00106$	$0.685$	$0.686$
(II)	${\Phi }_{\mathrm{iso}}$	$(0.08,1.21)$	$0.65$	$0.88$	$18.0$	$0.253$	$0.255$	$0.0676$	$0.0650$
(III)	${\Phi }_{\mathrm{Hern}}$	$(0.08,1.21)$	$0.65$	$0.88$	$11.0$	$1.15$	$1.15$	$0.192$	$0.195$

Table 2. Table 2: Properties of Orbits (IV)-(VII) integrated in the flattened power-law potential ( 56 ) with β = 1 / 2 𝛽 1 2 \beta=1/2 and q = 0.94 𝑞 0.94 q=0.94 (c.f. Figure 10 ). The initial conditions of Orbits (IV)-(VII) are identical except for the initial azimuthal velocity v ϕ subscript 𝑣 italic-ϕ v_{\phi} . We take initial ( R , v R , ϕ , Z , v Z ) = ( b ,  0.1 ​ G ​ M / b ,  0.1 ,  0.2 ​ b ,  0.5 ​ G ​ M / b ) 𝑅 subscript 𝑣 𝑅 italic-ϕ 𝑍 subscript 𝑣 𝑍 𝑏 0.1 𝐺 𝑀 𝑏 0.1 0.2 𝑏 0.5 𝐺 𝑀 𝑏 (R,v_{R},\phi,Z,v_{Z})=(b,\,0.1\sqrt{GM/b},\,0.1,\,0.2b,\,0.5\sqrt{GM/b}) , and initial v ϕ subscript 𝑣 italic-ϕ v_{\phi} is given below.

Orbit	Potential	$v_{\varphi }/\sqrt{GM/b}$	$A_{\mathrm{num}}^{\ast }$	${\Gamma }_{\mathrm{num}}$
(IV)	${\Phi }_{\mathrm{FPL}}$	$1.35$	$0.0332$	$0.243$
(V)	${\Phi }_{\mathrm{FPL}}$	$0.95$	$0.182$	$0.192$
(VI)	${\Phi }_{\mathrm{FPL}}$	$0.35$	$0.552$	$0.016$
(VII)	${\Phi }_{\mathrm{FPL}}$	$0.05$	$0.626$	$-0.085$

Table 3. Table 3: Properties of Orbits (VIII)-(XI) integrated in the Miyamoto-Nagai potential ( 57 ) with b h / b ℓ = 1 subscript 𝑏 ℎ subscript 𝑏 ℓ 1 b_{h}/b_{\ell}=1 (c.f. Figure 11 ). The initial conditions of Orbits (VII)-(XI) are identical except for the initial vertical coordinate Z 𝑍 Z . We take initial ( R , v R , ϕ , v ϕ , v Z ) = ( b ℓ ,  0 ,  0 ,  0.25 ​ G ​ M / b ℓ ,  0 ) 𝑅 subscript 𝑣 𝑅 italic-ϕ subscript 𝑣 italic-ϕ subscript 𝑣 𝑍 subscript 𝑏 ℓ  0  0 0.25 𝐺 𝑀 subscript 𝑏 ℓ  0 (R,v_{R},\phi,v_{\phi},v_{Z})=(b_{\ell},\,0,\,0,\,0.25\sqrt{GM/b_{\ell}},\,0) , and initial Z 𝑍 Z is given below.

Orbit	Potential	$Z/b_{\ell }$	$A_{\mathrm{num}}^{\ast }$	${\Gamma }_{\mathrm{num}}$
(VIII)	${\Phi }_{\mathrm{MN}},b_h/b_{\ell }=1$	$0.1$	$0.256$	$0.153$
(IX)	${\Phi }_{\mathrm{MN}},b_h/b_{\ell }=1$	$1.0$	$0.122$	$0.042$
(X)	${\Phi }_{\mathrm{MN}},b_h/b_{\ell }=1$	$2.0$	$0.0392$	$-0.163$
(XI)	${\Phi }_{\mathrm{MN}},b_h/b_{\ell }=1$	$3.0$	$0.0140$	$-0.384$

Table 4. Table 4: Summary of U 𝑈 U and Γ ​ ( U ) Γ 𝑈 \Gamma(U) ranges that are possible for orbits 𝐑 g subscript 𝐑 g \mathbf{R}_{\mathrm{g}} in different classes of potential Φ Φ \Phi .

Type of the potential/orbit ${𝐑}_g$	$U$ range	$\Gamma$ range
General axisymmetric potential	$-\infty \le U\le \infty$	$-\infty \le \Gamma \le \infty$
Spherical potential (assuming $d\rho /dr\le 0$ and finite mass)	$1/3\le U\le \infty$	$0\le \Gamma \le 1$
Midplane of a thin disk	$U=1$	$\Gamma =1/3$
Vertical cylindrical potential	$U=0$	$\Gamma =-1/3$
Axisymmetric harmonic potential	$0\le U\le 1$	$-1/3\le \Gamma \le 1/3$
Spherical harmonic potential	$U=1/3$	$\Gamma =0$
Keplerian potential	$U\to \infty$	$\Gamma \to 1$
Spherical cusp potential (density $\rho \propto r^{-\beta }$)	$U=1/\left(3-\beta \right)$	$\Gamma =\beta /\left[3\left(4-\beta \right)\right]$

Equations

$$\frac{{\mathrm{d}}^2({\mathbf{r}}_i+{\mathbf{R}}_{\mathrm{b}})}{\mathrm{d}{t}^2}=-(\nabla \Phi {)}_{{\mathbf{R}}_{\mathrm{b}}+{\mathbf{r}}_i}-\frac{G{m}_j}{|{\mathbf{r}}_i-{\mathbf{r}}_j{|}^3}({\mathbf{r}}_i-{\mathbf{r}}_j),$$

$$\frac{{\mathrm{d}}^2\mathbf{r}}{\mathrm{d}{t}^2}=-\left[(\nabla \Phi {)}_{{\mathbf{R}}_{\mathrm{b}}+{\mathbf{r}}_1}-(\nabla \Phi {)}_{{\mathbf{R}}_{\mathrm{b}}+{\mathbf{r}}_2}\right]-\frac{\mu \mathbf{r}}{r^3},$$

$$[(\nabla \Phi {)}_{{\mathbf{R}}_{\mathrm{b}}+{\mathbf{r}}_i}{]}_{\alpha }$$

$$+\sum _{\beta }{\left(\frac{{\partial }^2\Phi }{\partial R_{\alpha }\partial R_{\beta }}\right)}_{{\mathbf{R}}_{\mathrm{b}}}{r}_{i,\beta }+O(r_i^2),$$

$$\frac{{\mathrm{d}}^2\mathbf{r}}{\mathrm{d}{t}^2}=-\sum _{\beta }{r}_{\beta }\frac{\partial }{\partial R_{\beta }}{\left(\nabla \Phi \right)}_{{\mathbf{R}}_{\mathrm{b}}}-\frac{\mu \mathbf{r}}{r^3},$$

$$\frac{{\mathrm{d}}^2\mathbf{r}}{\mathrm{d}{t}^2}=-(\mathbf{r}\cdot \nabla ){\left(\nabla \Phi \right)}_{{\mathbf{R}}_{\mathrm{b}}}-\frac{\mu \mathbf{r}}{r^3}.$$

$$H=H_0+H_1,$$

$$H_0=$$

$$H_1=$$

$$+{\left(\frac{{\partial }^2\Phi }{\partial x\partial y}\right)}_{{\mathbf{R}}_{\mathrm{b}}}xy+{\left(\frac{{\partial }^2\Phi }{\partial x\partial z}\right)}_{{\mathbf{R}}_{\mathrm{b}}}xz+{\left(\frac{{\partial }^2\Phi }{\partial y\partial z}\right)}_{{\mathbf{R}}_{\mathrm{b}}}yz,$$

$$H_1=$$

$$\frac{{\mathrm{d}}^2{\mathbf{R}}_{\mathrm{b}}}{\mathrm{d}{t}^2}$$

$$=-(\nabla \Phi {)}_{{\mathbf{R}}_{\mathrm{b}}}\left[1+O\left(r^2/|{\mathbf{R}}_{\mathrm{b}}{|}^2\right)\right].$$

$${\mathrm{d}}^2{\mathbf{R}}_{\mathrm{g}}/\mathrm{d}{t}^2=-(\nabla \Phi {)}_{{\mathbf{R}}_{\mathrm{g}}}.$$

$$x=$$

$$-\cos i\sin \Omega \left[(\cos E-e)\sin \omega +\sqrt{1-e^2}\sin E\cos \omega \right]),$$

$$y=$$

$$+\cos i\cos \Omega \left[(\cos E-e)\sin \omega +\sqrt{1-e^2}\sin E\cos \omega \right]),$$

$$z=$$

$$L=\sqrt{\mu a};J=L\sqrt{1-e^2};J_z$$

$$H_0=-\frac{\mu }{2a}=-\frac{{\mu }^2}{2L^2},$$

$$⟨H{⟩}_M=H_0+⟨H_1{⟩}_M,$$

$$⟨\mathcal{F}{⟩}_M$$

$$⟨H_1{⟩}_M=\frac{1}{2}\sum _{\alpha \beta }{\Phi }_{\alpha \beta }⟨r_{\alpha }{r}_{\beta }{⟩}_M.$$

$$⟨H_1{⟩}_M=\frac{1}{2}\left[{\kappa }_x^2⟨x^2{⟩}_M+{\kappa }_y^2⟨y^2{⟩}_M+{\kappa }_z^2⟨z^2{⟩}_M\right].$$

$$⟨H_1{⟩}_M$$

$$\times \left[(2+3e^2)\left(1+\frac{{\kappa }_{-}^2}{{\kappa }_{+}^2}{\cos }^{2}i\right)+5\frac{{\kappa }_{-}^2}{{\kappa }_{+}^2}{e}^2{\sin }^{2}i\cos 2\omega \right],$$

$$⟨H_1{⟩}_M=$$

$$\dot{\omega }=\frac{\partial ⟨H_1{⟩}_M}{\partial J}=-\frac{3}{2}\frac{{\kappa }^2}{n_{\mathrm{K}}}\sqrt{1-e^2},$$

$${\overline{⟨H⟩}}_M=-\frac{{\mu }^2}{2L^2}+{\overline{⟨H_1⟩}}_M,$$

$${\overline{⟨H_1⟩}}_M=\frac{1}{2}\sum _{\alpha \beta }{\overline{\Phi }}_{\alpha \beta }⟨r_{\alpha }{r}_{\beta }{⟩}_M.$$

$${\Phi }_{xx}$$

$$+\frac{1}{2}\cos 2{\varphi }_{\mathrm{g}}\left[{\left(\frac{{\partial }^2\Phi }{\partial R^2}\right)}_{{\mathbf{R}}_{\mathrm{g}}}-{\left(\frac{1}{R}\frac{\partial \Phi }{\partial R}\right)}_{{\mathbf{R}}_{\mathrm{g}}}\right],$$

$${\Phi }_{yy}$$

$$-\frac{1}{2}\cos 2{\varphi }_{\mathrm{g}}\left[{\left(\frac{{\partial }^2\Phi }{\partial R^2}\right)}_{{\mathbf{R}}_{\mathrm{g}}}-{\left(\frac{1}{R}\frac{\partial \Phi }{\partial R}\right)}_{{\mathbf{R}}_{\mathrm{g}}}\right],$$

Full text

††pubyear: 2019††pagerange: Secular dynamics of binaries in stellar clusters I: general formulation and dependence on cluster potential–F

1 Introduction

Orbital motion of two bound point masses perturbed by external forces represents one of the oldest problems in celestial mechanics (Murray & Dermott, 1999). It naturally emerges in a variety of astrophysically relevant situations: classical secular Laplace-Largange evolution of planetary systems (Murray & Dermott, 1999), Lidov-Kozai oscillations forced by a distant companion in a triple system (Lidov, 1962; Kozai, 1962), evolution of Oort Cloud comets driven by the Galactic tide (Heisler & Tremaine, 1986), and so on. Due to the weakness of the external perturbation, the evolution is usually secular in nature and operates on long time scales.

Dense stellar systems — globular and open clusters, nuclear star clusters, galaxies themselves (hereafter collectively called ‘clusters’) — represent ideal environments for perturbing binaries. Historically, these perturbations were considered predominantly in the context of encounters with individual passing stars (Heggie, 1975; Hut & Bahcall, 1983; Hut, 1983), which should occur rather frequently in dense clusters (although not always in the perturbative limit). In particular, secular changes of binary orbital elements caused by a passage of a single star were calculated by Rasio & Heggie (1995), Heggie & Rasio (1996), and Hamers (2018).

At the same time, the gravitational effect of the smooth mass distribution of the full cluster on binary orbital evolution has been explored mainly in the context of Oort Cloud formation and evolution driven by the Galatic tide (Heisler & Tremaine 1986, hereafter HT86). In this case the ‘binary’ is the Sun-comet system, and Galactic tides can be at least as important as stellar encounters for the orbital evolution of Oort Cloud comets. Effects of the spatially smooth Galactic tide have also been studied in the context of wide stellar binaries (Jiang & Tremaine, 2010) and long-period planetary systems in the Galactic bulge (Veras & Evans, 2013a).

Interest in binaries in stellar clusters has been renewed recently by the discoveries of mergers of compact-object binaries by the LIGO collaboration (e.g. Abbott et al. 2016, Abbott et al. 2017). A number of channels for the production and orbital evolution of merging binaries in dense stellar systems have been proposed following these discoveries (Wen, 2003; O’Leary et al., 2006; Miller & Lauburg, 2009; Bae et al., 2014; Rodriguez et al., 2015; Rodriguez et al., 2016a; Rodriguez et al., 2016b; Antonini et al., 2016; Askar et al., 2017; Chatterjee et al., 2017; Petrovich & Antonini, 2017; Hamers et al., 2018; Hoang et al., 2018; Rodriguez et al., 2018; Samsing, 2018; Fragione & Kocsis, 2018). Some of them involve Lidov-Kozai (Lidov 1962; Kozai 1962; hereafter LK) coupling with a third distant companion111Silsbee & Tremaine (2017) and Antonini et al. (2017) also considered Lidov-Kozai cycles in isolated triples to explain LIGO events., which is either captured from the cluster environment (e.g. Antonini et al. 2016) or is effectively represented by a central super-massive black hole (Antonini & Perets, 2012; Petrovich & Antonini, 2017; Hoang et al., 2018). Under certain conditions, LK oscillations can naturally drive the binary orbit to become highly eccentric, boosting gravitational wave emission and substantially speeding up binary coalescence. Similar ideas (with different sources of dissipation at periapsis) have been previously considered for explaining the origin of other exotic objects typically found in stellar clusters, such as blue stragglers (e.g. Knigge et al. 2009; Perets & Fabrycky 2009).

In this work we consider the general problem of secular evolution of binary orbital elements driven by the tidal field that arises due to the smooth mass distribution of an axisymmetric cluster in which the binary moves. We do this for an arbitrary axisymmetric cluster potential and set no constraints on the type of orbit that the binary can have in the cluster. This allows us to formulate the most general framework for treating tide-driven secular evolution of binaries, applicable to a variety of astrophysical systems.

In doing this we neglect the stochastic effect of individual stellar passages on the orbital elements of the binary. Separating the cumulative effect of multiple stellar encounters from the smooth Galactic tide is a non-trivial exercise, as demonstrated previously by Collins & Sari (2010) and Jiang & Tremaine (2010). Nevertheless, for the purposes of clarity, we prefer to focus here on the effects of the mean tidal field due to the smooth mass distribution inside the cluster — effects of encounters with individual stars will be incorporated later.

This paper (the first in a series) is devoted to the derivation of the general Hamiltonian governing secular orbital evolution of a binary in an arbitrary axisymmetric cluster potential, as well as to exploring the dependence of some characteristics of this Hamiltonian on the properties of the cluster potential and the binary’s orbit within it. It is structured as follows. In §2 we derive the tidal Hamiltonian for the dynamical evolution of binary orbital elements due to any axisymmetric perturbation when expanded to quadrupole order. In §§3-5 we average the tidal potential over both the binary’s inner orbit and then over many orbits of the binary around the cluster, arriving at a simple doubly-averaged (secular) Hamiltonian which describes long-term evolution of the binary’s orbital elements. The coefficients entering this secular Hamiltonian depend on the potential of the host system and the binary’s barycentric orbit within this potential, and we explore this dependence in detail in §6. We verify the time-averaging procedure numerically in §7. In §8 we discuss the limitations of our theory, and show how our general results are connected with various special cases already explored by others (also in Appendices B & C).

In a companion paper (Hamilton & Rafikov 2019a, hereafter ‘Paper II’) we provide a complete study of the dynamics that result from the secular tidal Hamiltonian derived in this work, calculate the timescale and amplitude of eccentricity oscillations, and verify our theory with direct N-body integration.

2 Hamiltonian with cluster tides

Let us consider a binary system with semi-major axis $a$ and eccentricity $e$, consisting of point masses $m_{1}$ and $m_{2}$. Binary components gravitationally interact with each other and a fixed smooth background potential $\Phi$ of a much more massive system, which we will later take to be axisymmetric. The application we have most readily in mind is that of binary stars in the mean field potential of a globular or nuclear star cluster, and for this reason we will frequently refer to $m_{1}$ and $m_{2}$ as ‘stars’ and to the background system as ‘the cluster’. However it should be borne in mind that our analysis works for any system of two gravitationally bound objects (binary black holes, comet-Sun system, etc.) moving in any axisymmetric potential (galaxy, open cluster, young stellar cluster, etc).

Throughout this work we will refer to the binary’s orbit around the cluster as the ‘outer orbit’, while the orbit of the binary components about their common barycentre will be called the ‘inner orbit’, to coincide with the standard terminology in LK studies (e.g. Naoz 2016). To describe the outer and inner orbits we set up two coordinate systems — see Figure 1 for illustration.

The first, given by $\mathbf{R}=(X,Y,Z)$, has its origin at the centre of the cluster. In this coordinate system, the radius vector of the outer orbit, i.e. from the cluster centre to the barycentre of the binary is given by $\mathbf{R}_{\mathrm{b}}=(X_{\mathrm{b}},Y_{\mathrm{b}},Z_{\mathrm{b}})$. The second (non-inertial) coordinate system has its origin at $\mathbf{R}_{\mathrm{b}}$, and its axes are fixed to be aligned with those of the first system, so only its origin moves. The position of star $i=1,2$ in the non-inertial system is then given by $\mathbf{r}_{i}=(x_{i},y_{i},z_{i})$. The position of star $i$ relative to the centre of the cluster is $\mathbf{R}_{i}=\mathbf{R}_{\mathrm{b}}+\mathbf{r}_{i}=(X_{\mathrm{b}}+x_{i},Y_{\mathrm{b}}+y_{i},Z_{\mathrm{b}}+z_{i})$ and the barycentre is at $\mathbf{R}_{\mathrm{b}}=(m_{1}\mathbf{R}_{1}+m_{2}\mathbf{R}_{2})/(m_{1}+m_{2})$.

The equation of motion of star $i=1,2$ is then

[TABLE]

for $i\neq j$, where the subscript on derivatives means that we evaluate the derivative at $\mathbf{R}_{\mathrm{b}}+\mathbf{r}_{i}$.

Defining the relative position $\mathbf{r}=(x,y,z)\equiv\mathbf{r}_{1}-\mathbf{r}_{2}$, and $\mu\equiv G(m_{1}+m_{2})$, one obtains from (1)

[TABLE]

which is the general equation of relative motion of the binary components.

2.1 Tidal approximation

We now employ the tidal approximation, which means that in equation (2) we expand the potential $\Phi(\mathbf{R}_{\mathrm{b}}+\mathbf{r}_{i})$ around $\mathbf{R}_{\mathrm{b}}$. The Cartesian components of the vector $\nabla\Phi$ at position $\mathbf{R}_{\mathrm{b}}+\mathbf{r}_{i}$ are

[TABLE]

where $\alpha$, $\beta$ are the coordinate indices, so that e.g. $R_{\alpha}$ and $r_{i,\alpha}$ represent the components of $\mathbf{R}=(X,Y,Z)$ and $\mathbf{r}_{i}=(x_{i},y_{i},z_{i})$ respectively. We expect that terms $O(r_{i}^{2})$ will be subdominant because the distance to the centre of the cluster ($\sim|\mathbf{R}_{\mathrm{b}}|$) is much greater than the binary separation $r$. In this approximation we find

[TABLE]

with $r_{\alpha}$ the components of $\mathbf{r}=(x,y,z)$. Keeping in mind that, since the axes of the two coordinate systems are aligned, we may interchange $\partial/\partial R_{\alpha}$ with $\partial/\partial r_{\alpha}$, we can also write this as

[TABLE]

With $x,y,z$ as our canonical coordinates and $p_{x}=\dot{x},p_{y}=\dot{y},p_{z}=\dot{z}$ as the corresponding momenta, these equations of motion may be derived from the time-dependent Hamiltonian

[TABLE]

where

[TABLE]

and $H_{1}\ll H_{0}$. To compress the notation we write $(\partial^{2}\Phi/\partial r_{\alpha}\partial r_{\beta})_{\mathbf{R}_{\mathrm{b}}}\equiv\Phi_{\alpha\beta}$ so that the perturbing (‘tidal’) Hamiltonian reads

[TABLE]

where we sum over $\alpha,\beta=x,y,z$.

The dominant part of the Hamiltonian $H_{0}$ corresponds to the motion of an isolated binary star about its own barycentre, and has no explicit time dependence. The perturbing term $H_{1}$ takes into account the tidal effects of the external potential, which will drive the secular evolution of the binary orbital elements. It implicitly depends on time through $\mathbf{R}_{\mathrm{b}}(t)$, which we look at next.

According to equation (1) the evolution of $\mathbf{R}_{\mathrm{b}}$ is governed by

[TABLE]

The small correction terms on the right hand side of this equation mean that, in general, the motion of $\mathbf{R}_{\mathrm{b}}$ in the cluster does not coincide exactly with that of a test particle. However, at the level of accuracy needed in this work we can neglect this difference and assume that $\mathbf{R}_{\mathrm{b}}$ coincides with the ‘guide’ radius vector $\mathbf{R}_{\mathrm{g}}$, which evolves according to the equation of motion of a test particle in the cluster potential,

[TABLE]

In other words, in the following we set $\mathbf{R}_{\mathrm{b}}=\mathbf{R}_{\mathrm{g}}$ and calculate $\mathbf{R}_{\mathrm{g}}(t)$ using equation (11).

Our neglect of the terms quadratic and higher order in $r_{i}$ in equation (3) is equivalent to the so-called ‘quadrupole approximation’ in the hierarchical three-body problem. Keeping the next (quadratic) term in the expansion would correspond to the ‘octupole approximation’, and so on. In Appendix E we describe the extension of our tidal Hamiltonian to octupole order and provide a connection to the LK problem in the octupole approximation.

2.2 Orbital elements and Delaunay variables

We now introduce standard orbital elements in the frame of the binary. The reference direction is taken to be the $X$ direction and the reference plane the $(X,Y)$ plane (see Figure 1; we will later take the $Z$ axis to be the symmetry axis of the potential but the assumption of axisymmetry is not needed at the moment). We define binary argument of pericentre $\omega$, inclination $i$, longitude of ascending node $\Omega$ and mean anomaly $M$ relative to this reference plane and direction. When written in orbital elements the relative coordinates $\mathbf{r}=(x,y,z)$ become (Murray & Dermott, 1999):

[TABLE]

[TABLE]

[TABLE]

and $M=E-e\sin E$ where $E$ is the eccentric anomaly. It is important that the orbital elements are defined with respect to a reference frame with axis directions fixed in time (cf. Brasser 2001; Veras & Evans 2013a; Correa-Otto et al. 2017). In the limit of the cluster tide going to zero these orbital elements stay constant.

For dynamical studies it is often more convenient to use Delaunay variables, in which the actions

[TABLE]

are complemented by their conjugate angles $M$, $\omega$, $\Omega$. We will use them extensively in Paper II.

Since Delaunay variables are angle-action variables, we can easily identify the conserved quantities in the Hamiltonian. The dominant part of the Hamiltonian (7) reads

[TABLE]

while the perturbing Hamiltonian is given by equation (9) with $x,y$ and $z$ given by equations (12), (13) and (14) respectively (or their Delaunay equivalents).

3 Averaging the tidal Hamiltonian

Dynamics of binaries in stellar clusters benefits from a natural separation of scales. For example, a Solar mass binary with $a=20\,\mathrm{AU}$ has an inner orbital period of $\sim 100$ years, while its outer orbit around a globular cluster might have a period of $\sim 10^{5}$ years. As we show in Paper II, the resulting secular evolution of the binary’s orbital elements due to the tidal potential of the cluster may take $\sim 10^{8}$ years.

This naturally allows us to simplify our Hamiltonian (9), first by integrating out the fast evolution of the mean anomaly $M$ of the inner orbit (‘single-averaging’, see §3.1), and then by also integrating over many (outer) orbits of the binary around the cluster (‘double averaging’, see §3.2).

3.1 Singly-averaged Hamiltonian: averaging over the mean anomaly $M$

We begin by averaging over the shortest timescale in the problem, namely over the inner orbital motion of the binary components around their common barycentre. Our singly-averaged Hamiltonian is

[TABLE]

where the average of a quantity $\mathcal{F}$ over the mean anomaly $M$ is defined as

[TABLE]

The coefficients $\Phi_{\alpha\beta}$ depend on time only through $\mathbf{R}_{\mathrm{g}}(t)$, which is a ‘slow’ variable, so

[TABLE]

For reference, the full algebraic expressions for $\langle r_{\alpha}r_{\beta}\rangle_{M}$ in terms of orbital elements are given in Appendix A. The singly-averaged Hamiltonian (19) incorporating these expressions is completely general and can be used to describe orbital evolution of binaries moving in an arbitrary external potential.

Obviously we have eliminated the angle $M$, therefore the conjugate action $L=\sqrt{\mu a}$ is conserved, and so the binary’s semi-major axis $a$ is constant. When written in Delaunay variables, the singly-averaged Hamiltonian $\langle H_{1}\rangle_{M}$ is a function of $J,J_{z},\omega,\Omega$ and the time $t$ through the time-dependent coefficients $\Phi_{\alpha\beta}(\mathbf{R}_{\mathrm{g}}(t))$.

3.1.1 Example: orbits in a harmonic potential

For illustration, as well as to connect to subsequent results, we consider a binary orbiting in a globular cluster with a triaxial constant-density core. For orbits in this core the potential is that of a three-dimensional harmonic oscillator with frequencies $\kappa_{\alpha}$, namely $\Phi=\sum_{\alpha}\frac{1}{2}\kappa_{\alpha}^{2}R_{\alpha}^{2}$. Then $\Phi_{\alpha\beta}=\kappa_{\alpha}^{2}\delta_{\alpha\beta}$ so the singly-averaged Hamiltonian (19) becomes

[TABLE]

Let us now consider an axisymmetric core where the symmetry axis is the $Z$ axis. Then $\kappa_{x}=\kappa_{y}$ and the binary’s outer orbit fills a section of a cylindrical surface with elliptical cross-section (aligned with the $Z$ axis). Using equations (59)-(61) we end up with

[TABLE]

where $\kappa^{2}_{\pm}\equiv\kappa^{2}_{x}\pm\kappa^{2}_{z}$.

Note that the dependence on the longitude of the ascending node $\Omega$ has dropped out of this Hamiltonian, so the $z$-component of binary angular momentum $J_{z}$ is conserved. That is, in the reference frame of a binary orbiting an axisymmetric harmonic potential, the perturbation due to the tidal field of the cluster effectively becomes axisymmetric after averaging only over the inner orbit of the binary (single-averaging), not its outer orbit around the cluster (double-averaging). This is despite the fact that the outer orbit itself is not axisymmetric in general, even after averaging over a long time interval (its projection onto the $(X,Y)$ plane is an ellipse centred at the origin). This property does not hold for arbitrary potentials.

Things simplify further if we assume the core to be spherically symmetric. Without loss of generality we can then assume the outer orbit of the binary to be in the $Z=0$ plane, and put all frequencies equal to $\kappa$ (so that $\kappa_{-}=0$). We find that

[TABLE]

This singly-averaged Hamiltonian is now also independent of the argument of pericentre $\omega$. As a result, in this case there is no evolution of eccentricity or inclination of the binary. The only variation of its orbital elements is apsidal precession at the rate

[TABLE]

independent of the orientation of the binary orbit (i.e. $\omega$, $\Omega$, $i$). Here $n_{\rm K}=\sqrt{\mu/a^{3}}$ is the Keplerian mean motion of the binary.

3.2 Doubly-averaged Hamiltonian: averaging over time

As we already mentioned, binary orbital elements change significantly on timescales that are much longer than the outer orbital period of the binary around its host system. For that reason, it makes sense to average $\langle H\rangle_{M}$ over many outer orbits. Indicating such time-averages with an over bar, we write:

[TABLE]

where the doubly-averaged perturbing Hamiltonian $\overline{\langle H_{1}\rangle}_{M}$ differs from its singly-averaged predecessor $\langle H_{1}\rangle_{M}$ (equation (19)) only in that each $\Phi_{\alpha\beta}$ is now replaced by its time-averaged value $\overline{\Phi}_{\alpha\beta}$:

[TABLE]

This works because the outer orbit $\mathbf{R}_{\mathrm{g}}(t)$ only enters $\Phi_{\alpha\beta}$ and not $r_{\alpha}r_{\beta}$.

Equation (25) is the doubly-averaged perturbing Hamiltonian and is the main result of this section. It describes the secular evolution of the orbital elements of any binary perturbed by a smooth external potential $\Phi$. However in its current abstract form it is not of much use. In the following section we show how the time-averages $\overline{\Phi}_{\alpha\beta}$ may be calculated for cluster potentials possessing certain symmetries, culminating in the expressions (42), (43).

4 Time-averaging in axisymmetric potentials

So far we did not need to specify anything about the potential $\Phi$. However, we will now focus on binaries orbiting in fixed axisymmetric potentials (§§4.1-4.2). We describe the time-averaging procedure in spherical clusters (§4.3) and then extend it to general axisymmetric potentials (§4.4).

4.1 $\Phi_{\alpha\beta}$ in cylindrical coordinates

In an axisymmetric cluster we can choose the symmetry axis of the potential to be the $Z$ axis (like in Figure 1). Then it makes sense to write down the derivatives $\Phi_{\alpha\beta}$ in the cylindrical $(R,\phi,Z)$ coordinate system with origin at the centre of the cluster, where $R=\sqrt{X^{2}+Y^{2}}$ and $\phi=\tan^{-1}(Y/X)$. The axisymmetric potential may then be expressed as $\Phi(R,Z)$, and we find

[TABLE]

Here $\phi_{\mathrm{g}}$ is the azimuthal coordinate of $\mathbf{R}_{\mathrm{g}}$, namely $\tan^{-1}({Y_{\mathrm{g}}}/{X_{\mathrm{g}}})$, and again the subscripts on derivatives mean ‘evaluated at position $\mathbf{R}_{\mathrm{g}}(t)$’. The coefficients $\Phi_{\alpha\beta}$ have certain symmetry properties which will become important when we consider their time-averaged values in §4.2.

4.2 Orbit families and non-commensurable frequencies

We now consider which orbit families are possible in general axisymmetric potentials. Numerical orbit integration confirms that most orbits in axisymmetric potentials are regular and appear to respect a third integral of motion $I_{3}$ alongside energy $\mathcal{E}$ and the $Z$-component of angular momentum $\mathcal{L}_{Z}$ (Binney & Tremaine, 2008; Merritt, 2013). The vast majority of these regular orbits form a tube, or ‘torus’, around the symmetry axis: in an oblate potential they are short-axis tube orbits, while in a prolate potential they are (inner- or outer-) long axis tube orbits. Other possibilities include near-resonant regular orbits and irregular (chaotic) orbits, but both of these are typically scarce compared to the tubes.

We will ignore chaotic orbits in what follows since they are very rare in axisymmetric potentials (Regev, 2006). We are left with tube orbits and (near-)resonant non-tube orbits. The resonant family corresponds to $\mathbf{R}_{\mathrm{g}}(t)$ having commensurable frequencies. Mathematically, if we denote the frequencies of motion of $\mathbf{R}_{\mathrm{g}}(t)$ in each direction by the vector $\mathbf{\Omega}=(\Omega_{R},\Omega_{\phi},\Omega_{Z})$, we must consider whether there exists any triple of integers $\mathbf{n}=(n_{1},n_{2},n_{3})$ such that

[TABLE]

The role of commensurabilities and near-commensurabilities will be discussed in §§7-8.

If the frequencies are non-commensurable (i.e. relation (32) does not hold), then the outer orbit of the binary will trace out a non-repeating path around the cluster. Over time this path will densely fill a 3D axisymmetric torus whose symmetry axis is $Z$. We may therefore replace the time-average of a function following an orbit $\mathbf{R}_{\mathrm{g}}(t)$ with a weighted (by the time the orbit spends at a given point) volume-average over the torus.

Since the torus is axisymmetric, time-averaging over non-commensurable orbits inevitably involves integrating the expressions (4.1)-(31) over azimuthal angle $\phi_{\mathrm{g}}$. As a result, time-averages of $\Phi_{\alpha\beta}$ become

[TABLE]

see equations (4.1)-(31). In practice, vanishing of $\overline{\Phi}_{xy}$, $\overline{\Phi}_{xz}$, $\overline{\Phi}_{yz}$ typically requires averaging over many outer orbits — see §7.

4.3 Time-averages in spherical potentials

In spherical potentials the outer orbit of the binary remains in the same plane, which can be chosen to coincide with the $(X,Y)$ plane. In this plane the coefficients $\Phi_{xz}$ and $\Phi_{yz}$ vanish identically. In other words, equation (36) holds true even without averaging over the outer orbit. At the same time, $\overline{\Phi}_{xy}$ asymptotically tends to zero only upon averaging over many orbits, as we will see later in §7.1.

In the $(X,Y)$ plane the path of $\mathbf{R}_{\mathrm{g}}$ is a rosette, assuming it has non-commensurable radial and azimuthal frequencies; see Figure 7 for illustration. Over time the rosette densely fills an axisymmetric annulus with inner and outer radii corresponding to the pericentre $r_{\mathrm{p}}$ and apocentre $r_{\mathrm{a}}$ of the outer orbit $\mathbf{R}_{\mathrm{g}}(t)$. When discussing spherical potentials we will sometimes refer to this as the ‘axisymmetric annulus approximation’.

In this case it is easy to write down an analytical formula for the averages $\overline{\Phi}_{\alpha\beta}$ in terms of $r_{\mathrm{p}}$ and $r_{\mathrm{a}}$, as averaging over $\mathrm{d}t$ can be replaced with averaging over $\mathrm{d}R$ via $\mathrm{d}t=v_{R}^{-1}\mathrm{d}R$, where $v_{R}=[2(\mathcal{E}-\Phi(R))-\mathcal{L}^{2}/R^{2}]^{1/2}$ is the radial velocity. Specific energy $\mathcal{E}$ and angular momentum $\mathcal{L}$ of the outer orbit $\mathbf{R}_{\mathrm{g}}(t)$ in a spherical potential $\Phi$ can be explicitly expressed as function of $r_{\mathrm{p}}$ and $r_{\mathrm{a}}$ as follows:

[TABLE]

With this in mind, we can write the time-average of any radially-dependent function $\mathcal{F}(R)$ as

[TABLE]

Remembering that only the azimuthally-averaged versions of $\Phi_{\alpha\beta}$ provide non-zero contributions (see §4.2 and equations (33)-(36)), we see that in spherical potentials the time-averages $\overline{\Phi}_{\alpha\beta}$ can be calculated in a straightforward fashion via integration over radius $R$.

4.4 Time-averages in axisymmetric potentials

We would like to generalise the approach of §4.3 to axisymmetric potentials $\Phi(R,Z)$. This would involve averaging $\Phi_{\alpha\beta}$ over the $(R,Z)$ cross-section of an axisymmetric torus filled by the outer orbit of the binary — see the central columns of Figures 10 & 11 for examples of such cross-sections. However, there are several difficulties with this approach.

First, each $\mathrm{d}R\,\mathrm{d}Z$ element of the cross-section enters the averaging procedure with a certain weight proportional to the time the orbit spends in it. Calculating this weight is not trivial and involves the knowledge of the meridional velocity ($v_{R}$, $v_{Z}$) structure. For a general axisymmetric potential this calculation cannot be done analytically.

Second, even the shape of the outer boundary of the cross-section cannot be determined analytically for a general axisymmetric potential. The difficulty is that the knowledge of $\mathcal{E}$ and $Z$-component of angular momentum $\mathcal{L}_{Z}$ (which are integrals of motion in a general axisymmetric potential) is not sufficient to determine the shape of the meridional cross section of the torus: one also needs to know a third integral of motion $I_{3}$. The exact shape of the torus cross-section is known only for orbits in Staeckel potentials (Binney & Tremaine, 2008), since only for those do we have analytic expression for the third integral $I_{3}$. Even then, writing down a formula for the time-averaged coefficients $\overline{\Phi}_{\alpha\beta}$ is tiresome because of the awkward confocal-ellipsoidal coordinate system involved (Binney & Tremaine, 2008) and the complicated functional form of the third integral.

For these reasons, in this work we usually222There are special cases in certain axisymmetric potentials where we can compute $A,\Gamma$ (semi-)analytically, see §6.1 & §6.3.1. compute $A$ and $\Gamma$ in axisymmetric potentials by directly integrating the orbit of a guide $\mathbf{R}_{\mathrm{g}}(t)$ numerically using equation (11) for a given set of initial conditions $(R,v_{R},\phi,v_{\phi},Z,v_{Z})$, where $v_{R}$ is the velocity in the direction of increasing $R$, etc. This orbit is then used to carry out the time-average of $\Phi_{\alpha\beta}$ using a method outlined in Appendix F.

Note that, unlike in the spherically symmetric case, $\overline{\Phi}_{xz}$ and $\overline{\Phi}_{yz}$ no longer vanish identically due to a symmetry of the potential. Nevertheless, equations (35)-(36) are still fulfilled upon averaging over many outer orbits.

5 The secular Hamiltonian

Despite the fact that in general axisymmetric potentials we cannot write down a useful analytic expression for time-averages, we can still continue our derivation of the secular Hamiltonian owing to the symmetries of the $\overline{\Phi}_{\alpha\beta}$ coefficients (equations (33)-(36)). These symmetry properties allow us to greatly simplify the doubly-averaged perturbing Hamiltonian (25) so that it reads:

[TABLE]

Let us define the quanitities

[TABLE]

(note that the constants $A$ and $B$ are not related to the Oort constants!) and write $x,y,z$ in terms of orbital elements using equations (59)-(61). Then the Hamiltonian (40) becomes

[TABLE]

and $H_{1}^{*}$ is the ‘dimensionless Hamiltonian’

[TABLE]

Note that $H_{1}^{*}$ involves only a single dimensionless parameter $\Gamma$, while $C$ depends on $A$ (which has units of (frequency)2). In Paper II we will see that $\Gamma$ determines the phase space structure of the Hamiltonian while $A$ sets the timescale for secular evolution. All the information about the cluster properties and the characteristics of the (outer) orbit of the binary enter the Hamiltonian through these two parameters only. In §6 we investigate in detail how these key quantities depend on the form of the background potential and the outer orbit of the binary within the potential.

The dependence of the Hamiltonian upon the longitude of ascending node $\Omega$ has dropped out under time-averaging and so the $z$-component of angular momentum $J_{z}$ is conserved, as we would expect for an axisymmetric perturbation (which the cluster potential looks like from the binary frame upon averaging over many outer orbits). The dimensionless quantity $J_{z}/L=\sqrt{1-e^{2}}\cos i$ is an integral of motion, which implies that variations of eccentricity must come at the expense of changes in inclination and vice versa, just as in the LK mechanism (Lidov, 1962; Kozai, 1962; Naoz et al., 2013).

Finally, we note that the doubly-averaged perturbing Hamiltonian (42) appears very similar to the singly-averaged one derived for the example of an axisymmetric harmonic potential in §3.1.1 (equation (21)). Indeed, comparing equations (21) and (42) one might be tempted to say that axisymmetric harmonic potentials have $\Gamma=-\kappa_{-}^{2}/3\kappa_{+}^{2}$. However, this correspondence is a mathematical coincidence: the assumption of non-commensurability (§4.2) does not apply to harmonic potentials, for which all orbits are closed non-precessing ellipses. Despite their similarities the Hamiltonians (21) and (42) are different objects derived under different approximations.

5.1 Orbits in a Kepler potential: link to the Lidov-Kozai mechanism

Another example of such a mathematical coincidence is represented by the well known test particle quadrupole Lidov-Kozai (LK) problem (Lidov, 1962; Kozai, 1962). The Hamiltonian for this problem takes the form (43) if we were to set $\Gamma=1$. However, we have derived this Hamiltonian under the assumption that $\mathbf{R}_{\mathrm{g}}$ fills an axisymmetric torus, while in the LK case $\mathbf{R}_{\mathrm{g}}$ moves in a Keplerian ellipse, which in the test particle limit does not precess. Nevertheless, it is known that for elliptical orbits the time-averaged tidal Keplerian potential is exactly axisymmetric to quadrupole order (e.g. Katz et al. 2011; Naoz et al. 2011), and so (42) does in fact hold.

In Appendix B we show explicitly how the LK Hamiltonian is recovered in the ‘test particle quadrupole’ approximation from the singly-averaged equation (19) in the limit that the background potential $\Phi$ in which the binary orbits arises from a point mass at the origin. We recover the LK Hamiltonian exactly if we set $\Gamma=1$ in (43).

5.2 Epicyclic orbits in a disk: link to Heisler &

Tremaine (1986)

For wide binaries in the solar neighbourhood, the tidal potential of the Galactic disk can provide the dominant torque, as shown by Heisler & Tremaine (1986) for the Oort Cloud comets. Averaged over many orbits of the Sun around the Galaxy, the Galactic disk provides an axisymmetric tide onto the binary. In Appendix C we show how to calculate $\overline{\Phi}_{\alpha\beta}$ in the case where $\mathbf{R}_{\mathrm{g}}$ performs epicyclic motion in a disk. It is then easy to recover the tidal Hamiltonian of HT86 from (42). We reproduce the dimensionless version of HT86’s Hamiltonian by setting $\Gamma=1/3$ in (43).

6 Dependence of Hamiltonian coefficients $A$ and $\Gamma$ on the cluster potential and binary orbit

All of the information about the effect of the tidal potential on secular dynamics of the binary is contained in the two crucial quantities $A$ and $\Gamma$, which are constructed from the time-averages $\overline{\Phi}_{zz}$ and $\overline{\Phi}_{xx}$, see equation (41).

The quantity $A$ is a direct measure of the strength of the potential. Its influence on the dynamics is trivial: it enters the problem only as a proportionality constant of the Hamiltonian (equation (42)), and therefore merely sets the (inverse of the) secular timescale. In addition, $A$ is also a fairly intuitive quantity: if a binary is in a strong tidal potential we expect it will have a large $A$.

The meaning of $\Gamma$ is less intuitive than $A$ although its influence upon the system is quite profound. In Paper II we will see that the phase portrait of the Hamiltonian $H_{1}^{*}$ undergoes several bifurcations as we change the value of $\Gamma$, altering the dynamics completely. In particular, we show that there are four qualitatively different regimes — (i) $\Gamma>1/5$, (ii) $0<\Gamma\leq 1/5$, (iii) $-1/5<\Gamma\leq 0$, and (iv) $\Gamma\leq-1/5$. The value of $\Gamma$ is so important because, for instance, high-eccentricity excitation is ubiquitous for binaries in regime (i), whereas it is much less readily achieved by binaries in regime (ii). Hence, most effort in this section will be directed towards understanding which cluster potentials and outer binary orbits give rise to which values of $\Gamma$. So far we have seen that $\Gamma=1$ for any orbit in a Keplerian potential, and that $\Gamma=1/3$ for epicyclic orbits in a thin disk. In this section we explore in more detail how the values of $\Gamma$ (and $A$) depend on the form of the background potential $\Phi$ and the binary’s outer orbit $\mathbf{R}_{\mathrm{g}}$ within it.

We start by stating some general properties of $A$ and $\Gamma$ in §6.1. We then discuss the behavior of these parameters in certain spherical (§6.2) as well as axisymmetric (§6.3) potentials.

6.1 General properties of $A$ and $\Gamma$

Writing down Poisson’s equation in cylindrical coordinates

[TABLE]

and using equations (33)-(34) & (41), one can easily show that in a general axisymmetric potential

[TABLE]

where $\overline{\rho}$ is the cluster density in the vicinity of the binary time-averaged over many outer orbital periods. Then, defining the dimensionless ratio

[TABLE]

we can write $\Gamma$ quite generally as

[TABLE]

The function $\Gamma(U)$ is plotted in Figure 2.

In principle there are no limits on the values $U$ can take, although in practice, achieving values of $U$ less than $-1$ (and hence $\Gamma>1$) may require rather contrived orbits. An example of such an orbit with $U<-1$ and $\Gamma>1$ is given in Appendix D (see Figure 13). Note that $U<-1$ necessarily implies that $A<0$, see equation (46).

Somewhat stronger statements can be formulated for realistic spherically symmetric cluster potentials, as we show in Appendix D. In particular, one can demonstrate that in such potentials $A>0$, $B\geq 0$, and $0\leq\Gamma\leq 1$. In non-spherical potentials negative values of $\Gamma$ become possible for certain binary orbits as we will show in §6.3.

It is instructive to consider the values of $\Gamma$ for some specific potentials $\Phi$.

• •

In the case of a Keplerian cluster potential, i.e. a spherical point mass potential with vanishingly small density $\rho$ outside the centre, one has $\Phi_{zz}>0$, $\overline{\rho}\to 0$, $U\to+\infty$ and $\Gamma\to 1$ (see §5.1).

• •

In a spherical harmonic potential, symmetry dictates that $\Phi_{zz}=(1/3)\nabla^{2}\Phi=(4\pi/3)G\rho$ so that $U=1/3$ and $\Gamma=0$ (see §3.1.1).

• •

In a spherical cluster with a cusped density distribution $\rho\propto r^{-\beta}$ with $0<\beta<3$ (having finite mass at the centre) we have $\Gamma=\beta/[3(4-\beta)]$, see Appendix D.1.

• •

In a thin galactic disk, assuming that $\Phi_{zz}$ dominates over other spatial derivatives in Poisson’s equation, one has $\Phi_{zz}\approx 4\pi G\rho$; hence we find $U=1$ and $\Gamma=1/3$ (see §5.2).

• •

In the opposite limit of a ‘cylindrical’ (or highly prolate) potential $\Phi=\Phi(R)$ with no $Z$-dependence, one has $\Phi_{zz}=0$, $U=0$ and $\Gamma=-1/3$.

The values (or ranges) of $U$ and $\Gamma(U)$ for these and some other types of cluster potential are summarized in Table 4. We stress again that even though applying the ‘axisymmetric annulus’ approximation gives the correct results for Keplerian and spherical harmonic potentials, this is a mathematical coincidence unless the outer orbit of the binary in these potentials is circular (see §3.1.1 & 5.1).

6.2 Behavior of Hamiltonian characteristics in some spherical potentials

In spherical potentials the values of $\overline{\Phi}_{\alpha\beta}$ that enter $A$ and $\Gamma$ are computed using the analytic expression (39), which for a fixed potential depends only on the peri/apocentre $(r_{\mathrm{p}},r_{\mathrm{a}})$ of the binary’s outer orbit $\mathbf{R}_{\mathrm{g}}(t)$. We can define the outer orbit’s generalised semi-major axis $a_{\mathrm{g}}$ and generalised eccentricity $e_{\mathrm{g}}$ in terms of the peri/apocentre as

[TABLE]

These reduce to the usual orbital elements in the case of a Keplerian potential. These variables fully characterize the outer orbit of the binary in a given spherical potential.

In any spherical potential with scale radius $b$ and total mass $M$ we can also construct the dimensionless parameter $A^{*}\equiv A/(GM/b^{3})$; this normalization arises because $A$ is constructed from the second derivatives of the potential, which are of order333Note that $2\pi/\sqrt{GM/b^{3}}$ is roughly $\sim T_{\phi}$, the characteristic azimuthal orbital period of the binary around the cluster, so that $A\sim 4\pi^{2}A^{*}/T_{\phi}^{2}$. $GM/b^{3}$, see equation (41). This allows us to estimate

[TABLE]

Both $A^{*}$ and $\Gamma$ are dimensionless numbers which, for a given potential, depend only on $a_{\mathrm{g}}$ and $e_{\mathrm{g}}$. In the following we will explore the dependence of $A^{*}$ (rather than $A$, which also depends on the cluster mass and size) and $\Gamma$ on the shape of the potential and the binary orbit in it.

We use the following spherically symmetric potentials:

(i) the isochrone potential (which has a constant-density core and half mass radius $r_{\mathrm{h}}=3.06b$)

[TABLE]

(ii) the Plummer potential (also has a core and $r_{\mathrm{h}}=1.31b$)

[TABLE]

(iii) the Hernquist potential (has no core, and $r_{\mathrm{h}}=2.41b$)

[TABLE]

(iv) the NFW potential (has no core and has a divergent mass profile)

[TABLE]

The NFW potential arises from a density distribution

[TABLE]

In equation (53) the quantity $M$ is not the mass of the model (which is formally infinite), just a constant with units of mass.

For illustration, in Figure 7 (left panels) we show three examples of numerically integrated orbits in some of these potentials. (Each orbit was integrated for $100$ azimuthal periods; the first few periods are highlighted in red). The first two (‘I’ and ‘II’) orbit the isochrone potential (50), which has a constant density core for $r\lesssim b$. The third (‘III’) orbits the Hernquist potential (52), which is coreless. In Table 1 we list the peri/apocentre $r_{\mathrm{p/a}}$, semi-major axis $a_{\mathrm{g}}$, eccentricity $e_{\mathrm{g}}$, azimuthal period $T_{\phi}$, and the values of $A^{*}$ and $\Gamma$ calculated using equation (39). We also provide values of $A^{*}_{\mathrm{num}},\Gamma_{\mathrm{num}}$ obtained by direct averaging of $\Phi_{\alpha\beta}$ along each numerically integrated outer orbit (see Appendix F), to which we will return when discussing the validity of the axisymmetric averaging approximation in §7.

6.2.1 Behavior of $A^{*}$

In Figure 3 we plot $\log_{10}A^{*}$ in the $(a_{\mathrm{g}},e_{\mathrm{g}})$ plane for the potentials (50)-(53). We see that $A^{*}$ is a strong function of $a_{\mathrm{g}}$ but a weaker function of $e_{\mathrm{g}}$. The dependence on $e_{\mathrm{g}}$ emerges predominantly for orbits with $e_{\mathrm{g}}\gtrsim 0.5$; it is rather weak at all $e_{\mathrm{g}}$ for orbits with $a_{\mathrm{g}}\lesssim b$, where $b$ is the scale radius of the potential in question. The difference in radial $A^{*}$ behavior between different potentials is most pronounced for orbits with $a_{\mathrm{g}}\lesssim b$. In this region there is a sharp increase in $A^{*}$ in the uncored (Hernquist and NFW) potentials, but a much shallower gradient in the cored potentials (isochrone and Plummer).

We can make the comparison more quantitative by examining the radial profile of $A^{*}$ for circular outer orbits (of radius $a_{\mathrm{g}}$ and eccentricity $e_{\mathrm{g}}=0$). Then $A^{*}$ is a function of $a_{\mathrm{g}}/b$ only, and is plotted in Figure 4. We see that $A^{*}(a_{\mathrm{g}}/b)$ becomes significantly larger than $1$ for $a_{\mathrm{g}}\ll b$ in the case of non-cored potentials, but reaches a maximum of $0.5$ in the isochrone case. For those potentials with finite total mass $M$ we can construct the density-weighted average

[TABLE]

still assuming a circular outer orbit. We find $\langle A^{*}\rangle_{\rho}=0.0617,0.4234$ and $0.65$ for isochrone, Plummer and Hernquist potentials respectively. The isochrone model has by far the smallest $\langle A^{*}\rangle_{\rho}$.

6.2.2 Behavior of $\Gamma$

Figure 5 shows $\Gamma$ in the $(a_{\mathrm{g}},e_{\mathrm{g}})$ plane for the same four potentials. We see that $\Gamma\to 0$ for $a_{\mathrm{g}}\ll b$ in cored potentials. For the coreless potentials $\Gamma$ is always positive, as expected. We see that for $a_{\mathrm{g}}\gtrsim b$, the value of $\Gamma$ is quite sensitive to the outer orbit eccentricity $e_{\mathrm{g}}$ in all four potentials. At fixed $a_{\mathrm{g}}$, increasing $e_{\mathrm{g}}$ corresponds to a decrease in $\Gamma$. For the cored potentials this is because high-$e_{\mathrm{g}}$ orbits start probing the cluster core where the potential is roughly spherical harmonic (for which $\Gamma$ is effectively zero, see §3.1.1 and §6.1), which tends to lower $\Gamma$.

Meanwhile, increasing $a_{\mathrm{g}}$ at fixed $e_{\mathrm{g}}$ tends to increase $\Gamma$. At large $a_{\mathrm{g}}$ all finite mass potentials reduce to a Keplerian potential for which $\Gamma=1$. In the NFW potential, the $\Gamma$ profile is shallow because the potential decays slowly, namely as $\Phi_{\mathrm{NFW}}\sim r^{-1}\ln(r/b)$ for $r\gg b$. Hence it never becomes sufficiently Keplerian to reach $\Gamma\sim 1$.

To better illustrate this convergence at large $a_{\mathrm{g}}/b$, in Figure 6 we show $\Gamma(a_{\mathrm{g}}/b)$ for circular outer orbits $(e_{\mathrm{g}}=0)$ in the same four potentials as in Figure 5. We see that the $\Gamma\to 1$ convergence does occur for all potentials that are asymptotically Keplerian as $r\to\infty$, although in some cases one has to go to radii $a_{\mathrm{g}}\gtrsim 50b$ to observe it satisfactorily. Additionally, at very small radii the NFW density profile can be approximated as a power-law cusp, $\rho\propto r^{-1}$, see equation (54). Using the result listed in §6.1 (and derived in Appendix D.1) we expect to find $\Gamma=1/9$ as $a_{\mathrm{g}}\to 0$, and indeed this is reflected in Figure 6.

We note that some of the orbits in the potentials (50)-(53) will have commensurable (or almost commensurable) radial and azimuthal frequencies. For these orbits, i.e. at some points in $(a_{\mathrm{g}},e_{\mathrm{g}})$ space, equation (39) is not valid, because its derivation relies upon orbits densely filling an axisymmetric annulus, see §4.2. This is particularly true of potentials with a core at small $a_{\mathrm{g}}$, where the potential is close to harmonic (c.f. §7). Nevertheless, Figures 3 and 5 give a good idea of how $A$ and $\Gamma$ change as we consider different orbits within the cluster.

6.3 Behavior of Hamiltonian characteristics in axisymmetric potentials

For axisymmetric potentials it is difficult to make rigorous mathematical statements about $A$ and $\Gamma$. In Appendix D.2 we show that in this case, in principle, $\Gamma$ can take any value $\in(-\infty,\infty)$; however, to achieve extreme negative values of $\Gamma$, or $\Gamma>1$, may require very unusual orbits. (Two such examples are given in Appendix D.2). Meanwhile, in this section we focus on the most typical orbits in axisymmetric potentials via numerical examples. The $A^{*}$ and $\Gamma$ values in this section are calculated numerically using the procedure outlined in Appendix F, and are therefore denoted $A^{*}_{\mathrm{num}},\Gamma_{\mathrm{num}}$.

We will use two axisymmetric potentials in our numerical examples. The first is the flattened power-law potential (Evans, 1994):

[TABLE]

where $-\Phi_{0}$ is the central potential, $b$ is a core radius and $q$ is the oblateness parameter: $q<1$ corresponds to an oblate potential which can be used to model elliptical galaxies and galactic bulges (Evans, 1994). The natural definition of $A^{*}$ in this case is $A^{*}\equiv A/(|\Phi_{0}|/b^{2})$. We choose $\beta=1/2$ and $q=0.94$, meaning that this potential is only slightly flattened. One can derive a number of useful analytical results in such weakly non-spherical potentials; we defer this investigation to a future study. Here we simply demonstrate that even in the case of a weakly flattened potential large departures from the behaviour typical for purely spherical potentials described in §6.2 become possible.

The other potential we will use is the Miyamoto-Nagai potential (Miyamoto & Nagai, 1975):

[TABLE]

where $b_{\ell}$ is the scale length and $b_{h}$ is the scale height. As one changes the value of $b_{h}/b_{\ell}$, the Miyamoto-Nagai potential smoothly transitions from the Kuzmin potential of a razor thin disk ($b_{h}\ll b_{\ell}$) to the spherical Plummer potential frequently used to model globular clusters ($b_{h}\gg b_{\ell}$) (Binney & Tremaine, 2008). In Figure 8 we plot contours of constant density $\rho_{\mathrm{MN}}\equiv\nabla^{2}\Phi_{\mathrm{MN}}/(4\pi G)$ in the $(R,Z)$ plane for two Miyamoto-Nagai models used in this paper, namely $b_{h}/b_{\ell}=1$ and $b_{h}/b_{\ell}=0.05$. The natural definition of $A^{*}$ in this potential is $A^{*}\equiv A/(GM/b_{\ell}^{3})$.

6.3.1 Orbits in the midplane of an axisymmetric potential

The simplest non-spherical case to consider is when the binary’s outer orbit $\mathbf{R}_{\mathrm{g}}$ is confined to the $(X,Y)$ midplane of an axisymmetric potential. Then $\mathbf{R}_{\mathrm{g}}$ still traces a planar rosette with a fixed peri/apocentre $(r_{\mathrm{p}},r_{\mathrm{a}})$ just as in the spherical case, so we can easily compute $A$ and $\Gamma$ as in §4.3. In Appendix C we show how to compute $A,B,\Gamma$ in the case of epicyclic outer orbits in a disk-like potential. We find $A=B=\nu^{2}$, where $\nu$ is the vertical epicyclic frequency at the guiding radius; therefore $\Gamma=1/3$.

In fact, we already deduced in §6.1 that $\Gamma=1/3$ will hold for any orbit $\mathbf{R}_{\mathrm{g}}$ which is confined to the plane of a very thin axisymmetric disk. This follows from the fact that the curvature of the potential is by far greatest in the $Z$ direction at any given position in the disk, so that $\Phi_{zz}(R,0)\gg\Phi_{xx}(R,0)\,\,\forall\,\,R$. Then $A\approx B\approx\overline{\Phi}_{zz}$ and so $\Gamma\approx 1/3$. In Figure 9 we confirm this prediction using three very different orbits in the $(X,Y)$ midplane of a thin ($b_{h}/b_{\ell}=0.05$) Miyamoto-Nagai potential. The $\Gamma_{\mathrm{num}}$ values are (a) $0.321$, (b) $0.330$ and (c) $0.320$, all very close to $\Gamma=1/3$.

6.3.2 Orbits that are far from coplanar

As we show in Appendix D, we always have $\Gamma\geq 0$ in realistic, finite-mass spherical potentials. For $\Gamma$ to fall below zero the potential must be non-spherical, but also, according to definitions (41), the outer orbit of the binary must have $|\overline{\Phi}_{xx}|>|\overline{\Phi}_{zz}|$. Qualitatively, this implies that the average ‘radial curvature’ of the potential over the orbit needs to be greater than the average ‘vertical curvature’. This is not going to be the case while the orbit is confined near a single plane, as we have just seen. However, this situation is naturally realised in potentials that are highly prolate in the $Z$ direction (asymptotically ‘cylindrical’, with $\Phi(R,Z)=\Phi(R)$). In such potentials $\Gamma\approx-1/3$ (see §6.1). Also, to probe the negative $\Gamma$ regime we can consider orbits in non-spherical potentials that make large excursions ‘out of the plane’, i.e. in the $Z$ direction.

This is demonstrated in Figure 10, in which we plot four Orbits (‘IV’-‘VII’) in the flattened power-law potential (56) with $q=0.94$ and $\beta=1/2$. These four Orbits are initiated with exactly the same initial conditions except for their initial azimuthal velocity $v_{\phi}$; the full details of the initial conditions, as well as the resulting $A^{*}_{\mathrm{num}}$ and $\Gamma_{\mathrm{num}}$ values, are given in Table 2. In the top row of Figure 10 we have Orbit (IV), with initial $v_{\phi}=1.35\sqrt{GM/b}$. Orbit (IV) is certainly not planar but the typical excursions in $Z$ are fairly small compared to the excursions in $R$. As a result the $\Gamma_{\mathrm{num}}=0.243$ value is less than $1/3$ but still significantly greater than zero. As we move down the page we decrease the initial azimuthal velocity each time, so that Orbits (V)-(VII) initially have $v_{\phi}/\sqrt{GM/b}=0.95,\,0.35,\,0.05$ respectively (while keeping all other initial conditions the same). The radial excursions decrease as the initial azimuthal velocity decreases, until they become comparable to the vertical excursions. Eventually $\Gamma$ moves below zero in Orbit (VII), see Table 2. The $A^{*}$ values grow as we move down the page since the binary samples a stronger potential when it is closer to the origin.

For our final set of examples we compare four more Orbits (‘VIII’-‘XI’) in the Miyamoto-Nagai potential (57) with $b_{h}/b_{\ell}=1$, the intial conditions for which are given alongside their resulting $A^{*}_{\mathrm{num}},\Gamma_{\mathrm{num}}$ values in Table 3. At large distances this potential is significantly flatter than the $q=0.94$ flattened power-law potential. Orbits (VIII)-(XI) are plotted in the left hand column of Figure 11. Each Orbit has exactly the same initial conditions except we change the initial vertical coordinate $Z$, using $Z=0.1b_{\ell},\,b_{\ell},\,2b_{\ell}$ and $3b_{\ell}$ respectively. Increasing the initial $Z$ thickens the orbit. Orbit (VIII) (top row, initial $Z=0.1b_{\ell}$) is almost coplanar and has $\Gamma_{\mathrm{num}}=0.153$. The $\Gamma_{\mathrm{num}}$ value decreases as move down the page, reaching a minimum value of $\Gamma_{\mathrm{num}}=-0.384$ for Orbit (XI) (bottom row, $Z=3b_{\ell}$) which is thicker vertically than it is radially. Meanwhile, as we move from top to bottom the $A_{\mathrm{num}}^{*}$ value decreases because the orbit spends more time away from the midplane where the tidal potential is strongest. An even more extreme orbit in this potential, with the same initial conditions except $Z=4b_{\ell}$ (resulting in $\Gamma_{\mathrm{num}}=-1.4$), is presented in Figure 12.

Finally, note that Orbit (VIII) is very similar in appearance to the orbit in Figure 9a: both are roughly epicyclic, so they should obey equations (75) and (76). The difference between them is that in the case of Orbit (VIII) the potential felt by the binary is significantly less flattened, since this Orbit resides predominantly in the quasi-spherical core of the potential ($b_{h}/b_{\ell}=1$). As a result, $\Phi_{zz}$ does not dominate over $\Phi_{xx}$ and hence $\Gamma$ ends up being significantly less than $1/3$ (unlike the orbit in Figure 9a which explores much more flattened version of the Miyamoto-Nagai potential with $b_{h}/b_{\ell}=0.05$).

7 Validity of secular Hamiltonian

In §6 we focused on understanding the typical values of $A,\Gamma$ for various types of orbit in different potentials. However the Hamiltonian (42) is only valid if the symmetry conditions (33)-(36) for the time-averages $\overline{\Phi}_{\alpha\beta}$ are satisfied. In addition it is reasonable to require that $\overline{\Phi}_{\alpha\beta}$ converge to fixed values (say to within a few percent) on timescales significantly shorter than the timescale for secular evolution (which will be derived in Paper II). In this section we check the validity of these assumptions numerically. The procedure for calculating $\overline{\Phi}_{\alpha\beta}$ numerically is given in Appendix F.

7.1 Spherical potentials

In a spherical potential, orbits $\mathbf{R}_{\mathrm{g}}(t)$ that have non-commensurable frequencies densely fill an axisymmetric annulus. If this is true then the time-averaged coefficients $\overline{\Phi}_{\alpha\beta}$ should obey the symmetry properties (33), (35).

To verify this we use Orbits (I)-(III). Their initial conditions are given alongside their $A^{*},A^{*}_{\mathrm{num}},\Gamma,\Gamma_{\mathrm{num}}$ values in Table 1. The right panels in Figure 7 show the corresponding running average (from $t=0$ to current time) of numerically computed $\overline{\Phi}_{\alpha\beta}$. As the number of completed orbits grows, the time-averaged derivatives of the potential tend to converge towards fixed values.

Orbit (I) has rather large semi-major axis and small eccentricity, so that it stays far from the core at all times, filling its annulus densely. The ‘axisymmetric annulus’ approximation works very well in this case, so the (semi-)analytic and numerically computed values agree: $A^{*}=A^{*}_{\mathrm{num}}$ and $\Gamma=\Gamma_{\mathrm{num}}$ to within $1\%$ accuracy.

We have picked rather extreme examples in Orbits (II) and (III) in order to demonstrate behaviour of orbits $\mathbf{R}_{\mathrm{g}}$ that are both very radial and tightly bound near the centre of the cluster. Orbit (II) spends a lot of time in the isochrone potential’s constant density core where its frequencies are almost commensurable ($\Omega_{R}\approx 2\Omega_{\phi}$); as a result it precesses slowly, so that there are unfilled gaps left in its annulus even after $t=100T_{\phi}$. This issue does not arise in the uncored Hernquist potential, so Orbit (III) fills its annulus more efficiently than Orbit (II). Nevertheless, the axisymmetric approximation is still very successful in both cases, with a maximum discrepacy of $\sim 4\%$ arising between the $\Gamma$ and $\Gamma_{\mathrm{num}}$ values of Orbit (II).

However, we notice that while the converged symmetry properties of the $\overline{\Phi}_{\alpha\beta}$ (see equations (33), (35)) are well established after $\gtrsim 15T_{\phi}$ for Orbits (I) and (III), they are less well established for Orbit (II) even at $t\gtrsim 45T_{\phi}$. Again this is because Orbit (II) does not fill its annulus efficiently. This can be problematic because if the $\overline{\Phi}_{\alpha\beta}$ fail to converge on a timescale shorter than the secular evolution timescale, the doubly-averaged theory can break down, as we will see in Paper II.

7.2 Axisymmetric potentials

In a (non-spherical) axisymmetric potential, orbits $\mathbf{R}_{\mathrm{g}}(t)$ that have non-commensurable frequencies densely fill an axisymmetric torus and so the time-averaged coefficients $\overline{\Phi}_{\alpha\beta}$ should obey the symmetry properties (33), (35), and (36). To verify this numerically we use Orbits (IV)-(VII) in the flattened power-law potential (56) with $q=0.94$ and $\beta=1/2$ (see Table 2 and Figure 10), and Orbits (VIII)-(XI) in the Miyamoto-Nagai potential (57) with $b_{h}/b_{\ell}=1$ (see Table 3 and Figure 11).

Some features of the $\overline{\Phi}_{\alpha\beta}$ convergence plots are similar to the spherical case. For example, in Figure 10 the derivatives $\overline{\Phi}_{\alpha\beta}$ converge rather slowly in the bottom panel because the Orbit (VII) fills its torus rather sparsely, owing to the large fraction of time it spends in the almost-harmonic potential of the core.

Note that the rightmost columns in Figures 10 and 11 also show the convergence of $\Phi_{xz}$ and $\Phi_{yz}$. This is different from Figure 7 since now we are dealing with non-planar orbits so that these derivatives are no longer identically zero. Although the corresponding time-averages $\overline{\Phi}_{xz},\overline{\Phi}_{yz}$ do indeed converge to zero in all of our axisymmetric examples as expected, in most cases their convergence takes significantly longer than that of the other $\overline{\Phi}_{\alpha\beta}$ coefficients. This is what we would expect from looking at the $\phi_{\mathrm{g}}$ dependence of equations (4.1)-(31): the derivatives $\Phi_{xx}$, $\Phi_{yy}$ and $\Phi_{xy}$ fluctuate twice as rapidly with respect to $\phi_{\mathrm{g}}$ compared to $\Phi_{xz},\Phi_{yz}$. Slower convergence of $\overline{\Phi}_{xz}$ and $\overline{\Phi}_{yz}$ seems to be especially apparent for strongly non-coplanar orbits, i.e. orbits which make large excursions in the $Z$ direction. Orbits that inefficiently fill their torus (i.e. on timescales longer than the secular evolution timescale) can render the doubly-averaged theory inaccurate, as discussed in detail in §7 of Paper II.

8 Discussion

In deriving the secular Hamiltonian (equations (42), (43)) we relied on several approximations. First, we assumed that the outer orbit-averaging procedure used for computing potential derivatives $\overline{\Phi}_{\alpha\beta}$ converges rapidly when compared to the timescale for secular evolution. The rate of convergence of various $\overline{\Phi}_{\alpha\beta}$ components was explored in §7. In Paper II we will use direct numerical integrations of binaries orbiting in stellar cluster potentials to study how well this double averaging procedure works in practice.

Second, we truncated our expansion of the tidal Hamiltonian in §2.1 at the quadrupole order. However, studies of the LK mechanism have shown the importance of higher order — ‘octupole’ — terms for the dynamics of triples in certain situations (Lithwick & Naoz, 2011; Li et al., 2014). This raises a question of whether octupole terms can be important for the secular dynamics of binaries in external tidal fields. While we derive octupole-level corrections to our doubly-averaged Hamiltonian in Appendix E, in practice they are unlikely to be important. This is because in realistic situations the ratio of the semi-major axis of the inner binary orbit ($a\lesssim 100$ AU) is much smaller than the size of its outer orbit ($|\mathbf{R}_{\mathrm{g}}|\sim 1$ pc, comparable to the cluster size), rendering the timescale on which octupole-level effects may manifest themselves too long (we will see in Paper II that a characteristic timescale of secular evolution driven by quadrupole terms in a typical globular cluster is at least tens of Myrs).

Third, our calculation assumes a spatially smooth and time-invariant tidal potential. This approximation neglects the granularity and stochastic variability of the cluster potential caused by encounters with other stars, which are very important in dense environments of clusters (Heggie, 1975; Hut & Bahcall, 1983; Hut, 1983; Heggie & Rasio, 1996; Collins & Sari, 2008). The cumulative effect of a large number of such encounters is what eventually contributes to the smooth tidal field of the cluster (Collins & Sari, 2010); thus, one hopes that in the long run our framework should provide a qualitatively accurate picture of binary evolution in clusters. Nevertheless, in the future we plan to extend our calculation by including the effects of individual stellar encounters on evolution of the binary inner orbit (Weinberg et al., 1987). We will also explore the role of strong encounters (responsible for the formation and disruption of binaries in clusters, Heggie 1975) in resetting the whole course of secular evolution of the binary.

An effect that can modify the binary’s outer orbit in a stellar cluster is resonant relaxation. Rauch & Tremaine (1996) showed that in quasi-Keplerian systems (such as nuclear clusters dominated by a central super-massive black hole) angular momentum is efficiently exchanged between stellar orbits that have commensurable frequencies. When applied to a binary in a quasi-Keplerian cluster, precession of the binary’s outer orbit due to resonant interactions with other stars (so-called ‘vector resonant relaxation’) can alter its inclination relative to the inner orbital plane, potentially bringing an initially low-inclination binary into a high-inclination regime and triggering LK oscillations (Hamers et al., 2018). While this effect has not been explored for binaries in non-Keplerian potentials, vector resonant relaxation can indeed operate in non-Keplerian systems such as globular clusters (Meiron & Kocsis, 2018), as can more general forms of resonant relaxation allowing for exchange of both angular momentum and energy of the outer orbit (Hamilton et al., 2018).

Additionally, the fact that a binary is typically heavier than the average star in a cluster means that it will tend to sink towards the centre of the cluster via dynamical friction (Binney & Tremaine, 2008). Moreover, the global properties of the cluster itself may evolve as a result of two-body relaxation leading to core collapse. All of these effects can change the values of $A$ and $\Gamma$ for a given binary over long time intervals. We will explore their impact upon binary evolution in future work.

Finally, an important assumption that lies at the foundation of our time-averaging procedure is that different frequencies characterizing binary motion in the cluster are not commensurable with each other (see §4.2). If this condition is violated, the outer orbit $\mathbf{R}_{\mathrm{g}}$ no longer fills an axisymmetric torus inside the cluster uniformly in azimuth, rendering the equations (33)-(36) invalid. This issue is addressed in more detail next.

8.1 Commensurable frequencies

Orbits in realistic spherical potentials obey the following relation between the radial ($\Omega_{R}$) and azimuthal ($\Omega_{\phi}$) frequencies (Binney & Tremaine, 2008):

[TABLE]

Thus, in spherical potentials there are infinitely many rational values of $\Omega_{\phi}/\Omega_{R}$ in the interval $(1/2,1)$. However, even the orbits with rational values of $\Omega_{\phi}/\Omega_{R}=p/q$ should still be described (at least roughly) by the ‘filled annulus’ approximation as long as the integers $p,q\gg 1$. For these systems, equation (42) should be approximately valid.

Of course, a small number of low-order resonant points in the $(a_{\mathrm{g}},e_{\mathrm{g}})$ plane will have $\Omega_{\phi}/\Omega_{R}=p/q$ with integer $p,q\sim$ a few. Binaries on these outer orbits will not satisfy the axisymmetric averaging approximation. In particular, $\Omega_{\phi}/\Omega_{R}$ is always rational in the harmonic $(\Omega_{\phi}/\Omega_{R}\equiv 1/2)$ and Keplerian $(\Omega_{\phi}/\Omega_{R}\equiv 1)$ potentials. If the potential is purely harmonic we can treat it as in §3.1.1. Moreover, we show in Appendix B that Keplerian potentials are perfectly described by our doubly-averaged formalism with $\Gamma=1$. However it should be stressed that we have not used the axisymmetric averaging approximation in either of these cases: the harmonic potential just happens to be effectively axisymmetric after single-averaging, and the Keplerian potential is known to be axisymmetric under double averaging to quadrupole order. In neither case do orbits ‘fill their annulus’444The Keplerian potential is not axisymmetric to octupole order, see Appendix E..

Another problematic case is when the binary experiences a potential that is almost harmonic or almost Keplerian, so that the outer orbit precesses apsidally, but not quickly enough to fill a circular annulus on a secular timescale and thereby qualify for an axisymmetric treatment. For example, orbits that spend a lot of time in a constant-density core of the cluster potential experience an almost harmonic potential and so tend to fill their annulus very slowly (see Figure 7c).

8.2 Relation to previous work

Many previous studies have looked at secular evolution of binaries perturbed by external potentials. The effect of an arbitrary quadrupole perturbation upon a binary has been briefly considered by Mikkola & Nurmi (2006). In particular, their equation (20) gives the quadrupole potential experienced by a binary in a star cluster consisting of a large number of point masses $m_{k}$. Our perturbing Hamiltonian $H_{1}$ is recovered from their result in the mean-field limit (i.e. by replacing the exact potential of the cluster, $-G\sum_{k}m_{k}/|\mathbf{R}-\mathbf{R}_{k}|$, with the smooth potential $\Phi(\mathbf{R})$). However, Mikkola & Nurmi (2006) did not explicitly convert to orbital elements, perform any averaging, or develop any secular theory as we do here. In a similar vein, a short paper by Katz & Dong (2011) considered the secular dynamics of a binary perturbed by a generic quadratic potential and included axisymmetric potentials as a special case. They did convert to orbital elements but did not go much further; in particular they did not provide any prescription for computing the coefficients of the averaged perturbing potential.

Studies of tidal effects of the Galactic disk on wide binaries (Heisler & Tremaine, 1986; Byl, 1986; Yabushita, 1989) represent an important limit ($\Gamma\to 1/3$) of our general theory, see §5.2 and Appendix C. Since HT86, Galactic tides have been included in many studies of cometary orbits (e.g. Matese & Whitman 1989; Matese & Whitmire 1996; Breiter et al. 1996; Wiegert & Tremaine 1999; Brasser 2001; Fouchard 2004; Fouchard et al. 2005, 2006; Breiter et al. 2007), as well as planetesimal orbits (e.g. Higuchi et al. 2007).

Veras & Evans (2013a) considered a very general form of the perturbed two-body problem, allowing for both position- and velocity-dependent tidal forces to act upon the binary. Their equations (25)-(29) are more general versions of our singly-averaged equations (c.f. our singly-averaged Hamiltonian (19)), and our equations are recovered if one sets the velocity-dependent forces to zero. However they did not derive any analogues of our doubly-averaged equations. Veras & Evans (2013b) noted that Galactic forces may impact the evolution of exoplanetary systems around stars near the bulge of the Galaxy where the Galactic tide is much stronger than it is in the Solar neighbourhood.

Another interesting and obvious limit of our theory, $\Gamma=1$ — which is, however, rather distinct, see §5.1 and Appendix B — has been explored in numerous studies of Lidov-Kozai dynamics (Lidov, 1962; Kozai, 1962; Fabrycky & Tremaine, 2007; Naoz, 2016) and its extensions. One interesting extension to the LK problem was made by Petrovich & Antonini (2017), who considered the effect of an axisymmetric (non-spherical) nuclear cluster potential on compact-object binaries that are themselves orbiting a central super-massive black hole (SMBH). The non-spherical part of the cluster potential was considered to drive nodal precession of the binary’s quasi-Keplerian outer orbit around the SMBH (continuously changing the relative inclination in the triple system composed of the binary and SMBH, which is important for the operation of LK cycles in this sub-system), while the dominant spherical part drove apsidal precession of the outer orbit. Our doubly-averaged formalism covers this problem in the case where the characteristic timescales for nodal and apsidal precession of the outer orbit are much shorter than the secular timescale, so that the outer orbit fills its torus. Our singly-averaged equations cover it in all cases. However, unlike Petrovich & Antonini (2017), we also account for the direct effect of the tidal torque due to the potential of the cluster on the orbital elements of the inner orbit of the binary.

Several authors have considered the problem of a star in orbit around a SMBH in a nuclear cluster (e.g. Sridhar & Touma 1999; Ivanov et al. 2005; Löckmann et al. 2008; Šubr et al. 2009; Chang 2009; Haas et al. 2011; Merritt 2013; Li et al. 2015; Iwasa & Seto 2016, 2017). The SMBH-star system effectively forms a binary. The binary’s Keplerian orbital elements may then evolve on secular timescales due to some combination of (i) the mean field nuclear cluster potential, (ii) GR pericentre precession, (iii) an infalling massive black hole on a slowly decaying circular orbit, (iv) a circumnuclear ring of material, etc. While this class of problems is reminiscent of our work, they are not quite the same because in the SMBH-star case the barycentre of the binary does not move, and so there is no clean separation between single- and double-averaging. In some cases, e.g. for a binary that sits at the centre of a spherical cusp, there is even no well-defined tidal expansion of the potential. Instead, averaging of the potential is incorporated into the averaging over the stellar Keplerian orbit around the SMBH, which is different from our approach.

One of the most interesting recent applications of secular dynamics has been the possibility of substantial shrinking of binary orbits by LK cycles with dissipative effects. Such applications include the origin of hot Jupiters (Fabrycky & Tremaine, 2007; Naoz et al., 2011; Petrovich, 2015; Hamers, 2017), formation of blue stragglers (Perets & Fabrycky, 2009; Knigge et al., 2009), white-dwarf mergers (Thompson, 2011; Katz et al., 2011; Toonen et al., 2018), and compact-object binary mergers in globular or nuclear star clusters (Antognini et al., 2014; Rodriguez et al., 2015; Naoz, 2016; Silsbee & Tremaine, 2017; Petrovich & Antonini, 2017; Leigh et al., 2018). Binary evolution driven by cluster tides explored in our work represents a different evolutionary scenario that may lead to similar outcomes (without invoking a nearby third companion). This possibility is explored in a separate study (Hamilton & Rafikov, 2019b).

9 Summary

This work explores secular evolution of binary systems orbiting in axisymmetric stellar clusters. We derive a Hamiltonian describing this evolution for an arbitrary form of the smooth cluster potential, average it over the (inner) orbital motion of the binary, and then average it again over the (outer) orbit of the binary around the cluster assuming the potential is axisymmetric. Our results can be summarized as follows.

• •

When the doubly-averaged Hamiltonian is cast in dimensionless form, all the information about the tidal potential is contained in a single dimensionless parameter $\Gamma$, which depends on the background potential $\Phi$ and the orbit of the binary in the cluster. The value of this parameter determines the phase portrait of the binary evolution, which is explored in a companion study (Paper II).

• •

The timescale of secular evolution is set by another (dimensional) parameter $A$, which, like $\Gamma$, depends on the cluster potential and the binary’s outer orbit.

• •

In certain cases $A$ and $\Gamma$ can be calculated (semi-)analytically. Such cases include (a) orbits in spherical potentials, (b) orbits confined to the midplane of an axisymmetric potential, and (c) epicyclic orbits near the midplane of an axisymmetric potential. We demonstrate how our calculations reproduce the known results for Lidov-Kozai evolution, evolution of Oort Cloud comets due to the Galactic tide, and so on.

• •

We map out the behavior of $A$ and $\Gamma$ in different spherically symmetric potentials as a function of size and radial elongation of the binary orbit. We find that $\Gamma$ is small in the central regions of clusters with cored potentials, but tends to unity in clusters with finite mass as the orbit size increases. In general, $0\leq\Gamma\leq 1$ in realistic finite-mass spherical potentials.

• •

In general axisymmetric potentials, $\Gamma$ can easily attain negative values, in particular for highly inclined (i.e. non-coplanar) orbits.

• •

The accuracy with which our doubly-averaged Hamiltonian characterizes binary evolution deteriorates for highly non-coplanar orbits in axisymmetric potentials. Commeurability of orbital frequencies in the cluster potential may also present a problem for application of our theory at a quantitative level.

These results will be extensively used in Paper II, where we systematically explore the dynamics of binaries driven by the tidal field of clusters for different values of $\Gamma$. There we verify numerically the predictions of the secular theory based on our doubly-averaged Hamiltonian, and derive timescales for secular eccentricity oscillations. In future papers in this series, these calculations will be applied to understanding secular evolution of binaries in dense clusters and exploring their relevance for the formation of compact-object mergers (Hamilton & Rafikov, 2019b), blue stragglers, and so on.

Acknowledgements

We thank John Magorrian for pointing out some literature and Eugene Vasiliev, Fabio Antonini and Almog Yalinewich for useful comments. The numerical orbit integration was done with galpy (http://github.com/jobovy/galpy, Bovy 2015). CH is funded by a Science and Technology Facilities Council (STFC) studentship.

Appendix A The singly-averaged Hamiltonian in orbital elements

In this Appendix we give the full algebraic expressions in terms of orbital elements for the terms $\langle r_{\alpha}r_{\beta}\rangle_{M}$ that enter the singly-averaged Hamiltonian (19). They are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Appendix B Recovering the Lidov-Kozai quadrupole Hamiltonian

To derive the LK Hamiltonian we take equation (19) and average it over time using $\Phi(r)=-GM/r$, and with $\mathbf{R}_{\mathrm{g}}(t)$ describing a Keplerian ellipse with the focus at the origin.

First, it is obvious that for the this potential, $\Phi^{\prime\prime}(R_{\mathrm{g}})\propto\Phi^{\prime}(R_{\mathrm{g}})/R_{\mathrm{g}}\propto R_{\mathrm{g}}^{-3}$. Then we must average the right hand sides of equations (4.1)-(29) which requires that we average the quantities $R_{\mathrm{g}}^{-3},R_{\mathrm{g}}^{-3}\cos 2\phi_{\mathrm{g}}$ and $R_{\mathrm{g}}^{-3}\sin 2\phi_{\mathrm{g}}$.

Without loss of generality we may choose the argument of pericentre $\omega_{\mathrm{g}}$ of the ellipse to be zero, so that $\phi_{\mathrm{g}}=f_{\mathrm{g}}$, the true anomaly. Then $R_{\mathrm{g}}=a_{\mathrm{g}}(1-e_{\mathrm{g}}\cos E_{\mathrm{g}})$ with $E_{\mathrm{g}}$ being the eccentric anomaly, so for an arbitrary function $\mathcal{S}(\phi_{\mathrm{g}})$ we can write

[TABLE]

and we can convert between $f_{\mathrm{g}}$ and $E_{\mathrm{g}}$ using

[TABLE]

The answers are

[TABLE]

Using these identities, we find a remarkably simple relation between the time-averaged coefficients:

[TABLE]

and as expected $\overline{\Phi}_{xy}=\overline{\Phi}_{xz}=\overline{\Phi}_{yz}=0$. Note that the regime $\overline{\Phi}_{xx}=\overline{\Phi}_{yy}$ is exactly that of an axisymmetric perturbing potential (see e.g. Katz et al. (2011)). The resulting perturbing Hamiltonian will therefore be the same as (42), with the added simplification that $A=(GM/2)a_{\mathrm{g}}^{-3}(1-e_{\mathrm{g}}^{2})^{-3/2}$ and $\Gamma=1$. Making these substitutions we find

[TABLE]

This is precisely the dimensionless test particle quadrupole Lidov-Kozai Hamiltonian (Lidov, 1962; Kozai, 1962; Kinoshita & Nakai, 2007; Lithwick & Naoz, 2011; Antognini et al., 2014). It describes the secular evolution of a hierarchical triple system in which the outer orbit dominates the angular momentum budget.

Appendix C Epicyclic orbits

In this Appendix we look at the behavior of $A$ and $\Gamma$ in the case of a binary performing epicyclic motion in an axisymmetric disk, to connect with the results of Heisler & Tremaine (1986) (HT86), who calculated the secular effect of the Galactic tide on the Oort Cloud comets.

Let the guiding centre of the binary’s orbit be a circle of radius $R_{\mathrm{c}}$ in the $Z=0$ plane of the potential. The potential experienced by the binary can then be approximated as555We implicitly assume that the disk is symmetric about its midplane, $\Phi(R,Z)=\Phi(R,-Z)$, so that $\partial^{2}\Phi/\partial R\partial Z=0$ at $Z=0$. Otherwise the binary would not remain in the midplane anyway.

[TABLE]

Using this expression and equations (33)-(34) it is easy to show that

[TABLE]

while all other $\overline{\Phi}_{\alpha\beta}=0$; here

[TABLE]

is the angular frequency of the guiding centre, while

[TABLE]

are the radial and vertical epicyclic frequencies of $\mathbf{R}_{\mathrm{g}}$ respectively. Hence

[TABLE]

Near the midplane of a galactic disk, and in particular in the HT86 case of the solar neighbourhood of the Milky Way, it is almost always the case that $\Omega_{\mathrm{c}}\sim\kappa\ll\nu$. Thus to a very good approximation $A=\nu^{2}$ and $\Gamma=1/3$. Plugging these results into our doubly-averaged Hamiltonian (42), written in Delaunay variables, we find

[TABLE]

where we have eliminated $\nu$ in favour of the density in the Solar neighbourhood $\rho_{0}$ using Poisson’s equation $4\pi G\rho_{0}=(\nabla^{2}\Phi)_{\mathbf{R}_{\mathrm{g}}}\approx(\partial^{2}\Phi/\partial z^{2})_{\mathbf{R}_{\mathrm{g}}}\equiv\nu^{2}$. This is precisely the Hamiltonian arrived at by HT86 (c.f. their equation (14)) when considering the effect of the Galactic tide on the Oort Cloud comets.

Appendix D Signs and sizes of $A$ and $\Gamma$

Here we provide some technical details about the statements on the signs and values of $A$ and $\Gamma$ made in §6.1. Also, in Table 4 we summarize some information about these coefficients for certain potentials. We provide two examples of orbits with extreme values of $\Gamma$.

D.1 Spherical potentials

Consider a spherically symmetric potential $\Phi=\Phi(r)$, where $r\equiv\sqrt{R^{2}+Z^{2}}$ is the sphericl radius. According to our convention, the outer orbit of the binary always lies in $Z=0$ plane of the associated cylindrical $(R,\phi,Z)$ coordinate system. Then it is a simple matter to show that

[TABLE]

Using these conversions as well as equations (33), (34), (41) we find

[TABLE]

We can now prove that $A>0$ for (almost) any realistic spherical potential, and thereby show that for such systems $0\leq\Gamma\leq 1$. In a spherical potential $\Phi(r)$ we have

[TABLE]

where $M(r)$ is the cluster mass enclosed inside radius $r$. Also, Poisson’s equation reads

[TABLE]

allowing us to rewrite equation (80) as

[TABLE]

Since $\rho>0$ and $M>0$ at all radii, this inevitably results in $A>0$.

Finally, since

[TABLE]

equation (81) can be rewritten as

[TABLE]

For any spherical system in which the density is a non-increasing function of radius, $\mathrm{d}(M/r^{3})/\mathrm{d}r\leq 0$ for any $r$ and hence $B\geq 0$.

Let us now focus on spherical systems with $\mathrm{d}\rho/\mathrm{d}r\leq 0$. If the cluster has a constant density core, then $M(r)\sim r^{3}$ as $r\to 0$ and so $B\to 0$ (equation (86)). Hence if $\mathbf{R}_{\mathrm{g}}$ orbits entirely inside the constant density region, $\Gamma=0$. A potential without a core will always have a non-zero value of $\Gamma$ for orbits at small radii.

At the other extreme, as $r\to\infty$ we have $\rho\to 0$, and usually the enclosed mass $M(r)\to$ const (although see below for potentials arising from power-law cusp density profiles). Hence $\mathrm{d}(M(r)/r^{3})/\mathrm{d}r\to-3M/r^{4}$ where $M$ is the total mass of the cluster, and in this limit we get $B=3A$. Thus an orbit $\mathbf{R}_{\mathrm{g}}$ that spends its time exclusively at very large radii $r$ compared to the scale radius of the cluster will have $\Gamma\to 1$. This is precisely the Lidov-Kozai limit: for potentials that are Keplerian as $r\to\infty$ (i.e. those with finite mass), orbiting far from the core is equivalent to orbiting a point mass at the origin.

Finally, for any orbit in a spherical cluster with a power-law density cusp $\rho(r)\propto r^{-\beta}$ with $0<\beta<3$ (so that the mass is finite at the centre) one naturally has $M(r)=4\pi\rho r^{3}/(3-\beta)$. Then we find from equations (84), (86) that

[TABLE]

and so $U=1/(3-\beta)$ and $\Gamma=\beta/[3(4-\beta)]$.

D.2 Axisymmetric potentials

In a general axisymmetric potential there is no constraint on how negative the parameter $U$ (defined by equation (46) and plotted in Figure 2) can be. Non-spherical potentials naturally feature regions with $\Phi_{zz}<0$, especially near the poles; choosing a highly inclined (with respect to the equatorial plane of the potential) orbit with large radius so that $\overline{\rho}$ is vanishingly small, one can drive strongly negative $U$, thereby achieving extreme (positive or negative) values of $\Gamma$.

In Figure 12 we give an example of an orbit with $\Gamma_{\mathrm{num}}=-1.4$. We use the Miyamoto-Nagai potential with $b_{h}/b_{\ell}=1$ and precisely the same initial conditions as Orbits (VIII)-(XI) in the main text (see Table 3), except that we now take the initial $Z$ coordinate to be $5.0b_{\ell}$. All three panels show $100T_{\phi}$ of data. The large initial value of $Z$ means that the orbit spends a lot of time near the poles of the potential where $\Phi_{zz}<0$.

In Figure 13 we provide an example of a very polar orbit resulting in $\Gamma_{\mathrm{num}}=2.0$ ($U=-2.3$). We use the Miyamoto-Nagai potential $b_{h}/b_{\ell}=0.1$, and initial conditions

[TABLE]

In Figure 13a we display only the first $5$ vertical periods $T_{Z}$ of the orbit in the meridional $(R,Z)$ plane. When integrated for a long time, the orbit remains almost polar but precesses very slowly around the $Z$ axis until it eventually fills an axisymmetric torus after a few thousand $T_{Z}$. Figure 13a shows that the convergence of the $\overline{\Phi}_{\alpha\beta}$ coefficients in this case is very slow and takes $\gtrsim 5000T_{Z}$. In practice, unless the binary is very tight, secular theory is unlikely to work well for such an orbit. Indeed, for a relatively wide binary the secular evolution timescale is likely to be much shorter than $5000T_{Z}$, meaning that our assumption that the binary fills its torus (and hence the $\Phi_{\alpha\beta}$ converge) in much less than a secular timescale is violated. We will explore this issue in more detail in Paper II, §7.

Appendix E Octupole Hamiltonian

The Hamiltonian derived in §2 is correct to quadratic order in $a/|\mathbf{R}_{\mathrm{g}}|$, the so-called ‘quadrupole approximation’. We can attempt to derive a more accurate Hamiltonian by keeping the higher order terms in the series expansion of equation (3). The next (‘octupole’) term that we would include is

[TABLE]

We then use the fact that $\mathbf{r}\equiv\mathbf{r}_{1}-\mathbf{r}_{2}$ and $m_{1}\mathbf{r}_{1}+m_{2}\mathbf{r}_{2}=\mathbf{0}$ to write

[TABLE]

As a result, the equation for the relative motion $\mathrm{d}^{2}\mathbf{r}/\mathrm{d}t^{2}$ (equation (5)) can be written purely in terms of $\mathbf{r}$. The next order (‘octupole’) correction to the Hamiltonian (42) is

[TABLE]

Note that the octupole term vanishes for equal-mass binaries ($m_{1}=m_{2}$).

The corresponding doubly-averaged perturbing octupole term is then simply

[TABLE]

E.1 Time-averaging over an axisymmetric torus

As in the quadrupole case, it turns out that when time-averaged over an axisymmetric torus the coefficents $\overline{\Phi}_{\alpha\beta\gamma}$ satisfy various symmetry properties. After a little algebra one can show that

[TABLE]

and all other $\overline{\Phi}_{\alpha\beta\gamma}=0$ except for $\overline{\Phi}_{zzz}$. Hence

[TABLE]

Writing out the $\langle.\rangle_{M}$ factors in terms of orbital elements we have

[TABLE]

[TABLE]

Equations (94)-(96) provide a general framework for accounting for the octupole contribution to the tidal Hamiltonian in an arbitrary axisymmetric potential. Note there is no $\Omega$ dependence in the octupole Hamiltonian, so $J_{z}$ is still conserved to octupole order.

E.2 Link to the test particle octupole LK Hamiltonian

Note that one cannot recover the test particle octupole term of the doubly-averaged Lidov-Kozai Hamiltonian by putting $\Phi=-GM/r$ in (94), because equation (94) is derived under the axisymmetric approximation. The time-averaged potential of a perturber on a Keplerian orbit is only axisymmetric at the quadrupole level, and the symmetry is broken by octupole terms. Instead one must integrate over the outer Keplerian orbit exactly, as in Appendix B.

In general there are $10$ independent time-averaged coefficients $\overline{\Phi}_{\alpha\beta\gamma}$ to consider. We choose the outer orbit to be in the $Z=0$ plane so we can immediately eliminate four of these, $\Phi_{zzz}=\Phi_{xxz}=\Phi_{yyz}=\Phi_{xyz}=0$. We can also choose the pericentre of the outer orbit to be on the $X$ axis without loss of generality, so that the ellipse traced by the outer orbit is symmetric under $Y\to-Y$. Then all $\Phi_{\alpha\beta\gamma}$ that contain an odd number of $Y$ derivatives will be antisymmetric under $Y\to-Y$, so their time-averages over this ellipse will vanish: $\overline{\Phi}_{xxy}=\overline{\Phi}_{zzy}=\overline{\Phi}_{yyy}=0$. This leaves us with only three non-zero terms in the doubly-averaged octupole LK Hamiltonian:

[TABLE]

For reference we now write down the terms that make up equation (97). First we write down the necessary $\Phi_{\alpha\beta\gamma}$ coefficients in terms of cylindrical coordinates $R_{\mathrm{g}}$ and $\phi_{\mathrm{g}}$:

[TABLE]

The time-averages of these coefficients are

[TABLE]

The mean-anomaly averaged quantities are

[TABLE]

[TABLE]

[TABLE]

Plugging the results (E.2)-(101) in to (97), the resulting doubly-averaged test particle octupole Lidov-Kozai Hamiltonian is

[TABLE]

where $\theta\equiv\cos i$. Equation (107) is precisely the result found in standard LK literature (e.g. Ford et al. 2000; Lithwick & Naoz 2011; Naoz 2016).

Appendix F Numerical prescription for computing time-averages

To calculate the time-averages $\overline{\Phi}_{\alpha\beta}$ numerically we use the orbit integrator in galpy (Bovy, 2015). Given the initial position $\mathbf{R}_{\mathrm{g}}(0)$ and velocity of the binary’s outer orbit around the cluster, we integrate its equation of motion (11) numerically in the smooth cluster potential $\Phi$. We use a constant timestep $\Delta t$ so that after $k$ timesteps the time elapsed is $t_{k}=k\Delta t$. Then the running time-average of a quantity $\mathcal{F}(\mathbf{R}_{\mathrm{g}}(t))$ is

[TABLE]

In nearly all numerical examples shown in this paper we used $\Delta t\approx T_{\phi}/100$ where $T_{\phi}$ is the azimuthal period of $\mathbf{R}_{\mathrm{g}}$, and integrated the outer orbit for approximately $100T_{\phi}$. The exception is Figure 13, where we used $\Delta t\approx T_{Z}/10$ ($T_{Z}$ is the vertical period of $\mathbf{R}_{\mathrm{g}}$) and integrated the outer orbit for $5000T_{Z}$.