SECULAR RESONANCE SWEEPING OF THE MAIN ASTEROID BELT DURING PLANET MIGRATION

The dynamical structure of the Kuiper Belt suggests that the outer solar system experienced a phase of planetesimal-driven migration in its early history (Fernandez & Ip 1984; Malhotra 1993, 1995; Hahn & Malhotra 1999; Levison et al. 2008). Pluto and other Kuiper Belt objects that are trapped in mean motion resonances (MMRs) with Neptune are explained by the outward migration of Neptune due to interactions with a more massive primordial planetesimal disk in the outer regions of the solar system (Malhotra 1993, 1995). In addition, the so-called scattered disk of the Kuiper Belt can also be explained by the outward migration of Neptune (Hahn & Malhotra 2005) or by the effects of a high eccentricity phase of ice giant planet evolution during the outward migration of Neptune (Levison et al. 2008). The basic premise of planetesimal-driven migration is that the giant planets formed in a more compact configuration than we find them today, and that they were surrounded by a massive (∼50 M_⊕) disk of unaccreted icy planetesimals that was the progenitor of the currently observed Kuiper Belt (Hahn & Malhotra 1999). When planetesimals are preferentially scattered either inward (toward the Sun) or outward (away from the Sun), net orbital angular momentum is transferred between the disk and the large body, causing a drift in the large body's semimajor axis (Fernandez & Ip 1984; Kirsh et al. 2009). In many simulations of giant planet migration, icy planetesimals are preferentially scattered inward by each of the three outer giant planets (Saturn, Uranus, and Neptune) causing these planets to migrate outward. Due to Jupiter's large mass, planetesimals that encounter Jupiter are preferentially ejected out of the solar system, leading to a net loss of mass from the solar system and an inward migration of Jupiter.

Planetesimal-driven giant planet migration has been suggested as a cause of the Late Heavy Bombardment (LHB; Gomes et al. 2005; Strom et al. 2005); however, the link between these two events has yet to be definitively established (Chapman et al. 2007; Ćuk et al. 2010; Malhotra & Strom 2010). Such migration would have enhanced the impact flux of both asteroids and comets onto the terrestrial planets in two ways. First, many of the icy planetesimals scattered by the giant planets would have crossed the orbits of the terrestrial planets. Second, as the giant planets migrated, locations of mean motion and secular resonances would have swept across the asteroid belt, raising the eccentricities of asteroids to planet-crossing values.

Recently, Minton & Malhotra (2009) showed that the patterns of depletion observed in the asteroid belt are consistent with the effects of sweeping of resonances during the migration of the giant planets. The Jupiter-facing sides of some of the Kirkwood gaps (regions of the asteroid belt that are nearly empty due to strong Jovian MMRs) are depleted relative to the Sun-facing sides, as would be expected due to the inward migration of Jupiter and the associated inward sweeping of the Jovian MMRs. The region within the inner asteroid belt between semimajor axis range 2.1–2.5 AU also has excess depletion relative to a model asteroid belt that was uniformly populated and then subsequently sculpted by the gravitational perturbations of the planets over 4 Gyr, as would be expected due to the outward migration of Saturn and the associated inward sweeping of a strong secular resonance, the so-called ν₆ resonance, as explained below. In our 2009 study, we concluded that the semimajor axis distribution of asteroids in the main belt is consistent with the inward migration of Jupiter and outward migration of Saturn by amounts proposed in previous studies based on the Kuiper Belt resonance structure (e.g., Malhotra 1995). However, in that study the migration timescale was not strongly constrained, because only the relative depletion of asteroids in nearby semimajor axis bins could be determined, not their overall level of depletion. In the present paper, we explore in more detail the effect that planet migration would have had on the asteroid belt due to asteroid eccentricity excitation by the sweeping of the ν₆ secular resonance. From the observed eccentricity distribution of main belt asteroids, we find that it is possible to derive constraints on the secular resonance sweeping timescale, and hence on the migration timescale.

Secular resonances play an important role in the evolution of the main asteroid belt. The inner edge of the belt nearly coincides with the ν₆ secular resonance which is defined by g ≈ g₆, where g is the rate of precession of the longitude of pericenter, ϖ, of an asteroid and g₆ is the sixth eigenfrequency of the solar system planets (approximately the rate of precession of Saturn's longitude of pericenter). The ν₆ resonance is important for the delivery of near-Earth asteroids to the inner solar system (Scholl & Froeschle 1991). Williams & Faulkner (1981) showed that the location of the ν₆ resonance actually forms surfaces in a − e − sin i space, and Milani & Knezevic (1990) showed that those surfaces approximately define the "inner edge" of the main asteroid belt. However, as mentioned above, Minton & Malhotra (2009) found that, with the planets in their present configuration, planetary perturbations over the age of the solar system cannot fully account for the detailed orbital distribution of the asteroids in the inner asteroid belt. The pattern of excess depletion in inner asteroid belt noted by Minton & Malhotra (2009) is consistent with the effect of the inward sweeping of the ν₆ secular resonance. In general, the direction of motion of the ν₆ is anticorrelated with that of Saturn, so an inward sweeping of the ν₆ would be produced by an outwardly migrating Saturn.

Sweeping, or scanning, secular resonances have been analyzed in a number of previous works. Sweeping secular resonances due to the changing quadrupole moment of the Sun during solar spin-down have been explored as a possible mechanism for explaining the eccentricity and inclination of Mercury (Ward et al. 1976). Secular resonance sweeping due to the effects of the dissipating solar nebula just after planet formation has also been investigated as a possible mechanism for exciting the orbital eccentricities of Mars and of the asteroid belt (Heppenheimer 1980; Ward 1981). The dissipating massive gaseous solar nebula would have altered the secular frequencies of the solar system planets in a time-dependent way, causing locations of secular resonances to possibly sweep across the inner solar system, thereby exciting asteroids into the eccentric and inclined orbits that are observed today. This mechanism was revisited by Nagasawa et al. (2000), who incorporated a more sophisticated treatment of the nebular dispersal. However, O'Brien et al. (2007) have argued that the excitation (and clearing) of the primordial asteroid belt was unlikely due to secular resonance sweeping due to the dispersion of the solar nebula.

The special case of asteroids on initially circular orbits being swept by the ν₆ and ν₁₆ resonances has been investigated by Gomes (1997). In this paper, we consider the more general case of non-zero initial eccentricities; our analysis yields qualitatively new results and provides new insights into the dynamical history of the asteroid belt. This extends the work of Ward et al. (1976) and Gomes (1997) in developing analytical treatments of the effects of sweeping secular resonances on asteroid orbits. In doing so, we have developed an explicit relationship between the migration rate of the giant planets, the initial eccentricity of the asteroid and its initial longitude of perihelion, and the final eccentricity of the asteroid after the passage of the resonance. We show that for initially non-zero asteroid eccentricity, the sweeping of the ν₆ resonances can either increase or decrease asteroid eccentricities. Examining the orbits of observed main belt asteroids we find evidence for a double-peaked eccentricity distribution; this supports the case for a history of ν₆ sweeping. Quantitative comparison of our analytical theory with the semimajor axis and eccentricity distribution of asteroids yields new constraints on the timescale of planet migration.

We note that although our analysis is carried out in the specific context of the sweeping ν₆ resonance during the phase of planetesimal-driven migration of Jupiter and Saturn, the techniques developed here may be extended to other similar problems, for example, the sweeping of the inclination-node ν₁₆ resonance in the main asteroid belt, the ν₈ secular resonance in the Kuiper Belt, and farther afield, the sweeping of secular resonances in circumstellar or other astrophysical disks.

We adopt a simplified model in which a test particle (asteroid) is perturbed only by a single resonance, the ν₆ resonance. We use a system of units where the mass is in solar masses, the semimajor axis is in units of AU, and the unit of time is (2π)⁻¹ yr. With this choice, the gravitational constant, G, is unity. An asteroid's secular perturbations close to a secular resonance can be described by the following Hamiltonian function (Malhotra 1998):

$$\begin{equation} H_{{\rm sec}}=-g_0J+\varepsilon \sqrt{2J}\cos (w_p-\varpi), \end{equation} \tag{ 1 }$$

where w_p = g_pt + β_p describes the phase of the pth eigenmode of the linearized eccentricity–pericenter secular theory for the solar system planets (Murray & Dermott 1999), g_p is the associated eigenfrequency, ϖ is the asteroid's longitude of perihelion, $J=\sqrt{a}(1-\sqrt{1-e^2})$ is the canonical generalized momentum which is related to the asteroid's orbital semimajor axis a and eccentricity e, and −ϖ and J are the canonically conjugate pair of variables in this 1-degree-of-freedom Hamiltonian system. The coefficients g₀ and ε are given by

$$\begin{equation} g_0=\frac{1}{4a^{3/2}}\sum\nolimits _j\alpha _j^2b^{(1)}_{3/2}(\alpha _j)m_j, \end{equation} \tag{ 2 }$$

$$\begin{equation} \varepsilon =\frac{1}{4a^{5/4}}\sum\nolimits _j\alpha _j^2b^{(2)}_{3/2}(\alpha _j)m_jE_j^{(p)}, \end{equation} \tag{ 3 }$$

where the subscript j refers to a planet, E^(p)_j is the amplitude of the g_p mode in the jth planet's orbit, α_j = min {a/a_j, a_j/a}, m_j is the ratio of the mass of planet j to the Sun, and b⁽¹⁾_3/2(α_j) and b⁽²⁾_3/2(α_j) are Laplace coefficients; the sum is over all major planets. The summations in Equations (2) and (3) are over the eight major planets, for the greatest accuracy; however, we will adopt the simpler two-planet model of the Sun–Jupiter–Saturn in Section 3, in which case we sum over only the indices referring to Jupiter and Saturn; then g_p is an eigenfrequency of the secular equations of the two-planet system.

With fixed values of the planetary masses and semimajor axes, g₀, g_p, and ε are constant parameters in the Hamiltonian given in Equation (1). However, during the epoch of giant planet migration, the planets' semimajor axes change secularly with time, so that g₀, g_p, and ε become time-dependent parameters. In the analysis below, we neglect the time dependence of g₀ and ε, and adopt a simple prescription for the time dependence of g_p (see Equation (8) below). This approximation is physically motivated: the fractional variation of g₀ and ε for an individual asteroid is small compared to the effects of the "small divisor" g₀ − g_p during the ν₆ resonance sweeping event.

It is useful to make a canonical transformation to new variables (ϕ, P) defined by the following generating function:

$$\begin{equation} {\cal F}(-\varpi,P,t) = (w_p(t)-\varpi)P. \end{equation} \tag{ 4 }$$

Thus, $\phi = \partial {\cal F}/\partial P = (w_p(t)-\varpi)$ and $J=-\partial {\cal F}/\partial \varpi =P$. The new Hamiltonian function is $\tilde{H}_{{\rm sec}}=H_{{\rm sec}}+\partial {\cal F}/\partial t$,

$$\begin{equation} \tilde{H}_{{\rm sec}}=(\dot{w}_p(t)-g_0)J+\varepsilon \sqrt{2J}\cos \phi, \end{equation} \tag{ 5 }$$

where we have retained J to denote the canonical momentum, since P = J. It is useful to make a second canonical transformation to canonical eccentric variables,

$$\begin{equation} x = \sqrt{2J}\cos \phi, \qquad y=-\sqrt{2J}\sin \phi, \end{equation} \tag{ 6 }$$

where x is the canonical coordinate and y is the canonically conjugate momentum. The Hamiltonian expressed in these variables is

$$\begin{equation} \tilde{H}_{{\rm sec}}=(\dot{w}_p(t)-g_0)\frac{x^2+y^2}{2}+\varepsilon x. \end{equation} \tag{ 7 }$$

As discussed above, during planetary migration, the secular frequency g_p is a slowly varying function of time. We approximate its rate of change, $\dot{g}_p=2\lambda$, as a constant, so that

$$\begin{equation} \dot{w}_p(t) =g_{p,0}+2\lambda t. \end{equation} \tag{ 8 }$$

We define t = 0 as the epoch of exact resonance crossing, so that g_{p, 0} = g₀ (Ward et al. 1976). Then, $\dot{w}_p(t)-g_0=2\lambda t$, and the equations of motion from the Hamiltonian of Equation (7) can be written as

$$\begin{equation} \dot{x}=2\lambda t y, \end{equation} \tag{ 9 }$$

$$\begin{equation} \dot{y}={-}2\lambda t x-\varepsilon. \end{equation} \tag{ 10 }$$

These equations of motion form a system of linear, nonhomogenous differential equations, whose solution is a linear combination of a homogeneous and a particular solution. The homogeneous solution can be found by inspection, giving

$$\begin{equation} x_h(t)=c_1\cos \lambda t^2+c_2\sin \lambda t^2, \end{equation} \tag{ 11 }$$

$$\begin{equation} y_h(t)={-}c_1\sin \lambda t^2+c_2\cos \lambda t^2, \end{equation} \tag{ 12 }$$

where c₁ and c₂ are constant coefficients. We use the method of variation of parameters to find the particular solution. Accordingly, we replace the constants c₁ and c₂ in the homogeneous solution with functions A(t) and B(t), to seek the particular solution of the form

$$\begin{equation} x_p(t)=A(t)\cos \lambda t^2+B(t)\sin \lambda t^2, \end{equation} \tag{ 13 }$$

$$\begin{equation} y_p(t)={-}A(t)\sin \lambda t^2+B(t)\cos \lambda t^2. \end{equation} \tag{ 14 }$$

Substituting this into the equations of motion we now have

$$\begin{equation} \dot{A}\cos \lambda t^2+\dot{B}\sin \lambda t^2=0, \end{equation} \tag{ 15 }$$

$$\begin{equation} {-}\dot{A}\sin \lambda t^2+\dot{B}\cos \lambda t^2=-\varepsilon; \end{equation} \tag{ 16 }$$

therefore,

$$\begin{equation} \dot{A}=\varepsilon \sin \lambda t^2, \end{equation} \tag{ 17 }$$

$$\begin{equation} \dot{B}=-\varepsilon \cos \lambda t^2. \end{equation} \tag{ 18 }$$

Equations (17) and (18) do not have a simple closed-form solution, but their solution can be expressed in terms of Fresnel integrals (Zwillinger 1996). The Fresnel integrals are defined as follows:

$$\begin{equation} S(t)=\int _0^t\sin t^{\prime 2} dt^{\prime }, \end{equation} \tag{ 19 }$$

$$\begin{equation} C(t)=\int _0^t\cos t^{\prime 2} dt^{\prime }, \end{equation} \tag{ 20 }$$

and have the following properties:

$$\begin{equation} S(-t)=-S(t), \end{equation} \tag{ 21 }$$

$$\begin{equation} C(-t)=-C(t), \end{equation} \tag{ 22 }$$

$$\begin{equation} S(\infty)=C(\infty)=\sqrt{\frac{\pi }{8}}. \end{equation} \tag{ 23 }$$

Therefore,

$$\begin{equation} A(t)= \frac{\varepsilon }{\sqrt{|\lambda |}}S(t\sqrt{|\lambda |}), \end{equation} \tag{ 24 }$$

$$\begin{equation} B(t)=-\frac{\varepsilon }{\sqrt{|\lambda |}}C(t\sqrt{|\lambda |}). \end{equation} \tag{ 25 }$$

We denote initial conditions with a subscript i, and write the solution to Equations (9) and (10) as

$$\begin{eqnarray} x(t)&=&x_i\cos \big[\lambda \big(t^2-t_i^2\big)\big]+y_i\sin \big[\lambda \big(t^2-t_i^2\big)\big]\nonumber \\ && +\frac{\varepsilon }{\sqrt{|\lambda |}}[(S-S_i)\cos \lambda t^2-(C-C_i)\sin \lambda t^2],\qquad \end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray} y(t)&=&{-}x_i\sin \big[\lambda \big(t^2-t_i^2\big)\big]+y_i\cos \big[\lambda \big(t^2-t_i^2\big)\big]\nonumber \\ && -\frac{\varepsilon }{\sqrt{|\lambda |}}[(C-C_i)\cos \lambda t^2+(S-S_i)\sin \lambda t^2].\qquad \end{eqnarray} \tag{ 27 }$$

Because the asteroid is swept over by the secular resonance at time t = 0, we can calculate the changes in x, y by letting t_i = −t_f and evaluating the coefficients C_i, C_f, S_i, S_f far from resonance passage, i.e., for $t_f\sqrt{|\lambda |}\gg 1$, by use of Equation (23). Thus, we find

$$\begin{equation} x_f=x_i+\varepsilon \sqrt{\frac{\pi }{2|\lambda |}}\big[\cos \lambda t_i^2-\sin \lambda t_i^2\big], \end{equation} \tag{ 28 }$$

$$\begin{equation} y_f=y_i- \varepsilon \sqrt{\frac{\pi }{2|\lambda |}}\big[\cos \lambda t_i^2+\sin \lambda t_i^2\big]. \end{equation} \tag{ 29 }$$

The final value of J long after resonance passage is therefore given by

$$\begin{eqnarray} J_f&=&\frac{1}{2}\big(x_f^2+y_f^2\big)\nonumber \\ &=&\frac{1}{2}\big(x_i^2+y_i^2\big)+\frac{\pi \varepsilon ^2}{2|\lambda |} + \varepsilon \sqrt{\frac{\pi }{2|\lambda |}}\big[x_i\big(\cos \lambda t_i^2-\sin \lambda t_i^2\big)\nonumber \\ && -y_i\big(\cos \lambda t_i^2+\sin \lambda t_i^2\big)\big]\nonumber \\ &=&J_i+\frac{\pi \varepsilon ^2}{2|\lambda |}+\varepsilon \sqrt{\frac{2\pi J_i}{|\lambda |}}\cos \left(\phi _i-\lambda t_i^2-\frac{\pi }{4}\right). \end{eqnarray} \tag{ 30 }$$

With a judicious choice of the initial time, t_i, and without loss of generality, the cosine in the last term becomes cos ϖ_i, and therefore

$$\begin{equation} J_f=J_i+\frac{\pi \varepsilon ^2}{2|\lambda |}+\varepsilon \sqrt{\frac{2\pi J_i}{|\lambda |}}\cos \varpi _i. \end{equation} \tag{ 31 }$$

The asteroid's semimajor axis a is unchanged by the secular perturbations; thus, the changes in J reflect changes in the asteroid's eccentricity e. For asteroids with non-zero initial eccentricity, the phase dependence in Equation (31) means that secular resonance sweeping can potentially both excite and damp orbital eccentricities. We also note that the magnitude of eccentricity change is inversely related to the speed of planet migration. In linear secular theory (Equation (1)), eccentricity and inclination are decoupled, and therefore the effect of the sweeping ν₆ does not depend on the inclination. However, as Williams & Faulkner (1981) showed, the location of the ν₆ does depend on inclination, but the dependence is weak for typical inclinations of main belt objects. Nevertheless, there are populations of main belt asteroids at high inclination (such as the Hungaria and Phocaea families), and an analysis of secular resonance sweeping that incorporates coupling between eccentricity and inclination would be valuable for understanding the effects of planet migration on these populations; we leave this to a future investigation.

For small e, we can use the approximation $J\simeq {\frac{1}{2}}\sqrt{a}e^2$. Considering all possible values of cos ϖ_i ∈ { − 1, +1}, an asteroid with initial eccentricity e_i that is swept by the ν₆ resonance will have a final eccentricity in the range e_min to e_max, where

$$\begin{equation} e_{\rm min,max}\simeq \left|e_i\pm \delta _e\right| \end{equation} \tag{ 32 }$$

and

$$\begin{equation} \delta _e\equiv \left|\varepsilon \sqrt{\frac{\pi }{|\lambda |\sqrt{a}}}\right|. \end{equation} \tag{ 33 }$$

Equations (31)–(33) have the following implications.

• 1.  
  Initially circular orbits become eccentric, with a final eccentricity δ_e.
• 2.  
  An ensemble of orbits with the same a and initial non-zero e but uniform random orientations of pericenter are transformed into an ensemble that has eccentricities in the range e_min to e_max; this range is not uniformly distributed because of the cos ϖ_i dependence in Equation (31), rather the distribution peaks at the extreme values (see Figure 1 below).
• 3.  
  An ensemble of asteroids having an initial distribution of eccentricities which is a single-peaked Gaussian (and random orientations of pericenter) would be transformed into one with a double-peaked eccentricity distribution.

Zoom In Zoom Out Reset image size

In light of the above calculations, it is possible to conclude that if the asteroid belt were initially dynamically cold, that is, asteroids were on nearly circular orbits prior to secular resonance sweeping, then the asteroids would be nearly uniformly excited to a narrow range of final eccentricities, the value of which would be determined by the rate of resonance sweeping. Because asteroids having eccentricities above planet-crossing values would be unlikely to survive to the present day, it follows that an initially cold asteroid belt which is uniformly excited by the ν₆ sweeping will either lose all its asteroids or none. On the other hand, an initially excited asteroid belt, that is, a belt with asteroids that had non-zero eccentricities prior to the ν₆ sweeping, would have the asteroids' final eccentricities bounded by Equation (32), allowing for partial depletion and also broadening of its eccentricity distribution. In this section, we apply our theoretical analysis to the problem of the ν₆ resonance sweeping through the asteroid belt and compare the theoretical predictions with the observed eccentricity distribution of asteroids.

3.1. Parameters

In order to apply the theory, we must find the location of the ν₆ resonance as a function of the semimajor axes of the giant planets orbits, and also obtain values for the parameter ε (Equation (3)), for asteroids with semimajor axis values in the main asteroid belt. The location of the ν₆ resonance is defined as the semimajor axis, $a_{\nu _6}$, where the rate, g₀ (Equation (2)), of pericenter precession of a massless particle (or asteroid) is equal to the g₆ eigenfrequency of the solar system. In the current solar system, the g₆ frequency is associated with the secular mode with the most power in Saturn's eccentricity–pericenter variations. During the epoch of planetesimal-driven planet migration, Jupiter migrated by only a small amount but Saturn likely migrated significantly more (Fernandez & Ip 1984; Malhotra 1995; Tsiganis et al. 2005), so we expect that the variation in location of the ν₆ secular resonance is most sensitive to Saturn's semimajor axis. We therefore adopt a simple model of planet migration in which Jupiter is fixed at 5.2 AU and only Saturn migrates. We neglect the effects of the ice giants Uranus and Neptune, as well as secular effects due to the previously more massive Kuiper Belt and asteroid belt. With these simplifications, the g₆ frequency varies with time as Saturn migrates, so we parameterize g₆ as a function of Saturn's semimajor axis. In contrast with the variation of g₆, there is little variation of the asteroid's apsidal precession rate, g₀, as Saturn migrates. Thus, finding $a_{\nu _6}$ is reduced to calculating the dependence of g₆ on Saturn's semimajor axis.

To calculate g₆ as a function of Saturn's semimajor axis, we proceed as follows. For fixed planetary semimajor axes, the Laplace–Lagrange secular theory provides the secular frequencies and orbital element variations of the planets. This is a linear perturbation theory, in which the disturbing function is truncated to secular terms of second order in eccentricity and first order in mass (Murray & Dermott 1999). In the planar two-planet case, the secular perturbations of planet j, where j = 5 is Jupiter and j = 6 is Saturn, are described by the following disturbing function:

$$\begin{equation} R_{j}=\frac{n_j}{a_j^2}\left[\frac{1}{2}A_{jj}e_j^2+A_{jk}e_5e_6\cos (\varpi _5-\varpi _6)\right], \end{equation} \tag{ 34 }$$

where n is the mean motion and $\mathbf {A}$ is a matrix with elements

$$\begin{equation} A_{jj}=+n_j\frac{1}{4}\frac{m_k}{M_\odot +m_j}\alpha _{jk}\bar{\alpha }_{jk}b_{3/2}^{(1)}(\alpha _{jk}), \end{equation} \tag{ 35 }$$

$$\begin{equation} A_{jk}=-n_j\frac{1}{4}\frac{m_k}{M_\odot +m_j}\alpha _{jk}\bar{\alpha }_{jk}b_{3/2}^{(2)}(\alpha _{jk}), \end{equation} \tag{ 36 }$$

for j = 5, 6, k = 6, 5, and j ≠ k; α_jk = min {a_j/a_k, a_k/a_j}, and

$$\begin{equation} \bar{\alpha }_{jk}=\left\lbrace \begin{array}{@{}ll}1 & {:} a_j>a_k \\ a_j/a_k & {:} a_j<a_k. \\ \end{array} \right. \end{equation} \tag{ 37 }$$

The secular motion of the planets is then described by a set of linear differential equations for the eccentricity vectors, e_j(sin ϖ_j, cos ϖ_j) ≡ (h_j, k_j),

$$\begin{equation} \dot{h}_j=+\sum _{p=5}^6A_{pj}k_j,\;\;\;\;\dot{k}_j=-\sum _{p=5}^6A_{pj}h_j. \end{equation} \tag{ 38 }$$

For fixed planetary semimajor axes, the coefficients are constants, and the solution is given by a linear superposition of eigenmodes:

$$\begin{equation} \lbrace h_j,k_j\rbrace =\sum _pE_j^{(p)}\lbrace \cos (g_pt+\beta _p),\sin (g_pt+\beta _p)\rbrace, \end{equation} \tag{ 39 }$$

where g_p are the eigenfrequencies of the matrix $\mathbf {A}$ and E^(p)_j are the corresponding eigenvectors; the amplitudes of the eigenvectors and the phases β_p are determined by initial conditions. In our two-planet model, the secular frequencies g₅ and g₆ depend on the masses of Jupiter, Saturn, and the Sun and on the semimajor axes of Jupiter and Saturn.

For the current semimajor axes and eccentricities of Jupiter and Saturn the Laplace–Lagrange theory gives frequency values g₅ = 37 yr⁻¹ and g₆ = 223 yr⁻¹, which are lower than the more accurate values given by Brouwer & van Woerkom (1950) by 14% and 20%, respectively (Laskar 1988). Brouwer & van Woerkom (1950) achieved their more accurate solution by incorporating higher-order terms in the disturbing function involving 2λ₅–5λ₆, which arise due to Jupiter and Saturn's proximity to the 5:2 resonance (the so-called Great Inequality). By doing an accurate numerical analysis (described below), we found that the effect of the 5:2 resonance is only important over a very narrow range in Saturn's semimajor axis. More significant is the perturbation owing to the 2:1 near-resonance of Jupiter and Saturn. Malhotra et al. (1989) developed corrections to the Laplace–Lagrange theory to account for the perturbations from n + 1: n resonances in the context of the Uranian satellite system. Applying that approach to our problem, we find that the 2:1 near-resonance between Jupiter and Saturn leads to zeroth-order corrections to the elements of the $\mathbf {A}$ matrix.⁵ Including these corrections, we determined the secular frequencies for a range of values of Saturn's semimajor axis; the result for g₆ is shown in Figure 2(a) (dashed line).

Zoom In Zoom Out Reset image size

We also calculated values for the eccentricity–pericenter eigenfrequencies by direct numerical integration of the full equations of motion for the two-planet, planar solar system. In these integrations, Jupiter's initial semimajor axis was 5.2 AU, Saturn's semimajor axis, a₆, was one of 233 values in the range 7.3–10.45 AU, initial eccentricities of Jupiter and Saturn were 0.05, and initial inclinations were zero. The initial longitude of pericenter and mean anomalies of Jupiter were ϖ_{5, i} = 15° and λ_{5, i} = 92°, and Saturn were ϖ_{6, i} = 338° and λ_{6, i} = 625. In each case, the planets' orbits were integrated for 100 myr, and a Fourier transform of the time series of the {h_j, k_j} yields their spectrum of secular frequencies. For regular (non-chaotic) orbits, the spectral frequencies are well defined and are readily identified with the frequencies of the secular solution. The g₆ frequency as a function of Saturn's semimajor axis was obtained by this numerical analysis; the result is shown by the solid line in Figure 2 (a).

The comparison between the numerical analysis and the Laplace–Lagrange secular theory indicates that the linear secular theory, including the corrections due to the 2:1 near-resonance, is an adequate model for the variation in g₆ as a function of a₆. We adopted the latter for its convenience in the needed computations. The value of $a_{\nu _6}$ as a function of Saturn's semimajor axis was thus found by solving for the value of asteroid semimajor axis where g₀ = g₆; g₀ was calculated using Equation (2) and g₆ is the eigenfrequency associated with the p = 6 eigenmode (at each value of Saturn's semimajor axis). The result is shown in Figure 2(b).

We also used the analytical secular theory to calculate the eigenvector components E⁽⁶⁾_j in the secular solution of the two-planet system, for each value of Saturn's semimajor axis. We adopted the same values for the initial conditions of Jupiter and Saturn as in the direct numerical integrations discussed above. Finally, we computed the values of the parameter ε at each location $a_{\nu _6}$ of the secular resonance. The result is plotted in Figure 3. Despite the complexity of the computation, the result shown in Figure 3 is approximated well by a simple exponential curve, $\varepsilon \approx 3.5\times 10^{-9}\exp (2a_{\nu _6}/\mathrm{AU})$, in the semimajor axis range $2<a_{\nu _6}/\mathrm{AU}<4$.

Zoom In Zoom Out Reset image size

3.2. Four Test Cases

We checked the results of our analytical model against four full numerical simulations of the restricted four-body problem (the Sun–Jupiter–Saturn system with test particle asteroids) in which the test particles in the asteroid belt are subjected to the effects of a migrating Saturn. The numerical integration was performed with an implementation of a symplectic mapping (Wisdom & Holman 1991; Saha & Tremaine 1992), and the integration step size was 0.01 yr. Jupiter and Saturn were the only massive planets that were integrated, and their mutual gravitational influence was included. The asteroids were approximated as massless test particles. The current solar system values of the eccentricity of Jupiter and Saturn were adopted and all inclinations (planets and test particles) were set to zero. An external acceleration was applied to Saturn to cause it to migrate outward linearly and smoothly starting at 8.5 AU at the desired rate. As shown in Figure 2(b), the ν₆ resonance was located at 3 AU at the beginning of the simulation and swept inward past the current location of the inner asteroid belt in all simulations. In each of the four simulations, 30 test particles were placed at 2.3 AU and given different initial longitudes of pericenter spaced 12° apart. The semimajor axis value of 2.3 AU was chosen because it is far away from the complications arising due to strong MMRs. The only parameters varied between each of the four simulations were the initial osculating eccentricities of the test particles, e_i, and the migration speed of Saturn, $\dot{a}_{6}$. The parameters explored were

• 1.  
  e_i = 0.2, $\dot{a}_{6}=1.0\, \mathrm{AU}\,{\rm Myr}^{-1}$;
• 2.  
  e_i = 0.2, $\dot{a}_{6}=0.5\, \mathrm{AU}\, {\rm Myr}^{-1}$;
• 3.  
  e_i = 0.1, $\dot{a}_{6}=1.0\, \mathrm{AU}\,{\rm Myr}^{-1}$;
• 4.  
  e_i = 0.3, $\dot{a}_{6}=1.0\, \mathrm{AU}\, {\rm Myr}^{-1}$.

These values were chosen to illustrate the most relevant qualitative features. The migration rates of 0.5 AU Myr⁻¹ and 1.0 AU Myr⁻¹ are slow enough so that the change in eccentricity is substantial, but not so slow that the nonlinear effects at high eccentricity swamp the results. These test cases illustrate both how well the analytical model matches the numerical results, and where it breaks down.

Two aspects of the analytical model were checked. First, the perturbative equations of motion, Equations (9) and (10), were numerically integrated, and their numerical solution compared with that from the direct numerical integration of the full equations of motion. For the perturbative solution, we adopted values for λ that were approximately equivalent to the values of $\dot{a}_{6}$ in the full numerical integrations. Second, the eccentricity bounds predicted by the analytical theory, Equation (32), were compared with both numerical solutions. The results of these comparisons for the four test cases are shown in Figure 4. We find that the analytically predicted values of the maximum and minimum final eccentricities (shown as horizontal dashed lines) are in excellent agreement with the final values of the eccentricities found in the numerical solution of the perturbative equations, and in fairly good agreement with those found in the full numerical solution. Not surprisingly, we find that the test particles in the full numerical integrations exhibit somewhat more complicated behavior than the perturbative approximation, and Equation (32) somewhat underpredicts the maximum final eccentricity: this may be due to higher-order terms in the disturbing function that have been neglected in the perturbative analysis and which become more important at high eccentricity; effects due to close encounters with Jupiter also become important at the high eccentricities.

Zoom In Zoom Out Reset image size

3.3. Comparison with Observed Asteroid Eccentricities

Does the eccentricity distribution of main belt asteroids retain features corresponding to the effects of the ν₆ resonance sweeping? To answer this question, we need to know the main belt eccentricity distribution free of observational bias, and also relatively free of the effects of ∼4 Gyr of collisional evolution subsequent to the effects of planetary migration. We therefore obtained the proper elements of the observationally complete sample of asteroids with absolute magnitude H ⩽ 10.8 from the AstDys online data service (Knežević & Milani 2003); we excluded from this set the members of collisional families as identified by Nesvorný et al. (2006). These same criteria were adopted in Minton & Malhotra (2010) in a study of the long-term dynamical evolution of large asteroids. This sample of 931 main belt asteroids is a good approximation to a complete set of large asteroids that have been least perturbed by either dynamical evolution or collisional evolution since the epoch of the last major dynamical event that occurred in this region of the solar system; therefore, this sample likely preserves best the post-migration orbital distribution of the asteroid belt. The proper eccentricity distribution of these asteroids is shown in Figure 5. This distribution has usually been described in the literature by simply quoting its mean value (and sometimes a dispersion; Murray & Dermott 1999; O'Brien et al. 2007). Our best-fit single-Gaussian distribution to this data has a mean, μ_e and standard deviation, σ_e, given by μ_e = 0.135 ± 0.00013 and σ_e = 0.0716 ± 0.00022, and is plotted in Figure 5. However, we also note (by eye) a possible indication of a double-peak feature in the observed population. Our best-fit double-Gaussian distribution (modeled as two symmetrical Gaussians with the same standard deviation, but different mean values) to the same data has the following parameters:

$$\begin{eqnarray*} \mu ^{\prime }_{e,1}&=0.0846\pm 0.00011,\\ \mu ^{\prime }_{e,2}&=0.185\pm 0.00012,\\ \sigma _e^{\prime }&=0.0411\pm 0.00020. \end{eqnarray*}$$

More details of how these fits were obtained are described in Appendix A, where we also discuss goodness of fit of the single- and double-Gaussian distributions. A Kolmogorov–Smirnov (K-S) test shows that there is only a 4.5% probability that the observed eccentricities actually are consistent with the best-fit single-Gaussian distribution and a 73% probability that they are consistent with the best-fit double-Gaussian distribution; while the double Gaussian is apparently a better fit to the data compared to the single Gaussian, a K-S probability of only 73% is quite far from a statistically significant level of confidence. A dip test for multi-modality in the observed eccentricity distribution (Hartigan & Hartigan 1985) is also inconclusive; the likely reason for this is that our data sample is not very large, and the separation between the two putative peaks is too small in relation to the dispersion. Nevertheless, we bravely proceed in the next section with the implications of a double-Gaussian eccentricity distribution, with the caveat that some of the conclusions we reach are based on this statistically marginal result.

Zoom In Zoom Out Reset image size

By relating the g₆ secular frequency to the semimajor axis of Saturn, $\dot{g}_6$ can be related to the migration rate of Saturn, $\dot{a}_{6}$. In this section only the effects of the sweeping ν₆ secular resonance will be considered, and effects due to sweeping Jovian MMRs will be ignored. Because Jupiter is thought to have migrated inward a much smaller distance than Saturn migrated outward during planetesimal-driven migration, the effects due to migrating Jovian MMRs were likely confined to narrow regions near strong resonances (Minton & Malhotra 2009). In the inner asteroid belt between 2.1 and 2.8 AU, these would include the 3:1 and 5:2 resonances, currently located at approximately 2.5 AU and 2.7 AU, respectively. As shown in Figure 2(b), plausible parameters for the outward migration of Saturn would have allowed the ν₆ to sweep across the entire inner asteroid belt. Therefore, the ν₆ resonance would have been the major excitation mechanism across the 2.1–2.8 AU region of the main belt (and possibly across the entire main asteroid belt, depending on Saturn's pre-migration semimajor axis) during giant planet migration.

We used the results of our analytical model to set limits on the rate of migration of Saturn assuming a linear migration profile, with the caveat that many important effects are ignored, such as asteroid–Jupiter MMRs and Jupiter–Saturn MMRs (with the exception of the 2:1 resonance). We have confined our analysis to only the region of the main belt spanning 2.1–2.8 AU. Beyond 2.8 AU strong Jovian MMRs become more numerous. Due to the high probability that the icy planetesimals-driving planet migration would be ejected from the solar system by Jupiter, Jupiter likely migrated inward. The migration of Jupiter would have caused strong Jovian MMRs to sweep the asteroid belt, causing additional depletion beyond that of the sweeping ν₆ resonance (Minton & Malhotra 2009). A further complication is that sweeping Jovian MMRs may have also trapped icy planetesimals that entered the asteroid belt region from their source region beyond Neptune (Levison et al. 2009). The effects of these complications are reduced when we consider only the inner asteroid belt. From Figure 2(b), we find that the ν₆ would have swept the inner asteroid belt region between 2.1 and 2.8 AU when Saturn was between ∼8.5 and 9.2 AU. Therefore, the limits on $\dot{a}_6$ that we set using the inner asteroid belt as a constraint are only applicable for this particular portion of Saturn's migration history.

Our theoretically estimated final eccentricity as a function of initial asteroid semimajor axis and eccentricity is shown in Figure 6 for three different adopted migration rates of Saturn. The larger the initial asteroid eccentricities, the wider the bounds in their final eccentricities. If we adopt the reasonable criterion that an asteroid is lost from the main belt when it achieves a planet-crossing orbit (that is, crossing the orbits of either Jupiter or Mars) and that initial asteroid eccentricities were therefore confined to ≲ 0.4, then from Figure 6 Saturn's migration rate must have been $\dot{a}_6\gtrsim 0.15\, \mathrm{AU}\,{\rm Myr}^{-1}$ when the ν₆ resonance was sweeping through the inner asteroid belt. Our results indicate that if Saturn's migration rate had been slower than 0.15 AU Myr⁻¹ when it was migrating across ∼8.5–9.2 AU, then the inner asteroid belt would have been completely swept clear of asteroids by the ν₆ resonance.

Zoom In Zoom Out Reset image size

In light of our analysis and the observed dispersion of eccentricities in the asteroid belt (Figure 5), we can also immediately conclude that the pre-migration asteroid belt between 2.1 and 2.8 AU had significantly non-zero eccentricities. This is because, as discussed at the start of Section 3, an initially cold asteroid belt swept by the ν₆ resonance would either lose all its asteroids or none, and very low initial eccentricities would result in final asteroid eccentricities in a very narrow range of values (see Equation (32)), in contradiction with the fairly wide eccentricity dispersion that is observed. This conclusion supports recent results from studies of planetesimal accretion and asteroid and planet formation that the asteroids were modestly excited at the end of their formation (e.g., Petit et al. 2002).

In the Appendix we show that the double-Gaussian distribution is a slightly better fit to the main belt asteroid eccentricity distribution, but the statistical tests do not rule out a single-peaked distribution. We boldly proceed with considering the implications of the double-peaked eccentricity distribution to further constrain the migration rate of Saturn, with the caveat that these results can only be said to be consistent with the observations, rather than uniquely constrained by them.

If the pre-sweeping asteroid belt had a Gaussian eccentricity distribution, then the lower peak of the post-sweeping asteroid belt should be equal to the lower bound of Equation (32). We use the analytical theory to make a rough estimate of the parameter λ (and hence $\dot{a}_6$) that would yield a final distribution with lower peak near 0.09 and upper peak near 0.19 (which is similar to the best-fit double-Gaussian in Figure 5).

Applying Equation (32), we see that there are two possible solutions: 〈e_i〉 = 0.14,δ_e = 0.05 and 〈e_i〉 = 0.05,δ_e = 0.14. A corresponding migration rate of Saturn can be estimated from the value of δ_e using Equation (33), and the parameter relationships plotted in Figures 2 and 3. The former solution (δ_e = 0.05) requires a migration rate for Saturn of $\dot{a}_6=30\, \mathrm{AU}\, {\rm Myr}^{-1}$. We mention this implausible solution here for completeness, but we will not discuss it any further. The latter solution (δ_e = 0.14) requires a migration rate for Saturn of $\dot{a}_6=4\, \mathrm{AU}\, {\rm Myr}^{-1}$. We dub this solution the "cold belt" solution. This rate is comparable to the rates of planet migration found in the "Jumping Jupiter" scenario proposed by Brasser et al. (2009). A third solution exists if we consider that eccentricities in the main belt are restricted by the orbits of Mars and Jupiter on either side, such that stable asteroid orbits do not cross the planetary orbits. This limits asteroid eccentricities to values such that neither the aphelion of the asteroid crosses the perihelion distance of Jupiter, nor the perihelion of the asteroid crosses the aphelion distance of Mars. Maximum asteroid eccentricity is therefore a function of semimajor axis, where e_max = min (1 − Q_Mars/a − 1, q_Jupiter/a − 1), where Q and q are planet aphelion and perihelion, respectively, and a is the semimajor axis of the asteroid. In this case, an initial single-Gaussian eccentricity distribution with a mean greater than ∼0.3 would be severely truncated; therefore, we need only fit the lower peak of the double-Gaussian distribution at e = 0.09. Applying Equation (32), we find that δ_e = 0.21 provides a good fit. The corresponding migration rate of Saturn is $\dot{a}_6=0.8\, \mathrm{AU}\, {\rm Myr}^{-1}$. We dub this solution the "hot belt" solution.

We illustrate the two possible solutions for an ensemble of hypothetical asteroids having semimajor axes uniformly distributed randomly in the range 2.1–2.8 AU. In Figure 7(a) the initial eccentricity distribution is modeled as a Gaussian distribution with a mean 〈e_i〉 = 0.05 and a standard deviation of 0.01. This initial standard deviation was chosen so that the final standard deviation would be the same as that of the observed main belt. Figure 7(b) shows the eccentricity distribution after ν₆ resonance sweeping has occurred due to the migration of Saturn at a rate of 4 AU Myr⁻¹. The final distribution was calculated with Equation (31); we used values of ε shown in Figure 3, and the value of λ was calculated with the aid of Figure 2 which relates the value of g_p to the semimajor axis of Saturn. As expected, when an ensemble of asteroids with a single-peaked eccentricity distribution is subjected to the sweeping secular resonance, the result is a double-peaked eccentricity distribution. Because of the slight bias toward the upper limit of the eccentricity excitation band, proportionally more asteroids are found in the upper peak.

Zoom In Zoom Out Reset image size

In Figure 7(c) the initial eccentricity distribution is modeled as a truncated Gaussian: a Gaussian with mean 〈e_i〉 = 0.4 and standard deviation 0.1, but truncated at the semimajor axis-dependent Mars-crossing value. We used Equation (31) to calculate the eccentricity distribution of this hypothetical ensemble after ν₆ resonance sweeping with $\dot{a}_6=0.8\, \mathrm{AU}\, {\rm Myr}^{-1}$. Again, allowing that only those asteroids whose final eccentricities are below the Mars-crossing value will remain, the resulting post-migration eccentricity distribution is shown in Figure 7(d). In this case, we find the lower peak at the same eccentricity value as the lower peak in the observed main belt distribution (see Figure 5).

In both cases of possible solutions (initially cold main belt with 〈e_i〉 = 0.05 and $\dot{a}_6=4\, \mathrm{AU}\, {\rm Myr}^{-1}$; and initially hot main belt with 〈e_i〉 = 0.4 and $\dot{a}_6=0.8\, \mathrm{AU}\, {\rm Myr}^{-1}$), the theoretical models yield an excess of asteroids with eccentricities greater than 0.2 than in the observed main belt. However, as shown by Minton & Malhotra (2010), on gigayear timescales, the e ≳ 0.2 population of the asteroid belt is dynamically more unstable than the e ≲ 0.2 population. Thus, both solutions may be consistent with the observations, as post-sweeping dynamical erosion could result in a final eccentricity distribution resembling more closely the observed distribution.

The estimates of Saturn's migration rate quoted above depend strongly on the eccentricities of the giant planets during their migration. In deriving the above estimates, we adopted the present values of the giant planets' orbital eccentricities. The ν₆ resonance strength coefficient ε (Equation (3)) is proportional to the amplitude of the p = 6 mode, which is related to the eccentricities of the giant planets (namely, Saturn and Jupiter). From Equation (39), and the definition e_j(sin ϖ_j, cos ϖ_j) ≡ (h_j, k_j), the value of E^(p)_j is a linear combination of the eccentricities of the giant planets. Because Saturn is the planet with the largest amplitude of the p = 6 mode, from Equation (31) the relationship between the sweep rate and the value of Saturn's eccentricity is approximately λ_min∝e²₆. Therefore, to increase the limiting timescale by a factor of 10 would only require that the giant planets' eccentricities were ∼0.3 × their current value (i.e., e_{5, 6} ≈ 0.015).

Based on the existence of the inner asteroid belt, we conclude that Saturn's migration rate must have been ≳ 0.15 AU Myr^–1 as Saturn migrated from 8.5 to 9.2 AU (as the ν₆ resonance migrated from 2.8 to 2.1 AU). Migration rates lower than ∼0.15 AU Myr^–1 would be inconsistent with the survival of any asteroids in the inner main belt, as the ν₆ secular resonance would have excited asteroid eccentricities to planet-crossing values. This lower limit for the migration rate of Saturn assumes that Jupiter and Saturn had their current orbital eccentricities; the migration rate limit is inversely proportional to the square of the amplitude of the g₆ secular mode; if Jupiter and Saturn's eccentricities were ∼0.3 × their current value (i.e., e_{5, 6} ≈ 0.015) during the planet migration epoch, the limit on the migration rate decreases by a factor of ∼10. (This caveat also applies to the migration rate limits quoted below.)

Our analysis of secular resonance sweeping predicts that a single-peaked eccentricity distribution will be transformed into a double-peaked eccentricity by secular resonance sweeping. We find that the observed eccentricity distribution of asteroids may be consistent with a double-peaked distribution function, although we acknowledge that the statistics are poor. We used a double-Gaussian function to model the observed asteroid eccentricity distribution to set even tighter, albeit model-dependent constraints on the migration rate of Saturn.

We identified two possible migration rates that depend on the pre-migration dynamical state of the main asteroid belt. The first, the "cold primordial belt" solution, has an asteroid belt with an initial eccentricity distribution modeled as a Gaussian with 〈e_i〉 = 0.05; Saturn's migration rate of 4 AU Myr⁻¹ yields a final eccentricity distribution consistent with the observed asteroid belt. The second, the "hot primordial belt" solution, has an asteroid belt with an initial eccentricity distribution modeled as a Gaussian with 〈e_i〉 = 0.4, σ_e = 0.1, but truncated above the Mars-crossing value of eccentricity; in this case, Saturn's migration rate of 0.8 AU Myr⁻¹ is generally consistent with the observed asteroid belt.

Each of these solutions has very different implications for the primordial excitation and depletion of the main asteroid belt. The cold belt solution, with $\dot{a}_6=4\, \mathrm{AU}\, {\rm Myr}^{-1}$, would lead to little depletion of the asteroid belt during giant planet migration, as the ν₆ resonance would be unable to raise eccentricities to Mars-crossing values. This implies an initial dynamically quite cold asteroid belt with not more than ∼2 times the mass of the current main belt; the latter estimate comes from accounting for dynamical erosion over the age of the solar system (Minton & Malhotra 2010).

The hot belt solution, with $\dot{a}_6=0.8\, \mathrm{AU}\, {\rm Myr}^{-1}$, would lead to loss of asteroids directly due to excitation of eccentricities above planet-crossing values, by about a factor of ∼2. In this case, the main asteroid belt was more dynamically excited prior to resonance sweeping than we find it today. This implies much greater loss of asteroids prior to ν₆ sweeping, as the peak of the eccentricity distribution would be near the Mars-crossing value and subject to strong dynamical erosion (Minton & Malhotra 2010).

Each of these solutions has different implications for the model that the LHB of the inner solar system is linked to the epoch of planetesimal-driven giant planet migration (Gomes et al. 2005; Strom et al. 2005). These implications will be explored in a future paper.

We remind the reader that in order to elucidate the effects of ν₆ resonance sweeping, we have made a number of simplifying assumptions to arrive at an analytically tractable model. These simplifications include neglecting the effects of planets other than Jupiter and Saturn, the effects of sweeping Jovian MMRs on asteroids, the effects of a presumed massive Kuiper Belt during the epoch of planet migration, and the self-gravity and collisional interactions of a previously more massive asteroid belt. In addition, our analysis was carried out in the planar approximation, thereby neglecting any eccentricity–inclination coupling effects. These neglected effects can be expected to reduce somewhat the lower limit on Saturn's migration speed that we have derived, because in general they would reduce the effectiveness of the ν₆ in exciting asteroid eccentricities. Perhaps more importantly, giant planet migration would also lead to the sweeping of the main asteroid belt by the ν₁₆ inclination secular resonance (Williams & Faulkner 1981) whose effects could be used to infer additional constraints.

Recently, Morbidelli et al. (2010) found through numerical simulations with slow rates of planet migration (e-folding timescales exceeding 5 Myr) that the surviving asteroids in the main belt tend to be clumped around MMRs. The semimajor axis distribution of survivors is found very different from the observed distribution for asteroids, and also that a large proportion of survivors have inclinations above 20°, both inconsistent with the current main asteroid belt; comparisons made with the unbiased main belt asteroid distributions as described in Minton & Malhotra (2009). This is likely due to the effects of the sweeping ν₁₆ inclination longitude of ascending node secular resonance, analogous to the sweeping ν₆ eccentricity–pericenter secular resonance that we analyzed in the present study. While the effects of the sweeping ν₁₆ resonance are analogous to the ν₆, only affecting inclinations instead of eccentricities, a full analysis of the asteroid belt inclinations is beyond the scope of the present work, but will be explored in a future study.

A number of other studies have derived limits on the speed of planetesimal-driven giant planet migration. Murray-Clay & Chiang (2005) exclude an e-folding migration timescale τ ⩽ 1 Myr to 99.65% confidence based on the lack of a large observed asymmetry in the population of Kuiper Belt objects in the two libration centers of the 2:1 Neptune MMR. Boué et al. (2009) exclude τ ⩽ 7 Myr based on the observed obliquity of Saturn. The latter lower limit on the migration timescale is slightly incompatible with the lower limit on the rate of Saturn's migration of $\dot{a}_6>0.15\, \mathrm{AU}\, {\rm Myr}^{-1}$ we derive based on the existence of the inner asteroid belt. One way these can be reconciled is if Saturn's orbital eccentricity were a factor ∼2 smaller than its present value as it migrated from 8.5 to 9.2 AU, then some mechanism would need to have increased Saturn's eccentricity up to its present value by the time Saturn reached its present semimajor axis of ∼9.6 AU.

The authors thank the anonymous reviewer and the editor Eric Feigelson for useful comments. This research was supported in part by NSF grant no. AST-0806828 and NASA:NESSF grant no. NNX08AW25H. The work of David Minton was additionally partially supported by NASA NLSI/CLOE research grant no. NNA09DB32A

The binned eccentricity distribution may be modeled as a Gaussian probability distribution function, given by

$$\begin{equation} p(x)=\frac{1}{\sigma \sqrt{2\pi }}\exp \left[{-}\frac{(x-\mu)^2}{2\sigma ^2}\right], \end{equation} \tag{ A1 }$$

where σ is the standard deviation, μ is the mean, and x is the random variable; in our case x is the eccentricity. With an appropriate scaling factor, Equation (A1) can be used to model the number of asteroids per eccentricity bin. However, rather than fit the binned distribution, we instead perform a least-squares fit of the unbinned sample to the Gaussian cumulative distribution function given by

$$\begin{equation} P(x)=\frac{1}{2}+\frac{1}{2}\mathrm{erf}\left({-}\frac{x-\mu }{\sigma \sqrt{2}}\right). \end{equation} \tag{ A2 }$$

For the eccentricities of our sample of 931 main belt asteroids, the best-fit parameters are

$$\begin{eqnarray*} \mu _e&=0.135\pm 0.00013,\\ \sigma _e&=0.0716\pm 0.00022. \end{eqnarray*}$$

We also fit the data to a double-Gaussian distribution,

$$\begin{equation} p_2(x)=\frac{A^{\prime }}{\sigma ^{\prime }\sqrt{2\pi }}\left\lbrace \exp \left[-\frac{(x-\mu ^{\prime }_1)^2}{2\sigma ^{\prime 2}}\right]+\exp \left[{-}\frac{(x-\mu ^{\prime }_2)^2}{2\sigma ^{\prime 2}}\right]\right\rbrace. \end{equation} \tag{ A3 }$$

The cumulative distribution function for Equation (A3) is

$$\begin{equation} P_2(x)=\frac{1}{2}+\frac{1}{4}\left[\mathrm{erf}\left({-}\frac{x-\mu ^{\prime }_1}{\sigma \sqrt{2}}\right)+\mathrm{erf}\left({-}\frac{x-\mu ^{\prime }_2}{\sigma \sqrt{2}}\right)\right]. \end{equation} \tag{ A4 }$$

For the eccentricities of our sample of 931 main belt asteroids, we performed a least-squares fit to Equation (A4) and obtained the following best-fit parameters:

$$\begin{eqnarray*} \mu ^{\prime }_{e,1}&=0.0846\pm 0.00011,\\ \mu ^{\prime }_{e,2}&=0.185\pm 0.00012,\\ \sigma _e^{\prime }&=0.0411\pm 0.00020. \end{eqnarray*}$$

We evaluated the goodness of fit using the K-S test. This determines the probability that two distributions are the same, or in our case how well our model distributions fit the observed data (Press et al. 1992). The K-S test compares the cumulative distribution of the data against the model cumulative distribution function. We found that our asteroid sample has a probability of 4.5 × 10⁻² that it comes from the best-fit single Gaussian (Equation (A2)), but a probability of 0.73 that it comes from the double Gaussian (Equation (A4)). Therefore, the K-S tests indicate that the double Gaussian is a better fit to the data than the single Gaussian.

We performed Hartigan's dip test (Hartigan & Hartigan 1985) to test whether the observational data are consistent with a multi-peaked distribution. Hartigan's dip test calculates the probability that the distribution being tested has a single peak. Applying Hartigan's dip test to a given distribution yields in a test statistic; together with the sample size, the test statistic is matched to a p-value range in a precomputed table provided by Hartigan & Hartigan (1985). The p-value is a measure of the probability that the distribution actually has only one peak (the null hypothesis, for this problem). The smaller the calculated p-value, the less likely is the distribution to have a single peak and the more likely it is to have at least two significant peaks. A p-value of <0.05 indicates that the null hypothesis is very unlikely and that the given distribution has more than one peak. We applied the dip test to the eccentricity distribution of our sample of 931 main belt asteroids; we determined that the test statistic is 0.0107. From Hartigan & Hartigan (1985), this corresponds to a p-value range of 0.6 < p < 0.7 (based on a sample size of 1000). This indicates that the Hartigan dip test does not rule out the null hypothesis, i.e., a single-peaked distribution cannot be ruled out.

To further aid the interpretation of this test, we compare this result to that obtained by applying Hartigan's dip test to synthetic distributions that were explicitly double peaked by construction. Ten model distributions (of 1000 eccentricity values) were generated from the double-Gaussian function of Equation (A4), with parameter values μ₁ = 0.0846, μ₂ = 0.185, and σ = 0.0411 (same parameter values as the best fit for our sample of asteroids). Only the random seed was varied between each model distribution. Applying the Hartigan dip test, we find that the test statistic ranged between 0.00852 and 0.0126, corresponding to p-values between 0.98 and 0.3. This means that, according to the dip test, the null hypothesis—i.e., a single-peaked distribution—could not be ruled out for any of the model distributions (since none had p < 0.05), despite the fact that they were each drawn from an explicitly double-peaked distribution by construction. This illustrates that even with a sample size of nearly 1000, the dispersion in eccentricity around the two peaks is too large in comparison to the distance between the peaks and that the dip test is not sufficiently sensitive to detect the double-peaked underlying distribution. We interpret this to mean that the results of both the K-S test and Hartigan's dip test indicate that the main asteroid belt eccentricity distribution is consistent with being drawn from a double-peaked distribution, but that this cannot be definitely shown.