The theory of Type I migration has been widely used in many studies. Transiting multiplanet systems offer us the opportunity to examine the consistency between observation and theory especially for those systems harbouring planets in Mean Motion Resonance MMR. The displacement these resonant pairs show from exact commensurability provides us with information on their migration and eccentricitydamping histories. Here we adopt a probabilistic approach characterized by two distributions  appropriate for either the resonant or nonresonant planets  to fit the observed planet period ratio distribution. With the Markov chain Monte Carlo MCMC method we find that approx 15 rm per cent of exoplanets are in firstorder j  1 j MMRs the ratio of eccentricitytosemimajor axis damping is too high to allow overstable librations and that the results are byandlarge consistent with TypeI migration theory. In addition our modelling finds that a small fraction of resonant pairs is captured into resonance during migration implying late planet formation gaspoor. Most of the resonant pairs park themselves at the migration barrier indicating early planet formation gasrich. Furthermore after improving the criterion on twobody resonant trapping we obtain an upper limit of the disc surface density at the time the planets are locked in resonance.
['celestial mechanics', 'planets and satellites: dynamical evolution and stability', 'planets and satellites: formation', 'planet-disc interactions', 'Astrophysics - Earth and Planetary Astrophysics']
When where and how many planets end up in firstorder resonances
{'No Header': ', 1-18 (2022)', 'When, where, and how many planets end up in first-order resonances?': 'Shuo Huang ( 黄 硕 ), 1 GLYPH<150> 2 ★ Chris W. Ormel, 1 y  \n- 1 Department of Astronomy, Tsinghua University, 30 Shuangqing Rd, 100084 Beĳing, China\n- 2 Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands  \nAccepted XXX. Received YYY; in original form ZZZ', 'ABSTRACT': 'The theory of Type I migration has been widely used in many studies. Transiting multi-planet systems offer us the opportunity to examine the consistency between observation and theory, especially for those systems harbouring planets in Mean Motion Resonance (MMR). The displacement these resonant pairs show from exact commensurability provides us with information on their migration and eccentricity-damping histories. Here, we adopt a probabilistic approach, characterized by two distributions - appropriate for either the resonant or non-resonant planets - to fit the observed planet period ratio distribution. With the Markov chain Monte Carlo (MCMC) method, we find that GLYPH<25> 15% of exoplanets are in first order ( 𝑗 , 1: 𝑗 ) MMRs, the ratio of eccentricity-to-semi-major axis damping is too high to allow overstable librations and that the results are by-and-large consistent with Type-I migration theory. In addition, our modeling finds that a small fraction of resonant pairs is captured into resonance during migration, implying late planet formation (gas-poor). Most of the resonant pairs park themselves at the migration barrier, indicating early planet formation (gas-rich). Furthermore, after improving the criterion on two-body resonant trapping, we obtain an upper limit of the disc surface density at the time the planets are locked in resonance.  \nKey words: celestial mechanics - planet-disc interactions - planets and satellites: formation - planets and satellites: dynamical evolution and stability', '1 INTRODUCTION': 'Since the first discovery of exoplanets around solar-type stars (Mayor & Queloz 1995), the number of exoplanets has ballooned in the last three decades, exceeding 5 200 as of the present day. It is therefore appropriate to conduct population-level analyses to examine planet formation theories (Mordasini et al. 2015; Zhu & Dong 2021). When independent mass and radius measurements are available, planet bulk density and their composition can be inferred (Fortney et al. 2007; Seager et al. 2007; Piaulet et al. 2022), with which their mass accretion history and post-formation evolution can be constrained. The core accretion model successfully predicted the so-called "planet dessert" (Ida & Lin 2004), which refers to the paucity of planets with tens of Earth-mass within 3 au. The "radius valley" (Fulton et al. 2017) is manifested at a planet radius GLYPH<24> 2 𝑅 GLYPH<8> , which has been attributed to photoevaporation-driven mass loss (Owen & Wu 2013, 2017), planet formation location with respect to snow line (Luque & Pallé 2022; Izidoro et al. 2022), or core-powered mass loss (Ginzburg et al. 2018).  \nThe physical principles underlying migration of low-mass planets in gaseous discs (Type I migration) have long been established (Goldreich & Tremaine 1979; Lin & Papaloizou 1979). The total torque exerted on planets by the surrounding disc is typically negative, resulting in inward planet migration on time-scales shorter than the disc lifetimes (e.g., Ward 1997; Tanaka et al. 2002; Ribas et al. 2014; Winter et al. 2019). However, the direction of migration can be  \n- ★ E-mail: huangs20@mails.tsinghua.edu.cn\n- y E-mail: chrisormel@tsinghua.edu.cn  \nreversed at special locations where conditions materialize that render a net positive torque, resulting in migration traps. These locations include the region where the horseshoe saturates, i.e when the (positive) co-rotation torque compensates the (negative) Lindblad torque (Goldreich & Tremaine 1980; Ward 1991; Paardekooper et al. 2010, 2011), the disc inner edge (Liu et al. 2017; Romanova et al. 2019; Ataiee & Kley 2021) where the torque becomes one-sided, and the regions where the disc switches from optically thin to optically thick (0 GLYPH<149> 1 GLYPH<0> 1 au) (Masset et al. 2006). In addition, in the pebble accretion paradigm, the infalling dust can efficiently induce (positive) thermal torque onto the planets (Benítez-Llambay et al. 2015; Masset 2017; Guilera et al. 2019, 2021).  \nAlthough the present close-in positions of exoplanets indirectly hint at planet migration, it is hard to quantitatively test the theory based only on single-planet systems. Instead, multi-planet systems, especially those with planet pairs in Mean Motion Resonance (MMR) leave richer dynamical imprints against which the theory can be tested (Snellgrove et al. 2001; Papaloizou & Szuszkiewicz 2005). It is likely that such resonant architecture results from migration in a gas-rich environment (i.e., the disc) as energy dissipation is needed to trap planets in resonance (Terquem & Papaloizou 2007; Raymond et al. 2008; Rein 2012; Batygin 2015). One famous example is PDS 70, which harbours two directly imaged near-resonance planets in its protoplanetary disc (Bae et al. 2019; Benisty et al. 2021). Besides, a chain of planets in resonance might sculpt the asymmetry disc structure in HD 163296 (Isella et al. 2018; Garrido-Deutelmoser et al. 2023). In particular, there tends to be an excess of systems with planets\' period ratio just wide of commensurability (Fabrycky et al.  \n2014; Steffen & Hwang 2015), indicating that a certain fraction of planet pairs are truly in resonance.  \nFormally, two planets are said to be in \' 𝑗 , 1 \' : 𝑗 resonance if at least one of their resonance angles ( 𝜙 1 GLYPH<150> 2 = \' 𝑗 , 1 \' 𝜆 2 GLYPH<0> 𝑗𝜆 1 GLYPH<0> 𝜛 1 GLYPH<150> 2 , with 𝜆 𝑖 the mean longitudes and 𝜛 𝑖 the longitude of pericentres) librate around a fixed value. However, the values of the resonance angles are poorly constrained because it is hard to constrain 𝜛 𝑖 for near-circular orbits. We therefore turn our attention to their period ratios and define the dimensionless parameter  \nΔ = 𝑃 2 𝑃 1 GLYPH<0> 𝑗 , 1 𝑗 (1)  \nto measure the offset of the period ratio away from a first-order \' 𝑗 , 1 \' : 𝑗 commensurability. Here, 𝑃 1 is the period of the inner planet and 𝑃 2 that of the outer. If two planets are in a \' 𝑗 , 1 \' : 𝑗 resonance, the offset Δ must be close to zero. Xie (2014) and Ramos et al. (2017) emphasize that the exact value of Δ is determined by migration and eccentricity damping, linking the observed quantity Δ to planet migration (Charalambous et al. 2022). It offers us an opportunity to examine planet-disc interaction histories through planets in MMR in multi-planet systems. The migration history of such specific multi-planet systems like TRAPPIST-1 (Gillon et al. 2017; Luger et al. 2017; Huang & Ormel 2022), K2-24 (Petigura et al. 2018; Teyssandier & Libert 2020) and TOI-1136 (Dai et al. 2022), can be reconstructed.  \nYet, most exoplanets are obviously not in resonance because of their large offsets Δ . Various scenarios have been proposed to explain the overall observed non-resonant planetary architecture statistically. These include dynamical instability (Izidoro et al. 2017, 2021), disc winds (Ogihara et al. 2018), in situ formation of sub-Netunes (Dawson et al. 2015; Choksi & Chiang 2020), planetesimal scattering (Chatterjee & Ford 2015; Ghosh & Chatterjee 2022), stellar tides (Lithwick & Wu 2012; Delisle & Laskar 2014; Xie 2014; Sánchez et al. 2020), and stochastic forces (Rein & Papaloizou 2009; Goldberg &Batygin 2022). However, many of these works introduce additional free parameters in order to match the only observed quantity - the period ratio distribution. Overfitting may occur. One way to improve on this is to include more observational quantities. For example, Goldberg & Batygin (2022) and Choksi & Chiang (2022) find that by additionally accounting for the TTV signatures of hundreds of Kepler planets, a laminar disc alone cannot reproduce the observed TTV features. Either additional planets (perturbers) are required (Choksi &Chiang 2022) or their birth proto-discs are turbulent (Goldberg & Batygin 2022).  \nIn this work, we introduce a statistical approach to study all transiting exoplanets and try to answer when, where, and how many planets are captured in resonance. The observed planetary radii, orbital periods, and host stellar masses are taken into account. By assuming that the period ratio distribution is characterized by two distributions, representing both resonant and non-resonant planets, the probability that a planets pair is in resonance is evaluated. Whether the migration and eccentricity damping is consistent with migration theory is determined. In addition, our approach allows us to constrain the timing and the location, of the resonance trapping, which implies the pathway of planet formation.  \nThe paper is structured as follows: We first introduce our statistical model in Sect. 2. Using the MCMC method, we constrain the relation between eccentricity damping time-scale and migration time-scale in Sect. 4. Along with the resonance trapping criterion that we improve on (Sect. 3), we address when, where, and how many planets are in resonance Sect. 5. The discussion and conclusion of this study are presented in Sect. 6 and Sect. 7, respectively.', '2 METHODOLOGY': 'In this section, we describe our disc model and migration (Type I), and construct the likelihood function needed for the Monte Carlo Markov Chain (MCMC) simulations of Sect. 4. Since our migration model is linear with planet mass (Type I), we focus on the planets with relatively low planet-to-star mass ratios which are not likely to open a gap when they are in the protoplanet disc. We assume that the relevant disc quantities follow power-law distribution (Sect. 2.1). We describe the migration model in Sect. 2.2. The equilibrium dynamics of planets trapping in resonance are described in Sect. 2.3. The masses of observed planets are calculated from a mass-radius (M-R) relationship. Its prescription is given in Sect. 2.4. The total log-likelihood function we construct for the MCMC is detailed in Sect. 2.5.', '2.1 Disc model': "As most observed transiting planets are located within GLYPH<24> 1 au of their host star, we will describe the possible structures of the inner disc. We assume that the gas surface density always follows a power-law distribution:  \nΣ ' 𝑟 ' = Σ 1au GLYPH<16> 𝑟 1 au GLYPH<17> 𝑠 GLYPH<150> (2)  \nwhere Σ 1au is the gas surface density at 1 au and 𝑠 is its slope. The gas aspect ratio also follows a power-law distribution:  \nℎ ' 𝑟 ' = ℎ 1au GLYPH<16> 𝑟 1 au GLYPH<17> 𝑞 (3)  \nwhere ℎ 1au is the gas aspect ratio at 1 au and 𝑞 is its slope. Different assumptions about the disc structure, e.g., heating mechanisms, result in distinct values of 𝑠 , 𝑞 , Σ 1au , and ℎ 1au . Typically, the inner disc is optically thick and the main heating energy comes from viscous dissipation (Ruden & Lin 1986), while the outer disc is optically thin and stellar irradiation mainly heat the disc onto its surface layer (Chiang & Goldreich 1997).  \nFor discs dominated by stellar radiation, we make use of the disc structure from Liu et al. (2019). The gas surface density is:  \nΣ g GLYPH<150> irr = 250 / 𝑀 𝑔 10 GLYPH<0> 8 𝑀 GLYPH<12> GLYPH<157> yr ! GLYPH<18> 𝑀 ★ 𝑀 GLYPH<12> GLYPH<19> 9 14 GLYPH<18> 𝐿 ★ 𝐿 GLYPH<12> GLYPH<19> GLYPH<0> 2 7 GLYPH<16> 𝑟 au GLYPH<17> GLYPH<0> 15 14 g cm GLYPH<0> 2 (4)  \nwhere / 𝑀 𝑔 is stellar accretion rate and 𝐿 ★ is the star luminosity. The aspect ratio is:  \nℎ g GLYPH<150> irr = 0 GLYPH<149> 0245 GLYPH<18> 𝑀 ★ 1 𝑀 GLYPH<12> GLYPH<19> GLYPH<0> 4 GLYPH<157> 7 GLYPH<18> 𝐿 ★ 1 𝐿 GLYPH<12> GLYPH<19> 1 GLYPH<157> 7 GLYPH<16> 𝑟 1 au GLYPH<17> 2 GLYPH<157> 7 GLYPH<149> (5)  \nFor stars of mass between 0 GLYPH<149> 43 𝑀 GLYPH<12> and 2 𝑀 GLYPH<12> , which covers most of our star sample, the mass-luminosity relation is well represented by 𝐿 ★ GLYPH<157> 𝐿 GLYPH<12> = ' 𝑀 ★ GLYPH<157> 𝑀 GLYPH<12> ' 4 (Duric 2004). Therefore, the disc aspect ratio simplifies to:  \nℎ g GLYPH<150> irr = 0 GLYPH<149> 0245 GLYPH<16> 𝑟 1 au GLYPH<17> 2 GLYPH<157> 7 (6)  \nindependent of stellar mass.  \nIf the inner disc is dominated by viscous heating, its temperature structure is highly related to the viscous accretion rate and opacity. Following Liu et al. (2019), 𝑠 and 𝑞 are taken to be GLYPH<0> 3 GLYPH<157> 8 and GLYPH<0> 1 GLYPH<157> 16, while Σ 1au and ℎ 1au are not specified.", '2.2 Type I migration': "In the Type I migration regime, planet migration is the result of a net torque Γ net consisting of the Lindblad (Ward 1986, 1997), corotation (Goldreich & Tremaine 1979; Ward 1992) and thermal torques (Benítez-Llambay et al. 2015; Masset 2017; Guilera et al. 2019, 2021), etc. Usually, the net torque is negative and the planet migrates inward. In the Type I limit the migration speed is proportional to disc mass and planet mass. In the limit of a locally isothermal disc, which implies that temperature is a function of radius only, 𝑇 ' 𝑟 ' , the type I migration time-scale for the 𝑖 -th planet at distance 𝑟 𝑖 is:  \n𝜏 𝑎 𝑖 = 𝐿 𝑖 Γ net = 𝛾 𝐼 𝜏 𝑤 𝑖 ℎ ' 𝑟 𝑖 ' 2 GLYPH<150> (7)  \nwhere 𝐿 𝑖 is the angular momentum of the planet, 𝛾 𝐼 = 2 GLYPH<149> 7 GLYPH<0> 1 GLYPH<149> 1 𝑠 (Tanaka et al. 2002) is the Type I migration prefactor (D'Angelo & Lubow 2010), and ℎ ' 𝑟 𝑖 ' is the disc gas aspect ratio at 𝑟 𝑖 . The characteristic time of the orbital evolution (Tanaka & Ward 2004) is:  \n𝜏 𝑤 𝑖 = 1 𝜇 𝑖 𝑀 ★ Σ ' 𝑟 𝑖 ' 𝑟 2 𝑖 ℎ ' 𝑟 𝑖 ' 4 Ω 𝐾 ' 𝑟 𝑖 ' GLYPH<150> (8)  \nwhere 𝜇 𝑖 is the mass ratio of the 𝑖 -th planet over its host star and Ω 𝐾 ' 𝑟 𝑖 ' = √︃ 𝐺𝑀 ★ GLYPH<157> 𝑟 3 𝑖 is the Keplerian angular velocity at distance 𝑟 𝑖 . The eccentricity damping rate is proportional to the local surface density and planet mass. It is given by:  \n𝜏 𝑒 𝑖 = 𝐶 𝑒 𝜏 𝑤 𝑖 0 GLYPH<149> 78 = 𝐶 𝑒 0 GLYPH<149> 78 𝛾 𝐼 ℎ 2 𝜏 𝑎 GLYPH<150> (9)  \nwhere 𝐶 𝑒 stands for eccentricity damping efficiency. Although Cresswell & Nelson (2008) gives 𝐶 𝑒 GLYPH<25> 1, lower values are needed in other studies to reproduce specific systems. TRAPPIST-1 planets demand 𝐶 𝑒 GLYPH<25> 0 GLYPH<149> 1 (Huang & Ormel 2022) and K2-24 requires 𝐶 𝑒 GLYPH<25> 0 GLYPH<149> 28 (Teyssandier & Libert 2020). A recently discovered GLYPH<24> 100 Myr old exoplanet system TOI-1136, on the other hand, suggest 𝐶 𝑒 GLYPH<24> 10 (Dai et al. 2022). We are therefore agnostic about the value of 𝐶 𝑒 , which value we intend to constrain through our MCMC fitting.", '2.3 Dynamics of resonance trapping': "Resonance trapping is a natural outcome of convergent disc migration and eccentricity damping, especially for first-order resonances ' 𝑗 , 1 ' : 𝑗 . When two planets are in first-order resonance, their eccentricities and period ratio librate around their equilibrium values. Suchequilibrium has been studied by Goldreich & Schlichting (2014) and Terquem & Papaloizou (2019). They both give the equilibrium eccentricity (the eccentricity where tidal damping equals resonant excitation) of the inner planet:  \n𝑒 2 1 GLYPH<150> eq = 𝜏 𝑒 1 GLYPH<157> 𝜏 𝑎 2 GLYPH<0> 𝜏 𝑒 1 GLYPH<157> 𝜏 𝑎 1 2 ' 𝑗 , 1 ' GLYPH<16> 1 , 𝑗 𝑗 , 1 𝜇 1 𝛼𝜇 2 GLYPH<17> GLYPH<20> 1 , 𝜇 1 𝛼𝜇 2 GLYPH<16> 𝑗 𝑗 , 1 GLYPH<17> 2 GLYPH<16> 𝑓 0 2 𝑓 1 GLYPH<17> 2 𝜏 𝑒 1 𝜏 𝑒 2 GLYPH<21> (10)  \nwhere 𝛼 is ratio of the inner-to-outer semi-major axis . Here, and in the following, the subscription '1' stands for the inner planet and '2' for the outer planet. The relationship between the eccentricities of the inner and outer planets is:  \n𝑒 2 2 𝑒 2 1 = GLYPH<18> 𝜇 1 𝛼𝜇 2 𝑗 𝑗 , 1 𝑓 0 2 𝑓 1 GLYPH<19> 2 GLYPH<150> (11)  \nwhere 𝑓 1 and 𝑓 0 2 are coefficients tabulated in Terquem & Papaloizou (2019, their Table A1). The equilibrium value for the offset from  \nexact resonance is:  \nΔ eq = GLYPH<0> 𝑓 0 2 𝜇 1 1 𝑗 1 𝑒 2 GLYPH<150> eq GLYPH<149> (12)  \nAt this distance, the resonance repulsion equals the Type-I inward migration.  \nIn order to calculate the period ratios for planets in resonance, we need to know what their equilibrium eccentricities are. This calculation can be done only after the values of 𝜏 𝑒 1 GLYPH<157> 𝜏 𝑒 2 , 𝜏 𝑒 1 GLYPH<157> 𝜏 𝑎 2 and 𝜏 𝑒 1 GLYPH<157> 𝜏 𝑎 1 are known. Combining the disc model and migration model, we have  \n𝜏 𝑒 1 𝜏 𝑒 2 = 𝜇 2 𝜇 1 GLYPH<18> 𝑟 1 𝑟 2 GLYPH<19> 4 𝑞 GLYPH<0> 𝑠 GLYPH<0> 0 GLYPH<149> 5 (13)  \nand  \n𝜏 𝑒 1 𝜏 𝑎 2 = 𝜏 𝑒 2 𝜏 𝑎 2 𝜏 𝑒 1 𝜏 𝑒 2 = 𝐶 𝑒 ℎ ' 𝑟 2 ' 2 0 GLYPH<149> 78 𝛾 𝐼 𝜇 2 𝜇 1 GLYPH<18> 𝑟 1 𝑟 2 GLYPH<19> 4 𝑞 GLYPH<0> 𝑠 GLYPH<0> 0 GLYPH<149> 5 (14)  \nwhere 𝑠 and 𝑞 are gas surface density and aspect ratio gradient.  \nFor the eccentricity-to-semi-major axis damping of the inner planet, we distinguish it between two cases:  \n- (i) Migrating pair. If resonances are formed during migration and the ambient disc disperses before the planet pairs reach a migration barrier, planet migration and eccentricity damping follow Eq. (7) and Eq. (9). Therefore,  \n𝜏 𝑒 1 𝜏 𝑎 1 = 𝐶 𝑒 ℎ ' 𝑟 1 ' 2 0 GLYPH<149> 78 𝛾 𝐼 (15)  \nIt requires that the outer planet migrates faster than the inner planet to guarantee convergent migration.  \n(ii) Braking pair. On the other hand, if resonances are formed before/after the inner planet's migration is halted by a barrier (could be disc inner edge or the radius where reverse migration occurs), there is no net torque on the two planets. Angular momentum conservation gives 𝜏 𝑎 1 GLYPH<157> 𝜏 𝑎 2 = GLYPH<0> √︁ 𝜇 1 𝑟 1 GLYPH<157> 𝜇 2 𝑟 2 . Hence,  \n𝜏 𝑒 1 𝜏 𝑎 1 = GLYPH<0> 𝐶 𝑒 ℎ ' 𝑟 2 ' 2 0 GLYPH<149> 78 𝛾 𝐼 GLYPH<18> 𝜇 2 𝜇 1 GLYPH<19> 3 2 GLYPH<18> 𝑟 1 𝑟 2 GLYPH<19> 4 𝑞 GLYPH<0> 𝑠 GLYPH<0> 1 = GLYPH<0> 𝐶 𝑒 ℎ ' 𝑟 1 ' 2 0 GLYPH<149> 78 𝛾 𝐼 GLYPH<18> 𝜇 2 𝜇 1 GLYPH<19> 3 2 GLYPH<18> 𝑟 1 𝑟 2 GLYPH<19> 2 𝑞 GLYPH<0> 𝑠 GLYPH<0> 1 (16)  \nFor the outer planet 𝜏 𝑒 2 GLYPH<157> 𝜏 𝑎 2 = 𝐶 𝑒 ℎ ' 𝑟 2 ' 2 GLYPH<157> 0 GLYPH<149> 78 𝛾 𝐼 is always true. It is also apparent that the resonant equilibrium does not depend on disc mass ( Σ 1au ) but on power law indices and aspect ratio ( 𝑠 , 𝑞 , and ℎ 1au ).", '2.4 Mass-radius relations': 'Planets\' masses are also needed to calculate the equilibrium period ratios in resonance. However, most transiting planets have poorly constrained masses compared to their radii.  \nIf the planet mass is not yet constrained from e.g., Transit Timing Variation (TTV Agol et al. 2005) or Radial Velocity (RV), we then obtain the planet mass using a mass-radius relation. The planet sample for fitting the mass-radius relation is based on the data of transiting planets from the NASA Exoplanet Archive 1 . Planets with masses lower than 200 𝑀 GLYPH<8> , radii smaller than 20 𝑅 GLYPH<8> , or periods longer than 5 days (see Sect. 4.1) are selected while planets with periods shorter than 5 days are excluded.  \nFigure 1. Mass-radius (M-R) relationship for exoplanets. Planets are selected based on the data of transiting planets in the NASA Exoplanet Archive. Only planets with mass lower than 200 𝑀 GLYPH<8> , radii smaller than 20 𝑅 GLYPH<8> and periods longer than 5 days are included in the sample. Planets\' masses and radii are indicated by black dots, with 1 𝜎 error bars. The green colour indicates the fit result from Ramos et al. (2017), orange colour is the updated fit result by this work (Eq. (17)). The values within their 1 𝜎 dispersion ( 𝜎 m = 0 GLYPH<149> 374 for our fit) are contained within the light green and orange region respectively.  \n<!-- image -->  \nOur fitting approach is identical to Ramos et al. (2017). They use a broken power law expression, which fits two different power law relations for larger bodies ( 𝑅 ¡ 𝑅 crit ) and smaller bodies ( 𝑅 GLYPH<159> 𝑅 crit ).  \nlog 10 𝑀 𝑖 𝑀 GLYPH<8> ! = 8 > > < > > 𝑎 , 𝑏 log 10 GLYPH<16> 𝑅 𝑖 𝑅 GLYPH<8> GLYPH<17> if 𝑅 𝑖 GLYPH<20> 𝑅 crit 𝑐 , 𝑑 log 10 GLYPH<16> 𝑅 𝑖 𝑅 GLYPH<8> GLYPH<17> if 𝑅 𝑖 ¡ 𝑅 crit (17)  \n:  \nwhere 𝑀 𝑖 and 𝑅 𝑖 is the mass and radius of the 𝑖 -th planet. We make use of the Maximum Likelihood Estimator (MLE) to maximize a Gaussian likelihood centered at Eq. (17). The critical radius 𝑅 crit is also fitted. It estimates that: 𝑎 = 0 GLYPH<149> 69, 𝑏 = 0 GLYPH<149> 78, 𝑐 = 0 GLYPH<149> 11, 𝑑 = 1 GLYPH<149> 7, 𝑅 crit = 4 GLYPH<149> 23 𝑅 GLYPH<8> and the corresponding dispersion 𝜎 𝑚 = 0 GLYPH<149> 374. The best-fitting relation is shown in Fig. 1 in comparison with the relation fitted by Ramos et al. (2017). They do not differ significantly. The estimated value for 𝑅 crit is consistent with Teske et al. (2021) who suggest a single power law relation for planets with 𝑅 GLYPH<159> 3 GLYPH<149> 25 𝑅 GLYPH<8> .  \nEquation (17) allows us to calculate the planet mass and therefore the resonance offset (Eq. (12)). Moreover, the log-normal dispersion in planet mass enables us to check our model consistency. The reason is that any uncertainty in the mass, will propagate, through Eq. (10) and Eq. (12). Therefore, we expect that the obtained value for 𝜎 Δ from MCMC fitting is similar to, or exceeds, 𝜎 𝑚 .  \nThe log-normal dispersion in planet mass significantly simplifies our analysis. There are arguably more sophisticated forms of M-R relationships, e.g., the one given by Wolfgang et al. (2016) and improved by Teske et al. (2021). However, they assume planet mass follows a normal dispersion instead of a log-normal. In that case, the resulting Δ would follow a complicated form of the ratio distribution 2 .  \n2 If variables 𝑋 and 𝑌 follow a dependent (independent) normal distribution with nonzero mean values, the new variable 𝑍 = 𝑋 GLYPH<157> 𝑌 follows correlated (uncorrelated) non-central normal ratio distribution (Hinkley 1969; Hayya et al. 1975).  \nFigure 2. Classification of resonant and non-resonant planet pairs. In the pie chart, \'migrating\' means the pairs lock into resonance before reaching a migration barrier, whereas "at disc edge" indicates planets lock into resonance at a planet migration trap. "Inner slow" ("fast") means that the inner planet has a longer (shorter) migration time-scale than its outer planet. 𝑁 res is the number of resonant pairs. 𝑁 1, 𝑁 2, and 𝑁 5 are the numbers of Migrating , Braking , and non- resonant pairs in the pairs with the inner planets migrating slower than the outer. 𝑁 3 and 𝑁 4 are the numbers of Braking and non- resonant pairs among the pairs with inner planets migrating faster than the outer.  \n<!-- image -->', '2.5 A statistical model of resonant and non-resonant planets': 'We define the posterior distribution as: 𝑝 \' 𝜽 j Δ obs GLYPH<150> 𝑘 GLYPH<150> 𝑿 𝑘 \' = 𝑝 \' Δ obs GLYPH<150> 𝑘 j 𝜽 GLYPH<150> 𝑿 𝑘 \' 𝑝 \' 𝜽 j 𝑿 obs GLYPH<150> 𝑘 \' GLYPH<149> If the disc structure is not specified (without knowing the specific values of 𝑠 , 𝑞 in Eq. (2) and Eq. (3)), the unknown model parameters are 𝜽 = \' log 10 \' 𝐶 𝑒 ℎ 2 1au \' GLYPH<150> 𝜎 Δ GLYPH<150> 𝑠GLYPH<150> 𝑞 \' and the known parameters 𝑿 𝑘 = \' 𝑀 ★ GLYPH<150> 𝑀 1 GLYPH<150> 𝑀 2 GLYPH<150> 𝑟 1 GLYPH<150> 𝑟 2 GLYPH<150> 𝑗 \' . The index 𝑘 indicates the 𝑘 -th planet pair. If a disc structure is specified, the unknown model parameters are 𝜽 = \' log 10 𝐶 𝑒 GLYPH<150> 𝜎 Δ \' for the irradiation disc and 𝜽 = \' log 10 \' 𝐶 𝑒 ℎ 2 1au \' GLYPH<150> 𝜎 Δ \' for the viscous disc, and the known parameters are 𝑿 𝑘 = \' 𝑀 ★ GLYPH<150> 𝑀 1 GLYPH<150> 𝑀 2 GLYPH<150> 𝑟 1 GLYPH<150> 𝑟 2 GLYPH<150> 𝑗 GLYPH<150> ℎ 1au GLYPH<150> 𝑠GLYPH<150> 𝑞 \' . We assume that the prior 𝑝 \' 𝜽 j 𝑿 𝑘 \' follows a uniform distribution (Table 1). The resonance index 𝑗 is regarded as one of the known parameters such that we can analyse all pairs at once.  \nAs mentioned, resonance trapping naturally results from convergent migration: either two planets get trapped into resonance during migration (both inward) with the outer planet migrating faster than the inner one, or the inner planet reaches a migration barrier with the outer planet arriving at a later time. Following the discussion in Sect. 2.4, we assume that the period ratio of a planet pair that is in resonance obeys a log-normal distribution log 10 Δ GLYPH<24> N\' log 10 Δ m GLYPH<157> s \' 𝜽 GLYPH<150> 𝑿 𝑘 \' GLYPH<150> 𝜎 2 Δ \' :  \n𝑝 res GLYPH<150> m GLYPH<157> s \' Δ j 𝜽 GLYPH<150> 𝑿 𝑘 \' 𝑑 log 10 Δ = 1 √︃ 2 𝜋𝜎 2 Δ exp " GLYPH<0> » log 10 Δ GLYPH<0> log 10 Δ m GLYPH<157> s \' 𝜽 GLYPH<150> 𝑿 𝑘 \'… 2 2 𝜎 2 Δ # 𝑑 log 10 Δ GLYPH<150> (18)  \nwhere Δ m GLYPH<157> s indicates the resonance commensurability calculated by Eq. (12). The resonance offset Δ m indicates the value calculated for Migrating pairs using Eq. (15), while Δ s is for Braking pairs and is calculated using Eq. (16).  \nIn the Type I migration regime, if the inner planet is more massive (migrates faster), convergent migration can only occur when the inner planet has reached the migration barrier. On the other hand, if the inner planet\'s migration is slower than the outer, migration is always convergent. In that case, trapping can occur when both planets migrate inward or when the inner planet has reached the migration barrier. For ease of the likelihood calculation, we further divide the resonant planet pairs into three categories:  \nTable 1. Prior bounds or values we take for the disc parameters of the three different models. All priors follow a uniform distribution ( U ).  \nGroup 1: The inner planet migrated slower ( 𝜏 𝑎 1 ¡ 𝜏 𝑎 2 ) and it did not reach a migration barrier ( 𝑁 1 );  \nGroup 2: The inner planet migrated slower ( 𝜏 𝑎 1 ¡ 𝜏 𝑎 2 ) and it reached a migration barrier ( 𝑁 2 );  \nGroup 3: The inner planets migrated faster ( 𝜏 𝑎 1 GLYPH<20> 𝜏 𝑎 2 ) and it reached a migration barrier ( 𝑁 3 ).  \n𝑁 1 , 𝑁 2 and 𝑁 3 represent the number of resonant planet pairs corresponding to each type of resonance. In addition, 𝑁 4 ( 𝑁 5 ) represents the number of pairs that are not in resonance with the inner planet migrating faster (slower) than the outer planet. We provide a sketch to explain the five classes in Fig. 2.  \nIf the inner planets migrate faster ( 𝜏 𝑎 1 GLYPH<20> 𝜏 𝑎 2 ), those pairs in resonance must be Braking pairs. However, if the inner planets migrate slower, they could either be Migrating pairs or Braking pairs. The log 10 Δ distribution of planets in resonance is therefore (hereafter, 𝑝 res \' Δ \' represents 𝑝 res \' Δ obs j 𝜽 GLYPH<150> 𝑿 𝑘 \' , etc):  \n𝑝 res \' Δ \' = GLYPH<26> 𝑝 res GLYPH<150> s \' Δ \' if 𝜏 𝑎 1 GLYPH<20> 𝜏 𝑎 2 GLYPH<150> GLYPH<2> 𝑝 res GLYPH<150> s \' Δ \' , 𝑝 res GLYPH<150> m \' Δ \' GLYPH<3> if 𝜏 𝑎 1 ¡ 𝜏 𝑎 2 GLYPH<149> (19)  \nTheperiod ratio of planet pairs that are not in resonance is assumed to follow a uniform distribution:  \n𝑝 n GLYPH<0> res \' Δ \' 𝑑 log 10 Δ = \' ln 10 \' Δ Δ max 𝑑 log 10 Δ (20)  \nwhere Δ max is the value above which a planet pair with Δ obs is not considered in resonance. We use Δ max = \' 3 𝑗 , 2 \'GLYPH<157> \' 3 𝑗 GLYPH<0> 1 \' GLYPH<0> \' 𝑗 , 1 \'GLYPH<157> 𝑗 which is the distance from the first order \' 𝑗 , 1 \' : 𝑗 resonance to its closest external 3rd order \' 3 𝑗 , 2 \' : \' 3 𝑗 GLYPH<0> 1 \' resonance.  \nFinally, the total log-likelihood is written as:  \nln L = 𝑁 pairs ∑︁ 𝑘 = 1 ln GLYPH<2> 𝑝 res \' Δ obs GLYPH<150> 𝑘 \' , 𝑝 n GLYPH<0> res \' Δ obs GLYPH<150> 𝑘 \' GLYPH<3> GLYPH<150> (21)  \nwhere 𝑁 pairs is the number of planet pairs in our sample.', '3 RESONANCE TRAPPING CRITERION FOR THE RESTRICTED 3-BODY PROBLEM': 'In this section, we improve and numerically verify the two-body resonance trapping criterion. This new trapping condition will be used in Sect. 5.2 to further constrain the statistical results of Sect. 4.2. 3  \n3 Batygin & Petit (2023) have recently presented an analysis with a trapping condition also predicated on the equilibrium phase angle, Eq. (25), like in this Section. Their findings are consistent with ours.', '3.1 Theoretical Derivation': "In the restricted three-body problem, the outer planet is on a fixed circular orbit. The inner planet moves outward and its semi-major axis and eccentricity are damped on time-scales of 𝜏 𝑎 and 𝜏 𝑒 , respectively. Lagrange's planetary equation for the mean motion ( 𝑛 ) then reads:  \n/ 𝑛 1 = GLYPH<0> 3 𝑗 𝛼 𝑓 1 𝜇 2 𝑒 1 𝑛 2 1 sin 𝜙 1 GLYPH<0> 3 𝑛 1 2 𝜏 𝑎 , 𝑝𝑒 2 1 𝑛 1 𝜏 𝑒 GLYPH<150> (22)  \nwhere 𝛼 = 𝑎 1 GLYPH<157> 𝑎 2 , 𝜙 1 , 𝑒 1 and 𝑓 1 are the semi-major axis ratio, resonance angle, the eccentricity of the inner planet, and 𝑓 1 is a numerical factor that depends on the resonance index 𝑗 (Murray & Dermott 1999; Terquem & Papaloizou 2019). By definition, 𝑝 = 3 holds when the eccentricity damping operates at constant angular momentum (Teyssandier & Terquem 2014). Lagrange's planetary equation for eccentricity is:  \n/ 𝑒 1 = GLYPH<0> 𝛼 𝑓 1 𝜇 2 𝑛 1 sin 𝜙 1 GLYPH<0> 𝑒 1 𝜏 𝑒 1 GLYPH<150> (23)  \nWhen two planets are in resonance, the values of different orbital properties e.g., 𝑒 1 , 𝛼 and 𝜙 1 librate around their equilibrium values. The equilibrium eccentricity (Goldreich & Schlichting 2014; Terquem & Papaloizou 2019) is derived by putting / 𝑒 = / 𝑎 = 0 and eliminating sin 𝜙 1 in Eq. (22) and Eq. (23):  \n𝑒 1 GLYPH<150> eq = √︂ 𝜏 𝑒 1 2 ' 𝑗 , 1 ' 𝜏 𝑎 1 GLYPH<149> (24)  \nBy inserting the equilibrium eccentricity 𝑒 1 GLYPH<150> eq and 𝛼 eq GLYPH<25> » 𝑗 GLYPH<157> ' 𝑗 , 1 '… 2 GLYPH<157> 3 into Eq. (23), sin 𝜙 1 follows:  \nsin 𝜙 1 GLYPH<150> eq = GLYPH<0> 1 𝛼 𝑓 1 𝜇 2 𝑛 1 √︄ 1 2 ' 𝑗 , 1 ' 𝜏 𝑒 1 𝜏 𝑎 1 (25)  \nNaturally, its absolute value cannot exceed 1. Otherwise, for j sin 𝜙 1 j ¡ 1, no steady state exists and the planets will cross the resonance. Combining, / 𝑒 1 = 0, Eq. (23) and Eq. (24) we can write the resonance trapping condition:  \n𝜏 𝑎 1 𝜏 𝑒 1 GLYPH<21> 1 2 ' 𝑗 , 1 '' 𝛼 𝑓 1 𝜇 2 𝑛 1 ' 2 GLYPH<149> (26)  \nThe classical theory about the resonance trapping criterion is that the time for the planet to migrate across the libration width is shorter than the libration time-scale (Ogihara & Kobayashi 2013; Batygin 2015). For comparison, we also provide the criterion derived from the classical pendulum model (Murray & Dermott 1999; Ogihara & Kobayashi 2013; Huang & Ormel 2022):  \n𝜏 𝑎 1 𝜏 𝑒 1 GLYPH<21> 𝜋 ' 𝑗 , 1 ' 4 ' 𝛼 𝑓 1 𝜇 2 𝑛 1 ' 2 GLYPH<149> (27)  \nCompared to Eq. (27), the new criterion (Eq. (26)) has the same dependence on planet-to-star mass ratio 𝜇 2 and the orbital frequency 𝑛 1 but differs regarding the resonance index 𝑗 .", '3.2 Comparison with simulation': "We compare the new resonance trapping criterion above against the numerical simulation. The fiducial accelerations accounting for migration and eccentricity of planets in the simulations are expressed by:  \n𝒂 m = GLYPH<0> 𝒗 2 𝜏 𝑎 GLYPH<150> (28)  \n𝒂 e = GLYPH<0> 2 ' 𝒗 GLYPH<1> 𝒓 ' 𝒓 2 𝑟 2 𝜏 𝑒 (29)  \nFigure 3. Resonance trapping and crossing for the 2:1 resonance. For each point in the parameter space, we run five simulations with different initial conditions to obtain the probabilistic result of resonance trapping (green shading). The blue line is the resonance trapping criterion derived from the pendulum model (Eq. (27)) while the orange line is our newly derived criterion (Eq. (26)). The grey solid lines are the contours of sin 𝜙 1, with corresponding values labelled in white. The red dashed line (Eq. (30)) corresponds to the transition from stable resonance trapping (above) to overstable resonance (below). Below this line, the simulations progressively evolve into the overstable territory, see https://raw.githubusercontent.com/ shuohuangGIT/Infer-migration-history/main/q1.mp4 .  \n<!-- image -->  \n(Papaloizou & Larwood 2000; Cresswell & Nelson 2006, 2008). We make use of the WHfast integrator of the open-source N-body code REBOUND (Rein 2012). The migration and eccentricity damping on planets are implemented through REBOUNDx (Tamayo et al. 2020).  \nWe fix the outer planet on a circular orbit at »' 𝑗 , 1 'GLYPH<157> 𝑗 … 2 GLYPH<157> 3 au. The inner planet starts to migrate outward at 0 GLYPH<149> 8 au on a time-scale of 𝜏 𝑎 . Its eccentricity is damped on a time-scale of 𝜏 𝑒 . The planet mass is fixed at 10 𝑀 GLYPH<8> and the host mass is 1 𝑀 GLYPH<12> . We vary two parameters in the simulation: 𝜏 𝑎 from 10 4 yr (fast migration) to 10 7 yr (slow migration), and 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 ranging from 10 1 (inefficient eccentricity damping) to 10 4 (efficient eccentricity damping). Each parameter is sampled by 100 grid points evenly distributed in logspace. In order to capture the probabilistic behavior of resonance trapping, we run five simulations for each point in the 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 -𝜏 𝑎 parameter space, where we evenly sample the initial longitude of the inner planet.  \nWe conduct simulations for 𝑗 = 1 (2:1 resonance) and 𝑗 = 2 (3:2 resonance). We run the simulation until 𝑡 = 𝜏 𝑎 , but we take a snapshot at 𝑡 = 0 GLYPH<149> 2 𝜏 𝑎 . If the period ratio 𝑃 2 GLYPH<157> 𝑃 1 decreases below ' 𝑗 , 1 'GLYPH<157> 𝑗 , we classify the simulation outcome as a resonance crossing. The results are shown in Fig. 3 and Fig. 4 for 𝑗 = 1 and 𝑗 = 2 respectively. Resonance crossing cases are in white and resonance trapping cases are in green. The red dashed line indicates the boundary below which the trapping solution becomes overstable,  \nGLYPH<18> 𝜏 𝑎 𝜏 𝑒 GLYPH<19> overstable = 1 8 ' 𝑗 , 1 ' GLYPH<18> GLYPH<0> 3 𝑗 2 𝛼 𝑓 1 𝜇 2 GLYPH<19> 2 GLYPH<157> 3 (30)  \n(Goldreich & Schlichting 2014), which evaluates to 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 GLYPH<25> 170 for 𝜇 2 = 3 GLYPH<2> 10 GLYPH<0> 5 , both for 𝑗 = 1 and 2. Above this line, all systems are either trapped in resonance or not. The top-right corner indicates the parameter space where the two planets both get captured and permanently stay in resonance and the top-left indicates reso-  \nFigure 4. Same as Fig. 3, but for the 3:2 resonance. The video can be downloaded on Github: https://raw.githubusercontent.com/ shuohuangGIT/Infer-migration-history/main/q2.mp4 .  \n<!-- image -->  \nlow this line, some systems are still evolving and resonance trapping is only temporary.  \nWe indicate the trapping criterion derived from the pendulum model in blue and the improved trapping criterion (Eq. (26)) with the orange line. From Fig. 3 and 4 it is clear that the pendulum model criterion for resonance trapping (blue line) fails to quantitatively match the numerical simulations. Our new criterion (orange line), however, fits the simulations perfectly. The equilibrium value of sin 𝜙 1 for the simulation snapshots is calculated by averaging its value over a time span of 0 GLYPH<149> 1 𝜏 𝑎 before and after the snapshot time, e.g., 0 GLYPH<149> 1 GLYPH<0> 0 GLYPH<149> 3 𝜏 𝑎 for the snapshot at 𝑡 = 0 GLYPH<149> 2 𝜏 𝑎 . The values of sin 𝜙 1 GLYPH<150> eq is indicated by grey solid lines (contour) in Fig. 3 and 4. sin 𝜙 1 GLYPH<150> eq increases as getting closer to the orange line, which is also expected by Eq. (25). The picture of resonance trapping/crossing over the entire parameter space of the migration time-scale and eccentricity damping timescale has now been clarified. Migration plays a role in exciting the planet's eccentricity, while eccentricity damping reduces it. On one hand, if a planet pair is in resonance, the eccentricity damping balances its excitation and finally librates near the equilibrium value. If the migration speed is so fast that there is no steady state solution for the resonance angle 𝜙 1 GLYPH<150> eq (Eq. (26)), the resonance is crossed. Otherwise, resonance trapping is ensured. On the other hand, if eccentricity is excited to be high enough, planets can be captured into resonance, but only temporarily, because of the continuous increase of the resonance libration amplitude (overstability, cf. Goldreich & Schlichting 2014). Although not evident from the figures presented, all simulations located in the large green corner above the red dashed line exhibit permanent libration of sin 𝜙 1 , indicating that the planets are captured in resonance permanently. The amplitude of sin 𝜙 1 during libration increases as the ratio 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 decreases, and approaches the red dashed line denoting overstability in Fig. 3 and 4. The green region situated below the red dashed line represents simulations in which planets are temporarily captured in resonance, with their amplitude of sin 𝜙 1 increasing over time and circulating at the end of the simulation. In conclusion, both efficient eccentricity damping ( 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 is high) and slow migration are required for permanent resonance trapping.  \n2  \nFigure 5. Planet mass versus stellar mass in the sample. Planet pairs with period ratios larger than ' 𝑗 , 1 'GLYPH<157> 𝑗 and smaller than ' 3 𝑗 , 2 'GLYPH<157> ' 3 𝑗 GLYPH<0> 1 ' are chosen. We take 𝑗 equal to 1, 2, 3 and 4. Planet masses and their host masses are indicated by black dots, with 1 𝜎 error bar. If the planet mass is inferred from its radius through the M-R relation, its uncertainty of the masses is assigned to be 𝜎 m = 0 GLYPH<149> 374 dex. The grey dashed line indicates the position where the planet-to-star mass ratio equals 𝜇 = 10 GLYPH<0> 4 . We assume planets to the left of this line to follow Type I migration and otherwise Type II. Only Type I migrating planets are included in our analysis.  \n<!-- image -->", '4 MCMC ANALYSIS': 'In this section, we first discuss how we select our sample in Sect. 4.1. ThenweconductanMCMCfittingtoconstrainthemodelparameters.', '4.1 Sample selection': 'Our sample selection and all analysis are based on the NASA Exoplanet Archive. Our attention is drawn to planets detected through transits and TTV. As shown in Sect. 2.3, we need planets\' masses to calculate the equilibrium eccentricities and period ratios in resonance. If only the radius is available, the M-R relationship described in Sect. 2.4 is used to calculate the masses of those transiting planets. Additionally, the semi-major axes and stellar masses are extracted.  \nWe do not consider short-period planets in our sample. We take this simple step to reduce the effects of both photo-evaporation and stellar tides on the M-R relationship. First, photo-evaporation can alter the M-R relationship for planets with low-density atmospheres on time-scale of 1 Gyr (Fulton & Petigura 2018). It is believed to have triggered the so-called \'radius gap\' (Fulton et al. 2017). When planets get closer to their host stars, this effect is more obvious (Fulton & Petigura 2018). Second, stellar tides alter close-in planets\' orbital properties (Lithwick & Wu 2012; Batygin & Morbidelli 2013; Charalambous et al. 2018; Papaloizou et al. 2018) and blur the information of the planets inherited from their protoplanet disc. Both mechanism are very sensitive to planets\' semi-major axes. Excluding the planets with a cutoff period shorter than 5 days, though crude, can  \nsuppress the interference from stellar tides (Choksi & Chiang 2020) and photo-evaporation (Fulton & Petigura 2018) on our sample.  \nWe display the planet mass versus their host mass in Fig. 5. The average planet mass is 10 𝑀 GLYPH<8> and the average stellar mass is 1 𝑀 GLYPH<12> . The planets\' and stars\' mass uncertainties are indicated by error bars. If the planet mass is inferred from the M-R relation, the log-normal standard deviation is then 𝜎 𝑚 = 0 GLYPH<149> 374 (Eq. (17)). The migration speed of low-mass planets in the proto-planetary disc scales linearly with planet mass, as dictated by the Type I migration limit. As planets become massive enough to perturb their surrounding disc, their migration gradually switches to Type II (Kanagawa et al. 2018; Pichierri et al. 2022). We set the boundary between two types of migration as 𝜇 trans = 10 GLYPH<0> 4 . Since our interests focus on Type I migration only planets with 𝜇 GLYPH<159> 𝜇 trans are included in our sample.  \nFinally, the sample size is reduced to 371 and the period ratios for all planet pairs are given by Fig. 6. The selected planets come from systems with two and more planets, including those with resonance chains.  \nTo calculate the resonance offset, the resonance number 𝑗 is required. There are excesses of systems (\'peaks\') just wide of the integer period ratios, which is suggestive of resonances. We only consider first-order resonances: 2:1, 3:2, 4:3, and 5:4. We assume that planets with period ratios slightly larger than the integer ratios are potentially in resonance until the period ratio "hits" the third-order resonance, because the resonant interaction is weaker as they become further from exact commensurability. Planet pairs with period ratios larger than \' 𝑗 , 1 \'GLYPH<157> 𝑗 but smaller than \' 3 𝑗 , 2 \'GLYPH<157> \' 3 𝑗 GLYPH<0> 1 \' are possibly in \' 𝑗 , 1 \'GLYPH<157> 𝑗 resonance. Here, \' 3 𝑗 , 2 \'GLYPH<157> \' 3 𝑗 GLYPH<0> 1 \' is the location of the closest third-order resonance. The selection of a period ratio limit for identifying planets in 2:1 resonance may seem arbitrary, given that the period ratio can extend up to 2.5, where planet pairs are unlikely to be in resonance. However, a slightly smaller window for the 2:1 resonance would not affect our conclusions. Nonetheless, this choice is useful for identifying planets in 3:2, 4:3, and higher 𝑗 first-order resonances because planets located near these resonance locations are close to nearby higher-order resonances, and may therefore be more easily perturbed. The satisfied period ratio windows are highlighted in Fig. 10 (top panel) and the four lower panels zoom in on these four windows, for 𝑗 = 1, 2, 3, 4, where we instead show the distribution of the offset from resonance, Δ . Planets out of the windows may still be in first-order resonance, but their fraction must be very low and it is not covered by our analysis. We ignore other first-order resonances and all higher-order resonances.', '4.2 Implication on planet-disc interaction from MCMC': "Weuse emcee (Foreman-Mackey et al. 2013) to perform the MCMC analysis. We implement three different models:  \n- (i) General model. The disc structure is not specified and the MCMCis used to fit log 10 ' 𝐶 𝑒 ℎ 2 1au ' , 𝜎 Δ , 𝑠 , 𝑞 .\n- (ii) Irradiation model. Stellar irradiation is assumed to be the main heating source and the MCMC is used to fit log 10 𝐶 𝑒 and 𝜎 Δ .\n- (iii) Viscous model. Viscosity-driven accretion is assumed to be the main heating source and the MCMC is used to fit log 10 ' 𝐶 𝑒 ℎ 2 1au ' , 𝜎 Δ .  \nThe prior distribution of parameters log 10 ' 𝐶 𝑒 ℎ 2 1au ' , log 10 𝐶 𝑒 , 𝜎 Δ , 𝑠 , 𝑞 and the values we take for 𝑠 , 𝑞 and ℎ 1au for the three different models are shown in Table 1. The convergence of MCMC chains are checked. We make use of the criterion that MCMC converges if the autocorrection time is shorter than 1/50 times its chain length. We checked that our results all satisfied the convergence criterion.  \n<!-- image -->  \n3.00  \nFigure 6. Distribution of the period ratio of all observed planet pairs in our sample (green histograms). We highlight the four windows in the vicinity of four different first-order resonances: 2:1, 3:2, 4:3 and 5:4. The left boundary of each window is ' 𝑗 , 1 'GLYPH<157> 𝑗 and the right boundary is ' 3 𝑗 , 2 'GLYPH<157> ' 3 𝑗 GLYPH<0> 1 ' . The total number of pairs is 371 in the green histogram and 128 in the windows.Figure 7. Distribution of log 10 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 , the semi-major axis-to-eccentricity damping time-scale at the location of the inner planets averaged over all planet pairs, in General (blue), Irradiation (orange) and Viscous (green) model. We calculate this quantity based on the posterior distributions of the parameters in each model. Dashed lines indicate their median values.  \n<!-- image -->  \nFor the General model, we examine whether our method is capable to retrieve all the parameters in Appendix B. It turns out that almost all parameters are degenerate. Therefore, the fitted values for f log 10 ' 𝐶 𝑒 ℎ 2 1au ' GLYPH<150> 𝜎 Δ GLYPH<150> 𝑠GLYPH<150> 𝑞 g may not be reliable (Appendix B1). The result of the General model is shown and analysed in Appendix A. We also calculate the quantity log 10 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 , the semi-major axis-toeccentricity damping time-scale, at the location of the inner planets averaged over all planet pairs. This quantity, however, shows to be independent of the other parameters and can be reproduced within the 1 𝜎 error bar (Appendix B1).  \nWe calculate log 10 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 in all three models, and their distributions are shown in Fig. 7. Two key points can be made. First, different disc structures result in nearly identical distributions. The parameter log 10 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 is not sensitive to the assumed disc structure. Second, the value of log 10 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 - peaking at 4 and almost always larger than 3 - is high. The high semi-major axis-to-eccentricity damping time-scale ratio indicates that temporary capture (overstable libration) did not operate for the planets in our sample, which would require 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 GLYPH<25> 170 (Goldreich & Schlichting 2014) in Eq. (30).  \nBy specifying the disc structure - the Irradiation or Viscous  \nmodel - the parameters can be successfully retrieved within 1 𝜎 error bar (Appendix B2). We show the fit result of f log 10 𝐶 𝑒 GLYPH<150> 𝜎 Δ g and f log 10 ' 𝐶 𝑒 ℎ 2 1au ' , 𝜎 Δ g for the Irradiation and Viscous model in Fig. 8 and Fig. 9, respectively. The python package corner.py (ForemanMackey 2016) is used to generate the plots of the posterior distributions.  \nThe Irradiation model (Fig. 8) fits log 10 𝐶 𝑒 = 0 GLYPH<149> 24 , 0 GLYPH<149> 25 GLYPH<0> 0 GLYPH<149> 25 . For the viscous model, the disc aspect ratio is sensitive to the stellar accretion rate and disc opacity. Therefore, we fit the combination 𝐶 𝑒 ℎ 2 1au , and log 10 ' 𝐶 𝑒 ℎ 2 1au ' = GLYPH<0> 3 GLYPH<149> 60 , 0 GLYPH<149> 30 GLYPH<0> 0 GLYPH<149> 29 (Fig. 9). If we take ℎ 1au = 0 GLYPH<149> 0245, same value as the irradiation model, then log 10 𝐶 𝑒 = GLYPH<0> 0 GLYPH<149> 38 , 0 GLYPH<149> 30 GLYPH<0> 0 GLYPH<149> 29 . Increasing ℎ 1au (or 𝐿 ★ in Eq. (5)) would result in a smaller value of 𝐶 𝑒 ( 𝐶 𝑒 / ℎ GLYPH<0> 2 1au ). Theoretically, Tanaka et al. (2002) and Tanaka & Ward (2004) from the first principle calculate that log 10 𝐶 𝑒 = 0 for locally isothermal discs. The fitting outcomes from both models are consistent with it.  \nAdditionally, the fitted values for 𝜎 Δ are 0 GLYPH<149> 65 , 0 GLYPH<149> 10 GLYPH<0> 0 GLYPH<149> 08 and 0 GLYPH<149> 79 , 0 GLYPH<149> 10 0 GLYPH<149> 09 for the Irradiation and Viscous model, respectively. Their values are twice that of the mass dispersion. Indeed, we expect that the fitted 𝜎 Δ is of the same magnitude as 𝜎 m (Sect. 2.4). However, 𝜎 Δ is fitted to be slightly larger than our expectation. This could be an implication of turbulent discs (Rein & Papaloizou 2009; Goldberg & Batygin 2022) and/or post-disc perturbations (e.g. Lithwick & Wu 2012; Chatterjee & Tan 2014; Stock et al. 2020). We further run an MCMC fitting fixing 𝜎 Δ to 0.374, the resulting posterior distribution of log 10 𝐶 𝑒 or log 10 ' 𝐶 𝑒 ℎ 2 1au ' are not significantly different from what we present here. It gives log 10 ' 𝐶 𝑒 ℎ 2 1au ' = GLYPH<0> 3 GLYPH<149> 95 , 0 GLYPH<149> 29 GLYPH<0> 0 GLYPH<149> 27 for the Viscous model and log 10 𝐶 𝑒 = 0 GLYPH<149> 44 , 0 GLYPH<149> 19 GLYPH<0> 0 GLYPH<149> 20 for the Irradiation model.  \nIn summary, our MCMC model shows that eccentricity damping is highly effective (log 10 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 GLYPH<25> 4), making resonant over-stability unlikely. The observed period ratio excess of planets is consistent with predictions by Tanaka et al. (2002) and Tanaka & Ward (2004) ( 𝐶 𝑒 GLYPH<25> 1), irrespective of whether the disk structure is dominated by irradiation or viscous heating. However, the aspect ratio of a viscous inner disc depends on the disc opacity and stellar accretion rate (e.g. Liu et al. 2019), which limits our ability to constrain 𝐶 𝑒 .", '5 IMPLICATIONS FOR PLANET FORMATION': "In this section, we adopt the fit result from the Irradiation model and further study the implications of resulting resonant planets statistically. Ramos et al. (2017) and Charalambous et al. (2022) use similar  \nFigure 8. Corner plot of variables in the MCMC analysis (log 10 𝐶 𝑒 and 𝜎 Δ ) with 1 𝜎 , 2 𝜎 and 3 𝜎 confidence contours, for the Irradiation model, which fixes the surface density power law index 𝑠 = GLYPH<0> 15 GLYPH<157> 14 and the disc aspect ratio index 𝑞 = 2 GLYPH<157> 7. The 1 𝜎 uncertainty is labelled on the top of each column and indicated by left and right dashed lines. The middle dashed lines indicate their median values.  \n<!-- image -->  \nprescription for their disc structure. The reason why we choose the Irradiation model is the following. Even though resonance trapping can happen much earlier, planets' period ratios (offsets) are more evolved at the end of the disc lifetime when the migration and eccentricity damping time-scales are longer than the disc dispersal time-scale. The disc structure at this stage mostly determines what the corresponding mature planet system looks like. Because planet formation consumes solids and solids drift inward rapidly due to gas drag (Weidenschilling 1977; Andrews et al. 2012), the disc at this point becomes optically thin, rendering stellar irradiation the main heating source. The transition disc LkCa15 is arguably an example that low-mass planets can carve a large dust cavity (Leemker et al. 2022). Therefore, the Irradiation model is more applicable to transition discs.  \nThe best-fitting resonance offset distribution is plotted in Fig. 10 (blue lines), with upper and lower 3 𝜎 uncertainty (blue shaded regions). We assume that the offset of non-resonant pairs follows a uniform distribution, which is also indicated in Fig. 10 (grey dashed lines). The MCMC fits the observed distribution better than the uniform-only model because it fits more planets with small Δ and fewer planets with large Δ , just as observed. The complete sample with 128 pairs are fitted simultaneously. However, we use four panels to display the four near-resonance planets because Δ depends on the resonance number 𝑗 in a complex form. We cannot present one distribution of Δ to represent all 128 pairs while keeping the shape of log-normal profile.  \nFor each planet pair, we calculate its probability of being in resonance ( 𝑃 res GLYPH<150> 𝑘 , in Sect. 5.1). The total number of resonant pairs and planets' mean eccentricities in our sample is then obtained. The properties of their birthplace - the natal proto-planet disc - are then inferred, e.g., the upper limit of the surface density (Sect. 5.2) and the location of the migration barrier (Sect. 5.3).  \nFigure 9. Similar to Fig. 8, but for the Viscous model. The surface density power law index is fixed at 𝑠 = GLYPH<0> 3 GLYPH<157> 8 and the disc aspect ratio index at 𝑞 = GLYPH<0> 1 GLYPH<157> 16.  \n<!-- image -->", '5.1 Fraction of resonant pairs': "Given the fitted value for 𝐶 𝑒 and 𝜎 Δ , we can calculate the probability of each planet pair in resonance via two characteristic probability distribution functions (Eq. (19) and Eq. (20)): 𝑃 res GLYPH<150> 𝑘 = 𝑝 res ' Δ obs GLYPH<150> 𝑘 'GLYPH<157> » 𝑝 res ' Δ obs GLYPH<150> 𝑘 ' , 𝑝 n GLYPH<0> res ' Δ obs GLYPH<150> 𝑘 '… . In Fig. 10, we plot the histogram of period ratio and resonance offset distribution of planet pairs weighted by 𝑃 res GLYPH<150> 𝑘 (red hatches). It indicates the period ratio and resonance offset distribution of resonant pairs. The period ratio of resonant pairs peaks just wide of integer ratios, and, as the period ratio further increases, resonant pairs vanish. Not all pairs with period ratios close to integer ratio are in resonance.  \nThe total number of resonant pairs is the summation of resonant probability over all planet pairs:  \n𝑁 res = 𝑁 pairs ∑︁ 𝑘 = 1 𝑃 res GLYPH<150> 𝑘 GLYPH<149> (31)  \nWe label the average fraction of resonant pairs, 𝑓 res = 𝑁 res GLYPH<157> 𝑁 pairs , on the top right of four lower panels in Fig. 10. 𝑁 pairs is the number of pairs in the narrow period ratio windows.  \nWe also calculate the fraction of all resonant pairs among all pairs in our sample: 𝐹 res = 𝑁 res GLYPH<157> 371 = 14 GLYPH<149> 8 , 0 GLYPH<149> 5 GLYPH<0> 0 GLYPH<149> 7 %. The distribution is shown in Fig. 11 left panel. It is consistent with the crude estimation made by Wang & Ji (2014) (10%-20%). This number ignores higherorder resonances and first-order resonances with resonance numbers larger than 4 (5:4). The resonant fraction could therefore be higher.  \nWe split the resonant planets into three groups, each group has a number of resonant pairs 𝑁 1 , 𝑁 2 , 𝑁 3 , respectively, and 𝑁 res = 𝑁 1 , 𝑁 2 , 𝑁 3 (see Sect. 2.5 for detail and Fig. 2 for a sketch). Three different resonant fractions are calculated:  \n- GLYPH<15> 𝑓 res ' Mig j out ' = 𝑁 1 GLYPH<157> ' 𝑁 1 , 𝑁 2 , 𝑁 5 ' : fraction of Migrating pairs among the pairs where the inner planets migrate slower than the outer;\n- GLYPH<15> 𝑓 res ' Brk j out ' = 𝑁 2 GLYPH<157> ' 𝑁 1 , 𝑁 2 , 𝑁 5 ' : fraction of Braking pairs amongthepairswheretheinnerplanets migrate slower than the outer;  \nFigure 10. Distribution of period ratio (top panel, similar to Fig. 6) and resonance offset Δ (bottom panel) for the inferred resonant planet pairs (red hatches) versus all observed planet pairs in our sample (green histogram). The lower panels display the distribution of period ratio offset Δ within the highlighted four windows in the top panel. Grey dashed lines show the probability distribution of all planet pairs assuming they follow a uniform distribution (not in resonance). We compare it with the best MCMC fit result (blue curve, after bin size correction) assuming the Irradiated disc model. The shaded blue region indicates the 3 𝜎 uncertainty. We also label the fraction of inferred resonant planets, which is the ratio of the area between the red hatches and the green histograms (see Sect. 5.1 for details).  \n<!-- image -->  \nP  \nout/  \nP  \n- GLYPH<15> 𝑓 res ' Brk j inn ' = 𝑁 3 GLYPH<157> ' 𝑁 3 , 𝑁 4 ' : fraction of Braking pairs among the pairs where the inner planets migrate faster than the outer.  \nThis classification allows us to compare the fraction of resonant pairs under different physical conditions. The distributions of the three fractions are shown in Fig. 11, right panel. It shows that 𝑓 res ' Mig j out ' GLYPH<25> 0 GLYPH<149> 27. This implies that in some pairs the ambient gas disc disperses before the pair reaches a migration barrier, i.e., either planet migration is slow or the disc disperses rapidly following planet migration and formation. The gas-poor formation scenario for sub-Neptunes (Dawson et al. 2015; Choksi & Chiang 2020) would be consistent with this picture. However, a still larger fraction of resonant pairs, 𝑓 res ' Brk j out ' GLYPH<25> 0 GLYPH<149> 23, and 𝑓 res ' Brk j inn ' GLYPH<25> 0 GLYPH<149> 41 reach their migration barriers. It implies that gas-poor formation applies to some of the observed systems, but not all of them. Interestingly, 𝑓 res ' Brk j out ' is smaller than 𝑓 res ' Brk j inn ' , i.e., the inner planets tend to be more massive than the outers. A possible explanation would be that the inner edge of the disc - the location of the pressure maximum - would also be the place where pebbles accumulate. Pebble accretion at such locations can be very efficient (Chatterjee & Tan 2014; Jiang & Ormel 2023).  \nWealso plot the eccentricity distribution for the planets in resonant pairs, weighted by their probability of being in resonance ( 𝑃 res GLYPH<150> 𝑘 ), in Fig. 12. Generally, the outer planets have smaller eccentricities consistent with Eq. (11). Their values mostly fall between 10 GLYPH<0> 3 and 10 GLYPH<0> 2 . It suggests that if we observe a sub-Neptune planet pair with relatively high eccentricities ( 𝑒 GLYPH<24> 0 GLYPH<149> 1), they are not likely to be in  \nfirst-order resonance irrespective of their near-resonance period ratio. Post-disc perturbations (Choksi & Chiang 2022) could, however, excite eccentricities and change resonant pairs from apsidal antialignment to alignment (Laune et al. 2022). These apsidally aligned systems would have slightly larger eccentricities than what Fig. 12 predicts.", '5.2 Upper limit on the disc surface density at resonance trapping': "The constraints the MCMC model provides cannot be used to determine the absolute value of the natal disc surface density, as it cancels in the 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 expression. An upper limit for the natal disc surface density of resonant planets can, however, be deduced from the resonance trapping criterion. We use Eq. (26) in order to break the degeneracy and to find an upper limit for the disc surface density at the trapping location. It is assumed in the derivation of Eq. (26) that the outer planet is on a fixed circular orbit. Such an assumption is valid because the outer planets on average have lower eccentricities than their inner siblings (Fig. 12). From Sect. 4.2, we already found that 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 ¡ 10 3 for observed transiting planets, that is, overstability is not likely to occur. In this regime, Eq. (26) alone gives the resonance trapping condition. Therefore, we can get the critical migration time-scale, below which the resonance would be crossed. In the Type I migration regime, planet migration speed is proportional to gas disc surface density. We are therefore able to obtain the upper limit of the disc gas surface density.  \nFigure 11. Resulting distribution of resonant fraction for all pairs in our sample (left) and detailed resonant distribution for pairs with period ratio within the near integer ratio windows (right). Here 𝑓 res ' Mig j out ' is the fraction of planets that enter resonance during migration with the outer planet migrating faster. 𝑓 res ' Brk j out ' is the fraction of resonant planets stopping at the disc inner edge with the outer planet migrating faster. 𝑓 res ' Brk j inn ' is the fraction of resonant planets stopping at the disc inner edge with the inner planet migrating faster.  \n<!-- image -->  \nThe upper limit is reached by combining Eq. (7), Eq. (9), Eq. (16) and Eq. (26):  \nΣ max ' 𝑟 2 ' = 1 GLYPH<149> 6 j 𝑓 1 j ℎ 3 ' 𝑟 2 ' GLYPH<20> ' 𝑗 , 1 ' 𝐶 𝑒 𝛾 𝐼 𝜏 𝑒 1 𝜏 𝑒 2 GLYPH<21> 1 2 GLYPH<18> 𝑟 2 𝑀 1 𝑟 1 𝑀 2 GLYPH<19> 5 4 𝑀 ★ 𝑟 2 2 (32)  \ninserting a disc model and extrapolating to 1 au, we obtain:  \nΣ 1au GLYPH<150> max = 1 GLYPH<149> 6 j 𝑓 1 j ℎ 3 1au »' 𝑗 , 1 ' 𝐶 𝑒 𝛾 𝐼 … 1 GLYPH<157> 2 GLYPH<18> 𝑀 1 𝑀 2 GLYPH<19> 3 GLYPH<157> 4 GLYPH<2> GLYPH<18> 𝑟 1 𝑟 2 GLYPH<19> 2 𝑞 GLYPH<0> 𝑠 GLYPH<157> 2 GLYPH<0> 1 GLYPH<149> 5 GLYPH<16> 𝑟 2 1 𝑎𝑢 GLYPH<17> 3 𝑞 GLYPH<0> 𝑠 𝑀 ★ 𝑟 2 2 GLYPH<149> (33)  \nThe expression Σ 1au GLYPH<150> max / ℎ 3 1au is consistent with what Kajtazi et al. (2022) found in their simulations. In Fig. 13, we plot the Σ 1au GLYPH<150> max for each planet pair that is possible in resonance, versus their host mass. Thedeeperthecolour, the more likely it is that the pair is in resonance. The size of the symbol indicates the resonance index 𝑗 . The figure shows that high𝑗 resonances tend to be associated with high surface densities. This result is in line with Type-I migration theory. Massive discs result in faster migration, which allows the planets to cross the relatively strong resonance. In addition, Fig. 13 shows that Σ max increases with stellar mass. That is because the migration speeds in the Type I limit depend on the star-to-disc mass ratio. Higher surface density is required to migrate faster. We indicate the surface density of the Minimum mass solar nebula (i.e., MMSN, Hayashi 1981) and the Minimum mass extra-solar nebula (i.e., MMEN, Chiang & Laughlin 2013) in Fig. 13. The upper limits of nearly all disc surface densities are below that of the MMEN, while the inferred disc surface densities are centred around the MMSN value.  \nWe conclude this discussion with two final points. First, the value for Σ max we obtained refers to the time when the planets were locked into resonance, not the upper limit over the entire disc lifetime. Second, planets in higher 𝑗 resonances tend to provide higher upper limits on the disc surface density. However, they might alternatively have formed in close proximity to each other, avoiding crossing of lower𝑗 resonances. In that case, the true value of Σ 1au is likely to be less than Σ 1au GLYPH<150> max .", '5.3 Migration barrier reflects the disc inner rim': "Planets can migrate in the disc, but their migration is believed to be halted somewhere as otherwise all planets would be consumed by the host star. However, the location and mechanism of the migration barrier are under debate. There are two main explanations for the barriers: dust sublimation (Kama et al. 2009; Flock et al. 2019) and the stellar magnetosphere (Königl et al. 2011; Hartmann et al. 2016).  \nIf we assume a mass-luminosity relation for main sequence stars (Duric 2004), the dust sublimation radius becomes:  \n𝑟 subl = 𝑟 sub GLYPH<150> 0 GLYPH<18> 𝑀 ★ 𝑀 GLYPH<12> GLYPH<19> 2 GLYPH<150> (34)  \nwhere 𝑟 sub GLYPH<150> 0 is the silicate dust sublimation radius for solar mass stars:  \n𝑟 sub GLYPH<150> 0 = 0 GLYPH<149> 13 GLYPH<18> 𝐶 bw 1 GLYPH<19> 1 GLYPH<157> 2 GLYPH<18> 0 GLYPH<149> 1 𝜖 GLYPH<19> 1 GLYPH<157> 2 GLYPH<18> 1400K 𝑇 d GLYPH<19> 2 au GLYPH<150> (35)  \nwhere 𝐶 bw is the back-warming factor, 𝜖 is dust cooling efficiency and 𝑇 d is the dust sublimation temperature (Kama et al. 2009).  \nThe magnetospheric infall radius tends to expand and converge to the stellar corotation radius due to angular momentum locking (Long et al. 2005). The stellar corotation radius is expressed as:  \n𝑟 co = 0 GLYPH<149> 057 GLYPH<18> 𝑃 ★ 5 d GLYPH<19> 2 3 GLYPH<18> 𝑀 ★ 𝑀 GLYPH<12> GLYPH<19> 1 3 au (36)  \nThe rotation periods of T Tauri stars (e.g. Bouvier et al. 2007; Lee & Chiang 2017) and young star associations (e.g.,Upper Sco and NGC 2264 Roquette et al. 2021) are several days (1-10 days).  \nThe MCMC model is capable to identify (in a probabilistic sense) whether or not the planets have reached a migration barrier. We plot the location of the inner planet of those pairs - the presumed location of the migration barrier - in Fig. 13. The symbol size again represents the resonance index 𝑗 and the colour indicates how likely they are in resonance. From the plot, it can be seen that as the stellar mass increases the location of the migration barrier moves further away from the star. Motivated by the two theories about the migration barrier and aiming to figure out which radius is more consistent with observation, we fit the location of a migration barrier using the single  \nFigure 12. Eccentricity distribution of planet pairs that are in resonance, weighted by the probability of being in resonance. We compare the eccentricity distribution of the outer planets in all planet pairs with their inner planets.  \n<!-- image -->  \npower law expression:  \n𝑟 bar GLYPH<150> 𝑘 = 𝑟 0 GLYPH<18> 𝑀 ★GLYPH<150> 𝑘 𝑀 GLYPH<12> GLYPH<19> 𝑙 GLYPH<150> (37)  \nwhere 𝑟 0 is the location of the migration barrier for solar-mass stars. The Gaussian likelihood we construct is weighted by the resonance probability 𝑃 res GLYPH<150> 𝑘 :  \nln L = 𝑁 pairs ∑︁ 𝑘 = 1 ln h N' 𝑟 bar GLYPH<150> 𝑘 GLYPH<0> 𝑟 𝑘 GLYPH<150> 𝜎 2 𝑟 ' GLYPH<1> 𝑃 res GLYPH<150> 𝑘 i GLYPH<150> (38)  \nwhere 𝑟 bar GLYPH<150> 𝑘 , 𝑟 𝑘 and 𝜎 𝑟 are the model predicted migration barrier location for the 𝑘 -th pair (Eq. (37)), and the presumed location of the migration barrier and the standard deviation of the fitted migration barrier radius, respectively. If one pair has a larger 𝑃 res GLYPH<150> 𝑘 , it is more likely that the planets are formed in the protoplanet disc and undergo disc migration. The MLE fits 𝑟 0 = 0 GLYPH<149> 13, 𝑙 = 1 GLYPH<149> 49 and 𝜎 𝑟 is 0.079 au. The location of the migration barrier fitted by single power law relation is, however, shallower than the corresponding index of the dust sublimation radius but steeper than that of the magnetospheric radius. It may imply that there are two planet populations whose migration barriers are carved by either dust sublimation radius or magnetospheric radius. For this reason, we also fit a broken power law:  \n𝑟 bar GLYPH<150> 𝑘 = 8 > > > < > > > 𝑟 1 GLYPH<16> 𝑀 ★GLYPH<150> 𝑘 𝑀 GLYPH<12> GLYPH<17> 1 GLYPH<157> 3 if 𝑀 ★GLYPH<150> 𝑘 GLYPH<20> 𝑀 crit 𝑟 2 GLYPH<16> 𝑀 ★GLYPH<150> 𝑘 𝑀 GLYPH<12> GLYPH<17> 2 if 𝑀 ★GLYPH<150> 𝑘 ¡ 𝑀 crit GLYPH<150> (39)  \n:  \nwhere 𝑀 crit = ' 𝑟 1 GLYPH<157> 𝑟 2 ' 3 GLYPH<157> 5 𝑀 GLYPH<12> is the transition mass between the two different power laws. Similar to the single power law, we use MLE to fit the two parameters 𝑟 1 = 0 GLYPH<149> 08 au and 𝑟 2 = 0 GLYPH<149> 12 au and the transition mass 𝑀 crit is 0 GLYPH<149> 79 𝑀 GLYPH<12> (see Fig. 13). The low mass fit corresponds to a corotation period of 𝑃 ★ = 8 GLYPH<149> 76 days, which is consistent with the observed rotation period for young stars (Bouvier et al. 2007; Lee & Chiang 2017; Roquette et al. 2021). The value of 𝑟 2 agrees with Eq. (35). Around these high-mass stars, planets are trapped at the dust sublimation location as it exceeds 𝑟 co. However, the uncertainty 𝜎 𝑟 = 0 GLYPH<149> 078 au is only slightly smaller than the value given by the single power law, suggesting that the broken power law model is  \nonly marginally better. A larger sample will make for a more reliable analysis.", '6 DISCUSSION': 'In this work, we have constructed a model that connects the planet migration history to the observed values of the offset from integer period ratios ( Δ ). If the resulting planet pair is in resonance, Δ follows a log-normal distribution. On the other hand, if it is nonresonant, the corresponding Δ is assumed to follow a uniform distribution. Based on this, we have developed a statistical model that constrains the migration histories of the observed planets by conducting a Markov Chain Monte Carlo method (MCMC) analysis. We examine our MCMC method using self-generated mock data Sect. B and prove that it can indeed reproduce certain features.  \nOur model for resonance trapping is designed for two-planet systems and first-order resonances, but we also include observed systems with planets in a resonance chain. Kajtazi et al. (2022) have shown that the averaged properties of the resonance chain (involving three or more planets) still reflect the properties of the system as if there is only one resonant pair. Therefore, multi-planet systems do not significantly contaminate the results. In addition, we have not considered higher-order resonances. A more general model that applies to both first-order and higher-order resonances needs to be considered in the future. Finally, it is possible that the inner planets migrate across the inner disc edge and enter the disc cavity (Huang & Ormel 2022; Fitzmaurice et al. 2022), where our model would not be applicable. But those planets are plausibly massive enough to open a deep gap (Ataiee & Kley 2021; Chrenko et al. 2022), which are excluded by our sample selection.  \nA key assumption in the model is that the uncertainties in the planet masses and the ensuing Δ follow a log-normal distribution. We can then fit the excess of period ratio just wide of integer ratio with a log-normal profile, thus extracting pairs in resonance. If we would adopt different distributions for the mass, the resulting distributions for Δ would become far more complex and no longer allow us to express the corresponding likelihoods in closed form. However, certain post-disc dynamics, e.g., post-disc energy dissipation from planetesimal scattering (Chatterjee & Ford 2015; Ghosh & Chatterjee 2022), stellar tides (Lithwick & Wu 2012; Batygin & Morbidelli 2013) and stellar encounters (Cai et al. 2019; Stock et al. 2020), could slightly change the period ratios of planet pairs. Therefore, they may play roles in broadening, shifting, or even skewing the log-normal profile. What the resulting Δ distribution may look like needs to be investigated in future work. Once addressed, one may learn the postdisc perturbation histories the planets have experienced. However, this also requires a much larger sample size than what we have at present, as already in this work the MCMC is unable to break some model degeneracies.  \nOur model mainly applies to those small planets unable to open a gap (in Type I migration regime). Tanaka et al. (2002) and Tanaka &Ward (2004) showed that the semi-major axis damping time-scale ( 𝜏 𝑎 ) and eccentricity damping time-scale ( 𝜏 𝑒 ) have a relation for locally isothermal disc: 𝜏 𝑒 = 1 GLYPH<149> 28 𝐶 𝑒 ℎ 2 𝜏 𝑎 GLYPH<157> 𝛾 𝐼 and 𝐶 𝑒 = 1. If the planet partially opens a gap, the migration speed decreases linearly with the gas surface density in the gap (Kanagawa et al. 2018). The same holds for eccentricity damping (Pichierri et al. 2022). Therefore, the ratio 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 is independent of surface density and the constraints we obtain on it still hold when the planet opens a partial gap. Neither assuming a specific disc nor migration model, we obtain that 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 GLYPH<25> 10 4 , which is the most robust result of this  \nr [  \nau  \n]  \n0  \nFigure 13. Inferred maximum surface density at 1 au for the resonant planets (top panel) and location of the migration barrier (bottom panel) versus the host stellar mass, for the 128 planet pairs in the four windows in Fig. 10. Different sizes of the circles indicate different resonance index 𝑗 . The colour of the symbols indicates the probability of planets in resonance ( 𝑃 res GLYPH<150> 𝑘 ). In the top panel, error bars indicate the 1 𝜎 uncertainties in the stellar mass and Σ 1au GLYPH<150> max - the latter follows from the uncertainty in 𝐶 𝑒 . The left and right grey dashed lines indicate the surface density of the MMSN and MMEN, respectively. In the bottom panel, grey and green lines give the best-fitting location of the migration barrier (the location of the inner planet in resonant pairs) fitted by a single power law and a broken power law, respectively.  \n<!-- image -->  \nstudy. After adopting the irradiation disc model, we further obtain 𝐶 𝑒 GLYPH<25> 1 which is consistent with Tanaka et al. (2002) and Tanaka & Ward (2004). On the other hand, Charalambous et al. (2022) argue that 𝐶 𝑒 = 0 GLYPH<149> 1 because lower 𝐶 𝑒 increases the resonance offsets. A potential reason is that their simulations have assumed that trapping takes place at 1 au, where the disc aspect ratio is relatively large, while most transiting planets are found at GLYPH<24> 0 GLYPH<149> 1 au. Furthermore, by improving the analytical criterion for resonance trapping, we are able to constrain the upper limit on the natal disc surface density for those planets in resonance. The resulting maximum surface density is similar to that of the Minimum Mass Solar Nebular (MMSN) but smaller than that of the Minimum Mass Extra-solar Nebula (MMEN).  \nSeveral other migration prescriptions have been proposed, including more sophisticated ones such as those by Paardekooper et al. (2010, 2011). These prescriptions demonstrate that planets within a certain mass range can be naturally trapped at a location where the (positive) corotation torque exerted on the planet exceeds the (negative) Lindblad torque (Bitsch et al. 2013, 2014; Baruteau et al. 2014). To first order, for small planets, the migration behavior would still be  \npredominantly linear (with planet mass and disk mass) except near these trapping locations. In our model, such a scenario is naturally incorporated through the Braking pair. Regarding the dependence of the damping terms on planet eccentricity, non-linear correction terms have been proposed by e.g., Cresswell & Nelson (2006, 2008); Ida et al. (2020). However, at low eccentricity, these non-linear terms are irrelevant and we do not include these terms in our investigation. Finally, our findings on the ratio 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 are robust, irrespective of the specific disc migration prescriptions and non-linear terms. This is because this ratio has a one-to-one correspondence to the resonance offset Δ . Therefore, our conclusion regarding 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 remains valid even if we incorporate different migration models (disc structures) or non-linear terms.  \nThe imprint resonant trapping leaves behind tentatively allow us to assess where and when planets form in the discs. The evaluation of the probability of planet pairs in MMR is crucial, and this step also informs us of the number of resonant pairs. Since we employ different models for how planets get trapped in resonance (Sect. 2.5), it is possible to distinguish gas-poor formation scenarios (pairs trapped in  \nresonance when migrating) and gas-rich formation scenarios (pairs stopped by the migration barrier). The former could be identified with late formation, while the latter, which are more dominant, are connected to the early formation in gas-rich discs. The present orbits of these planets further hint at the location of migration barriers. However, due to our small sample, we cannot unambiguously identify the physical origin of the migration barrier; either the dust sublimation radius or the magnetospheric radius would fit the data. The situation is, however, expected to improve in the near future. The upcoming launch of the PLAnetary Transits and Oscillations of stars (PLATO) mission (e.g. Rauer et al. 2014) and The Earth 2.0 (ET) mission (Ge et al. 2022; Ye 2022) could drastically increase not only the number but also the precision of planet detections. It will provide us with a more precise analysis of the planet formation and migration histories reflected in the dynamical properties of resonant planets.', '7 CONCLUSIONS': 'We manage to construct a statistical model connecting planet migration theory to observed quantities of Kepler planets. Based on the inferred masses and resonance offsets, we conduct an MCMC analysis to extract the history of planet-disc interaction from planet-planet dynamics. The statistical approach provides us with the following findings:  \n- (i) The semi-major axis-to-eccentricity damping time-scale ratio can be constrained at log 10 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 GLYPH<25> 4 with a dispersion of GLYPH<25> 0 GLYPH<149> 3 dex, irrespective of the assumed disc model. The eccentricity damping is so efficient that overstable libration of resonances is unlikely to have occurred.\n- (ii) After assuming the isothermal irradiation disc model, we obtain that log 10 𝐶 𝑒 = 0 GLYPH<149> 24 , 0 GLYPH<149> 25 GLYPH<0> 0 GLYPH<149> 25 , which connects eccentricity damping and semi-major axis damping in the Type-I migration theory. It is consistent with Tanaka et al. (2002) and Tanaka & Ward (2004).\n- (iii) FromtheMCMCposterior,theprobability that a planet pair is in resonance follows. The fraction of transit planet pairs in first-order MMRamounts to 14 GLYPH<149> 8 , 0 GLYPH<149> 5 GLYPH<0> 0 GLYPH<149> 7 %.\n- (iv) Most of the inferred resonant planets are consistent with the scenario that they reached a migration barrier, indicative of early migration in a gas-rich disc. The location of the migration barrier could be the dust sublimation radius for massive stars ( 𝑀 ★ ¡ 0 GLYPH<149> 79 𝑀 GLYPH<12> ) and the magnetospheric radius for low mass stars ( 𝑀 ★ GLYPH<159> 0 GLYPH<149> 79 𝑀 GLYPH<12> ).\n- (v) Byevaluating the resonance strength of those inferred resonant planets, the upper limit of the proto-disc surface density during the planet formation era is obtained. Most systems have their Σ 1au below that of the Minimum Mass Extra-solar Nebula (MMEN) and half of them below that of the Minimum Mass Solar Nebula (MMSN).  \n(vi) It is found that the classical MMR trapping/crossing criterion based on the pendulum model does not match numerical simulation. We provide and numerically verify an improved criterion (Eq. (26)) based on the equilibrium resonance angle, which together with the overstability condition of Goldreich & Schlichting (2014) (Eq. (30)) fully describes the problem.  \nFuture work could feature a more complete model accounting for high-order resonances and resonance chains. Detection of more resonant planets and more precise measurements of planetary and stellar properties with upcoming missions will definitely improve the confidence of our analysis.', 'ACKNOWLEDGEMENTS': 'The authors appreciate the thoughtful comments of the referee, Dr. Carolina Charalambous. S.H. would like to thank Wei Zhu, Mario Flock, Simon Portegies-Zwart, Martĳn Wilhelm, and Remo Burn for their useful discussions. The authors acknowledge support by the National Natural Science Foundation of China (grant no. 12250610189). This work has made use of the NASA Exoplanet Archive. Software: emcee (Foreman-Mackey et al. 2013), corner.py (ForemanMackey 2016), Matplotlib (Hunter 2007; Caswell et al. 2021), REBOUND (Rein & Liu 2012) and REBOUNDx (Tamayo et al. 2020).', 'DATA AVAILABILITY': 'The data underlying this article will be shared on reasonable requests to the corresponding author.', 'APPENDIX A: RESULTS OF THE GENERAL DISC MODEL': "The full MCMC fitting corner plot of the General model is shown in Fig. A1. The fitted values for all four parameters are log 10 ' 𝐶 𝑒 ℎ 2 1au ' = GLYPH<0> 3 GLYPH<149> 55 , 0 GLYPH<149> 71 GLYPH<0> 0 GLYPH<149> 89 , 𝑠 = 1 GLYPH<149> 58 , 0 GLYPH<149> 58 GLYPH<0> 1 GLYPH<149> 05 , 𝑞 = 0 GLYPH<149> 31 , 0 GLYPH<149> 37 GLYPH<0> 0 GLYPH<149> 47 . However, these values do not necessarily refer to the real physical disc parameters, due to the degeneracy among them. Therefore, we generate mock planet period ratio data (with the same size as our sample used here) and examine this MCMC model without specifying the disc structure in Sect. B1. Those tests indeed suggest that 𝑠 and 𝑞 are not able to be retrieved.  \nFig. A1 reveals the degeneracy: 𝑠 negatively correlates with log 10 ' 𝐶 𝑒 ℎ 2 1au ' , 𝑞 correlates positively with log 10 ' 𝐶 𝑒 ℎ 2 1au ' and 𝑞 positively with 𝑠 . We found that such correlations are well represented by:  \n4 𝑞 GLYPH<0> 𝑠 = 𝐶 2 GLYPH<150>  \n(A1)  \nfrom Eq. (13) and  \nlog 10 ' 𝐶 𝑒 ℎ 2 1au ', 2 𝑞 log 10 𝑟 2 1 𝑎𝑢 GLYPH<0> log 10 ' 0 GLYPH<149> 78 𝛾 𝐼 ',' 4 𝑞 GLYPH<0> 𝑠 ' log 𝑟 1 𝑟 2 = 𝐶 1 GLYPH<150> (A2)  \nfrom Eq. (15), Eq. (14) and Eq. (16). 𝐶 1 and 𝐶 2 are two constants. Eq. (A2) and Eq. (A1) are indicated in Fig. A1 (green dashed lines), and they match the correlation. The fitted log-normal dispersion 𝜎 Δ , however, does not have any correlation with other parameters. It  \nFigure A1. Corner plot of variables in the MCMC analysis (log 10 𝐶 𝑒 , 𝜎 Δ , 𝑠 , and 𝑞 ) with 1 𝜎 , 2 𝜎 and 3 𝜎 confidence contours, w/o assuming the disc structure. We also calculate log 10 ' 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 ' accordingly and add this variable into the corner plot. The 1 𝜎 uncertainty is labelled on the top of each column and indicated by left and right dashed lines. The middle dashed lines indicate their median values. Green lines indicate the correlation between different variables, which is Eq. (A2) and Eq. (A1).  \n<!-- image -->  \ne  \nmakes sense because our model does not depend on 𝜎 Δ . We introduce another parameter:  \nlog 10 𝜏 𝑎 𝜏 𝑒 = 𝐶 𝑒 ℎ ' 𝑟 1 ' 2 0 GLYPH<149> 78 𝛾 𝐼 (A3)  \nand it can be calculated given f log 10 ' 𝐶 𝑒 ℎ 2 1au ' GLYPH<150> 𝜎 Δ GLYPH<150> 𝑠GLYPH<150> 𝑞 g . Here, 𝑟 1 is the observed semi-major axis of inner planet average over all planet pairs. This parameter is calculated and put in the corner plot. It shows that log 10 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 is a quantity independent of 𝑠 and 𝑞 , but slightly depends on log 10 ' 𝐶 𝑒 ℎ 2 1au ' .", 'APPENDIX B: MCMC PERFORMANCE EXAMINATION': 'Asdescribed in Sect. 4.1, we make use of planet masses, period ratios, semi-major axes and host masses from the NASA exoplanet database. Given these data as well as disc structure, we can then calculate the exact period ratios if planets are in resonance, and compare them to the observations. In this section, to examine the performance of the MCMC, we generate mock samples by replacing the actual period ratios with those assuming they are in or out of resonance. We randomly select a fraction of planet pairs to be in resonance and the  \nfraction is 𝑓 res. If they are not in resonance, the resulting period ratio follows a uniform distribution. If they are in resonance, the resulting period ratio is calculated via Eq. (12). We finally add log-normal noise to the planet masses with 𝜎 𝑚 , which enables us to compare it to the resulting 𝜎 Δ . In this way, the mock sample is generated, with known disc parameters. We then run the MCMC model to check whether we can reproduce the input parameters. The mock sample has the same size as the real sample we used in Sect. 4.', 'B1 Tests without assuming a disc structure': "We first examine whether the MCMC can reproduce log 10 ' 𝐶 𝑒 ℎ 2 1au ' , 𝑠 and 𝑞 , and whether the fitted 𝜎 Δ is comparable to (the input) 𝜎 𝑚 . Several sets of parameters are used for generating samples and model examination. The default parameters for generating the mock sample are log 10 ' 𝐶 𝑒 ℎ 2 1au ' = GLYPH<0> 2 GLYPH<149> 94, 𝜎 𝑚 = 0 GLYPH<149> 374, 𝑠 = GLYPH<0> 15 GLYPH<157> 14, 𝑞 = 2 GLYPH<157> 7 and 𝑓 res = 0 GLYPH<149> 5. We then test different values of 𝑠 , 𝑞 , and 𝑓 res. The resulting corner plots are shown in Fig. B1, Fig. B2 and Fig. B3. The true value of parameters for generating mock samples is labelled at the top right of each panel and indicated by blue lines.  \nUnfortunately, the true values of log 10 ' 𝐶 𝑒 ℎ 2 1au ' , 𝑠 and 𝑞 are not properly retrieved. The posterior distribution gives the expected val-  \nFigure B1. Similar to Fig. A1. Examination of MCMC model performance, without specifying the disc structure. In each panel, all parameters used are labelled in the upper right corner. We change the parameter 𝑠 (labelled after F ) in different panels to examine the performance of MCMC model and we change 𝑠 from GLYPH<0> 0 GLYPH<149> 5 and GLYPH<0> 2. In each corner plot, blue lines show the true values. The 1 𝜎 uncertainty is labelled on the top of each column and indicated by left and right dashed lines. The middle dashed lines indicate their median values.  \n<!-- image -->  \ne  \n<!-- image -->  \ne  \n0.24  \n<!-- image -->  \nFigure B2. Similar to Fig. B1, but with different 𝑞 . We change 𝑞 from 0 and 1.  \n<!-- image -->  \nues for log 10 ' 𝐶 𝑒 ℎ 2 1au ' , 𝑠 and 𝑞 , but they are not consistent with the true values. However, the correlation between the parameters is revealed. The corner plot shows that 𝑠 negatively correlates with log 10 ' 𝐶 𝑒 ℎ 2 1au ' , 𝑞 is positively correlated to log 10 ' 𝐶 𝑒 ℎ 2 1au ' and 𝑞 is positively correlated to 𝑠 . Their correlation is consistent with Eq. (A2) and Eq. (A1), indicated by green dashed lines.  \nSimilar to Fig. A1, we also plot the posterior distribution of log 10 ' 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 ' (Eq. (A3)), which is a rather independent variable. The  \nexamination result suggests that log 10 ' 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 ' can be fitted within 1.5 𝜎 error bar in all cases. The fitted 𝜎 Δ is always very close to the log-normal error 𝜎 𝑚 we impose for planet mass. It proves that the distribution of resonance offset Δ resulting from log-normal distributed planet masses also follows a log-normal distribution, when they are in resonance.  \nWe therefore conclude that our MCMC model is useful for fitting  \n0.19  \ne  \n<!-- image -->  \n0.17  \nFigure B3. Similar to Fig. B1, but with different 𝑓 res. We change fraction of pairs in resonance ( 𝑓 res) from 50% to 100%.  \n<!-- image -->  \nlog 10 ' 𝜏 𝑎 GLYPH<157> 𝜏 𝑒 ' but not log 10 ' 𝐶 𝑒 ℎ 2 1au ' , 𝑠 and 𝑞 if a disc structure has not been specified.", 'B2 Tests assuming a disc structure': 'Wehere examine whether the MCMC can reproduce log 10 𝐶 𝑒 , given a disc structure, and whether the fitted 𝜎 Δ is still comparable to 𝜎 𝑚 . We test several sets of parameters. The default parameter set is log 10 𝐶 𝑒 = 0 GLYPH<149> 28, 𝜎 𝑚 = 0 GLYPH<149> 374, 𝑠 = GLYPH<0> 15 GLYPH<157> 14, 𝑞 = 2 GLYPH<157> 7 and 𝑓 res = 0 GLYPH<149> 5. We then change each parameter to two other values, while keeping the other parameters the same.  \nThe MCMC results are shown in Fig. B4. Different panels show the fit result for the mock sample generated from different parameter sets, and the true values of the parameters are labelled on the top right and indicated by blue lines. The true value of log 10 𝐶 𝑒 can always be reproduced within 1.5 𝜎 when 𝑓 res ¡ 0 GLYPH<149> 5. The fitted 𝜎 Δ is always consistent with 𝜎 𝑚 as well. However, if we decrease the fraction of planets in resonance 𝑓 res to 0.1, the log10 𝐶 𝑒 can no longer be retrieved and therefore the MCMC is no longer valid. However, when fitting the observed data, we have 𝑓 res GLYPH<25> 0 GLYPH<149> 5 in most cases (Fig. 11 right panel). Therefore the MCMC results for the observed data are reliable.  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.  \ne  \nFigure B4. Examination of MCMC model performance, fixing the disc structure (the values of 𝑠 and 𝑞 ). The top panel shows the fiducial model parameters. In each panel, all parameters used are labelled in the upper right corner. We then change every single parameter separately (labelled after F ) in each panel in the lower half of the figure. We change log 10 𝐶 𝑒 to GLYPH<0> 1, 𝜎 2 𝑚 to 0 GLYPH<149> 07 and 0 GLYPH<149> 3, 𝑠 to GLYPH<0> 0 GLYPH<149> 5 and GLYPH<0> 2, 𝑞 to 0 and 1 and 𝑓 res to 10% and 100%. In each corner plot, blue lines show the true value of log 10 𝐶 𝑒 and 𝜎 m. The 1 𝜎 uncertainty is labelled on the top of each column and indicated by left and right dashed lines. The middle dashed lines indicate their median values.  \n<!-- image -->'}