Resonant Chains and the Convergent Migration of Planets in Protoplanetary Disks

An increasing number of compact planetary systems with multiple planets in a resonant chain have been detected. The resonant chain must be maintained by convergent migration of the planets due to planet–disk interactions if it is formed before the dispersal of the protoplanetary gas disk. For type I migration in an adiabatic disk, we show that an analytic criterion for convergent migration can be developed by requiring that any part of the resonant chain should be convergently migrating toward the remaining part. The criterion depends primarily on the logarithmic gradients α and β of the surface density and temperature profiles of the disk, respectively, and it is independent of the absolute values of the surface density and temperature. The analytic criterion is applied to the Kepler-60, Kepler-80, Kepler-223, TOI-178, and TRAPPIST-1 systems. Due to the variation of planetary masses within the resonant chains, we find that convergent migration typically requires rather extreme values of (α, β) that have little or no overlap with common disk models. Finally, we show that there is an empirical relationship between the distance of the innermost planet from the central star and the stellar mass for the observed resonant chain systems, which supports the idea that the resonant chains are formed and maintained by stalling the migration of the innermost planet near the inner edge of the disk truncated by the magnetic fields of the protostar.


INTRODUCTION
The first mean-motion resonance (MMR) in an extrasolar planetary system was discovered around the star GJ 876, with planets b and c in 2:1 resonance (Marcy et al. 2001).After two additional planets were discovered in this system, it also became the first system with a three-body resonance, with planets b, c, and e in a 4:2:1 Laplace resonance (Rivera et al. 2010;Batygin et al. 2015;Millholland et al. 2018).Many other systems with three or more planets in resonant chains have since been discovered, particularly in transit surveys such as Kepler and TESS.Some of the famous resonant chain systems are Kepler-60 (Steffen et al. 2013), Kepler-80 (MacDonald et al. 2016;Shallue & Vanderburg 2018), Kepler-223 (Mills et al. 2016), and TRAPPIST-1 (Gillon et al. 2017;Luger et al. 2017).
Convergent migration is required for MMR capture (Goldreich 1965;Murray & Dermott 1999).It is believed that planet-disk interaction is the dominant migration mechanism for assembling MMR and resonant chains of planets.Lee & Peale (2002) have shown that disk-induced migration is able to put GJ 876 b and c into the observed MMR with suitable eccentricity damping.For some of the observed resonant chains, migration simulations have also successfully assembled multiple planets into the observed resonant chains (e.g., Tamayo et al. 2017;Delisle 2017;MacDonald & Dawson 2018).However, as pointed out by Papaloizou et al. (2018), most of these simulations applied inward migration on the outermost planet only.Such a prescription ensures convergent migration and almost always results in a resonant chain system if the migration is sufficiently slow.The inner planets will be captured into MMR, and the resonant chain will be assembled from the outermost to the innermost planet.However, this may not reflect the real situation where all planets are embedded in the disk and are migrating individually.
Disk-induced planetary migration can be classified primarily into type I and type II, depending on whether the planetary body is massive enough to open a gap in the co-orbital region (Kley & Nelson 2012).Type I migration is valid for the Earth-mass planets involved in many observed resonant chains.Since type I torque depends on planetary mass and various disk properties (e.g., Paardekooper et al. 2010), we should examine what kind of protoplanetary disk would allow convergent migration.Batygin (2015) has derived a criterion for convergent type I migration of two planets in a classical Mestel disk, where m 1 and m 2 are the masses of the inner and outer planets, respectively, and a 1 and a 2 refer to their orbital semimajor axes. 1 This criterion assumes (1) a disk surface density profile Σ = Σ 0 a 0 /a, where Σ 0 is the surface density at the reference semimajor axis a 0 ; (2) a constant disk aspect ratio h = H/a throughout the disk; and (3) an isothermal equation of state.These assumptions make the criterion appear independent of the disk parameters.Equation (1) says that the outer planet has to be more massive than the inner planet for convergent inward migration and the subsequent MMR assembly to happen.Such criterion on the planetary masses is not satisfied in many observed resonant chain systems.Take TRAPPIST-1 as an example: the planet h has a mass of 0.3M ⊕ , which is significantly less massive than the inner planet g of 1.34M ⊕ (Grimm et al. 2018).Similar examples are seen in other resonant chain systems.This implies that these resonant chains are unlikely to be assembled in a classical Mestel disk.Regardless of how a resonant chain is assembled, if it is formed in the gas disk phase and coexists with the disk, the resonant chain must be maintained during migration, which requires convergent migration within the chain.In this study, we consider type I migration in an adiabatic disk and investigate what disk properties are required to maintain a resonant chain by convergent migration.The adopted disk and migration models are described in Section 2. In Section 3, we analytically derive the necessary conditions for the convergent migration of a pair of planets, which are shown to be accurate by comparison with numerical results.The criterion is generalized to the migration of multiple planets in a resonant chain in Section 4. The analytic criterion is tested numerically for the four-planet resonant chain system Kepler-223 and then applied to other observed resonant chains.We find that rather extreme disk parameters, far from the usual disk models, are required for convergent migration.The effects of uncertainties in planetary masses are discussed in Section 5. Terquem & Papaloizou (2007), Cossou et al. (2013), Brasser et al. (2018), and Huang & Ormel (2022) have shown that an inner disk edge can stall migration and facilitate the assembly of an MMR pair and resonant chain (see also Li 2014 andBatygin &Morbidelli 2020 for a similar scenario for the Laplace resonance of the Galilean satellites of Jupiter).In Section 6, we find a simple relationship between the mass of the host star and the orbital semimajor axis of the innermost planet in the observed resonant chains, which supports this idea.Our conclusions are summarized in Section 7.

DISK AND MIGRATION MODELS
We adopt an adiabatic disk model with surface density profile Σ = Σ 0 (a 0 /a) α and temperature profile T = T 0 (a 0 /a) β .Since the scale height H = c s /n, where the sound speed c s ∝ T 1/2 and the mean motion n ∝ a −3/2 , the reduced scale height h = H/a = h 0 (a/a 0 ) (−β+1)/2 .Thus, the disk has five parameters (α, β, γ, Σ 0 , h 0 ).Hereafter, we assume that the adiabatic index γ = 5/3.One can also think of α and β as the local power-law indices of the surface density and temperature profiles: Paardekooper et al. ( 2010) have derived a formula for the unsaturated nonlinear type I migration torque from an adiabatic disk on the planet: where and m and m 0 refer to the masses of the planet and the star, respectively.There are three components to the torque in Equation ( 2): (1) the first term with the coefficient −2.5−1.7β+0.1α is the Lindblad torque, (2) the second term with the coefficient 1.1(3/2 − α) is the barotropic term of the nonlinear horseshoe torque, and (3) the last term is the entropy-related term of the nonlinear horseshoe torque.Paardekooper et al. (2011) have studied the effects of viscous and thermal diffusion on the nonlinear corotation torque in Equations ( 2)-(4).Equations ( 2)-( 4) are valid if the viscous saturation parameter p ν and the thermal saturation parameter p χ are ∼ 0.5.In the limit p ν and p χ ≪ 1, the nonlinear horseshoe torque is replaced by the linear corotation torque.In the limit p ν and p χ ≫ 1, the corotation torque saturates, and only the Lindblad torque remains.The viscous saturation parameter p ν = 2 kx 3 s /3, where k = a 2 n/(2πν), x s = 1.1γ −1/4 (m/m 0 ) 1/2 h −1/2 is the dimensionless half-width of the horseshoe region, and ν is the viscosity.For α ν viscosity with ν = α ν c s H (Shakura & Sunyaev 1973) . For the planets in the resonant chain systems that we consider below, m/m 0 = 1.0-4.6 × 10 −5 and p ν = 0.27-0.85(αν /0.001) −1/2 (h/0.05)−7/4 .The thermal saturation parameter p χ = 3p ν ν/χ/2 = a 2 nx 3 s /(2πχ), where χ = 16γ(γ − 1)σ SB T 4 /(3κΣ 2 n 2 ) is the thermal diffusion coefficient, σ SB is the Stefan-Boltzmann constant, and κ is the opacity.The value of p χ depends on the actual values of a, n, Σ, and κ, not just the dimensionless parameters m/m 0 and h.If we assume that Σ 0 = 1700 g cm −2 and h 0 = 0.05 at a 0 = 1 au, similar to the minimum mass solar nebula (MMSN), and that κ = 1 cm 2 g −1 , p χ = 0.21-0.65 for the planets under consideration.Based on these estimates of p ν and p χ , the unsaturated nonlinear type I migration torque should be applicable to the planets considered in this paper.
The forced migration rate from the tidal torque on the planet is where the subscript f indicates forced migration.The correction from eccentricity (Goldreich & Schlichting 2014) is neglected to simplify derivation later, and it has little effect on the discussion in this paper on the convergent properties of migration.
If the expression in the square brackets in Equation ( 3) is positive (negative), there is a positive (negative) torque and outward (inward) migration.The boundary between inward and outward migration is simply deduced by setting Γ in Equation (3) to zero: This line is shown in Figure 1(a) in the (α, β) space with γ = 5/3.Inward (outward) migration occurs in the region above (below) the line.Note that the condition for inward/outward migration is independent of the mass m and location a of the planet.

CONDITIONS FOR CONVERGENT MIGRATION OF TWO PLANETS
For an inner planet of mass m 1 and an outer planet of mass m 2 , irrespective of whether the forced migration is inward or outward, the dividing line between the situation where the forced migration of the inner planet is faster than that of the outer planet and vice versa is 2)-( 5), the torque exerted on a planet of mass m at a is Γ ∝ Γ 0 ∝ m 2 h −2 Σa 4 n 2 ∝ m 2 a −α+β , and the forced migration rate ( ȧ/a) f ∝ ma −α+β−1/2 .Thus, the condition or where P 1 and P 2 are the orbital periods.This line is shown in Figure 1(b) in the (α, β) space for an example of two equal-mass planets, and it is independent of the period ratio because ln(m 2 /m 1 ) = 0.The forced migration of the inner (outer) planet is faster in the region above (below) the line.Note that Equation (7) for the Mestel disk with constant h (i.e., α = β = 1) agrees with Equation (1).Two planets undergo convergent migration if their migration rates satisfy for outward migration.Thus, the convergent migration zone is bounded by two lines, the disk constraint in Equation ( 6) and the planetary mass constraint in Equation ( 7) or (8).This is shown in Figure 1(c) for the example of two equal-mass planets, with the inward and outward convergent migration regions labeled.We stress again that convergent migration is only a necessary condition to maintain MMR.The capture probability also depends on the eccentricities and the migration and eccentricity damping rates (Mustill & Wyatt 2011;Goldreich & Schlichting 2014;Deck & Batygin 2015;Batygin & Petit 2023;Huang & Ormel 2023;Kajtazi et al. 2023).
A system with m 2 more massive than m 1 should favor convergent inward migration, while a system with a more massive m 1 should favor convergent outward migration.This intuition is verified with Equation (8) and demonstrated in Figure 2. In Figure 2, the dotted and solid lines have the same meaning as those in Figure 1 but now for P 2 /P 1 = 3/2 and m 2 /m 1 = 2 in panel (a) and m 2 /m 1 = 1/2 in panel (b).We see a large convergent inward migration region and no convergent outward migration region within the range of (α, β) shown in panel (a), and vice versa in panel (b).6) for an adiabatic disk with power-law indices α for surface density and β for temperature.(b) Boundary between faster forced migration of inner planet and faster forced migration of outer planet according to Equations ( 7) and ( 8) for m 1 = m 2 .(c) Combined plot, with the convergent outward and inward migration regions corresponding to the regions bounded by the red and blue lines in Figure 3.
We have tested the validity of the analytic conditions in Equations ( 6), (7), and (8) using numerical simulations.We start with a system with two planets already in 3:2 MMR, which is obtained from a prior migration simulation.For the masses, we adopt 1.125M ⊙ for the central star (which is the stellar mass of Kepler-223; see below) and 1M ⊕ for both planets.The initial orbital elements are listed in Table 1.The planets are initially in 3:2 MMR with low eccentricities, antialigned periapses, and both MMR angles in libration.For the disk parameters, we assume that Σ 0 = 1700 g cm −2 and  h 0 = 0.05 at a 0 = 1 au, similar to the MMSN.The disk parameters α and β are surveyed uniformly in a 30 × 30 grid from −3 to +3.The adopted masses, Σ 0 , and h 0 only affect the absolute migration rate but not the sign of the relative migration rate and the locations of the convergent and divergent migration zones in (α, β) space.We use the SyMBA integrator (Duncan et al. 1998) with imposed migration following Equations ( 2)-( 5).An eccentricity damping of ė/e = −100| ȧ/a| is applied on both planets, which avoids the instability caused by high eccentricity during convergent migration.
For each combination of α and β, the system is integrated for one migration timescale a 2 / ȧ2 of the outer planet.At the end of the simulation, we check whether the system is still in the 3:2 MMR by examining the libration of the resonant angles.
The simulation results are shown in Figure 3.The green dots represent the systems where the 3:2 MMR is broken during the simulation.The red dots indicate the systems where the planets move outward and the MMR is maintained.The blue dots indicate the systems where the planets move inward and the MMR is maintained.The solid lines are the same lines shown in Figure 1(c) from the analytic conditions in Equations ( 6), (7), and (8).The analytic theory and numerical simulations match quite well.Even systems just outside the zones bounded by the two lines break from 3:2 MMR, which means that a small relative divergent migration is enough to break the MMR within one migration timescale.

GENERALIZATION AND APPLICATION TO RESONANT CHAIN SYSTEMS
In this section, we generalize the analytic criterion for convergent migration of a pair of planets derived in the previous section to the case of multiple planets in a resonant chain.The Kepler-223 system is taken as an illustrative example to verify the analytic criterion with numerical simulations.We then apply the analytic criterion to the Kepler-60, Kepler-80, TOI-178, and TRAPPIST-1 systems.The results show that rather extreme values of (α, β) are required to maintain the resonant chains in these systems.

Convergent Migration Criterion for Multiple Planets in Resonant Chains
To construct the analytic criterion, one should distinguish between the actual migration rate of a resonant chain and the forced migration rate from the torque applied by the disk on each planet.If a resonant chain is maintained by convergent migration within the chain, the actual migration rate ȧ/a must be the same for all planets and the same as the comigration rate of the whole chain.For a resonant chain of N planets, the comigration rate ( ȧ/a) 12...N can be calculated by equating the work done by the torques on the individual planets to the total orbital energy change of the chain, with the assumption that the orbital period ratios are fixed: where ( ȧi /a i ) f is the forced migration rate of the planet i of mass m i and orbital period P i .Given a certain set of disk parameters (α, β, γ, Σ 0 , h 0 ), we can calculate the individual planet's forced migration rate ( ȧi /a i ) f from Equations (3)-( 5) and the comigration rate ( ȧ/a) 12...N from Equation (9) for any planet pair, triplet, or the whole resonant chain in a disk, under the assumption that they are locked in resonance.
It is straightforward to generalize the necessary criterion to maintain a resonant pair to a resonant chain.The idea is that any part of the resonant chain should be convergently migrating toward the remaining part.Taking a resonant chain of four planets as an example, the generalized criterion for convergent inward migration is This means that planet 1 should be migrating at a rate slower than that of the triplet 2-3-4.Similarly, the migration rate of the planet pair 1-2 should be slower than that of the planet pair 3-4, and the migration rate of the triplet 1-2-3 should be slower than that of planet 4. For outward migration, the inequalities are flipped, i.e.,  The set of inequalities for inward or outward migration have to be satisfied simultaneously because no part of the chain should be divergently moving away from the rest of the chain.
A similar criterion for a planet pair leads to Equations ( 7) and ( 8) for the combination α − β in the exponent.For a resonant chain, since the actual migration rate of two or more planets is given by the more complicated Equation (9), it is not possible to reduce the criterion to a similarly simple form.Each part of the resonant chain in Equations ( 10) and ( 11) results in an equation (e.g., | ȧ/a| 12 = | ȧ/a| 34 ) for α−β, which can be solved numerically and combined with the disk constraint in Equation ( 6) to determine the inward and outward convergent migration regions in the (α, β) space.The constraints from all parts of the resonant chain can then be combined to give the overlapping regions (if any) where the resonant chain can be maintained by convergent migration.

Application to Kepler-223
Let us use the 8:6:4:3 resonant chain in Kepler-223 as an example.We adopt a stellar mass of 1.125M ⊙ and the planetary masses as listed in Table 2 from Mills et al. (2016).Panels (a), (b), and (c) of Figure 4 show the convergent migration zones in the (α, β) space between planet b and triplet c-d-e, between pair b-c and pair d-e, and between triplet b-c-d and planet e, respectively.The inward convergent zone is shaded in blue, and the outward convergent zone is shaded in red.In each panel, one can see that the convergent zones have the same shapes as those for two planets in Figures 1 and 2. This is because one of the lines defining the convergent zones -the boundary between inward and outward migration in Equation ( 6) -is the same in all cases, while the other line from a comparison of the migration rate of a part of the resonant chain with the migration rate of the rest of the resonant chain (Equations ( 10) and ( 11)) is simply shifted vertically (i.e., a different value of α−β) for each part of the resonant chain.In panel (c), one can see that there is no convergent inward migration zone for triplet b-c-d and planet e in the region of (α, β) shown.This is because planet e has a small mass of 4.8M ⊕ , which makes it difficult for planet e to catch up with b-c-d in inward migration, unless extreme values of (α, β) are adopted.Panel (d) is the combined constraint, which shows the overlapping region from panels (a), (b), and (c).The analytic criterion indicates that the resonant chain in Kepler-223 cannot be maintained by inward convergent migration for any values of α and β between −3 and +3 and that it can be maintained by outward convergent migration in a narrow region of mostly negative values of α and β (i.e., both surface density and temperature increasing outward).Numerical simulations with imposed migration have been carried out to test the analytic criterion in Figure 4. We start with a system with four planets assembled into the desired 8:6:4:3 resonant chain by a prior migration simulation.The orbital parameters are listed in Table 2.The three-body resonance or Laplace angles θ L,in = −λ 1 + 2λ 2 − λ 3 and θ L,out = λ 2 − 3λ 3 + 2λ 4 (which are denoted as  Here λ i is the mean longitude of planet i.The disk parameters Σ 0 and h 0 , the values of α and β surveyed, and the eccentricity damping applied during the simulations are the same as those used in Section 3.Each simulation is integrated for one migration timescale of the outermost planet or 1 Myr, whichever is smaller. The simulation results are shown in Figure 5.As in Figure 3, the red and blue dots represent the systems where the resonant chain is maintained by convergent outward and inward migration, respectively, and the green dots represent the systems where the resonant chain is broken.For all the green dots, the period ratios evolve significantly from the original period ratios, and there are no ambiguous situations where the period ratios remain close to the original ones but the resonance angles are not librating.The region bounded by the two red solid lines is the convergent outward migration zone found analytically and shown in Figure 4(d).All of the simulations within this analytically derived region are indeed able to maintain the resonant chain by convergent outward migration.However, there are also some systems outside the analytic region with a sustained resonant chain.There are some red dots just below the line defined by the constraints in Equation ( 11), and there are five blue dots near the extension of the red line representing the boundary between inward and outward migration given by Equation (6).
Figure 6 shows the simulation with α = −2.586and β = −1.966,as an example of the simulations that are just below the analytic region but still able to maintain the resonant chain by outward convergent migration.The evolution of the orbital semimajor axes, eccentricities, and resonance  angles is shown in the upper left, lower left, and right panels, respectively.In addition to the three-body Laplace angles θ L,in and θ L,out , we also show the two-body MMR angles (with θ 12,in = 3λ 1 −4λ 2 +ϖ 1 and θ 12,out = 3λ 1 −4λ 2 +ϖ 2 , where ϖ i is the longitude of periapse of planet i, for planets 1 and 2, etc.).As the four planets migrate outward by ∼ 10%, the eccentricities are rather stable, with slightly decreasing amplitudes of variation toward the end of the integration, and all resonance angles remain in libration with no obvious change in the libration centers.It is not yet clear why this and other cases just below the analytic region are able to maintain the resonant chain during outward migration.One possibility is that the nonadjacent planets are also in first-order MMR in this system (both the b-d and c-e pairs are in 2:1 MMR), but further investigation is needed.
Figure 7 shows the evolution of the semimajor axes, eccentricities, and resonance angles of the simulation with α = 0.931 and β = 1.552.This is one of the blue dots near the red line representing the boundary between inward and outward migration, i.e., the line of no migration.The eccentricities decrease over the 1 Myr integration, and some of the resonance angles show an increase in libration amplitude around 0.4 Myr.But the migration is indeed very slow, so the differential migration is not significant enough to break the resonant chain.
The numerical results show that the analytic criterion for the resonant chain in Kepler-223 is not as accurate as in the case of two planets in MMR.Nevertheless, the region in (α, β) where the resonant chain is maintained during outward migration is only slightly larger, and the cases of sustained resonant chains during inward migration are special cases near the line of no migration with very slow migration.Therefore, the analytic criterion is still a good indicator of the types of disk that can maintain a resonant chain during migration.If the resonant chain in Kepler-223 migrated any significant amount after assembly, the migration should be inward because the observed planets are close to the star (with P b = 7.38 days).But our analysis shows that there is no region in (α, β) where the resonant chain is maintained by inward migration.Even if we allow for outward migration, the steady-state, constant α ν -viscosity disk models (dashed line in Figure 5) are outside the region where the resonant chain can be maintained by outward migration.

Application to Other Resonant Chain Systems
Now we apply the analytic criterion to several other resonant chain systems.For convenience, we label the planets in a resonant chain of N planets in numerical order from 1 for the innermost planet to N for the outermost planet.
For Kepler-60, the resonant chain consists of three planets in 5:4:3 resonance (Steffen et al. 2013).We adopt planetary masses of m 1 = 4.1M ⊕ , m 2 = 4.8M ⊕ , and m 3 = 3.8M ⊕ (Goździewski et al. 2016).Stellar mass is not needed to determine the locations of the convergent migration regions.Figure 8(a) shows the analytic result, which is the overlapping region from the convergent migration zones between pair 1-2 and planet 3 and between planet 1 and pair 2-3.There is no inward convergent migration zone that can maintain the resonant chain, as in Kepler-223.The outward convergent migration zone is slightly smaller than that for Kepler-223 and requires negative values of α and β.
Kepler-80 has six known planets, with all but the innermost planet f in a resonant chain.The outermost planet g cannot be included in our analysis, because its mass is not known (Shallue & Vanderburg 2018).The middle four planets are in 9:6:4:3 resonance.We adopt planetary masses of m 1 = 6.75M ⊕ , m 2 = 4.13M ⊕ , m 3 = 6.93M ⊕ , and m 4 = 6.74M ⊕ (MacDonald et al. 2016)

8(b)
shows that the outward convergent migration zone is similar to that of Kepler-223 and that there is a small region of inward convergent migration around α ∼ 2 and β ∼ 3 that can maintain the resonant chain.
TOI-178 is also a six-planet system with all but the innermost planet b in a resonant chain.The resonance is 18:9:6:4:3.We adopt planetary masses of m 1 = 4.77M ⊕ , m 2 = 3.01M ⊕ , m 3 = 3.86M ⊕ , m 4 = 7.72M ⊕ , and m 5 = 3.94M ⊕ (Leleu et al. 2021).Figure 8(c) shows that both inward and outward convergent migration zones are smaller than those of Kepler-80.Since the type I forced migration rate is proportional to planetary mass, the particularly massive m 4 in the middle of the chain makes it more difficult for the whole chain to maintain convergent migration.
For the seven-planet system TRAPPIST-1, we adopt planetary masses of m 1 = 1.374M ⊕ , m 2 = 1.308M ⊕ , m 3 = 0.388M ⊕ , m 4 = 0.692M ⊕ , m 5 = 1.0391M ⊕ , m 6 = 1.321M ⊕ , and m 7 = 0.326M ⊕ (Agol et al. 2021).Although the mean motions of the inner two planets b and c are nearly in the ratio 8:5:3 with that of planet d, it is not clear that they are part of the chain of first-order resonances (9:6:4:3:2) of the outer five planets (Gillon et al. 2017;Luger et al. 2017).Figure 8(d) shows the analytic result if we assume that all seven planets are in the resonant chain.Due to the innermost planet b (m 1 ) being the most massive planet in this system, the outward convergent migration zone is the largest for all the systems considered.There is no inward convergent migration zone that can maintain the resonant chain, because the outermost planet h (m 7 ) is significantly less massive than the rest of the planets (except m 3 ), which makes it difficult for m 7 to catch up with the rest of the planets in inward migration.If we assume that planets b and c are not in the resonant chain, we find that no region in (α, β) space can maintain the resonant chain in either inward or outward convergent migration, because both the innermost planet d (m 3 ) and the outermost planet h (m 7 ) of the resonant chain are now significantly less massive than the rest of the planets.
As shown in Section 4.2, the analytic criterion is a good indicator of the types of disk that can maintain a resonant chain during migration.even though the regions of inward and outward convergent migration in (α, β) may be slighter larger than the analytic ones.The analytic outward migration zone is the largest for TRAPPIST-1 if we assume that all seven planets are in the resonant chain.But even in this case, it just touches the steady-state, constant α ν -viscosity disk models (dashed line in Figure 8).In all cases, the resonant chain can be maintained by inward convergent migration in, at most, a small region of (α, β).As in the case of Kepler-223, the planets in the systems discussed in this section are close to the stars, and the difficulty with inward migration is particularly problematic.

EFFECTS OF UNCERTAINTIES IN PLANETARY MASSES
All of the resonant chain systems examined above were discovered by the transit method, and the uncertainties in the planetary masses can be as much as ∼ 40%.As the type I migration rate of a planet scales linearly with its mass, the mass uncertainties can significantly affect the estimation of convergent migration zones in (α, β).In this section, we examine how the convergent migration zones are changed if we adjust the planetary masses within 1σ.To enhance the inward convergent migration zone, we lower the mass of the innermost planet by 1σ and raise the mass of the outermost planet by 1σ.For the other planets in the chain, their masses are adjusted within 1σ to archive a relatively smooth mass variation.
For Kepler-223, the adjusted masses are m 1 = 6.3M ⊕ , m 2 = 6.8M ⊕ , m 3 = 6.7M ⊕ , and m 4 = 6.2M ⊕ .The analytic convergent migration zones with the adjusted masses are shown in Figure 9(a).Compared to the nominal mass result shown in Figure 4(d), the outward convergent migration zone is almost identical, but there is now a small inward convergent migration zone.For Kepler-60, we adjust the masses to m 1 = 3.4M ⊕ , m 2 = 3.8M ⊕ , and m 3 = 4.9M ⊕ .With m 3 > m 2 ∼ m 1 , the outward convergent migration zone in Figure 8(a) is replaced by a large inward migration zone in Figure 9(b).For Kepler-80, the adjusted masses are m 1 = 6.24M ⊕ , m 2 = 4.94M ⊕ , m 3 = 6.23M ⊕ , and m 4 = 7.97M ⊕ .Comparing the results with nominal masses in Figure 8(b) and adjusted masses in Figure 9(c), the inward migration zone is larger, while the outward migration zone is smaller.For TOI-178, we adjust the masses to m 1 = 4.09M ⊕ , m 2 = 3.01M ⊕ , m 3 = 3.86M ⊕ , m 4 = 6.20M ⊕ , and m 5 = 5.25M ⊕ .Comparing the uncertainties in the masses reported by Agol et al. (2021) are sufficiently small (6% or less) that the results are almost the same for any variation of the masses within ±1σ: the analytic convergent migration zone is almost the same as in Figure 8(d) if all seven planets are in the resonant chain, and no region in (α, β) space can maintain the resonant chain in either inward or outward convergent migration if planets b and c are not in the resonant chain.
These examples show that even if we take into account the observational uncertainties in the planetary masses to optimize the size of the (α, β) region where the resonant chain can be sustained by inward convergent migration, there is very little or no overlap between the inward convergent migration zone (shaded in blue in Figure 9) and the steady-state, constant α ν -viscosity disk models (dashed line with α + β = 3/2 in Figure 9) in most cases.In fact, within the observational uncertainties, there is no inward convergent migration zone for TRAPPIST-1 due to the small mass of the outermost planet h.The only exception is Kepler-60, where, within uncertainties, it is possible to have a substantial overlap between the inward convergent migration zone and the steady-state, constant α ν -viscosity disk models.Therefore, a more accurate determination of planetary masses in the resonant chain systems would further constrain their origin and migration in a disk.

FORMATION OF RESONANT CHAIN NEAR THE INNER DISK EDGE
We have demonstrated that the observed resonant chains are difficult to maintain in reasonable disk models if they migrate any significant amount after assembly.In this section, we show that a relationship between the orbital semimajor axis of the innermost planet and the stellar mass of the observed resonant chains supports the idea that they have a common origin in a stalling and assembly mechanism near the inner edge of the protoplanetary disk.
In a simple model of a disk with an inner edge, the surface density is zero at the inner edge at a tr , rises to a maximum at some distance from the inner edge, and then decreases away from the star.This means that the power-law index α = −d ln Σ/d ln a would change from highly negative near the inner edge to zero at the surface density maximum and then positive far from the inner edge.For any reasonable variation of the temperature power-law index β = −d ln T /d ln a, there is a zerotorque location at some distance from the inner edge, where α and β satisfy Equation ( 6), the type I migration torque is zero, and a planet would stop migrating.At distances closer to (farther from) the star, the torque would be positive (negative), and a planet would migrate outward (inward).As the innermost planet migrates inward, it would be trapped at the zero-torque location.The outer planets would continue their inward migration and would be naturally captured into resonance with the inner planets, regardless of the variation in planetary masses.As the resonant chain is assembled, the innermost planet would move slightly inward from the zero-torque location, so that its outward migration would balance the inward migration of the outer planets.
This idea of forming the resonant chain near the inner disk edge is also supported by the positions of the innermost planets in the observed resonant chain systems.In Figure 10, we show the orbital semimajor axis of the innermost planet versus the stellar mass for two populations of planets.The blue filled circles are the innermost planets in 11 confirmed and suspected resonant chain systems.In addition to and TRAPPIST-1, we also include HD 40307 (Díaz et al. 2016), Kepler-79 (Jontof-Hutter et al. 2014;Yoffe et al. 2021), K2-32 (Heller et al. 2019;Lillo-Box et al. 2020), K2-138 (Christiansen et al. 2018;Lopez et al. 2019), V1298 Tau (David et al. 2019;Suárez Mascareño et al. 2021), and YZ Ceti (Stock et al. 2020).Some of these are suspected but not confirmed resonant chains, and HD 40307 is suspected to be in a resonant chain in the past (Papaloizou & Terquem 2010).GJ 876 and HR 8799 are two resonant chain systems not included in this sample because of the presence of giant planets that undergo type II instead of type I migration.The open circles are the innermost planets in other systems with three or more planets, with red and green open circles for the planets discovered by transit and radial velocity, respectively.
For comparison, if the inner cavity of the protoplanetary disk is created by the magnetic fields of the protostar, the truncation radius is given by the balance between the stellar magnetic stress and the disk Reynolds stress (Starczewski et al. 2007;Chang et al. 2010): 10 -1 10 0 Stellar Mass (M ) 10 -2 10 -1 10 0 Semimajor Axis of Innermost Planet (AU) Figure 10.Orbital semimajor axis of innermost planet versus stellar mass for detected planetary systems with three or more planets, with red and green open circles for the planets discovered by transit and radial velocity, respectively, and blue filled circles for 11 confirmed and suspected resonant chain systems (see text for the list).The green dashed lines are the upper and lower estimates of the magnetic truncated radius a tr of the protoplanetary disk (Equation ( 12)), and the black dashed lines are the upper and lower estimates of 4a tr .The solid red line is a tr derived by Batygin et al. (2023).The solid blue line is the best-fit straight line to the blue filled circles.Data obtained from NASA Exoplanet Archive.
where η is a dimensionless factor of order unity, B * is the strength of the magnetic field at the stellar surface, R * and M * are the stellar radius and stellar mass, and Ṁdisk is the disk's accretion rate.If we adopt 0.5 < η < 1 for a star with an aligned dipole field, 500 G < B * < 1500 G, the stellar mass-radius relation R * /R ⊙ = 1.06(M * /M ⊙ ) 0.945 (Demircan & Kahraman 1991), and Ṁdisk = 10 −8 (M/M ⊙ ) 1.8 M ⊙ yr −1 (Muzerolle et al. 2003;Natta et al. 2006;Manara et al. 2016), we have a tr ∝ M 0.963 * and the estimated range of a tr is shown as the upper and lower green dashed lines in Figure 10.Since the zero-torque location is expected to be at some distance from the inner edge, we also show the estimated range of 4a tr as the black dashed lines.
As seen in Figure 10, the orbital semimajor axis of the innermost planet of the resonant chain systems (blue filled circles) shows a scaling relationship with the stellar mass that is roughly bounded by the two black dashed lines, which supports the idea that the resonant chains are assembled by halting the migration of the innermost planet at the zero-torque location near the inner disk edge.For the nonresonant chain multiplanet systems (open circles), many fall within or just outside the region bounded by the two black dashed lines, but some are significantly above the upper black dashed lines.This would be consistent with the former corresponding to systems where the innermost planet is stopped at the zero-torque location but the other planets in the system do not converge on it sufficiently to form a resonant chain and the latter corresponding to systems where the innermost planet stops migrating (due to, e.g., disk dispersal) before it reaches the zero-torque location.
There are significant uncertainties on how B * , R * , and Ṁdisk scale with M * .With different assumptions (B * ∝ M 1/3 * , R * ∝ M 1/3 * , and Ṁdisk ∝ M * ), Batygin et al. (2023) have recently argued that a tr ∝ M 1/3 * , which corresponds to an orbital period of ∼ 3.0 days independent of M * (red solid line in Figure 10).There is sufficient scatter in the resonant chain data shown in Figure 10 that they could be consistent with either a tr ∝ M 0.963 * or M 1/3 * .If we fit a straight line to the data, we find that a ∝ M µ * where µ = 0.71 ± 0.17 (solid blue line in Figure 10).In the future, as more and more resonant chain systems are discovered, if they continue to fall on the currently observed trend of orbital semimajor axis of the innermost planet versus stellar mass (a ∝ M µ * where µ ≈ 0.333-0.963),the case for a common origin near the inner disk edge would be strengthened.

CONCLUSIONS
Convergent migration is a necessary criterion to maintain planets in a resonant chain, if the resonant chain is formed before the dispersal of the protoplanetary gas disk and the planets continue to migrate due to planet-disk interactions.Therefore, we have investigated what kind of protoplanetary disk would allow such convergent migration.The planets are assumed to undergo type I migration in an adiabatic disk with power-law indices α = −d ln Σ/d ln a and β = −d ln T /d ln a.We have developed an analytic criterion to determine the convergent migration zone in the (α, β) parameter space.For a pair of planets in MMR, the convergent migration zone is bounded by two lines: the disk constraint (which is the boundary between inward and outward migration) and the planetary mass constraint (which is the boundary between faster inner planet migration and faster outer planet migration).The criterion was generalized to a resonant chain by requiring that any part of the resonant chain should be convergently migrating toward the remaining part.We have verified the analytic criterion with numerical simulations of a pair of 1M ⊕ planets in 3:2 MMR and the four-planet resonant chain in the Kepler-223 system.The analytic criterion was then applied to the Kepler-60, Kepler-80, TOI-178, and TRAPPIST-1 systems.In all cases, outward convergent migration requires rather extreme (mostly negative) values of (α, β), and there is little or no inward convergent migration zone.Even when we adjust the planetary masses within the uncertainties to maximize the inward convergent migration zone, there is little or no overlap between the inward convergent migration zone and common disk models (such as the steady-state, constant α ν disks) for all but one of the observed systems.
For TRAPPIST-1, the masses of the planets are well determined, and the biggest uncertainty in the existence and sizes of the convergent migration zone is whether the inner two planets b and c are (or were) in the resonant chain.For other resonant chains (in particular, Kepler-60), the uncertainties in the existence and sizes of the inward and outward convergent migration zones in (α, β) space can be improved by reducing the uncertainties in the planetary masses.In our analysis, we have assumed unsaturated nonlinear type I migration torque.The formulae developed by Paardekooper et al. (2011) can be used if the torques on the planets in a particular resonant chain system are significantly affected by the effects of viscous and thermal diffusion.However, the formulae depend on the viscous saturation parameter p ν and the thermal saturation parameter p χ , and it will not be possible to derive a simple analytic criterion for convergent migration.In fact, the criterion could vary as a function of position due to the functional form of p χ (see Section 2).
The innermost planets in the observed resonant chain systems are close to the stars, and an alternative scenario is that the resonant chain is formed and maintained by the inward migration of the outer planets toward the innermost planet after the migration of the latter is stalled near the inner disk edge.We have found support for this idea from a relationship between the orbital semimajor axis of the innermost planet and the stellar mass of the observed resonant chain systems.If resonant chains discovered in the future continue to fall on this relationship, the case for a common origin near the inner disk edge will be strengthened.

Figure 1 .
Figure 1.(a) Boundary between inward and outward migration according to Equation (6) for an adiabatic disk with power-law indices α for surface density and β for temperature.(b) Boundary between faster forced migration of inner planet and faster forced migration of outer planet according to Equations (7) and (8) for m 1 = m 2 .(c) Combined plot, with the convergent outward and inward migration regions corresponding to the regions bounded by the red and blue lines in Figure3.

Figure 3 .
Figure3.Numerical results in (α, β) space for the migration of a pair of 1M ⊕ planets in 3:2 MMR.Blue and red dots represent simulations where the MMR is maintained as the planets move inward and outward, respectively.Green dots represent the simulations where the MMR is broken.The blue and red solid lines are the analytic estimation of the convergent migration regions from Equations (6), (7), and (8).The black dashed line is α + β = 3/2 for the steady-state, constant α ν -viscosity disk models, with the star for the MMSN-like model and the other symbols for various regions of theGaraud & Lin (2007) disk model: pentagon for the weakly opaque region, square for the strongly opaque region, and triangle for the viscously heated region.

Figure 4 .
Figure 4. Analytic results on Kepler-223 from Equations (6), (10), and (11), with panels (a), (b), and (c) showing the regions in (α, β) space with convergent migration between planet b and triplet c-d-e, between pair b-c and pair d-e, and between triplet b-c-d and planet e, respectively, and panel (d) showing the combined result from the overlap of panels (a)-(c).Red indicates convergent outward migration, and blue indicates convergent inward migration.Panel (d) shows that the Kepler-223 resonant chain cannot be sustained by inward migration and that the convergent outward migration zone does not include the steady-state, constant α ν -viscosity disk models (dashed line).

Figure 6 .
Figure 6.Time evolution of semimajor axes (upper left panel), eccentricities (lower left panel), and twobody MMR and three-body angles (right panel) in a disk migration simulation of the Kepler-223 system with α = −2.586and β = −1.966.This (α, β) combination is just below the analytic convergent outward migration zone, but the resonant chain is maintained to the end of the simulation.

Figure 7 .
Figure7.Same as Figure6but for a simulation with α = 0.931 and β = 1.552.This (α, β) combination is near the line of no migration, and the resonant chain is not broken by the very slow inward migration after 1 Myr.

Figure 9 .
Figure 9. Analytic convergent migration zones in (α, β) similar to Figures 4(d) and 8 but for planetary masses adjusted within 1σ to enhance the inward convergent migration zone.

Table 1 .
Initial Conditions of Two Planets in 3:2 MMR for Disk-driven Migration Simulations • 197.60 • 296.6 • Note-The orbital elements are semimajor axis a, eccentricity e, inclination i, longitude of the ascending node Ω, argument of periapse ω, and mean anomaly M .The stellar mass is 1.125 M ⊙ .

Table 2 .
Initial Conditions of a Planetary System in 8:6:4:3 Resonance for Disk-driven Migration Simulations •Note-Stellar mass is 1.125 M ⊙ .