Secular chaotic dynamics in hierarchical quadruple systems, with applications to hot Jupiters in stellar binaries and triples

TL;DR

This paper investigates the complex secular chaotic dynamics in hierarchical quadruple star systems, revealing how precession can enhance eccentricity and influence hot Jupiter formation through high-eccentricity migration.

Contribution

It provides a detailed analysis of the secular dynamics in hierarchical quadruple systems, introducing a simplified model for eccentricity excitation due to outer orbit precession.

Findings

Eccentricity excitation is enhanced when precession period matches Lidov-Kozai timescale.

Outer orbit precession induces large mutual inclination, even from small initial inclinations.

A semi-analytic model explains the eccentricity amplification mechanism.

Abstract

Hierarchical quadruple systems arise naturally in stellar binaries and triples that harbour planets. Examples are hot Jupiters (HJs) in stellar triple systems, and planetary companions to HJs in stellar binaries. The secular dynamical evolution of these systems is generally complex, with secular chaotic motion possible in certain parameter regimes. The latter can lead to extremely high eccentricities and, therefore, strong interactions such as efficient tidal evolution. These interactions are believed to play an important role in the formation of HJs through high-eccentricity migration. Nevertheless, a deeper understanding of the secular dynamics of these systems is still lacking. Here we study in detail the secular dynamics of a special case of hierarchical quadruple systems in either the `2+2' or `3+1' configurations. We show how the equations of motion can be cast in a form representing a perturbed hierarchical three-body system, in which the outer orbital angular momentum vector is precessing steadily around a fixed axis. In this case, we show that eccentricity excitation can be significantly enhanced when the precession period is comparable to the Lidov-Kozai oscillation time-scale of the inner orbit. This arises from an induced large mutual inclination between the inner and outer orbits driven by the precession of the outer orbit, even if the initial mutual inclination is small. We present a simplified semi-analytic model that describes the latter phenomenon.

Equations

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

$$\frac{\mathrm{d}{\mathbf{e}}_1}{\mathrm{d}t}$$

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

$$\frac{\mathrm{d}{\mathbf{e}}_2}{\mathrm{d}t}$$

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

$$\frac{\mathrm{d}{\mathbf{e}}_3}{\mathrm{d}t}$$

$$+\frac{3}{8}{t}_{\mathrm{LK},2+2,23}^{-1}\frac{{\Lambda }_2}{{\Lambda }_3}\frac{1}{\sqrt{1-e_3^2}}\left[\left\{\left(1-6e_2^2\right)+25{\left({\mathbf{e}}_2\cdot {\hat{\mathbf{L}}}_3\right)}^2-5{\left({\mathbf{j}}_2\cdot {\hat{\mathbf{L}}}_3\right)}^2\right\}\left({\mathbf{e}}_3\times {\hat{\mathbf{L}}}_3\right)-10\left({\mathbf{e}}_2\cdot {\hat{\mathbf{L}}}_3\right)\left({\mathbf{e}}_3\times {\mathbf{e}}_2\right)+2\left({\mathbf{j}}_2\cdot {\hat{\mathbf{L}}}_3\right)\left({\mathbf{e}}_3\times {\mathbf{j}}_2\right)\right]$$

$$\equiv {\mathbf{e}}_3\times \mathbf{f}({\mathbf{e}}_1,{\mathbf{j}}_1,{\mathbf{e}}_2,{\mathbf{j}}_2,{\mathbf{j}}_3).$$

$$t_{\mathrm{LK},2+2,13}$$

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

$$\frac{\mathrm{d}{\hat{\mathbf{L}}}_3}{\mathrm{d}t}$$

$${\Omega }_3$$

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

$$\frac{\mathrm{d}{\mathbf{e}}_1}{\mathrm{d}t}$$

$$-5\left({\mathbf{e}}_1\cdot {\hat{\mathbf{L}}}_3\right)\left({\mathbf{j}}_1\times {\hat{\mathbf{L}}}_3\right)];$$

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

$$\frac{\mathrm{d}{\mathbf{e}}_2}{\mathrm{d}t}$$

$$+\frac{3}{4}{t}_{\mathrm{LK},23}^{-1}\left[\left({\mathbf{j}}_2\cdot {\hat{\mathbf{L}}}_3\right)\left({\mathbf{e}}_2\times {\hat{\mathbf{L}}}_3\right)+2\left({\mathbf{j}}_2\times {\mathbf{e}}_2\right)-5\left({\mathbf{e}}_2\cdot {\hat{\mathbf{L}}}_3\right)\left({\mathbf{j}}_2\times {\hat{\mathbf{L}}}_3\right)\right];$$

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

$$\frac{\mathrm{d}{\mathbf{e}}_3}{\mathrm{d}t}$$

$$+\frac{3}{8}{t}_{\mathrm{LK},3+1,13}^{-1}\frac{{\Lambda }_1}{{\Lambda }_3}\frac{1}{\sqrt{1-e_3^2}}\left[\left\{\left(1-6e_1^2\right)+25{\left({\mathbf{e}}_1\cdot {\hat{\mathbf{L}}}_3\right)}^2-5{\left({\mathbf{j}}_1\cdot {\hat{\mathbf{L}}}_3\right)}^2\right\}\left({\mathbf{e}}_3\times {\hat{\mathbf{L}}}_3\right)-10\left({\mathbf{e}}_1\cdot {\hat{\mathbf{L}}}_3\right)\left({\mathbf{e}}_3\times {\mathbf{e}}_1\right)+2\left({\mathbf{j}}_1\cdot {\hat{\mathbf{L}}}_3\right)\left({\mathbf{e}}_3\times {\mathbf{j}}_1\right)\right]$$

$$\equiv {\mathbf{e}}_3\times \mathbf{f}({\mathbf{e}}_1,{\mathbf{j}}_1,{\mathbf{e}}_2,{\mathbf{j}}_2,{\mathbf{j}}_3).$$

$$t_{\mathrm{LK},3+1,12}$$

$$t_{\mathrm{LK},3+1,13}$$

$$\frac{\mathrm{d}{\hat{\mathbf{L}}}_2}{\mathrm{d}t}$$

$${\Omega }_2$$

$$\frac{\mathrm{d}\mathbf{j}}{\mathrm{d}\tau }$$

$$\frac{\mathrm{d}\mathbf{e}}{\mathrm{d}\tau }$$

$$\frac{\mathrm{d}{\hat{\mathbf{L}}}_{\mathrm{out}}}{\mathrm{d}\tau }$$

$$t_{\mathrm{LK}}=\frac{m_0+m_1}{m_{\mathrm{out}}}\sqrt{\frac{a^3}{G(m_0+m_1)}}{\left(\frac{a_{\mathrm{out}}}{a}\right)}^3{\left(1-e_{\mathrm{out}}^2\right)}^{3/2}.$$

$$\beta \equiv {\Omega }_{\mathrm{out}}{t}_{\mathrm{LK}}.$$

$$H_{\mathrm{gen},\mathrm{rot}}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\left[\sqrt{1-e^2}\cos (i)\right]=-{\Omega }_{\mathrm{out}}\mathbf{j}\cdot \left(\hat{\mathbf{z}}\times {\hat{\mathbf{L}}}_{\mathrm{out}}\right).$$

$$e_{\max }=\sqrt{1-\frac{5}{3}{\cos }^2(i_0)},$$

$$0=\frac{\mathrm{d}e}{\mathrm{d}\tau }=-\frac{15}{4}\left(\mathbf{e}\cdot {\hat{\mathbf{L}}}_{\mathrm{out}}\right)\left[\hat{\mathbf{e}}\cdot \left(\mathbf{j}\times {\hat{\mathbf{L}}}_{\mathrm{out}}\right)\right].$$

$${\left(\hat{\mathbf{e}}\cdot {\hat{\mathbf{L}}}_{\mathrm{out}}\right)}^2=1-{\left(\hat{\mathbf{L}}\cdot {\hat{\mathbf{L}}}_{\mathrm{out}}\right)}^2,$$

$${\hat{\mathbf{L}}}_{\mathrm{out}}=\cos (\varphi )\sin (\alpha )\hat{\mathbf{x}}+\sin (\varphi )\sin (\alpha )\hat{\mathbf{y}}+\cos (\alpha )\hat{\mathbf{z}},$$

$$⟨(...)⟩\equiv \frac{1}{2\pi }{\int }_0^{2\pi }\mathrm{d}\varphi (...).$$

$$⟨\left[\hat{\mathbf{e}}\cdot \left(\mathbf{j}\times {\hat{\mathbf{L}}}_{\mathrm{out}}\right)\right]⟩$$

$$⟨{\left(\hat{\mathbf{e}}\cdot {\hat{\mathbf{L}}}_{\mathrm{out}}\right)}^2⟩$$

Full text

††pubyear: 2017††pagerange: Secular chaotic dynamics in hierarchical quadruple systems, with applications to hot Jupiters in stellar binaries and triples–10

1 Introduction

Approximately 1 per cent of all stellar FG dwarf systems are hierarchical quadruples (Tokovinin, 2014a, b). In addition to these purely stellar systems, hierarchical quadruple configurations also occur naturally in lower multiplicity stellar systems that harbour planets. For example, there are currently three hot Jupiters (HJs; Jupiter-like planets orbiting stars in several days) known in stellar triple systems, i.e., WASP-12b (Hebb et al., 2009; Bergfors et al., 2013; Bechter et al., 2014), HAT-P-8b (Latham et al., 2009; Bergfors et al., 2013; Bechter et al., 2014) and KELT-4Ab (Eastman et al., 2016). In these systems, the HJ and its host star are orbited by a stellar binary (see the left-hand panel of Fig. 1). Such a binary may have played a role in the formation and migration of the proto-HJ.

In particular, as shown by Hamers (2017a), the ‘binarity’ of the companion can introduce secular enhancement of the eccentricity of the proto-HJ orbit in a larger parameter space compared to the situation when the star+HJ system is orbited by a single star. In the latter case, the eccentricity excitation is driven by Lidov-Kozai (LK) oscillations (Lidov, 1962; Kozai, 1962) that arise in hierarchical three-body systems. Such an enhancement of the eccentricity excitation of proto-HJs in stellar triples compared to stellar binaries is relevant, because high-eccentricity migration models of HJs in stellar binaries (Wu & Murray, 2003; Fabrycky & Tremaine, 2007; Naoz et al., 2012; Petrovich, 2015a; Anderson et al., 2016; Muñoz et al., 2016) are faced with the problem that the predicted formation rates are about 5-10 times lower than observed. One of the reasons for the lower rates is that short-range force precession in the orbit of the proto-HJ suppresses secular excitation if the orbit of the stellar companion is relatively wide (Ngo et al., 2016).

Related to the above, the efficiency of high-eccentricity migration in stellar binaries could be enhanced if there are (currently undetected) massive planetary companions to HJs in stellar binaries, orbiting in-between the HJ and the stellar binary companion (see the right-hand panel of Fig. 1). If the planetary companion satisfies several constraints, it can mediate LK-like oscillations in the proto-HJ orbit induced by the stellar companion, even if the planetary companion was initially coplanar with respect to the proto-HJ. This would put further constraints on high-eccentricity migration if, in the future, such currently unseen companions to HJs in stellar binaries are found to be absent (Hamers, 2017b).

The two configurations of HJs discussed above can be classified as ‘2+2’ and ‘3+1’ quadruple systems (see Fig. 1). Although the orbit-averaged Hamiltonian and the equations of motion for these systems (and higher multiplicity systems) are known (Hamers et al. 2015; Vokrouhlický 2016; Hamers & Portegies Zwart 2016; see also, e.g., Petrovich 2015b; Liu et al. 2015 for vector-form equations for triples), a deeper understanding of the underlying mechanism for the enhanced eccentricity excitations, and the associated chaotic dynamics, is currently lacking.

In this paper, we study in detail the secular dynamics of low-mass objects (planets or HJs) in quadruple systems with three more massive bodies. We show that, with a number of additional assumptions, the equations of motion for the ‘2+2’ and ‘3+1’ configurations can be cast into a single, general, form that is mathematically identical to that of a hierarchical three-body system in which the outer orbital angular momentum vector is precessing steadily around a fixed axis (Section 2). Although the latter model is not amenable to analytic solutions, we show qualitatively how eccentricity excitation arises in this model (Section 3). In addition, we present another simplified model for which analytic results can be obtained, and which provides physical insight (Section 4). We conclude in Section 6.

Nearing the completion of this paper, we became aware of the simultaneous work of Petrovich & Antonini (2017), who discuss similar dynamics of perturbed hierarchical three-body systems in a different context. Petrovich & Antonini (2017) consider binaries embedded in a non-spherical nuclear star cluster, and find that extreme eccentricity excitation is possible if the LK time-scale associated with the torque of the central massive black hole is comparable to the nodal precession time-scale of the binary centre of mass associated with the nuclear star cluster.

2 Model

In this section, we present a simplified model for the secular dynamics of hierarchical quadruple systems in which a test particle is orbiting in a system with three more massive bodies. In hierarchical triple systems, a useful approximation is the quadrupole-order test particle limit in which the angular momentum of the inner orbit, $L_{\mathrm{in}}$, is negligible compared to the angular momentum of the outer orbit, $L_{\mathrm{out}}$. In this limit, the outer orbital angular momentum vector $\boldsymbol{L}_{\mathrm{out}}$ is constant (therefore, $e_{\mathrm{out}}$ is constant as well), and the system is completely integrable. This system gives rise to well-known LK oscillations, occurring if $\boldsymbol{L}_{\mathrm{in}}$ and $\boldsymbol{L}_{\mathrm{out}}$ are initially inclined by more than $i_{\mathrm{crit}}\equiv\mathrm{arccos}(\sqrt{3/5})\approx 39.2315^{\circ}$, and with a maximum eccentricity of $e_{\mathrm{max}}=[1-(5/3)(\hat{\boldsymbol{L}}_{\mathrm{in,init}}\cdot\hat{\boldsymbol{L}}_{\mathrm{out}})^{2}]^{1/2}$, where $\hat{\boldsymbol{L}}_{\mathrm{in,init}}$ is the initial unit inner orbital angular momentum vector (assuming a zero initial eccentricity).

Our model applies to hierarchical quadruple systems, but we will show that it can be considered as a perturbed hierarchical three-body problem in the test particle limit. In our case, $\boldsymbol{L}_{\mathrm{out}}$ (of the perturbed three-body problem) is no longer constant but precesses steadily around a fixed axis. With the introduction of this perturbation, the system is no longer integrable, and the evolution is generally more complex. In particular, chaotic secular behaviour can be induced. Below, we will show that this perturbed model applies to two types of restricted hierarchical quadruple systems (the ‘2+2’ and ‘3+1’ configurations). Both cases have direct applications to planetary systems.

2.1 2+2 quadruple systems

Consider hierarchical quadruple systems in the ‘2+2’ configuration, i.e., two binaries orbiting each other’s barycenter. The hierarchy and notation are indicated schematically in the left-hand panel of Fig. 1. To the quadruple-order, i.e., the second order in the ratios of the orbital separations, the orbit-averaged Hamiltonian was derived by Hamers et al. (2015) and Hamers & Portegies Zwart (2016), and consists of the hierarchical three-body Hamiltonian applied to the (1,3) orbit pair, plus the hierarchical three-body Hamiltonian applied to the (2,3) orbit pair. The equations of motion for the eccentricity $\boldsymbol{e}_{i}$ and angular-momentum $\boldsymbol{j}_{i}=\sqrt{1-e_{i}^{2}}\,\hat{\boldsymbol{L}}_{i}$ vectors for the three orbits follow from the Milankovitch equations (Milankovitch 1939; Musen 1961; Allan & Ward 1963; Allan & Cook 1964; Breiter & Ratajczak 2005; Tremaine et al. 2009), and read

[TABLE]

Here, $\Lambda_{i}$ is the angular momentum of orbit $i$ for a circular orbit which is constant in the secular approximation, i.e., $\Lambda_{i}=\mu_{i}\sqrt{GM_{i}a_{i}}$, where $\mu_{i}$ and $M_{i}$ are the reduced and the total mass, respectively, of binary $i$. The (two) LK time-scales are given by

[TABLE]

The function $\boldsymbol{f}$ in equation (1) is independent of $\boldsymbol{e}_{3}$, showing that $\boldsymbol{e}_{3}$ precesses around $\hat{\boldsymbol{f}}$ and the magnitude of $e_{3}$ remains constant. Therefore, the LK time-scales in equations (2) are constant as well.

The coupled equations (1) are generally not amenable to analytic solutions. We simplify them by assuming the test particle limit, $\Lambda_{1}\ll\Lambda_{2},\Lambda_{3}$, and setting $\boldsymbol{e}_{2}=\boldsymbol{0}$. Evidently, the torque of orbit 3 can excite $e_{2}$ if orbits 2 and 3 are initially sufficiently inclined. Here, we assume that $i_{\mathrm{23,init}}\lesssim 40^{\circ}$, such that the LK mechanism is not active for the orbit pair (2,3); therefore, $e_{2}$ is constant. The assumption on the $\Lambda_{i}$ implies that the first term in equation (1e), proportional to $\Lambda_{1}/\Lambda_{3}$, is negligible compared to the second term, which is proportional to $\Lambda_{2}/\Lambda_{3}$ (also assuming that $t_{\mathrm{LK,2+2,13}}$ and $t_{\mathrm{LK,2+2,23}}$ are not too distinct). Equation (1e) then reads

[TABLE]

Using that the total angular momentum vector, $\boldsymbol{L}_{\mathrm{tot}}=\boldsymbol{L}_{1}+\boldsymbol{L}_{2}+\boldsymbol{L}_{3}$, is conserved, equation (3) can be written in the form

[TABLE]

where $\boldsymbol{\Omega}_{3}$ is a constant vector with magnitude

[TABLE]

with $i_{23,\mathrm{init}}$ as the initial inclination between $\hat{\boldsymbol{L}}_{2}$ and $\hat{\boldsymbol{L}}_{3}$.

In summary, the secular dynamics of the restricted problem, $\Lambda_{1}\ll\Lambda_{2},\Lambda_{3}$ and $\boldsymbol{e}_{2}=\boldsymbol{0}$, are described by equations (1a), (1b) and (4). These restricted equations apply, e.g., to a planet orbiting a star (orbit 1) that is orbited by a more distant stellar binary (orbit 2), in a relatively wide orbit (orbit 3).

2.2 3+1 quadruple systems

Next, we consider hierarchical quadruple systems in the ‘3+1’ configuration, i.e., a triple orbited by a distant fourth body (see the second panel in Fig. 1). To quadrupole order, the orbit-averaged Hamiltonian is given by adding the relevant Hamiltonians from the hierarchical three-body Hamiltonian, i.e., the three Hamiltonians associated with the (1,2), (2,3) and (1,3) pairs (Hamers et al., 2015; Hamers & Portegies Zwart, 2016). The equations of motion read

[TABLE]

The (three) LK time-scales are now given by

[TABLE]

Again, the function $\boldsymbol{f}$ in equation (6f) is independent of $\boldsymbol{e}_{3}$, showing that $\boldsymbol{e}_{3}$ precesses around $\hat{\boldsymbol{f}}$ and $e_{3}$ is constant. In this case, the LK time-scale $t_{\mathrm{LK,3+1,12}}$ is generally not constant because $e_{2}$ is not generally constant.

We make the assumption $\Lambda_{1}\ll\Lambda_{2}\ll\Lambda_{3}$. Furthermore, we assume $\boldsymbol{e}_{2}=\boldsymbol{0}$. As before, we assume that $i_{\mathrm{23,init}}\lesssim 40^{\circ}$ such that the LK mechanism does not operate for the orbit pair (2,3). Equation (6c) can then be written in the form

[TABLE]

where $\hat{\boldsymbol{\Omega}}_{2}=\hat{\boldsymbol{L}}_{3}$ is constant, and

[TABLE]

In addition, for hierarchical systems ($a_{3}\gg a_{2}$) with a sufficiently large ratio $m_{2}/m_{3}$, $t_{\mathrm{LK,3+1,12}}\ll t_{\mathrm{LK,3+1,13}}$, i.e., the torque of orbit 3 on orbit 1 is negligible compared to the torque of orbit 2 on orbit 1. The restricted problem is then defined by equations (6a) and (6b), dropping terms proportional to $t_{\mathrm{LK,3+1,13}}^{-1}$, and equation (8). It applies, e.g., to a massive planetary companion to a proto-HJ in a stellar binary.

2.3 Generalized model

2.3.1 Equations of motion

Although they correspond to different hierarchical configurations, the restricted equations of motion for the ‘2+2’ and ‘3+1’ configurations are mathematically identical. They can be interpreted as representing a perturbed hierarchical three-body system in the test particle limit, where the outer angular momentum vector is precessing around a fixed axis with a fixed rate $\Omega_{\mathrm{out}}$. The model can be written in the scaled form

[TABLE]

Here, $\boldsymbol{e}$ and $\boldsymbol{j}$ correspond to the inner test-particle orbit, and $\tau\equiv t/t_{\mathrm{LK}}$ with the LK time-scale given by

[TABLE]

The semimajor axis of the test-particle orbit is denoted with $a$, and the subscript ‘out’ denotes properties of the ‘outer’ companion, which depends on the hierarchical configuration (in the ‘2+2’ case, $a_{\mathrm{out}}=a_{3}$, $e_{\mathrm{out}}=e_{3}$ and $m_{\mathrm{out}}=m_{2}+m_{3}$, whereas $a_{\mathrm{out}}=a_{2}$, $e_{\mathrm{out}}=e_{2}$ and $m_{\mathrm{out}}=m_{2}$ in the ‘3+1’ case). The dimensionless and constant parameter $\beta$ is defined as

[TABLE]

In equation (10c), the fixed axis around which $\hat{\boldsymbol{L}}_{\mathrm{out}}$ precesses is taken to be the $\hat{\boldsymbol{z}}$ axis. We denote the angle between $\hat{\boldsymbol{L}}_{\mathrm{out}}$ and $\hat{\boldsymbol{z}}$ with $\alpha$.

2.3.2 Hamiltonian and the ‘LK constant’

If $\beta=0$ in equations (10), then the (conserved) Hamiltonian is mathematically equivalent to the three-body test-particle Hamiltonian. For nonzero $\beta$, this Hamiltonian is not constant because of the fixed precession of $\hat{\boldsymbol{L}}_{\mathrm{out}}$. However, one can transform to a frame that is rotating around the $\hat{\boldsymbol{z}}$ axis with the same rate at which $\hat{\boldsymbol{L}}_{\mathrm{out}}$ is precessing around the $\hat{\boldsymbol{z}}$ axis. This gives (e.g., Kinoshita 1993)

[TABLE]

where $H_{0}$ is a constant.

In the test-particle three-body limit ($\beta=0$), the quantity $\sqrt{1-e^{2}}\cos(i)$, where $i$ is the mutual inclination between the inner and outer orbits, is a constant of the motion (known as the ‘LK’ or ‘Kozai’ constant). For $\beta\neq 0$, the ‘LK constant’ is no longer conserved, but evolves according to

[TABLE]

3 Eccentricity excitation in the generalized model

The secular three-body equations of motion in the quadrupole-order test particle limit are integrable, i.e., analytic solutions exist for the eccentricity as a function of time (Kinoshita & Nakai, 2007). However, they are no longer amenable to simple solutions if the outer angular momentum vector is precessing as in our model described by equations (10). Of course, it is possible to numerically integrate equations (10), which is done here.

The initial conditions for the numerical integrations are as follows. The vectors $\hat{\boldsymbol{L}}$, $\hat{\boldsymbol{L}}_{\mathrm{out}}$ and $\hat{\boldsymbol{z}}$ are initially all assumed to lie in the same plane; the latter plane is the $(x,z)$-plane. In addition, $\alpha$, the angle about which $\hat{\boldsymbol{L}}_{\mathrm{out}}$ precesses around $\hat{\boldsymbol{z}}$, is defined such that positive $\alpha$ corresponds to $\hat{\boldsymbol{L}}_{\mathrm{out}}$ having a positive component along the $x$ axis. The initial mutual inclination between $\hat{\boldsymbol{L}}$ and $\hat{\boldsymbol{L}}_{\mathrm{out}}$ is denoted with $i_{0}$. The initial mutual inclination between $\hat{\boldsymbol{L}}$ and $\hat{\boldsymbol{z}}$ is $i_{z,0}=i_{0}+\alpha$. The eccentricity vector $\boldsymbol{e}$ is assumed to be initially parallel to the $\hat{\boldsymbol{y}}$ axis, and its initial magnitude is set to $e=0.01$.

We carried out sets of 1600 integrations on a linear grid with $\beta\equiv\Omega_{\mathrm{out}}\,t_{\mathrm{LK}}$ running from 0 to 10, and the initial $i$, $i_{0}$, running from $0^{\circ}$ to $89^{\circ}$. Here, we set $\alpha=5^{\circ}$ or $30^{\circ}$ for each set. The duration of each integration was $\Delta\tau=1000$, i.e., corresponding to a physical time-span of $1000\,t_{\mathrm{LK}}$, and approximately 1000 LK oscillations if $\hat{\boldsymbol{L}}_{\mathrm{out}}$ were fixed. During the integrations, local maxima of $e$ were recorded, and the latter are shown as a function of $i_{0}$ in Fig. 2 for $\alpha=5^{\circ}$ (results for $\alpha=30^{\circ}$ are given in Appendix A). Each of the 20 panels corresponds to a different value of $\beta$ indicated therein; the number of bins in $i_{0}$ is 80.

Below, we discuss the maximum eccentricity behaviour as a function of $i_{0}$ for three regimes of $\beta$: small $\beta$ (Section 3.1), large $\beta$ (Section 3.2) and intermediate $\beta$ (Section 3.3). Short-range forces are discussed briefly in Section 3.4. In Section 4, we present a simplified model for the mutual inclination, and use this to approximately describe the behaviour in the regime of intermediate $\beta$.

3.1 Small $\beta$

In each panel of Fig. 2, the black dashed curve shows the canonical expression

[TABLE]

which applies to the ‘unperturbed’ problem with $\beta=0$, and assuming initially $e=0$ and $\cos(i_{0})<\sqrt{3/5}$. This result, which is described here in detail for further reference, can be obtained by dotting equation (10b) with $\hat{\boldsymbol{e}}$, giving the stationary points, i.e.,

[TABLE]

Assuming that at the stationary point $e\neq 0$, $e\neq 1$ and $\hat{\boldsymbol{e}}\cdot\hat{\boldsymbol{L}}_{\mathrm{out}}\neq 0$, this can be rewritten as the condition

[TABLE]

where we used a vector identity for the scalar product of two triple products. Using equation (17) to eliminate $\hat{\boldsymbol{e}}\cdot\hat{\boldsymbol{L}}_{\mathrm{out}}$, and using the LK constant to express $\hat{\boldsymbol{L}}\cdot\hat{\boldsymbol{L}}_{\mathrm{out}}$ in terms of the initial inclination, conservation of the Hamiltonian yields an algebraic equation for $e$ with the solution equation (15) for the (maximum) stationary point.

As shown in Fig. 2, for the smallest non-zero values of $\beta$, equation (15) still gives a good description as expected. For $\beta\gtrsim 0.2$ but still $\ll 1$, noticeable deviations start to occur; in particular, large eccentricity excitations, much larger than based on equation (15), are possible for inclinations near $90^{\circ}$.

3.2 Large $\beta$

If $\beta$ is large ($\beta\gg 1$), then $\hat{\boldsymbol{L}}_{\mathrm{out}}$ is precessing rapidly around the $\hat{\boldsymbol{z}}$ axis compared to the time-scale at which $\hat{\boldsymbol{L}}$ evolves (i.e., on a time-scale on the order of $t_{\mathrm{LK}}$). Intuitively, one might expect in this case that the rapid precession of $\hat{\boldsymbol{L}}_{\mathrm{out}}$ effectively means that $\hat{\boldsymbol{L}}_{\mathrm{out}}$ is pointing along $\hat{\boldsymbol{z}}$, although with a modified ‘effective length’ depending on $\alpha$, the angle at which $\hat{\boldsymbol{L}}_{\mathrm{out}}$ precesses around $\hat{\boldsymbol{z}}$.

In the $\beta\gg 1$ limit, we can average the Hamiltonian equation (13) over a precessional cycle of $\hat{\boldsymbol{L}}_{\mathrm{out}}$. To achieve this, we write

[TABLE]

where $\phi\in[0,2\pi)$ is the (rapidly changing) precessional phase angle. The precession average of a quantity $(...)$ is defined by

[TABLE]

A number of useful averages are

[TABLE]

Averaging equation (14) gives

[TABLE]

In other words, in the limit $\beta\gg 1$, the LK constant, $\boldsymbol{j}\cdot\hat{\boldsymbol{L}}_{\mathrm{out}}$, is replaced by $\boldsymbol{j}\cdot\hat{\boldsymbol{z}}$. Owing to the existence of this conserved quantity, the maximum eccentricity can be obtained analytically, as described below.

Following the same steps as in Section 3.1, the condition for a stationary eccentricity (still) reads $\left[\hat{\boldsymbol{e}}\cdot\left(\boldsymbol{j}\times\hat{\boldsymbol{L}}_{\mathrm{out}}\right)\right]=0$. From equation (20a), it follows that $\left[\hat{\boldsymbol{e}}\cdot\left(\boldsymbol{j}\times\hat{\boldsymbol{z}}\right)\right]=0$, which implies $\left(\hat{\boldsymbol{e}}\cdot\hat{\boldsymbol{z}}\right)^{2}=1-\left(\hat{\boldsymbol{L}}\cdot\hat{\boldsymbol{z}}\right)^{2}$.

Next, we substitute these results into the Hamiltonian equation (13). By equation (21), the last term $\propto\boldsymbol{j}\cdot\hat{\boldsymbol{z}}$ in the Hamiltonian is constant and can therefore be omitted. Dropping the latter term, using $\left(\hat{\boldsymbol{e}}\cdot\hat{\boldsymbol{z}}\right)^{2}=1-\left(\hat{\boldsymbol{L}}\cdot\hat{\boldsymbol{z}}\right)^{2}$ to eliminate $\hat{\boldsymbol{e}}\cdot\hat{\boldsymbol{z}}$ and equation (21) to relate $\hat{\boldsymbol{L}}\cdot\hat{\boldsymbol{z}}$ to the initial value, $(\hat{\boldsymbol{L}}\cdot\hat{\boldsymbol{z}})_{0}=\cos(i_{z,0})$, and applying equations (20b) and (20c), the Hamiltonian at the stationary eccentricity in the limit $\beta\gg 1$ reads

[TABLE]

Equating (3.2) to the initial averaged Hamiltonian with $e=0$,

[TABLE]

we can solve analytically for the maximum eccentricity as a function of $i_{z,0}$ and $\alpha$. The solution is

[TABLE]

which is, remarkably, independent of $\alpha$ (if expressed in terms of $i_{z,0}$). Equation (23) is simply equation (15) with $i_{0}$ replaced by $i_{z,0}$. In other words, in the limit of very rapid precession, one obtains the classical result for the maximum eccentricity with the initial mutual inclination now replaced by the initial inclination between $\hat{\boldsymbol{L}}$ and the $\hat{\boldsymbol{z}}$ axis.

In Fig. 2, we plot equation (23) with the solid black lines. The numerical integrations indeed agree well with equation (23) for large $\beta$ ($\beta\gtrsim 6$). Note that the largest $e_{\mathrm{max}}$ does not occur at $i_{0}=90^{\circ}$, but at $i_{0}=90^{\circ}-5^{\circ}=85^{\circ}$, since $\alpha=5^{\circ}$.

As an additional test of equation (23), we show in Fig. 3 with blue points the maximum eccentricities as a function of $\alpha$ determined by numerically integrating the equations of motion, where we set $i_{0}=60^{\circ}$ and $\beta=1000$. The numerical results are in good agreement with equation (23), shown with the solid red line.

3.3 Intermediate $\beta$

In the intermediate-$\beta$ regime ($0.2\lesssim\beta\lesssim 6$), $e_{\mathrm{max}}$ as a function of $i_{0}$ is much more complex. In particular, for $\beta$ around unity, extremely large eccentricity excitations are possible. Whereas in the limit $\beta=0$, $e_{\mathrm{max}}\rightarrow 1$ only for $i_{0}\rightarrow 90^{\circ}$, if $\beta\sim 1$, $e_{\mathrm{max}}\rightarrow 1$ for a much larger range of $i_{0}$ with $1-e_{\mathrm{max}}$ in the integrations reaching values lower than $10^{-6}$. The precise range of $i_{0}$ for extreme eccentricity excitation depends on $\beta$ (and $\alpha$; compare, Figures 2 and 10). For lower $\beta$, $0.2\lesssim\beta\lesssim 0.5$, and $\alpha=5^{\circ}$, there are specific inclinations $i_{0}$ for which $e_{\mathrm{max}}\rightarrow 1$, indicative of an overlap of resonances, and, therefore, chaos (Chirikov, 1979). For larger $\beta$, $0.5\lesssim\beta\lesssim 3$, there is a broad range of $i_{0}$ for which extreme eccentricities are reached. This range of local maximum eccentricities is manifested as a band in the $(i_{0},1-e_{\mathrm{max}})$ plane. The upper boundary on $e_{\mathrm{max}}$ of this band appears to be smooth function of $i_{0}$. A (qualitative) understanding of this boundary is useful for practical applications, since it describes the largest possible eccentricity (i.e., over all local values of $e_{\mathrm{max}}$) for a given $i_{0}$.

In the chaotic regime of intermediate $\beta$, the eccentricity generally evolves in a complicated way, with local maxima occurring at different values, and orbital flips occurring frequently. A number of examples of the inclination and eccentricity as a function of time for $\beta\approx 1.44$ and $\alpha=5^{\circ}$ are given in Fig. 4. In particular, for $i_{0}\approx 34^{\circ}$, there is significant eccentricity excitation even though $i_{0}$ is less than the critical LK angle. For $i_{0}\approx 79^{\circ}$, orbital flips occur, and they are associated with extremely high eccentricities.

In order to understand the underlying mechanism for enhanced eccentricity excitation for a given $i_{0}$, we explore the connection between the maximum inclination, $i_{\mathrm{max}}$, and the maximum eccentricity. This is motivated by the notion that, if the mutual inclination is initially small, a large mutual inclination can be generated between $\hat{\boldsymbol{L}}$ and $\hat{\boldsymbol{L}}_{\mathrm{out}}$ because of secular angular-momentum evolution. Subsequently, the LK mechanism can be triggered if the generated mutual inclinations are sufficiently large.

We determine from the numerical integrations the local maximum inclinations, which are shown as a function of $i_{0}$ for $\alpha=5^{\circ}$ in Fig. 5. For $\beta\ll 1$, $i_{\mathrm{max}}\approx i_{0}$; the latter is shown with the black solid lines. For $\beta\gg 1$, $i_{\mathrm{max}}\approx i_{0}+2\alpha$; the latter is shown with the green dotted lines. The large-$\beta$ behaviour can be explained by noting that if $\hat{\boldsymbol{L}}_{\mathrm{out}}$ is precessing rapidly, the dynamics of the inner binary are the same as if $\hat{\boldsymbol{L}}_{\mathrm{out}}$ were aligned along the $\hat{\boldsymbol{z}}$ axis (see Section 3.2). However, as a consequence of the initial geometry, the maximum inclination now occurs when the nodal angle of $\hat{\boldsymbol{L}}$ with respect to $\hat{\boldsymbol{L}}_{\mathrm{out}}$ is $\pi$, such that the angle between $\hat{\boldsymbol{L}}$ and $\hat{\boldsymbol{L}}_{\mathrm{out}}$ is $i_{0}+2\alpha$ .

In the intermediate $\beta$ regime, $i_{\mathrm{max}}$ can be much higher than either $i_{0}$ or $i_{0}+2\alpha$. For $0.6\lesssim\beta\lesssim 2$, $i_{\mathrm{max}}$ is larger than $i_{0}+2\alpha$ for a large range of $i_{0}$. In addition, depending on $\beta$, above a critical angle of $i_{0}$ there is a jump in $i_{\mathrm{max}}$. For example, for $\beta\approx 0.89$, this jump occurs at $i_{0}\approx 10^{\circ}$; for $\beta\approx 1.83$, it occurs at $i_{0}\approx 45^{\circ}$. These jumps are also reflected in the maximum eccentricities (see Fig. 2).

For a range of values $\beta$, orbital flips occur, i.e., $i_{\mathrm{max}}>90^{\circ}$ and the orbital orientation switches from prograde to retrograde. Note that such flips are not possible for $\beta=0$ (we are assuming the quadrupole-order approximation; flips do occur for $\beta=0$ at the octupole order, e.g., Naoz et al. 2011). For $0.21\lesssim\beta\lesssim 0.43$, flips only occur if $i_{0}$ is close to $90^{\circ}$. As expected, these flips are associated with high eccentricities, which can be seen when compared to Fig. 2. For $\beta$ closer to unity, the parameter space for flips is much larger, and a large region of retrograde orbits is populated with $i_{0}$ as low as $\approx 50^{\circ}$.

The maximum inclinations are generally strongly related to the maximum eccentricities. In Fig. 2, we show with solid green lines the ad hoc expression

[TABLE]

where $i_{\mathrm{max}}$ is the largest local maximum inclination determined from the numerical integrations. If $i_{\mathrm{max}}\geq 90^{\circ}$, we set $e_{\mathrm{max}}=1$, which is motivated by the expectation that if retrograde orbits are attained from initially prograde orbits, then $i=90^{\circ}$ at certain times in the evolution; moreover, orbital flips are typically associated with $e\rightarrow 1$.

Except in the limit of large $\beta$, the ad hoc expression equation (24) captures the maximum eccentricities. In particular, the shape of the maximum eccentricity envelopes is captured by equation (24); the same applies for the peaks or spikes, occurring for $0.26\lesssim\beta\lesssim 0.55$. This agreement supports the notion that high mutual inclinations, arising from precession of $\hat{\boldsymbol{L}}_{\mathrm{out}}$, can drive high eccentricities through the LK mechanism. Therefore, in order to understand the dependence of $e_{\mathrm{max}}$ on $i_{0}$, it is useful to understand how $i_{\mathrm{max}}$ depends on $i_{0}$. For this purpose, we discuss below, in Section 4, a simplified model that describes the inclination evolution only. This simplified model is more amenable to analytic treatment, and is therefore useful to gain more insight.

3.4 Short-range forces

We briefly discuss how the above results are modified if short-range forces are included. It is well known that short-range forces due to general relativity and tidal/ rotational bulges of the star/planet tend to suppress eccentricity excitation or limit the maximum eccentricity that can be achieved in LK oscillations (e.g., Wu & Murray 2003; see Fabrycky & Tremaine 2007; Liu et al. 2015 for analytic calculations). Here, we include general relativistic precession at the first post-Newtonian order. The additional precession breaks the scale invariance of the system; therefore, the absolute physical scales need to be specified. The latter are set as follows: the inner and outer semimajor axes and eccentricities are set to 1 and 100 AU, and 0.01 and 0.1, respectively. The assumed masses are $m_{0}=m_{2}=1\,\mathrm{M}_{\odot}$, and $m_{1}=1\,M_{\mathrm{J}}$. With these choices, the characteristic LK time-scale is $t_{\mathrm{LK}}\approx 0.16\,\mathrm{Myr}$, whereas the relativistic precession time-scale is

[TABLE]

where $P_{\mathrm{orb,1}}$ is the inner orbital period, and $r_{\mathrm{g}}=G(m_{0}+m_{1})/c^{2}$ is the gravitational radius. Note that for the assumed masses and semimajor axes, tidal effects play an important role (see Liu et al. 2015; Anderson et al. 2016), so our integrations in this section are for illustrative purposes only.

In Fig. 6, we show the maximum eccentricities as a function of $i_{0}$ for $\alpha=5^{\circ}$ similar to Fig. 2, now with the addition of general relativistic precession. The maximum eccentricities are generally the same compared to the situation without short-range forces, except when $1-e_{\mathrm{max}}$ is less than $\sim 1-10^{-3}$, for which there is a limit on $e_{\mathrm{max}}$. The horizontal green dashed lines show the smallest $1-e_{\mathrm{max}}$ reached for $\beta=0$ if $i_{0}=90^{\circ}$. The latter value can be computed semi-analytically at the quadrupole order (e.g., Blaes et al. 2002). The maximum eccentricity in the case $\beta=0$ is approximately the same compared to the case $\beta\neq 0$. Interestingly, the maximum eccentricities for $\beta\neq 0$ are (marginally) larger compared to the maximum $\beta=0$ value for $\beta\gtrsim 1$.

In summary, although the maximum eccentricities are limited by general relativistic precession, there is still significant enhancement of the eccentricity compared to $\beta=0$ for small inclinations.

4 A simplified model for the mutual inclination

In order to better understand the dynamics underlying the extreme eccentricity excitation observed in Section 3 in the intermediate $\beta$ regime, we explore a simplified model, henceforth referred to as the ‘inclination model’, in which we consider evolution of the angular-momentum vectors only, and disregard changes in the eccentricity vectors. This model is motivated by the observation made in Section 3.3 that the qualitative behaviour of $e_{\mathrm{max}}$ as a function of $i_{0}$ can be explained by determining the maximum inclinations from the numerical integrations, and applying to these the ad hoc expression equation (24). This indicates that if high inclinations can be attained due to interaction with the precessing $\hat{\boldsymbol{L}}_{\mathrm{out}}$, high eccentricities are triggered by the LK mechanism. Of course, our model is only an approximation because we consider the inclination and eccentricity evolution to be decoupled, i.e., the maximum inclination is driven by the interaction with $\hat{\boldsymbol{L}}_{\mathrm{out}}$ and the maximum eccentricity is driven by $i_{\mathrm{max}}$, the latter neglecting the fact that $\hat{\boldsymbol{L}}_{\mathrm{out}}$ is precessing. We note that the inclination model is mathematically very similar to the Hamiltonian model of Lai & Pu (2017).

4.1 Model description

The inclination model is obtained by setting $e=0$ in the model of Section 2, i.e.,

[TABLE]

The Hamiltonian is given by (see equation 13)

[TABLE]

To examine the evolution of $i$, the angle between $\hat{\boldsymbol{L}}$ and $\hat{\boldsymbol{L}}_{\mathrm{out}}$, we set up a rotating frame $(x^{\prime},y^{\prime}z^{\prime})$ such that $\hat{\boldsymbol{z}}^{\prime}=\hat{\boldsymbol{L}}_{\mathrm{out}}$ and $\hat{\boldsymbol{z}}=-\sin(\alpha)\,\hat{\boldsymbol{x}}^{\prime}+\cos(\alpha)\,\hat{\boldsymbol{z}}^{\prime}$. Let $\Omega$ be the nodal angle [measured from the $x^{\prime}$-axis in the $(x^{\prime},y^{\prime}$) plane] of the inner binary, such that $\hat{\boldsymbol{L}}=\cos(\Omega)\sin(i)\,\hat{\boldsymbol{x}}^{\prime}+\sin(\Omega)\sin(i)\,\hat{\boldsymbol{y}}^{\prime}+\cos(i)\,\hat{\boldsymbol{z}}^{\prime}$. Then, equation (27) reduces to

[TABLE]

Since $L\cos(i)$ and $\Omega$ are conjugate canonical variables (note that $L$ is related to $H_{0}$ according to $L=8H_{0}t_{\mathrm{LK}}$), the equations of motion for $(i,\Omega)$ are

[TABLE]

4.2 Maximum inclinations

In Fig. 5, the maximum inclinations obtained by numerically solving the equations of motion (29) are shown with the red dashed lines (below, in Section 4.3, we show how these maximum inclinations can be computed semi-analytically). The inclination model captures several behaviours of $i_{\mathrm{max}}$.

The limit $i_{\mathrm{max}}=i_{0}$ for small $\beta$ and $i_{\mathrm{max}}=i_{0}+2\alpha$ for large $\beta$. 2. 2.

The enhanced $i_{\mathrm{max}}$ relative to $i_{\mathrm{max}}=i_{0}+2\alpha$ in the intermediate $\beta$ regime. 3. 3.

The increase of $i_{\mathrm{max}}$ for decreasing $i_{0}$ for $0.13\lesssim\beta\lesssim 0.55$. 4. 4.

The presence of retrograde orbits for $i_{0}$ close to $90^{\circ}$; these flips are represented in Fig. 5 in the general model of Section 2 by the points near $140^{\circ}$, because the LK mechanism drives the inclination to $\approx 140^{\circ}$, which, in the integrations, is recorded as the maximum inclination.

In Fig. 2, the maximum eccentricities according to the ad hoc equation (24) computed with the maximum inclinations from the inclination model are shown with the red dashed lines (again, we set $e_{\mathrm{max}}=1$ if $i_{\mathrm{max}}>90^{\circ}$). Generally, the inclination model together with equation (24) captures the general trends of $e_{\mathrm{max}}$ with $\beta$. In particular, the envelopes for $0.55\lesssim\beta\lesssim 0.7$ and $2.98\lesssim\beta\lesssim 4.83$ are reproduced. Some of the features of high eccentricity at $0.13\lesssim\beta\lesssim 0.34$ are not consistent with the general model, and the maximum eccentricity envelope is not captured for $1.13\lesssim\beta\lesssim 2.34$. The latter can be expected: the inclination model does not reproduce the additional enhancement in $i_{\mathrm{max}}$ observed in the generalized model, and the $e_{\mathrm{max}}$-envelope appears to be associated with this enhancement. Evidently, the extra enhancement in eccentricity (compare the red dashed and green lines in Fig. 2) is due to a coupling between inclination induced in the inclination model, and the LK mechanism that becomes active at high inclinations.

The inclinations from the inclination model and the implied maximum eccentricities through equation (24) are shown with the red lines in the examples of Fig. 4. For $i_{0}\approx 23^{\circ}$, the inclination as a function of time is reproduced by the inclination model. For $i_{0}\approx 34^{\circ}$, the inclination model still captures some of the features of the inclination. In particular, the maximum inclination matches the value for the general model, and the implied eccentricities are similar.

4.3 Properties of the model

The maximum inclinations can be computed semi-analytically by using the Hamiltonian, equation (28). In Fig. 7, we show phase-space curves in the $(\sin\Omega,i)$ plane for different values of $\beta$ (with $\alpha=5^{\circ}$ fixed). Each panel corresponds to a different $i_{0}$. The following properties are revealed:

The inclination $i$ is always stationary at $\Omega=0$. This is also clear from equation (29a). 2. 2.

For small $\beta$ and $i_{0}$, $\Omega$ circulates, and the maximum inclination is $i_{\mathrm{max}}=i_{0}$, with the other stationary inclination $i_{\mathrm{min}}<i_{0}$. 3. 3.

At a critical value of $\beta$, $\beta_{\mathrm{crit}}$, the curves ‘flip’ over: the minimum inclination is now $i_{0}$, and the maximum inclination is $i_{\mathrm{max}}>i_{0}$ (note that this flip phenomenon is not an orbital flip). The value of $\beta_{\mathrm{crit}}$ depends on $i_{0}$. 4. 4.

For $\beta$ just above $\beta_{\mathrm{crit}}$, $i_{\mathrm{max}}$ increases with increasing $\beta$. As $\beta$ increases further, $i_{\mathrm{max}}$ decreases again.

Since the stationary inclinations (including the maximum inclinations) occur at $\Omega=\pi$, $i_{\mathrm{max}}$ can be found from energy conservation (equation 28),

[TABLE]

where we assumed that, initially, $\Omega=0$ (this is the case in all our numerical integrations).

As can be seen from Fig. 7 (particularly for small $i_{0}$), as $\beta$ increases from zero, at some point $\Omega$ no longer circulates but librates between two critical values of $\Omega$. At the critical $\beta$, the width of libration vanishes, i.e., $\Omega=0$ is constant, and $i=i_{0}$. Using the equations of motion (equation 29b), this implies that $\mathrm{d}\Omega/\mathrm{d}\tau=0$ for $\beta_{\mathrm{crit}}$, giving

[TABLE]

The value of $\beta_{\mathrm{crit}}$ is indicated in the top of each panel of Fig. 7.

In Fig. 8, we plot equation (31) as a function of $i_{0}$, for different values of $\alpha$. We distinguish between prograde ($i_{0}<90^{\circ}$; left-hand panel) and retrograde ($i_{0}>90^{\circ}$; right-hand panel) orbits. For $i_{0}>90^{\circ}$, $\beta_{\mathrm{crit}}$ becomes negative, i.e., the phase-space curve flips no longer occur. Depending on $\alpha$, $\beta_{\mathrm{crit}}$ becomes positive again for sufficiently retrograde orbits, showing that the flip phenomenon returns.

In Fig. 9, we show the maximum inclination as a function of $\beta$. The values of $\beta_{\mathrm{crit}}$ corresponding to each $i_{0}$ (and the chosen $\alpha$) are shown with the vertical dashed lines. This figure illustrates how much the inclination can increase from the initial value $i_{0}$, provided that $\beta$ lies within a certain range. For example, even if only $i_{0}=2^{\circ}$, $i_{\mathrm{max}}$ can reach values of $\approx 70^{\circ}$ if $\beta$ is near 0.5. The minimum value of $\beta$ for which large inclinations are attained is approximated by $\beta_{\mathrm{crit}}$, unless $i_{0}$ is small (i.e., $\lesssim 10^{\circ}$). After rapidly reaching the peak value of $i_{\mathrm{max}}$ as $\beta$ increases beyond $\beta_{\mathrm{crit}}$, $i_{\mathrm{max}}$ steadily decreases (unless the orbit is initially retrograde). In the limit of large $\beta$, the maximum inclination reaches the expected limit $i_{\mathrm{max}}=i_{0}+2\alpha$, indicated in Fig. 9 with the horizontal dashed lines.

In summary, the inclination model shows that high mutual inclinations can be induced, even if $i_{0}$ is small. The degree of inclination enhancement depends sensitively on $\beta$: there is a range in $\beta$ for which $i_{\mathrm{max}}$ is large which depends on $i_{0}$ (and $\alpha$).

5 Discussion

5.1 Applicability of the inclination model

We have shown that the inclination model (Section 4) qualitatively describes how high mutual inclinations can be achieved due to the nodal precession of the outer orbit, even for small initial mutual inclinations. We emphasize that the model, when combined with the ad hoc equation (24), does not quantitatively describe the intermediate-$\beta$ regime; in particular, the envelope of $e_{\mathrm{max}}$ with respect to $i_{0}$ in Fig. 2 is not as ‘deep’ compared to the numerical integrations of the equations of motion. Therefore, the inclination model should not be used for accurate predictions of $e_{\mathrm{max}}$ in the intermediate-$\beta$ regime. However, it can be used to estimate the minimum $i_{\mathrm{max}}$ for a given $i_{0}$ and $\beta$, thereby giving an indication whether or not high eccentricities are to be expected.

5.2 Applications to astrophysical systems

Amongst the astrophysical implications of the enhanced eccentricity oscillations discussed in this paper is the formation of HJs in stellar triple systems. We refer to section 3 of Hamers (2017a) for an example of the ‘2+2’ configuration, in which an HJ is formed through high-eccentricity migration in a stellar triple. In this example, migration is efficient, even though the initial mutual inclination between the proto-HJ and the orbit of the outer stars (i.e., orbit 3) is not large ($\approx 57^{\circ}$), and not high enough to induce migration in the case of a stellar binary (i.e., if orbit 2 were replaced by a point mass). In this example system, $\beta\approx 0.68$ showing that the system is in the regime where high mutual inclinations can be induced, driving high-eccentricity secular oscillations and tidal migration.

For reference, we give the explicit expressions for $\beta$ in terms of physical parameters111Equation (32) is the same as equations 10 of Hamers & Portegies Zwart (2016), and 1 of Hamers (2017a), modulo a factor of $(3/4)\cos(i_{23,\mathrm{init}})$. Equation (33) is the same as equations 13 of Hamers et al. (2015), 9 of Hamers & Portegies Zwart (2016), and 3 of Hamers (2017a), again modulo a factor of $(3/4)\cos(i_{23,\mathrm{init}})$. . For ‘2+2’ systems,

[TABLE]

Note that $a_{3}$ does not enter in this expression; it drops out in the ratio of the LK time-scales of orbit pairs (1,3) and (2,3). Of course, $a_{3}$ does set the absolute time-scale on which the secular evolution in the system occurs. For ‘3+1’ systems,

[TABLE]

Note that in this paper, we assumed $e_{2}=0$.

6 Conclusions

We have studied in detail the secular dynamics of hierarchical quadruple systems in either the ‘2+2’ or ‘3+1’ configurations (Fig. 1). The ‘2+2’ configuration applies, for example, to a planet orbiting a star, which is orbited by a more distant stellar binary (orbit 2), in a relatively wide orbit (orbit 3). An example of the ‘3+1’ configuration is a planet orbiting a star, which is orbited by two more distant objects (one more distant than the other, in hierarchical orbits); we assume that the eccentricity of the orbit of the first distant object is zero throughout the evolution (i.e., the outer objects are inclined by no more than $\approx 40^{\circ}$ such that the LK mechanism does not operate). In particular, the ‘2+2’ configuration applies to HJs in stellar triples (Hamers, 2017a), and the ‘3+1’ configuration applies to hypothetical planetary companions to HJs in stellar binaries (Hamers, 2017b). We have formulated the secular equations of motion for both configurations in one generalized model. In this model, the problem has been reduced to a hierarchical three-body problem with the perturbed outer orbital angular-momentum axis $\hat{\boldsymbol{L}}_{\mathrm{out}}$ precessing around a fixed axis, $\hat{\boldsymbol{z}}$, at a constant and prescribed (i.e., known) rate $\Omega_{\mathrm{out}}$ and an angle $\alpha$. Our conclusions are as follows.

1. Extremely high eccentricities, $e_{\mathrm{max}}\rightarrow 1$, can be attained in the inner orbit (i.e., the inner planetary orbit), even if the initial inclination of the inner orbit with respect to the outer orbit is small and less than the minimum LK angle (i.e., $<40^{\circ}$). These enhancements occur already at the quadrupole order.

2. The nature of the eccentricity enhancement depends sensitively on the (dimensionless) quantity $\beta\equiv\Omega_{\mathrm{out}}t_{\mathrm{LK}}$, where $t_{\mathrm{LK}}$ is the LK time-scale of the inner orbit driven by the outer orbit.

The case $\beta=0$ corresponds to the unperturbed three-body problem, for which the maximum eccentricity is given by the canonical LK expression, equation (15). High eccentricities can be attained, but only if $i_{0}$ is high.

In the limit of $\beta\gg 1$, we obtained an analytic result for $e_{\mathrm{max}}$ by averaging over the precessional cycle of the outer orbit (equation 23). Effectively, the classical LK result for the maximum eccentricity can be used, replacing the initial inclination between the inner and outer orbits by the initial inclination of the inner orbit with respect to the $\hat{\boldsymbol{z}}$ axis, i.e., the axis around which $\hat{\boldsymbol{L}}_{\mathrm{out}}$ is precessing.

Most importantly, if $\beta\sim 1$, then large $e_{\mathrm{max}}$ can be achieved for modest initial inclinations. In particular, $e_{\mathrm{max}}\rightarrow 1$ for $i_{0}\gtrsim 50^{\circ}$ if $\alpha=5^{\circ}$ (see Fig. 2). For $\alpha=5^{\circ}$, there is a complicated dependence of $e_{\mathrm{max}}$ on $i_{0}$, with chaotic ‘ridges’ occurring at specific $i_{0}$ and $0.26\lesssim\beta\lesssim 0.7$, and a broad envelope of maximum eccentricities appearing for $0.7\lesssim\beta\lesssim 3$. For larger $\alpha$, there are fewer ridges, and a large chaotic region appears for $\beta$ around 1 (see Fig. 10, in which $\alpha=30^{\circ}$).

3. The eccentricity excitations around $\beta\sim 1$ can be explained by the enhanced maximum inclination $i_{\mathrm{max}}$ as a result of the precessing $\hat{\boldsymbol{L}}_{\mathrm{out}}$, which subsequently drives strong LK evolution. This was illustrated in Fig. 2, in which the maximum eccentricities were computed numerically from the maximum inclinations using the ad hoc expression equation (24).

4. We briefly considered how the above is modified with the addition of short-range forces, in particular, general relativistic precession (Section 3.4). Although the maximum eccentricities are limited by general relativistic precession, there is still significant enhancement of the eccentricity compared to $\beta=0$ for small inclinations.

5. To further explore point (3), we considered a simplified model, the ‘inclination model’ (Section 4), in which we set $e=0$ for the inner orbit. This simple model can be studied semi-analytically, and correctly predicts that high mutual inclinations are attained even if $i_{0}$ is small, provided that $\beta$ is around unity. For example, even if $i_{0}$ is as small as $2^{\circ}$, $i_{\mathrm{max}}$ can reach values of $\approx 70^{\circ}$ if $\beta$ is near 0.5 and $\alpha=5^{\circ}$. For a large range of inclinations, a good estimate for the minimum value of $\beta$ for which high mutual inclinations can be achieved is given by equation (31). Typically, the inclination enhancement drops off rapidly as $\beta$ increases beyond $\beta_{\mathrm{crit}}$, with little enhancement occurring for $\beta\gtrsim 2$ (see Fig. 9). The behaviour of the eccentricity excitation around $\beta\sim 1$ for the general model can be crudely reproduced if the analytic ad hoc LK expression is used to compute $e_{\mathrm{max}}$ from the maximum inclination obtained in the semi-analytic inclination model.

Acknowledgements

We thank the referee, Cristobal Petrovich, for insightful and helpful comments. We also thank Scott Tremaine for comments on the manuscript, and Fabio Antonini, Evgeni Grishin and Ben Bar-Or for simulating discussions. ASH gratefully acknowledges support from the Institute for Advanced Study, and NASA grant NNX14AM24G. DL has been supported in part by NASA grants NNX14AG94G and NNX14AP31G, and a Simons Fellowship from the Simons Foundation.

Appendix A Maximum eccentricity as a function of the initial inclination for $\alpha=30^{\circ}$

In Fig. 10, we show a figure similar to Fig. 2, now with $\alpha=30^{\circ}$. The region in the parameter space for which large eccentricities are reached is now larger: $e_{\mathrm{max}}\rightarrow 1$ even for small $i_{0}$, in the range $0.43\lesssim\beta\lesssim 1.44$. The inclination model predicts $i_{\mathrm{max}}\geq 90^{\circ}$ in this regime, hence we set $e_{\mathrm{max}}=1$ (see the discussion in Section 3.3). Note that in the limit of large $\beta$, the inclination model correctly predicts $i_{\mathrm{max}}$ (compare the solid green and red dashed lines in Fig. 10); however, in this limit, the maximum eccentricity should be computed from $i_{z,0}$ (as shown with the solid black lines in Fig. 10; see Section 3.2 and equation 23).