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

Chiara E. Scardoni,1 Cathie J. Clarke,1 Giovanni P. Rosotti,2,3 Richard A.
Booth,4Richard D. Alexander,2 and Enrico Ragusa.2,5
1Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge
CB3 0HA, UK
2School of Physics and Astronomy, University of Leicester, University Road,
Leicester LE1 7RH, UK
3Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The
Netherlands
4Astrophysics Group, Imperial College London, Blackett Laboratory, Prince
Consort Road, London SW7 2AZ, UK
5Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique
de Lyon UMR5574, F-69230 Saint-Genis-Laval, France E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

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^{\prime}$ (which describes the depth of the gap carved by the planet in the
disc). We perform ‘live-planet’ simulations exploring a range of $K^{\prime}$
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^{\prime}$ 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.

###### keywords:

planet-disc interactions – accretion, accretion discs – protoplanetary discs –
hydrodynamics

††pubyear: 2015††pagerange: Inward and outward migration of massive planets:
moving towards a stalling radius–A

##  1 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, 1995; Ivanov et al., 1999), which can be further split
into two regimes, defined through the local disc to planet mass ratio

$B_{0}=\frac{4\pi\Sigma_{0}r^{2}}{m_{\rm p}}.$ (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 planet-dominated 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. (2020); Dempsey 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 $a_{\rm p}$ or in $e_{\rm
p}$, the migration rate for given system parameters is just proportional to
the disc mass; migration is parametrised in terms of $\Delta T/(\dot{M}l_{\rm
p})$, where $\Delta T$ is the torque acting on the planet, $\dot{M}$ is the
steady state accretion rate, $l_{\rm 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 $K^{\prime}=q^{2}/(\alpha
h^{3})$ (which controls the gap depth, Kanagawa et al. 2016), where $q=m_{\rm
p}/M_{*}$ is the planet to star mass ratio, $\alpha$ is the Shakura & Sunyaev
1973 viscosity parameter, and $h=H/r$ is the disc aspect ratio. The parameter
$K^{\prime}$ is a variant of the gap depth parameter $K=q^{2}/(\alpha h^{5})$,
which can be defined from the study of the torques’ balance (see for example
Fung et al., 2014; Kanagawa et al., 2015, 2018). Being proportional to the
planet’s gravitational torque ($\propto q^{2}/h^{3}$) and inversely
proportional to the disc’s viscous torque ($\propto\nu\propto\alpha h^{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 ($B>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 $B<1$ regime have been conducted with a
‘live’ (i.e. evolving) planet.

In this paper, we present a suite of 2D simulations lasting $\gtrsim 300$k
orbits, for a variety of aspect ratios $H/r$, planet masses $M_{p}$ and disc
masses $M_{\rm 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.

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. Name | $B_{0}$ | $m_{\rm p}[m_{\rm J}]$ | $h_{0}$ | $r_{\rm in}$ | $r_{\rm out}$ | $\alpha_{0}$ | $K^{\prime}_{0}$ | $K_{0}$ | $N_{\rm\rm orbits}$
---|---|---|---|---|---|---|---|---|---
L-m1-h036 | $0.046$ | $1$ | $0.036$ | $0.2$ | $15$ | $0.001$ | $21$ | $1.65\cdot 10^{4}$ | $6\cdot 10^{5}$
L-m3-h036 | $0.046$ | $3$ | $0.036$ | $0.2$ | $15$ | $0.001$ | $193$ | $1.49\cdot 10^{5}$ | $3\cdot 10^{5}$
L-m13-h036 | $0.046$ | $13$ | $0.036$ | $0.2$ | $15$ | $0.001$ | $3622$ | $2.9\cdot 10^{6}$ | $3\cdot 10^{5}$
L-m13-h06 | $0.046$ | $13$ | $0.06$ | $0.2$ | $15$ | $0.001$ | $782$ | $2.17\cdot 10^{5}$ | $3\cdot 10^{5}$
L-m13-h1 | $0.046$ | $13$ | $0.1$ | $0.2$ | $15$ | $0.001$ | $169$ | $1.69\cdot 10^{4}$ | $3\cdot 10^{5}$
L-m1-h06 | $0.046$ | $1$ | $0.06$ | $0.2$ | $15$ | $0.001$ | $5$ | $1.29\cdot 10^{3}$ | $6\cdot 10^{5}$
M-m1-h036 | $0.15$ | $1$ | $0.036$ | $0.2$ | $15$ | $0.001$ | $21$ | $1.65\cdot 10^{4}$ | $6\cdot 10^{5}$
M-m3-h036 | $0.15$ | $3$ | $0.036$ | $0.2$ | $15$ | $0.001$ | $193$ | $1.49\cdot 10^{5}$ | $6\cdot 10^{5}$
M-m13-h036 | $0.15$ | $13$ | $0.036$ | $0.2$ | $15$ | $0.001$ | $3622$ | $2.9\cdot 10^{6}$ | $3\cdot 10^{5}$
M-m13-h06 | $0.15$ | $13$ | $0.06$ | $0.2$ | $15$ | $0.001$ | $782$ | $2.17\cdot 10^{5}$ | $3\cdot 10^{5}$
M-m13-h1 | $0.15$ | $13$ | $0.1$ | $0.2$ | $15$ | $0.001$ | $169$ | $1.69\cdot 10^{4}$ | $3\cdot 10^{5}$
M-m1-h06 | $0.15$ | $1$ | $0.06$ | $0.2$ | $15$ | $0.001$ | $5$ | $1.29\cdot 10^{3}$ | $6\cdot 10^{5}$

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.

##  2 Simulations

###  2.1 Simulation parameters

We performed 12 2D hydrodynamical simulations of protoplanetary discs
containing a massive planet with the grid code Fargo 3D (Benítez-Llambay &
Masset, 2016), considering a cylindrical reference frame $(r,\varphi)$
centered in the star. We adopt dimensionless units $G=M_{*}=r_{0}=1$, where
$G$ is the gravitational constant, $M_{*}$ is the star mass and $r_{0}$ is the
planet’s initial location; the time unit is the inverse Keplerian frequency at
$r_{0}$, $\Omega_{\rm k}^{-1}$; this means that the planet at its initial
location $r_{0}$ requires $t=2\pi$ to complete an orbit. Each simulation is
run for $300$k orbits or $600$k, depending on the status of the steady state;
the requirement is that $\dot{M}(r)$ is within 10% of $\dot{M}(r_{\rm in})$ at
least until $r=2.5\ a_{\rm p}$.

The setup for the simulations is defined in analogy with the simulations by
Ragusa et al. (2018).111In turn, their choice of parameters was based on the
best fit for CI Tau disc by Rosotti et al. (2017) We consider logarithmically
spaced $N_{r}=430$ cells in the radial direction, from $r_{\rm in}=0.2$ to
$r_{\rm out}=15$, and linearly spaced $N_{\varphi}=580$ cells in the azimuthal
direction, extending from $0$ to $2\pi$. We assume a locally isothermal
equation of state set by Equation 3, and we parametrise the viscosity by
applying the $\alpha$ prescription by Shakura & Sunyaev (1973), $\nu=\alpha
c_{\rm s}H$, with $\alpha$ defined as a function of radius

$\alpha=\alpha_{0}r^{-0.63};$ (2)

in all the simulations we fix $\alpha_{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

$h=H/r=h_{0}r^{0.215},$ (3)

where $h_{0}$ is the aspect ratio at the initial planet location, for which we
consider 3 different values ($h_{0}=0.036,\ 0.6,\ 0.1$). Regarding the planet
mass, we explore 3 different values: $m_{\rm p}=1\ m_{\rm J}$, $m_{\rm p}=3\
m_{\rm J}$, and $m_{\rm p}=13\ m_{\rm J}$.

For each combination of parameters $h_{0}$ and $m_{\rm p}$, we perform 2
simulations, characterised by different values for the initial local disc to
planet mass ratio $B_{0}$ (see Equation 1). The initial disc surface density
at planet location $\Sigma_{0}$ is chosen to obtain the selected values for
the initial value $B_{0}$; we thus have ‘massive’ disc simulations, with
$B_{0}=0.15$, and ‘light’ disc simulations, with $B_{0}=0.046$. Note that in
physical units, $B_{0}=0.046$ and $B_{0}=0.15$ for a Jovian mass planet at the
radius of Jupiter corresponds to a local disc surface density of $1.3\ \rm
g/cm^{2}$ and $4.2\ \rm g/cm^{2}$, respectively. Apart from $m_{\rm p}$,
$h_{0}$ and $\Sigma_{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 $B_{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 $\times 2$ and $\times 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.

###  2.2 Initial and boundary conditions

The initial density profile is defined as

$\Sigma(r)=\Sigma_{0}r^{-0.3}\cdot{\rm e}^{\left(-{r}/{5}\right)^{1.7}},$ (4)

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 $r_{\rm out}=15$ by
setting both the velocity and the gas density to zero at $r_{\rm 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 $r_{\rm 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. 2020 and Dempsey et al. 2020; Dempsey et al. 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 ($\tau=\Omega_{\rm
k}/30$), we damp $v_{r}$ to the azimuthal average of the initial (viscous)
velocity $\langle v_{r,0}\rangle$ in the region from $r_{\rm in}$ to $r=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 $r_{\rm in}=0.1$, and
the damping zone extends in this case up to $r=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.

##  3 Results

###  3.1 Orbital parameters

Figure 1: Semi-major axis 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). 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).

In Figure 1 we show for the light disc simulations (left panel) and the
massive disc simulations (right panel) the evolution of the semi-major axes as
a function of evolutionary time, i.e. the physical time rescaled to the
viscous timescale at the initial planet position

$t_{\nu}=\frac{2}{3}\frac{r^{2}}{\nu},$ (5)

and we indicate with $t_{\nu,0}$ the viscous timescale evaluated at the
initial planet location. For reference, the initial viscous timescales in our
simulations are: $t_{\nu,0}\sim 9\cdot 10^{5}$ yr for $h=0.036$,
$t_{\nu,0}\sim 3\cdot 10^{5}$ yr for $h=0.06$, $t_{\nu,0}\sim 10^{5}$ yr for
$h=0.1$. After an initial adjustment (lasting $t\lesssim t_{\nu,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 $K^{\prime}=q^{2}/(\alpha h^{3})$ (where
$q=m_{\rm p}/M_{*}$), 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, 2016, 2017; Ragusa et al., 2018;
Teyssandier & Lai, 2019).

###  3.2 Torques

####  3.2.1 Dependence on initial disc mass

Figure 3: Torque acting on the planet normalised for the VSS accretion rate
times the local specific angular momentum as a function of $q$ (left panel)
and parameter $K^{\prime}=q^{2}/(\alpha h^{3})$ (right panel). The squares
indicate the torque acting computed as variation rate of the planet’s angular
momentum obtained from the simulations; the different colours correspond to
different simulations, as indicated in the legend.

In this section, we analyse the torque acting on the planet as a function of
the following parameter (Kanagawa et al., 2015)

$K^{\prime}=\frac{q^{2}}{\alpha h^{3}},$ (6)

which is a modification of the gap-opening parameter $K={q^{2}}/{(\alpha
h^{5})}$, i.e. a measure of the strength of the planet’s gravitational torque
($\propto q^{2}/h^{3}$) compared to the disc’s viscous torque
($\propto\nu\propto\alpha h^{2}$).

The parameter $K^{\prime}$ has been empirically studied by Kanagawa et al.
(2016); Kanagawa et al. (2018) and by Dempsey et al. (2020); Dempsey et al.
(2021), who noticed that $K^{\prime}$ correlates with the gap width and depth
(with the gap becoming wider and deeper for higher $K^{\prime}$ values).
Systems characterised by the same $K^{\prime}$ also appear to show similar
torque density profiles; thus Dempsey et al. (2020); Dempsey et al. (2021)
examined a number of fixed-planet222In fixed-planet simulations, the torque on
the planet is computed as integration of the torque arising from the disc,
from which planet migration is deduced; nonetheless, the planet orbital
parameters are kept fixed at their initial values (i.e. $a_{\rm p}=1$, $e_{\rm
p}=0$). simulations characterised by different $K^{\prime}$ values, and
suggested that $K^{\prime}$ is a good ‘ordering parameter’ for planet
migration. According to their analysis, in fact, planets in systems
characterised by $K^{\prime}\gtrsim 20$ migrate outward, whereas systems with
$K^{\prime}\lesssim 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 $\dot{M}l_{\rm p}$ (see their
Fig. 2). Note that for our convention of signs we have $\dot{M}<0$ for inward
migration; to avoid confusion we therefore plot $\Delta T/|\dot{M}l_{\rm 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

$\Delta T=\frac{{\rm d}J_{\rm p}}{{\rm d}t}=J_{\rm
p}\left(\frac{1}{2a}\cdot\dot{a}-\frac{e}{{1-e^{2}}}\cdot\dot{e}\right),$ (7)

being $J_{\rm p}=m_{\rm p}\sqrt{GM_{*}a_{\rm p}(1-e_{\rm p}^{2})}$. Thus, we
compute $\Delta T$ time-averaged over 30000 planet orbits,333 $30000$ orbits
was chosen as being several times the timescale for periodic exchange of
eccentricity between the planet and the disc, this being modulated on the beat
period between the precession of the aligned and anti-aligned eccentric modes
of the disc.  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 $K^{\prime}$, and negative torques for low
$K^{\prime}$ is confirmed, albeit with a suggestion of some decline in the
torque at values of $K^{\prime}>20-200$. Note, however, that the suite of
simulations performed so far does not explore the $K^{\prime}$ values
extensively enough to confirm that $K^{\prime}=20$ is the point of zero
torque, and further simulations are needed for a precise characterisation of
the $\Delta T-K^{\prime}$ relation.

We further notice an applicability limitation in the results obtained from
fixed planet simulations by Dempsey et al. (2020); Dempsey 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 $B_{0}\rightarrow 0$ (for which the planet is not
expected to modify its orbital parameters noticeably over the disc lifetime),
it cannot necessarily be used for $B_{0}\neq 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 $B_{0}$ values (i.e. ‘light disc’ vs
‘massive disc’), where we notice that the normalised torque properties change
for different $B_{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 $B_{0}$ and hence their ratio is independent of $B_{0}$;
therefore in Dempsey et al. (2020); Dempsey et al. (2021) simulations 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 $K^{\prime}$ ($>10^{3}$).

We find that the normalised torque values are indeed approximately independent
of $B_{0}$ in the case of the Jovian mass planet (see Figure 3) but that for
higher values of planet mass (and $K^{\prime}$) 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
$K^{\prime}$ value (the $13\ m_{\rm 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 $K^{\prime}$ 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 $e_{\rm p}$ rather than $a_{\rm p}$.

Figure 4: Azimuthal averaged density profiles $\Sigma/\Sigma_{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 $a_{\rm p}(1\pm e_{\rm p})$. The left panel refers to simulation
L-m13-h036; the right panel refers to simulation M-m13-h036.

The different evolutionary histories of the disc structure in the light and
massive 13 $m_{\rm 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 the two modes), resulting in the modulation of
planetary eccentricity and disc structure evident in Figures 2 and Figure
4.444In the light disc case with 13 $m_{\rm j}$ planet, the planet achieves
such a high eccentricity that the pericentre distance is only $\sim 2.5$ times
the inner disc edge, $r_{\rm in}$. In order to check that the location of the
inner boundary is not driving the evolution in this case, we have re-run the
light disc case with $r_{\rm in}$ reduced by a factor 2 to 0.1. We show in
Appendix A that the evolution of the planetary orbital elements is unaffected
by the reduction in $r_{\rm in}$.

We therefore conclude that the results based on non-migrating planet
simulations presented by Dempsey et al. (2020); Dempsey 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.

Figure 5: Azimuthal averaged density profiles $\Sigma/\Sigma_{0}$ (colour map)
at the end of the simulation time, plotted as a function of the radius
normalised to the planet location (both the density profile and the planet
location are taken at the final snapshot of each simulation). Each colour
refer to a different simulation, as indicated in the legend; the light and
massive disc simulation with the same planet mass and disc aspect ratio are
shown with the same colour with solid and dashed lines, respectively.

Finally we emphasise that Figure 3 is constructed after significant evolution
of the disc profiles (i.e. after $6\cdot 10^{5}$ planetary orbits in the case
of the massive and light simulations with $1\ m_{\rm j}$ planet, and in the
massive disc simulation with $3\ m_{\rm j}$ planet; and after $3\cdot 10^{5}$
planetary orbits in the remaining simulations) and use values of the torque
and accretion averaged over the preceding $3\cdot 10^{4}$ orbits. The result
that the simulations with higher $K^{\prime}$ values generally result in
normalised torque values that are $\sim 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 $B_{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 $B_{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 $K^{\prime}$ (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. (2020); Dempsey
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.

####  3.2.2 Torque components

Figure 6: Torques associated to the rate of change of semi-major axis
normalised for the VSS accretion rate times the local specific angular
momentum as a function of parameter $K^{\prime}=q^{2}/(\alpha h^{3})$. The
different colours correspond to different simulations, as indicated in the
legend.

Figure 7: Left panel: torque acting on the planet normalised for the VSS
accretion rate times the local specific angular momentum as a function of $K$.
Right panel: torques associated to the rate of change of semi-major axis
normalised for the VSS accretion rate times the local specific angular
momentum as a function of $K$. The different colours correspond to different
simulations, as indicated in the legend.

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 $K^{\prime}$

$\Delta T^{*}=\frac{m_{\rm p}}{2}\cdot\sqrt{\frac{GM_{*}}{a_{\rm p}}(1-e_{\rm
p}^{2})}\cdot\dot{a_{\rm p}}$ (8)

We notice that in this case there is little difference between the torques
associated to $\dot{a}_{\rm p}$ for the light and massive disc simulations
with 13 $m_{\rm j}$ and $h=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 $K^{\prime}$ (or $K$, see the next
section), the normalised $\Delta T^{*}$ is almost constant with $K^{\prime}$
(or $K$), whereas $\Delta T$ 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 $\Delta T^{*}$
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 $\Delta T^{*}$
value whose magnitude ($\sim~{}0.5$) is similar to that in lower mass planets.
This implies that the spread in normalised $\Delta T$ 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

$E_{\rm p}=-\frac{GM_{*}m_{\rm p}}{2a_{\rm p}},$ (9)

the rate of change of planetary orbital energy is given by

$\dot{E}_{\rm p}=\frac{GM_{*}m_{\rm p}}{2a_{\rm p}^{2}}\dot{a}_{\rm p}.$ (10)

Combining Equation 10 and the torque associated with the semi-major axis
variation $|\Delta T^{*}|$ (Equation 8) we deduce that

$\dot{E}_{\rm p}=\left(\frac{\Delta T^{*}}{\dot{M}l_{\rm
p}}\right)\cdot\left(\frac{GM_{*}}{a_{\rm p}}\right)\cdot\dot{M}.$ (11)

Since from our simulations the normalised $\Delta T^{*}$ 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, $\dot{a}_{\rm p}$, scales with
$\dot{M}/m_{\rm p}$ and, for simulations with the same value of $B_{0}$ and
$t_{\nu}$, 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 $\Delta T$ from the constant value at high $K$ (or $K^{\prime}$)
to the planet eccentricity excitation, and we noticed that the effect becomes
more and more important as $K$ ($K^{\prime}$) increases. We plan to further
explore the behaviour at high $K$ (or $K^{\prime}$) in future works.

####  3.2.3 Choice of ordering parameter

The torque analysis conducted so far is following Dempsey et al. (2020);
Dempsey et al. (2021) in assuming that $K^{\prime}$ is the appropriate
ordering parameter. In Figure 7 we instead plot the normalised torques as a
function of

$K=\frac{q^{2}}{\alpha h^{5}}$ (12)

for all our simulations, considering both the total torque (left panel) and
the torque associated with $\dot{a}_{\rm 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 $K=1.5\cdot 10^{4}$. We use this threshold in our ‘toy modeling’
presented in section 5. Remarkably the normalised torque associated with
$\dot{a}$ assumes a constant value of $\sim 0.5$ for higher values of $K$.

##  4 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 $r_{\rm
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 $r$ at different snapshots
for simulation M-m3-h036 with open boundary conditions.555In 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.

Figure 8: Viscous torque at different snapshots (showed in different colours,
as indicated by the legend) for simulation M-m3-h036 with viscous boundary
conditions.

Figure 9: Torque acting on the planet normalised for the VSS accretion rate
times the local specific angular momentum as a function of $q$ (left panel)
and parameter $K^{\prime}=q^{2}/(\alpha h^{3})$ (right panel). The diamonds
indicate the results from the open boundary condition simulations; the squares
illustrate the torque the results from the viscous boundary condition
simulations; the stars show the prediction for the open boundary condition
results, applying Equation 16 to the results obtained from viscous boundary
conditions.

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 $\Delta T^{\rm
in}/|\dot{M}l_{\rm p}|$, where $\Delta T^{\rm in}$ is the torque imparted to
the planet from 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
$r_{\rm in}$ to a few $r_{\rm p}$) reaches a viscous steady state after some
$t_{\nu,0}$ (see Figure 10). If we consider the area of the disc between
$r_{\rm in}$ and $r_{\rm p}$, the total angular momentum injected per unit
time into this region is due to advection and viscous torques at $r_{\rm in}$
and $r_{\rm p}$ is:

$\Delta T^{\rm in}=-\dot{M}(l_{\rm p}-l_{\rm in})+F_{\nu}(r_{\rm
in})-F_{\nu}(r_{\rm p})\sim-\dot{M}(l_{\rm p}-l_{\rm in})+F_{\nu}(r_{\rm
in}),$ (13)

where the approximation is valid for planets massive enough to create a deep
gap in the disc, so that $\Sigma(r_{\rm p})\sim 0$, and thus $F_{\nu}(r_{\rm
p})=3\pi\nu\langle\Sigma\rangle\langle l_{\rm k}\rangle\sim 0$. Since for open
boundary conditions $F_{\nu}(r_{\rm in})\sim 0$, the torque in this case is

$\Delta T_{\rm open\ BC}^{\rm in}=-\dot{M}(l_{\rm p}-l_{\rm in});$ (14)

whereas for viscous boundary conditions $F_{\nu}(r_{\rm in})=-\dot{M}l_{\rm
in}$, thus

$\Delta T_{\rm viscous\ BC}^{\rm in}=-\dot{M}l_{\rm p}.$ (15)

We can therefore estimate the difference in torque imparted by the inner disc

$\Delta T_{\rm open\ BC}^{\rm in}=\Delta T_{\rm viscous\ BC}^{\rm
in}+\dot{M}l_{\rm in},$ (16)

where $\dot{M}<0$ (for our convention of signs), thus $\Delta T_{\rm viscous\
BC}>\Delta T_{\rm open\ BC}$. Since in both cases, the quasi-steady state of
the inner disc requires that $\Delta T$ 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 10: Normalised accretion rate profiles time averaged over different
time intervals (in the legend we indicate the initial and final orbit of the
considered interval) for the massive disc simulation with $m_{\rm p}=3\ m_{\rm
j}$, $h=0.036$ (M-m3-h036). The black dashed lines show the $10\%$ limits on
the value of the normalised accretion rate, in order to consider the region in
steady state.

Figure 9 shows that for the 1 Jupiter mass simulations with $h=0.036$, the
difference in normalised torque is almost exactly that predicted by Equation
16 (with $r_{\rm in}=0.2$). Indeed, the outward migration is significantly
slowed when a zero torque boundary condition is applied at $r_{\rm 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 $13m_{\rm j}$ simulation with $h=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.

##  5 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) $\Delta
T/\dot{M}l_{\rm p}$ is a function of $K$ (Equation 12); (ii) $\Delta
T/\dot{M}l_{\rm p}$ is independent of disc mass.666As underlined in section
3.2.1, $\Delta T/\dot{M}l_{\rm p}$ does not depend on the disc mass only when
the planetary orbit maintains a low eccentricity; this is best satisfied for
the $1\rm m_{\rm J}$ case, which we therefore use in the toy model. Since the
disc is flared, we expect $K\propto h^{-5}$ to vary with radius as $K\propto
r^{-5f}$, where $f$ is the flaring index; in this section we assume constant
$\alpha$, but if we allowed it to vary with radius, it would have given a
further contribution to the dependence of $K$ on $r$.

We compute the planet migration time using the angular momentum definition
$M_{\rm p}l_{\rm p}=M_{\rm p}\sqrt{GM_{*}r_{\rm p}}$777Here we use $r_{\rm p}$
instead of $a_{\rm p}$ because we are assuming circular orbits.

$t_{\rm mig}=\frac{r_{\rm p}}{\dot{r}_{\rm p}}=\frac{1}{2}\cdot\frac{M_{\rm
p}l_{\rm p}}{\Delta T},$ (17)

where $t_{\rm mig}<0$ ($>0$) means inward (outward) migration. Using the
definitions $M_{\rm d}=4\pi r^{2}\Sigma$ and $\dot{M}=3\pi\nu\Sigma$, we can
rewrite the migration timescale as

$t_{\rm mig}=\frac{1}{2}\cdot\left(\frac{|\dot{M}|l_{\rm p}}{\Delta
T}\right)\cdot\left(\frac{M_{\rm p}}{M_{\rm d}}\right)\cdot\left(\frac{M_{\rm
d}}{|\dot{M}|}\right).$ (18)

Noticing that $M_{\rm p}/M_{\rm d}=1/B$ and $M_{\rm d}/\dot{M}=t_{\nu}$888We
define $t_{\nu}={4}{r^{2}}/3{\nu}$ this can be written as

$\frac{r_{\rm p}}{\dot{r}_{\rm p}}=\frac{1}{2}\cdot\left(\frac{|\dot{M}l_{\rm
p}|}{\Delta T}\right)\cdot\left(\frac{t_{\nu}}{B}\right).$ (19)

We now consider 3 different regimes:

  * (a)

classical Type ii migration with $B\gg 1$, for which $r_{\rm p}/\dot{r}_{\rm
p}=-t_{\nu}$, hence we deduce $\Delta T/\dot{M}l_{\rm p}=-0.5/B$ (see, for
example Lin & Papaloizou 1979; D’Angelo et al. 2005; Scardoni et al. 2020);

  * (b)

classical Type ii migration with $B\ll 1$, for which $r_{\rm p}/\dot{r}_{\rm
p}=-t{\nu}/B$, hence we deduce $\Delta T/\dot{M}l_{\rm p}=-0.5$ (see, for
example Syer & Clarke 1995; Ivanov et al. 1999);

  * (c)

migration at high $K$ where $\Delta T/\dot{M}l_{\rm p}\sim 0.5$ from our
numerical results.

We can combine the two limits (a) and (b) as $\Delta T/|\dot{M}l_{\rm
p}|=-0.5/(B+1)$; then we put all of these considerations together to obtain

$\frac{r_{\rm p}}{\dot{r}_{\rm p}}=\begin{dcases}-\frac{B+1}{B}\cdot
t_{\nu}\qquad&{\rm if}\ {K<K_{\rm lim}}\\\ \frac{1}{B}\cdot t_{\nu}\qquad&{\rm
if}\ {K>K_{\rm lim}}\end{dcases}$ (20)

where $K_{\rm lim}$ is the limiting value which determines the transition from
inward to outward migration. For the purpose of the toy model, we take $K_{\rm
lim}=1.5\cdot 10^{4}$ (as estimated from Figure 7).

Given all these considerations, we can define a function $g(K)$ such
that999Note that we define $g(K)$ in this way because we want it to be
proportional to $\Delta T/|\dot{M}l_{\rm p}|$, which is the quantity that we
have analysed in the previous sections.

$\frac{r_{\rm p}}{\dot{r}_{\rm p}}=\frac{1}{{g(K)}}\frac{B+1}{B}t_{\nu}.$ (21)

We therefore need a function with the following properties: for high $K$ and
low $B$, $g(K)\sim 1$, to obtain the lower case in Equation 20; for low $K$ we
want $g(K)\sim-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 $K$ and high $B$, because for the
typical disc parameters and plausible planet location, it is unlikely to
obtain this combination; consequently, for high values of $B$ 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 $g(K)$
as

$g(K)=\tanh\left(\frac{K-K_{\rm lim}}{W}\right),$ (22)

where $W$ controls the width of the transition in $K$ over which $g(K)$
changes sign. We then estimate the planet migration by solving the following
equation

${\dot{r}_{\rm p}=r_{\rm p}g(K)\frac{B}{B+1}\frac{1}{t_{\nu}},}$ (23)

where we underline that all the quantities on the right hand side are
functions of $r_{\rm p}$ (and hence time).

Since planets are expected to migrate inward (outward) for $K$ values smaller
(bigger) than $K_{\rm lim}=1.5\cdot 10^{4}$, and because $K$ increases during
inward migration and decreases during outward migration in the case of flaring
discs ($f>0$ in Equation 24 below), we expect them to always migrate towards
the location of the disc where $K=K_{\rm lim}$. We then call $r_{\rm lim}$ the
radius corresponding to $K_{\rm lim}=1.5\cdot 10^{4}$ which is therefore an
‘attractor’ for migrating planets

${r_{\rm lim}=\left(\frac{q^{2}}{K_{\rm
lim}\alpha_{0}h_{0}^{5}}\right)^{1/(a+5f)},}$ (24)

where $\alpha=\alpha_{0}r^{a}$ and $h=h_{0}r^{f}$.

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.

###  5.1 Toy model in a non-evolving disc

Figure 11: Toy model for a planet of mass $m_{\rm p}=1\ m_{\rm j}$ and a disc
characterised by $\alpha=0.001$, $h=0.02$ at $1\ \rm au$ and flaring index
$0.25$, producing a stalling radius located at $5\ \rm 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 $W$: $W=1$ for
the solid line; $W=1500$ for the dashed line; $W=3000$ for the dotted line. At
5 au, the disc to planet mass ratio is $B=0.1$, while the viscous timescale is
$t_{\nu}(5{\rm au})=1.6\ \rm Myr$

We first define a disc model characterised by $\alpha=0.001$ (constant
throughout the disc), thus $\nu\propto r^{2f+1/2}$; and
$\Sigma(r)=\Sigma(r_{\rm 1})(r/r_{\rm 1})^{-1}$, where $r_{1}$ is the value of
$1\ \rm au$ in code units, and $\Sigma(r_{\rm 1})$ is chosen to have $B_{5\rm
au}=0.1$ at 5 au,101010In the rest of the paper we will refer to the value of
$B$ at 5 au as $B_{5\rm au}$; while we will indicate as $B$ the value at the
planet location. corresponding to $M_{\rm disc}\sim 0.001\ \rm M_{\odot}$ for
a disc of 100 au. As the planet migrates (in either direction), the value of
$B$ varies as111111Note that the dependence of $B$ on $r_{\rm p}$ depends on
the chosen density profile: if $\Sigma\propto r^{-s}$, then $B\propto r_{\rm
p}^{2-s}$.

$B(t)=B_{5\rm au}\cdot\left(\frac{r_{\rm p}(t)}{5\ \rm au}\right).$ (25)

To design a disc model with $r_{\rm lim}=5\ \rm au$ (i.e. the location of
Jupiter), the disc aspect ratio at $1\ \rm 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 $r_{\rm lim}$ as an attractor, considering 3
different models: a planet initially located at $r_{\rm p}(t=0)=3\ \rm au\ <\
r_{\rm lim}$, which then migrates outwards until stalling at $5\ \rm au$ (red
lines); for $r_{\rm p}(t=0)=5\ \rm au\ =\ r_{\rm lim}$ the planet stalls at
its initial location (cyan line); for $r_{\rm p}(t=0)=7\ \rm au\ >\ r_{\rm
lim}$, instead, the planet migrate inwards (magenta lines). The solid, dashed
and dotted lines refer to values of parameter $W$ in Equation 22 equal to $1$,
$1500$, $3000$, respectively.

For the inward migrating planet (initially at $7\ \rm au$) the timescale
required to reach the stalling radius in all cases is less than $10\ \rm 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\ \rm au$, instead,
the planet migrates to $4\ \rm au$ in $10\ \rm Myr$, but requires $\sim 15\
\rm 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 $W$ in the function $g(K)$; this underlines the
importance of further simulations to explore the form of $g(K)$. 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.

###  5.2 Toy model in an evolving disc

Figure 12: Toy model for a Jupiter mass planet intially located at $3$ au
(left panel) and $7$ au (right panel), and migrating in a disc evolving
according to the Lynden-Bell & Pringle (1974) self-similar solution. The
different colours and line styles refer to disc models characterised by a
range of initial scale radii and $B_{5\rm au}$, whose values are indicated in
the legend.

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)

$\Sigma(R,t)=\frac{M_{\rm disc,0}}{2\pi R_{\rm s,0}}(2-b)\left(\frac{R}{R_{\rm
s,0}}\right)^{-b}\tau^{-\eta}\exp\left[-\frac{(R/R_{\rm
s,0})^{2-b}}{\tau}\right],$ (26)

where $M_{\rm disc,0}$ is the initial disc mass, defined by the choice of the
density profile and the initial disc scale radius; $R_{\rm s,0}$ is the
initial scale radius; $b$ 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
consistency with the previous section; $\eta=(5/2-b)/(2-b)$; $\tau$ is defined
as $\tau=1+t/t_{\nu}$, with $t_{\nu}=R_{\rm s}^{2}/[3(2-b)^{2}\nu(R_{\rm
s})]$. The surface density will thus evolve on the viscous timescale of the
instantaneous scale radius

$R_{\rm s}(t)=R_{\rm s,0}\left[\frac{\tau}{2(2-b)}\right]^{\frac{1}{2-b}}.$
(27)

We show in Figure 12 the evolution tracks for a Jupiter mass planet initially
located at $3\ \rm au$ (left panel) and $7\ \rm 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 $B_{0}$ (from 0.09 to 0.5, see the plot legend). For
reference, the initial viscous timescale at the planet location are
$t_{\nu}(3\ {\rm au})\sim 1\ \rm Myr$, and $t_{\nu}(7\ {\rm au})\sim 2.5\ \rm
Myr$; while the viscous timescales at the chosen scale-radii are $t_{\nu}(50\
{\rm au})=4.4\ \rm Myr$, $t_{\nu}(100\ {\rm au})=8.9\ \rm Myr$, $t_{\nu}(200\
{\rm au})=17.7\ \rm 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 $B$ to
decrease, making the migration timescale longer. This means that models
characterised by faster density evolution (i.e. those with lower $R_{\rm
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 $R_{\rm s,0}$ but same $B_{5\rm au}$).

Even more important is the effect of considering different $B_{5\rm au}$
values, with the migration timescale increasing while $B_{5\rm au}$ 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
$B$ 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 $B$, so that it can reach the stalling radius in $\lesssim 10\ \rm Myr$
even in the model characterised by the lowest $B_{5\rm au}$ and the smallest
$R_{\rm 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\ \rm au$ in our model, we need
to take either $B_{5\rm au}\gtrsim 0.4$, or we can decrease it to $B_{5\rm
au}\gtrsim 0.2$ if we take a disc with $R_{\rm s,0}\gtrsim 100\ \rm au$. For
lower $B_{5\rm au}$ values, even in the case of discs with large scale radii,
the planet fails to migrate to the stalling radius in $\lesssim 10\ \rm 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.

###  5.3 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 $r_{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 $T$ at 1 au as a
function of the $\alpha$ 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 $T\gtrsim 10$ K to $T\gtrsim 10^{2}$ K, when we choose values
for the $\alpha$ parameter in the range $\alpha\sim
10^{-4}-10^{-2}$.121212Note that in our model we consider a constant $\alpha$
value, for which it is common to consider $\alpha\sim 10^{-4}-10^{-3}$; if we
consider instead the fact that $\alpha$ depends on the specific location in
the disc, direct turbulence measurements (e.g., Carr et al. 2004) suggest that
at sub-au radii $\alpha\sim 10^{-2}$ is more typical. Since $T\gtrsim
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 $\sim 5$ au from the star. We
further notice that for higher disc temperature the planet stalling radius
decreases for fixed $\alpha$ because for a hotter disc, a larger region of the
disc is in the low $K$ (shallow gap) regime where inward migration is
expected.

If we consider more massive stars, the planet to star mass ratio decreases as
$q\propto M_{*}^{-1}$, whereas the disc aspect ratio
$h_{0}\propto\sqrt{T/M_{*}}$, with $T\propto M_{*}^{\xi}$ (the relation is
sub-linear, and the exact value of $\xi$ depends on the assumed stellar mass-
luminosity relation). Using these scale relations in Equation 24 we find

$r_{\rm lim}\propto M_{*}^{\frac{1-5\xi}{2(a+5f)}},$ (28)

therefore we expect the stalling radius value to increase (decrease) with
increasing stellar mass for $\xi<1/5$ ($\xi>1/5$). At fixed radius, the
temperature approximately scales as $T\propto M_{*}^{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.

Figure 13: Disc temperature at 1 au as a function of the viscosity parameter
$\alpha$, required to have a Jupiter mass planet stalling at $r_{0}$; the
magenta line correspond to $r_{0}=3\ \rm au$, the blue line correspond to
$r_{0}=5\ \rm au$, the cyan line correspond to $r_{0}=10\ \rm au$.

##  6 Implications

Planetary population synthesis models for giant planets are widely based on
the theory of Type II migration (Ida & Lin, 2004, 2008; Mordasini et al.,
2009a; Mordasini et al., 2009b; 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. 2020; Dempsey 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.
(2020); Dempsey 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 $\dot{M}l_{\rm 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 $\dot{M}l_{\rm 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 $h$), 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\cdot 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 $K$ 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
($t_{\nu_{0}}/B$) 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 $B$ values
(corresponding to fast migration timescale), and high $R_{\rm s,0}$ values
(corresponding to slow viscous evolution). Given these considerations, which
depend on poorly constrained parameters such as the value of $\alpha$ 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.131313Note 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 This is because it is implausible that the
disc gas surface density in the inner disc would ever be high enough to
maintain $B>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, $B=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,
2014), 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.

##  7 Conclusions

In this work we have presented a suite of long term ($300-600\rm k$ planet
orbits), 2D hydrodynamical simulations to test and expand recent findings by
Dempsey et al. (2020); Dempsey et al. (2021) of an empirical correlation
between the direction of planet migration and the value of the modified gap-
opening parameter $K^{\prime}$. 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 $K^{\prime}=q^{2}/(\alpha h^{3})$ (Equation 6) and $K=q^{2}/(\alpha
h^{5})$ (Equation 12). Our simulations suggest that $K$ 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
$K_{\rm lim}=1.5\cdot 10^{4}$.

We further notice that, regardless the choice of the ordering parameter, as
the planet mass increases some dependence on $B_{0}$ of the torques acting on
the planet is revealed; this effect could not be seen in fixed planet
simulations, where the normalised torque $\Delta T/|\dot{M}l_{\rm 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. (2020); Dempsey
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 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 $K_{\rm
lim}\sim 1.5\cdot 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 $K$ is a decreasing function of radius. We thus predict inward
(outward) migration for radii exterior (interior) to the location where $K$
has its limiting value.

We then model the migration behaviour of massive planets by describing the
dependence of migration on the parameter $K$. 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 $K<K_{\rm lim}$ ($K>K_{\rm
lim}$), they tend to go towards the location in the disc where $K=K_{\rm
lim}$. We thus propose the existence of a ‘stalling radius’ defined by the
location where $K=K_{\rm 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 $\alpha$ 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 $B$). 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 $B_{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 $K$ in the
vicinity of $K_{\rm 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.

## Acknowledgements

We would like to thank the referee for the careful reading of this paper and
for the useful comments. CES thanks Peterhouse Cambridge for a Ph.D.
studentship. RAB acknowledges the Royal Society University Research
Fellowship. GR acknowledges support from the Netherlands Organisation for
Scientific Research (NWO, program number 016.Veni.192.233) and from an STFC
Ernest Rutherford Fellowship (grant number ST/T003855/1). RDA and ER
gratefully acknowledges funding from the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation programme (grant
agreement no. 681601); RDA also acknowledges funding from STFC Consolidated
Grant ST/W000857/1. ER acknowledges funding from the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant agreement No 864965). This work has also been supported by
the European Union’s Horizon 2020 research and innovation programme under the
Marie Sklodowska-Curie grant agreement number 823823 (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 (Fargo3D,
Benítez-Llambay & Masset 2016) is available at http://fargo.in2p3.fr.

## References

  * Alcalá et al. (2014) Alcalá J. M., et al., 2014, A&A, 561, A2
  * Artymowicz (1993) Artymowicz P., 1993, ApJ, 419, 166
  * Baruteau et al. (2011) Baruteau C., Meru F., Paardekooper S. J., 2011, in Alecian G., Belkacem K., Samadi R., Valls-Gabaud D., eds, SF2A-2011: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics. pp 459–461
  * Baruteau et al. (2014) Baruteau C., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 667 (arXiv:1312.4293), doi:10.2458/azu_uapress_9780816531240-ch029
  * Benítez-Llambay & Masset (2016) Benítez-Llambay P., Masset F. S., 2016, ApJS, 223, 11
  * Bitsch et al. (2013) Bitsch B., Crida A., Morbidelli A., Kley W., Dobbs-Dixon I., 2013, A&A, 549, A124
  * Bitsch et al. (2019) Bitsch B., Izidoro A., Johansen A., Raymond S. N., Morbidelli A., Lambrechts M., Jacobson S. A., 2019, A&A, 623, A88
  * Carr et al. (2004) Carr J. S., Tokunaga A. T., Najita J., 2004, ApJ, 603, 213
  * Crida & Morbidelli (2007) Crida A., Morbidelli A., 2007, MNRAS, 377, 1324
  * Crida et al. (2006) Crida A., Morbidelli A., Masset F., 2006, Icarus, 181, 587
  * D’Angelo et al. (2005) D’Angelo G., Bate M. R., Lubow S. H., 2005, MNRAS, 358, 316
  * Dempsey et al. (2020) Dempsey A. M., Lee W.-K., Lithwick Y., 2020, ApJ, 891, 108
  * Dempsey et al. (2021) Dempsey A. M., Muñoz D. J., Lithwick Y., 2021, arXiv e-prints, p. arXiv:2105.05277
  * Donati et al. (2020) Donati J. F., et al., 2020, MNRAS, 491, 5660
  * Duffell & Chiang (2015) Duffell P. C., Chiang E., 2015, ApJ, 812, 94
  * Duffell et al. (2014) Duffell P. C., Haiman Z., MacFadyen A. I., D’Orazio D. J., Farris B. D., 2014, ApJ, 792, L10
  * Dürmann & Kley (2015) Dürmann C., Kley W., 2015, A&A, 574, A52
  * Fernandes et al. (2019) Fernandes R. B., Mulders G. D., Pascucci I., Mordasini C., Emsenhuber A., 2019, ApJ, 874, 81
  * Fung & Chiang (2016) Fung J., Chiang E., 2016, ApJ, 832, 105
  * Fung et al. (2014) Fung J., Shi J.-M., Chiang E., 2014, ApJ, 782, 88
  * Goldreich & Sari (2003) Goldreich P., Sari R., 2003, ApJ, 585, 1024
  * Goldreich & Tremaine (1979) Goldreich P., Tremaine S., 1979, ApJ, 233, 857
  * Goldreich & Tremaine (1980) Goldreich P., Tremaine S., 1980, ApJ, 241, 425
  * Hallam & Paardekooper (2018) Hallam P. D., Paardekooper S. J., 2018, MNRAS, 481, 1667
  * Hasegawa & Pudritz (2013) Hasegawa Y., Pudritz R. E., 2013, ApJ, 778, 78
  * Hasegawa & Pudritz (2014) Hasegawa Y., Pudritz R. E., 2014, ApJ, 794, 25
  * Ida & Lin (2004) Ida S., Lin D. N. C., 2004, ApJ, 604, 388
  * Ida & Lin (2008) Ida S., Lin D. N. C., 2008, ApJ, 673, 487
  * Ida et al. (2018) Ida S., Tanaka H., Johansen A., Kanagawa K. D., Tanigawa T., 2018, ApJ, 864, 77
  * Ivanov et al. (1999) Ivanov P. B., Papaloizou J. C. B., Polnarev A. G., 1999, MNRAS, 307, 79
  * Johns-Krull et al. (2016) Johns-Krull C. M., et al., 2016, ApJ, 826, 206
  * Kanagawa et al. (2015) Kanagawa K. D., Tanaka H., Muto T., Tanigawa T., Takeuchi T., 2015, MNRAS, 448, 994
  * Kanagawa et al. (2016) Kanagawa K. D., Muto T., Tanaka H., Tanigawa T., Takeuchi T., Tsukagoshi T., Momose M., 2016, PASJ, 68, 43
  * Kanagawa et al. (2018) Kanagawa K. D., Tanaka H., Szuszkiewicz E., 2018, ApJ, 861, 140
  * Kley & Dirksen (2006) Kley W., Dirksen G., 2006, A&A, 447, 369
  * Kley & Nelson (2012) Kley W., Nelson R. P., 2012, ARA&A, 50, 211
  * Lega et al. (2021) Lega E., et al., 2021, A&A, 646, A166
  * Lin & Papaloizou (1979) Lin D. N. C., Papaloizou J., 1979, MNRAS, 186, 799
  * Lin et al. (1993) Lin D. N. C., Papaloizou J. C. B., Kley W., 1993, ApJ, 416, 689
  * Lubow & D’Angelo (2006) Lubow S. H., D’Angelo G., 2006, ApJ, 641, 526
  * Lynden-Bell & Pringle (1974) Lynden-Bell D., Pringle J. E., 1974, MNRAS, 168, 603
  * Malik et al. (2015) Malik M., Meru F., Mayer L., Meyer M., 2015, ApJ, 802, 56
  * Manara et al. (2016) Manara C. F., et al., 2016, A&A, 591, L3
  * Mayor & Queloz (1995) Mayor M., Queloz D., 1995, Nature, 378, 355
  * Mordasini et al. (2009a) Mordasini C., Alibert Y., Benz W., 2009a, A&A, 501, 1139
  * Mordasini et al. (2009b) Mordasini C., Alibert Y., Benz W., Naef D., 2009b, A&A, 501, 1161
  * Nielsen et al. (2019) Nielsen E. L., et al., 2019, AJ, 158, 13
  * Paardekooper et al. (2022) Paardekooper S.-J., Dong R., Duffell P., Fung J., Masset F. S., Ogilvie G., Tanaka H., 2022, arXiv e-prints, p. arXiv:2203.09595
  * Papaloizou (2021) Papaloizou J. C. B., 2021, arXiv e-prints, p. arXiv:2107.07269
  * Papaloizou & Terquem (2006) Papaloizou J. C. B., Terquem C., 2006, Reports on Progress in Physics, 69, 119
  * Ragusa et al. (2018) Ragusa E., Rosotti G., Teyssandier J., Booth R., Clarke C. J., Lodato G., 2018, MNRAS, 474, 4460
  * Robert et al. (2018) Robert C. M. T., Crida A., Lega E., Méheut H., Morbidelli A., 2018, A&A, 617, A98
  * Rosotti et al. (2017) Rosotti G. P., Booth R. A., Clarke C. J., Teyssandier J., Facchini S., Mustill A. J., 2017, MNRAS, 464, L114
  * Rowther & Meru (2020) Rowther S., Meru F., 2020, MNRAS, 496, 1598
  * Scardoni et al. (2020) Scardoni C. E., Rosotti G. P., Lodato G., Clarke C. J., 2020, MNRAS, 492, 1318
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, in Bradt H., Giacconi R., eds, IAU Symposium Vol. 55, X- and Gamma-Ray Astronomy. p. 155
  * Sinclair et al. (2020) Sinclair C. A., Rosotti G. P., Juhasz A., Clarke C. J., 2020, MNRAS, 493, 3535
  * Stamatellos & Inutsuka (2018) Stamatellos D., Inutsuka S.-i., 2018, MNRAS, 477, 3110
  * Syer & Clarke (1995) Syer D., Clarke C. J., 1995, MNRAS, 277, 758
  * Teyssandier & Lai (2019) Teyssandier J., Lai D., 2019, MNRAS, 490, 4353
  * Teyssandier & Ogilvie (2016) Teyssandier J., Ogilvie G. I., 2016, MNRAS, 458, 3221
  * Teyssandier & Ogilvie (2017) Teyssandier J., Ogilvie G. I., 2017, MNRAS, 467, 4577
  * Thun et al. (2017) Thun D., Kley W., Picogna G., 2017, A&A, 604, A102
  * Veras & Armitage (2004) Veras D., Armitage P. J., 2004, MNRAS, 347, 613
  * Ward (1998) Ward W. R., 1998, in Woodward C. E., Shull J. M., Thronson Harley A. J., eds, Astronomical Society of the Pacific Conference Series Vol. 148, Origins. p. 338
  * de Val-Borro et al. (2006) de Val-Borro M., et al., 2006, MNRAS, 370, 529

##  Appendix A Test inner boundary effects

Figure 14: Semi-major axis (left panel) and eccentricity (right panel)
evolution in the light disc case with $13\ m_{\rm j}$ planet and $h=0.036$.
The magenta line shows the planet evolution in the original simulation run
L-m13-h=0.026, with $r_{\rm in}=0.2$; the black dashed line shows the result
from the simulation test, characterised by $r_{\rm in}=0.1$.

From Figure 2 we noticed that in simulation L-m13-h036 (i.e. the light disc
case with $13\ m_{\rm j}$ planet and $h=0.036$), the planet grows
significantly its eccentricity so that its pericentre distance is only $\sim
2.5$ times the inner disc edge. Since at the inner boundary we are applying
wave killing boundary conditions (see section 2), we must ensure that at the
in simulation L-m13-h036 the planet evolution is not affected by numerical
effects due to its closeness to the inner boundary when it is at the
pericentre. We thus present here the results from a simulation test,
characterised by all the same characteristics as those of simulation
L-m13-h036, apart from the inner radius and the damping region, which are both
taken to be half of their original values: we take $r_{\rm in}=0.1$, the
damping is taken from $r_{\rm in}$ to $r=0.15$; we also increase the number of
radial cells in order to obtain the same radial resolution as the original
simulation run.

In Figure 14 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 $r_{\rm 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.