11institutetext:  University Observatory, Faculty of Physics, Ludwig-
Maximilians-Universität München, Scheinerstr. 1, 81679 München, Germany
11email<EMAIL_ADDRESS>22institutetext: Exzellenzcluster ORIGINS, Boltzmannstr. 2, D-85748 Garching, Germany 33institutetext: Department of Physics and Astronomy, University of Leicester, University Road, Leicester LE1 7RH, UK 44institutetext: Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, the Netherlands 55institutetext: Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

# A Large Population Study of Protoplanetary Disks

Explaining the millimeter Size-Luminosity relation with or without sub-
structure

Apostolos Zormpas 11 Tilman Birnstiel 1122 Giovanni P. Rosotti 3344 Sean M.
Andrews 55

(Received ; accepted XXXX)

Recent sub-arcsecond resolution surveys of the dust continuum emission from
nearby protoplanetary disks showed a strong correlation between the sizes and
luminosities of the disks. We aim to explain the origin of the
(sub-)millimeter size-luminosity relation (SLR) between the $\mathrm{68\%}$
effective radius ($\mathrm{r_{eff}}$) of disks with their continuum luminosity
($\mathrm{L_{mm}}$), with models of gas and dust evolution in a simple viscous
accretion disk and radiative transfer calculations.

We use a large grid of models ($\mathrm{10^{5}}$ simulations) with and without
planetary gaps, varying the initial conditions of the key parameters. We
calculate the disk continuum emission and the effective radius for all models
as a function of time. By selecting those simulations that continuously follow
the SLR, we can derive constraints on the input parameters of the models.

We confirm previous results that models of smooth disks in the radial drift
regime are compatible with the observed SLR ($\mathrm{L_{mm}\propto
r_{eff}^{2}}$) but only smooth disks cannot be the reality. We show that the
SLR is more widely populated if planets are present. However they tend to
follow a different relation than smooth disks, potentially implying that a
mixture of smooth and sub-structured disks are present in the observed sample.
We derive a SLR ($\mathrm{L_{mm}\propto r_{eff}^{5/4}}$) for disks with strong
sub-structure. To be compatible with the SLR, models need to have an initially
high disk mass ($\mathrm{\geq 2.5\cdot 10^{-2}M_{\star}}$) and low turbulence-
parameter $\mathrm{\alpha}$ values ($\mathrm{\leq 10^{-3}}$). Furthermore, we
find that the grain composition and porosity drastically affects the evolution
of disks on the size-luminosity diagram where relatively compact grains that
include amorphous carbon are favoured. Moreover, a uniformly optically thick
disk with high albedo ($\rm{0.9}$) that follows the SLR cannot be formed from
an evolutionary procedure.

###### Key Words.:

Accretion disks – Planets and satellites: formation – Protoplanetary disks –
Circumstellar matter

##  1 Introduction

Planet formation is a process far from having a complete, robust and widely
accepted theory. Multiple theories aim to explore the way planets form (e.g.
Benz et al., 2014; Kratter & Lodato, 2016; Johansen & Lambrechts, 2017). To
understand planet formation, high resolution observations and large surveys of
the birth places of planets, the protoplanetary disks (PPDs) are essential. In
recent years, the Atacama Large Millimeter/Submillimeter Array (ALMA) has not
only provided a sizable number of highly resolved observations of PPDs but due
to its high sensitivity it also enabled several large, mid-resolution (of the
order of $100\text{\,}\mathrm{m}\mathrm{a}\mathrm{s}$) surveys of different
star-forming regions (for a review see Andrews, 2020, and references therein)
providing crucial information on population properties such as distributions
of disk sizes, fluxes, or spectral indices for disks across stars of different
masses and ages and across star-forming regions in different environments.

Highly resolved observations and large mid-resolution surveys are
complementary and are both essential for understanding the connections between
key properties of PPDs. One of the important diagnostics is the continuum
luminosity ($\mathrm{L_{mm}}$) at (sub-)millimeter (hereafter ”mm”)
wavelengths, since it can be a tracer of the disk mass that is produced by the
solid grains (Beckwith et al., 1990), i.e. the amount of material available to
form planets. Assuming a constant dust-to-gas ratio (usually $\rm{0.01}$ based
on observations of the interstellar medium, see Bohlin et al., 1978), the dust
mass can be converted to the total disk mass (dust and gas). Surveys that have
measured the disk dust mass ($\mathrm{M_{dust}}$), have been used to correlate
this property with the mass of the host star ($\mathrm{M_{\star}}$), and a
linear relation has been found between them (Andrews et al., 2013) that
appears to steepen with time (Pascucci et al., 2016; Ansdell et al., 2017;
Barenfeld et al., 2016). Furthermore, it has been observed a steeper than
linear relationship between $\mathrm{L_{mm}}$ and $\mathrm{M_{\star}}$
(Andrews et al., 2013; Pascucci et al., 2016; Ansdell et al., 2016). Recently,
theoretical studies have started to explain these correlations (e.g. Pascucci
et al., 2016; Stammler et al., 2019; Pinilla et al., 2020) using numerical
dust evolution models that include grain growth, radial drift, fragmentation
of dust particles, and particle traps (e.g. Birnstiel et al., 2010, 2012;
Krijt et al., 2016).

These observations of dust thermal continuum emission are crucial for the
characterization of disk evolution. However, the methods described above that
relate dust emission to physical quantities like total disk masses carry
uncertainty arising from the multiple assumptions such as the dust opacity and
the disk temperature that may vary for every disk (e.g. Hendler et al., 2017;
Ballering & Eisner, 2019). The grain opacity depends on the unknown particle
size distribution, composition, and particle structure (e.g. Birnstiel et al.,
2018, and references within), although considerable efforts have been
undertaken from modeling (e.g. Wada et al., 2008, 2009; Okuzumi et al., 2009;
Seizinger & Kley, 2013; Seizinger et al., 2013) and experimental studies (e.g.
Blum & Wurm, 2008; Güttler et al., 2010; Gundlach & Blum, 2015) of aggregation
and of the computation of optical properties of aggregates (e.g. Kataoka et
al., 2014; Min et al., 2016; Tazaki et al., 2016).

Another important property for the characterization of a disk population is
the disk size. Viscous theory (Lynden-Bell & Pringle, 1974) predicts that a
fraction of the disk mass keeps moving outwards, the so-called viscous
spreading, suggesting the disk size should increase with time. In principle a
measurement of disk size as function of time could measure this evolution to
test viscous theory and measure its efficiency. The most readily available
tracer of the disk size is the continuum as gas tracers suffer from uncertain
abundances (due to freeze-out, dissociation, and other) and sensitivity
constraints. Since the disk does not have a clear outer edge, we need to
introduce an effective radius and express the size as a function of the total
continuum emission. This size metric is called emission size or effective
radius ($\mathrm{r_{eff}}$) (Tripathi et al., 2017).

However, the dust component behaves differently than the gas, mainly due to an
effect termed radial drift (e.g. Whipple, 1972; Weidenschilling, 1977;
Takeuchi & Lin, 2002). The dust particles interact with the sub-Keplerian gas
disk via aerodynamic drag forces, leading them to migrate towards the star. As
an observational implication, the dust emission is less extended than the gas
emission (e.g. Andrews et al., 2012; Isella et al., 2012; Andrews et al.,
2016; Cleeves et al., 2016) as predicted by Birnstiel & Andrews (2014) (but
see Trapman et al., 2020). Radial drift is also heavily dependent on the grain
size, therefore grain growth (e.g. Birnstiel et al., 2012) has to be included
in the numerical studies that aim to use the dust disk radius. Rosotti et al.
(2019b) studied theoretically how the evolution of the disk dust radius
changes with time in a viscously evolving disk and addressed whether the
evolution of the dust disk radius is set by viscous spreading or by the dust
processes such as grain growth and radial drift. They found that viscous
spreading influences the dust and leads to the dust disk expanding with time.

Many surveys have been performed to explore the relation between the two
diagnostics (e.g. Andrews et al., 2010; Piétu et al., 2014; Hendler et al.,
2020). Recently, a sub-arcsecond resolution survey of 50 nearby protoplanetary
disks, conducted with the Submillimeter Array (SMA) by Tripathi et al. (2017),
showed a strong size-luminosity relation (hereafter ”SLR”) between the
observed population. The follow-up program in Andrews et al. (2018b), a
combined analysis of the Tripathi et al. (2017) data and ALMA data from the
Lupus disk sample (105 disks in total) confirmed the scaling relations between
the $\mathrm{r_{eff}}$ and the $\mathrm{L_{mm}}$, $\mathrm{M_{\star}}$.
However, not all studied star-forming regions show the same correlation but
appear to vary with the age of the region (Hendler et al., 2020).

In recent years, due to the unprecedented sensitivity and resolution, ALMA
provided a plethora of groundbreaking images of protoplanetary disks. Most of
these disks do not show a smooth and monotonically decreasing surface density
profile, but instead are composed of single or multiple symmetric annular sub-
structures instead e.g., HL Tauri (ALMA Partnership et al., 2015), TW Hya
(Andrews et al., 2016; Tsukagoshi et al., 2016), HD 163296 (Isella et al.,
2016), HD 169142 (Fedele et al., 2017), AS 209 (Fedele et al., 2018), HD
142527 (Casassus et al., 2013) and many more in the recent DSHARP survey
(Andrews et al., 2018a). Moreover, non-axisymmetric features like spiral arms
(e.g. Pérez et al., 2016; Huang et al., 2018b) and lopsided rings (e.g. van
der Marel et al., 2013) have been observed.

Many ideas tried to explore the origin of these ring-like sub-structures, but
one of the most favorable explanations is the formation of gaps in the gas
surface density, due to the presence of planets. A massive planet ($\rm{\geq
0.1M_{jup}}$) (Zhang et al., 2018) is able to open a gap in the surrounding
gaseous disk, thereby generating a pressure maximum. Later on, the dust
particles are migrating towards the local pressure bump due to radial drift
(Weidenschilling, 1977; Nakagawa et al., 1986) and consequently leading to the
annular shape (e.g Rice et al., 2006; Pinilla et al., 2012a). These narrow
rings may be optically thick or moderately optically thick, but in between
these features, the material is approximated as optically thin (Dullemond et
al., 2018). These rings however, can contain large amounts of dust that can
increase the total luminosity of a disk and its position with respect to the
SLR.

The SLR might contain crucial information about disk evolution and planet
formation theory. Our goal is to explore the physical origins of the SLR from
Tripathi et al. (2017) and Andrews et al. (2018b) by performing a large
population study of models with gas and dust evolution. We aim to characterize
the key properties of disks that reproduce the observational results. We
explore the differences in the SLR of disks that have a smooth surface density
profile and disks that contain weak and strong sub-structures. In section 2,
we discuss the methods we used to carry out our computational models. The
results of this analysis are presented in section 3, where we explain the
global effect of every parameter in the population of disks and we present the
general properties that disks should have to follow the SLR. In section 4 we
discuss the theoretical and observational implications of our results. We draw
our conclusions in section 5.

##  2 Methods

We carry out 1-D gas and dust evolution simulations using a slightly modified
version of the two-population model (two-pop-py) by Birnstiel et al. (2012,
2015), while we also mimic the presence of planets. As a post-processing step,
we calculate the intensity profile and the disk continuum emission. With the
purpose of running a population study, we use a large grid of parameters (see
Table 1), so that we can explore the differences that occur due to the
different initial conditions. In the next sections, we explain the procedure
in more detail.

###  2.1 Disk evolution

The gas is following the viscous evolution equation. For the disk evolution,
we use the turbulent effective viscosity as in Shakura & Sunyaev (1973),

$\mathrm{\nu=\alpha_{gas}\frac{c_{s}^{2}}{\Omega_{K}}}\,,$ (1)

and the dust diffusion coefficient as,

$\mathrm{D=\alpha_{dust}\frac{c_{s}^{2}}{\Omega_{K}}}\,,$ (2)

with $\mathrm{\alpha_{gas}}$ being the turbulence parameter, $\mathrm{c_{s}}$
the sound speed and $\mathrm{\Omega_{K}}$ the Keplerian frequency. The above
equation lacks the term $\rm{\frac{1}{1+St^{2}}}$, where $\rm{St}$ is the
Stokes number, but we can ignore it since the Stokes number is always
$\rm{<1}$ in our simulations (Youdin & Lithwick (2007)). We differentiate the
$\alpha$ parameter, in two different values. One for the gas
$\alpha_{\mathrm{gas}}$, and one for the dust $\alpha_{\mathrm{dust}}$, since
we later mimic planetary gaps by locally varying the viscosity (see subsection
2.3). In the smooth case $\alpha_{\mathrm{dust}}=\alpha_{\mathrm{gas}}$.

The dust is described by the two populations model of Birnstiel et al. (2012)
which evolves the dust surface density under the assumption that the small
dust is tightly coupled to the gas while the large particles can decouple from
the gas and drift inwards. The initial dust growth phase is using the current
dust-to-gas ratio, instead of the initial value as in Birnstiel et al. (2012).

We set the initial gas surface density according to the self-similar solution
of Lynden-Bell & Pringle (1974).

$\mathrm{\Sigma_{g}(r)}=\mathrm{\Sigma_{0}\left(\frac{r}{r_{c}}\right)^{-\gamma}Exp\left[-\left(\frac{r}{r_{c}}\right)^{2-\gamma}\right]}\,,$
(3)

where $\mathrm{\Sigma_{0}=(2-\gamma)M_{d}/2\pi r_{\mathrm{c}}^{2}}$ is the
normalization parameter which is set in every simulation by the disk mass
$\mathrm{M_{d}}$. The other parameters are $\mathrm{\gamma}=1$, which is the
viscosity exponent and it is not varied throughout our models and
$r_{\mathrm{c}}$ which is the characteristic radius of the disk (see Table 1).
When r $\ll$ $r_{\mathrm{c}}$, then $\mathrm{\Sigma_{g}}$ is a power law and
when $\mathrm{r\geq r_{\mathrm{c}}}$, $\mathrm{\Sigma_{g}}$ is dominated by
the exponential factor.

The initial dust distribution follows the gas distribution with a constant
dust-to-gas ratio of $\mathrm{\Sigma_{d}}$/$\mathrm{\Sigma_{g}}$=0.01. The
initial grain size (= monomer grain size) is
$a_{\mathrm{min}}=$0.1\text{\,}\mu\mathrm{m}$$. This monomer size stays
constant in time and space while the representative size for the large grains
increases with time as particles grow. The particle bulk density is
$\rm{\rho_{s}}=$1.7\text{\,}\mathrm{g}\mathrm{/}\mathrm{c}\mathrm{m}^{3}$$ for
the standard opacity model from Ricci et al. (2010) and
$\rm{\rho_{s}}=$1.675\text{\,}\mathrm{g}\mathrm{/}\mathrm{c}\mathrm{m}^{3}$$
for the DSHARP (Birnstiel et al., 2018) opacity but decreases for different
values of porosity (see subsection 2.4). We evolve the disks up to
$10\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ to study the long-term evolution
but in the following analysis we only show results from
$300\text{\,}\mathrm{k}\mathrm{y}\mathrm{r}$ to
$3\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ (see subsection 2.5).

Our $\mathrm{1-D}$ radial grid ranges from $\rm{0.05}$ to
$2000\text{\,}\mathrm{a}\mathrm{u}$ and the grid cells are spaced
logarithmically. We are using adaptive temperature that depends on the
luminosity of every star. Since the stellar mass changes in our grid, the
stellar luminosity is changing too. We follow the temperature profile:

$\mathrm{T=\left(\phi\frac{L_{\star}}{4\pi\sigma_{SB}r^{2}}+($10\text{\,}\mathrm{K}$)^{4}\right)^{1/4}}\,,$
(4)

as is Kenyon et al. (1996). In this equation, $\mathrm{L_{\star}}$ is the
stellar luminosity, $\mathrm{\phi}=0.05$ is the flaring angle, while
$\mathrm{\sigma_{SB}}$ is the Stefan-Boltzmann constant and $\mathrm{r}$ the
radius. The term $\mathrm{10^{4}}$ is a floor value so that we do not allow
the disk temperature to drop below $10\text{\,}\mathrm{K}$ at the outer parts
of the disk. We are using the evolutionary tracks of Siess et al. (2000) to
get the luminosity of an $1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ old star
of the given mass. The stellar luminosity and effective temperature is not
evolved in our simulations. However, the luminosity of a
$1\text{\,}\mathrm{M}_{\odot}$ star would decrease from
$\rm{\sim$2.4\text{\,}\mathrm{L}_{\odot}$}$ at
$1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ to
$\rm{\sim$1\text{\,}\mathrm{L}_{\odot}$}$ at
$3\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$. In subsection 4.8 we explore how a
change in stellar luminosity affects our results.

###  2.2 Population study

We use an extended parameter grid, by varying the initial values of the
turbulence parameter ($\mathrm{\alpha_{gas}}$), disk mass ($\mathrm{M_{d}}$),
stellar mass ($\mathrm{M_{\star}}$), characteristic radius ($r_{\mathrm{c}}$)
and fragmentation velocity ($\mathrm{v_{frag}}$). For every parameter, we pick
$\mathrm{10}$ values specified in Table 2, taking all the possible
combinations between them leading to a total of $\mathrm{100.000}$
simulations.

Table 1: Grid parameters of the model Parameter | Description | Value-Range
---|---|---
$\mathrm{\Sigma_{d}/\Sigma_{g}}$ | initial dust-to-gas ratio | $\rm{0.01}$
$\mathrm{\rho_{s}}$ [$\rm{g/cm^{3}}$] | particle bulk density | $\rm{1.7,1.675}$
$\mathrm{\gamma}$ | viscosity exponent | $\rm{1}$
$r$ [$\mathrm{au}$] | logarithmic grid extent | $\rm{0.05-2000}$
$\mathrm{n_{r}}$ [cells] | grid resolution | $\rm{400}$
$\mathrm{t}$ [years] | duration of each simulation | $\rm{10^{7}}$
Table 2: Variables of the model Parameter | Description | Values
---|---|---
$\mathrm{\alpha}$ | viscosity parameter | $\mathrm{\\{1,2.5,5,7.5\\}\cdot 10^{-4}}$
|  | $\mathrm{\\{1,2.5,5,7.5\\}\cdot 10^{-3}}$
|  | $\mathrm{\\{1,2.5\\}\cdot 10^{-2}}$
$\mathrm{M_{d}}$ [M⋆] | initial disk mass | $\mathrm{\\{1,2.5,5,7.5\\}\cdot 10^{-3}}$
|  | $\mathrm{\\{1,2.5,5,7.5\\}\cdot 10^{-2}}$
|  | $\mathrm{\\{1,2.5\\}\cdot 10^{-1}}$
$\mathrm{M_{\star}}$ [M⊙] | stellar mass | $\mathrm{0.2,0.4,0.6,0.8,1.0}$
|  | $\mathrm{1.2,1.4,1.6,1.8,2.0}$
$r_{\mathrm{c}}$ [$au$] | characteristic radius | $\mathrm{10,30,50,80,100}$
|  | $\mathrm{130,150,180,200,230}$
$\mathrm{v_{f}}$ [$cm/s$] | fragmentation velocity | $\mathrm{200,400,600,800,1000}$
|  | $\mathrm{1200,1400,1600,1800,2000}$
q | planet/star mass ratio |  $3\cdot 10^{-4}$, $10^{-3}$, $3\cdot 10^{-3}$
rp [rc] | planet position |  $1/3$, $2/3$

###  2.3 Planets

A large planet embedded in a disk produces a co-orbital gap in the gas
density. To mimick gap opening by planets in our simulations, we altered the
$\mathrm{\alpha_{gas}}$ turbulence parameter. Since in steady state
$\mathrm{\alpha_{gas}}$ $\mathrm{\cdot}$ $\mathrm{\Sigma_{g}}$ is constant,
the $\mathrm{\alpha}$ parameter and the surface density $\mathrm{\Sigma_{g}}$
are inversely proportional quantities, so a bump in the
$\mathrm{\alpha_{gas}}$ profile leads to a gap in the surface density profile.
The reason for the change in $\mathrm{\alpha_{gas}}$ and not in
$\mathrm{\Sigma_{g}}$ is that the surface density evolves according to the
Equation 3. By inserting the bump in the $\mathrm{\alpha_{gas}}$ the
$\mathrm{\Sigma_{g}}$ still evolves viscously and at the same time produces a
planetary gap shape.

Following the prescription from Kanagawa et al. (2016), we mimic the effect of
planets with different planet/star mass ratios $q$ (see Table 1). For
reference, $q=10^{-3}$ represents a Jupiter mass planet around a solar mass
star. This way, we can study the effect of planetary gaps and rings in the
observable properties of the disk and extract the key observables in a
computationally efficient way, avoiding the need to run expensive hydrodynamic
simulations for each combination of parameters.

Choosing the appropriate profile that mimics a planetary gap is tricky, so we
performed hydrodynamical simulations using FARGO-3D (Benítez Llambay & Masset,
2015), and we compared the effect on the observable quantities. The Kanagawa
et al. (2016) profile is an analytical approximation of the gap depth and
width but not necessarily represents the pressure bump that is caused by the
bump. Therefore we tested how strongly this assumption affects the properties
of the dust in the trap by comparing against proper hydrodynamical solutions
and disk evolution. We found that the depth of the gap is not as important to
the evolution of the disk on the SLR, but the width is the dominant factor. As
long as the planet is massive enough to create a strong pressure maximum and
thus stop the particles, the position of the pressure maximum is more
important than a precise value of the gap depth. So in summary the gap depth
is not what matters the most but the associated amplitude and location of the
pressure maximum. We should mention that the precise amount of trapping in the
bumps should still matter for e.g. planetesimal formation, but for our results
this is less relevant. We provide comparison plots and more details in
Appendix C.

We define the position of the planets in the disk $\mathrm{r_{p}}$, as a
function of the characteristic radius $r_{\mathrm{c}}$ (see Table 1). We
locate them either at $\mathrm{2/3}$, or at $\mathrm{1/3}$ of
$r_{\mathrm{c}}$. In our simulations we used none, one or two planets in these
positions. We refer the reader to subsubsection 3.1.1 for the effect of the
planet location and mass in the simulations.

###  2.4 Observables

Since the disk size is not one of the parameters that we measure directly,
using the characteristic radius $r_{\mathrm{c}}$ as a size metric is
problematic (Rosotti et al., 2019b). For this reason we define an observed
disk radius using the calculated surface brightness profile. Following
Tripathi et al. (2017) we adopt their approach to define an effective radius
($\mathrm{r_{eff}}$), as the radius which encloses a fixed fraction of the
total flux, $\mathrm{f_{\nu}}$ ($\mathrm{r_{eff}}$) = $\mathrm{xF_{v}}$. We
choose $\mathrm{x=68\%}$ of the total disk flux as a suitable intermediate
value to define $\mathrm{r_{eff}}$, as it is comparable to a standard
deviation in the approximation of a Gaussian profile.

We calculate the mean intensity $\rm{J_{\nu}}$ profile by using the scattering
solution from Miyake & Nakagawa (1993) of the radiative transfer equation

$\frac{J_{\nu}(\tau_{\nu})}{B_{\nu}(T(r))}=1-b\left(e^{-\sqrt{3\epsilon_{\nu}^{eff}}\left(\frac{1}{2}\Delta\tau-\tau_{\nu}\right)}+e^{-\sqrt{3\epsilon_{\nu}^{eff}}\left(\frac{1}{2}\Delta\tau+\tau_{\nu}\right)}\right)\,,$
(5)

with

$b=\left[\left(1-\sqrt{\epsilon_{\nu}^{eff}}\right)e^{-\sqrt{3\epsilon_{\nu}^{eff}}\Delta\tau}+1+\sqrt{\epsilon_{\nu}^{eff}}\right]^{-1}$
(6)

$\mathrm{B_{\nu}}$ is the Planck function and

$\tau_{\nu}=\left(\kappa_{\nu}^{\rm abs}+\kappa_{\nu}^{\rm
sca,eff}\right)\Sigma_{d}$ (7)

is the optical depth with $\mathrm{\kappa_{\nu}^{\rm abs}}$ the dust
absorption opacity and $\mathrm{\kappa_{\nu}^{\rm sca,eff}}$ the effective
scattering opacity which is obtained from Ricci et al. (2010) or Birnstiel et
al. (2018) (see below). As effective scattering opacity we refer to

$\kappa_{\nu}^{\rm sca,eff}=(1-g_{\nu})\kappa_{\nu}^{sca}$ (8)

where $\rm{g_{\nu}}$ is the forward-scattering parameter. $\rm{\Delta\tau}$ is

$\rm{\Delta\tau=\Sigma_{d}\kappa_{\nu}^{\rm tot}\Delta z}$ (9)

while

$\epsilon_{\nu}^{eff}=\frac{\kappa_{\nu}^{abs}}{\kappa_{\nu}^{abs}+\kappa_{\nu}^{sca,eff}}$
(10)

is the effective absorption probability. To calculate the intensity
$\rm{I_{\nu}^{out}}$ we follow the modified Eddington-Barbier approximation as
in Birnstiel et al. (2018):

$I_{\nu}^{out}\simeq\left(1-e^{-\Delta\tau/\mu}\right)S_{\nu}\left(\left(\frac{1}{2}\Delta\tau-\tau_{\nu}\right)/\mu=2/3\right)$
(11)

where $\rm{\mu=cos\theta}$

$S_{\nu}(\tau_{\nu})=\epsilon_{\nu}^{eff}B_{\nu}(T_{d})+(1-\epsilon_{\nu}^{eff})J_{\nu}(\tau_{\nu})$
(12)

is the source function.

Two-pop-py evolves only the dust and gas surface densities and the maximum
particle size111In the rest of the manuscript we often refer to ”grain size”
or ”particle size” rather than ”maximum grain size”, thereby implicitly
assuming a particle size distribution. The grain size at each radius is set by
either the maximum size possible in the fragmentation- or drift-limited
regimes, whichever is lower (see Birnstiel et al., 2012, for details). To
compute the optical properties of the dust, we therefore considered a
population of grains with a power-law size distribution, $\mathrm{n(a)}$
$\mathrm{\propto}$ $\mathrm{a^{-q}}$, with an exponent $\mathrm{q=2.5}$, for
$\mathrm{a_{min}}$ $\mathrm{\leq}$ $\mathrm{a}$ $\mathrm{\leq}$
$\mathrm{a_{max}}$. This choice follows Birnstiel et al. (2012) where the size
distribution is closer to $\rm{q=2.5}$ for disks that are in the drift limit
while for the fragmentation limited ones, a $\rm{q=3.5}$ choice would be more
suitable. Considering the disk mass is dominated by the large grains, the
choice of a smaller exponent does not alter our results significantly but it
matters for the details. Since the smooth simulations are mostly drift limited
the choice of $\rm{q=2.5}$ fits better these disks. Moreover, if a disk is
fragmentation limited then it is so mostly in the inner part but the main bulk
of the disk that defines the luminosity is in the outer part. Therefore the
luminosity will still depend mainly on the drift limited regime. The disks
with sub-structures can be fragmentation limited further out in the formed
rings but considering that these rings are mostly optically thick, the
difference between exponents is much smaller than for the smooth disks.

The grain composition consists of $\mathrm{10\%}$ silicates, $\mathrm{30\%}$
carbonaceous materials, and $\mathrm{60\%}$ water ice by volume. For a direct
comparison with observations (Tripathi et al., 2017; Andrews et al., 2018b),
we calculate the opacity in band 7 (i.e. at 850 $\mathrm{\mu m}$). Afterwards,
we use the absorption opacity to calculate the continuum intensity profile.

We also examined the effect of different opacity models and different grain
porosities. As a base model we used the composition from Ricci et al. (2010)
opacities but for compact grains (i.e. without porosity) as is Rosotti et al.
(2019a) (this model is noted as R10-0 throughout the paper). Furthermore we
used the DSHARP opacity model (Birnstiel et al., 2018) (DSHARP) and we altered
the grain porosity to $\mathrm{10\%}$ (little porous, DSHARP-10),
$\mathrm{50\%}$ (semi-porous, DSHARP-50) and $\mathrm{90\%}$ (very porous,
DSHARP-90). The particle bulk densities for the different porous grains are
$\rm{\rho_{s}}=$1.508\text{\,}\mathrm{g}\mathrm{/}\mathrm{c}\mathrm{m}^{3}$$,
$\rm{\rho_{s}}=$0.838\text{\,}\mathrm{g}\mathrm{/}\mathrm{c}\mathrm{m}^{3}$$
and
$\rm{\rho_{s}}=$0.168\text{\,}\mathrm{g}\mathrm{/}\mathrm{c}\mathrm{m}^{3}$$
respectively. An important feature of the opacity models that we used, is the
so called opacity cliff. As opacity cliff, we refer the sharp drop in the
opacity at $850\text{\,}\mu\mathrm{m}$, at a maximum particle size around
$0.1\text{\,}\mathrm{m}\mathrm{m}$, as defined in Rosotti et al. (2019a) (see
Figure 1). In all the figures that are shown in this paper, the opacity model
from Ricci et al. (2010) R10-0 is used, unless it is explicitly stated
otherwise.

Figure 1: Comparison between the opacity models that we used at 850
$\mathrm{\mu m}$ as a function of the maximum particle size and for a power-
law size distribution with an exponent of $\rm{q=2.5}$. We marked on the
figure the opacity cliff (where the opacity steeply drops by one order of
magnitude over a small range of variation in grain size) at a wavelength of
$850\text{\,}\mu\mathrm{m}$ (blue and orange shaded regions). The blue line
refers to the opacity model from Ricci et al. (2010) with compact grains
(marked as R10-0) and the orange to Birnstiel et al. (2018) with compact
grains (marked as DSHARP). Green, purple and red lines refer to
$\mathrm{10\%}$, $\mathrm{50\%}$ and $\mathrm{90\%}$ porous grains in the
DSHARP model. R10-0 and DSHARP opacity, differ by a factor of $\mathrm{\sim
8.5}$ at the position of the opacity cliff. As we increase the porosity of the
DSHARP model, we observe that the opacity cliff starts to flatten out, until
it completely disappears for very porous grains ($\mathrm{90\%}$). Moreover
the location of the cliff shifts to larger particle sizes as it diminished
with porosity. The black dashed line shows the particle size at
$\mathrm{1mm}$. The value for the R10-0 model at this size corresponds to
$\mathrm{\kappa_{\nu}^{R10-0}=$9.8\text{\,}\mathrm{c}\mathrm{m}^{2}\mathrm{/}\mathrm{g}$}$,
while for the DSHARP to
$\mathrm{\kappa_{\nu}^{DSHARP}=$4\text{\,}\mathrm{c}\mathrm{m}^{2}\mathrm{/}\mathrm{g}$}$.

###  2.5 Matching simulations

The behavior of the simulations on the size-luminosity diagram
($\mathrm{L_{mm}-r_{eff}}$ plane, hereafter SL Diagram) depends on the time
evolution of disks. According to Andrews et al. (2018b), the linear regression
of the joint data between Tripathi et al. (2017) and Andrews et al. (2018b),
gave a relation between the disk size and the
$340\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ luminosity. The effective radius
$\mathrm{r_{eff}}$ and the luminosity $\mathrm{L_{mm}}$ are correlated as

$\mathrm{log}\mathrm{r_{eff}}=(2.10^{+0.06}_{-0.03})+(0.49^{+0.05}_{-0.03})\mathrm{log}\mathrm{L_{mm}}\,,$
(13)

with a Gaussian scatter perpendicular to that scaling with a standard
deviation (1$\mathrm{\sigma}$) of $\mathrm{0.20^{+0.02}_{-0.01}}$ $\rm{dex}$
(where $\mathrm{r_{eff}}$ is in $\rm{au}$ and $\mathrm{L_{mm}}$ is in
$\rm{Jy}$ at a distance of $140\text{\,}\mathrm{p}\mathrm{c}$).

In Figure 2 (top left), we show an evolution track. It is the path that the
simulations follow on the luminosity-radius diagram in a chosen time-span.
Generically they move from the top-right (higher luminosity and radius)
regions to down-left as the disk evolves. We choose to plot the evolution
track from $300\text{\,}\mathrm{k}\mathrm{y}\mathrm{r}$ to
$3\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$. The thought behind this is that at
roughly $300\text{\,}\mathrm{k}\mathrm{y}\mathrm{r}$ our disks reach a quasi
steady state. In more detail, the drift speed in the drift limit is
$\mathrm{V_{r}=\epsilon V_{K}}$, with $\mathrm{\epsilon}$ being the dust-to-
gas ratio and $\mathrm{V_{k}}$ the Keplerian velocity. So after at least one
order of magnitude is lost in the dust, the evolutionary time becomes too
long. We refer the reader to section 4, where we show the evolution of disk
dust-to-gas ratio as a function of time for different cases. Longer
evolutionary times than the ones we are exploring here do not alter our
results significantly and we therefore exclude it to simplify the discussion.
It will be included as a topic of future research since at later stages, disks
are more strongly affected by dispersal. Moreover, our chosen time-span covers
the observed disks from the Andrews et al. (2018b) and Tripathi et al. (2017)
joint sample.

In order to filter our simulations, we divide them in categories. Every
simulation that at all times in our chosen time-span, lies within
$\mathrm{1-\sigma}$ (blue shaded region in Figure 2 (top left) of the SLR
(Equation 13) is considered as matching (see Figure 2, the green evolution
track). On the other hand, if at any time a simulation does not lie within the
area, it is considered as discrepant. The discrepant simulations can be
further divided into two sub-categories. One that is above the SLR (see Figure
2 in the top left, the purple and the yellow track) and one that is below (red
track in the same figure). With this classification, we can investigate what
are the main parameters that drive a simulation to be located on a certain
spot on the SL Diagram. It is worth mentioning that a fraction ($\mathrm{\sim
32\%}$) of the observational data points lie outside the 1-$\mathrm{\sigma}$
region by definition. This highlights that the definition of the matching
simulations is conservative with respect to the observational data.

Later on, we define the term matching fraction, as the percentage of the
matching simulations, against the total number of simulations performed with a
certain initial condition (see subsection 3.1).

##  3 Results

In this section we present the main results of this analysis. In subsection
3.1 we explain the effect of every parameter on the path of the disk on the SL
diagram. In subsection 3.2 we present the general properties that disks should
have to follow the SLR and we derive a theoretical SLR for disk with sub-
structures in subsubsection 3.2.2. In Appendix B we present an additional
analysis for the results discussed below.

###  3.1 Evolution tracks

In subsection 2.5, we explained what an evolution track is, while we showed
some examples in Figure 2 (top left). Every track is affected by the initial
conditions of the parameters chosen, and by the presence or not of a planet.
In the following sections we will explore the effect of the most important
parameters of our grid model and we will show how every parameter affects the
evolution track on the SL Diagram in Figure 2, Figure 3 and Figure 4.

Figure 2: Evolution tracks and explanation of matching and discrepant
simulations. Examples of a disk with the same initial conditions varying only
one parameter at a time. Top left: The SLR according to Andrews et al. (2018b)
and examples of evolution tracks. The black points correspond to the
observational data, and the black dashed line is the Equation 13, that shows
the relation between the luminosity and the effective radius. Finally, the
blue shaded region is the area within $1-\sigma$ of the SLR where we consider
our simulations as matching. The green track noted with the number 1 is
considered as a matching simulation, since it starts and ends inside the SLR.
On the other hand, the rest of the tracks are considered as discrepant. The
beginning of the track is where the empty bullet is (top right) and the end
where the number is printed (lower left). Top right: Varying only the presence
and the position of a Jupiter-mass planet. The red line (number 1) shows the
smooth disk, the green line (number 3), a disk with a Jupiter-mass planet at
2/3 of the $r_{\mathrm{c}}$, while the purple (number 2) corresponds to a
planet at 1/3 of the $r_{\mathrm{c}}$, and the yellow (number 4) to simulation
with two planets at the positions mentioned before. Bottom left: Varying only
the turbulence parameter $\mathrm{\alpha}$. Higher $\mathrm{\alpha}$ values
lead to higher luminosity. Bottom right: Varying only the characteristic
radius $r_{\mathrm{c}}$. Higher $r_{\mathrm{c}}$ values lead to larger and
more luminous disks.

Since the grid consists of 100.000 simulations, single evolution tracks do not
show the preferred initial conditions that allow a simulation to stay in the
SLR but only a representative case. In order to identify trends between the
initial conditions and the matching fraction of every disk, we construct
histograms, where on the y-axis we have the matching fraction (i.e. the
percentage of simulations that stay on the SLR for the chosen time-span) and
on the x-axis we have the value of each parameter in our grid model. Different
colors represent different simulation grids with or without planets and with
varying the planetary mass or positions (e.g. Figure 5). With black color we
show the simulations where we used a smooth surface density profile as in
Rosotti et al. (2019a). With green color we show the case where the
planet/star mass ratio is $\mathrm{q=1\cdot 10^{-3}}$ at a location of 1/3 of
the $r_{\mathrm{c}}$ (inner planet), with red the same planet at a distance of
2/3 of the $r_{\mathrm{c}}$ (outer planet) and with blue two planets of
$\mathrm{q=1\cdot 10^{-3}}$ at 1/3$r_{\mathrm{c}}$ and 2/3$r_{\mathrm{c}}$.
With the white hatched bars, we show the same cases but using a different
opacity model DSHARP (Birnstiel et al., 2018).

####  3.1.1 Effect of planetary parameters

In a smooth disk, the evolution track will evolve towards smaller radii and
smaller luminosity as the dust drifts inward, so the emission (and size)
decreases. Moreover, the opacity cliff moves further in as the radius in the
disk where the maximum particle size $\mathrm{a_{max}}$ is, the value for the
peak opacity decreases due to radial drift and grain growth as in Rosotti et
al. (2019b).

In contrast to this, in the case where a planet is present the pressure bump
that is formed stops the dust from drifting towards the host star delaying the
evolution of the disk on the SL Diagram, thus keeping the tracks on the SLR
for longer times. With this in mind, we expect a less extended evolution track
when we include planets, since planets are included early in the disk
evolution. In Figure 2 (top right), we show an example of the evolution tracks
of a disk with the same initial conditions, varying only the presence and the
position of a planet. The red line represents the evolution track of a disk
with a smooth surface density profile. If we include a planet with planet/star
mass ratio $\mathrm{q=10^{-3}}$ (Jupiter-mass in this case) at a location that
is close to the characteristic radius of the disk, in this case at 2/3 of the
$r_{\mathrm{c}}$, we see with the green line (number 3), that both effective
radius and the luminosity increase relative to the planet-less case, and the
evolution track is shorter at the same time-span as in all planet cases. This
is clear since the pressure bump is trapping particles and the dust mass is
retained. Therefore the luminosity is not decreasing as quickly. At the same
time, the fixed position of the pressure bump causes the effective radius to
remain the same. Together this means that the track in the SLR comes to a
halt.

On the other hand, if we place the planet close to the star, in this case at
1/3 of the $r_{\mathrm{c}}$, we observe a much longer track similar to the one
with the smooth profile as seen by the purple line. The size and the
luminosity of the disk change only slightly. The smaller radius compared to
the case where we have an outer planet is explained because the dust now stops
at the inner pressure bump and the luminosity is roughly the same for the two
cases. This could lead to the conclusion that a planet close to the star will
not affect dramatically its position on the SL Diagram, but this is not true
for all cases (see subsubsection 3.1.3). Instead, the evolution track here is
similar to the smooth one because the disk is too large and massive and the
inner planet cannot affect the evolution track too much.

As a last point, we included two planets, one at the 2/3 of the
$r_{\mathrm{c}}$ and one at 1/3 (yellow line). We observe a similar evolution
track as in the case where we have only one planet close to the outer radius,
with the disk being slightly more luminous. This is also explained by the fact
that the dust from the outer disk stops at the outer pressure bump, while the
dust that exists inside the outer planet stops at the pressure bump of the
inner planet. Since most of the dust mass initially resides outside the outer
planet, the dust that is trapped between the two planets contributes only
partially to the total luminosity. When two or more planets are present the
location of the outermost planet one is more dominant in the evolution track.

####  3.1.2 Effect of the turbulence $\mathrm{\alpha}$-parameter

The effect of the turbulence parameter $\mathrm{\alpha}$ on the evolution
tracks is straightforward, higher $\mathrm{\alpha}$ leads to higher
luminosity. In Figure 2 (bottom left), we show the evolution tracks of a disk
with the same initial conditions, varying only the
$\mathrm{\alpha}$-parameter. In this case, we choose a disk where we have also
inserted a Jupiter mass planet at the 2/3 of the characteristic radius
($r_{\mathrm{c}}$) since the effect is more prominent on these disks.

To understand the trends in Figure 2 (bottom left), it is instructive to
consider Figure 3, where we show an example of how different
$\mathrm{\alpha}$-values affect the efficiency of trapping. In the top panel
we show the dust mass local flow rate $\rm{\dot{M}_{acc,d}(r)}$
$[\mathrm{M_{\oplus}/yr}]$ as a function of radius. For low
$\mathrm{\alpha}$-values $\mathrm{1\cdot 10^{-4}}$ (red line), $\mathrm{5\cdot
10^{-4}}$ (blue) and $\mathrm{1\cdot 10^{-3}}$ (green), the mass is flowing
towards the bump and the local flow rate for $\mathrm{r<r_{p}}$ is small
($\mathrm{\leq 10^{-8}M_{\oplus}/yr}$), meaning that the trapping is efficient
enough to stop the dust from drifting towards the star. On the other hand for
the large $\mathrm{\alpha}$-values $\mathrm{5\cdot 10^{-3}}$ (yellow) and
$\mathrm{1\cdot 10^{-2}}$ (purple), the dust mass local flow rate stays almost
constant ($\mathrm{\leq 10^{-5}M_{\oplus}/yr}$) throughout the whole disk,
i.e. the bump is not trapping the particles, just locally slowing them down.
In the bottom panel we show the cumulative mass of the disk integrated from
inside out, as a function of radius. For the low $\rm{\alpha}$-values, most of
the mass is located in the bump, while for $\rm{\alpha}=\rm{5\cdot 10^{-3}}$
and $\rm{\alpha}=\rm{1\cdot 10^{-2}}$ it increases with the radius, meaning
that the bump allows more the grains to escape and therefore it does not
contain a significant fraction of the disk mass.

Returning to Figure 2 (bottom left), disks with low to medium
$\mathrm{\alpha}$-values ($10^{-4}$, $5\cdot 10^{-4}$, $10^{-3}$) are less
luminous than disks with higher $\mathrm{\alpha}$-values (red, blue and green
tracks). This is because the ring that is formed due to the pressure bump is
becoming too narrow and optically thick. The total flux emitted by an
optically thick ring of a given temperature is just a function of the emitting
area. So lower $\mathrm{\alpha}$ will lead to a narrower ring (Dullemond et
al., 2018) so less emitting area and therefore smaller luminosity,
independently of the amount of mass in the ring. On the other hand, high alpha
values work against trapping in various ways, leading to more luminous disks.
A higher $\mathrm{\alpha}$-value:

  * $\bullet$

decreases the particle size in the fragmentation limit. Smaller particles are
less efficiently trapped by radial drift (e.g. (Zhu et al., 2012)).

  * $\bullet$

increases the diffusivity which allows more grains to escape the bump.

  * $\bullet$

increases the viscosity in the same way, so more dust sizes are traveling with
the accreting gas.

  * $\bullet$

smears out the pressure peak, causing less efficient trapping by radial drift
(Pinilla et al., 2012b).

Moreover, with high $\alpha$-value, a higher dust-to-gas ratio is retained
because grain growth is impeded by fragmentation, hence radial drift is much
slower.

If we consider the two cases where the $\mathrm{\alpha}$-value is high
($5\cdot 10^{-3}$, $10^{-2}$), the planet cannot efficiently trap the dust and
the disk evolves further along the SLR to lower luminosities. As a matter of
fact, a disk that contains a planet with $\mathrm{\alpha=10^{-2}}$, behaves
the same as a smooth disk without a planet on the SL Diagram, implying that in
this case one would need a very massive planet (several Jupiter masses) to
significantly affect disk evolution. This results in the dust grains to become
smaller since the gas turbulent velocity is getting larger and the collisions
between them are more destructive. Consequently, the gap is becoming shallower
while there is also more diffusion.

Figure 3: Top panel: Local flow rate of the dust mass in
$\mathrm{M_{\oplus}/year}$ as a function of radius for a disk with a Jupiter
mass planet at $31\text{\,}\mathrm{a}\mathrm{u}$. For the low
$\mathrm{\alpha}$-values $\mathrm{1\cdot 10^{-4}}$ (red line), $\mathrm{5\cdot
10^{-4}}$ (blue) and $\mathrm{1\cdot 10^{-3}}$ (green), we observe that the
bump outside the planet location is large enough to hold the dust from
drifting towards the star, while for larger $\mathrm{\alpha}$-values
$\mathrm{5\cdot 10^{-3}}$ (yellow), $\mathrm{10^{-2}}$ (purple), there is only
a weak accumulation outside the planetary gap. Bottom panel: Cumulative dust
mass contained within a radius $r$, as function of radius. For
$\mathrm{\alpha}$-values $\mathrm{1\cdot 10^{-4},5\cdot 10^{-4},1\cdot
10^{-3}}$, roughly all the mass of the disk is inside the bump that is created
from the planet. For larger $\mathrm{\alpha}$-values $\mathrm{5\cdot
10^{-3},1\cdot 10^{-2}}$, the bump is minimal and not able to hold the dust
from drifting.

In Figure 5 (top left), we show the dependence of the
$\mathrm{\alpha}$-viscosity parameter to the matching fraction. As expected,
all simulations tend to favor low values of the turbulence parameter
$10^{-4}\leq\mathrm{\alpha}\leq 10^{-3}$. Smooth simulations show a clear
tendency towards low $\rm{\alpha}$-values because when they are drift
dominated they remain in the SLR. On the other hand sub-structured disks show
a preference towards $5\cdot 10^{-4}\leq\mathrm{\alpha}\leq 10^{-3}$. The dust
trapping is efficient enough and allows the disks to retain their mass, but
while moving to higher $\mathrm{\alpha}$-values $\mathrm{\alpha>2.5\cdot
10^{-3}}$, the trapping stops being efficient anymore leading to higher
chances that the evolution track will leave the SLR in the selected time-span.
If $\rm{\alpha}<5\cdot 10^{-4}$ the dust rings become narrow and the
luminosity is not large enough to place them in the SLR, therefore the
matching fraction decreases.

####  3.1.3 Effect of the characteristic radius - $r_{\mathrm{c}}$

In Figure 2 (bottom right), we plot the evolution tracks of a smooth disk with
the same initial conditions while varying only the characteristic radius
($\rm{r_{\mathrm{c}}=$10\text{\,}\mathrm{a}\mathrm{u}50\text{\,}\mathrm{a}\mathrm{u}100\text{\,}\mathrm{a}\mathrm{u}180\text{\,}\mathrm{a}\mathrm{u}$}$).
The same trend applies to disks where a planet is included, but for a smooth
disk the evolution tracks are longer and the effect is more easily visible.
The effect of the characteristic radius in the evolution tracks is
straightforward. The larger the $r_{\mathrm{c}}$ the more the evolution track
moves towards the top right of the plot. Meaning that for a larger disk we
expect larger luminosity. In more detail, increasing the characteristic radius
(going from red to green line) explains why the effective radius increases,
while the luminosity increases because the total disk mass remains fixed on
all these simulations. This result is consistent with the SLR from Andrews et
al. (2018b).

Taking a look at Figure 5 (middle left), we do not observe a continuous
pattern as in the other histograms. Smooth disks do not seem to depend on the
characteristic radius as disks of all sizes can reproduce the SLR. On the
other hand, sub-structured disks with small $r_{\mathrm{c}}$
($10\text{\,}\mathrm{a}\mathrm{u}$) are mostly above the correlation because
they have high luminosity relatively to their size (as it is explained by our
size luminosity estimate in subsection 4.3) and they are unable to enter the
correlation in time. For large radii ($>$150\text{\,}\mathrm{a}\mathrm{u}$$),
the disks can become too large but with low luminosity and can end up below
the correlation to the very right part of the SL-diagram (see Figure 6).

Therefore we observe a peak towards a specific characteristic radius, around
$\mathrm{80-$130\text{\,}\mathrm{a}\mathrm{u}$}$ when a planet is at a
location of 1/3 of the $r_{\mathrm{c}}$ (green color) and around
$\mathrm{30-$80\text{\,}\mathrm{a}\mathrm{u}$}$ for the case where a planet is
at a location of 2/3 of the $r_{\mathrm{c}}$ (red color). The inner planet
constrains the disk to a small size but with relatively high luminosity,
therefore placing it above the SLR before
$300\text{\,}\mathrm{k}\mathrm{y}\mathrm{r}$, while the opposite effect occurs
when an outer planet exists. When two planets are included we observe a mixed
situation of the single cases. The reason is that both pressure bumps compete
with each other and each one of them contributes with one of the ways
described above.

####  3.1.4 Effect of the disk mass - $\rm{M_{d}}$

In Figure 4 (top left), we plot the evolution tracks of a smooth disk with the
same initial conditions varying only the disk mass ($\mathrm{M_{d}}$). We
choose a smooth disk to show the effect more clearly, but the same principle
applies to the majority of disks.

The disk mass contributes to both luminosity and radius. Higher disk masses
lead to higher luminosities, both in the beginning and at the end of the
track. Since for a fixed $r_{\mathrm{c}}$, disks with higher $\mathrm{M_{d}}$
have higher $\mathrm{L_{mm}}$, it is only logical that more material will lead
to higher luminosity and vice versa. By the end of the evolution tracks, the
less massive disks have left the SLR. The dramatic curvature off the SLR for
the lowest $\mathrm{M_{d}}$ case ($\mathrm{M_{d}=5\cdot 10^{-3}M_{\star}}$) is
because all grain sizes become smaller than the opacity cliff. If we choose a
disk that contains a planet, the massive disks ($\mathrm{M_{d}\geq 5\cdot
10^{-2}}M_{\star}$) will still evolve towards lower radii and luminosities on
the SLR but the less massive ones will have shorter tracks. The pressure bump
will trap all the material outside of the planet position, so the emission and
the effective radius will both remain almost constant. The only case where the
track of a low mass disk can be long, is if the planet mass is small and the
pressure bump is not large enough to retain dust.

Figure 4: Evolution tracks with the same initial conditions varying only one
parameter at a time. Top left: Varying the disk mass of a smooth disk. Top
right: Varying the different opacity models and the porosity of the DSHARP
model of a smooth disk. Bottom left: Varying the stellar mass of disk that
contains a planet. Bottom right: Varying the fragmentation velocity of a
smooth disk.

In Figure 5 (top right), we plot the matching fraction of the disk mass
values. The tendency here is that higher disk mass leads to more simulation
inside the SLR. The reason for this is that high initial disk mass places the
disks above the SLR until they reach a stable state. While the dust is
drifting towards the host star the luminosity decreases, allowing them by our
chosen time-span, to reach the SLR and stay there for the remaining time.
Since most of the dust is in the trap at this point, the remaining evolution
time is set by the trap life time.

The difference is noticeable between the smooth and the planet(s) case. As we
see from the yellow bars, a disk with a smooth surface density profile must be
initially massive ($\rm{M_{d}\geq$0.025\text{\,}\mathrm{M}_{\star}$}$) to
remain in the SLR. The probability of a smooth simulation to match is even
larger than the cases where a planet is included.

Some of these results though are an effect of the opacity model used and the
chosen time-span, as we will discuss in subsubsection 3.1.5.

####  3.1.5 Effect of different opacity/grain porosity

In Figure 4 (top right), we are showing representatively the behavior of
several similar tracks varying only the opacity model. The simulations with
the Ricci et al. (2010) (R10-0, blue line) opacity model produce more luminous
and larger disks than the ones with the DSHARP (orange line) model due to the
higher value of the opacity. The R10-0 opacity is 8.5 times higher than the
DSHARP at the peak of the opacity cliff. If we use slightly porous grains
(DSHARP-10, green line) by altering the DSHARP opacity, we observe that the
effect is insignificant, as the shape of the opacity is roughly the same. On
the contrary, for semi and very porous grains (DSHARP-50 and DSHARP-90, purple
and red line) the opacity cliff starts to flatten out (Kataoka et al., 2014),
leading to a disk with low luminosity and no significant change in disk size.

On all the histogram figures (Figure 5) the same trend stands for either the
opacity from R10-0 Ricci et al. (2010) (solid color bars) or the DSHARP
(Birnstiel et al., 2018) (hatched bars). The difference is that more
simulations match when the R10-0 opacity model is used as opposed to the
DSHARP model. Disks with the DSHARP opacity are generally less bright because
they do not become optically thick in the rings and end up below the SLR.
Therefore it would need more dust (i.e. stronger traps) for them to be
luminous enough. Especially for the smooth case (yellow bars), there is only a
few simulations that match, hence the hatched bars are barely visible. For
smooth disks the total matching fraction is $\rm{29.6\%}$ with the R10-0,
while it is $\rm{0.8\%}$ with the DSHARP opacity. For disks with an inner
planet the matching fractions are $\rm{30.2\%}$ and $\rm{15.9\%}$
respectively. We refer the reader to subsection 4.3, where we explore the
overall impact of porous grains for the entire grid of models.

Figure 5: Histograms of the matching fraction for disk mass, characteristic
radius, stellar mass and fragmentation velocity. The matching fraction shows
the percentage of the simulations that remained on the SLR for the chosen
time-span
($\rm{$300\text{\,}\mathrm{k}\mathrm{y}\mathrm{r}$-$3\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$}$).
Top left: Dependence on the $\alpha$-value. There is a preference to low
$\alpha$-values $\rm{(10^{-4}\leq\alpha\leq 10^{-3})}$. Top right: Dependence
on the disk mass. There is a preference to high disk masses
$\mathrm{(0.025\leq\frac{M_{d}}{M_{\star}}\leq 0.25)}$. Middle left:
Dependence on the characteristic radius. Smooth disks do not depend on the
$r_{\mathrm{c}}$. Middle right: Dependence on the stellar mass. There is a
small preference to larger values. Bottom left: Dependence on the
fragmentation velocity. There is a small preference to larger values.

####  3.1.6 Effect of the stellar mass - $\mathrm{M_{\star}}$

In our models, the stellar mass is assumed to be directly correlated with the
disk mass. This is due to the fact that we varied the stellar mass, but we
always kept the disk-to-star mass ratio constant. Therefore, a larger star
mass implies a larger total mass of the dust, leading to larger continuum
luminosities.

In Figure 4 (bottom left), we plot the evolution tracks of a disk that
contains a planet with the same initial conditions, varying only the stellar
mass for
$\mathrm{M_{\star}=$0.2\text{\,}\mathrm{M}_{\odot}0.6\text{\,}\mathrm{M}_{\odot}1.0\text{\,}\mathrm{M}_{\odot}2.0\text{\,}\mathrm{M}_{\odot}$}$.
As expected, the largest value of $\mathrm{M_{\star}=2.0M_{\odot}}$ leads to
the largest and most luminous disk (green line), while the opposite is
happening for $\mathrm{M_{\star}=0.2M_{\odot}}$ (red line). There is similar
behavior for smooth disks, but in that case the radius of each disk will be
much smaller due to radial drift.

Furthermore, the stellar mass is the least important parameter on defining if
a simulation matches or not. The histogram in Figure 5 (middle right) confirms
this. Even though the trend shows that using higher stellar mass has higher
matching fraction, this is because the stellar mass scales to the disk mass
and higher disk mass leads to more matching simulations (see subsubsection
3.1.4).

This scaling implies that the luminosity ($\mathrm{L_{mm}}$) scales with the
stellar mass ($\mathrm{M_{\star}}$). Our models follow a relation that is not
as steep as the observed $\mathrm{L_{mm}\propto M_{\star}^{1.5}}$ in Andrews
et al. (2018b), because there is not a correlation between disk size and
stellar mass in our simulations. We refer the reader to subsection 4.5 and
Figure 10 where we explore further this relation.

####  3.1.7 Effect of the fragmentation velocity - $\mathrm{v_{frag}}$

In Figure 4 (bottom right), we plot the evolution tracks, of a smooth disk
with the same initial conditions, varying only the fragmentation velocity for
values of
$200\text{\,}\mathrm{c}\mathrm{m}\mathrm{/}\mathrm{s}600\text{\,}\mathrm{c}\mathrm{m}\mathrm{/}\mathrm{s}1000\text{\,}\mathrm{c}\mathrm{m}\mathrm{/}\mathrm{s}2000\text{\,}\mathrm{c}\mathrm{m}\mathrm{/}\mathrm{s}$.
We observe that for medium and large values of $\mathrm{v_{frag}}$ (in this
case for
$\mathrm{v_{frag}}\geq$600\text{\,}\mathrm{c}\mathrm{m}\mathrm{/}\mathrm{s}$$),
the evolution tracks overlap. Since most of our simulations are drift limited
we expect that no effect from the fragmentation velocity will take place for
these values since particles do not grow big enough to drift, so more mass
remains at large radii to produce more emission. Therefore, this effect only
arises when the fragmentation velocity is becoming too low leading for the
first snapshots to higher luminosity.

Moreover, if a disk is fragmentation limited then it is so mostly in the inner
part. Therefore, considering that the main bulk of the disk is in the outer
part, the emission that defines the luminosity will still depend on the drift
limited regime hence leading to the overlapping tracks. The effect of a planet
in these tracks is minimal and we expect a similar behavior to the case shown
here.

This can be validated in Figure 5 (bottom left), where we plot the matching
fraction compared to the fragmentation velocity and the tendency is to higher
fragmentation velocities for all cases. More specifically, if
$\mathrm{v_{frag}\geq 600cm/s}$ there is a large number of matching
simulations and it only gets larger with increasing fragmentation velocity. In
this range, the simulations are mostly drift limited. For low values of
$\mathrm{v_{frag}}$, most of the simulations are fragmentation limited and
they lose luminosity relatively quickly, moving them out of the SLR. Low
values of $\mathrm{v_{frag}}$ lead to smaller particles and less efficient
trapping. Therefore, those disks lose their solids too quickly. This is
analyzed in more detail in Appendix B.

For reference, following recent lab experiments (Blum, 2018),
$$100\text{\,}\mathrm{c}\mathrm{m}\mathrm{/}\mathrm{s}$\leq\mathrm{v_{frag}\leq$1000\text{\,}\mathrm{c}\mathrm{m}\mathrm{/}\mathrm{s}$}$
is considered as a value consistent with lab work and
$\mathrm{v_{frag}>$1000\text{\,}\mathrm{c}\mathrm{m}\mathrm{/}\mathrm{s}$}$ as
high fragmentation velocity.

###  3.2 Heat-maps

Single evolution tracks give us an idea of how a single simulation evolves on
the SLR but with a large sample of simulations as the one we have, it is not
easy to extract global results. Since there are too many tracks to over-plot
all of them on a single diagram, we treat the position of every simulation at
every snapshot as an independent sample and we plot them on a heat map. In
these figures, we plot the position of every simulation for a specific case
(smooth and planet in different locations), for three different snapshots
($300\text{\,}\mathrm{k}\mathrm{y}\mathrm{r}$,
$1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$,
$3\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$).

In Figure 6, we plot three different cases. In the first column we show the
smooth case, in the second one the one where we use a planet/star mass ratio
$\mathrm{q=10^{-3}}$ at 1/3 of the $r_{\mathrm{c}}$ and in the last one
$\mathrm{q=10^{-3}}$ at 2/3 of the $r_{\mathrm{c}}$. Each row (from top to
bottom) refers to the snapshots, $300\text{\,}\mathrm{k}\mathrm{y}\mathrm{r}$,
$1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ and
$3\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ respectively. The solid white line
refer to the SLR from Andrews et al. (2018b) and the dashed white lines show
the $\mathrm{1\sigma}$ deviation (same as the blue shaded region on the single
evolution tracks). The color-bar to the right shows the number of simulations
in a single cell. The red line is our prediction for the cases where we
include a planet (see subsubsection 3.2.2). Instead of following the relation
from Andrews et al. (2018b), they seem to follow a relation of
$\mathrm{L_{mm}\propto r_{eff}^{5/4}}$. In subsubsection 3.2.1 and
subsubsection 3.2.2 we perform a more detailed analysis on the this topic.

We observe that on all cases, most of the disks start inside and above the
correlation (first row at $300\text{\,}\mathrm{k}\mathrm{y}\mathrm{r}$). In
the smooth case, the disks lose a lot of their luminosity relatively quickly
and also shrink in size, making them move to the lower radii due to radial
drift. We end up with a large number ($\rm{29.6\%}$) of simulations occupying
and following the SLR. As we explain later subsubsection 3.2.2, the slope is
expected from Rosotti et al. (2019a) but the normalization depends on the
choice of opacities.

On the other hand, when we include a planet and the time increases the disks
are not decreasing in size, but mainly in luminosity. This is due to the
formation of the pressure bump that keeps the dust from drifting further in
keeping the same effective radius. The consequence is that if the disks leave
the SLR, it is due to luminosity decrease since they move vertically on the
diagram. The clustering in $\mathrm{r_{eff}}$ that is formed on the plots with
a planet are an artifact of our parameter grid. A randomly chosen value of
$r_{\mathrm{c}}$ and planet position would result in a continuous (non-
clustered) distribution. From the comparison of the two cases where a planet
is included, there are more matching simulations when a planet is in the inner
part of the disk ($\rm{30.2\%}$). Having a planet in the outer part leads the
disks to the right (large radii) and bottom part of the diagram and
consequently leaves them outside of the relation ($\rm{20.6\%}$ of the disks
match).

The difference becomes striking when we use the dust opacities from DSHARP in
Figure 7. Since the DSHARP opacity is lower than the R10-0 (see Figure 1),
many of the smooth disks start below the SLR. This leads the majority of them
outside of the relation by the last snapshot
($3\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$) and only ($\rm{0.8\%}$) of the
disks match. The same stands for the case where a planet is included. The
disks have lower luminosity in the first snapshot but the pressure bump formed
is big enough to keep them in the relation for the remaining time
($>\rm{11.1\%}$) of the disks match depending on the model.

Similar behaviour for the cases where we include strong sub-structures stands
for the opacity model with semi-porous grains, DSHARP-50, in Figure 8. Even
though there are less matching simulations in total, if the sub-structures are
large enough, they are able to keep a significant amount of simulation is the
clusters on the SLR ($>\rm{10.7\%}$). The same argument cannot be made for
simulations with the smooth surface density profile. With semi-porous grains
there is only a small fraction of matching simulations ($\rm{0.7\%}$). The
absence of a strong opacity cliff in the opacity profile leads to low
luminosities and consequently all the simulations below the SLR. An almost
identical point stands looking at the heat-map (Figure 19) for the case for
the case where very porous grains are used (DSHARP-90). The complete absence
of the opacity cliff does not allow a considerable fraction of smooth disks to
enter the SLR ($\rm{1.3\%}$) while a similar fraction of sub-structured disks
match as in the DSHARP-50 case ($>\rm{10.2\%}$). This heat-map is included in
Appendix D.

From these heat-maps we can extract three important results.

  * $\bullet$

Disks with strong traps (i.e. massive planets) follow a different SLR than
smooth disks, while smooth disks are more consistent with the SLR in terms of
the shape of the relation.

  * $\bullet$

Whether a smooth disk matches or not depends heavily on the opacity model.
Birnstiel et al. (2018) DSHARP opacities produce significantly less
simulations in the SLR than the Ricci et al. (2010) R10-0 model and only a
fraction of simulation will match with semi and very porous grains and the
model DSHARP-50 and DSHARP-90. Therefore, the porosity should be smaller than
$\mathrm{50}\%$ when the Birnstiel et al. (2018) opacities are used. However,
the distribution of simulations is significantly tighter than the observed
correlation for the smooth disks with the R10-0 opacity. As it is discussed
later (section 4), the observed correlation can be a mixture of smooth and
sub-structured disks that adds scatter to the simulated SLR.

  * $\bullet$

A bright disk (top right on the SL diagram) is more probable to remain in the
SLR if there is a pressure bump formed in the first $\mathrm{1Myr}$,
regardless of the opacity model.

Figure 6: Heat-maps of simulations with the Ricci et al. (2010) opacities. The
three different columns represent the smooth, planet at $\mathrm{1/3r_{c}}$
and a planet at $\mathrm{2/3r_{c}}$ respectively. The rows represent three
different snapshots at $\mathrm{300kyr}$, $\mathrm{1Myr}$ and $\mathrm{3Myr}$
respectively. The white solid is the SLR from Andrews et al. (2018b) and the
red solid line our fit for the cases where a planet is included. The color-bar
shows the number of simulations in a single cell. The blue dash-dotted line
shows the minimum limit
($\mathrm{r_{eff}\sim$10\text{\,}\mathrm{a}\mathrm{u}$}$) where observational
results are available. Figure 7: Heat-maps of simulations with the Birnstiel
et al. (2018) opacities. The three different columns represent the smooth,
planet at $\mathrm{1/3r_{c}}$ and a planet at $\mathrm{2/3r_{c}}$
respectively. The rows represent three different snapshots at
$\mathrm{300kyr}$, $\mathrm{1Myr}$ and $\mathrm{3Myr}$ respectively. The white
solid is the SLR from Andrews et al. (2018b) and the red solid line our fit
for the cases where a planet is included. The color-bar shows the number of
simulations in a single cell. The blue dash-dotted line shows the minimum
limit ($\mathrm{r_{eff}\sim$10\text{\,}\mathrm{a}\mathrm{u}$}$) where
observational results are available. Figure 8: Heat-maps of simulations with
the Birnstiel et al. (2018) D-50 opacities with 50% porostiy. The three
different columns represent the smooth, planet at $\mathrm{1/3r_{c}}$ and a
planet at $\mathrm{2/3r_{c}}$ respectively. The rows represent three different
snapshots at $\mathrm{300kyr}$, $\mathrm{1Myr}$ and $\mathrm{3Myr}$
respectively. The white solid is the SLR from Andrews et al. (2018b) and the
red solid line our fit for the cases where a planet is included. The color-bar
shows the number of simulations in a single cell. The blue dash-dotted line
shows the minimum limit
($\mathrm{r_{eff}\sim$10\text{\,}\mathrm{a}\mathrm{u}$}$) where observational
results are available.

####  3.2.1 Width of the pressure maxima

To understand the overall shape of the heat-map for the case with massive
planets (i.e. the red lines in Figure 6 and Figure 7), we will derive a
theoretical estimate in the following. This estimate depends on the width and
position of the pressure maximum formed outside the position of the gap-
opening planets. We therefore first derive a relation of the gas width vs
radius $r$, planet/star mass ratio $q$ and scale height $h$, using
hydrodynamical simulations of planet-disk interaction with the FARGO-3D code
(Benítez Llambay & Masset, 2015). In the subsequent subsubsection 3.2.2, we
estimate the SLR based on these empirically determined widths. For a complete
derivation of both sections, we refer the reader to subsection A.1 and
subsection A.2.

In addition to our two-pop-py models, we performed 24 hydrodynamical
simulations with FARGO-3D (Benítez Llambay & Masset, 2015) for different
planet/star mass ratios, planet locations and $\mathrm{\alpha}$-values (see
Table 3 in Appendix C). We used these simulations to calculate the width of
the outer pressure bump caused by the planet in the gas. The surface density
maximum is locally well fitted by a Gaussian which allows us to measure the
width (i.e. the standard deviation) using the curvature at the maximum. By
measuring all the widths of our hydrodynamical simulations we fit as a
multiple power-law to search how they scale with radius, scale height and the
$\mathrm{\alpha}$-parameter:

$\sigma_{g}=C\cdot h^{p}\cdot q^{k}\cdot\alpha^{l}\,$ (14)

where $\mathrm{C}$ is a constant, $h$ is the scale height, $q$ the mass ratio
and $\mathrm{\alpha}$ the turbulence parameter. We found that the width in the
measured range scales approximately as:

$\sigma_{g}\propto h^{0.81}\cdot q^{0.14}\cdot\alpha^{0.05}\,$ (15)

####  3.2.2 Size-luminosity relation of disks with companions

In Figure 6 we show a red line that scales differently from the SLR when we
include planets. We predict that this line fits with a correlation
$\mathrm{L_{mm}\propto r_{eff}^{5/4}}$. If we assume that all the luminosity
of a disk comes from rings that are approximately optically thick, we can
approximate:

$L\simeq A\cdot B_{\nu}\,$ (16)

where $\mathrm{A}$ is the area and $\mathrm{B_{\nu}}$ is the Planck function.
If we assume the Rayleigh Jeans approximation to approximate the Planck
function with the temperature, the equation becomes:

$L\simeq A\cdot T\,$ (17)

We make the assumption that the area of the pressure bump scales as
$\mathrm{A}\propto r\cdot\sigma_{d}$, where $\mathrm{r}$ is the radius and
$\mathrm{\sigma_{d}}$ is the width of the pressure bump in the dust and there
is linear scaling of $\mathrm{\sigma_{d}}$ with $\mathrm{h}$. The width of the
dust ring depends on the width of the gas

$\sigma_{d}\propto\sigma_{g}\cdot\sqrt{\frac{\alpha}{St}}\,$ (18)

as in (Dullemond et al., 2018), where $\mathrm{St}$ is the Stokes number,
while in the ring there is also effective dust trapping and diffusion. Using
the relation that previously measures the gas width subsubsection 3.2.1, we
find that the luminosity ($\mathrm{L_{mm}}$) scales with the radius as:

$L_{mm}\propto r_{eff}^{5/4}$ (19)

which is the relation that we did plot with the red line in Figure 6, Figure
7, Figure 8 and Figure 19. We find that this theoretical estimate well
explains the size luminosity scaling seen when strong sub-structure is
present. However, towards larger radii, this relation slightly over-predicts
the luminosity. For example, one could get a shallower slope looking at the
heat-map, because towards large radii our fitting line is above the main bulk
of the simulations. For a complete derivation of both sections, we refer the
reader to subsection A.1 and subsection A.2.

##  4 Discussion

To summarize the results discussed above, we have explored the observed trend
among the (sub-)mm disk continuum luminosity ($\mathrm{L_{mm}}$) and the
$\mathrm{68\%}$ effective radius ($\mathrm{r_{eff}}$) of protoplanetary disks.
Following the size-luminosity relation (SLR) obtained from Tripathi et al.
(2017) and Andrews et al. (2018b), $\mathrm{L_{mm}\propto r_{eff}^{2}}$, we
have shown which initial conditions are favorable for a disk to remain on the
SLR for a time-span of $300\text{\,}\mathrm{k}\mathrm{y}\mathrm{r}$ \-
$3\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$. We explored the effect of every
parameter on the disk evolution tracks on the SL Diagram, we got a visual
representation of how the disk population moves on the same diagram and we
found relations between the parameters (Appendix B). We present a different
correlation for disks that are dominated by strong sub-structures compared to
the disks that have a monotonically decreasing surface density profile.
Moreover we investigated the effect of different opacity models with compact
or porous grains and we conclude that it is a major factor for reproducing the
observational results. In the following sections we briefly recap these
results and we discuss in detail some of the implications.

###  4.1 Dominant parameters

In summary, our results imply that the most dominant parameters for the
evolution of disks are the viscosity parameter $\mathrm{\alpha}$, the initial
disk mass $\mathrm{M_{d}}$, the location of a giant planet if present, and the
opacity model that is used to derive the continuum intensity. The disks that
match the SLR are characterized by low turbulence ($\mathrm{\alpha}\leq
10^{-3}$), high disk mass ($\mathrm{M_{d}}\geq 2.5\cdot 10^{-2}M_{\star}$) and
are affected strongly by the existence of a strong trap (in this study caused
by a giant planet).

Turbulence-$\mathrm{\alpha}$ values larger than $10^{-3}$ lead to smaller
grains due to fragmentation and consequently to less luminous disks that do
not enter the SLR. Moreover, particles are diffused more efficiently, and due
to their size less efficiently trapped. Finally, the dust trap is not as
pronounced if the alpha viscosity is higher. All of this is acting in concert
to make dust trapping ineffective, and causes the disks to behave as if they
were smooth (see Figure 3 in section 3). If the fragmentation velocity is
large enough, there can be simulations that stay on the SLR, but we consider
these disks to have unrealistic initial conditions according to the known
literature.

High initial disk mass locates the disks initially either inside or above the
SLR (i.e. too bright for the given size) until they reach a quasi steady
state. This allows them to migrate to smaller luminosities while they evolve
up to $3\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ and still remain in the SLR.
This effect is aided by the right choice of opacity model and grain porosity.
Compact grains shift the position of disks to larger luminosity (see
subsection 4.3). On the other hand, most of the less massive and smaller disks
end up below the correlation at smaller luminosities, characterizing them as
discrepant.

Planets can alter the evolution path of the disk on the SL Diagram
significantly. An effectively trapping planet causes the disk to quickly
settle to a quasi-steady state on the SL Diagram, thereby leading to a shorter
track and thus delaying the evolution of disks towards lower luminosity and
radius. Disks with a massive planet in the inner part of the disk
($\mathrm{1/3r_{\mathrm{c}}}$) have more extended evolutionary track that are
overall less luminous. In contrast disks with a planet in the outer part
($\mathrm{2/3r_{\mathrm{c}}}$) have shorter and more luminous disks if all the
other parameters remain the same. This can be explained by the planet trapping
a large part of the disk solids at large radii. When two planets are included
then both of them contribute to the luminosity while the outer one defines the
effective radius of the disk. Overall, the existence of planets increases the
fraction of matching simulations on the SLR but this result is also a function
of the opacity (see subsection 4.3).

###  4.2 Position along the SLR

The position of a disk along the SLR is determined mainly by the disk mass
$\mathrm{M_{d}}$ and the disk size $r_{\mathrm{c}}$, as it can be seen in
Figure 2 and Figure 4. More massive and large disks are located on the top
right part of the SL diagram while small and less massive, to the middle and
left part. In Figure 9, we show the kernel-density-estimate (kde) of the
luminosity for all matching simulations for four different cases and three
different snapshots, while with black color we plot the observational kde from
Andrews et al. (2018b). The brightest disks that stay on the SLR are the ones
that contain planets which are located at the outer part of the disk (2/3 of
the $r_{\mathrm{c}}$, yellow and green line). A planet in the outer region
leads to larger / more luminous disks as explained in subsection 2.3 and
subsection 4.1. Massive planets at 1 and
$3\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ reproduce the peak of the observed
brightness distribution but overall produce too many bright and too few faint
disks. The peak of the luminosity distribution for smooth disks is generally
at much lower luminosities. Given these results, it is conceivable that the
observed sample consists of two distinct categories of disks: a brighter and
larger category due to massive outer planets that trap the dust while planets
in the second category are not massive enough to trap the dust effectively and
those disks evolve similar to a smooth disk.

Figure 9: Kernel-density-distribution of the luminosity for all matching
simulations using the Ricci et al. (2010) R10-0 opacity model, for four
different cases and three different snapshots, from $\mathrm{300kyr-3Myr}$.
The black line refers to the disks from the Andrews et al. (2018b) sample.
Disks with planets have larger luminosity, while smooth disks have low
luminosity at $\mathrm{3Myr}$. When two planets are included the luminosity is
larger than with a single planet.

###  4.3 Is there a preferable opacity model?

We used different opacity models for our study. As a base model, we made use
of the composition of Ricci et al. (2010) opacities but for compact grains
(i.e. no porosity) as in Rosotti et al. (2019a) (this model is noted as R10-0
throughout the paper). Moreover we used the non-porous Birnstiel et al. (2018)
opacities (DSHARP) and varied the porosity between $10\%$ DSHARP-10, $50\%$
DSHARP-50 and $90\%$ DSHARP-90.

Therefore, we find that independent of the model used, relatively compact
grains ($\mathrm{<50\%}$) are preferred instead of the highly porous grains.
When compact grains are included the initial position of the disks on the SL
Diagram is shifted towards higher luminosity giving it more time to evolve in
the SLR on our chosen time-span. Disks with the DSHARP opacity are generally
less bright and end up below the SLR. We remind the reader that the opacity at
our wavelength is a factor of $\mathrm{\sim 8.5}$ larger in the R10-0 case
compared to DSHARP at the opacity cliff location ($\mathrm{\sim 0.1-1mm}$),
with the difference mainly stemming from the choice of carbonaceous material
(Zubko et al., 1996 versus Henning & Stognienko, 1996, see the comparison in
Birnstiel et al., 2018). However, this point holds only for smooth disks and
disk with weak sub-structures (where a disk behaves as smooth). If this is the
case then only compact grains can explain the SLR while when sub-structures
are strong then any of the opacity models and the porosities tested in this
work can explain the sub-structured SLR. The latter applies because most of
the sub-structures become optically thick.

One can argue that alternative compositions that also exhibit a strong opacity
cliff and a high opacity would be equally suitable.

###  4.4 Which disks populate the SLR?

We have shown that when strong traps (i.e. massive planets) are included, then
disks follow a different SLR than the smooth ones. Based on measurements of
the width of the pressure maxima formed by the planet in hydrodynamical
simulations (see subsection A.1), we derived a theoretical prediction for
disks with sub-structures (subsubsection 3.2.2). Smooth disks with compact or
slightly porous grains seem to follow the Andrews et al. (2018b) relation
$\mathrm{L_{mm}\propto r_{eff}^{2}}$, while disks with massive planets follow
a relation $\mathrm{L_{mm}\propto r_{eff}^{5/4}}$.

This result does not imply that the observed disks from the Tripathi et al.
(2017) sample are all free of substructure, but that they might not show
strong enough and optically thick sub-structures as, e.g. AS209 or HD163296
(see Huang et al., 2018a). Figure 6 shows how smooth disks follow the SLR,
while disks with strong sub-structures follow a different relation but
intersecting the SLR at the bright end (top right part of the SLR). In
contrast, the less luminous disks follow the SLR if they are smooth, but disks
of the same effective size with substructure are too luminous (cf. bottom left
part of the SLR).

In Hendler et al. (2020, Fig. 6), disks seem to follow a universal relation
(close to the SLR) in all star-forming regions (Ophiuchus, Tau/Aur, Lupus,
Chal) except from USco, the oldest region ($\mathrm{\sim 10Myr}$). The
observed SLR therefore might flatten with age of the region. We could examine
these results since we evolve our simulations to $\rm{10Myr}$, but our models
do not include photo-evaporation and that would lead to uncertainties in the
results. For example, at the age of USco the detectable disk fraction is
$\rm{<20\%}$, while in our models it would be $\rm{100\%}$.

This raises a question: If a planet is not massive at early times, but around
$1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ has a planet/star mass ratio of
$\mathrm{q=10^{-3}}$, will the disk follow the observed SLR or the SLR with
strong sub-structures? According to the analysis on the evolution tracks in
subsection 3.1, most of the small and less luminous disks that do not
initially have a giant planet, drift towards lower radii and luminosities and
even below the SLR. Therefore strong sub-structures need to form in the first
$\sim$$1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ for the disk to follow the
$\mathrm{L_{mm}\propto r_{eff}^{5/4}}$ relation. This result might imply that
in most of the star-forming regions, strong sub-structures might not have
formed early enough for small disks, or that the sub-structure is weak. On the
other hand, bright and large disks can very well show strong sub-structures
and follow the SLR at the same time. This is indeed the case for the DSHARP
sample (Andrews et al., 2018a) that is biased towards bright disks and show
significant sub-structures in every source. The latter can be confirmed from
Figure 9 as in the previous section. The brightest disks that stay on the SLR
are the ones that contain planets that are located in the outer part of the
disk (yellow and green line).

The SLR can be explained if there is a mixture of both smooth and strong sub-
structured disks. Smooths disks always follow the SLR as shown already from
Figure 6, while the bright sub-structured disks populate the upper right part
of the SLR (Figures 6, 7, and 19). Disks with sub-structures that have large
$r_{\mathrm{c}}$ and low disk mass $\mathrm{M_{d}}$ populate the lower right
part of the plot under the SLR. These disks are not favored by Andrews et al.
(2018b) who finds a tentative positive correlation between the mass of the
star (or the disk) and the size of the disk. If massive, small disks are
excluded from the plot then the SLR could be reproduced by both sub-structured
disks that occupy the upper right part and smooth (or weakly sub-structured)
disks that occupy the lower left part of the SLR. Our results seem to be in
agreement with the observational classification from van der Marel & Mulders
(2021) where they suggest that all bright disks should have sub-structures
formed by giant planets. Moreover, the SLR for the sub-structured disks is
independent of the opacity model but it slightly over-predicts the luminosity
for the very large disks (see subsection 4.3 for more about the opacity).

###  4.5 $\mathrm{L_{mm}-M_{\star}}$ relation

In subsubsection 3.1.6 it is discussed that the stellar mass
($\mathrm{M_{\star}}$) is directly correlated with the disks mass
($\mathrm{M_{d}}$) and the disk temperature is only a weak function of the
$\mathrm{L_{\star}}$ (and therefore $\mathrm{M_{\star}}$). The fact that the
disk mass scales with the stellar mass implies that the luminosity
($\mathrm{L_{mm}}$) scales with the stellar mass. In Figure 10 the
$\mathrm{L_{mm}}-M_{\star}$ relation is shown for three different models for
all matching simulations. Smooth case (yellow lines), a planet with
planet/star mass ratio $\mathrm{q=10^{-3}}$ at $\mathrm{1/3r_{\mathrm{c}}}$
(green lines) and a planet with the same mass ratio at
$\mathrm{2/3r_{\mathrm{c}}}$ (red lines). The markers define the median value
of the luminosity at $1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ and the error
bars are the $\mathrm{75\%}$ percentile from the upper and lower value. The
blue line is the $\mathrm{L_{mm}\propto M_{\star}^{1.5}}$ correlation from
Andrews et al. (2018b), a correlation that is consistent with those found from
previous continuum surveys of comparable size and age (Andrews et al., 2013;
Ansdell et al., 2016; Pascucci et al., 2016). For any of our models, the
correlation is not as steep as the Andrews et al. (2018b), but the cases with
strong sub-structures have steeper profile than the smooth one. The reason is
that no correlation between disk size and stellar mass was imposed in the
parameter grid. If a size-mass correlation as inferred by Andrews et al.
(2018b) was imposed, the mass-luminosity relation is expected to steepen as
disks with optically thick sub-structures would be larger and therefore
brighter. However, reproducing the observed mass-luminosity trend will be part
of a future population synthesis study. A similar manifestation of the same
trend can be seen in Figure 13. In the panel $\mathrm{r_{\mathrm{c}}-M_{d}}$,
marked with white dots is the mean value of the characteristic radius for
every disk mass. In order for the correlation to reproduce the observations,
more massive disks should have initially been larger. In other words large,
low-luminosity disks would be expected in the lower right of the SL diagram,
but are not observed. Preliminary results indicate that these disks need to
have low disk mass $\rm{(M_{d}<10^{-2}M_{\star})}$ and be large in size
$(\rm{r_{\mathrm{c}}>$150\text{\,}\mathrm{a}\mathrm{u}$})$. Moreover the
turbulence parameter should be relatively small $\alpha$ $\rm{\leq 10^{-3}}$,
otherwise the disk will behave as smooth and will follow the SLR.

Figure 10: $\mathrm{L_{mm}}-M_{\star}$ relation at $\mathrm{1Myr}$ for three
different cases for the matching simulations. Smooth case (yellow lines), a
planet with planet/star mass ratio $\mathrm{q=10^{-3}}$ at
$\mathrm{1/3r_{\mathrm{c}}}$ (green lines) and a planet with the same ratio at
$\mathrm{2/3r_{\mathrm{c}}}$ (red lines). The points define the median value
of the luminosity at $\mathrm{1Myr}$ and the error bars are the
$\mathrm{75\%}$ percentile from the upper and lower value. The blue line is
the $\mathrm{L_{mm}\propto M_{\star}^{1.5}}$ correlation from Andrews et al.
(2018b).

###  4.6 Scattering

Scattering is included in our simulations as introduced in subsection 2.4.
Compared to the case where only the absorption opacity is used, the difference
is minimal and it can be observed only in few cases. With the inclusion of
scattering, the originally brightest disks (above the SLR) tend to move
towards lower luminosity (move down in the SL diagram). This happens for disks
that are optically thick hence for the ones that contain planets. This effect
favors the SLR and allows slightly more disks ($\rm{\sim 2\%}$) to enter in
the selected region. However, for moderately optically thick disks the
emission is larger and a small fraction of disks moves up (towards larger
luminosity) on the SL diagram (Figure 4 in Birnstiel et al. (2018)). This is
happening because the derived intensity (Equation 11) does not saturate to the
Planck function, but to a slightly smaller value for a non-zero albedo. This
is the well-known effect that scattering makes objects appear cooler than they
are in reality. On the other hand, for small optical depths ($\rm{\tau<<1}$)
the effect of scattering is insignificant because the intensity (Equation 11)
approaches
$\rm{I_{\nu}^{out}\longrightarrow\epsilon_{\nu}^{eff}B_{\nu}(T_{d})\Delta\tau/\mu}$
which is the identical solution as when $\rm{\kappa_{\nu}^{sca}}$ is set to
zero while $\rm{\kappa_{\nu}^{abs}}$ is kept unchanged as shown also in
Birnstiel et al. (2018).

The effect of scattering depends also on the albedo
($\rm{\eta_{\nu}=1-\epsilon_{\nu}^{eff}}$). For the compositions we use the
maximum effective albedo is $\rm{0.57}$ for R10-0 and $\rm{0.82}$ for DSHARP
while it can reach up to $\rm{\sim 0.97}$ for DSHARP-90. For these
compositions the effect of scattering is never more than a factor of $\rm{\sim
1.7}$ at a particle size of $\rm{1mm}$. The result of scattering is higher if
the we increase albedo but by using a plausible composition, the result is
effectively negligible.

However we obtain different results compared to Zhu et al. (2019). In that
work is discussed that a completely optically thick disk with high albedo
($\rm{0.9}$) can be constructed, which therefore lies along the SLR with the
right normalization (because the high albedo and high optical depth lowers the
luminosity). However our findings show that we cannot reach those results from
an evolutionary perspective. For smooth disks the dust drifts to the inner
part of the disk and the disk is no longer optically thick. On the other hand,
disks with sub-structures create only rings rather than disks that are
completely optically thick everywhere.

###  4.7 Prediction for longer wavelength

A recent study (Tazzari et al., 2021), showed a flatter SLR at $\rm{3.1mm}$
($\rm{Lmm\propto r_{eff}^{1.2}}$), confirming that emission at longer
wavelengths becomes increasingly optically thin.

We performed a series of simulations at $3.1\text{\,}\mathrm{m}\mathrm{m}$ to
compare with these results. Disks are fainter and smaller at
$3.1\text{\,}\mathrm{m}\mathrm{m}$. Two effects contribute to this fact. The
value of the opacity decreases and the opacity cliff moves to larger particle
sizes. This leads the disks to become optically thinner in comparison to the
$850\text{\,}\mu\mathrm{m}$ case. Secondly, the intensity at
$3.1\text{\,}\mathrm{m}\mathrm{m}$ is less according to the Plank’s spectrum.
Therefore the luminosity will be smaller.

In terms of the SLR, the slope for the smooth disks does not change since
these disks are never optically thick, therefore all disks simply move towards
lower luminosity and smaller radii. On the other hand, the sub-structured
disks that cover the SLR do not change in terms of slope but the large and
faint disks (right part of the heat-map in Figure 6) show a larger spread in
luminosity compared to the smaller wavelength. Disks that are very optically
thick and moderately optically thick have the same luminosity at
$\rm{\lambda=850\mu m}$, but at $3.1\text{\,}\mathrm{m}\mathrm{m}$ because of
the decrease in opacity the former category is still optically thick while the
latter no longer is, leading to a decrease in luminosity. Therefore the SLR
can become flatter if we can take into account these disks that do not belong
in the SLR.

With our models, the flatter relation from Tazzari et al. (2021), could be
explained by sub-structured and large smooth disks. The heat-map in Appendix D
confirms this. In this figure we plot the simulations at
$3.1\text{\,}\mathrm{m}\mathrm{m}$, using the R10-0 opacity model and we over-
plot the SLR from Tazzari et al. (2021) $\rm{L_{mm}\propto r_{eff}^{1.2}}$.
Sub-structured disks can very well explain this relation since it is similar
to the scaling relation we calculated for disks with strong substructures in
subsubsection 3.2.2. Small and smooth disks on the other hand cannot enter the
relation because they are too faint since the particles cannot grow to a size
where the opacity cliff is at $3.1\text{\,}\mathrm{m}\mathrm{m}$.

We have to mention the possibility that the flatter relation can be due to
observational bias towards large disks which tend to be sub-structured. If
small and faint disks are included in the sample, the observed SLR could be
steeper and closer to the SLR from Andrews et al. (2018b). Future
observational surveys should further investigate this possibility.

###  4.8 Limitations

It is important to keep in mind the limitations of this paper. The time-span
used for the simulations displayed in this paper and the figures is from
$300\text{\,}\mathrm{k}\mathrm{y}\mathrm{r}$ to
$3\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$. This does not exclude the
possibility that some disks with high disk mass might evolve a lot on the SLR
diagram for
$10\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$-$20\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$.
Disk dissipation has not been modeled in this paper and it is considered as
future work. In Figure 11 we show the kernel density estimate (kde) of the
global dust-to-gas ratio222The dust-to-gas ratio in the disk changes with
radius and time and this quantity is simply $\rm{M_{dust}/M_{gas}}$. for three
different snapshots between
$\mathrm{$300\text{\,}\mathrm{k}\mathrm{y}\mathrm{r}$-$1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$-$3\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$}$
and for three different cases. Smooth disks lose dust relatively quickly due
to radial drift, while disks with planets retain a much higher dust-to-gas
ratio because of the the strong trap. In the second panel there are cases
where the dust-to-gas ratio increases over the initial $\mathrm{0.01}$. These
are sub-structured disks with intermediate $\mathrm{\alpha}$-values, high
fragmentation velocity
($\mathrm{>$1000\text{\,}\mathrm{c}\mathrm{m}\mathrm{/}\mathrm{s}$}$) and
small size ($\mathrm{<$60\text{\,}\mathrm{a}\mathrm{u}$}$). The gas is removed
quicker than the dust, leading to a larger dust-to-gas ratio. Large
$\mathrm{\alpha}$ would lead to less trapping and the dust would drift as
usual while low $\mathrm{\alpha}$ would mean that the disk does not evolve
significantly.

As mentioned in section 2, the stellar luminosity is not evolved in the
simulation. If that were the case the disk luminosity would approximately
scale linearly with the stellar luminosity and further modulated by resulting
changes in the dust evolution. We therefore expect a general shift of the
disks towards lower luminosities but with the trends that have explored in
section 3 remaining the same. Since most of the simulations need to be
brighter to remain in the SLR, a change in the luminosity favors higher
$\rm{\alpha}-$values and smaller fragmentation velocities than the ones shown
before. An example of a heatmap in shown in Appendix D.

In our models the planets are already included at the beginning of the
simulations and they open a gap in the initially smooth surface density
profile relatively quickly. Realistically, the time-scales in the outer part
of the disk are much longer and the time-scale for planet formation changes
with the distance to the star (Johansen & Lambrechts, 2017). Therefore we
would expect that the inner planet should form first and the outer later as it
has been suggested by e.g. Pinilla et al. (2015). Since both planets start at
the same time, the inner one might trap more of the total disk mass and the
outer bump might be less bright than in our models in reality. The latter will
be included in a future work by including the outer planet later in the
simulation.

Figure 11: Evolution of the global disk dust-to-gas ratio of all matching
simulations with the Ricci et al. (2010) R10-0 opacity model, for three
different cases and three different snapshots, from
$\mathrm{$300\text{\,}\mathrm{k}\mathrm{y}\mathrm{r}$-$3\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$}$.
From left to right: Smooth case, a planet with planet/star mass ratio at
$\mathrm{1/3r_{\mathrm{c}}}$ and a planet with the same ratio at
$\mathrm{2/3r_{\mathrm{c}}}$. Different limits are used on the x-axis to
highlight the evolution of the dust-to-gas ratio. Initial dust-to-gas ratio is
$\mathrm{0.01}$. For the smooth case the dust-to-gas ratio is decreasing by
three orders of magnitude up to $3\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$.
When a planet is included the disk dust mass is retained and leads to a much
higher dust-to-gas ratio. In the case where the planet is the inner part of
the disk (middle column), there are cases at
$3\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ where the ratio is larger than
$\mathrm{0.01}$. The gas mass moves faster than the dust mass in this case.

##  5 Conclusions

In this paper we have performed a large population study of 1D models of gas
and dust evolution in protoplanetary disks, to study how the effective radius
and disk continuum emission evolves with time. We varied a range of initial
parameters and we included both smooth disks and disks that contain planets.
We compared our results with the observed trend between continuum sizes and
luminosities from Andrews et al. (2018b) and we managed to constrain the
initial conditions. Our findings are as follows.

  1. 1.

Disks with strong traps (i.e. massive planets) follow a different SLR than
smooth disks. Smooth disks follow the Andrews et al. (2018b) relation
$\mathrm{L_{mm}\propto r_{eff}^{2}}$ as shown by Rosotti et al. (2019a), while
disks with massive planets $\mathrm{L_{mm}\propto r_{eff}^{5/4}}$. This could
mean that not all disks in the Tripathi et al. (2017) and Andrews et al.
(2018b) joint sample have sub-structure as significant as e.g. HD163296. We
explained this result with a simple analytical derivation and we found that if
the gas width scales as we measured it from FARGO-3D and if the dust width
scales as we expect it from trapping and fragmentation, then theoretically the
luminosity scales as $\mathrm{L_{mm}\propto r_{eff}^{5/4}}$.

  2. 2.

If disks are following the SLR or not depends heavily on the opacity model.
When the DSHARP (Birnstiel et al., 2018) opacity is used, disks are not so
luminous in the first $300\text{\,}\mathrm{k}\mathrm{y}\mathrm{r}$ and the
majority of them end up below the SLR. Especially for smooth disks, the DSHARP
opacities produce a much lower number of simulations on the SLR compared to
models using the Ricci et al. (2010) R10-0 opacities ($\rm{0.8\%}$ with DSHARP
and $\rm{29.6\%}$ with R10-0). Therefore, with this opacity model, only disks
with sub-structures can populate the SLR. On the other hand, R10-0 opacities
can reproduce both disks with and without sub-structures since the absolute
value of the opacity at $850\text{\,}\mu\mathrm{m}$ is $\rm{\sim 8.5}$ times
larger than DSHARP, for particle sizes around
$\mathrm{\sim$0.1\text{\,}\mathrm{m}\mathrm{m}$}$ (position of the opacity
cliff) and the disks become luminous enough to enter the relation.

  3. 3.

The SLR is more widely populated when sub-structures are included in contrast
to a tight correlation for smooth disks. Sub-structured disks cover mostly the
upper right part (large and bright disks) of the SL diagram, while the lower
left (small and faint) is covered by smooth disks. This is an indication that
the SLR can be explained if there is a mixture of both smooth and strong sub-
structured disks.

  4. 4.

The grain porosity can drastically affect the evolution track of the disk.
Throughout our models, relatively compact grains ($\mathrm{<50\%}$ porosity)
are preferred for simulations that follow the SLR. If we use slightly porous
grains ($\mathrm{10\%}$) by altering the DSHARP opacity, the effect is
insignificant, as the shape of the opacity cliff is roughly the same. On the
contrary, for semi ($\mathrm{50\%}$) and porous grains ($\mathrm{90\%}$) the
opacity cliff flattens out, leading to disks with low luminosity. Only compact
grains can explain the SLR for smooth disks, while any porosity can explain it
when strong sub-structures are included.

  5. 5.

High initial disk mass gives a higher probability for a simulation to follow
the SLR. If this applies, the disk starts initially above the SLR (bright)
until it reaches a stable state at around
$\mathrm{\sim$300\text{\,}\mathrm{k}\mathrm{y}\mathrm{r}$}$. By this time it
will enter the relation and depending on the other initial conditions it will
either remain there and it will be considered as a matching simulation or
leave it.

  6. 6.

There is a preference towards low $\mathrm{\alpha}$-values (smaller than
$\mathrm{10^{-3}}$). This result is in line with other more direct methods of
determining $\mathrm{\alpha}$ (e.g. Flaherty et al. (2018)). There are
multiple reasons for this tendency. For $\mathrm{\alpha}\geq 2.5\cdot
10^{-3}$, disks tend to be more fragmentation dominated, the particle size
decreases and consequently are not trapped by the pressure bump (if any)
leading them outside the relation. Moreover the diffusivity increases and the
peak of the pressure bump smears out, leading to inefficient trapping. On the
other hand if $\mathrm{\alpha}$ is small, the ring that is formed is becoming
too narrow and disks tend to have smaller luminosity.

  7. 7.

The location of the planet as a function of the characteristic radius plays a
major role to the final outcome. If a planet is included in the inner part of
the disk ($\mathrm{1/3r_{c}}$), the disk has to be significantly larger in
order to retain the correct ratio of luminosity and effective radius to stay
in the SLR. Opposed to that, when an outer planet ($\mathrm{2/3r_{c}}$) is
included, the disk tends to be smaller in size. When two planets are included,
the location of the outermost one defines the size of the disk but a
combination of two defines the luminosity. These results are also affected by
the opacity model.

  8. 8.

We expect less extended evolution track when sub-structure is included. The
pressure bump halts the dust from drifting further in, constraining this way
the size of the disk and not allowing it to evolve further on the SLR.
Furthermore when two planets are included there is an indication that the
inner planet should form first otherwise there will not be a big enough
reservoir of material in order to form.

  9. 9.

We are not able to construct optically thick disks with high albedo
($\rm{0.9}$) that lie along the SLR with an evolutionary procedure as opposed
to Zhu et al. (2019). Smooth disks are not optically thick due to radial drift
while disks with sub-structure create only optically thick rings rather than a
uniform optically thick distribution.

  10. 10.

We chose different gap profiles based on Kanagawa et al. (2016) and compared
them again hydrodynamical simulations. We conclude that the depth of the gap
does not play an important role to the evolution of the disk on the SLR, as
long as the planet is big enough to stop the particles from drifting. The
width of the gap is the important parameter instead (see Appendix C, where we
compare the different profiles for different parameters.)

###### Acknowledgements.

T.B. acknowledges funding from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme under grant
agreement No 714769 and funding by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) - 325594231 under Ref no. FOR 2634/1. Furthermore,
this research was supported by the Excellence Cluster ORIGINS which is funded
by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany´s Excelence Strategy – EXC-2094-390783311. G.R. 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)

## References

  * ALMA Partnership et al. (2015) ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJ, 808, L3
  * Andrews (2020) Andrews, S. M. 2020, ARA&A, 58, 483
  * Andrews et al. (2018a) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018a, ApJ, 869, L41
  * Andrews et al. (2013) Andrews, S. M., Rosenfeld, K. A., Kraus, A. L., & Wilner, D. J. 2013, ApJ, 771, 129
  * Andrews et al. (2018b) Andrews, S. M., Terrell, M., Tripathi, A., et al. 2018b, ApJ, 865, 157
  * Andrews et al. (2010) Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2010, ApJ, 723, 1241
  * Andrews et al. (2012) Andrews, S. M., Wilner, D. J., Hughes, A. M., et al. 2012, ApJ, 744, 162
  * Andrews et al. (2016) Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJ, 820, L40
  * Ansdell et al. (2017) Ansdell, M., Williams, J. P., Manara, C. F., et al. 2017, AJ, 153, 240
  * Ansdell et al. (2016) Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46
  * Ballering & Eisner (2019) Ballering, N. P. & Eisner, J. A. 2019, AJ, 157, 144
  * Barenfeld et al. (2016) Barenfeld, S. A., Carpenter, J. M., Ricci, L., & Isella, A. 2016, ApJ, 827, 142
  * Beckwith et al. (1990) Beckwith, S. V. W., Sargent, A. I., Chini, R. S., & Guesten, R. 1990, AJ, 99, 924
  * Benítez Llambay & Masset (2015) Benítez Llambay, P. & Masset, F. 2015, FARGO3D: Hydrodynamics/magnetohydrodynamics code
  * Benz et al. (2014) Benz, W., Ida, S., Alibert, Y., Lin, D., & Mordasini, C. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 691
  * Birnstiel & Andrews (2014) Birnstiel, T. & Andrews, S. M. 2014, ApJ, 780, 153
  * Birnstiel et al. (2015) Birnstiel, T., Andrews, S. M., Pinilla, P., & Kama, M. 2015, ApJ, 813, L14
  * Birnstiel et al. (2010) Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79
  * Birnstiel et al. (2018) Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJ, 869, L45
  * Birnstiel et al. (2012) Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148
  * Blum (2018) Blum, J. 2018, Space Sci. Rev., 214, 52
  * Blum & Wurm (2008) Blum, J. & Wurm, G. 2008, ARA&A, 46, 21
  * Bohlin et al. (1978) Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132
  * Casassus et al. (2013) Casassus, S., van der Plas, G. M., Perez, S., et al. 2013, Nature, 493, 191
  * Cleeves et al. (2016) Cleeves, L. I., Öberg, K. I., Wilner, D. J., et al. 2016, ApJ, 832, 110
  * Dullemond et al. (2018) Dullemond, C. P., Birnstiel, T., Huang, J., et al. 2018, ApJ, 869, L46
  * Fedele et al. (2017) Fedele, D., Carney, M., Hogerheijde, M. R., et al. 2017, A&A, 600, A72
  * Fedele et al. (2018) Fedele, D., Tazzari, M., Booth, R., et al. 2018, A&A, 610, A24
  * Flaherty et al. (2018) Flaherty, K. M., Hughes, A. M., Teague, R., et al. 2018, ApJ, 856, 117
  * Gundlach & Blum (2015) Gundlach, B. & Blum, J. 2015, ApJ, 798, 34
  * Güttler et al. (2010) Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56
  * Hendler et al. (2020) Hendler, N., Pascucci, I., Pinilla, P., et al. 2020, ApJ, 895, 126
  * Hendler et al. (2017) Hendler, N. P., Mulders, G. D., Pascucci, I., et al. 2017, ApJ, 841, 116
  * Henning & Stognienko (1996) Henning, T. & Stognienko, R. 1996, A&A, 311, 291
  * Huang et al. (2018a) Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018a, ApJ, 869, L42
  * Huang et al. (2018b) Huang, J., Andrews, S. M., Pérez, L. M., et al. 2018b, ApJ, 869, L43
  * Isella et al. (2016) Isella, A., Guidi, G., Testi, L., et al. 2016, Phys. Rev. Lett., 117, 251101
  * Isella et al. (2012) Isella, A., Pérez, L. M., & Carpenter, J. M. 2012, ApJ, 747, 136
  * Johansen & Lambrechts (2017) Johansen, A. & Lambrechts, M. 2017, Annual Review of Earth and Planetary Sciences, 45, 359
  * Kanagawa et al. (2016) Kanagawa, K. D., Muto, T., Tanaka, H., et al. 2016, PASJ, 68, 43
  * Kataoka et al. (2014) Kataoka, A., Okuzumi, S., Tanaka, H., & Nomura, H. 2014, A&A, 568, A42
  * Kenyon et al. (1996) Kenyon, S. J., Yi, I., & Hartmann, L. 1996, ApJ, 462, 439
  * Kratter & Lodato (2016) Kratter, K. & Lodato, G. 2016, ARA&A, 54, 271
  * Krijt et al. (2016) Krijt, S., Ormel, C. W., Dominik, C., & Tielens, A. G. G. M. 2016, A&A, 586, A20
  * Lynden-Bell & Pringle (1974) Lynden-Bell, D. & Pringle, J. E. 1974, MNRAS, 168, 603
  * Min et al. (2016) Min, M., Rab, C., Woitke, P., Dominik, C., & Ménard, F. 2016, A&A, 585, A13
  * Miyake & Nakagawa (1993) Miyake, K. & Nakagawa, Y. 1993, Icarus, 106, 20
  * Nakagawa et al. (1986) Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375
  * Okuzumi et al. (2009) Okuzumi, S., Tanaka, H., & Sakagami, M.-a. 2009, ApJ, 707, 1247
  * Pascucci et al. (2016) Pascucci, I., Testi, L., Herczeg, G. J., et al. 2016, ApJ, 831, 125
  * Pérez et al. (2016) Pérez, L. M., Carpenter, J. M., Andrews, S. M., et al. 2016, Science, 353, 1519
  * Piétu et al. (2014) Piétu, V., Guilloteau, S., Di Folco, E., Dutrey, A., & Boehler, Y. 2014, A&A, 564, A95
  * Pinilla et al. (2012a) Pinilla, P., Benisty, M., & Birnstiel, T. 2012a, A&A, 545, A81
  * Pinilla et al. (2012b) Pinilla, P., Birnstiel, T., Ricci, L., et al. 2012b, A&A, 538, A114
  * Pinilla et al. (2015) Pinilla, P., Birnstiel, T., & Walsh, C. 2015, A&A, 580, A105
  * Pinilla et al. (2020) Pinilla, P., Pascucci, I., & Marino, S. 2020, A&A, 635, A105
  * Ricci et al. (2010) Ricci, L., Testi, L., Natta, A., et al. 2010, A&A, 512, A15
  * Rice et al. (2006) Rice, W. K. M., Armitage, P. J., Wood, K., & Lodato, G. 2006, MNRAS, 373, 1619
  * Rosotti et al. (2019a) Rosotti, G. P., Booth, R. A., Tazzari, M., et al. 2019a, MNRAS, 486, L63
  * Rosotti et al. (2019b) Rosotti, G. P., Tazzari, M., Booth, R. A., et al. 2019b, MNRAS, 486, 4829
  * Seizinger & Kley (2013) Seizinger, A. & Kley, W. 2013, A&A, 551, A65
  * Seizinger et al. (2013) Seizinger, A., Speith, R., & Kley, W. 2013, A&A, 559, A19
  * Shakura & Sunyaev (1973) Shakura, N. I. & Sunyaev, R. A. 1973, Astronomy and Astrophysics, 24, 337
  * Siess et al. (2000) Siess, L., Dufour, E., & Forestini, M. 2000, A&A, 358, 593
  * Stammler et al. (2019) Stammler, S. M., Drażkowska, J., Birnstiel, T., et al. 2019, ApJ, 884, L5
  * Takeuchi & Lin (2002) Takeuchi, T. & Lin, D. N. C. 2002, ApJ, 581, 1344
  * Tazaki et al. (2016) Tazaki, R., Tanaka, H., Okuzumi, S., Kataoka, A., & Nomura, H. 2016, ApJ, 823, 70
  * Tazzari et al. (2021) Tazzari, M., Clarke, C. J., Testi, L., et al. 2021, MNRAS, 506, 2804
  * Trapman et al. (2020) Trapman, L., Ansdell, M., Hogerheijde, M. R., et al. 2020, A&A, 638, A38
  * Tripathi et al. (2017) Tripathi, A., Andrews, S. M., Birnstiel, T., & Wilner, D. J. 2017, The Astrophysical Journal, 845, 44
  * Tsukagoshi et al. (2016) Tsukagoshi, T., Nomura, H., Muto, T., et al. 2016, ApJ, 829, L35
  * van der Marel & Mulders (2021) van der Marel, N. & Mulders, G. D. 2021, AJ, 162, 28
  * van der Marel et al. (2013) van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2013, Science, 340, 1199
  * Wada et al. (2008) Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2008, ApJ, 677, 1296
  * Wada et al. (2009) Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2009, ApJ, 702, 1490
  * Weidenschilling (1977) Weidenschilling, S. J. 1977, MNRAS, 180, 57
  * Whipple (1972) Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius, 211
  * Youdin & Lithwick (2007) Youdin, A. N. & Lithwick, Y. 2007, Icarus, 192, 588
  * Zhang et al. (2018) Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJ, 869, L47
  * Zhu et al. (2012) Zhu, Z., Nelson, R. P., Dong, R., Espaillat, C., & Hartmann, L. 2012, ApJ, 755, 6
  * Zhu et al. (2019) Zhu, Z., Zhang, S., Jiang, Y.-F., et al. 2019, ApJ, 877, L18
  * Zubko et al. (1996) Zubko, V. G., Mennella, V., Colangeli, L., & Bussoletti, E. 1996, MNRAS, 282, 1321

##  Appendix A Derivation of gap width and SLR

In this section we show the complete derivation of the width vs radius $r$,
planet/star mass ratio $q$ and scale scale height $h$ relation as explained in
subsubsection 3.2.1 and the derivation of the SLR of disks with strong sub-
structures as shown in subsubsection 3.2.2.

###  A.1 Derivation of the gap width

In addition to our two-pop-py models, we performed 24 hydrodynamical
simulations with FARGO-3D (Benítez Llambay & Masset 2015) for different
planet/star mass ratios, planet locations and $\mathrm{\alpha}$-values (see
Table 3 in Appendix C). We used these simulations to calculate the width of
the outer pressure bump caused by the planet in the gas. The surface density
maximum is locally well fitted by a Gaussian which allows us to measure the
width (i.e. the standard deviation) using the curvature at the maximum.

$y(x)=\frac{1}{\sqrt{2\pi}\sigma_{g}}e^{-\frac{(x-x_{0})^{2}}{2\sigma^{2}}}\,$
(20)

where $\mathrm{x_{0}}$ is the location of the surface density maximum. We took
the logarithmic derivative of the Gaussian,

$\frac{\partial lny}{\partial
lnx}=\frac{x}{y}\frac{1}{\sqrt{2\pi}\sigma_{g}}\left(-\frac{x-x_{0}}{\sigma^{2}}\right)e^{-\frac{(x-x_{0})^{2}}{2\sigma_{g}^{2}}}\,$
(21)

which gives us

$\frac{\partial lny}{\partial lnx}=-x\frac{x-x_{0}}{\sigma_{g}^{2}}$ (22)

If we are close to the pressure maximum we can approximate that
$\mathrm{x\xrightarrow{}x_{0}+\delta x}$, and the equation yields

$\frac{dlny}{dlnx}=-\frac{x_{0}\delta x}{\sigma_{g}^{2}}$ (23)

If we take the normal derivative of the above expression we end up with the
following expression

$\frac{\partial}{\partial\delta x}\frac{\partial lny}{\partial
lnx}=-\frac{x_{0}}{\sigma_{g}^{2}}$ (24)

Solving for $\mathrm{\sigma_{g}}$, we measure the width of the maximum by
calculating this derivative at the peak location for the azimuthally averaged
surface density profile derived from the FARGO3D runs.

$\sigma=\sqrt{-\frac{x_{0}}{\frac{\partial}{\partial\delta x}\frac{\partial
lny}{\partial lnx}}}$ (25)

We measured all the the widths of our hydrodynamical simulations and we
searched how do they scale with radius, scale height and the
$\mathrm{\alpha}$-parameter. We fit the width of the pressure bump as multiple
power-law as:

$\sigma_{g}=C\cdot h^{p}\cdot q^{k}\cdot q^{l}\,$ (26)

where $\mathrm{C}$ is a constant, $h$ is the scale height, $q$ the mass ratio
and $\mathrm{\alpha}$ the turbulence parameter. We found that the width in the
measured range scales approximately as:

$\sigma_{g}\propto h^{0.81}\cdot q^{0.14}\cdot\alpha^{0.05}\,$ (27)

###  A.2 Derivation of SLR of disks with companions

To understand the overall shape of the heat-map for the case with strong sub-
structures (i.e. the red lines in Figure 6 and Figure 7), we present here the
complete derivation of the theoretical estimate from subsubsection 3.2.2. This
estimate is based on the empirically determined widths of the previous
subsection A.1 as explained in subsubsection 3.2.1 and the position of the
pressure maximum formed outside the position of the gap-opening planets..

Assuming that all the luminosity of a disk comes from rings that are
approximately optically thick, we can approximate:

$L\simeq A\cdot B_{\nu}\,$ (28)

where $\mathrm{A}$ is the area and $\mathrm{B_{\nu}}$ is the Planck function.
If we assume the Rayleigh Jeans approximation, then we can substitute the
Planck function with the temperature $\mathrm{T}$ and the equation becomes

$L\propto A\cdot T$ (29)

We make the assumption that the area of the pressure bump scales as
$\mathrm{A}\propto r\cdot\sigma_{d}$, where $\mathrm{r}$ is the radius and
$\mathrm{\sigma_{d}}$ is the width of the pressure bump in the dust and there
is linear scaling of $\mathrm{\sigma_{d}}$ with $\mathrm{h}$. The temperature
is $\mathrm{T\propto r^{-1/2}}$ and $\mathrm{h\propto r^{5/4}}$. Combining the
above relations we obtain the theoretical expression for the luminosity,

$L_{mm}\propto r^{7/4}$ (30)

In the previous section we determined how the width of the pressure maximum
depends on the scale height $\mathrm{\sigma_{g}\propto h^{0.81}}$, so we need
to find a relation between the $\sigma_{g}$ and $\sigma_{d}$. In the pressure
maximum, we assume that the disk is optically thick and that it is
fragmentation limited since there is no radial drift. If the latter applies,
then,

$\sigma_{d}\propto\sigma_{g}\cdot\sqrt{\frac{\alpha}{St}}\,$ (31)

as in (Dullemond et al. 2018), where $\mathrm{St}$ is the Stokes number.
Moreover,

$St\propto\frac{1}{\alpha\cdot c_{s}^{2}}\,$ (32)

where $\mathrm{c_{s}}$ is the sound speed. Combining these two, we get that:

$h_{d}\propto h_{g}\cdot\alpha\cdot c_{s}$ (33)

The sound-speed is,

$c_{s}\propto r^{-1/4}$ (34)

Using all the above equations we conclude that,

$\sigma_{d}\propto\sigma_{g}\cdot r^{-1/4}$ (35)

or if we express the radius as a function of the scale height,

$\sigma_{d}\propto\sigma_{g}\cdot h^{-1/5}$ (36)

The scale height is $\mathrm{h=c_{s}/\Omega}$, where
$\mathrm{\Omega=\sqrt{GM_{\star}/r^{3}}}$ the Keplerian angular velocity,
leading to a dependency $\mathrm{h\propto r^{5/4}}$. Therefore we can use
Equation 27 and express the calculated $\mathrm{\sigma_{g}}$ as a function of
radius.

$\sigma_{g}\propto r^{1.01}$ (37)

Finally we can calculate the dust width $\mathrm{\sigma_{d}}$ as a function of
both the scale height and radius. Substituting into Equation 36, we get:

$\sigma_{d}\propto h^{0.61}$ (38)

and into Equation 35

$\sigma_{d}\propto r^{0.75}$ (39)

If we calculate again how the luminosity scales and replace the scaling of the
dust width $\mathrm{\sigma_{d}}\propto r^{0.75}$, we find that:

$L_{mm}\propto r^{5/4}$ (40)

which is the value that we did fit with the red line in Figure 6.

##  Appendix B Corner plots

To search for further correlations between the parameters we visualized the
results similar to a corner plot, however for a grid of parameters instead of
the usual sample density. The corner plots are projections of a multi-
dimensional parameter space to several two dimensional planes. In these
visualizations, every box of the plot is a 2D histogram of a parameter pair
while the color varies according to the number of simulations in every cell.
The color bar is given on the upper part of the plot. These visualizations
contain only the simulations that match the SLR as the histograms in Figure 5.

In Figure 12 we show the corner plot for the smooth case. The corner plots
overall confirm the results discussed in the previous sections. Starting from
the upper left part of the plot, we show the correlation between the
turbulence parameter $\mathrm{\alpha}$ and the fragmentation velocity. The
results show that all these plots are heavily dependent on whether the disks
are drift or fragmentation limited. The white dots confirm that a scaling
$\mathrm{\propto\sqrt{\alpha}}$ well explains the correlation: simulations
above this relation are strongly fragmentation limited and therefore do not
follow the SLR (see Rosotti et al. 2019b). The very low left cell
(corresponding to $\rm{\alpha}=10^{-4}$ and $\rm{v_{frag}}=$
$200\text{\,}\mathrm{c}\mathrm{m}\mathrm{/}\mathrm{s}$) is empty because it is
fragmentation limited and increasing the $\mathrm{\alpha}$-value only makes
the disks more fragmentation dominated. As the fragmentation velocity
increases then the allowed $\mathrm{\alpha}$-value increases too.

In Figure 13 and Figure 14, we show the corner plots of the cases where a
planet with a planet/star mass ratio $\mathrm{q=10^{-3}}$ is included at 1/3
and 2/3 of the $r_{\mathrm{c}}$ respectively. On both cases we obtain similar
results. On the second column a mild covariance can be seen between the
turbulence parameter $\mathrm{\alpha}$ and the disk mass. Most of the
simulations are towards low $\mathrm{\alpha}-$values and high disk masses
confirming once again the results from the 1D histograms but showing that
these two are the most dominant parameters.

In the bottom panel a strong correlation between $\mathrm{\alpha}$ and
$r_{\mathrm{c}}$ (for a limited $\mathrm{\alpha}$-range) is shown. With
vertical white dots, we mark the $\mathrm{\alpha=5\cdot 10^{-3}}$ value, where
there is a gradual boundary from matching to discrepant simulations. This
relation is visible on the second and third column of the Figure 6. In
summary, the very left clump of simulations above the correlation is the case
where $\mathrm{r_{c}}=$ $10\text{\,}\mathrm{a}\mathrm{u}$. This clump is
mostly above the correlation (unless $\mathrm{\alpha=10^{-4}}$) and it
explains the empty row for $\mathrm{r_{c}}=$ $10\text{\,}\mathrm{a}\mathrm{u}$
in the corner plots of Figure 13 and Figure 14.

With low $\mathrm{\alpha}$-values there is more trapping in the pressure bump
and the emission is optically thick, hence brighter disks. Since the disks
should not be too bright (because that would locate them above the SLR) or not
bright enough (this would locate them below the SLR) this yields only a range
of permitted $\mathrm{\alpha}$-values. This range depends on $r_{\mathrm{c}}$.
Moreover, when $\mathrm{\alpha}$ increases, the emission is not optically
thick anymore because the dust gets through the trap. More specifically, for
$\mathrm{\alpha}\geq 5\cdot 10^{-3}$ (as seen in Figure 3), the bump cannot
efficiently trap anymore and this is why this sharp transition from matching
to discrepant simulations occurs.

Therefore, in the case where a planet is included the strongest correlation is
between $\mathrm{\alpha}$-value and the characteristic radius
$r_{\mathrm{c}}$. The area of the parameter space with the highest fraction of
matching simulations is where we have small disks (around
$\rm{\sim$50\text{\,}\mathrm{a}\mathrm{u}$}$) and low
$\mathrm{\alpha}$-values. This allows once again to constrain the disk
characteristic radius to a value that is not larger than
$100\text{\,}\mathrm{a}\mathrm{u}$.

Figure 12: Corner plot of the smooth case. From left to right: Fragmentation
velocity $\mathrm{v_{frag}}[cm/s]$, turbulence $\mathrm{\alpha}-$value, disk
mass $\mathrm{M_{d}}[M_{\star}]$, characteristic radius $\mathrm{r_{c}}[au]$.
From top to bottom: $\mathrm{\alpha}$, $\mathrm{M_{d}}[M_{\star}]$,
$\mathrm{r_{c}}[au]$, $\mathrm{M_{\star}}[M_{\odot}]$. On the upper left plot,
the white dots show the curve where the distinction between the fragmentation
and the drift limited regime is. On the middle plot
($\mathrm{r_{\mathrm{c}}-\alpha}$) the white dots show the $\mathrm{\alpha}$
value where the pressure bump cannot hold the dust efficiently anymore. On the
panel $r_{\mathrm{c}}$\- $\mathrm{M_{d}}$, marked with white dots is the mean
value of the characteristic radius for every disk mass. Figure 13: Corner plot
with a planet with of a planet/star mass ratio $\mathrm{10^{-3}}$ at a
location of the 1/3$r_{\mathrm{c}}$. From left to right: Fragmentation
velocity $\mathrm{v_{frag}}[cm/s]$, turbulence $\mathrm{\alpha}-$value, disk
mass $\mathrm{M_{d}}[M_{\star}]$, characteristic radius $\mathrm{r_{c}}[au]$.
From top to bottom: $\mathrm{\alpha}$, $\mathrm{M_{d}}[M_{\star}]$,
$\mathrm{r_{c}}[au]$, $\mathrm{M_{\star}}[M_{\odot}]$. On the upper left plot,
the white dots show the curve where the distinction between the fragmentation
and the drift limited regime is. On the middle plot
($\mathrm{r_{\mathrm{c}}-\alpha}$) the white dots show the $\mathrm{\alpha}$
value where the pressure bump cannot hold the dust efficiently anymore. On the
panel $r_{\mathrm{c}}$\- $\mathrm{M_{d}}$, marked with white dots is the mean
value of the characteristic radius for every disk mass. Figure 14: Corner plot
where we have a Jupiter mass planet at a location of the 2/3$r_{\mathrm{c}}$.
From left to right: Fragmentation velocity $\mathrm{v_{frag}}[cm/s]$,
turbulence $\mathrm{\alpha}-$value, disk mass $\mathrm{M_{d}}[M_{\star}]$,
characteristic radius $\mathrm{r_{c}}[au]$. From top to bottom:
$\mathrm{\alpha}$, $\mathrm{M_{d}}[M_{\star}]$, $\mathrm{r_{c}}[au]$,
$\mathrm{M_{\star}}[M_{\odot}]$. On the upper left plot, the white dots show
the curve where the distinction between the fragmentation and the drift
limited regime is. On the middle plot ($\mathrm{r_{\mathrm{c}}-\alpha}$) the
white dots show the $\mathrm{\alpha}$ value where the pressure bump cannot
hold the dust efficiently anymore. On the panel $r_{\mathrm{c}}$\-
$\mathrm{M_{d}}$, marked with white dots is the mean value of the
characteristic radius for every disk mass.

##  Appendix C Gap profiles

As discussed in section 2, we compared the planet profile from Kanagawa et al.
(2016), against hydrodynamical simulations using FARGO (Benítez Llambay &
Masset 2015). We varied the turbulence parameter $\mathrm{\alpha}$, the
planet/star mass ratio $\mathrm{q}$, and the position of the planetary gap.
After obtaining the planet profile from the hydrodynamical simulation, we
evolved the disk with two-pop-py (Birnstiel et al. 2012) without gas
evolution, to avoid diffusion of the planetary gap. We obtained the disk
continuum luminosity and the effective radius, and we plotted the evolution
tracks on the SLR. We evolved the simulations for
$330\text{\,}\mathrm{k}\mathrm{y}\mathrm{r}$. The parameters that we used for
this test are summarized in the Table 3.

Table 3: Parameters for the gap comparison Parameter | Description | Values
---|---|---
$\alpha$ | viscosity parameter |  $10^{-4}$, $10^{-3}$, $10^{-3}$
rp [$au$] | planet position |  $30$, $60$, $100$
q | planet/star mass ratio |  $3\cdot 10^{-4}$, $10^{-3}$, $3\cdot 10^{-3}$

In Figures 15, 16, and 17, we show the comparison between the gap that is
obtained from the hydrodynamical simulation and the profile from Kanagawa et
al. (2016). In Figure 15 we compare the gap profile for a small planet with
planet/star mass ratio of $\mathrm{q=3\cdot 10^{-4}}$, for the three different
values of the $\mathrm{\alpha}$-parameter. From the low
$\mathrm{\alpha}$-value (first row, first column, $\mathrm{\alpha}=10^{-4}$),
we observe that the Kanagawa et al. (2016) profile fits reasonably well the
width of the gap but there is an offset on the depth. As we look on their
evolution tracks we do not observe any major difference. For
$\mathrm{\alpha}=10^{-3}$ (second row, first column) we have a good fit for
the depth and the width in Figure 16 and 17. On the evolution track, there is
a small fluctuation, but in the end of the simulation they converge to the
same point. For the high $\mathrm{\alpha}$-value (third row, first column),
the profile fits well the width and the depth of the gap. As a result, the
evolution tracks of the simulations (third row, second column) are almost
identical.

Similarly, for the other cases we can observe that even if the profile does
not fit perfectly the bottom of the gap the difference on the evolution tracks
is minimal. The track can be affected is only if the width of the gap is not
fitted properly. Therefore, we concluded that our chosen profile will not
affect the main point of the paper which is based on the observable
quantities.

Figure 15: Comparison between hydrodynamical and our dust evolution simulation
for different $\mathrm{\alpha}$-values and for a planet/star mass ratio
$\mathrm{q=3\cdot 10^{-4}}$ at $30\text{\,}\mathrm{a}\mathrm{u}$. On the left
column are the different gap profiles for different $\mathrm{\alpha}$-values
and on the right column the corresponding evolution tracks. We observe that
even though the profile from Kanagawa et al. (2016) is not overlapping with
the one obtained from the hydrodynamical simulation, the tracks on the right
column produce similar results. Figure 16: Comparison between hydrodynamical
and our dust evolution simulation for different $\mathrm{\alpha}$-values and
for a planet/star mass ratio $\mathrm{q=10^{-3}}$ at
$30\text{\,}\mathrm{a}\mathrm{u}$. On the left column are the different gap
profiles for different $\mathrm{\alpha}$-values and on the right column the
corresponding evolution tracks. We observe that even though the profile from
Kanagawa et al. (2016) is not overlapping with the one obtained from the
hydrodynamical simulation, the tracks on the right column produce similar
results. Figure 17: Comparison between hydrodynamical and our dust evolution
simulation for different $\mathrm{\alpha}$-values and for a planet/star mass
ratio $\mathrm{q=3\cdot 10^{-3}}$ at $100\text{\,}\mathrm{a}\mathrm{u}$. On
the left column are the different gap profiles for different
$\mathrm{\alpha}$-values and on the right column the corresponding evolution
tracks. We observe that even though the profile from Kanagawa et al. (2016) is