The impact of dynamic pressure bumps on the observational properties of protoplanetary disks

Context. Over the last years, large (sub-)millimetre surveys of protoplanetary disks in different star forming regions have well constrained the demographics of disks, such as their millimetre luminosities, spectral indices, and disk radii. Additionally, several high-resolution observations have revealed an abundance of substructures in the disk’s dust continuum. The most prominent are ring like structures, which are likely caused by pressure bumps trapping dust particles. The origins and characteristics of these pressure bumps, nevertheless, need to be further investigated. Aims. The purpose of this work is to study how dynamic pressure bumps affect observational properties of protoplanetary disks. We further aim to differentiate between the planetary-versus zonal flow-origin of pressure bumps. Methods. We perform one-dimensional 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 run radiative transfer calculations to obtain synthetic images, from which we obtain the different quantities of observations. Results. We find that the outermost pressure bump determines the disk’s dust size across different millimetre wavelengths and confirm that the observed dust masses of disks with optically thick inner bumps ( < 40au) are underestimated by up to an order of magnitude. Our modelled dust traps need to form early (<0.1Myr), fast (on viscous timescales), and must be long lived (>Myr) to obtain the observed high millimetre luminosities and low spectral indices of disks. While the planetary bump models can reproduce these observables irrespectively of the opacity prescription, the highest opacities are needed for the dynamic bump model, which mimics zonal flows in disks, to be in line with observations. Conclusions. Our findings favour the planetary-over the zonal flow-origin of pressure bumps and support the idea that planet formation already occurs in early class 0–1 stages of circumstellar disks. The determination of the disk’s effective size through its outermost pressure bump also delivers a possible answer to why disks in recent low-resolution surveys appear to have the same sizes across different millimetre wavelengths.


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 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(Alcalá et al. , 2017Manara et al. 2014Manara et al. , 2016Manara et al. , 2017. In addition, (sub-)millimetre observations with the Submillimetre Array (SMA) and the Atacama Large Millimetresubmillimetre 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. 2016Ansdell et al. , 2017Ansdell et al. , 2020Barenfeld 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. 2018c;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 ). 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, Fig. 7 & 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 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-)hydrodynamical 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 flowsorigin 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 Section 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 the Section 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 Section 6 we draw our main conclusions.

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): 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) For the kinematic viscosity ν, we use the so called α-prescription for viscous accretion disks from Shakura & Sunyaev (1973) that states 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 Ω k = GM * /r 3 being the Kepler frequency. The isothermal sound speed can be written as c s = k B T/µm 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.

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 with the gas midplane density ρ g = Σ g / √ 2πh 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 α 0 as where F i (t, r) for each bump i is given by 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 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 α 0 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 & 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 overdensities 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 Equation (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 then 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.

Dust evolution
The evolution of the dust is of particular importance, since its distribution determines the disk's appearance in the mmcontinuum 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 midplane, 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: 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): Stokes regime 1, where the mean-free path of the gas is given by λ mfp = 1/( √ 2 n σ H 2 ), with the gas' number density n = ρ g /µ m p and the cross-section for molecular hydrogen σ H 2 = 2 · 10 −15 cm 2 . 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): with η = − (1/2) (h g /r) 2 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 = ν/(1 + St 2 ) (Youdin & Lithwick 2007). The overall dust evolution can then be described by an advection-diffusion equation following Birnstiel et al. (2010): 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 backreaction 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).

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 Ω k St (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 d = h g √ α 0 /(α 0 + St) (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 also using the prescription of Birnstiel et al. (2010) for the dust scale height: 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 .

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 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 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 ) , where we neglect dust settling which can also set the maximum 10 1 10 2 r (au) grain size . 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 2. 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)).

Disk setup
We carried out our disk simulations with the code DustPy 1 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): with γ = 1 and the normalisation factor Σ 0 = M disk /(2πR 2 c ) 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 where the stellar effective temperature T eff = 3500 K, the stellar radius R * = 2 R 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 = 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 0 = 1 µm. We assume the grains to be compact spheres with a mass given by m = 4πa 3 ρ s /3 and their material density ρ s = 1.675 g/cm 3 , 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 1. Initial dust-to-gas ratio 0.01 v frag Fragmentation velocity 10 m/s

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 2. Initially, we include the bumps at T ini = 1000 yrs 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: since we induce them in α. 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 h g 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.

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

Radiative transfer calculations
To infer the observational appearance of the disk, we perform radiative transfer calculations with the code RADMC-3D 2 (Dullemond et al. 2012). As input, we take the 3D dust density distribu-tion (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 3 (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 7 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-spreadfunction 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. 2018b) 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, with respect to the maximum of the emission, of 1000. We choose this peak SNR, 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 SNR for the brightest disks, however for fainter disks this might not be achievable without too long integration times.
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): where d is the distance to the source, F 1.3mm is the flux taken at λ = 1.3 mm and B ν 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 κ ν = 2.3 cm 2 g −1 × (ν/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: with the absorption opacity κ i,abs ν calculated from the DSHARP opacities  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) where both the intensity I ν and the temperature of the blackbody function B ν get evaluated for λ = 1.3 mm. For the temperature T d , we use the profile introduced in equation (17).

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 parameters are kept constant. An overview of the different models with their corresponding parameters and observational properties can be found in Table 3.

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 of the Appendix.
In plot 2 a) 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 .

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.0mm drops more than the F 1.3mm , as longer wavelengths trace larger particles.
A last thing to point out in Fig. 3 a) 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.3mm by 20%, while F 3.0mm 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. 4 c)). 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 for larger growth times than F 1.3mm which will have a direct effect on the spectral indices, as we will discuss in the following subsection.

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.0 mm-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.3mm (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.0mm .
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.0mm compared to F 1.3mm , whose ratio in the end determines α mm .

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.1mm 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.1mm , 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.1mm , 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 4 b), 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.

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 τ peak abs , obtained through the dust surface density at the dust ring location of our DustPy models (Eq. (21)), to the observed optical depth τ peak obs , 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) show 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 cmsized 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 & 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.

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 of the Appendix.
The dynamic models, shown in plot 5 a) & 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 ⊕ 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 kyrs (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 5 b). 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 (τ peak abs = 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. 5 a). 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.

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 al-ready proposed by earlier studies 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. (2017Powell et al. ( , 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.0 mm (see Fig. A.4 in the Appendix). The maximal deviation for R 1.3mm d,68% /R 3.0mm d,68% is 2% and for R d,90% less than 5%. Again, one exception is the Two Bumps model, where R 1.3mm d,90% = 100 au traces the second bump, while R 3.0mm d,90% = 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.0mm /R 68,1.3mm = 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.0 mm, 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 R 68% 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 SNR of 200 (versus our 1000). Fig. 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 sizeluminosity 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 Article number, page 11 of 22  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 SNR 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 2 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 5 years depending on the disk's turbulent viscosity. For the case when the disk viscosity is high (10 −3 ), 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, ex-  Fig. 7: Size-luminosity relations and spectral indices dependent on disk sizes. In a) + 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 (R d,68% ) at λ =1.3 mm with 0.30" resolution, black dots show Lupus data also from Tazzari et al. (2021a). cept one observes the system at a 'lucky" time, while the structures are still present. However, when the viscosity is very low (10 −5 ), 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.

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 correlationline. The discrepancy of too low F 3.0mm is again especially pronounced for the models without efficient dust traps and also for the outer bump model. The low F 3.0mm fluxes of our disks translate directly into high millimetre spectral indices, as shown in the lower plot of Fig. 7  c). 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 form 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 7 c), 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 overdensities 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 .
About the effect of scattering on our results, Fig. 8 show 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.0mm . 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.0mm 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 mmsized 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.0 mm, 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.

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 3 versus ρ s,Dsharp = 1.68 g/cm 3 . 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 of the Appendix. Notably, the Ricci absorption opacities are an order of magnitude higher than the ones of Birnstiel et al. (2018). 2.5 NOTE -All observational properties have columns similarly labelled as columns in Table 3. For the Fiducial model τ peak obs could not be computed, as its peak intensity was larger than the blackbody emission.
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 mmfluxes 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.3mm , 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.0mm , 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. Fig. 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.

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, overplotted 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 , 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 10 1 10 0 F 1.3mm / cos(i) (Jy) 10 1  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 κ 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 α. Since we induce our Gaussian bumps onto the α-profile, the moderate value of α 0 =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 ∝ (α 0 Ω k ) −1 which sets the lower limit on how fast the gas can viscously adapt to changes. If for example, we take very low α 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 a 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 MHDsimulations, 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 r (au) 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 12M ⊕ 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.

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 section 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 the Appendix. As expected, we find the dust density to be decreased only at the location of the ring (by a factor of Σ fid d /Σ plt d ∼ 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 = 47M ⊕ , even more than what is stored in dust M 1Myr d = 42M ⊕ . This leads to decreased fluxes at 1.3 and 3.0 mm of 15% and 35%, respectively. As F 3.0mm drops stronger in emission than F 1.3mm , 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 τ plts abs = 1.4 versus τ fid abs = 3.5, a decrease of 60%. However, calculating the optical depth from the intensity peak of the synthetic observation, we only find τ plts obs = 0.55 versus τ fid obs = 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 τ plts abs = 14 versus τ fid abs = 36, still the planetesimal model with a decrease of over 60%. For the ob-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. 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.