Growing and Trapping Pebbles with Fragile Collisions of Particles in Protoplanetary Disks

[abridged] Recent laboratory experiments indicate that destructive collisions of icy dust particles occur with much lower velocities than previously thought. When these new velocities are considered from laboratory experiments in dust evolution models, a growth to pebble sizes in protoplanetary disks (PPDs) is difficult. This may contradict (sub-)mm observations and challenge the formation of planetesimals and planets. We investigate the conditions that are required in dust evolution models for growing and trapping pebbles in PPDs when the fragmentation speed is 1ms$^{-1}$ in the entire disk. We distinguish the parameters controlling the effects of turbulent velocities, vertical stirring, radial diffusion, and gas viscous evolution, always assuming that particles cannot diffuse faster (radially or vertically) than the gas. To form pebbles and produce effective particle trapping, the parameter that controls the particle turbulent velocities must be small ($\delta_t\lesssim10^{-4}$). In these cases, the vertical settling can limit the formation of pebbles, which also prevents particle trapping. Therefore the parameter that sets the vertical settling of the grains must be $\delta_z<10^{-3}$. Our results suggest that different combinations of the particle and gas diffusion parameters can lead to a large diversity of millimeter fluxes and dust-disk radii. When pebble formation occurs and trapping is efficient, gaps and rings have higher contrast at mm-emission than in the NIR. In the case of inefficient trapping, structures are also formed at the two wavelengths, producing deeper and wider gaps in the NIR. Our results highlight the importance of obtaining observational constraints of gas and particle diffusion parameters and the properties of gaps at short and long wavelengths to better understand basic features of PPDs and the origin of the structures that are observed in these objects.


Introduction
Planet formation occurs in protoplanetary disks, which are mainly composed by gas. Most of the information that we have about planets forming in disks comes from the dust that dominates the disk opacity. The initial properties of the dust in disks are inherited from the interstellar medium, such as its size and composition. These initial properties are expected to change by different physical processes that lead to collisions and transport of the particles (see, e.g., Birnstiel et al. 2016, for a review).
The growth of particles from (sub-)micron-sized dust to pebbles is fundamental for understanding the first steps of planet formation. From the observational point of view, we can have access to this process by observing the distribution of micron-sized dust to pebbles (submillimeter-to centimeter-sized particles) with observations from the optical to centimeter wavelengths, respectively. Solid bodies larger than centimeter-sized particles have very low opacities, and we are blind from the observational point of view about what occurs from the formation of pebbles to planets.
To connect current theories of planet formation with observations, we rely on connecting our observational knowledge of the small bodies (up to pebbles) to exoplanet statistics.
Recent deep observations of protoplanetary disks with high angular resolution have shown structures in these objects, such as rings, gaps, spiral arms, and lopsided shapes (e.g., Andrews 2020). These structures are believed to be regions of high pressure where the radial drift of pebbles is reduced or stopped completely (Whipple 1972), and where planetesimal formation can occur (e.g., Pinilla et al. 2012Pinilla et al. , 2020Dullemond et al. 2018;Lenz et al. 2019;Stammler et al. 2019;Klahr & Schreiber 2020). The current understanding of these structures from the observational and theoretical point of view is driving the field of planet-forming disks.
Another crucial part that helps to understand the first growth from interstellar medium dust to pebbles are laboratory experiments of particle collisions. They provide the potential outputs after collisions with a given initial relative velocity. Potential outputs include sticking, fragmen- In the standard models, they are all set to a value that depends on the α parameter that regulates the viscous evolution of the gas. We assumed all four parameters to be independent. tation, bouncing, or mass transfer (e.g., Windmark et al. 2012). They also define the limits for theoreticians of velocity thresholds for effective growth or destructive collisions (e.g., Blum & Wurm 2000Kimura et al. 2015Kimura et al. , 2020. These experiments are challenging, however, because the conditions must be set to be as close as possible to the physical properties of disks, which means low densities, pressure, and temperatures. It has commonly been thought that water-ice particles are more sticky than silicates (Blum & Wurm 2000;Wada et al. 2009Wada et al. , 2011Gundlach & Blum 2015;Musiolik et al. 2016), and the comparison has been made that this would be like snowballs versus sand. The difference in the fragmentation velocity between these two species is at least one order of magnitude. High fragmentation velocities for icy particles helps greatly in forming pebbles in the disk because the maximum grain size of pebbles in the fragmentation-dominated regime is proportional to the square of the threshold of the fragmentation velocity (e.g., Brauer et al. 2008;Birnstiel et al. 2010). Nonetheless, recent laboratory experiments have reproduced the temperatures of protoplanetary disks better, which changed our perspective on collisional growth because it was concluded that amorphous water-ice particles are as fragile as silicates or even more so (Gundlach et al. 2018;Musiolik & Wurm 2019;Steinpilz et al. 2019). The purpose of these experiments is to measure the sticking force when two particles are in contact. This force depends on the surface energy of the grains and therefore on the particle composition (e.g., Kimura et al. 2015). Most of the classical experiments used monolayer dust particles with an insignificant amount of water or no water at all. This limitation comes from keeping low values of temperature and pressure, and most of the information that we previously had about collision of (water) ices came from numerical experiments. More recent laboratory experiments were able to retain low temperatures in their collisions (<180 K) and showed that the surface energy of these grains decreases. This contradicts previous assumptions of collisions of water-ice particles.
The results of these new experiments suggest that the fragmentation velocities for icy particles is 1 m s −1 or even lower. This poses a problem not only for the formation of pebbles and hence for the formation of planets as a whole, but also makes the current interpretation of disk observations challenging. If the current structures are the byproduct of pressure bumps trapping the particles, grains have to grow at least to submillimeter sizes to efficiently slow down their drift at the pressure bump and be trapped. If particles remain small (typically 0.1 mm), radial drift is slow, and other effects such as turbulent mixing and radial diffusion may become more important, preventing the grains from being trapped. The maximum grain size that particles may reach with fragmentation velocities of 1 m s −1 may agree with the maximum grain size that Liu (2019) and Zhu et al. (2019) proposed to explain the optical depth of several disks recently observed with the Atacama Large Millimeter/submillimeter Array (ALMA) under self-scattering.
In this paper, we investigate the conditions that are required in dust evolution models to grow and trap pebbles in protoplanetary disks when the fragmentation speed is 1 m s −1 . The main questions that we consider are what values of turbulence, settling, and diffusion are needed to grow dust to pebbles and have particle trapping if the fragmentation velocities are as suggested from laboratory experiments. We define a reasonable value for the parameter that controls the viscous evolution of the gas versus dust evolution when we compare this with the observed properties of protoplanetary disks. To answer these questions, we investigate parameters that rule the turbulent velocities (δ t ), vertical stirring (δ z ), radial diffusion (δ z ), and viscous gas evolution (α), and consider that they are independent of each other. We distinguish the effect of each of these physical processes by always assuming that particles should not diffuse faster than the gas (i.e., δ r,z,t ≤ α). To test our models, we combine the results with radiative transfer simulations to compare them with observations at near-infrared and millimeter emission. This paper is organized as follows. In the next section, we briefly describe the models and the parameters that we assumed for the angular momentum transport of gas and dust (vertical and radial) diffusion. In this section we also describe our assumptions for the radiative transfer calculations. In Sect. 3 we present the results of our models. In Sect. 4 we discuss our results in the context of different observational diagnostics, such as millimeter fluxes, dust-disk radii, images, and spectral indices. Section 5 summarizes our conclusions.

Dust evolution models
We used the dustpy code developed by Stammler & Birnstiel (in prep) based on the dust evolution models from Birnstiel et al. (2010). This code solves the Smoluchowski coagulation equation (Smoluchowski 1916) simultaneously with the transport of the grains. Dust particles collide with each other and are transported in the disk by different mechanisms: thermal Brownian motion, vertical stirring and settling, turbulent mixing, and azimuthal and radial drift. All of these velocities were taken into account to calculate the relative velocities before collision. The output after collision (coagulation, fragmentation, or erosion) is obtained by comparing the relative velocities with the fragmentation speed v f , a parameter that we set to 1 m s −1 in the entire disk based on recent laboratory experiments.
Most of the transport and dynamics of the particles are regulated by the interaction of the dust particles with the gas. The dimensionless Stokes number quantifies the importance of the drag forces on the dynamics of the particles by comparing the stopping time of the particles to the Keplerian frequency at a given disk location. In the Epstein drag regime (Epstein 1924) and near to the disk midplane, the Stokes number is defined as (Brauer et al. 2008;Birnstiel et al. 2010) where a is the grain size, ρ s is the intrinsic volume density of the particles, which we set to 1.2 g cm −3 , and Σ g is the gas surface density. In our simulations we have the following general setup. We follow the dust evolution from 1000 years up to 1 Myr. We assume a disk around a solar-mass star with an initial disk mass of 0.05 M . Our radial grid is from 1 to 300 au with 300 cells that are logarithmically spaced. Initially, the gas-to-dust mass ratio is 100 and all the grains are one micron-sized particles. The particle size grid covers up to 2 meters with 180 cells, also logarithmically spaced. The temperature profile is a power law such that T = 160 K(r/au) −0.5 , but we assume a floor value of 10 K, implying that from ∼250 au the temperature becomes constant. We assume an exponentially tapered power-law function for the initial gas surface density (Lynden-Bell & Pringle 1974) with γ = 1 and a characteristic radius of R c = 80 au. The gas surface density is viscously evolving following the diffusion equation obtained from the continuity equation and angular momentum conservation of the gas in protoplanetary disks and assuming Keplerian frequency everywhere in the disk (Pringle 1981), that is, where ν is the disk kinematic viscosity. We use the parameterization by Shakura & Sunyaev (1973), such that where c s is the isothermal sound speed, h is the gas pressure scale height, and Ω is the Keplerian frequency. The parameter α is assumed to take values that range from 10 −5 to 10 −2 , the former being a value expected in the dead zones of protoplanetary disks (Nelson et al. 2013;Klahr & Hubbard 2014;Lyra 2014;Flock et al. 2015Flock et al. , 2020, and the latter being values expected in the magnetorotational instability active regions (e.g., Fromang & Nelson 2006;Davis et al. 2010;Simon et al. 2011;Flock et al. 2012;Turner et al. 2014).

Distinguishing the effects of gas diffusion versus dust (radial and vertical) diffusion
Current dust evolution models of protoplanetary disks assume a dependence of turbulent velocities, vertical stirring, and radial diffusion of dust particles on the parameter α that regulates the gas viscous evolution. In our simulations, we assume non-isotropic turbulence and thus the parameters that control these physical mechanism are independent of each other. We use the symbols δ t , δ z , and δ r , for the turbulence, settling, and radial diffusion, respectively. We always assume that particles cannot diffuse faster (radially or vertically) than the gas (i.e., δ r,z,t ≤ α).
Vertical stirring and settling. The velocities of very small particles ( 1 µm sized grains) are usually dominated by Brownian motion, but this becomes negligible as soon as particles grow to larger sizes (e.g., Birnstiel et al. 2016, see their Fig. 1). From micron-to sub-millimeter-sized dust grains, settling can be an import source for the velocity of the dust and their collisions. When the stirring by gas turbulence is neglected, the vertical velocity of grains that settle to the midplane (v z,settling ) is given by balancing the gas drag force with the vertical gravitational force, F grav = −mΩ 2 z. This gives the terminal velocity for particles coupled to the gas (St < 1) v z,settling = −zΩSt (Safronov 1969;Dullemond & Dominik 2004). This implies that the settling velocity decreases close to the midplane where the gas density increases, and it also increases for particles with higher Stokes number (for St < 1). Nonetheless, particles are also expected to be stirred up above the midplane by vertical turbulence. Vertical and radial turbulence may differ due to different instabilities, such as the Goldreich-Schubert Fricke instability (Nelson et al. 2013;Stoll & Kley 2014).
The relative velocities between particles due to settling depend on their height above the midplane, and we assume a Gaussian vertical mass distribution for which the scale height depends on the vertical stirring parameter δ z . The particle scale height h d is given by (Youdin & Lithwick 2007;Birnstiel et al. 2010) As demonstrated below, when the relative velocities due to turbulent motion and/or radial drift are low, relative velocities due to settling can dominate, in which case we use the settling limit obtained in Appendix A.
Although we do not include the vertical evolution of the particle density distribution, we calculate the total volume dust density as a function of h d and use these 2D (r,z) particle density distributions in the radiative transfer models. The volume dust density is assumed to follow where z = r cos(θ), with θ being a polar angle; and Σ d is the dust surface density.
Turbulent velocities. For particles larger than submillimeter sizes, turbulence and radial drift can dominate the relative velocities, and both of them depend on the Stokes number. We approximate the turbulent velocities by (Ormel & Cuzzi 2007) When we assume that turbulence dominates the relative velocities of particles and that sticking collisions occur below the fragmentation speed v f , the maximum grain size a f,turb is given by the so-called fragmentation barrier (Birnstiel et al. 2010, We assume values of δ t that are either equal to or lower than α. Physically, this would imply that what sets the gas viscous evolution may be detached from what dominates the turbulent motion of the grains. This is possible for instance in the weakly ionized regions (or dead zones), where the magneto-rotational instability (MRI) is not effective as expected in the dense midplane regions, while what drives the gas accretion and the angular momentum transport in the disk may be driven by other and more effective physical mechanism such as a MHD disk wind (e.g., Bai 2016).
Growth can also be limited by radial drift because the timescale of radial drift is short. When the growth timescales are balanced with the drift timescales (Klahr & Bodenheimer 2006) and when relative particle speeds are assumed to be dominated by turbulent motion, the drift barrier is given by where P is the disk gas pressure, and v K is the Keplerian speed.
The two-population model presented by Birnstiel et al. (2012) takes advantage of top-heavy size distributions and of the increase of the radial drift speed with size (for Stokes numbers below unity). Under these conditions, particle evolution can be simplified by assuming that a particle flux dominates large species and a small species that is radially transported mostly by gas motion and particle diffusion. Birnstiel et al. (2012) assumed in their model that fragmentation occurs as a result of relative turbulent motion or relative radial drift. For the drift limit, it is assumed that turbulent motion dominates. However, for low values of δ t , that is, lower relative turbulent velocities, and/or low particle fragmentation speeds v f , the dominating source of relative velocity can switch to relative vertical settling or relative radial drift in the fragmentation and drift regime. This means that in each regime turbulence, settling, or radial drift can also dominate, and thus each regime can have up to three different barriers, as explained in Appendix A. These limits can be used to further generalize the twopopulation model presented in Birnstiel et al. (2012), which allows us more freedom in choosing the viscosity parameter of the gas α, the particle turbulence parameter δ t , and the vertical (δ z ) and radial (δ r ) particle diffusion parameters separately.
Radial diffusion. For the dust transport, the advectiondiffusion equation of the dust surface density (Σ d ) is solved for each grain size, given by (e.g. Birnstiel et al. 2010) where D d is the dust radial diffusivity, which is assumed to be (Youdin & Lithwick 2007) In the classical models, δ r = α, but here we assumed independent values with the same motivation as for δ t or δ z , that is, the radial diffusion of the grains is independent of the physical mechanism that causes the viscous evolution of the gas.
2.1.2. Combination of the δ t , δ z , δ r , and α parameters Table 1 summarizes the combination of δ t , δ z , δ r , and α that we assumed for the models of this paper. Models 1 and 2 were the two extremes in our set of simulations. Model 1 considered a disk in which the parameter for the viscous evolution of the gas was set to the maximum value (10 −2 ) obtained in the active regions of disks from MRI Table 1. Summary of the models and results without pressure bumps (Fig. 2) Columns 2, 3, 4, and 5 show the values taken for the parameters that regulate the gas viscous evolution (α), radial diffusion (δr), turbulent mixing (δt), and vertical settling/stirring (δz), respectively. Column 6 shows a tick mark when pebbles (a >1 mm) have formed in most of the disk, or a cross otherwise. Column 7 summarizes the relative velocity that sets the maximum grain size, which is always in the fragmentation regime at 1 Myr of evolution. Columns 8 and 9 show the synthetic millimeter fluxes at 1.3 mm and 3.0 mm, respectively. Column 10 is the spectral index. Columns 11 and 12 show the dust radius that encloses 68% or 90% of the 1.3 mm flux. (1) (  calculations (e.g., Turner et al. 2014), and in which the remaining analyzed parameters had the same values. On the other hand, model 2 was the other extreme, in which the parameter for the viscous evolution of the gas and all the others were set to a very low value that is expected in a disk in wich MRI is inefficient (10 −5 ). The values for the parameters to explore in models 3 to 9 were such that they covered a transition from the two extremes of model 1 and 2. More specifically, models 7 and 9 also assumed that all the four parameters are equal, and the values were 10 −3 and 10 −4 , respectively. The remaining models (models 3, 4, 5, 6, and 8) mixed different values to resemble different radial/vertical diffusion and/or turbulent velocities of the particles.

Inclusion of a pressure bump
We followed the prescription introduced by Dullemond et al. (2018), which was also used by Stammler et al. (2019), to include a pressure bump in the disk. This assumes a bump that originates from variations of the viscosity parameter (α), modeled as where F (r) is given by From previous inspections of the results of models without bumps, we performed six simulations with this pressure bump, as summarized in Table 2. In these models, only α varied with radius, and the remaining parameters were set as before for comparison (see Table 2 for more details). We assumed the position of the bump at r p = 40 au, with an amplitude of f = 2 (as in Dullemond et al. 2018;Stammler et al. 2019), and a width of w gap = 5 au, which is approximately twice the local pressure scale height at r p = 40 au. However, the pressure bump itself formed farther out, at ∼ 52 au, where the pressure scale height is ∼3.5 au.

Radiative transfer calculations
To compute synthetic images and spectral energy distributions (SEDs) for each model, we considered the resulting vertical particle distribution (Eq. (6)) after 1 Myr of evolution as input to the radiative transfer calculations, for which we used RADMC-3D (Dullemond et al. 2012). We calculated the opacity of each grain size considering porous spheres with a dust mixture composed of astronomical silicates (Draine 2003), carbonaceous material (Zubko et al. 1996), and water ice Warren & Brandt (2008), as in Ricci et al. (2010). We assumed a blackbody radiation field from the central star as the radiation source and used 1 × 10 7 photons for our calculations.
We obtained the temperature profile for each grain size and calculated scattered-light images at 1.25 µm and 1.3 mm. We assumed a disk inclination of 20 • , a position angle of zero, and a distance to the source of 140 pc as the typical distances to nearby star-forming regions, such as Lupus, Taurus, and Ophiuchus. The synthetic scatteredlight images were computed including the full treatment of anisotropic scattering with polarization, as in Kataoka et al. (2015) and Pohl et al. (2016). Afterward, the images were convolved by a Gaussian point-spread function (with a size of 0.04" × 0.04") at the two wavelengths, which is similar to the angular resolution that can currently be obtained with Spectro-Polarimetric High-contrast Exoplanet REsearch at the Very Large Telescope (SPHERE/VLT) and ALMA. The intensity profiles created at 1.25 µm were then multiplied with the square of their distance to the star to compensate for the stellar illumination. From the intensity profiles at 1.3 mm we calculated the radius that enclosed either 68% or 90% of the emission by assuming that the observation sensitivity provides a signal-to-noise ratio with respect to the maximum of the emission (peak) of 1000. . We do not include the turbulence I regime or the azimuthal drift because they do not dominate the relative velocities in any case for the shown snapshot. The horizontal gray line shows a = 1 mm, and we call pebbles all particles that are larger than this.

Models without pressure bumps
The particle density distribution as a function of distance from the star and grain size of the models without any pressure bump (Table 1) is shown in Fig. 2. All panels show the results at 1 Myr of evolution (as the typical age of class II protoplanetary disks in nearby star forming regions like Taurus and Lupus) and in all cases the fragmentation speed is v f = 1m s −1 . Models 1 and 2 are the two extremes of our assumptions. For all the parameters that set the gas and dust diffusion as well as relative turbulent speeds, model 1 adopted the same value of 10 −2 , and model 2 corresponds to the case of 10 −5 . As expected in model 1, pebble (a > 1 mm) formation did not occur. The maximum grain size in the case of model 1 was set by high velocities because turbulence induced the fragmentation barrier (a f,turb , Eq. 8), and the maximum grain size was approximately 0.01 cm at 1 au and about 10 µm in the outer part of the disk (at about 100 au).
On the other hand, model 2 allowed the formation of pebbles in most of the disk, where the maximum grain size was set by fragmentation, for which the radial drift velocities dominate.
The results of model 1 and 2 are expected because a f,turb depends quadratically on the fragmentation velocity. When we assume a fragmentation velocity that is one order of magnitude lower than previously used, theory requires us to decrease the dust turbulence by at least two orders of magnitude if turbulence is the main source of relative velocities. Model 9 was the same as model 2, but with 10 −4 . It shows that the maximum grain size for r < 50 au is set by fragmentation by turbulence and by drift in the outer regions. In this case there is still some pebble formation. Model 7 assumed all gas and dust diffusion parameters to be equal to 10 −3 , and pebble formation did not occur in most of the disk, as in model 1. The first conclusion of these numerical experiments therefore is that with isotropic turbulence and v f = 1m s −1 , pebble formation can occur in However, the main motivation for the other experiments of this work was first that when α has values of 10 −5 −10 −4 , the disk accretion rate at million-year timescales are expected to be 10 −12 − 10 −11 M yr −1 , which is much lower than the typical accretion rates in different star-forming regions of different age (1-10 Myr, see, e.g., Manara et al. 2020). For values of α of 10 −3 − 10 −2 , the accretion rates values agree with the mean values from observations (i.e., 10 −9 − 10 −10 M yr −1 ). In the remaining numerical experiments we therefore kept α = 10 −2 or α = 10 −3 and varied the parameters responsible for the radial and vertical particle diffusion as well as the parameter that controls the collision speed due to turbulent motion.
The second motivation for these numerical experiments was that several protoplanetary disks appear to be flared when observed in scattered light (Avenhaus et al. 2018). Assuming that we trace micron-sized particles with these observations, and taking values of δ z = 10 −5 −10 −4 , micronsized particles would be mainly settled in the midplane, contradicting observations.
The third motivation is that the width of some rings in protoplanetary disks observed at very high resolution with ALMA have been resolved (e.g., Dullemond et al. 2018, with values of 4-20 au). This means that there must be radial diffusion inside pressure bumps to widen the distribution of millimeter-sized particles. All the observed rings would otherwise be "infinitely" narrow (dust grains would move to the location of the pressure maximum).
Our models 3, 4, 5, and 6 maintained high values of α, δ r , and δ z to guarantee accretion-rate values, width of rings, and dust aspect ratios that are in accordance with observations, while the parameter that controls the turbulent velocities (δ t ) of the particles was low. Model 8 is an experiment in which we also decreased δ z to 10 −4 because by previous inspections of the models, we realized that settling can set the maximum grain size and limit pebble formation. Table 1 shows when pebble formation occurred in our models and what set the maximum grain size in each case. At 1 Myr of evolution, all the cases are in the fragmentation regime, but the relative velocity that dominates and sets the maximum grain size is different.
From our numerical experiments in which we used different values for δ t , δ z , δ r , and α (i.e., models 3, 4, 5, 6, and 8), the maximum grain size was set either by settling or radial drift. In models 4 and 5, pebble formation is not possible in almost the entire disk (r > 2 au) when δ z = 10 −2 . The reason is the average settling velocities dominate the total relative velocities and set the maximum grain size. Models distribution, but the radial extension of the large grains is slightly larger for model 4, as we discuss in the Sect. 4 when we compare the results with observations. For models 3 and 6 with δ z = 10 −3 , the maximum grain size was mostly set by settling and was set by drift in only the farthest regions. In these two cases, pebble formation occurred for r < 10 au in model 3 and for r < 50 au for model 6. Model 6 produced pebbles in a wide range of the disk. This model is a good compromise for a comparison with observations, that is, an α = 10 −3 that gives accretion rates of 10 −9 M yr −1 , a δ z = 10 −3 that can lead to enough vertical stirring to produce a high dust-aspect ratio similar to those observed with SPHERE for T-Tauri and Herbig disks, and a δ r that can lead to radial diffusion and a wider distribution of pebbles (also for the potential presence of a pressure bump).
No direct observational constraint is available that could help us to fix δ t . Indirect measurements of the disk turbulence have been made using ALMA observations by searching for variations from Keplerian motion of the gas emission that may come from turbulence. Flaherty et al. (2015Flaherty et al. ( , 2017Flaherty et al. ( , 2020; Teague et al. (2016Teague et al. ( , 2018 found low turbulence values for four disks ( 10 −3 ). Flaherty et al. (2020) showed that DM Tau is the only system observed so far with high turbulence. Figure 3 shows the 2D (radial and vertical) particle density distribution assumed for the radiative transfer calculations and obtained using Eq. (6). We show the vertical and radial particle density distribution considering all the grain sizes. The greatest difference is seen between our two extremes model 1 (α, δ r,z,t = 10 −2 ) and model 2 (α, δ r,z,t = 10 −5 ), where model 2 clearly presents much less stirring of dust particles from the midplane, making the dusty disk thinner. Another clear result from Fig. 3 is that when the maximum grain size is set by the fragmentation limit due to radial drift, the disk becomes thinner in all grains. This occurs in all the disk for models 2 and 8 in which radial drift dominates relative velocities everywhere, and from the location where the barrier is set by radial drift and outward for models 3, 6, and 9. In addition, in models where δ z ≤ 10 −4 (models 2, 3, 8, and 9), the particle density distribution is much more dominated by the midplane than in the other models.

Models with a pressure bump in the outer disk
We chose our two extreme cases (model 1 and 2) for the inclusion of a pressure bump, together with the other two cases in which all the gas and dust diffusion parameters are constant (models 7 and 9). In addition to these cases, we included two other cases, one in which the maximum grain size was set by settling in the whole disk (model 4). The other case was model 6, in which pebble formation occurred, and the maximum grain size was set by a mix of settling and drift. The summary of the models that include a pressure bump as described in Eqns. (12) and (13) is listed in Table 2.
The particle density distribution as a function of distance from the star and grain size for the case of a gap at 40 au and a pressure bump at the edge of that gap is shown in Fig 4. Because this gap is formed in the α viscosity, any variation in the value of α changes the timescale within which this gap and bump are formed. Because all the plots are shown at 1Myr, in model 11 when α is the lowest, the gap is still opening and the amplitude of the bump is lowest. The gap is also shallower for model 15 (α = 10 −4 ) than in the rest of the models.
There is effective trapping, reflected by the significant particle accumulation at the pressure maximum, when pebble formation occurs (models 11, 13, and 15). In the other  Table 2. Summary of the models and results with a pressure bump (Fig. 4). In Column 3 a tick mark is given when efficient trapping occurs in the pressure bump. A cross is shown otherwise. (1) (2) cases, the strong radial mixing (both through the radial particle diffusion parameter and the radial gas velocity by which small particles are dragged along) does not lead to significant trapping and hence to no significant accumulation of particles at the pressure maximum (models 10, 12, and 14). When there is no efficient trapping (models 10, 12 and 14), there is a small decrease in the maximum particle size at the location of the gap because the fragmentation barrier is proportional to the gas surface density (regardless of whether it is set by turbulence or settling). When there is no trapping and the maximum grain size is set by settling (model 12), the decrease in particle size at the gap location is smaller.
For effective particle trapping, model 11 (α, δ r,z,t = 10 −5 ) shows a situation where grain growth is very efficient in the pressure maximum and particles grow to the maximum size that is allowed in our grid (2 m) at that location. The reason is that the relative velocity that dominates fragmentation is radial drift, which is zero at pressure maximum. This has been seen before in simulations with particle traps and low turbulence and diffusion of the grains (e.g., Pinilla et al. 2015). Models 13 (α, δ r,z = 10 −3 , and δ t = 10 −5 ) and 15 (α, δ r,z,t = 10 −4 ) where effective trap-ping occurs, settling or turbulence set the maximum grain size inside the pressure bump, respectively. The width of the dust accumulation inside the gap decreases for lower δ r , that is, model 11 (δ r = 10 −5 ) shows the narrowest accumulation, followed by model 15 (δ r = 10 −4 ), and then by model 13 (δ r = 10 −3 ). Figure 5 shows the vertical and radial distribution integrated over all the grain sizes for each model that includes a pressure bump. When there is no efficient dust trapping (models 10, 12, and 14), the decrease in particle density distribution is small, which is more evident in the upper layers of the disk, but it is still noticeable at the midplane.
For effective trapping (models 11, 13, and 15), the depletion of grains extends from the surface layers to the midplane, especially in cases 11 and 15, where the radial diffusion and vertical stirring are low. At the location of the trap, small grains clearly increase as a result of continuous fragmentation of dust particles, as seen before in Pinilla et al. (2012) and Pinilla et al. (2015), for example, but also in 2D dust evolution simulations by Drążkowska et al. (2019). Because of this local increase in small particles, their vertical distribution increases through vertical stirring. This is more clear in models 13 and 15, where the vertical stirring parameter is higher than in model 11. Model 13 has higher Article number, page 9 of 17 A&A proofs: manuscript no. MAIN radial diffusion, which causes a significant number of grains throughout the gap, more than in models 11 or 15.

Millimeter fluxes and spectral indices
Tables 1 and 2 provide the expected millimeter fluxes at 1.3 mm (ALMA Band 6) and 3.0 mm (ALMA Band 3) after taking the particle density distribution after 1 Myr of evolution. They then use radiative transfer simulations to create the SEDs. Model 1, in which the fragmentation barrier is the lowest and where pebble formation does not occur in the disk, has the lowest millimeter fluxes. In general, the cases where pebbles are not formed or they are formed in just a small region of the disk have the lower millimeter fluxes, except for model 7, in which all the gas and dust diffusion parameters were set to 10 −3 . None of the models show a correlation between the velocities that set the maximum grain size and the obtained fluxes. All of the fluxes are within the observed values of the disks in different star-forming regions around T-Tauri and Herbig stars (e.g., Ricci et al. 2012;Tripathi et al. 2017;Andrews et al. 2018b).
When we compare the models with and without a bump in the disk, in all of the cases (except for model 6 vs. model 13), the millimeter fluxes decrease when the bump is included. When trapping occurs, this originates from (i) grains growing to decimeter objects, which have very low opacities and contribute very little to the total fluxes (model 2 vs. model 11), or/and (ii) the emission inside the trap becoming optically thick (as in the case of model 15). In model 13 (model 6, but with a bump), grains in the trap have a maximum size of ∼1 cm and the emission is optically thin such that the millimeter fluxes slightly increase with the trapping, as expected because fewer grains are drifting toward the star. In the cases without efficient trapping, the fluxes decrease compared to the same cases without a bump. In all these cases, we have a local decrease of the maximum grain size at the gap location, which helps to decrease the total millimeter flux.
For all the models, we obtain the millimeter spectral index (α mm ) taking the fluxes at 1.3 mm and 3.0 mm, such that α mm = log 10 (F 1.3mm /F 3.0mm )/ log 10 (3.0/1.3). All the spectral indices in our work are much higher than observed (e.g., Ricci et al. 2012;Pinilla et al. 2017b). As shown in Pinilla et al. (2012), high spectral indices from dust evolution models are expected when radial drift is included. Even when pebble formation occurs and the millimeter fluxes agree with observations, the number of pebbles is not enough to have low spectral indices. Even for models with effective trapping, the values of the spectral index remain similar as in models without a bump, implying that this bump does not maintain enough pebbles to decrease the values of the spectral index (only 1-2% of the total dust mass is in pebbles). To solve the discrepancy, strong multiple bumps can assist in decreasing α mm (Pinilla et al. 2012). Alternatively, when only one single bump is included, it has to be strong enough, as the one expected from a giant planet, to decrease the spectral index (Pinilla et al. 2014). In this case, the spectral index depends on where the bump is located, becoming larger when the bump is located farther out. Another factor that can affect the final spectral indices is the lifetime of the pressure bumps because in previous works (e.g., Pinilla et al. 2012Pinilla et al. , 2013Pinilla et al. , 2014, pressure bumps were included from the beginning of the simulation and they are static, while in the present work the gap and bump grow on viscous timescales. Both completely static bumps and bumps that slowly grow with time have their limitations. A prescription of a bump in α allows us to have mass accretion onto the star as well as viscous spreading, but for very low values of α it can take more than 1 Myr to reach the final amplitude of the bump, which may not be the case of planets that cause the bump formation, for instance. A more realistic approach would include the viscous evolution and the grow, appearance, and disappearance of a number of bumps with a given frequency during the disk lifetime. To test whether multiple gaps as assumed in this paper versus static multiple bumps affect the spectral index, we performed four more models considering the α and δ r,z,t from model 9 (i.e., α, δ r,z,t = 10 −4 ), and we took five gaps or bumps with the same width of 2 local pressure scale heights. They were located at 20, 40, 60, 80, and 100 au. We list these models below.
-A model with five gaps as prescribed in this paper (Eq. 12 and Eq. 13) with f = 2 (name: smooth gaps in α). -As the previous model, but with f = 4 (name: strong gaps in α). -A model with five bumps as described in Pinilla et al. (2020), which were introduced as perturbations in the initial gas surface density profile. They remained static. Specifically, with The locations (r p,i ) of these perturbations were at 20, 40, 60, 80, and 100 au, and the width (w i ) of each of them was two local pressure scale heights. The amplitude (A i ) for all was 2 (name: smooth static bumps). -As the previous mdoel, but with an amplitude of 4 (name: strong static bumps).
The results of these tests are summarized in Fig 6, which shows the particle density distribution at 1 Myr of evolution. When the perturbations were introduced in the gas surface density (Eq. 14), the bumps were located at the position of the gaps in the other two cases. In addition, for static bumps, the perturbations overlapped in the outer part of the disk. This causes fewer bumps.
After we took these dust density distributions and performed the radiative transfer calculations, the spectral index for the models in which the gaps were introduced in α was 2.8, while the models in which the bumps were introduced in the gas surface density and remained static during the dust evolution was 2.5. The last values are similar to the observed values in different disks in several star-forming regions.
These results show that including more bumps help to decrease the spectral index. The fact that they live from the starting point of the simulations with a given amplitude decreases the values of the spectral index to the averaged values that have been observed. Because the bumps formed by variations in α evolve on viscous timescales, less trapping is expected than static pressure bumps if these times are longer than the growth timescale. The amplitude here does not seem to play an important role, probably because with this prescription of the bumps and gaps, there is not full filtration, or in other words, radial mixing still allows the grains to move out of the pressure bumps. A deeper investigation is required to confirm this.

Dust-disk radii
Using the intensity profile obtained from the synthetic maps at 1.3 mm, we calculated the radius that enclosed either 68% (R dust,68% ) or 90% (R dust,90% ) of the emission by assuming a signal-to-noise ratio (S/N) with respect to the maximum emission of 1000 in all cases. We decided to have an S/N compared to the peak instead of a constant S/N for all of the cases because most of the observations available with ALMA have a different sensitivity, all depending on the total flux of each disk and the observation goals. Currently, this S/N can be obtained for the bright disks with ALMA with a short (some minutes) integration time. This S/N is very challenging to achieve for faint disks, however, and therefore the following discussion is valid when observations with similar S/N with respect to the peak are compared. Tables 1 and 2 include these two dust radii for each simulation.
Our models without pressure bumps (except for model 1) span a wide range of values that agree with the values in different star-forming regions, such as Ophiuchus, Taurus, and Lupus (Hendler et al. 2020). The mean value of R dust,68% in these models (excluding model 1) is ∼ 53 au. In these models, the larger disks are also brighter. Our results demonstrate that a combination of gas and dust diffusion parameters can lead to very different particle density distributions, which directly affects the outer dust radius as obtained from millimeter observations. Model 1 is an outlier where the flux is the lowest we found, but this is the largest disk according to our models. This is a consequence of the high turbulence and radial diffusion because the high turbulence produces small grains that are easier to diffuse radially. However, this model has the lowest millimeter flux as well, meaning that we would need deeper observations to detect this outer radius. This result would imply that disks can appear large when observed at millimeter emission without the need of any pressure bump. Recent ALMA observations suggest that large disks remain bright because of substructures, and most of the small disks seem to be faint so far (Long et al. 2019). Current observations pose the question of whether the absence of structures is the reason for the small and faint disks. The current data have some observational bias because deeper observations are needed to determine the outer radius of fainter disks and higher resolution is required to detect structures in the very small disks.
The effect of radial diffusion can be more clearly seen when we compare model 4 and model 5, where the only parameter that is different between these two models is δ r . Model 5 has a δ r that is one order of magnitude lower and the disk appears to be smaller, both in the R dust,68% and R dust,90% , with a higher difference when R dust,90% is compared, which is twice as high for the disk with the higher radial diffusion (198 au for model 4 vs. 95 au for model 5). In these two models, the flux is very similar, and therefore this result is independent of the S/N of the real observations.
In our simulations without pressure bumps, R dust,90% can be 1.5 to 3.6 larger than R dust,68% , but this strongly depends on the assumption of the S/N. For instance, when we assume half of the S/N of our calculations, R dust,68% remains similar (variations within 5 au), while R dust,90% decreases such that R dust,90% can be only up twice R dust,68% .
When models with and without a pressure bump are compared, the dust radii increase for models with a bump, and this increase becomes higher when trapping is efficient (models 11, 13, and 15). The pressure bump is located at ∼52 au in all the models, and the outer radius is larger than this position, either R dust,68% or R dust,90% . Pinilla et al. (2020) showed that for strong and static pressure bumps, the outer radius traces the location of the farthest pressure bump in the disk. In this case, we do not see this trend, which is probably due to the nature of the pressure bump, which in this work develops with the gas viscous evolution and has a lower amplitude than in the models by Pinilla et al. (2020). The only model in which the outer dust radius does not increase is model 6, which might be the result of the high radial diffusion when compared to the other two cases with trapping. Figure 7 shows the radial intensity profile of the synthetic images (previously convolved with a Gaussian beam with a full width at half maximum, FWHM, of 0".04 in each case) after radiative transfer calculations assuming the particle density distributions without any pressure bumps from Fig. 3. The intensity profiles were obtained at 1.25 µm and 1.3 mm. The intensity profile at 1.25 µm was multiplied by r 2 to compensate for the stellar irradiation. The gray region shows the typical size of a coronograph or/and saturated inner pixels from typical observations at scattered light. after radiative transfer calculations assuming the particle density distributions from Fig. 3. The intensity profiles are obtained at 1.25 µm and 1.3 mm. The intensity profile at 1.25 µm is multiplied by r 2 to compensate for the stellar irradiation. The gray regions show the typical size of a coronograph or/and saturated inner pixels from typical observations at scattered light. Case of models without bumps.

Substructures
All the radial profiles of the intensity at 1.25 µm show a large gap, as reported by Birnstiel et al. (2015), even in the absence of any perturbation in the disk or pressure bump. This gap in scattered light originates from the removal of small grains by sweep-up collisions with the large particles. For this gap to form, there must be an equilibrium between vertical settling/stirring and drift, such that the large grains that are drifting inward are able to sweep-up the small grains in the disk surface. This gap forms in all our models, and the minimum of this gap is located at ∼ 25 au for all the cases. For most of the simulations, the width of the gap (extend up to ∼60-80 au) and its contrast (a factor of 10) are similar. Model 1 has the narrower and lowest contrast gap, while model 6 has the widest gap in our simulations (up to ∼100 au). Birnstiel et al. (2015) described the location of this gap based on the model assumptions.
The intensity at 1.3 mm of the synthetic observations of the models without any pressure bumps decreases smoothly with radius for all the models. In some cases, a local increase in intensity is observed close to the outer radial edge of our simulations. This structure should be taken with caution because it is an effect of assuming a floor value of 10 K for the temperature profile from ∼250 au. Because of this assumption, the radial profile of the pressure gradient becomes less steep and the large dust particles are locally enhanced. The radial profiles resemble a smooth decrease in the inner parts together with a sharp decrease in the outer disk. Similar to the dust-disk radii discussed in Sect 4.2, the location at which the transition from a smooth and sharp profile occurs depends on the gas surface density profile and its vacuous spread outward, the particle size, and on all the dust diffusion parameters that are involved, as discussed in Birnstiel & Andrews (2014).
For the models with a bump and in the cases without effective trapping (models 10, 12, and 14), there is a gap in the intensity profile at both wavelengths with a contrast equal to or higher than in scattered light at 1.25 µm than in the continuum emission at 1.3 mm. This is the opposite result: a long-lived and strong trap, from giant planet-disk interaction, where the millimeter gaps have higher contrast (e.g., de Juan Ovelar et al. 2013;Dong et al. 2015). This characteristic, that the gaps have higher contrast in scattered light than at the millimeter, resembles the results from dust-evolution models that include different disk icelines (Pinilla et al. 2017a), where trapping does not occur. Instead, there is a traffic-jam effect that is caused by the changing size of the dust particles at the location of the ice lines, which decreases or increases their radial drift. The gap in the absence of effective trapping appears to be due to the local decrease in maximum grain size at the gap location. When turbulence sets the maximum grain size (models 10 and 14), the contrast of the scattered light gap (measured as the ratio of the emission of the maximum at the outer edge of the gap and the minimum of the intensity at the gap location) is higher than at millimeter sizes, with a contrast of a factor of 10 in scattered light versus a contrast of a factor of 4-5 at millimeter emission. Without trapping and when the maximum grain size is set by settling (model 12), the contrasts of the gap (a factor of 2) and the width are the same at the two wavelengths.
For effective trapping (models 11, 13, and 15), all the gaps have higher contrast at 1.3 mm than at 1.25 µm, as expected from previous models (e.g., de Juan Ovelar et al. 2013). In model 11, the difference in the contrast between the two wavelengths is lowest, which may have two reasons. First, the large grains formed at pressure maximum, whose opacities are low and that do not significantly contribute to the 1.3 mm emission. Second, because the gap and bump are introduced as variations in α-viscosity and they evolve on a viscous timescale, the gap in model 11 (with α = 10 −5 ) takes longer to open and reach the final stage. At 1 Myr of evolution, the bump therefore does not have the same amplitude as in the other cases.
The contrast of the gap is higher in model 15 (α, δ r,z,t = 10 −4 ) than model 13 (α, δ r,z = 10 −3 , and δ t = 10 −5 ) because the radial diffusion is lower. Lower radial diffusion of the dust helps to keep more grains of different size inside the pressure bump, which contributes to increase the contrast in both wavelengths.
In all cases, the gap in scattered light seems to be as narrow or wider than at millimeter emission, which is the opposite from models of trapping with a strong and longlived pressure bump such as those expected from giant planets. The nature of this pressure bump forming with viscous timescales and having lower amplitudes leads to this difference. Therefore the final gap seen in the scattered light in the simulations with a pressure bump may result from a combination of a gap formed purely by dust evolution processes together with the effect of the bump in the distribution of the small dust particles. This characteristic of a wider gap in scattered light versus millimeter emission has been observed in several disks such as TW Hya and HD 163296 (see Fig. 6 in Pinilla et al. 2017a). However, a deeper investigation that studies how the lifetime of the pressure bumps affects different disk properties is required to test this conclusion.

Summary and conclusions
We investigate the conditions that are required in dust evolution models for growing and trapping pebbles in protoplanetary disks when the fragmentation speed is 1 m s −1 in the entire disk. In particular, we separated the parameters controlling the effects of turbulent velocities, vertical stirring, radial diffusion, and viscous evolution of the gas, always assuming that particles cannot diffuse faster (radially or vertically) than the gas. For some of our simulations, we included a smooth gap in the disk by variations of the α viscosity parameter, which led to a pressure bump in the outer disk (∼52 au). The simplified model by Birnstiel et al. (2012) does not include the effects of settling-dominated relative particle speeds. Furthermore, for the drift limit they did not consider radial drift-dominated growth. We extended the potential growth limits, including the effects of relative speeds dominated by settling, radial drift, and two different turbulence regimes. The calculations are shown in Appendix A. We compared our models with observations of protoplanetary disks in near-infrared and millimeter regimes. Our main conclusions are listed below.
-With isotropic turbulence, pebble formation can occur in disks when the gas and dust diffusion is low, with values of 10 −5 −10 −4 . Otherwise, pebble formation does not occur, and efficient accumulation of particles in a pressure bump is impossible. -When relative velocities due to turbulence or radial drift are low, settling can dominate the total relative velocities and can set the maximum grain size. In these cases, pebble formation is not possible when δ z > 10 −3 , which also limits the efficiency of trapping when a pressure bump is present. -In general, the millimeter fluxes are lower when pebbles are not formed. In none of the cases did we see a correlation between the processes that set the maximum grain size and the obtained fluxes. All of the fluxes are within the observed values of the disks in different star-forming regions around T-Tauri stars. The fluxes do not change significantly when a trap is included. -Most of the spectral indices in our work are much higher than observed because the number of pebbles is too low as a result of the effective drift. Even for models with one pressure bump and effective trapping, the spectral indices are higher than observed. To solve the discrepancy, multiple bumps assist to decrease α mm . The formation of the bumps and the formation time also affect the spectral index. Static multiple pressure bumps produce lower spectral indices than models in which the gaps and bumps are introduced as variations in α-viscosity. -Models without pressure bumps span a wide range of synthetic values of R dust,68% that agree with the values in different star-forming regions. When models with and without a pressure bump were compared, the dust radii increased for models with a bump, and this increased when trapping was efficient. -One of our simulations without a pressure bump (model 1, see Table 1) produced a large disk as a consequence of the high turbulence, high radial diffusion, and high α that causes the disk to spread and drag particles along. High turbulence limits the formation of pebbles, and small grains are continuously produced by fragmentation; they are easier to diffuse radially. This result would imply that disks can appear large when observed at millimeter emission without the need of any pressure bump. Deeper observations at high angular resolution are needed to test this idea and exclude that substructures are the reason for a potentially large size of a disk.
In the case of a bump and effective trapping, high radial diffusion also helps to widen the concentration of large particles in pressure maxima.
-As found in Birnstiel et al. (2015), even in the absence of any perturbation in the disk or pressure bump, a large gap is formed in the synthetic scattered-light images at 1.25 µm. For all our models, the minimum of this gap is located around ∼ 25 au with a similar width (extend up to ∼60-80 au) and contrast (a factor of 10). -When a pressure bump is included but there is no effective trapping, a gap opens in the intensity profile at short and long wavelengths, whose contrast is equal to or higher than in scattered light than at millimeter emissions. These results are similar to traffic-jam models, for example, to the dust evolution simulations that include the effect of different disk ice-lines. This is the opposite result of effective trapping by a long-lived and strong trap, for example, from giant planet-disk interaction, where the millimeter gaps have higher contrast.