SPEED LIMIT ON NEPTUNE MIGRATION IMPOSED BY SATURN TILTING

It is now well accepted that the solar system was more compact after the protoplanetary gas disk dissipated. Then planets migrated due to interactions with the primordial planetesimal belt. The Nice model (Gomes et al. 2005; Tsiganis et al. 2005; Morbidelli et al. 2005) gives a unified scenario of this planetary migration, but it is still not fully constrained. For example, the Nice model allowed two possible classes of late evolution (Nesvorný et al. 2007). In the first one, called "class MA," Neptune is scattered to 22–25 AU and reaches its final orbit by slowly migrating over more than 5 AU. In the second class, labeled "DE," Neptune is placed to its current orbital distance with large eccentricity ≈0.3 and then slowly circularizes. Besides, the Nice model does not constrain inclinations, and the timescale of this late evolution is uncertain. Nevertheless, constraints on the migration timescale were obtained from the distribution of the Kuiper belt on the one hand (Murray-Clay & Chiang 2005), and from the distribution of the main asteroid belt on the other (Minton & Malhotra 2009). They both assumed an MA migration type without long-term evolution of eccentricities and inclinations. The former obtained a migration timescale τ between 1 and 10 Myr, and the latter found τ ≲ 0.5 Myr.

In this Letter, we aim at giving new constraints based on Saturn tilting (Ward & Hamilton 2004; Hamilton & Ward 2004). According to Ward and Hamilton, Saturn's large obliquity, = 2673919 (Helled et al. 2009), is due to a resonance capture between its spin axis and Neptune's orbit. Given the large uncertainties on Saturn's precession rate $-0\farcs 75\pm 0\farcs 21\;{\rm {yr}}^{-1}$ (Ward & Hamilton 2004), the more accurate regression of Neptune's orbit plane s₈ = −0692 yr⁻¹ (Laskar et al. 2004) is indeed included in the errorbars. Ward & Hamilton (2004) assume that today the two frequencies are equal. In their scenario, the norm of the frequency of Neptune's ascending node was initially larger, and then it captured Saturn's spin axis as it decreased due to Neptune's migration and/or the dissipation of the planetesimal disk (Ward & Canup 2006). In their numerical model, they took a quasiperiodic model of the solar system and forced an exponential evolution of the frequency s₈. Here, we show that Saturn can tilt in both migration classes, and that it gives a lower limit on the migration timescale. This limit depends on Neptune's initial inclination.

A recent paper by Helled et al. (2009) seems to contradict Ward and Hamilton's scenario. It gives a new estimate of Saturn's precession rate −07542 ± 00002 yr⁻¹ that is incompatible with a resonance with s₈. We show that with this value, Saturn can still evolve to its current state but that it is very unlikely. We discuss this result in our conclusion.

Here we recall the equations of motion of a planet axis and give the current dynamical state of Saturn's spin axis. The evolution of the spin axis w of a planet in a fixed reference frame (i, j, k), where k is the direction of the total orbital angular momentum, is given by

$$\begin{equation} \frac{d\mathbf {w}}{dt} = -\alpha ({\mathbf {n}}\cdot {\mathbf {w}})\mathbf {n}\times \mathbf {w} \end{equation} \tag{ 1 }$$

where n = ^t(n_x, n_y, n_z) is the normal to the orbit, and α is the precession constant. Without planetary perturbations, n is fixed and the spin axis w precesses uniformly around n with constant obliquity cos  = n · w. However, in a multiplanetary system, n evolves due to secular interactions. The long-term evolution of n can be approximated by a quasiperiodic expression

$$\begin{equation} n_x+in_y = \sum _k I_k e^{i(\nu _k t+\varphi _k)}, \end{equation} \tag{ 2 }$$

where the ν_k (sorted with increasing amplitudes I_k) are combinations of the fundamental frequencies g_j, s_j (Laskar 1990). For Saturn, ν₂ = s₈ = −0692 yr⁻¹ and I₂ = 0064. As the other terms have only very weak effects on the behavior of Saturn's spin axis (Hamilton & Ward 2004), one can retain this single term in the orbital precession, which makes the problem integrable (Colombo 1966; Henrard & Murigande 1987). The associated autonomous Hamiltonian, written in a moving reference frame related to the orbital plane, reads

$$\begin{equation} H=-\frac{\alpha }{2}({\mathbf {n}}\cdot {\mathbf {w}})^2-\nu ({\mathbf {k}}\cdot {\mathbf {w}}) \end{equation} \tag{ 3 }$$

with now $\mathbf {k}={}^t{(I,0,\sqrt{1-I^2})}$ and n = ^t(0, 0, 1). The equation of motion, obtained from ${d\mathbf {w}}/{dt}={\boldsymbol{\nabla }}_{\mathbf {w}}{H}\times \mathbf {w}$, is

$$\begin{equation} \frac{d\mathbf {w}}{dt} = -\alpha ({\mathbf {n}}\cdot {\mathbf {w}})\mathbf {n}\times \mathbf {w}-\nu \mathbf {k}\times \mathbf {w}. \end{equation} \tag{ 4 }$$

This system possesses four relative equilibriums named Cassini states for which the three axes w, n, and k are collinear, and a separatrix delineating three zones in the phase space (see Figure 1). Hereafter, we label the three zones after the Cassini state they contain. Saturn's spin axis coordinates in the orbital frame are

Zoom In Zoom Out Reset image size

$$\begin{equation} \mathbf {w}=\left(\begin{array}{@{}c@{}}\sin \epsilon \cos \psi \\ \sin \epsilon \sin \psi \\ \cos \epsilon \end{array}\right) \end{equation} \tag{ 5 }$$

with = 2673919 (Helled et al. 2009) and ψ = −31° (Hamilton & Ward 2004). As ψ ≠ 0, the system is not in a Cassini state. Given the coordinates of w, if −αcos  ∈ [ − 0.730, − 0.666] ''yr⁻¹ (Ward & Hamilton 2004) then Saturn's spin axis is in resonance around Cassini state 2 with a libration amplitude larger than 31° (see Figure 1(b)), else it is in circulation around either Cassini state 1 (see Figure 1(a)) or Cassini state 3. Literature gives three different values of Saturn's precession rate. Two of them, −074 ± 07 yr⁻¹ (French et al. 1993) and −075 ± 021 yr⁻¹ (Ward & Hamilton 2004), are compatible with a libration in zone 2, whereas the third one, −07542 ± 00002 yr⁻¹ (Helled et al. 2009), constrains Saturn's axis to circulate in zone 1. In the following, we study these two cases. In Case I, we use the precession constant given by Helled et al. (2009), and in Case II we set α such that −αcos  = s₈ = −0692 yr⁻¹. In our numerical integrations detailed below, we take into account the dependences of α in Saturn's semimajor axis and eccentricity.

We integrate the secular equations of motion derived from the Hamiltonian of Laskar & Robutel (1995) written up to degree 4 in inclinations and eccentricities. In order to fit to the present value of s₈ (Laskar et al. 2004), a small constant offset δs₈ = −000342 yr⁻¹ is added in the model. This offset was obtained by frequency analysis (Laskar 1990) of our analytical model (Table 1). For the class "MA," we consider only the last 3 AU migration of Neptune. When Neptune was closer to the Sun, the frequency s₈ was too large to have any effect on Saturn's axis. Migration is simulated by an additional force leading to the following exponential law,

$$\begin{equation} a(t) = a_0+\Delta a(1-e^{-t/\tau }), \end{equation} \tag{ 6 }$$

with Δa = +0.1, − 0.3, − 1.3, − 3 AU respectively for Jupiter, Saturn, Uranus, and Neptune. It is scaled from Minton & Malhotra (2009) and it is in agreement with the full integration of Tsiganis et al. (2005). In the same way, for the class "DE" we apply an external force that gives a long-term exponential evolution of Neptune eccentricity starting at e₀ = 0.3 and finishing at its current value. For both classes, we did integrations with constant Neptune inclination, and others with an exponential damping with the same τ. For each value of τ, an integration in the past is done to obtain initial conditions for the orbital coordinates. Saturn's initial obliquity is then set to ₀ = 1.5°.

We now look at the effect of the dissipation of the remaining primordial planetesimal belt. Following a suggestion of A. Morbidelli (2009, private communication), the mass m_K of the planetesimal belt in the class "MA" is estimated by energy conservation as follows. Initially, planetesimals are distributed following Morbidelli et al. (2004, Figure 1) and during Neptune migration, planetesimals move from their initial position to Uranus' orbit. This leads to m_K = 1.7 ± 0.1M_⊕. During planet migration with a planetesimal disk, we force an exponential decrease of the planetesimal belt mass with the same timescale as the semimajor axis one. To first order, the averaged effect of a planetesimal of mass m_i and semimajor axis a_i on Neptune's nodal precession rate is

$$\begin{equation} \delta s_8(m_i,a_i) = -\frac{n_N}{4}\left(\frac{m_i}{m_0}\right) \left(\frac{a_N}{a_i}\right)^2 b_{3/2}^{(1)}(a_N/a_i) \end{equation} \tag{ 7 }$$

with n_N and a_N being, respectively, Neptune mean motion and semimajor axis, and m₀ the mass of the Sun. We model the planetesimal belt by a single annulus with semimajor axis a_K such that

$$\begin{equation} \delta s_8(m_k,a_k)=\sum _i \delta s_8(m_i,a_i). \end{equation} \tag{ 8 }$$

Using the mass distribution of Morbidelli et al. (2004), we found a_k = 60 ±  5AU. We run numerical integrations with and without a 2 Earth mass primordial planetesimal belt and found that the constraints on the migration timescales were unchanged. In the following, we give only the results of our integrations without a planetesimal disk.

In Case I, Saturn's spin axis circulates around Cassini state 1 with a large obliquity (Figure 1(a)). To show whether it is compatible with the Ward & Hamilton (2004) scenario or not, we did several numerical integrations without long-term evolution of Neptune inclination and enumerated those ending in zone 1 with an obliquity larger than or equal to 26.73919°. We first considered the class MA of the Nice model and varied the migration timescale τ from 100 Myr to 600 Myr every 10 Myr. Then, for each value of τ, we searched the range of the initial precession angle ψ for which the final state corresponds to our criterion. We found only three values for τ satisfying the criterion: τ ∈ {170, 180, 190} Myr. In each case, the range of possible values for ψ is extremely small $\Delta _\psi \lesssim 10^{-5}\deg$ (Figure 2, Ia, Ib, Ic). Thus assuming an equiprobable initial phase, the probability to find Saturn in its current state through this mechanism is less than 3 × 10⁻⁸ for any of the three selected τ. With a DE migration type, the widths of the initial longitude intervals Δ_ψ are identical. The only changes are in the values of the migration timescale τ leading to the large obliquity circulation state: τ ∈ {150, 260, 290, 310, 320} Myr. We discuss the implications of these results in the conclusion.

Zoom In Zoom Out Reset image size

In Case II, Saturn's spin axis is presently in resonance with Neptune's ascending node. In that case, planet migration must be slow enough for the capture to occur, but if it is too slow, then the evolution becomes adiabatic and the libration amplitude is too small (less than 31°; Figure 2, IIc). This latter constraint disappears if the precapture obliquity is larger than 45 (Ward & Hamilton 2004). We performed 2100 integrations for each of the two migration types MA and DE, τ going from 10 to 600 Myr every 10 Myr and ψ between 0° and 350° every 10°. The results are summarized in Figure 3, MA type in gray and DE type in black. Probabilities are now significant and reach 1 for a few timescales. We see clear lower limits, τ ⩾ 90 Myr (resp. τ ⩾ 170 Myr) for the MA (resp. DE) migration type. The difference in the results between the two Nice model classes comes mainly from the different dependence of the semimajor axis and the eccentricity on Neptune's regression frequency. In all these integrations, Neptune's inclination does not undergo long-term evolution. However, the amplitude I₂ of Saturn orbital quasiperiodic motion (2) is proportional to Neptune inclination. In Section 5, we show that the higher the inclination amplitude is, the faster a planet can be tilted. We thus studied the minimum timescale, for which Saturn's axis ends in zone 2 with a libration amplitude larger than 31°, as a function of Neptune's initial inclination (Figure 4, bold curves). In both migration classes, τ_min decreases rapidly to ≈20 Myr when Neptune's initial inclination increases to 4° and then it decreases slowly down to ≈7 Myr when Neptune's inclination goes to 10°.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this section, we compute analytically the minimal time required to tilt a planet as a function of its inclination I(t). We give also an algorithm to check whether Saturn can be tilted or not for a given migration.

We call θ the inclination of a planet equator relative to the invariant plane. As Saturn current inclination is small relative to its obliquity , the two angles θ and are similar. Let Φ and Ω be the longitude of the ascending node of the equator and of the orbit in the invariant plane. We have

$$\begin{eqnarray} &&\mathbf {w}=\left(\begin{array}{@{}c@{}}\sin \theta \sin \Phi\\ -\sin \theta \cos \Phi\\ \cos \theta\end{array}\right) _{(\mathbf {i},\mathbf {j},\mathbf {k})}\nonumber\\ &&\mathbf {n}=\left(\begin{array}{@{}c@{}} \sin I\sin \Omega\\ -\sin I\cos \Omega \\ \cos I\end{array}\right)_{(\mathbf {i},\mathbf {j},\mathbf {k})} \end{eqnarray} \tag{ 9 }$$

and

$$\begin{equation} \begin{array}{@{}l@{}} {\mathbf {n}}\cdot {\bf {w}} = \sin \theta \sin I\cos (\Phi -\Omega)+\cos \theta \cos I,\\ {(\bf {n}\times \bf {w})}\cdot {\bf {k}} = -\sin \theta \sin I\sin (\Phi -\Omega). \end{array} \end{equation} \tag{ 10 }$$

Thus, from (1),

$$\begin{equation} \frac{d\cos \theta }{dt} = \alpha \sin \theta \cos \theta \sin I\cos I [1+\gamma \cos (\Phi -\Omega)]\sin (\Phi -\Omega) \end{equation} \tag{ 11 }$$

where γ = tan θtan I can vary from 0 to infinity depending on the value of the obliquity. We now choose Ω that maximizes this time derivative as a function of I and θ. Doing so, we ensure that it is not possible to have a faster evolution of the equator inclination θ. This leads to

$$\begin{equation} \cos (\Phi -\Omega)=\frac{-1+\sqrt{1+8\gamma ^2}}{4\gamma }. \end{equation} \tag{ 12 }$$

Substituting this expression in (11) gives the maximal speed Θ_max such that dθ/dt ⩽ Θ_max

$$\begin{equation} {\Theta _{\rm max}}=\frac{\alpha \sqrt{2}}{16}\left|\cos \theta \sin 2I\right| \frac{\left(\sqrt{1+8\gamma ^2}+3\right)^{3/2}}{\left(\sqrt{1+8\gamma ^2}+1\right)^{1/2}}. \end{equation} \tag{ 13 }$$

After some calculus, it can be shown that Θ_max is an increasing function of tan I. Thus, if the only constraint on the orbit inclination amplitude is an upper limit I_max < π/2, the fastest evolution is obtained for I = I_max. In two asymptotic cases, the expressions of Θ_max are simpler. For I ≪ |π/2 − θ| or θ ≪ |π/2 − I|, we have |γ| ≪ 1 and thus

$$\begin{equation} {\Theta _{\rm max}}\approx \frac{\alpha }{2}\left|\cos \theta {\sin 2I}\right|. \end{equation} \tag{ 14 }$$

In the other case, if |θ − π/2| ≪ I or |I − π/2| ≪ θ, the parameter γ is arbitrarily large and (13) becomes

$$\begin{equation} {\Theta _{\rm max}}\approx \frac{\alpha }{2}\sin \theta \sin ^2I. \end{equation} \tag{ 15 }$$

Using the approximation for small angles |γ| ≪ 1 (14), the minimum time t_min required to bring θ from 0 to θ_end at constant inclination amplitude I is

$$\begin{equation} {t_{\rm min}} = -\frac{2}{\alpha \sin 2I} \ln \left|\tan \left(\frac{\pi }{4}-\frac{\theta _{\rm end}}{2}\right)\right|. \end{equation} \tag{ 16 }$$

In Saturn's case, the amplitude of the mode responsible for the tilt is I₂. Whenever Neptune's inclination is less than 10°, I₂ remains below 09, γ ⩽ 8 × 10⁻³ and Expressions (14), (16) are valid. This minimum time t_min decreases with the inclination amplitude I. For example, with I = I₂ = 0.064°, Equation (16) gives t_min = 105 Myr in Case I and t_min = 115 Myr in Case 2.

From this study, it is possible to check whether Saturn's axis can be tilted or not for a given migration. Let θ₂ be the value of θ at the Cassini state 2. θ₂ and its time derivative $\Theta _2=\dot{\theta }_2$ are functions of orbital parameters through I(t), ν(t) = s₈(t) and α(t). During a tilt, θ oscillates around the increasing θ₂ and $\Theta =\dot{\theta }$ reads

$$\begin{equation} \Theta = \Theta _2 + A{\omega _{\rm lib}}\sin ({\omega _{\rm lib}}t+\varphi) \end{equation} \tag{ 17 }$$

where A and φ are, respectively, the libration amplitude and a phase, and ω_lib is the libration amplitude given by Hamilton & Ward (2004)

$$\begin{equation} {\omega _{\rm lib}}= \sqrt{-\alpha \nu \sin \theta \sin I}. \end{equation} \tag{ 18 }$$

For a given migration, Saturn's axis can tilt if and only if there exist A and ϕ such that Θ ⩽ Θ_max (14 and 17) during all the evolution. Replacing θ by θ₂ in (14 and 18), one obtains a criterion that depends only on orbital parameters. We applied this criterion on the systems studied in Section 4. For each value of Neptune's initial inclination, we integrated once the system with a given τ. Then, we rescaled the derivatives Θ for different value of τ until the criterion is verified. The resulting values of τ_min are displayed in Figure (4, thin curves).

First of all, we see that the Helled et al. (2009) precession constant is incompatible with the Ward & Hamilton (2004) scenario. This is a robust result. Helled et al. (2009) obtained Saturn's precession constant from an empirical model of its internal structure. They used Saturn mass, radius, and gravitational coefficients J₂, J₄, and J₆ to fit a density profile represented by a sixth-degree polynomial. From this density profile they derived the normalized axial moment of inertia γ directly related to the precession constant. Our results suggest, rather, considering γ as an additional independent parameter to better constrain Saturn's interior. If Saturn is actually in libration in zone 2 then 0.2257 < γ < 0.2438 (Ward & Hamilton 2004).

Assuming the Hamilton & Ward (2004) precession constant, Saturn's spin axis is likely to evolve toward a libration in zone 2 whatever the migration class is as long as the timescale τ is sufficiently large. We found a strong dependence between the minimum timescale τ_min and Neptune's inclination. Thus, an external constraint on the speed limit of Neptune migration may also constrain its inclination. For instance, the upper boundary obtained by Murray-Clay & Chiang (2005) is τ ⩽ 10 Myr. In that case, our results show that under the hypothesis of Section 3, the initial inclination of Neptune's orbit must have been larger than 7°. On the other side, in all our studied cases, the minimum timescale must be at least greater than 7 Myr, whereas Minton & Malhotra (2009) found τ ≲ 0.5 Myr. This contradiction may be raised if one considers different evolution laws for the semimajor axes, eccentricities, and/or inclinations. In that scope, we have given in Section 5 an algorithm to know if Saturn can be tilted for any migration law.