© ESO, 2012

1. Introduction

The study of planet formation is an important field in astronomy with a large amount of research having been completed since the middle of the twentieth century, although countless unanswered questions remain. One of these questions is what is the cause of the observed grain growth to mm sized particles in the outer disk regions (Beckwith & Sargent 1991; Wilner et al. 2000, 2005; Testi et al. 2001, 2003; Andrews & Williams 2005; Rodmann et al. 2006; Natta et al. 2007; Isella et al. 2009; Lommen et al. 2009; Ricci et al. 2010a, 2011; Guilloteau et al. 2011), which suggests that there is a mechanism operating to prevent the rapid inward drift (Klahr & Henning 1997; Brauer et al. 2007; Johansen et al. 2007). Different efforts aim to explain theoretically the growth from small dust particles to planetesimals, which have led to the development of different numerical models, e.g. Nakagawa et al. (1981), Dullemond & Dominik (2005), Brauer et al. (2008), Zsom & Dullemond (2008), Okuzumi (2009), Birnstiel et al. (2010a). Since circumstellar disks exhibit a wide range of temperatures, they radiate from micron wavelengths to millimeter wavelengths, which is why they can be observed with infrared and radio telescopes. With the construction of different kinds of these telescopes, e.g. Spitzer, Herschel, SMA, EVLA, or ALMA, astronomers can observe in more detail the material inside accretion disks around young stars. The parallel development of theory and observations have allowed astrophysicists to study the different stages of planet formation, making this topic one of the most active fields in astronomy today.

In the first stage of planet formation, the growth from sub-micron sized particles to larger objects is a complex process that contains many physical challenges. In the case of a smooth disk with a radial gas pressure profile that monotonically decreases with radius, the dust particles drift inwards because the gas moves with a sub-Keplerian velocity due to the gas pressure gradient. Before a large object can be formed, the radial drift causes dust pebbles to move towards the star. Moreover, the large relative velocities produced by turbulence and radial drift cause the solid particles to reach velocities that lead to fragmentation collisions that prevent dust particles from forming larger bodies (Weidenschilling 1977; Brauer et al. 2008; Birnstiel et al. 2010b). The combination of these two problems is called a “meter-size barrier” because on timescales as short as 100 years, a one-meter-sized object at 1 AU moves towards the star owing to the radial drift, preventing any larger object from forming.

The observations in the inner regions of the disk, where planets such as Earth should form, are very difficult because these regions are so small on the sky that few telescopes can spatially resolve them. In addition, these regions are optically thick. However, a meter-size barrier in the inner few AU is equivalent to a “millimeter-size barrier” in the outer regions of the disk. These outer regions ( ≳ 50 AU) are much easier to spatially resolve and are optically thin. Moreover, one can use millimeter observations, which probe precisely the grain size range of the millimeter-size. Therefore, the study of dust growth in the outer disk regions may teach us something about the formation of planets in the inner disk regions.

Observations of protoplanetary disks at sub-millimeter and mm wavelengths show that the disks remain dust-rich over several million years with large particles in the outer regions (Natta et al. 2007; Ricci et al. 2010a). However, it remains unclear how to prevent the inward drift and how to explain theoretically that mm-sized particles are observed in the outer regions of the disk. Different mechanisms of planetesimal formation have been proposed to resist the rapid inward drift, such as: gravitational instabilities (Youdin & Shu 2002), the presence of zonal flows (Johansen et al. 2009, 2011; Uribe et al. 2011), or dead zones of viscously accumulated gas that form vortices (Varnière & Tagger 2006). With the model presented here, we wish to simulate mechanisms that allow long-lived pressure inhomogeneities to develop in protoplanetary disks, by artificially adding pressure bumps to a smooth density profile.

To confront the millimeter-size barrier, it is necessary to stop the radial drift by considering a radial gas pressure profile that does not monotonically decrease with radius. We instead take a pressure profile with local maxima adding a sinusoidal perturbation of the density profile. These perturbations directly influence the pressure, following a simple equation of state for the pressure in the disk. Depending on the size of the particle, the dust grains are perfectly trapped in the pressure peaks because a positive pressure gradient can cause the dust particles to move outwards. On the other hand, turbulence can mix part of the dust particles out of the bumps, so that overall there may still be some net radial inward drift. More importantly, dust fragmentation may convert part of the large particles into micron-size dust particles, which are less easily trapped and thus drift more readily inward.

Birnstiel et al. (2010b) compared the observed fluxes and mm spectral indices of the Taurus (Ricci et al. 2010a) and Ophiucus (Ricci et al. 2010b) star-forming regions with predicted fluxes and spectral indices at mm wavelengths. They neglected the radial drift,assuming that the dust particles to remained within the outer disk regions. They attempted to keep the spectral index at low values, which implied that the dust particles had millimeter sizes (Beckwith & Sargent 1991). However, they over predicted the fluxes. As an extension of their work, we model here the combination of three processes: the radial drift, the radial turbulent mixing, and the dust coagulation/fragmentation cycle in a bumpy surface density profile. Our principal aim is to find out how the presence of pressure bumps can help us to explain the retainment of dust pebbles in the outer regions of protoplanetary disks, while still allowing for moderate drift and achieving a closer match with the observed fluxes and mm spectral indeces. In addition, we show simulated images using different antenna configurations of the complete range of operations of ALMA, to study whether it is possible to detect these inhomogeneities with future ALMA observations.

This paper is ordered as follows: Sect. 2 describes the coagulation/fragmentation model and the sinusoidal perturbation that we apply to the initial conditions of the gas surface density. Section 3 describes the results of these simulations, the comparison between existing mm observations of young and forming disks, and the results of our model. We discuss whether the type of structures generated by our model will be detectable with future ALMA observations. In Sect. 4, we explore the relation of our model to the predictions of current simulations of the magnetorotational instability (MRI) (Balbus & Hawley 1991). Finally, Sect. 5 summarizes our results and the conclusions of this work.

Fig. 1 Comparison between the gas surface density (left plot) taken in this work (Eq. (1)) for two different values of the amplitude and constant width (dashed and dot-dashed lines). The Rossby wave instability (Regaly et al. 2012), and the presence of zonal flows caused by MHD instabilities (Uribe et al. 2011). Right plot shows the pressure gradient for each of the gas surface density profiles.

2. Dust evolution model

We used the model presented in Birnstiel et al. (2010a) to calculate the evolution of the dust surface density in a gaseous disk, the radial drift, and the amount of turbulent mixing. The dust size distribution evolves owing to the grain growth, cratering, and fragmentation. We accounted for the relative velocities caused by Brownian motion, turbulence, and both radial and azimuthal drift, as well as vertical settling are taken into account.

In this work, we do not consider the viscous evolution of the gas disk because the aim is to investigate how dust evolution is influenced by the stationary perturbations of an otherwise smooth gas surface density. The effects of a time-dependent perturbation and the evolution of the gas disk will be the subject of future work. For a comprehensive description of the numerical code, we refer to Birnstiel et al. (2010a).

To simulate the radial pressure maxima that allow the trapping of particles, we consider a perturbation of the gas surface density that is assumed for simplicity to be a sinusoidal perturbation such that ${\mathrm{\Sigma }}^{\mathrm{\prime }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\left(\mathrm{1}\mathrm{+}\mathit{A}\cos \left[\mathrm{2}\mathit{\pi }\frac{\mathit{r}}{\mathit{L}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\right]\right)\mathit{,}$(1)where the unperturbed gas surface density Σ(r) is given by the self similar solution of Lynden-Bell & Pringle (1974)$\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{0}}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{c}}}\right)}^{\mathrm{-}\mathit{\gamma }}\exp \left[\mathrm{-}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{c}}}\right)}^{\mathrm{2}\mathrm{-}\mathit{\gamma }}\right]\mathit{,}$(2)where r_c is the characteristic radius, taken to be 60 AU, and γ is the viscosity power index equal to 1, which are the median values found from high angular resolution imaging in the sub-mm of disks in the Ophiucus star-forming regions (Andrews et al. 2010). The wavelength L(r) of the sinusoidal perturbation depends on the vertical disk scale-height H(r) by a factor f as $\mathit{L}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathit{fH}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{,}$(3)with H(r) = c_sΩ^-1, where the isothermal sound speed c_s is defined as ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{B}}\mathit{T}}{\mathit{\mu }{\mathit{m}}_{\mathrm{p}}}\mathit{,}$(4)and the Keplerian angular velocity Ω is $\mathrm{\Omega }\mathrm{=}\sqrt{\frac{\mathit{G}{\mathit{M}}_{\mathit{\star }}}{{\mathit{r}}^{\mathrm{3}}}}\mathit{,}$(5)where k_B is the Boltzmann constant, m_p is the mass of the proton, and μ is the mean molecular mass, which in proton mass units is taken as μ = 2.3. For an ideal gas, the pressure is defined as $\mathit{P}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{=}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{\rho }\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathit{,}$(6)where ρ(r,z) is the gas density, such that ${\mathrm{\Sigma }}^{\mathrm{\prime }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{{}^{\mathrm{\int }}}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}\mathit{\rho }\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{d}\mathit{z}$. With the surface density described by Eq. (1), we can have pressure bumps such that the wavelength increases with radius. These bumps are static, which may not be entirely realistic. However, these can be a good approximation of long-lived, azimuthally extended pressure bumps,which might for instance be the result of MHD effects (see Johansen et al. 2009; Dzyurkevich et al. 2010). The influence of time-dependent pressure perturbations (e.g. Laughlin et al. 2004; Ogihara et al. 2007) on the dust growth process will be the topic of future work. The left plot of Fig. 1 (dashed lines) shows the behavior of the perturbed surface density for two values of the amplitude and fixed value of the width. The right plot of Fig. 1 shows the corresponding pressure gradient.

The very fine dust particles move with the gas because they are strongly coupled to the gas because the stopping time is very short. In the presence of a drag force, the stopping time is defined as the time that a particle, with a certain momentum, needs to become aligned to the gas velocity. However, when the particles are large enough and they do not move with the gas, they experience a head wind, because of the sub-Keplerian velocity of the gas, hence lose angular momentum and move inwards. In this case, the resulting drift velocity of the particles is given by Weidenschilling (1977)${\mathit{u}}_{\mathrm{drift}}\mathrm{=}\frac{\mathrm{1}}{{\mathrm{St}}^{-1}\mathrm{+}\mathrm{St}}\frac{{\mathit{\partial }}_{\mathit{r}}\mathit{P}}{\mathit{\rho }\mathrm{\Omega }}\mathrm{·}$(7)Comparing Eq. (7) with the expression for the drift velocity given by Birnstiel et al. (2010a, Eq. (19)), we find that the drag term for the radial movement of the gas is not taken here since we assume a stationary gas surface density. The Stokes number, denoted by St, describes the coupling of the particles to the gas. The Stokes number is defined as the ratio of the eddy turn-over time (1/Ω) to the stopping time. For larger bodies, the Stokes number is much greater than one, which implies that the particles are unaffected by the gas drag and consequently move on Keplerian orbits. When the particles are very small, St  ≪  1, they are strongly coupled to the gas. Since the gas orbits at sub-Keplerian velocity because of its pressure support, there is a relative velocity between the dust particles and the gas. The Stokes number is equal to unity when the particles are still influenced by the gas drag but not completely coupled to the gas, instead being marginally coupled and moving at speeds between the Keplerian and the sub-Keplerian gas velocity.

The retainment of dust particles owing to pressure bumps depends on the size of the particles. Since very small particles, with St  ≪  1, are tightly coupled to the gas, they do not drift inwards. However, radial drift becomes significant when the size of the particles increases and reaches its strongest value when St = 1 (see Eq. (7)). In the Epstein regime, where the ratio of the mean free path of the gas molecules λ_mfp to the particle size, denoted by a, satisfies λ_mfp/a ≥ 4/9, the Stokes number is given by $\mathrm{St}\mathrm{=}\frac{\mathit{a}{\mathit{\rho }}_{\mathrm{s}}}{{\mathrm{\Sigma }}_{\mathrm{g}}}\frac{\mathit{\pi }}{\mathrm{2}}\mathit{,}$(8)where ρ_s is the solid density of the dust grains, which is assumed to be constant (see Table 1). In this case, particles are small enough to be in this regime. Parameterizing the radial variation in the sound speed via ${\mathit{c}}_{\mathrm{s}}\mathrm{\propto }{\mathit{r}}^{\mathrm{-}\mathit{q}\mathit{/}\mathrm{2}}\mathit{,}$(9)where for a typical disk, the temperature is assumed to be a power law such that T ∝ r ^− q, which is an approximation to the temperature profile taken for this model. Therefore, the wavelength of the perturbation L(r) scales as: $\mathit{L}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathit{fH}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{\propto }\mathit{f}{\mathit{r}}^{\mathrm{(}\mathrm{-}\mathit{q}\mathrm{+}\mathrm{3}\mathrm{)}\mathit{/}\mathrm{2}}\mathit{.}$(10)The pressure bumps have the same amplitude A and wavelength L(r) as the density, because the over-pressure is induced by adding inhomogeneities to the gas surface density and parameterizing the temperature on the midplane by a power law (Nakamoto & Nakagawa 1994). The model taken here can artificially imitate e.g. the case of zonal flows in protoplanetary disks, where over-densities create pressure bumps. Zonal flows are formed by MRI, which depend on the degree of ionization in the disk, i.e. on the temperature of the disk and other factors such as the exposure to cosmic and stellar rays. The MRI appears to be the most probable source of turbulence, and if the turbulence is not uniform, there can be excitation of long-lived pressure fluctuations in the radial direction. For instance, Johansen et al. (2009) performed shearing box simulations of a MRI turbulent disk that exhibit large scale radial variations in Maxwell stresses of 10%. Dzyurkevich et al. (2010) presented three-dimensional (3D) global non-ideal MHD simulations including a dead zone that induces pressure maxima of 20–25%. Uribe et al. (2011) presented 3D global MHD simulations, leading to pressure bumps of around 25%. However, when the viscosity drops, the gas surface density changes causing a local inversion of the pressure gradient and an accumulation of dust particles. This matter accumulation causes Rossby wave instabilities (Lovelace et al. 1999) that lead to non-axisymmetric distributions in the disk, which we cannot exactly model at this moment since our dust evolution models are axisymmetric.

To constrain the values of the amplitude and wavelength of the perturbation, we take into account three different factors. First, it is important to analyze the necessary conditions to have a local outward movement of the dust particles. Second, we compare our assumptions with current studies of zonal flows (Johansen et al. 2009; Uribe et al. 2011). And thirdly, we only work with small-enough amplitude disturbances such that the disk has an angular momentum per unit mass that increases outwards, meaning that it is Rayleigh stable.

The Rayleigh criterion establishes that for a rotating fluid system, in the absence of viscosity, the arrangement of angular momentum per unit mass is stable if and only if it increases monotonically outward (Chandrasekhar 1961), implying that $\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}\mathit{r}{\mathit{v}}_{\mathit{\phi }}\mathrm{\right)}\mathit{>}\mathrm{0.}$(11)Since turbulence is necessary to ensure angular momentum transport, instabilities may occur if a magnetic field is present, and in those cases the angular velocity decreases with radius (MRI). For a typical α-turbulent disk, the MRI scale time is much longer than the time needed by the disk to recover the Rayleigh stability, which implies that the disk should remain quasi-stable at all times (Yang & Menou 2010). Any perturbation in the gas surface density that is Rayleigh unstable will almost instantly diminish sufficiently to ensure that the gas becomes Rayleigh stable again, thereby lowering its amplitude. This happens on a scale time that is much shorter than what MRI could ever counteract. The angular velocity of Eq. (11) is given by ${\mathit{v}}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{=}{\mathit{v}}_{\mathit{k}}^{\mathrm{2}}\mathrm{+}\frac{\mathit{r}}{\mathit{\rho }}\frac{\mathit{\partial P}}{\mathit{\partial r}}\mathrm{=}{\mathit{v}}_{\mathit{k}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{2}\mathit{\eta }\mathrm{)}$(12)where $\mathit{\eta }\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}\mathit{r}{\mathrm{\Omega }}^{\mathrm{2}}\mathit{\rho }}\frac{\mathit{\partial P}}{\mathit{\partial r}}\mathrm{·}$(13)The Rayleigh stability of the disk depends on the amplitude and the width of the bumps. In this case, we wish to study the influence of the amplitude of the perturbation, such that for this analysis we constrain the value of the wavelength of the perturbation, where f equals to unity, such that it remains consistent with the values expected from predictions of zonal flows models by Uribe et al. (2011, see Fig. 1.

Taking the perturbed gas density of Eq. (1) and f = 1.0, it is possible to find the upper limit of the amplitude to satisfy Eq. (11), i.e. the condition to remain Rayleigh stable at all times. This calculation allows the maximum value of the amplitude A to be at most  ~ 35% of the unperturbed density.

Various investigations have focused on the possibility of Rayleigh instabilities when disks have sharp peaks in their radial density profiles (see e.g. Papaloizou & Lin 1984; Li et al. 2001; and Yang & Menou 2010). These profiles can exist when the temperature in the midplane is insufficient to ionize the gas (Gammie 1996) and as a result the turbulence parameter α decreases. These regions are known as “dead zones” and these are possible locations for the formation of planet embryos. In the boundaries of these regions, Varnière & Tagger (2006) showed that it is possible to have a huge vortex with a local bump in the gas surface density. Lovelace et al. (1999) demonstrated that these perturbations create an accumulation of gas that leads to the disk being unstable to Rossby wave instability (RWI). As a comparison of the amplitudes generated by RWI vortices and the amplitudes of our perturbations, the left plot of Fig. 1 also shows the azimuthally average gas surface density of a large-scale anticyclonic vortex modeled by Regaly et al. (2012). In these cases, the equilibrium of the disk is affected in such a way that the disk becomes Rayleigh unstable. Since we focus on perturbations that allow the disks to remain Rayleigh stable at all evolutionary times, we do not consider our perturbed density to be similar to this kind of amplitudes.

Since the drift velocity is given by Eq. (7), to prevent an inward drift, the value of η from Eq. (13) must be negative, implying that the pressure gradient has to be positive in some regions of the disk. Doing this calculation for the condition that η < 0, we get that the values of the amplitude A have to be at least about 10%. In right plot of Fig. 1, we see that with an amplitude of 10% the pressure gradient barely reaches positive values in the inner regions of the disk ( ≲ 50 AU). Summarizing the upper and lower values of the amplitude, we should find that 0.10 ≲ A ≲ 0.35, when the width of the perturbation is taken to be f = 1.

Taking into account the growth-fragmentation cycle and the existence of pressure bumps, the radial drift efficiency can be lower if the bumps have appropriate values of the amplitude and the length scale. When the particles grow by coagulation, they reach a certain size with velocities high enough to cause fragmentation (fragmentation barrier). The two main contributors to the of relative velocities are the radial drift and the turbulence. In the bumps, the radial drift can be zero if the pressure gradient and the azimuthal relative velocities are high enough; however there are still non-zero relative velocities owing to the turbulence. Therefore, it is necessary to have this condition, to ensure that the particles do not reach the threshold where they fragment. The maximum turbulent relative velocity between particles, with St  ~  1, is given by Ormel & Cuzzi (2007)$\mathrm{\Delta }{\mathit{u}}_{\max }^{\mathrm{2}}\mathrm{\simeq }\frac{\mathrm{3}}{\mathrm{2}}\mathit{\alpha }\mathrm{St}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{,}$(14)where Eq. (14) is off by a factor of 2 for St  ≲ 0.1. Therefore, to break through the mm size barrier, we must ensure that Δu_max is smaller than the fragmentation velocities of the particles v_f. Collision experiments using silicates (Güttler et al. 2010) and numerical simulations (Zsom et al. 2010) show that there is an intermediate regime between fragmentation and sticking, where particles should bounce. In this present study, we did not take into account this regime because there are still many open questions in this field. For example, Wada et al. (2011) suggest that there is no bouncing regime for ice particles that may be present in the outer regions of the disks (see Schegerer & Wolf 2010). Laboratory experiments and theoretical work suggest that typical fragmentation velocities are on the order of few m s^-1 for silicate dust (see e.g. Blum & Wurm 2008). Outside the snow-line, the presence of ices affects the material properties, making the fragmentation velocities increase (Schäfer et al. 2007; Wada et al. 2009). Since in this work we consider a radial range from 1 AU to 300 AU, the fragmentation threshold velocity is assumed to be v_f = 10 m s^-1. All the parameters used for this model are summarized in Table  1.

For particles with St ≲ 1, taking the size at which the turbulent relative velocities Δu_max are as high as the fragmentation velocity v_f, we can find the maximum value of the grain size, which is ${\mathit{a}}_{\max }\mathrm{\simeq }\frac{\mathrm{4}{\mathrm{\Sigma }}_{\mathrm{g}}}{\mathrm{3}\mathit{\pi \alpha }{\mathit{\rho }}_{\mathrm{s}}}\frac{{\mathit{v}}_{\mathit{f}}^{\mathrm{2}}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}$(15)where a_max is valid only for St ≲ 1 because for larger bodies the turbulent relative velocities are smaller than the given in Eq. (14) (for a detailed discussion of the turbulent relative velocities see Ormel & Cuzzi 2007).

Fig. 2 Vertically integrated dust density distribution at 1 Myr for A = 0.1 and f = 1 (top) and A = 0.1 and f = 3 (bottom). The white line shows the particle size corresponding to a Stokes number of unity, which shows the same shape as the gas surface density Σ′ of Eq. (1) (see Eq. (8)). The blue line represents the maximum size of the particles before they reach fragmentation velocities (fragmentation barrier according to Eq. (15)).

Fig. 3 Ratio between final and the initial dust mass between 50 AU and 100 AU, at different times of evolution. Taking a constant value of the amplitude A = 0.1 and different values of the width of the perturbed density (Eq. (1)). For α = 10-3, we present our results in the left plot and for α = 10-4 in the right plot.

The dust grain distribution n(r,z,a) is the number of particles per cubic centimeter per gram interval in particle mass, which depends on the grain mass, the distance to the star r, and the height above the mid-plane z, such that $\mathit{\rho }\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathit{n}\mathrm{\left(}\mathit{r,z,a}\mathrm{\right)}\mathrm{·}\mathit{m}\mathrm{d}\mathit{m}$(16)is the total dust volume density. The quantity n(r,z,a) can change because of the grain growth and distribution of masses via fragmentation. The vertically integrated dust surface density distribution per logarithmic bin is defined to be $\mathit{\sigma }\mathrm{\left(}\mathit{r,a}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}\mathit{n}\mathrm{\left(}\mathit{r,z,a}\mathrm{\right)}\mathrm{·}\mathit{m}\mathrm{·}\mathit{a}\mathrm{d}\mathit{z}$(17)and the total dust surface density is then ${\mathrm{\Sigma }}_{\mathrm{d}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathit{\sigma }\mathrm{\left(}\mathit{r,a}\mathrm{\right)}\mathrm{d}\ln \mathit{a}\mathit{.}$(18)

3. Results

In this section we describe our results of dust evolution, the comparison to current mm observations of young disks and the observational signatures of our model that we will detect using ALMA.

3.1. Density distribution of dust particles

The simulations were performed for a disk of mass 0.05   M_⊙, with a surface density described by Eq. (1) from 1.0 AU to 300 AU, around a star with one solar mass. The turbulence parameter α is taken to be 10^-3, unless another value is specified. Figure 2 shows the vertically integrated dust density distribution, taking into account coagulation, radial mixing, radial drift, and fragmentation, after 1 Myr of the evolution of the protoplanetary disk. The solid white line shows the particle size corresponding to a Stokes number of unity. From Eq. (8), we can see that when St = 1, the particle size a is proportional to the gas surface density Σ′, and the solid line then reflects the shape of the surface density. The blue line of Fig. 2 represents the fragmentation barrier, which illustrates the maximum size of the particles before they reach velocities larger than the fragmentation velocity (see Eq. (15)). Hence, particles above the fragmentation barrier should fragment to smaller particles, which again contribute to the growth process.

Both plots of Fig. 2 have the same amplitude of the sinusoidal perturbation A = 0.1. The factor f, which describes the width of the perturbation, is taken to be f = 1 in the top plot, and f = 3 in the bottom plot of Fig. 2.

This result shows the following. First, the amplitude A = 0.1 of the perturbation is not high enough to have a positive pressure gradient in those regions (see right plot of Fig. 1) such that particles can be retained in the outer regions of the disk after some Myr. Instead the dust particles are still affected by radial drift and turbulence such that the particles do not grow to larger than mm size in the outer regions.

Second, taking a greater value of the factor f, at the same amplitude, implies that the retention of particles is even weaker. This is because with a wider perturbation is it harder to have positive pressure gradient. For a smaller value of f, the pressure gradient is expected to be higher because the surface profile should be steeper, hence any dust trapping much more efficient. However, the diffusion timescale τ_ν depends quadratically on the length scale ℓ. Therefore, when the wiggles of the perturbation are assumed to have a shorter wavelength, the diffusion times become much shorter, implying that the turbulence mixing ejects the dust particles from the bumps more rapidly, even when the surface density is steeper for narrow wiggles. More precisely, τ_ν ∝ ℓ²ν^-1 where the viscosity is defined as ν = αc_sh. For this reason, we note in Fig. 3-left plot that the trapping is more effective taking values of the width less than one, but only by a very small margin. As a result, the ratio of the final to the initial mass of the dust for r ∈  [50,100]  AU remains almost constant when the width is taken to be larger than 0.3. We conclude that for an amplitude of A = 0.1, the trapping after several million years does not become more effective when the wavelength of the perturbation is assumed to be shorter.

Only when the diffusion timescales become larger or equal to the drift timescales for a given pressure profile, turbulence parameter and Stokes number, can the dust particles be retained inside the bumps and therefore grow. From Eq. (7), we can deduce that the drift timescales can be given by τ_drift ∝ ℓ(∂_ℓP)^-1, where within the bumps ℓ is the width of the perturbations (which depends directly on f). As a consequence, after an equilibrium state is reached, both the drift and diffusion timescales are proportional to the square of the width. Accordingly, for a given value of f, the effects of turbulent mixing and radial drift cancel.

We note in Fig. 3 that for f = {0.3, 0.2} there is a small effect on the retention of particles for α = 10^-3 (left plot) and a significant effect for α = 10^-4 (right plot). This implies that for these values of f and A, the pressure gradient becomes positive enough in the outer regions (r ∈  [50,100]  AU) to trap the particles. Since A and f are the same in both cases, the pressure gradient is exactly the same. However, for α = 10^-3, the low efficiency that becomes evident when f is reduced, vanishes after two Myr, because turbulent mixing and radial drift cancel each other.

Nevertheless, when α is lowered by one order of magnitude (right panel of Fig. 3), the diffusion timescales become longer. In this case, the drift timescales become shorter than the diffusion timescales, hence the ratio of the final to the initial dust mass increases on average for each f. When f is small enough for there to be positive pressure gradient (f = {0.3, 0.2}), outward drift wins over turbulent mixing, and as a result the trapping of particles is a visibly effective. However, there is almost no difference between f = 0.3 and f = 0.2. On contrast, this counterbalance effect between radial drift and turbulence when f varies does not happen when the amplitude of the perturbation changes.

Fig. 4 Particle size corresponding to a Stokes number of unity for A = 0.3 and f = 1.0 and location of the fragmentation barrier for two different values of the turbulent parameter α. The dashed line corresponds to r = 50 AU to distinguish the maximum size particle in the outer regions of the disk for each case.

We fixed the value of the width of the perturbation to unity, because this value is consistent with current model predictions of zonal flows. The comparison between our assumption of the density inhomogeneities and the work of Uribe et al. (2011) is discussed in Sect. 4. In addition, for larger values of the width, we should have higher values of the amplitudes to maintain a positive pressure gradient. In this case, however the disk easily becomes Rayleigh unstable when the amplitude increases. These are the reasons why we fix the value of the width to unity and no higher.

Simulations of MRI-active disks suggest that the typical values of the turbulence parameter α are in the range of 10^-3 − 10^-2 (Johansen & Klahr 2005; Dzyurkevich et al. 2010). In this work, we focus on the results for α = 10^-3, because for greater turbulence the viscous timescales become shorter than the growth timescales of the dust, causing the particles to escape from the bumps and then drift radially inwards before any mm-sizes are reached. In addition, if α is taken to be one order of magnitude higher, the fragmentation barrier is lower by about one order of magnitude in grain size, implying that particles do not grow to grater than mm-sizes in the outer regions of the disk (see Eq. (14)). Figure 4 shows the location of the fragmentation barrier for the case of A = 0.3 and f = 1.0 and two different values of α. We note that the maximum value of the grain size for the case of α = 10^-2 is on the order of a few mm in the outer regions of the disk r > 50 AU, while for α = 10^-3 the grains even reach cm-sizes.

Figure 5 compares the surface density distribution for two different values of the amplitude of the perturbation A = 0.1 and A = 0.3, at different evolution times. Taking A = 0.3, we note that in the pressure bumps there is a high density of dust particles, even after 5 Myr of evolution for a maximum radius around 80 0AU. For A = 0.3 and r ≳ 100 AU, there is a small amount of particles above the fragmentation barrier for the different times of evolution. We note that the line of the fragmentation barrier given by Eq. (15) was calculated by taking into account only turbulent relative velocities, since radial and azimuthal turbulence relative velocities were assumed to be zero at the peaks of the bumps. Particles with St  > 1 were no longer strongly coupled to the gas, hence they are not affected by gas turbulence, and the relative velocities produced by turbulence were smaller, implying that they had been able to grow over the fragmentation barrier. Moreover, we can see in the right plot of Fig. 1 that the pressure gradient for A = 0.3 and r ≳ 100 AU is always negative, hence for those regions the pressure bumps may not reduce the efficacy of radial drift. Therefore, the total relative velocities for r ≳ 100 AU can decrease, leading some particles (with St  > 1) to grow beyond the fragmentation barrier. Only the particles with St  < 1 and that exceed the fragmentation barrier should eventually fragment to smaller particles. For regions r ≲ 100 AU, dust particles continuously grow to mm-size particles by coagulation because the collision velocities produced by turbulence are smaller that the assumed fragmentation velocity v_f.

Fig. 5 Vertically integrated dust density distribution with a fixed value of length scale of f = 1, for A = 0.1 (left column) and A = 0.3 (right column) at different times 0.5 Myr, 1.0 Myr, 3.0 Myr, and 5.0 Myr from top to bottom, respectively. The solid white line shows the particle size corresponding to a Stokes number of unity. The blue line represents the fragmentation barrier according to Eq. (15).

As we have mentioned before, the efficiency of the dust trapping depends on the amplitude of the pressure bumps. It is expected that for higher amplitudes there is a greater trapping of particles because the pressure gradient is also more positive (see right plot of Fig. 1). Taking the perturbed density of Eq. (1), we can see in Fig. 6 that between 50 AU and 100 AU from the star, the amount of dust grows significantly from A = 0.1 to A = 0.3. From A = 0.3 to A = 0.5, there is still a considerable growth, but the rate of growth is slower. From A = 0.5 to A = 0.7, the rate remains almost constant, reaching a threshold. When the amplitude is increased, the amount of dust particles retained in the bumps increases to a limit when the dust growth stops because the particles reach the maximum possible value where they start to fragment owing to the high relative velocities.

As we explained in Sect. 2, taking a fixed value of f = 1.0, the disk remains Rayleigh stable for 0.10 ≲ A ≲ 0.35, which means that for those values of the amplitude, these types of perturbations can be explained via MRI without any need to suggest that Rayleigh instability is present at any evolution time. The amplitude of A = 0.3 is the most consistent with current MHD simulations of zonal flows Uribe et al. (2011), where pressure reaches radial fluctuations of 25% (see left plot of Fig. 1). The amplitude of A = 0.1 is more consisted with the case of the pressure fluctuation of zonal flows of the order of 10% by Johansen et al. (2009).

Figure 7 shows the radial dependence of the dust-to-gas ratio for different times of the simulation, for two values of the perturbation amplitude and wavelength without the gas inward motion. For A = 0.1 and f = 1.0 (top-left plot of Fig. 7), we can see that the dust-to-gas ratio decreased significantly with time across the whole disk. This implies that the dust particles do not grow considerably with time, which is what we expected because with this amplitude the trapping of the particles into the pressure bumps is ineffective. Consequently, owing to turbulence the dust particles collide, fragment, and become even smaller, such that they mix and the retention of these small particles, with St  < 1, becomes more difficult. Thus, the radial drift is not reduced, the particles with St  ≲ 1 drift inward. The particles that survive are those that have very small St  ≪ 1, and are strongly coupled to the gas. Hence, the dust-to-gas ratio initially decreases quickly and then becomes almost constant with time, which implies that after several Myr only the very small dust particles remain. The top-right plot of Fig. 7 shows that taking the same amplitude but a shorter wavelength, the dust-to-gas ratio has the same behavior. This confirms that when α turbulence is constant, a decrease in the wavelength leads to shorter diffusion timescales. Therefore, the trapping is not more effective even if the surface density is steeper for narrow bumps.

Fig. 6 Ratio of the final to the initial dust mass between 50 AU and 100 AU, at different times of evolution. We assumed a constant value of the width f = 1.0 and different values of the amplitude of the perturbed density (Eq. (1)).

In contrast, owing to the strong over-pressures at A = 0.3 (bottom plots of Fig. 7), the dust-to-gas ratio remains almost constant with time for r < 100 AU, oscillating radially between  ~10^-3 to  ~10^-1. This oscillating behavior, even after 5 Myr of dust evolution, is possible because the particles are retained in the bumps and grow sufficient to increase the dust-to-gas ratio inside the bumps. Only around  ~100 AU from the star does the dust-to-gas ratio decrease slowly with time. This implies that for r < 100 AU, the drift is counteracted by the positive local pressure gradient when the timescales for the growth are comparable to the disk evolution times, i. e. the viscous timescales. Changing the width of the perturbation to f = 1.0 for the left-bottom plot and f = 0.7 for the right-bottom plot of Fig. 7, has only a minor effect on the dust-to-gas ratio as explained before.

3.2. Comparison to current data of young disks in the millimeter range

We compared the models predictions of the disk fluxes at millimeter wavelengths with observational data obtained for young disks in Class II young stellar objects (YSOs).

To do this, we calculated the time-dependent flux for the disk models described in Sect. 2. For the dust emissivity, we adopted the same dust model as Ricci et al. (2010a,b, 2011) and Birnstiel et al. (2010b), i.e. spherical composite porous grains made of silicates, carbonaceous materials, and water ice, with relative abundances from Semenov et al. (2003). At each stellocentric radius in the disk, the wavelength-dependent dust emissivity was calculated by considering the grain size number density n(r,z,a) derived from the dust evolution models at that radius, as described in Sect. 2.

Fig. 7 Dust-to-gas mass ratio in the disk for different times in the disk evolution and the parameters summarized in Table 1: A = 0.1 and f = 1 (top-left); A = 0.3 and f = 1 (bottom-left); A = 0.1 and f = 0.7 (top-right) and A = 0.3 and f = 0.7 (bottom-right).

In the millimeter wavelength regime, the opacity can be approximated by a power law (Miyake & Nakagawa 1993), which means that the flux can be approximated to F_ν ∝ ν^α_mm, where α_mm is known as the spectral index. The spectral index gives us information about the size distribution of the dust in the disk. Figure 8 shows the time-dependent predicted fluxes at  ~1 mm (F_1   mm) and the spectral index between  ~1 and 3 mm (α_{1 − 3   mm}) for a disk model with f = 1 and A = 0.1 (top panel) or A = 0.3 (bottom panel). In the same plot, we present mm-data for young disks in Taurus (Ricci et al. 2010a; and Ricci, priv. comm.), Ophiucus (Ricci et al. 2010b), and Orion Nubula Cluster (Ricci et al. 2011) are also shown.

As detailed in Birnstiel et al. (2010b), the F_1mm versus α_{1 − 3   mm} plot reflects some of the main properties of the dust population in the disk outer regions, which dominate the integrated flux at these long wavelengths. In particular, the 1mm-flux density is proportional to the total dust mass contained in the outer disk. The mm-spectral index is instead related to the sizes of grains: values lower than about three are caused by emission from grains larger than about  ~1 mm, whereas values around 3.5 are due to smaller grains (Natta et al. 2007).

In the case of the disk model with f = 1 and A = 0.1, the time evolution of the predicted mm-fluxes clearly reflects the main features of the dust evolution depicted in the top plot of Fig. 8. Grains as large as a few millimeters quickly form in the disk outer regions (R ≳ 50 AU), and most of them are initially retained in those regions ( ≲ 0.5 Myr). As a consequence, the mm-spectral index of the disk is slightly lower than 3 at these early stages. However, the radial drift of mm-sized pebbles becomes soon significant and, as already described in Sect. 3, perturbations with a length-scale of f = 1 and amplitude of A = 0.1 are inefficient in retaining mm-sized particles in the outer disk. The 1 mm-flux density significantly decreases because of the loss of dust from the outer regions, especially the mm-sized grains, which are efficient emitters at these wavelengths. Given that the spectral index is a proxy for the grain size, which is also affected by radial drift, its value increases with time because of the gradual loss of mm-sized pebbles in the outer disk. In this case, the under-predicted fluxes are inconsistent with observational data for disk ages  ≳ 1 Myr, i.e. with the mean estimated ages of PMS stars in the Taurus, Ophiucus, and Orion regions.

Interestingly, a disk with perturbations in the gas surface density with a larger amplitude of A = 0.3 shows different results. In this case, the trapping of particles in the pressure bumps is efficient enough to retain most of the large pebbles formed in the outer disk (see bottom plot of Fig. 8). Since radial drift is much less efficient in this case, the predicted 1 mm-flux density is less affected than in the A = 0.1 case and, more importantly, the spectral index levels off at a value of  ~2.5. This model provides a good match to the bulk of the mm-data for disk ages of a few Myr, as seen in the bottom panel of Fig. 8.

We also note that the match between the model presented here for A = 0.3 is closer than that obtained by Birnstiel et al. (2010b). In contrast to the present work, these authors completely ignored radial drift, thus restricted particles to remain artificially in the disk outer regions. Specifically, for a disk with the same unperturbed disk structure presented here, they found a higher 1 mm-flux density than we obtained in the A = 0.3 perturbation case at a few Myr. This indicates that to interpret the measured mm-fluxes of young disks we need to incorporate in our models both radial drift and a physical mechanism acting in the disk to trap, although not completely, mm-sized particles in the outer disk.

Fig. 8 Comparison of the observed fluxes at mm-wavelengths of young disks in Taurus (red dots; from Ricci et al. 2010a; and Ricci, priv. comm.), Ophiucus (blue dots; from Ricci et al. 2010b), and Orion Nebula Cluster (green dots; from Ricci et al. 2011) star-forming regions with the predictions of the disk models at different times in the disk evolution (star symbols). Disk ages are indicated by numbers, in Myr, above the star symbols. The predicted  ~1 mm-fluxes (x-axis) and spectral indices between  ~1 mm and 3 mm (y-axis) are for the disk models presented in Sect. 2 with perturbations characterized by f = 1 and either A = 0.1 (top panel) or A = 0.3 (bottom panel). The  ~1 mm-flux densities for the Orion disks have been scaled by a factor of (420 pc/140 pc)2 to account for the different distances estimated for the Orion Nebula Cluster (~420 pc, Menten et al. 2007) and Taurus and Ophiucus star-forming regions (~140 pc, Bertout et al. 1999; Wilking et al. 2008).

Fig. 9 Disk image at 2 Myr and observing wavelength of 0.45 mm, the amplitude of the perturbation is A = 0.3 and the factor f = 1 for: disk model with parameters of Table 1 (top) and a simulated image using full configuration of ALMA (bottom) with a maximum value of baseline of around 3 km and an observing time of 4 h. The contour plots are at  { 2,4,6,8 }  for the corresponding rms value (see Table 2).

Fig. 10 Comparison between the simulated images for an observing wavelength of 1 mm and 2 Myr of evolution, using the full antenna configuration of ALMA for two different values of the amplitude of the perturbation: A = 0.1 (top) and A = 0.3 (bottom). The contour plots are at  { 2,4,6,8 } , of the corresponding rms value (Table 2).

3.3. Future observations with ALMA

The Atacama Large Millimeter/sub-millimeter Array (ALMA) will provide an increase in sensitivity and resolution to observe in more detail the structure and evolution of protoplanetary disks. With a minimum beam diameter of  ~5 mas at 900 GHz, ALMA will offer a resolution down to 2 AU for disks observed in Orion and sub-AU for disk in Taurus-Auriga (Cossins et al. 2010). Using the Common Astronomy Software Applications (CASA) ALMA simulator (version 3.2.0), we run simulations to produce realistic ALMA observations of our model using an ALMA array of 50 antennas 12 m-each.

The selection of observing mode to obtain the images was chosen to have simultaneously the most favorable values for the resolution and sensitivity that should be available with ALMA. The spatial resolution depends on the observing frequency and the maximum baseline of the array. We did not consider the largest array because for very large baselines, the sensitivity could be not enough for the regions that we need to observe. Therefore, we used different antenna arrays depending on the observing frequency to achieve the highest possible resolution with enough sensitivity. The sensitivity depends on the number of antennas, the bandwidth (which is taken as Δν = 8 GHz for continuum observations), and the total observing time that was fixed to four hours for each simulation. The sensitivity also depends on the atmospheric conditions. ALMA is located in Llano de Chajnantor Observatory, where the precipitable water vapor (pwv) varies between 0.5 mm and 2.0 mm depending on the observable frequency. For the simulations, we assumed that the value of the pwv varies with frequency (see Table 2). The synthetic images are fully consistent with the opacity dust distribution discussed in Sect. 3.2.

Figure 9 shows a comparison between the model image and a simulated ALMA image using the full configuration of ALMA with a maximum baseline of around 3 km and an observation total time of four hours. This image is for an observing wavelength of 0.45 mm (band 9 of ALMA 620–750 GHz). We note that the simulated images take into account the atmospheric conditions and the expected receiver noise based on technical information of the antennas, but the residual noise after data calibrations and its uncertainties are not considered. We note (see Fig. 9-bottom plot) that with one of the full configurations (max. baseline  ~3 km), it is possible to distinguish some ring structures because the dust has drifted considerably into the rings relative to the gas.

In Fig. 10, we note again the importance of having a high value of the amplitude. If the gas surface density of the disk is A = 0.3, then the effects will be observable with ALMA. Both images of the figure have been computed with the complete antenna configuration of ALMA for an observing wavelength of 1 mm and 2 Myr of the disk evolution. We can see that for A = 0.1 it is impossible to detect any structures around the star even considering a perfect data calibration. We note that owing to the trapping of dust particles at the peaks of the pressure bumps, the contrasts between rings in the simulated images is very clear, around  ~20–25%, while the contrast for the gas is almost unrecognized.

Fig. 11 Disk simulated images with parameters of Table 1, A = 0.3 and f = 1 at 2 Myr of the disk evolution and for observing frequency of: 667 GHz with a maximum baseline of around 3 km (top left), 454 GHz (bottom left) with a maximum baseline of around 4 km, 300 GHz (top right) with a maximum baseline of around 7 km, and 100 GHz with a maximum baseline of 16 km (bottom right). The contour plots are at  { 2,4,6,8 }  the corresponding rms value (see Table 2).

Figure 11 compares the simulated images at different observing wavelengths using different antenna configurations of ALMA. The antenna configuration is chosen by CASA depending on the expected resolution. The highest quality image was obtained at 100 GHz and a maximum baseline of 16 km (most extended ALMA configuration), where it is possible to clearly detect the most external ring structure and some internal ring structures. Nevertheless, with more compact configurations at different frequencies it is still possible to detect some structures from the presence of the pressure bumps, which allow the formation of mm-sized particles. However, it is important to take into account that the simulated images of Figs. 9–12 assume a perfect data calibration after observations, and that for long baselines and high frequencies the calibration effects become more important.

Taking the ratio of the images at two different wavelengths, we evaluated the values of the spectral index α_{1 − 3   mm}, which indicates the location of mm-sized grains. For the full configuration of ALMA and a maximum value of the baseline of 12 km for both observing frequencies, some regions with large particles are distinguished that are regions of low spectral index α_{1 − 3   mm} ≲ 3, as explained in Sect. 3.2. In Fig. 12, we present is the spectral index of the model data (top plot) and the spectral index for two simulated images at 0.45 mm and 3.0 mm, and a total observing time of four hours.

4. Approach to zonal flows predictions

We have so far assumed ad-hoc models of pressure bumps. We have not, however, considered the processes that may cause such long-lived bumps in protoplanetary disks? We now examine whether zonal flows are the origin of long-lived pressure bumps.

One possible origin of pressure bumps is MRI turbulence. Hawley et al. (1995) and Brandenburg et al. (1995) presented the first attempts to simulate the nonlinear evolution of MRI in accretions disk, taking a box as a representation of a small part of the disk. More recent simulations have been performed for a higher resolution (see e.g. Johansen et al. 2009) and a more global setup (Flock et al. 2011).

In magnetorotational instability, “zonal flows” are excited as a result of the energy transportation from the MRI unstable medium scales, to the largest scales causing an inverse cascade of magnetic energy, and creating a large-scale variation in the Maxwell stress (Johansen et al. 2009).

Different 3D MHD simulations have shown that in the presence of zonal flows, pressure bumps can appear when there are drops in magnetic pressure throughout the disk (Johansen et al. 2009; Dzyurkevich et al. 2010; Uribe et al. 2011). Nevertheless, in the simulations of yet higher resolution performed by Flock et al. (2011), pressure bumps are not formed. There is no conclusion about how and whether these pressure bumps can be created via zonal flows.

An alternative for the origin of the pressure bumps produced by MRI is the change in the degree of ionization. The disk becomes MRI active when the degree of ionization is sufficiently high for the magnetic field to be strongly coupled to the gas. Variable degrees of ionization in the disk could cause local changes in the magnetic stress, which could induce structures in the density and pressure.

Fig. 12 Spectral index α1−3   mm of the model data (top plot) and the spectral index taking two simulated images at 0.45 mm and 3.0 mm, with a time observation of four hours and using the full configuration (maximum baseline of 12 km).

Fig. 13 Top plot: ratio of the surface density at two different azimuthal angles of of the disk from zonal flows simulation of Uribe et al. (2011). The azimuthal angles are chosen such that for a specific radius, the amplitude of the pressure bump has a maximum Σmax, and a minimum Σmin. Bottom plot: azimuthal velocity with respect to the Keplerian velocity for the azimuthal and time-averaged surface density of the midplane.

The question we now wish to answer is are the pressure bumps generated by zonal flows of MRI-turbulence, strong enough to trap the dust in a similar way to the models of Sect. 3? To find this out, we consider the 3D MHD simulations of MRI-turbulent protoplanetary disks by Uribe et al. (2011). These models have a resolution of (N_r,N_θ,N_φ) = (256,128,256). We note that the MRI dynamical timescale is around one orbital period (at one AU), while the dust growth timescale is longer than 100 orbital periods. It is currently infeasible to study the dust growth process self-consistently in the full time-dependent 3D model. This would require the development of a MHD model tens to hundreds of times longer than currently achievable. Hence, the strategy is to find first a quasi-steady state of the gas surface density from MRI evolution, in which structures in the pressure survive the entire simulation (around 1000 inner disk orbits). Afterwards, to do the coagulation/fragmentation simulation of the dust in 1D taking the gas surface density for a specific azimuthal angle in the midplane from the results from MHD simulations.

Fig. 14 Vertically integrated dust density distribution after 2 Myr and 4 Myr of dust evolution; taking the azimuthal and time-average profile of the gas surface density in the midplane from the MHD simulations (see Uribe et al. 2011, Fig. 3) without a planet. The solid white line shows the particle size corresponding to a Stokes number of unity, which has the same shape as the density profile.

Top plot of Fig. 13 shows the ratio of the time-averaged surface density at two different azimuthal angles where the amplitude of a pressure bump reaches a maximum and a minimum at two specific radii. We note, that the variations in the azimuthal angle are very uniform, around  ~5%. This is why we wish to work with the azimuthally averaged density.

For our simulations, we assumed that the pressure structure survives and is stable on dust growth and evolution timescales. The lifetime of these structures as determined by global disk simulations is still uncertain, but it has been found to be on the order of 10–100 local orbits (at the radial position of the bump) (Johansen et al. 2009; Uribe et al. 2011). It is still an open question whether this behavior can be directly re-scaled to apply to the outer parts of the disk and in any case, the structures should be eventually diffused on turbulent diffusion timescales. In the future, the continuous generation and evolution of these structures should be implemented alongside the dust evolution. However, for lack of a better model of this time-dependence at this stage, we assumed these structures to be static.

Since these MHD simulations use a radial domain where r ∈  [1;10] , we rescaled this grid logarithmically, such that the gas surface density was taken from 10 AU to 100 AU, and scaled the surface density such that the total disk mass was 0.05   M_⊙ (see Fig. 1-left plot (solid-line)). Comparing the gas surface density obtained from MHD simulations with the assumed perturbed density Σ′ (Eq. (1)), we could see that the amplitude of the surface density perturbation from zonal flows was around 25% and comparable with the amplitude of 30% of Σ′. The widths of the bumps from Uribe et al. (2011) are not uniform, but our assumption of f = 1 agrees well with some of those bumps.

The bottom plot of Fig. 13 shows the azimuthal velocity with respect to the Keplerian velocity for the azimuthal and time-averaged surface density of the midplane. We note that the azimuthal velocity exceeds the Keplerian velocity for some regions of the disk. This implies that for these regions, the presence of zonal flows allows us to have a positive pressure gradient leaving dust particles to move outwards. Therefore, the peaks of the pressure bumps created by zonal flows may be regions where dust particles can reach millimeter sizes.

Figure 14 shows the vertically integrated dust density distribution after 2 Myr and 4 Myr of dust evolution (using the model of Birnstiel et al. 2010a). We note that at that time of evolution the pressure bumps caused by zonal flows are able to retain mm and cm sized particles in the outer regions of the disk. This is why a high mass disk was considered in this case, in order to simulate large grains. Around 50–60 AU, there is clearly a region with a high vertically integrated dust density distribution for mm and cm-sized particles.We note that the peak around 100 AU is a result of the boundary condition.

5. Conclusions

Theoretical models of dust evolution in protoplanetary disks show that the growth from sub-micron sized particles to larger objects is prevented basically by two phenomena: radial drift and fragmentation. Nevertheless, infrared and radio observations show that millimeter-sized particles can survive under those circumstances in the outer regions of the disks. Therefore, various theoretical efforts have focused on explaining the survival of those bodies.

We have taken into account the strong inhomogeneities expected to persist in the gas density profile e.g. zonal flows, and used the coagulation/fragmentation and disk-structure models of Birnstiel et al. (2010a), to investigate how the presence of pressure bumps can cause a reduction in the radial drift, allowing the existence of millimeter-sized grains in agreement with observations. In this work, we assumed a sinusoidal function for the gas surface density to simulate pressure bumps. The amplitude and wavelength disturbances were chosen by considering the necessary conditions to have outward angular momentum transport in an α-turbulent type disk, outward radial drift of dust, and reasonable values compared to the predictions of studies of zonal flows (Uribe et al. 2011).

The results presented here suggest tha pressure bumps with a width of the order of the disk scale-height and an amplitude of 30% of the gas surface density of the disk, provide the necessary physical conditions for the survival of larger grains in a disk with properties summarized in Table 1. Comparisons between the observed fluxes of the Taurus, Ophiucus, and Orion Nebula Cluster star-forming regions with the results of the models ratify that the effect of the radial drift decreases allowing particles to grow. Figure 8 shows how models with these kind of disturbances reproduce more closely the mm-observations than models with full or without radial drift.

In addition, we have presented a comparison between the bumpy density profile assumed in this work and 3D MHD models of zonal flows that can cause long lived bumps in protoplanetary disks. We have shown that the pressure bumps produced by the zonal flows of Uribe et al. (2011) agree with the amplitudes and wavelengths used in this work. Therefore, considering these bumps, the survival of dust particles is possible in the outer regions after some Myr.

The simulated images using the CASA ALMA simulator (version 3.2.0) show that, with a different antenna configuration of the final ALMA stage, the ring structures, because of the pressure bumps, should be detectable. Future ALMA observations will have an important impact on our understanding the first stages of planet formation and will be very important in investigating if the grain growth and retetion can be explained by the presence of these kind of inhomogeneities in the gas density profile.

Acknowledgments

We acknowledge Francesco Trotta for his help with the code that we used in this work to derive the mm-fluxes. We would like to thank the referee, Wladimir Lyra, for his useful suggestions. This work was supported in part through the third funding line of German Excellence Initiative. T. Birnstiel acknowledges ESO Office for Science that provided funding for the visits to Garching. L. Ricci acknowledges the Ph.D. fellowship of the International Max-Planck-Research School. A. L. Uribe acknowledges the CPU time for running the simulations in the Bluegene/P supercomputer and the THEO cluster at the Rechenzentrum Garching (RZG) of the Max Planck Society. Finally, L. Testi acknowledges the allocation of an ASI contract to the INAF-Osservatorio Astrofisico di Arcetri.

References

All Tables

All Figures

	Fig. 1 Comparison between the gas surface density (left plot) taken in this work (Eq. (1)) for two different values of the amplitude and constant width (dashed and dot-dashed lines). The Rossby wave instability (Regaly et al. 2012), and the presence of zonal flows caused by MHD instabilities (Uribe et al. 2011). Right plot shows the pressure gradient for each of the gas surface density profiles.
In the text

Fig. 2 Vertically integrated dust density distribution at 1 Myr for A = 0.1 and f = 1 (top) and A = 0.1 and f = 3 (bottom). The white line shows the particle size corresponding to a Stokes number of unity, which shows the same shape as the gas surface density Σ′ of Eq. (1) (see Eq. (8)). The blue line represents the maximum size of the particles before they reach fragmentation velocities (fragmentation barrier according to Eq. (15)).

In the text

	Fig. 3 Ratio between final and the initial dust mass between 50 AU and 100 AU, at different times of evolution. Taking a constant value of the amplitude A = 0.1 and different values of the width of the perturbed density (Eq. (1)). For α = 10-3, we present our results in the left plot and for α = 10-4 in the right plot.
In the text

	Fig. 4 Particle size corresponding to a Stokes number of unity for A = 0.3 and f = 1.0 and location of the fragmentation barrier for two different values of the turbulent parameter α. The dashed line corresponds to r = 50 AU to distinguish the maximum size particle in the outer regions of the disk for each case.
In the text

Fig. 5 Vertically integrated dust density distribution with a fixed value of length scale of f = 1, for A = 0.1 (left column) and A = 0.3 (right column) at different times 0.5 Myr, 1.0 Myr, 3.0 Myr, and 5.0 Myr from top to bottom, respectively. The solid white line shows the particle size corresponding to a Stokes number of unity. The blue line represents the fragmentation barrier according to Eq. (15).

In the text

	Fig. 6 Ratio of the final to the initial dust mass between 50 AU and 100 AU, at different times of evolution. We assumed a constant value of the width f = 1.0 and different values of the amplitude of the perturbed density (Eq. (1)).
In the text

	Fig. 7 Dust-to-gas mass ratio in the disk for different times in the disk evolution and the parameters summarized in Table 1: A = 0.1 and f = 1 (top-left); A = 0.3 and f = 1 (bottom-left); A = 0.1 and f = 0.7 (top-right) and A = 0.3 and f = 0.7 (bottom-right).
In the text

Fig. 8 Comparison of the observed fluxes at mm-wavelengths of young disks in Taurus (red dots; from Ricci et al. 2010a; and Ricci, priv. comm.), Ophiucus (blue dots; from Ricci et al. 2010b), and Orion Nebula Cluster (green dots; from Ricci et al. 2011) star-forming regions with the predictions of the disk models at different times in the disk evolution (star symbols). Disk ages are indicated by numbers, in Myr, above the star symbols. The predicted  ~1 mm-fluxes (x-axis) and spectral indices between  ~1 mm and 3 mm (y-axis) are for the disk models presented in Sect. 2 with perturbations characterized by f = 1 and either A = 0.1 (top panel) or A = 0.3 (bottom panel). The  ~1 mm-flux densities for the Orion disks have been scaled by a factor of (420 pc/140 pc)2 to account for the different distances estimated for the Orion Nebula Cluster (~420 pc, Menten et al. 2007) and Taurus and Ophiucus star-forming regions (~140 pc, Bertout et al. 1999; Wilking et al. 2008).

In the text

Fig. 9 Disk image at 2 Myr and observing wavelength of 0.45 mm, the amplitude of the perturbation is A = 0.3 and the factor f = 1 for: disk model with parameters of Table 1 (top) and a simulated image using full configuration of ALMA (bottom) with a maximum value of baseline of around 3 km and an observing time of 4 h. The contour plots are at  { 2,4,6,8 }  for the corresponding rms value (see Table 2).

In the text

	Fig. 10 Comparison between the simulated images for an observing wavelength of 1 mm and 2 Myr of evolution, using the full antenna configuration of ALMA for two different values of the amplitude of the perturbation: A = 0.1 (top) and A = 0.3 (bottom). The contour plots are at  { 2,4,6,8 } , of the corresponding rms value (Table 2).
In the text

Fig. 11 Disk simulated images with parameters of Table 1, A = 0.3 and f = 1 at 2 Myr of the disk evolution and for observing frequency of: 667 GHz with a maximum baseline of around 3 km (top left), 454 GHz (bottom left) with a maximum baseline of around 4 km, 300 GHz (top right) with a maximum baseline of around 7 km, and 100 GHz with a maximum baseline of 16 km (bottom right). The contour plots are at  { 2,4,6,8 }  the corresponding rms value (see Table 2).

In the text

	Fig. 12 Spectral index α1−3   mm of the model data (top plot) and the spectral index taking two simulated images at 0.45 mm and 3.0 mm, with a time observation of four hours and using the full configuration (maximum baseline of 12 km).
In the text

Fig. 13 Top plot: ratio of the surface density at two different azimuthal angles of of the disk from zonal flows simulation of Uribe et al. (2011). The azimuthal angles are chosen such that for a specific radius, the amplitude of the pressure bump has a maximum Σmax, and a minimum Σmin. Bottom plot: azimuthal velocity with respect to the Keplerian velocity for the azimuthal and time-averaged surface density of the midplane.

In the text

	Fig. 14 Vertically integrated dust density distribution after 2 Myr and 4 Myr of dust evolution; taking the azimuthal and time-average profile of the gas surface density in the midplane from the MHD simulations (see Uribe et al. 2011, Fig. 3) without a planet. The solid white line shows the particle size corresponding to a Stokes number of unity, which has the same shape as the density profile.
In the text