Hints on the origins of particle traps in protoplanetary disks given by the Mdust − M? relation

Context. Demographic surveys of protoplanetary disks, carried out mainly with the Atacama Large Millimeter/submillimete Array, have provided access to a large range of disk dust masses (Mdust) around stars with different stellar types and in different star-forming regions. These surveys found a power-law relation between Mdust and M? that steepens in time, but which is also flatter for transition disks (TDs). Aims. We aim to study the effect of dust evolution in the Mdust − M? relation. In particular, we are interested in investigating the effect of particle traps on this relation. Methods. We performed dust evolution models, which included perturbations to the gas surface density with different amplitudes to investigate the effect of particle trapping on the Mdust −M? relation. These perturbations were aimed at mimicking pressure bumps that originated from planets. We focused on the effect caused by different stellar and disk masses based on exoplanet statistics that demonstrate a dependence of planet mass on stellar mass and metallicity. Results. Models of dust evolution can reproduce the observedMdust−M? relation in different star-forming regions when strong pressure bumps are included and when the disk mass scales with stellar mass (case of Mdisk = 0.05M? in our models). This result arises from dust trapping and dust growth beyond centimeter-sized grains inside pressure bumps. However, the flatter relation of Mdust −M? for TDs and disks with substructures cannot be reproduced by the models unless the formation of boulders is inhibited inside pressure bumps. Conclusions. In the context of pressure bumps originating from planets, our results agree with current exoplanet statistics on giant planet occurrence increasing with stellar mass, but we cannot draw a conclusion about the type of planets needed in the case of low-mass stars. This is attributed to the fact that for M? < 1M , the observed Mdust obtained from models is very low due to the efficient growth of dust particles beyond centimeter-sizes inside pressure bumps.


Introduction
Censuses performed on the thousands of exoplanets that have been discovered thus far have revealed a large diversity of planetary architectures, along with a variety of trends that may be the result of the common physical processes ruling the formation and evolution of planets in disks around young stars. Exoplanet statistics show that giant planet occurrence increases with stellar metallicity and mass (e.g., Fischer & Valenti 2005;Johnson et al. 2010), while sub-Neptune type planets are more common around M-dwarfs than around Sun-like stars (e.g., Mulders 2018). These observed trends stand as open challenges for current theories of planet formation, examining stellar and disk properties; in particular, questioning how the effects of the ratio between stellar mass and disk mass affect the final outcome of planet formation.
Current observations of protoplanetary disks with the Atacama Large Millimeter/submillimete Array (ALMA) provide demographic surveys of disks (Ansdell et al. 2016(Ansdell et al. , 2017(Ansdell et al. , 2018Barenfeld et al. 2016Barenfeld et al. , 2017Pascucci et al. 2016;Cazzoletti et al. 2019;Cieza et al. 2019;Long et al. 2019;van Terwisga et al. 2019) that aid in the understanding of the missing keys for current planet formation models. Using continuum observations of disks, we can currently access to different disk properties, including a large range of disk dust masses -an important property to understand the amount of material available in young disks to form different types of planets. The most common method to determine disk mass is to assume that the dust emission is optically thin over most of the disk volume and, therefore, that the detected flux is proportional to the dust disk mass (e.g., Beckwith et al. 1990). This dust mass is usually used as a tracer of the total disk mass by assuming a constant dust-to-gas ratio, conventionally the interstellar medium ratio, that is, 1/100.
Obtaining information about the amount of gas available in disks is challenging because faint optically thin molecular lines, such as 13 CO and C 18 O, must be observed and their intensity can be also affected by isotope selective processes (Bruderer et al. 2012;Miotello et al. 2014Miotello et al. , 2016 and the unknown carbon abundance relative to H 2 (e.g., Kama et al. 2016;Schwarz et al. 2016). In addition, even 13 CO and C 18 O may be optically thick inside the CO ice line at the disk midplane, resulting in an underestimation of the disk mass (Zhang et al. 2017), especially for warm disks around Herbig stars. In these cases, observations of rare CO isotopologues, such as 13 C 17 O are needed (although chal-Article number, page 1 of 16 arXiv:2001.11045v2 [astro-ph.EP] 7 Feb 2020 A&A proofs: manuscript no. pinilla_Must_Mstar lenging) for a more robust estimation of gas disk masses ) and, therefore, the best knowledge that we have so far about the disk material still comes from the dust emission.
Information obtained about the dusty material in disks from observations is, however, prone to uncertainties as well since it is based on several assumptions, such as the disk temperature (which is typically assumed to be 20 K but could vary depending on disk size; Hendler et al. 2017), and the dust opacity. The dust opacity is the major source of uncertainty because it depends on the grain size and composition, which are unknown in protoplanetary disks (e.g., Birnstiel et al. 2018). Nevertheless, when comparing millimeter-observations of protoplanetary disks with models of disk evolution, the same assumptions can be taken in order to facilitate the comparison and understand the effect of different crucial physical processes in the formation of planets, such as grain growth and drift.
Due to the fast radial drift that millimeter-and centimeter-sized particles undergo in the outer parts (>20 au) of protoplanetary disks, models of dust evolution at million-year timescales contradict millimeter observations of disks (Birnstiel et al. 2010;Pinilla et al. 2012). This is because the models predict a millimeter spectral index (an indicator of grain growth) that is higher than that which is observed for most protoplanetary disks in different star forming regions (e.g., Ricci et al. 2010;Testi et al. 2014). Reducing or completing the suppression of the radial drift of dust particles, for instance, with the presence of particle traps or pressure bumps, is necessary in order to keep millimeter grains in protoplanetary disks for million-year timescales (Pinilla et al. 2012).
Observational surveys that measure the disk dust mass found a power-law relation between M dust and M that steepens with time (e.g., Pascucci et al. 2016). This relation was recovered by dust evolution models that include the growth, fragmentation, and drift of particles from Krijt et al. (2016). Drift is crucial for reproducing the steepness of the M dust − M relation because it is expected to reduce the dust mass with time and because drift is more effective around low-mass stars (Pinilla et al. 2013). However, in these models, the disks are short in solids compared with observations taken after 1 Myr of evolution. When considering disks with large inner cavities resolved at millimeterwavelengths (the so-called transition disks; TDs), Pinilla et al. (2018) find that this relation is much flatter, suggesting that the flatness is due to an effective reduction of radial drift in pressure bumps.
In this paper, we aim to study the effect that dust evolution has on the M dust − M relation. We are particularly interested in investigating the effect of particle traps in the disk on this relation. Motivated by the fact that exoplanet statistics show a dependence on stellar mass and metallicity, we focus also on the effect of having different stellar and disk masses in the dust evolution models, assuming that the disk dust mass is a proxy of stellar metallicity. Furthermore, we investigate the effect of having pressure bumps of different amplitudes. In the context of planets being responsible for those pressure bumps, it is an open question of whether the M dust − M relations are, perhaps, tracing the presence of giant planets (hence stronger pressure bumps) around more massive stars, while sub-Neptune planets (hence weaker pressure bumps) around low-mass stars, in agreement with exoplanet statistics. In this pic-ture, planets have already formed at the very first million years and the millimeter disk mass is the leftover mass, which does not drift inward and can be still observed by ALMA.
This paper is organized as follows. In Sect. 2, we explain the dust evolution models and the assumptions made when comparing them with observations. In Sect. 3, we present the results of these dust evolution models and the comparison with current millimeter observations of planet-forming disks, in particular in the context of the M dust − M relation. In Sect. 4, we discuss our results and assumptions, and the main conclusions are given in Sect. 5.

Dust evolution models
To study the effect of pressure traps in the observed M dust − M relation, we performed dust evolution models that include the growth and dynamics of particles, based on the models from Birnstiel et al. (2010b).
In these models, all the particles are initially micronsized grains distributed as the gas with a dust-to-gas ratio of 1/100. We assumed a maximum grain size of 2 m, and the grain size grid was logarithmically spaced in 180 cells. For the growth process, we took into account the coagulation, fragmentation, and erosion of particles; while for the dynamics, particles are influenced by the gas drag and, hence, the radial and azimuthal drift, the turbulent mixing, vertical settling, and Brownian motion. Particles stick and grow when their relative velocities before collisions are below a threshold, which is set in the simulations to be 10 m s −1 . This value is based on numerical and laboratory experiments of icy particles, that show that above these velocities particles are expected to fragment (e.g., Blum & Wurm 2000;Blum 2018;Paszun & Dominik 2009;Gundlach & Blum 2015), although recent laboratory experiments show that ice particles may be as weak as silicate particles (e.g., Gundlach et al. 2018), with fragmentation velocities of ∼1 m s −1 . The dust diffusion is assumed to be as the turbulent gas viscosity (Youdin & Lithwick 2007) and the turbulent velocities are proportional to √ α, which is the effective viscosity parameter for the disk evolution (Shakura & Sunyaev 1973). In the models, this parameter α has an influence, therefore, on the relative velocity of particles (and their fragmentation) due to turbulence and settling, as well as on the diffusion or mixing of grains.
For the models, we assumed different stellar parameters and disk masses and fixed other disk parameters. This approach is aimed at achieving a better understanding of the effect of stellar and disk mass on the radial drift and trapping of particles and, therefore, on the final M dust that would be observed. Because radial drift is more efficient around low-mass stars (Pinilla et al. 2013;Zhu et al. 2018), we expect to obtain dust masses that decrease with decreasing stellar mass. To visualize this, we consider four different type of stars: (1) a very low-mass star with a mass of M = 0.1 M , a luminosity of L = 8 × 10 −2 L , and an effective temperature of 3,000 K. These properties are motivated by the CIDA 1, a low-mass star whose disk has a large inner cavity resolved at millimeter wavelengths ); (2) an intermediate-mass star with M = 0.3 M , with a luminosity of L = 0.2 L , and an effective temperature of 3,700 Kl (3) a solar-type star, that is, with M = 1 M , L = 2 L , and an effective temperature of 4,700 K. And finally, (4) a Herbig star with M = 2.0 M , a luminosity of L = 20 L , and an effective temperature of 9000 K. In the first set of simulations, we keep the disk mass the same independently of the the stellar mass, such that M disk = 0.05 M . These models aimed to investigate the effect of radial drift around different stellar masses only. In the second set of simulations, we assumed M disk = 0.05 M . According to protostellar disks formed in radiation hydrodynamical simulations of star-cluster formation, the disk mass is expected to be around 0.01-0.08 M for ages higher than 10000 years (Bate 2018) The temperature of the dust is assumed to be a power law that depends on the stellar luminosity, such that, where T 10 is taken to be 30 K (Andrews et al. 2013;Tripathi et al. 2017) and L is the stellar luminosity. Equation 1 is expected if the dust disk size is independent of stellar mass (Hendler et al. 2017).
Since the disk temperature is a function of stellar luminosity, we expect significant differences in some physical quantities across the range of stellar masses that we study. First, at a given distance, the sound speed is higher around more massive stars, making the disk scale height also higher (h(r) = c s Ω −1 , with c s being the sound speed and Ω the Keplerian frequency). Second, in dust evolution models, the relative velocities between particles are dominated by turbulence and radial drift. A change in temperature affects the maximum value of turbulent velocities, which are directly proportional to the sound speed, v turb ∝ √ α turb c s (Ormel & Cuzzi 2007). The disk surface density distribution is parametrized by an exponentially tapered power-law function, given by where γ = 1, R c = 80 au, and Σ 0 is taken such that the disk mass is either 0.05 M or 0.05 M . The radial grid is taken from 1 to 300 au and it is logarithmically spaced in 300 cells.
Particle traps. To include the pressure bumps in the disk, we assume Gaussian perturbations to the disk surface density given by where A is the amplitude, r p the center, and w the width of the Gaussian perturbation. We assume three pressure bumps in the disk motivated by the average amount of substructures observed in protoplanetary disks at high angular resolution Long et al. 2018) and claimed to be particle traps . Hence, the gas surface density of the perturbed disk is given by where B 1 , B 2 , and B 3 are given by Eq. 3. We consider two different types of perturbations. One type we refer to as strong pressure bumps, with A = 4; and another that we call weak pressure bumps, with A = 1. This simple parametric form aims to mimic the presence of multiple giant planets or (sub-) Neptune planets creating pressure bumps exterior to their orbits. We note that, in reality, A depends not only on the planet mass, but also on several disk parameters, such as viscosity and local scale height (e.g., Crida et al. 2006;Fung et al. 2014). When comparing the assumed gas surface density with the hydrodynamical simulations from Zhang et al. (2018), who studied three different disk scale heights and three different values of α for 4 different planet masses, the amplitude of our perturbations resemble planet masses of 33 M Earth and 0.3 M Jup in the case of α = 10 −3 and h/r = 0.05 when A = 4 and A = 1, respectively; or 0.3 M Jup and 1 M Jup when h = 0.07 − 0.1 (for A = 4 and A = 1, respectively). For this α value and planet masses, the gap width does not change significantly with the planet mass, so our assumption of keeping the same width independently of A is fair for this range of planet masses and disk scale heights. In both cases, the Gaussian perturbations are centered at 20, 40, and 80 au; and the width is the local scale height h(r). Because we assumed different stellar masses and luminosities (and, hence, different disk temperatures), the local scale height is different for each case. However, we set the width of the perturbation as the disk scale height that corresponds to the case of a very low mass star; this means that  for all cases, w has values of ∼2, ∼5, and, ∼11 au at 20, 40, and 80 au, respectively. The width of the perturbation is, therefore, larger than the local scale height for the other three stellar masses, becoming almost a factor of two for the Herbig disk. Although the Herbig disk is warmer and, hence, it is expected that the disk scale height is higher at a given location, the dependence on the stellar mass makes that the disk scale height increases for the low mass stars. Hence, the widths are larger or equal to the pressure scale height for all cases to ensure the stability of the pressure bumps (e.g., Pinilla et al. 2012;Dullemond et al. 2018). The width of the perturbation should not be confused with the width of the dust rings (see e.g., DSHARP) as the latter can be much smaller than the width of the pressure bump and of the local scale height. For the two types of pressure bumps, the perturbed gas surface density is normalized such that the disk mass remains as the unperturbed cases, that is, 0.05 M or 0.05 M . Figure 1 shows the gas surface density profile for the case without any traps (unperturbed density) and the cases with strong and weak pressure traps when M disk =0.05 M . Table 1 summarizes the model parameters considered in this work.

Synthetic dust mass
In order to compare the results of the models with observations, we took the dust density distribution at different times of evolution and calculated the vertical optical depth (τ ν ) at 0.87 mm, which corresponds to ALMA Band 7, the band adopted for most of the medium-resolution surveys of star-forming regions.
To calculate τ ν , we took the vertically integrated dust density distribution σ(r, a) at a given time of evolution and calculate τ ν = σ(r, a)κ ν / cos i, where κ ν is the opacity for each grain size and at a given frequency or wavelength (κ ν ). To obtain κ ν , we use Mie theory and we assume the same volume fractions and optical properties for the dust particles as in Ricci et al. (2010). Furthermore, we consider face-on disks only, that is, i = 0. This assumption has no effects for optically thin discs, but it can lead to overall higher disc fluxes and thus higher predicted masses if discs are close to being optically thick. With the optical depth, the intensity profile at a given wavelength is given by We obtain the total flux as where d is the distance to the source, which we set to 140 pc, and B λ (T (r)) is the Planck function, for which we take the temperature profile of Eq. 1.  To obtain the dust mass as is usually calculated from the millimeter fluxes, we assume optically thin emission (Hildebrand 1983) and, hence, For κ ν , we assume what is usually taken in disk surveys, which is a frequency-dependent relation given by (Andrews et al. 2013), and for B ν we assume a blackbody surface brightness corresponding to 20 K.

Stellar mass effect
In this section, we first describe the results when M disk = 0.05 M for all cases and focus on the effect of having different stellar mass. The top panels of Fig. 2 show the dust density distribution after 1 Myr of evolution and α = 10 −3 for the cases of no traps (left panel), weak traps (middle panel), and strong traps (right panel) for a stellar mass of 0.1 M . This figure shows the effect of particle trapping in pressure bumps and the effect of their amplitude.
Unperturbed density. In the absence of particle traps, grains grow from micron-sized particles to the maximum grain size that is seen in the case of 0.1 M set by fragmentation in most of the disk. When particles reach centimeter or millimeter sizes, they drift towards the star, depleting the disk of pebbles in million-year timescales. This is shown in the left panel of Fig. 3, which illustrates the effect of stellar mass on the evolution of the total mass in particles larger than 0.1 mm (hereafter, called mass in large dust). In this case, the maximum mass in large dust is quickly reached after few thousand of years of evolution and decreases with time until few large dust particles remain in the disk. The maximum of the total mass in large dust is around 50 M ⊕ (30% of the initial mass in solids) for the Herbig star and it is around 25 M ⊕ for the low-mass star, even when in all cases, the disk (and dust) mass was initially the same. This difference is a result of the growth and more efficient drift around low-mass stars. However, after ∼1 Myr of evolution, the mass in large dust is almost the same for the four stellar masses (∼ 20 M ⊕ ). At 5 Myr of evolution, the dust mass has decreased to less than 1 M ⊕ . Particle traps. When particle traps are included in the disk, particles can grow to millimeter or centimeter sizes and remain trapped in pressure maxima (middle and right panels of Fig. 2 and Fig. 3). In the case of weak traps (with A=1.0 in Eq. 3), the amplitude of the pressure bumps is enough to keep a high amount of large dust in the disk for the cases of a Solar-type star and a Herbig star (2 M ). The mass in large dust is around 130 M ⊕ and 110 M ⊕ for the 2 M stellar mass at 1 and 5 Myr of evolution, respectively (Fig. 3). For the Solar type star the mass in large dust is around 85 M ⊕ and 70 M ⊕ at 1 and 5 Myr of evolution, respectively. For the 0.1 M and 0.3 M stars, these pressure bumps are not enough to retain millimeter-or centimeter-sized particles and after around 0.5 Myr of evolution the mass in large dust decreases with time reaching values below 2 M ⊕ at 5 Myr of evolution. This is because of the more efficient radial drift around low-mass stars (Pinilla et al. 2013). The radial drift of particles is a consequence of the difference between the gas azimuthal velocity and the Keplerian speed, which is higher for disks around low-mass stars. Therefore, the weak traps cannot efficiently trap particles for the two low-mass stars that we consider (middle panel of Fig. 3).
When strong pressure bumps are assumed in the dust evolution models, their amplitude is high enough to retain most of the large dust independent of the stellar mass at million-year timescales (right panels of Figs. 2 and 3). From Due to the accumulation of dust particles in these strong pressure maxima, the disk emission becomes optically thick in several disk locations. In addition, dust grows as large as the maximum grain size that we allow in the simulations (right panels of Fig. 2). Both effects result on an underestimation of the dust disk masses when they are calculated from sub-millimeter fluxes, as discussed in Sect. 3.3.
Effect of the stellar mass on the fragmentation barrier. The maximum size that particles can reach before they fragment, or the fragmentation barrier, is given by ): When the disk surface density and α−viscosity are the same in the simulations, the effect of changing the stellar mass and luminosity is reflected in the sound speed. For a more luminous star, that is, with a warmer disk, the sound speed is higher, increasing the maximum turbulent velocity and therefore decreasing a frag . As a result, in a warmer disk, particles reach lower sizes. Figure 4 shows the case of strong pressure bumps in the case of a disk around a Solar type star. The maximum grain size that particles can reach, is lower than in the case of a 0.1 M stellar mass (top right panel of Fig. 2), an effect that is more clear inside the pressure bumps. In the case of a Solar-type star, more millimeter or centimeter grains remain in the disk at millionyear timescales because fragmentation keeps the particles smaller and closer to the sizes where they are more affected by the gas drag, which means that they are easier to trap in pressure bumps. This effect will influence the synthetic sub-millimeter fluxes that are derived from the models because when growth is very efficient (as in the case of strong bumps and low-mass stars), the very large grains ( 10 cm) that are formed inside the pressure maxima have very low opacities contributing very little to the millimeter fluxes.

Disk mass effect
There are two main effects of changing the disk mass to be a fraction of the stellar mass (0.05 M ) instead of assuming 0.05 M for all the cases. First of all, the fragmentation barrier (a frag , Eq. 8) is directly proportional to the gas surface density (Σ), implying that a frag is lower in low-mass disks even though the sounds speeds are lower. On the other hand, the radial drift of particles comes from the gas drag forces, which depend on dust properties such as their mass and cross-sectional area, but also on the gas disk via the gas density. The dimensionless Stokes number (St) measures the importance of drag force by comparing the stopping time of the particles to the Keplerian frequency at a given radius. Particles with the Stokes number equal unity are subject to the highest radial drift. At the midplane, the Stokes number is inversely proportional to the gas surface density and as a consequence, a millimeter-sized particle in the outer disk has a higher Stokes number in a lower mass disk. Therefore, the drift barrier is also lower in lower mass disks, which means that particles grow to smaller sizes in a low-mass disk, but also a particle with a given size drifts faster in a low-mass disk. In order to compare the results from dust evolution models with current millimeter observations of protoplanetary disks, in particular in the context of the M dust − M relation, we used the values for the slope and intercept for different star-forming regions from Pascucci et al. (2016) as stellar masses were homogeneously derived in this study, except for σ-Orionis, which was not included; hence, for this parameter, we took the values from Ansdell et al. (2017). In all regions, T = 20 K was assumed (see solid lines in different colors in Fig. 5). We also collected from the literature a sample of disks that have been observed with ALMA, which have either a large inner cavity or gap, defined as transition disks (TDs, empty circles), and disks with substructures, such as rings, gaps, or asymmetries, but still exhibiting substantial emission from the inner disk (full triangles). We selected the TDs that have been observed with ALMA and whose cavities are resolved, which implies that the TDs in our sample have large cavities ( 20 au). We gathered a total of 73 disks (43 TDs and 30 disks with substructures, see Table in the appendix). The sample of TDs double the number of targets presented previously in Pinilla et al. (2018). We placed the TDs and disks with substructures in the M dust − M relation and fit these points with the following relation: log 10 (M dust /M ⊕ ) = β log 10 (M /M ) + α. We obtained β = 0.72 +0.04 −0.04 and α = 1.86 +0.01 −0.01 , and a Pearson co-efficient of 0.50. The fit takes into account the uncertainties of the data, dominated by the 10% of uncertainty from flux calibration. When only the TDs are considered for the fit, the Pearson coefficient is higher (0.67) and, therefore, this correlation is stronger when only TDs are included. This is because for this correlation, we mainly have data for TDs around low-mass stars. If only disks with substructures are considered, the correlation is weak with a Pearson coefficient of 0.28. However, more observations in the low stellar mass regime that include the resolution of substructures are needed to confirm if there is indeed a correlation. The values for the slope and intercept change when including or excluding disks with substructures. When we only consider TDs, the values remain similar with β TD = 0.94 +0.06 −0.06 and α TD = 1.80 +0.02 −0.02 , but the relation is shallower when only disks with substructures are considered, with β SB = 0.42 +0.1 −0.1 and α SB = 1.90 +0.02 −0.02 . With the current data, there is not any substantial evidence that the population of TDs is different than the population of disks with substructures based on the dust disk mass. To compare the two samples, we calculate the Kolmogorov-Smirnov (KS) two-sided test, and we find a maximum deviation between the cumulative distributions of only 0.15 with a high probability (∼80%) that the two samples are similar. Overall, the M dust − M correlation of the whole sample is dominated by TDs, and it is shallower than seen from all disks in nearby star-forming as found in Pinilla et al. (2018), to which some of these disks actually belong to.
Next, we overplotted the M dust obtained from the models at 1 and 5 Myr of evolution with the observed correlations. We plotted M dust that were directly obtained from the dust evolution models for large dust without any postprocessing assumption (as e.g., the values in Fig. 3) and, after dust temperature and opacities were assumed, we calculated the dust mass expected from millimeter observations as explained in Sect. 2.2 (hereafter called synthetic dust mass). Figure 6 shows the results for all the models that assumed M disk = 0.05 M and Fig. 7 shows the results for the cases where M disk = 0.05 M .

Synthetic dust mass compared to mass in large dust from models
In general, the synthetic dust mass is lower than the mass in large dust that is obtained directly from the models (Fig. 6 and Fig. 7). This difference is because of two different reasons. First, some parts of the disk become optically thick.
In the case of an unperturbed density profile, the optical depth of the disk (τ ) reaches values of τ ∼ 0.1 in the case of low-mass stars (< M ), and τ ∼ 0.5 in the case of higher stellar masses at 1 Myr. The higher τ around more massive stars, even when the disk mass is assumed to be the same independently of the stellar mass, is due to the drift more efficient removal of large particles in the case of low-mass stars. At 5 Myr of evolution, the optical depth is much lower than unity in all the cases without particle traps, while the synthetic dust mass and the mass in large dust from the models are fairly in agreement. In the case of particle traps, the optical depth can reach values of between 1 to 10 inside the bumps. For the weak pressure bumps, these high values for τ are reached only in the innermost bump, while in the case of strong bumps, these values are reached at the peak of all the bumps. In the cases of particle traps, a second reason plays an important role to understand why the synthetic dust mass is lower than the mass in large dust from models, and it is due to the growth of solids beyond 10cm, hereafter referred to as boulders (in the Minimum Mass Solar Nebula (MMSN) model, the (sub-)meter-sized objects that are in between being strongly influenced by the gas and particles that are completely decoupled and moving with Keplerian speed are called "boulders".) Boulders have very low opacities and their contribution to the millimeter fluxes is insignificant, which means that their mass is hidden when calculated from the millimeterfluxes. This effect is clear when comparing the masses in the case of strong bumps. For low-mass stars (0.1 or 0.3 M ), boulders are formed inside the bumps as discussed in the previous section (Fig. 2), while in the case of disks around 1 or 2 M , the fragmentation barrier keeps the maximum grain size around millimeter and centimeter sizes in all the pressure bumps (Fig. 4). The formation of boulders reduces the optical depth inside the pressure bumps, but simultaneously, a lot of the solid material becomes invisible when observed at millimeter-wavelengths. As a consequence, while the synthetic dust mass is lower by a factor of ∼ 1.5 − 2 for the disks around stars M ≥ M , the synthetic dust mass is a factor ∼ 10 − 100 lower than the total mass in large dust for M < M . Hence, in the case of traps, the difference between the synthetic dust mass and the mass from the models originates from a combination of these two effects: optical thickness and boulder formation. For the weak traps, these two effects appear only in the innermost bump at 1 Myr of evolution. While for strong traps, these two effects appear in all the pressure traps at ≥ 1 Myr.
In summary, both optically thick regions and boulder formation lead to a large difference between the synthetic dust mass and the mass in large dust obtained directly from the models. Our result show that boulder formation has a more significant influence. Recently, Stammler et al. (2019) showed that to explain the optical depth between 0.2 and 0.5 inside the rings observed in the DSHARP sample (Andrews et al. 2018b), boulder and planetesimal formation must occur inside the pressure bump to reduce the optical depth. In our models, we do not reduce the mass of the dust inside the pressure bumps by assuming that boulders were formed when high dust-to-gas density ratio is reached, but rather particles naturally grow to sizes beyond decimeter when drift and turbulence are inefficient to fragment or keep the particles with millimeter to centimeter size.

Synthetic dust mass versus observations:
M disk = 0.05 M case In the less realistic case of disk masses independent of stellar masses, the observed slope in different star forming regions is not reproduced by the models, not when comparing directly with the models nor after post-processing to compare with the observations (see Fig. 6). Only in the case of weak traps, the slope of the sample of TDs and disks with substructures is reproduced when comparing the mass directly from the models at 1 Myr of evolution. When comparing it with the synthetic dust mass, the slope of the TDs and disks with substructures is similar to the models with and without traps, but values are lower than observed, meaning that the initial dust mass should be higher for all stellar masses to have an agreement with the relation of TDs and substructured disks, an assumption that may be unlikely since with the assumed disk mass of 0.05 M , the stellardisk mass ratio is already very high for the disks around 0.1 M (stellar-disk mass ratio of 0.5) and 0.3 M (stellardisk mass ratio of 0.17). Strong pressure bumps are needed to produce the high disk dust masses that are observed around stars with masses equal or higher than the Solar mass (when comparing with the synthetic dust mass). With weak or no pressure bumps, the millimeter fluxes and hence the dust disk masses are lower than observed already at 1 Myr of dust evolution. In the context of planets causing these pressure bumps, this implies that if the disk mass is initially the same for all disks around different stellar masses, only the strong pressure bumps (or giant planets) can explain the observed values of M dust for stellar masses higher than a solar mass. For low-mass stars no conclusion can be given besides that in all scenarios the synthetic dust mass is overestimated at 1 Myr. At 5Myr, the synthetic dust mass values are in agreement with the values in Lupus, Taurus and ChaI, although the averaged age of these regions is lower than 5 Myr. at 1 Myr when compared directly with the mass in large dust obtained from the models (see Fig. 7). When comparing with the synthetic dust mass predicted by our model, the slope at 1 Myr is slightly lower than in the actual observations. We find M dust values that are too low for the case of 1 and 2 M at 5 Myr (both dust mass: directly from models or the synthetic dust mass as inferred from the simulations).
In the case of weak traps, we obtained two interesting results. First, when comparing it with the mass in large dust directly obtained from the models, the M dust − M relation of TDs and disks with sub-structures is reproduced at 1 Myr of evolution. The sample of TDs and disks with substructures (Table A.1 in the appendix) is a combination of disks at different stages of evolution, from 0.5 to 10 Myr, but most of them are older than 1 Myr (Garufi et al. 2018). Second, due to the reasons explained above (optical thickness and boulder formation), when comparing it with the synthetic dust mass, we find that the model becomes steeper and closer to the slope of different star forming regions, especially at 1 Myr. At 5 Myr, the synthetic dust mass is lower than in any star forming region, in particular for disks around 1 and 2 M stellar masses.
The case of strong traps shows very interesting results but they must be interpreted with caution. First, the M dust − M relation observed in TDs and disks with substructures is fairly reproduced at any time of evolution when comparing with the mass in large dust obtained directly from the models. However, due to boulder formation, in particular for the 0.1 and 0.3 M , the synthetic dust mass is much lower reproducing again the steeper trends at any time of evolution. This implies that to explain the trend of TDs and disks with substructures, strong pressure bumps (perhaps due to giant planets) are needed, but in addition, boulder formation must be inhibited inside these pressure bumps, otherwise the observed fluxes and masses are too low, reproducing instead the steeper relations in different star-forming regions.
Therefore, in the context of planets causing these pressure bumps, we found similar results as in the case of M disk = 0.05 M for the disks around 1 and 2 M , meaning those for which only strong pressure bumps (or giant planets) can explain the observed values of dust mass at any time of evolution. Contrary to the case of M disk = 0.05 M , any of the three scenarios (no traps, weak or strong traps) can reproduce the M dust − M relations of different star formation regions for 0.1 and 0.3 M at 1 or 5 Myr of evolution. The agreement with the strong traps scenario is surprising, but it is due to boulder formation, while in the other two cases, it is due to the actual loss of particles attributed to the drift. forming regions (bottom right panel of Fig. 7) because of the trapping of dust particles and boulder formation inside pressure bumps, the latter in particular for M < 1 M . In contrast, the flatter relation observed for TDs and disks with substructures has not been reproduced by any model, specifically one based on post-processing for the purpose of comparison with the observations. In the framework of our models, one solution to reproduce the relation of TDs and disks with substructures is to consider that boulder formation is inhibited inside pressure bumps in the case of low-mass stars (M < 1 M ).
Our results seem to contradict those of Stammler et al. (2019), who suggests that dust must be transformed into planetesimals to obtain the optical depths values inferred inside the DSHARP rings (Andrews et al. 2018b), using dust evolution models that include trapping in a pressure bump. However, that work was done assuming a Herbig star, in particular to reproduce the observations of one of the brightest disk (HD 163296), while our conclusion about hindering boulder formation inside pressure bumps applies only to low-mass stars (M < 1 M ). In our case of a disk around a Herbig star, we do not have any natural production of boulders in pressure maxima and the small differences (a factor of ∼ 1.5 − 2) between the mass in large dust directly from the models and the observed total dust mass originate from the high optical depth (∼1-2) inside the pressure maxima. According to our results, inducing boulder formation in the pressure bumps in the Herbig case can help to reduce τ and be in agreement with the values obtained in the DSHARP sample. However, as found in our models, boulder formation reduces by a lager factor the synthetic dust mass (∼ 10 − 100, comparison between right panels of Fig. 7). Assuming a higher initial disk mass may help to have a better agreement with the synthetic dust mass as in Stammler et al. (2019), where the initial condition for the disk mass is 0.4 M (four times than what we assume in this work) and they do reproduce the intensity of one of the observed rings of HD 163296.
There are several possibilities that may hinder boulder formation in pressure maxima in disks around lowmass stars. The first possibility is that the velocity threshold for which particles fragment (fragmentation velocity) is lower than what we considered (10 ms −1 ). Laboratory experiments and numerical simulations show that these values are possible (e.g., Paszun & Dominik 2009;Gundlach & Blum 2015;Musiolik et al. 2016) for water ice grains, although recent experiments suggest that the fragmentation velocities for water ice grains may be as low or lower as for silicates, that is, 1 ms −1 (e.g., Musiolik & Wurm 2019). However, if this is the case, the fragmentation barrier would be two orders of magnitude lower (Eq. 8), which will imply that the maximum grain size would be set by fragmentation with typical values of ∼100µm in the outer disk (>20 au). These grain sizes are in agreement with the values proposed by Liu (2019) and Zhu et al. (2019), who suggest that the measured optical depth in several DSHARP disks can be explained by dust scattering when the optically thick dust has an albedo of ∼0.9 at 1.25 mm, which corresponds to a maximum grain size of ∼0.1-1 mm. Nevertheless, if this is the case, dust trapping may be inefficient in reproducing observations. Investigating the required conditions for particle trapping, and an agreement with observations when the fragmentation velocities are 1 ms −1 in the entire disk, is a subject of future research.
Other possibilities that may explain the inhibition of the formation of boulders include bouncing in the models or an increase in the disk α viscosity. Nevertheless, to include bouncing would keep the grains too small (Stokes numbers lower than 10 −4 , e.g., Zsom et al. 2010) such that they could not effectively be trapped in pressure maxima. The other alternative is to have a higher value of the α viscosity, which sets the fragmentation barrier and turbulence mixing strength; hence, determining how much grains can be diffused out from pressure bumps, which is discussed in the next subsection.

Effect of turbulence mixing strength
As discussed in the previous subsection, increasing the turbulence mixing strength may help to hinder boulder formation inside pressure bumps. In general, changing the value of the α viscosity in our models has a strong effect in the results. This effect has been studied and discussed in detail in de Juan Ovelar et al. (2016) and Bae et al. (2018), who show that by just varying the values of α and keeping the rest of the model parameters the same, the synthetic images from models can have very different outputs, such as: multiple rings and gaps, a cavity with a single ring, or compact disks.
In the case of planets, usually a higher mass planet is required to open a gap when the disk viscosity is higher. This balance between planet mass and viscosity has an effect in the potential trapping of particles ( Fig. 1 in de Juan Ovelar et al. 2016). The amplitude of A = 4 of our perturbation in Eq. 3 corresponds to planet masses similar to 0.1-1 M Jup for 10 −3 − 10 −2 α−viscosity. We note that increasing the amplitude of our perturbation, which would resemble gaps formed by more massive planets or less viscous disks, does not have a strong effect in our models, since with A = 4 effective trapping is already reached and therefore increasing A would only saturate the amount of dust trapped in pressure maxima. In our assumptions of the particle traps, the parametric perturbation does not depend on α, and hence the shape of the pressure bumps is independent of gas viscosity. As a consequence, the only effect of varying α in our models is reflected in the dust evolution.
Varying α directly influences the fragmentation barrier a frag ∝ 1/α, which means that for higher α, grains reach lower sizes when the maximum grain size is set by fragmentation. In addition, we assume dust diffusion D as Youdin & Lithwick (2007), which means D ∝ α/(1 + St 2 ) and, thus, higher dust diffusion for higher values of α. Increasing α implies lower a frag and higher D, which influence the potential trapping because if grains do not grow to values that attain a high Stokes number, they are more difficult to trap and they diffuse out of pressure traps easily. To see the effect of α on the mass in large dust, we simulated the case of M = 0.1 M and M disk = 0.05 M with three different values of α (10 −4 , 10 −3 , 10 −2 ). Figure 8 shows the main results of this numerical experiment.
When the disk viscosity is high (α = 10 −2 ) trapping is less efficient and it does not occur even for the case of strong traps (Fig. 8). In these cases, the peak of the mass in large dust is reached earlier in the simulations (∼ 10 5 years). On the other hand, when disk viscosity is low (α = 10 −4 ), the peak of the mass in large dust is reached later (∼1 Myr), and in this case trapping is very efficient even for weak traps. In this case, the fragmentation barrier is one order of magnitude higher than our standard models, which implies more efficient formation of boulders inside pressure bumps, creating more invisible solid material when observed at millimeter-emission and a higher discrepancy with the observed dust mass. Tripathi et al. (2017) and Andrews et al. (2018a) find a significant correlation between the disk size (R eff , usually defined as the radii that encloses 68% of the millimeter emission) and the continuum luminosity of the disk. Hendler et al. (2020) performed a homogeneous analysis of five star forming regions, finding that this relation flattens for Upper Sco, which is the oldest region in the sample (left panel of Fig. 9). The observed relation could originate from the growth and radial drift of dust particles, where the dust disk extension becomes smaller with time due to drift. Depending on the initial conditions, dust evolution models can re-produce the observed relations (Tripathi et al. 2017;Rosotti et al. 2019). It is, therefore, possible that the R eff − L mm flattens for disks in which radial drift is reduced or ineffective.

Disk size-luminosity relation
We checked the R eff − L mm relation for our sample of TDs and disks with substructures and found a significan correlation (Pearson coefficient of 0.57) with an intercept (A) equal to 2.29±0.02 and a slope (B) equal to 0.36±0.03, assuming log 10 (R eff /au) = B log 10 (L mm /mJy) + A. For this calculation, the millimeter luminosity is scaled to a 140 pc distance for all the disks. The values for R eff are the ones that enclose 90 or 95% of the millimeter flux in Band 6 or Band 7 with ALMA (the wavelength is between ∼ 0.8 and ∼1.3 mm). The relations found in Hendler et al. (2020), Tripathi et al. (2017) and Andrews et al. (2018a) used the radius that encloses the 68% of the flux, but Hendler et al. (2020) demonstrates that there is a linear relation close to 1-1 between the two values.
The values for R eff of the sample are larger than the averaged values found in different star forming regions, with a mean value of 116 au (see Table 2 from Hendler et al. (2020) for comparison). This is because most of the ALMA observations of TDs and disks with substructures are biased towards disks with larger radii. The slope of R eff − L mm in our sample is smaller than the values found in the youngest star forming regions studied in Hendler et al. (2020), that is, Ophiucus, Lupus, and Taurus/Auriga with ages between 1 and 3 Myr. The slope of our sample lies in between the values found for Chamaeleon I (2-3 Myr, slope of 0.4) and Upper-Sco (5-1 Myr, slope of 0.22), while the ages of the TDs and disks with substructures in our sample range between 0.1 Myr and 10 Myr.
Based on our dust evolution models, we calculate the radius that encloses 95% of the mass of the dust particles that are larger than 0.1 mm. The evolution of this radius is shown in the right panel of Fig. 9 for the models that assume strong pressure bumps, M disk = 0.05 M , and α = 10 −3 . Independent of the stellar mass, after ∼0.5 Myr of evolution the value of this radius converges to the location of the farther pressure maximum, that is, ∼ 75 au (similar values are found when the disk radius is calculated from the synthetic intensity profile). Based on these models, we expect a flat relation between R eff and L mm at any evolutionary stage. It is possible that the observed R eff − L mm relation for TDs and disks with substructures provides information about how far pressure bumps may form and how the less luminous disks that are around stars of lower mass may have the pressure bumps closer in. Exoplanet statistics suggest that low-mass stars have more planets with radii of 1-4 R ⊕ at the same semi-major axis (e.g., Mulders et al. 2015) but to test this idea, we require more detailed models to investigate different locations and the morphology of pressure bumps, along with a greater number of high-resolution observations around low-mass stars.

Conclusions
We investigated the effect of dust evolution and particle trapping in the observed M dust − M relation of different star-forming regions. We focused our attention on the effect of the stellar and disk mass. Both properties affect the mass in large dust at million-year timescales due to: (1) the more effective drift around low-mass stars; (2) the more efficient fragmentation for disks around stars of higher mass -which  are less efficient in more massive disks; and (3) the drift being more efficient in disks of lower mass. The combination of these effects influence the efficiency of trapping and the potential to form boulders in pressure bumps, which influence the final M dust . Our conclusions are: 1. Independent of the assumption for the disk mass, either M disk = 0.05 M or M disk = 0.05 M , strong traps are required to reproduce the observed M dust of disks around M ≥ 1 M at different time of evolution. In the context of planets as the origin of pressure bumps, this result agrees with current exoplanet statistics about giant planets, which are found to be more common around more massive stars. 2. When all disks have the same initial mass independent of the stellar mass (M disk = 0.05 M ), the dust evolution models produce a flatter M dust − M relation when compared to the observed ones from different star-forming regions. Therefore, we conclude that the initial disk mass has to scale with the stellar mass to reproduce the slope of the M dust − M relation, otherwise it is flatter than what is observed. In the case of M disk = 0.05 M for all cases, the slope obtained from the models is similar to the one observed for TDs and disks with substructures, however the modeled M dust is overall lower than the observed one. 3. When the disk mass is a fraction of the stellar mass (5%), our models reproduce the observed M dust − M relation from different star-forming regions. However, in the cases of an unperturbed density and weak traps, the trends are reproduced only at 1 Myr while at later times (5 Myr), M dust is too low when compared to observations for M ≥ 1M

5.
To explain the flatter relation of M dust −M for TDs and disks with substructures, strong pressure bumps (perhaps caused by giant planets) are needed. However, efficient growth must be inhibited inside these pressure bumps for low-mass stars, otherwise the observed fluxes and masses are too low, reproducing instead the steeper trends in different star-forming regions. Different possibilities can hinder very efficient growth inside pressure bumps, such as a lower fragmentation velocity or including the effect of bouncing. 6. Due to the efficient formation of decimenter bodies in pressure bumps for disks around low-mass stars M < 1 M , we cannot give a definitive conclusion on what type of traps are present in these disks. Indeed, all three cases we tested (unperturbed, weak, and strong traps) can reproduced the observed values in the more realistic case of disk mass scaling with stellar mass. 7. The R eff − L mm relation for TDs and disks with substructures is flatter than observed in Ophiucus, Lupus, and Taurus/Auriga star-forming regions and the slope value of our sample lies in the range found for older regions, that is, Chamaeleon I and Upper-Sco. This relation may flatten due to inefficient radial drift when pressure bumps exist in the disks, in which case R eff traces the location of the farther pressure bump at any time of evolution. Currently, high angular resolution observations of TDs and disks with sub-structures are biased toward the brightest and larger disks, so more observations of small and faint ones are required to test this idea.
Future multi-wavelength and high-resolution observations of disks around low-mass stars are needed to investigate if gaps or rings are present and how much grains have growth there in comparison to the rings and gaps in disks around more massive stars. Our results suggest that boulder formation must be inhibited inside pressure bumps in TDs and disks with sub-structures in order to explain the observed M dust . Further theoretical research is needed to investigate what conditions are required to explain the observed properties of this group of disks, in particular if boulder formation is hindered when the fragmentation velocities are lower than assumed in this work -as has been suggested by recent laboratory experiments.