© The Authors 2023

1 Introduction

The migration history of Jupiter and Saturn played a significant role in the shaping of our Solar System (Walsh et al. 2011; Tsiganis et al. 2005). To date, there are 163 giant exoplanets with masses between Saturn and five Jupiter masses and an orbital period larger than 100 days. Within this population of so-called warm Jupiters, there are 43 systems composed of at least two giant planets¹. Studying the dynamical history of such systems during their formation is a key point in understanding the variety of architectures of planetary systems (see Zhu & Dong 2021, for a review).

In the specific case of our Solar System, Masset & Snellgrove (2001) found that the progenitors of Jupiter and Saturn may have undergone convergent migration until Saturn was captured in the 3:2 mean motion resonance (MMR) with Jupiter. In such a compact orbital configuration, the two planets are observed to revert their direction of migration and start migrating outwards while remaining locked in the 3:2 resonance. The work of Masset and Snellgrove was extended by Morbidelli & Crida (2007), who explored a wider set of initial conditions and disc parameters.

Focussing on the specific case of our Solar System, Walsh et al. (2011) suggested that the Jupiter-Saturn pair did reverse their migration and that this phase of outward migration started when Jupiter was at a distance of about 1.5 au from the Sun. This scenario, labelled Grand Tack, would explain both the small mass of Mars with respect to the Earth and the variation of composition in the asteroid belt, which are key points in explaining the final architecture of our Solar System.

The reversal of migration, commonly called the Masset and Snellgrove mechanism, requires the following:

• (i)

  The mass of the interior planet must exceed that of the exterior one;

• (ii)

  The two planets must share a common gap so that the outer (negative) Lindblad torque exerted on the less massive planet does not overcome the (positive) inner Lindblad torque exerted on the more massive one;

• (iii)

  The orbital configuration of the pair of planets must be compact such that their horseshoe regions are close enough to overlap. In this way the gas that enters the horseshoe region outside the outer planet’s orbit may be transported in the disc region inside the orbit of the inner planet, passing through its horseshoe region. The pileup of gas at the inner edge of the gap then maintains the positive inner Lindblad torque that sustains the outward migration phase.

These three conditions are often met when considering a Jupiter-Saturn pair locked in 3:2 MMR. Furthermore, Crida et al. (2009b) showed that the outward migration can then be maintained on a very long range.

D’Angelo & Marzari (2012) also addressed the question of the migration of Jupiter and Saturn in a more complex thermo-dynamical model, where the gas disc is evolving and planets are accreting gas. They showed that there can be severe limitations to the outward migration conditions: the relative speed of the pair of planets may not be large enough for the exterior planet to cross the 2:1 resonance. In such a configuration, the conditions (ii) and (iii) are not satisfied. The authors also showed that the specific requirement on the planets’ masses (i.e. condition (i)) may easily be violated by planetary accretion. More recently, Chametla et al. (2020) have suggested that the conditions for the Grand Tack model depend on the precise initial orbit of Jupiter and Saturn and on their formation timescales. They find that capture in the 2:1 MMR and slow inward migration is favoured if Saturn enters in resonance with Jupiter before being fully grown, at about two-thirds of its final mass.

At odds with these results, Pierens et al. (2014) found that in relatively small mass, viscously heated discs (similar model to D’Angelo & Marzari 2012), Jupiter and Saturn may migrate outwards when locked in the 2:1 MMR. However, this result was obtained for viscosity values lower than the one used in D’Angelo & Marzari (2012). This implies that discs are characterised by an unusually small aspect ratio of the order of 0.02–0.03. As the authors noticed, discs are not expected to be so cold when adding stellar irradiation (although cold stellar irradiated discs have been recently discussed in Savvidou & Bitsch 2021), so that outward migration in the 2:1 MMR cannot be maintained on a very long range.

At the time of these cited studies, discs were believed to have a non-negligible viscosity to justify the high accretion rates of gas onto the central star. Viscosity was raised by the turbulence generated via the magneto rotational instability (MRI, Balbus & Hawley 1991). Therefore, the previously cited literature concerning the migration of a pair of giant planets mainly focusses on discs with a value for the α turbulence parameter (Shakura & Sunyaev 1973) in the interval [10⁻³, 10⁻²].

However, it has recently been shown, both theoretically and observationally, that protoplanetary discs are probably less viscous than previously thought. On the theoretical side, for example, it has been shown that the gas near the mid-plane has a level of ionisation that is too weak to sustain the MRI (Gammie 1996). And more recently, the inclusion of non-ideal magneto-hydrodynamics (MHD) effects, such as ambipolar diffusion, led to the conclusion that disc should not be turbulent even at the surface (see Turner et al. 2014 for a review). On the observational side, Pinte et al. (2016) derived a value of the turbulence parameter on HL Tau of the order of 10⁻⁴ by modelling dust settling. With the same method, Villenave et al. (2022) found a value of the order of 10⁻⁵. Another example of turbulent parameter <10⁻³ comes from observations of the Lupus and Upper Sco star-forming regions (Sellek et al. 2020), where such low viscosity values are required in order to reproduce the range of masses and accretion rates seen.

In this context, we consider it important to revisit the migration of pairs of giant planets in low-viscosity discs. In this paper, we consider two-dimensional discs. We are aware of the fact that low-viscosity discs do not provide the sufficient transport of gas required for a typical accretion rate observed in young stars (Hartmann et al. 1998; Manara et al. 2016). Magnetically driven disc winds have been proposed as a mechanism to remove angular momentum from thin ionised surface layers of low-viscosity protoplanetary discs (Suzuki & Inutsuka 2009; Bai & Stone 2013; Turner et al. 2014; Béthune et al. 2017). Angular momentum removal promotes, in these layers, a fast radial transport of gas towards the central star. We will mimic the effect of magnetically driven winds in three-dimensional hydrodynamical simulations (following the model in Lega et al. 2022) in a forthcoming paper. In this paper, we prefer to explore the dynamics in two-dimensional simulations, which have the advantage of providing a wide parameter exploration over long timescales in reasonable computation time.

The layout of this paper is as follows. We first present, in Sect. 2, how we conducted our study, including the description of the physical model and the numerical setup, as well as the structure of our simulations. We then present our results in Sect. 3, with each subsection detailing the parameter exploration of this study. In Sect. 4.1, we explain why Saturn never crosses the 2:1 MMR in low-viscosity discs using a semi-analytical method. We confirm the robustness of these results for a wider range of planetary masses in Sect. 4.2. Finally, we discuss our work and conclude in Sect. 5.

2 Model setup

2.1 Physical model

In this work, we consider a protoplanetary non-self-gravitating disc made of gas only. The gas follows the vertically integrated Navier-Stokes equations (Masset 2002). Similar to Lega et al. (2021), we used an adiabatic equation of state to which we added an exponential damping of the temperature perturbations. This damping was done on the timescale of the cooling time τ_c. Following Lega et al. (2021), we considered τ_c to be equal to one local orbital period.

The two-dimensional system is described in polar coordinates (r, ϕ), centred on the star, where r is the radial coordinate and ϕ the azimuthal coordinate. We adopted a flared aspect ratio given by Chiang & Goldreich (1997) such that h = h₀(r/r₀)^2/7 with h₀ = 0.05, unless otherwise stated, and r₀ is the unit length (see Sect. 2.2). The disc density is given by Σ = Σ₀(r/r₀)^−1/2 where Σ₀ = 6.76 × 10⁻⁴ in code units for the nominal case. This corresponds to 222 g cm⁻² for r₀ = 5.2 au and a Sun mass star. We used the a prescription (Shakura & Sunyaev 1973) for the kinematic viscosity, such that v = αc_sH, where c_s is the sound speed and H the scale height of the disc such that H = hr. Throughout this paper we use α = 10⁻⁴. This value is observationally and theoretically motivated as mentioned in the Introduction. Moreover from a numerical point of view, it has been shown that numerical convergence cannot be obtained for values of α lower than 5 × 10⁻⁵ (McNally et al. 2019). Numerical convergence with respect to the viscosity is discussed in Appendix A.

Two planets were added in the system at a different time. The inner planet is a Jupiter-mass planet with mass M_J/M_⋆ = q_J = 10⁻³. The outer planet, unless otherwise stated, is a Saturn-mass planet with mass M_S/M_⋆ = q_S = 2.9 × 10⁻⁴. For simplicity, the two planets are henceforth referred to as Jupiter and Saturn, respectively.

2.2 Units

The code units are such that G = M_⋆ = r₀ = 1, where r₀ is the unit length and the initial semi-major axis of Jupiter. Therefore, the orbital period at r₀ = 1 is 2π time code units. Planetary masses, M_p, are then given as a ratio to the stellar mass q_p = M_p/M_⋆.

We remark that the Navier–Stokes equations are invariant with respect to the mass of the system, even in the case of an adiabatic equation of state with our prescribed cooling (see Sect. 2.1), since the energy is directly proportional to the density. Only the mass ratios matter; it would not be the case if the cooling was given by the opacity of the disc, which depends on the amount of dust in physical units.

For example the three following sets of physical parameters correspond to one single simulation:

• (i)

  M_⋆ = M_⊙, M_p = M_J, ∑₀ = 1000 g cm⁻², r₀ = 5 au;

• (ii)

  M_⋆ = 2M_⊙, M_p = 2M_J, Σ₀ = 2000 g cm⁻², r₀ = 5 au;

• (iii)

  M_⋆ = M_⊙, M_p = M_J, Σ₀ = 250 g cm⁻², r₀ = 2.5 au.

Indeed, in the three cases, M_⋆ = 1, M_p = 10⁻³, and r₀ = 1 are in code units and therefore ${\text{Σ}}_0=3\times {10}^{-3}{M}_{\star }{r}_0^{-2}$. In the case of (ii), the mass unit is doubled with respect to case (i), while the distance unit is unchanged. In the case of (iii), the mass unit is unchanged while the distance unit is halved, such that the surface density is divided by four².

In other words, only the dimensionless quantities q, α, h₀ and the disc-to-star mass ratio µ = πΣ₀r₀²/M_✶ (see also Crida & Morbidelli 2007) are relevant to characterise a simulation. Similarly, the initial position of Saturn matters only in terms of the ratio to the semi-major axis of Jupiter. Thus, these are the parameters given in Table 1. The reader can then interpret the result choosing their favourite stellar mass and orbital dis-tance³. To give some sense of the timelines of the phenomenon that we observed in the frame of planetary formation, in this paper, we provide numbers in physical units based on M_✶ = M_⊙ and r₀ = 5.2 au, so that one orbital period at r₀ corresponds to 11.85 yr.

2.3 Numerical setup

The 2D simulations were carried out with the code FARGOCA, the code FARGO (Masset 2000) with Co-latitude Added (Lega et al. 2014)⁴. The code was parallelised using a hybrid combination of MPI between the nodes and OpenMP on shared memory multi-core processors.

The radial domain extends in the range r ∈ [0.2, 9]r₀, which corresponds to r ∈ [1.04, 46.8] au for r₀ = 5.2 au. The mesh was divided in N_rad = 568 radial cells with logarithmic spacing and N_sec = 940 azimuthal cells. The resolution is such that dr/r = dθ = 0.0067, so that a pressure scale length H is covered by about five cells at the inner edge of the grid. This resolution is proven to be sufficient in Appendix A. The boundary conditions used in the radial direction follow the prescription of the evanescent boundary condition (De Val-Borro et al. 2006).

We worked in the reference frame co-rotating with Jupiter so that Jupiter laid on the x-axis at all times. We recall that the system is centred on the star and, therefore, the star is accelerated by the gravity of the planets and the disc. Indirect forces thus arise and must be carefully taken into account. The planets feel the indirect forces that are due to their own gravity as well as that of the disc. The disc feels indirect forces from the planets. Since the self-gravity of the disc has not been accounted for, we did not account for the indirect forces of the gas onto itself (Crida et al. 2022, and in prep.). The gravitational potential of the planets on the disc was modelled as follows: $${\psi }_{\text{p}}=-\frac{G{M}_{\text{p}}}{\sqrt{d^2+{ϵ}^2}},$$(1)

where d is the distance between the planet and a grid cell centre. The smoothing parameter ϵ used to evaluate the potential in the vicinity of the planets scales with the scale height of the disc H such that ϵ = 0.6H as is usually done in 2D simulations (e.g. Crida et al. 2009a).

2.4 Simulation setup

The simulations we present in this paper proceeded in two phases. During the first phase, we let Jupiter evolve alone in the disc. We started by introducing the planet on a fixed circular orbit at a distance a_J = r₀ from the star, with the remaining orbital elements being initialised at zero. Its mass was increased analytically⁵ from zero up to q_J in 800 local orbits, following $$q\left(t\right)=q_{\text{J}}{\sin }^2\left(\frac{\pi }{2}\frac{t}{T_{\text{growth}}}\right),$$(2)

where $T_{\text{growth}}=2\pi \times 800{\left(\frac{a_{\text{J}}}{r_0}\right)}^{3/2}$ code time units. We waited an additional 400 orbits to let the disc stabilise before allowing Jupiter to migrate. This process avoids triggering instabilities that could arise from introducing a Jupiter mass planet in the disc directly (Hammer et al. 2017; Hallam & Paardekooper 2020). If Jupiter was allowed to grow and migrate at the same time, it would have a fast inward migration episode before reaching the type II regime. This would shift the positions of the planets to the inside of the disc and as explained in Sect. 2.2 the results would only have to be re-scaled. The initial phase of the migration of Jupiter alone in the nominal disc is shown in Fig. 1.

The second phase started when introducing the outer planet in the disc after 4000 initial orbital periods of Jupiter, corresponding to about 47 500 yr. We call this time T_0‚S. Saturn was also initialised on a circular orbit at a distance a_S determined by the ratio ${\left[\frac{a_{\text{S}}}{a_{\text{J}}}\right]}_{T_{0,\text{S}}}$. Adding the outer planet at this moment was motivated by the fact that planet formation is favoured to take place at pressure bumps located at the edge of planetary gaps (Eriksson et al. 2021, and references therein). We then increased the outer planet mass from zero to Saturn’s mass, with the same method as in Eq. (2) for 800 initial orbital periods of Saturn (corresponding to about 25 200 yr in our nominal case, ${\left[\frac{a_{\text{S}}}{a_{\text{J}}}\right]}_{T_{0,\text{S}}}=2$).

The planet was free to migrate during its growth. This means that the outer planet goes through a phase of type I migration before it reaches a mass significant enough to open a gap. This allowed the planet to rapidly migrate and possibly be caught in resonance with Jupiter. The case of Saturn growing on a fixed orbit and then being allowed to migrate has been studied, but it resulted in both planets opening separate gaps which do not overlap. The planets’ evolution is independent with very slow migrations.

Fig. 1 Migration of Jupiter alone in a low-viscosity disc, corresponding to the parameters of the nominal simulation N. The top panel shows the evolution of the semi-major axis of the planet and the bottom panel is its eccentricity. The blue Saturn-shaped marker in the top panel indicates the moment at which Saturn was added in the simulations presented in Sect. 3.

3 Simulation results

In the following section, we present the results of the migration of the pair of giant planets for different disc parameters and separations. Table 1 gathers the simulations that were run for this study, the parameters that were explored, and a summary of their results.

3.1 Nominal simulation

In our nominal simulation, Saturn was positioned outside of Jupiter’s gap at a distance of a_S/a_J = 2. Panel a of Fig. 2 shows the surface density of the disc at the moment of introduction of Saturn. The other panels of Fig. 2 show the positions of Jupiter and Saturn in the disc at different moments of the evolution of the system. Panel b shows Saturn, at 45% of its final mass, transitioning from a wake-driven migration to carving its own gap. The gap was fully opened in panel c once the planet was close to its final mass. At last, panel d of Fig. 2 shows the disc surface density after 118 000 yr, when Jupiter and Saturn were migrating together in a common gap.

We show the evolution of the semi-major axes and eccentricities of the planets in Fig. 3. Saturn reached its final mass around 72 600 yr. The four coloured circles in the figure represent the moments at which the surface density of the disc is shown in the four panels of Fig. 2. Figure 3 shows that while growing, Saturn migrated slowly inwards until the pair of planets were locked in the 2:1 MMR. Figure 4 confirms that the two resonant angles of the system are librating and therefore the planets are stably locked into the resonance. These angles are defined as $${\theta }_{\text{J,S}}=2{\lambda }_{\text{S}}-{\lambda }_{\text{J}}-{\varpi }_{\text{J,S}}$$

with λ_p and ϖ_p being the mean longitude and longitude of perihelion of planet p, respectively.

When the planets entered into resonance, the eccentricity of Jupiter was excited up to the order of ~0.1 while that of Saturn remained very low. This is consistent with θ_S having a large libration amplitude, while θ_J remains very close to 0; however, this is an uncommon result as shown by Michtchenko et al. (2008) regarding the stability of such a configuration. Furthermore, until now in the literature, Jupiter-Saturn migration studies had found that in most cases Saturn’s eccentricity is higher than that of Jupiter, and this has been seen for the planets locked in the 3:2 MMR (see e.g. Masset & Snellgrove 2001; D’Angelo & Marzari 2012) as well as in 2:1 MMR (Pierens et al. 2014). Readers can refer to Sect. 3.3 for further information about the eccentricities of the planets.

3.2 Disc thickness

In classical discs, for α ≥ 10⁻³, it has been shown that the disc thickness can influence the direction of migration of a pair of Jupiter–Saturn-like planets when they are locked in a 3:2 MMR, (Morbidelli & Crida 2007). Therefore, in this section, we show how the result obtained for the nominal simulation varies with the disc thickness.

In Fig. 5, we show the evolution of the orbital parameters of both planets for simulations with a colder, thinner disc, C with h₀ = 0.035 and with a hotter, thicker disc, H with h₀ = 0.06. In the thicker disc, as in the nominal case, the planets are locked into a 2:1 MMR. While in the thinner disc, the planets end up in a 5:2 MMR configuration.

We shall focus on case C and its unusual outcome. We see in Fig. 5 that while growing, Saturn migrates slowly inwards until the pair of planets gets locked in a 5:2 MMR. Figure 6 shows the libration of the resonant angles of the system and confirms that the planets are stably locked into a 5:2 resonance.

To understand further why Saturn’s migration slows down before the 5:2 MMR and is captured into this resonance, we show in Fig. 7 the evolution of the azimuthally averaged density profiles of the disc during the phase of inward migration and growth of Saturn. The profiles are centred on Jupiter’s position in order to highlight Saturn’s motion with respect to the 5:2 MMR. The position of Jupiter and Saturn are marked with filled and empty circles, respectively. The shape of the gap carved by Jupiter depends on the disc’s aspect ratio (Crida et al. 2006). Indeed, the thinner the disc, the steeper the density slope is at the edge of the gap and the wider the gap is. This can result in the creation of a density bump at the outer edge of Jupiter’s gap, beyond the 2:1 MMR. Besides, for a fixed α, the viscosity is proportional to h² so that it is twice as small in the C case than in the nominal N case. As a consequence, the bump is less quickly and less easily smoothed and eroded. While in case N, we did not observe any bump, in the thin disc case C there is a clear bump that can be seen from the lightest green curve in Fig. 7 corresponding to the disc profile at time T_0,S. Figure 7 shows that while migrating in the type I regime, Saturn gets trapped at the density maximum located at the edge of Jupiter’s gap. This density bump acts as a ‘planet trap’ where the steep positive gradient in density leads to a positive co-rotation torque that cancels out the negative total Lindblad torque (Masset et al. 2006). Once trapped, the planet continues to grow until it starts forming its own gap. The trapping of Saturn at the density peak slows down its migration considerably, allowing it to get locked in the 5:2 MMR. Finally, we note that the eccentricities of both planets are much higher in the thin disc case.

We find that the disc’s aspect ratio has an influence on the resonant configuration of the pair of planets, through the initial gap profile of Jupiter. The planet trap at the edge of Jupiter’s gap will determine whether Saturn will be locked in the 5:2 resonance or migrate further inwards towards the 2:1 MMR. We therefore expect that there is a threshold for the disc height below which the planets will end up being in the 5:2 MMR, or in the 2:1 MMR otherwise.

The overall behaviour of the pair of planets, however, does not differ significantly whether they are in a 5:2 or 2:1 resonance. Precisely, we notice the direction of the migration of the pair is inwards in both resonant configurations. In both cases the orbital configuration is not compact enough for the reversal of the migration direction to occur as is typically observed in a 3:2 MMR. A comment on the inward migration speed is subsequently provided in the discussion (see Sect. 5).

Fig. 2 Surface density of the disc at different times of the nominal simulation. The filled and empty circles mark the positions of Jupiter and Saturn, respectively. Panel a shows the disc at the moment of introduction of Saturn, T0‚S, while the outer planet still has zero mass. In panel b, the planet is at 45% of its total mass. In panel c, the planet is just about to reach its final mass. Panel d shows the disc and the planets closer to the end of the simulation. At this point both planets are in a large common gap and locked in the 2:1 MMR. The colours given to the planets correspond to the markers in Fig. 3.

Fig. 3 Orbital parameters’ evolution of Jupiter and Saturn in the nominal simulation. The top panel shows the semi-major axes of Jupiter in red and Saturn in blue. The shaded areas mark the positions of the peri-and apo-centre, q = a(1 − e) and Q = a(1 + e), of the planets, respectively. The black dashed line marks the position of the 2:1 resonance with Jupiter. The bottom panel shows the eccentricities of the planets. At the start of these curves, Jupiter is fully grown, while Saturn reaches its final mass around 72 700 yr. The coloured markers indicate the times at which the surface density is represented in the four panels of Fig. 2.

Fig. 4 Evolution of the resonant angles θJ‚S of both planets in the nominal simulation. The angle λp is the mean longitude of planet p and ϖp is its longitude for the perihelion.

Fig. 5 Orbital evolution of Jupiter and Saturn for the simulations C, N, and H (see Table 1). The top panel shows the semi-major axis of both planets (Jupiter in red, Saturn in blue), and the bottom panel shows their eccentricities. The markers indicate the time at which the planets get locked into a resonance: the triangle marks the 5:2 MMR, and the diamond marks the 2:1 MMR.

Fig. 6 Evolution of the four resonant angles in simulation C.

Fig. 7 Radial density profiles of simulation C in the time interval [45 500, 60 000] yr. The green colour gradient corresponds to a time evolution (and therefore also a mass evolution for Saturn) from light to dark. The filled and empty circles indicate the positions of Jupiter and Saturn.

3.3 Disc mass

We are now interested in studying the effects of the disc mass on the migration of the pair of planets. We have run the simulations N_Σ− and N_∑+ from Table 1 and have found that just as for the nominal simulation, the pair gets locked into a 2:1 MMR. We therefore conclude that the result obtained for the nominal simulation is robust when the disc mass changes. We are however interested in exploring the result obtained in the previous section for simulation C.

Therefore, we explored different disc masses in the case of the thin disc, and studied how robust the 5:2 MMR outcome is. Figure 8 shows the evolution of the orbital parameters comparing the thin disc with nominal mass C with a less massive, C_∑−, and a more massive disc, C_∑+. Since the shape and depth of Jupiter’s gap does not depend on the surface density, in all three cases Saturn starts at a similar position with respect to the local density distribution. In the more massive disc, C_∑+, Saturn’s migration starts similarly as in C and the planet initially gets captured in the 5:2 resonance before 60 000 yr. However, because the disc is more massive, the force it applies to the pair of planets is strong enough to force Saturn out of the 5:2 resonance and make the planet continue its migration inwards until it reaches the 2:1 resonance. Once in this resonance, the pair of planets migrates inwards, significantly faster than in the nominal disc mass case.

Since type I migration is proportional to the disc mass (Paardekooper et al. 2011), in the less massive disc C_Σ−, the migration speed of the Saturn embryo is slower than in case C. As a consequence, Saturn opens its gap before reaching the planet trap. Instead, the planet enters a regime of type II migration and crosses the 5:2 MMR unperturbed. This is shown in Fig. 9.

In the C_Σ− case, the eccentricities of Jupiter and Saturn both oscillate around 10%. This indicates that the high eccentricity of Jupiter in the other simulations is induced by the gas’ potential in the surroundings of the planets. We have tested this hypothesis in Appendix B.

Fig. 8 Same as Fig. 5, but for simulations CΣ−, C, and C∑+.

Fig. 9 Same as in Fig. 7, but for the time interval [45 500, 83 000] yr, for simulation CΣ−.

3.4 Distance between planets

The dependence of the outcome of the migration of Jupiter and Saturn on their initial distance separation has been shown in Chametla et al. (2020). In this section, we investigate how the initial position of Saturn with respect to the pressure bump present in thin disc case C may affect the final MMR configuration of the pair. As there is no pressure bump in nominal case N, we expect the planets to end in the 2:1 MMR in any case (provided it starts outside of the 2:1 and not too far away) and therefore we only explore this parameter in the cold disc case.

We ran simulations in which we placed Saturn at a distance a_S/a_J = 2.5 and a_S/a_J = 3. As explained in Sect. 2.2, the important quantity is the ratio of the semi-major axes (or equivalently of the orbital periods), not the difference⁶. We recall that the growth of Saturn is set to last 800 initial orbital periods, which entails the growth time of the planet being longer the further away it is placed in the disc. The consequence of this is that the planet remains longer (in physical time, not in orbital periods) in a fast migration regime compared to the nominal case. If the growth time is shorter, Saturn carves its gap before reaching Jupiter and both planets migrate independently. Figure 10 shows the results of these simulations.

Case C_2.5 is very similar to the C simulation; both planets are locked in the 5:2 resonance and their final eccentricities reach the same values. In the C₃ case, the pair of planets end in the 2:1 configuration; this is due to Saturn experiencing a fast inward migration episode between 75 000 and 80000 yr. Figure 11 shows the density profiles of this simulation at the time of this episode. The planet slowly migrates inwards within the low mass migration regime until it starts carving its gap near the edge of Jupiter’s gap. This pushes gas between the two planets creating an even stronger density peak than in the case of simulation C. That steep density gradient combined with a co-orbital mass deficit has positive feedback on the inward migration of Saturn, triggering a short episode of runaway migration. As a consequence, Saturn crosses the 5:2 MMR and continues its migration path until it encounters the 2:1 MMR. We recall that Saturn does not accrete gas from the disc while it grows. If this had been the case, this could have a consequence on the size of the density peak located between the two planets and therefore on the episode of the fast migration of Saturn.

Finally, if Saturn is formed further away in the disc (a_S ≳ 5a_J), it would not migrate fast enough to catch up to Jupiter and both planets would have independent migrations. We also assess in Appendix C the case in which Saturn would form on the inside of the 2:1 MMR, at an initial distance of a_S/a_J = 1.4. In this case only, we find Saturn and Jupiter locked in a 3:2 MMR and the pair migrates outwards. However this scenario requires that Saturn form inside Jupiter’s gap and this is therefore unlikely.

Fig. 10 Same as Fig. 5, but for simulations C, C2.5, and C3. Simulation C2.5 was stopped at about 110 000 yr as it overlaps almost perfectly with C in both orbital parameters.

Fig. 11 Same as in Fig. 7, but for simulation C3 in the time interval [75 000, 80 000] yr. It is important to note the jump of Saturn between masses 68% (r = 2.12aJ) and 71% (r = 1.76aJ) in 100 orbits at r0.

4 Why the 2:1 MMR is never crossed

We have seen so far that in low-viscosity discs, it seems that Saturn never crosses the 2:1 resonance. Indeed the planets are never in the 3:2 MMR and as a consequence the pair never migrates outwards. This is in contrast with results found in the literature that consider an α viscosity parameter no less than 10⁻³ (e.g. Masset & Snellgrove 2001; Morbidelli & Crida 2007; D’Angelo & Marzari 2012). Therefore, a relevant question to assess in this work is why the 2:1 resonance is never crossed at low viscosity. In this section, we answer this question using a semi-analytical method and generalise this result.

4.1 Resonance crossing criterion

The crossing of a resonance is dictated by the relative migration speed between the two bodies and the libration period of their resonant angle. This has been studied by D’Angelo & Marzari (2012) and further developed in Batygin (2015). The later paper has shown analytically that the outer body gets locked into a resonance if the libration period of the resonant angle is shorter than the characteristic timescale of convergence. Evidently this criterion for resonance locking depends on the migration rate of both planets.

If we know the analytical migration rate of the planets, we can derive a fully analytical method to understand the resonance crossing. In Batygin (2015), Sect, 3.1.2 (i), the author derives the criterion specifically applied to an inner planet migrating in a type II regime and an outer, less massive planet following a type I migration (cf. their Eq. (50)). However, we find that this derivation is based on assumptions that are too restrictive to be applied to our case. Firstly, the author assumes that most of the planetary mass is contained within the inner body and therefore neglects the mass of the outer planet in his derivation. In the case of Jupiter and Saturn with a mass ratio of M_J/M_S ≃ 3, it can be argued that this assumption is no longer valid. But most importantly, this criterion was derived using the analytical expression of type I migration from Tanaka et al. (2002) for the outer planet. However, type I migration is valid only up to a certain planetary mass. This transition mass depends on the viscosity of the disc.

To illustrate this point, we ran simulations of a planet growing from zero mass to M_p = M_S in 800 initial orbital periods for two viscosity parameters α = 10⁻⁴ and 10⁻³ as well as the aspect ratio of h₀ = 0.035 and we studied its migration as a function of its mass. The top panel of Fig. 12 shows the migration speed of the planet in absolute value, in addition to which we have superposed the analytical form of type I migration from Tanaka et al. (2002, which we recall is independent of the viscosity) given for a fixed r = r₀. With this figure, we show that the planet in a low-viscosity disc departs from the type I migration regime before reaching 10 M_⊕ (while not before 40 M_⊕ for α = 10⁻³). This is in agreement with the literature (McNally et al. 2019; Fung & Chiang 2017). As a consequence, using this analytical expression as the migration rate of the outer planet is not applicable in our configuration and we are therefore unable to use Eq. (50) of Batygin (2015).

Alternatively, in order to understand the locking of Saturn and Jupiter in the 2:1 resonance in our simulations, we used the general form of the criterion derived in Batygin (2015, cf. their Eq. (44)). We now assume that the relative migration rate between both planets is dominated by Saturn’s migration, thus neglecting Jupiter’s migration rate. We rewrite the criterion for resonance crossing as a critical migration timescale below which the outer planet would cross the resonance k : k − 1, expressed as $${\text{Ω}}_{\text{S}}{\tau }_{crit}\simeq \frac{2\pi k/\left(k-1\right)}{4{\left(3\right)}^{2/3}{\left[k{\left(k-1\right)}^2\right]}^{1/9}}{\left(\frac{f_{res}^{\left(1\right)}\left(M_{\text{J}}+M_{\text{S}}\right)}{M_{\star }}\right)}^{-4/3}$$(3)

where $\tau =\frac{a}{\left|\dot{a}\right|}$ is normalised by the angular velocity of the outer planet, and $f_{res}^{\left(1\right)}$ is a dimensionless constant (that can be found in Deck et al. 2013). We note that this criterion applies only to first order resonances.

Since we do not have an analytical expression for the migration of a planet in a low-viscosity disc, we used the migration timescales obtained by the simulations presented in Fig. 12 and compared it with the critical migration timescale for resonances 2:1 and 3:2 calculated from Eq. (3). We show in the bottom panel of Fig. 12 that, in the case of α = 10⁻³, the minimal migration timescale reached by the planet is close enough to the critical value to cross the 2:1 resonance, while it remains far from that of the 3:2 resonance. This illustrates the well-known result obtained for Jupiter and Saturn migration in ’classical’ viscous discs, that is the pair locked in the 3:2 MMR. In the low-viscosity case, however, it is clear that the migration timescale remains much higher than the timescale required to cross the 2:1 MMR. This seems to indicate that a planet growing and migrating in such a disc never migrates fast enough to cross the 2:1 MMR.

In the simulation presented in Fig. 12, the planet grows in an unperturbed disc. In fact, our study concerns the case of Saturn migrating in a disc where Jupiter has already formed a gap. Therefore, we ran a simulation in which we let Saturn migrate in the disc perturbed by the presence of Jupiter, but without it feeling the gravitational potential of Jupiter⁷. This allowed us to take into account the presence of Jupiter’s gap in the migration of Saturn as this one grew. Figure 13 shows the relative migration timescale (that is the variation timescale of the ratio of the semi-major axes a_S/a_J) in the perturbed disc as a function of the ratio of semi-major axes between the two planets. This figure confirms that the relative migration timescale at the time when the planets approach the 2:1 MMR is much higher than the critical value to cross that resonance. We recall that, in this simulation, Saturn does not feel the presence of Jupiter and therefore does not remain locked in the said resonance. We also checked that the same is true for our nominal aspect ratio of h₀ = 0. 05.

Incidentally, we also note in Fig. 13 that Saturn stops at the planet trap just outside the 5:2 resonance, as has already been discussed in Fig. 7. When it leaves the planet trap because it starts opening a gap, it approaches the 5:2 resonance very slowly. In this simulation, Jupiter is temporarily caught in this resonance with Saturn (Saturn does not feel Jupiter, but Jupiter feels Saturn such that we recovered the same response in Jupiter’s eccentricity as in simulation C). However, the disc eventually pushes Saturn across the 5:2 resonance and convergent migration resumes.

Fig. 12 Migration speed and timescale of a growing planet in different viscous discs. Top panel: migration speed in absolute value of a planet growing from zero to Saturn’s mass, MS, in 800 initial orbital periods. The dashed black line shows the type I migration speed from Tanaka et al. (2002, at r = r0 constant). Bottom panel: migration timescale τ = a/|ȧ| normalised by the angular velocity of the planet. The black dashed and dotted lines represent the critical migration timescale, calculated with Eq. (3), beyond which the planet would cross the 2:1 or 3:2 resonance, respectively. These simulations were run for viscosities corresponding to α = 10−4 and 10−3. The migration timescales plotted in this figure were smoothed by means of a sliding window average with a size of then and 50 points (which represents a time interval of about 100 and 600 yr) for the case α = 10−3 and α = 10−4, respectively.

Fig. 13 Relative migration timescale $\tau =\frac{a_{\text{S}}/a_{\text{J}}}{d\left(a_{\text{S}}/a_{\text{J}}\right)/dt}$ in a simulation where Saturn migrates in the cold disc perturbed by Jupiter, without feeling Jupiter’s gravity. The normalised timescale is plotted as a function of aS/ aJ, and the vertical black dashed line marks the location of the 2:1 MMR, while the grey dotted vertical line marks the position of the 5:2 MMR. The horizontal dashed line marks the critical migration timescale, calculated from Eq. (3) (for Mp = Ms), above which the planets should not cross the 2:1 resonance. The migration timescale was smoothed by means of a sliding window average with a size of 100 points (which represents a time interval of about 12 000 yr), in order to suitably reduce the noise.

4.2 Changing the outer planet’s mass

From the analysis performed in Sect. 4.1, we show that the outer planet does not cross the 2:1 MMR in the case of a low-viscosity disc independently of its mass, since the peak migration speed is reached way before Saturn’s mass. In this section, we further expand our study by considering a range of masses for the outer planet. The inner planet remains a Jovian mass planet.

We show the results of our simulations in Fig. 14 for the following masses of the outer planet M_p ∈ [M_J, 2M_J, 3M_J] and h₀ = 0. 05. As expected in all of these cases, the outer planet never crosses the 2:1 MMR with the inner planet. We note, however, that some of these systems are not stable in the long term due to close encounters between the planets.

In classically viscous discs, the outcome of the migration of pairs of giant planets often required some conditions on the mass ratio of the planets (e.g. Masset and Snellgrove mechanism). This work shows that in a low-viscosity disc, the pair gets mostly locked in a 2:1 MMR independently of the mass of the outer planet and this does not require the inner planet to be more massive.

Fig. 14 Same as Fig. 5, but for the simulations $N_{M_{\text{J}}},N_{2M_{\text{J}}}$, and $N_{3M_{\text{J}}}$. We note that $N_{2M_{\text{J}}}$ becomes unstable shortly after 130 000 yr.

5 Discussion and conclusions

In this paper, we have explored the dynamics of pairs of giant planets in low-viscosity discs. We ran two-dimensional hydrody-namical simulations, which allowed us to perform an extensive parameter exploration over long timescales.

Our results show that the planets primarily get locked in a 2:1 MMR. In some cases, we have found the planets in the higher order resonance 5:2. In both cases, we do not observe outward migration. In fact, we find that in most cases the planets migrate inwards, with Jupiter migrating at a speed of the order of a few astronomical units per million years. More precisely, in nominal simulation N, we find that the Jupiter–Saturn pair migrates inwards at a similar speed as Jupiter alone, with a speed of the order of 6 au/10⁶ yr. While in the thin disc case, C, once in resonance, Jupiter migrates at a speed of 2 au/10⁶ yr; in comparison, the migration speed of Jupiter alone in the cold disc is about 10 au/10⁶ yr. Our results therefore show that having a pair of giant planets does not stop inward migration, but it can in some cases slow down the migration speed with respect to a single planet.

We have shown that the planets never cross the 2:1 MMR, which is at odds with respect to higher viscosity discs. In order to explain our results, we have used an analytical criterion for resonance crossing based on the work of Batygin (2015). We show that a planet growing and migrating in a low-viscosity disc never reaches the migration speed required to cross the 2:1 MMR with Jupiter.

The only case in which outward migration is observed (and therefore is compatible with the Grand Tack model) is a scenario where Saturn would form inside the 2:1 resonance and get locked in the 3:2 MMR. However, owing to the large width of Jupiter’s gap, this scenario, implying that Saturn forms inside this gap, seems unlikely. Thus we conclude the Grand Tack to be almost impossible with α ≲ 10⁻⁴. This does not rule out a low-viscosity proto-solar nebula, though, since alternative explanations have been proposed for the small mass of Mars and the depletion of the asteroid belt (e.g. Drążkowska et al. 2016; Clement et al. 2018, 2019; Nesvorný et al. 2021; Izidoro et al. 2022; Morbidelli et al. 2022), so the Grand Tack may well not have happened after all.

As for the Nice model, our results show that if the proto-solar nebula had α ≲ 10⁻⁴, the Saturn/Jupiter period ratio would not be less than two after the gas disc phase, in contrast to most of the situations studied in the literature (e.g. Tsiganis et al. 2005; Morbidelli et al. 2007). However, the case of Jupiter and Saturn in the 2:1 MMR has been explored by Nesvorný & Morbidelli (2012) and Deienno et al. (2017), who found that this configuration is not incompatible with a global instability of the giant planets that would bring them to their current orbits, under some conditions. More recently, Clement et al. (2021b,a) argue that Jupiter and Saturn in the 2:1 MMR with significant eccentricities (more than 0.05 and up to 0.1 and 0.25) are favourable initial conditions for a global ’Nice-model-like’ instability. Our results support such a setup, although more work is needed to connect our final configuration to their initial conditions (see for instance Appendix B).

Outside of the Solar System, we remark that the only system with two giant planets observed in a protoplanetary disc, to date, namely PDS 70, occurs to be close to the 2:1 resonance (Haffert et al. 2019). Bae et al. (2019) have shown numerically that a PDS 70-like system initially positioned in the 2:1 MMR, in a classical viscous disc, remains dynamically stable over million-year timescales. However, their study does not address how the two planets have reached such a resonant configuration. Our results suggest that a low viscosity may explain their capture. Simulating the specific case of a PDS 70-like system is beyond the scope of this paper and will be the topic of a future work.

We also remark that the two-dimensional low-viscosity discs that we have considered in this paper do not have an accretion rate onto the star compatible with observed values. In a forthcoming study, we will use a new paradigm for the theoretical modelling of discs where the accretion onto the central star occurs in superficial layers while the mid-plane has a close to zero viscosity (Lega et al. 2022). In such discs, the migration of a single giant planet differs from the classical type II migration regime and depends on the thickness of the accretion layer. It is therefore interesting to extend this study to a pair of giant planets.

Acknowledgements

We thank the reviewer for comments that helped in improving the manuscript. We thank Frédéric Masset, Gabriele Pichierri and Michiel Lambrechts for enriching discussions. We acknowledge support by DFG-ANR supported GEPARD project (ANR-18-CE92-0044 DFG: KL 650/31-1). We also acknowledge HPC resources from GENCI DARI no. A0120407233 and from “Mesocentre SIGAMM” hosted by Observatoire de la Côte d’Azur. We wish to thank Alain Miniussi for maintenance and re-factorisation of the code FARGOCA.

Appendix A Viscosity and resolution

In this section, we first show the convergence of our results with higher resolution and secondly we show that, at lower viscosity values, such convergence cannot be achieved. We focus on the first phase of our simulations, that is when Jupiter is migrating alone in the disc. We recall that, in this paper, the resolution is such that in the nominal simulation we have about 7.4 cells per scale height at r₀, corresponding to N_rad × N_sec = 568 × 940. We show in Fig. A.1 that the simulation remains unchanged with the following grids: N_rad × N_sec = 852 × 1400 and N_rad × N_sec = 1136 × 1880.

In the case of lower viscosity, however, we show in Fig. A.2 that at higher resolution the migration of a single Jupiter-like planet is perturbed by small-scale instabilities in the discs. This confirms the findings from McNally et al. (2019), claiming that numerical convergence cannot be obtained for values of viscosity of about v < 10⁻⁷, corresponding in our case h=0.035 to about α ≃ 3 × 10⁻⁵ (and α ≃ 4 × 10⁻⁵ in the case h = 0.05).

Fig. A.1 Migration of Jupiter in a disc with viscosity α = 10−4 for higher resolutions.

Fig. A.2 Migration of Jupiter in a disc with viscosity α = 10−5 for higher resolutions.

Appendix B Removing gas potential: Effect on the eccentricities

In many simulations, we have found that Jupiter’s eccentricity is much higher than that of Saturn. We have suggested that these eccentricity values are maintained by the gas’ potential.

We have tested this hypothesis by gradually removing the force F(t) of the gas onto the planets with a taper function $$F\left(t\right)=F\left(t_0\right)e^{-\frac{\left(t-t_0\right)}{\tau }},$$(B.1)

with τ being the removal timescale and t₀ a suitable starting time. We show as an example in Figs. B.1 and B.2 the evolution of orbital parameters of simulation set C and C₃. At time t₀ ~ 120000 yr, we restarted the simulation removing the force of the gas onto the planets over a timescale τ = 10⁴ orbits at r₀ = 5. 2AU. We observed that:

• i)

  the pair of planets remains locked in their initial resonant configuration,

• ii)

  the eccentricity of Jupiter reduces and that of Saturn increases. This is consistent with the results obtained for C_Σ− and more generally expected from N-body simulations,

• iii)

  the libration angles θ have a non-zero amplitude indicating that the system is no longer at the stable equilibrium point. We note that the centre of libration can also change (e.g. in the case of the 2:1 MMR), as it is known to happen once the disc dissipates (see for example Siegel & Fabrycky (2021)).

Points i) and iii) are illustrated in Fig. B.3.

Fig. B.1 Orbital evolution of Jupiter and Saturn for simulation set C. The top and bottom panels show the semi-major axis and the eccentricity, respectively. The vertical dotted lines indicate the time t0 at which we started the simulation with gas removal, namely ’C + Gas removal’. The nominal simulation was continued for t > t0 until a significant difference appeared with respect to the gas removal.

Fig. B.2 Same as Fig. B.1, but for simulation set C3.

Fig. B.3 Resonant angles corresponding to the simulations shown in Fig. B.1 in the top panel, and those corresponding to Fig. B.2 in the bottom panel. We find that in both cases of 2:1 or 5:2 MMR, the planets remain locked into resonance after the gas removal. Although, their resonant angles either librate with larger amplitude or no longer around 0.

Appendix C Migration in the 3:2 MMR

We have shown that a Saturn mass planet is not able to cross the 2:1 MMR. However we wonder what the outcome of the migration of the pair of planets would be if Saturn were to form on the inside of the 2:1 MMR, so a_S/a_J < 1.58. In this case, we expect Jupiter and Saturn to be trapped in the 3:2 MMR, but the direction of their migration is an open question. Since the Masset and Snellgrove mechanism is independent of the viscous drift, we expect it to still be present in the low-viscosity case. However, in low-viscous discs, the gaps carved by planets are significantly deeper and it is not obvious that gas can flow from the outer to the inner disc, allowing the condition for outward migration (see (iii) in Sect. 1) to be satisfied.

Taking the same procedure as for the previous simulations, we placed Saturn at a distance of a_S = 1.4a_J, for h = 0.035. We show the result of this simulation in Fig. C.1. The planet moved inwards until reaching the 3:2 MMR in which the pair was stably locked for the rest of the simulation. Once the planets were in this configuration, they migrated outwards together as was observed in classically viscous discs (Masset & Snellgrove 2001; Morbidelli & Crida 2007). We checked that gas flows from the outward to the inner disc by placing gas tracers at the edge of the common gap. We indeed see that even in the low-viscous case, gas flows from the outer to the inner disc allowing the planets to migrate outwards.

We therefore find that in low-viscosity discs, it is possible to find the Masset and Snellgrove mechanism resulting in the pair Jupiter-Saturn migrating outwards, but only when the pair is locked in the 3:2 MMR. To have such a configuration of the planets requires that Saturn forms at a distance less than 1.58a_J. This is not likely if Jupiter is already formed as that would require Saturn to form in the gap of Jupiter. Another possibility would be that both planets form at the same time, allowing them to grow close enough to be in the 3:2 resonance. However, we have not studied this case as this is out of the scope of this paper.

Fig. C.1 Migration of Jupiter and Saturn if Saturn was to form on the inside of the of 2:1 MMR. The top panel shows the evolution of the semi-major axes of Jupiter, in red, and Saturn, in blue. The dashed black and grey lines represent the position of the 2:1 and 3:2 MMR with Jupiter, respectively. The bottom panel shows the eccentricities of both planets.

References

All Tables

All Figures

	Fig. 1 Migration of Jupiter alone in a low-viscosity disc, corresponding to the parameters of the nominal simulation N. The top panel shows the evolution of the semi-major axis of the planet and the bottom panel is its eccentricity. The blue Saturn-shaped marker in the top panel indicates the moment at which Saturn was added in the simulations presented in Sect. 3.
In the text

Fig. 2 Surface density of the disc at different times of the nominal simulation. The filled and empty circles mark the positions of Jupiter and Saturn, respectively. Panel a shows the disc at the moment of introduction of Saturn, T0‚S, while the outer planet still has zero mass. In panel b, the planet is at 45% of its total mass. In panel c, the planet is just about to reach its final mass. Panel d shows the disc and the planets closer to the end of the simulation. At this point both planets are in a large common gap and locked in the 2:1 MMR. The colours given to the planets correspond to the markers in Fig. 3.

In the text

Fig. 3 Orbital parameters’ evolution of Jupiter and Saturn in the nominal simulation. The top panel shows the semi-major axes of Jupiter in red and Saturn in blue. The shaded areas mark the positions of the peri-and apo-centre, q = a(1 − e) and Q = a(1 + e), of the planets, respectively. The black dashed line marks the position of the 2:1 resonance with Jupiter. The bottom panel shows the eccentricities of the planets. At the start of these curves, Jupiter is fully grown, while Saturn reaches its final mass around 72 700 yr. The coloured markers indicate the times at which the surface density is represented in the four panels of Fig. 2.

In the text

	Fig. 4 Evolution of the resonant angles θJ‚S of both planets in the nominal simulation. The angle λp is the mean longitude of planet p and ϖp is its longitude for the perihelion.
In the text

	Fig. 5 Orbital evolution of Jupiter and Saturn for the simulations C, N, and H (see Table 1). The top panel shows the semi-major axis of both planets (Jupiter in red, Saturn in blue), and the bottom panel shows their eccentricities. The markers indicate the time at which the planets get locked into a resonance: the triangle marks the 5:2 MMR, and the diamond marks the 2:1 MMR.
In the text

	Fig. 6 Evolution of the four resonant angles in simulation C.
In the text

	Fig. 7 Radial density profiles of simulation C in the time interval [45 500, 60 000] yr. The green colour gradient corresponds to a time evolution (and therefore also a mass evolution for Saturn) from light to dark. The filled and empty circles indicate the positions of Jupiter and Saturn.
In the text

	Fig. 8 Same as Fig. 5, but for simulations CΣ−, C, and C∑+.
In the text

	Fig. 9 Same as in Fig. 7, but for the time interval [45 500, 83 000] yr, for simulation CΣ−.
In the text

	Fig. 10 Same as Fig. 5, but for simulations C, C2.5, and C3. Simulation C2.5 was stopped at about 110 000 yr as it overlaps almost perfectly with C in both orbital parameters.
In the text

	Fig. 11 Same as in Fig. 7, but for simulation C3 in the time interval [75 000, 80 000] yr. It is important to note the jump of Saturn between masses 68% (r = 2.12aJ) and 71% (r = 1.76aJ) in 100 orbits at r0.
In the text

Fig. 12 Migration speed and timescale of a growing planet in different viscous discs. Top panel: migration speed in absolute value of a planet growing from zero to Saturn’s mass, MS, in 800 initial orbital periods. The dashed black line shows the type I migration speed from Tanaka et al. (2002, at r = r0 constant). Bottom panel: migration timescale τ = a/|ȧ| normalised by the angular velocity of the planet. The black dashed and dotted lines represent the critical migration timescale, calculated with Eq. (3), beyond which the planet would cross the 2:1 or 3:2 resonance, respectively. These simulations were run for viscosities corresponding to α = 10−4 and 10−3. The migration timescales plotted in this figure were smoothed by means of a sliding window average with a size of then and 50 points (which represents a time interval of about 100 and 600 yr) for the case α = 10−3 and α = 10−4, respectively.

In the text

Fig. 13 Relative migration timescale $\tau =\frac{a_{\text{S}}/a_{\text{J}}}{d\left(a_{\text{S}}/a_{\text{J}}\right)/dt}$ in a simulation where Saturn migrates in the cold disc perturbed by Jupiter, without feeling Jupiter’s gravity. The normalised timescale is plotted as a function of aS/ aJ, and the vertical black dashed line marks the location of the 2:1 MMR, while the grey dotted vertical line marks the position of the 5:2 MMR. The horizontal dashed line marks the critical migration timescale, calculated from Eq. (3) (for Mp = Ms), above which the planets should not cross the 2:1 resonance. The migration timescale was smoothed by means of a sliding window average with a size of 100 points (which represents a time interval of about 12 000 yr), in order to suitably reduce the noise.

In the text

	Fig. 14 Same as Fig. 5, but for the simulations $N_{M_{\text{J}}},N_{2M_{\text{J}}}$, and $N_{3M_{\text{J}}}$. We note that $N_{2M_{\text{J}}}$ becomes unstable shortly after 130 000 yr.
In the text

	Fig. A.1 Migration of Jupiter in a disc with viscosity α = 10−4 for higher resolutions.
In the text

	Fig. A.2 Migration of Jupiter in a disc with viscosity α = 10−5 for higher resolutions.
In the text

Fig. B.1 Orbital evolution of Jupiter and Saturn for simulation set C. The top and bottom panels show the semi-major axis and the eccentricity, respectively. The vertical dotted lines indicate the time t0 at which we started the simulation with gas removal, namely ’C + Gas removal’. The nominal simulation was continued for t > t0 until a significant difference appeared with respect to the gas removal.

In the text

	Fig. B.2 Same as Fig. B.1, but for simulation set C3.
In the text

	Fig. B.3 Resonant angles corresponding to the simulations shown in Fig. B.1 in the top panel, and those corresponding to Fig. B.2 in the bottom panel. We find that in both cases of 2:1 or 5:2 MMR, the planets remain locked into resonance after the gas removal. Although, their resonant angles either librate with larger amplitude or no longer around 0.
In the text

	Fig. C.1 Migration of Jupiter and Saturn if Saturn was to form on the inside of the of 2:1 MMR. The top panel shows the evolution of the semi-major axes of Jupiter, in red, and Saturn, in blue. The dashed black and grey lines represent the position of the 2:1 and 3:2 MMR with Jupiter, respectively. The bottom panel shows the eccentricities of both planets.
In the text