format: text

text:
# Secular Evolution of Hierarchical Planetary Systems

## 1 INTRODUCTION

Extrasolar planet searches using high-precision radial velocity observations to measure the reflex motion of the host stars along the line of sight have now yielded over $`100`$ extrasolar planets.111 See, e.g., http://www.obspm.fr/planets for a continuously updated catalog.  The discoveries include $`10`$ multiple planet systems with either two or, in the cases of $`\upsilon `$ And and 55 Cnc, three detected planets. There are also indications that about half of the stars with one known planet are likely to have additional detectable distant companions (Fischer et al. 2001). Many of the known multiple planet systems exhibit interesting dynamics. The two planets about GJ 876 are deep in three orbital resonances at the 2:1 mean-motion commensurability (Laughlin & Chambers 2001; Lee & Peale 2002), and there are possibly two other systems with mean-motion resonances: 2:1 for the two planets about HD 82943222 See http://obswww.unige.ch/$∼$udry/planet/hd82943syst.html.  (Goździewski & Maciejewski 2001) and 3:1 for the inner two planets about 55 Cnc (Marcy et al. 2002; Lee & Peale 2003). The outer two planets of the $`\upsilon `$ And system are apparently locked in a secular resonance with the lines of apsides of the two orbits being aligned on average (Rivera & Lissauer 2000; Lissauer & Rivera 2001; Chiang, Tabachnik, & Tremaine 2001), and as we shall demonstrate in this paper, the two planets of the HD 12661 system are also locked in a secular resonance, but with the lines of apsides being anti-aligned on average. The resonances are often vital in ensuring the long term stability of the multiple planet systems.

Although the origin, evolution, and stability of the orbital configurations of multiple planet systems can usually be analyzed with direct numerical integrations of the full equations of motion, a theoretical understanding of the dynamics often requires the development and application of theories with analytic approximations. The approximate theories also allow one to explore the parameter space much more rapidly than direct numerical integrations in the regions where the approximations are valid.

Because the orbital eccentricities and inclinations of the planets in the solar system are generally small while the ratios of the orbital semimajor axes of adjacent planets are generally large ($`a_j/a_{j+1}0.5`$ except for the Mars-Jupiter pair), the classical perturbation theory developed for the solar system is based on an expansion of the disturbing functions that describe the mutual gravitational interactions of the planets in powers of the eccentricities and inclinations. In particular, the Laplace-Lagrange secular solution for the evolution of a two-planet system that is not affected by mean-motion commensurabilities is based on retaining just the secular terms (i.e., those not involving the mean longitudes) in the disturbing functions up to second order in the eccentricities and inclinations (see, e.g., Murray & Dermott 1999). Since high orbital eccentricities are common among extrasolar planets and the distribution of the orbital eccentricities of the extrasolar planets in the known multiple planet systems is similar to that of the single planets (Fischer et al. 2003), the classical Laplace-Lagrange secular perturbation theory is not in general adequate for describing the secular evolution of the extrasolar planetary systems that are not affected by mean-motion commensurabilities.

On the other hand, many of the known multiple planet systems are hierarchical in the sense that the ratio(s) of the semimajor axes ($`a_1/a_2`$ for the two-planet systems and $`a_1/a_2`$ and/or $`a_2/a_3`$ for the three-planet systems) are small. There are $`12`$ pairs of adjacent planets in the $`10`$ known multiple planet systems (two of which have three planets). Among the $`9`$ pairs that are not known or suspected to be in mean-motion resonances, five have $`a_j/a_{j+1}0.1`$ and all but one (47 UMa) have $`a_j/a_{j+1}1/3`$. Thus a secular perturbation theory for a two-planet system that is based on an expansion in $`a_1/a_2`$ and valid for high eccentricities should provide an accurate description of the secular evolution of many extrasolar planetary systems.

It is useful to consider hierarchical two-planet systems in the context of the general hierarchical triple systems in which a third body orbits an inner binary on a much wider orbit. A hierarchical triple system can be treated as two binaries on slowly perturbed Kepler orbits by using Jacobi coordinates, where the position of the secondary of mass $`m_1`$ of the inner binary is relative to the primary of mass $`m_0`$ and the position of the third body of mass $`m_2`$ is relative to the center of mass of $`m_0`$ and $`m_1`$. A hierarchical two-planet system is simply a hierarchical triple system with $`m_1`$ and $`m_2`$ much smaller than $`m_0`$. In § 2 we derive the orbital parameters in Jacobi coordinates obtained by the observers from a two (or more generally multiple) Kepler fit to the radial velocity variations of a host star and show that Jacobi coordinates should be used in any analyses of hierarchical (and possibly other types of) planetary systems. Star-centered or astrocentric coordinates can introduce significant high-frequency variations in orbital elements that should be nearly constant on orbital timescales. The high-frequency variations can then lead to erroneous sensitivity of the evolution of the orbital elements to the starting epoch.

Secular perturbation theories based on an expansion in $`\alpha =a_1/a_2`$ have been developed for hierarchical triple systems. The expansion to order $`\alpha ^2`$, called the quadrupole approximation, was developed by Kozai (1962) and Harrington (1968). This quadrupole-level secular perturbation theory is successful in explaining the Kozai mechanism, whereby the perturbations between the inner binary and the third body can lead to large variations in the eccentricity of the inner binary and the mutual inclination angle between the inner and outer binaries if the initial mutual inclination is sufficiently high. However, the quadrupole approximation is not adequate for studying hierarchical two-planet systems. Although there is no direct information on the mutual inclination angles in the known hierarchical two-planet systems, the formation of planets from a common disk of materials surrounding the host star makes nearly coplanar orbits the most probable configuration. When the orbits are coplanar, the conservation of total angular momentum of the system and the secularly constant semimajor axes mean that the eccentricities of the two orbits are coupled and oscillate out of phase. The quadrupole term does not contribute to these eccentricity oscillations for coplanar orbits. Marchal (1990), Krymolowski & Mazeh (1999), and Ford, Kozinsky, & Rasio (2000) have extended the approximation to octupole (i.e., $`\alpha ^3`$) order. Blaes, Lee, & Socrates (2002) have recently applied this octupole-level secular perturbation theory, with a sign error corrected (see footnote 5 below) and with modifications to include the effects of general relativistic precession and gravitational radiation on the inner binary, to study the dynamical evolution of hierarchical triples of supermassive black holes. Unlike the quadrupole term, the octupole term does produce eccentricity oscillations for coplanar orbits.

In this paper we use both the octupole-level secular perturbation theory and direct numerical orbit integrations to investigate the dynamical evolution of coplanar hierarchical two-planet systems. The applicability of the octupole theory is limited by its use of an expansion in $`\alpha `$ and an averaging over the inner and outer orbital motions. In particular, the orbit averaging eliminates the effects of mean-motion commensurabilities and the possible development of instabilities. We establish the validity and limits of the octupole theory by comparison with direct numerical orbit integrations. In § 3 we summarize the derivation of the octupole theory, compare the octupole theory to the classical Laplace-Lagrange secular solution, and deduce some useful results from the octupole equations analytically. In § 4 we study the dynamical evolution of the HD 168443 and HD 12661 systems and their variants. The HD 168443 system is not in a secular resonance and its secular resonance variable $`\varpi _1\varpi _2`$ circulates, where $`\varpi _{1,2}`$ are the longitudes of periapse of the inner and outer orbits, respectively. For the HD 168443 system and a wide variety of systems with $`\alpha 0.1`$ (including some for which the octupole theory predicts rather unusual dynamical behaviors), we show that the octupole results are in excellent agreement with the direct integration results. Direct integrations of two-planet systems similar to HD 168443, but with different initial $`a_2`$ (and hence different initial $`\alpha `$), are used to show that systems with initial $`\alpha `$ above a critical value are generally unstable. As anticipated by the analytic results derived in § 3, the HD 12661 system is in a secular resonance with $`\varpi _1\varpi _2`$ librating about $`180^{}`$, and it shows large amplitude variations of both eccentricities. We show that the evolution of the HD 12661 system with the best-fit orbital parameters and $`\mathrm{sin}i=1`$ ($`i`$ is the inclination of the orbital plane from the plane of the sky) is affected by the proximity to the 11:2 mean-motion commensurability, but that small changes in the orbital period of the outer planet within the uncertainty can result in configurations that are not affected by mean-motion commensurabilities. For the latter type of configurations, we show that the octupole results are in reasonably good agreement with the direct integration results, even though $`\alpha `$($`0.32`$) is quite large. We also consider the effects of varying the inclination $`i`$ and show that the HD 12661 system is unstable if $`\mathrm{sin}i0.3`$. Our conclusions are summarized in § 5.

## 2 JACOBI ORBITAL PARAMETERS FROM MULTIPLE KEPLER FITS TO RADIAL VELOCITY OBSERVATIONS

Except for systems such as GJ 876, where the perturbations between the two planets are significant on orbital timescales and a dynamical fit to the stellar radial velocity variations is essential (Laughlin & Chambers 2001; Rivera & Lissauer 2001; Nauenberg 2002), it is often adequate to fit the radial velocity variations of a star with two (or more) planets over the time span of the available observations by assuming that the planets are on unperturbed Kepler orbits. Many authors have assumed that the orbital parameters obtained from the multiple Kepler fits are in astrocentric coordinates, but Lissauer & Rivera (2001; see also Note Added in Proof of Rivera & Lissauer 2000) have pointed out that Jacobi coordinates “better emulate” the assumption of unperturbed Kepler orbits. In this section we derive explicitly the orbital parameters in Jacobi coordinates obtained from a two-Kepler fit. It is straightforward to generalize the derivation to an $`N`$-Kepler fit. We show that especially for hierarchical systems such as HD 168443, the use of astrocentric coordinates can introduce erroneous features in the evolution of the orbital elements.

Let us consider a system consisting of a central star of mass $`m_0`$, an inner planet of mass $`m_1`$, and an outer planet of mass $`m_2`$, and use Jacobi coordinates, with $`𝒓_1`$ being the position of $`m_1`$ relative to $`m_0`$ and $`𝒓_2`$ being the position of $`m_2`$ relative to the center of mass of $`m_0`$ and $`m_1`$. We shall refer to the orbit of $`m_1`$ relative to $`m_0`$ as the inner orbit and the orbit of $`m_2`$ relative to the center of mass of $`m_0`$ and $`m_1`$ as the outer orbit.

If the inner orbit is an unperturbed Kepler orbit, the line-of-sight (LOS) component of the velocity of $`m_1`$ relative to $`m_0`$ is

$`V_{1,r}`$ $`=`$ $`\left[\dot{r}_1\mathrm{sin}(\omega _1+f_1)+r_1\dot{f}_1\mathrm{cos}(\omega _1+f_1)\right]\mathrm{sin}i_1`$ (1)
$`=`$ $`{\displaystyle \frac{2\pi a_1}{P_1\sqrt{1e_1^2}}}\left[\mathrm{cos}(\omega _1+f_1)+e_1\mathrm{cos}\omega _1\right]\mathrm{sin}i_1,`$
where a dot over a symbol denotes $`d/dt`$ and $`a_1`$, $`e_1`$, $`i_1`$, $`\omega _1`$, $`f_1`$, and

$$P_1=\frac{2\pi a_1^{3/2}}{\sqrt{G(m_0+m_1)}}$$
(2)
are, respectively, the semimajor axis, eccentricity, inclination, argument of periapse, true anomaly, and period of the inner orbit. In equation (2) $`G`$ is the gravitational constant. Note that the reference plane is the plane of the sky, and the planet approaches the observer at the ascending node, but the radial velocity of an approaching planet is negative. Then the LOS component of the velocity of $`m_0`$ relative to the center of mass of $`m_0`$ and $`m_1`$ is

$$V_{0,r}^{}=\frac{m_1}{m_0+m_1}V_{1,r}=K_1\left[\mathrm{cos}(\omega _1+f_1)+e_1\mathrm{cos}\omega _1\right],$$
(3)
where the amplitude

$$K_1=\left(\frac{2\pi G}{P_1}\right)^{1/3}\frac{m_1\mathrm{sin}i_1}{(m_0+m_1)^{2/3}}\frac{1}{\sqrt{1e_1^2}}.$$
(4)
Similarly, if the outer orbit is an unperturbed Kepler orbit, the LOS component of the velocity of $`m_2`$ relative to the center of mass of $`m_0`$ and $`m_1`$ is

$`V_{2,r}={\displaystyle \frac{2\pi a_2}{P_2\sqrt{1e_2^2}}}\left[\mathrm{cos}(\omega _2+f_2)+e_2\mathrm{cos}\omega _2\right]\mathrm{sin}i_2,`$ (5)
where $`a_2`$, $`e_2`$, $`i_2`$, $`\omega _2`$, $`f_2`$, and

$$P_2=\frac{2\pi a_2^{3/2}}{\sqrt{G(m_0+m_1+m_2)}}$$
(6)
are, respectively, the semimajor axis, eccentricity, inclination, argument of periapse, true anomaly, and period of the outer orbit. Then the LOS component of the velocity of the center of mass of $`m_0`$ and $`m_1`$ relative to the center of mass of the whole system is

$$V_{01,r}^{}=\frac{m_2}{m_0+m_1+m_2}V_{2,r}=K_2\left[\mathrm{cos}(\omega _2+f_2)+e_2\mathrm{cos}\omega _2\right],$$
(7)
where the amplitude

$$K_2=\left(\frac{2\pi G}{P_2}\right)^{1/3}\frac{m_2\mathrm{sin}i_2}{(m_0+m_1+m_2)^{2/3}}\frac{1}{\sqrt{1e_2^2}}.$$
(8)
Thus, if the orbits of the planets in Jacobi coordinates are unperturbed Kepler orbits, the radial velocity of the star $`m_0`$ is

$`V_r`$ $`=`$ $`V_{0,r}^{}+V_{01,r}^{}+V`$ (9)
$`=`$ $`K_1\left[\mathrm{cos}(\omega _1+f_1)+e_1\mathrm{cos}\omega _1\right]+K_2\left[\mathrm{cos}(\omega _2+f_2)+e_2\mathrm{cos}\omega _2\right]+V,`$
where $`V`$ is the LOS velocity of the center of mass of the whole system relative to the observer. Equation (9) is exactly the formula used by observers in two-Kepler fits, but with the amplitudes defined in equations (4) and (8) and the orbital periods defined in equations (2) and (6).333 Note that the planetary masses in the equations for the amplitudes and the orbital periods are the physical masses $`m_j`$ and not the Jacobi masses $`m_j_{k=0}^{j1}m_k/_{k=0}^jm_k`$ as stated by Lissauer & Rivera (2001).  Since the true anomaly $`f_j`$ depends on $`P_j`$, $`e_j`$, and the time of periapse passage $`T_{\mathrm{peri},j}`$, a two-Kepler fit directly yields five parameters $`P_j`$, $`K_j`$, $`e_j`$, $`\omega _j`$, and $`T_{\mathrm{peri},j}`$ for each orbit. Table 1 shows these parameters for the HD 168443 planets from Marcy et al. (2001) and the HD 12661 planets from D. A. Fischer (2002, private communication). If the stellar mass $`m_0`$ is known, equations (2), (4), (6), and (8) can be used to derive $`m_j`$ and $`a_j`$ for assumed $`\mathrm{sin}i_j`$. (A common but less accurate practice is to derive $`m_j\mathrm{sin}i_j`$ and $`a_j`$ by assuming that $`m_1`$ and $`m_2`$ are negligible compared to $`m_0`$.) Table 1 also shows $`m_j`$ (in units of Jupiter mass $`M_J`$) and $`a_j`$ for $`\mathrm{sin}i_j=1`$. Throughout this paper, we refer to the substellar companions of HD 168443 as planets, but it should be noted that their $`\mathrm{sin}i_j=1`$ (i.e., minimum) masses are about $`7.7`$ and $`17M_J`$, which are near or above the deuterium-burning limit.

The alternative assumption that the orbits of the planets in astrocentric coordinates are unperturbed Kepler orbits would also yield equation (9) for the radial velocity of the host star, but with $`K_j=[2\pi G(m_0+m_j)/P_j]^{1/3}m_j\mathrm{sin}i_j(m_0+m_1+m_2)^1(1e_j^2)^{1/2}`$ and $`P_j=2\pi a_j^{3/2}[G(m_0+m_j)]^{1/2}`$. However, as we show next, while the Jacobi orbits of the planets of a hierarchical system are nearly Keplerian on orbital timescales, the astrocentric orbit of the outer planet can deviate significantly from a Kepler orbit on orbital timescales.

Figure 1a shows the variations in the semimajor axes and eccentricities of the HD 168443 planets if we assume that the orbital parameters obtained from the two-Kepler fit are in astrocentric coordinates and plot the astrocentric $`a_j`$ and $`e_j`$. The planets are assumed to be on coplanar orbits with $`\mathrm{sin}i=1`$ (note that coplanar orbits have the same inclinations, $`i_1=i_2=i`$, and the same longitudes of ascending node, $`\mathrm{\Omega }_2=\mathrm{\Omega }_1`$, referenced to the plane of the sky). The solid and dotted lines are from direct numerical orbit integrations starting at $`T_{\mathrm{peri},1}`$ and $`T_{\mathrm{peri},2}`$ of Table 1, respectively. The direct integrations were performed using the Wisdom-Holman (1991) integrator contained in the SWIFT444 See http://www.boulder.swri.edu/$∼$hal/swift.html.  software package, with input and output in astrocentric coordinates. Figure 1b is similar to Figure 1a, but it shows the variations in the Jacobi $`a_j`$ and $`e_j`$, with the initial orbital parameters also in Jacobi coordinates. The direct integrations were performed using the modified Wisdom-Holman integrator described in § 4. In astrocentric coordinates (Fig. 1a), the evolution of the orbital elements is sensitive to the starting epoch; and while $`a_1`$ is nearly constant and $`e_1`$ oscillates only on the secular timescale, both $`a_2`$ and $`e_2`$ show significant fluctuations on the orbital timescales (note that we have reduced overcrowding in Fig. 1a by using a sampling interval of $`50\mathrm{yr}`$, which is long compared to the orbital timescales). Plots of the variations in the orbital elements of the HD 168443 planets in Marcy et al. (2001), Nagasawa, Lin, & Ida (2003), and Udry, Mayor, & Queloz (2003, from calculations by W. Benz) show similar behaviors and are probably due to the use of astrocentric coordinates. In contrast, in Jacobi coordinates (Fig. 1b), the evolution is not sensitive to the starting epoch (the solid and dotted lines in Fig. 1b are almost indistinguishable); both $`a_1`$ and $`a_2`$ are nearly constant, while both $`e_1`$ and $`e_2`$ oscillate only on the secular timescale, with the maximum in $`e_1`$ coinciding with the minimum in $`e_2`$ and vice versa. As we shall see, these behaviors can be understood with the octupole-level secular perturbation theory for coplanar hierarchical two-planet systems, which is also based on Jacobi coordinates.

From the facts that the Jacobi $`a_j`$ are nearly constant and that the Jacobi $`e_j`$ oscillate only on the secular timescale for hierarchical systems, it is easy to estimate the fractional fluctuation in the astrocentric $`a_2`$, which is primarily due to the velocity of the star $`m_0`$ relative to the center of mass of $`m_0`$ and $`m_1`$. It is approximately $`4(m_1/m_0)[(a_2/a_1)(1+e_{1,\mathrm{max}})/(1e_{1,\mathrm{max}})]^{1/2}`$ when the Jacobi $`e_1`$ is at its maximum $`e_{1,\mathrm{max}}`$ (and $`e_2`$ is at its minimum, which is assumed to be small and is neglected). The fluctuations in the astrocentric $`a_2`$ in Figure 1a are in good agreement with this estimate. Note that the fractional fluctuation is larger for smaller $`a_1/a_2`$, i.e., for systems that are more hierarchical.

In their study of the three-planet $`\upsilon `$ And system, Rivera & Lissauer (2000) have also reported that the use of Jacobi coordinates eliminates the high-frequency variations in the semimajor axes and eccentricities of the outer two planets and reduces the sensitivity of the evolution to the initial epoch. It is clear that Jacobi coordinates should be used in the analysis of hierarchical systems such as HD 168443 (where $`a_1/a_20.10`$) and $`\upsilon `$ And (where $`a_1/a_20.071`$ and $`a_2/a_30.33`$). It is likely that Jacobi coordinates should also be used in the analysis of other types of planetary systems, especially those for which multiple Kepler fits are adequate.

## 3 OCTUPOLE-LEVEL SECULAR PERTURBATION THEORY

As we mentioned in § 1, an approximate theory that describes the secular evolution of coplanar hierarchical two-planet systems in Jacobi coordinates, such as that shown in Figure 1b for the HD 168443 system, is the octupole-level secular perturbation theory. In this section we summarize the derivation of this octupole theory, compare it to the classical Laplace-Lagrange secular perturbation theory, and deduce some results from the octupole equations analytically.

### 3.1 Equations

As in §2, we consider a system consisting of a central star of mass $`m_0`$, an inner planet of mass $`m_1`$, and an outer planet of mass $`m_2`$, and use Jacobi coordinates, with $`𝒓_1`$ being the position of $`m_1`$ relative to $`m_0`$ and $`𝒓_2`$ being the position of $`m_2`$ relative to the center of mass of $`m_0`$ and $`m_1`$. The Hamiltonian of this system is

$$=\frac{Gm_0m_1}{2a_1}\frac{G(m_0+m_1)m_2}{2a_2}Gm_0m_2\left(\frac{1}{r_{02}}\frac{1}{r_2}\right)Gm_1m_2\left(\frac{1}{r_{12}}\frac{1}{r_2}\right),$$
(10)
where $`a_j`$ is the osculating semimajor axis of the $`j`$th orbit (with $`j=1`$ and $`2`$ for the inner and outer orbits, respectively), and $`r_{k2}`$ is the distance between $`m_k`$ and $`m_2`$. With $`r_1<r_2`$, both $`1/r_{02}`$ and $`1/r_{12}`$ can be expanded in powers of $`r_1/r_2`$, leading to the expanded Hamiltonian

$$=\frac{Gm_0m_1}{2a_1}\frac{G(m_0+m_1)m_2}{2a_2}\frac{G}{a_2}\underset{k=2}{\overset{\mathrm{}}{}}\alpha ^kM_k\left(\frac{r_1}{a_1}\right)^k\left(\frac{a_2}{r_2}\right)^{k+1}P_k(\mathrm{cos}\mathrm{\Phi }),$$
(11)
where $`\alpha =a_1/a_2`$,

$$M_k=m_0m_1m_2\frac{m_0^{k1}(m_1)^{k1}}{(m_0+m_1)^k},$$
(12)
$`P_k`$ is the Legendre polynomial of degree $`k`$, and $`\mathrm{\Phi }`$ is the angle between $`𝒓_1`$ and $`𝒓_2`$. The first two terms in equations (10) and (11) represent the independent Kepler motions of the inner and outer orbits, while the remaining terms represent the perturbations to the Kepler motions. For hierarchical systems with $`r_1r_2`$, the Jacobi decomposition leads to two slowly perturbed Kepler orbits, even if $`m_1`$ and $`m_2`$ are not much smaller than $`m_0`$.

In the general case where mutually inclined orbits are allowed, it is convenient to use the Delaunay variables with the invariable plane as the reference plane. The coordinates are the mean anomalies $`l_j`$, the arguments of periapse $`g_j=\omega _j`$, and the longitudes of ascending node $`h_j=\mathrm{\Omega }_j`$, and the conjugate momenta are

$`L_1`$ $`=`$ $`{\displaystyle \frac{m_0m_1}{m_0+m_1}}\sqrt{G(m_0+m_1)a_1},`$ (13)
$`L_2`$ $`=`$ $`{\displaystyle \frac{(m_0+m_1)m_2}{m_0+m_1+m_2}}\sqrt{G(m_0+m_1+m_2)a_2},`$ (14)
$`G_j`$ $`=`$ $`L_j\sqrt{1e_j^2},`$ (15)
$`H_j`$ $`=`$ $`G_j\mathrm{cos}i_j,`$ (16)
where $`e_j`$ and $`i_j`$ are the eccentricities and inclinations of the orbits. The momenta $`L_j`$, $`G_j`$, and $`H_j`$ are, respectively, the magnitude of the maximum possible angular momentum (if the orbit were circular), the magnitude of the angular momentum, and the $`z`$-component of the angular momentum of the $`j`$th orbit. By expressing $`\mathrm{cos}\mathrm{\Phi }=𝒓_1\mathbf{}𝒓_2/(r_1r_2)`$ in terms of the Delaunay variables, it is easy to show that the Hamiltonian contains $`h_1`$ and $`h_2`$ only in the combination $`h_1h_2`$, and $`h_1h_2=180^{}`$ if the invariable plane is used as the reference plane. Constant $`h_1h_2`$ and constant total angular momentum $`H_1+H_2`$ mean $`h_j`$ and $`H_j`$ can be eliminated from the problem — usually referred to as the elimination of nodes — and the system is reduced from $`6`$ to $`4`$ degrees of freedom (see, e.g., Marchal 1990).

Kozai (1962) studied the secular evolution of hierarchical triple systems with $`m_1m_2m_0`$ and $`e_2=0`$ by retaining just the lowest-order ($`\alpha ^2`$) quadrupole term in the series in equation (11) and averaging the Hamiltonian over the inner and outer orbital motions. Harrington (1968) generalized the quadrupole analysis to general hierarchical triple systems. Marchal (1990), Krymolowski & Mazeh (1999), and Ford et al. (2000) extended the approximation by retaining also the octupole term of order $`\alpha ^3`$ in the series in equation (11) and derived an octupole-level orbit-averaged Hamiltonian. The approach used in these studies to average the Hamiltonian is the von Zeipel method, which involves the determination of a canonical transformation such that the transformed Hamiltonian is independent of the transformed mean anomalies $`l_1`$ and $`l_2`$, and such that the original and transformed variables differ only by high frequency terms that are multiplied by powers of the small parameter $`\alpha `$. Since the transformed Hamiltonian is independent of the transformed $`l_j`$, the transformed $`L_j`$ and the corresponding semimajor axes are constant.

In this paper we focus on hierarchical two-planet systems with coplanar orbits, as the formation of planets from a common disk of materials surrounding the host star makes nearly coplanar orbits the most probable configuration. In the limit of coplanar orbits ($`i_1=i_2=0^{}`$), the Delaunay variables described above are not well-defined, since there are no ascending nodes with respect to the invariable plane, but the longitudes of periapse $`\varpi _j`$ are well-defined. Thus, for coplanar orbits, it is convenient to use the canonical variables $`l_j`$, $`\varpi _j`$, $`L_j`$, and $`G_j`$. In the limit of coplanar orbits, by noting that the expression for the angle $`\phi `$ between the directions of periapse reduces to $`\mathrm{cos}\phi =\mathrm{cos}(g_1g_2)=\mathrm{cos}(g_1g_2+180^{})=\mathrm{cos}(\varpi _1\varpi _2)`$ (since $`h_1h_2=180^{}`$), the doubly averaged octupole-level Hamiltonian derived by Marchal (1990), Krymolowski & Mazeh (1999), and Ford et al. (2000) can be written as555 There is a sign error common to both Krymolowski & Mazeh (1999) and Ford et al. (2000). We follow Ford et al. and denote the coefficients of the quadrupole and octupole terms by $`C_2`$ and $`C_3`$, respectively. Krymolowski & Mazeh denote these coefficients by $`C_1`$ and $`C_2`$, respectively. In the averaged Hamiltonian and the subsequent equations of motion of Ford et al., the coefficient $`C_3`$ should be replaced by $`C_3`$. Similarly, in the averaged Hamiltonian and the subsequent equations of motion of Krymolowski & Mazeh, the coefficient $`C_2`$ should be replaced by $`C_2`$.

$$\overline{}_{\mathrm{oct}}=\frac{G^2m_0^3m_1^3}{2(m_0+m_1)L_1^2}\frac{G^2(m_0+m_1)^3m_2^3}{2(m_0+m_1+m_2)L_2^2}2C_2(2+3e_1^2)+C_3e_1e_2(4+3e_1^2)\mathrm{cos}(\varpi _1\varpi _2),$$
(17)
where

$`C_2`$ $`=`$ $`{\displaystyle \frac{1}{16}}{\displaystyle \frac{G^2(m_0+m_1)^7m_2^7}{(m_0+m_1+m_2)^3(m_0m_1)^3}}{\displaystyle \frac{L_1^4}{L_2^3G_2^3}},`$ (18)
$`C_3`$ $`=`$ $`{\displaystyle \frac{15}{64}}{\displaystyle \frac{G^2(m_0+m_1)^9m_2^9(m_0m_1)}{(m_0+m_1+m_2)^4(m_0m_1)^5}}{\displaystyle \frac{L_1^6}{L_2^3G_2^5}},`$ (19)
and

$$e_j=\sqrt{1\left(\frac{G_j}{L_j}\right)^2}.$$
(20)
As in Ford et al. (2000), we do not include in equation (17) terms of order $`\alpha ^{7/2}`$ induced by the canonical transformation of the von Zeipel method, which were included partially by Marchal (1990) and fully by Krymolowski & Mazeh (1999).

As we mentioned above, the orbit averaging eliminates $`l_j`$ from $`\overline{}_{\mathrm{oct}}`$, and hence $`L_j`$ (and $`a_j`$) are constant. Furthermore, $`\overline{}_{\mathrm{oct}}`$ contains $`\varpi _1`$ and $`\varpi _2`$ only in the combination $`\varpi _1\varpi _2`$, and choosing this difference as a new variable in a canonical transformation simply reveals that the total angular momentum $`G_1+G_2`$ is an integral of motion. Thus, in the octupole approximation, the coplanar problem is reduced to one degree of freedom with the last two terms in equation (17) as an integral of motion. We can also see from equation (17) that the quadrupole approximation is not adequate for studying coplanar systems. If we drop the octupole term with coefficient $`C_3`$ in equation (17), the Hamiltonian is independent of $`\varpi _1`$ and $`\varpi _2`$, meaning that both $`e_1`$ and $`e_2`$ are constant, which is not generally true. (For mutually inclined orbits, only $`e_2`$ is constant in the quadrupole approximation; Harrington 1968.)

From the equations of motion

$$\frac{d\varpi _j}{dt}=\frac{\overline{}_{\mathrm{oct}}}{G_j}\mathrm{and}\frac{dG_j}{dt}=\frac{\overline{}_{\mathrm{oct}}}{\varpi _j},$$
(21)
we obtain the following equations for the variation of $`e_j`$ and $`\varpi _j`$:

$`{\displaystyle \frac{de_1}{dt}}`$ $`=`$ $`A_{12}e_2{\displaystyle \frac{(1e_1^2)^{1/2}\left(1+\frac{3}{4}e_1^2\right)}{(1e_2^2)^{5/2}}}\mathrm{sin}(\varpi _1\varpi _2),`$ (22)
$`{\displaystyle \frac{de_2}{dt}}`$ $`=`$ $`A_{21}e_1{\displaystyle \frac{\left(1+\frac{3}{4}e_1^2\right)}{(1e_2^2)^2}}\mathrm{sin}(\varpi _1\varpi _2),`$ (23)
$`{\displaystyle \frac{d\varpi _1}{dt}}`$ $`=`$ $`A_{11}{\displaystyle \frac{(1e_1^2)^{1/2}}{(1e_2^2)^{3/2}}}A_{12}\left({\displaystyle \frac{e_2}{e_1}}\right){\displaystyle \frac{(1e_1^2)^{1/2}\left(1+\frac{9}{4}e_1^2\right)}{(1e_2^2)^{5/2}}}\mathrm{cos}(\varpi _1\varpi _2),`$ (24)
$`{\displaystyle \frac{d\varpi _2}{dt}}`$ $`=`$ $`A_{22}{\displaystyle \frac{\left(1+\frac{3}{2}e_1^2\right)}{(1e_2^2)^2}}A_{21}\left({\displaystyle \frac{e_1}{e_2}}\right){\displaystyle \frac{(1+4e_2^2)\left(1+\frac{3}{4}e_1^2\right)}{(1e_2^2)^3}}\mathrm{cos}(\varpi _1\varpi _2),`$ (25)
where the constant coefficients

$`A_{12}`$ $`=`$ $`{\displaystyle \frac{4C_3}{L_1}}\left({\displaystyle \frac{G_2}{L_2}}\right)^5={\displaystyle \frac{15}{16}}n_1\left({\displaystyle \frac{m_2}{m_0+m_1}}\right)\left({\displaystyle \frac{m_0m_1}{m_0+m_1}}\right)\alpha ^4,`$ (26a)
$`A_{21}`$ $`=`$ $`{\displaystyle \frac{4C_3}{L_2}}\left({\displaystyle \frac{G_2}{L_2}}\right)^5={\displaystyle \frac{15}{16}}n_2{\displaystyle \frac{m_0m_1}{(m_0+m_1)^2}}\left({\displaystyle \frac{m_0m_1}{m_0+m_1}}\right)\alpha ^3,`$ (26b)
$`A_{11}`$ $`=`$ $`{\displaystyle \frac{12C_2}{L_1}}\left({\displaystyle \frac{G_2}{L_2}}\right)^3={\displaystyle \frac{3}{4}}n_1\left({\displaystyle \frac{m_2}{m_0+m_1}}\right)\alpha ^3,`$ (26c)
$`A_{22}`$ $`=`$ $`{\displaystyle \frac{12C_2}{L_2}}\left({\displaystyle \frac{G_2}{L_2}}\right)^3={\displaystyle \frac{3}{4}}n_2{\displaystyle \frac{m_0m_1}{(m_0+m_1)^2}}\alpha ^2,`$ (26d)
with $`\alpha =a_1/a_2`$ and the mean motions

$$n_1=\left[G(m_0+m_1)/a_1^3\right]^{1/2}\mathrm{and}n_2=\left[G(m_0+m_1+m_2)/a_2^3\right]^{1/2}.$$
(27)
Consistent with the fact that the system can be reduced to one degree of freedom, either equation (22) or (23) is redundant, since $`e_1`$ and $`e_2`$ are related by the conservation of total angular momentum, and equations (24) and (25) can be combined as

$`{\displaystyle \frac{d(\varpi _1\varpi _2)}{dt}}`$ $`=`$ $`A_{11}\left[{\displaystyle \frac{(1e_1^2)^{1/2}}{(1e_2^2)^{3/2}}}\left({\displaystyle \frac{L_1}{L_2}}\right){\displaystyle \frac{\left(1+\frac{3}{2}e_1^2\right)}{(1e_2^2)^2}}\right]`$ (28)
$`A_{12}\left[\left({\displaystyle \frac{e_2}{e_1}}\right){\displaystyle \frac{(1e_1^2)^{1/2}\left(1+\frac{9}{4}e_1^2\right)}{(1e_2^2)^{5/2}}}\left({\displaystyle \frac{L_1}{L_2}}\right)\left({\displaystyle \frac{e_1}{e_2}}\right){\displaystyle \frac{(1+4e_2^2)\left(1+\frac{3}{4}e_1^2\right)}{(1e_2^2)^3}}\right]`$
$`\times \mathrm{cos}(\varpi _1\varpi _2).`$
If we relax the assumption of coplanar orbits and allow a small mutual inclination angle $`i_{\mathrm{mu}}`$ between the inner and outer orbits, the number of degrees of freedom is increased from one to two, but an expansion of the orbit-averaged octupole-level Hamiltonian of Marchal (1990), Krymolowski & Mazeh (1999), and Ford et al. (2000) in powers of $`i_{\mathrm{mu}}`$ shows that equations (22)–(25) are only modified by terms of order $`i_{\mathrm{mu}}^2`$ and higher. One of the additional terms for $`de_1/dt`$ (which has an octupole term only in the coplanar limit; see eq. ) is a quadrupole term with coefficient $`C_2i_{\mathrm{mu}}^2`$. Thus we expect a small mutual inclination angle $`i_{\mathrm{mu}}\alpha ^{1/2}`$ to have little effect on the secular evolution of $`e_j`$ and $`\varpi _j`$ of a hierarchical system. This is confirmed by direct numerical orbit integrations in § 4.

It is important to note that the octupole-level secular perturbation theory is derived under the assumptions that there are no mean-motion commensurabilities and that $`\alpha `$ (or more precisely $`r_1/r_2`$) $`1`$. As we shall see in § 4, close proximity to even a rather high order commensurability (e.g., 11:2) can affect the evolution of a system. We shall also see in § 4 that although there is a limit to how large $`\alpha `$ can be for a system to be stable, the octupole theory can provide a reasonable description of the secular evolution for $`\alpha `$ as large as $`1/3`$ if the planets are not too massive.

### 3.2 Comparison with Classical Laplace-Lagrange Secular Perturbation Theory

As we mentioned in § 1, the classical secular perturbation theory developed for the solar system is based on an expansion in the eccentricities and inclinations. While it is valid to all orders in the ratio of the semimajor axes $`\alpha `$, it can be applied only to systems with small planetary masses on nearly circular and nearly coplanar orbits. In contrast, the octupole-level secular perturbation theory is based on an expansion in $`\alpha `$, and it can be applied to hierarchical systems with a wide range of masses, eccentricities and, in its most general form, inclinations. We now show that the octupole equations (22)–(25) agree with the equations for the classical Laplace-Lagrange secular solution for coplanar systems where the planetary masses, eccentricities, and $`\alpha `$ are all small.

To the lowest order in the eccentricities, equations (22)–(25) reduce to

$`{\displaystyle \frac{de_1}{dt}}`$ $`=`$ $`A_{12}e_2\mathrm{sin}(\varpi _1\varpi _2),`$ (29)
$`{\displaystyle \frac{de_2}{dt}}`$ $`=`$ $`A_{21}e_1\mathrm{sin}(\varpi _1\varpi _2),`$ (30)
$`{\displaystyle \frac{d\varpi _1}{dt}}`$ $`=`$ $`A_{11}A_{12}\left({\displaystyle \frac{e_2}{e_1}}\right)\mathrm{cos}(\varpi _1\varpi _2),`$ (31)
$`{\displaystyle \frac{d\varpi _2}{dt}}`$ $`=`$ $`A_{22}A_{21}\left({\displaystyle \frac{e_1}{e_2}}\right)\mathrm{cos}(\varpi _1\varpi _2),`$ (32)
where $`A_{12}=(15/16)(m_2/m_0)n_1\alpha ^4`$, $`A_{21}=(15/16)(m_1/m_0)n_2\alpha ^3`$, $`A_{11}=(3/4)(m_2/m_0)n_1\alpha ^3`$, and $`A_{22}=(3/4)(m_1/m_0)n_2\alpha ^2`$ in the limit $`m_1,m_2m_0`$.

The classical Laplace-Lagrange secular solution is based on retaining just the secular terms in the disturbing functions up to second order in the eccentricities and inclinations. For a coplanar planetary system, the equations for the variation of $`e_j`$ and $`\varpi _j`$ (see, e.g., Murray & Dermott 1999) are of the same form as equations (29)–(32), but with the $`A_{jk}`$ replaced by

$`A_{12}^{}`$ $`=`$ $`{\displaystyle \frac{1}{4}}n_1{\displaystyle \frac{m_2}{m_0+m_1}}\alpha ^2b_{3/2}^{(2)}(\alpha ),`$ (33a)
$`A_{21}^{}`$ $`=`$ $`{\displaystyle \frac{1}{4}}n_2{\displaystyle \frac{m_1}{m_0+m_2}}\alpha b_{3/2}^{(2)}(\alpha ),`$ (33b)
$`A_{11}^{}`$ $`=`$ $`{\displaystyle \frac{1}{4}}n_1{\displaystyle \frac{m_2}{m_0+m_1}}\alpha ^2b_{3/2}^{(1)}(\alpha ),`$ (33c)
$`A_{22}^{}`$ $`=`$ $`{\displaystyle \frac{1}{4}}n_2{\displaystyle \frac{m_1}{m_0+m_2}}\alpha b_{3/2}^{(1)}(\alpha ),`$ (33d)
where $`b_{3/2}^{(j)}(\alpha )`$ is the Laplace coefficient. (It should be noted that the classical theory uses astrocentric coordinates instead of Jacobi coordinates, but that the distinction between the two sets of coordinates vanishes in the limit $`m_1/m_0\alpha ^{1/2}/4`$; see §2.) Since $`b_{3/2}^{(2)}(\alpha )=(15/4)\alpha ^2`$ and $`b_{3/2}^{(1)}(\alpha )=3\alpha `$ to the lowest order in $`\alpha `$, the classical and octupole secular perturbation equations are identical in the limit of small planetary masses, eccentricities, and $`\alpha `$.

### 3.3 Some Analytic Results

Although it does not appear that the octupole-level secular perturbation equations can be solved analytically, some useful results can be deduced analytically. We begin by rewriting equations (22), (23), and (28) as

$`{\displaystyle \frac{de_1}{d\tau }}`$ $`=`$ $`\beta e_2{\displaystyle \frac{(1e_1^2)^{1/2}\left(1+\frac{3}{4}e_1^2\right)}{(1e_2^2)^{5/2}}}\mathrm{sin}(\varpi _1\varpi _2),`$ (34)
$`{\displaystyle \frac{de_2}{d\tau }}`$ $`=`$ $`\beta \lambda e_1{\displaystyle \frac{\left(1+\frac{3}{4}e_1^2\right)}{(1e_2^2)^2}}\mathrm{sin}(\varpi _1\varpi _2),`$ (35)
$`{\displaystyle \frac{d(\varpi _1\varpi _2)}{d\tau }}`$ $`=`$ $`\left[{\displaystyle \frac{(1e_1^2)^{1/2}}{(1e_2^2)^{3/2}}}\lambda {\displaystyle \frac{\left(1+\frac{3}{2}e_1^2\right)}{(1e_2^2)^2}}\right]`$ (36)
$`\beta \left[\left({\displaystyle \frac{e_2}{e_1}}\right){\displaystyle \frac{(1e_1^2)^{1/2}\left(1+\frac{9}{4}e_1^2\right)}{(1e_2^2)^{5/2}}}\lambda \left({\displaystyle \frac{e_1}{e_2}}\right){\displaystyle \frac{(1+4e_2^2)\left(1+\frac{3}{4}e_1^2\right)}{(1e_2^2)^3}}\right]`$
$`\times \mathrm{cos}(\varpi _1\varpi _2),`$
where

$`\beta `$ $`=`$ $`{\displaystyle \frac{A_{12}}{A_{11}}}={\displaystyle \frac{A_{21}}{A_{22}}}={\displaystyle \frac{5}{4}}\left({\displaystyle \frac{m_0m_1}{m_0+m_1}}\right)\alpha ,`$ (37)
$`\lambda `$ $`=`$ $`{\displaystyle \frac{A_{22}}{A_{11}}}={\displaystyle \frac{A_{21}}{A_{12}}}={\displaystyle \frac{L_1}{L_2}},`$ (38)
$`\tau `$ $`=`$ $`t/t_e,`$ (39)
and

$$t_e=\frac{1}{A_{11}}=\frac{4}{3\alpha ^3}\left(\frac{m_0+m_1}{m_2}\right)\frac{1}{n_1}.$$
(40)
Recall that either $`e_1`$ or $`e_2`$ can be eliminated from the above equations by the conservation of total angular momentum, which can be written in dimensionless form as $`\gamma =(G_1+G_2)/(L_1+L_2)=[\lambda (1e_1^2)^{1/2}+(1e_2^2)^{1/2}]/(\lambda +1)=`$ constant.

We can see from equations (34)–(36) that coplanar hierarchical systems with the same $`\beta `$ and $`\lambda `$ (which are constant) and the same initial $`e_1`$, $`e_2`$, and $`\varpi _1\varpi _2`$ have the same trajectory in the phase-space diagram of $`e_1`$ (or $`e_2`$) versus $`\varpi _1\varpi _2`$ and can differ only in the period of eccentricity oscillations. In particular, in the limit $`m_1,m_2m_0`$, since $`\beta (5/4)\alpha `$, $`\lambda (m_1/m_2)\alpha ^{1/2}`$, and $`t_e4(m_0/m_2)/(3\alpha ^3n_1)`$, systems with the same $`\alpha `$, $`m_1/m_2`$, and initial $`e_1`$, $`e_2`$ and $`\varpi _1\varpi _2`$, but different $`m_0`$, $`m_2`$ and/or $`a_1`$, differ only in the period of eccentricity oscillations, which is proportional to $`(m_0/m_2)/n_1`$. As we discussed in § 2, a two-Kepler fit yields $`e_j`$, $`\omega _j`$, and in the limit $`m_1,m_2m_0`$, $`m_j\mathrm{sin}i_j`$ and $`a_j`$, but not the inclinations $`i_j`$ of the orbital planes to the plane of the sky. (Hereafter, all inclinations are relative to the plane of the sky and not to the invariable plane.) If we assume that the orbits are coplanar, the above results mean that the trajectory in $`e_1`$ (or $`e_2`$) versus $`\varpi _1\varpi _2`$ and the amplitudes of eccentricity oscillations should be independent of the unknown inclination $`i`$ of both orbits, and that the period of eccentricity oscillations should be proportional to $`\mathrm{sin}i`$.

The phase-space structure of the secular evolution of coplanar hierarchical two-planet systems can be understood by plotting the trajectories of systems with the same $`\beta `$, $`\lambda `$, and $`\gamma `$ in a diagram of $`e_1`$ (or $`e_2`$) versus $`\varpi _1\varpi _2`$ (see, e.g., Figs. 2 and 6 below). The fixed points in the phase-space diagram can be found by solving $`de_1/d\tau =de_2/d\tau =d(\varpi _1\varpi _2)/d\tau =0`$. It is clear from equations (34) and (35) that $`de_1/d\tau =de_2/d\tau =0`$ requires $`\mathrm{sin}(\varpi _1\varpi _2)=0`$ or $`\varpi _1\varpi _2=0^{}`$ or $`180^{}`$. As we shall see, there are usually two libration islands, one about a fixed point at $`\varpi _1\varpi _2=0^{}`$ and another about a fixed point at $`\varpi _1\varpi _2=180^{}`$; additional fixed points are possible if the total angular momentum is low.

At $`\varpi _1\varpi _290^{}`$ and $`270^{}`$, since $`\mathrm{cos}(\varpi _1\varpi _2)0`$, equation (36) reduces to $`d(\varpi _1\varpi _2)/d\tau (1e_1^2)^{1/2}/(1e_2^2)^{3/2}\lambda (1+3e_1^2/2)/(1e_2^2)^2`$. Therefore, to the lowest order in the eccentricities, if $`\lambda 1`$, $`d(\varpi _1\varpi _2)/d\tau 0`$ and $`\varpi _1\varpi _2`$ should be nearly constant while both $`e_1`$ and $`e_2`$ change. According to equations (34) and (35), $`e_1`$ should be decreasing (from near its maximum possible value for the given total angular momentum to near zero) and $`e_2`$ should be increasing (from near zero to near its maximum possible value for the given total angular momentum) at $`\varpi _1\varpi _290^{}`$, where $`\mathrm{sin}(\varpi _1\varpi _2)1`$, and vice versa at $`\varpi _1\varpi _2270^{}`$, where $`\mathrm{sin}(\varpi _1\varpi _2)1`$. These behaviors are most consistent with a phase space that is dominated by two large libration islands, one about a fixed point at $`\varpi _1\varpi _2=0^{}`$ and another about a fixed point at $`\varpi _1\varpi _2=180^{}`$. Note also that for a given total angular momentum the maximum possible values of $`e_1`$ and $`e_2`$ are comparable if $`\lambda 1`$. Finite eccentricities change the condition for $`d(\varpi _1\varpi _2)/d\tau 0`$ at $`\varpi _1\varpi _290^{}`$ and $`270^{}`$ to $`\lambda (1e_1^2)^{1/2}(1e_2^2)^{1/2}/(1+3e_1^2/2)`$, which is slightly less than $`1`$ for moderate eccentricities. To estimate how much smaller than $`1`$ the critical value of $`\lambda `$ is for a given dimensionless total angular momentum $`\gamma `$, we substitute into the above condition the value of the eccentricities when they are equal \[$`e_1=e_2=(1\gamma ^2)^{1/2}`$\] and obtain $`\lambda _{\mathrm{crit}}=2\gamma ^2/(53\gamma ^2)`$. Note that $`\lambda _{\mathrm{crit}}=1`$ is recovered for $`\gamma =1`$. Therefore, large libration islands and large amplitude variations of both $`e_j`$ are likely if $`\lambda \lambda _{\mathrm{crit}}=2\gamma ^2/(53\gamma ^2)`$.

## 4 NUMERICAL RESULTS

In this section we study the dynamical evolution of the HD 168443 and HD 12661 systems and their variants and demonstrate the validity and limits of the octupole-level secular perturbation theory by comparison with direct numerical orbit integrations. Except for the direct integrations discussed in § 4.3, the planets are assumed to be on coplanar orbits. Since the octupole-level secular perturbation equations do not involve high frequency terms, numerical integrations of these equations are rapid. We integrated equations (35) and (36), with $`e_1`$ found from conservation of total angular momentum, using a Bulirsch-Stoer integrator. The direct numerical orbit integrations were performed using a modified version of the Wisdom-Holman (1991) integrator contained in the SWIFT software package. In addition to changing the input and output to Jacobi orbital elements, we divide the Hamiltonian into a part that describes the Kepler motions of the inner and outer orbits and a part that describes the perturbations to the Kepler motions using equation (10) instead of the division used by Wisdom & Holman (1991), which moves the term $`Gm_1m_2/r_2`$ from the perturbation Hamiltonian to the Kepler Hamiltonian of the outer orbit. This modified integrator can handle hierarchical systems where $`m_1`$ and $`m_2`$ are not small.

### 4.1 HD 168443

The best-fit orbital parameters of the HD 168443 planets from Marcy et al. (2001) and the inferred planetary masses and semimajor axes for $`\mathrm{sin}i=1`$ are listed in Table 1. This system has $`\alpha =a_1/a_2=0.102`$, $`\beta =5(m_0m_1)\alpha /[4(m_0+m_1)]=0.126`$, $`\lambda =L_1/L_2=0.143`$, and $`\gamma =(G_1+G_2)/(L_1+L_2)=0.963`$. Figure 1b shows the variations in the Jacobi $`a_j`$ and $`e_j`$, with the solid and dotted lines (almost indistinguishable) from direct numerical orbit integrations starting at $`T_{\mathrm{peri},1}`$ and $`T_{\mathrm{peri},2}`$ of Table 1, respectively, and the dashed lines in the lower two panels from the octupole-level secular perturbation theory. The trajectories in Figure 2 marked by the squares (which denote the current parameters of the HD 168443 system) are the trajectories of the $`\mathrm{sin}i=1`$ HD 168443 system in the phase-space diagrams of $`e_1`$ versus $`\varpi _1\varpi _2`$ (upper panels) and $`e_2`$ versus $`\varpi _1\varpi _2`$ (lower panels), with the left and right panels showing the results from direct integrations and the octupole theory, respectively (see below for a discussion of the other trajectories in Fig. 2). The direct integration results agree with the octupole theory in that the evolution is not sensitive to the starting epoch and that $`a_1`$ and $`a_2`$ are nearly constant. In addition, the trajectories in the diagram of $`e_1`$ (or $`e_2`$) versus $`\varpi _1\varpi _2`$ from the direct integrations are in excellent agreement with that from the octupole theory. The only noticeable difference is that the period of eccentricity oscillations predicted by the octupole theory is about $`3\%`$ longer that that found in the direct integrations ($`1.8\times 10^4\mathrm{yr}`$). Krymolowski & Mazeh (1999) have reported that the oscillation period predicted by the octupole theory can be improved by the inclusion of the terms of order $`\alpha ^{7/2}`$ induced by the canonical transformation of the von Zeipel method and neglected by us.

To study the effects of $`\mathrm{sin}i`$, we have performed two direct integrations of the coplanar HD 168443 system with $`\mathrm{sin}i=0.4`$, one starting at $`T_{\mathrm{peri},1}`$ and the other $`T_{\mathrm{peri},2}`$. (Although the two-Kepler fit to the radial velocity observations does not yield $`\mathrm{sin}i_j`$, the lack of evidence for stellar wobble in the Hipparcos astrometric data limits $`\mathrm{sin}i_20.4`$ for the outer planet of HD 168443; Marcy et al. 2001.) In both cases, the amplitudes of eccentricity oscillations and the trajectory in $`e_1`$ (or $`e_2`$) versus $`\varpi _1\varpi _2`$ are almost identical to those shown in Figures 1b and 2 for $`\mathrm{sin}i=1`$, and the factor by which the eccentricity oscillation period shortens agrees with $`\mathrm{sin}i=0.4`$ to better than $`\pm 0.001`$. These results are in good agreement with the analytic results derived in § 3.3 in the limit $`m_1,m_2m_0`$, even though $`m_2`$ is almost $`44M_J`$ (and $`m_2/m_00.042`$) when $`\mathrm{sin}i=0.4`$.

In addition to the trajectories of the $`\mathrm{sin}i=1`$ HD 168443 system, Figure 2 also shows the trajectories of systems with the same masses, initial semimajor axes, and total angular momentum as the $`\mathrm{sin}i=1`$ HD 168443 system. (Recall that the octupole results are also valid for other systems with the same $`\beta `$, $`\lambda `$, and $`\gamma `$.) There is excellent agreement between the direct-integration and octupole results in all cases. The HD 168443 system is not in a secular resonance and its $`\varpi _1\varpi _2`$ circulates. Indeed it is far from the relatively small libration islands about the fixed points at $`(\varpi _1\varpi _2,e_1)=(0^{},0.046)`$ and $`(180^{},0.702)`$ in the phase-space diagram of $`e_1`$ versus $`\varpi _1\varpi _2`$. The small libration islands and modest eccentricity variations in Figure 2 are consistent with the fact that these systems have $`\lambda =0.143`$ far from $`\lambda _{\mathrm{crit}}=0.836`$ for the dimensionless total angular momentum $`\gamma =0.963`$ of these systems.

The above $`e_1`$ values for the fixed points were obtained from the octupole theory as the roots of $`d(\varpi _1\varpi _2)/d\tau =0`$ at $`\varpi _1\varpi _2=0^{}`$ and $`180^{}`$, with $`d(\varpi _1\varpi _2)/d\tau `$ from equation (36) and $`e_2`$ from the conservation of total angular momentum. As we discussed in § 3.3, the octupole theory predicts that the fixed points must be at $`\varpi _1\varpi _2=0^{}`$ or $`180^{}`$. To demonstrate how the octupole theory can be used to explore the parameter space rapidly, we have used the same root finding procedure to determine the number and positions of fixed points as a function of the dimensionless total angular momentum $`\gamma =(G_1+G_2)/(L_1+L_2)`$ for systems with the same $`\beta `$ and $`\lambda `$ as the $`\mathrm{sin}i=1`$ HD 168443 system. The $`e_1`$ values for the fixed points at $`\varpi _1\varpi _2=0^{}`$ and $`180^{}`$ as a function of $`\gamma `$ are shown in the left and right panels of Figure 3, respectively. For all values of $`\gamma `$ shown in Figure 3, there is an elliptic fixed point at $`\varpi _1\varpi _2=0^{}`$ with relatively small $`e_1`$ (lower left panel) and an elliptic fixed point at $`\varpi _1\varpi _2=180^{}`$ with $`e_1`$ close to the maximum possible value indicated by the dashed line (right panel). These two fixed points are the ones seen above for systems with the same $`\gamma `$($`=0.963`$; dotted lines in Fig. 3) as the $`\mathrm{sin}i=1`$ HD 168443 system.

For $`0.872\gamma 0.8818`$, there are two additional fixed points at $`\varpi _1\varpi _2=0^{}`$ with $`e_1`$ very close to $`1`$ (upper left panel of Fig. 3). These additional fixed points emerge at $`\gamma 0.8818`$ with the same $`e_1`$ value (square in the upper left panel of Fig. 3), and the one with $`e_1`$ value that increases (decreases) with decreasing $`\gamma `$ is an elliptic (hyperbolic) fixed point. Direct numerical orbit integrations confirm that there are indeed stable librations about an elliptic fixed point close to the additional one predicted by the octupole theory. For example, for $`\gamma =0.88`$, the octupole theory predicts that the additional elliptic fixed point is at $`e_1=0.9948`$ (and $`e_2=0.1302`$). A direct integration of a system with these values of $`e_1`$, $`e_2`$, and $`\varpi _1\varpi _2`$ as initial conditions (and with the masses and semimajor axes of the $`\mathrm{sin}i=1`$ HD 168443 system and mean anomalies $`l_1=0^{}`$ and $`l_2=180^{}`$) shows that this system has $`\varpi _1\varpi _2`$ librating about $`0^{}`$ with an amplitude of about $`6^{}`$ and $`e_1`$ varying between $`0.9946`$ and $`0.9969`$. It should be noted that the treatment of this system with $`e_10.9946`$ and $`a_10.3\mathrm{AU}`$ as point masses interacting via Newtonian gravity is inadequate since the periapse distance of the inner planet is in fact less than the stellar radius. However, based on the octupole theory, we expect variants of this system with much larger $`a_1`$ (so that the periapse distance of the inner planet is well outside the stellar radius) but the same $`\alpha `$ to show similar libration behaviors. On the other hand, the inner planet could still come sufficiently close to the star for general relativistic precession and tidal effects to be important. For $`\gamma 0.872`$, the additional elliptic fixed point vanishes, leaving only the hyperbolic fixed point (upper left panel of Fig. 3). Both octupole and direct-integration calculations show that systems with $`\gamma 0.872`$, initial $`\varpi _1\varpi _2=0^{}`$ and initial $`e_1`$ above the hyperbolic fixed point have $`e_1`$ increasing to unity in a finite time. Again, we expect variants of these systems with larger $`a_1`$ (so that the periapse distance of the inner planet is well outside the stellar radius) but the same $`\alpha `$ to show similar increase in $`e_1`$. However, before $`e_1`$ reaches unity, either the effects of tides and general relativistic precession would eventually stop the increase in $`e_1`$, or the inner planet would collide with the star.

### 4.2 HD 12661

The best-fit orbital parameters of the HD 12661 planets from D. A. Fischer (2002, private communication) and the inferred planetary masses and semimajor axes for $`\mathrm{sin}i=1`$ are listed in Table 1. Figure 4 shows the variations in the Jacobi $`a_j`$ and $`e_j`$, with the solid and dotted lines from direct numerical orbit integrations starting at $`T_{\mathrm{peri},1}`$ and $`T_{\mathrm{peri},2}`$ of Table 1, respectively, and the dashed lines in the lower two panels from the octupole-level secular perturbation theory. The direct integration results are not consistent with the secular theory in that the evolution of the orbital elements is sensitive to the starting epoch. It turns out that a system with exactly the best-fit orbital parameters is very close to the 11:2 mean-motion commensurability ($`P_2/P_1=5.486=0.9975\times 11/2`$). The direct integration starting at $`T_{\mathrm{peri},1}`$ (solid lines in Fig. 4) shows irregular fluctuations in $`a_j`$ and $`e_j`$ (note, e.g., the irregular jumps in the mean values of $`a_1`$ and $`a_2`$ at successive minima of $`e_2`$) and is most likely chaotic. The direct integration starting at $`T_{\mathrm{peri},2}`$ (dotted lines in Fig. 4) does not show any obvious irregular jumps in $`a_j`$ or $`e_j`$, and may be either regular or very weakly chaotic. But its smaller amplitudes of eccentricity variations are due to the 11:2 commensurability. The chaos in one (and possibly both) of these calculations is due to the overlap of the resonances at the 11:2 commensurability (see Holman & Murray 1996 and Murray & Holman 1997 for a similar situation in the planar elliptic restricted three-body problem). For both of the direct integrations shown in Figure 4, we have examined the 10 eccentricity-type mean-motion resonance variables at the 11:2 commensurability and confirmed that some of them alternate between circulation and libration. We have extended the direct integrations shown in Figure 4 and found that both are stable for at least $`10^6\mathrm{yr}`$, but we cannot rule out the possibility that the chaos would lead to instability on longer timescales.

The orbital period $`P_2`$ (and the other orbital parameters) of the outer planet of HD 12661 are not currently known to high precision, because the time span of the available observations is comparable to $`P_2`$. To study the effects of varying $`P_2`$, we have performed two sets of direct integrations with different initial $`P_2`$, one starting at $`T_{\mathrm{peri},1}`$ and the other at $`T_{\mathrm{peri},2}`$. The initial $`P_2`$ were chosen such that $`(P_2/P_1)/(11/2)=0.98`$, $`0.99`$, $`1.01`$, $`1.02`$, and $`1.03`$. The initial values of the other orbital parameters that can be obtained from the two-Kepler fit ($`P_1`$, $`K_{1,2}`$, $`e_{1,2}`$, $`\omega _{1,2}`$, and $`T_{\mathrm{peri},1,2}`$) were fixed at the values listed in Table 1, and the planetary masses and initial semimajor axes were derived assuming that $`\mathrm{sin}i=1`$. We find that the cases with initial $`(P_2/P_1)/(11/2)=1.01`$ are also affected by the 11:2 mean-motion commensurability while those with initial $`(P_2/P_1)/(11/2)=0.98`$, $`0.99`$, and $`1.02`$ show regular secular evolution. Although the variations in $`a_j`$ for the cases with initial $`P_2/P_1=1.03\times 11/2=0.9997\times 17/3`$ indicate that these cases are probably affected by the 17:3 commensurability, the variations in $`e_j`$ and the trajectories in the diagram of $`e_1`$ versus $`\varpi _1\varpi _2`$ are qualitatively indistinguishable from those showing regular secular evolution (for at least $`2\times 10^5\mathrm{yr}`$).

In Figures 5 and 6 we show examples of the cases with regular secular evolution. Figure 5 shows the variations in $`a_j`$ and $`e_j`$ for the $`\mathrm{sin}i=1`$ HD 12661 system with $`P_2/P_1=0.99\times 11/2`$. The solid and dotted lines are from the direct integrations starting at $`T_{\mathrm{peri},1}`$ and $`T_{\mathrm{peri},2}`$, respectively, and the dashed lines are from the octupole-level secular perturbation theory. The trajectories of this system in the phase-space diagrams of $`e_1`$ versus $`\varpi _1\varpi _2`$ and $`e_2`$ versus $`\varpi _1\varpi _2`$ are the trajectories marked by the squares (which denote the current parameters of the HD 12661 system) in Figure 6. The direct integration results are not sensitive to the starting epoch. The amplitudes of eccentricity oscillations and the trajectories in the diagram of $`e_1`$ (or $`e_2`$) versus $`\varpi _1\varpi _2`$ predicted by the octupole theory are in reasonably good agreement with those from direct integrations, even though $`\alpha `$($`=0.323`$) is quite large. The main error of the octupole theory is in the period of eccentricity oscillations, with the predicted period ($`2.1\times 10^4\mathrm{yr}`$) about $`75\%`$ longer than that found in the direct integrations ($`1.2\times 10^4\mathrm{yr}`$).

Figures 5 and 6 show some interesting properties of the HD 12661 system that are shared by all of the cases studied above with different initial $`P_2`$, including those affected by the close proximity to a mean-motion commensurability. The HD 12661 system is in a secular resonance with $`\varpi _1\varpi _2`$ librating about $`180^{}`$. This means that the lines of apsides of the two orbits are on average anti-aligned. However, the amplitude of libration of $`\varpi _1\varpi _2`$ is large, and the orbital eccentricities of both planets exhibit large-amplitude variations. For the direct integrations shown in Figures 5 and 6, the amplitude of libration of $`\varpi _1\varpi _2`$ is $`56^{}`$, $`e_1`$ varies between $`0.09`$ and $`0.37`$, and $`e_2`$ varies between $`0.17`$ and $`0.37`$. HD 12661 is the first extrasolar planetary system found to have $`\varpi _1\varpi _2`$ librating about $`180^{}`$. The outer two planets of the $`\upsilon `$ And system (Rivera & Lissauer 2000; Lissauer & Rivera 2001; Chiang et al. 2001) and the two planets about GJ 876 (Laughlin & Chambers 2001; Lee & Peale 2002) are also likely in the secular resonance involving $`\varpi _1\varpi _2`$, but in both of these cases $`\varpi _1\varpi _2`$ librates about $`0^{}`$. (In the case of GJ 876, the secular resonance is associated with the simultaneous librations of both lowest order mean-motion resonance variables at the 2:1 commensurability.)

The libration of $`\varpi _1\varpi _2`$ and the large-amplitude variations of both $`e_j`$ in the HD 12661 system were anticipated by the analytic results derived in § 3.3. As we discussed in § 3.3, the octupole theory predicts that large libration islands and large-amplitude variations of both $`e_j`$ are likely if $`\lambda =L_1/L_2\lambda _{\mathrm{crit}}=2\gamma ^2/(53\gamma ^2)`$. The HD 12661 system has $`\lambda 0.83`$, which is almost identical to $`\lambda _{\mathrm{crit}}0.82`$ for the dimensionless total angular momentum $`\gamma 0.96`$ of this system. In addition to the trajectories of the $`\mathrm{sin}i=1`$ HD 12661 system with $`P_2/P_1=0.99\times 11/2`$, Figure 6 also shows the trajectories of systems with the same masses, initial semimajor axes, and total angular momentum. The direct-integration results (left panels) are in reasonably good agreement with the octupole results (right panels) and confirm that the phase space is dominated by large libration islands about $`\varpi _1\varpi _2=0^{}`$ and $`180^{}`$. The only discrepancy is a narrow region of circulation represented by a single trajectory in the direct-integration results. Thus libration of $`\varpi _1\varpi _2`$ is almost the only possibility. We can see in Figure 6 that $`e_1`$ decreases from near its maximum possible value to near zero and $`e_2`$ increases from near zero to near its maximum possible value at $`\varpi _1\varpi _290^{}`$, and vice versa at $`\varpi _1\varpi _2270^{}`$, again in agreement with the analysis in § 3.3.

Kiseleva-Eggleton et al. (2002) have recently reported that the $`\mathrm{sin}i=1`$ HD 12661 system with the best-fit orbital parameters in Table 1 is regular based on the MEGNO (Mean Exponential Growth of Nearby Orbits) technique. However, they assumed that the best-fit orbital parameters are in astrocentric coordinates and showed results from only one (unspecified) starting epoch. We have repeated the direct integrations shown in Figure 4, but with the best-fit orbital parameters assumed to be in astrocentric coordinates. Because of the larger fluctuations in the astrocentric $`a_2`$, the average ratio of the orbital periods is more sensitive to the starting epoch. The integration starting at $`T_{\mathrm{peri},1}`$ turns out to be sufficiently far from the 11:2 commensurability that it shows regular secular evolution similar to those in Figure 5. The integration starting at $`T_{\mathrm{peri},2}`$ is qualitatively similar to the integration in Jacobi coordinates starting at the same epoch (dotted lines in Figure 4) and is thus either regular or very weakly chaotic. These results again demonstrate the importance of interpreting the orbital parameters from multiple Kepler fits as orbital parameters in Jacobi coordinates.

The best-fit orbital parameters of the HD 12661 planets in Table 1 were provided by D. A. Fischer (2002, private communication), and some of them differ slightly from those published later in Fischer et al. (2003). Adopting the parameters in Fischer et al. (2003) does not qualitatively change the results presented here, although the smaller $`\varpi _1\varpi _2`$ for the parameters in Fischer et al. (2003) does lead to slightly larger amplitude of libration of $`\varpi _1\varpi _2`$ and slightly larger variations in $`e_j`$.

After this paper was completed, we learned of contemporaneous work on the HD 12661 system by Goździewski (2003) and Goździewski & Maciejewski (2003). Goździewski (2003) has examined a wide region of the parameter space near the HD 12661 best-fit solution in Table 1, using astrocentric coordinates and a combination of the MEGNO technique and direct integrations. His numerical results also show the proximity of the best-fit solution to the 11:2 commensurability and the robustness of the libration of $`\varpi _1\varpi _2`$. Goździewski & Maciejewski (2003) have subsequently reported a best-fit solution to the data set published by Fischer et al. (2003) with $`P_21660`$days or $`P_2/P_16.3`$. This illustrates the uncertainties mentioned above in the orbital parameters of the outer planet due to the time span of the available observations being comparable to $`P_2`$. Their best-fit solution also shows libration of $`\varpi _1\varpi _2`$ about $`180^{}`$ and large-amplitude variations of both $`e_j`$.

Mayor et al. (2003) have reported the detection of two planets around the star HD 83443. There is considerable doubt about the existence of the outer planet in this system (Butler et al. 2002), but it is interesting to note that the system as reported by Mayor et al. shares some of the properties of the HD 12661 system. The best-fit orbital periods of the HD 83443 planets are nearly in the ratio 10:1, and direct integrations of this system (neglecting the effects of tides and general relativistic precession) show that its evolution is also affected by the close proximity to the mean-motion commensurability. The HD 83443 system also shows large-amplitude variations of both eccentricities, although its $`\varpi _1\varpi _2`$ is circulating. The value of $`\lambda `$($`0.96`$) of the HD 83443 system is not as close to its value of $`\lambda _{\mathrm{crit}}`$($`0.79`$) as in the HD 12661 system, and there is a larger region of circulation in the phase space of systems with the same masses, initial semimajor axes, and total angular momentum as the HD 83443 system. It should be noted that the inner planet of the HD 83443 system is sufficiently close to the star ($`P3`$days) that tides and general relativistic precession can significantly affect the dynamics of the system (Wu & Goldreich 2002).

### 4.3 Small Mutual Inclination

To study the effects of a small mutual inclination angle $`i_{\mathrm{mu}}`$ between the orbits, we have repeated the direct integrations of the HD 168443 and HD 12661 systems shown in Figures 1b and 5 with initial $`i_{\mathrm{mu}}=1^{}`$ or $`2^{}`$. We assumed that the intersection of the orbital planes is initially along the line of sight, so that, initially, $`\mathrm{sin}i_1=\mathrm{sin}i_2=1`$ and the mutual inclination is the difference in the longitudes of ascending node referenced to the plane of the sky: $`i_{\mathrm{mu}}=\mathrm{\Omega }_2\mathrm{\Omega }_1`$. The results are nearly identical to those with coplanar orbits, confirming the analysis in § 3.1 based on an expansion of the more general form of the octupole theory in powers of $`i_{\mathrm{mu}}`$.

### 4.4 Stability

The applicability of the octupole theory is limited by its use of an expansion in $`r_1/r_2`$ and orbit averaging. In particular, the orbit averaging eliminates the effects of mean-motion commensurabilities and the possible development of instabilities. The direct integrations reported in § 4.1 show that the coplanar HD 168443 system is stable for $`\mathrm{sin}i0.4`$, and Marcy et al. (2001) have found that this system is stable for even smaller $`\mathrm{sin}i`$. The stability of the HD 168443 system is consistent with the empirical stability criteria for hierarchical triple systems derived most recently by Eggleton & Kiseleva (1995) and Mardling & Aarseth (2001). For example, the $`\mathrm{sin}i=1`$ HD 168443 system has $`\alpha =a_1/a_2=0.102`$, while the stability criterion of Eggleton & Kiseleva requires initial $`a_2(1e_2)/[a_1(1+e_1)]1.65`$ or $`\alpha 0.316`$ and that of Mardling & Aarseth requires initial $`a_2(1e_2)/a_13.17`$ or $`\alpha 0.252`$. The criterion of Eggleton & Kiseleva was tested over a wide range of mass ratios (including ones similar to those of the HD 168443 system) and is expected to be reliable to about $`20\%`$, while the criterion of Mardling & Aarseth does not have explicit dependence on $`m_1/m_0`$ and may not be very accurate for $`m_1/m_0`$ significantly different from unity. For a more detailed comparison with the above stability criteria, we have performed a set of direct integrations with different initial $`a_2`$: $`a_2=1.299`$$`0.738\mathrm{AU}`$ or $`\alpha =0.227`$$`0.400`$. The other initial conditions are identical to those of the $`\mathrm{sin}i=1`$ calculation shown as the solid lines in Figure 1b. The evolutions of $`\alpha `$ are shown in Figure 7. Except for the case with initial $`\alpha =0.333`$, which is apparently stabilized by the 5:1 mean-motion commensurability, there is a clear boundary between stable and unstable configurations at initial $`\alpha 0.30`$. We do not expect this stability boundary to move down significantly if we increase the integration length, since $`10^5\mathrm{yr}`$ span many eccentricity oscillation cycles (if the orbits are stable) and, as we can see in Figure 7, instabilities usually manifest themselves on much shorter timescales. The stability criterion of Eggleton & Kiseleva is a few percent larger than the stability boundary determined here, consistent with the $`20\%`$ accuracy of the former. Even the stability criterion of Mardling & Aarseth is only $`16\%`$ smaller than the stability boundary determined here (see, however, the next paragraph). Finally, we note that some of the configurations in Figure 7 that are stable for at least $`10^5\mathrm{yr}`$ appear to be chaotic (e.g., the case with initial $`\alpha =1/4`$). As in the case of the HD 12661 system, the chaos is due to the close proximity to a high-order mean-motion commensurability (8:1 for the case with initial $`\alpha =1/4`$).

In § 4.2 we have considered the coplanar HD 12661 system with $`\mathrm{sin}i=1`$ (i.e., minimum planetary masses) only. Unlike the HD 168443 system, which has $`\alpha =a_1/a_20.1`$ and is stable for any reasonable value of $`\mathrm{sin}i`$, the HD 12661 system has $`\alpha 0.32`$ (which is close to the stability boundary for the systems shown in Fig. 7 with the planetary masses of the $`\mathrm{sin}i=1`$ HD 168443 system) and could be unstable for moderately small $`\mathrm{sin}i`$. To study the effects of $`\mathrm{sin}i`$, we have performed two sets of direct integrations with different $`\mathrm{sin}i`$. The initial values of $`P_j`$, $`K_j`$, $`e_j`$, $`\omega _j`$, and $`T_{\mathrm{peri},j}`$ are the same as those of the HD 12661 system with $`P_2/P_1=0.99\times 11/2`$ shown in Figure 5, and the planetary masses and initial semimajor axes were derived for the assumed $`\mathrm{sin}i`$. The evolutions of $`\alpha `$ are shown in Figure 8 for $`\mathrm{sin}i=0.1`$, $`0.2,\mathrm{},1`$. The left and right panels are from the direct integrations starting at $`T_{\mathrm{peri},1}`$ and $`T_{\mathrm{peri},2}`$, respectively. The dashed line in each panel indicates the location of the 11:2 mean-motion commensurability. The influence of the 11:2 commensurability increases with decreasing $`\mathrm{sin}i`$ (which is expected since the width of a resonance increases with increasing planetary masses), but the cases with $`\mathrm{sin}i>0.3`$ are stable (for at least $`2\times 10^5\mathrm{yr}`$). All of the cases with $`\mathrm{sin}i0.3`$ are unstable. Are the stability criteria of Eggleton & Kiseleva (1995) and Mardling & Aarseth (2001) consistent with this transition from stable to unstable configurations at $`\mathrm{sin}i0.3`$? For the configurations shown in Figure 8, which have the same initial $`P_1`$ and $`P_2`$, the initial $`\alpha `$ decreases slightly from $`0.323`$ to $`0.322`$ as $`\mathrm{sin}i`$ decreases from $`1`$ to $`0.1`$ (and $`m_1`$ and $`m_2`$ increase; see eqs.  and ). For comparison, the stability boundary predicted by the criterion of Eggleton & Kiseleva decreases from initial $`\alpha =0.444`$ to $`0.348`$ as $`\mathrm{sin}i`$ decreases from $`1`$ to $`0.1`$, with initial $`\alpha =0.397`$ at $`\mathrm{sin}i=0.3`$. This stability boundary is consistent with the transition at $`\mathrm{sin}i0.3`$ if we simply move it down by $`23\%`$, which is compatible with the $`20\%`$ accuracy of the stability criterion of Eggleton & Kiseleva. On the other hand, the stability boundary predicted by the criterion of Mardling & Aarseth changes only slightly from initial $`\alpha =0.254`$ to $`0.253`$ as $`\mathrm{sin}i`$ decreases from $`1`$ to $`0.1`$. Thus, this stability boundary is not consistent with the transition at $`\mathrm{sin}i0.3`$, unless we allow adjustments that change with $`\mathrm{sin}i`$. This is not surprising since the stability criterion of Mardling & Aarseth does not have explicit dependence on $`m_1/m_0`$ and is not expected to be very accurate for $`m_1/m_0`$ significantly different from unity.

## 5 CONCLUSIONS

We have investigated the dynamical evolution of coplanar hierarchical two-planet systems where the ratio of the orbital semimajor axes $`\alpha =a_1/a_2`$ is small. Hierarchical two-planet systems are common among the known multiple planet systems and are likely to be common in the overall population of extrasolar planetary systems. We began by showing that the orbital parameters obtained from a two (or more generally multiple) Kepler fit to the radial velocity variations of a host star are best interpreted as Jacobi coordinates and that these coordinates should be used in any analyses of hierarchical (and possibly other types of) planetary systems. We used the HD 168443 system as an example to show that the use of astrocentric coordinates can introduce significant high-frequency variations in orbital elements that should be nearly constant on orbital timescales, leading to erroneous sensitivity of the evolution of the orbital elements on the starting epoch.

An approximate theory that can be applied to hierarchical two-planet systems with a wide range of masses, orbital eccentricities, and inclinations is the octupole-level secular perturbation theory, which was developed for general hierarchical triple systems by Marchal (1990), Krymolowski & Mazeh (1999), and Ford et al. (2000). The octupole theory is based on an expansion to order $`\alpha ^3`$ and orbit-averaging. We showed that the octupole approximation reduces the coplanar problem to one degree of freedom, with $`e_1`$ (or $`e_2`$) and $`\varpi _1\varpi _2`$ as the relevant phase-space variables. An analysis of the octupole equations yielded the following results: (1) the scaling properties of the equations imply that the amplitudes of eccentricity oscillations and the trajectory in the phase-space diagram of $`e_1`$ (or $`e_2`$) versus $`\varpi _1\varpi _2`$ should be independent of the inclination $`i`$ of both orbits from the plane of the sky (as long as $`m_1,m_2m_0`$) and that the period of eccentricity oscillations should be proportional to $`\mathrm{sin}i`$; (2) the fixed points in the phase-space diagram must be at $`\varpi _1\varpi _2=0^{}`$ or $`180^{}`$; and (3) if the ratio of the maximum orbital angular momenta, $`\lambda =L_1/L_2(m_1/m_2)\alpha ^{1/2}`$, for given semimajor axes is approximately equal to $`\lambda _{\mathrm{crit}}=2\gamma ^2/(53\gamma ^2)`$, where $`\gamma =(G_1+G_2)/(L_1+L_2)`$ is the dimensionless total angular momentum, then libration of $`\varpi _1\varpi _2`$ is almost certain, with possibly large amplitude variations of both eccentricities.

We used both the octupole-level secular perturbation theory and direct numerical orbit integrations to study the dynamical evolution of the HD 168443 and HD 12661 systems and their variants. The HD 168443 system is not in a secular resonance and its $`\varpi _1\varpi _2`$ circulates. For the family of systems with the same masses, initial semimajor axes, and total angular momentum as the $`\mathrm{sin}i=1`$ HD 168443 system, there are two relatively small libration islands in the phase-space diagram of $`e_1`$ (or $`e_2`$) versus $`\varpi _1\varpi _2`$, one about a fixed point at $`\varpi _1\varpi _2=0^{}`$ and the other about a fixed point at $`\varpi _1\varpi _2=180^{}`$, and the trajectory of the HD 168443 system is far from both. In all cases, the octupole results are in excellent agreement with the direct-integration results. Direct integrations of the HD 168443 system with $`\mathrm{sin}i=0.4`$ also confirmed the analytic octupole results on the effects of $`\mathrm{sin}i`$. To demonstrate how the octupole theory can be used to explore the parameter space rapidly, we used the octupole theory to determine the number and positions of fixed points as a function of the total angular momentum for systems with the same masses and semimajor axes as the $`\mathrm{sin}i=1`$ HD 168443 system. For sufficiently low total angular momentum, the existence of fixed point(s) in addition to those appropriate for the total angular momentum of the HD 168443 system leads to systems with unusual behaviors such as $`e_1`$ increasing to unity in a finite time or libration of $`\varpi _1\varpi _2`$ about $`0^{}`$ with $`e_1`$ very close to 1. These results were again confirmed by direct integrations. Direct integrations of systems similar to the $`\mathrm{sin}i=1`$ HD 168443 system, but with different initial $`a_2`$, showed that these systems with relatively massive planets \[$`(m_1+m_2)/m_00.024`$\] are unstable if the initial $`\alpha 0.30`$.

The HD 12661 system is the first extrasolar planetary system found to have $`\varpi _1\varpi _2`$ librating about $`180^{}`$. The secular resonance means that the lines of apsides of the two orbits are on average anti-aligned, although the amplitude of libration of $`\varpi _1\varpi _2`$ is large. The libration of $`\varpi _1\varpi _2`$ and the large-amplitude variations of both orbital eccentricities in the HD 12661 system, which has $`\lambda 0.83`$ and $`\lambda _{\mathrm{crit}}0.82`$, are consistent with the analytic results on systems with $`\lambda \lambda _{\mathrm{crit}}`$. Direct integrations showed that the evolution of the $`\mathrm{sin}i=1`$ HD 12661 system with the best-fit orbital parameters ($`P_2/P_1=0.9975\times 11/2`$) is affected by the close proximity to the 11:2 mean-motion commensurability, but that small changes in the orbital period of the outer planet within the uncertainty can result in configurations that are not affected by mean-motion commensurabilities. For the family of systems with the same masses, initial semimajor axes, and total angular momentum as the $`\mathrm{sin}i=1`$ HD 12661 system with $`P_2/P_1=0.99\times 11/2`$, which are not affected by mean-motion commensurabilities, the octupole results are in reasonably good agreement with the direct-integration results, even though $`\alpha `$($`0.32`$) is quite large. The phase space of this family of systems with $`\lambda \lambda _{\mathrm{crit}}`$ is dominated by large libration islands about $`\varpi _1\varpi _2=0^{}`$ and $`180^{}`$, which again confirmed the analytic results on systems with $`\lambda \lambda _{\mathrm{crit}}`$, and libration in the secular resonance is almost the only possibility. Finally, we showed that the HD 12661 system is unstable if $`\mathrm{sin}i0.3`$ \[or $`(m_1+m_2)/m_00.012`$\].

The comparisons with direct numerical orbit integrations established that the octupole-level secular perturbation theory is highly accurate for coplanar hierarchical two-planet systems with $`\alpha 0.1`$ and reasonably accurate even for systems with $`\alpha `$ as large as $`1/3`$, provided that $`\alpha `$ is not too close to a significant mean-motion commensurability or above the stability boundary. We have focused our study in this paper on coplanar systems in which the star and the planets can be treated as point masses interacting via Newtonian gravity, where a small mutual inclination between the orbits is shown to have little effect on the secular evolution. The more general form of the octupole theory is applicable to systems with arbitrary mutual inclinations. If the inner or both planets are very close to the star, tidal interactions and general relativistic precession can be important (see, e.g., Wu & Goldreich 2002). Nagasawa et al. (2003) have also discussed the effects of additional orbital precessions due to interactions with the protoplanetary disk during the epoch of disk depletion. It is straightforward to modify the octupole theory to include any of these effects (see, e.g., Blaes et al. 2002 for relativistic precession).

We thank Debra Fischer for furnishing the orbital parameters of the HD 12661 system before publication. We also thank O. Blaes, G. W. Marcy, N. Murray, E. J. Rivera, and Y. Wu for informative discussions and an anonymous referee for comments that greatly improved the presentation. This research was supported in part by NASA grants NAG5 11666 and NAG5 12087.