# Drifts of the substellar points of the TRAPPIST-1 planets

Alexandre Revol, Émeline Bolmont, Mariana Sastre, Gabriel Tobie, Anne-Sophie Libert, Mathilde Kervazo, Sergi Blanco-Cuaresma

► To cite this version:

Alexandre Revol, Émeline Bolmont, Mariana Sastre, Gabriel Tobie, Anne-Sophie Libert, et al. Drifts of the substellar points of the TRAPPIST-1 planets. Astronomy & Astrophysics - A&A, 2024, 691, 10.1051/0004-6361/202451532. insu-04841454

HAL Id: insu-04841454

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

Submitted on 16 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
LETTER TO THE EDITOR

# Drifts of the substellar points of the TRAPPIST-1 planets

Alexandre Revol<fnref n="1" />,<fnref n="2" />,*, Émeline Bolmont<fnref n="1" />,<fnref n="2" />,*, Mariana Sastre<fnref n="3" />, Gabriel Tobie<fnref n="4" />, Anne-Sophie Libert<fnref n="5" />,
Mathilde Kervazo<fnref n="4" />, and Sergi Blanco-Cuaresma<fnref n="6" />,<fnref n="7" />

Received 16 July 2024 / Accepted 18 September 2024

ABSTRACT

Accurate modeling of tidal interactions is crucial for interpreting recent JWST observations of the thermal emissions of TRAPPIST-1 b and c and for characterizing the surface conditions and potential habitability of the other planets in the system. Indeed, the rotation state of the planets, driven by tidal forces, significantly influences the heat redistribution regime. Due to their proximity to their host star and the estimated age of the system, the TRAPPIST-1 planets are commonly assumed to be in a synchronization state. In this work, we present the recent implementation of the co-planar tidal torque and forces equations within the formalism of Kaula in the *N*-body code Posidonius. This enables us to explore the hypothesis of synchronization using a tidal model well suited to rocky planets. We studied the rotational state of each planet by taking into account their multi-layer internal structure computed with the code Burnman. Our simulations show that the TRAPPIST-1 planets are not perfectly synchronized but oscillate around the synchronization state. Planet-planet interactions lead to strong variations on the mean motion and tides fail to keep the spin synchronized with respect to the mean motion. As a result, the substellar point of each planet experiences short oscillations and long-timescale drifts that lead the planets to achieve a synodic day with periods varying from 55 years to 290 years depending on the planet.

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

**1. Introduction**

The first results of JWST thermal emission observations of the rocky exoplanets of the TRAPPIST-1 system (Gillon et al. 2017) have been published; they focus on the characterization of the atmosphere of TRAPPIST-1 b (Greene et al. 2023; Ih et al. 2023) and TRAPPIST-1 c (Zieba et al. 2023; Lincowski et al. 2023). Uncertainties persist about the presence and/or composition of the atmosphere of TRAPPIST-1 b and c, and further observations and data analysis are needed. In this context, a good understanding of heat redistribution processes is crucial to accurately interpreting the data for atmospheric characterization. Given that the spin of the planets drives the heat redistribution regime and the tidal heating, and considering that the rotation states of the TRAPPIST-1 planets are driven by tides due to their close orbits to the star, accurately modeling these tidal interactions is essential for interpreting observations and inferring surface conditions.

We present recent developments made in the *N*-body
code Posidonius<fnref n="1" /> (Blanco-Cuaresma & Bolmont 2017; Bolmont

et al. 2020b). We implemented the Kaula (1964) formalism, which allowed us to study the spin evolution in a more relevant way for rocky planets compared to usual simple prescriptions of tides, such as a constant time lag (CTL) or a constant phase lag model (Hut 1981; Goldreich & Soter 1966; Goldreich 1966). These models do not reproduce the correct behavior for highly viscous objects, especially for the evolution of their rotation (Henning et al. 2009; Efroimsky & Makarov 2013). In contrast, we used an Andrade rheology (Andrade 1910) and modeled the multilayer internal structure of the TRAPPIST-1 planets using the code BurnMan (Cottaar et al. 2014; Myhill et al. 2021), assuming a rocky composition. We revisited the synchronized spin state of the TRAPPIST-1 planets and, in particular, the position of their substellar point. For each planet, we computed the maximal theoretical drift due to the mean motion perturbation assuming a spin rate changing only due to the tidal lag. Section 2.1 introduces the tidal formalism used (Kaula) and presents the equations implemented in the *N*-body code Posidonius. Internal structures used for the TRAPPIST-1 planets are described in Sect. 2.2. Section 3 presents the results of the simulations, as well as the variation in the spin states around the synchronization and the drift of the substellar point of each planet.

<footnote n="1">¹ Observatoire de Genève, Université de Genève, 51 Chemin Pegasi, 1290 Sauverny, Switzerland</footnote>

<footnote n="2">² Centre Vie dans l'Univers, Faculté des sciences, Université de Genève, Quai Ernest-Ansermet 30, 1211 Genève 4, Switzerland</footnote>

<footnote n="3">³ Kapteyn Astronomical Institute, University of Groningen, Groningen, The Netherlands</footnote>

<footnote n="4">⁴ Laboratoire de Planétologie et Géosciences, UMR-CNRS 6112, Nantes Université, 2 rue de la Houssinière, BP 92208, 44322 Nantes Cedex 3, France</footnote>

<footnote n="5">⁵ naXys, Department of Mathematics, University of Namur, 61 Rue de Bruxelles, 5000 Namur, Belgium</footnote>

<footnote n="6">⁶ Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="7">⁷ Laboratoire de Recherche en Neuroimagerie, University Hospital (CHUV), Lausanne, Switzerland</footnote>

<footnote>* Corresponding authors; alexandre.revol@unige.ch, emeline.bolmont@unige.ch</footnote>

<footnote n="1">¹ https://www.blancocuaresma.com/s/posidonius</footnote>
## 2. Numerical setup

We used Posidonius (Blanco-Cuaresma & Bolmont 2017; Bolmont et al. 2020b), an N-body code that integrates the dynamical evolution of planetary systems. The code accounts for additional effects, such as tidal interactions, rotational flattening, general relativity, the evolution of the central star, and a migration prescription in a protoplanetary disk. In the original version of Posidonius, an equilibrium tide model (Bolmont et al. 2015, 2020b) and a creep tide model (Gomes et al. 2021) were implemented. For this work, we only investigated the evolution of the system through the planetary tide. The code does not take the triaxiality of planets into account. We therefore neglected the net restoring force from the permanent deformation. This assumption is discussed in Sect. 4. For this study, we implemented a new model of tides in Posidonius following the formalism of Kaula (1964). The next section presents this model and the equations implemented in Posidonius. The code used in this study is available on GitHub<fnref n="2" /> and will be incorporated into the next version of Posidonius.

### 2.1. Tidal model

We implemented the tidal force and torque equations following the formalism of Kaula (1964). This formalism is based on the expansion of the tidal perturbing gravitational potential in Fourier harmonic modes. Each mode of the expansion is associated with an excitation frequency and raises an associated tidal bulge. This makes the formalism general enough to encapsulate the frequency-dependent response of a body subjected to stress and strain (Efroimsky & Makarov 2014; Boué et al. 2019). According to Kaula (1964), the tide-raising potential can be expressed in the Fourier domain in terms of Keplerian elements as

$$W = \sum_{l=2}^{\infty} \sum_{m=0}^{l} \sum_{p=0}^{l} \sum_{q \in \mathbb{Z}} W_{lmpq}(a, e, i, \omega_{lmpq}), \quad (1)$$

The indexes l, m, p, and q are indexes of harmonic modes, a, e, and i the semi-major axis, eccentricity, and orbital inclination, respectively, and $\omega_{lmpq}$ the excitation modes, defined as

$$\omega_{lmpq} = (l - 2p + q)\dot{M} + (l - 2p)\dot{\omega} + m(\dot{\Omega} - \dot{\theta}) \quad (2)$$

with $\dot{M}$ the mean anomaly derivative, $\dot{\omega}$ the argument of periastron derivative, $\dot{\Omega}$ the longitude of node derivative and $\dot{\omega}$ the spin rate of the planet (Efroimsky & Makarov 2014). The tidal potential raised by a deformed planet can be determined from the tidal-raising potential of Eq. (1) as (we omit the dependence on a, e, i, and $\omega_{lmpq}$ in the notation for simplicity)

$$U_{lmpq} = \left(\frac{R_p}{r}\right)^{l+1} \bar{k}_l(\omega_{lmpq}) W_{lmpq}|_{lag}, \quad (3)$$

with $R_p$ the radius of the planet, r the radial distance between the star and the planet, $\bar{k}_l(\omega_{lmpq})$ the complex Love number described in Sect. 2.2, and $U_{lmpq}|_{lag}$ the tidal-raising potential of Eq. (1), where each phase is shifted by an angle $\epsilon_{lmpq}$ associated with an excitation mode $\omega_{lmpq}$. The expression of the tidal force is then derived from the gravitational potential of Eq. (3). Neglecting the obliquity and orbital inclination, the tidal force

can be expressed as

$$F = - \frac{G M_{\star}^2}{a^6} \frac{R_p^5}{r} \sum_{q=-\infty}^{+\infty} \left[ \left( \frac{3}{4} G_{21q}(e)^2 \Re(\bar{k}_2) + \frac{9}{4} G_{20q}(e)^2 \Re(\bar{k}_2) \right) e_r \right] \quad (4)$$

$$+ \frac{3}{2} G_{20q}(e)^2 \Im(\bar{k}_2) e_{\varphi}, \quad (5)$$

with $G$ the gravitational constant and $M_{\star}$ the mass of the host star. $G_{21q}(e)$ and $G_{20q}(e)$ are the eccentricity functions $G_{lmpq}(e)$ with $l, m = 2, 0$ (see Cayley 1861 and Appendix A for details). $\Re(\bar{k}_2)$ and $\Im(\bar{k}_2)$ are the real and imaginary parts of the complex Love number, $\bar{k}_l$, at degree $l = 2$ (we have omitted the frequency dependence in the formula for clarity), and $e_r$ and $e_{\varphi}$ are the radial and azimuthal unit vectors. The tidal torque applied to the planet is computed directly in the code from the tidal force as

$$N = r \times F, \quad (6)$$

with r the vector $r e_r$.

### 2.2. Love numbers and internal structures of the Trappist-1 planets

The complex Love number, $\bar{k}_l(\omega_{lmpq})$, quantifies the response of a body subjected to tidal perturbations. The frequency dependence of the Love number associated with each planet is set as an input to the code, and interpolated in the code at each timestep for every excited frequency mode $\omega_{lmpq}$. We computed the frequency-dependent Love numbers using the internal structure model described in the following section.

#### 2.2.1. Internal structure

The complex Love number can be computed for any internal structure, accounting for the density, shear modulus, and viscosity profile. We computed its frequency dependence following the method described in Dumoulin et al. (2017), Tobie et al. (2005, 2019), and Bolmont et al. (2020a). We considered a multilayer model of planets and assumed a compressible Andrade rheology (Andrade 1914). The multilayer internal structures of the TRAPPIST-1 planets are built with the BurnMan code<fnref n="3" /> (Cottaar et al. 2014; Myhill et al. 2021). The BurnMan code computes the density, temperature, gravity, and pressure profiles from the planet's surface to its center. The viscosity profile in the mantle is estimated from the temperature profile, computed with the Rayleigh number with the value at the bottom of the mantle of the Earth ($R_a = 10^6$). The viscosity of the lithosphere is set to $10^{25}$ Pa.s. Internal structures are computed from possible compositions and core sizes compatible with the mass and radius estimations of TRAPPIST-1 planets from Agol et al. (2021). The masses and radii can be reproduced with a variety of internal structures and compositions. Thus, for each planet, we assumed a silicate mantle with pyrolytic composition and a liquid metallic core of 90% iron and 2% silicate, which corresponds to the Earth proportions (Hirose et al. 2021). With the composition fixed, we varied the relative radius of the core of each planet from 50 to 60% in order to match their total mass and radius. The masses and radii are given in Table 1. Figure 1 shows the temperature, density, viscosity, and rigidity profiles for each TRAPPIST-1 planet computed with the BurnMan code.

<footnote n="2">² https://github.com/DynaClim/posidonius</footnote>

<footnote n="3">³ https://geodynamics.github.io/burnman/</footnote>
<table><caption>**Table 1.** TRAPPIST-1 system parameters used as initial conditions for the simulations.</caption><thead><tr><th>Star</th><th>Mass M<sub>☉</sub></th><th>Radius R<sub>☉</sub></th><th>-</th><th>-</th><th>-</th><th>-</th></tr><tr><th>T1- A</th><th>0.0898 M<sub>☉</sub></th><th>0.1192 R<sub>☉</sub></th><th>-</th><th>-</th><th>-</th><th>-</th></tr><tr><th>Planets</th><th>Mass M<sub>⊕</sub></th><th>Radius R<sub>⊕</sub></th><th>a (10<sup>-2</sup> AU)</th><th>e (10<sup>-3</sup>)</th><th>ω (rad)</th><th>λ (rad)</th></tr></thead><tbody><tr><td>T1-b</td><td>1.374</td><td>1.116</td><td>1.154</td><td>6.694</td><td>3.9138</td><td>6.2807</td></tr><tr><td>T1-c</td><td>1.308</td><td>1.097</td><td>1.580</td><td>1.958</td><td>3.9568</td><td>3.2752</td></tr><tr><td>T1-d</td><td>0.388</td><td>0.788</td><td>2.227</td><td>7.221</td><td>3.6462</td><td>0.2445</td></tr><tr><td>T1-e</td><td>0.692</td><td>0.920</td><td>2.927</td><td>5.215</td><td>0.9017</td><td>1.4578</td></tr><tr><td>T1-f</td><td>1.039</td><td>1.045</td><td>3.851</td><td>9.502</td><td>3.1105</td><td>1.0048</td></tr><tr><td>T1-g</td><td>1.129</td><td>1.129</td><td>4.684</td><td>3.263</td><td>5.8092</td><td>1.4551</td></tr><tr><td>T1-h</td><td>0.755</td><td>0.755</td><td>6.192</td><td>4.403</td><td>3.2359</td><td>5.0763</td></tr></tbody></table>

**Notes.** The orbital parameters, i.e. the semi-major axis *a*, the eccentricity *e*, the longitude of ascending node *ω*, and the mean longitude *λ*, are taken from Agol et al. (2021).

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Internal profiles of the TRAPPIST-1 planets computed with the BurnMan code. The panels show, from left to right, the density profile, the temperature profile, the viscosity profile, and the shear modulus profile.</figcaption></figure>

### 2.2.2. Love numbers

Love numbers were computed following the method of Tobie et al. (2005, 2019) with the code LPG-Tide, which uses the radial profile of the density, viscosity, shear modulus, and seismic velocities. We used an Andrade rheology (Andrade 1910), which defines the complex compliance $\bar{J}$ as (Castillo-Rogez et al. 2011; Efroimsky 2012)

$$ \bar{J} = J + \beta(i\omega)^{-\alpha}\Gamma(1 + \alpha) - \frac{i}{\eta\omega}, \quad (7) $$

with $\eta$ the shear viscosity, $\omega$ the excitation frequency, $\beta$ and $\alpha$ two empirical parameters, and $\Gamma$ the Gamma function. Figure 2 shows the Love numbers' frequency dependences associated with the TRAPPIST-1 profiles presented in Sect. 2.2.1. The top panel of Fig. 2 represents the frequency dependence of the imaginary part of the Love number, and thus the dissipation within each body (Tobie et al. 2005; Dumoulin et al. 2017; Bolmont et al. 2020a).

## 3. TRAPPIST-1 rotational states

We performed N-body simulations of the TRAPPIST-1 system. The rotational state was assumed to be synchronous at the beginning of the simulation. The orbital inclinations and obliquities of the planets were assumed to be zero, given the system's high coplanar configuration and high tidal damping. The initial conditions were taken from Agol et al. (2021) and are listed in Table 1. The following section presents the evolution of the rotation state of each planet, their deviation from the synchronization state, and the variation in the position of the substellar points.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Frequency dependence of the imaginary part ($\mathfrak{J}(k_2)$; top panel) and real part ($\mathfrak{R}(k_2)$; bottom panel) of the Love number of the TRAPPIST-1 planets assuming Earth-like compositions, computed from the profiles presented in Sect. 2.2.1 and Fig. 1. The vertical dotted lines correspond to each planet's principal excitation frequency (i.e., the mean motion, *n*) if assumed to be perfectly synchronous.</figcaption></figure>

### 3.1. Non-synchronization

Figure 3 shows the variation of the rotation state, in terms of $\Omega/n$ (the spin rate over the mean motion), around the synchronization state for a period of 50 days. The simulation shows that
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Evolution of the rotation state of the TRAPPIST-1 planets in terms of $\Omega/n$ (the spin rate and the mean motion, respectively) around the synchronization state ($\Omega/n = 1$).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Spin rate and mean motion of planet e over 50 days.</figcaption></figure>

none of the planets is perfectly synchronized; each oscillates around the synchronization state (i.e., $\Omega/n = 1$) with different amplitudes and timescales. These variations are due to large mean motion variations caused by planet-planet interactions. As shown in Fig. 4 for planet e, the mean motion varies with such high amplitudes and on such short timescales that it prevents the synchronization of the spin.

## 3.2. Substellar point drift

We computed the substellar position of each planet by subtracting the rotation angle ($\theta$), computed from the integration of the spin rate, from the position angle ($\nu$) of the star in the sky of the planet (the true anomaly), as $\nu-\theta$. Substellar points were computed to be zero after 50 years of simulation.

Variations of the spin states around the synchronization led to short- and long-term variations of substellar points. Figure 5 shows short-term variations, which correspond to short variations of the mean motion, as shown in Fig. 4. Amplitudes of low librations vary between 1 and 2 degrees depending on the planet, with a timescale that varies from about 1.5 days for planet b to 18 days for the planet h, which approximately corresponds to the orbital period of each planet. Figure 6 shows long-term drifts of the substellar point of each planet for a period of 250 years. The time needed for the substellar points to perform a complete rotation (i.e., the duration of the synodic day) varies between 55 years and -290 years, depending on the planet, and is summarized in Table 2. The estimations of planets T1-e, g, and h (in

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Short-term librations of substellar points of the TRAPPIST-1 planets, in degrees, for a period of 1 year. Positive angles correspond to eastward drifts and negative angles to westward drifts. An eastward (respectively westward) rotation corresponds to a prograde (respectively retrograde) rotation of the planet with respect to its orbit.</figcaption></figure>

<table><caption>**Table 2.** Estimation of the length of the day of the TRAPPIST-1 planets according to the rotation period of their substellar point.</caption><thead><tr><td>Planet</td><td>Synodic days (years)</td></tr></thead><tbody><tr><td>T1-b</td><td>-69</td></tr><tr><td>T1-c</td><td>55</td></tr><tr><td>T1-d</td><td>160</td></tr><tr><td>T1-e</td><td>-1997</td></tr><tr><td>T1-f</td><td>-290</td></tr><tr><td>T1-g</td><td>2400</td></tr><tr><td>T1-h</td><td>1103</td></tr></tbody></table>

**Notes.** Numbers in gray are estimations from linear fits and should be used with caution.

gray) should be taken with caution, as they appear to be strongly tidally locked and have not completed a single synodic rotation in 250 yr. A positive drift corresponds to a retrograde rotation with respect to its orbit, while a negative drift corresponds to a prograde rotation. We must highlight that the formalism used does not take the triaxiality of the planets into account. Therefore, considering the absence of the restoring force in the permanent deformation, these estimations correspond to the maximum theoretical drift due to the mean motion variations. Accounting for the triaxiality should lead to a stronger tidal torque (Van Hoolst et al. 2013) and, thus, reduce the drift rate. The triaxiality, which is expected to be particularly pronounced for this close-in planet due to strong permanent tidal forces, should further limit the drift. Interestingly, the two inner planets have the strongest drifts of the system despite their shorter distances to the star, and they are drifting in the opposite direction. We find that the drift is linked to the position of the planets in the resonance chain and leads to strong planet-planet interactions. Planet b also shows numerical artifacts that make the computation of the substellar point challenging for a symplectic integrator. The combination of the very short orbital period and the shape of the Love number makes the dissipation in the low-frequency regime on the dissipation spectrum (see the top panel of Fig. 2).
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Long-term variation of substellar points of the TRAPPIST-1 planets for a period of 250 years. Horizontal dotted lines represent multiples of $2\pi$ rotation.</figcaption></figure>

## 4. Conclusion

This Letter presents recent developments made in the code Posidonius (Blanco-Cuaresma & Bolmont 2017; Bolmont et al. 2020a), i.e. the implementation of the Kaula formalism, that accounts for a realistic frequency dependence of the tidal response (Kaula 1964). We applied our work to the TRAPPIST-1 system and studied for the first time the variation of the spin of each planet around the synchronization state, accounting for the forcing of the tidal bulge and only assuming a realistic tidal model. We modeled the tidal response of each planet with the rheology of Andrade (1910) following the method of Bolmont et al. (2020a) with a multi-layered internal structure computed with the internal structure code BurnMan (Cottaar et al. 2014; Myhill et al. 2021).

We find that, despite the strong tidal interactions in the system, the spin of each planet is not perfectly synchronized (i.e., the spin rate of the planets is not strictly equal to their mean motion). Planet-planet interactions significantly perturb the mean motion of planets in compact configurations which prevents spin rates from remaining synchronous and causes the planets to oscillate around the synchronization state. As a consequence of these oscillations, planets experience variations of the position of their substellar point, leading to a slow and continuous drift. Continuous drifts lead the substellar points to perform complete rotations, which lead to day-night cycles, on retrograde rotation, with synodic day periods varying from 55 to 290 years depending on the planet.

We must highlight that the triaxiality of each planet is not taken into account in this study. It is still uncertain how the long-term shape of the planet evolves with respect to the spin-orbit evolution. On the one hand, if the spin and the orbit evolve on a timescale much shorter than the interior relaxation timescale, the planet would not have time to develop a permanent long-axis deformation in the star direction, and the planet should remain axisymmetric. On the other hand, if the spin orbit evolves on a much longer timescale than the relaxation timescale, a permanent long axis would develop and counteract any planet drift. Even under this condition, a drift is still possible, and the drift rate will be controlled by the rate at which the planet shape can readjust as the spin orbit changes. As we neglected the distorted shape of the planets and their long-term evolution, the estimations we obtain here should be considered as maximum synodic days from the spin oscillations around the synchronization state.

Correia & Delisle (2019) evaluated the damping of the triaxiality due to mean motion variations but assumed a fixed permanent deformation (assuming that the triaxiality of the planet has two components, the permanent deformation due to the intrinsic mass repartition in the planet and the tidally induced deformation; see also Appendix D of Delisle et al. 2017).

Nonsynchronous spin states of TRAPPIST-1 planets have already been suggested, the first time by Makarov et al. (2018) to be spin-orbit resonances and/or pseudo-synchronization from the eccentricities of the system. However, our findings show not only eccentricities but the overall non-Keplerian shapes of the orbits due to planet-planet interactions. These interactions have been studied by Vinson et al. (2019), Shakespeare & Steffen (2023), and Chen et al. (2023), but their results differ from our findings on several points.

In particular, Vinson et al. (2019) and Shakespeare & Steffen (2023) identified several regimes of rotation coming from the pendulum effect caused by the triaxiality of planets. But firstly, they did not compute the spin rate in a full consistent way with their *N*-body simulation. Secondly, the tidal model used was the common CTL model (Hut 1981; Eggleton et al. 1998), which is not well suited for studying the rotation of rocky planets subjected to tides (Makarov & Efroimsky 2013).

Our study is the first to analyze the spin state of rocky planets while accounting for their internal structure. The drifts of substellar points depend on the competition between mean motion variation and tidal damping timescales, which are influenced by the strength of the tides and, thus, by the internal structure and the potential presence of a surface magma or water ocean. Future work will further investigate the trend of the substellar drifts and provide a detailed analysis of the impact of tidal models used and the rheological properties and internal structures of TRAPPIST-1 planets. In this study, we assumed TRAPPIST-1 planets have a metallic core and a silicate mantle, which is particularly relevant for the inner planets (b, c, and d). However, the low bulk density of the outer planets (e, f, and g) suggests they have significant volatile content and may also be compatible with the presence of water or ice (e.g., Unterborn et al. 2018; Boldog et al. 2024). Exploring a broader range of internal structures is beyond the scope of this study but will be crucial for understanding the physical and thermal processes of these exoplanets, as well as for the interpretation of observational data and assessing the potential habitability of these planets.

**Acknowledgements.** AR and EB acknowledge the financial support of the SNSF (grant number: 200021_197176 and 200020_215760). The authors thank the referee, Michael Efroimsky, who provided good criticisms that help to improve our work. This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. The computations were performed at University of Geneva on the Baobab and Yggdrasil clusters. This research has made use of NASA's Astrophysics Data System.

## References

* Agol, E., Dorn, C., Grimm, S. L., et al. 2021, PSJ, 2, 1

* Andrade, E. N. D. C. 1910, Proc. R. Soc. London Ser. A, 84, 1

* Andrade, E. N. D. C. 1914, Proc. R. Soc. London Ser. A, 90, 329

* Blanco-Cuaresma, S., & Bolmont, E. 2017, EWASS Special Session 4 (2017): Star-planet Interactions (EWASS-SS4-2017)

* Boldog, Á., Dobos, V., Kiss, L. L., van der Perk, M., & Barr, A. C. 2024, A&A, 681, A109

* Bolmont, E., Raymond, S. N., Leconte, J., Hersant, F., & Correia, A. C. M. 2015, A&A, 583, A116

* Bolmont, E., Breton, S. N., Tobie, G., et al. 2020a, A&A, 644, A165

* Bolmont, E., Demory, B. O., Blanco-Cuaresma, S., et al. 2020b, A&A, 635, A117
Boué, G., Correia, A., & Laskar, J. 2019, ApJ, 82, 91

Castillo-Rogez, J. C., Efroimsky, M., & Lainey, V. 2011, J. Geophys. Res. (Planets), 116, E09008

Cayley, A. 1861, MmRAS, 29, 191

Chen, H., Li, G., Paradise, A., & Kopparapu, R. K. 2023, ApJ, 946, L32

Correia, A. C. M., & Delisle, J.-B. 2019, A&A, 630, A102

Cottaar, S., Heister, T., Rose, I., & Unterborn, C. 2014, Geochem. Geophys. Geosyst., 15, 1164

Delisle, J. B., Correia, A. C. M., Leleu, A., & Robutel, P. 2017, A&A, 605, A37

Dumoulin, C., Tobie, G., Verhoeven, O., Rosenblatt, P., & Rambaux, N. 2017, J. Geophys. Res. (Planets), 122, 1338

Efroimsky, M. 2012, Celest. Mech. Dyn. Astron., 112, 283

Efroimsky, M., & Makarov, V. V. 2013, ApJ, 764, 26

Efroimsky, M., & Makarov, V. V. 2014, ApJ, 795, 6

Eggleton, P. P., Kiseleva, L. G., & Hut, P. 1998, ApJ, 499, 853

Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542, 456

Goldreich, P. 1966, AJ, 71, 1

Goldreich, P., & Soter, S. 1966, Icarus, 5, 375

Gomes, G. O., Bolmont, E., & Blanco-Cuaresma, S. 2021, A&A, 651, A23

Greene, T. P., Bell, T. J., Ducrot, E., et al. 2023, Nature, 618, 39

Henning, W. G., O'Connell, R. J., & Sasselov, D. D. 2009, ApJ, 707, 1000

Hirose, K., Wood, B., & Vočadlo, L. 2021, Nat. Rev. Earth & Environ., 2, 645

Hut, P. 1981, A&A, 99, 126

Ih, J., Kempton, E. M. R., Whittaker, E. A., & Lessard, M. 2023, ApJ, 952, L4

Kaula, W. M. 1961, Geophys. J., 5, 104

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

Lincowski, A. P., Meadows, V. S., Zieba, S., et al. 2023, arXiv e-prints [arXiv:2308.05899]

Makarov, V. V., & Efroimsky, M. 2013, ApJ, 764, 27

Makarov, V. V., Berghea, C. T., & Efroimsky, M. 2018, ApJ, 857, 142

Myhill, R., Cottaar, S., Heister, T., Rose, I., & Unterborn, C. 2021, https://doi.org/10.5281/zenodo.5552756

Shakespeare, C. J., & Steffen, J. H. 2023, Bull. Am. Astron. Soc., 54, 203.04

Tisserand, F. 1889, Traité de mécanique céleste: Tome I, Perturbations des planètes d'après la méthode de la variation des constantes arbitraires (Gauthier-Villars)

Tobie, G., Mocquet, A., & Sotin, C. 2005, Icarus, 177, 534

Tobie, G., Grasset, O., Dumoulin, C., & Mocquet, A. 2019, A&A, 630, A70

Unterborn, C. T., Desch, S. J., Hinkel, N. R., & Lorenzo, A. 2018, Nat. Astron., 2, 297

Van Hoolst, T., Baland, R.-M., & Trinh, A. 2013, Icarus, 226, 299

Vinson, A. M., Tamayo, D., & Hansen, B. M. S. 2019, MNRAS, 488, 5739

Zieba, S., Kreidberg, L., Ducrot, E., et al. 2023, arXiv e-prints [arXiv:2306.10150]
# Appendix A: Table of the inclination and eccentricity polynomials

The eccentricity polynomials $G_{lpq}(e)$ are elliptic expansions given by Eqs. 23 and 24 of Kaula (1961) and can be computed from the Hansen functions $X_{l-2p+q}^{-(l-1),(l-2p)}$ (Cayley 1861; Tisserand 1889). They are presented in the Table A.1.

We considered eccentricities up to 0.3, which allowed us to study the eccentricity expansions up to order 7 (see the tables in Cayley 1861). As the tidal interactions are computed at the quadrupolar order $l=2$, the index $p$ is constrained between 0 and 2 and $q$ between -7 and 7. The elements $G_{22q}(e)$ can be computed with the symmetry rule $G_{2,2,q}(e) = G_{2,0,-q}(e)$.

<table><caption>**Table A.1.** Eccentricity polynomials, $G_{lpq}(e)$, from Cayley (1861), up to order 7 in eccentricity.</caption><thead><tr><th>l</th><th>p</th><th>q</th><th>G<sub>lpq</sub>(e)</th></tr></thead><tbody><tr><td>2</td><td>0</td><td>-7</td><td><sup>15625</sup>&frasl;<sub>129024</sub> e<sup>7</sup></td></tr><tr><td>2</td><td>0</td><td>-6</td><td><sup>4</sup>&frasl;<sub>45</sub> e<sup>6</sup></td></tr><tr><td>2</td><td>0</td><td>-5</td><td><sup>81</sup>&frasl;<sub>1280</sub> e<sup>5</sup> + <sup>81</sup>&frasl;<sub>2048</sub> e<sup>7</sup></td></tr><tr><td>2</td><td>0</td><td>-4</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>2</td><td>0</td><td>-3</td><td><sup>1</sup>&frasl;<sub>48</sub> e<sup>3</sup> + <sup>11</sup>&frasl;<sub>768</sub> e<sup>5</sup> + <sup>313</sup>&frasl;<sub>30720</sub> e<sup>7</sup></td></tr><tr><td>2</td><td>0</td><td>-2</td><td>0</td></tr><tr><td>2</td><td>0</td><td>-1</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> - <sup>143</sup>&frasl;<sub>18432</sub> e<sup>7</sup></td></tr><tr><td>2</td><td>0</td><td>0</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>2</td><td>0</td><td>1</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> - <sup>1763</sup>&frasl;<sub>2048</sub> e<sup>7</sup></td></tr><tr><td>2</td><td>0</td><td>2</td><td><sup>17</sup>&frasl;<sub>7</sub> e<sup>2</sup> - <sup>115</sup>&frasl;<sub>16</sub> e<sup>4</sup> + <sup>601</sup>&frasl;<sub>48</sub> e<sup>6</sup></td></tr><tr><td>2</td><td>0</td><td>3</td><td><sup>845</sup>&frasl;<sub>48</sub> e<sup>3</sup> - <sup>32525</sup>&frasl;<sub>16</sub> e<sup>5</sup> + <sup>208225</sup>&frasl;<sub>6144</sub> e<sup>7</sup></td></tr><tr><td>2</td><td>0</td><td>4</td><td><sup>333</sup>&frasl;<sub>16</sub> e<sup>4</sup> - <sup>13827</sup>&frasl;<sub>16</sub> e<sup>6</sup></td></tr><tr><td>2</td><td>0</td><td>5</td><td><sup>228347</sup>&frasl;<sub>3840</sub> e<sup>5</sup> - <sup>9071075</sup>&frasl;<sub>18432</sub> e<sup>7</sup></td></tr><tr><td>2</td><td>0</td><td>6</td><td><sup>73369</sup>&frasl;<sub>720</sub> e<sup>6</sup></td></tr><tr><td>2</td><td>0</td><td>7</td><td><sup>12144273</sup>&frasl;<sub>71680</sub> e<sup>7</sup></td></tr><tr><td>2</td><td>1</td><td>-7</td><td>G<sub>217</sub>(e)</td></tr><tr><td>2</td><td>1</td><td>-6</td><td>G<sub>216</sub>(e)</td></tr><tr><td>2</td><td>1</td><td>-5</td><td>G<sub>215</sub>(e)</td></tr><tr><td>2</td><td>1</td><td>-4</td><td>G<sub>214</sub>(e)</td></tr><tr><td>2</td><td>1</td><td>-3</td><td>G<sub>213</sub>(e)</td></tr><tr><td>2</td><td>1</td><td>-2</td><td>G<sub>212</sub>(e)</td></tr><tr><td>2</td><td>1</td><td>-1</td><td>G<sub>211</sub>(e)</td></tr><tr><td>2</td><td>1</td><td>0</td><td>(1 - e<sup>2</sup>)<sup>-3/2</sup> ≈ 1 + <sup>3</sup>&frasl;<sub>8</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> + O(e<sup>9</sup>)</td></tr><tr><td>2</td><td>1</td><td>1</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> + <sup>14309</sup>&frasl;<sub>6144</sub> e<sup>7</sup></td></tr><tr><td>2</td><td>1</td><td>2</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></tr><tr><td>2</td><td>1</td><td>3</td><td><sup>53</sup>&frasl;<sub>16</sub> e<sup>3</sup> + <sup>393</sup>&frasl;<sub>256</sub> e<sup>5</sup> + <sup>24753</sup>&frasl;<sub>10240</sub> e<sup>7</sup></td></tr><tr><td>2</td><td>1</td><td>4</td><td><sup>77</sup>&frasl;<sub>16</sub> e<sup>4</sup> + <sup>125</sup>&frasl;<sub>16</sub> e<sup>6</sup></td></tr><tr><td>2</td><td>1</td><td>5</td><td><sup>16</sup>&frasl;<sub>73</sub> e<sup>5</sup> - <sup>4987</sup>&frasl;<sub>6144</sub> e<sup>7</sup></td></tr><tr><td>2</td><td>1</td><td>6</td><td><sup>3167</sup>&frasl;<sub>320</sub> e<sup>6</sup></td></tr><tr><td>2</td><td>1</td><td>7</td><td><sup>432091</sup>&frasl;<sub>30720</sub> e<sup>7</sup></td></tr></tbody></table>

In practice, the order of the summation over $q$ is taken as a function of the value of the eccentricity. The maximum order considered is listed in Table A.2.

<table><caption>**Table A.2.** Maximum order of the summation over $q$ as a function of the eccentricity.</caption><tr><td>e</td><td>0</td><td>&lt;0.05</td><td>&lt;0.1</td><td>&lt;0.15</td><td>&lt;0.2</td><td>&lt;0.25</td><td>&lt;0.3</td><td>&gt;0.3</td></tr><tr><td>|q|<sub>max</sub></td><td>0</td><td>1</td><td>2</td><td>3</td><td>4</td><td>5</td><td>6</td><td>7</td></tr></table>