{'fraction_non_alphanumeric': 0.04637601092644599, 'fraction_numerical': 0.014747632859137855, 'mean_word_length': 4.130215714060778, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 32, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Theory of Secular Chaos and Mercury's Orbit

We study the chaotic orbital evolution of planetary systems, focusing on secular (i.e., orbit-averaged) interactions, because these often dominate on long timescales. We first focus on the evolution of a test particle that is forced by multiple massive planets. To linear order in eccentricity and inclination, its orbit precesses with constant frequencies. But nonlinearities can shift the frequencies into and out of secular resonance with the planets' eigenfrequencies, or with linear combinations of those frequencies. The overlap of these nonlinear secular resonances drive secular chaos in planetary systems. We quantify the resulting dynamics for the first time by calculating the locations and widths of nonlinear secular resonances. When results from both analytical calculations and numerical integrations are displayed together in a newly developed"map of the mean momenta"(MMM), the agreement is excellent. This map is particularly revealing for non-coplanar planetary systems and demonstrates graphically that chaos emerges from overlapping secular resonances. We then apply this newfound understanding to Mercury. Previous numerical simulations have established that Mercury's orbit is chaotic, and that Mercury might even collide with Venus or the Sun. We show that Mercury's chaos is primarily caused by the overlap between resonances that are combinations of four modes, the Jupiter-dominated eccentricity mode, the Venus-dominated inclination mode and Mercury's free eccentricity and inclination. Numerical integration of the Solar system confirms that a slew of these resonant angles alternately librate and circulate. We are able to calculate the threshold for Mercury to become chaotic: Jupiter and Venus must have eccentricity and inclination of a few percent. Mercury appears to be perched on the threshold for chaos.


INTRODUCTION
The question of the stability of planetary orbits in the Solar system has a long history, and has attracted the attention of some of the greatest scientists, including Newton, Laplace, Lagrange, Gauss, Poincaré, Kolmogorov, and Arnol'd. Newton thought that interplanetary perturbations are eventually destabilizing, and that divine intervention is required to restore the planets' orbits to their rightful places (Laskar 1996). Yet it is only over the last twenty years that the stability of the Solar system has been definitively settled, with the aid of computer simulations (Sussman & Wisdom 1988;Laskar 1989;Quinn et al. 1991;Wisdom & Holman 1991;Lecar et al. 2001;Laskar & Gastineau 2009). We now know that Newton was not far off: the Solar system is marginally stable: it is unstable, but on a timescale comparable to its age. In the inner Solar system, the planets' eccentricities chaotically diffuse on a billion-year timescale, with the two lightest planets, Mercury and Mars, experiencing particularly large variations. In fact, Mercury has roughly a 1% chance of colliding with Venus or the Sun within the next five billion years (Laskar & Gastineau 2009). By comparison, the giant planets in the outer Solar system are well-spaced, and their orbital elements undergo largely quasiperiodic variations, exhibiting chaotic diffusion only on extremely long timescales (Laskar 1996;Murray & Holman 1999).
The chaotic dynamics in the inner Solar system is pri-marily due to secular interactions (Laskar 2008). In general, interplanetary interactions can be decomposed into secular ones and MMR's (mean motion resonances-not to be confused with secular resonances). Secular interactions result from orbit-averaging the equations of motion.
Since averaging a Keplerian orbit produces an elliptical ring, secular evolution can be thought of as interactions between elliptical rings. Secular timescales are longthey are longer than the orbital time by at least the ratio of the star's mass to that of a planet. By contrast, interactions driven by MMR's depend on orbital phase, and typically occur on the orbital timescale, or longer if some of the planets' orbital periods are close to integer ratios. Intuitively, one would expect that the dynamics on long timescales can be treated by averaging over the fast orbital phase-i.e., that they are secular in nature. This is true in the inner Solar system. It is also true more generally for well-spaced planets that do not happen to lie near mean motion resonances. 3 Linear secular theory has been understood for hundreds of years, dating back to the famous solution of Laplace and Lagrange (see Murray & Dermott 2000). To linear order in the planets' eccentricities and inclinations, secular theory reduces to a simple eigenvalue problem, with two eigenmodes per planet-one for the eccentricity degree of freedom, and one for the inclination. Each eigenmode has a constant amplitude and a longitude that precesses uniformly in time. But linear secular theory is clearly incapable of describing the chaotic orbits of the 3 In the outer Solar system the dynamics is not mainly secular because the giant planets lie near a number of MMR's, such as the 5:2 between Jupiter and Saturn (the "Great Inequality"), and the 2:1 between Uranus and Neptune.

Solar system.
Despite the importance of secular chaos, there has been surprisingly little theoretical understanding of it (see, e.g., the review of Solar system chaos by Lecar et al. 2001). By contrast, chaos due to MMR's is wellunderstood, and accounts for the Kirkwood gaps in the asteroid belt (Wisdom 1983), and for the very weak chaos of the outer Solar system planets, which is due to 3-body MMR's (Murray & Holman 1999). In all cases that have been studied in the Solar system, chaos is caused by overlapping resonances (Chirikov 1979;Lecar et al. 2001). In linear secular theory, there can be secular resonances. And it is generally supposed that the chaos in the inner Solar system is caused by overlapping secular resonances. Yet thus far there has been little quantitative calculation. To our knowledge, the only previous theoretical work towards calculating secular chaos was by Sidlichovsky (1990), who considered the coplanar case, as we describe below ( §3).
Numerical attempts to identify the mechanism of chaos in the inner Solar system were made by Laskar (1990Laskar ( , 1992 and Sussman & Wisdom (1992). These authors found that the angle associated with the (secular) frequency (g mercury − g jupiter ) − (s mercury − s venus ) alternately librated and circulated in their simulations, where g is the apsidal precession rate, and s is the nodal precession rate (or, to be more precise, g and s here refer to the frequencies of the normal mode that is dominated by the corresponding planet). Laskar (1992) also found that two angles associated with Earth and Mars, corresponding to 2(g mars − g earth ) − (s mars − s earth ) and (g mars −g earth )−(s mars −s earth ), alternately librated, and conjectured that the overlap of those secular resonances was responsible for chaos. But, as Sussman & Wisdom (1992) note, Laskar's conjecture is not fully convincing, because there are too many unrelated angles that alternately circulate and librate, and it is not clear which are dynamically important. Furthermore, only one librating angle has been identified for Mercury. Yet chaos requires the overlap of at least two resonances, so why is Mercury chaotic (Lecar et al. 2001)? Without a theory for secular chaos, the dynamics remain obscure. For example, why does instability in the Solar system occur at such low values of eccentricity and inclination (∼ few percent)? What sets the timescale of the chaos? Can secular chaos shape the architecture of the inner Solar system (Laskar 1996)? And can it shape the architecture of extrasolar planetary systems (Wu & Lithwick 2010)? Without a theory for secular chaos, we will be forever at the mercy of computer simulations.
In this paper, we construct the theory for secular chaos of a test particle, and then apply the theory to Mercury. In §2, we present the test particle's equations of motion. In §3, we describe the coplanar solution, and in §4 we generalize to the case when bodies have non-zero inclinations. In §5, we apply the theory to N-body simulations of the real Mercury. We conclude in §6.

SECULAR EQUATIONS OF MOTION
We focus on the secular evolution of a massless test particle that is orbiting a star in the presence of multiple massive planets, assuming the planets' orbits are known.
The particle has six orbital elements, {a, e, i, λ, , Ω}, using standard notation (Murray & Dermott 2000). In secular theory, one averages over λ. As a consequence, a is a constant of motion, leaving only four orbital elements to be considered. The equations of motion for the particle's eccentricity and longitude of periapse (e and ) are given by Hamilton's equations for the Poincaré canonical variables Γ ≡ GM a 1 − √ 1 − e 2 and γ ≡ − (Murray & Dermott 2000). Since a is constant, it is simpler to choose the canonical momentum to be ∝ Γ/ GM a, so we introduce the momentum and its conjugate co-ordinate to be . Although this is a non-canonical transformation from Poincaré's variables, if we simultaneously re-scale the energy by defining as the Hamiltonian where E is the particle's energy per unit mass 4 then the equations of motion are Hamilton's equations, Therefore we may consider (p e , ) to be canonically conjugate. Similarly, for the inclination and longitude of node (i, Ω), we take the canonical variables to be related to the corresponding Poincaré variables in the same way by defining and taking its conjugate co-ordinate to be Ω. The equations of motion for (p e , , p i , Ω) are Hamilton's equations generated by the scaled H. These equations are exact as long as a is constant, which is the case for secular interactions. An alternative formulation of the equations of motion will also prove useful. For an arbitrary Hamiltonian H(p, q), one may define the complex canonical variable 5 Z ≡ √ pe iq . As may be directly verified, the equation of motion for Z is then dZ/dt = i∂H/∂Z * where H(Z, Z * ) = H(p, q). This complex equation of motion simultaneously encodes both of Hamilton's equations. Since our real canonical variables are (p e , , p i , Ω), we introduce the complex ones which, to leading order in e and i, are the usual complex eccentricity and inclination. 6 The equations of motion 4 The test particle's energy per unit mass E is given by equation (A1) in Appendix A for the case of a single external planet. 5 Our definition of the complex canonical variable differs from Ogilvie (2007) by a minus sign, and hence our Hamilton's equation also differs by a minus sign. 6 The symbol e denotes both the eccentricity and the exponential (Euler's constant), and the symbol i denotes both the inclina-for (z, ζ) are Throughout this paper, we freely switch between the set of real canonical variables and the set of complex ones, RCV : (p e , ; p i , Ω) (12) CCV : (z, ζ) Although Hamilton's equations of motion for the RCV and CCV are exact, we take an approximate form for the Hamiltonian by expanding H to fourth order in e and i, and to leading order in the ratios of semimajor axes, keeping only secular terms. The relevant terms are listed in Table 1 in Appendix A. Throughout the bulk of this paper, we focus on the c 1 , c 2 , c 4 , c 11 , c 12 , c 14 , and c 15 terms in that table. As we show in §4.5, the remaining terms in the table have a small effect in the parameter regime we focus on.

COPLANAR JUPITER AND SATURN
In this section, we consider the evolution of a test particle in the presence of two massive exterior planets ("Jupiter" and "Saturn"), with all bodies having zero inclinations. We assume that the massive planets' orbits are given by their linear Laplace-Lagrange solution, and evolve the test particle's equations of motion to leading nonlinear order (terms listed in Table 1 in Appendix A). This case was worked out by Sidlichovsky (1990). We describe it here in some detail because it sets the stage for the considerably more complicated case with non-zero inclinations ( §4).

Jupiter Only
We shall solve the coplanar case with a sequence of increasingly complicated sub-cases. Consider first the linear secular evolution of a test particle perturbed by a circular Jupiter. From Table 1 in the Appendix, the particle's Hamiltonian is where the constant is the particle's linear free precession rate induced by Jupiter (eq. [A4]); m J is Jupiter's mass, and α is the ratio of the particle's semimajor axis to Jupiter's. The equation of motion is dz/dt = i∂H/∂z * = iγz, with solution z = const × e iγt , an orbit with constant eccentricity that precesses at frequency γ.
For the second sub-case, we consider the test particle's linear evolution when Jupiter is assigned a constant eccentricity e J and precession rate g J , in which case the particle's Hamiltonian is tion and the imaginary unit. There should be no confusion because for the remainder of this paper we use as our dynamical variables either the RCV or CCV (eqs. [12]-[13]) in lieu of e and i.
where J = 5 4 αe J ; we drop nonlinear terms (i.e. the fourth-order terms in Table 1). Of course, if Jupiter were the only massive planet in the system, it would not precess (g J = 0). But we consider a finite g J in anticipation of what happens when Saturn is included.
Hamilton's equation is The solution is a sum of free and forced eccentricities: The forced solution diverges at secular resonance (γ = g J ). But this divergence is not physical. It is a consequence of dropping nonlinear terms. With the nonlinear term included, the test particle's Hamiltonian becomes 7 with equation of motion The nonlinear term reduces the free frequency from γ to g ≡ γ(1 − |z| 2 /2). We call g the nonlinear free frequency. It can also be introduced by first re-writing the Hamiltonian in terms of the real canonical variables p e and and defining Because g is a function of eccentricity, secular resonance can occur if the eccentricity is chosen so that g ≈ g J . This is a nonlinear secular resonance. Solutions of equation (21) are shown in Figure 1 for two cases, γ < g J and γ > g J . When γ < g J (left panel), the free frequency at small |z| is less than the forcing frequency, g| |z|→0 ≈ γ < g J . Increasing |z| decreases g, moving it further away from resonance. Hence the particle's evolution is similar to the linear case, but with a smaller precession rate. When γ > g J (right panel), the precession rate at small values of |z| exceeds g J . Increasing |z| causes the precession rate g to decrease until it nearly matches g J . A resonant island appears, within which the angle − g J t librates. Increasing |z| even further forces the free frequency to be less than g J , taking the particle out of secular resonance.
To calculate the width of the nonlinear secular resonance, we first write Hamiltonian (20) in terms of the real canonical variables p e and (eq. [8]), and then convert the Hamiltonian to a time-independent one by making the canonical transformation to (p e , − ), where  (21). Equivalently, they are level curves of Hamiltonian (24). The y-axis is approximately e 2 , and the x-axis is the phase, modulo 2π. Nonlinear secular resonance occurs for the case g J < γ.
yielding the new Hamiltonian This form for a Hamiltonian has been called the "second fundamental model for resonance," the first being the pendulum. Its properties have been extensively catalogued because it also applies to first-order mean motion resonances, for which the variables have a different interpretation (Henrard & Lemaitre 1983;Murray & Dermott 2000). The numerical integrations shown in Figure  1 trace level curves of this Hamiltonian. For parameters such that the resonant island both exists and is sufficiently far from |z| = 0 (as in the right panel of Fig.  1), one may approximate the √ p e multiplying the cosine term as constant, in which case the Hamiltonian is that of the pendulum. The center of the resonant island is where (d/dt) − = 0, i.e. at in this approximation; this is equivalent to g = g J . The width of the island may be found by first completing the square in the "kinetic" part of the Hamiltonian, i.e. setting −p 2 e /4 + ∆ · p e = − 1 4 (p e − p e * ) 2 +const. Since H − is constant, the half-width of the island is given by 3.2. Jupiter and Saturn We now add in the effect of a second coplanar planet, Saturn, and assume that Jupiter's and Saturn's evolution is described by their linear Laplace-Lagrange solution.
In that solution, Jupiter and Saturn participate in two normal modes, which we call the Jupiter-dominated and Saturn-dominated modes. We denote the eigenfrequencies of these two modes g J and g S . Jupiter's (complex) eccentricity is a sum of two terms, one for each mode, and may be written as e J,J e ig J t + e J,S e ig S t . Similarly, Saturn's eccentricity has one term ∝ e ig J t and another ∝ e ig S t . When all four terms are included, one arrives [27]): Section taken at times when e i(g J −g S )t = 1; the y-axis is approximately e 2 ; the x-axis is the phase, modulo 2π. In the left panel, the chosen parameters yield non-overlapping separatrices, and very little chaos is seen. In the right panel g J is slightly larger, yielding overlapping separatrices and a sea of chaos.
at the following form for the test particle's Hamiltonian (see eq. [20]): where γ is now re-interpreted to represent the test particle's precession rate due to both Jupiter and Saturn. The J term is due to the Jupiter-dominated mode, with J now a weighted sum of both Saturn's and Jupiter's eccentricity within that mode; similarly, S is for the Saturndominated mode. For the purposes of this section, rather than solving the linear Laplace-Lagrange equations, we consider the parameters γ, J , and S to be adjustable constants.
Hamiltonian (27) would also describe the test particle's evolution if Jupiter's and Saturn's eccentricities and precession rates were enforced by hand to be e J , e S , g J , and g S , with J = 5 4 α J e J and S = 5 4 α S e S . This reinterpretation of Hamiltonian (27), while unphysical as far as Jupiter and Saturn are concerned, can be helpful when considering the test particle's evolution under their influence. Figure 2 shows results of numerical integrations of Hamilton's equation, plotted as surfaces of section. Changing to real canonical variables (eq. [8]), we see that Hamiltonian (27) has two cosine terms. The J cosine term acting alone would yield a resonant island as long as γ g J . Similarly, the S term would yield an island if γ g S . The location and width of the islands are quantified in Figure 1. When acting together, there may be two resonant islands. The left panel of Figure 2 shows a case when the parameters have been chosen to yield two non-overlapping islands. The result is mostly regular motion. The right panel shows what happens when g J is increased, so that Jupiter's ∆ (eq. [19]) is reduced sufficiently that the islands overlap: the overlapping islands break up into a sea of chaos. This the well-known Chirikov resonance-overlap criterion for chaos (Chirikov 1979). From the widths and locations of the resonances as displayed in Figures 1 and 2, one deduces that the criterion for chaos is 2|g J − g S |/γ 4 1/2 J (2∆) 1/4 + 4 1/2 S (2∆ S ) 1/4 , where ∆ S ≡ (γ − g S )/γ (Sidlichovsky 1990).
It is instructive to consider the form of the surface of section shown in Figure 2 in more detail. This reasoning will also be helpful when we include a second degree of freedom (i.e. inclination) below. When the motion of z is regular, it can be written as a Fourier sum of terms with frequencies equal to the three fundamental frequencies of the problem (i.e. g, g J , g S , where g is the nonlinear free frequency), as well as integer combinations of these frequencies. But since only relative frequencies are physically meaningful, there are really only two fundamental frequencies, which may be chosen to be g − g J and g S − g J , i.e. relative to g J . In other words, if z is regular, thenz ≡ ze −ig J t is a doubly periodic function with periods 2π/(g − g J ) and 2π/(g S − g J ). In Figure 2, we choose to plot the amplitude versus phase ofz whenever the second period completed an integer number of cycles. As long as z is regular, the value ofz at those times is a singly periodic function, and hence appears in the plot as a connected curve. By contrast, when z is chaotic, it appears as scattered points.

Equations of Motion
In this section, we consider the evolution of a test particle that comes under the influence of a single planet ("Jupiter") that has fixed values of eccentricity e J , inclination i J , apsidal precession rate g J ≡˙ J , and nodal precession rate s J ≡Ω J . This is a model for the case when a particle comes under the influence of two planetary Laplace-Lagrange modes, one eccentric and one inclined. (See the discussion of Hamiltonian [27].) As we show below, in the real Solar system the main modes affecting Mercury are the Jupiter-dominated eccentricity mode and the Venusdominated inclination mode, with the Venus-dominated eccentricity mode also playing a role. Therefore for application to Mercury, e J and g J refer to the amplitude and frequency of the Jupiter-dominated eccentricity mode, while i J and s J refer to those of the Venus-dominated inclination mode. Nonetheless, for the purposes of the present section it is simplest to assume that Jupiter is the only planet, and that both of its precession rates are enforced by divine intervention.
We evolve the secular equations for the test particle to leading nonlinear order. Because the test particle now has two coupled degrees of freedom, its evolution is more complicated than before, and there are many more terms to include in its Hamiltonian. At first ( § §4.1-4.4), we include only the following terms from Table 1: where J = 5 4 αe J . Note that the effect of Jupiter's eccentricity is diluted by a factor ∼ α 1, whereas the effect of its inclination is undiluted by any such factor. In §4.5 we add in all the remaining terms from Table  1, and show that these additional terms have little effect in the parameter regime of interest. The terms in Hamiltonian (28) that are second order in eccentricity or inclination (i.e., the first two terms, and the bracketed terms on the second line) are responsible for linear evolution. We keep only three nonlinear terms, ∝ |z| 4 , |ζ| 4 and |z| 2 |ζ| 2 . As we show in this subsection, these are responsible for nonlinear frequency shifts. And as we show in subsequent subsections, nonlinear frequency shifts are crucial for resonance overlap and chaos. Even though the frequency shifts might be small (second order in eccentricity and inclination), they can still be sufficient to shift the frequency into and out of secular resonance. One of the terms we drop is the Kozai resonance, i.e., the term c 26 × z * 2 ζ 2 + c.c. in Table 1. That term has little effect on Mercury's evolution, because its phase is rapidly varying, and hence the term nearly averages to zero for parameters similar to Mercury's ( §4.5). By contrast, the frequency-changing terms can never average to zero. 8 The equations of motion are Expressed in terms of the real canonical variables (eq. [12]), Hamiltonian (28) has two cosine terms (i.e., primary resonances), which become important when their arguments vary slowly. We first focus on the term ∝ cos( − g J t). Defining the free nonlinear apsidal frequency one would expect this resonance to be important near where as may also be inferred by an inspection of equation (29).
To nonlinear order, g is a function of eccentricity and inclination, because of the nonlinear frequency-changing terms included in Hamiltonian (28). Hence by varying e and i, one can alter the free precession frequency and bring it into secular resonance. In the p e -p i plane, the resonance traces out a one dimensional curve-or, in fact, a straight line to leading nonlinear order. The second cosine term in the above Hamiltonian, cos(Ω − s J t), behaves similarly. Defining one would expect it to be important near where We shall make these considerations more precise in §4.3, where we also work out the resonant widths, and show that in addition to the primary resonances are a multitude of secondary resonances. The y-axis is the time-averaged |ζ| 2 ≈ i 2 , and the x-axis is the timeaveraged |z| 2 ≈ e 2 . The initial conditions are on a uniform grid in |ζ| 2 and |z| 2 . Resonant bands are clearly visible, as is the zone of chaos where the bands overlap. Because of the averaging, chaotic orbits that explore the full extent of the overlap zone give rise to points at the center of the zone. Chaotic orbits that only partially explore the overlap zone (as is true of the real Mercury today; see below) give rise to chaotic points at the edges of the overlap zone. The orbits are very regular at small pe, p i ; for comparison, the real Mercury currently has pe ∼ .05 and p i ∼ .02. Its motion would be regular under this weak forcing.

Simulations: Maps of the Mean Momenta
We run suites of simulations of Hamiltonian (28), i.e., of equations (29)-(30). There are only five parameters, γ, g J , s J , J , and i J . The linear frequency γ sets the overall timescale, and can be scaled out. We choose g J /γ = 0.72 and s J /γ = −1.26, since these are close to the true values in the Solar system for the Jupiterdominated eccentricity mode (relative to Mercury's free precession frequency), and for the Venus-dominated inclination mode (see §5). We also set the excitation amplitudes equal to each other, J = i J , and present sequences of simulations with various amplitudes. The true value in the Solar system for the corresponding modes is, very roughly, J ∼ i J ∼ 0.01 (see §5 for more precise values).
The dynamics is that of two nonlinearly coupled harmonic oscillators, each of which is also nonlinear and is forced periodically. We have attempted many different methods for visualizing the integration results, such as using catalogs of surfaces of section or Fourier transforms. However, most methods were too complicated, and obscured the underlying simplicity of the dynamics, i.e. that it is the overlapping of resonances that drive chaos. In the end, we invented a new method, the Map of the Mean Momenta (MMM). This method has many advantages over the usual surfaces of section. It is somewhat similar to the frequency map analysis of Laskar (1990) (see below). Figure 3 maps the results from around 50, 000 numerical integrations at very small excitation, J = i J = .0003, using the MMM. Each point in the plot is the timeaveraged value of |z| 2 ≈ e 2 and |ζ| 2 ≈ i 2 from a single simulation, averaged over a time span of 3 × 10 4 /γ. Before taking the time average, we filter with a Hanning filter (Laskar 1993), which leads to a faster convergence of the averages (when they do converge). The simulations were initialized with values of |z| 2 and |ζ| 2 that were equally spaced on a grid, with the spacing in |z| 2 twice that in |ζ| 2 ; the initial phases were = π/2 and Ω = −π/2. Three kinds of motion are readily apparent in the MMM: (i) regular and non-resonant, (ii) regular and resonant, and (iii) chaotic. Most of the figure is covered with a regular grid of points that nearly traces the initial conditions. Here, the values of z and ζ remain regular and non-resonant throughout the simulation. A few resonant bands also appear in the figure, where the motion is also regular. Note that if the initial conditions span a resonant island, and if the motion remains regular, then the time-averaged momenta (|z| 2 and |ζ| 2 ) exhibit a sharp discontinuity. Inside the island, they average to their values near the island center, whereas outside the island's separatrix they average to a value offset from the center by a finite amount, of order the resonance width (or, more accurately, around 1/4 of the resonance fullwidth). This can be seen clearly in Figure 3, where the regular points at the center of the resonant bands represent librating particles. Also apparent in the figure are the chaotic trajectories. These show up as the cluster of irregular points near where resonant bands intersect. We have checked the Lyapunov exponent, as well as surfaces of section (see below), to verify that points that appear on the figure to be chaotic are truly chaotic.
Figures 4-5 show the MMM for simulations with the same parameters as in Figure 3, but with J and i J increased first to 0.003, and then to 0.01. With increasing forcing, the locations of the resonant bands do not change, but they get wider, and higher order resonances become visible. As a result, the zone of chaos where the bands overlap expands. Surprisingly, even with the seemingly modest forcing of J = i J = 0.01-values that are comparable to those in the real Solar System (see below)-the zone of chaos approaches very low values of e and i, and close to the values for the real Mercury.
Our method for displaying results, the MMM, is similar in philosophy to frequency map analysis (FMA; Laskar 1990). But whereas in FMA one plots the frequencies of the co-ordinates, here we plot the averages of the momenta. We have also performed the FMA (not shown); and when we convert from frequencies to momenta via the inverse of equations (31) and (33), the resulting maps are almost identical to the MMM. For the purposes of the present paper, we prefer the MMM because its axes are approximately e 2 and i 2 , which are simpler to interpret than the precession frequencies.

Theory: Resonance Locations and Widths, and
Zone of Chaos To develop understanding of the behavior seen in the MMM's, we first re-write Hamiltonian (28) in terms of  the real canonical variables (eqs. [8]-[9]), and then transform from (p e , ) to (p e , − ≡ − g J t) and from (p i , Ω) to (p i , Ω − ≡ Ω − s J t), which transforms the Hamiltonian to a time-independent one: where ∆ is the linear apsidal frequency mismatch (eq. [19]), and is the mismatch for the nodal frequencies. In the absence of the coupling term (∝ p e p i ) both degrees of freedom would evolve independently according to the equations of the second fundamental model ( §3.1).
In the following, we determine the location and width of the two primary resonances: the eccentricity resonance ([1,0]) and the inclination resonance ([0,1]). To do so, we ignore inclination forcing (setting i J = 0) when studying the eccentricity resonance, and vice versa. The system is trivially integrable if either i J = 0 or J = 0. In the former case, p i is constant because the Hamiltonian does not depend on Ω − . Therefore H − is equivalent to the coplanar Hamiltonian (eq. [24]), but with ∆ → ∆ − 2p i . This implies that the center of the eccentricity resonance is located at (eq. [25]) and the island has half-width at fixed p i given by δp e = 4( J √ p e * ) 1/2 → [1, 0] (eq. [26]). The interior of this island is plotted in Figure  6 as a blue band, with parameter values as chosen for the simulations of Figure 3. Comparing the two figures shows that the above analytic expressions agree well with the result of the numerical integrations. Note that we plot the half-width in Figure 6 rather than the full-width because the average momentum of an orbit that lies just outside of a separatrix is approximately half-way to the edge of the resonance. Figure 6 also shows the inclination resonance as a red band. Reasoning as before, its center and half-width (at fixed p e ) are In addition to these primary resonances are an infinite number of secondary ones. Far from resonances, we may ignore the cosine forcing terms since they tend to average to zero. The test particle's free precession frequencies relative to Jupiter are then Resonances are important near where which passes close to the origin for our choices of ∆ and ∆ s . In general, the slope of a resonance line in the p e -p i plane is (4n + m)/(n − 4m), and all [m, n] resonant lines intersect at the point in the p e -p i plane where g − = s − = 0, i.e. at (p e * * , p i * * ), where Although there are an infinite number of secondary resonances, most are very weak, i.e. their widths are small. The most prominent resonances in Figure 3 are the primary resonances [1,0] and [0,1], whose widths have been worked out above. Next most prominent are the [1, -1] and [1, 1] resonances, whose widths may be understood qualitatively as follows (see Appendix B for a quantitative calculation). The linear solution for z is a sum of two terms, the free and forced complex eccentricities (eq. [18]). Similarly, to linear order ζ is a sum of free and forced complex inclinations. Therefore to leading nonlinear order, the coupling term in the Hamiltonian (p e p i = |z| 2 |ζ| 2 ) can be written as a sum of terms, one of which has the form z where the subscript f denotes forced and φ denotes free. This term has frequency corresponding to the [1,-1] resonance; hence the width of this resonance is ∼ |e φ i φ e f i f |.
after defining the free values as and the forced values as where the asterisk denotes values at resonant center, and we neglect here the small difference between the lowercase momenta and the upper-case ones used in the appendix. As in linear theory, the forced eccentricity scales as the inverse of the frequency detuning (eq. [18]), although now it is the nonlinear frequency detuning g − (eq. [41]) that is relevant; similarly, i f is inversely proportional to s − . The above widths for the [1,±1] resonances are shown in Figure 6. The [1, -1] resonance leaves p e + p i nearly constant, which produces trajectories in the p e -p i plane that have a slope of -1 (see Appendix B; we again neglect the small difference between lower-and upper-case momenta). Therefore for each value of (p e * , p i * ) at resonant center (i.e., where g − − s − = 0), the upper and lower envelopes of the resonant band are at (p e * ± δp e , p i * ∓ δp e ). Comparing with Figure 3 shows excellent agreement between theory and simulations.
The central grey region in Figure 6 is the overlap zone of the above four resonances. Its shape is peculiar because each resonance induces a trajectory with a particular orientation in the p e -p i plane, and we take the overlapping region to be wherever one resonance can induce motion into a second resonance. The zone of chaos seen in Figure 3 is about twice as large as that predicted in Figure 6. This difference is not surprising because we plot the half-widths in Figure 6, whereas one might expect that chaos would begin where the full-widths overlap. Furthermore, we shall show that higher order combinations of the four strongest resonances also play a role in the chaos. Nonetheless, the grey zone provides a reasonable estimate of the zone of chaos. Figure 7 plots the widths calculated above when the forcing J and i J take on the values of Figure 4. The latter figure is overlayed on Figure 7 , showing excellent agreement between theory and simulation in the zones of regular motion. And, as before, the zone of chaos from the simulations is around twice as large as the region where the resonant half-widths overlap.
The dimensions of the zone of chaos may be estimated analytically. These estimates will be used later to infer the threshold for chaos in the Solar system. In Figures  6-7, the horizontal and vertical spikes of the chaotic zone are due to the [1, 0] and [0, 1] resonances; the extent of the spikes is simply estimated by the resonance width for the horizontal and vertical spikes respectively, where p e * * and p i * * are the co-ordinates of the zone center (eq.
[45]). The extent of the zone caused by the overlap be- It might appear surprising that the width of the chaotic zone due to the overlap of a primary resonance with a secondary one (eq. [51]) is not much smaller than that due to the overlap between two primary resonances, despite the fact that the width of a primary resonance is first order in eccentricity or inclination (eqs. [38]&[40]), whereas the width of the [1,-1] is second order (eq. [46]). The reason for this is that the [1,-1] resonance is enhanced by the denominator that appears in the forced 9 Focusing on the region to the lower left of the zone center in Figure 6 or 7, the center of the [1,-1] resonance is displaced from zone center in the pe-p i plane by the vector −(1, 3/5)x (eq. [44]); we take x > 0. The orbit given by the half-width of the [1,-1] is therefore displaced from zone center by the vector −(1, 3/5)x + (−1, 1)δpe, where δpe 2(pe * * p i * * ) 1/4 √ J i J /(17x/10) (eq. [46]). Equating that vector to the displacement of the center of the [1,0] resonance, i.e. to (−1, 1/4)x (eq. [37]) yields equation (51). One could also proceed to calculate the widths of higher order resonances. Extending our reasoning from above, the width of an [m, n] resonance in the p e -p i plane should scale as where is comparable to the typical eccentricity, and i is comparable to the typical inclination. However, as we have seen for the [1, ±1] resonances, near-resonant denominators can make the widths significantly larger than this naive estimate.

Surfaces of Section and Libration of Resonant Angles
Here, instead of the global map (MMM), we investigate a few particular trajectories in detail to demonstrate the chaotic behavior, using both surfaces of section and resonant angles. Figure 8 shows a number of surfaces of section from Hamiltonian (28). Since the transformed Hamiltonian (eq. [35]) is time independent with two degrees of freedom, we follow the usual practice of taking a section whenever the phase of one of the degrees of freedom (here, − ) executes an integer number of cycles. At those times, we plot the amplitude versus phase of the second degree of freedom (p i vs. Ω − ). We may also understand this form for the surface of section as follows. Hamiltonian (28) has four fundamental frequencies, g J , s J , g, and s where g and s are the nonlinear free frequencies of z and ζ, respectively. Only relative frequencies are physically meaningful, and there are three of these, which we may choose to be g − g J , s − s J , and g − s. But Hamiltonian (28) does not depend on − Ω, and hence the frequency g − s does not enter. 10 Therefore there are two remaining fundamental frequencies, as for the coplanar case ( §3.2). And, as described there, to examine the characteristics of the motion, one may take a section whenever the phase corresponding to one of the fundamental frequencies (here, the phase − g J t which corresponds to g − g J ) executes an integer number of cycles.
The three upper left panels of Figure 8 show surfaces of section with values of g J and s J as before, and with various values of J and i J . All the surfaces of section shown have the same (constant) value of energy, H − = 0.0154γ. To map out phase space would require many different energy values, but for that purpose the MMM is more useful. The lower right panel in Figure 8 shows a "double section" of the (blue) chaotic trajectory in the upper right panel, i.e. it shows the two momenta wherever both − and Ω − have executed a half-integer number of cycles. At these times, the cosine terms in the transformed Hamiltonian vanish, and all trajectories with a fixed energy H − fall along the same branch of a hyperbola. Of course double sections from all the trajectories shown in Figure 8 would lie along the same branch of the hyperbola because they all have the same value of H − .
The surface of section in the upper left panel has a relatively low J and i J , and the motion is mostly regular for the value of energy chosen. This can also be seen in the MMM (Fig. 4) near the relevant hyperbola branch. The [1, -1] resonance is clearly evident in the top left panel of Figure 8. Its half-width is δp i = 0.0055, as compared to the prediction of 0.006 from equation (46). The upper right panel of Figure 8 shows the case with higher forcing. With this higher forcing, the [1, -1] resonance is wider, and the region near its separatrix has broken up into a wide zone of chaos. Some of the higher order resonances are visible in this section. In the lower left panel, the forcing has been raised further. Even though J and i J are still relatively small compared to unity, the zone of chaos is vast.
The blue chaotic trajectory in the upper-right panel of Figure 8 behaves qualitatively like the real Mercury, and the parameters are also similar ( §5). Therefore we investigate it in more detail. From its surface of section, we see that the separatrix of the [1, -1] resonance 10 The full fourth order Hamiltonian does depend on − Ω because of the Kozai term; see Table 1. Therefore one cannot take surfaces of section of the full fourth order Hamiltonian, but one can still plot its MMM (Fig. 10).

Full Fourth Order Hamiltonian
Thus far we have focused on the truncated fourth order Hamiltonian (eq. [28]). Figure 10 shows the MMM of the full fourth order Hamiltonian, expanded to leading order in α (i.e., including all terms in Table 1 in Appendix A). From the similarity of Figure 10 to the truncated integrations of Figure 5, we conclude that the terms dropped in the truncated Hamiltonian have little effect on the dynamics, particularly in the region of small e and i (lower left corner of the MMM).
The dropped terms have little effect because the only new fundamental relative frequency they introduce is g − s (see first paragraph of §4.4). This "Kozai frequency"  (Table 1 in Appendix A), rather than the truncated Hamiltonian (28). From the fact that the two figures are broadly similar, one can infer that the terms dropped from Hamiltonian (28) are of small importance in the regime of interest.
differs significantly from zero in the domain of Figure 10, and hence it can only combine with the other two relative frequencies (g − g J and s − s J ) to give resonances at high order. To be quantitative, the resonant line of the Kozai frequency is at g − s ≈ γ(2 + (3/2)p e − (5/2)p i ) ≈ 0, i.e., it traces the line p i ≈ (4 + 3p e )/5. Hence the Kozai resonance is at much larger p i than shown in Figure 10.
Aside from the Kozai term (c 26 in Table 1), all other terms dropped from the truncated Hamiltonian depend on Jupiter's eccentricity or inclination. These do not introduce new forcing frequencies because Jupiter's frequencies already appear in the test particle's orbit at linear order-in its forced e and i. While the dropped terms do change the amplitudes of the forcing terms, the change is small as long as Jupiter's e and i is smaller than the e and i it linearly forces in the test particle, as is true of Figure 10.
We suspect that the terms dropped from Hamiltonian (28) are quite often of secondary importance. This is largely true for the real Mercury ( §5). And we suspect that it is true more generally because if secular interactions between two planets are strong, then the forced e's and i's will typically (though not always) be larger than the forcing ones. Nonetheless, the dropped terms can be important in certain circumstances; for example, the Kozai term will play a role if a planet has a high inclination, and MMR's will be important for planets whose orbital periods are near integer ratios.

Fourier Transforms
In §5 we shall make the connection to the real Mercury. For that purpose, it will prove instructive to examine trajectories in Fourier space.
For a more exact comparison to Mercury, we consider here the Hamiltonian which differs from Hamiltonian (28) by the inclusion of a constantγ to allow the linear apsidal and nodal precession rates to differ from each other (see footnote 8).
Note that we also change notation so that i V and s V are the amplitude and precession rate of the Venus mode.
We focus on a one dimensional family of systems parameterized by κ, which scales all eccentricities and inclinations. More precisely, in this "κ-model" we choose the parameters J = i V = 0.008κ; and initial conditions |z| = 0.16κ, |ζ| = 0.07κ, = Ω = π/2. The remaining parameters areγ = 0.9, g J = 0.72γ, s V = −1.14γ. With these parameters, the center of the [1,0] resonance is at p e * = 2(0.28 − 2p i ), as before; and the center of the [0,1] resonance is at p i * = −2(0.24 − 2p e ), whereas before the constant was 0.26 rather than 0.24. This difference is of little consequence. For displaying the results of the integration, we shall choose γ = 5.87 /yr. Our rationale for choosing these particular numerical values will be explained in §5. Figure 11 shows the Fourier transforms of z and ζ for the κ-model trajectory that has κ = 0.75. We normalize the Fourier transform of z aŝ and similarly for ζ, where T is the duration of the Fourier transform, which we choose in the present subsection to be T = 4400/γ. With this normalization, if z has constant amplitude and frequency, i.e. if z(t) = k 0 e iω0t , then |ẑ| = |k 0 | at frequencies close to ω 0 . The trajectory used for Figure 11 is quasiperiodic-the peaks in the Fourier transform are simply spikes, whose widths become narrower for larger T . We call the two largest peaks in the top panel the forced and free z. The forced z is at frequency g J = 0.72γ = 4.23 /yr. The free z is at frequency g M = 5.78 /yr. Because of nonlinearities, the free frequency (g M ) differs from the linear free frequency (γ) by a small but non-negligible amount. Similarly, in the bottom panel the largest two peaks are the forced ζ at frequency s V = −1.14γ, and the free ζ at frequency s M = −5.45 /yr, which differs from the linear free frequency −γγ.
In addition to the free and forced z and ζ, there are a multitude of peaks in Figure 11 that are generated by nonlinear couplings. The peaks inẑ all fall at frequencies g M + m(g M − g J ) + n(s M − s V ) for integers m, n. Roughly speaking, the peak amplitudes become smaller for larger values of |m| and |n|. These amplitudes can be calculated perturbatively, as is sketched in the following. As before, we define the free and forced components as (z φ , z f , ζ φ , ζ f ), which have phases that rotate with frequencies (g M , g J , s M , s V ), respectively. To leading nonlinear order, the nonlinear terms in the equation for dz/dt are proportional to z|z| 2 = (z φ + z f )|z φ + z f | 2 and z|ζ| 2 = (z φ + z f )|ζ φ + ζ f | 2 (eq. [29]). These generate six new frequencies in z: Each frequency-generating term acts as a linear forcing on z. Together with the free and forced z, these account for eight of the peaks marked in the top panel of Figure 11; specifically, they account for the two highest peaks in each of the four left-most triplets. The other peaks are accounted for by higher order nonlinear terms. One of these other peaksthe one labelled [1,-1]-is quite large, even though one might naively have expected that it would be smaller because it enters at a higher nonlinear order. The reason for this is that its forcing frequency differs from g M by σ ≡ (g M − g J ) − (s M − s V ) which is quite small. Hence this near resonance amplifies the peak by g M /σ ∼ 20. We note parenthetically that the width of the [1,-1] resonance (as described in §4.3) is directly related to the amplitudes of the three peaks at g M and g M ± σ. The Fourier transform of ζ behaves similarly to that of z, with the frequency peaks at s M + m(g M − g J ) + n(s M − s V ). Figure 12 shows the Fourier transforms for the κ-model at higher values of κ. The left-most panels show the case κ = 1. The motion is largely quasiperiodic, but the nonlinearly generated peaks have increased significantly relative to the κ = 0.75 case. At κ = 1.3 the motion is chaotic, and at κ = 1.55 it is highly chaotic.

MERCURY
We integrate the eight Solar system planets with the SWIFT symplectic integrator (Levison & Duncan 1994), supplemented with a routine for Mercury's relativistic precession (see Wu & Lithwick 2010, for code details). We initialize the planets with their current orbits and use their actual masses. The integration timestep is 8 days. Fig. 12.-Fourier transforms of results from three κ-model simulations, with κ = 1, 1.3, 1.55. At κ = 1 (left-most panels), the motion is still largely quasiperiodic. The amplitudes of the nonlinearly generated peaks have risen significantly relative to Fig. 11 (κ = 0.75), even though the forced and free z and ζ have changed by a modest amount. At κ = 1.3, there is weak chaos-the peaks have widened, and neighboring peaks overlap. At κ = 1.55, the trajectory is highly chaotic.
One might suspect that Mercury's orbital evolution is more complicated than our toy model for a variety of reasons: its e and i are not too small, and hence the fourth order expansion is approximate; it is not massless, and hence backreacts onto the other planets (especially Venus); there are seven other planets that do not have constant orbital elements and frequencies but participate in the overall chaos of the Solar system; and Mercury can be affected by resonant terms. Despite these complications, we show that the chaotic behavior of Mercury is qualitatively similar to the Hamiltonian model. This is perhaps not too surprising, since nonlinear dynamics are largely driven by resonances and their overlap. Hence as long as a model roughly captures the locations and widths of the principal resonances, it should produce qualitatively correct behavior.

Fourier Transforms
To compare with the κ-model, we first consider cleaner cases by pre-multiplying the current eccentricities and inclinations of all planets by the reduction factor κ nbody . Figure 13 shows Mercury's Fourier transforms in a κ nbody = 0.75 integration lasting T = 150Myr. Comparing this with the κ-model at κ = 0.75 (Fig. 11) shows broad agreement. In truth, the parameters for the κmodel were chosen to match the free and forced z and ζ seen in Figure 13, i.e. the frequencies and heights of the four peaks marked g J , g M , s V , s M . Since there were eight quantities to match, we could do this by adjusting eight parameters in the κ-model: γ,γ, g J , s V , J , i V , as well as the initial values of |z| and |ζ|. Therefore it is not significant that the forced and free peaks in the two figures agree. What is significant is that the other peaks that are generated by nonlinear couplings of the forced and free peaks also largely agree. This indicates that the κ-model captures much of the nonlinearity as seen in the real Mercury.
There are, however, at least three differences of note. First, the peaks in Figure 13 are broader than those in Figure 11. This is because Figure 13 suffers from weak chaos. But it is remarkable how sharp the largest peaks are: even though the e's and i's of the Solar system have only been reduced by 25%, the resulting chaos is surprisingly weak. Note that the integration intervals in the two figures are the same, T = 4400/γ = 150Myr, and hence the finite width of the peaks in Figure 13 is not due to the finite T . A second difference between the two figures is that theẑ M peak at frequency g M + (g M − g J ) is significantly larger in the κ nbody integration. That peak is so large because it is overlapped by a peak forced by Venus's eccentricity mode, which has precession frequency g V ≈ 2g M −g J . In other words, for κ nbody = 0.75, Mercury is in a secular resonance with a librating angle that corresponds to frequency 2g M − g V − g J , and this largely hides the effect of g V in the Fourier transform of Figure 13. The third difference between the two figures is the peak inζ M that is caused by Uranus's inclination mode. But this peak appears to have little dynamical consequence for Mercury.
In Figure 14, the factor multiplying the initial e's and  Figure 13, but with κ nbody increased from 0.75 to 0.95. The motion is more chaotic here, and Venus's eccentricity forcing at frequency g V is distinct. The blue dashed arrows denote frequency spacings of g V − g M , and the blue dashed vertical lines denote peaks due to the g V mode and nonlinear couplings generated by that mode.
i's has been raised to κ nbody = 0.95. The resulting motion is more chaotic, as the widths of the peaks are wider than before, especially forζ M . In addition, the frequencies have been shifted sufficiently to break Mercury from the 2g M − g V − g J resonance, and the effect of Venus's eccentricity forcing is distinct. Even though the motion is more chaotic, the principal forcing peaks and their harmonics are still identifiable. Figure 15 shows the Fourier transform of the real Mercury (κ nbody = 1). Even though the initial e's and i's have been increased by only 5% relative to Figure 14, the motion is significantly more chaotic, and the nonlinearly generated peaks are less easily identifiable, especially those near Mercury's free frequencies g M and s M . Nonetheless, we conclude from the progression of Figures  13-15, that the κ-model captures much of the physics. In particular, two modes-the Jupiter eccentricity mode and the Venus inclination mode-are primarily responsible for driving Mercury's chaos. The most important element lacking from the κ-model appears to be the extra forcing by the Venus eccentricity mode. The fact that the κ-model becomes chaotic at a higher threshold than the real Mercury (κ ∼ 1.3 vs. κ nbody ∼ 1) is likely partly due to that extra forcing (i.e. that extra forcing makes the real Mercury more chaotic). An additional contributor to the discrepancy between the two critical κ's is that the simple κ-model does not accurately capture nonlinear frequency shifts, while the precise values of the frequencies are important for where the resonances overlap. Despite this, the difference between the critical κ's is not large, and this lends support to our claimed origin for Mercury's chaos.
We note parenthetically that while we only focus on a narrow range of frequencies in Figure 15, Mercury also has peaks at |ω| ∼ 20 /yr, due to forcing by Earth and Mars. However, these peaks have amplitudes 10 −3 , and appear to have little influence on Mercury's chaotic motion. Had they been important, one would have expected to see their influence in Figure 13, whereas all the main peaks in that figure have already been identified. Figure 16 shows results from a 2 Gyr SWIFT integration of the full Solar system (with no reduction of the initial e's and i's). The black curve in the top panel is Mercury's total |z M |, which is very nearly equal to its total eccentricity (eq. [8]), and illustrates the chaotic behavior of Mercury's orbit. The overplotted green curve is the absolute value of Mercury's free z M , which we define to be the part of its total z M that comes from the main peak in Figure 15, i.e., we first take the Fourier transform of z M , then set to zero all frequencies except those satisfying 4.9 /yr < ω < 6.5 /yr, and then take the inverse Fourier transform. By plotting the free z M , the shortterm variations are reduced, and long-term diffusion is clearer. The second panel in Figure 16 is the same as the top but for ζ M ; for the free ζ M , we filter out frequencies outside of the domain −6.3 /yr < ω < −4.7 /yr.

Resonant Angles
The bottom two panels of Figure 16 show the two fourangle combinations involving Mercury that were found to undergo libration episodes. The third panel shows the angle ( M − J ) − (Ω M − Ω V ), which is the angle that has frequency i.e., the [1,-1] angle (eq. [43]). (More precisely, we use the angles of the free elements; see below.) And the bottom panel shows the angle associated with the frequency Laskar (1992) has shown that the σ angle can change from libration to circulation. But our finding that the σ can as well is new. 11 We note that even though both the σ and σ angles undergo libration episodes, the four-angle combination that is the sum of the two, i.e. the angle associated with 2g M − g V − g J does not. That angle was found to librate in the κ nbody = 0.75 simulation (Fig.  13), and we suspect that it will eventually librate in a long enough integration of the full Solar system. The angles displayed in the bottom panels of Figure  16 were those of the free elements. For example, for M we first filtered the Fourier transform of z M as described above, and took the phase of the free part of z M . This filtering procedure is especially important for Ω V , because Venus's ζ V variations are dominated by forcings due to other modes, including the Mercury-, Earth-, and Mars-dominated modes (Laskar 1990), whereas we wish Ω V to denote the phase of the Venus-dominated mode. Therefore, we first filter the Fourier transform of ζ V , keeping only frequencies −7.7 /yr < ω < −6.3 /yr, and use for Ω V the phase of the filtered ζ V . Similarly, we filter z M , ζ M , z V , and z J with appropriate windows to obtain the other angles of interest. Our method of filtering for extracting mode angles differs from that of Laskar (1990), who extracts mode angles by projecting onto the numerically determined nonlinear "proper modes." We  [m, n]. The angles alternately librate, showing that the chaos is at least partly caused by the overlap of these resonances, as in the κ-model (Fig. 9). The plotted angles are the phases of the free orbital elements.
have experimented with a number of different methods, and also with changing the size of the filter window, and the duration of the integration, and found that our filtering method is simple to implement, is computationally efficient, and gives reliable results.
In addition to the two four-angle combinations of Figure 16, one might suspect that there are many more higher-order combinations that librate when both of those angles simultaneously circulate, as in the Hamiltonian model. In Figure 17, we zoom into the episode when the [1, −1] angle first undergoes an extended period of circulation, and plot some higher order combinations associated with the frequencies [m, n] ≡ m(g M − g J ) + n(s M − s V ). It can be seen that these angles librate in turn, just as in the Hamiltonian model (compare with Fig. 9). This provides another demonstration that the physics of the Hamiltonian model is similar to that of the real Mercury.

SUMMARY AND DISCUSSION
We have shown how secular chaos is driven by the overlap of secular resonances, both for a test particle modelled with a simplified Hamiltonian, and for the real Mercury. To linear order, secular frequencies are constant. But nonlinearities can shift planets into and out of secular resonance with each other, and when two resonances overlap, chaos results.
In § §2-4, we focused on the evolution of a test particle in the presence of multiple massive planets. The test particle was evolved to leading nonlinear order, and the e's, i's, and precession rates of the planets (or more properly of the planet modes) were taken to be constant. We first considered the simple case with zero inclinations, as was first worked out by Sidlichovsky (1990). In §4, we generalized to non-zero inclinations, when the test particle comes under the influence of one eccentric and one inclined planet mode. In that case, the particle has two free frequencies, its apsidal and nodal frequencies (g and s). Each of these is altered by the particle's e and i. Therefore each resonance traces out a one-dimensional curve in the particle's e-i plane, or equivalently in its p e -p i plane. A simple way to map out the dynamics is with the "mean momentum map" (MMM), whereby the particle's timeaveraged p e and p i are plotted against each other for different initial conditions. This shows where the resonances are, how wide they are, and how their overlap leads to chaos . We calculated analytically the locations and widths of the four strongest resonancesthe [1,0], [0,1], and [1,±1]-and showed that these agreed with the numerical MMM results (Figs. 6-7). Chaos in this case emerges from the overlap of resonances of the form [m, n] (eq. [43]), with typically n = m ± 1 and m a small integer. This may be seen in the MMM, in surfaces of section (Fig. 8), and also by explicitly tracing the chain of librating angles (Fig. 9). We also examined the test particle's trajectories in Fourier space .
In §5, we considered the orbital evolution of Mercury in N-body simulations. We showed that despite all the simplifications we made in the Hamiltonian models, the real Mercury behaved in a qualitatively similar manner. In particular: • Mercury's chaos is primarily driven by the s V and g J modes (i.e. the Venus-and Jupiter-dominated i and e modes), although the g V mode also plays a role. The nonlinear couplings between those modes and Mercury's own free modes (with frequencies g M and s M ) are primarily responsible for Mercury's chaos (Figs. 13-15).
• There are a slew of resonant angles that drive Mercury's chaos. Just as in the Hamiltonian model, a chain of resonant angles of the form m( M − J ) + n(Ω M − Ω V ) show sequential librations, for integers [m, n] (Fig. 17), where the angles refer to the phases of the free orbital elements. In addition, we identified a new four-angle combination, ( M − V ) + (Ω M − Ω V ), that can also undergo libration episodes (Fig. 16).
• Mercury is perched on the threshold of chaos. If one reduces the e's and i's of the planets by only 25%, Mercury's motion becomes nearly regular (Fig. 13). This behavior is also apparent in the Hamiltonian model (Fig. 5). We have also performed a κ nbody = 1.2 simulation, in which the planets' initial e's and i's were increased 20% (not shown). The result was violent instability, with Mercury ejected in ∼ 100 Myr.
Having identified the secular resonances responsible for Mercury's chaos, and calculated the widths and locations of those resonances, we can calculate the threshold for Mercury's chaos. We do that here in an approximate way. First, since Mercury's apsidal and nodal frequencies differ from g J and s V by ∼ 20%, the co-ordinates in the p e -p i plane where the two resonances overlap are around half that, or p e * * ∼ p i * * ∼ 0.1 (eq. [45]). Second, the width of the chaotic overlap zone between the [1,-1] and the [1,0] (or [0,1]) is ∼ 2( J i J ) 1/4 (p e * * p i * * ) 1/8 (eq. [51]; we include here an extra factor of 2 to account for the difference between the half-and full-width, as described in §4.3) Therefore if J ∼ i J p 3/2 e * * /4 ∼ 0.01, then the region of chaos will encroach upon the origin of the p e -p i plane. This explains why Mercury can be chaotic even though the eccentricities and inclinations in the Solar system are at the level of a few percent.
The work discussed in this paper can be extended in a number of directions. The theory can be extended to order-unity eccentricities and inclinations. Although that case will be more complicated, we suspect that the basic structure will remain, with resonant zones in the e-i plane whose overlap leads to chaos. One can also attempt to build a theory that includes long-term diffusion and massive planets, as well as incorporating MMR's.
A number of applications also come to mind, such as quantifying Mercury's chaotic diffusion and understanding how it came about that Mercury is perched on the threshold of chaos. The latter seems to be a clue for understanding how the Solar system arrived at its current marginally stable state. It would also be interesting to investigate the role of σ (eq. [56]). Is it an unnecessary coincidence for Mercury's chaos?
where the first two terms generate the identity transformation, and the functions k e and k i are first order in J and i J ; their form will be chosen to "kill" the cosine terms in the Hamiltonian to leading order. The von Zeipel generating function transforms variables as follows p e = P e + k e cos q e (B3) p i = P i + k i cos q i (B4) Q e = q e + ∂ Pe k e sin q e (B5) Inserting into the Hamiltonian yields H(P e , Q e ; P i , Q i ) = − P 2 e − P 2 i 4 + ∆P e + ∆ s P i − 2P e P i − k e k i cos(Q e − Q i ) − k e k i cos(Q e + Q i ) to second order, after setting To leading order, k e is twice the product of the free eccentricity ( √ P e ) with the forced eccentricity, where the forced eccentricity differs from the linear expression ( J /∆; see eq. [18]) by the terms P e /2 + 2P i in the denominator, which arise from the nonlinear shift of the frequency (eqs. [41]- [42]). Similarly k i is twice the product of the free and forced inclinations. Hence the strengths of the [1, -1] and [1,1] resonances are proportional to the products of the free and forced eccentricities and inclinations, as argued qualitatively in §4.
To determine the width of the [1, -1] resonance, we drop the last cosine term in the above Hamiltonian. Since P + ≡ P e + P i is an integral of the motion, we may re-write the Hamiltonian as H(P e , Q − ) = 2(P e − P * ) 2 − k e k i cos Q − , dropping a constant and defining Q − ≡ Q e − Q i and P * ≡ (5/8)P + − (1/4)(∆ g − ∆ s ). Therefore the half-width of the resonance is δP e = |k e k i | .
We take the amplitude of the cosine term to be fixed at its value at resonance center (e.g., Chirikov 1979). Since P e + P i is constant, the half-width in P i is the same, δP i = |k e k i |.