# ALMA constraints on assembly of Core Accretion planets

Sergei Nayakshin1, Vardan Elbakyan1,2 and Giovanni Rosotti1
1 School of Physics and Astronomy, University of Leicester, Leicester, LE1
7RH, UK.
2 Institute of Astronomy, Russian Academy of Sciences, Pyatnitskaya str. 48,
Moscow 119017, Russia<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Resolved dust continuum and CO line ALMA imaging, and in some cases detection
of H$\alpha$ emission, hint that young massive planets are abundant at wide
separations in protoplanetary discs. Here we show how these observations can
probe the runaway phase of planetary growth in the Core Accretion theory.
Planets in this phase have the right range of masses to account for the
predominantly moderate contrast gaps and rings seen in ALMA observations.
However, we find that these planets gain mass and migrate inward very rapidly.
As a result, the phase when they could produce gaps with properties similar to
those observed is very short, i.e., $t_{\rm gap}\lesssim 0.1$ Myr,
independently of the disc viscosity parameter. This would require many tens to
hundreds of gas giant planets to be born per ALMA system, violating the
available mass budget of solids in realistic discs. This also predicts
preponderance of discs with very wide gaps or complete inner disc holes, which
is not observed. We show that suppression of both planet accretion and
migration by a factor of at least ten is a possible solution to these serious
problems. Future population synthesis models of planet formation should aim to
address both exoplanetary data of older discless planetary systems and ALMA
discs with embedded planets in one framework.

###### keywords:

planet-disc interactions – protoplanetary discs – planets and satellites:
formation

††pubyear: 2022††pagerange: ALMA constraints on assembly of Core Accretion
planets–10

##  1 Introduction

A planet formation theory is a set of equations describing evolution of planet
variables, such as its mass $M_{\rm p}$, separation $a$, metallicity $z_{\rm
p}$, etc., inside the parent protoplanetary disc. Protoplanetary disc
evolution is itself a field of intense research (Williams & Cieza, 2011;
Alexander et al., 2014), with data indicating a tremendous variety in
protoplanetary disc properties even for stars of similar mass, age, and
environment (e.g., Manara et al., 2017; Andrews, 2020). Furthermore,
exoplanetary data are very strongly dominated by planets orbiting stars much
older than $\sim 10$ Myr (e.g., Winn & Fabrycky, 2015). These stars have lost
all signs of their primordial discs from which they grew. This makes it
impossible to ascertain the properties of the disc in which a given planet
formed. We also face similar uncertainties for the planets in the Solar
System. Therefore, in the last several decades, to compare planet formation
model predictions to observations theorists had to make do with statistical
assumptions about the disc structure and parameters, and then run their models
until disc dispersal (e.g., Pollack et al., 1996; Ida & Lin, 2004a; Mordasini
et al., 2012). This situation is less than ideal, leaving much wiggle room for
planet formation theories to contend with the data.

The recent protoplanetary disc imaging revolution (e.g., ALMA Partnership et
al., 2015) is a game changer. Instead of making statistical guesses about the
properties of protoplanetary discs we can now constrain these from dust disc
morphology, grain size, molecular line disc imaging, gas accretion rate onto
the star, etc. Unlike exoplanets orbiting discless stars, ALMA candidate
planets (e.g., ALMA Partnership et al., 2015; Andrews et al., 2018; Long et
al., 2018; Pinte et al., 2020) must be evolving right now, rather than in
their long forgotten past, promising more direct observational constraints on
the functions $F_{i}$.

In this paper we illustrate the power of this approach by focusing on a key
step in the Core Accretion (CA) scenario for planet formation (e.g., Safronov,
1969; Pollack et al., 1996). In all current flavours of CA, when a planetary
(solid) core reaches the critical mass of $\sim 10{\,{\rm M}_{\oplus}}$ while
still embedded in the gaseous protoplanetary disc, a very rapid – and thus
called runaway – gas accretion onto the planet sets in (e.g., Mizuno, 1980;
Stevenson, 1982; Ikoma et al., 2000). The runaway terminates when the planet
reaches $\sim$ Jovian mass and opens a deep gap in the surrounding gaseous
disc (Bate et al., 2003). For a population of growing planets, very few
planets are in the runaway phase at any given time; most are either too low
mass or already too massive to be accreting gas rapidly. Accordingly, CA
theories predict a pronounced planet desert in the runaway mass range, that
is, from $\sim$ twice Neptune ($0.1{\,{\rm M}_{\rm J}}$) to $\sim$ Jupiter
($1{\,{\rm M}_{\rm J}}$) masses111Note that the gap opening mass scales
roughly as $\propto(H/R)^{3}$. Since discs are geometrically thicker at larger
radii, the maximum planet mass in the CA scenario increases with separation
from the star. (Ida & Lin, 2004a, b), although the amplitude of the desert
depends on model detail (e.g., Mordasini et al., 2009a; Mordasini et al.,
2009b).

Early analysis of radial velocity CORALIE/HARPS exoplanetary data appeared to
confirm this theoretical prediction (Mayor et al., 2011). More recently,
microlensing observations (e.g., Suzuki et al., 2016, 2018; Jung et al., 2019)
and re-analysis of the CORALIE/HARPS sample (Bennett et al., 2021) indicated
instead a single power-law mass function from $\sim 10{\,{\rm M}_{\oplus}}$
into the Jovian regime. Furthermore, modern 3D numerical simulations find much
smaller gas accretion rates onto planets (Szulágyi et al., 2014; Ormel et al.,
2015; Fung & Chiang, 2016; Lambrechts & Lega, 2017; Lambrechts et al., 2019)
than previous 1D studies found.

Using the census of planets growing in their discs compiled by Lodato et al.
(2019), Nayakshin et al. (2019) showed that ALMA planets must be growing
$\sim$ an order of magnitude slower than CA scenario predicts. Nayakshin et
al. (2019) used idealised power-law gas-only models for ALMA discs and
neglected planet migration. In this paper we extend this work in several
significant ways. We include planet migration, and calculate the time-
dependent structure of discs via a 1D model that includes a self-consistent
gravitational torque exchange between the disc and the planet, including disc
gap opening. We also include gas accretion onto planets by explicitly removing
mass from the disc, thus combining constraints from planet accretion and
migration. Finally, dust dynamics in the disc is solved for, enabling analysis
of time evolution of planet-induced dust gaps, inner holes, and rings.

Our paper is organised as following. In §2 we introduce our methods, in §3 we
apply these methods to a well known protoplanetary disc in HD163296 to
illustrate how planets can be used to constrain CA scenario. In §4 we apply
these ideas to a sample of ALMA gaps, and in §5 we present a discussion of
main results of our paper.

##  2 Method

###  2.1 Gas, dust and planet dynamics

Due to numerical limitations, population synthesis modelling of planet
formation has so far only been performed in 1D (e.g., Ida & Lin, 2004a;
Mordasini et al., 2009a; Coleman & Nelson, 2014; Nayakshin & Fletcher, 2015;
Ndugu et al., 2019). For similar reasons, and to be directly comparable with
previous work, we follow this approach here, employing 1D geometrically thin
azimuthally averaged disc model.

We are interested in the outer regions of protoplanetary discs with accretion
rates usually well below $10^{-7}{\rm M_{\odot}}$ yr-1. Under such conditions
the discs are passively heated by irradiation from the central star, and the
midplane temperature scale with radius $R$ approximately as $T\propto
R^{-1/2}$ (D’Alessio et al., 2001). We therefore set

$T=T_{0}\left(\frac{R_{0}}{R}\right)^{1/2}$ (1)

where $R_{0}=100$ AU and $T_{0}$ is specific for each source we model. The
disc geometric aspect ratio is $H/R=h_{0}(R/R_{0})^{1/4}$, where
$h_{0}=(k_{b}T_{0}R_{0}/GM_{*}\mu)^{1/2}$, $k_{b}$ is the Boltzmann’s constant
and $\mu=2.4m_{p}$ is the mean molecular weight, and $M_{*}$ is the stellar
mass.

Time evolution of the gaseous disc surface density is described by the
standard viscous model, which reads

$\frac{\partial\Sigma}{\partial t}=\frac{3}{R}\frac{\partial}{\partial
R}\left[R^{1/2}\frac{\partial}{\partial
R}(R^{1/2}\nu\Sigma)\right]-\frac{1}{R}\frac{\partial}{\partial
R}\left(Rv_{\lambda}\Sigma\right)-\dot{\Sigma}_{-}\;.$ (2)

In this equation, the terms on the right are: the viscous diffusive mass
transfer, the radial mass flow due to the tidal torque from the planet, and
the mass loss term due to accretion onto embedded objects. The velocity
$v_{\lambda}=2\lambda_{t}/(R\Omega_{\rm K})$ may be thought as the velocity
imposed on the gas by the specific tidal torque $\lambda_{t}$ from the planet.
The integral of the tidal torque over the whole disc gives the total rate of
the angular momentum gain of the disc due to planet-disc interaction,
$\Lambda_{\rm t}$ (cf. eq. 22). We note in passing that $\Lambda_{\rm t}$ can
be positive (the planet migrates inward) or negative (the planet migrates
outward). The conservation of the angular momentum in the disc-planet system
demands that the planet-star separation $a$ evolves as

$\frac{da}{dt}=-\frac{2\Lambda_{\rm t}}{M_{\rm
p}}\left(\frac{a}{GM_{*}}\right)^{1/2}=-\frac{a}{t_{\rm mig}}\;,$ (3)

where we defined the planet migration time, $t_{\rm mig}$. We use the Crida et
al. (2006) parameter

$C_{\rm p}={3H\over 4R_{\rm H}}+50\alpha_{\rm v}\left({H\over
a}\right)^{2}{M_{*}\over M_{p}}$ (4)

to differentiate between type I planet migration (no deep gap in the gas
surface density at the location of the planet), $C_{\rm p}\geq 1$, and type II
regime (a deep gap is opened), $C_{\rm p}\leq 1$. In the type I regime we use
Paardekooper et al. (2011) expressions for $\Lambda_{\rm t}$, whereas in type
II regime the tidal torque is an integral (cf. eq. 22) over the specific
torque $\lambda_{t}$, given by the widely used expression in the literature,
(Lin & Papaloizou, 1986; Armitage & Bonnell, 2002; Lodato & Clarke, 2004;
Alexander et al., 2006; Lodato et al., 2009; Dipierro & Laibe, 2017):

$\displaystyle\lambda_{t}=$ $\displaystyle\displaystyle
f\frac{q^{2}}{2}\left(\displaystyle\frac{a}{\Delta
R}\right)^{4}\frac{GM_{*}}{a}$ $\displaystyle\text{ for }R>a$ (5)
$\displaystyle\lambda_{t}=$
$\displaystyle-f\displaystyle\frac{q^{2}}{2}\left(\displaystyle\frac{R}{\Delta
R}\right)^{4}\frac{GM_{*}}{a}$ $\displaystyle\text{ for }R<a\;,$

where $\Delta R=|R-a|$ and the factor $f$ is a constant set to $f=0.5$ here.
The torque term diverges at the planet location, and hence one smooths it for
$|\Delta R|<\max[H,R_{H}]$ (e.g., Syer & Clarke, 1995). As a number of authors
in the past we use a prescription to transition smoothly from type I to type
II planet migration (see appendix A for detail).

The protoplanetary disc viscosity is given by the Shakura & Sunyaev (1973)
model, $\nu=\alpha c_{s}H$, where $c_{s}=(P/\rho)^{1/2}$ is the isothermal
sound speed, with $P$ and $\rho$ being the disc midplane pressure and density.

We neglect dust back reaction on the gas and assume that dust radial velocity
$v_{\rm d}$ is given by the terminal velocity approximation plus the diffusion
term (see, e.g., Rosotti et al., 2016). The dust surface density then obeys
the mass continuity equation, modulo mass loss due to planet accretion (when
included):

$\frac{\partial\Sigma_{\rm d}}{\partial
t}=-\frac{1}{R}\frac{\partial}{\partial R}\left(Rv_{\rm d}\Sigma_{\rm
d}\right)-\dot{\Sigma}_{\rm d}\;,$ (6)

where

$\dot{\Sigma}_{\rm d}=\frac{\Sigma_{\rm d}}{\Sigma}\dot{\Sigma}_{-}\;.$ (7)

The latter equation merely states that dust is carried with the gas into the
planet during the runaway gas accretion phase onto it. This is appropriate
since our dust particles are usually in the small Stokes number regime. We do
not include pebble accretion onto our planets as this is important only for
planets before they enter the gas runaway phase (e.g., Ndugu et al., 2019) and
we are not concerned here with how exactly our planets gained enough solids
(via planetesimals or pebbles) to get into the runaway regime. We solve eqs. 2
and 6 in an explicit manner.

Protoplanetary disc surface density is often modelled via a power-law with an
exponential rollover,

$\Sigma=\Sigma_{0}\frac{R_{0}}{R}\exp\left[-\frac{R}{R_{c}}\right]\;,$ (8)

where $\Sigma_{0}$ is a normalisation constant, $R_{\rm c}$ is the critical
radius (e.g., in HD163296 it is 150 AU, see Dullemond et al., 2020). The
values of $R_{\rm c}$ for most of the observed discs are not well constrained,
but fortunately are not particularly important for our modelling here. Dust
continuum observations only constrain the outer edge of large dust particles
disc (e.g., Rosotti et al., 2019). The location of the outer edge of the gas
disc is harder to pinpoint but is likely to be $\sim$ a few times that of the
dust disc. Since the planets we are interested here are by definition within
the dust disc, $R_{\rm c}$ is probably much larger than $a$ in these systems.

For numerical convenience, and in line with a common approach in the field
(e.g., Rosotti et al., 2016; Zhang et al., 2018) we use the following initial
gas surface density profile:

$\Sigma=\Sigma_{0}\frac{R_{0}}{R}\;,$ (9)

sharply cut at $R=R_{\rm out}=200$ AU, where we place a reflecting boundary
condition. These simplifications do not affect our main conclusions but are
numerically convenient. Coupled with $T(R)$ given by eq. 1, this density
profile results in a radially constant gas accretion rate in the disc,
$\dot{M}$. To keep our disc in this state we further keep the gas and dust
surface densities at the outer boundary fixed, while allowing the dust and gas
to freely pass though the inner boundary (such boundary conditions are
sometimes called “evanescent”; Zhang et al., 2018).

###  2.2 Planet accretion

To model gas accretion term $\dot{\Sigma}_{-}$ onto our planets we follow the
well known Bern model for planet formation (e.g., Emsenhuber et al., 2021a,
b). In the Bern model, the gas capture radius of the planet is defined as

$R_{\rm gc}=\frac{k_{\rm acc}R_{\rm H}R_{\rm acc}}{R_{\rm acc}+k_{\rm
acc}R_{\rm H}}$ (10)

where $R_{\rm acc}=GM_{\rm p}/c_{s}^{2}$ and $k_{\rm acc}=1/4$. Gas accretion
proceeds in the 3D regime for low mass planets, when $R_{\rm gc}\leq H$, and
in 2D regime for high mass planets when $R_{\rm gc}>H$. The maximum (runaway)
2D gas accretion rate onto the planet in Emsenhuber et al. (2021a) model is

$\dot{M}_{2D}=2R_{\rm gc}\Sigma(a)v_{\rm rel}\;,$ (11)

where $v_{\rm rel}=$ max$[\Omega_{K}R_{\rm gc},c_{s}]$.

In the Bern model, the gap opening in the disc by the planet is not computed
from first principles since the planet tidal torque is not applied to the disc
while the planet is in the type I migration regime. The motivation for this
approach is the fact that the planet affects the gas disc weakly except very
close to it, so $\Sigma(a)$ in eq. 11 is approximately the unperturbed gas
surface density. The Bern model then assumes that a deep gap in the disc
emerges once $C_{\rm p}$ drops below unity. We note that this is a common
approximation in population synthesis for planet formation (e.g., Coleman &
Nelson, 2014).

Due to the aims of our paper, however, gap formation in the gas disc must be
explicitly resolved since dust is usually affected by the planets stronger
than gas (Dipierro et al., 2016), so that even planets with $C_{\rm p}>1$ can
open deep gaps in the discs (Dipierro & Laibe, 2017). Therefore we use the
actual $\Sigma(R)$ in eq. 11 rather than the unperturbed value. This reduces
the accretion rate onto the planet when it starts opening a gap, making the
main conclusions of our paper more conservative.

Since $\Sigma(R)$ may be changing rapidly near the planet location in our more
physically complete model, to find the gas accretion rate onto the planet,
$\dot{M}_{\rm acc}$, instead of using the unperturbed gas density in eq. 11,
we take the following integral:

$\dot{M}_{\rm acc}=f_{\rm geom}\int_{a-R_{\rm H}}^{a+R_{\rm
H}}dR\Sigma(R)v_{\rm gc}\;,$ (12)

where the gas capture velocity $v_{\rm gc}=$ max $[c_{s},\Omega_{K}|R-a|]$,
and $f_{\rm geom}=$ min $[R_{\rm gc}/H,1]$ is the geometrical factor that is
smaller than unity in the 3D accretion case. For the case of no gap opened in
the disc, within a factor of order unity, the result of eq. 12 is the same as
eq. (47) in Emsenhuber et al. (2021a). When the planet mass increases and the
gap begins to open, however, eq. 12 self-consistently reduces gas accretion
onto the planet, eventually terminating its growth and preventing the planet
from accreting the whole gas disc. We further comment that in the Bern model
the planets may be on eccentric orbits due to their interactions with other
planets, which increases the feeding zone of the planet and may lead to it
reaching higher terminal masses (cf. eqs. 49 & 50 in Emsenhuber et al.,
2021a). To be conservative in our conclusions, here we do not allow eccentric
orbits. Relaxing this assumption would make the need to suppress gas accretion
onto planets even stronger.

##  3 An Example: the planet at 86 AU in HD163296

To exemplify the range of results amid parameter dependencies we focus on one
particular case before turning our attention to the Lodato et al. (2019)
sample of candidate planets in ALMA discs. HD 136296 is one of the closest and
therefore one of the best studied Herbig Ae stars, with mass $M_{*}=2.04{\rm
M_{\odot}}$, luminosity $L_{*}=17\rm\,L_{\odot}$, and an estimated age of
$t_{*}\approx 5$ Myr (Andrews et al., 2018). Its protoplanetary disc is one of
a very few for which CO and DCO+ molecular line observations with ALMA are of
sufficiently high quality to set an upper limit to the disc turbulence
parameter $\alpha_{\rm v}\lesssim 3\times 10^{-3}$ (Flaherty et al., 2017). We
set $T_{0}=65$ K at $R_{0}=100$ AU (Fig. 6 in Rab et al., 2020). Note that
this temperature is found at the CO molecule emitting surface, and is larger
than $T\approx 25$ K temperature derived for the midplane of the disc in this
source by Dullemond et al. (2020). However, by picking the larger of the
derived temperatures here we obtain the most conservative results in the
context of our vertically and azimuthally averaged modelling. As we show in
Appendix A, the lower the disc temperature, the stronger the ALMA constraints
on the Core Accretion scenario. This is because at lower disc temperatures gas
accretion rate onto the planet is higher, the type I migration rate is faster
at the same disc mass, and a deep gap opening in the gas disc occurs sooner.

The putative planet we shall focus on produces a well defined gap at $a_{\rm
p}=86$ AU, characterised in Zhang et al. (2018) by $\Delta_{\rm gap}=\Delta
R/a=0.17$, where $\Delta R$ is the gap radial width, and the gap/ring
intensity contrast $C_{\rm gap}=0.15$ (for definitions see Dullemond et al.,
2018; Zhang et al., 2018, and §3.2 below) The case for this candidate planet
has also been strengthened by the detection of velocity perturbations both in
the vertical direction in the disc (Teague et al., 2019) and in the disc plane
direction (velocity kink) visible through the CO molecular line observations
(Pinte et al., 2018, 2020). The disc mass in this source has been estimated
through dust mass conversion to gas (assuming 100 for the gas-to-dust ratio)
and also through disc chemodynamical modelling at $M_{\rm disc}=(0.05-0.3){\rm
M_{\odot}}$ (Qi et al., 2011; Muro-Arena et al., 2018; Booth et al., 2019).
Further, Powell et al. (2019) used the dust radial drift arguments to estimate
the gas mass in HD163296 at $M_{\rm disc}=0.2{\rm M_{\odot}}$. We therefore
chose $\Sigma_{0}$ such that the gas mass of our disc is $0.2{\rm M_{\odot}}$.
We fix the inner and outer radii of our computational domain at 5 and 200 AU,
respectively.

###  3.1 Release from an artificial steady state

Figure 1: Evolution of the planet-disc system when the planet is allowed to
evolve according to one of the following: accrete; migrate; both accrete and
migrate (from left to right). The vertical lines show the radial location of
the planet at the times indicated in the legends. The initial planet location
and mass were found by searching for a good fit to the values of the gap width
and the ring/gap contrast observed for the D86 gap in the disc of HD163296.
The experiments here indicate that when the planet is allowed to evolve as
expected it leaves this desirable mass-location parameter space worryingly
quickly, in $\sim 0.03$ Myr.

Since multi-dimensional numerical simulations of dusty discs with embedded
planets are expensive, it is often assumed that planet mass and orbit remain
constant during the simulations (e.g., Zhang et al., 2018). One then runs the
gas-dust simulations for a range of planet masses, $M_{\rm p}$, disc viscosity
parameter $\alpha_{\rm v}$, maximum grain size $a_{\rm max}$, etc., to try and
zero in on model parameters that best reproduce the morphology of the ALMA
dust continuum emission. In this section we follow a similar approach, placing
the planet at $86$ AU.

We start by setting the disc viscosity parameter $\alpha_{\rm v}=10^{-3}$;
this parameter will be varied in the sections to follow. To reach a quasi
steady state the simulations are run for 0.5 Myr. By performing a number of
numerical experiments we found that grains with the maximum size $a_{\rm
max}=1.5$ mm and planet mass $M_{\rm p}\approx 1{\,{\rm M}_{\rm J}}$ reproduce
the observed values for the gap width and contrast in HD163296 well. This
planet mass is reasonably close to those favoured by previous authors (e.g.,
Liu et al., 2019). Once we determined the appropriate $M_{\rm p}$, we restart
the simulation from the quasi steady state that the system reached at $t=0.5$
Myr, now allowing the planet to do one of the following: (i) accrete gas and
dust as per §2 but not migrate; (ii) migrate but not accrete gas; (iii) both
accrete and migrate.

Fig. 1 shows the resulting gas (top panels) and dust (bottom panels) surface
density profiles at several different times. The time in the legend is counted
from restart, so $t=0$ corresponds to the quasi steady state. When only
accretion is enabled (left panels) the planet evolves towards a gas gap
opening planet at its (forced) location. This evolution is rapid, both in
terms of planet mass and the dust surface density profile. The planet nearly
doubles its mass in 0.05 Myr. The lower left panel shows that dust gap becomes
an order of magnitude deeper within just 0.1 Myr. Within 0.5 Myr the planet
grows in mass to about $5{\,{\rm M}_{\rm J}}$ and opens a complete dust gap
from the star to $\sim 100$ AU. Even the gas accretion onto the star is
affected by that point. Clearly such evolution turns the modest gap in this
system into a transition like disc with a very large inner hole, somewhat like
PDS70, except the planet is much further out in HD163296.

When only planet migration is allowed (middle panels) we observe that the dust
gap parameters change less significantly but the position of the gap moves
towards the star extremely rapidly. The planet-star separation shrinks by a
factor of two in about 0.03 Myr. This is slightly longer than analytical
estimates for a planet migrating in the type I planet regime; the migration
time in our disc is given by

$t_{\rm mig1}=\frac{1}{2\gamma\Omega}\frac{M_{*}^{2}}{M_{\rm p}\Sigma
R^{2}}\left(\frac{H}{R}\right)^{2}$ (13)

where $\Omega_{\rm K}=(GM_{*}/R^{3})^{1/2}$ is the Keplerian angular frequency
at the planet location. The dimensionless factor $\gamma=2.5$ (Paardekooper et
al., 2010) is set by the assumed disc properties. Numerically,

$t_{\rm mig1}=0.016\hbox{ Myr}\;\frac{R}{100\text{ AU}}\left(\frac{M_{*}}{{\rm
M_{\odot}}}\right)^{1/2}\left(\frac{1{\,{\rm M}_{\rm J}}}{M_{\rm
p}}\right)\frac{0.1{\rm M_{\odot}}}{M_{\rm disc}}\;.$ (14)

For HD163296 and $R=86$ AU this yields $t_{\rm mig}\approx 0.01$ Myr. The fact
that our planet migrates a factor of a few slower is probably due to its
opening a moderately deep gap in the disc.

The middle panel in Fig. 2 also shows a change in the dust morphology from the
“ring behind the planet” to “the ring inside the planet” one. These changes in
the dust surface density profile are consistent with the results of Meru et
al. (2019). We shall discuss this further in §3.4.

Finally, when both accretion and migration are permitted the system evolution
is even more rapid. As the planet gains mass its inward migration accelerates
since it is migrating in the type I regime in which the migration time is
proportional to the planet mass. For the parameters of this system the planet
stops migrating rapidly only when it reaches $\sim 10$ AU as it enters the
type II regime.

These experiments demonstrate that allowing the planet to evolve as predicted
by the runaway CA scenario leads to rapid and very significant changes in the
planet parameters and the dusty disc morphology. It appears that an
evolutionary time scale for this disc-planet system is of order $0.05$ Myr,
which is surprisingly short for a $\sim 5$ Myr old system.

###  3.2 Release at the onset of runaway

Figure 2: Snapshots of the disc gas and dust profiles (top panels, similar to
Fig. 1), and the history of the planet-disc system ( middle and bottom row of
panels). The middle row left and right panels show planet mass $M_{\rm p}$ and
separation $a$ versus time, respectively, whereas the bottom row right panel
presents $M_{\rm p}$ vs $a$. The bottom row left panel shows the evolution of
two key gap parameters versus time, the gap width $\Delta_{\rm gap}$ and the
gap contrast $C_{\rm gap}$ (eqs. 15 and 16, respectively). The different
background colours in the middle and the bottom rows are to help the reader to
locate the three characteristic regimes in the disc evolution, which often but
not always have a relation to the state of the planet. In particular, in the
earliest pre-runaway phase the planet is too low mass to open a significant
gap, so the disc looks featureless (white background); in the runaway phase
the planet gains mass quickly and the gap in the disc is moderate (green
background); in the detached/end of growth phase the planet reaches its final
mass, opens a deep gap in both gas and dust and migrates inward slowly in the
Type II regime (purple background).

Figure 3: Same as Fig. 2 but for a lower viscosity parameter, $\alpha_{\rm
v}=10^{-4}$. Note that the planet now reaches the type II migration regime
earlier and then migrates slower, so it spends a longer time at large
separations. However, the dusty disc morphology in this case evolves into the
observationally undesirable configuration with the gap too deep and too wide,
even sooner. Despite an order of magnitude lower value of $\alpha_{\rm v}$
than in Fig. 2, the duration of the moderate gap phase is similarly short,
$\sim 0.12$ Myr.

Figure 4: Same as Fig. 2 and Fig. 3 but for $\alpha_{\rm v}=10^{-2}$. In this
high viscosity case planet migrates inward too rapidly and does not open a gap
until it is ”too late”, when it reaches the inner 10 AU disc. This is the
reason why there is no purple coloured background in this figure. Instead we
use the cyan background to emphasise the times when there is no planet in the
outer $R>10$ AU disc since it migrated inside 10 AU. The outer disc will then
look featureless in ALMA observations. The duration of the moderate gap case
is again short, $0.08$ Myr.

Experiments presented in §3.1 are initialised assuming a non-migrating and non
accreting planet. In a realistic situation the planet spends a long time in
the relatively long low mass pre-runaway phase before embarking on the rapid
mass gain (e.g., Pollack et al., 1996; Mordasini et al., 2012). In this
section we start instead with planets at the initial mass of $M_{\rm
p}=0.1{\,{\rm M}_{\rm J}}$ which is high enough for the planet to enter the
runaway phase, especially at large distances (e.g., Piso & Youdin, 2014; Piso
et al., 2015). We do not model how the planet grows to the initial mass as
this is model dependent (e.g., planetesimal vs pebble accretion, see
Emsenhuber et al., 2021a; Ndugu et al., 2019). This important but model
dependent issue is fortunately mainly irrelevant to our paper since very low
mass planets do not produce significant gaps in dusty discs. This is certainly
true for HD163296 as its star is one of the most massive and bright in Lodato
et al. (2019) sample. Furthermore, below when we discuss the rest of the
sample (§4) we exclude from it all low mass planet candidates that may not yet
be in the runaway regime.

Our approach in this section is to treat the location of the planet at the
time it reached the runaway mass of $0.1{\,{\rm M}_{\rm J}}$ as a free
parameter. We keep the same stellar parameters, the disc setup, and the
maximum grain size as in §3.1. By trial and error we found that planets
starting at $a=114$ AU generally reach the desired mass, $M_{\rm p}\sim
0.5-1{\,{\rm M}_{\rm J}}$, at around the gap location, and hence we present
here only such experiments. We now present the simulations for three values of
the disc viscosity parameter: $\alpha_{\rm v}=10^{-2},10^{-3}$ and $10^{-4}$.

In order to relate to the observed parameteters of the $R=86$ AU gap in
HD163296, we define the gap width, $\Delta_{\rm gap}$ as

$\Delta_{\rm gap}=\frac{\Delta R_{+}+\Delta R_{-}}{a}\;,$ (15)

where $\Delta R_{+}$ is the distance between the planet location ($a$) and the
centre of the “ring” behind the planet. The ring position is defined as the
farthest local maximum in $\Sigma_{\rm dust}$ (this definition avoids local
peaks within the dusty gap which may be present for low values of $\alpha_{\rm
v}$). $\Delta R_{-}$ is defined as the radial extent of the region inward of
the planet location where the dust surface density is 10 times lower than the
unperturbed dust density. We also define the gap contrast simply as the ratio

$C_{\rm gap}=\frac{\Sigma_{\rm min}}{\Sigma_{\rm max}}\;,$ (16)

where $\Sigma_{\rm min}$ and $\Sigma_{\rm max}$ are the dust minimum and
maximum surface densities over the gap region. $\Sigma_{\rm max}$ in practice
is usually found at the centre of the ring behind the planet (e.g., see the
bottom left panel of fig. 1) or at the planet location for a rapidly migrating
planet (e.g., the cyan curve in the top right panel of fig. 2). This
definition is qualitatively similar to that used by Rosotti et al. (2016) if
the disc is optically thin and dust temperature does not vary appreciably over
the gap region. We note that there is a degree of arbitrariness in these
definitions, however they are conservative, and we found them to be more
robust when the planet is allowed to migrate than definitions introduced by
other authors in the past for static planets.

####  3.2.1 Viscosity parameter $\alpha_{\rm v}=10^{-3}$.

Fig. 2 presents a calculation in which $\alpha_{\rm v}=10^{-3}$. The
simulation was terminated at $t\approx 0.85$ Myr when the planet reached 10
AU, twice the inner boundary of our computational domain. The top row of
panels show several snapshots of the disc gas and dust surface density
profiles, similar to fig. 1. The middle and the bottom rows of panels detail
the system evolution. The middle left panel shows planet mass and Crida et al.
(2006) parameter $C_{\rm p}$, while the middle right shows planet separation,
$a$, versus time. The bottom row shows the dust gap width and contrast versus
time (left panel) and the planet mass versus position (right panel). The two
bottom rows are shaded with three different colours to emphasise the three
distinct states of the planet-disc system.

Further, to expose the pace of planet evolution at different times, the planet
mass-position diagram (bottom right panel) is sampled every 0.01 Myr. We see
that at the beginning of the run, the planet gains mass relatively slowly as
its “feeding zone” is small, and its migration is also relatively slow.
However, this accelerates strongly when the planet grows in mass by a factor
of a few. The planet crosses the radius of 86 AU (the centroid of the dusty
gap we study here) at $t=0.104$ Myr. The planet mass at that moment is $M_{\rm
p}=0.67{\,{\rm M}_{\rm J}}$. While this mass is smaller than $1{\,{\rm M}_{\rm
J}}$ we favoured in §3.1, we note that both the gap width and the gap contrast
(the bottom left panel) are not far off from those inferred from observations
of this gap ($\Delta_{\rm gap}=0.17$ and $C_{\rm gap}=0.15$; Zhang et al.,
2018). We place the green open circle on Fig. 2 to mark that moment.

Prior to $t\approx 0.07$ Myr the gap depth is too small (that is, the gap
contrast is $C_{\rm gap}\gtrsim 0.8$) to be compatible with the observations.
Since such a gentle gap may be unresolved (missed) in a disc less bright and
not so well studied as HD163296, we term this early time period “featureless”;
it has a white background in the figure. During the time marked with the green
background the parameters of the dusty gap are moderate, and are roughly
consistent with that seen in HD163296. After $t\approx 0.15$ Myr the planet
opens a gap too deep to be compatible with observations. During the phase that
comes after $t\approx 0.15$ Myr, coloured with purple background, the disc
morphology is that of a wide dust and gas gap transitioning into a complete
inner hole in dust (cf. the green dotted curve in the top right panel).

Note that the green region on the plot corresponds approximately to the
runaway phase of the planet growth, whereas the wide gap/inner hole disc state
corresponds quite well to the detached phase when the planet stops growing.
This could also be rephrased in terms of the two planet migration regimes. The
wide gap/inner hole disc state corresponds to the planets migrating in type II
regime, when $C_{\rm p}\leq 1$ (compare the left panels in the middle and the
bottom rows). This is unsurprising since planets open deep gaps in the gas
disc when $C_{\rm p}\leq 1$ (Crida et al., 2006); the gaps and holes in the
dust discs are always stronger than those in the gas discs (e.g., Dipierro et
al., 2016; Dipierro & Laibe, 2017). This correspondence between the dusty disc
morphology, on the one hand, and the planetary growth and migration phases, on
the other, is exciting: ALMA may be able to probe planetary growth and
migration in the CA scenario directly by the statistics of gaps/holes in
protoplanetary discs.

Quantitatively, the main conclusion from this numerical experiment confirms
what we found from the steady state start models in §3.1: the planet-disc
system spends a surprisingly short period of time in the corner of the
parameter space consistent with observations. The duration of this phase is no
longer than 0.07 Myr, which is comparable with the results we obtained from
the steady state start models.

####  3.2.2 Viscosity parameter $\alpha_{\rm v}=10^{-4}$.

Fig. 3 shows a numerical experiment fully analogous to that presented in Fig.
2, §3.2.1, but with a much lower viscosity parameter, $\alpha_{\rm
v}=10^{-4}$. The initial evolution of the planet is very similar to that in
fig. 2. However, due to lower viscosity, a gap in the disc is opened earlier,
and the planet switches to type II migration upon reaching the radius of 50
AU, rather than $\sim 30$ AU. The effects of the planet on the dust also
become noticeable much earlier, so that the featureless disc phase terminates
even sooner.

The largest difference with the $\alpha_{\rm v}=10^{-3}$ case is however the
much slower rate of planet migration in the Type II regime. Indeed, during the
1 Myr of this simulation the planet migrates only to about 33 AU. This has
major observational implications for the appearance of the dusty gap. The gap
in this case is situated farther out, and is much wider and deeper. Once
again, the gap contrast is formally zero after $t\sim 0.15$ Myr, which clearly
rules this scenario out as a reasonable model for the gap we study here.

####  3.2.3 Viscosity parameter $\alpha_{\rm v}=10^{-2}$.

Fig. 4 shows yet again the same numerical experiment but with a higher value
of the viscosity parameter, $\alpha_{\rm v}=10^{-2}$. In this case we observe
major qualitative differences in the disc-planet system evolution. The planet
does not manage to open a deep gap in the disc ($C_{\rm p}>1$) until it
reaches the inner boundary of our computational domain, exiting the grid.
Except for a very brief time, the gap contrast and the gap width are never
large enough to classify this system as a wide gap system. This is why there
is no region with purple background in Fig. 4. After the planet exits the
computational grid through the inner boundary at 5 AU, the outer disc looks
featureless. We use cyan background to highlight this state of the disc-planet
system. Note that our results would not be much different even if we extended
the grid to much smaller radii. Indeed, the type II migration in this case is
quite rapid and so the planet would migrate ”into the star” quickly. Even if
the planet were to hang about in the inner disc for an unspecified reason its
torques would not affect the disc at radii that ALMA can resolve in a typical
source. Finally, molecular line observations constrain $\alpha_{\rm v}\lesssim
3\times 10^{-3}$ (Flaherty et al., 2017). Concluding, this high viscosity
parameter experiment is also not promising as an explanation for the ALMA gap
at 86 AU in HD163296.

###  3.3 Observational implications for the Fiducial model of HD163296 gap at
86 AU

None of the three viscosity values considered above yielded long-lived
moderate gap discs, although for different reasons. To summarise, for the
lowest viscosity parameter, the planet runs away in mass to become a massive
gas giant at wide separations (tens of AU). It opens a deep gap/inner hole and
migrates inward very slowly in the type II regime. At the highest $\alpha_{\rm
v}$ explored, the planet does not fully reach its gas gap opening mass before
it reaches the inner $\sim 10$ AU disc, so it has a moderate gap while at wide
separations. However, the planet migrates from 100 AU to 10 AU in $0.08$ Myr.
The intermediate viscosity case, $\alpha_{\rm v}=10^{-3}$, shows elements of
both of these extreme cases, but the end result is nonetheless the same: the
moderate gap state lasts time $\Delta t_{\rm mod}\sim 0.08$ Myr.

A corollary of this rapid planet evolution is the fact that after just a brief
period of time the disc in HD163296 would look very different from the
observed one. Fig. 5 shows the gas and dust surface density profiles at
$t=0.2$ Myr for these three runs. At the lowest viscosity parameter, the gap
width and depth are too large compared with the observed values; at the
highest $\alpha_{\rm v}$ the planet is lost into the inner disc and so the
outer disc is featureless in both gas and dust; and at $\alpha_{\rm
v}=10^{-3}$ the planet both grew and migrated in too much. Furthermore, these
conclusions are robust to changes in the assumed temperature profile as shown
in Appendix A.

Figure 5: The gas and dust surface density profiles at time $t=0.2$ Myr for
the fiducial simulations with the three values of the viscosity parameter.
This shows that only $t=0.2$ Myr after the beginning of the runaway accretion
the dusty disc morphology is not right compared with the observed one.

The problem of too rapid planet evolution can also be understood via the
following argument. Let us consider an ensemble of protostars with stellar and
disc parameters as we assumed for HD163296, and let planets grow in these
discs according to the CA runaway scenario. We assume that the ages of these
systems is a uniformly distributed variable between 0 and $\sim 2t_{0}\sim 10$
Myr. If each system hatches one embryo massive enough to get into the runaway
regime, then we have a chance $\Delta t_{\rm mod}/2t_{*}\sim 0.01$ to catch
any given system in the moderate gap state. If we pick one system from this
ensemble in random, it would be very surprising if it were a disc with a
moderate gap. Further, if we observed a number of randomly selected objects in
this ensemble then we should expect that most of them would be in either
featureless (white background in Figs. 2, 3, and the cyan one in 4) or the
deep gap/inner hole, the low viscosity detached planet case (the purple
background in the figures). The featureless state corresponds to the times
before a massive core has grown and/or the planet migrated inside 10 AU. The
wide gap/inner hole state is realised at low viscosity parameter values when
the planet enters the type II regime.

This predicts a lack of systems with moderate gaps, and abundance of
featureless and wide gap discs. These theoretical predictions are rather
divergent from what we see for the observed sample of ALMA discs/planets
(Lodato et al., 2019).

###  3.4 Suppressing gas accretion and migration

Nayakshin et al. (2019) showed that decreasing the gas accretion rate onto
planets while they are on the runaway track of CA by $F_{\rm acc}\sim$ an
order of magnitude brings CA theory predictions much more in line with the
mass function of candidate ALMA planets. This finding is also consistent with
a number of recent simulations of gas accretion onto CA planets which suggests
that gas accretion may be much less efficient in the runaway phase than
previously believed (e.g., Szulágyi et al., 2014; Ormel et al., 2015; Fung &
Chiang, 2016; Lambrechts & Lega, 2017).

We therefore investigate how suppressing gas accretion by a factor $F_{\rm
acc}=10$ would affect our results. Nayakshin et al. (2019) held their planets
at fixed orbits artificially whereas our model includes planet migration.
Experiments presented in Fig. 1 demonstrated that the planets held at a
constant mass (which is equivalent to $F_{\rm acc}=\infty$) are still
challenged by the ALMA data as planet migration may also be too rapid.
Therefore we also introduce a factor $F_{\rm mig}\geq 1$ to reduce the rate of
planet migration. Similar approach was chosen in the past by many authors,
e.g., Ida & Lin (2004b).

Fig. 6 compares the gas and dust surface density profiles of the disc in the
fiducial $\alpha_{\rm v}=10^{-3}$ simulation presented in Fig. 2 with that
performed with slow-down factors $F_{\rm acc}=10$ and $F_{\rm mig}=1$ (slow
accretion but nominal migration), and $F_{\rm acc}=F_{\rm mig}=10$ (both
accretion and migration suppressed). As in the fiducial case, the planet
starts with mass $M_{\rm p}=0.1{\,{\rm M}_{\rm J}}$ and the initial separation
of $a_{\rm p}=114$ AU. In Fig. 6 we compare the disc surface density profiles
when the planet positions coincide, at $a\approx 60$ AU. This occurs at
different times in the simulations, and at different planet masses. We note
that our main conclusions do not depend on the planet position chosen for this
illustrative comparison.

Figure 6: Comparison of the disc surface density profiles of the nominally
growing planet, previously shown in fig. 2, with that of suppressed evolution
simulations. It is now plainly obvious that the dust surface density in the
fiducial simulation (black dot-dashed curves) is not an appropriate model to
account for ALMA observations because the bright ring is in front of the
planet, not behind. If only planet accretion is suppressed (cyan) then the
planet never grows massive enough to affect the disc significantly enough
before it migrates into the inner disc. The dust profile with both accretion
and migration suppressed, on the other hand, is qualitatively reasonable, and
is much longer lasting.

In the suppressed accretion but fiducial migration scenario (cyan curves in
the figure) the planet migrates inward so rapidly that it does not manage to
grow enough to affect the disc sufficiently. Its mass is only $0.2{\,{\rm
M}_{\rm J}}$ and with the chosen $\alpha_{\rm v}$; this is plainly too low to
open a detectable gas. This scenario looses planets into the inner disc too
rapidly to be compatible with observations of HD163296.

In the case when both accretion and migration are suppressed (red solid
curves), we see that the planet is much more massive by the time it reaches
$a=60$ AU. The planet mass is lower in the suppressed simulation than in the
fiducial simulation, but the gap opened by the planet is deeper in both gas
and dust. Additionally, in the fiducial simulation, there is a significant
radial offset between the positions of the gap and the planet. This is because
a very rapidly migrating planet may be said to push the material in front of
it, sort of snow ploughing it, creating the dusty ring in front of itself
rather than behind it (Meru et al., 2019). For the ring to be behind the
planet, the planet cannot migrate inward much more rapidly than the dust
drifts in. HD163295 is $\approx 5$ Myr old. Therefore a planet migrating on a
time scale of $\sim 0.1$ Myr could not possibly be providing a barrier to hold
the dust back. If this were the rule for protoplanetary discs then it would be
hard to invoke planets as the mechanism for slowing down dust radial drift (as
argued in, e.g., Pinilla et al., 2012), leaving additionally the conundrum of
dust survival in old discs unsolved (e.g., Birnstiel et al., 2009, 2012).

We also note the general success of the stationary (thus $F{\rm mig}\gg 1$)
planet models employed by Zhang et al. (2018) to explain the DSHARP sample of
ALMA gaps/rings. If the planets migrated rapidly in as predicted by the
standard theory, then the planetary scenario for the origin of gaps would not
be so convincing. Slowly migrating planets, such as the simulation with
$F_{\rm mig}=10$ in fig. 6 are clearly more promising in addressing this
issue.

The fiducial model also probably contradicts observations of the gas
kinematics in HD163296 as well. Pinte et al. (2020) show that the localized
velocity perturbation (kink) found in HD163296 near the gap at 86 AU falls
right into the middle of the gap after the deprojection of the CO emission
map. That is, the radial location of the planet coincides with the centroid of
the dusty gap. The fiducial model, however, shows a significant offset between
the location of the planet and the dusty gap position. The suppressed
migration and accretion model does not have such an offset and are thus more
promising.

Figure 7: Same as figs. 2 to 4 but for the simulations in which both migration
and accretion are suppressed by one order of magnitude, $F_{\rm acc}=F_{\rm
acc}=10$ (see §3.4) and for three different values of the disc viscosity
parameter $\alpha_{\rm v}$ as listed in the legend in the top right panel. The
times spent by the disc in the moderate gap configuration are shown in the
bottom right panel in the legend. Note that these are much longer than those
of fiducial simulations shown in Figs. 2 – 4.

Fig. 7 shows the time evolution of the planet-disc system for the suppressed
planet migration and accretion scenario, $F_{\rm acc}=F_{\rm mig}=10$, now for
three values of the viscosity parameter: $\alpha_{\rm v}=10^{-4}$, $10^{-3}$
and $10^{-2}$. The meaning of all of the curves is the same as in fig. 2.
While qualitatively this evolution may look similar to that shown in figs. 2
to 4, the discs spend much more time in the state of moderate gaps, here
defined as $C_{\rm gap}\leq 0.9$ and $\Delta_{\rm gap}\leq 0.4$. These times
are listed in the caption in the bottom right panel, and vary between 0.6 and
1 Myr. This is much more reasonable. If a few massive planet embryos formed in
the disc of HD163296 in the last few Myr, then we would have a decent
statistical chance to detect at least one of them in the state of moderate gap
now.

##  4 Results for the ALMA gap sample

We now use the sample of ALMA gaps/candidate planets investigated in Nayakshin
et al. (2019) (the wast majority of which is from Lodato et al., 2019) in
order to see how the lifetime of the observed gaps could constrain the CA
scenario for planet formation. To do this, we first remove from the sample any
gaps that are believed to be opened by planets less massive than the mass for
which the planet should be in the runaway CA regime, which we define here as
$0.1{\,{\rm M}_{\rm J}}$. We also chose only the gaps located beyond 10 AU. We
have 21 gaps remaining in the sample. The parameters of these systems are
presented in Table 1.

We then follow the “release from the steady state” approach from §3.1 rather
than the “release at the runaway” approach from §3.2. In the latter, the
initial location of the embryo is unknown and thus one needs to make guesses
and iterate to constrain it. The number of runs required to match our sample
of ALMA gaps and to then also study the implications for the CA scenario is
too large for such trial-and-error investigation. On the other hand, in §3 we
found that the simpler “release from the steady state” approach yields
quantitatively similar value for the duration that the gaps may remain in the
observed state.

For each object in the sample we set up a steady-state power-law disc model
with the same radial extent as before (5 to 200 AU). We determine the disc
mass via the following constraints. First of all, given the observed accretion
rate onto the star, $\dot{M}_{*}$, and its age, $t_{*}$, we require that there
must be at least $M_{\rm acc}=2\dot{M}_{*}t_{*}$ of gas to maintain this
accretion rate (e.g., Jones et al., 2012). Note that this minimum disc mass
estimate simply states that disc is the source of gas accretion onto the star
and must have fed the star for time $t=t_{*}$. The factor of 2 in the estimate
can be larger (up to $\sim 10$) and accounts for the fact that matter also
spreads outward (Jones et al., 2012). This estimate is conservative and is
likely to hold no matter the detail of the exact mass transfer mechanism from
the disc to the star. For example, the recent MHD disc wind formalism
developed by Tabone et al. (2021) will increase the disc mass estimate because
only a fraction of the disc mass reaches the star in this case.

Next, following the same arguments as Powell et al. (2017, 2019) we require
the drift timescale for 1 mm sized grains at the location of the planet to be
at least as long as $t_{*}$. While planets are often invoked as a solution to
the short dust radial drift times conundrum (e.g., Pinilla et al., 2012), most
of the observed ALMA gaps are quite shallow and hence may at best slow down
but not prevent the drift. We term this minimum disc mass estimate as $M_{\rm
drift}$.

We then take the disc mass to be a maximum of $M_{\rm acc}$ and $M_{\rm
drift}$. However, Hall et al. (2019) show that discs that exceed mass $M_{\rm
sg}=0.25M_{*}$ become very strongly self-gravitating. Spiral density arms in
such discs are likely to be detectable with ALMA, and transfer the mass onto
the star at very high rates, depleting the discs of mass in just $\sim 10^{4}$
years. Both of these effects strongly contradict what is observed for the ALMA
sample: the discs have annular symmetry and accretion rates onto the star some
$\sim 2-4$ orders of magnitude smaller. We therefore limit the disc mass by
$M_{\rm max}=0.2M_{*}$.

The disc surface temperature profile is set in this modelling by

$T(R)=\left[\phi\frac{L_{*}}{4\pi\sigma_{B}R^{2}}\right]^{1/4}\;,$ (17)

where $\phi=0.1$ and $L_{*}$ is the luminosity of the star (this model is
similar to that used by, e.g., Long et al., 2018). We set the disc viscosity
parameter to $\alpha_{\rm v}=10^{-3}$ in this section.

We then insert a planet with the mass inferred from the gap properties at the
location of the observed gap into the unperturbed disc. Planet growth and
migration are not allowed for 0.5 Myr to allow the disc-planet system to reach
a steady state. We then allow the planet to gain mass by accretion and
migrate. We define the planet migration and planet accretion time scales as

$t_{\rm mig}=-\frac{a}{\dot{a}}\;,$ (18)

and

$t_{\rm acc}=\frac{M_{\rm p}}{dM_{\rm p}/dt}$ (19)

respectively. Given the close association between the gap properties and the
planet mass and separation we define the gap lifetime as the following:

$t_{\rm gap}^{-1}=t_{\rm acc}^{-1}+t_{\rm mig}^{-1}$ (20)

Note that this definition is a smoothed out alternative to $t_{\rm
gap}=\min[t_{\rm acc},t_{\rm mig}]$.

Figure 8: Distribution of the observed gap lifetimes from the ALMA sample for
fiducial (top panel), accretion suppressed (middle panel), and both accretion
and migration suppressed scenarios (bottom). Only the last of the three cases
provides reasonably long gap lifetimes.

The top panel in Fig. 8 shows the histogram of $t_{\rm gap}$ evaluated with
the fiducial CA model for our 21 gaps in the sample. Table 1 shows $t_{\rm
gap}$ for each individual gap/planet. We observe in Fig. 8 that for most of
the gaps $t_{\rm gap}<0.1$ Myr, although there are some systems with longer
$t_{\rm gap}$. The mean value of $t_{\rm gap}$ for the sample is $0.15$ Myr.
These time scales are worryingly short compared with the system ages, $t_{*}$,
which are in the range of 1-10 Myr. More specifically, the mean of the ratio
$t_{\rm gap}/t_{*}$ is $0.06$ for the sample.

In the middle panel of Fig. 8 we test the suppressed gas accretion scenario,
setting $F_{\rm acc}=10$, but allowing planets to migrate at the fiducial rate
($F_{\rm mig}=1$). This extends the gap lifetime somewhat, especially for the
shortest $t_{\rm gap}$ bins. However, $t_{\rm gap}$ is still unacceptably
short on average.

Finally, in the bottom panel we investigate the scenario in which both gas
accretion and planet migration are suppressed by a factor of 10. In this case,
the typical gap lifetime is $\sim 0.5$ Myr, and the mean for the ratio $t_{\rm
gap}/t_{*}=0.27$. These results are much more reasonable, requiring a few gas
giant planets to be born per system in this sample.

Figure 9: Number of planets per system required for a probability of order
unity for ALMA to detect a planet in that system (cf. eq. 21).

In Fig. 9 we show the number of planets that has to be born in the system to
yield a probability of order unity for us to observe it,

$N_{\rm req}\sim\frac{t_{*}}{t_{\rm gap}}\;.$ (21)

The red open histogram shows this quantity for the fiducial model, that is,
for $F_{\rm acc}=F_{\rm mig}=1$. In this scenario, for gaps opened by planets
to be as common as implied by the observations we would need tens to a hundred
of gas giant planets to be born per star, which is clearly excessive. For
slower gas accretion onto planets, $F_{\rm acc}=10$, we obtain the green
shaded histogram in Fig. 9. This however relaxes the constraint only by a
factor of $\sim 2$. Finally, the model where both accretion and migration are
suppressed by an order of magnitude (the blue filled histogram) needs $\sim
10$ gas giant planets per star, although this is a significant dispersion
around this number, with many systems requiring just a few planets whereas for
some a few tens is needed.

The model with suppressed planet evolution, $F_{\rm acc}=F_{\rm mig}=10$ is
clearly favoured compared with the predictions of the other two scenarios.
While $N_{\rm req}$ is still somewhat uncomfortably large for a few systems in
the suppressed scenario, it is not yet clear if this constitutes a problem. It
is possible that future better observations may reduce $N_{\rm req}$. We also
need to remember that the present ALMA sample of resolved discs is dominated
by the brightest systems (Andrews et al., 2018; Long et al., 2018), and hence
they may be especially effective in spawning gas giant planets.

##  5 Discussion

###  5.1 Main results

In this paper we studied 1D models of dusty protoplanetary discs with embedded
gas giant planets growing according to the Core Accretion scenario. We found a
kinship between planet growth phases and the dusty disc morphology. As is well
known (e.g., Pollack et al., 1996), on the way to becoming gas giant planets,
CA planets go through the three main phases: (i) a low mass planetary core;
(ii) a moderate mass planet in the runaway gas accretion phase; (iii) a
massive gas giant planet that stopped growing by opening a deep gas gap in the
disc. We find that in the pre-runaway phase (i), the planet is too low mass to
open a detectable gap even in the dust disc. Barring other mechanisms creating
dust rings/gaps in discs, phase (i) hence corresponds to featureless, smooth,
protoplanetary discs. In the phase (ii) the planet mass is $\sim
0.1-0.5{\,{\rm M}_{\rm J}}$. Such a planet is massive enough to open gaps in
the dust but not gas (e.g., Dipierro et al., 2016; Rosotti et al., 2016). This
phase is probably the one most frequently seen in ALMA observations of bright
discs (e.g., Huang et al., 2018; Zhang et al., 2018; Long et al., 2018).
Finally, when the planet opens a deep gap in the gas disc, there is a very
wide and deep gap in the dust disc, evolving into an inner disc hole if the
mm-sized dust in the inner disc region is completely cutoff from the outer
disc (Rice et al., 2006).

Significantly, the theory predicts that phases (i) and (iii) are long lived,
whereas phase (ii) is very short. As expected based on classical works (e.g.,
Pollack et al., 1996; Ida & Lin, 2004a), we found that gas giant planets grow
through the mass range $M_{\rm p}\sim 0.1{\,{\rm M}_{\rm J}}$ to a (gas) gap
opening mass in a matter of $\lesssim 0.1$ Myr. This implies that in the CA
scenario the featureless discs, and the discs with very deep and wide holes,
should be the mode of protoplanetary discs, whereas the moderate dust-only
gaps should be a minority of discs transiting from phase (i) to phase (iii).
We compared this prediction with the statistics of the observed annular
structures in the samples of resolved protoplanetary discs (Andrews et al.,
2018; Long et al., 2018; Lodato et al., 2019) and uncovered a significant
mismatch between theoretical predictions and observations. On the quantitative
level, the gap opening mass is a function of the poorly constrained disc
viscosity parameter $\alpha_{\rm v}$; however, the predictions of the theory
disagree with observations for any $\alpha_{\rm v}$.

At low values of $\alpha_{\rm v}\ll 10^{-3}$, planets open deep gaps in the
disc quickly. Their growth therefore terminates at lower masses, typically
$M_{\rm p}\lesssim 1{\,{\rm M}_{\rm J}}$. They open a deep gap in the gas
disc, and migrate inward very slowly in the type II regime. This slow
migration could explain the surprisingly large abundance of massive gas giant
planets observed by ALMA; however, the morphology of dusty discs in this case
is deviant from what is observed. Planets that open deep gaps in gas discs
filter out large dust particles very efficiently (e.g., Rice et al., 2006),
leading to very wide gaps turning into complete inner holes in mm-sized dust
discs (see fig. 3). If this was the case then we would expect most ALMA discs
look like wide gap transition discs, such as the famous PDS70 system where two
massive gas giants opened an inner dust hole that extends from $\sim 2$ AU
from the star to $\sim 30$ AU (e.g., Haffert et al., 2019; Keppler et al.,
2019; Benisty et al., 2021). However, this system is unique amongst the ALMA
systems in many ways, with most other discs showing the signs of mild gaps.

At high values of the disc viscosity parameter, $\alpha_{\rm v}\gg 10^{-3}$,
on the other hand, planets are not able to open deep gaps in the gas disc
until they reach higher masses, $M_{\rm p}\sim$ a few ${\,{\rm M}_{\rm J}}$
typically. They therefore migrate in the type I regime for longer, and are
thus less likely to create discs with wide inner dust hole morphology.
However, these planets migrate into the inner disc too rapidly, making it hard
to understand the ubiquity of moderate gaps at large separations observed by
ALMA (unless they are not opened by planets).

At intermediate values of viscosity parameter, $\alpha_{\rm v}\sim 10^{-3}$,
models experience a mix of the two difficulties explained just above.
Additionally, we found that for any value of disc viscosity the phase at which
the planets produce moderate gaps in the disc, $t_{\rm gap}$, is too short,
$t_{\rm gap}\lesssim 0.1$ Myr (the phases shaded green in figs. 2 to 4). We
can estimate the probability of observing a phase of duration $t_{\rm gap}$ in
a system with age $t_{*}$ as $t_{\rm gap}/t_{*}$. Given that most bright ALMA
discs show dusty gaps, we require $N_{\rm req}\sim t_{*}/t_{\rm gap}$ gas
giant planets born per ALMA disc. This results in $N_{\rm req}$ in many tens
to a few hundred (cf. the red open histogram in fig. 9), which is clearly
excessive and would violate the available solid and (even) gas mass budget of
protoplanetary discs.

In §3.4 we therefore explored the “suppressed” planet evolution models. Early
one dimensional planet envelope contraction models (e.g., Mizuno, 1980;
Stevenson, 1982; Ikoma et al., 2000) arrived at very high gas accretion rates
in the runaway regime. However, 3D simulations by various authors in the
recent past (e.g., Szulágyi et al., 2014; Ormel et al., 2015; Lambrechts &
Lega, 2017) showed that gas accretion onto planets in the runaway regime may
be far less efficient. The key reason for this are complicated 3D gas flows
that carry gas both in and out of the Hill sphere of the planet, resulting in
a much smaller net planet accretion rate for Saturn-mass planets (e.g.,
Lambrechts et al., 2019). Nayakshin et al. (2019) found in simplified (gas
only, and no planet migration) models that suppressing gas accretion by a
factor $F_{\rm acc}\sim 10$ was required to account for the predominance of
sub-Jovian mass planets inferred to be present in ALMA discs (Zhang et al.,
2018; Lodato et al., 2009).

Our experiments here show that this conclusion also holds in models with
planet migration and dust dynamics included. Planets with gas accretion
suppressed by factor $F_{\rm acc}=10$ allow planets to spend more time in the
intermediate mass regime. Gaps opened by such planets in dusty discs are more
moderate (cf. the green dashed curve in fig. 6). However, at the full
unquenched migration speeds these planets disappear into the inner disc too
quickly, and the lifetime of the moderate gap phase, $t_{\rm gap}$, is
extended only by a factor of $\sim 2$ (cf. the middle panel in fig. 8).

We therefore also explored the models in which both gas accretion and planet
migration is suppressed. This provides far more promising results. First of
all, for slowly migrating planets the dust ring morphology switches from the
“inside ring” to the “outside ring” (see Meru et al., 2019, and fig. 6)
morphology. The latter is generally invoked as explanation for the observed
ALMA dust structures when one neglects planet migration (e.g., Dullemond et
al., 2018; Zhang et al., 2018). The models with unsuppressed migration result
instead in the “inside ring” dust morphology, which, in addition to all the
other problems already pointed out in this section, has so far not been
invoked for any of the observed ALMA protoplanetary discs. Secondly, models
with both growth and migration suppressed produce far longer lived moderate
gap structures (the bottom panel in fig. 8) that require a far smaller number
of gas giants born per ALMA disc (the blue filled histogram in fig. 9).

###  5.2 Caveats

Here we made an assumption that the annular features seen in the dusty
protoplanetary discs with ALMA are all due to embedded planets. This
assumption is clearly an oversimplification since other physical processes
were also shown to result in annular features, such as dead zones (Ruge et
al., 2016), condensation fronts for important chemical species (Zhang et al.,
2015), dust self-induced pile-ups (Gonzalez et al., 2015), effects of vortices
(e.g., Barge et al., 2017), etc. We however believe that our conclusions are
robust to at least a qualitative if not a quantitative level. We would expect
a major change to our conclusions only if the majority of gaps we included in
our sample are not due to planets. While this is in principle possible,
presence of gas giant planets in protoplanetary discs has now been confirmed
through independent means, such as localised velocity perturbations to gas
flows in vicinity of planets (e.g., Perez et al., 2015; Pérez et al., 2018;
Pinte et al., 2018; Casassus & Pérez, 2019; Teague et al., 2019; Pinte et al.,
2020), H-$\alpha$ imaging (Haffert et al., 2019), and even detection of
circum-planetary discs (Benisty et al., 2021). Nevertheless, more
observational and theoretical/simulation work is needed to disentangle
different mechanisms for formation of annular disc features.

In §3.4 we made the assumption that planet migration can be slowed down in the
type I regime with regards to the well known formulae from Paardekooper et al.
(2011). This allowed us to reduce the mismatch between the CA scenario
predictions and ALMA observations. It is not clear whether this is reasonable.
Nelson & Papaloizou (2004) showed that MHD turbulence may result in far slower
although chaotic type I migration of planets. However the values of the
turbulent viscosity parameter they considered were larger than the upper
limits on $\alpha_{\rm v}$ for the ALMA discs in which such measurements were
possible.

There is some recent work showing that luminous planets can migrate slower or
even outward due to “thermal torques” (Benítez-Llambay et al., 2015; Masset,
2017), but it appears that this effect is short lived (Guilera et al., 2021)
and is unlikely to reduce planet inward migration on time scales of $\gtrsim
1$ Myr; it also appears to affect planets only inside $10-20$ AU. Kimmig et
al. (2020) have demonstrated that planets migrating in the type II regime may
be driven outward if the main mechanism of angular momentum transfer is due to
MHD disc winds. However, we argued that we need type I migrating planets to
account for the majority of ALMA annular features as type II migrating planets
produce gaps that are too deep and wide.

The last point (on the potential importance of MHD driven disc winds, see,
e.g., Tabone et al., 2021) emphasises an important limitation of our study.
Here we assumed that we know the protoplanetary disc physics and thus we can
use ALMA observations to constrain CA scenario for planet formation. To be
specific, we assumed here that the angular momentum transport is viscous and
the disc is in a quasi steady state, which implies a minimum disc mass as
explained in §4. There can obviously be other complications such as non-steady
accretion. It is thus not impossible that the disc structure or its
parameters, such as its mass, is significantly off from what we currently
think it is. If that is true then perhaps the data can be made consistent with
the classical runaway accretion scenario of CA and instead require a
significant re-evaluation of our views of how the discs evolve.

##  6 Conclusions

In this paper we have shown that ALMA and other modern data of protoplanetary
discs with embedded planet candidates present us with a unique opportunity to
test planet formation models, supplementing constraints from the previous
exoplanetary data of older discless, planetary systems. We constrained
protoplanetary disc masses given their age, the presence of mm-sized grains,
and the observed gas accretion rates onto the star. This then allows us to
evaluate how planets embedded in these discs should evolve for a specific
planet formation model. As an example of this method, we examined the very
successful Bern population synthesis implementation (Emsenhuber et al., 2021a,
b) of the Core Accretion scenario. Our results suggest that if the observed
substructure is due to embedded planets, then both accretion and migration
will probably need to be substantially slower than expected. We also emphasise
that while we focused on the Bern model, our results would apply to any planet
formation model that incorporates planet accretion and migration at rates
similar to those in the Bern model. Finally, as noted in §5.2, another
possible interpretation of our results is that the data may require re-
evaluation of our disc models rather than planet formation physics. In any
event, we believe that our results indicate that the future of planet
formation population synthesis is in applying it to data of both ALMA
protoplanetary discs and the much older exoplanetary systems that have lost
all signatures of their primordial discs.

Table 1:  The sample of ALMA discs and candidate planets investigated in the paper, and the inferred gap lifetime $t_{\rm gap}$ in units of Myr for the three scenarios. See §4 for detail. Except for the total disc mass, $M_{\rm disc}$, all the entries to the left of the vertical bar in the Table are those derived in the previous literature through observations or simulations to fit dusty disc morphology (e.g., Zhang et al., 2018; Lodato et al., 2019). To the right of the vertical bar, the three columns list the results of our modelling, that is, the gap lifetimes for the Fiducial model, the one with gas accretion suppressed, and with both accretion and migration suppressed by factor of 10, respectively. Object | $M_{*}$ | $L_{*}$ | Age | $M_{\rm disc}$ | $R_{\rm gap}$ | $\Delta_{\rm gap}$ | $C_{\rm gap}$ | $M_{\rm p}$ | $t_{\rm gap}$ | $t_{\rm gap}$ | $t_{\rm gap}$
---|---|---|---|---|---|---|---|---|---|---|---
| [$M_{\odot}$] | [$L_{\odot}$] | [Myr] | [$M_{\odot}$] | [AU] | [AU] |  | [${\,{\rm M}_{\rm J}}$] | Fiducial | $F_{\rm a}=10,F_{\rm m}=1$ | $F_{\rm a}=F_{\rm m}=10$
DS TAU | 0.58 | 0.25 | 3.981 | 0.116 | 33 | 0.821 | 0.042 | 6.061 | 0.540 | 0.540 | 0.540
MWC480 | 1.9 | 17.4 | 6.310 | 0.382 | 73 | 0.454 | 0.013 | 2.300 | 0.166 | 0.194 | 0.241
GO TAU | 0.36 | 0.2 | 5.623 | 0.072 | 59 | 0.389 | 0.056 | 0.400 | 0.065 | 0.107 | 0.361
DL TAU | 0.98 | 0.65 | 1.000 | 0.059 | 67 | 0.207 | 0.157 | 0.163 | 0.032 | 0.225 | 1.035
DL TAU | 0.98 | 0.65 | 1.000 | 0.059 | 89 | 0.292 | 0.470 | 0.459 | 0.051 | 0.115 | 1.001
CI TAU | 0.89 | 0.8 | 1.585 | 0.104 | 14 | 0.636 | 0.455 | 0.750 | 0.239 | 0.239 | 0.239
CI TAU | 0.89 | 0.8 | 1.585 | 0.104 | 48 | 0.225 | 0.820 | 0.150 | 0.022 | 0.094 | 0.536
CI TAU | 0.89 | 0.8 | 1.585 | 0.104 | 119 | 0.186 | 0.600 | 0.400 | 0.028 | 0.077 | 0.705
HL TAU | 1.30 | 3.2 | 1.413 | 0.108 | 13 | 0.432 | 0.410 | 0.200 | 0.097 | 0.161 | 0.613
HL TAU | 1.30 | 3.2 | 1.413 | 0.108 | 32 | 0.152 | 0.670 | 0.270 | 0.041 | 0.097 | 0.636
HL TAU | 1.30 | 3.2 | 1.413 | 0.108 | 67 | 0.064 | 0.840 | 0.550 | 0.047 | 0.090 | 0.613
AS209 | 0.83 | 0.5 | 1.096 | 0.110 | 99 | 0.313 | 0.030 | 0.650 | 0.054 | 0.109 | 0.363
HD142666 | 1.58 | 9.1 | 12.589 | 0.316 | 16 | 0.219 | 0.730 | 0.300 | 0.051 | 0.123 | 0.230
HD143006 | 1.78 | 3.8 | 3.981 | 0.273 | 22 | 0.986 | 0.040 | 20.000 | 0.440 | 0.440 | 0.440
HD143006 | 1.78 | 3.8 | 3.981 | 0.273 | 51 | 0.251 | 0.530 | 0.330 | 0.027 | 0.068 | 0.257
HD163296 | 2.04 | 17.0 | 5.012 | 0.210 | 48 | 0.421 | 0.033 | 2.180 | 0.345 | 0.372 | 0.501
HD163296 | 2.04 | 17.0 | 5.012 | 0.210 | 86 | 0.188 | 0.150 | 1.000 | 0.045 | 0.090 | 0.379
SR4 | 0.68 | 1.2 | 0.80 | 0.136 | 11 | 0.573 | 0.230 | 2.160 | 0.235 | 0.235 | 0.235
HD169142 | 1.65 | 10.0 | 5.012 | 0.330 | 14 | 1.857 | 0.006 | 3.500 | 0.205 | 0.205 | 0.205
HD169142 | 1.65 | 10.0 | 5.012 | 0.330 | 40 | 0.750 | 0.012 | 0.700 | 0.061 | 0.097 | 0.142
PDS70 | 0.85 | 2.8 | 5.370 | 0.170 | 22 | 0.055 | 0.000 | 10.000 | 0.353 | 0.353 | 0.353

## Acknowledgement

Yinhao Wu is thanked for helping the authors to compare the migration rate of
planets in the 1D code with the 2D code FARGO3D (see Appendix). The authors
acknowledge funding from the UK Science and Technology Facilities Council
under the grant No. ST/S000453/1. GR acknowledges support from an STFC Ernest
Rutherford Fellowship (grant number ST/T003855/1). VE acknowledges support by
the Ministry of Science and Higher Education of the Russian Federation under
the grant 075-15-2020-780 (N13.1902.21.0039).

## Data availability statement

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * ALMA Partnership et al. (2015) ALMA Partnership et al., 2015, ApJ, 808, L3
  * Alexander et al. (2006) Alexander R. D., Clarke C. J., Pringle J. E., 2006, MNRAS, 369, 229
  * Alexander et al. (2014) Alexander R., Pascucci I., Andrews S., Armitage P., Cieza L., 2014, Protostars and Planets VI, pp 475–496
  * Andrews (2020) Andrews S. M., 2020, ARA&A, 58, 483
  * Andrews et al. (2018) Andrews S. M., et al., 2018, ApJ, 869, L41
  * Armitage & Bonnell (2002) Armitage P. J., Bonnell I. A., 2002, MNRAS, 330, L11
  * Barge et al. (2017) Barge P., Ricci L., Carilli C. L., Previn-Ratnasingam R., 2017, A&A, 605, A122
  * Bate et al. (2003) Bate M. R., Bonnell I. A., Bromm V., 2003, MNRAS, 339, 577
  * Benisty et al. (2021) Benisty M., et al., 2021, ApJ, 916, L2
  * Benítez-Llambay & Masset (2016) Benítez-Llambay P., Masset F. S., 2016, ApJS, 223, 11
  * Benítez-Llambay et al. (2015) Benítez-Llambay P., Masset F., Koenigsberger G., Szulágyi J., 2015, Nature, 520, 63
  * Bennett et al. (2021) Bennett D. P., Ranc C., Fernandes R. B., 2021, AJ, 162, 243
  * Birnstiel et al. (2009) Birnstiel T., Dullemond C. P., Brauer F., 2009, A&A, 503, L5
  * Birnstiel et al. (2012) Birnstiel T., Klahr H., Ercolano B., 2012, A&A, 539, A148
  * Booth et al. (2019) Booth A. S., Walsh C., Ilee J. D., Notsu S., Qi C., Nomura H., Akiyama E., 2019, ApJ, 882, L31
  * Casassus & Pérez (2019) Casassus S., Pérez S., 2019, ApJ, 883, L41
  * Coleman & Nelson (2014) Coleman G. A. L., Nelson R. P., 2014, MNRAS, 445, 479
  * Coleman & Nelson (2016) Coleman G. A. L., Nelson R. P., 2016, MNRAS, 460, 2779
  * Crida & Bitsch (2017) Crida A., Bitsch B., 2017, Icarus, 285, 145
  * 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’Alessio et al. (2001) D’Alessio P., Calvet N., Hartmann L., 2001, ApJ, 553, 321
  * Dipierro & Laibe (2017) Dipierro G., Laibe G., 2017, MNRAS, 469, 1932
  * Dipierro et al. (2016) Dipierro G., Laibe G., Price D. J., Lodato G., 2016, MNRAS,
  * Dullemond et al. (2018) Dullemond C. P., et al., 2018, ApJ, 869, L46
  * Dullemond et al. (2020) Dullemond C. P., Isella A., Andrews S. M., Skobleva I., Dzyurkevich N., 2020, A&A, 633, A137
  * Emsenhuber et al. (2021a) Emsenhuber A., Mordasini C., Burn R., Alibert Y., Benz W., Asphaug E., 2021a, A&A, 656, A69
  * Emsenhuber et al. (2021b) Emsenhuber A., Mordasini C., Burn R., Alibert Y., Benz W., Asphaug E., 2021b, A&A, 656, A70
  * Flaherty et al. (2017) Flaherty K. M., et al., 2017, ApJ, 843, 150
  * Fung & Chiang (2016) Fung J., Chiang E., 2016, ApJ, 832, 105
  * Gonzalez et al. (2015) Gonzalez J. F., Laibe G., Maddison S. T., Pinte C., Ménard F., 2015, MNRAS, 454, L36
  * Guilera et al. (2021) Guilera O. M., Miller Bertolami M. M., Masset F., Cuadra J., Venturini J., Ronco M. P., 2021, MNRAS, 507, 3638
  * Haffert et al. (2019) Haffert S. Y., Bohn A. J., de Boer J., Snellen I. A. G., Brinchmann J., Girard J. H., Keller C. U., Bacon R., 2019, Nature Astronomy, 3, 749
  * Hall et al. (2019) Hall C., Dong R., Rice K., Harries T. J., Najita J., Alexander R., Brittain S., 2019, ApJ, 871, 228
  * Huang et al. (2018) Huang J., et al., 2018, ApJ, 869, L42
  * Ida & Lin (2004a) Ida S., Lin D. N. C., 2004a, ApJ, 604, 388
  * Ida & Lin (2004b) Ida S., Lin D. N. C., 2004b, ApJ, 616, 567
  * Ikoma et al. (2000) Ikoma M., Nakazawa K., Emori H., 2000, ApJ, 537, 1013
  * Jones et al. (2012) Jones M. G., Pringle J. E., Alexander R. D., 2012, MNRAS, 419, 925
  * Jung et al. (2019) Jung Y. K., et al., 2019, AJ, 157, 72
  * Keppler et al. (2019) Keppler M., et al., 2019, A&A, 625, A118
  * Kimmig et al. (2020) Kimmig C. N., Dullemond C. P., Kley W., 2020, A&A, 633, A4
  * Lambrechts & Lega (2017) Lambrechts M., Lega E., 2017, A&A, 606, A146
  * Lambrechts et al. (2019) Lambrechts M., Lega E., Nelson R. P., Crida A., Morbidelli A., 2019, A&A, 630, A82
  * Lin & Papaloizou (1986) Lin D. N. C., Papaloizou J., 1986, ApJ, 309, 846
  * Liu et al. (2019) Liu Y., et al., 2019, A&A, 622, A75
  * Lodato & Clarke (2004) Lodato G., Clarke C. J., 2004, MNRAS, 353, 841
  * Lodato et al. (2009) Lodato G., Nayakshin S., King A. R., Pringle J. E., 2009, MNRAS, 398, 1392
  * Lodato et al. (2019) Lodato G., et al., 2019, MNRAS, 486, 453
  * Long et al. (2018) Long F., et al., 2018, ApJ, 869, 17
  * Manara et al. (2017) Manara C. F., et al., 2017, A&A, 604, A127
  * Masset (2017) Masset F. S., 2017, MNRAS, 472, 4204
  * Mayor et al. (2011) Mayor M., et al., 2011, ArXiv e-prints (astro-ph 1109.2497),
  * Meru et al. (2019) Meru F., Rosotti G. P., Booth R. A., Nazari P., Clarke C. J., 2019, MNRAS, 482, 3678
  * Mizuno (1980) Mizuno H., 1980, Progress of Theoretical Physics, 64, 544
  * 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
  * Mordasini et al. (2012) Mordasini C., Alibert Y., Benz W., Klahr H., Henning T., 2012, A&A, 541, A97
  * Muro-Arena et al. (2018) Muro-Arena G. A., et al., 2018, A&A, 614, A24
  * Nayakshin & Fletcher (2015) Nayakshin S., Fletcher M., 2015, MNRAS, 452, 1654
  * Nayakshin et al. (2019) Nayakshin S., Dipierro G., Szulágyi J., 2019, MNRAS, 488, L12
  * Ndugu et al. (2019) Ndugu N., Bitsch B., Jurua E., 2019, MNRAS, 488, 3625
  * Nelson & Papaloizou (2004) Nelson R. P., Papaloizou J. C. B., 2004, MNRAS, 350, 849
  * Ormel et al. (2015) Ormel C. W., Shi J.-M., Kuiper R., 2015, MNRAS, 447, 3512
  * Paardekooper et al. (2010) Paardekooper S.-J., Baruteau C., Crida A., Kley W., 2010, MNRAS, 401, 1950
  * Paardekooper et al. (2011) Paardekooper S. J., Baruteau C., Kley W., 2011, MNRAS, 410, 293
  * Perez et al. (2015) Perez S., Dunhill A., Casassus S., Roman P., Szulágyi J., Flores C., Marino S., Montesinos M., 2015, ApJ, 811, L5
  * Pérez et al. (2018) Pérez S., Casassus S., Benítez-Llambay P., 2018, MNRAS, 480, L12
  * Pinilla et al. (2012) Pinilla P., Birnstiel T., Ricci L., Dullemond C. P., Uribe A. L., Testi L., Natta A., 2012, A&A, 538, A114
  * Pinte et al. (2018) Pinte C., et al., 2018, ApJ, 860, L13
  * Pinte et al. (2020) Pinte C., et al., 2020, ApJ, 890, L9
  * Piso & Youdin (2014) Piso A.-M. A., Youdin A. N., 2014, ApJ, 786, 21
  * Piso et al. (2015) Piso A.-M. A., Youdin A. N., Murray-Clay R. A., 2015, ApJ, 800, 82
  * Pollack et al. (1996) Pollack J. B., Hubickyj O., Bodenheimer P., Lissauer J. J., Podolak M., Greenzweig Y., 1996, Icarus, 124, 62
  * Powell et al. (2017) Powell D., Murray-Clay R., Schlichting H. E., 2017, ApJ, 840, 93
  * Powell et al. (2019) Powell D., Murray-Clay R., Pérez L. M., Schlichting H. E., Rosenthal M., 2019, ApJ, 878, 116
  * Qi et al. (2011) Qi C., D’Alessio P., Öberg K. I., Wilner D. J., Hughes A. M., Andrews S. M., Ayala S., 2011, ApJ, 740, 84
  * Rab et al. (2020) Rab C., Kamp I., Dominik C., Ginski C., Muro-Arena G. A., Thi W. F., Waters L. B. F. M., Woitke P., 2020, A&A, 642, A165
  * Rice et al. (2006) Rice W. K. M., Armitage P. J., Wood K., Lodato G., 2006, MNRAS, 373, 1619
  * Rosotti et al. (2016) Rosotti G. P., Juhasz A., Booth R. A., Clarke C. J., 2016, MNRAS, 459, 2790
  * Rosotti et al. (2019) Rosotti G. P., Tazzari M., Booth R. A., Testi L., Lodato G., Clarke C., 2019, MNRAS, 486, 4829
  * Ruge et al. (2016) Ruge J. P., Flock M., Wolf S., Dzyurkevich N., Fromang S., Henning T., Klahr H., Meheut H., 2016, A&A, 590, A17
  * Safronov (1969) Safronov V. S., 1969, Evoliutsiia doplanetnogo oblaka.
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Stevenson (1982) Stevenson D. J., 1982, P&SS, 30, 755
  * Suzuki et al. (2016) Suzuki D., et al., 2016, ApJ, 833, 145
  * Suzuki et al. (2018) Suzuki D., et al., 2018, ApJ, 869, L34
  * Syer & Clarke (1995) Syer D., Clarke C. J., 1995, MNRAS, 277, 758
  * Szulágyi et al. (2014) Szulágyi J., Morbidelli A., Crida A., Masset F., 2014, ApJ, 782, 65
  * Tabone et al. (2021) Tabone B., Rosotti G. P., Cridland A. J., Armitage P. J., Lodato G., 2021, MNRAS,
  * Teague et al. (2019) Teague R., Bae J., Bergin E. A., 2019, Nature, 574, 378
  * Williams & Cieza (2011) Williams J. P., Cieza L. A., 2011, ARA&A, 49, 67
  * Winn & Fabrycky (2015) Winn J. N., Fabrycky D. C., 2015, ARA&A, 53, 409
  * Zhang et al. (2015) Zhang K., Blake G. A., Bergin E. A., 2015, ApJ, 806, L7
  * Zhang et al. (2018) Zhang S., et al., 2018, ApJ, 869, L47

##  Appendix A Gaseous disc evolution and planet migration

Here we provide some finer technical detail of planet-disc angular momentum
exchange implementation in our code. The total torque $\Lambda_{t}$ from the
planet onto the disc is the integral over the entire disc,

$\Lambda_{t}=\int_{R_{\rm in}}^{R_{\rm out}}\lambda_{t}(R)\Sigma(R)2\pi
RdR\;.$ (22)

Strictly speaking, the widely used type II torque implementation given by eq.
5 is well justified for massive planets that open deep gaps in the discs (Type
II migration). For less massive planets, type I migration occurs, for which
the corotation torques are important (e.g., Paardekooper et al., 2010). Eq. 5
then usually under-estimates the migration rate of the planet. A frequent
choice for $\Lambda_{t}$ in this regime are the expressions from Paardekooper
et al. (2011). In planet formation population synthesis one then also neglects
the reverse (planetary) torque on the disc, setting $\lambda_{t}=0$, arguing
that in type I the planet is not very massive and hence can be thought as a
test particle. Such an approach is used, for example, by Coleman & Nelson
(2014, 2016); Emsenhuber et al. (2021a, b).

Unfortunately, the abrupt transition from using eq. 5 for $\lambda_{t}$ in
type II regime to setting it to zero in type I is not appropriate for us here
since we aim to model both planet migration and tidal torque effects on the
gas and dust components of the disc. The latter is more strongly influenced by
the planet for Stokes numbers not too small (e.g., Dipierro & Laibe, 2017).

Therefore, instead of neglecting the planet torques on the disc in type I
regime, we blend the two approaches used in the literature previously. In the
type II regime, that is, when $C_{\rm p}\leq 1$, we use eq. 5. In Type I
regime, when $C_{\rm p}>1$, instead of setting $\lambda_{t}$ to zero we write

$\displaystyle\lambda_{t}=g\lambda_{\rm out}$ $\displaystyle\text{ for }R>a$
(23) $\displaystyle\lambda_{t}=g^{-1}\lambda_{\rm in}$ $\displaystyle\text{
for }R<a\;,$

where $\lambda_{\rm out}=\lambda_{t}(R>a)$, $\lambda_{\rm
in}=\lambda_{t}(R<a)$ as given in eq. 5, and $g$ is a factor of order unity
introduced to re-normalise the tidal torque on both sides of the planet to
obtain the expected planet migration rate. We find $g$ by computing the
integral tidal torque, $\Lambda_{t}$ from the planet on the disc, and
demanding that it equals to that given by the Type I torque $\Lambda_{1}$
Paardekooper et al. (2011), with the minus sign,

$g\Lambda_{\rm out}+g^{-1}\Lambda_{\rm in}=-\Lambda_{1}\;,$ (24)

where $\Lambda_{\rm in}$ and $\Lambda_{\rm out}$ are given by the integral in
eq. 22, except only inside and outside $a$, respectively. In practice $g$ is
usually between $1$ and $2$. Note that one may think that a symmetric re-
normalisation of eq. 5 is a more logical idea, however we found that in
general such an approach requires $g$ that can be either much larger or much
smaller than unity. Furthermore, $g$ may even have negative values if type I
migration direction is opposite to that predicted to type II.

Comparing our 1D planet migration treatment with 2D simulations with the
public code FARGO3D (Benítez-Llambay & Masset, 2016) for a particular disc
setup, and a range of disc masses, we found a reasonable agreement. The 1D
planet migration time was a factor $f_{1to2}=$1-3 times shorter than the 2D
result. The variation in $f_{1to2}$ is mainly due to fact that the 1D
treatment cannot capture the type III runaway migration that may occur on the
boundary between type I and type II regimes (e.g., Crida & Morbidelli, 2007;
Crida & Bitsch, 2017).

##  Appendix B A colder disc in HD163296

As remarked in §3, we chose a rather hot disc to study in HD163296 by setting
$T_{0}=65$ K at $R_{0}=100$ AU following Rab et al. (2020). This temperature
is however found at the CO molecule emitting surface, and is larger than
$T\approx 25$ K temperature derived for the midplane of the disc in this
source by Dullemond et al. (2020). Here we show that choosing this lower
temperature results in even faster planet growth and hence makes our main
conclusions even stronger.

Fig. 10 shows calculation of the disc-planet system evolution repeated in
exactly the same way as shown in Fig. 2, but now for a colder disc, with
$T_{0}=25$ K. Since gas accretion and gas capture radii are larger in colder
discs (cf. eq. 10), gas accretion rate is higher in a colder disc too.
Therefore, planet runaway to their final masses more rapidly, and hence the
time that the gap properties are commensurate with a moderately massive planet
becomes yet shorter, that is, less than $0.1$ Myr for the $\alpha_{\rm
v}=10^{-3}$ case presented in Fig. 10. For the same reason the planet opens a
deep gap in the gas and switches into the type II migration regime earlier and
remains at large separations for longer. While retention of the planet at
large separation for a few Myr is a welcome result, the gas and dust disc
morphology are very different from those observed in HD163296 when the planet
is in type II migration regime. Therefore we conclude that our main results
are robust with respect to changes in the assumed disc temperature.

Figure 10: Same as Fig. 2 but for a colder disc with $T_{0}=25$ K. Note that
the planet accretion is faster in this colder disc, therefore it reaches its
final mass faster. For the same reason it opens a deep gap and switches to
type II migration more rapidly than in the hotter disc studied in Fig. 2.