PARTICLE TRANSPORT IN EVOLVING PROTOPLANETARY DISKS: IMPLICATIONS FOR RESULTS FROM STARDUST

We use the thin-disk model for the viscous evolution of the gas-disk surface density (Pringle 1981). We also include loss of disk gas due to photoevaporation by the central star, so that the surface-density profile evolves according to

$$\begin{equation} \frac{\partial \Sigma _{{g}}}{\partial t} = \frac{3}{R} \frac{\partial }{\partial R} \left[R^{1/2} \frac{\partial }{\partial R} \left(\nu \Sigma _{{g}} R^{1/2}\right) \right] - \dot{\Sigma }_\mathrm{g,wind}\left(R,t\right)\,, \end{equation} \tag{ 1 }$$

where Σ_g is the surface density, t is time, R is the radial distance from the central star, and ν is the local disk viscosity. The $\dot{\Sigma }_{\rm g,wind}$ term is the loss due to photoevaporation by EUV (Hollenbach et al. 1994) or X-rays (Owen et al. 2010). Here, we use the Alexander & Armitage (2007) parameterization of the results of Font et al. (2004), which calculates photoevaporative mass loss due to a 10⁴² s⁻¹ photon flux of EUV from the central star. We do not study here the evolution of particles during photoevaporation (long after the presumed epoch when crystalline material must have been incorporated into comets or their progenitor bodies); we include photoevaporation only because it is essential for providing reasonable model-disk lifetimes.

We evolve the disk numerically using a standard time-explicit-finite-differencing scheme similar to that of Pringle et al. (1986). We use 600 grid cells spaced logarithmically. Innermost and outermost grid points are placed at 0.1 and 15,000 AU, respectively.

We assume that the disk temperature and viscosity are both fixed and vertically uniform. We use a temperature profile of T(R) = T₀R^−1/2 normalized to a temperature of 280 K at 1 AU. This is appropriate for a passive, flared disk after most of the gas infall has occurred and the initial rapid phase of accretion has slowed. Note that very early disk temperatures are expected to be much higher (Kenyon & Hartmann 1987; Bell et al. 1997; Tscharnuter et al. 2009). However, a more precise treatment of disk temperature would require models or assumptions involving, for example, disk chemistry and disk-forming infall rates. For the disk viscosity, we use the α-prescription of Shakura & Sunyaev (1973), which gives

$$\begin{equation} \nu \left(R\right) = \alpha c_\mathrm{s} H_{{g}} \,, \end{equation} \tag{ 2 }$$

where α is a non-dimensional scaling parameter, c_s is the local sound speed, and H_g is the local-disk scaleheight. Having assumed T ∝ R^1/2, Equation (2) leads to ν ∝ R. We adopt α = 10⁻². This value falls within the range of values estimated observationally (Hartmann et al. 1998; Hueso & Guillot 2005; King et al. 2007) and allows dissipation of our model disks in less than 10 Myr. The parameterization in Equation (2) assumes that the source of the viscosity is turbulence within the gas disk. Equation (2) is also important for establishing the turbulent diffusion properties for the particle ensemble as discussed in Section 3.5.

The initial surface density distribution of the gas disk depends on the mass and angular momentum distribution of the parent cloud that collapsed to form the star–disk system (Lin & Pringle 1990; Hueso & Guillot 2005). Our evolving-disk models do not use a steady-disk profile, but rather assume a finite disk, so that at t = 0, the surface density is

$$\begin{equation} \Sigma _{{g}}\left(R, t=0 \right) = \frac{\dot{M}_0}{3 \pi \nu } \left(1 - \sqrt{\frac{R_\mathrm{in}}{R}}\right) \exp {\left(-\frac{R}{R_{d}}\right)} \,, \end{equation} \tag{ 3 }$$

where $\dot{M}_0$ is the t = 0 accretion rate onto the central star, R_in is the inner-disk boundary (set equal to the inner-grid boundary, R_1/2 = 0.099 AU), and R_d is a variable controlling the compactness of the t = 0 profile. In our nominal-disk model, R_d = 20 AU, and the initial disk mass of 0.03 M_☉ dictates $\dot{M}_0 = 1.8 \times 10^{-7} \,M_\odot$ yr⁻¹. The evolution of this model is shown in Figure 1. As the disk evolves, the surface density drops and the disk spreads outward before dissipation via photoevaporation occurs after ∼5.5 Myr. That disks may form in initially compact configurations has been discussed by several authors (Lin & Pringle 1990; Hueso & Guillot 2005; Dullemond et al. 2006) and may be indirectly supported by the steep disk surface densities calculated by Weidenschilling (1977b) for estimates of the minimum-mass solar nebula, and, more recently, by Desch (2007). Depending on the strength of the disk viscosity, a compact initial configuration may also be necessary to process expected disk masses within observed disk lifetimes.

Zoom In Zoom Out Reset image size

In Section 5.4, we vary R_d (retaining M_d,0 = 0.03 M_☉) to consider the effects of varied disk-forming conditions on the transport of particles. To facilitate a comparison with other models, we also consider two models of dust motion within a static disk, including a steady-disk model where R_d = ∞ (retaining $\dot{M}_0 = 1.8 \times 10^{-7} \,M_\odot$ yr⁻¹). As shown in Figure 1, this steady-disk model results in a disk profile that is matched to our t = 0 nominal-disk profile in the inner-disk regions, but maintains a shallower profile and much larger gas surface densities in the outer disk.

The mean radial flow of gas experienced by solid particles within the disk depends on aspects of disk turbulence that cannot yet be predicted from first principles. We maintain this uncertainty explicitly by considering two bounding cases for the gas flow, the "accretion-flow" case and the "midplane-flow" case. For the accretion-flow case, we use the vertically averaged gas velocity, v_acc (Pringle 1981)

$$\begin{equation} v_\mathrm{acc} = - \frac{3}{\Sigma _{{g}} R^{1/2}} \frac{\partial }{\partial R}\left(\nu \Sigma _{{g}} R^{1/2}\right). \end{equation} \tag{ 4 }$$

This velocity can be derived from the 1D fluid equations and results in a gas flow that is predominantly inward. Importantly, however, there is a region of outward-flowing gas in the outer disk where the disk is expanding. The boundary of this region moves outward as the disk evolves and spreads.

In a real 2D disk, the radial flow may vary with height above the disk midplane, and particles that settle toward the midplane may experience a very different radial-gas velocity than that predicted by Equation (4). In particular, 2D (R, z) disks in which the viscosity is given by a generalized α model with no vertical mixing of angular momentum yield an outward flow at the midplane at most disk radii (Urpin 1984; Rozyczka et al. 1994). This flow is highly conducive to the transport of inner-disk grains to the outer disk. For the midplane-flow case, we follow Takeuchi & Lin (2002) and use the azimuthally symmetric Navier–Stokes equation to calculate v_merid, the gas velocity at the midplane:

$$\begin{eqnarray} v_\mathrm{merid} &=& \Bigg[\frac{\partial }{\partial R}\big(R^2 \Omega _{{g}}\big) \Bigg]^{-1} \Bigg[ \frac{1}{R \rho _{{g}}}\frac{\partial }{\partial R} \Bigg(R^3 \nu \rho _{{g}} \frac{\partial \Omega _{{g}}}{\partial R} \Bigg)\nonumber\\ &&+ \frac{R^2 \nu }{\rho _{{g}}}\frac{\partial }{\partial z} \Bigg(\rho _{{g}} \frac{\partial \Omega _{{g}}}{\partial z} \Bigg) - R^2 \frac{\partial \Omega _{{g}}}{\partial t} \Bigg]. \end{eqnarray} \tag{ 5 }$$

Here, ρ_g is the local-gas volume density and z is height above the midplane. Ω_g is the local orbital angular velocity of the gas, determined via force balance in the radial direction (Takeuchi & Lin 2002),

$$\begin{equation} R\Omega _{{g}}^2 - \frac{{\it G M}_\star R}{(R^2 + z^2)^{3/2}} - \frac{1}{\rho _{{g}}}\frac{\partial p}{\partial R} = 0. \end{equation} \tag{ 6 }$$

In this equation, M_⋆ is the mass of the central star (equal to 1 M_☉) and p is the local gas pressure. Note that while Takeuchi & Lin (2002) assumed a power law for the gas distribution and used that to derive a simplified expression for the flow structure of the gas, we have used the input of the disk surface density from our disk-evolution model to solve Equation (5) numerically by assuming vertically uniform temperature and viscosity.

Figure 2 shows the two gas-velocity cases at t = 0 and at t = 1 Myr. Note that the outward expansion of the disk dominates both velocity cases at large disk radii.

Zoom In Zoom Out Reset image size

Our model for particle transport tracks an ensemble of particles within a 1D gas-disk environment and subjects the particles to two non-gravitational forces: (1) aerodynamic drag against the mean gas flow and (2) turbulent diffusion of the ensemble within the disk gas. Coagulation is neglected. Our particle trajectories therefore consist of two components: (1) mean radial motion assuming Epstein drag against the mean-disk flow and (2) a random-walk component due to turbulent diffusion.

These two physical effects are represented numerically by two components to the dust radial velocity: v_srd, the mean radial velocity induced by Epstein drag and v_turb, the turbulent random walk velocity. These velocities are calculated at each evolutionary time step of the gas disk model for each grid point of the model disk; values between grid points are linearly interpolated as necessary. Together, these velocities are used to integrate the dust-grain trajectory as

$$\begin{equation} r_d\left(t+\Delta t\right) = r_d\left(t\right) + \Delta t \times \left(v_{\rm srd} \pm v_\mathrm{turb}\right),\ \end{equation} \tag{ 7 }$$

where r_d is the dust-particle radial position in the disk and Δt is the time step. The time step is chosen locally, to insure accurate integration of particle trajectories, but is capped at the global disk evolution time step. The fidelity of this method in reproducing analytical results for trace-particle diffusion and transport in a static disk is discussed in the Appendix. The relative motions of the dust and gas can vary from tightly coupled to divergent depending on the surface-area-to-mass ratio of the simulation particles and the local density of the disk gas.

We solve for the mean radial velocity of the dust grains, v_srd, subject to Epstein drag against the mean gas flow. In the Epstein-drag regime, the size of the particle is less than the mean free path of a gas molecule and the drag scales with both the particle velocity and the thermal velocity of the gas as (Weidenschilling 1977a)

$$\begin{equation} F_D = m_{d}\frac{C_R}{3}\rho _{{g}}v_\mathrm{therm} (v_{{g}}-v_{{d}}) \,, \end{equation} \tag{ 8 }$$

where C_R and m_d are the surface-area-to-mass ratio and mass of the dust grain, respectively, ρ_g is the local gas density, $v_{\rm therm} = \sqrt{8/\pi } c_{\rm s}$ is the local thermal velocity of the gas, and v_g and v_d are the gas and dust velocities. In our model, we assume that $\rho _{{g}} = \rho _{\rm g,mid} = \Sigma _{{g}} / \sqrt{2\pi } H_{{g}}$, the gas density at the disk midplane. In the accretion-flow case, inputting the midplane density to Equation (8) is an approximation that selects for the tightest possible coupling between the gas and the particle ensemble at a given R; it leads to the most conservative estimate of outward mixing in this gas-flow case. However, the midplane density is the most appropriate choice for use in the midplane-flow case, which presupposes that the entire particle ensemble has settled to the disk midplane.

To solve for the steady mean radial velocity of the dust grains, we must include both the drag force and the orbital motions in a polar coordinate system. (See the Appendix for the explicit force-balance equations of particle motion.) We follow Takeuchi & Lin (2002) and simplify these forces by assuming that the radial acceleration of the grain is zero, dv_r,d/dt ≈ 0, and that the azimuthal acceleration corresponds to a change in the Keplerian velocity with R, dv_ϕ,d/dt ≈ −v_Kv_r,d/2R. The equations for the radial and azimuthal motions of a grain now produce two equations for the steady mean radial velocity, v_srd, as a function of its steady azimuthal velocity, v_sϕd:

$$\begin{eqnarray} v_{\rm srd} &=& v_{r,{{g}}} + \frac{3\left(v_{s\phi d}^2-v_\mathrm{K}^2\right)}{C_R R \rho _{{g}} v_\mathrm{therm}} \nonumber \\ v_{\rm srd} &=& -\, \frac{1}{3}C_R\,R\rho _{{g}}v_\mathrm{therm} \frac{\left(v_{s\phi d}-v_{\phi,{{g}}}\right)}{\left(v_{s\phi d}-v_\mathrm{K}/2\right)} \,, \end{eqnarray} \tag{ 9 }$$

which we solve iteratively for v_srd using the gas-disk bulk flow from our disk model at each radial grid point. In the Appendix, we demonstrate that the trajectories produced using the gridded-v_srd values agree well with direct numeric integrations of the force-balance equations of grain motion for particles in disks with Epstein-drag stopping times ranging from less than 10⁻⁴ Kepler times to greater than 10 Kepler times.

Figure 2 shows how the mean radial velocities of the dust compare to the radial velocity of the gas. Very small grains are well coupled to the gas flow, but the crystalline silicates larger than 10 μm associated with the Stardust collection are large enough that at large distances from the star (tens of AU) and at late times in the disk evolution, the dust motion significantly departs from the gas flow.

To model the turbulent diffusion of the particle ensemble, we add a random walk to the individual grain motion via the addition of a turbulent velocity component, v_turb. This turbulent velocity is set by the dust particle diffusivity, D_p, and by the time step of the random walk motion; it is given by

$$\begin{equation} v_\mathrm{turb} = \pm \sqrt{\frac{2D_\mathrm{p}}{\Delta t}} \,. \end{equation} \tag{ 10 }$$

In our particle-transport model, the time step of the random walk is limited by the time step of the gas-disk-evolution model, but is also allowed to be no greater than the local Kepler time, 1/Ω_K (Ω_K is the Keplerian angular velocity), and no greater than either ΔR₅/v_srd or ΔR₅/v_turb, where ΔR₅ is the radial distance across five grid cells in the direction of particle motion.

We follow Youdin & Lithwick (2007) in calculating the radial diffusion coefficient of the dust grains via

$$\begin{equation} D_\mathrm{p} = D_{{g}} \, \frac{1+4\tau _\mathrm{s}^2}{\left(1+\tau _\mathrm{s}^2\right)^2}, \end{equation} \tag{ 11 }$$

where D_g is the gas diffusion coefficient and τ_s = Ω_Kt_stop. Here, t_stop is the exponential-stopping time. In the Epstein-drag regime, it is given by t_stop = 3/(C_Rρ_gv_therm). The dust and gas diffusion coefficients are equal throughout most of the disk, except in the outermost regions where the gas surface density is very low and the dust diffusion coefficient drops to zero. In our fiducial models, we assume that because the disk viscosity is derived from turbulence within the disk, D_g = ν. In Section 5.3, we vary Sc = ν/D_g to evaluate the impact of this assumption on the potential for the outward mixing of inner-disk grains.

While Equation (10) provides an accurate representation of the particle diffusivity, it is not sufficient for producing accurate diffusion of the particle ensemble in our model gas disk. Among other things, it does not account for the mass distribution of the disk gas. In a real random walk, the time and distance of each step vary according to the local properties of the gas environment. In our model, however, the time step is fixed by other, mostly numerical, considerations. Therefore, to allow our particle ensemble to diffuse according to the actual gas distribution, we must calculate an appropriate probability for whether a given particle in the ensemble will step radially inward (p_in) or outward (p_out).

The proper weighting between p_in and p_out is a function of the imposed time step. Instantaneously (Δt → 0), we know that p_in/p_out = 1 because Δt = 0 represents only a single random walk encounter. However, for Δt > 0, the probability that multiple encounters will occur within that time step becomes nonzero. Furthermore, additional encounters will occur not at the initial position being considered, but either inward or outward of that location depending on the outcome of past encounters. Therefore, the properties of the diffusive medium (the gas) outside the local point of interest play a role in determining p_in/p_out for Δt > 0. The spatial extent of this region of interest is set by how far particles may travel for a given Δt. For the particle ensemble, this is the rms of the displacement of the diffusing particles, or

$$\begin{equation} \Delta R = \sqrt{\left<\left(\Delta x\right)^2\right>} = v_\mathrm{turb}\Delta t = \sqrt{ 2 D_\mathrm{p}\Delta t } \,\,. \end{equation} \tag{ 12 }$$

The weighting for p_in is set by the gas properties between R − ΔR and R; and for p_out by the properties between R and R + ΔR. Asymptotically, the diffusion must yield everywhere a uniform particles(dust)-to-gas ratio. Therefore, after an infinite time, the region with more gas mass must also have proportionally more dust particles. However, because regions of higher diffusivity reach the steady state of uniform concentration more quickly, for a finite time step, proportionally more particles will also mix into a region of higher diffusivity than into one of lower diffusivity. Therefore,

$$\begin{equation} \frac{p_\mathrm{in}}{p_\mathrm{out}} = \frac{\left(M_{{g}}\times D_\mathrm{p}\right)_{\rm in}}{\left(M_{{g}}\times D_\mathrm{p}\right)_{\rm out}} = \frac{\displaystyle \int _{R-\Delta R}^R D_\mathrm{p}\Sigma _{{g}} R\,dR}{\displaystyle \int _R^{R+\Delta R} D_\mathrm{p}\Sigma _{{g}} R\,dR} \,. \end{equation} \tag{ 13 }$$

To evaluate Equation (13) within our model gas disk, we assume uniform, mean values of the gas surface density and particle diffusivity between each grid point. Because we calculate p_in once at each disk-evolution time step, but Δt of the particle motions may vary locally for the individual grains, we calculate p_in at each grid point for both Δt_evolve and $\frac{1}{2}\Delta t_{\rm evolve}$, then interpolate parabolicly for other time-step sizes as needed, using p_in(Δt = 0) = 0.5.

We have verified that this random-walk method does a good job of reproducing the expected diffusion profile of a contaminant within a 1D gas disk. Figure 3 presents a comparison of the distribution of particles using our particle-transport model compared to the expected distribution for an analytical example derived by Clarke & Pringle (1988). A more detailed discussion of the fidelity of our particle-diffusion model may be found in the Appendix.

Zoom In Zoom Out Reset image size

Figure 6 shows the baseline results for particle mixing if the effective gas velocity seen by the particles is that of the accretion-flow case. We plot the fraction of 20 μm particles beyond 25 AU that originate from each of the source regions, together with the same quantity for perfectly coupled particles (i.e., of negligible size). We also show the gas fraction beyond 25 AU. As expected from the reduced models, diffusion is fast enough to mix some particles upstream into the comet-forming region. Indeed, the first particles pass 25 AU in less than 20,000 years. However, loss of particles by accretion onto the star is also rapid, especially close to the inner edge of the disk. After 10⁵ years only ∼30% of all 20 μm particles remain in the disk, and most of those originated in the outer half of our total source distribution (the peak fractions of particles in the comet-forming region are 2.6% from the inner-half zone, but 11% from the outer-half zone). The lifetime of particles that do manage to attain large radii is also limited, first by the fact that the transition radius that divides gas inflow from outflow itself moves out as the disk evolves, and second because declining gas densities in the outer disk eventually preclude retention of 20 μm particles. These effects mean that by t = 10⁶ years, less than 1% of 20 μm particles exist in the disk beyond 25 AU, even if they started in the outer-half source region. By t = 2.1 Myr, all have been lost inward onto the parent star.

Zoom In Zoom Out Reset image size

To quantify the effect of these mixing patterns on the abundance of hot, inner-disk grains in the comet-forming regions, we plot in Figure 7 the normalized concentration, C_N, of source-region particles beyond 25 AU as a function of time. We define the normalized concentration as the number of particles in a region divided by the gas mass in that region, then normalized by the initial number of source-region particles divided by the t = 0 gas mass in that source region:

$$\begin{equation} C_N = \frac{\left(\mathrm{number \, particles}\right)_\mathrm{outer \, disk} / \left(\mathrm{gas \, mass}\right)_\mathrm{outer \, disk}}{\left(\mathrm{number \, particles}\right)_{\mathrm{source},\,t=0} / \left(\mathrm{gas \, mass}\right)_{\mathrm{source},\,t=0}}. \end{equation} \tag{ 14 }$$

With this definition, C_N is like a scaled dust-to-gas ratio, where C_N = 1 means that the ratio of particles-of-interest to gas in the region of interest (the outer disk) is the same as where those particles originated at the time that they formed. We note that while the denominator of Equation (14) is almost constant across our source region (i.e., the initial dust to gas ratio is a fixed value to within 10%), the normalized concentration beyond 25 AU is an average over a broad region. The local values can vary widely, depending upon the parameters of the disk. For the accretion-flow case simulations, the highest local values of the normalized concentration typically lie just outside of 25 AU, with a lower local C_N in the regions hundreds of AU from the parent star.

Zoom In Zoom Out Reset image size

From Figure 7, we see that the peak values of the normalized concentration—and thus the maximum extent to which the outer disk can be contaminated by particles from the hot inner regions—is relatively modest. For 20 μm particles, the peak in C_N from the all source region is approximately 6%. This is actually higher than the value obtained from the outer-half source region, reflecting the fact that although a greater fraction of particles mix outward from this region, there is less mass there to start with. We also observe that particle size strongly influences the lifetime of particles in the outer disk. For well-coupled particles, C_N declines only modestly from its peak at a few hundred thousand years out to several Myr. This reflects the fact that well-coupled particles that survive the initial evolution have thoroughly mixed with most of the entire disk gas mass, so that their subsequent loss is largely in proportion with the disk-gas accretion rate. By contrast, 20 μm particles are lost much faster and only significantly contaminate the outer disk between ∼5 × 10⁴ and 8 × 10⁵ years.

Dramatically higher efficiencies of particle retention are obtained if, instead of experiencing the mean gas flow, the particles settle and experience an outward flow at the disk midplane (recall that this requires little or no vertical mixing of angular momentum, so that the steep radial density gradient at the midplane results in outflow). Results for this midplane-flow case are plotted in Figure 8. In this limit, few 20 μm particles are initially lost inward onto the parent star. Around 98% remain in the disk after 10⁵ years, and around 90% of particles from each source region are transported to beyond 25 AU, with only a slight selection against particles from the innermost regions of the disk. This means that the normalized concentration curves for the different source regions are all nearly the same, only scaled relative to a given source region's innate ability to contaminate the outer disk—the relative fraction of t = 0 disk mass in that source region.

Zoom In Zoom Out Reset image size

Because 20 μm and well-coupled particles are swept to large distances and remain there for a long time, there is a period when the change of normalized concentration beyond 25 AU is a function solely of the loss of disk gas mass via accretion. Without the concurrent loss of accreting particles, the particle-to-gas ratio reaches values even higher than in the original 0.5–10 AU source regions. This is particularly so for the well-coupled particles, which never leave the outer disk once they have entered it. The normalized concentration of 20 μm particles drops to half-maximum at around t ∼ 2.6 Myr, about three times the contamination time window given in the accretion-flow simulations.

Note, however, that our model cannot account for the radial mixing that would occur in the presence of vertical mixing and a vertically stratified gas flow. If the gas flow is rapidly inward far from the midplane and outward near the midplane, then vertical mixing will lead to radial mixing, even in the absence of radial diffusion. Ciesla (2007, 2009) discusses in detail how a given degree of vertical mixing can lead to enhanced mixing of grain populations or to the partial segregation of inward-flowing and outward-flowing populations. Therefore, the normalized concentration for the midplane-flow case merely represents the extreme upper limit on the extent and duration of the contamination of the outer disk by inner-disk material.

The results imply that the efficacy of the fiducial disk model (R_d = 20 AU, Sc = 1) for outward transport of 20 μm grains depends upon both the location (and extent) of the source region for crystalline silicates and upon the radial gas flow experienced by the particles. In the accretion flow case, the peak of contamination of the outer disk (beyond 25 AU) would result in a crystalline-to-amorphous ratio of ∼6%. Attaining this limit requires that at t = 0 all silicate grains in the entire inner disk out to 10 AU were crystalline (and likewise assumes that all beyond that were amorphous). If, instead, the source of crystalline silicates is confined to inward of 2.5 AU then the peak C_N beyond 25 AU is only about 0.2%. Moreover, we find that even these modest degrees of contamination are relatively short-lived. Around 20 μm sized grains (matched to the upper-size end of the Stardust grains) only have a contamination lifetime in the outer disk of ≲1 Myr. Unless the timing of the epoch when grains became assembled into cometary material happened to coincide with the peak of the contamination, the crystalline fraction would be diluted, either by grains grown locally in the cold comet-forming region or by those raining inward from the even-colder outer regions of the disk. Finally, our results for the accretion flow case might be yet further reduced by the fact that some crystalline silicate formation mechanisms may not yield 100% crystalline fractions out to large radii. For example, if the primary source of crystalline silicates is shocks out to ∼10 AU, but a sizable fraction of silicates in that region remain amorphous, then that cuts into the potential maximum in the crystalline-to-amorphous ratio produced in the outer disk.

In the midplane-flow case, on the other hand, substantial outward transport can occur almost irrespective of the location of the source region. If there is an outward flow of gas at the disk midplane, then grains that are sufficiently well settled will readily reach large distances (Ciesla 2009). An outward-midplane flow also allows material that has been transported outward to persist there for a longer period of time. Our midplane-flow simulations produce an upper-limit maximum C_N in the outer disk of 30% for the inner-quarter source region. Even if vertical mixing means that only a fraction of the maximum can be achieved, it is still plausible to believe that the CAI-like grains captured during the Stardust mission may have reached the outer disk via mixing within the disk gas alone. This scenario may also be consistent with the observations of Watson et al. (2009), who find higher crystalline-to-amorphous silicate mass ratios in disks where the grains appear more settled toward the disk midplane, and of Olofsson et al. (2009), who find that possibly a large fraction of crystalline silicates exist in the outer regions (≲10 AU versus ≲1 AU) of some disks.

In Section 4, we demonstrated that non-negligible particle size bars particles from the outermost regions of the disk because at large distances, the gas is so tenuous that the particles decouple from the gas motions, and inward-pointing head-wind drag dominates. Figure 9 plots the dust-mean-radial velocities for several different particle sizes in the nominal-disk model at t = 0 for the accretion-flow case. (In the midplane-flow case, the curves are similar, only shifted toward outward velocities as in Figure 2.) Though turbulent diffusion can send particles upstream of their average radial-velocity flow, grains are effectively barred from the regions of the disk beyond where their average-radial velocity falls to large, inward-pointing values. It is clear from Figure 9 that the range within the disk that dust grains may occupy depends strongly on the particle size, and that the mixing achievable out to the comet-forming regions will become sharply limited as one approaches millimeter particle sizes. Also, where this cutoff in grain size occurs changes, not only in time as the disk thins, but also from disk to disk depending on the total mass and mass distribution of each system.

Zoom In Zoom Out Reset image size

We parameterize the outward mixing of grains as a function of size by plotting in Figure 10 the peak values of the normalized concentration of grains beyond 25 AU as a function of grain size. The error bars show statistical errors based on the number of particles beyond 25 AU in each run. We plot the maximum values of C_N from each of the different source regions, as well as the upper-limit values provided by the midplane-flow simulations (dashed lines with black markers). While the magnitude of outward mixing varies by orders of magnitude across these populations, the trends with grain size are the same across all sets: as seen previously, the highest normalized concentrations of inner-disk particles in the outer disk occur for the smallest grain sizes, up to a few tens of microns, and the peak-C_N values drop off sharply at around millimeter grain sizes. The results of Ciesla (2007, 2009) for a vertically stratified disk flow suggest that the loss of larger grain sizes may be partially mitigated by settling toward the midplane, thereby increasing the relative outward transport and retention of these particles. Still, the concurrent drop in our midplane-flow upper-limit concentrations suggests that this can only hold back the loss of large grains from the outer disk to a point. For 2 mm sized particles, for example, even the midplane-flow runs yield a peak normalized concentration beyond 25 AU for all simulation particles of only ∼3.4%.

Zoom In Zoom Out Reset image size

The other important variable to consider is the timing of these mixing events. In Figure 11, we depict, as a function of particle size, the time frame over which the normalized concentration of inner-disk particles in the comet-forming region is at a value of half-maximum or greater. For each source region (for the accretion-flow case simulations), we plot the latest time of half-maximum concentration prior to the peak and the earliest time of half-maximum after the peak. (Due to limited particle statistics, some of these quantities are quite noisy.) We also include the boundary of post-peak half-maximum concentration in the midplane-flow simulations to illustrate when headwind-drag inspiral will overwhelm other physical effects even in the most extreme gas-flow scenario.

Zoom In Zoom Out Reset image size

Figure 11 demonstrates that the smallest inner-disk particles remain mixed into the outer disk for millions of years and that even inner-disk grains as large as 20 μm may be noticeably present at large distances for ∼1–2 Myr. However, while the mixing of all grains outward is relatively rapid in these simulations, not only do the larger particles experience relatively less outward mixing, but that mixing is similarly short lived. Even in the midplane-flow simulations, post-peak half-maximum concentration is reached by t ∼ 6.2 × 10⁵ years for 0.2 mm particles and as soon as t ∼ 1.1 × 10⁵ years for 2 mm particles. The fact that larger inner-disk particles are short-lived in the outer disk is demonstrated more bluntly in Figure 12, where we plot the normalized concentration as a function of grain size at t = 1 Myr. While the mixing effects are relatively uniform for small particles at this time, no particles 0.2 mm and larger remain in the outer disk in the accretion-flow simulation, and the upper-limits of the midplane-flow simulations suggest that even in that extreme those large-grain populations are in rapid decline, if not already disappeared.

Zoom In Zoom Out Reset image size

These results imply that it is the upper envelope of particle sizes recovered by Stardust that provides the most stringent constraints on disk evolution and particle transport. In our simulations, particles at least as large as 2 μm appear to behave as though well mixed with the disk gas. Therefore, the presence of high-temperature particles in this size range within a Jupiter-family comet can largely be explained by outward mixing of disk material aided by some disk expansion from an initially compact state. However, the particles recovered by Stardust include crystalline silicates as large as 20 μm. While the outward transport of 20 μm particles to the comet-forming region does occur in our models, these particles do not remain well coupled to the disk gas indefinitely. This suggests that either (1) mixing of particles and planetesimal formation in the solar nebula occurred in a disk massive enough for 20 μm particles to remain well coupled to the gas motions, or that (2) the formation of comets or cometesimals occurred within the first million years or so of disk evolution. This second constraint may be in agreement with Stardust results. Early studies found no aqueous minerals in the Stardust materials, suggesting that the outward transport and incorporation of these materials into the 81P/Wild 2 cometesimals occurred very early, before the primary accretion and fragmentation of the chondritic-asteroid parent bodies (Wooden 2008). However, the possible identification of igneous materials in the Stardust samples may contradict this scenario or else put interesting time-line constraints on the formation time scales of planetesimals and cometesimals across the range of solar system radii (Stodolna et al. 2010; Joswiak et al. 2010). Another constraint linked closely with time is the large-end grain size for which our model predicts measurable outward mixing. Not only do the peak (and midplane-flow–upper-limit peak) values of C_N beyond 25 AU drop off markedly for grains a few millimeters in size, but the time frame for that outward transport is also restricted (to barely more than 10⁵ years for 2 mm sized grains), so that comet formation would have to be quite rapid to capture such large inner-disk grains if they did appear in the comet-forming regions.

The transport of different-sized particles is also of interest for models that include the effects of grain growth. Some observational studies of disks suggest a link between disk crystallinity and grain growth (van Boekel et al. 2005; Olofsson et al. 2009), with characteristic crystalline grain sizes of a few microns (distinctly larger than typical ISM grains). At the very least, grain growth and crystallization likely occur on similar time scales. The results presented here suggest that as dust grains grow, they should become less likely to enter and more likely to leave the outer regions of the disk. Large grains become mostly confined to small disk radii where the disk gas is still dense enough to support them. We know that growth at centimeter and meter scales poses a constraint to disk models and planetesimal formation because particles of this size are expected to fall inward toward the parent star very rapidly, depriving the disk of these solids. However, growth to millimeter scales may also place constraints on cometesimal formation if the outer regions of the disk cannot support particles of even these sizes. At the very least, it could restrict the composition of bodies formed in the outer disk, if most of the particles available are those falling inward from large AU.

Next, we explore how the outward mixing of inner-disk grains depends on the Schmidt number, which is the ratio of the disk viscosity to the gaseous diffusivity (so that a lower Schmidt number means more relative diffusivity.) In Pavlyuchenkov & Dullemond (2007), the authors argue for an Sc lower-limit of 1/3. Therefore, we have run simulations of particle transport varying the Schmidt number by a factor of 3, and in Figure 13 we plot the values of maximum normalized concentration beyond 25 AU for four particle sizes and two source regions as a function of Sc. The primary effect shown in Figure 13 is that a higher diffusivity can lead to a substantially higher degree of outward transport in the general accretion-flow case of our particle-transport simulations. Here we see variations of an order-of-magnitude or more in the peak outer-disk normalized concentrations between the highest (Sc = 1/3) and lowest (Sc = 3) diffusivity simulations across the range of particle sizes considered. This is less than the variation predicted in the model scenario of Pavlyuchenkov & Dullemond (2007; σ_{Sc = 1/3}/σ_{Sc = 3} ∝ (R/R_source)⁴), but is still quite a substantial effect considering that our normalized concentrations are reported for the whole of the outer disk and that our pool of contaminants is limited to those present in the source regions at t = 0 only.

Zoom In Zoom Out Reset image size

More generally though, the Schmidt number controls the relative importance of advection versus diffusion of particles. When the diffusivity is low, advection dominates, so that more grains reach the outer disk if the bulk flow is outward, and fewer do so if their bulk flow is inward (including the case of 2 mm particles in the outward-flowing-gas midplane-flow case whose radial advection is dominated by headwind drag). Conversely, high diffusivity pushes the system in the direction of a flow-independent state where diffusion, both within the disk gas distribution and inward onto the parent star, is the dominant term. The diffusivity can also have an important effect on outer-disk contamination from the innermost source regions. For 20 μm sized particles from the inner-quarter source region, Sc = 1/3 simulations produce a peak normalized concentration in the comet-forming region of 1.7%, up from 0.2% in the baseline Sc = 1 simulations. In the low-diffusivity simulations, inner-quarter source grains of all sizes are completely absent from the comet-forming region unless directly advected there by outward-flowing gas in the midplane-flow simulations.

Finally, to some degree, the diffusivity also affects the timing of outward mixing. This is plotted for 20 μm sized particles in Figure 14. The higher diffusivity simulations allow higher normalized concentrations in the outer disk to occur sooner and to last longer (∼4 × 10⁴–1.1 × 10⁶ years half-max to half-max), whereas for lower diffusivity the outer-disk normalized concentration peaks later and more briefly (∼9 × 10⁴–5.5 × 10⁵ years).

Zoom In Zoom Out Reset image size

Our results for the variation in transport efficiency with Schmidt number suggest that two distinctly different disk models could be consistent with observations. One possibility is that the diffusivity is high, as suggested by Pavlyuchenkov & Dullemond (2007). In this case, it seems plausible that a significant fraction of grains, formed in high temperature regions close to the star, can end up in the comet-forming region, even if the disk is highly turbulent and no significant settling occurs. The recovery of CAI-type grains by Stardust—whose condensation requires especially high temperatures—poses the strongest constraint on this scenario, and we have not demonstrated explicitly that it is possible. However, given that the early disk was certainly hotter than the disk which we have modeled, it is plausible that enough material from the innermost disk could reach large enough radii. Hot outflows near the disk midplane, which Tscharnuter et al. (2009) find very early in the star/disk formation process, would assist outward transport from the innermost disk.

A qualitatively distinct disk model with low diffusivity is also conceivable. Our results suggest that if Sc is significantly larger than unity, there is negligible transport of crystalline silicates subject to the accretion flow to the region where Jupiter family comets form. This is true even if silicates can form at radii as large as 10 AU. A low diffusivity disk, however, could potentially be favorable to the establishment of an outward midplane flow, which is able to move settled particles outward efficiently. Ciesla (2009) showed that a transport scenario dominated by high-altitude–inward and low-altitude–outward advection can lead to segregated grain populations with different processing histories. This is possibly compatible with disk observations of Olofsson et al. (2009) that suggest some disks are more crystalline at larger distances.

Finally, we assess the role of the disk mass distribution and its evolution on outward transport. First, we contrast the results of runs assuming a static disk with those of the evolving disk, and then move on to exploring the effects of the initial compactness of the disk in the evolving-disk model. Of our two static-disk models, the first uses the t = 0 profile of the nominal-disk held static (the static0-disk model). The second is the steady-disk limit of the t = 0 nominal-disk profile, where the disk surface density follows Equation (3), with R_d = ∞ (the steady-disk model).

Results from the two static-disk models are shown in Figure 15, which plots the normalized concentration of 20 μm particles beyond 25 AU with time for two different source regions. The mass per AU of the gas disk does not drop off in the outer disk for the steady-disk model but is instead continuous out to the edge of our simulation space. Therefore, the normalized concentrations beyond 25 AU for this disk model are calculated using the disk gas mass only out to 100 AU (for a total of 0-to-100 AU disk mass of 0.162 M_☉; 0.03 M_☉ of gas is contained within the first 21 AU in the steady-disk model). Figure 15 includes the results for the evolving nominal-disk model for comparison. The most obvious result from Figure 15 is that the two static-disk models, together with the nominal evolving-disk model, form a hierarchy in outward mixing that is based on disk structure. In the static0-disk model, the region of outward-flowing accretion flow remains fixed relatively close to the central star and those outward velocities are fairly rapid. Therefore, more than twice the peak fraction of nominal-disk particles move to the outer disk (and stay there, as per the static-disk-model tendencies seen in Section 4). The surface-density structure in the steady-disk model means that in the accretion-flow case the gas radial-velocity is everywhere inward, and hence the number of 20 μm particles beyond 25 AU peaks at less than 40% of the value for the corresponding nominal evolving-disk simulations, and those particles eventually all fall back inward onto the parent star—and more rapidly than in the evolving-disk simulations. Roughly 23% remain in the disk after 10⁵ years (when 30% remain in the nominal-disk accretion-flow simulations).

Zoom In Zoom Out Reset image size

However, in terms of the normalized concentration of inner-disk particles in the outer disk, the disk-gas mass distribution also really sets the steady-disk simulations apart from the rest. With so much more mass in the outer disk than the other two models, the steady-disk model is that much less intrinsically capable of contaminating the outer disk with inner-disk grains than are the disks with a more compact disk mass distribution. Finally, for the upper limits placed by the midplane-flow simulations both static-disk models send more than 90% of the simulation particles from all source regions to beyond 25 AU; the outward advection is only more rapid in the static0-disk model. Therefore, the upper bounds on the mixing in these models are set only by their relative source region and comet-forming region disk-gas mass ratios.

Although some trends are the same for static as compared to evolving-disk models, the results imply that the evolution of the disk—particularly at early times—has a substantial impact on particle transport. Steady-disk models tend to underestimate the degree of outward particle transport possible in protoplanetary disks in the early stages of disk evolution, when the most processing of high-temperature materials occurs. The outward transport of particles is substantially more efficient at early times, consistent with the results of Bockelée-Morvan et al. (2002) and Ciesla (2009), even here assuming a fixed diffusivity. How strong the effects of evolution are will depend, of course, on the details of the initial conditions for the disk. To examine this dependence, we vary the initial compactness of our evolving-disk model (parameterized by R_d in Equation (3)). We consider R_d = 5, 10, and 40 AU, in contrast to R_d = 20 AU in the nominal-disk model. We retain 0.03 M_☉ as the starting mass of the disk, so that the lifetimes of the model disks vary between 3.6 and 6.9 Myr, and the t = 0 mass accretion rates between 8.3 × 10⁻⁷ and 8.6 × 10⁻⁸ M_☉ yr⁻¹. We do not consider possible variations in heating for disks forming in more or less compact configurations, retaining the static temperature profile described in Section 3.1.⁵

Reducing the compactness of the model disk has two consequences: the disk initially expands more rapidly and it does so from an initial state where a greater fraction of the mass is at small radii (and thus hot and potentially able to form crystalline material). The combined magnitude of these effects is shown in Figure 16, where we plot the maximum normalized concentration of 20 μm particles beyond 25 AU as a function of R_d. Variations of approximately one order of magnitude in R_d result in nearly two orders of magnitude difference in the peak concentrations for each source region. For the smallest disk—with R_d = 5 AU—the normalized concentration of all simulation particles peaks at ∼ 68%, more than an order of magnitude in excess of the value for the nominal simulations. A substantial part of this increase is due to the larger reservoir of mass that a compact disk has at small radii, which can potentially contaminate the (proportionally less massive) outer disk. That this mass-distribution effect is important is made clear by the fact that even the upper-limit concentrations from the midplane-flow simulations show a strong (order-of-magnitude) trend with R_d.

Zoom In Zoom Out Reset image size

The fact that the outer disks of the initially more-compact-disk models are less massive (and more tenuous) does place some constraints on the outward mixing of large grains in our particle-transport simulations. In Figure 17, we plot the peak normalized concentration of 2 mm particles beyond 25 AU as a function of R_d. All of the disks, of course, have less outward transport of these larger grains, but in our simulations, the peak in the fraction of particles sent to beyond 25 AU occurs for R_d = 10 AU, which sends up to 0.25% of all 2 mm simulation particles beyond 25 AU (4% in the midplane-flow case) compared to less than 0.1% for R_d = 5 AU (less than 2% in the midplane-flow case). Because of the larger source region mass for R_d = 5 AU, this drop in fractional outward transport constitutes a leveling off in the maximum-achievable normalized concentrations of inner-disk particles in the outer disk.

Zoom In Zoom Out Reset image size

Comparing these results to observations of crystalline silicates in other disks allows us to construct a possible scenario for the formation and transport of high-temperature minerals in those disks. In the Watson et al. (2009) survey, the disks observed were around T Tauri stars ∼1–2 Myr old and showed a slight correlation between the crystalline-silicate mass fraction and the measured accretion rate onto the star. As shown in Figure 18 for our model disks, t = 1 Myr corresponds to a time when those disks that were initially most compact (with high t = 0 accretion rates) now have lower accretion rates than the initially more-extended disks. Therefore, small grains that formed at small AU near t = 0 should be more broadly distributed in disks that now have lower accretion rates. However, lower accretion rates also mean disks that are more tenuous. The initially most compact disks lose mass the fastest, causing their outward transport efficiency to drop rapidly with time. Therefore, grains that formed at small AU after t ∼ 2 × 10⁵ yr may be better distributed in the initially more-extended disks, which now have higher accretion rates.

Zoom In Zoom Out Reset image size

These properties could potentially explain the opposing trends between crystallinity and observed accretion rate reported by Watson et al. (2009). Pyroxene is favored over olivine when the minerals are formed by condensation (Gail 2004; Wooden et al. 2005). Therefore, a possible scenario is that pyroxene formed primarily at very early times when the disks were hottest and became most broadly distributed throughout the initially most-compact disks. Most of the olivine, on the other hand, may have formed over a longer period of time. It might, therefore, have become more broadly distributed within disks that retained significant surface densities for longer times (the initially more-extended disks). This hypothesis does not, however, explain why the Watson et al. (2009) survey found no correlation between crystalline mass fraction and disk mass.

Our particle-transport simulations use an approximation technique to calculate particle trajectories. Here we compare the trajectories produced using this technique to precision numerical integrations of those trajectories. The equations of motions for a particle in a laminar disk-gas flow in the Epstein-drag regime are (Takeuchi & Lin 2002)

$$\begin{equation} \frac{d v_{r,{d}}}{dt} = \frac{v_{\phi,{d}}^2-v_\mathrm{K}^2}{R} - \frac{1}{3}C_R\,\rho _{{g}} v_{\mathrm{therm}} (v_{r,{d}}-v_{r,{{g}}}) \,, \end{equation} \tag{ A1 }$$

$$\begin{equation} \frac{dv_{\phi,{d}}}{dt} = -\frac{v_{\phi,{d}}v_{r,{d}}}{R} - \frac{1}{3}C_R\,\rho _{{g}} v_\mathrm{therm} (v_{\phi,{d}}-v_{\phi,{{g}}}) \,, \end{equation} \tag{ A2 }$$

where v_r,d and v_ϕ,d are the radial and azimuthal velocities of the dust particle, v_r,g and v_ϕ,g are the local radial and azimuthal velocities of the disk gas, v_K is the local Keplerian velocity, C_R is the particle surface-area-to-mass ratio, and ρ_g and v_therm are the local-gas volume density and thermal velocity, respectively.

In the case of laminar-gas flow with no turbulent forcing of the particle motions, we can numerically solve for each particle trajectory directly. To do this, we use a Runge–Kutta-style integrator, where Equations (A1) and (A2) are coupled with

$$\begin{equation} \frac{dr_{{d}}}{dt} = v_{r,{d}} \end{equation} \tag{ A3 }$$

to define the particle motion radially within the disk. The integrator uses step-size adjustment to reach a relative accuracy of 10⁻⁵–10⁻⁶ in R_d, v_r,d, and v_ϕ,d. Note, however, that the integrator becomes numerically expensive when the individual terms in Equations (A1) and (A2) are large but sum to values that are small. Therefore, if the initially computed acceleration (dv_r,d/dt or dv_ϕ,d/dt or both) is very small (<10⁻¹⁴ cm s⁻²), the integrator sets the small acceleration equal to zero. Also, if v_ϕ,d ∼ v_K (10⁻¹⁰ relative difference) then the v²_ϕ,d − v²_K term in Equation (A1) is set equal to zero, but only if $\frac{1}{3}C_R\rho _{{g}}v_{\rm therm} \ge \frac{v_{\rm K}}{R}$. If v_ϕ,d ∼ v_K but the second condition is not met, then the v²_ϕ,d − v²_K term is still important and the integrator instead calculates the particle trajectory to the lower relative accuracy of 10⁻⁵.

Most of the comparison simulations presented below are run using the midplane-flow-velocity case, as defined in Section 3.2 of this paper. We begin by presenting particle trajectories calculated within an initially tenuous disk (R_d = 500 AU, Σ_g(1 AU) ∼ 1 g cm⁻² at t = 0) that is photoevaporatively cleared within ≲10⁵ years. Radial particle trajectories in such disks can be complex and so provide good tests of our particle-transport model. Figure A1 plots, for both numeric integrations and particle-transport simulations, the trajectories and radial velocities of 0.2 mm sized particles initiated at three different radial locations within this tenuous disk. Here, such large particles are only loosely coupled to the gas flow, with exponential-stopping times (T_s) at t = 0 of a few to greater than 10 Kepler times. Nevertheless, there is good agreement between our particle-transport simulations and the more precise numeric integrations.

Zoom In Zoom Out Reset image size

However, our particle-transport simulations do lack certain details; as shown in Figure A2, they do not account for orbital eccentricity. Thus, the particle-transport simulations are limited in their depiction of complex gas–particle interactions. In Figure A2, these interactions pertain to a particle disengaging from and outward-sweeping photoevaporation front. However, the results of the particle-transport simulations presented in this paper are not sensitive to such a high level of detail. Instead, we are most interested in the broad effects of disk mass and particle size. For the tenuous disk, we ran simulations for particles of sizes 0.2 μm, 2 μm, 20 μm, and 0.2 mm from 20 starting locations, 0.5–500 AU. The exponential-stopping times in the simulations ranged from less than 10⁻⁴ to greater than 10 Kepler times, and in general, the relative radial positions of the simulated trajectories agreed with the numeric integrations to within 2% or better.

Zoom In Zoom Out Reset image size

Where our particle-transport simulations deviate the most from the numeric integrations is for particles experiencing a long period of steady infall. This is true regardless of how well coupled the particles are to the gas motions and is an effect of the particle-transport model repeatedly overshooting the estimate of the infalling particle trajectory. Figure A3(a) plots trajectory comparisons for particles initiated at 500 AU in a massive, but initially compact disk (M_D = 0.05 M_☉ and R_d = 10 AU). In this scenario, the particles begin very decoupled from the gas (T_s > 10³ Kepler times), and after their initial infall toward the main disk they become moderately well coupled (T_s ∼ 10⁻² and smaller). In Figure A3(b), we plot the relative difference in the particle radial positions between the particle-transport model trajectories and the numerical integrations. Here we can see that for each infall, these two trajectory calculations diverge as the simulated trajectories outpace those of the numeric integrations. Therefore, exaggerated, rapid infall of particles may be a drawback of our particle-transport model. However, inspection of all of our simulated trajectories suggests that this exaggerated infall should have a negligible effect on our general understanding of the timing and trends of particle transport. As a check, we have also run comparison integrations using the accretion-flow case in our nominal-disk model (Section 5.1). We ran four particle sizes (0.2 μm–0.2 mm) initiated at 15 AU and found positional differences of less than 2% for all but the very final stages of infall.

Zoom In Zoom Out Reset image size

Therefore, in the absence of turbulence, the simulated trajectories produced by our particle-transport model appear to be a sufficiently accurate representation of radial particle transport in an evolving disk under the influence of gas drag, both from the pressure-supported–azimuthal and the radial gas flow.

Here we demonstrate the fidelity of the random-walk method in our particle-transport model in reproducing accurate diffusion of the particle ensemble. Using the analytical solutions of Clarke & Pringle (1988) we compare the time-dependent particle distributions from our particle-transport simulations to the analytically expected contaminant distributions and discuss the effects of our choice of simulation time step on the results. All solutions and simulations presented below use static-steady-disk–surface-density profiles.

In Clarke & Pringle (1988), the authors solve for the time-dependent concentration, C, of a contaminant initiated in a ring at a specified radius: C(t = 0) = C₀δ(R − R₀), where C is the ratio of the local contaminant mass over the local disk mass, R is the distance from the central star, and R₀ is the initial position of the contaminant. The disk that they consider has a surface-density profile given by Σ = Σ₀R^−a, where Σ₀ and a are constants; it is also a steady disk, so that for disk viscosity, ν = ν₀R^b (b and ν₀ constants), b = a, and the gas-accretion velocity may be written as v_R = v₀R^b−1 where v₀ = −3ν₀/2.

Assuming azimuthal symmetry and that the gas and contaminant are vertically well mixed, the diffusion equation in polar coordinates may be written as (Gail 2001)

$$\begin{eqnarray} \frac{\partial }{\partial t}\left(\Sigma C\right) &&+ \frac{1}{R}\frac{\partial }{\partial R}\left(R\Sigma C v_R\right)\nonumber\\ &&- \frac{1}{R}\frac{\partial }{\partial R} \left[R D_{ig}\Sigma \left(\frac{\partial C}{\partial R}\right)\right] = P_i\left(R,t\right) \,, \end{eqnarray} \tag{ A4 }$$

where D_ig is the coefficient of diffusion of the contaminant within the gas and P_i corresponds to the rate of production of the contaminant per radial increment within the disk. (Clarke & Pringle 1988, and the setup of our simulations) assume zero production of the contaminant, P_i(R, t) = 0, and that the diffusivity scales with the disk viscosity, D_ig = ζν, where ζ is a constant. For these particle-transport simulations using a static disk profile, Equation (A4) is further constrained by ∂Σ/∂t = 0 and so v_R refers specifically to the set dust-mean-radial velocity. Using the profiles for the disk surface density and viscosity above, the full diffusion equation applicable to these simulations is

$$\begin{eqnarray} &&\frac{1}{\nu _0}\frac{\partial C}{\partial t} + \frac{v_0}{\nu _0}R^{b-1}\frac{\partial C}{\partial R} + \left(b-a\right)\frac{v_0}{\nu _0}R^{b-2}C\nonumber\\ &&\quad = R^{a-1}\frac{\partial }{\partial R} \left[\zeta R^{1+b-a}\left(\frac{\partial C}{\partial R}\right)\right] \,, \end{eqnarray} \tag{ A5 }$$

which simplifies to

$$\begin{equation} \frac{1}{\nu _0}\frac{\partial C}{\partial t} + \frac{v_0}{\nu _0}R^{a-1}\frac{\partial C}{\partial R} = R^{a-1}\frac{\partial }{\partial R} \left[\zeta R\left(\frac{\partial C}{\partial R}\right)\right] \, \end{equation} \tag{ A6 }$$

in the steady-disk (a = b) scenario of Clarke & Pringle (1988). The scenario of Clarke & Pringle (1988) assumes an evolving (accreting) gas disk, but because it is steady (has an unchanging relative-surface-density profile) and the diffusion equations do not depend on the actual value of the gas surface density, Equation (A6) fits both that scenario and our static steady-disk scenario.

The general solution given in Clarke & Pringle (1988) to Equation (A6) using the initial contaminant distribution given above is

$$\begin{eqnarray} C\left(R,t\right) &=& \frac{|q|C_0 R_0^{1-a}}{2 \zeta \nu _0 t} \left(\frac{R}{R_0}\right)^{v_0/2\zeta \nu _0} I_\beta \left(\frac{q^2 R_0^{1/q}R^{1/q}}{2\zeta \nu _0 t}\right)\nonumber\\ &&\times \exp \left[-q^2\frac{\left(R_0^{2-a}+R^{2-a}\right)}{4\zeta \nu _0 t} \right] \,, \nonumber \\ q &=& \frac{2}{\left(2-a\right)} \,, \nonumber \\ \beta &=& \frac{|v_0/\zeta \nu _0|}{|2-a|} \,, \end{eqnarray} \tag{ A7 }$$

where I_β(x) is the modified Bessel function of the first kind of order β.

However, Equation (A7) is not valid for the case of a = b = 2, for which Clarke & Pringle (1988) provide the solution

$$\begin{equation} C\left(R,t\right) = \frac{C_0}{R_0}\left(4\pi \zeta \nu _0 t\right)^{-1/2} \exp \left(-\frac{\left[\ln \left(R/R_0\right)- v_0 t\right]^2}{4\zeta \nu _0 t}\right) \,. \end{equation} \tag{ A8 }$$

Note, however, that Equation (A8) does not match Equation (3.1.3) in Clarke & Pringle (1988) because of a typo in that paper, which accidentally presents the solution for twice the background velocity, v_R → 2v_R. Also, when a = b = 2, the natural space to solve the equation is logarithmic space. This happens to be the space used for the width of our model-disk grid cells and may be related to the only odd behavior that we find in our particle transport simulations (see below).

From our particle-transport simulations, we output the number of particles per radial grid space. Therefore, to compare the solutions provided by Equations (A7) and (A8) to our simulations, we must first convert the analytic C(R, t) into an expected fractional-mass distribution. The concentration can be expressed as the ratio of the contaminant-to-gas surface densities, C = σ/Σ, where σ is the surface density of the contaminant. Because these test cases assume a static-disk–surface-density profile, it is then relatively simple to solve for the fractional mass of the contaminant m_i, in a given grid space, where i denotes the grid space, which we will normalize by m_tot, the total mass of the contaminant in the disk at t = 0, where

$$\begin{equation*} m_\mathrm{tot} = 2\pi \int _0^\infty R \Sigma C_0 \,\delta \left(R-R_0\right) dR = 2\pi R_0 \Sigma \left(R_0\right) C_0 \,. \end{equation*}$$

We write

$$\begin{equation} \left(\frac{m_i}{m_\mathrm{tot}}\right) = \left[ R_0 \Sigma \left(R_0\right) C_0 \right]^{-1} \int _{R=R_{i-1/2}}^{R_{i+1/2}} R \Sigma \left(R\right) C\left(R,t\right) dR \,, \end{equation} \tag{ A9 }$$

where R_i−1/2 and R_{i + 1/2} are the inner and outer boundaries of the grid space, respectively. Combining Equations (A7) and (A9), and using a variable substitution, s = R^1/q, we can write the general expression for the expected mass distribution as

$$\begin{eqnarray} \left(\frac{m_i}{m_\mathrm{tot}}\right) &=& \frac{q|q|R_0^{-v_0/2\zeta \nu _0}}{2\zeta \nu _0 t} \exp \left[-\frac{q^2R_0^{2-a}}{4\zeta \nu _0 t}\right]\nonumber\\ &&\times \int _{R=R_{i-1/2}}^{R_{i+1/2}} s^\gamma I_\beta \left(\frac{q^2s_0 s}{2\zeta \nu _0 t}\right) \exp \left(-\frac{q^2s^2}{2\zeta \nu _0 t}\right) ds \nonumber \\ \gamma &=& \frac{\left(2-a+v_0/\zeta \nu _0\right)}{\left(2-a\right)} \nonumber \\ s_0 &=& R_0^{1/q} \,. \end{eqnarray} \tag{ A10 }$$

We solve the integral numerically using a Simpson integrator and the Numerical Recipes functions for the modified Bessel functions (Press et al. 1992, Sections 6.6 and 6.7).

For the case of a = b = 2 given in Equation (A8) the expected mass distribution in a static disk is given by

$$\begin{equation} \left(\frac{m_i}{m_\mathrm{tot}}\right) = \left(4\pi \zeta \nu _0 t\right)^{-1/2} \int _{R=R_{i-1/2}}^{R_{i+1/2}} \exp \left[-\frac{\left(Y-v_0 t\right)^2}{4\zeta \nu _0 t} \right] dY \,, \end{equation} \tag{ A11 }$$

where Y = ln(R/R₀), and we again, we use a Simpson integrator to solve for the expected (m_i/m_tot) for these cases.

To compare our particle-transport model to the analytic test cases described in Equations (A10) and (A11), we ran simulations of 10,000 particles initiated at 40 AU; we include some simulations with zero background velocity, for which we substitute v₀ = 0 into the solution equations above (instead of the accretion case v₀ = −3ν₀/2). We vary the value of a, but all simulations use ν(1 AU) = 4.941 × 10¹⁴ cm² s⁻¹ and ζ = 1. In these test-case simulations, the values of Σ₀ and the particle size are arbitrary, since the background velocity of the particles, v_R, is prescribed.

Panel (a) of Figure A4 plots a comparison between our numerical simulations of particle transport and the analytically expected contaminant-mass distribution for the case of a = 1.5 and zero background velocity at two different times. It shows that the overall, evolving distribution of the particle ensemble is in good agreement with the analytical solution. However, the scatter produced by the discreet nature of the particle simulations makes it difficult to judge the accuracy with which the simulations reproduce this analytical test case. Therefore, panel (b) of Figure A4 plots an alternate comparison of the simulated and analytical mass distributions. Here the values plotted correspond to a summing of the fractional-mass per grid space either from R = 0 outward, or from R = R₀ = 40 AU inward (for R < R₀) and outward (for R > R₀). Again, the simulations and the analytical solutions are in good agreement.

Zoom In Zoom Out Reset image size

Simulating particle diffusion with zero background velocity in a steady disk for a range of a-values, we find, in general, that our simulations agree well with the analytical solutions. The exceptions occur for very steep (a = 3), or very shallow (a = 1) disk profiles and are caused by the finite nature of our simulation space. A viscosity (diffusivity) profile that goes as R³ leads to rapid transport of particles over very large distances at large R. Therefore, the outer edge of our simulation space acts as a sink for particles. At early times, before a significant fraction of particles have been lost, the simulations and analytical solutions are in good agreement, as shown in Figure A5; they deviate from each other at later times. In the a = 1 case, it is loss past the inner boundary of the simulation space that causes the simulations to deviate from the analytical solutions.

Zoom In Zoom Out Reset image size

Loss of simulation particles past the simulation-space boundaries is not an issue when we include the background-accretion velocity. In Figure A6, we plot the mass-distribution profiles at t = 10⁴ years for a range of a-values with the accretion velocity included. Because v_R is inward, it keeps particles away from the outer boundary and loss past the inner boundary is a part of the analytical solution as well as the simulations. All of these cases show good agreement between the simulations and the analytical solutions.

Zoom In Zoom Out Reset image size

Finally, we would like to be sure that our particle-transport model correctly simulates diffusion of the particle ensemble for the global time steps used in our static-disk models (Δt = 2000 years), and for the time steps required by the disk-evolution model (Δt ∼ 0.05 years). As was described in Section 3.3, maximum limits are placed on the particle time stepping to ensure that disk conditions local to each particle are obeyed. However, smaller time steps mean more steps to simulate a given amount of time in the disk; therefore, the probabilities of a particle stepping inward versus outward are closer and closer to equal. Cumulatively, however, all time-step sizes used must evolve the distribution of the particle ensemble equally. Therefore, we have run a few test-case simulations using time-step sizes ranging from 10⁻³ years to 100 years. For the case of a = 1 and v₀ = 0, the resultant distributions versus the analytical solutions are plotted in Figure A7. All of these simulations agree well with the analytical solution and the scatter produced for simulations with different time-step sizes are all comparable.

Zoom In Zoom Out Reset image size

While the agreement between simulations using different time-step sizes is good for almost all of the simulations we tested, one set of simulations shows a marked change in the scatter depending on the size of the time step used. Figure A8 plots simulations run for a range of time-step sizes for the case of a = b = 2 and v₀ = −3ν₀/2 compared to the analytical solutions at two times. Looking at the summed-mass distributions, it is clear that all of these simulations do distribute the particle ensemble appropriately for this diffusion test case. However, the mass-fraction-per-bin profiles show significant, resonant-like scatter in the Δt = 10 years and Δt = 1 year simulations; the scatter also appears a bit higher than usual in the Δt = 10⁻¹ year simulation. While it is not clear what causes this large scatter for these time steps, it is suggestive that the natural space for diffusion in this regime ($a=b=2 \longrightarrow \ln \left(R\right)$-space) is the same as the space wherein our disk-model grid cells are equal width. Furthermore, the resonant-like-scatter region in the Δt = 10 years simulation is restricted to outward of about 15 AU, and this region moves outward with time. Inward of 15 AU in these simulations, the Kepler time (1/Ω_K) is less than 10 years and decreases toward smaller R; the time-step sizes of individual particles at small AU are restricted to smaller than the global 10 year time step. Therefore, the resonant-like-scatter behavior is only observable when the same local time-step size is used for the entire region and when the source of all the particles is a single point in R. Nonresonant-like scatter at small AU allows the particle distribution to spread out randomly, following the local diffusivity and removing the single-source-point condition starting at small AU and moving outward.

Zoom In Zoom Out Reset image size

While this resonant-like-scatter behavior is curious and somewhat vexing, it is not a problem for the simulations presented in this paper. Aside from the fact that our model disks (Section 3.1) use ν = ν₀R (b = 1), the time-step sizes for our steady-disk and static0-disk models are large enough that the local time-step restrictions of the particle-transport model apply almost everywhere. The evolving-disk models use a time-step size that is small enough (Δt ∼ 0.05 years) that resonant-like scatter is not expected to be an issue. Note, though, that resonant-like-scatter behavior may be something to watch out for when using this method of particle-transport modeling.