Inward and outward migration of massive planets: moving towards a stalling radius

Recent studies on the planet-dominated regime of Type II migration showed that, contrary to the conventional wisdom, massive planets can migrate outwards. Using `fixed-planet' simulations these studies found a correlation between the sign of the torques acting on the planet and the parameter $K'$ (which describes the depth of the gap carved by the planet in the disc). We perform `live-planet' simulations exploring a range of $K'$ and disc mass values to test and extend these results. The excitation of planet eccentricity in live-planet simulations breaks the direct dependence of migration rate (rate of change of semi-major axis) on the torques imposed, an effect that `fixed-planet' simulations cannot treat. By disentangling the contribution to the torque due to the semi-major axis evolution from that due to the eccentricity evolution, we recover the relation between the magnitude and sign of migration and $K'$ and argue that this relation may be better expressed in terms of the related gap depth parameter $K$. We present a toy model in which the sign of planetary migration changes at a limiting value of $K$, through which we explore planets' migration in viscously evolving discs. The existence of the torque reversal shapes the planetary system's architecture by accumulating planets either at the stalling radius or in a band around it (defined by the interplay between the planet migration and the disc evolution). In either case, planets pile up in the area $1-10$ au, disfavouring hot Jupiter formation through Type II migration in the planet-dominated regime.


INTRODUCTION
Planets form from protoplanetary disc material and continue to interact with this material by exchanging orbital energy and angular momentum via tidal torques. This interaction shapes both planetary orbital architecture and disc structure. Planet-disc interactions have been studied over many decades, long before the detection of the first exoplanet by Mayor & Queloz (1995): see for example Lin & Papaloizou (1979); Goldreich & Tremaine (1979, 1980. Recent reviews of the topic include Papaloizou & Terquem (2006); Kley & Nelson (2012); Baruteau et al. (2014); Papaloizou (2021); Paardekooper et al. (2022).
Once a planet has formed in a protoplanetary disc, depending on its mass it might be able to open a gap in the disc surface density, or it may remain embedded in the disc. On this basis, planet migration can be classified (e.g. Artymowicz 1993;Ward 1998;Kley & Nelson 2012) as type I migration (Ida & Lin 2008;Bitsch et al. 2013), for light planets that migrate embedded in the disc; or Type II migration, for massive planets that open a gap in the disc. In this paper, we are interested in Type II migration (e.g. Lin et al. 1993; Syer & Clarke ★ E-mail: ces204@cam.ac.uk 1995; Ivanov et al. 1999), which can be further split into two regimes, defined through the local disc to planet mass ratio (1) The 'disc-dominated' regime is defined as the regime where the local disc mass is higher than the planet mass; while the 'planet-dominated' regime, conversely, requires the planet to be more massive than the local disc. The classical theory of the disc-dominated regime assumes that the planet is locked in the disc gap, with no possible flux of material crossing its orbit, and it migrates following the disc viscous evolution (e.g. Lin & Papaloizou 1979). In recent years, this simple picture has been questioned by hydrodynamics simulations' results; for example, Duffell et al. (2014) and Dürmann & Kley (2015) suggested that planet migration is unlocked from the disc viscous evolution and the planet can migrate faster than the disc material; other works (e.g. Lubow & D'Angelo 2006) modelled the gas flow through the planetary gap. A number of recent studies investigated the problem (e.g. Robert et al. 2018;Kanagawa et al. 2018;Scardoni et al. 2020;Lega et al. 2021), and Scardoni et al. (2020) proposed that the hydrodynamic simulations mostly undergo a transient phase and over longer viscous timescales the planet migration conforms to the usual Type II picture; however, no definitive conclusion has been reached yet.
Strongly related to the migration problem, is the study of the gap's properties, since most of the angular momentum is exchanged between the planet and the disc in proximity of the gap edges (Goldreich & Tremaine 1979, 1980; some important works analysing the gap shape and width are for example those by Crida et al. (2006), who provided an analytical formula linking the gap's and the planet-disc system properties; Fung et al. (2014), who analysed the expected density in the gap; and by Fung & Chiang (2016), who demonstrated that in 3D simulations, massive planets carve gaps whose properties are consistent with those produced in 2D simulations.
Regardless of the details of migration velocity in this regime, Type II migration is often suggested as an explanation for the existence of massive exoplanets characterised by small semi-major axes (the so-called 'hot Jupiters'). Moreover Type II migration in the planetdominated regime -the main focus of this paper -might help in explaining the population of Jupiter-like planets located at larger radii. Indeed, in the latter regime, the planet's inertia is supposed to have a crucial impact on migration, slowing it down significantly.
The very first papers exploring this problem (Syer & Clarke 1995;Ivanov et al. 1999), studied the problem in 1D under the assumption that no material can cross the planet's location. Since the planet migrates more slowly than the local viscous velocity and the material cannot cross the gap, the fate of the inner disc is to disappear, as it is rapidly accreted by the central star. For the same reason the gas located in the outer disc tends to accumulate at the outer gap edge. As the surface density of the gas at the outer gap edge increases, the torque pushing the planet inwards is enhanced, until reaching an equilibrium, where the planet moves at a rate matching the inward motion of the density maximum.
More recent studies have analysed the planet-dominated regime in different conditions, finding that in some circumstances the planet migration is not only slowed down, but even reversed. Crida & Morbidelli (2007), for example, studied the problem in the presence of a gap with a non-negligible amount of material. Adding the corotational torque in their computations, they found that it can have a significant influence on planet migration: if the Reynolds number is low enough the planet migrates outwards instead of inwards towards the star. Under the assumption that the planet migration is locked to the viscous evolution, another situation for potential outward migration is the case that the planet is initially located in the disc area which is viscously expanding (Veras & Armitage 2004). Hallam & Paardekooper (2018) instead showed that, under certain circumstances, the illumination of the outer gap edge by the central star's radiation can lower the torque exerted by the outer disc on the planet; they thus suggested that this is a possible mechanism to slow down or even reverse massive planet migration.
Further investigation on the planet migration direction has been recently conducted by Dempsey et al. (2020Dempsey et al. ( , 2021, who used a set of 2D systems in the planet-dominated regime, in a stable steady state condition. To study planet migration they model the situation where the planet's orbital parameters (semi-major axis and eccentricity) are not allowed to evolve; they then calculate planet migration from the torques acting on the fixed planet. Since the planet cannot react to the disc either in p or in p , the migration rate for given system parameters is just proportional to the disc mass; migration is parametrised in terms of Δ /( p ), where Δ is the torque acting on the planet, is the steady state accretion rate, p is the specific angular momentum of the planet. Through their study, they found that for typical disc parameters, Jupiter-like (or lighter) planets migrate inwards, while super-Jupiter planets migrate outwards; interestingly, the sign of the torque (and thus direction of migration) seems to be related to the gap parameter = 2 /( ℎ 3 ) (which controls the gap depth, Kanagawa et al. 2016), where = p / * is the planet to star mass ratio, is the Shakura & Sunyaev 1973 viscosity parameter, and ℎ = / is the disc aspect ratio. The parameter is a variant of the gap depth parameter = 2 /( ℎ 5 ), which can be defined from the study of the torques' balance (see for example Fung et al. 2014;Kanagawa et al. 2015Kanagawa et al. , 2018. Being proportional to the planet's gravitational torque (∝ 2 /ℎ 3 ) and inversely proportional to the disc's viscous torque (∝ ∝ ℎ 2 ), it can be interpreted as a measure of the relative strength of the torques.
These interesting results, however, are limited by the assumption of a fixed planet. Albeit in a different migration regime ( > 1, see Equation 1), Scardoni et al. (2020) showed that when the planet is allowed to change its orbital parameters, the disc density readjusts to the presence of the moving planet, which modifies the disc-planet interaction. Although the two migration regimes cannot be directly compared, this work highlights the potential importance of including the evolution of the planet orbital parameters in the computation of disc-planet torques. The importance of the gap-adjustment due to planet migration in determining the evolution of the planet orbital parameters has recently also been underlined also by Lega et al. (2021). Furthermore, in live planet simulations, the torque exerted on the planet by the disc modifies both the semi-major axis and the eccentricity (e.g. Kley & Dirksen 2006;Duffell & Chiang 2015;Teyssandier & Ogilvie 2016;Rosotti et al. 2017;Ragusa et al. 2018), potentially producing a non-trivial relation between the torque and planet migration.
Another fundamental aspect of planet-disc interaction is the long timescale over which the disc re-adjusts itself to the presence of the planet; in fact, although the timescale of tidal interaction is short, the disc structure change might happen over the significantly longer viscous timescale (e.g. Ward 1998;Lin & Papaloizou 1979). The importance of considering the long term evolution of planet-disc systems has been demonstrated by Ragusa et al. (2018), who showed that the long term evolution obtained from hydrodynamic simulations can differ from the trend at the first stages of evolution. Apart from their work -which was mainly focused on the eccentricity evolution -no long term migration studies in the < 1 regime have been conducted with a 'live' (i.e. evolving) planet.
In this paper, we present a suite of 2D simulations lasting 300k orbits, for a variety of aspect ratios / , planet masses and disc masses disc . These simulations are intended to be complementary to those by Dempsey et al. (2021), and allow us to investigate both the importance or otherwise of employing a 'live' planet, as well as the role of different boundary conditions. Since these simulations are extremely expensive from the computational point of view, such simulations have never been performed before; consequently this is a new area of exploration for the disc-planet interaction problem, and the present simulation set inevitably leaves some questions unanswered. Nevertheless, they provide an intriguing starting point for further exploration in the field.
The paper is organised as follows: in section 2 we describe the long numerical simulations performed to study planet-disc interaction; section 3 contains the description of the planets' orbital parameters, and the analysis of the torques resulting from the disc to planet interaction; in section 4 we discuss the influence of boundary conditions on the obtained torques on the planet; a toy model exploring the secular evolution of planets is developed in section 5, and in section 6 we discuss the implications of this model for planetary demographics; finally, in section 7 we draw our conclusions. Table 1. Simulation parameters. The name of each simulation is chosen in the following way: 'L' or 'M' to indicate a light or massive disc, respectively; 'm' followed by a number to indicate the planet mass (measured in Jupiter masses); 'h' followed by a number to indicate the aspect ratio.

Simulation parameters
We performed 12 2D hydrodynamical simulations of protoplanetary discs containing a massive planet with the grid code F 3D (Benítez-Llambay & Masset 2016), considering a cylindrical reference frame ( , ) centered in the star. We adopt dimensionless units = * = 0 = 1, where is the gravitational constant, * is the star mass and 0 is the planet's initial location; the time unit is the inverse Keplerian frequency at 0 , Ω −1 k ; this means that the planet at its initial location 0 requires = 2 to complete an orbit. Each simulation is run for 300k orbits or 600k, depending on the status of the steady state; the requirement is that ( ) is within 10% of ( in ) at least until = 2.5 p .
The setup for the simulations is defined in analogy with the simulations by Ragusa et al. (2018). 1 We consider logarithmically spaced = 430 cells in the radial direction, from in = 0.2 to out = 15, and linearly spaced = 580 cells in the azimuthal direction, extending from 0 to 2 . We assume a locally isothermal equation of state set by Equation 3, and we parametrise the viscosity by applying the prescription by Shakura & Sunyaev (1973), = s , with defined as a function of radius in all the simulations we fix 0 = 0.001. In our simulations we consider a variety of disc aspect ratios and planet masses. Specifically, we consider flared discs, defining the aspect ratio as follows where ℎ 0 is the aspect ratio at the initial planet location, for which we consider 3 different values (ℎ 0 = 0.036, 0.6, 0.1). Regarding the planet mass, we explore 3 different values: p = 1 J , p = 3 J , and p = 13 J . For each combination of parameters ℎ 0 and p , we perform 2 simulations, characterised by different values for the initial local disc to planet mass ratio 0 (see Equation 1). The initial disc surface density at planet location Σ 0 is chosen to obtain the selected values for the initial value 0 ; we thus have 'massive' disc simulations, with 0 = 0.15, and 'light' disc simulations, with 0 = 0.046. Note 1 In turn, their choice of parameters was based on the best fit for CI Tau disc by Rosotti et al. (2017) that in physical units, 0 = 0.046 and 0 = 0.15 for a Jovian mass planet at the radius of Jupiter corresponds to a local disc surface density of 1.3 g/cm 2 and 4.2 g/cm 2 , respectively. Apart from p , ℎ 0 and Σ 0 , all the other parameters are kept fixed among all the simulations. Being interested in the planet-dominated regime of Type II migration, we only considered values 0 < 1; for planet migration with high B values we refer to our previous study (Scardoni et al. 2020), specifically focused on the disc-dominated regime of Type II migration In addition to the simulations outlined in Table 1 we performed two high resolution simulations as convergence test: we performed simulations L-m13-h036 and M-m13-h036 with ×2 and ×4 the standard resolution used in this paper, and verified that the obtained orbital parameters are the same as those obtained using the standard resolution.

Initial and boundary conditions
The initial density profile is defined as The planet mass is gradually increased during the first 50 orbits of the simulation, while its orbital parameters are kept fixed; once the planet reaches its full mass, it is allowed to evolve (i.e. to migrate and develop eccentricity) under the action of torques arising from planet-disc interaction. We apply closed boundary conditions at the outer edge out = 15 by setting both the velocity and the gas density to zero at out , i.e. no material is added to the simulation, in order to minimise the effect of the outer boundary condition. This generates a zero flux at out , which physically corresponds to the case of a disc truncated by a wide binary. At the inner boundary, instead, we apply 'viscous' boundary conditions, i.e. we enforce the inner material velocity to match the viscous velocity (as in Scardoni et al. 2020and Dempsey et al. 2020, 2021. This choice ensures that the inner boundary has a net zero effect on the total angular momentum budget of the disc, i.e. the torque supplied exactly matches the rate of angular momentum advected by accretion through the inner boundary. In practice, the simulations are run for a small fraction of the viscous timescale at the outer edge so that the form of the outer boundary should not be critical to the calculations. On the contrary, boundary effects can disturb the inner disc and thus the region around the planet on much shorter timescales, due to spurious wave propagation. We therefore employ the wave damping method by de Val-Borro et al. (2006) for the radial velocity at the inner boundary. This means that over the damping timescale ( = Ω k /30), we damp to the azimuthal average of the initial (viscous) velocity ,0 in the region from in to = 0.3.
We also perform an additional simulation characterised by the same parameters as simulation L-m13-h036, except for the inner radius and damping zone, which are taken to be half of the default values; thus we take in = 0.1, and the damping zone extends in this case up to = 0.15. This is to ensure that in the cases where the planet develops high eccentricity (and thus gets close to the inner boundary), the material in the damping zone does not affect the migration owing to numerical modification to the physical torque exerted on the planet.

Orbital parameters
In Figure 1 we show for the light disc simulations (left panel) and the massive disc simulations (right panel) the evolution of the semimajor axes as a function of evolutionary time, i.e. the physical time rescaled to the viscous timescale at the initial planet position and we indicate with ,0 the viscous timescale evaluated at the initial planet location. For reference, the initial viscous timescales in our simulations are: ,0 ∼ 9 · 10 5 yr for ℎ = 0.036, ,0 ∼ 3 · 10 5 yr for ℎ = 0.06, ,0 ∼ 10 5 yr for ℎ = 0.1. After an initial adjustment (lasting ,0 ), we can notice a variety of behaviours: inward migration (e.g. L-m1-h06, M-m1-h06), outward migration (e.g. L-m13-h036, M-m3-h036), or even stalling (e.g. L-m13-h1, M-m13-h06). This interesting behaviour can be related to the value of the gap-opening parameter = 2 /( ℎ 3 ) (where = p / * ), whose analysis is deepened in section 3.2. Figure 2, instead, illustrates the planet eccentricity evolution. This plot shows a variety of behaviours, with some simulations (especially at higher planet masses) exhibiting significant eccentricity growth. It is also worth noticing that, in many cases, the planet eccentricity presents oscillations (as already observed, for example, by Duffell & Chiang 2015; Thun et al. 2017;Rosotti et al. 2017;Ragusa et al. 2018). These oscillations are due to the disc eccentricity vector evolving as a superposition of two rigidly precessing normal modes (Teyssandier & Ogilvie 2016Ragusa et al. 2018;Teyssandier & Lai 2019).

Dependence on initial disc mass
In this section, we analyse the torque acting on the planet as a function of the following parameter (Kanagawa et al. 2015 which is a modification of the gap-opening parameter = 2 /( ℎ 5 ), i.e. a measure of the strength of the planet's gravitational torque (∝ 2 /ℎ 3 ) compared to the disc's viscous torque (∝ ∝ ℎ 2 ).
The parameter has been empirically studied by Kanagawa et al. (2016Kanagawa et al. ( , 2018 and by Dempsey et al. (2020Dempsey et al. ( , 2021, who noticed that correlates with the gap width and depth (with the gap becoming wider and deeper for higher values). Systems characterised by the same also appear to show similar torque density profiles; thus Dempsey et al. (2020Dempsey et al. ( , 2021) examined a number of fixed-planet 2 simulations characterised by different values, and suggested that is a good 'ordering parameter' for planet migration. According to their analysis, in fact, planets in systems characterised by 20 migrate outward, whereas systems with 20 present inward planet migration. Since their finding relies on simulations characterised by an evolving disc and fixed planets, our goal is to test whether the torque imposed on a fixed planet on a circular orbit is the same that applies in the case of a 'live' planet which can respond to the torque during the evolution by changing its orbital elements.
For this purpose, we follow the approach of Dempsey et al. (2021), and plot the total torque on the planet normalised to p (see their Fig. 2). Note that for our convention of signs we have < 0 for inward migration; to avoid confusion we therefore plot Δ /| p | whose sign only depends on the sign of the torque exerted by the disc on the planet (positive for outward migration and negative for inward migration).
The presence of a live planet, enables us to compute the torque acting on the planet directly as the planet's angular momentum variation rate being p = p √︃ * p (1 − 2 p ). Thus, we compute Δ timeaveraged over 30000 planet orbits, 3 taken in the last stages of the simulation, when the inner disc has reached a steady state and the planet migration has stabilised. The results are shown in Figure 3, from which we can notice that the general trend of positive torques for high , and negative torques for low is confirmed, albeit with a suggestion of some decline in the torque at values of > 20−200. Note, however, that the suite of simulations performed so far does not explore the values extensively enough to confirm that = 20 is the point of zero torque, and further simulations are needed for a precise characterisation of the Δ − relation. We further notice an applicability limitation in the results obtained from fixed planet simulations by Dempsey et al. (2020Dempsey et al. ( , 2021. Those simulations do not conserve the total angular momentum of the system, because they do not allow the planet to evolve its orbital parameters. Although this approach is completely adequate to represent the limiting behaviour in the case that 0 → 0 (for which the planet is not expected to modify its orbital parameters noticeably over the disc lifetime), it cannot necessarily be used for 0 ≠ 0 studies, where the planet might change its orbital parameters over the simulation timescale, and in the process affects the interaction with the disc. This limitation becomes evident when we compare our results for different 0 values (i.e. 'light disc' vs 'massive disc'), where we notice that the normalised torque properties change for different 0 values in some cases. In the approach adopted by Dempsey et al. (2021), by contrast, both the torque and the accretion rate are assumed to vary linearly with 0 and hence their ratio is independent of 0 ; therefore in Dempsey et al. (2020Dempsey et al. ( , 2021 simulations  Table 1 for the simulations' name definition).

Figure 2.
Eccentricity as a function of the evolutionary time for light disc simulations (left panel) and massive disc simulations (right panel). The different colours correspond to different simulations, as indicated in the plot legend (see Table 1 for the simulations' name definition).
the torque applied to the planet is, by construction, proportional to the disc mass. Our simulations demonstrate instead that although the assumption of torques independent of the disc mass is adequate for low mass planets, it breaks down in the case of high (> 10 3 ).
We find that the normalised torque values are indeed approximately independent of 0 in the case of the Jovian mass planet (see Figure 3) but that for higher values of planet mass (and ) the normalised torque values become increasingly divergent between the 'light' and 'massive' simulations. This is particularly marked in the case of the simulation with the largest value (the 13 J planet with the lowest disc aspect ratio). These results can be readily understood from Figure 2 where it can be seen that the planet develops significant eccentricity during the course of the simulation. As discussed by Ragusa et al. (2018); Teyssandier & Lai (2019), the eccentricity evolution of the planet depends on the interplay between eccentric modes of the disc whose structure depends on the disc to planet mass ratio. Thus we see in Figure 2 that not only do the simulation with largest develop large eccentricities but that the eccentricity evolution is markedly different in the 'light' and 'massive' cases.
It is therefore unsurprising that the torques on the planet evolve differently; indeed, in the light case, the eccentricity grows to the point where although the semi-major axis increases, the torque on the planet is actually negative (Figure 3), a result that is reconciled by the large planetary eccentricity in this case. We could thus consider simulation L-m13-h036 as an 'outlier' in Figure 3 and Figure 7 because it is the extreme case where the torque mainly affects p rather than p .
The different evolutionary histories of the disc structure in the light and massive 13 j case are illustrated in Figure 4 (in the left and right panel, respectively). The strong growth of planetary eccentricity in the light case is associated with the system settling into a single eigenmode of the system, i.e. the slow mode in which the apsidal precession of the disc and planet are aligned. In the massive case, conversely, the system exists in a state of superposition of aligned and anti-aligned modes (see Ragusa et al. 2018, for more discussion on  . Azimuthal averaged density profiles Σ/Σ 0 (colour map) as a function of the radius (x axis) and time normalised to the viscous timescale at the initial planet location (y axis). The magenta line shows the time evolution of the planet's semi-major axis; the white lines show the time evolution of p (1 ± p ). The left panel refers to simulation L-m13-h036; the right panel refers to simulation M-m13-h036. the two modes), resulting in the modulation of planetary eccentricity and disc structure evident in Figures 2 and Figure 4. 4 We therefore conclude that the results based on non-migrating planet simulations presented by Dempsey et al. (2020Dempsey et al. ( , 2021) are a good approximation for light planets; nonetheless a more complex model, accounting for the planet's eccentricity growth and its dependence on disc mass, is needed when higher mass planets are considered.
Finally we emphasise that Figure 3 is constructed after significant evolution of the disc profiles (i.e. after 6 · 10 5 planetary orbits in the case of the massive and light simulations with 1 j planet, and in the massive disc simulation with 3 j planet; and after 3 · 10 5 planetary orbits in the remaining simulations) and use values of the torque and accretion averaged over the preceding 3 · 10 4 orbits. The result that the simulations with higher values generally result in normalised torque values that are ∼ 0.5 is all the more remarkable given the substantial variety in the surface density evolution that occurs in the various simulations. This is illustrated in Figure 5 which shows snapshots of the surface density profiles (as a function of radius normalised to the current semi-major axis of the planet) in each of the simulations at the time that the normalised torques shown in Figure 3 are evaluated. In each case the solid and dashed versions of the curves relate to corresponding pairs of simulations with different values of 0 . The density is normalised to the initial surface density at the location of the planet in the unperturbed disc and therefore the relative surface density in each simulation pair can be obtained by scaling with the relevant values of 0 for the light and massive disc simulations. Figure 5 demonstrates that it is only in the case of the simulation pair (L-m1-h06 and M-m1-h06) with lowest (i.e. narrowest gap) that the surface density profiles are simply scaled versions of each other. This can be traced to the fact that the planetary eccentricity is not excited in this case (see Figure 2) as can be expected in the case of a narrow gap where corotation torques dominate over Lindblad resonances (Goldreich & Tremaine 1979, 1980. At the other extreme, the simulations where the planetary eccentricity is most excited (i.e. the m13-h036 and m13-h1 pairs) show that the inner edge of the planet carved cavity is close to the 3:1 inner eccentric Lindblad resonance which is associated with strong driving of orbital eccentricity (Lin et al. 1993;Goldreich & Sari 2003).
We note that, in contrast to the simulations of Dempsey et al. (2020Dempsey et al. ( , 2021, we do not drive the accretion rate at a prescribed rate in these simulations; the accretion rate at the inner edge (which is applied to the calculation of the normalised torque in Figure 3) evolves rapidly during the first viscous timescale at the planetary radius, thereafter setting up a quasi-steady state within a few times the planetary radius (e.g. see Figure 10). Thereafter, the magnitude of the accretion rate declines on the viscous timescale at the half mass radius of the disc, as expected. It is this slow decline in the accretion rate in the inner disc (and the associated normalisation of the surface density) that drives the leveling off in the evolution of the semi-major axis seen in Figure 1. In section 5 we will explore toy models which use the normalised torque values extracted from our simulations to calculate the secular evolution of planets in an evolving disc.

Torque components
So far we have evaluated the total torque on the planet, including contributions from both rate of change of semi-major axis and rate of change of eccentricity, as in Equation 7. In Figure 6 we plot instead the torque component associated with the rate of change of only the semi-major axis against the parameter We notice that in this case there is little difference between the torques associated to p for the light and massive disc simulations with 13 j and ℎ = 0.036. This suggests that the differences in total torque between these two simulations do not result from interactions that transfer significant orbital energy to the planet, but are instead associated with interactions that primarily affect the planetary eccentricity.
We furthermore notice that at large (or , see the next section), the normalised Δ * is almost constant with (or ), whereas Δ seems to decrease. This behaviour suggests that as the planet mass increases, the resonances that excite/damp the planet's eccentricity are strengthened/ weakened; while the resonances responsible for migration seem less affected.
This result has interesting implications on the rate of energy transfer to the planet. Given that in this case the variation of normalised Δ * between simulations is much less than if the total normalised torque is considered, even the L-m13-h036 simulation (where the total normalised torque is large and negative: see Figure 2) has a positive normalised Δ * value whose magnitude (∼ 0.5) is similar to that in lower mass planets. This implies that the spread in normalised Δ values for the most massive planets is driven by differences in eccentricity evolution which are not associated with significant contributions to the transfer of energy between disc and planet. Given that the rate of change of planetary orbital energy is given by Combining Equation 10 and the torque associated with the semimajor axis variation |Δ * | (Equation 8) we deduce that Since from our simulations the normalised Δ * is roughly constant (right panel in Figure 7), Equation 11 implies that the energy transfer to the planet at a given location depends only on the accretion rate through the disc and is independent of planet mass. Moreover we note that the migration rate, p , scales with / p and, for simulations with the same value of 0 and , should be independent of planet mass. This can be seen from the asymptotic gradients of the migration tracks shown in Figure 1.
The analysis of the torque components allowed us to attribute the deviation of the torques Δ from the constant value at high (or ) to the planet eccentricity excitation, and we noticed that the effect becomes more and more important as ( ) increases. We plan to further explore the behaviour at high (or ) in future works.

Choice of ordering parameter
The torque analysis conducted so far is following Dempsey et al. (2020Dempsey et al. ( , 2021 in assuming that is the appropriate ordering parameter. In Figure 7 we instead plot the normalised torques as a function of for all our simulations, considering both the total torque (left panel) and the torque associated with p (right panel). While the number of simulations is too small to make a definitive judgement on the parameter that best captures the transition between narrow gaps with inward migration and broad gaps with outward migration, Figure 7 suggests that this transition occurs at = 1.5 · 10 4 . We use this threshold in our 'toy modeling' presented in section 5. Remarkably the normalised torque associated with assumes a constant value of ∼ 0.5 for higher values of .

ROLE OF INNER BOUNDARY CONDITIONS
In this section, we discuss the role of boundary conditions in planet migration simulations. For this purpose, we performed a set of additional simulations, characterised by the same parameters as those presented above, but with 'open' boundary conditions at the inner boundary. The open boundary conditions allow the material to leave the grid at the inner edge at its own radial velocity (numerically, we set both the velocity and the density in the inner ghost cell equal to the last active cell). As a consequence, at in = 0.2 we have a zero torque boundary condition, as illustrated in Figure 8, where the viscous torque is shown as a function of at different snapshots for simulation M-m3-h036 with open boundary conditions. 5 From the physical point of view, the open boundary conditions produce a strong depletion of the inner disc; we can therefore apply them to study the physical case of a disc with a zero-torque inner cavity; for example, it may represent a case where the angular velocity is subject to a turning point (e.g. the classical boundary layer), or any case where material is removed from the cavity without injecting angular momentum to exterior material (e.g. photoevaporation or magnetospheric accretion).
The different physical configuration naturally produces differences in the migration properties which can be analysed in terms of Δ in /| p |, where Δ in is the torque imparted to the planet from 5 In this section, we use simulation M-m3-h036 as an example to illustrate the viscous torque and the accretion rate profile through the disc. Note, however, that the same behaviour is shown by all the other simulations in our sample. the disc interior to the planet. In order to analyse the different torques acting on the planet in the two cases, we commence by noticing that for both the boundary condition choices, the inner disc (from in to a few p ) reaches a viscous steady state after some ,0 (see Figure 10). If we consider the area of the disc between in and p , the total angular momentum injected per unit time into this region is due to advection and viscous torques at in and p is: where the approximation is valid for planets massive enough to create a deep gap in the disc, so that Σ( p ) ∼ 0, and thus whereas for viscous boundary conditions ( in ) = − in , thus We can therefore estimate the difference in torque imparted by the inner disc where < 0 (for our convention of signs), thus Δ viscous BC > Δ open BC . Since in both cases, the quasi-steady state of the inner disc requires that Δ is balanced by the torque applied by the planet, it follows that the form of the inner boundary condition should alone determine the normalised torque of the inner disc on the planet. The overall physical consequence is that planets in discs with an inner cavity are expected to be less prone to outward migration than planets in discs where the angular momentum lost by advection is balanced by the viscous torque of inward lying material.
Using Equation 16 we can predict the open boundary conditions' results from the values obtained for viscous boundary conditions (and viceversa). We illustrate the extent to which this simple prediction of how the inner boundary condition affects the planet migration through the plot in Figure 9. In that plot, we show the torque acting on the planet (computed as time derivative of the planet angular momentum) for both the viscous boundary conditions (squares) and the open boundary conditions (diamonds). The stars, instead, show the prediction for the open boundary conditions' results, based on the viscous boundary conditions' results, using Equation 16, and assuming that the inner boundary condition does not affect the torque on the planet from the outer disc. Each colour corresponds to a different simulation, which is indicated in the legend. Figure 9 shows that for the 1 Jupiter mass simulations with ℎ = 0.036, the difference in normalised torque is almost exactly that predicted by Equation 16 (with in = 0.2). Indeed, the outward migration is significantly slowed when a zero torque boundary condition is applied at in = 0.2, because of the diminished torque delivered by the depleted inner disc. The prediction of Equation 16 is appropriate as long as (a) the inner disc is in steady state interior to the planet (as assumed in the derivation of Equation 16), and (b) the modification of the structure of the inner disc due to the changed boundary condition does not modify the torque exerted by the outer disc. This explains the quantitative agreement between the predictions of Equation 16 and the change in total torque experienced by the planet shown in Figure 9 for simulations characterised by low mass planets, where both the conditions are satisfied. We have verified, through examination of radial profiles of accretion rate (see Figure 10), that the condition (a) is satisfied in all simulations. We therefore attribute the fact that the quantitative agreement with Equation 16 declines as the planet mass increases to the breakdown of assumption (b). This is particularly evident in the case of the light 13 j simulation with ℎ = 0.036 where Equation 16 fails to replicate the dependence of torque on boundary condition both in magnitude and in sign. It is notable that these parameters result in significant planetary eccentricity in the case of both boundary conditions but that the eccentricity growth is weaker in the case of the zero torque (cavity) boundary condition (compare Figure 2 with Figure 6 of Ragusa et al. 2018). The modification of the outer disc structure in response to the different orbital evolution of the planet means that the difference in total torque is not simply related to the effect of the boundary condition on the inner disc alone.

EXPLORATION OF SECULAR EVOLUTION USING A TOY MODEL
In this section we define a toy model which provides a simplified prediction of massive planets' migration, under the following key assumption: (i) Δ / p is a function of (Equation 12); (ii) Δ / p is independent of disc mass. 6 Since the disc is flared, we expect ∝ ℎ −5 to vary with radius as ∝ −5 , where is the flaring index; in this section we assume constant , but if we allowed it to vary with radius, it would have given a further contribution to the dependence of on .
We compute the planet migration time using the angular momentum definition p p = p √︁ * p 7 where mig < 0 (> 0) means inward (outward) migration. Using the definitions d = 4 2 Σ and = 3 Σ, we can rewrite the 6 As underlined in section 3.2.1, Δ / p does not depend on the disc mass only when the planetary orbit maintains a low eccentricity; this is best satisfied for the 1m J case, which we therefore use in the toy model. 7 Here we use p instead of p because we are assuming circular orbits.
We can combine the two limits (a) and (b) as Δ /| p | = −0.5/( + 1); then we put all of these considerations together to obtain where lim is the limiting value which determines the transition from inward to outward migration. For the purpose of the toy model, we take lim = 1.5 · 10 4 (as estimated from Figure 7). Given all these considerations, we can define a function ( ) such that 9 We therefore need a function with the following properties: for high and low , ( ) ∼ 1, to obtain the lower case in Equation 20; for low we want ( ) ∼ −1, to recover the classical prediction for Type II migration (upper case in Equation 20). Note that in deriving the formula in Equation 21 we have not considered the case with high and high , because for the typical disc parameters and plausible planet location, it is unlikely to obtain this combination; consequently, for high values of we expect the planet to migrate inward following the 'classical' Type II migration rate (as analysed in Scardoni et al. 2020). We thus model the functional form of ( ) as where controls the width of the transition in over which ( ) changes sign. We then estimate the planet migration by solving the following equation where we underline that all the quantities on the right hand side are functions of p (and hence time).
Since planets are expected to migrate inward (outward) for values smaller (bigger) than lim = 1.5·10 4 , and because increases during inward migration and decreases during outward migration in the case of flaring discs ( > 0 in Equation 24 below), we expect them to always migrate towards the location of the disc where = lim . We then call lim the radius corresponding to lim = 1.5 · 10 4 which is therefore an 'attractor' for migrating planets where = 0 and ℎ = ℎ 0 . Note that in the toy model we do not include the effect of eccentricity growth, because we develop the toy model in the case of Jovian mass planets where eccentricity growth is very modest. We neither include the effect of the planet on the surface density profile, because by comparing the evolution of the disc with and without a Jupiter mass planet we found that the two systems behave similarly over the long timescales.

Toy model in a non-evolving disc
We first define a disc model characterised by = 0.001 (constant throughout the disc), thus ∝ 2 +1/2 ; and Σ( ) = Σ( 1 )( / 1 ) −1 , where 1 is the value of 1 au in code units, and Σ( 1 ) is chosen to have 5au = 0.1 at 5 au, 10 corresponding to disc ∼ 0.001 M for a disc of 100 au. As the planet migrates (in either direction), the value of varies as 11 To design a disc model with lim = 5 au (i.e. the location of Jupiter), the disc aspect ratio at 1 au is taken equal to 0.02, while the flaring index is 0.25. We then insert a Jupiter mass planet and we solve Equation 23. In Figure 11 we illustrate the behaviour of lim as an attractor, considering 3 different models: a planet initially located at p ( = 0) = 3 au < r lim , which then migrates outwards until stalling at 5 au (red lines); for p ( = 0) = 5 au = r lim the planet 10 In the rest of the paper we will refer to the value of at 5 au as 5au ; while we will indicate as the value at the planet location. 11 Note that the dependence of on p depends on the chosen density profile: if Σ ∝ − , then ∝ 2− p . Figure 11. Toy model for a planet of mass p = 1 j and a disc characterised by = 0.001, ℎ = 0.02 at 1 au and flaring index 0.25, producing a stalling radius located at 5 au. The different colours show the planet migration tracks for an initial planet location equal to 3 au (red lines), 5 au (cyan lines), and 7 au (magenta lines). The different line styles refer to different choices of parameter : = 1 for the solid line; = 1500 for the dashed line; = 3000 for the dotted line. At 5 au, the disc to planet mass ratio is = 0.1, while the viscous timescale is (5au) = 1.6 Myr stalls at its initial location (cyan line); for p ( = 0) = 7 au > r lim , instead, the planet migrate inwards (magenta lines). The solid, dashed and dotted lines refer to values of parameter in Equation 22 equal to 1, 1500, 3000, respectively.
For the inward migrating planet (initially at 7 au) the timescale required to reach the stalling radius in all cases is less than 10 Myr, which means that the planet can reach the stalling location during the disc lifetime. In the case of the planet initially located at 3 au, instead, the planet migrates to 4 au in 10 Myr, but requires ∼ 15 Myr to reach the stalling radius. We caution, however, that the timescale on which the planet makes its final approach to the stalling radius is very sensitive to the choice of in the function ( ); this underlines the importance of further simulations to explore the form of ( ). We furthermore notice that as the obtained migration timescale is comparable to the disc lifetime, we would expect the disc density to evolve in time over the planet's migration and model this possibility in the following section.

Toy model in an evolving disc
Since the planet migration timescale expected from the toy model is comparable to the disc lifetime, we define a variant of the toy model to take into account of the density evolution during migration.
For this purpose, while solving Equation 23 we evolve the density profile at each timestep according to the self-similar solution by Lynden-Bell & Pringle (1974) where disc,0 is the initial disc mass, defined by the choice of the density profile and the initial disc scale radius; s,0 is the initial scale radius; is a parameter of the model which is the power law exponent for the radial variation of the viscosity that we fix to 1 for We show in Figure 12 the evolution tracks for a Jupiter mass planet initially located at 3 au (left panel) and 7 au (right panel), in a disc whose density is evolving as described above. In both cases, we consider different models, characterised by a range of initial scale radii (50 au with the solid lines, 100 au with the dashed lines, and 200 au with the dotted lines), and a range of 0 (from 0.09 to 0.5, see the plot legend). For reference, the initial viscous timescale at the planet location are (3 au) ∼ 1 Myr, and (7 au) ∼ 2.5 Myr; while the viscous timescales at the chosen scale-radii are (50 au) = 4.4 Myr, (100 au) = 8.9 Myr, (200 au) = 17.7 Myr. Since the density is reduced while the disc evolves, the migration timescale becomes longer in both the cases when we consider an evolving disc (compare the red lines to Figure 11). In fact, a decrease in the surface density causes the instantaneous value of the parameter to decrease, making the migration timescale longer. This means that models characterised by faster density evolution (i.e. those with lower s,0 ) slow down the migration more effectively than those characterised by slower density evolution; this behaviour can be easily seen in Figure 12 by comparing the dotted, dashed, and solid lines with the same colour (i.e. with different s,0 but same 5au ).
Even more important is the effect of considering different 5au values, with the migration timescale increasing while 5au decreases. Furthermore we notice that in all the models the planet starting from 3 au is more affected by the surface density time evolution, because the value of at its initial location is lower, thus it is more sensitive to further reductions in the local density value. In contrast, the planet starting from 7 au reaches the stalling radius before decreasing significantly the value of , so that it can reach the stalling radius in 10 Myr even in the model characterised by the lowest 5au and the smallest s,0 (see the red solid line in the right panel).
This result has interesting implication for the properties of a system characterised by outward migrating planets. From Figure 12, we can deduce that to have a Jupiter mass planet migrating to 5 au in our model, we need to take either 5au 0.4, or we can decrease it to 5au 0.2 if we take a disc with s,0 100 au. For lower 5au values, even in the case of discs with large scale radii, the planet fails to migrate to the stalling radius in 10 Myr; nonetheless, we still expect some outward migration for those planets, with a final radius which is determined by the disc lifetime rather than by the value of the stalling radius.

Stalling radius and disc properties
In this section we focus on how the disc properties affect the location of the stalling radius. To investigate this problem we consider a Jupiter mass planet, and we compute the disc's properties required to have the planet stalling at a given radius 0 . In the left panel of Figure 13, we consider 3 different stalling radii -3 au (magenta line), 5 au (blue line), 10 au (cyan line) -and for each we plot the disc temperature at 1 au as a function of the parameter. Focusing first on the blue line, we notice that for a Jupiter mass planet that stalls 5 au the disc temperature at 1 au may vary from 10 K to 10 2 K, when we choose values for the parameter in the range ∼ 10 −4 − 10 −2 . 12 Since 10 − 10 2 K at 1 au is a sensible temperature range for protoplanetary discs, this model is in agreement with the findings by Fernandes et al. (2019) and Nielsen et al. (2019) who suggest that around solar-like stars there is a peak in the number of giant planets located at ∼ 5 au from the star. We further notice that for higher disc temperature the planet stalling radius decreases for fixed because Figure 13. Disc temperature at 1 au as a function of the viscosity parameter , required to have a Jupiter mass planet stalling at 0 ; the magenta line correspond to 0 = 3 au, the blue line correspond to 0 = 5 au, the cyan line correspond to 0 = 10 au.
for a hotter disc, a larger region of the disc is in the low (shallow gap) regime where inward migration is expected.
If we consider more massive stars, the planet to star mass ratio decreases as ∝ −1 * , whereas the disc aspect ratio ℎ 0 ∝ √︁ / * , with ∝ * (the relation is sub-linear, and the exact value of depends on the assumed stellar mass-luminosity relation). Using these scale relations in Equation 24 we find therefore we expect the stalling radius value to increase (decrease) with increasing stellar mass for < 1/5 ( > 1/5). At fixed radius, the temperature approximately scales as ∝ 0.15 * (Sinclair et al. 2020), we therefore expect that for the typical disc the location of the stalling radius increases for higher stellar mass; this is in agreement with the results from Nielsen et al. (2019) that higher mass stars should have the peak of the Jupiter mass planet distribution at higher radius with respect to solar mass stars.

IMPLICATIONS
Planetary population synthesis models for giant planets are widely based on the theory of Type II migration (Ida & Lin 2004, 2008Mordasini et al. 2009a,b;Ida et al. 2018;Bitsch et al. 2019). In its classic form this involves inward migration on a viscous timescale which decelerates at the point that the planet mass becomes comparable with the local disc mass. However, since the viscous timescale decreases with decreasing orbital radius, the net effect is that the planet arrives at the disc inner edge in a finite time (Syer & Clarke 1995;Ivanov et al. 1999). The continued inward driving by the disc is related to the assumption of zero flow past the planet so that material can always accumulate exterior to the planetary orbit and drive continued inward migration.
The present paper has demonstrated that even in calculations where the planetary orbital elements are allowed to evolve, giant planets do not arrive at the disc inner edge in a finite time. This is because as the planet migrates in to regions of the disc where the local disc mass dominates the planet mass, the disc establishes a quasi-steady state flow past the planet (cf Dempsey et al. 2020Dempsey et al. , 2021, preventing the inexorable extraction of angular momentum of the planet by material accumulating outside the planetary orbit. Moreover, as previously analysed by Dempsey et al. (2020Dempsey et al. ( , 2021 the sign of the torque depends on the quasi-steady state structure of the gas in the vicinity of the planet. A key point is that the disc interior to the planet is always in a steady state and so, in order to ensure net zero accumulation of angular momentum in this region, must impart a spin up torque to the planet of magnitude p . On the other hand, the disc outside the planet is never globally in a steady state and thus the angular momentum that it extracts from the planet is not constrained to be p , its magnitude depending on the structure of the gap exterior to the planet. Thus for wider gaps (either due to higher planet mass or lower disc aspect ratio ℎ), the spin down torque on the planet is reduced as a result of the substantial clearing external to the planetary orbit: Figure 5 demonstrates that the asymmetry between surface density interior and exterior to the planet becomes more pronounced for higher planet masses and colder (geometrically thinner) discs. We show that even in cases where the planet acquires a significant eccentricity as a result of interaction in the disc (i.e. those cases where the inner edge of the gap extends to close to the 3:1 mean motion resonance), the evolution of the planetary semi-major axis is well described in terms of a simple switch at a value of K (Equation 22) of around 1.5 · 10 4 (Figure 7).
We illustrate in section 5 the implementation of some toy models that incorporate this evolutionary phenomenology. We examined the possibility that the location of zero torque might be imprinted on the architecture of planetary systems. As expected, the planetary evolutionary tracks demonstrate a convergence towards this attractor (Figure 11), located at 5 au given the planet mass and normalisation of the disc temperature profile (see Figure 13 for the sensitivity of the attractor location to disc parameters for a Jovian mass planet). However, the evolution may be too slow for planets to necessarily get close to the attractor location. First of all, the attainment of high values (associated with planet stalling) is favoured by relatively low alpha values where the viscous timescale is relatively long. Secondly, the speed of convergence upon the attractor location depends on the range of disc radii over which the torque magnitude undergoes a sign switch. Clearly, if this occurs relatively gradually, there will be a broad radial range where the torque values are low and therefore the evolution towards the attractor becomes very slow. Our present simulation set (right hand panel of Figure 7) does not allow us to tightly constrain the behaviour near the point of zero torque. Finally, because the planet is not damming up the disc upstream of the planetary orbit, the disc's effect on the planet is progressively weakening on a timescale set by the viscous timescale of the outer disc. Thus whether the planet can be driven to the attractor depends on the migration timescale ( 0 / ) compared with the viscous timescale of the outer disc: as illustrated in Figure 12, efficient driving of the planet to its attractor location during the observed lifetime of protoplanetary discs (e.g. Alcalá et al. 2014;Manara et al. 2016) is achieved for relatively high values (corresponding to fast migration timescale), and high s,0 values (corresponding to slow viscous evolution). Given these considerations, which depend on poorly constrained parameters such as the value of and the form of the torque's dependence on gap shape in the region close to the location of zero net torque, it would be premature to argue that this process should impose a strong signature by piling up planets at a specific orbital location. Nevertheless, it is expected that Jovian mass planets should accumulate in the 1 − 10 au range by this process.
What is certainly clear, however, is that these results preclude the production of hot Jupiters by disc mediated migration, at least in a non-self gravitating disc. 13 This is because it is implausible that the disc gas surface density in the inner disc would ever be high enough to maintain > 1 (where roughly viscous timescale inward migration is expected) right down to the inner edge of the disc. For example, for a Jupiter mass planet at 0.1 au, = 1 corresponds to such a high surface density that the total disc mass would likely exceed a solar mass within a few au. If the candidate hot Jupiter in the Classical T Tauri star CI Tau (Johns-Krull et al. 2016) is confirmed (see counterarguments by Donati et al. 2020) it would then present a puzzle concerning how such a massive planet would have arrived in the innermost disc by an age of a few Myr. Rapid inward migration of giant planets during the earliest self-gravitating disc phase has been proposed by Baruteau et al. (2011), Malik et al. (2015, though more recent works suggest that even then it may not be possible to migrate into the innermost disc (Stamatellos & Inutsuka 2018;Rowther & Meru 2020).
It is also interesting to notice that the existence of planet traps have already been suggested in a series of papers by Hasegawa & Pudritz (2013, for lower mass planets (or planetary cores). They showed that low mass planets can be trapped at zero torque locations, potentially produced at disc radii characterised by significant density/thermal gradients, such as dead zone boundaries and ice lines.
We finally caution that this work relies on the assumption that only one planet is formed in the disc. In the case that multiple planets are formed, simulations including more migrating planets are needed.

CONCLUSIONS
In this work we have presented a suite of long term (300 − 600k planet orbits), 2D hydrodynamical simulations to test and expand recent findings by Dempsey et al. (2020Dempsey et al. ( , 2021 of an empirical correlation between the direction of planet migration and the value of the modified gap-opening parameter . We extend the torque analysis to systems with finite disc mass by considering a live planet allowed to modify its orbital parameters. Our live planet simulations confirm that there is a switch between inward and outward migration which is associated with the creation of deep gaps in the disc. Gap depth is an increasing function of planet mass and a decreasing function of disc aspect ratio and can alternatively be parametrised by the quantities = 2 /( ℎ 3 ) (Equation 6) and = 2 /( ℎ 5 ) (Equation 12). Our simulations suggest that may be the better ordering parameter for describing planetary migration (see right hand panel of Figure 6) and that the switch from inward to outward migration occurs at lim = 1.5 · 10 4 .
We further notice that, regardless the choice of the ordering parameter, as the planet mass increases some dependence on 0 of the torques acting on the planet is revealed; this effect could not be seen in fixed planet simulations, where the normalised torque Δ /| p | values are disc mass independent by construction. This means that if we consider low mass planets, the fixed-planet simulations by Dempsey et al. (2020Dempsey et al. ( , 2021 are a good approximation of the planets' migration; if we consider higher mass planets, the results are modified by the growth of the planet's eccentricity in a way that depends on the disc mass (Ragusa et al. 2018;Teyssandier & Lai 2019). Nonetheless 13 Note that the contrary conclusion by Rosotti et al. (2017) was a result of these authors not adopting viscous boundary conditions so that the loss of the inner disc led to a net negative torque on the planet: the Rosotti et al. (2017) simulation corresponds to the diamond for the L-m13-h036 simulation in Figure 9 we showed that, by disentangling the contribution to the torque due to the semi-major axis variation from the contribution due to the eccentricity evolution, the massive planet migration is well described by the change of sign of the disentangled torque at lim ∼ 1.5 · 10 4 . It is a general feature of protoplanetary discs with realistic heating that the aspect ratio of the disc increases with radius and thus that is a decreasing function of radius. We thus predict inward (outward) migration for radii exterior (interior) to the location where has its limiting value.
We then model the migration behaviour of massive planets by describing the dependence of migration on the parameter . We evaluate this migration in the context of the secular evolution of a viscous disc. This model allows us to obtain the following results: (1) since planets migrate inwards (outwards) for < lim ( > lim ), they tend to go towards the location in the disc where = lim . We thus propose the existence of a 'stalling radius' defined by the location where = lim (i.e. the location of zero torque on the planet); (2) we study the dependence of the stalling radius on the disc parameters (temperature and viscosity parameter), finding that typical disc parameters enable stalling radii in the range 3-10 au, in agreement with the peak in the Jupiter distribution at a few au (Fernandes et al. 2019;Nielsen et al. 2019); (3) when we include the effect of disc density evolution in the model, the migration is slowed down, due to the density reduction with time (which reduces the parameter ). As a consequence, the planet migration towards the stalling radius might be limited by the disc lifetime in rapidly evolving systems characterised by relatively low 0 .
The toy model suggests that while planets should migrate towards a stalling radius set by the planet mass and disc aspect ratio, whether or not they attain their stalling radii depends both on the initial location of the planet and the over-all radial extent of the disc, since the latter determines the rate at which the disc surface density declines. Thus it is likely that this effect does not imprint a strong pile up of planets at their respective stalling radii but rather causes them to occupy a broad band of radii in the range 3 − 10 au. Future quantification of the torque dependence on in the vicinity of lim will help to constrain this further. In any case we do not expect planets to be able to move in from this band by disc mediated migration, thus posing a difficulty for hot Jupiter formation at early times.
(DUSTBUSTERS). This work was in part performed using the Cambridge Service for Data Driven Discovery (CSD3), part of which is operated by the University of Cambridge Research Computing on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC component of CSD3 was funded by BEIS capital funding via STFC capital grants ST/P002307/1 and ST/R002452/1 and STFC operations grant ST/R00689X/1. This work was in part performed using the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure.

DATA AVAILABILITY
The code used to perform the simulations contained in this paper (F 3D, Benítez-Llambay & Masset 2016) is available at http: //fargo.in2p3.fr. In Figure A1 we show the semi-major axis (left panel) and eccentricity (right panel) as a function of orbit, for both the original L-m13-h036 run (magenta line) and the simulation test with in = 0.1 (black dashed line). We notice that the semi-major and eccentricity evolution is essentially the same in both the runs, confirming that the planet evolution in simulation L-m13-h036 is not affected by boundary numerical effects. This paper has been typeset from a T E X/L A T E X file prepared by the author.