Spin and orbital dynamics of planets undergoing thermal atmospheric tides using a vectorial approach

Earth-mass planets are expected to have atmospheres and experience thermal tides raised by the host star. These tides transfer energy to the planet that can counter the dissipation from bodily tides. Indeed, even a relatively thin atmosphere can drive the rotation of these planets away from the synchronous state. Here we revisit the dynamical evolution of planets undergoing thermal atmospheric tides. We use a novel approach based on a vectorial formalism, which is frame independent and valid for any configuration of the system, including any eccentricity and obliquity values. We provide the secular equations of motion after averaging over the mean anomaly and the argument of the pericenter, which are suitable to model the long-term spin and orbital evolution of the planet.


Introduction
In general, Earth-mass planets are believed to be composed of a large rocky mantle covered by a thin atmospheric layer (e.g., Komacek & Abbot 2019;Wordsworth & Kreidberg 2022), as for the Earth, Venus, and Mars in the Solar System.The molecules that make up the atmospheres can absorb part of the radiation they receive from the host star, giving rise to local temperature gradients, which in turn create pressure variations.As a result, large-scale periodic mass redistribution inside the atmosphere occurs, while it attempts to return to an equilibrium state, known as thermal atmospheric tides.
Observations on Earth and Mars show that the pressure variations can be essentially decomposed in a diurnal and a semidiurnal tidal wave (e.g., Withers et al. 2011;Auclair-Desrotour et al. 2017).These oscillations are one of the most regular meteorological phenomena on Earth, and they are easily detectable by any station around the world (e.g., Chapman & Lindzen 1970).The thermal inertia of the atmosphere introduces a delay between the stellar heating and the thermal response, which creates an asymmetry in the mass redistribution with respect to the substellar point.As a consequence, the gravitational pull of the star exerts a torque on the atmosphere.The diurnal wave has no net torque, but the semidiurnal wave gives rise to angular momentum exchanges within the atmosphere (e.g., Correia & Laskar 2003a).
Thermal atmospheric tides garnered special interest with the discovery of the retrograde rotation of Venus (e.g., Carpenter 1964).Because the atmosphere and the mantle are usually well coupled by friction at the surface, the variations in the angular momentum of the atmosphere are then transferred to the mantle of the planet, modifying its spin.To explain the peculiar rotation of Venus, Gold & Soter (1969) thus proposed that it can be the result of a balance between bodily tides (i.e., gravitational tides raised in the mantle), which drive the planet toward synchronous rotation, and thermal tides, which drive it away.This effect is negligible on Earth, because it is too distant from the Sun; however, for Earth-like planets in the habitable zone of Kdwarf stars, thermal tides can lead to asynchronous equilibria similar to that of Venus (Leconte et al. 2015).
Formation studies predict that Earth-mass planets are common around main sequence dwarf stars (e.g., Schlecker et al. 2021).Indeed, there is already a large population of detected low-mass exoplanets, and their number is likely to grow rapidly (e.g., Winn 2018).These planets are found in a wide range of orbital configurations, from very eccentric to compact multibody resonant orbits (e.g., Borucki et al. 2013).We thus need a correct modeling of their spins and orbital dynamics.The estimates for the tidal evolution of a planet are based on a general formulation of the tidal potential (e.g., Kaula 1964).The classical expansion of this potential in elliptical elements depends on the chosen frame and introduces multiple index summations, which can lead to confusion and errors in the equations of motion.Moreover, mistakes such as neglecting energy or momentum conservation considerations are more easily done (Boué & Efroimsky 2019).For bodily tides, it has been shown that the equations of motion are more easily expressed in terms of angular momentum vectors (e.g., Correia & Valente 2022).Therefore, in this paper, we aim to also use these vectors to study the dynamics of a planet undergoing thermal tides.This formalism is independent of the reference frame and allows us to simply add the contributions of multiple bodies in the system.
In Sect.2, we obtain the tidal potential of a planet whose atmosphere is excited through thermal forcing.In Sect.3, we derive the equations of motion to study the spin and orbital evolution of a planet-star system under the action of thermal tides.In Sect.4, we average the equations of motion over the mean anomaly and the argument of the pericenter, which provide simpler expressions that are easy to implement and suitable for longterm evolution studies.In Sect. 5 and Sect.6, we provide the equations of motion in the simplified cases of small eccentric-

Thermal atmospheric tides
We consider a planet with total mass m, which is composed of a perfectly spherical mantle with radius R, and a thin atmosphere with mass M. We assume that the mantle is completely rigid, while the atmosphere can be deformed owing to the thermal forcing from the radiation of a nearby star.

Atmospheric tidal potential
The gravitational potential in a generic point of the space, r, at a given time, t, generated by all the particles that compose the atmosphere is given by (e.g., Goldstein 1950) where r ′ is the position of an atmosphere mass element dM.For a frame centered in the planet, we can expand Eq. ( 1) in Legendre polynomials, P ℓ (x), as with where r = r/r is the unit vector and r = ||r|| is the norm, and we assume r > r ′ .For ℓ = 0 and ℓ = 1, we have P 0 (r • r′ ) = 1 and P 1 (r • r′ ) = r • r′ , respectively, and thus we can rewrite where is the center of mass of the atmosphere.Furthermore, V 0 is the potential of a spherical symmetric atmosphere, which can be absorbed in the total potential of a point-mass planet.For simplicity, we can set r cm = 0, and thus V 1 = 0.The resulting potential of a perturbed atmosphere, also known as tidal potential, is then given by the differential potential Moreover, if we assume r ≫ r ′ , we can keep only the term in ℓ = 2 (quadrupolar approximation), and finally get

Deformation of the atmosphere
For a frame attached to the planet, we can express r = (r, θ, ϕ) and r ′ = (r ′ , θ ′ , ϕ ′ ) in spherical coordinates.Then where ρ a (r ′ , t) is the local density of the atmosphere, and is the differential solid angle.Assuming a constant radius for the planet, R, we have with For terrestrial planets, the height of the atmosphere, H, is usually negligible compared to the radius of the planet1 , that is Thus, we can take a thin layer approximation and get We further assume that the self-gravity fluctuations in the atmosphere can also be neglected (e.g., Cowling 1941).Then, we use the hydrostatic equilibrium equation, where p(r ′ , t) is the local pressure and g(r ′ ) is the local gravity, to eliminate the integral in height and where p s (θ ′ , ϕ ′ , t) = p(R, θ ′ , ϕ ′ , t) is the surface pressure.We also assume that p(R + H, θ ′ , ϕ ′ , t) = 0, and that g(r ′ ) = g(R) = g is constant.In a very general formulation, the surface pressure depends on the amount of energy per unit area received from the star, that is, the insolation with F being an operator that depends on the composition and the physical properties of the atmosphere.The insolation W can be seen as a forcing function.For a measurement at the top of the atmosphere, the insolation is given by (e.g., Ward 1974) where S ⋆ = L ⋆ /(4πa 2 ) is the "solar" constant, L ⋆ is the luminosity of the star, r ⋆ is the position of the star with respect to the center of mass of the planet (which is a function of the time), a is the semimajor axis of the star-planet orbit, and We can expand expression (18) in Legendre polynomials, as with and The insolation variations can also be expressed in the frequency domain by performing a Fourier transform with Since the Legendre polynomials form an orthogonal basis, any other quantity with the same symmetry as the insolation, such as the surface pressure variations (Eq.( 17)), can also be expanded in a similar way.Therefore, we write where pn are the coefficients of order n of the decomposition of the surface pressure in Legendre polynomials.They can also be expressed in terms of spherical harmonics Y m ℓ and their complex conjugate Ȳm ℓ (e.g., Abramowitz & Stegun 1972) Then, by replacing ( 25) and (26) in expression ( 16) we get and, using the orthogonality formula 2π where δ n,2 is the Kronecker delta, we finally get The mass distribution inside the planet is usually better characterized by its inertia tensor.Then, we can rearrange Eq. (29) as with and where T denotes the transpose, I is the identity matrix, I(t) is the atmosphere inertia tensor that accounts for the departure of the mass distribution from a sphere (Eq.( 57)), while Λ(t) is a forcing tensor that depends on time through r ⋆ .The subscript B is there to remind us that Eqs. ( 30)−(33) were obtained in the body frame, that is, in a frame attached to the planet.

Deformation in an arbitrary frame
It is also possible to compute the deformation of the atmosphere in an arbitrary frame, I(t).Following Correia & Valente (2022), we let S be the rotation matrix that allows us to convert any vector r B in a frame attached to the planet into another frame r, such that r = S r B , and Then, the atmospheric tidal potential (Eq. ( 30)) becomes Similarly, we can also compute the forcing tensor (Eq.( 33)) in the new frame as We note that the expression of I B (t) (Eq.( 31)) also corresponds to the convolution product between p2 (t) and Λ B (t), where Then, in the frequency domain we have .Reference frame (p, q, s) and angles.We note that k is a unit vector that is normal to the orbital plane, s is a unit vector that is normal to the equatorial plane of the planet, and p is a unit vector along the line of nodes of the two planes.The obliquity ε is the angle between the two planes, ϖ is the argument of the pericenter, and υ is the true anomaly.
(39) Since S(t) is a rotation matrix, we note that S(t ′ +τ) = S(τ)S(t ′ ), and thus with At this stage, we need to adopt a specific frame to proceed, in order to express S(τ).We adopt here a Cartesian reference frame ( p, q, s), such that s is the spin axis, where k is a unit vector that is normal to the orbital plane of the star, p is aligned with the line of nodes between the equator of the planet and the orbital plane, and ε is the angle between these two planes, also known as the obliquity (Fig. 1).This particular frame is very useful to obtain the secular equations of motion (Sect.4.2), as it allows us to decouple the rotational perturbations from the orbital ones (Correia 2006).For a planet rotating about the s axis with angular velocity ω, we thus have The result for the product S(τ) Λ(σ) S(τ) T is given in Table 1.
For instance, for the Î23 (σ) coefficient of the inertia tensor Î(σ) (Eq.( 57)), we then finally have Similarly, for the remaining coefficients of Î(σ), we obtain with and

Surface pressure variations
The deformation of the atmosphere is given by the inertia tensor, Î(σ), which is obtained from the second harmonic of the surface pressure, p2 (σ), combined with the forcing tensor, Λ(σ) (Eq.( 40)).While the forcing tensor is well determined, as it only depends on the position of the star (Eq.( 36)), the surface pressure is subject to many uncertainties, as it depends on the composition and the physical properties of the atmosphere.
In order to compute p2 (σ), we need to adopt some dynamical model for the atmosphere.The study of thermal tides has been initiated by the pioneer work of Siebert (1961) and Chapman & Lindzen (1970).They assumed that the atmosphere of the planet corresponds to an ideal gas obeying the perfect gas law, together with the conservation of mass and the Navier-Stokes Table 1.Coefficients of the product matrix S(τ) Λ(σ) S(τ) T (Eq.( 40)).
equation.These equations were then linearized around the equilibrium state.Chapman & Lindzen (1970) studied the case of a fast rotating planet, such as the Earth, but for long-term evolution studies, the planet may encounter slow rotation regimes near synchronization, such as Venus (e.g., Correia & Laskar 2001, 2003b).In the later case, it is important to take into account the effect of radiative losses (Auclair-Desrotour et al. 2017).Assuming a slowly rotating convective atmosphere, the hydrodynamic equations can be simplified drastically, since the Coriolis acceleration can be neglected.The thin atmospheric layer close to the surface has a strongly negative temperature gradient (Seiff et al. 1980), which subjects this layer to convective instability (Baker et al. 2000).As a result, gravity waves cannot propagate in the unstable region above the surface and we can set the Brunt-Väisälä frequency to approximately zero, which corresponds to the adiabatic temperature gradient.Following Dobrovolskis & Ingersoll (1980), thermal tides are thus solely considered to be generated by the average heat at the ground, Ĵ0 , which finally gives the following for the second harmonic of the surface pressure variations (Auclair-Desrotour et al. 2017, Eq. ( 166)): with κ = 1 − γ −1 , where γ is the adiabatic exponent, ρ 0 = ρ a (0) is the mean surface density of the atmosphere, and σ 0 is the radiative frequency of the atmosphere, which depends on its thermal capacity.This theoretical estimation agrees incredibly well with empirical estimations derived from generic global climate circulation models (Leconte et al. 2015, Fig. 1).Interestingly, it is also similar to the Love number of a Maxwell tidal model (e.g., Correia et al. 2014;Auclair-Desrotour et al. 2019).The minus sign in the expression of p2 (σ) causes the pressure variations to lead the star, a phenomenon well documented for the Earth (e.g., Chapman & Lindzen 1970).The p2 (σ) is a complex number, whose modulus gives the amplitude of the pressure variations and the argument gives the phase lag between the substellar point and the maximal deformation of the atmosphere.It can be decomposed in its real and imaginary parts as which is very useful when we write the secular equations of motion (see Sect. 4.2), because the imaginary part characterizes the atmosphere's viscous response.We thus have It is also important to note that a(σ) is always an even function and b(σ) is always an odd function.

Two-body problem with thermal tides
We consider a system composed by a planet and a star with masses m and m ⋆ , respectively, on a Keplerian orbit.The orbital angular momentum is given by where a is the semimajor axis, e is the eccentricity, and k is the unit vector along the direction of G, which is normal to the orbit.The star is a point mass object with luminosity L ⋆ .The planet is composed by a completely rigid spherical mantle with radius R, and surrounded by a thin atmosphere that can be deformed under the action of thermal tides.The planet rotates with angular velocity ω = ω s, where s is the unit vector along the direction of the spin axis.The rotational angular momentum of the planet is given by where C is the moment of inertia of the mantle together with the moment of inertia of an unperturbed atmosphere and is the inertia tensor that accounts for the departure of the mass distribution in the atmosphere from a sphere (Eq.( 34)).In the absence of the stellar heating, we have I = 0.In general, the deformations in the atmosphere are very small with respect to the radius of the spherical mantle, R, and yield periodic changes in the moments of inertia, such that İω ≪ C ω (e.g., Frouard & Efroimsky 2018).Therefore, for simplicity, we can assume that I i j ≪ C (i, j = 1, 2, 3), and thus

Tidal force and torque
The atmospheric tidal potential (Eq. ( 35)) creates a differential gravitational field around the planet given by a The star, with mass m ⋆ and located at r = r ⋆ , interacts with this field, with a resulting tidal force with and where r⋆ = r ⋆ /r ⋆ = ( x1 , x2 , x3 ) is the unit vector.We decompose F and all the following vectorial quantities in the frame ( p, q, s) (Fig. 1).There is no loss of generality with this frame choice, because the vectors can always be expressed in another basis.However, the frame (p, q, s) greatly facilitates the computation of the inertia tensor I (Sect.3.2) since the position of the star does not depend on the planet's rotation rate: and where n = µ/a 3 is the mean motion, υ is the true anomaly, and ϖ is the argument of the pericenter (Fig. 1).
We follow the evolution of the system in an inertial frame because we can always project an inertial vector in a noninertial coordinate system.For the orbital evolution, we thus obtain The first term is responsible for the Keplerian motion, while the second term corresponds to the orbital correction introduced by the tidal force, which is responsible for the modifications in the orbit and spin of the planet.The evolution of the angular momentum vectors are computed from the tidal torque, and, due to the conservation of the total angular momentum, with

Expansion in Hansen coefficients
The coefficients of the inertia tensor Î(σ) (Eqs.( 44)−( 49)) are given in the frequency domain.In order to use them in the equations of motion (Sect.3.1), we need to return to the time domain using an inverse Fourier transform (Eq.( 39)): In general, the perturbations introduced by the forcing tensor Λ(t) are quasi-periodic (Eq.( 36)).Then, only a discrete number of frequencies exist, and we can express I(t) as a series: In the frame (p, q, s), the position of the star only depends on the orbital motion (Eq.( 63)), and so the only forcing frequencies are the orbital mean motion, n, and its harmonics.Thus, we can express Λ(t) and I(t) through the Hansen coefficients, X ℓ,m k , with In Table 2, we provide the expression of the Hansen coefficients used in this work expanded up to e 6 .For instance, for the Λ 23 coefficient (Eq.( 36)), we get from expression ( 63) where Similarly, for the Λ 13 coefficient (Eq.( 36)), we get Then, for the I 23 coefficient we finally have (Eqs.( 44) and ( 70)) For the remaining coefficients of I(t), we get (Eqs.( 45)−(49)) −k (e).The exact expression of these coefficients is given by expression (72). with The orbital and spin evolution of the planet under the action of thermal atmospheric tides is completely described by the set of equations ( 65), ( 67) and ( 75)-(80).

Secular dynamics
In general, the thermal atmospheric tides slowly modify the spin and the orbit of the planet, in a timescale much longer than the orbital and precession periods of the system.Therefore, we can average the equations of motion (section 3.1) over the mean anomaly and the argument of the pericenter, and obtain the equations of motion for the secular evolution of the system.

Averaging process
Following Correia & Valente (2022), to average the equations of motion, we first expand them in Hansen coefficients (Eq.( 71)), similarly to what we have done with the inertia tensor (Sect 3.2).For instance, for the last term in the s component of the tidal torque we have (Eq.( 68)) Then, we replace the expression of I 23 also expanded in Hansen coefficients (Eq.( 75)), and average over the mean anomaly, M, which is equivalent to retaining only the terms with k ′ = −k, and over the argument of the pericenter Finally, we decompose the surface pressure variations p2 (kn±ω) in its real and imaginary parts (Eq.( 53)), make use of their parity properties, and use the simplification This last arrangement is very useful, because we are able to combine terms in b(ω ± kn) in a single term b(ω − kn), which considerably simplifies the expression of the equations of motion.
In this particular case, we also note that there is no longer the contribution from a(ω − kn).
Another simplification is that we do not need to follow the evolution of the position vector anymore (Eq.( 65)).Indeed, as we average over the mean anomaly, M, the position of the planet in the orbit is no longer defined, and as we average over ϖ, the position of the pericenter is also no longer defined.Therefore, in the secular case, the equations of motion can be given by the evolution of the orbital angular momentum (Eq.( 55)) together with the evolution of the eccentricity.Alternatively, we prefer to use the evolution of the orbital energy, because it can be directly obtained from the power of the tidal force (Eq.( 60)) as and thus provides a simpler expression, from which we can later easily derive the evolution of the eccentricity.

Tidal torque and power
The orbital and spin evolution of the planet are completely described by the evolution of the angular momentum vectors, G and L, given by the torque (Eqs.( 66) and ( 67)), and by the evolution of the orbital energy, E, given by the power (Eqs.( 85)).
The secular torque (Eq.( 68)) can be written as For zero obliquity (ε = 0), the vectors p and q of the basis are not defined (Eq.( 42)).There are no singularities in this problem because both T p , T q ∝ sin ε; however, to avoid numerical issues, it is easier to express the torque simply using the unit vectors of the angular momentum quantities, k and s, as where with and where ρ = 3g/(4πGR) is the mean density of the planet.These expressions were obtained using the algebraic manipulators Maxima (2022) and TRIP (Gastineau & Laskar 2011).
Finally, for the secular power (Eq.( 85)) we get

Orbital and spin evolution
The set of equations ( 87) and ( 93) allows us to track the secular evolution of the system using the angular momentum vectors and the orbital energy.For a more intuitive description of the orbital and spin evolution, we can relate these quantities to the orbital elements and the rotation of the planet.The semimajor axis is directly given from the orbital energy (Eq.( 84)) the eccentricity from the orbital angular momentum (Eq.( 55)) and the rotation rate from the rotational angular momentum (Eq.( 58)), The angle between the orbital and equatorial planes (also known as obliquity or inclination), can be obtained from both angular momentum vectors as (Fig. 1) For a better comparison with previous studies, we can also directly obtain the evolution of all these quantities.The semimajor axis evolution is given from expression (93), while the eccentricity evolution can be computed from expressions (87), ( 95) and (98) as that is, using expressions ( 89), (90), and ( 93), For the rotation rate evolution, we have the following from Eqs. ( 67) and ( 96): that is, using expressions ( 89) and ( 90) The obliquity (or inclination) evolution is given from expressions (66), ( 67) and ( 97) In general, for a planet around a star, we have |L| ≪ |G|, and so the evolution is dominated by the first term.Finally, we can also obtain the evolution of the precession angles, that is, the angular velocity of the longitude of the node, Ω, and the precession speed of the spin axis, ψ.The line of nodes is aligned with the vector p (Fig. 1) and thus

Energy balance
The average total energy transferred to the planet due to thermal atmospheric tides is given by where E is the orbital energy (Eq.( 84)), is the rotational energy, and where ω is given by expression (101).
The total energy is then obtained by combining expressions ( 93) and ( 102) as We note that regardless of the orbital and spin parameters, the total energy transferred is always positive since b(σ) is an odd function and b(|σ|) < 0 (Eq.( 54)).

Expansion up to e 2
Most Earth-mass planets are observed in multiplanet systems, whose eccentricities are relatively small for stability reasons.Therefore, to simplify the equations of motion, we can truncate the series expansion in Hansen coefficients, and retain only terms in e 2 or smaller.

Tidal torque and power
For the average tidal torque (Eq.( 87)), we have and for the average power (Eq.( 93)) (113)

Orbital and spin evolution
The semimajor axis evolution can be directly obtained replacing Eq. ( 113) in Eq. ( 98).For the eccentricity, we have (Eq.( 100)) The evolution of the obliquity is directly obtained by replacing Eq. ( 110) in Eq. ( 103), and for the rotation rate (Eq.( 102)), This last expression for e = 0 is equivalent to Eq. ( 11) in Dobrovolskis (1980) and Eq. ( 7) in Correia & Laskar (2003a).It should also match Eq. ( 19) in Cunha et al. (2015), but they performed a development of the atmospheric tidal potential in r −5 ⋆ , while here the potential only depends on r −3 ⋆ (Eq.( 35)), and the extra factor in r −2 ⋆ is included in the forcing function (Eq.( 36)).

Energy balance
The total energy is obtained from expression (109) as (116)

Planar case
The final outcome of tidal evolution is the alignment of the spin axis with the normal to the orbit (see Correia et al. 2003).Therefore, to simplify the equations of motion, many works assume that this alignment is always present, that is, the motion is planar (ε = 0).In this case, we have x = 1 and s = k (Fig. 1).

Tidal torque and power
Using the simplification s = k yields (Eqs.( 86) and ( 87)) for the average tidal torque, with For the average power, we get the following from Eq. ( 93 The evolution of the obliquity is simply given by ε = 0 (Eq.( 103)) since sin ε = 0, that is, the motion remains planar.

Energy balance
The total energy transferred to the planet due to tides is obtained from expression ( 109

Conclusion
In this paper, we have revisited the spin and orbital evolution of a planet with a dense atmosphere under the action of thermal tides.We derive the secular equations of motion using a vectorial formalism, where the basis only depends on the unit vectors of the spin and orbital angular momenta.These vectors are related to the spin and orbital quantities, thus they are easy to obtain and independent of the chosen frame.The equations only depend on series of Hansen coefficients, which are widely used in celestial mechanics.They are obtained after averaging over the mean anomaly and over the argument of the pericenter, because thermal tides are expected to modify the orbital elements on a timescale much longer than the evolution of these two angles.
In some problems, where the pericenter evolves slowly, we can also perform a single average over the mean anomaly following the method presented in Sect. 3 of Correia & Valente (2022) for bodily tides.
The expression of the second harmonic of the surface pressure variations (Eq.( 52)) was obtained following the approximations in the work by Auclair-Desrotour et al. (2017), which is in very good agreement with the results obtained with general circulation models (Leconte et al. 2015;Auclair-Desrotour et al. 2019).Nevertheless, for more refined atmospheric models (e.g., Wu et al. 2023), we also expect that we only need to correct the expression for the surface pressure variations (Eq.( 52)) and that the dynamical equations derived in Sect. 4 do not change.
The vectorial formalism presented in this paper is well suited to study the long-term evolution of Earth-mass rocky planets.In addition to thermal atmospheric tides, we generally need to consider the bodily tides (e.g., Correia & Valente 2022), rotational deformation and general relativity corrections (e.g., Correia et al. 2011).For multiplanet systems, the secular interactions can be obtained either by developing the perturbing functions in terms of Laplace coefficients, suited for nonresonant compact systems (e.g., Boué & Fabrycky 2014), or in terms of Legendre polynomials, suited for hierarchical systems (e.g., Correia et al. 2016).All these previous studies adopted the same kind of reference vectors and thus each individual effect can be simply added to get the complete dynamical evolution of the planet.
Fig.1.Reference frame (p, q, s) and angles.We note that k is a unit vector that is normal to the orbital plane, s is a unit vector that is normal to the equatorial plane of the planet, and p is a unit vector along the line of nodes of the two planes.The obliquity ε is the angle between the two planes, ϖ is the argument of the pericenter, and υ is the true anomaly.

Table 2 .
Hansen coefficients up to e 6 .