Secular dynamics of binaries in stellar clusters II: dynamical evolution

TL;DR

This paper explores the secular evolution of binaries in stellar clusters using a Hamiltonian framework, revealing four dynamical regimes influenced by the cluster's tidal field parameter, with implications for exotic stellar objects.

Contribution

It provides a detailed phase-space analysis of binary dynamics in clusters, identifying bifurcations and regimes, and validates the formalism with numerical simulations.

Findings

Phase-space structure varies with parameter $\Gamma$

Four distinct dynamical regimes identified

Formalism accurately predicts secular evolution timescales

Abstract

Dense stellar clusters are natural sites for the origin and evolution of exotic objects such as relativistic binaries (potential gravitational wave sources), blue stragglers, etc. We investigate the secular dynamics of a binary system driven by the global tidal field of an axisymmetric stellar cluster in which the binary orbits. In a companion paper (Hamilton & Rafikov 2019a) we developed a general Hamiltonian framework describing such systems. The effective (doubly-averaged) Hamiltonian derived there encapsulates all information about the tidal potential experienced by the binary in its orbit around the cluster in a single parameter $\Gamma$. Here we provide a thorough exploration of the phase-space of the corresponding secular problem as $\Gamma$ is varied. We find that for $\Gamma > 1/5$ the phase-space structure and the evolution of binary orbital element are qualitatively similar to the Lidov-Kozai problem. However, this is only one of four possible regimes, because the dynamics are qualitatively changed by bifurcations at $\Gamma = 1/5,0,-1/5$. We show how the dynamics are altered in each regime and calculate characteristics such as secular evolution timescale, maximum possible eccentricity, etc. We verify the predictions of our doubly-averaged formalism numerically and find it to be very accurate when its underlying assumptions are fulfilled, typically meaning that the secular timescale should exceed the period of the binary around the cluster by $\gtrsim 10-10^2$ (depending on the cluster potential and binary orbit). Our results may be relevant for understanding the nature of a variety of exotic systems harboured by stellar clusters.

Figures

Tables

Table 1. Table 1: Summary of ( ω , j ) 𝜔 𝑗 (\omega,j) locations and values of minima and maxima of j 𝑗 j , values of Δ Δ \Delta , and Θ Θ \Theta and D 𝐷 D ranges, for different types of orbit in each Γ Γ \Gamma regime. These determine minimum and maximum eccentricities as well as timescales of secular eccentricity oscillations via equation ( 33 ).

Type of orbit	Regime	$(\omega ,j_{\min })$	$(\omega ,j_{\max })$	$\Delta$	$\Theta \in$	$D\in$
Librating, $\Gamma >1/5$		$(\pm \pi /2,j_{-})$	$(\pm \pi /2,j_{+})$	$j_0^2-j_{-}^2$	$(0,\frac{1+5\Gamma }{10\Gamma })$	$(D_{-},0)$
Circulating, $\Gamma >1/5$		$(\pm \pi /2,j_{-})$	$(0,j_0)$	$j_{+}^2-j_{-}^2$	$(0,1)$	$(0,1-\Theta )$
Librating,		$(\pm \pi /2,j_{-})$	$(\pm \pi /2,j_{+})$	$j_{+}^2-j_0^2$	$(0,\frac{10\Gamma }{1+5\Gamma })$	$(1-\Theta ,D_{-})$
Circulating,		$(0,j_0)$	$(\pm \pi /2,j_{+})$	$j_{+}^2-j_{-}^2$	$(0,1)$	$(0,1-\Theta )$
		$(\pm \pi /2,j_{+})$	$(0,j_0)$	$j_0^2-j_{-}^2$	$(0,1)$	$(0,1-\Theta )$
Librating, $\Gamma \le -1/5$		$(\pm \pi /2,j_{+})$	$(\pm \pi /2,j_{-})$	$j_0^2-j_{+}^2$	$(0,\frac{1+5\Gamma }{10\Gamma })$	$(D_{+},0)$
Circulating, $\Gamma \le -1/5$		$(\pm \pi /2,j_{+})$	$(0,j_0)$	$j_{-}^2-j_{+}^2$	$(0,1)$	$(0,1-\Theta )$

Table 2. Table 2: Summary of ( ω , j ) 𝜔 𝑗 (\omega,j) locations and values of minima and maxima of H 1 ∗ superscript subscript 𝐻 1 H_{1}^{*} , depending on Γ Γ \Gamma and whether fixed points exist.

Category	$H_1^{\ast }$ min location	$H_{1,\min }^{\ast }$	$H_1^{\ast }$ max location	$H_{1,\max }^{\ast }$
$\Gamma >1/5$, fixed points exist	$j=\sqrt{\Theta }$	$\left(5-3\Theta \right)\left(1-3\Gamma \right)$	$(\pm \pi /2,j_f)$	$H_{-}$
$\Gamma >1/5$, fixed points do not exist	$j=\sqrt{\Theta }$	$\left(5-3\Theta \right)\left(1-3\Gamma \right)$	$j=1$	$2\left(1-3\Theta \Gamma \right)$
, fixed points exist	$j=1$	$2\left(1-3\Theta \Gamma \right)$	$(\pm \pi /2,j_f)$	$H_{-}$
, fixed points do not exist	$j=1$	$2\left(1-3\Theta \Gamma \right)$	$j=\sqrt{\Theta }$	$\left(5-3\Theta \right)\left(1-3\Gamma \right)$
	$j=1$	$2\left(1-3\Theta \Gamma \right)$	$j=\sqrt{\Theta }$	$\left(5-3\Theta \right)\left(1-3\Gamma \right)$
$\Gamma \le -1/5$, fixed points exist	$(\pm \pi /2,j_f)$	$H_{+}$	$j=\sqrt{\Theta }$	$\left(5-3\Theta \right)\left(1-3\Gamma \right)$
$\Gamma \le -1/5$, fixed points do not exist	$j=1$	$2\left(1-3\Theta \Gamma \right)$	$j=\sqrt{\Theta }$	$\left(5-3\Theta \right)\left(1-3\Gamma \right)$

Equations

$${\overline{⟨H_1⟩}}_M=C{H}_1^{\ast },\mathrm{where}C=\frac{A{a}^2}{8}$$

$$H_1^{\ast }$$

$$\Gamma$$

$$0<$$

$$-1/5<$$

$$\Gamma$$

$$\Theta \equiv J_z^2/L^2=(1-e^2){\cos }^{2}i,j\equiv J/L=\sqrt{1-e^2},$$

$$H_1^{\ast }=$$

$$-15\Gamma (j^2-\Theta )(1-j^2)\cos 2\omega ].$$

$$0\le e\le \sqrt{1-\Theta },{\Theta }^{1/2}\le j\le 1,$$

$$\frac{\mathrm{d}\omega }{\mathrm{d}t}$$

$$\frac{\mathrm{d}j}{\mathrm{d}t}$$

$$j_{\mathrm{f}}(\Theta ,\Gamma )={\left(\frac{\Theta }{\Lambda }\right)}^{1/4},\Lambda (\Gamma )=\frac{5\Gamma +1}{10\Gamma }.$$

$$\Theta <\min \left(\Lambda ,{\Lambda }^{-1}\right).$$

$$j^2(\omega =0)=\frac{5/3-2\Gamma \Theta -5\Gamma -H_1^{\ast }/3}{1-5\Gamma }.$$

$$D\equiv 1-j^2(\omega =0)=\frac{H_1^{\ast }/3-2/3+2\Gamma \Theta }{1-5\Gamma }.$$

$$D=e^2\left(1+\frac{10\Gamma }{1-5\Gamma }{\sin }^{2}i{\sin }^2\omega \right).$$

$$j_0\equiv j(\omega =0)=\sqrt{1-D}.$$

$$j_{\pm }\equiv j(\omega =\pi /2)=\sqrt{\frac{\Sigma \pm \sqrt{{\Sigma }^2-10\Gamma \Theta \left(1+5\Gamma \right)}}{1+5\Gamma }},$$

$$\Sigma$$

$$=\frac{1+5\Gamma }{2}(1-e^2)+5\Gamma \left({\cos }^{2}i+e^2{\sin }^{2}i{\cos }^2\omega \right).$$

$$H_1^{\ast }(j=1)=2(1-3\Gamma \Theta ),$$

$$H_1^{\ast }(j=\Theta )=(5-3\Theta )(1-3\Gamma ),$$

$$H_1^{\ast }(\omega =\pi /2,j=j_{\mathrm{f}})=\{\begin{array}{l}{H}_{-},\Gamma >0,\\ H_{+},\Gamma \le -1/5,\end{array}$$

$$H_{\pm }=5(1+3\Gamma )+24\Gamma \Theta \pm 6\sqrt{10\Gamma \Theta (1+5\Gamma )}.$$

$$D(j=1)=0,$$

$$D(j=\Theta )=1-\Theta ,$$

$$D(\omega =\pi /2,j=j_{\mathrm{f}})=\{\begin{array}{l}{D}_{-},\Gamma >0,\\ D_{+},\Gamma \le -1/5,\end{array}$$

$$D_{\pm }=\frac{1}{(1-5\Gamma )}\left[1+5\Gamma +10\Gamma \Theta \pm 2\sqrt{10\Gamma \Theta \left(1+5\Gamma \right)}\right].$$

$$\cos 2\omega =\frac{(j^2-3\Gamma \Theta )(5-3j^2)-j^2{H}_1^{\ast }}{15\Gamma (j^2-\Theta )(1-j^2)}.$$

$$\frac{\mathrm{d}{j}^2}{\mathrm{d}t}=\pm \frac{12C}{L}\sqrt{(25{\Gamma }^2-1)(j_0^2-j^2)(j_{+}^2-j^2)(j^2-j_{-}^2)}.$$

$$t_{\sec }=2{\int }_{j_{\min }^2}^{j_{\max }^2}{\left(\frac{\mathrm{d}{j}^2}{\mathrm{d}t}\right)}^{-1}\mathrm{d}{j}^2,$$

$$\Delta \equiv \max [j_{+}^2,j_{-}^2,j_0^2]-\min [j_{+}^2,j_{-}^2,j_0^2],$$

$$t_{\sec }=\frac{t_1}{\sqrt{|1-25{\Gamma }^2|\Delta }}K\left(\sqrt{\frac{j_{\max }^2-j_{\min }^2}{\Delta }}\right),$$

$$t_1$$

Full text

††pubyear: 2018††pagerange: Secular dynamics of binaries in stellar clusters II: dynamical evolution–B.4

1 Introduction

Orbital evolution of binary systems in dense stellar clusters is a very rich problem. Historically, the focus in this area has been on the cumulative effect of multiple encounters of the binary with the other stellar members of the cluster (Heggie, 1975). In particular, the excitation of binary eccentricity via this mechanism was explored (Heggie & Rasio, 1996) as a means of understanding the properties of binary pulsars in globular clusters (Phinney, 1992; Phinney & Kulkarni, 1994). At the same time, the effect of the smooth component of the cluster gravitational field on the internal dynamics of the binary was largely overlooked in the past (although it was routinely accounted for in the studies of the Oort Cloud dynamics and wide binaries driven by the Galactic tide, see e.g. Heisler & Tremaine 1986; Brasser 2001; Veras & Evans 2013).

To fill this gap, in Hamilton & Rafikov (2019a) — hereafter ‘Paper I’ — we explored secular evolution of a binary system driven by the smooth tidal field of a cluster in which the binary orbits. We developed a general Hamiltonian framework which describes the evolution of the orbital elements of a binary fully accounting for the details of the (axisymmetric or spherical) cluster potential and the binary orbit in it. This was done by expanding the cluster potential near the barycentre111To sufficient accuracy, the binary’s barycentre follows the trajectory of a test particle, $\mathrm{d}^{2}\mathbf{R}_{\mathrm{g}}/\mathrm{d}t^{2}=-\nabla\Phi$ where $\Phi$ is the smooth cluster potential. $\mathbf{R}_{\mathrm{g}}$ of the binary consisting of two point masses $m_{1}$ and $m_{2}$ and performing two averages of the resultant tidal perturbation: first over the binary’s ‘inner’ orbital period (i.e. over its mean anomaly), and second over many ‘outer’ orbits of the binary itself around the cluster. As a result, defining the binary orbital elements — semi-major axis $a$, eccentricity $e$, inclination $i$, longitude of the ascending node $\Omega$, argument of pericentre $\omega$ and mean anomaly $M$ — with respect to a frame tied to the equatorial plane of the axisymmetric cluster (or the plane of the outer orbit of the binary in the case of a spherical cluster), Paper I arrived at a secular (‘doubly-averaged’) perturbing Hamiltonian

[TABLE]

is a constant with dimensions of energy per unit mass and $H_{1}^{*}$ is the ‘dimensionless Hamiltonian’

[TABLE]

The quantities $A$ (with units of $($frequency$)^{2}$) and $\Gamma$ (dimensionless) emerging in these expressions are constants, which encode all of the information about the time-averaged tidal potential experienced by the binary as it moves around the cluster. Paper I explored in detail the dependence of $A$ and $\Gamma$ upon the shape of the background potential and the binary’s orbit within it.

The quantity $A$ is an overall proportionality constant and hence sets the timescale for the secular eccentricity oscillations. The dimensionless part of the Hamiltonian, $H_{1}^{*}$, depends on the potential and the binary’s outer orbit only through the constant $\Gamma$ (see equation 2). The value of $\Gamma$ determines the phase-space morphology for the dynamics of the inner orbit of the binary, as we will see in this work. In particular, in Paper I we found that $\Gamma\to 1$ is reached when the cluster potential is Keplerian, in which case $H_{1}^{*}$ reduces to the test particle quadrupole Lidov-Kozai (hereafter ‘LK’) Hamiltonian (Lidov, 1962; Kozai, 1962; Antognini, 2015; Naoz, 2016). In addition, the Hamiltonian for binaries orbiting in a thin Galactic disk (first derived by Heisler & Tremaine 1986) was recovered from equation (2) when $\Gamma=1/3$. It was also shown that for binaries in realistic spherical potentials $0\leq\Gamma\leq 1$, while in non-spherical potentials $\Gamma$ can be (but is not necessarily) negative.

The goal of the present paper is to systematically explore the dynamics that result from the general secular theory based on the Hamiltonian (1)-(2). It turns out that the dynamics in the range $\Gamma>1/5$ are qualitatively very similar to those arising in the LK case, but that bifurcations occur when $\Gamma=1/5$, [math] and $-1/5$, which change the picture significantly. As a result, in this work we separately treat four distinct dynamical regimes:

[TABLE]

We recognise that our ‘clean’, idealised scenario in which binaries are torqued by an entirely smooth, time-independent potential ignores various effects which must be included if our theory is to be of astrophysical relevance. The most obvious of these is perturbations of the binary by stellar flybys. However for clarity we choose to focus exclusively on the idealised problem for most of the paper, before discussing non-ideal effects at the end, including flybys (§9.3).

In §2 we derive general results that hold for all $\Gamma$: conditions for the existence of fixed points (§2.2), the criteria for phase-space trajectories to librate or circulate (§2.3), the values of maximum and minimum eccentricities (§2.4), the timescale of eccentricity oscillations (§2.6), etc. Then, in §§3-6, we explore the details of each of the $\Gamma$ regimes (3)-(6) separately. The validity of the doubly-averaged (secular) theory is verified numerically in §7. The impact of general relativistic (GR) pericentre precession on the cluster tide-driven secular evolution is explored in §8. We collect our results in §9, discuss them in light of the existing literature and comment on the applicability of our formalism. We summarise our findings and discuss future applications in §10.

2 General aspects of secular dynamics

Our goal is to understand evolution of the orbital elements of the inner orbit of the binary — eccentricity, inclination, etc. — as it moves in the cluster potential. We do this by carefully investigating the phase-space of the dimensionless Hamiltonian $H_{1}^{*}$ (equation (2)).

It will prove useful to express the Hamiltonian in terms of Delaunay variables. We therefore define the actions $L=\sqrt{G(m_{1}+m_{2})a},J=L\sqrt{1-e^{2}}$, and $J_{z}=J\cos i$, their corresponding angles being the mean anomaly $M$, the argument of pericentre $\omega$ and the longitude of ascending node $\Omega$ respectively. (As in Paper I, the symmetry axis of the axisymmetric cluster is chosen to be the $Z$ axis. In the case of a spherical cluster, the $Z$ axis is chosen perpendicular to the plane of the binary’s outer orbit).

The Hamiltonian (1) does not depend on the mean anomaly $M$, so the action $L$ is conserved. Hence we can choose to work with dimensionless versions of our variables $(J,J_{z})$; following the notation of Antognini (2015) we define

[TABLE]

and we clearly have $0\leq\Theta\leq j^{2}\leq 1$. Obviously $j$ is just the dimensionless angular momentum. Then

[TABLE]

Since $H_{1}^{*}$ is independent of the longitude of the ascending node $\Omega$, the dimensionless quantity $\Theta$ is also an integral of motion (the analog of the ‘Kozai constant’). This simply reflects conservation of the $Z$-component of angular momentum since the doubly-averaged potential is axisymmetric. Since $\Theta$ is conserved, we can always infer the time evolution of binary inclination from the behavior of its eccentricity via $\cos^{2}i=(1-e^{2})^{-1}\Theta$. Finally, definitions (7) imply that $e$ and $j$ must obey

[TABLE]

to be physically meaningful for a fixed $\Theta$.

Given these considerations, the secular Hamiltonian is a function of the dimensionless angle-action variables $\omega$ and $j$. The equations of motion fully describing their evolution are

[TABLE]

Our subsequent analysis of binary dynamics is largely based on these evolution equations.

2.1 Phase portrait

Since the dimensionless doubly-averaged Hamiltonian (8) ends up being a function of $j=\sqrt{1-e^{2}}$ and $\omega$, one can get a good perception of the secular dynamics by plotting contours of $H_{1}^{*}$ in the $(\omega,e)$ plane. We do this in Figures 4, 5, 6, & 7 for the $\Gamma$ ranges (3)-(6) respectively, and for varying $\Theta$. In each panel the limiting eccentricity $e_{\mathrm{lim}}=\sqrt{1-\Theta}$ is represented by a dashed black line. The direction in which orbits traverse their trajectories is indicated with green arrows in Figures 4a, 5d, 6a & 7d.

We will explain the features of these phase portraits as we discuss each of them individually in §§3-6. An observation that we would like to make now is that all phase-space trajectories are split into two families: librating and circulating. The librating orbits are closed contours of $H_{1}^{*}$ which loop around a fixed point (always located at $\omega=\pm\pi/2$ as explained in §2.2), whereas circulating orbits run over all $\omega\in(-\pi,\pi)$. The separatrices between families of librating and circulating orbits are indicated with red dashed lines. Depending on the values of $\Gamma$ and $\Theta$, one could have either only circulating orbits (typically the case for large $\Theta$, and always true for $-1/5<\Gamma\leq 0$), only librating orbits (a rare case realized for $\Gamma=1/5$, see Figure 5), or a mix of both (for low enough $\Theta$ in most dynamical regimes).

2.2 Fixed points and orbit families

We start by exploring characteristics of the fixed points of the system, around which phase-space orbits librate. As these points are extrema of $H_{1}^{*}$ in $(\omega,j)$ space, they must be solutions of $\mathrm{d}\omega/\mathrm{d}t=\mathrm{d}j/\mathrm{d}t=0$. With a small amount of algebra, we find from equations (10)-(11) two possible formal solutions for the non-trivial fixed points in our phase-space222Fixed points also exist for all $\omega$ along the lines $j=1$ and $j=\sqrt{\Theta}$, but these are trivial in the sense that they can never be reached by orbits which do not start on those lines. However, they are still important because they bound the phase-space and are often the locations of maximum/minimum $H_{1}^{*}$, see §2.5., namely $(\omega,j)=(\pm\pi/2,j_{\mathrm{f}})$ where

[TABLE]

The value of $\omega$ corresponding to fixed points agrees with the phase portraits discussed in §2.1.

For a given $\Gamma$, fixed points can exist in the $(\omega,j)$ phase-space as long as $j_{\mathrm{f}}(\Theta,\Gamma)$ satisfies the condition (9). Solving the inequality $\sqrt{\Theta}<j_{\mathrm{f}}<1$ gives the following constraint on $\Theta$ for the existence of fixed points:

[TABLE]

Depending on the value of $\Gamma$, this constraint may or may not have meaningful solutions. The functions $\Lambda(\Gamma)$ and $\Lambda^{-1}(\Gamma)$ are plotted in Figure 1. We can see from this plot that there are four distinct $\Gamma$ regimes, given by the ranges (3)-(6). We will return to Figure 1 when discussing the existence of fixed points in §§3-6.

2.3 Does a given orbit librate or circulate?

Next we work out whether a given phase-space orbit with specified $H_{1}^{*}$, $\Theta$, $\Gamma$ librates or circulates. We do this by considering the behaviour at $\omega=0$: from the morphology of the phase portraits (Figures 4, 5, 6, & 7), it is clear that if constant $H_{1}^{*}=H_{1}^{*}(\Theta,\Gamma,\omega,j)$ has a physical solution for $j$ at $\omega=0$ then the orbit circulates; if not, it has to librate about one of the fixed points.

We find from equation (8) a formal solution

[TABLE]

The trajectory is circulating whenever the condition (9) is satisfied for $j(\omega=0)$. If $j(\omega=0)$ does not obey this inequality, then it does not represent a physically meaningful solution. As a result, the orbit must librate around one of the fixed points and never reach $\omega=0$.

Let us define the quantity

[TABLE]

It represents $e^{2}(\omega=0)$ and is an integral of motion since it depends on $H_{1}^{*}$ and $\Theta$. It will prove useful to eliminate $\Theta$ and $H_{1}^{*}$ in the above expression in favour of $e,i,\omega$. We find after some algebra

[TABLE]

Equations (9) and (15) imply that for $0<D<1-\Theta$ the trajectory is circulating, whereas for $D<0$ or $D>1-\Theta$ it librates around a fixed point. If both families exist then the separatrix between them corresponds to either $D=0$ or $D=1-\Theta$ (see the upcoming Figure 3).

2.4 Maximum and minimum eccentricities

We now find the minimum and maximum eccentricities that a binary reaches as it evolves along its phase-space trajectory. We do this by finding the extrema of $j$, which we call $j_{\mathrm{min/max}}$.

From the phase portraits in Figures 4, 5, 6, & 7 it is clear that librating orbits, whenever they exist, have both their minimum and maximum eccentricities at $\omega=\pm\pi/2$. On the other hand, for circulating orbits the maximum and minimum eccentricities can be at either $\omega=0$ or $\omega=\pm\pi/2$ depending on the value of $\Gamma$; see Figures 4 and 5. To find $j_{\mathrm{min/max}}$, we therefore separately plug $\omega=0,\pi/2$ into $H_{1}^{*}$ (equation (8)) and solve for $j$.

The solution for $\omega=0$ is simply the square root of equation (14), and we denote it $j_{0}$:

[TABLE]

For $\omega=\pi/2$ there are two solutions, which we call $j_{\pm}$:

[TABLE]

where

[TABLE]

We would like to stress here that although we use the notation $j_{\pm}^{2},j_{0}^{2}$, these quantities are nothing more than possible solutions to the equations $\mathrm{d}j/\mathrm{d}t=\mathrm{d}\omega/\mathrm{d}t=0$, and should not therefore be interpreted as always positive. Indeed, depending on the $\Gamma$ regime either two or none of $j_{\pm}^{2},j_{0}^{2}$ will lie in the allowed physical range $(\Theta,1)$ (equation (9)) — the remaining one or three will lie outside this range and can even be negative, but are physically irrelevant.

This is demonstrated in Figure 2, which shows $j_{\pm}^{2},j_{0}^{2}$ as functions of $\Gamma$ for various points in $(D,\Theta)$ space. Depending on the dynamical regime (determined by the value of $\Gamma$) each of these $j_{\pm}$, $j_{0}$ solutions, if they exist, can correspond to either the minimum or maximum of $j$, as we describe in further detail in §§3-6.

2.5 Range of parameter values

We would first like to determine the range of values that the integral of motion $H_{1}^{*}$ can take. To do this we need to find the extrema of our Hamiltonian $H_{1}^{*}$ in the $(\omega,e)$ (or equivalently $(\omega,j)$) phase-space. It is clear from the phase portraits (Figures 4, 5, 6 & 7) that extrema of $H_{1}^{*}$ can only occur in three distinct locations: fixed points333Evaluating at $\omega=-\pi/2$ gives the same answers as evaluating at $\omega=\pi/2$. Fixed points exist in both locations, but to avoid confusion with other upcoming $\pm$ signs we prefer just to consider $\omega=\pi/2$ from here on. $(\omega=\pi/2,j=j_{\mathrm{f}})$, the line $j=1$, and the line $j=\sqrt{\Theta}$. Evaluating $H_{1}^{*}$ in these locations, we find

[TABLE]

where444Note that the functions $H_{\pm}$, unlike $j_{\pm}$, are not two independent solutions. The $\pm$ sign is due merely to an algebraic peculiarity that arises when evaluating the single-valued function $H_{1}^{*}$ at $(\pi/2,j_{\mathrm{f}})$, depending on the $\Gamma$ range. For a given $\Gamma$, only one of $H_{\pm}$ is correct. The same consideration holds for $D_{\pm}$ in the next paragraph.

[TABLE]

(No fixed points exist in the range $-1/5<\Gamma\leq 0$). One must then investigate each $\Gamma$ regime independently to work out which of the above corresponds to a maximum or minimum. We will not pursue the details here but we state the results for each $\Gamma$ regime in §§3-6, and summarise them in Table 2.

Rather than $H_{1}^{*},\Theta$, it is often convenient to take $D$ and $\Theta$ to be our primary integrals of motion, so that for fixed $\Gamma$ each phase-space trajectory of a binary corresponds to a point in the $(D,\Theta)$ plane (see the upcoming Figure 3). It is obvious that $\Theta$ can always run between [math] and $1$ for circulating orbits, and is bounded by equation (13) for librating orbits. We would like to know which values $D$ can take for a given $\Theta,\Gamma$. From equation (15) we see that extrema of $D$ are also extrema of $H_{1}^{*}$, so we must evaluate $D$ at the same three locations as $H_{1}^{*}$, see equations (21)-(23). We find

[TABLE]

where

[TABLE]

In each $\Gamma$ regime and for each type of orbit the maximum and minimum $D$ will correspond to some combination of $D_{\pm}$, $D=0$ and $D=1-\Theta$. These limits give rise to distinctive morphologies of the physically allowed regions in $(D,\Theta)$ plane — see Figure 3. They are summarised in Table 1 and discussed in §§3-6.

2.6 Timescales of eccentricity oscillations

We can also derive a general expression for the timescale of secular eccentricity oscillations, for any binary and for any given $\Gamma$.

Since $H_{1}^{*}$ is a conserved quantity we can rearrange (8) to get $\omega$ explicitly in terms of $j$:

[TABLE]

Plugging this into equation (11) allows us to eliminate $\omega$ from $\mathrm{d}j/\mathrm{d}t$, turning it into an equation for $j$ only. Factorizing the result and multiplying both sides by $j$ we arrive at a simple equation for the rate of change of $j^{2}$:

[TABLE]

The square root here is well defined because for $|\Gamma|>1/5$ and $<1/5$ the signs of bracketed terms change in such a way that the whole expression under the square root is positive555For $\Gamma=\pm 1/5$ we have $\mathrm{d}j^{2}/\mathrm{d}t=0$. This reflects the fact that secular evolution stalls and $t_{\rm sec}\to\infty$ for $\Gamma=\pm 1/5$ — see equation (33) and Figure 3., as can be checked using the results collected in Table 1 (also see Figure 2).

The maximum and minimum $j$ reached by a given phase-space orbit and satisfying the constraint (9) are denoted $j_{\mathrm{min/max}}$. They correspond to two of the three possible roots $j_{\pm},j_{0}$ depending on the orbit type and the value of $\Gamma$, as we will see in §§3-6. Regardless of their precise values, an entire oscillation runs from $j_{\mathrm{min}}$ to $j_{\mathrm{max}}$ and back again, so an entire secular oscillation takes

[TABLE]

which is expressible in terms of complete elliptic integrals of the first kind $K(k)=\int_{0}^{\pi/2}\mathrm{d}\alpha/\sqrt{1-k^{2}\sin^{2}\alpha}$ (Gradshteyn & Ryzhik, 2014). Defining666Note that, in general, $\mathrm{min}[j_{+}^{2},j_{-}^{2},j_{0}^{2}]\neq j_{\mathrm{min}}^{2}$. This is because $\mathrm{min}[j_{+}^{2},j_{-}^{2},j_{0}^{2}]$ can take any value (including negative values, see Figure 2), whereas $j_{\mathrm{min}}$ is the physical minimum angular momentum reached by a given binary, which must lie between $\sqrt{\Theta}$ and $1$. An analogous statement holds for $j_{\mathrm{max}}$.

[TABLE]

we find, in general, that

[TABLE]

where

[TABLE]

is the characteristic secular timescale. In the numerical estimate we assumed the binary orbits a spherical cluster with scale radius $b$ and total mass $M$, and defined $A^{*}\equiv A/(GM/b^{3})$. Maps of $A^{*}$ for different cluster potentials and binary orbits are presented in Paper I.

In Paper I we also noted that $A\sim 4\pi^{2}A^{*}T_{\phi}^{-2}$, where $T_{\phi}=2\pi(GM/b^{3})^{-1/2}$ is the characteristic azimuthal period of the outer orbit of the binary around the cluster. Introducing the period of the inner orbit of the binary $T_{b}=2\pi[G(m_{1}+m_{2})/a^{3}]^{-1/2}$, we can then use equation (34) to estimate $t_{1}$ as

[TABLE]

A similar estimate of the characteristic secular timescale — ratio of the square of the outer orbital period to the inner orbital period — is known to hold for the LK problem (Naoz, 2016). Although equation (36) was derived for a spherical cluster potential, we also expect the same scaling to hold in (non-spherical) axisymmetric potentials.

The ratio $t_{\mathrm{sec}}/t_{1}$ is plotted in the $(D,\Theta)$ plane for various fixed $\Gamma$ in Figure 3. In each panel, circulating orbits fill the triangle $0<D<1-\Theta$ while regions of librating orbits are shown with white hashing. In §2.2 we noted that the four distinct $\Gamma$ regimes (equations (3)-(6)) would give rise to different phase-space behaviours. We see in Figure 3 that they also give rise to different morphologies of the allowed regions in the $(D,\Theta)$ plane, to be discussed in §§3-6.

The analytic derivation of eccentricity oscillation period (equation (33)) was previously done in the LK limit ($\Gamma=1$) by Vashkov’yak (1999) and Kinoshita & Nakai (2007), and in the $\Gamma=1/3$ limit by Brasser (2001). Our expression (33) generalizes their results for arbitrary external tidal potentials of the type explored in Paper I.

Much of §§3-6 will be focused on deriving the values of $j_{\mathrm{min/max}}$, $\Delta$, and the bounds on $\Theta$, $D$ and $H_{1}^{*}$ appropriate to each of the distinct $\Gamma$ regimes (3)-(6). For ease of reference all of the results are collected in Tables 1 and 2, which will be discussed in more detail in §9.

3 Secular dynamics in the case $\Gamma>1/5$

In this section we focus exclusively on the case $\Gamma>1/5$. In short, we find the dynamics in this regime to be qualitatively similar to (but quantitatively different from) the ‘test particle quadrupole’ LK problem. A similar investigation in the LK limit ($\Gamma=1$) was carried out by Antognini (2015). This regime also covers the case of $\Gamma=1/3$ relevant for Oort Cloud comets perturbed by the Galactic tide (Heisler & Tremaine, 1986).

In Figure 4 we plot contours of the dimensionless doubly-averaged perturbing Hamiltonian $H_{1}^{*}$ (equation (2)) in the $(\omega,e)$ plane, with $\Gamma=1,\,0.7,\,0.4$ from top to bottom and $\Theta=0.1,\,0.5,\,0.8$ from left to right.

Fixed points exist in the left and centre columns, but not in the rightmost column (large $\Theta$) where there are only circulating orbits. Circulating orbits show prograde pericentre precession $\dot{\omega}>0$, while librating orbits traverse clockwise loops in the $(\omega,e)$ plane. As we increase $\Theta$ (i.e. move from left to right), the maximum eccentricity reached by the average orbit sharply decreases.

We see that whenever a fixed point is present, the circulating orbits run ‘over the top’ of the librating orbits. As a result, the eccentricity of a fixed point provides a lower bound on the maximum eccentricity reached by any orbit in the phase-space777This is not true for general $\Gamma$ — c.f. §4.. Furthermore, the eccentricity of the fixed point increases slightly as we decrease $\Gamma$ (move down the page) — see equation (12). Thus for $\Gamma$ close to (but greater than) $0.2$ and $\Theta\to 0$ more binaries may reach very high eccentricities than for $\Gamma\sim 1$.

We now proceed to explain these features quantitatively.

3.1 Fixed points and orbit families

Looking at Figure 1, we see that whenever $\Gamma>1/5$, for fixed points to exist (shaded regions) $\Theta$ must be less than $\Lambda$. However, $\Lambda$ has a minimum value of $1/2$ in this $\Gamma$ range (namely as $\Gamma\to\infty$); hence fixed points always exist for $\Gamma>1/5$ provided we choose $\Theta$ small enough, see equation (13). The precise requirement is

[TABLE]

The range of $\Theta$ for which fixed points exist increases as $\Gamma$ decreases.

This result allows us to understand the lack of librating orbits in panels (c), (f) and (i) in Figure 4: their $\Theta$ value is too large for the range (37). Physically, at large $\Theta$ the inclination of the binary with respect to the symmetry plane of the cluster is too low for its tidal field to efficiently torque the binary to high eccentricities. In the LK case $\Gamma=1$ we recover the classic result that the critical $\Theta$ range for fixed points to exist is $\Theta\in(0,3/5)$. For initially circular binaries the resulting minimum inclination in then $\cos^{-1}(\sqrt{3/5})\approx 39.2^{\circ}$ (see also §9.1).

We can convert (12) into an eccentricity via $e=\sqrt{1-j^{2}}$:

[TABLE]

This helps us to understand why in Figure 4, the maximum eccentricity is largest for small $\Theta$, and only weakly dependent on $\Gamma$: since each trajectory reaches a maximum eccentricity which is at least $e_{\mathrm{f}}$, decreasing $\Theta$ will increase that maximum. And since $\Lambda$ is a weak function of $\Gamma$ in this range (taking values $\Lambda\in(1/2,1)$ — see Figure 1), dependence on $\Gamma$ is not very strong.

3.2 Range of parameter values

In §2.5 we mentioned that $\Theta$ is bounded by equation (13) for librating orbits and runs between [math] and $1$ for circulating orbits. From equations (25)-(27) we know that in the regime $\Gamma>1/5$, the extrema of $D$ correspond to some combination of $D=0$, $D=1-\Theta$ and $D_{-}$. It can be shown that $D_{-}$ is negative all for $\Theta$ when $\Gamma>1/5$, while $D=1-\Theta$ is obviously positive. Hence the bounds on $D$ are:

[TABLE]

These ranges dictate the morphology of the $(D,\Theta)$ plane in the top row of Figure 3.

It is easy to verify that for $\Gamma>1/5$ the minimum of $H_{1}^{*}$ is attained at $j^{2}=\Theta$ (i.e. along the black dashed lines of limiting eccentricity in Figure 4), so $H^{*}_{1,\mathrm{min}}$ is given by equation (22). As for $H^{*}_{1,\mathrm{max}}$, if fixed points exist for a given $\Gamma,\Theta$, then the Hamiltonian is maximised at the fixed point and $H^{*}_{1,\mathrm{max}}=H_{-}$, see equation (23). If fixed points do not exist, then $H^{*}_{1,\mathrm{max}}$ is attained at $j^{2}=1$ (i.e. along the line of zero eccentricity), and hence is given by equation (21).

3.3 Maximum and minimum eccentricities

For $\Gamma>1/5$, librating orbits, if they exist, will have minimum angular momentum $j_{\mathrm{min}}=j_{-}$ and maximum angular momentum $j_{\mathrm{max}}=j_{+}$, since these are the two solutions at $\omega=\pi/2$, and $j_{+}>j_{-}$ (see Figure 4 and equation (18)).

Meanwhile, circulating orbits also reach minimum angular momentum at $\omega=\pi/2$ (Figure 4), so that their $j_{\mathrm{min}}=j_{-}$. At the same time, $j_{\mathrm{max}}=j_{0}$ for circulating orbits since the maximum value of their angular momentum (lowest eccentricity) is reached at $\omega=0$.

Maximum and minimum eccentricites are then given by $e_{\mathrm{max/min}}=(1-j_{\mathrm{min/max}}^{2})^{1/2}$ for both types of phase-space trajectories.

3.4 Timescales of eccentricity oscillations

To find the timescale $t_{\mathrm{sec}}$ using equation (31) we need to know the values of $j_{\mathrm{min/max}}$ (§3.3) and $\Delta$ for each orbit family. First of all, in the $\Gamma>1/5$ regime it is clear (equation (18)) that we always have $j_{+}^{2}>j_{-}^{2}\geq 0$ (see Figure 2 for an illustration).

Librating orbits in this regime have $D<0$ (see Figure 3), so we see from equation (14) that $j_{0}^{2}>1$. Since $j_{\pm}^{2}$ provide upper and lower bounds on the true angular momentum $j^{2}$, it must be the case that $j_{-}^{2}<j^{2}<j_{+}^{2}<j_{0}^{2}$, from which we read off that librating orbits have $\Delta=j_{0}^{2}-j_{-}^{2}$.

Meanwhile for circulating orbits $D>0$ and we know that $j^{2}$ is bounded from above by $j_{0}^{2}<1$, so we must have $j_{-}^{2}<j^{2}<j_{0}^{2}<j_{+}^{2}$, see Figure 2a. Hence $\Delta=j_{+}^{2}-j_{-}^{2}$ for circulating orbits.

Using these results we plot the ratio $\log_{10}(t_{\mathrm{sec}}/t_{1})$ in $(D,\Theta)$ space for $\Gamma=1,0.5,0.25$ in the top row of Figure 3. The triangle $0<D<1-\Theta$ contains the circulating orbits in each case; the orbits outside of this triangle librate. The bounds on $\Theta$ and $D$ are given by equations (37) and (39) respectively. The timescale for oscillations is seen to depend primarily on the proximity to the separatrix between librating and circulating orbits at $D=0$. Along the separatrix the timescale for secular oscillations diverges888To see this note that as $D\to 0$, we have $j_{+}\to j_{0}=1$, and the function $K(k)$ diverges logarithmically as $k\to 1$.. As we decrease $\Gamma$ from $1$, the region containing librating orbits gets larger, though of course the triangle of circulating orbits is unchanged. The value $t_{\mathrm{sec}}/t_{1}$ is amplified when we decrease $\Gamma$, so that the timescale for oscillations increases as we decrease $\Gamma$ (at fixed $A$).

For any $\Gamma$ in the approximate range $0.25\leq\Gamma\lesssim 0.5$, and sufficiently far from the separatrix, $t_{1}$ provides a good estimate of $t_{\mathrm{sec}}$. For $\Gamma\gtrsim 0.5$, large portions of the $(D,\Theta)$ space have secular timescales that are shorter than $t_{1}$ by a factor of a few. As $\Gamma\to 0.2$, the timescale diverges everywhere in $(D,\Theta)$ space (see equation (33)).

All of the results arrived at in this section will change when we leave the regime $\Gamma>1/5$.

4 The case $0<\Gamma\leq 1/5$

We now turn to the second regime, $0<\Gamma\leq 1/5$, which is realised quite naturally for example by binaries orbiting close to the core of a spherical cluster (see Paper I).

We begin as in §3 by showing the phase portraits as one varies $\Gamma$ and $\Theta$. In Figure 5 we plot contours of $H_{1}^{*}$, with $\Gamma=0.2,0.1,0.01$ from top to bottom and $\Theta=0.1,\,0.5,\,0.8$ from left to right. The green arrows in panel (d) show the sense in which orbits traverse their trajectories. We immediately note qualitative differences between the plots with $0<\Gamma\leq 1/5$ (Figure 5) and those for $\Gamma>1/5$ (Figure 4). For $\Gamma=1/5=0.2$ only librating orbits exist, as we can see in the top row of Figure 5. In panels (d) and (e) we again have fixed points at $\omega=\pm\pi/2$ and librating orbits surrounding them, but note that the circulating orbits which exist for $\Gamma=0.1,0.01$ now have eccentricity minima at $\omega=\pm\pi/2$ and maxima at $\omega=0$, which is the opposite of the $\Gamma>1/5$ case. In the $(\omega,e)$ phase plane, circulating orbits now run ‘underneath’ the librating orbits, whereas for $\Gamma>1/5$ they ran ‘over the top’. As a result, fixed points no longer provide a lower bound on the maximum eccentricity. The librating orbits still run clockwise but the circulating orbits now display retrograde precession, $\dot{\omega}<0$, whereas in the $\Gamma>1/5$ case we had $\dot{\omega}>0$.

4.1 Fixed points and orbit families

Figure 1 shows that for $\Gamma=1/5$ we have $\Lambda=\Lambda^{-1}=1$. This means that all $\Theta\in(0,1)$ accommodate librating trajectories, and nothing circulates. Moving to lower $\Gamma$ in Figure 1 one finds $\Lambda^{-1}(\Gamma)<\Lambda(\Gamma)$, so that for fixed points and librating trajectories to exist $\Theta$ must now now be less than $\Lambda^{-1}$ (i.e. the shaded region is now bounded by the blue curve):

[TABLE]

The range of $\Theta$ for which fixed points exist diminishes as $\Gamma\to 0$.

The fact that circulating orbits have changed their sense of pericentre precession from prograde to retrograde is easily explained by calculating $\mathrm{d}\omega/\mathrm{d}t$ (equation (10)) at $\omega=0$. The result is proportional to $(5\Gamma-1)$, so that $\dot{\omega}$ is positive when $\Gamma>1/5$ and negative when $\Gamma<1/5$.

4.2 Range of parameter values

We again want to derive the bounds on the $(D,\Theta)$ plane and to find the extrema of $H_{1}^{*}$. We begin as in §2.5 by considering the limits on $D$.

First of all, according to its definition (17), $D$ diverges when $\Gamma=1/5$, which is in agreement with the absence of circulating orbits in the top row of Figure 5 (circulating orbits require $0<D<1-\Theta$). For $0<\Gamma<1/5$ the quantity $10\Gamma/(1-5\Gamma)$ is positive; then it follows from equation (16) that the minimum value of $D$ for any fixed $\Theta$ is zero (attained for $e=0$). Hence, $D=0$ is the lower bound. The other possible bounds on $D$ are $D_{-}$ and $1-\Theta$ (see equations (25)-(27)); it turns out that $D_{-}\geq 1-\Theta$ for all $\Theta$ in this $\Gamma$ range, so we conclude that:

[TABLE]

Looking at the timescale plots in the second row of Figure 3, we see that the $(D,\Theta)$ plane morphology has completely changed compared to $\Gamma>1/5$ (top row).

This time the minimum of $H_{1}^{*}$ is situated at $j^{2}=1$ (i.e. along the line of zero eccentricity in Figure 5), so $H_{1,\mathrm{min}}^{*}$ is given by equation (21). If fixed points exist for a given ($\Gamma,\Theta$), then $H_{1,\mathrm{max}}^{*}=H_{-}$ is found at the fixed point, see equation (23). Otherwise $H_{1,\mathrm{max}}^{*}$ is found on the line $j^{2}=\Theta$ (the line of limiting eccentricity), and its value is given by equation (22).

4.3 Maximum and minimum eccentricities

For $\Gamma>1/5$ it was easy to determine for example that librating orbits have $j_{-}^{2}<j^{2}<j_{+}^{2}<j_{0}^{2}$. From this we were able to instantly read off $j_{\mathrm{min/max}}=j_{\pm}$ and $\Delta=j_{0}^{2}-j_{-}^{2}$ (§3.4). When $\Gamma\leq 1/5$ things become more complicated. While it is possible to calculate the analagous results in each $\Gamma$ regime algebraically, it is easier and more informative to look at Figure 2, in which we plot $j_{\pm}^{2},j_{0}^{2}$ as a function of $\Gamma$ for various fixed $(D,\Theta)$.

Let us consider only Figure 2a to begin with, which is for $\Theta=0.15,D=0.15$ and therefore sits inside the triangle of circulating orbits ($0<D<1-\Theta$) for any $\Gamma\neq 1/5$. The properties of circulating orbits in each $\Gamma$ regime can be read off from Figure 2a. In the current regime $0<\Gamma<1/5$ we see that $j_{-}^{2}<j_{0}^{2}<j_{+}^{2}$. Moreover, the horizontal dotted lines represent $j^{2}=\Theta$ and $j^{2}=1$, so physical solutions $j^{2}$ must lie between these two lines. Hence $j^{2}$ must be bounded by the green and blue lines (that is by $j_{0}^{2}$ and $j_{+}^{2}$). We therefore deduce that $j_{\mathrm{min}}=j_{0}$, $j_{\mathrm{max}}=j_{+}$, and $\Delta=j_{+}^{2}-j_{-}^{2}$.

This exercise can be repeated with panels (b) and (c) to find the corresponding results for librating orbits. The only subtlety is that as one changes $\Gamma$, the region of the $(D,\Theta)$ plane corresponding to librating orbits changes shape and moves (this is most clearly seen by comparing different panels along the same row in Figure 3). The $(D,\Theta)$ points we have chosen for panels (b) and (c) in Figure 2 do not quite fall inside the region of librating orbits for some $\Gamma$ because that region becomes too small, which is why for example the red and blue curves in panel (b) are not defined in a small range near $\Gamma=-0.2$. However these plots are good enough for anticipating the results of the algebraic calculations. When $0<\Gamma<1/5$ the $j^{2}$ ranges for librating orbits can be read off from Figure 2c: for $j^{2}$ to lie between the upper and lower horizontal dotted lines we must have $j_{0}^{2}<j_{-}^{2}<j^{2}<j_{+}^{2}$; as a result, $j_{\mathrm{min}}=j_{-}$, $j_{\mathrm{max}}=j_{+}$, and $\Delta=j_{+}^{2}-j_{0}^{2}$.

4.4 Timescales of eccentricity oscillations

In the second row of Figure 3 we present contour plots of $\log_{10}(t_{\mathrm{sec}}/t_{1})$ in $(D,\Theta)$ space for $\Gamma=0.15,0.1,0.01$. This time the bounds on $\Theta$ and $D$ are given by equations (40) and (41) respectively. The triangle $0<D<1-\Theta$ again contains the circulating orbits, but now the librating orbits have moved to the right of this triangle, and the separatrix corresponds to $D=1-\Theta$. Along the separatrix the timescale for secular oscillations again diverges. The timescale is also infinite everywhere in the $(D,\Theta)$ plane in the special case $\Gamma=1/5$ (see equation (33)).

For $0\leq\Gamma\lesssim 0.15$, we see that $t_{1}$ again provides a fair estimate of $t_{\mathrm{sec}}$, although it is really a lower bound on $t_{\mathrm{sec}}$ in much of the $(D,\Theta)$ space, whereas for $\Gamma>1/5$ it was typically an effective upper bound.

5 The case $-1/5<\Gamma\leq 0$

The two regimes $-1/5<\Gamma\leq 0$ and $\Gamma\leq-1/5$ cannot be realised by binaries orbiting spherical potentials. In fact, they typically require rather extreme orbits in strongly aspherical potentials. For this reason, here and in §6 we simply summarise the qualitative results for each regime — the details are given in Appendices A and B respectively.

The regime $-1/5<\Gamma\leq 0$ is typically realised when the binary’s outer orbit makes large excursions in the $Z$ direction in an oblate potential, i.e. is highly inclined with respect to the potential’s symmetry plane (see Paper I). Orbits in prolate potentials can also result in this range of $\Gamma$.

In Figure 6 we plot contours of constant $H_{1}^{*}$ for $\Gamma=-0.01,-0.1,-0.18$ from top to bottom and $\Theta=0.1,\,0.5,\,0.8$ from left to right. The phase portrait has undergone another bifurcation as we passed through $\Gamma=0$. We see from Figure 6 that only circulating orbits exist in this $\Gamma$ regime, all with retrograde precession ($\dot{\omega}<0$). Fixed points and librating trajectories do not emerge at all in this regime, which is explained in Appendix A along with some more technical details.

6 The case $\Gamma\leq-1/5$

Situations where $\Gamma\leq-1/5$ may arise, for example, from orbits that are highly inclined with respect to the symmetry plane of a strongly flattened potential. Cylindrical potentials, i.e. ones in which $\Phi(R,Z)=\Phi(R)$ (no $Z$-dependence, extremely prolate configurations), also fall into this regime as they always have $\Gamma=-1/3$, see Paper I.

In Figure 7 we plot contours of $H_{1}^{*}$ with $\Gamma=-0.2,-0.4,1$ from top to bottom and $\Theta=0.1,\,0.5,\,0.8$ from left to right. A final bifurcation has occured as we moved below $\Gamma=-1/5$, such that the phase portrait morphology now looks similar to the case $\Gamma>1/5$, with fixed points at $\omega=\pm\pi/2$ and circulating orbits running ‘over the top’ of librating orbits in the $(\omega,e)$ plane. However, the direction of $\omega$ precession is the opposite to the $\Gamma>1/5$ case: circulating orbits in the $\Gamma\leq-1/5$ regime have retrograde precession ($\dot{\omega}<0$) while librating orbits run anticlockwise. Appendix B provides additional details on this $\Gamma$ regime.

This completes our detailed exploration of the dynamical regimes corresponding to the different ranges of $\Gamma$ (equations (3)-(6)) considered in this work.

7 Accuracy of the doubly-averaged approximation

In Paper I we developed three successive levels of approximation for the evolution of the binary’s inner orbital elements.

First, we had the set of six time-dependent equations for the relative position and velocity of the binary components, subject to the gravitational force of each other and the smooth background potential of the cluster. No tidal approximation had been made at this stage. These six equations can only be solved by direct numerical (‘N-body’) integration.

Second, we had a set of four ‘singly-averaged’ (SA) equations, which were obtained by tidally expanding the six N-body equations, recasting them in Hamiltonian form, and averaging over the binary’s mean anomaly (i.e. over the ‘inner orbit’). The singly-averaged equations are still explicitly time-dependent through the external potential $\Phi(\mathbf{R}_{\mathrm{g}}(t))$, where $\mathbf{R}_{\mathrm{g}}(t)$ is the barycentric position of the binary (assumed to move as a test particle in the cluster potential).

Finally, we derived a system of two ‘doubly-averaged’ (DA) time-independent equations (10)-(11) resulting from the secular Hamiltonian (1), which was itself obtained by averaging the singly-averaged Hamiltonian over many ‘outer orbits’ of the binary around the cluster. Our time-averaging procedure relied on the assumption that $\mathbf{R}_{\mathrm{g}}(t)$ densely fills an axisymmetric torus whose symmetry axis coincides with the symmetry ($Z$) axis of the potential (in the case of a spherical potential, the ‘torus’ is a two-dimensional annulus perpendicular to $Z$). More technically, we required that the time-averages of the derivatives $\Phi_{\alpha\beta}$ of the potential (where $\Phi_{xy}=\partial^{2}\Phi/\partial X\partial Y$, etc.) converge to constant values — in particular, $\overline{\Phi}_{xy}\equiv\overline{\Phi}_{xz}=\overline{\Phi}_{yz}=0$. This condition is almost always satisfied for orbits in any axisymmetric potential after a sufficient number of outer orbital periods. Then, the time-dependent torque on the binary can be replaced with a converged time-independent torque that arises from an axisymmetric mass distribution as seen from the binary’s frame. In §§2-6 of the present paper we have examined the dynamics that arise from the doubly-averaged equations (10)-(11).

However, it is to be expected that the DA theory will break down under certain circumstances. The goal of this section is to explore the validity of the DA approximation for computing inner binary dynamics in the different $\Gamma$ regimes covered in §§3-6. To do so, we present several examples comparing numerical integrations of the N-body999Note the term ‘N-body’ is not meant to imply that we integrate an entire cluster of, say, $10^{5}$ particles. Instead we integrate the exact two-body equations of motion in the presence of the time-dependent external potential calculated via numerical orbit integration of $\mathbf{R}_{\mathrm{g}}(t)$., SA and DA equations.

We will see that the secular approximation is good as long as the ratio of the secular timescale to the outer orbital period, $t_{\mathrm{sec}}/T_{\phi}$, is large enough. Given that $t_{\mathrm{sec}}\sim T_{\phi}^{2}/T_{\mathrm{b}}$ (see equations (33) and (36)), this is equivalent to the statement that $T_{\phi}/T_{\mathrm{b}}$ be sufficiently large. Hence, alongside each case where the DA theory fails we present another example with identical initial conditions except for a smaller binary semi-major axis. The effect of this is to decrease the inner binary orbital period $T_{\mathrm{b}}$, rendering the DA theory valid.

7.1 Method

We use two particular forms of the background potential $\Phi$ in our examples. The first is the spherical isochrone potential

[TABLE]

where $r=\sqrt{X^{2}+Y^{2}+Z^{2}}$ is the spherical radius. This is a model potential for a spherical star cluster with total mass $M$ and scale radius $b$ (Binney & Tremaine, 2008). Since the potential (42) is spherical we can always choose the plane of the binary’s outer orbit $\mathbf{R}_{\mathrm{g}}$ to be the $(X,Y)$ plane. The other potential we will use is the Miyamoto-Nagai potential (Miyamoto & Nagai, 1975):

[TABLE]

where $R=\sqrt{X^{2}+Y^{2}}$ is the usual cylindrical radius. Here $M$ is the total mass of the model, $b_{\ell}$ is the scale length and $b_{h}$ is the scale height. By changing the value of $b_{h}/b_{\ell}$, the Miyamoto-Nagai potential interpolates between the Kuzmin potential of a razor thin disk ($b_{h}\ll b_{\ell}$) and the spherical Plummer potential ($b_{h}\gg b_{\ell}$) frequently used to model globular clusters (Binney & Tremaine, 2008).

In either case the binary’s outer orbit $\mathbf{R}_{\mathrm{g}}(t)$ is stipulated via its initial conditions $(R,v_{\mathrm{R}},Z,v_{\mathrm{Z}},\phi,v_{\mathrm{\phi}})$, where $\phi=\tan^{-1}({Y/X})$ is the azimuthal angle in cylindrical coordinates, $v_{\mathrm{R}}$ is the velocity of the binary in the direction of increasing $R$, etc. We integrate the orbit $\mathbf{R}_{\mathrm{g}}(t)$ in this potential numerically using galpy (Bovy 2015; see Appendix F of Paper I for details), and then feed the resulting time series into the SA and N-body equations101010Note that in some examples where there is an extremely large separation between the secular timescale and inner orbital period it becomes prohibitively expensive to integrate the N-body equations of motion, so we just show the DA and SA results.. The DA equations (10)-(11) are integrated using the theoretical values of $A$ and $\Gamma$ when they are available (i.e. in the spherical isochrone case); otherwise we use the numerical prescription outlined in Appendix F of Paper I and denote them by111111In spherical potentials, $A$ and $\Gamma$ can be calculated theoretically by stipulating the outer orbit’s peri/apocentre $r_{\mathrm{p/a}}$. See Paper I for more information about calculating $A$, $A_{\mathrm{num}},\Gamma$, $\Gamma_{\mathrm{num}}$ and possible small discrepancies between the theoretical and numerical values. $A_{\mathrm{num}},\Gamma_{\mathrm{num}}$. In all numerical examples the binary has constituent masses $m_{1}=m_{2}=0.5M_{\odot}$.

7.2 Accuracy of the doubly-averaged approximation for $\Gamma>1/5$.

Here we give two examples where the N-body, SA and DA integrations are in very good agreement, both in the regime $\Gamma>1/5$ (explored in §3). In the first the binary orbits the spherical isochrone potential, and in the second it orbits the Miyamoto-Nagai potential. The details of each example are given in the following paragraphs, and the results are shown in Figures 8 and 9 respectively.

[Figure 8: Isochrone potential, $\Gamma=0.4$.]

We consider a binary orbiting a spherical isochrone cluster (42) with scale radius $b=1\mathrm{pc}$ and total mass $M=10^{5}M_{\odot}$. The initial conditions for the outer orbit $\mathbf{R}_{\mathrm{g}}$ are as follows:

[TABLE]

It is easy to show that this choice of initial conditions is equivalent to a choice of peri/apocentre $(r_{\mathrm{p}}/b,r_{\mathrm{a}}/b)=(2.40,3.03)$. The theoretical $A,\Gamma$ values that result are $A=0.0206(GM/b^{3})$ and $\Gamma=0.4$. The outer orbit is shown in Figure 8a. In Figure 8 all panels show $100T_{\phi}$ worth of data.

The (rather soft) binary has initial orbital elements $(a,e,i,\Omega,M)=(10^{3}\mathrm{AU},0.5,70^{\circ},17.188^{\circ},161.36^{\circ})$ with three different initial values of $\omega$, namely (A) $91.67^{\circ}$ (librating phase-space orbit), (B) $5.73^{\circ}$ (circulating orbit), and (C) $34.14^{\circ}$ (an orbit very close to the separatrix). From Figures 8b,c we see that the agreement is extremely good between N-body, SA and DA calculations for phase-space trajectories (A) and (B), even though their secular timescales are longer than $T_{\phi}$ by just several tens. Although initially the agreement in trajectory (C) is also excellent, there is a divergence between N-body, SA and DA calculations when we reach low eccentricity because (C) was chosen to be so close to the separatrix (where the secular timescale formally is infinite, see Figure 3). The DA result circulates while the others librate.

[Figure 9: Miyamoto-Nagai potential, $\Gamma_{\mathrm{num}}=0.370$.]

This time the binary orbits the Miyamoto-Nagai potential (43) with $b_{\ell}=b_{h}=3.5\mathrm{kpc}$ and total mass $M=10^{11}M_{\odot}$. The initial conditions of the outer orbit are

[TABLE]

The projections of the outer orbit onto the $(X,Y)$ and $(R,Z)$ planes shown in Figures 9a,b. Again all panels show the first $100T_{\phi}$ of integration time. Since the outer orbit is not planar we do not have theoretical $A,\Gamma$ values; numerically we find $A_{\mathrm{num}}=0.0149(GM/b_{\ell}^{3})$ and $\Gamma_{\mathrm{num}}=0.370$ (not too different from $\Gamma=1/3$ corresponding to binaries near the midplane of a thin disc, see Paper I).

The binary has initial orbital elements $(a,e,i,\Omega,M)=(5\times 10^{4}\,\mathrm{AU},0.5,70^{\circ},17.188^{\circ},161.36^{\circ})$, and we consider three initial values of $\omega$, namely (A) $91.67^{\circ}$ (librating orbit), (B) $5.73^{\circ}$ (circulating orbit), and (C) $30.65^{\circ}$ (an orbit very close to the separatrix). It is again evident (Figures 9c,d) that the agreement is extremely good between N-body, SA and DA calculations for trajectories (A) and (B). The agreement in the eccentricity timeseries for the separatrix trajectory (C) is also good over the first half-oscillation, but there is then once again a discrepancy in the phase portrait as to whether the orbit ought to librate or circulate. Since the semi-major axis used in this example is typical for Oort Cloud comets, we conclude that the DA approximation should work well for characterising the secular evolution of the long-period comets (Heisler & Tremaine, 1986).

In both of these examples the DA theory is very accurate, except for describing phase-space trajectories that lie extremely close to the separatrix between librating and circulating orbits. Note that in each case, $100$ outer orbital periods was enough time for at least two secular cycles to take place, so $t_{\mathrm{sec}}/T_{\phi}$ was at most $\sim 50$ and usually smaller. We will see in the upcoming sections that in some circumstances a much greater timescale separation $t_{\mathrm{sec}}/T_{\phi}$ is required for the secular approximation to be valid.

7.3 Accuracy of the doubly-averaged approximation for $0<\Gamma\leq 1/5$.

We now consider the case $0<\Gamma\leq 1/5$ studied in §4. We again provide one example in the spherical isochrone potential and one in the non-spherical Miyamoto-Nagai potential.

This time, in each instance we use two different initial semi-major axes to show how the secular approximation improves for smaller $a$ (shorter $T_{\mathrm{b}}$).

[Figure 10: Isochrone potential, $\Gamma=0.162$.]

The potential is exactly as in the $\Gamma=0.4$ case discussed in §7.2 (Figure 8), but now we choose different initial conditions for the outer orbit, namely

[TABLE]

This corresponds to $(r_{\mathrm{p}}/b,r_{\mathrm{a}}/b)=(0.9,1.6)$, giving the theoretical values $\Gamma=0.162$ and $A=0.124(GM/b^{3})$. Figure 10a shows the outer orbit for the first $100T_{\phi}$ of integration time.

In Figures 10b,c the binary’s initial orbital elements are also unchanged from those in Figure 8 (in particular, we still use $a=10^{3}\mathrm{AU}$). We again integrate for $100T_{\phi}$. We see that the agreement between N-body, SA and DA calculations is good for trajectory (A) and reasonable for (B). Trajectory (C), which is very near the separatrix, does not show very good agreement between the DA and other approximations. The reason that the DA theory is so much less accurate in this example than for $\Gamma=0.4$ (Figure 8) — despite having similar ratios of $t_{\mathrm{sec}}/T_{\phi}$ — is that this time the binary does not fill its annulus fast enough for the time-averaged potential coefficients $\overline{\Phi}_{\alpha\beta}$ to converge rapidly.

To show that the secular approximation can be improved, we rerun integrations of the same three trajectories with exactly the same initial conditions except we use a smaller semi-major axis, $a=100\mathrm{AU}$ (lowering $T_{b}$ by $\sim 30$ and, correspondingly, increasing $t_{\rm sec}$ by the same factor). We integrate for $3000T_{\phi}$. The results are shown in Figures 10d,e. The secular timescales are much longer now (from $t_{\mathrm{sec}}/T_{\phi}\sim 650$ for trajectory (A), to $\sim 1300$ for (C)). As a result, the binary fills its annulus many times per secular period and the agreement between DA and SA integrations is almost perfect.

[Figure 11: Miyamoto-Nagai potential, $\Gamma_{\mathrm{num}}=0.045$.]

The binary orbits the Miyamoto-Nagai potential (43) with $b_{\ell}=3.5\mathrm{kpc}$ and $b_{h}/b_{\ell}=5$, which is a less flattened potential than the $b_{h}/b_{\ell}=1$ example from Figure 9. The total mass is again $M=10^{11}M_{\odot}$. The initial conditions of the outer orbit are

[TABLE]

From Figures 11a,b (which both show the first $100T_{\phi}$ of integration time) we see that $\mathbf{R}_{\mathrm{g}}$ makes large excursions in the $Z$ direction: the binary’s barycentric orbit is about as thick vertically as it is radially. Numerically we find $A_{\mathrm{num}}=0.00826(GM/b_{\ell}^{3})$ and $\Gamma_{\mathrm{num}}=0.045$.

In Figures 11c,d the binary has initial orbital elements $(a,e,i,\Omega,M)=(5\times 10^{4}\mathrm{AU},0.9,70^{\circ},17.188^{\circ},161.36^{\circ})$ and three initial values of $\omega$, namely (A) $90^{\circ}$ (librating orbit), (B) $24.9^{\circ}$ (circulating orbit), and (C) $40.1^{\circ}$ (separatrix). We integrate the equations of motion for $100T_{\phi}$. Although the N-body and SA integrations agree extremely well (the red and green lines in Figures 11c,d are almost indistinguishable), the agreement with the DA integration is only reasonable for trajectories (A) and (B) and is poor for trajectory (C). Moreover, the smooth periodicity of the DA solution has disappeared in the SA and N-body integrations, which show chaotic small-scale oscillations. This is unsurprising — despite the fact that $t_{\mathrm{sec}}/T_{\phi}\gtrsim 50$ in each case, the binary does not fill its torus densely enough over sufficiently few $T_{\phi}$ to render the axisymmetric secular approximation valid.

The DA theory fares much better in Figures 11e,f, in which we rerun the same three trajectories except with a smaller semi-major axis, $a=10^{4}\mathrm{AU}$ (i.e. lowering $T_{b}$ by $\sim 10$ and increasing $t_{\rm sec}$ by the same amount). We integrate for $1000T_{\phi}$, so that $t_{\mathrm{sec}}/T_{\phi}\gtrsim 500$. Of course, the timescales involved here are much longer than the age of the universe and so are not relevant in practice, but this example demonstrates how the secular approximation becomes more accurate when the binary’s outer orbit has a better chance to fill its axisymmetric torus.

7.4 Accuracy of the doubly-averaged approximation in the cases $-1/5<\Gamma\leq 0$ and $\Gamma\leq-1/5$

In this section we present one numerical example in each of the regimes $-1/5<\Gamma\leq 0$ and $\Gamma\leq-1/5$ (explored in §5 and §6 respectively). In both cases we use a Miyamoto-Nagai potential with $b_{h}/b_{\ell}=1$; each time we give an example in which the DA theory works very poorly, and one in which it works well.

[Figure 12: Miyamoto-Nagai potential, $\Gamma_{\mathrm{num}}=-0.163$.]

The potential and outer orbit initial conditions are exactly as in Figure 9, except that the initial vertical coordinate is $Z=2b_{\ell}$. The resulting (vertically extended) outer orbit (shown in Figures 12a,b for the first $100T_{\phi}$ of integration time) results in a negative $\Gamma$ value: we find $\Gamma_{\mathrm{num}}=-0.163$ and $A_{\mathrm{num}}=0.0392(GM/b_{\ell}^{3})$.

For clarity we only show a single phase-space orbit (‘B’) in Figure 12c (integrated for $100T_{\phi}$), which obviously circulates since there are no fixed points in this $\Gamma$ regime (§5). The initial orbital elements for trajectory (B) are identical to those in Figure 9 — in particular, we again use $a=5\times 10^{4}\mathrm{AU}$. From the $(\omega,e)$ phase-space portrait (Figure 12c) we see that the N-body and SA integrations agree nicely, but the agreement with the DA theory is very poor. Comparing with Figure 9, we note that we have caused the DA theory to fail simply by changing a single parameter, the initial $Z$ coordinate of the outer orbit. This is true because the outer orbit is now much more vertically extended and so takes many more outer orbital periods to fill its torus, delaying the convergence of the time averages of the potential derivatives $\Phi_{\alpha\beta}$ that enter the DA Hamiltonian.

Much better agreement is found in Figure 12d, in which we use a semi-major axis of $a=5\times 10^{3}\mathrm{AU}$ (increasing the secular timescale by a factor $\sim 30$) and integrate for $2000T_{\phi}$. Other than semi-major axis, trajectory (B) in this figure has the same initial conditions as trajectory (B) in Figure 12c. Trajectory (A) differs from (B) only in that it has initial $\omega=91.67^{\circ}.$

[Figure 13: Miyamoto-Nagai potential, $\Gamma_{\mathrm{num}}=-0.384$.]

The potential and the initial conditions of the outer orbit are exactly the same as in Figure 12, except the initial $Z$ coordinate is now even larger, $Z=3b_{\ell}$, making the orbit very highly inclined. This thickens the outer orbit’s torus further (the first $100T_{\phi}$ are shown in Figures 13a,b) and results in $\Gamma_{\mathrm{num}}=-0.384$ and $A_{\mathrm{num}}=0.0142(GM/b_{\ell}^{3})$. From panel (a) it is clear that the outer orbit has not come close to filling its torus after $100T_{\phi}$.

Figure 13c shows a trajectory with exactly the same initial conditions as trajectory (B) in Figure 12c (i.e. $a=5\times 10^{4}$ AU), integrated for $100T_{\phi}$. The DA theory is seen to be very inaccurate here. For a better example, in Figure 13d we again lower the semi-major axis to $a=5\times 10^{3}\mathrm{AU}$ and integrate for $2000T_{\phi}$. We see that this time (B) is captured almost perfectly by the DA theory, as are two other trajectories, namely (A) (identical initial conditions except $\omega=91.67^{\circ}$) and (C) (with $\omega=67.6^{\circ}$).

7.5 Discussion: validity of the doubly-averaged secular theory

The examples presented in §7.2-7.4 illustrate two possible ways in which the DA theory can be in error:

(i) The secular approximation is only a good one if the timescale for evolution of the inner orbital elements is much longer than the time for the time-averages of $\Phi_{\alpha\beta}$ to converge (i.e. for the binary to ‘fill its torus’). Otherwise, the secular approximation can break down and the DA equations can fail completely to describe the evolution.

(ii) The torque experienced by the binary fluctuates on the timescale of its outer orbital period, leading to small fluctuations in the orbital elements on this timescale. Even when the DA equations provide a good description of the dynamics on average, they will always fail to resolve these short-timescale fluctuations (which are fully captured at the level of the SA approximation).

So far in this section we primarily explored the validity of the secular approximation in the sense of $\Phi_{\alpha\beta}$ convergence (effect (i)). However, in practice both effects (i) and (ii) are typically present in our calculations and distinguishing them is important.

The best illustration of the two effects operating simultaneously is provided by Figures 12c and 13c, which show that the secular (DA) approximation completely fails to follow the dynamics quantitatively over long timescales. First, the overall shape of each trajectory (the ‘carrier signal’) computed with both the N-body and SA theory does not line up with the DA prediction. This mismatch is entirely the consequence of effect (i). Additionally, the N-body and SA curves also exhibit short-timescale ‘wiggles’ on top of the smoother ‘carrier’ phase curve. These wiggles are caused by effect (ii), i.e. torque variations along the outer orbit.

Regardless of how good is the secular approximation by the measure of $\Phi_{\alpha\beta}$ convergence, there are always small short-timescale fluctuations due to effect (ii) that the DA theory cannot capture (e.g. see Figures 9c,d, 10b,c, & 11c,d). In many applications of the (quadrupole) LK theory effect (i) is largely absent because there the timescale for a binary to effectively ‘fill its torus’ is simply equal to its outer orbital period. The primary deviations from the DA prediction in the LK case are then due to short-timescale fluctuations (effect (ii)), which has been studied widely in this setup (Ivanov et al. 2005; Katz & Dong 2012; Antonini & Perets 2012; Bode & Wegg 2014; Antonini et al. 2014; Antognini et al. 2014; Luo et al. 2016; Grishin et al. 2018). While we largely pass over them from now on, the short-period SA fluctuations can become vitally important when eccentricities get very close to $1$ (see LK references above). Their magnitude depends strongly on the shape of the outer orbit. Accounting for such fluctuations in a systematic way in the general $\Gamma$ case is a topic for future work.

Focusing now on the issue of the $\Phi_{\alpha\beta}$ convergence, the two examples presented in §7.2 (Figures 8 and 9 — both in the regime $\Gamma>1/5$) showed very good agreement between DA theory and direct numerical integration, even for binaries whose secular evolution timescales were significantly shorter than $100T_{\phi}$ (for example, trajectory (A) in Figure 8 had $t_{\mathrm{sec}}/T_{\phi}\sim 25$). However, all other examples required a much larger ratio of $t_{\mathrm{sec}}/T_{\phi}$ for the DA theory to be rendered accurate (typically $\sim$ a few $\times 100$). This is because the secular approximation is valid only when the timescale for secular evolution is much longer than the time for the binary to fill its torus densely. The number of outer orbital periods required to fill the torus densely depends strongly on the form of the potential and the choice of outer orbit.

The DA approximation is often most easily satisfied (i.e. it works for relatively small values of $t_{\mathrm{sec}}/T_{\phi}$) in spherical potentials, because then the ‘torus’ reduces to a two-dimensional annulus (e.g. Figure 8). Not only does this decrease the volume that must be filled, but also the derivatives $\Phi_{xz},\Phi_{yz}$ automatically vanish and so pose no problem to the convergence. Circular outer orbits and orbits that avoid any central core tend to fill their annuli particularly efficiently (and typically correspond to $\Gamma>1/5$; see Paper I).

However, many spherical cluster potentials are cored (such as the isochrone and Plummer models), and so binaries that spend significant time near the centre of these clusters (i.e. those with small $r_{\mathrm{p}}$) experience an almost-harmonic potential. Since orbits in a harmonic potential are closed ellipses, the apsidal precession of such outer orbits can be very slow; this frequently leads to unfilled gaps being left in the annulus even after $\sim 100T_{\phi}$ (see Paper I for an example). As a result, the secular approximation may require relatively large values of $t_{\mathrm{sec}}/T_{\phi}$ to be valid, as was the case in Figure 10. In spherical potentials, small $r_{\mathrm{p}}$ tends to correspond to small (but always positive) $\Gamma$, so this issue will arise most often for binaries in the regime $0<\Gamma\leq 1/5$.

In non-spherical potentials the situation is often worse simply because there is a third dimension of the torus for the outer orbit to fill. In addition, the derivatives $\Phi_{xz},\Phi_{yz}$ are no longer identically zero in general, so we must wait for them to converge, and this typically takes longer than for the other $\Phi_{\alpha\beta}$ (see §7.2 of Paper I). In these potentials the secular approximation is most easily satisfied by binaries on outer orbits that are coplanar or nearly coplanar — then the torus is small in volume, and, if the potential is strongly flattened, the vertical oscillations tend to be very rapid, so the torus is filled efficiently (a good example is Figure 9). Orbits confined near the midplane of a strongly flattened potential have $\Gamma\approx 1/3$.

On the other hand, when the outer orbit has a large vertical extent, filling a torus takes many more outer orbital periods and hence very large values of $t_{\mathrm{sec}}/T_{\phi}$ are required for the secular theory to be accurate (e.g. Figures 12 & 13). This in turn implies that for $\Gamma<0$, the DA theory may be of limited practical use in certain cases (when $t_{\rm sec}/T_{\phi}$ is not large enough) because achieving negative $\Gamma$ typically requires outer orbits that make large excursions in the $Z$ direction.

8 Effect of general relativistic precession on the cluster-tide driven evolution

So far our secular theory considered only the gravitational tidal effect of a stellar cluster on binary evolution. However in a realistic astrophysical situation there could be other, short-range forces which must be taken into account, particularly at high eccentricity when the pericentre distance becomes small (see §9.2). Depending on the type of binary (i.e the masses and sizes of its components) and its semi-major axis these could include (i) prograde precession of the argument of pericentre $\omega$ due to general relativity (GR), (ii) precession due to the oblateness or tidal distortions of the binary components, (iii) loss of energy and angular momentum due to gravitational wave emission, (iv) tidal dissipation within the components of the binary, etc. The first two effects are conservative in that they do not change the energy of the system and preserve binary semi-major axis, while the latter two lead to energy losses and tend to shrink the binary orbit.

In this section we will only consider (i), namely GR pericentre precession. We can include GR precession in our doubly-averaged theory by adding the following term to the Hamiltonian (1) (Fabrycky & Tremaine, 2007):

[TABLE]

The angle brackets remind us that this term is derived by averaging over the binary’s mean anomaly. The Hamiltonian $\langle H_{\mathrm{GR}}\rangle_{M}$ is independent of the longitude of the ascending node $\Omega$, so the $z$-component of angular momentum $J_{z}$ is conserved; hence, $\Theta=(1-e^{2})\cos^{2}i$ remains an integral of motion. Another integral of motion is the full perturbation energy $\overline{\langle H_{1}\rangle}_{M}+\langle H_{\mathrm{GR}}\rangle_{M}$ — the sum of the cluster tide and GR Hamiltonians, equations (1) and (48).

We put (48) into dimensionless form by dividing by $C=Aa^{2}/8$ (see equation (8)). Then we must compare the perturbation $H_{1}^{*}$ due to the potential of the cluster to the corresponding dimensionless perturbation due to GR:

[TABLE]

where the relative strength of GR precession compared to the cluster tide is measured by the (not necessarily small) parameter

[TABLE]

In the numerical estimate we have again assumed that the binary is orbiting a spherical cluster with scale radius $b$ and total mass $M$, and typical values of $A^{*}$ are given in Paper I. The parameter $\epsilon_{\mathrm{GR}}$ represents, up to a constant factor, the ratio of the GR apsidal precession rate for a circular binary to the binary precession rate due to the cluster tide.

The prograde pericentre precession rate induced by GR is

[TABLE]

With this we can evaluate the ratio of $\dot{\omega}_{\mathrm{GR}}$ to the precession rate $\dot{\omega}_{1}$ due to the background cluster tide alone (equation (10)):

[TABLE]

We expect GR effects to be most important at high eccentricity, when $j^{2}\sim\Theta\ll 1$. Also, in most cases of interest $j_{\mathrm{min}}$ occurs at $\omega=\pi/2$. Plugging these relations into (53) we can evaluate

[TABLE]

We can ignore GR precession only when $|\dot{\omega}_{\mathrm{GR}}/\dot{\omega}_{1}|\lesssim 1$, which requires rather small $\epsilon_{\mathrm{GR}}\lesssim 60\Gamma\Theta^{1/2}(j/\Theta^{1/2})^{-1}$.

Since $H^{*}_{\mathrm{GR}}$ is independent of $\omega$, the equation for $\mathrm{d}j/\mathrm{d}t$ is unchanged from (11). As a result the maximum eccentricity reached by a binary perturbed by the Hamiltonian $H^{*}\equiv H_{1}^{*}+H_{\mathrm{GR}}^{*}$ will be found at the same value of $\omega$ as for a binary perturbed only by $H_{1}^{*}$. Following the derivation of equation (30), we can express $\omega$ through $j$ and $H^{*}$ using equations (8) and (49) and plug the result into equation (11). After some manipulations we find that

[TABLE]

where the definitions of $j_{\pm},j_{0}$ now involve $H^{*}$ replacing $H_{1}^{*}$ everywhere. In the limit $\epsilon_{\mathrm{GR}}\to 0$ this equation reduces to (30). We have used equation (55) to find maximum eccentricities when calculating merger rates of compact object binaries in globular and nuclear star clusters (Hamilton & Rafikov, 2019b).

In principle, the minima and maxima of $j$ (and hence $e$) in the presence of GR precession can be calculated analytically. Indeed, the 8th order polynomial inside the radical on the right-hand side of equation (55) can be factorized into the product of $j$, a depressed cubic polynomial, and a quartic polynomial; the latter two have analytic roots, which correspond to the extrema of $j$ (since $\mathrm{d}j^{2}/\mathrm{d}t=0$ when $j$ equals one of these roots). However, the symbolic expressions for these roots are too complicated to be worth presenting here.

One can draw qualitative conclusions regarding the role of GR precession in our problem. Typically its effect is to increase the rate at which $\omega$ precesses, i.e. the rate at which the binary’s apsidal line changes its orientation with respect to the fixed cluster axes. Then the cluster has less time to coherently torque the binary per secular cycle, resulting in a reduced maximum eccentricity (see e.g. Fabrycky & Tremaine 2007; Liu et al. 2014). However, there are subtleties in this problem. In the LK case the maximum eccentricity can actually be increased by GR precession if one includes octupole terms in the perturbing Hamiltonian (Ford et al. 2000; Naoz et al. 2013). In our (strictly quadrupole) case it turns out that increasing the strength of GR precession can increase $e_{\mathrm{max}}$ for some binaries if $\Gamma<1/5$, because of the way that GR precession modifies the phase space morphology. We defer a more careful exploration of the secular problem with the GR precession in all $\Gamma$ regimes to a future study.

9 Discussion

The main result of this work is the unveiling of a variety of new dynamical regimes that characterise the orbital evolution of a binary system subject to an external gravitational tidal field (‘cluster’). While the results of §2 are completely general, we found that we need to investigate dynamics in four separate regimes, corresponding to certain ranges of the parameter $\Gamma$ characterising the external tidal field of the cluster and the binary orbit in it. For $\Gamma>1/5$ (§3), the results were found to be qualitatively very similar to those previously derived in the test particle quadrupole Lidov-Kozai problem, which is recovered exactly by taking $\Gamma=1$ (Vashkov’yak, 1999; Kinoshita & Nakai, 2007; Antognini, 2015).

However, when leaving the regime $\Gamma>1/5$ several qualitative differences emerge. The condition for the existence of fixed points changes, as do locations of minimum and maximum eccentricities in the $(\omega,e)$ phase-space and the morphology of the $(D,\Theta)$ plane; even the very existence of the fixed points and orbits librating around them changes with $\Gamma$. The first bifurcation in the qualitative dynamics happens at $\Gamma=1/5$ but there are others at $\Gamma=0$ and $\Gamma=-1/5$, so we separately treated the regimes $0<\Gamma\leq 1/5$ (described in §4), $-1/5<\Gamma\leq 0$ (in §5), and $\Gamma\leq-1/5$ (in §6).

In Tables 1 & 2 we collect the results of §§3-6. In Table 1 we provide the locations and values of minimum/maximum $j$, the values of $\Delta$ which enter the secular timescale (33), and the allowed ranges of the constants of motion $D$ and $\Theta$, for each family of phase-space orbit and in each $\Gamma$ regime. Table 2 collects the locations and values of the extrema of the dimensionless Hamiltonian $H_{1}^{*}$, depending on the $\Gamma$ regime and whether or not fixed points exist.

In this work we never explicitly considered the possibility of $\Gamma>1$ — all our examples were given for $\Gamma\leq 1$. Situations in which $\Gamma$ exceeds unity are possible. However, we found in Paper I that this regime is realised only for rather extreme binary orbits inside the cluster, e.g. close to polar, which justifies our overall neglect of the $\Gamma>1$ possibility. Also, we should note that none of the results obtained in the $\Gamma>1/5$ regime (§3) explicitly assumed $\Gamma\leq 1$; they also hold when $\Gamma>1$. The only substantial difference in this case would be that $A$ becomes negative for $\Gamma>1$. As a result, maxima and minima of the (dimensional) perturbing Hamiltonian would swap their locations in the $(\omega,e)$ phase-space, and the phase-space trajectories would be traversed in the opposite direction compared to the $1/5<\Gamma\leq 1$ case.

For our doubly-averaged theory to properly characterise binary orbital evolution certain conditions should be met. We already saw in §7 that the description based on the doubly-averaged Hamiltonian (1) may fail when the secular timescale $t_{\rm sec}$ is not much longer than the period of the outer orbit of the binary $T_{\phi}$. In such cases one should resort to using the singly-averaged framework described in Paper I, which always works very well. Other phenomena that may affect the orbital evolution of binaries in clusters (on top of the smooth cluster tide-driven secular evolution) are discussed in §9.3.

9.1 Critical inclination for the existence of fixed points

A classic result of LK theory is the value of the critical initial inclination $i_{0}=i_{\mathrm{c}}$, above which fixed points appear in the $(\omega,e)$ phase-space. This critical angle marks the onset of large eccentricity oscillations, and a qualitative departure from classical Laplace-Lagrange dynamics. It provides a constraint on which binary orientations can lead to large eccentricity excursions.

Assuming an initial binary eccentricity of zero, we can calculate $i_{\mathrm{c}}$ for general $\Gamma$ using the conservation of $\Theta=\cos^{2}i_{0}$. The upper bounds on $\Theta$ for fixed points to exist ($=3/5$ in the test particle quadrupole LK case) are given by equations (37), (40) and (63); it then easily follows that the existence of fixed points requires $i_{0}>i_{\mathrm{c}}$ where

[TABLE]

with $\Lambda(\Gamma)$ given in equation (12). There are no fixed points for any initial inclination if $-1/5<\Gamma\leq 0$, see §5. If the initial eccentricity of the binary $e_{0}$ is non-zero, the argument of the $\cos^{-1}$ needs to be additionally divided by $\sqrt{1-e_{0}^{2}}$, further lowering $i_{\rm c}$.

We plot $i_{\mathrm{c}}$ as a function of $\Gamma$ in Figure 14. In the LK limit $\Gamma=1$ we recover the classic result $i_{\mathrm{c}}=\cos^{-1}\sqrt{3/5}\approx 39.2^{\circ}$. As we decrease $\Gamma$ from $1$, fixed points exist for ever smaller initial inclinations, until $i_{\mathrm{c}}$ reaches zero at $\Gamma=1/5$; note however that the secular timescale diverges as $|\Gamma|\to 1/5$. Asymptotes at $\Gamma=-1/5$ and $\Gamma=0$ ensure that fixed points never exist between those values. As $\Gamma\to-\infty$ we find $\Lambda\to 1/2$ and so $i_{\mathrm{c}}\to 45^{\circ}$.

9.2 High eccentricity behaviour

One is often interested in the high-eccentricity behaviour of binaries undergoing LK-like cycles, because it is at small pericentre distances that exotic effects like GR precession (8), mass transfer, gravitational wave emission and tidal circularisation become important. To explore these possibilities we consider orbits that are capable of reaching $e\to 1$ or $j\to 0$, which necessarily requires

[TABLE]

(see definitions 7) and study their behavior as they evolve through highest eccentricity. We will focus on orbits that start at eccentricities that are not too close to unity, and ignore the effects of GR precession (§8).

We base our discussion on equation (30). One of the roots ($j_{\pm},j_{0}$) corresponds to the smallest value of $j$ satisfying the constraint (9), which we have called $j_{\rm min}$. Normally, when (57) is true, the other two roots are much larger in magnitude. With this in mind equation (30) can be approximately integrated in the vicinity of this root $j_{\rm min}$ as

[TABLE]

where we defined

[TABLE]

(see equation (34)) and $j_{1}$ and $j_{2}$ are the other two roots of equation (30), i.e not $j_{\rm min}$ (which we normalized by its smallest possible value $\Theta^{1/2}$); $|j_{1}|$,$|j_{2}|\gg j_{\rm min}$. Time $t$ is counted from the point of reaching highest eccentricity, i.e. $j(t=0)=j_{\rm min}$, while $t_{\rm min}$ is the characteristic evolution timescale in the vicinity of $j_{\rm min}$ — the time it takes for $j$ to change from $j_{\rm min}$ (at $t=0$) to $2^{1/2}j_{\rm min}$. Note also that $j(t)\propto t$ when $j\gtrsim j_{\rm min}$.

The fact that the time spent near maximum eccentricity is of order $t_{\mathrm{min}}\propto j_{\mathrm{min}}=\sqrt{1-e_{\mathrm{max}}^{2}}$ is a familiar result from the Lidov-Kozai case (e.g. Miller & Hamilton 2002). It has been used to characterise the timescale for gravitational-wave induced mergers of binaries driven to high eccentricity through the LK effect (e.g. Thompson 2011; Antonini & Perets 2012; Bode & Wegg 2014; Liu & Lai 2018; Grishin et al. 2018; Randall & Xianyu 2018). We will employ the more general result (58)-(59) applicable for arbitrary axisymmetric perturbation (not just that of a point mass companion) when exploring the rate of compact-object binary mergers in stellar clusters predicted by our theory in future work.

Finally, phase space morphology (as determined by the value of $\Gamma$) is a crucial factor in determining how many phase space orbits are able to reach $e\to 1$. For example, provided (57) is satisfied, all orbits in the regime $\Gamma>1/5$ will reach very high eccentricities, whereas this ceases to be true for $0<\Gamma<1/5$ (compare Figures 4a,d,g with Figures 5d,g). This effect is important when calculating merger rates of compact object binaries in stellar clusters (Hamilton & Rafikov, 2019b).

9.3 Stellar scattering and other non-ideal effects

In our work we have assumed that the gravitational field of the cluster can be adequately approximated as time-independent. This is of course not true in general. For example, globular clusters in the Milky Way can be shocked and tidally stripped as they move through the Galactic disk. Also, globular clusters undergo secular evolution on $\gtrsim\mathrm{Gyr}$ timescales which eventually leads to core collapse. Both of these effects would directly modify the mean field potential $\Phi$ and could also alter a binary’s outer orbit dramatically.

Additionally, we have assumed that the cluster’s potential is perfectly smooth. However, one must remember that the cluster’s true potential is in fact the sum of the potentials of the many individual stars comprising it. As a result the true potential felt by the binary is both granular in space and stochastic in time. In practice, these effects can be explored by considering the impact of individual stellar passages in the vicinity of a binary on its orbital elements. The issue of binaries undergoing flyby encounters has been studied widely. Heggie & Rasio (1996) first considered the case of ‘secular encounters’, where the scattering event takes much longer than the orbital period of the inner binary (Hamers, 2018). This regime is appropriate if one is studying perturbations to the orbits of relatively tight systems (hard binaries), such as millisecond pulsars or hot Jupiters. On the other hand, Collins & Sari (2008) (see also Collins & Sari 2010) considered the opposite regime in which the timescale for the flyby interaction is much shorter than the inner binary period, so that the encounter can be treated in the impulse approximation. This is the correct description when studying the dynamics of the Oort Cloud comets in the Galaxy or very soft binaries in clusters. Finally, when the approach distance and velocity of the external perturber are comparable to the semi-major axis and the orbital speeds of the binary components, the binary changes its orbital elements in a dramatic fashion on a short (non-secular) timescale, with a high chance of being disrupted (Heggie, 1975; Goodman & Hut, 1993). This would completely reset the course of the smooth secular evolution of the binary orbit explored in this work. Thus, the prescription needed for estimating the effects of stellar scattering depends on the physical problem one wishes to address.

We defer a careful study of the coupling between the effects of stochastic stellar encounters and the smooth cluster tide-driven evolution of binaries to a future work. Here we simply estimate the characteristic time between close encounters of a binary with field stars. Assuming that all perturbers have mass $m$ and can be drawn from an homogeneous, isotropic Maxwellian distribution with number density $n$ and velocity disperision $\sigma$, we can estimate the typical time elapsed before the binary experiences a collision with impact parameter $q_{\mathrm{coll}}=a/2$ as (Binney & Tremaine, 2008):

[TABLE]

where

[TABLE]

is a measure of gravitational focusing, and we have used typical values of $n$ and $\sigma$ for a globular cluster (although as the binary moves through the cluster, the velocity dispersion and number density of field stars it experiences may change dramatically). One can see that depending on cluster mass, number density, binary semi-major axis, etc., $t_{\mathrm{coll}}$ can be larger or smaller than the secular timescale due to cluster tides (equations (33) and (34)). Moreover, weaker (secular) encounters (Heggie & Rasio, 1996) which cause slow random walk of the binary orbital elements would occur more frequently.

Thus, it is usually very important to take into account the effect of stellar flybys. However, we do not believe that this diminishes the astrophysical relevance of cluster tides, for several reasons. First, tidal effects can be important even in the outskirts of clusters where $n$ is low and stellar encounters are rare. In fact, tides can drive compact-object binaries to merge out to cluster-centric distances of several parsecs (Hamilton & Rafikov, 2019b). Second, in massive centrally cusped clusters (such as nuclear clusters with or without a central massive black hole), secular timescales due to tides can be as short as $\sim 10^{5}\mathrm{yr}$, potentially leading to interesting effects before close encounters occur. Third, while dense stellar environments can lead to frequent disruption of binaries they can also result in efficient binary formation, i.e. they can act as a source of new binaries that can then undergo tidal evolution.

9.4 Relation to previous work

The secular dynamics of binaries presented in this paper have been investigated thoroughly by other authors in the LK ($\Gamma=1$) limit (the ‘test particle quadrupole’ LK problem). In particular, Vashkov’yak (1999) and Kinoshita & Nakai (2007) derived analytically the maximum and minimum eccentricities and the timescale of LK oscillations. Antognini (2015) rederived the same results and provided an approximate fitting formula for the timescale.

A study of the phase-space portrait of binaries perturbed by the Galactic tide — a problem investigated by Heisler & Tremaine (1986) and many others, see Paper I — has been performed by Brasser (2001). Keeping only the $\partial^{2}\Phi/\partial Z^{2}$ contribution in the tidal expansion of the potential (equivalent to $\Gamma=1/3$), they derived the fixed points, secular timescale, and criteria for circulation and libration in $(\omega,e)$ space.

Petrovich & Antonini (2017) considered an extension to the LK problem in which a binary orbits a supermassive black hole (SMBH), and its (outer) orbit is perturbed by a non-spherical nuclear cluster potential (the inner orbit is assumed to be unperturbed by the cluster). Unlike our study, Petrovich & Antonini (2017) only looked at the effect of the cluster potential on the outer orbit of the binary and completely ignored the direct effect of the cluster potential on the secular dynamics of the inner orbital elements. Relevant for this work, part of their paper involves an investigation of the $(\omega,e)$ phase portrait of the inner binary in the quadrupole approximation, assuming (a) the outer orbit is almost circular and (b) the cluster potential is only weakly flattened. However, our doubly-averaged formalism does not cover this part of Petrovich & Antonini (2017)’s paper, because in this particular limit the outer orbit’s nodal precession timescale is long compared to the secular evolution time, so the perturbing potential cannot be considered axisymmetric (the situation here is similar to that described in §7.4).

10 Summary

We considered the secular dynamics of binaries arising from the general doubly-averaged tidal Hamiltonian derived in Paper I. Our study focused on exploring the phase portraits describing the evolution of binaries perturbed by the tidal field of a host cluster. We unravelled a number of new dynamical regimes previously not accounted for in studies of binary evolution problem, and provided their full classification. Our results can be briefly summarized as follows.

• •

We find that that under a wide range of initial conditions, a generic axisymmetric potential (including spherical potentials) can generate a sufficient tidal torque on a binary to allow it to perform large-amplitude secular eccentricity oscillations reminiscent of the LK mechanism.

• •

The morphology of the binary evolution in the phase-space of its orbital elements (e.g. $\omega$ and $e$) is uniquely set by the value of a single dimensionless parameter $\Gamma$, which encodes all information about the shape of the cluster potential and the binary orbit within it. We mapped out different dynamical behaviours of the binary as a function of $\Gamma$.

• •

Although the dynamics are qualitatively similar to the LK mechanism for $\Gamma>1/5$, there are bifurcations in the phase-space portrait when $\Gamma=\pm 1/5$ and [math] such that the dynamics become drastically different from LK case. We provide detailed descriptions of the binary evolution in each of the corresponding dynamical regimes.

• •

We numerically verify our theoretical predictions and find that they work well when the timescale for secular evolution is much longer than the time for the binary’s outer orbit to fill an axisymmetric torus inside the cluster. Such circumstances may be rare when $\Gamma<0$, because this regime typically requires strongly non-coplanar outer orbits that may take large number (several hundred) of orbital periods to fill a torus.

• •

General relativistic pericentre precession typically acts to quench secular eccentricity oscillations. Its effect can be easily included in our general doubly-averaged formalism.

• •

While the LK mechanism is efficient at driving high eccentricity oscillations, it requires the presence of a long-term distant companion to a binary. In contrast, every binary in a cluster feels the cluster potential, just as every comet feels the Galactic tide. As a result, the effect considered in this work, while possibly weaker than in the standard LK scenario, should be more ubiquitous in nature since it is available to any binary bound to an axisymmetric host system.

The theory we have developed in Papers I and II could be of importance to various astrophysical problems, such as formation of blue stragglers, X-ray binaries, hot Jupiters, and compact-object mergers (i.e. gravitational wave sources) in globular and nuclear star clusters. The possibility of the cluster tide-driven evolution explored here presents an interesting alternative to the well-explored hierarchical triple (LK) scenario.

In a third paper in this series we apply our formalism to the problem of binary compact-object mergers in globular and nuclear star clusters (Hamilton & Rafikov, 2019b). In subsequent work we will investigate the role of stochastic orbital element perturbations due to stellar flybys, among other things.

Acknowledgements

We thank Henrik Latter, Eugene Vasiliev and Antranik Sefilian for useful conversations. The orbit integration in §7 was done using galpy (http://github.com/jobovy/galpy, Bovy 2015). The N-body calculations were done with MERCURY (Chambers, 1999). CH is funded by a Science and Technology Facilities Council (STFC) studentship.

Appendix A Detailed characteristics of the $-1/5<\Gamma\leq 0$ regime.

Here we provide more details about the properties of the dynamical regime $-1/5<\Gamma\leq 0$, see §5.

A.1 Fixed points and orbit families

When $\Gamma=0$ the dimensionless Hamiltonian is simply $H_{1}^{*}=2+3e^{2}$, so the phase portrait consists of straight horizontal lines: all orbits circulate with $e(t)=e(0)$.

When $-1/5<\Gamma<0$, we use the constraint (13) to explore the possibility of fixed points. Since both $\Lambda$ and $\Lambda^{-1}$ are negative (see Figure 1), while $\Theta>0$, we conclude there are no fixed points in this $\Gamma$ regime. As a result, all orbits circulate, in agreement with Figure 6.

A.2 Range of parameter values

In the absence of fixed points in this $\Gamma$ regime, the only possible bounds on $D$ are $D=0$ and $D=1-\Theta$. Hence the $(D,\Theta)$ plane consists simply of a triangle of circulating orbits (see the third row of Figure 3):

[TABLE]

It is easy to show that the Hamiltonian is maximised at $j^{2}=\Theta$ (i.e. at the upper limit on eccentricity in Figure 6), so $H_{1,\mathrm{max}}^{*}$ obeys equation (22). Similarly it is minimised along the line of zero eccentricity ($j^{2}=1$), so $H_{1,\mathrm{min}}^{*}$ is given by equation (21).

A.3 Maximum and minimum eccentricities

Figure 6 shows that circulating orbits’ maximum and minimum eccentricities are back at $\omega=\pm\pi/2$ and $\omega=0$ respectively, as they were in the $\Gamma>1/5$ case (§3).

To understand the ordering of $j_{\pm}^{2},j_{0}^{2}$ we only need to consider panel (a) of Figure 2, because we have only circulating orbits in this $\Gamma$ regime. For $-1/5<\Gamma\leq 0$ the ordering of $j_{0}^{2}$ and $j_{+}^{2}$ has flipped compared to $0<\Gamma\leq 1/5$, while $j_{-}^{2}$ still lies outside of the physical region (i.e. it is not bounded by the horizontal dotted lines $j^{2}=\Theta,1$). Hence we must have $j_{-}^{2}<j_{+}^{2}<j^{2}<j_{0}^{2}$, so that $j_{\mathrm{min}}=j_{+}$, $j_{\mathrm{max}}=j_{0}$, and $\Delta=j_{0}^{2}-j_{-}^{2}$.

A.4 Timescales of eccentricity oscillations

The timescale $\log_{10}(t_{\mathrm{sec}}/t_{1})$ is plotted in the third row of Figure 3 for $\Gamma=-0.01,-0.1,-0.18$. We have only a triangle of circulating orbits, and their secular timescale is rather well approximated by $t_{1}$.

Appendix B Detailed characteristics of the $\Gamma\leq-1/5$ regime.

Here we provide more details about the properties of the dynamical regime $\Gamma\leq-1/5$, see §6.

B.1 Fixed points and orbit families

When $\Gamma=-1/5$ there are still no librating orbits, because $\Lambda(-1/5)=0$. However, librating orbits emerge as we decrease $\Gamma$ further. In terms of the constraint on $\Theta$ for fixed points and librating orbits to exist, the regime $\Gamma<-1/5$ mirrors the first $(\Gamma>1/5)$ regime in that we again require $\Theta<\Lambda$ (in Figure 1, the shaded region is bounded by the red curve):

[TABLE]

B.2 Range of parameter values

It is perhaps unsurprising from the morphology of the phase portraits (compare Figures 4 and 7) that the $(D,\Theta)$ plane for $\Gamma\leq-1/5$ (bottom row of Figure 3) looks similar to the $\Gamma>1/5$ case (top row of Figure 3). However, in this case the librating orbits are bounded to the left by $D_{+}$ (see equation (27)):

[TABLE]

As for the extrema of $H_{1}^{*}$, the only change from the case $-1/5<\Gamma\leq 0$ is that we now have fixed points available. If they exist then the Hamiltonian is minimised at the fixed point and so $H_{1,\mathrm{min}}^{*}$ obeys equation (23); if not, it is minimised at the zero eccentricity line $j^{2}=1$ (equation (21)). The maximum $H_{1,\mathrm{max}}^{*}$ is always found at $j=\sqrt{\Theta}$ and is therefore given by equation (22).

B.3 Maximum and minimum eccentricities

For circulating orbits we may again inspect Figure 2a. When $\Gamma\leq-1/5$ the physical solutions $j^{2}$ run from $j_{+}^{2}$ to $j_{0}^{2}$ as in the $-1/5<\Gamma<0$ case, but $j_{-}^{2}>1$ is suddenly larger than the others, so $j_{+}^{2}<j^{2}<j_{0}^{2}<j_{-}^{2}$. Thus $j_{\mathrm{min}}=j_{+}$, $j_{\mathrm{max}}=j_{0}$ and $\Delta=j_{-}^{2}-j_{+}^{2}$.

For librating orbits we read off from Figure 2b that $j_{+}^{2}<j^{2}<j_{-}^{2}<j_{0}^{2}$ (with $j_{0}>1$), giving $j_{\mathrm{min}}=j_{+}$, $j_{\mathrm{max}}=j_{-}$ and $\Delta=j_{0}^{2}-j_{+}^{2}$.

B.4 Timescales of eccentricity oscillations

We plot $\log_{10}(t_{\mathrm{sec}}/t_{1})$ for $\Gamma=-0.25,-0.5,-0.8$ in the bottom row of Figure 3. Bounds on $\Theta$ and $D$ are given by equations (63) and (64) respectively. The separatrix lies along $D=0$.

Along the separatrix the timescale for secular oscillations once again diverges. The timescale also diverges everywhere in $(D,\Theta)$ space for $\Gamma=-1/5$, see equation (33). However, as $\Gamma$ is lowered, one can see that $t_{\rm sec}$ becomes substantially smaller than $t_{1}$, just as in the case of $\Gamma\to 1$ considered in §3.4.