Gas and dust evolution in protoplanetary disks
T. Birnstiel  C. P. Dullemond  F. Brauer^{}
Junior Research Group at the MaxPlanckInstitut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany
Received 24 November 2009 / Accepted 26 January 2010
Abstract
Context. Current models of the size and radial evolution of
dust in protoplanetary disks generally oversimplify either the radial
evolution of the disk (by focussing at one single radius or by
using steady state disk models) or they assume particle growth to
proceed monodispersely or without fragmentation. Further studies of
protoplanetary disks  such as observations, disk chemistry and
structure calculations or planet population synthesis models 
depend on the distribution of dust as a function of grain size and
radial position in the disk.
Aims. We attempt to improve upon current models to be able to
investigate how the initial conditions, the buildup phase, and the
evolution of the protoplanetary disk influence growth and transport
of dust.
Methods. We introduce a new model similar to Brauer et al.
(2008, A&A, 480, 859) in which we now include the timedependent
viscous evolution of the gas disk, and in which more advanced input
physics and numerical integration methods are implemented.
Results. We show that grain properties, the gas pressure
gradient, and the amount of turbulence are much more influencing the
evolution of dust than the initial conditions or the buildup phase of
the protoplanetary disk. We quantify which conditions or environments
are favorable for growth beyond the meter size barrier. High gas
surface densities or zonal flows may help to overcome the problem of
radial drift, however already a small amount of turbulence poses a much
stronger obstacle for grain growth.
Key words: accretion, accretion disks  circumstellar matter  stars: formation  stars: premain sequence  infrared: stars
1 Introduction
The question of how planets form is one of the key questions in modern astronomy today. While it has been studied for centuries, the problem is still far from being solved. The agglomeration of small dust particles to larger ones is believed to be at least the first stage of planet formation. Both laboratory experiments (Blum et al. 2000) and observations of the 10 m spectral region (Bouwman et al. 2001; van Boekel et al. 2003) conclude that grain growth must take place in circumstellar disks. The growth from submicron size particles to many thousand kilometer sized planets covers 13 orders of magnitude in spatial scale and over 40 orders of magnitude in mass. To assemble a single 1 km diameter planetesimal requires the agglomeration of about 10^{27} dust particles. These dynamic ranges are so phenomenal that one has to resort to special methods to study the growth of particles though aggregation in the context of planet(esimal) formation.
A commonly used method for this purpose makes use of particle size distribution functions. The time dependent evolution of these particle size distribution functions has been studied by Weidenschilling (1980), Nakagawa et al. (1981), Dullemond & Dominik (2005), Brauer et al. (2008a, hereafter (B08a)) and others. It was concluded that dust growth by coagulation can be very quick initially (in the order of thousand years to grow to centimeter sized aggregates at 1 AU in the solar nebula), but it stalls around decimeter to meter size due to what is known as the ``meter size barrier'': a size range within which particles achieve large enough velocities to undergo destructive collisions and fast radial inward drift toward the central star (Weidenschilling 1977; Nakagawa et al. 1986).
While the existence of this meter size barrier (ranging in fact from a couple of centimeters to tens of meters at 1 AU) has been known for over 30 years, a thorough study of this barrier, including all known mechanisms that induce motions of dust grains in protoplanetary disks, and at all regions in the disk, for various conditions in the disk and for different properties of the dust (such as sticking efficiency and critical fragmentation velocity), has been only undertaken recently (see B08a). It was concluded that the barrier is indeed a very strong limiting factor in the formation of planetesimals from dust, and that special particle trapping mechanisms are likely necessary to overcome the barrier.
However, this work was based on a static, nonevolving gas disk model. It is known that over the duration of the planet formation process the disk itself also evolves dramatically (LyndenBell & Pringle 1974; Hueso & Guillot 2005; Hartmann et al. 1998), which may influence the processes of dust coagulation and fragmentation. It is the goal of this paper to introduce a combined diskevolution and dustevolution model which also includes additional physics: we include relative azimuthal velocities, radial dependence of fragmentation critical velocities and the Stokesdrag regime for small Reynolds numbers.
The aim is to find out what the effect of disk formation and evolution is on the process of dust growth, how the initial conditions affect the final outcome, and whether certain observable signatures of the disk (for instance, its degree of dustiness at a given time) can tell us something about the physics of dust growth.
Moreover, this model will serve as a supporting model for complementary modeling efforts such as the modeling of radiative transfer in protoplanetary disks (which needs information about the dust properties for the opacities) and modeling of the chemistry in disks (which needs information about the total amount of dust surface area available for surface chemistry). In this paper we describe our model in quite some detail, and thus provide a basis for future work that will be based on this model.
Furthermore, additional physics, such photoevaporation or layered accretion can be easily included, which we aim to do in the near future.
As outlined above, this model includes already many processes which influence the evolution of the dust and the gas disk. However, there are several aspects we do not include such as backreactions by the dust through opacity or collective effects (Weidenschilling 1997), porosity effects (Ormel et al. 2007), grain charging (Okuzumi 2009) or the ``bouncing barrier'' (Zsom et al. 2010).
This paper is organized as follows: Sect. 2 will describe all components of the model including the radial evolution of gas and dust, as well as the temperature and vertical structure of the disk and the physics of grain growth and fragmentation. In Sect. 3 we will compare the results of our simulations with previous steadystate disk simulations and review the aforementioned growth barriers. As an application, we show how different material properties inside and outside the snow line can cause a strong enhancement of dust within the snow line. Section 4 summarizes our findings. A detailed description of the numerical method along with results for selected test cases can be found in the Appendix.
2 Model
The model presented in this work combines a 1D viscous gas disk evolution code and a dust evolution code, taking effects of radial drift, turbulent mixing, coagulation and fragmentation of the dust into account. We model the evolution of gas and dust in a vertically integrated way. The gas disk is viscously evolving after being built up by infalling material from a collapsing molecular cloud core.
The radial distribution of grains is subject to gas drag, radial drift, and turbulent mixing. To which extend each effect contributes, depends on the grain/gas coupling of each grain size. By simultaneously modeling about 100200 different grain sizes, we are able to follow the detailed evolution of the dust subdisk being the superposition of all sizes of grain distributions.
So far, the evolution of the dust distribution depends on the evolving gas disk but not vice versa. A completely self consistently coupled code is a conceptually and numerically challenging task which will be the subject of future work.
2.1 Evolution of gas surface density
Our description of the viscous evolution of the gas disk follows closely the models described by Nakamoto & Nakagawa (1994) and Hueso & Guillot (2005). In this paper we shall therefore be relatively brief and put emphasis on differences between those models and ours.
The viscous evolution of the gas disk can be described by the continuity equation,
where the gas radial velocity is given by (see LyndenBell & Pringle 1974)
is the gas surface density, r the radius along the disk midplane and the gas viscosity. The source term on the right hand side of Eq. (1), denoted by can be either infall of material onto the disk or outflow.
The molecular viscosity of the gas is too small to account for relevant accretion onto the star, the time scale of viscous evolution would be in the order of several billion years. Observed accretion rates and disk lifetimes can only be explained if turbulent viscosity drives the evolution of circumstellar disks. Therefore Shakura & Sunyaev (1973) parameterized the unknown viscosity as the product of a velocity scale and a length scale. The largest reasonable values for these scales in the disk are the pressure scale height
(3) 
and the sound speed . Therefore the viscosity is written as
(4) 
where is the turbulence parameter and .
Today it is generally believed that disks transport angular momentum by turbulence, however the source of this turbulence is still debated. The magnetorotational instability is the most commonly accepted candidate for source of turbulence (Balbus & Hawley 1991). Values of are expected to be in the range of 10^{3} to some 10^{2} (see Johansen & Klahr 2005; Dzyurkevich et al. 2010). Observations confirm this range with higher probability for larger values (see Andrews et al. 2009).
If the disk becomes gravitationally unstable, large scale mechanisms of angular momentum transport such as through the formation of spiral arms come into play. The stability of the disk can be described in terms of the Toomre parameter (Toomre 1964)
(5) 
Values below a critical value of describe a weakly unstable disk, which forms nonaxisymmetric instabilities like spiral arms. Q values below unity lead to fragmentation of the disk.
The effect of these nonaxisymmetric structures is to transport angular momentum outward and rearranging the surface density in the disk so as to counteract the unstable configuration. This mechanism is therefore to a certain extent selflimiting. Values above are not influenced by instabilities, values below form instabilities which increase Q until the disk is marginally stable again. Our model approximates this mechanism by increasing the turbulence parameter according to the recipe of Armitage et al. (2001),
where is a free parameter of the model which is taken to be 10^{3}, unless otherwise noted.
During the time of infall onto the disk, we use a constant, high value of to mimic the effective redistribution of surface density during the infall phase which also increases the overall stability of the disk. Once the infall stops, we gradually decrease the turbulence parameter on a timescale of 10 000 years until it reaches its input value.
Equation (1) is a diffusion equation, which is stiff. This means, one faces the problem that the numerical step of an explicit integration scheme goes (where is the radial grid step size) which would make the simulation very slow. One possible solution to this problem is using the method of implicit integration. This scheme keeps the small time scales of diffusion i.e. the fast modes in check. We are not interested in these high frequency modes, but they would become unstable if we used a large time step. With an implicit integration scheme (see Sect. A.1) the time step can be chosen larger without causing numerical instabilities, thus increasing the speed of the computation.
2.2 Radial evolution of dust
If the average dusttogas ratio in protoplanetary disks is in the order of 10^{2} (as found in the ISM), then the dust does not dynamically influence the gas, while the gas strongly affects the dynamics of the dust.
Thermally, however, the dust has potentially a massive influence on the gas disk evolution through its opacity, but we will not include this in this paper. Therefore the evolution of the gas disk can, in our approximation, be done without knowledge of the dust evolution, while the dust evolution itself does need knowledge of the gas evolution.
There might be regions, where dust accumulates (such as the midplane of the disk or deadzones and snowlines) and its influence becomes significant or even dominant but we will not include feedback of such dust enhancements onto the disk evolution in this paper.
For now, we want to focus on the equations of motion of dust
particles under the assumption that gas is the dominant material by
mass. The interplay between dust and gas can then be described in terms
of a dimensionless coupling constant, the Stokes number which is defined as
where is the eddy turnover time and is the stopping time.
The stopping time of a particle is defined as the ratio of the momentum
of a particle divided by the drag force acting on it. There are
four different regimes for the drag force which determine the
dusttogas coupling (see Whipple 1972; Weidenschilling 1977). Which regime applies to a certain particle, depends on the ratio between mean free path
of the gas molecules and the dust particle size a (i.e. the Knudsen number) and also on the particle Reynoldsnumber
with
being the gas molecular viscosity
(8) 
the mean thermal velocity. The mean free path is taken to be
(9) 
where n denotes the midplane number density and cm^{2}.
The different regimes^{} are
Here u denotes the velocity of the dust particle with respect to the gas, denotes the mean thermal velocity of the gas molecules, is the solid density of the particles and is the local gas density.
For now, we will focus on the Epstein regime. To calculate the
Stokes number, we need to know the eddy turnover time. As noted
before, our description of viscosity comes from a dimensional analysis.
We use a characteristic length scale
and a characteristic velocity scale
of the eddies. This prescription is ambiguous in a sense that it does
not specify if angular momentum is transported by large, slow moving
eddies or by small, fast moving eddies, that is
(11) 
This is rather irrelevant for the viscous evolution of the gas disk, since all values of q lead to the same viscosity, but if we calculate the turnovereddy time, we get
(12) 
The Stokes number and therefore the dusttogas coupling as well as the relative particle velocities strongly depend on the eddy turnover time and therefore on q. In this work q is taken to be 0.5 (following Schräpler & Henning 2004; Cuzzi et al. 2001) which leads to a turnovereddy time of
(13) 
The Stokes number then becomes
The overall radial movement of dust surface density can now be described by an advectiondiffusion equation,
(15) 
where the total Flux, has contributions from a diffusive and an advective flux. The diffusive part comes from the fact that the gas is turbulent and the dust couples to the gas. The dust is therefore turbulently mixed by the gas. Mixing counteracts gradients in concentration, in this case it is the dusttogas ratio of each size that is being smoothed out by the turbulence. The diffusive flux can therefore be written as
(16) 
The ratio of gas diffusivity over dust diffusivity is called the Schmidt number. We follow Youdin & Lithwick (2007), who derived
(17) 
and assume the gas diffusivity to be equal to the turbulent gas viscosity .
The second contribution to the dust flux is the advective flux,
(18) 
where is the radial velocity of the dust. There are two contributions to it,
The first term is a drag term which comes from the radial movement of the gas which moves with a radial velocity of , given by Eq. (2). Since the dust is coupled to the gas to a certain extend, the radially moving gas is able to partially drag the dust along.
The second term in Eq. (19)
is the radial drift velocity with respect to the gas. The gas in a
Keplerian disk does in fact move subkeplerian, since it feels the
force of its own pressure gradient which is usually pointing inwards.
Larger dust grains do not feel this pressure and move on a keplerian
orbit. Therefore, from a point of view of a (larger) dust particle,
there exists a constant head wind, which causes the particle to loose
angular momentum and to drift inwards. This depends on the coupling
between the gas and the particle and is described by the second term in
Eq. (19). u_{n} denotes the maximum drift velocity of a particle,
which has been derived by Weidenschilling (1977). Here, we introduced a radial drift efficiency parameter . This parameter describes how efficient the radial drift actually is, as several mechanisms such as zonal or meridional flows might slow down radial drift. This will be investigated in Sect. 3.5.
Putting all together, the time dependent equation for the surface density of one dust species of mass m_{i} is given by
where we have included a source term on the right hand side which can be positive in the case of infall or recondensation of grains and negative in the case of dust evaporation or outflows. This source term does not include the sources caused by coagulation and fragmentation processes (see Sect. 2.6.1). All sources will be combined into one equation later which is implicitly integrated in an unsplit scheme (see Sect. A.3).
Note that both, the diffusion coefficient and the radial velocity depend on the Stokes number and therefore on the size of the particle.
Figure 1: Total amount of infallen surface density as function of radius according to the ShuUlrich infall model (see Sect. 2.5) assuming a centrifugal radius of 8 AU (solid line), 33 AU (dashed line), and 100 AU (dashdotted line). 

Open with DEXTER 
2.3 Temperature and vertical structure
The vertical structure can be assumed as being in hydrostatic equilibrium at all times if the disk is geometrically thin (
)
and the vertical sound crossing time is much shorter than the radial
drift time scale of the gas. The isothermal vertical density structure
is then given by
(22) 
where
The viscous heating is given by Nakamoto & Nakagawa (1994)
(24) 
where denotes the turbulent gas viscosity and the Kepler frequency.
Nakamoto & Nakagawa (1994) give a solution for the midplane temperature with an optically thick and an optically thin contribution,
where we used and the approximation . and are Rosseland and Planck mean opacities which will be discussed in the next section.
contains contributions due to stellar or external irradiation. Here, we use a formula derived by Ruden & Pollack (see Ruden & Pollack 1991, Appendix B),
(26) 
with a fixed , following Hueso & Guillot (2005).
The main source of opacity is the dust. Due to viscous heating, the
temperature will increase with surface density. If the temperature
rises above
K,
the dust (i.e. the source of opacity) will evaporate, which stops
the disk from further heating until all dust is vaporized.
To simulate this behavior in our model, we calculate a gas
temperature (assuming a small, constant value for gas opacity) in the
case where the dust temperature rises above the evaporation
temperature. Then
is approximated by
(27) 
only if from Eq. (25) would be larger than .
This is a thermostat effect: if T rises above 1500 K, dust will evaporate, the opacity will drop and the temperature stabilizes at T=1500 K. Once even the very small gas opacity is enough to get temperatures >1500 K, all the dust is evaporated and the temperature rises further.
2.4 Opacity
In the calculation of the midplane temperature we use Rosseland and Planck mean opacities which are being calculated from a given frequency dependent opacity table. The results are stored in a look up table and interpolated during the calculations. The opacity table is for a mixture of 50% silicates and 50% carbonaceous grains.
Since these are dust opacities, we convert them from cm^{2}/g dust to cm^{2}/g gas by multiplying the values with the dusttogas ratio , which is a fixed parameter in our model, taken to be the canonical value of 0.01. This assumes that the mean opacity of the gas is very small and that the dusttogas ratio does not change with time. To calculate the opacity selfconsistently, the total mass of dust and the distribution of grain sizes has to be taken into account, meaning that the dust evolution has a back reaction on the gas by determining the opacity. For now, our model does not take backreactions from dust to gas evolution into account. Only in the case where the temperature rises above 1500 K, the drop of opacity due to dust evaporation is considered, as described above.
2.5 Initial infall phase: cloud collapse
The evolution of the protoplanetary disk also depends on the initial infall phase which builds up the disk from the collapse of a cold molecular cloud core. This process is still not well understood. First similarity solutions for a collapsing sphere were calculated by Larson (1969) and Penston (1969). Shu (1977) found a self similar solution for a singular isothermal sphere. The Larson & Penston solution predicts much larger infall rates compared to the insideout collapse of Shu ( and respectively).
More recent work has shown that the infall rates are not constant over time, but develop a peak of high accretion rates and drop to smaller accretion rates at later times. The maximum accretion rate is about if opacity effects are included (see Larson 2003). Analytical, pressurefree collapse calculations of Myers (2005) show similar behavior but with a smaller peak accretion rate of . They also argue that the effects of pressure and magnetic fields will further increase the time scales of cloud collapse.
This initial infall phase is important since it provides the initial condition of the disk and also influences the whole simulation by providing a source of small grains and gas at larger distances to the star during later times of evolution.
It should be noted that several groups perform 3D hydrodynamic simulations of collapsing cloud cores which show more complicated evolution (e.g., Banerjee & Pudritz 2006; Whitehouse & Bate 2006). However, to be able to study general trends of the infall phase, we use the Shumodel since it is adjustable by a few parameters whose influences are easy to understand. In this model the collapse proceeds with an infall rate of which stays constant throughout the collapse.
We assume the singular isothermal sphere of mass
to be in solid body rotation
.
If infalling material is conserving its angular momentum, all matter from a shell of radius
will fall onto the star and disk system (with mass
)
within the so called centrifugal radius,
where G is the gravitational constant and . The path of every parcel of gas can then be described by a ballistic orbit until it reaches the equatorial plane. The resulting flow onto the disk is
(29) 
where
(30) 
and
(31) 
as described in Ulrich (1976). Here, is given by
(32) 
The centrifugal radius can therefore be approximated by (cf. Hueso & Guillot 2005)
(33) 
We admit that this recipe for the formation of a protoplanetary disk is perhaps oversimplified. Firstly, as shown by Visser et al. (2009), the geometrical thickness of the disk changes the radial distribution of infalling matter onto the disk surface, because the finite thickness may ``capture'' an infalling gas parcel before it could reach the midplane. Secondly, star formation is likely to be messier than our simple singlestar axisymmetric infall model. And even in such a simplified scenario, the Shu insideout collapse model is often criticized as being unrealistic. However, it would be far beyond the scope of this paper to include a better infall model. Here we just want to get a feeling for the effect of initial conditions on the dust growth, and we leave more detailed modeling to future work.
2.6 Grain growth and fragmentation
The goal of the model described in this paper is to trace the evolution of gas and dust during the whole lifetime of a protoplanetary disk, including the grain growth, radial drift and turbulent mixing.
Here, the problem arises that radial drift and coagulation ``counteract'' each other in the regime of particles: coagulation of smaller sizes restores the population around , whereas radial drift preferentially removes these particles. To be able to properly model this behavior, the time step has to be chosen small enough if the method of operator splitting is used.
The upper limit for this time step can be very small. If larger steps are used the solution will ``flipflop'' back and forth between the two splitted substeps of motion and coagulation, and the results become unreliable. A method to allow the choice of large time steps while preserving the accuracy is to use a nonsplitted scheme in which the integration is done ``implicitly''. (B08a) already use this technique to avoid flipflopping between coagulation and fragmentation of grains, and we refer to that paper for a description of the general method. What is new in the current paper is that this implicit integration scheme is extended to also include the radial motion of the particles. So now radial motion, coagulation and fragmentation are done all within a single implicit integration scheme. See Appendix A.3 for details.
2.6.1 Smoluchowski equation
The dust grain distribution n(m,r,z), which is a function of mass m, distance to the star r and height above the midplane z, describes the number of particles per cm^{3} per gram interval in particle mass. This means that the total dust density in g cm^{3} is given by
(34) 
With this definition of n(m,r,z), the coagulation/fragmentation at one point in the disk can be described by a general twobody process,
where M(m,m',m'') is called the kernel. In the case of coagulation and fragmentation, this is given by
For better readability, the dependency of M on radius and height above the midplane was omitted here. The combined coagulation/fragmentation kernel consists of four terms containing the coagulation kernel K, the fragmentation kernel L and the distribution of fragments after a collision S.
Figure 2: Sources of relative particle velocities considered in this model (Brownian motion velocities are not plotted) at a distance of 10 AU from the star. The turbulence parameter in this simulation was 10^{3}. It should be noted that relative azimuthal velocities do not vanish for very large and very small particles. 

Open with DEXTER 
The first two terms in Eq. (36) correspond to gain (masses m' and mm' coagulate) and loss (m coagulates with m') due to grain growth.
The third term describes the fragmentation of masses m and m', governed by the fragmentation kernel L and the fourth term describes the fact that when masses m' and m'' fragment, they distribute some of their mass via fragments to smaller sizes.
The coagulation and fragmentation kernels will be described in Sect. 2.6.3, the distribution of fragments, S, will be described in the next section.
To be able to trace the size and radial evolution of dust in a combined way, we need to express all contributing processes in terms of the same quantity. Hence, we will formulate the coagulation/fragmentation equation in a vertically integrated way. The vertical integration can be done numerically (as in B08a), however coagulation processes are most important near the midplane, which allows to approximate the kernels as being vertically constant (using the values at the midplane). If the vertical dust density distribution of each grain size is taken to be a Gaussian (see Sect. 2.6.4, Eq. (51)), then the vertical integration can be done analytically, as discussed in Appendix A.2. Unlike the steadystate disk models of (B08a) which have fixed surface density and temperature profiles, we need to recompute the coagulation and fragmentation kernels (which are functions of surface density and temperature) at every time step. Therefore this analytical integration also saves significant amounts of computational time.
We therefore define the vertically integrated dust surface density distribution per logarithmic bin of grain radius,
as
where n(r,z,a) and n(r,z,m) are related through . The total dust surface density at any radius is then given by
(38) 
Defining as in Eq. (37) makes it a gridindependent density unlike the mass integrated over each numerical bin. This way, all plots of are meaningful without knowledge of the size grid that was used. Numerically, however we use the discretized values as defined in the Appendix.
In our description of growth and fragmentation of grains, we always assume the dust particles to have a constant volume density meaning that we do not trace the evolution of porosity of the particles as this is currently computationally too expensive with a statistical code as used in this work. This can be achieved with MonteCarlo methods as in Ormel et al. (2007) or Zsom & Dullemond (2008), however these models do not yet include the radial motion of dust and therefore cannot trace the global evolution of the dust disk.
2.6.2 Distribution of fragments
The distribution of fragments after a collision,
S(m,m',m''), is commonly described by a power law,
The value has been investigated both experimentally and theoretically. Typical values have been found in the range between 1 and 2, by both experimental (e.g., Blum & Muench 1993; Davis & Ryan 1990) and theoretical studies (Ormel et al. 2009). Unless otherwise noted, we will follow (B08a) by using the value of .
In the case where masses of the colliding particles differ by orders of
magnitude, a complete fragmentation of both particles is an
unrealistic scenario. More likely, cratering will occur (Sirono 2004),
meaning that the smaller body will excavate a certain amount of mass
from the larger one. The amount of mass removed from the larger one is
parameterized in units of the smaller body ,
(40) 
The mass of the smaller particle plus the mass excavated from the larger one is then distributed to masses smaller than according to the distribution described by Eq. (39). In this work, we follow (B08a) by assuming to be unity, unless otherwise noted.
2.6.3 Coagulation and fragmentation kernels
The coagulation kernel
K(m_{1},m_{2}) can be factorized into three parts,
(41) 
and, similarly, the fragmentation kernel can be written as
(42) 
Here, denotes the relative velocity of the two particles, is the geometrical cross section of the collision and and are the coagulation and fragmentation probabilities which sum up to unity. In general, all these factors can also depend on other material properties such as porosity, however we always assume the dust grains to have a volume density of g cm^{3}.
The fragmentation probability is still not well known and subject of both theoretical (Paszun & Dominik 2009;
Wada et al. 2008) and experimental research (see Güttler et al. 2010; Blum & Wurm 2008). In this work, we adopt the simple recipe
(43) 
with a transition width and the fragmentation speed as free parameter which is assumed to be 1 m s^{1}, following experimental work of Blum & Muench (1993) and theoretical studies of Leinhardt & Stewart (2009).
2.6.4 Relative particle velocities
The different sources of relative velocities considered here are Brownian motion, relative radial and azimuthal velocities, turbulent relative velocities and differential settling. These contributions will be described in the following, an example of the most important contributions is shown in Fig. 2.
Brownian motion, the thermal movement of particles, dependents
on the mass of the particle. Hence, particles of different mass have an
average velocity relative to each other which is given by
(44) 
Particle growth due to Brownian relative motion is most effective for small particles.
Radial drift, as described in Sect. 2.2
also induces relative velocities since particles of different size are
differently coupled to the gas. The relative velocity is then
(45) 
where the radial velocity of the dust, is given by Eq. (19).
Azimuthal relative velocities are induced by gas drag in a
similar way as radial drift. However while only particles
(plus/minus 2 orders of magnitude) around
are significantly drifting, relative azimuthal velocities do not vanish
for encounters between very large and very small particles
(see Fig. 2). Consequently, large particles are constantly suffering high velocity impacts of smaller ones. According to Weidenschilling (1977) and Nakagawa et al. (1986), the relative azimuthal velocities for gasdominated drag are
where u_{n} is defined by Eq. (20).
Turbulent motion as source of relative velocities is discussed in detail in Ormel & Cuzzi (2007). They also derive closed form expressions for the turbulent relative velocities which we use in this work.
Differential settling is the fifth process we consider contributing to relative particle velocities. Dullemond & Dominik (2004)
constructed detailed models of vertical disk structure describing the
depletion of grains in the upper layers of the disk. They show that the
equilibrium settling speed for particles in the Epstein regime is given
by
which can be derived by equating the frictional force
and the vertical component of the gravity force,
.
To limit the settling speed to velocities smaller than half the vertically projected Kepler velocity, we use
for calculating the relative velocities.
Since we do not resolve the detailed vertical distribution of
particles, we take the scale height of each dust size as average height
above the midplane, which gives
(48) 
Table 1: Parameters of the fiducial model.
The dust scale height is calculated by equating the time scale for settling,
(49) 
and the time scale for stirring (Dullemond & Dominik 2004; Schräpler & Henning 2004; Dubrulle et al. 1995),
(50) 
By limiting the vertical settling velocity as in Eq. (47) and by constraining the dust scale height to be at most equal to the gas scale height, one can derive the dust scale height to be
This prescription is only accurate for the dust close to the midplane, however most of the dust (and hence most of the coagulation/fragmentation processes) take place near the midplane, therefore this approximation is accurate enough for our purposes.
3 Results
3.1 Viscous evolution of the gas disk
We will now focus on the evolution of a disk around a T Tauri like star. We assume the rotation rate of the parent cloud core to be s^{1}, which corresponds to 0.06 times the break up rotation rate of the core. The disk is being builtup from inside out due to the ShuUlrich infall model, with the centrifugal radius being 8 AU. The parameters of our fiducial model are summarized in Table 1.
Figure 3: Evolution of disk surface density distribution ( top) and midplane temperature ( bottom) of the fiducial model described in 3.1. 

Open with DEXTER 
Figure 4: Evolution of disk mass and stellar mass (solid and dashed line on left scale respectively) and accretion rate onto the star (dashdotted line on right scale). Adapted from Fig. 5 in Hueso & Guillot (2005). 

Open with DEXTER 
Figure 3 shows how the gas surface density and the midplane temperature of this model evolve as the disk gets built up, viscously spreads and accretes onto the star. It can be seen that viscous heating leads to a strong increase of temperature at small radii. This effect becomes stronger as the disk surface density increases during the infall phase. After the infall has ceased, the surface density and therefore also the amount of viscous heating falls off.
This effect also influences the position of the dust evaporation radius, which is assumed to be at the radius where the dust temperature exceeds 1500 K. This position moves outwards during the infall (because of the stronger viscous heating described above). Once the infall stops, the evaporation radius moves back to smaller radii as the large surface densities are being accreted onto the star.
Figure 4 shows the evolution of accretion rate onto the star, stellar mass and disk mass. The infall phase lasts until about 150 000 years. At this point, the disk looses its source of gas and smallgrained dust and the disk mass drops off quickly until the disk has adjusted itself to the new condition. This also explains the sharp drop of the accretion rate. The slight increase in the accretion rate afterwards comes from the change in after the infall stops (see Sect. 2.1). Hueso & Guillot (2005) find a steeper, powerlaw decline of the accretion rate after the end of the infall phase because their model does not take the effects of gravitational instabilities into account.
3.2 Fiducial model without fragmentation
We will now focus on the dust evolution of the disk. This fiducial simulation includes only grain growth without fragmentation or erosion. All other parameters as well as the evolution of the gas surface density and midplane temperature are the same as discussed in the previous section. The evolution of this model is visualized in Fig. 5.
The contour levels in Fig. 5 show the vertically integrated dust surface density distribution per logarithmic bin of grain radius, , as defined in Eq. (37). The left column of Fig. 5 shows the results of dust evolution for a steady state (i.e. not viscously evolving) gas disk as described in (B08a).
The right column shows the evolution of the dust density distribution of the fiducial model, in which the gas disk is gradually built up through infall from the parent molecular cloud core, and the gas disk viscously spreads and accretes matter onto the star. The solid line marks the grain size corresponding to at the given radius. This grain size will be called hereafter. In the Epstein regime, is proportional to the gas surface density of the disk, which can be seen from Eq. (14).
There are several differences to the simulations of grain growth in steadystate disks, presented in (B08a): viscously evolving disks accrete onto the star by transporting mass inwards and angular momentum outwards. This leads to the fact that some gas has to be moving to larger radii because it is ``absorbing'' the angular momentum of the accreting gas. This can be seen in numerical simulations, but also in the self similar solutions of LyndenBell & Pringle (1974). They show that there is a radius which divides the disk between inward and outward moving gas. This radius itself is constantly moving outwards, depending on the radial profile of the viscosity.
Figure 5: Snapshots of the vertically integrated dust density distributions (defined in Eq. (37)) of a steady state disk as in (B08a) (left column) and of an evolving disk (fiducial model, right column). No coagulation is calculated within the evaporation radius (denoted by the dashdotted line), fragmentation is not taken into account in both simulations. The solid line shows the particle size corresponding to a Stokes number of unity. Since (see Eq. (14)) this curve in fact has the same shape as , so it reflects, as a ``bonus'', what the gas disk looks like. The radius dividing the evolving disk into accreting and expanding regions is marked by the dotted line and the arrows. Particles which are located below the dashed line have a positive flux in the radial direction due to coupling to the expanding gas disk and turbulent mixing (particles within the closed contour in the upper right plot have an inward pointing flux). 

Open with DEXTER 
The radius in the fiducial model was found to move from around 20 AU at the end of the infall to 42 AU at 1 Myr and is denoted by the dotted line in Fig. 5. Important here is that small particles are well enough coupled to the gas to be transported along with the outward moving gas while larger particles decouple from the gas and drift inwards. Those parts of the dust distribution which lie below the dotted line in Fig. 5 have positive flux in the radial direction due to the gascoupling.
Figure 6: Comparison of the radial dependence of the dusttogas ratio for the stationary simulations (thick lines) and the evolving disk simulations (thin lines). The four panels compare the results after 1 Myr of evolution with/without fragmentation (left/right column) and with/without effects of nonaxisymmetric instabilities (top/bottom row). It can be seen that the dusttogas ratio of the evolving disk simulations is almost everywhere lower than the one of the stationary simulations. See Sect. 3.4 for more details. 

Open with DEXTER 
One might expect that the dotted and dashed lines always coincide for small grains as they are well coupled to the gas, however, it can be seen that this is not the case. The reason for this is that turbulent mixing also contributes to the total flux of dust of each grain size. The smallest particles in between the dotted line and the dashed line in the lower two panels of Fig. 5 do have positive velocities, but due to a gradient in concentration of these dust particles, the flux is still negative.
During the early phases of its evolution, a disk which is built up from inside out quickly grows large particles at small radii which are lost to the star by radial drift. During further evolution, growth timescales become larger and larger (since the dusttogas ratio is constantly decreasing) while only the small particles are mixed out to large radii.
The radial dependence of the dusttogas ratio after 1 Myr is shown in Fig. 6. These simulations show that the initial conditions of the stationary disk models (such as shown in (B08a) and in the left column of Fig. 5) are too optimistic since they assume a constant dusttogas ratio at the start of their simulation throughout the disk which is not possible if the disk is being builtup from inside out unless the centrifugal radius is very large (in which case the disk would probably fragment gravitationally) and grain fragmentation is included to prevent the grains from becoming large and start drifting strongly. Removal of larger grains and outward transport of small grains lead to the fact that the dusttogas ratio is reduced by 0.51.5 orders of magnitude compared to a stationary model. This effect is also discussed in Sect. 3.4.
3.3 Fiducial model with fragmentation
The situation changes significantly, if we take grain fragmentation into account. As discussed in Sect. 2.6.2, we consider two different kinds of outcomes for graingrain collisions with relative velocities larger than the fragmentation velocity : cratering (if the masses differ by less than one order of magnitude) and complete fragmentation (otherwise).
Figure 7: Evolution of the dust density distribution of the fiducial model as in Fig. 5 but with fragmentation included. The dashed contour line (in the lower two panels) around the upper end of the size distribution and around small particles at >60 AU marks those parts of the distribution which have a positive (outward pointing) fluxes. 

Open with DEXTER 
For particles larger than about
,
relative velocities are dominated by turbulent motion (and to a
lesser extend by vertical settling). Since the relative velocities
increase with Stokes number (and therefore with grain size), we
can calculate the approximate position of the fragmentation barrier by
equating the assumed fragmentation velocity
with the approximate relative velocities of the particles,
Particles larger than this size are subject to highvelocity collisions which will erode or completely fragment those particles. This is only a rough estimate as the relative velocities also depend on the size of the smaller particle and radial drift also influence the grain distribution. However Eq. (52) reproduces well the upper end of the size distribution in most of our simulations and therefore helps understanding the influence of various parameters on the outcome of these simulations.
The evolution of the grain size distribution including fragmentation is depicted in Fig. 7. The initial condition is quickly forgotten since particles grow on very short timescales to sizes at which they start to fragment. The resulting fragments contribute again to the growth process until a semiequilibrium of growth and fragmentation is reached.
It can be seen that particles stay much smaller than in the model without fragmentation. This means that they are less affected by radial drift on the one hand and better transported along with the expanding gas disk on the other hand. Consequently, considerable amounts of dust can reach radii of the order of 100 AU, as seen in Fig. 8.
The approximate maximum Stokes number, defined in Eq. (52), is inversely proportional to the temperature (since relative velocities are proportional to ), which means that in regions with lower temperature, particles can reach larger Stokes numbers. By equating drift and drag velocities of the particles (cf. Eq. (19)), it can be shown that the radial velocities of particles with Stokes numbers larger than about , are being dominated by radial drift.
Due to the high temperatures below 5 AU (caused by viscous heating), stays below this value which prevents any significant radial drift within this radius, particles inside 5 AU are therefore only transported along with the accreting gas. Particles at larger radii and lower temperatures can drift (although only slightly), which means that there is a continuous transport of dust from the outer regions into the inner regions. Once these particles arrive in the hot region, they get ``stuck'' because their Stokes number drops below . The gas within about 5 AU is therefore enriched in dust, as seen in Fig. 8. The dusttogas ratio at 1 AU after 1 Myr is increased by 25% over the value of infalling matter, which is taken to be 0.01.
Figure 8 also shows a relatively sharp decrease in the dust to gas ratio at a few hundred AU. At these radii, the gas densities become so small that even the smallest dust particles decouple from the gas and start to drift inwards.
The thick line in Fig. 8 shows as comparison the dusttogas ratio of the stationary disk model (cf. left column of Fig. 5) after 1 Myr, which starts with a radially constant initial dusttogas ratio of 0.01.
Figure 9 shows the semiequilibrium grain surface density distribution at 1, 10 and 100 AU in the fiducial model with fragmentation at 1 Myr. The exact shape of these distributions depends very much on the prescription of fragmentation and cratering. In general the overall shape of these semiequilibrium distributions is always the same: a powerlaw or nearly constant distribution (in ) for small grains and a peak at some grain size , beyond which the distribution drops dramatically. The peak near the upper end of the distribution is caused by cratering. This can be understood by looking at the collision velocities: the relative velocity of two particles increases with the grain size but it is lower for equalsized collisions than for collisions with particles of very different sizes (see Fig. 2). The largest particles in the distribution have relative velocities with similar sized particles which lie just below the fragmentation velocity (otherwise they would fragment). This means that the relative velocities with much smaller particles (which are too small to fragment the bigger particles but can still damage them via cratering) are above this critical velocity. This inhibits the further growth of the big particles beyond , causing a ``traffic jam'' close to the fragmentation barrier. The peak in the distribution represents this traffic jam.
Figure 8: Evolution of the radial dependence of the dusttogas ratio in the fiducial model including fragmentation with the times corresponding to the snapshots shown in Fig. 7. The initial dusttogas ratio is taken to be 0.01. The thick dashed curve represents the result at 1 Myr of the static disk model for comparison. 

Open with DEXTER 
Figure 9: Vertically integrated (cf. Eq. (37)) grain surface density distributions as function of grain radius at a distance of 1 AU (solid), 10 AU (dashed) and 100 AU (dotdashed) from the star. These curves represent slices through the bottom panel of Fig. 7. 

Open with DEXTER 
3.4 Influences of the infall model
In the fiducial model without fragmentation, continuous resupply of material by infall is the cause why the disk has much more small grains than compared to the stationary disk model (cf. Fig. 5), which relatively quickly consumes all available micrometer sized dust. The effect has already been found in Dominik & Dullemond (2008): if all grains start to grow at the same time, then the bulk of the mass grows in a relatively thin peak to larger sizes (see Fig. 6 in (B08a)). However if the bulk of the mass already resides in particles of larger size, then additional supply of small grains by infall is not swept up effectively because of the following: firstly, the number density of large particles is small (they may be dominating the mass, but not necessarily the number density distribution) and secondly, they only reside in a thin midplane layer while the scale height of small particles equals the gas scale height.
We studied how much the disk evolution depends on changes in the infall model.
For a given cloud mass, the so called centrifugal radius , which was defined in Eq. (28), depends on the temperature and the angular velocity of the cloud. Both can be varied resulting in a large range of possible centrifugal radii reaching from a few to several hundred AU. Since the centrifugal radius is the relevant parameter, we varied only the rotation rate of the cloud core. We performed simulations with three different rotation rates which correspond to centrifugal radii of about 8 (fiducial model), 33 AU, and 100 AU. For each centrifugal radius, we performed two simulations: one which includes effects of gravitational instabilities (GI)  i.e. increased during infall and according to Eq. (6)  and one which neglects them.
However for a centrifugal radius of 100 AU, too much matter is loaded onto the cold outer parts of the disk and consequently, the disk would fragment through gravitational instability. We cannot treat this in our simulations, hence, for the case of 100 AU, we show only results which neglect all GI effects.
The resulting dusttogas ratios are being shown in Fig. 6.
Two general aspects change in the case of higher rotation rates: firstly, more of the initial cloud mass has to be accreted onto the star by going through the outer parts of the disk. Consequently, the disk is more extended and more massive than compared to the case of low rotation rate.
Secondly, as aforementioned the high surface densities in the colder regions at larger radii cause the disk to become less gravitationally stable.
If grain fragmentation is not taken into account in the simulations, both effects cause more dust mass to be transported to larger radii. Growth and drift timescales are increasing with radius and the dust disk with centrifugal radius of 33 AU (100 AU) can stay 5 (35) times more massive than in the low angular momentum case after 1 Myr if GI effects are neglected.
If GI effects are included, matter is even more effectively transported outward, the dusttodisk mass ratio for 8 and 33 AU is increased from 5 to 8.
However if fragmentation is included, it does not matter so significantly, where the dust mass is deposited onto the disk since grains stay so small during the buildup phase of the disk (due to grain fragmentation by turbulent velocities) that they are well coupled to the gas. Outwards of AU (without GI effects) or of a few hundred AU (if GI effects are included), the gas densities become so small that even the smallest grains start do decouple from the gas. They are therefore not as effectively transported outwards. In these regions, the amount of dust depends on the final centrifugal radius while at smaller radii, turbulent mixing quickly evens out all differences in the dusttogas ratio (see left column of Fig. 6).
It can be seen, that in all simulations, the dusttogas ratio is lower than in the stationary disk model. The trend in the upper right panel in Fig. 6 suggests that for a centrifugal radius of 100 AU and the enhanced radial transport by GI effects might have a higher dusttogas ratio than the stationary disk model. However in this case, the disk would become extremely unstable and would therefore fragment.
The reason for this is the following: to be able to compare both simulations, the total mass of the diskstar system is the same as in the stationary disk models. How the total mass is distributed onto disk and star depends on the prescription of infall. If a centrifugal radius of 100 AU is used, the disk becomes so massive that it significantly exceeds the stability criterion .
3.5 The radial drift barrier revisited
According to the current understanding of planet formation, several mechanisms seem to prevent the formation of large bodies via coagulation quite rigorously. The most famous ones  radial drift and fragmentation  have already been introduced above. Radial drift has first been discussed by Weidenschilling (1977), while the importance of the fragmentation barrier (which may prevent grain growth at even smaller sizes) was discussed in (B08a). In the following, we want to test some ideas about how to weaken or to overcome these barriers apart from those already studied in (B08a).
(B08a) has quantified the radial drift barrier by equating the
timescales of growth and radial drift which leads to the condition
where is the dusttogas ratio and and are the drift and growth timescales respectively. The parameter describes how much more efficient growth around must be, so that the particles are not removed by radial drift. To overcome the drift barrier, obviously either particle growth must be accelerated, or the drift efficiency has to be decreased. (B08a) have numerically measured the parameter to be around 12. In other words, this means that the growth timescales have to be decreased (e.g. by an increased dusttogas ratio) until the condition in Eq. (53) is fulfilled.
However, there are other ways of breaking through the drift barrier. Firstly, the drift timescale for
particles also depends on the temperature (via the pressure gradient). A simple approximation from Eq. (53) with a 0.5
star and a dusttogas ratio of 0.01 gives
(54) 
which means that particles should be able to break through the drift barrier at 1 AU if the temperature is below 100 K (or 200 K for a solar mass star). Dullemond et al. (2002) have constructed vertical structure models of passively irradiated circumstellar disks using full frequency and angledependent radiative transfer. They show that the midplane temperature of such a T Tauri like system at 1 AU can be as low as 60 K. Reducing the temperature by some factor reduces the drift time scale by the a factor of similar size which we will call the radial drift efficiency (cf. Eq. (20)).
Zonal flows as presented in Johansen et al. (2006) could be an alternative way of decreasing the efficiency of radial drift averaged over typical time scales of particle growth. Johansen et al. (2006) found a reduction of the radial drift velocity of up to 40% for metersized particles.
Meridional flows (e.g., Kley & Lin 1992; Urpin 1984) might also seem interesting in this context, however they do not directly influence the radial drift efficiency but rather reverse the gasdrag effect. This might be important for small particles (which, however are not strongly settling to the midplane) but for particles, would have to be extremely high to have significant influence: even would result in a reduction of the particles radial velocity by approximately only a few percent.
Another possibility to weaken the drift barrier is changing its radial dependence. The reason for this is the following: particle radial drift is only a barrier if it prevents particles to cross the size . Since particles grow while drifting, the particle size corresponding to needs to increase as well, to be a barrier. Otherwise drifting particles would grow (at least partly) through the barrier while they are drifting. If is decreasing in the direction toward the star, then a particle that drifts inwards would have an increasing Stokes number even if the particle does not grow at all.
Figure 10: Evolution of the dust surface density distribution of the fiducial model at 200 000 years for different drift efficiencies , without fragmentation. The solid line denotes the grain size of particles with a Stokes number of unity. Gas outside of the radius denoted by the dotted line as well as particles below the dashed line have positive radial velocities. See Sect. 3.2 for more discussion. 

Open with DEXTER 
Figure 11: Dust grain surface density distribution as in Fig. 10 at 200 000 years but including the Stokes drag regime. The drift efficiency is set to and fragmentation is not taken into account. It can be seen that (solid line) increases with radius until about 1 AU, which facilitates the break through the drift barrier. 

Open with DEXTER 
In the Epstein regime, the size corresponding to
is proportional to the gas surface density
(55) 
meaning that a relatively flat profile of surface density (or even a profile with positive slope) is needed to allow particles to grow through the barrier. However, our simulations of the viscous gas disk evolution do not yield surface density profiles with positive slopes outside the dust evaporation radius.
To quantify the arguments above, we have performed simulations with varying drift efficiency to test how much the radial drift has to be reduced to allow break through. We have additionally included the first Stokes drag regime to see how the radial drift of particles is influenced by it.
Figure 10 shows the grain surface density distribution after 200 000 years of evolution for three different drag efficiencies. The most obvious changes can be seen in the region where the line (solid line) is relatively flat: the grain distribution is shifted towards larger Stokes numbers. As explained above, particles can grow while drifting, which can be seen in the case of . The Stokes number and size of the largest particles is significantly increasing towards smaller radii. However the radial drift efficiency has to be reduced by 80% to produce particles which are large enough to escape the drift regime and are therefore not lost to the star.
The situation changes, if the Stokes drag is taken into
account: if gas surface densities are high enough for the dust
particles to get into a different drag regime, a change in the
radial dependency of
can be achieved. The Epstein drag regime for particles with Stokes number of unity is only valid if
otherwise, drag forces have to be calculated according to the Stokes drag law since the Knudsen number becomes smaller than 4/9 (see Weidenschilling 1977). The Stokes number is then given by
(57) 
with being the cross section of molecular hydrogen. Interestingly, is independent of and proportional to the square root of the pressure scale height which decreases towards smaller radii. This means that  as long as the surface density is high enough  it does not depend on the radial profile of the surface density. In this regime the radial drift itself could move particles over the drift barrier since drifting inwards increases the Stokes number of a particle without increasing its size.
Results of simulations which include the Stokes drag are shown in Fig. 11. In this case, particles can already break through the drift barrier if . This value and the position of the breakthrough depends on where the drag law changes from Epstein to Stokes regime and therefore on the disk surface density. As noted above, larger surface densities generally shift the position of regime change towards larger radii making it easier for particles to break through the drift barrier.
It should be noted that the physical way to avoid the Stokes drag regime is to decrease the surface densities, however we chose the same initial conditions for both cases  with and without Stokes drag  and just neglected the Stokes drag in the latter computations to be able to compare the efficiency factors independent of other parameters such as disk mass or temperature.
3.6 The fragmentation barrier revisited
In the previous section, we have shown that several mechanisms could allow particles to break through the radial drift barrier, however fragmentation puts even stronger constraints on the formation of planetesimals.
Figure 12: Dust surface density distributions (top row) and the according solidtogas ratio (bottom row) for the case of radiusdependent fragmentation velocity after years of evolution. In the upper row, the vertical dashed line denotes the position of the snow line at the midplane, the solid line corresponds to and the dotdashed line shows the approximate location of the fragmentation barrier according to Eq. (52). The snow line on the right plot lies further outside since viscous heating is stronger if is larger. In the bottom row, the solid line denotes the vertically integrated dusttogas ratio while the dashed line denotes the dusttogas ratio at the disk midplane. The icy dust grains outside the snow line are assumed to fragment at a critical velocity of 10 m s^{1} while particles inside the snow line fragment at 1 m s^{1}. The plots on the left and right hand side differ in the amount of turbulence in the disk ( and 10^{2}, respectively). 

Open with DEXTER 
As shown by Ormel & Cuzzi (2007), the largest relative velocities are of the order of
(58) 
If particles should be able to break through the fragmentation barrier, then they need to survive these large relative velocities, meaning that has to be smaller than the fragmentation velocity of the particles, or
(59) 
This condition is hard to fulfill with reasonable fragmentation velocities, unless is very small. E.g., for and a temperature of 200 to 250 K, the fragmentation velocity needs to be higher than 4 m s^{1}, which could already be seen in the simulations by Brauer et al. (2008b), who have simulated particle growth near the snowline.
Even in the case of very low turbulence, relative azimuthal velocities of large ( ) and small grains ( ) are of the order of 30 m s^{1}, which means that large particles are constantly being ``sandblasted'' by small grains. The only way of reducing these velocities significantly is decreasing the pressure gradient (see Eqs. (20) and (46)).
Another possibility to overcome this problem would be if highvelocity impacts of smaller particles would cause net growth, as has been found experimentally by Wurm et al. (2005) and Teiser & Wurm (2009).
Taken together, these facts make environments as the inner edge of dead zones ideal places for grain growth (see Kretke & Lin 2007; Brauer et al. 2008b): shutting of MRI leads to low values of , which are needed to reduce turbulent relative velocities and the low pressure gradients prevent radial drift and high azimuthal relative velocities.
3.7 Dust enhancement inside the snow line
As already noted in a previous paper (Birnstiel et al. 2009), significant loss of dust by radial drift can be prevented by assuring that particles stay small enough and are therefore not influenced by radial drift. For typical values of ( 10^{3}10^{2}), this means that the fragmentation velocity must be smaller than about 0.55 m s^{1}. If particles have higher tensile strength, they can grow to larger sizes which are again affected by radial drift.
Typical fragmentation velocities for silicate grains determined both theoretically and experimentally are of the order of a few m s^{1} (for a review, see Blum & Wurm 2008). The composition of particles outside the snowline is expected to change due to the presence of ices. This can influence material properties and increase the fragmentation velocity (see Schäfer et al. 2007; Wada et al. 2009).
We have performed simulations with a radially varying fragmentation speed. We assume the fragmentation speed inside the snowline to be 1 m s^{1}, outside the snowline to be 10 m s^{1}. It should be noted that we do not follow the abundance of water in the disk or the composition of grains, we only assume particles outside the snow line to have larger tensile strength due to the presence of ice. To be able to compare both simulations, we used the same 1 Myr old 0.09 gas disk around a solar mass star as initial condition. The gas surface density profile of this disk was derived by a separate run of the disk evolution code. We used this gas surface density profile and a radially constant dusttogas ratio as initial condition for the simulations presented in this subsection. Apart from the level of turbulence, the input for both simulations is identical, the results are therefore completely independent of uncertainties caused by the choice of the infall model.
Results of the simulations are shown in Fig. 12. A one order of magnitude higher fragmentation velocity causes the maximum grain size to be about two orders of magnitude larger, which follows from Eq. (52) since in the Epstein regime (all particles in these simulations are small enough to be in the Epstein regime). This effect can be seen in Fig. 12. Particles outside the snowline are therefore more strongly drifting inwards (because they reach larger Stokes numbers) where they are being pulverized as soon as they enter the region within the snowline.
Strong drift outside the snow line and weaker radial drift inside the snow line cause the dusttogas ratio within the snow line to increase significantly (see bottom row of Fig. 12): in the case of , the dusttogas ratio reaches values between 0.39 and 0.10 in the region from 0.2 to 1.9 AU.
Simulations for a less massive star (0.5 ) show the same behavior, however the increase in dusttogas ratio is not as high as for a solar mass star (dusttogas ratio of 0.270.20 from 0.24 AU).
This effect is not as significant in the case of stronger turbulence, where the maximum dusttogas ratio is around 0.027. The evolution of the dusttogas ratio at a distance of 1 AU from the star is plotted for both cases in Fig. 13 (the minor bump is an artifact of the initial condition).
The reason for this difference lies in the locations of the drift and fragmentation barriers. The approximate position of the fragmentation barrier is represented by the dotdashed line in Fig. 12. The radial drift barrier cannot be defined as sharply, however radial drift is strongest at a Stokes number of unity, which corresponds to the solid line in Fig. 12. An increase of by one order of magnitude lowers the fragmentation barrier by about one order of magnitude in grain size.
In the lower turbulence case, the fragmentation barrier lies close to . Most particles are therefore drifting inwards before they are large enough to experience the fragmentation barrier. Hence, the outer parts of the disk are significantly depleted in small grains.
In contrast to this case, fragmentation is the stronger barrier for grain growth throughout the disk in the high turbulence simulation (right column in Fig. 12). It can be seen that particles of smaller sizes are constantly being replenished by fragmentation.
With these results in mind, the evolution of the disk mass (bottom panel of Fig. 14) seems counterintuitive: the mass of the high turbulence dust disk is decaying faster than in the low turbulence case. This can be understood by looking at the global dusttogas ratio of the disks (top panel of Fig. 14) which does not differ much in both cases. This means that the increased dust mass loss in the high turbulence disk is due to the underlying evolution of the gas disk. Particles in the high turbulence disk may have smaller Stokes numbers (causing drift to be less efficient), however the inward dragging by the accreting gas is stronger in this case.
To show how much the dust evolution depends on the fragmentation velocity, we included the case of a lower fragmentation velocity throughout the disk in Fig. 14. It can be seen that the dust mass is retained at its initial value for much longer timescales.
Figure 13: Dusttogas ratio at a distance of 1 AU from the central star as a function of time for the case of low ( , solid line) and high ( , dashed line) turbulence. 

Open with DEXTER 
Figure 14: Evolution of the global dusttogas ratio ( top panel) and the dust disk mass ( bottom panel) for the simulations shown in Fig. 12. The solid and dashed lines correspond to the low and high turbulence case, respectively. The dashdotted line shows the evolution of the disk if a low fragmentation velocity is assumed throughout the disk. 

Open with DEXTER 
4 Discussion and conclusions
We constructed a new model for growth and fragmentation of dust in circumstellar disks. We combined the (size and radial) evolution of dust of (B08a) with a viscous gas evolution code which takes into account the spreading and accretion, irradiation and viscous heating of the gas disk. The dust model includes the growth/fragmentation, radial drift/drag and radial mixing of the dust. We reimplemented and substantially improved the numerical treatment of the Smoluchowski equation of (B08a) to solve for the combined size and radial evolution of dust in a fully implicit, unsplit scheme (see Appendix A). In addition to that, we also included more physics such as relative azimuthal velocities, radial dependence of fragmentation critical velocities and the Stokes drag regime. The code has been tested extensively and was found to be very accurate and massconserving (see Appendix B).
We compared our results of grain growth in evolving protoplanetary disks to those of steady state disk simulations, similar to (B08a). In spite of many differences in details, we confirm the most general result of (B08a): radial drift and particle fragmentation set strong barriers to particle growth. If fragmentation is included in the calculations, then it poses the strongest obstacle for the formation of planetesimals. Very low turbulence ( ) and fragmentation velocities of more than a few m s^{1} are needed to be able to overcome the fragmentation barrier in the case of turbulent relative velocities.
This model includes also the initial buildup phase of the disk, which is still a very poorly understood phase of disk evolution. We use the ShuUlrich infall model which represents a strong simplification. However, the following novel findings of this work do not or only weakly depend on the buildup phase of the disk:
 Apart from an increased dusttogas ratio (as in B08a), other mechanisms such as streaming instabilities or a decreased temperature may be able to weaken this barrier by decreasing the efficiency of radial drift. We found that in simulations without fragmentation the radial drift efficiency needs to be reduced by 80% to produce particles which crossed the metersize barrier and are large enough to resist radial drift.
 If the gas surface density is above a certain threshold (in our simulations about 140 g cm^{2} at 1 AU or 330 g cm^{2} at 5 AU, see Eq. (56)) then the drag force which acts on the dust particles has to be calculated according to the Stokes drag law, instead of the Epstein drag law. The drift barrier in this drag regime is shifting to smaller sizes for smaller radii (independent on the radial profile of the surface density) which means that pure radial drift can already transport dust grains over the drift barrier or at least to larger Stokes numbers even without simultaneous grain growth. In this case, the drift efficiency needs to be reduced only by about 25% to produce large bodies.
 If relative azimuthal velocities are included, then grains with are constantly ``sandblasted'' by small grains (if they are present) which (in our prescription of fragmentation) causes erosion and stops graingrowth even in the case of low turbulence. Only decreasing the radial pressure gradient significantly weakens both relative azimuthal and radial velocities. Low turbulence and a small radial pressure gradient together are needed to allow larger bodies to form. These conditions may be fulfilled at the inner edge of dead zones (Brauer et al. 2008b; Kretke & Lin 2007; see also Dzyurkevich et al. 2010). Future work needs to investigate the disk evolution and grain growth of disks with dead zones. Our prescription of fragmentation and erosion may also need rethinking since Wurm et al. (2005) and Teiser & Wurm (2009) find netgrowth by high velocity impacts of small particles onto larger bodies.
 Higher tensile strengths of particles outside the snowline allows particles to grow to larger sizes, which are more strongly affected by radial drift. Particles therefore drift from outside the snowline to smaller radii where they fragment and almost stop drifting (since their radial velocity is decreased by almost two orders of magnitude). This can cause an increase the dusttogas ratio inside the snowline by more than 1.5 orders of magnitude.
 The critical fragmentation velocity and its radial dependence (and to a lesser extent the level of turbulence) is a very important parameter determining if the dust disk is drift or fragmentation dominated. A drift dominated disk is significantly depleted in small grains and only a fragmentation dominated disk can retain a significant amount of dust for millions of years as is observed in T Tauri disks.
 Disk spreading causes small particles ( m) to be transported outward at radii of 60190 AU even in 1 Myr old disks.
 Small particles provided by infalling material are not effectively swept up by large grains if the bulk of the dust mass has already grown to larger sizes.
 In an insideout buildup of circumstellar disks, grains growth is very fast (timescales of some 100 years) because densities are high and orbital timescales are small. Large grains are quickly lost due to drift towards the star if fragmentation is neglected. Fragmentation is firstly needed to keep grains small enough to be able to transport a significant amount of dust to large radii by disk spreading and secondly to retain it in the disk by preventing strong radial inward drift.
We like to thank Thomas Henning, Hubert Klahr, Chris Ormel, and Andras Zsom for insightful discussions. We also thank the referee, Hidekazu Tanaka for his fast and insightful review which helped to improve this paper.
Appendix A: Algorithm
In the next two sections, we will first discuss how the equations of radial evolution of gas (Eq. (1)) and dust (Eq. (21)) as well as the coagulation/fragmentation of dust (Eq. (35)) are solved separately. In Sect. A.3 we will then describe how this model treats the radial and the size evolution of dust in an unsplit, fully implicit way.
A.1 Advectiondiffusion algorithm
To be able to also model both, the evolution of dust and gas
implicitly, we constructed a scheme which solves a general form of an
advectiondiffusion equation,
which can be adapted to both, Eqs. (1) and (21) by proper choice of parameters u, , g, h, K and L.
We use a fluxconserving donorcell scheme which is implicit in
.
The time derivative in Eq. (A.1), written in a discretized way becomes
(A.2) 
where i denotes timedimension and n denotes spacedimension.
The advective part is discretized as
where and are the future flux and the surface between cell n and cell n+1 and V_{n} is the volume of cell n. The advective interface fluxes can then be written as
=  (A.4)  
=  (A.5) 
The diffusive interface flux becomes
=  (A.6)  
=  (A.7) 
Working out these equations and separating the values of leads to
with the coefficients
Equation (A.8) can now be solved by any matrixsolver, but since it is a tridiagonal matrix, the fastest analytical way to solve it is by forward elimination/backward substitution.
It should be noted that Eq. (A.1) is implicit only in which means that the equations we solve are only implicit in the surface density. In the case of the viscous accretion disk, described by Eq. (1), we face the problem that also the turbulent gas viscosity depends on the temperature which (in the case of viscous heating) depends on the surface density. This can cause numerical instabilities to develop.
To stabilize the code, we use a scheme which estimates the temperature in several predictor steps. The actual time step is then done with the predicted temperature.
A.2 Coagulation algorithm
Discretizing Eq. (35) on a mass grid m_{i} gives
(A.10) 
where the dust particle number density is
(A.11) 
If we assume that the coagulation and fragmentation kernels are constant in z and that the vertical distribution of grains is a Gaussian with a scale height according to Eq. (51),
(A.12) 
we can vertically integrate the coagulation/fragmentation equation.
where
(A.14) 
Discretizing Eq. (A.13) also in radial direction and rewriting it in terms of the quantity
(A.15) 
yields
where the vector is the source function for each of the m mass bins.
The numerical change of within a time step is . The timediscretized version of Eq. (A.16) then becomes (omitting second order terms)
(A.17) 
Since the first term on the right hand side is the explicit source function, we can write
(A.18) 
Using the vectors and , this can be rewritten in matrix notation,
Where
(A.20) 
denotes the Jacobian of the source function and the unity matrix. The solution for the future values can now be derived by inverting the matrix in Eq. (A.19),
(A.21) 
A.3 Fully implicit scheme for radial motion and coagulation
To be able to solve the radial motion and the Smoluchowski equation fully implicitly, we rewrite Eq. (A.8) as
where M is the tridiagonal matrix with entries A,B and C and source term which are given by Eq. (A.9). M is now rewritten by separating off the diagonal terms and by dividing by ,
(A.23) 
which brings Eq. (A.22) in a form similar to Eq. (A.19),
(A.24) 
The coagulation/fragmentation is now determined by the matrix J and the source vector and similarly, the radial motion is determined by matrix and source vector
(A.25) 
This allows us to combine both operators into one matrix equation,
In principle, to solve for the vector , the matrix on the left hand side of Eq. (A.26) has to be inverted. This is a numerically challenging task since the inverse matrix in our simulations would have about 150500 million elements. The matrix on the left hand side of Eq. (A.26), however is a sparse matrix (schematically depicted in Fig. A.1).
We can therefore solve Eq. (A.26) by an iterative incomplete LU decomposition solver for sparse matrices provided by the Sixpack library^{} of S. E. Norris.
Figure A.1: Pictographic representation of the matrix on the left hand side of Eq. (A.26) with six radial and five mass grid points. White elements are zero, dark grey elements contain contributions from radial transport of dust ( J), light gray elements contain contributions from the coagulation/fragmentation ( ), and black elements contain both contributions. The upper left and lower 5 rows represent the boundary conditions where coagulation/fragmentation is not taken into account. The matrix in the simulations would typically have a size of 15 000^{2}. 

Open with DEXTER 
Figure A.2: Comparison between simulation result and analytical solution for growth of equal sized particles. The solid lines denote the position of the peak of the grain size distribution. In the top panel, only Brownian motion is considered as source of relative particle velocities, in the bottom panel, turbulent relative velocities are considered as well. The parameters of these simulation are T=196 K, = 1.6 , , and . 

Open with DEXTER 
Appendix B: Test cases
To check if the numerical implementation presented above accurately solves Eq. (A.26),
we compare results of the simulation to some analytical solutions: the
growth rates of particles can be approximated if we assume the grain
size distribution to be a delta function and the sticking probability
to be unity. The increase of mass per collision is then given by the
mass of the particle, m divided by the time between two collisions, ,
which can be written as
(B.1) 
Using and , we derive
where is the relative velocity, the dust density, the solid density of the dust particles and the cross section of the collision.
With this formula, we can estimate the growth time scales if the
relative velocities of the dust grains is given. For equal sized
particles and Brownian motion, we get
(B.3) 
for turbulent motion, the relative velocities are (see Ormel & Cuzzi 2007)
(B.4) 
Integrating Eq. (B.2) gives
for Brownian motion and
(B.6) 
for relative velocities due to turbulent motion.
A comparison between analytical solution and simulation result for Brownian motion growth is shown in the top panel of Fig. A.2. It can be seen that the position of the peak of the grain size distribution (solid line in Fig. A.2) follows the analytical result of Eq. (B.5).
Figure B.1: Test case for only radial drift, without coagulation. The solid line denotes the massaveraged position of the radial distribution of grains for each grain size after 10^{3} years. The dashed line is the expected solution for a single particle, the dotted line denotes the initial radius of the particles. 

Open with DEXTER 
A similar comparison for the case of relative velocities due to turbulent motion is not as straightforward since both the turbulent relative velocities and the dust scale height (Eq. (51)) are subdivided into several regimes. We therefore integrated Eq. (B.2) numerically for the case of Brownian motion and turbulent relative velocities; the results are shown in the bottom panel of Fig. A.2. As before, we see that  after the initial condition is overcome  the simulation result follows closely the monodisperse approximation. For grains larger than , the analytical solution is also plotted in the bottom panel of Fig. A.2, almost coinciding with the simulation results.
The radial motion of dust particle was tested in a similar fashion: starting from a grain distribution at a radius of 5.5 AU, we let the particles drift (taking the gas drag and the radial drift into account, see Eq. (19)) without coagulation. We compare the results to results of a numerical integration of the equation of motion for a single particle. The results are shown in Fig. B.1.
We find that the size distribution behaves as expected: small particles are well coupled to the gas, they almost retain their initial position since the radial motion due to gas drag is in the order of 0.01 AU. Larger particles, having a larger Stokes number drift towards the star on shorter timescales. Particles of a few centimeters (corresponding to ) are already lost after about 700 years.
The mass in all test cases was found to be conserved on the order of 10^{11}% of the initial value.
References
 Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2009, ApJ, 700, 1502 [NASA ADS] [CrossRef] [Google Scholar]
 Armitage, P. J., Livio, M., & Pringle, J. E. 2001, MNRAS, 324, 705 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Banerjee, R., & Pudritz, R. E. 2006, ApJ, 641, 949 [NASA ADS] [CrossRef] [Google Scholar]
 Birnstiel, T., Dullemond, C. P., & Brauer, F. 2009, A&A, 503, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blum, J., & Muench, M. 1993, Icarus, 106, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blum, J., Wurm, G., Kempf, S., et al. 2000, Phys. Rev. Lett., 85, 2426 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Bouwman, J., Meeus, G., de Koter, A., et al. 2001, A&A, 375, 950 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brauer, F., Dullemond, C. P., & Henning, T. 2008a, A&A, 480, 859 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brauer, F., Henning, T., & Dullemond, C. P. 2008b, A&A, 487, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cuzzi, J. N., Hogan, R. C., Paque, J. M., & Dobrovolskis, A. R. 2001, ApJ, 546, 496 [NASA ADS] [CrossRef] [Google Scholar]
 Davis, D. R., & Ryan, E. V. 1990, Icarus, 83, 156 [NASA ADS] [CrossRef] [Google Scholar]
 Dominik, C., & Dullemond, C. P. 2008, A&A, 491, 663 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dubrulle, B., Morfill, G. E., & Sterzik, M. F. 1995, Icarus, 114, 237 [NASA ADS] [CrossRef] [Google Scholar]
 Dullemond, C. P., & Dominik, C. 2004, A&A, 421, 1075 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dullemond, C. P., van Zadelhoff, G. J., & Natta, A. 2002, A&A, 389, 464 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dzyurkevich, N. S., Flock, M., Turner, N. J., Klahr, H., & Henning, Th. 2010, A&A, in press [Google Scholar]
 Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hartmann, L., Calvet, N., Gullbring, E., & D'Alessio, P. 1998, ApJ, 495, 385 [NASA ADS] [CrossRef] [Google Scholar]
 Hueso, R., & Guillot, T. 2005, A&A, 442, 703 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Johansen, A., & Klahr, H. 2005, ApJ, 634, 1353 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Klahr, H., & Henning, T. 2006, ApJ, 636, 1121 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Kley, W., & Lin, D. N. C. 1992, ApJ, 397, 600 [NASA ADS] [CrossRef] [Google Scholar]
 Kretke, K. A., & Lin, D. N. C. 2007, ApJ, 664, L55 [NASA ADS] [CrossRef] [Google Scholar]
 Larson, R. B. 1969, MNRAS, 145, 271 [NASA ADS] [CrossRef] [Google Scholar]
 Larson, R. B. 2003, Rep. Prog. Phys, 66, 1651 [NASA ADS] [CrossRef] [Google Scholar]
 Leinhardt, Z. M., & Stewart, S. T. 2009, Icarus, 199, 542 [NASA ADS] [CrossRef] [Google Scholar]
 LyndenBell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [NASA ADS] [CrossRef] [Google Scholar]
 Myers, P. C. 2005, ApJ, 623, 280 [NASA ADS] [CrossRef] [Google Scholar]
 Nakagawa, Y., Nakazawa, K., & Hayashi, C. 1981, Icarus, 45, 517 [NASA ADS] [CrossRef] [Google Scholar]
 Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Nakamoto, T., & Nakagawa, Y. 1994, ApJ, 421, 640 [NASA ADS] [CrossRef] [Google Scholar]
 Okuzumi, S. 2009, ApJ, 698, 1122 [NASA ADS] [CrossRef] [Google Scholar]
 Ormel, C. W., & Cuzzi, J. N. 2007, A&A, 466, 413 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ormel, C. W., Spaans, M., & Tielens, A. G. G. M. 2007, A&A, 461, 215 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ormel, C. W., Paszun, D., Dominik, C., & Tielens, A. G. G. M. 2009, A&A, 502, 845 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Paszun, D., & Dominik, C. 2009, A&A, 507, 1023 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Penston, M. V. 1969, MNRAS, 145, 457 [NASA ADS] [Google Scholar]
 Ruden, S. P., & Pollack, J. B. 1991, ApJ, 375, 740 [NASA ADS] [CrossRef] [Google Scholar]
 Schäfer, C., Speith, R., & Kley, W. 2007, A&A, 470, 733 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schräpler, R., & Henning, T. 2004, ApJ, 614, 960 [NASA ADS] [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Shu, F. H. 1977, ApJ, 214, 488 [NASA ADS] [CrossRef] [Google Scholar]
 Sirono, S.I. 2004, Icarus, 167, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Teiser, J., & Wurm, G. 2009, MNRAS, 393, 1584 [NASA ADS] [CrossRef] [Google Scholar]
 Toomre, A. 1964, ApJ, 139, 1217 [NASA ADS] [CrossRef] [Google Scholar]
 Ulrich, R. K. 1976, ApJ, 210, 377 [NASA ADS] [CrossRef] [Google Scholar]
 Urpin, V. A. 1984, Soviet Astron., 28, 50 [NASA ADS] [Google Scholar]
 van Boekel, R., Waters, L. B. F. M., Dominik, C., et al. 2003, A&A, 400, L21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Visser, R., Dishoeck, E. F. V., Doty, S. D., & Dullemond, C. P. 2009, A&A, 495, 881 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2008, ApJ, 677, 1296 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2009, ApJ, 702, 1490 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J. 1977, MNRAS, 180, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J. 1980, Icarus, 44, 172 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J. 1997, Icarus, 127, 290 [NASA ADS] [CrossRef] [Google Scholar]
 Whipple, F. L. 1972, From Plasma to Planet, 211 [Google Scholar]
 Whitehouse, S. C., & Bate, M. R. 2006, MNRAS, 367, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Wurm, G., Paraskov, G., & Krauss, O. 2005, Icarus, 178, 253 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588 [NASA ADS] [CrossRef] [Google Scholar]
 Zsom, A., & Dullemond, C. P. 2008, A&A, 489, 931 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zsom, A., Omel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Footnotes
 ... F. Brauer^{}
 Deceased September 2009.
 ... regimes^{}
 It should be noted that ``Stokes regime'' refers to the regime where the drag force on a particle is described by the Stokes law  this is not directly related to the Stokes number.
 ... library^{}
 Available from www.engineers.auckland.ac.nz/ snor007
All Tables
Table 1: Parameters of the fiducial model.
All Figures
Figure 1: Total amount of infallen surface density as function of radius according to the ShuUlrich infall model (see Sect. 2.5) assuming a centrifugal radius of 8 AU (solid line), 33 AU (dashed line), and 100 AU (dashdotted line). 

Open with DEXTER  
In the text 
Figure 2: Sources of relative particle velocities considered in this model (Brownian motion velocities are not plotted) at a distance of 10 AU from the star. The turbulence parameter in this simulation was 10^{3}. It should be noted that relative azimuthal velocities do not vanish for very large and very small particles. 

Open with DEXTER  
In the text 
Figure 3: Evolution of disk surface density distribution ( top) and midplane temperature ( bottom) of the fiducial model described in 3.1. 

Open with DEXTER  
In the text 
Figure 4: Evolution of disk mass and stellar mass (solid and dashed line on left scale respectively) and accretion rate onto the star (dashdotted line on right scale). Adapted from Fig. 5 in Hueso & Guillot (2005). 

Open with DEXTER  
In the text 
Figure 5: Snapshots of the vertically integrated dust density distributions (defined in Eq. (37)) of a steady state disk as in (B08a) (left column) and of an evolving disk (fiducial model, right column). No coagulation is calculated within the evaporation radius (denoted by the dashdotted line), fragmentation is not taken into account in both simulations. The solid line shows the particle size corresponding to a Stokes number of unity. Since (see Eq. (14)) this curve in fact has the same shape as , so it reflects, as a ``bonus'', what the gas disk looks like. The radius dividing the evolving disk into accreting and expanding regions is marked by the dotted line and the arrows. Particles which are located below the dashed line have a positive flux in the radial direction due to coupling to the expanding gas disk and turbulent mixing (particles within the closed contour in the upper right plot have an inward pointing flux). 

Open with DEXTER  
In the text 
Figure 6: Comparison of the radial dependence of the dusttogas ratio for the stationary simulations (thick lines) and the evolving disk simulations (thin lines). The four panels compare the results after 1 Myr of evolution with/without fragmentation (left/right column) and with/without effects of nonaxisymmetric instabilities (top/bottom row). It can be seen that the dusttogas ratio of the evolving disk simulations is almost everywhere lower than the one of the stationary simulations. See Sect. 3.4 for more details. 

Open with DEXTER  
In the text 
Figure 7: Evolution of the dust density distribution of the fiducial model as in Fig. 5 but with fragmentation included. The dashed contour line (in the lower two panels) around the upper end of the size distribution and around small particles at >60 AU marks those parts of the distribution which have a positive (outward pointing) fluxes. 

Open with DEXTER  
In the text 
Figure 8: Evolution of the radial dependence of the dusttogas ratio in the fiducial model including fragmentation with the times corresponding to the snapshots shown in Fig. 7. The initial dusttogas ratio is taken to be 0.01. The thick dashed curve represents the result at 1 Myr of the static disk model for comparison. 

Open with DEXTER  
In the text 
Figure 9: Vertically integrated (cf. Eq. (37)) grain surface density distributions as function of grain radius at a distance of 1 AU (solid), 10 AU (dashed) and 100 AU (dotdashed) from the star. These curves represent slices through the bottom panel of Fig. 7. 

Open with DEXTER  
In the text 
Figure 10: Evolution of the dust surface density distribution of the fiducial model at 200 000 years for different drift efficiencies , without fragmentation. The solid line denotes the grain size of particles with a Stokes number of unity. Gas outside of the radius denoted by the dotted line as well as particles below the dashed line have positive radial velocities. See Sect. 3.2 for more discussion. 

Open with DEXTER  
In the text 
Figure 11: Dust grain surface density distribution as in Fig. 10 at 200 000 years but including the Stokes drag regime. The drift efficiency is set to and fragmentation is not taken into account. It can be seen that (solid line) increases with radius until about 1 AU, which facilitates the break through the drift barrier. 

Open with DEXTER  
In the text 
Figure 12: Dust surface density distributions (top row) and the according solidtogas ratio (bottom row) for the case of radiusdependent fragmentation velocity after years of evolution. In the upper row, the vertical dashed line denotes the position of the snow line at the midplane, the solid line corresponds to and the dotdashed line shows the approximate location of the fragmentation barrier according to Eq. (52). The snow line on the right plot lies further outside since viscous heating is stronger if is larger. In the bottom row, the solid line denotes the vertically integrated dusttogas ratio while the dashed line denotes the dusttogas ratio at the disk midplane. The icy dust grains outside the snow line are assumed to fragment at a critical velocity of 10 m s^{1} while particles inside the snow line fragment at 1 m s^{1}. The plots on the left and right hand side differ in the amount of turbulence in the disk ( and 10^{2}, respectively). 

Open with DEXTER  
In the text 
Figure 13: Dusttogas ratio at a distance of 1 AU from the central star as a function of time for the case of low ( , solid line) and high ( , dashed line) turbulence. 

Open with DEXTER  
In the text 
Figure 14: Evolution of the global dusttogas ratio ( top panel) and the dust disk mass ( bottom panel) for the simulations shown in Fig. 12. The solid and dashed lines correspond to the low and high turbulence case, respectively. The dashdotted line shows the evolution of the disk if a low fragmentation velocity is assumed throughout the disk. 

Open with DEXTER  
In the text 
Figure A.1: Pictographic representation of the matrix on the left hand side of Eq. (A.26) with six radial and five mass grid points. White elements are zero, dark grey elements contain contributions from radial transport of dust ( J), light gray elements contain contributions from the coagulation/fragmentation ( ), and black elements contain both contributions. The upper left and lower 5 rows represent the boundary conditions where coagulation/fragmentation is not taken into account. The matrix in the simulations would typically have a size of 15 000^{2}. 

Open with DEXTER  
In the text 
Figure A.2: Comparison between simulation result and analytical solution for growth of equal sized particles. The solid lines denote the position of the peak of the grain size distribution. In the top panel, only Brownian motion is considered as source of relative particle velocities, in the bottom panel, turbulent relative velocities are considered as well. The parameters of these simulation are T=196 K, = 1.6 , , and . 

Open with DEXTER  
In the text 
Figure B.1: Test case for only radial drift, without coagulation. The solid line denotes the massaveraged position of the radial distribution of grains for each grain size after 10^{3} years. The dashed line is the expected solution for a single particle, the dotted line denotes the initial radius of the particles. 

Open with DEXTER  
In the text 
Copyright ESO 2010