Particle evolution in protoplanetary discs

Global variation of the dust-to-gas ratio in evolving protoplanetary discs


Recent theories suggest planetesimal formation via streaming and/or gravitational instabilities may be triggered by localized enhancements in the dust-to-gas ratio, and one hypothesis is that sufficient enhancements may be produced in the pile-up of small solid particles inspiralling under aerodynamic drag from the large mass reservoir in the outer disc. Studies of particle pile-up in static gas discs have provided partial support for this hypothesis. Here, we study the radial and temporal evolution of the dust-to-gas ratio in turbulent discs, that evolve under the action of viscosity and photoevaporation. We find that particle pile-ups do not generically occur within evolving discs, particularly if the introduction of large grains is restricted to the inner, dense regions of a disc. Instead, radial drift results in depletion of solids from the outer disc, while the inner disc maintains a dust-to-gas ratio that is within a factor of of the initial value. We attribute this result to the short time-scales for turbulent diffusion and radial advection (with the mean gas flow) in the inner disc. We show that the qualitative evolution of the dust-to-gas ratio depends only weakly upon the parameters of the disc model (the disc mass, size, viscosity, and value of the Schmidt number), and discuss the implications for planetesimal formation via collective instabilities. Our results suggest that in discs where there is a significant level of midplane turbulence and accretion, planetesimal formation would need to be possible in the absence of large-scale enhancements. Instead, trapping and concentration of particles within local turbulent structures may be required as a first stage of planetesimal formation.

accretion, accretion discs — planets and satellites: formation — protoplanetary discs — stars: pre-main-sequence — stars: variables: T Tauri, Herbig Ae/Be

1 Introduction

Uncovering the pathway by which sub-micron dust grains grow into gravitationally bound bodies suitable for building planets is a long-standing problem. Pairwise collisions between dust particles in protoplanetary discs are thought to drive rapid growth up to mm sizes (Dullemond & Dominik, 2005). Beyond this size scale, continued growth up to the km-scale of planetesimals is potentially frustrated by two distinct physical barriers. First, the measured material properties of solid aggregates, which favor sticking upon collision for small sizes, switch to favour bouncing (Güttler et al., 2010; Zsom et al., 2010), or collisional fragmentation due to high, turbulence-driven velocity dispersions (e.g. Birnstiel, Dullemond & Brauer, 2010) for . Even if these barriers could be surmounted, however, direct collisional growth faces a second obstacle from headwind drag (Weidenschilling, 1977), which becomes extremely strong for particles in the cm to m-size range. The loss of substantial quantities of solids into the star is not in itself deleterious – indeed the relatively small solid mass in the Minimum Mass Solar Nebula (Weidenschilling, 1977b) indicates loss probably occurred – but the residence time of critically coupled particles is so short that it is hard to see how catastrophic depletion could be averted if bodies grow steadily through m-scales.

These considerations have led to renewed interest in collective mechanisms for planetesimal formation, which appeal (typically) to gravitational instability in a dense particle layer (Safronov, 1972; Goldreich & Ward, 1973) to effect a dynamical time scale jump from mm or cm-sizes to km-scales (for a review see Chiang & Youdin, 2010). Several precursor mechanisms have been proposed that might act to create conditions locally favourable for gravitational instability. They include vertical settling (Youdin & Shu, 2002), turbulent concentration (Cuzzi, Hogan & Shariff, 2008; Johansen et al., 2007), and streaming instabilities (Johansen & Youdin, 2007; Johansen, Youdin & Mac Low, 2009b; Bai & Stone, 2010a). There is some evidence that these mechanisms could form surprisingly large first-generation planetesimals (), of a size consistent with that derived — by combining the measured size distribution with the results of collisional-grinding simulations — for the primordial asteroid belt (Morbidelli et al., 2009).

Which, if any, of these mechanisms is responsible for planetesimal formation is hard to say, since their feasibility depends in part on the uncertain properties of the intrinsic turbulence within the disc (Armitage, 2011). Some theoretical constraints are nonetheless possible, since most of the formation-via-collapse models require local enhancement of the height-averaged dust-to-gas ratio (the local disc metallicity) above solar values in order to exceed the threshold for gravitational instability 1. The magnitude of the required enhancement varies both with the model, and with the properties of the gas disc, but typical values quoted vary between roughly 2 and 10 times standard metallicity (Youdin & Shu, 2002; Johansen, Youdin & Mac Low, 2009b; Bai & Stone, 2010a). Either a local build-up of solids or a local depletion of gas (Throop & Bally, 2005; Alexander & Armitage, 2007) would be equally effective.

Further clues come from observations. To match Solar System constraints, conditions for planetesimal formation must extend over a significant radial range (because planetesimals formed with substantially distinct compositions) and persist for a large fraction of the disc lifetime. The evidence for an extended era of planetesimal formation comes from age dating of primitive meteorites, and in particular from the finding that the chondrules within those meteorites differ in age by as much as 2.5 Myr (Amelin et al., 2002). Furthermore, recent age-dating of iron-core meteorites suggests the era of planetesimal formation extended nearly as far back as the beginning of Solar-system history, as defined by the formation of calcium-aluminium-rich inclusions (CAIs) (Kleine et al., 2005; Bottke et al., 2006).

The replenishment of cold dust in debris discs (Wyatt, 2008) implies that planetesimal populations around other stars are not solely confined to the inner disc, but there is no observational requirement that planetesimal formation spans the entire range of radii where gas was once present. Given this, one obvious way to create an enhancement of the dust-to-gas ratio at small to intermediate disc radii is to appeal to the large reservoir of mass typically present further out. If small particles spiralling inward pile-up as they migrate radially, then the “problem” of headwind drag becomes instead a necessary ingredient for planetesimal formation. In a static (inviscid) disc model, headwind drag acting on an initially uniform dust-to-gas ratio does in fact result in just such a pile-up, with peak values of the enhancement exceeding those needed to produce gravitational collapse via vertical sedimentation (Youdin & Shu, 2002; Youdin & Chiang, 2004). Although not the focus of the original work, radial pile-ups would also result in highly favorable conditions for triggering streaming instabilities.

Inviscid-disc models, which assume that the gas has neither a net radial flow nor local turbulent motions, include only the minimal set of potential aerodynamic effects that can modify the radial particle distribution. Such models also provide a poor representation of the global disc environment, which is expected (and observed, e.g. Hartmann et al., 1998) to change substantially over the time period of interest for planetesimal formation. Our goal in this paper is to quantify the variation in the dust-to-gas ratio if the gas disc, instead of being static, evolves with time under the action of local turbulent transport of angular momentum and photoevaporation. Headwind drag will still occur in an evolving disc, but the propensity of particles to pile-up will additionally be affected by several new processes,

  • The global evolution of the gas disc, which thins with time and whose outer edge must expand radially to conserve angular momentum as gas is accreted on to the star.

  • Radial advection with the flow of disc gas, which may counter the slow-down of drift that leads to pileup in the inner disc.

  • Turbulent diffusion.

  • Disc-clearing photoevaporation.

The first three of these processes (also explored in, e.g., Stepinski & Valageas, 1996) are ultimately byproducts of the same turbulence within the disc, so the inclusion of evolution introduces fewer independent parameters into a calculation of particle transport than one might at first fear. Indeed, the primary conclusion of our work is that a range of disc models all produce roughly the same result: as the disc evolves, the outer disc becomes depleted of solids, but the lost material does not pile-up significantly in the inner disc. Instead, the inner disc maintains roughly its primordial dust-to-gas ratio as the disc drains away. As time progresses, the thinning gas supports successively smaller particle sizes across the range of Solar System radii.

The plan of this paper is as follows. In §2 we outline the one dimensional disc evolution and particle transport models that we employ. These models immediately specify the radial dependence of several critical time scales relevant to particle transport, and in §3 we discuss the elementary consequences of those time scales for dust evolution. In §4, we show the dust distributions resulting from our simulations, which we compare in §5 to the dust-to-gas enhancements thought necessary to lead to precipitation (formation via vertical sedimentation as in Youdin & Shu (2002)) of large bodies, and, briefly in §6, to observations of dust in discs around other stars. We discuss the implications of our results for planetesimal formation models in §7.

2 Methods

2.1 Gas disc evolution

Our simulations assume a 1D, viscously evolving disc model, whose surface density profile is at small radii, with a characteristic exponential fall-off from ,


where is the surface density of the disc gas, is the initial disc mass, is radial distance from the central star, and  AU is the inner boundary of the disc-model grid. In the fiducial model, and  AU, but simulations are also run in discs with and or  AU and 40 AU.

Viscous evolution of the disc follows the thin-disc model of disc evolution (Pringle, 1981) on a 1D grid of 600 points spaced logarithmically in from 0.1 to 15,000 AU. We use the Shakura & Sunyaev (1973) alpha-parameterization for the disc viscosity,


where is the viscosity, is the local gas sound speed, is the local Keplerian angular velocity, and , is the vertically-isothermal disc-gas scale-height. is Boltzman’s constant, and is the average mass of a gas particle, set equal to 2.34 hydrogen atoms. In the fiducial model, , but simulations are also run in discs using and . The disc-evolution model also includes photoevaporation at the end of the disc lifetime due to photons s of EUV radiation, as in Hughes & Armitage (2010), Equation (1).

We assume that the viscosity in Equation (2) derives from a local turbulent stress (Balbus & Papaloizou, 1999). In this case, both the evolution of the gas disc profile, and the magnitude of the particle diffusivity, , share a common scaling with the midplane disc temperature, . Due to gas infall and accretional heating, disc temperatures at early times are much hotter than those of later discs (Cassen, 2001). To evaluate the evolution of the midplane temperature, we assume that local heating by viscous dissipation and stellar irradiation is balanced by radiative cooling. We compute these terms using standard methods (detailed in Appendix A), assuming a constant dust-to-gas ratio (i.e., ignoring coupling to the evolving dust distribution). For computational economy, when running our dust-transport simulations, we model the gas-disc evolution using a fit to the computed energy-balanced temperature evolution. At a given time,


where  K is the background temperature of the star-forming cloud, and and follow an exponential decay in time from an initially hotter, steeper distribution ( K, ) to a cooler, shallower distribution appropriate for a passively-heated, flared disc ( K, ). To compute the stellar irradiation, which dominates at late times and large radii, we assume a stellar luminosity of and a disc-starlight intersection angle set to everywhere. We have included the time-evolving nature of the disc-temperature distribution for completeness in considering the fundamental differences between an evolving versus a static disc-mass distribution. However, at the time-scales / masses and accretion rates considered in our models, disc surface-density evolution is nearly identical using either a static or evolving temperature distribution. We have also verified that using a fit to the evolving temperature distribution, rather than calculating it in each dust-transport simulation at every , does not affect our dust-distribution results.

The disc-gas radial velocity, , is calculated using the 1D accretional velocity derived due to the viscous evolution of the disc surface-density profile.


This accretion-flow velocity is predominantly inward but includes outward flow at the outer edge of the expanding disc. The inclusion of the accretion velocity of the gas marks one of the primary differences between our evolving-disc model and the inviscid-disc model of Youdin & Shu (2002). It is important because, for a typical gas surface density profile (), does not decrease toward smaller disc radii, unlike the headwind-drag velocity. We note, however, that it is possible for discs to evolve without there being a significant accretion flow at the disc midplane. This situation occurs in dead-zone models (Gammie, 1996; Armitage, 2011), in which vertical variation in the efficiency of angular momentum transport results in gas accretion that occurs primarily along the disc surface. Qualitatively, we expect that particle transport within dead-zone disc models would differ substantially depending on whether the scale height of the particle disc was large enough for the particles to experience the (rapid) surface layer flow. However, more detailed study of these models is left to others.

The gas azimuthal velocity, , effecting headwind-drag on the simulation grains, is calculated at the disc midplane assuming pressure balance within the vertically-isothermal disc model (Takeuchi & Lin, 2002). For the numerical calculation, we use


where is the gas angular velocity, is the power-law index of the temperature distribution used in Equation (3), and is the locally-calculated power-law index of the disc surface-density profile.

2.2 Particle Transport Setup

Particle Motions

We use the particle-transport code outlined in Hughes & Armitage (2010) to evolve the dust distribution, which is modeled as an ensemble of particles. Each particle is subject to aerodynamic advection and random-walk turbulent diffusion within the disc. Particle trajectories are integrated following


where is a particle’s radial position, the time-step, is a particle’s radial velocity due to gas-drag and advection within the mean gas flow, and is the mean velocity for effecting random-walk turbulent diffusion of the particle ensemble. The random walk is asymmetric, with the weighting between the inward and outward steps calculated from the local disc conditions so as to recover the evolution of a fluid model (Hughes & Armitage, 2010).

Gas drag is calculated assuming that particles experience Epstein drag appropriate to their size. The Epstein drag regime is appropriate for the vast majority of the disc’s radial and temporal extent, though the largest grains considered (cm-sized) technically fall within the Stokes-drag regime in the innermost disc ( AU) at very early times. The radial velocity of a grain is calculated at each radial grid-point by iteratively solving for the steady velocity of particles that have attained terminal velocity (Takeuchi & Lin, 2002; Hughes & Armitage, 2010),


Here is the dust azimuthal velocity, is the local Kepler velocity, is the dust-grain surface-area-to-mass ratio, is the local gas density (set equal to the value at the midplane in a vertically isothermal disc), and is the gas thermal velocity. We assume spherical grains of internal density 1 g cm, and vary the grain radius, . We consider grain sizes of m – 2 cm, which is meant to represent the general range of grain sizes produced by coagulative grain-growth processes. Above some critical size, laboratory experiments and coagulation models show that bouncing and collisional fragmentation represent barriers to further growth (Brauer, Dullemond & Henning, 2008; Zsom et al., 2010). The precise value of the critical size depends upon the gas disc conditions, and is not precisely known. We assume 2 cm as the upper-limit grain size for tracking a global dust distribution. Except in the case of instantaneous grain growth (discussed in § 4.2), our simulation results are qualitatively identical for maximum assumed grain sizes as small as 0.2 mm.

The local velocity for turbulent diffusion is,


where is the local diffusivity of the particle ensemble. We use the Youdin & Lithwick (2007) scaling for the particle radial diffusivity,


where is the gas-particle diffusivity, and is the normalized gas-drag stopping time. In general we assume a Schmidt number of , but also run simulations in the fiducial-disc model using and 2. Equation (10) has been numerically validated by Carballido, Bai & Cuzzi (2011). However, note that for , it deviates from the earlier common expression given by Cuzzi, Dobrovolskis & Champney (1993): , where , and is a normalized eddy turnover timescale. We have run a comparison simulation of the evolving dust distribution within the fiducial disc model using this other dust-diffusivity parameterization (with commonly used ) and have found that the results are qualitatively identical to those using the parameterization in Equation (10). As discussed in §3, headwind-drag generally confines grains to regions of the disc with .

Initial Particle Distribution and Dust Mass Allocation

We assume that, at , the entire mass of disc solids is in small particles (0.2 m) whose radial distribution matches that of the disc gas. Since coagulation occurs rapidly for such small-size particles, to a first approximation we expect that the grain-size distribution will rapidly approach an equilibrium dictated by collision/fragmentation processes (e.g., Birnstiel, Ormel & Dullemond (2011)). In more detail, however, the outer regions of a disc will experience a time-delay in the production of large grains relative to the inner, denser regions. This should translate into a delay before which rapid inspiral of large grains can occur from the outer disc.

The grain trajectories calculated in our simulations are non-coagulative, considering a single grain-size for the full lifetime of a given run. In order to account for the delay in the appearance of larger grain sizes, we employ a pair of simplified grain-growth models (described in Appendix B) to define a radially dependent timescale for the appearance of larger grains. Larger grains are then initialized at time within ten radial zones for which grain-growth time scales (at characteristic radii) have been calculated. This initialization of all grains larger than 0.2 m diameter occurs only once per zone, explicitly neglecting the collisional or coagulative regeneration and resupply of any grain-size population that may be removed due to transport effects. However, as discussed in §4, limiting the appearance of large grains (especially in the inner disc where they would be most likely to regenerate) has the primary effect of retaining that dust mass against loss inward onto the parent star. This simplification therefore results in the highest dust-to-gas ratio in the inner disc achievable with our disc-evolution and dust-transport models.

Our simplified grain-growth models consider coagulation via differential settling (in the rapid-upper-limit case in the absence of turbulence) and random motions. We do not consider growth via differential-radial drift, which, e.g., Stepinski & Valageas (1996), found to produce more rapid growth time-scales, nor directly via turbulent velocity dispersions. Our assumed dust distributions, therefore, may represent a lower-limit on the rapidity and radial extent of large-grain initialization. In §4.2, we consider the upper limit cases where large-grain formation is a) instantaneous, and b) independent of disc location (not confined to the inner regions of the disc) in order to define the boundaries of the behavior of our evolving dust/gas distributions within the grain-growth–model parameter space.

Figure 1

Figure 1: Timing for the initiation of larger grains within our fiducial simulations for each of our radial zones. These correspond to the more rapid of the two growth timescales produced by our cartoon grain-growth models, discussed in Appendix B. Points not plotted beyond a given radial position indicate that grains are not grown to that size within the disc lifetime.

plots the grain-growth timing employed within our fiducial disc simulations. We run our two cartoon grain-growth models independently and initiate larger grain sizes according the most rapid growth timescale indicated. Due to disc evolution and depletion, sometimes only one (or neither) model will produce grains of a given size within a given zone. Our ten zones have characteristic radii spaced logarithmically from 0.5 to 1000 AU, and zone boundaries placed logarithmically between. We simulate transport of 0.2 m–2cm-sized grains, spaced per decade in grain size, running an ensemble of grains for each grain-size transport simulation. Grains initialized within a zone are placed randomly following the gas-mass distribution within that zone.

Each simulation particle tracks an assigned gas-equivalent mass, , of solid material within the disc. The gas-equivalent mass can be used to directly measure the local change in the dust-to-gas ratio (the dust enhancement factor, ) within the disc. By adding up all gas-equivalent mass within a radial bin, , and comparing that to the disc-gas mass within that bin, , we can calculate the locally-measured dust-to-gas enhancement factor within our simulations.


for indicates that the local disc metallicity is unchanged from its initial (solar) value, whereas indicates an enhancement in the local dust-to-gas ratio.

Our simulations do not explicitly track refractory and icy material separately or include evaporation/condensation effects for grains crossing the snow-line. Instead, the composition (and therefore relative mass of grains) is assumed to be that appropriate for the local disc-temperature conditions. While the enhancement factor is, therefore, the fundamental output of our dust-transport simulations, it can be translated into a local value for the dust-grain surface density, , separate from the disc gas. To do this, we solve for at a given binned using those assumed compositions and the more general definition of the enhancement factor,


where refers to the initial (solar) dust-gas composition of the disc, the baseline metallicity. used in this work is taken from Lodders (2003) and is broken into two components:


where is the total fraction of condensable material thought to be present in a solar-composition gas, and is the fraction of that material believed to actually be condensed based on the local gas temperature. The values of used are

with transitions corresponding to the condensation of water ice (the snow line) and to the condensation of methane (and ammonia) ice.

Because initially all disc solids are contained within 0.2 m grains, at all 0.2 m grains are assigned . When a new grain size, , is initialized within a radial zone, all values for grains 0.2 m – presently within that zone are reassigned to follow an updated mass allocation for the grain-size distribution. For this mass re-allocation, we assume a power-law distribution of grain sizes,


where is the number density of grains of a given radius , and is the power-law index of the distribution. For the fiducial simulations, we take , corresponding to the size distribution of grains measured for the ISM (Mathis, Rumpl & Nordsieck, 1977). However, observations linked to grain-growth within discs often predict flatter values of (e.g. Ricci et al., 2011), and in §4 we also report results that include and -4, the later commonly used to depict a purely collisional size distribution, e.g., Bai & Stone (2010a).

To calculate the fraction of the total dust mass within a zone to be reallocated into a given grain-size bin, we follow Garaud (2007), where, for a general value of ,


where is the local volume density of all solids, is the volume density of solids within size bin , and are the lower and upper bounds, respectively, on the entire size distribution, and and are the bounds on a given size bin. For ,


We define size bins logarithmically, so that, e.g., redistributing mass across 0.2 m and 2 m size bins ( and 1 m) has bin boundaries at , 0.316, and 3.16 m. Once we know the fraction of gas-equivalent mass to be portioned into each size bin, we divide it evenly among the particles of that size currently within the zone being updated.

Note that because our dust-transport trajectories are non-coagulative, tracking only a single grain size per model run, grain growth in our simulations in each radial zone occurs only when our grain-growth models indicate that the next largest grain size should be initiated. At that time, all of the dust mass currently within that zone is summed and redistributed according to Equations (15) or (16), thereby updating the values of all particles (including those newly initiated) currently within that zone. We therefore maintain a roughly equilibriated grain-size distribution within our disc during the early period of growth up to large-grain sizes. However, as radial drift proceeds and preferentially depletes mass from the larger grain-size bins, this equilibration is not maintained. Furthermore, we assume that larger grains will be built from grains only slightly smaller than themselves. Therefore, in terms of mass reallocation, we do not initiate new, larger grains in a zone if radial drift has removed all grains of the next smallest size from that zone.

Without a mechanism to replenish the large-dust population, our simulations therefore produce a conservative estimate of the degree of dust-mass infall experienced within an evolving disk. As (discussed in §4) grain growth and the infall of large grains primarily leads to the draining of dust-mass from the disc, we expect that our simulations represent rough upper-limits on local dust-to-gas enhancements in the inner/main disc at times greater than the earliest disc-evolution timescales ( Myr.) For possible alternatives at earlier times, see §4.2.

Finally, we report the results for our dust-distribution simulations as a map of dust-to-gas enhancement, , within radial and temporal space, . While our simulations report mass values for each radial grid-point (600 total) of our model disc and at each output time-step, usually  years, we present our results at lower resolution in order to reduce noise from our particle-ensemble representation of the dust distribution. We present results using mass radially summed across three grid cells for 200 radial grid-points of results output, and averaged across usually 250 time-steps for a results  years.

3 Relevant Timescales and Disc-model Constraints

We begin by considering several timescales relevant to grain transport, which are plotted in Figure 2

Figure 2: Comparison of dust-transport and grain-growth timescales for 0.2 mm grains within the fiducial disc model. Timescales for inward advection, inward drift, and diffusion of grains are plotted for  years of the fiducial disc model. Timescales for grain-growth are from and include disc evolution. There is a break in the plotted advection timescale where the dust flow is radially outward. The dashed line references the disc lifetime.

for 0.2 mm grains at early times within the fiducial disc model. The first is the time-scale to grow grains to a given size, for which we plot the initiation timing taken from the models shown in Figure 1. Next is the time-scale for grains to spiral inward onto the parent star. This can be approximated as the distance to the star divided by the local inward-drift/drag velocity. However, there are actually two inward velocities to consider. The first is the head-wind drag velocity by itself, which is commonly used in such loss estimates applied to static disc models. The second is the drag-advection velocity on the grain, calculated as from Equations (7) & (8) above, which includes the (mostly) inward gas flow resultant from disc evolution. Within our model, the drift-only velocity, , may be calculated with the same method as only setting . The two inward-loss timescales are therefore defined as


The final timescale to consider is the timescale for turbulent/diffusive redistribution of the particle ensemble, which for Figure 2 we define as


From the timescales plotted in Figure 2, one sees that there are several important trends to consider.

  1. The shortest time-scale is that for the appearance of these 0.2 mm grains within the disc, out to at least a few tens of AU. At any one radius, then, it would be reasonable to assume that the initial phases of grain growth are effectively instantaneous. However, when we consider the disc as a whole, we find that even our lower limit on the growth time scale in the outer disc exceeds other time scales in the inner disc. This implies that by the time sub-mm grains have formed at tens of AU, we expect the major fraction of similar grains to have been lost from the inner few AU.

  2. There is a large difference between the time-scales and . The inward drift time scale, which would represent the flow of particles in an inviscid disc, is relatively long and almost constant with . These features are consistent with the idea that radial drift results in a pile-up of solids in the inner disc. In an evolving disc, however, the advection time scale is substantially shorter everywhere inside 10 AU. Solids can therefore be lost via advection along with the gas, without having had time to pile-up.

  3. Finally, the advection and diffusion timescales are similar, except in the outer disc where headwind-drag infall dominates. Turbulent redistribution of grains can therefore matter in the inner disc. Since turbulence represents a diffusive process, it will act to limit sharp radial gradients in particle concentration.

For planetesimal formation, two basic properties of the background disc model are particularly important: the solid-gas coupling strength for particles of a given size, and the magnitude of the deviation from Keplerian rotation due to gas pressure gradients. These quantities enter into models for planetesimal formation via either streaming instabilities, or vertical precipitation followed by collapse. It is therefore useful to discuss how they vary in the gas disc model that we assume.

Local dust-gas coupling is parameterized using the local gas-drag stopping time on particles within the disc. This depends on the size of the particle considered, but for characterizing the disc, we can plot the local particle-size corresponding to , which is a stopping time of particular interest to some theories of planetesimal formation. In the Epstein-drag regime


where and are the dust-grain radius and internal density, is the local gas density, and is the thermal velocity of the gas particles. At the midplane in a vertically-isothermal disc model, Equation (20) translates to


In Figure 3,

Figure 3: grain-size contours (m–2 m) overlying our simulated distribution of 0.2 mm grains within the fiducial disc model. Enhancement values of the distribution are scaled to the total simulated dust-grain distribution, of which 0.2 mm grains make up some fraction. Note that contours follow contours in a vertically isothermal disc (Eq. 21).

we plot the contours for  g cm grains within the fiducial disc model, overlain on the 0.2 mm component of our simulated dust-grain distribution. We have highlighted the contour corresponding to for 0.2 mm grains, and one can see that the dust distribution within the disc is confined well interior () to this contour. This is consistent with the time-scales for radial drift presented in Figure 2 ( occurs near 200 AU at the peak in the head-wind-drag effect), and also with the simulations of Brauer, Dullemond & Henning (2008), though radial confinement in our fiducial simulations is about an order-of-magnitude greater than theirs, perhaps due in part to our approximate treatment of grain growth. We find identical behavior (confinement to regions of the disc) for the simulated distributions of each dust-size component we have run, including only minor variation if the diffusivity () is varied by a factor of two. Headwind drag therefore confines larger dust solids to the inner/main disc.

The precipitation mechanism of planetesimal formation explored by Youdin & Shu (2002) is formally independent of the local dust size distribution (though grains do need to be large enough to settle toward the disc midplane, particularly if intrinsic disc turbulence is important for the vertical dust profile). However, streaming-instability models for planetesimal formation find that particle clumping and collapse occurs specifically for particles that are marginally coupled to the gas motions (Johansen et al., 2007; Johansen & Youdin, 2007; Bai & Stone, 2010a). This means particles with normalized stopping times, , preferentially near 1 (though Bai & Stone (2010a) report that grains as small as still participate in the streaming instability). Because the headwind-drag effect dictates the rapid removal of grains from their disc locations, it is hard to reconcile planetesimal formation models that require with aerodynamic grain-transport theory. Observationally, however, chondrules (mm in size) do make up a large fraction of most meteorites (Cuzzi, Davis & Dobrovolskis, 2003). Planetesimals in the inner region of the Solar Nebula may then have been largely built from mid-to-large sized grains (see, e.g., Cuzzi, Hogan & Bottke (2010) for a discussion), and planetesimal formation models using such grain sizes may be well justified.

Figure 4: contours (red) overlying grain-size contours (black) for the fiducial-disc model. contours are for values of 0.05, 0.1, 0.25, 0.5, and 1. grain-size contours are spaced in decades 0.2 m–2 m.
Figure 5: Composite map of dust/gas enhancement factor (top panel) within the fiducial disc model and maps of some of the contributing grain-size bins (bottom panels). Regions colored in black have a solids/gas enhancement factor of or greater. The black contours trace the evolving gas surface density, , spaced by orders of magnitude between and g cm. The grey regions mask out bins for which the particle statistics were too low to report dust/gas enhancement values.

Finally, a requirement of planetesimal formation common to both precipitation and streaming-instability models (and consistent with the model presented in Cuzzi, Hogan & Bottke (2010)) is that the local radial pressure gradient within the disc is small. The radial pressure gradient can be parameterized via,


since it is the source of the offset between gas and Keplerian orbital velocities. The dependence of the Youdin & Shu (2002) precipitation mechanism on low is explicitly outlined in §5. Some papers reporting simulations of streaming-instability clumping of particles instead use , or


and in Figure 4 we plot contours in red overlying the grain-size contours for the fiducial-disc model. In Bai & Stone (2010b), the authors report simulations of particle clumping via the streaming-instability in which they vary between 0.025, 0.05, and 0.1. They find clumping requires for most particles and super-solar metallicity when , and a monotonic increase in required metallicity for increasing . Comparing these results to the spatial distribution of values within our disc model, it is clear that even the purely gaseous properties of the disc work to disfavor planetesimal formation at large radii. It remains to be explored whether potentially significant dust-to-gas enhancements above solar can be produced closer in.

4 Results: The Global Distribution of Solids

4.1 Fiducial Results

Figure 5 shows the simulated solids enhancement (change in the dust-to-gas ratio) as a function of and within the model disc. It is overlain with contours showing the evolution of the gas surface density. The bottom panels of the figure plot contributions from the individual grain-size components that make up the total solid mass. (Note that the green, scallop-shaped enhancement contours in the inner/main disc in the composite map are artifacts of the size-binning, matching the terminal-infall boundaries of successively larger grain sizes.) In these maps, green areas represent no change from the dust-to-gas ratio, red areas represent depletion of dust relative to the gas, and black areas represent enhancement of the dust or greater.

From the figure, we see that the fiducial-disc simulations do not produce large radial concentrations of dust in the main and inner discs until near the end of the disc lifetime. In §5, we will compare these measured enhancement values to predicted requirements for the precipitation model of planetesimal formation. Recall, however, that a common minimum requirement quoted to aid planetesimal formation is a doubling of the dust-to-gas ratio above solar metallicity, and that, as discussed in §1, it is the first few million years of disc evolution that are of greatest interest in terms of bulk planetesimal formation. Potentially significant enhancements are not produced by our models during this phase of disc evolution.

Figure 5 illustrates several other characteristics of the simulations. At early times, a noticeable portion of the small-grain population is advected outward as the disc expands. However, inward drift of particles quickly sets in, and leads to the depletion of dust relative to the gas in the outer regions of the disc. In the inner disc, this influx of mass is in turn offset by concurrent loss (due to further inspiral) of most of the larger-grain population. The net result is a significant loss of dust mass from the evolving disc.

The significance of the mass contained within the larger grain sizes is depicted in Figure 6

Figure 6: Measured enhancement of local dust/gas for the fiducial-disc model simulations at Myr varying the upper-limit on the grain-size distribution. For no grain growth (0.2 m grains only) the dust enhancement factor remains everywhere near unity. Allowing larger grain sizes allows inward drift and some pile-up of the dust population.

where, for Myr, we plot the enhancement factor produced when the results of our fiducial simulations are compiled (using sub-sets of the simulation output) for grain-size distributions truncated at successively smaller particle sizes. Here, we see that grain-growth to roughly 20 m sizes results in inward drift that depletes the outer disc of dust, and it is for grain-size truncation at m that our simulations produce peak enhancement ratios within the inner/main disc. Growth to larger sizes, however, subsequently depletes the inner disc of this dust-mass excess, bringing enhancement values back toward unity.

Finally, we can solve for a representation of the local surface-density of dust particles, , using Equations (12) and (13) in order to gain a sense of the diminution of solids available in an evolving, accreting disc as time progresses. In Figure 7,

Figure 7: Time evolution of (a) the dust surface-density profile, and (b) the dust mass-per-radius distributions for the fiducial-disc simulations.

we plot the time evolution of and of the local dust-mass per AU within the fiducial-disc simulations. The snow-line is clearly marked by jumps in the distributions at around a few AU. At , our fiducial disc has just over 100 Earth-masses of solid material spread throughout the entire disc, with the greatest concentration just beyond the snow line at around 5–10 AU. As the disc evolves, the solids distribution spreads outward with the disc and decreases in total quantity, to 40 at years, and 10 at Myr.

4.2 Varying Disc Model Parameters

To quantify the role of radial dependence on dust/gas enhancement results through both the local maximum of grain sizes considered and the time delay in the formation of the largest grains, we present two simulations using alternative grain-growth model assumptions. In the first, the initial particle-size distribution does not vary as a function of disc radius. All particles of all sizes (0.2 m – 2 cm) are instead initiated at evenly distributed throughout the entire disc-gas-mass distribution. This is similar to the radially independent size distribution used as an initial condition in the Youdin & Shu (2002) calculations of dust pile-up via radial drift, except that we retain the simple power-law form of the size distribution with throughout, while Youdin & Shu (2002) considered a two-component distribution representing chondrules and matrix material.

The results of our simulations using a radially independent grain-size distribution are plotted in Figure 8

Figure 8: Enhancement map for a simulation using a radially independent initial grain-size distribution.

and are distinct from those of our fiducial simulations. Due to the large fraction of total dust mass now residing in large grain sizes, most of the dust-mass in the disc rapidly falls inward to small AU. It is thus quickly lost onto the parent star, but not before creating a period of distinct enhancement in the dust-to-gas ratio in the inner disc. The peak in this enhancement occurs at  years at about solar (initial) metallicity. (If we instead assume a maximum grain size of only 2 mm, the peak in the enhancement occurs at  years.) After this peak, disc metallicity declines swiftly, and by years the metallicity across the whole disc is at a fraction its value. While this extreme-case grain-size distribution does produce significant enhancements in the dust-to-gas ratio early in the disc lifetime, the time to reach this enhancement peak is shorter even than the expected time-scale for disc assembly ( years). It is thus hard to envisage a disc formation scenario in which the pile-up seen in this instant-grain-growth-everywhere simulation could be physically realized. Moreover, such a short window supplied by this simulation for planetesimal formation is in direct contradiction with meteoritical evidence implying large-body formation in the Solar System persisted for millions of years.

In the second alternative simulation, we again assume that grain growth is instantaneous, initiating all particles at . However, we now return to the assumption that the largest grains will grow only in the inner, dense regions of the disc. The results depicted in Figure 9

Figure 9: Enhancement map for a simulation with all grains initiated at , but with maximum grain sizes set by radial zone using our cartoon grain-growth model.

use the same simulation output as those of Figure 8 but have been compiled selecting only those grains initiated in the zones appropriate for their size as defined by our cartoon grain-growth model, e.g., 2 mm grains initiated at  AU. 2 This simulation produces results similar to those of our fiducial simulations but with two notable exceptions. First, the small-grain population is depleted relative to the fiducial simulations, producing a somewhat greater depletion of dust relative to the gas at large distances and at late times. In the fiducial simulations, small grains have time to advect and diffuse into the massive outer disc before dust mass at smaller distances is converted into larger sized grains.

Second, there is a peak in the enhancement factor of just over solar metallicity at small radii around  years. This peak corresponds to the terminal infall of the 0.2 mm sized grains, which are carrying a larger portion of the dust mass than in the fiducial simulations. Assuming instantaneous grain growth also allows infall of these larger grains without delay from the further regions of the disc (out to  AU). While the 2 mm and 2 cm grain-size bins are also holding more dust mass than in the fiducial simulation and infall more swiftly than the 0.2 mm grains, their presence is too radially confined to produce an earlier peak in the dust/gas enhancement like that seen in Figure 8. The enhancement peak produced in this simulation is timed more favorably to assist in planetesimal formation than the fully radially independent simulation, particularly under the potentially lower enhancement requirements of streaming-instability–type planetesimal formation. However, its duration is still noticeably shorter than the meteoritical evidence would seem to indicate necessary.

The results of these simulations therefore seem to indicate that faster grain growth leads to a greater probability of dust pile-up in the inner regions of the disc, but that this effect is limited by the mass supplied in large grains initially a large distances from the star. Furthermore, faster growth may still be insufficient to supply an extended period of dust/gas enhancement. For the remainder of our simulations, we return to initiating particles radially and temporally based on our cartoon models of grain growth.

Figure 10

Figure 10: Enhancement maps from simulations run with and 2. We use the fluid-dynamics convention that , while retaining the local scaling for particle diffusivity based on grain size using Equation (10). These simulations were run with output at lower resolution than in the fiducial setup and therefore compiled using six radial-grid-points per radial-bin and averaging 25 time-points for every years time-bin.

plots the enhancement distributions produced varying the Schmidt number, (inversely varying the gas-disc scaling for the particle diffusivity), by a factor of two in either direction. We find the precise value of used to be unimportant to our results. A lower diffusivity leads to greater segregation of the particle-size distribution and emphasizes the discrete particle sizes used for these simulations. However, a real disc should have a mostly continuous particle-size distribution, and therefore the results for the global distribution of solids cannot be said to be qualitatively different from the fiducial model. The higher-diffusivity case allows large grains to remain in the disc slightly longer than usual (all 20 m grains are lost after 4.1 Myr rather than by 3.8 Myr as in the fiducial simulations), but again, this is insufficient to substantially alter the global distribution of solids.

In Figure 11,

Figure 11: Enhancement maps for simulations run in several different disc models, separately varying (horizontal), (vertical), and (diagonal). Fiducial-model map (, AU, ) shown in the center panel. Note that vertical scales are each set to the same years-per-inch to illustrate the varying rates of disk evolution.

we plot the composite enhancement maps for simulations run varying three disc-model parameters: the initial disc mass, ; the initial disc compactness, , and the magnitude of the disc viscosity, . Qualitatively, the patterns of evolving global enhancement are roughly constant across the different disc-model setups. The infall of the smaller grain sizes causes depletion in the dust relative to the gas at the outer edge of the disc, while in the main disc enhancement values remain near unity for somewhere between 40 and 60% of the total disc lifetime. The initially more compact disc shows somewhat greater outer-disc depletion because a larger fraction of grains are processed to larger sizes at early times, while the disc using shows some of the enhanced segregation of our discrete particle sizes typical of the lower-diffusivity simulation in Figure 10 (as, in our transport model, lower means lower disc viscosity and diffusivity). However, the basic shape of dust-to-gas values within the disc as controlled by aerodynamic forces appears fairly universal. The primary result, of these simulations then, is that a smoothly-defined, azimuthally-symmetric, evolving disc model tends to lose dust-mass smoothly inward onto the parent star at a rate that generally keeps pace with the loss of accreting disc-gas mass, and that variety in basic disc properties/parameters does not tend to produce special-case discs with respect to this process.

5 Comparison with Required Enhancement Factors

In this section, we compare our simulated enhancement factors to those required by the Youdin & Shu (2002) precipitation model of planetesimal formation. In the precipitation model, planetesimal formation is accomplished by the gravitational collapse of solids settled out to the disc midplane. As outlined in Sekiya (1998), it assumes that the disc is quiescent, with no global turbulence to stir the particles upward, and that grains will settle out to the limit provided by the Kelvin-Helmholz shearing instability. For midplane dust densities in this settled state above a certain mass threshold, the gravity of the dust sub-disc should dominate and lead to the collapse of solids at the midplane. The criteria for this collapse, given in Youdin & Shu (2002) Equation (15), is


where is the critical surface density of the solids for collapse, is the critical Richardson number defining the degree of dust settling prior to the onset of the Kelvin-Helmholtz instability, and is a correction factor accounting for the self-gravity of the gas. Youdin & Shu (2002) define




For a gas density of , where in a vertically isothermal disc model:


using the conventions for and defined in §2.2.2. One expects that any degree of turbulence beyond the Kelvin-Helmholz effect should stir the particle distribution to greater vertical heights than those assumed in this scenario; certainly, the used in our fiducial disc-evolution models suggests that a significant fraction of the small dust grains should be vertically well-mixed with the gas. However, assuming that grains could settle out to their Kelvin-Helmholtz limit allows us to calculate the absolute minimum enhancement-factor required for precipitation of planetesimals to proceed in that assumed non-turbulent environment.

The precipitation-collapse criterion is separate from the simpler gravitational-collapse criterion using the Toomre- parameter, . From Youdin & Shu (2002),


The minimum enhancement criterion for -collapse is then


which tends to dominate in the outer regions of massive, cold discs.

While the critical Richardson number is commonly set to , recent numerical simulations by Lee et al. (2010) find the onset of Kelvin-Helmholtz instability occurring at


where is the current dust-to-gas volume-density ratio at the disc midplane. can be written as and, if we continue to pursue to limit of maximum settling so that follows from in Equation-set (28), we can write from Equation (30) in terms of the current dust-to-gas enhancement and a value defining the current degree of settling. Absolute maximum settling should therefore occur to (Eq30)(Eq28), or


Below, we compare our simulated dust enhancements to precipitation/collapse requirements using , so that Equations (27) & (29) become




This choice results in a slight decrease (up to a factor of two less) in minimum dust enhancements required in the inner/main disc to produce planetesimal formation relative to the Youdin & Shu (2002) calculations, which use .

Figure 12

Figure 12: Required Enhancements for midplane precipitation using the results of Lee et al. (2010) and Equation (31) to calculate for the discs of Youdin & Shu (2002) and our fiducial-disc model. In our fiducial disc at , the snow-line minimum value is .

plots the required enhancement factors for the disc models discussed in Youdin & Shu (2002) (panel-a) and for our fiducial disc at several times (panel-b). In Table 1,

model (g cm) q T (K) q

H 1700 3/2 280 1/2
A 1700 3/2 170 0.63
B 1700 3/2 100 3/4
Af 1700 1 170 0.63
Af5 8500 1 170 0.63

fiducial 1180 1 280 0.60.5
Table 1: Parameters for disc models used in Youdin & Shu (2002) versus rough equivalents for our fiducial disc near . In general, and . ”H” designates the minimum-Solar-Nebula model of Hayashi-1981 (see Armitage (2010) pp 4–5).

we list the main parameters defining the different Youdin & Shu (2002) disc models, as well as the values roughly corresponding to the inner regions of our fiducial disc at early times. (Youdin & Shu (2002) use simple power-law disc models while our models have an exponential tail-off of surface-density in the outer disc.)

Figure 12 demonstrates several interesting features of these enhancement requirements.

  • First, note that our fiducial model calls for remarkably high in the outer disc, particularly at early times. This is due to the strong dependence of on . When the radial pressure gradient is steep, as near the outwardly expanding disc edge, shearing between the dust and gas is strong, requiring large dust-to-gas ratios to be overcome. This limiting requirement is shared by some other proposed models of planetesimal formation, such as the streaming-instability model, as discussed in §3.

  • Second, while the Youdin & Shu (2002) paper de-emphasizes the roles of the disc gravity, represented by , the time-series for of our fiducial model, particularly in the inner disc, demonstrates the importance of the local disc mass for meeting dust/gas-precipitation requirements. (and therefore ) becomes large, when is small, which may occur for either a low local-disc surface density, or for a locally hot disc.

  • Finally, the curves for both sets of disc models emphasize the importance of the ice lines, particularly the snow line, as places of potential local minimum in metallicity requirements for planetesimal formation. While our particle-transport model does not treat evaporation and vapor-transport effects, the simulations of Ciesla & Cuzzi (2006) do show an additional peak in the dust population just outside the snow line due, at least partially, to local concentration of diffused and recondensed vapor. This also corresponds to a peak in the disc opacity and heating, however, so it is not clear what the degree of advantage the region just outside the snowline lends to planetesimal formation via gravitational precipitation.

In Figure 13,

Figure 13: Measured enhancement values for simulations in the fiducial-disc model at several times and using several particle-size distributions. Compared to the enhancement factor required for precipitation/collapse as given by Equations (32) & (33). The dashed curves plot the simulated distribution using the , -independent grain-size distribution discussed in §4.2, which in the left-top-most panel is plotted at the peak in enhancement, years.

we compare minimum enhancement factors required for planetesimal precipitation to enhancement results from our simulations compiled using three different assumed particle-size distributions. From this figure, it is clear that the size distribution used to compile our simulation results has little impact on the conclusions of the fiducial runs. Radial drift of particles does not produce enhancement sufficient to lead to planetesimal formation via precipitation for this model disc. The exception occurs at the very end of the disc lifetime when photoevaporative clearing of the disc leads to an outward-sweeping pressure-maximum point at large AU where goes to zero and drops to small values.

The first two panels of Figure 13 also include the results of the simulation discussed in §4.2 where all particle sizes were initiated at and distributed throughout the entire disc gas mass, with no radially-dependent grain-growth restrictions. While this simulation yields a distinct peak in the dust-to-gas enhancement early on in the inner disc, it is still less than the minimum enhancement required by the precipitation model.

Finally, while the simulated enhancement distributions vary little between different evolving-disc models (as shown in Figure 11) the enhancement factors required by the precipitation model should vary according to factors such as disc mass. In Figure 14,

Figure 14: values near the beginning ( years) and end (before the onset of photoevaporation) of the disc lifetime for the various evolving-disc models employed. In the fiducial disc ,  AU, and . Note that the onset of photoevaporation and disc clearing is controlled by the disc accretion rate.

we plot the enhancement criteria for precipitation near the beginning and end of the disc lifetime for the various discs used in generating Figure 11. While more massive discs, in particular, produce lower minimum enhancement requirements, these requirements remain generally high relative to the dust-to-gas enhancements produced in our simulations.

6 Comparison to Disc Observations

One of the most general results from our simulated dust distributions is that we predict the outer-most edges of a disc should be depleted in dust relative to the gas. Furthermore, we predict that the dust that does reside in the outer reaches of the disc should be fine-grained in order to remain mixed with such tenuous gas. It is not yet clear from observations whether the outermost regions of a disc tend to be gas rich, though our prediction is consistent with some observations that find smaller disc radii when examining the dust versus the gas signatures of a disc (Hughes et al., 2008). As good observations at the outer edges of a disc require high sensitivity, a full test of outer-edge disc metallicities awaits more sophisticated facilities, such as ALMA.

However, there are observations of dust grain-size distributions within discs around some stars that are inconsistent with the dust-distribution simulations we have presented. Our simulations predict a global loss of mm-sized and larger grains from the disc within half a million years (even in our simulation of §4.2 where grains of all sizes were placed out to the largest distances within the disc). However, observations in the (sub-)mm bands see evidence for mm-sized grains, not just within discs, but out to 100 AU or farther within those discs (Testi et al., 2003; Ricci et al., 2010b). The cartoon grain-growth model we use to constrain the radial extent of simulation particles does not allow mm-sized grains to be initiated beyond 22 AU an any of the disc models we’ve run, and the inward flow from both accretion and headwind drag certainly precludes any simulated mm grains from reaching 100 AU. Furthermore, even in the simulations of Hughes & Armitage (2010) that considered grain transport (using the same fiducial disc model) within an outward flow of disc gas at the midplane, less than 6% of mm-sized grains initiated between 0.5 and 10 AU reached even 25 AU, and then only for a brief period of time before headwind drag forced them inward once more.

That headwind drag tends to bar large grains from the outer disc, in contradiction with observations, is well-known (e.g. Weidenschilling, 1977; Takeuchi & Lin, 2005; Brauer et al., 2007). Brauer et al. (2007) specifically focus on attempting to solve the radial-drift problem by including increasingly sophisticated physics into their transport calculations, including dust settling and collective effects, and angular momentum exchange between a dust sub-layer and the main gas disc. The authors find these effects in general insufficient to slow radial drift of mm-sized grains to velocities compatible with their observed persistence in discs that are generally 1 Myr old or older.

The apparent inability of dust transport models to match observations poses a puzzle, since the basic one-dimensional description of the gas disc evolution appears to be a reasonably good fit to statistical measures of the population such as disc sizes, lifetimes, and accretion rates. It is possible that the axisymmetric gas disc model, while globally descriptive, is qualitatively inadequate as a background model for studies of particle transport. One way in which this could happen is if local or non-axisymmetric structures are integral to accurately describing the large-scale distribution and transport of large dust grains. Otherwise, explaining outer-disc, large-grain populations with current aerodynamic disc theory seems to require such extreme scenarios as a collisional population of outer-disc planetesimals (of unknown origin), or else very massive outer gas discs of very low metallicity (due to the sustained depletion of large grains observed to be in residence) and currently unknown impact on global disc evolution.

7 Discussion and Conclusions

In this paper, we have presented simulations of the global (re)distribution of dust within protoplanetary discs. Our models include evolution of the gas-disc surface-density profile, aerodynamic advection and diffusion of the dust-particle ensemble within the gas, and simulated grain-growth constraints confining the appearance of the largest grain sizes to the inner/main regions of the disc. We find that the global distribution of solids within a 1D evolving disc model follows a fairly common pattern, and does not depend strongly on the values of parameters characterizing the initial disc, or on those describing the turbulence within it. Specifically,

  1. Large, mm or cm-sized particles are lost rapidly onto the the parent star (within 0.5 Myr).

  2. Growth of dust grains up to tens-of-microns size and larger leads to the depletion of dust relative to the gas toward the outer regions of the disc (beyond 100 AU).

  3. The dust-to-gas ratio within the main and inner discs remains near solar (its initial value) for at least the first 40–60% of the disc lifetime.

In general, grains of a given size are confined to regions of the disc interior to at least and experience accelerated infall as () contours near the outer boundary of that dust-population distribution evolve inward.

In Youdin & Shu (2002), the authors present calculations of dust-to-gas enhancement factors resulting from radial drift of dust grains within a static disc, and find a steady increase in dust-to-gas ratios at small disc radii with time. This is in marked contrast to our own simulation results, which do not support the radial-drift hypothesis as a means to increase the metallicity in the inner disc in time for bulk planetesimal formation. The difference between the two calculations is primarily the difference between a purely static- and an evolving-disc scenario. As outlined by the drift/advection timescales presented in Figure 2, a static disc produces ever-slower inward drift as grains approach the central star, providing a natural environment for grain pileup. An evolving disc, however, will advect material steadily inward onto the parent star. Furthermore, as the disc evolves and thins, drift due to headwind drag is accelerated, thus limiting the slowing of inward drift possible, even presuming a zone of the disc were to be free of an accretion/advection flow.

While our simulations do not paint a picture of large-scale increases to disc metallicity within an evolving protoplanetary disc, they do suggest a state of long-term continuity of conditions, perhaps for the first several million years of disc lifetime. This is in keeping with observations that suggest the processing of solids in discs starts early (Kwon et al., 2009; Ricci et al., 2010a) but extends late, as evidenced by the spread in ages measured for chondrules embedded within meteorites (Amelin et al., 2002). Furthermore, that our models predict similar conditions across a range in disc-model parameter space is in keeping with the apparent ubiquity of large-body processing within discs around other stars. The fraction of stars observed to have circumstellar debris discs (with scattered-light levels around a thousand times brighter than that of our own zodiacal dust and Kuiper belt) peaks at around 50% for 20 Myr old B and A stars (Currie, Plavchan & Kenyon, 2008).

We have failed to find circumstances in which global redistribution of solids, due to radial drift, results in substantial enhancements in the dust-to-gas ratio anywhere within an evolving disc. If large enhancements are in fact a prerequisite for planetesimal formation, it may instead be necessary to appeal to local or non-axisymmetric structures to provide them, although the effect of such structures may also be constrained by the headwind-drag produced lack of particles, discussed in §3. Local disc structures that could prove important for dust concentration and/or transport include disc-opacity transitions, such as at the snow-line, spiral arms produced within gravitationally unstable discs (Rice et al., 2006), and vortices within the flow of disc gas. In discs subject to the magnetorotational instability, some simulations suggest that MRI turbulence produces local structure sufficiently long lived to concentrate dust particles (Johansen, Youdin & Klahr, 2009a). Determining which, if any, of these processes is important for particle transport and planetesimal formation will require an understanding not just of the strength of disc turbulence, but also of its detailed physical nature.


We would like to thank Til Birnstiel and an anonymous reviewer for comments which lead to a much clearer and generally improved manuscript. Our work was supported by NASA, under award NNX09AB90G and NNX11AE12G from the Origins of Solar Systems and Astrophysics Theory programs, and by the NSF under award AST-0807471.

Appendix A Calculating the Evolving Disc-temperature Distribution

The model-disc temperature is calculated using a power-law parameterization of the energy-balanced temperature at the disc midplane. To do this, the disc is first evolved once using an iterative solver to calculate the energy-balanced temperature at each grid point and each time-step due to passive and viscous heating and radiative cooling. We then make a power-law fit from the outputted temperature history.

For the energy-balanced temperature, we calculate three components to the disc heating. The flux of energy due to viscous dissipation near the disc midplane is given by


where the viscosity, , is taken from the alpha-parameterization described in §2.1. The flux due to heating via stellar illumination is


where is the Stefan-Boltzman constant, is the passive temperature due to stellar radiation, and the radiative efficiency is assumed to be . For our disc evolution models we use and set the incidence angle to the disc at everywhere. Finally, heating due to background radiation is given by


where K is the background temperature of the cloud environment.

These are used to solve for the midplane temperature of the disc following the approximation derived by Nakamoto & Nakagawa (1994):


where is the temperature at the disc midplane, is the Rosseland mean optical depth, and is the Planck mean optical depth. For the opacity, , we use the piecewise parameterization provided by Bell & Lin (1994) in which depends on both temperature and the mean gas density (for which we use the density at the disc midplane, ). For simplicity we assume , but note that the specific ratio chosen has negligible impact on the energy-balanced temperature calculated. The iterative solver seeks a convergence in according to Equation (37), updating the corresponding values for the midplane density (for a vertically isothermal disc profile), opacity, and disc viscosity at each pass.

The primary focus of this paper is the non-uniform nature of the dust/gas distribution in an evolving disc, an effect which is disregarded for the purposes of the energy-balanced disc-temperature calculation. In this model of disc temperature, the effects of the dust distribution are felt only in the amount of energy deposited via disc accretion. If growth of the dust grains were taken into account, dust opacity should drop for larger grain sizes, thereby allowing for less disc heating via accretion. However, at the time-scales / masses and accretion rates considered in our models, the evolution of disc surface-density contours is reasonably convergent using either a static (cooler, passive disc) or our evolving–energy-balanced temperature distribution. If, on the other hand, massive infall of the dust population led to significant dust enhancement in the inner disc, dust opacities in the inner disc should increase. However, at early times, such infall is comprised solely of large dust grains and so we expect that the effect should be negligible. Infall of the small-grain population occurs only when the disc is too thin for accretional heating to be a significant contributor to the disc temperature.

Next, we fit the temperature history of a given disc model to a parameterization following the form


First and are calculated at each time-step using the intersection at and 200 AU with the outputted, energy-balanced temperatures. Next, a quasi-exponential decay in time is fit separately to both and following the form


with and representative of the initial and final states of the temperature evolution, and and fitting constants chosen roughly by inspection. The parameterization constants used for the temperature evolution of each of the disc models of this paper are given below in Table 2.

() (AU) (K) (K) (yr)
0.01 20 458 280 1.03
-0.59 -0.5 1.09
0.03 5 500 280 0.966
-0.61 -0.5 1.01
0.03 20 500 280 0.917
-0.61 -0.5 0.959
0.03 20 10 500 280 3.510 0.925
-0.61 -0.5 4.010 0.966
0.03 20 500 280 0.907
-0.61 -0.5 0.943
0.03 40 472 280 1.01
-0.6 -0.5 1.07
0.09 20 600 280 1.00
-0.65 -0.5 1.03
Table 2: Disc parameters and power-law temperature fitting constants (Equations (38) & (39) ) for the disc models used in this chapter. The temperature is lower if there is initially much less disc mass at small AU. The time-scaling constants, and give a sense of the relative lifetimes of the various disc models, though the lifetimes are more than an order-of-magnitude longer than either of the scaling-constants. The parameters of the fiducial model are in bold.

These parameterizations are used to calculate disc temperature (and hence, the surface-density evolution of the disc gas) during each of the dust-transport runs (one per grain size per disc model) used to build our dust/gas-evolution simulations.

Appendix B Two Cartoon Grain-growth Models

Within the particle-transport simulations, we use two simple-case grain-growth models to place constraints on when and where grains of a given size will appear within the model disc. These models consider growth of a particle due to differential settling from some height within the disc (called the raindrop model because it mirrors the growth of raindrops falling through the atmosphere), and particle growth due to random, thermal motions at the disc midplane. The basic schematic for these two models is given in Figure 15.

Figure 15: Depictions of grain growth in two cartoon models. The modeled, growing grain is highlighted in red, while the background grains from which it grows are shown in blue with black outlines. Arrows depict the relative magnitude and direction (random in the random-motions model) of grain motions. Growth in the raindrop model is dependent on grain height, , above the disc midplane, while the grain growing in the random-motions model is assumed to exist always at .

The raindrop model is described in Dullemond & Dominik (2005). It considers growth under the extreme condition that only the particle of interest is settling toward the disc midplane and that all other particles remain small and suspended within the disc gas. In this limit


where is the mass of the particle, is the assumed dust-to-gas mass ratio within the disc, set to for all runs of our grain-growth models, is the gas density local to height above the midplane, is the sound speed of the gas, and is the collisional cross-section of the growing particle, in this model set to , where is the particle radius. The settling velocity of the grain, , is defined for the case of Epstein drag. While the final, size of a grain in this model depends upon the height it is initiated, growth saturates for . In our simulations, we therefore initiate grains from this height for the raindrop-growth calculation.

For the random-motions growth model, we retain the assumption that all background grains remain small (m), allowing only the particle-of-interest to grow. Growth in this regime is controlled by the thermal velocity dispersion, , between grains, and is prescribed by


where is the mass of a background particle, and the growth model is run at the disc midplane where densities are highest and growth rates will be fastest. For this model we use the formal definition for the collisional cross-section: .

Note that formally, growth via thermal velocity dispersion is referred to as Brownian grain growth. However, in growth via Brownian motion, the entire small-grain population grows in size concurrently, so that growth is initially fast, but slows quickly as the thermal velocity drops respective to the larger grain sizes. However, in turbulent models of grain-growth, growth via random motions can remain fast as the turbulent velocity dispersion increases for the larger size grains. See, e.g., the static-disc test cases presented in the appendix of Birnstiel, Dullemond & Brauer (2010). By fixing the background grain population to small sizes in our random-motions grain-growth model, we retain some of the faster grain-growth experienced by larger grains in a turbulent environment, albeit, not to the extent that it occurs in formal turbulent-growth simulations and independent of the specific level of turbulence used within our disc-evolution simulations.

In Figure 16,

Figure 16: Grain growth using the two cartoon-grain-growth models. Blue: the raindrop growth model. Red: the random-motions growth model. g cm. Growth here is modeled within the evolving fiducial disc model.

we plot raindrop growth and random-motions growth within our fiducial disc model at three disc radii. Because the inner regions of a disc are hotter and denser, particles there grow fastest and to the largest possible sizes. Particles of a given size are initiated within the particle-transport simulations as soon as either of these grain-growth models grows a particle to that size.


  1. The Cuzzi, Hogan & Shariff (2008) model of planetesimal formation may be an exception as it considers only height-local dust/gas enhancement, possibly achievable by vertical settling alone.
  2. Note: To compare to Figure 1, 9.6 AU is the outer boundary of the radial zone with characteristic radius defined at 6.3 AU.


  1. Alexander, R. D., Armitage, P. J., 2007, MNRAS, 375, 500
  2. Amelin, Y., Krot, A. N., Hutcheon, I. D., Ulyanov, A. A., 2002, Sci, 297, 1678
  3. Armitage, P.J., 2010, Astrophysics of Planet Formation. Cambridge University Press, Cambridge.
  4. Armitage, P. J., 2011, ARA&A, 49, 195
  5. Bai, X.-N., Stone, J. M., 2010a, ApJ, 722, 1437
  6. Bai, X.-N., Stone, J. M., 2010b, ApJ, 722, L220
  7. Balbus, S. A., Papaloizou, J. C. B., 1999, ApJ, 521, 650
  8. Bell, K. R., Lin, D. N. C., 1994, ApJ, 427, 987
  9. Birnstiel, T., Dullemond, C. P., Brauer, F., 2010, A&A, 513, A79
  10. Birnstiel, C., Ormel, C. W., Dullemond, C. P., 2011, A&A, 525, A11
  11. Bottke, W. F., Nesvorný, D., Grimm, R. E., Morbidelli, A., O’Brien, D. P., 2006, Nat, 439, 821
  12. Brauer, F., Dullemond, C. P., Johansen, A., Henning, Th., Klahr, H., Natta, A., 2007, A&A, 469, 1169
  13. Brauer, F., Dullemond, C. P., Henning, Th., 2008, A&A, 480, 859
  14. Carballido, A., Bai, X.-N., Cuzzi, J. N., 2011, MNRAS, 415, 93
  15. Cassen, P., 2001, Meteoritics & Planetary Science, 36, 671
  16. Chiang, E., Youdin, A. N., 2010, Annu. Rev. Earth Planet. Sci., 38, 493
  17. Ciesla, F. F., Cuzzi, J. N., 2006, Icarus, 181, 178
  18. Currie, T., Plavchan, P., Kenyon, S. J., 2008, ApJ, 688, 597
  19. Cuzzi, J. N., Dobrovolskis, A. R., Champney, J. M., 1993, Icarus, 106, 102
  20. Cuzzi, J. N., Davis, S. S., Dobrovolskis, A. R., 2003, Icarus, 166, 385
  21. Cuzzi, J. N., Hogan, R. C., Shariff, K., 2008, ApJ, 687, 1432
  22. Cuzzi, J. N., Hogan, R. C., Bottke, W. F., 2010, Icarus, 208, 518
  23. Dullemond, C. P., Dominik, C., 2005, A&A, 434, 971
  24. Gammie, C. F., 1996, ApJ, 457, 355
  25. Garaud, P., 2007, ApJ, 671, 2091
  26. Goldreich, P., Ward, W. R., 1973, ApJ, 183, 1051
  27. Güttler, C., Blum, J., Zsom, A., Ormel, C. W., Dullemond, C. P., 2010, A&A, 513, A56
  28. Hartmann, L., Calvet, N., Gullbring, E., D’Alessio, P., 1998, ApJ, 495, 385
  29. Hughes, A. L. H., Armitage, P. J., 2010, ApJ, 719, 1633
  30. Hughes, A. M., Wilner, D. J., Qi, C., Hogerheijde, M. R., 2008, ApJ, 678, 1119
  31. Johansen, A., Youdin, A., 2007, ApJ, 662, 627
  32. Johansen, A., Oishi, J. S., Mac Low, M.-M., Klahr, H., Henning, T., Youdin, A., 2007, Nat, 448, 1022
  33. Johansen, A., Youdin, A., Klahr, H., 2009a, ApJ, 697, 1269
  34. Johansen, A., Youdin, A., Mac Low, M.-M., 2009b, ApJ, 704, L75
  35. Kleine, T., Mezger, K., Palme, H., Scherer, E., 2005, 36th Lunar and Planetary Sciences Conference, abstract 1431
  36. Kwon, W., Looney, L. W., Mundy, L. G., Chiang, H.-F., Kemball, A. J., 2009, ApJ, 696, 841
  37. Lee, A. T., Chiang, E., Asay-Davis, X., Barranco, J., 2010, ApJ, 718, 1367
  38. Lodders, K., 2003, ApJ, 591, 1220
  39. Mathis, J. S., Rumpl, W., Nordsieck, K. H., 1977, ApJ, 217, 425
  40. Morbidelli, A., Bottke, W. F., Nesvorný, D., Levison, H. F., 2009, Icarus, 204, 558
  41. Nakamoto, T., Nakagawa, Y., 1994, ApJ, 421, 640
  42. Pringle, J. E., 1981, ARA&A, 19, 137
  43. Ricci, L., Testi, L., Natta, A., Neri, R., Cabrit, S., Herczeg, G. J., 2010a, A&A, 512, A15
  44. Ricci, L., Testi, L., Natta, A., Brooks, K. J., 2010b, A&A, 521, A66
  45. Ricci, L., Mann, R. K., Testi, L., Williams, J. P., Isella, A., Robberto, M., Natta, A., Brooks, K. J., 2011, A&A, 525, A81
  46. Rice, W. K. M., Lodato, G., Pringle, J. E., Armitage, P. J., Bonnell, I. A., 2006, MNRAS, 372, L9
  47. Safronov, V. S., 1972, Evolution of the protoplanetary cloud and formation of the earth and planets. Translation: Israel Program for Scientific Translations, Keter Publishing House, 212 p.
  48. Sekiya, M., 1998, Icarus, 133, 298
  49. Shakura, N. I., Sunyaev, R. A., 1973, A&A, 24, 337
  50. Stepinski, T. F., Valageas, P., 1996, A&A, 309, 301
  51. Takeuchi, T., Lin, D. N. C., 2002, ApJ, 581, 1344
  52. Takeuchi, T., Lin, D.N.C., 2005, ApJ, 623, 482
  53. Testi, L., Natta, A., Shepherd, D. S., Wilner, D. J., 2003, A&A, 403, 323
  54. Throop, H. B., Bally, J., 2005, ApJ, 623, L149
  55. Weidenschilling, S. J., 1977, MNRAS, 180, 57
  56. Weidenschilling, S. J., 1977b, Ap&SS, 51, 153
  57. Wyatt, M. C., 2008, ARA&A, 46, 339
  58. Youdin, A. N., Chiang, E. I., 2004, ApJ, 601, 1109
  59. Youdin, A. N., Lithwick, Y., 2007, Icarus 192, 588
  60. Youdin, A. N., Shu, F. H., 2002, ApJ, 580, 494
  61. Zsom, A., Ormel, C. W., Güttler, C., Blum, J., Dullemond, C. P., 2010, A&A, 513, A57
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description