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

Adrian S. Hamers and Dong Lai
Institute for Advanced Study, School of Natural Sciences, Einstein Drive, Princeton, NJ 08540, USA
Cornell Center for Astrophysics and Planetary Science, Department of Astronomy, Cornell University, Ithaca, NY 14853, USA
Institute for Theory and Computation, Harvard-Smithsonian Center for Astrophysics, Cambridge, MA 02138, USA
E-mail: hamers@ias.edu
Accepted 2017 May 23. Received 2017 May 20; in original form 2017 April 26
###### 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.

###### keywords:
planets and satellites: dynamical evolution and stability – planet-star interactions – gravitation
pubyear: 2017pagerange: Secular chaotic dynamics in hierarchical quadruple systems, with applications to hot Jupiters in stellar binaries and triplesA

## 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, , is negligible compared to the angular momentum of the outer orbit, . In this limit, the outer orbital angular momentum vector is constant (therefore, is constant as well), and the system is completely integrable. This system gives rise to well-known LK oscillations, occurring if and are initially inclined by more than , and with a maximum eccentricity of , where 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, (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.

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 and angular-momentum 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

 dj1dt =34t−1LK,2+2,13[(j1⋅^L3)(j1×^L3)−5(e1⋅^L3)(e1×^L3)]; (1a) de1dt =34t−1LK,2+2,13[(j1⋅^L3)(e1×^L3)+2(j1×e1)−5(e1⋅^L3)(j1×^L3)]; (1b) dj2dt =34t−1LK,2+2,23[(j2⋅^L3)(j2×^L3)−5(e2⋅^L3)(e2×^L3)]; (1c) de2dt =34t−1LK,2+2,23[(j2⋅^L3)(e2×^L3)+2(j2×e2)−5(e2⋅^L3)(j2×^L3)]; (1d) dj3dt =34t−1LK,2+2,13Λ1Λ3[5(e1⋅^L3)(e1×^L3)−(j1⋅^L3)(j1×^L3)]+34t−1LK,2+2,23Λ2Λ3[5(e2⋅^L3)(e2×^L3)−(j2⋅^L3)(j2×^L3)]; (1e) de3dt =38t−1LK,2+2,13Λ1Λ31√1−e23[{(1−6e21)+25(e1⋅^L3)2−5(j1⋅^L3)2}(e3×^L3)−10(e1⋅^L3)(e3×e1)+2(j1⋅^L3)(e3×j1)] +38t−1LK,2+2,23Λ2Λ31√1−e23[{(1−6e22)+25(e2⋅^L3)2−5(j2⋅^L3)2}(e3×^L3)−10(e2⋅^L3)(e3×e2)+2(j2⋅^L3)(e3×j2)] ≡e3×f(e1,j1,e2,j2,j3). (1f)

Here, is the angular momentum of orbit for a circular orbit which is constant in the secular approximation, i.e., , where and are the reduced and the total mass, respectively, of binary . The (two) LK time-scales are given by

 tLK,2+2,13 =m0+m1m2+m3 ⎷a31G(m0+m1)(a3a1)3(1−e23)3/2;tLK,2+2,23=m2+m3m0+m1 ⎷a32G(m2+m3)(a3a2)3(1−e23)3/2. (2a)

The function in equation (1) is independent of , showing that precesses around and the magnitude of 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, , and setting . Evidently, the torque of orbit 3 can excite if orbits 2 and 3 are initially sufficiently inclined. Here, we assume that , such that the LK mechanism is not active for the orbit pair (2,3); therefore, is constant. The assumption on the implies that the first term in equation (1e), proportional to , is negligible compared to the second term, which is proportional to (also assuming that and are not too distinct). Equation (1e) then reads

 dj3dt ≈−34t−1LK,2+2,23Λ2Λ3(^L2⋅^L3)(^L2×^L3). (3)

Using that the total angular momentum vector, , is conserved, equation (3) can be written in the form

 d^L3dt ≈−Ω3×^L3, (4)

where is a constant vector with magnitude

 Ω3 =34t−1LK,2+2,23cos(i23,init), (5)

with as the initial inclination between and .

In summary, the secular dynamics of the restricted problem, and , 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).

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

 dj1dt =34t−1LK,3+1,12[(j1⋅^L2)(j1×^L2)−5(e1⋅^L2)(e1×^L2)]+34t−1LK,3+1,13[(j1⋅^L3)(j1×^L3)−5(e1⋅^L3)(e1×^L3)]; (6a) de1dt =34t−1LK,3+1,12[(j1⋅^L2)(e1×^L2)+2(j1×e1)−5(e1⋅^L2)(j1×^L2)]+34t−1LK,3+1,13[(j1⋅^L3)(e1×^L3)+2(j1×e1) −5(e1⋅^L3)(j1×^L3)]; (6b) dj2dt =34t−1LK,3+1,12Λ1Λ2[5(e1⋅^L2)(e1×^L2)−(j1⋅^L2)(j1×^L2)]−34t−1LK,3+1,23[5(e2⋅^L3)(e2×^L3)−(j2⋅^L3)(j2×^L3)]; (6c) de2dt =38t−1LK,3+1,12Λ1Λ21√1−e22[{(1−6e21)+25(e1⋅^L2)2−5(j1⋅^L2)2}(e2×^L2)−10(e1⋅^L2)(e2×e1)+2(j1⋅^L2)(e2×j1)] +34t−1LK,23[(j2⋅^L3)(e2×^L3)+2(j2×e2)−5(e2⋅^L3)(j2×^L3)]; (6d) dj3dt =34t−1LK,23Λ2Λ3[5(e2⋅^L3)(e2×^L3)−(j2⋅^L3)(j2×^L3)]+34t−1LK,13Λ1Λ3[5(e1⋅^L3)(e1×^L3)−(j1⋅^L3)(j1×^L3)]; (6e) de3dt =38t−1LK,3+1,23Λ2Λ31√1−e23[{(1−6e22)+25(e2⋅^L3)2−5(j2⋅^L3)2}(e3×^L3)−10(e2⋅^L3)(e3×e2)+2(j2⋅^L3)(e3×j2)] (6f) +38t−1LK,3+1,13Λ1Λ31√1−e23[{(1−6e21)+25(e1⋅^L3)2−5(j1⋅^L3)2}(e3×^L3)−10(e1⋅^L3)(e3×e1)+2(j1⋅^L3)(e3×j1)] (6g) ≡e3×f(e1,j1,e2,j2,j3). (6h)

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

 tLK,3+1,12 =m0+m1m2 ⎷a31G(m0+m1)(a2a1)3(1−e22)3/2;tLK,3+1,23=m0+m1+m2m3 ⎷a32G(m0+m1+m2)(a3a2)3(1−e23)3/2; (7a) tLK,3+1,13 =m0+m1m3 ⎷a31G(m0+m1)(a3a1)3(1−e23)3/2. (7b)

Again, the function in equation (6f) is independent of , showing that precesses around and is constant. In this case, the LK time-scale is generally not constant because is not generally constant.

We make the assumption . Furthermore, we assume . As before, we assume that such that the LK mechanism does not operate for the orbit pair (2,3). Equation (6c) can then be written in the form

 d^L2dt ≈−Ω2×^L2, (8)

where is constant, and

 Ω2 =34t−1LK,3+1,23cos(i23,init). (9)

In addition, for hierarchical systems () with a sufficiently large ratio , , 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 , 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 . The model can be written in the scaled form

 djdτ =34[(j⋅^Lout)(j×^Lout)−5(e⋅^Lout)(e×^Lout)]; (10a) dedτ =34[(j⋅^Lout)(e×^Lout)+2(j×e)−5(e⋅^Lout)(j×^Lout)]; (10b) d^Loutdτ =−β(^z×^Lout). (10c)

Here, and correspond to the inner test-particle orbit, and with the LK time-scale given by

 tLK=m0+m1mout√a3G(m0+m1)(aouta)3(1−e2out)3/2. (11)

The semimajor axis of the test-particle orbit is denoted with , and the subscript ‘out’ denotes properties of the ‘outer’ companion, which depends on the hierarchical configuration (in the ‘2+2’ case, , and , whereas , and in the ‘3+1’ case). The dimensionless and constant parameter is defined as

 β≡ΩouttLK. (12)

In equation (10c), the fixed axis around which precesses is taken to be the axis. We denote the angle between and with .

#### 2.3.2 Hamiltonian and the ‘LK constant’

If in equations (10), then the (conserved) Hamiltonian is mathematically equivalent to the three-body test-particle Hamiltonian. For nonzero , this Hamiltonian is not constant because of the fixed precession of . However, one can transform to a frame that is rotating around the axis with the same rate at which is precessing around the axis. This gives (e.g., Kinoshita 1993)

 Hgen,rot =H0[(1−6e2)+15(e⋅^Lout)2−3(j⋅^Lout)2+8β(j⋅^z)], (13)

where is a constant.

In the test-particle three-body limit (), the quantity , where is the mutual inclination between the inner and outer orbits, is a constant of the motion (known as the ‘LK’ or ‘Kozai’ constant). For , the ‘LK constant’ is no longer conserved, but evolves according to

 ddt[√1−e2cos(i)]=−Ωoutj⋅(^z×^Lout). (14)

## 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 , and are initially all assumed to lie in the same plane; the latter plane is the -plane. In addition, , the angle about which precesses around , is defined such that positive corresponds to having a positive component along the axis. The initial mutual inclination between and is denoted with . The initial mutual inclination between and is . The eccentricity vector is assumed to be initially parallel to the axis, and its initial magnitude is set to .

We carried out sets of 1600 integrations on a linear grid with running from 0 to 10, and the initial , , running from to . Here, we set or for each set. The duration of each integration was , i.e., corresponding to a physical time-span of , and approximately 1000 LK oscillations if were fixed. During the integrations, local maxima of were recorded, and the latter are shown as a function of in Fig. 2 for (results for are given in Appendix A). Each of the 20 panels corresponds to a different value of indicated therein; the number of bins in is 80.

Below, we discuss the maximum eccentricity behaviour as a function of for three regimes of : small (Section 3.1), large (Section 3.2) and intermediate (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 .

### 3.1 Small β

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

 emax=√1−53cos2(i0), (15)

which applies to the ‘unperturbed’ problem with , and assuming initially and . This result, which is described here in detail for further reference, can be obtained by dotting equation (10b) with , giving the stationary points, i.e.,

 0=dedτ=−154(e⋅^Lout)[^e⋅(j×^Lout)]. (16)

Assuming that at the stationary point , and , this can be rewritten as the condition

 (^e⋅^Lout)2=1−(^L⋅^Lout)2, (17)

where we used a vector identity for the scalar product of two triple products. Using equation (17) to eliminate , and using the LK constant to express in terms of the initial inclination, conservation of the Hamiltonian yields an algebraic equation for with the solution equation (15) for the (maximum) stationary point.

As shown in Fig. 2, for the smallest non-zero values of , equation (15) still gives a good description as expected. For but still , noticeable deviations start to occur; in particular, large eccentricity excitations, much larger than based on equation (15), are possible for inclinations near .

### 3.2 Large β

If is large (), then is precessing rapidly around the axis compared to the time-scale at which evolves (i.e., on a time-scale on the order of ). Intuitively, one might expect in this case that the rapid precession of effectively means that is pointing along , although with a modified ‘effective length’ depending on , the angle at which precesses around .

In the limit, we can average the Hamiltonian equation (13) over a precessional cycle of . To achieve this, we write

 ^Lout=cos(ϕ)sin(α)^x+sin(ϕ)sin(α)^y+cos(α)^z, (18)

where is the (rapidly changing) precessional phase angle. The precession average of a quantity is defined by

 ⟨(...)⟩≡12π∫2π0dϕ(...). (19)

A number of useful averages are

 ⟨[^e⋅(j×^Lout)]⟩ =[^e⋅(j×^z)]cos(α); (20a) ⟨(^e⋅^Lout)2⟩ =14[1+(^e⋅^z)2−cos(2α){1−3(^e⋅^z)2}]; (20b) ⟨(^L⋅^Lout)2⟩ =14[1+(^L⋅^z)2−cos(2α){1−3(^L⋅^z)2}]. (20c)

Averaging equation (14) gives

 ddt(j⋅^z)≈0. (21)

In other words, in the limit , the LK constant, , is replaced by . 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 . From equation (20a), it follows that , which implies .

Next, we substitute these results into the Hamiltonian equation (13). By equation (21), the last term in the Hamiltonian is constant and can therefore be omitted. Dropping the latter term, using to eliminate and equation (21) to relate to the initial value, , and applying equations (20b) and (20c), the Hamiltonian at the stationary eccentricity in the limit reads

 ¯¯¯¯¯Hgen,rot,β≫1,stat −34(1−e2){1+cos2(iz,0)1−e2−cos(2α)(1−3cos2(iz,0)1−e2)}]. (22)

Equating (3.2) to the initial averaged Hamiltonian with ,

 ¯¯¯¯¯Hgen,rot,β≫1,init=H0[1−34{1+cos2(iz,0)−cos(2α)(1−3cos2(iz,0))}],

we can solve analytically for the maximum eccentricity as a function of and . The solution is

 emax =√1−53cos2(iz,0), (23)

which is, remarkably, independent of (if expressed in terms of ). Equation (23) is simply equation (15) with replaced by . 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 and the axis.

In Fig. 2, we plot equation (23) with the solid black lines. The numerical integrations indeed agree well with equation (23) for large (). Note that the largest does not occur at , but at , since .

As an additional test of equation (23), we show in Fig. 3 with blue points the maximum eccentricities as a function of determined by numerically integrating the equations of motion, where we set and . The numerical results are in good agreement with equation (23), shown with the solid red line.

### 3.3 Intermediate β

In the intermediate- regime (), as a function of is much more complex. In particular, for around unity, extremely large eccentricity excitations are possible. Whereas in the limit , only for , if , for a much larger range of with in the integrations reaching values lower than . The precise range of for extreme eccentricity excitation depends on (and ; compare, Figures 2 and 10). For lower , , and , there are specific inclinations for which , indicative of an overlap of resonances, and, therefore, chaos (Chirikov, 1979). For larger , , there is a broad range of for which extreme eccentricities are reached. This range of local maximum eccentricities is manifested as a band in the plane. The upper boundary on of this band appears to be smooth function of . 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 ) for a given .

In the chaotic regime of intermediate , 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 and are given in Fig. 4. In particular, for , there is significant eccentricity excitation even though is less than the critical LK angle. For , 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 , we explore the connection between the maximum inclination, , 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 and 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 for in Fig. 5. For , ; the latter is shown with the black solid lines. For , ; the latter is shown with the green dotted lines. The large- behaviour can be explained by noting that if is precessing rapidly, the dynamics of the inner binary are the same as if were aligned along the axis (see Section 3.2). However, as a consequence of the initial geometry, the maximum inclination now occurs when the nodal angle of with respect to is , such that the angle between and is .

In the intermediate regime, can be much higher than either or . For , is larger than for a large range of . In addition, depending on , above a critical angle of there is a jump in . For example, for , this jump occurs at ; for , it occurs at . These jumps are also reflected in the maximum eccentricities (see Fig. 2).

For a range of values , orbital flips occur, i.e., and the orbital orientation switches from prograde to retrograde. Note that such flips are not possible for (we are assuming the quadrupole-order approximation; flips do occur for at the octupole order, e.g., Naoz et al. 2011). For , flips only occur if is close to . As expected, these flips are associated with high eccentricities, which can be seen when compared to Fig. 2. For closer to unity, the parameter space for flips is much larger, and a large region of retrograde orbits is populated with as low as .

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

 emax=√1−53cos2(imax), (24)

where is the largest local maximum inclination determined from the numerical integrations. If , we set , which is motivated by the expectation that if retrograde orbits are attained from initially prograde orbits, then at certain times in the evolution; moreover, orbital flips are typically associated with .

Except in the limit of large , 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 . This agreement supports the notion that high mutual inclinations, arising from precession of , can drive high eccentricities through the LK mechanism. Therefore, in order to understand the dependence of on , it is useful to understand how depends on . 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