© ESO 2020

1 Introduction

The measurement of transit timing variations (TTVs) in the context of multi-transiting planet systems can be a very efficient method for deriving the dynamical parameters of a planetary system, such as mass and eccentricity (see Agol & Fabrycky 2018, for a review). The TRAPPIST-1 system has been intensely monitored by TRAPPIST, K2, and Spitzer, which led to estimates of the masses of the planets by Grimm et al. (2018). Recently additional Spitzer observations were obtained thanks to the Spitzer proposal no. 14223 (Agol et al. 2019).

In most studies on TTVs the model used is an N-body model assuming point-mass and Newtonian dynamics, sometimes with a prescription for general relativity (as in Grimm et al. 2018; Jordán & Bakos 2008; Pál & Kocsis 2008). Theoretical studies have considered the possible impact of tides and quadrupole distortion on transit times (Miralda-Escudé 2002; Heyl & Gladman 2007). However, the influence of tides has never been consistently taken into account in a multi-planet context.

Some studies do take into account tidal decay (e.g., Maciejewski et al. 2018), but decay typically occurs on time- scales that are much longer than the typical duration of the observations available for TRAPPIST-1. However, tidal forces are not only a dissipative effect (which drives migration and spin evolution), there is also a non-dissipative effect which depends on the real part of the complex Love number of degree-2, k₂, which quantifies the shape of the tidal deformation. This deformation can lead to a precession of the orbit, which can lead to TTVs. In addition, for fast rotating planets, rotational flattening can also drive a precession of the orbit, which in turn can lead to TTVs. These effects have been considered in systems with a single hot-Jupiter planet (see Ragozzine & Wolf 2009, for a comparative study of each effect). The precession of the orbit leads to observable TTVs, which then can inform the internal structure of the planet through the determination of the Love number.

However, these effects are usually never taken into account when investigating the TTVs of multi-planetary systems. We show that in the context of TRAPPIST-1 (Gillon et al. 2017, 2016), the inclusion of tidal forces may lead to an observable TTV signal. In contrast with TRAPPIST-1 d–h, planets b and c are in proximity to a higher-order resonance (increasing the frequency of the TTV pattern modulation; see Agol et al. 2005) and exhibit small TTV amplitudes (2–5 min); both effects inflate the uncertainties on the masses and eccentricities, as shown in Grimm et al. (2018). Interestingly, a recent TTV analysis using data obtained in the past two years yields a 34% increase in mass for TRAPPIST-1b and a 13% increase in mass for TRAPPIST-1c (in prep.) compared to Grimm et al. (2018). These mass increases of the two inner planets drew our attention to physical processes that could impact the planet physical and orbital parameters on secular timescales. As the parameters for the other planets have remained relatively insensitive to the addition of new data, two hypotheses remain. A first possibility is that the changing masses are due to an incomplete sampling of the TTV pattern that should resolve as new data are included. A second possibility is that dynamical models are missing physical processes that impact the close-in planets more strongly, such as tides and rotational flattening.

We show in this paper that the precession caused by general relativity, by tidal deformation, and by rotational flattening could lead to significantly different TTVs for the two inner planets of TRAPPIST-1.

2 Simulation setup

We use POSIDONIUS¹ v2019.07.30 (Blanco-Cuaresma & Bolmont 2017; Blanco-Cuaresma & Bolmont, in prep.), an N-body code which allows users to take into account additional forces and torques: tidal forces and torques, rotational flattening forces and torques, and general relativity (Bolmont et al. 2015). As in Bolmont et al. (2015), tides are computed between a planet and the star independently of the other planets, and the planet–planet tides are not taken into account (which is justified, see Hay & Matsuyama 2019). In Posidonius, we use the integrator IAS15 (Rein & Spiegel 2015) to compute the evolution of the system for 1500 days, which is approximately the time range available from all the observations collected from the system, and we fix the maximum time step allowed at 0.01 day = 14 min. We tested the convergence of our code with time steps of 0.005 day and 0.001 day, and find that the transit timings are stable to a precision of better than 10⁻⁶ s.

2.1 Tidal model

Posidonius accounts for equilibrium tides following the prescription of Bolmont et al. (2015), which is an implementation of the constant time-lag model (Mignard 1979; Hut 1981; Eggleton et al. 1998). The equilibrium tide is the result of the hydrostatic adjustment of a body; instead, the dynamical tide is the tidal response corresponding to the propagation of waves, for example inertial waves in the convective region of stars (see Zahn 1975) or gravito-inertial waves in a planetary liquid layer (see Auclair-Desrotour et al. 2019).

We review here the expressions for the tidal force and torques. Let us consider a star, defined by its mass M_⋆, its radius R_⋆, its degree-2 potential Love number k_2,⋆, its (constant) time lag Δτ_⋆, and its spin vector Ω_⋆. Let us consider one planet, j, orbiting the star at a distance r_j. The planet is defined by its mass $M_{{\text{p}}_{\text{j}}}$, its radius $R_{{\text{p}}_{\text{j}}}$, its degree-2 potential Love number $k_{2,{\text{p}}_{\text{j}}}$, its (constant) time lag $\Delta {\tau }_{{\text{p}}_{\text{j}}}$, and its spin vector ${\Omega }_{{\text{p}}_{\text{j}}}$.

Let us define $F_{\text{diss},{\text{p}}_{\text{j}}}$ and $F_{\text{nodiss},{\text{p}}_{\text{j}}}$ as the dissipative part and the non-dissipative part, respectively, of the force exerted on planet j due to the planetary tide as (Bolmont et al. 2015) $$\begin{array}{l}{F}_{\text{nodiss},{\text{p}}_{\text{j}}}=\frac{-3\mathcal{G}}{r_{\text{j}}{}^7}{M}_{\star }{}^2{k}_{2,{\text{p}}_{\text{j}}}{R}_{{\text{p}}_{\text{j}}}{}^5{e}_{r_{\text{j}}},\\ F_{\text{diss},{\text{p}}_{\text{j}}}=-9\mathcal{G}\frac{\stackrel{.}{r_{\text{j}}}}{r_{\text{j}}{}^8}{M}_{\star }{}^2{R}_{{\text{p}}_{\text{j}}}{}^5{k}_{2,{\text{p}}_{\text{j}}}\Delta {\tau }_{{\text{p}}_{\text{j}}}{e}_{r_{\text{j}}}\\ +3\mathcal{G}\frac{M_{\star }{}^2{R}_{{\text{p}}_{\text{j}}}{}^5}{r_{\text{j}}{}^7}{k}_{2,{\text{p}}_{\text{j}}}\Delta {\tau }_{{\text{p}}_{\text{j}}}\left({\Omega }_{{\text{p}}_{\text{j}}}-\stackrel{.}{{\theta }_{\text{j}}}\right)\times e_{r_{\text{j}}},\end{array}$$

where $e_{r_{\text{j}}}$ is the unit vector r_j∕r_j and $\stackrel{.}{{\theta }_{\text{j}}}$ is a vector collinear with the orbital angular momentum of planet j, the norm of which is equal to the time derivative of the true anomaly. Let us define F_diss,⋆ and F_nodiss,⋆ as the dissipative part and the non-dissipative part, respectively, of the force exerted on planet j due to the stellar tide as $$\begin{array}{l}{F}_{\text{nodiss},\star }=\frac{-3\mathcal{G}}{r_{\text{j}}{}^7}{M}_{{\text{p}}_{\text{j}}}{}^2{k}_{2,\star }{R}_{\star }{}^5{e}_{r_{\text{j}}},\\ F_{\text{diss},\star }=-9\mathcal{G}\frac{\stackrel{.}{r_{\text{j}}}}{r_{\text{j}}{}^8}{M}_{{\text{p}}_{\text{j}}}{}^2{R}_{\star }{}^5{k}_{2,\star }\Delta {\tau }_{\star }{e}_{r_{\text{j}}}\\ +3\mathcal{G}\frac{M_{{\text{p}}_{\text{j}}}{}^2{R}_{\star }{}^5}{r_{\text{j}}{}^7}{k}_{2,\star }\Delta {\tau }_{\star }\left({\Omega }_{\star }-\stackrel{.}{{\theta }_{\text{j}}}\right)\times e_{r_{\text{j}}}.\end{array}$$

The total force as a result of the tides acting on a planet j is therefore given by the sum of these contributions (Bolmont et al. 2015) $$\begin{array}{ll}{F}_{{\text{p}}_{\text{j}}}^T\hfill & =F_{\text{diss},{\text{p}}_{\text{j}}}+F_{\text{nodiss},{\text{p}}_{\text{j}}}+F_{\text{diss},\star }+F_{\text{nodiss},\star }.\hfill \end{array}$$(5)

2.2 Rotational flattening model

To account for rotational flattening, we also follow here the prescription of Bolmont et al. (2015), which assumes that the deformation due to the rotational flattening results in a triaxial ellipsoid symmetric with respect to the rotation axis (Murray & Dermott 1999). This deformation is quantified by a parameter, J₂, which depends on the radius, mass, and spin of the body and on the potential Love number of degree-2 for a perfectly fluid body (which we call here the fluid Love number; Correia & Rodríguez 2013). We define this parameter for planet j and the star as $$\begin{array}{ll}{J}_{2,{\text{p}}_{\text{j}}}\hfill & =k_{2f,{\text{p}}_{\text{j}}}\frac{{\Omega }_{{\text{p}}_{\text{j}}}{}^2{R}_{{\text{p}}_{\text{j}}}{}^3}{3\mathcal{G}{M}_{{\text{p}}_{\text{j}}}},\hfill \\ J_{2,\star }\hfill & =k_{2f,\star }\frac{{\Omega }_{\star }{}^2{R}_{\star }{}^3}{3\mathcal{G}{M}_{\star }}.\hfill \end{array}$$

Let us define $F_{\text{rot},{\text{p}}_{\text{j}}}$ the force exerted on planet j due to the rotational flattening of planet j and F_rot,⋆ the force exerted on planet j due to the rotational flattening of the star as (Murray & Dermott 1999; Correia et al. 2011) $$\begin{array}{ll}{F}_{\text{rot},{\text{p}}_{\text{j}}}\hfill & =\left(-\frac{3}{r_{\text{j}}{}^5}{C}_{{\text{p}}_{\text{j}}}+\frac{15}{r_{\text{j}}{}^7}{C}_{{\text{p}}_{\text{j}}}\frac{{\left(r_{\text{j}}.{\Omega }_{{\text{p}}_{\text{j}}}\right)}^2}{{\Omega }_{{\text{p}}_{\text{j}}}{}^2}\right)r_{\text{j}}-\frac{6}{r_{\text{j}}{}^5}{C}_{{\text{p}}_{\text{j}}}\frac{r_{\text{j}}.{\Omega }_{{\text{p}}_{\text{j}}}}{{\Omega }_{{\text{p}}_{\text{j}}}{}^2}{\Omega }_{{\text{p}}_{\text{j}}},\hfill \\ F_{\text{rot},\star }\hfill & =\left(-\frac{3}{r_{\text{j}}{}^5}{C}_{\star }+\frac{15}{r_{\text{j}}{}^7}{C}_{\star }\frac{{\left(r_{\text{j}}.{\Omega }_{\star }\right)}^2}{{\Omega }_{\star }{}^2}\right)r_{\text{j}}-\frac{6}{r_{\text{j}}{}^5}{C}_{\star }\frac{r_{\text{j}}.{\Omega }_{\star }}{{\Omega }_{\star }{}^2}{\Omega }_{\star },\hfill \end{array}$$

where C_⋆ and $C_{{\text{p}}_{\text{j}}}$ are defined as follows: $$\begin{array}{ll}{C}_{{\text{p}}_{\text{j}}}\hfill & =\frac{1}{2}\mathcal{G}{M}_{{\text{p}}_{\text{j}}}{M}_{\star }{J}_{2,{\text{p}}_{\text{j}}}{R}_{{\text{p}}_{\text{j}}}{}^2,\hfill \\ C_{\star }\hfill & =\frac{1}{2}\mathcal{G}{M}_{{\text{p}}_{\text{j}}}{M}_{\star }{J}_{2,\star }{R}_{\star }{}^2.\hfill \end{array}$$

The resulting force on planet j due to the rotational deformation of both the star and planet j is the sum of the two contributions: $$F_{{\text{p}}_{\text{j}}}^R=F_{\text{rot},{\text{p}}_{\text{j}}}+F_{\text{rot},\star }.$$(12)

2.3 General relativity

We use three different prescriptions for general relativity: Kidder (1995), which is the one used in Mercury-T (Bolmont et al. 2015; Anderson et al. 1975; and Newhall et al. 1983).

The prescription of Kidder (1995) was designed for two bodies, and Posidonius takes into account the post-Newtonian, the spin-orbit, and the post²-Newtonian contributions to the total acceleration (Eqs. (2.2b)–(2.2d), respectively), as well as the spin precession equations for both bodies (Eqs. (2.4a) and (2.4b)).

The prescription of Anderson et al. (1975) accounts for the post-Newtonian acceleration of two bodies. We refer the reader to Eq. (12) in Anderson et al. (1975) where the expression of this acceleration is given.

The prescription of Newhall et al. (1983) is more complete insofar as it accounts for the post-Newtonian effect between all bodies. We refer to Eq. (1) in Newhall et al. (1983), which gives the point-mass acceleration. Posidonius accounts for this acceleration, except for the last term which accounts for the perturbation of five solar system asteroids.

3 Transit timing variations

We performed simulations of the TRAPPIST-1 system switching on and off these various effects: the effect of the planetary tide (by varying the Love number k_2,p and the time lag τ_p), the effect of the stellar tide (by varying the Love number k_2,⋆), the effect of the rotational flattening of the planets (by varying the fluid Love number k_2f,p), the effect of the rotational flattening of the star (by varying the fluid Love number k_2f,⋆), and the effect of general relativity. We list in Table 1 the reference values of the parameters varied here, and we refer the reader to Appendix A for the parameters that remain constant in our simulations, including the initial orbital elements for the planets. We tested the three different prescriptions of the general relativity introduced in Sect. 2.3. They gave very similar results so that in the following we compare the other effects with respect to the simulations performed using the prescription of Kidder (1995).

The planetary reference values were taken to be representative of the Earth; in particular, the quantity k_2,p Δτ_p is equal to 213 s (Neron de Surgy & Laskar 1997). The stellar reference values were chosen to be representative of fully convective M dwarfs (Bolmont et al. 2015).

For all the simulations we performed, we calculated the transit timing variations (Agol & Fabrycky 2018) as follows: (i) for each transit we find the time of the transit mid-time by performing an interpolation to find the precise time a given planet crosses a reference direction. This corresponds to the “Observed transit time” O; (ii) we evaluate the “Calculated transit times” C by performing a linear fit² of the transit times calculated in step i) over the total number of transits; (iii) we calculate the difference O–C to obtain the TTVs as a function of the epoch (or transit number).

To quantify the impact of each additional effect on the simulated TTVs, we compute the difference between the TTVs calculated taking into account an additional effect and the TTVs obtained for a Newtonian N-body integration.

4 Influence of each effect on the TTVs

We performed a set of six simulations of the TRAPPIST-1 system to test the impact of the additional effects listed in Sect. 2, and we compare each with a Newtonian N-body simulation. One after another, we explored the effect of each parameter using the reference values in Table 1and general relativity. We tested the influence of the dissipative part of the planetary tide by assuming $\Delta {\tau }_{\text{p}}=1\times \Delta {\tau }_{\text{p}}{}^{\text{ref}}$ for all planets, with all the other parameters set to zero. We repeated the operation for the non-dissipative part of the planetary tide (through k_2,p, equal for all planets), the rotational flattening of the planets (through k_2f,p, equal for all planets), the non-dissipative part of the stellar tide (through k_2,⋆), the rotational flattening of the star (through k_2f,⋆), and for general relativity.

Figure 1 shows the results for planets b to d. The top panels of Fig. 1 show the transit timing variations for the three planets for the seven simulations, and the bottom panels show the difference between the TTVs and the TTVs corresponding to the pure N-body simulation. The different additional effects have a very limited impact on the shape of the TTVs, but computing the difference with the pure N-body case reveals the amplitude of each effect.

For TRAPPIST-1b (T-1b), the dominant effects are the non-dissipative part of the planetary tide (green in Fig. 1a) and general relativity (pink), respectively accounting for a difference in TTVs of about −0.63 s and 0.56 s after 1500 days. The effect of the rotational flattening of the planets (red) plays a smaller role, but still accounts formore than half the amplitude due to the non-dissipative part of the planetary tide with a difference of −0.33 s. The effect of the dissipative part of the planetary tide (orange) and the non-dissipative part of the stellar tide (purple) are completely negligible (accounting for a difference ~1 × 10⁻³ and ~ 1 ×10⁻⁴ s, respectively), which is in agreement with Ragozzine & Wolf (2009). The effect of the rotational flattening of the star (brown) is much smaller (accounting for a difference of −0.057 s), but might contribute to a lesser extent. Accounting for all effects (gray) leads to an absolute difference of −0.45 s at the end of the 1500-day simulation. The effects of general relativity and the non-dissipative part of the planetary tide almost cancel each other out, while the amplitude is determined by the effect of the rotational flattening of the planet (red curve) and of the star (brown curve).

For TRAPPIST-1c (T-1c), the dominant effects are the non-dissipative part of the planetary tide (accounting for a difference of 1 s, in green in Fig.1b), followed by the effect of the rotational flattening of the planet (0.53 s, in red), followed by the effect of general relativity (−0.44 s, pink). The effect of the rotational flattening of the star accounts for 0.067 s (purple) and the effect of the dissipative part of the planetary tide remains negligible (0.013 s, orange). As with T-1b, to reproduce the difference observed when all effects are taken into account (in gray in Fig.1b), we need to account for the non-dissipative part of the planetary tide, the rotational flattening of the planet, general relativity, and the rotational flattening of the star to a lesser extent. We note that the precession of the orbits due to the rotational flattening depends on the square of the spin frequency of the considered body (Ragozzine & Wolf 2009). Here we use a rotation period of 3.3 days for TRAPPIST-1 (Luger et al. 2017). It is possible that the rotation is slower (as the period distribution of nearby late M dwarfs shows, Newton et al. 2016), in which case the contribution of the rotational flattening of the star would be even less important.

For TRAPPIST-1d (T-1d), Fig. 1c shows that the dominant effect is general relativity (−2.03 s, in pink). The effect of the non-dissipative part of the planetary tide (accounting for −0.81 s, green) and the effect of the rotational flattening of the planet (accounting for −0.43 s, red) should probably also be taken into account. This is also true for all the external planets: the amplitude due to general relativity is much higher, but at the same time not accounting for at least the non-dissipative part of the planetary tide and the rotational flattening of the planet leads to small offsets (see Appendix B, Fig B.1).

Fig. 1 Impact of various additional effects on the TTVs of (a) planet b, (b) planet c, and (c) planet d. Top panel: transit timing variations for a pure N-body simulation (blue); for a simulation using $\Delta {\tau }_{\text{p}}=1\times \Delta {\tau }_{\text{p}}{}^{\text{ref}}$ for all planets (orange), using $k_{2,\text{p}}=1\times k_{2,\text{p}}{}^{\text{ref}}$ for all planets (green), using $k_{2f,\text{p}}=1\times k_{2f,\text{p}}{}^{\text{ref}}$ for all planets (red), using $k_{2,\star }=1\times k_{2,\star }{}^{\text{ref}}$ (purple), using $k_{2f,\star }=1\times k_{2f,\star }{}^{\text{ref}}$ (brown), and where we only consider the general relativity (pink); and for a simulation where all effects are taken into account. Bottom panel: corresponding TTVs differences with the Newtonian N-body case.

5 Potential observable effects

We performed simulations for which we varied the potential Love number k_2,p and the fluid Love number k_2f,p over a wide range. We first treat these parameters as free parameters with no limitations on their value and then we discuss the validity of this approach in Sect. 6.

As in Sect. 4, we always assume the same value of the Love numbers for all planets. We vary the parameter k_2,p from $1\times k_{2,\text{p}}{}^{\text{ref}}=0.299$ to $50\times k_{2,\text{p}}{}^{\text{ref}}=14.95$; the impact on the TTV differences with the pure N-body case can be seen in the top panel of Fig. 2a for T-1b and in the top panel of Fig. 2b for T-1c.

Considering the highest value of the tidal Love number leads to a difference in TTVs after 1500 days of −31.7 s for T-1b and 50.7 s for T-1c. The amplitude of these effects are then comparable to the precision achievable today on the observed TTVs of the two inner planets. If the tidal Love number could reach such high values, the effect of the non-dissipative part of the planetary tide (the tidal deformation of T-1b and T-1c) would therefore be detectable.

Similarly, we vary the parameter k_2f,p from $1\times k_{2f,\text{p}}{}^{\text{ref}}=0.9532$ to $50\times k_{2f,\text{p}}{}^{\text{ref}}=47.66$, and the impact on the TTV differences with the pure N-body case can be seen in the bottom panel of Fig. 2a for T-1b and in the bottom panel of Fig. 2b for T-1c. Assuming the highest value of the fluid Love number leads to a difference inTTVs after 1500 days of -17.3 s for T-1b and 27.7 s for T-1c.

Considering all the important effects given in Sect. 4 and relaxing the range of possible Love numbers might therefore be a way to settle the question of the increasing masses of T-1b and T-1c and to settle the two hypotheses given in the introduction: is it a sampling problem or are we missing dynamical processes? The answer to this question depends upon the potential degeneracy of these effects with varying the N-body parameters, and on the duration and precision of the transit timing measurements. We do not explore these effects here.

6 Discussion

We showed that for systems like TRAPPIST-1 the effect of the tidal deformation of the planets (through the planetary tidal Love number k_2,p), the effect of the rotational deformation of the planets (through the planetary fluid Love number k_2f,p), and the effect of the rotational flattening of the star (through the stellar fluid Love number k_2f,⋆) can impact the TTVs of the two inner planets at the same order of magnitude as the general relativity, if we assume “standard” values for these parameters. We also showed that the tidal dissipation (responsible for the misalignment which drives long term tidal evolution) does not significantly impact the TTVs of the system over the short observation time that we simulated. By relaxing the assumptions on the planetary tidal and fluid Love numbers, we also showed that a high Love number can lead to differences in TTV on the order of ~ 10 s. This difference is potentially observable with the current precision we have on the transit timings, unless there is significant degeneracy with other N-body parameters.

However, it is commonly accepted that a tidal Love number cannot exceed 1.5, which corresponds to a homogeneous body. This means that the physical range of our study should encompass maximum values, which are $5\times k_{2,\text{p}}{}^{\text{ref}}=5\times 0.299$. Limiting ourselves to this value would entail a difference in TTVs for T-1b of less than 2.5 s, which is below the precision we can achieve today. On the other hand, it is known that if a planet has a liquid layer (liquid water ocean, or liquid magma ocean), the response of the body becomes more complex; in particular, it becomes highly dependent on the excitation frequency. Specifically, if a frequency excites a resonant mode of the ocean, the tidal response can be much higher than a homogeneous rocky-planet model would predict (see, e.g., Auclair-Desrotour et al. 2019). Investigating this aspect consistently will require us to generalize the tidal formalism used here to account for the frequency dependence of the dynamical tide (e.g., using the formalism of Kaula 1961).

That is why we think that we might need to perform a TTV analysis of the TRAPPIST-1 system accounting for the various physical processes described here with no particular preconception about the values of the parameters for the planetary Love numbers. If the TTVs are reproduced by having a T-1b planet with a high Love number, this could be a sign for a liquid layer on the planet, possibly a magma ocean given the flux it receives and the tidal heat flux it might generate (e.g., Turbet et al. 2018; Makarov et al. 2018).

While a difference of a few ~10 s is potentially observable, it could be interpreted as a system with slightly different planetary masses and periods by a classical TTV retrieval code. Our group is thus currently investigating if these effects could be picked up with such a retrieval code and, if so, under what conditions (e.g., duration of the observations, precision of the timings). We are also working on implementing these effects in the TTV analysis pipeline and plan to revisit the analysis with the additional parameters mentioned earlier (Grimm et al., in prep.).

Fig. 2 Differences in TTVs with the pure N-body case for (a) planet b and (b) planet c. Top panel: fixing the planetary fluid Love number and the dissipation to their reference values (see Table 1), the potential Love number k2,p is varied between 1 and 50$\times k_{2,\text{p}}{}^{\text{ref}}~15$. Bottom panel: fixing the planetary Love number and the dissipation to their reference values, the potential Love number k2f,p is varied between 1 and 50$\times k_{2f,\text{p}}{}^{\text{ref}}~48$.

Acknowledgements

The authors would like to thank the anonymous referee for helping improving the manuscript. This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation. B.-O.D. acknowledges support from the Swiss National Science Foundation (PP00P2_163967). P.A.-D. acknowledges financial support from the European Research Council via the Consolidator grant EXOKLEIN (grant number 771620). F.S. acknowledges support from the CNRS/INSU PNP (Programme National de Planétologie). This research has made use of NASA’s Astrophysics Data System.

Appendix A Initial conditions for the TRAPPIST-1 simulations

To ensure the reproducibility of our simulations, we give here the exact initial conditions we took for the system. Table A.1 gives the stellar parameters used for the integration of the system. The stellar mass and radius come from Gillon et al. (2017) and the rotation comes from Luger et al. (2017). The value of the radius of gyration squared $r{g}_{\star }^2$ (Hut 1981) comes from Bolmont et al. (2015) and should be typical of a fully convective dwarf.

Table A.2 gives the masses and radii of the planets as well as the initial orbital elements. We consider that all planets have the same radius of gyration squared $r{g}_{\text{p}}^2=0.3308$ (where this quantity is related to the moment of inertia $I_{\text{p}}=M_{\text{p}}{(R_{\text{p}}r{g}_{\text{p}})}^2$). We consider that all planets have a zero obliquity (angle between the direction perpendicular to the orbital plane and the rotation axis of the planet) and that they are tidally locked (see discussion in Luger et al. 2017).

To perform the integration of the system, we used POSIDONIUS v2019.07.30³. This versionwas slightly altered to be able to fix a maximum time step size (0.01 day). The initial conditions can be found in Bolmont et al. (2020)⁴.

Appendix B Transit timing variations for the four outer planets of TRAPPIST-1

As in Figs. 1 and B.1 shows the difference in TTVs with the N-body case for seven different simulations (N-body and simulations with additional effects) for TRAPPIST-1e to TRAPPIST-1h. General relativity is the dominant effect, but the non-dissipative part of the planetary tidal force (via k_2,p) and the rotational flattening of the planets (via k_2f,p) are still contributing marginally. Only for TRAPPIST-1h is general relativity the only relevant process to account for.

Fig. B.1 As the bottom panels of Fig. 1, but for (a) TRAPPIST-1e, (b) TRAPPIST-1f, (c) TRAPPIST-1g, and (d) TRAPPIST-1h. The general relativity is the dominant effect, but the planetary deformation (due to tides or rotation) is not quite completely negligible, except for TRAPPIST-1h.

References

All Tables

All Figures

Fig. 1 Impact of various additional effects on the TTVs of (a) planet b, (b) planet c, and (c) planet d. Top panel: transit timing variations for a pure N-body simulation (blue); for a simulation using $\Delta {\tau }_{\text{p}}=1\times \Delta {\tau }_{\text{p}}{}^{\text{ref}}$ for all planets (orange), using $k_{2,\text{p}}=1\times k_{2,\text{p}}{}^{\text{ref}}$ for all planets (green), using $k_{2f,\text{p}}=1\times k_{2f,\text{p}}{}^{\text{ref}}$ for all planets (red), using $k_{2,\star }=1\times k_{2,\star }{}^{\text{ref}}$ (purple), using $k_{2f,\star }=1\times k_{2f,\star }{}^{\text{ref}}$ (brown), and where we only consider the general relativity (pink); and for a simulation where all effects are taken into account. Bottom panel: corresponding TTVs differences with the Newtonian N-body case.

In the text

Fig. 2 Differences in TTVs with the pure N-body case for (a) planet b and (b) planet c. Top panel: fixing the planetary fluid Love number and the dissipation to their reference values (see Table 1), the potential Love number k2,p is varied between 1 and 50$\times k_{2,\text{p}}{}^{\text{ref}}~15$. Bottom panel: fixing the planetary Love number and the dissipation to their reference values, the potential Love number k2f,p is varied between 1 and 50$\times k_{2f,\text{p}}{}^{\text{ref}}~48$.

In the text

	Fig. B.1 As the bottom panels of Fig. 1, but for (a) TRAPPIST-1e, (b) TRAPPIST-1f, (c) TRAPPIST-1g, and (d) TRAPPIST-1h. The general relativity is the dominant effect, but the planetary deformation (due to tides or rotation) is not quite completely negligible, except for TRAPPIST-1h.
In the text