© ESO 2022

1 Introduction

The TRAPPIST-1 planetary system consists of at least seven Earth-size planets around an ultra-cool red dwarf star of mass M_*≈ 0.09 M_⊙ (Gillon et al. 2016, 2017; Luger et al. 2017; Delrez et al. 2018; Van Grootel et al. 2018; Ducrot et al. 2018; Agol et al. 2021). One remarkable feature of TRAPPIST-1 is that the planets are possibly in the longest resonant chain known to date, with period ratios between adjacent planets close to the following ratios (starting with the innermost pair): 8:5, 5:3, 3:2, 3:2, 4:3, and 3:2 (Luger et al. 2017).

Resonant configurations have been observed in other planetary systems, although they are not the most common orbital configuration (Fabrycky et al. 2014). The formation of mean motion resonances (hereafter MMRs) is often attributed to convergent migration of two planets in a gaseous disc (Snellgrove et al. 2001; Lee & Peale 2002). Furthermore, several exoplanetary systems have been observed in chains of resonances involving three or more planets (see, e.g., Rivera et al. 2010; Goździewski et al. 2016; MacDonald et al. 2016; Mills et al. 2016; Christiansen et al. 2018; Leleu et al. 2021, and the aforementioned TRAPPIST-1 system). Three-body resonances are often referred to as Laplace resonances, in reference to the orbital configuration of the Jovian satellites Io, Europa, and Ganymede. Hydrodynamical simulations have long predicted that resonant chains of low-mass planets could form via type-I migration in gaseous discs (McNeil et al. 2005; Cresswell & Nelson 2006; Terquem & Papaloizou 2007; Ogihara & Ida 2009; Ida & Lin 2010).

Using formation models based either on pebble accretion or planetesimal accretion, several studies have shown that resonant chains such as the one seen in TRAPPIST-1 can naturally arise during the early stages of formation and migration in a disc (Ormel et al. 2017; Schoonenberg et al. 2019; Coleman et al. 2019; Huang & Ormel 2022). As noted by Coleman et al. (2019), the majority of pairs end up in first-order resonances (mainly the 2:1, 3:2 and 4:3 MMRs), with a few in second-order resonance (such as the 5:3 MMR which is relevant for this work), but no pair ends up in third-order resonance (which may be relevant for TRAPPIST-1 since the inner pair has a period ratio of ~1.6, and therefore lies near the 8:5 MMR). Using a simpler approach where planets are assumed to have already formed, and are initially placed just outside their current positions, and where migration forces are only applied to the outermost planet, Tamayo et al. (2017) showed that disc migration allows the capture of the TRAPPIST-1 planets into their current configuration. Such a configuration is stable on very long timescales, as opposed to the early N-body simulations by Gillon et al. (2017) which typically showed instability on a timescale of 0.5 Myr. Finally Papaloizou et al. (2018, see also Brasser et al. 2019; Huang & Ormel 2022) studied the migration and tidal evolution of the system, arguing that the planets may have migrated in two distinct groups, re-assembling later on, thanks to outward migration of the inner planets driven by tidal interactions with the central star.

The TRAPPIST-1 system presents strong TTVs (i.e., deviations of the transit times with respect to the strictly periodic Keplerian case). These variations are useful in order to infer the presence of unseen planets (Miralda-Escudé 2002; Schneider 2004; Agol et al. 2005; Holman & Murray 2005; Cochran et al. 2011; Ford et al. 2012; Steffen et al. 2012) and to provide a dynamical measurement of planetary masses and eccentricities (Nesvorný & Morbidelli 2008; Lithwick et al. 2012; Hadden & Lithwick 2017). This is especially the case for systems harboring planets in or near MMR for which the variations are strongly enhanced. Based on transit timing data presented in Ducrot et al. (2020), a refined TTV model of the TRAPPIST-1 system was recentlypublished by Agol et al. (2021), allowing for an improvement of the mass estimates of the seven planets.

In this paper, we use the new analysis from Agol et al. (2021) to better constrain the formation and dynamics of the TRAPPIST-1 system. Whether all the planet pairs are indeed in resonances, and how did these resonances form are still open questions. In Sect. 2, we carry out a frequency analysis of the best-fit solution from Agol et al. (2021) and characterize the resonant nature of the solution. In Sect. 3, we explore whether this resonant configuration can be produced bydisc migration. Our results are discussed in Sect. 4 and summarized in Sect. 5.

2 Dynamics of TRAPPIST-1

In this section, we study in detail the resonant chain of TRAPPIST-1 based on the latest model inference obtained by Agol et al. (2021). The masses and orbital properties of the system are listed in Table 1. By means of frequency analysis, we also provide a deep analysis of the periods shown by the TTVs to explore how they are connected to the dynamical evolution of the system.

2.1 Resonant chain

The seven planets of the recent best-fit solution of Agol et al. (2021) are all influenced by 2- and/or three-body resonances. The last column of Table 1 gives the orbital period ratios between adjacent planets, as well as the closest low-order two-body MMRs. The recent best-fit TTV model of Agol et al. (2021) (see their Fig. 2), based on 4 yr of observation, is reproduced in the left panel of Fig. 1. Since the planets revolve about the ultra-cool red dwarf star with very short orbital periods (from 1.5 to 19 days), a survey of TRAPPIST-1 on 4 yr was sufficient to clearly see the imprint of the different two-body MMRs between adjacent planets in the TTVs. Indeed, the characteristic timescale of the TTVs (i.e., the period of the sinusoidal TTV signal) for two planets near a j + k :j MMR is (see, e.g, Lithwick et al. 2012, for first order resonances): $$P_{\text{TTV}}=\frac{1}{|j/P_1-(j+k)/P_2|},$$(1)

with P₁ and P₂ the orbital periods of the inner and outer planets, respectively. This gives a characteristic period of ~ 1.3 yr for the first-order resonances amongst the five outer planets of TRAPPIST-1, which is clearly visible in the TTVs of planets d to h in the left panel of Fig. 1. For pair b–c, which is near a third-order resonance, the formula gives a period of 0.45 yr. For pair c–d, which is near a second-order resonance, it gives a period of 0.72 yr. However, the TTVs also present several additional shorter and longer periods which we aim to identify precisely in the following based on the analysis of the dynamical evolution of TRAPPIST-1.

Given the data of Table 1, the relevant two-body resonant angles are: $$\begin{array}{ll}{\theta }_{4:3}^{(1)}\hfill & =4{\lambda }_2-3{\lambda }_1-{\varpi }_1,\hfill \\ {\theta }_{4:3}^{(2)}\hfill & =4{\lambda }_2-3{\lambda }_1-{\varpi }_2,\hfill \\ {\theta }_{3:2}^{(1)}\hfill & =3{\lambda }_2-2{\lambda }_1-{\varpi }_1,\hfill \\ {\theta }_{3:2}^{(2)}\hfill & =3{\lambda }_2-2{\lambda }_1-{\varpi }_2,\hfill \\ {\theta }_{8:5}^{(1)}\hfill & =8{\lambda }_2-5{\lambda }_1-3{\varpi }_1,\hfill \\ {\theta }_{8:5}^{(2)}\hfill & =8{\lambda }_2-5{\lambda }_1-3{\varpi }_2,\hfill \\ {\theta }_{8:5}^{(3)}\hfill & =8{\lambda }_2-5{\lambda }_1-2{\varpi }_1-{\varpi }_2,\hfill \\ {\theta }_{8:5}^{(4)}\hfill & =8{\lambda }_2-5{\lambda }_1-{\varpi }_1-2{\varpi }_2,\hfill \\ {\theta }_{5:3}^{(1)}\hfill & =5{\lambda }_2-3{\lambda }_1-2{\varpi }_1,\hfill \\ {\theta }_{5:3}^{(2)}\hfill & =5{\lambda }_2-3{\lambda }_1-2{\varpi }_2,\hfill \\ {\theta }_{5:3}^{(3)}\hfill & =5{\lambda }_2-3{\lambda }_1-{\varpi }_1-{\varpi }_2,\hfill \end{array}$$

for two adjacent planets labeled 1 and 2, with orbital periods P₁ <P₂, and with ϖ the longitude of pericenter and λ the mean longitude. These angles can be combined to form the following three-body resonant angles for the successive triplets of planets: $$\begin{array}{l}{\theta }_{10:25:15}=10{\lambda }_1-25{\lambda }_2+15{\lambda }_3,\\ {\theta }_{3:9:6}=3{\lambda }_1-9{\lambda }_2+6{\lambda }_3,\\ {\theta }_{2:5:3}=2{\lambda }_1-5{\lambda }_2+3{\lambda }_3,\\ {\theta }_{2:6:4}=2{\lambda }_1-6{\lambda }_2+4{\lambda }_3,\\ {\theta }_{3:6:3}=3{\lambda }_1-6{\lambda }_2+3{\lambda }_3.\end{array}$$

In order to study further the resonant nature of TRAPPIST-1, we integrate it forward in time, starting from the initial conditions given by the best fit model of Agol et al. (2021). Throughout this paper, simulations made use of the REBOUND N-body code (Rein & Liu 2012). The simulations were integrated using WHFast, a symplectic Wisdom-Holman integrator (Rein & Tamayo 2015; Wisdom & Holman 1991). Our integrator includes relativistic corrections as implemented in REBOUNDx (Tamayo et al. 2020), but neglects any tidal interactions. In Fig. 2, we show the orbital evolution (the eccentricities and the pericenter arguments) of the 7 planets over the same 1600 days baseline as Agol et al. (2021). The eccentricities clearly exhibit some periodic variations, with a distinctive pattern for planets b and c on one hand, and all the other planets on the other hand. This is more apparent when looking at the arguments of pericenters which have a nearly constant average value for planets b and c, and are circulating for the other planets. There is also a rapid oscillation for planets b and c which occurs on a timescale of ~ 11.76 days which is driven by the 3:2 resonant term; this has a small amplitude given that their period ratio is far from 3/2.

In the left panel of Fig. 3, we computed the two-body resonant angles based on the proximity to the given MMR shown in Table 1. The best-fit solution shows a resonant behavior, with most two-body angles librating, apart from the angles associated with the 8:5 MMR between the innermost two planets which are circulating. We have also checked the behavior of the 3:2 resonant angles of this innermost pair, and did not observe a resonant behavior. However, while not in MMR, the pair b–c presents an alignment of the pericenters, as shown by the libration of the difference of the pericenter arguments around 0, with a period of 40 yr (bottom panel of Fig. 2).

Some similarities in periodicity are apparent between Figs. 1, 2, and 3. In particular, the two-planet resonant angles between all pairs except the innermost one show variations with the characteristic period of about 1.3 yr previously discussed, the same period that is observed in the TTVs of planets d to h.

The three-body resonant angles quoted in Eqs. (13)–(17) are also shown in Fig. 3, on the same timescale of 1600 days as in Agol et al. (2021) in the middle panel and on a longer timescale of 100 yr in the right panel. As previously noted by Libert & Renner (2013), the dynamics of three-planet (or Laplace) resonances occurs on a timescale longer than that of two-planet MMRs. The resonant behavior shown in Fig. 3 suggests that the system is dominated by two-body resonances over the observed baseline of 1600 days, while it is dominated by three-body resonances over longer timescales of over a decade. Moreover, we observe that all the three-body resonant angles of the best-fit solution librate. As expected, they appear nearly constant on the 4 yr of observation, while they present interesting periods on their long-term evolution.

In light of the above, we computed the TTVs of the best-fit solution of TRAPPIST-1 on a timescale of 100 yr in the right-hand side of Fig. 1. We will further explore the periodicity of the “short-term” and “long-term” TTV signals in the next section.

Fig. 1 TTVs of the best-fit solution for TRAPPIST-1 from Agol et al. (2021). The left panel shows the TTVs over the 1600 days ofobservation as in Fig. 2 of Agol et al. (2021), while the right panel shows the TTVs over 100 yr.

Fig. 2 Orbital evolution of the Agol et al. (2021) best fit: eccentricity in the top panel, argument of pericenter in the middle panel, and difference of arguments of pericenters of planets b and c in the bottom panel (note the different timescale onthe bottom panel).

Fig. 3 Time evolution of the resonant angles of the Agol et al. (2021) best fit. The left column shows the two-body resonant angles from Eqs. (2) to (12). The color scheme is as follows: blue is first angle θ(1), orange θ(2), and if present, green and red θ(3) and θ(4), respectively. The middle and right columns show the three-body resonant angles from Eqs. (13) to (17). In the middle column they are plotted over the same 1600 days timescale as the two-body angles. On the right column they are plotted on a longer timeframe of 100 yr.

2.2 Frequency analysis

In order to further investigate the dynamics of the system, we carried out a frequency analysis of the best-fit solution. The frequency analysis of the TTVs and resonant angles was done using the method outlined in Laskar (1993). This method approximates a periodic signal as a finite sum of trigonometric functions, and allows for a very accurate computation of the amplitude, frequency, and phase of each component of the approximate signal.

Since the three-body resonances impact the TTVs on timescales longer than two-body resonances (Libert & Renner 2013), we conducted the frequency analysis on two datasets: the “short-term” dataset of 1600 days, and the “long-term” dataset of 100 yr. By carrying out a frequency analysis on the two datasets, we can extract the relevant periods of the TTV signal of each planet. In total, nine periods were unveiled, most of them being present in the TTV signal of each planet but with different amplitudes. They are listed in the first column of Table 2.

Regarding the short-term dynamics of the TTVs, the signals of the outer five planets are dominated by a periodic modulation of 1.3 yr, as can be seen on the left-hand side column of Fig. 1. An additional very short period of 0.09 yr is also present, and particularly visible in planets e–h. The importance of these two periods will be further discussed in Sect. 4. Planets b and c, on the other hand, are mostly dominated by a short period variation which is three times faster, at 0.032 yr (corresponding to the ~ 11.76 days period seen in the eccentricities and arguments of pericenter of planets b and c), as well as an additional period variation of 0.45 yr.

The long-term TTVs mainly vary on a period of 31.5 yr for all the planets, as visible on the right-hand side column of Fig. 1. Planets d–g possess a common period of 3.3 yr. A common period of 5.1 yr is shared across all the planets. As also visible in the last panel of Fig. 1, a period of 12.3 yr is also reported by the frequency analysis.

Table 2 shows the main periods found in the time series of the resonant angles, for the two-body resonances between adjacent pairs in Cols. 2–7 and for the three-body resonances in Cols. 8–12. A tick indicates that the period is present in the signal with a large amplitude. A circled bold tick indicates the periodic component with the largest amplitude. This table allows one to immediately visualize two groups: two-body angles are associated with short periods and three-body angles to long periods. In particular, the period of 1.3 yr is clearly visible in the eccentricity, TTVs and two-body resonances, pointing towards a chain of two-body resonances as the origin of this signal. On the other hand, periods of 31.5, 12.3, 5.1, and 3.3 yr modulate the signal on long timescales, both in TTVs and three-body resonant angles.

Finally, we note that we explored 100 draws from the posterior distribution of Agol et al. (2021). They all exhibit the same resonant structure, and in all of them we recovered the main periods of 1.3 and 31.5 yr. The angles associated with the 8:5 MMR between planets b and c was also circulating in all of the simulations.

Our main conclusion is that on short timescales, the main frequencies of the TTVs derived from the best-fit solution arise from two-body resonances, while they arise from three-body resonances on longer timescales. Our analysis is relevant for follow-up observational campaigns of TRAPPIST-1. In the coming decade, enough data could be collected to refine the best-fit solution. Our analysis suggest that decade-long observational span could give enough time to see the 3.3 and 5.1 yr frequencies arise. More measurements would also help regarding the peculiar dynamics of the inner pair, since its currently observed period ratio of 1.6 is hard to reproduce through disc migration (see Sect. 3). Hence, in addition to being relevant to future observations, our analysis puts constraints on the physical processes that took place during the formation and migration of TRAPPIST-1. In the next section, we explore whether the resonant chain of TRAPPIST-1 can be formed via smooth disc migration.

3 Formation of TRAPPIST-1 by disc migration

In this section, we study the late-stage formation of TRAPPIST-1, in particular the migration phase. We aim to see whether the resonant chain of TRAPPIST-1 discussed in the previous section can be reproduced by disc migration.

3.1 Disc model

Given the masses of the TRAPPIST-1 planets, they are unlikely to be in the gap-opening (i.e., Type-II) migration regime. The planet-to-star mass ratios in the TRAPPIST-1 system range from ~ 10⁻⁵ to ~ 5 × 10⁻⁵. Around a Sun-like star, this would correspond to a super-Earth to Neptune size planet. In this mass regime, we follow Papaloizou et al. (2018) and apply semi-major axis and eccentricity damping forces corresponding to Type-I migration. We use the forces prescribed by Cresswell & Nelson (2008) which we implemented in REBOUNDx (Tamayo et al. 2020). We assume that the disc mass remains constant throughout the simulations.

Let us consider a planet of mass M_p orbiting a star of mass M_*, at a distance r_p and with a Keplerian frequency Ω_p. The planet is embedded in a disc of surface density Σ and aspect ratio H∕r. We define β = −dlnΣ∕dlnr the slope of the surface density. Following Tanaka et al. (2002), we define the wave timescale as: $${\tau }_{\text{w}}=\frac{M_{*}}{M_{\text{p}}}\frac{M_{*}}{\Sigma (a_{\text{p}})a_{\text{p}}^2}{\left(\frac{H}{r}\right)}^4{\Omega }_{\text{p}}^{-1}.$$(18)

Hydrodynamical simulations by Cresswell & Nelson (2008) indicate that the semi-major axis then decays at a rate given by $${\tau }_a=\frac{{\tau }_{\text{w}}}{2.7+1.1\beta }{\left(\frac{H}{r}\right)}^{-2}P(e_{\text{p}})$$(19)

with $$P(e)=\frac{1+{\left(\frac{e}{2.25(H/r)}\right)}^{1.2}+{\left(\frac{e}{2.84(H/r)}\right)}^6}{1-{\left(\frac{e}{2.02(H/r)}\right)}^4}.$$(20)

Note that if the surface density increases sharply with radius (as would for instance be the case at the inner edge of the disc), one would have β < 0. This raises the possibility that the term 2.7 + 1.1β is negative, hence halting or even reverting the direction of migration of the planet as it reaches the inner edge (see, e.g., Kretke & Lin 2012). Similarly, this formula indicates that the migration can also be outward for e > 2.02H∕r (Cresswell & Nelson 2008). Finally, the eccentricity is damped at a rate: $${\tau }_e=\frac{{\tau }_{\text{w}}}{0.78}F(e_{\text{p}})$$(21)

with $$F(e)=1-0.14{\left(\frac{e}{(H/r)}\right)}^2+0.06{\left(\frac{e}{(H/r)}\right)}^3.$$(22)

3.2 Halting Type-I migration at the inner edge

The Type-I migration timescale defined in Eq. (19) is inversely proportional to the planet mass. All other things being equal, more massive planets migrate faster. Resonant capture of two migrating planets requires a convergent migration, that is, the outer planet needs to migrate faster than the inner one. Planets b and c have very similar masses (see Table 1), and their relative migration speed is slow. Although planet b is slightly more massive, intricate details of the formation process and disc structure allow for the possibility of convergent migration. More interestingly, planet d is significantly less massive than planets b and c. Assuming a smooth disc described by a power-law all the way to the inner edge, Eq. (19) indicates that planet d never migrates fast enough to catch up with planets b and c and capture planet c in the 5:3 MMR. This can be seen on the left-hand side panel of Fig. 4, where we plot the migration time of all seven planets, assuming they migrate in a disc with surface density Σ = Σ₀r^−1∕2, where Σ₀ is calculated assuming the disc mass is M_d = 10⁻⁴ M_⊙ between an inner edge r_in = 0.01 au and outer edge r_out = 5 au. It is clear that planet d migrates much slower than planets b and c. However, planet e, f, and g all migrate faster than the planet immediately interior to it. This raises the possibility that planets d, e, f and g were all captured together in a resonant chain, which would explain why their TTVs share so many common periods, as shown in Sect. 2.2. Finally, planet h being of much lower mass, it would have arrived later.

A possible scenario is therefore that planets b and c migrate together until the inner edge of the disc where they stopped migrating. Meanwhile (or perhaps even at a later stage) further out in the disc, planets d to g migrated, with migration speed increasing for further out planets, causing them to form one long chain (or two small chains). They finally encountered planets b and c whosemigration was already halted at the inner edge, captured them in a resonance, and settled at the inner edge of the disc. Eventually planet h came in and captures planet g in a resonance. Note that Papaloizou et al. (2018) also examined the possibility that TRAPPIST-1 migrated as two subsystems. In their case, it was an outward migration caused by tidal interactions with the central star which allowed the two sub-systems to “re-attach” (see also Huang & Ormel 2022, for a more complete scenario invoking both disc migration and stellar tides).

In order for the scenario described above to work, one needs a mechanism to halt migration at the inner edge. As previously noted, the term 2.7 +1.1β in the migration timescale given by Eq. (19) allows for the migration to be slowed down or even reversed for β < 0, i.e. for positive surface density slopes. A sharp edge at the inner end of the disc can provide such a positive slop. We therefore model the disc surface density profile as a combination of a power-law and an inner edge: $$\Sigma (r)={\Sigma }_0{\left(\frac{r}{r_{\text{in}}}\right)}^{-s}\tanh {\left(\frac{r-0.7r_{\text{in}}}{r_{\text{in}}}\right)}^6.$$(23)

This surface density profile is plotted in Fig. 5 for s = 1∕2. The corresponding migration timescales (zoomed in on the inner 0.1 au) are shown on the right-hand side of Fig. 4. As the planets approach the inner edge, their migration speed becomes slower, and can even be reversed (as indicated by the negative migration times). This provides a way to halt migration at the inner edge (see, e.g., Masset et al. 2006; Ogihara et al. 2018).

Fig. 4 Migration time of the TRAPPIST-1 planets in a simple power-law disc (left), and in a disc presenting an inner edge (right) such as the one prescribed by Eq. (23). The presence of an inner edge slows down migration, and can even reverse it (corresponding to negative migration times).

3.3 Initial conditions

We assume that the planets have fully grown to their presently observed mass at the start of the simulation, although the planets might have continued experiencing mass growth during the later stages of migration (see, e.g., Coleman et al. 2019). We assume the planetary system to be coplanar, and we initiate the planets with very low eccentricities (10⁻⁴). The initial orbital periods of the planets are selected using the following relation $$P_{\text{k},\text{0}}=a^k{b}_k{P}_{k,\text{obs}}$$(24)

where P_k,0 is the initial period of planet k, P_k,obs its currently observed period (see Table 1), a is a random number which determines how wide the planets are initiated away from their current period ratio (Tamayo et al. 2017 used a = 1.02), and b_k is a random number between 1.2 and 1.9. Overall, this law causes the planets to start with periods between 1.2 and 2.5 larger than their observed ones. The longitudes of pericenter, longitudes of ascending nodes, and mean anomalies are all randomly selected between 0 and 2π. Relativistic precession is included, as implemented in REBOUNDx.

Regarding our disc model, we chose the disc mass randomly from a uniform distribution between M_d = 10⁻⁵ M_⊙ and M_d = 10⁻⁴ M_⊙, with a radial extension between an inner edge r_in = 0.01 and an outer edge r_out = 5 au, following the estimate provided by Papaloizou et al. (2018). We also assume that the index of the power-law part of the surface density in Eq. (23) is either s = 1∕2 or s = 1. We pick H∕r from a uniform distribution between 0.025 and 0.035. As noted by Cresswell & Nelson (2006), it is sometimes necessary to multiply the eccentricity damping timescale given by Eq. (21) by a numerical factor Q_e to provide a better match between the linear theory and the simulations. They suggested Q_e = 0.1. In preliminary N-body simulations, we found that in order to obtain eccentricities similar to the ones observed in TRAPPIST-1, low values of Q_e were necessary, suggesting a strong eccentricity damping. We chose Q_e from a uniform distribution between 0.02 and 0.1.

Fig. 5 Surface density profile given by Eq. (23) used in our simulations in order to stop migration at the inner edge.

3.4 Results

In this section, we present the results of the simulations of migrating planets in a disc, aimed at reproducing the architecture of TRAPPIST-1. In Fig. 6, we display the outcomes of an ensemble of 1000 simulations carried out using the initial conditions described in Sect. 3.3. We plot the final period ratios of adjacent pairs of planets, as well as the observed period ratios. Apart from pairs b-c and g-h, all pairs of planets are routinely found in configurations matching their observed ones (dashed lines in Fig. 6). The planets c-d are regularly captured in the 5:3 MMR (although the most frequent capture is in the 3:2 MMR), d–e and e–f both in the 3:2 MMR, and f–g in the 4:3 MMR (although many pairs also end in the 3:2 MMR), as found in the observations. The slow migration of planet h, as detailed in Sect. 3.2, combined with our finite integration time, leads to many simulation outcomes where planet h is not captured in any resonance, although we observe pile-ups at the 2:1 and 3:2 MMRs with planet g.

Regarding the inner planets b and c, most pairs are trapped in the 3:2 or the 5:3 MMR, with additional pairs captured in the 2:1 MMR. Seven simulations showed pairs b–c having a period ratio of 1.6. Interestingly, in these seven simulations, the pair c–d was also captured in the observed 5:3 MMR, suggesting that three-body resonances play a significant role during the migration (see also Charalambous et al. 2018). In addition, all the simulations that led to planets b and c having a period ratio of 1.6 had a similar disc mass of M_d ≈ 3.5 × 10⁻⁵ M_⊙.

Fig. 6 Final period ratio found in the 792 simulations of disc migration for the TRAPPIST-1 planets, achieved using the initial conditions described in Sect. 3.3. The dashed vertical line in each panel represents the currently observed period ratio of the given pair (see Table 1).

Fig. 7 Final orbital configuration of one of the simulations described in Sect. 3.3. From top to bottom: semi-major axes, period ratios, eccentricities, arguments of pericenters, and difference of arguments of pericenters between planets b and c.

3.5 Analysis of one simulation

We now focus on one particular simulation where the period ratio between planets b and c was locked around 1.6. For this simulation, the disc parameters are H∕r = 0.0344, s = 0.5, Σ₀ =2.386 × 10⁻⁶M_⊙ au⁻², and Q_e =0.0134. In this particular simulation, planet h was not captured in the 3:2 MMR with planet g, so we restarted a series of simulations with the same initial conditions but a slightly different initial semi-major axis for planet h. We present in Fig. 7 one such simulation where planet h was captured in the 3:2 MMR with planet g. As can be seen on Fig. 7, the final semi-major axes and period ratios are very close to the observed ones. The eccentricities are also in the range of observed values, albeit with smaller amplitude compared to the best-fit solution presented in Fig. 2. We also point out that the angular momentum deficit (AMD, see Laskar 1997) of this simulated system differs only by 5% compared to the AMD of the observed TRAPPIST-1 system¹. One major difference with the best-fit solution shown in Fig. 2 arises when considering the arguments of pericenter. In the best-fit case, planets b and c precess at much lower rates than the rest of the system, and their difference in arguments of pericenters oscillates around 0 with a period of about 40 yr. Regarding the system arising from disc migration, the arguments of pericenter are now all circulating at the same rate, but the difference of arguments of pericenters of the pair b-c is librating around π.

We now examine the resonant state of the simulated system, shown in Fig. 8. Regarding planets d–h, all the two-body and three-body resonant angles show very similar trends, compared to the best-fit solution shown in Fig. 3. Although the amplitudes of libration slightly differ, the main frequencies appearing in the two-body MMRs in Fig. 2 are also visible in Fig. 8. The evolutions of the three-body angles are similar on a short timescale (note that the center of libration of the angles associated with triplets d-e-f and e-f-g are of the reversed sign), while on longer timescales they reveal a main period of 40 yr, slightly longer than the 31.5 yr period of the best-fit solution.

Differences arise when considering the two-body angles associated with the 8:5 MMR between planets b and c, which are in a state of clear circulation in the best-fit case. In the disc-migration simulation, although they spend a large fraction of their time around 0 or π, they are not strictly librating, but rather circulating. Similarly, in the system formed by disc migration, planets c and d do not appear to be as deep in the 5:3 MMR as in the best-fit case. In both cases, the angle ${\theta }_{5:3}^{(2)}=5{\lambda }_{\text{d}}-3{\lambda }_{\text{c}}-2{\varpi }_{\text{d}}$ is librating. However the two other angles behave differently in both sets of simulations due to the different precession rate of ϖ_c. In Fig. 9, we compare the evolution of the angles associated with the 3:2 MMR for pairs b-c and c-d in both sets of simulations, as well as the associated Laplace angle. While the 3:2 angles are not librating in the best-fit case, they are librating with large amplitude in the disc-migration simulation. This libration is only made possible because the arguments of pericenters of planets b and c are precessing much faster in the disc-migration case.

Finally, the TTVs associated with this simulation are shown in Fig. 10. The TTVs of the system formed by smooth migration in a disc reproduce very well the ones of the observed TRAPPIST-1 system. Comparing Figs. 10 and 1, we recover the very fast variations in the TTVs of planets b and c, as well as the 0.09 and 1.3 yr signals in the outer planets. In addition, the amplitudes of the TTVs are similar. As was previously highlighted in Teyssandier & Libert (2020) for the K2-24 planetary system, we have shown that it is possible to recover most of the orbital architecture and TTV signal of TRAPPIST-1 from smooth disc migration.

There is a good agreement between the TTV signals of the best-fit solution and the disc-migration simulation, despite the fact that the resonant dynamics of the pairs b–c and c–d differ between the two cases. In order to understand this discrepancy, we examine the trajectories of the planets in the e cos ω–esinω plane (see Fig. 11). It is clear from this figure that the architecture of the outer system (planets d to h) is very similar between the best-fit solution and the migration simulation. Although the centers are slightly offset, in both cases the trajectories follow annuli of similar radius and thickness. The forced eccentricities of the outer planets are driven by their captures into first-order resonances. Their free eccentricities have been damped during disc migration and are therefore expected to be small. As was previously noted, the main difference is in the trajectories of planets b and c. In the migration simulation, they follow small thick annuli centered on (0, 0), meaning that their free eccentricities are very small compared to their forced eccentricities. In the best-fit case, the precession of planets b and c is very slow, and both planets have a larger value of free eccentricity than forced eccentricity. Thus, over the 4 yr window represented in Fig. 11, the trajectories of these two planets are off-centered from the origin. This raises the question of why the TTVs are similar in both cases. As noted by Lithwick et al. (2012) and Linial et al. (2018), TTVs are mostly influenced by the forced eccentricities and a linear combination of free eccentricities, which are similar in both sets of simulations. In the best-fit case, the free eccentricities of planets b and c are offset by a comparable amount, so that the trajectories of their eccentricity vectors maintain similar behaviors on long timescales. Thus, it is difficult to constrain the absolute values of the free eccentricities of each planet, while the relative values are better constrained. Hence the planets still experience the same proximity at conjunctions, and thus feel the same kicks, resulting in the same TTV pattern.

Fig. 8 Resonant state of the simulation shown in Fig. 7. The columns are the same as in Fig. 3.

Fig. 9 Time evolution of the resonant angles related to the 3:2 MMR between planets b and c (top panel), c and d (middle panel), and their combination as a Laplace angle (bottom panel), for the best-fit solution (left column) and the illustrative simulation of Fig. 7 (right column).

Fig. 10 TTVs computed at the end of the simulation shown in Fig. 7.

Fig. 11 Dynamical evolution of the TRAPPIST-1 system represented in the e cos ω–esin ω plane. Left: evolution of the best-fit solution from Agol et al. (2021). Right: evolution of the system obtained through disc migration. In both panels, the system is represented over the course of 4 yr.

4 Discussion

The numerical experiments presented in Sect. 3.4 show that under generic initial conditions, forming a system like TRAPPIST-1 through disc migration alone is a rare event. Since only a handful of our simulations gave the correct period ratios for pair b-c, it is impossible to estimate the occurrence rate of systems where the angles associated with its 8:5 MMR librate or circulate. We therefore ran 200 additional simulations similar to the ones presented in Tamayo et al. (2017): all planets pairs start within 2% of their observed period ratios, and only the outermost planet migrates with a fixed timescale of 3 × 10⁴ yr (eccentricitydamping is applied to all planets). For planets b-c, we find that 92% of the simulations show a libration of the angles associated with the 3:2 MMR (even though their period ratio is ~ 1.6). The rest show libration of the 8:5 resonant angles. Among those, some even show a very slow variation of the arguments of pericenter of planets b and c (as is observed in the best-fit solution). Regarding planets c–d, 46% show libration of the 3:2 MMR angles, while 54% show libration of the 5:3 MMR angles.

In order to efficiently estimate how close these systems are to the observed one, let us first consider the orbital periods given in Table 1. The orbital period ratio for each pair of planets suggests that they are close to the ratio between two integers from the following series of integers: n = [2, 3, 4, 6, 9, 15, 24]; that is, the orbital frequencies correlate tightly with this integer sequence (the first integer n = 2 is associated with planet h, while n = 24 is associated with planet b). In Fig. 12, we plot the planets’ orbital frequencies versus this series of integers. As expected, the orbital frequencies nearly lie on the same line, whose slope corresponds to a period of P₀ =36.1 d, which matches the 0.09 yr period found in our frequency analysis (see Table 2). Once this linear variation is subtracted, one can see on the bottom panel that the orbital frequencies are all offset from the period ratios by almost exactly the same frequency, with very small residuals (≲ 10⁻⁴ day⁻¹). This frequency of ~ 0.002 day⁻¹ corresponds to a period of P_TTV = 492 days, which matches the main TTV period of 1.3 yr (as seen in Table 2). This suggests that in a frame rotating retrograde to the planets’ orbits with a period 492 days, the orbital periods are all close to harmonics of a same period of 36.1 d. In comparison, we show in Fig. 13 the same analysis for our best-case migration simulation. Using the same sequence of integers, the slope corresponds to a similar period, P₀ =36.2 days, and an offset period P_TTV = 504 days (with residuals ≲1.5 × 10⁻⁴ d⁻¹), which is also similar to the observed frequency offset.

The two parameters, P₀ and P_TTV, allow us to estimate whether or not simulations are close to the observed system. For the disc migration simulations presented in Sect. 3.4, as well as the simplified simulations presented in this section, we fitted the same series of integers to the final periods, and deduced the values of P₀ and P_TTV. We discarded systems with a scatter in the residuals which is larger than 2 × 10⁻⁴ day⁻¹. The results are shown in Fig. 14. The simulations presented in this section (where only planet h is migrating) all have a similar P₀ since they started very close to their observed period ratios, but present a large scatter in P_TTV, and none of them match the observed P_TTV. Among the simulations presented in Sect. 3.4, the particular simulation that we illustrated in Sect. 3.5 (blue circle in Fig. 14) is the closest to the observed system (orange cross). This method can quickly classify systems depending on whether their period ratios are close to a resonant chain, and identify those simulated systems which resemble the observed system.

Fig. 12 Top panel: The blue points (planets h to b, left to right) are the observed orbital frequencies P−1 (where the P’s are the planets’ orbital periods, see Table 1) as a function of the series of integers n = [2, 3, 4, 6, 9, 15, 24]. When applying a linear fit (green line), the slope corresponds to a period P0 = 0.09 yr. Bottom panel: plot of $P^{-1}-n{P}_0^{-1}$ as a functionof n, showing a residual frequency corresponding to a period of ~1.3 yr.

Fig. 13 Same as Fig. 12, but with the orbital frequencies obtained via disc migration (see Sect. 3.5).

Fig. 14 Relevant periods of simulated systems: P0 and PTTV (see Sect. 4). Green points are for the simulations presented in Sect. 3.4 and grey points for the simulations presented in Sect. 4. The blue circle represents our best-case disc migration simulation (Sect. 3.5), and the orange cross represents the best-fit solution of Agol et al. (2021).

5 Conclusion

In this work, we have studied the dynamics and formation of the best-fit solution to the TRAPPIST-1 system given by Agol et al. (2021). In the first part, we focused on the dynamics acting on timescales accessible by future observations. In particular we focus on a short timescale of 4 yr corresponding to the extent of the existing dataset published by Agol et al. (2021), and on a longer timescale of 100 yr which reveals dynamical effects visible on this time frame. Our frequency analysis has revealed an intricate chain of resonant interactions that shapes the dynamics of the system. We find that the dynamics is dominated by two-body resonances on short timescales, and by three-body resonances on long timescales. This is in agreement with the analysis of Libert & Renner (2013) which shows that three-body resonances impact the TTV signal on timescales much longer than two-body resonances. We show that all the planets are simultaneously in two-body and three-body resonances, except for the innermost pair.

In the second part, we carried out ~1000 N-body simulations of planets migrating in a disc in an attempt to reproduce the resonant chain of the TRAPPIST-1 system. For planets migrating in discs with a steep inner edge and a strong eccentricity damping (Q_e of 0.02–0.1), we found that the outer system was well reproduced by our migration scenario, while the dynamics of the inner system (the eccentricity vectors of planets b and c in particular) did not match that of the best fit solution. Most of the two-body and three-body resonances between the TRAPPIST-1 planets are easily formed, apart from the 8:5 MMR between planets b and c. With regardto these two planets, we found that despite their proximity to the 8:5 MMR, it is the angles associated with the 3:2 MMR which are librating. Papaloizou et al. (2018) argued that this commensurability could be obtained through tidal expansion of the inner pair away from the 3:2 MMR. Since we do not include tidal interactions with the central star, we chose to focus on the few simulations in which disc migration lead planets b and c to stop near their observed period ratio of 1.6. We presented a particular simulation in which thesystem at the end of the disc phase was similar to the observed one. We also computed the TTVs at the end of this simulation and showed that they are similar to the ones given by the best-fit solution of Agol et al. (2021). This exercise shows that it is possible that the resonant chain of TRAPPIST-1 formed through smooth disc migration. However, owing for the difficulty of locking the innermost pair b-c near its observed period ratio, it appears to be a low-probability event with our choice of initial conditions. Whether tidal effects would change the long-term dynamics of this pair is beyond the scope of this paper, but we note that tidal interactions would damp the free eccentricities, further reducing the amplitude of the corresponding TTVs (Lithwick & Wu 2012).

In addition to being relevant to future observations, our dynamical analysis puts constraints on the physical processes that took place during the formation and migration of TRAPPIST-1. The TTVs collected these last four years of observation hold key information about the two-body resonant dynamics of TRAPPIST-1. In the coming years, further monitoring of TRAPPIST-1 will supply extended TTVs in which we predict the signature of three-body resonances will become visible. In particular, enough data could be collected soon to recover the periods of 3.3 and 5.1 yr, as well as to give more precise information on the peculiar dynamics of the inner pair of TRAPPIST-1, which may prove to be insightful for formation scenarios.

Acknowledgements

The work of J.T. is supported by a Fonds de la Recherche Scientifique – FNRS Postdoctoral Research Fellowship. The work of A.S.L. is supported by the Fonds de la Recherche Scientifique – FNRS under Grant No. F.4523.20 (DYNAMITE MIS-project). E.A. acknowledges support from NASA’s NExSS Virtual Planetary Laboratory, funded under NASA Astrobiology Institute Cooperative Agreement Number NNA13AA93A, and the NASA Astrobiology Program grant 80NSSC18K0829. Computational resources have been provided by the Consortium des Équipements de Calcul Intensif (CÉCI), supported by the FNRS-FRFC, the Walloon Region, and the University of Namur (Conventions No. 2.5020.11, GEQ U.G006.15, 1610468 and RW/GEQ2016). The research done in this project made use of the SciPy stack (Jones et al. 2001), including NumPy (Oliphant 2006) and Matplotlib (Hunter 2007), as well as Astropy, (http://www.astropy.org) a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018). Simulations in this paper made use of the REBOUND code which is freely available at http://github.com/hannorein/rebound.

References

All Tables

All Figures

	Fig. 1 TTVs of the best-fit solution for TRAPPIST-1 from Agol et al. (2021). The left panel shows the TTVs over the 1600 days ofobservation as in Fig. 2 of Agol et al. (2021), while the right panel shows the TTVs over 100 yr.
In the text

	Fig. 2 Orbital evolution of the Agol et al. (2021) best fit: eccentricity in the top panel, argument of pericenter in the middle panel, and difference of arguments of pericenters of planets b and c in the bottom panel (note the different timescale onthe bottom panel).
In the text

Fig. 3 Time evolution of the resonant angles of the Agol et al. (2021) best fit. The left column shows the two-body resonant angles from Eqs. (2) to (12). The color scheme is as follows: blue is first angle θ(1), orange θ(2), and if present, green and red θ(3) and θ(4), respectively. The middle and right columns show the three-body resonant angles from Eqs. (13) to (17). In the middle column they are plotted over the same 1600 days timescale as the two-body angles. On the right column they are plotted on a longer timeframe of 100 yr.

In the text

	Fig. 4 Migration time of the TRAPPIST-1 planets in a simple power-law disc (left), and in a disc presenting an inner edge (right) such as the one prescribed by Eq. (23). The presence of an inner edge slows down migration, and can even reverse it (corresponding to negative migration times).
In the text

	Fig. 5 Surface density profile given by Eq. (23) used in our simulations in order to stop migration at the inner edge.
In the text

	Fig. 6 Final period ratio found in the 792 simulations of disc migration for the TRAPPIST-1 planets, achieved using the initial conditions described in Sect. 3.3. The dashed vertical line in each panel represents the currently observed period ratio of the given pair (see Table 1).
In the text

	Fig. 7 Final orbital configuration of one of the simulations described in Sect. 3.3. From top to bottom: semi-major axes, period ratios, eccentricities, arguments of pericenters, and difference of arguments of pericenters between planets b and c.
In the text

	Fig. 8 Resonant state of the simulation shown in Fig. 7. The columns are the same as in Fig. 3.
In the text

	Fig. 9 Time evolution of the resonant angles related to the 3:2 MMR between planets b and c (top panel), c and d (middle panel), and their combination as a Laplace angle (bottom panel), for the best-fit solution (left column) and the illustrative simulation of Fig. 7 (right column).
In the text

	Fig. 10 TTVs computed at the end of the simulation shown in Fig. 7.
In the text

	Fig. 11 Dynamical evolution of the TRAPPIST-1 system represented in the e cos ω–esin ω plane. Left: evolution of the best-fit solution from Agol et al. (2021). Right: evolution of the system obtained through disc migration. In both panels, the system is represented over the course of 4 yr.
In the text

Fig. 12 Top panel: The blue points (planets h to b, left to right) are the observed orbital frequencies P−1 (where the P’s are the planets’ orbital periods, see Table 1) as a function of the series of integers n = [2, 3, 4, 6, 9, 15, 24]. When applying a linear fit (green line), the slope corresponds to a period P0 = 0.09 yr. Bottom panel: plot of $P^{-1}-n{P}_0^{-1}$ as a functionof n, showing a residual frequency corresponding to a period of ~1.3 yr.

In the text

	Fig. 13 Same as Fig. 12, but with the orbital frequencies obtained via disc migration (see Sect. 3.5).
In the text

	Fig. 14 Relevant periods of simulated systems: P0 and PTTV (see Sect. 4). Green points are for the simulations presented in Sect. 3.4 and grey points for the simulations presented in Sect. 4. The blue circle represents our best-case disc migration simulation (Sect. 3.5), and the orange cross represents the best-fit solution of Agol et al. (2021).
In the text