Evolution of Planetary Obliquity: The Eccentric Kozai–Lidov Mechanism Coupled with Tide

Planetary obliquity plays a significant role in determining the physical properties of planetary surfaces and climate. As direct detection is constrained due to the present observation accuracy, kinetic theories are helpful for predicting the evolution of planetary obliquity. Here the coupling effect between the eccentric Kozai–Lidov effect and the equilibrium tide is extensively investigated; the planetary obliquity is observed to follow two kinds of secular evolution paths, based on the conservation of total angular momentum. The equilibrium timescale of the planetary obliquity t eq varies along with r t , which is defined as the initial timescale ratio of the tidal dissipation and secular perturbation. We numerically derive the linear relationship between t eq and r t with the maximum likelihood method. The spin-axis orientation of S-type terrestrials orbiting M-dwarfs reverses over 90° when r t > 100, then enters the quasi-equilibrium state between 40° and 60°, while the maximum obliquity can reach 130° when r t > 104. Numerical simulations show that the maximum obliquity increases with the semimajor axis ratio a 1/a 2, but is not so sensitive to the eccentricity e 2. The likelihood of an obliquity flip for S-type terrestrials in general systems with a 2 < 45 au is closely related to m 1. The observed potentially oblique S-type planets HD 42936 b, GJ 86 Ab, and τ Boo Ab are explored and found to have a great possibility of rotating head-down over the secular evolution of spin.


INTRODUCTION
The planetary obliquity is simply defined as the angle between the planetary spin and normal orientation of the orbital plane.In our solar system, Uranus and Venus are well-known as retrograde planets with extremely high planetary obliquity, where the obliquity of Uranus is 97 • , whereas that of Venus is approximate to 178 • (Goldreich & Peale 1970;Dobrovolskis 1980).The primordial planetary obliquity provides a key clue to the understanding of evolution of our solar system, for example, the four inner terrestrial planets may have experienced secular chaotic scenario of spin-axis orientation (Laskar & Robutel 1993).For exoplanets, the evolution of planetary obliquity will directly affect the temperature distribution on planetary surface, thereby leading to a variable seasonality of terrestrial planets (Gaidos & Williams 2004).Moreover, as to hot-Jupiters, high planetary obliquity can give rise to a severe atmosphere loss rate (Nikolov & Sainsbury-Martinez 2015).The obliquity can greatly influence the atmospheric escape rate of Earth-like planets orbiting late-type M-dwarfs (Dong et al. 2019), which plays a vital role in regulating surficial habitability.
As of today, 2M0122b is the only derived oblique exoplanet, which is believed to be a planetary mass companion around the 120 Myr host star at an orbital distance of 50 AU.Thus, the diversity of planetary obliquity should be clearly understood due to planetary for-mation.Bryan et al. (2020) measured three inclinations referring to the line-of sight, where i o for the planet orbit, i p for the planetary spin and i * for the stellar spin.The longitude of ascending node of the planet's equatorial plane on its orbital plane is defined as λ p .According to the law of cosines, the most possible planet obliquity is estimated to be greater than 50 • .So naturally arises a question: how does such high obliquity origin?
Aside from indirect calculations through the spatial geometric configurations, a wide variety of observational methods were carried out to reveal the planetary obliquity.Kawahara (2016) analyzed the frequency modulation of the scattered light curve of a directly imaged exoplanet, which is induced by the planetary obliquity and orbital inclination.The modulation factor changes with difference attitudes of the spin rotations.Schwartz et al. (2016) further predicted that the planet's spin axis affects the time-resolved photometry, and east-west albedo contrast plays a vital role in constraining obliquity.Nikolov & Sainsbury-Martinez (2015) investigated the Rossiter-McLaughlin effect during secondary eclipse (RMse), and showed that for transiting exoplanets, the ratio of the ingress-to-egress radial velocity amplitudes are subjected to the planetary rotational rate and axial tilt.The synthetic near-infrared spectroscopic data were utilized to estimate the sky-projected spin axis orientation and equatorial velocity of the planet.Currently, detection of planetary obliquity is limitedly known due to the spectroscopy accuracy.RMse is only effective when the host stars are brighter than K ∼ 6 mag and large aperture telescopes (i.e., ∼ 40 m) (Nikolov & Sainsbury-Martinez 2015).In this sense, kinetic theories are very helpful to understand the evolution of planetary obliquity in advance and its predictions will be further examined by future high-precision measurements.Colombo (1966) proposed to describe the evolution of obliquity with the existence of a companion in the system, when the spin precession rate and orbital precession rate become comparable, a resonance can occur to excite a larger obliquity.The spin axis motion is characterized as precessing around its orbital axis and the total angular momentum vector.The system configuration is specified by Euler angles (θ, ϕ, ψ).The equilibrium solutions of the canonical equations of motion (Equation (6) in Fabrycky et al. (2007)) are defined as Cassini states (Colombo 1966;Peale 1969), where the planetary equatorial plane maintains constant inclination to the plane of the ecliptic, and the spin axis remains coplanar with the normal to its orbit and the normal to the ecliptic plane.Total four Cassini states may exist for a given system, each with a preferred obliquity.
Cassini states was adopted to explain moderately oblique exoplanet induced by planet-disk interaction (Su & Lai 2020) and the planet-planet interaction in multi-planetary systems (Su & Lai 2022a).The planetary obliquities can be maintained between 60 • and 80 • throughout the entire lifetime of Earth-like planets in the habitable zone of M-dwarf stars (Valente & Correia 2022), which provides friendly temperate environment and conditions for habitability.
In addition to the axial tilt of exoplanet, the spin-orbit alignment between planet and star plays a crucial part in the dynamical evolution of system.Storch & Lai (2014) intensively explored spin-orbit coupling effect between stellar spin and planet orbit, which leads to chaotic evolution of the stellar spin axis during Kozai cycles.Storch & Lai (2015) identified the origin of the chaos to be secular spin-orbit resonances overlaps.Furthermore, Campante et al. (2016) empirically constrained the angle between a planet's orbital axis and its parent star's spin axis with the asteroseismology.Vervoort et al. (2022) indicated that the period and amplitude of the tilt of a planet's rotational axis relates directly to long-term habitability of Earth-like planets.
In this work, we primarily aim to investigate the planetary obliquity evolution of S-type terrestrials orbiting M-dwarfs in binary systems.The planets that simply revolve a single star in binary are referred to S-type orbits, whereas those move around both stars are P-type orbits.As of May 2022, 154 binary star systems with 217 planets were discovered, among which more than 70 are S-type planets.For S-type systems, Kozai-Lidov model (Kozai 1962;Lidov 1962;von Zeipel 1910) is suitable for a hierarchical three-body system with more intensive orbital precession when the perturbing object is more massive than the planet.In the secular evolution of a hierarchical triple body system, the perturbation of the third body exerts on a much longer timescale than its orbital period.The Kozai-Lidov mechanism for hierarchical triple systems was studied under a wide variety of circumstances (Harrington 1968;Lee & Peale 2003;Naoz et al. 2013;Teyssandier et al. 2013;Tan et al. 2020).Li et al. (2014) explored Kozai-Lidov mechanism and characterized the parameter space that allowed large amplitude oscillations in eccentricity and inclination under the test particle limit.For a system with an eccentric outer orbit, the octupole level terms in the perturbation Hamiltonian become significant, thereby triggering EKL effect.The planetary eccentricity increases approximately one, and the flips of orbital inclination can lead to a retrograde orbit (Naoz et al. 2013).Furthermore, the orbital flips of S-type planets were found to origin from EKL (Huang & Ji 2022).Here we will explore similar flips of planetary spin-axis due to EKL.
Tidal dissipation and the dynamical torque from the adjacent planets or disk both play critical roles in the spin evolution.Wang et al. (2019) adopted a comprehensive scenario in the restricted three-body system, including the equilibrium tide and the third body perturbation, and showed that the obliquity librates for a long time or decays slowly down to zero.The planetary angular momentum, along with a positive dynamical loop consisting of orbital migration and obliquity tides, can well produce Ultra-Short-Period planets (Millholland & Spalding 2020).Formation of oblique super-Earths and several hot Jupiters can be further explained with tidal dissipation and spin-orbit capture (Su & Lai 2022b).
In the real situation of a hierarchical system when the stellar companion is more massive than the perturbed body and the secondary's eccentricity is non-negligible, the non-restricted EKL mechanism will work to trigger the planet's evolution.Based on the conservation of angular momentum, in this work the mixed scenario between EKL and equilibrium tide is adopted to reveal a diverse secular evolution paths for planetary obliquity.The comprehensive evolution timescale of the obliquity depends on the masses, orbital elements of all bodies, and the tidal factors of the planet throughout the secular evolution.The angular momentum transferred from the distant companion to the perturbed body can sustain or excite the planetary obliquity and contend with the tidal decay effect.
To better understand the nature of planetary obliquity of exoplanets, we utilize N-body package MERCURY-T (Bolmont et al. 2015) to conduct more than 2000 numerical simulations that each run evolves at a timescale of roughly 10 7 yr to obtain two typical evolution paths of planetary obliquity corresponding two regions in the Hamiltonian level curves.Subsequently, we derive the relationship of equilibrium timescale versus the planetary mass m 1 , the initial orbit a 1 and the timescale ratio r t .In the parameter space of semi-major axis, we present the diagram of maximum obliquity to estimate the flip ratio of obliquity.Here we further investigate reversal conditions of planetary spin axis involved in the EKL and tide effect simultaneously, which are employed to address the evolution of S-type planets, such as HD 42936 Ab, GJ 86 Ab and τ Boot Ab.
This work is structured as follows.In Section 2, the theory of the mixed mechanism is described under the frame of central body equatorial coordinate, including theoretical analysis of the relationship between initial timescale ratio r t and the equilibrium timescale.Section 3 describes evolution paths for general S-type terrestri-als, and provides the relationship between the maximum planetary obliquity, the equilibrium timescale and r t .The obliquity flip ratio is further calculated for a wide variety of initial conditions.Section 4 presents the secular evolution of several potential oblique S-type planets, and predict the possible value of the equilibrium obliquity.In final, we summarize major outcomes.

THE MIXED SCENARIOS
In the secular evolution of triple body systems, the perturbation arising from the third body acts on a time scale much longer than the orbital period.If the perturbation of a circular outer orbit works, the Hamiltonian of this system can be expanded to the quadrupole level, which is referred to Kozai-Lidov mechanism as aforementioned.Here we adopt the Kozai-Lidov mechanism and the equilibrium tide to characterize the total evolution for S-type planets.It is interesting that the rotation angular momentum, the planetary orbital angular momentum and the stellar companion's angular momentum are all participating to transfer to each other.The geometry configuration of the system is plotted in Figure 1, where the reference plane is the equatorial plane of the host star m 0 , m 1 and m 2 are the perturbed planet and the outer stellar companion.− → G 1 and − → G 2 are angular momentum vectors of the inner and the outer orbits.− → G is the total angular momentum vector of the orbital motion.
− → N p is the angular momentum of planetary spin.

Quadrupole Force vs Tide
The classical Kozai-Lidov mechanism shows that the inclination and eccentricity of inner test particle oscillate over secular evolution in a hierarchical system (von Zeipel 1910;Kozai 1962;Lidov 1962).With the secular approximation, the inner and the outer orbits only exchange angular momentum, thus the semi-major axis of orbits do not change.When the SMA ratio α r = a 1 /a 2 is a small parameter, the perturbation term of full Hamiltonian can be expanded as a power series in α r (Naoz 2016): where k 2 is the gravitational constant (with the mass unit of M ⊕ and the length unit of au), r 1 is the distance between m 0 and m 1 , r 2 is the distance between the center of mass of the inner binary and m 2 , P j is the Legendre polynomial, Φ is the angle between the vectors r 1 and r 2 , and the subscripts j = 1, 2 represent the inner and outer orbits, respectively.
In the non-test particle approximation under the classical Kozai-Lidov mechanism, the eccentricity and the inclination oscillate regularly in a well-defined timescale t kl (Antognini 2015): This relationship was derived under the consideration of the equation of motion of ω, argument of perihelion, by integrating between the maximum and minimum eccentricities.
Meanwhile, we employ the constant time-lag equilibrium tide at arbitrary rotation (Hut 1981;Leconte et al. 2011;Wang et al. 2019) from the start of integration in the present work.Because the planetary eccentricity will be triggered to an extremely large value by Kozai-Lidov effect, which may result in unaccepable outcomes for the low-order expansion of eccentricity in the constant phase delay model.When the orbit plane does not coincide with the equatorial plane of the planet, Wang et al. (2019) used the Oxyz coordinate system to describe the equilibrium tide at arbitrary rotation and the equations of motion for the constant time-lag tidal model.
We choose a central body coordinate system, and the stellar equatorial plane is treated as the reference plane.The origin is the host star, where x-y plane runs parallel to the equatorial plane of the star, z-axis is positively aligned with the rotation axis of the star.The Oxyz coordinate system is in line with the frame illustrated in Figure 1.When the tidal potential function expands to order 2 (Kozai 1965), the tidal acceleration is derived as (Wang et al. 2019): where ŵ is the unit vector of orbital angular momentum, r is the vector from the star to the center of mass of the planet, and θ is the unit vector along the orbital velocity.G is the gravitational constant, r is the distance from the planet to the center of the star, k 2 is the 2th order Love number, θ is the instantaneous orbital angular velocity.τ denotes the time delay factor in the constant time lag tidal model (Hut 1981).Ω is the longitude of ascending node of the planetary orbit, and i is the orbital inclination of the planet's orbit relative to the Oxyz system.Ω p is the planetary rotational velocities.I and Θ, refer to the inclination and the longitude of ascending node of the planet's equatorial plane with respect to the Oxyz system, respectively.We can calculate equations of rotational motion based on the total angular momentum conservation: where ⃗ C is the total angular momentum vector, and the moment of inertia I p = m 1 R 2 p r 2 g , with r 2 g the square of radius of gyration (constant).Total angular momentum ⃗ C and orbital angular momentum H 1 , H 2 can be calculated by the initial orbital elements: a 0 , e 0 , i 0 , Ω 0 , I p , Ω p,0 , I 0 , Θ 0 .In the evolution, the rotational velocity of planet Ω p can be estimated with orbital parameters according to angular momentum conservation.Hence, I and Θ can be further derived through three components of the planetary spin vector.
The evolutionary timescale of I is defined as (Wang et al. 2019): where I p is the momentum inertia of the planet, T t is a typical tidal dissipation timescale, which keeps constant and can be calculated by: τ ⊕ is the time delay factor of Earth.In the following, we define the planetary obliquity as θ p , which is the angle between the orbital angular momentum vector and the rotational angular momentum vector (Figure 1).Note that Equation (1) and the secular perturbation timescale t kl in Equation ( 2) are defined in the Jacobian coordinate system (Naoz 2016).For a host star with m 0 = 0.1 M ⊙ and the terrestrial planet with a mass of 1 M ⊕ , the center of mass of them will be inside the host star if the semi-major axis a 1 < 0.2 AU.The difference of coordinates between the Jacobian frame and the central body frame can be neglected.Here we employ the astrocentric coordinate as the system frame of reference.
According to Colombo (1966), the spin precession rate and orbital precession rate is comparable to excite planetary obliquity.As the tidal dissipation timescale t I is much greater than t kl induced by the Kozai-Lidov mechanism.Blue and red solid lines in Figure 2 simply focus on the comparison of t kl and t I for a wide variety of initials m 1 , e 2 , a 1 and a 2 according to the theoretical equations, which is also considered to be a pair of initial conditions.Figure 2 indicates that when a 1 increases, the timescale t I can rise but t kl decreases.By contrast, t kl grows when a 2 increases.For a 1 = 0.02 -0.04 AU, the line of t I crosses t kl at roughly 10 4 yr.The tidal dissipation timescale goes up along with the planetary mass, and t kl is not so sensitive to the eccentricity e 2 of the secondary.In addition, Figure 2 exhibits that t kl remains larger when e 2 = 0.2, t I and t kl are equivalent (the intersections of two curves) at 5 × 10 3 -3 × 10 4 yr.When e 2 = 0.4, the EKL timescale is smaller, and the intersection range of t I and t kl is 4 × 10 3 -3 × 10 4 yr.For m 1 = 1 M ⊕ , the timescale of tide is larger in comparison with that of m 1 = 0.2 M ⊕ , and the crossing point of t I and t kl varies from 7 × 10 3 to 5 × 10 4 yr.

Octupole Force vs Tide
Unlike the improved tidal triple body model (Wang et al. 2019), here we consider the combined effects of the octupole level perturbation and the equilibrium tide for non-restricted hierarchical systems.In such model, the total angular momentum is conserved, which consists of the orbital angular momentum of planet, that of the secondary star and the angular momentum of planet' spin.The three components of angular momentum are mutually coupling over secular evolution.However, for the restricted model, the total angular momentum of the planet is assumed to be constant.
As Naoz (2016) pointed out, the EKL secular evolution timescale is sensitive to e 2 of the secondary star, which is difficult to be directly calculated.Thus we numerically investigate secular evolution of the systems with N-body simulations.
In Figure 2, red dots exhibit the comprehensive evolution timescale of planetary obliquity corresponding to different comparison between t I and t kl .Vertical coordinates of the red dots represent the mean evolution timescale of the planetary obliquity during the first integration output step from t = 0 to t = δt.In the runs, the output step δt is set to be 100 yr.When the vertical coordinates of the dots fall down to be negative (gray lines) , the planetary obliquity goes down or evolves to be retrograde.For a relatively distant secondary star, θ p will decrease initially with a 2 > 10 AU and a 1 = 0.02 -0.08 AU.

SECULAR EVOLUTION OF THE OBLIQUITY
Here the principal goal for this work is mainly to explore the spin evolution of terrestrial planets, as the planetary obliquity plays a crucial role in habitability, therefore our study will cast light on spin evolution for habitable planets.Over two hundred S-type planetary systems were discovered so far, however, there are only 10 binaries that own eccentricities, which those of the secondary stars are detected to be greater than 0.2.Note that the octupole force frequency is proportional to the eccentricity e 2 of the perturber under EKL scenario.

Numerical setup
To extensively investigate the spin evolution for terrestrial planets, here we adopt the standard equilibrium tidal model and constant time lag model given by MERCURY-T (Bolmont et al. 2015) with numerical simulations.Here all runs refer to the N-body simulations.The constant time lag model consists on the assumption that the bodies under consideration are made of a weakly viscous fluid (Alexander 1973).MERCURY-T models the tidal forces between the star and the planets but we neglect the tidal interaction between planets, then adds the forces in the acceleration of orbital motion and the spin angular momentum.The tidal torque exerted by planet and the star is explicitly written in Bolmont et al. (2015).
Here host stars are considered to be non-evolutionary in the simulations.Mass of the secondary star is set to be 0.2 M ⊙ .The initials of a 1 are in the range 0.02 -0.13 AU, which keeps the orbital periastron of a rocky planet outside of the Roche limit of a host star (Liu et al. 2013).
For rocky planets, the lower value of Roche limit is 1.58 R ⊙ .The separation of the binaries a 2 are given in [2,5,10,15,20] AU.According to t I and t kl in Equation ( 2) and ( 7), t I can always be comparable with t kl for each set of a 2 .To maintain the system stability, we choose e 1 = 0.2, e 2 = [ 0.2, 0.4] for terrestrial planets orbiting M-dwarfs, and set e 2 = [0.6,0.8] for solar-type host star, respectively.The initial inclinations are given to be i 1 = 50 Recently, the stellar spin was proved to has large implications on the dynamics of planetary systems (Becker et al. 2020;Chen et al. 2022;Faridani et al. 2023), our calculations demonstrate that the radius of M-dwarf with age of 700 Myr and a mass of 0.1 M ⊙ will decrease from 0.1238 R ⊙ to 0.1236 R ⊙ .The terrestrial planet's orbital angular momentum variation due to the stellar spin evolution over 2 × 10 7 yr is less than 4 × 10 −6 , and the final planetary obliquity differs from that of the non-evolving model by ∼ 1 • .Thus we adopt the nonevolving M-dwarf host model, which will not bring about significant influence to our results.
In the non-evolving host body model, the time delay factor in the constant time lag tidal model is set as τ ⊕ = 698 s (Bolmont et al. 2015).For terrestrials, we select the love number k 2 = 0.305, the radius of gyration r 2 g = 0.3308, and the planetary radius R p = R ⊕ .The initial planet rotation period is set to be 24 hr, and the initial value of θ p is assumed to be 10 • in the simulations.
The former study revealed that general relativity (GR) may have influence on triples' evolution, which may result in resonant-like dynamics (Naoz et al. 2013;Liu et al. 2015) or destabilize the existing resonance (Hansen & Naoz 2020).We conduct additional runs by considering GR effect with the same initials as Figure 4.For a terrestrial planet at a 1 = 0.04 AU with the GR and non-GR models, the difference in final rotation period is 0.5 hours and the final planetary orbital angular momentum differs by 6 × 10 −5 .To speed up the integration, here we simply adopt the non-GR model in the simulations.
In this work, we conduct over 2000 simulations that each evolves for 10 Myr (∼ 10 4 t kl ) with MERCURY-T.When all runs reach 10 7 yr, the simulation is ceased once the equilibrium state of planetary obliquity occurs after 10 7 yr, which is the stopping condition.

Typical Evolution Paths
According to the origin of Cassini States described in this work and Equation (8) in Su & Lai (2022b), we first plot the Hamiltonian level curves in the phase space of cos θ p −ϕ to theoretically describe the spin evolution feature, ϕ is the angle to represent the precessional phase of spin vector about the orbital angular momentum vector.We find the dynamical structures similar to Su & Lai (2022b) as shown in Figure 3, including up to three kinds of typical oscillation modes, which refer to three kinds of Cassini States.For region I in Figure 3, the obliquity suffers excitation below 90 • .θ p experiences a larger oscillation amplitude and the spin-axis flip with θ p > 90 • in region II.After flip, θ p stably librates between 40 • and 60 • .In region III, the planetary obliquity is initially stirred up to θ p > 90 • , finally evolves to the stable retrograde spin orientation.
The Hamiltonian level curves indicate that the planet spin evolution dynamics is greatly affected by the initial conditions, the equilibrium obliquity corresponding to Cassini State II that decreases with the precession ratio of the orbit and spin-axis.Based on the theoretical results from the Hamiltonian level curves, we then mainly investigate the typical oscillation cases of region I and II, to summarize the relationship between initial conditions and the final evolution feature.
Figure 4 presents a typical set of obliquity evolution with respect to diverse initial orbits for terrestrials, where a 1 of the planet varies from 0.02 to 0.2 AU.Panel (a)-(d) follow the evolution path of region I in Figure 3, and panel (e)-(f) follow the region II.We also perform the time-series evolution of inclination and eccentricity in Figure 5 to show the effect of the EKL strength, and the excitation effect of EKL is greater when the planet and the secondary star are closer.The eccentricity can be excited to 0.7 and the orbital inclination decreases to lower than 40 • when a 1 = 0.1 AU.The oscillation amplitude of e 1 and i 1 are increasing with a 1 .The variation of semi-major axis is too little during the evolution to capture the orbit decay, thus we could not see the noticeable circulation.
The evolution results can be further distinguished by r t , as shown in Figure 6.Here we perform the numerical results with respect to four groups of m 1 = 0.2 M ⊕ , 1.0 M ⊕ and e 2 = 0.2, 0.4.With the increase of r t , the maximum obliquity also increases.Figure 6 shows that the planetary obliquity without flip when r t < 1, while the obliquity of Earth-mass planet could flip when r t ranges from 1 to 10 3 , and θ p will always flip when r t spans from 10 3 to 10 6 .
It is noteworthy that when 10 < r t < 10 3 , θ max goes up to a local bulge and then declines to a valley.It is speculated that the extreme planetary eccentricity originates from EKL mechanism when r t gets larger, and the tidal effect is dramatically enhanced owing to a much closer pericenter to the host star.Under the combined scenarios, the planetary obliquity undergoes an inhibited growth.On the other hand, we further compare  the obliquity evolution due to e 2 .When e 2 rises to 0.4, the more obvious obliquity excitation occurs, and the maximum θ p enters the flip region with smaller logr t .For r t > 10 4 , the maximum obliquity converges to ∼ 130 • .In total numerical simulations, we calculate the evolution paths for terrestrial planets with m 1 = 0.2 ∼ 8 M ⊕ , where the general evolution trend is consistent.When m 1 grows to 1 M ⊕ , the local maximum obliquity of the bulge ascends to 100 • with a 2 = 20 AU.Additional simulation results reveal that when m 1 increases to 4 M ⊕ and 8 M ⊕ , the local maximum obliquity can arrive at 105 • and 110 • , respectively.It can be inferred that more runs may locate at the region of flip with larger planetary masses for terrestrials around M-dwarfs.
We carry out additional N-body simulations with MERCURY-T to study more general S-type systems with an extremely higher e 2 .Currently, there are over five S-type binary systems that are observed to have an eccentric outer companion with e 2 > 0.5, thus we can further investigate the evolution of θ p in these cases.The simulations show that the terrestrials orbiting Mdwarfs can be scattered or captured by the secondary with e 2 = 0.8.Furthermore, we conduct a couple of runs to explore terrestrials orbiting solar-type stars in comparison with those about M-dwarfs, where 55% of them can survive under the comprehensive effect of EKL and tidal dissipation.
Aside from typical oscillation cases presented in Figure 4, we also find chaotic evolution behaviour.Storch & Lai (2014) demonstrated the stellar spin-orbit angle can have a chaotic region, and similar chaotic behaviour of the planetary obliquity exist in our work.These chaotic cases always occur when a 2 = 2 AU and e 2 > 0.6, the extremely excitation of eccentricity leads to the instability of the system.Comparing with the stellar spin-orbit angle chaotic behaviour, the chaotic results of planetary obliquity in our work origin from spin-orbit coupling effect between the planet spin and the planet orbit, rather than between the stellar spin and the planet orbit.

The Equilibrium Timescale
In Section 3.2, θ p remains a quasi-equilibrium state for over 10 7 yr prior to a decrease to zero, and the timescale of entering the equilibrium obliquity is denoted as t eq .We then calculate the ratio between the orbital period and rotation period P r over the entire evolution.For all terrestrial cases with m 1 = 1 M ⊕ and a 1 = 0.1 ∼ 0.2 AU, P r will fall down to 6 at the time of t eq , which seems to be in a high-order spin-orbit resonant-like state.The total angular momentum of the planet, the summation of − → G 1 and − → N p , also drops down to a constant.We infer that the planetary angular momentum is trapped in the high-order spin-orbit resonant-like state.Su & Lai (2022a) showed that the equilibrium obliquity ranges from −90 • to 90 • for Cassini State I, II and IV, whereas the obliquity varies from 160 • to 180 • for Cassini State III.The related precession ratio of spin vector and orbital angular momentum vector is in (10 −2 -10 1 ).Here the equilibrium θ p performs the oscillation state consisting of Cassini State II-evolves between 40 • and 60 • , which agrees with the typically evolution in our simulations.Winn & Holman (2005) addressed that Cassini state II is the most favorable configuration to maintain a significant obliquity.In this study, we adopt a non-restricted model, the initial planetary mass is set to be 0.2 -8 M ⊕ , then we perform in a larger parameter space of r t spanning from 10 −2 to 10 7 .Our findings unveil a new equilibrium value of the maximum obliquity θ max close to 130 • when r t > 10 4 (Figure 6).
In the equilibrium stage, the initial planetary orbit residing in the habitable zone of M-dwarf is not shrunk to approach the Roche-limit, thus the spin-axis orientation remains relatively stable, thereby helping maintaining the planetary habitability.
Figure 7 shows the contours of equilibrium timescale t eq in the parameter space a 1 -a 2 , where t eq evenly distributes in the dimension of a 1 but varies with a 2 and e 2 .Next, we observe that t eq is always less than 10 5 yr for a 1 = 0.09 AU (see the dashed lines in each panel).For a 2 < 20 AU, t eq is approximately proportional to the initial orbit of planet, which primarily governs tidal timescale.For a 1 = 0.02 AU and a 1 > 0.1 AU, the obliquity equilibrium timescale keeps unchanged for the variations a 2 at same a 1 .When the semi-major axis meets 0.04 < a 1 < 0.08 AU, a 2 = 20 AU, the tidal timescale is equivalent to the secular perturbation timescale.
We further investigate the relationship between t eq and the timescale ratio r t as shown in Figure 8, where Figure 8 (a)-(b) shows the planetary mass of m 1 = 0.2 M ⊕ while m 1 = 1M ⊕ in Figure 8 (c)-(d), with respect to e 2 = 0.2 and 0.4, respectively, for each panel.We then utilize a maximum likelihood (ML) function and numerically optimize the relation.The fitting model is logt eq = m M L × (logr t ) + b M L .The fitting results of t eq and t I /t kl are provided in the panel legend.When a 2 = 2 or 5 AU, the tidal torque and the EKL effect are intensive for (t I /t kl ) < 10 2 , the equilibrium timescale of the obliquity much arises from the tidal dissipation, thus the blue dots deviate to the lower end.

The Obliquity Flip Ratio R flip
In Section 3.1, we describe two typical evolution paths of obliquity, and cases follow region II can always trig-Figure 5. Evolution of the eccentricity and inclination for S-type terrestrials with a mass of 1 M⊕ and semi-major axis smaller than 0.1 AU, other initials same as Figure 4.The eccentricity can be excited to 0.7 and the inclination decreases to below 40 • when a1 = 0.1 AU.The oscillation amplitude of e1 and i1 are increasing with a1.When a1 > 0.1 AU, the upper and lower limit e1 and i1 in the evolution are not sensitive to a1 any more.
ger the flip of spin-axis.We extend the distance of the second star to 45 AU in this section, and investigate the obliquity flip ratio in the parameter space of a 1 -a 2 .In Figure 9, over 700 simulations are conducted for S-type terrestrials with m 1 = 0.2 ∼ 8 M ⊕ , accompanied with a constant perturber of 0.2 solar mass, in which each panel is plotted with a grid resolution of 10 × 18.The flip ratio is calculated by the fraction of the simulations that undergoes flip.
The flip ratio in Figure 9(a)-(b) are 0.49, 0.55, 0.67, 0.72, respectively.We conclude that the reverse of the spin axis can be very common for S-type terrestrials with a distance companion.When the planetary mass increases, the flip ratio increases with the same initial semi-major axis, thus the flips get easier with higher m 1 /m 2 .This will provide the clues to predict what kind of exoplanets can evolve to be head-down, like Venus in our solar system.We also change the eccentricity of the companion, but the flip region of obliquity is not much sensitive to e 2 .For the Earth-mass planet, the flip ratio spans from 0.52 to 0.56 when e 2 = 0.2 -0.6.In this section, several S-type planets are reported to figure out whether they rotate with the oblique attitude under the combined scenarios of EKL and equilibrium tide.Considering the timescale of secular evolution, we select the terrestrial planet HD 42936 b and the Jupiterlike GJ 86 Ab and τ Boot Ab as case study.In Table 1 and 2, we summarize essential parameters for three potential oblique S-type systems.The primary star HD 42936 is a solar-like star with mass of 0.87 M ⊙ at a distance of 48.9± 0.6 pc, while the semimajor axis of the secondary star is simply 1.22 ± 0.02 AU, and the eccentricity is relatively high with e 2 = 0.594.Such a compact and hierarchical system naturally brings secular perturbation scenario and tidal evolution into our mind.Barnes et al. (2020) concluded that the EKL mechanism, as one of the most likely formation models, can drive HD 42936 b move to current orbital configuration, in which the eccentricity of planet can be initially stirred up, and then shrink to a circular orbit when the tidal timescale is comparable to or shorter than t kl .Hence, HD 42936 is a good example to explore the spin evolution with our mixed mechanism and reveal the dependence of the spin-axis evolution tendency on the initial conditions.HD 42936 b was discovered by radial velocity, probably moving at an inclined orbit.Thus in this work, the initial semi-major axis of this planet is assumed to be 0.06 -0.12 AU, and the initial inclination ranges from 10 • to 170 • .When 40 • < i p < 140 • , HD 42936 b will collide with the primary star or be scattered out of the system, before that the obliquity θ p can be excited to 150 • .For the survivals, θ p can be triggered to the maximum value of 80 • .Figure 10 shows the obliquity evolution for HD 42936 b with the initial prograde orbits, while Figure 11 stands for those results for the retrograde orbits at beginning.

GJ 86 Ab
Aside from terrestrials in binary systems, we again investigate the spin evolution of hot-Jupiters with the mass ranging from 0.3 to 8 M J .The simulation results show that the obliquity can be triggered to retrograde, and the maximum obliquity can arrive at 160 • , being indicative of that the equilibrium obliquity may be sensitive to the planetary mass.
GJ 86 is a nearby S-type system at a distance of 10.9 pc from the solar system, and hosts a giant planet GJ 86 Ab (Queloz et al. 2000) at the orbital distance of 0.11 AU, the outer companion is a white dwarf in the eccentric orbit with a 2 = 21 AU, e 2 = 0.4 (Mugrauer & Neuhäuser 2005;Lagrange et al. 2006).The orbital elements and masses of both companions are measured (Zeng et al. 2022)   2016; Brandt 2018).The current distance between the binary is 21 AU, whereas the mutual closest approach may reach 9 AU (Zeng et al. 2022), being indicative of that the gas-giant may origin from the disk truncation.
The secular perturbation also participates in the obliquity evolution due to close separation of the binary.Our simulations suggest that GJ 86 Ab have a large possibility to follow obliquity evolution path II, and maintain the head-down rotation attitude.GJ 86 is selected as one of the candidates for CHES (Closeby Habitable Exoplanet Survey) mission (Ji et al. 2022), which will observe nearby solar-type stars to hunt for terrestrial planets in the habitable zones at an ultrahigh resolution via astrometry.The future observations will extensively provide the orbital elements and true mass of GJ 86 Ab, which helps place constraints on the dynamics of spin evolution.

τ Boot Ab
The giant planet orbiting τ Boot is also one of the best-known exoplanets in the nearby stars from our solar system.The planet τ Boot Ab was first discovered in 1996 (Butler et al. 1997).Brogi et al. (2012) utilized radial velocity measurements to determine an orbital inclination of i p = 44.5 ± 1.5 • and a true planetary mass of 5.95 ± 0.28 M J .They presented the atmospheric characterization of this non-transiting planet and showed the planetary temperature distribution de- creases towards higher altitudes.Table 1 and 2 summarize the essential parameters for τ Boot.For a 1 = 0.046 AU, a AB = 45 AU, we still believe the perturbation from the outer companion can excite the orbital eccentricity and spin obliquity of τ Boot Ab, as the eccentricity e 2 of the secondary remains extremely high.Our numerical results indicate that θ p follows the evolution path II, and the maximum value can reach above 100 • , then falls to the equilibrium stage before 100 Myr (Figure 12).

CONCLUSIONS AND DISCUSSION
In this work, we explore the spin evolution of terrestrial planets as the planetary obliquity plays a vital role in habitability.As the binary systems of stars that harbor planets are very common, thus secular perturbation from the secondary may greatly lead to a diverse plan-etary evolution.Here we adopt a comprehensive model composed of the octupole level secular perturbation and equilibrium tide to study dynamical evolution of S-type planets' obliquity, which may provide essential clues for elucidating spin evolution for habitable planets.
We first calculate the comprehensive evolution timescale of θ p , which relates very much to the Kozai-Lidov timescale t kl and tidal timescale t I .We further plot the Hamiltonian level curves of the spin dynamics to theoretically derive the evolution paths of terrestrials' planetary obliquity, and also conduct numerical simulations to extensively explore the evolution of planetary obliquity induced by the coupling effect due to EKL resonance and the equilibrium tide, the planetary obliquity equilibrium state varies as a function of r t .In addition, the maximum θ p during evolution can be distinguished  According to the numerical results in the parameter space of semi-major axis a 1 and a 2 , we present diagrams of the equilibrium timescale and maximum obliquities.The equilibrium timescale t eq evenly distributes in the dimension of a 1 , and the timescale is always smaller than 10 5 yr when a 1 < 0.09 AU.For a 2 < 20 AU, t eq is approximately proportional to the initial orbit of the planet, which mainly determines the tidal timescale.Then we derive the linear relationship between t eq and the timescale ratio r t with the maximum likelihood method, logt eq is approximately proportional to logr t and the slope of fitted lines rises with m 1 .
Moreover, we estimate the ratio of obliquity flip in the extensive parameter space a 2 < 45 AU.We conclude that the reverse of spin axis appear to be common for S-type terrestrial planets accompanied with a distant massive companion.When the planetary mass ascends, the flip ratio grows with the same initial semi-major axis, thus the flips become easier for a larger m 1 /m 2 .To employ our comprehensive model in the observed systems and predict the evolution of potential oblique planets, then we further numerically investigate three cases of HD 42936 b, GJ 86 Ab and τ Boot Ab, to understand their evolution paths of θ p .
Aside from the planetary spin-orbit flip, the flip of the star-spin-orbit angle induced by the EKL mechanism were noted in the literature previously (Naoz et al. 2011(Naoz et al. , 2012;;Petrovich 2015a;Lai et al. 2018;Stephan et al. 2018), and the original high stellar spin-orbit obliquity will be damped for the cooler star (Lai 2012;Dawson 2014).As above-mentioned, Storch & Lai (2014) intensively studied the spin-orbit coupling effect between the stellar spin and the planet orbit, which gives rise to the chaotic evolution of stellar spin axis during Kozai cycles.
Effects of tidal dissipation and stellar spin-down will also influence the final distribution of spin-orbit misalignment angles of hot-Jupiters.Storch & Lai (2015) explored the origin of this chaotic stellar spin behaviour, they identified secular spin-orbit resonances and the res-onance overlaps are responsible for the chaos.Key parameters including the adiabaticity parameter, the ratio of the Kozai-Lidov nodal precession rate and the stellar spin precession rate are proved to dominate the degree of chaos.In our research, similar chaotic behavior of the planetary obliquity also exist in several simulation cases.It was also demonstrated that before being destroyed by their stars either Roche-limit crossing or engulfment during stellar expansion, there could be diverse results on the final stellar-spin orbit angle (Petrovich 2015a,b;Stephan et al. 2018;Angelo et al. 2022).
It was recently highlighted that chaotic tide plays a key role in the evolution of planetary systems (Wu 2018;Vick et al. 2019).In this case, in eccentric orbits, the planet's orbital energy is converted into internal fluid energy via tidal stretching and compression at periastron passage, causing the orbit to circularize over time.For Earth-like planets, a multilayered internal structure would increase the efficiency of tidal dissipation and affect the climate and habitability of the planet.For Jupiter-mass planets, increasing the rate of planetary orbital energy loss accelerates the tidal dissipation process and faster entry into the proposed equilibrium state.
With the improvement of observational accuracy, the space-borne telescopes, such as JWST (James Webb Space Telescope) (Gardner et al. 2006), TESS (Transiting Exoplanet Survey Satellite) (Ricker et al. 2015), PLATO (Planetary Transits and Oscillations of stars) (Rauer et al. 2014), ARIEL (Atmospheric remotesensing infrared exoplanet large survey) (Tinetti et al. 2021) and CHES (Closeby Habitable Exoplanet Survey) (Ji et al. 2022), will release a bulk of detailed information about the planetary three dimension orbital configuration, the planetary mass, the presence of atmosphere, and spin orientation, thereby further identifying the stable oblique planet in habitable zone.There are prospects for constraining exoplanet obliquities in the coming years, such as using high-resolution spectroscopy to obtain the projected rotational radial-velocity of the planet (Snellen et al. 2014;Bryan et al. 2018) based on the Rossiter-Mclaughlin effect during the second transit eclipsing.The planetary obliquity is also an important indicator for the seasonality of habitable terrestrials and the atmosphere loss rate of hot-Jupiters.This can be further examined by upcoming observations.In future work, we will further explore the origin of these chaotic behaviors, the stellar spin evolution, and the connection between planetary obliquity and stellar obliquity as well.

Figure 1 .
Figure 1.Geometry of the hierarchical system.−→ G1 and − → G2 are angular momentum vectors of the inner and the outer orbits.− → G is the total angular momentum vector of the orbital motion.−→ N p is the rotation angular momentum.

Figure 2 .
Figure 2. The relationship between tI and t kl obtained from the theoretical equations.In panel (a)-(b), m1=0.2M⊕, in (c)-(d), m1=1 M⊕. e2 is set to be 0.2 in panels (a)-(c) and e2 =0.4 in (b)-(d).Dotted lines are the comprehensive evolution timescale of the planetary obliquity during the first output step of simulations.The comprehensive evolution timescale of the obliquity locally decreases when tI > t kl for a2 > 10 AU and a1 = 0.02 -0.08 AU under EKL.

Figure 3 .
Figure 3. Three examples of Hamiltonian level curves for the spin evolution of terrestrial planets around the M-dwarf.The planetary semi-major axis ranges from 0.05 AU to 0.1 AU, and the precession frequency ratio noted as η by (Su & Lai 2020) is selected as [0.12, 0.45, 1.2] respectively.The different Cassini states are marked out with labels I, II, III.

Figure 6 .
Figure 6.Distribution of the maximum obliquity θmax in the dimension of timescale ratio rt.The semi-major axis a2 ranges from 2 to 20 AU.Regions with or without flip are also filled with pink and green respectively, with 90 • as horizontal boundaries.
with the observations from radial velocity and high angular resolution imaging along with the absolute astrometry data from Hipparcos (Perryman et al. 1997) and Gaia (Gaia Collaboration et al.

Table 1 .
Parameters for binary stars in three potential oblique S-type systems

Table 2 .
Physical parameters for planets in three potential oblique S-type systems.