scispace - formally typeset
SciSpace - Your AI assistant to discover and understand research papers | Product Hunt

Journal ArticleDOI

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

AbstractIn order to understand how the chemical and isotopic compositions of dust grains in a gaseous turbulent protoplanetary disk are altered during their journey in the disk, it is important to determine their individual trajectories. We study here the dust-diffusive transport using lagrangian numerical simulations using the the popular "turbulent diffusion" formalism. However it is naturally expressed in an Eulerian form, which does not allow the trajectories of individual particles to be studied. We present a simple stochastic and physically justified procedure for modeling turbulent diffusion in a Lagrangian form that overcomes these difficulties. We show that a net diffusive flux F of the dust appears and that it is proportional to the gas density ({\rho}) gradient and the dust diffusion coefficient Dd: (F=Dd/{\rho}\timesgrad({\rho})). It induces an inward transport of dust in the disk's midplane, while favoring outward transport in the disk's upper layers. We present tests and applications comparing dust diffusion in the midplane and upper layers as well as sample trajectories of particles with different sizes. We also discuss potential applications for cosmochemistry and SPH codes.

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

...read more

Content maybe subject to copyright    Report

HAL Id: in2p3-00623484
http://hal.in2p3.fr/in2p3-00623484
Submitted on 20 May 2020
HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entic research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.
L’archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diusion de documents
scientiques de niveau recherche, publiés ou non,
émanant des établissements d’enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
Three-dimensional Lagrangian Turbulent Diusion of
Dust Grains in a Protoplanetary Disk: Method and
First Applications
S. Charnoz, L. Fouchet, J. Aléon, Manuel Moreira
To cite this version:
S. Charnoz, L. Fouchet, J. Aléon, Manuel Moreira. Three-dimensional Lagrangian Turbulent Diusion
of Dust Grains in a Protoplanetary Disk: Method and First Applications. The Astrophysical Journal,
American Astronomical Society, 2011, 737, pp.33. �10.1088/0004-637X/737/1/33�. �in2p3-00623484�

The Astrophysical Journal, 737:33 (17pp), 2011 August 10 doi:10.1088/0004-637X/737/1/33
C
2011. The American Astronomical Society. All rights reserved. Printed in the U.S.A.
THREE-DIMENSIONAL LAGRANGIAN TURBULENT DIFFUSION OF DUST GRAINS IN A
PROTOPLANETARY DISK: METHOD AND FIRST APPLICATIONS
S
´
ebastien Charnoz
1,5
, Laure Fouchet
2
,J
´
er
ˆ
ome Aleon
3
, and Manuel Moreira
4
1
Laboratoire AIM, Universit
´
e Paris Diderot/CEA/CNRS, UMR 7158, 91191 Gif sur Yvette cedex, France; charnoz@cea.fr
2
Physikalisches Institut, Universit
¨
at Bern, CH-3012 Bern, Switzerland
3
Centre de Spectrom
´
etrie Nucl
´
eaire et de Spectrom
´
etrie de Masse, CNRS/IN2P3, Universit
´
e Paris Sud 11, B
ˆ
atiment 104, 91405 Orsay Campus, France
4
Institut de Physique du Globe de Paris, Universit
´
e Paris-Diderot, UMR CNRS 7154, 1 rue Jussieu, 75238 Paris cedex 05, France
Received 2010 December 20; accepted 2011 May 17; published 2011 July 26
ABSTRACT
In order to understand how the chemical and isotopic compositions of dust grains in a gaseous turbulent
protoplanetary disk are altered during their journey in the disk, it is important to determine their individual
trajectories. We study here the dust-diffusive transport using Lagrangian numerical simulations using the popular
“turbulent diffusion” formalism. However, it is naturally expressed in an Eulerian form, which does not allow the
trajectories of individual particles to be studied. We present a simple stochastic and physically justified procedure
for modeling turbulent diffusion in a Lagrangian form that overcomes these difficulties. We show that a net diffusive
flux F of the dust appears and that it is proportional to the gas density (ρ) gradient and the dust diffusion coefficient
D
d
: (F = D
d
× grad(ρ)). It induces an inward transport of dust in the disk’s midplane, while favoring outward
transport in the disk’s upper layers. We present tests and applications comparing dust diffusion in the midplane and
upper layers as well as sample trajectories of particles with different sizes. We also discuss potential applications
for cosmochemistry and smoothed particle hydrodynamic codes.
Key words: methods: numerical protoplanetary disks
Online-only material: color figures
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 pres-
sure. Modeling these processes is important and necessary as
an increasing amount of observational data testifies that, both
in our own solar system and in protoplanetary disks orbit-
ing distant stars, large-scale transport of solids occurs over
tenths of astronomical units (AU). Samples from the comet
81P/Wild2 (Brownlee et al. 2006; Zolensky et al. 2006;
Westphal et al. 2009) and observations of comet 9P/Tempel
1 (Lisse et al. 2006) have revealed the presence of crystalline
silicates, which may have originated within 2 AU of the Sun and
then been transported outward by some mechanism. In this re-
spect, perhaps the most spectacular result of the Stardust mission
to comet Wild2 is the discovery of refractory inclusions formed
at high temperatures (1500 K) in the innermost regions of
the solar protoplanetary disk (1 AU) in a comet originating in
the Kuiper Belt (e.g., Zolensky et al. 2006; Simon et al. 2008;
Matzel et al. 2010). Similar observations made with the Spitzer
Space Telescope have revealed that crystalline silicates are also
found in T Tauri disks with crystalline-to-amorphous silicate
ratios varying greatly from one disk to another (Van Boekel
et al. 2005; Olofsson et al. 2009;Watsonetal.2009). These
results suggest that global dust transport processes are indeed
active, but that their nature and efficiency may differ signifi-
cantly from one disk to another. Since dust settling is expected
to modify the photometric appearance of a disk (Dullemond &
Dominik 2004), observations can reveal the disk structure and
signs of transport. For example, in the GG Tauri circumbinary
disk (Duch
ˆ
ene et al. 2004; Pinte et al. 2007), multi-wavelength
5
Author to whom any correspondence should be addressed.
observations have revealed a radial size sorting of dust particles
at least qualitatively consistent with the radial migration induced
by gas drag.
Many studies have addressed different physical aspects of
dust transport (for radiation pressure, see, e.g., Vinkovi
´
c 2009;
photophoresis, see, e.g., Krauss & Wurm 2005; stellar wind,
see, e.g., Shu et al. 2001; turbulent diffusion, see, e.g., Gail
2001;Ciesla2009; photoevaporation, see, e.g., Alexander &
Armitage 2009). These studies use various tools and formalisms,
either Eulerian or Lagrangian, which are most often designed to
study one specific aspect of transport. For example, multi-fluid
simulations are appropriate for describing the turbulent transport
of fine dust (see, e.g., Fromang & Nelson 2009), but they are very
computationally demanding. At the opposite end, populations
of large particles (which are less collisional) are not accurately
described by a fluid approach and Lagrangian approaches seem
to be more adapted given that they track individual trajectories
(see, e.g., Johansen & Youdin 2007; Johansen et al. 2007).
It would thus be useful to build a three-dimensional (3D)
modular code that allows for a simple coupling of the different
physical processes of dust transport. As a first step, the present
paper describes a 3D Lagrangian dust transport code in which
the equation of motion is numerically solved with as few
approximations as possible. While the motion of small dust
particles that are tightly coupled to gas is analytically well
known, the motion of particles loosely coupled to gas requires
direct integration. Here we focus on the motion of dust particles
under gas drag and turbulent diffusion. Other processes will
be incorporated into future work. Because the effect of gas
drag in a laminar disk using a Lagrangian description is largely
documented (see, e.g., Barri
`
ere-Fouchet et al. 2005), this paper
focuses mainly on the most difficult aspect: the inclusion of
turbulent diffusion in a Lagrangian code. We describe below the
various features of this code.
1

The Astrophysical Journal, 737:33 (17pp), 2011 August 10 Charnoz et al.
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. One very attractive advantage of this ap-
proach is that it affords the possibility of tracking individual
particle trajectories, and thus, of reconstructing the thermody-
namical history of particles in a protoplanetary disk environ-
ment. This should have many applications for cosmochemistry.
For instance, chondrules and refractory inclusions, which are
the main millimeter- to-centimeter-sized components of prim-
itive chondritic meteorites, are thought to have undergone a
complex multi-stage evolution during the first 2–3 Myr of so-
lar system history, before being incorporated into asteroidal or
cometary planetesimals at various heliocentric distances. Chon-
drules (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 ir-
radiation 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, track-
ing 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 so-
lar protoplanetary disk and to unraveling the origin of planetary
volatiles, such as water or organic molecules with prebiotic po-
tential.
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. A direct approach would be to couple
a particle-based dust transport code with a 3D magnetohydro-
dynamic (MHD) simulation of a turbulent gas disk (as in Fro-
mang & 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 rea-
son, 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 environmen-
tal 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 turbu-
lent 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 differ-
ence 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 exam-
ples in Section 4.2). We will see, in particular, that proper im-
plementation 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 Thermo-
dynamics. 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.
1.3. Three-dimensional System
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;Ciesla2007; Tscharnuter
&Gail2007) 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.
2

The Astrophysical Journal, 737:33 (17pp), 2011 August 10 Charnoz et al.
2. NUMERICAL IMPLEMENTATION OF
TURBULENT DIFFUSION
2.1. Basic Concepts and Definitions
In a Lagrangian code, the motion of an individual dust particle
is described using the classical Newtonian formalism:
d v
dt
=
F
m
v −v
g
τ
, (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:
τ =
s
ρC
s
, (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
3
)
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)isrewritten:
d v
dt
=
F
m
v −v
g
τ
+
δv
g
τ
, (3)
where δv
g
is the acceleration of dust induced by the turbulent
gas velocity field that induces a random walk (or turbulent dif-
fusion) of the particle, consistent with the “turbulent diffusion”
model. We assume that v
g
can be identified with the unper-
turbed 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 diffu-
sion 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):
δv = δ v
M
+ δ v
T
(4)
δr = δr
M
+ δr
T
, (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.Thetermδr
M
is obtained by direct
numerical integration of the equation:
δr
M
=
t+dt
t
v(t) dt (6)
with
d v
dt
=
F
m
v −v
g
τ
. (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
7
and 10
3
, 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
8
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.
2.2. The “Good Mixing” Problem
A diffusion process is usually described by Fick’s law, which
relates the diffusive flux of material F
T
to the density gradient:
F
T
=−D ·
grad(ρ), (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
.Tomimica
3

The Astrophysical Journal, 737:33 (17pp), 2011 August 10 Charnoz et al.
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 σ
2
r
: the variance of δv
T
as σ
2
v
,
δr
T
and δr
T
being their respective mean. In the “position”
representation, the kick on position is a random Gaussian with
δr
T
=
δr
T
=0
σ
2
r
= 2Ddt
. (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., (σ
2
r
)
1/2
= dt × (σ
2
V
)
1/2
):
δv
T
=
δv
T
=0
σ
2
v
=
2D
dt
. (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
s
H. 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) =
ρ
0
exp(z
2
/2H
2
), 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.
2.3. Good Diffusion with a Constant Diffusion Coefficient
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
dX = Vdt +(2Ddt)
1/2
W, (11)
where W is a Gaussian random variable with mean 0 and
variance σ
2
= 1 (such that (2Ddt)
1/2
× W is a random
variable with mean 0 and variance 2Ddt, as in Equation (9)). A
basic result of Einstein–Brownian diffusion is that the resulting
motion has an average position X so that dX/dt = V, and
that the standard deviation (X−X)
2
evolves according to
d/dt((X−X)
2
) = 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:
∂ρ
∂t
+
∂ρV
∂x
=
∂x
D
∂ρ
d
∂x
, (12)
also written as
∂ρ
∂t
+
∂x
ρV −D
∂ρ
∂x
= 0. (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
4

Figures (11)
Citations
More filters

Journal ArticleDOI
Abstract: The Disk Substructures at High Angular Resolution Project (DSHARP) provides a large sample of protoplanetary disks with substructures that could be induced by young forming planets. To explore the properties of planets that may be responsible for these substructures, we systematically carry out a grid of 2D hydrodynamical simulations, including both gas and dust components. We present the resulting gas structures, including the relationship between the planet mass, as well as (1) the gaseous gap depth/width and (2) the sub/super-Keplerian motion across the gap. We then compute dust continuum intensity maps at the frequency of the DSHARP observations. We provide the relationship between the planet mass, as well as (1) the depth/width of the gaps at millimeter intensity maps, (2) the gap edge ellipticity and asymmetry, and (3) the position of secondary gaps induced by the planet. With these relationships, we lay out the procedure to constrain the planet mass using gap properties, and study the potential planets in the DSHARP disks. We highlight the excellent agreement between observations and simulations for AS 209 and the detectability of the young solar system analog. Finally, under the assumption that the detected gaps are induced by young planets, we characterize the young planet population in the planet mass-semimajor axis diagram. We find that the occurrence rate for >5 M-J planets beyond 5-10 au is consistent with direct imaging constraints. Disk substructures allow us to probe a wide-orbit planet population (Neptune to Jupiter mass planets beyond 10 au) that is not accessible to other planet searching techniques.

256 citations


Journal ArticleDOI
27 Apr 2012-Science
TL;DR: It is found that icy grains originating in the outer disk, where temperatures were less than 30 kelvin, experienced ultraviolet irradiation exposures and thermal warming similar to that which has been shown to produce complex organics in laboratory experiments, implying that organic compounds are natural by-products of protoplanetary disk evolution and should be important ingredients in the formation of all planetary systems, including the authors' own.
Abstract: Complex organic compounds, including many important to life on Earth, are commonly found in meteoritic and cometary samples, though their origins remain a mystery. We examined whether such molecules could be produced within the solar nebula by tracking the dynamical evolution of ice grains in the nebula and recording the environments to which they were exposed. We found that icy grains originating in the outer disk, where temperatures were less than 30 kelvin, experienced ultraviolet irradiation exposures and thermal warming similar to that which has been shown to produce complex organics in laboratory experiments. These results imply that organic compounds are natural by-products of protoplanetary disk evolution and should be important ingredients in the formation of all planetary systems, including our own.

183 citations


Journal ArticleDOI
Abstract: The Atacama Large Millimeter Array (ALMA) has been returning images of transitional disks in which large asymmetries are seen in the distribution of mm-sized dust in the outer disk. The explanation in vogue borrows from the vortex literature by suggesting that these asymmetries are the result of dust trapping in giant vortices, excited via Rossby wave instability (RWI) at planetary gap edges. Due to the drag force, dust trapped in vortices will accumulate in the center, and diffusion is needed to maintain a steady state over the lifetime of the disk. While previous work derived semi-analytical models of the process, in this paper we provide analytical steady-state solutions. Exact solutions exist for certain vortex models. The solution is determined by the vortex rotation profile, the gas scale height, the vortex aspect ratio, and the ratio of dust diffusion to gas-dust friction. In principle, all these quantities can be derived from observations, which would give validation of the model, also giving constrains on the strength of the turbulence inside the vortex core. Based on our solution, we derive quantities such as the gas-dust contrast, the trapped dust mass, and the dust contrast at the same orbital location. We apply our model to the recently imaged Oph IRS 48 system, finding values within the range of the observational uncertainties.

157 citations


Journal ArticleDOI
Abstract: The Atacama Large Millimeter Array has returned images of transitional disks in which large asymmetries are seen in the distribution of millimeter sized dust in the outer disk. The explanation in vogue borrows from the vortex literature and suggests that these asymmetries are the result of dust trapping in giant vortices, excited via Rossby wave instabilities at planetary gap edges. Due to the drag force, dust trapped in vortices will accumulate in the center and diffusion is needed to maintain a steady state over the lifetime of the disk. While previous work derived semi-analytical models of the process, in this paper we provide analytical steady-steady solutions. Exact solutions exist for certain vortex models. The solution is determined by the vortex rotation profile, the gas scale height, the vortex aspect ratio, and the ratio of dust diffusion to gas-dust friction. In principle, all of these quantities can be derived from observations, which would validate the model and also provide constrains on the strength of the turbulence inside the vortex core. Based on our solution, we derive quantities such as the gas-dust contrast, the trapped dust mass, and the dust contrast at the same orbital location. We apply our model to the recently imaged Oph IRS 48 system, finding values within the range of the observational uncertainties.

147 citations


Journal ArticleDOI
Abstract: Context. When a planet becomes massive enough, it gradually carves a partial gap around its orbit in the protoplanetary disk. A pressure maximum can be formed outside the gap where solids that are loosely coupled to the gas, typically in the pebble size range, can be trapped. The minimum planet mass for building such a trap, which is called the pebble isolation mass (PIM), is important for two reasons: it marks the end of planetary growth by pebble accretion, and the trapped dust forms a ring that may be observed with millimetre observations.Aims. We study the effect of disk turbulence on the PIM and find its dependence on the gas turbulent viscosity, aspect ratio, and particles Stokes number.Methods. By means of 2D gas hydrodynamical simulations, we found the minimum planet mass to form a radial pressure maximum beyond the orbit of the planet, which is the necessary condition to trap pebbles. We then carried out 2D gas plus dust hydrodynamical simulations to examine how dust turbulent diffusion impacts particles trapping at the pressure maximum. We finally provide a semi-analytical calculation of the PIM based on comparing the radial drift velocity of solids and the root mean square turbulent velocity fluctuations around the pressure maximum.Results. From our results of gas simulations, we provide an expression for the PIM vs. disk aspect ratio and turbulent viscosity. Our gas plus dust simulations show that the effective PIM can be nearly an order of magnitude larger in high-viscosity disks because turbulence diffuse particles out of the pressure maximum. This is quantified by our semi-analytical calculation, which gives an explicit dependence of the PIM with Stokes number of particles.Conclusions. Disk turbulence can significantly alter the PIM, depending on the level of turbulence in regions of planet formation.

91 citations


Frequently Asked Questions (2)
Q1. What are the contributions mentioned in the paper "Three-dimensional lagrangian turbulent diffusion of dust grains in a protoplanetary disk: method and first applications" ?

The authors study here the dust-diffusive transport using Lagrangian numerical simulations using the popular “ turbulent diffusion ” formalism. However, it is naturally expressed in an Eulerian form, which does not allow the trajectories of individual particles to be studied. The authors present a simple stochastic and physically justified procedure for modeling turbulent diffusion in a Lagrangian form that overcomes these difficulties. The authors show that a net diffusive flux F of the dust appears and that it is proportional to the gas density ( ρ ) gradient and the dust diffusion coefficient Dd: ( F = Dd/ρ × grad ( ρ ) ). The authors present tests and applications comparing dust diffusion in the midplane and upper layers as well as sample trajectories of particles with different sizes. The authors also discuss potential applications for cosmochemistry and smoothed particle hydrodynamic codes. 

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.