Spin Dynamics of Extrasolar Giant Planets in Planet–Planet Scattering

The Mercury code is modified to include planet oblateness by modifying the force term generated by the gravitational field around the planet. The gravitational field has a quadruple moment term J₂ = 0.0147, consistent with Jupiter and slightly less than Saturn's. All simulated planets adopt this oblateness. The J₂ moment modifies the spherically symmetric force field around a point mass, to an axisymmetric force field around an extended mass. The code assumes zero reaction time of planet spin to external perturbations, but in reality there might be a time lag due to viscous dissipation for a fluid body (such as a giant planet) to react to a torque. In general, the internal dynamics of a precessing fluid body is rather complicated, potentially involving turbulence and strong dissipation (Papaloizou & Pringle 1982; Barker 2016; Lin & Ogilvie 2017). However, since the internal dissipation conserves angular momentum, we expect that as far as forced precession is concerned, the giant planet effectively behaves as a rigid body to a leading-order approximation. In any case, as in all previous works on planetary obliquities, we will ignore any complications associated with internal fluid motions in this paper.

The spin of an oblate planet evolves under the torque from the central star and other planets. In planetary close encounters, the orbital parameters of a planet change rapidly, and the spin of a planet can receive large torques from another planet, in a shorter period of time than their orbital periods. Therefore, this work adopts an instantaneous equation of motion for the spin, instead of the secular one commonly used.

For an extended mass under the gravitational influence of a point mass, different parts of it receive differing amount of gravitational attraction. The difference in gravitational attraction in different parts of the extended mass causes a torque, leading to spin evolution. The force on a point mass M from an extended mass m is

$$\begin{eqnarray}&&\begin{array}{l}\vec{F}=\displaystyle \frac{-\,G\,M\,m}{{r}^{2}}\left\{\hat{r}-\displaystyle \frac{3{J}_{2}}{2}{\left(\displaystyle \frac{R}{r}\right)}^{2}\,[(5{(\hat{r}\cdot \hat{s})}^{2}-1)\hat{r}\right.\\ \qquad \qquad \left.-2(\hat{r}\cdot \hat{s})\hat{s}]+\cdots \right\},\end{array}\end{eqnarray} \tag{ 1 }$$

where G is the gravitational constant, r the distance between the centers of m and M, $\hat{r}$ the unit vector pointing from m to M, R the radius of m, and $\hat{s}$ the normalized spin vector of m (e.g., Hilton 1991).

The reaction torque per moment of inertia on the extended body is

$$\begin{eqnarray}&&\displaystyle \frac{-{\boldsymbol{r}}\times {\boldsymbol{F}}}{I}=\displaystyle \frac{d}{{dt}}\left({\rm{\Omega }}\,\hat{s}\right),\end{eqnarray} \tag{ 2 }$$

where Ω is the rotation rate of the planet, which is set to be 2π/(10 hr). The rotation period is close to that of Jupiter. For simplicity, the rotation rate is assumed to remain constant throughout the simulations.

Combining the two equations above, we can obtain the instantaneous equation of motion for the spin components of mass m:

$$\begin{eqnarray}&&\displaystyle \frac{d\hat{s}}{{dt}}=\displaystyle \frac{3\,G\,M}{{r}^{3}}\displaystyle \frac{{J}_{2}}{\lambda \,{\rm{\Omega }}}\left(\hat{r}\cdot \hat{s}\right)\left(\hat{r}\times \hat{s}\right),\end{eqnarray} \tag{ 3 }$$

where $\lambda =\tfrac{I}{{{mR}}^{2}}$ is the normalized moment of inertia. λ is set to 0.25, close to that of Jupiter, in all simulations.

The sources of torque on planet spin include the central star and other planets. The implementation of the above equation of motion into Mercury is verified with the past evolution of Mars's obliquity. The simulation is checked against the result in Ward (1973) using the initial conditions provided in Quinn et al. (1991), with all planets and Pluto included in the simulation. An initial obliquity of 2518 is applied to Mars. Although his approach was based on secular equations, our result for the obliquity evolution closely resembles Ward (1973) (Figure 1)—this is expected since no close encounters or strong scatterings are involved in the evolution. The integration error limit of the simulations is set at 10⁻¹² per time step, and the initial time step is 1 day. The integrator automatically adjusts the time step in accordance with the set error limit.

Zoom In Zoom Out Reset image size

Before presenting our numerical result, we first summarize spin dynamics as numerous previous works on planetary/stellar obliquities are based on secular theory (e.g., Ward 1973; Laskar & Robutel 1993; Ward & Hamilton 2004; Storch et al. 2014; Storch & Lai 2015; Anderson et al. 2016; Millholland & Batygin 2019; Millholland & Spalding 2020; Su & Lai 2020).

In secular systems, when a planet's own spin and orbital angular momenta are under perturbation, $| \tfrac{d\hat{s}}{{dt}}|$ and $| \tfrac{d\hat{L}}{{dt}}|$ evolve. The obliquity evolution can be categorized into three different regimes. When $| \tfrac{d\hat{s}}{{dt}}|$ (rate of change of the direction of spin momentum) is much greater than $| \tfrac{d\hat{L}}{{dt}}|$ (rate of change of the direction of orbital angular momentum), the spin axis of the planet will follow its own orbit normal, and the planet will retain a constant obliquity. If $| \tfrac{d\hat{s}}{{dt}}|$ is much smaller than $| \tfrac{d\hat{L}}{{dt}}|$, the spin of the planet does not follow its own orbit normal. Instead, it precesses with a constant angular separation around the normal to system's invariable plane. In two-planet systems, the spin of a planet precesses around the total orbital angular momentum axis of the two planets. When $| \tfrac{d\hat{s}}{{dt}}| \sim | \tfrac{d\hat{L}}{{dt}}|$, resonance crossing may lead to significant spin evolution, and can generate high obliquity planets. Eccentricity evolution in planet–planet scattering can also alter the rate of spin evolution, because it alters r-dependent terms in Equation (3). Therefore, eccentricity evolution in planet–planet scattering can also contribute to resonance crossings (e.g., Storch et al. 2014).

While the secular equations do not strictly apply to our problem, they can be used to provide a qualitative understanding of how a planet may attain a significant obliquity during the scattering process. Consider planet m₁ on an approximately circular orbit (with semimajor axis a₁) while planet m₂ on an eccentric orbit (with a₂ and e₂) interact with m₁ near its pericenter (thus a₂(1 − e₂) ∼ a₁). The spin axis of m₁ precesses around its orbital axis ${\hat{L}}_{1}$ (driven by the torque from the host star) at a rate given by

$$\begin{eqnarray}&&{\omega }_{s}=\displaystyle \frac{3{k}_{q}}{2\lambda }\displaystyle \frac{{M}_{\star }}{{m}_{1}}{\left(\displaystyle \frac{{R}_{1}}{{a}_{1}}\right)}^{3}{{\rm{\Omega }}}_{1},\end{eqnarray} \tag{ 4 }$$

where ${k}_{q}={J}_{2}({{Gm}}_{1}/{R}_{1}^{3})/{{\rm{\Omega }}}_{1}^{2}$ (with R₁ and Ω₁ the radius and rotation rate of m₁). On the other hand, the orbital axis ${\hat{L}}_{1}$ precesses around the total angular momentum axis ${\vec{L}}_{1}+{\vec{L}}_{2}$ at the rate

$$\begin{eqnarray}&&{\omega }_{l}\simeq \displaystyle \frac{3{m}_{2}}{4{M}_{\star }}{\left(\displaystyle \frac{{a}_{1}}{{a}_{2}}\right)}^{3}\displaystyle \frac{{n}_{1}}{{(1-{e}_{2}^{2})}^{3/2}},\end{eqnarray} \tag{ 5 }$$

where n₁ is the mean motion of m₁. Equation (5) assumes a₂(1 − e₂) ≫ a₁ (e.g., Liu et al. 2014); when this is not satisfied, the actual ω_l can be larger by a factor of a few. The ratio of the two precession frequencies is

$$\begin{eqnarray}\begin{array}{rcl}\eta & = & \displaystyle \frac{{\omega }_{l}}{{\omega }_{s}}\simeq \displaystyle \frac{3.8\lambda }{{k}_{q}}\left(\displaystyle \frac{{\rho }_{1}}{{\rm{g}}\,{\mathrm{cm}}^{-3}}\right){\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{-3/2}\left(\displaystyle \frac{{m}_{2}}{{M}_{J}}\right)\left(\displaystyle \frac{{P}_{1}}{10\,\mathrm{hr}}\right)\\ & & \times \,{\left[\displaystyle \frac{{a}_{1}}{{a}_{2}(1-{e}_{2})}\right]}^{3}{\left(\displaystyle \frac{{a}_{1}}{\mathrm{au}}\right)}^{3/2}{\left(\displaystyle \frac{1-{e}_{2}}{1+{e}_{2}}\right)}^{3/2},\end{array}\end{eqnarray} \tag{ 6 }$$

where P₁ = 2π/Ω₁ and ρ₁ are the rotation period and mean density of of the planet, and λ ∼ k_q for giant planets.

During planet–planet scatterings that lead to planet ejection, a₂ and e₂ increase (while a₂(1 − e₂) remains close to a₁) as a result of energy exchange between m₁ and m₂, in a process analogous to diffusion (see Pu & Lai 2021). Thus the ratio η can pass through ∼1 in a dynamical way. This "transient" resonance crossing may potentially give rise to appreciable obliquity excitation of the planet.

In planet–planet scattering, the spin of a planet evolves under torques from the host star and other planets in the system. Because of the short timescale evolution of planetary orbits, and very short duration of planetary close encounters, we use the nonsecular equation for the planet's spin evolution (Equation (1)). Here we discuss how planets can experience significant spin evolution and gain high obliquity in planet–planet scattering. During the planet–planet scattering phase, close planetary encounters, as well as strong mutual perturbation outside of planetary close encounters, can change the planet's orbital elements by large amounts over short periods, which in turn can change $| \tfrac{d\hat{L}}{{dt}}|$ and $| \tfrac{d\hat{s}}{{dt}}|$. This can drive the system close to a state where $| \tfrac{d\hat{L}}{{dt}}| \simeq | \tfrac{d\hat{s}}{{dt}}|$, and thus lead to large variations in the obliquity due to stellar torques.

Figures 2 and 3 show an example of the spin and orbit evolution of a two-planet system in the base simulation set for the first two million years. Figure 3 shows typical orbital evolution of planets during the planet–planet scattering phase. The spin of the planets evolves chaotically. The outer planet's spin axis evolves beyond 80° from its initial orientation at 2 million years. The inner planet's spin evolution is relatively stable. Figure 2 also shows the corresponding evolution of inclination and obliquity. The dramatic spin evolution of the outer planet leads to its high obliquity (max = 90°), under a relatively modest inclination. The inner planet experiences very little spin evolution, with a maximum displacement of 14. The contrast in the scale of spin evolution between the outer and inner planets can be explained by considering the three regimes laid out in the last subsection. When the planet's orbital precession rate $| \tfrac{d\hat{L}}{{dt}}|$ is close to $| \tfrac{d\hat{s}}{{dt}}|$, resonance effects can cause the spin to significantly evolve. In the top panel of Figure 4, the ratio between $| \tfrac{d\hat{L}}{{dt}}|$ and $| \tfrac{d\hat{s}}{{dt}}|$ for the outer planet lingers around 1 very frequently (${\mathrm{log}}_{10}| \tfrac{{\rm{d}}\hat{s}}{\mathrm{dt}}| /| \tfrac{{\rm{d}}\hat{L}}{\mathrm{dt}}| \sim 0$), so its spin evolves significantly. $| \tfrac{d\hat{L}}{{dt}}|$ and $| \tfrac{d\hat{s}}{{dt}}|$ have chances to be close to each other because in planet–planet scattering, both can evolve significantly, as shown in the middle and bottom panels of Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the inner planet, $| \tfrac{d\hat{L}}{{dt}}| /| \tfrac{d\hat{s}}{{dt}}|$ stays much higher than 1 for the most part, so its spin evolves very little in the example depicted in Figure 2. However, the inner (remnant) planet does have a chance to attain larger obliquities, as revealed in our ensemble of simulations (see Section 4.3).

This section discusses the results of planet–planet scattering in systems with two giant planets in four simulation sets (Table 1). The statistics use only simulations with one surviving planet after the other planet is ejected via planet–planet scattering (23.26% of all simulations; see Table 1), which is around 5500 simulations in each set. In two-planet systems, the instability boundary is very narrow. If both planets survive, its very likely they have none or very little interactions, so it is not very interesting to study them for our purpose. Simulations with planet–star (2.71% of all simulations) or planet–planet collisions (73.61% of all simulations; see Table 1) are excluded from the statistics, as obliquity evolution for the collision of fluid bodies is beyond the scope of this work (see Li & Lai 2020; Li et al. 2020). Simulations with energy error $| \tfrac{\delta E}{E}|$ greater than ${ \mathcal O }({10}^{-4})$ are also excluded (0.36% of all simulations), which is an adequate threshold for multiplanet systems (Barnes & Quinn 2004; Raymond et al. 2010). A negligible fraction of systems (0.03% of all simulations) that do not experience planet destruction in 10⁸ yr are also excluded. Figures 5 and 6 show the distribution of the surviving planet's final spin orientation and obliquity for simulation sets 1–3. The spin inclination is measured relative to the z-axis, which is aligned with the inner planet's initial spin and orbit normal. The final spin and obliquity distribution is quite broad compared with the orbital inclination distribution (Figure 7). Overall, the final spin and obliquity distributions (Table 3) are broader than the inclination distributions, with 19.51% planets' obliquity beyond 40°, and 8.77% of the surviving planets are on retrograde obliquity (i.e., the rotation is in reversed direction relative to the orbit). The inclination distribution only has 0.12% planets beyond 40°. The rate of highly oblique planets (>40°) are 17.86%, 31.98%, and 10.43% in sets 1, 2, and 3. There are 5.07% surviving planets on retrograde obliquity in the base set, 21.46% in set 2, and 3.21% in set 3. The distribution of obliquity for set 3 is more clustered at lower regions compared to the other sets. From the rate of high and retrograde obliquity in set 2 (the inner planet ends up at 0.19 au after the outer planet is ejected), we may expect a high rate of highly oblique planets close to the central star, and some oblique planets farther out (Figures 5 and 6). A moderate peak in Figure 6 for close-in planets in sets 1 and 2 between obliquities 170°–180° implies that Venus-like obliquity is possible. The simulation results imply that planets closer to the star receive more torque from it and become more evolved in obliquity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The final spin distributions of the three simulation sets all fit best with a log-normal distribution (Figure 8):

$$\begin{eqnarray}&&f(s,\sigma )=\displaystyle \frac{1}{\sqrt{2\pi }\sigma }{e}^{-\tfrac{1}{2}{\left(\tfrac{{ln}(s)-\mu }{\sigma }\right)}^{2}},\end{eqnarray} \tag{ 7 }$$

σ is the standard deviation, and μ the mean of the log-normal distribution. We do not have an explanation for the log-normal distribution, but we find it provides a good fit—this is a finding from our numerical experiments. Table 4 shows the statistics of the log-normal fit.

Zoom In Zoom Out Reset image size

For completeness, the Kolmogorov–Smirnov (K-S) test shows that the final eccentricity distributions between the three simulation sets are similar, with K-S statistics 0.08–0.10 and p-values 0.96–1.0 (Figure 9). The low value of K-S statistics implies the distributions are similar, and the p-value is above our set 95% confidence level, so one cannot reject the null hypothesis that the distributions are the same. This result agrees with previous findings (Li et al. 2020). For the total number of close encounters, close-in planets in set 2 experience more than the other 2 sets, and set 1 and set 3 have nearly identical distributions as shown by the K-S test. We have also found that swapping the masses of the inner planet with the outer planet does not produce a significant difference in the distributions of final obliquities and spin inclinations, as tested for simulation set 3.

Zoom In Zoom Out Reset image size

Simulation set 4 has the same initial setting as the base set, but torques between planets are turned off. Figure 10 shows the histograms of final spin inclinations for the base set and set 4. The distributions have a K-S statistic of 0.023 and p-value of 1.00. The statistics cannot reject the null hypothesis that the two distributions are the same, and the statistical evidence is very strong that planet–planet torque is very insignificant. As shown in Figure 11, we can see that the stellar torque is greater than planet–planet torque most of the time. Although the torque between planets is very large during close encounters (at the first blue peak ∼20 yr, the torque from the inner planet is ∼10⁵ greater than the stellar torque), due to their short duration and small mass of planets relative to the star, statistically the summed effect of planet–planet torque is not very important. In the simulation shown in Figure 11, the summed torque on the outer planet from the star is 3 × 10⁴ greater than from the inner planet in the first 2 million years of simulation time, so planet–planet torque is rather insignificant compared to stellar torque. The results show that planet–planet torque does not statistically affect simulation results with regards to spin or obliquity (Lee et al. 2007), and helps conclude that the torque from the central star is the main driver for the spin and obliquity evolution of the planets.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size