EDP Sciences
Free access
Volume 502, Number 1, July IV 2009
Page(s) 385 - 389
Section Planets and planetary systems
DOI http://dx.doi.org/10.1051/0004-6361/200911865
Published online 15 June 2009

On the pressure of collisionless particle fluids
(Research Note)

The case of solids settling in disks

F. Hersant1,2

1 - Université de Bordeaux, Laboratoire d'Astrophysique de Bordeaux (LAB), 2 rue de l'Observatoire, BP 89, 33271 Floirac Cedex, France
2 - CNRS/INSU - UMR5804, BP 89, 33271 Floirac Cedex, France

Received 17 February 2009 / Accepted 4 May 2009

Aims. Collections of dust, grains, and planetesimals are often treated as a pressureless fluid. We study the validity of neglecting the pressure of such a fluid by computing it exactly for the case of particles settling in a disk.
Methods. We solve a modified collisionless Boltzmann equation for the particles and compute the corresponding moments of the phase space distribution: density, momentum, and pressure.
Results. We find that whenever the Stokes number, defined as the ratio of the gas drag timescale to the orbital timescale, is more than 1/2, the particle fluid cannot be considered as pressureless. While we show it only in the simple case of particles settling in a laminar disk, this property is likely to remain true for most flows, including turbulent flows.

Key words: hydrodynamics - methods: analytical - planetary systems: formation

1 Introduction

The interaction of solid particles with gas in protoplanetary disks is one of the key mechanisms for planetesimals and planet formation. Moreover, this fundamental process gets more and more attention for other astrophysical media (interstellar medium, planetary atmospheres, etc.) as grains can dampen fluid motions and be the targets of adsorption and condensation of molecules and can be the base for complex chemistry, e.g., the formation of the molecule H2 in the interstellar medium. With the advance of observational facilities and techniques, dust dynamics is understood more and more as a crucial process for most media, while, due to the discrete nature of grains, modelling their dynamics remains tricky.

There are two main methods for computing the dynamics of particles embedded in a fluid. One can either study a large number of individual particles using N-body techniques or consider the particles as a fluid, which is treated with a modified Navier-Stokes/Euler equation. While the first method is reliable and precise, it can only be applied to a finite (and in fact limited) number of particles; therefore, macroscopic quantities reconstructed from the N-body simulations outputs are sometimes inaccurate due to a low sampling of the phase space. The second method enables one to consider an almost infinite ensemble of particles by computing their macroscopic properties (density, velocity field, etc.). This makes the ``two-fluid'' approach (gas and solids) the preferred method of treating the dynamics and growth of solid particles in a fluid, especially in the planet formation community. Consequently, a dense and rich literature is based on this approach, and many fundamental mechanisms of planetary formation have been unraveled (e.g. Inaba et al. 2005; Stepinski & Valageas 1996; Fromang & Papaloizou 2006; Barrière-Fouchet et al. 2005; Cuzzi et al. 1993; Johansen et al. 2006; Dubrulle et al. 1995).

The two-fluid approach is usually implemented by introducing two Euler equations, coupled through the gas drag force. Moreover, in the fluid equation for particles, the pressure term is usually ignored. This approximation is often justified by the statement that collisions between particles and gas are much more probable than collisions between particles making the pressure term negligible compared to the gas drag force.

The purpose of the present paper is to question the validity of ignoring the pressure of the particle fluid by computing it in the very simple case of particles settling in a disk. The vertical dynamics of a solid body in a disk can be understood as a damped oscillator. In a Keplerian disk, if gas is absent or very tenuous, the particle oscillates around the midplane at a frequency equal to the orbital frequency (the oscillation comes from inclined orbits). The friction with the gas damps this oscillation, as the particle feels a ``headwind'' when it crosses the gas disk. Like any damped oscillator, there are two regimes: an underdamped regime for low damping and an overdamped regime, where the particle has no time to oscillate for a full period. This occurs for high damping rates.

We can compute the pressure of the solid fluid embedded in a laminar gas disk by solving exactly the collisionless Boltzmann (or Vlasov) equation for the phase space density of particles. As we show, as soon as the particles are larger than some critical size, the pressure cannot be neglected. This work extends the work of Garaud et al. (2004), who comes to similar conclusions based on asymptotic calculations.

In Sect. 2 is developed the general formalism. Section 3 is devoted to a very simple test case that can be used to check the validity of more general solutions and intuit how the pressure varies. In Sect. 4, the exact general solution of the Vlasov equation is given. In Sect. 5, this solution is applied to a special case that retains most of the expected general properties. A surprisingly simple solution is found for the pressure and the scaling laws, and asymptotics are discussed. In Sect. 6, we discuss our results and their implications for more complex and turbulent flows and give our conclusions.

2 General formalism

We consider the equation for the evolution of the phase space density in 1+1 dimensions (one in space, one in velocity) under the action of a drag force from an inert fluid. This is not exactly the common collisionless Boltzmann or Vlasov equation, as the drag force is not conservative. Instead, the general Vlasov equation with non-conservative forces takes the form (see e.g. Youdin & Lithwick 2007; Lutsko 1997; Carballido et al. 2006):

$\displaystyle \frac{\partial f}{\partial t}+u\frac{\partial f}{\partial z} + \frac{\partial}{\partial
u}\left(\frac{F}{m}f\right)=0$     (1)

where z is the space coordinate, u is the velocity coordinate, F/m the non-conservative force per mass (acceleration), and f the phase space density (or probability density function). When f is known, the macroscopic quantities (density, velocity, energy) are moments of the distribution f. In particular, the zeroth order moment is the density n(z,t):

\begin{displaymath}n(z,t)=\int_{-\infty}^{+\infty} f(z,u,t) {\rm d}u,
\end{displaymath} (2)

the first-order moment is the momentum (the product of the density by the macroscopic velocity $\tilde{u}$)

\begin{displaymath}n(z,t)\ \tilde{u}(z,t)=\int_{-\infty}^{+\infty} u f(z,u,t) {\rm d}u,
\end{displaymath} (3)

and the second-order moment is related to the pressure (a scalar in 1D):
                      P(z,t) = $\displaystyle \int_{-\infty}^{+\infty} (u-\tilde{u})^2 f(z,u,t) {\rm d}u$  
  = $\displaystyle \int_{-\infty}^{+\infty} u^2 f(z,u,t)
{\rm d}u - n\tilde{u}^2.$ (4)

In the case of particles settling in a disk, there are two forces: gravity and gas drag. When the disk is geometrically thin, the vertical gravity can be written as $\frac{F}{m}=-\omega^2 z$, where $\omega$ is the orbital frequency, while the drag force, in the absence of gas motions, can be written as $\frac{F}{m}=-u/\tau$ where $\tau$ is the so-called stopping time (Whipple 1972). This is the characteristic timescale for a particle to stop in a fluid. Here, we consider $\tau$ as a free parameter, which would in general depend upon the particle size, morphology, density, the density of the gas, and the turbulent intensity in the grain trail. Here, we consider it as independent of both space and velocity. This is an important simplification, corresponding to the so-called Epstein regime, where grains are smaller than the mean free path of the underlying gas. It should, however, be noted that the mean free path $l_{\rm mfp}$ in a circumstellar disk is very large, typically varying like $l_{\rm mfp}=10^8 \left(\frac{R}{100\ {\rm AU}}\right)^{5/2}$ cm, assuming that the gas number density $n_{\rm g}$ varies like $n_{\rm g}=10^8 \left(\frac{R}{100\ {\rm AU}}\right)^{-5/2}$ cm-3(Hersant et al. 2005; Guilloteau & Dutrey 1998).

Under these conditions, Eq. (1) becomes

\begin{displaymath}\partial_t f +u\partial_z f - \partial_u \left( \omega^2 z +\frac{u}{\tau} \right) f=0
\end{displaymath} (5)

where z is now the vertical space variable and u the vertical velocity. To simplify the notations, we use the convention $\partial_x=\frac{\partial}{\partial x}$. This equation can be developed into pure advective terms:

\begin{displaymath}\partial_t f +u\partial_z f - \left( \omega^2 z +\frac{u}{\tau} \right) \partial_u f -\frac{f}{\tau}=0.
\end{displaymath} (6)

The last term on the lefthand side ($-f/\tau$) is the difference with a Vlasov equation for conservative forces and deserves some comments. This term induces an exponential growth of f. Without solving the equation, one can understand that the conservation of the number of particles (following from Eq. (1)), requires that the surface of the phase space field of density f has to collapse in one direction so that, at infinite times, the field of f tends to a discontinuous distribution, consisting of a family of curves, but nothing like a continuous surface. This property of the equation remains valid even in 3+3 dimensions and with gas turbulence.

At this point, it is convenient to adimension the time, introducing two new variables $t'=\omega t$ and $u'=u/\omega$. Using these new variables and omitting the primes for clarity, we get

\begin{displaymath}\partial_t f +u\partial_z f - \left(z +\frac{u}{S} \right) \partial_u f -\frac{f}{S}=0
\end{displaymath} (7)

where $S=\omega \tau$ is a Stokes number of the problem. For small (large) particles, S is small (large) due to the strong (weak) coupling of the particles with the fluid.

We now can artificially get rid of the exponential growth term by defining

\begin{displaymath}\psi=f\ {\rm e}^{-t/S}
\end{displaymath} (8)

and finally get

\begin{displaymath}\partial_t \psi +u\partial_z \psi - \left( z +\frac{u}{S} \right)\partial_u \psi=0.
\end{displaymath} (9)

The behavior of this equation is rather easy to foresee. For infinite S, the remaining terms correspond to a rotation in phase space. This rotation is a direct consequence of energy conservation. In real space, particles oscillate around the midplane with zero velocity at the highest absolute altitude and maximum velocity when particles cross the midplane. The gas-drag term induces a collapse along the velocity axis. The net outcome of both the rotation and the collapse along the velocity axis (which corresponds in the real space to the dampening by gas drag of the inclination of the particle orbit) is a global collapse towards the center of the phase space (z=0, u=0), but not a spherical collapse.

3 Simple case

Although the general solution is given in the next section, it is useful to solve Eq. (9) in the simplest case: infinitely large or massive particles ($S=+\infty$). In this case, the equation admits a simple stationary solution with separated variables:

\begin{displaymath}\psi=K{\rm e}^{-\frac{z^2+u^2}{2\sigma^2}}.
\end{displaymath} (10)

Using the definition of the moments, the pressure is straightforward:

\begin{displaymath}P(z)=n(z) \sigma^2,
\end{displaymath} (11)

where the density is

\begin{displaymath}n(z)=K\sigma \sqrt{2\pi}{\rm e}^{-\frac{z^2}{2\sigma^2}}.
\end{displaymath} (12)

Equation (11) is similar to an equation of state for isothermal fluids, where $\sigma$ plays the role of a (constant) sound speed. As we will see in Sect. 5, this behavior remains valid even in more general cases. With the equation of state (11), solution (12) is the result of the hydrostatic equilibrium of an isothermal gas, the pressure force balancing the vertical gravity. The stationary solution (10) being even in u, it corresponds to a zero macroscopic velocity. In this case, where particles have an infinite inertia, the pressure term is thus dominant in the structure of the particle disk.

4 General solution

Equation (9) can be solved in the general case, using standard techniques (the method of characteristics for example). We will not develop this here as it is easier to check that the following solution is the general solution of Eq. (9):

\end{displaymath} (13)

where $\phi$ is the initial condition and (z0, u0) is the initial phase space position of a particle located in (z, u) at a time t. The initial positions are given by

z_0 = f_1(t) z\ {\rm e}^{\alpha t} ...
...^{\alpha t} + f_3(t) u\ {\rm e}^{\alpha t},
\end{displaymath} (14)

where f1, f2 and f3 are time dependent functions given by

f_1(t)={\rm ch}(\beta t)-\frac{\alp...
...ta t)+\frac{\alpha}{\beta}{\rm sh}(\beta t)
\end{displaymath} (15)

and $\alpha$ and $\beta $ are parameters depending on the Stokes number S:

\alpha=\frac{1}{2S}\\ [2mm]
\end{displaymath} (16)

Here, $\alpha$ is always real and positive, while $\beta $ is either a positive real number (for small Stokes numbers), or a purely imaginary number (for large Stokes numbers). However, whatever $\beta $, the time functions f1, f2, and f3 are always real.

5 Particular case for the initial conditions

Computing the moments of the phase space distribution cannot be done without the initial conditions. Here, we choose an initial condition composed of Gaussians, with one Gaussian per dimension, which is relevant to most physically meaningful cases:

\begin{displaymath}\phi(z,u)=K{\rm e}^{-\frac{z^2}{2\sigma_z^2}}{\rm e}^{-\frac{u^2}{2\sigma_u^2}}
\end{displaymath} (17)

where $\sigma_z$ is the initial thickness of the particle sub-disk and $\sigma_u$ is the initial (presumably small) velocity dispersion. The normalization constant K is related to the surface density $\Sigma$ by $\Sigma=2\pi K \sigma_z
\sigma_u$. This initial condition assumes that the particles initially have a vanishing mean velocity, a velocity dispersion and are randomly distributed vertically, following a Gaussian distribution in this case. Other choices would lead to similar results unless the symmetry with respect to the disk midplane is broken (particles initially located only in the upper half of the disk, for example).

5.1 The corresponding moments

Using definition (2), the density corresponding to the initial conditions (17) can be written as

\begin{displaymath}n(z,t)=K{\rm e}^{2\alpha
t}{\rm exp}\left(-\frac{z^2f_1^2(t)}...
...rm exp}\left(-\frac{z^2f_2^2(t)}{2\sigma_u'^2(t)}\right)L(z,t)
\end{displaymath} (18)

where the exponential dependencies can be included into the Gaussian widths:

\sigma_z'=\sigma_z {\rm e}^{-\alpha...
\sigma_u'=\sigma_u {\rm e}^{-\alpha t}
\end{displaymath} (19)

and L(z,t) is a function defined as

{\rm exp}\left(-\frac{u^2f_2^...
...frac{u^2f_3^2(t)}{2\sigma_u'^2(t)}\right){\rm e}^{-pu}{\rm d}u
\end{displaymath} (20)


\end{displaymath} (21)

Let us define

\end{displaymath} (22)


\begin{displaymath}L_p=\int_{-\infty}^{+\infty}{\rm e}^{-\frac{u^2}{2\lambda^2}}{\rm e}^{-pu}{\rm d}u.
\end{displaymath} (23)

Here we keep the p-dependence explicitly as it will be useful in the following. We have L(z,t)=Lp if ptakes the value given by Eq. (21). This expression can be now written as

Lp=F(p)+F(-p) (24)

where F(p) is the Laplace transform of a Gaussian:

\begin{displaymath}F(p)=\int_{0}^{+\infty}{\rm e}^{-\frac{u^2}{2\lambda^2}}{\rm ...
...p^2}{2}}{\rm erfc}\left(\frac{\lambda
\end{displaymath} (25)

Therefore, using the identity ${\rm erfc}(x)+{\rm erfc}(-x)=2$, we have

\begin{displaymath}L_p=\lambda\sqrt{2\pi}{\rm e}^{\frac{\lambda^2p^2}{2}}.
\end{displaymath} (26)

Using the definition (3) and a similar procedure, we get for the first-order moment:

\begin{displaymath}n\tilde{u}=K{\rm e}^{2\alpha
t}{\rm exp}\left(-\frac{z^2f_1^2...
... exp}\left(-\frac{z^2f_2^2(t)}{2\sigma_u'^2(t)}\right)H(z,t),
\end{displaymath} (27)


\begin{displaymath}H(z,t)=\int_{-\infty}^{+\infty} u\ {\rm e}^{-\frac{u^2}{2\lambda^2}}{\rm e}^{-pu}{\rm d}u.
\end{displaymath} (28)

Now Hp can be written as a function of Fp, using basic properties of the Laplace transform:

\begin{displaymath}H_p=\partial_p F(-p)-\partial_p F(p)=-\partial_p L_p.
\end{displaymath} (29)


\begin{displaymath}\tilde{u}=-\frac{\partial_p L_p}{L_p}\cdot
\end{displaymath} (30)

Applying the same procedure to definition (4), we get the following expression for the particle fluid pressure:

\end{displaymath} (31)

and using Eq. (26), we finally get the simple expression for the pressure:

\end{displaymath} (32)

Here $\lambda$ can be understood as the velocity dispersion of the particles, and its expression 22 generalizes the asymptotics of Garaud et al. (2004).

It is interesting to note that the general pressure only depends on z through the density. The Euler equation for the particle fluid involves only gradients of the pressure so this gradient is in fact directly proportional to the density gradient. In other words, the particles seem to follow a barotropic equation of state. This behavior seems closely related to the assumption of constant stopping time (Garaud et al. 2004).

5.2 Scaling and asymptotics

5.2.1 Case $\beta $ real (small particles)

When $\beta $ is real (and positive), $\alpha$ is large, it is easy to check that $\lambda$ is a decreasing function of time. Its initial value is

\end{displaymath} (33)

and the large time asymptotics are

\begin{displaymath}\lambda(t \to +\infty) \sim \frac{\beta
...ght)^2\left(\alpha+\beta\right)^2}}{\rm e}^{-(\alpha+\beta)t}.
\end{displaymath} (34)

Starting from a presumably very low value (the initial velocity dispersion of particles), $\lambda$ decreases with time, asymptotically like an exponential. The location where the pressure decreases the least is the midplane, since the density tends to a Dirac function when times goes to infinity,

\begin{displaymath}n(z=0,t)=K \lambda\sqrt{2\pi}{\rm e}^{2\alpha t},
\end{displaymath} (35)

$\displaystyle P(z=0,t)=\lambda^3 K \sqrt{2\pi} {\rm e}^{2\alpha t}
{\rm e}^{-(\alpha+3\beta)t}.$     (36)

Even in the midplane, the pressure quickly drops with time, starting from a low value. In this case, as has long been known (see e.g. Garaud et al. 2004; Cuzzi et al. 1993; Dubrulle et al. 1995), it is safe to neglect the pressure of the particle fluid. Consequently, particles steadily sediment towards the disk midplane at a velocity close to the terminal velocity $w^\infty=\omega^2z\tau$ (defined as the asymptotic velocity of any body free-falling in a fluid), which is the asymptotic velocity of both the particle fluid and individual particles here.

5.2.2 Case $\beta $ imaginary (large particles)

We define $\gamma$ such that $\beta={\rm i} \gamma$. In such a case, $\lambda$ can be written as

\begin{displaymath}\lambda={\rm e}^{-\alpha t}\frac{\sigma_u}{\sqrt{\left( {\rm ...
...igma_z}\right)^2\frac{1}{\gamma^2}{\rm sin}^2\ \gamma t}}\cdot
\end{displaymath} (37)

For small initial velocity dispersions, the location of the maxima of $\lambda$ corresponds to

\begin{displaymath}{\rm tan}\ \gamma t=-\frac{\alpha}{\gamma}
\end{displaymath} (38)

and the amplitude of these maxima is

\begin{displaymath}\lambda_{\max}(t)=\sigma_z \frac{\gamma}{\alpha}\sqrt{\alpha^2+\gamma^2}\ {\rm e}^{-\alpha t}.
\end{displaymath} (39)

It is crucial to note here that the amplitude of the maxima of the time-dependent velocity dispersion $\lambda$ is independent of the initial velocity dispersion $\sigma_u$. The velocity dispersion is created and maintained by the vertical gravity.

Since the midplane density is still given by Eq. (35), we get the amplitude of the midplane pressure maxima:

\begin{displaymath}P_{\max}(t)=K\sqrt{2 \pi}\left(\sigma_z \frac{\gamma}{\alpha}\sqrt{\alpha^2+\gamma^2}
\right)^3 {\rm e}^{-\alpha t}.
\end{displaymath} (40)

These pressure maxima have a much lower decreasing rate than in the non-oscillating case, both because it involves only $\alpha$ here and because $\alpha$ has a lower value.

The midplane pressure as a function of time for different values of the Stokes number S is displayed in Fig. 1. The transition between the case with oscillation S > 1/2 and the overdamped case S < 1/2 appears very clearly. The pressure is computed for a ratio $\frac{\sigma_u}{\sigma_z}=0.1$. This ratio defines the pressure at t=0 and the width of the maxima, but neither the amplitude of the maxima nor the decreasing rate.

\end{figure} Figure 1:

Pressure of the particle fluid in the midplane as a function of the dimensionless time ($\omega t$), for different values of the Stokes number S

Open with DEXTER

Figure 1 illustrates that, while the pressure of the solids fluid can be safely neglected for small values of the Stokes number S, the situation is drastically different when the Stokes number reaches S=1/2. When particles oscillate around the midplane, the fluid macroscopic velocity is different from the terminal velocity $w^\infty$ as the macroscopic pressure force slows down the vertical contraction of the particle fluid. In other words, the velocity of the particle fluid is not that of individual particles and the fluid has a nonvanishing internal energy (contrary to the case for small particles).

6 Discussion and conclusions

We have shown that, in general, the pressure of the particle fluid cannot be neglected as soon as their Stokes number is more than 0.5. The particles seen as a fluid can have an internal energy, and the pressure force can partly balance gravity. In the extreme case where particles have an infinite inertia, the structure of the particle fluid is in fact hydrostatic. The existence of pressure is completely independent of the presence or absence of collisions between solid particles. Instead, the pressure is the result of a velocity dispersion. While velocity dispersion can be the result of collisions, it can also be created dynamically when external forces excite the velocity dispersion (gravity in the present case) or when the flow itself is converging locally. Even though our results concern only a very simple and specific problem, they have implications for a much broader variety of astrophysical and non-astrophysical flows.

The Stokes number can be defined in a more general way as the ratio of a dynamical timescale, creating velocity differences between the solids and the underlying fluid (gas or liquid), and the stopping time, the characteristic timescale for the particles to reach the fluid velocity. There are many processes able to create velocity differences between the gas and solids, gravity being only one of them. For example, when two rivers join, particles embedded in the rivers do not see the confluence like the water does and this can create a velocity difference between water and the particles, and a velocity dispersion in the particle flow at the confluence.

For particles embedded in a turbulent flow, it is common to define a scale-dependent Stokes number by $S_\eta=
\tau/t_\eta$, where $t_\eta$ is the turnover time of a turbulent eddy of size $\eta$. This Stokes number depends on the size of the turbulent eddy and scales as $S_\eta \varpropto \eta^{2/3}$ for Kolmogorov-like turbulence (see e.g. Cuzzi et al. 2001). The smallest reachable eddy is dependent on the strength of turbulence, quantified by the Reynolds number of the flow $Re=\frac{UL}{\nu}$, where U and L are characteristic velocities and length scale, respectively, and $\nu$ is the molecular kinematic viscosity of the fluid. For Kolmogorov turbulence, the smallest turbulent scale $\eta_{\min}$ of the flow scales as $\eta_{\min} \varpropto
R{\rm e}^{-3/4}$. For flows with large enough Reynolds numbers (many astrophysical and geophysical flows reach Reynolds numbers over 1010), the smallest turbulent scales have $S_\eta$ smaller than 1, even for small grains. This is a well-known behavior of particle-laden flows. Particles are centrifuged out of the turbulent eddies whentheir scale Stokes number reaches unity (Squires & Eaton 1991). This creates a converging flow of particles towards regions of low vorticity. Consequently, this induces a velocity dispersion inside the particle flow and, for similar reasons to those described in this paper, the particle fluid develops a pressure. As a side note, when turbulence modelling is considered, turbulent pressure terms can arise from fluctuating or subgrid velocity dispersions (Dobrovolskis et al. 1999). However, turbulent pressure terms are a consequence of Reynolds averaging or spatial filtering and disappear when turbulence is directly simulated instead of modelled.

This is hopefully a good illustration of the care one has to take whenever a fluid equation is applied to a collection of solid particles in an astrophysical flow. Whenever any dynamical timescale becomes smaller than the characteristic coupling timescale between the particles and the underlying fluid, the particles cannot be considered as a pressureless fluid, whatever the collision rate between particles.

Part of this work was done during the author's stay at Jeremiah Horrocks Institute for Astrophysics & Supercomputing, University of Central Lancashire, whose hospitality is gratefully acknowledged. I am grateful to J. Braine, S. Courty, B. Dubrulle, and J.-M. Huré for comments on the manuscript, and to the anonymous referee and T. Guillot for suggestions that helped me to significantly improve the presentation of the paper.


All Figures

\end{figure} Figure 1:

Pressure of the particle fluid in the midplane as a function of the dimensionless time ($\omega t$), for different values of the Stokes number S

Open with DEXTER
In the text

Copyright ESO 2009