Tilting Uranus: Collisions versus Spin–Orbit Resonance

Gravitational torques from the Sun on an oblate planet cause the planet's spin axis to precess backwards, or regress, about the normal to its orbital plane (Colombo 1966). Similarly, gravitational perturbations cause a planet's inclined orbit to regress around the Sun. A match between these two precession frequencies results in a secular spin–orbit resonance. In this case, the spin axis remains fixed relative to the planet's orbital pole, and the two vectors precess about the normal to the invariable plane. The longitudes of the two axial vectors, ϕ_α and ϕ_g , are measured from a reference polar direction to projections onto the invariable plane, and the resonance argument is given as (Hamilton & Ward 2004)

$$\begin{eqnarray}&&{\rm{\Psi }}={\phi }_{\alpha }-{\phi }_{g}.\end{eqnarray} \tag{ 1 }$$

The precession rate of Uranus's spin axis can be derived from first principles by considering the torques of the Sun and the Uranian moons on the planet's equatorial bulge. Following Colombo (1966), if $\hat{\sigma }$ is a unit vector that points in the direction of the total angular momentum of the Uranian system, then

$$\begin{eqnarray}&&\displaystyle \frac{d\hat{\sigma }}{{dt}}=\alpha (\hat{\sigma }\times \hat{n})(\hat{\sigma }\cdot \hat{n}),\end{eqnarray} \tag{ 2 }$$

where $\hat{n}$ is a unit vector pointing in the direction of Uranus's orbital angular momentum, α is the spin precession rate near zero degree tilts, and t is time. Uranus's axial precession period is therefore

$$\begin{eqnarray}&&{T}_{\alpha }=\displaystyle \frac{2\pi }{\alpha \cos \epsilon },\end{eqnarray} \tag{ 3 }$$

where is the obliquity and $\cos \epsilon =\hat{\sigma }\cdot \hat{n}$. The precession frequency near zero obliquity, α, incorporates the torques from the Sun and the planet's moons on the central body (Tremaine 1991):

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{3{n}^{2}}{2}\displaystyle \frac{{J}_{2}(1-\tfrac{3}{2}{\sin }^{2}{\theta }_{p})+q}{K\omega \cos {\theta }_{p}+l}.\end{eqnarray} \tag{ 4 }$$

Here, $n={({{GM}}_{\odot }/{r}_{p}^{3})}^{1/2}$ is the orbital angular speed of the planet, G is the gravitational constant, M_⊙ is the Sun's mass, r_p is the Sun–planet distance, ω is the planet's spin angular speed, J₂ is its quadrupole gravitational moment, and K is its moment of inertia normalized by ${M}_{p}{R}_{p}^{2}$. For Uranus today, M_p = 14.5 M_⊕, R_p = 2.56 × 10⁹ cm, K = 0.225, and J₂ = 0.00334343. ¹ The parameter $q\equiv \tfrac{1}{2}{\sum }_{i}{({M}_{i}/{M}_{p})({a}_{i}/{R}_{p})}^{2}(1-\tfrac{3}{2}{\sin }^{2}{\theta }_{i})$ is the effective quadrupole coefficient of the satellite system, and $l\equiv {R}_{p}^{-2}{\sum }_{i}{({M}_{i}/{M}_{p})({{GM}}_{p}{a}_{i})}^{\tfrac{1}{2}}\cos {\theta }_{i}$ is the angular momentum of the satellite system divided by ${M}_{p}{R}_{p}^{2}$. The masses and semimajor axes of the satellites are M_i and a_i , $\cos {\theta }_{p}=\hat{s}\cdot \hat{\sigma }$, and $\cos {\theta }_{i}=\hat{{l}_{i}}\cdot \hat{\sigma }$, where $\hat{s}$ is the direction of the spin angular momenta of the central body, and $\hat{{l}_{i}}$ is the normal to the satellite's orbit (Tremaine 1991). Note that M_i ≪ M_p , where M_p is the mass of the planet, and because the satellite orbits are nearly equatorial, we can take θ_p = θ_i = 0.

Torques from the main Uranian satellites on the planet contribute significantly to its precessional motion, while those from other planets and satellites can be ignored. We therefore limit ourselves to Uranus's major moons—Oberon, Titania, Umbriel, Ariel, and Miranda. We find q = 0.01558, which is about 4.7 times larger than Uranus's J₂, and l = 2.41 × 10⁻⁷, which is smaller than K ω by about a factor of 100. So from Equation (4), the effective quadrupole coefficient of the satellite system plays a much more significant role in the planet's precession period than the angular momentum of the satellite system. At its current obliquity, = 98°, Uranus's precession period is about 210 million years (or α = 0.0062 arcsec yr⁻¹), and reducing Uranus's obliquity to 0° results in a precession period 7.2 times faster: 29 million years (or α = 0.045 arcsec yr⁻¹). This pole precession rate is much longer than any of the giant planets' fundamental frequencies (Murray & Dermott 1999), but it can be sped up to ≈2 Myr by placing Uranus at around 7 au. This is just fast enough for Uranus to resonate with a similar planet—Neptune—located beyond Saturn.

Placing Uranus's orbit between those of Jupiter and Saturn is not entirely ad hoc. Thommes et al. (1999, 2002, 2003) argue that at least the ice giants' cores might have formed between Jupiter and Saturn (4–10 au), as the timescales there for the accretion of planetesimals through an oligarchic growth model, when the large bodies in the planetary disk dominate the accretion of surrounding planetesimals, are more favorable than farther away. The Nice model (Gomes et al. 2005; Morbidelli et al. 2005; Tsiganis et al. 2005) places Uranus closer to the Sun but beyond Saturn for similar reasons; however, having the ice giants form between Jupiter and Saturn is not inconsistent with the Nice model. If Uranus and Neptune were indeed formed between Jupiter and Saturn and later ejected sequentially, then a secular spin–orbit resonance between Uranus and Neptune is possible. Note that such close encounters would not yield any significant obliquity excitations because the perturbing torque is too weak as it depends on the planet's gravitational quadrupole moment (Lee et al. 2007). A related possibility that is also sufficient for our purposes is if the planets were formed from pebble accretion, as the pebble isolation mass can be similar everywhere in the outer solar system (Lambrechts et al. 2014), allowing Neptune to be initially formed beyond Saturn and Uranus between Jupiter and Saturn. In the following, we assume that Uranus is fully formed with its satellites located near their current configurations to derive the spin axis precession rate. We also include only the known giant planets because adding a third or fourth ice giant, as suggested by the Nice model to better reproduce the solar system (Nesvorný 2011; Batygin et al. 2012; Nesvorný & Morbidelli 2012), would increase the planet's orbital precession rates and make it more difficult for Uranus to obtain a spin–orbit resonance. Only if we succeed to tilt Uranus reasonably under these conditions would we consider introducing more giant planets to the model.

Calculating Uranus's obliquity evolution requires tracking the planets' orbits while also appropriately tuning Neptune's nodal precession rate. We use the HNBody Symplectic Integration package (Rauch & Hamilton 2002) to track the motion of bodies orbiting a central massive object using symplectic integration techniques based on two-body Keplerian motion, and we move Neptune radially with an artificial drag force oriented along the velocity vector using the package HNDrag. These packages do not follow spins, so we have written an integrator that uses a fifth-order Runge–Kutta algorithm (Press et al. 1992) and reads in HNBody data to calculate Uranus's axial orientation due to torques applied from the Sun (Equation (2)). For every time step, the integrator requires the distance between the Sun and Uranus. Because HNBody outputs the positions and velocities at a given time frequency different from the adaptive step that our precession integrator uses, calculating the precessional motion requires interpolation. To minimize interpolation errors, we use a torque averaged over an orbital period that is proportional to $\langle {r}_{p}^{-3}\rangle ={a}_{p}^{-3}{(1-{e}_{p}^{2})}^{-\tfrac{3}{2}}$, where a_p is the planet's semimajor axis and e_p is its eccentricity. This is an excellent approximation because Uranus's orbital period is 10⁵–10⁶ times shorter than its precession period. We tested the code for a two-body system consisting of just the Sun and Uranus, and recovered the analytic result for the precession of the spin axis (Figure 1).

Zoom In Zoom Out Reset image size

For our simulations, we place Jupiter and Saturn near their current locations (5 au and 9 au, respectively), Uranus at 7 au, and Neptune well beyond Saturn at 17 au. Leaving Uranus in between the two gas giants for more than a few million years is unstable (Lecar & Franklin 1973; Franklin et al. 1989; Gladman & Duncan 1990; Holman & Wisdom 1993), but eccentricity dampening from remnant planetesimals can delay the instability. Scattering between Uranus and the planetesimals provides a dissipative force that temporarily prevents Uranus from being ejected, and we mimic this effect by applying an artificial force to damp Uranus's eccentricity. We apply the force in the orbital plane and perpendicular to the orbital velocity to damp the eccentricity while preventing changes to the semimajor axis (Danby 1992). With Uranus's orbit relatively stable, we then seek a secular resonance between its spin and Neptune's orbit.

Capturing into a spin–orbit resonance also requires the two angular momentum vectors, the planet's spin axis and an orbital pole, and the normal to the invariable plane be coplanar. Equilibria about which the resonance angle librates are called "Cassini states" (Colombo 1966; Peale 1969; Ward 1975; Ward & Hamilton 2004), and there are multiple vector orientations that can yield a spin–orbit resonance. In our case, the resonance angle, Ψ, librates about Cassini state 2 because Uranus's spin axis and Neptune's orbital pole precess on opposite sides of the normal to the invariable plane.

As Neptune migrates outwards away from the Sun, its nodal precession frequency slows until a resonance is reached with Uranus's spin precession rate. If the consequence of the resonance is that Uranus's obliquity increases (Ward 1974), then its spin precession period increases as well (Equation (3)), and the resonance can persist. The time evolution of the resonance angle and obliquity are given by Hamilton & Ward (2004):

$$\begin{eqnarray}&&\dot{{\rm{\Psi }}}=-\alpha \cos \epsilon -g\cos I,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\dot{\epsilon }=g\sin I\sin {\rm{\Psi }},\end{eqnarray} \tag{ 6 }$$

where g is the negative nodal precession rate, and I is the amplitude of the inclination induced by Neptune's perturbation on Uranus's orbit. If Neptune migrates outward slowly enough, then $\dot{{\rm{\Psi }}}$ is small and the two planets can remain in resonance nearly indefinitely.

Figure 2 shows Uranus undergoing capture into a spin–orbit resonance when Neptune crosses ∼24 au en route to its current location at 30 au. Here we have set Neptune's migration rate to ∼0.05 au Myr⁻¹, which is within the adiabatic limit—the fastest possible rate to generate a capture with _i ≈ 0°. The adiabatic limit occurs when Neptune's migration takes it across the resonance width in about a libration time, which is just 2π/w_lib with ${w}_{\mathrm{lib}}=\sqrt{-\alpha g\sin \epsilon \sin I}$ (Hamilton & Ward 2004). Just as slow changes to the support of a swinging pendulum do not alter the pendulum's motion, gradual changes to Neptune's orbit do not change the behavior of the libration. However, if Neptune's migration speed exceeds the adiabatic limit, then resonance cannot be established. The top panel of Figure 2 shows Uranus tilting to 60° in 150 Myr when Neptune reaches its current location, and all the way to 90° in 600 Myr if we allow Neptune to continue outwards. Planets migrate by scattering planetesimals, which can decrease inclinations; accordingly, we optimistically assumed an initial value for Neptune's inclination at twice its current value. Because we have increased Neptune's inclination and moved Neptune out as fast as possible and yet still allowed capture, 150 million years represents a rough lower limit to the time needed to tilt Uranus substantially.

Zoom In Zoom Out Reset image size

The bottom panel of Figures 2 and 3 both show the evolution of the resonance angle, and the angle oscillates with a libration period of about 30 Myr about the equilibrium point. The libration period increases as increases in accordance with the predictions of w_lib. The noticeable offset of the equilibrium below Ψ = 0° in Figures 2 and 3 is due to the rapid migration of Neptune (Hamilton & Ward 2004):

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{\mathrm{eq}}=\displaystyle \frac{\dot{\alpha }\cos \epsilon \,+\,\dot{g}\cos I}{\alpha g\sin \epsilon \sin I}.\end{eqnarray} \tag{ 7 }$$

Recall that g, the nodal precession frequency, is negative, α is positive, and as Neptune migrates away from the Sun, $\dot{g}$ is positive. Because α is constant, $\dot{\alpha }=0$, and so Ψ_eq is slightly negative in agreement with Figure 2. We conclude that although a spin–orbit resonance with Neptune can tilt Uranus over, the model requires that Uranus be pinned between Jupiter and Saturn for a few hundred million years, yet close encounters with either gas giant makes it unstable to leave Uranus there for more than a million years (Lecar & Franklin 1973; Franklin et al. 1989; Gladman & Duncan 1990; Holman & Wisdom 1993). Is there any room for improvement?

Zoom In Zoom Out Reset image size

Solar system evolution models calculate planetary migration timescales on the order of 10⁶–10⁷ yr (Hahn & Malhotra 1999; Thommes et al. 1999, 2002, 2003; Gomes et al. 2005; Morbidelli et al. 2005; Hahn & Malhotra 2005; Tsiganis et al. 2005; de Sousa et al. 2020). This is incompatible with this resonance capture scenario, which requires at least 10⁸ yr. Speeding up the tilting timescale significantly would require a stronger resonance. The strength of this resonance is proportional to the migrating planet's inclination and it sets the maximum speed at which a capture can occur (Hamilton 1994). Although Neptune's initial orbital inclination angle is unknown, a dramatic reduction in the tilting timescale is implausible.

Another possibility is that the gas giants were once closer to the Sun where tidal forces are stronger. Some evidence for this comes from the fact that the giant planets probably formed closer to the snow line (Ciesla & Cuzzi 2006), where volatiles were cold enough to condense into solid particles. Shrinking the planets' semimajor axes by a factor of 10% decreases the resonance location by about 3 au and reduces the obliquity evolution timescale by about 15%. Although this is an improvement, a timescale on the order of 10⁸ yr seems to be the fundamental limit on the speed at which a significant obliquity can be reached (Rogoszinski & Hamilton 2016; Quillen et al. 2018).

Less critical than the timescale problem but still important is the inability of the obliquity to exceed 90° (Figure 2). The reason for this follows from Equation (3), which shows that Uranus's precession period approaches infinity as approaches 90°. Neptune's migration speed then is faster than the libration timescale, and the resonance ceases. This effect is more apparent in Figure 3, in which the libration period also increases with obliquity. The resonance breaks when the resonance angle stops librating about an equilibrium point and instead circulates a full 2π radians. Quillen et al. (2018) show that a related resonance that occurs when the planets are also close to a mean-motion resonance could tilt the planet past 90°, but this, like the resonance considered here, is probably too weak. Keeping Uranus between Jupiter and Saturn for 10⁸ yr is as implausible as the planet having once had a massive distant moon (Boué & Laskar 2010).

Recall that the leading hypothesis for Uranus's tilt is a single Earth-mass impactor striking the planet's polar region (e.g., Slattery et al. 1992; Ida et al. 2020), but that Morbidelli et al. (2012) argue for two or more collisions if the satellites were formed primordially from a circumplanetary disk (Szulágyi et al. 2018). In this section, we consider each of these scenarios and derive the resulting probability distributions for such impacts. To do this, we designed a collisional code that builds up a planet by summing the angular momenta of impactors to determine the planet's final obliquity and spin rate under various circumstances, and we typically run this for half a million randomized instances. Our assumptions are that the impactors originate within the protoplanetary disk, they approach a random location on the planet on trajectories that parallel its orbital plane, and all the mass is absorbed upon impact. Because nearly every object in the solar system orbits in roughly the same direction, the impactors' relative speed would be at most several tens of percent of Uranus's orbital speed (6.8 km s⁻¹). Because we expect most impactors to follow orbits with lower eccentricities, we sample relative velocities between 0 and 0.4 times Uranus's circular speed.

Considering that the impactor's relative speed is small compared to the planet's escape speed (21.4 km s⁻¹), gravitational focusing is also important. For cases where gravitational focusing is strong, the impact cross section is large and the impactor is focused to a hyperbolic trajectory aimed more closely toward the planet's center. Because head-on collisions do not impart any angular momentum, we expect the planet's spin state to be more difficult to change when focusing is included. The impact parameter for this effect is given by b with

$$\begin{eqnarray}&&{b}^{2}={R}_{p}^{2}(1+{\left({V}_{\mathrm{esc}}/{V}_{\mathrm{rel}}\right)}^{2}).\end{eqnarray} \tag{ 8 }$$

Also, because we do not know how the density profile changes between impacts, we maintain the dimensionless moment of inertia at $K\equiv \tfrac{I}{{M}_{p}{R}_{p}^{2}}=0.225$ but vary the planet's radius as the cube root of the total mass. Although these assumptions are mildly inconsistent, we find that even large impacts incident on a mostly formed Uranus yield just small changes in radius and that the final spin rates changes by only about 10% for other mass–radius relations. Finally, Podolak & Helled (2012) suggest a maximum impact boundary of around 0.95 R_p as beyond this, the impactor simply grazes the planet's atmosphere and departs almost unaffected. For simplicity, and in the spirit of approximation, we ignore this subtlety.

In Figures 11(a) and (b), we assume that the planet's initial spin rate was low to highlight the angular momentum imparted by impacts. Because ${V}_{\mathrm{esc}}^{2}=2{{GM}}_{p}/{R}_{p}$, the impact cross section b² ∝ R_p for V_rel ≪ V_esc (Equation (8)). The corresponding probability density distribution of impact locations is $d(\pi {b}^{2})/{{dR}}_{p}$, which is constant; therefore, the spin distribution induced from a single collision is flat (Figure 11(a)). However, if the impactor's relative speed is instead much greater than the planet's escape speed, then gravitational focusing is weak, and the impactors will be traveling on nearly straight lines. In this case, a single collision produces a spin distribution that increases linearly, as there is an equal chance of striking anywhere on the planet's surface. But because gravitational focusing only varies the radial concentration of impacts on a planet's surface, the obliquity distribution for a single impact onto an initially nonspinning planet with or without gravitational focusing is uniform. A Uranian core formed from the accretion of many small objects, by contrast, would likely have a very low spin rate (Lissauer & Kary 1991; Dones & Tremaine 1993a, 1993b; Agnor et al. 1999) because each successive strike likely cancels out at least some of the angular momentum imparted from the previous impact (Figures 11(c) and (d)). The planet would also have a narrower range of likely obliquities because the phase space available for low tilts is small.

Zoom In Zoom Out Reset image size

The calculation for the planet's final spin state for many impacts behaves similarly to a random walk, so from the central limit theorem, each directional component of the imparted angular momentum can be described by a normal distribution. The theoretical curve of Figure 11(c) is given by the probability distribution f_L (l), which describes the probability that L, the magnitude of the planet's spin angular momentum $L=\sqrt{{L}_{X}^{2}+{L}_{Y}^{2}+{L}_{Z}^{2}}$, takes the value l:

$$\begin{eqnarray}&&{f}_{L}(l)=\displaystyle \frac{2{l}^{2}{e}^{-{l}^{2}/2{\sigma }^{2}}}{\sqrt{2\pi }\,{\sigma }^{2}{\sigma }_{z}}\,{\rm{\Phi }}(0.5;1.5;-\beta {l}^{2})\end{eqnarray} \tag{ 9 }$$

(Dones & Tremaine 1993a, Equation (109)). Here, σ is the standard deviation for the components of the planet's spin angular momentum that lie in the orbital plane, σ_z is the standard deviation for the component perpendicular to the orbital plane, and $\beta =({\sigma }^{2}-{\sigma }_{z}^{2})/2{\sigma }^{2}{\sigma }_{z}^{2}$. The angular momentum imparted is always perpendicular to the plane of the impactor's trajectory. After multiple impacts, standard deviations are related by ${\sigma }_{z}\approx \sqrt{2}\sigma$, so β < 0. Finally, Φ(0.5; 1.5; β l²) is the confluent hypergeometric function of the first kind. The corresponding obliquity probability distribution is

$$\begin{eqnarray}&&{f}_{\epsilon }(\varepsilon )=\left|\displaystyle \frac{1}{4\sqrt{2}\,{\sigma }^{2}{\sigma }_{z}}\displaystyle \frac{\tan (\varepsilon )}{{\cos }^{2}(\varepsilon )}{\left(\displaystyle \frac{{\tan }^{2}(\varepsilon )}{2{\sigma }^{2}}+\displaystyle \frac{1}{2{\sigma }_{z}^{2}}\right)}^{-3/2}\right|\end{eqnarray} \tag{ 10 }$$

(Dones & Tremaine 1993a, Equation (111)); we provide derivations of these two equations in the Appendix. Notice how well these calculations agree with the numerical result for many impacts (Figures 11(c) and (d)). Consequentially, keeping the total mass imparted constant and increasing the number of impactors in Figure 11(c) from 100 to 1000 would shift the peak to slower spin rates by a factor of $\sqrt{10}$. Because Uranus's spin period is quite fast, its spin state could not have simply been a by-product of myriad small collisions.

Accordingly, we will now consider the intermediary cases with only a few impactors incident on a nonspinning planet. Figure 12 shows the product of two equal-sized hits, and the resulting distributions already resemble the limit of multiple collisions. If the masses of the two impactors differ significantly, however, the corresponding spin and obliquity distributions are more similar to those of the single impact case (Figure 13). Therefore, while the planet's obliquity distribution may be more or less flat, its spin rate strongly depends on both the number of strikes and the total mass in impactors.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1 shows a range of possible collisions onto a nonspinning planet. Here we show that the smallest amount of mass necessary to push Uranus toward its observed spin state is about 0.4 M_⊕, regardless of the number of impacts. The odds of this happening decrease for each additional collision because each impact needs to hit at exactly the right location. We also provide statistics for impactors much greater than an Earth-mass in the last section of Table 1. Impactors this massive would likely violate our no mass-loss assumption, yet the odds of generating Uranus's current spin state are still low. A more detailed analysis of these impacts is beyond the scope of this paper; however, see Kegerreis et al. (2018, 2019) for a smooth particle hydrodynamics analysis on the effects that impacts have on Uranus's rotation rate and internal structure.

We also explored cases with multiple unequal-sized impactors and discovered that the order of the impacts does not matter, as expected, and that the odds are improved for similar-sized impactors. An example of this can be seen in Figures 12(a) and 13(a), where for the same total mass, the spin distribution for two equally sized impactors is concentrated near Uranus's current spin state, whereas the distribution is flatter for two unequal-sized impacts. We conclude that a small number of equal impacts totaling to about 1 M_⊕ is the most likely explanation for Uranus's spin state if the planet was initially nonspinning.

Gas accretion almost certainly provides a significant source of angular momentum, so much so that we might expect the giant planets to be spinning at near break-up velocities if they accreted gas from an inviscid thin circumplanetary disk (Bodenheimer & Pollack 1986; Lissauer et al. 2009; Ward & Canup 2010). Instead, we observe the gas giants to be spinning several times slower, so there must have been some process for removing excess angular momentum. This mechanism may be a combination of multiple effects: magnetic braking caused by the coupling between a magnetized planet and an ionized disk (Lovelace et al. 2011; Batygin 2018), vertical gas flow into the planet's polar regions and additional midplane outflows from a thick circumplanetary disk (Tanigawa et al. 2012; Szulágyi et al. 2014), and magnetically driven outflows (Quillen & Trilling 1998; Fendt 2003; Lubow & Martin 2012; Gressel et al. 2013). Because both Uranus and Neptune spin at about the same rates and have likely harbored circumplanetary disks of their own (Szulágyi et al. 2018), we suspect that gas accretion is responsible, though pebble accretion may also contribute a significant amount of prograde spin (Visser et al. 2020). As such, the planet's initial obliquities should be near 0° as the angular momentum imparted by gas accretion is normal to the planet's orbital plane.

First, we explore cases where the planet initially spins slowly (T_i ≪ 17.2 hr). In Figure 14, we have Uranus's initial spin period four times slower than its current value, tilted to 40°, and the planet was struck by two Earth-mass impactors. In this case, even if Uranus was tilted initially by another method, the odds of generating Uranus's current spin state are about the same as if the planet was untilted. This is shown in Table 2, and the entries show similar likelihoods to the nonspinning case. However, both the nonspinning and slow-spinning cases are improbable for two reasons. First, the mechanism responsible for removing excess angular momentum during gas accretion needs to be extremely efficient. And second, the odds that both Uranus and Neptune were spun up similarly by impacts require significant fine-tuning.

Zoom In Zoom Out Reset image size

Table 2. An Initially Slow-rotating Uranus

N	Mi	MT	i	Probability (lU )	Normalized Probability
1	1.0	1.0	0°	4.5 × 10−3	0.90
2	0.2	0.5	0°	5.4 × 10−4	0.11
2	0.5	1.0	0°	1.0 × 10−2	2.00
2	1.0	2.0	0°	4.7 × 10−3	0.94
2	1.5	3.0	0°	2.5 × 10−3	0.50
1	1.0	1.0	40°	4.7 × 10−3	0.94
2	0.25	0.5	40°	9.0 × 10−4	0.18
2	0.5	1.0	40°	1.0 × 10−2	2.00
2	1.0	2.0	40°	5.0 × 10−3	1.00
2	1.5	3.0	40°	2.7 × 10−3	0.54
1	1.0	1.0	70°	4.8 × 10−3	0.96
2	0.25	0.5	70°	1.7 × 10−3	0.34
2	0.5	1.0	70°	1.0 × 10−2	2.00
2	1.0	2.0	70°	5.0 × 10−3	1.00
2	1.5	3.0	70°	2.7 × 10−3	0.54

Note. This table shows the same calculations as in Table 1, but with the planet having an initial spin period of 68.8 hr. _i is Uranus's initial obliquity. The normalized probability column divides the probability by 5 × 10⁻³ as in Table 1.

Download table as:  ASCIITypeset image

Accordingly, we investigate the effects of gas accretion by considering impacts onto an untilted fast-spinning (T_i = 17.2 hr) Uranus. Note that because we are adding angular momentum vectors, the order does not matter; therefore, striking Uranus with a giant impactor before the planet accretes gas will yield the same probability distributions as the reverse case considered here. For an initial spin period near Uranus's current value, the minimum impactor mass increases by $\sqrt{2}$ from ∼ 0.4 M_⊕ to 0.55 M_⊕ over the nonspinning case because the planet already has the correct ∣ L ∣, which must be rotated by ∼90° by the impact. However, while the slowly spinning cases have a relatively flat obliquity distribution, a fast-spinning planet is more resistant to change. For example, striking this planet with a 1 M_⊕ object will most likely yield little to no change to the planet's spin state (see Figure 7(a) in Rogoszinski & Hamilton 2020). Introducing more impactors does not change this conclusion appreciably; the planet still tends to remain with a low tilt and similar spin period. Figure 15 demonstrates this with the most favorable case of two 1 M_⊕ strikes onto an untilted planet already spinning with a 17.2 hr period. Additional cases are reported in Table 3.

Zoom In Zoom Out Reset image size

Table 3. An Initially Fast-rotating Uranus

N	Mi	MT	i	Probability (lU )	Normalized Probability
1	1.0	1.0	0°	3.4 × 10−3	0.68
2	0.25	0.5	0°	0	0
2	0.5	1.0	0°	3.7 × 10−3	0.74
2	1.0	2.0	0°	4.1 × 10−3	0.82
2	1.5	3.0	0°	2.6 × 10−3	0.52
5	0.6	3.0	0°	6.1 × 10−3	1.22
10	0.3	3.0	0°	7.5 × 10−3	1.50
15	0.2	3.0	0°	6.0 × 10−3	1.20
1	1.0	1.0	40°	4.5 × 10−3	0.90
2	0.25	0.5	40°	1.3 × 10−3	0.26
2	0.5	1.0	40°	7.4 × 10−3	1.48
2	1.0	2.0	40°	4.7 × 10−3	0.94
2	1.5	3.0	40°	2.6 × 10−3	0.52
1	1.0	1.0	70°	8.3 × 10−3	1.66
2	0.25	0.5	70°	2.6 × 10−2	5.20
2	0.5	1.0	70°	1.4 × 10−2	2.80
2	1.0	2.0	70°	5.7 × 10−3	1.14
2	1.5	3.0	70°	2.7 × 10−3	0.54

Note. This table shows the same calculations as in Table 1, but with the planet having an initial spin period of 17.2 hr. _i is Uranus's initial obliquity. The final column normalizes the probability column by 5 × 10⁻³ as in Table 1.

Download table as:  ASCIITypeset image

If Uranus was initially tilted by a 40° resonance kick, its rapid rotation ensures that its spin state will tend to remain relatively unaffected by subsequent impacts. This can be seen in Figure 16(a) with a 1 M_⊕ strike, where the probability of tilting Uranus to 98° is only 4.5 × 10⁻³. The odds do improve if the number of impacts increases (Figure 16(b)), but they are not better than the nonspinning case. However, if Uranus was initially tilted by 70° via a spin–orbit resonance (Rogoszinski & Hamilton 2020), then two 0.5 M_⊕ strikes generate a favorable result (Figure 16(c)). Also, only in this case will two 0.25 M_⊕ strikes yield even better likelihoods (see Figure 8 in Rogoszinski & Hamilton 2020). Therefore, if Uranus's and Neptune's current spin rates were a by-product of gas accretion, then a large resonance kick can significantly reduce the mass needed in later impacts.

Zoom In Zoom Out Reset image size

Finally, the mechanism that removes angular momentum during gas accretion could have been very weak, and Uranus would have been initially spinning very fast (T_i ≫ 17.2 hr). In this case, slowing down Uranus's spin rate and tilting the planet over would require very massive impacts. As discussed in the previous subsection, changing the planet's spin state with many impactors requires more impacting mass to compensate for the partial cancellations of impact effects. Table 4 shows that 10 impacts totaling to 4 M_⊕ produce plausible outcomes. While their obliquity distributions peak at around 30°, which favors a Neptune formation scenario, the planets would still likely be spinning twice as fast as they are today (Figure 17). We would therefore not expect the required massive impactors to spin down both Uranus and Neptune similarly. Additionally, 10 independent strikes is less probable than only 2, while also requiring the solar system to have been populated with many massive rogue planetary cores.

Zoom In Zoom Out Reset image size

Table 4. An Initially Very Fast-rotating Uranus

N	Mi	MT	i	Probability (lU )	Normalized Probability
1	1.0	1.0	0°	2.3 × 10−3	0.46
2	0.25	0.5	0°	0	0.00
2	0.5	1.0	0°	2.6 × 10−4	0.05
2	1.0	2.0	0°	2.7 × 10−3	0.54
2	1.5	3.0	0°	2.0 × 10−3	0.40
1	1.0	1.0	40°	4.1 × 10−3	0.82
2	0.25	0.5	40°	0	0
2	0.5	1.0	40°	2.0 × 10−3	0.40
2	1.0	2.0	40°	4.1 × 10−3	0.82
2	1.5	3.0	40°	2.5 × 10−3	0.50
1	1.0	1.0	70°	2.1 × 10−3	0.42
2	0.25	0.5	70°	1.2 × 10−4	0.02
2	0.5	1.0	70°	3.3 × 10−3	0.66
2	1.0	2.0	70°	3.0 × 10−3	0.60
2	1.5	3.0	70°	2.4 × 10−3	0.48
5	0.8	4.0	0°	3.4 × 10−3	0.68
10	0.4	4.0	0°	5.0 × 10−3	1.00
15	0.266 7	4.0	0°	4.4 × 10−3	0.88

Note. This table shows the same calculations as the previous tables, but the planet is spinning with a period of 8.6 hr. The final column has been normalized as in Table 1.

Download table as:  ASCIITypeset image

The calculation for the angular momentum distribution of a planet from multiple strikes at random locations on the planet's surface is a random walk scenario. We start with the magnitude of the spin angular momentum of a planet:

$$\begin{eqnarray}&&L=\sqrt{{L}_{X}^{2}+{L}_{Y}^{2}+{L}_{Z}^{2}},\end{eqnarray} \tag{ A1 }$$

where the probability distribution (${f}_{{L}_{k}}({l}_{k})$) of each component (L_k ) of the angular momentum vector is described by a normal distribution as a by-product of the central limit theorem:

$$\begin{eqnarray}&&{f}_{{L}_{k}}({l}_{k})=\displaystyle \frac{1}{{\sigma }_{k}\sqrt{2\pi }}{e}^{-{l}_{k}^{2}/2{\sigma }_{k}^{2}}.\end{eqnarray} \tag{ A2 }$$

As such, to find the distribution of the magnitude of the angular momentum, we will first need to determine the square of each distribution and then the sum of three squares, and finally, to take the square root of the sum as seen in Equation (A1).

The distribution of the square of each component (L_k ²) can be calculated by assuming that X and Y are continuous random variables (i.e., "variates" as depicted in upper case), with x and y as specific elements in the ranges of their corresponding variates (i.e., also called "quantiles," depicted here in lower case; Grinstead & Snell 2006). X and Y have cumulative distribution functions F_X and F_Y , and Y is described by a strictly increasing function of X: Y = ϕ(X). F_Y (y) = P(Y ≤ y), where the right-hand side describes the probability that the variate Y is less than or equal to a number y, which is equal to P(ϕ(X) ≤ y) = P(X ≤ ϕ⁻¹(y)) = F_X (ϕ⁻¹(y)).

So, for the variate X² and its corresponding quantile x²:

$$\begin{eqnarray}&&{F}_{{X}^{2}}({x}^{2})=P({X}^{2}\leqslant {x}^{2})=P(-x\leqslant X\leqslant x).\end{eqnarray} \tag{ A3 }$$

The right-hand side can be rearranged accordingly:

$$\begin{eqnarray}&&P(-x\leqslant X\leqslant x)=P(X\leqslant x)-P(X\leqslant -x),\end{eqnarray} \tag{ A4 }$$

so that

$$\begin{eqnarray}&&{F}_{{X}^{2}}({x}^{2})={F}_{X}(x)-{F}_{X}(-x).\end{eqnarray} \tag{ A5 }$$

The corresponding density distribution function for an arbitrary variate Y is ${f}_{Y}(y)=\tfrac{d}{{dy}}{F}_{Y}(y)$. Starting with F_Y (y) = F_X (ϕ⁻¹(y)), we take the derivative of each side and employ the chain rule to obtain ${f}_{Y}(y)={f}_{X}({\phi }^{-1}(y))\tfrac{d}{{dy}}{\phi }^{-1}(y)$. So,

$$\begin{eqnarray}&&{f}_{{X}^{2}}({x}^{2})=\displaystyle \frac{{f}_{X}(x)+{f}_{X}(-x)}{2x}.\end{eqnarray} \tag{ A6 }$$

Because the normal distribution is centered at zero and is symmetric, the density distribution for L_k ² is then

$$\begin{eqnarray}&&{f}_{{L}_{k}^{2}}({l}^{2})=\displaystyle \frac{1}{{\sigma }_{k}l\sqrt{2\pi }}{e}^{-{l}^{2}/2{\sigma }_{k}^{2}},\end{eqnarray} \tag{ A7 }$$

which is the distribution for a chi squared with one degree of freedom.

Next, the density distribution of the sum of two independent random variables is their convolution. Let ${L}_{{XY}}^{2}={L}_{x}^{2}+{L}_{y}^{2}$ and its corresponding density distribution:

$$\begin{eqnarray}&&{f}_{{L}_{{XY}}^{2}}({l}_{{xy}}^{2})={\int }_{0}^{{l}_{{xy}}^{2}}{f}_{{L}_{X}^{2}}({l}_{{xy}}^{2}-{l}_{y}^{2}){f}_{{L}_{Y}^{2}}({l}_{y}^{2}){{dl}}_{y}^{2},\end{eqnarray} \tag{ A8 }$$

where L_Y ² ranges from 0 to L_XY ². Note that the standard deviations for both ${f}_{{L}_{X}}$ and ${f}_{{L}_{Y}}$ are equal with σ = σ_x = σ_y . Thus, combining Equations (A7) and (A8),

$$\begin{eqnarray}\begin{array}{rcl}{f}_{{L}_{{XY}}^{2}}({l}_{{xy}}^{2}) & = & \displaystyle \frac{1}{2\pi {\sigma }^{2}}{\displaystyle \int }_{0}^{{l}_{{xy}}^{2}}\left({e}^{-({l}_{{xy}}^{2}-{l}_{y}^{2})/2{\sigma }^{2}}{\left({l}_{{xy}}^{2}-{l}_{y}^{2}\right)}^{-0.5}\right)\\ & & \times \left({e}^{-{l}_{y}^{2}/2{\sigma }^{2}}{\left({l}_{y}^{2}\right)}^{-0.5}\right){{dl}}_{y}^{2}=\displaystyle \frac{{e}^{-{l}_{{xy}}^{2}/2{\sigma }^{2}}}{2{\sigma }^{2}}.\end{array}\end{eqnarray} \tag{ A9 }$$

Now, let ${L}^{2}={L}_{{XY}}^{2}+{L}_{Z}^{2}$ and repeat the above process. The probability distribution ${f}_{{L}^{2}}({l}^{2})$ describes the probability that L² takes the value l², and ${f}_{{L}_{Z}^{2}}({l}_{z}^{2})$ describes the probability that L_Z ² takes the value l_z ². We explicitly treat the general case σ ≠ σ_z .

The density distribution for L² is

$$\begin{eqnarray}\begin{array}{rcl}{f}_{{L}^{2}}({l}^{2}) & = & {\displaystyle \int }_{0}^{{l}^{2}}{f}_{{L}_{{XY}}^{2}}({l}^{2}-{l}_{z}^{2}){f}_{{L}_{Z}^{2}}({l}_{z}^{2}){{dl}}_{z}^{2}\\ & = & \displaystyle \frac{1}{2\sqrt{2\pi }\,{\sigma }^{2}{\sigma }_{z}}{\displaystyle \int }_{0}^{{l}^{2}}{e}^{-({l}^{2}-{l}_{z}^{2})/2{\sigma }^{2}}{e}^{-{l}_{z}^{2}/2{\sigma }_{z}^{2}}\,{\left({l}_{z}^{2}\right)}^{-0.5}\,{{dl}}_{z}^{2},\end{array}\end{eqnarray} \tag{ A10 }$$

let $\beta =({\sigma }^{2}-{\sigma }_{z}^{2})/2{\sigma }^{2}{\sigma }_{z}^{2}$, $\gamma =\beta {l}_{z}^{2}$, and $d\gamma =\beta {{dl}}_{z}^{2}$, and so

$$\begin{eqnarray}&&{f}_{{L}^{2}}({l}^{2})=\displaystyle \frac{{e}^{-{l}^{2}/2{\sigma }^{2}}}{2\sqrt{2\pi }\,{\sigma }^{2}{\sigma }_{z}}\displaystyle \frac{1}{\sqrt{\beta }}{\int }_{0}^{\beta {l}^{2}}{e}^{-\gamma }\,{\gamma }^{-0.5}\,d\gamma .\end{eqnarray} \tag{ A11 }$$

Equation (A11) is of similar form to Equation (109) found in Dones & Tremaine (1993a). Applying Equation (A6) to ${f}_{{L}^{2}}$ and noting that because L is the magnitude of the planet's angular momentum, f_L ( − l) = 0. We find ${f}_{L}(l)={f}_{{L}^{2}}(l)\cdot 2l$. The probability distribution describing the angular momentum of the planet for β > 0, or σ_x = σ_y > σ_z is then

$$\begin{eqnarray}&&{f}_{L}(l)=\displaystyle \frac{{{le}}^{-{l}^{2}/2{\sigma }^{2}}}{\sqrt{2\pi }\,{\sigma }^{2}{\sigma }_{z}}\displaystyle \frac{1}{\sqrt{\beta }}\,\gamma (0.5,\beta {l}^{2}),\end{eqnarray} \tag{ A12 }$$

where γ(0.5, β l²) is the lower incomplete gamma function. For β < 0 (σ_x = σ_y < σ_z ):

$$\begin{eqnarray}&&{f}_{L}(l)=\displaystyle \frac{{{le}}^{-{l}^{2}/2{\sigma }^{2}}}{\sqrt{2\pi }\,{\sigma }^{2}{\sigma }_{z}}\displaystyle \frac{1}{\sqrt{-\beta }}(2l\sqrt{-\beta }){\rm{\Phi }}(0.5;1.5;-\beta {l}^{2})\end{eqnarray} \tag{ A13 }$$

where Φ(0.5; 1.5; − β l²) is the confluent hypergeometric function of the first kind. For β = 0, where σ = σ_x = σ_y = σ_z (isotropic case), the form is particularly simple:

$$\begin{eqnarray}&&{f}_{L}(l)=\displaystyle \frac{2{l}^{2}{e}^{-{l}^{2}/2{\sigma }^{2}}}{\sqrt{2\pi }\,{\sigma }^{3}}.\end{eqnarray} \tag{ A14 }$$

The obliquity angle () is defined by $\tan (\epsilon )=\sqrt{{L}_{x}^{2}+{L}_{y}^{2}}/{L}_{z}={L}_{{XY}}/{L}_{z}$. To find the distribution of the quotient of two independent variants, we let Q = X/Y, where X and Y are independent random variables. Then, F_Q (q) = P(Q ≤ q) = P(X/Y ≤ q). If Y > 0, then X ≤ yq, while if Y < 0, then X ≥ yq. Therefore, P(X/Y ≤ q) = P(X ≤ yq, Y > 0) + P(X ≥ yq, Y < 0). These constraints determine the integral limits in the corresponding cumulative distribution:

$$\begin{eqnarray}\begin{array}{rcl}{F}_{Q}(q) & = & {\displaystyle \int }_{y=0}^{\infty }{\displaystyle \int }_{x=-\infty }^{{yq}}{f}_{{XY}}(x,y){dxdy}\\ & & +{\displaystyle \int }_{y=-\infty }^{0}{\displaystyle \int }_{x={yq}}^{\infty }{f}_{{XY}}(x,y){dxdy},\end{array}\end{eqnarray} \tag{ A15 }$$

and density distribution:

$$\begin{eqnarray}&&{f}_{Q}(q)={\int }_{0}^{\infty }{{yf}}_{{XY}}({yq},y){dy}+{\int }_{-\infty }^{0}(-y){f}_{{XY}}({yq},y){dy}.\end{eqnarray} \tag{ A16 }$$

So, to calculate the obliquity distribution, let σ ≠ σ_z , and $U=\tan (\epsilon )$ with $u=\tan (\varepsilon )$ as the corresponding quantile. Thus,

$$\begin{eqnarray}\begin{array}{rcl}{f}_{U}(u) & = & {\displaystyle \int }_{0}^{\infty }{l}_{z}{f}_{{L}_{{XY}}}({{ul}}_{z}){f}_{{L}_{z}}({l}_{z}){{dl}}_{z}\\ & & +{\displaystyle \int }_{-\infty }^{0}-{l}_{z}{f}_{{L}_{{XY}}}({{ul}}_{z}){f}_{{L}_{z}}({l}_{z}){{dl}}_{z},\end{array}\end{eqnarray} \tag{ A17 }$$

which becomes

$$\begin{eqnarray}&&{f}_{U}(u)=2{\int }_{0}^{\infty }\displaystyle \frac{| u| {l}_{z}^{2}}{{\sigma }^{2}{\sigma }_{z}\sqrt{2\pi }}{e}^{-{l}_{z}^{2}{u}^{2}/(2{\sigma }^{2})}{e}^{-{l}_{z}^{2}/(2{\sigma }_{z}^{2})}{{dl}}_{z}.\end{eqnarray} \tag{ A18 }$$

If we let $\alpha =\tfrac{{u}^{2}}{2{\sigma }^{2}}+\tfrac{1}{2{\sigma }_{z}^{2}}$, then the equation is now of the form

$$\begin{eqnarray}&&{\int }_{0}^{\infty }{t}^{2}{e}^{-\alpha {t}^{2}}{dt}=\displaystyle \frac{\sqrt{\pi }}{4{\alpha }^{1.5}},\end{eqnarray} \tag{ A19 }$$

and so when normalized:

$$\begin{eqnarray}&&{f}_{U}(u)=\left|\displaystyle \frac{u}{4\sqrt{2}\,{\sigma }^{2}{\sigma }_{z}{\alpha }^{1.5}}\right|.\end{eqnarray} \tag{ A20 }$$

We can change variables to obliquity () by setting ${f}_{\epsilon }(\varepsilon )\,=\tfrac{{du}}{d\varepsilon }{f}_{U}(u)$ where $\tfrac{{du}}{d\varepsilon }={\sec }^{2}(\varepsilon )$. We find

$$\begin{eqnarray}&&{f}_{\epsilon }(\varepsilon )=\left|\displaystyle \frac{1}{4\sqrt{2}\,{\sigma }^{2}{\sigma }_{z}}\displaystyle \frac{\tan (\varepsilon )}{{\cos }^{2}(\varepsilon )}{\left(\displaystyle \frac{{\tan }^{2}(\varepsilon )}{2{\sigma }^{2}}+\displaystyle \frac{1}{2{\sigma }_{z}^{2}}\right)}^{-3/2}\right|.\end{eqnarray} \tag{ A21 }$$

This is equivalent to the obliquity distribution given in Dones & Tremaine (1993a, Equation (111)). For the isotropic case, σ_z = σ, the distribution reduces to

$$\begin{eqnarray}&&{f}_{\epsilon }(\varepsilon )=\left|\displaystyle \frac{1}{2}\sin (\varepsilon )\right|.\end{eqnarray} \tag{ A22 }$$