A Algorithm

Gas- and dust evolution in protoplanetary disks

Key Words.:
accretion, accretion disks – circumstellar matter – stars: formation, pre-main-sequence – infrared: stars


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 build-up phase, and the evolution of the protoplanetary disk influence growth and transport of dust.

Methods:We introduce a new version of the model of Brauer et al. (2008) in which we now include the time-dependent 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 build-up 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.


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 sub-micron 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 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][]Brauer:2008p215. 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, non-evolving gas disk model. It is known that over the duration of the planet formation process the disk itself also evolves dramatically (Lynden-Bell & Pringle 1974; Hartmann et al. 1998; Hueso & Guillot 2005), which may influence the processes of dust coagulation and fragmentation. It is the goal of this paper to introduce a combined disk-evolution and dust-evolution model which also includes additional physics: we include relative azimuthal velocities, radial dependence of fragmentation critical velocities and the Stokes-drag 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 back-reactions 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., in press).

This paper is organized as follows: Section 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 Section 3 we will compare the results of our simulations with previous steady-state 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 in-falling 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 100–200 different grain sizes, we are able to follow the detailed evolution of the dust sub-disk 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 Lynden-Bell & Pringle 1974)


is the gas surface density, the radius along the disk mid-plane 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


and the sound speed . Therefore the viscosity is written as


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 magneto-rotational instability is the most commonly accepted candidate for source of turbulence (Balbus & Hawley 1991). Values of are expected to be in the range of to some see Johansen & Klahr 2005; Dzyurkevich et al., in press. Observations confirm this range with higher probability for larger values (see Andrews & Williams 2007).

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)


Values below a critical value of describe a weakly unstable disk, which forms non-axisymmetric instabilities like spiral arms. values below unity lead to fragmentation of the disk.

The effect of these non-axisymmetric 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 self-limiting. Values above are not influenced by instabilities, values below form instabilities which increase 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 , 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.

Eq. 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 Section 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 dust-to-gas ratio in protoplanetary disks is in the order of (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 mid-plane of the disk or dead-zones and snow-lines) 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 turn-over 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 dust-to-gas 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 (i.e. the Knudsen number) and also on the particle Reynolds-number with being the gas molecular viscosity


the mean thermal velocity. The mean free path is taken to be


where denotes the mid-plane number density and  cm.

The different regimes1 are


Here 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 turn-over 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


This is rather irrelevant for the viscous evolution of the gas disk, since all values of lead to the same viscosity, but if we calculate the turn-over-eddy time, we get


The Stokes number and therefore the dust-to-gas coupling as well as the relative particle velocities strongly depend on the eddy turnover time and therefore on . In this work is taken to be 0.5 (following Cuzzi et al. 2001; Schräpler & Henning 2004) which leads to a turn-over-eddy time of


The Stokes number then becomes


The overall radial movement of dust surface density can now be described by an advection-diffusion equation,


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 dust-to-gas ratio of each size that is being smoothed out by the turbulence. The diffusive flux can therefore be written as


The ratio of gas diffusivity over dust diffusivity is called the Schmidt number. We follow Youdin & Lithwick (2007), who derived


and assume the gas diffusivity to be equal to the turbulent gas viscosity .

The second contribution to the dust flux is the advective flux,


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 sub-keplerian, 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. 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 Section 3.5.

Putting all together, the time dependent equation for the surface density of one dust species of mass 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 re-condensation 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 Section 2.6.1). All sources will be combined into one equation later which is implicitly integrated in an un-split scheme (see Section 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 in-fallen surface density as function of radius according to the Shu-Ulrich infall model (see Section 2.5) assuming a centrifugal radius of 8 AU (solid line), 33 AU (dashed line), and 100 AU (dash-dotted line).

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




The viscous heating is given by Nakamoto & Nakagawa (1994)


were denotes the turbulent gas viscosity and the Kepler frequency.

Nakamoto & Nakagawa (1994) give a solution for the mid-plane temperature with an optically thick and an optical 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, App. B),


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


only if from Eq. 25 would be larger than .

This is a thermostat effect: if rises above 1500 K, dust will evaporate, the opacity will drop and the temperature stabilizes at K. Once even the very small gas opacity is enough to get temperatures K, all the dust is evaporated and the temperature rises further.

2.4 Opacity

In the calculation of the mid-plane 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/g dust to cm/g gas by multiplying the values with the dust-to-gas 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 dust-to-gas ratio does not change with time. To calculate the opacity self-consistently, 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 back-reactions 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 inside-out 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, pressure-free 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 Shu-model 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 in-falling 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 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






as described in Ulrich (1976). Here, is given by


The centrifugal radius can therefore be approximated by (cf. Hueso & Guillot (2005))


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 in-falling matter onto the disk surface, because the finite thickness may “capture” an in-falling gas parcel before it could reach the midplane. Secondly, star formation is likely to be messier than our simple single-star axisymmetric infall model. And even in such a simplified scenario, the Shu inside-out 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

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 . It should be noted that relative azimuthal velocities do not vanish for very large and very small particles.

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 “flip-flop” back and forth between the two splitted sub-steps 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 non-splitted scheme in which the integration is done “implicitly”. B08a already use this technique to avoid flip-flopping 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.

Smoluchowski equation

The dust grain distribution , which is a function of mass , distance to the star and height above the mid-plane , describes the number of particles per cm per gram interval in particle mass. This means that the total dust density in g cm is given by


With this definition of , the coagulation/fragmentation at one point in the disk can be described by a general two-body process,


where is called the kernel. In the case of coagulation and fragmentation, this is given by


For better readability, the dependency of on radius and height above the mid-plane was omitted here. The combined coagulation/fragmentation kernel consists of four terms containing the coagulation kernel , the fragmentation kernel and the distribution of fragments after a collision .

The first two terms in Eq. 36 correspond to gain (masses and coagulate) and loss ( coagulates with ) due to grain growth.

The third term describes the fragmentation of masses and , governed by the fragmentation kernel and the fourth term describes the fact that when masses and fragment, they distribute some of their mass via fragments to smaller sizes.

The coagulation and fragmentation kernels will be described in section 2.6.3, the distribution of fragments, , 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][]Brauer:2008p215, however coagulation processes are most important near the mid-plane, which allows to approximate the kernels as being vertically constant (using the values at the mid-plane). If the vertical dust density distribution of each grain size is taken to be a Gaussian (see Section 2.6.4, Eq. 51), then the vertical integration can be done analytically, as discussed in Appendix A.2. Unlike the steady-state 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 and are related through . The total dust surface density at any radius is then given by


Defining as in Eq. 37 makes it a grid-independent 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 Monte-Carlo methods as in Ormel et al. (2007) or Zsom & Dullemond (2008), however these models have do not yet include the radial motion of dust and therefore cannot trace the global evolution of the dust disk.

Distribution of fragments

The distribution of fragments after a collision, , 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 ,


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.

Coagulation and fragmentation kernels

The coagulation kernel can be factorized into three parts,


and, similarly, the fragmentation kernel can be written as


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.

The fragmentation probability is still not well known and subject of both theoretical (Paszun & Dominik 2009; Wada et al. 2008) and experimental research (see Blum & Wurm 2008; Güttler et al. 2009). In this work, we adopt the simple recipe


with a transition width and the fragmentation speed as free parameter which is assumed to be 1 m s, following experimental work of Blum & Muench (1993) and theoretical studies of Leinhardt & Stewart (2009).

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 Figure 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


Particle growth due to Brownian relative motion is most effective for small particles.

Radial drift, as described in section 2.2 also induces relative velocities since particles of different size are differently coupled to the gas. The relative velocity is then


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 St=1 are significantly drifting, relative azimuthal velocities do not vanish for encounters between very large and vary small particles (see Figure 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 gas-dominated drag are


where 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 mid-plane, which gives

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

The dust scale height is calculated by equating the time scale for settling,


and the time scale for stirring (Dubrulle et al. 1995; Schräpler & Henning 2004; Dullemond & Dominik 2004),


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 mid-plane, however most of the dust (and hence most of the coagulation/fragmentation processes) take place near the mid-plane, 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, which corresponds to 0.06 times the break up rotation rate of the core. The disk is being built-up from inside out due to the Shu-Ulrich infall model, with the centrifugal radius being 8 AU. The parameters of our fiducial model are summarized in Table 1.

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

Figure 3 shows how the gas surface density and the mid-plane 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 small-grained 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 Section 2.1). Hueso & Guillot (2005) find a steeper, power-law 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

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 dash-dotted 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).

Table 1: Parameters of the fiducial model.

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 mid-plane temperature are the same as discussed in the previous section. The evolution of this model is visualized in Figure 5.

The contour levels in Figure 5 show the vertically integrated dust surface density distribution per logarithmic bin of grain radius, , as defined in Eq. 37. The left column of Figure 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 St=1 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 steady-state 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 Lynden-Bell & 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.

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 Figure 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 Figure 5 have positive flux in the radial direction due to the gas-coupling.

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 Figure 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 dust-to-gas ratio is constantly decreasing) while only the small particles are mixed out to large radii.

The radial dependence of the dust-to-gas ratio after 1 Myr is shown in Figure 6. These simulations show that the initial conditions of the stationary disk models (such as shown in B08a and in the left column of Figure 5) are too optimistic since they assume a constant dust-to-gas ratio at the start of their simulation throughout the disk which is not possible if the disk is being built-up 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 dust-to-gas ratio is reduced by 0.5–1.5 orders of magnitude compared to a stationary model. This effect is also discussed in Section 3.4.

3.3 Fiducial model with fragmentation

The situation changes significantly, if we take grain fragmentation into account. As discussed in Section 2.6.2, we consider two different kinds of outcomes for grain-grain 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 6: Comparison of the radial dependence of the dust-to-gas 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 non-axisymmetric instabilities (top/bottom row). It can be seen that the dust-to-gas ratio of the evolving disk simulations is almost everywhere lower than the one of the stationary simulations. See Section 3.4 for more details.
Figure 7: Evolution of the dust density distribution of the fiducial model as in Figure 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  AU marks those parts of the distribution which have a positive (outward pointing) fluxes.

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 high-velocity 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 Figure 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 semi-equilibrium 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 Figure 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 Figure 8. The dust-to-gas ratio at 1 AU after 1 Myr is increased by 25% over the value of in-falling 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 Figure 8 shows as comparison the dust-to-gas ratio of the stationary disk model (cf. left column of Figure 5) after 1 Myr, which starts with a radially constant initial dust-to-gas ratio of 0.01.

Figure 9 shows the semi-equilibrium 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 semi-equilibrium distributions is always the same: a power-law 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 equal-sized collisions than for collisions with particles of very different sizes (see Figure 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 dust-to-gas ratio in the fiducial model including fragmentation with the times corresponding to the snapshots shown in Figure 7. The initial dust-to-gas ratio is taken to be 0.01. The thick dashed curve represents the result at 1 Myr of the static disk model for comparison.
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 (dot-dashed) from the star. These curves represent slices through the bottom panel of Figure 7.

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. Figure 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 Figure 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 mid-plane 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 dust-to-gas ratios are being shown in Figure 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 dust-to-disk 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 build-up 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 dust-to-gas ratio (see left column of Figure 6).

It can be seen, that in all simulations, the dust-to-gas ratio is lower than in the stationary disk model. The trend in the upper right panel in Figure 6 suggests that for a centrifugal radius of 100 AU and the enhanced radial transport by GI effects might have a higher dust-to-gas 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 disk-star 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 .

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 section 3.2 for more discussion.

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 dust-to-gas ratio and and are the drift and growth timescales respectively. The parameter describes how much more efficient growth around St=1 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 dust-to-gas 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 St=1 particles also depends on the temperature (via the pressure gradient). A simple approximation from Eq. 53 with a 0.5  star and a dust-to-gas ratio of 0.01 gives


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 angle-dependent radiative transfer. They show that the mid-plane 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 meter-sized particles.

Meridional flows (e.g., Urpin 1984; Kley & Lin 1992) might also seem interesting in this context, however they do not directly influence the radial drift efficiency but rather reverse the gas-drag effect. This might be important for small particles (which, however are not strongly settling to the mid-plane) but for St=1 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 St=1 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 11: Dust grain surface density distribution as in Figure 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.

In the Epstein regime, the size corresponding to St=1 is proportional to the gas surface density


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

Figure 12: Dust surface density distributions (top row) and the according solid-to-gas ratio (bottom row) for the case of radius-dependent fragmentation velocity after 2 years of evolution. In the upper row, the vertical dashed line denotes the position of the snow line at the mid-plane, the solid line corresponds to and the dot-dashed 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 dust-to-gas ratio while the dashed line denotes the dust-to-gas ratio at the disk mid-plane. The icy dust grains outside the snow line are assumed to fragment at a critical velocity of 10 m s while particles inside the snow line fragment at 1 m s. The plots on the left and right hand side differ in the amount of turbulence in the disk ( and , respectively).

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


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

As shown by Ormel & Cuzzi (2007), the largest relative velocities are of the order of


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


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, which could already seen in the simulations by Brauer et al. (2008b), who have simulated particle growth near the snow-line.

Even in the case of very low turbulence, relative azimuthal velocities of large () and small grains () are of the order of 30 m s, which means that large particles are constantly being ’sand-blasted’ by small grains. The only way of reducing these velocities significantly is decreasing the pressure gradient (see Equations 20 and 46).

Another possibility to overcome this problem would be if high-velocity 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 Brauer et al. 2008b; Kretke & Lin 2007): 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 (), this means that the fragmentation velocity must be smaller than about 0.5–5 m s. 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 (for a review, see Blum & Wurm 2008). The composition of particles outside the snow-line 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 snow-line to be 1 m s, outside the snow-line to be 10 m s. 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 dust-to-gas 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 Figure 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 Figure 12. Particles outside the snow-line 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 snow-line.

Strong drift outside the snow line and weaker radial drift inside the snow line cause the dust-to-gas ratio within the snow line to increase significantly (see bottom row of Figure 12): in the case of , the dust-to-gas 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 dust-to-gas ratio is not as high as for a solar mass star (dust-to-gas ratio of 0.27–0.20 from 0.2–4 AU).

This effect is not as significant in the case of stronger turbulence, where the maximum dust-to-gas ratio is around 0.027. The evolution of the dust-to-gas ratio at a distance of 1 AU from the star is plotted for both cases in Figure 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 dot-dashed line in Figure 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 Figure 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 Figure 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 Figure 14) seems counter-intuitive: 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 dust-to-gas ratio of the disks (top panel of Figure 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 Figure 14. It can be seen that the dust mass is retained at its initial value for much longer timescales.

Figure 13: Dust-to-gas 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.
Figure 14: Evolution of the global dust-to-gas ratio (top panel) and the dust disk mass (bottom panel) for the simulations shown in Figure 12. The solid and dashed lines correspond to the low and high turbulence case, respectively. The dash-dotted line shows the evolution of the disk if a low fragmentation velocity is assumed throughout the disk.

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 re-implemented 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, un-split 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 mass-conserving (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 are needed to be able to overcome the fragmentation barrier in the case of turbulent relative velocities.

This model includes also the initial build-up phase of the disk, which is still a very poorly understood phase of disk evolution. We use the Shu-Ulrich infall model which represents a strong simplification. However, the following novel findings of this work do not or only weakly depend on the build-up phase of the disk:

  • Apart from an increased dust-to-gas ratio (?)as in][]Brauer:2008p215, 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 meter-size 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 at 1 AU or 330 g cm 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 ’sand-blasted’ by small grains (if they are present) which (in our prescription of fragmentation) causes erosion and stops grain-growth 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., in press). 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 net-growth by high velocity impacts of small particles onto larger bodies.

  • Higher tensile strengths of particles outside the snow-line allows particles to grow to larger sizes, which are more strongly affected by radial drift. Particles therefore drift from outside the snow-line 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 dust-to-gas ratio inside the snow-line 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.

The following results depend on the build-up phase of the disk. However unless the collapse of the parent cloud is not inside-out or so fast to cause disk fragmentation, we expect only slight alteration of the results:

  • Disk spreading causes small particles ( m) to be transported outward at radii of 60–190 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 inside-out build-up 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 Section A.3 we will then describe how this model treats the radial and the size evolution of dust in an un-split, fully implicit way.

a.1 Advection-diffusion 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 advection-diffusion equation,


which can be adapted to both, Eq. 1 and Eq. 21 by proper choice of parameters , , , , and .

We use a flux-conserving donor-cell scheme which is implicit in . The time derivative in Eq. 60, written in a discretized way becomes


where denotes time-dimension and denotes space-dimension.

The advective part is discretized as


where and are the future flux and the surface between cell and cell and is the volume of cell . The advective interface fluxes can then be written as


The diffusive interface flux becomes


Working out these equations and separating the values of leads to


with the coefficients


Eq. 67 can now be solved by any matrix-solver, but since it is a tri-diagonal matrix, the fastest analytical way to solve it is by forward elimination/backward substitution.

It should be noted that Eq. 60 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 gives


where the dust particle number density is


If we assume that the coagulation and fragmentation kernels are constant in and that the vertical distribution of grains is a Gaussian with a scale height according to Eq. 51,


we can vertically integrate the coagulation/fragmentation equation.




Discretizing Eq. 72 also in radial direction and rewriting it in terms of the quantity