# The dynamics of discrete particles in turbulent flows: open issues and
current challenges in statistical modeling

Jean-Pierre Minier<EMAIL_ADDRESS>EDF R$\&$D, Mécanique des
Fluides, Energie et Environnement, 6 quai Watier, 78400 Chatou, France
Christophe Henry<EMAIL_ADDRESS>Université Côte d’Azur, Inria,
CNRS, Sophia-Antipolis, France

###### Abstract

This article is an invitation. It is, first, an invitation to consider as a
subject worthy of attention the wide range of situations where small discrete
elements, either bubbles, droplets or solid particles, are embedded in
turbulent flows. Occurring often at a human scale and in our daily
environments, these turbulent dispersed two-phase flows display complex
behavior due to the interplay of two fundamental interactions, the fluid-
particle and particle-particle interactions, compounded by the turbulence of
the carrier flow. This is not a domain where the basic laws are unknown but
where the huge number of degrees of freedom involved call for reduced, or
coarse-grained, statistical descriptions to be developed. Since we are
considering transport and collision phenomena or relaxation processes, it
would seem that they can be handled by kinetic theory. In the general case of
non-fully resolved turbulent flows, we are however dealing with particles
influenced by random media with non-zero time and space correlations. The
second invitation is therefore to recognize the limitations of kinetic-based
descriptions and to address the challenges driving us to extend the classical
framework, for fluid-particle as well as particle-particle interactions.
Taking the standpoint provided by the modern formulation of stochastic
processes and focusing on the description of the particle phase, this review
proposes a step-by-step pedagogical presentation of current models while
pointing out new directions and remaining uncharted territories. This is done
to provide answers to the question ‘why?’ as much as ‘how?’ and to try to
kindle interest into these open and fascinating issues.

††preprint: Proposal submitted to RMP

###### Contents

  1. I The physics of dispersed turbulent two-phase flows
    1. I.1 A rich tapestry of applications
      1. I.1.1 Colloid suspensions and river deltas
      2. I.1.2 Plume dispersion and bubbly flows
      3. I.1.3 The physics of sand dunes
    2. I.2 An intricate range of phenomena
      1. I.2.1 The two fundamental interactions: particle-fluid and particle-particle
      2. I.2.2 Addressing complexity
    3. I.3 The governing equations
      1. I.3.1 The Navier-Stokes equations for fluid flows
      2. I.3.2 The equations of motion for particles
    4. I.4 The issue to address in statistical physics
  2. II Reduced statistical descriptions and the probabilistic framework
    1. II.1 From microscopic to macroscopic descriptions
      1. II.1.1 Mass density functions and mean fields
      2. II.1.2 A unified approach to statistical operators
    2. II.2 Stochastic processes
      1. II.2.1 Definition and key characteristics
      2. II.2.2 Markovian processes
      3. II.2.3 The two building blocks: the Wiener and Poisson processes
    3. II.3 Jump-diffusion processes and extended Fokker-Planck equations
    4. II.4 A double hierarchy: one- or N-particle PDF and particle-attached variables
  3. III Statistical modeling of particle transport in turbulent flows
    1. III.1 Shortcomings of kinetic-like descriptions
    2. III.2 Slow and fast variables and the Markovian approach
      1. III.2.1 Colored or white noises and well-posed PDF equations
      2. III.2.2 The Markovian approach
    3. III.3 The Kolmogorov theory and Lagrangian models
  4. IV Modeling the velocity of the fluid seen
    1. IV.1 The physics of particle dispersion
      1. IV.1.1 The crossing-trajectory effect
      2. IV.1.2 Stochastic diffusion models for the velocity of the fluid seen
      3. IV.1.3 Remarks on the particle preferential concentration effect
    2. IV.2 Stochastic models for fluid particles
    3. IV.3 Langevin models for the velocity of the fluid seen
      1. IV.3.1 State-of-the-art two-phase SLM
      2. IV.3.2 Applications and extensions in LES
      3. IV.3.3 Proposals for general two-phase GLM
      4. IV.3.4 Consistency issues and limitations
    4. IV.4 New macroscopic approaches and microscopic models
      1. IV.4.1 Extended Kolmogorov hypotheses for Lagrangian timescales
      2. IV.4.2 Capturing space and time correlations through time changes
      3. IV.4.3 Accounting for the statistical signature of coherent structures
  5. V Similarities and differences with complex fluids
    1. V.1 Neither complete order nor disorder: thermal versus turbulence noise
    2. V.2 Fast-variable elimination and the Brownian limit
      1. V.2.1 Accounting for Brownian motion
      2. V.2.2 Brownian dynamics in turbulent flows
      3. V.2.3 Non-local and local closures
    3. V.3 Summary and two-way influences
  6. VI Statistical modeling of particle collisions in turbulent flows
    1. VI.1 Developments in sample and physical spaces
    2. VI.2 Limitations of present collision kernels
    3. VI.3 The molecular chaos assumption
    4. VI.4 Jump-diffusion processes as modeling tools
  7. VII Conclusions and perspectives

##  I The physics of dispersed turbulent two-phase flows

What is a dispersed two-phase flow and what are its main characteristics? To
provide answers to these questions, it is best to let Nature do the talking
through three examples among a wide range of applications. It is actually
difficult to do justice to the richness and variety of applications involving
dispersed two-phase flows and we have retained situations corresponding mainly
to natural phenomena, leaving out a number of interesting industrial processes
in which these flows play a key role.

These examples are helpful to introduce the fundamental interactions (Sec.
I.2), the governing equations (Sec. I.3), and the statistical issue which is
the main theme of this review in Sec. I.4. The probabilistic framework is
outlined in Sec. II and from then on we follow the modeling road. The issues
related to particle transport are analyzed in Sec. III and we present current
models as well as recent developments for the velocity of the fluid seen in
Sec. IV. Connections with soft matter are discussed in Sec. V, before
addressing particle collisions in Sec. VI and proposing some conclusions on
the road ahead in Sec. VII.

###  I.1 A rich tapestry of applications

####  I.1.1 Colloid suspensions and river deltas

The first illustration concerns colloidal suspension stability (detailed
accounts can easily be found in textbooks on colloids [29, 48, 49], among a
vast literature on the subject). As sketched in Fig. 1(a), colloids are small
particles sensitive to Brownian effects and which do not sediment (cf. Sec.
I.3.2 for a more precise definition). In a liquid solution, these colloids are
also subject to forces acting between them. These interactions are well
captured by the DLVO theory (see extensive descriptions in [49, 53] on the
formulation of the DLVO theory), which combines attractive forces (typically
due to van der Waals forces) when colloids are very close one to another and
repulsive ones (typically due to the repulsion between overlapping double
layers, which result from an excess of counter-ions in the vicinity of a
charged surface). The range of these repulsive forces corresponds to the Debye
length, which is a function of the local ion concentration (through the ionic
strength of the solvent). When the chemical conditions are such that the Debye
length is large, colloids repel one another and the suspension is stable. When
the chemical conditions are such that the Debye length is small, colloids can
get close enough for the attractive forces to take over. In that case,
aggregates begin to form until their size is too large for their weight to be
balanced by Brownian effects and the suspension is unstable.

(a) Stability of colloidal suspensions.

(b) Dynamics of a river delta.

Figure 1: The role of dispersed two-phase flows in the formation of river
deltas: (a) The formation of aggregates leads to unstable colloidal
suspensions (source: Wikipedia); (b) The resulting deposition forms delta
beds, as for the Copper River in the Gulf of Alaska (source: NASA Earth
Observatory)).

This is the much-simplified but general picture partly explaining the
formation of river deltas, as illustrated in Fig. 1(b) [1]. When colloids are
embedded in rivers with a small salt content (thus, a large Debye length), the
suspension is stable and colloids are basically carried along by the river
flow. When they reach more salty water from the ocean (thus, inducing a small
Debye length), colloids start to agglomerate and deposit on the river beds
forming river deltas over time. The rate at which this aggregation and
sedimentation take place is governed by these chemical conditions but also by
hydrodynamical ones, such as currents, recirculations, local turbulence
effects or eddies created by the river delta geography. It is also interesting
to note that this physical situation implies direct links between microscopic
scales (what happens to a colloidal suspension) and macroscopic ones (the
formation of river delta and river shores).

####  I.1.2 Plume dispersion and bubbly flows

For all its familiarity, plume dispersion is one of the best example to
exhibit the effects of fluid turbulence. As seen in Fig. 2(a), a stream of gas
material (which is often steam and not necessarily a noxious substance) is
released out of a chimney and transported by the wind [100]. The plume is
rapidly spreading due to the local turbulent fluctuations of the wind flow and
the initially-compact form of the released cloud is progressively torn apart
by the wind random motions so that it blends in the surrounding flow. Far less
familiar, however, is another example of plume dispersion, namely hydrothermal
vents through fissures on seabeds as shown in Fig. 2(b). Discovered in the
late 20th century in connection with tectonic plate theory, they are one of
the hallmarks of the volcanic activity near mid-ocean ridges and thus related
to the formation of oceanic crust (oceanic lithosphere) as well as to under-
water life [69]. They involve superheated water, heated by hot rocks in the
mantle, rushing through the sea floor. When this hot water mixes with cold sea
water, minerals precipitate giving them various colors depending on the type
of minerals (the name of black smokers illustrated here is when the plume is
black).

(a) Smoke plume dispersing from an industrial chimney.

(b) Submarine hydrothermal vent in the Rainbow site.

(c) Rising bubbles in a glass.

Figure 2: Illustrations showing the dispersion of two-phase flows in marine
and atmospheric flows: (a) A plume of smoke from the Dunbar Cement Works
chimney dispersing in the air (source: Walter Baxter, geograph.org.uk); (b)
Hydrothermal vent releasing a black smoke due to the precipitation of
dissolved minerals in cold seawater (source: Ifremer); (c) Bubble of Champagne
rising in a glass (Reprinted from [66] with permission from American Chemical
Society).

At first sight, it is not obvious that we are dealing with a dispersed two-
phase flow. However, it is not only possible but appropriate in such
situations to describe the flow of released materials as a set of fluid-like
individual elements. This corresponds to the important tracer-particle limit
obtained by considering particles having vanishing inertia, but still buoyancy
forces if there is a temperature difference with the surrounding fluid as in
the case of hydrothermal vents. Another similar point source dispersion
problem is illustrated in Fig. 2(c) where a stream of bubbles is rising in a
liquid [66], clearly showing this time the dispersed two-phase flow nature of
the flow. Note that, due to the adsorption of dissolved carbon dioxide, these
bubbles grow as they rise from the bottom of the glass, which impacts their
dynamics and rising velocities, but remain as individual bubbles up to the
free surface.

####  I.1.3 The physics of sand dunes

With the publication in 1941 of his book on desert dunes, reprinted in [3], R.
A. Bagnold took a much-to-be-admired initiative and can be considered as the
founding father of a new topic in physics, namely ‘the physics of blown sand’.
This is a fascinating subject not only for its aesthetic appeal (cf. Fig.
3(a)) but also for the diversity of the physical processes involved [61].
Under the actions of strong-enough winds, sand particles are blown off the
edge of sand dunes, as manifested by the layers streaming away in Fig. 3(a).
This corresponds to the process of particle resuspension by turbulent flows,
which is of wide applicability but still contains puzzling issues [44].
Contrary to colloids, sand particles have non-negligible inertia. Therefore,
if the wind velocity diminishes or if they were entrained by sudden gusts in
an intermittent wind regime, they are transported over some distance but can
deposit again. They are then often referred to as saltating particles. When
they come in contact with another dune, the kinetic energy imparted by the
wind flows means that they hit the surface of previously-deposited layers of
sand particles with some force and this can create ‘splashing effects’,
inducing a complex interplay between the two mirror processes of resuspension
and deposition. This leads also to an ever-changing landscape of sand dunes
whose formation is tightly coupled to hydrodynamical conditions. Indeed, when
a sand dune is formed or when its slope is modified, there is a back-effect on
the strength as well as the orientation of the wind velocity which, in turn,
changes the way sand particles are being blown off, and so on.

(a) Sand being blown off dunes by the wind.

(b) Sahara dust over the Atlantic Ocean.

Figure 3: Illustrations showing how tiny solid particles can be transported
across large scales: (a) Picture of sand being blown off dunes (left, source:
Flickr); (b) Satellite observation of Sahara dust transported in the
atmosphere above the Atlantic Ocean (right, source: NASA Earth Observatory).

Sand particles have different morphologies and sizes which can cover a range
of values. We are therefore dealing with poly-disperse turbulent two-phase
flows and the consequence of the difference in sand particle inertia is that
their responses to wind solicitations vary greatly. Above the surface of
desert dunes, sand particles interact and collective effects must be
considered. On the other hand, once airborne, the tiniest sand particles can
be carried over very large distances, in which case the main driving force is
the wind currents that disperse them over regions at the scale of a continent
or an ocean, as shown in Fig. 3.

###  I.2 An intricate range of phenomena

As it transpires from these examples, dispersed two-phase flows involve two
non-miscible thermodynamical phases with one phase, the dispersed one, present
as a set of discrete elements (bubbles, droplets, solid particles) embedded in
the other phase, the continuous one (also referred to as the carrier phase),
which can be a gas or a liquid flow. Dispersed two-phase flows represent an
important subclass within the general class of two-phase flows. Note that
particle-laden flows are always in that category whereas the topology of the
interface between two fluid phases in the dispersed flow regime implies that,
for the same volume, the surface of contact is increased compared to the
separated one. When bubbles or droplets are present, internal processes (such
as recirculation, surface tension) can take place but we consider that, once
formed, bubbles and droplets are treated as single entities and we are
essentially concerned with their dynamics in turbulent flows. In that sense,
apart from their different densities or similar properties, we do not
distinguish between solid particles, droplets or bubbles, and refer to them as
‘discrete particles’ or simply ‘particles’.

####  I.2.1 The two fundamental interactions: particle-fluid and particle-
particle

These examples are also useful to bring out that there are two fundamental
interactions involved: the fluid-particle and the particle-particle
interactions. The fluid-particle interaction is at play in particle transport
as can be seen from colloids carried by rivers, fluid-like elements dispersing
in the atmosphere, etc. In the vast majority of cases where flows are
turbulent, we are therefore talking about turbulent convective transport. Yet,
even when a fluid is at rest, Brownian motion can induce diffusive transport
and we can describe both effects as fluid-induced transport. The particle-
particle interaction is at play in the formation (or breakup) of colloid
aggregates or when saltating particles hit layers of deposited sand particles.
By particle-particle interaction, we are referring to short-range
interactions, such as the forces occurring between colloids whose surfaces are
separated by distances of the order of the Debye length, and not just to
particles bouncing off one another. Nevertheless, compared to hydrodynamical
scales, the Debye length is negligible and these particle-particle
interactions can be regarded as near-contact ones.

It follows that we are at the crossroads between several scientific domains
since these fundamental interactions involve aspects pertaining to fluid
mechanics (in particular, turbulence) but also to contact mechanics and to
interface chemistry. A further remark is to note that, for dispersed two-phase
flows, we are essentially dealing with a two-step mechanism. The first one is
the transport step involved in particle dispersion but also responsible for
bringing discrete particles nearby. The second step, which is the collision
step (whether it is actual collisions, agglomeration, etc.), can only take
place when particles are in the immediate vicinity of one another, as the
result of the transport step. This successive mechanism is the reason why
issues related to the transport step are given more emphasis in this review
and also to pave the way for specific analysis of the open questions related
to particle-particle interactions in turbulent flows.

####  I.2.2 Addressing complexity

Through these issues, we are touching upon the general question of ‘how to
address complexity?’ This is far too vast a subject to be discussed here in
details. Nevertheless, this question lies in the background and is connected
to some points discussed in later developments. For instance, it will be seen
that the selection of the variables retained to characterize particles as
mechanical systems is related to our choice to describe as ‘noise’ or
‘disorder’ some external effects. One of the main themes of this review is to
emphasize that this is a key point if we are to devise complete and well-based
descriptions so that these aspects should be carefully weighted. To put it in
more philosophical terms, we should not be too quick in qualifying as noise
some effects, since what we may call ‘chaos’ is perhaps an order whose reading
we have not yet learned. A second example concerns the modeling approach to
follow when trying to model one variable which is displaying ‘complex
behavior’. Should we retain the idea of devising a single model, which is then
likely to become complex? Or should we consider that complexity is perhaps
best captured by considering random alternation between several models, each
of which being simple? This debate reappears in Sec. IV.4.3.

###  I.3 The governing equations

A dispersed two-phase flow is a composite system in which the continuous phase
is described in terms of fields whereas the dispersed phase is described in
terms of particles. In the following, we provide the basic governing equations
for the two phases, with more emphasis on the particle one since their
statistical treatment is the main theme of the present review, as indicated in
Sec. I.4.

####  I.3.1 The Navier-Stokes equations for fluid flows

At the continuum level of description, the governing equations of the fluid
phase are the continuity, the Navier-Stokes (NS) and the transport equations
for a set of $N_{s}+1$ scalar fields which gathers the $N_{s}$ relevant
species mass fractions
$\bm{\phi}_{f}(t,\mathbf{x})=(\phi_{f,l})_{l=1,\ldots,N_{s}}$ (if chemical
reactions are involved) to which the fluid enthalphy (or another
thermodynamical potential) is added, along with an equation of state for
compressible flows. In the case of constant-property flows, these equations
are [89]

$\displaystyle\frac{\partial U_{{\rm f},k}}{\partial x_{k}}=0~{},$ (1a)
$\displaystyle\frac{\partial U_{{\rm f},i}}{\partial t}+U_{{\rm
f},k}\,\frac{\partial U_{{\rm f},i}}{\partial x_{k}}=-\frac{1}{\rho_{\rm
f}}\frac{\partial P_{\rm f}}{\partial x_{i}}+\nu_{\rm
f}\,\frac{\partial^{2}U_{{\rm f},i}}{\partial x_{k}\partial x_{k}}~{},$ (1b)
$\displaystyle\frac{\partial\phi_{{\rm f},l}}{\partial t}+U_{{\rm
f},k}\,\frac{\partial\phi_{{\rm f},l}}{\partial x_{k}}=\Gamma_{\rm
f}\,\frac{\partial^{2}\phi_{{\rm f},l}}{\partial x_{k}\partial x_{k}}+S_{{\rm
f},l}~{},$ (1c)

where $\mathbf{U}_{\rm f}(t,\mathbf{x})$ and $P_{\rm f}(t,\mathbf{x})$ are the
fluid velocity and pressure fields respectively, $\rho_{\rm f}$ the fluid
density, $\nu_{\rm f}$ its kinematic viscosity and $\Gamma_{\rm f}$ the scalar
diffusivity. In Eq. (1c), the last term on the right-hand side (rhs) is a
reactive source term $S_{{\rm f},l}=\hat{S}_{{\rm f},l}(\bm{\phi}_{\rm
f}(t,\mathbf{x}))$ where $\hat{S}_{{\rm f},l}$ are known functions. Eqs. (1)
are the traditional field equations for fluid mechanics. They are generally
obtained directly at the hydrodynamical level of description by writing
balance equations for fluid mass, momentum and energy densities in which
constitutive relations are introduced using the notion of local thermodynamic
equilibrium [12] (a much-recommended new presentation is given in [108]). They
can also be derived from the kinetic theory (see [86, chapter 7] on the
subject, among several references), which is more in line with the theme of
this review (cf. discussions in Sec. V (in particular, in Sec. V.2.3) as well
as in Sec. VI).

In the present work, we are not interested in reactive flows and we do not
need to specify the expressions of reactive source terms. The important point
is, however, that these terms are obtained from the instantaneous scalar
variables $\bm{\phi}_{\rm f}$ through a (usually) non-linear but known
relation. This is relevant when modeling reactive flows with similar source
terms and the PDF approach to single-phase reactive flows has the key
advantage of treating them without approximation (see extensive discussions in
[90, 41]). Note that the Navier-Stokes equation in Eq. (1b) is written at
locations $\mathbf{x}$ where no particles are present at time $t$. The
equation for scalars, Eq. (1c), is given for the sake of completeness and is
only referred to (without reactive source terms) when discussing the diffusive
limit for passive scalars in Sec. V.2. From now on, we limit ourselves to
dynamical aspects of incompressible constant-property flows represented by
Eqs. (1a)-(1b).

####  I.3.2 The equations of motion for particles

To describe the particle phase, it is natural to adopt a Lagrangian point of
view. This approach consists in tracking explicitly each particle by solving
the evolution equations for the state vector $\mathbf{Z}_{\rm
p}=(\mathbf{X}_{\rm p},\mathbf{U}_{\rm p},\bm{\Omega}_{\rm p})$ where
$\mathbf{X}_{\rm p}$ is the location of the particle center of mass while
$\mathbf{U}_{\rm p}$ and $\bm{\Omega}_{\rm p}$ are the translational and
rotational velocities, respectively (more details on how to select the
relevant variables entering the state vector are provided in Sec. II).
Applying the fundamental laws of classical mechanics, a general system of
evolution equations for the particle state vector are the following ordinary
differential equations (ODE):

$\displaystyle\frac{\mathrm{d}\mathbf{X}_{\rm p}}{\mathrm{d}t}$
$\displaystyle=\mathbf{U}_{\rm p}~{},$ (2a) $\displaystyle m_{\rm
p}\frac{\mathrm{d}\mathbf{U}_{\rm p}}{\mathrm{d}t}$
$\displaystyle=\mathbf{F}_{\rm f\to p}+\mathbf{F}_{\rm p\to p}+\mathbf{F}_{\rm
ext}~{},$ (2b) $\displaystyle I_{\rm p}\frac{\mathrm{d}\bm{\Omega}_{\rm
p}}{\mathrm{d}t}$ $\displaystyle=\mathbf{M}_{\rm f\to p}+\mathbf{M}_{\rm p\to
p}~{},$ (2c)

where $m_{\rm p}$ the mass of the particle and $I_{\rm p}$ its moment of
inertia. In Eqs. (2b)-(2c), $\mathbf{F}_{\rm f\to p}$ and $\mathbf{M}_{\rm
f\to p}$ represent the forces and torques exerted by the fluid on each
particle, $\mathbf{F}_{\rm p\to p}$ and $\mathbf{M}_{\rm p\to p}$ the forces
and torques due to particle-particle interactions and $\mathbf{F}_{\rm ext}$
forces due to external fields (such as gravity).

The derivation of approximate expressions for hydrodynamic forces acting on a
particle embedded in a fluid flow has a long history, going back to the 19th
century with the Basset-Boussinesq-Oseen (BBO) expression. However, in spite
of numerous studies, it cannot be said that this issue has been completely
solved, though a state-of-the-art formulation has emerged. Present
formulations retain the drag, lift, added-mass and Basset forces [21, 36, 70,
62, 15] which are added to the gravity force, but the expressions of these
forces vary according to different authors. The lift force is a special case
since several expressions are proposed but each one obtained in a very
specific or asymptotic situation. The result is a catalog of widely different
formulations [21, 43] so that there is, at the moment, no well-established
consensus about a unique expression (the ‘optimal lift force’ in McLaughlin’s
works [71, 113] is perhaps the most general expression available). Approximate
expressions for Basset force are also still debated and, consequently, present
models for these two forces cannot be regarded as being reliable enough. In
the present work, we are more interested in a general and well-accepted form
to develop probabilistic approaches and, for these reasons, the lift and
Basset forces are left out. When lift forces are not considered, the particle
rotation $\bm{\Omega}_{p}$ does not play an explicit role in the dynamics of
spherical particles and can also be left out of the particle state vector
which is then reduced to $\mathbf{Z}_{\rm p}=(\mathbf{X}_{\rm
p},\mathbf{U}_{\rm p})$. Similarly, we do not consider electro- or thermo-
phoresis phenomena which, if present, are accounted for by the addition of
corresponding formulations in the particle momentum equation (see discussion
in [72, section 2.3]). Nevertheless, it is important to realize that the
framework described below is an open one and can easily be extended to account
for these effects. Hence, present restrictions should not be seen as
limitations but rather as an invitation to widen the range of applications by
introducing future developments related to these additional physical
phenomena.

For small spherical particles with diameters $d_{\rm p}\sim\eta_{\rm K}$,
where $\eta_{\rm K}$ is the Kolmogorov length scale that represents the
smallest length scale for fluid motions in a turbulent flow (to be defined in
Sec. III.3), the reference formulation of the particle momentum equation is
derived in [36] and writes

$m_{\rm p}\frac{\mathrm{d}\mathbf{U}_{\rm p}}{\mathrm{d}t}=\frac{\pi d_{\rm
p}^{3}}{6}\rho_{\rm f}\frac{\mathrm{D}\mathbf{U}^{\mathcal{V}_{\rm p}}_{\rm
s}}{\mathrm{D}t}+\frac{\pi d_{\rm p}^{3}}{6}(\rho_{\rm p}-\rho_{\rm
f})\mathbf{g}\\\ +\frac{1}{2}\frac{\pi d_{\rm p}^{2}}{4}\rho_{\rm f}C_{\rm
D}|\,\mathbf{U}^{\mathcal{S}_{\rm p}}_{\rm s}-\mathbf{U}_{\rm
p}\,|(\mathbf{U}^{\mathcal{S}_{\rm p}}_{\rm s}-\mathbf{U}_{\rm p})+\frac{\pi
d_{p}^{3}}{6}\rho_{\rm f}C_{\rm
A}\left(\frac{\mathrm{D}\mathbf{U}^{\mathcal{V}_{\rm p}}_{\rm
s}}{\mathrm{D}t}-\frac{\mathrm{d}\mathbf{U}_{\rm p}}{\mathrm{d}t}\right),$ (3)

where $\mathrm{D}/\mathrm{D}t$ refers to the derivative along a fluid particle
trajectory

$\frac{\mathrm{D}missing}{\mathrm{D}t}=\frac{\partial}{\partial t}+U_{{\rm
f},k}\frac{\partial}{\partial x_{k}}.$ (4)

In Eq. (3), the first two terms on the rhs are the fluid acceleration and the
buoyancy force (involving the difference between the particle density
$\rho_{\rm p}$ and the fluid one), the third term is the general form of the
drag force written with the drag coefficient $C_{\rm D}$ while the fourth term
is the added-mass force written with a coefficient $C_{\rm A}$ (usually taken
as $C_{\rm A}=1/2$). The added-mass force depends on the difference between
the acceleration of the fluid seen (along its own trajectory) and the particle
acceleration while the drag force is expressed in terms of the velocity slip
$\mathbf{U}_{\rm s}^{\mathcal{S}_{\rm p}}-\mathbf{U}_{\rm p}$. The last two
terms involve $\mathbf{U}^{\mathcal{S}_{\rm p}}_{\rm s}$ and
$\mathbf{U}^{\mathcal{V}_{\rm p}}_{\rm s}$, which are the fluid velocities
averaged over the surface $\mathcal{S}_{\rm p}$ and the volume
$\mathcal{V}_{\rm p}$ of the particle, respectively, i.e.

$\displaystyle\mathbf{U}^{\mathcal{S}_{\rm p}}_{\rm s}$
$\displaystyle=\frac{1}{\mathcal{S}_{\rm p}}\int_{\mathcal{S}_{\rm
p}}\mathbf{U}_{\rm f}(t,\mathbf{r})\,d\mathbf{r}~{},$ (5)
$\displaystyle\mathbf{U}^{\mathcal{V}_{\rm p}}_{\rm s}$
$\displaystyle=\frac{1}{\mathcal{V}_{\rm p}}\int_{\mathcal{V}_{\rm
p}}\mathbf{U}_{\rm f}(t,\mathbf{r})\,d\mathbf{r}~{}.$ (6)

These velocities are expressed by a series expansion around the ‘velocity’ at
the particle center, which introduces the notion of the velocity of the fluid
seen, $\widehat{\mathbf{U}}_{\rm s}(t)=\mathbf{U}^{\rm ud}_{\rm
f}(t,\mathbf{X}_{\rm p}(t))$

$\displaystyle\mathbf{U}^{\mathcal{S}_{\rm p}}_{\rm s}$
$\displaystyle\simeq\widehat{\mathbf{U}}_{\rm s}+\frac{d_{\rm
p}^{2}}{24}\left(\nabla^{2}\mathbf{U}^{\rm ud}_{\rm
f}\right)(t,\mathbf{X}_{\rm p}(t))~{},$ (7)
$\displaystyle\mathbf{U}^{\mathcal{V}_{\rm p}}_{\rm s}$
$\displaystyle\simeq\widehat{\mathbf{U}}_{\rm s}+\frac{d_{\rm
p}^{2}}{40}\left(\nabla^{2}\mathbf{U}^{\rm ud}_{\rm
f}\right)(t,\mathbf{X}_{\rm p}(t))~{}.$ (8)

The notation $\mathbf{U}^{\rm ud}_{\rm f}(t,\mathbf{X}_{\rm p}(t))$ designates
the ‘undisturbed velocity field’ that would exist if the particle were not
present at the location $\mathbf{x}=\mathbf{X}_{\rm p}(t)$ at time $t$ [36,
70] and the last terms added to $\widehat{\mathbf{U}}_{\rm s}(t)$ on the rhs
of Eqs. (7) and (8) are the Faxen terms [36, 70, 21].

For small particles (say $d_{p}\leq\eta_{\rm K}$), the Faxen terms are small
enough to be neglected which means that, for the expression of these
hydrodynamical forces, particles are considered as mere points. This
corresponds to the point-wise approximation for particles and a classical form
of the particle momentum equation is

$m_{\rm p}\frac{\mathrm{d}\mathbf{U}_{\rm p}}{\mathrm{d}t}=\frac{\pi d_{\rm
p}^{3}}{6}\rho_{\rm f}\frac{\mathrm{D}\mathbf{U}_{\rm
s}}{\mathrm{D}t}+\frac{\pi d_{\rm p}^{3}}{6}\left(\rho_{\rm p}-\rho_{\rm
f}\right)\mathbf{g}\\\ +\frac{1}{2}\frac{\pi d_{\rm p}^{2}}{4}\rho_{\rm
f}\,C_{\rm D}|\,\mathbf{U}_{\rm s}-\mathbf{U}_{\rm p}\,|(\mathbf{U}_{\rm
s}-\mathbf{U}_{\rm p})+\frac{\pi d_{\rm p}^{3}}{6}\,C_{\rm A}\,\rho_{\rm
f}\left(\frac{\mathrm{D}\mathbf{U}_{\rm
s}}{\mathrm{D}t}-\frac{\mathrm{d}\mathbf{U}_{\rm p}}{\mathrm{d}t}\right)~{}.$
(9)

The velocity of the fluid seen is now given as the local instantaneous value
of the fluid velocity at the same time and at the particle position, i.e.
$\mathbf{U}_{\rm s}(t)=\mathbf{U}_{\rm f}(t,\mathbf{X}_{\rm p}(t))$, where
$\mathbf{U}_{\rm f}(t,\mathbf{x})$ represents the fluid velocity field without
having to distinguish undisturbed and disturbed fluid velocity fields anymore.
For particles heavier than the fluid $\rho_{p}\gg\rho_{f}$, it can be shown
that the drag force is the dominant one [70, 36, 76] and the particle momentum
equation is then further reduced to

$\frac{\mathrm{d}\mathbf{U}_{\rm p}}{\mathrm{d}t}=\frac{\mathbf{U}_{\rm
s}-\mathbf{U}_{\rm p}}{\tau_{\rm p}}+\mathbf{g}$ (10)

where the drag force is written so as to bring out the particle relaxation
timescale $\tau_{\rm p}$

$\tau_{\rm p}=\frac{\rho_{\rm p}}{\rho_{\rm f}}\frac{4\,d_{\rm p}}{3\,C_{\rm
D}|\mathbf{U}_{\rm s}-\mathbf{U}_{\rm p}|}.$ (11)

This timescale, which is a measure of particle inertia, is the key notion for
particle transport. More precisely, $\tau_{\rm p}$ represents the timescale
needed for particle velocities to adjust to the local fluid velocity seen. In
the Stokes regime, which is valid when ${\rm Re}_{\rm p}\leq 1$, with ${\rm
Re}_{\rm p}=|\,\mathbf{U}_{\rm s}-\mathbf{U}_{\rm p}\,|\,d_{\rm p}/\nu_{\rm
f}$ the particle Reynolds number, the drag coefficient is $C_{\rm D}=24/{\rm
Re}_{\rm p}$. In that case, we retrieve the classical form $\mathbf{F}_{\rm
f\to p}^{\rm Drag}=3\pi\rho_{\rm f}\,\nu_{\rm f}\,d_{\rm
p}\left(\mathbf{U}_{\rm s}-\mathbf{U}_{\rm p}\right)$ and the particle
relaxation timescale is given by the well-known expression

$\tau_{\rm p}=\frac{\rho_{\rm p}}{\rho_{\rm f}}\frac{d_{\rm p}^{2}}{18\nu_{\rm
f}}.$ (12)

For general values of ${\rm Re}_{\rm p}$, the drag coefficient is obtained
through empirical correlations and an often-retained formula is [21, 16]

$C_{\rm D}=\begin{cases}\dfrac{24}{{\rm Re}_{\rm p}}\left[\,1+0.15\,{\rm
Re}_{\rm p}^{0.687}\,\right]&\text{if}\;{\rm Re}_{\rm p}\leq 1000,\\\
0.44&\text{if}\;{\rm Re}_{\rm p}\geq 1000.\end{cases}$ (13)

This correlation is for isolated particles or when particle concentration is
not high enough for collective effects, representing hydrodynamical influences
of each particle on its neighbors (the ‘wake effect’), to be significant. If
needed, particle wake effects can be accounted for with modified correlations
based, for example, on the local particle volumetric fraction [21].

In most cases involving particles or droplets in a turbulent flow, the
preceding expressions provide a satisfactory picture and are sufficient to
describe the influence of fluid flows on particle dynamics. For diameters
larger than the Kolmogorov length scale, it is seen from the general form of
the particle momentum equation, cf. Eq. (3), and the expressions given in Eqs.
(5)-(6) that a non-negligible particle size induces a filtering effect and
that fluid velocity fluctuations with length scales smaller than $d_{\rm p}$
tend to be smoothed out and act as an underlying noise. Specific developments
can then be considered (see for example interesting proposals in [31, 30]
among other ideas), but since we concentrate essentially on turbulent
dispersion and collisional effects, the particle point-wise approximation is
retained in the present review from now on (see further comments in Sec. I.4).

Even when added-mass forces are neglected, the acceleration of the fluid seen
(i.e., $\mathrm{D}\mathbf{U}_{\rm s}/\mathrm{D}t$) is sometimes retained in
the particle velocity equation which is then written as

$\frac{\mathrm{d}\mathbf{U}_{\rm p}}{\mathrm{d}t}=-\frac{1}{\rho_{\rm
p}}\nabla P_{\rm f}+\frac{\mathbf{U}_{\rm s}-\mathbf{U}_{\rm p}}{\tau_{\rm
p}}+\mathbf{g}~{},$ (14)

where we have used the Euler form of the NS equations by neglecting viscosity
to relate the fluid pressure gradient and acceleration as

$\frac{\mathrm{D}\mathbf{U}_{\rm s}}{\mathrm{D}t}=-\frac{1}{\rho_{\rm
f}}\nabla P_{\rm f}+\mathbf{g}~{}.$ (15)

Eq. (14), or variations of it when $C_{\rm A}\neq 0$, is useful for small
sediments or bubbles, as discussed below.

##### Brownian effects and particle collisions

So far, we have only considered forces arising from the continuum description
of fluid flows. However, small enough particles are sensitive to Brownian
effects due to the molecular nature of the fluid. This vague notion of ‘small
enough’ is quantified by introducing the criterion that, in the absence of a
fluid flow, the particle settling velocity is counterbalanced by the random
velocities imparted by collisions with the molecules of the fluid, so that
these particles do not sediment, as the colloids discussed in Sec. I.1.1.
Using a one-dimensional formulation in the direction aligned with gravity and
the equipartition theorem of statistical physics, the diameter $d_{\rm p}$ of
these particles can be estimated by

$\sqrt{\frac{k_{\rm B}\,\Theta_{\rm f}}{m_{\rm p}}}\sim\tau_{\rm
p}\,g\Longrightarrow d_{\rm p}\simeq\left(\frac{\rho_{\rm f}^{2}\nu_{\rm
f}^{2}}{\rho_{\rm p}^{3}\,g^{2}}\,k_{\rm B}\,\Theta_{\rm f}\right)^{1/7}~{},$
(16)

where $k_{\rm B}$ is the Boltzmann constant and $\Theta_{f}$ the fluid
temperature. This defines more precisely colloidal particles which have
therefore a diameter of the order of a few microns ($d_{\rm p}\leq 1-2\mu m$)
while, for $d_{\rm p}\geq 5-10\mu m$, particles are called inertial and
Brownian effects become negligible. In practice, it is best to include
Brownian motions in the particle momentum equation whatever the particle
diameter so that its effects diminish continuously when considering increasing
particle diameters without having to introduce an artificial cut-off between
colloidal and inertial particles. Although implementing Brownian effects is
straightforward, it is connected to interesting questions and is addressed in
more details in Sec. V.2. We also need to account for particle-particle
interactions, typically collisions, in the particle momentum equation. This
issue, which is more involved than Brownian effects, is given specific
attention in Sec. VI. Moreover, Brownian motions are typically expressed by
Wiener processes while the effects of random collisions are accounted for with
Poisson processes. Since these two fundamental stochastic processes are only
introduced in Sec. II, it is best to address the issues related to Brownian
effects and particle collisions in Secs. V.2 and VI, respectively, after
having detailed stochastic models related to turbulent dispersion in Sec. IV.

##### Remarks on bubbly flows

Compared to solid particles and droplets in a gas, the case of small bubbles
in a liquid raises specific concerns and a few comments are in order. First,
bubbles tend to have larger sizes (due to surface tension and the resulting
inner pressure) which implies that the volume fraction occupied by the bubble
phase becomes quickly appreciable. In practice, this volumetric back effect
from bubbles to the liquid phase is important and the liquid-phase equations
need to be modified accordingly (through some interface tracking methods).
Second, even if the point-wise approximation is less valid for bubbles, we can
still follow their center-of-mass motion. Yet, since bubble density is much
smaller than the liquid one ($\rho_{\rm p}\ll\rho_{\rm f}$), we cannot neglect
the added-mass force anymore. This leads to the bubble momentum equation

$\frac{\mathrm{d}\mathbf{U}_{\rm p}}{\mathrm{d}t}=-\frac{1}{\rho_{\rm
p}}\left(\frac{1+C_{\rm A}}{1+C_{\rm A}\rho_{\rm f}/\rho_{\rm p}}\right)\nabla
P_{\rm f}+\frac{\mathbf{U}_{\rm s}-\mathbf{U}_{\rm p}}{\tau_{\rm
p}}+\mathbf{g}~{}.$ (17)

Assuming that we are in the Stokes regime, the relaxation timescale in Eq.
(17) is transformed to become

$\tau_{\rm p}=\left(\frac{\rho_{\rm p}+C_{\rm A}\rho_{\rm f}}{\rho_{\rm
f}}\right)\frac{d_{\rm p}^{2}}{18\nu_{\rm f}}.$ (18)

As a consequence, we must know not only $\mathbf{U}_{\rm s}$ but also
$\mathrm{D}\mathbf{U}_{\rm s}/\mathrm{D}t$. To slightly anticipate on
developments to come, it will be seen, however, that models based on a two-
step formulation (cf. Sec. IV.4.2) generate the increments of $\mathbf{U}_{\rm
s}$ (which corresponds to $\mathrm{d}\mathbf{U}_{\rm s}/\mathrm{d}t$, the
acceleration along bubble trajectories) as well as the velocity increments
along a fluid particle trajectory in the so-called Lagrangian step (which
corresponds to $\mathrm{D}\mathbf{U}_{\rm s}/\mathrm{D}t$). This means that
new models for $\mathbf{U}_{\rm s}$ are able to provide the information needed
to account for the added-mass force and we can focus on what is at stake from
a statistical point of view.

###  I.4 The issue to address in statistical physics

With these clarifications, we consider discrete particle dynamics expressed by
the system

$\displaystyle\frac{\mathrm{d}\mathbf{X}_{\rm p}}{\mathrm{d}t}$
$\displaystyle=\mathbf{U}_{\rm p}~{},$ (19a)
$\displaystyle\frac{\mathrm{d}\mathbf{U}_{\rm p}}{\mathrm{d}t}$
$\displaystyle=\frac{\mathbf{U}_{\rm s}-\mathbf{U}_{\rm p}}{\tau_{\rm
p}}+\mathbf{g}~{},$ (19b)

where $\tau_{\rm p}$ is given by Eqs. (11) and (13). As explained for bubbles,
we have in fact only slightly reduced the applicability of the dynamical
formulation. However, by doing so, we have stripped down the problem to its
bare statistical essential which is the treatment of the velocity of the fluid
seen. Indeed, the advantage of Eq. (19b) is to bring out that the real issue
is how to model the ‘driving force’, $\mathbf{U}_{\rm s}(t)$, more than the
various forms taken by the particle momentum equation depending on which
forces are retained. This can be rephrased using an analogy where the particle
momentum equation is regarded as a (non-linear) filter from which we derive
statistics of interest which are outputs of the filter. However, the real
issue is to model the input (or the driving force) of this filter, which is
the velocity of the fluid seen $\mathbf{U}_{\rm s}(t)$. In the rare situations
where the fluid flow is fully resolved and where we have access to the
instantaneous velocity field $\mathbf{U}_{\rm f}(t,\mathbf{x})$, then
$\mathbf{U}_{\rm s}(t)$ is exactly known and we are dealing with a fully-
characterized deterministic time signal as the input of the particle-momentum
filter, or a series of such time functions if we are considering a number of
discrete particles. In that case, we can proceed with the particle kinetic
variables $(\mathbf{X}_{\rm p},\mathbf{U}_{\rm p})$ as in classical mechanics.
However, for the vast majority of turbulent flows the number of degrees of
freedom is so high (cf. estimations in Sec. III.3) that we have only access to
a limited information, usually in the form of the first few moments of the
fluid velocity field. The loss of information implied by such reduced, or
coarse-grained, descriptions is then reflected by the representation of
$\mathbf{U}_{\rm s}(t)$ as a stochastic process, which means that we enter the
world of probabilistic formulations.

##  II Reduced statistical descriptions and the probabilistic framework

There are basically two steps in dispersed two-phase flow modeling. The first
one deals with the selection of the variables gathered in the particle state
vector and with the construction of the stochastic processes used to model
particle dynamics. This step encompasses most of the issues related to physics
and answers the questions: How do we describe a mechanical system? How do we
represent its dynamical evolution? The second step concerns the probabilistic
framework guiding us from particle stochastic models to statistics of interest
in practical situations. This step is more mathematical in nature and answers
the questions: What are the main stochastic processes? How do we handle them
to stay clear of mathematical pitfalls?

To concentrate on the physical issues in the following sections while relying
on a safe probabilistic framework, we start with a mere outline of the
mathematical background. This is precisely what is not suggested to newcomers
in this subject who are encouraged to study in-depth the mathematical aspects
of stochastic processes to ‘learn their trade’. For that purpose, a highly-
recommended textbook is [85] which does a very good job at presenting
stochastic processes and, in particular, the interest of formulations in terms
of trajectories instead of considering only PDF equations. It is telling that
half of this book on polymer models is dedicated to the mathematical
properties of stochastic processes.

###  II.1 From microscopic to macroscopic descriptions

####  II.1.1 Mass density functions and mean fields

The PDF approach to single-phase and dispersed two-phase turbulent flows has
reached a mature level and detailed accounts can be found in [90, 89, 76, 41,
72, 73]. The formalism bears similarities with the classical formulation in
terms of distribution functions but there are also differences. First, the
normalization is in terms of the mass of the fluid or particle phases instead
of the number of particles. Second, this PDF framework grew out of the
probabilistic description of turbulent single-phase flows where the
‘microscopic level of description’ (i.e., the Navier-Stokes equations) is made
up by fields and, for dispersed two-phase flows, a combination of particles
and fields. This is in contrast with classical statistical physics where the
microscopic level involves always discrete particles (typically
atoms/molecules) whereas the macroscopic one corresponds to fields (the
hydrodynamical level). Consequently, a specific distinction is made between
Lagrangian and Eulerian PDFs to avoid confusions between descriptions where
the position is a variable (Lagrangian) or a parameter (Eulerian).

The methodology can be developed for general $N$-particle PDF descriptions and
state vectors (whose selections are addressed later on), but to keep simple
notations we consider essentially the one-particle PDF approach. This allows
to handle a large number of particles as samples instead of a large number of
pairs of particles for two-particle PDF models, etc. Adopting therefore a
Lagrangian standpoint and focusing on the particle phase, we introduce the
decomposition of the particle state vector $\mathbf{Z}=(\mathbf{X}_{\rm
p},\mathbf{Z}_{\rm c})$ where the particle location is always present and
where $\mathbf{Z}_{\rm c}$ is the complementary part of the state vector. In
sample space, the corresponding variables are noted
$\mathbf{z}=(\mathbf{y}_{\rm p},\mathbf{z}_{\rm c})$. Starting from the
unconditional Lagrangian PDF $p^{\rm L}(t;\mathbf{y}_{\rm p},\mathbf{z}_{\rm
c})$ and using $M_{\rm p}$ for the total mass of the discrete particles in the
domain, the Lagrangian and Eulerian mass density functions (MDFs) for the
disperse phase are defined by

$\displaystyle F_{\rm p}^{\rm L}(t;\mathbf{y}_{\rm p},\mathbf{z}_{\rm c})$
$\displaystyle=M_{\rm p}\,p^{\rm L}(t;\mathbf{y}_{\rm p},\mathbf{z}_{\rm c}),$
(20a) $\displaystyle F_{\rm p}^{E}(t,\mathbf{x};\mathbf{z}_{\rm c})$
$\displaystyle=F_{\rm p}^{L}(t;\mathbf{y}_{\rm p}=\mathbf{x},\mathbf{z}_{\rm
c})=\int F_{\rm p}^{\rm L}(t;\mathbf{y}_{\rm p},\mathbf{z}_{\rm
c})\delta(\mathbf{y}_{\rm p}-\mathbf{x})\,\mathrm{d}\mathbf{y}_{\rm p}~{}.$
(20b)

The key relation between Eulerian and Lagrangian descriptions is [5, 76]

$F_{\rm p}^{E}(t,\mathbf{x};\mathbf{z}_{\rm c})=\int p^{\rm
L}(t;\mathbf{x},\mathbf{z}_{\rm c}\,|\,t_{0};\mathbf{x}_{0},\mathbf{z}_{{\rm
c},0})\,F_{\rm p}^{E}(t_{0},\mathbf{x}_{0};\mathbf{z}_{{\rm
c},0})\mathrm{d}\mathbf{x}_{0}\,\mathrm{d}\mathbf{z}_{{\rm c},0}~{},$ (21)

which shows that the Lagrangian transition PDF is the propagator of the
information, here the Eulerian MDF from the initial state
$(t_{0};\mathbf{x}_{0},\mathbf{z}_{{\rm c},0})$ to the present one
$(t;\mathbf{x},\mathbf{z}_{\rm c})$. This central property explains why most
of the physics is contained in modeling the Lagrangian transition PDF. Once
the Eulerian MDF is defined, particle mean field properties are extracted. For
a particle variable written as $H_{\rm p}(t;\mathbf{Z}_{\rm c})$, its average
$\langle H_{\rm p}\rangle$ (which is a field variable), is defined as

$\alpha_{\rm p}(t,\mathbf{x})\,\rho_{\rm p}\langle H_{\rm
p}\rangle(t,\mathbf{x})=\int H_{\rm p}(t;\mathbf{z}_{\rm c})F_{\rm p}^{\rm
E}(t,\mathbf{x};\mathbf{z}_{\rm c})\,\mathrm{d}\mathbf{z}_{\rm c}~{},$ (22)

where $\alpha_{\rm p}(t,\mathbf{x})$ is the mean particle volumetric fraction.
In a discrete sense, with $N_{\rm p}$ stochastic particles, the definitions of
Lagrangian and Eulerian MDFs become

$\displaystyle F_{{\rm p},N_{\rm p}}^{\rm L}(t\,;\,\mathbf{y}_{\rm
p},\mathbf{z}_{\rm c})$ $\displaystyle=\sum_{i=1}^{N_{\rm p}}m_{\rm
p}^{[i]}\delta(\mathbf{y}_{\rm p}-\mathbf{x}_{\rm
p}^{[i]})\,\delta(\mathbf{z}_{\rm c}-\mathbf{Z}_{\rm c}^{[i]})$ (23)
$\displaystyle F_{{\rm p},N_{\rm p}}^{\rm E}(t,\mathbf{x}\,;\,\mathbf{z}_{\rm
c})$ $\displaystyle=F_{{\rm p},N_{\rm p}}^{\rm L}(t\,;\,\mathbf{y}_{\rm
p}=\mathbf{x},\mathbf{z}_{\rm c})$ (24)

where $m_{\rm p}^{[i]}$ is the mass of the particle labeled $[i]$. This shows
that in a small volume around a point $\mathbf{x}$ where averages are obtained
from Monte Carlo estimations, that is as the ensemble averages over the
$N_{\rm p}^{(\mathbf{x})}$ particles present in that volume, we get the
equivalent of Favre, or mass-weighted, averages

$\langle H_{\rm p}\rangle(t,\mathbf{x})\simeq\langle H_{\rm p}\rangle_{N_{\rm
p}}=\frac{\sum_{i=1}^{N_{\rm p}^{(\mathbf{x})}}m_{\rm p}^{[i]}H_{\rm
p}(t;\mathbf{z}_{\rm c}^{[i]}(t))}{\sum_{i=1}^{N_{\rm p}^{(\mathbf{x})}}m_{\rm
p}^{[i]}}~{}.$ (25)

####  II.1.2 A unified approach to statistical operators

In Eq. (25), we retrieved the traditional ensemble averaging applied in the
Reynolds-averaged Navier-Stokes (RANS) equations in turbulence modeling, using
$F_{\rm p}^{\rm L}=\langle F_{{\rm p},N_{\rm p}}^{\rm L}\rangle\simeq F_{{\rm
p},N_{\rm p}}^{\rm L}$ for large $N_{\rm p}$. This is an echo of the all-
important role of the fine-grained PDF defined as
$\delta(\mathbf{Z}(t)-\mathbf{z})$, since we have
$p(t;\mathbf{z})=\langle\,\delta(\mathbf{Z}(t)-\mathbf{z})\,\rangle$. In
recent decades, another commonly encountered approach is the Large-Eddy
Simulation (LES) method in which a spatial filtering operation is applied
instead of an averaging one. For single-phase turbulent flows where we handle
fields, e.g. $Q(t,\mathbf{x})$, this means that we consider the statistical
operator noted $\langle\cdot\rangle_{\rm ls}$

$\langle Q\rangle_{\rm ls}(t,\mathbf{x})=\int
Q(t,\mathbf{x}^{\prime})\,G(\mathbf{x}^{\prime}-\mathbf{x})\,\mathrm{d}\mathbf{x}^{\prime}~{},$
(26)

where the spatial filter $G$ is taken as a spatially and temporally invariant
positive function with a compact support and such that $\int
G(\mathbf{x})\,\mathrm{d}\mathbf{x}=1$. In a series of papers, Pope and co-
workers [37, 101, 102] demonstrated that the PDF road can still be followed
provided that we consider the filtered density function (FDF) defined as
$p_{\rm
ls}(t;\mathbf{z})=\langle\,\delta(\mathbf{Z}(t)-\mathbf{z})\,\rangle_{\rm ls}$
(note that this makes sense since the position is always present in the state
vector and the decomposition $\mathbf{Z}=(\mathbf{X}_{\rm p},\mathbf{Z}_{\rm
c})$ is quite relevant here). This is referred to as the FDF approach. Its
original formulation consisted in handling field quantities only but, given
the upstream role of the Lagrangian MDF, this suggests to extend the FDF
approach to a system made up by either fields or particles by manipulating the
discrete Lagrangian MDF. In [72, section 7], this idea to derive FDF models
for LES from a Lagrangian description in terms of particles by starting from
the discrete Lagrangian Mass Density Function (LMDF) was developed.

We can go one step further and combine the RANS and LES descriptions in a
single formulation. Indeed, if we consider the discrete LMDF in Eq. (23) with
the same decomposition of the state vector,

$\widetilde{F}_{\rm p}^{\rm L}(t;\mathbf{y}_{\rm p},\mathbf{z}_{\rm
c})=\sum_{i=1}^{N_{\rm p}}m_{\rm p}^{[i]}\,\delta(\mathbf{y}_{\rm
p}-\mathbf{x}^{[i]}_{\rm p}(t))\,\delta(\mathbf{z}_{\rm
c}-\mathbf{Z}^{[i]}_{\rm c}(t))~{},$ (27)

we remark that the LMDF needed for RANS formulations and the Lagrangian
Filtered Density Function (LFDF) needed for LES are retrieved by applying the
different statistical operators, $\langle\cdot\rangle$ and
$\langle\cdot\rangle_{ls}$ respectively. This general approach based on the
discrete LMDF as the ‘parent function’ is represented in Fig. 4.

$\widetilde{F}_{\rm p}^{\rm L}(t;\mathbf{y}_{\rm p},\mathbf{z}_{\rm
c})$$F_{\rm p}^{\rm L}(t;\mathbf{y}_{\rm p},\mathbf{z}_{\rm
c})=\langle\widetilde{F}_{\rm p}^{\rm L}\rangle(t;\mathbf{y}_{\rm
p},\mathbf{z}_{\rm c})$$F_{\rm p,ls}^{\rm L}(t;\mathbf{y}_{\rm
p},\mathbf{z}_{\rm c})=\langle\widetilde{F}_{\rm p}^{\rm L}\rangle_{\rm
ls}(t;\mathbf{y}_{\rm p},\mathbf{z}_{\rm c})$$F_{\rm p}^{\rm
E}(t,\mathbf{x}_{\rm p};\mathbf{z}_{\rm c})=F_{\rm p}^{\rm
L}(t,\mathbf{y}_{\rm p}=\mathbf{x}_{\rm p},\mathbf{z}_{\rm c})$$F_{\rm
p,ls}^{\rm E}(t,\mathbf{x}_{\rm p};\mathbf{z}_{\rm c})=F_{\rm p,ls}^{\rm
L}(t,\mathbf{y}_{\rm p}=\mathbf{x}_{\rm p},\mathbf{z}_{\rm c})$RANS: mean-
field equationsLES: filtered-field
equations$\langle\cdot\rangle$$\langle\cdot\rangle_{\rm ls}$ Figure 4: Sketch
of the central role played by the discrete LMDF $\widetilde{F}_{\rm p}^{\rm
L}$ from which the application of different statistical operators yields the
Lagrangian and Eulerian MDFs and the resulting macroscopic descriptions.

The great benefit from such a unified formulation is that we do not have to
bother anymore about which framework we are evolving in (RANS/LES) or which
statistical averages to apply (ensemble or spatial filter). We can concentrate
on how the particle state vector $\mathbf{Z}_{\rm p}$ is represented by a
stochastic processes.

###  II.2 Stochastic processes

####  II.2.1 Definition and key characteristics

##### Definition and key characteristics

A stochastic process, noted as $\mathbf{Z}$ or $(\mathbf{Z}(t),t\geq 0)$, is a
family of random variables indexed by a parameter which is usually time (with
$t\in[0\,;\,T_{\text{max}}]$, where $T_{\text{max}}$ is an end time with the
possible value $T_{\text{max}}=+\infty$). In the following, we consider
vector-valued stochastic processes in $\Re^{d}$ (where $d$ is the dimension)
and denote the values in the corresponding sample space by $\mathbf{z}$. To
capture the key characteristics of a stochastic process, it is necessary to
resort to the mathematical definition according to which $\mathbf{Z}$ is a
family of measurable functions on an underlying probability space $\Omega$
equipped with a $\sigma$-algebra $\mathcal{F}$ and a probability measure
$\mathbb{P}$. This gives that $\mathbf{Z}$ as a measurable function of two
variables:

$\begin{gathered}[0\,;\,T_{\text{max}}]\times(\Omega,\mathcal{F},\mathbb{P})\longrightarrow(\Re^{d},\mathcal{\mathbf{B}}^{d})\\\
(t,\omega)\longmapsto\mathbf{Z}(t,\omega)~{},\end{gathered}$ (28)

where $\mathcal{\mathbf{B}}^{d}$ is the Borel $\sigma$-algebra for $\Re^{d}$
obtained as the $d$-tensorial-product of the Borel $\sigma$-algebra for $\Re$.

In some presentations, this rigorous definition is skipped and stochastic
processes are introduced directly in the image space
$(\Re^{d},\mathcal{\mathbf{B}}^{d})$. This turns out to be unfortunate. A
first reason is that the introduction of $\sigma$-algebras allows to quantify
the intuitive notions of the information content and of fine- or coarse-
grained descriptions which correspond to different embedded $\sigma$-algebras,
or sub-$\sigma$-algebras, embodying the information resolved by one
description. A second reason is that the correspondence between two ways to
characterize a stochastic process may be overlooked. On the one hand, we can
consider the family of vector-valued random variables $\mathbf{Z}(t)$ at fixed
times $t$. This corresponds to the ‘PDF point of view’ where the aim is to
derive the PDF equation satisfied by $p(t;\mathbf{z})$ in sample space. On the
other hand, we can consider a family of ‘elementary events’, represented by
different values of $\omega$ and handle a (large) number of time functions
$t\longmapsto\mathbf{Z}_{\omega}(t)$. This corresponds to the ‘trajectory
point of view’ where the objective is to write the time-evolution equations of
these ‘particles’ (samples of the PDF to be used in Monte Carlo estimations).
For a stochastic process, there is more information contained in its
trajectories than in its PDF. However, since we are interested in
approximating statistics from a stochastic process, we can refer to an
equivalence between the PDF and trajectory points of view in a weak sense.
Said otherwise, considering for instance Langevin type of equations in
physical space for a number of trajectories has the same status as the PDF
equation in sample space. This is why the trajectory point of view is often
adopted in the rest of the article.

##### The law of a stochastic process

Knowing a stochastic process is not equivalent to knowing the family of its
one-time PDFs $(p(t;\mathbf{z}),t\geq 0)$ since it implies knowing also the
joint laws of joint random variables $\mathbf{Z}(t)$ and
$\mathbf{Z}(t^{\prime})$ at any times $t$ and $t^{\prime}$, etc. Using a
discrete time setting for the sake of simplicity, this means that knowing the
law of a stochastic process is equivalent to the knowledge of the joint PDF
$p(t_{1},\mathbf{z}_{1};t_{2},\mathbf{z}_{2},\ldots;t_{N},\mathbf{z}_{N})$ for
any set of $N$ times and for any values of the chosen times
$(t_{1};t_{2},\ldots;t_{N})$. It is thus clear that the amount of information
required is huge and, in particular, much larger than the sole access to the
one-time PDF functions $p(t,\mathbf{z})$. In practice, we need to restrict
ourselves to more tractable sub-classes of stochastic processes, such as
Markov processes.

####  II.2.2 Markovian processes

In the world of ODEs, knowledge of an initial condition (at a time $t_{0}$)
and of the rate of change of a deterministic dynamical system is enough to
predict its future. A classical example is Lagrangian/Hamiltonian analytical
mechanics. The counterpart of this notion for probabilistic approaches leads
to define Markov processes for which knowledge of the present (in a
probabilistic sense) and of the law of conditional increments is enough to
predict the future (still in a probabilistic sense). To capture this idea, it
is sufficient to translate the Markov property in a discrete time setting,
which gives

$p(t_{n+1},\mathbf{z}_{n+1}\,|\,\left(t_{n},\mathbf{z}_{n};t_{n-1},\mathbf{z}_{n-1};\ldots;t_{1},\mathbf{z}_{1}\right))=p(t_{n+1},\mathbf{z}_{n+1}\,|\,t_{n},\mathbf{z}_{n})$
(29)

where $t_{n}$ represents the present time, $t_{n+1}$ the future and
$(t_{n-1},\ldots,t_{1})$ the past, while $\mathbf{z}_{i}$ is the value of the
process at $t=t_{i}$ (i.e. $\mathbf{Z}(t_{i})=\mathbf{z}_{i}$). An important
element is that the condition entering the conditional PDF, written as
$(t_{n},\mathbf{z}_{n})$ in Eq. (29), represents in fact the information known
at the present time $t=t_{n}$.

The fundamental property of Markov processes is that the law of the stochastic
process is determined from the knowledge of only two functions: the initial
PDF condition $p(t_{0},\mathbf{z}_{0})$ and the transition PDF
$p(t,\mathbf{z}\,|\,t_{0},\mathbf{z}_{0})$ from a state $\mathbf{z}_{0}$ at
$t=t_{0}$ to a state $\mathbf{z}$ at a later time $t$. Indeed, all the
$N$-time PDFs are reconstructed from the chain-rule [85, 35, 97]

$p(t_{n};\mathbf{z}_{n};t_{n-1},\mathbf{z}_{n-1};\cdots;t_{1},\mathbf{z}_{1};t_{0},\mathbf{z}_{0})=p(t_{n},\mathbf{z}_{n}\,|\,t_{n-1},\mathbf{z}_{n-1})\cdots
p(t_{1},\mathbf{z}_{1}\,|\,t_{0},\mathbf{z}_{0})p(t_{0},\mathbf{z}_{0})$ (30)

demonstrating that information on the complete law of the process is derived.
In particular, we obtain the Chapman-Kolmogorov equation

$p(t,\mathbf{z}\,|\,t_{0},\mathbf{z}_{0})=\int
p(t,\mathbf{z}\,|\,t_{1},\mathbf{z}_{1})\,p(t_{1},\mathbf{z}_{1}\,|\,t_{0},\mathbf{z}_{0})\,d\mathbf{z}_{1}~{},$
(31)

which forms the basis of path-integral formulations [114].

In a weak sense, Markov processes are characterized by the infinitesimal
operator

$\mathcal{L}_{t}g(\mathbf{z})=\lim_{dt\to
0}\frac{\langle\left(g(\mathbf{Z}_{t+dt})|\mathbf{Z}_{t}=\mathbf{z}\right)\rangle-g(\mathbf{z})}{dt}~{},$
(32)

which measures the effect of a conditional increment of the Markov process
over a test function $g$. It can then be shown [2, 35, 85] that the
transitional PDF $p(t,\mathbf{z}\,|\,t_{0},\mathbf{z}_{0})$ is the solution of
two equations depending on whether we fix the end condition $(t,\mathbf{z})$
or the initial one $(t_{0},\mathbf{z}_{0})$. In the first case, this
corresponds to the Kolmogorov backward equation written as

$\left\\{\begin{aligned}
&\frac{\partial\,p(t,\mathbf{z}\,|\,t_{0},\mathbf{z}_{0})}{\partial
t_{0}}+\mathcal{L}_{t_{0}}\left[\,p(t,\mathbf{z}\,|\,t_{0},\mathbf{z}_{0})\,\right]=0~{},\\\
&p(t,\mathbf{z}\,|\,t_{0},\mathbf{z}_{0})=\delta(\mathbf{z}-\mathbf{z}_{0})\quad
t_{0}\to t~{}.\end{aligned}\right.$ (33)

In the second case, the evolution equation corresponds to the Kolmogorov
forward equation

$\left\\{\begin{aligned}
&\frac{\partial\,p(t,\mathbf{z}\,|\,t_{0},\mathbf{z}_{0})}{\partial
t}=\mathcal{L}_{t}^{\bot}\left[\,p(t,\mathbf{z}\,|\,t_{0},\mathbf{z}_{0})\,\right]~{},\\\
&p(t,\mathbf{z}\,|\,t_{0},\mathbf{z}_{0})=\delta(\mathbf{z}-\mathbf{z}_{0})\quad
t\to t_{0}~{},\end{aligned}\right.$ (34)

where $\mathcal{L}_{t}^{\bot}$ is the adjoint operator of $\mathcal{L}_{t}$.
From a physical point of view, we often consider the time evolution of a
dynamical system from an observed initial state at a given time $t_{0}$ and we
are therefore mostly concerned with Kolmogorov forward equations but the
backward one is also of great interest (cf. the celebrated Feynman-Kac
formula). When dealing with Markov processes, these equations are central
since all the needed information is generated by the transition PDF. Note that
the same Kolmogorov forward equation is obtained for the one-time PDF
$p(t,\mathbf{z})$ by integrating over all initial conditions but actual
information is in the transition PDF. In practice, the subclass of Markov
stochastic processes is the only one for which we can develop a complete
description. For non-Markov processes, it is still possible to consider PDF
equations for the one-time PDF but it is already clear that the description of
the stochastic process is incomplete. It is shown in Sec. III that more
serious troubles lie ahead for non-Markov processes.

####  II.2.3 The two building blocks: the Wiener and Poisson processes

##### The Poisson process

For random discrete events, the reference model is the Poisson process where
there is a very small probability to have one event but each one implying a
change represented by a jump (multiple jumps at one time are not possible).
The trajectories of the Poisson process are therefore piecewise constant with
jumps (having a step of one unity in the standard Poisson process) occurring
at random times [35, 60] as illustrated in Fig. 5.

Figure 5: Some properties of the Poisson process $N(t)$: example of one
trajectory of the process jumping at random times $T_{i}$ which follows an
exponential distribution and the resulting pdf of $N(t)$ at each time $t$
being a Poisson distribution. Reprinted from [19] with permission from
Springer.

A Poisson process is characterized by its intensity $\lambda$ which is the
mean value of the number of events per unit time. The number of events in a
time interval $[t\,,\,t+\Delta t\,]$, represented by the stationary and
independent increments $\Delta N(t)=N(t+\Delta t)-N(t)$, is a Poisson random
variable

$\mathbb{P}[\,\Delta N(t)=k\,]=\frac{(\lambda\Delta
t)^{k}}{k!}e^{-\lambda\Delta t}$ (35)

from which it follows that the mean and variance of $\Delta N(t)$ are linear
with respect to $\Delta t$

$\langle\Delta N(t)\rangle=\lambda\Delta t\qquad\langle(\Delta
N(t)-\langle\Delta N(t)\rangle)^{2}\rangle=\lambda\Delta t.$ (36)

In any finite interval, the times at which the Poisson process jumps are
uniformly distributed. In practice, this property explains that deviations of
measured particle distributions from the Poisson law are regarded as
manifesting an underlying order. There are typically two ways to simulate a
Poisson process. The first method is based on waiting times (the time
intervals between successive random jumps) and, since values are constant
between jumps, the idea is to go directly from one event to the next one.
These waiting times are random variables following an exponential distribution
with the same intensity $\lambda$:

$\mathbb{P}[\,T_{i}=t\,]=\lambda e^{-\lambda t}\quad\Rightarrow\langle
T_{i}\rangle=\frac{1}{\lambda}~{},$ (37)

and are generated before applying the unit-one jump or more general events
when considering generalized Poisson processes (see below in section II.3).
This method is often used to simulate collision and agglomeration events in
molecular dynamics (cf. Sec. VI.1) but imposes variable time steps. The second
way is to generate the number of events occurring in fixed time intervals (cf.
Sec. VI.4). When the time interval $\Delta t$ is much smaller than the
relevant timescale $1/\lambda$ (i.e., $\lambda\,\Delta t\ll 1$), the
statistics of the increments, which follow the Poisson distribution given in
Eq. (35), have the simplified expression:

$\mathbb{P}[\,\Delta N(t)=0\,]\simeq 1-\lambda\Delta
t~{},\quad\mathbb{P}[\,\Delta N(t)=1\,]\simeq\lambda\Delta
t~{},\quad\mathbb{P}[\,\Delta N(t)=k\,]\simeq 0\;(k\geq 2)~{}.$ (38)

##### The Wiener process

The Wiener process is the canonical model of Brownian motion [35, 85] and can
be defined by the following three properties [55]:

  1. (i)

The process has independent increments, i.e. $(W(t_{3})-W(t_{2}))$ and
$(W(t_{1})-W(t_{0}))$ are independent when $t_{0}<t_{1}<t_{2}<t_{3}$;

  2. (ii)

The trajectories of the process are continuous functions (almost everywhere);

  3. (iii)

The increments of the Wiener process $(W(t_{2})-W(t_{1}))$ are centered
Gaussian random variables and with a variance equal to $(t_{2}-t_{1})$.

From these characteristics, it results that the Wiener process $W$ is a
Gaussian, Markov process, with zero mean and a covariance equal to $\langle
W(t)W(t^{\prime})\rangle=\min(t,t^{\prime})$. Its trajectories are continuous
and represent continuous evolutions where there is a near-one chance that
modifications happen but with small changes. The most important properties are
that the trajectories of the Wiener process are non-differentiable at any
point and have even unbounded total variations in any interval. Furthermore,
the increments $dW=W(t+dt)-W(t)$ are stationary and independent with
successive moments

$\langle
dW\rangle=0,\;\;\langle(dW)^{2}\rangle=dt,\;\;\langle(dW)^{n}\rangle=o(dt)\quad\forall
n>2~{}.$ (39)

The linear variation of $\langle(dW)^{2}\rangle$ with respect to $dt$ is of
paramount importance in the development of stochastic calculus [2, 60, 83].
One trajectory of a Wiener process is displayed in Fig. 6 illustrating the
wild variations at any scale.

Figure 6: One trajectory, or realization, of the Wiener process showing a
continuous but ragged profile. The window, with a zoom of the trajectory,
reveals the self-similar nature and the fast fluctuations of the trajectories
of the Wiener process at all scales.

###  II.3 Jump-diffusion processes and extended Fokker-Planck equations

To explain the shift from ODEs to Stochastic Differential Equations (SDEs), we
consider a one-dimensional dynamical system $Z(t)$ under the influence of
Gaussian white-noise $\zeta(t)$

$\frac{dZ(t)}{dt}=A(t,Z(t))+B(t,Z(t))\,\zeta(t)~{}.$ (40)

To arrive at a well-posed mathematical definition, the idea is to benefit from
the smoothing properties of the integration using the identification of
$\int_{0}^{t}\zeta(s)\,ds$ with a white-noise process, i.e.
$\zeta_{t}\,dt=dW(t)$. Instead of Eq. (40), we therefore try to give a meaning
to

$Z(t)=Z(t_{0})+\underbrace{\int_{t_{0}}^{t}A(s,Z(s))\,ds}_{\text{regular
integral}}+\underbrace{\int_{t_{0}}^{t}B(s,Z(s))\,dW(s)}_{\text{stochastic
integral}}.$ (41)

Due to the infinite total variation in any finite interval of the trajectories
of $W$, we cannot obtain the stochastic integral as the limit of classical
Riemann-Stieltjes sums. For a non-anticipating process $Z$, the stochastic
integral is then defined in the Ito sense as

$\int_{t0}^{t}B(s,Z(s))\,dW_{s}=\text{ms-}\lim_{N\to\infty}\sum_{i=1}^{N}B(t_{i},Z(t_{i}))\left(W(t_{i+1})-W(t_{i})\right)$
(42)

where $(t_{i})_{i=1,N}$ represents a partition of the interval $[0\,;\,t]$
(with $t_{1}=0$ and $t_{N+1}=t$) and where the limit is to be understood in
the mean-square sense [2, 85, 83] and not as a convergence trajectory by
trajectory. In the physics literature, these equations are written in
incremental form as a short-hand notation and are referred to as ‘Langevin
equations’

$dZ(t)=A(t,\,Z(t))\,dt+B(t,\,Z(t))\,dW(t)~{},$ (43)

while their corresponding PDF equation in sample space is the Fokker-Planck
(FP) equation

$\frac{\partial p(t,z\,|\,t_{0},z_{0})}{\partial
t}=-\frac{\partial[\,A(t,z)\,p(t,z\,|\,t_{0},z_{0})\,]}{\partial
z}+\frac{1}{2}\frac{\partial^{2}[\,B^{2}(t,z)\,p(t,z\,|\,t_{0},z_{0})\,]}{\partial
z^{2}}~{}.$ (44)

The FP equation is a forward Kolmogorov equation, Eq. (34), with the
infinitesimal operator $\mathcal{L}_{t}$ being

$\mathcal{L}_{t}[\cdot]=A(t,z)\frac{\partial\,[\,\cdot\,]}{\partial
z}+\frac{1}{2}B^{2}(t,z)\frac{\partial^{2}\,[\,\cdot\,]}{\partial z^{2}}~{}.$
(45)

Since the FP equation appears as a convection-diffusion equation, we speak of
stochastic diffusion processes and the functions $A$ and $B$ are referred to
as the drift and diffusion coefficients, respectively. Note that, regardless
of the sign of $B$, the diffusion coefficient in the PDF equation is $B^{2}$
and is always positive or null.

The correspondence between Langevin and Fokker-Planck equations is easily
extended to the multi-dimensional case. For instance, the SDEs written for a
$d$-dimensional process $\mathbf{Z}=(Z_{1},\ldots,Z_{d})$ are

$dZ_{i}(t)=A_{i}(t,\mathbf{Z}(t))\,dt+B_{ij}(t,\mathbf{Z}(t))\,dW_{j}(t)~{},$
(46)

with $(W_{j})_{j=1,d}$ a set of $d$ independent Wiener processes. In Eq. (46),
the drift $\mathbf{A}=(A_{i})$ is now a vector while the diffusion
$\mathbf{B}=(B_{ij})$ is a matrix. The corresponding Fokker-Planck equation is

$\frac{\partial p(t,\mathbf{z}\,|\,t_{0},\mathbf{z}_{0})}{\partial
t}=-\frac{\partial\left[\,A_{i}(t,\mathbf{z})\,p(t,\mathbf{z}\,|\,t_{0},\mathbf{z}_{0})\,\right]}{\partial
z_{i}}\\\
+\frac{1}{2}\frac{\partial^{2}\left[\,D_{ij}(t,\mathbf{z})\,p(t,\mathbf{z}\,|\,t_{0},\mathbf{z}_{0})\,\right]}{\partial
z_{i}\partial z_{j}}~{},$ (47)

where $\mathbf{D}=\mathbf{B}\mathbf{B}^{\bot}$ (or $D_{ij}=B_{ik}B_{jk}$) is a
symmetric definite-positive matrix. As in the one-dimensional case, the
positivity of $\mathbf{D}$ justifies the reference to a diffusive nature
induced by the rapidly-varying terms, $B_{ij}(t,\mathbf{Z}(t))\,dW_{j}(t)$, in
Eq. (46). The fact that several matrices $\mathbf{B}$ can correspond to the
same diffusion matrix $\mathbf{D}$ translates the equivalence between the
trajectory and PDF points of view in a weak sense.

The physical meaning of $A$ and $B$ in the Langevin and Fokker-Planck
equations are revealed by considering the statistics of the conditional
increments (using a one-dimensional version)

$\displaystyle\langle\Delta Z\,|\,Z(t)=z\rangle$ $\displaystyle=A(t,z)\,\Delta
t~{},$ (48a) $\displaystyle\langle(\Delta Z)^{2}\,|\,Z(t)=z\rangle$
$\displaystyle=B^{2}(t,z)\,\Delta t~{}.$ (48b)

As shown in Fig. 7, the drift term, $A(t,z)$, governs the mean evolution of
the conditional increments while the diffusion coefficient, $B(t,z)$, governs
the spread of the conditional increments around its mean value.

Figure 7: Statistics over small time increments on the conditional
trajectories of a stochastic diffusion process illustrating the physical
meaning of the drift and diffusion coefficients.

There is regularly some confusion about the Gaussian hypothesis in SDEs. As
expressed by Eqs. (48), the Gaussian hypothesis applies only to the
conditional increments and, in the general situation when $B(t,z)$ is not
constant, the resulting process $Z$ can deviate from Gaussianity. Though this
point has been clarified [76, 72], it is worth emphasizing that the
observation of non-Gaussian distributions does not invalid models based on
stochastic diffusion processes.

Apart from the case of Brownian motion itself (where $A=0$ and $B=1$), the
simplest example of a stochastic diffusion process is the Ornstein-Uhlenbeck
(OU) process whose trajectories are

$dZ(t)=-\frac{Z(t)}{\tau}\,dt+K\,dW(t)~{},$ (49)

where the timescale $\tau$ and the diffusion coefficient $K$ are constant.
From Ito calculus [35], it is trivial to show that $K$, $\tau$ and $\langle
Z^{2}\rangle$ are related through $K^{2}=2\langle Z^{2}\rangle/\tau$ which
corresponds to a fluctuation-dissipation theorem [86]. Another immediate
property is that the auto-correlation function of the stationary process
$\rho(t^{\prime}-t)=\langle Z(t)Z(t^{\prime})\rangle/\langle Z^{2}\rangle$ is
$\rho(t^{\prime}-t)=\exp(-(t^{\prime}-t)/\tau)$, revealing that the relaxation
timescale $\tau$ is also the integral timescale since
$\tau=\int_{0}^{+\infty}\rho(s)\,\mathrm{d}s$. This is the first hint that
timescales play a central role when formulating stochastic models, as
exemplified in later sections, in particular Sec. IV.

So far, we have built SDEs based on the Wiener process for diffusion processes
with continuous trajectories. Similar steps can be made with the Poisson
process to introduce sudden jumps. Still using the trajectory point of view,
the SDEs for a jump-diffusion process are written for a one-dimensional
process

$dZ(t)=A(t,Z(t))\,dt+B(t,Z(t))\,dW(t)+C(t,Z(t))\,dN(t)~{},$ (50)

where $N$ a Poisson process with intensity $\lambda$ and $C(t,Z(t))$
represents the amplitude of the jumps. The jumps contribute to the statistics
of the conditional increments over a time step $\Delta t$ and Eqs. (48) become

$\displaystyle\langle\Delta Z\,|\,Z(t)=z\rangle=A(t,z)\,\Delta
t+C(t,z)\lambda\,\Delta t~{},$ (51a) $\displaystyle\langle(\Delta
Z)^{2}\,|\,Z(t)=z\rangle=B^{2}(t,z)\,\Delta t+C^{2}(t,z)\lambda\,\Delta t~{}.$
(51b)

The introduction of jumps is illustrated in Fig. 8, revealing the
discontinuous nature of the trajectories and their diffusive behavior between
successive jumps.

Figure 8: Trajectories of a jump-diffusion process with their diffusive
behavior in between sudden jumps.

From the properties of the Poisson process, we know that, over an
infinitesimal time increment $dt$ and retaining only the first order
contributions, the increments of $N(t)$ can take only two possible values

$\mathbb{P}[\,dN(t)=k\,]=(1-\lambda
dt)\,\delta_{k,0}+\lambda\,dt\,\delta_{k,1}\Rightarrow\langle(dN(t))^{m}\rangle=\lambda\,dt\quad\forall
m.$ (52)

The Poisson jumps contribute to any order in $dt$ (whereas the increments of
the Wiener process contributes only to the second order in $dt$). The SDEs of
a jump-diffusion process can be generalized by considering random amplitudes
for the jumps (this is referred to as a compound Poisson process), which gives

$dZ(t)=A(t,Z(t))\,dt+B(t,Z(t))\,dW(t)+C(t,Z(t),Q)\,dN(t)~{},$ (53)

where $Q$ is an independent random variable. The law of the random jumps
becomes then

$W(y\,|\,t,z)=W(y\,|\,Z(t)=z)=\mathbb{P}[\,C=y\,|\,Z(t)=z\,]\,\lambda~{},$
(54)

where the probability $\mathbb{P}$ refers to the law of the independent random
variable $Q$. In sample space, the equation for the transitional PDF
$p(t,z\,|\,t_{0},z_{0})$ is

$\frac{\partial p(t,z\,|\,t_{0},z_{0})}{\partial
t}=-\frac{\partial[\,A(t,z)\,p(t,z\,|\,t_{0},z_{0})\,]}{\partial
x}+\frac{1}{2}\frac{\partial^{2}[\,B^{2}(t,z)\,p(t,z\,|\,t_{0},z_{0})\,]}{\partial
x^{2}}\\\
+\int[\,W(z\,|\,t,y)\,p(t,y\,|\,t_{0},z_{0})-W(y\,|\,t,z)\,p(t,z\,|\,t_{0},z_{0})\,]\,dy.$
(55)

This form is referred to as the Chapman-Kolmogorov PDF equation in [35]. Jump-
diffusion processes resurface in Sec. VI.4 to account for diffusive effects
due to random continuous motions of a turbulent flow and sudden velocity jumps
due to particle collisions.

###  II.4 A double hierarchy: one- or N-particle PDF and particle-attached
variables

Having set the probabilistic framework, we return to the first step mentioned
at the beginning of this section. Yet, before considering their modeled
evolution, we need to define the variables in the state vector. This includes:
(a) choosing a one- or a $N$-particle PDF formulation; (b) determining the
relevant variables to describe each particle. Point (a) is a very classical
issue in statistical physics in relation with the well-known BBGKY hierarchy
for any system with (at least) two-particle interactions. At the moment, two-
particle PDF models (not to mention $N$-particle ones, with $N\geq 3$) are not
developed enough for general inhomogeneous turbulent flows, especially wall-
bounded ones, and are therefore far less tractable [91, 89]. For these reasons
and though references are made to such approaches (cf. Secs. IV and VI), we
retain the one-particle PDF level of description and, correspondingly, a one-
point Eulerian PDF description of fields. It is then essential to be aware of
the inherent loss of information since, for instance, two-point correlations
cannot be extracted from one-particle PDF models and particle-particle
interactions, such as actions at distance or even collisions, have to be
modeled through mean-field formulations. Point (b) turns out to be a challenge
to traditional statistical views and requires specific analysis. This is
addressed in Sec. III.

##  III Statistical modeling of particle transport in turbulent flows

###  III.1 Shortcomings of kinetic-like descriptions

The Boltzmann picture is dominating the description of transport phenomena to
such an extent that the kinetic framework is rarely questioned. Yet, whether
kinetic variables are necessarily the only ones to retain in the particle
state vector is not a point over which to pass too quickly. For discrete
particles in turbulent flows, there are typically three situations. The first
one corresponds to the case of fully-resolved turbulent flows, where the fluid
velocity field is known at every point and every time. This means that the
velocity of the fluid seen is similar to an external deterministic force
field. It is then accounted for without approximation in Boltzmann-like PDF
models, i.e. in PDF formulations based on the kinetic state vector
$\mathbf{Z}_{\rm p}^{\rm r}=(\mathbf{X}_{\rm p},\mathbf{U}_{\rm p})$ where we
handle a PDF noted $p^{\rm r}(t;\mathbf{z}_{\rm p}^{\rm r})=p^{\rm
r}(t;\mathbf{y}_{\rm p},\mathbf{V}_{\rm p})$ in sample space. In these
notations, the superscript r indicates a reduced description compared to
$\mathbf{Z}_{\rm p}=(\mathbf{X}_{\rm p},\mathbf{U}_{\rm p},\mathbf{U}_{\rm
s})$ and the PDF $p^{\rm r}$ is the marginal of $p(t;\mathbf{z}_{\rm
p})=p(t;\mathbf{y}_{\rm p},\mathbf{V}_{\rm p},\mathbf{V}_{\rm s})$. The second
situation corresponds to high-inertia particles for which the underlying fluid
can be regarded as white-noise, leading to a Langevin equation for the
discrete particle velocity and a resulting FP equation for $p^{\rm r}$. The
third situation corresponds to the usual case of non-fully-resolved turbulent
flows, where $\mathbf{U}_{\rm s}$ is neither deterministic nor white-noise but
exhibits memory effects. This is the situation addressed by the kinetic-PDF
model.

By manipulating the fine-grained PDF $\mathcal{P}(t;\mathbf{y}_{\rm
p},\mathbf{V}_{\rm p})=\delta(\mathbf{X}_{\rm p}(t)-\mathbf{y}_{\rm
p})\delta(\mathbf{U}_{\rm p}(t)-\mathbf{V}_{\rm p})$ and using standard
techniques [89] with $p^{\rm r}(t;\mathbf{y}_{\rm p},\mathbf{V}_{\rm
p})=\langle\mathcal{P}(t;\mathbf{y}_{\rm p},\mathbf{V}_{\rm p})\rangle$, the
kinetic PDF equation is derived from the particle equations keeping only the
drag force which is the central point. The unclosed PDF equation is

$\frac{\partial p^{\rm r}}{\partial t}+\frac{\partial\left[\,V_{{\rm
p},i}\,p^{\rm r}\,\right]}{\partial y_{{\rm
p},i}}=\left[\,\left\langle\frac{V_{{\rm p},i}}{\tau_{\rm
p}}\,|\,(\mathbf{y}_{\rm p},\mathbf{V}_{\rm p})\right\rangle p^{\rm
r}\,\right]-\frac{\partial}{\partial V_{{\rm
p},i}}\left[\,\left\langle\frac{U_{{\rm s},i}}{\tau_{\rm
p}}\,|\,(\mathbf{y}_{\rm p},\mathbf{V}_{\rm p})\right\rangle p^{\rm
r}\,\right]~{}.$ (56)

It first appears that the complete expression of the particle relaxation time
$\tau_{\rm p}$ cannot be accounted for at this level of description, since it
is usually a function of $(\mathbf{U}_{\rm p},\mathbf{U}_{\rm s})$. In the
kinetic approach, it is assumed that we are only dealing with particles in the
Stokes regime so that $\tau_{\rm p}\simeq\tau_{\rm p}^{\rm St}$ becomes
independent of $(\mathbf{U}_{\rm p},\mathbf{U}_{\rm s})$, which is already a
limitation. With this approximation, the unclosed kinetic PDF equation is
written as

$\frac{\partial p^{\rm r}}{\partial t}+\frac{\partial\left[\,V_{{\rm
p},i}\;p^{\rm r}\,\right]}{\partial y_{{\rm p},i}}=\frac{\partial}{\partial
V_{{\rm p},i}}\left[\,\frac{V_{{\rm p},i}}{\tau_{p}^{\rm St}}\,p^{\rm
r}\,\right]-\frac{\partial}{\partial V_{{\rm
p},i}}\left[\,\frac{1}{\tau_{p}^{\rm St}}\langle U_{{\rm f},i}\rangle\,p^{\rm
r}\,\right]-\frac{\partial}{\partial V_{{\rm
p},i}}\left[\,\frac{1}{\tau_{p}^{\rm St}}\langle u^{{}^{\prime}}_{{\rm
s},i}\,|\,(\mathbf{y}_{\rm p},\mathbf{V}_{\rm p})\rangle\,p^{\rm r}\,\right]$
(57)

where the velocity of the fluid seen has been decomposed as $\mathbf{U}_{\rm
s}=\langle\mathbf{U}_{\rm f}\rangle(t,\mathbf{X}_{\rm p}(t))+\mathbf{u}_{\rm
s}^{{}^{\prime}}$. This does not imply, however, that $\langle\mathbf{u}_{\rm
s}^{{}^{\prime}}\rangle=0$ since the set of velocities of the fluid seen
represents only a subset of fluid particle velocities at a given location
(unless the so-called well-mixed condition is satisfied), and its mean value
is called the drift velocity $\mathbf{U}_{\rm d}$, which is thus equal to
$\mathbf{U}_{\rm d}=\langle\mathbf{U}_{\rm s}\rangle-\langle\mathbf{U}_{\rm
f}\rangle(t,\mathbf{X}_{\rm p}(t))$.

Closure of the open flux in sample space is derived from Furutsu-Novikov-
Donsker (FND) relation [34, 82, 23], which expresses the correlation between a
Gaussian field with zero mean and an arbitrary functional of that field. It is
applied by assuming that the fluctuating velocity field
$\mathbf{u}_{f}(t,\mathbf{x})$ is a random Gaussian field and writes for a
functional $F[t;\mathbf{u}_{f}]$ of $\mathbf{u}_{f}(t,\mathbf{x})$

$\langle u_{{\rm f},i}(t,\mathbf{x})\,F[t;\mathbf{u}_{\rm
f}]\rangle=\int_{0}^{t}\int_{\mathbf{x}^{\prime}}R_{ik}(t,\mathbf{x};t^{\prime},\mathbf{x}^{\prime})\left\langle\frac{\delta
F[t;\mathbf{u}_{\rm f}]}{\delta u_{{\rm
f},k}(t^{\prime},\mathbf{x}^{\prime})}\right\rangle\mathrm{d}\mathbf{x}^{\prime}\,\mathrm{d}t^{\prime}$
(58)

where $R_{ik}(t,\mathbf{x};t^{\prime},\mathbf{x}^{\prime})=\langle u_{{\rm
f},i}(t,\mathbf{x})\,u_{{\rm f},k}(t^{\prime},\mathbf{x}^{\prime})\rangle$ is
the fluid two-point two-time correlation. Then, applying the FND relation to
the fine-grained PDF $\mathcal{P}$ yields the expression of the flux closure
as [14]

$\frac{1}{\tau_{p}^{\rm St}}\langle U_{{\rm s},i}\,|\,\mathbf{z}^{r}_{\rm
p}\rangle\,p^{r}=\frac{1}{\tau_{\rm p}^{\rm St}}\langle U_{{\rm
f},i}\rangle(t,\mathbf{y}_{\rm p})\,p^{r}+\frac{1}{\tau_{\rm p}^{\rm
St}}U_{{\rm d},i}(t,\mathbf{y}_{\rm
p})\,p^{r}-\lambda_{ij}\,\frac{\partial\left[\,p^{r}\,\right]}{\partial
y_{{\rm p},j}}-\mu_{ij}\,\frac{\partial\left[\,p^{r}\,\right]}{\partial
V_{{\rm p},j}}$ (59)

with $\lambda_{ij}$ and $\mu_{ij}$ given by

$\displaystyle\lambda_{ij}$ $\displaystyle=\frac{1}{\tau_{\rm p}^{\rm
St}}\int_{0}^{t}\langle\Gamma_{jk}(t,t^{\prime})\,R_{ik}(t,\mathbf{y}_{\rm
p};t^{\prime},\mathbf{X}_{\rm p}(t^{\prime}))\rangle_{(\mathbf{y}_{\rm
p},\mathbf{V}_{\rm p})}\,dt^{\prime}$ (60a) $\displaystyle\mu_{ij}$
$\displaystyle=\frac{1}{\tau_{\rm p}^{\rm
St}}\int_{0}^{t}\langle\dot{\Gamma}_{jk}(t,t^{\prime})\,R_{ik}(t,\mathbf{y}_{\rm
p};t^{\prime},\mathbf{X}_{\rm p}(t^{\prime}))\rangle_{(\mathbf{y}_{\rm
p},\mathbf{V}_{\rm p})}\,dt^{\prime}$ (60b)

where the notation $\langle.\rangle_{(\mathbf{y}_{\rm p},\mathbf{V}_{\rm p})}$
indicates the averaged value conditioned on the particle trajectory that
‘arrives’ at $(\mathbf{y}_{\rm p},\mathbf{V}_{\rm p})$ at time $t$, which
explains that the dispersion tensors are functions of $(\mathbf{y}_{\rm
p},\mathbf{V}_{\rm p})$ even if the eventual dependence on $\mathbf{V}_{\rm
p}$ is often neglected. In these equations, $\Gamma_{jk}(t,t^{\prime})$ stands
for the response function

$\Gamma_{jk}(t,t^{\prime})=\frac{\delta X_{{\rm p},j}(t)}{\delta u_{{\rm
f},k}(t^{\prime},\mathbf{X}_{\rm p}(t^{\prime}))}$ (61)

that measures the effect of a perturbation of the fluctuating fluid velocity
seen at an earlier time $t^{\prime}$ on the particle position $\mathbf{X}_{\rm
p}(t)$ at time $t$, and $\dot{\Gamma}=\partial\Gamma/\partial t$. With these
expressions, the flux can also be written as

$\frac{1}{\tau_{\rm p}^{\rm St}}\langle U_{{\rm s},i}\,|\,\mathbf{z}_{\rm
p}\rangle\,p^{r}=\frac{1}{\tau_{\rm p}^{\rm St}}\langle U_{{\rm
f},i}\rangle(t,\mathbf{y}_{\rm
p})\,p^{r}+\kappa_{i}\,p^{r}-\frac{\partial\left[\,\lambda_{ij}\,p^{r}\,\right]}{\partial
y_{{\rm p},j}}-\frac{\partial\left[\,\mu_{ij}\,p^{r}\,\right]}{\partial
V_{{\rm p},j}}$ (62)

with

$\kappa_{i}(t;\mathbf{z}_{\rm p})=\frac{1}{\tau_{\rm p}^{\rm St}}U_{{\rm
d},i}(t,\mathbf{y}_{\rm p})+\frac{\partial\lambda_{ij}}{\partial y_{{\rm
p},j}}+\frac{\partial\mu_{ij}}{\partial V_{{\rm p},j}}~{}.$ (63)

Inserting Eq. (62) into Eq. (57), we obtain the kinetic PDF equation which, in
a compact form, is

$\frac{\partial p^{\rm r}}{\partial t}=-\frac{\partial}{\partial z^{\rm
r}_{l}}\left[A^{\rm KE}_{l}\,p^{\rm
r}\right]+\frac{1}{2}\frac{\partial^{2}}{\partial z^{\rm r}_{l}\partial z^{\rm
r}_{m}}\left[B^{\rm KE}_{lm}\,p^{\rm r}\right]~{},$ (64)

where the components of the drift vector $(A^{\rm KE}_{l})_{l=1,6}$ are

$A^{\rm KE}_{l}=\begin{cases}V_{p,l}&\;l=1,3,\\\ \dfrac{\langle U_{{\rm
f},l-3}\rangle-V_{{\rm p},l-3}}{\tau_{\rm p}^{\rm
St}}+\kappa_{l-3}&\;l=4,6~{},\end{cases}$ (65)

and where the symmetrical matrix $B^{\rm KE}_{lm}$ is (using a bloc notation
with $i,j=1,3$)

$\mathbf{B}^{\rm KE}=\begin{pmatrix}0&\vline&(\lambda_{ij})\\\ \--&&\---\\\
(\lambda_{ji})&\vline&(\mu_{ij})+(\mu_{ji})\end{pmatrix}~{}.$ (66)

As such, Eq. (64) looks like a classical convection-diffusion equation with
$B^{\rm KE}_{lm}$ in the second-order derivative appearing as a diffusion
matrix. However, it is immediate to show that this symmetrical matrix $B^{\rm
KE}_{lm}$ is negative definite. Indeed, its determinant is
$det(\mathbf{B}^{\rm KE})=-\left(det(\bm{\lambda})\right)^{2}<0$ and,
consequently, the matrix $\mathbf{B}^{\rm KE}$ has always at least one
negative eigenvalue. This point was first put forward in [94] where it was
associated with an ‘anti-diffusive’ behavior and the non-Markovian nature of
the reduced state vector $\mathbf{Z}^{\rm r}_{\rm p}$. It was later discussed
in a comprehensive analysis of kinetic- and dynamic-PDF models in [78] to
which readers are referred to for further details. The unavoidable consequence
is that this kinetic-PDF equation cannot be solved, unless for very special
initial conditions, and, therefore, does not qualify as an acceptable PDF
model.

Is the failure of kinetic descriptions only due to specific closures or to the
sole situation of small discrete particles in non-fully-resolved turbulent
flows? Or is it related to a poor choice of the variables entering the
particle state vector and the loss of Markovianity? In relation to these
queries is the remark that, when considering homogeneous situations (for the
sake of avoiding more complex notations in general flows, see in-depth
discussions of these relations in [78]), $\lambda_{ij}$ and $\mu_{ij}$ are
simply expressed as the correlations between particle positions and velocities
and the variable that is eliminated from the reduced PDF description, that is
the velocity of the fluid seen:

$\displaystyle\lambda_{ij}$ $\displaystyle=\frac{1}{\tau_{\rm p}^{\rm
St}}\left\langle u_{{\rm s},i}(t)\,X_{{\rm p},j}(t)\right\rangle~{},$ (67a)
$\displaystyle\mu_{ij}$ $\displaystyle=\frac{1}{\tau_{\rm p}^{\rm
St}}\left\langle u_{{\rm s},i}(t)\,U_{{\rm p},j}(t)\right\rangle~{}.$ (67b)

This raises the question if leaving out $\mathbf{U}_{\rm s}$ was appropriate.
These issues are now investigated.

###  III.2 Slow and fast variables and the Markovian approach

####  III.2.1 Colored or white noises and well-posed PDF equations

The failure of the kinetic description for discrete particles in non-fully-
resolved turbulent flows can be cast in the more general framework of
dynamical systems under the influence of ‘external noises’ [39, 107, 54, 106].
To shed light on this issue, we consider a dynamical system characterized by a
state vector $\mathbf{Z}(t)=(Z_{1}(t),Z_{2}(t),\ldots,Z_{n}(t))$ with
evolution equations

$\frac{dZ_{i}}{dt}=A_{i}(t,\mathbf{Z})+B_{ij}(t,\mathbf{Z})\,\xi_{j}$ (68)

where $\bm{\xi}(t)=(\xi_{j}(t))_{j=1,n}$ represents an ‘external noise’. For
the sake of simplicity since it does not affect results to come, the
dependence of $\mathbf{A}$ and $\mathbf{B}$ is simply written as
$\mathbf{A}(t,\mathbf{Z})$ and $\mathbf{B}(t,\mathbf{Z})$ instead of a more
general dependence including statistics of the stochastic process (such as its
mean, variance, etc.).

When the external noise $\bm{\xi}$ is a Gaussian process with independent
values, we are evolving within the well-established framework of Ito SDEs and
Fokker-Planck equations, cf. Sec. II.3. The situation is different for
‘colored noises’, that is when $\bm{\xi}$ has a non-zero correlation
timescale. This is found, for example, for stationary Gaussian processes
simulated as a set of independent Ornstein-Uhlenbeck (OU) processes

$d\xi_{k}=-\frac{\xi_{k}}{\tau}\,dt+\sigma\,dW_{k}~{},$ (69)

where $\tau$ is the timescale of $\xi_{k}$ and $\sigma$ a constant equal to
$\sigma^{2}=2\langle\xi_{k}^{2}\rangle/\tau$ from the classical fluctuation-
dissipation theorem [35]. For the sake of simplicity, $\tau$ and $\sigma$ are
retained for each component $\xi_{k}$ (differences can be accounted for
through the choice of $\mathbf{B}$). The auto-correlation of this process is
an exponential function
$\langle\xi_{k}(t)\xi_{k}(s)\rangle=\langle\xi^{2}\rangle\exp(-|t-s|/\tau)$
and is therefore not delta-correlated when $\tau\neq 0$. An important element
is that for such external noises, the process $\mathbf{Z}$ in Eq. (68) is no
longer Markovian [54, 106, 97] although a classical remark is to note that
Markovianity is retrieved by considered the extended process
$(\mathbf{Z},\bm{\xi})$ [54, 106, 97]. Even in a non-Markovian case, we can
still consider the equation satisfied by the one-time PDF of the process
$p(t;\mathbf{z})$, from which a set of PDEs are derived for some relevant
statistical moments. It is therefore essential that this equation be
mathematically well-posed for the resulting macroscopic descriptions to be
regarded as acceptable, regardless of the choices of $\mathbf{A}$ and
$\mathbf{B}$ in Eq. (68).

##### PDF equation and the well-posed criterion

The PDF equation for Gaussian colored noise is derived by standard techniques
[89] from the fine-grained PDF
$\mathcal{P}(t;\mathbf{z})=\delta(\mathbf{Z}(t)-\mathbf{z})$ with
$p(t;\mathbf{z})=\langle\mathcal{P}(t;\mathbf{z})\rangle$. The exact but
unclosed PDF equation for $p(t;\mathbf{z})$ is

$\frac{\partial p}{\partial t}=-\frac{\partial}{\partial
z_{i}}\left[A_{i}(t,\mathbf{z})\,p\,\right]-\frac{\partial}{\partial
z_{i}}\left[B_{ik}(t,\mathbf{z})\langle\xi_{k}\,|\,(t,\mathbf{z})\rangle\,p\,\right]$
(70)

which, using $\langle X\mathcal{P}\rangle=\langle X|(t,\mathbf{z})\rangle
p(t,\mathbf{z})$, can be written as

$\frac{\partial p}{\partial t}=-\frac{\partial}{\partial
z_{i}}\left[A_{i}(t,\mathbf{z})\,p\,\right]-\frac{\partial}{\partial
z_{i}}\left[B_{ik}(t,\mathbf{z})\langle\xi_{k}\mathcal{P}\rangle\right]~{}.$
(71)

In Eq. (71), the open flux $\langle\xi_{k}\mathcal{P}\rangle$ can be written
as $\langle\xi_{k}\mathcal{F}[t;\bm{\xi}]\rangle$, where $\mathcal{F}[\cdot]$
stands for a functional dependence since $\mathcal{P}$ can be seen as a
functional of the Gaussian centered process $\bm{\xi}$. As for derivation of
the kinetic PDF, we can apply the Furutsu-Novikov-Donsker (FND) relation [34,
82, 23]

$\langle\xi_{k}(t)F[t;\bm{\xi}]\rangle=\int_{0}^{t}\langle\xi_{k}(t)\xi_{l}(t^{\prime})\rangle\left\langle\frac{\delta
F[t;\bm{\xi}]}{\delta\xi_{l}(t^{\prime})}\right\rangle dt^{\prime}.$ (72)

Then, using

$\displaystyle\frac{\delta\mathcal{P}(t;\mathbf{z})}{\delta\xi_{l}(t^{\prime})}$
$\displaystyle=-\frac{\delta
Z_{j}(t)}{\delta\xi_{l}(t^{\prime})}\frac{\partial\mathcal{P}(t;\mathbf{z})}{\partial
z_{j}}$ (73a) $\displaystyle=-\frac{\partial}{\partial
z_{j}}\left[\frac{\delta
Z_{j}(t)}{\delta\xi_{l}(t^{\prime})}\mathcal{P}\right]+\left(\frac{\partial}{\partial
z_{j}}\left[\frac{\delta
Z_{j}(t)}{\delta\xi_{l}(t^{\prime})}\right]\right)\mathcal{P}$ (73b)

and applying the averaging operator, we obtain

$\langle\xi_{k}\mathcal{P}\rangle=\alpha_{k}(t,\mathbf{z})\,p(t;\mathbf{z})-\frac{\partial\left[\,\lambda_{kj}(t,\mathbf{z})\,p(t;\mathbf{z})\,\right]}{\partial
z_{j}}$ (74)

with

$\alpha_{k}(t,\mathbf{z})=\int_{0}^{t}\langle\xi_{k}(t)\xi_{l}(t^{\prime})\rangle\left\langle\frac{\partial}{\partial
z_{j}}\left[\frac{\delta
Z_{j}(t)}{\delta\xi_{l}(t^{\prime})}\right]|(t,\mathbf{z})\right\rangle
dt^{\prime}$ (75)

and where $\lambda_{kj}$ is given by

$\lambda_{kj}(t,\mathbf{z})=\int_{0}^{t}\langle\xi_{k}(t)\xi_{l}(t^{\prime})\rangle\left\langle\frac{\delta
Z_{j}(t)}{\delta\xi_{l}(t^{\prime})}|(t,\mathbf{z})\right\rangle
dt^{\prime}~{}.$ (76)

On the other hand, applying directly the FND relation to $Z_{j}(t)$ in Eq.
(72) leads to

$\langle\xi_{k}Z_{j}(t)\rangle=\int\lambda_{kj}(t,\mathbf{z})p(t;\mathbf{z})d\mathbf{z}=\int_{0}^{t}\langle\xi_{k}(t)\xi_{l}(t^{\prime})\rangle\left\langle\frac{\delta
Z_{j}(t)}{\delta\xi_{l}(t^{\prime})}\right\rangle dt^{\prime}~{}.$ (77)

Combining Eqs. (74) and Eq. (71) gives the closed form for the one-time PDF
equation

$\frac{\partial p}{\partial t}=-\frac{\partial}{\partial
z_{i}}\left[\,(A_{i}+B_{ik}\alpha_{k})\,p\,\right]+\frac{\partial}{\partial
z_{i}}\left[B_{ik}\frac{\partial\left(\lambda_{kj}\,p\right)}{\partial
z_{j}}\right]~{},$ (78)

which can be re-arranged as

$\frac{\partial p}{\partial t}=-\frac{\partial}{\partial
z_{i}}\left[\widetilde{A}_{i}\,p\right]+\frac{1}{2}\frac{\partial^{2}}{\partial
z_{i}\partial z_{j}}\left[\widetilde{D}_{ij}\,p\right]~{},$ (79)

with $\widetilde{A}_{i}=A_{i}+B_{ik}\alpha_{k}+\lambda_{kj}\partial
B_{ik}/\partial z_{j}$ and $\widetilde{D}_{ij}$ the symmetrical matrix

$\widetilde{D}_{ij}(t,\mathbf{z})=B_{ik}(t,\mathbf{z})\lambda_{kj}(t,\mathbf{z})+B_{jk}(t,\mathbf{z})\lambda_{ki}(t,\mathbf{z})~{}.$
(80)

The well-posed nature of the PDF equation, Eq. (79), is determined by the
matrix $\widetilde{D}_{ij}$ in Eq. (80). Indeed, if $\widetilde{D}_{ij}$ is
not positive definite (there is at least one negative eigenvalue), this
implies that a marginal of the one-time PDF appears as the ‘solution’ of an
‘anti-diffusion’ PDE which are ill-posed in the sense that they can only be
solved for very special initial conditions. Hence, the well-posed criterion is
to require that $\widetilde{D}_{ij}(t,\mathbf{z})$ have only positive (or
null) eigenvalues, $\forall\mathbf{z}$.

##### Analysis of the linear case

In [73, section 9.2], the well-posed property of such PDF equations is
investigated, especially for a dynamical system with a general structure
containing kinetic descriptions, and it is shown that the positive nature of
$D_{ij}$ is only obtained when taking the white-noise limit. In the present
context, it is sufficient to consider a simpler situation where the drift
vector is linear in $\mathbf{Z}(t)$ (or linearized around a given point in
sample space). Then, Eq. (68) becomes

$\frac{dZ_{i}}{dt}=-G_{ik}Z_{k}+B_{ik}\xi_{k}$ (81)

where $\mathbf{G}$ is a constant matrix representing return-to-equilibrium
effects and where the colored noise $\bm{\xi}$ is a set of independent
stationary OU processes as in Eq. (69). In the linear case, the ‘response
functions’ $\delta Z_{i}(t)/\delta\xi_{l}(t^{\prime})$ in Eq. (76) are
independent of the sample space value $\mathbf{z}$, showing that
$\lambda_{ji}=\langle Z_{i}\,\xi_{j}\rangle$. In the stationary state where
$\lambda_{ji}$ reach constant values, the correlations are easily derived
through

$\frac{d\lambda_{ji}}{dt}=0\Longrightarrow\left(\delta_{ik}+\tau
G_{ik}\right)\langle Z_{k}\xi_{j}\rangle=B_{ij}\tau\langle\xi^{2}\rangle~{},$
(82)

which can be inverted to give
$\lambda_{ji}=\tau\langle\xi^{2}\rangle\widetilde{G}^{-1}_{ik}B_{kj}$ with
$\widetilde{\mathbf{G}}=\mathbb{1}+\tau\mathbf{G}$.

To show that the positive-definite property of $\widetilde{D}_{ij}$ is not
automatically satisfied, it is sufficient to consider a specific counter-
example. Taking for instance a simple two-dimensional situation where
$\mathbf{Z}=(Z_{1},Z_{2})$ with an isotropic noise term, i.e.
$\mathbf{B}=B\mathbb{1}$, and with a return-to-equilibrium matrix of the form

$\mathbf{G}=\begin{pmatrix}1&-\kappa\\\ -\kappa&1\end{pmatrix}~{}.$ (83)

The evolution equations for this system are therefore

$\displaystyle\frac{\mathrm{d}Z_{1}}{\mathrm{d}t}$
$\displaystyle=-Z_{1}+\kappa Z_{2}+B\xi_{1}~{},$ (84a)
$\displaystyle\frac{\mathrm{d}Z_{2}}{\mathrm{d}t}$
$\displaystyle=-Z_{2}+\kappa Z_{1}+B\xi_{2}~{},$ (84b)

and from Eq. (80) we obtain that the determinant of the $(2\times 2)$ matrix
$\widetilde{D}_{ij}$ is

$det(\widetilde{\mathbf{D}})=\displaystyle\frac{4B^{2}\left(\tau\langle\xi^{2}\rangle\right)^{2}}{(1+\tau)^{2}-\tau^{2}\kappa^{2}}~{}.$
(85)

As soon as $\tau\neq 0$, this determinant is negative for large values of
$\kappa$ and, hence, such a system is ill-posed. In other words, as soon as
the two components $Z_{1}$ and $Z_{2}$ are strongly coupled, the resulting
probabilistic description becomes ill-based even for such a trivial system.
Note for small values of $\tau$, the inverse matrix can be approximated by
$\widetilde{\mathbf{G}}^{-1}=\mathbb{1}-\tau\mathbf{G}+\mathcal{O}(\tau^{2})$,
which shows that, for the general formulation in Eq. (81), the correlations
$\lambda_{ji}$ can be written as

$\lambda_{ji}=\tau\langle\xi^{2}\rangle B_{ij}-\tau^{2}\langle\xi^{2}\rangle
G_{ik}B_{kj}+\mathcal{O}(\tau^{2})~{}.$ (86)

Using this approximation in Eq. (80), the symmetrical matrix
$\widetilde{D}_{ij}$ is obtained as

$\widetilde{D}_{ij}=2\tau\langle\xi^{2}\rangle(BB^{\bot})_{ij}-\tau^{2}\langle\xi^{2}\rangle
C_{ij}+\mathcal{O}(\tau^{2})~{},$ (87)

where $\mathbf{C}=(BB^{\bot})G^{\bot}+G(BB^{\bot})$ is a symmetrical matrix.
The first term on the rhs of Eq. (87) constitutes a positive-definite matrix
but the second one is unsure. The important point is that this second term is
explicitly dependent upon $\mathbf{G}$. The only possibility to ensure that
the resulting matrix $\widetilde{D}_{ij}$ remains definite positive whatever
the choice of $\mathbf{G}$ is to take the limit $\tau\to 0$. Yet, in order to
retain a non-zero $\widetilde{D}_{ij}$ matrix, this limit must be taken as:
$\tau\to 0$ with $\langle\xi^{2}\rangle\to\infty$, such that
$\lim\limits_{\tau\to 0}\left(\tau\langle\xi^{2}\rangle\right)=K$ where $K$ is
a positive constant. This corresponds to the Markovian approximation which is
now introduced.

####  III.2.2 The Markovian approach

It follows from the preceding discussions that there are two main modeling
philosophies. The first one consists in keeping the same set of variables, for
instance kinetic variables $\mathbf{Z}_{\rm p}^{\rm r}=(\mathbf{X}_{\rm
p},\mathbf{U}_{\rm p})$, whatever the context. We have then to handle colored
noises when ‘external forcing’ involves non-zero time or space correlations
with the risk of ending up with ill-posed PDF formulations, as in the case for
discrete particles in non-fully-resolved turbulent flows. The second modeling
philosophy consists in adjusting the particle state vector by including
additional variables until the eliminated degrees of freedom and/or the
‘external forcing’ can be treated as white-noise terms on the now-extended
mechanical system under consideration. The technical derivations of the
eigenvalues of the diffusion matrix should not hide the physical issues at
stake. In a thermodynamic formulation, the existence of a negative eigenvalue
indicates that one is trying to describe a system whose contact with the
‘external world’ cannot be treated as a contact with a heat bath since it
contains a (negative) correlation and, thus, an underlying order that needs to
be taken into account. On the other hand, positive eigenvalues of the second-
order matrix $\widetilde{D}_{ij}$ means that the corresponding effects can be
regarded as the sum of uncorrelated ‘pure noise’ perturbations, leading to
real diffusive actions on the system. In this second situation, we can regard
the extended system as being in contact with the equivalent of ‘heat bathes’
and we are now dealing with Markovian systems. To underline the importance of
this notion, we quote [86, page 213] who gave an excellent presentation of
this principle: “One needs to keep sufficiently many variables and the
appropriate non-linearities for achieving a realistic description of a system
by Markovian time-evolution equations. We here insist on avoiding explicit
memory effects and take the standpoint that memory effects always indicate the
existence of unrecognized variables relevant to the definition of a proper
system for understanding certain phenomena of interest.”

The application of this principle consists in classifying the degrees of
freedom of a system as slow and fast variables with respect to an observation
time $\Delta t$ which needs to be introduced (in a discrete time version, this
observation time interval corresponds to the time step). Variables whose auto-
correlation timescale are larger than $\Delta t$ are defined as slow variables
while variables with an auto-correlation timescale smaller than $\Delta t$ are
defined as fast variables. As such, this is not sufficient since the
eliminated fast variables would appear as colored noises in the evolution
equations of the retained slow variables. The search is therefore of a scale
separation which allows to treat fast variables as white-noise effects while
slow variables have not changed appreciably (this is the essence of the
‘slaving principle’).

To exemplify this notion, it is instructive to consider a toy model involving
a variable $X_{\rm slow}$ whose time-rate-of-change is $X_{\rm fast}$, i.e.
$\mathrm{d}X_{\rm slow}/\mathrm{d}t=X_{\rm fast}$. In the spirit of the
slaving principle, we consider that $X_{\rm fast}$ is a centered process which
has reached its equilibrium distribution conditioned on a given value of the
slow variable $X_{\rm slow}=x$. Therefore, $X_{\rm fast}$ is a stationary
process with variance $\langle X^{2}_{\rm fast}\rangle$ and an auto-
correlation function $R_{X_{\rm fast}}$, $\langle X^{2}_{\rm fast}\rangle
R_{X_{\rm fast}}(t^{\prime}-t)=\langle X_{\rm fast}(t)X_{\rm
fast}(t^{\prime})\rangle$, which depends only on the time difference
$t^{\prime}-t$ and whose integral timescale is $T_{X_{\rm fast}}$.
Straightforward calculations show that for $X_{\rm slow}(0)=0$

$\frac{\mathrm{d}\langle X^{2}_{\rm slow}\rangle}{\mathrm{d}t}=2\langle
X^{2}_{\rm fast}\rangle\int_{0}^{t}R_{X_{\rm fast}}(s)\,\mathrm{d}s~{}.$ (88)

Since by definition $T_{X_{\rm fast}}=\int_{0}^{+\infty}R_{X_{\rm
fast}}(s)\,\mathrm{d}s$, we get that for ‘long-enough time lapses’

$t\gg T_{X_{\rm fast}}\;\Longrightarrow\;\langle X^{2}_{\rm
slow}\rangle(t)\simeq\left(2\langle X^{2}_{\rm fast}\rangle T_{X_{\rm
fast}}\right)t~{},$ (89)

which is the linear behavior for the second-order moment of $X_{\rm slow}$
characterizing the diffusive regime. If we wish to obtain the same result at
any time $t$ (of the order of the observation time), we need to take the limit
of vanishing timescale $T_{X_{\rm fast}}$ and, to retain a non-zero diffusive
coefficient for the evolution of $\langle X^{2}_{\rm slow}\rangle$ in Eq.
(89), we are led to assume that we have

$T_{X_{\rm fast}}\to 0\;\text{and}\;\langle X^{2}_{\rm
fast}\rangle\to+\infty\quad\text{such that}\quad\langle X^{2}_{\rm
fast}\rangle T_{X_{\rm fast}}\to\mathcal{D}~{},$ (90)

where $\mathcal{D}$ appears as the diffusion coefficient for $X_{\rm slow}$.
This corresponds to the white-noise limit for $X_{\rm fast}$ and, in terms of
the equation for the trajectories of $X_{\rm slow}$, we have replaced an ODE
by a SDE

$\mathrm{d}X_{\rm slow}=X_{\rm
fast}\,\mathrm{d}t\;\Longrightarrow\;\mathrm{d}X_{\rm
slow}=\sqrt{2\mathcal{D}}\,\mathrm{d}W~{},$ (91)

where $W$ is a Wiener process. Note that the resulting diffusion coefficient
$\mathcal{D}$ can be written as

$\mathcal{D}=\int_{0}^{\tau}\langle X_{\rm fast}(t)X_{\rm
fast}(0)\rangle\mathrm{d}t~{},$ (92)

where $\tau$ is an intermediate timescale separating $X_{\rm slow}$ from the
rapidly-varying part of its time derivative, which is here $X_{\rm fast}$ (as
long as $\tau\gg T_{X_{\rm fast}}$, the upper limit of the integral does not
modify significantly the integral compared to $\tau=+\infty$). This is
actually a Green-Kubo expression [86]. Following the slaving principle, an
heuristic formulation consists in keeping the same results conditioned on a
given value $x$ of $X_{\rm slow}$, so that we obtain $\mathcal{D}(x)$ as in
Eq. (92) with the averaging operator being $\langle\cdot\,|\,X_{\rm
slow}=x\rangle$. This turns out to be valid for the diffusion coefficient
resulting from the elimination of fast variables but the shift from an ODE to
a SDE, as in Eq. (91), may involve additional drift terms when $T_{X_{\rm
fast}}(x)$ is an explicit function of the slow variable. In the course of the
following sections, these questions resurface regularly and they are addressed
more rigorously in Sec. V.2. Going back to the well-posed criterion for PDF
equations involving colored noises, it is seen that the conclusion reached to
guarantee the positive nature of the diffusion matrix $D_{ij}$ in Eq. (87) is
the same one with $\tau$, $\langle\xi^{2}\rangle$ and $K$ being replaced here
by $T_{X_{\rm fast}}$, $\langle X^{2}_{\rm fast}\rangle$ and $\mathcal{D}$,
respectively. The conclusion is that well-posed formulations are obtained when
taking the white-noise limit based on a scale separation to distinguish
between slow and (very) fast variables.

For discrete particles in non-fully-resolved turbulent flows, the Markovian
approach leads to include the velocity of the fluid seen in the particle state
vector (having similar timescale as $\mathbf{U}_{\rm p}$ for relatively low-
inertia particles, $\mathbf{U}_{\rm s}$ is clearly a slow variable). However,
do we need additional variables with, for instance, the time derivative of
$\mathbf{U}_{\rm s}$? Answers to this question are provided by the Kolmogorov
theory.

###  III.3 The Kolmogorov theory and Lagrangian models

In spite of some limitations, the Kolmogorov description of turbulent flows
remains the reference theory in turbulence modeling. The first theory was
presented in 1941 (the K41 theory) and later refined to account for
intermittency in 1962 (the refined or K62 theory) [81, 33, 47]. Interestingly,
A. N. Kolmogorov, who was one of the most brilliant mathematician of the 20th
century, chose a rather qualitative approach based on the image of an energy
cascade from which statistical predictions are derived. As in Richardson’s
first pictorial description in 1922 (“Big whorls have little whorls, which
feed on their velocity, and little whorls have lesser whorls, and so on until
viscosity”), the energy cascade corresponds to a description in terms of what
is loosely defined as ‘an eddy’ (or a velocity fluctuating ‘component’),
characterized by its ‘size’ $l$, velocity scale $\delta u_{\rm f}(l)$ and
timescale $\tau_{\rm f}(l)$: energy is produced at the large scales imposed by
the geometry of the flow domain or the boundary conditions and is transferred
through the inertial range (eddies for which energy is neither created nor
dissipated) until it is dissipated at the smallest scales by viscous motions,
see Fig. 9(a). Given its central role, the Kolmogorov theory has made its way
in classical textbooks on turbulence and extensive accounts are available [81,
33, 89]. Consequently, we only give a brief outline of its main
characteristics with a view towards its application for Lagrangian stochastic
models

(a) The energy cascade

(b) Energy spectra as a function of $Re$

Figure 9: The Kolmogorov picture of turbulence: (a) Energy is produced at the
large scales and is transferred through the inertial range before being
dissipated by viscous motions at the same rate $\langle\epsilon_{\rm
f}\rangle$; (b) Evolution of the energy spectra when the Reynolds number based
on the Taylor length scale ${\rm Re}_{\lambda}\sim{\rm Re}^{1/2}$ increases
(reprinted from [89] with permission from Cambridge University Press).

Since traditional presentations tend to concentrate on spatial correlations,
it is worth recalling that the fundamental K41 theory is more general and
proceeds from a Lagrangian vision. In that sense, the most comprehensive
description remains the one given in [81, chapter 8]. It defines the notion of
locally isotropic turbulence by considering the fields relative to a chosen
fluid particle in a small space-time region around that moving particle. For
velocity statistics, this means that we consider the relative field,
$\delta\mathbf{v}_{\rm f}(\tau,\mathbf{r})=\mathbf{U}_{\rm
f}(t_{0}+\tau,\mathbf{x})-\mathbf{U}_{\rm f}(t_{0},\mathbf{x}_{0})$, in the
reference frame moving with the velocity of a chosen fluid particle,
$\mathbf{U}_{\rm f}(t_{0},x_{0})$, and where the space coordinate $\mathbf{r}$
is $\mathbf{r}=\mathbf{x}-\mathbf{x}_{0}-\mathbf{U}_{\rm
f}(t_{0},\mathbf{x}_{0})\,\tau$. Based on this description, the K41 theory
states that, for high-Reynolds-number turbulent flows and for small enough
$r=|\mathbf{r}|$ and $\tau$, turbulence is locally isotropic (in the small
space-time region defined by $\mathbf{r}$ and $\tau$ around the moving fluid
particle). Then, the K41 first similarity hypothesis is that statistics within
that small space-time region are uniquely determined by $\langle\epsilon_{\rm
f}\rangle$ (which represents also the rate of energy transfer in the inertial
range), $\nu_{\rm f}$ (the fluid viscosity) and the space or time coordinates
$\mathbf{r}$ and $\tau$, but not by the specific velocity of the ‘observation
fluid particle ’ $\mathbf{U}_{\rm f}(t_{0},\mathbf{x}_{0})$.

The first outcome of the Kolmogorov theory is the expression of the length,
velocity and time scales of the smallest scales of turbulence

$\eta_{\rm K}=\left(\frac{\nu_{\rm f}^{3}}{\langle\epsilon_{\rm
f}\rangle}\right)^{3/4},\quad u_{\rm K}=\left(\nu_{\rm f}\langle\epsilon_{\rm
f}\rangle\right)^{1/4},\quad\tau_{\rm K}=\left(\frac{\nu_{\rm
f}}{\langle\epsilon_{\rm f}\rangle}\right)^{1/2}~{}.$ (93)

If we introduce $L_{\rm f}$ and $u_{\rm f}$ the length and velocity of the
large-scale motions and use the estimation $\langle\epsilon_{\rm f}\rangle\sim
u_{\rm f}^{3}/L_{\rm f}$, we obtain the ratios between the length and time
scales of the largest to the smallest scales

$\frac{L_{\rm f}}{\eta_{\rm K}}\sim{\rm Re}^{3/4}~{},\quad\frac{T_{\rm
f}}{\tau_{\rm K}}\sim{\rm Re}^{1/2}~{},$ (94)

with $T_{\rm f}=L_{\rm f}/u_{\rm f}$ and ${\rm Re}=u_{\rm f}L_{\rm f}/\nu_{\rm
f}$ the Reynolds number based on the large scales. This provides a way to
assess the complexity involved in turbulence. Indeed, a complete spatial
resolution of the velocity field (i.e., capturing all the turbulent eddies) in
a three-dimensional flow implies that we must handle a number of degrees of
freedom that scales as ${\rm Re}^{9/4}$ which have, furthermore, to be tracked
in time (just to simulate one large-eddy turnover time $T_{\rm f}$ we must
adopt a time resolution of the order of $\tau_{\rm K}$ implying another ${\rm
Re}^{1/2}$ factor in the effort required). Given that many flows corresponds
to Reynolds numbers of the order of $10^{6}-10^{8}$ (for example, in the
atmospheric boundary layer), this estimation demonstrates that we are dealing
with huge numbers of degrees of freedom so that statistical reduced
descriptions are unavoidable.

The Kolmogorov scales allow to properly define the inertial range as the space
and time domain where $\eta_{\rm K}\ll r\ll L_{\rm f}$ and $\tau_{\rm
K}\ll\tau\ll T_{\rm f}$. Then, the second Kolmogorov similarity hypothesis
assumes that, in the inertial range, statistics do not depend anymore on the
fluid viscosity $\nu_{\rm f}$. This yields scaling relations for the eddy
velocity and time scales in the inertial range

$\delta u_{\rm f}(l)=\left(\langle\epsilon_{\rm f}\rangle\,l\right)^{1/3}\sim
u_{\rm f}\left(l/L_{\rm f}\right)^{1/3}~{},\quad\tau_{\rm
f}(l)=\left(l^{2}/\langle\epsilon_{\rm f}\rangle\right)^{1/3}\sim\frac{L_{\rm
f}}{u_{\rm f}}\left(l/L_{\rm f}\right)^{2/3}~{}.$ (95)

##### Eulerian statistics of velocity differences

The most usual application is for the Eulerian fluid velocity correlations
which are characterized by the tensor

$\displaystyle D_{{\rm f},ij}(t,\mathbf{x},\mathbf{r})$
$\displaystyle=\langle\,\delta v_{{\rm f},i}(0,\mathbf{r})\,\delta v_{{\rm
f},j}(0,\mathbf{r})\,\rangle~{},$ (96a)
$\displaystyle=\langle\,\left[\,U_{{\rm f},i}(t,\mathbf{x}+\mathbf{r})-U_{{\rm
f},i}((t,\mathbf{x})\,\right]\,\left[\,U_{{\rm
f},j}(t,\mathbf{x}+\mathbf{r})-U_{{\rm
f},j}((t,\mathbf{x})\,\right]\,\rangle.$ (96b)

When isotropy prevails, this velocity-structure tensor does not depend on
$\mathbf{x}$ anymore and is written as [81, 89]

$D_{{\rm f},ij}(t,\mathbf{r})=D_{\rm f,NN}(t,r)\delta_{ij}+\left[D_{\rm
f,LL}(t,r)-D_{\rm f,NN}(t,r)\right]\frac{r_{i}r_{j}}{r^{2}}~{},$ (97)

where the two scalar functions $D_{\rm f,LL}(t,r)$ and $D_{\rm f,NN}(t,r)$ are
the longitudinal and transverse structure functions, respectively. Moreover,
the continuity constraint implies that $D_{\rm f,NN}$ is uniquely determined
by $D_{\rm f,LL}$ and $D_{{\rm f},ij}(t,\mathbf{r})$ is therefore fully
characterized by one scalar function $D_{\rm f,LL}(t,r)$. In the inertial
range, the second Kolmogorov similarity hypothesis yields then that

$D_{\rm f,LL}(t,r)=C_{2}\left(\langle\epsilon_{\rm
f}\rangle\,r\right)^{2/3}~{},$ (98)

where $C_{2}$ is a constant. This is the same result as the one given in Eqs.
(95) since $D_{\rm f,LL}$ represents $\left(\delta u_{\rm f}(l)\right)^{2}$
(using $\langle\epsilon_{\rm f}\rangle\sim u_{\rm f}^{3}/L_{\rm f}$). The
Fourier transform of $D_{\rm f,LL}(t,r)$ gives the (longitudinal) kinetic
energy spectrum with the well-known $-5/3$ variation in the inertia range (see
Fig. 9(b)). This point has been the subject of numerous experimental and
numerical studies ever since the 1960s and detailed discussions can be found
in the references mentioned above. When $\nu_{\rm f}\to 0$ (with fixed $u_{\rm
f}$ and $L_{\rm f}$), ${\rm Re}\to+\infty$ but the kinetic energy dissipation
rate $\langle\epsilon_{\rm f}\rangle$ tends towards a finite value which is
the only remaining trace of the vanishing viscosity. Smaller and smaller
scales are generated and the energy spectrum is stretched to larger and larger
wave numbers but with the same slope, as shown in Fig. 9(b).

##### Statistics of temporal velocity increments

The Eulerian velocity difference at a fixed point $\mathbf{x}_{0}$ and at two
instants $t_{0}$ and $t_{0}+\tau$ can be expressed in terms of the velocity
field $\delta\mathbf{v}_{\rm f}$ defined in the reference system moving with
the velocity $\mathbf{U}_{\rm f}(t_{0},\mathbf{x}_{0})$ as

$\delta\mathbf{U}^{(t_{0},\mathbf{x}_{0})}_{\rm
f}(\tau)\simeq\delta\mathbf{v}_{\rm f}(\tau,-\mathbf{U}_{\rm
f}(t_{0},\mathbf{x}_{0})\tau).$ (99)

However, the situation is more complicated than for the Eulerian increments
since $\delta\mathbf{U}^{(t_{0},\mathbf{x}_{0})}_{\rm f}(\tau)$ depends
explicitly on the reference velocity $\mathbf{U}_{{\rm f},0}=\mathbf{U}_{\rm
f}(t_{0},\mathbf{x}_{0})$. Therefore, we can only conclude from the Kolmogorov
theory that there is a conditional probability distribution for
$\delta\mathbf{U}^{(t_{0},\mathbf{x}_{0})}_{\rm f}(\tau)$. For a given value
of $\mathbf{U}_{{\rm f},0}$, the tensor $D_{{\rm
f},ij}^{(t_{0},\mathbf{x}_{0})}(\tau)=\langle\delta
U^{(t_{0},\mathbf{x}_{0})}_{{\rm f},i}(\tau)\,\delta
U^{(t_{0},\mathbf{x}_{0})}_{{\rm f},j}(\tau)\rangle$ depends on the two
functions $D^{(t_{0},\mathbf{x}_{0})}_{{\rm f},||}$ and
$D^{(t_{0},\mathbf{x}_{0})}_{{\rm f},\bot}$ which correspond to the
longitudinal and transverse directions, as in Eq. (97). These two functions
depend on $\tau$ as well as on $\mathbf{r}=\mathbf{U}_{{\rm f},0}\tau$ and in
the inertial range we have

$D^{(t_{0},\mathbf{x}_{0})}_{{\rm f},||}=\langle\epsilon_{\rm
f}\rangle\,\tau\,\alpha_{||}\left(\frac{|\mathbf{U}_{{\rm
f},0}|^{2}}{\langle\epsilon_{\rm f}\rangle\tau}\right),\quad
D^{(t_{0},\mathbf{x}_{0})}_{{\rm f},\bot}=\langle\epsilon_{\rm
f}\rangle\,\tau\,\alpha_{\bot}\left(\frac{|\mathbf{U}_{{\rm
f},0}|^{2}}{\langle\epsilon_{\rm f}\rangle\tau}\right),$ (100)

where $\alpha_{||}$ and $\alpha_{\bot}$ are universal functions which have to
be specified. The dependence on $\mathbf{U}_{{\rm f},0}$ and the conditional
nature of the previous results can be removed by resorting to the Taylor, or
frozen turbulence, hypothesis. According to this hypothesis, the turbulent
fluctuations are much smaller than the mean velocity, i.e. $\mathbf{U}_{{\rm
f},0}\simeq\langle\mathbf{U}_{{\rm f},0}\rangle$. Over a small time interval
$\tau$, the turbulent fluctuations at a fixed point are regarded as being
transported without modification at a constant velocity
$\langle\mathbf{U}_{{\rm f},0}\rangle$. This frozen-turbulence assumption
removes the dependence on $\mathbf{U}_{{\rm f},0}$ and allows to write the
statistics of $\delta\mathbf{U}^{(t_{0},\mathbf{x}_{0})}_{\rm f}(\tau)$ in
terms of the velocity-structure functions derived for the Eulerian space
increments by replacing $\mathbf{r}$ by $\langle\mathbf{U}_{{\rm
f},0}\rangle\tau$. This leads to

$D^{(t_{0},\mathbf{x}_{0})}_{{\rm f},||}=D_{||}(\langle\mathbf{U}_{{\rm
f},0}\rangle\tau),\quad D^{(t_{0},\mathbf{x}_{0})}_{{\rm
f},\bot}=D_{\bot}(\langle\mathbf{U}_{{\rm f},0}\rangle\tau),$ (101)

and thus in the inertial range to

$D^{(t_{0},\mathbf{x}_{0})}_{{\rm f},||}=C\left(\langle\epsilon_{\rm
f}\rangle\langle\mathbf{U}_{{\rm f},0}\rangle\,\tau\right)^{2/3},\quad
D^{(t_{0},\mathbf{x}_{0})}_{{\rm
f},\bot}=\frac{4}{3}\,C\left(\langle\epsilon_{\rm
f}\rangle\langle\mathbf{U}_{{\rm f},0}\rangle\,\tau\right)^{2/3}.$ (102)

##### Lagrangian statistics for particle velocity increments

The Kolmogorov theory can be directly applied to the velocity increments of a
fluid particle $\mathrm{d}\mathbf{U}_{\rm f}(\tau)=\mathbf{U}_{\rm
f}(t+\tau)-\mathbf{U}_{\rm f}(t)$ which, by choosing this particle as the
reference one in the Kolmogorov approach, are written as
$\mathrm{d}\mathbf{U}_{\rm f}(\tau)=\delta\mathbf{v}_{\rm f}(\tau,0)$. The
locally isotropic nature of turbulent flows implies that, for $\tau\ll T_{\rm
f}$, the Lagrangian second-order velocity structure function
$\langle\mathrm{d}U_{{\rm f},i}(\tau)\,\mathrm{d}U_{{\rm f},j}(\tau)\rangle$
has an isotropic form

$\langle\mathrm{d}U_{{\rm f},i}(\tau)\,\mathrm{d}U_{{\rm
f},j}(\tau)\rangle=D_{\rm f}^{\rm L}(\tau)\,\delta_{ij}~{}.$ (103)

For time differences in the inertial range, i.e. $\tau_{\rm K}\ll\tau\ll
T_{\rm f}$, the second Kolmogorov hypothesis gives that

$D_{\rm f}^{\rm L}(\tau)=C_{0}\,\langle\epsilon_{\rm f}\rangle\,\tau~{},$
(104)

with $C_{0}$ a constant. This is a significant result since the linear
variation in $\tau$ of the second-order moment of velocity increments is a
signature of white-noise effects and of a diffusive behavior of fluid particle
velocities.

While there is no well-marked separation in terms of length scales, it is
important to note that there is one in terms of timescales. More precisely, it
is seen that $\delta u_{\rm f}(l)$ diminishes as $l$ becomes smaller but
without any sharp decrease between two comparable scales (the decrease is
continuous) whereas there is a clear-cut distinction between the correlation
timescales of a fluid particle velocity $\mathbf{U}_{\rm f}$ and its
acceleration $\mathbf{A}_{\rm f}$. Fluid particle velocities are governed by
the large-scale motions of a turbulent flow and scale as $(\mathbf{U}_{\rm
f})^{2}\sim u_{\rm f}^{2}$ with a timescale $T_{\rm L}\sim T_{\rm f}=u_{\rm
f}^{2}/\langle\epsilon_{\rm f}\rangle$, while fluid particle accelerations are
governed by the small-scale motions and scale as $(\mathbf{A}_{\rm
f})^{2}\sim\langle\epsilon_{\rm f}\rangle/\tau_{\rm K}$ with a timescale
$\tau_{\rm A}$ which is of the order of the Kolmogorov timescale $\tau_{\rm
A}\simeq\tau_{\rm K}$. This situation plays a key role in the analysis of
stochastic models. Indeed, we can use the result obtained with the toy model
presented above, with $X_{\rm slow}$ being the fluid particle velocity $U_{\rm
f}$ and $X_{\rm fast}$ its acceleration $A_{\rm f}$ (in a one-dimensional
notation). In high-Reynolds-number flows where $\tau_{\rm A}\ll T_{\rm L}$ and
with $\langle\epsilon_{\rm f}\rangle$ remaining finite as indicated above, we
retrieve the white-noise limit since we have

$\tau_{\rm A}\to 0\;\text{and}\;\langle A_{\rm
f}^{2}\rangle\to+\infty~{},\quad\text{such that}\quad\tau_{\rm A}\langle
A_{\rm f}^{2}\rangle\to\langle\epsilon_{\rm f}\rangle~{}.$ (105)

This is in line with the prediction already given in Eq. (104) and this
indicates that, in the inertial range, we expect fluid particle velocity
increments to be described by a stochastic model containing a white-noise term
such as $\sqrt{C_{0}\,\langle\epsilon_{\rm f}\rangle}\,\mathrm{d}\mathbf{W}$,
as will be seen in more detail in Sec. IV.

##### Intermittency and refined hypothesis

So far, this outline of the Kolmogorov theory has essentially followed the K41
picture in that the rate of kinetic transfer and dissipation
$\langle\epsilon_{\rm f}\rangle$ is taken as a parameter rather than as a
random variable. This is related to the intermittency of the flow produced by
the inherent transfer processes of turbulence, referred to as ‘inner
intermittency’ (to be distinguished from external intermittency due, for
example, to large-scale mixing between laminar and turbulent flows) and we
have assumed that $\epsilon_{\rm f}\simeq\langle\epsilon_{\rm f}\rangle$.
While the above picture can be kept substituting $\langle\epsilon_{\rm
f}\rangle$ with $\epsilon_{\rm f}$ , especially for second-order moments where
the impact of intermittency is small, a refined description based on a log-
normal distribution for the dissipation of kinetic energy was developed in the
K62 theory. For this question, the basic presentation of the refined
similarity hypothesis is [81, section 25.2] and this question has been
investigated in several studies, in particular in [33]. Without addressing in
detail this issue, related comments are proposed in connection with coherent
structures in Sec. IV.4.3.

##### Particle state vectors and Lagrangian models

From these results, it appears that the Kolmogorov theory provides clear
indications as to which state vectors should be selected to follow the
Markovian approach. For tracer particles, the relevant state vector is
$\mathbf{Z}_{\rm f}=(\mathbf{X}_{\rm f},\mathbf{U}_{\rm f})$ while for
discrete particles the velocity of the fluid seen should be included, leading
to $\mathbf{Z}_{\rm p}=(\mathbf{X}_{\rm p},\mathbf{U}_{\rm p},\mathbf{U}_{\rm
s})$, since we can expect the ‘acceleration’ of the fluid seen (or its time-
rate-of-change) to be also a fast variable. In terms of the trajectories of
the stochastic process $\mathbf{Z}_{\rm p}$, the general structure of one-
particle PDF models for discrete particles in turbulent flows is

$\displaystyle\mathrm{d}\mathbf{X}_{\rm p}$ $\displaystyle=\mathbf{U}_{\rm
p}\,\mathrm{d}t~{},$ (106a) $\displaystyle\mathrm{d}\mathbf{U}_{\rm p}$
$\displaystyle=\frac{\mathbf{U}_{\rm s}-\mathbf{U}_{\rm p}}{\tau_{\rm
p}}\,\mathrm{d}t+\mathbf{F}_{\rm f\to p}\,\mathrm{d}t~{},$ (106b)
$\displaystyle\mathrm{d}\mathbf{U}_{\rm s}$ $\displaystyle=(\text{stochastic
model})~{},$ (106c)

where $\mathbf{F}_{\rm f\to p}$ in Eq. (106b) stands for additional forces due
to the fluid acting on particles (i.e. other than drag). The task ahead is to
work out appropriate stochastic models for $\mathbf{U}_{\rm s}$ in Eq. (106c).

##  IV Modeling the velocity of the fluid seen

###  IV.1 The physics of particle dispersion

####  IV.1.1 The crossing-trajectory effect

In a nutshell, particle dispersion in turbulent flows is due to the
fluctuating fluid velocities encountered by discrete particles and is, thus,
governed by the statistical properties of the velocity of the fluid seen
$\mathbf{U}_{\rm s}(t)$. Our purpose, in this section, is to bring out the
main physical phenomena at play and the relevant variables that need to be
accounted for when modeling $\mathbf{U}_{\rm s}(t)$.

The problem can be described as follows (see also Fig. 10). At time $t$, we
consider a discrete particle (P) located at a position $\mathbf{X}_{\rm
p}(t)$, with a velocity $\mathbf{U}_{\rm p}(t)$ and a velocity of the fluid
seen $\mathbf{U}_{\rm s}(t)$. Then, after a small time interval $\Delta t$,
the discrete particle has a probability to move to a downstream location
$\mathbf{X}_{\rm p}(t+\Delta t)$, while a fluid particle (F) starting from
$\mathbf{X}_{\rm p}(t)$ with the velocity $\mathbf{U}_{\rm s}(t)$ at time $t$
has a probability to move to another position $\mathbf{X}_{\rm f}(t+\Delta
t)$. The modeling issue is to estimate $\mathbf{U}_{\rm s}(t+\Delta t)$, which
is the velocity of another fluid particle than (F), say (F’), at time
$t+\Delta t$. In most cases, due to particle inertia and/or external forces
inducing mean velocity slips $\langle\mathbf{U}_{\rm r}\rangle$ (with
$\mathbf{U}_{\rm r}=\mathbf{U}_{\rm s}-\mathbf{U}_{\rm p}$ the relative
velocity), the trajectories of the discrete particle (P) and of the fluid
particle (F) located at the same point at time $t$ separate during the time
interval $\Delta t$. This is referred to as the crossing-trajectory effect
(CTE). As such, the CTE is related to particle inertia and the various effects
mentioned above but its precise definition in the present modeling context
will be narrowed down below. When we re-express the velocity of the fluid seen
in terms of the instantaneous fluid velocity field $\mathbf{U}_{\rm
f}(t,\mathbf{x})$, it appears that the modeling issue is to simulate
$\mathbf{U}_{\rm s}(t+\Delta t)=\mathbf{U}_{\rm f}(t+\Delta t,\mathbf{X}_{\rm
p}(t+\Delta t))$ knowing $\mathbf{U}_{\rm s}(t)=\mathbf{U}_{\rm
f}(t,\mathbf{X}_{\rm p}(t))$. This clearly involves a two-time two-point
conditional fluid velocity correlation, which reveals that it can only be
treated without approximation at the level of two-particle PDF models, i.e.
when following pairs of fluid particles or when generating spatial information
about the fluid velocity field in the vicinity of every discrete particle
being tracked. At the level of one-particle PDF models, this constitutes an
open issue and additional assumptions have to be made.

Figure 10: Illustration of the CTE and of its modeled representation. The CTE
occurs as soon as there is a relative velocity between discrete particles and
the fluid, due to particle inertia and mean velocity drifts (left). In the
present formulations, it is accounted for only due to mean relative velocities
(right).