THREE-DIMENSIONAL LAGRANGIAN TURBULENT DIFFUSION OF DUST GRAINS IN A PROTOPLANETARY DISK: METHOD AND FIRST APPLICATIONS

A Lagrangian code treats each particle individually, which means that it is not well suited for representing a fluid system where the mean free path is short. It is, however, well designed for tracking pointlike particles in a gas disk that have very few pairwise interactions. One very attractive advantage of this approach is that it affords the possibility of tracking individual particle trajectories, and thus, of reconstructing the thermodynamical history of particles in a protoplanetary disk environment. This should have many applications for cosmochemistry. For instance, chondrules and refractory inclusions, which are the main millimeter- to-centimeter-sized components of primitive chondritic meteorites, are thought to have undergone a complex multi-stage evolution during the first 2–3 Myr of solar system history, before being incorporated into asteroidal or cometary planetesimals at various heliocentric distances. Chondrules (Mg–Fe-rich) and refractory inclusions (also known as calcium–aluminum inclusions, CAIs hereafter) were formed in the inner solar system and experienced multiple heating and irradiation events before, or during, their transport to the region where they were finally incorporated into a meteorite. As they are the major building blocks of rocky planetesimals, tracking their thermodynamical history provides key information for deciphering the physics and chemistry of planetary formation. Similarly, the history of frozen volatile-rich dust grains from outer solar system regions is essential to understanding the solar protoplanetary disk and to unraveling the origin of planetary volatiles, such as water or organic molecules with prebiotic potential.

The disk is expected to be magnetorotational instability (MRI) turbulent, inducing an efficient mixing of the dust component and the gas component. A direct approach would be to couple a particle-based dust transport code with a 3D magnetohydrodynamic (MHD) simulation of a turbulent gas disk (as in Fromang & Nelson 2005; Johansen & Youdin 2007; Johansen et al. 2007). However, while conceptually simple, MHD simulations remain very computationally demanding, and consequently are currently limited to a few thousand orbits at most. For this reason, we turn our attention to a simplified turbulence model, the popular turbulent diffusion model, which mimics turbulent transport as a diffusive process through a Brownian motion with an efficiency parameter, the dust diffusion coefficient D_d. This diffusion coefficient is thought to be comparable in magnitude to the turbulent viscosity coefficient. More recently, Fromang & Nelson (2009) have shown that D_d may increase with the distance from the midplane. The description of turbulence as a diffusive process, though not fully accurate, is widely used and underpinned by vast literature. We have therefore built on earlier studies to develop an efficient Lagrangian code. For example, Lagrangian diffusion has been extensively used in environmental studies for the transport of air and oceanborne pollutants (see Wilson & Sawford 1996 for a review). In the planetary science literature, several models couple gas drag and turbulent diffusion within a Lagrangian framework, but this is either in a form not adapted to large particles (which decouple from the gas and undergo oscillations) or limited to some specific prescription of the gas density field. For example, in Ciesla (2010) and in Hughes & Armitage (2010), a one-dimensional (1D) stochastic diffusion model (partly similar to ours; see Section 2) is applied, but the treatment of gas density variations is different. Our method is readily applicable to any 3D gas disk and thus much more general than that used in Ciesla (2010). It is also more physically justified than in Hughes & Armitage's (2010) model, which is mainly empirical. Another major difference is our use of an implicit solver to integrate the dust motion, which allows us to accurately follow any range of particle size, from the most tightly gas-coupled particles, such as polycyclic aromatic hydrocarbons (PAHs), to those totally decoupled from the gas, such as kilometer-sized planetesimals (see the examples in Section 4.2). We will see, in particular, that proper implementation of turbulent diffusion is not straightforward, with one frequently encountered problem being that of satisfying the "good mixing condition," i.e., reaching an asymptotic state in which dust is well mixed with the gas for any gas spatial density distribution, in accordance with the Second Law of Thermodynamics. This requires a special treatment of diffusion in a Lagrangian approach and constitutes the core of this paper. The case of a non-constant diffusion is also addressed.

Another useful physical aspect of our code is that we consider dust motion in 3D. Transport of dust to altitudes high above the midplane is expected due to the efficient vertical diffusion induced by turbulence. Most studies treat radial mixing in the protoplanetary disk through a one dimensional (1D) approach (see, e.g., Gail 2001; Brauer et al. 2008; Hughes & Armitage 2010). Yet, two-dimensional (2D) models that explicitly treat the vertical motion of dust particles (see, e.g., Takeuchi & Lin 2002; Dullemond & Dominik 2004; Ciesla 2007; Tscharnuter & Gail 2007) show that the disk is stratified, which may have an impact on global transport. For example, Takeuchi & Lin (2002) show that the radial velocity of dust depends quadratically on Z inducing an outward gas drag in the disk's upper layers (for above ∼1.5 pressure scale heights, see Takeuchi & Lin 2002, Equations (11) and (17)). This depends on the gas density profile only, and is thus a robust result. We shall also see that radial diffusion is more effective at altitudes far from the midplane, due to the shallower slope of the radial density gradient for Z > 0 (see Section 3.2). The importance of including the vertical dimension is also emphasized by Bai & Stone (2010) in the context of dust transport in a dead zone. Thus, in order to ensure that the approach remains as general as possible, it is important to incorporate the vertical dimension into the system and integrate the motion of each particle with the fewest approximations possible. In the present paper, the azimuthal direction is also explicitly included. However, as there is no planet for the moment, the disk remains azimuthally symmetric, while the system is intrinsically evolved in 3D.

For practical use, this code will be referred to as LIDT3D (for Lagrangian implicit dust transport in 3D).

The paper is organized as follows: in Section 2, we describe the procedure used to introduce turbulent diffusion into a Lagrangian form. It should be noted that the gas disk considered in this paper is a simple and non-evolving parameterized gaseous disk (as in Takeuchi & Lin 2002) in order to test the dust transport algorithm. This algorithm, however, is independent of the choice of disk and can be readily extended to any disk sampled on a 3D grid. In Section 3, we present several tests aimed at reproducing known results about turbulent diffusion in a gaseous disk, and in Section 4 we discuss some applications to compare dust diffusion in 2D (at the disk midplane only) and in 3D and to show individual trajectories of different-sized dust grains extracted from the simulations.

In a Lagrangian code, the motion of an individual dust particle is described using the classical Newtonian formalism:

$$\begin{equation} \frac{{d\vec v}}{{dt}} = \frac{{\skew9\overrightarrow F _* }}{m} - \frac{{{\vec v} - {\vec v} _g }}{\tau }, \end{equation} \tag{ 1 }$$

where F_* is the gravitational force of the central star, the second term is the gas drag force, v is the particle's velocity, v_g is the gas velocity, and m is the particle's mass. The dust stopping time τ is in the Epstein regime:

$$\begin{equation} \tau = \frac{{a\rho _s }}{{\rho C_s }}, \end{equation} \tag{ 2 }$$

where a is the dust grain radius, ρ_s is the dust material density, ρ is the gas density, and C_s is the local sound velocity. Particles with sizes larger than the mean free path of gas molecules may be in the "Stokes regime," where the exact expression for τ depends on the gas Reynolds number (see Equation (10) of Birnstiel et al. 2010). However, our discussion and the method described here do not closely depend on the expression for the stopping time τ. We introduce the Stokes number St, which is the particle coupling time τ divided by the eddy turnover time τ_ed (St = τ/τ_ed). τ_ed is about 1/Ω_k (see, e.g., Fromang & Papaloizou 2006), with Ω_k representing the local Keplerian frequency (Ω_k = (GM_*/r³)^1/2), where M_*, G, and r are the star's mass, the universal gravitational constant, and the distance to the star projected along the disk midplane, respectively, such that St ∼ τΩ_k. In the laminar case, once the gas velocity field is known by applying, for example, a numerical method like that used by Tscharnuter & Gail (2007), or an analytical model like that in Takeuchi & Lin (2002; as is used here), Equation (1) is easy to solve numerically in the absence of turbulence. However, since gas is turbulent, the gas velocity field is highly complex, then it is more convenient to introduce a turbulent diffusion model into this equation to take the stochastic motion of the gas into account (see e.g., Hinze 1959; Youdin & Lithwick 2007). We first transform Equation (1) by decomposing v_g into a mean-field contribution 〈v_g〉 plus a turbulent fluctuation δv_g, such that v_g = 〈v_g 〉+ δv_g, so that Equation (1) is rewritten:

$$\begin{equation} \frac{{d\vec v}}{{dt}} = \frac{{\skew9\overrightarrow F _* }}{m} - \frac{{{\vec v} - \langle {\vec v} _g \rangle }}{\tau } + \frac{{\delta {\vec v} _g }}{\tau }, \end{equation} \tag{ 3 }$$

where δv_g/τ is the acceleration of dust induced by the turbulent gas velocity field that induces a random walk (or turbulent diffusion) of the particle, consistent with the "turbulent diffusion" model. We assume that 〈v_g〉 can be identified with the unperturbed laminar flow. This gas flow can either be numerically computed as in Ciesla (2009) or Tscharnuter & Gail (2007) or taken from analytical models as in Takeuchi & Lin (2002). Whereas the former method is more versatile, we choose to adopt the latter in the present paper for the sake of simplicity and also to enable us to test our dust transport model against analytical results. Yet, the method for including turbulent diffusion that we present below is not dependent on this choice and a gas disk sampled on a 2D or 3D grid can be easily implemented to provide the 〈v_g〉 field.

We introduce δv_T and δr_T (kicks in velocity and position) defined as δr_T = dt × δv_Tv, where dt is the time step. They are decomposed into the following elements by integrating Equation (3):

$$\begin{equation} \delta {\vec v} = \delta {\vec v} _M + \delta {\vec v} _T \end{equation} \tag{ 4 }$$

$$\begin{equation} \delta {\vec r} = \delta {\vec r} _M + \delta {\vec r} _T, \end{equation} \tag{ 5 }$$

where δv_M and δv_T are velocity increments due to the central star gravity and gas drag with the disk mean velocity field (δv_M) and gas drag with the turbulent velocity field (δv_T). As regards positions, δr_M is position increments due to initial velocity plus the star's gravity field and gas drag with the mean flow. δr_T corresponds to the position increment due to turbulent diffusion that induces a particle's random walk. The turbulent random walk of dust is entirely contained in the δr_T and δv_T terms, as discussed in Section 2.2. The term δr_M is obtained by direct numerical integration of the equation:

$$\begin{equation} \overrightarrow {\delta r_M } = \int_t^{t + dt} {v(t)\,dt} \end{equation} \tag{ 6 }$$

with

$$\begin{equation} \frac{{d\vec v}}{{dt}} = \frac{{\skew9\overrightarrow F _* }}{m} - \frac{{{\vec v} - \langle {\vec v} _g \rangle }}{\tau }. \end{equation} \tag{ 7 }$$

Equations (6) and (7) can be solved using a variety of implicit or explicit variable-order methods (see, e.g., Press et al. 1992). An implicit method is highly recommended here given that Equation (1) is a stiff equation due to the two very different timescales at play: the stopping timescale τ and the Keplerian orbital timescale T_k. The smallest particles have a very short value of τ (meaning that they are tightly coupled to the gas) that may be much smaller than T_k: for example, 0.1 μm and 1 mm particles at 1 AU in the minimum-mass solar nebula have τ/T_k about 10⁻⁷ and 10⁻³, respectively. If an explicit integrator such as the popular explicit fourth-order Runge–Kutta scheme is used, the integration time step will be limited to a fraction of the gas-coupling timescale τ. It will then be impossible to include the smallest particles (with a small Stokes number) or to integrate them over long timescales.

In the present work. we use a Bulirsch–Stoer scheme with a semi-implicit solver (Bader & Deuflhard 1983) as described in Press et al. (1992, chapter 16.6). This yields excellent results in terms of accuracy and rapidity, at least down to a Stokes number of about 10⁻⁸ when we use the Bulirsch–Stoer extrapolation method up to the eighth order in the case of a laminar flow. However, the taking into account of turbulent diffusion (see the following section) is only first-order accurate due to operator split. In this case, the Bulirsch–Stoer scheme is used with a second-order integrator. Due to the adaptive time step scheme, it yields excellent results in terms of stability whereas the intrinsic accuracy is only first order. However, extensive tests (see Section 4) have shown that our results precisely match analytical expectations even when turbulence is included.

A diffusion process is usually described by Fick's law, which relates the diffusive flux of material F_T to the density gradient:

$$\begin{equation} \skew9\overrightarrow F _T = - D \cdot \overrightarrow {{\rm grad}} (\rho), \end{equation} \tag{ 8 }$$

where D is the diffusion coefficient and ρ is the local material density. For the specific case of dust in the solar nebula, D will be written hereafter as D_d and ρ is written as ρ_d. To mimic a random walk of the particles and obtain a diffusion flux obeying Equation (8), a well-known method is to express δr_T or δv_T as Gaussian random variables with 0 mean and a variance dependant on the diffusion coefficient and the time step (see, e.g., Wilson & Sawford 1996; Hughes & Armitage 2010; and Appendix A of this paper). We now go on to consider a 1D variable only, but the procedure is readily generalized to 3D. We express the variance of δr_T as σ²_r: the variance of δv_T as σ²_v, 〈δr_T〉 and 〈δr_T〉 being their respective mean. In the "position" representation, the kick on position is a random Gaussian with

$$\begin{equation} \delta r_T = \left\{ \begin{array}{@{}l} \dsty \langle \delta r_T \rangle = 0 \\ \ms \dsty \sigma _r^2 = 2D_{} dt \end{array} \right.\!\!. \end{equation} \tag{ 9 }$$

In the "velocity" representation, the kick on velocity is a random Gaussian variable δv_T constructed so that, after time integration, it induces the same average mean dispersion on positions as δr_T (i.e., (σ²_r)^1/2 = dt × (σ²_V)^1/2):

$$\begin{equation} \delta v_T = \left\{ \begin{array}{@{}l} \dsty \langle \delta v_T \rangle = 0 \\ \ms \dsty \sigma _v^2 = \frac{{2D}}{{dt}} \end{array} \right.\!\!. \end{equation} \tag{ 10 }$$

Both methods are possible: one of the two representations must be chosen.

These kinds of procedures are known to yield good results for the random walk of a solute particle (i.e., dust) with constant D_d inside a solvent (i.e., gas) of uniform density. It has been used several times in environmental studies, notably to study the diffusion of pollutants in the atmosphere or ocean (see, e.g., Wilson & Sawford 1996). For transport of dust in the protoplanetary disk, Youdin & Lithwick (2007) and Hughes & Armitage (2010) use similar procedures in which δv_T is a discrete random variable allowed to take two values, +(2D/dt)^1/2 or −(2D/dt)^1/2.

However, problems arise when the solvent (i.e., gas) has a non-uniform density in space. In this case, the simple procedure described above is no longer valid. It can be easily verified that this method does not satisfy the "good mixing condition": at steady state, a solute (i.e., dust) has a uniform concentration throughout the solvent (i.e., gas) in order to reach maximum entropy. In other words, it must have the same spatial density as the solvent times a constant multiplicative factor. The procedure described above (Equation (9) or Equation (10)) will lead inexorably to a uniform density distribution of the solute (i.e., dust) throughout the whole space, even though the solvent (i.e., gas) has a non-uniform density. For our astrophysical case, this means that dust would diffuse everywhere in space, out of the gaseous disk itself, in the absence of central star gravity. This problem has long been identified in environmental studies and several solutions exist with different derivations (see, e.g., Sothl & Thomson 1999 or Ermak & Nasstrom 2000). Yet most of them are either expressed in terms of gas velocity fluctuations (which are not known or not "well documented" in the protoplanetary disk literature) or deal only with the self-diffusion of non-inertial material rather than being expressed in terms of diffusion coefficient, as is required here. For the case of protoplanetary disks, we wish to use only the diffusion coefficient and gas macroscopic properties in order to compare our results with previous analytical studies. To cure this problem of "good mixing," it could be tempting to introduce a spatial dependence in the diffusion coefficient D to enforce a uniform concentration at steady state. But this would run counter to the usual evaluations of the dust diffusion coefficient D_d in a turbulent disk: in an α-turbulent isothermal disk, D_d is assumed to be close to the turbulent viscosity D_d ∼ αC_sH. Since C_s and H depend only on r in an isothermal disk, D_d is a constant function of Z. Hughes & Armitage (2010) identified this problem of dust-to-gas ratio and propose an empirical method in which δv_T is a non-Gaussian discrete random variable. It is designed so that there is higher probability of a dust particle migrating toward higher density regions with weights depending on the local density profile. However, although this procedure is successful in the framework of their study, it is performed in the radial direction only and its extension to 3D systems is somewhat unclear since the number of possible directions becomes infinite. Ciesla (2010) presents a physically motivated derivation of the dust displacement algorithm, but treats the problem in the vertical direction only, for a disk density in the form ρ(z) = ρ₀exp(−z²/2H²), and within the limit of very small particles that are strongly coupled to the gas and thus follow the same path as the gas molecules in the absence of turbulence. We present below a heuristic derivation to properly account for gas density variations. While sharing some similarities with Ciesla (2010), this derivation is more general and not constrained by all the previously mentioned limitations or assumptions.

For the sake of clarity, we restrict ourselves to the case of a constant diffusion coefficient in the current section. The more general case of Lagrangian diffusion with a varying diffusion coefficient is treated in Section 2.4. We first recall some basic principles of a 1D Gaussian random walk of a solute particle in a uniform and static solvent. We call X the position of a Lagrangian particle. For a particle with Lagrangian velocity V and constant diffusion coefficient D, the spatial position increment dX due to random walk as described in Equation (9) is

$$\begin{equation} dX = Vdt + (2D\,dt)^{1/2} W, \end{equation} \tag{ 11 }$$

where W is a Gaussian random variable with mean 0 and variance σ² = 1 (such that (2D dt)^1/2 × W is a random variable with mean 0 and variance 2D dt, as in Equation (9)). A basic result of Einstein–Brownian diffusion is that the resulting motion has an average position 〈X〉 so that d〈X〉/dt = V, and that the standard deviation 〈(X−〈X〉)²〉 evolves according to d/dt(〈(X−〈X〉)²〉) = 2D. In a Lagrangian simulation involving a large number of test particles with positions that evolve according to Equation (11), the local ensemble average is a solution to the Eulerian advection diffusion equation of a solute in a solvent with uniform volume density, which reads in Cartesian coordinates:

$$\begin{equation} \frac{{\partial \rho }}{{\partial t}} + \frac{{\partial \rho \langle V\rangle }}{{\partial x}} = \frac{\partial }{{\partial x}}\left({D\frac{{\partial \rho _d }}{{\partial x}}} \right), \end{equation} \tag{ 12 }$$

also written as

$$\begin{equation} \frac{{\partial \rho }}{{\partial t}} + \frac{\partial }{{\partial x}}\left({\rho \langle V\rangle - D\frac{{\partial \rho }}{{\partial x}}} \right) = 0. \end{equation} \tag{ 13 }$$

Bearing in mind that the mass conservation equation reads ∂ρ/∂t+∂F/∂x = 0, where F is the flux, the term ρ〈V〉 in the parentheses of Equation (13) is the advective flux and −D∂ρ/∂x is the diffusive flux. Unfortunately, since the solvent (gas) density is not uniform, the transport equation of the solute (dust) does not have exactly the same form as Equation (12). Thus, it cannot be solved using the simple method of Equation (11) (equivalent to Equation (9)). In a gaseous protoplanetary disk, the transport equation of dust in the gas disk is given rather by (see, e.g., Dubrulle et al. 1995; Takeuchi & Lin 2002)

$$\begin{equation} \frac{{\partial \rho _d }}{{\partial t}} + \frac{\partial }{{\partial x}}\left({\rho _d V_d - \rho _g D_d \frac{\partial }{{\partial x}}\left({\frac{{\rho _d }}{{\rho _g }}} \right)} \right) = 0, \end{equation} \tag{ 14 }$$

where ρ_g is the density of gas, V_d, D_d, and ρ_d are the Eulerian velocity, the diffusion coefficient, and the density of dust in the disk, respectively. The difference between Equation (13) and Equation (14) is that the term D∂ρ/∂x is replaced by ρ_gD_d∂/∂x (ρ_d/ρ_g), which accounts for gas density variations. Note that when ρ_g is constant, Equation (14) reduces to Equation (13), as expected.

To build a Lagrangian approach, we wish to rewrite Equation (14) into a form functionally close to Equation (13), in which the diffusion term depends only on ∂ρ_d/∂x. This is simply done by developing the diffusion term of Equation (14):

$$\begin{equation} \frac{{\partial \rho _d }}{{\partial t}} + \frac{\partial }{{\partial x}}\left({\rho _d \left({V_d + \frac{{D_d }}{{\rho _g }}\frac{{\partial \rho _g }}{{\partial x}}} \right) - D_d \frac{{\partial \rho _d }}{{\partial x}}} \right) = 0. \end{equation} \tag{ 15 }$$

By identifying Equation (15) with Equation (13) we can recover Equation (13) simply by adding a correction term to the mean velocity of particles:

$$\begin{equation} \langle V\rangle = V_d + \frac{{D_d }}{{\rho _g }}\frac{{\partial \rho _g }}{{\partial x}}. \end{equation} \tag{ 16 }$$

This means simply that a gradient in gas density induces a net diffusive flux D/ρ_g × grad(ρ_g) directed toward higher density regions. This correction term is a systematic component of the diffusive term arising from non-homogeneous diffusion (see, e.g., Sothl & Thomson 1999; Ermak & Nasstrom 2000; Van Milligen et al. 2005). This is physically required in order to have a well-mixed final steady state.

We now come back to the mathematical implementation of dust diffusion in our Lagrangian code, including correction for the gas density gradient. In the "position" representation, the kick on position is a random Gaussian variable δr_T with mean 〈δr_T〉 and variance σ_r² now given by

$$\begin{equation} \delta r_T = \left\{ \begin{array}{@{}l} \dsty \langle \delta r_T \rangle = \frac{{D_d }}{{\rho _g }}\frac{{\partial \rho _g }}{{\partial x}}dt \\ \ms \dsty \sigma _r^2 = 2D_d\, dt \end{array} \right.\!\!. \end{equation} \tag{ 17 }$$

So the displacement during one time step is δr_T = 〈δr_T〉+ Wσ_r (W being a normal random variable). In the "velocity" representation, the kick on velocity is a random Gaussian variable δv_T with mean 〈δv_T〉 and variance σ_v² linked to δr_T through the relation δr_T = δv_T × dt, so 〈δv_T〉 = 〈δr_T〉/dt+Wσ_r/dt or in other words,

$$\begin{equation} \delta v_T = \left\{ \begin{array}{@{}l} \dsty \langle \delta v_T \rangle = \frac{{D_d }}{{\rho _g }}\frac{{\partial \rho _g }}{{\partial x}} \\ \ms \dsty \sigma _v^2 = \frac{{2D_d }}{{dt}} \end{array} \right.\!\!. \end{equation} \tag{ 18 }$$

In the velocity representation, we first integrate the motion considering only the gas drag with the mean flow and with the star's gravity. At the end of the time step, the particle's velocity is modified by adding the velocity kick as defined in Equation (18).

To illustrate the validity of these results, we show in Figure 1 simple tests of 1D diffusion of dust in a non-uniform gas medium and with periodic boundary conditions. Gas density is given a sinusoidal density profile and all test particles are released at a same starting location, X = 0.5. As there is no net transport, particles are subject to turbulent diffusion only. The diffusion coefficient is arbitrarily set to 0.001. In a first set of runs (Figure 1, left column), the density correction term is neglected such that δr_T (or δv_T) for each particle is drawn from a random distribution according to Equation (9) (kicks on positions). In Figure 1, we see that after 1000 time steps a steady state is reached in which the absolute spatial density of particles is uniform (Figure 1(a)), whereas the gas is not uniform. As a result, the final dust concentration profile is not uniform (Figure 1(c)), which is not physical.

Zoom In Zoom Out Reset image size

In a second run, we include the correction term, thus following Equation (17). The time evolution of the system is shown in the right column of Figure 1. We see that the system tends toward a state where dust concentration in the gas is asymptotically uniform (Figure 1(b)) and takes on a volume density profile similar to the profile of the gas, up to a constant multiplicative factor (Figure 1(d)), which is physical. These results illustrate that our numerical procedure naturally allows diffusion to tend toward a state of uniform concentration, in agreement with the "good mixing condition." This will be tested further in "real conditions" in Section 3.

We now consider a coefficient of diffusion that varies in space, so D is now written as D(x). Such a behavior may be encountered in many physical situations: for example, numerical simulations (see, e.g. Fromang & Nelson 2009; Turner et al. 2010) show that MRI turbulence is not uniform vertically in the disk and that the velocity fluctuations δV_z tend to increase with Z. In consequence, the effective gas diffusion coefficient increases according to D_g ∼ δV_z²/Ω. This may be especially important in the case a "dead zone" (a laminar region) is present in the disk's midplane under an active upper layer. In this case, D_g may vary by several orders of magnitudes on distances of a few scale heights only.

It can be shown that when D(x) has a non-zero gradient, the classic procedure of a random walk whose step is a random variable with 0 mean and 2D dt standard deviation (like in Equation (11)) is not anymore an accurate solution to the equation ∂X/∂t = ∂²Dx/∂t². The right random variable to consider in this case has been studied by our colleagues in environmental studies (see, e.g., Ermak & Nasstrom 2000), and is given by (see Ermak & Nasstrom 2000 and Appendix A) in the position representation:

$$\begin{equation} \delta r_T = \left\{ \begin{array}{@{}l} \dsty \langle \delta r_T \rangle = \frac{{\partial D(x)}}{{\partial x}}\,dt \\ \ms \dsty \sigma _r^2 = 2D(x)\,dt + \left({\frac{{\partial D(x)}}{{\partial x}}dt} \right)^2 \end{array} \right. \end{equation} \tag{ 19 }$$

and in the velocity representation

$$\begin{equation} \delta v_T = \left\{ \begin{array}{@{}l} \dsty \langle \delta v_T \rangle = \frac{{\partial D(x)}}{{\partial x}} \\ \ms \dsty \sigma _v^2 = \frac{{2D(x)}}{{dt}} + \left({\frac{{\partial D(x)}}{{\partial x}}} \right)^2 \end{array} \right.\!\!. \end{equation} \tag{ 20 }$$

Ermak & Nasstrom (2000) suggest considering a distribution with a non-zero skewness (the skewness is the third moment of a distribution, and is 0 for a symmetric distribution) to make the result better. In the current paper, we found excellent result considering only symmetric distributions (i.e., a Gaussian) as described by Equation (19) or Equation (20). We now need to put all things together to treat the most general case.

We now treat the most general case of dust in the protoplanetary disk, with a varying dust diffusion coefficient D_d and also with a varying gas density (the solvent). The good random variable for the dust random walk is simply obtained by considering that in the absence of a gas-density gradient the random walk is described by Equation (19) and that when a gas-density gradient is present, an additional term appears on 〈δr_t〉 only given by Equation (17). So in the position representation, δr_t is now

$$\begin{equation} \delta r_T = \left\{ \begin{array}{@{}l} \dsty \langle \delta r_T \rangle = \frac{{D_d }}{{\rho _g }}\frac{{\partial \rho _g }}{{\partial x}}dt + \frac{{\partial D_d (x)}}{{\partial x}}\,dt \\ \ms \dsty \sigma _r^2 = 2D_d (x)\,dt + \left({\frac{{\partial D_d (x)}}{{\partial x}}\,dt} \right)^2 \end{array} \right.\!, \end{equation} \tag{ 21 }$$

and in the velocity representation

$$\begin{equation} \delta v_T = \left\{ \begin{array}{@{}l} \dsty \langle \delta v_T \rangle = \frac{{D_d }}{{\rho _g }}\frac{{\partial \rho _g }}{{\partial x}} + \frac{{\partial D_d (x)}}{{\partial x}} \\ \ms \dsty \sigma _v^2 = \frac{{2D_d (x)}}{{dt}} + \left({\frac{{\partial D_d (x)}}{{\partial x}}} \right)^2 \end{array} \right.\!\!. \end{equation} \tag{ 22 }$$

Equations (21) and (22) are the core result of the present paper.

We have seen above that two methods are possible: either a kick on positions or a kick on velocities. For example, Ciesla (2010) uses a kick on position, whereas Hughes & Armitage (2010) and Youdin & Lithwick (2007) use a kick on velocities. Which is the better choice? Both methods have their own caveats and suffer from the fact that Brownian motion is still a crude physical model that hides unphysical infinite velocities and/or accelerations.

• 1.  
  A kick on position is an instantaneous transport that depends on the time step (Equation (21)). It induces a discrepancy between the velocity field and the position field because of the sheared nature of a Keplerian disk. This may, however, be partially cured by applying the kick at the intermediate time step and by self-consistently correcting velocity at the end of the time step.
• 2.  
  A kick on velocity depends on the time step (Equation (22)), inducing again velocity dispersion among particles that are time step dependant.

After testing both procedures, it turns out that a kick on positions seems to yield somewhat better results for two reasons:

• 1.  
  In the limit of dt → 0, which is often desirable for numerical accuracy and stability, the magnitude of kicks on positions converges toward 0 (Equation (21)), which is numerically stable, whereas kicks on velocity diverge (Equation (22)), which is numerically unstable.
• 2.  
  Since the velocity of small particles is efficiently damped by gas drag, velocity kicks may induce a strong gas drag that brakes accelerations a_b, the magnitude of which is ∼δv_T/τ ∼ (2D_d/τ²dt)^1/2. Generally, the time step, dt, is a small fraction of the orbital period, i.e., dt = ε/Ω_k with ε ∼ 0.1%–10%. Since D_d is of the order of the turbulent viscosity D_d ∼ αH²Ω_k, a_b is proportional to Ω_k/(τε^1/2). Small particles with a very small Stokes number will thus be subject to very high accelerations requiring very small integration time steps to ensure accuracy and/or stability. Reducing the time-step will induce a decrease of ε and will not solve the problem as it increases the velocity kick. It may then be impossible to integrate very small particles in the limit St → 0.

For these reasons, integration with velocity kicks appeared to be about five times slower than with position kicks (for particles with St > 10⁻⁵) using an implicit Bulirsch–Stoer solver with adaptive time step control (Press et al. 1992). We thus chose to integrate particle motion using kicks on positions. This yields good results for all particle sizes, as shown in Section 3.

Various expressions exist in the literature for D_d. In a turbulent protoplanetary disk it is generally thought that turbulence is the main process driving the outward transport of angular momentum (and thus mass accretion) and diffusion of dust. Numerical simulations (see, e.g., Fromang & Nelson 2009) have shown that the effective gas diffusion coefficient, D_g, is close to the turbulent viscosity that describes the transport of angular momentum. The diffusion coefficient D_d is linked to the turbulent viscosity through the Schmidt number Sc, defined as the ratio of the gas diffusion coefficient (D_g) to the dust diffusion coefficient (D_d). In other words, Sc = D_g/D_d. D_g is approximately αC_sH as measured in numerical simulations of MRI turbulence (where α is the turbulent viscosity parameter and H is the local pressure scale height of the gas disk; see, e.g., Fromang & Papaloizou 2006) this gives

$$\begin{equation} D_d \approx \frac{{\alpha C_s H}}{{{\rm Sc}}}. \end{equation} \tag{ 23 }$$

Sc depends on the particle size. It tends to 1 in the limit of very small particles and increases to infinity for very big particles. This means that very big particles no longer "feel" the turbulence because of their inertia. The expression of Sc is a matter of debate and, in the literature, depends on the value given to St. Following the work of Youdin & Lithwick (2007) we assume the correlation time to be the eddy time (τ_{c =} τ_edd). The eddy turnover time is τ_edd ∼ 1/(α^2γ-1Ω_k). The value of γ depends on the structure of the turbulence: if turbulent diffusion is driven by large slow-moving turbulent eddies then γ → 1, but if driven by small fast-moving turbulent eddies γ→0 (Brauer et al. 2008). However, several studies (Cuzzi et al. 2001; Schräpler & Henning 2004) report γ = 1/2, such that τ_edd ∼ 1/Ω_k (Ω_k representing the local Keplerian frequency) in agreement with Fromang & Papaloizou (2006), who measure a τ_c ∼ 0.15 orbital period. In LIDT3D, the time correlation in the kicks was introduced stochastically, similarly to Youdin & Lithwick (2007): at each time step and for every particle a random number is generated, and the probability of applying a kick is the ratio dt/τ_c. Thus, particles close to the central star will be subject to kicks much more frequently than those at greater distances. This is illustrated in Section 4.2.

The expression of the Schmidt number is a matter of some debate. Dullemond & Dominik (2004), for example, use Sc = 1+St. Recently, however, in their paper on diffusion driven by anisotropic turbulence, Youdin & Lithwick (2007) suggest that Sc = (1+Ω_k²τ²)²/(1+4 Ω_kτ), for radial diffusion only, on the basis of an analytical study coupling sedimentation and orbital oscillation. This is in agreement with the numerical simulations of Carballido et al. (2006). This expression will be used in our code unless otherwise specified. In addition to these analytical expressions, some recent papers (Fromang & Papaloizou 2006; Fromang & Nelson 2009) have attempted to characterize the diffusion of small dust particles using direct 3D MHD simulations. They show that the vertical distribution of dust does not fit a Gaussian distribution (as often assumed in analytical models), with major discrepancies for the upper layers (Z/H > 1). As explained by the authors, this may be due to larger velocity fluctuation (ΔV_z) of the gas for increasing values of Z/H. As a result, D increases quadratically with Z and has a value of ∼0.002C_sH at the midplane. They provide a basic fit to their results in Equation (26) of Fromang & Nelson (2009). Non-uniform dust diffusion will be presented in Section 4.1.

In the rest of this paper, a 3D analytical model of a locally isothermal gas disk is used for the sake of simplicity. It provides gas density, temperature, and velocity field. However, the method presented in Section 2 is easily applied to any disk sampled on a 3D grid provided that the density gradient can be computed. Here we use the disk structure described in Takeuchi & Lin (2002). The 3D gas density as a function of R and Z is

$$\begin{equation} \rho _g (r,z) = \rho {}_0r^p e^{\frac{{ - Z^2 }}{{2H(r)^2 }}}, \end{equation} \tag{ 24 }$$

where H²(r) = H₀²r^q+3 and ρ₀ is the gas density at unit radius, and p and q are two constants. With these notations, the exponent of the surface density profiles is s = p+(q+3)/2 (Takeuchi & Lin 2002). We use the fiducial disk model of Brauer et al. (2008) with indices p = −2.25 and q = −0.5, corresponding to a surface density with exponent s = −1, a gas surface density of 200 kg m⁻² and a temperature of 204 K at 1 AU with a central star mass of 0.5 M_☉. Unless otherwise specified, the turbulent parameter is α = 0.01 and uniform throughout the disk. In order to isolate transport effects due to diffusion and sedimentation from those caused by advection in the gas flow, we fix the radial and vertical gas velocity at zero. Consequently, the nebula considered in the following tests is simply a modified minimum-mass nebula. In a forthcoming paper, this simple (but convenient) model will be replaced by a 2D time-evolving disk, as in the work of Ciesla (2009).

Turbulent diffusion opposes dust settling in the vertical direction. The dust thus attains a vertical steady-state distribution, which in general is not a Gaussian (strictly speaking). For small dust particles that reach terminal velocity rapidly (in less than one orbital period), the equation governing the evolution of equilibrium dust density along the Z direction is (Dullemond & Dominik 2004; Dubrulle et al. 1995)

$$\begin{equation} \frac{{\partial \rho _d }}{{\partial t}} = \frac{\partial }{{\partial z}}\left({\rho _d \tau \Omega _k^2 z} \right) + \frac{\partial }{{\partial z}}\left({D_d \rho _g \frac{\partial }{{\partial z}}\left({\frac{{\rho _d }}{{\rho _g }}} \right)} \right). \end{equation} \tag{ 25 }$$

The term ρ_dτΩ_k² is the particle's terminal velocity in the gas, which is assumed to be reached in a short time (in other words, the stopping timescale is much shorter than the diffusion and orbital timescales, i.e., St ≪1), which means that Equation (25) applies to small particles only. When D_d is constant, this equation can be analytically solved to determine the steady-state vertical distribution of tightly coupled particles (Dubrulle et al. 1995; Fromang & Nelson 2009):

$$\begin{equation} \rho _d = \rho _{d,{\rm mid}} \exp \left({\frac{{(\Omega _k \tau)_{{\rm mid}} C_s H}}{{D_d }}\big({e^{\frac{{z^2 }}{{2H^2 }}} - 1} \big) - \frac{{z^2 }}{{2H^2 }}} \right). \end{equation} \tag{ 26 }$$

For particles with stopping times much larger than their orbital period and which are loosely coupled to the gas (⇔ Ω_kτ ⩾ 1 ⇔ St ⩾ 1 assuming that the eddy turnover time is ∼1/Ω_k), vertical motion is mainly driven by their orbital oscillation. The scale height H_p of these large particles can be estimated using the argument of Youdin & Lithwick (2007). The time they need to return to the midplane after a vertical perturbation is given by their stopping time τ. Thus, H_p is easily found by equating the diffusion time over H_p (∼H_p²/D_d) to the stopping time (τ) implying H_p = (D_dτ)^1/2 (equivalent to H_p = (D_g/τΩ_k²)^1/2 reported in Equation (6) of Youdin & Lithwick 2007 since D_d ∼ D_g/(τΩ_k)² according to their Equation (4)). This is more naturally expressed in terms of the Stokes number (τ ∼ St/Ω_k), which gives H_p = (D_d St/Ω_k)^1/2. Note that D_d depends on the particle's size through the Schmidt number. H_p is equivalent to Equation (6) of Youdin & Lithwick (2007) but, in our expression, D_d is not explicitly given in terms of the particle stopping time. Assuming that particles are vertically distributed along a Gaussian with scale height H_p, it is possible to derive a simple expression for the vertical distribution of particles in the regime St ⩾ 1 (loose coupling):

$$\begin{equation} \rho _d = \rho _{d,{\rm mid}} e^{ - \frac{{ - z^2 \Omega _k }}{{2(D_d \cdot {\rm St})}}}. \end{equation} \tag{ 27 }$$

We now verify below whether the analytical distributions of dust in both the tightly coupled regime (Equation (26)) and the loosely coupled regime (Equation (27)) are reproduced with the method described in Section 2.5, which is the core of the LIDT3D code.

We first run a simulation with particles orbiting in an annulus extending from r_min = 1 AU to r_max = 1.01 AU in a gaseous disk as described in Section 3.1. We use periodic radial boundary conditions (i.e., when a particle is found at position r_max+ dr, dr > 0, it is shifted to position r_min+dr modulo r_max−r_min). The particle's velocity is also corrected to account for the Keplerian shear. Initially 70,000 different-sized particles (10,000 particles per size bin: 0.1, 1, 10, ..., 10⁵ μm) are uniformly distributed vertically from Z = −2H to Z = +2H. Particles evolve under the action of gas drag and turbulent diffusion, using the algorithm described in Section 2.5 (note here that D_d is constant with height so that all terms in ∂D_d/∂Z vanish). The Schmidt number is assumed to be a constant, Sc = 1.5. In Figure 2, particle location is plotted after 20,000 orbits and clearly shows that particles are vertically sorted according to their sizes, with the smallest dust particles floating in the upper layers and the largest orbiting close to the midplane.

Zoom In Zoom Out Reset image size

Vertical distributions of same-sized particles are plotted in Figure 3 (solid line) and compared to Equation (26) (dashed line) and Equation (27) (dotted line in Figure 3). The vertical distribution obtained with LIDT3D closely matches the analytical estimates as shown hereafter. We see in Figure 3 that, for Ω_kτ ranging from 10⁻⁶ to 10⁻², LIDT3D produces an excellent solution, even though the integration time step is about 10⁴ times the particle stopping time owing to the efficiency of the implicit solver. For Ω_kτ > 0.1, the assumption of tight gas–particle coupling breaks down. This explains the discrepancy between the dashed and solid lines in Figure 3(a). To explore the loose-coupling regime, an additional simulation is run for particles with Ω_kτ = 2.5 and 25, with the assumption that the Schmidt number increases with the particle stopping time Sc ∼ 1 + (Ω_kτ)² (valid for radial diffusion, as stated in Youdin & Lithwick 2007, whereas it is used here as an isotropic diffusion coefficient, as in Birnstiel et al. 2010). The result is plotted in Figure 4 against the analytical models for tight coupling (Equation (26)) and loose coupling (Equation (27)). The distributions produced by LIDT3D are closely matched by the analytical distribution for loose coupling (dotted line).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

These tests show that for both small and large particles, LIDT3D produces a good physical solution. An additional test is also presented in Appendix B for the case of a vertically varying diffusion coefficient. Any particle size can thus be treated self-consistently with a single tool, from submicron dust to macroscopic boulders. However, note that for dust particles at St ≫1, large-scale correlated motion may also arise in MRI simulations (i.e., velocity fluctuations with long correlation distances; S. Fromang 2010, private communication). We have not yet taken these into consideration. However, as this is not an intrinsic limitation of the method, we plan to introduce such spatial correlations into the system in the near future.

To test radial diffusion against simple models, we run a purely 2D simulation of dust transport in the midplane (Z = 0). We consider only very small particles (0.1 μm), with Stokes numbers in the range 10⁻⁷–10⁻⁸ between 1 and 10 AU. Since the particles are tightly coupled to the gas, which is not moving radially (V_r = 0), their radial motion results uniquely from turbulent diffusion (with α = 0.01). Ten thousand particles were released at 10 AU and evolved using the procedure described in Section 2.5. The resulting particle volume density as a function of distance is reported in Figure 5 (colored lines). As expected, the density distribution widens with time, adopting a somewhat asymmetric shape owing to the tendency for particles to mix well with the gas. We compare the distributions obtained with our stochastic algorithm with the numerical integration of the advection–diffusion equation in the radial direction given by

$$\begin{equation} \frac{{dC}}{{dt}} = \frac{1}{{\rho _g r}}\frac{\partial }{{\partial r}}\left({ - rV_r \rho _g C + r\rho _g D_d \frac{{\partial C}}{{\partial r}}} \right), \end{equation} \tag{ 28 }$$

where C is the concentration of dust in the gas, D_d is the diffusion coefficient, and V_r is the radial velocity of the gas (= 0 here). Equation (28) is numerically integrated using a second-order derivative operator (centered difference) and a first-order time scheme with a time step of about 0.1 yr. The result from this Eulerian model is shown by the black solid line in Figure 5. The agreement with LIDT3D (colored lines) is excellent, apart from the regions below 2 AU at T > 10,000 yr due to edge effects (all particles passing below 1 AU are eliminated from our particle simulation in this example).

Zoom In Zoom Out Reset image size

In the examples considered in the present paper (Sections 3 and 4), numerous tests showed that ∂D_d/∂x appeared to play solely a minor role as far as D_d is defined as a constant fraction of C_sH (like in the standard formulation D_d ∼ αC_sH with α = constant) which is the standard case of the current paper. However, for disks containing a central "dead zone" in their midplane, D_d may be strongly varying with Z. Thus the gradient of D_d may even dominate in some cases. These results will be presented in a forthcoming paper concentrating on the subject of diffusion in dead zones.

To conclude this section, we have shown that our method for calculating Lagrangian diffusion fulfills the "good mixing" constraint for both vertical and radial diffusion; it is in close agreement with analytical models for vertical mixing (in regimes of both tight and loose coupling to the gas), and in excellent agreement with the numerical integration of the diffusion–advection equation for radial diffusion.

In many studies, both analytical and numerical, dust diffusion is studied along the radial direction only, so as to make computation tractable. In order to do this, the authors solve the radial transport equation (Equation (28)) by applying one of the following three assumptions: either they average dust density in the vertical direction and assume perfect mixing (Gail 2001), or consider diffusion only at the disk's midplane (Hughes & Armitage 2010), or vertically integrate the diffusion equation assuming a Gaussian vertical distribution of dust (Brauer et al. 2008; Birnstiel et al. 2010). Moreover, in most of the aforementioned studies, particle dynamics is cast in the simplified form of asymptotic dynamics (at times much larger than the stopping time), and is most often computed at the midplane only. Whereas these approaches are very useful for understanding the large-scale properties of dust diffusion, the fact that they do not explicitly resolve the Z direction may lead to uncertainties or biases. None of the aforementioned approximations are really physically justified as, indeed, dust particles never reach perfect vertical mixing with the gas (see Section 3.2) or, due to turbulence affects, remain confined to the midplane. They are always in a somewhat intermediate state. At this point, we would emphasize the need to track dust motion in the R and Z directions simultaneously and provide some comparisons of 3D simulations against purely midplane cases. It is to be expected that dust diffusion is more effective in the disk's upper layers for the following reasons: For a gas disk structure such as in Equation (24), the magnitude of the gas density gradient is smaller in the disk's upper layers than at the midplane (see the top panel of Figure 6). This has two important consequences:

• 1.  
  As described in Takeuchi & Lin (2002), for Z > 1.5H the gas rotation profiles become super-Keplerian, inducing an outward gas drag force (see Section 3.1 of Takeuchi & Lin 2002), which favors an outward migration of dust particles. This effect cannot be captured unless the Z direction is explicitly taken into account, as is the case here (see below).
• 2.  
  As dust particles need to reach a perfect mixing with gas in a Lagrangian description, they are subject to a diffusion velocity that is proportional to D_d/ρ_g × grad(ρ_g), which is strongly directed inward for Z = 0. However, as Z increases and the density gradient decreases, this inward diffusion velocity becomes increasingly weaker (in magnitude) and can even become positive for some values of Z/R (see the example in the bottom panel of Figure 6). The net effect is again an increased efficiency of diffusion in the disk's upper layers.

Zoom In Zoom Out Reset image size

All these elements favor a non-homogeneous diffusion in the vertical direction, in which diffusion efficiency is an increasing function of altitude. We present below some illustrations of this, comparing them to the case of dust dynamics in the midplane only. It is far beyond the scope of the current paper to quantify these effects precisely, but we wish to provide examples in a disk with a surface density decreasing as r⁻¹ and with α = 0.01.

As a qualitative example, we first run a simulation using 5000 particles, with sizes ranging from 0.1 μm to 1 cm, starting at 10 AU. We prevent radial migration for 500 yr so as to reach a vertical steady state and then allow the system to evolve radially (Figure 7). The Schmidt number here is Sc = (1+Ω_k²τ²)²/(1+4Ω_kτ) (valid for diffusion in the radial direction; Youdin & Lithwick 2007). As expected there is a clear gradient of particles with Z. The particle cloud does not spread symmetrically around the release location. In fact, particles with Z > H spread radially more rapidly (Figure 7 at 3000 and 5500 yr) than those with Z < H.

Zoom In Zoom Out Reset image size

To observe the spatial and time evolution of dust distribution and to emphasize the differences with the case of midplane diffusion only, we run two simulations with 0.1 μm particles: one fully three dimensional (3D cases) and the other two dimensional confined to the disk's midplane (2D case). In both cases, 50,000 particles are initially released at 5 AU from the central star. The resulting surface density of dust as a function of time is plotted in Figure 8. After a time, particles spread in the disk and try to reach a "good mix" with the gas disk, and thus rapidly adopt a power-law surface density profile decreasing with R. The difference between the 2D and 3D cases is striking: in the 3D case, particles spread outward much more rapidly than in the 2D case, In addition, the particle surface density is systematically higher in the 3D case than in the 2D case for regions further than the initial release point. This result is easily explained: in the 2D case, due to the strong negative density gradient of gas in the midplane, the diffusive flux is strongly directed inward due to the diffusive flux in D/ρ_g × grad(ρ_g), in such a way that outward diffusion is to a certain extent prevented. Conversely, in the 3D case, particles can explore the upper layers of the disk where the gas density gradient is much shallower, or even positive, which thus favors outward diffusion. Note that a zero radial velocity of the gas is used in order to isolate the diffusion effect.

Zoom In Zoom Out Reset image size

The present example illustrates simply that dust transport may be very significantly different when we explicitly consider the vertical direction of the disk, and that simulation of dust diffusion in the midplane only, or when perfect vertical mixing is assumed, may lead to large uncertainties, or even errors.

Our second example involves the direct tracking of dust motion in protoplanetary disks. One of the strengths of a Lagrangian description is that individual particle paths in the disk can be extracted from simulations. This can be useful for deciphering the thermochemical histories of dust and grains that will be incorporated into meteorites or planets. In Figures 9–11, we present different examples of the trajectories of particles with increasing sizes. Initially, all the particles were released in the disk's midplane at 5 AU from the central star. Different behaviors can be readily identified. Small micrometer-sized particles have very stochastic paths in the R/Z plane, with their (R, Z) location increasing or decreasing stochastically (see, e.g., Figures 9 and 10). At larger distances, we also note that the frequency of turbulent kicks decreases while their magnitude increases: this is due to the eddy correlation time (see Section 3.3), which increases with the distance to the star. As a result, turbulent kicks are less frequent for increasing values of R, while the magnitude of the kick (the diffusion coefficient), such as the turbulent viscosity, increases with R (∝ R^q+3/2; q = −0.5 represents the radial exponent index of the gas disk temperature). Our last example involves the dynamics of a 10 cm radius particle (Figure 11), which shows significant differences from the previous cases, as it drifts rapidly toward the central star due to strong gas drag. As is clearly visible in Figure 11, it is also subject to vertical oscillations around the disk's midplane, the amplitude of which diminishes as the particle gets closer to the star.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

During their journey in the disk, the dust particles go through different episodes of heating or cooling (see the bottom panels of Figures 9–11) which may affect their chemical and isotopical composition due to multiple adsorptions and desorptions of the surrounding gas. Similarly, refractory inclusions in meteorites and comets have undergone multiple heating events in gaseous environments with drastically different oxygen isotopic compositions (¹⁶O-enriched solar composition and ¹⁶O-depleted planetary composition) and oxygen fugacity (reducing and oxidizing). Thus tracking their trajectories together with their thermodynamical histories may reveal critical clues on how oxygen reservoirs were distributed in the solar protoplanetary disk. Because oxygen is the major rock-forming element, understanding how its isotopic composition evolved from solar to planetary is key to understand the physics of planetary formation.

These trajectories will be used in a forthcoming paper to study the chemical evolution of dust during its transport in the protoplanetary disk.