Contents

Abstract

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. 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. We introduce a new model similar to Brauer et al. (2008a) 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. 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.

Abstract

Grains in circumstellar disks are believed to grow by mutual collisions and subsequent sticking due to surface forces. Many fields of research involving circumstellar disks, such as radiative transfer calculations, disk chemistry, magneto-hydrodynamic simulations and others largely depend on the resulting grain size distribution.

As detailed calculations of grain growth and fragmentation are both numerically challenging and computationally expensive, we aim to find simple recipes and analytical solutions for the grain size distribution in circumstellar disks for a scenario in which grain growth is limited by fragmentation and radial drift can be neglected.

We generalize previous analytical work on self-similar steady-state grain distributions. Numerical simulations are carried out to identify under which conditions the grain size distributions can be understood in terms of a combination of power-law distributions. A physically motivated fitting formula for grain size distributions is derived using our analytical predictions and numerical simulations. We find good agreement between analytical results and numerical solutions of the Smoluchowski equation for simple shapes of the kernel function. The results for more complicated and realistic cases can be fitted with a physically motivated “black box” recipe presented in this chapter.

Abstract

Protoplanetary disks are observed to remain dust-rich for up to several million years. Theoretical modeling, on the other hand, raises several questions. Firstly, dust coagulation occurs so rapidly, that if the small dust grains are not replenished by collisional fragmentation of dust aggregates, most disks should be observed to be dust poor, which is not the case. Secondly, if dust aggregates grow to sizes of the order of centimeters to meters, they drift so fast inwards, that they are quickly lost. We attempt to verify if collisional fragmentation of dust aggregates is effective enough to keep disks ‘dusty’ by replenishing the population of small grains and by preventing excessive radial drift. With a new and sophisticated implicitly integrated coagulation and fragmentation modeling code, we solve the combined problem of coagulation, fragmentation, turbulent mixing and radial drift and at the same time solve for the 1-D viscous gas disk evolution. We find that for a critical collision velocity of 1 m/s, as suggested by laboratory experiments, the fragmentation is so effective, that at all times the dust is in the form of relatively small particles. This means that radial drift is small and that large amounts of small dust particles remain present for a few million years, as observed. For a critical velocity of 10 m/s, we find that particles grow about two orders of magnitude larger, which leads again to significant dust loss since larger particles are more strongly affected by radial drift.

Abstract

Observations at sub-millimeter and mm wavelengths will in the near future be able to resolve the radial dependence of the mm spectral slope in circumstellar disks with a resolution of around a few AU at the distance of the closest star-forming regions. We aim to constrain physical models of grain growth and fragmentation by a large sample of (sub-)mm observations of disks around pre-main sequence stars in the Taurus-Auriga and Ophiuchus star-forming regions. State-of-the-art coagulation/fragmentation and disk-structure codes are coupled to produce steady-state grain size distributions and to predict the spectral slopes at (sub-)mm wavelengths. This work presents the first calculations predicting the mm spectral slope based on a physical model of grain growth. Our models can quite naturally reproduce the observed mm-slopes, but a simultaneous match to the observed range of flux levels can only be reached by a reduction of the dust mass by a factor of a few up to about 30 while keeping the gas mass of the disk the same. This dust reduction can either be caused by radial drift at a reduced rate or during an earlier evolutionary time (otherwise the predicted fluxes would become too low) or due to efficient conversion of dust into larger, unseen bodies.

\chapterstyle

veelo \setsecnumdepthsubsection Dissertation submitted to the Combined Faculties for the Natural Sciences and for Mathematics of the Ruperto-Carola University of Heidelberg, Germany for the degree of Doctor of Natural Sciences

Put forward by

[0.5cm] M.Sc. Tilman Birnstiel Born in Miltenberg am Main Date of oral examination 18.10.2010

\movetooddpage

The Evolution of Gas and Dust in Protoplanetary Accretion Disks Referees: Prof. Dr. Cornelis P. Dullemond Prof. Dr. Ralf S. Klessen

\movetooddpage

Zusammenfassung

Obwohl Staub nur etwa ein Prozent der Masse einer protoplanetaren Scheibe ausmacht, ist er doch ein wichtiger Bestandteil von chemischen Modellen, Modellen zur Planetenentstehung oder in der Modellierung von Strahlungstransport und Beobachtungen. Der anfängliche Wachstumsprozess von sub-m großen Staubpartikeln bis hin zu Planetesimalen, sowie der radiale Transport von Staub in den Gasscheiben um junge Sterne ist das Thema dieser Arbeit. Radiale Drift, Sedimentation, turbulenter Transport sowie Koagulation, Fragmentation und Erosion bestimmen die zeitliche Entwicklung von zirkumstellarem Staub.

Wir gehen dieses Problem von drei verschiedenen Richtungen an: analytische Berechnungen, numerische Simulationen und Vergleich zu Beobachtungen. Wir beschreiben die physikalischen und numerischen Konzepte, mit denen unser Modell die Entwicklung des Staubes über Millionen von Jahren in einer zeitabhängigen, viskosen Gasscheibe simuliert. Wir vergleichen die simulierten Staubgrößenverteilungen mit unseren analytischen Vorhersagen und leiten daraus ein einfaches Rezept zur Bestimmung von stationären Größenverteilungen ab. Mit dem vorgestellten Modell ist es uns möglich zu zeigen, dass Fragmentation erklären kann, warum zirkumstellare Scheiben für mehrere Millionen Jahre reich an Staub sein können. Des Weiteren befassen wir uns mit dem Problem, vor das uns Beobachtungen stellen, nämlich dass Staubpartikel von der Größe einiger Millimeter in großen Abständen von ihrem Zentralstern nachgewiesen wurden. Unter der Annahme, dass radiale Drift ineffektiv ist, können wir einige beobachtete spektrale Indices und Flüsse im mm-Wellenlängenbereich reproduzieren. Lichtschwächere Quellen deuten darauf hin, dass das Staub-zu-Gas Verhältnis oder die Opazitäten geringer sind als weithin angenommen. \movetooddpage

Abstract

Dust constitutes only about one percent of the mass of circumstellar disks, yet it is of crucial importance for the modeling of planet formation, disk chemistry, radiative transfer and observations. The initial growth of dust from sub-m sized grains to planetesimals and also the radial transport of dust in disks around young stars is the topic of this thesis. Circumstellar dust is subject to radial drift, vertical settling, turbulent mixing, collisional growth, fragmentation and erosion.

We approach this subject from three directions: analytical calculations, numerical simulations, and comparison to observations. We describe the physical and numerical concepts that go into a model which is able to simulate the radial and size evolution of dust in a gas disk which is viscously evolving over several million years. The resulting dust size distributions are compared to our analytical predictions and a simple recipe for obtaining steady-state dust size distributions is derived. With the numerical model at hand, we show that grain fragmentation can explain the fact that circumstellar disks are observed to be dust-rich for several million years. Finally, we investigate the challenges that observations present to the theory of grain evolution, namely that grains of millimeter sizes are observed at large distances from the star. We have found that under the assumption that radial drift is ineffective, we can reproduce some of the observed spectral indices and fluxes. Fainter objects point towards a reduced dust-to-gas ratio or lower dust opacities. \movetooddpage\movetooddpage \epigraphtextpositionflushleftright \epigraphThe light dove, cleaving in free flight the thin air, whose resistance it feels, might imagine, that her movements would be far more free and rapid in airless space. Critique of pure reason (1787), Introduction, III.
– Immanuel Kant

Contents

Chapter \thechapter Introduction

The average density in the universe today is about one atom per cubic meter (Schneider 2008). Nevertheless, we are living on an over-density of 29 orders of magnitude, known to us as the Earth. Even though the question how our world was created is almost as old as mankind itself, we are still far away from being able to answer it. Today, astronomers are able to draw a coarse picture of how stars and their planetary systems are “produced” but at the same time, many more questions have been discovered which have not been posed before.

This chapter tries to give a very brief overview over the historical aspects and the nowadays widely believed scenario of star and planet formation, which sets the scene for the following chapters of this thesis.

1 Historical perspective

The question about the origin of the Earth and our solar system has always been particularly attractive as it is one of the keys to constraining the uniqueness of our own existence. The earth was believed to be the center of the universe, until 400 years ago in January 1610, Galileo Galilei pointed his telescope to Jupiter (Drake 2003). By the discoveries of the moons of Jupiter and the phases of Venus, he did not only become the father of modern astronomy but also did he degrade the Earth to be just one planet amongst others, orbiting the Sun. This hypothesis was soon accepted by other thinkers of his time as evidence was accumulating.

Figure 1: The nebula hypothesis, first suggested by Laplace. In this picture of solar system formation, the initial cloud contracts (a) the contraction causes the forming disk to spin up due to angular momentum conservation (b). Once the centrifugal force equals the gravitational force of the star, the disk fragments locally into rings (c) and finally, planets condense out of the rings (d).

While previous theories about the world were mostly mystic or philosophical, the Age of Enlightenment and the observations of Galilei and Kepler led to more mechanical and more realistic philosophies. One of the first of such philosophies about the solar system was put forward by Descartes in 1644. He imagined the solar system and the universe in total as a system of vortices, thus implying that each star itself could host a planetary system like our own. While the idea that the solar system was born out of a large vortex was not too far from what we know today, his theory still lacked a mathematical basis. He writes (see, Gaukroger 2006):

“Finally, we see that, although these whirlpools always attempt a circular motion, they practically never describe perfect circles, but sometimes become too great in width or in length. Thus we can easily imagine that all the same things happen to the planets; and this is all we need to explain all their remaining phenomena.”

Nevertheless, things were a bit more complicated than Descartes imagined and his ideas were quickly swept away by the Newtonian/Keplerian system which seemed to explain perfectly the motions of the planets by a mathematical formalism. Yet it did not explain the origin of the solar system.

During the 18th century, two scenarios were developed, the first of which was later known as the tidal theory, first put forward by the french naturalist Georges Buffon, the other one was called the nebula hypothesis which was proposed by Pierre Simon de Laplace in 1796 (Murdin 2001). In Buffons scenario, a close encounter between the sun and a comet should have stripped off material which later condenses to form the earth and the other planets. This theory implied that the formation of the solar system was a particularly rare event which was consistent with the fact that no evidence for other planetary systems was known at that time.

For Laplace, on the other hand, the very uniform motion of the planets and their satellites proved that the system was created by a more general law of nature. Thus he proposed that the solar atmosphere once extended much further out through the whole solar system, initially being hot and luminous, and then cooled and contracted (see Fig. 1). In his view, matter would spin up due to angular momentum conservation until it condenses locally, forming rings around the sun out of which in turn the planets condense. Other defenders of this scenario include Immanuel Kant and Friedrich Wilhelm Herschel who already observed nebulae which seemed to resemble the nebulae predicted by this theory. In this scenario, planet formation is a natural, universal by-product of star formation which would mean that each star could in principle have a planetary system.

This theory became more and more accepted and few people at that time still considered theories of interaction. In 1905, T. C. Chamberlin and F. R. Moulton successfully revived Buffons one-century-old dualistic theory (see Brush 1996), arguing that the Earth could not have formed out of hot gas because its atmosphere would have evaporated. They suggested that the solar system formed through a close encounter of the sun with another star which induces tidal bulges which are then thrown off into the surrounding space. However, unlike in Buffons scenario, the planets did not “condense” out of this matter, in their view, but were assembled out of small cold particles which they called planetesimals (from ’infinitesimal planets’). But also the tidal theory had its flaws: Henry N. Russel showed that the expelled gas would be far too hot to condense into planets and that it is impossible to expel enough material to large enough radii to explain the existence of the gas giants (Murdin 2001).

During the first and the beginning of the second half of the twentieth century, many different scenarios and theories were discussed, amongst them the externally-triggered cloud collapse model of A. G. W. Cameron, the theory of magnetic breaking of Hannes Alfvén and Fred Hoyle or the Bondi-Hoyle-Lyttleton cloud capture scenario. One of the most influencing ones was by Carl Friedrich von Weizsäcker who revisited the idea of Descartes’ vortices on a more solid basis. He was also the first to realize that a turbulent disk would transport mass inside and angular momentum outside (Bodenheimer 2006).

The subject of star and planet formation experienced an ever increasing impetus since the second half of the last century. Larson (1969), Penston (1969), and Shu (1977) derived first solutions of the collapse of self gravitating spheres, while the development of computers enabled astronomers to carry out simulations of the proposed scenarios which could not be computed analytically. Today, it is generally accepted that stars form out of molecular clouds with turbulence being both triggering as well as regulating the formation of stars (see, e.g. McKee & Ostriker 2007).

If a cloud core exceeds a certain critical mass, called the Jeans mass1, it becomes gravitationally unstable and will collapse, unless other stabilizing effects, such as turbulence, rotational energy or magnetic fields are strong enough to stabilize the cloud. These effects also gain importance as the cloud collapses: rotational energy can cause binary formation and/or the the formation of a circumstellar disk (Larson 2003). Additionally, magnetic fields play a role by driving outflows and transporting angular momentum.

These facts show how diverse and complicated the formation of stars is, already on larger scales. Things will not become simpler if we focus on smaller scales now. Large parts of our current basic understanding of disk evolution and planet formation has to be contributed to a few influential papers: Lynden-Bell & Pringle (1974) and Shakura & Sunyaev (1973) pioneered the theory of accretion disks. These disks are an important aspect of star formation and at the same time the birth places of planets. Safronov (1969) and independently Goldreich & Ward (1973) envisioned a scenario of planetesimal formation which is (with modifications) still discussed today.

2 Theoretical concepts

Circumstellar disks are the birth places of planets, thus understanding the mechanism which determine their shape and structure is the key to finding out how planets can form. The evolution of the gas and the dust in circumstellar disks, the growth and fragmentation of the precursors of planet formation and some observational implications thereof are the topic of this thesis. In the following, we will discuss the most important theoretical concepts of circumstellar disks and planet formation and also some observational constrains which together provide the framework of this thesis.

2.1 Structure and evolution of circumstellar disks

Already the most basic observations of our solar system contain valuable hints about the formation of stars and planets: most of the mass (about 99.87%) of the solar system is contained in the sun, which consists almost entirely of hydrogen and helium. In contrast to this, the angular momentum of Jupiter already exceeds the angular momentum of the sun by about two orders of magnitude (e.g., Armitage 2009). The typical angular momentum of a one solar mass molecular cloud core (see Goodman et al. 1993) is another 3-4 orders of magnitude larger than the total angular momentum of the solar system (assuming a sphere of constant density in solid-body rotation). This tells us that which ever mechanism produced our solar system out of such a cloud core needs to transport angular momentum outside while most mass in transfered inwards. This is known as the angular momentum problem of star formation.

A collapse of a cloud with non-vanishing angular momentum will form a disk. This can be understood if one imagines each parcel of gas to orbit around the core center maintaining its angular momentum. If the total angular momentum of the previously mentioned sphere of constant density is about g cm s (assuming a density of cm and the ratio of rotational over gravitational energy ), then the parcel of gas with an average specific angular momentum will hit the equator at a radius of about  AU (Ulrich 1976). At the mid-plane, it will collide with a symmetrically moving parcel of gas. This will produce a so called accretion shock where the kinetic energy is dissipated and thus a disk is formed. In reality, the physics of a cloud core collapse is much more complicated because the rotational support can produce multiple stars/disks and turbulence and magnetic fields also play an important role during the collapse. The formation of a disk is however a well established outcome, also of more realistic simulations of cloud collapse (see, Inutsuka et al. 2010, for state-of-the-art simulations). The disks around new-born stars are usually called circumstellar disks, accretion disks (implying an inward flux of mass) or protoplanetary disks (implying that planets will form inside of them). Due to the geometry of a disk, the vertical sound-crossing time scale is much shorter than the radial one, thus it is usually a good approximation to treat the vertical structure to be in an equilibrium, while the radial evolution is happening on larger time scales.

A first order vertical structure of the gas disk can be derived by assuming the disk to be isothermal (vertical temperature profile const.) and thin (i.e., the height above the mid-plane is small compared to the distance to the star ). Then, the vertical component of the stellar gravity needs to be balanced by a vertical pressure gradient,

(1)

where is the isothermal sound speed, the constant of gravity, the gas density, and the mass of the central star. From this we derive by integration a Gaussian density profile

(2)

where we denoted the mid-plane density as and defined the pressure scale height to be the ratio of the sound speed over the Kepler frequency .

Having a first expression for the vertical structure, we can now focus on the radial structure and its evolution. Depending on the desired resolution/complexity, this can be done in 1, 2 or 3 dimensions. However, in this thesis, we are interested in (parameter studies of) the long-time evolution of these disks and multidimensional calculations are computationally too expensive to simulate several million years of disk evolution (especially if grain growth is considered as we will discuss later on). Therefore – and also to get an analytical feeling of the disk evolution – we will now consider the 1+1D approach, where an axisymmetric disk is assumed to be in a vertical equilibrium as discussed above. We are thus only interested in the time evolution of the surface density, which is the vertically integrated column density, defined as

(3)

If there is viscosity acting within the disk (we will cover this question further below), then two neighboring sheets of gas will interact through friction due to the differential rotation. The gas is approximately orbiting in Keplerian rotation. Therefore, an inner sheet of gas is rotating faster than an outer one and viscous friction thus accelerates the outer sheet while the inner sheet is decelerated. This mechanism obviously transports angular momentum outside and at the same time allows matter to flow inwards.

The mathematical description of this process is called the theory of accretion disks and was first derived by Lynden-Bell & Pringle (1974). They derived the radial velocity of the gas due to mass and angular momentum conservation in the presence of a viscosity to be

(4)

The evolution of the surface density is then given by the vertically integrated mass conservation equation, and can be written as

(5)

(for a derivation, see also Pringle 1981). This equation can now be solved numerically for a given viscosity and initial condition (see Chapter id1 and id1), but it is instructive to discuss some general properties of this equation analytically, which we will do in the following.

Substituting and in Eq. 5, we can derive a diffusion equation

(6)

with diffusion constant (assuming a constant ). The diffusion constant defines the time on which diffusion acts to smooth out gradients in the concentration (in this case: ). A viscous disk with size will therefore evolve on the viscous timescale given by

(7)

From the observed live time of circumstellar disks, which is about 3 million years (see Section 3), we can now get constraints on the viscosity of these disks. It turns out that the viscous time scale of a 100 AU disk (with typical values for the temperature and surface density) is of the order of years if we consider only molecular viscosity. This is much longer than the age of the universe. To explain the observed disk lifetimes, another mechanism of angular momentum transport is needed and the anomalous viscosity provided by turbulence is typically believed to be solution (although gravitational instability and disk winds may also play a role).

The source of the turbulence is still unknown. Keplerian disks are linearly stable to axisymmetric perturbations (e.g., Safronov 1969), but a number of other possible sources have been discussed, such as the convective instability (today discarded due to the fact that it would transport angular momentum inside instead of outside, see Stone & Balbus 1996), the gravitational instability (Toomre 1964; Gammie 2001), or the baroclinic instability (Klahr & Bodenheimer 2003; Lesur & Papaloizou 2010) to name a few of them. However none of these has been as successful as the magnetorotational instability (MRI), discovered by Balbus & Hawley (1991), because it draws its energy directly from the Keplerian shear flow (Balbus & Hawley 1998). Therefore it has been shown to be a very reliable and effective mechanism of angular momentum transport which is effective enough to explain the observed accretion rates. Still, it requires two conditions in order to operate in circumstellar disks: outwardly decreasing differential rotation and an ionization fraction of at least (see Balbus & Hawley 1998).

While the former condition is always the case in (nearly) Keplerian disks, the latter condition may be questionable. A certain level of ionization is provided by cosmic rays, high energy irradiation of the central star, and by decay of radioactive nuclei within the disk. However dense enough regions and/or the large scale magnetic field line structure may shield some parts of the disk (e.g., Gammie 1996; Sano et al. 2000; Turner & Drake 2009). In this case, there will be a zone where MRI shuts off (called dead zone) surrounded by an MRI-active layer. Since the vertically averaged angular momentum transport in layered (i.e., partly-dead) regions is not as effective as in active regions, material will accumulate in the dead parts of the disk (Gammie 1996). These layered disks can become unstable and produce outbursts (thought to explain the FU Orionis objects, e.g. Armitage et al. 2001) but stationary solutions are also possible (Terquem 2008).

Coming back to the evolution of the disks surface density evolution, we note that there are possible mechanisms providing the necessary anomalous viscosity, but it is still unknown which mechanism operates and how effective it really transports angular momentum. It has therefore been proven very helpful to parametrize our ignorance in the following way, proposed by Shakura & Sunyaev (1973): turbulent viscosity is provided by eddies with a certain size and velocity . The largest reasonable values in a circumstellar disk are the scale height and the sound speed , respectively. The viscosity can therefore be written as a fraction of a product of these quantities,

(8)

where the dimensionless parameter specifies the amount of turbulence in the disk and is therefore usually bound between . Typical values of from simulations of MRI-active disks are in the range of to a few times (e.g., Johansen & Klahr 2005; Dzyurkevich et al. 2010). For a given value of , we can now calculate the viscosity if the temperature profile is known. Realistic models need to take into account the effects of external irradiation, viscous heating, optical depth, and shadowing. Thus multidimensional radiative transfer is in principle necessary to derive a realistic temperature structure, however we are only interested in some scaling relations of the surface density in this introductory chapter and will therefore assume a power-law profile for the temperature

(9)

and also for the surface density

(10)

Postulating a steady-state of the surface density () and using the power-laws for temperature and surface density in Eq. 5, we derive

(11)

This means that for a given temperature profile, a disk will evolve towards an equilibrium power-law with index . The time scale of this evolution is given by Eq. 7. From this result, we can also derive the accretion rate of the disk: the product rule applied to Eq. 4 gives

(12)

Inserting this equation into the definition of the mass accretion rate,

(13)

results in two contributions to the mass flux,

(14)

We can see now that this mass accretion rate is constant in a steady state disk since the product of is constant with radius and that the mass in a non-stationary disk will be rearranged to approach a steady state. Eq. 5 can also be solved analytically in a less hand-waving way, which can be found in Lynden-Bell & Pringle (1974). Their self-similar solution of the disk evolution shows similar features. In general, a viscous disk (with power-laws for temperature and viscosity) evolves in the following way: the disk profile in the inner parts of the disk (defined by a characteristic radius, see Hartmann et al. 1998, for example) follows the same power-law derived above, while in the outer parts, the surface density decreases exponentially with radius. The mass accretion rate causes the mass of the disk to gradually decrease. In order to achieve this, a small fraction of the mass needs to “absorb” the angular momentum. The outer parts of the disk therefore spreads while the inner parts accrete onto the star. The final outcome of this evolution is the state of lowest energy: almost all the mass is in the central object while almost all angular momentum is in a small fraction of the mass at a large radius. This arrangement is similar to what is observed in our solar system.

2.2 The dynamics of dust

So far we have only talked about the gas disk evolution. The canonical value of the dust-to-gas ratio in the interstellar medium used in the literature is only 1%, thus one might be tempted to neglect it. However, not only is this small fraction of the total mass the material out of which terrestrial planets, asteroids and the cores of giant planets are made, it is also an important ingredient to many aspects of circumstellar disks. Almost all the opacity of the disk is provided by grains. They do therefore determine most observable properties of disk. Additionally, they are important for chemical surface reaction (e.g., formation of water) or for photoelectric heating which strongly influences the temperature in the upper layers of the disk (and also the FUV photoevaporation, see Gorti et al. 2009, for example) and also for the ionization fraction, which detemins the coupling of the disk to the magnetic field and thus its MRI activity (e.g., Sano et al. 2000; Turner & Drake 2009). All these facts tell us how important the radial and also the size distribution of dust in circumstellar disks is. In this subsection we will introduce the most important aspects wich determine the evolution of circumstellar dust.

As noted above, the dust mass in circumstellar disk is believed to be (initially) only about 1% of the gas mass since the disk material stems from the interstellar medium. Therefore, the gas is the dynamically dominating phase. This can only change if some transport mechanism effectively accumulates the dust relative to the gas until the mass fraction approaches unity. We will therefore begin by discussing the transport of dust in both the vertical and the radial direction. Analyzing the interplay between gas and dust motion is beyond the scope of this chapter (for details, see Chiang & Youdin 2010, and references therein). We will rather discuss the steady-state solutions of the dust velocity without considering back reactions to the gas.

While the gas reacts to pressure forces, the dust only feels the gravitational forces acting on it and the drag force by which it is coupled to the gas motion. The drag force is given by

(15)

where is the mass of the particle, is the particle velocity relative to the gas and is the stopping time, which will be discussed in more detail in Chapter id1. This coupling to the gas causes dust particles to move towards the region of highest pressure. If we consider a particle at a height above the mid-plane, then its vertical motion would be an oscillation around the mid-plane if it were not coupled to the gas. The drag force damps this motion and thus the particle will settle towards the mid-plane. Equating the vertical component of the gravitational force and the drag force (Eq. 15) yields the terminal settling velocity

(16)

It turns out to be very practical to define a dimensionless stopping-time of the particle called the Stokes number St. It describes the coupling of the particle to the gas motion and is in our context defined as the ratio of the stopping time and the orbital time,

(17)

This way, particles with different structure or composition but the same St will behave dynamically entirely the same.

As noted before, the gas disk is most probably turbulent. The distribution of dust is thus mixed along with the turbulent gas motion, again depending on the coupling to the gas. The vertical distribution of dust is now determined by the two counteracting mechanisms of turbulent mixing and vertical settling. We can derive a steady-state distribution from the diffusion-advection equation of the dust concentration (e.g., Dubrulle et al. 1995)

(18)

Postulating a steady-state (), using Eq. 16 and integration results in a Gaussian distribution of the dust concentration with a scale height

(19)

This leads to the question what the diffusion constant is. Youdin & Lithwick (2007) showed that the dust diffusion coefficient is in good approximation given by

(20)

It is then typically assumed that the gas diffusivity is equal to the viscosity . Johansen & Klahr (2005) found that in MRI turbulence, radial diffusion is as strong as turbulent viscosity, while the vertical diffusivity is slightly weaker.

The dust-gas coupling also influences the radial distribution of the dust, again, by turbulent mixing and “radial settling” towards the star but also by co-moving it along with the radial gas flux. These topics will extensively be discussed in Chapter id1. Especially the radial drift has important implications for planet formation, as we will see in the following: larger bodies (i.e., ) are not significantly affected by the gas drag (because of their small surface to mass ratio) and thus move on Keplerian orbits. In contrast to this, very small particles () are so tightly coupled to the gas that they are basically frozen-in to the gas. In between these two limits, there is a critical regime around a Stokes number of unity: these particles are marginally coupled, thus they are not frozen-in to the gas but are orbiting at Keplerian velocity. However unlike the much larger bodies, they are still significantly influenced by the gas drag. The fact that the gas is orbiting at sub-Keplerian velocity due to its pressure support induces a relative velocity between the dust particles and the gas. The dust particles loose angular momentum because of this head wind and consequently spiral inwards. The resulting radial velocity can be written as

(21)

as derived by Weidenschilling (1977) and Nakagawa et al. (1986) (here, is the gas pressure). For typical disk models, the time scale of this drift can be very short, in the order of only a few local orbits (e.g., Brauer et al. 2007). This drift poses a fundamental problem for all theories which assume hierarchical growth from small dust particles to planetesimals and will be discussed in further detail in Chapter id1. The basic ideas of grain growth and planet(-esimal) formation will be introduced in the next section.

2.3 From dust to planets

The growth from sub-m sized monomers to larger aggregates is inevitably the first stage of the formation of planets. Small dust aggregates are sticky because of their large surface to mass ratio and inter-molecular binding forces (e.g. van der Waals forces or hydrogen bonds). This fact is proven by both experimental (e.g., Heim et al. 1999; Poppe et al. 2000) and theoretical work (e.g., Dominik & Tielens 1997). In order to assemble larger aggregates, these sub-m sized dust particles need to collide. There are several mechanisms which could provide relative motion of the dust particles: Brownian motion, the random motion caused thy the thermal motion of the gas, leads to random velocities between small particles. The smaller the particle, the more it is influenced by Brownian motion. Radial, azimuthal, and vertical motion of the dust (as discussed in the previous section) is also size dependent, thus, particles with different sizes will have different velocities which leads to relative motion between grains of different sizes. Different sized dust grains couple to different sized eddies of the turbulent motion of the gas. The relative particle velocities induced by turbulence have been derived by Ormel & Cuzzi (2007).

The relative velocities together with the efficient sticking of small dust aggregates leads to a fast growth of dust; sizes of around centimeters can be reached within less than 1000 years at a few AU, as was shown by many authors (Weidenschilling 1980, 1995; Nakagawa et al. 1981; Tanaka et al. 2005; Dullemond & Dominik 2005; Brauer et al. 2008a, to name a few). However, two problems arise: particles of decimeter sizes at a few AU (or millimeter sizes at around 100 AU) have a Stokes number of unity and, once they are formed, are quickly lost towards the star. Additionally, the relative velocities due to drift and turbulence also reach their maximum at these sizes and are large enough to fragment or erode the particles (e.g., Brauer et al. 2008a; Blum & Wurm 2008). Additionally, it is questionable if the efficient sticking assumption which is often used is still valid in this size regime (Youdin 2010). Aggregates of mm size do not seem to stick at any velocity, but rather bounce or fragment (Blum & Wurm 2008). This obstacle might be overcome by the collisions of aggregates of very different sizes (e.g., Wurm et al. 2005; Teiser & Wurm 2009).

The gravitational instability of the dust is an attractive scenario which basically “jumps” over this problematic size regime. In the original idea, developed by Safronov (1969) and independently by Goldreich & Ward (1973), growth and vertical settling of dust particles forms a dust layer at the mid-plane of the disk which becomes unstable to its own gravity and thus collapses to form planetesimals. Weidenschilling (1980) showed that this scenario does not hold since the thin dust layer is unstable to the Kelvin-Helmholtz instability because of the shear at the boundary between the dust and the gas. It turns out that mid-plane turbulence is both a curse and a cure (Youdin 2010): it may on average stir up concentrations of solids, but it also does create random over-densities. The streaming instability (Youdin & Goodman 2005) is particularly effective in creating dust over-densities, it can even create gravitationally bound clumps, as shown by Johansen et al. (2007). However, this scenario still requires particles of around decimeter sizes, underlining the importance of the initial “sticking” growth of dust.

Leaving these problems aside, a new mode of growth is reached, once particles have grown to sizes where gravity dominates the accretion of solid bodies (however, gas drag can still influence the accretion process, see Ormel & Klahr 2010). The gravitational attraction effectively increases the cross section (Lissauer 1993)

(22)

where is the relative velocity of the planetesimals and is the escape velocity of the two colliding bodies. The accretion rate in a system where gravitational focusing is important () scales superlinearly with mass (Ormel et al. 2010a). Under these conditions, it can be shown (e.g. Armitage 2009) that the ratio of the masses of two bodies growth with time. This is called runaway growth and will continue until the runaway bodies have reached sizes of a few 100 km (Wetherill & Stewart 1989; Thommes et al. 2003; Ormel et al. 2010b). At this point, the largest body gravitationally stirs up the distribution of smaller bodies and the accretion rate then becomes again sub-linearly with mass (e.g. Ida & Makino 1993; Kokubo & Ida 2000; Thommes et al. 2003). Ultimately, the growth of the largest body (called the oligarch) is stalled once it has accreted all material within its gravitational reach, the so-called feeding zone. The resulting protoplanets typically have about the mass of Mars in the terrestrial zone (Thommes & Duncan 2006).

Unless the disk mass is scaled up by a factor of 10, this scenario fails to explain the formation of Earth and Venus and also the formation of the cores of gas giant planets in situ: according to Pollack et al. (1996), a core of about 10  is needed to start runaway gas accretion and thus produce gas giants. Possible solutions to this include the migration of small bodies by aerodynamic drag as well as the migration of the planets and cores. Scattering between the cores is also an attractive scenario as it could possibly explain the formation of Uranus and Neptune which are too massive to be explained by core accretion at their current location in the disk (Thommes et al. 2003).

An entirely different approach to planet formation is the disk instability scenario: while the disk is typically gravitationally stable in the inner parts, this might be debatable in the outer parts where the disk is colder and the orbital frequency is lower. This can be seen from Toomre’s stability criterion (Toomre 1964) which states that disk becomes gravitationally unstable if the stability parameter

(23)

becomes smaller than unity. Whether the disk really fragments can however not be seen just from this parameter. Before collapse sets in (i.e., ), the disk develops spiral arms which tend to stabilize the disk by rearranging the surface density (Gammie 2001). Additionally, the disk must be able to cool quickly enough so that a contracting clump is not sheared apart by the Keplerian rotation. Several authors (e.g., Boss 1997, 2000, 2007; Rice et al. 2005; Meru & Bate 2010) have produced very different results which depend highly on the efficiency of cooling in the disk (and thus also on the numerical scheme). Today it seems unlikely that disk instability alone can explain giant planet formation since the simulations typically produce planets at the high end of the mass scale (Rafikov 2005), however some of the planets detected by direct imaging (see Sect. 3.3) are candidates for this formation scenario.

3 Observational constraints

After having introduced the basic concepts theoretically, it is time to provide some proof from the observational side which led to our current understanding (and also the problems) of disk evolution and planet formation. In this section, we will discuss the most important constraints, derived from observations, ordered by size: from observations of gas and tiny dust grains, to the meteorites in our own system, and finally to the observed population of planets in and beyond our solar system.

Figure 2: Origin of the spectral energy distribution: at larger wavelength, the disk emission outweighs the stellar spectrum. Generally, the shorter wavelength emission comes from hotter, inner regions of the disk. Large wavelength (sub-millimeter) probe the dust mass in the outer, optically thin regions. Figure adapted from Dullemond et al. (2007).

3.1 Disk observations

Circumstellar disks can be detected by several techniques, each tracing specific features, temperatures and/or regions of the disk. Combining this knowledge enables us not only to proof the existence of disks but also to probe the physics which is happening inside of them.

Most of the disk mass is expected to reside in gas (because disks are made out of interstellar material), however this component is the hardest to detect as it is often frozen out, especially in the outer, cold parts of the disk and in the mid-plane.

Another gas feature is a strong H excess, which is is believed to be due to an accretion flow onto the stellar surface (see Bouvier et al. 2007). Traditionally, T Tauri stars are divided into classical T Tauri stars (CTTS) and weak-lined T Tauri stars (WTTS), depending on the strength of the emission lines. The accreting material has to come from circumstellar material and thus CTTS are thought to be younger objects which still have an accretion disk around them. Indeed, the accretion rate (typically of the order of solar masses per year) and the fraction of young stars which have detectable disks (see below) seem to correlate (in the same way) with the age of the system (e.g., Calvet et al. 2000; Haisch et al. 2001; Sicilia-Aguilar et al. 2006a, b), indicating a typical lifetime of a gaseous circumstellar disk to be around 3 million years. Only few gaseous disks have been detected around stars older than 6 million years but dust can continue to exist in debris disks, which are probably the remnants of planet(-esimal) formation, similar to the asteroid belts in our own solar system.

The accretion signature and the observed lifetimes of disks show that some mechanism must efficiently transport angular momentum outwards and mass inwards. In general, the transition from accreting to non-accreting seems to be very rapid, in the order of years (Simon & Prato 1995; Wolk & Walter 1996; Andrews & Williams 2005). This is in contradiction with simple viscous evolution which predicts a gradually fading disk (e.g., Hueso & Guillot 2005). Photoevaporative dissipation of the disk (Hollenbach et al. 1994; Clarke et al. 2001; Alexander et al. 2006) or planets blocking the accretion of the disk (e.g., Calvet et al. 2002) are some of the frequently suggested scenarios. However at every age some sources do not show any signs of accretion. The reason for this is yet unknown.

Much easier to detect than the gas is the dust which re-radiates the stellar irradiation in the infrared (IR). This emission occurs as a long-wavelength excess above the stellar photospheric flux, where different wavelength correspond to different temperatures and thus different regions in the disk. The disk geometry and the according spectral energy distribution (SED) is schematically shown in Fig. 2. In fact this IR excess was the first observational evidence for disks around young stars (e.g., Strom et al. 1989; Beckwith et al. 1990).

The emission of warm dust in the inner disk can be used as a tracer for grain growth at smaller sizes: the silicate 10 m emission feature is prominent for small dust sizes (around sub-m) but disappears as the grain size is increased beyond a few m. van Boekel et al. (2003) found evidence for grain growth at these aggregate sizes as one would expect from both theory and laboratory experiments (Blum & Wurm 2000; Blum et al. 2000).

While the near- and mid-IR traces mostly the warm/hot inner regions and the surface layers of the optically thick dust disk, observations at millimeter wavelength probe the cold outer regions of the disk. Radiation at these wavelength comes from the far end of the disk SED, which is close to the Rayleigh-Jeans limit. The disk SED can be approximated by (e.g., Beckwith et al. 1990; Andrews & Williams 2007)

(24)

if the disk is optically thin and the emission is from an isothermal region. Here denotes the frequency, the opacity (per gram of dust), the Planck function at a characteristic temperature , the dust mass, and the distance to the source. Eq. 24 shows that under these conditions, the SED is proportional to the product of dust mass and opacity. For an assumed opacity law, the dust mass can be determined. The total mass of the disk can then roughly be estimated of one assumes the dust-to-gas ratio to be the same as in the interstellar case. Disk masses were found to be ranging from 0.005 to 0.14 solar masses spread out over sizes of between 20 to several hundred AU (e.g. Andrews et al. 2009). The resolution of observations at these wavelength has reached a state where the emission can be radially resolved (e.g. Isella et al. 2010). The Atacama Large Millimeter Array (ALMA) will be able to resolve the emission down to a few AU for the most nearby disks.

In the framework of viscous disks and -turbulence (see Eqs. 5 and 8), spatially resolved observations can be used to constrain important parameters of the theoretical disk models. Results indicate for example that the turbulence parameter should be in the range of –Ð (Andrews et al. 2009), which compares well to values predicted by MRI simulations.

Additionally, the mm-wavelength regime is also an important probe of grain growth: Eq. 24 shows that the spectral dependency of the SED is given by the dust opacity and the Planck function. If we know the mid-plane temperature (or are sufficiently deep in the Rayleigh-Jeans regime), measuring the slope of the SED (by observations at two or more different wavelengths) provides a direct measure of the spectral index of the dust opacity law . depends on the dust size distribution (and the opacity model). For a power-law dust size distribution,

(25)

we can now compare which parameters and are possible for a certain measured value of (see Fig. 3). For example, a value of unity can only be reached by a distribution which has a maximum grain size of around a few millimeters.

Several authors (e.g. Testi et al. 2003; Wilner et al. 2005; Rodmann et al. 2006; Ricci et al. 2010) detected emission which indicates mm- or cm-sized dust particles at several tens of AU. On the other hand, Brauer et al. (2007) has shown that (for reasonable disk parameters) this distribution of grains should drift towards the star on a short time scale. A solution to this problem is yet unknown, but particle clumping in streaming instabilities or spiral density waves could significantly decrease the drift time scale.

Figure 3: Dust opacity (upper panel) and the resulting opacity index for different dust size distributions. The parameters of the dust size distribution are the maximum grain sizes and the (negative) power-law index (see Eq. 25).

3.2 Meteorites

Additional constrains are coming from a quite different direction, namely meteorites: they provide ample information with a precision unmatched by any astronomical observation about the formation time scales of solar system bodies. The largest class of meteorites, the chondrites are extraordinarily diverse, but they are undifferentiated which means that their parent body has not been molten which would have differentiated it into different phases. While the oldest known rocks on earth are “only” about million years old (O’neil et al. 2008) (some zircon minerals are as old as million years, see Wilde et al. 2001), meteorites can be significantly older: the age of calcium-aluminium rich inclusions (CAIs), which can be found in chondrites, was determined to be 4568.3 million years, with an uncertainty of only 700 000 years (Burkhardt et al. 2008). This way, meteorites are the “fossils of the planet formation process”, conserving hints about the structure and the chemical composition of the early solar system.

Most chondrites consist of three main components: Firstly, the CAIs, up to cm-sized particles, which are considered to be the first condensates in the solar system. Secondly, chondrules, small, mostly spherical silicate crystals with diameters in the range of 50 m to several mm (Scott 2007) comprise more than 70% by volume of the ordinary chondrites (Youdin 2010). Thirdly, the space between CAIs and chondrules is filled with the matrix, which consists of sub-m sized grains of similar minerals. A cut through a fragment of the famous Allende meteorite is shown in Fig. 4.

The decay of radioactive nuclei can be used to infer both the ages and the formation time-scales of these components (for a review, see Russell et al. 2006). Results indicate that CAIs are the oldest components and were formed during a time interval of about 0.25 million years. Chondrules were formed during an extended period of around 2 million years, starting 1–2 million years after the formation of CAIs (Kleine et al. 2009). Today, it is believed that chondrules also formed earlier than this, but these early chondrules have been incorporated into parent bodies which differentiated due to the heat from Al-decay. Dating of iron meteorites indicate that they have formed within million years after CAI formation, they are therefore often thought to be core material of the early parent bodies which were differentiated by the Al-heating.

Despite the ample laboratory data (or perhaps due to it) there is no satisfying theory which explains the formation of these bodies. Still two major points that have been learned from meteoritic data should be mentioned: firstly, apart from some very volatile materials, the composition of chondritic meteorites is strikingly similar to the composition of the sun. This tells us that meteorites (and thus the planets) formed from the same material as the sun itself. Secondly, the constraints on the lifetime of the solar nebula are in very good agreement with the lifetimes of circumstellar disks.

Figure 4: A fragment of the Allende meteorite.

3.3 Exoplanets

The most important observational constraint is that planets do exist, not only in our own solar system. The basic concepts behind the most successful planet-searching techniques are rather simple: the radial velocity method measures the “wobble” of the star induced by the common motion around the center of mass, while the transit method uses the dip in the light curve to proof the presence of an extrasolar planet. Together, these two methods are responsible for more than 93% of all exoplanet detections. Still, it is the high precision of the measurements that is needed, which makes detections of exoplanets such a challenging task. Therefore it has taken 40 years after the first proposal of the radial velocity method by Struve (1952) that an extrasolar planetary companion (around a neutron star) could be detected by Wolszczan & Frail (1992). Detections of planets around solar-type stars followed a few years later (Mayor & Queloz 1995).

Today, the number of extrasolar planets is ever increasing (unlike the number of planets in the solar system, which has decreased), proving that planet formation is indeed a common by-product of star formation. Extrapolation of recent discoveries suggest that about 20% of all stars have planets orbiting within 20 AU (Cumming et al. 2008). Population synthesis models predict an even higher rate: about 30-40% of all FGK stars should host planets (Mordasini et al. 2009).

Fig. 5 shows the mass of the currently 2 known planets as function of their semimajor axis. The known population of exoplanets suggests that stars with higher metallicity are more likely to host planets (Fischer & Valenti 2005), which is in general in agreement with the core accretion scenario.

Recently, the first direct imaging detections of extrasolar planets have been reported (e.g., Kalas et al. 2008; Marois et al. 2008), which have typically several Jupiter masses and are located at large distances from the star (hundreds of AU). Some of these are thought to have formed by disk instability (Meru & Bate 2010).

Figure 5: The mass of the currently known planets and exoplanets as function of their semimajor axis. The solid black circles (•) denote transiting planets, the empty circles (◦) denote radial velocity detections, and all other detections (pulsar timing, microlensing, and direct imaging) are displayed as black dots (). The red symbols represent the solar system planets.

4 Summary

Summarizing all of this, problems with planet formation seem to come in all sizes: smaller particles, as discussed above, face the problem that growth can be halted by charging of the grains (see Okuzumi 2009), by fragmentation and erosion (see Chapter id1), or by bouncing (see Zsom et al. 2010). If larger grains are formed, they seem to be lost quickly to the star due to the radial drift. And even if one can explain the formation of planetary cores, these bodies still face the problem of migration due to interaction with the disk. Additionally, planet formation needs to happen relatively quickly (i.e., a few million years), otherwise the gas disk, which is needed for the formation of gas giant planets, will have disappeared. Thus, we can see that we are able to understand the formation of the solar system and the planets far better than Descartes, however we are still far away from a final, consistent answer.

5 About this thesis

This thesis focuses on the very first stages of planet formation, namely the growth of sub-m particles to planetesimals along with the evolution of the radial distribution of these particles. Dust is an important ingredient in circumstellar disks, not only because terrestrial planets and the cores of giant planets are made out of it, but also it is the main source of opacity in disks (thus, important for interpreting observations), it influences the coupling of the disk to magnetic fields by influencing the ionization fraction (important for the viscous evolution of the whole disk) and it is an important catalyst in disk chemistry, especially for the formation of water. Therefore, we build a numerical code which is able to trace the radial and size evolution of dust suspended in a gaseous disk which itself is evolving due to viscous transport of angular momentum.

Chapter id1

describes the physics of growth and radial transport of dust as well as the physics of the gas disk which dominates the dynamics of the dust. This is the first work which takes into account the disk build-up phase and subsequent gas disk evolution in calculations of the dust evolution. We investigate how the dust-to-gas ratio (which is an important parameter for both observations of disks and the theory of planetesimal formation) evolves over time if coagulation and fragmentation are taken into account. We also revisit the barriers to planet formation by sticking collisions, namely the radial drift and the fragmentation barrier. We show that radial variations in the tensile strength of particles (caused by compositional changes due to the presence of ices) can cause a strong enhancement of dust inside the snow line. The fragmentation velocity turns out to be one of the most important parameters determining the particle shape as well as the dust disk shape and lifetime.

The outward spreading of the gas disk is able to carry dust along. Whether and to which extent this effect works strongly depends upon the maximum particle size and thus also on the critical fragmentation velocity of the gains. There is always a radius where even the smallest grains decouple from the spreading gas disk. Future observations (e.g., with ALMA) might be able to detect the different outer radii of the gas and the dust disks.

Chapter id1

discusses the numerical implementation of the model. The radial transport of the gas and the evolution (in size and radial position) of dust are solved using implicit integration schemes. We also show test cases which validate the results of our numerical simulations by comparison to analytic estimates.

Chapter id1

focuses on the grain distribution in the scenario of an equilibrium between coagulation and fragmentation. We generalize previous analytical work on the slope of steady-state grain size distribution by taking into account both fragmentation and coagulation at the same time. These results not only validate the distributions we derived in the previous chapters but they also give us the opportunity to understand the influence of parameters like the temperature or the gas surface density on the size distribution of grains in circumstellar disks. We also derive a fitting recipe which can readily be used for further modeling such as grain surface chemistry in disks, dust-dependent MRI simulations or radiative transfer modeling.

Chapter id1

investigates the the fact that disks are generally observed to be dusty for several million years and how this can be understood if grain fragmentation is considered. We show that grain fragmentation at impact velocities similar to what is found in experimental studies causes grains to be continuously fragmented to small enough sizes to avoid the radial drift regime and to keep the optical depth as high as it is observed.

Chapter id1

is an application of our previously derived steady state dust size distributions to the field of millimeter observations. A grid of different input parameters is used to calculate the dust distribution in steady-state. From this, we can predict the spectral indices and flux levels at millimeter wavelength and compare this to an existing sample of observations. We show that typical values can explain the low spectral indices and also the flux levels of a subset of the observations. We propose several possibilities which further reduce the fluxes while keeping the spectral indices low enough.

Chapter id1

summarizes our findings and discusses future prospects of this thesis which are planned or already in progress.

Chapter \thechapter Evolution of gas and dust in protoplanetary disks

\chapterprecishere

From Birnstiel, Dullemond, & Brauer 2010a, A&A, 513, 79

6 Introduction

The question of how planets form is one of the key questions in modern astronomy today. While it has been studied for centuries, the problem is still far from being solved. The agglomeration of small dust particles to larger ones is believed to be at least the first stage of planet formation. Both laboratory experiments (Blum et al. 2000) and observations of the 10 m spectral region (Bouwman et al. 2001; van Boekel et al. 2003) conclude that grain growth must take place in circumstellar disks. The growth from sub-micron size particles to many thousand kilometer sized planets covers 13 orders of magnitude in spatial scale and over 40 orders of magnitude in mass. To assemble a single 1 km diameter planetesimal requires the agglomeration of about dust particles. These dynamic ranges are so phenomenal that one has to resort to special methods to study the growth of particles though aggregation in the context of planet(esimal) formation.

A commonly used method for this purpose makes use of particle size distribution functions. The time dependent evolution of these particle size distribution functions has been studied by Weidenschilling (1980), Nakagawa et al. (1981), Dullemond & Dominik (2005), Brauer et al. (2008a) (hereafter B08a) and others. It was concluded that dust growth by coagulation can be very quick initially (in the order of thousand years to grow to centimeter sized aggregates at 1 AU in the solar nebula), but it stalls around decimeter to meter size due to what is known as the “meter size barrier”: a size range within which particles achieve large enough velocities to undergo destructive collisions and fast radial inward drift toward the central star (Weidenschilling 1977; Nakagawa et al. 1986).

While the existence of this meter size barrier (ranging in fact from a couple of centimeters to tens of meters at 1 AU) has been known for over 30 years, a thorough study of this barrier, including all known mechanisms that induce motions of dust grains in protoplanetary disks, and at all regions in the disk, for various conditions in the disk and for different properties of the dust (such as sticking efficiency and critical fragmentation velocity), has been only undertaken recently (see 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; Hartmann et al. 1998; Hueso & Guillot 2005), which may influence the processes of dust coagulation and fragmentation. It is the goal of this work 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 chapter 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 chapter is organized as follows: Section 7 will describe all components of the model including the radial evolution of gas and dust, as well as the temperature and vertical structure of the disk and the physics of grain growth and fragmentation. In Section 8 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 9 summarizes our findings. A detailed description of the numerical method along with results for selected test cases can be found in Chapter id1.

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

7.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 thesis 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,

(26)

where the gas radial velocity is given by (see Lynden-Bell & Pringle 1974)

(27)

is the gas surface density, the radius along the disk mid-plane and the gas viscosity. The source term on the right hand side of Eq. 26, denoted by can be either infall of material onto the disk or outflow.

The molecular viscosity of the gas is too small to account for relevant accretion onto the star, the time scale of viscous evolution would be in the order of several billion years. Observed accretion rates and disk lifetimes can only be explained if turbulent viscosity drives the evolution of circumstellar disks. Therefore Shakura & Sunyaev (1973) parameterized the unknown viscosity as the product of a velocity scale and a length scale. The largest reasonable values for these scales in the disk are the pressure scale height

(28)

and the sound speed . Therefore the viscosity is written as

(29)

where is the turbulence parameter and .

Today it is generally believed that disks transport angular momentum by turbulence, however the source of this turbulence is still debated. The magneto-rotational instability is the most commonly accepted candidate for source of turbulence (Balbus & Hawley 1991). Values of are expected to be in the range of to some (see Johansen & Klahr 2005; Dzyurkevich et al. 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)

(30)

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

The effect of these non-axisymmetric structures is to transport angular momentum outward and rearranging the surface density in the disk so as to counteract the unstable configuration. This mechanism is therefore to a certain extent self-limiting. Values above are not influenced by instabilities, values below form instabilities which increase until the disk is marginally stable again. Our model approximates this mechanism by increasing the turbulence parameter according to the recipe of Armitage et al. (2001),

(31)

where is a free parameter of the model which is taken to be , unless otherwise noted.

During the time of infall onto the disk, we use a constant, high value of to mimic the effective redistribution of surface density during the infall phase which also increases the overall stability of the disk. Once the infall stops, we gradually decrease the turbulence parameter on a timescale of 10 000 years until it reaches its input value.

Eq. 26 is a diffusion equation, which is stiff. This means, one faces the problem that the numerical step of an explicit integration scheme goes (where is the radial grid step size) which would make the simulation very slow. One possible solution to this problem is using the method of implicit integration. This scheme keeps the small time scales of diffusion i.e. the fast modes in check. We are not interested in these high frequency modes, but they would become unstable if we used a large time step. With an implicit integration scheme (see Section 10.1) the time step can be chosen larger without causing numerical instabilities, thus increasing the speed of the computation.

7.2 Radial evolution of dust

If the average dust-to-gas ratio in protoplanetary disks is in the order of (as found in the ISM), then the dust does not dynamically influence the gas, while the gas strongly affects the dynamics of the dust.

Thermally, however, the dust has potentially a massive influence on the gas disk evolution through its opacity, but we will not include this in this work. 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 work.

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

(32)

where is the eddy turn-over time and is the stopping time.

The stopping time of a particle is defined as the ratio of the momentum of a particle divided by the drag force acting on it. There are four different regimes for the drag force which determine the dust-to-gas coupling (see Whipple 1972; Weidenschilling 1977). Which regime applies to a certain particle, depends on the ratio between mean free path of the gas molecules and the dust particle size (i.e. the Knudsen number) and also on the particle Reynolds-number with being the gas molecular viscosity

(33)

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

(34)

where denotes the mid-plane number density and  cm.

The different regimes3 are

(35)

Here denotes the velocity of the dust particle with respect to the gas, denotes the mean thermal velocity of the gas molecules, is the solid density of the particles and is the local gas density.

For now, we will focus on the Epstein regime. To calculate the Stokes number, we need to know the eddy turn-over time. As noted before, our description of viscosity comes from a dimensional analysis. We use a characteristic length scale and a characteristic velocity scale of the eddies. This prescription is ambiguous in a sense that it does not specify if angular momentum is transported by large, slow moving eddies or by small, fast moving eddies, that is

(36)

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

(37)

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

(38)

The Stokes number then becomes

(39)

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

(40)

where the total Flux, has contributions from a diffusive and an advective flux. The diffusive part comes from the fact that the gas is turbulent and the dust couples to the gas. The dust is therefore turbulently mixed by the gas. Mixing counteracts gradients in concentration, in this case it is the dust-to-gas ratio of each size that is being smoothed out by the turbulence. The diffusive flux can therefore be written as

(41)

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

(42)

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

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

(43)

where is the radial velocity of the dust. There are two contributions to it,

(44)

The first term is a drag term which comes from the radial movement of the gas which moves with a radial velocity of , given by Eq. 27. 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. 44 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. 44. denotes the maximum drift velocity of a particle,

(45)

which has been derived by Weidenschilling (1977). Here, we introduced a radial drift efficiency parameter . This parameter describes how efficient the radial drift actually is, as several mechanisms such as zonal or meridional flows might slow down radial drift. This will be investigated in Section 8.5.

Putting all together, the time dependent equation for the surface density of one dust species of mass is given by

(46)

where we have included a source term on the right hand side which can be positive in the case of infall or re-condensation of grains and negative in the case of dust evaporation or outflows. This source term does not include the sources caused by coagulation and fragmentation processes (see Section 7.6.1). All sources will be combined into one equation later which is implicitly integrated in an un-split scheme (see Chapter id1).

Note that both, the diffusion coefficient and the radial velocity depend on the Stokes number and therefore on the size of the particle.

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

7.3 Temperature and vertical structure

The vertical structure can be assumed as being in hydrostatic equilibrium at all times if the disk is geometrically thin () and the vertical sound crossing time is much shorter than the radial drift time scale of the gas. The isothermal vertical density structure is then given by

(47)

where

(48)

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

(49)

were denotes the turbulent gas viscosity and the Kepler frequency.

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

(50)

where we used and the approximation . and are Rosseland and Planck mean opacities which will be discussed in the next section.

contains contributions due to stellar or external irradiation. Here, we use a formula derived by Ruden & Pollack (see Ruden & Pollack 1991, App. B),

(51)

with a fixed , following Hueso & Guillot (2005).

The main source of opacity is the dust. Due to viscous heating, the temperature will increase with surface density. If the temperature rises above K, the dust (i.e. the source of opacity) will evaporate, which stops the disk from further heating until all dust is vaporized. To simulate this behavior in our model, we calculate a gas temperature (assuming a small, constant value for gas opacity) in the case where the dust temperature rises above the evaporation temperature. Then is approximated by

(52)

only if from Eq. 50 would be larger than .

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

7.4 Opacity

In the calculation of the mid-plane temperature we use Rosseland and Planck mean opacities which are being calculated from a given frequency dependent opacity table. The results are stored in a look up table and interpolated during the calculations. The opacity table is for a mixture of 50% silicates and 50% carbonaceous grains.

Since these are dust opacities, we convert them from cm/g dust to cm/g gas by multiplying the values with the dust-to-gas ratio , which is a fixed parameter in our model, taken to be the canonical value of 0.01. This assumes that the mean opacity of the gas is very small and that the dust-to-gas ratio does not change with time. To calculate the opacity self-consistently, the total mass of dust and the distribution of grain sizes has to be taken into account, meaning that the dust evolution has a back reaction on the gas by determining the opacity. For now, our model does not take back-reactions from dust to gas evolution into account. Only in the case where the temperature rises above 1500 K, the drop of opacity due to dust evaporation is considered, as described above.

7.5 Initial infall phase: cloud collapse

The evolution of the protoplanetary disk also depends on the initial infall phase which builds up the disk from the collapse of a cold molecular cloud core. This process is still not well understood. First similarity solutions for a collapsing sphere were calculated by Larson (1969) and Penston (1969). Shu (1977) found a self similar solution for a singular isothermal sphere. The Larson & Penston solution predicts much larger infall rates compared to the inside-out collapse of Shu ( and respectively).

More recent work has shown that the infall rates are not constant over time, but develop a peak of high accretion rates and drop to smaller accretion rates at later times. The maximum accretion rate is about if opacity effects are included (see Larson 2003). Analytical, pressure-free collapse calculations of Myers (2005) show similar behavior but with a smaller peak accretion rate of . They also argue that the effects of pressure and magnetic fields will further increase the time scales of cloud collapse.

This initial infall phase is important since it provides the initial condition of the disk and also influences the whole simulation by providing a source of small grains and gas at larger distances to the star during later times of evolution.

It should be noted that several groups perform 3D hydrodynamic simulations of collapsing cloud cores which show more complicated evolution (e.g., Banerjee & Pudritz 2006; Whitehouse & Bate 2006). However, to be able to study general trends of the infall phase, we use the Shu-model since it is adjustable by a few parameters whose influences are easy to understand. In this model the collapse proceeds with an infall rate of which stays constant throughout the collapse.

We assume the singular isothermal sphere of mass to be in solid body rotation . If in-falling material is conserving its angular momentum, all matter from a shell of radius will fall onto the star and disk system (with mass ) within the so called centrifugal radius,

(53)

where is the gravitational constant and . The path of every parcel of gas can then be described by a ballistic orbit until it reaches the equatorial plane. The resulting flow onto the disk is

(54)

where

(55)

and

(56)

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

(57)

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

(58)

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

7.6 Grain growth and fragmentation

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

The goal of the model described in this thesis is to trace the evolution of gas and dust during the whole lifetime of a protoplanetary disk, including the grain growth, radial drift and turbulent mixing.

Here, the problem arises that radial drift and coagulation “counteract” each other in the regime of particles: coagulation of smaller sizes restores the population around , whereas radial drift preferentially removes these particles. To be able to properly model this behavior, the time step has to be chosen small enough if the method of operator splitting is used.

The upper limit for this time step can be very small. If larger steps are used the solution will “flip-flop” back and forth between the two splitted sub-steps of motion and coagulation, and the results become unreliable. A method to allow the choice of large time steps while preserving the accuracy is to use a non-splitted scheme in which the integration is done “implicitly”. B08a already use this technique to avoid flip-flopping between coagulation and fragmentation of grains, and we refer to that paper for a description of the general method. What is new in the current work 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 Chapter id1 for details.

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

(59)

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

(60)

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

(61)

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

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

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

The coagulation and fragmentation kernels will be described in section 7.6.3, the distribution of fragments, , will be described in the next section.

To be able to trace the size and radial evolution of dust in a combined way, we need to express all contributing processes in terms of the same quantity. Hence, we will formulate the coagulation/fragmentation equation in a vertically integrated way. The vertical integration can be done numerically (as in 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 Section 7.6.4, Eq. 76), then the vertical integration can be done analytically, as discussed in Section 10.2. Unlike the steady-state disk models of B08a which have fixed surface density and temperature profiles, we need to recompute the coagulation and fragmentation kernels (which are functions of surface density and temperature) at every time step. Therefore this analytical integration also saves significant amounts of computational time.

We therefore define the vertically integrated dust surface density distribution per logarithmic bin of grain radius, as

(62)

where and are related through . The total dust surface density at any radius is then given by

(63)

Defining as in Eq. 62 makes it a grid-independent density unlike the mass integrated over each numerical bin. This way, all plots of are meaningful without knowledge of the size grid that was used. Numerically, however we use the discretized values as defined in Chapter id1.

In our description of growth and fragmentation of grains, we always assume the dust particles to have a constant volume density meaning that we do not trace the evolution of porosity of the particles as this is currently computationally too expensive with a statistical code as used in this work. This can be achieved with Monte-Carlo methods as in Ormel et al. (2007) or Zsom & Dullemond (2008), however these models have do not yet include the radial motion of dust and therefore cannot trace the global evolution of the dust disk.

Distribution of fragments The distribution of fragments after a collision, , is commonly described by a power law,

(64)

The value has been investigated both experimentally and theoretically. Typical values have been found in the range between 1 and 2, by both experimental (e.g., Blum & Muench 1993; Davis & Ryan 1990) and theoretical studies (Ormel et al. 2009). Unless otherwise noted, we will follow B08a by using the value of .

In the case where masses of the colliding particles differ by orders of magnitude, a complete fragmentation of both particles is an unrealistic scenario. More likely, cratering will occur (Sirono 2004), meaning that the smaller body will excavate a certain amount of mass from the larger one. The amount of mass removed from the larger one is parameterized in units of the smaller body ,

(65)

The mass of the smaller particle plus the mass excavated from the larger one is then distributed to masses smaller than according to the distribution described by Eq. 64. In this work, we follow B08a by assuming to be unity, unless otherwise noted.

Coagulation and fragmentation kernels The coagulation kernel can be factorized into three parts,

(66)

and, similarly, the fragmentation kernel can be written as

(67)

Here, denotes the relative velocity of the two particles, is the geometrical cross section of the collision and and are the coagulation and fragmentation probabilities which sum up to unity. In general, all these factors can also depend on other material properties such as porosity, however we always assume the dust grains to have a volume density of  g cm.

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

(68)

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

Relative particle velocities The different sources of relative velocities considered here are Brownian motion, relative radial and azimuthal velocities, turbulent relative velocities and differential settling. These contributions will be described in the following, an example of the most important contributions is shown in Figure 7.

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

(69)

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

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

(70)

where the radial velocity of the dust, is given by Eq. 44.

Azimuthal relative velocities are induced by gas drag in a similar way as radial drift. However while only particles (plus/minus 2 orders of magnitude) around St =1 are significantly drifting, relative azimuthal velocities do not vanish for encounters between very large and vary small particles (see Figure 7). 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

(71)

where is defined by Eq. 45.

Turbulent motion as source of relative velocities is discussed in detail in Ormel & Cuzzi (2007). They also derive closed form expressions for the turbulent relative velocities which we use in this work.

Differential settling is the fifth process we consider contributing to relative particle velocities. Dullemond & Dominik (2004) constructed detailed models of vertical disk structure describing the depletion of grains in the upper layers of the disk. They show that the equilibrium settling speed for particles in the Epstein regime is given by which can be derived by equating the frictional force and the vertical component of the gravity force, . To limit the settling speed to velocities smaller than half the vertically projected Kepler velocity, we use

(72)

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

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

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

(74)

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

(75)

By limiting the vertical settling velocity as in Eq. 72 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

(76)

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.

8 Results

8.1 Viscous evolution of the gas disk

We will now focus on the evolution of a disk around a T Tauri like star. We assume the rotation rate of the parent cloud core to be  s, which corresponds to 0.06 times the break up rotation rate of the core. The disk is being built-up from inside out due to the Shu-Ulrich infall model, with the centrifugal radius being 8 AU. The parameters of our fiducial model are summarized in Table 1.

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

Figure 8 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 9 shows the evolution of accretion rate onto the star, stellar mass and disk mass. The infall phase lasts until about 150 000 years. At this point, the disk looses its source of gas and small-grained dust and the disk mass drops off quickly until the disk has adjusted itself to the new condition. This also explains the sharp drop of the accretion rate. The slight increase in the accretion rate afterwards comes from the change in after the infall stops (see Section 7.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.

8.2 Fiducial model without fragmentation

Figure 10: Snapshots of the vertically integrated dust density distributions (defined in Eq. 62) of a steady state disk as in B08a (left column) and of an evolving disk (fiducial model, right column). No coagulation is calculated within the evaporation radius (denoted by the dash-dotted line), fragmentation is not taken into account in both simulations. The solid line shows the particle size corresponding to a Stokes number of unity. Since (see Eq. 39) this curve in fact has the same shape as , so it reflects, as a “bonus”, what the gas disk looks like. The radius dividing the evolving disk into accreting and expanding regions is marked by the dotted line and the arrows. Particles which are located below the dashed line have a positive flux in the radial direction due to coupling to the expanding gas disk and turbulent mixing (particles within the closed contour in the upper right plot have an inward pointing flux).
parameter symbol value unit
turbulence parameter -
irradiation angle -
cloud mass
effect. speed of sound in core cm s
rotation rate of cloud core s
solid density of dust grains g/cm
stellar radius
stellar temperature K
Table 1: Parameters of the fiducial model.

We will now focus on the dust evolution of the disk. This fiducial simulation includes only grain growth without fragmentation or erosion. All other parameters as well as the evolution of the gas surface density and mid-plane temperature are the same as discussed in the previous section. The evolution of this model is visualized in Figure 10.

The contour levels in Figure 10 show the vertically integrated dust surface density distribution per logarithmic bin of grain radius, , as defined in Eq. 62. The left column of Figure 10 shows the results of dust evolution for a steady state (i.e. not viscously evolving) gas disk as described in B08a.

The right column shows the evolution of the dust density distribution of the fiducial model, in which the gas disk is gradually built up through infall from the parent molecular cloud core, and the gas disk viscously spreads and accretes matter onto the star. The solid line marks the grain size corresponding to St =1 at the given radius. This grain size will be called hereafter. In the Epstein regime, is proportional to the gas surface density of the disk, which can be seen from Eq. 39.

There are several differences to the simulations of grain growth in steady-state disks, presented in B08a: viscously evolving disks accrete onto the star by transporting mass inwards and angular momentum outwards. This leads to the fact that some gas has to be moving to larger radii because it is ’absorbing’ the angular momentum of the accreting gas. This can be seen in numerical simulations, but also in the self similar solutions of Lynden-Bell & Pringle (1974). They show that there is a radius which divides the disk between inward and outward moving gas. This radius itself is constantly moving outwards, depending on the radial profile of the viscosity.

The radius in the fiducial model was found to move from around 20 AU at the end of the infall to 42 AU at 1 Myr and is denoted by the dotted line in Figure 10. Important here is that small particles are well enough coupled to the gas to be transported along with the outward moving gas while larger particles decouple from the gas and drift inwards. Those parts of the dust distribution which lie below the dotted line in Figure 10 have positive flux in the radial direction due to the gas-coupling.

One might expect that the dotted and dashed lines always coincide for small grains as they are well coupled to the gas, however, it can be seen that this is not the case. The reason for this is that turbulent mixing also contributes to the total flux of dust of each grain size. The smallest particles in between the dotted line and the dashed line in the lower two panels of Figure 10 do have positive velocities, but due to a gradient in concentration of these dust particles, the flux is still negative.

During the early phases of its evolution, a disk which is built up from inside out quickly grows large particles at small radii which are lost to the star by radial drift. During further evolution, growth timescales become larger and larger (since the dust-to-gas ratio is constantly decreasing) while only the small particles are mixed out to large radii.

The radial dependence of the dust-to-gas ratio after 1 Myr is shown in Figure 11. These simulations show that the initial conditions of the stationary disk models (such as shown in B08a and in the left column of Figure 10) are too optimistic since they assume a constant dust-to-gas ratio at the start of their simulation throughout the disk which is not possible if the disk is being built-up from inside out unless the centrifugal radius is very large (in which case the disk would probably fragment gravitationally) and grain fragmentation is included to prevent the grains from becoming large and start drifting strongly. Removal of larger grains and outward transport of small grains lead to the fact that the dust-to-gas ratio is reduced by 0.5–1.5 orders of magnitude compared to a stationary model. This effect is also discussed in Section 8.4.

8.3 Fiducial model with fragmentation

The situation changes significantly, if we take grain fragmentation into account. As discussed in Section 7.6.2, we consider two different kinds of outcomes for grain-grain collisions with relative velocities larger than the fragmentation velocity : cratering (if the masses differ by less than one order of magnitude) and complete fragmentation (otherwise).

Figure 11: Comparison of the radial dependence of the dust-to-gas ratio for the stationary simulations (thick lines) and the evolving disk simulations (thin lines). The four panels compare the results after 1 Myr of evolution with/without fragmentation (left/right column) and with/without effects of non-axisymmetric instabilities (top/bottom row). It can be seen that the dust-to-gas ratio of the evolving disk simulations is almost everywhere lower than the one of the stationary simulations. See Section 8.4 for more details.
Figure 12: Evolution of the dust density distribution of the fiducial model as in Figure 10 but with fragmentation included. The dashed contour line (in the lower two panels) around the upper end of the size distribution and around small particles at  AU marks those parts of the distribution which have a positive (outward pointing) fluxes.

For particles larger than about , relative velocities are dominated by turbulent motion (and to a lesser extend by vertical settling). Since the relative velocities increase with Stokes number (and therefore with grain size), we can calculate the approximate position of the fragmentation barrier by equating the assumed fragmentation velocity with the approximate relative velocities of the particles,

(76)

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. 76 reproduces well the upper end of the size distribution in most of our simulations and therefore helps understanding the influence of various parameters on the outcome of these simulations.

The evolution of the grain size distribution including fragmentation is depicted in Figure 12. The initial condition is quickly forgotten since particles grow on very short timescales to sizes at which they start to fragment. The resulting fragments contribute again to the growth process until a semi-equilibrium of growth and fragmentation is reached.

It can be seen that particles stay much smaller than in the model without fragmentation. This means that they are less affected by radial drift on the one hand and better transported along with the expanding gas disk on the other hand. Consequently, considerable amounts of dust can reach radii of the order of 100 AU, as seen in Figure 13.

The approximate maximum Stokes number, defined in Eq. 76, is inversely proportional to the temperature (since relative velocities are proportional to ), which means that in regions with lower temperature, particles can reach larger Stokes numbers. By equating drift and drag velocities of the particles (cf. Eq. 44), it can be shown that the radial velocities of particles with Stokes numbers larger than about , are being dominated by radial drift.

Due to the high temperatures below 5 AU (caused by viscous heating), stays below this value which prevents any significant radial drift within this radius, particles inside 5 AU are therefore only transported along with the accreting gas. Particles at larger radii and lower temperatures can drift (although only slightly), which means that there is a continuous transport of dust from the outer regions into the inner regions. Once these particles arrive in the hot region, they get “stuck” because their Stokes number drops below . The gas within about 5 AU is therefore enriched in dust, as seen in Figure 13. 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 13 also shows a relatively sharp decrease in the dust to gas ratio at a few hundred AU. At these radii, the gas densities become so small that even the smallest dust particles decouple from the gas and start to drift inwards.

The thick line in Figure 13 shows as comparison the dust-to-gas ratio of the stationary disk model (cf. left column of Figure 10) after 1 Myr, which starts with a radially constant initial dust-to-gas ratio of 0.01.

Figure 14 shows the semi-equilibrium grain surface density distribution at 1, 10 and 100 AU in the fiducial model with fragmentation at 1 Myr. The exact shape of these distributions depends very much on the prescription of fragmentation and cratering. In general the overall shape of these semi-equilibrium distributions is always the same: a power-law or nearly constant distribution (in ) for small grains and a peak at some grain size , beyond which the distribution drops dramatically. The peak near the upper end of the distribution is caused by cratering. This can be understood by looking at the collision velocities: the relative velocity of two particles increases with the grain size but it is lower for equal-sized collisions than for collisions with particles of very different sizes (see Figure 7). The largest particles in the distribution have relative velocities with similar sized particles which lie just below the fragmentation velocity (otherwise they would fragment). This means that the relative velocities with much smaller particles (which are too small to fragment the bigger particles but can still damage them via cratering) are above this critical velocity. This inhibits the further growth of the big particles beyond , causing a “traffic jam” close to the fragmentation barrier. The peak in the distribution represents this traffic jam.

Figure 13: Evolution of the radial dependence of the dust-to-gas ratio in the fiducial model including fragmentation with the times corresponding to the snapshots shown in Figure 12. The initial dust-to-gas ratio is taken to be 0.01. The thick dashed curve represents the result at 1 Myr of the static disk model for comparison.
Figure 14: Vertically integrated (cf. Eq. 62) grain surface density distributions as function of grain radius at a distance of 1 AU (solid), 10 AU (dashed) and 100 AU (dot-dashed) from the star. These curves represent slices through the bottom panel of Figure 12.

8.4 Influences of the infall model

In the fiducial model without fragmentation, continuous resupply of material by infall is the cause why the disk has much more small grains than compared to the stationary disk model (cf. Figure 10), which relatively quickly consumes all available micrometer sized dust. The effect has already been found in Dominik & Dullemond (2008): if all grains start to grow at the same time, then the bulk of the mass grows in a relatively thin peak to larger sizes (see Figure 6 in B08a). However if the bulk of the mass already resides in particles of larger size, then additional supply of small grains by infall is not swept up effectively because of the following: firstly, the number density of large particles is small (they may be dominating the mass, but not necessarily the number density distribution) and secondly, they only reside in a thin mid-plane layer while the scale height of small particles equals the gas scale height.

We studied how much the disk evolution depends on changes in the infall model.

For a given cloud mass, the so called centrifugal radius , which was defined in Eq. 53, depends on the temperature and the angular velocity of the cloud. Both can be varied resulting in a large range of possible centrifugal radii reaching from a few to several hundred AU. Since the centrifugal radius is the relevant parameter, we varied only the rotation rate of the cloud core. We performed simulations with three different rotation rates which correspond to centrifugal radii of about 8 (fiducial model), 33 AU, and 100 AU. For each centrifugal radius, we performed two simulations: one which includes effects of gravitational instabilities (GI) – i.e. increased during infall and according to Eq. 31 – and one which neglects them.

However for a centrifugal radius of 100 AU, too much matter is loaded onto the cold outer parts of the disk and consequently, the disk would fragment through gravitational instability. We cannot treat this in our simulations, hence, for the case of 100 AU, we show only results which neglect all GI effects.

The resulting dust-to-gas ratios are being shown in Figure 11.

Two general aspects change in the case of higher rotation rates: firstly, more of the initial cloud mass has to be accreted onto the star by going through the outer parts of the disk. Consequently, the disk is more extended and more massive than compared to the case of low rotation rate.

Secondly, as aforementioned the high surface densities in the colder regions at larger radii cause the disk to become less gravitationally stable.

If grain fragmentation is not taken into account in the simulations, both effects cause more dust mass to be transported to larger radii. Growth and drift timescales are increasing with radius and the dust disk with centrifugal radius of 33 AU (100 AU) can stay 5 (35) times more massive than in the low angular momentum case after 1 Myr if GI effects are neglected.

If GI effects are included, matter is even more effectively transported outward, the dust-to-disk mass ratio for 8 and 33 AU is increased from 5 to 8.

However if fragmentation is included, it does not matter so significantly, where the dust mass is deposited onto the disk since grains stay so small during the build-up phase of the disk (due to grain fragmentation by turbulent velocities) that they are well coupled to the gas. Outwards of  AU (without GI effects) or of a few hundred AU (if GI effects are included), the gas densities become so small that even the smallest grains start do decouple from the gas. They are therefore not as effectively transported outwards. In these regions, the amount of dust depends on the final centrifugal radius while at smaller radii, turbulent mixing quickly evens out all differences in the dust-to-gas ratio (see left column of Figure 11).

It can be seen, that in all simulations, the dust-to-gas ratio is lower than in the stationary disk model. The trend in the upper right panel in Figure 11 suggests that for a centrifugal radius of 100 AU and the enhanced radial transport by GI effects might have a higher dust-to-gas ratio than the stationary disk model. However in this case, the disk would become extremely unstable and would therefore fragment.

The reason for this is the following: to be able to compare both simulations, the total mass of the disk-star system is the same as in the stationary disk models. How the total mass is distributed onto disk and star depends on the prescription of infall. If a centrifugal radius of 100 AU is used, the disk becomes so massive that it significantly exceeds the stability criterion .

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

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

(77)

where is the dust-to-gas ratio and and are the drift and growth timescales respectively. The parameter describes how much more efficient growth around St =1 must be, so that the particles are not removed by radial drift. To overcome the drift barrier, obviously either particle growth must be accelerated, or the drift efficiency has to be decreased. B08a have numerically measured the parameter to be around 12. In other words, this means that the growth timescales have to be decreased (e.g. by an increased dust-to-gas ratio) until the condition in Eq. 77 is fulfilled.

However, there are other ways of breaking through the drift barrier. Firstly, the drift timescale for St =1 particles also depends on the temperature (via the pressure gradient). A simple approximation from Eq. 77 with a 0.5  star and a dust-to-gas ratio of 0.01 gives

(78)

which means that particles should be able to break through the drift barrier at 1 AU if the temperature is below 100 K (or 200 K for a solar mass star). Dullemond et al. (2002) have constructed vertical structure models of passively irradiated circumstellar disks using full frequency- and angle-dependent radiative transfer. They show that the mid-plane temperature of such a T Tauri like system at 1 AU can be as low as 60 K. Reducing the temperature by some factor reduces the drift time scale by the a factor of similar size which we will call the radial drift efficiency (cf. Eq. 45).

Zonal flows as presented in Johansen et al. (2006) could be an alternative way of decreasing the efficiency of radial drift averaged over typical time scales of particle growth. Johansen et al. (2006) found a reduction of the radial drift velocity of up to 40% for meter-sized particles.

Meridional flows (e.g., Urpin 1984; Kley & Lin 1992) might also seem interesting in this context, however they do not directly influence the radial drift efficiency but rather reverse the gas-drag effect. This might be important for small particles (which, however are not strongly settling to the mid-plane) but for St =1 particles, would have to be extremely high to have significant influence: even would result in a reduction of the particles radial velocity by approximately only a few percent.

Another possibility to weaken the drift barrier is changing its radial dependence. The reason for this is the following: particle radial drift is only a barrier if it prevents particles to cross the size . Since particles grow while drifting, the particle size corresponding to St =1 needs to increase as well, to be a barrier. Otherwise drifting particles would grow (at least partly) through the barrier while they are drifting. If is decreasing in the direction toward the star, then a particle that drifts inwards would have an increasing Stokes number even if the particle does not grow at all.

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

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

(79)

meaning that a relatively flat profile of surface density (or even a profile with positive slope) is needed to allow particles to grow through the barrier. However, our simulations of the viscous gas disk evolution does not yield surface density profiles with positive slopes outside the dust evaporation radius.

To quantify the arguments above, we have performed simulations with varying drift efficiency to test how much the radial drift has to be reduced to allow break through. We have additionally included the first Stokes drag regime to see how the radial drift of particles is influenced by it.

Figure 15 shows the grain surface density distribution after 200 000 years of evolution for three different drag efficiencies. The most obvious changes can be seen in the region where the line (solid line) is relatively flat: the grain distribution is shifted towards larger Stokes numbers. As explained above, particles can grow while drifting, which can be seen in the case of . The Stokes number and size of the largest particles is significantly increasing towards smaller radii. However the radial drift efficiency has to be reduced by 80% to produce particles which are large enough to escape the drift regime and are therefore not lost to the star.