TRAPPIST1 is an 0.09 MSUBSUB star which harbours a system of seven Earthsized planets. Two main features stand out i all planets have similar radii masses and compositions and ii all planets are in resonance. Previous works have outlined a pebbledriven formation scenario where planets of similar composition form sequentially at the HSUB2SUBO snowline 0.1 au for this lowmass star. It was hypothesized that the subsequent formation and migration led to the current resonant configuration. Here we investigate whether the sequential planet formation model is indeed capable to produce the presentday resonant configuration characterized by its twobody and threebody mean motion resonances structure. We carry out Nbody simulations accounting for typeI migration stellar tidal damping disc eccentricitydamping and featuring a migration barrier located at the discs inner edge. Due to the sequential migration planets naturally form a chain of firstorder resonances. But to explain the period ratios of the bcdsystem which are presently in higher order resonances we find that planets b and c must have marched across the migration barrier into the gasfree cavity before the disc has dispersed. We investigate both an early and late cavity infall scenario and find that the early infall model best matches the constraints as well as being more probable. After the dispersal of the gaseous disc stellar tidal torque also contributes towards a modest separation of the inner system. We outline how the insights obtained in this work can be applied to aid the understanding of other compact resonant planet systems.
['celestial mechanics', 'planets and satellites: dynamical evolution and stability', 'planets and satellites: formation', 'planet-disc interactions', 'planet-star interactions', 'protoplanetary discs', 'Astrophysics - Earth and Planetary Astrophysics']
The dynamics of the TRAPPIST1 system in the context of its formation
{'Shuo Huang ( 黄 硕 ), 1 Chris W. Ormel, 1 ★': '1 Department of Astronomy, Tsinghua University, 100084 Beĳing, China  \nAccepted XXX. Received YYY; in original form ZZZ', 'ABSTRACT': "TRAPPIST-1 is an 0.09 𝑀 GLYPH<12> star, which harbours a system of seven Earth-sized planets. Two main features stand out: (i) all planets have similar radii, masses, and compositions; and (ii) all planets are in resonance. Previous works have outlined a pebble-driven formation scenario where planets of similar composition form sequentially at the H2O snowline ( GLYPH<24> 0 GLYPH<149> 1 au for this low-mass star). It was hypothesized that the subsequent formation and migration led to the current resonant configuration. Here, we investigate whether the sequential planet formation model is indeed capable to produce the present-day resonant configuration, characterized by its two-body and three-body mean motion resonances structure. We carry out N-body simulations, accounting for type-I migration, stellar tidal damping, disc eccentricity-damping, and featuring a migration barrier located at the disc's inner edge. Due to the sequential migration, planets naturally form a chain of first-order resonances. But to explain the period ratios of the b/c/d-system, which are presently in higher-order resonances, we find that planets b and c must have marched across the migration barrier, into the gas-free cavity, before the disc has dispersed. We investigate both an early and late cavity infall scenario and find that the early infall model best matches the constraints, as well as being more probable. After the dispersal of the gaseous disc, stellar tidal torque also contributes towards a modest separation of the inner system. We outline how the insights obtained in this work can be applied to aid the understanding of other compact resonant planet systems.  \nKey words: celestial mechanics - planet-disc interactions - planet-star interactions - protoplanetary discs - planets and satellites: formation - planets and satellites: dynamical evolution and stability", '1 INTRODUCTION': "TRAPPIST-1 is the first M8-dwarf star found to harbour seven transiting planets (Gillon et al. 2016, 2017). While the orbits are located within 0.1 au from the star in a compact (resonant) configuration, TRAPPIST-1's low luminosity renders planets d, e, and f to reside at equilibrium temperatures similar to Earth (Luger et al. 2017). Therefore, the system is at the forefronts of characterization efforts and habitability studies (Lincowski et al. 2018) and a prime target of the James Webb Space Telescope (Gillon et al. 2020; Meadows et al. 2021; Lim et al. 2021; Lin et al. 2021a). A particular intriguing question is whether planet systems like TRAPPIST-1 are common. Lying within the 1,000 closest stars from the Sun, it is natural to conjecture that planet systems like TRAPPIST-1 must be common in the galaxy. However, Sestovic & Demory (2020) and Sagear et al. (2020) argue that the seven planet TRAPPIST-1 discovery has been serendipitous, based on a study of over 500 low mass targets observed by the K2 spacecraft. Still, due to the challenging nature of detecting small planets around low-mass stars, the occurrence of TRAPPIST-1-like multi-planet systems is poorly constrained.  \nTo characterize the planets of the TRAPPIST-1 system, Delrez et al. (2018) and Grimm et al. (2018) evaluate the planetary radii and masses via transit depths and Transit Timing Variation (TTV) from Spitzer data. Recently, Ducrot et al. (2020) present all transits records for TRAPPIST-1 from Spitzer. Agol et al. (2021) combine  \nSpitzer observations (Delrez et al. 2018; Ducrot et al. 2020), groundbased observations (Gillon et al. 2016, 2017; Ducrot et al. 2018; Burdanov et al. 2019) and other space (K2 and HST) observations (Luger et al. 2017; Grimm et al. 2018), and refine the physical and dynamical properties of the TRAPPIST-1 planets. From these studies it follows that all seven planets in the TRAPPIST-1 system have roughly the same densities (Dorn et al. 2018; Agol et al. 2021), and similar low eccentricities ( 𝑒 glyph[lessorsimilar] 0 GLYPH<149> 01) and inclinations ( 𝑖 glyph[lessorsimilar] 0 GLYPH<149> 01) (Agol et al. 2021). The resonant configuration is one of the most astonishing features in the TRAPPIST-1 system. From inner-to-outer, the period ratio of every adjacent planet pair is near 8:5, 5:3, 3:2, 3:2, 4:3, 3:2. The observed near-integer period ratios likely reflect corresponding Mean Motion Resonances (MMRs). In addition, every adjacent planet triplet is connected by various three-body resonances (3BRs) (Luger et al. 2017).  \nThe similarity among the TRAPPIST-1 planets may reflect a similar formation origin. Because of its compact nature, it is natural that the TRAPPIST-1 planets (or their building blocks) experienced significant radial migration. Planetesimals could originate within the snowline driven by the traffic jam effect of sublimated icy particles (Drążkowska et al. 2016; Ida & Guillot 2016; Schoonenberg et al. 2018). On the other hand, re-condensation from outward diffusion of water vapour enhances the dust-to-gas ratio exterior the snowline, thus triggering streaming instability to form planetesimals (Schoonenberg & Ormel 2017; Ormel et al. 2017; Hyodo et al. 2019). In both scenarios, planets are formed at the same location, which explains the common composition. Population synthesis models where planets are initialized at random locations - also claim to  \nFigure 1. Schematic of the different stages in the dynamical evolution of the TRAPPIST-1 planet system. Time is oriented from left to right with the mass of the disc indicated by the vertical width of the line. Different boxes indicate different events (stages) that occur over specific time intervals. In the Early Cavity Infall model, Stage I and II are contemporaneous.  \n<!-- image -->  \nreproduce the TRAPPIST-1-like properties (Burn et al. 2021). After reaching a certain mass, planets migrate inwards, explaining the similar size of the planets, independent of whether growth is driven by planetesimal accretion or pebble accretion (Schoonenberg et al. 2019; Coleman et al. 2019; Burn et al. 2021; Zawadzki et al. 2021). An intriguing question is how much water these planets were born with and how much still remains. According to Agol et al. (2021) the present-day water mass fraction are much lower than prediction by Ormel et al. (2017) and Schoonenberg et al. (2019). However, 26 Al isotropic heating can desiccate terrestrial planets (Lichtenberg et al. 2019). In addition, Raymond et al. (2021) recently constrains that any water that (had been) present on the TRAPPIST-1 planets is likely to be accreted in the protoplanet disc phase.  \nMost of the above formation models did not address the dynamical evolution of the formed planets. In the TRAPPIST-1 system, 3BRs bind planets together to render the planets long-term dynamically stable (Tamayo et al. 2017). At the same time, Mah (2018) and Brasser et al. (2019) find that the libration centre of the 3BR angles may jump to different values after tens of thousand years. Until now, only several planetary systems are known to harbour 3BRs: GJ 876 (Rivera et al. 2010), Kepler-60 (Goździewski et al. 2016), Kepler-80 (Xie 2013; MacDonald et al. 2016, 2021), Kepler-223 (Mills et al. 2016), K2-138 (Christiansen et al. 2018; Lopez et al. 2019), HD-158259 (Hara et al. 2020), and TOI-178 (Leleu et al. 2021). Recently, one 3BR is confirmed in Kepler-221 (Goldberg & Batygin 2021) seven years after its discovery (Rowe et al. 2014). As hypothesized in this work, these resonances are fossils of the formation and evolution of planetary systems, which we can use to our advantage to reconstruct the planet formation history.  \nThe convergent migration of planets in quasi-circular orbits generally results in first-order MMRs (Lee & Peale 2002; Correia et al. 2018). However, in the TRAPPIST-1 system, the innermost two planet pairs are near second or even third order MMRs. Conceivably, the innermost two planet pairs were near 3:2 MMR just after planet formation because it is the first order MMR closest to the current period ratio, whereafter the inner three planets separated from each other. Lin et al. (2021b) presented an analytical model in which Earth-mass planets in discs around low mass stars end up in first order MMRs chain. However, they did not investigate the deviation from first order MMRs for the inner three planets. Coleman et al. (2019), using a population synthesis approach where planets form at arbitrary locations, observe several cases that have separated inner planets after planet formation and disc dispersal. However, the corresponding 3BRs are inconsistent with the observation. It is likely that a more specific, subtle mechanism operated to shape the dynamical structure of the inner planets.  \nThere are several ways to separate planets away from first-order MMR. In two-planets systems, divergent migration in the disc can drive planets away from MMRs, e.g., by magnetospheric rebound (Liu et al. 2017) or by adjusting the zero net torque location using a migration map approach (Wang et al. 2021). In the absence of a gaseous disc, stellar tidal dissipation can reduce the eccentricities of planets, and divergently migrate planets away from exact MMR (Lithwick & Wu 2012; Batygin & Morbidelli 2013). In multi-planet systems connected by 3BR, divergent migration, as well as tidal dissipation, can separate planets with each other maintaining 3BR relationships (Charalambous et al. 2018; Papaloizou et al. 2018). In particular, Papaloizou et al. (2018) obtain the observed configuration of TRAPPIST-1 planets purely using stellar tidal dissipation. However, their initial conditions arbitrarily lie close to the present state, raising the question of how their initial state did materialize. Teyssandier et al. (2021) find that the observed high order resonances can be reproduced only with a smooth inner disc edge, but even then it is a rare simulation outcome.  \nThe goal of this paper is to connect the present long-term dynamical behaviour of the TRAPPIST-1 planets to their formation in the planet-forming disc. Our model consists of two steps. In the first step, a chain of first-order MMRs is set up in the proto-planet disc. Planets sequentially form and migrate inward one after the other. Migration ceases at the disc's inner edge. The second step involves the separation of the inner b/c/d subsystem until the observed period ratios are reached. For this stage, we distinghuish between late infall, where planets b and c fall in the cavity after the first step has completed, and early infall, in which planets b and c cross the disc edge upon their arrival. After the disc disappears, long-term tidal dissipation further expands the planet system modestly, helping to finally achieve the observed resonant configuration.  \nUsing N-body simulations, we first examine several simple models that fail to establish a more separated inner subsystem. We conclude that planets b and c had to migrate interior to the disc truncation radius, which allowed them to separate already during the disc phase. Finally, accounting for long-term stellar tidal damping, we manage to reproduce the dynamical configuration of TRAPPIST-1 planets including their period ratios, eccentricities, and 3BRs.  \nThis paper is structured as follows. In Sect. 2, we describe the terminologies and our physical method. We then present the results from simple models and assess the conditions for meeting Objective I in Sect. 3. Turning to Objective II, we carry out a parameter study to search for the conditions that can successfully reproduce our objectives (in Sect. 4). We then examine the Early Infall model in Sect. 5. Analytical supports for our findings is given in Sect. 6. Finally, we  \nTable 1. Radii, masses, semi-major axes and eccentricities of the TRAPPIST-1 planets. The values listed in the first three rows are used in the simulations. We compare the eccentricities in our simulation results to the eccentricities listed here from Agol et al. (2021). We also list the observed value of the 3BR angles in the TRAPPIST-1 system (Luger et al. 2017).  \nTable 2. Definition of quantities, symbols and units. The last five parameters separated by horizontal line are the free parameters of the simulations.  \nassess our models and the significance of our findings in Sect. 7 and conclude in Sect. 8.", '2.1 Preliminaries and Terminology': '2.1.1 Simulation stages  \nWe conduct N-body simulations to integrate the planets equation of motion. In this work we follow the dynamical evolution of planets in three stages as shown in Fig. 1:  \n- (I) the stage after their formation, where planets migrates inwards in the gas disc and end up in a chain of first-order resonances;\n- (II) the stage where planets b and c enter the gas-free cavity;\n- (III) the post-disc stage where stellar tides operate on GLYPH<24> Gyr timescales.  \nA key simplification is to simulate each of these phases separately, instead of running one simulation that contains all these stages. Not every model contains each of these three stages. Usually Stage I, II and III occur in chronological order, but in the Early Infall model in Sect. 5, Stage I and Stage II happen simultanously.', '2.1.2 Parameters': 'All masses, semi-major axes and radii of the TRAPPIST-1 planets are adopted from Agol et al. (2021) and are listed in Table 1. The  \nvalues of the default parameters and the range of free parameters appearing in the equations of this paper are listed in Table 2.', '2.1.3 Resonance angles': "For two planets 𝑖 and 𝑖 , 1, a two-body 𝑗 GLYPH<0> 𝑜 𝑖 : 𝑗 MMRangle can be defined  \n𝜙 𝑖GLYPH<150> 𝑖 , 1; 𝑖 , 𝑋 = ' 𝑗 GLYPH<0> 𝑜 𝑖 ' 𝜆 𝑖 GLYPH<0> 𝑗𝜆 𝑖 , 1 , 𝑜 𝑖 𝜛 𝑖 , 𝑋 (1)  \nwhere 𝑜 𝑖 is the order of the resonance, 𝜆 𝑖 and 𝜆 𝑖 , 1 are the mean longitudes of planets 𝑖 and 𝑖 , 1, and 𝜛 𝑖 , 𝑋 is either the longitude of periapsis of planet 𝑖 ( 𝑋 = 0) or 𝑖 , 1 ( 𝑋 = 1). A similar MMR angle can be defined for the ' 𝑘 GLYPH<0> 𝑜 𝑖 , 1 ' : 𝑘 MMR of planets 𝑖 , 1 and 𝑖 , 2. Then we can define the 3BR angle such that it does not involve the longitudes of pericentre:  \nΦ 𝑝GLYPH<150> 𝑝 , 𝑞GLYPH<150> 𝑞 𝑖GLYPH<150> 𝑖 , 1 GLYPH<150> 𝑖 , 2 = 𝑝𝜆 𝑖 GLYPH<0> ' 𝑝 , 𝑞 ' 𝜆 𝑖 , 1 , 𝑞𝜆 𝑖 , 2 (2)  \nwhere 𝑝 = 𝑜 𝑖 , 1 ' 𝑗 GLYPH<0> 𝑜 𝑖 ' and 𝑞 = 𝑜 𝑖 𝑘 . Observationally, it is hard to constrain the longitudes of pericentre when, as is typical, the eccentricity is low. The three-body angles can however be constrained more easily. TRAPPIST-1's c/d/e planets are currently near 3:5 and 2:3 MMRs, corresponding to 𝑝 = 3 and 𝑞 = 6, and the 3BR angle Φ 3 GLYPH<150> 9 GLYPH<150> 6 c GLYPH<150> d GLYPH<150> e = 3 Φ 1 GLYPH<150> 3 GLYPH<150> 2 c GLYPH<150> d GLYPH<150> e . We list the five main 3BR angles of the TRAPPIST1 system in Table 1.", '2.2 Numerical integration': "Every planet obeys the equation of motion which includes the two body forces from the star and other planets, the angular momentum loss and orbital circularization induced by the gas disc and damping from tides with the star:  \n𝑭 i = 𝑁 ∑︁ 𝑗 ≠ 𝑖 𝐺𝑚 𝑗 GLYPH<12> GLYPH<12> 𝑟 𝑗𝑖 GLYPH<12> GLYPH<12> 3 𝒓 ji , 𝑭 disc GLYPH<150> 𝑖 , 𝑭 star GLYPH<150> 𝑖 GLYPH<149> (3)  \nIn Eq. (3) the first term on the right-hand side is the standard gravitational interaction N-body force, where 𝒓 𝑗𝑖 = 𝒓 𝑗 GLYPH<0> 𝒓 𝑖 and the index 𝑗 represent the planets other than planet 𝑖 . The term 𝑭 disc GLYPH<150> 𝑖 indicates the gas disc tidal force and 𝑭 star GLYPH<150> 𝑖 indicates the stellar tidal damping force acting on planet 𝑖 .  \nSimulations were conducted using the WHFast integration method in REBOUND (Rein & Liu 2012) and our additional forces were implemented via REBOUNDx (Tamayo et al. 2020). In Eq. (3), only the 𝑭 disc term is proportional to the disc gas density. During the disc dispersing, we add an exponential decay factor exp »GLYPH<0> ' 𝑡 GLYPH<0> 𝑡 0 'GLYPH<157> 𝜏 d … to 𝑭 disc . In addition, the disc dispersal time-scale 𝜏 d is taken to be 10 5 years. We ignore the disc torque term in the post-disc environment of stage III. The time-step during the simulations is less than 5% of the orbital period of the innermost planet.", '2.3 Planet-disc interaction': 'The planet-disc interaction can cause planet migration (Goldreich & Tremaine 1980; Ward 1997) and eccentricity-damping (Ward 1988; Artymowicz 1993; Tanaka & Ward 2004). Following Papaloizou & Larwood (2000) we write the disc tidal force as:  \n𝑭 disc GLYPH<150> 𝑖 = GLYPH<0> 𝒗 i 2 𝜏 a GLYPH<150> 𝑖 GLYPH<0> 2 \' 𝒗 i GLYPH<1> 𝒓 i \' 𝒓 i j 𝑟 𝑖 j 2 𝜏 d e GLYPH<150> 𝑖 (4)  \nwhere 𝜏 a GLYPH<150> 𝑖 and 𝜏 d e GLYPH<150> 𝑖 indicate the semi-major axis and eccentricitydamping time-scales of planet 𝑖 due to migration.  \nWe assume the gas disc is truncated at a radius 𝑟 c by the magnetic field from the star (Pringle & Rees 1972; D\'Angelo & Spruit 2010). We assume the inner disc is viscously relaxed (Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974) within the snowline 𝑟 ice GLYPH<24> 0 GLYPH<149> 1 au (Schoonenberg et al. 2019; Lin et al. 2021b), implying constant / 𝑀 𝑔 = 3 𝜋𝜈 Σ where 𝜈 is the viscosity and / 𝑀 𝑔 is taken to be 10 GLYPH<0> 10 𝑀 GLYPH<12> yr GLYPH<0> 1 (Manara et al. 2015). The aspect ratio is fixed at ℎ = 0 GLYPH<149> 03 following Ormel et al. (2017), whose value is motivated by viscous heating and lamppost heating (Rafikov & De Colle 2006) and the spectral energy distribution fitting (Mulders & Dominik 2012). In an alphadisc model, 𝜈 = 𝛼ℎ 2 Ω 𝑟 2 with constant 𝛼 parameter (Shakura & Sunyaev 1973), the disc surface density is therefore:  \nΣ g = Σ 0 GLYPH<18> 𝑟 𝑟 0 GLYPH<19> GLYPH<0> 0 GLYPH<149> 5 (5)  \nwhere Σ 0 is the disc surface density at distance 𝑟 0 ( GLYPH<159> 𝑟 ice ). Planets\' semi-major axes and eccentricities are damped on time-scales:  \n𝜏 a GLYPH<150> 𝑖 = ℎ 2 GLYPH<157> 𝑓 a 2 𝛾 𝐼 𝑞 𝑖 𝑞 gas Ω K = 𝜏 a GLYPH<150> GLYPH<8> GLYPH<18> 𝑞 𝑖 𝑞 GLYPH<8> GLYPH<19> GLYPH<0> 1 GLYPH<18> 1 𝑓 a GLYPH<19> (6)  \n𝜏 e GLYPH<150> 𝑖 = 𝐶 𝑒 ℎ 2 𝜏 a GLYPH<150> 𝑖 𝑓 a 𝑓 e GLYPH<150> (7)  \nwhere 𝛾 𝐼 is the prefactor in the type I migration torque expression and 𝑞 𝑖 = 𝑚 𝑖 GLYPH<157> 𝑀 ★ (Tanaka et al. 2002; Kley & Nelson 2012). The semi-major axis damping time-scale is half of the angular momentum damping time-scale (Teyssandier & Terquem 2014). 𝜏 a GLYPH<150> GLYPH<8> is the nominal Type I semi-major axis damping time-scale:  \n𝜏 a GLYPH<150> GLYPH<8> = ℎ 2 2 𝛾 𝐼 𝑞 GLYPH<8> 𝑞 gas Ω K (8)  \n= 1 GLYPH<149> 2 GLYPH<2> 10 4 GLYPH<16> 𝛾 𝐼 2 GLYPH<17> GLYPH<0> 1 / 𝑀 𝑔 10 GLYPH<0> 10 𝑀 GLYPH<12> yr GLYPH<0> 1 ! GLYPH<18> 𝛼 10 GLYPH<0> 3 GLYPH<19> yr GLYPH<150> (9)  \nwhere 𝑞 GLYPH<8> stand for the mass ratio of an Earth mass planet over the stellar mass and 𝑞 gas = Σ g 𝑟 2 GLYPH<157> 𝑀 ★ . Therefore 𝜏 𝑎GLYPH<150> GLYPH<8> and 𝐶 𝑒 are the parameters that control the damping of the planets in the N-body simulations.  \nIn Eq. (7), 𝐶 𝑒 expresses the rate of eccentricity-damping mechanism over semi-major axis damping. In the linear theory its value is 1.28 (Tanaka & Ward 2004), which is consistent with the 3D hydrodynamical simulations (Cresswell & Nelson 2008). However, Cresswell & Nelson (2006) suggest that 𝐶 𝑒 GLYPH<24> 0 GLYPH<149> 1 is sometimes needed to fit their hydro-dynamical simulations while generally this parameter encapsulates the uncertainty in (and variety of) disc migration mechanisms, which amount to considerable uncertainty in the type-I migration prefactor 𝛾 𝐼 (Paardekooper et al. 2011; Baruteau et al. 2014; Benítez-Llambay et al. 2015; McNally et al. 2017). When the disc surface density gradient ( GLYPH<0> d ln Σ 𝑔 GLYPH<157> d ln 𝑟 ) and temperature gradient ( GLYPH<0> d ln 𝑇 GLYPH<157> d ln 𝑟 ) become small, Type I migration slows down  \nFigure 2. Sketch of 𝑓 a and 𝑓 e. The 𝑥 -axis denotes the distance measured from 𝑟 c, the 𝑦 -axis the value of 𝑓 a or 𝑓 e. The horizontal black-dashed lines indicate 𝐴 𝑎 and 𝐴 𝑒 . The vertical dashed dotted line represents the location of 𝑟 c.  \n<!-- image -->  \ndistance from the star  \n( 𝛾 𝐼 decreases) and 𝐶 𝑒 goes down. In this study we take an agnostic view and examine 𝐶 𝑒 values in the range of 0.05 to 1.  \nThe parameters 𝑓 a and 𝑓 e determine the behaviour of semi-major axis and eccentricity-damping near the disc inner edge. When the location of one planet 𝑟 GLYPH<29> 𝑟 c, 𝑓 a = 1. If the planet reaches the disc\'s inner edge, only Lindblad torques from the exterior disc remain. Similarly, only the upper horseshoe motion and associated co-orbital torque is present. In other words, the torques have become one-sided. For low-mass planets, the (positive) one-sided co-rotation torque will win out over the (negative) Lindblad torques, halting migration (Lin & Papaloizou 1979, 1993; Paardekooper & Papaloizou 2009; Liu et al. 2017). Therefore, 𝑓 a must change sign (become negative) at the disc\'s inner edge. We adopt an expression for 𝑓 a in Eq. (6) as:  \n𝑓 a = 8 > > > > < > > > > : 0 \' 𝑟 GLYPH<159> 𝑟 c GLYPH<0> Δ \' 𝐴 𝑎 𝑟 c GLYPH<0> Δ GLYPH<0> 𝑟 Δ \' 𝑟 c GLYPH<0> Δ GLYPH<159> 𝑟 GLYPH<159> 𝑟 c \' 𝐴 𝑎 𝑟 GLYPH<0> 𝑟 c GLYPH<0> Δ Δ \' 𝑟 c GLYPH<159> 𝑟 GLYPH<159> 𝑟 c , Δ , 1 GLYPH<157> 𝐴 𝑎 \' 1 \' 𝑟 ¡ 𝑟 c , Δ , 1 GLYPH<157> 𝐴 𝑎 \' (10)  \nwhere Δ and 𝐴 𝑎 determine the width and strength of the trap. If sharply truncated at the disc\'s inner edge, Δ is comparable to the half-width of the horseshoe region so Δ GLYPH<24> 𝐻 c. For simplicity, we take Δ = 2 𝐻 c, where 𝐻 c is the scale height of gas disc at truncation radius. According to Liu et al. (2017), 𝐴 𝑎 follows as the ratio of the (positive) one-side corotation torque (eq. 11 in Liu et al. 2017) over the Type I migration torque (two times the reciprocal of Eq. (6) times angular momentum)  \n𝐴 𝑎 = Γ c GLYPH<150> 1s Γ I = 𝐶 hs 𝛾 𝐼 ℎ GLYPH<18> 𝑞 𝑝 ℎ 3 GLYPH<19> GLYPH<0> 1 GLYPH<157> 2 = 135 𝛾 𝐼 GLYPH<18> 𝑞 𝑝 10 GLYPH<0> 5 GLYPH<19> GLYPH<0> 1 GLYPH<157> 2 GLYPH<18> ℎ 0 GLYPH<149> 03 GLYPH<19> 1 GLYPH<157> 2 (11)  \nwhere 𝐶 hs = 2 GLYPH<149> 46. For a planet-to-star mass ratio of 𝑞 = 10 GLYPH<0> 5 and ℎ = 0 GLYPH<149> 03 this gives 𝐴 𝑎 = 68. Therefore we take 𝐴 𝑎 around this value. Equation (10) is similar in shape to the "deep drop" model Ataiee & Kley (2021) adopt in their N-body simulations.  \nIn our simulations, we find that it is sometimes advantageous to enhance the eccentricity damping of planets at the disc edge ( 𝑓 𝑒 ¡ 1). Eccentricity damping of planets could be more efficient near the disc\'s inner edge for the following reasons:  \n- GLYPH<15> Material may pile up significantly near 𝑟 c if 𝑟 c is larger than the corotation radius of the star (D\'Angelo & Spruit 2010);\n- GLYPH<15> The absence of inner Lindblad resonances (they lie in the cavity), which would excite planets\' eccentricity (Ward 1988);  \nTable 3. Summary of our model designs.  \n- GLYPH<15> hydro-dynamical simulation indicate that eccentrities for planets at the disc inner edge are more strongly damped (see Fig.15 of Ataiee & Kley 2021).  \nTherefore we express 𝑓 e in Eq. (4) with a similar expression as 𝑓 a:  \n𝑓 e = 8 > > > > > < > > > > > : 0 \' 𝑟 GLYPH<159> 𝑟 c GLYPH<0> Δ \' 𝐴 𝑒 𝑟 GLYPH<0> 𝑟 c , Δ Δ \' 𝑟 c GLYPH<0> Δ GLYPH<159> 𝑟 GLYPH<159> 𝑟 c \' \' 𝐴 𝑒 GLYPH<0> 1 \' 𝑟 c , Δ , 1 GLYPH<157> 𝐴 𝑎 GLYPH<0> 𝑟 Δ , 1 GLYPH<157> 𝐴 𝑎 , 1 \' 𝑟 c GLYPH<159> 𝑟 GLYPH<159> 𝑟 c , Δ , 1 GLYPH<157> 𝐴 𝑎 \' 1 \' 𝑟 ¡ 𝑟 c , Δ , 1 GLYPH<157> 𝐴 𝑎 \' (12)  \nwhere 𝐴 𝑒 is the enhancement of the eccentricity-damping near the truncation radius. Not every simulations feature this anomalous eccentricity damping at the disc edge.  \nWeplot 𝑓 a and 𝑓 e in Fig. 2. Near the truncation radius, 𝑓 a decreases sharply and becomes negative, thus pushing planets in this region outward. The outward torque peaks at 𝑟 c. Within 𝑟 c, the outward migration torque decreases to zero because there is no gas in the disc cavity (See Sect. 5.1, however, for the effects of a distant disc torque on planets inside the cavity).', '2.4 Tidal dissipation': "For a compact planet system in which planets are close to their host star, such as in TRAPPIST-1, tidal dissipation can be significant. The stellar tidal damping force term in Eq. (3) is:  \n𝑭 star GLYPH<150> 𝑖 = GLYPH<0> 2 ' 𝒗 i GLYPH<1> 𝒓 i ' 𝒓 i j 𝑟 𝑖 j 2 𝜏 s e GLYPH<150> 𝑖 GLYPH<150> (13)  \nwhere 𝜏 s e GLYPH<150> 𝑖 refers to the eccentricity-damping time-scale induced by tidal interaction with the host star:  \n𝜏 s e GLYPH<150> 𝑖 = 7 GLYPH<149> 6 GLYPH<2> 10 5 𝑄 sim GLYPH<18> 𝑚 𝑖 𝑚 GLYPH<8> GLYPH<19> GLYPH<18> 𝑀 GLYPH<12> 𝑀 ★ GLYPH<19> 1 GLYPH<149> 5 GLYPH<18> 𝑅 GLYPH<8> 𝑅 𝑖 GLYPH<19> 5 GLYPH<16> 𝑎 𝑖 0 GLYPH<149> 05 au GLYPH<17> 6 GLYPH<149> 5 yr GLYPH<150> (14)  \nwhere 𝑄 sim = 3 𝑄 GLYPH<157> ' 2 𝑘 2 ' , 𝑄 is the tidal dissipation function and 𝑘 2 is the Love number (Goldreich & Soter 1966; Papaloizou et al. 2018). For solar system planets in the terrestrial mass range, 𝑄 sim is estimated between 50-2500. In this paper we take 𝑄 sim = 100 in most cases, similar to the terrestrial planets in our solar system. In Stage III, following Papaloizou et al. (2018), we take 𝑄 sim = 0 GLYPH<149> 1 to accelerate the simulation by a factor 1,000.", '3 PRELIMINARY SIMULATIONS': 'In this section, we summarize our model design and the principal outcomes. This work concentrates on the dynamic evolution after the planet formation process. We insert planets sequentially at the beginning of the simulation with fixed masses following Table 1. Ormel et al. (2017), Schoonenberg et al. (2019) and Lin et al. (2021b) hypothesize that the TRAPPIST-1 planets are formed at the snowline ( GLYPH<25> 0 GLYPH<149> 1 au) and quickly accreted material, while Coleman et al. (2019) and Burn et al. (2021) initialize many planetary embryos in the disc to grow them to TRAPPIST-1 planets. In the context of this work planets could also have formed further out; but we simply fix the mass of every planets during the simulation.  \nDue to type I migration (Tanaka et al. 2002), planets migrate towards the host star when they grow massive enough and stop at the disc truncation radius. The migration speed is so fast that planets tend to be trapped in first order MMR while avoiding the higherorder MMR. Therefore, we initialize every new planet just beyond the 2:1 MMR location of the outermost planet by default. Still, by virtue of the short formation time (Lin et al. 2021b) planets may be born near 3:2 and 4:3 period ratios with the inner planet. Therefore, it is possible that planets are born at locations within the 3:2 MMR. We investigate whether this imprint of planet formation is needed to reproduce the orbital configuration.  \nTherefore, our first step is to park all planets in a chain of first order MMR, i.e., 3:2, 3:2, 3:2, 3:2, 4:3, 3:2 (hereafter referred to as Objective I). However, this contrasts the present-day dynamical state where the period ratio of planets c/b and d/c are 8/5 and 5/3, respectively, i.e., close to higher-order MMRs. Our second step is to evolve the first order MMR chain to the 8:5, 5:3, 3:2, 3:2, 4:3, 3:2 period ratios chain (hereafter Objective II). We separate the simulation into three stages as introduced in Fig. 1.  \nThe set of our models is listed in Table 3. Different columns express model characteristics. In bC model, we set the cavity radius at 𝑟 c = 0 GLYPH<149> 013 au, which coincides with the observed location of planet b. In the bCcgI model, we insert planet g interior to the 3:2 MMR location. In the dCcgI model, we set 𝑟 c = 0 GLYPH<149> 023 au, which coincides with the observed location of planet d, and insert planet g interior to 3:2. In the dCcgIcT model, we additionally switch on the one-sided repulsive torque from the gas disc on planet c in the cavity, while in dCcgIcL we use a particular choice for this repulsive force. The first three models will be discussed in this section. Model dCcgIcT and dCcgIcL will be discussed in Sect. 4 and Sect. 5, respectively.  \nFigure 3. One-by-one migration of the TRAPPIST-1 planets after their sequential formation in the proto-planetary disc for a nominal migration time of 𝜏 a GLYPH<150> GLYPH<8> = 5 GLYPH<2> 10 3 yr and 𝐶 𝑒 = 0 GLYPH<149> 1. The eccentricity enhancement at the disc inner edge is 𝐴 𝑒 = 1. Solid lines with different colours indicate the semimajor axes of different planets. The horizontal black dashed line indicates the position of the migration barrier 𝑟 c. Vertical black and red dashed segments between every two planets trajectories indicate the moment planets converge to and depart from 2:1 and 3:2 period ratios.  \n<!-- image -->', '3.1 Planet b at the cavity edge ( bC )': "Objective I - the parking of planets in a chain of first order MMRs close to the observed period ratios - can most straightforwardly be achieved by the planet's sequential migration and pileup near 𝑟 𝑐 . Therefore, we let 𝑟 c coincide with the location of planet b. We vary 𝜏 a GLYPH<150> GLYPH<8> and 𝐶 𝑒 , and run a set of Stage I simulations; 𝐴 𝑒 is taken to be unity. In the simulations, we insert a new planet every GLYPH<24> 10 kyr (long enough for the inner planets to stabilize). For simplicity, We initialize planet b near 𝑟 c, and then the other planets outside the 2:1 period ratio.  \nOne of our typical outcomes is shown in Fig. 3. We adopt 𝜏 a GLYPH<150> GLYPH<8> = 5 GLYPH<2> 10 3 yr and 𝐶 𝑒 = 0 GLYPH<149> 1. Planet c gets trapped by the 2:1 MMR with planet b. Similarly, planet d gets trapped by the 2:1 MMR with planet c. Planet e directly crosses the 2:1 MMR because planets d & e are less massive and further from the star than planet b, c, and d. Therefore, the Type-I migration can overcome the resonant repulsion and it gets trapped by the 3:2 MMR with planet d. At the same time, planet d crosses the 2:1 MMR because of the pushing from resonant interaction with planet e. Afterwards, the later planets are trapped by the 3:2 MMR with their former peers. The truncation radius acts to prevent planets migrating into the cavity. At the end of this simulation, the period ratios of every two adjacent planets are about 3:2, except planets b & c, which is about 2:1. However, Objective I requires planets b & c to park in 3:2 and planets f & g in 4:3.  \nNaively, closer-spaced MMR can be obtained with a faster migration rate (lower 𝜏 𝑎GLYPH<150> GLYPH<8> ). However, if we take 𝜏 a GLYPH<150> GLYPH<8> GLYPH<159> 5 GLYPH<2> 10 3 yr to migrate planets more rapidly, it results in the breaking of the 3:2 MMRofplanet d & e. Therefore, the value of the nominal migration time-scale 𝜏 a GLYPH<150> GLYPH<8> must exceed 5 GLYPH<2> 10 3 yr when the cavity radius is fixed at 𝑟 c = 0 GLYPH<149> 013 au. Consequently, either b/c and f/g end up in too wide MMRord/e ends up in a too close MMR, thus failing to accomplish Objective I.  \nFigure 4. Similar to Fig. 3, but with planet c starting out inside the 2:1 period ratio and planet g initially inside the 3:2 period ratio. We take 𝜏 a GLYPH<150> GLYPH<8> = 5 GLYPH<2> 10 3 yr, 𝐶 𝑒 = 0 GLYPH<149> 1 and 𝐴 𝑒 = 1.  \n<!-- image -->", '3.2 Formation imprints for planet c and g ( bCcgI )': 'Instead of accelerating planet migration, the closer MMR of planet b & c and planet f & g can be obtained by reducing the time/space span between their formation - a condition that we refer to an imprint of planet formation. In the bCcgI model, the parameters in the simulation are the same as in the bC model. We make the additional assumption that planets c and g are born within 2:1 and 3:2 period ratios in Stage I. One of our typical outcomes is shown in Fig. 4. Planet c gets trapped in the 3:2 MMR with planet b from the outset. Because planets f & g are inserted simultaneously and more massive planet migrate faster Eq. (6), planet g catches up with planet f and gets trapped in the 4:3 MMR. At the end of this simulation, the period ratios of every adjacent two planets are in 3:2 except planet f & g (4:3), fulfilling the requirement of Objective I.  \nObjective II is to separate the b/c/d subsystem from the 3:2 and 3:2 period ratios to 8:5 and 5:3 period ratios. As a first attempt, we disperse the disc to see if we are able to evolve the planet system naturally. During disc dispersal, the semi-major axis damping and eccentricity-damping forces decrease along with the profile of the disc exponentially on a time-scale 𝜏 d = 10 5 yr. After 20 𝜏 d we run the post-disc simulation (Stage III) where only the stellar tidal force operates on the planets.  \nFigure 5 shows the post-disc simulation. The eccentricities of the planets decrease with time because of tidal dissipation during the simulation, while the 3BR angles keep librating around values near 𝜋 . The libration centres of the 3BR phase angles of the last adjacent planet triplets lie around 150 GLYPH<14> due to the perturbation from non-adjacent resonances (Siegel & Fabrycky 2021). Due to tidal dissipation, the gravitational potential energy of this system is decreasing while angular momentum is conserved. Therefore, angular momentum transfers outwards, from inner planets to outer planets and the period ratios of every adjacent planet pair increase together (Batygin & Morbidelli 2013; Papaloizou et al. 2018; Goldberg & Batygin 2021).  \nHowever, the simulation does not result in the observed TRAPPIST-1 planetary orbital configuration. The 3BRs that connect the planets are inconsistent with the observation, i.e., the angle Φ 2 GLYPH<150> 5 GLYPH<150> 3 c GLYPH<150> d GLYPH<150> e (see the notation in Sect. 2) does not librate in the observation. Stellar tides increase the period ratios especially those of the outer planet pairs. Additionally, this simulation outcome results in  \nFigure 5. Period ratios (top panel), eccentricities (bottom-left panel) and 3BR angles (bottom-right panel) during tidal dissipation (with 𝑄 sim = 0 GLYPH<149> 1). In the period ratios plot, dashed lines indicate the 3:2 and 4:3 period ratios. Arrows indicate the observed value of the period ratios. In the eccentricity plot, we label the name of each planets on the left and indicate the present-day eccentricities by error bars with different colours (the x-coordinate of the error bars does not have meaning). In the 3BR angles plots, we label the expression of each angle on the upper left of each panel. The angle labelled in red does not librate at present.  \n<!-- image -->  \neccentricities of planet b & c too high compared to the other planets, which is inconsistent with TTV analysis (Agol et al. 2021).  \nTo sum up, bCcgI model succeeds in Objective I but fails in Objective II. Rather than separating the entire planet system, only the separation of the b/c/d subsystem after Stage I is required. Specifically, this requires that Φ 2 GLYPH<150> 5 GLYPH<150> 3 c GLYPH<150> d GLYPH<150> e needs to be broken before the disc disappears because stellar tidal torques alone cannot. Since disc migration alone rarely reproduces the resonant chain in the TRAPPIST1 system (Teyssandier et al. 2021), there must be special events that drive the expansion of the b/c/d subsystem before disc dispersal.', '3.3 Planet d at the cavity edge ( dCcgI )': 'Ormel et al. (2017) discussed one scenario to decouple the innermost two planets with the disc as the disc truncation radius moves outward due to magnetospheric rebound (Liu et al. 2017), allowing the 3:2 MMRof the inner two planet pairs to break. Here, we aim to obtain the same result by migrating planet b & c across the cavity radius during disc dispersal, motivated by that planet b and c are three times more massive than planet d (Table 1), therefore are easier to open a gap in the disc. Then, stellar tides might increase the period ratios of c/b and d/c Papaloizou et al. (2018). Therefore, 𝑟 c is set to be 0 GLYPH<149> 023 au, coinciding with the observed location of planet d.  \nOne of our typical results for Stage I is shown in Fig. 6. Similar to bCcgI model in Fig. 4, all seven planets park themselves at the desired first order MMRs. To obtain the configuration that planets b and c are in the cavity while the outer planets are in the disc, we switch off the disc torque term (in Eq. (3)) on the planet at the disc inner edge, ignoring for simplicity the details of how the disc profile evolves (Stage II). Figure 7 shows the simulation results for Stage II. At the beginning of this simulation, all planets are beyond the  \nFigure 6. Similar to Fig. 4, but with 𝑟 c = 0 GLYPH<149> 023 au. We take 𝜏 a GLYPH<150> GLYPH<8> = 1 GLYPH<149> 4 GLYPH<2> 10 4 yr, 𝐶 𝑒 = 0 GLYPH<149> 1 and 𝐴 𝑒 = 1.  \n<!-- image -->  \ntruncation radius 𝑟 c. After first switching off the disc torque on planet b, the outer planets push planet b further into the cavity until planet c has reached the cavity radius. Because its mass is small, planet h migrates slowest and is left behind the other planets. Therefore, the period ratio of planet h/g increases until the migration of planet g is stopped by the resonant interaction with the other planets in the disc. Then planet h catches up with its inner planets and the 3:2 MMR reforms. Meanwhile, stellar tides separate planet b in the cavity from exact 3:2 resonance and the period ratio of planet c/b increases to GLYPH<25> 1 GLYPH<149> 56. Next, we switch off the disc torque on planet c, after which, the period ratio of planet h/g does not increase so much and the 3-body resonance angle Φ 3 GLYPH<150> 6 GLYPH<150> 3 f GLYPH<150> g GLYPH<150> h = 3 𝜆 f GLYPH<0> 6 𝜆 g , 3 𝜆 h librates.  \nHowever, after planet d stops migrating at the truncation radius,  \nFigure 7. Stage II of model dCcgI . In the semi-major axis plot, arrows refer to the time that we switch off the disc force on the corresponding planet. In the 3BR angle plots, we label the expression of each angle in the upper left of each panel. The label in red means that the 3-body angle does not presently librate.  \n<!-- image -->  \nstellar tides are unable to sufficiently expand the b/c/d subsystem, i.e. to increase the period ratios of planet c/b and d/c to their observed values. As can be seen from Fig. 7, the 3BR angle Φ 2 GLYPH<150> 5 GLYPH<150> 3 c GLYPH<150> d GLYPH<150> e = 2 𝜆 c GLYPH<0> 5 𝜆 d , 3 𝜆 e still librates, preventing planet c from moving further inward. Therefore, we need an additional mechanisms to break the 3BRs, specifically Φ c GLYPH<150> d GLYPH<150> e , such that the planet b/c/d subsystem is allowed to expand until the observed period ratios (Objective II). This breaking of the 3-body resonance must occur before the disc disappears as otherwise all planets would undergo separation by stellar tides, similar to what is going on in Fig. 5. Therefore, a more sophisticated model is required to expand the b/c/d subsystem.', '3.4 Parameter variation for Objective I ( dCcgI )': 'To assess the likelihood of fulfilling Objective I, we conduct a sensitivity study of the model parameters during Stage I of the dCcgI model, assuming planets all start migrating outside the 2:1 MMRs of their inner planets, except for planet g, which starts just inside 3:2. The goal is to investigate whether (and which) additional planets feature such formation imprints, i.e., whether it is necessary to initialize planets inside 3:2 MMR or 2:1 MMR.  \nWe run a set of Stage I simulations, varying three parameters: the migration time-scale 𝜏 a GLYPH<150> GLYPH<8> , the eccentricity-damping prefactor 𝐶 𝑒 and the eccentricity-damping enhancement at the truncation radius 𝐴 𝑒 . The parameter 𝐶 𝑒 is taken in the range of [0.1, 1], 𝜏 a GLYPH<150> GLYPH<8> in the range of [5, 20] kyr and 𝐴 𝑒 within [1, 40]. We always initialize planet g inside the 3:2 MMR of planet f and other planets wide of the 2:1 MMRof their inner planets.  \nThe simulation results are shown in Fig. 8. We divide the simulation outcomes into four groups according to:  \n(i) the final orbital configuration of planet b & c: circles indicate that planet b & c are connected by 3:2 MMR and triangles indicate that planet b & c are connected by 2:1 MMR.  \n(ii) the final configuration of the outer planets other than planet b &c: green symbols indicate that all outer planet pairs are consistent with Objective I and blue symbols that at least one planet pair is connected by a wider MMR than envisioned by Objective I (and no pair in closer MMRs).  \n(iii) red crosses indicate that there is at least one planet pair connected by a closer MMR than stipulated in Objective I.  \nIn panel (a), more green and blue triangles are present at higher 𝜏 𝑎GLYPH<150> GLYPH<8> since planets are easier trapped into wider MMRs - i.e., the 2:1 - under conditions of slower migration. As 𝐶 𝑒 increases, the green and blue triangles move to the left. Higher 𝐶 𝑒 results in higher planet eccentricities, decreasing the critical time-scale to cross the resonance (consistent with the analysis in Sect. 6.1). Other panels show similar patterns but the higher the 𝐴 𝑒 , the easier for planet c to cross the 2:1 MMR with planet b. Therefore, more circles appear at high 𝜏 a GLYPH<150> GLYPH<8> as 𝐴 𝑒 increases. The analytical explanation can be found in Sect. 6.1. All green circles satisfy the requirement of Objective I. For runs labelled by blue markers, it is mostly planet h that fails to cross the 2:1 MMR. Therefore, we conclude that a formation imprint may be necessary for planets c, g and h for Objective I to be fulfilled - its outcome is surely not universal but does also not require extremely tuned conditions.  \nSince one of our assumptions is that the disc is viscously relaxed (as described in Sect. 2), the 𝛼 -viscosity parameter can be constrained using Eq. (9). In Fig. 8, green circles are centred near 𝜏 a GLYPH<150> GLYPH<8> = 12 kyr, corresponding to a value of 𝛼 GLYPH<24> 10 GLYPH<0> 3 . 𝛼 GLYPH<24> 10 GLYPH<0> 3 is a typical value seen in observational studies, e.g., on the disc gas radii (Trapman et al. 2020) and stellar accretion rates (Hartmann et al. 2016), as well as theoretical studies, e.g. on modelling the dust size distribution (Birnstiel et al. 2018). A higher value 10 GLYPH<0> 2 GLYPH<0> 10 GLYPH<0> 3 is suggested in the close-in region in the disc (Gammie 1996; Carr et al. 2004). Note that our modelling constrains the migration time-scale parameter 𝜏 𝑎GLYPH<150> GLYPH<8> ;  \nFigure 8. Sensitivity study of the dCcgI model. In all five panels the 𝑥 -axis denotes the migration time-scale 𝜏 a GLYPH<150> GLYPH<8> and the 𝑦 -axis the eccentricity-damping 𝐶 𝑒 , but in each panel the eccentricity-damping enhancement at the truncation radius parameter, 𝐴 𝑒 , differs. We label the value of 𝐴 𝑒 at the top of each panel. Different symbols represent different final orbital configuration in the simulation. Circles represent that planet b & c are connected by 3:2 MMR and triangles represent that planet b & c are in 2:1 MMR. Red crosses indicate there is at least one planet pair in a closer MMR than the desired configuration of Objective I, green indicates that all outer planet pairs are consistent with Objective I and blue indicates that at least one planet pair is in an MMR wider than the desired resonance chain.  \n<!-- image -->  \nwhenconverting to the disc proprieties (like 𝛼 ) additional uncertainty arises from other model-dependent parameters, i.e., / 𝑀 ★ , 𝛾 𝐼 .', '4 TRAPPING AND ESCAPE OF PLANETS IN 3BRS': "This section examines how the TRAPPIST-1 planets evolve from Objective I to Objective II, in which the inner three planets attain their final period ratios. We postulate that this is achieved after planets b and c enter the cavity and experience a phase of inward migration. Specifically, we envision that planet c is being pushed inwards, from an initial period ratios of 3:2 (with planet d) to the final ratio of 5:3. The rationale behind this scenario is that if the Φ 2 GLYPH<150> 5 GLYPH<150> 3 b GLYPH<150> c GLYPH<150> d 3BR is maintained during this process, the inner planets will naturally end up near their observed period ratios.  \nHowever, while (gently) pushing planet c inwards, it needs to overcome other 3BR with the outer planets. We list four 3BRs that may halt the inward migration of planet c and therefore the expansion of the b/c/d subsystem. In order of increasing period ratio 𝑃 d GLYPH<157> 𝑃 c, as illustrated in Fig. 9:  \n- GLYPH<15> Φ 1 GLYPH<17> Φ 2 GLYPH<150> 5 GLYPH<150> 3 c GLYPH<150> d GLYPH<150> e with period ratio 𝑃 d GLYPH<157> 𝑃 c = 1 GLYPH<149> 5;\n- GLYPH<15> Φ 2 GLYPH<17> Φ 3 GLYPH<150> 15 GLYPH<150> 12 c GLYPH<150> e GLYPH<150> f (non-adjacent) with period ratio 𝑃 d GLYPH<157> 𝑃 c = 1 GLYPH<149> 556;\n- GLYPH<15> Φ 3 GLYPH<17> Φ 5 GLYPH<150> 14 GLYPH<150> 9 c GLYPH<150> d GLYPH<150> e with period ratio 𝑃 d GLYPH<157> 𝑃 c = 1 GLYPH<149> 6;\n- GLYPH<15> Φ 4 GLYPH<17> Φ 3 GLYPH<150> 9 GLYPH<150> 6 c GLYPH<150> d GLYPH<150> e with period ratio 𝑃 d GLYPH<157> 𝑃 c = 1 GLYPH<149> 667.  \nThe last 3BR is the observed period ratio while the first three need to be crossed. We therefore search for parameter combinations that can  \nFigure 9. Four 3BR angles in which planet c could be trapped when it migrates from its 3:2 (top) period ratio to its observed value of 5:3 (bottom). The period ratios in blue are consistent with the present observation while those in red are inconsistent. Several three-body resonances need to be passed. Γ 1s (or its corresponding time-scale 𝜏 1s) represents the torque from the exterior remnant disc.  \n<!-- image -->  \nFigure 10. Results of the Late Infall model parameter study, featuring a one-sided torque Γ 1s (here expressed in terms of a time-scale 𝜏 1s) acting on planet c. Different panels refer to different values of the cavity eccentricity enhancement factor 𝐴 𝑒 . In each panel, the 𝑥 -axis denotes the semi-major axis damping time-scale acting on planet c in the cavity. The 𝑦 -axis denotes the eccentricity-damping efficiency (Eq. (7)). The nominal migration time 𝜏 a GLYPH<150> GLYPH<8> is taken to be 16 kyr. The grey dashed line is 𝜏 1s GLYPH<2> 𝐶 𝑒 = 2 GLYPH<149> 5 GLYPH<2> 10 5 yr. Different markers represent different final configurations. The markers '1', '2', '3' and the green dot represent respectively the Φ 1, Φ 2, Φ 3 and Φ 4 end stages, as shown in the four panels in Fig. 9. Furthermore, grey crosses mean that planet c crosses Φ 4, red crosses that at least one of outer planet pair breaks the desired MMRs, i.e., Objective I is not yet achieved, and red squares mean that planet d crosses the disc inner edge. Finally, light-green dots indicate at least one configuration ends up in Φ 4 in simulations with larger 𝜏 a GLYPH<150> GLYPH<8> .  \n<!-- image -->  \nbreak the first three 3BRs but not Φ 4 . As a comprehensive (analytical) investigation of 3BR is beyond the scope of this work (but see some comments in Sect. 6), we instead solve the problem numerically in Sect. 4.1.That is, we vary the torque on planet c, as well as the disc parameters, to assess the conditions for which planet c crosses Φ 1 , Φ 2 and Φ 3 but stays in Φ 4 . Based on these findings, in Sect. 4.2, we implement a complete scenario for the emergence of the TRAPPIST1 dynamical configuration. In Sect. 5.1, we comment on the nature of the inward torque that planet c experiences.", '4.1 Crossing and trapping of the three-body resonances ( dCcgIcT )': "We conduct a parameter study to assess the conditions necessary for the above described scenario, varying the torque on planet c 𝜏 GLYPH<0> 1 1s , the nominal migration time 𝜏 𝑎GLYPH<150> GLYPH<8> (see Eq. (6)), the eccentricityto-semi-major axis damping parameter 𝐶 𝑒 (see Eq. (7)), and the eccentricity-damping boost at the disc inner edge ( 𝐴 𝑒 ; see Eq. (12)). The numerical simulations are conducted with all seven planets: five in the disc (d, e, f, g, and h) and two (b and c) in the cavity. All planets are initialized with a configuration that is consistent with Objective I. Next, we employ a negative torque on planet c in the cavity to break the 3BRs. We record the outcome of the simulation (see below). We conduct a parameter study varying 𝜏 1s in the range » 10 5 GLYPH<150> 10 6 … yr, 𝜏 a GLYPH<150> GLYPH<8> in the » 12 GLYPH<150> 36 … kyr range, 1 GLYPH<20> 𝐴 𝑒 GLYPH<20> 40 and 0 GLYPH<149> 1 GLYPH<20> 𝐶 𝑒 GLYPH<20> 1. Each simulation is integrated for GLYPH<25> 𝜏 1s GLYPH<157> 3, such that it is long enough for 𝑃 d GLYPH<157> 𝑃 c to increase from 1.5 to 1.667 if planet c is not trapped.  \nEscape from 3BR is qualitatively different than escape from twobody MMR. In the two-body case, under the divergent situation described above, planet c will always move away from planet d; we have verified that the critical torque is smaller than 10 GLYPH<0> 7 yr GLYPH<0> 1 . However, the presence of a third resonant body (here planet e or f) qualitatively changes the picture. Even though the torque acting on planet c is divergent, it does not necessarily break the resonance. For example, planet c, d, and e are connected by the 3BR Φ 1 , corresponding to the  \n3:2 resonance locations of both planet c & d and d & e. Therefore, an increase in the period ratio 𝑃 d GLYPH<157> 𝑃 c, due to a negative (divergent) torque acting on planet c, will likewise result in an increase of 𝑃 e GLYPH<157> 𝑃 d , if Φ 1 does not break. In the scenario that planet d and e are in the disc, Type I (convergent) migration, however, prevents the 𝑃 e GLYPH<157> 𝑃 d expansion. Therefore, 𝑃 d GLYPH<157> 𝑃 c tends to stay near 1.5.  \nThe results of our parameter study are displayed in Fig. 10. We present the final configurations from the simulations in which 𝜏 a GLYPH<150> GLYPH<8> = 16 kyr, and vary 𝐴 𝑒 (panels), 𝜏 1s ( 𝑥 -axis) and 𝐶 𝑒 ( 𝑦 -axis). Empirically, as 𝐶 𝑒 increases from 0.1 to 1.0, Φ 2 becomes stronger and is more capable to trap planet c in the cavity (MMRs hold similar characteristics, which is analyzed in Sect. 6.1). It is found that planet c will always strand in either Φ 1 or Φ 2 when the combination 𝜏 1s GLYPH<2> 𝐶 𝑒 = 𝜏 c exceeds a certain value. Therefore, we limit 𝜏 c = 2 GLYPH<149> 5 GLYPH<2> 10 5 yr, which is indicated by the grey dashed line in Fig. 10.  \nDifferent markers represent different final configurations in the simulations. At the end of each simulation, we first check whether the outer planets are still in the desired MMRs; if not we mark them by red crosses. Second, we check whether planet d enters the disc cavity; if so we mark them by red squares. In both scenarios, Objective I has failed to materialize. Most red crosses and red squares are in the panel where 𝐴 𝑒 = 1 and 𝐶 𝑒 is relatively high, because the eccentricity of planet d is not sufficiently damped. Then, planet d will cross the disc inner edge when its eccentricity becomes as large as the migration barrier width ( 𝑒 ¡ GLYPH<24> Δ ), or planet e crosses the 3:2 MMR due to perturbation from eccentric planet d. If planet c is trapped by Φ 1 , Φ 2 , Φ 3 or Φ 4 , we mark the corresponding simulations with '1', '2', '3' or a green dot respectively. Otherwise, if planet c crosses all four 3BR listed here, we mark the outcome by a grey cross.  \nFigure 10 provides some insights into the physical properties of 3BRs. In each panel, planet c tends to be trapped in Φ 1 (marker '1') at high 𝐶 𝑒 . As 𝐶 𝑒 decrease down to 0.1, Φ 2 (marker '2') replaces it in most simulations. It indicates that Φ 1 becomes weaker if more efficient eccentricity damping is adopted. When increasing the one-  \nFigure 11. The period ratio and 3BR angles evolution in one of the systems drawn from Fig. 10, where 𝐴 𝑒 = 15, 𝐶 𝑒 = 0 GLYPH<149> 6 and 𝜏 1s = 0 GLYPH<149> 1 Myr. In the left panel, the period ratios 1.66, 1.60, 1.5, 1.33 are indicated by the four dashed horizontal lines from the upper to the bottom. In the middle panel, we label the expression of each 3BR angle. In the right panel, we zoom in on the period ratios 𝑃 c GLYPH<157> 𝑃 b and 𝑃 d GLYPH<157> 𝑃 c around the time where the b/c/d 3BR breaks.  \n<!-- image -->  \nplanet c in the cavity, planet c breaks through Φ 2 and gets trapped in Φ 4 (green dots) in several simulations. These simulations reveal two points. First, the boundary between trapping in Φ 2 and Φ 4 is not so clear, hinting that the trapping process of planet c is stochastic. Second, planet c is either trapped by Φ 4 or just crosses Φ 4 after Φ 2 , rather than be trapped by Φ 3 . It implies that Φ 3 is much weaker than Φ 2 and Φ 4 . Furthermore, Φ 2 and Φ 4 may have comparable strengths to trap planet c . Otherwise, more green dots would have popped up at the locations now occupied by grey crosses.  \nWe give one possible reason for the stochastic nature of 3BR trapping/escaping. Apart from Φ 1 , Φ 2 , Φ 3 and Φ 4 , there are 9 GLYPH<2> 4 other 3BRs involving the outer five planets that could trap planet c at the four period ratios 1 GLYPH<149> 5, 1 GLYPH<149> 556, 1 GLYPH<149> 6, 1 GLYPH<149> 667. For instance, Φ 3 GLYPH<150> 6 GLYPH<150> 3 c GLYPH<150> d GLYPH<150> g is able to trap planet c at 𝑃 d GLYPH<157> 𝑃 c = 1 GLYPH<149> 667. Some of these 3BRs feature several different libration centers when differing the initial conditions slightly (Tamayo et al. 2017). All these different 3BRs or the same 3BR with different libration centers may have distinct strengths to trap planet c in the cavity. If the initial condition is different, planets may end up in a different configuration, thus contributing to the stochastic behavior.  \nAt the beginning of all of these simulations, 𝑃 c GLYPH<157> 𝑃 b increases in sync with 𝑃 d GLYPH<157> 𝑃 c because 3BR 'bind' planet b/c/d together. However, in most of the simulations where planet c is trapped by Φ 4 as a final state, 𝑃 c GLYPH<157> 𝑃 b is slightly smaller than 8/5. We present an example in Fig. 11. It shows that when 𝑃 d GLYPH<157> 𝑃 c is about to arrive at 5/3, planet c undergoes an inward resonant repulsion (Lithwick & Wu 2012) which breaks the b/c/d resonance. Thereafter, 𝑃 c GLYPH<157> 𝑃 b decreases to 1 GLYPH<149> 595 as shown on the right panel in Fig. 11. This intermediate configuration is a natural output from our model, which is regarded as an initial configuration by Papaloizou et al. (2018). Long term tidal dissipation then expand 𝑃 c GLYPH<157> 𝑃 b back to GLYPH<25> 8 GLYPH<157> 5(Sect. 4.2, cf. Papaloizou et al. 2018), which is consistent with the current observation.  \nDisc evolution may also help to move the TRAPPIST-1 planets into their observed configuration. We run five sets of simulations similar to Fig. 10, but with 𝜏 a GLYPH<150> GLYPH<8> equal to 20 kyr, 24 kyr, 28 kyr, 32 kyr and 36 kyr (not shown). In several of these simulations with higher 𝜏 a GLYPH<150> GLYPH<8> planet c becomes trapped by Φ 4 and we mark them by light-green dots in Fig. 10 together. Therefore, if disc evolution operates during the migration planet c may avoid trapping in Φ 1 and Φ 2 to end up in Φ 4 .  \nIn this model, we have successfully reproduced the orbital configuration of TRAPPIST-1 planets as well as their resonant relationship. However, most of the green dots feature low 𝐶 𝑒 ( GLYPH<159> 0 GLYPH<149> 7) and ex-  \ny high 𝐴 𝑒 ( glyph[greaterorsimilar] 10). Besides, the one-sided torque on planet c in the cavity also needs to be tuned. Considering all these caveats, the probability of TRAPPIST-1 analogs produced by this scenario is not very compelling.", '4.2 Long term tidal dissipation ( dCcgIcT )': "The breaking of 3BR of planets b, c, and d is a natural outcome in dCcgIcT model, as described in Fig. 11. The 3BR of the innermost three planets is believed to librate (Luger et al. 2017; Agol et al. 2021; Teyssandier et al. 2021). To reform this 3BR, long-term tidal dissipation is needed.  \nWe therefore run a long-term simulation, similar to Fig. 5 for one particular parameter set (the same as the one in Fig. 11). We first disperse its disc on a e-folding time-scale of 10 5 yr and run the simulation for 20 𝜏 d . Using the final snapshot from the disc dispersing simulation as an initial condition, we start long term tidal dissipation process. We take 𝑄 sim = 0 GLYPH<149> 1 for all planets to accelerate the simulation by 10 3 times (Papaloizou et al. 2018).  \nSince the disc dispersal does not modify the orbital properties of the planets in this system, we only present the evolution of the planets after removing the disc in Fig. 12. The period ratio 𝑃 c GLYPH<157> 𝑃 b increases with time due to tidal damping during the first GLYPH<25> 0 GLYPH<149> 2 Myr simulation time ( GLYPH<25> 0 GLYPH<149> 2 Gyr physical time). After that, the period ratio 𝑃 c GLYPH<157> 𝑃 b reaches 1.6 and Φ 2 GLYPH<150> 5 GLYPH<150> 3 b GLYPH<150> c GLYPH<150> d reforms. Then, the period ratio increase of the inner planets slows down because all planets are linked by a chain of 3BRs. At 𝑡 GLYPH<25> 2 GLYPH<149> 5 Myr simulation time ( 𝑡 GLYPH<25> 2 GLYPH<149> 5 Gyr physical time), the eccentricities and period ratios of all planets are consistent with the inferred values by Agol et al. (2021) indicated by the 1 𝜎 error bars in Fig. 12. The libration centre of all five adjacent 3BRs lie within the observationally-inferred values as well.  \nOur simulation results can constrain the tidal quality factor of TRAPPIST-1 b. The estimated age of TRAPPIST-1 is 7 GLYPH<149> 6 Gyr (Burgasser & Mamajek 2017), suggesting that the simulation in Fig. 12 evolves GLYPH<25> 3 times faster than the real TRAPPIST-1 system. We also run Stage III simulations that apply tidal torques on (i) only planet b; (ii) only planet b and c. The results are consistent with the simulation that applied the tidal torque on all seven planets. The tidally induced expansion is therefore dominated by the tidal dissipation of planet b (cf. Papaloizou et al. 2018). Since 𝑄 sim = 3 𝑄 GLYPH<157> ' 2 𝑘 2 ' , the chosen simulation constrains the tidal quality factor of planet b as 𝑄 b GLYPH<25> 200 𝑘 2 . This value falls within the upper limit of 5 GLYPH<150> 000 𝑘 2 calculated by Brasser et al. (2019). Note that the bulk planetary eccentricities just after removing the surrounding disc are proportional to  \nFigure 12. Long-term tidal damping simulations of the post-disc phase, illustrating how, after the completion of Objective II, the system evolves towards its present configuration. We present period ratios (top panel), eccentricities (bottom-left) and 3BR angles (bottom-right) during tidal dissipation (with 𝑄 sim = 0 GLYPH<149> 1), after dispersing the gas disc of the system in Fig. 11. The figure style is the same as Fig. 5. We label the observed range of 3BR angles (Luger et al. 2017) on the lower right panel with blue error bars. Tidal damping reforms the Φ 2 GLYPH<150> 5 GLYPH<150> 3 b GLYPH<150> c GLYPH<150> d 3BR and damps the eccentricities towards the observationally-inferred values.  \n<!-- image -->  \np 𝐶 𝑒 (Goldreich & Schlichting 2014; Teyssandier & Terquem 2014; Terquem & Papaloizou 2019). The time one system spends in the tidal dissipation simulation to fit the observed TRAPPIST-1 system is positively correlated with the initial bulk eccentricities. Our model features that 𝐶 𝑒 glyph[lessorsimilar] 0 GLYPH<149> 6, and therefore the corresponding constraint on the tidal quality factor is 𝑄 b glyph[greaterorsimilar] 200 𝑘 2 .  \nIn Fig. 12, long term stellar tidal damping does not alter the period ratios very much compared to Fig. 5. Since the initial eccentricities of planets b and c are ten times smaller than in Fig. 5, stellar tidal forces on planets b and c are also smaller (Eq. (13)) hence dissipating the energy of the system on a longer time-scale. Therefore, the departures from exact commensurability increase more slowly in Fig. 12.", '5 EARLY CAVITY INFALL ( dCcgIcL )': 'In the preceding, we have assumed that planets b and c crossed the cavity only after the emergence of a first order resonance chain, separating Objective I and II. Objective II is the most challenging to meet, as this requires planet c to cross the Φ 1 -Φ 3 3BRs but not to overshoot Φ 4 . However, in reality these stages may be intertwined, such that planets b and c may enter the cavity before the outer resonance chain has established, that is, before (some of the) Φ 𝑖 even exist. For this scenario, we hypothesize that planet b and c, being larger, were able to open partial gaps in the disc that caused them to be disconnected from the disc (Ataiee & Kley 2021), whereas planet d, being smaller, did not suffer this fate and stayed at the cavity edge. While in the cavity we assume, as before, that planet c experiences an inward migration torque, which causes it to drift away slowly from planet d.  \nWe only concentrate on how likely it is to evolve to Objective II from Objective I, because our analysis in Sect. 3.4 (see Fig. 8) is  \nstill roughly applicable. The model is illustrated in Fig. 13. At 𝑡 = 0, we assume that planets b and c and d are initially in 3:2 MMR with 𝑑 at the cavity edge. Subsequently, planets e, f, g, and h enter the resonant chain one after another, each separated by an interval time 𝜏 pl - the planet formation time. Meanwhile planet c experiences a negative torque, expressed in terms of a migration time-scale 𝜏 1s . We now introduce a specific model for the nature of this "one-sided" torque (Sect. 5.1), linking it to the the scaling of nominal Type I semi-major axis damping time-scale 𝜏 𝑎GLYPH<150> GLYPH<8> , i.e., 𝜏 1s is no longer a free parameter. Furthermore, we no longer vary the eccentricity-damping enhancement parameter at the cavity edge, 𝐴 𝑒 = 1, leaving us to explore 𝐶 𝑒 , 𝐴 𝑎 , 𝜏 𝑎GLYPH<150> GLYPH<8> , and 𝜏 pl in Sect. 5.2.', '5.1 Disc repulsion on planets in the cavity': "The Lindblad torque is capable to drive planets in the cavity inward. It arises due to gravitational exchange at the Lindblad resonance locations in the disc (Ward 1997). In the case that the planet is located in the inner disc cavity, only those outermost Lindblad resonances, which still lie in the disc will contribute. Because the nominal migration rate (for an 1 Earth mass planet) and the distant migration rate both scale with disc mass, the rates are linked through:  \n𝜏 GLYPH<0> 1 1s GLYPH<150> GLYPH<8> = 𝜏 GLYPH<0> 1 a GLYPH<150> GLYPH<8> ∑︁ 𝑚 in disc 𝑓 𝑚 = 0 GLYPH<149> 0054 𝜏 GLYPH<0> 1 𝑎GLYPH<150> GLYPH<8> (15)  \nwhere 𝑓 𝑚 is the torque, arising from the 𝑚 th order Lindblad resonance, normalized to the nominal Type I migration torque ( 𝐿 GLYPH<157> 2 𝜏 𝑎GLYPH<150> GLYPH<8> , where 𝐿 is the angular momentum of the planet), and the last expression assumes only 𝑚 = 1 contributes (see below). The Lindblad resonance condition for an orbit of 𝑚 th order torque reads:  \n𝑚 » Ω ' 𝑟 ' GLYPH<0> Ω p … = GLYPH<6> 𝜅 √︁ 1 , 𝑚 2 ℎ 2 GLYPH<150> (16)  \nFigure 13. Schematic of the Early Cavity Infall model ( dCcgIcL ). At 𝑡 = 0 (upper panel), planets b and c have already entered the cavity and planet d has migrated to the inner disc edge. Planet b, c, and d are in 3:2 MMRs. Planet e is still forming or migrating (on time-scale 𝜏 𝑎 ) far out in the disc. Planet c is pushed by the external disc inward on time-scale 𝜏 1s. At 𝑡 = 𝜏 pl (middle panel; 𝜏 pl is the planet formation interval time) planet e enters the resonant chain and is inserted inside the 2:1 MMR of planet d. At this point, planet f is forming or migrating in the outer disc. At 𝑡 = 2 𝜏 pl (bottom panel), planet f ends up inside the 2:1 MMR of planet e. Similarly, planet g and h are inserted after every 𝜏 pl , inside 3:2 MMR and 2:1 MMR respectively. In each panel, we zoom in on the region around planet c and indicate the four 3BRs introduced in Sect. 4 that can trap planet c by solid vertical lines. At early times some of the Φ 𝑖 may not yet exist (dashed vertical lines; unlabelled), allowing planet c to cross these locations.  \n<!-- image -->  \nwhere Ω ' 𝑟 ' is the Keplerian velocity at distance 𝑟 , Ω p is the Keplerian velocity on the planet orbit and 𝜅 is the epicycle frequency. We follow the calculation by Ward (1997) and plot the value of 𝑓 𝑚 for 𝑚 from 1 to 50 in Appendix A.  \nIn the case of the TRAPPIST-1 planets, after planets b and c have entered the cavity, planet d resides near the cavity radius location. As the period ratios of planets d and b are larger than 2:1, there is no one-sided Lindblad torque from the outer disc on planet b. On the other hand, the period ratio of planets c and d exceeds 1.5 but is less than 2. Therefore, only the outer 2:1 Lindblad resonance of planet c falls in the disc and the relevant order is 𝑚 = 1. The value of 𝑓 1 for the outer Lindblad resonance is 0.0054. For instance, for 𝜏 𝑎 = 10 kyr we arrive at 𝜏 1s = 1 GLYPH<149> 85 Myr. Although 𝜏 1s is long, the time to cross a distance from -say- the 3:2 to the 5:3 resonance location is lower by about a factor GLYPH<24> 10. Therefore 𝜏 1s may compete with the planet formation interval time-scale 𝜏 pl .", '5.2 Parameter study ( dCcgIcL )': "We conduct a parameter study to assess the conditions needed to obtain the desired resonance configuration. Apart from varying the interval time 𝜏 pl from 20 to 200 kyr, we also vary the eccentricityto-migration damping parameter 𝐶 𝑒 from 0.05 to 1, the migration threshold barrier parameter 𝐴 𝑎 from 50 to 150, and the nominal planet migration parameter 𝜏 𝑎GLYPH<150> GLYPH<8> from 4 kyr to 22 kyr. We let the system evolve until GLYPH<24> 𝜏 1s (Eq. (15)). The results are shown in Fig. 14. We mark the final configurations following the same style as Fig. 10. Wealso mark the configurations where planet c stays in Φ 4 but where the period ratio has deviated by some amount from 8:5 ( 𝑃 c GLYPH<157> 𝑃 b GLYPH<159> 1 GLYPH<149> 53 or 𝑃 c GLYPH<157> 𝑃 b ¡ 1 GLYPH<149> 63) with grey dots.  \nIn this Early Cavity Infall model the new parameter 𝜏 pl sensitively determines the resonant structure. Although we assume 𝜏 pl is the  \nsame for all outer planets, what matters chiefly is the arrival time of planets e and f. If 𝜏 pl is too long, planet c would migrate so far such that 𝑃 d GLYPH<157> 𝑃 c ¡ 5 GLYPH<157> 3 before planet e enters the resonant chain. Our study does not cover such long 𝜏 pl . In contrast, if 𝜏 pl is too short, planet c would not migrate very far with 𝑃 d GLYPH<157> 𝑃 c GLYPH<159> 14 GLYPH<157> 9 after the arrival of planet f, which bring about the existence of Φ 2 . Then it is unlikely for planet c to overcome Φ 2 . We are hence searching for the sweet spot where planet c finally stays at Φ 4 .  \nStarting from the bottom-right corner of each panel in Fig. 14, the one-side torque is relatively weak and 𝜏 pl is short. Planet c mostly stays in Φ 1 or Φ 2 . With increasing 𝜏 pl and decreasing 𝜏 𝑎GLYPH<150> GLYPH<8> , planet c can migrate over larger distances in the cavity before the outer planets appear. Hence, planet c increasingly ends in Φ 3 and Φ 4 . Interestingly, in the left two columns green dots are more frequent than marker '3' and there is usually a clear boundary between the two markers. The corresponding one-sided torque at the boundary therefore represents the strength of Φ 3 . It implies that Φ 4 is sturdier than Φ 3 , which provides a wide range of parameter space to overcome Φ 3 but stay in Φ 4 . For this reason, we colour '3' light green - if (for some reason) planet c manages to escape Φ 3 , it is still likely to end up in Φ 4 .  \nWhen 𝜏 𝑎GLYPH<150> GLYPH<8> becomes too small, some planets in the disc will break their desired MMRs. As shown in Fig. 14 when 𝐶 𝑒 GLYPH<159> 0 GLYPH<149> 2, there are some red crosses on the left and a clear boundary between red crosses and green dots. The corresponding disc migration torque at the boundary indicates the critical torque, above which the desired configuration of outer planets (d, e, f, g, and h) will break. As 𝐶 𝑒 increases from 0.05 to 0.1, the desired configuration of outer planets becomes more stable (the boundary between red crosses and green dots moves to lower 𝑡 𝑎 ). This trend is consistent with Fig. 8, and is analysed in Sect. 6.1.  \nHowever, when increasing 𝐶 𝑒 to higher values, the boundary between red crosses and green dots becomes less sharp and the ap-  \nFigure 14. Final orbital configurations produced by the Early Cavity Infall model ( dCcgIcL ). Different columns and different rows refer to different eccentricitydamping parameter 𝐶 𝑒 and migration threshold height 𝐴 𝑎 . In each panel, the 𝑥 -axis denotes the Type I semi-major axis damping time-scale 𝜏 𝑎GLYPH<150> GLYPH<8> , the 𝑦 -axis denotes the planet formation span 𝜏 pl . Different markers represent different final configurations, following the same style as Fig. 10. We additionally mark the configurations in which planet c stays in Φ 4 with 𝑃 c GLYPH<157> 𝑃 b GLYPH<159> 1 GLYPH<149> 53 or 𝑃 c GLYPH<157> 𝑃 b ¡ 1 GLYPH<149> 63 by grey dots.  \n<!-- image -->  \nce of red crosses becomes stochastic. The reason is that with higher 𝐶 𝑒 planet eccentricities are higher, which renders phenomena like resonance crossing more stochastic (Ogihara & Kobayashi 2013). Therefore, we see more and more red crosses distributed in the right columns of Fig. 14 (higher 𝐶 𝑒 ).  \nWe test the influence of the migration threshold barrier 𝐴 𝑎 by three different values, 50, 100, 150 in Fig. 14. If 𝐴 𝑎 is small, it shows more red squares, where planet d enters the cavity, when 𝐶 𝑒 is high. When increasing 𝐴 𝑎 , the results of 𝐴 𝑎 = 100 is largely similar to 𝐴 𝑎 = 150, signifying a threshold above which planet d stays at the disc edge.  \nWefinally comment on the meaning of the grey dots in Fig. 14. As a rule of thumb we see that, similar to Fig. 11, when planet c reaches Φ 4 the 3BR Φ 2 GLYPH<150> 5 GLYPH<150> 3 b GLYPH<150> c GLYPH<150> d breaks, whereafter it will reform due to long term tidal dissipation (Fig. 12). However, it is found that in some of the simulations, the 3BR Φ 2 GLYPH<150> 5 GLYPH<150> 3 b GLYPH<150> c GLYPH<150> d already breaks early. Two physically different mechanisms result in grey dots: (a) If 𝜏 pl is so long that Φ 4 is reached before planets f and g appear, the outer planets f and g will further push planet c slightly inward. This case applies to the top left two panels in Fig. 14, where 𝜏 pl is long and 𝐴 𝑎 is small. (b) The 3BR Φ 2 GLYPH<150> 5 GLYPH<150> 3 b GLYPH<150> c GLYPH<150> d breaks early when Φ 3 is crossed or when a new planet enter the resonant chain. Since planet c is still repulsed inward by the Lindblad torque, 𝑃 c GLYPH<157> 𝑃 b decreases towards GLYPH<25> 1 GLYPH<149> 5. We have run longterm tidal dissipation simulation starting with the analog satisfying 𝑃 c GLYPH<157> 𝑃 b GLYPH<25> 1 GLYPH<149> 53, but found that it fails to reform the 3BR Φ 2 GLYPH<150> 5 GLYPH<150> 3 b GLYPH<150> c GLYPH<150> d .  \nWhen adopting low 𝐶 𝑒 ( glyph[lessorsimilar] 0 GLYPH<149> 20) the low eccentricities of the outer planets render the system more stable, preventing resonance crossing.  \nInterestingly, there is another green dots island on the panel where 𝐶 𝑒 = 1 and 𝐴 𝑎 = 100. We do not have a clear explanation for the occurrence of so many green dots in this panel. Possibly, the strength of Φ 3 is non-linearly dependent on 𝐶 𝑒 , such that it is strongest at 𝐶 𝑒 GLYPH<25> 0 GLYPH<149> 5 and breaks again when 𝐶 𝑒 = 1.  \nTo achieve the observed TRAPPIST-1 configuration, the early cavity infall model compares favourably to the late infall model of Sect. 4.1. The early cavity infall model requires less parameter tuning - we fix 𝐴 𝑒 = 1 and link 𝜏 1 𝑠 to 𝜏 𝑎GLYPH<150> GLYPH<8> - and the values for the 𝐴 𝑎 and 𝜏 1s parameters are natural or in line with formation models (see Sect. 7.1). In particular for low 𝐶 𝑒 , the likelihood of a successful outcome is quite high.", '6.1 Trapping under convergent migration in MMRs': 'Through numerical simulations, we find that more efficient eccentricity-damping promotes resonance crossing (see Fig. 8). The trend can be understood using the pendulum model (Murray & Dermott 1999) for internal MMRs. Assuming small deviations from the equilibrium value for the resonance angle, the libration width of a 𝑗 : 𝑗 GLYPH<0> 1 MMR can be written as:  \nΔ 𝑎 = GLYPH<18> 16 j 𝐶 𝑗 j 𝑒 2 𝑛 GLYPH<19> 1 GLYPH<157> 2 𝑎GLYPH<150> (17)  \nwhere 𝑛 , 𝑒 and 𝑎 are the mean motion, eccentricity and semi-major axis of the inner planet, respectively. The constant 𝐶 𝑗 arises from the  \nresonant part of the disturbing function (Murray & Dermott 1999):  \n𝐶 𝑗 = 𝑚 0 𝑀 ★ 𝑛𝛼 𝑗 𝑓 𝑗 (18)  \nwhich is proportional to the inner planet mass 𝑚 0 . In Eq. (18), 𝛼 𝑗 and 𝑓 𝑗 are the inner-to-outer planet semi-major axis ratio and the direct term in the expansion of the disturbing function. The value of the most commonly used 𝑓 𝑗 can be seen on page 334 in (Murray & Dermott 1999). The libration time-scale is:  \n𝜏 lib = GLYPH<18> 2 𝜋 3 j 𝐶 𝑗 j 𝑛𝑒 GLYPH<19> 1 GLYPH<157> 2 GLYPH<149> (19)  \nFollowing Ogihara & Kobayashi (2013), we write the critical migration time-scale 𝜏 a GLYPH<150> crit as:  \n𝜏 a GLYPH<150> crit = 𝜏 lib 𝑎 Δ 𝑎 = 2 𝜋 16 𝐶 2 𝑗 𝑒 2 ! 1 GLYPH<157> 2 GLYPH<150> (20)  \nbelow which the corresponding MMR will break. In the case of convergent migration of two planets, the equilibrium eccentricities of planets balancing eccentricity-damping from the disc and excitation by resonant interaction reads:  \n𝑒 2 GLYPH<24> 𝜏 e 𝜏 a (21)  \n(Goldreich & Schlichting 2014; Teyssandier & Terquem 2014; Terquem & Papaloizou 2019). Combining Eq. (20) and Eq. (21) we obtain:  \n𝜏 a GLYPH<150> crit GLYPH<24> 2 𝜋 16 𝐶 2 𝑗 𝜏 e GLYPH<149> (22)  \nThis expression implies that we can either employ efficient eccentricity-damping (low 𝜏 e) or use smaller planet mass (low 𝐶 𝑗 ), in order to enlarge 𝜏 a GLYPH<150> crit and cross the MMR. Note that although the above expressions are derived for an internal MMRs (where the outer planet is fixed), those for the external MMRs are almost the same (Murray & Dermott 1999).  \nThe above analytic results is consistent with Fig. 8. The value of 𝜏 𝑎GLYPH<150> GLYPH<8> at the boundary between triangles and circles indicate the critical migration time-scale 𝜏 𝑎GLYPH<150> GLYPH<8> of the 2:1 MMR of planet b & c. As 𝐶 𝑒 decreases or 𝐴 𝑒 increases, the boundary moves to the right, implying longer 𝜏 𝑎GLYPH<150> GLYPH<8> .', '6.2 Eccentricity-period relation in 3BRs': "In a three-planet system connected by one 3BR with each of the adjacent two planets connected by 𝑗 GLYPH<0> 1: 𝑗 and 𝑘 GLYPH<0> 1: 𝑘 MMRs, the eccentricities and period ratios are related to each other. We define the distance to exact commensurability of planet 𝑖 and 𝑖 , 1 as:  \n𝛿 𝑖GLYPH<150> 𝑖 , 1 = ' 𝑗 GLYPH<0> 1 ' 𝑃 𝑖 , 1 𝑗 𝑃 𝑖 GLYPH<0> 1 GLYPH<150> (23)  \nwhere 𝑃 𝑖 is the orbital period of planet 𝑖 . The MMR angles in a planet pair near 𝑗 GLYPH<0> 1: 𝑗 MMRare:  \n𝜙 𝑖GLYPH<150> 𝑖 , 1; 𝑖 , 𝑋 = ' 𝑗 GLYPH<0> 1 ' 𝜆 𝑖 GLYPH<0> 𝑗𝜆 𝑖 , 1 , 𝜛 𝑖 , 𝑋 GLYPH<149> (24)  \nwhere 𝑋 = 0 or 𝑋 = 1. The time derivatives / 𝜙 𝑖GLYPH<150> 𝑖 , 1; 𝑖 , 𝑋 average to zero if 𝜙 𝑖GLYPH<150> 𝑖 , 1; 𝑖 , 𝑋 is librating. Taking the time derivative of Eq. (24) and combining it with Eq. (23), we have:  \n𝛿 𝑖GLYPH<150> 𝑖 , 1 = GLYPH<0> / 𝜛 𝑖 , 𝑋 𝑃 𝑖 , 1 2 𝜋𝑗 GLYPH<149> (25)  \nA similar relation holds for planet 𝑖 , 1 & 𝑖 , 2. Without external perturbation the three planets are locked in common apsidal precession, / 𝜛 𝑖 = / 𝜛 𝑖 , 1 = / 𝜛 𝑖 , 2 . Then, the ratio of the distance to exact commensurability of two planet pairs is:  \n𝛿 𝑖 , 1 GLYPH<150> 𝑖 , 2 𝛿 𝑖GLYPH<150> 𝑖 , 1 = 𝑗 𝑃 𝑖 , 2 𝑘𝑃 𝑖 , 1 GLYPH<25> 𝑗 𝑘 GLYPH<0> 1 GLYPH<149> (26)  \nIn Laplace resonances where 𝑗 = 𝑘 , the right hand term of Eq. (26) is larger than unity. Ramos et al. (2017) and Teyssandier & Libert (2020) insert Lagrange's planetary equation, / 𝜛 = 𝐶 𝑗 cos 𝜙 GLYPH<157> 𝑒 , into Eq. (25) and get:  \n𝛿 𝑖GLYPH<150> 𝑖 , 1 = 1 𝑗 𝛼 GLYPH<0> 1 GLYPH<157> 2 𝑗 𝑓 𝑗 𝑚 𝑖 , 1 𝑀 ★ 1 𝑒 𝑖 (27)  \nCombining Eq. (26) and Eq. (27), we arrive at:  \n𝑒 𝑖 𝑒 𝑖 , 1 = 𝑚 𝑖 , 1 𝑓 𝑗 𝑚 𝑖 , 2 𝑓 𝑘 GLYPH<18> 𝛼 𝑘 𝛼 𝑗 GLYPH<19> 1 GLYPH<157> 2 𝑃 𝑖 , 2 𝑃 𝑖 , 1 GLYPH<25> 𝑚 𝑖 , 1 𝑓 𝑗 𝑚 𝑖 , 2 𝑓 𝑘 GLYPH<20> 𝑘 ' 𝑗 GLYPH<0> 1 ' 𝑗 ' 𝑘 GLYPH<0> 1 ' GLYPH<21> 1 GLYPH<157> 3 𝑘 𝑘 GLYPH<0> 1 GLYPH<149> (28)  \nIn Laplace resonances, the right hand term of Eq. (28) is larger than unity when planets have nearly the same mass. Therefore, outer planets tend to further depart from exact commensurability (Eq. (26)) and have smaller eccentricities (Eq. (28)). The eccentricities and the period ratios relationships of model bC , in which planets stay in first-order MMR, is revealed by Fig. 5. The deviations from exact commensurability increase and the eccentricities decrease from planet b to planet h. In the final stage of dCcgIcT (Stage III in Fig. 12, planet b & c and c & d are near MMRs of order higher than one. Because the analysis above only holds for first-order resonances, the relationship between the eccentricity and period ratios can be seen for outer planets but not for planet b/c/d.", '7.1 Model Assessment': "In this work, we have aimed to connect the present-day dynamical properties of the TRAPPIST-1 planets to their formation. We have assumed that the planets are formed at locations further away from their present orbits, were subject to inward (Type I) migration, which was halted at the disc inner edge. Subsequently, the planets formed a chain of first-order MMRs, predominantly of 3:2 commensurability, and a plethora of three-body resonances (Stage I). In the next Stage II, which is either separated in time from Stage I (Sect. 4) or occurs contemporaneously with Stage I (Sect. 5), the two inner-most planets fell into the cavity, allowing their period ratios to expand due to the (weak and distant) torque from the disc. Finally, stellar tidal interactions after the demise of the gas disc (Stage III) provided the final act in moving the planets into a configuration consistent with their present-day orbital periods, period ratios, eccentricities, and three-body resonance angles.  \nThe first successful model matching the observed TRAPPIST-1 configuration has Stage I, II, and III occurring separately (Sect. 4). However, the model has several drawbacks: (i) The probability of a successful case is low. (ii) It requires rather high additional eccentricity-damping on the planet at the disc inner edge ( 𝐴 𝑒 GLYPH<29> 10). Intriguingly, such enhanced eccentricity-damping at the disc edge has recently been observed (but not quite understood) in hydrodynamical simulations of planet-disc interaction (Ataiee & Kley 2021). Magnetospheric effects may also contribute towards a pileup of material near the cavity radius (D'Angelo & Spruit 2010), which would  \nbe the most natural explanation for the enhanced damping. Therefore, this model demands a very high surface density at the edge of the disc or some unknown mechanisms that would give rise to the anomalous efficient eccentricity-damping around the truncation radius. (iii) The physical origin of the one-sided negative torque on the planet in the disc cavity is not motivated. The torque required is usually GLYPH<24> 10 times stronger than the Lindblad torque calculated in Sect. 5.1. Hence, for this scenario to work, the surface density of the gas needs to be boosted by a factor GLYPH<24> 10 near the cavity region. Even then, the probability of an outcome consistent with the TRAPPIST-1 system remains low.  \nFinally, we have investigated the Early Cavity Infall model (Sect. 5) which features planets b and c entering the cavity before the outer planets park in resonance near the inner disc, combining Stage I and Stage II. In this scenario, we adopt a value for the Lindblad torque acting on planet c consistent with the disc mass. The key drawback of this model is that it introduces another parameter - the time interval between planet formation ( 𝜏 pl ) - which needs to be tuned such that it is at least one-tenth of the time-scale over which planet c is being pushed in. 1 Only then can planet c migrate sufficiently far inwards to cross the locations of the 3BRs before their appearance. In addition, success is more likely for low values of the eccentricity damping parameter 𝐶 𝑒 . Nevertheless, this scenario is attractive for the following reasons: (i) the planet formation interval time ( 𝜏 pl ) of GLYPH<24> 10 5 yr is in line with the pebble-driven growth scenario that we have previously explored (Ormel et al. 2017; Schoonenberg et al. 2019). (ii) Other parameters are also natural or in line with the constraints from Objective I. No anomalous eccentricity damping is present ( 𝐴 𝑒 = 1), the value of 𝐴 𝑎 is in line with theory (Liu et al. 2017), and the values of 𝜏 𝑎GLYPH<150> GLYPH<8> that are required are overall consistent with the Objective I results. Hence, we favour this scenario to explain the architecture of the the TRAPPIST-1 system.  \nIn the Early Cavity Infall model (Fig. 14), efficient eccentricitydamping is preferred to prevent the outer planets, especially planet d, from crossing the desired resonances due to the stochastic nature of resonance trapping at high eccentricity. For the TRAPPIST-1 system, such high eccentricity-over-semi-major axis damping (low 𝐶 𝑒 ) is also used in other literature studies, e.g., Teyssandier et al. (2021) and MacDonald et al. (2021). The most straightforward explanation is that Type I disc migration, being a complex phenomenon, operated less efficient in TRAPPIST-1, i.e., that 𝛾 𝐼 is smaller than what linear theory prescribes. In addition, subdued Type-1 migration would also ensure 𝐴 𝑎 will be above the required threshold. A low 𝐶 𝑒 parameter is also required in the K2-24 planetary system to reproduce its observed TTV signals (Teyssandier & Libert 2020).", '7.2 Resonances Establishing and Breaking': 'Three-body resonances play a crucial role in shaping these scenarios. A multi-planet resonant system as TRAPPIST-1 is characterized by many adjacent and non-adjacent three-body resonances, with some of the Φ 𝑖 librating, while others do not. This dynamical imprint strongly constrains the formation history of the TRAPPIST-1 system. In 3BRs, resonance locking limits the ability for period expansion. In the literature, certain mechanisms that have been invoked to separate two-planet systems by divergent migration, e.g., stellar tidal damping (Lithwick & Wu 2012; Batygin & Morbidelli 2013) or magnetospheric rebound (Liu et al. 2017). When, however, the  \nplanets are instead connected to other planets through three-body resonances, all planets participate in the expansion and the effects on -say- the inner two planets are potentially limited. Specifically, for TRAPPIST-1 the Φ 2 GLYPH<150> 5 GLYPH<150> 3 c GLYPH<150> d GLYPH<150> e resonance, corresponding to the 3:2 and 3:2 MMR of planets c/d/e, naturally forms during the migration phase (Stage I), but greatly hinders the further period expansion of the inner planet system towards their observed period ratios. As this angle is currently not librating, we conclude that it must have been broken early, which would be achieved when planets b and c fell into the cavity and were repulsed by the gas disc. This allows the subsequent expansion of the inner planet system b/c/d towards their observed orbital period ratios.  \nThis (temporary) decoupling of the inner planet system from the other TRAPPIST-1 planets was already employed in the study by Papaloizou et al. (2018), where initially two dynamically separate subsystems (b/c/d and e/f/g/h), each locked in their respective threebody resonances, dynamically converge into the present configuration under the action of stellar tidal forces only. As in that scenario stellar tides would primarily operate on the inner system, it naturally explains its expansion. However, Papaloizou et al. (2018) do not motivate the origin of the two subsystems, their initial separations and the initial period ratios within each system. Indeed, the initial setup of these simulations follows from the formation stage. We find that during the expansion, planets may well enter the wrong three-body resonance (especially Φ 2 GLYPH<150> 5 GLYPH<150> 3 c GLYPH<150> d GLYPH<150> e , Φ 3 GLYPH<150> 15 GLYPH<150> 12 c GLYPH<150> e GLYPH<150> f and Φ 5 GLYPH<150> 14 GLYPH<150> 9 c GLYPH<150> d GLYPH<150> e ; see Fig. 9, Fig. 10 and Fig. 14), which stellar tides alone is too weak to break. This supports our scenario where the inner-most planets were pushed into the disc cavity.  \nTo achieve a first-order resonant chain, we need to assume that planets c and g start in a 3:2 and 4:3 resonance (Fig. 8). Coleman et al. (2019) and Burn et al. (2021) also initialize planet embryos close to each other such to avoid trapping in wider MMR like 2:1 and 3:2. A natural explanation is that these planets formed rapidly oneafter-another, which is conceivable in the pebble-driven formation model (Schoonenberg et al. 2019; Lin et al. 2021b).  \nIn this work, we simplify the disc model without complicating discplanets interactions via introducing, e.g., radial gravity (Nagasawa et al. 2003; Pan et al. 2020) and stochastic torque (Rein & Papaloizou 2009; Paardekooper et al. 2013) induced by the disc and raised on planets. How much difference those mechanisms can make to 3BR trapping and escaping need to be addressed in the future.', '7.3 Application to other systems': 'In our model, the planet in the disc cavity can be repulsed inward slowly by Lindblad resonance from the external disc. Three-body resonance with two planets in the external disc can then halt the inward migration in the disc cavity. This mechanims therefore enables a configuration with the inner planets more separated from each other than the outer planets. Although three-body resonances are not ubiquitous in the Kepler data (Goldberg & Batygin 2021), we expect that configurations with separated inner planets near higerorder resonances may be connected by three-body resonances.  \nWe give three examples. HD 158259 is a six planets system (Hara et al. 2020). The period ratios are 𝑃 c GLYPH<157> 𝑃 b = 1 GLYPH<149> 57, 𝑃 c GLYPH<157> 𝑃 b = 1 GLYPH<149> 51, 𝑃 c GLYPH<157> 𝑃 b = 1 GLYPH<149> 53, 𝑃 c GLYPH<157> 𝑃 b = 1 GLYPH<149> 51 and 𝑃 c GLYPH<157> 𝑃 b = 1 GLYPH<149> 44, where planet b and c are more separated than other adjacent plant pairs. Interestingly, planet c, d and e in HD 158259 are near 3:2 mean-motion resonance, similar to planet d, e and f in TRAPPIST-1. Based on our analysis in Sect. 4, we predict that the three-body resonance angle Φ b GLYPH<150> d GLYPH<150> e =  \n3 𝜆 b GLYPH<0> 15 𝜆 d , 12 𝜆 𝑒 in HD 158259 system is librating. This angle is equivalent to Φ 2 in TRAPPIST-1 (Sect. 4 and Sect. 5).  \nOther systems, while not sharing the same period ratios configuration as TRAPPIST-1, do feature wider-spaced inner planets. Recently, Leleu et al. (2021) found that the five outermost planets of the six-planet TOI-178 system constitute a 2:4:6:9:12 chain of Laplace resonances. The period ratio of the second-to-innermost planet is 1 GLYPH<149> 69. Since this period ratio is close to 5:3, it may indicate that planet b is connected by the outer planets through three-body resonance. In the Kepler-80 system, the period ratios of the six planets are 3.1, 1.511, 1.518, 1.35, and 1.538 from inner to outer respectively (Xie 2013; MacDonald et al. 2021). Although the period ratio of 3.1 cannot be explained by the Lindblad torque mechanism, it would still be worthwhile to investigate whether a 3BR trapping may occur at this location. To address what role 3BRs in these systems have played in shaping the configuration of these systems, the first step is to evaluate their strength by numerical simulations such as done in this work.', '8 CONCLUSIONS': 'We conducted N-body simulations to reconstruct the origin of the dynamical structure of planets in the TRAPPIST-1 system. Our model covers the stages directly after planet formation ( GLYPH<24> 10 5 yr) in the proto-planet disc, until the present age. Two objectives were pursued: the formation of a chain of first-order MMR (Objective I) and the evolution of the inner-most planets towards their present dynamical configuration (Objective II). In the Late Infall Model planets b and c entered the disc cavity at some time after all planets entering the resonant chain, and were repulsed by the one-sided Lindblad torque exerted by the exterior disc. Together with stellar tidal damping, this provided for the expansion of the inner sub-system towards its present configuration. In the preferred Early Infall Model, planet b and c fall in the cavity before the outer planets entering the resonant chain. Then, planet c can migrate across the resonant locations in the cavity before their existence.  \nWe list our main findings from this study:  \n- (i) The 4:3 MMR of planet f & g stands out as the closest resonance. In simple disc models planets f & g fail to cross the 3:2 resonance without preventing other planets from obtaining a more compact configuration. The 4:3 MMR could be explained as an imprint of rapid sequential planet formation, whereafter f and g migrated convergently.\n- (ii) Simple scenarios that feature only disc and stellar tidal damping fail to meet Objective II. Due to three-body resonance locking, the inner subsystem cannot expand independently from the rest of the system.\n- (iii) Objective II can be reached when planets b and c entered the disc cavity. The disc further repulsed the planets in the cavity by virtue of the order 𝑚 = 1 outer Lindblad torque whose resonant location still resides in the disc. This repulsion is a crucial mechanism in aligning the inner system towards the observed 3BRs and period ratios.\n- (iv) Efficient eccentricity-damping ( 𝐶 𝑒 GLYPH<24> 0 GLYPH<149> 1) and a strong disc inner barrier ( 𝐴 𝑎 GLYPH<24> 100) were found to be advantageous to obtain the desired resonance configuration. These imply a small value for the Type I migration prefactor 𝛾 𝐼 .\n- (v) In the Late Infall model, planet b and c enter the disc cavity after outer planets enter the resonant chain. To reproduce TRAPPIST1, it is necessary to enhance the eccentricity damping of planet d at the cavity by large amounts ( 𝐴 𝑒 glyph[greaterorsimilar] 10).  \n(vi) In our preferred scenario, planet b and planet c enter the disc cavity early, before the outer planets have joined the resonance chain. To allow planet c to migrate across the barrier resonances, the planet formation interval time 𝜏 pl is about 10 5 yr.  \n(vii) Tidal dissipation decreases the eccentricities and increases the period ratio among planets in a three-body resonance chain. After all planets have been connected by three-body resonances, it is the period ratio of the outermost planet pair (h and g) that increases most, while the rate of the expansion is determined by the tidal dissipation of the innermost planet b. For our model, we constrain the tidal quality factor of planet b 𝑄 𝑏 glyph[greaterorsimilar] 200 𝑘 2 .  \n(viii) Three-body resonances play an important role in shaping the dynamical structure of the TRAPPIST-1 planets. In contrast to two-body resonances, 3BRs do not necessarily break when exposed to a divergent force. In this study, a numerical approach was used to assess the strength of 3BR and we have highlighted the role of non-adjacent 3BR. Due to their significance in shaping scenarios containing divergent migration, a more comprehensive understanding of 3BR is in need.', 'ACKNOWLEDGEMENTS': 'We thank the referee for their constructive suggestions and comments. We thank Eric Agol for providing us with an early version of the transit timing and photometric analysis of TRAPPIST1 planets (Agol et al. 2021) and corrections. We also thank Ramon Brasser, Beibei Liu and Jean Teyssandier for useful email discussions. Software: Matplotlib (Hunter 2007; Caswell et al. 2021), REBOUND (Rein & Liu 2012) and REBOUNDx (Tamayo et al. 2020). We also use Pylaplace to calculate the Laplace coefficients, https://pypi.org/project/pylaplace/ .', 'DATA AVAILABILITY': 'The data underlying this article will be shared on reasonable requests to the corresponding author.', 'REFERENCES': 'Agol E., et al., 2021, The Planetary Science Journal, 2, 1 Artymowicz P., 1993, ApJ, 419, 166  \nAtaiee S., Kley W., 2021, arXiv e-prints, p. arXiv:2102.08612  \nBaruteau C., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 667 ( arXiv:1312.4293 ), doi:10.2458/azu\\_uapress\\_9780816531240-ch029  \nBatygin K., Morbidelli A., 2013, AJ, 145, 1  \nBenítez-Llambay P., Masset F., Koenigsberger G., Szulágyi J., 2015, Nature,  \n520, 63  \nBirnstiel T., et al., 2018, ApJ, 869, L45  \nBrasser R., Barr A. C., Dobos V., 2019, MNRAS, 487, 34  \nBurdanov A. Y., et al., 2019, MNRAS, 487, 1634  \nBurgasser A. J., Mamajek E. E., 2017, ApJ, 845, 110  \nBurn R., Schlecker M., Mordasini C., Emsenhuber A., Alibert Y., Henning  \n- T., Klahr H., Benz W., 2021, arXiv e-prints, p. arXiv:2105.04596  \nCarr J. S., Tokunaga A. T., Najita J., 2004, ApJ, 603, 213  \n- Caswell T. A., et al., 2021, matplotlib/matplotlib: REL: v3.4.2,  \ndoi:10.5281/zenodo.592536  \nCharalambous C., Martí J. G., Beaugé C., Ramos X. S., 2018, MNRAS, 477,  \n1414  \nChristiansen J. L., et al., 2018, AJ, 155, 57  \nColeman G. A. L., Leleu A., Alibert Y., Benz W., 2019, A&A, 631, A7  \nFigure A1. Torque magnitude as a function of order, m, normalized to the Type I migration time-scale for scale height ℎ = 0 GLYPH<149> 03 and disc surface density gradient 𝑘 = 0 GLYPH<149> 5. Left panel present the order m from 1 to 50. Right panel zoom in m from 1 to 5.  \n<!-- image -->', 'APPENDIX A: LINDBLAD TORQUE': 'Following Ward (1997), we calculate the magnitude of the Lindblad torque at different order ( 𝑚 ) of Lindblad resonance. We insert disc scale height ℎ = 0 GLYPH<149> 03, disc power-law surface density gradient 𝑘 = 𝑑 ln Σ g GLYPH<157> 𝑑 ln 𝑟 = 0 GLYPH<149> 5 and disc temperature gradient 𝑙 = d ln 𝑇 GLYPH<157> d ln 𝑟 = 1. The left panel of Fig. A1 shows the Lindblad torque scale by Type I migration torque at the first 50 orders Lindblad resonance location and the right panel only shows the first five orders. For m=1 only an outer Lindblad torque is present, which is used in this work. Its value is 0.0054.  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.'}