© ESO 2022

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 has also enabled several large intermediate resolution (on the order of 100 mas) 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 intermediate 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 (L_mm) at (sub-)millimeter 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 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 measure the disk dust mass (M_dust) have been used to correlate this property with the mass of the host star (M_*), 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, a steeper than linear relationship has been observed between L_mm and M_* (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 outward, known as viscous spreading, suggesting the disk size should increase with time. In principle a measurement of disk size as a 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 and dissociation, among others) 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 (r_eff) (Tripathi et al. 2017).

However, the dust component does not behave in the same way as 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 toward the star. As an observational implication, the dust emission is less extended than the gas emission (e.g., Andrews et al. 2012, 2016; Isella et al. 2012; Cleeves et al. 2016) 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 subarcsecond resolution survey of 50 nearby protoplanetary disks, conducted with the Submillimeter Array (SMA) by Tripathi et al. (2017), showed a strong size-luminosity relation (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 the ALMA data from the Lupus disk sample (105 disks in total), confirmed the scaling relations between the r_eff and the L_mm, M_*. 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 has 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 substructures, for example HL Tauri (ALMA Partnership 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 for the origin of these ring-like substructures have been explored 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 (≥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 migrate toward the local pressure bump, due to radial drift (Weidenschilling 1977; Nakagawa et al. 1986), and consequently lead 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 between these features the material is approximated as optically thin (Dullemond et al. 2018). However, these rings 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 substructures. In Sect. 2, we discuss the methods we used to carry out our computational models. The results of this analysis are presented in Sect. 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 Sect. 4, we discuss the theoretical and observational implications of our results. We draw our conclusions in Sect. 5.

2 Methods

We carry out 1D 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 follows the viscous evolution equation. For the disk evolution, we use the turbulent effective viscosity as in Shakura & Sunyaev (1973), $$v={\alpha }_{\text{gas}}\frac{c_{\text{s}}^2}{{\text{Ω}}_{\text{K}}},$$(1)

and the dust diffusion coefficient as $$\text{D}={\alpha }_{\text{dust}}\frac{c_{\text{s}}^2}{{\text{Ω}}_{\text{K}}},$$(2)

with α_gas being the turbulence parameter, c_s the sound speed, and Ω_K the Keplerian frequency. The above equation lacks the term $\frac{1}{1+{\text{St}}^2}$, where St is the Stokes number, but we can ignore it since the Stokes number is always <1 in our simulations (Youdin & Lithwick 2007). We differentiate the a parameter in two different values. One for the gas α_gas and one for the dust α_dust, since we later mimic planetary gaps by locally varying the viscosity (see Sect. 2.3). In the smooth case α_dust = α_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 inward. 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), $${\displaystyle \sum {}_{\text{g}}\left(r\right)=}{\displaystyle \sum {}_0}{\left(\frac{r}{r_c}\right)}^{-\gamma }\text{EXP}\left[-{\left(\frac{r}{r_c}\right)}^{2-\gamma }\right],$$(3)

where ${\sum }_0=(2-\gamma )M_{\text{d}}/2\pi r_{\text{c}}^2$ is the normalization parameter, which is set in every simulation by the disk mass M_d. The other parameters are γ = 1, which is the viscosity exponent and it is not varied throughout our models and rc which is the characteristic radius of the disk (see Table 1). When r ≪ r_c, then Σ_g is a power law and when r ≥ r_c, Σ_g is dominated by the exponential factor.

The initial dust distribution follows the gas distribution with a constant dust-to-gas ratio of Σ_d/Σ_g = 0.01. The initial grain size (= monomer grain size) is a_min = 0.1 µm. This monomer size stays constant in time and space, while the representative size for the large grains increases with time as the particles grow. The particle bulk density is ρ_s = 1.7 g cm⁻³ for the standard opacity model from Ricci et al. (2010) and ρ_s = 1.675 g cm⁻³ for the DSHARP (Birnstiel et al. 2018) opacity, but decreases for different values of porosity (see Sect. 2.4). We evolve the disks up to 10 Myr to study the long-term evolution, but in the following analysis we only show results from 300 kyr to 3 Myr (see Sect. 2.5).

Our 1D radial grid ranges from 0.05 to 2000 au and the grid cells are spaced logarithmically. We use adaptive temperature, which depends on the luminosity of every star. Since the stellar mass changes in our grid, the stellar luminosity changes as well. We follow the temperature profile, $$T={\left(\varphi \frac{L_{★}}{4\pi {\sigma }_{\text{SB}}{\text{r}}^2}+{\left(10\text{K}\right)}^4\right)}^{1/4},$$(4)

as in Kenyon et al. (1996). In this equation, L_* is the stellar luminosity, ϕ = 0.05 is the flaring angle, σ_SB is the Stefan-Boltzmann constant, and r is the radius. The term 10⁴ is a lower limit so that we do not allow the disk temperature to drop below 10 K at the outer parts of the disk. We use the evolutionary tracks of Siess et al. (2000) to get the luminosity of a 1 Myr old star of the given mass. The stellar luminosity and effective temperature is not evolved in our simulations. However, the luminosity of a 1 M_⊙ star would decrease from ~2.4 L_⊙ at 1 Myr to ~ 1 L_⊙ at 3 Myr. In Sect. 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 (α_gas), disk mass (M_d), stellar mass (M_*), characteristic radius (r_c), and fragmentation velocity (υ_frag). For every parameter, we pick the 10 values specified in Table 2, taking all the possible combinations between them leading to a total of 100 000 simulations.

2.3 Planets

A large planet embedded in a disk produces a co-orbital gap in the gas density. To mimic gap opening by planets in our simulations, we altered the α_gas turbulence parameter. Since in steady state α_gas · Σ_g is constant, the α parameter and the surface density Σ_g are inversely proportional quantities, so a bump in the α_gas profile leads to a gap in the surface density profile. The reason for the change in α_gas and not in Σ_g is that the surface density evolves according to Eq. (3). By inserting the bump in the α_gas, the Σ_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⁻³ represents a Jupitermass 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 does not necessarily represent 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 them 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. 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 example in 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 r_p, as a function of the characteristic radius r_c (see Table 1). We locate them either at 2/3 or at 1/3 of r_c. In our simulations we used zero, one, or two planets in these positions. We refer the reader to Sect. 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_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 defining an effective radius (r_eff) as the radius that encloses a fixed fraction of the total flux, f_v (r_eff) = xF_v. We choose x = 68% of the total disk flux as a suitable intermediate value to define r_eff as it is comparable to a standard deviation in the approximation of a Gaussian profile.

We calculate the mean intensity J_v profile by using the scattering solution from Miyake & Nakagawa (1993) of the radiative transfer equation $$\frac{J_v\left({\tau }_v\right)}{B_v\left(T\left(r\right)\right)}=1-b\left(e^{-\sqrt{3{ϵ}_v^{\text{eff}}}\left(\frac{1}{2}\text{Δ}\tau -{\tau }_v\right)}+e^{-\sqrt{3{ϵ}_v^{\text{eff}}}\left(\frac{1}{2}\text{Δ}\tau +{\tau }_v\right)}\right),$$(5)

with $$b={\left[\left(1-\sqrt{{ϵ}_v^{\text{eff}}}\right)e^{-\sqrt{3{ϵ}_v^{\text{eff}}}\text{Δ}\tau }+1\sqrt{{ϵ}_v^{\text{eff}}}\right]}^{-1},$$(6)

where B_v is the Planck function and $${\tau }_v=\left({\kappa }_v^{\text{abs}}+{\kappa }_v^{\text{sca,eff}}\right){\displaystyle \sum {}_{\text{d}}}$$(7)

is the optical depth with ${\kappa }_{\nu }^{\text{abs}}$ the dust absorption opacity and ${\kappa }_{\nu }^{\text{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 }_v^{\text{sca,eff}}=\left(1-g_v\right){\kappa }_v^{\text{sca}},$$(8)

where g_v is the forward-scattering parameter, Δτ is $$\text{Δ}\tau ={\displaystyle \sum {}_{\text{d}}}{\kappa }_v^{\text{tot}}\text{Δ}z,$$(9)

while $${ϵ}_v^{\text{eff}}=\frac{{\kappa }_v^{\text{abs}}}{{\kappa }_v^{\text{abs}}+{\kappa }_v^{\text{sca,eff}}}$$(10)

is the effective absorption probability. To calculate the intensity $I_{\nu }^{\text{out}}$ we follow the modified Eddington–Barbier approximation as in Birnstiel et al. (2018), $$I_v^{\text{out}}\simeq \left(1-{\text{e}}^{-\text{Δ}\tau /\mu }\right)S_v\left(\left(\frac{1}{2}\text{Δ}\tau -{\tau }_v\right)/\mu =2/3\right),$$(11)

where µ = cos θ and $$S_v\left({\tau }_v\right)={ϵ}_v^{\text{eff}}{B}_v\left(T_{\text{d}}\right)+\left(1-{ϵ}_v^{\text{eff}}\right)J_v\left({\tau }_v\right)$$(12)

is the source function.

Two-pop-py evolves only the dust and gas surface densities and the maximum particle 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, n(a) ∝ a^−q, with an exponent q = 2.5, for a_min ≤ a ≤ a_max. This choice follows Birnstiel et al. (2012) where the size distribution is closer to q = 2.5 for disks that are in the drift limit, while for the fragmentation-limited ones, a choice of q = 3.5 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 q = 2.5 fits these disks better. Moreover, if a disk is fragmentation-limited then it is so mostly in the inner part, but the 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 substructures can be fragmentation-limited farther 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 10% silicates, 30% carbonaceous materials, and 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 µm). Afterward, 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 the Ricci et al. (2010) opacities (hereafter denoted R10-0), but for compact grains (i.e., without porosity) as in the model of Rosotti et al. (2019a). Furthermore, we used the DSHARP opacity model (Birnstiel et al. 2018) and we altered the grain porosity to 10% (little porous, DSHARP-10), 50% (semi-porous, DSHARP-50), and 90% (very porous, DSHARP-90). The particle bulk densities for the different porous grains are ρ_s = 1.508 g cm⁻³, ρ_s = 0.838 g cm⁻³, and ρ_s = 0.168 g cm⁻³, respectively. An important feature of the opacity models that we used is called the opacity cliff, which refers to the sharp drop in the opacity at 850 µm, at a maximum particle size around 0.1 mm, as defined in Rosotti et al. (2019a) (see Fig. 1). In all the figures that are shown in this paper the R10-0 opacity model from Ricci et al. (2010) is used, unless it is explicitly stated otherwise.

2.5 Matching simulations

The behavior of the simulations on the size-luminosity diagram (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) gives a relation between the disk size and the 340 GHz luminosity. The effective radius r_eff and the luminosity L_mm are correlated as $$\log r_{\text{eff}}=\left({2.10}_{-0.03}^{+0.06}\right)+\left({0.49}_{-0.03}^{+0.05}\right)\log L_{\text{mm}},$$(13)

with a Gaussian scatter perpendicular to that scaling with a standard deviation (1σ) of ${0.20}_{-0.01}^{+0.02}$ dex (where r_eff is in au and L_mm is in Jy at a distance of 140 pc).

In Fig. 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. They move from the top right (higher luminosity and radius) to the bottom left as the disk evolves. We plot the evolution track from 300 kyr to 3 Myr. The thought behind this is that at roughly 300 kyr our disks reach a quasi-steady state. Specifically, the drift speed in the drift limit is V_r = ϵVK, with ϵ being the dust-to-gas ratio and 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 Sect. 4, where we show the evolution of disk dust-to-gas ratio as a function of time for different cases. Longer evolutionary times than we are exploring here do not alter our results significantly, and are not included to simplify the discussion. They 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 divided them in categories. A simulation that lies within 1σ (blue shaded region in Fig. 2, top left, of the SLR) (Eq. (13)) at all times in our chosen time span is considered matching (see Fig. 2, green evolution track). On the other hand, if at any time a simulation does not lie within the area, it is considered discrepant. The discrepant simulations can be further divided into two subcategories. One that is above the SLR (see Fig. 2, purple and yellow tracks, top left) and one below (red track). With this classification we can identify the main parameters that drive a simulation to be located in a certain spot on the SL Diagram. It is worth mentioning that a fraction (~32%) of the observational data points lie outside the 1er 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 Sect. 3.1).

Fig. 1 Comparison between the opacity models that we used at 850 µm as a function of the maximum particle size and for a power-law size distribution with an exponent of q = 2.5. Shown is the opacity cliff (see text for details) at a wavelength of 850 µm (blue and orange shaded regions). The blue line refers to the opacity model from Ricci et al. (2010) with compact grains (labeled R10-0) and the orange line to that of Birnstiel et al. (2018) with compact grains (labeled DSHARP). The green, purple, and red lines refer to 10%, 50%, and 90% porous grains in the DSHARP model. The R10-0 and DSHARP opacity values differ by a factor of ~8.5 at the position of the opacity cliff. As the porosity of the DSHARP model increases, the opacity cliff starts to flatten out, until it completely disappears for very porous grains (90%). The location of the cliff shifts to larger particle sizes as it diminishes in porosity. The black dashed line shows the particle size at 1 mm. The value for the R10-0 model at this size corresponds to ${\kappa }_{\nu }^{\text{R10}-0}=9.8{\text{cm}}^2{\text{g}}^{-1}$, while for the DSHARP to ${\kappa }_{\nu }^{\text{DSHARP}}=4{\text{cm}}^2{\text{g}}^{-1}$.

3 Results

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

3.1 Evolution tracks

In Sect. 2.5, we explain what an evolution track is, while we show some examples in Fig. 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 explore the effect of the most important parameters of our grid model and we show how every parameter affects the evolution track on the SL Diagram in Fig. 2–4.

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 constructed 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., Fig. 5). In black we show the simulations where we used a smooth surface density profile, as in Rosotti et al. (2019a). In green we show the case where the planet/star mass ratio is q = 1 × 10⁻³ at a location of 1/3 of the r_c (inner planet), in red the same planet at a distance of 2/3 of the r_c (outer planet), and in blue two planets of q = 1 × 10⁻³ at l/3r_c and 2/3r_c. The white hatched bars 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 evolves toward smaller radii and lower luminosity as the dust drifts inward, so the emission (and size) decreases. Moreover, the opacity cliff moves farther in as the radius in the disk where the maximum particle size a_max is (i.e., the value for the peak opacity) decreases, due to radial drift and grain growth, as in Rosotti et al. (2019b).

In contrast, in the case where a planet is present the pressure bump that is formed stops the dust from drifting toward 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 Fig. 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 q = 10⁻³, 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_c, we see from the green line (number 3) that both the 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 does not decrease 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_c, we observe a much longer track similar to the one with the smooth profile (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 Sect. 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_c and one at 1/3. 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 is more dominant in the evolution track.

Fig. 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: 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 Eq. (13), which shows the relation between the luminosity and the effective radius. Finally, the blue shaded region is the area within 1σ of the SLR where the simulations are considered to be matching. The green track labeled number 1 is considered a matching simulation since it starts and ends inside the SLR. The other tracks are considered discrepant. The beginning of the track is where the empty bullet is (top right) and the end is where the number is found (bottom left). Top right: varying only the presence and the position of a Jupiter-mass planet. The red line (number 1) is for a smooth disk and the green line (number 3) for a disk with a Jupiter-mass planet at 2/3 of the rc; the purple line (number 2) corresponds to a planet at 1/3 of the rc and the yellow line (number 4) to a simulation with two planets at the positions mentioned above. Bottom left: varying only the turbulence parameter α. Higher α values lead to higher luminosity. Bottom right: varying only the characteristic radius rc. Higher rc values lead to larger and more luminous disks.

3.1.2 Effect of the turbulence α-parameter

The effect of the turbulence parameter α on the evolution track is straightforward: higher α leads to higher luminosity. In Fig. 2 (bottom left), we show the evolution tracks of a disk with the same initial conditions, varying only the α-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_c) since the effect is more prominent on these disks.

To understand the trends in Fig. 2 (bottom left), it is instructive to consider Fig. 3, where we show an example of how different α-values affect the efficiency of trapping. In the top panel, we show the dust mass local flow rate $M_{\text{acc,d}}(r)\left[M_{\oplus }{\text{yr}}^{-1}\right]$ as a function of radius. For low α-values 1 × 10⁻⁴ (red line), 5 × 10⁻⁴ (blue), and 1 × 10⁻³ (green), the mass flows toward the bump and the local flow rate for r < r_p is low (≤10⁻⁸ M_⊕ yr⁻¹), meaning that the trapping is efficient enough to stop the dust from drifting toward the star. On the other hand, for the high α-values 5 × 10⁻³ (yellow) and 1 × 10⁻² (purple), the dust mass local flow rate stays almost constant (≤10⁻⁵ M_⊕ yr⁻¹) throughout the whole disk, meaning that the bump does not trap the particles, just locally slows 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 α-values, most of the mass is located in the bump, while for α = 5 × 10⁻³ and α = 1 × 10⁻² it increases with the radius, meaning that the bump allows more grains to escape and therefore it does not contain a significant fraction of the disk mass.

Returning to Fig. 2 (bottom left), disks with low to medium α-values (10⁻⁴, 5 × 10⁻⁴, 10⁻³) are less luminous than disks with higher α-values. This is so 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. Lower α will lead to a narrower ring (Dullemond et al. 2018), and thus less emitting area and therefore lower luminosity, independently of the amount of mass in the ring. On the other hand, high α-values work against trapping in various ways, leading to more luminous disks. A higher α-value decreases the particle size in the fragmentation limit, which are less efficiently trapped by radial drift (e.g., Zhu et al. 2012). It increases the diffusivity, which allows more grains to escape the bump and it increases the viscosity in the same way, so more dust sizes are traveling with the accreting gas. Furthermore it smears out the pressure peak, causing less efficient trapping by radial drift (Pinilla et al. 2012b). Moreover, with high α-values 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 α-value is high (5 × 10⁻³, 10⁻²), the planet cannot efficiently trap the dust and the disk evolves farther along the SLR to lower luminosities. A disk that contains a planet with α = 10⁻² behaves the same as a smooth disk without a planet on the SL Diagram, implying that in this case a very massive planet (several Jupiter masses) would be needed to significantly affect disk evolution. This results in the dust grains becoming smaller since the gas turbulent velocity is increasing and the collisions between them are more destructive. Consequently, the gap is becoming shallower while there is also more diffusion.

In Fig. 5 (top left), we show the dependence of the α-viscosity parameter on the matching fraction. As expected, all simulations tend to favor low values of the turbulence parameter 10⁻⁴ ≤ α ≤ 10⁻³. Smooth simulations show a clear tendency toward low α-values because when they are drift dominated they remain in the SLR. On the other hand substructured disks show a preference toward 5 × 10⁻⁴ ≤ α ≤ 10⁻³. The dust trapping is efficient enough and allows the disks to retain their mass, but while moving to higher α-values α > 2.5 × 10⁻³, the trapping stops being efficient and leads to higher chances that the evolution track will leave the SLR in the selected time span. If α < 5 × 10⁻⁴ the dust rings become narrow and the luminosity is not large enough to place them in the SLR, and therefore the matching fraction decreases.

Fig. 3 Top panel: local flow rate of the dust mass in M⊕ yr−1 as a function of radius for a disk with a Jupiter mass planet at 31 au. For the low α-values 1 × 10−4 (red line), 5 × 10−4 (blue) and 1 × 10−3 (green), we observe that the bump outside the planet location is large enough to hold the dust from drifting toward the star, while for larger α-values 5 × 10−3 (yellow), 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 α-values 1 × 10−4, 5 × 10−4, 1 × 10−3, roughly all the mass of the disk is inside the bump that is created from the planet. For larger α-values 5 × 10−3, 1 × 10−2, the bump is minimal and not able to hold the dust from drifting.

3.1.3 Effect of the characteristic radius r_c

In Fig. 2 (bottom right), we plot the evolution tracks of a smooth disk with the same initial conditions while varying only the characteristic radius (r_c = 10, 50, and 100 au, and 180 au). 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 on the evolution tracks is straightforward. The larger the r_c the more the evolution track moves toward the top right of the plot, meaning that for a larger disk we expect higher luminosity. In more detail, increasing the characteristic radius (going from the 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 Fig. 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, substructured disks with small r_c (10 au) are mostly above the correlation because they have high luminosity relative to their size (as explained by our size-luminosity estimate in Sect. 4.3) and they are unable to enter the correlation in time. For large radii (>150 au) 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 Fig. 6).

Therefore, we observe a peak toward a specific characteristic radius, around 80–130 au when a planet is at a location of 1/3 of the r_c and around 30–80 au for the case where a planet is at a location of 2/3 of the r_c. The inner planet constrains the disk to a small size but with relatively high luminosity, therefore placing it above the SLR before 300 kyr, 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 the two pressure bumps compete with each other and each contributes in one of the ways described above.

Fig. 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 a disk that contains a planet. Bottom right: varying the fragmentation velocity of a smooth disk.

3.1.4 Effect of the disk mass - M_d

In Fig. 4 (top left), we plot the evolution tracks of a smooth disk with the same initial conditions varying only the disk mass (M_d). We choose a smooth disk to show the effect more clearly, but the same principle applies to most disks.

The disk mass contributes to both luminosity and radius. Higher disk masses lead to higher luminosities, both at the beginning and at the end of the track. Since for a fixed r_c, disks with higher M_d have higher 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 of the SLR for the lowest M_d case (M_d = 5 × 10⁻³ M_*) occurs because all grain sizes become smaller than the opacity cliff. If we choose a disk that contains a planet, the massive disks (M_d ≥ 5 × 10⁻² M_*) will still evolve toward 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.

In Fig. 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 toward the host star the luminosity decreases, allowing them in 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) cases. As we see from the yellow bars, a disk with a smooth surface density profile must be initially massive (M_d ≥ 0.025 M_*) to remain in the SLR. The probability of a smooth simulation to match is even greater 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 discuss in Sect. 3.1.5.

Fig. 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 (300 kyr − 3 Myr). Top left: dependence on the α-value. There is a preference to low α-values (10−4 ≤ α ≤ 10−3). Top right: dependence on the disk mass. There is a preference to high disk masses $\left(0.025\le \frac{{\text{M}}_{\text{d}}}{{\text{M}}_{\star }}\le 0.25\right)$. Middle left: dependence on the characteristic radius. Smooth disks do not depend on the rc. Middle right: dependence on the stellar mass. There is a slight preference to higher values. Bottom left: dependence on the fragmentation velocity. There is a slight preference to higher values.

Fig. 6 Heat maps of simulations with the R10-0 opacities. From left to right, the three columns represent the smooth case, a planet at 1/3rc, and a planet at 2/3rc. From top to bottom, the rows represent three different snapshots at 300 kyr, 1 Myr, and 3 Myr. The white solid line 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 (reff ~ 10 au) where observational results are available.

3.1.5 Effect of different opacity and grain porosity

In Fig. 4 (top right), we represent the behavior of several similar tracks varying only the opacity model. The simulations with the Ricci et al. (2010) opacity model (R10-0, blue line) produce more luminous and larger disks than those with the DSHARP model (orange line) 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-porous and very porous grains (DSHARP-50 and DSHARP-90, purple and red line, respectively) 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.

In all the histogram figures (Fig. 5), the same trend stands for either the opacity from R10-0 (solid color bars) or 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 are only a few simulations that match, hence the hatched bars are barely visible. For smooth disks the total matching fraction is 29.6% with the R10-0 opacity, while it is 0.8% with the DSHARP value. For disks with an inner planet the matching fractions are 30.2% and 15.9%, respectively (see Sect. 4.3, where we explore the overall impact of porous grains for the entire grid of models).

3.1.6 Effect of the stellar mass - M_*

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

In Fig. 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 0.2 M_⊙, 0.6 M_⊙, 1.0 M_⊙, and 2.0 M_⊙. As expected, the highest value of M_* = 2.0 M_⊙ leads to the largest and most luminous disk (green line), while the opposite is true for M_* = 0.2 M_⊙ (red line). There is similar behavior for smooth disks, but in that case the radius of each disk is much smaller due to radial drift.

Furthermore, the stellar mass is the least important parameter when defining whether a simulation matches. The histogram in Fig. 5 (middle right) confirms this. Even though the trend shows that using higher stellar mass has a greater matching fraction, this occurs because the stellar mass scales to the disk mass, and higher disk mass leads to more matching simulations (see Sect. 3.1.4).

This scaling implies that the luminosity (L_mm) scales with the stellar mass (M_*). Our models follow a relation that is not as steep as the observed $L_{\text{mm}}\propto M_{\ast }^{1.5}$ in Andrews et al. (2018b), because there is not a correlation between disk size and stellar mass in our simulations (see Sect. 4.5 and Fig. 10 where we explore further this relation).

3.1.7 Effect of the fragmentation velocity - v_frag

In Fig. 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 cm s⁻¹, 600 cm s⁻¹, 1000 cm s⁻¹, and 2000 cm s⁻¹. We observe that for medium and high values of υ_frag (in this case for υ_frag ≥ 600 cm 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 becomes too low, which leads to a higher luminosity in the first snapshots.

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 Fig. 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 υ_frag ≥ 600 cm 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 υ_frag, most of the simulations are fragmentation-limited and they lose luminosity relatively quickly, moving them out of the SLR. Low values of υ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, in the recent lab experiment of Blum (2018), 100 cm s⁻¹ ≤ υ_frag ≤ 1000 cm s⁻¹ is considered a value consistent with lab work and υ_frag > 1000 cm s⁻¹ a 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 like the one we have it is not easy to extract global results. Since there are too many tracks to plot all of them in 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 kyr, 1 Myr, 3 Myr).

In Fig. 6, we plot three different cases. In Col. 1, we show the smooth case, in the Col. 2 the case where we use a planet/star mass ratio q = 10⁻³ at 1/3 of the r_c, and in the last column q = 10⁻³ at 2/3 of the r_c. The snapshots are at three times as shown, as are the SLR and the standard deviation. The red line is our prediction for the cases where we include a planet (see Sect. 3.2.2). Instead of following the relation from Andrews et al. (2018b), they seem to follow a relation of $L_{\text{mm}}\propto r_{\text{eff}}^{5/4}$. In Sects. 3.2.1 and 3.2.2, we perform a more detailed analysis on the this topic.

We observe that most of the disks start inside and above the correlation (first row at 300 kyr). In the smooth case the disks lose a great deal 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 (29.6%) of simulations occupying and following the SLR. As we explain in Sect. 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 do not decrease in size, but mainly in luminosity. This is due to the formation of the pressure bump that keeps the dust from drifting farther in, thus 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 in the diagram. The clustering in r_eff that is formed in the plots with a planet are an artifact of our parameter grid. A randomly chosen value of r_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 (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 (20.6% of the disks match).

The difference becomes striking when we use the dust opacities from DSHARP in Fig. 7. Since the DSHARP opacity is lower than the R10-0 (see Fig. 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 Myr) and only (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 (> 11.1%) of the disks match depending on the model.

A similar behavior for the cases where we include strong substructures is seen for the opacity model with semi-porous grains, DSHARP-50, in Fig. 8. Even though there are fewer matching simulations in total, if the substructures are large enough, they are able to keep a significant number of simulations in the clusters on the SLR (>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 (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 is true looking at the heat map (Fig. D.2) 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 (1.3%), while a similar fraction of sub-structured disks match, as in the DSHARP-50 case (>10.2%). This heat map is included in Appendix D.

From these heat maps we can extract three important results:

• 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;

• Whether a smooth disk matches or not depends heavily on the opacity model. Birnstiel et al. (2018) DSHARP opacities produce significantly fewer simulations in the SLR than the Ricci et al. (2010) R10-0 model and only a fraction of the simulations match the semi-porous and very porous grains and the model DSHARP-50 and DSHARP-90. Therefore, the porosity should be lower than 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 discussed in Sect. 4, the observed correlation can be a mixture of smooth and substructured disks that adds scatter to the simulated SLR;

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

Fig. 7 Heat maps of simulations with the Birnstiel et al. (2018) opacities. From left to right, the three different columns represent the smooth case, a planet at 1 /3rc, and a planet at 2/3rc. From top to bottom, the rows represent three different snapshots at 300 kyr, 1 Myr, and 3 Myr. The white solid line 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 (reff ~ 10 au) 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 Figs. 6 and 7), we 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 versus 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 Sect. 3.2.2, we estimate the SLR based on these empirically determined widths. For a complete derivation of the two sections, we refer the reader to Appendices A.1 and 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 a-values (see Appendix C.1). 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 α-parameter, $${\sigma }_{\text{g}}=C\cdot h^{\text{p}}\cdot q^k\cdot {\alpha }^l,$$(14)

where C is a constant, h is the scale height, q the mass ratio, and α the turbulence parameter. We find that the width in the measured range scales approximately as $${\sigma }_{\text{g}}\propto h^{\text{0}\text{.81}}\cdot q^{0.14}\cdot {\alpha }^{0.05}.$$(15)

Fig. 8 Heat maps of simulations with the Birnstiel et al. (2018) D-50 opacities with 50% porostiy. From left to right, the three columns represent the smooth case, a planet at 1 /3rc, and a planet at 2/3rc. From top to bottom, the rows represent three different snapshots at 300 kyr, 1 Myr, and 3 Myr. The white solid line 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 (reff ~ 10 au) where observational results are available.

3.2.2 Size-luminosity relation of disks with companions

In Fig. 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 $L_{\text{mm}}\propto r_{\text{eff}}^{5/4}$. If we assume that all the luminosity of a disk comes from rings that are approximately optically thick, we can approximate as $$L\simeq A\cdot B_v,$$(16)

where A is the area and B_v 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 A ∝ r · σ_d, where r is the radius and σ_d is the width of the pressure bump in the dust and there is linear scaling of σ_d with h. The width of the dust ring depends on the width of the gas, $${\sigma }_{\text{d}}\propto {\sigma }_{\text{g}}\cdot \sqrt{\frac{\alpha }{\text{St}}},$$(18)

as in Dullemond et al. (2018), where 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 in Sect. 3.2.1, we find that the luminosity (L_mm) scales with the radius as $$L_{\text{mm}}\propto r_{\text{eff}}^{5/4}$$(19)

which is the relation that we plot with the red line in Figs. 6–8, and D.2. We find that this theoretical estimate nicely explains the size-luminosity scaling seen when strong substructure is present. However, toward larger radii, this relation slightly overpredicts the luminosity. For example, we could get a shallower slope looking at the heat map because toward large radii our fitting line is above the main bulk of the simulations. For a complete derivation of the two sections, we refer the reader to Appendices A.1 and A.2.

4 Discussion

To summarize the results discussed above, we explored the observed trend among the (sub-)mm disk continuum luminosity (L_mm) and the 68% effective radius (r_eff) of protoplanetary disks. Following the size-luminosity relation (SLR) obtained from Tripathi et al. (2017) and Andrews et al. (2018b), $L_{\text{mm}}\propto r_{\text{eff}}^2$, we showed which initial conditions are favorable for a disk to remain on the SLR for a time span of 300 kyr - 3 Myr. 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 substructures 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 in 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 a, the initial disk mass 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 (α ≤ 10⁻³) and high disk mass (M_d ≥ 2.5 × 10⁻² M_*), and they are affected strongly by the existence of a strong trap (in this study caused by a giant planet).

Turbulence-α values greater than 10⁻³ 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, are less efficiently trapped. Finally, the dust trap is not as pronounced if the alpha viscosity is higher. All of this acts in concert to make dust trapping ineffective, and causes the disks to behave as if they were smooth (see Fig. 3 in Sect. 3). If the fragmentation velocity is high 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 lower luminosities while they evolve up to 3 Myr 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 higher luminosity (see Sect. 4.3). On the other hand, most of the less massive and smaller disks end up below the correlation at lower 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 toward lower luminosity and radius. Disks with a massive planet in the inner part of the disk (1/3 r_c) have more extended evolutionary tracks that are overall less luminous. In contrast, disks with a planet in the outer part (2/3 r_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 presence of planets increases the fraction of matching simulations on the SLR, but this result is also a function of the opacity (see Sect. 4.3).

4.2 Position along the SLR

The position of a disk along the SLR is determined mainly by the disk mass M_d and the disk size r_c, as can be seen in Fig. 2 and Fig. 4. More massive and large disks are located in the top right part of the SL diagram, while small and less massive are in the middle and left parts. In Fig. 9 we show the kernel density estimate (kde) of the luminosity for all matching simulations for four different cases and three different snapshots, and also plot the observational kde from Andrews et al. (2018b). The brightest disks that stay on the SLR are those that contain planets that are located at the outer part of the disk (2/3 of the r_c, yellow and green lines). A planet in the outer region leads to larger, more luminous disks as explained in Sects. 2.3 and 4.1. Massive planets at 1 and 3 Myr 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 these disks evolve similarly to a smooth disk.

4.3 Opacity model preferences

We used different opacity models for our study. As a base model, we made use of the composition of Ricci et al. (2010) opacities (R10-0), but for compact grains (i.e., no porosity) as in Rosotti et al. (2019a). 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 (<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 toward 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 recall that the opacity at our wavelength is a factor of ~8.5 higher in the R10-0 case compared to DSHARP at the opacity cliff location (~0.1–1 mm), 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 disks with weak substructures (where a disk behaves as smooth). If this is the case then only compact grains can explain the SLR; instead, when substructures are strong then any of the opacity models and the porosities tested in this work can explain the substructured SLR. The latter applies because most of the substructures become optically thick.

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

Fig. 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 300kyr – 3 Myr. The black line refers to the disks from the Andrews et al. (2018b) sample. Disks with planets have higher luminosity, while smooth disks have low luminosity at 3 Myr. When two planets are included the luminosity is higher than with a single planet.

4.4 Types of disks on the SLR

We show 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 Appendix A.1), we derived a theoretical prediction for disks with substructures (Sect. 3.2.2). Smooth disks with compact or slightly porous grains seem to follow the Andrews et al. (2018b) relation $L_{\text{mm}}\propto r_{\text{eff}}^2$, while disks with massive planets follow a relation $L_{\text{mm}}\propto r_{\text{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 they might not show strong enough and optically thick substructures, such as AS209 or HD163296 (see Huang et al. 2018a). Figure 6 shows how smooth disks follow the SLR, while disks with strong substructures follow a different relation but intersect 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 for USco, the oldest region (~ 10 Myr). The observed SLR therefore might flatten with age of the region. We could examine these results since we evolve our simulations to 10 Myr, 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 <20%, while in our models it would be 100%.

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

The SLR can be explained if there is a mixture of both smooth and strong substructured disks. Smooths disks always follow the SLR, as shown in Fig. 6, while the bright substructured disks populate the upper right part of the SLR (Figs. 6, 7, and D.2). Disks with substructures that have large r_c and low disk mass M_d populate the lower right part of the plot under the SLR. These disks are not favored by Andrews et al. (2018b), who find 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 substructured disks that occupy the upper right part and smooth (or weakly substructured) disks that occupy the lower left part of the SLR. Our results seem to be in agreement with the observational classification of van der Marel & Mulders (2021) who suggest that all bright disks should have substructures formed by giant planets. Moreover, the SLR for the substructured disks is independent of the opacity model, but it slightly overpredicts the luminosity for the very large disks (see Sect. 4.3 for more about the opacity).

4.5 L_mm - M_* relation

In Sect. 3.1.6 we discuss that the stellar mass (M_*) is directly correlated with the disks mass (M_d) and that the disk temperature is only a weak function of (and therefore M_*). The fact that the disk mass scales with the stellar mass implies that the luminosity (L_mm) scales with the stellar mass. In Fig. 10, the L_mm − relation is shown for three different models for all matching simulations: the smooth case (yellow lines), a planet with planet/star mass ratio q = 10⁻³ at 1/3r_c (green lines), and a planet with the same mass ratio at 2/3r_c (red lines). The markers define the median value of the luminosity at 1 Myr and the error bars are the 75% percentile from the upper and lower value. The blue line is the $L_{\text{mm}}\propto r_{\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 for the Andrews et al. (2018b) model, but the cases with strong substructures have steeper profiles than the smooth cases. 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) were imposed, the mass-luminosity relation is expected to steepen as disks with optically thick substructures 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 Fig. B.2. In the panel r_c − M_d, shown as 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 (M_d < 10⁻² M_*) and be large in size (r_c > 150 au). Moreover the turbulence parameter should be relatively small α ≤ 10⁻³, otherwise the disk would behave as smooth and would follow the SLR.

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

4.6 Scattering

Scattering is included in our simulations, as introduced in Sect. 2.4. Compared to the case where only the absorption opacity is used, the difference is minimal and it can be observed only in a few cases. With the inclusion of scattering, the originally brightest disks (above the SLR) tend to move toward lower luminosity (move down in the SL diagram). This happens for disks that are optically thick, hence for those that contain planets. This effect favors the SLR and allows slightly more disks (~2%) to enter the selected region. However, for moderately optically thick disks the emission is larger and a small fraction of disks move up (toward higher luminosity) on the SL diagram (Fig. 4 in Birnstiel et al. 2018). This happens because the derived intensity (Eq. (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 (τ ≪ 1) the effect of scattering is insignificant because the intensity (Eq. (11)) approaches $I_{\nu }^{\text{out}}\to {ϵ}_{\nu }^{\text{eff}}{B}_{\nu }(T_{\text{d}})\Delta \tau /\mu$, which is the identical solution to when ${\kappa }_{\nu }^{\text{out}}$ is set to zero while ${\kappa }_{\nu }^{\text{abs}}$ is kept unchanged, as also shown in Birnstiel et al. (2018).

The effect of scattering also depends on the albedo $({\eta }_{\nu }=1-{ϵ}_{\nu }^{\text{eff}})$. For the compositions we use, the maximum effective albedo is 0.57 for R10-0 and 0.82 for DSHARP, while it can reach ~0.97 for DSHARP-90. For these compositions the effect of scattering is never more than a factor of ~ 1.7 at a particle size of 1 mm. The result of scattering is higher if we increase the 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 the authors discuss that a completely optically thick disk with high albedo (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 substructures create only rings rather than disks that are completely optically thick everywhere.

4.7 Predictions for longer wavelengths

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

We performed a series of simulations at 3.1 mm for a comparison with these results. The disks are fainter and smaller at 3.1 mm. Two effects contribute to this. First, 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 µm case. Second, the intensity at 3.1 mm is less according to Plank’s spectrum. Therefore, the luminosity will be lower.

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 toward lower luminosities and smaller radii. On the other hand, the substructured 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 Fig. 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 λ = 850 µm, but at 3.1 mm 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 substructured and large smooth disks. The heat map in Appendix D confirms this. In this figure we plot the simulations at 3.1 mm, using the R10-0 opacity model and we overplot the SLR from Tazzari et al. (2021): $L_{\text{mm}}\propto r_{\text{eff}}^{1.2}$. Substructured disks can explain this relation very well since it is similar to the scaling relation we calculated for disks with strong substructures in Sect. 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 mm.

We should mention the possibility that the flatter relation can be due to observational bias toward large disks, which tend to be substructured. 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 investigate this possibility further.

Fig. 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 300 kyr-3 Myr. From left to right, the smooth case, a planet with planet/star mass ratio at 1/3rc, and a planet with the same ratio at 2/3rc. Different limits are used on the x-axis to highlight the evolution of the dust-to-gas ratio. The initial dust-to-gas ratio is 0.01. For the smooth case the dust-to-gas ratio decreases by three orders of magnitude up to 3 Myr. 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 Myr where the ratio is higher than 0.01. The gas mass moves faster than the dust mass in this case.

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 kyr to 3 Myr. This does not exclude the possibility that some disks with high disk mass might evolve a great deal on the SLR diagram for 10 Myr–20 Myr. Disk dissipation has not been modeled in this paper, but it will be considered for future work. In Fig. 11 we show the kde of the global dust-to-gas ratio² for three different snapshots between 300 kyr – 1 Myr – 3 Myr 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 strong trap. In the second panel there are cases where the dust-to-gas ratio increases over the initial 0.01. These are sub-structured disks with intermediate α-values, high fragmentation velocities (>1000 cm s⁻¹), and small sizes (<60 au). The gas is removed more quickly than the dust, leading to a higher dust-togas ratio. Large values of α would lead to less trapping and the dust would drift as usual, while low α would mean that the disk does not evolve significantly.

As mentioned in Sect. 2, the stellar luminosity is not evolved in the simulation. If this were the case the disk luminosity would scale approximately linearly with the stellar luminosity and would be further modulated by resulting changes in the dust evolution. We therefore expect a general shift of the disks toward lower luminosities, but with the trends that have explored in Sect. 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 α-values and lower fragmentation velocities than the values shown before. An example of a heat map in shown in Appendix D.1.

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 timescales in the outer part of the disk are much longer and the timescale 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 planet later, as has been suggested (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.

5 Conclusions

In this paper we 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. Disks with strong traps (i.e., massive planets) follow a different SLR than smooth disks. Smooth disks follow the Andrews et al. (2018b) relation, $L_{\text{mm}}\propto r_{\text{eff}}^2$, as shown by Rosotti et al. (2019a), while disks with massive planets follow $L_{\text{mm}}\propto r_{\text{eff}}^{5/4}$. This could mean that not all disks in the Tripathi et al. (2017) and Andrews et al. (2018b) joint sample have a substructure as significant as HD163296, for example. We explained this result with a simple analytical derivation and we find 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 $L_{\text{mm}}\propto r_{\text{eff}}^2$

2. Whether disks follow the SLR depends heavily on the opacity model. When the DSHARP (Birnstiel et al. 2018) opacity is used, disks are not as luminous in the first 300 kyr 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 (0.8% with DSHARP and 29.6% with R10-0). Therefore, with this opacity model, only disks with substructures can populate the SLR. On the other hand, R10-0 opacities can reproduce disks both with and without substructures since the absolute value of the opacity at 850 µm is ~8.5 times higher than DSHARP for particle sizes around ~0.1 mm (position of the opacity cliff) and the disks become luminous enough to enter the relation;

3. The SLR is more widely populated when substructures are included in contrast to a tight correlation for smooth disks. Substructured 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 substructured disks;

4. The grain porosity can drastically affect the evolution track of the disk. Throughout our models, relatively compact grains (<50% porosity) are preferred for simulations that follow the SLR. If we use slightly porous grains (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-porous (50%) and porous grains (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 substructures are included;

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 ~300 kyr. By this time it will enter the relation, and depending on the other initial conditions it will either remain there and be considered a matching simulation or will leave it;

6. There is a preference toward low a-values (lower than 10⁻³). This result is in line with other more direct methods of determining a (e.g., Flaherty et al. 2018). There are multiple reasons for this tendency. For α ≥ 2.5 × 10⁻³, disks tend to be more fragmentation dominated, the particle size decreases, and consequently they 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 a is low, the ring that is formed becomes too narrow and disks tend to have lower luminosity;

7. The location of the planet as a function of the characteristic radius plays a major role in the final outcome. If a planet is included in the inner part of the disk (1/3rc), the disk has to be significantly larger in order to retain the correct ratio of luminosity to effective radius to stay in the SLR. Instead, when an outer planet (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 the two defines the luminosity. These results are also affected by the opacity model;

8. We expect a less extended evolution track when substructure is included. The pressure bump halts the dust from drifting farther 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. We are not able to construct optically thick disks with high albedo (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 substructure create only optically thick rings rather than a uniform optically thick distribution.

10. We chose different gap profiles based on Kanagawa et al. (2016) and compared them again to 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. Instead, the width of the gap is the important parameter (see Appendix C, where we compare the different profiles for different parameters.)

This study shows how the combination of observed and simulated populations allows us to put constraints on crucial unknowns such as the disk turbulence, or the dust opacity. Future work is required to also investigate the effects of disk buildup and dissipation as well as planet migration and planetesimal formation.

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).

Appendix A Derivation of gap width and SLR

In this section we show the complete derivation of the width versus radius r, planet/star mass ratio q, and scale height h relation, as explained in Sect. 3.2.1, and the derivation of the SLR of disks with strong substructures, as shown in Sect. 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 α-values (see Appendix C.1). 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\left(x\right)=\frac{1}{\sqrt{2\pi }{\sigma }_{\text{g}}}{e}^{-\frac{{\left(x-x_0\right)}^2}{2{\sigma }^2}},$$(A.1)

where x₀ 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 }_{\text{g}}}\left(-\frac{x-x0}{{\sigma }^2}\right)e^{-\frac{{\left(x-x_0\right)}^2}{2{\sigma }_{\text{g}}^2}},$$(A.2)

which gives us $$\frac{\partial lny}{\partial lnx}=-x\frac{x-x0}{{\sigma }_{\text{g}}^2}.$$(A.3)

If we are close to the pressure maximum, we can approximate that x → x₀ + δx, and the equation yields $$\frac{dlny}{dlnx}=-\frac{x_0\delta x}{{\sigma }_{\text{g}}^2}.$$(A.4)

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 }_{\text{g}}^2}.$$(A.5)

Solving for σ_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}}}.$$(A.6)

We measured all the widths of our hydrodynamical simulations, and we investigated how they scale with radius, scale height, and the α-parameter. We fit the width of the pressure bump as multiple power-law as $${\sigma }_{\text{g}}=C\cdot h^p\cdot q^k\cdot {\alpha }^l,$$(A.7)

where C is a constant, h is the scale height, q the mass ratio, and α the turbulence parameter. We find that the width in the measured range scales approximately as $${\sigma }_{\text{g}}\propto h^{\text{0}\text{.81}}\cdot q^{0.14}\cdot {\alpha }^{0.05}.$$(A.8)

A.2 Derivation of SLR of disks with companions

To understand the overall shape of the heat map for the case with strong substructures (i.e., the red lines in Fig. 6 and Fig. 7), we present here the complete derivation of the theoretical estimate from Sect. 3.2.2. This estimate is based on the empirically determined widths (see Appendix A.1), as explained in Sect. 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_v,$$(A.9)

where A is the area and B_v is the Planck function. If we assume the Rayleigh–Jeans approximation, then we can substitute the Planck function with the temperature T and the equation becomes $$L\propto A\cdot T.$$(A.10)

We make the assumption that the area of the pressure bump scales as A ∝ r · σ_d, where r is the radius and σ_d is the width of the pressure bump in the dust and there is linear scaling of σ_d with h. The temperature is T ∝ r^−1/2 and h ∝ r^5/4. Combining the above relations, we obtain the theoretical expression for the luminosity: $$L_{\text{mm}}\propto r^{7/4}.$$(A.11)

In the previous section we determined how the width of the pressure maximum depends on the scale height σ_g ∝ h^0.81, so we need to find a relation between the σ_g and σ_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 }_{\text{d}}\propto {\sigma }_{\text{g}}\cdot \sqrt{\frac{\alpha }{\text{St}}},$$(A.12)

as in (Dullemond et al. 2018), where St is the Stokes number. Moreover, $$St\propto \frac{1}{\alpha \cdot c_{\text{s}}^2},$$(A.13)

where c_s is the sound speed. Combining these two, we obtain that $$h_{\text{d}}\propto h_{\text{g}}\cdot \alpha \cdot c_{\text{s}}.$$(A.14)

The sound speed is $$c_{\text{s}}\propto r^{-1/4}.$$(A.15)

Using all the above equations we conclude that $${\sigma }_{\text{d}}\propto {\sigma }_{\text{g}}\cdot r^{-1/4},$$(A.16)

or, if we express the radius as a function of the scale height, $${\sigma }_{\text{d}}\propto {\sigma }_{\text{g}}\cdot h^{-1/5}.$$(A.17)

The scale height is h = c_s/Ω, where $\Omega =\sqrt{G{M}_{\star }/r^3}$ the Keplerian angular velocity, leading to a dependency h ∝ r^5/4. Therefore, we can use Eq. (A.8) and express the calculated σ_g as a function of radius: $${\sigma }_{\text{g}}\propto r^{1.01}.$$(A.18)

Finally, we can calculate the dust width σ_d as a function of both the scale height and radius. Substituting into Eq. (A.17) we obtain $${\sigma }_{\text{d}}\propto h^{0.61},$$(A.19)

and into Eq. (A.16) we obtain $${\sigma }_{\text{d}}\propto r^{0.75}.$$(A.20)

If we calculate again how the luminosity scales and replace the scaling of the dust width σ_d ∝ r^0.75, we find that $$L_{\text{mm}}\propto r^{5/4},$$(A.21)

which is the value that we fit with the red line in Fig. 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 in the upper part of the plot. These visualizations contain only the simulations that match the SLR as the histograms in Fig. 5.

In Fig. B.1 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 α 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 ∝ $\propto \sqrt{\alpha }$ nicely explains the correlation: simulations above this relation are strongly fragmentation-limited, and therefore do not follow the SLR (see Rosotti et al. 2019b). The bottom left cell (corresponding to α = 10⁻⁴ and υ_frag = 200 cm/s) is empty because it is fragmentation-limited and increasing the α-value only makes the disks more fragmentation-dominated. As the fragmentation velocity increases then the allowed α-value increases as well.

In Fig. B.2 and Fig. B.3, we show the corner plots of the cases where a planet with a planet/star mass ratio q = 10⁻³ is included at 1/3 and 2/3 of the r_c, respectively. In both cases we obtain similar results. In the second column a mild covariance can be seen between the turbulence parameter α and the disk mass. Most of the simulations are toward low α-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 α and r_c (for a limited α-range) is shown. With vertical white dots we indicate the α = 5 × 10⁻³ value, where there is a gradual boundary from matching to discrepant simulations. This relation is visible in the second and third columns of Fig. 6. In summary, the leftmost clump of simulations above the correlation is the case where r_c = 10 au. This clump is mostly above the correlation (unless α = 10⁻⁴), and it explains the empty row for r_c = 10 au in the corner plots of Fig. B.2 and Fig. B.3.

With low α-values there is more trapping in the pressure bump and the emission is optically thick, hence the disks are brighter. 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 α-values. This range depends on r_c. Moreover, when α increases, the emission is not optically thick anymore because the dust gets through the trap. More specifically, for α ≥ 5 × 10⁻³ (as seen in Fig. 3), the bump can no longer efficiently trap and this is why this transition from matching to discrepant simulations occurs.

Therefore, in the case where a planet is included the strongest correlation is between α-value and the characteristic radius r_c. The area of the parameter space with the highest fraction of matching simulations is where we have small disks (around ~50 au) and low α-values. This allows us once again to constrain the disk characteristic radius to a value that is not larger than 100 au.

Fig. B1 Corner plot of the smooth case. From left to right: Fragmentation velocity υfrag[cm/s], turbulence α-value, disk mass Md[M*], characteristic radius rc[au]. From top to bottom: α, Md[M*], rc[au], M*[M⊙]. In the upper left plot the white dots show the track where the distinction between the fragmentation- and the drift-limited regime is. In the middle plot (rc − α) the white dots show the α-value where the pressure bump cannot hold the dust efficiently anymore. In the panel rc - Md, the mean value of the characteristic radius for each disk mass is shown as a white dot.

Fig. B.2 Corner plot with a planet with of a planet/star mass ratio 10−3 at a location of the 1/3rc. From left to right: Fragmentation velocity υfrag[cm/s], turbulence a-value, disk mass Md[M*], characteristic radius rc[au]. From top to bottom: α, Md[M*], rc[au], M*[M⊙]. In the upper left plot the white dots show the track where the distinction between the fragmentation and the drift limited regime is. In the middle plot (rc − α) the white dots show the α-value where the pressure bump cannot hold the dust efficiently anymore. In the panel rc- Md, shown as white dots, is the mean value of the characteristic radius for every disk mass.

Fig. B.3 Corner plot where we have a Jupiter-mass planet at a location of 2/3rc. From left to right: Fragmentation velocity υfrag[cm/s], turbulence α-value, disk mass Md[M*], characteristic radius rc[au]. From top to bottom: α, Md[M*], rc[au], M*[M⊙]. In the upper left plot the white dots show the track where the distinction between the fragmentation- and the drift-limited regime is. In the middle plot (rc − α) the white dots show the α-value where the pressure bump cannot hold the dust efficiently anymore. In the panel rc - Md, shown as white dots, is the mean value of the characteristic radius for every disk mass.

Appendix C Gap profiles

As discussed in Sect. 2, we compared the planet profile from Kanagawa et al. (2016) to hydrodynamical simulations using FARGO (Benítez Llambay & Masset 2015). We varied the turbulence parameter α, the planet/star mass ratio q, and the position of the planetary gap r_p. 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 kyr. The parameters that we used for this test are summarized in the following table.

In Figures C.1, C.2, and C.3 we show the comparison between the gap that is obtained from the hydrodynamical simulation and the profile from Kanagawa et al. (2016). In Fig. C.1 we compare the gap profile for a small planet with planet/star mass ratio of q = 3 × 10⁻⁴ for the three different values of the α-parameter. From the low α-value (first row, first column, α = 10⁻⁴), we observe that the Kanagawa et al. (2016) profile fits the width of the gap reasonably well, but there is an offset for the depth. As we look at their evolution tracks we do not observe any major difference. For α = 10⁻³ (second row, first column) we have a good fit for the depth and the width in Fig. C.2 and C.3. On the evolution track there is a small fluctuation, but at the end of the simulation they converge to the same point. For the high α-value (third row, first column) the profile nicely fits 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 perfectly fit the bottom of the gap, the difference in the evolution tracks is minimal. The track can be affected only if the width of the gap is not fitted properly. Therefore, we conclude that our chosen profile will not affect the main point of the paper, which is based on the observable quantities.

Fig. C.1 Comparison between a hydrodynamical and our dust evolution simulation for different α-values and for a planet/star mass ratio q = 3 × 10−4 at 30 au. In the left column are the different gap profiles for different α-values and in the right column the corresponding evolution tracks. Even though the profile from Kanagawa et al. (2016) does not overlapping with that obtained from the hydrodynamical simulation, the tracks in the right column produce similar results.

Fig. C.2 Comparison between a hydrodynamical and our dust evolution simulation for different α-values and for a planet/star mass ratio q = 10−3 at 30 au. In the left column are the different gap profiles for different α-values and in the right column the corresponding evolution tracks. Even though the profile from Kanagawa et al. (2016) does not overlapping with that obtained from the hydrodynamical simulation, the tracks in the right column produce similar results.

Fig. C.3 Comparison between a hydrodynamical and our dust evolution simulation for different α-values and for a planet/star mass ratio q = 3 × 10−3 at 100 au. In the left column are the different gap profiles for different α-values and in the right column the corresponding evolution tracks. Even though the profile from Kanagawa et al. (2016) does not overlapping with that obtained from the hydrodynamical simulation, the tracks in the right column produce similar results.

Appendix D Additional heat maps

In this section we add the additional heat maps where the 10% porosity (DSHARP-10) and 90% porosity (DSHARP-90) opacity is used. We plot the position of every simulation for three different snapshots (300 kyr, 1 Myr, 3 Myr). For DSHARP-10 (Fig. D.1) we obtain similar results to the DSHARP case Fig. 7 (with no porosity) as the two opacities are similar to one another (see Fig. 1). For the DSHARP-90 case (Fig. D.2), the complete absence of the opacity cliff does not allow a considerable amount of smooth disks to enter the SLR, while the same fraction of substructured disks match, as in the DSHARP-50 Fig. 8 case.

D.1 Changing the stellar luminosity

In Sect. 4.8 we performed a test where we used the luminosity of the star at 3 Myr instead of 1 Myr, which is used throughout this work. Since the stellar luminosity decreases from 1 Myr to 3 Myr, there is an overall movement of disks to lower luminosities that is shown in Fig. D.4 in the right panel. In this example we plot the heat map of a case where there is a planet at 1/3r_c at a snapshot of 1 Myr.

Fig. D.1 Heat maps of representative simulations with the Birnstiel et al. (2018) D-10 opacities with 10% porosity. From left to right, the three columns represent the smooth case, a planet at 1 /3rc, and a planet at 2/3rc. From top to bottom, the rows represent three different snapshots at 300 kyr, 1 Myr, and 3 Myr. The white solid line is the SLR from Andrews et al. (2018b) and the red solid line our fit for the cases where we include a planet. The color bar shows the number of simulations in a single cell. The blue dash-dotted line shows the minimum limit (reff ~ 10 au) where observational results are available.

Fig. D.2 Heat maps of representative simulations with the Birnstiel et al. (2018) D-90 opacities with 90% porosity. From left to right, the three columns represent the smooth case, a planet at 1 /3rc, and a planet at 2/3rc. From top to bottom, the rows represent three different snapshots at 300 kyr, 1 Myr, and 3 Myr. The white solid line is the SLR from Andrews et al. (2018b) and the red solid line our fit for the cases where we include a planet. The color bar shows the number of simulations in a single cell. The blue dash-dotted line shows the minimum limit (reff ~ 10 au) where observational results are available.

Fig. D.3 Heat maps of representative simulations with the R10-0 opacity at 3.1 mm for a direct comparison with Tazzari et al. (2021). From left to right, the three columns represent the smooth case, a planet at 1/3rc, and a planet at 2/3rc. From top to bottom, the rows represent three different snapshots at 300 kyr, 1 Myr, and 3 Myr. The white solid line is the SLR from Tazzari et al. (2021). The color bar shows the number of simulations in a single cell. The blue dash-dotted line shows the minimum limit (reff ~ 10 au) where observational results are available.

Fig. D.4 Heat map comparison between disks where the stellar luminosity at 3 Myr (right panel) is used instead of 1 Myr (left panel). There is an overall shift toward lower luminosities, but the trend remains the same.

References

All Tables

All Figures

Fig. 1 Comparison between the opacity models that we used at 850 µm as a function of the maximum particle size and for a power-law size distribution with an exponent of q = 2.5. Shown is the opacity cliff (see text for details) at a wavelength of 850 µm (blue and orange shaded regions). The blue line refers to the opacity model from Ricci et al. (2010) with compact grains (labeled R10-0) and the orange line to that of Birnstiel et al. (2018) with compact grains (labeled DSHARP). The green, purple, and red lines refer to 10%, 50%, and 90% porous grains in the DSHARP model. The R10-0 and DSHARP opacity values differ by a factor of ~8.5 at the position of the opacity cliff. As the porosity of the DSHARP model increases, the opacity cliff starts to flatten out, until it completely disappears for very porous grains (90%). The location of the cliff shifts to larger particle sizes as it diminishes in porosity. The black dashed line shows the particle size at 1 mm. The value for the R10-0 model at this size corresponds to ${\kappa }_{\nu }^{\text{R10}-0}=9.8{\text{cm}}^2{\text{g}}^{-1}$, while for the DSHARP to ${\kappa }_{\nu }^{\text{DSHARP}}=4{\text{cm}}^2{\text{g}}^{-1}$.

In the text

Fig. 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: 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 Eq. (13), which shows the relation between the luminosity and the effective radius. Finally, the blue shaded region is the area within 1σ of the SLR where the simulations are considered to be matching. The green track labeled number 1 is considered a matching simulation since it starts and ends inside the SLR. The other tracks are considered discrepant. The beginning of the track is where the empty bullet is (top right) and the end is where the number is found (bottom left). Top right: varying only the presence and the position of a Jupiter-mass planet. The red line (number 1) is for a smooth disk and the green line (number 3) for a disk with a Jupiter-mass planet at 2/3 of the rc; the purple line (number 2) corresponds to a planet at 1/3 of the rc and the yellow line (number 4) to a simulation with two planets at the positions mentioned above. Bottom left: varying only the turbulence parameter α. Higher α values lead to higher luminosity. Bottom right: varying only the characteristic radius rc. Higher rc values lead to larger and more luminous disks.

In the text

Fig. 3 Top panel: local flow rate of the dust mass in M⊕ yr−1 as a function of radius for a disk with a Jupiter mass planet at 31 au. For the low α-values 1 × 10−4 (red line), 5 × 10−4 (blue) and 1 × 10−3 (green), we observe that the bump outside the planet location is large enough to hold the dust from drifting toward the star, while for larger α-values 5 × 10−3 (yellow), 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 α-values 1 × 10−4, 5 × 10−4, 1 × 10−3, roughly all the mass of the disk is inside the bump that is created from the planet. For larger α-values 5 × 10−3, 1 × 10−2, the bump is minimal and not able to hold the dust from drifting.

In the text

	Fig. 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 a disk that contains a planet. Bottom right: varying the fragmentation velocity of a smooth disk.
In the text

Fig. 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 (300 kyr − 3 Myr). Top left: dependence on the α-value. There is a preference to low α-values (10−4 ≤ α ≤ 10−3). Top right: dependence on the disk mass. There is a preference to high disk masses $\left(0.025\le \frac{{\text{M}}_{\text{d}}}{{\text{M}}_{\star }}\le 0.25\right)$. Middle left: dependence on the characteristic radius. Smooth disks do not depend on the rc. Middle right: dependence on the stellar mass. There is a slight preference to higher values. Bottom left: dependence on the fragmentation velocity. There is a slight preference to higher values.

In the text

Fig. 6 Heat maps of simulations with the R10-0 opacities. From left to right, the three columns represent the smooth case, a planet at 1/3rc, and a planet at 2/3rc. From top to bottom, the rows represent three different snapshots at 300 kyr, 1 Myr, and 3 Myr. The white solid line 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 (reff ~ 10 au) where observational results are available.

In the text

Fig. 7 Heat maps of simulations with the Birnstiel et al. (2018) opacities. From left to right, the three different columns represent the smooth case, a planet at 1 /3rc, and a planet at 2/3rc. From top to bottom, the rows represent three different snapshots at 300 kyr, 1 Myr, and 3 Myr. The white solid line 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 (reff ~ 10 au) where observational results are available.

In the text

Fig. 8 Heat maps of simulations with the Birnstiel et al. (2018) D-50 opacities with 50% porostiy. From left to right, the three columns represent the smooth case, a planet at 1 /3rc, and a planet at 2/3rc. From top to bottom, the rows represent three different snapshots at 300 kyr, 1 Myr, and 3 Myr. The white solid line 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 (reff ~ 10 au) where observational results are available.

In the text

Fig. 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 300kyr – 3 Myr. The black line refers to the disks from the Andrews et al. (2018b) sample. Disks with planets have higher luminosity, while smooth disks have low luminosity at 3 Myr. When two planets are included the luminosity is higher than with a single planet.

In the text

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

In the text

Fig. 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 300 kyr-3 Myr. From left to right, the smooth case, a planet with planet/star mass ratio at 1/3rc, and a planet with the same ratio at 2/3rc. Different limits are used on the x-axis to highlight the evolution of the dust-to-gas ratio. The initial dust-to-gas ratio is 0.01. For the smooth case the dust-to-gas ratio decreases by three orders of magnitude up to 3 Myr. 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 Myr where the ratio is higher than 0.01. The gas mass moves faster than the dust mass in this case.

In the text

Fig. B1 Corner plot of the smooth case. From left to right: Fragmentation velocity υfrag[cm/s], turbulence α-value, disk mass Md[M*], characteristic radius rc[au]. From top to bottom: α, Md[M*], rc[au], M*[M⊙]. In the upper left plot the white dots show the track where the distinction between the fragmentation- and the drift-limited regime is. In the middle plot (rc − α) the white dots show the α-value where the pressure bump cannot hold the dust efficiently anymore. In the panel rc - Md, the mean value of the characteristic radius for each disk mass is shown as a white dot.

In the text

Fig. B.2 Corner plot with a planet with of a planet/star mass ratio 10−3 at a location of the 1/3rc. From left to right: Fragmentation velocity υfrag[cm/s], turbulence a-value, disk mass Md[M*], characteristic radius rc[au]. From top to bottom: α, Md[M*], rc[au], M*[M⊙]. In the upper left plot the white dots show the track where the distinction between the fragmentation and the drift limited regime is. In the middle plot (rc − α) the white dots show the α-value where the pressure bump cannot hold the dust efficiently anymore. In the panel rc- Md, shown as white dots, is the mean value of the characteristic radius for every disk mass.

In the text

Fig. B.3 Corner plot where we have a Jupiter-mass planet at a location of 2/3rc. From left to right: Fragmentation velocity υfrag[cm/s], turbulence α-value, disk mass Md[M*], characteristic radius rc[au]. From top to bottom: α, Md[M*], rc[au], M*[M⊙]. In the upper left plot the white dots show the track where the distinction between the fragmentation- and the drift-limited regime is. In the middle plot (rc − α) the white dots show the α-value where the pressure bump cannot hold the dust efficiently anymore. In the panel rc - Md, shown as white dots, is the mean value of the characteristic radius for every disk mass.

In the text

Fig. C.1 Comparison between a hydrodynamical and our dust evolution simulation for different α-values and for a planet/star mass ratio q = 3 × 10−4 at 30 au. In the left column are the different gap profiles for different α-values and in the right column the corresponding evolution tracks. Even though the profile from Kanagawa et al. (2016) does not overlapping with that obtained from the hydrodynamical simulation, the tracks in the right column produce similar results.

In the text

Fig. C.2 Comparison between a hydrodynamical and our dust evolution simulation for different α-values and for a planet/star mass ratio q = 10−3 at 30 au. In the left column are the different gap profiles for different α-values and in the right column the corresponding evolution tracks. Even though the profile from Kanagawa et al. (2016) does not overlapping with that obtained from the hydrodynamical simulation, the tracks in the right column produce similar results.

In the text

Fig. C.3 Comparison between a hydrodynamical and our dust evolution simulation for different α-values and for a planet/star mass ratio q = 3 × 10−3 at 100 au. In the left column are the different gap profiles for different α-values and in the right column the corresponding evolution tracks. Even though the profile from Kanagawa et al. (2016) does not overlapping with that obtained from the hydrodynamical simulation, the tracks in the right column produce similar results.

In the text

Fig. D.1 Heat maps of representative simulations with the Birnstiel et al. (2018) D-10 opacities with 10% porosity. From left to right, the three columns represent the smooth case, a planet at 1 /3rc, and a planet at 2/3rc. From top to bottom, the rows represent three different snapshots at 300 kyr, 1 Myr, and 3 Myr. The white solid line is the SLR from Andrews et al. (2018b) and the red solid line our fit for the cases where we include a planet. The color bar shows the number of simulations in a single cell. The blue dash-dotted line shows the minimum limit (reff ~ 10 au) where observational results are available.

In the text

Fig. D.2 Heat maps of representative simulations with the Birnstiel et al. (2018) D-90 opacities with 90% porosity. From left to right, the three columns represent the smooth case, a planet at 1 /3rc, and a planet at 2/3rc. From top to bottom, the rows represent three different snapshots at 300 kyr, 1 Myr, and 3 Myr. The white solid line is the SLR from Andrews et al. (2018b) and the red solid line our fit for the cases where we include a planet. The color bar shows the number of simulations in a single cell. The blue dash-dotted line shows the minimum limit (reff ~ 10 au) where observational results are available.

In the text

Fig. D.3 Heat maps of representative simulations with the R10-0 opacity at 3.1 mm for a direct comparison with Tazzari et al. (2021). From left to right, the three columns represent the smooth case, a planet at 1/3rc, and a planet at 2/3rc. From top to bottom, the rows represent three different snapshots at 300 kyr, 1 Myr, and 3 Myr. The white solid line is the SLR from Tazzari et al. (2021). The color bar shows the number of simulations in a single cell. The blue dash-dotted line shows the minimum limit (reff ~ 10 au) where observational results are available.

In the text

	Fig. D.4 Heat map comparison between disks where the stellar luminosity at 3 Myr (right panel) is used instead of 1 Myr (left panel). There is an overall shift toward lower luminosities, but the trend remains the same.
In the text