A simple model for the evolution of the dust population in protoplanetary disks
Key Words.:
accretion, accretion disks – protoplanetary disks – stars: premainsequence, circumstellar matter – planets and satellites: formationAbstract
Context:The global size and spatial distribution of dust is an important ingredient in the structure and evolution of protoplanetary disks and in the formation of larger bodies, such as planetesimals.
Aims:We aim to derive simple equations that explain the global evolution of the dust surface density profile and the upper limit of the grain size distribution and which can readily be used for further modeling or for interpreting of observational data.
Methods:We have developed a simple model that follows the upper end of the dust size distribution and the evolution of the dust surface density profile. This model is calibrated with stateoftheart simulations of dust evolution, which treat dust growth, fragmentation, and transport in viscously evolving gas disks.
Results:We find very good agreement between the full dustevolution code and the toy model presented in this paper. We derive analytical profiles that describe the dusttogas ratios and the dust surface density profiles well in protoplanetary disks, as well as the radial flux by solid material “rain out”, which is crucial for triggering any gravity assisted formation of planetesimals.We show that fragmentation is the dominating effect in the inner regions of the disk leading to a dust surface density exponent of , while the outer regions at later times can become driftdominated, yielding a dust surface density exponent of . Our results show that radial drift is not efficient in fragmenting dust grains. This supports the theory that small dust grains are resupplied by fragmentation due to the turbulent state of the disk.
Conclusions:
1 Introduction
The dust content is fundamental for many aspects of the evolution, structure, and observations of protoplanetary disks. Not only is this the material out of which terrestrial planets and the cores of giant planets form, but it also dominates the opacity, thereby determining the temperature structure and the observational appearance of protoplanetary disks. The grain surfaces also provide the ground for chemical surface reactions (e.g., Aikawa & Nomura 2006) and electron capture, so that dust strongly influences the ionization state and possibly the turbulent state of the disk (Sano et al. 2000; Bai & Goodman 2009; Turner et al. 2010).
While there have been studies of gas line emission (e.g., Dutrey et al. 2007; Panić et al. 2008), most disk mass estimates are based on observations of dust emission, assuming a constant dusttogas ratio (see Andrews & Williams 2005, 2007). However, it is well known that the dusttogas ratio in disks should be an evolving quantity (e.g., Weidenschilling 1977a; Nakagawa et al. 1981; Keller & Gail 2004; Brauer et al. 2007; Ciesla 2009). The problem with our understanding of dust transport is that it depends on the structure of the disk, which is uncertain, and also on the size distribution of the grains, which is also poorly understood.
Several authors have investigated simplified models of dust transport and size evolution: for example, Garaud (2007), Kornet et al. (2001), Rózyczka et al. (2004), Klahr & Bodenheimer (2006), or Hughes & Armitage (2010). A more complete picture is given by global simulations that selfconsistently treat both the transport and collisional evolution, such as the works by Weidenschilling (1997), Brauer et al. (2008b), Birnstiel et al. (2010), or Okuzumi et al. (2011). The last models (which we will call full models as opposed to the simplified models, such as the one presented in this paper), however, are often so complicated that first it is difficult to derive a physical understanding from them, and second they are computationally expensive.
We therefore put forward a model that is motivated by and based on the full simulations mentioned before (in this case, Birnstiel et al. 2010); however, it is very simplified so easier to understand, easier to implement, and computationally much less intensive. The aim of this work is to identify a few equations that govern the spatial evolution of dust as found in the full simulations. We derive upper limits for the size distribution from which we can estimate the drift velocity of the largest particles. The size distribution is then represented by just two grain sizes, one size representing the small dust that is coupled to the gas and another size representing the largest particles, which can be drifting inwards at a much higher velocity. While mathematically similar, the basic approach of this work is thus quite different from other simplified models (e.g. Kornet et al. 2001; Garaud 2007): we do not estimate the outcome of the physical processes, but instead identify the important concepts in the full models to construct a simplified model that is calibrated to reproduce the outcome of the full model.
Klahr & Bodenheimer (2006) have estimated how much solid material can be concentrated in decimetersized radially drifting particles in local dust traps. Once these overdensities reach Toomre unstable values and exceed the local Roche density, planetesimal formation via self gravity in the dust layer is the consequence (Johansen et al. 2006). The latter mechanism has so far only be demonstrated for ad hoc local dust to gas ratios and size distributions. What is missing is a dust evolution model, as the one presented in this work, which predicts a local dust flux and the size distribution for models of gravoturbulent planetesimal formation (Klahr & Johansen 2008).
The structure of this paper is as follows: in Section 2, we introduce a few concepts and formulae that are instrumental for the rest of this work. Section 3 presents some fiducial simulation results of the complete model that are used in Section 4 to derive and test a simple model of dust transport. Section 5 makes use of this model to discuss and explain some of the implication of it, such as the observational significance of the model, which profiles of the dust surface density or the dusttogas ratio are to be expected under certain physical conditions, as well as the resulting dust mass flux. Our findings are summarized in Section 6.
2 Background
Dust particles in protoplanetary disks are subject to radial drift (which depends on the pressure structure of the disk, see Weidenschilling 1977a), gas drag (i.e., particles following the gas accretion flow) and turbulent mixing. These effects depend on the size of the particles, which is also an evolving quantity because grains grow by sticking collisions. However there are several mechanisms that can limit collisional growth of dust particles: charging effects (Okuzumi 2009), bouncing (Güttler et al. 2010; Zsom et al. 2010), fragmentation and radial drift (Weidenschilling 1977a; Brauer et al. 2008a; Birnstiel et al. 2010).
The charging and bouncing barriers are very much dependent on the material, porosity, or fractal dimension of the grains, and are still strongly debated topics (e.g., Wada et al. 2011). In one way or another these barriers need to be overcome in order to explain the existence of larger bodies such as asteroids or planets. The fact that disks are observed to be rich in small dust particles (for a recent review, see Williams & Cieza 2011) tells us that fragmentation must be effective in wide regions of the disk, otherwise coagulation would quickly render the disk optically thin even in the infrared wavelength regime (see Dullemond & Dominik 2004; Birnstiel et al. 2009). We take this observational fact as motivation for only considering fragmentation/cratering and radial drift as growthlimiting mechanisms (Dominik & Dullemond 2008). For describing both of these effects, it is convenient to define the Stokes number St, which is the ratio of the stopping time of a particle and the turnover time of the largest eddy in the turbulent environment,
(1) 
which under typical assumptions (Epstein drag law, compact spherical particles, isothermal gas density profile, and , see Cuzzi et al. 2001) simplifies for particles near the disk midplane to
(2) 
Here is the radius and the internal density of the dust aggregate (taken to be 1.6 g cm) and the gas surface density. The Stokes number is a dimensionless quantity that describes the aerodynamic properties of a dust particle and thus its coupling to the gas flow, i.e. particles with different sizes, densities, or shapes but the same Stokes number are aerodynamically identical.
Fragmentation of dust particles stops further growth because relative velocities tend to increase with the Stokes number of the particle (for ). They can reach values well above the fragmentation threshold velocity (Brauer et al. 2008a). Lab experiments measured threshold velocities of around 1 m s for the onset of fragmentation of silicate dust grains (e.g., Blum & Wurm 2008). It has been suspected that icy particles fragment only at higher velocities due to the 10 times higher surface energies (Wada et al. 2008; Gundlach et al. 2011). Recent numerical studies of Wada et al. (2009) found that the fragmentation velocity for 10m sized icy aggregates could be as high as 50 m s. However, experiments with silicate dust grains find a fragmentation velocity that decreases with grain size (Beitz et al. 2011), which can at least partly be attributed to the growing importance of inhomogeneities in larger grains (Geretshauser et al. 2011). In the following we discuss mostly results for a threshold velocity of 10 m s, since we are interested in the global evolution of dust, i.e. the regions beyond the ice line. The Stokes number, at which the typical impact velocity of similarsized grains equals the fragmentation threshold velocity can be approximated by (see Birnstiel et al. 2009)
(3) 
where is the turbulence parameter (Shakura & Sunyaev 1973) and is the isothermal sound speed and we assume . The value of is still largely unknown, but observations of accretion signatures or disk lifetimes suggest values between and . Theoretical models of the magnetorotational instability tend to reproduce these ranges of values (e.g. Johansen & Klahr 2005; Dzyurkevich et al. 2010), but global simulations also suggest a radial dependence of (see Sano et al. 2000; Flock et al. 2011). The simulations shown in this paper were carried out with radially constant values of , but the model introduced in this work is able to treat also radial variations of .
Apart from fragmentation, radial drift can also be a limiting factor for grain growth. Radial inward drift of dust particles occurs because the head wind induced by the subKeplerian gas disk removes angular momentum from the dust particles. The drift velocity (Weidenschilling 1977a)
(4) 
thus depends on the gas pressure gradient since is the difference between the orbital velocity of the dust grain and the gas (i.e., the head wind velocity), which is given by (see Weidenschilling 1977a; Nakagawa et al. 1986)
(5) 
Here denotes the gas pressure, the midplane gas density, and the Keplerian angular velocity. It is straightforward to also include reduction factors of the radial drift speed in Eq. 4, as found in Johansen et al. (2006), however this is not the subject of this work.
3 Typical simulation results
In this section we review some typical simulation results of the complete model that treats dust growth, fragmentation and cratering as well as radial transport due to gas drag, radial drift and turbulent mixing. These results are used in Section 4 to derive a simple model that reproduces the overall evolution of the dust surface density. The gas surface density is viscously evolving using the typical turbulent viscosity prescription (Shakura & Sunyaev 1973). The temperature is estimated taking passive irradiation and viscous heating into account. A detailed description of the model can be found in Birnstiel et al. (2010). Figures 1 and 2 show snapshots from simulations of a 0.1 disk around a solar mass star. The contours denote , the dust surface density distribution per logarithmic size bin, which is defined such that
(6) 
is the total dust surface density. Throughout this paper we define as the cylindrical distance to the central star. The initial gas surface density profile for the simulations is close to the selfsimilar solutions from LyndenBell & Pringle (1974)
(7) 
for a viscosity powerlaw index of unity and a characteristic radius of 60 AU (see Fig. 3). The only significant deviation is in the inner parts of the disk where the viscous heating leads to a steeper temperature profile and consequently to a flatter surface density profile. Outside of about 3 AU, there is no deviation % between the simulations and the selfsimilar solution. The evolution of the gas surface density for both simulations is shown in Fig. 3.
Figures 1 and 2 differ only in the turbulence parameter , which is in the former and in the latter case. There are a few fundamental features that are important for the further understanding of this work:

in the left panels of Fig. 1 and 2 it can be seen that grains in the inner regions quickly grow and reach a maximum size, the socalled fragmentation barrier, which is discussed in Section 4.1.1. The size limit set by fragmentation is shown as black solid lines
^{1} . A steady state is reached in which fragmentation and coagulation balance each other. We call these regions fragmentation limited. As we go to larger radii in the first panel, we see that grains at about 20 AU and beyond are still at the state of initial growth, as discussed in Section 4.3. In the central panel of Fig. 2 we see a similar behavior, but growth does happen on slightly shorter time scales because the growth time scale of small particles () is shorter due to the higher turbulent velocities. It is only for larger particles that the growth time scale does not depend on the turbulence strength (Brauer et al. 2008a). 
In the central panels of Figs. 1 and 2, we can see that grains also in the outer regions have grown significantly. The size distribution in the inner regions are basically unchanged compared to earlier times since they are still in coagulationfragmentation equilibrium. For reasons which are discussed later, the barriers to growth, shown as dashed and solid black lines in Figs. 1 and 2 also evolve with time. Most significant is the fact that the drift limit (dashed lines, discussed in Section 4.1.2) is moving down to smaller sizes.

In the outer regions at later times, the drift barrier (dashed line) has moved below the fragmentation barrier (solid line) as can be seen in the central and right panels of Fig. 1. This means that dust particles in those regions are no longer limited by fragmentation, but instead they are more rapidly drifting away than being replenished by growth from smaller sizes. We call those regions drift limited.
4 The twopopulation model
In this section, we develop an analytical understanding of dust evolution in size and spatial dimension. The simulation results of the comprehensive model, which have been introduced in the previous section are used to derive simple equations that characterize the size distribution. This enables us to develop a simple model that can reproduce the evolution of the dust surface density of the fullscale simulations well. While the following section comprehensively explains the derivation of the model, Appendix B summarizes the model in a concise algorithmic form.
4.1 Size limits
We crudely simplify the grain size distribution by only two representative sizes, and , where is the monomer size and is the representative size of the large grain population (see also the sketch in Fig 4). While is taken to be constant in time and space, increases with time as particles grow and the disk structure evolves. This growth is not without limits. As explained in Section 2, we consider the two most important barriers towards further growth, namely grain fragmentation and radial drift. We find that can also decrease with time if the growth limits drop to smaller sizes.
Turbulent fragmentation
In the case where the fragmentation of dust grains due to turbulent relative motion is the limiting effect, analytical models of the grain size distribution have been derived in Birnstiel et al. (2011). In these stationary and local grain size distributions, a significant fraction of the dust mass is concentrated in the largest sizes, slightly below the limiting Stokes number . We parametrize this offset by an orderofunity constant such that the representative size for the largest grains in a fragmentationdominated distribution can be stated as (using Eqns. 1 and 3)
(8) 
Radial drift
Radial drift can also set a limit to the local size of dust grains if the time within which particles are removed by drift is similar or less than the time scale at which the particles are replenished, i.e. the growth time scale. To derive the growth time scale
(9) 
we follow Brauer et al. (2008a) by using the growth rate of monodisperse coagulation
(10) 
and the approximate relative velocities between the similarsized dust grains for turbulent motion (see Ormel & Cuzzi 2007)
(11) 
Assuming , the Epstein drag law, and a scale height of the dust distribution of (Youdin & Lithwick 2007)
(12) 
the growth time scale can be written as
(13) 
where is the vertically integrated dusttogas mass ratio and we neglected a factor of order unity.
It should be noted that in this case (i.e. turbulent relative velocities), the growth time scale is independent of the particle size because the rate of change of the particle radius becomes directly proportional to the particle radius itself. The drift time scale
(14) 
depends on the drift velocity (cf. Eq. 4) and thus the Stokes number of the particle. It is given by
(15) 
where is the Keplerian velocity and
(16) 
is the absolute value of the powerlaw index of the gas pressure profile. By equating the growth time scale and the size dependent drift time scale, we derive an estimate of the maximum Stokes number that can be reached before radial drift removes the dust particles,
(17) 
or in terms of particle size (again assuming Epstein drag law),
(18) 
This boundary is not as sharp as the fragmentation barrier, as can be seen in Fig. 1. We have introduced another orderofunity factor , which we calibrate later with the detailed numerical simulations.
Fragmentation by relative drift velocities
So far we have only considered fragmentation by turbulent relative velocities. However in the case of very low turbulent velocities, one might expect that the relative drift speed of the particles, either in the vertical or in the radial dimension, is still high enough to cause fragmentation. Dust grains in the tenuous disk atmosphere settle towards the midplane with velocities that can exceed the radial ones (Nakagawa et al. 1986). In a turbulent disk, mixing counteracts this downward motion and an equilibrium state is reached on time scales less than years, which means that the vertical motion is not significant throughout the lifetime of the disk. In contrast to this, radial drift motion does not vanish towards the midplane and is therefore a systematic velocity that stays relevant throughout the lifetime of the disk. We therefore only focus on the fragmenting effects of the radial dust motion. We show in Section 5.1 that this effect is irrelevant for the simulations of a turbulent disk. For laminar disks, this effect can become relevant.
To derive a size limit for this case, we need an estimate of the relative velocities induced by radial drift. From Eq. 4, we can derive the relative drift speed between two particles for Stokes numbers below unity
(19) 
Expressing the Stokes number of the smaller particle as , and equating Eq. 19 with the fragmentation threshold velocity , we derive
(20) 
which is the largest Stokes number that can be reached before the relative drift speed exceeds the fragmentation threshold. The most important collision types determining the mass gain/loss of the largest bodies is typically collisions with similar sized bodies. is therefore a reasonable assumption.
4.2 Mass fractions & transport equation
In order to represent the size distribution of the dust with only two representative sizes, we also need to know the ratio of mass contained in large and small particles. We represent the mass ratio between the two populations by another yet unknown factor ,
(21)  
(22) 
where is the surface density in the small grains while is the surface density in large grains. This representation is schematically shown in Fig. 4. In the case of a dust size distribution which is in coagulation/fragmentation equilibrium, small dust is replenished. In the case of a driftlimited size distribution, fragmentation is not efficient enough since particles drift inwards before reaching sizes at which they fragment.
For particles with a Stokes number below unity, this allows us to evolve the dust distribution by a single advectiondiffusion equation
(23) 
using the mass weighted velocity
(24) 
and the gas diffusivity . Here, and denote the radial velocity of the dust sizes and , respectively. The radial velocity of dust grains is the sum of the drift velocity given by Eq. 4 and the gas drag induced velocity
(25) 
where is the gas radial velocity. A more thorough justification of Eq. 23 is given in Appendix A.
4.3 Particle growth
The mass flux and thus the whole dust evolution is strongly influenced by the growth time scale. This is due to the fact that the dusttogas coupling and thus the radial velocity strongly depends on the particle size. Thus, particles in the outer parts, which grow less quickly, start to drift at later times, once they have reached sizes large enough to be influenced by drift. From the growth time scale (i.e. the sizedoubling time) given in Eq. 13 we can estimate the time it takes to grow from the smallest size to a maximum size
(26) 
where
(27) 
In order to take this delayed drift into account, we define a time dependence of the larger representative size
(28) 
This is an oversimplified estimate since different growth mechanisms are at play at different sizes, however we found a good agreement between our estimate and the simulation result. This can be seen in Figs. 1 and 2, where the dashdotted curve represents the time dependent upper size limit according to Eq. 28. Only at early times (left panels of Figs. 1 and 2), differs from the fragmentation or drift barrier.
It can also be seen that this estimate is slightly too low, i.e. it is slightly below the upper size of the simulated dust distribution at all radii AU in the left panel of Fig. 2. This is due to the fact that the growth time scale of Eq. 13 implicitly assumes growth to be driven only by equal sized collisions from turbulent motion according to Eq. 11. However, the initial stages of particle growth involve different turbulent regimes (see Ormel & Cuzzi 2007), which lead to faster growth for higher values of .
parameter  drift limited  fragmentation 

case  limited case  
–  0.37  
0.55  –  
0.97  0.75 
4.4 Calibration and test cases
In the semianalytic model described above, we have introduced three factors , , and . These factors are necessary since our model is not derived in a rigorously analytical way, but includes several assumptions and estimates. In order to still achieve a good agreement with the detailed numerical simulations, we need to calibrate these orderofunity factors by comparison to a grid of the full simulation results.
In order to do this, we divide the grain size distribution at each radius in two populations: the large population, which are all grains whose drift velocity is higher than the gas drag induced velocity, while everything that is smaller belongs to the small population. The gas radial velocity can be written as (see Takeuchi & Lin 2002)
(29) 
where is defined as and as , and is taken to be constant. Thus, the gas drag velocity is lower than the drift velocity for particles with a Stokes number larger than
(30) 
which for typical values (, , ) yields .
As described in Section 4.1, the fragmentation is a rather strict upper limit to the size distribution. The dust mass flow is therefore not dominated by grain sizes at the fragmentation barrier but slightly below it. We therefore calculated a fluxaveraged grain size for the large population
(31) 
where
(32) 
is the grain size corresponding to a Stokes number of and is the total radial velocity of grain size . We compared the outcome of the simplified model to a grid of 39 full simulations.
For the case of a drift limited size distribution, the size limit is not as strict as in the fragmentation limited case, i.e. there can exist grains slightly larger than . The same method as above yields therefore a value of of 0.55, which is higher than in the case of a fragmentation limited size distribution. These values of and should thus be universal, as long as the collision outcome model is not changed.
The dashed and solid lines in Fig. 1 denote the representative grain size for a drift limited and a fragmentation limited maximum grain size, respectively. It can be seen that these sizes, together with the time dependence introduced in Section 4.3, describe the limits of the grain size distribution well.
In order to properly capture the radial transport of dust with this simple model, we need also to calibrate the parameter , which describes how the mass is distributed between the two populations. Typical size distributions can be seen in the top panel of Fig. 5, denoted by the solid lines. Both the drift and the fragmentation limited distributions have most of the mass in the large grain population. However fragmentation limited distributions (e.g., the solid black line in the top panel of Fig. 5) have a larger fraction of mass in the small sizes compared to drift limited distributions (solid gray lines). The value of thus depends on whether drift or fragmentation is the limiting factor. We found typical values of around 0.75 for the fragmentation limited case and around 0.97 for the drift limited case. All values are summarized in Table 1.
It should be noted that the flux is typically determined by the large population. This is shown in Fig. 6. The plotted quantity is defined such that integration over (or equally ) yields the total dust flux. This means that a onepopulation model, which neglects the small population can describe the dust evolution in the drift limited case well. However in the case of a fragmentationdominated distribution or an outwardspreading gas disk, the small population can contain enough mass or have a significant velocity to be relevant and generally better agreement is obtained if this population is taken into account, while the additional complications and computational costs of the twopopulation model are minimal.
Another important test case is whether the total dust mass flux is reproduced. In Fig. 7 we compare the total dust flux from the full simulations (solid line) to our twopopulation model (dashed line). The total dust flux is calculated as
(33) 
The dashdotted line in Fig. 7 represents the gas accretion rate multiplied by the initial dusttogas ratio, i.e. the dust mass flux that would be found if all dust particles were perfectly coupled to the gas accretion flow. The error of the toy model is always within a factor of two, which is a surprisingly good result considering the simplicity of the model. The results are roughly consistent with the orderofmagnitude estimates of Klahr & Bodenheimer (2006).
Finally, Fig. 8 compares the profiles of the dusttogas ratio. It can be seen that the twopopulation model reproduces the full simulation results reasonabl well.
5 Discussion
The twopopulation model presented in Section 4 was shown to reproduce the time evolution of several important quantities, such as the dusttogas ratio, the dust mass flux and also the size of the drifting particles. In this section we discuss some of the applications of this model.
5.1 Fragmentation by drift velocities
In the following subsection we investigate whether the relative velocities due to radial drift of the particles can become high enough to play a role. For this to be the case, the size limit of Eq. 20 must be smaller than the limits set by both turbulent fragmentation (cf., Eq. 3) and by transport (cf., Eq. 17). These two condition can be rewritten as
(34) 
and
(35) 
respectively. Considering a typical disk where and assuming , we find that for strong turbulence (), the first condition is never fulfilled even for high fragmentation velocities of m s. For weaker turbulence such as , the regions beyond AU can fulfill Eq. 34. Eq. 35 sets constraints from the opposite direction, such that the region between about 10 and 60 AU could be influenced by drift induced fragmentation.
To quantify the error we do by considering only the size limits by drift and fragmentation, we consider the ratio between and , which is given by
(36) 
It is important to note that this deviation only needs to be considered in regions where it is larger than unity and turbulent fragmentation does not play a role. Since the regions further inside are dominated by turbulent fragmentation, the largest error we do is overestimating the upper end of the distribution by a factor of about 2.5. As we go to 60 AU the error is decreasing while the regions beyond 60 AU are in any case limited by particle drift and not by fragmentation. Decreasing the fragmentation threshold velocity seems to increase the error, however in this case the region where drift induced fragmentation applies, then disappears because turbulent fragmentation becomes more important.
Most importantly, we see that Eq. 36 depends on the dusttogas ratio. Since the region beyond 10 AU is strongly drifting, the dusttogas ratio is quickly decreasing as can be seen in Fig. 8. Therefore, we can safely ignore the size limit set by Eq. 20 if we are only concerned about the evolution of the dust surface density or the upper end of the dust size distribution.
We also carried out simulations with a very small amount of turbulence (). As expected, particles grow to large Stokes numbers until fragmentation due to radial drift sets in. However even in this case, the dusttogas ratio drops so quickly that after 1 Myr the drift limit (Eq. 17) becomes everywhere smaller than the limit set by drift induced fragmentation (Eq. 20), such that grains do not fragment anymore. In terms of the presence of small dust, this effect does matter: it means that a driftlimited size distribution cannot sustain efficient fragmentation on long timescales and thus contains much smaller amounts of small dust grains (see right panel of Fig. 1). The observational fact that disks are observed to be rich in small dust for several millions of years thus seems to be in contradiction with a completely driftdominated disk unless there is an external source of small dust (see Dominik & Dullemond 2008) or small dust is mixed and fragmented into surface layers with stronger turbulence, an effect that we cannot take into account in a vertically integrated dust evolution code. We have shown that radial drift alone cannot sustain the relative velocities to keep fragmentation active for the observed lifetimes of protoplanetary disks. Fragmentation together with vertical mixing is therefore likely driven by turbulent motion.
5.2 Analytical dust surface density profiles
The evolution of the dusttogas ratio is important for theories of planetesimal formation (e.g., Johansen et al. 2007, 2011) and planet formation (see Klahr & Bodenheimer 2006; Lissauer & Stevenson 2007; Mordasini et al. 2009). Also a constant dusttogas ratio is usually assumed for deriving disk masses from dust continuum observations (e.g., Andrews & Williams 2005, 2007), which is one of the biggest sources of error in disk mass estimates. In this section, we investigate how the dust surface density and thus the dusttogas ratio changes with time for a typical model of a circumstellar disk, starting with a constant dusttogas ratio of 0.01. The following discussion is based on two important aspects:

As discussed in previous sections, the dust mass is concentrated towards the upper end of the size spectrum, which is also where the drift velocity is the highest. This means that the dust flux is governed by the upper end of the size distribution.
A semisteady dust surface density is then set by the rate at which the outer regions provide dust which is inward drifting with a velocity of . Mass conservation then dictates that the dust accretion rate
(37) 
has to be constant for all . This yields a dust surface density profile
(38) 
For a driftdominated distribution, we can derive the drift velocity of the representative size from Eqns. 17 and 4, which results in a surface density profile given by
(39) 
which for a gas surface density profile with index is proportional to .
For the case of a fragmentationdominated distribution, such a result cannot be obtained as uniquely, because very small particles which are smaller than (i.e., they have a Stokes number of , see Eqns. 30 and 32) are coupled to the gas, while for somewhat larger particles, radial drift starts to play a role. For a dust surface density distribution with an upper size limit larger than (cf. Eq. 32), we can derive a stationary surface density profile as we did before but using Eq. 8 as upper size limit. This yields
(40) 
which for constant values of and is proportional to , as the Minimum Mass Solar Nebula profile (Weidenschilling 1977b; Hayashi 1981).
A dust distribution that consists entirely of small dust grains ( for the largest grains) should be so well coupled to the gas that the dusttogas ratio is constant, i.e.
(41) 
because the dust is just comoving with the gas flow.
We have now derived three different shapes of the dust surface density distribution. In Fig. 9, we compare our simple estimates to the previously discussed simulation results of the full model. Both simulations are for a 0.1 disk with AU around a solar mass star assuming a fragmentation velocity of 10 m s. The turbulence in the simulation of the upper panel is weaker () than in the bottom panel (). It can be seen that the simulation in the upper panel becomes driftdominated in the outer parts (i.e. ). The drift estimate (cf., Eq. 39) agrees well with the simulation result in the regions between about 2 and 200 AU, while the inner parts follow the fragmentation limited estimate, Eq. 40. The estimates do not capture the dropoff in the outermost parts of the disks because they assume a constant accretion rate, while in reality, the surface density and thus the accretion rate must decline towards the outer end of the disk. The stronger turbulence in the bottom panel causes the upper limit of the dust size distribution to decrease. The fragmentationdominated region now reaches out to larger radii at all times while grains in the innermost regions become so small that they are well coupled to the gas and thus can be described by the mixing estimate of Eq. 41. Only the regions beyond about 60 AU stay driftdominated. For a lower fragmentation threshold velocity of 1 m s, the entire disk was found to be fragmentation or mixing dominated.
The results shown in Fig. 9 are for 3 Myrs of dust evolution. At earlier times, for example 1 Myr (not shown), the dusttogas ratio in the outer parts is still higher and thus the drift limited region is smaller, covering only the region beyond 4 and 150 AU for the low and high turbulence case, respectively.
Andrews et al. (2012) have recently found a discrepancy in the sizes and shapes of the dust and gas surface densities in the disk of TW Hya, based on observations of 870 m dust emission and CO 32 emission. In Fig. 10, we plot these observationally derived profiles along with the analytical, drift limited dust surface density (using the observational gas surface density in Eq. 39). The analytical profile has been scaled to fit the value of the dust profile from Andrews et al. (2012) at 10 AU. Such a profile would be expected for a moderate level of turbulence at late stages of disk evolution (the age of TW Hya is estimated to be about 10 Myrs, see Kastner et al. 1997; Webb et al. 1999). Fig. 10 shows that the shape of the analytical solution agrees very well with the observations. One notable difference is the shape of the outer edge. Andrews et al. (2012) note that the remnant mismatch between their bestfit model and the 870 m emission is likely related to the shape of the outer edge, which in their case is a sharp cutoff. In contrast to that, our model assumes a constant dust mass accretion rate through the disk. As stated earlier, the mass accretion rate has to drop off at some radius as the dust reservoir in the outer regions gets depleted. The two dust profile in Fig. 10 thus most probably represent two extreme cases, i.e. a too sharp and a too smooth outer edge.
5.3 Evolution of the dusttogas ratio and dust accretion rate
The results presented in the previous section can directly be translated into a profile of the dusttogas ratio. The absolute value of this profile however depends on the remnant dust accretion rate, which is initially higher due to radial drift, but more quickly decreasing than the gas accretion rate, as shown in Fig. 7. The fact that the dust accretion rate is higher than the gas accretion rate multiplied by the initial dusttogas ratio is important for theories of planetesimal or planetary core formation because this allows a dust trap such as a vortex or zonal flow to achieve local dust overdensities. Depending on the trapping efficiency and the rate at which dust is accumulating in the trap, the dust overdensities can reach critical values to trigger the streaming instability and gravoturbulent formation of planetesimals, as discussed in Klahr & Bodenheimer (2006) and Johansen et al. (2006).
In the low turbulence case (top panels of Fig. 8) it can be seen that the disk regions inward of 4 AU get strongly enriched in dust because grains quickly grow to larger sizes and drift. This leads to a pileup effect similar to Youdin & Chiang (2004). This behavior can be understood by dividing Eq. 40 by the gas surface density, which results in a dusttogas ratio proportional to . Since steady state disks obey , we can also write
(42) 
Thus the steep increase in temperature due to viscous heating causes the increase in the dusttogas ratio through the temperature dependence of the gas disk profile.
In the high turbulent case (bottom panels of Fig. 8), the dust enhancement is much less effective because of the more effective fragmentation of dust grains due to higher turbulent velocities. This causes the drift velocities to be lower than in the lowturbulent case and for a given accretion rate, this yields only a modest enhancement of the dusttogas ratio.
It is important to note that so far we have considered the vertically integrated dusttogas ratio, which is what future observations may be able to constrain. For theoretical works, the midplane dusttogas ratio is usually more important. Due to settling, this will be even larger than the integrated value, if grains are large enough to efficiently settle down to the midplane. For particles with Stokes numbers , Eq. 12 describes the dust scale height as function of Stokes number reasonably well. For the low turbulence case, the particles at the drift limit (Eq. 17) reach Stokes number of up to 0.1 at early times, later dropping 0.01 (see also Fig. 5). According to Eq. 12, this means that the dust scale height and thus the midplane dusttogas ratio is enhanced by a factor of 3 to 10 compared to the vertically integrated values shown in Fig. 8. We find that at early times – up to years – the midplane dusttogas ratio in the whole disk is increased at least by a factor of 2 over the canonical values. After that time, the loss of dust mass counteracts the settling effects. This result obviously depends on the initial condition of the dusttogas ratio, which in our case was taken to be a constant value of 0.01.
6 Summary and conclusions
In this paper, we have derived a simple model that describes the radial evolution of the dust surface density under the combined influence of growth and fragmentation of compact grains as well as radial transport mechanisms. This has been achieved by finding upper limits to the grain size distribution, which are functions of time and the physical conditions at each radius. Good agreement between this very simple numerical model and a fullfledged, computationally intensive dust evolution code was found.
Important parameters are the fragmentation threshold velocity , the level of turbulence , the initial dusttogas ratio, and the temperature and density profile of the gas disk. The effective dust transport velocity can be estimated by representing the dust distribution by only two characteristic grain sizes, the small and the large population. This allows us to derive analytical solutions for the dust surface density profiles in protoplanetary disks.
Our findings are summarized in the following:

The simple twopopulation model presented in this work describes the evolution of the dustsurface density and the evolution of the largest grain sizes well. Good agreement between this simplified model and a fullfledged dust evolution code was found. Additionally, the dusttogas ratio, the dust flux, and the size of the drifting particles can be derived from it.

The upper end of the grain size distribution can be described as limited by turbulent fragmentation (cf., Eq. 8) or by radial drift (cf., Eq. 18). Fragmentation due to relative radial drift alone (Eq. 20) is found to be ineffective in nonlaminar disks (see Section 5.1). This supports the theory that disks are turbulent because radial drift alone cannot cause efficient fragmentation over long enough timescales to agree with observations.

We derived three different analytical profiles of the dustsurface density for different physical conditions: firstly, for a strongly drifting dust distribution (Eq. 39), secondly for a fragmentation limited distribution (40), and thirdly for a distribution where all grains are so small that they are well coupled to the gas (Eq. 41). The free parameter of the profiles is the dust accretion rate provided by the outer regions. The analytical profiles were found to fit to the simulation results of the full dust evolution code of Birnstiel et al. (2010) very well.

We found that at late times of their evolution, disks can become driftdominated, which for typical gas disk profiles (, ) leads to a dust surface density profile proportional to . These results agree with the bestfit models for dust and gas profiles of the 10 Myr old TW Hya disk, as recently found by Andrews et al. (2012).

Similar to Youdin & Chiang (2004), the vertically integrated dusttogas ratio is strongly enhanced in the innermost regions if the dust is significantly drifting in the outer region of the disk. This pileup is aided by the combined effects of drift and fragmentation.

An enhancement of the midplane dusttogas ratio was found in the lowturbulent simulation (). This is because particles can reach large enough Stokes numbers (up to ) to efficiently settle to the midplane. In the case of efficient fragmentation, particles remain too small to significantly enhance the midplane dusttogas ratio.

As an important parameter for models of planetesimal or planetary core formation the radial mass flux in solids has been determined as a time and space dependent function. Depending on the disk parameters, the dust accretion rates can lie anywhere from a factor of a few up to two orders of magnitude above the value expected from the gas accretion rate, scaled by the dusttogas ratio. These values allow a dust trap such as a vortex or zonal flow to achieve a locally critical overdensity on the order of only a few tens of orbits and could trigger the streaming instability and gravoturbulent formation of planetesimals, as discussed in Klahr & Bodenheimer (2006) and Johansen et al. (2006).
Acknowledgements.
We like to thank Satoshi Okuzumi, Chris Ormel, Jürgen Blum, Anna Hughes, Phil Armitage, Sean Andrews, Luca Ricci for interesting discussions. We also thank the anonymous referee for a thorough and insightful report, which helped to improve this paper.Appendix A Simplified dust transport equation
Assuming the gas to be the dynamically dominant medium, the radial transport of each trace species can be described by an advectiondiffusion equation,
(43) 
where is the surface density of the tracer, the gas surface density and and are the tracer velocity and diffusivity, respectively.
In our case, the two trace species considered are the small and the large representative grain sizes, as described in Section 4. The fraction of the total mass for each species is defined in Eq. 22. Thus, the evolution of the total dust surface density can be written as
(44)  
Using Eq. 22, we can simplify this to an advectiondiffusion equation for the total dust surface density
(45) 
where
(46) 
(47) 
and is the mass weighted velocity given by Eq. 24. The dust diffusivities are (see Youdin & Lithwick 2007)
(48) 
For Stokes numbers smaller than one, which is always fulfilled in this paper, the diffusivities are basically equal to the gas diffusivity and thus and . This allows us to simulate the radial evolution of the total dust surface density by using just the gas diffusivity and the mass weighted velocity, as in Eq. 23.
Appendix B Algorithm of the twopopulation model
In this section, we summarize the algorithm of the twopopulation model.

calculate the representative size for a fragmentation limited size distribution given by Eq. 8.

calculate the representative size for a drift limited size distribution given by Eq. 18.

calculate the representative size in the case of drift induced fragmentation
(49) where given by Eq. 20.

the representative size of the large population is now given by the smallest of the size limits
(50) 
the initial phase where particles grow from the smallest size to is approximated by
(51) where
(52) 
the velocities and of the two populations are then given by
(53) where
(54) is the drift velocity, the gas pressure at the midplane
(55) and St the Stokes number is defined by
(56) 
the effective dust transport velocity is the mass averaged velocity of both populations
(57) where the mass fraction depends on the size limiting mechanism
(58) 
Having calculated the effective transport velocity , the evolution of the dust surface density can be simulated by numerically solving the advectiondiffusion equation, Eq. 23.
Footnotes
 More precisely, the solid line is the representative size in a fragmentation limited distribution that is by definition slightly below the upper end of the size distribution, as can be seen in Fig. 1.
 The values in this table represent the values of the model parameters that best reproduce the results of the full simulations.
 We carried out 39 simulations using as initial condition either selfsimilar solutions (Eq. 7) with AU, , disk masses of 0.01 and 0.1 , fragmentation velocities of 1, 3, and 10 m s or powerlaw models with gas surface density exponents of 1 with outer radii of 250 AU.
References
 Aikawa, Y. & Nomura, H. 2006, ApJ, 642, 1152
 Andrews, M. R. & Williams, J. P. 2007, ApJ, 671, 1800
 Andrews, S. M. & Williams, J. P. 2005, ApJ, 631, 1134
 Andrews, S. M., Wilner, D. J., Hughes, A. M., et al. 2012, ApJ, 744, 162
 Bai, X.N. & Goodman, J. 2009, ApJ, 701, 737
 Beitz, E., Güttler, C., Blum, J., et al. 2011, ApJ, 736, 34
 Birnstiel, T., Dullemond, C. P., & Brauer, F. 2009, A&A, 503, L5
 Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, 79
 Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2011, A&A, 525, 11
 Blum, J. & Wurm, G. 2008, ARA&A, 46, 21
 Brauer, F., Dullemond, C. P., & Henning, T. 2008a, A&A, 480, 859
 Brauer, F., Dullemond, C. P., Johansen, A., et al. 2007, A&A, 469, 1169
 Brauer, F., Henning, T., & Dullemond, C. P. 2008b, A&A, 487, L1
 Ciesla, F. J. 2009, Icarus, 200, 655
 Cuzzi, J. N., Hogan, R. C., Paque, J. M., & Dobrovolskis, A. R. 2001, ApJ, 546, 496
 Dominik, C. & Dullemond, C. P. 2008, A&A, 491, 663
 Dullemond, C. P. & Dominik, C. 2004, A&A, 421, 1075
 Dutrey, A., Henning, T., Guilloteau, S., et al. 2007, A&A, 464, 615
 Dzyurkevich, N., Flock, M., Turner, N. J., Klahr, H., & Henning, T. 2010, A&A, 515, 70
 Flock, M., Dzyurkevich, N., Klahr, H., Turner, N. J., & Henning, T. 2011, ApJ, 735, 122
 Garaud, P. 2007, ApJ, 671, 2091
 Geretshauser, R. J., Speith, R., & Kley, W. 2011, A&A, 536, 104
 Gundlach, B., Kilias, S., Beitz, E., & Blum, J. 2011, Icarus, 214, 717
 Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, 56
 Hayashi, C. 1981, Prog. Theor. Phys. Suppl., 70, 35
 Hughes, A. L. H. & Armitage, P. J. 2010, ApJ, 719, 1633
 Johansen, A. & Klahr, H. 2005, ApJ, 634, 1353
 Johansen, A., Klahr, H., & Henning, T. 2006, ApJ, 636, 1121
 Johansen, A., Klahr, H., & Henning, T. 2011, A&A, 529, 62
 Johansen, A., Oishi, J. S., Low, M.M. M., et al. 2007, Nature, 448, 1022
 Kastner, J. H., Zuckerman, B., Weintraub, D. A., & Forveille, T. 1997, Science, 277, 67
 Keller, C. & Gail, H.P. 2004, A&A, 415, 1177
 Klahr, H. & Bodenheimer, P. 2006, ApJ, 639, 432
 Klahr, H. & Johansen, A. 2008, Physica Scripta, 130, 4018
 Kornet, K., Różyczka, M., & Stepinski, T. F. 2004, A&A, 417, 151
 Kornet, K., Stepinski, T. F., & Różyczka, M. 2001, A&A, 378, 180
 Lissauer, J. J. & Stevenson, D. J. 2007, PP V, 591
 LyndenBell, D. & Pringle, J. E. 1974, MNRAS, 168, 603
 Mordasini, C., Alibert, Y., & Benz, W. 2009, A&A, 501, 1139
 Nakagawa, Y., Nakazawa, K., & Hayashi, C. 1981, Icarus, 45, 517
 Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375
 Okuzumi, S. 2009, ApJ, 698, 1122
 Okuzumi, S., Tanaka, H., Takeuchi, T., & Sakagami, M. 2011, ApJ, 731, 96
 Ormel, C. W. & Cuzzi, J. N. 2007, A&A, 466, 413
 Panić, O., Hogerheijde, M. R., Wilner, D., & Qi, C. 2008, A&A, 491, 219
 Rózyczka, M., Kornet, K., Bodenheimer, P., & Stepinski, T. 2004, Rev. Mex. Astron. Astrof. Conf. Ser., 22, 91
 Sano, T., Miyama, S. M., Umebayashi, T., & Nakano, T. 2000, ApJ, 543, 486
 Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337
 Takeuchi, T. & Lin, D. N. C. 2002, ApJ, 581, 1344
 Turner, N. J., Carballido, A., & Sano, T. 2010, ApJ, 708, 188
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2008, ApJ, 677, 1296
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2009, ApJ, 702, 1490
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2011, ApJ, 737, 36
 Webb, R. A., Zuckerman, B., Platais, I., et al. 1999, ApJ, 512, L63
 Weidenschilling, S. J. 1977a, MNRAS, 180, 57
 Weidenschilling, S. J. 1977b, Ap&SS, 51, 153
 Weidenschilling, S. J. 1997, Icarus, 127, 290
 Williams, J. P. & Cieza, L. A. 2011, ARA&A, 49, 67
 Youdin, A. N. & Chiang, E. I. 2004, ApJ, 601, 1109
 Youdin, A. N. & Lithwick, Y. 2007, Icarus, 192, 588
 Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, 57