TILTING JUPITER (A BIT) AND SATURN (A LOT) DURING PLANETARY MIGRATION

Consider a planet revolving about the Sun and rotating with angular velocity ω about an instantaneous spin axis characterized by a unit vector ${\boldsymbol{s}}$. With solar gravitational torques applied on the planet, ω remains constant, but ${\boldsymbol{s}}$ evolves in the inertial space according to (e.g., Colombo 1966)

$$\begin{eqnarray}&&\frac{d{\boldsymbol{s}} }{dt}=-\alpha ({\boldsymbol{c}} \cdot {\boldsymbol{s}} )({\boldsymbol{c}} \times {\boldsymbol{s}} ),\end{eqnarray} \tag{ 1 }$$

where ${\boldsymbol{c}}$ denotes a unit vector along the orbital angular momentum, and

$$\begin{eqnarray}&&\alpha =\frac{3}{2}\frac{G{{M}_{\odot }}}{\omega \;{{b}^{3}}}\frac{{{J}_{2}}+q}{\lambda +l}\end{eqnarray} \tag{ 2 }$$

is the precession constant of the planet. Here, G is the gravitational constant, ${{M}_{\odot }}$ is the mass of the Sun, $b=a\sqrt{1-{{e}^{2}}}$, where a and e are the orbital semimajor axis and eccentricity, J₂ is the quadrupole coefficient of the planetary gravity field, and $\lambda =C/M{{R}^{2}}$ is the planetary moment of inertia C normalized by a standard factor MR², where M is planet's mass and R its reference radius (to be also used in the definition of J₂). The term $q=\frac{1}{2}{{\sum }_{j}}({{m}_{j}}/M){{({{a}_{j}}/R)}^{2}}$ is an effective, long-term contribution to the quadrupole coefficient due to the massive, close-in regular satellites with masses m_j and planetocentric distances a_j, and $l={{\sum }_{j}}({{m}_{j}}a_{j}^{2}{{n}_{j}})/(M{{R}^{2}}\omega )$ is the angular momentum content of the satellite orbits (n_j denotes their planetocentric mean motion) normalized to the characteristic value of the planetary rotational angular momentum.

For both Jupiter and Saturn, q is slightly dominant over J₂ in the numerator of the last term in Equation (2), while l is negligible in comparison to λ in the denominator. Spacecraft observations have been used to accurately determine the values of J₂, q, and l. On the other hand, λ cannot be derived in a straightforward manner from observations, because it depends on the structure of planetary interior. Using models of planetary interior, Helled et al. (2011) determined that Jupiter's ${{\lambda }_{{\rm J}}}$ is somewhere in the range between 0.2513 and 0.2529 (here rescaled for the equatorial radius R = 71,492 km of the planet). This would imply Jupiter's precession constant ${{\alpha }_{{\rm J}}}$ to be in the range between 2.754 and 2.772 arcsec yr⁻¹. These values are smaller than the ones considered in Ward & Canup (2006), if their proposed small angular distance between Jupiter's pole and Cassini state C₂ is correct.

Similarly, Helled et al. (2009) determined Saturn's precession constant ${{\alpha }_{{\rm S}}}$ to be in the range between 0.8443 and 0.8447 arcsec yr⁻¹. Again, these values differ from those inferred by Ward & Hamilton (2004) and Hamilton & Ward (2004), if a resonant confinement of the Saturn's spin axis in their proposed scenario is true. This difference has been noted and discussed in Boué et al. (2009). If Helled's values are correct, Saturn's spin axis cannot be presently locked in the resonance with s₈. However, it seems possible that of ${{\alpha }_{{\rm S}}}$ derived in Helled et al. (2009) may have somewhat larger uncertainty than reflected by the formal range of the inferred λ values. Also, the interpretation is complicated by the past orbital evolution of planetary satellites which may have also contributed to changes of ${{\alpha }_{{\rm S}}}$ (and ${{\alpha }_{{\rm J}}}$). For these reason, and in the spirit of previous studies (e.g., Ward & Hamilton 2004; Ward & Canup 2006; Boué et al. 2009), here we consider a wider range of the precession constant values α for both Jupiter and Saturn.

Assuming α constant for a moment, a difficult element preventing a simple solution of Equation (1) is the time evolution of ${\boldsymbol{c}}$. This is because mutual planetary interactions make their orbits precess in space on a characteristic timescale of tens of thousands of years and longer. In addition, during the early phase of planetary evolution, the precession rates of planetary orbits may have been faster due to the gravitational torques from a massive population of planetesimals in the trans-Neptunian disk. In the Keplerian parametrization of orbits, the unit vector ${\boldsymbol{c}}$ depends on the inclination I and longitude of node Ω, such that ${{{\boldsymbol{c}} }^{T}}=({\rm sin} I{\rm sin} {\Omega },-{\rm sin} I{\rm cos} {\Omega },{\rm cos} I)$. Traditionally, the difficulties with the time dependence in Equation (1) are resolved using a transformation to a reference frame fixed on an orbit, where ${{{\boldsymbol{c}} }^{T}}=(0,0,1)$.

The transformation from the inertial to orbit coordinate frames is achieved by applying a 3-1-3 rotation sequence with the Eulerian angles $({\Omega },I,-{\Omega })$. This transforms Equation (1) to the following form

$$\begin{eqnarray}&&\frac{d{\boldsymbol{s}} }{dt}=-[\alpha ({\boldsymbol{c}} \cdot {\boldsymbol{s}} ){\boldsymbol{c}} +{\boldsymbol{h}} ]\times {\boldsymbol{s}} ,\end{eqnarray} \tag{ 3 }$$

where now the planetary spin vector ${\boldsymbol{s}}$ is expressed with respect to the orbit frame, and ${\boldsymbol{c}}$ is now a unit vector along the z-axis. In effect, the time dependence has been moved to the vector quantity ${{{\boldsymbol{h}} }^{T}}=({\mathcal{A}},{\mathcal{B}},-2{\mathcal{C}})$ with

$$\begin{eqnarray}&&{\mathcal{A}}={\rm cos} {\Omega }\;\dot{I}-{\rm sin} I{\rm sin} {\rm \ \Omega\ \dot{\ \Omega\ }},\end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray}&&{\mathcal{B}}={\rm sin} {\Omega }\;\dot{I}+{\rm sin} I{\rm cos} {\Omega }\;{\rm \dot{\ \Omega\ }},\end{eqnarray} \tag{ 4b }$$

$$\begin{eqnarray}&&{\mathcal{C}}={{{\rm sin} }^{2}}I/2\;{\rm \dot{\ \Omega\ }},\end{eqnarray} \tag{ 4c }$$

and the over-dots mean time derivative.

A further development consists in introducing complex and non-singular parameter $\zeta ={\rm sin} (I/2){\rm exp} (\imath {\Omega })$ ($\imath =\sqrt{-1}$ is the complex unit) that replaces I and Ω. First order perturbation theory for quasi-circular and near-coplanar orbits indicates ζ for each planet can be expressed using a finite number of the Fourier terms with the s_i frequencies uniquely dependent on the orbital semimajor axes and masses of the planets, and amplitudes set by the initial conditions (e.g., Brouwer & Clemence 1961). In the models from nonlinear theories or numerical integrations, ζ can still be represented by the Fourier expansion $\zeta =\sum {{A}_{i}}{\rm exp} [\imath ({{s}_{i}}t+{{\phi }_{i}})]$ (e.g., Applegate et al. 1986; Laskar 1988), with the linear terms having typically the largest amplitudes A_i.

As discussed in Section 1, terms with present frequencies ${{s}_{7}}\simeq -2.985$ and ${{s}_{8}}\simeq -0.692$ arcsec yr⁻¹ are of a particular importance in this work. The s₇-term is the largest in Uranus' ζ representation, and the s₈-term is the largest in Neptune's ζ representation, and these terms also appear, though with smaller amplitudes, in the ζ variable of Jupiter and Saturn, because mutual planetary interactions enforce all fundamental frequency terms to appear in all planetary orbits. In terms of ζ, Equations (4a)–(4c) become

$$\begin{eqnarray}&&{\mathcal{A}}+\imath \;{\mathcal{B}}=\frac{2}{\sqrt{1-\zeta \bar{\zeta }}}\left( \frac{d\zeta }{dt}-\imath \zeta \;{\mathcal{C}} \right),\end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray}&&{\mathcal{C}}=\frac{1}{2\imath }\left( \bar{\zeta }\;\frac{d\zeta }{dt}-\zeta \;\frac{d\bar{\zeta }}{dt} \right),\end{eqnarray} \tag{ 5b }$$

where $\bar{\zeta }$ is complex conjugate to ζ. For small inclinations, relevant to this work, we therefore find that ${\mathcal{A}}+\imath \;{\mathcal{B}}\simeq 2(d\zeta /dt)$, and ${\mathcal{C}}\simeq 0$.

Another important aspect is of the problem is that Equation (3) derives from a Hamiltonian

$$\begin{eqnarray}&&{\mathcal{H}}({\boldsymbol{s}} ;t)=\frac{\alpha }{2}{{({\boldsymbol{c}} \cdot {\boldsymbol{s}} )}^{2}}+{\boldsymbol{h}} \cdot {\boldsymbol{s}} ,\end{eqnarray} \tag{ 6 }$$

such that

$$\begin{eqnarray}&&\frac{d{\boldsymbol{s}} }{dt}=-{{\nabla }_{{\boldsymbol{s}} }}{\mathcal{H}}\times {\boldsymbol{s}} .\end{eqnarray} \tag{ 7 }$$

This allowed Breiter et al. (2005) to construct an efficient Lie-Poisson integrator for a fast propagation of the secular evolution of planetary spins. In Section 3, we will use the leap-frog variant of the Hamiltonian's LP2 splitting from Breiter et al. (2005). To propagate ${\boldsymbol{s}}$ through a single integration step, Breiter et al. (2005) method requires that the orbital semimajor axis, eccentricity and ${\boldsymbol{c}}$ are provided at times corresponding to the beginning and end of the step. These values can be supplied from an analytic model of orbit evolution, or be directly obtained from a previous numerical integrations of orbits where a, e and ${\boldsymbol{c}}$ were recorded with a conveniently dense sampling.

The Hamiltonian formulation in Equation (6) allows us to introduce several important concepts of the Cassini dynamics. A long tradition in astronomy is to represent ${\boldsymbol{s}}$ with obliquity and precession angle ψ such that ${{{\boldsymbol{s}} }^{T}}=({\rm sin} \varepsilon {\rm sin} \psi ,{\rm sin} \varepsilon {\rm cos} \psi ,{\rm cos} \varepsilon )$. The benefit of this parametrization is that the unit spin vector is expressed using only two variables. The drawback is that resulting equations are singular when $\varepsilon =0$.

The conjugated momentum to ψ is $X={\rm cos} \varepsilon$. The Hamiltonian is then (e.g., Laskar & Robutel 1993)

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {\mathcal{H}}(X,\psi ;t) & = & \frac{\alpha }{2}{{X}^{2}}-2{\mathcal{C}}X \\ {} & {} & +\sqrt{1-{{X}^{2}}}({\mathcal{A}}{\rm sin} \psi +{\mathcal{B}}{\rm cos} \psi ). \\ \end{array}\end{eqnarray} \tag{ 8 }$$

For the low-inclination orbits, ${\mathcal{C}}$ is negligible, while ${\mathcal{A}}$ and ${\mathcal{B}}$ are expanded in the Fourier series with the same frequency terms as those appearing in ζ itself. A model of fundamental importance, introduced by G. Colombo (Colombo 1966; Henrard & Murigande 1987; Ward & Hamilton 2004), is obtained when only one Fourier term in ${\mathcal{A}}$ and ${\mathcal{B}}$ is considered. This Colombo model obviously serves only as an approximation of the complete spin axis evolution, since all other Fourier terms in ${\mathcal{A}}$ and ${\mathcal{B}}$ act as a perturbation. Nevertheless, the Colombo model allows us to introduce several important concepts that are the basis of the discussion in Section 3.

In the Colombo model, the orbital inclination is fixed and the node precesses with a constant frequency. Put in a compact way, $\zeta =A{\rm exp} [\imath (st+\phi )]$ is the single Fourier term, and $A={\rm sin} I/2$ and ${\Omega }=st+\phi$. Transformation to new canonical variables $X^{\prime} =-X$ and $\varphi =-(\psi +{\Omega })$, and scaling with the nodal frequency s, results in a time-independent Hamiltonian

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {\mathcal{H}}(X^{\prime} ,\varphi ) & = & \kappa X{{^{\prime} }^{2}}-{\rm cos} IX^{\prime} \\ {} & {} & +{\rm sin} I\sqrt{1-X{{^{\prime} }^{2}}}{\rm cos} \varphi , \\ \end{array}\end{eqnarray} \tag{ 9 }$$

where $\kappa =\alpha /(2s)$ is a dimensionless parameter. Note that the orbit-plane angle φ is measured from a reference direction that is 90° ahead of the ascending node. The general structure of the phase flow of solutions ${\mathcal{H}}(X^{\prime} ,\varphi )=C$, with C constant, derives from the location of the stationary points. Depending on the parameter values $(\kappa ,I)$, there are either two ($|\kappa |\lt {{\kappa }_{\star }}$) or four ($|\kappa |\gt {{\kappa }_{\star }}$) such stationary solutions (called the Cassini states). The critical value of κ reads (e.g., Henrard & Murigande 1987)

$$\begin{eqnarray}&&{{\kappa }_{\star }}(I)=\frac{1}{2}{{\left( {{{\rm sin} }^{2/3}}I+{{{\rm cos} }^{2/3}}I \right)}^{3/2}}.\end{eqnarray} \tag{ 10 }$$

Therefore, for small I, ${{\kappa }_{\star }}\simeq \frac{1}{2}$. The stationary solutions are located at $\varphi =0{}^\circ$ or $\varphi =180{}^\circ$ meridians in the orbital frame, and have obliquity values given by a transcendental equation

$$\begin{eqnarray}&&\kappa {\rm sin} 2\varepsilon =-{\rm sin} (\varepsilon \mp I),\end{eqnarray} \tag{ 11 }$$

with the upper sign for $\varphi =0{}^\circ$, and the lower sign for $\varphi =180{}^\circ$.

In the present work, we are mainly interested in the Cassini state C₂ located at $\varphi =0{}^\circ$. For Jupiter, the C₂ state related to frequency $s={{s}_{7}}$ is subcritical since $|\kappa |\lt \frac{1}{2}$ for all estimates of Jupiter's precession constant found in the literature. Only two Cassini states exist in this regime, and ${\boldsymbol{s}}$ must circulate about C₂. In the case of Saturn, $|\kappa |\gt {{\kappa }_{\star }}\simeq \frac{1}{2}$ for $s={{s}_{8}}$, four Cassini states exist in this situation, and ${\boldsymbol{s}}$ was suggested to librate in the resonant zone about C₂. The configuration of vectors in C₂ can be inferred from Equation (11). If the inclination is significantly smaller than the obliquity , we have $\alpha {\rm cos} \varepsilon \simeq -s$. Since the term on the left-hand side of this relation is the precession frequency of the planet's spin (see Equation (1)), we find that the spin and orbit vectors will co-precess with the same rate about the normal vector to the inertial frame. Small resonant librations about C₂ would reveal themselves by small departures of the spin vector from this ideal state. The maximum width ${\Delta }\varepsilon$ of the resonant zone in obliquity at the $\varphi =0{}^\circ$ meridian can be obtained using an analytic formula (e.g., Henrard & Murigande 1987)

$$\begin{eqnarray}&&{\rm sin} \frac{{\Delta }\varepsilon }{2}=\frac{1}{|\kappa |}\sqrt{\frac{{\rm sin} 2I}{{\rm sin} 2{{\varepsilon }_{4}}}},\end{eqnarray} \tag{ 12 }$$

where ${{\varepsilon }_{4}}$ is the obliquity of the unstable Cassini state C₄ (a solution of Equation (11) at $\varphi =180{}^\circ$ having the intermediate value of the obliquity).

Figure 1, top panel, shows how the location of the Cassini states and the resonance width ${\Delta }\varepsilon$ depend on κ, which is the fundamental parameter that changes during planetary migration. For sake of this example we assumed orbital inclination $I=0\buildrel{\circ}\over{.} 5$ (note that the overall structure remains similar for even smaller inclination values considered in the next section, but would be only less apparent in the figure). The C₁ and C₄ stationary solutions bifurcate when $\kappa ={{\kappa }_{\star }}$ at a non-zero critical obliquity value ${{\varepsilon }_{\star }}={\rm atan}({{{\rm tan} }^{1/3}}I)$ (e.g., Henrard & Murigande 1987; Ward & Hamilton 2004). Note that ${\Delta }\varepsilon$ is significant in spite of a very small value of the inclination, which manifests through its dependence on a square root of ${\rm sin} 2I$ in (12). The bottom panels show examples of the phase portraits ${\mathcal{H}}(X^{\prime} ,\varphi )=C$ for both sub-critical $|\kappa |\lt {{\kappa }_{\star }}$ and super-critical $|\kappa |\gt {{\kappa }_{\star }}$ cases.

Zoom In Zoom Out Reset image size

NM12 reported the results of nearly 10⁴ numerical integrations of planetary instability, starting from hundreds of different initial configurations of planets that were obtained from previous hydrodynamical and N-body calculations. The initial configurations with the 3:2 Jupiter–Saturn mean motion resonance were given special attention, because Jupiter and Saturn, radially migrating in the gas disk before its dispersal, should have become trapped into their mutual 3:2 resonance (e.g., Masset & Snellgrove 2001; Morbidelli & Crida 2007; Pierens & Nelson 2008). They considered cases with four, five and six initial planets, where the additional planets were placed onto resonant orbits between Saturn and the inner ice giant, or beyond the orbit of the outer ice giant. The integrations included the effects of the transplanetary planetesimal disk. NM12 experimented with different disk masses, density profiles, instability triggers, etc., in an attempt to find solutions that satisfy several constraints, such as the orbital configuration of the outer planets, survival of the terrestrial planet, and the distribution of orbits in the asteroid belt.

NM12 found the dynamical evolution in the four planet case was typically too violent if Jupiter and Saturn start in the 3:2 resonance, leading to ejection of at least one ice giant from the Solar System. Planet ejection can be avoided if the mass of the transplanetary disk of planetesimals was large (${{M}_{{\rm disk}}}\gtrsim$ $50\;{{M}_{{\rm Earth}}}$, where ${{M}_{{\rm Earth}}}$ is the Earth mass), but such a massive disk would lead to excessive dynamical damping (e.g., the outer planet orbits become more circular then they are in reality) and to smooth migration that violates constraints from the survival of the terrestrial planets, and the asteroid belt. Better results were obtained when the Solar System was assumed to have five giant planets initially and one ice giant, with mass comparable to that of Uranus or Neptune, was ejected into interstellar space by Jupiter (Nesvorný 2011; Batygin et al. 2012). The best results were obtained when the ejected planet was placed into the external 3:2 or 4:3 resonances with Saturn and ${{M}_{{\rm disk}}}\simeq 20$ ${{M}_{{\rm Earth}}}$. The range of possible outcomes was rather broad in this case, indicating that the present Solar System is neither a typical nor expected result for a given initial state.

The most relevant feature of the NM12 models for this work is that the planetary migration happens in two stages (see Figure 2). During the first stage, that is before the instability happens, Neptune migrates into the outer disk at $\simeq 20-30$ AU. The migration is relatively fast during this stage, because the outer disk still has a relatively large mass. We analyzed several simulations from NM12 and found that Neptune's migration can be approximately described by an exponential with the e-folding timescale $\tau \simeq 10$ Myr for ${{M}_{{\rm disk}}}=20$ ${{M}_{{\rm Earth}}}$ and $\tau \simeq 20$ Myr for ${{M}_{{\rm disk}}}=15$ ${{M}_{{\rm Earth}}}$. The instability typically happens in the NM12 models when Neptune reaches $\simeq 28$ AU. The main characteristic of the instability is that planetary encounters occur, mainly between the extra ice giant and all other planets. The instability typically lasts $\sim {{10}^{5}}$ years and terminates when the extra ice giant is ejected from the solar system by Jupiter. The second stage of migration starts after that. The migration of Neptune is much slower during this period, because the outer disk is now very much depleted. From simulations in NM12 we find that $\tau \simeq 30$ Myr for ${{M}_{{\rm disk}}}=20$ ${{M}_{{\rm Earth}}}$ and $\tau \simeq 50$ Myr for ${{M}_{{\rm disk}}}=15$ ${{M}_{{\rm Earth}}}$. Moreover, rather then being precisely exponential, the migration slows down relative to an exponential with fixed τ, such as, effectively, the very late stages have larger τ values ($\tau \sim 100$ Myr) than the ones immediately following the instability. Uranus accompanies the migration of Neptune on timescales similar to those mentioned above.

Zoom In Zoom Out Reset image size

The frequencies s₇ and s₈, which are the most relevant for this work, are initially somewhat higher than the precession constants of Jupiter and Saturn, mainly because of the torques from the outer disk. The extra term from the third ice giant initially located at $\simeq 10$ AU has much faster frequency than the precession constants and does not interfere with the obliquity of Jupiter and Saturn during the subsequent evolution. The s₇ and s₈ frequencies slowly decrease during both stages. Their e-folding timescales may slightly differ from the migration e-folding timescales mentioned above, due to the nonlinearity of the dependence of the secular frequencies on the semimajor axis of planets. Our tests show that they are about (90–95)% of the e-folding timescales of planetary semimajor axes.

From analyzing the behavior of frequencies in different simulations we found that s₈ should cross the value of Jupiter's precession constant during the first migration stage, that is before the instability. The main characteristic of this crossing is that the planetary orbits are very nearly coplanar during this stage. The amplitude I in Equation (9) should thus be very small. It is not known exactly, however, how small. In Section 3.2, we consider amplitudes down to 0005 (about 10 times smaller than the present value of I₅₈, where I₅₈ is the amplitude of the 8th frequency term in Jupiter's orbit), and show that the effects on Jupiter's obliquity are negligible if the amplitudes were even lower. The second characteristic of the first migration stage is that the evolution of s₈ happens on a characteristic timescale of $\simeq 10-20$ Myr. Since the total change of s₈ during this interval is several arcsec yr⁻¹, and the first stage typically lasts ∼10–20 Myr, the average rate of change is very roughly, as an order of magnitude estimate, $d{{s}_{8}}/dt\sim 0.1$ arcsec yr⁻¹ Myr⁻¹. The actual value of $d{{s}_{8}}/dt$ during crossing depends on several unknowns, including when exactly the crossing happens during the first stage. Also, the changes of s₈ could have been slower if the first stage lasted longer than in the NM12 simulations, as required if the instability occurred at the time of the Late Heavy Bombardment (e.g., Gomes et al. 2005). In Section 3.2, we will consider values in the range $0.005\lt d{{s}_{8}}/dt\lt 0.05$ arcsec yr⁻¹ Myr⁻¹, and show that the obliquity of Jupiter cannot be pumped up to its current value if $d{{s}_{8}}/dt\gt 0.05$ arcsec yr⁻¹ Myr⁻¹ (assuming that ${{I}_{58}}\lesssim 0\buildrel{\circ}\over{.} 05$).

Interesting effects on obliquities should happen during the second migration stage. First, the s₈ frequency reaches the value of the precession constant of Saturn. There are several differences with respect to the s₈ crossing of Jupiter's precession constant during the first stage (discussed above). The orbital inclinations of planets were presumably excited to their current values during the instability. Therefore, the amplitude I₆₈ should be comparable to its current value, ${{I}_{68}}\simeq 0\buildrel{\circ}\over{.} 064$, during the second stage. We see this happening in the NM12 simulations. First, there is a brief period during the instability, when the inclinations of all planets are excited by encounters with the ejected ice giant. The inclination of Neptune is modest, at most $\simeq 2{}^\circ$, and is rather quickly damped by the planetesimal disk. Also, the invariant plane of the solar system changes by ≃05 when the third ice giant is ejected during the instability. The final inclinations are of this order. The current amplitudes are ${{I}_{58}}\simeq 0\buildrel{\circ}\over{.} 066$ (s₈ term in Jupiter's orbit) and ${{I}_{68}}\simeq 0\buildrel{\circ}\over{.} 064$ (s₈ term in Saturn's orbit; see e.g., Laskar 1988).

Another difference with respect the first stage is that the evolution of s₈ is much slower during the second stage. If, as indicated by the NM12 integration, s₈ changes by ∼1 arcsec yr⁻¹ in 100 Myr, then the average rate of change is very roughly $d{{s}_{8}}/dt\sim 0.01$ arcsec yr⁻¹ Myr⁻¹. The actual rate of change can be considerably lower than this during the very late times, when the effective τ was lower than during the initial stages. Finally, during the very last gasp of migration, the s₇ frequency should have approached the precession constant of Jupiter. We study this case in an adiabatic approximation when the rate of change of s₇ is much slower than any other relevant timescale. We find that the present obliquity of Jupiter can be excited by the interaction with the s₇ term only if the precession constant of Jupiter is somewhat larger than inferred by Helled et al. (2011), in accord with the results of Ward & Canup (2006).

Since s₈ remains larger than ${{\alpha }_{{\rm S}}}$ during the first stage, we do not expect any important effects on Saturn's obliquity during this stage. If Saturn's obliquity was low initially, it should have remained low in all times before the instability. We therefore focus on the case of Jupiter in this section. From the analysis of the NM12 numerical simulation in Section 3.1, we infer that the s₈ frequency crossed ${{\alpha }_{{\rm J}}}$ during the first stage. The values of $d{{s}_{8}}/dt$ and I₅₈ during crossing are not known exactly from the NM12 simulations, because they depend on details of the initial conditions. We therefore consider a range of values and determine how Jupiter's obliquity excitation depends on them. The results can be used to constrain future simulations of the planetary instability/migration.

We consider Colombo's model with only one Fourier term in ζ of Jupiter, namely that of the s₈ frequency.⁴ The inclination term I₅₈ in Jupiter's orbit was treated as a free parameter. The range of values was set to be between zero and roughly the current value of ≃0066. As we discussed in Section 3.1, it is reasonable to assume that I₅₈ was smaller than the current value, because the orbital inclination of planets should have been low in times before the instability. The value of ${{\alpha }_{{\rm J}}}$ was obtained by rescaling the present value to the the semimajor axis of Jupiter before the instability (${{a}_{{\rm J}}}\simeq 5.45$ AU). To do so we used Equation (2) and assumed ${{\alpha }_{{\rm J}}}\propto a_{{\rm J}}^{-3}$. No additional modeling of possible past changes of ${{\alpha }_{{\rm J}}}$, for instance due to satellite system evolution or planetary contraction, was implemented. The s₈ frequency was slowly decreased from a value larger than ${{\alpha }_{{\rm J}}}$ to a value smaller than ${{\alpha }_{{\rm J}}}$.

Motivated by the numerical simulations of the instability discussed in Section 3.1, we assumed the initial s₈ value of −4 arcsec yr⁻¹ and let it decrease to −1.2 arcsec yr⁻¹ by the end of each test. The rate of change, $d{{s}_{8}}/dt$, was treated as a free parameter. The integrations were carried for several tens of Myr for the highest assumed rates and up to several hundreds of Myr for the lowest rates. We recorded Jupiter's obliquity during the last 5 Myr of each run, and computed the mean value ${{\varepsilon }_{{\rm fin}}}$. Jupiter's initial spin axis was oriented toward the pole of the Laplacian plane. (The value of ${{\varepsilon }_{{\rm fin}}}$ reported in Figure 3 was averaged over all possible phases of the initial spin axis on the ${\mathcal{H}}(X^{\prime} ,\varphi )=C$ level curve, with C defined by ${\boldsymbol{s}}$ oriented toward the pole of the Laplacian plane.)

Zoom In Zoom Out Reset image size

Figure 3 shows the results. For most parameter combinations shown here the s₈ resonance swept over ${{\alpha }_{{\rm J}}}$ without having the ability to capture Jupiter's spin vector in the resonance. This happened because $d{{s}_{8}}/dt$ was relatively large and I₅₈ was relatively small, thus implying that the s₈ frequency crossed the resonant zone in a time interval that was shorter than the libration period. Captures occurred only in extreme cases (largest I₅₈ and smallest $d{{s}_{8}}/dt$). These cases ended up generating very large obliquity values of Jupiter and are clearly implausible. The plausible values of $d{{s}_{8}}/dt$ and I₅₈ correspond to the cases where Jupiter's obliquity was not excited at all, thus assuming that Jupiter obtained its present obliquity later, or was excited by up to $\simeq 3{}^\circ$. To obtain ${{\varepsilon }_{{\rm J}}}\simeq 3{}^\circ$, $d{{s}_{8}}/dt$ and I₅₈ would need to have values along the bold line labeled 3 in Figure 3, which extends diagonally in $d{{s}_{8}}/dt$ and I₅₈ space. An example of a case where the obliquity of Jupiter was excited to near $\simeq 3{}^\circ$ value is shown in Figure 4.

Zoom In Zoom Out Reset image size

We now turn our attention to the second stage, when the migration slowed down and the obliquities of Jupiter and Saturn should have suffered additional perturbations. At the beginning of stage, that is just after the time of the instability, the s₈ frequency is already lower than ${{\alpha }_{{\rm J}}}$, but still higher than ${{\alpha }_{{\rm S}}}$, while the s₇ frequency is higher than ${{\alpha }_{{\rm J}}}$. Since the s₇ and s₈ frequencies are slowly decreasing during the second stage, a possibility arises that Jupiter's obliquity was (slightly) excited when s₇ approached ${{\alpha }_{{\rm J}}}$ (e.g., Ward & Canup 2006), and that Saturn's obliquity was strongly excited by capture into the spin–orbit resonance with s₈ (e.g., Hamilton & Ward 2004; Ward & Hamilton 2004; Boué et al. 2009).

3.3.1. Jupiter

Ward & Canup (2006) suggested a possibility that Jupiter's present obliquity may be explained by the proximity of ${{\alpha }_{{\rm J}}}$ to the current value of the s₇ frequency. They showed that, if $\kappa ={{\alpha }_{{\rm J}}}/(2{{s}_{7}})$ is sufficiently close to the critical value from Equation (10), namely $\simeq \frac{1}{2}$ for small inclinations, the obliquity of the Cassini state C₂ may be significant. Thus, as s₇ adiabatically approached to ${{\alpha }_{{\rm J}}}$, Jupiter's obliquity may have been excited along. As a supportive argument for this scenario they pointed out that ${{\varphi }_{{\rm J}}}\simeq 0{}^\circ$, where the Cassini state 2 is located.

To test this possibility we run a suite of simulations, assuming an exponential convergence to the current value ${{s}_{7}}\simeq -2.985$ arcsec yr⁻¹. Specifically, we set ${{s}_{7}}(t)={{s}_{7}}\;+[{{s}_{7}}(0)-{{s}_{7}}]{\rm exp} (-t/\tau )$, where τ and ${{s}_{7}}(0)$ are parameters. The initial value ${{s}_{7}}(0)$ at the beginning of the second stage was obtained from the numerical simulation discussed in Section 3.1. Here we chose to use ${{s}_{7}}(0)=-3.5$ arcsec yr⁻¹, however we verified that the results are insensitive to this choice. The e-folding timescale τ depends on how slow or fast planets migrate. Given that the planetary migration is slow during the second stage, we chose $\tau =100-200$ Myr. This assures that the approach of ${{s}_{7}}(t)$ to ${{\alpha }_{{\rm J}}}$ is adiabatical. The amplitude I₅₇ is assumed to be constant and equal to its current value (${{I}_{57}}\simeq 0\buildrel{\circ}\over{.} 055$).

Two additional parameters need to be specified: (i) Jupiter's initial spin state, and (ii) ${{\alpha }_{{\rm J}}}$. As for (i), the results discussed in Section 3.2 indicate that Jupiter's obliquity may have remained near zero during the first stage, if I₅₈ was too small and/or $d{{s}_{8}}/dt$ was too fast, or could have been potentially excited to $\simeq 3{}^\circ$, if I₅₈ and $d{{s}_{8}}/dt$ combined in the right way. Therefore, here we treat the obliquity of Jupiter at the beginning of stage 2 as a free parameter. As for (ii), as we discussed in Section 1, the present value of ${{\alpha }_{{\rm J}}}$ somewhat uncertain. We therefore performed various simulations, where ${{\alpha }_{{\rm J}}}$ takes on a number of different values between 2.75 and 2.95 arcsec yr⁻¹. A similar approach has also been adopted by Ward & Canup (2006).

Figure 5 reports the results. The top panel shows how the obliquity ${{\varepsilon }_{{\rm C}2}}$ of the Cassini state C₂ depends on the assumed value of ${{\alpha }_{{\rm J}}}$. This is calculated from Equation (11). The trend is that ${{\varepsilon }_{{\rm C}2}}$ increases with ${{\alpha }_{{\rm J}}}$, because the larger values of ${{\alpha }_{{\rm J}}}$ correspond to a situation where the system is closer to the exact resonance with s₇. If ${{\alpha }_{{\rm J}}}\lt 2.8$ arcsec yr⁻¹, as inferred from models in Helled et al. (2011), ${{\varepsilon }_{{\rm C}2}}$ is too small to significantly contribute to ${{\varepsilon }_{{\rm J}}}$. This case would imply that Jupiter's present obliquity had to be acquired during the earlier stages and is possibly related to the non-adiabatic s₈ resonance crossing discussed in Section 3.2. If, on the other hand, ${{\alpha }_{{\rm J}}}\simeq 2.92-2.94$ arcsec yr⁻¹, Jupiter would owe its present obliquity to the proximity between ${{\alpha }_{{\rm J}}}$ and s₇. This would imply that the obliquity excitation during the first stage of planetary migration must have been minimal. Figures 3 and 5 express the joint constraint on the planetary migration also for the intermediate cases, where the present obliquity of Jupiter arose as a combination of both effects discussed here.

Zoom In Zoom Out Reset image size

Figure 6 illustrates the two limiting cases discussed above. In panel (a), we assumed that the parameters during the first migration stage were such that the obliquity of Jupiter was excited to its current value during the s₈ crossing (such as shown in Figure 4). Also, we set ${{\alpha }_{{\rm J}}}=2.77$ arcsec yr⁻¹, corresponding to the best theoretical value of Helled et al. (2011). The Cassini state C₂ corresponding to the s₇ term is only slightly displaced from the center of the plot, and does not significantly contribute to the present obliquity value. In panel (b) of Figure 6, we set ${{\alpha }_{{\rm J}}}=2.93$ arcsec yr⁻¹. This value implies that ${{\varepsilon }_{{\rm C}2}}\simeq 2\buildrel{\circ}\over{.} 6$. The present obliquity of Jupiter would then be in large part due to the "forced" term arising from the proximity of s₇. The initial excitation of Jupiter's obliquity during the first stage would have to be minor in this case.

Zoom In Zoom Out Reset image size

3.3.2. Saturn

In the case of Saturn, all action is expected to take place during the second stage of planetary migration. Insights gleaned from the numerical simulations discussed in Section 3.1 show that the s₈ frequency should have very slowly approached ${{\alpha }_{{\rm S}}}$, thus providing a conceptual basis for capture of Saturn's spin vector in a resonance with s₈ (e.g., Hamilton & Ward 2004; Ward & Hamilton 2004; Boué et al. 2009). To study this possibility, we assume that the I₆₈ amplitude was excited to its current value during the instability, and remained nearly constant during the second stage of planetary migration. This choice is motivated by the NM12 simulations, where the inclination of Neptune is never too large. Note that Boué et al. (2009) investigated the opposite case where Neptune's inclination was substantially excited during the instability and remained high when the resonance with s₈ occurred. This type of strong inclination excitation does not happen in the NM12 models.

Here we assume that the planetary migration was very slow during the second stage and parametrize ${{s}_{8}}(t)$ as ${{s}_{8}}(t)={{s}_{8}}+[{{s}_{8}}(0)-{{s}_{8}}]{\rm exp} (-t/\tau )$ with $\tau \geqslant 80$ Myr. The initial frequency value at the beginning of stage 2, ${{s}_{8}}(0)$, is estimated from the NM12 simulations. We find that ${{s}_{8}}(0)\simeq -1.3$ arcsec yr⁻¹, and use this value to set up the evolution of ${{s}_{8}}(t)$. We also assume a range of ${{\alpha }_{{\rm S}}}$ values. This has the following significance. As already pointed out by Boué et al. (2009), the best-modeling values of ${{\alpha }_{{\rm S}}}$ from Helled et al. (2009) are not compatible with a resonant location of Saturn's spin axis. This is because the Cassini state C₂ would be moved to a significantly larger obliquity value ($\geqslant 34{}^\circ$). So these values of ${{\alpha }_{{\rm S}}}$ would imply that Saturn's spin circulates about the Cassini state C₁. On the other hand, the significant obliquity of Saturn requires an increase when the s₈ value was crossing ${{\alpha }_{{\rm S}}}$ value, as schematically shown in the left panel of Figure 4 for Jupiter's obliquity during the phase 1. Boué et al. (2009) tested this scenario using numerical simulations and found it extremely unlikely: initial data of an insignificant measure have led this way to the current spin state of Saturn. Indeed, here we recover the same result with a less extensive set of numerical simulations.

Given the arguments discussed above we therefore tend to believe that the precession constant of Saturn may be somewhat smaller than the one determined by Helled et al. (2009). For instance, R. A. Jacobson (2015, personal communication) determined the Saturn precession from the Saturn's ring observations. The mean precession rate obtained by him is 0.725 arcsec yr⁻¹ (formal uncertainty of about 6%). This value would indicate ${{\alpha }_{{\rm S}}}$ in the range between 0.769 and 0.864 arcsec yr⁻¹. Both Ward & Hamilton (2004) and Boué et al. (2009) report other observational constraints of Saturn's pole precession that have comparably large uncertainty. We therefore sampled a larger interval of the ${{\alpha }_{{\rm S}}}$ values to make sure that all interesting possibilities are accounted for.

Our numerical simulations thus spanned a grid of two parameters: (i) ${{\alpha }_{{\rm S}}}$ discussed above, and (ii) τ, the e-folding timescale of the s₈ frequency that slowly changes due to residual migration of Neptune and depletion of the outer disk. The amplitude related to the s₈ term in the inclination vector ζ of Saturn is kept constant, namely ${{I}_{68}}=0\buildrel{\circ}\over{.} 064$. To keep number of tested free parameters low, we assumed initial orientation of Saturn's spin axis ${\boldsymbol{s}}$ to be near the pole of the invariable plane. Specifically, we set its obliquity to 01 in the reference frame of the s₈-frequency Fourier term in ζ. To prevent fluke results, we sampled 36 values of longitude φ of the Saturn's pole in the same reference frame, incrementing it from 0° by 10°. Each of the simulations covered a 1 Gyr timespan. We recorded Saturn's pole orientation during the last 150 Myr time interval. We specifically analyzed if it passes close to the current location of Saturn's pole, namely, ${{\varepsilon }_{8}}\simeq 27\buildrel{\circ}\over{.} 4$ and ${{\varphi }_{8}}\simeq -31\buildrel{\circ}\over{.} 4$ in the s₈-frequency reference frame (see Table 2 in Ward & Hamilton 2004). A numerical run was considered successful if the simulated Saturn's pole passed through a box of $\pm 0\buildrel{\circ}\over{.} 2$ in obliquity and $\pm 3{}^\circ$ in longitude φ around the planet's values $({{\varepsilon }_{8}},{{\varphi }_{8}})$ during the recorded 150 Myr time interval. Note that the libration period of Saturn's pole around Cassini state C₂, if captured in the spin–orbit resonance, is ∼(50–100) Myr, depending on the libration amplitude. This set our requirement for the timespan over which we monitored Saturn's pole position.

Figure 7 shows the results from this suite of runs. The shaded region shows correlated ${{\alpha }_{{\rm S}}}$-τ pairs that provided successful match to the Saturn's pole position. We note that no successful solutions were obtained for ${{\alpha }_{{\rm S}}}\gt 0.812$ arcsec yr⁻¹ and all successful solutions correspond to the capture in the resonance zone around the Cassini state C₂. The absence of low-probability solutions in which Saturn's pole would circulate about the Cassini state C₁ for larger ${{\alpha }_{{\rm S}}}$ may be related to the limited number of simulations performed. No solutions were also obtained for $\tau \gt 215$ Myr. This is because for such long e-folding timescales the resonant capture process would be strictly adiabatic and the simulated spin would not meet the condition of at least 30° libration amplitude (see discussion in Hamilton & Ward 2004; Ward & Hamilton 2004). The area occupied by successful solutions splits into two branches for ${{\alpha }_{{\rm S}}}\simeq 0.78$ arcsec yr⁻¹. This is because the obliquity of the Cassini state C₂ is ${{\varepsilon }_{{\rm C}2}}\simeq 27\buildrel{\circ}\over{.} 4$ for the critical value of ${{\alpha }_{{\rm S}}}$, and the solutions have a minimum libration amplitude of $\simeq 31{}^\circ$.

Zoom In Zoom Out Reset image size

Figure 8 shows two evolution paths of Saturn's spin vector obtained in two different simulations. As mentioned above, in both cases Saturn's spin state was captured into the spin–orbit resonance s₈, and remained in the resonance for the full length of the simulation. The final states of both simulations are a good proxy for the Saturn's present spin state. These two cases differ from each other principally because the path in panel (b) shows librations with a larger amplitude than the path in panel (a). Note some of the solutions, such as (b) here, may attain a significant librations amplitude. This is because with the corresponding values of $\tau \simeq 100$ Myr the evolution is not adiabatic and the librations amplitude is excited immediately after capture. Therefore, the complicated evolution histories proposed in Hamilton & Ward (2004) may not be needed.

Zoom In Zoom Out Reset image size