Exploring Tidal Obliquity Variations with SMERCURY-T

Steven M. Kreyche; Jason W. Barnes; Billy L. Quarles; John E. Chambers

doi:10.3847/PSJ/ac1ce6

1. Introduction

Planetary obliquity is a fundamental parameter that factors into the assessment of a planet's habitability. The obliquity of a planet, or axial tilt, is the angle between the planet's spin axis and orbit normal. Together, the value of the obliquity and its evolution over time determine the nature of the planet's seasons (Williams & Kasting 1997; Williams & Pollard 2003; Kane & Torres 2017; Kilic et al. 2017; Kang 2019), atmospheric and oceanic processes (Kilic et al. 2018; Dong et al. 2019; Guendelman & Kaspi 2019; Olson et al. 2020), and ultimately, its climatic stability (Milanković 1998; Armstrong et al. 2014; Spiegl et al. 2015; Deitrick et al. 2018b; Colose et al. 2019). Various studies investigated the obliquity evolution of bodies in the solar system, as well as exoplanets, for purposes such as gaining insights into their dynamical history, making predictions for the future, or even judging their potential for habitability.

Many processes can contribute to alter a planet's obliquity over time. Constant torques exerted on a planet's equatorial bulge by its planetary neighbors in addition to the regular precessional motion of the planet's orbit plane set a range of benign obliquity variations. On the other hand, a planet can experience large obliquity variations of many tens of degrees if it enters into a secular spin–orbit resonance. Such a resonance occurs when a planet's spin-axis precession becomes commensurate with a driving eigenfrequency of its orbital precession (Colombo 1966; Peale 1969; Henrard & Murigande 1987). The resultant obliquity variations of this resonance depend on the amplitude of the orbital mode in addition to the planet's proximity to the resonance center. Provided that a planet resides within a spin–orbit resonance, its resonant angle, σ = ψ + ϕ, will typically librate about 0° (stable Cassini states 1 and 2 or unstable Cassini state 4) or 180° (Cassini state 2, stable), where ψ is the planet's spin-precession angle that defines the azimuthal orientation of the planet's spin vector and ϕ is the phase angle of the orbital precession mode (Deitrick et al. 2018a; Saillenfest et al. 2019; Su & Lai 2020). Furthermore, in the case where the widths of multiple spin–orbit resonances overlap, a planet's spin axis can be subject to chaotic evolution in which its obliquity swings across a wide range of values. Mars is the poster child for this chaotic behavior, as its obliquity swings from ∼0° to 60° (Laskar & Robutel 1993b; Touma & Wisdom 1993; Laskar et al. 2004b; Li & Batygin 2014). The Earth would suffer a similar fate if not for the Moon's help in speeding up the Earth's spin-axis precession enough so that the Earth peacefully resides outside of the chaotic region (Laskar et al. 1993a, 2004a; Brasser & Walsh 2011).

Considering the working parts that contribute to this resonance, rhythmic gravitational tugs from neighboring planets drive orbital precession, while torques from the planet's star and moons fuel its spin precession. The rate of a planet's spin precession obeys

$\begin{eqnarray}&&\dot{\psi }=\displaystyle \frac{\alpha \cos \epsilon }{{\left(1-{e}^{2}\right)}^{\tfrac{3}{2}}},\end{eqnarray} \tag{ 1 }$

where is the planet's obliquity, and e is its orbital eccentricity. In the case of a moonless planet, the precession constant, α, goes as

$\begin{eqnarray}&&\alpha \approx \displaystyle \frac{3{n}^{2}}{2{\omega }_{p}}\displaystyle \frac{{J}_{2}}{\overline{C}}.\end{eqnarray} \tag{ 2 }$

Here n is the planet's mean motion, ω_p is its spin rate, J₂ is its oblateness coefficient, and $\overline{C}$ is its normalized polar moment of inertia (this expression is approximate, based on the assumption that the planet's dynamical ellipticity is roughly proportional to ${{\omega }_{p}}^{2}$ ; Laskar & Robutel 1993b; Li & Batygin 2014; Quarles et al. 2020).

The status of a planet's presence within a spin–orbit resonance can change with time. For instance, although the Earth rests easy now, it may have crossed a high-order resonance just a few million years ago (Laskar et al. 1993a, 2004a). Change can be brought on by a variety of processes, such as those that shift the location and amplitudes of potentially resonant orbital frequencies. Planetary orbital migration mechanisms are one such way to achieve these shifts (Brasser & Lee 2015; Vokrouhlický & Nesvorný 2015; Millholland & Laughlin 2019). In a similar vein, the works of Saillenfest et al. (2020, 2021a, 2021b) studied the past and future influence that satellite migration may have had on the obliquity evolution of Jupiter and Saturn. Millholland & Batygin (2019) and Su & Lai (2020) studied another mechanism in which a shrinking protoplanetary disk could cause sweeping spin–orbit resonances to influence the obliquity distribution of newly formed protoplanets. Another study applied this concept to Uranus's circumplanetary disk to attempt to infer its dynamical history (Rogoszinski & Hamilton 2020).

Alternatively, a planet can be driven into or out of a spin–orbit resonance by processes that alter its spin-precession rate (Equation (1)). Tidal spin torques are one mechanism that can dynamically shape the fate of a planet's spin state over time. The tidal friction generated from a system's star pulling on the equatorial bulges of a planet will slowly drive it toward tidal equilibrium, in which its obliquity will typically increase followed by a gradual decline to small values, while its spin rate undergoes rotation braking until its spin period becomes synchronous with its orbital period (or pseudosynchronous for eccentric orbits). This journey amounts to a constant decline in a planet's precession constant and, consequently, its spin-precession rate over time, allowing for opportunities to encounter spin–orbit resonances as time goes on. Tidal spin torques can drive a planet's spin axis to traverse parameter space and undergo scenarios such as spin–orbit resonance crossings, captures, or even traverses of chaotic resonance regions. Since the magnitude of the tidal force decays roughly as r⁻⁶ with distance, this process is likely especially important for planets in close proximity to their star. The Earth will likely undergo chaotic obliquity variations itself in a few billion years after entering the chaotic zone thanks to its decaying rotation rate (Neron de Surgy & Laskar 1997).

We introduce SMERCURY-T (Kreyche et al. 2021), a joining of the existing code SMERCURY (Lissauer et al. 2012) with elements from Mercury-T (Bolmont et al. 2016), so that we can investigate these dynamic tidal spin phenomena and build a more complete picture of planetary spin dynamics. While SMERCURY utilizes its implementation of the spin Hamiltonian from Touma & Wisdom (1994) to accurately compute the orbital and rotational dynamics of a planetary system due to N-body interactions, it is not equipped to compute tidal evolution. On the other hand, Mercury-T is capable of computing the tidal evolution of a system but limits the study of planetary spin evolution by neglecting planet–planet cross-terms in its spin calculations. In addition, Mercury-T lacks a scheme to efficiently explore parameter space like the "ghost planet" scheme featured in the SMERCURY-T code (discussed in Section 2.1). Similar open-source codes, such as Posidonius (Blanco-Cuaresma & Bolmont 2017) and VPLanet (Barnes et al. 2020), are also guilty of these charges. We bridge this gap by branching SMERCURY-T from SMERCURY to include modifications that allow for the complete simultaneous computation of the orbital and spin evolution of a planet while under the effects of stellar tidal spin torques.

In this paper, we begin with Section 2, in which we describe the SMERCURY-T code and discuss our implementation of its new modules. Following that, in Section 3, we perform several experiments that demonstrate the phenomena of spin–orbit resonance crossings, captures, and chaotic crossings that express the versatility of SMERCURY-T as a tool. We conclude in Section 4.

2. Methods

2.1. The Base Code

We present our new code, SMERCURY-T, as a modified version of its predecessor, SMERCURY (Lissauer et al. 2012; Barnes et al. 2016; Kreyche et al. 2020; Quarles et al. 2020). SMERCURY is a mixed-variable symplectic N-body integrator that is itself a modified spin-tracking version of the original MERCURY package (Chambers 1999). Just like the code that SMERCURY-T derives from, it is capable of computing the orbital evolution of multiplanet systems while tracking the spin state of one planet and its "ghosts." The ghost planet scheme entails the simultaneous simulation of massless copies of the spin-tracked planet, in which the user assigns a unique initial spin state to each ghost planet that will evolve while its orbital evolution mimics that of the original planet. This scheme makes SMERCURY-T advantageous over similar codes by allowing for efficient exploration of parameter space; however, its present implementation restricts the user to only being able to simulate the spin evolution of one system's planet at a time. SMERCURY-T is publicly available (Kreyche et al. 2021).

SMERCURY-T treats the spin computation by considering the spin-tracked planet to be an axisymmetric rigid body. Gravitational torques from the planet's star and neighboring planets exerted on the spin-tracked planet's equatorial bulge modify its obliquity. When preparing a simulation, we express the planet's spin rate by assigning it a value for its zonal gravity coefficient J₂ according to the Darwin–Radau relation (Hubbard 1984; Murray & Dermott 1999). This takes the form

$\begin{eqnarray}&&{J}_{2}=\displaystyle \frac{{\omega }_{p}^{2}{R}^{3}}{3{{Gm}}_{p}}\left(\displaystyle \frac{5}{D}-1\right),\end{eqnarray} \tag{ 3 }$

with ω_p, R, and m_p as the planet's spin rate, mass, and equatorial radius, respectively. Here G is the universal gravitational constant, and the placeholder D is short for

$\begin{eqnarray}&&D=\displaystyle \frac{25}{4}{\left(\displaystyle \frac{3}{2}\overline{C}-1\right)}^{2}+1,\end{eqnarray} \tag{ 4 }$

where $\overline{C}$ is the planet's normalized polar moment of inertia. We calculate R according to

$\begin{eqnarray}&&R={\left[\displaystyle \frac{{{Gm}}_{p}D}{{\omega }_{p}^{2}}\left(\displaystyle \frac{2-\sqrt{4-\tfrac{30{\omega }_{p}^{2}}{{GD}\pi \rho }}}{10}\right)\right]}^{\displaystyle \frac{1}{3}},\end{eqnarray} \tag{ 5 }$

with ρ as the average density of the planet. We refer to Lissauer et al. (2012) for a more thorough description of this technique.

2.2. Code Modifications

In this section, we describe the modifications that we include within the SMERCURY-T package based on the alterations that we made to the base SMERCURY code. Here we include the details and implementation of the tidal spin torque module and then discuss the results of a test experiment to verify its accuracy. We also discuss our additional inclusion of a general relativity module.

2.2.1. Tidal Spin Torque Module

The most significant inclusion in the SMERCURY-T code is the tidal spin torque module. We follow the methods of Bolmont et al. (2015) and adopt the constant time lag tidal formulation expressed by Leconte et al. (2010) and others (Mignard 1979; Hut 1981; Eggleton et al. 1998), which is valid for arbitrary values of obliquity, spin, and orbital eccentricity. Contrary to the rigid-body considerations that we employ for the N-body torque computations, this model assumes that the body is made of a weakly viscous fluid (Alexander 1973). In the future, SMERCURY-T's modular design will make it simple to implement and explore alternative tidal models. Our approach differs from the aforementioned works in that we neglect orbital tides and consider only the effects of tidal spin torque on the spin-tracked planet from its star. This choice saves valuable computation time while still allowing for sufficient study of the behavior of planetary spin evolution in most planetary systems, except the most compact ones. Readers can refer to Table 2 of Mardling & Lin (2004) or perform quick calculations with handy equations such as those found in Rodríguez & Ferraz-Mello (2010) to estimate the importance of orbital tides in their desired system. We justify our study of the systems that we consider in Section 3 by estimating the semimajor axis damping timescale ( $a/\dot{a}$ ) and eccentricity damping timescale ( $e/\dot{e}$ ) to be practically infinite for the Earth-mass planet under study compared to our integration time of 4 Gyr. Another consequence of our methodology is that it is technically in violation of the law of conservation of angular momentum, since the spin-tracked planet's spin angular momentum decays while the system gets nothing back in return. The choice to neglect the effects of this feedback is acceptable for the study of planetary systems, since the total angular momentum of the system will always vastly dwarf the spin-tracked planet's spin angular momentum.

The task to modify the spin routine of SMERCURY-T to include tidal spin torques is a relatively straightforward one. Since SMERCURY-T already integrates the summation of all of the torques felt by the spin-tracked planet's rotational bulge at each time step, all that remains is calculating the additional tidal torque within the routine. From Bolmont et al. (2015), the nonaveraged tidal torque felt by a planet from its star goes as

$\begin{eqnarray}&&{{\boldsymbol{N}}}_{{\boldsymbol{T}}}=3G\displaystyle \frac{{{m}_{s}}^{2}{R}^{5}}{{r}^{7}}{k}_{2}\tau (r{{\boldsymbol{\omega }}}_{p}-({\boldsymbol{r}}\cdot {{\boldsymbol{\omega }}}_{p}){{\boldsymbol{e}}}_{r}-{{\boldsymbol{e}}}_{r}\times {\boldsymbol{v}}),\end{eqnarray} \tag{ 6 }$

where r is the radial vector between the planet and its star, ω _p is the planet's spin vector, e _r is the radial unit vector, v is the planet's velocity vector, k₂ is its potential Love number of degree 2, τ is the constant time lag, m_s is the mass of the star, and G is the universal gravitational constant. This calculation treats the star as a point mass, while the planet is made up of a reduced central point mass with two point-mass bulges. Figure 1 from Bolmont et al. (2015) is a useful diagram that demonstrates some of these elements. Converting this result to heliocentric coordinates allows for a final expression that conforms with the requirements of SMERCURY-T's spin algorithm. Taking the planet's spin angular momentum as S , the adjusted tidal torque is

$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}{\boldsymbol{S}}=-\displaystyle \frac{{m}_{s}}{{m}_{s}+{m}_{p}}{{\boldsymbol{N}}}_{T}.\end{eqnarray} \tag{ 7 }$

As the spin-tracked planet's spin rate damps over time due to the inclusion of the tidal torque, it is necessary to update the planet's J₂ and equatorial radius at each time step to represent the reduction of its equatorial bulge. This is important not only due to the dependence of the strength of the tidal torque on R but even more so due to the sensitive nature of the N-body spin dynamics computations to the value of J₂ as the planet's spin state evolves. We develop a scheme to update these parameters and provide the details in the Appendix. This scheme involves instructing SMERCURY-T to use the value of the planet's spin angular momentum to solve a cubic equation to ultimately determine new values for R and J₂ using Equations (A6) and (3).

One possible concern with the inclusion of the tidal spin torque module is that there may be scenarios in which the magnitude of the computed tidal torque is so small that it underflows machine precision, causing the planet's tidal evolution to cease entirely. We solve this issue by implementing a "tidal tolerance" input parameter that the user can specify prior to running a simulation. This parameter is a fractional value that determines a certain threshold after being multiplied by the total spin angular momentum at the current time step. If the tidal torque components exceed this threshold, then they are integrated and added to the spin vector. However, if they fall below the threshold, then SMERCURY-T saves their values and adds them to a running total to be checked the next time around. A typical value for the tidal tolerance parameter that we use in this paper is 10⁻¹².

2.3. Tidal Spin Torque Module Verification

We verify the implementation of the SMERCURY-T tidal spin torque module by simulating and comparing the results of a test system's tidal evolution to that of a secular code based on the framework of Leconte et al. (2010). The secular code is capable of solving the averaged equations of one planet's orbital and rotational tidal evolution. However, since SMERCURY-T only accounts for the effects of tidal spin torques, for a fair comparison, we hold the orbital elements constant and evolve the secular code solely with the spin equations. We implement the same scheme from Section 2.2.1, detailed in the Appendix, into the secular code to similarly update the planet's equatorial radius at each time step.

The test system that we use here is a one-planet system that consists of the Sun and an Earth-mass planet with a semimajor axis of 0.2 au and an orbital eccentricity of 0.3. This choice ensures a short integration time and demonstrates the capability of SMERCURY-T to handle large eccentricities. For both simulations, we assign this planet a polar moment of inertia, average density, initial spin period, initial equatorial radius, potential Love number of degree 2, constant time lag, and initial obliquity according to Table 1. We integrate the SMERCURY-T simulation with a time step of 1.6 days (∼5% of the innermost planet's orbital period, our usual rule of thumb) and sample at an interval of 1000 yr. The time step for the averaged equations of the secular code does not need to be as small, so we proceed with 1000 yr increments. We run both simulations for 10 Myr, which is more than enough time for the planet to reach tidal equilibrium.

Table 1. Planet Physical Parameters

$\overline{C}$	ρ (g cm⁻³)	P₀ (hr)	R₀ (km)	k₂	τ (s)	₀ (deg)
0.3296108	5.513 6	24	6378	0.305	698	45

Note. Physical values for the test planet that we describe in Section 2.3. Here we list the planet's moment of inertia ( $\overline{C}$ ), average density (ρ), initial spin period (P₀), initial equatorial radius (R₀), potential Love number of degree 2 (k₂), constant time lag (τ), and initial obliquity (₀).

Download table as: ASCII Typeset image

Looking to Figure 1, we show a comparison between the results of the two codes for the spin evolution of the Earth-mass planet due to the effects of tidal spin torques. Examining the evolution of the planet's obliquity, spin period, equatorial radius, J₂ (initially calculated with Equation (3)), and spin angular momentum, we show that SMERCURY-T is in good agreement with the secular code. Although there is a slight discrepancy between the extent of the planet's obliquity values near the peaks of the plot in Figure 1's top left panel, each parameter ultimately tracks and converges to the same end value with little difference between the two codes. In summary, we see that tidal spin torques cause the planet to initially experience an increase in its obliquity until its spin rate decays such that the altitudinal component of its spin angular momentum overtakes the azimuthal component, as shown in the top and bottom right panels of Figure 1. From there, its spin period continues to increase while its obliquity declines until it reaches tidal equilibrium when its obliquity is near zero and its spin period is that of its pseudosynchronous period (∼21 days here). The planet's equatorial radius shrinks by about 7 km, while its J₂ reduces considerably.

2.3.1. General Relativity Module

We also incorporate the general relativistic (GR) force due to the post-Newtonian potential into SMERCURY-T as an optional module following Section 2.2 of Bolmont et al. (2016). The GR force not only plays an important role with regard to the stability of a planetary system, but we reason that it can also have an indirect influence on planetary obliquity evolution. This influence stems from the GR effect that causes the precession of periapsis of planetary orbits to speed up. Bearing this in mind, GR can indirectly influence planetary obliquity by shifting the potentially resonant orbital eigenfrequencies. While these effects are important to a degree in our own solar system, they are likely especially important for compact exoplanet systems whose planets reside deep within their star's gravitational well.

In the spirit of one famous confirmation of Einstein's general theory of relativity (Einstein 1916), we verify the implementation of our GR module by checking the contribution of GR to the precession of periapsis of our planetary neighbor Mercury. We compare the results of our SMERCURY-T simulation to that of the measured rate from the MESSENGER spacecraft reported by Park et al. (2017). We obtain an average value of 42 farcs 9790 per Julian century over the course of our 1 Myr simulation, in which we used a system consisting only of the Sun and Mercury to isolate the effect of GR. Specifically, this value is the difference between Mercury's precession rate when the GR module was enabled and another simulation when it was not. Here we set the integration time step to be 5% of Mercury's orbital period and sampled at 100 yr intervals. This value is in good agreement and lies within the error bars of the Park et al. (2017) result.

3. Experiments

In this section, we apply the code modifications discussed in Section 2.2 to demonstrate the usefulness and versatility of SMERCURY-T as a tool to study planetary spin dynamics. Here we showcase examples of three different scenarios driven by the effects of tidal spin torques: a spin–orbit resonance crossing, a spin–orbit resonance capture, and a chaotic resonance crossing. We enable the general relativity and tidal spin torque modules ("tidal tolerance" set to 10⁻¹²) for these experiments, with the exception of the initial frequency analysis runs detailed in Sections 3.1.1 and 3.3.1.

3.1. Spin–Orbit Resonance Crossing

3.1.1. Initial Conditions and Preanalysis

We design an experiment to exemplify the scenario in which tidal spin torques drive a planet to undergo a spin–orbit resonance crossing as its increasing obliquity and spin period cause its precession constant to decline over time. We choose to use a simple toy system consisting of the Sun, an Earth-mass planet initially placed at 0.5 au, and a Jupiter-mass planet initially placed at 2.5 au. This choice allows us to study this process appropriately while minimizing the necessary integration time thanks to the stronger presence of tides closer to the Sun. We display the masses and orbital elements of these bodies in Table 2 listed under system E-J, which shows that the Earth- and Jupiter-mass planet share an initial mutual inclination of 1°.

Table 2. System Orbital Parameters

System	Planet	m_p (M_⊙)	a (au)	e	I (deg)	ω (deg)	Ω (deg)	M (deg)
E-J	Earth	3.003 × 10⁻⁶	0.5	0	1	0	0	0
	Jupiter	9.542 44 × 10⁻⁴	2.5	0	0	0	0	0

E-E-J	Inner Earth	3.003 × 10⁻⁶	0.5	0	5	0	0	0
	Outer Earth	3.003 × 10⁻⁶	0.7	0	2.5	0	0	0
	Jupiter	9.542 44 × 10⁻⁴	3.5	0	0	0	0	0

Note. Initial values for the planetary mass (m_p), semimajor axis (a), eccentricity (e), inclination (I), argument of periapsis (ω), longitude of ascending node (Ω), and mean anomaly (M) for the Earth–Jupiter (E-J) and Earth–Earth–Jupiter (E-E-J) systems that we test in this paper. We take the mass of the Sun to be 1.98911 × 10³⁰ kg.

Download table as: ASCII Typeset image

We perform an initial simulation of the E-J system over the course of 10 Myr, sampling on an interval of 1000 yr with a numerical time step of 6.46 days (∼5% of the Earth-mass planet's orbital period). For this part, we are only interested in the system's orbital evolution, so we leave the GR module enabled but disable the tidal spin torque module in the interest of shaving off some computation time. Taking these results, we perform a frequency analysis routine by applying a frequency-modified Fourier transform (Šidlichovský & Nesvorný 1996) to the Earth-mass planet's inclination and eccentricity vectors ([ $\sin \tfrac{I}{2}\cos {\rm{\Omega }}$ , $\sin \tfrac{I}{2}\sin {\rm{\Omega }}]$ and $[e\cos \varpi$ , $e\sin \varpi ]$ respectively). This process extracts the frequencies, amplitudes, and initial phases of the precessional modes that could enter resonance with the Earth-mass planet. First-order resonances appear directly in the inclination series. We display these values in Table 3 listed under system E-J. This information allows us to compute both the location of a resonance center and its width enclosed by the separatrix as a function of a planet's precession constant and obliquity, following the procedure of Saillenfest et al. (2019). The simplicity of the E-J system results in the single dominant ν₁ mode in the inclination series at ≈−22 farcs 69 yr⁻¹. This frequency sets the target resonance for the Earth-mass planet in the experiments we describe in Sections 3.1.2 and 3.2.2.

Table 3. Secular Orbital Frequencies

			Inclination Vector Series			Eccentricity Vector Series
System	j	ν_j (arcsec yr^–1)	S_j × 10⁸	${\phi }_{j}^{(0)}$ (deg)	μ_j (arcsec yr^–1)	E_j × 10⁸	${\theta }_{j}^{(0)}$ (deg)
E-J	1	−22.69	871,406	359.94	22.89	804.86	179.80
	2	0.00	1221	0.00	0.03	431	183.31
	3	−22.47	183	15.31	−9.11	15	41.66
	4	−22.90	170	149.56	−6.23	14	346.47
	5	−22.30	102	6.21	−7.45	13	235.11

E-E-J	1	−10.67	3,736,801	359.99	30.74	704	180.09
	2	−32.25	617,977	359.84	0.03	617	180.22
	3	−3.97	8169	0.09	12.59	53	2.33
	4	−53.83	2096	179.43	9.15	14	191.46
	5	10.91	786	358.52	−98.64	8	77.37

Note. The results from our frequency-modified Fourier transform analysis of the inclination ([ $\sin \tfrac{I}{2}\cos {\rm{\Omega }}$ , $\sin \tfrac{I}{2}\sin {\rm{\Omega }}$ ]) and eccentricity ([ $e\cos \varpi$ , $e\sin \varpi$ ]) vectors of the Earth-mass planet in the E-J system and the inner Earth-mass planet in the E-E-J system, where I is the orbital inclination, Ω is the longitude of the ascending node, e is the eccentricity, and ϖ is the longitude of periapsis. We show the top five values of the frequencies (ν_j), amplitudes (S_j), and phases ( ${\phi }_{j}^{(0)}$ ) of the inclination vector, as well as the top five values of the frequencies (μ_j), amplitudes (E_j), and phases ( ${\theta }_{j}^{(0)}$ ) of the eccentricity vector.

Download table as: ASCII Typeset image

Targeting the ν₁ mode, we initialize the Earth-mass planet so that it will begin its tidal evolution just outside of the resonance and cross paths with the resonance center some time later. Therefore, we assign it an initial obliquity of 50° and a precession constant of 60'' yr⁻¹ (using Equation (2) with a spin period of ≈54.7 hr and J₂ = 2.03779 × 10⁻⁴). These planetary characteristics are not unrealistic and could be the aftermath of various primordial scenarios, such as giant impacts. We set its initial precession angle to be 0°, while we use the same values for its average density, normalized polar moment of inertia, potential Love number of degree 2, and constant time lag as found in Table 1. We run this simulation for 4 Gyr and sample every 10,000 yr with the same 6.46 day time step.

3.1.2. Results and Discussion

We explore the results of the spin–orbit resonance crossing experiment by looking to Figure 2. We see that the Earth-mass planet begins with an obliquity of 50° and experiences benign obliquity variations of a few degrees that are primarily due to the precessional motion of its orbit. The planet's precession constant declines, while tidal spin torques cause its obliquity to gradually increase as it approaches the resonance. Upon entering the resonance region, the planet's obliquity swings wildly in the range of ∼50°–60°. Then, about 98 Myr into the simulation, the planet's obliquity suddenly shifts 10° as it crosses the resonance center and swings into a new range of ∼60°–70°. The planet then proceeds to exit the resonance and eventually converges toward tidal equilibrium with the behavior of Cassini state 1, as its obliquity draws near zero and its spin period approaches the synchronous period of ∼129 days. The nature of these short-term obliquity variations as the planet crossed through this resonance during its tidal evolution would likely be harmful toward its prospects for climatic stability and habitability. However, a true assessment of these prospects would require a more complete study that would explore the application of a climate model. We include an accessory video that shows a time lapse of the crossing event in Figure 3.

Figure 3. The real-time duration of this video is 30 s, during which we show a time lapse of the results from the upper panel of Figure 2; refer to the caption of Figure 2 for a description of the chart elements. Here we represent the Earth-mass planet's parameter position with a black dot and display an associated time stamp for each frame in the lower right corner. From the range 93–103 Myr into the simulation, tidal spin torques from the star cause the Earth-mass planet's precession constant to decline, while it experiences strong obliquity variations as it traverses the resonance region. The planet's obliquity receives a prompt kick upon crossing the resonance center from the left to the right side and then begins to move away from the resonance center afterward.

(An animation of this figure is available.)

Download figure:

Video Standard image High-resolution image

Further examination of the resonance crossing event leads us to reason that our choice of the planet's initial precession angle ensured that it crossed through the resonance's hyperbolic point (Saillenfest et al. 2021a). In other words, this behavior resulted due to the chance timing and positioning of the planet's spin axis as its spin precession wound down while the planet crossed the resonance center. Looking to Figure 4, we examine the evolution of the angle σ = ψ − Ω centered around the time of the resonance center crossing event at ∼98 Myr, where ψ is its spin-precession angle and Ω is the planet's longitude of ascending node. We treat this angle as a close approximation of the resonant angle, in which it is diagnostic of the nature of a resonance event. This angle tends to librate about 180° when a planet is captured in resonance in the vicinity of Cassini state 2; otherwise, it usually circulates (Deitrick et al. 2018a; Quarles et al. 2020). Here we show that the resonant angle circulates while the planet approaches and moves away from the resonance center while an inflection point appears at the time of the crossing. This behavior reinforces that the planet was never captured into the resonance at any point in time and indeed crossed through its hyperbolic point.

**Figure 4.** We show the evolution of the Earth-mass planet's (from the E-J system described in Table 2) resonant angle over time centered about the time that the planet crossed the ν₁ resonance center. In this case, the planet's initial precession angle, ψ, was 0°. For the resonant angle, we use σ = ψ − Ω, with Ω as the planet's longitude of ascending node. This angle is diagnostic of the presence of a planet in a spin–orbit resonance. Here the angle circulates with a brief inflection point at the time of the crossing event (∼98 Myr), which indicates that the planet was not captured into the resonance and instead crossed through its hyperbolic point (Saillenfest et al. 2021a).
Download figure:
Standard image High-resolution image

3.2. Spin–Orbit Resonance Capture

3.2.1. Initial Conditions and Preanalysis

Following the spin–orbit resonance crossing experiment described in Section 3.1.1, we now consider an alternative fate for the Earth-mass planet in which it could instead become captured into the spin–orbit resonance. Here we use the exact same initial conditions as the former experiment, save one: the initial precession angle. Since the outcome is probabilistic based on the value of the precession angle upon its arrival at the resonance (Henrard 1982; Saillenfest et al. 2021a), we consider a slew of initial precession angles but handpick the case of the Earth-mass planet with an initial precession angle of 180°. We stress here that although this serves our example, this is an atypical pathway toward resonance capture, since the planet crosses the resonance "from above" (α/ν₁ proceeds toward unity from greater to smaller values; Hamilton & Ward 2004; Ward & Hamilton 2004). We simulate this system for 4 Gyr with the same sample rate and time step as the experiment described in Section 3.1.1.

3.2.2. Results and Discussion

We show the results of the spin–orbit resonance capture experiment by looking to Figure 5. While the Earth-mass planet begins on a track very similar to the crossing experiment described in Section 3.1.2, its spin evolution takes a very different path once it enters the resonance. About 75 Myr into the simulation, the planet crosses the resonance center and is subsequently captured into the ν₁ resonance. From here, its obliquity experiences large variations of order ∼20° as it is bound to the follow the track of the resonance's center as its precession constant declines. Eventually, at around the 350 Myr mark, the separatrix of the resonance disappears, and the planet exits the resonance. Since the planet's ride within the resonance led to an inefficient decay of its spin rate, the planet's obliquity then trends to greater values until it peaks at about 50° and then reverses while the planet heads toward tidal equilibrium. The path that the planet's obliquity ultimately took from start to finish spanned a wide range of values with frequently large obliquity variations over short timescales. Similar to the resonance crossing experiment from Section 3.1.2, this behavior would likely be detrimental toward the planet's prospects for habitability. Look to Figure 6 for a video that shows a time lapse of the initial resonance capture event.

Figure 6. The real-time duration of this video is 30 s, during which we show a time lapse of the results from the top panel of Figure 5; refer to the caption of Figure 2 for a description of the chart elements. Similar to Figure 3, we represent the Earth-mass planet's parameter position with a black dot and display an associated time stamp for each frame in the lower right corner. From the range 73–83 Myr into the simulation, tidal spin torques from the star cause the Earth-mass planet's precession constant to decline, while it experiences strong obliquity variations as it traverses the resonance region. The planet becomes captured upon encountering the resonance center, after which it maintains large obliquity variations about the resonance center while its precession constant continues to gradually decline.

(An animation of this figure is available.)

Download figure:

Video Standard image High-resolution image

Why did this experiment lead to the Earth-mass planet getting captured into the resonance while the other did not? Recalling that we set this planet to have an initial precession angle of 180° rather than the former experiment's 0°, we show the temporal evolution of its resonant angle centered around the event in Figure 7. Here we show that prior to the planet's encounter with the resonance center its resonant angle circulates. However, following the event, the resonant angle then librates about 180°, albeit very widely, which is confirmation of the planet's capture into the resonance. This time, our choice of initial precession angle ensured that the conditions were ripe for capture into the resonance at the time of crossing.

**Figure 7.** Similar to Figure 4, here we show an alternative evolution of the Earth-mass planet's (from the E-J system described in Table 2) resonant angle over time centered about the time that the planet crossed the ν₁ resonance center. In this case, the planet's initial precession angle was 180°. We see that this time, the resonant angle transitions from circulating to librating about 180° upon the encounter, which is a strong indicator that it was captured into the resonance.
Download figure:
Standard image High-resolution image

3.3. Chaotic Crossing

3.3.1. Initial Conditions and Preanalysis

For our final experiment, we design a scenario that exemplifies a chaotic spin–orbit resonance crossing. The resonance overlap criterion specifies that chaotic regions of parameter space exist in the case where two first-order resonances overlap. Therefore, we now consider a system with an additional planet to introduce an additional forced inclination frequency (Saillenfest et al. 2019). We present this system in Table 2 listed under system E-E-J, which shows that our new system consists of the Sun, an Earth-mass planet at 0.5 au with an inclination of 5°, an Earth-mass planet at 0.7 au with an inclination of 2 fdg 5, and a Jupiter-mass planet at 3.5 au whose orbit sits flat with an inclination of 0°.

We conduct the same frequency analysis as described in Section 3.1.1 to extract the information on the relevant modes of this system. We display this information in Table 3 listed under system E-E-J. This time, our inclusion of the additional Earth-mass planet raised the "temperature" of the system, where the inner Earth shows two prominent driving frequencies in its inclination series. The first mode, ν₁, has a frequency of ≈−10 farcs 67 yr⁻¹ and is the primary driver of the inner Earth-mass planet's orbital precession. The second mode, ν₂, has a frequency of ≈−32 farcs 25 yr⁻¹ with a considerably smaller amplitude. These two frequencies have resonant widths that overlap to form a chaotic region.

We set the inner Earth's obliquity to 60° and its precession constant to 110'' yr⁻¹ (using Equation (2) with a spin period of ≈29.88 hr and J₂ = 6.83897 × 10⁻⁴) so that it will begin its tidal evolution just outside of the ν₂ resonance region and tidally evolve to cross into the chaotic region. The inner Earth has the same average density, normalized polar moment of inertia, potential Love number of degree 2, and constant time lag found in Table 1. We simulate this system over the course of 4 Gyr by sampling on 10,000 yr intervals with a 6.46 day time step.

3.3.2. Results and Discussion

We show the results of the chaotic spin–orbit resonance crossing experiment in Figure 8. These results share similarities with Figures 4(a) and (b) of Neron de Surgy & Laskar (1997), who studied possible chaotic futures for the real Earth. Here the inner Earth begins with an obliquity of 60° and experiences variations of order 5°–10° that are on par with benign variations resulting primarily from the precessional motion of its orbit. As tidal spin torques dampen the planet's precession constant and push its obliquity toward higher values, the planet crosses into the ν₂ resonance after about 150 Myr. The resonance excites larger obliquity variations as the planet first crosses the center of the ν₂ resonance and then into the chaotic overlap region of the ν₂ and ν₁ resonances. Chaos ensues, in which the planet's obliquity violently swings as much as ∼60° in an unpredictable fashion within the range 30°–90°. This behavior spans 200–600 Myr into the simulation until the planet finally transitions to crossing through just the ν₁ resonance and eventually exits to proceed toward tidal equilibrium, where its obliquity trends toward zero and its spin period approaches the synchronous period of ∼129 days. Although this planetary system is fictional, its example of chaotic spin evolution implicates the probable catastrophic effects that this behavior could have on the planet's well-being. Refer to Figure 9 for a video that shows a time lapse of a period of this chaotic behavior.

Figure 9. The real-time duration of this video is 30 s, during which we show a time lapse of the results from the top panel of Figure 8; refer to the caption of Figure 2 for a description of the chart elements. Similar to Figure 3, we represent the Earth-mass planet's parameter position with a black dot and display an associated time stamp for each frame in the lower right corner. From the range 210–220 Myr into the simulation, tidal spin torques from the star cause the Earth-mass planet to cross through the chaotic overlap region as its precession constant declines, causing its obliquity to undergo large and chaotic swings.

(An animation of this figure is available.)

Download figure:

Video Standard image High-resolution image

4. Conclusion

In this paper, we presented our new code, SMERCURY-T, a numerical N-body integrator that is capable of computing the orbital and rotational evolution of a planet under the influence of tidal spin torques. Additionally, we included a module that adds the effects of the GR force to planetary orbital evolution. We verified the proper implementation of these inclusions and showed the potential of SMERCURY-T as a powerful and versatile tool that is capable of tackling a variety of problems.

We used SMERCURY-T to perform a series of experiments to investigate different phenomena that a planet could experience due to the effects of tidal spin torques. These processes are likely important to the long-term dynamics of many planets. The spin–orbit resonance crossing and capture experiments demonstrated the possibility that a planet can take vastly different pathways that are determined purely by its initial conditions. While the resonance crossing experiment showed an example in which tidal spin torques drew a planet to encounter a spin–orbit resonance followed by an abrupt kick to its obliquity, the resonance capture experiment instead resulted in the planet engaging in a stable libration within the resonance while it was transported toward lower obliquity for a period of time. We also showed the results of the crossing of a chaotic region formed by the overlap of two spin–orbit resonances. The large chaotic swings of the planet's obliquity are familiar and consistent with the dynamics of some of the bodies in our own solar system.

Although the planetary systems that we considered in this paper were fictional, our general study of the nature of tidal obliquity evolution demonstrates possible scenarios that may have implications for the past, present, and future dynamics of planetary obliquity evolution. Knowing that the climatic stability and, consequently, habitability of the planet is strongly related to the nature of its obliquity evolution, understanding the dynamic role that tidal spin torques can play is critical knowledge. Large and rapid obliquity variations like those that we observed in our experiments would likely be detrimental for a planet's climatic stability and therefore its prospects for habitability. Although we did not investigate this claim in this paper, an interesting avenue for future work could involve the application of a climate model to explore the potentially destructive consequences of these evolutionary pathways.

For future work, SMERCURY-T can serve as a tool to further explore tidal obliquity evolution and seek answers to questions about the planets in our solar system and beyond. The effects of tidal spin torques are especially relevant to planets that are in close proximity to their star, such our own system's inner terrestrial planets or those in compact exoplanet systems such as TRAPPIST-1. Future characterization of exoplanet obliquities will allow for the opportunity to perform more complete assessments of their dynamical nature.

We thank the anonymous referees for their insights and suggestions that helped to improve this manuscript. This work was supported by the NASA Habitable Worlds program, grant No. 80NSSC19K0312.

Appendix: Scheme to Update Planetary Radius and J₂

In this appendix, we discuss the methodology that SMERCURY-T uses to update the planetary radius and zonal gravity coefficient of a planet throughout a simulation with tidal spin torques enabled. The inclusion of the tidal spin torque module that we discuss in Section 2.2.1 calls for the need to update the values of the spin-tracked planet's equatorial radius, R, and zonal gravity coefficient, J₂, at each time step within a SMERCURY-T simulation (or at an interval set by the "tidal tolerance" parameter). This is necessary because the effects of the tidal torque on the planet's equatorial bulge cause its spin rate, and therefore its equatorial radius and J₂, to decay over time. These parameters are critical for SMERCURY-T to compute the magnitude of the tidal torque following Equation (6), as well as to populate the planet's inertia tensor for the N-body torque calculations following Lissauer et al. (2012).

Sticking with the rigid-body considerations described in Section 2.1, we begin an expression for the planet's equatorial radius (like that of Equation (5)) as

$\begin{eqnarray}&&{R}^{3}=\displaystyle \frac{{{Gm}}_{p}D}{{\omega }_{p}^{2}}\left(\displaystyle \frac{2-\sqrt{4-\tfrac{30{\omega }_{p}^{2}}{{GD}\pi \rho }}}{10}\right),\end{eqnarray} \tag{ A1 }$

where m_p, ω_p, and ρ are the planet's mass, spin rate, and average density, respectively. Here G is the universal gravitational constant, and

$\begin{eqnarray}&&D=\displaystyle \frac{25}{4}{\left(\displaystyle \frac{3}{2}\overline{C}-1\right)}^{2}+1,\end{eqnarray} \tag{ A2 }$

where $\overline{C}$ is the planet's normalized polar moment of inertia. Now we can shuffle the equation S = I_p ω_p for the planet's rotational angular momentum to get

$\begin{eqnarray}&&{\omega }_{p}=\displaystyle \frac{S}{{I}_{p}}=\displaystyle \frac{S}{\overline{C}{m}_{p}{R}^{2}}.\end{eqnarray} \tag{ A3 }$

Here S is the magnitude of the planet's spin angular momentum, and I_p is its moment of inertia. We plug Equation (A3) into Equation (A1) to work toward a solution for R. This gives

$\begin{eqnarray}&&\begin{array}{l}{R}^{3}=\displaystyle \frac{G{{m}_{p}}^{3}D{\overline{C}}^{2}{R}^{4}}{{S}^{2}}\\ \quad \left(\displaystyle \frac{2-\sqrt{4-\tfrac{30{S}^{2}}{G{{m}_{p}}^{2}D{\overline{C}}^{2}\pi \rho {R}^{4}}}}{10}\right).\end{array}\end{eqnarray} \tag{ A4 }$

Shifting some things around, squaring both sides, and isolating the R terms to one side yields

$\begin{eqnarray}&&\begin{array}{l}-\displaystyle \frac{40{S}^{2}}{G{{m}_{p}}^{3}D{\overline{C}}^{2}}{R}^{3}+\displaystyle \frac{100{S}^{4}}{{G}^{2}{{m}_{p}}^{6}{D}^{2}{\overline{C}}^{4}}{R}^{2}\\ \,+\,\displaystyle \frac{30{S}^{2}}{G{{m}_{p}}^{2}D{\overline{C}}^{2}\pi \rho }=0.\end{array}\end{eqnarray} \tag{ A5 }$

This reduces to

$\begin{eqnarray}&&-\displaystyle \frac{4}{{m}_{p}}{R}^{3}+\displaystyle \frac{10{S}^{2}}{G{{m}_{p}}^{4}D{\overline{C}}^{2}}{R}^{2}+\displaystyle \frac{3}{\pi \rho }=0.\end{eqnarray} \tag{ A6 }$

This equation takes the form ax³ + bx² + cx + d = 0, for which we solve for R with the use of Cardano's general cubic formula (Cardano 1993). We obtain the corresponding value of ω_p by plugging the solution for R back into Equation (A3). Finally, we calculate the planet's J₂ with Equation (3).

Exploring Tidal Obliquity Variations with SMERCURY-T

Article metrics

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction