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

Ema F. S. Valente, Alexandre C. M. Correia

► To cite this version:

Ema F. S. Valente, Alexandre C. M. Correia. Spin and orbital dynamics of planets undergoing thermal atmospheric tides using a vectorial approach. Astronomy & Astrophysics - A&A, 2023, 679, 10.1051/0004-6361/202245443. insu-04849151

HAL Id: insu-04849151

https://insu.hal.science/insu-04849151v1

Submitted on 19 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Spin and orbital dynamics of planets undergoing thermal atmospheric tides using a vectorial approach

Ema F. S. Valente<fnref n="1" /> and Alexandre C. M. Correia<fnref n="1" />,<fnref n="2" />

Received 11 November 2022 / Accepted 6 March 2023

## ABSTRACT

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.

**Key words.** planet-star interactions – planets and satellites: atmospheres – planets and satellites: dynamical evolution and stability – planets and satellites: terrestrial planets – celestial mechanics

## 1. 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 K-dwarf 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

<footnote n="1">¹ CFisUC, Departamento de Física, Universidade de Coimbra, 3004-516 Coimbra, Portugal
e-mail: acor@uc.pt</footnote>

<footnote n="2">² IMCCE, Observatoire de Paris, PSL Université, 77 Av. Denfert-Rochereau, 75014 Paris, France</footnote>
provide simpler expressions that are easy to implement and suitable for long-term evolution studies. In Sects. 5 and 6, we provide the equations of motion in the simplified cases of small eccentricity, and planar motion, respectively. Finally, the last section is devoted to the conclusions.

## 2. 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.

### 2.1. Atmospheric tidal potential

The gravitational potential in a generic point of the space, $\mathbf{r}$, at a given time, $t$, generated by all the particles that compose the atmosphere is given by (e.g., Goldstein 1950)

$$V(\mathbf{r}, t) = -G \int_M \frac{dM}{|\mathbf{r} - \mathbf{r}'|}, \qquad (1)$$

where $\mathbf{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_l(x)$, as

$$V(\mathbf{r}, t) = \sum_{\ell=0}^{\infty} V_{\ell}(\mathbf{r}, t), \qquad (2)$$

with

$$V_{\ell}(\mathbf{r}, t) = -\frac{G}{r} \int_{M} \left(\frac{r'}{r}\right)^{\ell} P_{\ell}(\hat{\mathbf{r}} \cdot \hat{\mathbf{r}}') dM, \qquad (3)$$

where $\hat{\mathbf{r}} = \mathbf{r}/r$ is the unit vector and $r = \|\mathbf{r}\|$ is the norm, and we assume $r > r'$. For $\ell = 0$ and $\ell = 1$, we have $P_0(\hat{\mathbf{r}} \cdot \hat{\mathbf{r}}') = 1$ and $P_1(\hat{\mathbf{r}} \cdot \hat{\mathbf{r}}') = \hat{\mathbf{r}} \cdot \hat{\mathbf{r}}'$, respectively, and thus we can rewrite

$$V_0(\mathbf{r}, t) = -\frac{GM}{r} \qquad (4)$$

and

$$V_1(\mathbf{r}, t) = -\frac{GM}{r^2} \hat{\mathbf{r}} \cdot \mathbf{r}_{\text{cm}}, \qquad (5)$$

where

$$\mathbf{r}_{\text{cm}} = \frac{1}{M} \int_{M} \mathbf{r}' dM, \qquad (6)$$

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 $\mathbf{r}_{\text{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

$$\Delta V(\mathbf{r}, t) = V(\mathbf{r}, t) - V_0(\mathbf{r}, t) - V_1(\mathbf{r}, t) = \sum_{\ell=2}^{\infty} V_{\ell}(\mathbf{r}, t). \qquad (7)$$

Moreover, if we assume $r \gg r'$, we can keep only the term in $\ell = 2$ (quadrupolar approximation), and finally get

$$\Delta V(\mathbf{r}, t) \approx V_2(\mathbf{r}, t) = -\frac{G}{r^3} \int_M r'^2 P_2(\hat{\mathbf{r}} \cdot \hat{\mathbf{r}}') dM. \qquad (8)$$

### 2.2. Deformation of the atmosphere

For a frame attached to the planet, we can express $\mathbf{r} = (r, \theta, \phi)$ and $\mathbf{r}' = (r', \theta', \phi')$ in spherical coordinates. Then

$$dM = \rho_a(\mathbf{r}', t) r'^2 dr'd\Omega', \qquad (9)$$

where $\rho_a(\mathbf{r}', t)$ is the local density of the atmosphere, and

$$d\Omega' = \sin\theta' d\theta' d\phi' \qquad (10)$$

is the differential solid angle. Assuming a constant radius for the planet, $R$, we have

$$\Delta V(\mathbf{r}, t) = -\frac{G}{r^3} \int_{r'=R}^{+\infty} \int_{\phi'=0}^{2\pi} \int_{\theta'=0}^{\pi} \rho_a(\mathbf{r}', t) r'^4 P_2(\hat{\mathbf{r}} \cdot \hat{\mathbf{r}}') dr' d\Omega', \qquad (11)$$

with

$$\hat{\mathbf{r}} \cdot \hat{\mathbf{r}}' = \cos\theta \cos\theta' + \sin\theta \sin\theta' \cos(\phi - \phi'). \qquad (12)$$

For terrestrial planets, the height of the atmosphere, $H$, is usually negligible compared to the radius of the planet<fnref n="1" />, that is

$$R \le r' \le R + H, \quad \text{and} \quad \frac{H}{R} \ll 1. \qquad (13)$$

Thus, we can take a thin layer approximation and get

$$\Delta V(\mathbf{r}, t) = -GR \left(\frac{R}{r}\right)^3 \int_{r'=R}^{R+H} \int_{\phi'=0}^{2\pi} \int_{\theta'=0}^{\pi} \rho_a(\mathbf{r}', t) P_2(\hat{\mathbf{r}} \cdot \hat{\mathbf{r}}') dr' d\Omega'. \qquad (14)$$

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,

$$dp = -\rho_a(\mathbf{r}', t)g(\mathbf{r}') dr', \qquad (15)$$

where $p(\mathbf{r}', t)$ is the local pressure and $g(\mathbf{r}')$ is the local gravity, to eliminate the integral in height

$$\Delta V(\mathbf{r}, t) = -\frac{GR}{g} \left(\frac{R}{r}\right)^3 \int_{\phi'=0}^{2\pi} \int_{\theta'=0}^{\pi} p_s(\theta', \phi', t) P_2(\hat{\mathbf{r}} \cdot \hat{\mathbf{r}}') d\Omega', \qquad (16)$$

and where $p_s(\theta', \phi', t) = p(R, \theta', \phi', t)$ is the surface pressure. We also assume that $p(R+H, \theta', \phi', t) = 0$, and that $g(\mathbf{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 $W(\theta', \phi', t)$,

$$p_s(\theta', \phi', t) = F[W(\theta', \phi', t)], \qquad (17)$$

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)

$$W(\theta', \phi', t) = S_{\star} \left(\frac{a}{r_{\star}}\right)^2 \begin{cases} \hat{\mathbf{r}}' \cdot \hat{\mathbf{r}}_{\star}, & \hat{\mathbf{r}}' \cdot \hat{\mathbf{r}}_{\star} > 0 \\ 0, & \hat{\mathbf{r}}' \cdot \hat{\mathbf{r}}_{\star} \le 0 \end{cases}, \qquad (18)$$

<footnote n="1">¹ More generally, we have $p(R+z) \approx p(R)$ exp($-z/H_0$), where $H_0$ is the “scale height”. In the case of the Earth, $H_0 \sim 8$ km, and because the pressure decreases exponentially, 95% of the mass of the gas is contained in $H \sim 3H_0 \sim 24$ km $\ll R_{\oplus}$ (e.g., Chapman & Lindzen 1970).</footnote>
where $S_\star = L_\star/(4\pi a^2)$ is the “solar” constant, $L_\star$ is the luminosity of the star, $\mathbf{r}_\star$ 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

$$ \hat{\mathbf{r}}' \cdot \hat{\mathbf{r}}_\star = \cos \theta' \cos \theta_\star + \sin \theta' \sin \theta_\star \cos(\phi' - \phi_\star). \quad (19) $$

We can expand expression (18) in Legendre polynomials, as

$$ W(\theta', \phi', t) = \sum_{n=0}^{\infty} W_n \Theta_n(\mathbf{r}', t), \quad (20) $$

with

$$ \Theta_n(\mathbf{r}, t) = \left(\frac{a}{r_\star}\right)^2 P_n(\hat{\mathbf{r}} \cdot \hat{\mathbf{r}}_\star), \quad (21) $$

and

$$ W_n = S_\star \frac{2n+1}{2} \int_0^1 x P_n(x) \, dx. \quad (22) $$

The insolation variations can also be expressed in the frequency domain by performing a Fourier transform

$$ W(\theta', \phi', t) = \sum_{n=0}^{\infty} \int W_n \hat{\Theta}_n(\mathbf{r}', \sigma) e^{i\sigma t} d\sigma, \quad (23) $$

with

$$ \hat{\Theta}_n(\mathbf{r}, \sigma) = \int \Theta_n(\mathbf{r}, t) e^{-i\sigma t} dt. \quad (24) $$

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

$$ p_s(\theta', \phi', t) = \sum_{n=0}^{\infty} \int \hat{p}_n(\sigma) \hat{\Theta}_n(\mathbf{r}', \sigma) e^{i\sigma t} d\sigma, \quad (25) $$

where $\hat{p}_n$ 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_\ell^m$ and their complex conjugate $\bar{Y}_\ell^m$ (e.g., Abramowitz & Stegun 1972)

$$ P_\ell(\hat{\mathbf{r}} \cdot \hat{\mathbf{r}}') = \frac{4\pi}{2\ell+1} \sum_{m=-\ell}^{\ell} Y_\ell^m(\theta, \phi) \bar{Y}_\ell^m(\theta', \phi'). \quad (26) $$

Then, by replacing (25) and (26) in expression (16) we get

$$ \Delta V(\mathbf{r}, t) = -\frac{GR}{g} \left(\frac{R}{r}\right)^3 \frac{4\pi}{5} \sum_{m=-2}^{2} Y_2^m(\theta, \phi) \sum_{n=0}^{\infty} \int \hat{p}_n(\sigma) \\ \times \int_{\phi'=0}^{2\pi} \int_{\theta'=0}^{\pi} \hat{\Theta}_n(\mathbf{r}', \sigma) \bar{Y}_2^m(\theta', \phi') \, d\Omega' e^{i\sigma t} \, d\sigma, \quad (27) $$

and, using the orthogonality formula

$$ \int_{\phi'=0}^{2\pi} \int_{\theta'=0}^{\pi} P_n(\hat{\mathbf{r}}' \cdot \hat{\mathbf{r}}_\star) \bar{Y}_2^m(\theta', \phi') \, d\Omega' = \frac{4\pi}{5} \bar{Y}_2^m(\theta_\star, \phi_\star) \delta_{n,2}, \quad (28) $$

where $\delta_{n,2}$ is the Kronecker delta, we finally get

$$ \Delta V(\mathbf{r}, t) = -\frac{GR}{g} \left(\frac{R}{r}\right)^3 \frac{4\pi}{5} \int \hat{p}_2(\sigma) \hat{\Theta}_2(\mathbf{r}, \sigma) e^{i\sigma t} d\sigma. \quad (29) $$

The mass distribution inside the planet is usually better characterized by its inertia tensor. Then, we can rearrange Eq. (29) as

$$ \Delta V(\mathbf{r}_B, t) = \frac{3G}{2r^3} \hat{\mathbf{r}}_B \cdot \mathcal{I}_B(t) \cdot \hat{\mathbf{r}}_B, \quad (30) $$

with

$$ \mathcal{I}_B(t) = \int \hat{p}_2(\sigma) \hat{\Lambda}_B(\sigma) e^{i\sigma t} d\sigma, \quad (31) $$

$$ \hat{\Lambda}_B(\sigma) = \int \Lambda_B(t) e^{-i\sigma t} dt, \quad (32) $$

and

$$ \Lambda_B(t) = -\frac{4\pi R^4}{5g} \left(\frac{a}{r_\star}\right)^2 \left(\hat{\mathbf{r}}_{\star B} \hat{\mathbf{r}}_{\star B}^T - \frac{1}{3} \mathbb{I}\right), \quad (33) $$

where $^T$ denotes the transpose, $\mathbb{I}$ is the identity matrix, $\mathcal{I}(t)$ is the atmosphere inertia tensor that accounts for the departure of the mass distribution from a sphere (Eq. (56)), while $\Lambda(t)$ is a forcing tensor that depends on time through $\mathbf{r}_\star$. 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.

## 2.3. Deformation in an arbitrary frame

It is also possible to compute the deformation of the atmosphere in an arbitrary frame, $\mathcal{I}(t)$. Following Correia & Valente (2022), we let $S$ be the rotation matrix that allows us to convert any vector $\mathbf{r}_B$ in a frame attached to the planet into another frame $\mathbf{r}$, such that $\mathbf{r} = S\mathbf{r}_B$, and

$$ \mathcal{I}(t) = S(t) \mathcal{I}_B(t) S(t)^T. \quad (34) $$

Then, the atmospheric tidal potential (Eq. (30)) becomes

$$ \Delta V(\mathbf{r}, t) = \frac{3G}{2r^3} \hat{\mathbf{r}} \cdot \mathcal{I}(t) \cdot \hat{\mathbf{r}}. \quad (35) $$

Similarly, we can also compute the forcing tensor (Eq. (33)) in the new frame as

$$ \Lambda(t) = S(t) \Lambda_B(t) S(t)^T = -\frac{4\pi R^4}{5g} \left(\frac{a}{r_\star}\right)^2 \left(\hat{\mathbf{r}}_\star \hat{\mathbf{r}}_\star^T - \frac{1}{3} \mathbb{I}\right). \quad (36) $$

We note that the expression of $\mathcal{I}_B(t)$ (Eq. (31)) also corresponds to the convolution product between $\tilde{p}_2(t)$ and $\Lambda_B(t)$,

$$ \mathcal{I}_B(t) = [\tilde{p}_2 * \Lambda_B](t) = \int \tilde{p}_2(t - t') \Lambda_B(t') dt', \quad (37) $$

where

$$ \tilde{p}_2(t) = \int \hat{p}_2(\sigma) e^{i\sigma t} d\sigma. \quad (38) $$
<table><caption>**Table 1.** Coefficients of the product matrix $S(\tau) \hat{\Lambda}(\sigma) S(\tau)^T$ (Eq. (40)).</caption><tr><td>-ΔΛ(σ) cos 2ωτ - Λ̂₁₂(σ) sin 2ωτ - ½Λ̂₃₃(σ)</td><td>-ΔΛ(σ) sin 2ωτ + Λ̂₁₂(σ) cos 2ωτ</td><td>Λ̂₁₃(σ) cos ωτ - Λ̂₂₃(σ) sin ωτ</td></tr><tr><td>-ΔΛ(σ) sin 2ωτ + Λ̂₁₂(σ) cos 2ωτ</td><td>ΔΛ(σ) cos 2ωτ + Λ̂₁₂(σ) sin 2ωτ - ½Λ̂₃₃(σ)</td><td>Λ̂₁₃(σ) sin ωτ + Λ̂₂₃(σ) cos ωτ</td></tr><tr><td>Λ̂₁₃(σ) cos ωτ - Λ̂₂₃(σ) sin ωτ</td><td>Λ̂₁₃(σ) sin ωτ + Λ̂₂₃(σ) cos ωτ</td><td>Λ̂₃₃(σ)</td></tr></table>

Then, in the frequency domain we have

$$ 
\begin{aligned}
\hat{I}(\sigma) &= \int \mathcal{I}(t) e^{-i\sigma t} dt = \int S(t) \mathcal{I}_B(t) S(t)^T e^{-i\sigma t} dt \\
&= \iint \tilde{p}_2(t-t') S(t) \Lambda_B(t') S(t)^T e^{-i\sigma t} dt'dt \quad (39) \\
&= \iint \tilde{p}_2(\tau) S(t'+\tau) \Lambda_B(t') S(t'+\tau)^T e^{-i\sigma(t'+\tau)} dt'd\tau.
\end{aligned}
$$

Since $S(t)$ is a rotation matrix, we note that $S(t' + \tau) = S(\tau)S(t')$, and thus

$$ 
\begin{aligned}
\hat{I}(\sigma) &= \iint \tilde{p}_2(\tau) S(\tau) \Lambda(t') S(\tau)^T e^{-i\sigma(t'+\tau)} dt'd\tau \quad (40) \\
&= \int \tilde{p}_2(\tau) S(\tau) \hat{\Lambda}(\sigma) S(\tau)^T e^{-i\sigma\tau} d\tau,
\end{aligned}
$$

with

$$ \hat{\Lambda}(\sigma) = \int \Lambda(t) e^{-i\sigma t} dt. \quad (41) $$

At this stage, we need to adopt a specific frame to proceed, in order to express $S(\tau)$. We adopt here a Cartesian reference frame $(\mathbf{p}, \mathbf{q}, \mathbf{s})$, such that $\mathbf{s}$ is the spin axis,

$$ \mathbf{p} = \frac{k \times s}{\sin \varepsilon}, \quad \mathbf{q} = \frac{k - \cos \varepsilon s}{\sin \varepsilon}, \quad \cos \varepsilon = \mathbf{k} \cdot \mathbf{s}, \quad (42) $$

where $\mathbf{k}$ is a unit vector that is normal to the orbital plane of the star, $\mathbf{p}$ is aligned with the line of nodes between the equator of the planet and the orbital plane, and $\varepsilon$ 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 $\mathbf{s}$ axis with angular velocity $\omega$, we thus have

$$ S(\tau) = \begin{bmatrix} \cos \omega\tau & -\sin \omega\tau & 0 \\ \sin \omega\tau & \cos \omega\tau & 0 \\ 0 & 0 & 1 \end{bmatrix}. \quad (43) $$

The result for the product $S(\tau) \hat{\Lambda}(\sigma) S(\tau)^T$ is given in Table 1. For instance, for the $\hat{I}_{23}(\sigma)$ coefficient of the inertia tensor $\hat{I}(\sigma)$ (Eq. (57)), we then finally have

$$ 
\begin{aligned}
\hat{I}_{23}(\sigma) &= \int \tilde{p}_2(\tau) [\hat{\Lambda}_{13}(\sigma) \sin \omega\tau + \hat{\Lambda}_{23}(\sigma) \cos \omega\tau] e^{-i\sigma\tau} d\tau \\
&= \hat{\Lambda}_{13}(\sigma) \int \tilde{p}_2(\tau) \sin \omega\tau e^{-i\sigma\tau} d\tau \\
&\quad + \hat{\Lambda}_{23}(\sigma) \int \tilde{p}_2(\tau) \cos \omega\tau e^{-i\sigma\tau} d\tau \quad (44) \\
&= \frac{1}{2} \hat{p}_2(\sigma - \omega) [\hat{\Lambda}_{23}(\sigma) - i \hat{\Lambda}_{13}(\sigma)] \\
&\quad + \frac{1}{2} \hat{p}_2(\sigma + \omega) [\hat{\Lambda}_{23}(\sigma) + i \hat{\Lambda}_{13}(\sigma)].
\end{aligned}
$$

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Reference frame $(\mathbf{p}, \mathbf{q}, \mathbf{s})$ and angles. We note that $\mathbf{k}$ is a unit vector that is normal to the orbital plane, $\mathbf{s}$ is a unit vector that is normal to the equatorial plane of the planet, and $\mathbf{p}$ is a unit vector along the line of nodes of the two planes. The obliquity $\varepsilon$ is the angle between the two planes, $\omega$ is the argument of the pericenter, and $\nu$ is the true anomaly.</figcaption></figure>

Similarly, for the remaining coefficients of $\hat{I}(\sigma)$, we obtain

$$ 
\begin{aligned}
\hat{I}_{13}(\sigma) &= \frac{1}{2} \hat{p}_2(\sigma - \omega) [\hat{\Lambda}_{13}(\sigma) + i \hat{\Lambda}_{23}(\sigma)] \\
&\quad + \frac{1}{2} \hat{p}_2(\sigma + \omega) [\hat{\Lambda}_{13}(\sigma) - i \hat{\Lambda}_{23}(\sigma)],
\end{aligned}
\quad (45) $$

$$ 
\begin{aligned}
\hat{I}_{12}(\sigma) &= \frac{1}{2} \hat{p}_2(\sigma - 2\omega) [\hat{\Lambda}_{12}(\sigma) + i \Delta\hat{\Lambda}(\sigma)] \\
&\quad + \frac{1}{2} \hat{p}_2(\sigma + 2\omega) [\hat{\Lambda}_{12}(\sigma) - i \Delta\hat{\Lambda}(\sigma)],
\end{aligned}
\quad (46) $$

$$ \hat{I}_{33}(\sigma) = \hat{p}_2(\sigma)\hat{\Lambda}_{33}(\sigma), \quad (47) $$

$$ \hat{I}_{22}(\sigma) = \Delta\hat{I}(\sigma) - \frac{1}{2}\hat{I}_{33}(\sigma), \quad (48) $$

$$ \hat{I}_{11}(\sigma) = -\Delta\hat{I}(\sigma) - \frac{1}{2}\hat{I}_{33}(\sigma) \quad (49) $$

with

$$ 
\begin{aligned}
\Delta\hat{I}(\sigma) &= \frac{1}{2} \hat{p}_2(\sigma - 2\omega) [\Delta\hat{\Lambda}(\sigma) - i \hat{\Lambda}_{12}(\sigma)] \\
&\quad + \frac{1}{2} \hat{p}_2(\sigma + 2\omega) [\Delta\hat{\Lambda}(\sigma) + i \hat{\Lambda}_{12}(\sigma)],
\end{aligned}
\quad (50) $$

and

$$ \Delta\hat{\Lambda}(\sigma) = \frac{1}{2} [\hat{\Lambda}_{22}(\sigma) - \hat{\Lambda}_{11}(\sigma)]. \quad (51) $$
## 2.4. Surface pressure variations

The deformation of the atmosphere is given by the inertia tensor, $\hat{\mathcal{I}}(\sigma)$, which is obtained from the second harmonic of the surface pressure, $\hat{p}_2(\sigma)$, combined with the forcing tensor, $\hat{\Lambda}(\sigma)$ (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 $\hat{p}_2(\sigma)$, 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 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, $\hat{J}_0$, which finally gives the following for the second harmonic of the surface pressure variations (Auclair-Desrotour et al. 2017, Eq. (166)):

$$ \hat{p}_2(\sigma) = -\frac{\kappa \rho_0 \hat{J}_0}{\sigma_0 + i\sigma}, \qquad (52) $$

with $\kappa = 1 - \gamma^{-1}$, where $\gamma$ is the adiabatic exponent, $\rho_0 = \rho_a(0)$ is the mean surface density of the atmosphere, and $\sigma_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 $\hat{p}_2(\sigma)$ causes the pressure variations to lead the star, a phenomenon well documented for the Earth (e.g., Chapman & Lindzen 1970).

The $\hat{p}_2(\sigma)$ 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

$$ \hat{p}_2(\sigma) = a(\sigma) - ib(\sigma), \qquad (53) $$

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

$$ a(\sigma) = -\frac{\kappa \rho_0 \hat{J}_0 \sigma_0}{\sigma_0^2 + \sigma^2}, \quad \text{and} \quad b(\sigma) = -\frac{\kappa \rho_0 \hat{J}_0 \sigma}{\sigma_0^2 + \sigma^2}. \qquad (54) $$

It is also important to note that $a(\sigma)$ is always an even function and $b(\sigma)$ is always an odd function.

## 3. Two-body problem with thermal tides

We consider a system composed by a planet and a star with masses $m$ and $m_\star$, respectively, on a Keplerian orbit. The orbital angular momentum is given by

$$ \mathcal{G} = \beta \sqrt{\mu a(1 - e^2) k}, \qquad (55) $$

where $a$ is the semimajor axis, $e$ is the eccentricity, $\beta = m_\star m/(m_\star + m)$, $\mu = \mathcal{G}(m_\star + m)$, $\mathcal{G}$ is the gravitational constant, and $k$ is the unit vector along the direction of $\mathcal{G}$, which is normal to the orbit. The star is a point mass object with luminosity $L_\star$. 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 $\omega = \omega s$, where $s$ is the unit vector along the direction of the spin axis. The rotational angular momentum of the planet is given by

$$ L = C\omega + \mathcal{I} \cdot \omega, \qquad (56) $$

where $C$ is the moment of inertia of the mantle together with the moment of inertia of an unperturbed atmosphere and

$$ \mathcal{I} = \begin{bmatrix} I_{11} & I_{12} & I_{13} \\ I_{12} & I_{22} & I_{23} \\ I_{13} & I_{23} & I_{33} \end{bmatrix} \qquad (57) $$

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 $\mathcal{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 $\mathcal{I}\omega \ll C\dot{\omega}$ (e.g., Frouard & Efroimsky 2018). Therefore, for simplicity, we can assume that $I_{ij} \ll C$ ($i, j = 1, 2, 3$), and thus

$$ L \approx C\omega = C\omega s. \qquad (58) $$

### 3.1. Tidal force and torque

The atmospheric tidal potential (Eq. (35)) creates a differential gravitational field around the planet given by

$$ a(\mathbf{r}) = -\nabla_r [\Delta V(\mathbf{r})]. \qquad (59) $$

The star, with mass $m_\star$ and located at $\mathbf{r} = \mathbf{r}_\star$, interacts with this field, with a resulting tidal force

$$ F = m_{\star} a(\mathbf{r}_{\star}) = F_{1} + F_{2}, \qquad (60) $$

with

$$ F_1 = \frac{15 G m_\star}{r_\star^4} \left[ \frac{I_{22} - I_{11}}{2} (\hat{x}_2^2 - \frac{1}{5}) + \frac{I_{33} - I_{11}}{2} (\hat{x}_3^2 - \frac{1}{5}) \right. \\ \left. + I_{12} \hat{x}_1 \hat{x}_2 + I_{13} \hat{x}_1 \hat{x}_3 + I_{23} \hat{x}_2 \hat{x}_3 \right] \hat{r}_\star \qquad (61) $$

and

$$ F_2 = -\frac{3Gm_\star}{r_\star^4} \left[ (I_{22} - I_{11})\hat{x}_2 q + (I_{33} - I_{11})\hat{x}_3 s \right. \\ \left. + I_{12}(\hat{x}_1 q + \hat{x}_2 p) + I_{13}(\hat{x}_1 s + \hat{x}_3 p) + I_{23}(\hat{x}_2 s + \hat{x}_3 q) \right], \qquad (62) $$
where $\hat{\mathbf{r}}_{\star} = \mathbf{r}_{\star}/r_{\star} = (\hat{x}_{1}, \hat{x}_{2}, \hat{x}_{3})$ is the unit vector. We decompose $\mathbf{F}$ and all the following vectorial quantities in the frame $(\mathbf{p}, \mathbf{q}, \mathbf{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 $(\mathbf{p}, \mathbf{q}, \mathbf{s})$ greatly facilitates the computation of the inertia tensor $\mathcal{I}$ (Sect. 3.2) since the position of the star does not depend on the planet's rotation rate:

$$ 
\begin{aligned}
\hat{x}_1 &= \hat{\mathbf{r}}_{\star} \cdot \mathbf{p} = \cos(\varpi + \upsilon), \\
\hat{x}_2 &= \hat{\mathbf{r}}_{\star} \cdot \mathbf{q} = \cos \varepsilon \sin(\varpi + \upsilon), \\
\hat{x}_3 &= \hat{\mathbf{r}}_{\star} \cdot \mathbf{s} = -\sin \varepsilon \sin(\varpi + \upsilon),
\end{aligned}
\qquad (63) $$

and

$$ 
\begin{aligned}
\dot{\mathbf{r}}_{\star} \cdot \mathbf{p} &= \frac{-na}{\sqrt{1-e^2}} (\sin(\varpi + \upsilon) + e \sin \varpi), \\
\dot{\mathbf{r}}_{\star} \cdot \mathbf{q} &= \frac{na}{\sqrt{1-e^2}} \cos \varepsilon (\cos(\varpi + \upsilon) + e \cos \varpi), \\
\dot{\mathbf{r}}_{\star} \cdot \mathbf{s} &= \frac{-na}{\sqrt{1-e^2}} \sin \varepsilon (\cos(\varpi + \upsilon) + e \cos \varpi),
\end{aligned}
\qquad (64) $$

where $n = \sqrt{\mu/a^3}$ is the mean motion, $\upsilon$ is the true anomaly, and $\varpi$ 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

$$ \ddot{\mathbf{r}}_{\star} = -\frac{\mu}{r_{\star}^2} \hat{\mathbf{r}}_{\star} + \frac{\mathbf{F}}{\beta}. \qquad (65) $$

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,

$$ \dot{\mathbf{G}} = \mathbf{T} = \mathbf{r}_{\star} \times \mathbf{F} = \mathbf{r}_{\star} \times \mathbf{F}_{2}, \qquad (66) $$

and, due to the conservation of the total angular momentum,

$$ \dot{\mathbf{L}} = -\dot{\mathbf{G}} = -\mathbf{T}, \qquad (67) $$

with

$$ \mathbf{T} = -\frac{3\mathbf{G}m_{\star}}{r_{\star}^3} \begin{bmatrix} (I_{33} - I_{22})\hat{x}_2\hat{x}_3 - I_{12}\hat{x}_1\hat{x}_3 + I_{13}\hat{x}_1\hat{x}_2 + I_{23}(\hat{x}_2^2 - \hat{x}_3^2) \\ (I_{11} - I_{33})\hat{x}_1\hat{x}_3 + I_{12}\hat{x}_2\hat{x}_3 + I_{13}(\hat{x}_3^2 - \hat{x}_1^2) - I_{23}\hat{x}_1\hat{x}_2 \\ (I_{22} - I_{11})\hat{x}_1\hat{x}_2 + I_{12}(\hat{x}_1^2 - \hat{x}_2^2) - I_{13}\hat{x}_2\hat{x}_3 + I_{23}\hat{x}_1\hat{x}_3 \end{bmatrix}. \qquad (68) $$

### 3.2. Expansion in Hansen coefficients

The coefficients of the inertia tensor $\hat{\mathcal{I}}(\sigma)$ (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)):

$$ \mathcal{I}(t) = \int \hat{\mathcal{I}}(\sigma) e^{i\sigma t} d\sigma. \qquad (69) $$

In general, the perturbations introduced by the forcing tensor $\Lambda(t)$ are quasi-periodic (Eq. (36)). Then, only a discrete number of frequencies exist, and we can express $\mathcal{I}(t)$ as a series:

$$ \mathcal{I}(t) = \sum_k \hat{\mathcal{I}}(\sigma_k) e^{i\sigma_k t}. \qquad (70) $$

In the frame $(\mathbf{p}, \mathbf{q}, \mathbf{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 $\Lambda(t)$ and $\mathcal{I}(t)$ through the Hansen coefficients, $X_k^{\ell,m}$,

$$ \left(\frac{r_{\star}}{a}\right)^{\ell} e^{im\upsilon} = \sum_{k=-\infty}^{+\infty} X_{k}^{\ell,m}(e) e^{ikM}, \qquad (71) $$

with

$$ 
\begin{aligned}
X_k^{\ell,m}(e) &= \frac{1}{2\pi} \int_{-\pi}^{\pi} \left(\frac{r_{\star}}{a}\right)^{\ell} e^{i(m\upsilon - kM)} dM \\
&= \frac{1}{\pi} \int_{0}^{\pi} \left(\frac{1 - e^2}{1 + e \cos \upsilon}\right)^{\ell} \cos(m\upsilon - kM) dM.
\end{aligned}
\qquad (72) $$

In Table 2, we provide the expression of the Hansen coefficients used in this work expanded up to $e^6$.

For instance, for the $\Lambda_{23}$ coefficient (Eq. (36)), we get from expression (63)

$$ 
\begin{aligned}
\Lambda_{23} &= -\frac{4\pi R^4}{5g} \left(\frac{a}{r_{\star}}\right)^2 \hat{x}_2 \hat{x}_3 = \frac{4\pi R^4}{5g} \left(\frac{a}{r_{\star}}\right)^2 \sin \varepsilon \cos \varepsilon \sin^2(\varpi + \upsilon) \\
&= \frac{\pi R^4}{5g} \sin \varepsilon \cos \varepsilon \left(\frac{a}{r_{\star}}\right)^2 (2 - e^{i2(\varpi + \upsilon)} - e^{-i2(\varpi + \upsilon)}) \\
&= \frac{\pi R^4}{5g} \sin \varepsilon \sum_{k=-\infty}^{+\infty} x \left(2X_k^{-2,0} - e^{i2\varpi} X_k^{-2,2} - e^{-i2\varpi} X_k^{-2,-2}\right) e^{ikM},
\end{aligned}
\qquad (73) $$

where $x = \cos \varepsilon = \mathbf{k} \cdot \mathbf{s}$.

Similarly, for the $\Lambda_{13}$ coefficient (Eq. (36)), we get

$$ \Lambda_{13} = i \frac{\pi R^4}{5g} \sin \varepsilon \sum_{k=-\infty}^{+\infty} \left( e^{-i2\varpi} X_k^{-2,-2} - e^{i2\varpi} X_k^{-2,2} \right) e^{ikM}. \qquad (74) $$

Then, for the $I_{23}$ coefficient we finally have (Eqs. (44) and (70))

$$ I_{23} = \frac{\pi R^4}{10g} \sin \varepsilon \sum_{k=-\infty}^{+\infty} \left\{ \hat{p}_2(kn - \omega) \left[ 2xX_k^{-2,0} - (1+x) e^{i2\varpi} X_k^{-2,2} \right. \right. \\
\left. \left. \qquad + (1-x) e^{-i2\varpi} X_k^{-2,-2} \right] \right. \\
\left. \qquad + \hat{p}_2(kn + \omega) \left[ 2xX_k^{-2,0} + (1-x) e^{i2\varpi} X_k^{-2,2} \right. \right. \\
\left. \left. \qquad - (1+x) e^{-i2\varpi} X_k^{-2,-2} \right] \right\} e^{ikM}. \qquad (75) $$

For the remaining coefficients of $\mathcal{I}(t)$, we get (Eqs. (45)–(49))

$$ I_{13} = i \frac{\pi R^4}{10g} \sin \varepsilon \sum_{k=-\infty}^{+\infty} \left\{ \hat{p}_2(kn - \omega) \left[ 2xX_k^{-2,0} - (1+x) e^{i2\varpi} X_k^{-2,2} \right. \right. \\
\left. \left. \qquad + (1-x) e^{-i2\varpi} X_k^{-2,-2} \right] \right. \\
\left. \qquad - \hat{p}_2(kn + \omega) \left[ 2xX_k^{-2,0} + (1-x) e^{i2\varpi} X_k^{-2,2} \right. \right. \\
\left. \left. \qquad - (1+x) e^{-i2\varpi} X_k^{-2,-2} \right] \right\} e^{ikM}, \qquad (76) $$
<table><caption>**Table 2.** Hansen coefficients up to $e^6$.</caption><thead><tr><th>k</th><th>X<sub>k</sub><sup>-2,0</sup>(e)</th><th>X<sub>k</sub><sup>-2,2</sup>(e)</th><th>X<sub>k</sub><sup>-3,0</sup>(e)</th><th>X<sub>k</sub><sup>-3,2</sup>(e)</th></tr></thead><tbody><tr><td>-6</td><td><sup>1223</sup>&frasl;<sub>320</sub> e<sup>6</sup></td><td></td><td><sup>3167</sup>&frasl;<sub>320</sub> e<sup>6</sup></td><td></td></tr><tr><td>-5</td><td><sup>1097</sup>&frasl;<sub>384</sub> e<sup>5</sup></td><td></td><td><sup>1773</sup>&frasl;<sub>256</sub> e<sup>5</sup></td><td></td></tr><tr><td>-4</td><td><sup>103</sup>&frasl;<sub>48</sub> e<sup>4</sup> - <sup>129</sup>&frasl;<sub>160</sub> e<sup>6</sup></td><td>-<sup>2</sup>&frasl;<sub>45</sub> e<sup>6</sup></td><td><sup>77</sup>&frasl;<sub>16</sub> e<sup>4</sup> + <sup>129</sup>&frasl;<sub>160</sub> e<sup>6</sup></td><td><sup>4</sup>&frasl;<sub>45</sub> e<sup>6</sup></td></tr><tr><td>-3</td><td><sup>13</sup>&frasl;<sub>8</sub> e<sup>3</sup> - <sup>25</sup>&frasl;<sub>128</sub> e<sup>5</sup></td><td>-<sup>27</sup>&frasl;<sub>640</sub> e<sup>5</sup></td><td><sup>53</sup>&frasl;<sub>16</sub> e<sup>3</sup> + <sup>393</sup>&frasl;<sub>256</sub> e<sup>5</sup></td><td><sup>81</sup>&frasl;<sub>1280</sub> e<sup>5</sup></td></tr><tr><td>-2</td><td><sup>5</sup>&frasl;<sub>4</sub> e<sup>2</sup> + <sup>1</sup>&frasl;<sub>6</sub> e<sup>4</sup> + <sup>21</sup>&frasl;<sub>64</sub> e<sup>6</sup></td><td>-<sup>1</sup>&frasl;<sub>24</sub> e<sup>4</sup> - <sup>3</sup>&frasl;<sub>80</sub> e<sup>6</sup></td><td><sup>9</sup>&frasl;<sub>4</sub> e<sup>2</sup> + <sup>7</sup>&frasl;<sub>4</sub> e<sup>4</sup> + <sup>141</sup>&frasl;<sub>64</sub> e<sup>6</sup></td><td><sup>1</sup>&frasl;<sub>24</sub> e<sup>4</sup> + <sup>7</sup>&frasl;<sub>240</sub> e<sup>6</sup></td></tr><tr><td>-1</td><td>e + <sup>3</sup>&frasl;<sub>8</sub> e<sup>3</sup> + <sup>65</sup>&frasl;<sub>192</sub> e<sup>5</sup></td><td>-<sup>1</sup>&frasl;<sub>24</sub> e<sup>3</sup> - <sup>13</sup>&frasl;<sub>384</sub> e<sup>5</sup></td><td><sup>3</sup>&frasl;<sub>2</sub> e + <sup>27</sup>&frasl;<sub>16</sub> e<sup>3</sup> + <sup>261</sup>&frasl;<sub>128</sub> e<sup>5</sup></td><td><sup>1</sup>&frasl;<sub>48</sub> e<sup>3</sup> + <sup>11</sup>&frasl;<sub>768</sub> e<sup>5</sup></td></tr><tr><td>0</td><td>1 + <sup>1</sup>&frasl;<sub>2</sub> e<sup>2</sup> + <sup>3</sup>&frasl;<sub>8</sub> e<sup>4</sup> + <sup>5</sup>&frasl;<sub>16</sub> e<sup>6</sup></td><td></td><td>1 + <sup>3</sup>&frasl;<sub>2</sub> e<sup>2</sup> + <sup>15</sup>&frasl;<sub>8</sub> e<sup>4</sup> + <sup>35</sup>&frasl;<sub>16</sub> e<sup>6</sup></td><td></td></tr><tr><td>1</td><td>e + <sup>3</sup>&frasl;<sub>8</sub> e<sup>3</sup> + <sup>65</sup>&frasl;<sub>192</sub> e<sup>5</sup></td><td>-e + <sup>3</sup>&frasl;<sub>8</sub> e<sup>3</sup> + <sup>7</sup>&frasl;<sub>192</sub> e<sup>5</sup></td><td><sup>3</sup>&frasl;<sub>2</sub> e + <sup>27</sup>&frasl;<sub>16</sub> e<sup>3</sup> + <sup>261</sup>&frasl;<sub>128</sub> e<sup>5</sup></td><td>-<sup>1</sup>&frasl;<sub>2</sub> e + <sup>1</sup>&frasl;<sub>16</sub> e<sup>3</sup> - <sup>5</sup>&frasl;<sub>384</sub> e<sup>5</sup></td></tr><tr><td>2</td><td><sup>5</sup>&frasl;<sub>4</sub> e<sup>2</sup> + <sup>1</sup>&frasl;<sub>6</sub> e<sup>4</sup> + <sup>21</sup>&frasl;<sub>64</sub> e<sup>6</sup></td><td>1 - <sup>7</sup>&frasl;<sub>2</sub> e<sup>2</sup> + <sup>29</sup>&frasl;<sub>16</sub> e<sup>4</sup> - <sup>53</sup>&frasl;<sub>288</sub> e<sup>6</sup></td><td><sup>9</sup>&frasl;<sub>4</sub> e<sup>2</sup> + <sup>7</sup>&frasl;<sub>4</sub> e<sup>4</sup> + <sup>141</sup>&frasl;<sub>64</sub> e<sup>6</sup></td><td>1 - <sup>5</sup>&frasl;<sub>2</sub> e<sup>2</sup> + <sup>13</sup>&frasl;<sub>16</sub> e<sup>4</sup> - <sup>35</sup>&frasl;<sub>288</sub> e<sup>6</sup></td></tr><tr><td>3</td><td><sup>13</sup>&frasl;<sub>8</sub> e<sup>3</sup> - <sup>25</sup>&frasl;<sub>128</sub> e<sup>5</sup></td><td>3e - <sup>69</sup>&frasl;<sub>8</sub> e<sup>3</sup> + <sup>369</sup>&frasl;<sub>64</sub> e<sup>5</sup></td><td><sup>53</sup>&frasl;<sub>16</sub> e<sup>3</sup> + <sup>393</sup>&frasl;<sub>256</sub> e<sup>5</sup></td><td><sup>7</sup>&frasl;<sub>2</sub> e - <sup>123</sup>&frasl;<sub>16</sub> e<sup>3</sup> + <sup>489</sup>&frasl;<sub>128</sub> e<sup>5</sup></td></tr><tr><td>4</td><td><sup>103</sup>&frasl;<sub>48</sub> e<sup>4</sup> - <sup>129</sup>&frasl;<sub>160</sub> e<sup>6</sup></td><td><sup>13</sup>&frasl;<sub>2</sub> e<sup>2</sup> - <sup>55</sup>&frasl;<sub>3</sub> e<sup>4</sup> + <sup>239</sup>&frasl;<sub>16</sub> e<sup>6</sup></td><td><sup>77</sup>&frasl;<sub>16</sub> e<sup>4</sup> + <sup>129</sup>&frasl;<sub>160</sub> e<sup>6</sup></td><td><sup>17</sup>&frasl;<sub>2</sub> e<sup>2</sup> - <sup>115</sup>&frasl;<sub>6</sub> e<sup>4</sup> + <sup>601</sup>&frasl;<sub>48</sub> e<sup>6</sup></td></tr><tr><td>5</td><td><sup>1097</sup>&frasl;<sub>384</sub> e<sup>5</sup></td><td><sup>295</sup>&frasl;<sub>24</sub> e<sup>3</sup> - <sup>13745</sup>&frasl;<sub>384</sub> e<sup>5</sup></td><td><sup>1773</sup>&frasl;<sub>256</sub> e<sup>5</sup></td><td><sup>845</sup>&frasl;<sub>48</sub> e<sup>3</sup> - <sup>32525</sup>&frasl;<sub>768</sub> e<sup>5</sup></td></tr><tr><td>6</td><td><sup>1223</sup>&frasl;<sub>320</sub> e<sup>6</sup></td><td><sup>345</sup>&frasl;<sub>16</sub> e<sup>4</sup> - <sup>10569</sup>&frasl;<sub>160</sub> e<sup>6</sup></td><td><sup>3167</sup>&frasl;<sub>320</sub> e<sup>6</sup></td><td><sup>533</sup>&frasl;<sub>16</sub> e<sup>4</sup> - <sup>13827</sup>&frasl;<sub>160</sub> e<sup>6</sup></td></tr><tr><td>7</td><td></td><td><sup>69251</sup>&frasl;<sub>1920</sub> e<sup>5</sup></td><td></td><td><sup>228347</sup>&frasl;<sub>3840</sub> e<sup>5</sup></td></tr><tr><td>8</td><td></td><td><sup>42037</sup>&frasl;<sub>720</sub> e<sup>6</sup></td><td></td><td><sup>73369</sup>&frasl;<sub>720</sub> e<sup>6</sup></td></tr></tbody></table>

**Notes.** $X_k^{l,-2}(e) = X_{-k}^{l,2}(e)$. The exact expression of these coefficients is given by expression (72).

$$I_{12} = \frac{\pi R^4}{20g} \sum_{k=-\infty}^{+\infty} \left\{ \hat{p}_2(kn - 2\omega) \left[ 2(1 - x^2) X_k^{-2,0} + (1 + x)^2 e^{i2\varpi} X_k^{-2,2} + (1 - x)^2 e^{-i2\varpi} X_k^{-2,-2} \right] - \hat{p}_2(kn + 2\omega) \left[ (1 + x)^2 e^{-i2\varpi} X_k^{-2,-2} + (1 - x)^2 e^{i2\varpi} X_k^{-2,2} + 2(1 - x^2) X_k^{-2,0} \right] \right\} e^{ikM}, \quad (77)$$

$$I_{33} = \frac{\pi R^4}{5g} \sum_{k=-\infty}^{+\infty} \left\{ \hat{p}_2(kn) \left[ 2 \left(x^2 - \frac{1}{3}\right) X_k^{-2,0} + (1 - x^2) \left( e^{i2\varpi} X_k^{-2,2} + e^{-i2\varpi} X_k^{-2,-2} \right) \right] \right\} e^{ikM}, \quad (78)$$

$$I_{22} = \Delta I - \frac{I_{33}}{2} \quad \text{and} \quad I_{11} = -\Delta I - \frac{I_{33}}{2}, \quad (79)$$

with

$$\Delta I = \frac{\pi R^4}{20g} \sum_{k=-\infty}^{+\infty} \left\{ \hat{p}_2(kn - 2\omega) \left[ 2(1 - x^2) X_k^{-2,0} + (1 + x)^2 e^{i2\varpi} X_k^{-2,2} + (1 - x)^2 e^{-i2\varpi} X_k^{-2,-2} \right] + \hat{p}_2(kn + 2\omega) \left[ (1 + x)^2 e^{-i2\varpi} X_k^{-2,2} + (1 - x)^2 e^{i2\varpi} X_k^{-2,0} \right] \right\} e^{ikM}. \quad (80)$$

The orbital and spin evolution of the planet under the action of thermal atmospheric tides is completely described by the set of Eqs. (65), (67) and (75)–(80).

## 4. 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 (Sect. 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.

### 4.1. 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))

$$-\frac{3Gm_{\star}}{r^3} I_{23} \hat{x}_1 \hat{x}_3 = \frac{3Gm_{\star}}{2r^3} I_{23} \sin \varepsilon \sin(2\varpi + 2\nu) \\ = \frac{3Gm_{\star}}{4ia^3} I_{23} \sin \varepsilon \sum_{k'=-\infty}^{+\infty} \left( e^{i2\varpi} X_{k'}^{-3,2} - e^{-i2\varpi} X_{k'}^{-3,-2} \right) e^{ik'M}. \quad (81)$$

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

$$ \begin{aligned} & \left\langle -\frac{3Gm_{\star}}{r^3} I_{23} \hat{x}_1 \hat{x}_3 \right\rangle_{M,\omega} \\ &= \frac{3\pi Gm_{\star}R^4}{40ia^3g} \sin^2\epsilon \sum_{k=-\infty}^{+\infty} \left\{ \hat{p}_2(kn - \omega) \left[ (1+x)X_k^{-2,2}X_{-k}^{-3,-2} \right. \right. \\ & \qquad \left. \left. + (1-x)X_k^{-2,-2}X_{-k}^{-3,2} \right] \right. \\ & \qquad \left. - \hat{p}_2(kn + \omega) \left[ (1-x)X_k^{-2,2}X_{-k}^{-3,-2} \right. \right. \\ & \qquad \left. \left. + (1+x)X_k^{-2,-2}X_{-k}^{-3,2} \right] \right\}. \end{aligned} \quad (82) $$

Finally, we decompose the surface pressure variations $\hat{p}_2(kn \pm \omega)$ in its real and imaginary parts (Eq. (53)), make use of their parity properties, and use the simplification $X_{-k}^{\ell,m} = X_k^{\ell,-m}$ to write

$$ \left\langle -\frac{3Gm_{\star}}{r^3} I_{23} \hat{x}_1 \hat{x}_3 \right\rangle_{M,\omega} = \frac{3\pi Gm_{\star}R^4}{20a^3g} \sin^2\epsilon \times \sum_{k=-\infty}^{+\infty} b(\omega - kn) \left[ (1+x)X_k^{-2,2}X_k^{-3,2} + (1-x)X_k^{-2,-2}X_k^{-3,-2} \right]. \quad (83) $$

This last arrangement is very useful, because we are able to combine terms in $b(\omega \pm kn)$ in a single term $b(\omega - 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(\omega - 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 $\omega$, 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,

$$ E = -\frac{Gm_{\star}m}{2a}, \quad (84) $$

because it can be directly obtained from the power of the tidal force (Eq. (60)) as

$$ \dot{E} = \dot{r}_{\star} \cdot \mathbf{F}, \quad (85) $$

and thus provides a simpler expression, from which we can later easily derive the evolution of the eccentricity.

## 4.2. 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

$$ \langle \mathbf{T} \rangle_{M,\omega} = T_p \mathbf{p} + T_q \mathbf{q} + T_s \mathbf{s}. \quad (86) $$

For zero obliquity ($\epsilon = 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 \propto \sin \epsilon$; 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

$$ \langle \mathbf{T} \rangle_{M,\omega} = \bar{T}_1 \mathbf{k} + \bar{T}_2 \mathbf{s} + \bar{T}_3 \mathbf{k} \times \mathbf{s}, \quad (87) $$

where

$$ \bar{T}_1 = \frac{T_q}{\sin \epsilon}, \quad \bar{T}_2 = T_s - \cos \epsilon \frac{T_q}{\sin \epsilon}, \quad \bar{T}_3 = \frac{T_p}{\sin \epsilon}, \quad (88) $$

with

$$ \begin{aligned} \bar{T}_1 = & -\frac{3\mathcal{K}}{32} \sum_{k=-\infty}^{+\infty} \left\{ 3b(kn) \left[ (1-x^2)(X_k^{-2,2}X_k^{-3,2} - X_k^{-2,-2}X_k^{-3,-2}) \right] \right. \\ & \left. + 2b(\omega-kn) \left[ (1-x)^2(2+x)X_k^{-2,-2}X_k^{-3,-2} \right. \right. \\ & \left. \left. + 4x^3X_k^{-2,0}X_k^{-3,0} - (1+x)^2(2-x)X_k^{-2,2}X_k^{-3,2} \right] \right. \\ & \left. + b(2\omega-kn) \left[ 4x(1-x^2)X_k^{-2,0}X_k^{-3,0} \right. \right. \\ & \left. \left. + (1-x)^3X_k^{-2,-2}X_k^{-3,-2} - (1+x)^3X_k^{-2,2}X_k^{-3,2} \right] \right\}, \end{aligned} \quad (89) $$

$$ \begin{aligned} \bar{T}_2 = & -\frac{3\mathcal{K}}{32} \sum_{k=-\infty}^{+\infty} \left\{ 3b(kn) \left[ x(1-x^2)(X_k^{-2,-2}X_k^{-3,-2} - X_k^{-2,2}X_k^{-3,2}) \right] \right. \\ & \left. - 2b(\omega-kn) \left[ (1-x)^2(1+2x)X_k^{-2,-2}X_k^{-3,-2} \right. \right. \\ & \left. \left. + 4x^2X_k^{-2,0}X_k^{-3,0} + (1+x)^2(1-2x)X_k^{-2,2}X_k^{-3,2} \right] \right. \\ & \left. - b(2\omega-kn) \left[ (1-x)^3X_k^{-2,-2}X_k^{-3,-2} \right. \right. \\ & \left. \left. + 4(1-x^2)X_k^{-2,0}X_k^{-3,0} + (1+x)^3X_k^{-2,2}X_k^{-3,2} \right] \right\}, \end{aligned} \quad (90) $$

$$ \begin{aligned} \bar{T}_3 = & \frac{3\mathcal{K}}{32} \sum_{k=-\infty}^{+\infty} \left\{ xa(kn) \left[ -4(1-3x^2)X_k^{-2,0}X_k^{-3,0} \right. \right. \\ & \left. \left. - 3(1-x^2)(X_k^{-2,-2}X_k^{-3,-2} + X_k^{-2,2}X_k^{-3,2}) \right] \right. \\ & \left. + 2a(\omega-kn) \left[ -(1-x)^2(1+2x)X_k^{-2,-2}X_k^{-3,-2} \right. \right. \\ & \left. \left. + 4x(1-2x^2)X_k^{-2,0}X_k^{-3,0} + (1+x)^2(1-2x)X_k^{-2,2}X_k^{-3,2} \right] \right. \\ & \left. - a(2\omega-kn) \left[ 4x(1-x^2)X_k^{-2,0}X_k^{-3,0} \right. \right. \\ & \left. \left. + (1-x)^3X_k^{-2,-2}X_k^{-3,-2} - (1+x)^3X_k^{-2,2}X_k^{-3,2} \right] \right\}, \end{aligned} \quad (91) $$

and

$$ \mathcal{K} = \frac{4\pi Gm_{\star}R^4}{5a^3g} = \frac{3m_{\star}}{5\bar{\rho}} \left(\frac{R}{a}\right)^3, \quad (92) $$

where $\bar{\rho} = 3g/(4\pi G\mathcal{R})$ is the mean density of the planet. These expressions were obtained using the algebraic manipulators Maxima<fnref n="2" /> and TRIP (Gastineau & Laskar 2011).

<footnote n="2">² Maxima, a Computer Algebra System (Version 5.46.0), https://maxima.sourceforge.io/</footnote>
Finally, for the secular power (Eq. (85)) we get

$$ \langle \dot{E} \rangle_{M,\omega} = -\frac{\mathcal{K}}{64} \sum_{k=-\infty}^{+\infty} kn \left\{ b(kn) \left[ 4(1-3x^2)^2 X_k^{-2,0} X_k^{-3,0} \right. \right. \\ 
+ 9(1-x^2)^2 (X_k^{-2,-2} X_k^{-3,-2} + X_k^{-2,2} X_k^{-3,2}) \left. \right] \\ 
- 12 b(\omega-kn)(1-x^2) \left[ 4x^2 X_k^{-2,0} X_k^{-3,0} \right. \\ 
+ (1-x)^2 X_k^{-2,-2} X_k^{-3,-2} + (1+x)^2 X_k^{-2,2} X_k^{-3,2} \left. \right] \\ 
- 3 b(2\omega-kn) \left[ 4(1-x^2)^2 X_k^{-2,0} X_k^{-3,0} \right. \\ 
+ (1-x)^4 X_k^{-2,-2} X_k^{-3,-2} + (1+x)^4 X_k^{-2,2} X_k^{-3,2} \left. \right] \right\}. \quad (93) $$

### 4.3. Orbital and spin evolution

The set of Eqs. (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))

$$ a = -\frac{\beta\mu}{2E}, \quad (94) $$

the eccentricity from the orbital angular momentum (Eq. (55))

$$ e = \sqrt{1 - \frac{(\mathbf{G} \cdot \mathbf{k})^2}{\beta^2 \mu a}} = \sqrt{1 - \frac{\mathbf{G} \cdot \mathbf{G}}{\beta^2 \mu a}}, \quad (95) $$

and the rotation rate from the rotational angular momentum (Eq. (58)),

$$ \omega = \frac{\mathbf{L} \cdot \mathbf{s}}{C} = \frac{\sqrt{\mathbf{L} \cdot \mathbf{L}}}{C}. \quad (96) $$

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)

$$ \cos \varepsilon = \mathbf{k} \cdot \mathbf{s} = \frac{\mathbf{G} \cdot \mathbf{L}}{\sqrt{(\mathbf{G} \cdot \mathbf{G})(\mathbf{L} \cdot \mathbf{L})}}. \quad (97) $$

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),

$$ \dot{a} = \frac{2a^2}{\beta\mu}\dot{E}, \quad (98) $$

while the eccentricity evolution can be computed from expressions (87), (95) and (98) as

$$ \dot{e} = \frac{1 - e^2}{2ae} \dot{a} - \frac{\mathbf{G} \cdot \mathbf{T}}{\beta^2 \mu a e} = \frac{\sqrt{1 - e^2}}{\beta na^2 e} \left( \frac{\sqrt{1 - e^2}}{n} \dot{E} - \bar{T}_1 - \bar{T}_2 x \right), \quad (99) $$

that is, using expressions (89), (90), and (93),

$$ \dot{e} = -\frac{\mathcal{K}}{64} \frac{\sqrt{1-e^2}}{\beta na^2 e} \sum_{k=-\infty}^{+\infty} \left\{ b(kn) \left[ 4(1-3x^2)^2 X_k^{-2,0} X_k^{-3,0} k \sqrt{1-e^2} \right. \right. \\ 
+ 9(1-x^2)^2 X_k^{-2,-2} X_k^{-3,-2} (2+k\sqrt{1-e^2}) \\ 
- 9(1-x^2)^2 X_k^{-2,2} X_k^{-3,2} (2-k\sqrt{1-e^2}) \left. \right] \\ 
- 12 b(\omega-kn)(1-x^2) \left[ 4x^2 X_k^{-2,0} X_k^{-3,0} k \sqrt{1-e^2} \right. \\ 
+ (1-x)^2 X_k^{-2,-2} X_k^{-3,-2} (2+k\sqrt{1-e^2}) \\ 
- (1+x)^2 X_k^{-2,2} X_k^{-3,2} (2-k\sqrt{1-e^2}) \left. \right] \\ 
- 3 b(2\omega-kn) \left[ 4(1-x^2)^2 X_k^{-2,0} X_k^{-3,0} k \sqrt{1-e^2} \right. \\ 
+ (1-x)^4 X_k^{-2,-2} X_k^{-3,-2} (2+k\sqrt{1-e^2}) \\ 
- (1+x)^4 X_k^{-2,2} X_k^{-3,2} (2-k\sqrt{1-e^2}) \left. \right] \right\}. \quad (100) $$

For the rotation rate evolution, we have the following from Eqs. (67) and (96):

$$ \dot{\omega} = -\frac{\mathbf{T} \cdot \mathbf{s}}{C} = -\frac{\bar{T}_1 x + \bar{T}_2}{C}, \quad (101) $$

that is, using expressions (89) and (90)

$$ \dot{\omega} = -\frac{3\mathcal{K}}{32C} \sum_{k=-\infty}^{+\infty} \left\{ 2b(\omega - kn)(1 - x^2) \left[ 4x^2 X_k^{-2,0} X_k^{-3,0} \right. \right. \\ 
+ (1 + x)^2 X_k^{-2,2} X_k^{-3,2} + (1 - x)^2 X_k^{-2,-2} X_k^{-3,-2} \left. \right] \\ 
+ b(2\omega - kn) \left[ 4(1 - x^2)^2 X_k^{-2,0} X_k^{-3,0} \right. \\ 
+ (1 + x)^4 X_k^{-2,2} X_k^{-3,2} + (1 - x)^4 X_k^{-2,-2} X_k^{-3,-2} \left. \right] \right\}. \quad (102) $$

The obliquity (or inclination) evolution is given from expressions (66), (67) and (97),

$$ \dot{e} = \frac{\mathbf{T} \cdot \mathbf{k} - \mathbf{T} \cdot \mathbf{s} \cos \varepsilon}{\sin \varepsilon \sqrt{\mathbf{L} \cdot \mathbf{L}}} - \frac{\mathbf{T} \cdot \mathbf{s} - \mathbf{T} \cdot \mathbf{k} \cos \varepsilon}{\sin \varepsilon \sqrt{\mathbf{G} \cdot \mathbf{G}}} \\ 
= \left( \frac{\bar{T}_1}{C\omega} - \frac{\bar{T}_2}{\beta \sqrt{\mu a(1 - e^2)}} \right) \sin \varepsilon \approx \frac{\bar{T}_1}{C\omega} \sin \varepsilon. \quad (103) $$

In general, for a planet around a star, we have $|\mathbf{L}| \ll |\mathbf{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, $\dot{\Omega}$, and the precession speed of the spin axis, $\dot{\psi}$. The line of nodes is aligned with the vector $\mathbf{p}$ (Fig. 1) and thus

$$ \dot{\Omega} = \frac{\dot{G}}{|\mathbf{G}|} \cdot \mathbf{p} = \frac{\mathbf{T} \cdot (\mathbf{k} \times \mathbf{s})}{|\mathbf{G}| \sin \varepsilon} = \frac{\bar{T}_3 \sin \varepsilon}{\beta \sqrt{\mu a(1 - e^2)}}, \quad (104) $$

$$ \dot{\psi} = \frac{\dot{L}}{|\mathbf{L}|} \cdot \mathbf{p} = - \frac{\mathbf{T} \cdot (\mathbf{k} \times \mathbf{s})}{|\mathbf{L}| \sin \varepsilon} = - \frac{\bar{T}_3 \sin \varepsilon}{C\omega}. \quad (105) $$
### 4.4. Energy balance

The average total energy transferred to the planet due to thermal atmospheric tides is given by

$$ \dot{U} = \langle \dot{E} + \dot{K} \rangle_{M,\omega}, \qquad (106) $$

where E is the orbital energy (Eq. (84)),

$$ K = \frac{\omega \cdot L}{2} \qquad (107) $$

is the rotational energy, and

$$ \dot{K} = C\omega\dot{\omega}, \qquad (108) $$

where $\dot{\omega}$ is given by expression (101).

The total energy is then obtained by combining expressions (93) and (102) as

$$ \begin{aligned} \dot{U} = -\frac{\mathcal{K}}{64} \sum_{k=-\infty}^{+\infty} & \left\{ kn b(kn) \left[ 4(1-3x^2)^2 X_k^{-2,0} X_k^{-3,0} \right. \right. \\ & \left. \left. \qquad + 9(1-x^2)^2 (X_k^{-2,-2} X_k^{-3,-2} + X_k^{-2,2} X_k^{-3,2}) \right] \right. \\ & \left. \qquad + 12(\omega-kn)b(\omega-kn)(1-x^2) \left[ 4x^2 X_k^{-2,0} X_k^{-3,0} \right. \right. \\ & \left. \left. \qquad + (1-x)^2 X_k^{-2,-2} X_k^{-3,-2} + (1+x)^2 X_k^{-2,2} X_k^{-3,2} \right] \right. \\ & \left. \qquad + 3(2\omega-kn)b(2\omega-kn) \left[ 4(1-x^2)^2 X_k^{-2,0} X_k^{-3,0} \right. \right. \\ & \left. \left. \qquad + (1-x)^4 X_k^{-2,-2} X_k^{-3,-2} + (1+x)^4 X_k^{-2,2} X_k^{-3,2} \right] \right\}. \end{aligned} \qquad (109) $$

We note that regardless of the orbital and spin parameters, the total energy transferred is always positive since $b(\sigma)$ is an odd function and $b(|\sigma|) < 0$ (Eq. (54)).

## 5. 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.

### 5.1. Tidal torque and power

For the average tidal torque (Eq. (87)), we have

$$ \begin{aligned} \bar{T}_1 = -\mathcal{K} & \left\{ b(3n) \frac{189}{32} (1-x^2) e^2 \right. \\ & + b(2n) \frac{9}{16} (1-x^2) (1-6e^2) \\ & + b(n) \frac{9}{32} (1-x^2) e^2 \\ & + b(\omega+3n) \frac{63}{32} (1-x)^2 (2+x) e^2 \\ & + b(\omega+2n) \frac{3}{16} (1-x)^2 (2+x) (1-6e^2) \\ & \left. + b(\omega+n) \frac{3}{32} (2-3x+13x^3) e^2 \right\} \end{aligned} $$

$$ \begin{aligned} & + b(\omega) \frac{3}{4} x^3 (1 + 2e^2) \\ & - b(\omega - n) \frac{3}{32} (2 + 3x - 13x^3) e^2 \\ & - b(\omega - 2n) \frac{3}{16} (1 + x)^2 (2 - x) (1 - 6e^2) \\ & - b(\omega - 3n) \frac{63}{32} (1 + x)^2 (2 - x) e^2 \\ & + b(2\omega + 3n) \frac{63}{64} (1 - x)^3 e^2 \\ & + b(2\omega + 2n) \frac{3}{32} (1 - x)^3 (1 - 6e^2) \\ & + b(2\omega + n) \frac{3}{64} (1 - x) (1 + 10x + 13x^2) e^2 \\ & + b(2\omega) \frac{3}{8} x (1 - x^2) (1 + 2e^2) \\ & - b(2\omega - n) \frac{3}{64} (1 + x) (1 - 10x + 13x^2) e^2 \\ & - b(2\omega - 2n) \frac{3}{32} (1 + x)^3 (1 - 6e^2) \\ & - b(2\omega - 3n) \frac{63}{64} (1 + x)^3 e^2 \end{aligned} \qquad (110) $$

$$ \begin{aligned} \bar{T}_2 = \mathcal{K} & \left\{ b(3n) \frac{189}{32} x (1 - x^2) e^2 \right. \\ & + b(2n) \frac{9}{16} x (1 - x^2) (1 - 6e^2) \\ & + b(n) \frac{9}{32} x (1 - x^2) e^2 \\ & + b(\omega + 3n) \frac{63}{32} (1 - x)^2 (1 + 2x) e^2 \\ & + b(\omega + 2n) \frac{3}{16} (1 - x)^2 (1 + 2x) (1 - 6e^2) \\ & + b(\omega + n) \frac{3}{32} (1 + 9x^2 + 2x^3) e^2 \\ & + b(\omega) \frac{3}{4} x^2 (1 + 2e^2) \\ & + b(\omega - n) \frac{3}{32} (1 + 9x^2 - 2x^3) e^2 \\ & + b(\omega - 2n) \frac{3}{16} (1 + x)^2 (1 - 2x) (1 - 6e^2) \\ & + b(\omega - 3n) \frac{63}{32} (1 + x)^2 (1 - 2x) e^2 \\ & + b(2\omega + 3n) \frac{63}{64} (1 - x)^3 e^2 \\ & + b(2\omega + 2n) \frac{3}{32} (1 - x)^3 (1 - 6e^2) \\ & + b(2\omega + n) \frac{3}{64} (1 - x) (13 + 10x + x^2) e^2 \\ & + b(2\omega) \frac{3}{8} (1 - x^2) (1 + 2e^2) \\ & + b(2\omega - n) \frac{3}{64} (1 + x) (13 - 10x + x^2) e^2 \\ & + b(2\omega - 2n) \frac{3}{32} (1 + x)^3 (1 - 6e^2) \\ & + b(2\omega - 3n) \frac{63}{64} (1 + x)^3 e^2 \end{aligned} \qquad (111) $$
$$
\begin{equation}
\begin{aligned}
\bar{T}_3 = -\mathcal{K} & \left\{ a(3n) \frac{189}{32} x(1-x^2)e^2 \right. \\
& + a(2n) \frac{9}{16} x(1-x^2)(1-6e^2) \\
& + a(n) \frac{9}{32} x(5-13x^2)e^2 \\
& + a(0) \frac{3}{8} x(1-3x^2)(1+2e^2) \\
& + a(\omega+3n) \frac{63}{32} (1-x)^2 (1+2x)e^2 \\
& + a(\omega+2n) \frac{3}{16} (1-x)^2 (1+2x)(1-6e^2) \\
& + a(\omega+n) \frac{3}{32} (1-12x-3x^2+26x^3)e^2 \\
& - a(\omega) \frac{3}{4} x(1-2x^2)(1+2e^2) \\
& - a(\omega-n) \frac{3}{32} (1+12x-3x^2-26x^3)e^2 \\
& - a(\omega-2n) \frac{3}{16} (1+x)^2 (1-2x)(1-6e^2) \\
& - a(\omega-3n) \frac{63}{32} (1+x)^2 (1-2x)e^2 \\
& + a(2\omega+3n) \frac{63}{64} (1-x)^3 e^2 \\
& + a(2\omega+2n) \frac{3}{32} (1-x)^3 (1-6e^2) \\
& + a(2\omega+n) \frac{3}{64} (1-x)(1+10x+13x^2)e^2 \\
& + a(2\omega) \frac{3}{8} x(1-x^2)(1+2e^2) \\
& - a(2\omega-n) \frac{3}{64} (1+x)(1-10x+13x^2)e^2 \\
& - a(2\omega-2n) \frac{3}{32} (1+x)^3 (1-6e^2) \\
& - a(2\omega-3n) \frac{63}{64} (1+x)^3 e^2 \}
\end{aligned}
\tag{112}
\end{equation}
$$

$$
- a(\omega - 2n) \frac{3}{16} (1 + x)^2 (1 - 2x) (1 - 6e^2)
$$

$$
- a(\omega - 3n) \frac{63}{32} (1 + x)^2 (1 - 2x) e^2 \\
+ a(2\omega + 3n) \frac{63}{64} (1 - x)^3 e^2
$$

$$
+ a(2\omega + 2n) \frac{3}{32} (1 - x)^3 (1 - 6e^2)
$$

$$
+ a(2\omega + n) \frac{3}{64} (1 - x)(1 + 10x + 13x^2)e^2
$$

$$
+ a(2\omega) \frac{3}{8} x(1 - x^2)(1 + 2e^2)
$$

$$
- a(2\omega - n) \frac{3}{64} (1 + x)(1 - 10x + 13x^2)e^2
$$

$$
- a(2\omega - 2n) \frac{3}{32} (1 + x)^3 (1 - 6e^2)
$$

$$
- a(2\omega - 3n) \frac{63}{64} (1 + x)^3 e^2 \}
$$

and for the average power (Eq. (93))

$$
\begin{equation}
\begin{aligned}
\langle \dot{E} \rangle_{M,\bar{\omega}} = -n\mathcal{K} & \left\{ b(3n) \frac{567}{64} (1-x^2)^2 e^2 \right. \\
& + b(2n) \frac{9}{16} (1-x^2)^2 (1-6e^2) \\
& + b(n) \frac{3}{64} (7-30x^2+39x^4) e^2 \\
& + b(\omega+3n) \frac{189}{32} (1-x^2)(1-x)^2 e^2 \\
& + b(\omega+2n) \frac{3}{8} (1-x^2)(1-x)^2 (1-6e^2) \\
& + b(\omega+n) \frac{3}{32} (1-x^2)(1-2x+13x^2) e^2 \\
& - b(\omega-n) \frac{3}{32} (1-x^2)(1+2x+13x^2) e^2
\end{aligned}
\tag{113}
\end{equation}
$$

5.2. 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))

$$
\begin{equation}
\begin{aligned}
\dot{e} = -\frac{Ke}{\beta na^2} & \left\{ b(3n) \frac{189}{64} (1-x^2)^2 \right. \\
& - b(2n) \frac{9}{32} (1-x^2)^2 \\
& + b(n) \frac{3}{64} (1-18x^2+33x^4) \\
& + b(\omega+3n) \frac{63}{32} (1-x^2)(1-x)^2 \\
& - b(\omega+2n) \frac{3}{16} (1-x^2)(1-x)^2 \\
& - b(\omega+n) \frac{3}{32} (1-x^2)(1-2x-11x^2) \\
& + b(\omega-n) \frac{3}{32} (1-x^2)(1+2x-11x^2) \\
& + b(\omega-2n) \frac{3}{16} (1-x^2)(1+x)^2 \\
& - b(\omega-3n) \frac{63}{32} (1-x^2)(1+x)^2 \\
& + b(2\omega+3n) \frac{63}{128} (1-x)^4 \\
& - b(2\omega+2n) \frac{3}{64} (1-x)^4 \\
& + b(2\omega+n) \frac{3}{128} (1-x)^2 (11+26x+11x^2) \\
& - b(2\omega-n) \frac{3}{128} (1+x)^2 (11-26x+11x^2) \\
& + b(2\omega-2n) \frac{3}{64} (1+x)^4 \\
& - b(2\omega-3n) \frac{63}{128} (1+x)^4 \}
\end{aligned}
\tag{114}
\end{equation}
$$
The evolution of the obliquity is directly obtained by replacing Eq. (110) in Eq. (103), and for the rotation rate (Eq. (102)),

$$ \dot{\omega} = -\frac{\mathcal{K}}{C} \left\{ \begin{aligned} & b(\omega + 3n) \frac{63}{32} (1-x^2)(1-x)^2 e^2 \\ & + b(\omega + 2n) \frac{3}{16} (1-x^2)(1-x)^2(1-6e^2) \\ & + b(\omega + n) \frac{3}{32} (1-x^2)(1-2x+13x^2)e^2 \\ & + b(\omega) \frac{3}{4} x^2 (1-x^2)(1+2e^2) \\ & + b(\omega - n) \frac{3}{32} (1-x^2)(1+2x+13x^2)e^2 \\ & + b(\omega - 2n) \frac{3}{16} (1-x^2)(1+x)^2(1-6e^2) \\ & + b(\omega - 3n) \frac{63}{32} (1-x^2)(1+x)^2 e^2 \\ & + b(2\omega + 3n) \frac{63}{64} (1-x)^4 e^2 \\ & + b(2\omega + 2n) \frac{3}{32} (1-x)^4(1-6e^2) \\ & + b(2\omega + n) \frac{3}{64} (1-x)^2(13+22x+13x^2)e^2 \\ & + b(2\omega) \frac{3}{8} (1-x^2)^2(1+2e^2) \\ & + b(2\omega - n) \frac{3}{64} (1+x)^2(13-22x+13x^2)e^2 \\ & + b(2\omega - 2n) \frac{3}{32} (1+x)^4(1-6e^2) \\ & + b(2\omega - 3n) \frac{63}{64} (1+x)^4 e^2 \end{aligned} \right\}. \quad (115) $$

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_{\star}^{-5}$, while here the potential only depends on $r_{\star}^{-3}$ (Eq. (35)), and the extra factor in $r_{\star}^{-2}$ is included in the forcing function (Eq. (36)).

### 5.3. Energy balance

The total energy is obtained from expression (109) as

$$ \begin{aligned} \dot{U} = -\mathcal{K} \left\{ \begin{aligned} & (3n) b(3n) \frac{189}{64} (1-x^2)^2 e^2 \\ & + (2n) b(2n) \frac{9}{32} (1-x^2)^2 (1-6e^2) \\ & + (n) b(n) \frac{3}{64} (7-30x^2+39x^4)e^2 \\ & + (\omega + 3n) b(\omega + 3n) \frac{63}{32} (1-x^2)(1-x)^2 e^2 \\ & + (\omega + 2n) b(\omega + 2n) \frac{3}{16} (1-x^2)(1-x)^2(1-6e^2) \\ & + (\omega + n) b(\omega + n) \frac{3}{32} (1-x^2)(1-2x+13x^2)e^2 \\ & + (\omega) b(\omega) \frac{3}{4} x^2 (1-x^2)(1+2e^2) \\ & + (\omega - n) b(\omega - n) \frac{3}{32} (1-x^2)(1+2x+13x^2)e^2 \end{aligned} \right\} \end{aligned} $$

$$ \begin{aligned} & + (\omega - 2n) b(\omega - 2n) \frac{3}{16} (1 - x^2)(1 + x)^2 (1 - 6e^2) \\ & + (\omega - 3n) b(\omega - 3n) \frac{63}{32} (1 - x^2)(1 + x)^2 e^2 \\ & + (2\omega + 3n) b(2\omega + 3n) \frac{63}{128} (1 - x)^4 e^2 \\ & + (2\omega + 2n) b(2\omega + 2n) \frac{3}{64} (1 - x)^4 (1 - 6e^2) \\ & + (2\omega + n) b(2\omega + n) \frac{3}{128} (1 - x)^2 (13 + 22x + 13x^2)e^2 \\ & + (2\omega) b(2\omega) \frac{3}{16} (1 - x^2)^2 (1 + 2e^2) \\ & + (2\omega - n) b(2\omega - n) \frac{3}{128} (1 + x)^2 (13 - 22x + 13x^2)e^2 \\ & + (2\omega - 2n) b(2\omega - 2n) \frac{3}{64} (1 + x)^4 (1 - 6e^2) \\ & + (2\omega - 3n) b(2\omega - 3n) \frac{63}{128} (1 + x)^4 e^2 \end{aligned} \quad (116) $$

## 6. 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 ($\epsilon = 0$). In this case, we have $x = 1$ and $s = k$ (Fig. 1).

### 6.1. Tidal torque and power

Using the simplification $s = k$ yields (Eqs. (86) and (87))

$$ \langle \mathbf{T} \rangle_{M,\omega} = T_s \mathbf{k} = (\overline{T}_1 + \overline{T}_2) \mathbf{k}, \quad (117) $$

for the average tidal torque, with

$$ T_s = \frac{3\mathcal{K}}{2} \sum_{k=-\infty}^{+\infty} b(2\omega - kn) X_k^{-2,2} X_k^{-3,2}. \quad (118) $$

For the average power, we get the following from Eq. (93) with $x = 1$:

$$ \langle \dot{E} \rangle_{M,\omega} = -\frac{\mathcal{K}}{4} \sum_{k=-\infty}^{+\infty} kn \left[ b(kn) X_k^{-2,0} X_k^{-3,0} - 3 b(2\omega - kn) X_k^{-2,2} X_k^{-3,2} \right]. \quad (119) $$

### 6.2. Orbital and spin evolution

For the semimajor axis, we have (Eqs. (98) and (119))

$$ \frac{\dot{a}}{a} = -\frac{\mathcal{K}}{\beta na^2} \sum_{k=-\infty}^{+\infty} \frac{k}{2} \left[ b(kn) X_k^{-2,0} X_k^{-3,0} - 3 b(2\omega - kn) X_k^{-2,2} X_k^{-3,2} \right], \quad (120) $$

or, up to $e^2$,

$$ \frac{\dot{a}}{a} = \frac{3\mathcal{K}}{\beta na^2} \left\{ b(2\omega - 2n) (1 - 6e^2) - b(n) \frac{1}{2} e^2 + b(2\omega - n) \frac{1}{4} e^2 + b(2\omega - 3n) \frac{63}{4} e^2 \right\}. \quad (121) $$
Since $x = 1$, for the eccentricity we get (Eq. (100)),

$$ \dot{e} = -\frac{\mathcal{K}}{4} \frac{\sqrt{1-e^2}}{\beta na^2 e} \sum_{k=-\infty}^{+\infty} \left\{ b(kn) X_k^{-2,0} X_k^{-3,0} k \sqrt{1-e^2} + 3 b(2\omega - kn) X_k^{-2,2} X_k^{-3,2} (2-k\sqrt{1-e^2}) \right\}, \quad (122) $$

or, up to $e^2$,

$$ \dot{e} = -\frac{3\mathcal{K}e}{4\beta na^2} \left\{ b(n) + b(2\omega - 2n) + \frac{1}{2}b(2\omega - n) - \frac{21}{2}b(2\omega - 3n) \right\}, \quad (123) $$

and for the rotation rate (Eq. (102))

$$ \dot{\omega} = -\frac{3\mathcal{K}}{2C} \sum_{k=-\infty}^{+\infty} b(2\omega - kn) X_k^{-2,2} X_k^{-3,2}, \quad (124) $$

or, up to $e^2$,

$$ \dot{\omega} = -\frac{3\mathcal{K}}{2C} \left\{ b(2\omega - 2n)(1 - 6e^2) + b(2\omega - n)\frac{1}{2}e^2 + b(2\omega - 3n)\frac{21}{2}e^2 \right\}. \quad (125) $$

The evolution of the obliquity is simply given by $\dot{\epsilon} = 0$ (Eq. (103)) since $\sin \epsilon = 0$, that is, the motion remains planar.

### 6.3. Energy balance

The total energy transferred to the planet due to tides is obtained from expression (109) with $x = 1$,

$$ \dot{U} = -\frac{\mathcal{K}}{4} \sum_{k=-\infty}^{+\infty} \left\{ kn b(kn) X_k^{-2,0} X_k^{-3,0} + 3(2\omega - kn) b(2\omega - kn) X_k^{-2,2} X_k^{-3,2} \right\}. \quad (126) $$

## 7. 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.

*Acknowledgements.* This work was supported by grant SFRH/BD/137958/2018, and by projects CFisUC (UIDB/04564/2020 and UIDP/04564/2020), GRAVITY (PTDC/FIS-AST/7002/2020), PHOBOS (POCI-01-0145-FEDER-029932), and ENGAGE SKA (POCI-01-0145-FEDER-022217), funded by COMPETE 2020 and FCT, Portugal.

## References

Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions (New York: Dover)

Auclair-Desrotour, P., Laskar, J., & Mathis, S. 2017, A&A, 603, A107

Auclair-Desrotour, P., Leconte, J., & Mergny, C. 2019, A&A, 624, A17

Baker, R. D., Schubert, G., & Jones, P. W. 2000, J. Atmos. Sci., 57, 184

Borucki, W. J., Agol, E., Fressin, F., et al. 2013, Science, 340, 587

Boué, G., & Efroimsky, M. 2019, Celest. Mech. Dyn. Astron., 131, 30

Boué, G., & Fabrycky, D. C. 2014, ApJ, 789, 110

Carpenter, R. L. 1964, AJ, 69, 2

Chapman, S., & Lindzen, R. 1970, Atmospheric Tides. Thermal and Gravitational (Dordrecht: Reidel)

Correia, A. C. M. 2006, Earth Planet. Sci. Lett., 252, 398

Correia, A. C. M., & Laskar, J. 2001, Nature, 411, 767

Correia, A. C. M., & Laskar, J. 2003a, J. Geophys. Res. Planets, 108, 5123

Correia, A. C. M., & Laskar, J. 2003b, Icarus, 163, 24

Correia, A. C. M., & Valente, E. F. S. 2022, Celest. Mech. Dyn. Astron., 134, 24

Correia, A. C. M., Laskar, J., & Néron de Surgy, O. 2003, Icarus, 163, 1

Correia, A. C. M., Laskar, J., Farago, F., & Boué, G. 2011, Celest. Mech. Dyn. Astron., 111, 105

Correia, A. C. M., Boué, G., Laskar, J., & Rodríguez, A. 2014, A&A, 571, A50

Correia, A. C. M., Boué, G., & Laskar, J. 2016, Celest. Mech. Dyn. Astron. 126, 189

Cowling, T. G. 1941, MNRAS, 101, 367

Cunha, D., Correia, A. C. M., & Laskar, J. 2015, Int. J. Astrobiol., 14, 233

Dobrovolskis, A. R. 1980, Icarus, 41, 18

Dobrovolskis, A. R., & Ingersoll, A. P. 1980, Icarus, 41, 1

Frouard, J., & Efroimsky, M. 2018, MNRAS, 473, 728

Gastineau, M., & Laskar, J. 2011, ACM Commun. Comput. Algebra, 44, 194

Gold, T., & Soter, S. 1969, Icarus, 11, 356

Goldstein, H. 1950, Classical Mechanics (Reading: Addison-Wesley)

Kaula, W. M. 1964, Rev. Geophys., 2, 661

Komacek, T. D., & Abbot, D. S. 2019, ApJ, 871, 245

Leconte, J., Wu, H., Menou, K., & Murray, N. 2015, Science, 347, 632

Schlecker, M., Pham, D., Burn, R., et al. 2021, A&A, 656, A73

Seiff, A., Kirk, D. B., Young, R. E., et al. 1980, J. Geophys. Res., 85, 7903

Siebert, M. 1961, Adv. Geophys., 7, 105

Ward, W. R. 1974, J. Geophys. Res., 79, 3375

Winn, J. N. 2018, in Handbook of Exoplanets, eds. H. J. Deeg, & J. A. Belmonte (Berlin: Springer), 195

Withers, P., Pratt, R., Bertaux, J. L., & Montmessin, F. 2011, J. Geophys. Res. Planets, 116, E11005

Wordsworth, R., & Kreidberg, L. 2022, ARA&A, 60, 159

Wu, H., Murray, N., Menou, K., Lee, C., & Leconte, J. 2023, Sci. Adv., 9, 27