© The Authors 2024

1 Introduction

Two key questions raised by the rich body of exoplanet data provided by the Kepler mission were why most exoplanets in multi-planet systems are not in mean motion resonance (MMR) and why there are pileups just wide of exact commensurabilities (Fabrycky et al. 2014). Identifying physical processes that could explain the rarity of exoplanets at exact MMRs when there is a sudden increase in the planets’ number at slightly larger period ratios sparked numerous investigations (e.g. Baruteau & Papaloizou 2013; Petrovich et al. 2013; Delisle & Laskar 2014; Goldreich & Schlichting 2014; Wang & Ji 2017; Charalambous et al. 2022). One pioneering study was carried out by Goldreich & Schlichting (2014), who suggested the majority of observed exoplanets cannot permanently stay in an MMR because their masses are in a regime such that even if they are captured in an MMR, their orbits become overstable in a time comparable to their eccentricity damping timescale (τ_e) and the resonance is eventually broken¹. They conclude that since τ_e is typically much shorter than the migration timescale, τ_a (e.g. Goldreich & Tremaine 1980), the planets spend more time between the resonances rather than in them. This finding triggered many follow-up studies.

Goldreich & Schlichting (2014) analytically investigated the stability of a two-planet system undergoing convergent migration in the presence of planet eccentricity damping by the disc. Assuming the outer planet is more massive than the inner one (a restricted three-body problem) and that their eccentricities are very small, they found three possible fates for the system. Depending on the mass of the outer planet, the system can (a) be captured permanently in resonance with a constant resonance angle and period ratio (stability), (b) stay in resonance while the resonance angle and the period ratio oscillate with a constant amplitude (limit cycle), or (c) be caught in resonance for a time of about τ_e but eventually become overstable. Overstability happens when the amplitudes of the oscillating eccentricities and resonance angles grow until the planetary orbits circulate and the planets come out of the resonance. Goldreich & Schlichting (2014) analytically calculated the criterion for overstable librations and show that it depends on the eccentricity damping timescale, migration timescale, and planet-to-star mass ratio when the restricted three-body problem applies. Delisle et al. (2015) and Deck & Batygin (2015) improved upon the study of Goldreich & Schlichting (2014) by taking the masses of both planets into account. They show that the ratio of the two planets’ masses is a key parameter when determining the outcome of a system with convergent migration and the time spent in resonance when overstable librations lead to escape. Deck & Batygin (2015) categorised the outcome of the system into three regimes, similar to Goldreich & Schlichting (2014) but allowed for a finite, and in fact often small, mass ratio between the two planets. They also studied the long-term evolution of the systems and show that for a planetary pair undergoing overstability, the escape time from the resonance can be 5–10 times longer than τ_e. Moreover, an overstable system goes through successive resonance capturing, with several consequences. Firstly, because the subsequent resonances are more compact, if the system is captured in a stable resonance, the final spacing of the system would be very closely packed. However, if the next resonances are overstable, the system cascades until it reaches the chaotic zone where various MMRs overlap, and the planets are scattered. Secondly, because an overstable system stays in each subsequent resonance for a time longer than τ_e, the system would spend most of its evolutionary time in resonances rather than between them. These results hint that using more realistic models probably leads to even longer staying times in resonances for the system.

One step towards more realistic models is to perform N - body simulations so that the gravitational interaction between the objects can be calculated with fewer simplifications. Xiang-Gruess & Papaloizou (2015) performed a numerical survey to study the long-term evolution of co-planar low-mass two-planet systems near the disc inner edge, which guarantees the convergent migration of the planets. They included the effect of disc-planet tidal interaction with two acceleration terms in the equation of motion of the planets. These terms depend on τ_e and τ_a, which were considered constant during the whole evolution time unless the planets reached the disc inner edge. They find that only a few percent of their models ended in second-order resonances, some of which were overstable. Later, Xu et al. (2018) studied first-order MMR capturing using N -body simulations and taking the eccentricity-related non-linear terms into account when determining τ_e and τ_a (Cresswell & Nelson 2008). They conclude that, compared to models with constant timescales, the eccentricity of the captured planets is larger and the resonant systems tend to be more stable. Recently, Nesvorný et al. (2022) conducted a study on the exoplanetary system TOI-216, finding that a large eccentricity and a 2:1 resonance libration amplitude of the inner planet suggest the presence of a limit cycle during the disc phase. These results highlight the need to investigate the overstability of resonant systems using hydrodynamical simulations.

Hydrodynamical modelling of low-mass resonant planetary systems is very computationally expensive because: (a) a very long running time is needed to bring the planets into resonance and to follow the stability of the system afterwards, and (b) such simulations need high spatial resolution to properly resolve the disc torques on the planets. Among the attempts at modelling the low-mass resonant systems using hydrodynamical simulations (e.g. Papaloizou & Szuszkiewicz 2005; Pierens & Nelson 2008; Pierens et al. 2011; Paardekooper et al. 2013; Hands & Alexander 2018; Ataiee & Kley 2021), overstability has only been reported in Hands & Alexander (2018) and Ataiee & Kley (2021). In the former, the orbital evolution of two moderate-mass planets in a 2:1 MMR in a locally isothermal disc was studied for various disc viscosities using the PLUTO code², as well as the FARGO3D code³ for one test model. In all but one of their models, the planets cross the resonance; the exception is a model run with FARGO3D and described in their appendix. In all the other models, while the Goldreich & Schlichting (2014) overstability condition is marginally satisfied, the planets’ eccentricities are excited as the planets cross the 2:1 MMR and are damped afterwards. However, in the model in their appendix, the planets stay in the resonance for more than 2000 yr; the oscillation amplitude of the eccentricities slowly grows until the resonance is broken. This is the typical behaviour of an overstable system (e.g. Goldreich & Schlichting 2014; Deck & Batygin 2015). Recently, Ataiee & Kley (2021) found some overstable systems in their study of moderate-mass resonant multi-planet systems close to a disc inner edge. They performed locally isothermal hydrodynamical simulations using the FARGO code⁴ and modelled the inner disc edge as either smooth or sharp as the surface density drops. Only for a smooth inner edge with a small disc aspect ratio of 0.03 did they find overstable systems that clearly cascade into higher-order packed resonances. These cases suggest that overstability can indeed happen in more realistic models. The notable point about these two hydrodynamical studies is that the planets are not low-mass planets that undergo type-I migration without opening gaps but instead are moderate-mass planets that open partial gaps.

In hydrodynamical simulations, the planets’ eccentricity damping, migration, and disc structure are all interrelated, in contrast to analytical studies and N-body simulations. When the planets open a (partial) gap, their eccentricity damping and migration timescales are modified. Moreover, the planetary gaps can merge as the planets’ orbits approach each other. This can create asymmetries in the torques on the planets (Masset & Snellgrove 2001) and cause eccentricity damping. The question of how the eccentricity of a gap-opening planet is damped is answered in Pichierri et al. (2023). Another important issue that is automatically included in hydrodynamical simulations is the time evolution of τ_a and τ_e. They are two crucial parameters in forming overstable systems and are undoubtedly time-dependent, especially if the disc viscosity is low. As a planet migrates and/or opens a (partial) gap, the torque on the planet that determines these two timescales changes (Paardekooper 2014; McNally et al. 2019). Whether the criteria suggested by the analytical calculations are valid and how gap opening affects overstability are questions that need further investigation.

In this study, we aim to inspect the overstability at 2:1 MMR for a system with two moderate-mass planets. Because of the complexity of the models in Ataiee & Kley (2021) due to the inner-edge surface density evolution, we started with the model in the appendix of Hands & Alexander (2018), determined the proper initial condition, and performed a parameter study on the disc parameters. We find various outcomes of overstability, limit cycles, and even divergent migration for different disc parameters when the planetary masses remain untouched.

The manuscript is ordered as follows. In Sect. 2, we introduce our nomenclature to avoid any ambiguity. Then we detail our method, initial conditions, and setup. The results are described in Sect. 3. Thereafter, we discuss our findings and compare them with previous work in Sect. 4. Finally, the summary and conclusions come in Sect. 5.

2 Models

In this section, we define some frequently used key terms in the manuscript. Then, we determine the appropriate physical model and initial conditions for our fiducial model and our parameter study. Finally, we introduce the hydrodynamic codes utilised for running the simulations, along with details about the setup.

2.1 Nomenclature

We define several commonly used terms in the manuscript as follows:

• Commensurability. A commensurability is where the period (or mean motion) ratio of two planets equals the ratio of two integers.

• Mean motion resonance (MMR). Two planets are in MMR when, in addition to being in a commensurability, at least one resonance angle stays limited to a span smaller than 2π.

• Stability. We call a system stable in an MMR if the planets’ period ratio, eccentricities, and resonance angle(s) remain bound around an equilibrium value until the end of the simulation. It means that even if these quantities oscillate but their oscillation amplitudes remain constant over time, we consider the system stable.

• Limit cycle. A stable system is in limit cycle if the planets’ period ratio, eccentricities, and resonance angle(s) oscillate around an equilibrium value with constant amplitude.

• Overstability. A system becomes overstable if during the captured time in MMR, the planets’ period ratio, eccentricities, and resonance angle(s) oscillate around an equilibrium value with growing amplitudes until the planets’ orbits circulate and the system comes out of the MMR. Such a system can be subsequently captured in either the next stable or unstable MMRs.

• Resonance passage. When the planets pass an MMR location without being captured, the planets’ eccentricities increase suddenly and the pattern of resonance angle changes slightly but the system does not remain in this configuration and crosses it quickly. We call this state ‘resonance passage’.

2.2 Physical model and initial conditions

Power laws are employed to describe the radial dependence of surface density (Σ), temperature (T), and scale height (H) of our gaseous disc. These profiles are as follows:

$$\begin{array}{l}\text{Σ}(r)={\text{Σ}}_0{\left(\frac{r}{r_0}\right)}^{-{\alpha }_{\text{Σ}}},\hfill \\ T(r)=T_0{\left(\frac{r}{r_0}\right)}^{-\beta },\hfill \\ H(r)=rh(r)=r{h}_0{\left(\frac{r}{r_0}\right)}^f\text{,}\hfill \end{array}$$(1)

where r₀ represents the unit of length chosen to be 1 au in this study. The variables Σ₀, T₀, and h₀ denote the surface density, temperature, and aspect ratio at r₀, respectively. The time of one orbit at r₀ was defined as 1 yr. We assumed that the protoplanetary disc is locally isothermal. Since sound speed c_s is proportional to the square of temperature and aspect ratio is defined as h(r) = c_s/v_K with v_K being the Keplerian velocity, the relation between disc flaring index and temperature slope becomes β = −2f + 1. In this study, the prescription for viscosity in the disc is the α-viscosity model, which is expressed as v = α_vc_sH (Shakura & Sunyaev 1973). To ensure the viscous equilibrium of the system, we selected the initial condition such that the mass accretion rate through the disc remains constant, denoted as Ṁ = 3πvΣ. At a unit length of r₀ = 1 au, its value is approximately Ṁ(r = r₀) ≈ 1.7 × 10⁻¹⁰ M_★ yr⁻¹. This condition requires that α_Σ = 2f + 1/2 in the context of an a-viscosity model. In all models in this study, we adopted an α_Σ value of 1.0, which dictates a flaring index of f = 0.25 for a disc in viscous equilibrium.

In the fiducial model, we set Σ₀ to $1.13\times {10}^{-4}{M}_{\star }/r_0^2$, which corresponds to 1000 g cm⁻² for a host star with a solar mass, M_★ = 1 M_⊙. In addition, we set h₀ to 0.05 and used 10⁻⁵ for α_v.

In our parameter study, we investigated the occurrence of overstability by varying the parameters of the disc. Specifically, we made changes to Σ₀, α_V, and h₀. We adjusted the value of Σ₀ to [1/4,1/3,1/2,2,3] times the value in the fiducial model. We also examined the values of α_v ∈ [10⁻⁶,10⁻⁴,10⁻³]. For h₀, we considered two smaller and one larger values as 0.03, 0.04, and 0.06.

Inspired by the study conducted by Hands & Alexander (2018), our models consist of two super-Earth mass planets initially in circular orbits with no inclination. Accretion onto the planets is not considered. The mass ratio of the inner planet to the star is q₁ = 1.5 × 10⁻⁵, and the mass ratio of the outer planet to the star is q₂ = 3 × 10⁻⁵. This corresponds to 5 and 10 times the Earth’s mass (M_⊕), respectively, assuming the host star has a solar mass. These planets are positioned near 2:1 resonance at orbital radii of a₁ = 1.23 r₀ and a₂ = 2.0 r₀ to save the computational cost. The planets immediately begin migrating, and depending on the specific model, they may eventually reach the 2:1 resonance.

2.3 Codes and setup

In our numerical study, we utilised two different codes to conduct our research. The main code we employed was FARGO3D, developed by Benítez-Llambay & Masset (2016). This code allows us to solve the hydrodynamics equations with a vertically integrated approach. To directly compare our results with a previous study by Hands & Alexander (2018), we also performed simulations of the fiducial model using the PLUTO code (Mignone et al. 2007).

We adopted a 2D cylindrical coordinate system, specifically constructing the (r, ϕ) coordinates to configure the grid for the disc. The azimuthal direction (ϕ) was bounded by [0,2π] to ensure complete coverage of the entire azimuthal domain. For the radial direction (r), the computational domain is limited from [r_in, r_out] = [0.2,7.0] r₀. To determine the optimal value for r_out and ensure the boundary does not impact results while saving computation time, we explored various outer boundaries ranging from 3.1 r₀ (Hands & Alexander 2018, hereafter HA18) to 10 r₀. In Fig. 1, we compare the azimuthally averaged surface density perturbation, ΔΣ_ave/Σ₀, for two identical models with different r_out. Furthermore, both models contain almost the same number of cells in the horseshoe region. The behaviour of models remains identical until approximately 2.3 r₀. However, beyond this point, slight differences emerge due to the proximity of the boundary region to the initial position of the outer planet. To protect against the influence of boundary and damping regions on the torque acting on the outer planet and the overall system evolution, we set r_out to 7 r₀.

To ensure a proper resolution across the disc, we used a computational grid with evenly spaced azimuthal and logarithmically spaced radial directions. We created a computational domain containing N_ϕ = 1024 cells in the azimuthal direction and N_r = 700 cells in the radial direction. The resulting cells are approximately square with this grid configuration. We selected the resolution of our simulation to achieve a scale of 10 cells per scale height at unit of length, resulting in approximately 4 cells within half of the horseshoe region. According to the study of Paardekooper et al. (2010), this number of cells is considered sufficient for resolving the horseshoe region. The chosen resolution constructs a balance between accurately capturing the dynamics of the horseshoe region and computational efficiency. However, we conducted a resolution study in Sect. 3.3.1.

To handle wave reflection in the radial direction near the edges of the domain, we implemented the wave-killing prescription proposed by de Val-Borro et al. (2006). This technique effectively suppresses wave reflections. The damping zones were placed between [0.2–0.22]r₀ and [6.38–7]r₀, providing a buffer region where the waves are attenuated without interfering with the main computational domain. However, in the study by HA18, they damped the surface density to zero within wave-killing zones in their PLUTO simulations.

In order to prevent numerical singularities in the evolution of potential gravity near the planets, we chose a smoothing length of ϵ = 0.6 H_p, as recommended in the study by Müller et al. (2012). This smoothing length ensures smooth and stable computations in the vicinity of the planets. To compare, it is noteworthy that HA18 employed a different value, ϵ = 0.4 H_p, in their study.

In our models, we considered the indirect-term resulting from the planets’ gravity. However, we did not incorporate the gas-indirect term (GIT)⁵ into our analysis. Nevertheless, in Fig. 2, we illustrate the impact of GIT on the occurrence of over-stability. Our findings indicate that the GIT does not significantly affect the occurrence of overstability.

In all models, we implemented the torque correction method presented by Baruteau & Masset (2008, hereafter BM08) to manage the elimination of disc self-gravity. This correction involves subtracting the azimuthally averaged surface density of the disc before calculating the torque exerted on each planet. Additionally, we investigated whether overstability occurs when the BM08 correction is disabled. Figure 2 illustrates that overstability indeed occurs in this case. However, we observe that the planets in this setup reach resonance earlier compared to our fiducial model and they spend less time in resonance before breaking it. We summarise the setup in Table 1.

Fig. 1 Azimuthally averaged surface density perturbation at t = 104 yr for the fiducial model with rout = 7r0 (solid blue line) and 3.1r0, as in HA18 (dotted magenta). The vertical lines show the damping radius. The damping zones were placed between [0.2–0.22]r0 and [2.82–3.1]r0 for a model with rout = 3.1r0.

Fig. 2 Impact of the GIT and BM08 correction on overstability occurrence. Evolution of the semi-major axis of planets (top), the period ratio of the outer planet to the inner one (middle), and the eccentricity of the inner planet (larger eccentricities) and the outer one (shorter eccentricities; bottom). In the top panel, the solid lines represent the inner planet and the dashed lines the outer planet.

3 Results

This section presents the results of our fiducial model, followed by an examination of the impact of viscosity, disc mass, and aspect ratio on the fiducial model. We observed (partial) gaps in some of our models, which had been reported in the recent study by Ataiee & Kley (2021).

Fig. 3 Time evolution of the planets’ orbital properties. From top to bottom: the planets’ semi-major axes, the orbital period ratio of the outer planet to the inner one, and the eccentricity of the inner and outer planets. In the top panel, the solid lines represent the inner planet and the dashed lines the outer planet. The two planets are initially located close to a 2:1 resonance. The horizontal dashed line, r = 0.22, corresponds to the inner damping radius. The vertical dark and light grey lines represent the time when the planets are in resonance. In our fiducial model, the dark grey line corresponds to t = 4840 yr and the light grey line to t = 11 660 yr.

3.1 Fiducial model

We began by examining the fiducial model, which is characterised by an initial surface density of Σ₀ = 1.13 × 10⁻⁴ M_★/r₀², a viscous parameter of α_v = 10⁻⁵, and an aspect ratio of h₀ = 0.05. In this model, two planets are initially placed near 2:1 commensurability, with one planet located at 1.23 and the other at 2.0. After about 5000 yr of convergent migration, the planets enter into a 2:1 resonance that becomes overstable within approximately 6000 yr. Then, after this period, the resonance breaks. Figure 3 provides an overview of the orbital evolution of the two planets in the fiducial model. In the top panel, we plot the semi-major axes of both planets. The planets migrate inwards, and as the outer planet is more massive, it migrates faster compared to the inner one. In the second panel, we illustrate the period ratio of the outer planet to the inner one. After approximately 5000 yr, the planets enter a 2:1 resonance configuration. They remain in this resonant phase for about 6000 yr before eventually escaping the resonance. In the last panel, we display the eccentricity evolution of the planets. Initially, the planets have nearly circular orbits, and as they approach the resonant phase, their eccentricities are excited. Then they start to oscillate around an equilibrium value during the resonant period. The amplitude of these oscillations grows over time until the eccentricity is eventually damped back to near zero as they exit the resonance phase. This behaviour of eccentricity is one of the indicators of overstability in the system.

In Fig. 4, the phase space of the eccentricity and resonance angle of the inner planet is presented. Before the planet’s eccentricity is excited, denoted by the blue spot, the eccentricity was initially almost zero, resulting in esin(ϕ) and ecos(ϕ) values of zero. Once the planet enters the resonance phase, both the eccentricity and resonance angle start oscillating around their respective equilibrium values, e_eq and ϕ_eq. These oscillations are visually represented by the repeating patterns observed in the plot. Over time, the amplitude of these oscillations increases, leading to the size expansion of the patterns in the phase space plot. Upon the planet’s exit from the resonance phase, the eccentricity is damped. Simultaneously, the resonance angle shifts between zero and 2π, causing the appearance of circular patterns centred around zero (the reddish part of the curve). The characteristics of this phase space exhibit similarities to the phase space established by Goldreich & Schlichting (2014), both qualitatively and in terms of the patterns observed. This agreement suggests that the system can be classified as overstable based on the definition provided in Sect. 2.1. We leave the investigation of the equilibrium resonance angle to Sect. 4.

In Fig. 5, we illustrate the azimuthally averaged surface density perturbation caused by planets at various time intervals for our fiducial model. Both planets carve partial gaps (~30%) in the disc, thereby disrupting its initially smooth structure. This is consistent with findings of Crida et al. (2006) and Duffell & MacFadyen (2013), who suggest a higher possibility of gap opening in low-viscosity discs. A detailed analysis and discussion regarding gap opening is presented in Sect. 4. Additionally, we performed a test for the impact of a disc with opened gaps on overstability occurrence as established in Appendix A.

Fig. 4 Phase space of the inner planet for the fiducial model. The colour indicates the time.

Fig. 5 Azimuthally averaged surface density perturbations before, during, and after resonance capture. The dotted dark grey line shows the average surface density perturbation at t = 3000 yr, the solid blue at t = 7500 yr, and the dashed light grey at t = 12 500 yr.

Fig. 6 Influence of viscosity on the evolution of the planetary system for αv = 10−6,10−5,10−4, and 10−3, where 10−5 represents our fiducial model. The figure displays the semi-major axes, the period ratio of planets, and the eccentricity of inner and outer planets from top to bottom. The first row is divided into right and left panels, illustrating the semi-major axes of the inner and outer planets, respectively.

3.2 Parameter study

As part of this study, we investigated the impact of disc parameters, including initial surface density, viscosity, and aspect ratio, on the overstability occurrence of the planetary system. By varying these parameters, we aimed to gain an understanding of how they influence the occurrence and characteristics of the overstability phenomenon in the planetary system.

3.2.1 The effect of viscosity

We explored the impact of viscosity on the overstability of planetary systems by conducting simulations with various values of the viscosity parameter, α_v, equal to [10⁻⁶,10⁻⁴,10⁻³].

Figure 6 presents the results of the models with different viscosity. The two panels in the first row illustrate the evolution of the semi-major axis for the inner and outer planets, respectively, from left to right. The results indicated that increasing the disc viscosity decreases the migration speed of both planets in the disc. In the lower-viscosity models, the evolution of the semi-major axes of both planets overlaps. However, for the model with the highest viscosity in this study (α_v = 10⁻³), the planets migrated more slowly than the others. This happens because the positive co-rotation torque is larger for this value of viscosity. We discuss this point further in Sect. 4.2.

In the second panel of Fig. 6, which shows the period ratio of the two planets, we observed that in all models, planets entered the 2:1 resonance. In the low-viscosity models, we observed overstable behaviour. However, in the highest-viscosity model in this study (α_v = 10⁻³), the system remained in resonance. For α_v < 10⁻³, higher viscosities led to earlier captures into the 2:1 resonance and a shorter duration of the resonance phase. In the α_v = 10⁻³ model, where planets only captured into resonance without showing signs of overstability, resonance capture occurred very late, and the planets remained in the resonance phase throughout the entire simulation period.

In the last two panels of Fig. 6, we demonstrate the effect of viscosity on the eccentricity evolution of the planets. When the viscosity parameter, α_v, is less than 10⁻³, both planets display typical behaviour of overstable planetary systems. This behaviour is characterised by eccentricity enhancement followed by oscillations around an equilibrium value during the resonance phase. Before the resonance breaking and eccentricity damping, the amplitude of these oscillations grew significantly. The average eccentricity value of the inner planet was approximately ten times higher than that of the outer planet. In the model with α = 10⁻³, the eccentricity remained at an equilibrium value, approximately half of the value observed in other models, but it did not dampen. To investigate this behaviour, we extended the simulation time of the α = 10⁻³ model up to three times the typical simulation duration used for our fiducial model.

Figure 7 illustrates the time-extended simulation of the α_v = 10⁻³ model. The eccentricity exhibited oscillations around an equilibrium value. The amplitude of oscillations initially increased and then decreased. Eventually, the eccentricity stabilised at a constant value without experiencing damping. This phenomenon is known as a limit cycle, where the system’s behaviour formed a repetitive pattern over time.

To further confirm the occurrence of a limit cycle in the α_v = 10⁻³ model, we present its phase space in Fig. 8. This plot illustrates the excitation of eccentricity from zero, followed by oscillations in both the resonance angle and eccentricity around an equilibrium value. Although the amplitude of these oscillations initially increased and then decreased, the eccentricity and resonance angle consistently remained at the equilibrium value, never damping to zero. The phase space plot clearly confirms the presence of the limit cycle.

Fig. 7 Period ratio and eccentricity evolution of planets for the model with α = 10−3. To track the evolution of the α = 10−3 model, we extended the simulation time up to 45 000 yr. The presence of a limit cycle is evident in the eccentricity evolution, denoted by persisting eccentricity oscillations without damping to zero.

3.2.2 The effect of disc surface density

In this section we explore various disc masses, specifically [$\frac{1}{4}$, $\frac{1}{3}$, $\frac{1}{2}$, 2, 3] times Σ₀. We aimed to understand the influence of disc mass variations on the orbital properties of the planetary system.

In Fig. 9, we compare the fiducial model with models with initial surface densities of $\left[\frac{1}{4},\frac{1}{3},\frac{1}{2}\right]{\text{\hspace{0.17em}Σ}}_0$. The first panel, displaying the semi-major axes of the inner and outer planets. In the second panel, we show the period ratio of the outer planet to the inner one. For $\left[\frac{1}{2},\frac{1}{3},\frac{1}{4}\right]{\text{\hspace{0.17em}Σ}}_0$ models, the planet enters 2:1 resonance at approximately t ≈ [8800, 13 200, 17 800] yr. As a result, planets embedded in discs with higher surface density are captured into the 2:1 resonance more rapidly. They stay in the resonance phase for a shorter duration, for $\left[\frac{1}{2},\frac{1}{3},\frac{1}{4}\right]{\text{\hspace{0.17em}Σ}}_0$ models, the residence time in the resonance phase is about ∆t ≈ [16 200, 69 200, 81200] yr. In the $\frac{1}{2}{\text{Σ}}_0$ model, the migration was converging similarly to the fiducial model, while for the $\left[\frac{1}{3},\frac{1}{4}\right]{\text{\hspace{0.17em}Σ}}_0$ models, the migration diverge after t ≈ [82 400, 99 000] yr. According to Kanagawa & Szuszkiewicz (2020), this divergence can be attributed to the profound gaps created by the planets, estimated to be approximately 70 to 80% of the disc’s initial surface density for the outer planet in each model.

In the last two panels, we illustrate the eccentricity evolution of the inner and outer planets. These panels demonstrate the occurrence of overstability in planetary systems with surface density lower than Σ₀. As shown in the pioneering studies (e.g. Tanaka et al. 2002), torque is directly proportional to surface density. Torque is also inversely proportional to the migration timescale, ${\tau }_{\text{a}}^{-1}$. Therefore, a reduction in surface density results in a smaller torque and an increased migration timescale. Our results show that the duration of resonance is longer in lower-surface-density models.

The time at which planets exit overstability differs from the time they exit 2:1 resonance. For $\left[\frac{1}{3},\frac{1}{4}\right]{\text{\hspace{0.17em}Σ}}_0$ models, the duration that planets remain in overstability is about ∆t ≈ [28 800, 49 200]. After planets exit overstability, the resonance angle is not restricted to a limited value. It changes between [0, 2π], as illustrated in Figs. B.1 and B.2. Thus, the $\left[\frac{1}{3},\frac{1}{4}\right]{\text{Σ}}_0$ models exhibit overstability until approximately t ≈ [42 000, 67 000] yr. After this period, the resonance angle circulates and no longer oscillates around its equilibrium value. Consequently, the planets are no longer overstable, they simply remain in 2:1 commensurability. After a considerable time, they eventually exit this commensurability.

In Fig. 10, we present the results of increasing the initial surface density of the disc and its impact on the occurrence of overstability in planetary systems. The first panel illustrates the planets’ inward migration, demonstrating the expected linear relationship between surface density and migration rate. As evidenced in the second panel, all models reach the 2:1 resonance. However, contrary to other models, the 3Σ₀ model diverges upon exiting from resonance. This divergence may be attributed to the wake-planet interaction, which is effective when planets open partial gaps (Baruteau & Papaloizou 2013).

The last two panels depict the eccentricity evolution of the inner and outer planets’ orbits. Overstability was evident in the 2Σ₀ model, whereas in the 3Σ₀ model, oscillations did not grow around a specific equilibrium value, indicating the absence of overstability. According to the eccentricity evolution in these panels, it is evident that increasing the surface density leads to a decrease in the amplitude of eccentricity oscillations. This indicates that a higher surface density enhances the efficiency of eccentricity damping. In the 3Σ₀ model, the eccentricity damping efficiency lies within a moderate range. It is not efficient enough to cause a permanent resonance and not low enough to allow the amplitude of oscillations to grow, resulting in overstability.

As depicted in Figs. 9 and 10, the variation in the surface density of the protoplanetary disc has a significant impact on the occurrence of overstability, the stability of planetary systems, and their overall structure. In Sect. 4 we explore in detail how variations in surface density influence the architecture of planetary systems.

Fig. 8 Phase space of the inner planet for the α = 10−3 model. The colour indicates the evolutionary time. The presence of a limit cycle is clearly evident in the phase space, denoted by the repeating oscillations of the inner planet’s orbit around a stable equilibrium value.

Fig. 9 Impact of a decreasing surface density of the disc on the orbital properties of planets. From top to bottom: evolution of planets semi-major axes, the period ratio of planets, and the eccentricity of the inner and outer planets. The dashed lines represent the outer planet, while the solid lines represent the inner planet. The dotted line, which corresponds to r = 0.22, represents the inner damping radius.

Fig. 10 Impact of an increasing surface density of the disc on the orbital properties of planets. From top to bottom: evolution of planets semi-major axes, period ratio of planets, and eccentricity of the inner and outer planets. The dashed lines represent the outer planet, while the solid lines represent the inner planet. The dotted line, r = 0.22, represents the inner damping radius.

Fig. 11 Impact of increasing the aspect ratio on the evolution of the planetary system. From top to bottom: evolution of planets semi-major axes, the period ratio of planets, and the eccentricity of the planets. The solid lines represent the inner planet, while the dashed lines represent the outer planet.

3.2.3 The effect of the aspect ratio

To explore the influence of disc thickness and temperature on the occurrence of overstability, we performed simulations with h₀ = 0.03,0.04, and 0.06. Figure 11 illustrates the impact of increasing the aspect ratio to h₀ = 0.06 on the planetary system. In the first panel, we observe that the migration speed slowed with increasing disc thickness, as expected (e.g. Paardekooper et al. 2010). The second panel shows that the planets reached the 2:1 resonance. However, the capture into resonance occurred later than in the fiducial model, and the planets remained in resonance for a longer period. To fully capture the system’s evolution, we extended the simulation runtime. In the last panel, we illustrate the eccentricity evolution of the inner and outer planets. The occurrence of overstability was evident in these plots. Initially, the eccentricity was excited, leading to the growth of oscillations. However, with time, the eccentricity starts to dampen, and eventually, the eccentric orbits turn to circular orbits.

In Fig. 12, we present the influence of decreasing the aspect ratio of the disc on the evolution of the planetary system. The first panel displays the inward migration of both planets. Contrary to our expectations from type-I migration, decreasing the aspect ratio did not accelerate migration. Specifically, for the h₀ = 0.03 model, the migration speed decreased significantly, while for the h₀ = 0.04 model, it remained comparable to our fiducial model. In the second panel, we present the period ratio of the planets. Initially, for the first few years of the simulation, there was an agreement between the period ratio of h₀ = 0.04 and the fiducial model. However, after a few years, the migration diverged, and the planets did not even capture into resonance when the aspect ratio was decreased for both models. Instead, they remained in nearly circular orbits throughout the entire simulation without experiencing any overstability. These results indicated that decreasing the aspect ratio of the disc has a significant impact on the migration behaviour and resonance capture of the planets. Gap openings in these models surely led to these results, as discussed further in Sect. 4.

In Table 2, we provide a simple overview of the outcomes of our parameter study, indicating whether planet captures into resonance, overstability, or limit cycle occurred.

Fig. 12 Influence of a decreasing aspect ratio on the evolution of the planetary system. The first panel shows the evolution of the semi-major axis and the second panel the period ratio of the planets.

3.3 The effect of numerics

In this section, we explore the impact of numerical parameters on the occurrence of overstability. Initially, we conducted a resolution study to assess how spadai resolution affects the occurrence of overstability. Following this, we compared our results with similar simulations carried out with the PLUTO code. While recognising the importance of numerical parameters, we observe the occurrence of overstability across the majority of cases.

Fig. 13 Effect of resolution on overstability occurrence.

3.3.1 Resolution study

To understand how resolution impacts overstability occurrence, we performed simulations at various resolutions while maintaining other parameters consistent with our fiducial model. These resolutions ranged from half to twice that of our fiducial model. The results of these simulations are summarised in Fig. 13 and Table 3.

We observe overstability occurrence for cpx_s = [3.26,4,5,6], where cpx_s stands for cells per horseshoe half-width. For cpx_s = [2,3], the overstability did not occur. However, for cpx_s = [3], the planets reached 2:1 resonance. In cpx_s = 7.8, which is our highest resolution, planets reached the 2:1 resonance, but overstability did not happen during the simulation time, which was about 30kyr. However, from the eccentricity evolution of the planets (and the resonance angle libration, which for the sake of brevity is not shown here), we expect that this model will eventually show overstable behaviour. The difference can be due to how properly the torque on the planet is resolved. It is known that the lowest cpx_s for resolving the horseshoe drag is 4 (Paardekooper et al. 2010), but we are not aware of how increasing the resolution changes the horseshoe drag on the planet. Further investigation is needed to properly study how the resolution change the torque on the planet that is beyond the scope of this paper.

3.3.2 Comparison with PLUTO

To corroborate our results further, we compared our fiducial model against the numerical hydrodynamics code PLUTO v4.4 (Mignone et al. 2007). By default, PLUTO uses a Godunov-type, finite-volume approach with a second-order space- and time-accurate scheme to advance the hydrodynamics equations. To speed up calculations, and to facilitate a fair comparison to FARG03D, we also used the FARGO module implemented in PLUTO by Mignone et al. (2012).

For our fiducial setup, we chose the HLLC Riemann solver (Toro et al. 1994; a method in computational fluid dynamics for efficiently solving the Riemann problem, incorporating contact waves to capture flow discontinuities and complex phenomena accurately), and second-order reconstruction and time-stepping schemes (LINEAR and RK2, respectively). The planets and star were integrated with an RK4 scheme similar to Thun & Kley (2018, see Appendix A therein), and the back-reaction of the disc on the star was included. Since we neglected disc self-gravity, we also accounted for the correction to the planetary torques by Baruteau & Masset (2008).

Our findings are shown in Fig. 14. While our fiducial model shows an eccentricity growth for both planets similar to FARG03D, the planets are unable to break the 2:1 resonance. Increasing the resolution by a factor of 1.4 in both directions (for a resolution of N_r × N_ϕ = 990 × 1448 cells), however, allows the planets to break the resonance, indicating that the overstability is operating here. Alternatively, opting for a third-order, weighted essentially non-oscillatory (WENO) reconstruction (Yamaleev & Carpenter 2009) at the fiducial resolution also reproduces the overstability, while also matching the amplitude of the eccentricity growth for both planets.

While not shown in the figure, we tested various other combinations of numerical options. We find that, at the fiducial resolution and regardless of numerical setup, the WENO reconstruction is necessary for the planets to break the resonance. At a higher resolution, however, the default numerical setup is sufficient. This suggests that the numerical scheme used by PLUTO is not able to capture the overstability at low resolution, but that the overstability is not a numerical artefact. This sensitivity to the numerical method between PLUTO and FARGO has also been observed in the context of buoyancy waves (Ziampras et al. 2023), but showed numerical convergence at sufficient resolution.

We note that the self-gravity correction by Baruteau & Masset (2008) is essential to trapping the two planets in a 2:1 resonance in the first place. Without the correction, the planets pass through the 2:1 resonance with only a brief excitation to their eccentricity similar to Hands & Alexander (2018), but can still be later captured into 3:2 resonance after around ~8000 yr.

Fig. 14 Time evolution of the eccentricity (top) and period ratio (bottom) of the two planets for various configurations of PLUTO, compared to the fiducial FARGO3D setup presented in Table 1. The two planets do not break out of resonance in the fiducial PLUTO setup, but a higher resolution or a third-order reconstruction (WENO) both recover the behaviour found by FARG03D.

4 Discussion

4.1 Equilibrium eccentricity and resonant angle

The resonant angle and eccentricity oscillate around an equilibrium value, as defined in Sect. 2.1. Figure 15 illustrates the evolution of the resonant angle, with the colour scheme indicating changes in the eccentricity of the inner planet. In all of our overstable models, we can estimate the times of entry into and exit from resonance when the eccentricity of the inner planet is approximately e_in = 0.015. Consequently, by calculating the average of resonant angle in the resonance phase, we determined that ϕ_eq = 22.54° represents the equilibrium value of resonant angle during planetary resonance for our fiducial model.

In Fig. 16, we present the equilibrium values of the resonant angle, the inner and outer planet eccentricities, and the eccentricity ratio for all of our models in 2:1 MMR, regardless of whether they exhibit overstable behaviour. Vertical lines separate different model categories. The left segment comprises viscous models, the middle contains the disc surface density models, and the right segment includes the aspect ratio model.

According to earlier literature (e.g. Tanaka et al. 2002), there exists a simple linear relationship between the eccentricity damping timescale and the migration timescale, τ_e = h²(r)τ_a. Additionally, based on the study by Goldreich & Schlichting (2014), the equilibrium resonance angle is given by sin ϕ_eq ∝ 1/τ_e. Thus, a decrease in τ_e leads to a reduction in the equilibrium resonance angle. In viscous models, the equilibrium resonant angle is nearly identical across all three overstable models. As illustrated in Fig. 6, the overstable viscous models exhibit approximately similar migration rates, resulting in almost identical migration timescales. In contrast, the equilibrium resonant angle value for the limit cycle model is half that of the fiducial model. This discrepancy is reasonable because of the low migration rate associated with the limit cycle model. In models investigating surface density, the equilibrium resonant angle decreases as the migration timescale decreases with a reduction in the disc’s surface density. According to the equilibrium value of our non-overstable model, we can conclude that models with ϕ_eq ≈ 50° and higher do not experience overstability. Models examining the aspect ratio reveal a notable decrease in the equilibrium resonant angle and a decreasing migration timescale with an increasing aspect ratio. Therefore, our results align with earlier studies in these cases.

According to the equilibrium eccentricity relation established in Goldreich & Schlichting (2014), e_eq ∝ (τ_e/τ_a)^1/2, for our fiducial model, the predicted value of e_eq is approximately 0.02, which deviates from the simulated value of about 0.03. Similar discrepancies were observed in other models when comparing the e_eq values predicted by Goldreich & Schlichting (2014) to the simulation results. These contrasts highlight the need for modifications to the e_eq formulation for application in hydrodynamic simulations. Additionally, while eccentricity damping rates and migration rates are given values in N-body simulations and analytical calculations, determining these rates in hydrodynamics involves a more complex process. In hydro-dynamic simulations, such as those discussed in this study, eccentricity damping and migration are indeed computed self-consistently through applying the forces from the disc and other planets at every time step.

In Goldreich & Schlichting (2014), a criterion is proposed that if $e_{\text{eq}}\ge q_{\text{in}}^{1/3}$, the resonance is considered temporary. In our study, where $q_{\text{in}}^{1/3}\approx 0.024$, this agrees with our simulations. All our simulations, except the limit cycle and non-overstable models, exhibit equilibrium eccentricity values exceeding the specified threshold. While the limit cycle model seems to be permanently in resonance, the non-overstable model does not exhibit permanent resonance. This suggests that the criterion may be specifically indicative of the presence of overstability in models.

The bottom panel in Fig. 16 illustrates the ratio of the inner planet’s eccentricity to the outer one. With the exceptions of the limit cycle model, non-overstable model, and models with α_v = 10⁻⁴ and Σ(r₀) = 2Σ₀, which closely match the value of e_in/e_out ≈ 7 from Goldreich & Schlichting (2014), other models exhibit variations. This suggests that eccentricity and capturing in an exact MMR depend not solely on the planet mass ratio and position but also on disc properties and the interaction between the disc and the planet.

Fig. 15 Time evolution of resonant angle for the fiducial model. The colour represents the eccentricity of the inner planet. The horizontal line represents the equilibrium value of the resonant angle during the resonance period of planets, which is ϕeq = 22.54°.

Fig. 16 Top to bottom: equilibrium values for the resonant angle, eccentricities of the inner and outer planets, and the ratio of the planet’s eccentricity for the models that captured them in a 2:1 MMR. The viscosity models, initial surface density models, and aspect ratio models are depicted from left to right, separated by vertical lines. The square, diamond, cross, and circle markers represent the fiducial model, the overstable models, the non-overstable model captured in the 2:1 MMR, and the limit cycle model, respectively. Additionally, Σ1 = Σ(r0).

4.2 The role of co-rotation torque

In Fig. 6, we see a considerable trend in the migration rates of models with varying viscosity. We find that the migration rates are nearly identical for models with α_v < 10⁻³. However, for the model with α_v = 10⁻³, the migration rate is remarkably low compared to the others. To understand the origins of this behaviour, we examined the torque exerted on planets. Our analysis was based on equations found in the studies for Lindblad torques, as well as unsaturated and saturated co-rotation torques (e.g. Goldreich & Tremarne 1980; Paardekooper et al. 2010, 2011). The scaled Lindblad torque, Γ_lin/Γ₀ ≈ −2.437, remains unaffected by variations in viscosity, aligning with expectations. In contrast, the co-rotation torque can become saturated with varying viscosity. For α_v = [10⁻⁶, 10⁻⁵, 10⁻⁴, 10⁻³], the scaled saturated co-rotation torque is approximately [0.003, 0.022, 0.141, 0.466]. Consequently, it becomes crucial to investigate the magnitude of the saturated torque. The value of the unsaturated co-rotation torque remains constant across all viscous models, equal to 2.997. The positive value of the saturated co-rotation torque reduces the absolute value of the total torque to [2.434, 2.415, 2.293, 1.971], resulting in a slower migration of planets, especially in the model with α = 10⁻³.

Fig. 17 Time evolution of gaps’ positions and depths. Top: time evolution of the semi-major axis and the position at which surface density has its minimum value. Bottom: time evolution of the surface density in the absence of the planet (disc’s power-law surface density) in blue, and the surface density in the position of planets in dark and light red for inner and outer planets.

4.3 Gap opening

In our study, we observe that planets create partial gaps around their orbits. In this section we focus on examining these gaps in some of our models, including the fiducial model, the one with Σ(r₀) = Σ₀/4, and h₀ = 0.03.

Figure 17 depicts the evolution of gaps in our fiducial model. The top panel illustrates the time evolution of the planet’s semi-major axes and the position of the bottom of the gaps. We determined the gap’s bottom position using an algorithm that identifies a local minimum in the azimuthally averaged ID profile of surface density around the planet’s orbit. The outer gap, situated beyond the orbit of the outer planet, corresponds to the study by Kanagawa et al. (2020). Initially, the first gap formed inside the orbit of the inner planet and gradually moved closer to it. However, the gap’s position shifted outside the inner planet’s orbit around 10 000 yr. The bottom panel illustrates the evolution of the disc surface density in the absence (following the power-law relation) and the presence of planets at their positions. The surface density at the initial point, where no gap exists (Σ_n), increases over time. This increase is attributed to inward migration, leading to movement towards areas of higher surface density. In contrast, in the presence of planets, the surface density at the planet’s position (Σ_pl) decreases due to gap formation.

In Fig. 18, we display the surface density of gaps and its perturbation, expressed in terms of the position of the bottom of the gap. The colour signifies time. The outer planet’s gap is consistently located beyond its orbit. Conversely, the inner planet’s gap initially forms within its orbit but gradually moves outwards over time, passing the planet’s position.

The inner planet creates a deeper gap compared to the outer one. According to earlier studies (e.g. Crida et al. 2006; Duffell & MacFadyen 2013), massive planets, low viscosities, and low aspect ratios form deeper gaps. Despite the lower mass of the inner planet, its location in a region with a lower viscosity and a lower aspect ratio leads to the formation of a deeper gap.

To continue our investigation, we explored the reason behind the placement of the inner planet’s gap inside planet’s orbit, which differs from the finding by Kanagawa et al. (2020). Figure 19 illustrates the parameter $\frac{\text{Σ}}{{\text{Σ}}_0}$, where Σ and Σ₀ denote the values of disc surface density at a specific time point before planetary resonance (t = 3000 yr) and the initial time (t = 0), respectively, for our fiducial model. There are two additional models in this comparative analysis, each with a single planet. One model includes a planet with the mass of the inner planet in our fiducial model, while the other incorporates a planet with the mass of the outer planet in the fiducial model. Both planets in these models are initially positioned in the same locations as those in the fiducial model. In Fig. 19, it is evident from the single-planet simulations that the gap for both planets lies outside their respective orbits, in agreement with the findings of Kanagawa et al. (2020). However, in the fiducial model with two planets in resonance, we noticed a distinct behaviour, concerning the gap of the inner planet. This deviation may arise due to the presence of the outer planet. The outer planet, with its notable mass, behaves independently of the inner planet, resulting in its gap opening resembling that of a single planet. On the other hand, the gap opening for the inner planet is different in the models with single and double planets. The gold line shows the gap profile for a test simulation with identical planet locations and evolution time but null eccentricity for the planets. It shows that the eccentricity of the planets are not responsible for the change in the gap profile.

After analysing the gap structure of the fiducial model, we shifted our focus to the gap configurations of the Σ(r₀) = Σ₀/4 model. Figure 20 presents the surface density within the gap and its perturbation for the Σ(r₀) = Σ₀/4 model. From previous studies (e.g. Goldreich & Tremaine 1980), we know that decreasing the initial surface density of the disc leads to slow migration. Figure 9 highlights that a slower migration speed extends the overstability process. Consequently, planets have more time to deepen the gaps, as seen in Fig. 20. Contrary to our initial expectation that different surface densities would merely scale the overstability, the gaps deepen due to the longer time of evolution and a low-viscosity parameter, as discussed by Duffell & MacFadyen (2013). Similar to the fiducial model, the outer planet’s gap is located beyond the planet, whereas the inner gap transitions from inside the orbit to the outside.

In the following, we investigate the gap structure within the non-resonance h₀ = 0.03 model to understand its impact on migration speed. According to early studies (e.g. Goldreich & Tremaine 1980), the scaled torque is proportional to the inverse square of the aspect ratio, Γ₀ ∝ h⁻². Consequently, a planet in a disc with a lower aspect ratio exhibits a higher migration rate. However, our simulations reveal the opposite trend. In Fig. 21, we present a 2D snapshot of the disc’s surface density at different time points. Even within the initial few thousand years of simulation, planets in this model create notably more profound gaps compared to gaps in the fiducial model, confirming the results from the literature for low aspect ratios (e.g. Crida et al. 2006). These substantial gaps are responsible for the observed slow migration in discs with lower aspect ratios. This is essentially transiting from Type-I to Type-II migration.

Fig. 18 Evolution of the surface density of gaps and its perturbation in terms of positions where surface density has its minimum value. The colour corresponds to time.

Fig. 19 Azimuthal average of disc surface density for multiple models at t = 3000 yr. We include our fiducial model (FM), two other models with single planets (one with a low-mass planet and the other with a large-mass planet), and another model identical to the FM but with zero eccentricity. The orange circles represent planets in the single-planet models, while the red circles depict planets in the FM.

Fig. 20 Evolution of gap surface density and its perturbations in terms of the positions of the gap bottoms for the Σ(r0) = Σ0/4 model. The stars indicate the times when planets enter resonance.

Fig. 21 Snapshots of the disc’s surface density over time at 2000, 20 000, and 30 000 yr for the h0 = 0.03 model. The position of gaps is marked by gold circles, and planets by green circles.

4.4 Comparison with pure type-l migration

Since the migration torques play a key role in the resonance evolution and occurrence of overstability, we assess in this section the proximity of our planets’ migration rate to pure type-I migration, which is a situation that no gap formed in disc. Given that we have observed gap formation in our planetary systems, we sought to understand how the torque, as calculated for systems with a gap (e.g. Kanagawa et al. 2015, 2018; Duffell 2015), influences the migration of our planets. For this purpose, we compared our fiducial model’s orbital evolution with the study by Paardekooper et al. (2010), where no gaps are formed, and with the study by Kanagawa et al. (2018), accounting for gap formation in the disc (see Kanagawa et al. 2018, Eq. (27)). It is worth mentioning that we used the same smoothing length as in our hydrodynamic simulations when evaluating Type-I migration. As the co-rotation torque is nearly saturated in this model, the effect of gap opening would be mostly on the Lindblad torque.

In top panel of Fig. 22, we compare the orbital evolution of both the inner and outer planets obtained from our simulations with the orbital evolution of a planet if the torques on them are calculated according to the prescription in Paardekooper et al. (2010) and Kanagawa et al. (2018). Our simulations agree with the semi-major axes evolution according to pure type-I migration up to a certain point. However, a significant deviation is observed when applying a gap by Kanagawa et al. (2018) formula. This suggests that in torques correction by Kanagawa et al. the depth of the gap for the moderate planetary masses are over-estimated. Concerning the inner planet, its orbital evolution closely follows pure type-I migration until the planet exits resonance. Beyond this point, the differences between the two become apparent. For the outer planet, our simulation diverges from pure type-I migration once the planets reach resonance.

To investigate the difference between our fiducial model and pure type-I migration, we performed two additional simulations with only a single planet in each run, with one simulation for the inner low-mass planet (solid line) and one for the outer high-mass planet (dashed line), in red in top panel of Fig. 22.

Comparing these single-planet simulations shows a closer agreement with our fiducial model, suggesting that gap formation has minor effect on the migration rate of the planet.

To further assess this difference, we analysed the evolution of the quantity (a_p – a_sim)/a_p, where a_sim and a_p represent the semi-major axis of the planet in our simulations and pure type-I migration. The bottom panel of Fig. 22, indicates that the discrepancy between our fiducial model and pure type-I migration is mainly attributed to resonance and gravitational interactions to a minimal extent.

In the following analysis, we aim to assess whether the model with $\text{Σ}\left(r_0\right)=\frac{{\text{Σ}}_0}{4}$, characterised by the formation of a deep gap as depicted in Fig. 20, agrees with the migration rate predicted by Kanagawa et al., for discs exhibiting significant gaps. Figure 23 presents a comparison between the semi-major axis of this model, pure type-I migration, and migration with gap opening. A substantial deviation between the simulation and the Kanagawa et al. theory is observed, particularly in the early stages. This suggests that the long-term evolution of the system led to the formation of these pronounced gaps (e.g. Duffell & MacFadyen 2013), and the deviation from pure type-I migration in later times indicates that type-I migration in this model is no longer pure.

In the subsequent analysis, we compare the simulation results of the h = 0.03 model with semi-analytical predictions. As illustrated in Fig. 24, a remarkable agreement is observed between the semi-major axis of the h = 0.03 model and the prediction presented by Kanagawa et al., for discs featuring significant gaps. This close correspondence strongly supports the presence of a deep gap in the structure of this disc. Despite the planets having a moderate mass, the migration behaviour closely resembles type-II migration. The figure indicates that the semi-major axis, that predicted by Kanagawa et al., exhibits almost no migration.

This is attributed to the Kanagawa et al. expectation that $\frac{{\text{Σ}}_{\min }}{{\text{Σ}}_0}$ equals 0.03 for h₀ = 0.03, signifying a gap with a depth of 97 percent. In such a deep gap, the planets essentially orbit around a fixed semi-major axis.

Fig. 22 Comparison of semi-major axes obtained through simulation with different theoretical models. Top: orbital evolution of our FM, along with two semi-theoretical models from Paardekooper et al. (2010, without a gap opening) and Kanagawa et al. (2018, with a gap opening), alongside two additional simulations using the same disc properties as the FM but focussing on individual planets. Solid lines represent the semi-major axis of the inner planet, while dashed lines represent the outer planet. Bottom: evolution of the quantity (ap – asim)/ap, where ap represents the semi-major axis according to pure type-I migration and asim denotes the semi-major axis of both the inner and outer planets in our FM, as well as the single small and big planets in our single-planet simulations. The period of resonance is indicated between the dark and light grey time intervals.

Fig. 23 Same as the top panel of Fig. 22 but for the Σ(r0) = Σ0/4 model.

Fig. 24 Same as the top panel of Fig. 22 but for the h0 = 0.03 model.

5 Conclusion

We have investigated the overstability in a 2:1 resonant planetary system, considering two planets with moderate masses of 5 and 10 M_⊕, using hydrodynamic simulations. Through an organised parameter study, we found that disc properties, especially mass and viscosity, play a key role in the occurrence of overstability. Low-viscosity discs exhibited overstability, while a high viscosity led to limit cycle phenomena. Variations in the disc mass had a profound impact, with low-mass discs consistently displaying overstability. In these models, planets have enough time to open deeper gaps compared to the higher-surface-density models. For lower-aspect-ratio discs (e.g. h = 0.03), the planets formed very deep gaps and underwent type-II migration. Our study takes a step closer to reality and shows that overstability is not merely a theoretical construct, especially in low-mass and low-viscosity discs: it can aid in forming the packed resonant system. We would like to emphasise that our results are particularly relevant for planets undergoing convergent migration and crossing the 2:1 MMR at approximately one astronomical unit separations from the star.

In a future study, we aim to calculate migration and eccentricity timescales in our hydrodynamic simulations and compare the results with theoretical calculations of previous studies carried out using our hydrodynamical simulations.

This study pushes us closer to understanding the complex interaction between planetary migration and disc properties. Overstability, previously restricted to theoretical domains, is now on the threshold of empirical confirmation.

Acknowledgements

This work was supported by Ferdowsi University of Mashhad under Grant No. 3/58699 (1401/11/30). The authors acknowledge support by the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant no. INST 37/935-1 FUGG. The authors thank Cornelis Dullemond and Gabriele Pichierri for their helpful discussions. In heartfelt remembrance, ZA expresses her deep gratitude to Willy Kley, whose unwavering guidance, support, and encouragement have deeply influenced her life journey. We would like to thank the referee O. Chrenko for his careful comments which helps to improve the manuscript.

Appendix A Impact of a disc with opened gaps on overstability occurrence

To examine the influence of opened gaps on overstability occurrence, we conducted three simulations using different models, including our fiducial model, the α_v = 10⁻⁶ model, and the h₀ = 0.03 model. In these simulations, we initially allowed the system to evolve for 5000 years without migration. Subsequently, after the planets opened gaps, we enabled migration for the rest of the simulation time. The results are depicted in Fig. A.1. It is apparent that for models experiencing overstability simultaneously with gap formation, the overstability reoccurred, albeit with slight variations in resonance braking time. Additionally, in the model where overstability did not occur together with gap opening, we did not observe overstability again.

Fig. A.1 Impact of the disc containing opened gaps on the occurrence of overstability. The solid lines show the models that the evolution of the system is together with gap formation. The dotted lines illustrate the impact of the disc with opened gaps on the system.

Appendix B Resonance angle of low-surface-density discs

Fig. B.1 Resonance angle evolution for the $\text{Σ}\left(r_0\right)=\frac{1}{3}{\text{Σ}}_0$ model. The colour indicates the eccentricity of the inner planet.

Fig. B.2 Resonance angle evolution for the $\text{Σ}\left(r_0\right)=\frac{1}{4}{\text{Σ}}_0$ model. The colour indicates the eccentricity of the inner planet.

References

All Tables

All Figures

	Fig. 1 Azimuthally averaged surface density perturbation at t = 104 yr for the fiducial model with rout = 7r0 (solid blue line) and 3.1r0, as in HA18 (dotted magenta). The vertical lines show the damping radius. The damping zones were placed between [0.2–0.22]r0 and [2.82–3.1]r0 for a model with rout = 3.1r0.
In the text

Fig. 2 Impact of the GIT and BM08 correction on overstability occurrence. Evolution of the semi-major axis of planets (top), the period ratio of the outer planet to the inner one (middle), and the eccentricity of the inner planet (larger eccentricities) and the outer one (shorter eccentricities; bottom). In the top panel, the solid lines represent the inner planet and the dashed lines the outer planet.

In the text

Fig. 3 Time evolution of the planets’ orbital properties. From top to bottom: the planets’ semi-major axes, the orbital period ratio of the outer planet to the inner one, and the eccentricity of the inner and outer planets. In the top panel, the solid lines represent the inner planet and the dashed lines the outer planet. The two planets are initially located close to a 2:1 resonance. The horizontal dashed line, r = 0.22, corresponds to the inner damping radius. The vertical dark and light grey lines represent the time when the planets are in resonance. In our fiducial model, the dark grey line corresponds to t = 4840 yr and the light grey line to t = 11 660 yr.

In the text

	Fig. 4 Phase space of the inner planet for the fiducial model. The colour indicates the time.
In the text

	Fig. 5 Azimuthally averaged surface density perturbations before, during, and after resonance capture. The dotted dark grey line shows the average surface density perturbation at t = 3000 yr, the solid blue at t = 7500 yr, and the dashed light grey at t = 12 500 yr.
In the text

Fig. 6 Influence of viscosity on the evolution of the planetary system for αv = 10−6,10−5,10−4, and 10−3, where 10−5 represents our fiducial model. The figure displays the semi-major axes, the period ratio of planets, and the eccentricity of inner and outer planets from top to bottom. The first row is divided into right and left panels, illustrating the semi-major axes of the inner and outer planets, respectively.

In the text

	Fig. 7 Period ratio and eccentricity evolution of planets for the model with α = 10−3. To track the evolution of the α = 10−3 model, we extended the simulation time up to 45 000 yr. The presence of a limit cycle is evident in the eccentricity evolution, denoted by persisting eccentricity oscillations without damping to zero.
In the text

	Fig. 8 Phase space of the inner planet for the α = 10−3 model. The colour indicates the evolutionary time. The presence of a limit cycle is clearly evident in the phase space, denoted by the repeating oscillations of the inner planet’s orbit around a stable equilibrium value.
In the text

Fig. 9 Impact of a decreasing surface density of the disc on the orbital properties of planets. From top to bottom: evolution of planets semi-major axes, the period ratio of planets, and the eccentricity of the inner and outer planets. The dashed lines represent the outer planet, while the solid lines represent the inner planet. The dotted line, which corresponds to r = 0.22, represents the inner damping radius.

In the text

	Fig. 10 Impact of an increasing surface density of the disc on the orbital properties of planets. From top to bottom: evolution of planets semi-major axes, period ratio of planets, and eccentricity of the inner and outer planets. The dashed lines represent the outer planet, while the solid lines represent the inner planet. The dotted line, r = 0.22, represents the inner damping radius.
In the text

	Fig. 11 Impact of increasing the aspect ratio on the evolution of the planetary system. From top to bottom: evolution of planets semi-major axes, the period ratio of planets, and the eccentricity of the planets. The solid lines represent the inner planet, while the dashed lines represent the outer planet.
In the text

	Fig. 12 Influence of a decreasing aspect ratio on the evolution of the planetary system. The first panel shows the evolution of the semi-major axis and the second panel the period ratio of the planets.
In the text

	Fig. 13 Effect of resolution on overstability occurrence.
In the text

	Fig. 14 Time evolution of the eccentricity (top) and period ratio (bottom) of the two planets for various configurations of PLUTO, compared to the fiducial FARGO3D setup presented in Table 1. The two planets do not break out of resonance in the fiducial PLUTO setup, but a higher resolution or a third-order reconstruction (WENO) both recover the behaviour found by FARG03D.
In the text

	Fig. 15 Time evolution of resonant angle for the fiducial model. The colour represents the eccentricity of the inner planet. The horizontal line represents the equilibrium value of the resonant angle during the resonance period of planets, which is ϕeq = 22.54°.
In the text

Fig. 16 Top to bottom: equilibrium values for the resonant angle, eccentricities of the inner and outer planets, and the ratio of the planet’s eccentricity for the models that captured them in a 2:1 MMR. The viscosity models, initial surface density models, and aspect ratio models are depicted from left to right, separated by vertical lines. The square, diamond, cross, and circle markers represent the fiducial model, the overstable models, the non-overstable model captured in the 2:1 MMR, and the limit cycle model, respectively. Additionally, Σ1 = Σ(r0).

In the text

	Fig. 17 Time evolution of gaps’ positions and depths. Top: time evolution of the semi-major axis and the position at which surface density has its minimum value. Bottom: time evolution of the surface density in the absence of the planet (disc’s power-law surface density) in blue, and the surface density in the position of planets in dark and light red for inner and outer planets.
In the text

	Fig. 18 Evolution of the surface density of gaps and its perturbation in terms of positions where surface density has its minimum value. The colour corresponds to time.
In the text

Fig. 19 Azimuthal average of disc surface density for multiple models at t = 3000 yr. We include our fiducial model (FM), two other models with single planets (one with a low-mass planet and the other with a large-mass planet), and another model identical to the FM but with zero eccentricity. The orange circles represent planets in the single-planet models, while the red circles depict planets in the FM.

In the text

	Fig. 20 Evolution of gap surface density and its perturbations in terms of the positions of the gap bottoms for the Σ(r0) = Σ0/4 model. The stars indicate the times when planets enter resonance.
In the text

	Fig. 21 Snapshots of the disc’s surface density over time at 2000, 20 000, and 30 000 yr for the h0 = 0.03 model. The position of gaps is marked by gold circles, and planets by green circles.
In the text

Fig. 22 Comparison of semi-major axes obtained through simulation with different theoretical models. Top: orbital evolution of our FM, along with two semi-theoretical models from Paardekooper et al. (2010, without a gap opening) and Kanagawa et al. (2018, with a gap opening), alongside two additional simulations using the same disc properties as the FM but focussing on individual planets. Solid lines represent the semi-major axis of the inner planet, while dashed lines represent the outer planet. Bottom: evolution of the quantity (ap – asim)/ap, where ap represents the semi-major axis according to pure type-I migration and asim denotes the semi-major axis of both the inner and outer planets in our FM, as well as the single small and big planets in our single-planet simulations. The period of resonance is indicated between the dark and light grey time intervals.

In the text

	Fig. 23 Same as the top panel of Fig. 22 but for the Σ(r0) = Σ0/4 model.
In the text

	Fig. 24 Same as the top panel of Fig. 22 but for the h0 = 0.03 model.
In the text

	Fig. A.1 Impact of the disc containing opened gaps on the occurrence of overstability. The solid lines show the models that the evolution of the system is together with gap formation. The dotted lines illustrate the impact of the disc with opened gaps on the system.
In the text

	Fig. B.1 Resonance angle evolution for the $\text{Σ}\left(r_0\right)=\frac{1}{3}{\text{Σ}}_0$ model. The colour indicates the eccentricity of the inner planet.
In the text

	Fig. B.2 Resonance angle evolution for the $\text{Σ}\left(r_0\right)=\frac{1}{4}{\text{Σ}}_0$ model. The colour indicates the eccentricity of the inner planet.
In the text