#  Liouville models of particle-laden flow

Daniel Domínguez-Vázquez, Gustaaf B. Jacobs
Department of Aerospace Engineering
San Diego State University
San Diego, CA 92182, USA
&Daniel M. Tartakovsky
Department of Energy Science and Engineering
Stanford University
Stanford, CA 94305, USA
<EMAIL_ADDRESS>to whom correspondence should be addressed:
gjacobs<EMAIL_ADDRESS>
###### Abstract

Langevin (stochastic differential) equations are routinely used to describe
particle-laden flows. They predict Gaussian probability density functions
(PDFs) of a particle’s trajectory and velocity, even though experimentally
observed dynamics might be highly non-Gaussian. Our Liouville approach
overcomes this dichotomy by replacing the Wiener process in the Langevin
models with a (small) set of random variables, whose distributions are tuned
to match the observed statistics. This strategy gives rise to an exact
(deterministic, first-order, hyperbolic) Liouville equation that describes the
evolution of a joint PDF in the augmented phase-space spanned by the random
variables and the particle position and velocity. Analytical PDF solutions for
canonical models of particle-laden flows serve to establish the relationship
between the Langevin and Liouville approaches. Finally, our framework is used
to derive a new analytical PDF model for fluidized homogeneous heating
systems.

##  1 Introduction

Dispersion of particles in chaotic flows is of fundamental importance to a
plethora of natural environments and engineering technologies. For example,
mixing of evaporating fuel droplets in a combustor determines the fuel-to-air
ratio, homogeneity of the mixture, and ultimately the efficiency of the
combustor; upwelling of surface dust particles by shocked, unsteady nozzle-jet
flows can damage launch vehicles; etc. An accurate model prediction of such
flows at large scales is necessary to elucidate fundamental physics and
improve engineering design.

A typical model of turbulent particle-laden flows involves the spatial scales
ranging from tens of meters (the application scale) to microns (the smallest
scales of turbulence and particulates). Consequently, a predictive simulation
based on “first principles” is not practical and a portion of the smaller
scales has to be treated empirically through stochastic “subgrid-scale”
models, as it is done in, e.g., large-eddy simulations (LES) or Reynolds-
averaged Navier-Stokes (RANS) models. In the context of particle-laden flows,
this stochasticity stems from either turbulence of a carrier gas/liquid
(subgrid velocity fluctuations) or from random particle kinetics.

Classical kinetic theory [1, 2, 3] deals with stochasticity on the particle
phase through an equation for the joint probability density function (PDF) in
the six–dimensional phase space composed by particle position and velocity.
This equation requires a closure approximation, which can be derived either
analytically or numerically [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
18]. The high–dimensionality of this formulation leads to the so–called curse
of dimensionality, which has inspired reduced formulations based on a number
of moments of the joint PDF [19, 20, 21, 22, 23, 24, 25, 26].

Equivalent stochastic differential equations (SDEs), which in this field are
referred to as a generalized Langevin model [18], can be used to reduce the
dimensionality of this PDF equation. The model relies on Wiener increments and
random walks to account for stochasticity in particle trajectories [27, 28];
more recent versions of such stochastic models of the dispersed particle phase
have been proposed in the literature [29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
10, 39, 40, 41, 42, 43, 44, 45]. A representative formulation of this Langevin
approach, in $d$ spatial dimensions, describes the temporal evolution of an
inertial point particle’s location,
$\boldsymbol{X}_{\text{p}}(t)\in\mathbb{R}^{d}$, and velocity,
$\boldsymbol{U}_{\text{p}}(t)\in\mathbb{R}^{d}$, by adding Wiener increments
to the particle position and velocity to Newton’s second law:

$\displaystyle\text{d}{\boldsymbol{X}_{\text{p}}}$
$\displaystyle={\boldsymbol{U}_{\text{p}}}\text{d}t+\boldsymbol{b}_{x}(\boldsymbol{X}_{\text{p}},\boldsymbol{U}_{\text{p}},t)\,\text{d}\boldsymbol{W}_{x},$
(1a) $\displaystyle\text{d}{\boldsymbol{U}_{\text{p}}}$
$\displaystyle=\boldsymbol{\mu}_{a}(\boldsymbol{X}_{\text{p}},\boldsymbol{U}_{\text{p}},t)\text{d}t+\boldsymbol{b}_{u}(\boldsymbol{X}_{\text{p}},\boldsymbol{U}_{\text{p}},t)\,\text{d}\boldsymbol{W}_{u}.$
(1b)

The particle’s slow-varying average acceleration,
$\boldsymbol{\mu}_{a}\in\mathbb{R}^{d}$, might depend in general on
$\boldsymbol{U}_{\text{p}}$ and $\boldsymbol{X}_{\text{p}}$ through the
evaluation of the bulk-flow velocity at the particle position,
$\boldsymbol{u}(\boldsymbol{X}_{\text{p}},t)\in\mathbb{R}^{d}$ and the
relative velocity of the two phases, time, and other model parameters. The
particle’s fast-varying dynamics is represented by the (random) Wiener process
$\boldsymbol{W}_{u}(t)\in\mathbb{R}^{m}$, whose strength is encapsulated in
the (deterministic) diffusion tensor $\boldsymbol{b}_{u}\in\mathbb{R}^{d\times
m}$ that, in general, varies with $\boldsymbol{X}_{\text{p}}$,
$\boldsymbol{U}_{\text{p}}$, and time $t$. Similarly, Wiener increments
$\boldsymbol{W}_{x}(t)\in\mathbb{R}^{m}$ are usually added to the particle
position with diffusion tensor $\boldsymbol{b}_{x}\in\mathbb{R}^{d\times m}$
(see for example Refs.[41] and [46]). The joint PDF of the particle position
and velocity,
$f_{\boldsymbol{X}\boldsymbol{U}}(\boldsymbol{x}_{\text{p}},\boldsymbol{u}_{\text{p}};t)$,
in the case of Gaussian white-noise in (1), satisfies the $(2d)$-dimensional
Fokker-Planck equation

$\displaystyle\begin{split}\frac{\partial
f_{\boldsymbol{X}\boldsymbol{U}}}{\partial
t}&+\nabla_{\boldsymbol{x}_{\text{p}}}\cdot\left(\boldsymbol{u}_{\text{p}}f_{\boldsymbol{X}\boldsymbol{U}}\right)+\nabla_{\boldsymbol{u}_{\text{p}}}\cdot\left(\boldsymbol{\mu}_{a}(\boldsymbol{x}_{\text{p}},\boldsymbol{u}_{\text{p}},t)f_{\boldsymbol{X}\boldsymbol{U}}\right)=\nabla\cdot(\boldsymbol{\mathcal{D}}\nabla
f_{\boldsymbol{X}\boldsymbol{U}}),\end{split}$ (2)

with the diffusion tensor
$\boldsymbol{\mathcal{D}}=\boldsymbol{b}\boldsymbol{b}^{\top}/2\in\mathbb{R}^{2d\times
2d}$ and
$\nabla=(\nabla_{\boldsymbol{x}_{\text{p}}},\nabla_{\boldsymbol{u}_{\text{p}}})^{\top}\in\mathbb{R}^{2d}$.

At least two major challenges hinder the development of consistent Langevin
models for particle-laden flows. First, since the Wiener process is Gaussian
and white noise, it gives rise to the Gaussian PDF
$f_{\boldsymbol{X}\boldsymbol{U}}(\boldsymbol{x}_{\text{p}},\boldsymbol{u}_{\text{p}};t)$.
Yet, empirically observed (and numerically simulated) single-point statistics
of particles are rarely Gaussian [47, 48, 49, 50, 51], which suggests that
models like (1) provide ad hoc (and possibly inaccurate) representations of
subgrid particle dynamics. Second, closure models for diffusion coefficients
in turbulent flows, let alone particle-laden turbulent flows, remain one of
the open problems in fluid mechanics. In addition to having these foundational
issues, the Langevin approach is also computationally expensive when Monte
Carlo methods are used to estimate the tails of joint PDFs.

A way to overcome these limitations is to replace the Wiener increments in (1)
with random forcings, which are represented via series expansions involving
proper sets of random variables [52, 53, 51, 54, 55]. This approach replaces
(stochastic) Langevin equations (1) with a coupled system of $2d$ ordinary
differential equations with random coefficients,

$\displaystyle\frac{\text{d}{\boldsymbol{X}_{\text{p}}}}{\text{d}t}$
$\displaystyle={\boldsymbol{U}}_{\text{p}}+\sum_{i=1}^{N}\boldsymbol{Z}_{i}\circ\boldsymbol{\varphi}_{i}(\boldsymbol{X}_{\text{p}},\boldsymbol{U}_{\text{p}},t),$
(3a) $\displaystyle\frac{\text{d}{\boldsymbol{U}_{\text{p}}}}{\text{d}t}$
$\displaystyle=\boldsymbol{\mu}_{a}(\boldsymbol{X}_{\text{p}},\boldsymbol{U}_{\text{p}},t)+\sum_{i=1}^{N}\boldsymbol{\Xi}_{i}\circ\boldsymbol{\phi}_{i}(\boldsymbol{X}_{\text{p}},\boldsymbol{U}_{\text{p}},t).$
(3b)

Here, $\circ$ is the Hadamard product, $\boldsymbol{Z}_{i}\in\mathbb{R}^{d}$
and $\boldsymbol{\Xi}_{i}\in\mathbb{R}^{d}$ are the vectors of random
coefficients in the $N$th-order polynomials with orthogonal bases
$\boldsymbol{\varphi}_{i}\in\mathbb{R}^{d}$ and
$\boldsymbol{\phi}_{i}\in\mathbb{R}^{d}$ respectively, defined
deterministically. While the Wiener process, $\boldsymbol{W}_{x}(t)$ or
$\boldsymbol{W}_{u}(t)$ in (1), is a random _function_ , an infinite-
dimensional object that necessitates the use of stochastic calculus, the
parameters $\boldsymbol{Z}_{i}$ and $\boldsymbol{\Xi}_{i}$ are random
_variables_ ; hence, for any realization $\boldsymbol{z}_{i}$ and
$\boldsymbol{\xi}_{i}$ of $\boldsymbol{Z}_{i}$ and $\boldsymbol{\Xi}_{i}$, (3)
are deterministic differential equations amenable to a standard numerical or
analytical treatment. Moreover, the joint PDF of the random inputs and outputs
in (3) defined in an augmented phase space that includes the random
coefficients,
$f_{\boldsymbol{X}\boldsymbol{U}\boldsymbol{\Xi}\boldsymbol{Z}}(\boldsymbol{x}_{\text{p}},\boldsymbol{u}_{\text{p}},\boldsymbol{\xi},\boldsymbol{z};t)$,
satisfies the Liouville equation

$\displaystyle\begin{split}&\frac{\partial
f_{\boldsymbol{X}\boldsymbol{U}\boldsymbol{\Xi}\boldsymbol{Z}}}{\partial
t}+\nabla_{\boldsymbol{x}_{\text{p}}}\cdot[(\boldsymbol{u}_{\text{p}}+\sum_{i=1}^{N}\boldsymbol{z}_{i}\circ\boldsymbol{\varphi}_{i}(\boldsymbol{x}_{\text{p}},\boldsymbol{u}_{\text{p}},t))f_{\boldsymbol{X}\boldsymbol{U}\boldsymbol{\Xi}\boldsymbol{Z}}]\\\
&+\nabla_{\boldsymbol{u}_{\text{p}}}\cdot[(\boldsymbol{\mu}_{a}(\boldsymbol{x}_{\text{p}},\boldsymbol{u}_{\text{p}},t)+\sum_{i=1}^{N}\boldsymbol{\xi}_{i}\circ\boldsymbol{\phi}_{i}(\boldsymbol{x}_{\text{p}},\boldsymbol{u}_{\text{p}},t))f_{\boldsymbol{X}\boldsymbol{U}\boldsymbol{\Xi}\boldsymbol{Z}}]=0,\end{split}$
(4)

a hyperbolic partial-differential equation that admits a solution via the
method of characteristics (MoC) [55]. This solution might or might not be
Gaussian, depending on, e.g., the PDF of the initial data or the random
coefficients. Alternatively, (3)–(4) can be used to learn the series in (3)
from data, making the Liouville approach adequate for developing stochastic
models using data assimilation techniques [56]. We point out that the
functions $\boldsymbol{\varphi}_{i}$ and $\boldsymbol{\phi}_{i}$ do not need
to act as basis functions necessarily for the formulation to be valid.

Our analysis pursues two goals. The first is to elucidate that the Liouville
model of particle-laden flows, i.e., Eq. (4), amounts to a generalization of
the Langevin approach. We show it in two canonical examples in which the
former reduces to the latter when the random variables $\boldsymbol{Z}_{i}$
and $\boldsymbol{\Xi}_{i}$ are Gaussian and the source terms are found to fit
the solution of the Fokker-Planck equivalent model. An analytical
generalization for non-Gaussian systems is provided with the Liouville model.
The second goal is to show the potential of the Liouville approach for
developing stochastic models of the particle phase. By adjusting the random
coefficients to fit data, the Liouville approach can be framed in a data-
driven manner. We propose an analytical model for fluidized homogeneous
heating systems (FHHS) based on the Liouville approach. The accuracy of this
model is established via comparison with the predictions of both particle-
resolved direct numerical simulations (PR-DNS) and previous models proposed in
the literature [57, 39, 41, 58, 43].

Section 2 contains a brief formulation of both Langevin and Liouville
frameworks for modeling the particle dynamics in multiphase flows. In Sections
3, we derive analytical solutions of the Fokker-Planck equation and the
Liouville equation for canonical-flow examples and establish the equivalency
between the Langevin and Liouville approaches. In Section 4, our Liouville
approach is used to derive a new model of homogeneous fluidized heated
systems. Main conclusions drawn from this study are collated in Section 5.

##  2 Stochastic Models of Particle Dynamics

If (the sum of) forces
$\boldsymbol{g}(\boldsymbol{x}_{\text{p}},\boldsymbol{u}_{\text{p}},t)$ acting
on a particle were known with certainty, then its trajectory
$\boldsymbol{x}_{\text{p}}(t)$ and velocity $\boldsymbol{u}_{\text{p}}(t)$
would be deterministic and computable from Newton’s second law,

$\displaystyle\frac{\text{d}\boldsymbol{x}_{\text{p}}}{\text{d}t}=\boldsymbol{u}_{\text{p}},\qquad\frac{\text{d}\boldsymbol{u}_{\text{p}}}{\text{d}t}=\boldsymbol{g}.$
(5)

In particle-laden flows, these forces are unresolved/uncertain, e.g., due to
multi-particle interactions or drag on particles of non-canonical shapes. When
treated probabilistically, this uncertainty gives rise to stochastic dynamics,
such that the deterministic functions $\boldsymbol{x}_{\text{p}}(t)$ and
velocity $\boldsymbol{u}_{\text{p}}(t)$ are replaced with the random processes
$\boldsymbol{X}_{\text{p}}(t)$ and velocity $\boldsymbol{U}_{\text{p}}(t)$.

###  2.1 Particle trajectory

A reduced-complexity model of particle dynamics superimposes stochastic
fluctuations onto a given deterministic drift velocity
$\boldsymbol{u}_{\text{p}}(t)$. This model is used for example in the
literature to compute Lagrangian Coherent Structures (LCS) of fluid particles
for uncertain fluid velocity field [59], also denoted particle Langevin model
[41] for inertial particles with unresolved drag.

####  2.1.1 Langevin approach

For deterministic drift velocity $\boldsymbol{u}_{\text{p}}(t)$ and dispersion
tensor $\boldsymbol{b}_{x}$, the Langevin approach is to replace (5) with
stochastic ODEs [59, 41]

$\displaystyle\text{d}{\boldsymbol{X}_{\text{p}}}=\boldsymbol{u}_{\text{p}}\text{d}t+\boldsymbol{b}_{x}({\boldsymbol{X}_{\text{p}}},t)\text{d}\boldsymbol{W}_{x}.$
(6)

A single-point PDF of the stochastic trajectory $\boldsymbol{X}_{\text{p}}(t)$
in the Langevin equation (6),
$f_{\boldsymbol{X}}(\boldsymbol{x}_{\text{p}};t)$, satisfies the Fokker-Planck
equation [60],

$\displaystyle\frac{\partial f_{\boldsymbol{X}}}{\partial
t}+\nabla_{\boldsymbol{x}_{\text{p}}}\cdot\left(\boldsymbol{u}_{\text{p}}f_{\boldsymbol{X}}\right)=\nabla_{\boldsymbol{x}_{\text{p}}}\cdot(\boldsymbol{\mathcal{D}}\nabla_{\boldsymbol{x}_{\text{p}}}f_{\boldsymbol{X}}),\quad\boldsymbol{\mathcal{D}}=\frac{\boldsymbol{b}_{x}\,\boldsymbol{b}_{x}^{\top}}{2}.$
(7)

####  2.1.2 Liouville approach

Instead of using the Wiener process $\boldsymbol{W}_{x}(t)$, the Liouville
approach accounts for stochasticity of the particle dynamics using time-
invariant random variables $\boldsymbol{Z}_{i}$, characterized by the joint
PDF $f_{\boldsymbol{Z}}(\boldsymbol{z})$. Specifically, the stochastic ODE (6)
is replaced with an ODE with random coefficients,

$\displaystyle\frac{\text{d}\boldsymbol{X}_{\text{p}}}{\text{d}t}=\boldsymbol{u}_{\text{p}}+\sum_{i=1}^{N}\boldsymbol{Z}_{i}\circ\boldsymbol{\varphi}_{i}(\boldsymbol{X}_{\text{p}},t),$
(8)

where $\boldsymbol{\varphi}_{i}(\boldsymbol{X}_{\text{p}},t)$ is a prescribed
vector function conforming a basis, for example Chebyshev polynomials. While
the derivation of the PDF equation for
$f_{\boldsymbol{X}}(\boldsymbol{x}_{\text{p}};t)$ requires a closure
approximation [61], the equation for the joint PDF
$f_{\boldsymbol{X}\boldsymbol{Z}}(\boldsymbol{x}_{\text{p}},\boldsymbol{z};t)$
of the model input $\boldsymbol{Z}_{i}$ and output
$\boldsymbol{X}_{\text{p}}(t)$ is exact:

$\displaystyle\frac{\partial f_{\boldsymbol{X}\boldsymbol{Z}}}{\partial
t}+\nabla_{{\boldsymbol{x}}_{\text{p}}}\cdot\left[\left(\boldsymbol{u}_{\text{p}}+\sum_{i=1}^{N}\boldsymbol{z}_{i}\circ\boldsymbol{\varphi}_{i}(\boldsymbol{x}_{\text{p}},t)\right)f_{\boldsymbol{X}\boldsymbol{Z}}\right]=0.$
(9)

This result can be derived via either the Liouville theorem [62, 63] or the
method of distributions [61, 51, 55].

Since $f_{\boldsymbol{X}}(\boldsymbol{x}_{\text{p}};t)$ is the marginal of
$f_{\boldsymbol{X}\boldsymbol{Z}}(\boldsymbol{x}_{\text{p}},\boldsymbol{z};t)$,
i.e.,

$f_{\boldsymbol{X}}(\boldsymbol{x}_{\text{p}};t)=\int
f_{\boldsymbol{X}\boldsymbol{Z}}(\boldsymbol{x}_{\text{p}},\boldsymbol{z};t)\text{d}\boldsymbol{z},$

the integration of (9) yields

$\displaystyle\frac{\partial f_{\boldsymbol{X}}}{\partial
t}+\nabla_{\boldsymbol{x}_{\text{p}}}\cdot\left(\boldsymbol{u}_{\text{p}}f_{\boldsymbol{X}}\right)=-\nabla_{\boldsymbol{x}_{\text{p}}}\cdot\left(\int\sum_{i=1}^{N}\boldsymbol{z}_{i}\circ\boldsymbol{\varphi}_{i}f_{\boldsymbol{X}\boldsymbol{Z}}\,\text{d}\boldsymbol{z}\right).$
(10)

As mentioned above, the derivation of a PDF equation for $f_{\boldsymbol{X}}$
requires a closure to approximate the right hand side term. Comparison of (10)
and (7) indicates that the Langevin and Liouville approaches yield the
identical predictions of the PDF
$f_{\boldsymbol{X}}(\boldsymbol{x}_{\text{p}};t)$ if

$\displaystyle\nabla_{\boldsymbol{x}_{\text{p}}}\cdot\left(\boldsymbol{\mathcal{D}}\nabla_{\boldsymbol{x}_{\text{p}}}f_{\boldsymbol{X}}+\int\sum_{i=1}^{N}\boldsymbol{z}_{i}\circ\boldsymbol{\varphi}_{i}f_{\boldsymbol{X}\boldsymbol{Z}}\,\text{d}\boldsymbol{z}\right)=0.$
(11)

If a solution of a Fokker-Planck equation is known, the above identity may be
used to develop an equivalent Liouville model. Taking moments of (7) and (9),
similar relations can be found for each moment, adding constrains to the
Liouville approach to match an equivalent Langevin model. A moment model based
on the Liouville approach is provided in Appendix A.

###  2.2 Particle velocity

A class of models of particle-laden flows takes $\mathbf{g}$ in (5) to
represent the drag force acting on a particle (see Ref. [64] and references
therein). This force is proportional to the difference between the carrier
flow velocity at the point $\boldsymbol{x}_{\text{p}}$ occupied by the
particle, $\boldsymbol{u}(\boldsymbol{x}_{\text{p}},t)$, and the particle
velocity, $\boldsymbol{u}_{\text{p}}(t)$, so that (5) is replaced with

$\displaystyle\frac{\text{d}\boldsymbol{x}_{\text{p}}}{\text{d}t}=\boldsymbol{u}_{\text{p}},\qquad\frac{\text{d}\boldsymbol{u}_{\text{p}}}{\text{d}t}=\frac{\mathcal{F}}{\tau_{\text{p}}}\left(\boldsymbol{u}-\boldsymbol{u}_{\text{p}}\right).$
(12)

The coefficient of proportionality comprises the correction factor
$\mathcal{F}$, which accounts for a possible deviation from the Stokes drag
(due to, e.g., finite particle-Reynolds number, particle Mach number, volume
fraction of the particles), and the particle’s response time
$\tau_{\text{p}}$.

####  2.2.1 Langevin approach

A probabilistic treatment of (12), which is referred to velocity Langevin
model [41], uses Wiener increments to account for velocity fluctuations:

$\displaystyle\text{d}{\boldsymbol{X}_{\text{p}}}$
$\displaystyle={\boldsymbol{U}_{\text{p}}}\text{d}t,$ (13a)
$\displaystyle\text{d}{\boldsymbol{U}_{\text{p}}}$
$\displaystyle=\frac{\mathcal{F}}{\tau_{\text{p}}}\left[\boldsymbol{u}({\boldsymbol{X}_{\text{p}}},t)-{\boldsymbol{U}_{\text{p}}}\right]\text{d}t+\boldsymbol{b}_{u}({\boldsymbol{X}_{\text{p}}},{\boldsymbol{U}_{\text{p}}},t)\,\text{d}\boldsymbol{W}_{u}.$
(13b)

The Fokker-Planck equation for the joint PDF
$f_{\boldsymbol{X}\boldsymbol{U}}\left({\boldsymbol{x}_{\text{p}}},{\boldsymbol{u}_{\text{p}}};t\right)$
is given by

$\displaystyle\begin{split}\frac{\partial
f_{\boldsymbol{X}\boldsymbol{U}}}{\partial
t}&+\nabla_{\boldsymbol{x}_{\text{p}}}\cdot(\boldsymbol{u}_{\text{p}}f_{\boldsymbol{X}\boldsymbol{U}})+\nabla_{\boldsymbol{u}_{\text{p}}}\cdot\left[\frac{\mathcal{F}}{\tau_{\text{p}}}(\boldsymbol{u}-\boldsymbol{u}_{\text{p}})f_{\boldsymbol{X}\boldsymbol{U}}\right]=\nabla_{\boldsymbol{u}_{\text{p}}}\cdot(\boldsymbol{\mathcal{D}}\nabla_{\boldsymbol{u}_{\text{p}}}f_{\boldsymbol{X}\boldsymbol{U}}),\end{split}$
(14)

with prescribed functions
$\boldsymbol{u}=\boldsymbol{u}(\boldsymbol{x}_{\text{p}},t)$ and
$\boldsymbol{\mathcal{D}}(\boldsymbol{x}_{\text{p}},\boldsymbol{u}_{\text{p}},t)=\boldsymbol{b}_{u}\boldsymbol{b}_{u}^{\top}/2$.

####  2.2.2 Liouville approach

Our Liouville model relies on random variables $\boldsymbol{Z}_{i}$ and
$\boldsymbol{\Xi}_{i}$, with PDFs $f_{\boldsymbol{Z}}(\boldsymbol{z})$ and
$f_{\boldsymbol{\Xi}}(\boldsymbol{\xi})$, and basis functions given by
$\boldsymbol{\phi}_{i}(\boldsymbol{X}_{\text{p}},\boldsymbol{U}_{\text{p}},t)$
and
$\boldsymbol{\varphi}_{i}(\boldsymbol{X}_{\text{p}},\boldsymbol{U}_{\text{p}},t)$
to represent the stochastic particle dynamics:

$\displaystyle\frac{\text{d}\boldsymbol{X}_{\text{p}}}{\text{d}t}$
$\displaystyle=\boldsymbol{U}_{\text{p}}+\sum_{i=1}^{N}\boldsymbol{Z}_{i}\circ\boldsymbol{\varphi}_{i}(\boldsymbol{X}_{\text{p}},\boldsymbol{U}_{\text{p}},t),$
(15a) $\displaystyle\frac{\text{d}\boldsymbol{U}_{\text{p}}}{\text{d}t}$
$\displaystyle=\frac{\mathcal{F}}{\tau_{\text{p}}}\left[\boldsymbol{u}({\boldsymbol{X}_{\text{p}}},t)-{\boldsymbol{U}_{\text{p}}}\right]+\sum_{i=1}^{N}\boldsymbol{\Xi}_{i}\circ\boldsymbol{\phi}_{i}(\boldsymbol{X}_{\text{p}},\boldsymbol{U}_{\text{p}},t).$
(15b)

The corresponding Liouville equation for the joint PDF
$f_{\boldsymbol{X}\boldsymbol{U}\boldsymbol{Z}\boldsymbol{\Xi}}(\boldsymbol{x}_{\text{p}},\boldsymbol{u}_{\text{p}},\boldsymbol{z},\boldsymbol{\xi};t)$
is

$\displaystyle\begin{split}\frac{\partial
f_{\boldsymbol{X}\boldsymbol{U}\boldsymbol{Z}\boldsymbol{\Xi}}}{\partial
t}&+\nabla_{\boldsymbol{x}_{\text{p}}}\cdot\left[\left(\boldsymbol{u}_{\text{p}}+\sum_{i=1}^{N}\boldsymbol{z}_{i}\circ\boldsymbol{\varphi}_{i}(\boldsymbol{x}_{\text{p}},\boldsymbol{u}_{\text{p}},t)\right)f_{\boldsymbol{X}\boldsymbol{U}\boldsymbol{Z}\boldsymbol{\Xi}}\right]\\\
&+\nabla_{\boldsymbol{u}_{\text{p}}}\cdot\left[\left(\frac{\mathcal{F}}{\tau_{\text{p}}}\left(\boldsymbol{u}(\boldsymbol{x}_{\text{p}},t)-\boldsymbol{u}_{\text{p}}\right)+\sum_{i=1}^{N}\boldsymbol{\xi}_{i}\circ\boldsymbol{\phi}_{i}(\boldsymbol{x}_{\text{p}},\boldsymbol{u}_{\text{p}},t)\right)f_{\boldsymbol{X}\boldsymbol{U}\boldsymbol{Z}\boldsymbol{\Xi}}\right]=0.\end{split}$
(16)

Since $f_{\boldsymbol{X}\boldsymbol{U}}$ is the marginal of
$f_{\boldsymbol{X}\boldsymbol{U}\boldsymbol{Z}\boldsymbol{\Xi}}$, the
integration of (16) over the coefficients yields an unclosed PDF equation for
$f_{\boldsymbol{X}\boldsymbol{U}}$. Comparing that equation with the Fokker-
Planck equation (14) yields a compatibility condition,

$\displaystyle\begin{split}&\nabla_{\boldsymbol{x}_{\text{p}}}\cdot\left[\int\sum_{i=1}^{N}\boldsymbol{z}_{i}\circ\boldsymbol{\varphi}_{i}\left(\int
f_{\boldsymbol{X}\boldsymbol{U}\boldsymbol{Z}\boldsymbol{\Xi}}\,\text{d}\boldsymbol{\xi}\right)\text{d}\boldsymbol{z}\right]+\nabla_{\boldsymbol{u}_{\text{p}}}\cdot\left[\boldsymbol{\mathcal{D}}\nabla_{\boldsymbol{u}_{\text{p}}}f_{\boldsymbol{X}\boldsymbol{U}}+\int\sum_{i=1}^{N}\boldsymbol{\xi}_{i}\circ\boldsymbol{\phi}_{i}\left(\int
f_{\boldsymbol{X}\boldsymbol{U}\boldsymbol{Z}\boldsymbol{\Xi}}\,\text{d}\boldsymbol{z}\right)\text{d}\boldsymbol{\xi}\right]=0,\end{split}$
(17)

for the Langevin and Liouville approaches. Relation (17) can be used to
compare both Fokker-Planck and Liouville approaches. We note that the
Liouville approach has additional degrees as compared to the Fokker-Planck
approach. By comparing moments of both approaches, additional relations can be
found.

##  3 Analytical solutions to Liouville models

###  3.1 Particle trajectory

Lattanzi et al. [41] considered a one-dimensional version of the Fokker-Planck
equation (7) with constant velocity $u_{\text{p}}$ and diffusion tensor
${b_{x}}_{ij}=\sqrt{2D}\delta_{ij}$,

$\displaystyle\frac{\partial f_{X}}{\partial t}+u_{\text{p}}\frac{\partial
f_{X}}{\partial x_{\text{p}}}=D\frac{\partial^{2}f_{X}}{\partial
x_{\text{p}}^{2}},$ (18)

subject to the deterministic initial condition,
$f_{X}(x_{\text{p}};0)=\delta(x_{\text{p}})$. Its solution,

$\displaystyle f_{X}(x_{\text{p}};t)=\frac{1}{\sqrt{4\pi
Dt}}\exp\left(-\frac{\left(x_{\text{p}}-u_{\text{p}}t\right)^{2}}{4Dt}\right),$
(19)

provides a probabilistic prediction of the particle trajectory,
$X_{\text{p}}(t)$, within the Langevin framework. As expected, this solution
is Gaussian.

In the corresponding one-dimensional Liouville approach, we search for
relations to define the random forcings such that the solution of the Fokker-
Planck equation and the Liouville equation coincide. Starting from (18) and
(19), we take the first moment of the Fokker-Plank solution and compare that
with the first moment of the Liouville approach. In particular, if we employ
the law of the unconscious statistician (LOTUS) [65] to the characteristics of
the Liouville equation (9),

$\displaystyle\frac{\text{d}\hat{x}_{\text{p}}}{\text{d}t}=u_{\text{p}}+\sum_{i=1}^{N}z_{i}\varphi_{i}(t),$
(20)

to find the first moment of the particle position which according to (19) is
$\bar{X}_{\text{p}}=u_{\text{p}}t$, we find

$\displaystyle\bar{X}_{\text{p}}=u_{\text{p}}t+\int\left(\sum_{i=1}^{N}z_{i}\int_{0}^{t}\varphi_{i}(t^{\prime})\text{d}t^{\prime}\right)f_{\boldsymbol{Z}}(\boldsymbol{z})\text{d}\boldsymbol{z},$
(21)

where we have applied the deterministic initial condition $X_{\text{p}}(0)=0$,
according to the heat kernel solution (19). The second term in the right hand
side of (21) must be zero to match the solution of the Fokker-Planck equation.
We therefore find a condition to match the first moment

$\displaystyle\int\left(\sum_{i=1}^{N}z_{i}\int_{0}^{t}\varphi_{i}(t^{\prime})\text{d}t^{\prime}\right)f_{\boldsymbol{Z}}(\boldsymbol{z})\text{d}\boldsymbol{z}=0.$
(22)

Then, looking at the second moment of the solution (19),
$\overline{{X_{\text{p}}^{\prime}}^{2}}=2Dt$, using the same procedure we find
the condition

$\displaystyle\int[u_{\text{p}}t+\sum_{i=1}^{N}{z}_{i}\int_{0}^{t}\varphi_{i}(t^{\prime})\text{d}t^{\prime}-\bar{X}_{\text{p}}^{2}f_{\boldsymbol{Z}}(\boldsymbol{z})\text{d}\boldsymbol{z}=2Dt.$
(23)

Conditions (22) and (23) can be fulfilled by considering a single random
coefficient $\Xi\equiv Z_{i}$, for $i=1,\dots,N$, which when used in (22)
leads to

$\displaystyle\bar{\Xi}\sum_{i=1}^{N}\int_{0}^{t}\varphi_{i}(t^{\prime})\text{d}t^{\prime}=0,\qquad\bar{\Xi}=0.$
(24)

Then, applying this to also (23) gives

$\displaystyle\sigma_{\Xi}\sum_{i=1}^{N}\int_{0}^{t}\varphi_{i}(t^{\prime})\text{d}t^{\prime}=\sqrt{2Dt},$
(25)

that can be fulfilled by considering a single function
$\varphi\equiv\varphi_{i}$ such that

$\displaystyle\varphi=\sqrt{\frac{D}{2t}},\qquad\sigma_{\Xi}=1.$ (26)

Note that in general, basis function may be used to approximate the source by
a finite number of terms $\sum_{i=1}^{N}\varphi_{i}(t)\approx\sqrt{D/(2t)}$.
The random source term can be simply described by a single term
$\Xi\,\varphi(t)$ such that (9) reduces to

$\displaystyle\frac{\partial f_{X\Xi}}{\partial
t}+\left(u_{\text{p}}+\xi\sqrt{\frac{D}{2t}}\right)\frac{\partial
f_{X\Xi}}{\partial x_{\text{p}}}=0.$ (27)

Since the initial position of the particle, $X_{p}(0)=0$ is known with
certainty, i.e., deterministic, this Liouville equation is subject to the
initial condition
$f_{X\Xi}(x_{\text{p}},\xi;0)=f_{\Xi}(\xi)\delta(x_{\text{p}})$. Its solution,
obtained via the method of characteristics in Appendix C, is

$\displaystyle
f_{X\Xi}(x_{\text{p}},\xi;t)=f_{\Xi}(\xi)\delta(x_{\text{p}}-u_{\text{p}}t-\xi\sqrt{2Dt}).$
(28)

The PDF of the particle position, $f_{X}(x_{\text{p}};t)$, is the marginal of
$f_{X\Xi}$ computed by integrating (28) over $\xi$, while accounting for the
properties of the the Dirac delta function $\delta(-x)=\delta(x)$ and
$|\gamma|\delta(\gamma x)=\delta(x)$ for a given scalar $\gamma$,

$\displaystyle\begin{split}f_{X}(x_{\text{p}};t)=\frac{1}{\sqrt{2Dt}}f_{\Xi}\left(\frac{x_{\text{p}}-u_{\text{p}}t}{\sqrt{2Dt}}\right).\end{split}$
(29)

Figure 1: Particle position PDF $f_{X}$ at times $t=[0,\ 0.05,\ 0.8]$ given by
(29) with $D=2.5$ and $u_{\text{p}}=5$ for normal and triangular distributions
of the parameter $\Xi$. The uniform distribution has been omitted for clarity
(see Fig. 2b).

If $f_{\Xi}(\xi)$ is standard normal, $\Xi\sim\mathcal{N}(0,1)$, then the
Liouville solution (29) reduces to the Gaussian solution (19) predicted by the
Langevin approach. Otherwise, the PDF $f_{X}$ in (29) is non-Gaussian (Fig.
1). The temporal evolution of the particle-trajectory PDF
$f_{X}(x_{\text{p}};t)$ is depicted in Figure 2 for the normal,
$\mathcal{N}(0,1)$; uniform, $\mathcal{U}(-\sqrt{3},\sqrt{3})$; and
triangular, $\mathcal{T}(-2\sqrt{2},\sqrt{2})$; distributions of the parameter
$\Xi$. These three distributions have the same mean, 0, and variance, 1, but
different higher moments. The normal distribution gives rise to the Gaussian
$f_{X}(x_{\text{p}};t)$, associated with the Wiener process. On the other
hand, the triangular distribution translates into a skewed (non-Gaussian) PDF
of the particle trajectory. The uniform distribution provides the same
statistics than the normal up to the second moments, despite both PDFs differ
at any given time (Figs. 2a and 2b).

(a)

(b)

(c)

Figure 2: Temporal evolution of $f_{X}(x_{\text{p}};t)$, PDF of the particle
trajectory $X_{\text{p}}(t)$, predicted by the Liouville solution (29) for (a)
normal, (b) uniform and (c) triangular distributions of the parameter $\Xi$.
The deterministic parameters are set to $D=2.5$ and $u_{\text{p}}=5$.

The compatibility condition (11) may be employed to obtain the diffusion
tensor $\mathcal{D}$ for the Fokker-Planck equation from the results obtained
with the Liouville approach. In fact, the substitution of (28) and (29) into
(18) proves the identity, showing consistency between the Langevin and
Liouville approaches.

###  3.2 Particle velocity

The so-called velocity Langevin model [41] deals with a one-dimensional
version of (13),

$\displaystyle\text{d}{{X}_{\text{p}}}=\tau_{\text{p}}{{U}_{\text{p}}}\text{d}t,\qquad\text{d}{{U}_{\text{p}}}=-{{U}_{\text{p}}}\text{d}t+\sqrt{2D}\text{d}{W},$
(30)

in a reference frame moving at the average particle velocity [41]. The
corresponding Fokker-Planck equation for the joint PDF
$f_{XU}(x_{\text{p}},u_{\text{p}};t)$ of the particle trajectory
$X_{\text{p}}(t)$ and velocity $U_{\text{p}}(t)$ is

$\displaystyle\frac{\partial f_{XU}}{\partial
t}+\tau_{\text{p}}{u}_{\text{p}}\frac{\partial f_{XU}}{\partial
x_{\text{p}}}-\frac{\partial}{\partial
u_{\text{p}}}\left(u_{\text{p}}f_{XU}\right)=D\frac{\partial^{2}f_{XU}}{\partial
u_{\text{p}}^{2}}.$ (31)

If the initial particle position and velocity are known with certainty,
$X_{\text{p}}(0)=0$ and $U_{\text{p}}(0)=v_{0}$, then (31) is subject to the
initial condition
$f_{XU}(x_{\text{p}},u_{\text{p}};0)=\delta(x_{\text{p}})\delta(u_{\text{p}}-v_{0})$.
Integration of this equation over $x_{\text{p}}$ yields the Fokker-Planck
equation for the marginal $f_{U}(u_{\text{p}};t)$.

$\displaystyle\frac{\partial f_{U}}{\partial t}-\frac{\partial}{\partial
u_{\text{p}}}\left(u_{\text{p}}f_{U}\right)=D\frac{\partial^{2}f_{U}}{\partial
u_{\text{p}}^{2}},$ (32)

subject to the initial condition
$f_{U}(u_{\text{p}};0)=\delta(u_{\text{p}}-v_{0})$. This problem admits the
Gaussian solution [60, 66, 67, 41]

$\displaystyle f_{U}(u_{\text{p}};t)=\frac{1}{\sqrt{2\pi
D(1-\text{e}^{-2t})}}\exp\left[-\frac{(u_{\text{p}}-v_{0}\text{e}^{-t})^{2}}{2D(1-\text{e}^{-2t})}\right].$
(33)

Our Liouville approach replaces the Langevin model (30) with the one-
dimensional version of (15). To find the random forcings that make coincide
both approaches, we proceed similar as in the previous case (see Section 3.1).
Equivalently to the use of the LOTUS with the characteristics of the Liouville
equation to obtain moments, the model developed in Appendix A, can also be
used to find expressions for the moments using the Liouville approach. In
essence, we compare the moments of the solution of the Fokker-Planck equation
(31) to find the random forcings in the Liouville approach. As in the previous
example, we assume a single random coefficient can be used
$\Xi\equiv\Xi_{i}=Z_{i}$ and the basis functions reduce to a deterministic
function of time, $\varphi(t)$ and $\phi(t)$ for the position and velocity
equations respectively. Then, using the moment equations (63) derived in
Appendix B.2, one can find the closed-form functions

$\displaystyle\begin{split}\varphi=\frac{\text{d}\sigma_{X}}{\text{d}t}-\tau_{p}\sigma_{U}=\tau_{\text{p}}\sqrt{D}\left(\frac{1-2\text{e}^{-t}+\text{e}^{-2t}}{\sqrt{2t-3+4\text{e}^{-t}-\text{e}^{-2t}}}-\sqrt{1-\text{e}^{-2t}}\right),\end{split}$
(34a) $\displaystyle\phi$
$\displaystyle=\frac{\text{d}\sigma_{U}}{\text{d}t}+\sigma_{U}=\sqrt{\frac{D}{1-\text{e}^{-2t}}},$
(34b)

where we have taken as in the previous example, $\bar{\Xi}=0$, and
$\sigma_{\Xi}=1$. The functions $\sigma_{X}(t)$ and $\sigma_{U}(t)$ are the
standard deviations of $X_{\text{p}}(t)$ and $U_{\text{p}}(t)$ according to
the Fokker-Planck equation (31), given by[41]

$\displaystyle\sigma_{X}$
$\displaystyle=\sqrt{\tau_{\text{p}}^{2}D\left(2t-3+4\text{e}^{-t}-\text{e}^{-2t}\right)},$
(35a) $\displaystyle\sigma_{U}$
$\displaystyle=\sqrt{D\left(1-\text{e}^{-2t}\right)}.$ (35b)

In general, data-driven techniques can be used to adjust the random forcings
in (3) by defining an optimization problem. See for example Ref. [56] where
the PDF of the random coefficients are learned from DNS with an adjoint-based
data assimilation procedure using Fourier series as basis functions.

With the random forcings defined by the random coefficient $\Xi$ and the
functions in (34), our Liouville approach replaces the Langevin model (30)
with

$\displaystyle\frac{\text{d}{X}_{\text{p}}}{\text{d}t}=\tau_{\text{p}}U_{\text{p}}+\Xi\varphi,\qquad\frac{\text{d}{U}_{\text{p}}}{\text{d}t}=-U_{\text{p}}+\Xi\phi,$
(36)

The Liouville equation for the joint PDF
$f_{XU\Xi}(x_{\text{p}},u_{\text{p}},\xi;t)$ is

$\displaystyle\begin{split}\frac{\partial f_{XU\Xi}}{\partial
t}&+\frac{\partial}{\partial
x_{\text{p}}}\left[\left(\tau_{\text{p}}u_{\text{p}}+\xi\varphi\right)f_{XU\Xi}\right]+\frac{\partial}{\partial
u_{\text{p}}}\left[\left(-u_{\text{p}}+\xi\phi\right)f_{XU\Xi}\right]=0,\end{split}$
(37)

subject to the initial condition
$f_{XU\Xi}(x_{\text{p}},u_{\text{p}},\xi;0)=\delta(x_{\text{p}})\delta(u_{\text{p}}-v_{0})f_{\Xi}(\xi)$,
which gives rise to the solution (Appendix C)

$\displaystyle\begin{split}f_{XU\Xi}=\text{e}^{t}f_{\Xi}(\xi)\delta[x_{\text{p}}-\tau_{\text{p}}v_{0}(1-\text{e}^{-t})-\xi\sigma_{X}]\,\delta[\text{e}^{t}(u_{\text{p}}-\xi\sigma_{U})-v_{0}].\end{split}$
(38)

The marginals of (38) are the PDFs of the particle position and velocity,

$\displaystyle
f_{X}(x_{\text{p}};t)=\frac{1}{\sigma_{X}}f_{\Xi}\left(\frac{x_{\text{p}}-\bar{X}_{\text{p}}}{\sigma_{X}}\right),$
(39)

and

$\displaystyle f_{U}(u_{\text{p}};t)$
$\displaystyle=\frac{1}{\sigma_{U}}f_{\Xi}\left(\frac{u_{\text{p}}-\bar{U}_{\text{p}}}{\sigma_{U}}\right),$
(40)

where $\bar{X}_{\text{p}}(t)$ and $\bar{U}_{\text{p}}(t)$ are the means of
$X_{\text{p}}(t)$ and $U_{\text{p}}(t)$, respectively (given in (64)).
Solutions (39) and (40) are generalizations of the velocity Langevin model
[41], with (40) reducing to (33) when the random variable $\Xi$ is standard
normal, $\Xi\sim\mathcal{N}(0,1)$. The temporal evolution of the PDFs
$f_{X}(x_{\text{p}};t)$ and $f_{U}(u_{\text{p}};t)$ is depicted in Figure 3,
for the normal, uniform and triangular distributions of $\Xi$. The spatial
variations of the PDFs of particle positions and velocity related to non-
Gaussian effects (Figs. 3b– 3c and 3e– 3f) can not be represented by the
Fokker-Plank model (31) (Figs. 3a and 3d). The Liouville approach may be used
then to formulate non-Gaussian stochastic models.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: Temporal evolution of $f_{X}(x_{\text{p}};t)$ and
$f_{U}(u_{\text{p}};t)$, PDFs of the particle trajectory $X_{\text{p}}(t)$ and
velocity $U_{\text{p}}(t)$, predicted by the Liouville solutions (39) and (40)
for normal (left column), uniform (middle column) and triangular (right
column) distributions of the parameter $\Xi$. The deterministic parameters are
set to $\tau_{\text{p}}=10$, $v_{0}=1$ and $D=1/50$.

When applied to (37), the method of characteristics provides a geometric
interpretation of the dynamics of the joint PDF
$f_{XU\Xi}(x_{\text{p}},u_{\text{p}},\xi;t)$ and its marginals. The joint PDF
$f_{XU\Xi}$ is transported along the characteristics
$\boldsymbol{y}(t)=[\hat{x}_{\text{p}}(t),\hat{u}_{\text{p}}(t),\xi]^{\top}$,
defined by (69). This process defines a flow map [55] and is akin to mapping
the initial state $f_{XU\Xi}(x_{\text{p}},u_{\text{p}},\xi;0)$ onto its
counterpart at time $t$ via the method of transformations [63]:

$\displaystyle
f_{XU\Xi}(x_{\text{p}},u_{\text{p}},\xi;0)=J(t)f_{XU\Xi}(x_{\text{p}},u_{\text{p}},\xi;t),\quad
J=\left|\frac{\partial\,{\boldsymbol{y}}(t)}{\partial\,{\boldsymbol{y}}(0)}\right|,$
(41)

where $J(t)$ the determinant of the transformation Jacobian. Its direct
evaluation from the characteristics equation (69a) yields $J(t)=\exp(-t)$,
which is consistent with (38) and (69b). Numerical evaluation of the Jacobian
and the use of the method of transformations, i.e., using (41), provides an
alternative to computing the solution of (69b).

##  4 A Liouville model for fluidized homogeneous heating systems

In a fluidized homogeneous heating system (FHHS), particles with diameter
$d_{\text{p}}$ are released into a current of carrier-fluid (with density
$\rho_{\text{f}}$ and viscosity $\mu_{\text{f}}$). As time increases, the
particle velocity $\boldsymbol{U}_{\text{p}}(t)$ deviates from the fluid
velocity in such a way that the variances of the particle velocity along each
component reach a steady value. The quantities

$\displaystyle
T=\frac{1}{3}\overline{\boldsymbol{U}_{\text{p}}^{\prime}\cdot\boldsymbol{U}_{\text{p}}^{\prime}}\quad\text{and}\quad\mathit{Re}_{T}=\frac{\rho_{\text{f}}\,d_{\text{p}}}{\mu_{\text{f}}}\sqrt{T}$
(42)

are referred to as the granular temperature and the thermal Reynolds number,
respectively [43]. Notice that a Reynolds decomposition is used for
$\boldsymbol{U}_{\text{p}}=\bar{\boldsymbol{U}}_{\text{p}}+\boldsymbol{U}_{\text{p}}^{\prime}$.
They have been studied analytically via kinetic theory [57] or empirical
treatment of Langevin models [58] and numerically via PR-DNS [39] or Langevin
equations [43]. We use these results to ascertain the veracity of a new model
based on our Liouville approach.

In a reference frame moving with the average particle velocity [41, 43], a
one-dimensional Liouville model of the random particle velocity
$U_{\text{p}}(t)$ is

$\displaystyle\frac{\text{d}U_{\text{p}}}{\text{d}t}$
$\displaystyle=-\frac{1}{\tau_{\text{p}}}U_{\text{p}}+\Xi\phi,$ (43)

such that $T=\sigma_{U}^{2}/3$. The random variable $\Xi$ has zero mean,
$\bar{\Xi}=0$, and standard deviation
$\sigma_{\Xi}=\sigma_{\Xi}(\mathit{Re}_{\text{m}})$, which depends on the mean
Reynolds number in the PR-DNS [43],

$\displaystyle\mathit{Re}_{\text{m}}=(1-\omega)\frac{\rho_{\text{f}}d_{\text{p}}}{\mu_{\text{f}}}|\boldsymbol{V}|,$

where $\boldsymbol{V}$ is the average relative velocity, and $\omega$ is the
average particle solid volume fraction. The deterministic function,

$\displaystyle\phi=\frac{1}{\tau_{\text{p}}}\left(1-\text{e}^{-\mathcal{C}_{1}t/\tau_{\text{p}}}\right)^{\mathcal{C}_{2}-1}\left[1+\left(\mathcal{C}_{1}\mathcal{C}_{2}-1\right)\text{e}^{-\mathcal{C}_{1}t/\tau_{\text{p}}}\right]$
(44)

is chosen phenomenologically, with the constants $\mathcal{C}_{1}$ and
$\mathcal{C}_{2}$ characterizing the behavior of the granular temperature $T$
at early and late times, respectively. Both constants depend on
$\mathit{Re}_{\text{m}}$, $\rho_{\text{p}}/\rho_{\text{f}}$ and $\omega$, and
are adjusted to match the PR-DNS results [39, 58, 43].

The Liouville equation for the joint PDF $f_{U\Xi}(u_{\text{p}},\xi;t)$ is

$\displaystyle\frac{\partial f_{U\Xi}}{\partial t}+\frac{\partial}{\partial
u_{\text{p}}}\left[\left(-\frac{1}{\tau_{\text{p}}}u_{\text{p}}+\xi\phi\right)f_{U\Xi}\right]=0.$
(45)

Its solution is (Appendix C)

$\displaystyle\begin{split}f_{U\Xi}=f_{U\Xi}^{0}\text{e}^{t/\tau_{\text{p}}}=\text{e}^{t/\tau_{\text{p}}}f_{\Xi}(\xi)\delta[(u_{\text{p}}-\xi\eta)\text{e}^{t/\tau_{\text{p}}}]=\frac{1}{\eta}f_{\Xi}(\xi)\delta\left(\frac{u_{\text{p}}}{\eta}-\xi\right),\end{split}$
(46a) where
$\displaystyle\eta(t)=(1-\text{e}^{-\mathcal{C}_{1}t/\tau_{\text{p}}})^{\mathcal{C}_{2}}.$
(46b)

The marginalization along $\xi$ gives the PDF of the particle velocity
$U_{\text{p}}(t)$ in FHHS,

$\displaystyle
f_{U}(u_{\text{p}};t)=\int\frac{1}{\eta}f_{\Xi}(\xi)\delta\left(\frac{u_{\text{p}}}{\eta}-\xi\right)\text{d}\xi=\frac{1}{\eta}f_{\Xi}\left(\frac{u_{\text{p}}}{\eta}\right).$
(47)

In accordance with (42), the variance of this PDF, $\sigma_{U}^{2}(t)$, is
proportional to granular temperature $T(t)$. If $\Xi$ were Gaussian,
$f_{\Xi}\sim\mathcal{N}(0,\sigma_{\Xi}^{2})$, then

$\displaystyle\bar{U}_{\text{p}}(t)={\bar{U}_{\text{p}}}(0)\text{e}^{-t/\tau_{\text{p}}}\quad\text{and}\quad\sigma_{U}^{2}(t)=\sigma_{\Xi}^{2}(1-\text{e}^{-\mathcal{C}_{1}t/\tau_{\text{p}}})^{2\mathcal{C}_{2}}.$
(48)

Figure 4: Temporal evolution of the thermal Reynolds number $Re_{T}(t)$
alternatively predicted with our model (48) (solid lines with square symbol),
the model of Ref.[58] (dotted lines with triangle symbol) and the Langevin
approach of Ref.[43] (dashed lines with diamond symbol), for several values of
the mean Reynolds number $\mathit{Re}_{\text{m}}$. The open circles indicate
the PR-DNS data [39, 43], to which
$\sigma_{\Xi}=\sigma_{\Xi}(\mathit{Re}_{\text{m}})$ and
$\mathcal{C}_{1}=\mathcal{C}_{1}(\mathit{Re}_{\text{m}})$ are fitted (Fig. 5).
The average volume fraction and density ratio are $\omega=0.1$ and
$\rho_{\text{p}}/\rho_{\text{f}}=100$ respectively. For all cases
$\mathcal{C}_{2}=1.2$ and $\tau_{\text{p}}=0.14$.

Figure 4 exhibits the granular temperature predicted with this solution,
$T(t)=\sigma_{U}^{2}/3$, for $\tau_{\text{p}}=0.14$,
$\rho_{\text{p}}/\rho_{\text{f}}=100$ and $\omega=0.1$ [43]. The functional
dependencies $\sigma_{\Xi}=\sigma_{\Xi}(\mathit{Re}_{\text{m}})$ and
$\mathcal{C}_{1}=\mathcal{C}_{1}(\mathit{Re}_{\text{m}})$ are obtained by
fitting the predictions of $T(t)$ from (48) to the PR-DNS data from [39]. This
fitting reveals that both $\sigma_{\Xi}$ and $\mathcal{C}_{1}$ depend linearly
in a log-log scale on $\mathit{Re}_{\text{m}}$ (Fig. 5), whereas
$\mathcal{C}_{2}$ is constant ($\mathcal{C}_{2}=1.2$). Our model matches the
PR-DNS data better than its Langevin alternatives [43, 58]. It captures both
the early growth and asymptotic values of $T(t)$, while the models of
Refs.[58] and [43] underestimate the steady values of $T(t)$ and predict the
unphysical/unobserved long-term fluctuations in $T(t)$, respectively.

Figure 5: Functional dependencies of the model parameters $\sigma_{\Xi}$ and
$\mathcal{C}_{1}$ on the mean Reynolds number $\mathit{Re}_{\text{m}}$. The
squares and circles indicate the values of $\sigma_{\Xi}$ and
$\mathcal{C}_{1}$ obtained via fitting (48) to the PR-DNS data [39] for
$\rho_{\text{p}}/\rho_{\text{f}}=100$ and $\omega=0.1$. These data are
represented by
$\log(\sigma_{\Xi})=0.06258\log(\mathit{Re}_{\text{m}})+\log(0.7866)$ and
$\log(\mathcal{C}_{1})=2.446\log(\mathit{Re}_{\text{m}})+\log(0.5411)$, with
the coefficient of determination $R^{2}=0.998$.

###  4.1 Phenomenology of granular-temperature dynamics

It has been argued that granular temperature in a FHHS, $T(t)$, exhibits the
source-and-sink dynamics [57, 39, 58],

$\displaystyle\frac{\text{d}T}{\text{d}t}=S(T)-\Gamma(T).$ (49)

In Ref.[57], the source is taken to be $S\sim T^{-1/2}$, which is singular at
$T=0$, i.e., at initial time $t=0$. In contrast, the ensemble averaging of our
Liouville model (43) gives rise to the source $S(T)$ and sink $\Gamma(T)$ (see
Appendix A for detail),

$\displaystyle S=\frac{2}{\sqrt{3}}\sigma_{\Xi}\phi
T^{1/2}\quad\text{and}\quad\Gamma=\frac{2}{\tau_{\text{p}}}T,$ (50)

that behave well at $T=0$. Figure 6 exhibits the curves $S(T)$ and $\Gamma(T)$
given by (50) and by the model in Ref.[57], in comparison with PR-DNS data
[39]. While the mismatch between both models and the data is noticeable, our
model exhibit the correct qualitative behavior, while the model in Ref.[57]
does not.

Figure 6: The source, $S(T)$, and sink, $\Gamma(T)$, terms in the granular-
temperature model (49). Our model (50) yields the solid and dashed green lines
for $S(T)$ and $\Gamma(T)$, respectively. The model of Ref.[57] has $S(T)$ and
$\Gamma(T)$ depicted by the dash-dotted and dotted black lines, respectively.
The PR-DNS estimates of $S(T)$ and $\Gamma(T)$[39] are represented by the non-
filled and filled red triangles, respectively. The parameter values are
$\mathit{Re}_{\text{m}}=20$, $\rho_{\text{p}}/\rho_{\text{f}}=100$ and
$\omega=0.1$.

###  4.2 Particle-velocity distributions

A Langevin counterpart to our Liouville model (43) is a stochastic ODE
$\text{d}{U_{\text{p}}}=-(1/\tau_{\text{p}}){U_{\text{p}}}\text{d}t+\sqrt{2D}\text{d}W$,
which gives rise to the Fokker-Planck equation for the particle velocity PDF,
$f_{U}(u_{\text{p}};t)$,

$\displaystyle\frac{\partial f_{U}}{\partial t}-\frac{\partial}{\partial
u_{\text{p}}}\left(\frac{1}{\tau_{\text{p}}}u_{\text{p}}f_{U}\right)=\frac{\partial}{\partial
u_{\text{p}}}\left(D\frac{\partial f_{U}}{\partial u_{\text{p}}}\right).$ (51)

Many studies focus on the numerical determination of the diffusion coefficient
$D(u_{\text{p}},t)$ using the standard deviation of particle acceleration and
its integral time scale [43] or directly fitting a correlation from PR-DNS
data [39]. Our Liouville approach provides an alternative means to accomplish
this task analytically. Substituting $f_{U}(u_{\text{p}};t)$ from (47) into
(51) yields a second-order ODE for $D(u_{\text{p}},\cdot)$, whose solution is

$\displaystyle
D=\eta\phi\sigma_{\Xi}^{2}+\frac{K}{u_{\text{p}}}\exp\left(\frac{u_{\text{p}}^{2}}{2\eta^{2}\sigma_{\Xi}^{2}}\right),$
(52)

with $K$ an integration constant whose value is not constrain by the
compatibility condition (see (11) and (17)), and it is an extra degree of
freedom of the model. Here we choose $K=0$, such that the diffusion
coefficient is simply $D=\eta\phi\sigma_{\Xi}^{2}$. Using the analytical
relations derived previously, the diffusion tensor can be expressed in terms
of the granular temperature as

$\displaystyle
D=\frac{3}{\tau_{\text{p}}}\left[\left(1-\mathcal{C}_{1}\mathcal{C}_{2}\right)T+\frac{\mathcal{C}_{1}\mathcal{C}_{2}\sigma_{\Xi}^{1/\mathcal{C}_{2}}}{3^{1/(2\mathcal{C}_{2})}}T^{1-1/(2\mathcal{C}_{2})}\right].$
(53)

Models in the literature [39, 43] and, in general, numerical approximations of
the diffusion coefficient in the Langevin approach are not consistently linked
to Fokker-Planck equations and an equivalent PDF formulation is seldom
investigated. Moreover, the full statistical description is provided by a PDF
description (46a). The model presented here with the Liouville approach leads
to equivalent Langevin and Fokker-Planck formulations. In Fig 7 we compare the
derived diffusion coefficient (53) with the model proposed in Ref.[39],
showing agreement at lower values of the granular temperature, linked to the
evolution at the early times of the FHHS. The model in Ref.[39] is based on a
correlation that contains $8$ constants fitted from PR-DNS data.

Figure 7: Diffusion coefficient in terms of the granular temperature in the
model presented in Ref.[39] and the one developed here, equation (53), for
$\mathit{Re}_{m}=20$, $\omega=0.1$, $\rho_{\text{p}}/\rho_{\text{f}}=100$ and
$\tau_{\text{p}}=0.14$.

###  4.3 Miscellaneous musings

The advantage of the Liouville approach is obvious, the source and sink terms
for the evolution of the granular temperature can be proposed starting from
defining the standard deviations of the particle velocity (48) based on
phenomenology. Then, we proceed to find the source term $\phi$ going back to
the characteristics that can be analytically described. The model also leads a
diffusion coefficient that can be computed accordingly by comparing the
Fokker-Planck equation with the marginalized Liouville equation. Trying to
formulate the problem from the Langevin perspective without previous knowledge
of the diffusion coefficient and arrive at the expression in (52) is a hard
task. The procedure can be used to develop a model based on the time evolution
of $T$ reported in PR-DNS results by adjusting the parameters $\sigma_{\Xi}$,
$\mathcal{C}_{1}$ and $\mathcal{C}_{2}$.

The temporal evolution of the granular temperature in FHHS is characterized by
a source and sink behavior [39] stabilizing at a constant value after an
initial growth. To this purpose, the parameter $\sigma_{\Xi}$ is responsible
for the steady value, and the coefficients $\mathcal{C}_{1}$ and
$\mathcal{C}_{2}$ characterize the growth of the granular temperature at two
distinguished time scales, at early and later times respectively. The Langevin
models in Ref.[41] correspond to ODE systems with a random coefficient in the
Liouville approach with source terms not defined at the initial time, making
the solution singular at $t=0$, corresponding to an initial infinite impulse
of the particles (events moving in phase space). See the previous Sections 3.1
and 3.2 that reproduce results reported in Ref.[41]. Another model where the
solution is not defined at the initial time is the one presented in Ref.[57],
where the granular temperature grows infinitely at the initial time. In our
model, the initial impulse in the granular temperature is bounded and governed
by the parameter $\mathcal{C}_{2}$, leading to an ODE system whose solution is
bounded at $t=0$.

##  5 Conclusions

A framework to formulate stochastic descriptions of particle-laden flows based
on a Liouville PDF governing equation introduced. A connection between
Langevin dynamical systems and the new Liouville approach allows to
analytically treat the solution with the method of characteristics and
circumvent numerical difficulties related to Itô calculus and the numerical
resolution of the Fokker-Planck equation. The Liouville approach is based on
expressing stochastic processes with random coefficients constants in time and
additional source terms. The problem is then governed by the Liouville
equation that can be derived with the method of distributions in an augmented
dimensionality that includes the added random coefficients, the particle
position and particle velocity. With proper selection of the random
coefficients and source terms, both Langevin and Liouville approaches
coincide. This implies that the Liouville equation marginalized along the
extra dimension of the random coefficients leads to the Fokker-Planck equation
of an equivalent Langevin system.

The Liouville approach is compared with canonical Eulerian-Lagrangian models
of particle-laden flows formulated with Langevin equations used to obtain the
time evolution of particle dispersion and granular temperature in the
literature. We develop a full description of the statistics with the Liouville
approach, that analytically describes exact results of the PDFs of the
particle position and velocity as well as their moments.

The potential of the Liouville approach to develop modeling of particle
statistics is assessed by proposing a phenomenological model to FHHS. Based on
the analytical character of the Liouville approach, the sink and source terms
that govern the time evolution of the granular temperature can be represented
in the Liouville approach by selecting a source term that depends on the non-
dimensional parameters of the problem. This leads to matching of PR-DNS
results more accurately than previous proposed numerical and analytical models
based on the Langevin approach. The FHHS model presented here also eliminates
the singularity at the initial time in the granular temperature equation, that
leads to an ill-posed problem that contains an initial infinite impulse.
Previous models in the literature did not offer descriptions of the granular
temperature without such singularity. This singularity at the initial time can
be interpreted as an initial infinite impulse by virtue of the Liouville
approach that links events with possible particles. Such events are associated
with deterministic characteristic equations that are derived from the
Liouville equation with the method of characteristics. The hyperbolicity of
the Liouville equation is therefore advantageous as opposed to the advection-
diffusion character of the Fokker-Planck equation. The model presented here is
well posed at the initial time and describes the early evolution of the
granular temperature accurately as compared with PR-DNS. Additionally, the
Liouville approach may be used to find the diffusion coefficient and with it
other descriptions of the same problem based on the Fokker-Planck or Langevin
equations.

The Liouville approach may be used to describe other Langevin dynamical
systems in general, making its numerical treatment simpler. It is not bounded
by Gaussianity and does not rely on sampling nor requires closure. It opens
new possibilities in proposing reduced models that capture statistics obtained
with PR-DNS. Moreover, because events and their probability to occur are
traced deterministically, each event has a deterministic equation given by the
characteristics that can be physically interpreted as a possible particle in
the system. Then each event is associated with a physical particle and how
likely to occur is that particle in the system. Therefore it provides a
connection between the statistics and the particle physics.

Additionally, moment equations provide a reduced description, which may lead
to closure challenges depending on the non-linearity of the particle
equations. In the results presented here, the moment equations are found in
closed form.

Extensions of this work include accounting for mass and energy exchange
between the two phases, as well as the development of stochastic models for
particle statistics in flows with particle-shock interactions [68, 69, 70].

## Funding

This work was supported by the Air Force Office of Scientific Research under
award numbers FA9550-19-1-0387 and FA9550-21-1-0381, and by San Diego State
University Graduate Fellowship.

## Declaration of interests

The authors report no conflict of interest.

## Data Availability Statement

The data that support the findings of this study are available upon reasonable
request.

##  Appendix A Moment Models

With the Liouville approach, descriptions of moments of the particle phase
variables may also be derived using the systems (8) or (15). These include the
variance of particle positions (particle dispersion), the variance of particle
velocity (granular temperature) and higher moments as the third and fourth
moment of the particle position less reported in the literature. To write a
general moment model based on the systems of ODEs with a random coefficient
described previously, we first rewrite the ODE systems in terms of the vector
of state variables $\boldsymbol{q}$. Then, both systems (8) and (15) can be
written as

$\displaystyle\frac{\text{d}\boldsymbol{q}}{\text{d}t}=\boldsymbol{g}(\boldsymbol{q}),$
(54)

where $\boldsymbol{q}=(\boldsymbol{x}_{\text{p}},\ \Xi)\in\mathbb{R}^{4}$ for
the model in (8) and $\boldsymbol{q}=(\boldsymbol{x}_{\text{p}},\
\boldsymbol{u}_{\text{p}},\ \Xi)\in\mathbb{R}^{7}$ for (15) with the
corresponding function with right hand side terms here taken as
$\boldsymbol{g}$, and considering a single random coefficient $\Xi$. Then, we
employ a Reynolds decomposition of the random variables
$(\cdot)=\overline{(\cdot)}+(\cdot)^{\prime}$. Substituting in (54) and
averaging, we arrive to the first two moment equations

$\displaystyle\frac{\text{d}\overline{q}_{i}}{\text{d}t}$
$\displaystyle=\overline{g}_{i},$ (55a)
$\displaystyle\frac{\text{d}{\overline{{q_{i}^{\prime}}{q_{j}^{\prime}}}}}{\text{d}t}$
$\displaystyle=\overline{q_{j}^{\prime}g_{i}^{\prime}}+\overline{q_{i}^{\prime}g_{j}^{\prime}},$
(55b)

where unclosed terms appear as the average of the non-linear function
$\boldsymbol{g}$ is in principle unknown, as well as second order moments
involving it. Closures for such terms may be based on a priori computations
[71] or expansions [72, 54]. If the right hand side of the ODE system is
composed by a product of non-linear functions of the state variables,
$\boldsymbol{q}$, statistical moments of order higher appear in each moment
equation, defining a classic closure problem. Take for exmaple
$\boldsymbol{g}(\boldsymbol{q})=\boldsymbol{r}(\boldsymbol{q})\circ\boldsymbol{h}(\boldsymbol{q})$,
then the system (55) is rewritten as

$\displaystyle\frac{\text{d}\overline{q}_{i}}{\text{d}t}$
$\displaystyle=\overline{r}_{i}\overline{h}_{i}+\overline{r_{i}^{\prime}h_{i}^{\prime}},$
(56a)
$\displaystyle\frac{\text{d}{\overline{{q_{i}^{\prime}}{q_{j}^{\prime}}}}}{\text{d}t}$
$\displaystyle=\overline{q_{j}^{\prime}h_{i}^{\prime}}\overline{r}_{i}+\overline{q_{j}^{\prime}r_{i}^{\prime}}\
\overline{h}_{i}+\overline{q_{i}^{\prime}h_{j}^{\prime}}\overline{r}_{j}+\overline{q_{i}^{\prime}r_{j}^{\prime}}\
\overline{h}_{j}+\overline{q_{j}^{\prime}h_{i}^{\prime}r_{i}^{\prime}}+\overline{q_{i}^{\prime}h_{j}^{\prime}r_{j}^{\prime}},$
(56b)

with $\boldsymbol{r}$ and $\boldsymbol{h}$ two non-linear functions. This is
for example is the case if a correction of the Stokes drag $\mathcal{F}$ based
on the relative velocity is used in (15). By combining closures based on
expansions with truncation of higher order terms and a splitting algorithm, in
Ref.[54], the first two moments were computed in closed form.

When applied to (43), the method of moments yields

$\displaystyle\frac{\text{d}T}{\text{d}t}$
$\displaystyle=-\frac{2}{\tau_{\text{p}}}T+\frac{2}{3}\overline{\Xi^{\prime}U_{\text{p}}^{\prime}}\phi,$
(57a)
$\displaystyle\frac{\text{d}\overline{\Xi^{\prime}U_{\text{p}}^{\prime}}}{\text{d}t}$
$\displaystyle=-\frac{1}{\tau_{\text{p}}}\overline{\Xi^{\prime}U_{\text{p}}^{\prime}}+\sigma_{\Xi}^{2}\phi.$
(57b)

The solution of the second ODE is

$\displaystyle\overline{\Xi^{\prime}U_{\text{p}}^{\prime}}=\sigma_{\Xi}^{2}(1-\text{e}^{-\mathcal{C}_{1}t/\tau_{\text{p}}})^{\mathcal{C}_{2}}.$
(58a)

Since $T=\sigma_{U}^{2}/3$, with $\sigma_{U}^{2}$ given by (48), substituting
this result into (57a), yields the source $S$ and sink $\Gamma$ terms in (50).

##  Appendix B Analytical solutions based on the moment model

###  B.1 Particle trajectory

The moment equations of the one-dimensional version of (8) with a single
random variable $Z=\Xi$ read as

$\displaystyle\frac{\text{d}\overline{X}_{\text{p}}}{\text{d}t}$
$\displaystyle=\ u_{\text{p}}+\overline{\Xi}\varphi,$ (59a)
$\displaystyle\frac{\text{d}\overline{{X_{\text{p}}^{\prime}}^{2}}}{\text{d}t}$
$\displaystyle=2\overline{X_{\text{p}}^{\prime}U_{\text{p}}^{\prime}}+2\overline{\Xi^{\prime}X_{\text{p}}^{\prime}}\varphi,$
(59b)
$\displaystyle\frac{\text{d}\overline{\Xi^{\prime}X_{\text{p}}^{\prime}}}{\text{d}t}$
$\displaystyle=\overline{{\Xi^{\prime}}^{2}}\varphi.$ (59c)

Considering constant and deterministic velocity $u_{\text{p}}$, one has
$u_{\text{p}}=\overline{u}_{\text{p}}$ and $u_{\text{p}}^{\prime}=0$. The
impulse initial condition, given by a Dirac delta function located at zero,
gives
$\overline{X}_{\text{p}_{0}}=\overline{{X_{\text{p}}^{\prime}}_{0}^{2}}=0$.
Also, initially the random coefficient and the particle locations are
statistically independent such that
$\overline{\Xi^{\prime}X_{\text{p}}^{\prime}}_{0}=\sigma_{\Xi}\sigma_{{X_{\text{p}}}_{0}}=0$.
With this initial condition, the system (59) leads to the following solution

$\displaystyle\overline{X}_{\text{p}}$
$\displaystyle=u_{\text{p}}t+\overline{\Xi}\sqrt{2Dt},$ (60a)
$\displaystyle\overline{{{X_{\text{p}}^{\prime}}^{2}}}$
$\displaystyle=2\overline{{\Xi^{\prime}}^{2}}Dt,$ (60b)
$\displaystyle\overline{\Xi^{\prime}X_{\text{p}}^{\prime}}$
$\displaystyle=\overline{{\Xi^{\prime}}^{2}}\sqrt{2Dt},$ (60c)

where if $\overline{\Xi}=0$ and $\overline{{\Xi^{\prime}}^{2}}=1$, in
concordance with the classic heat kernel solution where
$\Xi\sim\mathcal{N}(0,1)$, the solution (60) leads to the classic result for
the particle dispersion $\overline{{X_{\text{p}}^{\prime}}^{2}}=2Dt$. The
average position grows linear in time $\overline{X}_{\text{p}}=u_{\text{p}}t$,
and the correlation of position with the random coefficient is
$\overline{\Xi^{\prime}X_{\text{p}}^{\prime}}=\sigma_{X_{\text{p}}}$. The
characteristic in equation (67a) can be recasted in the form
$\hat{x}_{\text{p}}=\overline{X}_{\text{p}}+\xi\sigma_{X_{\text{p}}}$
according to (60). This linear relation indicates that for each value of
$\xi$, a quantity $\xi\sigma_{X_{\text{p}}}$ is added to the average path and
causes that corresponding trajectory to deviate from the average.

The general expression for the average particle position can be also obtained
by averaging of the characteristic equation $\hat{x}_{\text{p}}(\xi,t)$ with
respect to $\xi$, which corresponds to the application of the law of the
unconscious statistician (LOTUS) [65]. Making use of (67a), with the impulse
located at zero $\hat{x}_{\text{p}}^{0}=0$, one has

$\displaystyle\begin{split}\overline{X}_{\text{p}}&=\int\hat{x}_{\text{p}}(\xi,t)f_{\Xi}(\xi)d\xi=\int\left(\hat{x}_{\text{p}}^{0}+u_{\text{p}}t+\xi\sqrt{2Dt}\right)f_{\Xi}(\xi)d\xi=u_{\text{p}}t+\overline{\Xi}\sqrt{2Dt}\end{split}$
(61)

and for the $n$-th central moment

$\displaystyle\overline{{X_{\text{p}}^{\prime}}^{n}}$
$\displaystyle=\int(\hat{x}_{\text{p}}(\xi,t)-\overline{X}_{\text{p}})^{n}f_{\Xi}(\xi)d\xi=\overline{{\Xi^{\prime}}^{n}}\left(2Dt\right)^{n/2},$
(62)

for $n\geq 2$. Notice that the even moments are non-zero for any of the chosen
distributions for $\Xi$ ($\mathcal{N}$, $\mathcal{U}$ or $\mathcal{T}$)
whereas the odds moments are zero for symmetric distributions. These result
generalize those reported in Ref.[41] for non-Gaussian white noise.

###  B.2 Particle velocity

The moment equations of the system (36) read as

$\displaystyle\frac{\text{d}\overline{X}_{\text{p}}}{\text{d}t}$
$\displaystyle=\tau_{\text{p}}\overline{U}_{\text{p}}+\overline{\Xi}\varphi,$
(63a) $\displaystyle\frac{\text{d}\overline{U}_{\text{p}}}{\text{d}t}$
$\displaystyle=-\overline{U}_{\text{p}}+\overline{\Xi}\phi,$ (63b)
$\displaystyle\frac{\text{d}\overline{{X^{\prime}_{\text{p}}}^{2}}}{\text{d}t}$
$\displaystyle=2\tau_{\text{p}}\overline{X_{\text{p}}^{\prime}U_{\text{p}}^{\prime}}+2\varphi\overline{\Xi^{\prime}X_{\text{p}}^{\prime}},$
(63c)
$\displaystyle\frac{\text{d}\overline{X^{\prime}_{\text{p}}U_{\text{p}}^{\prime}}}{\text{d}t}$
$\displaystyle=\tau_{\text{p}}\overline{{U_{\text{p}}^{\prime}}^{2}}-\overline{X_{\text{p}}^{\prime}U_{\text{p}}^{\prime}}+\varphi\overline{\Xi^{\prime}U_{\text{p}}^{\prime}}+\phi\overline{\Xi^{\prime}X_{\text{p}}^{\prime}},$
(63d)
$\displaystyle\frac{\text{d}\overline{{U_{\text{p}}^{\prime}}^{2}}}{\text{d}t}$
$\displaystyle=-2\overline{{U_{\text{p}}^{\prime}}^{2}}+2\phi\overline{\Xi^{\prime}U_{\text{p}}^{\prime}},$
(63e)
$\displaystyle\frac{\text{d}\overline{\Xi^{\prime}X_{\text{p}}^{\prime}}}{\text{d}t}$
$\displaystyle=\tau_{\text{p}}\overline{\Xi^{\prime}U_{\text{p}}^{\prime}}+\overline{{\Xi^{\prime}}^{2}}\varphi,$
(63f)
$\displaystyle\frac{\text{d}\overline{\Xi^{\prime}U_{\text{p}}^{\prime}}}{\text{d}t}$
$\displaystyle=-\overline{\Xi^{\prime}U_{\text{p}}^{\prime}}+\overline{{\Xi^{\prime}}^{2}}\phi,$
(63g)

where the equations are closed because in this reference frame the flow
velocity does not appear in the average velocity equation and the Stokes drag
has been assumed, $\mathcal{F}=1$. The solution to the system (63) considering
$\overline{\Xi}=0$ and $\sigma_{\Xi}=1$, according to
$\Xi\sim\mathcal{N}(0,1)$, and all other first and second moments zero at the
initial time except ${\overline{U}_{\text{p}}}_{0}=v_{0}$, as well as
statistical independence of the variables initially, is given by

$\displaystyle\overline{X}_{\text{p}}$
$\displaystyle=\tau_{\text{p}}v_{0}\left(1-e^{-t}\right),$ (64a)
$\displaystyle\overline{U}_{\text{p}}$ $\displaystyle=v_{0}e^{-t}$ (64b)
$\displaystyle\overline{{X^{\prime}_{\text{p}}}^{2}}$
$\displaystyle=\tau_{\text{p}}^{2}D\left(2t-3+4e^{-t}-e^{-2t}\right),$ (64c)
$\displaystyle\overline{X^{\prime}_{\text{p}}U_{\text{p}}^{\prime}}$
$\displaystyle=\tau_{\text{p}}D\sqrt{\left(1-e^{-2t}\right)\left(2t-3+4e^{-t}-e^{-2t}\right)},$
(64d) $\displaystyle\overline{{U_{\text{p}}^{\prime}}^{2}}$
$\displaystyle=D\left(1-e^{-2t}\right),$ (64e)
$\displaystyle\overline{\Xi^{\prime}X_{\text{p}}^{\prime}}$
$\displaystyle=\sqrt{\tau_{\text{p}}^{2}D\left(2t-3+4e^{-t}-e^{-2t}\right)},$
(64f) $\displaystyle\overline{\Xi^{\prime}U_{\text{p}}^{\prime}}$
$\displaystyle=\sqrt{D\left(1-e^{-2t}\right)}.$ (64g)

The averages particle position and velocity as well as the variances of
particle positions $\overline{{X_{\text{p}}^{\prime}}^{2}}$ and velocities
$\overline{{U_{\text{p}}^{\prime}}^{2}}$ is exact to the results reported in
Ref.[41] for the velocity Langevin model (see Figure 8a). In addition, the
solution (64) includes the second moments of variables combined such as the
correlation of particle positions and velocities,
$\overline{X_{\text{p}}^{\prime}U_{\text{p}}^{\prime}}$, which have not been
previously reported in the literature for this particular model. For
completeness, we show in Figure 8b the rest of the second moments of the
particle phase which are related to the standard deviation of the particle
position and velocity by
$\overline{X_{\text{p}}^{\prime}U_{\text{p}}^{\prime}}=\sigma_{X}\sigma_{U}$,
$\overline{\Xi^{\prime}X_{\text{p}}^{\prime}}=\sigma_{\Xi}\sigma_{X}$ and
$\overline{\Xi^{\prime}U_{\text{p}}^{\prime}}=\sigma_{\Xi}\sigma_{U}$.

The general system (63), may be used for any other initial conditions or
correlations of the fluctuations by selecting $\overline{\Xi}$ and
$\overline{{\Xi^{\prime}}^{2}}$ differently. For example, if the particle
velocity is initialized with a Maxwellian distribution $\mathcal{N}(0,D)$,
such that $\overline{{U_{\text{p}}^{\prime}}^{2}_{0}}=D$, the system (63)
reproduces the classic ballistic–diffusive result of Taylor [73] where
$\overline{{X_{\text{p}}^{\prime}}^{2}}=2\tau_{\text{p}}^{2}D\left(t-1+e^{-t}\right)$.
This is also for $\overline{\Xi}=0$ and $\overline{{\Xi^{\prime}}^{2}}=1$.

Higher moments are also found with this moment model. Using the
characteristics equations (69a) and applying the LOTUS, one may find the
following general expressions for the $n$-th central moments

$\displaystyle\overline{{X_{\text{p}}^{\prime}}^{n}}$
$\displaystyle=\overline{{\Xi^{\prime}}^{n}}\left[\tau_{\text{p}}^{2}D\left(2t-3+4e^{-t}+e^{-2t}\right)\right]^{n/2},$
(65a) $\displaystyle\overline{{U_{\text{p}}^{\prime}}^{n}}$
$\displaystyle=\overline{{\Xi^{\prime}}^{n}}\left[D\left(1-e^{-2t}\right)\right]^{n/2}.$
(65b)

(a)

(b)

Figure 8: Moments of the particle position and velocities according to the
solutions (64). In (a) solutions of the particle position and velocity
variance which coincides with solutions reported in Ref.[41] and (b)
additional moments correlating particle positions with velocities and the
random coefficient with particle variables. The cases correspond to particle
time constants $\tau_{\text{p}}=[1,\ 10,\ 100]$ and $(\tau_{p}D)^{-1}=5$.

##  Appendix C Derivation of analytical solutions

The method of characteristics, i.e., the mapping of the fixed (Eulerian)
coordinate system in the phase space, $x_{\text{p}}$, onto the moving
(Lagrangian) coordinate system in the phase space, $\hat{x}_{\text{p}}(t)$,
transforms the PDE (27) into a system of ODEs

$\displaystyle\frac{\text{d}\hat{x}_{\text{p}}}{\text{d}t}=$
$\displaystyle\;u_{\text{p}}+\xi\sqrt{\frac{D}{2t}},\qquad\hat{x}_{\text{p}}(0)=\hat{x}_{\text{p}}^{0};$
(66a) $\displaystyle\frac{\text{d}f_{X\Xi}}{\text{d}t}=$
$\displaystyle\;0,\qquad
f_{X\Xi}(\hat{x}_{\text{p}}(0),\xi;0)=f_{\Xi}(\xi)\delta(\hat{x}_{\text{p}}^{0}).$
(66b)

The solution to the characteristics (66) is

$\displaystyle\hat{x}_{\text{p}}$
$\displaystyle=\hat{x}_{\text{p}}^{0}+u_{\text{p}}t+\xi\sqrt{2Dt},$ (67a)
$\displaystyle f_{X\Xi}$ $\displaystyle=f_{X\Xi}^{0}.$ (67b)

After substituting $\hat{x}_{\text{p}}^{0}$ from (67a) we arrive to (28).

The characteristics of the Liouville equation (37) are given by

$\displaystyle\frac{\text{d}\hat{x}_{\text{p}}}{\text{d}t}$
$\displaystyle=\tau_{\text{p}}\hat{u}_{\text{p}}+\xi\varphi,\quad\hat{x}_{\text{p}}(0)=\hat{x}_{\text{p}}^{0};$
(68a) $\displaystyle\frac{\text{d}\hat{u}_{\text{p}}}{\text{d}t}$
$\displaystyle=-\hat{u}_{\text{p}}+\xi\phi,\quad\hat{u}_{\text{p}}(0)=\hat{u}_{\text{p}}^{0};$
(68b) and
$\displaystyle\begin{split}\frac{\text{d}{f}_{XU\Xi}}{\text{d}t}&={f}_{XU\Xi},\\\
f_{XU\Xi}(&\hat{x}_{\text{p}}(0),\hat{u}_{\text{p}}(0),\xi;0)=\delta(\hat{x}_{\text{p}}^{0})\delta(\hat{u}_{\text{p}}^{0}-v_{0})f_{\Xi}(\xi).\end{split}$
(68c)

The joint PDF changes along time because the non-conservative version of the
equation (37) is non-homogeneous. This occurs because the third term in (37)
has an explicit dependency with the particle velocity which leads to a non
zero right hand side in the ODE for the joint PDF (68c). The solution to (68)
is

$\displaystyle\hat{x}_{\text{p}}(t)$
$\displaystyle=\hat{x}_{\text{p}}^{0}+\tau_{\text{p}}\hat{u}_{\text{p}}^{0}(1-\text{e}^{-t})+\xi\sigma_{X},\quad\hat{u}_{\text{p}}(t)=\hat{u}_{\text{p}}^{0}\text{e}^{-t}+\xi\sigma_{U},$
(69a) $\displaystyle f_{XU\Xi}$ $\displaystyle=f_{XU\Xi}^{0}\text{e}^{t}.$
(69b)

where $\sigma_{X}$ and $\sigma_{U}$ are defined by (35). They are derived with
the moment model (55). The general expression for the joint PDF
$f_{XU\Xi}(x_{\text{p}},u_{\text{p}},\xi;t)$ is found by obtaining
$\hat{x}_{\text{p}}^{0}$ and $\hat{u}_{\text{p}}^{0}$ from the characteristics
(69a) and substituting them into (69b):

$\displaystyle
f_{XU\Xi}=e^{t}f_{\Xi}(\xi)\delta(\hat{x}_{\text{p}}^{0})\delta(\hat{u}_{\text{p}}^{0}-v_{0}),$
(70)

which leads to (38).

The characteristics of (45) are

$\displaystyle\frac{\text{d}\hat{u}_{\text{p}}}{\text{d}t}$
$\displaystyle=-\frac{1}{\tau_{\text{p}}}\hat{u}_{\text{p}}+\xi\phi,$ (71)
$\displaystyle\frac{\text{d}{f}_{U\Xi}}{\text{d}t}$
$\displaystyle=\frac{1}{\tau_{\text{p}}}{f}_{U\Xi},$ (72)

that one can be solve analytically, obtaining

$\displaystyle\hat{u}_{\text{p}}$
$\displaystyle=\hat{u}_{\text{p}}^{0}\text{e}^{-t/\tau_{\text{p}}}+\xi(1-\text{e}^{-\mathcal{C}_{1}t/\tau_{\text{p}}})^{\mathcal{C}_{2}},$
(73a) $\displaystyle f_{U\Xi}$
$\displaystyle=f_{U\Xi}^{0}\,\text{e}^{t/\tau_{\text{p}}},$ (73b)

which leads to (46).

## References

  * Buyevich [1971] Y. A. Buyevich.  Statistical hydromechanics of disperse systems Part 1. Physical background and general equations.  _Journal of Fluid Mechanics_ , 49(3):489–507, 1971.
  * Buyevich [1972a] Y. A. Buyevich.  Statistical hydromechanics of disperse systems. Part 2. Solution of the kinetic equation for suspended particles.  _Journal of Fluid Mechanics_ , 52(2):345–355, 1972a.
  * Buyevich [1972b] Y. A. Buyevich.  Statistical hydromechanics of disperse systems. Part 3. Pseudo-turbulent structure of homogeneous suspensions.  _Journal of Fluid Mechanics_ , 56(2):313–336, 1972b.
  * Reeks [1980] M. W. Reeks.  Eulerian direct interaction applied to the statistical motion of particles in a turbulent fluid.  _Journal of Fluid Mechanics_ , 97(3):569–590, 1980.
  * Reeks [1983] M. W. Reeks.  The transport of discrete particles in inhomogeneous turbulence.  _Journal of Aerosol Science_ , 14(6):729–739, 1983.
  * Reeks [1991] M. W. Reeks.  On a kinetic equation for the transport of particles in turbulent flows.  _Physics of Fluids A: Fluid Dynamics_ , 3(3):446–456, 1991.
  * Reeks [1992] M. W. Reeks.  On the continuum equations for dispersed particles in nonuniform flows.  _Physics of Fluids A: Fluid Dynamics_ , 4(6):1290–1303, 1992.
  * Swailes and Darbyshire [1997] D. C. Swailes and K. F. F. Darbyshire.  A generalized Fokker-Planck equation for particle transport in random media.  _Physica A: Statistical Mechanics and its Applications_ , 242(1-2):38–48, 1997.
  * Garzó and Dufty [1999] V. Garzó and J. W. Dufty.  Dense fluid transport for inelastic hard spheres.  _Physical Review E_ , 59(5):5895, 1999.
  * Garzó et al. [2012] V. Garzó, S. Tenneti, S. Subramaniam, and C. M. Hrenya.  Enskog kinetic theory for monodisperse gas–solid flows.  _Journal of Fluid Mechanics_ , 712:129–168, 2012.
  * Bragg et al. [2012] A. Bragg, D. C. Swailes, and R. Skartlien.  Drift-free kinetic equations for turbulent dispersion.  _Physical Review E_ , 86(5):056306, 2012.
  * Swailes et al. [1998] D. C. Swailes, Y. A. Sergeev, and A. Parker.  Chapman-Enskog closure approximation in the kinetic theory of dilute turbulent gas-particulate suspensions.  _Physica A: Statistical Mechanics and its Applications_ , 254(3-4):517–547, 1998.
  * Zaichik and Alipchenkov [2007] L. I. Zaichik and V. M. Alipchenkov.  Refinement of the probability density function model for preferential concentration of aerosol particles in isotropic turbulence.  _Physics of Fluids_ , 19(11):113308, 2007.
  * Alipchenkov and Zaichik [2005] V. M. Alipchenkov and L. I. Zaichik.  Dispersion and clustering of bidisperse particles in isotropic turbulence.  _Fluid Dynamics_ , 40(1):83–94, 2005.
  * Pandya and Mashayek [2001] R. V. R. Pandya and F. Mashayek.  Probability density function modeling of evaporating droplets dispersed in isotropic turbulence.  _AIAA Journal_ , 39(10):1909–1915, 2001.
  * Pandya and Mashayek [2002] R. V. R. Pandya and F. Mashayek.  Turbulent thermal diffusion and barodiffusion of passive scalar and dispersed phase of particles in turbulent flows.  _Physical Review Letters_ , 88(4):044501, 2002.
  * Pandya and Mashayek [2003] R. V. R. Pandya and F. Mashayek.  Non-isothermal dispersed phase of particles in turbulent flow.  _Journal of Fluid Mechanics_ , 475:205–245, 2003.
  * Reeks [2021] M. W. Reeks.  The development and application of a kinetic theory for modeling dispersed particle flows.  _Journal of Fluids Engineering_ , 143(8), 2021.
  * Marchisio et al. [2002a] D. Marchisio, A. Barresi, G. Baldi, and R. O. Fox.  Comparison between the Classes Method and the Quadrature Method of Moments for multiphase systems.  In _Proceedings 8th International Conference on Multiphase Flows in Industrial Plants_ , pages 283–300. ANIMP Servizi srl, 2002a.
  * Marchisio et al. [2002b] D. L. Marchisio, J. Pikturna, L. Wang, R. D. Vigil, and R. O. Fox.  Quadratic method of moments for population balances in CFD applications: Comparison with experimental data.  _Chem. Engg. Trans_ , 1:305–310, 2002b.
  * Marchisio et al. [2003a] D. L. Marchisio, R. D. Vigil, and R. O. Fox.  Quadrature method of moments for aggregation-breakage processes.  _Journal of Colloid and Interface Science_ , 258(2):322–334, 2003a.
  * Marchisio et al. [2003b] D. L. Marchisio, J. T. Pikturna, R. O. Fox, R. D. Vigil, and A. A. Barresi.  Quadrature method of moments for population-balance equations.  _AIChE Journal_ , 49(5):1266–1276, 2003b.
  * Heylmun et al. [2019] J. C. Heylmun, B. Kong, A. Passalacqua, and R. O. Fox.  A quadrature-based moment method for polydisperse bubbly flows.  _Computer Physics Communications_ , 244:187–204, 2019.
  * Bryngelson et al. [2020] S. H. Bryngelson, A. Charalampopoulos, T. P. Sapsis, and T. Colonius.  A Gaussian moment method and its augmentation via LSTM recurrent neural networks for the statistics of cavitating bubble populations.  _International Journal of Multiphase Flow_ , 127:103262, 2020.
  * Charalampopoulos et al. [2022] A. Charalampopoulos, S. H. Bryngelson, T. Colonius, and T. P. Sapsis.  Hybrid quadrature moment method for accurate and stable representation of non-Gaussian processes applied to bubble dynamics.  _Philosophical Transactions of the Royal Society A_ , 380(2229):20210209, 2022.
  * Bryngelson et al. [2023] S. H. Bryngelson, R. O. Fox, and T. Colonius.  Conditional moment methods for polydisperse cavitating flows.  _Journal of Computational Physics_ , page 111917, 2023.
  * Pope [1985] S. B. Pope.  PDF methods for turbulent reactive flows.  _Progress in Energy and Combustion Science_ , 11(2):119–192, 1985.
  * Haworth and Pope [1986] D. C. Haworth and S. B. Pope.  A generalized Langevin model for turbulent flows.  _Physics of Fluids_ , 29(2):387–405, 1986.
  * Iliopoulos et al. [2003] I. Iliopoulos, Y. Mito, and T. J. Hanratty.  A stochastic model for solid particle dispersion in a nonhomogeneous turbulent field.  _International Journal of Multiphase Flow_ , 29(3):375–394, 2003.
  * Gao and Mashayek [2004a] Z. Gao and F. Mashayek.  Stochastic modeling of evaporating droplets polydispered in turbulent flows.  _International Journal of Heat and Mass Transfer_ , 47(19-20):4339–4348, 2004a.
  * Gao and Mashayek [2004b] Z. Gao and F. Mashayek.  Stochastic model for nonisothermal droplet-laden turbulent flows.  _AIAA Journal_ , 42(2):255–260, 2004b.
  * Gao and Mashayek [2004c] Z. Gao and F. Mashayek.  A stochastic model for gravity effects in particle-laden turbulent flows.  _J. Fluids Eng._ , 126(4):620–625, 2004c.
  * Shotorban and Mashayek [2005] B. Shotorban and F. Mashayek.  Modeling subgrid-scale effects on particles by approximate deconvolution.  _Physics of Fluids_ , 17(8):081701, 2005.
  * Shotorban and Mashayek [2006a] B. Shotorban and F. Mashayek.  A stochastic model for particle motion in large-eddy simulation.  _Journal of Turbulence_ , 7:N18, 2006a.
  * Shotorban and Mashayek [2006b] B. Shotorban and F. Mashayek.  On stochastic modeling of heavy particle dispersion in large-eddy simulation of two-phase turbulent flow.  In _IUTAM Symposium on Computational Approaches to Multiphase Flow: Proceedings of an IUTAM Symposium held at Argonne National Laboratory, October 4–7, 2004_ , pages 373–380. Springer, 2006b.
  * Sengupta et al. [2009] K. Sengupta, B. Shotorban, G. B. Jacobs, and F. Mashayek.  Spectral-based simulations of particle-laden turbulent flows.  _International Journal of Multiphase Flow_ , 35(9):811–826, 2009.
  * Pozorski and Apte [2009] J. Pozorski and S. V. Apte.  Filtered particle tracking in isotropic turbulence and stochastic modeling of subgrid-scale dispersion.  _International Journal of Multiphase Flow_ , 35(2):118–128, 2009.
  * Pai and Subramaniam [2012] M. G. Pai and S. Subramaniam.  Two-way coupled stochastic model for dispersion of inertial particles in turbulence.  _Journal of Fluid Mechanics_ , 700:29–62, 2012.
  * Tenneti et al. [2016] S. Tenneti, M. Mehrabadi, and S. Subramaniam.  Stochastic Lagrangian model for hydrodynamic acceleration of inertial particles in gas–solid suspensions.  _Journal of Fluid Mechanics_ , 788:695–729, 2016.
  * Esteghamatian et al. [2018] A. Esteghamatian, F. Euzenat, A. H., M. Lance, and A. Wachs.  A stochastic formulation for the drag force based on multiscale numerical simulation of fluidized beds.  _International Journal of Multiphase Flow_ , 99:363–382, 2018.
  * Lattanzi et al. [2020] A. M. Lattanzi, V. Tavanashad, S. Subramaniam, and J. Capecelatro.  Stochastic models for capturing dispersion in particle-laden flows.  _Journal of Fluid Mechanics_ , 903, 2020.
  * Knorps and Pozorski [2021] M. Knorps and J. Pozorski.  Stochastic modeling for subgrid–scale particle dispersion in large–eddy simulation of inhomogeneous turbulence.  _Physics of Fluids_ , 33(4):043323, 2021.
  * Lattanzi et al. [2022a] A. M. Lattanzi, V. Tavanashad, S. Subramaniam, and J. Capecelatro.  Stochastic model for the hydrodynamic force in Euler–Lagrange simulations of particle-laden flows.  _Physical Review Fluids_ , 7(1):014301, 2022a.
  * Friedrich et al. [2022] J. Friedrich, B. Viggiano, M. Bourgoin, R. B. Cal, and L. Chevillard.  Single inertial particle statistics in turbulent flows from Lagrangian velocity models.  _Physical Review Fluids_ , 7(1):014303, 2022.
  * Pietrzyk et al. [2022] K. Pietrzyk, J. A. K. Horwitz, F. M. Najjar, and R. W. Minich.  On analysis and stochastic modeling of the particle kinetic energy equation in particle-laden isotropic turbulent flows.  _Physics of Fluids_ , 34(1):013316, 2022.
  * Sheikhi et al. [2007] M. R. H. Sheikhi, P. Givi, and S. B. Pope.  Velocity-scalar filtered mass density function for large eddy simulation of turbulent reacting flows.  _Physics of fluids_ , 19(9), 2007.
  * Gheorghiu et al. [2003] S Gheorghiu, JR Van Ommen, and M-O Coppens.  Power-law distribution of pressure fluctuations in multiphase flow.  _Physical Review E_ , 67(4):041305, 2003.
  * Biferale et al. [2012] L Biferale, P Perlekar, M Sbragaglia, and FJPRL Toschi.  Convection in multiphase fluid flows using lattice boltzmann methods.  _Physical Review Letters_ , 108(10):104502, 2012.
  * Capecelatro et al. [2015] J. Capecelatro, O. Desjardins, and R. O. Fox.  On fluid–particle dynamics in fully developed cluster-induced turbulence.  _Journal of Fluid Mechanics_ , 780:578–635, 2015.
  * Forgues et al. [2019] F. Forgues, L. Ivan, A. Trottier, and J. G. McDonald.  A gaussian moment method for polydisperse multiphase flow modelling.  _Journal of Computational Physics_ , 398:108839, 2019.
  * Domínguez-Vázquez et al. [2021] D. Domínguez-Vázquez, G. B. Jacobs, and D. M. Tartakovsky.  Lagrangian models of particle-laden flows with stochastic forcing: Monte Carlo, moment equations, and method of distributions analyses.  _Physics of Fluids_ , 33(3):033326, 2021.
  * Venturi et al. [2013] D. Venturi, D. M. Tartakovsky, A. M. Tartakovsky, and G. E. Karniadakis.  Exact PDF equations and closure approximations for advective-reactive transport.  _Journal of Computational Physics_ , 243:323–343, 2013.
  * Cho et al. [2014] H. Cho, D. Venturi, and G. E. Karniadakis.  Statistical analysis and simulation of random shocks in stochastic Burgers equation.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 470(2171):20140080, 2014.
  * Domínguez-Vázquez and Jacobs [2024] D. Domínguez-Vázquez and G. B. Jacobs.  SPARSE–R: A point-cloud tracer with random forcing.  _International Journal of Multiphase Flow_ , 170:104653, 2024.  ISSN 0301-9322.  doi: https://doi.org/10.1016/j.ijmultiphaseflow.2023.104653.
  * Dominguez-Vazquez et al. [2023] D. Dominguez-Vazquez, S. A. Castiblanco-Ballesteros, G. B. Jacobs, and D. M. Tartakovsky.  Lagrangian Liouville models of multiphase flows with randomly forced inertial particles.  _arXiv preprint arXiv:2312.07713_ , 2023.
  * Domínguez-Vázquez et al. [2022] D. Domínguez-Vázquez, Q. Wang, and G. B. Jacobs.  Adjoint-based particle forcing reconstruction and uncertainty quantification.  _arXiv preprint arXiv:2211.10849_ , 2022.
  * Koch and Sangani [1999] D. L. Koch and A. S. Sangani.  Particle pressure and marginal stability limits for a homogeneous monodisperse gas-fluidized bed: kinetic theory and numerical simulations.  _Journal of Fluid Mechanics_ , 400:229–263, 1999.
  * Lattanzi et al. [2022b] A. M. Lattanzi, V. Tavanashad, S. Subramaniam, and J. Capecelatro.  Fluid-mediated sources of granular temperature at finite Reynolds numbers.  _Journal of Fluid Mechanics_ , 942:A7, 2022b.
  * Schneider et al. [2011] D. Schneider, J. Fuhrmann, W. Reich, and G. Scheuermann.  A variance based FTLE-like method for unsteady uncertain vector fields.  In _Topological Methods in Data Analysis and Visualization II: Theory, Algorithms, and Applications_ , pages 255–268. Springer, 2011.
  * Risken and Frank [1996] H. Risken and T. Frank.  _The Fokker-Planck Equation: Methods of Solution and Applications_.  Springer, Berlin, 1996.
  * Tartakovsky and Gremaud [2016] D. M. Tartakovsky and P. A. Gremaud.  Method of distributions for uncertainty quantification.  In R. Ghanem, D. Higdon, and H. Owhadi, editors, _Handbook of Uncertainty Quantification_. Springer, Cham., 2016.  doi: 10.1007/978-3-319-11259-6_27-1.
  * Moyal [1949] J. E. Moyal.  Stochastic processes and statistical physics.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , 11(2):150–210, 1949.
  * Soong and Bogdanoff [1974] T. T. Soong and J. L. Bogdanoff.  _Random Differential Equations in Science and Engineering_.  Academic Press, New York, 1974.
  * Mashayek and Pandya [2003] F. Mashayek and R. V. R. Pandya.  Analytical description of particle/droplet-laden turbulent flows.  _Progress in Energy and Combustion Science_ , 29(4):329–378, 2003.
  * Flury [2013] B. Flury.  _A first course in multivariate statistics_.  Springer Science & Business Media, 2013.
  * Pope [2000] S. B. Pope.  _Turbulent flows_.  Cambridge University Press, 2000.
  * Gardiner [2009] C. Gardiner.  _Stochastic methods_.  Springer, Berlin, 2009.
  * Sen et al. [2018] O. Sen, N. J. Gaul, K. K. Choi, G. B. Jacobs, and H. S. Udaykumar.  Evaluation of multifidelity surrogate modeling techniques to construct closure laws for drag in shock–particle interactions.  _Journal of Computational Physics_ , 371:434–451, 2018.
  * Das et al. [2018a] P. Das, O. Sen, K.K . Choi, G. Jacobs, and H. S. Udaykumar.  Strategies for efficient machine learning of surrogate drag models from three-dimensional mesoscale computations of shocked particulate flows.  _International journal of multiphase flow_ , 108:51–68, 2018a.
  * Das et al. [2018b] P. Das, O. Sen, G. Jacobs, and H. S. Udaykumar.  Metamodels for interphase heat transfer from mesoscale simulations of shock–cylinder interactions.  _AIAA journal_ , 56(10):3975–3987, 2018b.
  * Davis et al. [2017] S. L. Davis, G. B. Jacobs, O. Sen, and H. S. Udaykumar.  SPARSE—A subgrid particle averaged reynolds stress equivalent model: testing with a priori closure.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 473(2199):20160769, 2017.
  * Domínguez-Vázquez et al. [2023] D. Domínguez-Vázquez, B. F. Klose, and G. B. Jacobs.  Closed SPARSE—A predictive particle cloud tracer.  _International Journal of Multiphase Flow_ , 161:104375, 2023.
  * Taylor [1922] G. I. Taylor.  Diffusion by continuous movements.  _Proceedings of the London Mathematical Society_ , 2(1):196–212, 1922.