EDP Sciences
Open Access
Issue
A&A
Volume 513, April 2010
Article Number A79
Number of page(s) 21
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/200913731
Published online 30 April 2010
A&A 513, A79 (2010)

Gas- and dust evolution in protoplanetary disks

T. Birnstiel - C. P. Dullemond - F. Brauer[*]

Junior Research Group at the Max-Planck-Institut 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 build-up 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 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.

Key words: accretion, accretion disks - circumstellar matter - stars: formation - stars: pre-main 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 $\mu$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 1027 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, 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; 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 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. 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 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,

\begin{displaymath}%
\frac{\partial\ensuremath{\Sigma_{\rm g}}\ }{\partial t} - ...
...\ensuremath{\Sigma_{\rm g}}\ r ~ u_{\rm g}\right) = S_{\rm g},
\end{displaymath} (1)

where the gas radial velocity $u_{\rm g}$ is given by (see Lynden-Bell & Pringle 1974)

\begin{displaymath}%
u_{\rm g} = - \frac{3}{\ensuremath{\Sigma_{\rm g}}\ \sqrt{r...
...ft( \ensuremath{\Sigma_{\rm g}}\ \nu_{\rm g} \sqrt{r} \right).
\end{displaymath} (2)

$\ensuremath{\Sigma_{\rm g}} = \int_{-\infty}^{\infty} \ensuremath{\rho_{\rm g}} (z) \ensuremath{{\rm d}} z$ is the gas surface density, r the radius along the disk mid-plane and \ensuremath{\nu_{\rm g}} the gas viscosity. The source term on the right hand side of Eq. (1), denoted by $S_{\rm g}$ 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  $\ensuremath{H_{\rm p}} $

\begin{displaymath}%
\ensuremath{H_{\rm p}} = \frac{\ensuremath{c_{\rm s}} }{\ensuremath{\Omega_{\rm k}} }
\end{displaymath} (3)

and the sound speed $\ensuremath{c_{\rm s}} $. Therefore the viscosity is written as

\begin{displaymath}%
\ensuremath{\nu_{\rm g}} = \alpha \: \ensuremath{c_{\rm s}}\: \ensuremath{H_{\rm p}} ,
\end{displaymath} (4)

where $\alpha $ is the turbulence parameter and $\alpha \leq 1$.

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 $\alpha $ 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)

\begin{displaymath}%
Q = \frac{\ensuremath{c_{\rm s}}\ensuremath{\Omega_{\rm k}} }{\pi ~G ~ \ensuremath{\Sigma_{\rm g}} }\cdot
\end{displaymath} (5)

Values below a critical value of $Q_{\rm cr} = 2$ describe a weakly unstable disk, which forms non-axisymmetric instabilities like spiral arms. Q 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 $Q_{\rm cr}$ are not influenced by instabilities, values below  $Q_{\rm cr}$ form instabilities which increase Q until the disk is marginally stable again. Our model approximates this mechanism by increasing the turbulence parameter $\alpha $ according to the recipe of Armitage et al. (2001),

\begin{displaymath}%
\alpha(r) = \alpha_0 + 0.01 \left( \left(\frac{Q_{\rm cr}}{\min(Q(r),Q_{\rm cr})}\right)^2 - 1 \right),
\end{displaymath} (6)

where $\alpha_0$ 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 $\alpha = 0.1$ 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 ${\propto}\Delta r^2$ (where $\Delta r$ 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 dust-to-gas 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 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

\begin{displaymath}%
{\rm St}= \frac{\tau_{\rm s}}{\tau_{\rm ed}},
\end{displaymath} (7)

where $\tau_{\rm ed}$ is the eddy turn-over time and $\tau_{\rm s}$ 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  $\lambda_{\rm mfp}$ of the gas molecules and the dust particle size a (i.e. the Knudsen number) and also on the particle Reynolds-number $Re = {2 a u/\nu_{\rm mol}}$ with  $\nu_{\rm mol}$ being the gas molecular viscosity

\begin{displaymath}%
\nu_{\rm mol} = \frac{1}{2}~ \bar u ~ \lambda_{\rm mfp},
\end{displaymath} (8)

$\bar u$ the mean thermal velocity. The mean free path is taken to be

\begin{displaymath}%
\lambda_{\rm mfp} = \frac{1}{n\:\sigma_{{\rm H}_2}}
\end{displaymath} (9)

where n denotes the mid-plane number density and $\sigma_{{\rm H}_2} =2~\ensuremath{\times 10^{-15}} $ cm2.

The different regimes[*] are

\begin{displaymath}%
\tau_{\rm s} =
\left\{
\begin{array}{lll}
\frac{\rho_{\rm ...
...{\rm for }~Re>800& {\rm Stokes~regime~3}\\
\end{array}\right.
\end{displaymath} (10)

Here u denotes the velocity of the dust particle with respect to the gas, $\bar u = c_{\rm s} \sqrt{{\pi}/{8}}$ denotes the mean thermal velocity of the gas molecules, $\rho _{\rm s}$ is the solid density of the particles and  $\rho_{\rm g}$ 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 $L_{\rm c}$ and a characteristic velocity scale $V_{\rm c}$ 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

\begin{displaymath}%
\ensuremath{\nu_{\rm g}} = (\alpha^{q} V_{\rm c}) \cdot (\alpha^{1-q} L_{\rm c}).
\end{displaymath} (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 turn-over-eddy time, we get

\begin{displaymath}%
\tau_{\rm ed} = \frac{2\pi L_{\rm c}}{V_{\rm c}} = \alpha^{1-2q}\: \frac{1}{\ensuremath{\Omega_{\rm k}} }\cdot
\end{displaymath} (12)

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 q. In this work q is taken to be 0.5 (following Schräpler & Henning 2004; Cuzzi et al. 2001) which leads to a turn-over-eddy time of

\begin{displaymath}%
\tau_{\rm ed} = \frac{1}{\ensuremath{\Omega_{\rm k}} }\cdot
\end{displaymath} (13)

The Stokes number then becomes

\begin{displaymath}%
{\rm St}= \frac{\rho_{\rm s}\cdot a}{\ensuremath{\Sigma_{\r...
...frac{\pi}{2} \qquad {\rm for }~a<\frac{9}{4}\lambda_{\rm mfp}.
\end{displaymath} (14)

The overall radial movement of dust surface density  $\ensuremath{\Sigma_{\rm d}} $ can now be described by an advection-diffusion equation,

\begin{displaymath}%
\frac{\partial \ensuremath{\Sigma_{\rm d}} }{\partial t} + ...
...frac{\partial }{\partial r} \Bigl( r ~ F_{\rm tot} \Bigr) = 0,
\end{displaymath} (15)

where the total Flux, $F_{\rm tot}$ 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

\begin{displaymath}%
F_{\rm diff} = - D_{\rm d} \: \frac{\partial }{\partial r} ...
...h{\Sigma_{\rm g}} }\right) \cdot \ensuremath{\Sigma_{\rm g}} .
\end{displaymath} (16)

The ratio of gas diffusivity $D_{\rm g}$ over dust diffusivity $D_{\rm d}$ is called the Schmidt number. We follow Youdin & Lithwick (2007), who derived

\begin{displaymath}%
{\rm Sc}\equiv \frac{D_{\rm g}}{D_{\rm d}} = {1+{\rm St}^2},
\end{displaymath} (17)

and assume the gas diffusivity to be equal to the turbulent gas viscosity  \ensuremath{\nu_{\rm g}}.

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

\begin{displaymath}%
F_{\rm adv} = \ensuremath{\Sigma_{\rm d}}\cdot u_{\rm r},
\end{displaymath} (18)

where $u_{\rm r}$ is the radial velocity of the dust. There are two contributions to it,

\begin{displaymath}%
u_{\rm r} = \frac{u_{\rm g}}{1 + {\rm St}^2} - \frac{2 u_{n}}{{\rm St}+ {\rm St}^{-1}}\cdot
\end{displaymath} (19)

The first term is a drag term which comes from the radial movement of the gas which moves with a radial velocity of $u_{\rm g}$, 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). un denotes the maximum drift velocity of a particle,

\begin{displaymath}%
u_{n} = - E_{\rm d} \cdot \frac{ \frac{\partial P_{\rm g}}{\partial r}}{2 \: \rho_{\rm g} \: \ensuremath{\Omega_{\rm k}} },
\end{displaymath} (20)

which has been derived by Weidenschilling (1977). Here, we introduced a radial drift efficiency parameter $E_{\rm d}$. 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 mi is given by

\begin{displaymath}%
\frac{\partial \ensuremath{\Sigma_{\rm d}} ^i}{\partial t} ...
...ensuremath{\Sigma_{\rm g}}\right] \right\rbrace = S_{\rm d}^i,
\end{displaymath} (21)

where we have included a source term $S^i_{\rm d}$ 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 Sect. 2.6.1). All sources will be combined into one equation later which is implicitly integrated in an un-split 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.

\begin{figure}
\par\includegraphics[width=8.5cm,clip]{plots/13731fg01}
\end{figure} Figure 1:

Total amount of in-fallen surface density as function of radius according to the Shu-Ulrich infall model (see Sect. 2.5) assuming a centrifugal radius of 8 AU (solid line), 33 AU (dashed line), and 100 AU (dash-dotted 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 ( $\ensuremath{H_{\rm p}} ~(r)/r\ll r$) 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

\begin{displaymath}%
\ensuremath{\rho_{\rm g}} (z) = \rho_0 \: \exp\left( - \frac{1}{2}\: \frac{z^2}{\ensuremath{H_{\rm p}} ^2} \right),
\end{displaymath} (22)

where

\begin{displaymath}%
\rho_0 = \frac{\ensuremath{\Sigma_{\rm g}} }{\sqrt{2\:\pi} \: \ensuremath{H_{\rm p}} }\cdot
\end{displaymath} (23)

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

\begin{displaymath}%
Q_+ = \ensuremath{\Sigma_{\rm g}}\: \ensuremath{\nu_{\rm g}...
...{\partial \ensuremath{\Omega_{\rm k}} }{\partial r} \right)^2,
\end{displaymath} (24)

where $\ensuremath{\nu_{\rm g}} $ denotes the turbulent gas viscosity and $\ensuremath{\Omega_{\rm k}} $ the Kepler frequency.

Nakamoto & Nakagawa (1994) give a solution for the mid-plane temperature with an optically thick and an optically thin contribution,

\begin{displaymath}%
{\ensuremath{T_{\rm mid}} ^4} = \frac{9}{8 \ensuremath{\sig...
...th{c_{\rm s}} ^2\: \ensuremath{\Omega_{\rm k}} + T_{\rm irr}^4
\end{displaymath} (25)

where we used $\ensuremath{\nu_{\rm g}} = \alpha~ \ensuremath{c_{\rm s}} ~ \ensuremath{H_{\rm p}} $ and the approximation $\tau_{\rm R/P} = \kappa_{\rm R/P} ~\ensuremath{\Sigma_{\rm g}} /2$. \ensuremath{\kappa_{\rm R}} and \ensuremath{\kappa_{\rm P}} are Rosseland and Planck mean opacities which will be discussed in the next section.

$T_{\rm irr}$ contains contributions due to stellar or external irradiation. Here, we use a formula derived by Ruden & Pollack (see Ruden & Pollack 1991, Appendix B),

\begin{displaymath}%
T_{\rm irr} = T_\star \cdot \left[ \frac{2}{3 \pi} \left(\f...
..._{\rm p}} }{{\rm d}\!\ln r} - 1 \right) \right]^{\frac{1}{4}},
\end{displaymath} (26)

with a fixed ${\rm d ln}\ensuremath{H_{\rm p}} /{\rm d ln} r = 9/7$, 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 ${\sim}1500$ 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  $T_{\rm mid}$ is approximated by

\begin{displaymath}%
T_{\rm mid} = {\rm max}(T_{\rm gas},T_{\rm evap}),
\end{displaymath} (27)

only if $T_{\rm mid}$ from Eq. (25) would be larger than  $T_{\rm evap}$.

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 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 cm2/g dust to cm2/g gas by multiplying the values with the dust-to-gas ratio  $\epsilon_0$, 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 ( $\dot m_{\rm in} \approx 47\: c_{\rm s}^3/G$ and $0.975\: c_{\rm s}^3/G$ 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 $13\: c_{\rm s}^3/G$ 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 $\dot m_{\rm in} = 7.07 \ensuremath{c_{\rm s}} ^3/G$. 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 $\dot m_{\rm in}=0.975 \: c_{\rm s}^3/G$ which stays constant throughout the collapse.

We assume the singular isothermal sphere of mass $M_{\rm cloud}$ to be in solid body rotation  $\Omega_{\rm s}$. If in-falling material is conserving its angular momentum, all matter from a shell of radius $r_{\rm s}$ will fall onto the star and disk system (with mass  $m_{\rm cent}(t)$) within the so called centrifugal radius,

\begin{displaymath}%
r_{\rm centr}(t) = \frac{\Omega_{\rm s}^2\: r_{\rm s}^4}{G \: m_{\rm cent}(t)},
\end{displaymath} (28)

where G is the gravitational constant and $r_{\rm s} = 0.975\cdot c_{\rm s}\:t /2$. 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

\begin{displaymath}%
\dot\Sigma_{\rm d}(r,t) = 2\;\rho_1(r,t) \cdot u_{z}(r,t),
\end{displaymath} (29)

where

\begin{displaymath}%
u_{z}(r,t) = \sqrt{\frac{G \: m_{\rm cent}(t)}{r}} \cdot \mu,
\end{displaymath} (30)

and

\begin{displaymath}%
\rho_1(r,t) = \frac{\dot m_{\rm in}}{8 \pi \sqrt{G \: m_{\r...
...r^3}} \cdot \frac{r}{r_{\rm centr}(r,t)} \cdot \frac{1}{\mu^2}
\end{displaymath} (31)

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

\begin{displaymath}%
\mu = \sqrt{1-r/r_{\rm centr}(r,t)}.
\end{displaymath} (32)

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

\begin{displaymath}%
r_{\rm centri}(t) \!\simeq\!
1.4\! \left( \frac{\Omega_{\r...
... s}} }{3\!\times\! 10^4\;{\rm cm~s}^{-1}}\right)^{-8}{\rm AU}.
\end{displaymath} (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 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

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 ${\rm St}=1$ particles: coagulation of smaller sizes restores the population around ${\rm St}=1$, 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.

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 mid-plane z, describes the number of particles per cm3 per gram interval in particle mass. This means that the total dust density in g cm-3 is given by

\begin{displaymath}%
\rho(r,z)=\int_0^\infty n(m,r,z) \cdot m \; \ensuremath{{\rm d}} {m}.
\end{displaymath} (34)

With this definition of n(m,r,z), the coagulation/fragmentation at one point in the disk can be described by a general two-body process,
                     $\displaystyle %
\frac{\partial}{\partial t} n(m,r,z)$ = $\displaystyle \iint_0^\infty M(m,m',m'')$  
    $\displaystyle \times~ n(m',r,z) \cdot n(m'',r,z) ~\ensuremath{{\rm d}} m' ~\ensuremath{{\rm d}} m'',$ (35)

where M(m,m',m'') is called the kernel. In the case of coagulation and fragmentation, this is given by
             M(m,m',m'') = $\displaystyle \frac{1}{2}K(m',m'') \cdot \delta(m'+m''-m)$  
    $\displaystyle - ~ K(m',m'') \cdot \delta(m''-m)$  
    $\displaystyle +~\frac{1}{2}L(m',m'') \cdot S(m,m',m'')$  
    $\displaystyle -~ L(m',m'') \cdot \delta(m-m'').$ (36)

For better readability, the dependency of M on radius and height above the mid-plane 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.

\begin{figure}
\par\includegraphics[width=16cm,clip]{plots/13731fg02}
\end{figure} 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 $\alpha $ 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 m-m' 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 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 Sect. 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,  $\sigma(r,a)$ as

\begin{displaymath}%
\sigma(r,a) = \int_{-\infty}^\infty n(r,z,a) \cdot m\cdot a\; \ensuremath{{\rm d}} {z},
\end{displaymath} (37)

where n(r,z,a) and n(r,z,m) are related through $m=4\pi/3 \rho_{\rm s}a^3$. The total dust surface density at any radius is then given by

\begin{displaymath}%
\Sigma_{\rm d}(r) = \int_0^\infty \sigma(r,a)\; {\rm dln}a.
\end{displaymath} (38)

Defining $\sigma(r,a)$ as in Eq. (37) makes it a grid-independent density unlike the mass integrated over each numerical bin. This way, all plots of  $\sigma(r,a)$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 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,

\begin{displaymath}%
n(m) {\rm d}m \propto m^{-\xi} {\rm d}m.
\end{displaymath} (39)

The value $\xi$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 $\xi=1.83$.

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 $m_{\rm s}$,

\begin{displaymath}%
m_{\rm out} = \chi \: m_{\rm s}.
\end{displaymath} (40)

The mass of the smaller particle plus the mass excavated from the larger one is then distributed to masses smaller than $m_{\rm s}$ according to the distribution described by Eq. (39). In this work, we follow (B08a) by assuming $\chi$ to be unity, unless otherwise noted.

2.6.3 Coagulation and fragmentation kernels

The coagulation kernel K(m1,m2) can be factorized into three parts,

\begin{displaymath}%
K(m_1,m_2) = \Delta u(m_1,m_2) \: \sigma_{\rm geo}(m_1,m_2) \: p_{\rm c},
\end{displaymath} (41)

and, similarly, the fragmentation kernel can be written as

\begin{displaymath}%
L(m_1,m_2) = \Delta u(m_1,m_2) \: \sigma_{\rm geo}(m_1,m_2) \: p_{\rm f}.
\end{displaymath} (42)

Here, $\Delta u(m_1,m_2)$ denotes the relative velocity of the two particles, $\sigma_{\rm geo}(m_1,m_2)$ is the geometrical cross section of the collision and $p_{\rm c}$ and $p_{\rm f}$ 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 $\rho_{\rm s}=1.6$ 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

\begin{displaymath}%
p_{\rm f} = \left\{
\begin{array}{ll}
0& {\rm if }~\Delta u...
...u_{\rm f}} -\Delta u}{\delta u}& {\rm else}
\end{array}\right.
\end{displaymath} (43)

with a transition width $\delta u$ and the fragmentation speed  \ensuremath{u_{\rm f}} 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

\begin{displaymath}%
\Delta u_{\rm BM} = \sqrt{\frac{8 k_{\rm B} \: T (m_1 + m_2)}{\pi \: m_1 \: m_2 }}\cdot
\end{displaymath} (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

\begin{displaymath}%
\Delta u_{\rm RD} = \left\vert u_{\rm r}(m_1) - u_{\rm r}(m_2) \right\vert,
\end{displaymath} (45)

where the radial velocity of the dust, $u_{\rm r}$ 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 ${\rm St}=1$ 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 gas-dominated drag are

\begin{displaymath}%
u_\varphi = \left\vert u_{n} \cdot \left( \frac{1}{1+{\rm St}_1^2} - \frac{1}{1+{\rm St}_2^2} \right) \right\vert,
\end{displaymath} (46)

where un 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 $u_{\rm sett} = - z \: \ensuremath{\Omega_{\rm k}}\: {\rm St}$ which can be derived by equating the frictional force $F_{\rm fric} = - m \: u / t_{\rm fric}$ and the vertical component of the gravity force, $F_{\rm grav} = - m \: z \: \ensuremath{\Omega_{\rm k}} ^2$. To limit the settling speed to velocities smaller than half the vertically projected Kepler velocity, we use

\begin{displaymath}%
u_{\rm sett} = - z \: \ensuremath{\Omega_{\rm k}}\:\min\left({\rm St},\frac{1}{2}\right)
\end{displaymath} (47)

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

\begin{displaymath}%
\Delta u_{\rm DS} = \left\vert h_i \cdot \min({\rm St}_i,1/...
... St}_j,1/2) \right\vert \: \cdot \ensuremath{\Omega_{\rm k}} .
\end{displaymath} (48)

Table 1:   Parameters of the fiducial model.

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

\begin{displaymath}%
t_{\rm sett} = \frac{z}{u_{\rm sett}}
\end{displaymath} (49)

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

\begin{displaymath}%
t_{\rm stir} = \frac{z^2}{D_{\rm d}}\cdot
\end{displaymath} (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

\begin{displaymath}%
h_{\rm d} = \ensuremath{H_{\rm p}}\cdot \min\left(1, \sqrt{\frac{\alpha}{\min({\rm St},1/2)(1+{\rm St}^2)}} \right)\cdot
\end{displaymath} (51)

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 $7\ensuremath{\times 10^{-14}} $ s-1, 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.

\begin{figure}
\par\includegraphics[width=7.5cm,clip]{plots/13731fg03}
\end{figure} Figure 3:

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

Open with DEXTER

\begin{figure}
\par\includegraphics[width=8.5cm,clip]{plots/13731fg04}
\end{figure} 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 Fig. 5 in Hueso & Guillot (2005).

Open with DEXTER

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 $\alpha $ after the infall stops (see Sect. 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

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 Fig. 5.

The contour levels in Fig. 5 show the vertically integrated dust surface density distribution per logarithmic bin of grain radius,  $\sigma(r,a)$, 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 ${\rm St}=1$ at the given radius. This grain size will be called $a_{{\rm St}=1}$ hereafter. In the Epstein regime, $a_{{\rm St}=1}$ 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 $R_\pm$ 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.

\begin{figure}
\par\includegraphics[width=7.8cm,clip]{plots/13731fg05a}\includegraphics[width=7.8cm,clip]{plots/13731fg05b}
\end{figure} 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 $a_{\rm St=1} \propto \ensuremath {\Sigma _{\rm g}} $ (see Eq. (14)) this curve in fact has the same shape as  $\ensuremath {\Sigma _{\rm g}} (r)$, 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 $R_\pm$ 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 gas-coupling.

\begin{figure}
\par\includegraphics[width=16cm,clip]{plots/13731fg06}
\end{figure} 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 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 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 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 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 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 grain-grain collisions with relative velocities larger than the fragmentation velocity  \ensuremath{u_{\rm f}}: cratering (if the masses differ by less than one order of magnitude) and complete fragmentation (otherwise).

\begin{figure}
\par\includegraphics[width=7.5cm,clip]{plots/13731fg07}
\end{figure} 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 ${\rm St}\approx 10^{-3}$, 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  \ensuremath{u_{\rm f}} with the approximate relative velocities of the particles,

\begin{displaymath}%
{\rm St}_{\rm max} \simeq \frac{\ensuremath{u_{\rm f}} ^2}{2\alpha~\ensuremath{c_{\rm s}} ^2}\cdot
\end{displaymath} (52)

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 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 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 Fig. 8.

The approximate maximum Stokes number, defined in Eq. (52), is inversely proportional to the temperature (since relative velocities are proportional to  \ensuremath{c_{\rm s}}), 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 $\alpha/2$, are being dominated by radial drift.

Due to the high temperatures below $\sim$5 AU (caused by viscous heating), ${\rm St}_{\rm max}$ 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 $\alpha/2$. The gas within about 5 AU is therefore enriched in dust, as seen in Fig. 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 Fig. 8 shows as comparison the dust-to-gas ratio of the stationary disk model (cf. left column of Fig. 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 $\sigma$) for small grains and a peak at some grain size  $a_{\rm max}$, 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 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  $a_{\rm max}$, causing a ``traffic jam'' close to the fragmentation barrier. The peak in the distribution represents this traffic jam.

\begin{figure}
\par\includegraphics[width=8.5cm,clip]{plots/13731fg08}
\end{figure} 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 Fig. 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.

Open with DEXTER

\begin{figure}
\par\includegraphics[width=8.5cm,clip]{plots/13731fg09}
\end{figure} 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 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 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  $r_{\rm centr}$, 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 $\alpha $ 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 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 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 ${\sim}10$ 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 Fig. 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 Fig. 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 $M_{\rm disk}/M_\star \la 0.1$.

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

\begin{displaymath}%
\frac{\tau_{\rm g}}{\tau_{\rm d}} = \frac{1}{\epsilon_0} \l...
...c{\ensuremath{H_{\rm p}} }{r} \right)^2 \leq \frac{1}{\gamma},
\end{displaymath} (53)

where $\epsilon_0$ is the dust-to-gas ratio and $\tau_{\rm d}$ and $\tau_{\rm g}$ are the drift and growth timescales respectively. The parameter $\gamma$ describes how much more efficient growth around ${\rm 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 $\gamma$ 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 ${\rm St}=1$ particles also depends on the temperature (via the pressure gradient). A simple approximation from Eq. (53) with a 0.5 $M_\odot$ star and a dust-to-gas ratio of 0.01 gives

\begin{displaymath}%
T < 103~{\rm K } \left(\frac{r}{{\rm AU}}\right)^{-1},
\end{displaymath} (54)

which means that particles should be able to break through the drift barrier at 1 AU if the temperature is below $\sim$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 $E_{\rm d}$ (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., 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 gas-drag effect. This might be important for small particles (which, however are not strongly settling to the mid-plane) but for ${\rm St}=1$ particles, $\alpha $ would have to be extremely high to have significant influence: even $\alpha = 0.1$ 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 $a_{{\rm St}=1}$. Since particles grow while drifting, the particle size corresponding to ${\rm 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  $a_{{\rm St}=1}$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.

\begin{figure}
\par\includegraphics[width=8.5cm,clip]{plots/13731fg10}
\end{figure} Figure 10:

Evolution of the dust surface density distribution of the fiducial model at 200 000 years for different drift efficiencies $E_{\rm d}$, without fragmentation. The solid line denotes the grain size $a_{{\rm St}=1}$ 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

\begin{figure}
\par\includegraphics[width=8.5cm,clip]{plots/13731fg11}
\end{figure} 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 $E_{\rm d} = 0.75$ and fragmentation is not taken into account. It can be seen that $a_{\rm {\rm St}=1}$ (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 ${\rm St}=1$ is proportional to the gas surface density

\begin{displaymath}%
a_{{\rm St}=1} = \frac{2 \ensuremath{\Sigma_{\rm g}} }{\pi \rho_{\rm s}},
\end{displaymath} (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 $E_{\rm d}$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 $a_{{\rm St}=1}$ 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 $E_{\rm d}=0.5$. 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 $a_{\rm {\rm St}=1}$ can be achieved. The Epstein drag regime for particles with Stokes number of unity is only valid if

\begin{displaymath}%
\ensuremath{\Sigma_{\rm g}}\la 108\frac{{\rm g}}{{\rm cm}^2...
..._{\rm s}}{1.6\frac{{\rm g}}{{\rm cm}^3}}\right)^{\frac{1}{2}},
\end{displaymath} (56)

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

\begin{displaymath}%
{\rm St}= \frac{\sqrt{2\pi}}{9} \frac{\rho_{\rm s}\; \sigma...
...mu\; \ensuremath{m_{\rm p}} } \: \ensuremath{H_{\rm p}} ^{-1},
\end{displaymath} (57)

with $\sigma_{{\rm H}_2}$ being the cross section of molecular hydrogen. Interestingly, $a_{{\rm St}=1}$ is independent of  $\ensuremath{\Sigma_{\rm g}} $ 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 $E_{\rm d} \la 0.75$. 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.

\begin{figure}
\par\includegraphics[width=18cm,clip]{plots/13731fg12}
\vspace*{4mm}
\end{figure} 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\ensuremath {\times 10^{5}} $ 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 $a_{{\rm St}=1}$ 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 $\alpha $ 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-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 ( $\alpha = 10^{-3}$ and 10-2, respectively).

Open with DEXTER

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

\begin{displaymath}%
\Delta u_{\rm max} \simeq \sqrt{2 ~ \alpha} ~ \ensuremath{c_{\rm s}} .
\end{displaymath} (58)

If particles should be able to break through the fragmentation barrier, then they need to survive these large relative velocities, meaning that  $\Delta u_{\rm max}$ has to be smaller than the fragmentation velocity of the particles, or

\begin{displaymath}%
\frac{\ensuremath{u_{\rm f}} }{\ensuremath{c_{\rm s}} } \ga \sqrt{2 ~ \alpha}.
\end{displaymath} (59)

This condition is hard to fulfill with reasonable fragmentation velocities, unless $\alpha $ is very small. E.g., for $\alpha = 10^{-5}$ 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 snow-line.

Even in the case of very low turbulence, relative azimuthal velocities of large ( ${\rm St}\ga 1$) and small grains ( ${\rm St}\ll 1$) are of the order of 30 m s-1, 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 Eqs. (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 Kretke & Lin 2007; Brauer et al. 2008b): shutting of MRI leads to low values of $\alpha $, 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 $\alpha $ ( 10-3-10-2), this means that the fragmentation velocity must be smaller than about 0.5-5 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 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-1, outside the snow-line 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 $M_\odot$ 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 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 ${\rm St}_{\rm max} \propto a_{\rm max}$ 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 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 Fig. 12): in the case of $\alpha = 10^{-3}$, 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 $M_\odot$) 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 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 dot-dashed 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 $\alpha $ 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 $a_{{\rm St}=1}$. 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 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 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.

\begin{figure}
\par\includegraphics[width=8.5cm,clip]{plots/13731fg13}
\end{figure} 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 ( $\alpha = 10^{-3}$, solid line) and high ( $\alpha =10^{-2}$, dashed line) turbulence.

Open with DEXTER

\begin{figure}
\par\includegraphics[width=8.5cm,clip]{plots/13731fg14}
\end{figure} Figure 14:

Evolution of the global dust-to-gas 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 dash-dotted 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 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 ( $\alpha\la 10^{-5}$) 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 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 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 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-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 ${\rm St}>1$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. 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 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 (${\la}10$ $\mu$m) to be transported outward at radii of $\sim$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.

Acknowledgements
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 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,

\begin{displaymath}%
\frac{\partial \ensuremath{\mathcal{N}} }{\partial t} + \fr...
...{N}} }{h}\right) \right] = K + L \cdot \ensuremath{\mathcal{N}}\end{displaymath} (A.1)

which can be adapted to both, Eqs. (1) and (21) by proper choice of parameters u, $D_{\rm d}$, g, h, K and L.

We use a flux-conserving donor-cell scheme which is implicit in  $\ensuremath{\mathcal{N}} $. The time derivative in Eq. (A.1), written in a discretized way becomes

\begin{displaymath}%
\frac{\partial N}{\partial t} \hat{=} \frac{\ensuremath{\ma...
...N}} ^{i+1}_n - \ensuremath{\mathcal{N}} ^{i}_n}{t_{i+1}-t_{i}}
\end{displaymath} (A.2)

where i denotes time-dimension and n denotes space-dimension.

The advective part is discretized as

\begin{displaymath}%
\frac{\partial }{\partial x} \left( \ensuremath{\mathcal{N}...
... - \frac{F_{n-\frac{1}{2}}^{i+1} \cdot S_{n-\frac{1}{2}}}{V_n}
\end{displaymath} (A.3)

where $F^{i+1}_{n+\frac{1}{2}}$ and $S_{n+\frac{1}{2}}$ are the future flux and the surface between cell n and cell n+1 and Vn is the volume of cell n. The advective interface fluxes can then be written as
                          $\displaystyle %
F_{n+\frac{1}{2}}^{i+1}$ = $\displaystyle \ensuremath{\mathcal{N}} _n \cdot {\rm max}(0,u_{n+\frac{1}{2}}) + \ensuremath{\mathcal{N}} _{n+1} \cdot \min(0,u_{n+\frac{1}{2}})$ (A.4)
$\displaystyle F_{n-\frac{1}{2}}^{i+1}$ = $\displaystyle \ensuremath{\mathcal{N}} _{n-1} \cdot {\rm max}(0,u_{n-\frac{1}{2}}) + \ensuremath{\mathcal{N}} _n \cdot \min(0,u_{n-\frac{1}{2}})$ (A.5)

The diffusive interface flux becomes
                     $\displaystyle %
F_{{\rm d},n+\frac{1}{2}}^{i+1}$ = $\displaystyle D_{{\rm d},n+\frac{1}{2}} \: h_{n+\frac{1}{2}} \: \frac{\frac{g_{...
...l{N}} _{n+1} - \frac{g_{n}}{h_{n}}\ensuremath{\mathcal{N}} _{n}}{x_{n+1}-x_{n}}$ (A.6)
$\displaystyle F_{{\rm d},n-\frac{1}{2}}^{i+1}$ = $\displaystyle D_{{\rm d},n-\frac{1}{2}} \: h_{n-\frac{1}{2}} \: \frac{\frac{g_{...
...} _{n} - \frac{g_{n-1}}{h_{n-1}}\ensuremath{\mathcal{N}} _{n-1}}{x_{n}-x_{n-1}}$ (A.7)

Working out these equations and separating the values of $\ensuremath{\mathcal{N}} $ leads to

\begin{displaymath}%
\ensuremath{\mathcal{N}} _n^i = A_n \cdot \ensuremath{\math...
...^{i+1} + C_n \cdot \ensuremath{\mathcal{N}} _{n+1}^{i+1} + D_n
\end{displaymath} (A.8)

with the coefficients

                                     
  $\displaystyle A_n = -\frac{\Delta t}{V_n} \left( {\rm max}(0,u_{n-\frac{1}{2}})...
...ac{1}{2}} \cdot S_{n-\frac{1}{2}} \cdot g_{n-1}}{(x_n-x_{n-1}) h_{n-1}} \right)$  
    $\displaystyle B_n = 1- \Delta t L_n + \frac{\Delta t}{V_n} \Bigg( {\rm max}(0,u_{n+\frac{1}{2}}) \cdot S_{n+\frac{1}{2}}$  
    $\displaystyle \qquad~ - ~ \min(0,u_{n-\frac{1}{2}})\cdot S_{n-\frac{1}{2}}$  
    $\displaystyle \qquad~ +~ \frac{D_{{\rm d},n+\frac{1}{2}}\cdot h_{n+\frac{1}{2}}...
...t h_{n-\frac{1}{2}}\cdot S_{n-\frac{1}{2}} \cdot g_n}{(x_n-x_{n-1}) h_n} \Bigg)$  
    $\displaystyle C_n = \frac{\Delta t}{V_n} \left( \min(0,u_{n+\frac{1}{2}})\cdot ...
...frac{1}{2}}\cdot S_{n+\frac{1}{2}}\cdot g_{n+1}}{(x_{n+1}-x_n) h_{n+1}} \right)$  
    $\displaystyle D_n = -\Delta t \cdot K_n.$ (A.9)

Equation (A.8) 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. (A.1) is implicit only in  $\ensuremath{\mathcal{N}} $ 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  $\ensuremath{\nu_{\rm g}} $ 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 mi gives

\begin{displaymath}%
\frac{\partial }{\partial t} n_k(r,z) = \sum_{ij} M_{ijk} ~n_i(r,z)~ n_j(r,z),
\end{displaymath} (A.10)

where the dust particle number density is

\begin{displaymath}%
n_i(r,z) = \int_{m_{i-1/2}}^{m_{i+1/2}} n(m,r,z)~ \ensuremath{{\rm d}} m.
\end{displaymath} (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),

\begin{displaymath}%
n_k(r,z) = \frac{N_k(r)}{\sqrt{2\pi}h_k(r)} \cdot \exp\left(-\frac{z^2}{2h_k(r)^2} \right),
\end{displaymath} (A.12)

we can vertically integrate the coagulation/fragmentation equation.
                           $\displaystyle %
\frac{\partial }{\partial t} N_k(r)$ = $\displaystyle \int_{-\infty}^\infty \dot n_k \: {\rm d}z
= \sum_{ij} M_{ijk} \: N_i(r) \: N_j(r) \: \frac{1}{2\: \pi \: h_i \: h_j}$  
    $\displaystyle \times ~\int_{-\infty}^\infty \exp{\left(-\frac{h_i^2+h_j^2}{2h_i^2h_j^2} \cdot z^2\right)} {\rm d}z$  
  = $\displaystyle \sum_{ij} \widetilde M_{ijk} \: N_i(r) \: N_j(r),$ (A.13)

where

\begin{displaymath}%
\widetilde M_{ijk} = \frac{1}{\sqrt{2\:\pi\:\left(h_i^2+h_j^2\right)}}\: M_{ijk}.
\end{displaymath} (A.14)

Discretizing Eq. (A.13) also in radial direction and rewriting it in terms of the quantity

\begin{displaymath}%
\ensuremath{\mathcal{N}} _{kl} = N_k(r_l)\cdot r_l
\end{displaymath} (A.15)

yields

\begin{displaymath}%
\frac{\partial }{\partial t} \ensuremath{\mathcal{N}} _{kl}...
...cal{N}} _{il}\cdot \ensuremath{\mathcal{N}} _{jl} := {S}_{kl}.
\end{displaymath} (A.16)

where the vector $\vec{S}=\{S_{kl}\}_{k=1,m}$ is the source function for each of the m mass bins.

The numerical change of $\ensuremath{\mathcal{N}} _{kl}$ within a time step $\Delta t = t_i-t_{i-1}$ is $\Delta \ensuremath{\mathcal{N}} _{kl} = \ensuremath{\mathcal{N}} _{kl}^{i+1}-\ensuremath{\mathcal{N}} _{kl}^{i}$. The time-discretized version of Eq. (A.16) then becomes (omitting second order terms)

  $\displaystyle \frac{\Delta \ensuremath{\mathcal{N}} _{kl}}{\Delta t} = \sum_{ij...
...t(\ensuremath{\mathcal{N}} _{jl} + \Delta \ensuremath{\mathcal{N}} _{jl}\right)$  
  $\displaystyle \; = \sum_{ij} \frac{\widetilde M_{ijk}}{r_l} \: \left( \ensurema...
...a \ensuremath{\mathcal{N}} _{il} \Delta \ensuremath{\mathcal{N}} _{jl}}\right).$ (A.17)

Since the first term on the right hand side is the explicit source function, we can write

\begin{displaymath}%
\frac{\Delta \ensuremath{\mathcal{N}} _{kl}}{\Delta t} = S_...
...mathcal{N}} _{jl} \cdot \Delta \ensuremath{\mathcal{N}} _{il}.
\end{displaymath} (A.18)

Using the vectors $\Delta \vec{N}= \{\Delta \ensuremath{\mathcal{N}}\}_{k=1,n_m}$ and $\vec{N}= \{\ensuremath{\mathcal{N}}\}_{k=1,n_m}$, this can be rewritten in matrix notation,

\begin{displaymath}%
\left( \frac{\mathbbm{1}}{\Delta t} -\mathbf{J} \right) \Delta \vec{N} = \vec{S}.
\end{displaymath} (A.19)

Where

\begin{displaymath}%
J_{ki} = \sum_{j} \frac{\widetilde{M}_{ijk} + \widetilde{M}_{jik} }{r_l}\:\ensuremath{\mathcal{N}} _{jl}
\end{displaymath} (A.20)

denotes the Jacobian of the source function and  $\mathbbm{1}$ the unity matrix. The solution for the future values can now be derived by inverting the matrix in Eq. (A.19),

\begin{displaymath}%
\vec{N}^{i+1} = \vec{N}^{i}+\Delta \vec{N}
= \vec{N}^{i}+ ...
...\mathbbm{1}}{\Delta t} -\mathbf{J} \right)^{-1} \cdot \vec{S}.
\end{displaymath} (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

\begin{displaymath}%
\mathbf{M} \cdot \vec{N}^{i+1}= \vec{N}^i- \vec{D},
\end{displaymath} (A.22)

where M is the tri-diagonal matrix with entries A,B and C and source term $\vec{D}$ which are given by Eq. (A.9). M is now rewritten by separating off the diagonal terms and by dividing by $\Delta t$,

\begin{displaymath}%
\mathbf{M} = \Delta t\cdot \mathbf{\hat M} + \mathbbm{1}
\end{displaymath} (A.23)

which brings Eq. (A.22) in a form similar to Eq. (A.19),

\begin{displaymath}%
\left( \frac{\mathbbm{1}}{\Delta t} + \mathbf{\hat M} \righ...
...\vec{N} = -\mathbf{\hat M} \cdot \vec{N}^i - \vec{D}/\Delta t.
\end{displaymath} (A.24)

The coagulation/fragmentation is now determined by the matrix  J and the source vector $\vec{S}$ and similarly, the radial motion is determined by matrix  $\mathbf{\hat M}$ and source vector

\begin{displaymath}%
{\bf R} = -{\mathbf{\hat M} \cdot \vec{N}^i - \vec{D}}/\Delta t.
\end{displaymath} (A.25)

This allows us to combine both operators into one matrix equation,

\begin{displaymath}%
\left( \frac{\mathbbm{1}}{\Delta t} + \mathbf{\hat M} -\mathbf{J} \right) \cdot \Delta \vec{N} = \bf {R} + \vec{S}.
\end{displaymath} (A.26)

In principle, to solve for the vector  $\Delta \vec{N}$, 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 150-500 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.

\begin{figure}
\par\includegraphics[width=7cm,clip]{plots/13731fg15}
\vspace*{3mm}
\end{figure} 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 ( $\mathbf{\hat M}$), 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 0002.

Open with DEXTER

\begin{figure}
\par\includegraphics[width=8.5cm,clip]{plots/13731fg16}
\end{figure} 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, $\rho _{\rm s}$ = 1.6  ${\rm g/cm}^3$, $\ensuremath {\Sigma _{\rm g}} =18$  ${\rm g/cm}^2$, $\ensuremath{\Sigma_{\rm d}} =0.18~{\rm g/cm}^2$ and $\alpha = 10^{-3}$.

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, $\tau$, which can be written as

\begin{displaymath}%
\frac{{\rm d}m}{{\rm d}t} = \frac{m}{\tau}\cdot
\end{displaymath} (B.1)

Using ${\rm d}m={4\pi} \rho_{\rm s}a^2 {\rm d}a$ and $\tau = m/(\ensuremath{\rho_{\rm d}}\;\sigma \; \Delta u)$, we derive

\begin{displaymath}%
\frac{{\rm d}a}{{\rm d}t} = \frac{\ensuremath{\rho_{\rm d}} }{\rho_{\rm s}} \; \Delta u,
\end{displaymath} (B.2)

where $\Delta u$ is the relative velocity, $\ensuremath{\rho_{\rm d}} $ the dust density, $\rho _{\rm s}$ the solid density of the dust particles and $\sigma = 4\pi a^2$ 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

\begin{displaymath}%
\Delta u_{\rm BM} = \sqrt{\frac{16 \; \ensuremath{k_{\rm b}}\; T}{\pi \; m}},
\end{displaymath} (B.3)

for turbulent motion, the relative velocities are (see Ormel & Cuzzi 2007)

\begin{displaymath}%
\Delta u_{\rm TM} \approx \left\{ \begin{array}{ll}
\ensure...
...}{{\rm St}}} & {\rm for }~{\rm St}\gg 1.\\
\end{array}\right.
\end{displaymath} (B.4)

Integrating Eq. (B.2) gives

\begin{displaymath}%
a(t) = \left( \frac{5}{2} \frac{\rho_{\rm d}}{\rho_{\rm s} ...
...left( t - t_0\right) + a_0^{\frac{5}{2}} \right)^{\frac{2}{5}}
\end{displaymath} (B.5)

for Brownian motion and

\begin{displaymath}%
a(t) \approx \left\{\begin{array}{ll}
a_0\cdot\exp\left(\fr...
...(t - t_0) + a_0 & {\rm for }~{\rm St}> 1\\
\end{array}\right.
\end{displaymath} (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).

\begin{figure}
\par\includegraphics[width=8.5cm,clip]{plots/13731fg17}
\vspace*{6mm}
\end{figure} Figure B.1:

Test case for only radial drift, without coagulation. The solid line denotes the mass-averaged position of the radial distribution of grains for each grain size after 103 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 straight-forward since both the turbulent relative velocities and the dust scale height (Eq. (51)) are sub-divided 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 mono-disperse approximation. For grains larger than ${\rm St}=1$, 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  ${\rm St}=1$) 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

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

  \begin{figure}
\par\includegraphics[width=8.5cm,clip]{plots/13731fg01}
\end{figure} Figure 1:

Total amount of in-fallen surface density as function of radius according to the Shu-Ulrich infall model (see Sect. 2.5) assuming a centrifugal radius of 8 AU (solid line), 33 AU (dashed line), and 100 AU (dash-dotted line).

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=16cm,clip]{plots/13731fg02}
\end{figure} 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 $\alpha $ 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

  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{plots/13731fg03}
\end{figure} 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

  \begin{figure}
\par\includegraphics[width=8.5cm,clip]{plots/13731fg04}
\end{figure} 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 Fig. 5 in Hueso & Guillot (2005).

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=7.8cm,clip]{plots/13731fg05a}\includegraphics[width=7.8cm,clip]{plots/13731fg05b}
\end{figure} 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 $a_{\rm St=1} \propto \ensuremath {\Sigma _{\rm g}} $ (see Eq. (14)) this curve in fact has the same shape as  $\ensuremath {\Sigma _{\rm g}} (r)$, 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

  \begin{figure}
\par\includegraphics[width=16cm,clip]{plots/13731fg06}
\end{figure} 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 Sect. 3.4 for more details.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{plots/13731fg07}
\end{figure} 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

  \begin{figure}
\par\includegraphics[width=8.5cm,clip]{plots/13731fg08}
\end{figure} 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 Fig. 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.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8.5cm,clip]{plots/13731fg09}
\end{figure} 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 Fig. 7.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8.5cm,clip]{plots/13731fg10}
\end{figure} Figure 10:

Evolution of the dust surface density distribution of the fiducial model at 200 000 years for different drift efficiencies $E_{\rm d}$, without fragmentation. The solid line denotes the grain size $a_{{\rm St}=1}$ 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

  \begin{figure}
\par\includegraphics[width=8.5cm,clip]{plots/13731fg11}
\end{figure} 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 $E_{\rm d} = 0.75$ and fragmentation is not taken into account. It can be seen that $a_{\rm {\rm St}=1}$ (solid line) increases with radius until about 1 AU, which facilitates the break through the drift barrier.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=18cm,clip]{plots/13731fg12}
\vspace*{4mm}
\end{figure} 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\ensuremath {\times 10^{5}} $ 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 $a_{{\rm St}=1}$ 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 $\alpha $ 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-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 ( $\alpha = 10^{-3}$ and 10-2, respectively).

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8.5cm,clip]{plots/13731fg13}
\end{figure} 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 ( $\alpha = 10^{-3}$, solid line) and high ( $\alpha =10^{-2}$, dashed line) turbulence.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8.5cm,clip]{plots/13731fg14}
\end{figure} Figure 14:

Evolution of the global dust-to-gas 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 dash-dotted line shows the evolution of the disk if a low fragmentation velocity is assumed throughout the disk.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=7cm,clip]{plots/13731fg15}
\vspace*{3mm}
\end{figure} 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 ( $\mathbf{\hat M}$), 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 0002.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8.5cm,clip]{plots/13731fg16}
\end{figure} 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, $\rho _{\rm s}$ = 1.6  ${\rm g/cm}^3$, $\ensuremath {\Sigma _{\rm g}} =18$  ${\rm g/cm}^2$, $\ensuremath{\Sigma_{\rm d}} =0.18~{\rm g/cm}^2$ and $\alpha = 10^{-3}$.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8.5cm,clip]{plots/13731fg17}
\vspace*{6mm}
\end{figure} Figure B.1:

Test case for only radial drift, without coagulation. The solid line denotes the mass-averaged position of the radial distribution of grains for each grain size after 103 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

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.