Issue 
A&A
Volume 646, February 2021



Article Number  A166  
Number of page(s)  18  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202039520  
Published online  26 February 2021 
Migration of Jupitermass planets in lowviscosity discs
^{1}
Laboratoire Lagrange, UMR7293, Université de Nice SophiaAntipolis, CNRS, Observatoire de la Côte d’Azur, Boulevard de l’Observatoire,
06304
Nice Cedex 4,
France
email: elena.lega@oca.eu
^{2}
Astronomy Unit, School of Physics and Astronomy, Queen Mary University of London,
London
E1 4NS,
UK
^{3}
Institut für Astronomie und Astrophysik, Universität Tübingen, Auf der Morgenstelle 10,
72076
Tübingen,
Germany
Received:
24
September
2020
Accepted:
21
December
2020
Context. TypeII migration of giant planets has a speed proportional to the disc’s viscosity for values of the α viscosity parameter larger than 10^{−4}. Previous studies based on twodimensional simulations, have shown that, at even lower viscosities, migration can be very chaotic and is often characterised by phases of fast migration. The reason is that vortices appear in lowviscosity discs due to the Rossbywave instability at the edges of the gap opened by the planet. Migration is then determined by vortexplanet interactions.
Aims. Our goal is to study giant planet migration in lowviscosity discs with 3D simulations. In 3D, vortices are more complex than the simple vertical extension of their 2D counterparts; their impact on planet migration is therefore not obvious.
Methods. We performed numerical simulations using two gridbased codes: FARGOCA for threedimensional simulations and FARGOADSG for the two dimensional case. Twodimensional simulations were used mainly for preliminary tests to check the impact of selfgravity on vortex formation and on vortexdisc dynamics. After selecting disc masses for which selfgravity is not important at the planet location, threedimensional simulations without selfgravity can be safely used. We have considered an adiabatic equation of state with exponential damping of temperature perturbations in order to avoid the development of the vertical shear instability. In our nominal simulation, we set α = 0 so that only numerical viscosity is present. We then performed simulations with nonzero α values to assess the threshold of prescribed viscosity below which the new migration processes appear.
Results. We show that for α ≲ 10^{−5} two migration modes are possible, which differ from classical TypeII migration in the sense that they are not proportional to the disc’s viscosity. The first occurs when the gap opened by the planet is not very deep. This occurs in 3D simulations and/or when a big vortex forms at the outer edge of the planetary gap, diffusing material into the gap. The desaturation of coorbital and corotation resonances keeps the planet’s eccentricity low. Inward planet migration then occurs as long as the disc can refill the gap left behind by the migrating planet, either due to diffusion caused by the presence of the vortex or to the inward migration of the vortex itself due to its interaction with the disc. We call this type of migration ‘vortexdriven migration’, which differs from ‘vortexinduced’ migration described in Lin & Papaloizou (2010, MNRAS, 405, 1473, and 2011a, MNRAS, 415, 1445). This migration is very slow and cannot continue indefinitely because eventually the vortex dissolves. The second migration mode occurs when the gap is deep so that the planet’s eccentricity grows to a value e ~ 0.2 due to inefficient eccentricity damping by corotation resonances. Once the planet is on an eccentric orbit, gas can pass through the gap and planet migration unlocks from the disc’s viscous evolution. This second, faster migration mode appears to be typical of twodimensional models in discs with slower damping of temperature perturbations.
Conclusions. Vortexdriven migration in lowviscosity discs can be very slow and eventually reverses and stops, offering an interesting mechanism to explain the existence of the coldJupiter population, even if these planets originally started growing at the disc’s snowline.
Key words: methods: numerical / planets and satellites: dynamical evolution and stability / protoplanetary disks / planetdisk interactions
© E. Lega et al. 2021
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.
1 Introduction
The origins of giant planets remain elusive. Radial velocity surveys have found giant planets to exist around roughly 10% of Sunlike stars (Mayor et al. 2011; Cumming et al. 2008). However, only ~1% of Sunlike stars have hot Jupiters on very shortperiod orbits (Howard et al. 2010). Very few stars have warm Jupiters with orbital radii as large as 0.5–1 au (Butler et al. 2006; Udry et al. 2007). Instead, when considering the unbiased distribution, most giant planets are found between 1 and several au (Butler et al. 2006; Udry et al. 2007; Cumming et al. 2008; Howard et al. 2010; Mayor et al. 2011), while there are hints that their number decreases again farther out (Mayor et al. 2011; Fernandes et al. 2019). In our Solar System, of course, there are no giant planets within 5 au from the Sun, although Jupiter may have been at ~2 au in the past (Walsh et al. 2011). The preference for giant planets to orbit relatively far from the parent star, in contrast to superEarths for example, is puzzling because giant planets must have formed in the presence of gas in the protoplanetary disc and, consequently, they should have migrated towards the central star due to planetdisc interactions.
Giant planet migration, also known as TypeII migration, has been extensively studied in the literature since the pioneering work by Lin & Papaloizou (1986). In essence, a giant planet opens a gap in the gas around its orbit. Any migration of the planet has to occur in concert with the readjustment of the gas in the disc for the gap to migrate together with the planet. In a viscous accretion disc, this can happen only on a viscous timescale. Hence TypeII migration is expected to be inwards, with a radial speed of the order of v_{r} ~ 1.5ν∕r, where ν is the viscosity, that is to say, with the same radial speed at which the gas viscously accretes towards the central star (Ward 1997). Young stars typically accrete gas at a rateof ~10^{−8}M_{⊙} y^{−1} – with a large, order of magnitude scatter around this value (Hartmann et al. 1998; Manara et al. 2016) – for a density of gas in the disc at 1 au comparable to that of the Minimum Mass Solar Nebula model (2 × 10^{−4}M_{⊙} au^{−2}; Weidenschilling 1977; Hayashi 1981); this implies a radial speed of ~ 8 au My^{−1}. This means that, at TypeII migration speeds in viscous discs, giant planets should move towards their host stars on timescales much shorter than disc lifetimes. If this were true, the presence of giant planets at several au from the parent star would require these planets to have formed very far away (i.e. beyond ~ 20 au) and/or quite late in the history of the protoplanetary disc (Coleman & Nelson 2014; Bitsch et al. 2015).
However, this simple solution is not without problems. The sweet spot for the rapid formation of massive cores capable of accreting gas and becoming giant planets is the snowline (Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017). If the giant planets we observe today orbiting beyond 1 au formed at distances ≥ 20 au, where are the planets that should have formed faster and in larger numbers near the snowline? A popular idea was that these planets have fallen onto the central star and the planets that we observe now are the last of the Mohicans (Lin 1997; Laughlin & Adams 1997). Today this idea has been mostly discarded because of the pileup of extrasolar planets on orbits with periods of 3–10 days, and modern models of interactions between stellar magnetospheres and protoplanetary discs (Matt & Pudritz 2005; Mohanty & Shu 2008; Adams 2012) suggest that discs are truncated at a few 10^{−2} au, which should prevent planet migration continuing all the way to the stellar surface.
If the giant planets we observe formed at the snowline, they must have migrated only a few astronomical units during the lifetime of the disc. Because the rate of TypeII migration is in principle proportional to the disc viscosity, this may suggest very low viscosities in protoplanetary discs. This idea is supported by modern studies on turbulent viscosity in discs. Turbulence was originally expected to arise from the magnetorotational instability (MRI; Balbus & Hawley 1991). However, it was later understood that the ionisation of the gas near the midplane of the disc is too weak to sustain the MRI (Gammie 1996; Stone et al. 1996), introducing the concept of the dead zone. Even more recently, the inclusion of nonideal MHD effects, such as ambipolar diffusion, led to the conclusion that the coupling between the magnetic field and the gas should not make the disc turbulent even at its surface (see Turner et al. (2014) for a review). Another often considered source of turbulence, the vertical shear instability (VSI hereafter; Nelson et al. 2013; Stoll & Kley 2014) should also not be active at the disc midplane within a few astronomical units from the star because the disc’s cooling rates are too slow (however, see Pfeil & Klahr 2020 for a different view). Thus, the idea of formation and migration of giant planets in low viscosity discs may be appealing.
Unfortunately, the dependence of TypeII migration on viscosity in the limit of small viscosity is far from clear. Duffell et al. (2014) claimed that the paradigm of TypeII migration is flawed because the gas can pass through the planet’s orbit, from one side of the gap to the other, so that the planet is not locked in the viscous evolution of the disc. For this reason, the planet’s migration speed can be different from − 1.5ν∕r where ν is the viscosity. Dürmann & Kley (2015) confirmed this result, but Robert et al. (2018) showed that the passage of gas through the gap is inhibited at small viscosity because the gap is much larger than the planet’s horseshoe region; consequently the planet’s migration rate remains proportional to ν although not necessarily equal to − 1.5ν∕r. However, all these studies have been performed using twodimensional simulations and, more importantly, have considered viscosities that were not very small for numerical reasons. Adopting the usual prescription ν = αh^{2}Ω (Shakura & Sunyaev 1973), where h is the pressure scaleheight of the disc and Ω is the orbital frequency, the viscosities considered in these works corresponded to α > 10^{−4}. For these values of α, TypeII migration, even if proportional to ν, is still too fast to explain the current location of most giant planets, if these planets formed at the snowline. At face value, α should be ~ 10^{−5} for TypeII migration to be slow enough, if the proportionality between migration speed and α (or ν) holds also for such a low viscosity.
Attempts to simulate planet migration in inviscid discs have led to very chaotic and erratic planetary evolutions, often characterized by phases of very fast migration (Lin & Papaloizou 2010, 2011a; McNally et al. 2019). The reason is that in low viscosity discs with embedded giant planets large scale vortices can appear due to edge instabilities of planetary gaps and Rossby wave instabilities at pressure bumps (Lovelace et al. 1999; Koller et al. 2003; Lin & Papaloizou 2010). The interaction of the planet with a vortex can lead to several effects. For instance, the passage of the vortex from one side of the gap to the other can accelerate the planet’s migration in the opposite direction (called ‘vortexinduced’ migration in Lin & Papaloizou 2010). A further complication is that gapopening planets in inviscid discs may undergo eccentricity growth because of (i) their interaction with the vortices and (ii) the saturation of corotation resonances that are normally responsible for damping eccentricity (Goldreich & Tremaine 1980; Goldreich & Sari 2003); in turn the planet eccentricity can have a strong feedback on planet migration.
However, it has been shown that the formation and evolution of vortices depends sensitively on the simulation setup. In 2D simulations, Zhu & Baruteau (2016) showed that vortices generated by the Rossby wave instability are less pronounced if the selfgravity of the disc is taken into account. Using twodimensional (2D hereafter) models Lin & Papaloizou (2011a,b) showed that, when considering selfgravity in discs with suitable mass, the vortexinduced migration observed in Lin & Papaloizou (2010) is considerably delayed. In threedimensional (3D hereafter) models, vortices are more complex than the vertical extension of their 2D counterparts (Barranco & Marcus 2005; Meheut et al. 2010). Thus the interplay between planets and vortices and the effects of this interplay on migration are far from clear.
In this context, the goal of this paper is to provide a comprehensive exploration of giant planet migration in lowviscosity discs. We approach this problem with a suite of 2D simulations with and without disc selfgravity and 3D simulations without selfgravity. More specifically, after a brief presentation of our physical models (Sect. 2) and of the simulation setups (Sect. 3), we present in Sect. 4 the structure of an inviscid disc under the effect of a Jupitermass planet kept on a fixed circular orbit at 5.2 au, with emphasis on the generation of vortices and their evolution. Then, in Sect. 5 we investigate the planet migration in inviscid discs using both 2D and 3D simulations. We also provide the range of validity of our results by finding the transition to classical type II migration. We show in Sect. 6 that results depend on the planet formation site. Numerical tests such as: convergence with respect to resolution and gravitational smoothing of the potential in the planet vicinity as well as a discussion of dependence of results on the choice of the equation of state are discussed in the appendix. The conclusions and a global discussion of the results are reported in Sect. 7.
In a subsequent paper, we will continue investigating giant planet migration in lowviscosity 3D discs in two directions: (i) taking into account selfgravity in 3D discs for the study of migration of giant planets forming at large distances from the star (~ 20 au) (ii) modeling the effect of the radial advection of gas in the disc due to angular momentum removal by magnetised disc winds (Bai & Stone 2013; Turner et al. 2014; Gressel et al. 2015; Béthune et al. 2017). Advection is required in a realistic model of low viscosity discs to provide a mass flux to the star comparable to the observed stellar accretion rates. The effect of advection on giant planets migration will be the object of a forthcoming paper. Nevertheless, a deep understanding of planet migration in lowviscosity discs with no advection, as that developed in this paper, is required in order to understand giant planet migration in a more realistic scenario.
2 Physical models
We briefly describe our 2D and 3D models in this section.
2.1 The 3D model
The protoplanetary disc is treated as a non selfgravitating gas whose motion is described by the NavierStokes equations. We use spherical coordinates (r, φ, θ) where r is the radial distance from the star (which is at the origin of the coordinate system), φ is the azimuthal coordinate measured from the xaxis and θ the polar angle measured from the zaxis (the colatitude). The midplane of the disc is located at the equator . We work in a coordinate system which rotates with angular velocity:
where M_{⋆} is the mass of the central star, G is the gravitational constant, and r_{p}(0) is the initial distance to the star from a planet of mass m_{p}, assumed to be on a circular orbit. The gravitational influence of the planet on the disk is modelled as in Kley et al. (2009) using the full gravitational potential for disc elements having distance d from the planet larger than a fraction ϵ of the Hill radius, and a smoothed potential for disc elements with d < ϵ.
We integrate the NavierStokes equations taking into account indirect forces that account for the acceleration of the star by the disc and planet (Masset 2002). As shown in Zhu & Baruteau (2016), indirect forces have an impact on shaping vortices that is proportional to the disc mass. We add an equation for the internal energy e = ρc_{v}T to the NavierStokes equations, where ρ and T are the volume density and the temperature of the disc gas and c_{v}v is the specific heat at constant volume: (1)
where τ_{c} is the cooling time and T_{0} is the initial temperature, defined as , with h_{0} being the disc aspect ratio, μ the mean molecular weight (μ = 2.3 g mol^{−1} for a standard solar mixture) and R_{gas} is the idealgas constant.
In short, we use an adiabatic EoS, on top of which we exert an exponential damping of the temperature perturbations. We do not use the simpler locally isothermal EoS in order to avoid disc instabilities like the VSI. In the quoted paper it is shown that a disc with anequation of state taking into account thermal relaxation like in Eq. (1) with suitable values of τ_{c} is not prone to develop such an instability. According to the same paper we will consider τ_{c} equal to 1 orbital period at the planet location. A selfconsistent simulation of the disc’s evolution indeed shows that no significant VSI should develop at Jupiter’s distance (Ziampras et al., in prep., but see Pfeil & Klahr 2020 for results with a different simulation setup.). The results for different cooling times are reported in Appendix A.2.
2.2 The 2D model
Twodimensional models have been commonly used in the literature for studies of planetary migration. Evolution over long timescales (up to hundreds of thousands of planetary orbits) is possible at reasonable computational cost using 2D simulations, which makes them useful for exploring long term behaviour. The obvious disadvantage is that possible genuine 3D effects cannot be observed.
In the twodimensional model we solve the verticallyintegrated NavierStokes equations using polar coordinates (r, φ) and, similar to the 3D model described above, we also solve Eq. (1) with thermal relaxation, where in 2D the internal energy density is defined by e = Σc_{v}T, and Σ is the surface density. We also consider a cooling time equal to 1 orbital period at the planet position, and we will discuss in Appendix A.1 the sensitivity of results to the choice of the cooling time.
We undertake 2D simulations that either include or neglect the effects of the disc selfgravity. Vortices are an inevitable consequenceof gap opening by planets in inviscid discs, and their evolution is expected to depend on the relative influences of the indirect term and selfgravity (Zhu & Baruteau 2016). Hence, it is important to test the effects of selfgravity and to determine under whichconditions it is important or can be neglected.
The 2D simulations are performed in a frame of reference in which the star is located at the origin, and which rotates at a rate that corresponds to the instantaneous angular velocity of the planet (i.e. corotating with the planet even when the planet migrates).
3 Setup of numerical simulations
Our 3D simulations are done with the code FARGOCA (FARGO with Colatitude Added; Lega et al. 2014)^{1}. The code is based on the FARGO code (Masset 2000) extended to three dimensions. The fluid equations are solved using a second order upwind scheme with a timeexplicitimplicit multistep procedure. The code is parallelised using a hybrid combination of MPI between the nodes and OpenMP on shared memory multicore processors. The two dimensional simulations are done with FARGOADSG, a version of the FARGO code which implements an energy equation and disc selfgravity as explained in Baruteau and Masset (2008).
The code units are G = M_{*} = 1, and the unit of distance r_{1} = 1 is arbitrary when expressed in au. The unit of time is therefore . For the simulations in this paper we adopt the SunJupiter distance as the unit of length in au: r_{1} = 5.2. When presenting simulation results distances are expressed in au and time in years.
3.1 Disc parameters
Unless otherwise stated, all models have a radial domain defined by r_{min} ≤ r ≤ r_{max} with r_{min} = 1.04 au and r_{max} = 46.8 au. Such a wide radial extension is quite unusual in these kind of studies. The paper focuses on the case of a planet initially at 5.2 au but, in order to test the dependence of the results on the planetstar distance, we also considered the case of a planet initially at 13 au (see Sect. 6), and in this case we found that an extended radial domain is required to ensure that results are not affected by the outer boundary. In the radial direction we use the classical prescription (de ValBorro et al. 2006) of evanescent boundary condition.
Importantly, in our nominal simulations, we set α = 0, namely we consider inviscid discs. Of course, there will still be some resolutiondependent numerical diffusivity in the code. This issue will be discussed in Sect. 5.4, but we anticipate here that, for our nominal resolution, the numerical viscosity will not be larger than the one in a viscous disc with α = 10^{−5}. Moreover, in all the simulations (2D and 3D) we consider artificial viscosity to stabilize shocks (precisely we use von NeumannRichtmyer viscosity and corresponding heating terms as described by Stone & Normann 1992).
In the 3D models, the meridional domain extends from the midplane to 12° above the midplane (θ = 78°–90°), about four disc scaleheights. We do not study inclined planets and therefore we do not need to extend the domain below the midplane. Mirror boundary conditions are applied at the midplane as in Kley et al. (2009) and reflecting boundaries are applied at the disc surface. The initial temperature profile is the same for all the 3D discs: isothermal in the vertical direction, while in the radial direction the temperature scales as 1∕r. This gives rise to a model in which the disc aspect ratio, h, is constant with radius, and we adopt the value h = h_{0} = 0.05.
In the 2D models we adopt a radial profile for the disc aspect ratio given by Chiang & Goldreich (1997) for a passively heated disc: with h_{0} = 0.05. This choice was made because in earlier test calculations undertaken with h = 0.05 throughout the disc, the value of the Toomre Q parameter approaches unity at the edge of the disc models, such that the effects of selfgravity on global disc evolution became very obvious in the form of global spiral waves being excited. Adopting a variable h removed this problem.
We remark that the aspect ratio is the same in 2D and 3D models at the planet’s location, 5.2 au, and is slowly diverging away from this radius. The aspect ratio does not change during the simulations^{2}. This is true also for simulations with non zero viscosity (see Sect. 5.4). Precisely, we do not include the term corresponding to the viscous heating in Eq. (1) in order to have the same vertical structure for all the simulations and make them comparable.
The disc surface density in 2D and 3D models is given by , with Σ_{0} = 6.76 × 10^{−4} in code units for our nominal disc (corresponding to 222 g cm^{−2} at 5.2 au). We will call the corresponding simulations ’nominal’ or M simulations (see Table 1). Classic TypeII migration is not dependent on the disc mass, except in the inertial limit where the local disc mass is much smaller than the planet mass and migration slows down (Quillen et al. 2004). However, in a low viscosity disc, in which vortices can form at the edge of the planetinduced gap, the mass of the disc may play an important role. In particular, Zhu & Baruteau (2016) showed that the strength and evolution of vortices depends on both the gravitational acceleration that they exert on the central star (as expressed through the indirect term) and the selfgravity of the disc. Therefore, we also run simulations with the nominal disc mass divided by some integer n, we will call these simulations M∕n or Nominal/n.
Achieving long run times for 3D simulations is challenging in terms of computational cost, therefore we have used a moderate resolution of (N_{r}, N_{θ}, N_{φ}) = (568, 16, 360) with uniform radial grid spacing. At this resolution we consider a smoothing length ϵ = 0.8R_{H}, that is four gridcells in the smoothing length. Although this value is quite large, it has been shown in Fung et al. (2017) that convergence in the torque measurements in non viscous disks requires at least three grid cells per smoothing length. We provide a test in the Appendix with smoothing length ϵ = 0.4 R_{H} at the same resolution and show that results are converged, and we further reduce the smoothing length (ϵ = 0.2 R_{H}) in the case where we double the resolution (see Appendix B). We note that in our comparison between 2D selfgravitating and nonselfgravitating discs, we adopt the same number of cells in the radial and azimuthal directions as in the 3D runs, but the radial grid spacing is logarithmic, as is required by the selfgravity solver in FARGOADSG. Table 1 lists the principal parameters and simulation names to which we will refer in the following. We notice that viscosity is labeled with the α parameter. Since the vertical structure remains constant in the simulations α gives directly the kinematic viscosity ν through the relation ν = αh^{2}Ω.
Simulation parameters for 3D and 2D discs.
3.2 Planet growth on fixed orbits
In the 2D simulations the planet grows to its final mass over 800 orbits while being held on a fixed circular orbit. In the migration runs, the planet is released after the system has been evolved for a further 400 orbits. In the 3D runs theplanet grows to its final mass over about 200 orbits, and the system is evolved for a further 600 orbits before the planet is released in the migration runs. In this way we avoid the excitation of instabilities that would arise if the planet were initialised with its final mass (see also Hammer, Kratter, & Lin 2017; Hallam & Paardekooper 2020).
Fig. 1
Toomre Q values for 2D selfgravitating disc models. We focus on the nominal or M disc, but we also run control simulations with the mass of the disc divided by 2, 4 and 10. The Q = 1 value is indicated by the dotdashed line, and the locations of the planets at r = 5.2 and r = 13 au are indicated by the filled red circles. 

Open with DEXTER 
4 Planets on fixed orbits: comparing disc structures for simulations in 2D (with and without selfgravity) and in 3D
As discussed in Sect. 1, Zhu & Baruteau (2016) have shown that the evolution of a protoplanetary disc with a vortex at a pressure bump depends on whether the indirect term and selfgravity are accounted for. Neglecting both of them is appropriate only for a disc with very low mass and results in a single vortex that remains at the location of the pressure bump. For more massive discs, the inclusion of the indirect term, without selfgravity, results in a large, radiallyextended vortex that migrates inwards because of the excitation of strong spiral density waves. This leads to significant global restructuring of the disc. Including also selfgravity results in a more stable evolution, with one or more vortices forming at the pressure bump and remaining there for the duration of the simulation.
The formation of a gapforming planet in an inviscid protoplanetary disc inevitably leads to the formation of a vortex at the pressure bump located at the outer edge of the gap, and the subsequent evolution is expected to depend on the disc mass because of the influence of the indirect term and selfgravity. In this section, we first present results from a suite of 2D simulations that specifically address the question: under which conditions does selfgravity become unimportant? Then, for the conditions for which selfgravity is not important we present 3D simulations without selfgravity. We note that all simulations include the indirect term.
Zhu & Baruteau (2016) suggested that a controlling factor is the value of the Toomre Q parameter at the vortex location. The radial profile of Q is shown in Fig. 1 for different disc models: M, M∕2, M∕4 and M∕10, where the disc mass is nominal or divided by 2, 4 and 10 respectively. Because Q depends on radius, we consider models in which the giant planet is located at either r_{p} = 5.2 or r_{p} = 13 au. We discuss in the following the case of a planet at r_{p} = 5.2 and we dedicate Sect. 6 to a discussion about planets forming farther out from the star (r_{p} = 13 au). Our main suite of 2D and 3D simulations adopt a local cooling time of τ_{c} = 1 orbital period. We have found, however, that qualitative changes in the results (especially in 2D) arise when the cooling time changes, and hence we present a brief analysis of how the cooling time affects the results in Appendix A.
4.1 Planet at r_{p} = 5.2 au, 2D model
We consider the disc evolution with the planet orbiting at r_{p} = 5.2 au using the nominal parameters described in Sect. 3.1. We note the local Toomre parameter Q(r = 5.2 au) = 23.5, the aspect ratio h(r = 5.2 au) = 0.05 and the local cooling time τ_{c} = 1 orbital period.
Our discussion will focus on the part of the disc that orbits exterior to the planet. During the early stages, the simulations show the formation of a gap and of three or four vortices that quickly merge into a single vortex at the outer edge of the gap.
Results for the nominal M_{2D}sg and M_{2D} simulations are shown in Fig. 2. The selfgravitating simulation produces a vortex at the gap outer edge, which digs a shallow secondary gap at about 10 au by the emission of spiral waves, and weakens over time. Figure 3 shows that the size of the pressure bump slowly decreases over time. The nonselfgravitating simulation produces a stronger vortex and a deeper secondary gap than in M_{2D}sg, but the vortex is also observed to weaken over time and in the rightmost panels in Fig. 2 we see that the final discs have very similar structures. Comparing selfgravitating and nonselfgravitating models, the radial profiles of the azimuthally averaged surface densities are seen to differ in detail in Fig. 3, mostly concerning the depth of the secondary gap at ~ 10 au, but are qualitatively similar. In both the selfgravitating and the nonselfgravitating simulations, the outer edge of the main gap is found to become Rayleigh unstable after approximately 800 planet orbits (~ 9500 yr) and this coincides with the dissipation of the vortex, suggesting a causal link.
We have checked (Fig. 3, middle and right panels) that, when we reduce the disc mass by factors of 2 and 4, the selfgravitating and nonselfgravitating models show rapid convergence in their behaviour even concerning small details, such as the depth of the secondary gap. We conclude that for the chosen parameters, these models all show good convergence in their behaviour regardless of whether selfgravity is included or not.
4.2 Planet at r_{p} = 5.2 au, 3D model
The 2D approximation is commonly used since protoplanetary discs are vertically thin and very often 2D results are shown to be consistent with results obtained with 3D models. If this is true in viscous discs, caution has to be taken for inviscid discs. Results for the M_{3D} simulation are presented in Fig. 2, bottom panels. We notice that a vortex forms at the outer edge of the gap also in 3D, and that it is radially larger and more massive than the vortex formed in the 2D models. Vortex exchange of angular momentum with the disc results in the evolution of the gap’s outer edge. On the timescale of Fig. 2, the vortex appears to be persistent over time, though slightly smoothing out. The vortex exerts a pressure wave on the disc which superposes and interferes (sometimes non linearly) with the wake launched by the planet, which explains why the wave pattern observed in Fig. 2 (bottom panels) is more complex than that observed in 2D simulations. Through the wave the vortex exerts a torque on the disc which is effective enough to open a secondary gap (Fig. 2, bottom right panel) that is more pronounced than in the equivalent 2D simulation.
When looking at the radial surface density profiles in Fig. 4, we see that the planetinduced gaps appear very different. In 2D the gap becomes deeper and wider with time. In 3D the gap looks almost stationary over time; the gap is even partially refilled and, at time t = 23 720 yr, contains three orders of magnitude more gas than in the 2D simulations. This difference between 2D and 3D gaps in low viscosity discs was already shown in Morbidelli et al. (2014). In 3D the gap’s depth is regulated by the gas meridional circulation which continuously refills the gap, so that the gap depth saturates at a much larger value than in 2D simulations. Consequently the disc does not satisfy Rayleigh’s criterion for instability and the vortex is not smoothed out by the diffusion that this instability generates. In turn, the vortex exerts a negative torque on the outer edge of the gap, which contributes significantly to limiting the gap’s depth. The difference between 2D and 3D simulations persists also for smaller disk masses. However, Appendix A will show that the results of 2D simulations will approach those of 3D simulations when shorter cooling times are assumed, e.g. τ_{c} = 0.01 orbit.
Fig. 2
Top panels: contours of surface density for the M_{2D}sg simulation of the nominal disc, for times shown (in years) at the top of each panel. Middle panels: same as top panels but for the M_{2D} simulation (i.e. without selfgravity). Bottom panels: same as top panels but for the M_{3D} simulation. The planet is located at r = 5.2 au. 

Open with DEXTER 
5 Planet migration
We investigate the migration of the planet initially at r = 5.2 au, for which we have shown that selfgravity has essentially no impact even for the nominal disc mass and τ_{c} = 1 orbit. We have shown that 2D and 3D discs behave very differently for these model parameters; it is therefore instructive to study planet migration in both cases in order to assess the range of behaviours that can occur in inviscid disc models even though, in principle, the 3D case should be more realistic than the 2D one.
In order to better understand what follows, remember that a giant planet opening a gap in the disc can migrate only at the speed at which the disc can readjust itself in order for planet and gap to migrate together (Ward 1997). In viscous discs, the gas can displace radially on the viscous timescale, which sets the planet migration speed to be proportional to ν∕r. In the inviscid case, the gap can move radially only in two ways. One possibility is that the planet displaces gas from the inner edge to the outer edge of the gap (Duffell et al. 2014; Dürmann & Kley 2015). Robert et al. (2018) showed that this is not possible in the limit of low viscosity if the planet is on a circular orbit, because the gap is too wide with respect to the location of the separatrices of the planet’s horseshoeshaped coorbital region. But, if the planet acquires a sufficiently eccentric orbit, it can transfer gas from the inner to the outer disc. Then it can migrate inwards, irrespective of the viscous timescale. The other possibility is that the vortex forming at the outer edge of the gap pushes gas into the gap or migrates itself inwards due to the pressure waves that it generates in the disc (Paardekooper et al. 2010; Zhu & Baruteau 2016). In this case the outer disc can spread inwards on a timescale that is not related to viscous relaxation but is that of vortexdisc interactions (unrelated to viscosity). Thus, the planet can migrate on the same timescale, while pushing the inner disc inward. As we show below, 2D and 3D simulations display these two mechanisms.
Fig. 3
Azimuthallyaveraged radial surface density profiles versus time for 2D runs with and without selfgravity, and a Jupitermass planet at 5.2 au. Moving from left to right: simulations M, M∕2 and M∕4 (i.e. decreasing disc masses). Note that the curves corresponding to different times have been vertically offset for clarity. 

Open with DEXTER 
Fig. 4
Azimuthallyaveraged radial surface density profiles at different times for nominal disc mass (M) simulations with a Jupitermass planet at 5.2 au. Left panel: same as the left panel of Fig. 3 for simulation M_{2D}, except that the scale on the yaxis is logarithmic, to appreciate the gap’s depth and its evolution with time. Right panel: same for the M_{3D} simulation. The gaps in 2D and 3D inviscid discs (with adiabatic EoS and a cooling timescale equal to one planet’s orbital period) appear very different; see text for discussion. 

Open with DEXTER 
Fig. 5
Left panels: evolution of semimajor axes and eccentricities for planets in the simulations M_{2D} sg, M_{2D} ∕2sg and M_{2D}∕4sg. Right panels: same as left panels except for the corresponding nonselfgravitating disc models. 

Open with DEXTER 
5.1 2D model
The evolution of the semimajor axes and the eccentricities are shown in Fig. 5 for the selfgravitating and nonselfgravitating M, M∕2 and M∕4 models. We begin by noting that the evolution of these quantities is similar in the selfgravitating and nonselfgravitating runs, which supports the conclusion of Sect. 4.1 that selfgravity is unimportant in this disc model. From now on we just discuss in detail the results of the nonselfgravitating nominal model in this section.
Figure 5 shows that after an initial adjustment phase when the planet is released from its fixed circular orbit, the planet migrates inwards at an approximately steady rate of ~ 5.55 au Myr^{−1}. During this time the eccentricity rises steadily up to e_{p} ~ 0.25. As discussed in Sect. 4.1, a vortex forms at the outer gap edge as the planet grows and opens a gap, but the vortex is short lived and dissipates shortly before the planet is released. Hence, the planet is embedded in a disc in which there is essentially no internal angular momentum transport, except due to the density waves excited by the planet itself.
A gapforming planet in a viscous disc can experience eccentricity growth when the gap becomes deep enough thateccentricity driving by external Lindblad resonances (ELRs) dominates the combined eccentricity damping due to coorbital Lindblad resonances (CLRs, which sit at the location of the planets emimajor axis, a_{p}) and corotation resonances (CORs; Goldreich & Tremaine 1980). Diagrams showing the locations of the various competing resonances are shown in Goldreich & Sari (2003). To first order in the planet eccentricity, e_{p}, and considering the outer disc only, we note that the 1:3 ELR (located at r = 2.08 a_{p}) is the outermost eccentricity driving resonance, and the 1:2 COR (located at 1.58 a_{p}) is the outermost damping resonance. Growth from an initially very small (formally infinitesimal) eccentricity can occur when the gap is wide enoughthat it extends beyond the 1:2 COR, such that the 1:3 ELR can drive up the eccentricity without competition from CORs (Papaloizou et al. 2001). Alternatively, eccentricity growth can occur when the initial eccentricity has a significant finite amplitude, because the CORs can then become partially saturated, and hence unable to dominate the ELRs (Goldreich & Sari 2003; Duffel & Chiang 2015).
In an inviscid disc, with no internal transport of angular momentum, CORs saturate completely and the formation of a deep gap renders the CLRs ineffective. Thus eccentricity growth should be expected, as shown in the simulations. Figure 6 depicts the surface density profile in the disc at two different times; also indicated are the planet semimajor axes and the locations of the 1:2 COR and the 1:3 ELR. At t = 39138 yr, when the eccentricity has already undergone significant growth, the gap is insufficiently wide to deplete the gas at the 1:2 resonance, suggesting this resonance is saturated due to the absence of viscosity or any other source of angular momentum transport such as Reynolds stress caused by a vortex.
The growth of the eccentricity gives rise to the higheccentricity mode of migration discussed above. It is clear from comparing Figs. 6 and 4 that the gap structure evolves very differently when the planet’s orbit is free to change versus when it is kept fixed and circular. In particular, the gap becomes less deep but wider as the planet moves away from the original outer edge of the gap while at the same time becoming eccentric. The eccentric orbit induces a flow of gas from the inner to the outer disc, and this appears to be sufficient to allow the planet to migrate inwards at the rate indicated by Fig. 5 (meaning ~ 6.7 au per Myr) over large radial scales and long time scales. It is noteworthy that the maximum eccentricity of e ~ 0.25 obtained by the planet is much larger than that reported by Duffel & Chiang (2015), who demonstrated that a Jovianmass planet in a viscous disc with initial eccentricity e = 0.04 could experience growth only up e = 0.07 before interaction with the gap edges stalls the eccentricity growth. Here, it appears that a Jovian mass planet in an inviscid disc can open a wider gap because torques exerted on the disc at ELRs as the eccentricity grows operate unopposed by viscosity. Hence, a larger eccentricity can develop.
Fig. 6
Surface density profile at the times shown in the legend for the run M_{2D}. Also shown are the semimajor axes of the planet (filled circles), and the locations of the 2:1 corotation resonances and of the 3:1 external Lindblad resonances at the times indicated. 

Open with DEXTER 
Fig. 7
Evolution of semimajor axis and eccentricity for a migrating planet in simulation M_{3D}. Notice the outward migration after 45 000 yr. The colored dots give the values of (a, e) at times corresponding to the panels of Fig. 8. 

Open with DEXTER 
5.2 3D model
The planet is released after 800 orbits at r = 5.2 au, i.e. at the time corresponding approximately to the middle bottom panel in Fig. 2. The planet remains on an orbit with small eccentricity (e ~ 0.01) because (i) the gap is not very deep and (ii) the CORs remain unsaturated because of the meridional circulation and the presence of the vortex, which generates an effective viscosity by exciting Reynolds stresses. Because of the small planet eccentricity, the passage of gas through the gap is virtually null. Thus, the migration mode observed in the previous 2D simulations is not observed. Nevertheless, the planet migrates inwards during the first ~ 40 000 yr (Fig. 7). Migration is possible because the vortex itself tends to migrate inwards, thus refilling the gap as the planet moves away. Vortex migration is due to its ability to exchange angular momentum with the disc, via the generation of spiral density waves (Paardekooper et al. 2010). Vortex migration is inwards because the vortex is located at the outer edge of the gap, so that the disc inside the vortex’s orbit is strongly depleted and the interaction with the outer disc dominates. By the actionreaction principle, the vortex transfers angular momentum to the disc beyond its orbit, carving a secondary gap (Fig. 8). Once this gap is formed, the vortex slows down, spreads radially and eventually dissipates completely, leaving the outer disc partially depleted near the planet’s orbit.
The depletion of the outer disc relative to the inner disc creates an imbalance in the torques exerted on the planet, with the positive torque exerted by the inner disc becoming the (slightly) dominant one (Fig. 9, top panel). As a consequence, the direction of migration is reversed and the planet slowly migrates outwards (Fig. 7, top). However, the process of outward migration cannot continue indefinitely. The torque contribution from the inner disc decreases when the planet moves outwards (because its distance from the inner disc increases), eventually balancing out with the negative torque exerted by the outer disc. Consequently, migration slowly stops and the planet remains on a quasicircular orbit embedded in a gap that is much wider than the initial one, resulting from the merging of the planet’s main gap with the secondary gap carved by the nowdisappeared vortex (Fig. 8, rightmost panel). Incidentally, this breaks the relationship usually used to deduce the planet’s mass from the gap’s width in observed disks. Figure 9, central panel shows that when the planet is at rest the gap becomes much deeper than the one during the migration phase. Nevertheless, enough gas remains in the gap and the CORs are sufficiently desaturated by the meridional circulation to damp the eccentricity of the planet once it stops migrating (Fig. 7, bottom).
In order to study more precisely the contributions of different regions of the disc – and of the vortex in particular – to planet migration, we show in Fig. 9 (bottom panel) the specific torque γ(r) exerted on the planet by a ring of the disc placed at r with width dr which is: (2)
where the subscript z indicates the vertical component of the vector enclosed in the brackets. The integral of γ(r) over the radial coordinate provides the total specific torque Γ, plotted in the top panel of Fig. 9: (3)
The vortex provides periodically positive to negative contributions to γ depending on the planetvortex azimuthal distance that average to zero over a planetvortex synodic period. To highlight the vortex location, its radial migration and its fading over time, it is convenient to select snapshots at which the planet and the vortex have the same relative shift in azimuth. We did this by sampling the hydrodynamical quantities every 1/10 of orbit and choosing snapshots with the vortex always ahead of the planet’s location as shown in Fig. 8, corresponding to the maximal value of γ at the vortex radial location during the vortex’s synodic period. The fact that γ is positive in these snapshots does not imply that the vortex exerts a net positive torque on the planet (remember that averaged over a synodic period the torque is zero). With this procedure, the bottom panel of Fig. 9 reveals that the position of the maximum in γ, corresponding to the vortex’s radial position, is located at the gap’s outer edge (middle panel). Both vortex and gap’s edge shift inwards with time, until t = 35 557 yr. Moreover the maximum value of γ decreases with time, revealing the weakening of the vortex. For times larger than t = 88 000 yr the vortex has been completely dissipated (2 rightmost panels of Fig. 8) and accordingly there is no positive contribution to γ.
In summary, after an initial short (~40 000 yr) phase of relatively fast inward migration, during which the planet migrates less than 0.5 au (Fig. 10), the migration of the planet slows down, then reverses direction, leading to an outward migration that naturally damps out at 4.85 au. A test of convergence with respect to the resolution of the simulation is provided in the Appendix.
When we decrease the mass of the disc, the vortex is less massive, roughly proportionally to the disc surface density. However, the indirect term also scales with the disc mass and therefore we cannot expect a simple linear decrease of the migration speed with the disc mass. In the M_{3D}∕2 and M_{3D}∕4 simulations the planet migrates slowly inward at a speed of about 1 au Myr^{−1} (slightly faster for M_{3D}∕2 than for M_{3D}∕4). In both cases we observe the vortex spreading radially and dissipating with time while opening a small dip beyond the vortex orbit. The torque imbalance between inner and outer disc appears to be not enough for the reversal of the migration direction.
Fig. 8
Contours of surface density for the M_{3D} simulation and a migrating Jupitermass planet. The panels correspond to different times, reported on top of each panel, corresponding to the dots in Fig. 7. It appears clearly that a secondary gap is carved by vortexdisc interaction (second and third panels) and that the vortex weakens by spreading radially, leaving the outer disc partially depleted (fourth and fifth panels). 

Open with DEXTER 
Fig. 9
Top panel: evolution of specific torques from the inner (blue curve) and outer disc (red curve) and of the total specific torque (cyan) in simulation M_{3D}. A runningaverage has been applied over 500 orbits to smooth out the oscillations of the torque exerted by the outer disc. Such oscillations are due the synodic period of the vortex with respect to the planet. At times corresponding to the dots in the top panel (the same times selected in Figs. 7 and 8) we plot: (middle panel) surface density profiles; (bottom panel) specific radial torque density γ(r) (see Eq. (2) in the text). 

Open with DEXTER 
Fig. 10
Migration speed as a function of time, for discs with different masses in 3D simulations. The planet always starts at 5.2 au. The migration speed decreases, as expected, with the disc mass. The direction of migration is reversed at t ≃ 45 000 yr in the nominal mass disc case (simulation M_{3D}). 

Open with DEXTER 
5.3 Definition of vortexdriven migration
The result illustrated above shows that, in inviscid disks, migration of giant planets on quasicircular orbit is sustained because the migration of the vortex allows to refill the gap leftbehind by the migrating planet. We call this process ‘vortexdriven’ migration, to distinguish it from the ‘vortexinduced’ migration of Lin & Papaloizou (2010, 2011a) which was due to the passage of the vortex through the planet’s coorbital region. As we mentioned in the introduction, vortex migration is much less prominent if selfgravity is taken into account (Zhu & Baruteau 2016). Indeed, in Appendix A.1 we will show a 2D simulation with selfgravity and fast disc cooling that will behave similarly to our M_{3D} simulation, i.e. with a pronounced vortex and a slowly migrating planet on a quasicircular orbit, but showing no vortex migration. In these cases, planet migration is possible because the vortex, although not moving, is able to push gas inwards, refilling the gap leftbehind by the planet. We include also this case in the ‘vortexdriven’ migration category. In both cases, the migration of the planet occurs on the timescale characterising the vortexdisc interaction (either manifested by vortex migration or vortexdriven disc spreading), which is unrelated to the disc viscosity. Thus vortexdriven migration is conceptually different from TypeII migration. Vortexdriven migration is dominated by TypeII migration if the disc viscosity is large enough (as quantified in the next section) or superseded by the eccentric migration mode (as seen in Sect. 5.1) if the planet’s eccentricity becomes large enough to allow the passage of gas across the gap, unlocking the planet from the disc evolution.
In summary, vortexdriven migration is a slow, inward and smooth migration due to the presence of a vortex beyond the orbit of the planet. It is very different from vortexinduced migration, which is a violent process inducing a phase of fast migration. Importantly, vortexdriven migration does not seem to be indefinite. By modifying the disk around it (i.e. spreading the disk inwards, opening a secondary gap) the vortex eventually fades away and its effect stops. When this happens, a transient phase of outward migration is even possible if the disappearance of the vortex leaves a gap that is asymmetric relative to the position of the planet.
Vortexdriven migration, being slow and transient, offers a new and interesting explanation of the existence of the coldJupiters, without requiring (as in Coleman & Nelson 2014; Bitsch et al. 2015) that these planets formed much beyond the snowline.
5.4 Migration in viscous discs
This sectionaddresses the question: how low does the viscosity need to be for the vortexdriven migration to be dominant? For this purpose we ran 3D simulations with values of the α viscosity parameter in the range [10^{−5}, 10^{−3}]. In viscous discs the equation of state should consistently include a term corresponding to viscous heating. However, the viscous heating would also change the vertical disc structure (Kley et al. 2009) making the comparison less clear. In order to compare equivalent discs we null the viscous heating contribution.
Figure 11 shows the migration speed as a function of the planetstar distance. Planets monotonically migrate towards the star for α down to 10^{−4}, with a speed which scales with the viscosity, as expected in the TypeII migration regime. For α =10^{−5} a vortex forms very similarly to the nonviscous case and vortexdriven migration is observed with speed slightly larger than in the nonviscous case and persisting for a longer time, thus driving the planet closer to the star (to ~ 4.55 au instead of ~ 4.75 au at the moment of migration reversal). Therefore, we conclude that the transition between TypeII migration and vortexdriven migration occurs for values of α in the interval [10^{−5}, 10^{−4}]. This makes vortexdriven migration an interesting process acting at viscosities that are not exceedingly small, but could be typical of realistic discs (Turner et al. 2014).
When considering nonviscous discs in numerical simulations there is always a nonzero numerical viscosity. The fact that we see a different evolution for decreasing values of α down to 10^{−5} (Fig. 11) and that thesimulation with α = 10^{−5} is different from the inviscid simulation suggests that the numerical viscosity for our inviscid 3D simulation M_{3D} is not larger than the one modeledassuming α = 10^{−5}.
Fig. 11
Migration speed as a function of planet’s distance from the star in 3D simulations of discs with different viscosities (i.e. simulations M_{3D}, M_{3D α5}, M_{3D α4} and M_{3Dα3}; see Table 1). A running average procedure has been used to reduce fluctuations, with the consequence that all but the green curve start at about 5.1 au instead of the released distance of 5.2 au. In the α = 10^{−3} and 10^{−4} cases migration is always inwards (negative radial speed). Notice that we have divided by 10 the migration speed of M_{3D α3} in order to scale it proportionally to viscosity relatively to the α = 10^{−4} case. The migration pattern in the inviscid case does not respect this scaling and turns positive after a short phase of inward migration. The same happens in the M_{3D α5} simulation. 

Open with DEXTER 
6 Distant planets
In some planetevolution models, giant planets are postulated to start forming at large distances from the parent star. From the results of this paper, one could be tempted to conclude that in lowviscosity discs these planets remain at large distances throughout the disc’s lifetime and may correspond to the giant planets observed by direct imaging surveys. Here, we warn that the results are not invariant with the distance of the planet from the central star. This is because the role of selfgravity becomes more prominent farther out.
6.1 Fixed planet at r_{p} = 13 au, 2D model
We start this discussion by considering the M_{2D13au} and M_{2D13au}sg simulations. The local Toomre parameter Q(r = 13 au) = 7.74 and the aspect ratio h(r = 13 au) = 0.065 for this model (see Fig. 1). We plot in Fig. 12 contours of the surface density.
As a consequence of the relatively low Q value at the planet’s location, a behaviour similar to that described in Zhu & Baruteau (2016) is replicated in our simulations. The selfgravitating disc (M_{2D13au}sg) forms a single vortex that remains stably orbiting at the planetinduced pressure bump throughout the simulation and, although a detailed inspection of the surface density evolution shows some intermittency in the structure of the vortex and the pressure bump, the evolution is in general smooth and involves the vortex weakening with time. On longer time scales, however, the disc develops a strong m = 1 eccentric mode, that eventually becomes so strong that it destroys the disc (remember that the planet is kept on a circular orbit, so that the disc’s eccentricity cannot be limited by transferring the disc’s angular momentum deficit to the planet). The nonselfgravitating run M_{2D13au} displays very different behaviour. A large vortex forms and migrates into the gap region and disperses, smoothing the profile of the gap outer edge. Over long run times (i.e. >200 000 yr, not shown in Fig. 12), we observe that nonaxisymmetric structures form again in the disc.
To reach genuine convergence between selfgravitating and nonselfgravitating disc in the case of distant planets we need to reduce the disc mass by a factor 10 (see Fig. 13), which brings Q to ~ 70, comparable to the value in the M∕2 simulations for r_{p} = 5.2 au. In this case the disc forms a single vortex at the outer edge of the gap that remains there over the simulation run time. The vortex weakens over time, and the only significant difference between the selfgravitating and the nonselfgravitating runs is the length of time required for vortex weakening to occur. Because of the lack of convergence in the disc behaviour for selfgravity and noselfgravity runs we do not report a detailed analysis of migration, but it is worth mentioning that in the M_{2D13au}sg simulation the planet becomes very eccentric (unsurprisingly, given the eccentricity acquired by the disc), which makes the planet go out of the gap at perihelion and aphelion, triggering very fast inward migration.
6.2 Fixed Planet at r_{p} = 13 au, 3D model
We remind the reader that our 3D code does not feature selfgravity. Thus, it is not a surprise that the disc surface density (Fig. 12, bottom panels) appears more similar to the 2D nonselfgravitating model (middle panels) than to the selfgravitating case (top panels). The vortex interacts with the disc and weakens while forming a secondary gap. At the outer edge of the secondary gap a second vortex forms while the inner vortex gets weaker. Reducing the disc mass by a factor 10, we have obtained results quite similar to the 2D models, with a weak vortex that appears to evolve very slowly over time.
It is then clear that in order to study giant planet migration realistically at large distances in nominal discs one needs to use a 3D code with selfgravity, which to our knowledge has never be done before. Without selfgravity, only migration in lowdensity discs can be studied and we find a similar behavior as that discussed in Sect. 5.2.
7 Discussion and conclusions
In this paper we have examined giant planet migration in low viscosity discs with 2D and 3D numerical simulations. First, we have considered a Jupitermass planet at 5.2 au in a nominal disc with density at 5.2 au of 220 g cm^{−2} or smaller. The presence of a giant planet in a low viscosity disc opens a gap and triggers vortex formation at the outer edge of the gap. The presence of the vortex raises the question of the importance of the selfgravity in the disc (Zhu & Baruteau 2016). Thus, we compared 2D simulations with and without disc selfgravity to conclude that, apart for minor and temporary differences, the evolution observed with a disc cooling time τ_{c} = 1 orbit are equivalent to each other. Then we moved to 3D simulations of the planetdisc interaction. 3D simulations are in principle more realistic than 2D simulations because they can capture genuine 3D effects such as the meridional circulation of gas in the vicinity of gaps (Morbidelli et al. 2014), but our code FARGOCA does not allow to take selfgravity into account, so that its use is valid only where selfgravity is not relevant, as shown by the aforementioned 2D simulations.
For the nominal model with τ_{c} = 1 orbit, the results of 3D simulations are quite different from those of 2D simulations: the gap is shallower and the vortex at its outer edge is much larger and persistent. The two aspects are linked in a positive feedback. Because the gap is shallower, the disc never becomes Rayleigh unstable, unlike in 2D simulations. Thus the vortex is not eroded by the diffusivity generated by this instability; in turn the presence of the vortex exerts a torque on the disc which helps to push gas into the gap. This effect, together with the meridional circulation of gas that exists only in 3D, makes the 3D gap shallower.
As a consequence of the sharp differences in the gap’s depth and vortex properties, we observe a big difference in the orbital evolution of the planet when the latter is free to migrate. In 2D simulations the planet rapidly acquires a significantly eccentric orbit with e ~ 0.25, apparently because the depth of the gap and the lack of effective viscosity operating in the disc allows coorbital and corotation resonances, that usually damp the eccentricity, to saturate (Goldreich & Tremaine 1980; Goldreich & Sari 2003). In their paper, Goldreich and Sari conjectured that this mechanism can explain the large orbital eccentricities observed for extrasolar coldJupiters. We find that this is unlikely. In fact, as soon as the planet’s orbit becomes eccentric, the planet undergoes a relatively fast migration (~ 6.7 au Myr^{−1} which, although slower than TypeII migration in an α = 10^{−3} disk, impliesfull orbital decay in less than a Myr relative to the planet’s initial position at 5.2 au). Thus, we suggest it is unlikely that an eccentric giant planet could remain in the coldJupiter region at the end of the disc’s lifetime, given the speed of migration. Moreover, we do not see such an eccentricity excitation in the 3D simulations, as discussed below. Thus, dynamical instabilities after the removal of the protoplanetary disc (Ford & Rasio 2008; Chatterjee et al. 2008; Jurić & Tremaine 2008) seem to provide a much more robust mechanism for the eccentricity excitation of giant planets.
In 3D simulations the planet remains on a quasicircular orbit, probably because there is still enough gas in the gap to damp continuously its eccentricity. The meridional circulation may also play a role in unsaturating corotation resonances. The vortex tends to migrate inwards by exerting a torque on the outer disc and therefore can accompany the planet in its inward migration. Thus, initially planet migration is inward and occurs at the speed of vortex migration, which is unrelated to the viscosity of the disc (Paardekooper et al. 2010). We call this ’vortexdriven’ migration (see Sect. 5.3). However, because the disc is inviscid, the torque exerted by the vortex opens a secondary gap just beyond the vortex orbit. The vortex then spreads radially and becomes less pronounced. These processes partially deplete the disc just beyond the planet’s gap and, consequently, the planet starts to feel a positive torque from the inner disc that exceeds the negative torque from the outer disc: its migration is reversed. By being pushed out, eventually the planet adjusts its position relative to the inner and outer discs to balance the two torques and its migration stops. If the mass of the disc is reduced, the relative strength of the vortex is also reduced. Consequently, no significant secondary gap is opened; migration is not reversed, but it slows down to a rate smaller than 1 au My^{−1}; it should eventually stop on the longterm.
In both cases, these results can potentially explain the lack of largescale migration of giant planets. If inward migration is shortranged and eventually slows down, or reverses and stops, the presence of numerous giant planets beyond 1 au from the central star (i.e. coldJupiters) becomes consistent with the idea that the sweetspot for giant planet formation is the disc’s snowline, which typically is in the same region or only slightly beyond. This result is very promising to solve the giant planet migration problem, but it needs to be confirmed using more realistic discs. In fact, the disc we have adopted with very low viscosity and therefore negligible radial transport would be incompatible with the mass accretion rates that are typically observed for young stars. If the protoplanetary discs have low viscosities because of the lack of magnetorotational instability, most likely the accretion of gas towards the star occurs due to laminar stresses where angular momentum is removed in disc winds (see Turner et al. 2014 for a review). In a future work we will study giant planet migration in these discs. Nevertheless, the results presented in this paper open a new perspective and will serve as a basis for comparison.
In the final part of the paper we have also stressed that the result of slow, or reversed, giantplanet migration for low viscosities holds only in the inner part of the disc. The outer part of the disc is closer to the gravitationally instability limit and therefore can respond very differently to the presence of a giant planet. Unfortunately, our 3D code cannot account for selfgravity in its present form. A future work will be devoted to performing 3D selfgravitating simulations. Here, by using 2D selfgravitating simulations we observed that the disc may develop a strong m = 1 eccentric mode, which in turn excites the eccentricity of the planet, triggering fast inward migration. This result suggests that giant planet migration may be fast in the outer disc and slow in the inner disc, so that giant planets, wherever they form, should all pileup in the coldJupiter region. If this is true, we would then predict a strong dropoff in the giant planet radial distribution, with only the most massive planets – itive to planetdisc interactions – remaining at large distances where they can be probed by direct imaging surveys. The stakes are high, so it is important to confirm or improve on these preliminary conclusions with 3D selfgravitating simulations.
Fig. 12
Same as Fig. 2, but for a planet located at r_{p} = 13 au. 

Open with DEXTER 
Fig. 13
Azimuthallyaveraged radial surface density profiles versus time for 2D runs with and without selfgravity, and a Jupitermass planet at 13 au. Moving from left to right: simulations M∕2, M∕4 and M∕10. Note that the curves corresponding to different times have been vertically offset for clarity. 

Open with DEXTER 
Acknowledgements
We acknowledge support by DFGANR supported GEPARD project (ANR18CE920044 DFG: KL 650/311). We also acknowledge HPC resources from GENCI DARI n.A0080407233 and from “Mesocentre SIGAMM” hosted by Observatoire de la Côte d’Azur. L.E. wish to thank Alain Miniussi for maintainance and refactorization of the code FARGOCA. R.P.N. acknowledges support from the STFC through the Consolidated Grants ST/M001202/1 and ST/P000592/. This research utilised Queen Mary’s Apocrita HPC facility, supported by QMUL ResearchIT. T.R. acknowledges support by the DFG research group FOR2634: “Planet Formation Witnesses and Probes: transition Disks” under grant KL650/291.
Appendix A Sensitivity of results to the cooling time
A.1 2D simulations
In addition to running simulations for cooling time τ_{c} = 1 orbit, we also ran simulation suites for selfgravitating and nonselfgravitating models with τ_{c} = 0.01 and τ_{c} = 10 orbits. We considered the disc masses of the M, M∕2 and M∕4 simulations. Our results indicate that for 2D simulations, the detailed outcome depends on the cooling time in a rather complicated manner.
Fig. A.1
Left panels: evolution of semimajor axes and eccentricities for planets in simulations M_{2D}, M_{2D} ∕2 and M_{2D}∕4 with cooling time τ_{c} = 0.01 orbit. Right panels: same as left panels except for the corresponding nonselfgravitating disc models. 

Open with DEXTER 
A.1.1 τ_{c} = 0.01 orbits
Figure A.2 shows the surface density evolution for the M_{2D} simulation with cooling time of 0.01 orbits. The top panels are for the selfgravitating run and the bottom panels are for the nonselfgravitating disc. This figure should be compared with Fig. 3 in the main text, which shows the results for τ_{c} = 1 orbit. The shorter cooling time gives rise to different qualitative behaviour of the disc and the embedded planet, which is released from a fixed orbit at the time corresponding to the leftmost panels in Fig. A.2. The selfgravitating run shows the formation of a strong vortex at the outer edge of the gap, and this remains relatively narrow in its radial extent and stays at the gap edge without migrating. Over time the vortex weakens and dissipates, and appearsto have almost completely dispersed in the rightmost panel of Fig. A.2. Analysis of the disc at late times, extending beyond the time of the rightmost panel of Fig. A.2 to ~ 500 000 yr, indicates that a vortex remains at the pressure bump at the outer edge, but with a strength that waxes and wanes intermittently throughout the run. The nonselfgravitating run results in the formation of a more radially extended vortex, and the evolution here is much more similar to the 3D results described in Sect. 5.2. The radial width of the vortex is large enough for the excitation of spiral waves to drive its migration into the gap, while at the same time the vortex creates a secondary gap (lower panels, second from left). The third panel shows that the formation of a secondary gap also leads to the formation of a secondary pressure bump, and this forms a second vortex via the Rossby wave instability. Eventually, the primary vortex disperses, leaving the weaker secondary vortex orbiting at the gap edge (fourth panel). This vortex persists for the whole simulation. Hence, we see that in these selfgravitating and nonselfgravitating runs, the initial evolution of the system is quite different because of the formation ofa large, migrating vortex in the latter simulation, but the long term behaviour is actually very similar.
Fig. A.2
Contours of surface density for simulation M_{2D} with τ_{c} = 0.01 orbits and a migrating Jupitermass planet. The panels correspond to different times, reported on top of each panel. Top panels: selfgravitating discs. Bottom panels: nonselfgravitating discs. 

Open with DEXTER 
We do not show figures for the M_{2D}∕2 and M_{2D}∕4 runs, but instead comment that the qualitative differences between selfgravitating and nonselfgravitating models are much smaller than shown in Fig. A.2. Reducing the disc mass reduces the importance of the indirect term, and hence the radial width of the vortex that forms is smaller in these cases. Hence, the outer gapedge vortex does not migrate quickly into the gap, but remains at the gap edge and slowly dissipates while generating a shallow secondary gap. The main difference between the selfgravitating and nonselfgravitating models is that the vortex weakens more quickly in the selfgravitating models. The vortex in the nonselfgravitating models is persistent until the ends of the simulations, providing a source of weak effective viscosity via the generation of Reynolds stresses, whereas the vortex in the selfgravitating discs appears to be more intermittent in its behaviour.
Fig. A.3
Contours of surface density for simulations M_{2D} with τ_{c} = 10 orbits and a migrating Jupitermass planet. The panels correspond to different times, reported on top of each panel. Top panels: selfgravitating discs. Bottom panels: nonselfgravitating discs. 

Open with DEXTER 
Fig. A.4
Left panels: evolution of semimajor axes and eccentricities for planets in selfgravitating disc simulations M_{2D}, M_{2D} ∕2 and M_{2D}∕4 where τ_{c} = 10 orbit. Right panels: same as left panels except for the corresponding nonselfgravitating disc models. 

Open with DEXTER 
Figure A.1 shows the evolution of semimajor axes and eccentricities for the M_{2D}, M_{2D}∕2 and M_{2D}∕4 runs with τ_{c} = 0.01. The selfgravitating models show the planet initially adjusts its orbital location in the gap due to an imbalance of torques, and experiences a phase of eccentricity growth due to direct interaction with the vortex. The eccentricity, however, is eventually damped by interaction with the disc as the vortex weakens, and remains at a very low value (~ 5 × 10^{−3}) for the remainder of the run. During this low eccentricity phase, we see that the migration rate is very slow, corresponding to migration over a distance of 1 au in 4 × 10^{6} yr (i.e. a migration speed of 0.25 au per Myr) for the nominal model. The orbital evolution during the M_{2D} simulation is initially more extreme because of the strong interaction with the migrating vortex, whereas the evolution in the M_{2D}∕2 and M_{2D}∕4 simulations are similar to the selfgravitating models. Over the long term, however, we see that in these models the planet also exhibits low eccentricities and very slow migration rates. Hence, the τ_{c} = 0.01 selfgravitating and nonselfgravitating models during longterm evolution exhibit the low eccentricity mode of migration observed in the M_{3D} simulations described in Sect. 5.2. This maintenance of low eccentricity in an inviscid disc is not expected, since the corotation resonances that are responsible for damping eccentricity are expected to saturate and switch off. The persistence of vortices at the outer gap edge, however, seems to provide a source of effective viscosity that allows these resonances to actively damp the eccentricity.
A.1.2 τ_{c} = 10 orbits
Figure A.3 shows the surface density evolution for the M_{2D}sg (top panels) and M_{2D} models (bottom panels). The selfgravitating disc shows similar behaviour to that seen previously for both τ_{c} = 0.01 and τ_{c} = 1, namely the formation of a radiallynarrow vortex that remains at the gap edge while weakening and eventually dissipating due to the onset of the Rayleigh instability. For τ_{c} = 1 and 10 the vortex dissipates much more quickly than in the τ_{c} = 0.01 run. The nonselfgravitating run shows the formation of a radiallyextended vortex which forms a significant secondary gap and migrates into the planetinduced gap before dispersing. This occurs relatively quickly in the simulation, and leaves behind a smooth outer disc without any evidence for additional vortices being present (unlike in the τ_{c} = 0.01 run). The evolution of the discs for the M_{2D}∕2 and M_{2D}∕4 models show much stronger convergence with respect to whether selfgravity is included or not. Here the vortex at the outer edge is radially narrow and does not migrate into the gap and strongly interacts with the planet. Instead it sits at the gap edge in all cases, weakening with time and eventually dissipating.
Fig. A.5
Contours of surface density for the M_{3D} simulations with the locally isothermal equation of state and for increasing values of τ_{c} (in orbits). All the panels correspond to t =8296 yr (the same time as for the {middle panels} of Fig. 2, for comparison). 

Open with DEXTER 
Figure A.4 shows the evolution of semimajor axes and eccentricities for all runs with τ_{c} = 10. Apart from the initial phase of evolution of the M_{2D}sg case, where the planet interacts strongly with the vortex, we see that these models all display the higheccentricity mode of migration described in the main text for the models with τ_{c} = 1. The mean migration rates vary somewhat between the runs, but correspond to migration over a distance of 1 au in 150 000–400 000 yr. Hence, Jovian mass planets in these models can migrate to the star within typical disc lifetimes.
A.2 3D simulations
In addition to the Nominal 3D case with τ_{c} = 1 we ran simulations with fixed planet only for the locally isothermal equation of state and for τ_{c} = 10 and τ_{c} = 50 orbits.
We can consider the locally isothermal case equivalent to having an extremely short cooling time, so that in Fig. A.5 we show results (from left to right panels) from short cooling times to high values such as τ_{c} = 50 orbits.
In all cases a big vortex forms at the outer edge of the gap and qualitatively all the panels are very similar. It appears that 3D models are much less sensitive to cooling times than 2D ones. We remark that for large values of the cooling time, the vertical temperature profile is not damped efficiently toward the initial condition: h_{0} = 0.05. Actually, for τ_{c} = 50 orbital periods the aspect ratio fluctuates around a mean value of 0.057 at the time snapshot of Fig. A.5. Qualitatively there is no influence on vortex formation but a more detailed study of the role of the thermal structure on vortex formation and planet migration goes beyond the scope of the present paper.
Due to the expensive computation required for 3D migration, we decided to limit our analysis of the sensitivity to the cooling time to fixed planets. However, the similarity in vortex formation indicates that the migration histories would probably be very similar to the case τ_{c} = 1 orbit described in Sect. 5.
Appendix B Numerics
We have seen above that vortexdriven migration occurs if gaps are only moderately deep, so it is important to check the dependence of results with respect to the choice of the smoothing of the planetary potential. Moreover, vortex formation and evolution may depend on numerical resolution. It is therefore important to check also the validity of results with respect to the grid stepsize. We present these tests below.
In all the simulations we smooth the potential of the planets for disc elements with distance d from the planet: d < ϵ. We have considered ϵ = 0.8 R_{H} in our nominal M_{3D} simulations. In order to check the sensitivity to the smoothing length we have run a simulation at nominal resolution with ϵ = 0.4 R_{H} for 800 orbits.
Figure B.1 shows that the surface density profiles for the nominal simulation (blue curve) and the new one with reduced smoothing length (red curve) at the time of planet release. As one can see, the two curves are essentially identical. Thus, we can conclude that the limited depth of the gap is not due to an excessive smoothing but corresponds to the real dynamic of the gas.
When the planet is released (Fig. B.2, red curve) the migration is very similar, although slightly faster than the nominal simulation. Importantly, we see the same reversal of migration direction, although this happens slightly earlier with the planet being slightly closer to the Sun.
Testing the effect of resolution is hard because of the long integration times. Thus we have restarted the simulation with ϵ = 0.4 R_{H} with double resolution at t = 10^{4} yr when the planet is still on a fixed orbit. We recall that the nominal M_{3D} simulation had a resolution (N_{r}, N_{θ}, N_{φ}) = (568, 16, 360). In the high resolution simulations we reduced further the smoothing length to ϵ = 0.2 R_{H}. The surface density profile obtained after 100 orbits is shown by the green curve in Fig. B.1 and is extremely similar to those of the two other simulations. We then released the planet. On a timescale of about 800 orbits the migration rate is indistinguishable from that of the previous simulations (Fig. B.2). Then, at time ~ 19 000 yr we observe an episode of rapid inward migration (and eccentricity increase up to 0.07), before turning to slow outward migration as in the nominal resolution case. The outward migration phase starts at t = ~23 000 yr when the vortex is weakened by radial spreading and the eccentricity is damped to ~ 0.001. Finally, planet’s migration stops at a position in which the torques from the inner and outer discs balance. From this test we conclude that in the non viscous case the dynamics described in the main paper is not an artefact of resolution, although the evolution of the planet can change at the quantitative level using different resolutions. To obtain a quantitativeconvergence of results with respect to resolution and smoothing length a low prescribed viscosity may be required.
To do this we have also run the double resolution test by considering the low viscosity α =10^{−5} simulation M_{3Dα5} for which we have observed vortexdriven migration (see Fig. 11). We have restarted the simulation with α =10^{−5} by doubling the resolution at t ~ 16 000 yr when theplanet is already in the migration phase. This choice is motivated on the one hand by computational cost (about 1 orbit at 5.2 au per hour using 240 cores with Hybrid OpenMP+MPI parallelization), on the other hand by the results obtained with double resolution in the inviscid case. In fact, as stated above, up to t ~ 19 000 yr migration doesn’t show any dependency on resolution (Fig. B.2).
Fig. B.1
Surface density profiles at t = 10000 yr for standard resolution non viscous 3D simulations with different smoothing length and for an additional run with double resolution and smoothing length ϵ = 0.2 Hill radii. We call “SR” the standard resolution simulations and we use the label “HR” for the high resolution one. The "HR" simulation is a restart of the “SR” case with ϵ = 0.4 R_{H} (red curve) at t = 10^{4} yr. The green curve corresponds to 100 orbits after the restart. 

Open with DEXTER 
Figure B.3 shows that the migration rate is not affected by the change in the resolution (we did not observe any increase in eccentricity either). We have integrated this case up to t ~ 23 000 yr. We considerthis integration time (corresponding to the end of the fast inward migration phase observed for the non viscous case) long enough to prove quantitative robustness of results with respect to resolution and smoothing length.
Fig. B.2
Evolution of the planet semimajor axis for the SR M_{3D} simulation with different smoothing lengths and for the additional HR simulation with smoothing length ϵ = 0.2 Hill radii. 

Open with DEXTER 
Fig. B.3
Evolution of the planet semimajor axis for the SR M_{3D} simulation with viscosity α = 10^{−5} and for the additional HR simulation with smoothing length ϵ = 0.2 Hill radii. 

Open with DEXTER 
We conclude that a small prescribed viscosity (α = 10^{−5}) is necessary to show convergence of results with numerical resolution also at the quantitative level, confirming the results of the paper.
References
 Adams, F. C., & Gregory, S. G. 2012, ApJ, 744, 55 [Google Scholar]
 Bai, X., & Stone, J. M. 2013, ApJ, 769, 76 [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [Google Scholar]
 Barranco, J. A., & Marcus, P. S. 2005, ApJ, 623, 1157 [Google Scholar]
 Baruteau, C., & Masset, F., 2008, ApJ. 678, 483–97 [Google Scholar]
 Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., Lambrechts, M., & Johansen, A., 2015, A&A, 582, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Butler, R. P., Wright, J. T., Marcy, G. W., et al. 2006, ApJ, 646, 505 [Google Scholar]
 Chatterjee, S., Ford, E. B., Matsumura, S., et al. 2008, ApJ, 686, 580 [Google Scholar]
 Chiang, E. I., & Goldreich, P., 1997, ApJ., 490, 368 [Google Scholar]
 Coleman, G. A. L., & Nelson, R. P., 2014, MNRAS, 445, 479 [Google Scholar]
 Cumming, A., Butler, R. P., Marcy, G. W., et al. 2008, PASP, 120, 531 [CrossRef] [Google Scholar]
 de ValBorro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529 [Google Scholar]
 Drążkowska, J., & Alibert, Y. 2017, A&A, 608, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Duffel, P. C., & Chiang, E. 2015, ApJ, 812, 94 [Google Scholar]
 Duffell, P. C., Haiman, Z., MacFadyen, A. I., et al. 2014, ApJ, 792, L10 [Google Scholar]
 Dürmann, C., & Kley, W. 2015, A&A, 574, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fernandes, R. B., Mulders, G. D., Pascucci, I., et al. 2019, ApJ, 874, 81 [Google Scholar]
 Ford, E. B., & Rasio, F. A. 2008, ApJ, 686, 621 [Google Scholar]
 Fung, J., & Chiang, E. 2017, ApJ, 839, 100 [Google Scholar]
 Fung, J., Masset, F., Lega, E., & Velasco, D. 2017, AJ, 153, 124 [Google Scholar]
 Gammie, C. F. 1996, ApJ, 457, 355 [Google Scholar]
 Goldreich, P., & Sari, R. 2003, ApJ, 585, 1024 [Google Scholar]
 Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425 [Google Scholar]
 Gressel, O., Turner, N. J., Nelson, R. P., et al. 2015, ApJ, 801, 84 [Google Scholar]
 Hallam, P. D., & Paardekooper, S.J. 2020, MNRAS, 491, 5759 [Google Scholar]
 Hammer, M., Kratter, K. M., Lin, M.K. 2017, MNRAS, 466, 3533. [Google Scholar]
 Hartmann, L., Calvet, N., Gullbring, E., et al. 1998, ApJ, 495, 385 [Google Scholar]
 Hayashi, C. 1981, PThPS, 70, 35 [Google Scholar]
 Howard, A. W., Marcy, G. W., Johnson, J. A., et al. 2010, Science, 330, 653 [Google Scholar]
 Jurić, M., & Tremaine, S. 2008, ApJ, 686, 603 [Google Scholar]
 Kanagawa, K. D., Tanaka, H., Muto, T., et al. 2015, MNRAS, 448, 994 [Google Scholar]
 Kley, W., Bitsch, B., & Klahr, H. 2009, A&A, 506, 971 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Koller, J., Li, H., & Lin, D. N. C. 2003, ApJ, 596, L91 [Google Scholar]
 Laughlin, G., & Adams, F. C. 1997, ApJ, 491, L51 [Google Scholar]
 Lega, E., Crida, A., Bitsch, B., & Morbidelli, A. 2014, MNRAS, 440, 683 [Google Scholar]
 Lin, D. N. C. 1997, IAU Colloq., 121, 321 [Google Scholar]
 Lin, D. N. C., & Papaloizou, J. 1986, ApJ, 309, 846 [Google Scholar]
 Lin, M.K., & Papaloizou, J. C. B. 2010, MNRAS, 405, 1473 [NASA ADS] [Google Scholar]
 Lin, M.K., & Papaloizou, J. C. B. 2011a, MNRAS, 415, 1426 [Google Scholar]
 Lin, M.K., & Papaloizou, J. C. B. 2011b, MNRAS, 415, 1445 [Google Scholar]
 Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ 513, 85 [Google Scholar]
 Manara, C. F., Rosotti, G., Testi, L., et al. 2016, A&A, 591, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Masset, F. 2000, A&AS, 141, 165 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Masset, F. 2002, A&A, 387, 605 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Matt, S., & Pudritz, R. E. 2005, ApJ, 632, L135 [Google Scholar]
 Mayor, M., Marmier, M., Lovis, C., et al. 2011, ArXiv eprints [arXiv:1109.2497] [Google Scholar]
 McNally, C. P., Nelson, R. P., Paardekooper, S.J., & BenítezLlambay, P. 2019, MNRAS, 484, 728 [Google Scholar]
 Meheut, H., Casse, F., Varniere, P., & Tagger, M. 2010, A&A, 516, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mohanty, S., & Shu, F. H. 2008, ApJ, 687, 1323 [Google Scholar]
 Morbidelli, A., Szulágyi, J., Crida, A., et al. 2014, Icarus, 232, 266 [Google Scholar]
 Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610 [Google Scholar]
 Paardekooper, S.J., Lesur, G., & Papaloizou, J. C. B. 2010, ApJ, 725, 146 [Google Scholar]
 Papaloizou, J. C. B. Nelson, R. P., & Masset, F. 2001, A&A, 336, 263 [Google Scholar]
 Pfeil, T., & Klahr, H. 2020, ApJ, submitted [arXiv:2008.11195] [Google Scholar]
 Quillen, A. C., Blackman, E. G., Frank, A., et al. 2004, ApJ, 612, L137 [Google Scholar]
 Robert, C. M. T., Crida, A., Lega, E., Méheut, H., & Morbidelli, A. 2018, A&A, 617, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schoonenberg,D., & Ormel, C. W. 2017, A&A, 602, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33 [NASA ADS] [Google Scholar]
 Stoll, M. H. R., & Kley, W. 2014, A&A, 572, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stone, J. M., & Norman, M. L 1992, ApJS 80, 753 [Google Scholar]
 Stone, J. M., Hawley, J. F., Gammie, C. F., et al. 1996, ApJ, 463, 656 [Google Scholar]
 Turner, N. J., Fromang, S., Gammie, C., et al. 2014, Protostars and Planets VI (Tucson, AZ: University of Arizona Press), 411 [Google Scholar]
 Udry, S., Fischer, D., & Queloz, D. 2007, Protostars and Planets V (Tucson, AZ: University of Arizona Press), 685 [Google Scholar]
 Walsh, K. J., Morbidelli, A., Raymond, S. N., et al. 2011, Nature, 475, 206 [Google Scholar]
 Ward, W. R. 1997, Icarus, 126, 261 [Google Scholar]
 Weidenschilling, S. J. 1977, Ap&SS, 51, 153 [Google Scholar]
 Zhu, Z., & Baruteau, C. 2016, MNRAS, 458, 3918 [Google Scholar]
The simulations presented in this paper have been obtained with a recently refactorized version of the code that can be found at: https://disc.pages.oca.eu/fargOCA/public/
All Tables
All Figures
Fig. 1
Toomre Q values for 2D selfgravitating disc models. We focus on the nominal or M disc, but we also run control simulations with the mass of the disc divided by 2, 4 and 10. The Q = 1 value is indicated by the dotdashed line, and the locations of the planets at r = 5.2 and r = 13 au are indicated by the filled red circles. 

Open with DEXTER  
In the text 
Fig. 2
Top panels: contours of surface density for the M_{2D}sg simulation of the nominal disc, for times shown (in years) at the top of each panel. Middle panels: same as top panels but for the M_{2D} simulation (i.e. without selfgravity). Bottom panels: same as top panels but for the M_{3D} simulation. The planet is located at r = 5.2 au. 

Open with DEXTER  
In the text 
Fig. 3
Azimuthallyaveraged radial surface density profiles versus time for 2D runs with and without selfgravity, and a Jupitermass planet at 5.2 au. Moving from left to right: simulations M, M∕2 and M∕4 (i.e. decreasing disc masses). Note that the curves corresponding to different times have been vertically offset for clarity. 

Open with DEXTER  
In the text 
Fig. 4
Azimuthallyaveraged radial surface density profiles at different times for nominal disc mass (M) simulations with a Jupitermass planet at 5.2 au. Left panel: same as the left panel of Fig. 3 for simulation M_{2D}, except that the scale on the yaxis is logarithmic, to appreciate the gap’s depth and its evolution with time. Right panel: same for the M_{3D} simulation. The gaps in 2D and 3D inviscid discs (with adiabatic EoS and a cooling timescale equal to one planet’s orbital period) appear very different; see text for discussion. 

Open with DEXTER  
In the text 
Fig. 5
Left panels: evolution of semimajor axes and eccentricities for planets in the simulations M_{2D} sg, M_{2D} ∕2sg and M_{2D}∕4sg. Right panels: same as left panels except for the corresponding nonselfgravitating disc models. 

Open with DEXTER  
In the text 
Fig. 6
Surface density profile at the times shown in the legend for the run M_{2D}. Also shown are the semimajor axes of the planet (filled circles), and the locations of the 2:1 corotation resonances and of the 3:1 external Lindblad resonances at the times indicated. 

Open with DEXTER  
In the text 
Fig. 7
Evolution of semimajor axis and eccentricity for a migrating planet in simulation M_{3D}. Notice the outward migration after 45 000 yr. The colored dots give the values of (a, e) at times corresponding to the panels of Fig. 8. 

Open with DEXTER  
In the text 
Fig. 8
Contours of surface density for the M_{3D} simulation and a migrating Jupitermass planet. The panels correspond to different times, reported on top of each panel, corresponding to the dots in Fig. 7. It appears clearly that a secondary gap is carved by vortexdisc interaction (second and third panels) and that the vortex weakens by spreading radially, leaving the outer disc partially depleted (fourth and fifth panels). 

Open with DEXTER  
In the text 
Fig. 9
Top panel: evolution of specific torques from the inner (blue curve) and outer disc (red curve) and of the total specific torque (cyan) in simulation M_{3D}. A runningaverage has been applied over 500 orbits to smooth out the oscillations of the torque exerted by the outer disc. Such oscillations are due the synodic period of the vortex with respect to the planet. At times corresponding to the dots in the top panel (the same times selected in Figs. 7 and 8) we plot: (middle panel) surface density profiles; (bottom panel) specific radial torque density γ(r) (see Eq. (2) in the text). 

Open with DEXTER  
In the text 
Fig. 10
Migration speed as a function of time, for discs with different masses in 3D simulations. The planet always starts at 5.2 au. The migration speed decreases, as expected, with the disc mass. The direction of migration is reversed at t ≃ 45 000 yr in the nominal mass disc case (simulation M_{3D}). 

Open with DEXTER  
In the text 
Fig. 11
Migration speed as a function of planet’s distance from the star in 3D simulations of discs with different viscosities (i.e. simulations M_{3D}, M_{3D α5}, M_{3D α4} and M_{3Dα3}; see Table 1). A running average procedure has been used to reduce fluctuations, with the consequence that all but the green curve start at about 5.1 au instead of the released distance of 5.2 au. In the α = 10^{−3} and 10^{−4} cases migration is always inwards (negative radial speed). Notice that we have divided by 10 the migration speed of M_{3D α3} in order to scale it proportionally to viscosity relatively to the α = 10^{−4} case. The migration pattern in the inviscid case does not respect this scaling and turns positive after a short phase of inward migration. The same happens in the M_{3D α5} simulation. 

Open with DEXTER  
In the text 
Fig. 12
Same as Fig. 2, but for a planet located at r_{p} = 13 au. 

Open with DEXTER  
In the text 
Fig. 13
Azimuthallyaveraged radial surface density profiles versus time for 2D runs with and without selfgravity, and a Jupitermass planet at 13 au. Moving from left to right: simulations M∕2, M∕4 and M∕10. Note that the curves corresponding to different times have been vertically offset for clarity. 

Open with DEXTER  
In the text 
Fig. A.1
Left panels: evolution of semimajor axes and eccentricities for planets in simulations M_{2D}, M_{2D} ∕2 and M_{2D}∕4 with cooling time τ_{c} = 0.01 orbit. Right panels: same as left panels except for the corresponding nonselfgravitating disc models. 

Open with DEXTER  
In the text 
Fig. A.2
Contours of surface density for simulation M_{2D} with τ_{c} = 0.01 orbits and a migrating Jupitermass planet. The panels correspond to different times, reported on top of each panel. Top panels: selfgravitating discs. Bottom panels: nonselfgravitating discs. 

Open with DEXTER  
In the text 
Fig. A.3
Contours of surface density for simulations M_{2D} with τ_{c} = 10 orbits and a migrating Jupitermass planet. The panels correspond to different times, reported on top of each panel. Top panels: selfgravitating discs. Bottom panels: nonselfgravitating discs. 

Open with DEXTER  
In the text 
Fig. A.4
Left panels: evolution of semimajor axes and eccentricities for planets in selfgravitating disc simulations M_{2D}, M_{2D} ∕2 and M_{2D}∕4 where τ_{c} = 10 orbit. Right panels: same as left panels except for the corresponding nonselfgravitating disc models. 

Open with DEXTER  
In the text 
Fig. A.5
Contours of surface density for the M_{3D} simulations with the locally isothermal equation of state and for increasing values of τ_{c} (in orbits). All the panels correspond to t =8296 yr (the same time as for the {middle panels} of Fig. 2, for comparison). 

Open with DEXTER  
In the text 
Fig. B.1
Surface density profiles at t = 10000 yr for standard resolution non viscous 3D simulations with different smoothing length and for an additional run with double resolution and smoothing length ϵ = 0.2 Hill radii. We call “SR” the standard resolution simulations and we use the label “HR” for the high resolution one. The "HR" simulation is a restart of the “SR” case with ϵ = 0.4 R_{H} (red curve) at t = 10^{4} yr. The green curve corresponds to 100 orbits after the restart. 

Open with DEXTER  
In the text 
Fig. B.2
Evolution of the planet semimajor axis for the SR M_{3D} simulation with different smoothing lengths and for the additional HR simulation with smoothing length ϵ = 0.2 Hill radii. 

Open with DEXTER  
In the text 
Fig. B.3
Evolution of the planet semimajor axis for the SR M_{3D} simulation with viscosity α = 10^{−5} and for the additional HR simulation with smoothing length ϵ = 0.2 Hill radii. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.