Migration of pairs of giant planets in low-viscosity discs

When considering the migration of Jupiter and Saturn, a classical result is to find the planets migrating outwards and locked in the 3:2 mean motion resonance (MMR). These results were obtained in the framework of viscously accreting discs, in which the observed stellar accretion rates constrained the viscosity values. However, it has recently been shown observationally and theoretically that discs are probably less viscous than previously thought. Therefore, in this paper, we explore the dynamics of pairs of giant planets in low-viscosity discs. We performed two-dimensional hydrodynamical simulations using the grid-based code FARGOCA. In contrast to classical viscous discs, we find that the outer planet never crosses the 2:1 resonance and the pair does not migrate outwards. After a wide parameter exploration, including the mass of the outer planet, we find that the planets are primarily locked in the 2:1 MMR and in some cases in the 5:2 MMR. We explain semi-analytically why it is not possible for the outer planet to cross the 2:1 MMR in a low-viscosity disc. We find that pairs of giant planets migrate inwards in low-viscosity discs. Although, in some cases, having a pair of giant planets can slow down the migration speed with respect to a single planet. Such pairs of slowly migrating planets may be located, at the end of the disc phase, in the population of exoplanets of 'warm Jupiters'. However, the planets never migrate outwards. These results could have strong implications on the Solar System's formation scenarios if the Sun's protoplanetary disc had a low viscosity.


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 1 .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. 1 From the Open Exoplanet Catalogue by Rein (2012).
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 thermodynamical 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) 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 magnetohydrodynamics (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 −4 by modelling dust settling.With the same method, Villenave et al. (2022) found a value of the order of 10 −5 .Another example of turbulent parameter <10 −3 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 lowviscosity 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.

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 0 (r/r 0 ) 2/7 with h 0 = 0.05, unless otherwise stated, and r 0 is the unit length (see Sect. 2.2).The disc density is given by Σ = Σ 0 (r/r 0 ) −1/2 where Σ 0 = 6.76 × 10 −4 in code units for the nominal case.This corresponds to 222 g cm −2 for r 0 = 5.2 au and a Sun mass star.We used the α prescription (Shakura & Sunyaev 1973) for the kinematic viscosity, such that ν = αc s H, 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 −4 .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 −5 (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 −3 .The outer planet, unless otherwise stated, is a Saturnmass planet with mass M S /M = q S = 2.9 × 10 −4 .For simplicity, the two planets are henceforth referred to as Jupiter and Saturn, respectively.

Units
The code units are such that G = M = r 0 = 1, where r 0 is the unit length and the initial semi-major axis of Jupiter.Therefore, the orbital period at r 0 = 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: Indeed, in the three cases, M = 1, M p = 10 −3 , and r 0 = 1 are in code units and therefore Σ 0 = 3 × 10 −3 M r −2 0 .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 2 .In other words, only the dimensionless quantities q, α, h 0 and the disc-to-star mass ratio µ = πΣ 0 r 0 2 /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 distance 3 .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 0 = 5.2 au, so that one orbital period at r 0 corresponds to 11.85 yr.

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)  4 .The code was parallelised using a hybrid combination of MPI between the nodes and OpenMP on shared memory multi-core processors.
2 More detailed calculations: 3 Only for a system governed by the equations that we simulated and in which no other physical phenomena have significant effects (e.g.UV photo-evaporation by the central star). 4The simulations presented in this paper have been obtained with a recently re-factorised version of the code that can be found at https: //disc.pages.oca.eu/fargOCA/public/The radial domain extends in the range r ∈ [0.2, 9]r 0 , which corresponds to r ∈ [1.04, 46.8] au for r 0 = 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: 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).

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 0 from the star, with the remaining orbital elements being initialised at zero.Its mass was increased analytically 5 from zero up to q J in 800 local orbits, following where T growth = 2π × 800 a J r 0 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 a S a J T 0,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, a S a J T 0,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.

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.

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 θ J,S = 2λ S − λ J − 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.

Disc thickness
In classical discs, for α ≥ 10 −3 , 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 = 0.035 and with a hotter, thicker disc, H with h 0 = 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 2 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.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.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).

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.

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 difference6 .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   7, but for simulation C 3 in the time interval [75 000, 80 000] yr.It is important to note the jump of Saturn between masses 68% (r = 2.12a J ) and 71% (r = 1.76aJ ) in 100 orbits at r 0 .
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 3 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 80 000 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.

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 −3 (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.

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 −4 and 10 −3 as well as the aspect ratio of h 0 = 0.035 and we studied its migration as a function of its mass.The top panel of Fig. 12  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 where τ = a |ȧ| is normalised by the angular velocity of the outer planet, and f (1)  res 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 −3 , 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 Fig. 13.Relative migration timescale τ = a S /a J d(a S /a J )/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 a S /a J , 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 M p = M S ), 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.
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 Jupiter7 .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 semimajor 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 = 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.

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 Fig. 14.Same as Fig. 5, but for the simulations N M J , N 2M J , and N 3M J .We note that N 2M J becomes unstable shortly after 130 000 yr. 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 = 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.

Discussion and conclusions
In this paper, we have explored the dynamics of pairs of giant planets in low-viscosity discs.We ran two-dimensional hydrodynamical 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 6 yr.While in the thin disc case, C, once in resonance, Jupiter migrates at a speed of 2 au/10 6 yr; in comparison, the migration speed of Jupiter alone in the cold disc is about 10 au/10 6 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 −4 .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. 2018Clement et al. , 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 −4 , 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.
A190, page 1 of 13 Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0),which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.This article is published in open access under the Subscribe to Open model.Subscribe to A&A to support open access publication.
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.

Fig. 2 .
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, T 0,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 .
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 periand 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 .
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 .
Fig. 5. Orbital evolution of Jupiter and Saturn for the simulations C, N, and H (see Table1).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 .
Fig. 6.Evolution of the four resonant angles in simulation C.

Fig. 7 .
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.

Fig. 10 .
Fig.10.Same as Fig.5, but for simulations C, C 2.5 , and C 3 .Simulation C 2.5 was stopped at about 110 000 yr as it overlaps almost perfectly with C in both orbital parameters.
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, M S , in 800 initial orbital periods.The dashed black line shows the type I migration speed from Tanaka et al. (2002, at r = r 0 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.

Table 1 .
Summary of the simulation parameters explored in this study.
Notes.Column 1 -simulation name.Column 2 -parameters: mass and starting position of the outer planet with respect to the position of the inner planet, disc aspect ratio, and disc surface density.Column 3mean motion resonance (MMR) of the planets at the end of the simulation.Bold parameters indicate the changed values with respect to the nominal case N. ( * ) Unstable configuration after 130 kyr.