TILTING SATURN WITHOUT TILTING JUPITER: CONSTRAINTS ON GIANT PLANET MIGRATION

In the last two decades, considerable progress has been made in regards to understanding the processes that shaped the dynamical evolution of the solar system, and it seems the migration of the giant planets is a critical chapter in this story. Angular momentum exchange with a remnant planetesimal disk induces planetary migration, moving Jupiter slightly inward and Saturn, Uranus, and Neptune outwards (Fernandez & Ip 1984; Hahn & Malhotra 1999). The existence of many Kuiper Belt objects in orbital mean-motion resonances (MMRs) with Neptune can be explained by the outward migration of Neptune. If the migration of the planets was smooth, the orbital eccentricities of the resonant Kuiper Belt objects and the asymmetry in the resonant angle distribution of objects in the 2:1 resonance can be used to infer that Neptune has migrated by $\sim 7\;{\rm{AU}}$ on a timescale $\gtrsim 1\;{\rm{Myr}}$ and that the initial mass of the planetesimal disk is $\sim 50\;{M}_{\oplus },$ where ${M}_{\oplus }$ is the mass of the Earth (Malhotra 1995; Hahn & Malhotra 1999; Murray-Clay & Chiang 2005). However, this smooth migration model has difficulty reproducing several features of the solar system, and alternative scenarios have been invoked.

Thommes et al. (1999) proposed a scenario in which Uranus and Neptune grew as cores between Jupiter and Saturn, were scattered close to their current locations due to planet–planet gravitational interactions, and had their random velocities damped by interactions with a remnant planetesimal disk. In this way the challenge of forming the ice giants at their current locations, with long orbital timescales and low (estimated) solid surface density at the time of formation, could be neatly side-stepped. Their work provided the first strong numerical evidence that the outer solar system could have undergone considerable rearrangement from the time that the giant planets achieved their current masses.

The successor of this model for the configuration of the giant planets after the protoplanetary gas disk has dissipated was presented by Tsiganis et al. (2005) and is known as the Nice (or "classic Nice") model. It proposes that the giant planets were once in a considerably more compact configuration, in which Saturn started just inside the 2:1 MMR with Jupiter, and the ice giants were located at 11–13 and 13.5–$17\;{\rm{AU}}.$ When the 2:1 resonance was crossed by the divergent migration of Jupiter and Saturn caused by the gravitational interactions with planetesimals, the system became dynamically unstable, exciting eccentricities and inclinations, and scattering Uranus and Neptune to close to their present locations (in some cases, after switching their order). This scenario has many nice properties, including possibly explaining the Late Heavy Bombardment (Gomes et al. 2005) and Jupiter's Trojan population (Morbidelli et al. 2005), while typically preserving the regular satellite populations during the planetary encounters.

However, even granting the classic Nice model its successes, there were obvious questions about the history: how did the system arrive at such an initial configuration? The timescales for migration induced by interactions with the protoplanetary gas disk to drive the giant planets into the inner system are shorter than the formation times for the object (e.g., Tanaka et al. 2002). The recent revisions to the type I migration rate from a more sophisticated treatment of thermodynamics (e.g., Paardekooper et al. 2010) have complicated the situation, and it remains a challenge. Building upon work done by Masset & Snellgrove (2001) and Morbidelli & Crida (2007), a variant of the classic Nice model was proposed by Morbidelli et al. (2007). This resonant Nice model attempts to bridge the gap between the gas-rich phase (typically studied with one or two planets embedded in a gas disk modeled using a hydrodynamics code) and the gas-free phase (studied with many planets and planetesimals using an N-body code). The authors use the observation that convergent migration in the gas-rich phase can result in capture into MMR and that a 3:2 MMR between Jupiter and Saturn can prevent Jupiter from migrating inwards. They construct a number of stable configurations in which each giant planet is in resonance with its neighbors (with Saturn and the inner ice giant in the 3:2 MMR, and the ice giants in the 4:3 or 5:4 MMR). Subsequent migration due to interactions with the planetesimals breaks the resonances, and an instability resembling that of classic Nice occurs.

Recent studies on giant planet migration have tried to determine the actual evolution based on constraints from the secular architecture of the giant planets (Morbidelli et al. 2009), the current orbits of the terrestrial planets (Brasser et al. 2009) and the asteroid belt (Morbidelli et al. 2010). In essence, it is required that Jupiter scattered one of the ice giants outwards and thus migrated on a very short timescale to avoid both making the orbits of the terrestrial planets too eccentric and causing secular resonances to sweep through the asteroid belt. Unfortunately this scenario had a very low probability of success, which led Nesvorný (2011) and Batygin et al. (2012) to propose a scenario in which the solar system began with five giant planets locked in MMRs. Nesvorný & Morbidelli (2012) have studied a large number of five planet systems to try to reproduce the current solar system and match all the constraints, which they did with a probability of about 5% in the best cases.

One thing that has not yet been considered are the effects of migration on planetary obliquities (the angles between the spin axes of the planets and their respective orbit normals) as the tilt angles of the giant planets preserve information about migration and encounter history. The case of Saturn is particularly interesting because of a coincidence between the spin precession rate of Saturn and the vertical secular eigenfrequency associated with Neptune, as well as the near match of Saturn's 27° obliquity to that of Cassini state 2 (see Section 2) of the secular spin–orbit resonance. Ward & Hamilton (2004) and Hamilton & Ward (2004) have proposed an elegant scenario, in which the secular spin–orbit resonance is responsible for tilting Saturn (but see Helled et al. 2009 for a potential problem with this scenario). Hamilton & Ward (2004) also note that, while they concentrate on changes in secular eigenfrequencies due to the changing gravitational potential of a gradually decaying Kuiper Belt, all that matters is that the frequencies change and not the underlying mechanism, and so a change in frequency due to migration should also work. At the same time, even though Jupiter's spin precession rate is close to the vertical secular eigenfrequency associated with Uranus (Ward & Canup 2006), its 3° obliquity means that its obliquity evolution was very different from Saturn's.

Whether this spin–orbit resonance mechanism can work in the various migration scenarios is therefore a natural question. Here we perform secular and particle N-body integrations and we initially consider smooth migration, Nice-like, and five planet models for the giant planets, and study the resulting obliquities. The obliquities of Jupiter and Saturn, in combination with the survival and orbital properties of four giant planets, provide strong constraints on the migration and encounter histories of the giant planets. In Section 2 we consider the conditions needed to tilt Saturn by spin–orbit resonance. In Section 3 we describe our numerical methods, and in Section 4 our initial configurations. We present our results in Section 5, with the secular and N-body simulations in Sections 5.1 and 5.2. The results are discussed in Section 6 and summarized in Section 7.

We construct a fixed reference frame with unit vectors $({\boldsymbol{i}},{\boldsymbol{j}},{\boldsymbol{k}})$ where ${\boldsymbol{k}}$ is parallel to the total orbital angular momentum of the system (and hence perpendicular to the invariable plane) and ${\boldsymbol{i}}$ and ${\boldsymbol{j}}$ are in the invariable plane. Consider a planet whose unit vector of the spin direction ${\boldsymbol{s}}$ is tilted at an angle from the unit orbit normal ${\boldsymbol{n}},$ i.e., $\mathrm{cos}\varepsilon ={\boldsymbol{s}}\cdot {\boldsymbol{n}}.$ The orbit has an inclination $\mathrm{cos}i={\boldsymbol{n}}\cdot {\boldsymbol{k}}.$ If there are no perturbations on the orbit and thus ${\boldsymbol{n}}$ is fixed in space, the torque from the Sun on the oblate figure of the planet (and on any satellites whose orbits are locked to the planet's equator plane) causes ${\boldsymbol{s}}$ to precess uniformly about ${\boldsymbol{n}}$ at the frequency $\alpha {(1-{e}^{2})}^{-3/2}\mathrm{cos}\varepsilon .$ Here e is the orbital eccentricity and α is the precession constant, given by

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{3G\;{M}_{\odot }}{2\nu {a}^{3}}\left(\displaystyle \frac{{J}_{2}+q^{\prime} }{{\mathcal{C}}+{\ell }}\right),\end{eqnarray} \tag{ 1 }$$

where G is the gravitational constant, ${M}_{\odot }$ is the solar mass, a is the orbital semimajor axis of the planet, ν is its rotation frequency, J₂ is the second zonal harmonic coefficient of its gravity field, ${\mathcal{C}}$ is its moment of inertia normalized to mR², m and R are the mass and radius of the planet, ℓ is the orbital angular momentum (normalized to ${{mR}}^{2}\nu$) of the satellites whose orbits are locked to the planet's equator, and ${q}^{\prime }/{J}_{2}$ is the ratio of the torque on the satellites to that directly exerted on the planet (see Ward & Hamilton 2004 for explicit definitions of ℓ and q').

Secular perturbations from other planets force ${\boldsymbol{n}}$ to precess about ${\boldsymbol{k}}.$ In the simple case where ${\boldsymbol{n}}$ precesses at a constant rate s, there are either two or four locations where ${\boldsymbol{s}}$ remains coplanar with ${\boldsymbol{n}}$ and ${\boldsymbol{k}},$ and ${\boldsymbol{s}}$ and ${\boldsymbol{n}}$ precess at the same rate s about ${\boldsymbol{k}}.$ These secular spin–orbit resonance configurations are called the Cassini states (Colombo 1966; Peale 1969, 1974).

Ward & Hamilton (2004) and Hamilton & Ward (2004) suggest that Saturn's obliquity of nearly 27° may be the result of Saturn being caught in Cassini state 2 with the vertical secular eigenfrequency s₈ associated with Neptune.³ They show that Saturn's obliquity ${\varepsilon }_{{\rm{S}}}$ can be increased from an initially small value to 27° by capture into Cassini state 2 if $| {s}_{8}/{\alpha }_{{\rm{S}}}| \gt 1$ initially and decreases to $\lt 1.$ They also show that the more complicated nodal regression of Saturn caused by perturbations from Jupiter and Uranus do not affect the capture into spin–orbit resonance and the associated growth in Saturn's obliquity.

One way to decrease $| {s}_{8}|$ (and hence $| {s}_{8}/{\alpha }_{{\rm{S}}}|$) is through the outward migration of Neptune. The timescale ${\tau }_{s}={\alpha }_{{\rm{S}}}/{\dot{s}}_{8}$ for the decrease in $| {s}_{8}|$ must be adequately long for Saturn's spin to stay in resonance. The critical timescale ${\tau }_{s,{\rm{min}}}$ can be estimated from requiring the time that it takes for the obliquity of Cassini state 2 to move by an angle comparable to the width of the separatrix to be longer than the libration period about Cassini state 2. With the obliquity of Cassini state 2 given by $\mathrm{cos}{\varepsilon }_{{\rm{S}}}\approx -({s}_{8}/{\alpha }_{{\rm{S}}})\mathrm{cos}i$ (when ${\varepsilon }_{{\rm{S}}}$ is large), the half width of the separatrix ${\rm{\Delta }}{\varepsilon }_{{\rm{S}}}=2{(\mathrm{tan}i/\mathrm{tan}{\varepsilon }_{{\rm{S}}})}^{1/2},$ and the libration frequency ${\omega }_{{\rm{lib}}}={(-{\alpha }_{{\rm{S}}}{s}_{8}\mathrm{sin}{\varepsilon }_{{\rm{S}}}\mathrm{sin}i)}^{1/2}$ for small-amplitude libration about Cassini state 2 (Henrard & Murigande 1987; Hamilton & Ward 2004; Ward & Hamilton 2004). Then $2{\rm{\Delta }}{\varepsilon }_{{\rm{S}}}/{\dot{\varepsilon }}_{{\rm{S}}}\approx 2\pi /{\omega }_{{\rm{lib}}}$ gives

$$\begin{eqnarray}&&{\tau }_{s,{\rm{min}}}\approx \displaystyle \frac{\pi }{2{\alpha }_{{\rm{S}}}\mathrm{sin}i\mathrm{sin}{\varepsilon }_{{\rm{S}}}}.\end{eqnarray} \tag{ 2 }$$

The condition ${\dot{s}}_{8}\approx {\omega }_{{\rm{lib}}}^{2}$ found by Hamilton & Ward (2004) gives a similar ${\tau }_{s,{\rm{min}}}$ (but without the coefficient $\pi /2$). Boué et al. (2009) have also shown using a different analytic argument that the critical timescale ${\tau }_{s,{\rm{min}}}\propto 1/({\alpha }_{{\rm{S}}}i).$

We have verified Equation (2) with numerical simulations using the symplectic integrator SyMBA, which was modified to include the evolution of the spin axis (Lee et al. 2007; see also Section 3). We integrated Saturn alone with ${\alpha }_{{\rm{S}}}=0\buildrel{\prime\prime}\over{.} 79\;{{\rm{yr}}}^{-1}$ (Tremaine 1991), $e={10}^{-3}$ and imposed orbital nodal precession rate s that changes linearly from $-5{\alpha }_{{\rm{S}}}$ to $-0.88{\alpha }_{{\rm{S}}}$ according to $s(t)=(-5+t/{\tau }_{s}){\alpha }_{{\rm{S}}}.$ The spin axis is initially parallel to ${\boldsymbol{k}}.$ Figure 1(a) shows the evolution of Saturn's obliquity in a set of simulations with $i=0\buildrel{\circ}\over{.} 1$ and different ${\tau }_{s}.$ Saturn is captured into Cassini state 2 and remains there until the end of the simulation (when ${\varepsilon }_{{\rm{S}}}$ is slightly larger than 27°) if ${\tau }_{s}\gt {\tau }_{s,{\rm{min}}}\approx 700\;{\rm{Myr}}.$ For ${\tau }_{s}\lt {\tau }_{s,{\rm{min}}},$ Saturn is initially captured into Cassini state 2, but it escapes from the resonance when the obliquity becomes too large; the final obliquity increases with increasing ${\tau }_{s}$ and can be more than 10° for ${\tau }_{s}\gtrsim \displaystyle \frac{1}{2}{\tau }_{s,{\rm{min}}}.$ From similar simulations with $i=0\buildrel{\circ}\over{.} 01-10^\circ ,$ we have determined ${\tau }_{s,{\rm{min}}}$ as a function of i, and the results are shown in Figure 1(b) for Saturn, as well as Jupiter with ${\alpha }_{{\rm{J}}}=2\buildrel{\prime\prime}\over{.} 67\;{{\rm{yr}}}^{-1}$. We confirm Equation (2) over a wide range in i and find that for Saturn ${\tau }_{s,{\rm{min}}}\sim 20\;{\rm{Myr}}$ for i of a few degrees. We note that $20\;{\rm{Myr}}$ is roughly the upper limit of the typical migration timescale of the giant planets caused by planetesimal scattering (Morbidelli et al. 2010). Also, in our earlier simulations of the classic Nice model (Lee et al. 2007), we found that the orbital inclination of Saturn can reach $\sim 5^\circ$ (and those of the ice giants 10°–15°) during the encounter phase in some simulations (see, e.g., Figure 2 of Lee et al. 2007).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Summarizing, it appears possible to tilt Saturn to its current obliquity if $| {s}_{8}/{\alpha }_{{\rm{S}}}|$ changes slowly enough. In the next sections we investigate whether or not this can be realized in the various migration scenarios.

To study the obliquity evolution in the various migration scenarios, we performed several types of numerical simulations of the Sun and the giant planets. One set of simulations consisted of integrating the secular evolution of the migrating giant planets without encounters and MMRs. The second set of simulations integrated the full N-body problem with planetesimals. These integrations used the Kepler-adapted symplectic N-body code SyMBA (Duncan et al. 1998), a descendant of the original techniques of Wisdom & Holman (1991) and Kinoshita et al. (1991), as modified by Lee et al. (2007) to incorporate spin evolution. There are some additional modifications for both the secular and N-body simulations that we now describe.

3.1. Secular Simulations

Since we expect the obliquity evolution to be driven by resonances between the nodal regression and spin precession frequencies, both rather small compared to the orbital frequencies, a useful approximation to the full N-body behavior for exploratory simulations can be obtained by modeling the purely secular evolution of the system. We have adapted the spin-SyMBA code of Lee et al. (2007) to remove mean-motion effects and close encounters when desired. The method is the following.

We consider a system made up of the Sun and some number of planets. The planets are assumed to undergo principle axis rotation, and the effects of the oblateness of the planets on the orbital evolution are neglected. The integration step has the form

$$\begin{eqnarray}&&{M}^{\;\tau /2}\;{R}^{\;\tau /2}\;{S}^{\;\tau /2}\;{D}^{\;\tau }\;{S}^{\;\tau /2}\;{R}^{\;\tau /2}\;{M}^{\;\tau /2}\end{eqnarray} \tag{ 3 }$$

where τ is the time step and each term is an operator. The operator D advances the planets along their osculating Kepler orbits, S advances the spin axes according to the mutual torques between planets and the torques due to the Sun, R handles the secular interactions between the planets, and M generates radial migration (whose recipe is discussed later in this subsection). Unlike with the full code described by Lee et al. (2007), the integration of the secular system is performed in heliocentric coordinates rather than canonical heliocentric, and so there is no linear drift term to correct for the barycentric momentum of the Sun (Duncan et al. 1998). Since we do not treat close encounters in the secular approximation, we do not need the recursive subdivision of the time step of SyMBA. The Kepler drift operator is straightforward, and the spin torque operator follows Lee et al. (2007), but the secular operator replacing the interaction term is new.

The secular interaction between the planets can be treated in several ways. Here we forgo the eigenvalue approach and instead directly integrate Lagrange's planetary equations with the first-order secular potential (see, e.g., Murray & Dermott 1999). Let a, e, i be the semimajor axis, eccentricity, and inclination of a planet, ϖ the longitude of perihelion, Ω the longitude of the ascending node, and n the mean motion. To apply the R operator, for each pair of planets (with indices $j,k;j\ne k$) we write the disturbing function on planet j due to k as

$$\begin{eqnarray}{R}_{{jk}} & = & {n}_{j}\;{a}_{j}^{2}\left[\displaystyle \frac{1}{2}\;{A}_{{jj}}\;{e}_{j}^{2}+{A}_{{jk}}\;{e}_{j}\;{e}_{k}\;\mathrm{cos}({\varpi }_{j}-{\varpi }_{k})\right.\\ & & \left.+\displaystyle \frac{1}{2}\;{B}_{{jj}}\;{i}_{j}^{2}+{B}_{{jk}}\;{i}_{j}\;{i}_{k}\;\mathrm{cos}({{\rm{\Omega }}}_{j}-{{\rm{\Omega }}}_{k})\right].\end{eqnarray} \tag{ 4 }$$

Here the matrix elements A and B are functions of the masses and semimajor axes and they contain information about the coupling between the planets (Brouwer & Woerkom 1950). We can then compute the time derivatives

$$\begin{eqnarray}&&{\dot{e}}_{{jk}}=-\displaystyle \frac{1}{{n}_{j}\;{a}_{j}^{2}\;{e}_{j}}\displaystyle \frac{\partial {R}_{{jk}}}{\partial {\varpi }_{j}},\;\quad {\dot{i}}_{{jk}}=-\displaystyle \frac{1}{{n}_{j}\;{a}_{j}^{2}\;{i}_{j}}\displaystyle \frac{\partial {R}_{{jk}}}{\partial {{\rm{\Omega }}}_{j}},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\dot{\varpi }}_{{jk}}=\displaystyle \frac{1}{{n}_{j}\;{a}_{j}^{2}\;{e}_{j}}\displaystyle \frac{\partial {R}_{{jk}}}{\partial {e}_{j}},\;\quad {\dot{{\rm{\Omega }}}}_{{jk}}=\displaystyle \frac{1}{{n}_{j}\;{a}_{j}^{2}\;{i}_{j}}\displaystyle \frac{\partial {R}_{{jk}}}{\partial {i}_{j}}.\end{eqnarray} \tag{ 6 }$$

The total rate of change in eccentricity of planet j is given by

$$\begin{eqnarray}&&{\dot{e}}_{j}=\displaystyle \sum _{k\ne j}{\dot{e}}_{{jk}}\end{eqnarray} \tag{ 7 }$$

and similarly for the changes in the other elements ${\dot{i}}_{j},{\dot{\varpi }}_{j},$ and ${\dot{{\rm{\Omega }}}}_{j}.$ These equations are integrated using a fourth-order Runge–Kutta scheme.

It may seem odd to evolve what is fundamentally a secular ring problem using point particles, as this requires repeatedly solving Kepler's equation during the Kepler drifts, which should be unnecessary in a secular model. However, in this way the well-tested spin vector infrastructure of the Lee et al. (2007) modifications to the SyMBA code could be reused. We have confirmed by comparing with a different secular code (one which solves the corresponding eigenproblem instead) and with the full N-body integration that this approach recovers the expected evolution.

Note that the secular step itself is rather time-consuming compared to a simple force evaluation due to the need to compute multiple trigonometric functions and the Laplace coefficients. The separation between the orbital timescale and the secular, spin, and migration timescales provides several opportunities for optimization. Advancing the secular term R every (orbital timescale) step is unnecessary, and applying it every 10 steps serves to make the runtime comparable to the N-body run, without changing the resulting evolution. Another improvement can be obtained by dividing each of the secular, spin, and migration timescales by some constant factor (making sure to keep them several orders of magnitude longer than the orbital timescale) during the simulation and then scaling the resulting output. This is possible because the behavior of the spin–orbit resonance depends not on the absolute values of these timescales but on their ratios. In our regime, using a factor of 10–20 produces negligible errors relative to the original simulation and provides almost all of the possible speed benefits. We have verified that our results are insensitive to these approximations, which are used only in the secular integrations.

In the standard picture, the final configuration of the giant planets should evolve mostly due to mutual interactions and interactions with an external remnant planetesimal disk (such as the Kuiper Belt). Accordingly, with the possible exception of objects which get scattered out into the Kuiper Belt and Scattered Disk (Duncan & Levison 1997), imposing a prescription for migration and random velocity damping during the secular simulations should be a reasonable approximation to the effect of the planetesimal disk. Let t be the time since the start of the simulation, and τ the migration timescale. We follow Lee et al. (2007), and take as migration expressions

$$\begin{eqnarray}&&{\dot{a}}_{k}=\displaystyle \frac{{\rm{\Delta }}{a}_{k}}{\tau }\mathrm{exp}(-t/\tau ),\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&\quad {\dot{a}}_{k}=\displaystyle \frac{{\rm{\Delta }}{a}_{k}}{\tau }\left(\displaystyle \frac{t}{\tau }\right)\mathrm{exp}[-{t}^{2}/(2{\tau }^{2})],\end{eqnarray} \tag{ 9 }$$

where ${\rm{\Delta }}{a}_{k}$ is a parameter measuring the amount of distance the planet should travel due to the migration to end up on its current orbit. In the secular problem, ${\rm{\Delta }}{a}_{k}$ is simply the difference between the current value and the starting value. Equation (8) results in faster migration at earlier times than Equation (9), but the two resulting trajectories meet again at $t=2\tau$ when the distance travelled is $\sim 0.86{\rm{\Delta }}{a}_{k},$ after which Equation (9) has faster migration. Equation (8) has a continually decreasing $\dot{a},$ whereas Equation (9) has $\dot{a}$ increasing until a maximum at τ (corresponding to a maximum speed of $\sim 0.6{\rm{\Delta }}{a}_{k}/\tau$) beyond which it decays to zero. The M operator of Equation (3) simply applies the changes in a to the planets. As in Lee et al. (2007), we evolve all planets under Equation (8) in our integrations of the smooth migration model, and we use Equation (8) for Jupiter and Saturn and Equation (9) for Uranus and Neptune in our integrations of the Nice model. Note that we do not apply any eccentricity or inclination damping to the secular simulations: as there is no stirring mechanism, e and i would damp too quickly to serve as a useful model for the spin–orbit resonance crossing. We evolve every simulation to $10\tau .$

This concludes our description of the secular simulations. We now move on to discuss the full N-body simulations.

3.2. N-body Simulations

We restrict ourselves to four classes of initial conditions, namely a version of the smooth migration scenario of Hahn & Malhotra (1999) with four planets, a variant of the resonant Nice model (Morbidelli et al. 2007) with four planets, and two five planet models from Nesvorný & Morbidelli (2012). We integrate the Hahn & Malhotra (1999) initial conditions with the secular equations only, the resonant Nice with the secular equations and full N-body model, and the five planet cases with full N-body only.

4.1. Smooth Migration

4.2. Resonant Nice Configuration

In this case, the four giant planets are all initially in resonance. Saturn and Jupiter are in the 3:2, as are Saturn and the inner ice giant. There are two distinct stable configurations for Uranus and Neptune found by Morbidelli et al. (2007), namely the 4:3 and the 5:4; we study only the 4:3, which has the outer ice giant at 11.6 AU. To generate the initial conditions, we used an approach similar to that of Morbidelli et al. (2007), but instead of a hydrodynamic code, we employed the N-body code SyMBA, with forced migration and eccentricity damping to model the effects of the proto-planetary gas disk. We began with Jupiter and Saturn and placed them just outside the 3:2 MMR. Jupiter was forced to migrate outwards and Saturn inwards, and the migration rates were chosen so that they would be nearly stationary after capture into the 3:2 resonance. The eccentricity damping rate is specified by Lee et al. (2007)

$$\begin{eqnarray}&&{\dot{e}}_{k}/{e}_{k}=-{K}_{e}| {\dot{a}}_{k}/{a}_{k}| .\end{eqnarray} \tag{ 10 }$$

We adopted K_e = 100 for Jupiter and adjusted K_e for Saturn to match the equilibrium eccentricities reported by Morbidelli et al. (2007). After the Jupiter–Saturn pair had reached equilibrium in the 3:2 resonance, we added an ice giant just outside the 3:2 resonance with Saturn and forced it to migrate inwards until Saturn and the ice giant reached equilibrium in the 3:2 resonance. This process was repeated for the 4:3 resonance between the ice giants. For the ice giants, which should be undergoing type I migration, we assumed that the values of K_e are identical and adjusted K_e until we matched the eccentricities reported by Morbidelli et al. (2007) at each stage of this process. To model the dispersal of the gas disk, we continued the integration with an exponential decay in the migration and eccentricity damping rates.

4.3. Five-planet Configurations

In this section we present the results of our secular and N-body simulations. We discuss the smooth migration case first to explain the underlying mechanism of how to tilt the gas giants, followed by the resonant Nice and five planet simulations.

5.1. Secular Simulations

To understand how the secular spin–orbit resonance scenario works in the full planetary migration context, but without the complicating factors present in noisy N-body simulations, we apply the secular spin approach of Section 3.1 as an intermediate scheme between the one-planet simulations in Section 2 and the full N-body runs to two of the histories (smooth migration and resonant Nice) described in Section 4.

5.1.1. Smooth Migration

The simplest model to consider is a purely secular model of the smooth migration scenario, as it behaves as the analytic theory would predict. Figure 2 presents an overview of a successful run that reproduces the current obliquities of the gas giants. In this example $\tau =32\;{\rm{Myr}}$ and initial $i=2^\circ$ for all planets. The upper panels describe the evolution of the orbits, with the migration in semimajor axes imposed by Equation (8). The inclinations oscillate around the forced values due to secular interactions. The lower panels describe the evolution of the system spins. The obliquities of all four planets remain low until approximately $20\;{\rm{Myr}},$ during which time the evolution of is driven entirely by changes in the inclination, as can be seen by comparing the bottom left and top right panels. However, after about $20\;{\rm{Myr}}$ the obliquity of Saturn rises to a mean of $\sim 33^\circ$ with a superimposed oscillation of about 4^◦. That this obliquity rise is due to capture into a secular spin–orbit resonance is clear from the bottom right panel, which compares the spin precession frequencies α of the planets with the three non-degenerate vertical secular eigenfrequencies $| {s}_{6}| ,| {s}_{7}| ,$ and $| {s}_{8}|$ of the system. At the beginning of the simulation, no eigenfrequency is particularly close to any spin precession frequency of a planet, with the slowest eigenfrequency $| {s}_{8}|$ (which is that dominated by Neptune's contribution) situated between the spin precession frequencies of Jupiter and Saturn. However, as the semimajor axes change due to migration, so do both the eigenfrequencies and the spin precession frequencies. After $20\;{\rm{Myr}},$ when Saturn is observed to begin tilting over, we see that $| {s}_{8}|$ matches the spin precession frequency of Saturn, and it is this secular spin–orbit resonance which could explain the current tilt.

Not every set of parameters results in such a Saturnian obliquity: as discussed in Section 2, both τ and i play a role such that $\tau \times i$ must exceed a certain value. Figure 3 shows the resulting Jovian and Saturnian obliquities from simulations with $\tau =1-100\;{\rm{Myr}}$ and initial $i=0\buildrel{\circ}\over{.} 01-2^\circ .$ In general, as the timescale increases so does the resulting Saturnian obliquity, up to a maximum of $\sim 30^\circ .$ This is very close to the observed 27^◦, thus providing a satisfying explanation for Saturn's obliquity. Moreover, in this migration timescale regime, at least an inclination of $1^\circ -2^\circ$ is required, which also gives a Jovian obliquity compatible with the current value. Thus smooth migration can work at large $\tau \times i$ (≳30 Myr deg): the minimum $\tau \times i$ producing a simulation with final Saturnian obliquity $\gt 15^\circ$ is 32 Myr deg, although some simulations with $\tau \times i$ above this limit have final Saturnian obliquity $\lt 15^\circ .$ This joint constraint will prove to be useful.

Zoom In Zoom Out Reset image size

Although the smooth migration model is able to reproduce Saturn's obliquity if $\tau \times i$ is sufficiently large, we want to emphasize that it has great difficulties in explaining the current secular architecture of the giant planets (Morbidelli et al. 2009) and the terrestrial planets (Brasser et al. 2009), as well as the orbital properties of the main belt asteroids (Morbidelli et al. 2010), and should not be considered as a proxy for the past evolution of the solar system.

5.1.2. Resonant Nice

In principle the situation with the resonant Nice model should be very similar. The outward migration of Saturn and the ice giants will lower the spin precession frequency of Saturn slightly and the eigenfrequencies considerably, so the same resonance crossing should occur, and it does. The difficulty becomes apparent when we examine an equivalent to the smooth migration run shown in Figure 2.

Figure 4 shows the evolution of a resonant Nice model with $\tau =10\;{\rm{Myr}}$ and initial ${i}_{{\rm{N}}}=4^\circ$ (the initial inclinations of the other planets are zero). This relatively high inclination of Neptune was chosen to mimic the typical inclination value during the encounter phase from the N-body simulations. The Nice initial configuration is significantly more compact, which has the effect of increasing the eigenfrequencies relative to the spin precession frequencies. In particular, the same frequency that successfully tilted Saturn in the smooth migration scenario, $| {s}_{8}| ,$ is now larger than Jupiter's spin precession frequency and must cross it to reach its final value near Saturn's. The situation for Jupiter is exactly comparable to the case of tilting Saturn, and the only question is whether the product of the migration timescale and the inclination is large enough at the time of the resonance crossing near $10\;{\rm{Myr}}$ to produce a significant Jovian obliquity. (Recall that Jupiter can be tilted by a faster decrease in $| {s}_{8}|$ than Saturn; see Figure 1.) Unfortunately, the answer is yes. In this run, at the end we have ${\varepsilon }_{{\rm{J}}}\simeq 6^\circ ,$ which is much greater than the observed 3^◦, even though Jupiter does not stay in resonance until the end. Moreover, Saturn's obliquity of $\sim 8^\circ$ is not substantially larger. To increase ${\varepsilon }_{{\rm{S}}}$ we require an increase in either ${i}_{{\rm{N}}}$ or τ or both, but since the same resonance is at work and the crossing times are only separated by $\sim \tau ,$ improving the Saturn value will make Jupiter's obliquity unrealistically large. Figure 5 summarizes the outcomes from a series of resonant Nice runs with τ varying from 1 to $10\;{\rm{Myr}}$ and the initial i_N from 1° to 10°. As feared, ${\varepsilon }_{{\rm{J}}}\sim {\varepsilon }_{{\rm{S}}}:$ the mean and standard deviation of ${\varepsilon }_{{\rm{J}}}/{\varepsilon }_{{\rm{S}}}=0.67\pm 0.12,$ whereas the observed ratio is 0.11. The only runs which show ${\varepsilon }_{{\rm{S}}}\geqslant 15^\circ$ have $10^\circ \leqslant {\varepsilon }_{{\rm{J}}}\leqslant 20^\circ ,$ and have long migration times (8 or $10\;{\rm{Myr}}$).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

However, the secular models leave out many relevant dynamical effects. For example, we apply semimajor axis migration via our prescriptions, and the real system will undergo a more sudden transition during resonance crossing or breaking and will be far more stochastic. Perhaps Jupiter's spin–orbit resonance crossing happens much more quickly than Saturn's, preventing too large a tilt. In the secular simulations we also let the inclinations undergo the usual evolution. Some experiments (not presented here) introducing inclination damping proved unfruitful as the Jupiter crossing occurs first, and lowering the inclination enough to avoid tilting Jupiter results in a far too small Saturnian obliquity. In the real system, one could imagine that the inclinations are low during the Jupiter spin–orbit resonance crossing and higher due to encounters when Saturn crosses, recovering the smooth migration model behavior. To determine whether these possibilities are realistic scenarios we must move beyond the secular models to more realistic N-body models, and to those we now turn.

5.2. N-body Simulations

We saw from the secular simulations that in the resonant Nice model the obliquities of Jupiter and Saturn are similar because of $| {s}_{8}|$ crossing the spin precession frequencies of both planets. This did not occur in the smooth migration setup because of the wider spacing of the planets.⁵ We now present the results of our full N-body simulations. In what follows we need to quantify what we consider to be a successful simulation that adequately reproduces the secular properties and spacing of the outer solar system, in addition to the obliquities of Jupiter and Saturn. To that end we identify several conditions that the system must adhere to. These are similar to the requirements of Nesvorný & Morbidelli (2012) but with some modifications. First, we need to have four giant planets (Jupiter, Saturn, and two ice giants) at the end of the simulation—criterion A of Nesvorný & Morbidelli (2012). Second, the planets need to be within 10% of their current orbits, i.e., $0.9\leqslant {a}_{f}/{a}_{c}\leqslant 1.1$ where a_f is the final semimajor axis of each planet at the end of the simulation, and a_c is its current semimajor axis—similar to criterion B of Nesvorný & Morbidelli (2012) which imposes a margin of 20%. Third, we impose a restriction on the difference in the longitudes of perihelion, ${\rm{\Delta }}{\varpi }_{{\rm{JS}}}={\varpi }_{{\rm{J}}}-{\varpi }_{{\rm{S}}}.$ The current secular evolution of the eccentricities of Jupiter and Saturn is such that ${\rm{\Delta }}{\varpi }_{{\rm{JS}}}$ circulates through 360° and that the eccentricity of Jupiter (Saturn) reaches a maximum when ${\rm{\Delta }}{\varpi }_{{\rm{JS}}}=180^\circ$ (0°). Morbidelli et al. (2009) have shown that, at the very least, encounters between an ice giant and Saturn are needed to produce this secular behavior. Thus, we require that ${\rm{\Delta }}{\varpi }_{{\rm{JS}}}$ circulates. We used Fourier analysis over the last 5 Myr of the simulation to establish the behavior of ${\rm{\Delta }}{\varpi }_{{\rm{JS}}}.$ This approach differs slightly from Nesvorný & Morbidelli (2012) who keep track of the amplitude of the eccentricity eigenmode of Jupiter (their criterion C). For the final obliquities, we require that ${\varepsilon }_{{\rm{J}}}\lt 5^\circ$ and ${\varepsilon }_{{\rm{S}}}\gt 15^\circ .$

During the course of this investigation, we became aware of Jupiter often being too close to the Sun at the end of our N-body simulations. For example, the mean and standard deviation of the final ${a}_{{\rm{J}}}=4.98\pm 0.21\;{\rm{AU}}$ for the simulations of the resonant Nice model with four planets at the end. This affects not only whether the planets are within 10% of their current orbits, but also potentially the obliquities of Jupiter and Saturn, because the secular spin–orbit resonance problem is not scale free. The precession constants scale with the semimajor axis of Jupiter as $\alpha \propto {a}_{{\rm{J}}}^{-3}$ while the secular eigenfrequencies scale as $s\propto {a}_{{\rm{J}}}^{-3/2},$ if we fix the semimajor axis ratios. So the final $| s/\alpha |$ would be smaller by $\approx (3/2){\rm{\Delta }}{a}_{{\rm{J}}}/{a}_{{\rm{J}}}$ if the final ${a}_{{\rm{J}}}$ is closer than the actual value $5.20\;{\rm{AU}}$ by ${\rm{\Delta }}{a}_{{\rm{J}}}.$ Since the initial ${a}_{{\rm{J}}}$ is in some sense arbitrary, we rescale each simulation so that the final ${a}_{{\rm{J}}}=5.20\;\mathrm{AU}.$ If the final $| {s}_{8}/{\alpha }_{{\rm{S}}}| \lt 1$ (i.e., there is resonance crossing) before rescaling and $| {s}_{8}/{\alpha }_{{\rm{S}}}| \gt 1$ (i.e., there is no resonance crossing) after rescaling, we use the obliquity of Saturn before the crossing in the original (unscaled) simulation as its final value. We apply the same procedure to $| {s}_{7}/{\alpha }_{{\rm{J}}}|$ for Jupiter.

5.2.1. Final Orbits and Obliquities

First we analyze the outcome from the N-body simulations of the resonant Nice model, with the planets residing initially in the 3:2, 3:2, 4:3 resonant chain and a planetesimal disk of $50\;{M}_{\oplus }.$ In Figure 6 we plot the final ${\varepsilon }_{{\rm{J}}}$ versus ${\varepsilon }_{{\rm{S}}}.$ The values of both ${\varepsilon }_{{\rm{J}}}$ and ${\varepsilon }_{{\rm{S}}}$ are averages taken over the last 50 Myr in order to wash out any large-amplitude, long-period oscillations, except, as noted earlier, when there is secular spin–orbit resonance crossing before rescaling but no crossing after rescaling ($| {s}_{7}/{\alpha }_{{\rm{J}}}|$ for Jupiter and $| {s}_{8}/{\alpha }_{{\rm{S}}}|$ for Saturn), in which case the obliquity is taken before the crossing. The left panel is a log–log version of the right panel. The small black dots depict the final obliquities for simulations that have four planets (but the planets may or may not be near their correct positions). The larger blue dots refer to simulations where the planets are also near the correct positions (but ${\rm{\Delta }}{\varpi }_{{\rm{JS}}}$ may librate or circulate). Finally, the red dots denote simulations that match all of our criteria in orbital properties for success. The dotted line denotes ${\varepsilon }_{{\rm{J}}}={\varepsilon }_{{\rm{S}}}.$

Zoom In Zoom Out Reset image size

There are several visible trends in Figure 6. First, the obliquities of both Saturn and Jupiter take on values ranging from under 1° to over 30°. They generally do not follow the line ${\varepsilon }_{{\rm{J}}}={\varepsilon }_{{\rm{S}}}.$ Second, there are roughly equal numbers of cases with ${\varepsilon }_{{\rm{J}}}\gt {\varepsilon }_{{\rm{S}}}$ and ${\varepsilon }_{{\rm{J}}}\lt {\varepsilon }_{{\rm{S}}}.$ We find that ${\varepsilon }_{{\rm{J}}}\gt {\varepsilon }_{{\rm{S}}}$ 50% of the time for all simulations that produce four planets (irrespective of their final orbits), and ${\varepsilon }_{{\rm{J}}}\gt {\varepsilon }_{{\rm{S}}}$ 51% of the time for the cases when the planets are close to their current orbits. However, the mean and standard deviation of ${\varepsilon }_{{\rm{J}}}/{\varepsilon }_{{\rm{S}}}=1.86\pm 2.75$ and 2.24 ± 3.40 for each sample, indicating that the distribution of ${\varepsilon }_{{\rm{J}}}$/${\varepsilon }_{{\rm{S}}}$ has a longer tail at large values. Compared to the secular simulations of the same model in Section 5.1.2, the mean ${\varepsilon }_{{\rm{J}}}$/${\varepsilon }_{{\rm{S}}}$ is even further from the observed value of 0.11, but there is a much larger spread in ${\varepsilon }_{{\rm{J}}}$/${\varepsilon }_{{\rm{S}}}$.

Unfortunately only very few simulations matched our criteria for being considered successful. A total of 84 runs out of 1280 (6.6%) had the planets in the right positions. Of these only 13 have ${\rm{\Delta }}{\varpi }_{{\rm{JS}}}$ circulate. Of these 13 only one has ${\varepsilon }_{{\rm{S}}}\gt 15^\circ$ and ${\varepsilon }_{{\rm{J}}}\lt 5^\circ$ (see also Table 1). However, it is not clear how essential of a constraint the circulation of ${\rm{\Delta }}{\varpi }_{{\rm{JS}}}$ is. Nesvorný & Morbidelli (2012) state that requiring the circulation of ${\rm{\Delta }}{\varpi }_{{\rm{JS}}}$ is the real bottleneck in matching the current orbital architecture. Removing this restriction yields two additional good cases matching the planetary architecture and the obliquities. In summary, it appears that the four-planet model has some difficulty reproducing the current solar system. This fact was already established by Nesvorný & Morbidelli (2012), who reported similar statistics for the orbital constraints, but the probability is even lower (∼0.08%–0.23%) when we consider the obliquities.

Table 1.  Statistics of N-body Simulation Results

Model	${M}_{{\rm{disk}}}$	Initial ${a}_{{\rm{N}}}$	f4	fa	${f}_{\varpi }$	${f}_{\varepsilon }$
3:2 3:2 4:3	50 M⊕	11.6 AU	35%	6.6%	1.0%	0.08%–0.23%
3:2 3:2 4:3 4:3	35 M⊕	14.2 AU	23%	6.7%	1.9%	0.16%–0.47%
3:2 3:2 2:1 3:2	20 M⊕	22.2 AU	27%	10%	3.3%	<0.08%

Download table as:  ASCIITypeset image

Next we report on the outcome of the N-body simulations with five initial planets. Figure 7 is similar to Figure 6 but now we focus on the compact five planet case with the initial configuration 3:2, 3:2, 4:3, 4:3 and a planetesimal disk of $35\;{M}_{\oplus }.$ The outcome is similar to, but slightly better than, the resonant Nice case with four planets: we have six cases where the giant planets are in the right positions and have their current obliquities, but in only two of these ${\rm{\Delta }}{\varpi }_{{\rm{JS}}}$ circulates. There are fewer cases where Jupiter is tilted more than Saturn: ${\varepsilon }_{{\rm{J}}}\gt {\varepsilon }_{{\rm{S}}}$ 37% of the time for all simulations that produce four planets, and ${\varepsilon }_{{\rm{J}}}\gt {\varepsilon }_{{\rm{S}}}$ 33% of the time for the cases when the planets are close to their current orbits. But $\langle {\varepsilon }_{{\rm{J}}}/{\varepsilon }_{{\rm{S}}}\rangle =1.44\pm 2.12$ for all final cases with four planets, and 1.45 ± 2.02 for when the final planets are close to their current orbits, again indicating that there is a longer tail at large values of ${\varepsilon }_{{\rm{J}}}$/${\varepsilon }_{{\rm{S}}}$.

Zoom In Zoom Out Reset image size

Finally, we turn to the loose five planet system, whose initial configuration is 3:2 3:2 2:1 3:2 and the outermost ice giant resides at $22.2\;{\rm{AU}}.$ the mass of the planetesimal disk is $20\;{M}_{\oplus }.$ we plot the final obliquities of Jupiter and Saturn in Figure 8. the main difference with the same figures for the resonant nice case (Figure 6) and the compact five planet configuration (Figure 7) is an absence of high ${\varepsilon }_{{\rm{S}}}$ at low ${\varepsilon }_{{\rm{J}}}.$ only one simulation with four planets too far from their correct positions at the end has final ${\varepsilon }_{{\rm{S}}}\gt 15^\circ$ and ${\varepsilon }_{{\rm{J}}}\lt 5^\circ .$ the fraction of simulations with ${\varepsilon }_{{\rm{J}}}\gt {\varepsilon }_{{\rm{S}}}$ is 45% for all runs that have four planets at the end and 61% for runs where the planets are near their current orbits. we also find that $\langle {\varepsilon }_{{\rm{J}}}/{\varepsilon }_{{\rm{S}}}\rangle =3.98\pm 3.45$ and 4.39 ± 7.24 for both samples, respectively, so that this wider initial configuration of the planets leaves Jupiter's obliquity higher compared to Saturn's than the compact one.

Zoom In Zoom Out Reset image size

5.2.2. Evolution Examples

We now illustrate the mechanisms behind the tilting of Jupiter and Saturn through some figures. As was discussed earlier, tilting Saturn requires that ${\alpha }_{{\rm{S}}}$ resonates with s₈, or passes through this resonance. For Jupiter, on the other hand, there are two ways to tilt it: passage through resonance with s₈ or s₇. Figure 9 shows the evolution of the system that tilts both Jupiter and Saturn to their current obliquities in a simulation of the compact five planet model. The top-left panel shows the perihelion distance, semimajor axis and aphelion distance of all planets as a function of time. To guide the eye each planet is identified by a different color. The top-right panel depicts the period ratio between Jupiter and Saturn. The bottom-right panel depicts ${\varepsilon }_{{\rm{J}}}$ (black) and ${\varepsilon }_{{\rm{S}}}$ (blue) and the bottom-left panel shows the evolution of ${\alpha }_{{\rm{J}}}$ (magenta), ${\alpha }_{{\rm{S}}}$ (black), $| {s}_{7}|$ (red) and $| {s}_{8}|$ (blue) as a function of time. When starting with three ice giants there is an extra vertical eigenfrequency that disappears when the third ice giant is ejected. The frequencies decrease with increasing semimajor axis and this property allowed us to keep track of which frequency corresponded to which planet and ultimately identify s₇ and s₈ belonging to the two surviving ice giants.

Zoom In Zoom Out Reset image size

After approximately 3 Myr of evolution, the system becomes unstable. The innermost ice giant becomes Neptune, and the outermost ice giant evolves little and ends up being Uranus. The period ratio between the gas giants increases rapidly when Saturn scatters the lost ice giant and then evolves more slowly. The two remaining ice giants eventually settle into stable orbits at t ≈ 7 Myr and migrate toward their current orbits without undergoing encounters. The frequency $| {s}_{8}|$ crosses ${\alpha }_{{\rm{J}}}$ shortly after the system goes unstable but the crossing occurs during the jump and is too quick to allow capture in resonance. It then steadily approaches ${\alpha }_{{\rm{S}}}$ and after about 100 Myr it is close enough for Saturn to be caught in resonance. Saturn's obliquity increases and it stays in resonance until the end of the simulation. The resonant state is marked by the long-period oscillations in the obliquity of Saturn. At the same time $| {s}_{7}|$ approaches ${\alpha }_{{\rm{J}}}$ but does not come close enough for capture in resonance because Uranus stays fairly close to the Sun.

From Section 2 we know that a slow migration of Neptune is needed to tilt Saturn. From Fourier analysis we obtain for the above simulation ${I}_{68}=0\buildrel{\circ}\over{.} 065$ and ${\alpha }_{{\rm{S}}}=0\buildrel{\prime\prime}\over{.} 693\;{{\rm{yr}}}^{-1},$ where I₆₈ is the forced inclination in Saturn associated with s₈. Thus ${\tau }_{s,{\rm{min}}}\sim {10}^{9}\;\mathrm{year}$, according to Figure 1(b) and Equation (2). In the simulation Neptune migrates outwards by 0.4 AU between t = 80 Myr and t = 350 Myr, and then essentially stops. This suggests that $\dot{a}\sim {\rm{\Delta }}a/{\rm{\Delta }}t\sim 2\times {10}^{-9}\;{\rm{AU}}\;{{\rm{yr}}}^{-1}.$ The timescale for migration is then ${t}_{{\rm{mig}}}=a/\dot{a}\gg {10}^{9}\;\mathrm{year}$ and thus Saturn gets captured in resonance.

Unfortunately we also have many cases where Jupiter gets tilted significantly more than Saturn. In Figure 10 we plot the evolution of a simulation of the compact five planet model where Jupiter gets tilted but Saturn does not. The panels are the same as for Figure 9. The outermost ice giant is ejected when the system becomes unstable at t ≈ 9 Myr. From the bottom-left panel we see that ${\alpha }_{{\rm{J}}}$ gets caught in resonance with s₇. Thus it appears that cases of high obliquity of Jupiter are often caused by capture or passage through resonance with s₇ rather than s₈ because the latter occurs early in the violent phase when the planets migrate fast. Applying Figure 1(b) to Jupiter, we have ${I}_{57}=0\buildrel{\circ}\over{.} 072$ and ${\alpha }_{{\rm{J}}}=3\buildrel{\prime\prime}\over{.} 24\;{{\rm{yr}}}^{-1}$ and thus ${\tau }_{s,{\rm{min}}}\sim 200\;\mathrm{Myr}$. Uranus migrates outward by 0.3 AU from t = 50 Myr until t = 200 Myr so that ${t}_{{\rm{mig}}}\gg 200\;\mathrm{Myr}$. Thus Jupiter gets caught in the resonance.

Zoom In Zoom Out Reset image size

5.2.3. Spin–Orbit Resonance Crossing

Figure 11 shows the cumulative distribution of Jupiter's final obliquity, compared with the cumulative distribution of Jupiter's obliquity before any $| {s}_{7}| /{\alpha }_{{\rm{J}}}$ crossing (if one occurs), for the simulations with four planets in the right positions at the end (blue dots in Figures 6–8). Except for one simulation in the compact five-planet model, all simulations have ${\varepsilon }_{{\rm{J}}}\lt 10^\circ$ before any $| {s}_{7}| /{\alpha }_{{\rm{J}}}$ crossing occurs, but at least 18% of the simulations have final ${\varepsilon }_{{\rm{J}}}\gt 10^\circ .$ This confirms that the cases of high Jupiter obliquity are often caused by capture or passage through resonance with s₇ instead of s₈.

Zoom In Zoom Out Reset image size

To understand the differences in the distributions of the final obliquities of Jupiter and Saturn between the three models, we compare in Figures 12–14 the final values of $| {s}_{7}| ,| {s}_{8}| ,{\alpha }_{{\rm{J}}},$ and ${\alpha }_{{\rm{S}}}$ for the N-body simulations with four planets in the right positions at the end. In each figure, the filled circles and error bars show the mean values and three standard deviations ($3\sigma$) of the final values of $| {s}_{7}|$ versus ${a}_{{\rm{U}}}$ and $| {s}_{8}|$ versus a_N, with s₇ and s₈ from the Laplace–Lagrange secular theory. The red line shows the final ${\alpha }_{{\rm{J}}},$ which is fixed after rescaling the final ${a}_{{\rm{J}}}$ to $5.20\;{\rm{AU}},$ and the cyan band shows the $\pm 3\sigma$ range of the final values of ${\alpha }_{{\rm{S}}}.$ The vertical magenta dashed lines indicate the ±10% range of the current ${a}_{{\rm{U}}}$ and ${a}_{{\rm{N}}}.$

Zoom In Zoom Out Reset image size

From Figure 12 for the resonant Nice model, we see that $| {s}_{7}|$ is about $0.7\sigma$ above ${\alpha }_{{\rm{J}}},$ which means that $| {s}_{7}|$ crosses ${\alpha }_{{\rm{J}}}$ in about 24% of the simulations with four planets in the right positions, if we assume a normal distribution. This agrees with the $\sim 24\%$ of simulations with final ${\varepsilon }_{{\rm{J}}}\gt 10^\circ$ (see Figure 11), showing again that the cases where Jupiter has a high obliquity are most likely caused by a resonance crossing with s₇ instead of s₈. For the compact five planet model (Figure 13), $| {s}_{7}|$ is further ($\sim 1\sigma$) above ${\alpha }_{{\rm{J}}},$ and there are fewer ($\sim 20\%$) simulations with ${\varepsilon }_{{\rm{J}}}\gt 10^\circ .$ On the other hand, although more than half of the loose five planet simulations have final $| {s}_{7}| \lt {\alpha }_{{\rm{J}}}$ (see Figure 14), only $\sim 18\%$ of the simulations have ${\varepsilon }_{{\rm{J}}}\gt 10^\circ .$ This suggests that, in this model, the $| {s}_{7}| /{\alpha }_{{\rm{J}}}$ resonance crossing often occurs when Uranus is migrating too quickly, which is consistent with the fact that Uranus is on average farther from the Sun than in the other two models.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Unfortunately, this problem of resonance crossing occurring too quickly also prevents the tilting of Saturn in many cases. For the resonant Nice model (Figure 12), more than half of the simulations have $| {s}_{8}| \lt {\alpha }_{{\rm{S}}},$ but only $\sim 8\%$ of the simulations have ${\varepsilon }_{{\rm{S}}}\gt 10^\circ ,$ indicating that the $| {s}_{8}| /{\alpha }_{{\rm{S}}}$ resonance crossing often occurs when Neptune is migrating too quickly. Although the compact five planet model has a similar mean value of ${a}_{{\rm{N}}},$ it has a narrower range in ${a}_{{\rm{N}}}$ and more overlap between $| {s}_{8}|$ and ${\alpha }_{{\rm{S}}}$ (see Figure 13), which results in a higher fraction ($\sim 19\%$) of simulations with ${\varepsilon }_{{\rm{S}}}\gt 10^\circ .$ Finally, in the loose five planet model (Figure 14), Neptune is on average too far from the Sun, and almost all of the simulations have $| {s}_{8}| \lt {\alpha }_{{\rm{S}}},$ but only $\sim 3\%$ of the simulations have ${\varepsilon }_{{\rm{S}}}\gt 10^\circ .$

5.2.4. Summary

We summarize the outcome of our N-body simulations in Table 1. The column f₄ is the probability of the simulation having four giant planets after 500 Myr. The column f_a is the probability that the system has four planets with each planet within 10% of its current orbit. The next column, f_ϖ, is the probability of having four planets close to their current orbits and ${\rm{\Delta }}{\varpi }_{{\rm{JS}}}$ circulating. The last column, f, is the probability that the system also has ${\varepsilon }_{{\rm{J}}}\lt 5^\circ$ and ${\varepsilon }_{{\rm{S}}}\gt 15^\circ .$ For the resonant Nice and compact five planet models, we list two values for f, with the lower value requiring the circulation of ${\rm{\Delta }}{\varpi }_{{\rm{JS}}}$ and the higher value removing this restriction on the state of ${\rm{\Delta }}{\varpi }_{{\rm{JS}}}.$ For the loose five planet model, none of the simulations with four planets close to their current orbits has ${\varepsilon }_{{\rm{J}}}\lt 5^\circ$ and ${\varepsilon }_{{\rm{S}}}\gt 15^\circ ,$ and there is only an upper limit on f.

Nesvorný & Morbidelli (2012) have examined the statistics of orbital properties for a large number of models, each with only 30–100 N-body simulations. In this paper, we have examined three representative models with 1280 simulations each. Our results on the orbital statistics are much more accurate, and they are consistent with those found by Nesvorný & Morbidelli (2012). In particular, the loose five planet model is the most successful in reproducing the positions of the giant planets and the circulation of ${\rm{\Delta }}{\varpi }_{{\rm{JS}}}.$ However, it is clear from Table 1 that the loose five planet model fails to reproduce the obliquities of Jupiter and Saturn (with the probability ${f}_{\varepsilon }\lt 0.08\%$) with the setup we have chosen. Thus the obliquities of Jupiter and Saturn provide a strong constraint that can alter the relative merits of different models and initial conditions.

In summary, it appears that both the resonant Nice and compact five planet models are able to reproduce simultaneously the current orbital architecture of the giant planets and the obliquities of the gas giants, with the compact five planet model faring slightly better. However, the probability is less than $0.5\%,$ even if we lift the restriction on the state of ${\rm{\Delta }}{\varpi }_{{\rm{JS}}}.$ The loose five planet model fails with a probability of less than $0.08\%.$

Even though the results presented in the previous section are extensive, they are by no means complete. In this section we discuss a few outstanding issues.

First, we did not impose any condition on the migration speed of the giant planets. Nesvorný & Morbidelli (2012) imposed having the Saturn–Jupiter period ratio evolve from below 2.1 to beyond 2.3 in shorter than 1 Myr (their criterion D). We verified that all nine cases in our N-body simulations that reproduce the current obliquities of the gas giant have some sort of jump in the period ratio, although only one case from the compact five planet model had a jump in the right range, and is depicted in Figure 9. This would suggest that the compact five planet case is better overall at matching all constraints than the resonant Nice model.

Second, we only tested a single four planet and two five planet configurations because we opted to bracket the two end points of the initial spacing of the giant planets to determine whether the final outcomes favor one case over the other. The answer is clearly yes. We have had to run many more simulations per set of initial conditions than Nesvorný & Morbidelli (2012) because we are trying to satisfy an additional constraint, namely reproducing the obliquities of the gas giants. Given that we found at most a few cases out of 1280 that match all constraints from one set of initial conditions, it becomes evident why we ran as many simulations as we did for each set of initial conditions, and why it is infeasible to test a larger sample of initial conditions.

Third, we generally did not change the mass of the planetesimal disk when running the simulations for one configuration. One could argue that if we increased the disk mass from 35 M_⊕ to 50 M_⊕ in the compact five planet case we may get a higher probability matching all constraints. However, increasing the disk mass from 35 M_⊕ to 50 M_⊕ would most likely damp the eccentricities of Jupiter and Saturn too much and make them inconsistent with their current values—see Nesvorný & Morbidelli (2012) for a more detailed discussion. We ran the same compact five planet model with a 50 M_⊕ disk. Even though the probabilities f₄ to ${f}_{\varepsilon }$ show no improvement and are nearly identical to those of the compact five planet model with the 35 M_⊕ disk, we found that the eccentricities of Jupiter and Saturn were inconsistent with their current values ($\gt 2\sigma$ result). Thus, increasing the disk mass does not help our case and we do not report the results here. In addition, our current disk mass of 35 M_⊕ typically increases the period ratio of Jupiter and Saturn by 0.2 after the jump, which Nesvorný & Morbidelli (2012) state may indicate the disk is too massive.

Another issue that we have not considered is the uncertainties in the spin precession constant ${\alpha }_{{\rm{J}}}$ and ${\alpha }_{{\rm{S}}}.$ Vokrouhlický & Nesvorný (2015) have published an analysis based on one-planet simulations similar to those in Section 2 (i.e., integration of spin evolution for a planet whose orbital nodal precession rate $| s|$ is decreased on timescale ${\tau }_{s}$). Their analysis differs from ours in assuming that the final values of s₇ and s₈ are fixed at the observed solar system ones, but they considered a range of values for both ${\alpha }_{{\rm{J}}}$ and ${\alpha }_{{\rm{S}}}.$ For the $| {s}_{8}| /{\alpha }_{{\rm{J}}}$ resonance crossing, they find constraints on I₅₈ and ${\tau }_{s}$ for the resonance crossing to be too fast for capture and Jupiter's obliquity to be $\leqslant 3^\circ$ after this crossing (which agrees with the open symbols in Figure 11). Similarly, for the final approach of s₇ to its current value, they find constraints on ${\alpha }_{{\rm{S}}}$ and ${\tau }_{s}$ that yield Jupiter's final obliquity to be $\leqslant 3^\circ .$ For the final approach of s₈ to its current value, they considered different values of ${\alpha }_{{\rm{S}}}$ and ${\tau }_{s}$ and find that capture into the $| {s}_{8}| /{\alpha }_{{\rm{S}}}$ resonance requires large ${\tau }_{s}$ and low ${\alpha }_{{\rm{S}}}$ (including the value of ${\alpha }_{{\rm{S}}}$ adopted in the present paper). However, Vokrouhlický & Nesvorný (2015) did not assess the probability that such conditions on the migration rates are satisfied in the various models of Nesvorný & Morbidelli (2012), nor did they assess the probability of these models giving the right values of s₇ and s₈, which we have done here for three specific configurations. Thus, while they suggest that the mechanism behind tilting Saturn is the resonance crossing between s₈ and ${\alpha }_{S},$ we show here how well this mechanism performs with full N-body simulations in several specific settings.

A last point that requires elaboration is whether the obliquities of the gas giants can be used to constrain the migration of the giant planets at all. At issue is the following. The current solar system has ${\varepsilon }_{{\rm{J}}}\sim 3^\circ$ and ${\varepsilon }_{{\rm{S}}}\sim 27^\circ$ and for both planets ${\alpha }_{{\rm{J}}}\mathrm{cos}{\varepsilon }_{{\rm{J}}}\sim | {s}_{7}|$ and ${\alpha }_{{\rm{S}}}\mathrm{cos}{\varepsilon }_{{\rm{S}}}\sim | {s}_{8}| .$ With a 10% margin in the final semimajor axes of the planet one can compute the range of expected final obliquities of the gas giants assuming that during a resonance crossing there is perfect capture i.e., $\mathrm{cos}{\varepsilon }_{{\rm{J}}}^{\prime }=| {s}_{7}^{\prime }/{\alpha }_{{\rm{J}}}|$ if $| {s}_{7}^{\prime }/{\alpha }_{{\rm{J}}}| \leqslant 1,$ and $\mathrm{cos}{\varepsilon }_{{\rm{J}}}^{\prime }=1$ otherwise. Here primed quantities are the values at the end of the simulation. Doing a series expansion of the Laplace–Lagrange theory around the current semimajor axes of the ice giants we find that ${\varepsilon }_{{\rm{J}}}^{\prime }\in (0,44^\circ )$ and ${\varepsilon }_{{\rm{S}}}^{\prime }\in (0,53^\circ ).$ With such a large range of final obliquities over a small range of final semimajor axes of the giant planets, one might argue that the final obliquities of the gas giants cannot be used to constrain the migration.

However, this argument requires one assumption that does not stand up to scrutiny. Indeed, if there was perfect capture then in a plot of $| s/\alpha |$ versus the data points should be clustered near $\varepsilon =\mathrm{arccos}(| s/\alpha | )$ when $| s/\alpha | \leqslant 1$ and near $\varepsilon =0$ otherwise. Figures 15–17 depict the final obliquities of the gas giants as a function of $| s/\alpha | .$ It is clear that this trend is only partly visible. We see some clustering of high ${\varepsilon }_{{\rm{J}}}$ near $| {s}_{7}/{\alpha }_{{\rm{J}}}| \sim 1$—and in some cases near $| {s}_{7}/{\alpha }_{{\rm{J}}}| \gt 1$ because of the inherent uncertainty in the Laplace–Lagrange model. For Jupiter the trend is more obvious than it is for Saturn, and the trend for Saturn is almost absent in the loose five planet model. This indicates that when a crossing does occur, capture often does not. In other words, the assumption that during a crossing there is perfect capture is clearly untrue.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

So what are we to make of this? The answer is that for Saturn's obliquity to be where it is, two things need to be satisfied simultaneously: (a) Neptune's migration needs to be slow enough so that $\tau \times i\gtrsim 30$ Myr deg, and (b) the spacing of the planets needs to be such that $| {s}_{8}/{\alpha }_{{\rm{S}}}| \approx 1.$ We showed in Figures 12–14 that the final semimajor axes of Uranus and Neptune (and indirectly, also the one of Saturn) are not the same for each set of simulations. Consequently, the final values of s₇ and s₈ are not the same either, despite the requirement of the final orbits being within 10% of the current ones. As an example, for the resonant Nice model, the average and standard deviation of the final semimajor axes of Saturn, Uranus and Neptune in the cases where the planets are near their current orbits, are $\langle {a}_{{\rm{S}}}\rangle =8.98\pm 0.352\;{\rm{AU}},\;$ $\langle {a}_{{\rm{U}}}\rangle =18.69\pm 0.99\;{\rm{AU}}$ and $\langle {a}_{{\rm{N}}}\rangle =30.14\pm 1.70\;{\rm{AU}}.$ In contrast, for the compact five-planet model the values are $\langle {a}_{{\rm{S}}}\rangle =9.33\pm 0.451\;{\rm{AU}},$ $\langle {a}_{{\rm{U}}}\rangle =18.69\pm 0.956\;{\rm{AU}}$ and $\langle {a}_{{\rm{N}}}\rangle =29.74\pm 1.08\;{\rm{AU}}.$ Even though these are all within 10% of the current values, the current spacing of the giant planets is better matched for the compact five planet case than for the resonant Nice case. It is the difference in the final spacing, coupled with the tail of Neptune's migration, that determines the obliquities of the gas giants and whether one model matches the observed constraints better than another. In light of this, we conclude that the obliquities of the gas giants can indeed be used to constrain the evolution and initial conditions of the giant planets.

We have performed secular and N-body simulations tracking planetary obliquity to study constraints on migration histories of the outer solar system from both spin and orbit, covering the smooth migration scenario, the resonant Nice model and two five planet scenarios. The obliquities of Jupiter and Saturn provide a strong constraint on the models. The secular results leave us in the following situation: in the very simplest case, that of smooth migration, the Ward & Hamilton (2004) secular spin–orbit resonance mechanism for tilting Saturn's spin axis appears to work nicely if the product of the migration timescale and the orbital inclinations is sufficiently large ($\tau \times i\gtrsim 30$ Myr deg). On the other hand, the resonant Nice model, which is preferable on many grounds (such as explaining the orbital properties of the giant planets, terrestrial planets, and main belt asteroids), is more problematic, because the secular eigenfrequency s₈ is initially above the spin precession frequency of Jupiter, ${\alpha }_{{\rm{J}}}$. The fundamental problem is that there are incompatible constraints. We need to migrate slowly ($\tau \gtrsim 8\;{\rm{Myr}}$) at typical inclinations to consistently tilt Saturn's spin axis by $\gt 15^\circ ,$ but for those same inclinations we need to migrate quickly (say, $\tau \lesssim 2\;{\rm{Myr}}$) to avoid tilting Jupiter by more than the observed 3°. We find that on average ${\varepsilon }_{{\rm{J}}}\sim {\varepsilon }_{{\rm{S}}}$ in our secular simulations of the resonant Nice model.

At the same time, the N-body simulations appear to tell a different story. The resonant Nice model is able to reproduce the current orbital architecture of the giant planets and the obliquities of Jupiter and Saturn, but only with a small probability. The compact five planet model fares slightly better, but the loose five planet model has great difficulty reproducing the obliquities, even though it is the most successful in reproducing the orbital properties. However, it is possible that a slight change in the initial conditions or the outer edge of the planetesimal disk could improve the outcome.

Ultimately the following needs to happen: (1) There needs to be fast migration during the encounter phase to avoid tilting Jupiter through resonance passage with s₈. (2) Then a late, slow migration of Neptune to its current location could complete the task by tilting Saturn through the resonance with s₈. (3) At the same time Uranus must stay close enough to the Sun to avoid tilting Jupiter through the resonance with s₇. Condition (1) is almost always satisfied in the N-body simulations. But the chances of satisfying conditions (2) and (3) are limited, as the $| {s}_{8}| /{\alpha }_{{\rm{S}}}$ resonance crossing often occurs when Neptune is migrating too fast, especially in the loose five-planet model, while Jupiter sometimes crosses the $| {s}_{7}| /{\alpha }_{{\rm{J}}}$ resonance. Thus the main obstacle encountered to reproduce the obliquities of both Jupiter and Saturn is the final orbital spacing of the giant planets, coupled with the tail of Neptune's migration. Both the resonant Nice and compact five planet models appear to be able to overcome the obstacle, but with low probability.

The authors are grateful for the support of Hong Kong RGC Grant HKU 7030/11P. We thank an anonymous reviewer for his/her valuable comments. We thank D. Nesvorný for providing one of his five planet configurations and fruitful discussions, A. Morbidelli for comments on an earlier version, S. J. Peale for useful discussions during the early stages of this work, and D. S. McNeil for earlier work on this project. We are grateful for being able to use the TIARA grid computing cluster at the Institute for Astronomy and Astrophysics, Academia Sinica, Taiwan.