© The Authors 2026

1 Introduction

Planet orbital migration driven by planet-disk interactions is a fundamental process shaping planetary system architectures. While migration in smooth disks has been extensively studied (Kley & Nelson 2012; Baruteau & Papaloizou 2013; Baruteau et al. 2014; Paardekooper et al. 2023), the dynamics near disk boundaries represent a distinct and crucial regime that is far less explored (Masset et al. 2006; Ogihara et al. 2010; Liu et al. 2017; Brasser et al. 2018; Liu et al. 2025).

When the gas surface density gradient at a cavity edge becomes sufficiently steep, the corotation torque can become strongly positive and dominate over the Lindblad torque, reversing the direction of inward migration (Masset et al. 2006; Liu et al. 2017; Romanova et al. 2026). This mechanism was identified in the context of magnetospheric cavities and termed “rebound” migration (Liu et al. 2017). Recent works have extended this concept to photoevaporative cavities, which open at larger radii (i.e., typically of a few astronomical units, au) and expand outward as the disk disperses (Liu et al. 2022). Hydrodynamical simulations have verified that a single planet near such an expanding edge can be pushed outward, with the rebound rate depending sensitively on planet mass and disk conditions (Liu et al. 2025, hereafter Paper I).

In multiplanet systems, the picture becomes more complex due to mutual gravitational interactions. Convergent migration often leads systems to be captured in mean-motion resonances (MMRs). However, the analytical work carried out by Liu et al. (2017) proposed that an expanding cavity edge can drive resonant pairs outward and potentially pull them apart. The observed prevalence of nearly resonant or fully nonresonant planet pairs in transit and radial velocity surveys (Fabrycky et al. 2014; Dai et al. 2024; Hu et al. 2025) indicates the need for mechanisms that alter the resonant configurations during and/or after late-stage gas disk evolution (Liu & Ormel 2017; Pichierri et al. 2024; Pan et al. 2025; Hansen et al. 2025); also see Papaloizou (2011); Baruteau & Papaloizou (2013); Izidoro et al. (2017); Wu et al. (2024); Lin et al. (2025, 2026).

Paper I provided the first hydrodynamic demonstration of a single planet near an expanding photoevaporative cavity edge undergoing outward rebound migration. In this work, we extend these investigations to multiplanet systems at the same type of disk edge. Using high-resolution simulations, here we study how planet pairs with varying mass ratios and initial orbital configurations evolve under the combined effects of disk torques and mutual interactions. We examine whether such planet pairs remain in or escape from the resonant configurations. Our results provide key inputs for future dedicated N-body studies aiming to explain the orbital architectures of observed exoplanet systems.

The paper is structured as follows. In Section 2, we introduce our numerical setup and model parameters. Sections 3 and 4 describe the migration outcomes for two-planet pairs and multiplanetary systems with various masses and initial orbital configurations. We present our discussion and conclusion in Sections 5 and 6, respectively.

2 Method

Following Paper I, we used the 2D grid-based hydrodynamical code Dusty FARGO-ADSG to investigate multiplanet migration in a photo-evaporating disk during its gas dispersal phase. Our simulations only model the gas content of the disk and the disk self-gravity is discarded. The governing equations are expressed as $$\begin{array}{rl}& \frac{\mathrm{\partial }\mathrm{\Sigma }}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\mathrm{\Sigma }\mathit{v})={\dot{\mathrm{\Sigma }}}_{\mathrm{P}\mathrm{E}},\\ & \frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+\mathit{v}\cdot \mathrm{\nabla }\mathit{v}=-\frac{\mathrm{\nabla }P}{\mathrm{\Sigma }}-\mathrm{\nabla }\mathrm{\Phi }+{\mathit{f}}_{\nu },\end{array}$$(1)

where Σ, P, v, and f_ν represent the gas surface density, pressure, velocity and viscous force per unit mass, and Φ is the total gravitational potential, which includes the direct contributions from the star and the planet, as well as the indirect contribution from the planet due to the acceleration of the coordinate frame (see Eq. (2) of Paper I). A locally isothermal equation of state was assumed. We adopted the mass-loss rate profile from Owen et al. (2012), where ${\dot{\mathrm{\Sigma }}}_{\mathrm{P}\mathrm{E}}$ is peaked in the inner disk region of ~2–3 au. This code allows us to self-consistently model the co-evolution of the planets and disk under the influence of turbulent viscous accretion and stellar X-ray photo-evaporation (PE). The code units for length, mass, and time are defined as r₀ = 10 au, M₀=M_⋆=1M_⊙, and $t_0={\mathrm{\Omega }}_0^{-1}=\sqrt{r_0^3/G{M}_0}$ (see Table 1), where Ω₀ is the Keplerian angular frequency at r₀ and G is the gravitational constant. We also defined P₀=2π/Ω₀, ≈31 yr as the canonical time unit for figures.

The gas turbulent viscosity is parameterized using the Shakura–Sunyaev α-prescription: ν=αc_sH, where α is set to 10⁻³, c_s is the sound speed, H=c_s/Ω_K is the disk’s pressure scale height, h=H/r is the disk aspect ratio, and ${\mathrm{\Omega }}_{\mathrm{K}}=\sqrt{G{M}_{\star }/r^3}$ is the Keplerian angular frequency at radial distance r. The initial disk surface density and aspect ratio are given by Σ₀(r)=Σ₀(r/r₀)^s and h(r)=h₀(r/r₀)^k, with default parameters listed in Table 1.

In this work, we chose a set of fiducial parameters that differs slightly from that of Paper I to reduce computational cost. While the parameters in Paper I were already computationally expensive for single-planet simulations, our new configuration, featuring a photo-evaporation rate 1.5 times higher and a slightly lower disk aspect ratio of h₀=0.05, keeps computational demands modest enough for simulating the rebound migration of multiple planets throughout the disk dispersal phase.

As in Paper I, grid cells are sampled logarithmically in radius between 0.2 r₀ and 10 r₀, and linearly in azimuth over the full 2π. Two resolution grids are used: a standard resolution (N_r=480, N_ϕ=770) for most cases and a high resolution (N_r=1200, N_ϕ=1900) specifically for the planet pairs with masses of 5 M_⊕ and 20 M_⊕. This choice of N_r ensures that the full width of the planet’s horseshoe region (where $\approx 2.4r_{\mathrm{p}}\sqrt{q_{\mathrm{p}}/h},q_{\mathrm{p}}$ is the planet-to-star mass ratio) can be resolved by approximately ten radial cells, which is desirable for accurate corotation torque calculations. The corresponding N_ϕ values yield approximately square grid cells for better numerical accuracy.

To optimize the computational efficiency, we first evolved the disk in 1D (N_ϕ = 1) without planets until an inner cavity began to form, defined as the gas surface density dropping below 10⁻⁵ g/cm². We then restarted the simulations in 2D, introducing planets at their full masses. The 1D simulations used standard open boundary conditions at both radial edges, while the 2D simulations employed wave-damping zones at these edges to prevent reflections of planet-induced wakes. We refer to Paper I for further numerical details. We considered five planet masses: 5 M_⊕, 10 M_⊕, 20 M_⊕, 100 M_⊕, and 318 M_⊕, representing super-Earths with two characteristic masses, Neptune-like planets, and gas giants similar to Saturn (0.3 M_J) and Jupiter (1 M_J). Various combinations of these five planet masses are adopted for the inner and outer planets.

We consider the planets are initially placed close to the first-order MMRs either the 3:2 or 2:1 resonances. These two states are most easily captured during long-range convergent disk migration (Quillen 2006; Kajtazi et al. 2023; Lin et al. 2025) and are also prevalent among observed multiplanet systems (Fabrycky et al. 2014). The innermost planet begins at r₀=10 au in our 2D simulations, while the outer planet’s orbit is placed slightly exterior to the exact resonant location (e.g., at 16.5 au for the 2:1 resonance).

We also calculated the disk torque exerted onto the planet Γ via $$\mathrm{\Gamma }={\int }_r{\int }_{\varphi }\mathrm{\Sigma }\frac{\mathrm{\partial }{\mathrm{\Phi }}_{\mathrm{p}}}{\mathrm{\partial }\varphi }rdrd\varphi =\sum _{r}T(r),$$(2)

where Φ_p is the planet’s gravitational potential, T(r) is the torque density that represents the ϕ-averaged torque across the radial grid cells.

3 Migration of two-planet systems

In this section, we present the results of simulations for two-planet pairs, varying both their masses and mass ordering. The simulation parameters are summarized in Table 2. We classify the results into three categories based on the planet masses. In all cases, both planets are inserted into the disk at t = 5 × 10³ P₀ when the inner cavity has just formed and the code is switched to 2D evolution. The specific evolution for each pair is detailed in the following subsections.

Fig. 1 Evolution of disk surface density, planet semimajor axes, period ratios, and resonant angles for a two-planet system consisting of 5 M⊕ and 20 M⊕ planets. The top (bottom) row shows the configuration in which the inner planet is less (more) massive. In the second column, the semimajor axes of the inner and outer planets are given by the red and blue curves, respectively, with shaded regions representing their perihelion and aphelion distances. The black dashed line marks the inner edge of the expanding cavity, defined as the radius where the surface density falls below 0.1 g/cm2. The fourth column plots the resonant angles 2λ1 − 3λ2 + ϖ2 (blue) and ϖ1 − ϖ2 (red), where the subscripts 1 and 2 denote the inner and outer planets, respectively.

3.1 Low-mass planet pairs

The top panels of Figure 1 illustrate the evolution of the disk surface density and the semimajor axes of planets for a system comprising an inner 5 M_⊕ planet and an outer 20 M_⊕ planet (fid_M02_R1.5 in Table 2). Initially positioned with a period ratio slightly above 1.5, the planets migrate inward rapidly and are captured into a 3:2 MMR within a few hundred P₀ (Figures 1b and 1c). As the inner planet approaches the cavity edge at t ≃ 1.5 × 10⁴ P₀, its inward migration halts as it enters the low-density cavity. Subsequently, as the cavity edge sweeps across the outer planet’s orbit, this inner planet reverses its migration direction, moving outward between t ≃ 2 × 10⁴ P₀ and 4 × 10⁴ P₀. This divergent migration forces the pair out of their initial resonant configuration, as confirmed by the resonant angle of 2λ₁ − 3λ₂ + ϖ₂ (λ and ϖ are mean longitude and longitude of periapsis, and the subscripts 1 and 2 denote the inner and outer planet) in Figure 1d. Ultimately, it librates until t ≃ 2 × 10⁴ P₀ before transitioning to circulation. The outer planet eventually stalls at 16 au, resulting in a final period ratio of 2.85.

The bottom panels of Figure 1 show the case with the mass ordering of the two planets reversed (fid_M20_R1.5). Here, the less massive outer planet migrates more slowly than its inner counterpart. The inner planet first migrates inward to the cavity and rebounds as the cavity expands. This evolution causes their period ratio to increase within the first 10⁴ P₀ before declining afterward. The two planets migrate convergently and are temporarily captured into the 3:2 MMR at 2 × 10⁴ P₀ (Figure 1g and h). The pair leaves the resonance as the cavity sweeps further outward, causing both planets to sequentially drop into the cavity. However, because the outer planet is less massive, its rebound migration is weaker when the cavity edge sweeps by. Consequently, the planets settle with a final period ratio of 1.7, only slightly higher than the original resonant value.

Figure 2 shows the snapshots of the surface density log(Σ) (left), and radial profiles of the individual planet’s torque density and total torque of the system (right) in fid_M02_R1.5 and fid_M20_R1.5. In the right panels, the red and blue solid lines show the torque density of the inner and outer planets, respectively, while the black dashed line represents their total torque (Γ_in + Γ_out). The planet pair is close to cavity edge in the top row of Figure 2 at t=1.5 × 10⁴ P₀. In this case, the outer 20 M_⊕ planet features an asymmetric torque density distribution, with the inner peak 10% higher than the outer peak. The torque onto the inner planet is very modest, since it is less massive and very close to the cavity where the density is strongly depleted. Thus, the net torque of the system is positive and driven by the outer planet (Figure 2b). This divergent migration causes the planet pair to move out of the resonance (Figure 1b). As the cavity further retreats, the inner planet does not migrate outward fast enough, leaving only the outer planet to rebound. This can also be seen in the middle row of Figure 2. The outer planet experiences a stronger, asymmetric positive torque as the cavity approaches it at t=4 × 10⁴ P₀.

On the other hand, in the case with an inner massive planet at t=1.5 × 10⁴ P₀ (bottom row of Figure 2), the inner massive planet located at the disk edge experiences a much stronger positive torque. For the outer low-mass planet, its torque contributions from the inner and outer peaks nearly cancel out. The total torque for the system is positive. In this circumstance the migration is convergent, which differs from that in the top row of Figure 2.

It is important to note that this work adopts a broad definition of rebound migration, encompassing any scenario where a planet migrates outward in tandem with an expanding disk cavity. Previous studies have shown that for low-mass planets, this outward migration is dominated by the corotation torque, driven either by a positive density gradient near the boundary (Masset et al. 2006), or a sharp boundary transition when only upper U-turn is present within the planetary horseshoe region (Liu et al. 2017). In the former case, the corotation torque is twosided and driven by a positive vortensity gradient, whereas in the latter limit, the angular momentum exchange in the corotation region comes solely from the upper horseshoe. Alternatively, for gap-opening planets, the required positive torque can originate from asymmetric Lindblad torques induced by disk eccentricity (Dempsey et al. 2021; Liu et al. 2025). The aim of this work is not to specify which underlying mechanism drives outward migration in each specific case, but highlight this resonance breaking channel through rebound.

Fig. 2 Snapshots from the data in Figure 1 of surface density log(Σ) (left), and radial distribution of the individual planet’s torque density T(r) and their total torque (right). Top two rows show results for fid_M02_R1.5 at t=1.5 × 104 P0 and 4 × 104 P0. Bottom row gives the same for fid_M20_R1.5 at t=1.5 × 104 P0. In the right column, the red and blue solid lines are the torque densities of inner and outer planets, whereas the black dashed line represents the total net torque summed over the two planets. The torque profiles are averaged over 20 snapshots within 1 P0. In all panels, red and blue circles indicate the positions of the inner and outer planets, respectively. The video is available online and on https://github.com/bbliu-astro/movies/blob/main/hydrodynamic_rebound/fiducial32_p0p2.gif.

Fig. 3 Same as Figure 1, but for a two-planet system consisting of 20 M⊕ and 0.3 MJ planets.

3.2 Low-mass and massive planet pairs

We also performed simulations with a low-mass planet of 20 M_⊕ and a giant planet of 0.3 M_J. The top panels of Figure 3 show the results from fid_M23_R2.0, where the inner planet mass is significantly lower than that of the outer planet. Although the pair starts with a period ratio slightly beyond 2, the outer massive planet opens a partial gap and undergoes rapid type III migration. Consequently, the planets easily bypass the 2:1 MMR and become trapped in the 3:2 MMR. Both planets migrate inward and fall into the disk cavity at t=9 × 10³ P₀ and 1.1 × 10⁴ P₀, respectively.

Their resonant state survives through the inside-out disk clearing phase. Notably, as the disk dissipates, the resonant state transitions from stable to overstable (Goldreich & Schlichting 2014; Lin et al. 2025). This is evidenced in Figure 3d by the resonant angle switching from a low amplitude libration around 180° to an increased amplitude of 60°, accompanied by an elevated eccentricity of the inner super-Earth planet from <0.01 to 0.12 (Figure 3b).

We observe an entirely different evolutionary pattern when the mass ordering of the planet pair is reversed (bottom panels of Figure 3). The planets initially migrate divergently until the inner giant planet enters the cavity. The outer super-Earth continues inward migration and is temporarily trapped in a 2:1 resonance. However, as the cavity edge sweeps outward, it pushes the super-Earth planet away, leaving the final system substantially displaced from the resonance.

The torque density and cumulative distributions are further analyzed in Figure 4. In the configuration with a less massive inner planet (top row of Figure 4), at t=9 × 10³ P₀, the inner planet lies within the cavity; whereas the outer planet opens a partial gap that overlaps with the gas clearing, driven by photo-evaporation. Under these conditions, the outer planet experiences a overall negative torque, causing both planets to eventually migrate into the inner disk cavity.

Conversely, in the case with a more massive inner planet at t=9 × 10³ P₀, the inner planet experiences a net positive torque at the disk edge; whereas the outer planet experiences a nearly zero net torque farther from the edge (middle row of Figure 4). As can be seen from Figures 3e and f, the cavity quickly sweeps through the inner planet’s orbit from t=8 × 10³ P₀ to 10⁴ P₀. Although the total torque for the pair is positive, the migration speed induced by the net torque in this case is lower than the cavity expansion rate. Consequently, the inner planet cannot keep pace with the outward-moving cavity and falls into it. At later times the disk only carries the outer, less massive planet, which continues to surf the expanding cavity until t~4 × 10⁴ P₀. For example, at t=3.5 × 10⁴ P₀, the positive torque peak on the outer planet is nearly twice as strong as the negative torque peak (bottom row of Figure 4).

3.3 Massive planet pairs

In the simulations with an inner planet of 0.3 M_J and an outer planet of 1 M_J, the pair is quickly captured into the 3:2 MMR. Their period ratio remains locked at 1.5 while their eccentricities are excited to ~0.1–0.2 during disk dispersal (Figures 5b and d). The pair remains in the resonant state but with a highamplitude libration of the resonant angles, which is a signature of overstability (Lin et al. 2025). This fragile nature, associated with large-amplitude libration and elevated eccentricities, indicates that the system might be unstable in the long term (Hu et al. 2025). However, as shown in Figure 5f, when the two planet masses are reversed, both planets undergo outward migration from the beginning. They reside in the 5:3 MMR until t ~ 2.5 × 10⁴ P₀. As the disk edge continues to retreat, the inner Jupiter-mass planet enters the cavity, leaving only the outer Saturn-mass planet to migrate outward and further increasing their period ratio (Figure 5g).

In contrast to the low-mass cases discussed above, the presence of two massive planets induces strong perturbations in the disk gas, as seen in Figure 6. Following the method of paper I, we computed the azimuthally averaged radial eccentricity profile of the disk. Treating the gas as pressure-less particles on Keplerian orbits around the central star¹, we first calculated the eccentricity for each gas cell from its local position and velocity. The eccentricity of the gas ring at radius, r, was then obtained by azimuthally averaging these individual cell values. Figure 7 illustrates the disk eccentricity for various planet pair configurations. Notably, the disk becomes significantly eccentric with e_g~0.15 for a Saturn-Jupiter pair in fid_M34_R1.5 (purple dots in Figure 7). Consequently, the resulting torque density profile departs from the classical two-lobed peak structure, instead exhibiting oscillations driven by large-amplitude non-Keplerian gas motions (top row of Figure 6). At t=7 × 10³ P₀, the torque originates primarily from the outer disk, yielding a net negative Lindblad torque that drives both planets inward.

When the Jupiter-Saturn pair are in nearly circular orbits (red dots in Figure 7), the principal term contributes to the disk eccentricity evolution (e.g., see Eq. (13) of Teyssandier & Ogilvie 2016 with l=m). For massive planets, gas in the horseshoe region is significantly depleted, eliminating effects from the corotation resonance. Generally, the outer Lindblad resonance excites eccentricity, while the inner Lindblad resonances damp it down (Goldreich & Sari 2003; Masset 2008; Teyssandier & Ogilvie 2016). Since the inner disk is quickly removed through photoevaporation, only the outer disk remains, which causes disk eccentricity to increase (purple line in Figure 7).

When the inner planet is more massive, however, the inner disk is not fully depleted before t=10⁴ P₀ (Figure 5e) and the disk still remains largely circular (red line in Figure 7). The two planets open a common gap, and the torques onto the two planets within the gap region nearly cancel out. The net torque on the system is thus mainly determined by the disk gas on both sides of the planets. Because the inner planet is more massive than the outer one, the inner residual disk exerts a stronger positive integrated torque onto the Jupiter-mass planet that overcomes the negative torque from the outer disk onto the Saturn-mass planet (see black dashed line in Fig. 6d). As a result, the system experiences a net positive torque, driving the outward migration of both planets. This migration pattern was first identified by Masset & Snellgrove (2001). We confirm that the similar physics also operates in the photo-evaporating disk with an expanding cavity. These two planets migrate outward together up to t=4 × 10⁴ P₀. Afterward, the Jupiter-mass planet falls into the disk cavity, while the Saturn-mass planet retains a net positive torque (bottom row of Fig. 6).

Fig. 4 Snapshots from the data in Figure 3 of surface density log(Σ) on the left and radial distribution of the individual planet’s torque density T(r) and their total torque on the right. Top row shows results for fid_M23_R2.0 at t=9 × 103 P0 and the bottom two rows give the same for fid_M32_R2.0 at 9 × 103 P0 and 3.5 × 104 P0, respectively.

Fig. 5 Same as Figure 1, but for a two-planet system consisting of 0.3 MJ and 1 MJ planets.

Fig. 6 Snapshots from the data in Figure 5 of surface density log(Σ) on the left and radial distribution of the individual planet’s torque density T(r) and their total torque on the right. Top row shows results for fid_M34_R1.5 at t=7 × 103 P0 and the bottom two rows give those for fid_M43_R1.5 at 7 × 103 P0 and 4 × 104 P0, respectively.

Fig. 7 Azimuthally averaged radial distribution of the disk eccentricity. Filled circles show the planet’s location in each simulation. Dashed part of the curves indicate the gas-depleted cavity region, where the calculation of the gas eccentricity is strongly affected by the reset in the gas density. The Saturn-Jupiter planet pair excites a large gas eccentricity, in contrast to the other three cases.

4 Migration of multiple planets

In this section, we discuss the simulations for three-planet systems. The masses of the planets from the innermost to the outermost are 10 M_⊕, 20 M_⊕, and 0.3 M_J. We explore their migration with two gas disk configuration: a fiducial model (Σ₀=44 g/cm²) and a model with lower gas density (Σ₀=13 g/cm²). Initially, all neighboring planet pairs are placed close to their 3:2 MMRs.

The top panels of Figure 8 illustrate the evolution in the fiducial disk model. The outer two planets are captured into a 3:2 resonance and migrate rapidly inward. This leads to a close encounter that scatters the innermost planet into the outer disk region even in the gas-rich disk phase (also see Marzari et al. 2010). The inner two planets of 20 M_⊕ and 0.3 M_J stay locked in the resonance as the disk cavity retreats outward. When the cavity approaches the outermost planet, it rebounds farther out. As seen in this proof-of-concept simulation, the orbital architecture is altered during disk dispersal, and the system evolves from initial compact configuration into final one consisting of a resonant pair and an isolated low-mass planet with wider orbital separations.

The bottom panels of Figure 8 show results for the same planet configuration but in a lower mass disk. In this case, disk dispersal proceeds much faster. Convergent migration results in planets sequential trapping into 3:2 resonances. Because the cavity expands too rapidly, no rebound occurs during gas dissipation, leaving the resonant chain intact. Combined with the analysis of Paper I, these results demonstrate once again that rebound migration requires specific conditions. This does not occur if the cavity expands either too quickly or too slowly, where this rate depends on the gas disk mass and/or the strength of stellar photoevaporation (see Figure 7 of Paper I).

We summarize all 20 simulations in Table 2 for systems with both two or three planets. The left panel of Figure 9 presents simulated systems with their masses and period ratios, while the right panel shows a selection of observed exoplanet systems for comparison. Since our simulations are focused on solar-mass stellar hosts and the photoevaporation mainly alters the planets whose orbits are on the order of 1–10 au, we selected the observed systems around solar-like stars where the innermost planet has an orbital period longer than 100 days. From inner to outer, planets are color-coded orange, blue, and purple, respectively, with the symbol size corresponding to planet mass. The solid horizontal line marks systems that remain in resonance, whereas the dashed line indicates period ratios that have deviated from resonance. Overall, these cases demonstrate the potential of rebound migration to reshape system architectures. Further studies are warranted to explore in depth how the period ratios of systems depend on planet masses and their ordering.

5 Discussion

Our simulations assume a locally isothermal equation of state, meaning the disk temperature varies only with radial distance and remains constant over time. Under this approximation, the corotation torque exerted on the planets is driven entirely by the vortensity gradient. While this approximation is widely used in planet–disk interaction studies, it neglects the thermal evolution of the gas and the entropy-related effects that can contribute to the corotation torque (Paardekooper et al. 2010, 2011). In reality, protoplanetary disks possess finite cooling timescales. The inclusion of a full energy equation would introduce entropy-related corotation torques, which could modify the net torque balance and alter the planets’ migration rates.

The current study relies solely on an X-ray photoevaporation prescription, adopting the mass-loss rate profile from Owen et al. (2012). However, alternative dispersal mechanisms such as far-ultraviolet (FUV) photoevaporation (Gorti et al. 2015) or magnetically driven disk winds (Bai et al. 2016; Suzuki et al. 2016) could produce fundamentally different cavity structures, edge gradients, and overall dispersal timescales (Pascucci et al. 2023). Because the efficacy of rebound migration is highly sensitive to the speed at which the cavity expands, a systematic study of the rebound efficiency under various dispersal models is an important direction for future works.

Furthermore, throughout our 2D simulations, planets are introduced at their full masses without gas accretion. As recently demonstrated by Ida et al. (2026), gas accretion onto a planet can significantly alter the surrounding gas flow and its resulting migration torque. For our massive planet cases in particular, incorporating realistic gas accretion rates could further deplete the local gas surface density, potentially shifting the delicate balance between Lindblad and corotation torques during the critical rebound phase. Although simulations that self-consistently include gas accretion onto planets are beyond the scope of this paper, such models are needed to fully capture the migration physics of massive planets.

Finally, for computational feasibility, we restricted our simulations to a fixed turbulent viscosity parameter of α=10⁻³ and disk aspect ratio of h₀=0.05. The rebound mechanism, as demonstrated in Paper I, is known to depend on these parameters. Moreover, we have only considered two initial period ratios (3:2 and 2:1 resonances) and a limited set of planet mass combinations. A more comprehensive parameter study is beyond the computational scope of the present work, but we caution that our results may not be fully representative of the entire diversity of protoplanetary disks and planetary systems. In particular, the sensitivity of rebound migration to the disk dispersal timescale (illustrated by the comparison between the fiducial and low-mass disk models in Section 4) suggests that a careful exploration of the cavity expansion rate, which depends on both the photoevaporation strength and the initial disk mass, is essential for robust statistical predictions.

Fig. 8 Evolution of disk surface density, planet semimajor axes, period ratios, and resonant angles for a three-planet system consisting of 10 M⊕, 20 M⊕, and 0.3 MJ planets. Top and bottom rows show the configurations of the fiducial disk model and less massive disk model. In the second column, the semimajor axes of the planets from the inner to outer are given by the red, blue and purple curves, respectively, with shaded regions representing their perihelion and aphelion distances. The black dashed line marks the inner disk edge, defined as the radius where the surface density falls below 0.1 g/cm2. In the third column, the blue and red lines give the outer-to-inner planet period ratio for the inner and outer pairs, respectively. The fourth column plots the resonant angles of the inner (red) and outer (blue) planet pairs.

Fig. 9 Period ratios of all simulated systems (left) and a selection of observed extra-solar systems (right). Planets are ordered from inner to outer, color-coded as orange, blue, and purple, respectively, with symbol size representing planet mass. The solid horizontal line marks systems that remain in resonance, while the dashed line indicates period ratios that deviate from resonance. The vertical gray shaded regions denote the 3:2 and 2:1 MMRs. The observed systems are around sun-like stars, with the innermost planet having an orbital period longer than 100 days.

6 Conclusions

In this work, we extended our 2D hydrodynamical investigation of disk-planet interactions to explore how rebound migration operates in multiplanet systems near the expanding cavity of a photoevaporating disk. Using the Dusty FARGO-ADSG code with a source term in the continuity equation that mimics the mass loss rate obtained from Owen et al. (2011, 2012), we simulated systems of two or three planets spanning a wide mass range, from super-Earths to Jupiters, under varying disk conditions. Our simulations confirm that rebound migration, driven by the strong positive corotation torque at the cavity edge, remains effective in multiplanet environments, where it can significantly reshape system architectures. The outcome depends critically on planet masses, mass ordering, and the disk parameters. Our key conclusions are as follows:

• In low-mass pairs (e.g., 5–20 M_⊕), divergent migration driven by the expanding cavity can pull planets out of resonance, leading to period ratios significantly larger than the initial resonant values. Rebound is most pronounced for the outer planet when it is more massive (Figures 1 and 2);

• In mixed-mass pairs (e.g., 20 M_⊕ and 0.3 M_J), the rapid migration of a massive outer planet can bypass wider resonances and establish a stable 3 : 2 MMR that is capable of surviving the disk clearing. However, in cases where the inner planet is more massive, divergent migration caused by cavity expansion can cause the system to become displaced from resonance (Figures 3 and 4);

• In massive pairs (e.g., 0.3 M_J and 1 M_J), strong disk perturbations and eccentricity excitation can modify torque profiles. However, outward migration is still possible in cases where the two planets open a common gap and the inner planet is sufficiently massive to maintain a net positive torque (Figures 5, 6 and 7).

In three-planet systems, we find that rebound migration can alter the system’s architecture during disk dispersal. In our fiducial disk, due to the rapid inward migration of the outer two massive planets, the inner planet undergoes orbital crossing and later gets scattered outward. The outermost planet later rebounds as the cavity sweeps past it, yielding a final architecture of a close-in resonant pair and a distant, lower mass planet. In contrast, in a lower mass disk where photoevaporation drives faster cavity expansion, rebound is suppressed. Here, convergent migration results in the sequential capture and maintenance of a resonant chain without instability, highlighting the sensitivity of system outcomes to disk dispersal timescales (Figure 8).

This study provides the first hydrodynamical evidence that rebound migration can exert a critical influence on the orbital evolution and final architectures of multiplanet systems during the late stages of disk dispersal. By breaking resonances and widening period ratios through divergent migration, rebound migration offers a robust mechanism to transition systems from tightly packed, primordial resonant chains into the nonresonant configurations that characterize much of the observed exoplanet population (Figures 9). Therefore, the late-stage disk dispersal phase is a critical epoch for defining the final observational signatures of multiplanet systems. While a detailed statistical comparison with observational data is beyond the scope of this study, future works should further explore the parameter space of planet masses, disk properties, and photoevaporation models, along with the long-term stability of systems shaped by rebound migration.

Data availability

Movie associated with Fig. 2 is available at https://www.aanda.org.

Acknowledgements

We thank the referee for careful reading and helpful suggestions that improved the manuscript. This work is supported by the National Key R&D Program of China (2024YFA1611803). BL acknowledges support from the National Natural Science Foundation of China (Nos. 12222303 and 12173035) and the start-up grant of the Bairen program from Zhejiang University. The simulations and analysis presented in this article were carried out on the SilkRiver Supercomputer of Zhejiang University. Y.-P.L. is supported in part by the Natural Science Foundation of China (Grants Nos. 12373070 and 12192223), and the Natural Science Foundation of Shanghai (Grant No. 23ZR1473700).

References

All Tables

All Figures

Fig. 1 Evolution of disk surface density, planet semimajor axes, period ratios, and resonant angles for a two-planet system consisting of 5 M⊕ and 20 M⊕ planets. The top (bottom) row shows the configuration in which the inner planet is less (more) massive. In the second column, the semimajor axes of the inner and outer planets are given by the red and blue curves, respectively, with shaded regions representing their perihelion and aphelion distances. The black dashed line marks the inner edge of the expanding cavity, defined as the radius where the surface density falls below 0.1 g/cm2. The fourth column plots the resonant angles 2λ1 − 3λ2 + ϖ2 (blue) and ϖ1 − ϖ2 (red), where the subscripts 1 and 2 denote the inner and outer planets, respectively.

In the text

Fig. 2 Snapshots from the data in Figure 1 of surface density log(Σ) (left), and radial distribution of the individual planet’s torque density T(r) and their total torque (right). Top two rows show results for fid_M02_R1.5 at t=1.5 × 104 P0 and 4 × 104 P0. Bottom row gives the same for fid_M20_R1.5 at t=1.5 × 104 P0. In the right column, the red and blue solid lines are the torque densities of inner and outer planets, whereas the black dashed line represents the total net torque summed over the two planets. The torque profiles are averaged over 20 snapshots within 1 P0. In all panels, red and blue circles indicate the positions of the inner and outer planets, respectively. The video is available online and on https://github.com/bbliu-astro/movies/blob/main/hydrodynamic_rebound/fiducial32_p0p2.gif.

In the text

	Fig. 3 Same as Figure 1, but for a two-planet system consisting of 20 M⊕ and 0.3 MJ planets.
In the text

	Fig. 4 Snapshots from the data in Figure 3 of surface density log(Σ) on the left and radial distribution of the individual planet’s torque density T(r) and their total torque on the right. Top row shows results for fid_M23_R2.0 at t=9 × 103 P0 and the bottom two rows give the same for fid_M32_R2.0 at 9 × 103 P0 and 3.5 × 104 P0, respectively.
In the text

	Fig. 5 Same as Figure 1, but for a two-planet system consisting of 0.3 MJ and 1 MJ planets.
In the text

	Fig. 6 Snapshots from the data in Figure 5 of surface density log(Σ) on the left and radial distribution of the individual planet’s torque density T(r) and their total torque on the right. Top row shows results for fid_M34_R1.5 at t=7 × 103 P0 and the bottom two rows give those for fid_M43_R1.5 at 7 × 103 P0 and 4 × 104 P0, respectively.
In the text

Fig. 7 Azimuthally averaged radial distribution of the disk eccentricity. Filled circles show the planet’s location in each simulation. Dashed part of the curves indicate the gas-depleted cavity region, where the calculation of the gas eccentricity is strongly affected by the reset in the gas density. The Saturn-Jupiter planet pair excites a large gas eccentricity, in contrast to the other three cases.

In the text

Fig. 8 Evolution of disk surface density, planet semimajor axes, period ratios, and resonant angles for a three-planet system consisting of 10 M⊕, 20 M⊕, and 0.3 MJ planets. Top and bottom rows show the configurations of the fiducial disk model and less massive disk model. In the second column, the semimajor axes of the planets from the inner to outer are given by the red, blue and purple curves, respectively, with shaded regions representing their perihelion and aphelion distances. The black dashed line marks the inner disk edge, defined as the radius where the surface density falls below 0.1 g/cm2. In the third column, the blue and red lines give the outer-to-inner planet period ratio for the inner and outer pairs, respectively. The fourth column plots the resonant angles of the inner (red) and outer (blue) planet pairs.

In the text

Fig. 9 Period ratios of all simulated systems (left) and a selection of observed extra-solar systems (right). Planets are ordered from inner to outer, color-coded as orange, blue, and purple, respectively, with symbol size representing planet mass. The solid horizontal line marks systems that remain in resonance, while the dashed line indicates period ratios that deviate from resonance. The vertical gray shaded regions denote the 3:2 and 2:1 MMRs. The observed systems are around sun-like stars, with the innermost planet having an orbital period longer than 100 days.

In the text