© The Authors 2025

1 Introduction

In the past two decades, with various transit and radial-velocity surveys, the number of multi-exoplanetary systems detected has exceeded 900¹. Among these, there is an excess of systems with adjacent planets in period ratios close to integer (Fabrycky et al. 2014; Steffen & Hwang 2015; Huang & Ormel 2023; Hamer & Schlaufman 2024; Dai et al. 2024), which are often identified with mean-motion resonances (MMRs). The MMRs in these multi-planetary systems are believed to form early, in the gasrich disk phase, through Type-I convergent migration (Goldreich & Tremaine 1979; Lin & Papaloizou 1979), with the exact resonant state related to the disk and planet properties (Ogihara & Kobayashi 2013; Kajtazi et al. 2023; Huang & Ormel 2023; Batygin & Petit 2023). Therefore, resonant systems preserve relics of formation processes and can unveil the footprint of the dynamical history of multi-planetary systems (Snellgrove et al. 2001; Papaloizou & Szuszkiewicz 2005; Teyssandier & Libert 2020; Huang & Ormel 2022; Pichierri et al. 2024).

The formal dynamical identification of a two-body resonance (2BR) is in terms of an angle (Murray & Dermott 1999), which relies on a precise knowledge of the argument of pericenter and is observationally hard to assess. However, it is much easier to recognize three planets in a three-body resonance (3BR) because the corresponding resonance angle only depends on the angular positions of the planets (their mean longitudes). Multiple exoplanetary systems have been confirmed to host planets in 3BRs, including GJ 876 (Rivera et al. 2010), Kepler-80. (Xie 2013; MacDonald et al. 2016, 2021), Kepler-221 (Rowe et al. 2014; Goldberg & Batygin 2021), Kepler-60 (Goździewski et al. 2016), Kepler-223 (Mills et al. 2016), TRAPPIST-1 (Gillon et al. 2017; Luger et al. 2017; Agol et al. 2021; Huang & Ormel 2022), K2-138 (Christiansen et al. 2018; Leleu et al. 2019, 2021; Cerioni & Beaugé 2023), HD-158259 (Hara et al. 2020), TOI178 (Leleu et al. 2021), TOI-1136 (Dai et al. 2023), HD-110067 (Luque et al. 2023; Lammers & Winn 2024), and Kepler-402 (see Section 6.4).

Among the planetary systems with a 3BR, Kepler-221 is unique because it hosts a non-adjacent 3BR containing a secondorder 2BR (Goldberg & Batygin 2021). Kepler-221 is a G5-type star (R_s = 0.82 R_⊙) hosting four planets (b, c, d, and e) at a distance of 385 pc, with the system information listed in Table 1 (Rowe et al. 2014; Berger et al. 2018). Kepler-221 system is likely to be relatively young, with its large lithium abundance suggesting that the system is younger than the Hyades with an age estimated around 650 Myr (Berger et al. 2018; Goldberg & Batygin 2021). In the Kepler-221 system, the radii of the planets are well constrained, but the masses of the planets are unknown due to weak TTVs (Berger et al. 2018; Goldberg & Batygin 2021). In the system, planets b, c, and e are in a 6:3:1 3BR with the period ratios relatively far away from the integer period ratio (the period ratio between planets b and c is 2.035, and the period ratio between c and e is 3.228 ). However, planet d is not in any resonance with other planets.

If the resonance chain in the Kepler- 221 planetary system formed due to convergent migration during the disk phase, it is generally expected that the planets will be locked into a firstorder resonance chain (Lee & Peale 2002). However, except for planets b and c, this is not what is observed in the Kepler-221 system. Two observations stand out. First, planets c and e are close to (or in) 3:1 resonance, which is weaker than a first-order resonance. Second, planet d is clearly not in resonance with other planets (not even close to any). These two features cannot be explained by a simple disk migration model.

To explain the unique dynamical configuration of the Kepler221 system, we propose a new dynamical model in which there were originally five planets in a first-order resonance chain, as would be expected from migration theory. The resonance breaks in the post-disk phase due to dynamical instability from infalling materials from a potential debris disk (Izidoro et al. 2017, 2021; Liu et al. 2022; Nagpal et al. 2024) or planetesimal scattering (Raymond et al. 2022; Ghosh & Chatterjee 2023; Griveaud et al. 2024; Wu et al. 2024). Therefore, the system became unstable, and two planets merged into the current planet d, while the three-body resonance between the other three planets (b, c, and e) reformed. Over long timescales, tidal dissipation fueled the expansion of these planets to their current period ratios. To investigate this scenario quantitatively, a modular approach was adopted, where we investigated the viability of each of these steps separately. The advantage of this approach is that the resulting conclusions are self-contained. In particular, when planet (b, c, e) 3BR expands toward the observed period ratios, we find that planet d puts additional constraints on the planet’s tidal parameters and masses, which are independent of the detailed formation scenario.

This paper is structured as follows. In Section 2, we present the dynamical configuration of the Kepler-221 system and the new four-phased dynamical model. In Section 3, we present a consistent simulation throughout all phases that successfully reproduced the observation configuration. Then we use N-body simulations to verify our model and conduct parameter studies following a chronological order. In Section 4, the evolution of the system up to the point of three-body reformation is discussed. Section 5 investigates the orbital expansion phase, which puts constraints on the mass ratio of planets in resonance. According to the simulation results, we assess the five-planet formation model in Section 6. Finally, we summarize our results and present them in Section 7.

2 Model

2.1 Resonance in the Kepler-221 planetary system

Planets b, c, and e in Kepler-221 are very likely in resonance (Goldberg & Batygin 2021). The most prominent feature of two planets in resonance is that their period ratio is close to integer ratios. The more dynamically accurate description of resonance amounts to following the positions of the planet conjunctions in the orbit. This is encapsulated in terms of the resonance angle.

The resonance angle can be expressed as: $${\varphi }_{1,2,X}=(j+o)l_1-j{l}_2-o{\varpi }_X,$$(1)

where l₁ and l₂ stand for the mean longitude for the inner and outer planet, ϖ_X stands for either the longitude of periapsis for planet 1 (X = 1) or planet 2 (X = 2), which, respectively, correspond to the inner and outer resonance angle. In Equation (1), j and o are integers called the resonance number and resonance order, respectively. For two-body resonances, first-order resonances are more stable (o = 1) than higher-order resonances (o > 1) (Murray & Dermott 1999). For two planets in resonance, either one or both inner and outer resonance angles should librate about a fixed value. For planets not in resonance, the resonance angle instead circulates from 0 to 2π.

In a system with three or more planets, we can calculate the three-body resonance (3BR) angle by combining outer and inner two-body resonance angles and eliminating the longitude of periapsis of the middle planet. The expression of the 3BR angle is: $${\varphi }_{3BR}=p{\lambda }_1-r{\lambda }_2+q{\lambda }_3,$$(2)

where λ₁, λ₂, and λ₃ are the mean longitude for the three planets. The quantity (p + q − r) is called the 3BR order. Planets are locked in 3BR if such an angle librates around a fixed value. The zeroth-order 3BR and first-order 3BR are generally stable and hard to break (Petit 2021). For the Kepler-221 system, the period ratio between planets b and c is 2.035, and the period ratio between c and e is 3.228, which is close to 2:1 and 3:1, respectively. For simplicity, a librating φ_3BR = pλ₁−r λ₂+q λ₃ is referred to by using the (p, q, r) 3BR notation in the following. Therefore, planets b, c, and e might be in the (2, 3, 5) zeroth-order 3BR. The B-value of a 3BR is defined by taking the time derivative of Equation (2): $$B={\dot{\varphi }}_{3BR}=p{n}_1-r{n}_2+q{n}_3,$$(3)

where n₁, n₂, and n₃ are the mean motion for the three planets and B is the B -value. For planets in 3BR, because their 3BR angle librates around a certain value, the B-value averaged over time should be close to zero.

The normalized B-value for a 3BR can be defined as (Fabrycky et al. 2014; Goldberg & Batygin 2021): $$B_{\text{norm}}=|B|/⟨n⟩,$$(4)

where ⟨n⟩ is the average of the mean motion of the planets. For planets b, c, and e of Kepler-221, B_norm = 1.03 × 10⁻⁵ when (p, q, r) = (2, 3, 5). Compared with other confirmed 3BR pairs such as K2-138 (Christiansen et al. 2018), Kepler-80 (MacDonald et al. 2016), and Kepler-60 (Goździewski et al. 2016), the normalized B-value for Kepler-221 is statistically closer to zero, which implies that planets b, c, and e are in (2, 3, 5) zeroth-order 3BR (Goldberg & Batygin 2021).

2.2 The formation of the Kepler-221 planetary system

To solve the puzzle of the formation of the (b, c, e) three-body resonance chain, we hypothesize that Kepler-221 originally harbored five planets, but two of them (d₁ and d₂) merged to form planet d. The model that we envision consists of four phases, as illustrated in Figure 1. In chronological order:

1. In the disk phase, five planets formed and then migrated inward due to Type-I migration. All planets are assumed to stay beyond the inner edge (Liu et al. 2017). This leads to convergent migration and the adjacent planets are locked into a stable first-order two-body resonance chain: 2:1, 3:2, 4:3, and 3:2. The disk then disperses and the stable resonance chain remains, as shown in panel b in Figure 1.

2. In the collision phase, planets d₁ and d₂ collide and merge into d. After the disk dispersal, instability first breaks the resonance between d₁ and d₂, which triggers system-wide instability that causes all resonances in the system to break (similar to what Raymond et al. 2022 investigate for the TRAPPIST-1 system). As planets d₁ and d₂ are the closest and lightest, they merge into one single planet d (Matsumoto et al. 2012).

3. In the third phase – the post-collision phase – the (b, c, e) 3BR reforms (panels d and e of Figure 1). Planets b, c, and e are initially not in resonance after the collision, but they stay very close to the zeroth-order 6:3:1 resonance position. Tidal damping in the system will damp the planets’ eccentricities. Convergent migration due to tidal damping results in the reformation of the b, c, and e 3BR.

4. After reforming the zeroth-order 3BR, the system expands to the current period ratio – the orbital expansion phase (see panel f of Figure 1). In this last phase, tidal damping will continue to expand the system while preserving the 3-body resonance (Goldberg & Batygin 2021). This is because tidal damping preserves the angular momentum but dissipates the energy. As a result, the period ratio between planets in resonance in the system expands and the system reaches the observed position at present.

In this scenario, the first three phases are essential to explain the presence of the 3-body (b, c, e) resonance chain in Kepler221. The last phase explains the departure of period ratios from exact integer ratios through (b, c, e) 3BR expansion. This phase is independent of the details of the previous three phases; it only requires that planets b, c, and e are in resonance.

2.3 Numerical implementation

To verify the above scenario, we conduct a suite of N-body simulations using the REBOUND package (Rein & Liu 2012). We use the WHFast integrator (Rein & Tamayo 2015) in the disk phase with a timestep of 0.2% of the orbital period of the innermost planet. In the collision, post-collision, and expansion phases, we use the Mercurius integrator (Rein et al. 2019). The dissipation forces are implemented using the REBOUNDx package (Tamayo et al. 2020).

In the disk phase, planets migrate inward and resonances naturally form. The disk tidal force on a certain planet i following (Papaloizou & Larwood 2000) is: $${\mathit{F}}_{\text{disk},\text{i}}=-\frac{{\mathbf{v}}_{\text{i}}}{2t_{\text{a},\text{i}}}-\frac{2({\mathbf{v}}_{\text{i}}·{\mathbf{r}}_{\text{i}}){\mathbf{v}}_{\text{i}}}{|r_{\text{i}}{|}^2{t}_{\text{e},\text{i}}},$$(5)

where t_{a, i} and t_{e, i} correspond to the semi-major axis damping timescale and eccentricity damping timescale on planet i, respectively. After the disk has dispersed, stellar tides damp the planet eccentricities. While conserving angular momentum, the eccentricity of planets is damped on a timescale of (Papaloizou et al. 2018): $$t_e=7.6\times {10}^5{Q}_{\mathrm{p}\mathrm{h}\mathrm{y}}\frac{m_i}{m_{\oplus }}{\left(\frac{m_{\odot }}{m_{\star }}\right)}^{1.5}{\left(\frac{r_{\oplus }}{R_i}\right)}^5{\left(\frac{a_i}{0.05\mathrm{a}\mathrm{u}}\right)}^{6.5}\mathrm{y}\mathrm{r},$$(6)

where m_i, a_i is the mass and semi-major axis of planet i and m_⋆ is the mass of the star. For simplicity, the simulations are conducted in the plane. In our simulations, we only apply eccentricity damping according to Equation (6). However, it must be understood that the effective Q_phy accounts for additional damping mechanisms, such as obliquity tides. Therefore, Q_phy should not necessarily be identified with the tidal quality factor of a body. For example, when obliquity tides operate the equivalent tidal damping parameter, accounting for both tidal dissipation and obliquity tides, is expressed as $Q_{\text{phy}}^{-1}=Q_{\text{d}}^{-1}+Q_{\text{o}}^{-1}$, where Q_o is related to the obliquity ϵ of the planet (see Section 6.3), Q_d = 3Q/2 k₂, Q is the tidal dissipation function of the planet, and k₂ is the Love number (Papaloizou et al. 2018).

In addition, once planets are locked in 3BR in Phase IV, we accelerate the simulation by 100 times faster by adopting a Q_sim in the simulation that is much smaller than the effective Q_phy:Q_sim = Q_phy/100. This will speed up the computations. A similar approach was used in previous literature (Papaloizou et al. 2018; Huang & Ormel 2022).

The masses (see Table 2) of the planets in the Kepler-221 system have not been directly inferred by observations due to weak TTV signals (Goldberg & Batygin 2021). Therefore, we work with three different mass models. In model 1 “mass-radius” it is assumed that the planets follow the mass-radius relationship derived by Chen & Kipping (2017). In model 2 “peas-in-a-pod” it is assumed that the planets are of equal mass (Weiss & Petigura 2020; Weiss et al. 2023). Finally, Model 2a and 2d are variations of Model 2 with a slightly more massive planet c and a twice as massive planet d.

3 Charting the course: a successful simulation

In this section, we prove the feasibility of the four-phase model by presenting a successful N-body simulation covering all four phases discussed in Section 2. The successful simulation is shown in Figure 2 with the initial mass of planets according to model M2a in Table 2, consistent with the mass constraint of the system (discussed in Section 5.2.2). These and other parameters are discussed in detail in the later section; here we focus on the narrative. The general setup of the simulations in this work is listed in Table 3.

The simulation starts with the disk phase as shown in panels a and e in Figure 2. In the disk phase, five planets (planets b, c, d₁, d₂, and e) are initialized a little beyond the first-order resonance locations and semi-major axis damping and eccentricity damping are added on all planets except b, as shown in panel a. Planet b is assumed to stay at the disk inner edge, where positive co-rotation torques prevent planets from migrating further inward (Paardekooper & Papaloizou 2009; Liu et al. 2017). This convergent migration locks the planets in a chain of firstorder resonances around 100 yr. The planets are quickly captured into two- and three-body resonances, including a 3BR between b, c, and e (panel e). The formation of resonances also excites the eccentricity of the planets (panel a). The eccentricities of planets after disk dispersal are determined by the disk damping parameters (See Section 4.1).

After disk dispersal (end of Phase I), the planets in the simulation stay close to the resonance position (integer period ratio), which is typical for young resonant planetary systems (Hamer & Schlaufman 2024; Dai et al. 2024). Due to the adopted damping parameters the planets end up with a relatively high eccentricity of around 0.01–0.1 (Teyssandier & Terquem 2014; Papaloizou et al. 2018; Yang & Li 2024).

In the collision phase, dynamical instability is assumed to have taken place shortly after disk dispersal. The situation is illustrated in panels b and f of Figure 2. After t = 20 yr, a kick is applied to planet d₂, which results in the breaking of the entire resonance chain. Because the breaking of the resonance chain allows for close encounters, the system becomes dynamically unstable which further excites the planet eccentricities. Because of the high (initial) eccentricities, orbit crossing follows rapidly (Zhou et al. 2007), around t = 100 yr after resonance breaking, and merging of planets d₁ and d₂ is achieved after 800 yr. Because of the short orbital instability phase, the other planets are still close to their resonant locations (see Section 4.1). In the simulation, we assume, for simplicity, that two planets will experience perfect merging if their distance becomes closer to the sum of their radii. In reality, even if planets d1 and d2 did not have a perfect merger and experienced a hit-and-run collision, the “runner” is expected to return for a second giant impact resulting in planet merging (Asphaug et al. 2021).

Panels c and g represent the post-collision phase in Figure 2. After the merging of d₁ and d₂, the eccentricity of planet d decreases because planets d₁ and d₂ merged at their respective perihelion and aphelion distances, leading to a more circular orbit for planet d after merging (Kokubo & Ida 1995). Specifically, the eccentricity of planet d₁ and d₂ is around 0.1 before the merger while the eccentricity of planet d is around 0.05 after merging. Tidal damping also continuously decreases the eccentricity. In the post-collision phase, planets b and c undergo convergent migration due to the outward torque that is applied to planet b (see Section 5.2.1). Such an outward torque on planet b is essential to reform the resonance chain (see Section 4.3). The 2:1 two-body resonance between planet b and c reforms first (after 1.7 Myr; the two-body angle is not shown in Figure 2). Thereafter, the c and e two-body 3:1 resonance and the b, c, and e 3BR reform simultaneously at around 5.7 Myr. The conditions are now in place for the planet period ratios P_c/P_b and P_e/P_c to expand along the zeroth-order 3BR line (blue-dashed line in Figure 2.

Panels d and h in Figure 2 show the evolution of the (b, c, e) and (c, d, e) period ratios during the orbital expansion phase (Phase IV). After the (b, c, e) 3BR reforms, tidal damping operates to expand the period ratio of the planets in resonance, migrating planets b and c inward and planet e outward. Initially, the planets are close to their corresponding 2:1 and 3:1 period ratios. The dashed lines in Figure 2 indicate different 3BRs characterized by p, q, and r in Equation (3): $$\frac{P_3}{P_2}={(\frac{r}{q}-\frac{p{P}_2}{q{P}_1})}^{-1},$$(7)

where P₁, P₂, P₃ correspond to the inner, middle, and outer planets (i.e., (b, c, e) in Figure 2d and (c, d, e) in Figure 2h). In panel d of Figure 2, tidal dissipation moves the planets along the blue dashed line corresponding to the (2, 3, 5) zeroth-order 3BR between b, c, and e (Charalambous et al. 2018; Papaloizou et al. 2018). At the same time, the period ratio of planets c, d, and e also evolves along a line that avoids interaction with certain (c, d, e) 3BRs as shown in Figure 2h. The slope of the expansion seen in Figure 2h is related to the masses of the planets, as is discussed in Section 5.1. After about 2 Gyr of expansion, the period ratios of the planets are consistent with their observed values, implying a successful reproduction of the dynamical configuration of the Kepler-221 system.

In this section, a particular simulation has been selected to demonstrate the four-phase model. The successful result in Figure 2 is not guaranteed and the success rate of each phase is discussed in detail in Section 4.4. Obviously, the successful outcome is somewhat tuned and parameter-dependent. Specifically, several key milestones can be identified. These include: the high eccentricity the planets acquire post-disk phase; the occurrence of a collision between planets d 1 and d 2; the reformation of the b/c two-body resonance, followed by the 3BR between (b, c, e); the ability for sustained orbital expansion of the (b, c, e) system toward the current period ratios; and the inability of planet to dislodge these planets out of the 3BR they are in today. However, the successful result in Figure 2 serves the purpose of charting the course in this section. In the following, we discuss in detail the physical conditions that must be satisfied to pass each successive milestone and quantify the overall success rate of each model step.

4 Early resonance formation, dynamical instability, and 3BR reformation (Phase I-III)

In the previous section, we introduced a consistent simulation throughout all phases which successfully reproduced the observed configuration. In this section, the evolution from Phase I to III until the (b, c, e) resonance reformation will be investigated. The masses of planets in the simulations follow Model M2a in Table 2 in this section. The focus of this section will be on the conditions for which planet d₁ and d₂ manage to merge and the b, c, and e 3BR to successfully reform after the assumed planetary collision in Phase II. The range of parameters for the simulation from Phase I to III are listed in Table 3.

4.1 Evolution until dynamical instability

In the disk phase, the five planets (planets b, c₁, d₁, d₂, and e) form a stable first-order resonance chain by disk migration with planets b, c, and e in 6:3:1 resonance, as shown in panels a and e in Figure 2. The mass of planet d is split between d₁ and d₂ such that its center of mass lies close to the observed position of planet d after the orbital expansion phase. Therefore the mass ratio of m_d1/m_d2 is assumed to be in the range of [1, 1.04]. The formation of the first-order resonance chain is guaranteed because the planets are initialized close to their corresponding resonance position. To test the stability of the resonance chain after resonance formation, we apply a disk dispersal timescale t_d in the range of [50 yr, 10⁴ yr]. The planets remain in resonance throughout disk dispersal. This is because the number of planets in the system is lower than the critical number in a firstorder resonance chain (Matsumoto & Ogihara 2020). Therefore, dynamical instability is required to break the resonance chain.

It has been suggested that dynamical instability (Izidoro et al. 2017) or planetesimal scattering (Raymond et al. 2022) could perturb the system stability and break the resonance chain. Here, we do not model the specifics of the resonance breaking process, but simply enforce it by perturbing planet d₂ with a “kick” just strong enough to break the d₁ and d₂ resonance. In the simulation, we apply a small instantaneous torque on planet d₂ that decreases its semi-major axis (Δa/a) in the range of [0.5%, 1.2%] within 10⁻³ yr. The simulations show a Δa/a as small as 0.7% breaks the entire resonance chain (Zhou et al. 2007; Pichierri & Morbidelli 2020) and results in the rapid merging of planets d₁ and d₂, as shown in panels b and f in Figure 2.

In the collision phase, a quick merging between planets d₁ and d₂ is preferred. In that case, planets (b, c, e) would not experience much perturbation of their semi-major axes and would stay close to the 3BR resonance position, which is beneficial for reformation (see below). On the other hand, if the collision phase is long, the planets would undergo a long time of instability that would cause the positions of planets (b, c, e) to deviate significantly from the exact integer period ratio. This is disadvantageous for the subsequent 3BR reformation process. How rapid d₁ and d₂ merge after resonance breaking is determined by the eccentricity of planets d₁ and d₂. If the initial eccentricity is high (for instance, around 0.1), orbit crossing and merging of planets d₁ and d₂ will quickly follow the breaking of resonance. On the other hand, if the eccentricity remains low, the orbit of planets d₁ and d₂ would only cross after the planet eccentricity excites and the time required for the planet merging will be much longer (Tamayo et al. 2021). Also, a higher eccentricity of planets in the resonance chains implies that they are deeper in resonance with the period ratio closer to integer (Huang & Ormel 2023). This is also advantageous for the 3BR reformation afterward with the initial position of the planet closer to the (b, c, e) 3BR position. The ratio between the semi-major axis damping timescale and eccentricity damping timescale t_{a, i}/t_{e, i} determines the eccentricity of the planets after disk dispersal (Teyssandier & Terquem 2014; Papaloizou et al. 2018) with the eccentricity being proportional to the square-root of the damping ratio parameter. For example, in 45% of the simulations planets d₁ and d₂ would merge within 10³ yr with t_{a, i}/t_{e, i} = 250 while the merging rate within 10³ yr is lower at 10% when t_{a, i}/t_{e, i} = 450. Therefore, in 70% of the simulations with t_{a, i}/t_{e, i} = 450, planets (b, c, e) exhibit a large deviation from 3BR with B_norm > 0.015 after the merging of planets d₁ and d₂ (see Equation (4)), making the 3BR reformation afterward almost impossible.

The success rate of mass model M2a and M2d (see Table 2 is tested in the simulation with the parameters in the range of Table 3. With mass model M2a, in 51% of the simulations planets d₁ and d₂ merge within 10⁴ yr. In the other 49% of the simulation, planets (b, c, e) already deviate a lot from the integer period ratio, making the future 3BR reformation impossible, so the simulation is truncated. On the other hand, the merging of other planets is possible when mass model M2d is assumed because of more massive planets d₁ and d₂. In 47% of the simulations still d₁ and d₂ merge within 10⁴ yr. In 6% of the simulations planets d₂ and e merge while in 11% of the simulation planets c and d₁ merge. The other simulations are truncated since planets (b, c, e) already deviate from the 3BR position.

4.2 Resonance reformation after dynamical Instability

Panels c and f of Figure 2 represent the post-collision phase after the merging of planets d₁ and d₂. While the success rate in Phase I and II is relatively high, the reformation of the (b, c, e) 3BR after dynamical instability is much more challenging. Figure 2 features in fact a successful condition in which the b, c, and e 3BR reforms because of an outward (positive) torque applied on planet b. This ensures that planets b and c convergently approach each other, facilitating 2BR formation (see Section 4.3). Otherwise, the reformation of the (2, 3, 5) zeroth-order 3BR cannot be accomplished (Petit 2021). Specifically, the desired scenario is that planets b and c first reform the 2:1 two-body resonance under convergent migration, which is relatively straightforward. Then b, c, and e 3BR form at the same time with the c and e in 3:1 resonance. However, in the scenario where Q_phy is the same for all planets, planets b and c tend to experience divergent migration. As a result, there is only a slight chance to reform the b, c, and e 3BR.

We run over 50 simulations with the setup applying the same Q_phy (see Section 2.3) on all planets but with slightly different initial positions and eccentricities of the planets for each simulation. In all these simulations, planets b, c, and e fail to reform the 3BR. An illustrative example is given by panels a and d in Figure 3. Instead of forming the desired 3BR, planet (b, c, e) either form a first-order 3BR ( 40%; as in Figure 3a) or the formation of the 3BR fails altogether ( 60% ). The reason why a first-order 3BR forms but not the zeroth-order is because the formation of a first-order 3BR is not preconditioned on the requirement of being close to exact resonance (see Appendix A) (Petit 2021). Therefore, in the absence of convergent motions, the planets would either fail to be locked in resonance or form the first-order 3BR and expand along it (red dashed line in Figure 3).

4.3 Potential mechanisms facilitating resonance reformation

In the previous paragraph, we outlined that convergent migration between planets b and c is essential to reform the two-body and three-body resonances between planets b, c, and e. Here, we envision two such mechanisms that can decrease P_c/P_b: i) a positive torque on planet b; ii) an enhanced tidal dissipation on planet c.

The first mechanism is a potential outward torque operating on planet b, which promotes the convergent migration between b and c and reforms the resonance. Panels b and e in Figure 3 show the 3BR reformation process corresponding to an exponentially decaying outward torque on planet b. The origin of the torque could include mass loss on planet b (Wang & Lin 2023) or dynamical tides (Ahuir et al. 2021). In the simulation, we take the form of the exponentially decaying torque: $$\mathrm{\Gamma }={\mathrm{\Gamma }}_0\exp (-\frac{t}{t_{\mathrm{\Gamma }}}),$$(8)

where t is the simulation time after the merging of planets d1 and d2, Γ₀ is the value of the torque at t = 0 and t_Γ is the timescale on which the torque decays. The initial value of the torque can equivalently be parameterized in terms of an orbital decay timescale, for instance, if it acts on planet b: $$t_{a0}^{\mathrm{b}}\equiv \frac{L_{\mathrm{b}}}{2{\mathrm{\Gamma }}_0},$$(9)

where L_b is the current total angular momentum of planet b. This is the form given in Table 3.

In panel c and g of Figure 2, the torque applied to planet b amounts to a 2% mass loss in 7.5 Myr during the post-collision phase, which corresponds to an initial torque of Γ₀ = 8 × 10²⁴N · m. Expressed in terms of the parameters defined above, t_{a0, III}^b = 2 × 10⁸ yr and t_{Γ, III}^b = 7.5 Myr. Mass loss of planet b would also excite the eccentricity of the planet on a timescale of t_e = 2m/γṁ (Correia et al. 2020). At the stated mass loss rates the corresponding eccentricity excitation timescale therefore amounts to a minimum eccentricity excitation timescale (with γ = 1) of t_e = 10⁸ yr (such timescale is even longer in Phase IV, see Section 5.2.1). Because tidal damping also operates on the planets with a timescale of 10⁶ yr (see Equation (6)), it is safe to neglect the excitation of eccentricity due to mass loss.

Panel b of Figure 3 presents the period ratio diagram of a simulation with the imposed torque operating on planet b. It can be seen that the planet period ratio first oscillates around the resonance position relatively far from the integer period ratio. Then the positive torque on b decreases P_c/P_b while quickly trapping planets b and c in two-body resonance. Planets b, c, and e also quickly get captured in 3BR, as can be seen from the librating resonance angle shown in panel e of Figure 3. Subsequently, the period ratio of the planets expands along the zeroth-order resonance line.

The second mechanism to promote resonance reformation is more efficient damping on planet c. Planet c could potentially have a relatively small Q_phy (with Q_{c, phy}/Q_{b, phy} = 0.1) due to efficient obliquity tides (Louden et al. 2023; Goldberg & Batygin 2021) while the effective tidal Q parameter for other planets is much larger (see Section 6.3). Tidal damping on planet c would result in its inward migration and contribute to the convergent migration between planets b and c, necessary to the formation of the b and c 2BR and the b, c, and e 3BR. Panel c and f in Figure 3 correspond to such a simulation. Here, damping on planet c is stronger, with a smaller tidal Q parameter for planet c,(Q_c/Q_b = 0.1). As shown in the figure, planets b, c, and e first oscillate around the zeroth-order resonance (black dashed line) with large amplitude. Then tidal damping on the system migrates planets b and c inward, decreasing P_c/P_b while increasing P_e/P_c. The subsequent convergent migration first forms the b and c 2BR. Afterward, tidal damping on planets b and c migrates c outward, leading to the convergent migration between planets c and e and the reformation of (b, c, e) 3BR.

The reformation of the (b, c, e) 3BR requires convergent migration between planets b and c. Otherwise, planets b, c, and e will be trapped in the first-order 3BR or will not get trapped in any 3BR. Here, we have demonstrated the viability of resonance reformation for two physically plausible scenarios of an outward torque on planet b and more effective tides on planet c. Nevertheless, the reformation of the (b, c, e) 3BR is stochastic and cannot be guaranteed. This is discussed in more detail in Section 4.4.

4.4 Success rate analysis

To understand the success of resonance reformation postcollision, we run 200 simulations with different parameters (see Table 3). In 101 of them planets d₁ and d₂ collide, as shown in Figure 4, while in the other simulations planets c and d₁ collide, or some planets are ejected out of the system. The overall success rate in the post-collision phase is around 10% ( 9 successful in 101 simulations). In the simulation, an outward torque is applied on planet b corresponding to a mass loss of 2% to promote the convergent migration between planets. Tidal damping is applied on all planets with the tidal Q_phy parameter the same on each planet and varies from 3 to 10. A simulation is categorized as a successful simulation if planet (b, c, e) successfully gets trapped in the zeroth-order 3BR.

The resonance reformation process is a critical step of the model. To obtain meaningful statistics on its viability we conduct a parameter study of Phases II and III. Altogether 200 simulations are run varying parameters such as tidal damping strength and additional torque on planet b. As it is already shown that convergent migration is critical, an outward torque is applied on planet b with t_{0, III}^b between [2.5, 7.5] Myr and t_{a, III}^b in the range of [6 × 10⁷, 2 × 10⁸] yr (see above, Section 4.3). Tidal damping is applied on all planets with the tidal Q_phy parameter the same for each planet and varying from 3 to 10. A simulation is categorized as successful if planet (b, c, e) successfully gets trapped in the zeroth-order 3BR. Based on these simulations, we identify milestones that must be taken to guarantee resonance reformation. The results are presented graphically in Figure 4 in the form of a decision tree. Although the bulk success rate is low, less than 10%, Figure 4 shows that it can be significant when certain conditions and Ansatzes are met.

Reformation of the 3BR is only possible when P_c/P_b > 2 and P_e/P_c > 3 immediately post-collision (this happens in 40% of the simulations). Here, the period ratio is measured as the average in the first 100 yr after merging of planets d₁ and d₂. Period ratios above exact resonance are needed in order to set up conditions for convergently migrating planets into resonance (see Section 4.2). In 68% of these simulations, planets b and c undergo a relatively slow convergent migration (−0.018 Myr⁻¹ < Δ(P_c/P_b)/Δ t < 0) and the b/c 2BR successfully reforms. Due to the outward torque on b and the tidal damping on the planets, planet c would migrate outward and undergo convergent migration with planet e after the formation of b/c 2BR. The trapping of planets (b, c, e) prefer slow migration of c (with −0.015 Myr⁻¹<Δ(P_e/P_c)/Δ t < 0, occurring in 78% of the simulation). Otherwise, if P_e/P_c decreases too rapidly, the period ratio between planets c and e will quickly drop below three, and the trapping of c/e 2BR and (b, c, e) 3BR becomes impossible. When all these migration criteria are satisfied, the (b, c, e) 3BR has a chance of forming in 43% of the simulation ( 9 out of 21 ).

In reality, the situation could belong in this category because we omitted the higher Q_phy parameters from our parameter study due to computational constraints. Specifically, we have enforced a maximum simulation time of 10 Myr, within which the resonance reformation must take place. Therefore, a more gentle convergent migration during the formation of b/c 2BR and (b, c, e) 3BR would naturally satisfy conditions 3) and 4) in Figure 4 with a higher tidal Q parameter (Wu 2005; Jackson et al. 2008). On the other hand, it is impossible to reform the (b, c, e) 3BR with Q_phy < 3 due to too rapid migration. Therefore, given the Ansatz that a collision took place, the (maximum) success rate of the (b, c, e) resonance reformation stands at about 17% (=0.4 × 0.43).

Finally, we run a similar success rate analysis for mass model M2d in Table 2 assuming planet d is about twice as massive as other planets. The success rate of 3BR reformation in Phase III is around 10%, similar to mass model M2a. This implies that a massive planet d, potentially due to the merging of planets d1 and d2, has minimal impact on the reformation of (b, c, e) 3BR in Phase III. In conclusion, the reformation of 3BR in the post-collision phase appears viable under plausible physical conditions.

5 Orbital expansion (Phase IV)

After the formation of (b, c, e) 3BR, discussed in Section 4, the period ratios between b, c, and e are still close to integer. On the other hand, planets b, c, and e in the Kepler-221 system are currently locked into three-body resonance, relatively far away from integer period ratio, with P_c/P_b = 2.035 and P_e/P_c = 3.228. This implies that Phase IV – the orbital expansion phase must have happened in the history of the Kepler-221 system (Goldberg & Batygin 2021), regardless of the formation origin of the 3BR. The 3BR expansion excluding planet d in the Kepler-221 system was already investigated in Goldberg & Batygin (2021). In this section, we run the expansion phase including planet d, and demonstrate under which conditions tidal damping acting on the b, c, and e three-body resonance is capable of expanding these planets to their current observed period ratios. We will show that planet d can significantly affect the orbital expansion phase due to (b, c, d) 3BR crossing, which could put constraints on the system parameters including the mass ratio of planets.

5.1 Orbital expansion: Failure and success

In the orbital expansion phase, a successful simulation matching the observed period ratio has two requirements: i) planets (b, c, e) should stay in resonance and expand to the current period ratio P_c/P_b = 2.035 and P_e/P_c = 3.228; ii) when planets (b, c, e) reach the current period ratio, planet d should also match the observed period ratio with P_d/P_c = 1.765 and P_e/P_d = 1.829. A successful expansion matching these two conditions, as in panels d and g in Figure 2, is not guaranteed.

Because the orbital expansion phase is not related to the formation process of (b, c, e) 3BR, we used an optimized model to form the (b, c, e) 3BR to simplify the parameter setup as the initial masses of the planets. We initialize five planets (b, c, d₁, d₂ and e) in the first-order resonance chain consistent with Phase I and adiabatically decrease the masses of d₁ and d₂ to zero to ensure that the remaining planet b, c and e are still in 3BR close to integer period ratio. Then we add planet d back. The detailed setup of the optimized model is discussed in Appendix B.

Figure 5 shows the outcome of two unsuccessful simulations with the same masses of planets following M1 in Table 2 and different assumptions on the initial value of P_d/P_c. Initially, the planets are close to their corresponding 2:1 and 3:1 period ratios, which is a condition for the trapping in the zeroth-order 3BR (see Section 4.2). The dashed lines in Figure 5 indicate different 3BRs characterized by p, q and r in Equation (7). Tidal dissipation moves planets along (2, 3, 5) zeroth-order 3BR between b, c, and e (Charalambous et al. 2018; Papaloizou et al. 2018), indicated in the left panels of Figure 5 by the blue dashed line.

While the (b, c, e) system starts from near-exact resonance locations, there is some freedom to choose the initial starting point of planet d. In Figure 5, panels a and b present a simulation where planet d is positioned below the (p, q, r) = (3, 5, 8) 3BR resonance line. In addition, only tidal damping operates on the planets and all planets have their Q_phys = 4.6, which is consistent with the successful simulation in Figure 2. These low Q-values are solely adopted for computational expediency since they speed up the simulation. However, as is discussed in more detail in Section 6.3 the young age demands a very efficient damping mechanism. The adopted value for Q_phys in the orbital expansion phase simulations does not affect the simulation outcome. A larger Q_phys will not change the simulation outcome, but it will make the duration of the orbital expansion to the observed period ratio longer. Therefore, the age of the system (younger than 650 Myr ) might be inconsistent with the simulation, as explained in Section 6.3. During tidal expansion the periods of planets b, c, and e change, following (by definition) the (2, 3, 5) 3BR line (panel a), while P_d remains largely constant. As a result, the increase of P_e and decrease of P_c lead to the increase of both P_e/P_d and P_d/P_c and the system moves away from the (3, 5, 8) resonance line (panel b). Therefore, during the orbital expansion planets c, d, and e would not encounter this 3BR, which could dislodge it from the (b, c, e) resonance (see below), and planets b, c, and e successfully expand to the observed period ratio (panel a). Nevertheless, this scenario fails as the present location of planet d with respect to c and e (the black cross in Figure 5b) cannot be reached.

The slope in the P_e/P_d versus P_d/P_c plane (Figure 5b) is the result of the conservation of the 3BR resonance angle and angular momentum. Specifically, without considering the potentially disturbing effects of planet d, the slope in this plane can be analytically solved (see Appendix C): $$a_{\mathrm{c}\mathrm{d}\mathrm{e}}=\frac{\mathrm{\Delta }\left(P_{\mathrm{e}}/P_{\mathrm{d}}\right)}{\mathrm{\Delta }\left(P_{\mathrm{d}}/P_{\mathrm{c}}\right)}=\frac{4.90m_{\mathrm{b}}+4.94m_{\mathrm{c}}}{7.27m_{\mathrm{e}}-m_{\mathrm{b}}}.$$(10)

From Equation (10), it is obvious that the slope of the expansion increases with a more massive planet b and c or a less massive e. In Figure 5b, the effect of the b, c, and e 3BR expansion on the evolution of P_d/P_c and P_e/P_d is represented by the black line, which has a slope of 0.89. The simulation including planet d, represented by the colored line in panel b in Figure 5, shows little difference from the analytical solution represented by the black line. The difference stems from a minor outward migration of planet d due to secular interactions with the other planets (see Appendix D). Therefore, the value of P_d/P_c in the simulation (colored line) is a little larger with respect to the analytical prediction (black line).

A possible solution to match the low a_cde is to decrease the initial semi-major axis of planet d, and thereby P_d/P_c, in such a way that it could reach the correct observed position. One example is shown in panels c and d of Figure 5. However, this simulation is also unsuccessful, because planets c, d, and e would encounter the (3, 5, 8) 3BR during the orbital expansion, as shown in panel d. The resonance overlapping experienced at this point would break the b, c, and e 3BR and stop the orbital expansion of these planets (Petit et al. 2020). We find that the (b, c, e) 3BR are dislodged in more than 95% simulations, similar to what happened in panel c.

Therefore, the system likely started its orbital expansion right of the (3, 5, 8) 3BR line in the P_e−P_d−P_c period plane. Reaching the observed period ratio is only possible with: $$a_{\text{cde}}>\mathrm{1.56.}$$(11)

Thus, the expansion of P_d/P_c must feature some degree of slowdown, while that of P_e/P_d must not. We envision two possible mechanisms to achieve this, which are shown in Figure 6.

1. A (temporary) slowdown in the expansion of P_d/P_c due to outward migration of planet c, for example through a positive torque acting on planet b. Because planets b, c, and e are linked by 3BR, the outward migration of b would also lead to the increase in the semi-major axis of c and e. An example of a scenario that includes positive torques is shown in panels a and b of Figure 6. Specifically, a positive torque, with t_{a0, IV}^b = 2.1 × 10¹⁰ yr and t_Γ^b = 1.4 Gyr, is added to planet b (see below, Section 5.2.1). The torque slows down the expansion of P_d/P_c with respect to P_e/P_d. Its decaying nature results in an S-shaped trajectory in the (c, d, e) period ratio plane (Figure 5b). As a result, the planets successfully reach the observed period ratio near the end of the simulation (see Section 5.2.1).

2. A combination of planet masses that results in sufficiently steep a_cde. A corresponding simulation is shown in panels c and d of Figure 6. As shown in Equation (10), increasing the slope in the P_e−P_d−P_c plane can be achieved with a more massive bor cor a less massive e. An example is mass model M2a in Table 2, which assumes a slightly more massive planet c with respect to the other planets but is still in line with the peas-in-a-pod mass model (Weiss & Petigura 2020; Weiss et al. 2023). The orbital expansion resulting from these planet masses is shown in panels e and f in Table 2. Due to the increased slope in the P_e/P_d−P_d/P_c plane, the initial position of planet d can be chosen to avoid the (3, 5, 8) 3BR. This scenario also ensures a successful orbital expansion to the observed period ratio.

Goldberg & Batygin (2021) expanded b, c, and e 3BR, but planet d was neglected because its current position is far from a 3BR with any other planet. We stress, however, that planet d is crucial in the orbital expansion phase. Although planet d is far from 3BR currently, when planets c, d, and e hit a 3BR during the orbital evolution, it could break the 3BR between b, c, and e and result in a different architecture from the observation. Therefore, it is necessary to include planet d in the orbital expansion simulation. The initial, relative positions of planet d and the (c, d, e) 3BR are crucial for the subsequent orbital expansion of the system.

5.2 Orbital expansion: Specific scenarios

In this subsection, we conduct parameter studies for the two mechanisms outlined above: (i) outward migration of the (b, c, e) system; or (ii) a change in the masses of the planets.

5.2.1 Outward-moving planet b

The first mechanism is the outward migration of planet b. If the masses of the planets follow M1 in Table 2, the planets could not reach the observed period ratio, as shown in Figure 5. However, the planets could be secured into the observed position if the outward migration of planet b slows down the increase of P_d/P_c. The origin of the torque could be due to multiple reasons, including mass loss (Carroll-Nellenback et al. 2017; Wang & Lin 2023; Vazan et al. 2024) or dynamical tides (Bolmont & Mathis 2016; Benbakoura et al. 2019; Ahuir et al. 2021). Importantly, these torques are time-dependent and operate only at early times.

For isotropic mass loss, the total angular momentum is conserved and the planet migrates to a higher orbit upon losing mass. This process effectively amounts to a torque operating on the planet. If the mass loss fraction is A and the exponential decaying mass loss timescale is τ_m, the corresponding torque on planet b according to Wang & Lin(2023), Γ_ML is $${\mathrm{\Gamma }}_{\mathrm{M}\mathrm{L}}=\frac{L_{\mathrm{b}}}{{\tau }_{\mathrm{m}}}\frac{A\exp (-t/{\tau }_{\mathrm{m}})}{1+A\exp (-t/{\tau }_{\mathrm{m}})}\approx {\mathrm{\Gamma }}_{\mathrm{M}\mathrm{L},0}\exp (-t/{\tau }_{\mathrm{m}}),$$(12)

where Γ_{ML, 0} = AL_b/τ_m, t is the simulation time after the 3BR reformation, L_b is the current total angular momentum of planet b. Such torque is equivalent to a semi-major axis expansion timescale t_{a0, IV}^b following Equation (8) with t_{Γ, IV}^b = τ_m. In the orbital expansion phase, if the masses of the planets follow model M1 in Table 2, the torque required to expand the planet to the observed position corresponds to t_{a0, IV}^b = 2.1 × 10¹⁰ yr and t_{Γ, IV}^b = 1.4 Gyr. Notice that t_{Γ, IV}^b is much longer than the corresponding timescale in Phase III (t_{Γ, III}^b, see Table 3) because the mass loss operates over a longer timescale in Phase IV.

On the other hand, dynamical tides are better described with a power-law decay model. We fit the results obtained by Ahuir et al. (2021), where a rapid decline of the torque is seen to occur after a time τ_d: $${\mathrm{\Gamma }}_{\mathrm{D}\mathrm{T}}={\mathrm{\Gamma }}_{\mathrm{D}\mathrm{T},0}{\left(\frac{t+{\tau }_{\mathrm{d}}}{{\tau }_{\mathrm{d}}}\right)}^{-2},$$(13)

where Γ_{DT, 0} is the initial dynamical torque and τ_d = t_Γ is the time after which the dynamical tides rapidly decay. For solar-likes star such as Kepler-221 we put τ_d = 460 Myr (Ahuir et al. 2021). Figure 7 plots these expressions for Γ_ML and Γ_DT by the solid green and blue lines, respectively. The torque due to the dynamical tide on planet b is represented by the blue solid line (Ahuir et al. 2021), which is much stronger than the equilibrium tides represented by the black line in Figure 7 (assuming Q_⋆′ = 10 ⁶).

We have conducted simulations with the above timedependent torque expressions acting on planet b, where we vary Γ_{ML, 0} and Γ_{DT, 0}. We find that if the masses of the planets follow M1 in Table 2, torque expressions that are approximately a factor 2 greater than the reference expressions could obtain the desired configuration where planet d ends up at the observed period ratio. These curves are shown with dashed lines in Figure 7. From the figure, it is clear that the torque required in the simulation is of the same order of magnitude as the torque generated by mass loss or dynamical tides. After 500 Myr, the increase in P_d/P_c slows down, allowing the system to execute an S-curve bend as illustrated in Figure 6 d. Finally, the planet period ratios increase to the observed values while avoiding the (3, 5, 8) 3BR line.

5.2.2 Mass constraint on b, c, and e

A more direct way to avoid crossing the (c, d, e) (3, 5, 8) 3BR is to assume different masses of the planets. The masses of planets b, c, and e determine the slope of expansion of P_d/P_c and P_e/P_d defined as a_cde in Equation (10). To reach the observed period ratio at the end of the orbital expansion phase, the slope should satisfy a_cde > 1.56 (see Section 5 and Appendix C). Therefore, if there are no additional mechanisms (for instance Section 5.2.1) acting on the planets, we can constrain the masses and densities of the planets with the current period ratio in the orbital expansion phase.

Figure 8 shows the constraints on the mass and density ratios in the orbital expansion phase in which the migration of planet d in the orbital expansion phase has been neglected. Successful orbital expansion is only possible if the combination of mass or density ratios falls in the green region, which corresponds to a_cde > 1.56. In general, this requires a more massive planet c than planet e (for instance, mass model M2a or M2d in Table 2, green dot in Figure 8). Because mass model M2a and M2d have the same mass for planets b, c, and e, only one dot is presented in Figure 8. Initial masses according to the mass-radius relationship in Chen & Kipping (2017) would result in too small a_cde for successful orbital expansion (M1 in Table 2, blue dot in Figure 8) and require additional mechanisms (for instance, stronger damping on planet d and outward torque on planet b, see Section 5).

The peas-in-a-pod mass model (M2 in Table 2; orange dot in Figure 8) assumes equal masses for all planets, resulting in a_cde = 1.56, just sufficient according to Equation (11). Yet, it would be difficult to reach the observed period ratio, because planet d would also migrate slightly outward due to secular interactions, resulting in a deviation from the observed period ratio during expansion (see Appendix D). Therefore, a modest adjustment of the masses – still within the peas-in-a-pod framework – is necessary. For example, mass model M2a or M2d with a 25% more massive planet c would ensure a successful orbital expansion phase.

6 Discussion

6.1 Assessment of the model

In this work, we have investigated several mechanisms to explain the observational constraints that the Kepler-221 system presents to us. A summary of the comments for each of these mechanisms is listed in Table 4 in chronological order. In Phase I and II, the merging of planets d₁ and d₂ is a natural result of the first-order resonance breaking after disk dispersal, and a quick merging within hundreds of years is achieved with high eccentricities for planets after disk dispersal (see Section 4.1).

Perhaps the most critical phase for the collision scenario is the re-establishment of the 6:3:1 nonadjacent 3BR between planets b, c, and e. We found that reformation is not spontaneous (Petit 2021), but instead requires either a stronger damping on planet c or a positive torque on planet b. The latter scenario is favored because the outward migration of planet b also promotes a successful subsequent orbital expansion phase, which ensures that planets c, d, and e can avoid 3BR encounters during Phase IV (see Section 5.2.1). Still, resonance reformation is not guaranteed, but stochastic, and successful only in ∼17% of our simulations (see Section 4.3).

It is likely that planets b, c, and e have experienced significant orbital expansion while locked in the 3BR configuration. If the system is young, as suggested by Berger et al. (2018), eccentricity tides alone are unlikely to be the sole driver responsible for the 3BR expansion. Therefore, additional damping mechanisms, such as obliquity tides (Goldberg & Batygin 2021), are required to accelerate the expansion. However, a too-strong dissipation would hinder the reformation of (b, c, e) 3BR in Phase III. This tension between Phase III (which favors weak damping) and Phase IV (which favors robust damping) is discussed below in Section 6.3.

Also, we found in this work that planet d plays an instrumental role in this expansion as it has the potential to kick the planets out of resonance (see Section 5.1). This outcome could be avoided by tuning the parameters for the outward migration of planet b, but a more straightforward explanation is that prolonged expansion was made possible by the planet’s mass ratios.

6.2 Mass constraints

Our modeling provides constraints on the yet undetermined masses of the Kepler-221 planets (Berger et al. 2018). The first clue comes from the formation scenario (Phase II-III), where it was postulated that planet d arises from the merger of two planets. The merging could have attributed to significant atmosphere loss but little mass loss (Ghosh et al. 2024; Dou et al. 2024). If so, it is expected that the density of planet d is higher, which would make planet d the most massive planet in the current system. The reformation of 3BR in the post-collision phase is still viable with planet d twice as massive as other planets (see Section 4.4). Furthermore, once the (b, c, e) resonance is established, the long-term evolution of the system (Phase IV) is mostly independent of planet d’s mass (see Appendix D). Therefore, a more massive planet d remains a possibility.

The second stronger clue comes from the orbital expansion phase. This clue is stronger as it is independent of the formation model. In the orbital expansion phase, a slope of a_cde > 1.56 in the (c, d, e) orbital plane (see Section 5.1) is required, which constrains the masses and densities of planets b, c, and e in the absence of additional dissipative forces or torques on planets (see Section 5.2.2). From this, similar masses with a slightly more massive planet c are expected, such as in model M2a, consistent with the peas-in-a-pod model (Weiss & Petigura 2020; Weiss et al. 2023). This model renders the densities of the planets progressively decrease with distance, ρ_b > ρ_c > ρ_e (see Figure 8). A possible origin of the similar masses and the decreasing density of the planets in the system could be due to photoevaporation. Inflated low-mass planets with orbits of a few tenths of an AU could experience a boil-off phase after disk dispersal, in which the planets lose their envelope and contract quickly with only a few of their original envelope left after the disk dispersal (Owen & Wu 2016). Afterwards, an X-ray and EUV-dominated photoevaporation could further erode the planet atmosphere with a timescale of approximately 100 Myr (Owen & Wu 2013; Lopez & Fortney 2013). One example of systems with photoevaporation-dominated planets is Kepler-36, which hosts two close-in planets (closer than 0.15 au ) with similar atmosphere fractions when they are born and the inner planet has lost almost all of its atmosphere (Carter et al. 2012; Lopez & Fortney 2013; Owen 2019). A similar photoevaporation process could have happened in the Kepler-221 system for planet b, resulting in a planet b of terrestrial density and puffier outer planets. An additional advantage of photoevaporation on planet b is that the associated rapid mass loss would contribute to the outward migration of planets b and c, which facilitates success in Phase III and Phase IV of our model (see Sects. 5.2.1 and 4.3).

Assuming model M2a of Table 2, Figure 9 plots the planets in the Kepler-221 system along with the planets of the exoplanet sample for which masses and radius are available². The plot shows that assuming similar masses following the peas-in-a-pod model (Weiss & Petigura 2020; Weiss et al. 2023) in M2a in Table 2 results in planet c and e puffier than Neptune. However, such an assumption is not unreasonable because similar super-puff planets exist in the exoplanet sample, especially in resonance chains (red background dots in Figure 9). For example, TOI-178 hosts six planets in a first-order resonance except for the innermost planets (Leleu et al. 2021). Among these planets, TOI-178 d and g are super-puff planets with radii similar to and density lower than the planets c and e estimated in M2a in Table 2 (Leleu et al. 2021; Delrez et al. 2023). Kepler-51 hosts three planets in a 3:2:1 resonance chain, and two of the planets are super-puff planets with masses similar to and radii much larger than the Kepler-221 estimation for planets c and e in mass models M2a. Therefore, the outer planets in the Kepler-221 system could potentially be another sample of super-puff planets in resonance.

It is worthwhile to conduct RV measurements on the Kepler-221 system with facilities such as Keck (Vogt et al. 1994). The significance of the RV detection is as follows: (i) if planet d turns out to be more massive than the other planets, it suggests a collisional origin. (ii) If the mass ratios of the planets b, c, and e fall in the green region in Figure 8, it offers the most natural scenario to explain the sustained orbital expansion phase (see Section 5). And (iii) it can provide evidence that planets c and e belong to the super-puff planets, while planet b is on the other side of the radius valley (Fulton et al. 2017) (see Figure 9).

6.3 Age constraint

Kepler-221 is a young system with an age below 650 Myr inferred from its large lithium abundance (Berger et al. 2018; Goldberg & Batygin 2021). Therefore, planets b, c, and e in 3BR need to undergo a sufficiently rapid orbital expansion to reach their present period ratios.

We define the time required to expand the b, c, and e 3BR from the integer period ratio to the observed period ratio as t_exp.. In our simulation, t_exp mainly depends on the tidal damping strength on the planets (determined by Q_phy) and is hardly affected by other parameters such as the mass ratios of the planets. For example, panels e and f in Figure 6 show a successful expansion with masses of the planets according to mass model M2a in Table 2 and the same Q_phy for all planets. From the orbital expansion phase simulations (see Section 5), the relation between t_exp and Q_phy can be expressed as: $$t_{\text{exp}}\approx 380Q_{\text{phy}}\text{Myr},$$(14)

which follows from the linearity of the damping rate with Q. Because the expansion duration time t_exp should not exceed the age of the system, we can conclude Q_phy < 2 corresponding to the maximum age of 650 Myr for the Kepler-221 system.

Such results conflict with the formation model of Kepler-221 in two aspects. Firstly, The value of Q_phy is over two orders of magnitudes smaller than the typical value of the tidal quality factor Q for terrestrial planets (Lee et al. 2013; Silburt & Rein 2015), indicating the existence of other much more efficient dissipation mechanisms to speed up the simulation. Secondly, the reformation process of 3BR requires gentle convergent migration between planets (b, c, e) corresponding to relatively large Q_phy (in Phase III we apply Q_phy > 3, see Table 3 and Section 4.4). This conflicts with the Q_phy < 2 constraint derived from the age of the system.

Goldberg & Batygin (2021) point out that the timescale puzzle could be solved with some additional dissipation mechanism required (for instance, obliquity tides; Millholland & Laughlin 2019). Obliquity tides occur when planets have a large axial tilt (obliquity) (Millholland & Laughlin 2019) and have been invoked to speed up the orbital expansion in the Kepler-221 system (Goldberg & Batygin 2021). A non-zero obliquity for planets can be maintained if the planet is in Cassini state where it reaches a secular spin-orbit resonance with synchronized precession of planetary spin and orbital angular momentum. The total energy dissipation rate for a planet with non-zero eccentricity and obliquity can be expressed as (Millholland & Laughlin 2019): $$\dot{E}(e,ϵ)=\frac{2K}{1+{\cos }^2ϵ}[{\sin }^2ϵ+e^2(7+16{\sin }^2ϵ)],$$(15)

where e is the eccentricity of the planet and ϵ is the obliquity of the planets. Here, 7Ke² is the tidal dissipation in the absence of obliquity (ϵ = 0), where K depends on the usual stellar and planet properties (Millholland & Laughlin 2019).

Planets are likely to have low obliquity during the disk phase as their angular momentum and spin vectors are aligned with the disk. But during the brief phase of dynamical instability and the resonance reformation afterward, planet b or c could have acquired a mutual inclination with respect to the other planets, which is a requirement for trapping in a Cassini state. Additionally, planet c dominates the tidal expansion process at the same level as planet b due to the large radius of c compared to b (see Equation (6)). Equation (15) demonstrates that even modest obliquity (for example, on planet c, consistent with Section 4.3) can dominate the energy dissipation, and drive the tidal expansion at rates much higher than what can be obtained from tidal damping alone.

However, the drawback of a small Q_phy is that such strong tidal damping would also lead to a rapid migration during the 3BR reformation phase. This makes the 3BR reformation unlikely (see Section 4.4). This implies either that the age of the system is actually older, or that some unknown mechanisms have sped up the orbital expansion in Phase IV. For example, Figure 2 is a successful simulation throughout all phases, which corresponds to a system age of around 2 Gyr. On the other hand, the orbital expansion configuration in Phase IV is unrelated to the expansion speed. Therefore, the constraints on the masses of the planets are solid regardless of the system age (see Section 5.2.2).

6.4 Application to other systems

Other exo-planetary systems with one planet outside the firstorder resonance chain may also be described with ideas introduced in this work. In particular, we highlight the K2-138 and Kepler-402 systems with the information of these systems listed in Tables 5 and 6. K2-138 is a K1-type star with six super-Earths. The period ratios between the inner five planets are close to 3:2 and all planets in K2-138 are in three-body resonance, with the innermost three pairs in (p, q, r) = (2,3,5) 3BR. The outermost three planets (planets e, f, and g ) are locked in a peculiar first-order (2, 3, 4) 3BR, which is the first pure first-order 3BR detected in a multi-planetary system (Cerioni & Beaugé 2023). Kepler-402 is an F2-type star with four super-Earths (planets b, c, d, and e). Similar to Kepler-221, planets b, c, and e are in a (7, 8, 15) 3BR, with the inner period ratio close to 3:2 and the outer period ratio close to 16:9. Planet d is not in resonance with any other planet.

The 3BRs in the K2-138 system are listed in Table 5. The inner three 3BRs pairs are close to (2, 3, 5) 3BR, with the period ratio increasing from inner pairs to outer pairs. Therefore, it is reasonable to assume that the inner five planets are originally locked into 3:2 two-body resonances with each adjacent triple also locked in 3BRs. At 1.544 and 3.290 the period ratios of planets e, f, and g are relatively far from integer period ratios, which could be the result of orbital expansion along (2, 3, 4) first-order 3BR from an inner 3:2 and outer 3:1 two-body resonance. The formation of first-order 3BRs was also seen in our simulations for the Kepler-221 system when we attempted to reform the (b, c, e) 3BR in the post-collision phase (see Section 4 and Appendix A). This demonstrates that the formation of pure first-order 3BR is possible in systems similar to Kepler-221. Compared to zeroth-order 3BRs, a first-order 3BR does not originate from two two-body resonances and could form relatively far from exact 2-body commensurabilities (Petit 2021) (see Appendix A). A possible formation scenario of the K2-138 could be that the inner five planets are first locked in a chain of 3:2 two-body resonances. Then tidal damping expands the period ratio of the inner planet, and the outer planet g coincidentally comes close to the (2, 3, 4) first-order 3BR and becomes trapped, similar to the resonance reformation process discussed in Section 4.3 and Appendix A.

The period ratios of the Kepler-402 planetary system are shown in Table 6 with the corresponding closest 3BR. It is clear from the figure that the adjacent planets are far from 3BR with larger normalized B -value but planets b, c, and e are in the (7, 8, 15) pure 3BR relatively far from integer period ratio, similar to the Kepler-221 system. This 3BR between planets b, c, and e might expand from an inner 3:2 and outer 16:9 two-body resonance. However, it is unlikely that such a high-order resonance can be formed with simple convergent migration in the disk phase. Also, the middle planet d makes the resonance formation between c and e harder (see Section 4). Therefore, the (7, 8, 15) 3BR for planets b, c, and e probably does not originate from two-body resonance between c and e. Possibly, the Kepler-402 system formed similarly to Kepler-221 in the sense that that there were originally two planets between planets c and e (d₁ and d₂) with all five planets in a first-order resonance chain. The resonance number for planets d₁ and d₂ two-body resonance is very high (7:6 or 8:7 resonance, facilitating the resonances to break quickly after disk dispersal (Matsumoto et al. 2012), and resulting in the merging of planets d₁ and d₂ into planet d. Afterward, tidal damping reforms the 3BR between planets b, c, and e, and the system period ratios expand to the observed value.

7 Conclusion

We have proposed a multi-phase formation model for the dynamical history of the Kepler-221 system, in line with its present architecture, of which the most peculiar feature is the out-ofresonance intermediate planet d. The envisioned scenario relies on two Ansatzes. First, there were originally five planets (planets b, c, d₁, d₂, and e) in a chain of first-order resonances. Second, after disk dispersal, the system experienced an instability, which broke all resonances and caused planets d₁ and d₂ to merge into planet d with its position consistent with its observed location. Under dissipative processes such as tidal damping, the 6:3:1 three-body resonance between planets b, c, and e reformed. On evolutionary timescales, the period ratio of these three planets expanded under tidal dissipation until they reached the observed period ratio.

To investigate the feasibility of the model we have carried out a detailed parameter study, in a modular fashion, where we detail the conditions necessary to overcome each milestone. The conclusions that emerge from this study are the following:

1. Immediately after the collision of d₁ and d₂, convergent migration between planets b and c is essential to reform the (b, c, e) 3BR. This implies a positive torque operating on planet b or an anomalously small tidal-Q of planet c. Stronger damping on planet c (possibly due to strong obliquity tides (Goldberg & Batygin 2021)) would migrate planet c inward while a positive torque on planet b would migrate planet b outward. If the appropriate conditions are satisfied, the success rate amounts to 17% (see Figure 4);

2. If these conditions are not present, resonance reformation would fail or the (b, c, e) planets would end up in a (1, 3, 3) first-order three-body resonance;

3. The properties of planet d are crucial for the outcome of the orbital expansion of the (b, c, e) sub-system. In order for the planets’ period ratios to evolve toward their current values, planets c, d, and e must have avoided the (3, 5, 8) resonance, which would break the b, c, and e 3BR with a probability of 95%;

4. Therefore, planets (b, c, e) should have evolved along a line of slope a_cde>1.56 in the period ratio plane of P_e/P_d vs P_d/P_c. Conservation of angular momentum and the (b, c, e) 3BR angle gives an analytical expression for a_cde (see Equation (10)), which only depends on the mass ratios of planets b, c, and e;

5. The condition a_cde>1.56 constrains the mass ratios of planets b, c, and e in the Kepler-221 system (see Section 6.2) to values consistent with the peas-in-a-pod scenario. Satisfying this mass constraint, a sustained expansion of the system is guaranteed. In addition, a positive torque on planet b could have increased a_cde to ensure successful orbital expansion when the mass constraint is not satisfied. This mechanism is particularly attractive as it also promotes resonance reformation;

6. Resonance reformation requires the planet to gentle converge onto each other at rates that correspond to effective damping parameter Q_phy > 3. On the other hand, the presumably young age of the system would require effective Q_phy < 2. This could imply that either the system actually is older, or that some unknown mechanisms have sped up the orbital expansion process in Phase IV;

7. The dynamical model proposed for the Kepler-221 system can be applied to other planetary systems in resonance, in particular, K2-138, which features three planets in a pure first-order 3BR, and Kepler-402, which, as with Kepler-221, features an intermediate planet not part of the resonance chain.

Data availability

The data underlying this article will be shared on reasonable requests to the corresponding author.

Acknowledgements

We thank the referee for their constructive comments, which significantly improved the presentation of this work. This work is supported by the National Natural Science Foundation of China under grant No. 12250610189 12233004 and 12473065. The authors gratefully thank Sharon Xuesong Wang, Fei Dai, Haochang Jiang, Yu Wang and Helong Huang for their useful discussions. Software: matplotlib (Hunter 2007; Caswell et al. 2021), REBOUND (Rein & Liu 2012), REBOUNDx (Tamayo et al. 2020).

Appendix A Trapping in first-order 3BR

In the post-collision phase, planets b, c, and e can escape from zeroth-order 3BR and instead get trapped in first-order 3BR, as shown in panels a and d in Figure 3. Although both originate from b/c 2:1 and c/e 3:1 2BR, the definitions of zeroth-order and first-order 3BR for planets b, c, and e are different. The expression for b/c 2:1 outer and c/e 3:1 inner 2BR angle are: $${\varphi }_{\mathrm{b}\mathrm{c},\mathrm{o}}={\lambda }_{\mathrm{b}}-2{\lambda }_{\mathrm{c}}+{\varpi }_{\mathrm{c}},$$(A.1) $${\varphi }_{\mathrm{c}\mathrm{e},\mathrm{i}}={\lambda }_{\mathrm{c}}-3{\lambda }_{\mathrm{e}}+2{\varpi }_{\mathrm{c}},$$(A.2)

where λ_i is for mean longitude for the planet i and ϖ_c represents the longitude of periapsis of planet c. For the definition of zeroth-order 3BR, ϖ_c is eliminated: $${\varphi }_{\mathrm{b}\mathrm{c}\mathrm{e},0}=2{\varphi }_{\mathrm{b}\mathrm{c}}-{\varphi }_{\mathrm{c}\mathrm{e}}=2{\lambda }_{\mathrm{b}}-5{\lambda }_{\mathrm{c}}+3{\lambda }_{\mathrm{e}},$$(A.3)

The definition of a first-order 3BR contains the argument of periapsis of one planet and requires that p + q − r = 1 (see Equation (2)). One example of first-order 3BR for planets b, c, and e is: $${\varphi }_{\mathrm{b}\mathrm{c}\mathrm{e},1}={\varphi }_{\mathrm{b}\mathrm{c}}-{\varphi }_{\mathrm{c}\mathrm{e}}={\lambda }_{\mathrm{b}}-3{\lambda }_{\mathrm{c}}+3{\lambda }_{\mathrm{e}}-{\varpi }_{\mathrm{c}},$$(A.4)

If planet b/c and c/e are in a 2BR resonance (i.e., φ_{bc, o} and φ_{ce, i} librate), both φ_{bce, 0} and φ_{bce, 1} librate because these two 3BR angles are linear combination of the 2BR angle. In this case, planets b, c, and e are trapped in (2, 3, 5) zeroth-order 3BR (see Equation (2)). If φ_{bee, 1} liberates, but not the zeroth-order 3BR angle and 2BR angles, then planets b, c, and e are trapped in (1, 3, 3) first-order 3BR. The strength of 0th and first-order 3BRs are roughly determined by p + q, with the smaller value representing stronger resonance (Petit 2021).

Generally, the formation of zeroth-order 3BR is an inevitable result of 2BR formation, so planets should first be captured when their period ratios are close to integer ratios. If planets are far away from integer ratio, then they are more likely to be captured into the stronger first-order 3BR with smaller p + q values for planets b, c, and e. Figure A.1 illustrates the process of resonance formation in the post-collision phase (see Section 4). The left panels show the simulation where no additional torque is applied on planet b, and the tidal Q parameter is the same for all planets. Planet b, c, and e skip the zeroth-order 3BR and get trapped into first-order 3BR. This is because planets b and c did not first reform 2BR because they do not undergo convergent migration (Petit 2021). So the planets period ratio expands skipping the zeroth-order resonance (black dashed line in panel a), and forms the stronger first-order 3BR along the red dashed line. Therefore, only the first-order resonance angle librates (see panels cand din Figure A.1).

The right panel in Figure A.1 corresponds to a simulation with an additional positive torque on planet b (see Section 4.3), which results in the outward migration of planet b. This decreases the period ratio between planets b and c (panel b) and leads to convergent migration between these two planets, and first reform b/c 2BR near integer ratio around 1.5 Myr (panel f). This promotes the reformation of the zeroth-order 3BR, which forms simultaneously with the c/e 2BR (panel d and f) around 1.7 Myr.

The simulations show that the formation of pure first-order 3BR is possible if planets are close to strong first-order 3BR position due to resonance breaking event. Similar scenarios could be applied to the K2-138 system, which hosts the only pure first-order 3BR resonance identified in the exoplanet census (Cerioni & Beaugé 2023). In the K2-138 system, the outermost three planets are trapped in a pure first-order 3BR. Possibly, the formation of K2-138 first-order 3BR is similar to the resonance formation process in Kepler-221, as shown in Figure A.1 (see Section 6.4).

Appendix B Optimized collision phase

In order to obtain the proper initial conditions for Phase IV in Section 5, we optimized Phase II and III, in which planets d₁ and d₂ are removed and planet d is added back adiabatically. Although artificial, this step ensures that planets b, c, and e remain in resonance and allow us to focus on the ensuing orbital expansion phase. We integrate the simulations of this orbital expansion phase until the point where the resonance breaks and the expansion process is the same as what we observe applying the initial condition directly from Phase III (see Figure 2), which proves the viability of this optimized model.

The optimized Phase II and III are advantageous to preserve the 3BR while merging planets d₁ and d₂. In this way, we arrive at a compact version of Kepler- 221 with planets b, c, and e in 3BR, and it can now be examined whether this configuration can evolve into the present configuration of Kepler-221.

Following the model introduced in Section 2, we proceed with simulation in chronological order with the optimized model, starting from the collision phase. Panels a and c in Figure B.1 show an example simulation of the optimized collision phase. We initialize the five planets in first-order resonance with masses according to Model 1 in Table 2 (Chen & Kipping 2017). In panel a, the (b, c, e) 3BR angle is first on a fixed value different from 0 or π because multiple 3BR are present for the five planets in the system. In the optimized collision phase, planets d₁ and d₂ are removed by decreasing their masses adiabatically with a timescale of 1 Myr until their masses reach zero. When d₁ and d₂ are removed, only (b, c, e) 3BR is left and it goes to the isolated fixed point librating around π (Delisle 2017). After this, planet d is inserted at its current position (panel a) and its mass adiabatically increases until m_d reaches the value shown in Table 2 with a timescale of 0.3 Myr, which ensures that the 3BR between planets b, c, and e do not break due to the mass change of planet d.

Appendix C Mass dependence in Phase IV

In the orbital expansion phase, the slope of expansion of P_d/P_c and P_e/P_d is crucial because it determines whether planet d can reach the observed position, as shown in panel b of Figure 5. In this section, we derive the analytical expression of the slope (Equation (10)) under the assumption that planet d is stationary.

We first consider the simplified model for the orbital expansion phase, as shown in Figure C.1. There are two constraints for the starting period ratio in Figure C.1. Firstly, the starting period ratio of P_d/P_c and P_e/P_d should be exterior to (to the right of) the (3, 5, 8) 3BR between planets c, d, and e (the blue dashed line); otherwise the 3BR between planets b, c, and e would break due to the encounter of the c, d, and e 3BR (as explained in Section 5). Secondly, the initial period ratio of c and e should satisfy P_e/P_c = 3 (right of the red dashed line in Figure C.1). Based on these two constraints, the lower limit of the slope of P_d/P_c and P_e/P_d for successful expansion to observed period ratio corresponds to the black solid line in Figure C.1, which starts at the intersection of the two constraint lines (black dot in Figure C.1) and ends at the observed period ratio, with a slope of 1.56. If we assume planet d does not migrate during the orbital expansion (P_d is a constant), then the slope can be expressed as: $$a_{\mathrm{c}\mathrm{d}\mathrm{e}}=\frac{\mathrm{\Delta }\left(P_{\mathrm{e}}/P_{\mathrm{d}}\right)}{\mathrm{\Delta }\left(P_{\mathrm{d}}/P_{\mathrm{c}}\right)}=-\frac{P_{\mathrm{c}}^2}{P_{\mathrm{d}}^2}\frac{\mathrm{\Delta }{P}_{\mathrm{e}}}{\mathrm{\Delta }{P}_{\mathrm{c}}},$$(C.1)

Therefore, the constraint of the slope for successful orbital expansion is a_cde > 1.56. Next, we analytically solve for ΔP_e/ΔP_c under the assumption of invariance of the resonance angle and angular momentum conservation.

When three planets with orbital period P₁, P₂ and P₃ are in a (p, q, r) 3BR the period ratio of the three planets satisfies: $$\frac{P_3}{P_2}=\frac{q{P}_1}{r{P}_1-p{P}_2},$$(C.2)

The total angular momentum for the three planets in a 3BR is conserved and can be expressed as: $$L_{\mathrm{t}\mathrm{o}\mathrm{t}}={\left(\frac{G^2{M}^2}{2\pi }\right)}^{1/3}(m_1{P}_1^{1/3}+m_2{P}_2^{1/3}+m_3{P}_3^{1/3}),$$(C.3)

Combining angular momentum conservation (Equation (C.3)) and the resonance constraint, we can eliminate P₁ to obtain: $$L_{\text{tot}}={\left(\frac{G^2{M}^2}{2\pi }\right)}^{1/3}[m_1{\left(\frac{p{P}_2{P}_3}{r{P}_3-q{P}_2}\right)}^{1/3}+m_2{P}_2^{1/3}+m_3{P}_3^{1/3}],$$(C.4)

where the total angular momentum L_tot is conserved. Implicitly differentiating this expression towards P₂ gives: $$\frac{m_2}{P_2^{2/3}}+\frac{m_3}{P_3^{2/3}}{Y}_{32}+\frac{p(r{P}_3^2-q{P}_2^2{Y}_{32})m_1}{{(q{P}_2-r{P}_3)}^2{\left(\frac{p{P}_2{P}_3}{r{P}_3-q{P}_2}\right)}^{2/3}}=0,$$(C.5)

where Y₃₂ = ΔP₃/ΔP₂ – the slope that we are seeking. Here we want to evaluate ΔP₃/ΔP₂ at a certain period ratio P₃/P₂ = P₃₂. Solving this equation for Y₃₂ we obtain the slope of the expansion: $$Y_{32}=\frac{p^{1/3}r{P}_{32}^2{m}_1+{(r{P}_{32}-q)}^{4/3}{P}_{32}^{2/3}{m}_2}{p^{1/3}q{m}_1-{(r{P}_{32}-q)}^{4/3}{m}_3},$$(C.6)

Similarly, by expressing P₃ in Equation (C.3) with P₁ and P₂ and defining P₂₁ = P₂/P₁, the slope in the period ratio between planets 1 and 2 can be obtained: $$Y_{21}=\frac{q^{1/3}p{P}_{21}^2{m}_3-{(r-p{P}_{21})}^{4/3}{P}_{21}^{2/3}{m}_2}{q^{1/3}r{m}_3+{(r-p{P}_{21})}^{4/3}{m}_2},$$(C.7)

where Y₂₁ = ΔP₂/ΔP₁.

For the Kepler-221 system, (p, q, r) = (2, 3, 5), ΔP_e/ΔP_c = Y₃₂, and P₃₂ = 3. Relabelling (1, 2, 3) → (b, c, e) and substituting (p, q, r) →(2, 3, 5) Equation (C.5) evaluates to: $$\frac{\mathrm{\Delta }{P}_{\mathrm{e}}}{\mathrm{\Delta }{P}_{\mathrm{c}}}=\frac{15m_{\mathrm{b}}+15.12m_{\mathrm{c}}}{m_{\mathrm{b}}-7.27m_{\mathrm{e}}},$$(C.8)

Therefore we obtain the dependence of a_cde on the mass of planets (also shown in Equation (10)): $$a_{\text{cde}}=\frac{4.90m_{\mathrm{b}}+4.94m_{\mathrm{c}}}{7.27m_{\mathrm{e}}-m_{\mathrm{b}}},$$(C.9)

Equation (C.9) shows the slope is steeper with more massive planets b and c and shallower with a more massive e. The slope in the simulation is consistent with this analytical expression although marginal differences arise due to the migration of planet d during the orbital expansion phase, which will be further discussed in Appendix D.

Appendix D Origin of the migration of planet d in Phase IV

In panel b of Figure 5, the evolution of P_d/P_c and P_e/P_d in the simulation (colored line) show minor deviation from the analytical evolution (black line) due to the outward migration of planet d. To explain such migration, this section investigates secular interaction between planets (Murray & Dermott 1999). Secular interaction contributes to the minimum value of the eccentricity of the planets without resonance interaction. In the orbital expansion phase, tidal dissipation expands the b, c, and e 3BR and conserves angular momentum, while planet d would migrate outward slightly due to angular momentum change. The relative tidal dissipation strength of resonant system expansion depends on the initial eccentricity of planets in resonance (i.e. b, c, and e).

Therefore, we estimate the eccentricity of planets c and e due to their secular interactions with planet d with the Laplace-Lagrange theory (Murray & Dermott 1999), in which the evolution of the eccentricity of the planets due to mutual secular interaction is analytically expressed by: $$h_j=e_{j1}\sin (g_{1}t+{\beta }_1)+e_{j2}\sin (g_{2}t+{\beta }_2)$$(D.1) $$k_j=e_{j1}\cos (g_{1}t+{\beta }_1)+e_{j2}\cos (g_{2}t+{\beta }_2)$$(D.2) $$e_j(t)={(h_j^2+k_j^2)}^{1/2}$$(D.3)

where j = 1, 2. e_j(t) shows the time evolution of the eccentricity of the two planets. Parameters in the formula including e_j1, e_{j 2}, g₁, g₂, β₁ and β₂ depend on the initial mass, semimajor axis, eccentricity, and orbital position of the two planets. We calculate the mean and standard deviation of the eccentricity of planets c and e due to the secular perturbation of planet d averaged over random initial conditions at the start of the orbital expansion phase, and compare it with the mean and average of the eccentricity in the simulation. The result is shown in Figure D.1, in which the blue line represents the simulation result and the orange line represents the analytical result. From the figure, it is clear that the two curves have similar values, indicating that the secular perturbation of planet d contributes most of the eccentricity of planets c and e. The difference between the simulation and analytical value is an outcome of secular and resonance perturbation from other planets. Also, the eccentricity for planets c and e is larger with a more massive planet d.

The dependence of the eccentricity of planets b, c, and e on the mass of d helps explain the evolution of planet d in the post-collision phase, as shown in Figure D.2, which shows the outward or inward migration speed of planet d with different initial conditions.

For the blue line in Figure D.2, planet d migrates neither outward nor inward because no tidal damping is applied on planet d. For the yellow line, tidal damping is added on planets b, c, and e but not on planet d, and planet d would migrate outward at a faster speed with a more massive planet d. This is because the tidal damping on planets b, c, and e would expand the 3BR, increasing the relative period ratio. During this process, tidal dissipation conserves angular momentum, planets b and c would migrate inward and planets d and e would migrate outward slightly due to the angular momentum exchange. A more massive d planet will increase the eccentricity of the other planets, increasing their relative tidal dissipation and the outward migration speed of planet d. Such effect of outward migration explains the difference between the black line and the colored line in panel b of Figure 5, in which the simulation (colored line) deviates from the analytical result (black line) because planet d also moves outward during the orbital expansion phase.

Tidal damping is only applied on planet d for the green line in Figure D.2, thus planet d migrates inward. This is because the eccentricity of planet d is already at the minimum value allowed by the secular interaction from other planets, so the tidal damping on planet d would not change its eccentricity but still migrate the planet inward. The tidal dissipation strength is independent of the mass of planet d, and less massive planets are affected more by the dissipation, leading to the mass dependence that the speed of inward migration is faster for the less massive planet d.

The red line in Figure D.2 corresponds to the case where tidal damping is added on all the planets. The value of the red line is similar to the addition of the yellow and blue line combining the two effects, with the less massive planet d migrating inward and the more massive planet d migrating outward.

References

---

All Tables

All Figures