# 3D Lagrangian turbulent diffusion of dust grains in a protoplanetary disk: method and first applications

## Summary (5 min read)

### 1. INTRODUCTION

- The transport of solids in turbulent, gaseous protoplanetary disks involves highly diverse physical processes such as gas drag, turbulence, photophoretic gas pressure and radiation pressure.
- Here the authors focus on the motion of dust particles under gas drag and turbulent diffusion.

### 1.1. A Lagrangian Approach

- 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.
- This should have many applications for cosmochemistry.
- 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.
- 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.

### 1.2. Turbulent Diffusion

- The disk is expected to be magnetorotational instability (MRI) turbulent, inducing an efficient mixing of the dust component and the gas component.
- Lagrangian diffusion has been extensively used in environmental studies for the transport of air and oceanborne pollutants (see Wilson & Sawford 1996 for a review).
- It is also more physically justified than in Hughes & Armitage’s (2010) model, which is mainly empirical.
- Another major difference is their 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).
- The authors 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.

### 1.3. Three-dimensional System

- Another useful physical aspect of their code is that the authors consider dust motion in 3D.
- This depends on the gas density profile only, and is thus a robust result.
- 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 Section 3, the authors present several tests aimed at reproducing known results about turbulent diffusion in a gaseous disk, and in Section 4 they 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.

### 2.1. Basic Concepts and Definitions

- 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).
- The authors introduce the Stokes number St, which is the particle coupling time τ divided by the eddy turnover time τ ed (St = τ/τ ed).
- As regards positions, δrM is position increments due to initial velocity plus the star’s gravity field and gas drag with the mean flow.
- The turbulent random walk of dust is entirely contained in the δrT and δvT terms, as discussed in Section 2.2. (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).
- The Bulirsch–Stoer scheme is used with a second-order integrator.

### 2.2. The “Good Mixing” Problem

- To mimic a random walk of the particles and obtain a diffusion flux obeying Equation (8), a well-known method is to express δrT or δvT as Gaussian random variables with 0 mean and a variance dependant on the diffusion coefficient and the time step .
- (10) Both methods are possible: one of the two representations must be chosen.
- 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.
- 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 selfdiffusion of non-inertial material rather than being expressed in terms of diffusion coefficient, as is required here.

### 2.3. Good Diffusion with a Constant Diffusion Coefficient

- The more general case of Lagrangian diffusion with a varying diffusion coefficient is treated in Section 2.4.
- Note that when ρg is constant, Equation (14) reduces to Equation (13), as expected.
- The diffusion coefficient is arbitrarily set to 0.001.
- In a first set of runs , the density correction term is neglected such that δrT (or δvT ) for each particle is drawn from a random distribution according to Equation (9) (kicks on positions).
- The authors see that the system tends toward a state where dust concentration in the gas is asymptotically uniform ) and takes on a volume density profile similar to the profile of the gas, up to a constant multiplicative factor ), which is physical.

### 2.4. Good Diffusion with a Varying Diffusion Coefficient

- 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 δVz tend to increase with 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 the current paper, the authors found excellent result considering only symmetric distributions (i.e., a Gaussian) as described by Equation (19) or Equation (20).
- The authors now need to put all things together to treat the most general case.

### 2.5. Good Diffusion: Putting All Things Together

- The authors now treat the most general case of dust in the protoplanetary disk, with a varying dust diffusion coefficient Dd 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 〈δrt〉 only given by Equation (17).

### 2.6. Numerical Considerations: Kicks on Position or on Velocity?

- 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.
- It induces a discrepancy between the velocity field and the position field because of the sheared nature of a Keplerian disk.
- Reducing the time-step will induce a decrease of ε and will not solve the problem as it increases the velocity kick.
- For these reasons, integration with velocity kicks appeared to be about five times slower than with position kicks (for particles with St > 10−5) using an implicit Bulirsch–Stoer solver with adaptive time step control (Press et al. 1992).

### 2.7. Expressions of the Diffusion Coefficient

- Various expressions exist in the literature for Dd. Numerical simulations (see, e.g., Fromang & Nelson 2009) have shown that the effective gas diffusion coefficient, Dg, is close to the turbulent viscosity that describes the transport of angular momentum.
- 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.
- This expression will be used in their code unless otherwise specified.
- As explained by the authors, this may be due to larger velocity fluctuation (ΔVz) of the gas for increasing values of Z/H. Non-uniform dust diffusion will be presented in Section 4.1.

### 3.1. Gas Disk Model

- It provides gas density, temperature, and velocity field.
- Here the authors use the disk structure described in Takeuchi & Lin (2002).
- 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, the authors fix the radial and vertical gas velocity at zero.
- The nebula considered in the following tests is simply a modified minimum-mass nebula.

### 3.2. Vertical Diffusion of Particles: Testing Small and Large Particles

- Turbulent diffusion opposes dust settling in the vertical direction.
- Vertical distributions of same-sized particles are plotted in Figure 3 (solid line) and compared to Equation (26) (dashed line) and Equation (27) .
- 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 +.
- The authors have not yet taken these into consideration.

### 3.3. Radial Diffusion in the Midplane

- Since the particles are tightly coupled to the gas, which is not moving radially (Vr = 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.
- 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.
- For disks containing a central “dead zone” in their midplane, Dd may be strongly varying with Z.

### 4.1. Comparing Midplane and Out-of-plane Diffusion

- In many studies, both analytical and numerical, dust diffusion is studied along the radial direction only, so as to make computation tractable.
- This effect cannot be captured unless the Z direction is explicitly taken into account, as is the case here (see below).
- The particle cloud does not spread symmetrically around the release location.
- 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.
- The present example illustrates simply that dust transport may be very significantly different when the authors 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.

### 4.2. Tracking Particle Paths in Protoplanetary Disks: A Bridge Toward Cosmochemistry

- The authors second example involves the direct tracking of dust motion in protoplanetary disks.
- In Figures 9–11, the authors present different examples of the trajectories of particles with increasing sizes.
- At larger distances, the authors 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.
- These trajectories will be used in a forthcoming paper to study the chemical evolution of dust during its transport in the protoplanetary disk.

### 5. CONCLUSION

- The gas disk is treated separately from the dust and no retroaction of dust on gas is considered for the moment.
- In SPH simulations of gas and dust (as in Barrière-Fouchet et al. 2005 or in Fouchet et al. 2010 which includes a planet), the method presented in Section 2 provides a direct way to introduce a perturbation into dust–particle motion that mimics the effect of turbulence without the need for an explicit MHD code algorithm.
- The authors are indebted to an anonymous referee whose careful review, patience, and tolerance resulted in a much improved paper.

Did you find this useful? Give us your feedback

##### Citations

^{1}, Harvard University

^{2}, Pontifical Catholic University of Chile

^{3}, Ludwig Maximilian University of Munich

^{4}, Heidelberg University

^{5}, Rice University

^{6}, University of Chile

^{7}, University of Grenoble

^{8}, Centre national de la recherche scientifique

^{9}, Tsinghua University

^{10}, California State University, Northridge

^{11}

337 citations

196 citations

157 citations

157 citations

105 citations

##### Related Papers (5)

##### Frequently Asked Questions (2)

###### Q2. What have the authors stated for future works in "Three-dimensional lagrangian turbulent diffusion of dust grains in a protoplanetary disk: method and first applications" ?

In this paper, the authors have presented a numerical tool ( LIDT3D ) for simulating 3D transport of dust in a turbulent gaseous protoplanetary disk using an implicit Lagrangian approach, in order to allow tracking of individual particles in gas. In the first example, the authors show that diffusive transport in the upper layers of the disk is more efficient than in the disk ’ s midplane and may significantly increase the radial transport of dust. In particular, the authors show that the density of particles diffusing in the disk ’ s midplane only is about 5–10 times less than for the 3D case, in which particles are allowed to move vertically in the disk ( Figure 8 ). The authors then plan to use this to compute the equilibrium distribution of dust in observed protoplanetary disks and couple the results with a radiative transfer code to compare them to observations.