© J. Stadler et al. 2022

1 Introduction

The discovery of exoplanets shows that planet formation is an efficient process in protoplanetary disks. Observing the circumstellar material around young stars in different star forming regions is crucial to understand the physical processes that regulate the evolution of protoplanetary disks, and hence planet formation.

Surveys at infrared wavelengths with Spitzer and Herschel provided information about the disk fraction around young stars and provided hints of undergoing disk evolution (Oliveira et al. 2010; Fedele et al. 2010; Kim et al. 2016). An example of this evolution are the disks that have low infrared emission compared to the averaged disks, suggesting an inner disk depleted of hot dust (Strom et al. 1989; Calvet et al. 2002; Harvey et al. 2007). Follow-up surveys of protoplanetary disks with state-of-the art telescopes and instruments, such as the Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE) and X-shooter at the Very Large Telescope (VLT) have provided more information about the disk structures of the hot upper layers of the disks (e.g. Avenhaus et al. 2018), the stellar properties and accretion rates (Alcalá et al. 2014, 2017; Manara et al. 2014, 2016, 2017, respectively).

In addition, (sub-)millimetre observations with the Sub-millimetre Array (SMA) and the Atacama Large Millimetre-submillimetre Array (ALMA) have provided data of more than thousands of protoplanetary disks in different star forming regions (e.g. Andrews et al. 2013; Pascucci et al. 2016; Ansdell et al. 2016, 2017, 2020; Barenfeld et al. 2016; Eisner et al. 2018; van Terwisga et al. 2020), giving demographics of the disk properties and the effect of different physical environments (as for example irradiation from massive stars present in some of the observed regions). Some of the most studied properties of these disks are their millimetre luminosities (which are used to infer the dust mass, although with a high level of uncertainty), their dust disk sizes, and their spectral indices.

These observations evidence that the millimetre-flux decreases with age (e.g. Grant et al. 2021), with a couple of exceptions in Corona-Australis and Ophiuchus (Cazzoletti et al. 2019; Williams et al. 2019, respectively). In addition, a correlation between the size of the disk inferred from millimetre observations and the disk luminosity has been found in disks from different regions (Tripathi et al. 2017; Andrews et al. 2018b; Hendler et al. 2020), which seems to flatten in older systems. This correlation may arise from the radial drift of dust particles that is expected to make the disks smaller and less luminous with time, when observed at mm-wavelengths (Rosotti et al. 2019). However, Tazzari et al. (2021a) showed that the outer dust disk radii for a sample of disks in Lupus do not change between wavelengths (0.9, 1.3, and 3.1 mm), which contradicts the expectations of radial drift (see also, Long et al. 2020, Figs. 7 and 8).

Moreover, most of the high resolution observations of protoplanetary disks show that disks are not smooth, but instead they display substructures (e.g. Andrews et al. 2018a; Long et al. 2018; Cieza et al. 2021). These substructures may be the product of dust trapping in pressure maxima (e.g. Pinilla et al. 2012; Dullemond et al. 2018) which can reduce or fully suppress the radial drift. Hence, if most of the disks host pressure bumps it is natural that the observed correlations contradict current models where radial drift is very efficient. However, a recent study by Jiang et al. (2022) also questions this link between continuum rings and pressure bumps, as they find no statistically significant evidence for a universal correlation between line emission and continuum substructures.

Pressure maxima can have different origins in the disk and they may have different behaviours and leave different imprints in the disk depending on the origin (Pinilla & Youdin 2017, for a review). These origins can be divided in planetary origin, (magneto-) hydro-dynamical (Flock et al. 2015), or dust evolution processes (e.g. Okuzumi & Inutsuka 2015, for variations of the dust properties and evolution near the ice lines). Currently, it is debatable if most of the observed structures are due to planets, partially because planets have been only confirmed in PDS 70 (Keppler et al. 2018; Müller et al. 2018) and very recently suggested in AB Aur (Currie et al. 2022; Zhou et al. 2022). One of the main differences that are expected from planetary pressure bumps versus (magneto-) hydro-dynamical effects is how long these bumps live and how dynamic they are. On one hand, once a planet is formed in the disk and it is massive enough to perturb the pressure profile of the disk or even the gas azimuth velocities to create traffic jams, the expected structures are quite stable in time (they do not appear and disappear). Still, the bump might jointly migrate (inwards) with the planet (e.g. Meru et al. 2019; Pérez et al. 2019; Weber et al. 2019), as long as it does not undergo rapid runaway migration, for which multiple dust traps can be formed (e.g. Wafflard-Fernandez & Baruteau 2020). On the other hand, pressure bumps that are from (magneto-)hydro-dynamical origin, also known as zonal flows, may appear and disappear in different disk locations (e.g. Uribe et al. 2011; Simon et al. 2013; Flock et al. 2015) with different properties (e.g. amplitudes, widths).

In this paper, we aim to understand if current observations enable us to distinguish between planetary- versus zonal flows-origin of the pressure bumps. We investigate the effect that different pressure bumps properties have on observational quantities such as the disk luminosity (including the estimated dust mass), the spectral index, the disk dust outer radius, and the optical depth. So far, theoretical comparison of observed relations, such as the correlation between the disk dust size the disk luminosity has mostly been studied in detail in the absence of pressure bumps (Powell et al. 2019; Toci et al. 2021). Only recently, Zormpas et al. (2022) find that the inclusion of strong pressure bumps (planets) flatten the slope of the size-luminosity relation and that the location of the outermost planet defines the dust size of the disk. In this work, we focus our analysis of dust evolution on different properties of the pressure bumps, including: their amplitude and location, the number of bumps, and the growth timescale of the bumps. Additionally, we also perform simulations with variable pressure bumps, in order to mimic zonal-flow type bumps.

This paper is structured as follows: In the upcoming section, the models for the disk and dust evolution and pressure bumps are described. The third section introduces the models’ simulation setups with the disk’s initial conditions, variable pressure bump parameters and the assumptions for the radiative transfer calculations. In Sect. 4, the results are presented. We start with an extensive parameter study of the pressure bump variables, where in each subsection we present the parameter’s influence on the disks’ observables. The following subsection will focus on the results of the dynamic bumps model, which tries to mimic zonal flows in disks. In Sect. 5, we discuss our results in the context of the various disk observables, such as the dust disk radii, millimetre fluxes, spectral indices, and dust masses. Additionally, we investigate the influence of the dust opacity prescription and planetesimal formation on our results. Finally, in Sect. 6 we draw our main conclusions.

2 Models

2.1 Disk evolution

In our model, we consider a 1D axi-symmetric circumstellar disk whose viscous evolution can be described by the following continuity equation (Lüst 1952; Lynden-Bell & Pringle 1974; Pringle 1981):

$$\frac{\partial }{\partial t}{\text{Σg}}_{}=-\frac{1}{r}\frac{\partial }{\partial r}\left({\text{Σg}}_{}r{v}_{\text{g}}\right),$$(1)

where Σ_g is the gas surface density, r is the radial distance to the star and v_g is the gas radial velocity due to the viscosity which is given by (Lynden-Bell & Pringle 1974)

$$v_{\text{g}}=-\frac{3}{{\text{Σg}}_{}\sqrt{r}}\frac{\partial \left({\text{Σg}}_{}v\sqrt{r}\right)}{\partial r}.$$(2)

For the kinematic viscosity v, we use the so called α-prescription for viscous accretion disks from Shakura & Sunyaev (1973) that states

$$v=\alpha c_{\text{s}}{h}_{\text{g}}.$$(3)

Here, α describes a dimensionless quantity that measures the efficiency of the angular momentum transport in the disk, h_g = c_s/Ω_k is the gas pressure scale height with ${\text{Ωk}}_{}=\sqrt{G{M}_{*}\text{/}{r}^3}$ being the Kepler frequency. The isothermal sound speed can be written as $c_{\text{s}}=\sqrt{k_{\text{B}}T\text{/}\mu m_{\text{p}}}$, where k_B is the Boltzmann constant, T the disk’s temperature, m_p the proton mass and µ = 2.3 the mean molecular weight.

2.2 Pressure bump model

The focus of this work is on the properties of pressure bumps and how these affect observational properties of disks. Hence, we want to create bumps in the pressure profile of the gas, where the disk pressure in our vertically isothermal assumption is given by

$$P={\rho }_{\text{g}}{c}_{\text{s}}^2=\frac{{\text{Σg}}_{}{c}_{\text{s}}^2}{\sqrt{2\pi }{h}_{\text{g}}}$$(4)

with the gas midplane density ${\rho }_{\text{g}}={\Sigma }_{\text{g}}\text{/}\sqrt{2\pi }{h}_{\text{g}}$, where most of the mass can be expected.

Physically, a planetary gap forms when a planet is massive enough for its wakes to create shocks that progressively expel gas out of its horseshoe region. In our model, we implement this by locally increasing the turbulent viscosity, hence the gas radial velocity v_g has a strong increase in the gap, so the gap remains even if the disk is viscously evolved. The gap in the gas surface density then naturally produces a local pressure maximum or bump at the outer gap edge, since P(r) declines with increasing r. This prescription has similarly been used in several previous studies (e.g. Dullemond et al. 2018; Stammler et al. 2019; Pinilla et al. 2021).

To this end, we modify the turbulent viscosity parameter by imposing a Gaussian factor on its constant background value α₀ as

$$\alpha \left(t,r\right)={\alpha }_0·\left(1+{\displaystyle \sum _i{F}_i\left(t,r\right)}\right),$$(5)

where F_i (t, r) for each bump i is given by

$$F_i\left(t,r\right)=A\left(t\right)·exp\left(\frac{-{\left(r-r_{\text{gap,}i}\right)}^2}{2w_{\text{gap,}i}^2}\right),$$(6)

with r_gap the radial gap location, w_gap the Gaussian width and A the time dependent amplitude of the bump. We parameterize the latter as

$$A\left(t\right)=\{\begin{array}{ll}0\hfill & t<T_{\text{ini}}\hfill \\ A_{\max }\hfill & t\ge T_{\text{fin}}\hfill \\ A_{\max }·\sin \left(\frac{\pi }{2}\left(\frac{t-T_{\text{ini}}}{T_{\text{fin}}-T_{\text{ini}}}\right)\right)\hfill & \text{otherwise,}\hfill \end{array}$$(7)

with A_max the maximum amplitude, T_ini the initial time where the bump starts to appear in the disk and T_fin the final time of the growth phase, after which the bump reached its full amplitude A_max. In contrast to these long-lived ‘planetary’ bumps, we will also explore models, in which the bump will disappear again. For those, Eq. (7) will be adjusted by A_max = 0 for t ≥ T_fin and the sin-function is normalised by π instead of π/2.

We note that the above prescribed bumps have fixed radial locations throughout the simulation, therewith neglecting any planetary migration, a simplification which will be discussed in Sect. 5. Furthermore, the bumpy α-profile gets only imposed on the kinematic viscosity, while keeping the global constant value of α₀ for the turbulent velocities of the dust particles, as well as for the radial and vertical diffusion of the dust, which will be introduced in the upcoming sections. Physically, this implies that what sets the gas viscous evolution may be detached from the above mentioned processes of the dust evolution, which has been recently studied in detail by Pinilla et al. (2021).

Dynamic bumps A special case of the above described disappearing bumps is our Dynamic Bumps model, which mimics zonal flows. First, we divide the disk up into several zones. The zones centres are located in the disk where r/(5 × h_g) are integers and their width is such that, this ratio spans ±0.5 around the integers, for example the boundaries of the first zone are at 2.5 and 7.5 h_g, respectively. Within each of these zones, pressure bumps can appear and disappear at random positions. However, only one bump can form within each zone and the bumps need to have at least a distance of one pressure scale height h_g to each other. As in the above ‘planetary’ model, the gaps in the gas surface density get created through a Gaussian bump, imposed only on the turbulent α-viscosity profile (Eq. (6)), affecting only the gas kinematics but not those of the dust.

Before the bumps actually form and get imposed on the α-profile a certain formation time T_form needs to pass. This is similar to MHD-simulations of the Rossby wave instability (RWI) that creates vortices, where it also takes some time till over-densities accumulate in the disk, which are then able to trap dust (e.g. Meheut et al. 2012; Flock et al. 2015; Cui & Bai 2021). After the formation time passed, the bumps have a given lifetime T_life, in which their amplitude evolves, like the sinusoidal function defined in Eq. (7), and after T_life elapsed they disappear again. We set both T_form and T_life proportional to the local orbit time t_orb = 2π/Ω_k at the corresponding radial location of the bump. Once a bump has lived for its given lifetime, it disappears and a new random position within the zone will be drawn for a new bump and the cycle begins anew.

It is important to state, that this model does not create bumps in the gas pressure based on physical conditions, but is rather forced to act in a similar way as seen in numerical simulations (e.g. Dittrich et al. 2013; Lenz et al. 2019). As zonal flows are strictly speaking over-densities in the gas density, we also ran a comparison model inducing gaps rather than bumps in α, leading to bumps in Σ_g. However, we found no major differences between the two prescriptions, so we choose to implement the pressure bumps as bumps in α similar to the ‘planetary’ bump model.

2.3 Dust evolution

The evolution of the dust is of particular importance, since its distribution determines the disk’s appearance in the mm-continuum that are observed with interferometers like ALMA and which we want to model in this work. It’s important to emphasise that the dust dynamics differs from that of the gas, as dust particles drift towards the star or concentrate at the mid-plane, which will be explained in the following paragraphs. This varying behaviour of dust and gas was predicted decades ago by theory (Whipple 1972; Weidenschilling 1977) and has been repeatedly seen in observations (e.g. HD 163296, Isella et al. 2018; Varga et al. 2021; Öberg et al. 2021; Teague et al. 2021).

The most important quantity to track the dust dynamical behaviour is the Stokes number St = τ_stopΩ_k, where the stopping time τ_stop describes the timescale on which dust particles couple to the gas motion. The strength of the (de)coupling, in turn, depends on the difference between gas and dust velocities (v_g and v_d, respectively), as well as on the particle’s Stokes number through the drag force, which gets exerted onto the dust by the gas (per unit mass) as:

$$f_{\text{D}}=-\frac{v_{\text{d}}-v_{\text{g}}}{{\tau }_{\text{stop}}}.$$(8)

As the global pressure maximum lies towards the star, the gas usually flows at sub-Keplerian velocity due to its own pressure support. Small enough dust particles (St ≪ 1) are well coupled to the gas, meaning they are dragged along with it. Larger boulders (St ≫ 1), however, can decouple from the gas motion and are not affected by it. A special case are mid-sized pebbles (St ≈ 1), which centrifugal force that acts against gravity is most insufficient, causing them to drift (inwards) the fastest.

Assuming spherical dust grains with radii a and a material density ρ_s, the two prevalent regimes for the Stokes number at the disk’s midplane are (Brauer et al. 2008; Birnstiel et al. 2010):

$$\text{St}=\frac{\pi }{2}\frac{a{\rho }_{\text{s}}}{{\text{Σg}}_{}}\{\begin{array}{lll}1\hfill & a<9/4·{\lambda }_{\text{mfp}}\hfill & \text{Epstein\hspace{0.17em}regime}\hfill \\ \frac{4}{9}\frac{a}{{\lambda }_{\text{mfp}}}\hfill & \text{otherwise}\hfill & \text{Stokes\hspace{0.17em}regime\hspace{0.17em}1,}\hfill \end{array}$$(9)

where the mean-free path of the gas is given by ${\lambda }_{\text{mfp}}=1/\left(\sqrt{2}n{\sigma }_{{\text{H}}_2}\right)$, with the gas’ number density n = ρ_g/µ m_p and the cross-section for molecular hydrogen ${\sigma }_{{\text{H}}_2}=2\times {10}^{-15}$ cm². Based on the above equation, we can note that the Stokes number is essentially a ‘dynamical grain size’ that tracks the dust grain dynamics. We further define the dust’s radial velocity following Nakagawa et al. (1986) and Takeuchi & Lin (2002):

$$v_{\text{d}}=\frac{1}{1+{\text{St}}^2}{v}_{\text{g}}-\frac{2\text{St}}{1+{\text{St}}^2}\eta v_{\text{k}},$$(10)

with η = − (1/2)(h_g/r)² dlnP/dln r and the Keplerian orbital velocity v_k.

Additionally dust particles also diffuse due to the gas’ turbulent motion, which counter acts concentration gradients. The dust diffusivity is hereby given through D_d = v/(1 + St²) (Youdin & Lithwick 2007). The overall dust evolution can then be described by an advection–diffusion equation following Birnstiel et al. (2010):

$$\frac{\partial {\text{Σd}}_{}}{\partial t}+\frac{1}{r}\frac{\partial \left(r{\text{Σd}}_{}{v}_{\text{d}}\right)}{\partial r}-\frac{1}{r}\frac{\partial }{\partial r}\left(r{D}_d{\text{Σg}}_{}\frac{\partial }{\partial r}\left(\frac{{\text{Σd}}_{}}{{\text{Σg}}_{}}\right)\right)=0,$$(11)

which is numerically being solved for the surface density Σ_d of each dust species, i.e. size bin.

We note that we do not consider the effects of the dust back-reaction onto the gas in this work, which can safely be done, since this effect becomes negligible for low dust-to-gas ratios (є ≤ 0.01) and long timescales (~Myr), as shown in Gárate et al. (2020).

2.3.1 Dust vertical settling

The vertical motion of dust particles is important when considering collisions of dust particles and, therewith associated, the growth of grains to larger sizes (which will be explained in more detail in the next subsection). The settling velocity of grains well coupled to the gas (St < 1) is given by v_Z,Settl = −z Ω_kSt (Safronov 1969; Dullemond & Dominik 2004). Thus, larger particles with high Stokes numbers (for St<1) quickly settle towards the midplane, where v_Z,Settl then decreases. Dust settling in protoplanetary disks is nowadays confirmed by high resolution observations of highly inclined disks (e.g. Villenave et al. 2020).

The height above the midplane which influences the relative velocities between grains and thus their growth, is approximated by the dust scale height $h_{\text{d}}=h_{\text{g}}\sqrt{{\alpha }_0\text{/}\left({\alpha }_0+\text{St}\right)}$ (Dubrulle et al. 1995). While this work does not treat the vertical evolution of the dust density distribution, we do compute the 3D total volume density of the dust ρ_d, which will later serve as an input for the radiative transfer calculations. Following Pinilla et al. (2021), we define it such that

$${\rho }_{\text{d}}\left(r,\phi ,z,\text{St}\right)=\frac{{\text{Σd}}_{}}{\sqrt{2\pi }{h}_{\text{d}}}\exp \left(-\frac{z^2}{2h_{\text{d}}^2}\right),$$(12)

also using the prescription of Birnstiel et al. (2010) for the dust scale height:

$$h_{\text{d}}\left(\text{St}\right)=h_{\text{g}}·\min \left(1,\sqrt{\frac{{\alpha }_0}{\min \left(\text{St},1/2\right)\left(1+{\text{St}}^2\right)}}\right),$$(13)

where it is ensured that h_d does not become larger than the gas scale height h_g, which would otherwise be unphysical. We note that we found no differences in the disk’s observables of our Fiducial setup, whether using the Birnstiel et al. (2010) or Dubrulle et al. (1995) prescription of h_d.

2.3.2 Dust growth

In this work, we simultaneously describe the evolution of dust and gas, including the growth and fragmentation of the dust grains (Brauer et al. 2008; Birnstiel et al. 2010).

In general, the dust phase consists of a distribution of different particle sizes, whose dynamics are mainly determined by their Stokes number (Eq. (9)). The particles interact and collide with each other, where different collisional outcomes are possible, such as sticking, fragmentation or erosion of grains. The Smoluchowski coagulation equation (Smoluchowski 1916) is solved to model the dust growth. Whether particles grow or fragment in this process mainly depends on their respective sizes or masses and their relative velocities with respect to each other, as described in Birnstiel et al. (2010).

While the velocities of small particles (≤ 1 µm) are governed by Brownian motion, the picture changes for larger grains (see Fig. 1 of Birnstiel et al. 2016). Then, there are different transport mechanisms that become the dominant contributors of the particles’ relative velocities before collision, namely turbulent mixing, vertical stirring or settling and azimuthal or radial drift. If grains grow, fragment, or erode upon collision is thus determined by comparing the computed relative velocities with the fragmentation speed v_frag, a parameter that we take to be constant throughout the disk (Birnstiel et al. 2009). In the fragmentation limit, the maximum Stokes number or size a particle can reach is

$${\text{St}}_{\text{frag}}=\frac{1}{2}\frac{v_{\text{frag}}^2}{{\alpha }_0{c}_{\text{s}}^2}.$$(14)

The other process that locally limits the growth of the particle distribution is radial drift (Brauer et al. 2008; Birnstiel et al. 2010). It dominates when the growth timescale exceeds the drift timescale, so particles drift faster towards the star than they can grow to larger sizes. In this case, the maximum Stokes number is given by

$${\text{St}}_{\text{drift}}={\left|\frac{\text{d}\ln P}{\text{d}\ln r}\right|}^{-1}\frac{v_k^2}{c_{\text{s}}^2}\epsilon ,$$(15)

with є = Σ_d/Σ_g being the column density dust-to-gas ratio. The above equations then ultimately set the maximum size a particle can grow to St_max = min(St_frag, St_drift) (Birnstiel et al. 2012), where we neglect dust settling which can also set the maximum grain size (Pinilla et al. 2021). A display of the both size limits (orange and cyan lines) can be seen in Fig. 1 that shows the dust size distribution for our Fiducial model at 1 Myr of evolution, whose parameters can be found in Table 1. Inside the pressure bumps the maximum grain size is limited by fragmentation because radial drift is reduced at these locations. In the rest of the disk, the maximum grain size is usually limited by drift at Myr timescales. The white solid line in the plot corresponds to the particle size of St = 1, which is proportional to the gas surface density (see Eq. (9)).

Fig. 1 Dust size distribution for the Fiducial model at 1 Myr. The contours show the growth limits for drift (orange) and fragmentation (cyan), as well as the particle sizes corresponding to St = 1 (white) and 1 mm (grey).

3 Simulation setup

3.1 Disk setup

We carried out our disk simulations with the code DustPy¹ by Stammler & Birnstiel (2022) which builds on the dust evolution model of Birnstiel et al. (2010). This code solves both the gas and dust evolution equations described in the previous section. We assume the disk to be around a Solar mass star and initially set up the gas surface density of the disk following the profile of Lynden-Bell & Pringle (1974):

$${\text{Σ}}_g\left(r,t_0\right)={\text{Σ}}_0{\left(\frac{r}{R_{\text{c}}}\right)}^{-\gamma }\exp \left[-{\left(\frac{r}{R_{\text{c}}}\right)}^{2-\gamma }\right],$$(16)

with γ = 1 and the normalisation factor ${\text{Σ}}_0=M_{\text{disk}}/\left(2\pi R_{\text{c}}^2\right)$ is defined through the initial disk mass, which we set to be M_disk = 0.05M_⊙. We take the characteristic radius R_c = 60 au, from where the profile gets exponentially tapered. The radial grid contains n_r = 300 logarithmically spaced cells and spans from 5 to 1000 au, as we are only interested in the outer parts of the disk.

The midplane temperature profile is described by a power law of a flaring passively irradiated disk with

$$T\left(r\right)=T_{\text{eff}}{\left(\frac{0.5{\phi }_{\text{irr}}{R}_{*}^2}{r^2}\right)}^{1/4}\approx 160\text{K}{\left(\frac{r}{1\text{\hspace{0.17em}au}}\right)}^{-1/2},$$(17)

where the stellar effective temperature T_eff = 3500 K, the stellar radius R_* = 2R_⊙ and the irradiation angle φ_irr = 0.05. We choose this lower than solar value of T_eff, to shrink the width of our imposed gaps, which are determined by the pressure scale height and therewith the disk’s temperature (see next section). Additionally, we assume the dust and gas to have the same temperature and the disk to be vertically isothermal.

For the mass distribution of the dust particles, we also take a logarithmically spaced grid with n_m = 141 cells, covering sizes from 0.5 µm up to 2 m. The initial dust-to-gas ratio is set to be є₀ = 0.01 and initially the dust particles will have a size distribution n(a) ∝ a^−3.5, following the so-called MRN distribution of interstellar grains by Mathis et al. (1977) with a maximum grain-size of a₀ = 1 µm. We assume the grains to be compact spheres with a mass given by m = 4πα³ρ_s/3 and their material density ρ_s = 1.675 g cm⁻³, in line with our opacity composition. Finally, we assume the fragmentation velocity to be v_frag = 10 m s⁻¹ throughout the whole disk (Gundlach et al. 2011; Gundlach & Blum 2015). An overview of the disk model parameters are listed in Table 2.

3.2 Pressure bump variables

To study the nature of pressure bumps and their observables, we want to explore a variety of possible parameters in our pressure bump models. The key parameters to be varied are the gap location r_gap and the time-dependent amplitude of the Gaussian bump A(t). The latter is defined through its maximum amplitude A_max and the growth time T_grow = T_fin – T_ini (see Eq. (7)). To this end, we vary each of these parameters with respect to a Fiducial model, while the other parameters are kept constant. The parameters of the Fiducial model, which later will be taken as a reference point, are r_gap = 30 au, A_max = 2 and T_fin = 0.10 Myr. An overview of the parameter space for the ‘planetary’ pressure bump models can be found in Table 1.

Initially, we include the bumps at T_ini = 1000 yr and then let them grow over several thousand years. However, how fast bumps can form is locally governed by the viscous timescale (Armitage 2010). Inserting ∆r = w_gap and simplifying, we obtain:

$${\tau }_{\text{visc}}\left(r\right)\approx \frac{{\left(\text{Δ}r\right)}^2}{v}=\frac{w_{\text{gap}}^2}{v}=\frac{4}{{\alpha }_0{\text{Ωk}}_{}},$$(18)

since we induce them in a. For the Gaussian width of the bump, in the above equation, we take w_gap = 2 h_g(r_gap) at the gap location and also inserted the kinematic viscosity. The width needs to be larger (or equal) than hg to ensure the stability of the bumps (Li et al. 2000; Ono et al. 2016). This timescale is inversely proportional to the Keplerian frequency, so outer disk regions have longer viscous time scales, and accordingly bumps in the outer parts will take longer to grow, which naturally sets a lower limit of T_grow ≥ τ_visc(r_gaP).

In order to actually create pressure bumps, where dust particles can be trapped, it is essential that there is a radius where the pressure gradient changes sign, meaning there exists a local pressure maximum in the disk (Pinilla et al. 2012). For radially increasing gap locations, the pressure scale height, thus the Gaussian width w_gap and the disk’s aspect ratio (h_g/r ∝ r^1/4) increase, leading to a smaller Gaussian factor in Eq. (5) and a shallower pressure gradient. This puts a constraint on the amplitude of the bump A_max which needs to be large enough to fulfil the above condition. To this end, we take a higher value A_max = 3 at the gap location of r_gap = 80 au, to ensure to create a particle trap.

3.3 Dynamic bumps model

For the Dynamic Bumps model, we do not vary the parameters for certain bumps in terms of their amplitude and location, but rather want to focus on which effect the bumps’ lifetimes have on the disk’s observational properties. As described in Sect. 2.2 the disk will be divided up into several zones. Since we are only interested in the outer parts of the disk, we restrict the regions where bumps can form to be between 8–100 au. We set the Gaussian bumps width to be w_gap = 2h_g, in line with the radial extent of vortices found in global MHD-simulations (e.g. Flock et al. 2015; Cui & Bai 2021). As explained in the previous section, we take A_max = 3, as to create efficient dust traps across the whole disk. However, this leads to the fact that dust trapping will be more efficient in the inner zones, where the pressure gradient will be steeper due to the flared nature of our disk (h_g/r ∝ r^1/4). We assume that all bumps start forming at T_form = 500 t_orb local orbit times. We then run three simulations varying the lifetimes of the dynamic bumps by T_life = (200, 400, 1000) t_orb.

We want to compare these dynamic bumps models to a static (‘planetary’) one. To accurately do that, we set up another simulation where we induce three bumps at r_gap = [10, 30, 80] au into the disk. They take T_form = 500 t_orb to appear and then another T_grow = 200 t_orb to reach their final amplitude of A_max = 3, after which they remain unchanged in the disk.

3.4 Radiative transfer calculations

To infer the observational appearance of the disk, we perform radiative transfer calculations with the code RADMC-3D² (Dullemond et al. 2012). As input, we take the 3D dust density distribution (see Eq. (12)) for each grain size obtained from our DustPy simulations at an evolution time of 1 Myr.

We assume the composition of the grains as stated in the Disk Substructures at High Angular Resolution Project (DSHARP) paper of Birnstiel et al. (2018), consisting of water ice (Warren & Brandt 2008), astronomical silicates (Draine 2006), as well as troilite and refractory organics (Henning & Stognienko 1996).

We then calculate the opacities for each grain size using the program Optool³ (Dominik et al. 2021), where we further assume the grains to be compact spheres. The scattering phase function for grains that are comparable in size to the modelled wavelength has a strong forward peak. This poses a problem, due to the limited resolution of our radiative transfer model. Therefore, we set the scattering opacity of the grains to an upper limit within the first 5° of the phase function, which thereafter needs to be renormalised again. This means that within those first degrees forward scattering will effectively be treated as no scattering at all.

In accordance with our disk model, we set up a Solar mass star with an effective temperature of T_eff = 3500 K which is treated as a blackbody point source and emits 10⁷ photon packages for our radiative transfer calculations. The spherical grid of the model has N_r = 300 radial, N_ϕ = 64 azimuthal and N_θ = 180 polar cells with the star at its centre. The latter corresponds to three or six cells per pressure scale height at 10 and 100 au, respectively. The wavelength grid, which is needed to compute the dust temperature, spans from 0.1 µm to 1 cm and contains 150 cells.

First of all, the 3D dust equilibrium temperature gets computed, which is then used to calculate synthetic observations of the disk in the 1.3 mm and 3.0 mm dust continuum. Hereby, we use the full treatment of anisotropic-scattering of dust grains with polarisation (Kataoka et al. 2015) of the code. Even though, it is computationally expensive, we later want to compare this to models where we do not treat scattering and investigate the differences.

For all calculations, we assume the disk to have an inclination of i = 20° and to be at a distance of d = 140 pc, which is the mean value of the nearest star forming regions. In a final step, the synthetic images get convolved by a Gaussian point-spread-function with a full-width at half-maximum (FWHM) of 0.04″ × 0.04″ to mimic the current angular resolutions obtained by interferometers like ALMA (e.g. Andrews et al. 2018a) and to make sure to resolve the disk’s substructures.

From the convolved intensity profiles of the disks, we can calculate their millimetre fluxes at the two wavelengths, as well as their dust disk radii. The latter is defined through the radius where 68% or 90% of the emission at 1.3 mm (or additionally 3.0 mm) is enclosed, by assuming that the sensitivity of the observations has a peak signal-to-noise ratio (S/N), with respect to the maximum of the emission, of 1000. We choose this peak S/N, instead of a constant value, following Pinilla et al. (2021). The reason behind this is that each observation has a different sensitivity, depending on the total flux of the disk, as well as on the limitations of the integration time. Current observations conducted with ALMA are able to achieve this S/N for the brightest disks, however for fainter disks this might not be achievable without too long integration times.

Furthermore, we calculate the disk’s radially integrated millimetre spectral index taking its fluxes at 1.3 mm and 3.0 mm as

$${\alpha }_{\text{mm}}=\frac{{\log }_{10}\left(F_{1.3\text{\hspace{0.17em}mm}}\text{/}{F}_{3.0\text{\hspace{0.17em}mm}}\right)}{{\log }_{10}\left(3.0\text{/}1.3\right)}.$$(19)

Another quantity, which we want to infer from our synthetic observations is the dust mass, that can be calculated through the optically thin approximation of Hildebrand (1983):

$$M_{\text{d},ob\text{s}}\simeq \frac{d^2{F}_{1.3\text{\hspace{0.17em}mm}}}{{\kappa }_v{B}_v\left(T\left(r\right)\right)},$$(20)

where d is the distance to the source, F_1,3mm is the flux taken at λ = 1.3 mm and B_v is the Planck function, for which mostly a blackbody temperature of T = 20 K is assumed. For the opacity, we take a frequency-dependent prescription, that is frequently used in disk surveys and is given by κ_v = 2.3 cm² g⁻¹ × (v/230 GHz)^0.4 (Andrews et al. 2013). It is important to stress that this relation only holds as long as most of the disk has an optical depth of τ ≤ 1, otherwise the optically thin approximation breaks down and the dust mass cannot accurately be traced anymore.

Lastly, we actually want to infer the optical depth from our synthetic observations and compare this to the absorption optical depth of our models:

$${\tau }_{v,\text{abs}}\left(r\right)={\displaystyle \sum _i^N{\text{Σd}}_i^{}\left(r\right){\kappa }_v^{i,\text{abs}},}$$(21)

with the absorption opacity ${\kappa }_v^{i,abs}$ calculated from the DSHARP opacities (Birnstiel et al. 2018) and where the sum index corresponds to the different dust species. The highest emission of our bump models is naturally measured at the peak location of dust ring r_peak, as can be seen in Fig. 6. Therefore, we also want to focus on computing the optical depth at this peak location, since there we might expect the emission to become optically thick. Observationally, this is usually calculated by solving the radiative transfer equation for the optical depth (see Eq. (11) of Dullemond et al. 2018)

$$\begin{array}{c}{I}_v\left(r\right)=B_v\left(T_{\text{d}}\left(r\right)\right)\left(1-\exp \left[-{\tau }_v\left(r\right)\right]\right)\\ \iff {\tau }_{\text{obs}}=-\ln \left(1-\frac{I_{1.3\text{\hspace{0.17em}mm}}\left(r_{\text{peak}}\right)}{B_{1.3\text{\hspace{0.17em}mm}}\left(T_{\text{d}}\left(r_{\text{peak}}\right)\right)}\right),\end{array}$$(22)

where both the intensity I_v and the temperature of the blackbody function B_v get evaluated for λ = 1.3 mm. For the temperature T_d, we use the profile introduced in Eq. (17).

4 Results

4.1 Parameter effects on observables

In this subsection, we will present the parameter study of our ‘planetary’ pressure bump model and how the three main variables, namely the gap location r_gap, the maximum amplitude A_max and the growth time of the bump T_grow, affect the observational properties of our disk. To this end, we vary each of the parameters with respect to our Fiducial model, while the other parameters are kept constant. An overview of the different models with their corresponding parameters and observational properties can be found in Table 3.

4.1.1 Dust disk radii

The main result of this subsection is that the (farthest out) pressure bump determines the disk dust radius. We show this, investigating how the dust disk radius changes with our different bump models. In Fig. 2, we plotted dust disk radius enclosing 68% and 90% of the cumulative flux at 1.3 mm over each model parameter. An overview plot of all model’s millimetre intensity profiles with the location of each emission radius can be found in Fig. A.4.

In plot 2a, we find that the disk radius shows a linear relationship to the gap location of our model. Hence, the (farthest out) pressure bump determines the disk dust radius. Interestingly, the R_d,68% and R_d,90% (green and blue crosses) of the Two Bumps model, plotted in the figure, each trace its corresponding bumps at 30 and 80 au, while the rest of the models have no big difference between their 68% and 90% radii. All the data points lie above the 1:1 line, since the cumulative flux traces the emission from the dust rings, which lie exterior to the corresponding gap location.

The amplitude of the bump has only a minor effect on the disk radius, as long as a particle trap is created (A_max ≳ 2). This can be seen in (b), where higher amplitudes create deeper gaps that move the local pressure maximum further out, therewith increasing the disk radius by a small factor.

In subplot (c), we see that the growth time of the bump has no effect on the dust disk radius at all, which just stays constant with varying T_grow.

4.1.2 Millimetre fluxes

Small gap locations (r_gap < R_c) and short growth times of the bumps are needed to observe high millimetre fluxes. This can be seen in Fig. 3, where we study the fluxes at λ = 1.3 mm and λ = 3.0 mm dependent on the bump parameters.

In the first plot (a), we see that the Fiducial model, with the gap at 30 au, displays the highest millimetre flux. Although, the model with the gap further in traps more dust, its millimetre flux is lower since the width of the dust ring is narrower (w_gap = 2 h_g(r_gap)) which leads also to stronger dust trapping around the pressure maximum. This increases the optical depth of the Inner Bump model, thereby shielding more dust from our view and reducing its emission. We will explain the influence of the optical depth on our results in more detail in a coming subsection.

For the Outer Bump model, with its gap further out, the flux is lower as it traps less dust, because r_gap > R_c and most of the dust mass is located within the characteristic radius. Hence, the bump cannot trap enough dust mass from drifting inwards, in order to grow sufficient millimetre pebbles for high millimetre fluxes. This also explains, why the F_{3.0 mm} drops more than the F_{1.3 mm}, as longer wavelengths trace larger particles.

A last thing to point out in Fig. 3a is that the fluxes for Two Bumps model (cross markers) very well coincide with the Fiducial one, even though it has a second bump at 80 au. The model’s additional bump only increases F_{1.3 mm} by 20%, while F_{3.0 mm} decreases compared to the Fiducial setup.

In subplot (b), we investigate the influence of the bump’s amplitude on the millimetre fluxes. Here, the only major change in fluxes happens when the amplitude A_max increases from 1 to 2, because an amplitude factor of 1 is insufficient to actually trap particles. This bump only slows down the particle drift and is not able to halt a large amount of dust mass from accreting onto the star (see also Fig. 4c). The minor decrease in the fluxes for A_max > 2 is due to the fact, that for larger bump amplitudes, we increase the pressure gradient. Thus, the grains accumulate closer around the local pressure maximum, leading again to higher optical depths.

The third subplot (c) displays the flux dependency on the growth time of the pressure bumps. We find the general trend that both millimetre fluxes decrease for longer growth times. This is straightforward to understand, since the earlier the pressure bump is in the disk, the more dust it can trap from drifting inwards and also the more time dust particles have to grow to larger mm-sizes. We further note, that F_{3.0 mm} decreases faster for larger growth times than F_{1.3 mm} which will have a direct effect on the spectral indices, as we will discuss in the following subsection.

Fig. 2 Dust disk radius parameter plots. Dust disk radii enclosing 68% (green) or 90% (blue) of the 1.3 mm cumulative flux plotted over the gap location, bump amplitude and bump growth time (from a-c). In sub-plot (a), the cross markers correspond to the Two Bumps model and the dashed-dotted line shows the 1:1 -line.

Fig. 3 Flux and spectral index parameter plots. (Upper) Disk millimetre fluxes at 1.3 mm (blue) and 3.0 mm (green) and according millimetre spectral indices (lower) plotted over the gap location, bump amplitude and bump growth time (from a-c).

4.1.3 Spectral indices

Strong and fast-forming dust traps in the disk interior are needed to produce small spectral indices of disks. They are determined through the logarithmic ratio of the mm-fluxes of our disk models and we show their dependence on the bump parameters in Fig. 3.

The lower left subplot shows, that the spectral index increases for bumps located further out in the disk, which is a consequence of less mm-dust-grains prevalent in the bumps, leading to lower 3.0mm-fluxes, as well as to lower optical depths at the dust ring location. Curiously, the Two Bumps model has a higher spectral index than the Fiducial model, due to its higher F_{1.3 mm} (see Fig. 3a). This happens as some of the dust grains that were trapped in the inner, optically thick regions in the Fiducial model, are now retained in the optically thin outer regions.

In plot (b), we see the spectral index’ dependency on the bump amplitude which decreases for higher values of A_max. Since the Smaller A model does not trap particles sufficiently, it has a quite high spectral index. For larger amplitudes, α_mm then only decreases a little which is due to the stronger trapping around the pressure maxima, as explained in the earlier section.

The lower plot (c) in Fig. 3 shows that the spectral index increases for longer growth times of the bump. As already mentioned in the last subsection and which can be seen in plot (c) above, F_3.0mm decreases faster for long growth times than F_1.3mm. Thus, for higher T_grow less dust mass is trapped to grow sufficiently large mm-sized pebbles with high opacities that particularly contribute to F_{3.0 mm}.

As already mentioned in the last subsection, it understandably needs longer times and higher accumulations of dust masses to grow sufficiently large mm-sized pebbles with high opacities that especially contribute to a higher F_{3.0 mm} compared to F_{1.3 mm}, whose ratio in the end determines α_mm.

Fig. 4 Dust mass and optical depth parameter plots. (Upper) Total dust mass (blue), dust mass for grains larger than 0.1 mm (orange) and observed dust mass (green) plotted over the gap location, bump amplitude and bump growth time (from a-c). (Lower) Absorption and observed peak optical depth also in dependence on bump model parameters. The Inner Bump model reaches τabs = 35, which is why we used a logarithmic scale in subplot (a).

4.1.4 Dust masses

The observed dust masses of models with optically thick inner bumps are underestimated by up to an order of magnitude in comparison with the actual dust mass. This can be seen by the bump model’s influence on the various inferred dust masses, at an evolution time of our disk of 1 Myr. In Fig. 4, we plot the total dust mass of our simulation M_d,tot, the dust mass contained in grains larger than 0.1 mm M_{d>0.1 mm} and the observed dust mass M_d,obs over each model parameter. The latter is computed through the optically thin approximation and therefore solely mirrors the 1.3 mm flux of the models. We investigate M_{d>0.1 mm}, since at an observing wavelength of a few millimetre we do not trace µm-sized grains, so it is more insightful to compare M_d,obs to M_{d>0.1 mm} mm, instead of M_d,tot.

In subplot (a), we find our simulation dust masses decreasing for larger gap radii, since the closer in the bump is, the more dust it can halt from accreting onto the star. Comparing the observed dust mass with the actual ones, we find that the largest discrepancy occurs for the Inner Bump model, where M_d,obs is an order of magnitude smaller than M_d,tot. This is not surprising, since this model has the highest of all dust masses, mainly contained within the dust trap, therefore also having the highest optical depth at its peak location. Consequently, with the optically thin approximation we can only trace dust down to optical depths of roughly τ ≈ 1, meaning the largest fraction of the dust mass is hidden from our view. The Fiducial model and the Outer Bump model have both similar observed dust masses, though the former underestimates the real mass, while the latter traces its millimetre mass pretty well.

Looking at the Two Bumps model, we see that it has a higher observed dust mass than the Fiducial one, despite having a lower mm-sized dust mass. This is due to the second outer bump of the model. Instead of trapping more particles in the inner bump, therewith just increasing its optical depth, we now have an otherwise hidden fraction of the (mm) dust mass trapped in the outer bump which contributes to its 1.3 mm flux.

Moving on to plot 4b, we find the best match between the observed and actual dust masses for the Small A bump model. The models with an increased amplitude, on the one hand, increase their model dust masses. On the other hand, their observed dust masses decrease a small amount, which can be traced back to stronger trapping around the pressure maxima, leading to higher optical depths of the dust ring.

The last subplot (c), shows the dust mass dependency on the growth time of the bumps. As elaborated in previous sections, the model dust masses decrease for longer growth times, since the trapping is less efficient, if the bumps need longer times before they fully develop their local pressure maximum. This leads to a convergence of the two lines with M_{d, obs}, which does not drop as steeply.

4.1.5 Optical depths

The main result of this last subsection is that the observed optical depths of all dust traps in the inner regions of the disk are underestimated. We find this by comparing the absorption optical depth ${\tau }_{\text{abs}}^{\text{peak}}$, obtained through the dust surface density at the dust ring location of our DustPy models (Eq. (21)), to the observed optical depth ${\tau }_{\text{obs}}^{\text{peak}}$, calculated through the ring peak intensities of our synthetic observations (Eq. (22)). The corresponding plots are shown in Fig. 4.

First of all, we find all values of the observed optical depth to lie between τ_obs ≈ 0.2–0.7 for the models with bumps. On the other hand, the absorption optical depth usually spans values up to τ_abs ≤ 4 (lower panels). However, in the case of the Inner Bump it reaches up to τ_abs = 35. The absorption optical depth generally increases for larger amplitudes, shorter growth times and smaller gap radii.

The general trend of the curves in subplots (a)-(c) shows the same dependency on the bump parameters for the absorption and observed optical depths, as already seen in the upper panel of Fig. 4, between the total and observed dust masses, respectively. This is not surprising since they trace the same quantities: τ_abs and M_d,tot trace the dust surface density and τ_obs and M_d,obs the 1.3 mm flux density, both at the corresponding dust ring location.

In the lower left subplot (a), we additionally plotted the observed optical depths τ_obs,noscat (orange) for our models where we did not treat scattering in the radiative transfer calculations. As can be seen, this increases the peak flux and therewith the optical depth of our rings in the inner parts of the disk by roughly 40%, and slightly decreases it for the outer bump by 15% (see also Fig. 8). In the inner rings, we find a large amount of mm- to cm-sized dust pebbles which have a high scattering opacity (albedo). Furthermore, the scattering phase function of such grains has a strong forward peak, as already mentioned in Sect. 3.4. So, if we treat scattering in these optically thick models, a considerable amount of radiation gets scattered out of the line-of-sight leading to a decreased emission and optical depth compared to the models without scattering. In other words, scattering reduces the (optical) depth in our rings from where photons can escape (Zhu et al. 2019). This effect does not play a role in the optically thin (τ_abs < 1 ) outer parts of the disk, where the non-treatment of scattering even leads to a small increase in the observed optical depth. To sum up, we find that for nearly all models the observed optical depth is underestimated. Only for models without dust trapping (No Bump and Small A) and for the Outer Bump model, can τ_abs be well traced by τ_obs.

A final remark on the parameter effects: We also ran models where ‘planetary’ bumps form and then disappear after a finite lifetime of 0.5–1 Myr. However, those disappearing bumps display no lasting effect on the disk’s dust evolution in terms of observability. After a few hundred thousand years they cannot be differentiated anymore from the model without any bump, regardless of their lifetime and bump location.

4.2 Dynamic bumps

The main result of this section is that dynamic bumps, in contrast to ‘planetary’ ones, cannot efficiently halt the drift of pebbles in order to display high millimetre fluxes and low spectral indices observed in protoplanetary disks. We present the observational properties of our Dynamic Bumps models, which are set to mimic zonal flows in protoplanetary disks, in Table 4. In Fig. 5, we show the gas and dust surface densities at the evaluation time of the observations at 1 Myr, as well as the total dust masses over the models’ disk lifetime. The detailed dust size distribution of the models at 1 Myr can be found in Fig. A.2.

The dynamic models, shown in plots 5a and b, with bump lifetimes of 200 (Std), 400 (Long) and 1000 local orbital times, show quite low dust surface densities and nearly lost all of their initial dust mass already after 1 Myr of the disk’s evolution. After 3 Myr, less than 2 M_⊕ Me worth of dust are retained in the disk. In (b), it is worth noting that the doubled lifetime of the Long model does not slow down the exponential loss of the dust mass, compared to the Std setup. So it comes as no surprise that these models also show very low millimetre fluxes, in fact, they even have the lowest of all models in this work, even lower than the model with no bump, shown as a comparison in the plots. Thus, their respective millimetre spectral index of α_mm = 3.6 lies only slightly below the value of α_ISM ≃ 3.7, which is typically found in the optically thin emission of interstellar medium (ISM) grains (Tazzari et al. 2021b), where no dust growth is expected to happen.

On the contrary, the Static bumps model displays the highest dust mass and millimetre luminosity of all models in this work. The reason for this is its strong dust trap (inner bump) at 10 au, which is already fully developed after ≃20 kyr (700 orbits), due to the short orbital timescales here. Once the trap has formed, it efficiently halts the dust drift and therefore the disk’s dust mass does not significantly decrease anymore, as can be seen in the blue line of plot 5b. This is the important difference between the dynamic and the static model. If the dynamic ones would also have efficient dust traps in the inner parts of their disk, they would retain their dust mass and thereby reaching higher luminosities.

A last thing to point out about the Static model is that it has the highest total dust mass of all models, however, its observed dust mass is underestimated by a factor of five. Here, most of the disk’s dust mass is contained in the inner bump thus making it highly optically thick (${\tau }_{\text{abs}}^{\text{peak}}=15$). Hence, the optically thin approximation for measuring the dust mass is not applicable anymore.

To conclude this section, we shortly want to touch on the measured disk radii of the three different models. Looking at the Table 4 above, we find the R_d,68%, of the Static model to be the largest, but not significantly larger than for the the dynamic ones. Interestingly however, is the fact that R_d,90% is the lowest for the Static model, even though it displays a dust ring in Σ_d at ~100 au of Fig. 5a. This ring indeed, does not substantially contribute to overall flux anymore which is why we do not accurately trace it within the 90% of the emission and as a result, cannot trace the outer edge of the disk anymore.

Fig. 5 Surface density and dust masses plots of dynamic bumps models. (a) Gas (dashed line) and dust (solid line) surface densities at 1 Myr of evolution for the different dynamic bumps models and the model without any bumps for comparison. (b) Total dust mass of each model over the disk lifetime.

5 Discussion

5.1 Disk sizes

We find nearly all of the measured disk sizes (R_d,68%) coincide with the location of the farthest out pressure bump, as was already proposed by earlier studies (Pinilla et al. 2020; Zormpas et al. 2022). The exception are the models with multiple bumps, where the measured sizes (R_d,68%,90%) may indicate the location of the inner dust traps (see Fig. A.4, Two Bumps and Static models).

In a recent study, Franceschi et al. (2022) shows that in the presence of efficient dust traps, the dust disk’s outer edge (or ‘dust line’) only traces the trap location across different wavelengths. Therefore, the ‘dust line’ cannot be used as a tracer of radial drift anymore to determine the disk’s mass, as proposed by earlier studies of Powell et al. (2017, 2019). The results of Franceschi et al. (2022) are in agreement with the findings of our dust trapping models, when we calculate the emission radii for our intensity maps at 1.3 mm and compare them to the disk sizes at 3.0mm (see Fig. A.4). The maximal deviation for $R_{\text{d},68\%}^{1.3\text{\hspace{0.17em}mm}}\text{/}{R}_{\text{d},68\%}^{3.0\text{\hspace{0.17em}mm}}$ is 2% and for R_d,90% less than 5%. Again, one exception is the Two Bumps model, where $R_{\text{d},90\%}^{1.3\text{\hspace{0.17em}mm}}=100$ au traces the second bump, while $R_{\text{d},90\%}^{3.0\text{\hspace{0.17em}mm}}=41$ au coincides with the first.

A recent multi-wavelength survey of the Lupus star-forming region by Tazzari et al. (2021a) finds that their observed disk sizes are nearly the same across different millimetre wavelengths. The authors find a ratio between their sample’s mean millimetre dust disk radii of R_{68,3.0 mm}/R_{68,1.3 mm} = 0.94 ± 0.23, which matches our results well. This finding is however in contrast to dust evolution models (e.g. Birnstiel et al. 2012; Long et al. 2020) of smooth disks which predict much smaller radii at 1.3 mm and 3.0mm, than were observed (see their Fig. 4). As their survey was conducted at low resolution (0.3″), it is not possible to resolve hidden substructures within their disks. However, some of their (larger) sources overlap with the DSHARP sample, that resolved them and found substructures, which could explain why R68% is constant with wavelength at least for some disks.

To compare our results to the survey of Tazzari et al. (2021a), we cautiously note the caveat that we only consider a solar-type star with a fixed disk setup. Nevertheless to be as accurate as possible, we convolve our synthetic images with the survey’s larger beam size of 0.30″ (versus our 0.04″) and further assume a lower S/N of 200 (versus our 1000). Figure 6 displays an overview of the synthetic images of our models at high and low angular resolution.

We show our derived disk sizes at low and high angular resolution, in dependence of their 1.3 and 3.0 millimetre flux in the upper plots of Fig. 7 and then compare our models to the size–luminosity relation(s) found in Lupus by Tazzari et al. (2021a). In the lower plot (c) of Fig. 7, we compare our (low resolution) disk sizes to each of the Lupus sample and find them to span the same range of values. The discussion about the spectral indices and the luminosities is in the next subsection.

The decrease in resolution used to match the sample of Tazzari et al. (2021a) increases the measured disk sizes by a small factor of a few percent for most models. Only the Inner Bump model shows a significant increase in its measured disk size at the lower resolution (doubling its 68% emission radius at 1.3 mm, from 13 au to 24 au), which occurs because both the pressure bump located at 10 au and the disk’s physical extent are unresolved, leading to a spread of the emission over a wider area for the larger beam size.

Nevertheless, both measured radii of the Inner Bump model still lie within the correlation found in the Lupus sample, as shown in the plots of Fig. 7, both for the 1.3 mm radius and for the matching low spectral index. If, the smallest disks of the sample, that are comparable in size to the beam (R_d < 40 au), exhibit substructures, they cannot be resolved at this low angular resolution (see Fig. 6, Inner Bump and Fiducial models). Indeed, it might even be favourable that small disks do contain unresolved dust traps, as they also show the lowest spectral indices of the Lupus sample. The presence of pressure bumps can help to obtain such low values of α_mm, as has been shown by Pinilla et al. (2012) and is in agreement with our findings.

Building on the assumption of unresolved dust traps, we might find the smallest (Lupus) disks to be even smaller, when observed at higher angular resolution. This correlation has also been observed for some disks in Taurus by (Hendler et al. 2020, see Fig. 12), supporting our results. These points lead us to an important conclusion regarding the size-luminosity relation of disks. Since the measured luminosities and the sizes of larger disks, do not change or strongly decrease with higher angular resolution, we expect the slope of the relation to become steeper as the smallest disks are resolved. However, high angular resolution observations with high S/N are especially challenging to obtain for small disks, as most of them appear to be faint (e.g. Long et al. 2019; Kurtovic et al. 2021). Still, they do find their sample’s size-luminosity slope to steepen by including the resolved disks of DSHARP into their sample.

To conclude this section, we want to discuss the effect of neglecting planetary migration in our models, i.e. that we assume a fixed radial location of the bumps throughout the simulation. Looking at the corresponding planetary masses of our bump amplitudes in Table 1 and taking into account our disk setup, our planets might be subject to either type II or type III migration (e.g. Masset & Papaloizou 2003). In the former case, we would expect the bump to jointly migrate inwards with the planet, thus just decreasing the disk’s effective size. In the case of rapid runaway (type III) migration, however, the planet is able to form multiple pressure bumps on its discontinuous inward migration (see Fig. 5 of Wafflard-Fernandez & Baruteau 2020). The lifetime of such pressure bumps can be between several thousand years to 10⁵ yr depending on the disk’s turbulent viscosity. For the case when the disk viscosity is high (10⁻³), models of type III migration are similar to the models with disappearing bumps that we have investigated. These cases show no lasting effect on the disk’s dust evolution in terms of observability, except one observes the system at a ‘lucky’ time, while the structures are still present. However, when the viscosity is very low (10⁻⁵), the pressure bumps formed by planets in type III migration can be long-lived and reproduce several observable properties of disks as expected from static bumps.

We lastly note that the other bump model parameters, such as the growth times and amplitudes, do not have a major influence on the measured dust disk radii. The latter only if the amplitude factor is not large enough to act as an efficient dust trap, which usually sets the disks outer edge.

Fig. 6 Synthetic observations of models at low and high angular resolution. Images for λ = 1.3 mm, observed with high (0.04″) and low (0.30″) angular resolution, the according beams are indicated in the lower left corner of the plots. The black and white arcs in the images indicate the gap locations of the models.

5.2 Millimetre fluxes and spectral indices

In this section, we discuss which model parameters influence the millimetre fluxes of our disks and their according spectral index.

First of all, we compare how our measured millimetre fluxes fit into the size–luminosity relations shown in Fig. 7. We find all models with efficient dust traps in the inner regions of the disk (<40 au) have high enough intensities (and sizes) to fit within the confidence interval of the 1.3 mm relation. Yet, the models with outer bumps display too low 1.3 mm fluxes, compared to their large sizes, since their bumps lie too far out (r_gap > R_c) to trap enough dust. Similarly, the model without any dust trap, is also too faint for its medium disk size. Our models with more than one bump fit the size-luminosity relation well. Having more than one pressure bump in their disks, preferably in the inner parts, reduces their measured disk radii while simultaneously increasing their millimetre intensity.

In subplot (b) of Fig. 7, we also compared our models’ 3.0 mm flux (and radii) to the size-luminosity relation found at this wavelength by Tazzari et al. (2021a). They fit the correlation worse than for 1.3 mm in plot (a). While the disk sizes, at least for the models with bumps are constant with wavelength (as can be seen comparing the two plots), all of our models show too insufficient emission at 3.0 mm. Only the Static model, with its three bumps, still lies within the confidence interval of the correlation-line. The discrepancy of too low F_{3,0 mm} is again especially pronounced for the models without efficient dust traps and also for the outer bump model.

The low F_{3,0 mm} fluxes of our disks translate directly into high millimetre spectral indices, as shown in the lower plot of Fig. 7c. Here, we compare α_mm of our models to the values found in the Lupus sample by Tazzari et al. (2021a), which distribution is statistically indistinguishable from the star forming regions of Taurus and Ophiuchus (Tazzari et al. 2021b). We first note that all of our derived spectral indices are higher than what was observed for disks with similar sizes in Taurus. The only models that still fit within the given errors are the ones with fast-forming efficient dust traps at 10 and 30 au, respectively. Models with longer bump growth times, no or less efficient dust traps, and also the ones without a close-in bump, show too high spectral indices of α_mm ≥ 3.0. Some of their computed α_mm are even close to the value of 3.7, found in the interstellar medium (dashed line; Goldsmith et al. 1997). Nevertheless, we still observe the same general trend in subplot 7c, namely of an increased α_mm with larger disk size of our models.

These results suggest that in order to obtain the low spectral indices observed, we need disks to include long-lived and efficient dust traps, especially in their inner disk parts to halt the dust drift. This finding, together with low millimetre luminosities of our Dynamic Bumps model, makes us disfavour the zonal-flow origin of pressure bumps, in which the bumps induced by over-densities in the disk have a short lifetime. On the other hand, this then supports the planetary-origin of the dust traps, in which case the pressure bumps are indeed long-lived and thus can trap the large amounts of dust that we observe with high resolution observations of disks (Dullemond et al. 2018).

About the effect of scattering on our results, Fig. 8 shows the radial flux profiles of our three models with varying gap location, either treating scattering in the radiative transfer or not. In the left subplot (a), we find that the inclusion of scattering into our calculation decreases the 1.3 mm flux of our models, at the location of the optically thick inner dust rings. However, in the optically thin regions scattering increases the flux, as it is the case for the whole disk of the Outer Bump model. The same holds true for the 3.0 mm flux profile in subplot (b). The only differences here is that the dust ring of the Fiducial model is not optically thick anymore at this wavelengths, hence its peak emission flux does not decrease when scattering is included.

The implications of the (non-)treatment of dust scattering in models with optically thick regions or substructures has been studied in detail by Zhu et al. (2019). They find that dust scattering can significantly reduce the emission of optically thick disk regions, in line with our findings. They further state that the spectral index can be used to determine the optical thickness of disks. Using the DSHARP-opacities, they expect a jump of α_mm in the radial direction from optically thick (α_mm < 2.5) to optically thin regions (α_mm > 2.5) of the disks. This is exactly what we find, looking at the radial profile of the spectral index in subplot (c) of Fig. 8, where especially the Inner Bump model lies below 2.5 at its location of the optically thick dust ring. Zhu et al. (2019) also state that in the optically thick regime, the spectral index depends on the albedo rather than the opacity, as in the optically thin regime. Therefore, we include scattering in our models to be sure to correctly determine the emission of the optically thick dust rings.

Even though our optically thick substructures display regions of low α_mm in the radial profile, the overall values are still too high compared to profiles obtained by multi-wavelengths observations (e.g. see Tazzari et al. 2021a; Huang et al. 2018), where α_mm ≤ 3.5 in most cases, even in the outer disk regions. These generally unobserved high spectral indices of our models arise from the low flux values at 3.0 mm, since they overall contain too few mm-sized dust grains with high opacities. There are multiple issues regarding the low F_{3.0 mm}. Firstly, most of the dust mass in the disk is trapped in the dust rings. Due to optical depth effects, it can only contribute a certain portion to the emission. Large grains that contribute most to F_{3.0 mm} are trapped even closer around the pressure maximum, compared to grains of lower opacity. This leads to higher optical depths and lower emission.

Secondly, the bulk of the dust mass is stored in the largest centimetre-sized grains for all of our models, also for the one(s) with a bump at 80 au. Within the dust traps, the particle size is determined through the fragmentation limit, as displayed in Fig. 1. Within this regime, the maximum size a particle can grow to through collisions, is proportional to the square of the fragmentation velocity (see Eq. (14)). Moderately decreasing v_frag to values lower than 10 m s⁻¹ helps to grow a larger amount of mm-sized pebbles, whose drift is reduced and which usually have opacities an order of magnitude higher than grains of cm sizes (see Fig. A.3), therewith obtaining higher disk fluxes. This is indeed the case for a simulation we ran with v_frag = 3 m s⁻¹ which displays 62% and 24% higher fluxes at 1.3 mm and 3.0mm, respectively, compared to our Fiducial setup. However, the lower v_frag also leads to an increase of α_mm = 3.0 versus the original 2.7.

Laboratory experiments, that have revisited the sticking properties of dust aggregates in protoplanetary disks even support this idea of a lower fragmentation velocity (Gundlach et al. 2018; Musiolik & Wurm 2019; Steinpilz et al. 2019). Recently, Pinilla et al. (2021) studied disk models with such a fragmentation velocity of 1 m s⁻¹ , combined with different values of the turbulence parameters in disks. Even though, they have a similar disk setup and also conduct radiative transfer calculations with scattering, they obtain disk millimetre fluxes that are nearly an order of magnitude higher than the values we found. However, they use amorphous grains with a different opacity prescription, namely the one stated in Ricci et al. (2010), whose generally higher opacity values for mm-sized grains translate into their larger millimetre fluxes. This stresses the fact that the opacities of grains are still the largest uncertainty we face in modelling the disks millimetre appearances as discussed in the next subsection.

Fig. 7 Size-luminosity relations and spectral indices dependent on disk sizes. In (a) and (b), the fluxes are divided by the cosine of the inclination (i = 20 Å). The plotted 68% emission disk radii are found to decrease from low (0.30″) to high (0.04″) resolution for all models irrespective of wavelength which is highlighted by lines connecting both radii. The grey dashed lines show the relations for Lupus at λ = 1.3 mm and λ = 3.0 mm found by Tazzari et al. (2021a), where the grey area represents the 68% confidence interval around the median. (c) Spectral index (α1.3–3.0mm) over disk effective size (Rd,68%) at λ = 1.3mm with 0.30″ resolution, black dots show Lupus data also from Tazzari et al. (2021a).

Fig. 8 Radial flux and spectral index profiles for the three models with varying gap location, with and without scattering. (a and b) Radial flux profiles at λ = 1.3 and 3.0 mm. (c) The models’ radial profile of the spectral index. In all plots, the solid lines treat scattering in the radiative transfer, the dashed lines not.

5.3 The effect of the dust opacity

The composition of interstellar dust grains is still not well known. We initially chose the DSHARP opacities (which considers compact grains), since most of the studies we compare our results to, used this prescription as well (as Dullemond et al. 2018; Stammler et al. 2019; Zhu et al. 2019; Franceschi et al. 2022). Nevertheless, it is important to discuss the effect a different dust opacity model has, when modelling the observational properties of protoplanetary disks.

For this subsection, we calculated the opacity of the grains following the widely used prescription of Ricci et al. (2010). In contrast to the DSHARP opacity model by Birnstiel et al. (2018), the grains are assumed to be porous spheres with a vacuum volume fraction of 40%, composed of astronomical silicates (Draine 2003), carbonaceous material (Zubko et al. 1996) and water ice (Warren & Brandt 2008). First, this leads to a decreased dust material density of ρ_s,Ric = 0.96 g cm⁻³ versus ρ_s,Dsharp = 1.68 g cm⁻³. However, we are not rerunning our dust evolution models with this lower density, since we only investigate the change of opacity here.

The resulting millimetre absorption and scattering opacities for the two prescriptions can be compared in Fig. A.3. Notably, the Ricci absorption opacities are an order of magnitude higher than the ones of Birnstiel et al. (2018). Additionally, they also display no major difference between the two mm-wavelengths for larger than millimetre grains, whereas the 3.0 mm DSHARP opacities are lower than the 1.3 mm ones. These differences have a major impact on the calculated mm-fluxes and spectral indices of our modelled disks. The resulting observational properties for some chosen models, calculated using the Ricci opacities are shown in Table 5.

As it can be seen in the table, all models increase their millimetre fluxes by at least a factor of three, in the case of the No Bump and Std Dynamic models even by over an order of magnitude. Through the increases of F_{1.3 mm}, we also obtain larger observed dust masses which are in most cases now even higher than the actual model masses at 1 Myr. The highest increase is experienced for F_{3.0 mm}, as the Ricci absorption opacities are over an order of magnitude larger than the DSHARP prescription. This has a direct consequence for the calculated spectral index that is for all chosen models now α_mm ≤ 2.8. In turn, α_mm just mirrors the difference in the opacity values between 1.3 mm and 3.0 mm, as in nearly all models we grow a sufficient amount of pebbles which in both opacity prescriptions have high opacity values.

Figure 9 shows how the changed opacity prescription affects the models’ size-luminosity relations and spectral indices when compared to the Lupus sample of Tazzari et al. (2021a). As discussed above, the mm-luminosities for all models increase, however, their disk sizes stay roughly the same. This leads to the point that (smooth) models, which did not fit the relation in Fig. 7, fit it now the best, while the Fiducial and Static Bumps model, display now too high luminosities compared to their sizes. Zormpas et al. (2022) also finds that whether disks lie on the size-luminosity relation heavily depends on the chosen opacity model and that for DSHARP opacities preferably disks with sub-structures populate the relation. Looking at the spectral indices, we now find all models to lie within the range of the Lupus sample, even the ones which before did not.

The above discussed substantial differences in the observational properties highlight the need for more studies investigating the opacities of protoplanetary disks, in order to weigh up the usage of the various opacity prescriptions when calculating the disks’ observables. For instance, dust particles may grow more fluffy than assumed in this paper, until they grow and compact to planetesimals (Kataoka et al. 2013). The absorption and scattering mass opacity are highly influenced by the filling factor of the dust particles, which observationally can be distinguished by deep and high angular resolution at multiple wavelengths and polarisation observations of protoplanetary disks (Kataoka et al. 2014). However, if the dust particles are very fluffy, it is expected that they do not drift towards regions of high pressure in contraction to the observed substructures.

5.4 Optical depths and dust masses

The discrepancy between the observed and modelled dust masses is directly related to the optical depth at the peak location of our dust rings. This has already been seen by comparing the according variables in Fig. 4 that show the same dependencies on the bump model parameters. If the dust rings have τ_abs ≥ 1, the optically thin approximation breaks down and we find that the bulk of the dust mass, which is located in the rings, cannot accurately be traced anymore. This is shown in Fig. 10, where we display the particle density distribution for our models, over-plotted by the emission surface corresponding to an optical depth of τ = 1. Above this surface radiation or photons can escape the rings unimpeded, without being absorbed or scattered. Therefore, the dust mass that lies below the white contours is hidden from our view and can thus not be traced by observations anymore.

It is important to note that the height of the emission surface changes whether we include scattering in our radiative transfer calculations or not. Ignoring dust scattering leads to an underestimation τ, especially in the inner regions of the disks (<50 au) (Zhu et al. 2019). Similar to their results, we find that dust scattering reduces the emission from optically thick regions, as shown in Fig. 8, but increases it for optically thin ones. The latter is the case for the Outer Bump model where it leads to a further decrease in its optical depth, hence the model without scattering not even shows an emission surface of τ = 1. On the contrary, the rings of the Inner Bump and Fiducial model raise their observed optical depths to τ_obs = 0.58 and τ_obs = 1.02, respectively, when scattering is not included. This is an increase of over 40% for both models compared to the values including scattering, and is due to their increased emission at the peak location of the ring.

The general case is that we do not observe the dust rings with τ_abs > 1 to be optically thick, as all of our τ_obs < 1. However, the way of calculating the observed optical depth (through Eq. (22)) is already assuming the thermal emission of the dust to be optically thin (Dullemond et al. 2018), so it is not useful to measure the optical depth of optically thick rings with this formula. The argument of Dullemond et al. (2018), to use the optically thin assumption for their DSHARP sample is that they do not observe the source’s radial ring profiles to be flat-topped, which should be the case for highly optically thick rings. Yet, we do not observe a flat top of our radial ring emission profiles for our (highly) optically thick models in Fig. 8. This is likely due to the fact that such a flat-topped region is very narrow, such that it cannot be resolved by observations. The inclusion of scattering further reduces and broadens the surface, as well as the convolution with the Gaussian beam which additionally adds to the smoothing of the peaked emission profile and thus also to a decreased observed optical depth. In the end, it remains challenging to thoroughly determine the optical depths of observed disks with substructures.

To conclude the discussion on the optical depths of the rings, we want to point out that there are still large sources of uncertainty in calculating these values. Not only the temperature prescription for τ_obs, but also the grain composition and the thereof derived opacities used to calculate τ_abs are still not well known. Furthermore, we assume τ_abs to be the ‘actual’ optical depth of our models. This is not entirely correct, since we ignore the scattering opacity ${\kappa }_{\lambda }^{\text{scat}}$ of the grains which would significantly contribute to the then summed up extinction optical depth τ_ext = τ_abs + τ_sca. However, as scattering itself is angle dependent, this summation is not straightforward to do, so in first order it is reasonable to compare τ_obs to τ_abs, instead of τ_ext.

About dust disk masses, we first want to point out that we confirm the findings of Binkert et al. (2021), which also find the dust masses of optically thick inner bumps are underestimated by up to an order of magnitude. While we treat dust evolution in our one dimensional disk models, they carried out three dimensional two-fluid (gas + dust) hydrodynamic simulations to obtain their results.

Regarding the effects the initial conditions have on the obtained dust masses, the most important parameters are of course the initial disk mass and dust-to-gas ratio. Doubling the initial disk’s mass of the Fiducial model, we also found the total dust mass at 1 Myr to double. However, due to optical thickness of the dust ring, where the bulk of the mass is located, this leads only to an increase in the mm-fluxes (and M_d,obs) of about 35%.

Another very important parameter, we have not yet discussed, is the turbulent viscosity a. Since we induce our Gaussian bumps onto the a-profile, the moderate value of α₀ = 1e–3 makes sure that we do not experience too long growth times of our pressure bumps within the disks. The key parameter is the viscous timescale τ_visc ∝ (α₀Ω_k)⁻¹ which sets the lower limit on how fast the gas can viscously adapt to changes. If for example, we take very low a values, it may take up to 1 Myr for a bump in the outer region of the disk to reach its final amplitude. An increased turbulence value on the other hand, would lead to less dust trapping within the rings, as the higher dust diffusivity would smooth out the density gradients, eventually decreasing the overall dust mass.

Another initial bump model parameter that effects the disk masses is the time, after which we initially induce the bumps into our viscosity profile, namely at T_ini = 1 kyr. Most of studies are already including their pressure bumps at the beginning of the simulation and then let the viscous timescale set their according growth times (e.g. Dullemond et al. 2018; Pinilla et al. 2021). Recently, Kuznetsova et al. (2022) also argued for an early formation of pressure bumps in disks, which in their work are generated through anisotropic streams of infalling material. As we are particularly interested in the growth times of pressure bumps, we ran two models, similar to our Fiducial setup with a bump at 30 au, but where we start to grow the bump later than 1 kyr. For the first model, we induced the bump only at T_ini = 0.1 Myr of evolution which led to a decrease in total dust mass of 35% compared to the Fiducial model. In a second setup, we even increased T_ini to 0.5 Myr and further let the bump grow over the same period of time, in comparison to the fiducial value of T_grow = 0.1 Myr at 30 au. For this model, the according dust masses at 1 and 3 Myr of evolution, even decreased by more than a factor of eight. Both models show significantly less dust mass, than the Fiducial setup that would lead to even lower observed mm-fluxes and higher spectral indices for these models.

Regarding the formation time of zonal flows in global MHD-simulations, they can develop as fast as a few tens to hundreds of orbits, depending on the magnetisation and ionisation of the gas. They then remain stable in the disk for several hundreds of orbits. In a recent study by Riols et al. (2020) their zonal flow structures even remains steady over their entire simulation time (~1000 orbits). We ran another simulation taking this value for T_life together with the formation time of their zonal flows (~200 orbits). It turned out that the increased lifetime does not help in halting the dust loss (see Fig. 5). The important point is, that the inner-most zone passes ~25 bump life cycles till 1 Myr in (each 38 kyr at 10 au), which sets the bottleneck of how much dust the disk looses towards the star. Therefore, we only obtain a total dust mass of 12 M_⊕ for the above model, still an order of magnitude less than in the Static one.

To conclude, the observed dust masses of disks are still hard to measure, especially as we cannot be sure whether observed substructures are optically thick or not, therewith hiding a bulk of the dust mass contained within disks. However, optically thick substructures might ease the discrepancy between the observed dust masses in mature class II disks, compared to the solid material actually needed to reproduce observed exoplanetary systems (Tychoniec et al. 2020). Another possibility is that part of the dust mass could be contained within larger boulders of metre to kilometre sized objects.

Fig. 9 Size-luminosity relations and spectral indices dependent on disk sizes. As Fig. 7 but using the Ricci et al. (2010) opacities.

Fig. 10 Vertical and radial dust density distribution. Shown are the models with varying gap locations after 1 Myr of evolution, over-plotted by the contours of the optical depth emission surface corresponding to τ = 1 with either treatment of full (solid) or no (dashed) scattering.

5.5 Planetesimal formation

In this last subsection, we shortly want to discuss the impact planetesimal formation has on our Fiducial model, where it might reduce the dust mass stored within the ring, converting it into ‘invisible’ planetesimals and therewith reducing its optical depth.

Dust rings may be ideal places to form metre to kilometre sized objects, as the local dust-to-gas ratio є is strongly enhanced here. This might trigger the so-called ‘streaming instability’ (e.g. Youdin & Goodman 2005) that can lead to the process of gravoturbulent planetesimal formation (Johansen et al. 2006). We refer to Sect. 7.3 of Dullemond et al. (2018) for a more detailed discussion about the mechanisms that can lead to the streaming instability in dust rings.

To this end, we simply follow the setup of Stammler et al. (2019) that builds on the planetesimal model of Schoonenberg et al. (2018). Here, planetesimal formation is triggered as soon as a characteristic value є_crit = 1.0 of the dust-to-gas ratio at the midplane is reached. If this is the case, a fraction ζ = 0.1 of the dust surface density per settling time scale of each dust mass bin will be converted into the surface density of planetesimals. Similarly, we do not further evolve the planetesimals’ surface density, even though planetesimal-planetesimal collisions could lead to dust reproduction and hence a lower planetesimal mass.

We take our Fiducial model, with the bump at 30 au and now add the formation of planetesimals to the model. The dust size distribution of the according model, can be found in Fig. A.1. As expected, we find the dust density to be decreased only at the location of the ring (by a factor of ${\text{Σd}}_{\text{fid}}^{}{\text{/Σd}}_{\text{pld}}^{}~3$), compared to the size distribution of the Fiducial model. Apart from the density at the ring location, there are no differences in the dust size distribution between the models.

In Fig. 11, we plotted the total masses of the models stored in dust and planetesimals over time. First, we see that the difference in dust mass between the models does not increase further after about 1.5 Myr of evolution. Already at 1 Myr, we find over half of the total dust mass of the Fiducial model to be converted into planetesimals, namely M_pltes = 47 M_⊕, even more than what is stored in dust $M_{\text{d}}^{1\text{Myr}}=42M_{\oplus }$. This leads to decreased fluxes at 1.3 and 3.0mm of 15% and 35%, respectively. As F_{3.0 mm} drops stronger in emission than F_{1.3 mm}, this directly impacts the spectral index which increased to α_mm = 3.0.

But how about the optical depths of the model? Checking again the absorption optical depth at the peak location of the dust rings, we find ${\tau }_{\text{abs}}^{\text{plts}}|,1.4$ versus ${\tau }_{\text{abs}}^{\text{fid}}=3.5$, a decrease of 60%. However, calculating the optical depth from the intensity peak of the synthetic observation, we only find ${\tau }_{\text{abs}}^{\text{plts}}=0.55$ versus ${\tau }_{\text{obs}}^{\text{fid}}=0.70$, a small reduction of ~20%.

The same trend can be found using the higher Ricci et al. (2010) opacities (see Table 5 for obs. prop.), here we find both rings highly optically thick ${\tau }_{\text{abs}}^{\text{plts}}=14$ versus ${\tau }_{\text{abs}}^{\text{fid}}=36$, still the planetesimal model with a decrease of over 60%. For the observed optical depths the picture is harder to interpret, as for our Fiducial model, we were not able to determine ${\tau }_{\text{obs}}^{\text{fid}}$ anymore, since the peak emission exceeds the prescribed blackbody emission. Therefore, we increase the blackbody temperature prescription of Eq. (17) by 5 K to 165 K at 1 au, to obtain ${\tau }_{\text{obs}}^{\text{plts}}=2.2$ versus ${\tau }_{\text{obs}}^{\text{fid}}=\text{4.2}$. On contrary to the DSHARP-opacities, which were also used by Stammler et al. (2019), we now find both rings to be optically thick, even though planetesimal formation lowered τ_obs by nearly a factor of two. Still, the problem persists that τ_obs underestimates the ‘actual’ absorption optical depth τ_abs of the rings, for the Ricci-prescription even greater than for the DSHARP-opacities.

Stammler et al. (2019) proposed that planetesimal formation provides a mechanism to explain the low observed optical depths of dust rings in the DSHARP-sample (Dullemond et al. 2018), as it naturally regulates the dust density in the rings. In our adaption, we similarly find the absorption optical depth of our dust ring to decrease, but only at the cost of lower millimetre fluxes and more importantly an increased spectral index of the disk (similarly using the DSHARP-opacities), which is not supported by observations (though if we use the high Ricci et al. 2010 opacities the α_mm remains unchanged). Additionally, Stammler et al. (2019) calculated their intensity profile using Eq. (22) without the treatment of scattering. We on the contrary, performed full radiative transfer calculations with an-isotropic scattering to obtain synthetic observations. From the convolved intensity profiles, we further obtain the observed optical depths of the ring(s) similar to Dullemond et al. (2018) and find them in both cases to be optically thin.

To conclude this section, we first note that it is challenging to obtain correct optical depth values from observations, as we measure them already assuming optical thinness and additionally, as the peak of the emission and thus derived optical depth is also depending on the beam-size, decreasing with larger beams. Secondly, we should be cautious as to compare these observed values to the ones of models without scattering, as scattering can significantly reduce the emission and therewith observed optical depths, at least for those of optically thick dust rings in inner disk regions (≤40 au), as we have shown in the above sections. These points make it hard to conclude whether planetesimal formation is actually at work in dust rings, regulating their optical depths, or if the fine-tuned values of τ are just a result of the error-prone determination of the observed optical depths. A recent study by Carrera & Simon (2022) also presents evidence that the streaming instability is indeed unlikely to form planetesimals from mm-grains inside axisymmetric pressure bumps.

Fig. 11 Dust mass evolution of the model with planetesimal formation. Total mass of the fiducial and planetesimal formation model in dust (blue and green) and in planetesimals (brown) over the disk lifetime.

6 Conclusions

In this work, we investigated how dynamic pressure bumps affect the observational properties of protoplanetary disks. Furthermore, we aimed to differentiate between the planetary- versus zonal flow-origin of pressure bumps. To this end, we performed gas and dust evolution simulations, setting up models with varying pressure bump features, including their amplitude and location, growth time and number of bumps. We subsequently ran radiative transfer calculations, including the full treatment of anisotropic scattering, to create synthetic images from which we inferred the disk’s observed properties. Finally, we compared our results to multi-wavelength observations of protoplanetary disks in the millimetre dust continuum.

The main conclusions of this paper are summarised as follows:

1. The location of the farthest out pressure bump determines the dust disk effective size, calculated through its emission radius enclosing (68% or) 90% of the flux. This delivers a possible answer to why disks in low-resolution surveys appear to have the same sizes across different wavelengths (Tazzari et al. 2021a);

2. Dust traps need to form early (<0.1 Myr), fast (on viscous timescales) and must be long lived (>Myr) in order to obtain low spectral indices of disks. This favours the planetary-origin of pressure bumps and supports the idea that planet formation may already occur in the class 0-I phases of circumstellar disks (e.g. Tychoniec et al. 2020) and that their migration is type II during most of the disk’s lifetime;

3. Unlike the planetary bumps, our dynamic bump model with dis- and re-appearing pressure bumps of finite lifetimes can only reproduce the disk’s high fluxes and low spectral indices observed with the highest prescription of opacities that we investigated, thus disfavouring the zonal flow-origin of most of the observed disk substructures;

4. The wavelength-independent sizes of unresolved small disks (R_dust ≤ beam-size) with low spectral indices may indicate that they contain unresolved dust traps. By resolving them, we could measure the disks to be significantly smaller in size. This would lead to a steepening of the slope for the size-luminosity relation found in disks, since the determination of the luminosity (in general) and the sizes of larger disks, are less affected by the chosen resolution;

5. We find a large discrepancy between the observed and modelled optical depths of dust rings within our disks, especially for those in inner regions (<40 au), independently of the opacity model. For the DSHARP opacities, the observed values all lie between τ_obs ≈ 0.2–0.7 for the models with bumps, even though some of them are optically thick, reaching optical depths in the models of τ_abs ≈ 2–35. This difference is because the determination of τ_obs implicitly assumes the emission to be optically thin. Additionally, the non-treatment of scattering underestimates the optical depth of dust rings, particularly in the inner part of the disks;

6. The observed dust masses of disks with optical thick inner bumps (≤40 au) are underestimated by up to an order of magnitude. This might ease the discrepancy between observed dust masses of mature class II disks and the solid material actually needed to reproduce the observed exoplanetary systems (Tychoniec et al. 2020). However, these results are particularly affected by the chosen opacity model;

7. Planetesimal formation within dust rings can help to reduce the modelled optical depths of rings, but only at the cost of lower millimetre luminosities and high spectral indices, when using the DSHARP-opacities. The prevalent determination of observed optical depths (assuming optical thinness), however, makes it hard to conclude whether planetisimal formation is actually at work in dust rings, or if the observed fine-tuned values of τ are just a result of the procedure determining them in the first place;

8. Most of the inferred observational disk properties are sensitive to the chosen opacity prescription (Birnstiel et al. 2018). Using another prescription (Ricci et al. 2010) leads to higher fluxes, dust masses and optical depths, as well as lower spectral indices (α_mm ≤ 2.8). The dust disk effective sizes for models with long-lived pressure bumps are less affected by the choice of opacity.

In summary, our findings favour the planetary- over the zonal flow-origin of ring and gap substructures in the dust continuum of protoplanetary disks. Furthermore, this supports the idea that planet formation may already occur in early Class 0-I stages of circumstellar disks. It remains an open question if planetesimals can be formed in our zonal flow model, who then can deliver the initial seeds to form planetary cores in early times, depending on their formation timescales.

Acknowledgements

We would like to thank referee Clément Baruteau for his very constructive feedback, as well as Mario Flock for the fruitful discussion regarding zonal flows. The authors acknowledge funding from the Alexander von Humboldt Foundation in the framework of the Sofja Kovalevskaja Award endowed by the Federal Ministry of Education and Research, 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) under grants 361140270, 325594231, and Germany’s Excellence Strategy – EXC-2094 – 390783311. Software: DustPy (Stammler & Birnstiel 2022), Matplotlib (Hunter 2007), Numpy (van der Walt et al. 2011), OpTool (Dominik et al. 2021), RADMC-3D (Dullemond et al. 2012), Scipy (Virtanen et al. 2020).

Appendix A Supporting figures

Fig. A.1 Dust size distributions for the ‘planetary’ bump models at 1 Myr. The contours show the growth limits for drift (orange) and fragmentation (cyan), as well as the the particle sizes corresponding to St = 1 (white) and 1 mm (grey). The parameters for the displayed models can be found in Table 3.

Fig. A.2 Dust size distributions dynamic bump models. As Fig. A.l but for the dynamic bumps models.

Fig. A.3 Millimetre opacities over dust grainsize. Absorption and scattering opacities over dust grainsize for the observation wavelengths of λ = 1.3 (blue) & 3.0mm (red). The composition of the grains leading to the displayed opacities are adopted from the DSHARP paper of Birnstiel et al. (2018) and the paper of Ricci et al. (2010)

Fig. A.4 Radial millimetre intensity profiles. Profiles for synthetic images at 1.3 mm (blue) and 3.0 mm (red), both convolved with a Gaussian beam with a FWHM of 0.04″×0.04″. The dashed and dotted-dashed lines show the 68% and 90% emission disk radii of the models, respectively.

References

All Tables

All Figures

	Fig. 1 Dust size distribution for the Fiducial model at 1 Myr. The contours show the growth limits for drift (orange) and fragmentation (cyan), as well as the particle sizes corresponding to St = 1 (white) and 1 mm (grey).
In the text

	Fig. 2 Dust disk radius parameter plots. Dust disk radii enclosing 68% (green) or 90% (blue) of the 1.3 mm cumulative flux plotted over the gap location, bump amplitude and bump growth time (from a-c). In sub-plot (a), the cross markers correspond to the Two Bumps model and the dashed-dotted line shows the 1:1 -line.
In the text

	Fig. 3 Flux and spectral index parameter plots. (Upper) Disk millimetre fluxes at 1.3 mm (blue) and 3.0 mm (green) and according millimetre spectral indices (lower) plotted over the gap location, bump amplitude and bump growth time (from a-c).
In the text

Fig. 4 Dust mass and optical depth parameter plots. (Upper) Total dust mass (blue), dust mass for grains larger than 0.1 mm (orange) and observed dust mass (green) plotted over the gap location, bump amplitude and bump growth time (from a-c). (Lower) Absorption and observed peak optical depth also in dependence on bump model parameters. The Inner Bump model reaches τabs = 35, which is why we used a logarithmic scale in subplot (a).

In the text

	Fig. 5 Surface density and dust masses plots of dynamic bumps models. (a) Gas (dashed line) and dust (solid line) surface densities at 1 Myr of evolution for the different dynamic bumps models and the model without any bumps for comparison. (b) Total dust mass of each model over the disk lifetime.
In the text

	Fig. 6 Synthetic observations of models at low and high angular resolution. Images for λ = 1.3 mm, observed with high (0.04″) and low (0.30″) angular resolution, the according beams are indicated in the lower left corner of the plots. The black and white arcs in the images indicate the gap locations of the models.
In the text

Fig. 7 Size-luminosity relations and spectral indices dependent on disk sizes. In (a) and (b), the fluxes are divided by the cosine of the inclination (i = 20 Å). The plotted 68% emission disk radii are found to decrease from low (0.30″) to high (0.04″) resolution for all models irrespective of wavelength which is highlighted by lines connecting both radii. The grey dashed lines show the relations for Lupus at λ = 1.3 mm and λ = 3.0 mm found by Tazzari et al. (2021a), where the grey area represents the 68% confidence interval around the median. (c) Spectral index (α1.3–3.0mm) over disk effective size (Rd,68%) at λ = 1.3mm with 0.30″ resolution, black dots show Lupus data also from Tazzari et al. (2021a).

In the text

	Fig. 8 Radial flux and spectral index profiles for the three models with varying gap location, with and without scattering. (a and b) Radial flux profiles at λ = 1.3 and 3.0 mm. (c) The models’ radial profile of the spectral index. In all plots, the solid lines treat scattering in the radiative transfer, the dashed lines not.
In the text

	Fig. 9 Size-luminosity relations and spectral indices dependent on disk sizes. As Fig. 7 but using the Ricci et al. (2010) opacities.
In the text

	Fig. 10 Vertical and radial dust density distribution. Shown are the models with varying gap locations after 1 Myr of evolution, over-plotted by the contours of the optical depth emission surface corresponding to τ = 1 with either treatment of full (solid) or no (dashed) scattering.
In the text

	Fig. 11 Dust mass evolution of the model with planetesimal formation. Total mass of the fiducial and planetesimal formation model in dust (blue and green) and in planetesimals (brown) over the disk lifetime.
In the text

	Fig. A.1 Dust size distributions for the ‘planetary’ bump models at 1 Myr. The contours show the growth limits for drift (orange) and fragmentation (cyan), as well as the the particle sizes corresponding to St = 1 (white) and 1 mm (grey). The parameters for the displayed models can be found in Table 3.
In the text

	Fig. A.2 Dust size distributions dynamic bump models. As Fig. A.l but for the dynamic bumps models.
In the text

	Fig. A.3 Millimetre opacities over dust grainsize. Absorption and scattering opacities over dust grainsize for the observation wavelengths of λ = 1.3 (blue) & 3.0mm (red). The composition of the grains leading to the displayed opacities are adopted from the DSHARP paper of Birnstiel et al. (2018) and the paper of Ricci et al. (2010)
In the text

	Fig. A.4 Radial millimetre intensity profiles. Profiles for synthetic images at 1.3 mm (blue) and 3.0 mm (red), both convolved with a Gaussian beam with a FWHM of 0.04″×0.04″. The dashed and dotted-dashed lines show the 68% and 90% emission disk radii of the models, respectively.
In the text