© P. Pinilla et al. 2021

1 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. 2012, 2020; Dullemond 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, fragmentation, 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 2000, 2008; Kimura et al. 2015, 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. 2009, 2011; Gundlach & Blum 2015; Musiolik et al. 2016), and the comparison has been made that this would be like snowballsversus 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 ofpebbles 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. 2010a). 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 temperatureand 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⁻¹ 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⁻¹ 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⁻¹. 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.

Fig. 1 Schematic view of turbulent velocities, vertical stirring and settling, and radial diffusion of dust. 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 (see Sect. 2.1.1).

2 Models

2.1 Dust evolution models

We used the dustpy code developed by Stammler & Birnstiel (in prep.) based on the dust evolution models from Birnstiel et al. (2010a). 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 Brownianmotion, 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⁻¹ in the entiredisk 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. 2010a) $$\text{St}=\frac{a{\rho }_{\text{s}}}{{\Sigma }_{\text{g}}}\frac{\pi }{2},$$(1)

where a is the grain size, ρ_s is the intrinsic volume density of the particles, which we set to 1.2 g cm⁻³, and Σ_g is the gas surface density.

In our simulations we have the following general setup. We follow the dust evolution from 1000 yr 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\text{K}{(r/\text{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) $${\Sigma }_{\text{g}}(r,t_0)={\Sigma }_0{\left(\frac{r}{R_{\text{c}}}\right)}^{-\gamma }\exp \left[-{\left(\frac{r}{R_{\text{c}}}\right)}^{2-\gamma }\right]$$(2)

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, $$\frac{\partial {\Sigma }_{\text{g}}(r,t)}{\partial t}=\frac{3}{r}\frac{\partial }{\partial r}\left[r^{1/2}\frac{\partial }{\partial r}(\nu {\Sigma }_{\text{g}}{r}^{1/2})\right],$$(3)

where ν is the disk kinematic viscosity. We use the parameterization by Shakura & Sunyaev (1973), such that $$\nu =\alpha c_{\text{s}}h\text{with}h=\frac{c_{\text{s}}}{\Omega },$$(4)

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⁻⁵ to 10⁻², 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. 2015, 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).

2.1.1 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Ω²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. 2010a) $$h_{\text{d}}(\text{St})=h\times \text{min}\left(1,\sqrt{\frac{{\delta }_z}{\text{min}(\text{St},1/2)(1+{\text{St}}^2)}}\right).$$(5)

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 $${\rho }_{\text{d}}(r,\phi ,z,\text{St})=\frac{{\Sigma }_{\text{d}}(r,\text{St})}{\sqrt{2\pi }{h}_{\text{d}}(r,\text{St})}\exp \left(-\frac{z^2}{2h_{\text{d}}^2(r,\text{St})}\right),$$(6)

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) $$v_{\text{f}}{}^2~3{\delta }_{\text{t}}\text{St} c_s^2.$$(7)

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. 2010a, 2012), $$a_{\text{f},\text{turb}}=\frac{2}{3\pi }\frac{{\Sigma }_{\text{g}}}{{\rho }_s{\delta }_{\text{t}}}\frac{v_{\text{f}}{}^2}{c_{\text{s}}^2}.$$(8)

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 (Birnstiel et al. 2012) $$a_{\text{drift}}=\frac{2{\Sigma }_{\text{d}}}{\pi {\rho }_{\text{s}}}\frac{v_{\text{K}}^2}{c_{\text{s}}^2}{\left|\frac{\text{d}\ln P}{\text{d}\ln r}\right|}^{-1},$$(9)

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 radialdrift speed with size (for Stokes numbers below unity). Under these conditions, particle evolution can be simplified by assumingthat 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 two-population 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 advection-diffusion equation of the dust surface density (Σ_d) is solved for each grain size, given by (e.g. Birnstiel et al. 2010a) $$\frac{\partial {\Sigma }_{\text{d}}}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left(r{\Sigma }_{\text{d}}{v}_{\text{r},\text{d}}\right)-\frac{1}{r}\frac{\partial }{\partial r}\left[r{\Sigma }_{\text{g}}{D}_{\text{d}}\frac{\partial }{\partial r}\left(\frac{{\Sigma }_{\text{d}}}{{\Sigma }_{\text{g}}}\right)\right]=0,$$(10)

where D_d is the dust radial diffusivity, which is assumed to be (Youdin & Lithwick 2007) $$D_{\text{d}}=\frac{{\delta }_{\text{r}}{c}_{\text{s}}h}{1+{\text{St}}^2}.$$(11)

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⁻²) obtained in the active regions of disks from MRI 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⁻⁵). The values for the parameters to explore in models 3–9 were such that they covered a transition from the two extremes of model 1 and 2. Morespecifically, models 7 and 9 also assumed that all the four parameters are equal, and the values were 10⁻³ and 10⁻⁴, 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.

2.1.3 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 $$\alpha \prime \left(r\right)=\frac{\alpha }{F\left(r\right)},$$(12)

where $F\left(r\right)$ is given by $$F\left(r\right)=\exp \left[-f\exp \left(-\frac{{\left(r-r_{\text{p}}\right)}^2}{2w_{\text{gap}}^2}\right)\right].$$(13)

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.

2.2 Radiative transfer calculations

To compute synthetic images and spectral energy distributions (SEDs) for each model, we considered the resulting vertical particledistribution (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⁷ 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 scattered-light 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 (S/N) with respect to the maximum of the emission (peak) of 1000.

3 Results

3.1 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 = 1 m s⁻¹.

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⁻², and model 2 corresponds to the case of 10⁻⁵. 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⁻⁴. 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⁻³, 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 = 1 m s⁻¹, pebble formation can occur in disks when the gas and dust diffusion are low, with values of 10⁻⁵−10⁻⁴.

However, the main motivation for the other experiments of this work was first that when α has values of 10⁻⁵−10⁻⁴, the disk accretion rate at million-year timescales are expected to be 10⁻¹²−10⁻¹¹ M_⊙ yr⁻¹, 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⁻³−10⁻², the accretion rates values agree with the mean values from observations (i.e., 10⁻⁹−10⁻¹⁰ M_⊙ yr⁻¹). In the remaining numerical experiments we therefore kept α = 10⁻² or α = 10⁻³ 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⁻⁵−10⁻⁴, micron-sized 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 ALMAhave 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⁻⁴ 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⁻². The reason is the average settling velocities dominate the total relative velocities and set the maximum grain size. Models 4 and 5 are identical except for the radial diffusion, which is lower in model 5. This creates a small effect in the particle 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⁻³, 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⁻³ that gives accretion rates of  10⁻⁹ M_⊙ yr⁻¹, a δ_z = 10⁻³ 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. (2015, 2017, 2020); Teague et al. (2016, 2018) found low turbulence values for four disks (≲10⁻³). 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⁻²) and model 2 (α, δ_r,z,t = 10⁻⁵), 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⁻⁴ (models 2, 3, 8, and 9), the particle density distribution is much more dominated by the midplane than in the other models.

Fig. 2 Particle density distribution after 1 Myr of evolution as a function of distance from the star and grain size of the models without any pressure bump (Table 1). The solid white line represents a Stokes number of unity. Cyan and red show the fragmentation (Eq. (A.17)) and drift (Eq. (A.32)) limit. The style of the line shows what dominates the relative velocity: turbulence (Eqs. (A.16) and (A.27)), radial drift (Eqs. (A.10), and (A.24)), and settling (Eqs. (A.15) and (A.25)). 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.

3.2 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 Eqs. (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 1 Myr, 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⁻⁴) 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 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⁻⁵) 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⁻³, and δ_t = 10⁻⁵) and 15 (α, δ_r,z,t = 10⁻⁴) where effective trapping 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⁻⁵) shows the narrowest accumulation, followed by model 15 (δ_r = 10⁻⁴), and then by model 13 (δ_r = 10⁻³).

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 incases 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 radial diffusion, which causes a significant number of grains throughout the gap, more than in models 11 or 15.

Fig. 3 Vertical and radial particle density distribution for all grain sizes of models without any bumps after 1 Myr of evolution.

4 Discussion

4.1 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⁻³. 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. 2018).

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₁₀(F_{1.3 mm}∕F_{3.0 mm})∕log₁₀(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. 2012, 2013, 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 tohave 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⁻⁴), 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 (Eqs. (12) and (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,

$${\Sigma }_{\text{g}}^{\prime }(r)={\Sigma }_{\text{g}}(r)\times \left[1+{\sum }_i{B}_i\right],$$(14)

with $$B_i(r)=A_i\exp \left(-\frac{{(r-r_{\text{p},i})}^2}{2w_i^2}\right).$$(15)

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.

Fig. 4 Particle density distribution of the models with a pressure bump in the outer part of the disk (dip located at 40 au) after 1 Myr of evolution. We follow the same procedure as in Fig. 2 for the barriers.

Fig. 5 Vertical and radial particle density distribution of all grain sizes of the models in which a pressure bump is located in the outer disk.

Fig. 6 As Fig. 2, but with multiple pressure bumps induced either by variations in α and different gap depths (top panels), or introduced in the gas surface density that remains constant over time (bottom panels). The growth barriers follow the same convention as in Fig. 2. In these cases, α, δr,z,t = 10−4 for all four simulations.

4.2 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 emissionby assuming a 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 consequenceof the high turbulence and radial diffusion because the high turbulence produces small grains that are easier todiffuse 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–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.

Fig. 7 Radial intensity profile of the synthetic images (previously convolved with a Gaussian with an FWHM of 0′′.04 in each case) 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 r2 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.

4.3 Substructures

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

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 theintensity 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 ice-lines (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⁻⁵) 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⁻⁴) than model 13 (α, δ_r,z = 10⁻³, and δ_t = 10⁻⁵) 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 long-lived 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.

Fig. 8 As Fig. 7, but when a pressure bump is formed in the outer part of the disk.

5 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⁻¹ 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⁻⁵−10⁻⁴. 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⁻³, 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 isincluded.

• Most of the spectral indices in our work are much higher than observed because the number of pebbles is too lowas 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 foundin 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.

Acknowledgements

We are thankful to the referee for the constructive report. P.P. acknowledges support provided by the Alexander von Humboldt Foundation in the framework of the Sofja Kovalevskaja Award endowed by the Federal Ministry of Education and Research. C.T.L. was in parts funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) as part of the Schwerpunktprogramm (SPP, Priority Program) SPP 1833 “Building a Habitable Earth”. S.M.S. has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No 714769. S.M.S. acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Ref no. FOR 2634/1. We are thankful to A. Pohl for her insightful help with the coupling of dust evolution models and RADMC-3D. We thank H. Klahr and R. Latka for discussions.

Appendix A Growth barrier estimates

In this appendix we provide analytical prescriptions for different growth barriers, which depend on the relative speed among the largest grains. We assumed mono-disperse growth for nonvanishing relative velocities for equal-sized particles, and duo-disperse growth where relative velocities of particles of the same size is zero. In the latter case, we model growth assuming that the mass-dominating species grows through collisions of particles that are N times smaller. We modeled the relative speed due toradial drift (Δv_r), azimuthal drift (Δv_φ), vertical settling (Δv_z), and turbulent motion (Δv_{turb II}) as $$\begin{array}{l}\Delta v_{\text{r}}\approx \frac{c_{\text{s}}^2}{v_{\text{K}}}\left|\frac{\partial \ln P}{\partial \ln r}\right|(1-N)\text{St},\\ \Delta v_{\phi }\approx \frac{1}{2}\frac{c_{\text{s}}^2}{v_{\text{K}}}\left|\frac{\partial \ln P}{\partial \ln r}\right|\left|\frac{1}{{\text{St}}^2+1}-\frac{1}{N^2{\text{St}}^2+1}\right|,\\ \Delta v_{\text{z}}\approx h_{\text{d}}\Omega (1-N)\text{St}=c_{\text{s}}(1-N)\frac{\text{St}}{\sqrt{1+\text{St}/{\delta }_z}},\\ \Delta v_{\text{turb}\text{II}}\approx \sqrt{3{\delta }_{\text{t}}\text{St}}{c}_{\text{s}}.\end{array}$$

These are approximations for the different contributions to the relative speed if St < 1 and if the azimuthal pressure gradient is zero. For the second turbulent regime, particles were required to have Stokes numbers larger than $1/\sqrt{\text{Re}}$, where the Reynolds number is given by $$\begin{array}{ll}\text{Re}\hfill & =\frac{{\nu }_{\text{turb}}}{{\nu }_{\text{mol}}}\approx \frac{\sqrt{\pi }{\delta }_{\text{t}}{c}_{\text{s}}}{\Omega \sqrt{2}\lambda }\hfill \\ \hfill & =\frac{\sqrt{\pi }{\delta }_{\text{t}}n{\sigma }_{{\text{H}}_2}{c}_{\text{s}}}{\Omega }\stackrel{z=0}{\approx }\frac{{\delta }_{\text{t}}{\Sigma }_{\text{g}}{\sigma }_{{\text{H}}_2}}{\sqrt{2}{m}_{\text{g}}}.\hfill \end{array}$$(A.5)

Here, we adopted the cross section of molecular hydrogen, $${\sigma }_{{\text{H}}_2}\approx 2E-5{\text{cm}}^2,$$(A.6)

and assumed that the mean molecular weight m_g is about 2.3 proton masses, which is roughly $$m_{\text{proton}}\approx 1.67E-24\text{g}.$$(A.7)

The mean free path of gas molecules following the Maxwell-Boltzmann distribution is given by $$\lambda =\frac{1}{\sqrt{2}}\frac{1}{n{\sigma }_{{\text{H}}_2}}$$(A.8)

with number density of gas molecules n. We assumed that the collisional behavior can be modeled by collision partners with a factor of N ≈ 0.5 difference in Stokes number (however, see comments at the end of this appendix). The factor that depends on the Stokes number in the expression of the the relative azimuthal drift Δv_φ needs simplifications to find simple analytical expressions for the growth barriers for which the relative speed among the largest grains is dominated by this effect. This expression can be simplified by neglecting higher order terms as follows: $$\left|\frac{1}{{\text{St}}^2+1}-\frac{1}{N^2{\text{St}}^2+1}\right|=\frac{(1-N^2){\text{St}}^2}{\cancel{N^2{\text{St}}^4}+(1+N^2){\text{St}}^2+1.}$$(A.9)

The fragmentation limit can be estimated by setting the dominating relative speed equal to the critical fragmentation speed v_f (Birnstiel et al. 2010b, 2012). For the case of relative speeds dominated by radial drift, this yields (Birnstiel et al. 2012) $${\text{St}}_{\text{frag},r}=\frac{v_{\text{f}}{v}_{\text{K}}}{c_{\text{s}}^2}{\left|\frac{\partial \ln P}{\partial \ln r}\right|}^{-1}\frac{1}{1-N}.$$(A.10)

By neglecting the $\mathcal{O}({\text{St}}^2)$ term in the denominator on the right side of Eq. (A.9), the fragmentation limit in the case of relative speeds dominated by azimuthal drift can be estimated to be $${\text{nn St}}_{\text{frag},\phi }=\sqrt{2\frac{v_{\text{f}}{v}_{\text{K}}}{c_{\text{s}}{}^2}{\left|\frac{\partial \ln P}{\partial \ln r}\right|}^{-1}\frac{1}{1-N^2}}.$$(A.11)

However, the relative velocities due to azimuthal drift of the largest grains with similar-sized grains for St < 1 are expected to be lower than those caused by relative radial drift: $${\text{relvel}}_{\phi }\approx \frac{1-N^2}{2(1-N)}{\text{Strelvel}}_r<{\text{relvel}}_r.$$(A.12)

For settling-dominated relative speeds, a quadratic equation of the following form is achieved: $${\text{St}}^2-\text{St}\frac{A}{{\delta }_z}-A=0,$$(A.13)

where $$A:=\frac{v_{\text{f}}{}^2}{c_{\text{s}}{}^2{(1-N)}^2}.$$(A.14)

The physical solution is found to be $${\text{St}}_{\text{frag},z}=\frac{A}{2{\delta }_z}+\sqrt{\frac{A^2}{4{\delta }_z^2}+A}.$$(A.15)

The fragmentation limit in the turbulence-dominated regime is given by (Birnstiel et al. 2010a, 2012) $${\text{St}}_{\text{frag},\text{turb}}=\frac{1}{3}\frac{v_{\text{f}}{}^2}{{\delta }_{\text{t}}{c}_{\text{s}}{}^2}.$$(A.16)

Because higher relative velocities lead to critical fragmentation speed at smaller sizes, the overall fragmentation limit can be estimated through $${\text{St}}_{\text{frag}}=\min \left\{{\text{ St}}_{\text{frag},\text{turb}};{\text{St}}_{\text{frag},r};{\text{St}}_{\text{frag},z}\right\}.$$(A.17)

The drift limit can be obtained by setting the growth timescale $${\tau }_{\text{grow}}=2\sqrt{2\pi }\frac{{\rho }_{\text{s}}{h}_{\text{g}}}{{\Sigma }_{\text{d}}}\frac{a}{\sqrt{1+\text{St}/{\delta }_{\text{z}}}\Delta v},$$(A.18)

and drift timescale $${\tau }_{\text{drift}}=\frac{r}{\left|v_{\text{drift}}\right|}=\frac{\cancel{{\text{St}}^2}+1}{\text{St}}\frac{r{v}_{\text{K}}}{c_{\text{s}}{}^2}{\left|\frac{\partial \ln P}{\partial \ln r}\right|}^{-1},$$(A.19)

equal to each other (Klahr & Bodenheimer 2006; Birnstiel et al. 2012). In the derivation of this growth timescale, we assumed a growth rate (particle radius per time) that is half as high as in canonical derivations from Stepinski & Valageas (1997) and Kornet et al. (2001). When the Smoluchowki equation is used, a factor of 0.5 naturally shows up in the growth rate (e.g., Brauer 2008). This growth timescale is an estimate for mono-disperse growth only. However, it is still good enough to estimate the growth barriers, which can be fine-tuned with fudge factors or through the N parameter.

With help of the Stokes number in the midplane in the Epstein-drag regime Eq. (1), the growth timescale can be written in a form independentof particle radius a: $${\tau }_{\text{grow}}=4\sqrt{\frac{2}{\pi }}\frac{h_{\text{g}}}{\Delta v}\frac{1}{Z}\frac{\text{St}}{\sqrt{1+\text{St}/{\delta }_{\text{z}}}},$$(A.20)

where Z : = Σ_d∕Σ_g is the column density ratio of all particles to gas. By neglecting second-order terms in the Stokes number, the relative speed due to radial drift can be written as $${\text{nn relvel}}_r\approx (1-N)\left|v_{\text{drift}}\right|.$$(A.21)

Setting the growth and drift timescales equal leads to the quadratic equation $${\text{St}}^2-\frac{B}{{\delta }_{\text{z}}}\text{St}-B=0,$$(A.22)

where $$B:={\left(\frac{\pi (1-N)}{4\sqrt{2\pi }}\frac{v_{\text{K}}}{c_{\text{s}}}Z\right)}^2.$$(A.23)

The physical solution yields the drift limit in the radial drift-dominated growth regime $${\text{St}}_{\text{drift},r}=\frac{B}{2{\delta }_{\text{z}}}+\sqrt{\frac{B^2}{4{\delta }_{\text{z}}{}^2}+B}.$$(A.24)

By assuming settling-dominated relative velocities, we obtain $${\text{St}}_{\text{drift},z}=\frac{\sqrt{\pi }}{4\sqrt{2}}(1-N)Z{\left(\frac{v_{\text{K}}}{c_{\text{s}}}\right)}^2{\left|\frac{\partial \ln P}{\partial \ln r}\right|}^{-1}.$$(A.25)

The case of dominating relative azimuthal drift is not interesting because of Eq. (A.12). For turbulent motion as the dominating source of relative velocities, equating the growth and drift timescale after squaring both sides yields $$\frac{{\text{St}}^2}{\cancel{{\delta }_{\text{t}}/\text{St}}+{\delta }_{\text{t}}/{\delta }_{\text{z}}}=\sqrt{\frac{3\pi }{32}}Z{\left(\frac{v_{\text{K}}}{c_{\text{s}}}\right)}^2{\left|\frac{\partial \ln P}{\partial \ln r}\right|}^{-1},$$(A.26)

which gives (for St ≫ δ_z) $${\text{St}}_{\text{drift},\text{turb}\text{II}}=\sqrt{\frac{3\pi }{32}}\sqrt{\frac{{\delta }_{\text{t}}}{{\delta }_{\text{z}}}}Z{\left(\frac{v_{\text{K}}}{c_{\text{s}}}\right)}^2{\left|\frac{\partial \ln P}{\partial \ln r}\right|}^{-1}.$$(A.27)

According to Ormel & Cuzzi (2007), particles with St < Re^−1∕2 are in the first turbulent regime. This regime can be reached in the very outer parts of circumstellar disks. It has a relative speed $$\begin{array}{ll}{\text{nn nn relvel}}_{\text{turb}\text{I}}\hfill & \approx \sqrt{{\delta }_{\text{t}}}{c}_{\text{s}}\sqrt{\frac{1-N}{1+N}}\hfill \\ \hfill & \text{hphantom}\approx \times \underset{\approx \sqrt{(1-N^2)}{\text{Re}}^{1/4}\text{St}}{\underbrace{\sqrt{\frac{{\text{St}}^2}{\text{St}+{\text{Re}}^{-1/2}}-\frac{N^2{\text{St}}^2}{N\text{St}+{\text{Re}}^{-1/2}}}}}\hfill \\ \hfill & \approx \sqrt{{\delta }_{\text{t}}}{c}_{\text{s}}(1-N){\text{Re}}^{1/4}\text{St}.\hfill \end{array}$$(A.28)

Following the same steps as above to obtain the drift limit, the same quadratic equation as Eq. (A.22) with the solution Eq. (A.24) is obtained, but for B → C with

$$C:=\frac{\pi }{32}{\left(\frac{v_{\text{K}}}{c_{\text{s}}}\right)}^4{Z}^2{\left|\frac{\partial \ln P}{\partial \ln r}\right|}^{-2}{\delta }_{\text{t}}{(1-N)}^2\sqrt{\text{Re}}.$$(A.29)

The solution is given by $${\text{St}}_{\text{drift},\text{turb}\text{I}}=\frac{C}{2{\delta }_{\text{z}}}+\sqrt{\frac{C^2}{4{\delta }_{\text{z}}{}^2}+C}.$$(A.30)

Because of the different scaling of the relative velocities, this drift limit St_{drift,turb I} fits the simulation results better for N = 0.1. The overall drift limit in the turbulence-dominated case is estimated as $${\text{nn St}}_{\text{drift},\text{turb}}=\{\begin{array}{ll}{\text{nn St}}_{\text{drift},\text{turb}\text{II}},\hfill & {\text{ St}}_{\text{drift},\text{turb}\text{I}}\ge {\text{Re}}^{-1/2}\hfill \\ {\text{nn St}}_{\text{drift},\text{turb}\text{I}},\hfill & {\text{ St}}_{\text{drift},\text{turb}\text{I}}<{\text{Re}}^{-1/2}\hfill \end{array}.$$(A.31)

Because faster growth caused by higher relative velocities leads to shorter growth timescales, the overall drift limit is given by $${\text{St}}_{\text{drift}}=\max \left\{{\text{ St}}_{\text{drift},\text{turb}};{\text{St}}_{\text{drift},r};{\text{St}}_{\text{drift},z}\right\}.$$(A.32)

Both the fragmentation limit Eq. (A.17) and the drift limit Eq. (A.32) can be converted into particle radii using the St -a relation in the Epstein regime Eq. (1).

The N parameter can potentially be different for the fragmentation limit or drift limit because this depends on the size distribution, which is different in these limits. While the size distribution in drift-limited regions is very top heavy, fragmentation-limited regions show shallower size distributions. In the most cases, N = 0.3 for fragmentation limits and N = 0.8 for drift limits mimic the evolution of the maximum grain size quite well. We summarize the values used for N in different cases in Table A.1. Especially the fragmentation limits work very well. The drift limit at late times and in the most outer regions of the disk can significantly deviate from the data of numerical simulations. This could be an effect of violated conditions for the analytical estimates (e.g., relative speeds are not dominated by one effect, or the size distribution does not allow analmost mono-disperse or duo-disperse approach).

References

All Tables

All Figures

	Fig. 1 Schematic view of turbulent velocities, vertical stirring and settling, and radial diffusion of dust. 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 (see Sect. 2.1.1).
In the text

Fig. 2 Particle density distribution after 1 Myr of evolution as a function of distance from the star and grain size of the models without any pressure bump (Table 1). The solid white line represents a Stokes number of unity. Cyan and red show the fragmentation (Eq. (A.17)) and drift (Eq. (A.32)) limit. The style of the line shows what dominates the relative velocity: turbulence (Eqs. (A.16) and (A.27)), radial drift (Eqs. (A.10), and (A.24)), and settling (Eqs. (A.15) and (A.25)). 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.

In the text

	Fig. 3 Vertical and radial particle density distribution for all grain sizes of models without any bumps after 1 Myr of evolution.
In the text

	Fig. 4 Particle density distribution of the models with a pressure bump in the outer part of the disk (dip located at 40 au) after 1 Myr of evolution. We follow the same procedure as in Fig. 2 for the barriers.
In the text

	Fig. 5 Vertical and radial particle density distribution of all grain sizes of the models in which a pressure bump is located in the outer disk.
In the text

	Fig. 6 As Fig. 2, but with multiple pressure bumps induced either by variations in α and different gap depths (top panels), or introduced in the gas surface density that remains constant over time (bottom panels). The growth barriers follow the same convention as in Fig. 2. In these cases, α, δr,z,t = 10−4 for all four simulations.
In the text

Fig. 7 Radial intensity profile of the synthetic images (previously convolved with a Gaussian with an FWHM of 0′′.04 in each case) 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 r2 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.

In the text

	Fig. 8 As Fig. 7, but when a pressure bump is formed in the outer part of the disk.
In the text