# Towards a lattice-Fokker-Planck-Boltzmann method for simulating thermal
noise in compressible phase-changing fluids

K.J. Petersen<EMAIL_ADDRESS>University of British Columbia, Okanagan,
Canada  J.R. Brinkerhoff  University of British Columbia, Okanagan, Canada

###### Abstract

Molecular fluctuations are known to affect the macroscopic and spatio-temporal
evolution of a host of physical phenomena central to the study of biological
systems, turbulence, and reactive mixtures, among others. In phase-changing
fluids metastability and nucleation rates of embryos are known to be non-
trivially affected by thermal noise stemming from molecules' random velocity
fluctuations, which ultimately determine the long-term growth, morphology, and
decay of macroscopic bubbles in cavitation and boiling. We herein present the
mathematical groundwork for a lattice-based solution of the combined
stochastic Fokker-Planck and Boltzmann equations that by proxy solve the
compressible fluctuating Navier-Stokes-Fourier equations. We present the
derivation of the kinetic lattice-Fokker-Planck-Boltzmann equations
facilitated by means of Gauss-Hermite quadrature, and show by multi-scale
asymptotic analysis that the non-equilibrium dynamics in velocity space
inherent to the Fokker-Planck equation manifest as fluctuating stresses
analogous to the conventional viscous-stress tensor. The key implication of
this analysis is a direct correlation between the thermal white noise of the
Fokker-Planck equation and the effective viscosities $\mu,\zeta$ that in
concert with a cubic equation of state enable simulations of phase-changing
fluids perturbed by stochastic thermal fluctuations.

Usage

Our theoretical development is inspired by the lattice-Fokker-Planck
development in [Moroni  _et al._ Phys. Rev. E 73, 066707] and employs
Particles-on-Demand for kinetic theory presented in [Reyhanian  _et al._ Phys.
Rev. E 102, 020103(R)] to mitigate numerical instability issues associated
with high compressibility levels and temperature gradients.

lattice-Fokker-Planck-Boltzmann method (LFPBM), Particles-on-Demand (PonD),
Navier-Stokes-Fourier (NSF) equations, kinetic theory of liquids, phase
change, thermal noise, metastability

††preprint: APS/123-QED

##  I Introduction

###  I.1 Fluctuations, nucleation, and phase transitions

Nucleation is responsible for initiating homogeneous as well as heterogeneous
phase transitions in fluid systems, and occurs in a myriad of natural and
engineered processes. Phase transitions are commonly described in literature
to be either pressure-driven along isotherms (cavitating) or thermally-driven
along isobars (boiling) but in reality evolve in a blended manner
simultaneously along both the $p$ and $T$ dimensions [3]. Cavitation dominates
in highly advective flows with large pressure gradients causing local
pressures to fall below the saturation pressure, whereas boiling occurs in
diabetic flows where local temperatures exceed the saturation temperature. The
relevance of either of the modes is exemplified by technological applications
such as medical drug delivery and acoustic therapies facilitated by cavitation
[4], turbomachinery erosion and efficiency loss from cavitation [5], pool
boiling in cryogenic spills [3], as well as natural occurrences in skeletal-
joint cavitation [6], cavitation in the xylem of drought-stricken trees [7],
cavitation as a predation tool in mantis shrimp [8], to mention a few exotic
cases.

Inherent to cavitation and boiling is the property that phase change may occur
either homogeneously or heterogeneously, and that the thermodynamics can be
metastable. In heterogeneous phase change, the presence of nucleation
sites—e.g. on the surface of suspended solid particles and in undissolved
gases in cavitation, or in microscale crevices in the surface topology of a
solid boiling substrate—can non-trivially alter the effective liquid tensile
strength by an order of magnitude in comparison to the homogeneous counterpart
occurring in pure fluids [9, 10, 11], thus increasing the probability of
nucleation events occurring.

Commonly, phase transitions are thought to occur instantaneously when the
pressure falls below the saturation state, and conversely when the temperature
exceeds the saturation temperature, when they in reality are relaxed processes
in time. This means that the thermodynamic state of a stretched or superheated
liquid in between its binodal and spinodal curves may remain in its original
(liquid) state for a finite period of time until it is perturbed [12, 13]
where the phase stability is fundamentally dictated by entropy generation and
the entropy-extremum principle [14]. Such perturbations can by exemplification
be pressure discontinuities in shock waves propagating on the macroscales,
and/or thermal noise below and around the microscales where large ensembles of
molecules can produce fluctuations in the density and temperature of a fluid
thence provoking a phase transition [15, 16].

In addition to the above complexity, phase transitions statistically evolve
across scales; new thermodynamic phases originate from the molecular scale and
develops into the macroscale [17, 9], and the size distributions of impurities
span the nm-mm scales. Consequently, the nuclei density (number of nucleation
sites/m3), and the rate at which vapour phases grow during phase change
markedly impacts the thermo-hydrodynamics of many flows including the cloud-
shedding dynamics in cavitating hydrofoils [18], or hypothetically the Taylor-
Helmholtz instability and conjugate heat transfer in pool boiling [3].

Capturing the heterogeneous and metastable nature of phase-transitioning
fluids in numerical simulations is rather difficult, and in computational,
macroscopic simulations of industrial relevance it is frequently assumed that
phase transitions occur homogeneously and instantaneously, as found in our
recent literature reviews [3, 19]. Due to the span of length scales in
heterogeneous phase transitions, resolving nucleation sites in macroscale
simulations is computationally intractable with contemporary computational
resources. Thus, a more viable approach to the problem is to involve sub-grid
modelling of unresolved nuclei that can act as nucleation sites. In cavitation
and boiling simulations such heterogeneity in the form of undissolved
microbubbles and solid surfaces, rigorously coupled with the physics of phase
transformations and capillary effects of surface imperfections, is yet to be
fully represented in computational fluid dynamics. Nevertheless, there do
exist strategies for introducing sub-grid nuclei into simulations. In
cavitation, works such as [20] evolved sub-grid microbubbles with a Rayleigh-
Plesset equation including phase-change mass transfer. In the case of boiling,
nucleation-site density models for pool boiling [21, 22] have proven useful,
but still rely on empirically informed activation models. Otherwise, the
computational fluid dynamics community predominantly rely on lower-fidelity,
empirical models for predicting latent heat transfer in cavitation and boiling
simulations [3], such as those of Kunz  _et al._ [23], Schnerr and Sauer [24],
and Zwart  _et al._ [25], which account for nucleation sites via a constant
nucleation-site volume fraction, but not for non-equilibrium thermodynamic
effects nor the impact of micro-scale thermal fluctuations on nucleation.

Recently, Menzl  _et al._ [9] showed that in homogeneously cavitating water,
the classical nucleation theory of Volmer and Weber [26], Farkas [27], Becker
and Döring [28] combined with a macroscopic, thermal-noise-augmented Rayleigh-
Plesset equation was able to reproduce nucleation rates from molecular-
dynamics simulations demonstrating the non-trivial diffusive effects of
molecular fluctuations in phase transitions. The concept of including thermal-
noise in macroscopic simulations is the fundament of fluctuating-hydrodynamics
theory. In recent works by Magaletti  _et al._ [13], Gallo  _et al._ [29, 17],
Gallo [15] the authors combined the theory with Navier-Stokes simuations of
boiling and cavitation. From a similar perspective as the authors' Gaussian
white-noise enforced viscous-stress and heat-flux components, we aim to
investigate the fluid physics of thermal noise and its effect on phase
transitions across the micro and macroscales.

The Fokker-Planck equation (FPE) by Fokker [30] and Planck [31] has been
suggested to be a suitable model of liquids [32] compared to the popularly
employed Boltzmann equation in lattice-Boltzmann methods, which arguably is a
powerful tool in studying phase-changing flows [19]. This is owed to the FPE
being directly derived from the Langevin equation, which models the velocity
space of particles including their chaotic fluctuations that arise due to the
continuous interaction between the electro-static potentials of neighbouring
molecules and other external forces. The FPE thereby mimics the stochastic
behaviour of the Langevin equation, from the underlying molecular dynamics
that manifests as thermal noise in the FPE. We are interested in investigating
the suitability of the FPE as a versatile model of liquids; as the continuum
density, momentum and total energy are recovered from the moments of the
primitive variables of the FPE, the inherited Langevin noise directly projects
into the macroscopic variables as additional noisy components
$\widetilde{\bm{u}},\widetilde{T},\widetilde{p}$, etc. Moreover, entropy
generated by thermal noise can be directly derived from the FPE [33] and as
entropy maxima dictate phase stability [14], the FPE framework may present an
avenue for rigorously studying the non-equilibrium residence times of fluids
in metastable regimes. The infinitesimal oscillations in pressure and
temperature fields are presumed to enable the perturbation of (respectively)
stretched and superheated metastable fluids, and thereby spontaneous
nucleation.

The FPE has already enjoyed success in investigating the turbulence-energy
cascade and its associated entropy generation [34], modelling rarefied-gas
dynamics [35, 36, 37, 38], ion transport in clays [39], homogeneous and
heterogeneous condensation [40], among others. Thus, in a series of articles,
we wish to assess the potential numerical and theoretical capabilities of the
FPE tailored to nucleation and first-order phase transitions. As such, with
this article we initially seek a lattice-based solution framework of the FPE
for a simple, pure van der Waals fluid based on a cubic equation of state
(EOS). We acknowledge that real noise in fluids may not be white, but
attributed with a coloured spectrum, pending further investigation in a
dedicated literature review. Nevertheless, prior to exploring coloured noise,
we first wish to establish the lattice solution with idealized white noise.
Eventually the work can be extended to accommodate more unique spectra for
various fluids, more sophisticated EOS, as well as multi-component kinetic
models. To understand how nucleation is mediated by thermal noise, section
§I.2 lays out the Langevin foundation of the FPE and explores its benefits as
a kinetic model of liquids [41].

###  I.2 The Fokker-Planck equation

First-order phase transitions occur when a critically large population of
molecules experience sudden large fluctuations in the momenta that manifest a
perturbation in entropy and hence a transition between stable states
associated with the concavity of entropy [14]. Since the 1970's the Fokker-
Planck equation has proven a powerful tool for studying problems that involve
stochastic fluctuations (noise), and how fluctuations affect systems around
critical transition points [42]. In describing the FPE, we start with the 1D
linear Langevin equation for Brownian acceleration of a particle with mass $m$
and speed $\xi$ as,

$\displaystyle\eta=\partial_{t}\xi=-\gamma\xi+\widetilde{\Gamma}(t),$ (1)

in which $\gamma\doteq\alpha/m=1/\tau$ is a friction factor that correlates
with a Stokes-like damping-force coefficient $\alpha$, and a relaxation time
$\tau$ [42]. For simplicity we only included the Langevin equation for a
single velocity component, but it generalizes to a molecular-velocity vector
$\bm{v}=\bigl{[}\xi_{x},\xi_{y},\xi_{z}\bigr{]}^{\dagger}$ in three
dimensions. In the article we use bold-font for vector quantities and for the
viscous stress tensor $\bm{\tau}$. The stochasticity in the equation is
facilitated by the per-unit mass Langevin fluctuation force term
$\widetilde{\Gamma}=\widetilde{F}/m$. Assuming that $\widetilde{\Gamma}(t)$ is
Gaussian and $\delta$-correlated, a set of key properties arise from the
fluctuation-dissipation theorem (FDT); namely that the average over an
ensemble is zero $\bigl{\langle}\widetilde{\Gamma}(t)\bigr{\rangle}=0$, and
that the ensemble-average of a product of two Langevin force-terms at
different times is zero if $t-t^{\prime}$ is greater than the correlation time
$\tau_{c}$, i.e.
$\bigl{\langle}\widetilde{\Gamma}(t)\widetilde{\Gamma}(t^{\prime})\bigr{\rangle}=0,|t-t^{\prime}|\geq\tau_{c}$.
Conversely, when we are interested in timescales smaller than the correlation
time, the ensemble average is Dirac $\delta$-correlated,

$\displaystyle\bigl{\langle}\widetilde{\Gamma}(t)\widetilde{\Gamma}(t^{\prime})\bigr{\rangle}=q\delta(t-t^{\prime}),$
(2)

where the noise-strength can be shown to be $q=2\gamma k_{B}T/m$, where
$k_{B}$ is Boltzmann's constant and $T$ is the absolute temperature. In a
three dimensional case $\widetilde{\Gamma}(t)$ would be a vector quantity with
mutually independent components for each of the spatial dimensions.
Conversely, if the noise components were mutually dependent the strength $q$
would be a matrix quantity. The noise is ``white'' characterized by the
spectral density of the Langevin equation (1) being independent of frequency.
On the other hand, for non-$\delta$-correlated terms, the spectral density is
frequency-dependent and hence the noise is attributed as being ``coloured''.
We here clarify that we are not interested in directly solving the Langevin
equation for a finite number of particles, but instead wish to solve an
equivalent equation of motion for a distribution $f(\bm{v};\bm{x},t)$ of
Langevin particles. With this notation of the arguments that we generally
adopt in the article, we imply that $f$ is inherently dependent on the
molecular velocities being the quadrature nodes of the population, whereas the
spatio-temporal dependence is indirectly owed to the time-integration of the
FPE, as well will detail later. The bridge between the FP and Langevin
equations is the Kramers-Moyal (KM) expansion [43, 44] around a single or
larger set of fluctuating variables $\bm{v}\in\\{\chi\\}$, where we
specifically consider the molecular velocity vector $\bm{v}$ in one and two
spatial dimensions. If the Langevin equation with Gaussian $\delta$-correlated
noise governs $\bm{v}$ it can be shown [42] that the KM expansion,

$\displaystyle\partial_{t}f\bigl{(}\bm{v};\bm{x},t\bigr{)}=\sum_{n=1}^{\infty}\bigl{(}-\partial_{\bm{v}}\bigr{)}^{n}\mathcal{D}^{(n)}f\bigl{(}\bm{v};\bm{x},t\bigr{)},$
(3)

reduces to the FPE, with vanishing coefficients $\mathcal{D}^{(n)}$ for the
higher-order terms with $n\geq 3$. In this case the expansion comprises the
deterministic drift coefficient $\mathcal{D}^{(1)}\ (n=1)$, and the diffusion
coefficients $\mathcal{D}^{(2)}\ (n=2)$ accounting for the fluctuations in the
variable. Truncation to any finite order $n>2$ results in introducing
vanishingly small noise contributions from the Langevin equation, guaranteeing
that the process is statistically continuous except for some non-vanishing
cases [45] where discrete, anomalous events play a significant role [34]. The
drift and diffusion terms read,

$\displaystyle\mathcal{D}^{(1)}$ $\displaystyle=-\gamma\bm{v}-\bm{\eta},$ (4)
$\displaystyle\mathcal{D}^{(2)}$ $\displaystyle=1/2q=\gamma k_{B}T/m,$ (5)

which in concert with (3) for $f(\bm{v};\bm{x},t)$ in position and velocity-
space produces the FP equation [42],

$\displaystyle\bigl{(}\partial_{t}+v_{\alpha}\partial_{\alpha}+\eta_{\alpha}\partial_{v_{\alpha}}\bigr{)}f=\widehat{\Omega}^{\textrm{(FP)}}\circ
f\doteq\gamma\partial_{v_{\alpha}}\bigl{(}v_{\alpha}+v_{T}^{2}\partial_{v_{\alpha}}\bigr{)}f,$
(6)

where $v_{T}=\sqrt{k_{B}Tm^{-1}}$ is the thermal velocity. We further appended
the LHS body-forcing term with acceleration $\bm{F}_{B}/m=\bm{\eta}$, which is
not to be confused with the Langevin acceleration. We employ the notation for
the FP operator $\widehat{\Omega}^{\textrm{(FP)}}\circ f$ imposed on the
populations, and apply the Einstein notation for which repeated subscripts
denote summation over pertinent indices. At this point, we emphasize that the
FP operator effectively models liquid states [32, 41] and so we also need to
consider a hard-sphere binary collision model for gaseous states. To that end,
the Boltzmann operator
$\Omega^{\textrm{(B)}}\doteq\tau_{c}^{-1}\bigl{(}f^{\textrm{(eq)}}-f\bigr{)}$
with the Maxwell-Boltzmann (MB) equilibrium distribution,

$\displaystyle f^{\textrm{(eq)}}=\frac{\rho}{\bigl{(}2\pi
RT\bigr{)}^{D/2}}\exp\left(-\frac{\bigl{(}\bm{v}-\bm{u}\bigr{)}^{2}}{2RT}\right),$
(7)

introduced into the RHS of (6) yields a Fokker-Planck-Boltzmann (FPB) type
equation,

$\displaystyle\bigl{(}\partial_{t}+v_{\alpha}\partial_{\alpha}\bigr{)}f=\add{-\eta_{\alpha}\partial_{v_{\alpha}}f}$
$\displaystyle+\gamma\partial_{v_{\alpha}}(v_{\alpha}+v_{T}^{2}\partial_{v_{\alpha}})f$
$\displaystyle+\frac{1}{\tau_{c}}\bigl{(}f^{\textrm{(eq)}}-f\bigr{)},$ (8)

where the FP and Boltzmann operators share the same MB equilibrium [32], and
the acceleration term has been moved from the advection side. In the
equilibrium (7) $\rho$ is the mass density, $R$ the gas constant, $D$ the
number of spatial dimensions, $\bm{v}$ the molecular velocity, and $\bm{u}$
the continuum velocity. We are interested in achieving a unique, transient
solution of the FPB equation in seven degrees of freedom, i.e.
$(\bm{v},\bm{x},t)$, and its moment equations enabling simulation of large-
scale fluid-dynamics problems with phase transitions. The FP operator is
discretized on the lattice following Moroni  _et al._ [1] to derive kinetic
equations for phase-changing flows of pure van der Waals fluids. The Boltzmann
operator is modelled following the conventional lattice-Boltzmann method with
the Bhatnagar-Gross-Krook (BGK) single-relaxation time collision model. As
such, this paper revolves around the derivation and nature of such a lattice
solution, and is organized as follows: In §II we seek a phase-space
discretized form of (I.2) via a Hermite-series expansion resulting in a set of
explicit kinetic lattice equations with a discretized velocity space in
accordance with conventional $DmQn$ lattice models, exemplified by the $D2Q9$
variant in particular. Then, in §III a Chapman-Enskog multiscale analysis is
carried out to prove that the continuum conservation laws are satisfied by the
lattice equations. A discussion in §IV revolving around thermodynamic
consistency, as well as strengths and weaknesses of the model, is followed by
a conclusion §V giving an outlook on planned applications.

##  II Kinetic equations

###  II.1 Advection

The FPBE (I.2) comprises two processes that need to be treated separately in
prospective simulations. Going forward we attribute relaxation due to the BGK
operator with ``collisions'', and that due to the FP operator with
``interactions''. The LHS advection terms represent the advection of
populations prior to particle collisions and interactions, and it is
responsible for the spatio-temporal discretization of $f$ in
$(\bm{x},t)$-space. It is well known that the stochastic sample path of a
Langevin particle can be predicted by integration of the (exemplary 1D) Itô
stochastic-differential equation (SDE),

$\displaystyle
d\xi(t)=\mathcal{D}^{(1)}\bigl{(}\xi(t),t\bigr{)}dt+\sqrt{\mathcal{D}^{(2)}\bigl{(}t\bigr{)}}dW(t),$
(9)

in which the first term prefixed by $\mathcal{D}^{(1)}$ accounts for the
deterministic effects of the Langevin equation and the second term prefixed by
$\mathcal{D}^{(2)}$ the stochastic dynamics of the Langevin equation as
prescribed by the Weiner process $W(t)$ with the difference
$dW(t)=\widetilde{\Gamma}(t)dt$ [46]. The equivalence between the Itô
formalism and the FPE is drift and diffusion coefficients that occur in the
SDE and (I.2).

As the stochastic sample path is inherently realized by the velocity space
$\bm{v}(\bm{x},t)$ of the FP operator during collisions, we need to account
for the associated non-equilibrium effects that arise in the populations
$f(\bm{v};\bm{x},t)$. We do this by adapting the lattice-discrete velocity
space to the local thermo-hydrodynamics by means of the recently proposed
Particles-on-Demand method in the form [47],

$\displaystyle\bm{v}_{i}=\sqrt{\theta}\bm{c}_{i}+\bm{u},$ (10a) such that the
inertial reference frame is shifted by the continuum velocity
$\bm{u}(\bm{x},t)$ rendering $\bm{c}_{i}$ the constant relative velocities as
given in the conventional LBM, and the discrete velocities are further scaled
by the local thermodynamics through the reduced pressure $\theta(\bm{x},t)$,
that are predicted by the FPB populations. The subscript $i$ corresponds to
each of the $n$ quadrature vectors of the pertinent $DmQn$ model. We adopt the
real-gas formulation, $\displaystyle\theta=\frac{p}{\rho RT_{L}},$ (10b) where
$T_{L}$ is the lattice temperature of a corresponding isothermal, ideal-gas
realization of the LFPBE, represented by a known constant for every $DmQn$
lattice model [48, 49].

The recast velocities (10) enable simulation of high-Mach and thermal flows
owing to the fact that the Hermite-expansion of the equilibrium populations
shared by both the FP and Boltzmann equations become exact and velocity-error
free [47, 50, 49],

$\displaystyle f_{i}^{\textrm{(eq)}}$ $\displaystyle=\rho w_{i},$ (11)
$\displaystyle g_{i}^{\textrm{(eq)}}$ $\displaystyle=\rho
w_{i}\left(2e-D\frac{p}{\rho}+v_{i}^{2}\right),$ (12)

informed by $\rho(\bm{x},t),D$, the quadrature weights $w_{i}$, internal
energy $e(\bm{x},t)$, and real-gas pressure $p(\bm{x},t)$.
$v_{i}^{2}=\bm{v}_{i}\cdot\bm{v}_{i}$ is understood to be the dot product.
Mitigating velocity errors in a prospective simulation with a lattice-FPBE is
especially attractive as previous lattice-solutions [32] of the FPE have
suffered from significantly reduced stability limits with increasing Reynolds
numbers in comparison to conventional lattice-Boltzmann models.

The benefits of the scheme are contrasted by the complexity of the required
solution methods for the considered lattice equations; in the conventional LBM
the discrete velocity sets ensure exact streaming between lattice nodes, but
this is not the case with (10) that is more likely to incur populations
streaming between off-lattice locations and lattice nodes. Thus, establishing
the populations at off-lattice locations demands a predictor-corrector scheme
and a reconstruction procedure on the basis of the adjacent known on-lattice
populations. Moreover, the off-lattice nature of the discrete velocities
requires that spatio-temporal integration of (I.2) is carried out backwards in
time along the characteristic lines. We adopt the second-order accurate
trapezoidal integration resulting in the kinetic lattice equations,

$\displaystyle f_{i}(\bm{x},t)-f_{i}(\bm{x}-\bm{v}_{i}\delta t,t-\delta t)$
$\displaystyle=\omega^{\textrm{(FP)}}\Omega^{(\textrm{FP},\eta)}_{i}(f_{i})$
$\displaystyle+\omega^{\textrm{(B)}}\Omega^{\textrm{(B)}}_{i}(f_{i}),$ (13a)
$\displaystyle g_{i}(\bm{x},t)-g_{i}(\bm{x}-\bm{v}_{i}\delta t,t-\delta t)$
$\displaystyle=\omega^{\textrm{(FP)}}\Omega^{(\textrm{FP},\eta)}_{i}(g_{i})$
$\displaystyle+\omega^{\textrm{(B)}}\Omega^{\textrm{(B)}}_{i}(g_{i}),$ (13b)
where the populations $f$ and $g$ recovering density and total energy,
respectively, are relaxed by the relaxation frequencies
$\omega^{\textrm{(FP)}}$ and $\omega^{\textrm{(B)}}$ that are imposed on the
lattice-discretized FP-interaction $\Omega^{\textrm{(FP)}}_{i}$ and BGK-
collision $\Omega^{\textrm{(B)}}_{i}$ operators.

In the temporal discretization of (13), the streaming operator results in,

$\displaystyle\int_{t-\delta t}^{t}\Omega^{(\textrm{S})}\circ f(t)\ dt$
$\displaystyle=f^{\prime}-f$ $\displaystyle=\frac{\delta
t}{2}\Bigl{(}\Omega\circ f^{\prime}+\Omega\circ f\Bigr{)},$ (14)

where we adopt the notation for the advected population is
$f^{\prime}=f(\bm{x},t)$, and the initial $f=f(\bm{x}-\bm{v}\delta t,t-\delta
t)$. However, the purpose of the integration is to be able to predict the
post-advected state $f^{\prime}$, and thus it is unfeasible to compute
$\Omega\circ f^{\prime}$ unless a transformation is made to avoid the
implicitness of the equation. An explicit form can be procured analogously to
the BGK model with the population remapping,

$\displaystyle\tilde{f}\doteq f-\frac{\delta t}{2}\Omega\circ f,$ (15)

which in the BGK case yields the kinetic equation [51],

$\displaystyle\tilde{f}^{\prime}-\tilde{f}$
$\displaystyle=-\omega^{\textrm{(B)}}\Omega^{\textrm{(B)}}\circ f$
$\displaystyle=-\frac{\tau^{-1}\delta t}{1+\tau^{-1}\delta
t/2}\Omega^{\textrm{(B)}}\circ\tilde{f}.$ (16)

It was previously shown [52] that the trapezoidal integration of the FP
operator yields a similar relaxation frequency, although instead mediated by
the friction factor $\gamma$,

$\displaystyle\omega^{\textrm{(FP)}}=\frac{\gamma\delta t}{1+\gamma\delta
t/2}.$ (17)

In the BGK case, the zeroth-through-second moments coincide for both
$f,\tilde{f}$, but this is not the case for the FP operator in which only the
zeroth moments coincide. As such, sampling to compute the higher-order moments
$\sum_{i}f_{i}\bm{v}_{i},\sum_{i}f_{i}v_{i,\alpha}v_{i,\beta}$ must be done on
the basis of $f_{i}$ (and similarly for $g_{i}$). This is relevant as the
moments are necessary for carrying out the Hermite-series expansion, as we
will present in the following.

###  II.2 Collision and interaction

Following Reyhanian  _et al._ [50], we adopt the exact equilibria in (11,12)
in the BGK collision model. However, we still need to establish a lattice
model for the FP operator, which will be facilitated with a Hermite expansion.
In the following we repeat some of the work already documented in the article
by Moroni  _et al._ [1], but do so for the sake of completeness. However, as
we will see later, our explicit formulations for the operator is different
from that reported by Moroni  _et al._ owing to the characteristic eigenvalue
property of operator modified by the discrete velocity set (10) of PonD. We
adopt a notation for the operators in which circumflex-accented (
$\widehat{\phantom{a}}$ ) operators denote that the operators are to be
considered only as the operational expressions excluding the populations. This
notation is important in the derivation of the expansions. The discretization
of the operators as well as the populations $f,g$ is facilitated by a Hermite-
polynomial tensor $\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})$ basis in
$D$ dimensions. Just as these tensors can be exploited to transform gradient
terms in the continuous FPE (6), the continuous population
$f(\bm{v};\bm{x},t)$ in phase-space can be expanded in terms of its Hermite
coefficients $\mathcal{F}_{\underline{\alpha}}^{(l)}(\bm{x},t)$ as [1]

$\displaystyle
f(\bm{v};\bm{x},t)=w(\bm{v})\sum_{l=0}^{\infty}\frac{1}{v_{T}^{2l}l!}\mathcal{F}_{\underline{\alpha}}^{(l)}\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v}),$
(18)

where the populations are proportional to the Gaussian weight function [1],

$\displaystyle w(\bm{v})=\frac{1}{\bigl{(}2\pi
v_{T}^{2}\bigr{)}^{D/2}}\exp\left[-\frac{(\bm{v}-\bm{u})^{2}}{2v_{T}^{2}}\right].$
(19)

The subscript $\underline{\alpha}=\alpha_{1}\cdots\alpha_{l}$ in the Hermite
coefficients and polynomials contracts the $l$ indices of the $l$-rank Hermite
tensor $\mathcal{H}_{\underline{\alpha}}^{(l)}$. To that end, the zeroth-
through third-rank tensors are [53, 54],

$\displaystyle\mathcal{H}^{(0)}(\bm{v})$ $\displaystyle=1,$ (20a)
$\displaystyle\mathcal{H}_{\alpha}^{(1)}(\bm{v})$ $\displaystyle=v_{\alpha},$
(20b) $\displaystyle\mathcal{H}_{\alpha\beta}^{(2)}(\bm{v})$
$\displaystyle=v_{\alpha}v_{\beta}-v_{T}^{2}\delta_{\alpha\beta},$ (20c)
$\displaystyle\mathcal{H}_{\alpha\beta\gamma}^{(3)}(\bm{v})$
$\displaystyle=v_{\alpha}v_{\beta}v_{\gamma}-v_{T}^{2}\bigl{(}v_{\alpha}\delta_{\beta\gamma}+v_{\beta}\delta_{\alpha\gamma}+v_{\gamma}\delta_{\alpha\beta}\bigr{)},$
(20d)

and further detailed by Moroni  _et al._ [1]. The Hermite coefficients for the
populations are the velocity integrals which occupy $\mathbb{R}^{D}$,

$\displaystyle\mathcal{F}_{\underline{\alpha}}^{(l)}(\bm{x},t)=\int
d\bm{v}~{}f(\bm{v};\bm{x},t)\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v}).$
(21)

To obtain a tenable series expansion, we truncate (18) at an order $K$ over
$l$ ranks and discard those exceeding $K$, such that we approximate
$f(\bm{v};\bm{x},t)$ as [1],

$\displaystyle
f(\bm{v};\bm{x},t)=w(\bm{v})\sum_{l=0}^{K}\frac{1}{v_{T}^{2l}l!}\mathcal{F}_{\underline{\alpha}}^{(l)}\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v}).$
(22)

By the Gauss-Hermite quadratures, we seek to carry out the velocity integrals
at a finite set of velocity quadrature-nodes, or abscissae, such that the
continuous integrals become $Q$ easily-digestible discrete sums over the index
$i=0,...,Q-1$. Crucially, given any polynomial $p(\bm{v})$ with a lower order
than $2K$, the continuous integral of its product with a Gaussian function
becomes the quadrature [1],

$\displaystyle\int
d\bm{v}~{}w(\bm{v})p(\bm{v})=\sum_{i=0}^{Q-1}w_{i}p(\bm{v}_{i}),$ (23)

accompanied by $\bm{v}_{i}\in\mathbb{R}^{D}$ and $w_{i}$ identified as the
quadrature nodes (10) and weights, respectively. In (22), $f/w$ can be
characterized as a maximum $K$-order polynomial, which conveniently can be
exploited to define the familiar discrete populations $f_{i},g_{i}$. Starting
with the conventional $m$-order moments of $f$,

$\displaystyle M^{(m)}_{\underline{\alpha}}\doteq\int
d\bm{v}~{}f(\bm{v};\bm{x},t)v_{\alpha_{1}}\cdots v_{\alpha_{m}},$ (24)

can be rewritten as [1],

$\displaystyle\int d\bm{v}~{}$
$\displaystyle\frac{w(\bm{v})}{w(\bm{v})}f(\bm{v};\bm{x},t)v_{\alpha_{1}}\cdots
v_{\alpha_{m}}$
$\displaystyle=\sum_{i=0}^{Q-1}\frac{w_{i}}{w(\bm{v}_{i})}f_{i}(\bm{v}_{i};\bm{x},t)v_{i,\alpha_{1}}\cdots
v_{i,\alpha_{m}}$
$\displaystyle\doteq\sum_{i=0}^{Q-1}f_{i}v_{i,\alpha_{1}}\cdots
v_{i,\alpha_{m}},$ (25)

where we have arrived at the familiar discrete population $f_{i}$ by the
ratio,

$\displaystyle\frac{f_{i}(\bm{x},t)}{f_{i}(\bm{v}_{i};\bm{x},t)}=\frac{w_{i}}{w(\bm{v}_{i})}.$
(26)

$K\geq 2$ is required to reproduce the correct hydrodynamic behaviour governed
by the Navier-Stokes-Fourier (NSF) equations. The series up to these
$K$-orders require the conventional quadratures of the $f_{i},g_{i}$
populations which represent mass and total energy, respectively [51, 50],

$\displaystyle\rho$ $\displaystyle\doteq\int
d\bm{v}~{}f(\bm{v};\bm{x},t)=\sum_{i=0}^{Q-1}f_{i},$ (27a) $\displaystyle
j_{\alpha}$ $\displaystyle\doteq\int
d\bm{v}~{}f(\bm{v};\bm{x},t)v_{\alpha}=\sum_{i=0}^{Q-1}f_{i}v_{i,\alpha},$
(27b) $\displaystyle P_{\alpha\beta}$ $\displaystyle\doteq\int
d\bm{v}~{}f(\bm{v};\bm{x},t)v_{\alpha}v_{\beta}=\sum_{i=0}^{Q-1}f_{i}v_{i,\alpha}v_{i,\beta},$
(27c) $\displaystyle Q_{\alpha\beta\gamma}$ $\displaystyle\doteq\int
d\bm{v}~{}f(\bm{v};\bm{x},t)v_{\alpha}v_{\beta}v_{\gamma}=\sum_{i=0}^{Q-1}f_{i}v_{i,\alpha}v_{i,\beta}v_{i,\gamma},$
(27d)

and

$\displaystyle\mathcal{E}$ $\displaystyle\doteq\int
d\bm{v}~{}g(\bm{v};\bm{x},t)=\sum_{i=0}^{Q-1}g_{i},$ (28a) $\displaystyle
q_{\alpha}$ $\displaystyle\doteq\int
d\bm{v}~{}g(\bm{v};\bm{x},t)v_{\alpha}=\sum_{i=0}^{Q-1}g_{i}v_{i,\alpha},$
(28b) $\displaystyle R_{\alpha\beta}$ $\displaystyle\doteq\int
d\bm{v}~{}g(\bm{v};\bm{x},t)v_{\alpha}v_{\beta}=\sum_{i=0}^{Q-1}g_{i}v_{i,\alpha}v_{i,\beta},$
(28c)

Here, we emphasize that $\mathcal{E}=2\rho E$, and
$\rho,u_{\alpha}=j_{\alpha}/\rho,E=e+u_{\alpha}u_{\alpha}/2$ are the
macroscopic observables density, continuum velocity, and total energy, which
are of the interest in prospective simulations, and recoverable as in any
other conventional lattice-Boltzmann model. They also participate in the
hydrodynamic limit derived in §III. Furthermore, we can observe that the
Hermite coefficients
$\mathcal{F}_{\underline{\alpha}}^{(l)},\mathcal{G}_{\underline{\alpha}}^{(l)}$
are merely linear combinations of these, and read,

$\displaystyle\mathcal{F}^{(0)}$ $\displaystyle=\rho,$ (29a)
$\displaystyle\mathcal{F}_{\alpha}^{(1)}$ $\displaystyle=j_{\alpha},$ (29b)
$\displaystyle\mathcal{F}_{\alpha\beta}^{(2)}$
$\displaystyle=P_{\alpha\beta}-v_{T}^{2}\rho\delta_{\alpha\beta},$ (29c)
$\displaystyle\mathcal{F}_{\alpha\beta\gamma}^{(3)}$
$\displaystyle=Q_{\alpha\beta\gamma}-v_{T}^{2}\bigl{(}\delta_{\alpha\beta}j_{\gamma}+\delta_{\alpha\gamma}j_{\beta}+\delta_{\beta\gamma}j_{\alpha}\bigr{)},$
(29d)

for $f_{i}$ and,

$\displaystyle\mathcal{G}^{(0)}$ $\displaystyle=\mathcal{E}=2\rho E,$ (30a)
$\displaystyle\mathcal{G}_{\alpha}^{(1)}$ $\displaystyle=q_{\alpha},$ (30b)
$\displaystyle\mathcal{G}_{\alpha\beta}^{(2)}$
$\displaystyle=R_{\alpha\beta}-v_{T}^{2}\bigl{(}2\rho
E\bigr{)}\delta_{\alpha\beta},$ (30c)

for $g_{i}$. In the following we will establish the Hermite series for
$\Omega^{\textrm{(FP)}}_{i}$ and $\Omega^{(\eta)}_{i}$. As these reside in the
exact same subspace as (22), we can utilize the quadrature of
$\mathcal{F}_{\underline{\alpha}}^{(l)}$ and
$\mathcal{G}_{\underline{\alpha}}^{(l)}$ in the series.

####  II.2.1 The lattice Fokker-Planck operator

We are now prepared to apply the above methodology on the operators
$\widehat{\Omega}^{\textrm{(FP)}}=\gamma\partial_{v_{\alpha}}\bigl{(}v_{\alpha}+v_{T}^{2}\partial_{v_{\alpha}}\bigr{)}$
and $\widehat{\Omega}^{(\eta)}=-\eta_{\alpha}\partial_{v_{\alpha}}$ in order
to derive explicit functionals for
$\Omega^{\textrm{(FP)}}_{i}\bigl{(}f_{i};g_{i}\bigr{)}$ and
$\Omega^{(\eta)}_{i}\bigl{(}f_{i};g_{i}\bigr{)}$. We start with the FP
operator. Extending the function $\widehat{\Omega}^{\textrm{(FP)}}\circ f$ by
$w_{i}/w(\bm{v}_{i})$ enables us to write the action of the Hermite expansion
on that function [1],

$\displaystyle\Omega^{\textrm{(FP)}}_{i}=w_{i}\sum_{l=0}^{K}\frac{1}{v_{T}^{2l}l!}\Omega^{\textrm{(FP)}}_{\underline{\alpha}}\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v}_{i}),$
(31a) where the Hermite coefficient for the operator is computed as,
$\displaystyle\Omega^{\textrm{(FP)}}_{\underline{\alpha}}=\int
d\bm{v}~{}\widehat{\Omega}^{\textrm{(FP)}}\circ
f\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v}).$ (31b)

We are interested in finding an operational expression that can be used to
rewrite $\widehat{\Omega}^{\textrm{(FP)}}$ in a more convenient, explicit
form, to which end we need to exploit the eigenvalue property of the
characteristic Gaussian $w(\bm{v})$. The product
$\widehat{\Omega}^{\textrm{(FP)}}\bigl{[}w(\bm{v})\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\bigr{]}$,
which is identifiable in (31) by insertion of (22), can be constructed and
will comprise a set of gradient and divergence terms up to the second order,
which can be simplified by this eigenvalue property. In App. A we use the
property to expand the product, where the former is found to be,

$\displaystyle\partial_{v_{\beta}}w(\bm{v})=-\frac{\sqrt{\theta}c_{\beta}}{v_{T}^{2}}w(\bm{v}),$
(32)

which indeed is an eigenfunction with the eigenvector
$-\sqrt{\theta}c_{\beta}/v_{T}^{2}$. This result accounts for the recast
velocities (10), whereas the corresponding eigenvalue $-c_{\beta}/v_{T}^{2}$
found in the analysis of Moroni  _et al._ [1] considers the conventional
peculiar velocity set $\bm{c}_{i}\in\mathbb{R}^{Q}$ used by all LB models with
exact streaming. As we demonstrate in App. B, the operational expression
derived from
$\widehat{\Omega}^{\textrm{(FP)}}\bigl{[}w(\bm{v})\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\bigr{]}$
can be used to compute the Hermite coefficients
$\Omega^{\textrm{(FP)}}_{\underline{\alpha}}$, melded with the previously
found Hermite coefficients for $f(\bm{v};\bm{x},t)$ and
$\mathcal{F}_{\underline{\alpha}}^{(l)}$. With the eigenvalue property, we
find that the product expands to,

$\displaystyle\widehat{\Omega}^{\textrm{(FP)}}\bigl{[}w(\bm{v})\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\bigr{]}=\gamma\left(1+u_{\beta}\frac{v_{i,\beta}+u_{\beta}}{v_{T}^{2}}-l\right)\times$
$\displaystyle\bigl{[}w(\bm{v})\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\bigr{]}-\gamma\frac{2u_{\beta}}{v_{T}^{2}}\bigl{[}w(\bm{v})\mathcal{H}_{\beta\underline{\alpha}}^{(l+1)}(\bm{v})\bigr{]},$
(33)

whereas Moroni  _et al._ [1] reported the result
$\widehat{\Omega}^{\textrm{(FP)}}\Bigl{[}w(\bm{v})\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\Bigr{]}=-\gamma
l\Bigl{[}w(\bm{v})\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\Bigr{]}$ with
the peculiar velocity set. In (II.2.1) $v_{i,\beta}$ are the $Q$ recast
velocity vectors (10), and $\mathcal{H}_{\beta\underline{\alpha}}^{(l+1)}$ is
the $(l+1)$–rank tensor relative to the current rank being considered for
$\mathcal{H}_{\underline{\alpha}}^{(l)}$. Thus, the local thermo-hydrodynamics
imposed by $\bm{v}_{i}$ directly affect the FP dynamics. As (II.2.1) is fully
explicit and $\mathcal{F}_{\underline{\alpha}}^{(l)}$ (29) are merely linear
combinations of the population moments (24), we have all the ingredients for
computing the Hermite coefficients
$\Omega^{\textrm{(FP)}}_{\underline{\alpha}}$. The first three are computed in
Tab. 1 using the explicit form of the coefficients (B). To compute the
integral in the coefficients $\widehat{\Omega}^{\textrm{(FP)}}$ we have
utilized the Hermite orthonormality relation as is also showcased in App. B.
These coefficients, combined with (31a), yield the series expanded up to
$l\leq K=2$, for both $f_{i},g_{i}$:

$\displaystyle\Omega^{\textrm{(FP)}}_{i}(f_{i})$
$\displaystyle=w_{i}\gamma\Biggl{[}\rho\left(1+u_{\gamma}\frac{v_{i,\gamma}+u_{\gamma}}{v_{T}^{2}}\right)+\frac{v_{i,\alpha}}{v_{T}^{2}}\left(u_{\gamma}\frac{v_{i,\gamma}+u_{\gamma}}{v_{T}^{2}}-2\right)j_{\alpha}$
$\displaystyle+\frac{v_{i,\alpha}v_{i,\beta}-v_{T}^{2}\delta_{\alpha\beta}}{2v_{T}^{4}}\left(\left\\{u_{\gamma}\frac{v_{i,\gamma}+u_{\gamma}}{v_{T}^{2}}-1\right\\}\left\\{P_{\alpha\beta}-v_{T}^{2}\rho\delta_{\alpha\beta}\right\\}-2\left\\{u_{\alpha}j_{\beta}+u_{\beta}j_{\alpha}\right\\}\right)\Biggr{]},$
(34a) $\displaystyle\Omega^{\textrm{(FP)}}_{i}(g_{i})$
$\displaystyle=w_{i}\gamma\Biggl{[}\mathcal{E}\left(1+u_{\gamma}\frac{v_{i,\gamma}+u_{\gamma}}{v_{T}^{2}}\right)+\frac{v_{i,\alpha}}{v_{T}^{2}}\left(u_{\gamma}\frac{v_{i,\gamma}+u_{\gamma}}{v_{T}^{2}}-2\right)q_{\alpha}$
$\displaystyle+\frac{v_{i,\alpha}v_{i,\beta}-v_{T}^{2}\delta_{\alpha\beta}}{2v_{T}^{4}}\left(\left\\{u_{\gamma}\frac{v_{i,\gamma}+u_{\gamma}}{v_{T}^{2}}-1\right\\}\left\\{R_{\alpha\beta}-v_{T}^{2}(2\rho
E)\delta_{\alpha\beta}\right\\}-2\left\\{u_{\alpha}q_{\beta}+u_{\beta}q_{\alpha}\right\\}\right)\Biggr{]}.$
(34b)

A similar analysis can be carried out for $\widehat{\Omega}^{(\eta)}$ via its
corresponding
$\widehat{\Omega}^{(\eta)}\Bigl{[}w(\bm{v})\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\Bigr{]}$
function.

Table 1: Hermite coefficients for the FP and acceleration lattice operators evaluated with (B, B) as well as the Hermite coefficients for the populations $\mathcal{F}_{\underline{\alpha}}^{(l)}(f_{i})$ (29) and $\mathcal{G}_{\underline{\alpha}}^{(l)}(g_{i})$ (30). The operator coefficients are ultimately used to complete the series expansion yielding $\Omega^{\textrm{(FP)}}_{i}$ and $\Omega^{(\eta)}_{i}$. $l\quad$ | $\Omega^{\textrm{(FP)}}_{\underline{\alpha}}(l)\circ\\{f_{i};g_{i}\\}$ | $\Omega^{(\eta)}_{\underline{\alpha}}(l)\circ\\{f_{i};g_{i}\\}$
---|---|---
0 | $\gamma\left[1+u_{\gamma}\frac{v_{i,\gamma}+u_{\gamma}}{v_{T}^{2}}\right]\Bigl{\\{}\mathcal{F}^{(0)};\mathcal{G}^{(0)}\Bigr{\\}}$ | $-\left[\frac{\eta_{\gamma}u_{\gamma}}{v_{T}^{2}}\right]\Bigl{\\{}\mathcal{F}^{(0)};\mathcal{G}^{(0)}\Bigr{\\}}$
1 | $\gamma\left[u_{\gamma}\frac{v_{i,\gamma}+u_{\gamma}}{v_{T}^{2}}-2\right]\Bigl{\\{}\mathcal{F}_{\alpha}^{(1)};\mathcal{G}_{\alpha}^{(1)}\Bigr{\\}}$ | $\left[\eta_{\alpha}\Bigl{\\{}\mathcal{F}^{(0)};\mathcal{G}^{(0)}\Bigr{\\}}-\frac{\eta_{\gamma}u_{\gamma}}{v_{T}^{2}}\right]\Bigl{\\{}\mathcal{F}_{\alpha}^{(1)};\mathcal{G}_{\alpha}^{(1)}\Bigr{\\}}$
2 | $\gamma\Bigl{[}\left(u_{\gamma}\frac{v_{i,\gamma}+u_{\gamma}}{v_{T}^{2}}-1\right)\Bigl{\\{}\mathcal{F}_{\alpha\beta}^{(2)};\mathcal{G}_{\alpha\beta}^{(2)}\Bigr{\\}}\quad$ | $\Bigl{[}\eta_{\alpha}\Bigl{\\{}\mathcal{F}_{\beta}^{(1)};\mathcal{G}_{\beta}^{(1)}\Bigr{\\}}+\eta_{\beta}\Bigl{\\{}\mathcal{F}_{\alpha}^{(1)};\mathcal{G}_{\alpha}^{(1)}\Bigr{\\}}-\frac{\eta_{\gamma}u_{\gamma}}{v_{T}^{2}}\Bigl{\\{}\mathcal{F}_{\alpha\beta}^{(2)};\mathcal{G}_{\alpha\beta}^{(2)}\Bigr{\\}}\Bigr{]}$
| $\qquad-2\Bigl{(}u_{\alpha}\Bigl{\\{}\mathcal{F}_{\beta}^{(1)};\mathcal{G}_{\beta}^{(1)}\Bigr{\\}}+u_{\beta}\Bigl{\\{}\mathcal{F}_{\alpha}^{(1)};\mathcal{G}_{\alpha}^{(1)}\Bigr{\\}}\Bigr{)}\Bigr{]}$ |

####  II.2.2 The lattice acceleration operator

Moroni  _et al._ [1] found the operational expression
$\widehat{\Omega}^{(\eta)}\Bigl{[}w(\bm{v})\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\Bigr{]}=\frac{\eta_{\beta}}{v_{T}^{2}}\Bigl{[}w(\bm{v})\mathcal{H}_{\beta\underline{\alpha}}^{(l+1)}(\bm{v})\Bigr{]}$
with the peculiar velocity set. With the recast velocities, we instead report
the result,

$\displaystyle\widehat{\Omega}^{(\eta)}\Bigl{[}w(\bm{v})\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\Bigr{]}=\frac{\eta_{\beta}}{v_{T}^{2}}\Bigl{[}w(\bm{v})\bigl{(}\mathcal{H}_{\beta\underline{\alpha}}^{(l+1)}(\bm{v})-u_{\beta}\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\bigr{)}\Bigr{]}.$
(35)

This result can again be exploited in the same Hermite series (31a) as for the
FP operator, recognizing that we are aiming at obtaining the combined lattice
operator
$\Omega^{(\textrm{FP},\eta)}_{i}=\Omega^{\textrm{(FP)}}_{i}+\Omega^{(\eta)}_{i}$,

$\displaystyle\Omega^{(\eta)}_{i}=w_{i}\sum_{l=0}^{K}\frac{1}{v_{T}^{2l}l!}\Omega^{(\eta)}_{\underline{\alpha}}\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v}_{i}),$
(36a) where the Hermite coefficient for the operator is computed as,
$\displaystyle\Omega^{(\eta)}_{\underline{\alpha}}=\int
d\bm{v}~{}\widehat{\Omega}^{(\eta)}[f]\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v}).$
(36b)

Using (35), we get the relation for the Hermite coefficients (B) derived in
App. B, which we use to obtain the first three Hermite coefficients in Tab. 1.
These coefficients, in concert with (36a) yield the full $(l\leq K=2)$–kinetic
operators,

$\displaystyle\Omega^{(\eta)}_{i}(f_{i})$
$\displaystyle=w_{i}\gamma\Biggl{[}-\frac{v_{\gamma}^{(\eta)}u_{\gamma}}{v_{T}^{2}}\rho+\frac{v_{i,\alpha}}{v_{T}^{2}}\left(v_{\alpha}^{(\eta)}\rho-\frac{v_{\gamma}^{(\eta)}u_{\gamma}}{v_{T}^{2}}j_{\alpha}\right)$
$\displaystyle\qquad+\frac{v_{i,\alpha}v_{i,\beta}-v_{T}^{2}\delta_{\alpha\beta}}{2v_{T}^{4}}\biggl{(}\left\\{v_{\alpha}^{(\eta)}j_{\beta}+v_{\beta}^{(\eta)}j_{\alpha}\right\\}$
$\displaystyle\qquad-\frac{v_{\gamma}^{(\eta)}u_{\gamma}}{v_{T}^{2}}\left\\{P_{\alpha\beta}-v_{T}^{2}\rho\delta_{\alpha\beta}\right\\}\biggr{)}\Biggr{]},$
(37a) $\displaystyle\Omega^{(\eta)}_{i}(g_{i})$
$\displaystyle=w_{i}\gamma\Biggl{[}-\frac{v_{\gamma}^{(\eta)}u_{\gamma}}{v_{T}^{2}}\mathcal{E}+\frac{v_{i,\alpha}}{v_{T}^{2}}\left(v_{\alpha}^{(\eta)}\mathcal{E}-\frac{v_{\gamma}^{(\eta)}u_{\gamma}}{v_{T}^{2}}q_{\alpha}\right)$
$\displaystyle\qquad+\frac{v_{i,\alpha}v_{i,\beta}-v_{T}^{2}\delta_{\alpha\beta}}{2v_{T}^{4}}\biggl{(}\left\\{v_{\alpha}^{(\eta)}q_{\beta}+v_{\beta}^{(\eta)}q_{\alpha}\right\\}$
$\displaystyle\qquad-\frac{v_{\gamma}^{(\eta)}u_{\gamma}}{v_{T}^{2}}\left\\{R_{\alpha\beta}-v_{T}^{2}\mathcal{E}\delta_{\alpha\beta}\right\\}\biggr{)}\Biggr{]},$
(37b)

where we adopted $v_{\alpha}^{(\eta)}=\eta_{\alpha}\gamma^{-1}$ as a pseudo
velocity designed by the dimensionality of $[\gamma]=s^{-1}$. Having rewritten
the acceleration term this way, enables us to factorize terms with
corresponding ranks of
$\mathcal{F}_{\underline{\alpha}}^{(l)},\mathcal{G}_{\underline{\alpha}}^{(l)}$
into the same groups as all coefficients now are prefixed by $\gamma$. Having
discretized both operators on the lattice we can combine them in more compact,
convenient forms,

$\displaystyle\Omega^{(\textrm{FP},\eta)}_{i}(f_{i})$
$\displaystyle=\Omega^{(\eta)}_{i}+\Omega^{\textrm{(FP)}}_{i}$
$\displaystyle=w_{i}\gamma\left[\overline{\rho}+\frac{v_{i,\alpha}}{v_{T}^{2}}\overline{\jmath}_{\alpha}+\frac{v_{i,\alpha}v_{i,\beta}-v_{T}^{2}\delta_{\alpha\beta}}{2v_{T}^{4}}\overline{P}_{\alpha\beta}\right],$
(38a) $\displaystyle\Omega^{(\textrm{FP},\eta)}_{i}(g_{i})$
$\displaystyle=\Omega^{(\eta)}_{i}+\Omega^{\textrm{(FP)}}_{i}$
$\displaystyle=w_{i}\gamma\left[\overline{\mathcal{E}}+\frac{v_{i,\alpha}}{v_{T}^{2}}\overline{q}_{\alpha}+\frac{v_{i,\alpha}v_{i,\beta}-v_{T}^{2}\delta_{\alpha\beta}}{2v_{T}^{4}}\overline{R}_{\alpha\beta}\right],$
(38b) where we recasted the $f$-moments yielding density, momentum flux, and
pressure and higher-order tensors into new functionals,
$\displaystyle\quad\overline{\rho}=\left(1+u_{\gamma}\frac{v_{i,\gamma}+u_{\gamma}-v_{\gamma}^{(\eta)}}{v_{T}^{2}}\right)\rho,$
(38c)
$\displaystyle\quad\overline{\jmath}_{\alpha}=j_{\alpha}^{(\eta)}+\left(u_{\gamma}\frac{v_{i,\gamma}+u_{\gamma}-v_{\gamma}^{(\eta)}}{v_{T}^{2}}-2\right)j_{\alpha},$
(38d)
$\displaystyle\quad\overline{P}_{\alpha\beta}=\left(v_{\alpha}^{(\eta)}-2u_{\alpha}\right)j_{\beta}+\left(v_{\beta}^{(\eta)}-2u_{\beta}\right)j_{\alpha}$
$\displaystyle\qquad+\left(u_{\gamma}\frac{v_{i,\gamma}+u_{\gamma}-v_{\gamma}^{(\eta)}}{v_{T}^{2}}-1\right)\bigl{(}P_{\alpha\beta}-v_{T}^{2}\rho\delta_{\alpha\beta}\bigr{)},$
(38e) where $j_{\alpha}^{(\eta)}=\rho v_{\alpha}^{(\eta)}$ is a pseudo
momentum. We have omitted the quadrature-node index notation in the recast
moments and note that their dependence on $v_{i,\gamma}$ is implied going
forward. Evidently, $v_{\gamma}^{(\eta)}$ originates from the redefinition of
$\eta_{\gamma}$. Similarly the recast functionals for the $g_{i}$-population
moments read,
$\displaystyle\quad\overline{\mathcal{E}}=\left(1+u_{\gamma}\frac{v_{i,\gamma}+u_{\gamma}-v_{\gamma}^{(\eta)}}{v_{T}^{2}}\right)\mathcal{E},$
(38f)
$\displaystyle\quad\overline{q}_{\alpha}=u_{\gamma}\frac{v_{i,\gamma}+u_{\gamma}-v_{\gamma}^{(\eta)}}{v_{T}^{2}}q_{\alpha}+\Bigl{(}v_{\alpha}^{(\eta)}-2u_{\alpha}\Bigr{)}\mathcal{E},$
(38g)
$\displaystyle\quad\overline{R}_{\alpha\beta}=\Bigl{(}v_{\alpha}^{(\eta)}-2u_{\alpha}\Bigr{)}q_{\beta}+\Bigl{(}v_{\beta}^{(\eta)}-2u_{\beta}\Bigr{)}q_{\alpha}$
$\displaystyle\qquad+\left(u_{\gamma}\frac{v_{i,\gamma}+u_{\gamma}-v_{\gamma}^{(\eta)}}{v_{T}^{2}}-1\right)\Bigl{(}R_{\alpha\beta}-v_{T}^{2}\mathcal{E}\delta_{\alpha\beta}\Bigr{)}.$
(38h)

The lattice operators (38,38) are the final results for the kinetic scheme.
They define the ``interaction'' process and can be implemented after the semi-
Lagrangian advection step to successively simulate the spatio-temporal
evolution of $\rho,\bm{u},E$ and other macroscopic quantities of interest.

###  II.3 Thermodynamics

Phase transformations include non-equilibrium thermodynamics and large
variations in viscosity and density that reside on a multitude of spatio-
temporal scales that historically have been difficult to treat in the
continuum without resorting to limiting assumptions that appear in many
empirical mass-transfer models for multi-component mixtures, as we have
documented in a recent review [3].

Whereas mixture models account for the latent heat transfer between species
$\dot{m}h_{\textrm{lv}}$, our current single-component implementation evolves
$T,p$ and phase transitions under the saturation curve in tandem with the
empirical van der Waals (vdW) equation of state (EOS). Even though the EOS is
empirical its simple non-monotonic topology with a distinct energy minimum
enables nucleation in kinetic theories as was shown in [2] that employ the
same EOS. Moreover, the cubic vdW EOS is the simplest EOS that retains the
ability to predict spinodals [55] and is inherently mechanically unstable in
the spinodal region and metastable in the binodal region. When a metastable
fluid is perturbed the fluid state will transition towards the binodal curve
where the local free energy is minimized [56]. Reyhanian  _et al._ presented
the governing thermodynamic equations, but for conciseness we reproduce them
here for later reference in our multi-scale analysis (§III). We duly note that
the kinetic equations conserve the density (27a), momentum (27b), and total
energy (28a), where the latter predicts the internal energy $e$ via,

$\displaystyle E=e+\frac{u_{\alpha}u_{\alpha}}{2}.$ (39)

Thus, given $e$, the temperature can be computed by reorganizing the EOS,

$\displaystyle e=c_{v}T-a\rho-\frac{p_{0}}{\rho},$ (40)

in which the term $p_{0}/\rho$ is included to retain positive temperatures in
the subcritical region $T/T_{\textrm{cr}}<0.84375$ [50]. We can then proceed
with defining the thermodynamic pressure as,

$\displaystyle p=\frac{\rho RT}{1-b\rho}-a\rho^{2},$ (41)

which contributes the non-ideal energy contribution $T\bigl{(}\partial
p/\partial T\bigr{)}_{v}\partial_{\alpha}u_{\alpha}$ in the Fourier equation,
as we will document later in the Chapman-Enskog multiscale analysis in
§III.2.2. This analysis further necessitates deriving the isobaric specific
heat and speed-of-sound, to arrive at the thermo-hydrodynamic limit. The
former reads,

$\displaystyle c_{p}\doteq\biggl{(}\frac{\partial h}{\partial
T}\biggr{)}_{p},$ (42)

and the real-gas speed-of-sound,

$\displaystyle\varsigma^{2}\doteq\biggl{(}\frac{\partial
p}{\partial\rho}\biggr{)}_{s},$ (43)

and can be rewritten using the cyclic and Maxwell relations [57] as,

$\displaystyle c_{p}$ $\displaystyle=\biggl{(}\frac{\partial e}{\partial
T}\biggr{)}_{p}+p\biggl{(}\frac{\partial v}{\partial T}\biggr{)}_{p},$ (44)
$\displaystyle\varsigma^{2}$ $\displaystyle=\biggl{(}\frac{\partial
p}{\partial\rho}\biggr{)}_{T}+\frac{T}{\rho^{2}c_{v}}\biggl{(}\frac{\partial
p}{\partial T}\biggr{)}_{v}^{2}.$ (45)

##  III Chapman-Enskog multiscale analysis

As it is our goal to utilize the kinetic model (13) in mesoscale simulations
where we recover macroscopic quantities (27, 28) from microscopic populations,
it is crucial to investigate how the equations behave across scales.
Specifically, we conduct a Chapman-Enskog multiscale analysis [51]. Initially,
the post-advection state is approximated by a Taylor series, after which the
modified kinetic equations are expanded in terms of perturbation series of all
the variables. Subsequently, all terms in the equations can be separated by
their representative orders forming a hierarchy of perturbation relations for
the advection and interaction/collision sides of the equations. Then, the
zero-through-second moments can be taken of these relations to form moment
equations, i.e. scale-dependent forms of the continuum conservation laws. This
step is crucial as it comprises the moments stemming from non-equilibrium,
higher-order contributions of the perturbation series that can be chosen to
either be included or excluded with the aim of assessing their mathematical
formalism and physical implications in the continuum, thus creating a direct
avenue between the microscopic and macroscopic scales. Finally, the moment
equations can be recombined by inverting the perturbation series which results
in the Navier-Stokes-Fourier (NSF) equations. By including selected non-
equilibrium populations in the moment equations it is possible to elucidate
the asymptotic behaviour of different parts of the NSF hierarchy, to which end
we initially focus on the shear-stress tensor.

###  III.1 Perturbation series

The mesoscale behaviour, especially that stemming from the fluctuations of the
FPE, evolves at different characteristics length and time scales. Analysis of
the scale-effects as they manifest in the continuum can be treated with
multiscale methods in which we seek to separate derivative terms, populations,
and macroscopic observables into contributions that scale with a ``smallness''
parameter $\epsilon$ where $\epsilon\ll 1$ is associated with the
thermohydrodynamic limit [41]. By exemplification, the viscous effects of FP
may reside on longer time scales whereas shocks propagate on shorter time
scales. Specifically, we decompose our distributions $f_{i},g_{i}$, spatial
and temporal gradients $\partial_{\alpha},\partial_{t}$, and macroscopic
moments into a family of perturbation series around $\epsilon$ raised to
different exponent values. Firstly, the populations expand to,

$\displaystyle f_{i}$ $\displaystyle=f_{i}^{(0)}+\epsilon
f_{i}^{(1)}+\epsilon^{2}f_{i}^{(2)}+\mathcal{O}\bigl{(}\epsilon^{3}\bigr{)},$
(46a) $\displaystyle g_{i}$ $\displaystyle=g_{i}^{(0)}+\epsilon
g_{i}^{(1)}+\epsilon^{2}g_{i}^{(2)}+\mathcal{O}\bigl{(}\epsilon^{3}\bigr{)},$
(46b) where the zeroth-order contributions are the known equilibrium
distributions, such that the higher orders act as unknown, stochastic,
fluctuating components. Propagation and diffusion phenomena are captured by
decomposing the time derivative into two time variables,
$\displaystyle\partial_{t}=\epsilon\partial_{t}^{(1)}+\epsilon^{2}\partial_{t}^{(2)}+\mathcal{O}\bigl{(}\epsilon^{3}\bigr{)}.$
(46c) As for the spatial derivatives, these are evolved on the hydrodynamic
scale that is appropriately resolved as,
$\displaystyle\partial_{\alpha}=\epsilon\partial_{\alpha}^{(1)}+\mathcal{O}\bigl{(}\epsilon^{2}\bigr{)}.$
(46d) Moreover, we are interested in the zeroth-through-second-order
perturbations in all of the macroscopic observables in order to recover the
correct continuum-conservation laws, $\displaystyle\rho$
$\displaystyle=\rho^{(0)}+\epsilon\rho^{(1)}+\epsilon^{2}\rho^{(2)},$ (46e)
$\displaystyle\mathcal{E}$
$\displaystyle=\mathcal{E}^{(0)}+\epsilon\mathcal{E}^{(1)}+\epsilon^{2}\mathcal{E}^{(2)},$
(46f) $\displaystyle j_{\alpha}$ $\displaystyle=j_{\alpha}^{(0)}+\epsilon
j_{\alpha}^{(1)}+\epsilon^{2}j_{\alpha}^{(2)},$ (46g) $\displaystyle
q_{\alpha}$ $\displaystyle=q_{\alpha}^{(0)}+\epsilon
q_{\alpha}^{(1)}+\epsilon^{2}q_{\alpha}^{(2)},$ (46h) $\displaystyle
P_{\alpha\beta}$ $\displaystyle=P_{\alpha\beta}^{(0)}+\epsilon
P_{\alpha\beta}^{(1)}+\epsilon^{2}P_{\alpha\beta}^{(2)},$ (46i) $\displaystyle
R_{\alpha\beta}$ $\displaystyle=R_{\alpha\beta}^{(0)}+\epsilon
R_{\alpha\beta}^{(1)}+\epsilon^{2}R_{\alpha\beta}^{(2)},$ (46j) and similarly
for the linear recombinations of the above,
$\overline{\rho}=\overline{\rho}^{(0)}+\epsilon\overline{\rho}^{(1)}+\epsilon^{2}\overline{\rho}^{(2)}$
and so forth.

We interpret the zeroth orders as ``equilibrium'' contributions that dictate
the phenomena that evolve on the macroscales, whereas increasing orders
represent perturbations on diminishing scales. Consequently, we can imagine
that microscopic fluctuations due to thermal noise may be observable on the
$\epsilon^{1},\epsilon^{2}$ scales, where the coarser fluctuations associated
with $\epsilon^{1}$ are likely more relevant to analyze in the context of
mesoscopic nucleation.

###  III.2 Expansions: Taylor series, perturbations, separation, moments, and
recombination

####  III.2.1 Perturbation-series expansion

We expand the previously derived lattice equations (13) to prove their
consistency in the thermohydrodynamic limit. In the following we adopt the
notation for the pre-advection populations
$f_{i}=f_{i}(\bm{x}-\bm{v}_{i}\delta t,t-\delta t)$ and post-advection
populations $f_{i}^{\prime}=f_{i}(\bm{x},t)$. First we treat the advection
side of the equations, and seek to approximate the advected state
$f_{i}^{\prime}$ by its Taylor expansion around the datum $f_{i}$ where we
exploit the common ansatz that third and higher order terms $n>2$ are very
small and do not significantly affect the macroscopic behaviour [51]. As such,
to find the Navier-Stokes-Fourier equations, we only retain the two lowest
orders in the Taylor-series expansion,

$\displaystyle f_{i}^{\prime}=f_{i}$ $\displaystyle+\delta
t\bigl{(}\partial_{t}+v_{i,\delta}\partial_{\delta}\bigr{)}f_{i}$
$\displaystyle+\frac{\delta
t^{2}}{2}\bigl{(}\partial_{t}+v_{i,\delta}\partial_{\delta}\bigr{)}\bigl{(}\partial_{t}+v_{i,\varepsilon}\partial_{\varepsilon}\bigr{)}f_{i}+\mathcal{O}\bigl{(}\delta
t^{3}\bigr{)},$ (47)

and similarly for $g_{i}$. It follows that the lattice equations (13) can be
rewritten to,

$\displaystyle\delta t\bigl{(}\partial_{t}$
$\displaystyle+v_{i,\delta}\partial_{\delta}\bigr{)}f_{i}+\frac{\delta
t^{2}}{2}\bigl{(}\partial_{t}+v_{i,\delta}\partial_{\delta}\bigr{)}\bigl{(}\partial_{t}+v_{i,\varepsilon}\partial_{\varepsilon}\bigr{)}f_{i}$
$\displaystyle=\omega^{\textrm{(FP)}}\Omega^{\textrm{(FP)}}_{i}+\omega^{\textrm{(B)}}\Omega^{\textrm{(B)}}_{i},$
(48)

where we have neglected the higher-order $\mathcal{O}\bigl{(}\delta
t^{3}\bigr{)}$ terms. Before treating the collision side of the equation, we
make the initial perturbation analysis on the advection terms and retain them
for further analysis. The analysis start with another common ansatz: only the
two lowest orders in $Kn$ are required to obtain the continuum NSEs where only
the coarse non-equilibrium effects in $f_{i}^{(1)}$ are considered in addition
to the equilibrium $f_{i}^{(0)}$ [51]. Substituting in the pertinent variables
from (46) yields the expanded advection side,

$\displaystyle\biggl{\\{}\delta
t\Bigl{[}\bigl{(}\epsilon\partial_{t}^{(1)}+\epsilon^{2}\partial_{t}^{(2)}\bigr{)}+\epsilon
v_{i,\alpha}\partial_{\alpha}^{(1)}\Bigr{]}+\frac{\delta
t^{2}}{2}\Bigl{[}\bigl{(}\epsilon\partial_{t}^{(1)}+\epsilon^{2}\partial_{t}^{(2)}\bigr{)}+\epsilon
v_{i,\alpha}\partial_{\alpha}^{(1)}\Bigr{]}\Bigl{[}\bigl{(}\epsilon\partial_{t}^{(1)}+\epsilon^{2}\partial_{t}^{(2)}\bigr{)}+\epsilon
v_{i,\beta}\partial_{\beta}^{(1)}\Bigr{]}\biggr{\\}}\bigl{(}f_{i}^{(0)}+\epsilon
f_{i}^{(1)}+\epsilon^{2}f_{i}^{(2)}\bigr{)},$ (49)

where contributions separated by orders in $\epsilon$ are,

$\displaystyle\epsilon^{0}:\qquad 0,$ (50a)
$\displaystyle\epsilon^{1}:\qquad\delta
t\bigl{(}\partial_{t}^{(1)}+v_{i,\alpha}\partial_{\alpha}^{(1)}\bigr{)}f_{i}^{(0)},$
(50b) $\displaystyle\epsilon^{2}:\qquad\Bigl{[}\delta
t\partial_{t}^{(2)}+\frac{\delta
t^{2}}{2}\bigl{(}\partial_{t}^{(1)}+v_{i,\alpha}\partial_{\alpha}^{(1)}\bigr{)}\bigl{(}\partial_{t}^{(1)}+v_{i,\beta}\partial_{\beta}^{(1)}\bigr{)}\Bigr{]}f_{i}^{(0)}$
$\displaystyle\qquad\qquad+\delta
t\Bigl{(}\partial_{t}^{(1)}+v_{i,\alpha}\partial_{\alpha}^{(1)}\Bigr{)}f_{i}^{(1)}.$
(50c)

These contributions will eventually equate with the corresponding-order
contributions from the analogously expanded collisions operators. To that end,
the BGK operator is expanded with the perturbation series yielding the
separated orders,

$\displaystyle\epsilon^{0}:\qquad\omega^{\textrm{(B)}}\bigl{(}f_{i}^{\textrm{(eq)}}-f_{i}^{(0)}\bigr{)},$
(51a) $\displaystyle\epsilon^{1}:\qquad-\omega^{\textrm{(B)}}f_{i}^{(1)},$
(51b) $\displaystyle\epsilon^{2}:\qquad-\omega^{\textrm{(B)}}f_{i}^{(2)}.$
(51c)

Proceeding with the FP operator, we similarly expand it by parts in
$\Omega^{\textrm{(FP)}}_{i}=\Omega^{\textrm{(FP,0)}}_{i}+\epsilon\Omega^{\textrm{(FP,1)}}_{i}+\epsilon^{2}\Omega^{\textrm{(FP,2)}}_{i}$
(wherein the perturbation order $n_{\epsilon}$ should not be confused with
those of the Hermite series), in which each of the orders are further
decomposed by $n_{\epsilon}=\\{0,1,2\\}$,

$\displaystyle\Omega^{(\textrm{FP},n_{\epsilon})}_{i}=w_{i}\left[\overline{\rho}^{(n_{\epsilon})}+\frac{v_{i,\alpha}}{v_{T}^{2}}\overline{\jmath}_{\alpha}^{(n_{\epsilon})}+\frac{v_{i,\alpha}v_{i,\beta}-v_{T}^{2}\delta_{\alpha\beta}}{2v_{T}^{4}}\overline{P}_{\alpha\beta}^{(n_{\epsilon})}\right].$
(52)

Thus, separation of orders simply yields,

$\displaystyle\epsilon^{0}:\qquad\omega^{\textrm{(FP)}}\Omega^{\textrm{(FP,0)}}_{i},$
(53a)
$\displaystyle\epsilon^{1}:\qquad\omega^{\textrm{(FP)}}\Omega^{\textrm{(FP,1)}}_{i},$
(53b)
$\displaystyle\epsilon^{2}:\qquad\omega^{\textrm{(FP)}}\Omega^{\textrm{(FP,2)}}_{i},$
(53c)

where the various functionals (38c)–(38) are linear combinations of the
perturbation series of the macroscopic variables (46e)–(46j), such that the
corresponding $\epsilon$-orders are also linearly retained.

Now, the next step is to recombine the different orders of the perturbation
series of the advection terms, and the BGK and FP operators, resulting in,

$\displaystyle\epsilon^{0}:\qquad
0=\omega^{\textrm{(B)}}\bigl{(}f_{i}^{\textrm{(eq)}}-f_{i}^{(0)}\bigr{)}$
$\displaystyle\qquad+\omega^{\textrm{(FP)}}w_{i}\left[\overline{\rho}^{(0)}+\frac{v_{i,\alpha}}{v_{T}^{2}}\overline{\jmath}_{\alpha}^{(0)}+\frac{v_{i,\alpha}v_{i,\beta}-v_{T}^{2}\delta_{\alpha\beta}}{2v_{T}^{4}}\overline{P}_{\alpha\beta}^{(0)}\right],$
(54a)
$\displaystyle\epsilon^{1}:\qquad\bigl{(}\partial_{t}^{(1)}+v_{i,\delta}\partial_{\delta}^{(1)}\bigr{)}f_{i}^{(0)}=-\frac{\omega^{\textrm{(B)}}}{\delta
t}f_{i}^{(1)}$ $\displaystyle\qquad+\frac{\omega^{\textrm{(FP)}}}{\delta
t}w_{i}\left[\overline{\rho}^{(1)}+\frac{v_{i,\alpha}}{v_{T}^{2}}\overline{\jmath}_{\alpha}^{(1)}+\frac{v_{i,\alpha}v_{i,\beta}-v_{T}^{2}\delta_{\alpha\beta}}{2v_{T}^{4}}\overline{P}_{\alpha\beta}^{(1)}\right],$
(54b) $\displaystyle\epsilon^{2}:\qquad\Bigl{[}\partial_{t}^{(2)}+\frac{\delta
t}{2}\bigl{(}\partial_{t}^{(1)}+v_{i,\delta}\partial_{\delta}^{(1)}\bigr{)}\bigl{(}\partial_{t}^{(1)}+v_{i,\varepsilon}\partial_{\varepsilon}^{(1)}\bigr{)}\Bigr{]}f_{i}^{(0)}$
$\displaystyle\qquad+\Bigl{(}\partial_{t}^{(1)}+v_{i,\delta}\partial_{\delta}^{(1)}\Bigr{)}f_{i}^{(1)}=-\frac{\omega^{\textrm{(B)}}}{\delta
t}f_{i}^{(2)}$ $\displaystyle+\frac{\omega^{\textrm{(FP)}}}{\delta
t}w_{i}\left[\overline{\rho}^{(2)}+\frac{v_{i,\alpha}}{v_{T}^{2}}\overline{\jmath}_{\alpha}^{(2)}+\frac{v_{i,\alpha}v_{i,\beta}-v_{T}^{2}\delta_{\alpha\beta}}{2v_{T}^{4}}\overline{P}_{\alpha\beta}^{(2)}\right],$
(54c)

to which end we make a third ansatz: there exists multiple roots in the
zeroth-order equality (54) where $f_{i}^{(0)}=f_{i}^{\textrm{(eq)}}$ is a
candidate as is conventionally the case in multi-scale analyses of the
Boltzmann equation. We can further infer that for the equality to hold, we
must require that,

$\displaystyle\overline{\rho}^{(0)}$ $\displaystyle=0,$ (55)
$\displaystyle\overline{\jmath}_{\alpha}^{(0)}$ $\displaystyle=0_{\alpha},$
(56) $\displaystyle\overline{P}_{\alpha\beta}^{(0)}$
$\displaystyle=0_{\alpha\beta}.$ (57)

The zero equalities are the only information we can extract from the zeroth-
order relation for now. Prior to taking the moments of the perturbation
relations, we can simplify the first and second-order relations. If we
factorize $f_{i}^{(1)}$ out of $\Omega^{\textrm{(FP,1)}}_{i}$ in (54) and
adopt the notation $\widehat{\Omega}_{i}^{\textrm{(FP,1)}}\circ f_{i}^{(1)}$
as the FP operational expression operating on the first-order populations, the
material derivative of the equilibrium populations can be written as,

$\displaystyle\epsilon^{1}:\qquad
D_{t}^{(\bm{v},1)}f_{i}^{(0)}=f_{i}^{(1)}\left(\frac{\omega^{\textrm{(FP)}}}{\delta
t}\widehat{\Omega}_{i}^{\textrm{(FP)}}-\frac{\omega^{\textrm{(B)}}}{\delta
t}\right),$ (58)

where we employ the notation
$D_{t}^{(\bm{v},1)}=\partial_{t}^{(1)}+v_{i,\delta}\partial_{\delta}^{(1)}$
for conciseness, and recognize that $\widehat{\Omega}_{i}^{\textrm{(FP)}}$ is
in fact independent of the perturbation order. This result is useful for
eliminating the first-order derivatives of $f_{i}^{(0)}$ in (54c), and
combining all prefactors of $f_{i}^{(1)}$, by substituting in this
reformulated $\epsilon^{1}$ contribution yielding,

$\displaystyle\epsilon^{2}:\qquad$
$\displaystyle\partial_{t}^{(2)}f_{i}^{(0)}+D_{t}^{(\bm{v},1)}\left(1+\frac{\omega^{\textrm{(FP)}}\widehat{\Omega}_{i}^{\textrm{(FP)}}}{2}-\frac{\omega^{\textrm{(B)}}}{2}\right)f_{i}^{(1)}$
$\displaystyle=-\frac{\omega^{\textrm{(B)}}}{\delta
t}f_{i}^{(2)}+\frac{\omega^{\textrm{(FP)}}}{\delta
t}\Omega^{\textrm{(FP,2)}}_{i}.$ (59)

This allows us to again in a reverse manner recast the relation with
$f_{i}^{(1)}$ isolated in (58), such that the second-order contributions now
obey,

$\displaystyle\epsilon^{2}:\qquad$
$\displaystyle\partial_{t}^{(2)}f_{i}^{(0)}+\Bigl{(}\partial_{t}^{(1)}+v_{i,\delta}\partial_{\delta}^{(1)}\Bigr{)}\Bigl{(}\partial_{t}^{(1)}+v_{i,\varepsilon}\partial_{\varepsilon}^{(1)}\Bigr{)}\delta
t$
$\displaystyle\times\left(1+\frac{2-\omega^{\textrm{(B)}}}{2\omega^{\textrm{(FP)}}\widehat{\Omega}_{i}^{\textrm{(FP)}}}-\frac{2+\omega^{\textrm{(FP)}}\widehat{\Omega}_{i}^{\textrm{(FP)}}}{2\omega^{\textrm{(B)}}}\right)f_{i}^{(0)}$
$\displaystyle\qquad=-\frac{\omega^{\textrm{(B)}}}{\delta
t}f_{i}^{(2)}+\frac{\omega^{\textrm{(FP)}}}{\delta
t}\Omega^{\textrm{(FP,2)}}_{i},$ (60)

where we have obtained an equilibrium differential equation with $f_{i}^{(0)}$
instead of $f_{i}^{(1)}$, where the non-equilibrium effects are administered
through the $\omega^{\textrm{(B)}},\omega^{\textrm{(FP)}}$ coefficients of
$f_{i}^{(0)}$. The perturbation equations (54), (58), (III.2.1) are now at a
stage where we can further the analysis by forming the moment equations.

####  III.2.2 Moment equations

Initially we can take the zero, second, and third moments of (54), from which
the continuity and momentum equation are obtained. The zeroth-moment equation
is obtained by summing both sides of the relation over $i$, the first-moment
equation by summing over $i$ multiplied by $v_{i,\varepsilon}$, the second-
moment equation by $v_{i,\varepsilon}v_{i,\zeta}$, and so forth. Thereto, we
can investigate the effects of perturbations on the different macroscopic
observables by including or excluding various moments. In the simplest
analysis, we can make the ansatz that the non-equilibrium moments [58],

$\displaystyle\sum_{i=0}^{Q-1}\bigl{\\{}f_{i}^{\textrm{(neq)}},g_{i}^{\textrm{(neq)}}\bigr{\\}}$
$\displaystyle=\bigl{\\{}\rho^{\textrm{(neq)}},\mathcal{E}^{\textrm{(neq)}}\bigr{\\}}=0,$
(61a)
$\displaystyle\sum_{i=0}^{Q-1}\bigl{\\{}f_{i}^{\textrm{(neq)}},g_{i}^{\textrm{(neq)}}\bigr{\\}}v_{i,\delta}$
$\displaystyle=\bigl{\\{}j_{\delta}^{\textrm{(neq)}},q_{\delta}^{\textrm{(neq)}}\bigr{\\}}=0_{\delta},$
(61b)
$\displaystyle\sum_{i=0}^{Q-1}\bigl{\\{}f_{i}^{\textrm{(neq)}},g_{i}^{\textrm{(neq)}}\bigr{\\}}v_{i,\delta}v_{i,\varepsilon}$
$\displaystyle=\bigl{\\{}P_{\delta\varepsilon}^{\textrm{(neq)}},R_{\delta\varepsilon}^{\textrm{(neq)}}\bigr{\\}}=0_{\delta\varepsilon}.$
(61c)

are all zero for the non-equilibrium states in $n_{\epsilon}\geq 1$. This
conventionally renders the continuity and Euler equations, and does not give
any insight into non-equilibrium effects. Now, as we are interested in
investigating the non-equilibrium behaviour of the FP equation on smaller
scales it may be worthwhile to contest the conventional zero-assumptions of
(61a–61c) and include the coarse first-order perturbations in the analysis. As
we know that the entropy-extremum principle corresponds to a minimization of
energy [14], we can include the contributions from
$\mathcal{E}^{(1)}=\sum_{i}g_{i}^{(1)}$. Additionally, for thermally-driven
phase-change processes such as pool boiling it is apt to consider the
fluctuations in the heat flux
$q_{\delta}^{(1)}=\sum_{i}g_{i}^{(1)}v_{i,\delta}$, and for cavitation which
is a stress-induced phenomenon, the momentum flux
$P_{\delta\varepsilon}^{(1)}=\sum_{i}f_{i}^{(1)}v_{i,\delta}v_{i,\varepsilon}$.
As the inclusion of $P_{\delta\varepsilon}^{(1)}$ in the second-moment
equation of the $\epsilon$-relation based on the Boltzmann equation results in
the viscous-stress tensor $\tau_{\delta\varepsilon}$ with the dynamic
viscosity $\mu=\bigl{(}1/\omega^{\textrm{(B)}}-1/2\bigr{)}p\delta t$ [2], one
could suspect that $\omega^{\textrm{(FP)}}$ (17) could be attributed to a
pseudo-viscosity from the $\gamma$-damping of thermal fluctuations—analogous
to the eddy-viscosity in turbulence models. Thus, we are interested in the
same non-equilibrium moments and how they are affected by the FP dynamics, and
as such assume that the values and correlations of all non-equilibrium
moments, except for that unknown and non-zero $P_{\delta\varepsilon}^{(1)}$
moment, are zero,

$\displaystyle\overline{\rho}^{(2)}$ $\displaystyle\doteq 0,$ (62a)
$\displaystyle
j_{\alpha}^{(1)}=j_{\alpha}^{(2)}=\overline{\jmath}_{\alpha}^{(1)}=\overline{\jmath}_{\alpha}^{(2)}$
$\displaystyle\doteq 0_{\alpha},$ (62b) $\displaystyle
P_{\alpha\beta}^{(2)}=\overline{P}_{\alpha\beta}^{(2)}$ $\displaystyle\doteq
0_{\alpha\beta},$ (62c)
$\displaystyle\mathcal{E}^{(1)}=\mathcal{E}^{(2)}=\overline{\mathcal{E}}^{(1)}=\overline{\mathcal{E}}^{(2)}$
$\displaystyle\doteq 0,$ (62d) $\displaystyle
q_{\alpha}^{(2)}=\overline{q}_{\alpha}^{(2)}$ $\displaystyle\doteq
0_{\alpha},$ (62e) $\displaystyle
R_{\alpha\beta}^{(1)}=R_{\alpha\beta}^{(2)}=\overline{R}_{\alpha\beta}^{(1)}=\overline{R}_{\alpha\beta}^{(2)}$
$\displaystyle\doteq 0_{\alpha\beta},$ (62f)

in addition to the non-zero counterpart,

$\displaystyle P_{\alpha\beta}^{(1)}\neq 0.$ (63a)

This ultimately allows us to take the zero-through-second moments of (54,
III.2.1), which usually are necessary to recover the NS equations. Moreover,
for the FP contributions, we are primarily interested in the non-equilibrium
effects to the momentum and energy equations, and as such we omit the density-
fluctuations in the zero-moment continuity equation. As for the Boltzmann
collision operator, we rely on the usual methodology of omitting non-
equilibrium contributions from the zero and first moment equations of the
perturbation equations [51]. From the moments of the first-order perturbation
expansion we get the (first-order) continuity, Euler, and a higher-order
equation,

$\displaystyle\partial_{t}^{(1)}\rho^{\textrm{(eq)}}+\partial_{\delta}^{(1)}j_{\delta}^{\textrm{(eq)}}$
$\displaystyle=0,$ (64a)
$\displaystyle\partial_{t}^{(1)}j_{\delta}^{\textrm{(eq)}}+\partial_{\varepsilon}^{(1)}P_{\delta\varepsilon}^{\textrm{(eq)}}$
$\displaystyle=0,$ (64b)
$\displaystyle\partial_{t}^{(1)}P_{\delta\varepsilon}^{\textrm{(eq)}}+\partial_{\zeta}^{(1)}Q_{\delta\varepsilon\zeta}^{\textrm{(eq)}}$
$\displaystyle=-\frac{\omega^{\textrm{(B)}}}{\delta
t}P_{\delta\varepsilon}^{(1)}+\frac{\omega^{\textrm{(FP)}}}{\delta
t}\widehat{\Omega}^{\textrm{(FP)}}P_{\delta\varepsilon}^{(1)},$ (64c) in which
we exploited the zero-non-equilibrium moments of $f_{i}^{(1)}$, as well as
factorized out $\sum_{i=0}^{Q-1}f_{i}$ from $\Omega^{\textrm{(FP)}}$ leaving
the moment of the operational expression itself
$\widehat{\Omega}^{\textrm{(FP)}}$. We will show later that the last moment
equation, as in the LBM, recovers the viscous-stress tensor, albeit with
additional stresses manifested by
$\widehat{\Omega}^{\textrm{(FP)}}P_{\delta\varepsilon}^{(1)}$. To that end,
the explicit formalism of $P_{\delta\varepsilon}^{(1)}$ is not recovered
directly from (64c) but rather appears in the first-moment of the second-order
perturbation expansion (III.2.1). Lastly, the equilibrium moments are defined
from the equilibrium populations (11) as, $\displaystyle\rho^{\textrm{(eq)}}$
$\displaystyle=\sum_{i=0}^{Q-1}f_{i}^{\textrm{(eq)}}=\rho,$ (64d)
$\displaystyle j_{\delta}^{\textrm{(eq)}}$
$\displaystyle=\sum_{i=0}^{Q-1}f_{i}^{\textrm{(eq)}}v_{i,\delta}=\rho
u_{\delta},$ (64e) $\displaystyle P_{\delta\varepsilon}^{\textrm{(eq)}}$
$\displaystyle=\sum_{i=0}^{Q-1}f_{i}^{\textrm{(eq)}}v_{i,\delta}v_{i,\varepsilon}=\rho
u_{\delta}u_{\varepsilon}+p\delta_{\delta\varepsilon},$ (64f) $\displaystyle
Q_{\delta\varepsilon\zeta}^{\textrm{(eq)}}$
$\displaystyle=\sum_{i=0}^{Q-1}f_{i}^{\textrm{(eq)}}v_{i,\delta}v_{i,\varepsilon}v_{i,\zeta}=\rho
u_{\delta}u_{\varepsilon}u_{\zeta}+p\bigl{[}u\delta\bigr{]}_{\delta\varepsilon\zeta},$
(64g)

where
$\bigl{[}u\delta\bigr{]}_{\delta\varepsilon\zeta}=u_{\delta}\delta_{\varepsilon\zeta}+u_{\varepsilon}\delta_{\delta\zeta}+u_{\zeta}\delta_{\delta\varepsilon}$.
Applying the product rule to $\partial_{\delta}^{(1)}(\rho u_{\delta})$ in the
continuity equation (64a) enables rewriting the material derivative as,

$\displaystyle D_{t}^{(1)}\rho=-\rho\partial_{\delta}^{(1)}u_{\delta}.$ (65)

This material derivative appears in the product-rule expanded momentum
equation considering the equilibrium moments $j_{\delta}^{\textrm{(eq)}}$ and
$P_{\delta\varepsilon}^{\textrm{(eq)}}$, which after some algebra allows us to
rewrite the momentum equation as,

$\displaystyle
D_{t}^{(1)}u_{\delta}=-\frac{1}{\rho}\partial_{\varepsilon}^{(1)}p\delta_{\delta\varepsilon}.$
(66)

Noting that the same perturbation series apply to $g_{i}$, its corresponding
moment equation hierarchy is,

$\displaystyle\partial_{t}^{(1)}\mathcal{E}^{\textrm{(eq)}}+\partial_{\delta}^{(1)}q_{\delta}^{\textrm{(eq)}}$
$\displaystyle=0,$ (67a)
$\displaystyle\partial_{t}^{(1)}q_{\delta}^{\textrm{(eq)}}+\partial_{\varepsilon}^{(1)}R_{\delta\varepsilon}^{\textrm{(eq)}}$
$\displaystyle=-\frac{\omega^{\textrm{(B)}}}{\delta
t}q_{\delta}^{(1)}+\frac{\omega^{\textrm{(FP)}}}{\delta
t}\widehat{\Omega}^{\textrm{(FP)}}q_{\delta}^{(1)},$ (67b) where the moments
of the equilibrium populations (12) read,
$\displaystyle\mathcal{E}^{\textrm{(eq)}}$
$\displaystyle=\sum_{i=0}^{Q-1}g_{i}^{\textrm{(eq)}}=2\rho E,$ (67c)
$\displaystyle q_{\delta}^{\textrm{(eq)}}$
$\displaystyle=\sum_{i=0}^{Q-1}g_{i}^{\textrm{(eq)}}v_{i,\delta}=2\rho
u_{\delta}H,$ (67d) $\displaystyle R_{\delta\varepsilon}^{\textrm{(eq)}}$
$\displaystyle=\sum_{i=0}^{Q-1}g_{i}^{\textrm{(eq)}}v_{i,\delta}v_{i,\varepsilon}=2\rho
u_{\delta}u_{\varepsilon}\bigl{(}H+p/\rho\bigr{)}+2pH\delta_{\delta\varepsilon},$
(67e) and $H=E+p/\rho=e+u^{2}/2+p/\rho$ is the total enthalpy. All of the
equilibrium moments except $Q_{\delta\varepsilon\zeta}^{\textrm{(eq)}}$ are
symmetric tensors. To arrive at the NSF equations, the second $g_{i}$-moment
equation is not required, and thus we have excluded it.

After some expansion with the product rule, and substituting the van der Waals
internal energy (40), the derivative relation $T\bigl{(}\partial p/\partial
T\bigr{)}_{v}=a\rho^{2}+p$ on a vdW-EOS basis, and the rewritten momentum
equation (65), it can be shown that the first-order contribution to the
temperature equation is,

$\displaystyle D_{t}^{(1)}T=-\frac{T}{\rho C_{v}}\biggl{(}\frac{\partial
p}{\partial T}\biggr{)}_{v}\partial_{\delta}^{(1)}u_{\delta}.$ (68)

Furthermore, an analogous pressure equation can be established by considering
that $p=p(\rho,T)$ where applying the corollary chain rule,

$\displaystyle D_{t}^{(1)}p=\biggl{(}\frac{\partial
p}{\partial\rho}\biggr{)}_{T}D_{t}^{(1)}\rho+\biggl{(}\frac{\partial
p}{\partial T}\biggr{)}_{\rho}D_{t}^{(1)}T,$ (69)

can be rewritten with the first-order continuity and energy equation
contributions (65,68), as well as the speed-of-sound (45), yielding,

$\displaystyle
D_{t}^{(1)}p=-\rho\varsigma^{2}\partial_{\delta}^{(1)}u_{\delta}.$ (70)

This relation will be exploited later to derive the viscous stresses in the
moment equations. As of now we have obtained the necessary zero and first-
moment equations of $f_{i},g_{i}$, and are now left with deriving an
operational expression for the moment $P_{\delta\varepsilon}^{(1)}$ terms in
the second-moment equation (64c). This last part of the analysis revolves
around the second-order perturbations (III.2.1) from which the zero moment
becomes,

$\displaystyle\partial_{t}^{(2)}\rho$ $\displaystyle=0,$ (71)

and consequently we can conclude that the density predominantly evolves on the
$\epsilon$ scales via (65). Formulating the first-moment equation of the
$\epsilon^{2}$ scales necessitates more algebraic work. First we note that in
the expanded first-moment equation stemming from (III.2.1) successive
instances of the momentum equation (66) appear and significantly simplify the
result to,

$\displaystyle\partial_{t}^{(2)}\rho
u_{\delta}=\partial_{\varepsilon}^{(1)}\left(\partial_{t}^{(1)}P_{\delta\varepsilon}^{\textrm{(eq)}}+\partial_{\zeta}^{(1)}Q_{\delta\varepsilon\zeta}^{\textrm{(eq)}}\right)\delta
t\biggl{[}\left(\frac{1}{\omega^{\textrm{(B)}}}-1\right)$
$\displaystyle+\left(\frac{\omega^{\textrm{(FP)}}}{2\omega^{\textrm{(B)}}}\widehat{\Omega}^{\textrm{(FP)}}-\frac{1-\omega^{\textrm{(B)}}/2}{\omega^{\textrm{(FP)}}}\frac{1}{\widehat{\Omega}^{\textrm{(FP)}}}\right)\biggr{]},$
(72)

where the moment $\widehat{\Omega}^{\textrm{(FP)}}$ of the operator is
$\sum_{i=0}^{Q-1}\widehat{\Omega}_{i}^{\textrm{(FP)}}$. Nevertheless, given
that
$\widehat{\Omega}^{\textrm{(FP)}}=\widehat{\Omega}^{\textrm{(FP)}}(\bm{x},\bm{u},t)$
we cannot immediately exclude the relaxation-frequency-dependent terms, except
for $\bigl{(}1/\omega^{\textrm{(B)}}-1\bigr{)}$, in the derivative terms, and
thus need to dissect the equation further. After some algebra it follows that
the second-order contributions to the momentum equation take the form,

$\displaystyle\partial_{t}^{(2)}u_{\delta}=\frac{\delta
t}{\rho}\biggl{\\{}\left(\frac{1}{\omega^{\textrm{(B)}}}-1\right)\partial_{\varepsilon}^{(1)}\left(\partial_{t}^{(1)}P_{\delta\varepsilon}^{\textrm{(eq)}}+\partial_{\zeta}^{(1)}Q_{\delta\varepsilon\zeta}^{\textrm{(eq)}}\right)$
$\displaystyle\quad+\frac{\omega^{\textrm{(FP)}}}{2\omega^{\textrm{(B)}}}\partial_{\varepsilon}^{(1)}\Bigl{[}\widehat{\Omega}^{\textrm{(FP)}}\left(\partial_{t}^{(1)}P_{\delta\varepsilon}^{\textrm{(eq)}}+\partial_{\zeta}^{(1)}Q_{\delta\varepsilon\zeta}^{\textrm{(eq)}}\right)+\left(P_{\delta\varepsilon}^{\textrm{(eq)}}\partial_{t}^{(1)}+Q_{\delta\varepsilon\zeta}^{\textrm{(eq)}}\partial_{\zeta}^{(1)}\right)\widehat{\Omega}^{\textrm{(FP)}}\Bigr{]}$
$\displaystyle\quad-\frac{1-\omega^{\textrm{(B)}}/2}{\omega^{\textrm{(FP)}}}\partial_{\varepsilon}^{(1)}\biggl{[}\left(\widehat{\Omega}^{\textrm{(FP)}}\right)^{-1}\left(\partial_{t}^{(1)}P_{\delta\varepsilon}^{\textrm{(eq)}}+\partial_{\zeta}^{(1)}Q_{\delta\varepsilon\zeta}^{\textrm{(eq)}}\right)+\left(P_{\delta\varepsilon}^{\textrm{(eq)}}\partial_{t}^{(1)}+Q_{\delta\varepsilon\zeta}^{\textrm{(eq)}}\partial_{\zeta}^{(1)}\right)\left(\widehat{\Omega}^{\textrm{(FP)}}\right)^{-1}\biggr{]}\biggr{\\}},$
(73)

in which the superposition of the three prefactors of
$\partial_{\varepsilon}^{(1)}(\ \cdots)$ recover the dynamic viscosity in the
thermohydrodynamic limit that together with
$\partial_{t}^{(1)}P_{\delta\varepsilon}^{\textrm{(eq)}}+\partial_{\zeta}^{(1)}Q_{\delta\varepsilon\zeta}^{\textrm{(eq)}}$
constitute the viscous-stress tensor, as we will show in the following. In the
result we can observe three major components to the stresses: the ordinary
Boltzmann-equation contributions by
$\bigl{(}1/\omega^{\textrm{(B)}}-1\bigr{)}$ as well as contributions from the
FP operator with both the material derivative of itself (either in the
continuous (6) or derived Hermite-space formalism) linearly dependent on the
moments, as well as more FPE-analogous material derivatives of the moments
operated on by the operator.

In an analogous analysis of the perturbation expansion of $g_{i}$, only the
zero-moment is required for recovering the Fourier equation. Considering the
perturbation expansion of not $f_{i}$ but $g_{i}$ in (III.2.1) the resulting
moment equation simplifies to,

$\displaystyle\partial_{t}^{(2)}T=\frac{1}{2\rho
C_{v}}\biggl{\\{}\partial_{\delta}^{(1)}\left(\partial_{t}^{(1)}q_{\delta}^{\textrm{(eq)}}+\partial_{\varepsilon}^{(1)}R_{\delta\varepsilon}^{\textrm{(eq)}}\right)\delta
t\biggl{[}\left(\frac{1}{\omega^{\textrm{(B)}}}-1\right)$
$\displaystyle\quad+\left(\frac{\omega^{\textrm{(FP)}}}{2\omega^{\textrm{(B)}}}\widehat{\Omega}^{\textrm{(FP)}}-\frac{1-\omega^{\textrm{(B)}}/2}{\omega^{\textrm{(FP)}}}\frac{1}{\widehat{\Omega}^{\textrm{(FP)}}}\right)\biggr{]}-2\rho
u_{\delta}\partial_{t}^{(2)}u_{\delta}\biggr{\\}},$ (74)

which requires exploiting the thermodynamic pressure (41), the internal energy
(40), the continuity equation (65), and the aforementioned pressure equation
(70). Finally,
$\widehat{\Omega}^{\textrm{(FP)}}=\widehat{\Omega}^{\textrm{(FP)}}(\bm{x},\bm{u},t)$,
necessitates further expansion yielding,

$\displaystyle\partial_{t}^{(2)}T=\frac{1}{2\rho C_{v}}\biggl{\\{}\delta
t\left(\frac{1}{\omega^{\textrm{(B)}}}-1\right)\partial_{\delta}^{(1)}\left(\partial_{t}^{(1)}q_{\delta}^{\textrm{(eq)}}+\partial_{\varepsilon}^{(1)}R_{\delta\varepsilon}^{\textrm{(eq)}}\right)-2\rho
u_{\delta}\partial_{t}^{(2)}u_{\delta}$ $\displaystyle\quad+\delta
t\left(\frac{\omega^{\textrm{(FP)}}}{2\omega^{\textrm{(B)}}}\right)\partial_{\delta}^{(1)}\left[\widehat{\Omega}^{\textrm{(FP)}}\left(\partial_{t}^{(1)}q_{\delta}^{\textrm{(eq)}}+\partial_{\varepsilon}^{(1)}R_{\delta\varepsilon}^{\textrm{(eq)}}\right)+\left(q_{\delta}^{\textrm{(eq)}}\partial_{t}^{(1)}+R_{\delta\varepsilon}^{\textrm{(eq)}}\partial_{\varepsilon}^{(1)}\right)\widehat{\Omega}^{\textrm{(FP)}}\right]$
$\displaystyle\quad-\delta
t\left(\frac{1-\omega^{\textrm{(B)}}/2}{\omega^{\textrm{(FP)}}}\right)\partial_{\delta}^{(1)}\left[\left(\widehat{\Omega}^{\textrm{(FP)}}\right)^{-1}\left(\partial_{t}^{(1)}q_{\delta}^{\textrm{(eq)}}+\partial_{\varepsilon}^{(1)}R_{\delta\varepsilon}^{\textrm{(eq)}}\right)+\left(q_{\delta}^{\textrm{(eq)}}\partial_{t}^{(1)}+R_{\delta\varepsilon}^{\textrm{(eq)}}\partial_{\varepsilon}^{(1)}\right)\left(\widehat{\Omega}^{\textrm{(FP)}}\right)^{-1}\right]\biggr{\\}}.$
(75)

This hierarchy of equations—specifically the continuity (65,71), the momentum
(66,73), and the energy (68,75) equations—form the foundation for recovering
the NSF equations by recombination of the perturbation components into the
native series formulation (46).

####  III.2.3 Thermo-hydrodynamic limit: compressible, fluctuating Navier-
Stokes-Fourier equations

At this point, we recover the fluctuating Navier-Stokes-Fourier equations from
the hierarchy of perturbed moment equations by recombining each of the orders
in $\epsilon$. By merging (65,71), and recognizing that
$\epsilon\partial_{t}^{(1)}+\epsilon^{2}\partial_{t}^{(2)}=\partial_{t}$ and
$\epsilon\partial_{\delta}^{(1)}=\partial_{\delta}$ we get the continuity
equation,

$\displaystyle\partial_{t}\rho+u_{\delta}\partial_{\delta}\rho=-\rho\partial_{\delta}u_{\delta}.$
(76)

Similarly, the momentum equation becomes,

$\displaystyle\rho D_{t}u_{\delta}$
$\displaystyle=-\partial_{\varepsilon}p\delta_{\delta\varepsilon}-\partial_{\varepsilon}\bigl{(}\tau_{\delta\varepsilon}+\widetilde{\tau}_{\delta\varepsilon}\bigr{)},$
(77a) where the viscous-stress tensor stemming exclusively from Boltzmann-type
collisions is,
$\displaystyle\tau_{\delta\varepsilon}=-\mu^{\textrm{(B)}}\left(\partial_{\delta}u_{\varepsilon}+\partial_{\varepsilon}u_{\delta}-\frac{2}{D}\bigl{(}\partial_{\zeta}u_{\zeta}\bigr{)}\delta_{\delta\varepsilon}\right)$
$\displaystyle-\zeta^{\textrm{(B)}}\bigl{(}\partial_{\zeta}u_{\zeta}\bigr{)}\delta_{\delta\varepsilon},$
(77b) with the shear and volume viscosities respectively as,
$\displaystyle\mu^{\textrm{(B)}}$
$\displaystyle=\left(\frac{1}{\omega^{\textrm{(B)}}}-1\right)p\delta t,$ (77c)
$\displaystyle\zeta^{\textrm{(B)}}$
$\displaystyle=\left(\frac{1}{\omega^{\textrm{(B)}}}-1\right)\left(\frac{D+2}{D}-\frac{\rho\varsigma^{2}}{p}\right)p\delta
t.$ (77d) In addition to this ``traditional'' stress tensor, we also get
``noisy'' stress tensor components perturbed directly by the FP operator and
its inverse, separated by complementing perturbative components (denoted by
``$+$'' and ``$-$'' subscripts) that concurrently amplify and attenuate
stresses, $\displaystyle\widetilde{\tau}_{\delta\varepsilon}$
$\displaystyle=\overbrace{\widehat{\Omega}^{\textrm{(FP)}}\left[-\mu^{\textrm{(FPB)}}_{+}\left(\partial_{\delta}u_{\varepsilon}+\partial_{\varepsilon}u_{\delta}-\frac{2}{D}\bigl{(}\partial_{\zeta}u_{\zeta}\bigr{)}\delta_{\delta\varepsilon}\right)-\zeta^{\textrm{(FPB)}}_{+}\bigl{(}\partial_{\zeta}u_{\zeta}\bigr{)}\delta_{\delta\varepsilon}\right]}^{(\textrm{I})}$
$\displaystyle+\frac{1}{\widehat{\Omega}^{\textrm{(FP)}}}\left[-\mu^{\textrm{(FPB)}}_{-}\left(\partial_{\delta}u_{\varepsilon}+\partial_{\varepsilon}u_{\delta}-\frac{2}{D}\bigl{(}\partial_{\zeta}u_{\zeta}\bigr{)}\delta_{\delta\varepsilon}\right)-\zeta^{\textrm{(FPB)}}_{-}\bigl{(}\partial_{\zeta}u_{\zeta}\bigr{)}\delta_{\delta\varepsilon}\right]\hskip
159.3356pt\Biggr{\\}}{\footnotesize(\textrm{II})}$
$\displaystyle+\mu^{\textrm{(FPB)}}_{+}\frac{P_{\delta\varepsilon}^{\textrm{(eq)}}}{p}\left\\{\left[\left(\partial_{v_{\zeta}^{(\eta)}}\widehat{\Omega}^{\textrm{(FP)}}\right)\frac{A_{\zeta}\rho\varsigma^{2}}{\gamma
m}+\left(\partial_{v_{T}^{2}}\widehat{\Omega}^{\textrm{(FP)}}\right)\frac{k_{B}T}{m\rho
c_{v}}\biggl{(}\frac{\partial p}{\partial
T}\biggr{)}_{v}\right]\bigl{(}\partial_{\zeta}u_{\zeta}\bigr{)}+\left(\partial_{u_{\zeta}}\widehat{\Omega}^{\textrm{(FP)}}\right)\frac{1}{\rho}\bigl{(}\partial_{\delta}p\bigr{)}\delta_{\delta\zeta}\right\\}\qquad\Biggr{\\}}{\footnotesize(\textrm{III})}$
$\displaystyle+\underbrace{\mu^{\textrm{(FPB)}}_{-}\frac{P_{\delta\varepsilon}^{\textrm{(eq)}}}{p}\left\\{\left[\left(\partial_{v_{\zeta}^{(\eta)}}\frac{1}{\widehat{\Omega}^{\textrm{(FP)}}}\right)\frac{A_{\zeta}\rho\varsigma^{2}}{\gamma
m}+\left(\partial_{v_{T}^{2}}\frac{1}{\widehat{\Omega}^{\textrm{(FP)}}}\right)\frac{k_{B}T}{m\rho
c_{v}}\biggl{(}\frac{\partial p}{\partial
T}\biggr{)}_{v}\right]\bigl{(}\partial_{\zeta}u_{\zeta}\bigr{)}+\left(\partial_{u_{\zeta}}\frac{1}{\widehat{\Omega}^{\textrm{(FP)}}}\right)\frac{1}{\rho}\bigl{(}\partial_{\delta}p\bigr{)}\delta_{\delta\zeta}\right\\}}_{(\textrm{IV})},$
(77e) in which the viscosities manifest due to both the Boltzmann and FP
relaxation frequencies,
$\displaystyle\mu^{\textrm{(FPB)}}_{+}=\left(\frac{\omega^{\textrm{(FP)}}}{2\omega^{\textrm{(B)}}}\right)p\delta
t,$
$\displaystyle\mu^{\textrm{(FPB)}}_{-}=\left(-\frac{1-\omega^{\textrm{(B)}}/2}{\omega^{\textrm{(FP)}}}\right)p\delta
t,$
$\displaystyle\zeta^{\textrm{(FPB)}}_{+}=\left(\frac{\omega^{\textrm{(FP)}}}{2\omega^{\textrm{(B)}}}\right)\left(\frac{D+2}{D}-\frac{\rho\varsigma^{2}}{p}\right)p\delta
t,$
$\displaystyle\zeta^{\textrm{(FPB)}}_{-}=\left(-\frac{1-\omega^{\textrm{(B)}}/2}{\omega^{\textrm{(FP)}}}\right)\left(\frac{D+2}{D}-\frac{\rho\varsigma^{2}}{p}\right)p\delta
t.$ (77f)

Thus, we get four pair-wise complementing ``fluctuating'' stress tensors in
(77): two traditional stress tensor components (I,II) pre-factored by FP-
operator $\widehat{\Omega}^{\textrm{(FP)}}$ and its inverse counterpart
$1/\widehat{\Omega}^{\textrm{(FP)}}$. The remaining two viscous-stress
components (III,IV) are due to the spatial and temporal derivatives of
$\widehat{\Omega}^{\textrm{(FP)}}$ in (77), and account for the diffusion
contributed from the FP operator itself as it varies due to external forces
(via $\bm{v}^{(\eta)}$ in the $\partial_{v_{\zeta}^{(\eta)}}$ terms), thermal
fluctuations (via $v_{T}^{2}$ in the $\partial_{v_{T}^{2}}$ terms), and
continuum velocity ($\bm{u}$ in the $\partial_{u_{\zeta}}$ terms),
respectively. We duly note that by fluctuating we imply that the stochastic
dynamics of the Langevin equation is captured by the dynamics in the velocity
space of the FPE rather than computed directly as discrete, stochastic events
sampled at random from a characteristic Gaussian white noise distribution, and
that this noise permeates the stresses.

In obtaining the velocity gradient and divergence terms in the stresses, i.e.
$\left(\partial_{\delta}u_{\varepsilon}+\partial_{\varepsilon}u_{\delta}-2/D\bigl{(}\partial_{\zeta}u_{\zeta}\bigr{)}\delta_{\delta\varepsilon}\right)$
and $\bigl{(}\partial_{\zeta}u_{\zeta}\bigr{)}\delta_{\delta\varepsilon}$, we
have exploited the expanded formulation,

$\displaystyle\partial_{t}^{(1)}P_{\delta\varepsilon}^{\textrm{(eq)}}$
$\displaystyle+\partial_{\zeta}^{(1)}Q_{\delta\varepsilon\zeta}^{\textrm{(eq)}}=p\left(\partial_{\delta}^{(1)}u_{\varepsilon}+\partial_{\varepsilon}^{(1)}u_{\delta}\right)$
$\displaystyle+\left(p-\rho\varsigma^{2}\right)\partial_{\zeta}^{(1)}u_{\zeta}\delta_{\delta\varepsilon},$
(78)

previously identified in (73). The noisy stresses arise analogously to the
stresses by the Boltzmann operator. During the recombination we obtain the
auxiliary components with the derivatives
$\left(P_{\delta\varepsilon}^{\textrm{(eq)}}\partial_{t}^{(1)}+Q_{\delta\varepsilon\zeta}^{\textrm{(eq)}}\partial_{\zeta}^{(1)}\right)\Bigl{\\{}\widehat{\Omega}^{\textrm{(FP)}};1/\widehat{\Omega}^{\textrm{(FP)}}\Bigr{\\}}$
imposed on the FP operational expressions and its inverse counterpart,
respectively. This form can be recast to a convenient material-derivative
basis as
$Q_{\delta\varepsilon\zeta}^{\textrm{(eq)}}=P_{\delta\varepsilon}^{\textrm{(eq)}}u_{\zeta}$
and thus,

$\displaystyle\left(P_{\delta\varepsilon}^{\textrm{(eq)}}\partial_{t}^{(1)}+Q_{\delta\varepsilon\zeta}^{\textrm{(eq)}}\partial_{\zeta}^{(1)}\right)=P_{\delta\varepsilon}^{\textrm{(eq)}}D_{t}^{(1)}.$
(79)

This trick enables us to reuse the chain rule to effectively rewrite
$P_{\delta\varepsilon}^{\textrm{(eq)}}D_{t}^{(1)}\Bigl{\\{}\widehat{\Omega}^{\textrm{(FP)}},1/\widehat{\Omega}^{\textrm{(FP)}}\Bigr{\\}}$
where we note that
$\widehat{\Omega}^{\textrm{(FP)}}=\widehat{\Omega}^{\textrm{(FP)}}\bigl{(}\bm{u}^{n},\bm{v}^{(\eta)},v_{T}^{n}(T)\bigr{)}$
is a function of velocity polynomials of various even orders $n$,

$\displaystyle D_{t}^{(1)}\widehat{\Omega}^{\textrm{(FP)}}$
$\displaystyle=\left(\frac{\partial\widehat{\Omega}^{\textrm{(FP)}}}{\partial
u_{\zeta}}\right)_{v_{\zeta}^{(\eta)},v_{T}^{2}}D_{t}^{(1)}u_{\zeta}$
$\displaystyle+\left(\frac{\partial\widehat{\Omega}^{\textrm{(FP)}}}{\partial
v_{\zeta}^{(\eta)}}\right)_{v_{T}^{2},u_{\zeta}}D_{t}^{(1)}v_{\zeta}^{(\eta)}$
$\displaystyle+\left(\frac{\partial\widehat{\Omega}^{\textrm{(FP)}}}{\partial
v_{T}^{2}}\right)_{u_{\zeta},v_{\zeta}^{(\eta)}}D_{t}^{(1)}v_{T}^{2}.$ (80)

Therein, we in addition to (66) can identify the following material
derivatives derived previously by reformulating
$v_{T}^{2},v_{\zeta}^{(\eta)}$. Firstly, we can exploit the previously derived
material derivative (68) that combined with the direction-independent, root-
mean-square formulation of the average thermal molecular speed
$v_{T}^{2}=k_{B}T/m$ for molecules of mass $m$ yields,

$\displaystyle D_{t}^{(1)}v_{T}^{2}$
$\displaystyle=\frac{k_{B}}{m}D_{t}^{(1)}T$
$\displaystyle=-\frac{k_{B}T}{m\rho c_{v}}\biggl{(}\frac{\partial p}{\partial
T}\biggr{)}_{v}\partial_{\zeta}^{(1)}u_{\zeta}.$ (81)

As we have not specified a particular force model in this study we rely on
Newton's second law of motion $\eta_{\zeta}=F_{\zeta}/m$ and (70) to formulate
the material derivative of $v_{\zeta}^{(\eta)}$,

$\displaystyle D_{t}^{(1)}v_{\zeta}^{(\eta)}$
$\displaystyle=-\frac{A_{\zeta}\rho\varsigma^{2}}{\gamma
m}\partial_{\zeta}^{(1)}u_{\zeta},$ (82)

where we have further introduced the unit area $A_{\zeta}$ as a consequence of
converting the force into a pressure $p=F_{\zeta}/A_{\zeta}$, such that
$A_{\zeta}$ in a simulation scenario would be the grid-cell cross-sectional
area. The non-ideal speed of sound $\varsigma^{2}$ was previously found in
(45). Consequently, we obtain the explicit result for,

$\displaystyle
P_{\delta\varepsilon}^{\textrm{(eq)}}D_{t}^{(1)}\widehat{\Omega}^{\textrm{(FP)}}=-P_{\delta\varepsilon}^{\textrm{(eq)}}\Biggl{\\{}\biggl{[}\left(\partial_{v_{\zeta}^{(\eta)}}\widehat{\Omega}^{\textrm{(FP)}}\right)\frac{A_{\zeta}\rho\varsigma^{2}}{\gamma
m}$
$\displaystyle\quad+\left(\partial_{v_{T}^{2}}\widehat{\Omega}^{\textrm{(FP)}}\right)\frac{k_{B}T}{m\rho
c_{v}}\biggl{(}\frac{\partial p}{\partial
T}\biggr{)}_{v}\biggr{]}\bigl{(}\partial_{\zeta}u_{\zeta}\bigr{)}$
$\displaystyle\quad+\left(\partial_{u_{\zeta}}\widehat{\Omega}^{\textrm{(FP)}}\right)\frac{1}{\rho}\bigl{(}\partial_{\delta}p\bigr{)}\delta_{\delta\zeta}\Biggr{\\}},$
(83)

and similarly for the $1/\widehat{\Omega}^{\textrm{(FP)}}$ destructive
components with $\mu^{\textrm{(FPB)}}_{-}$, where we in both cases recombine
$\partial_{\zeta}^{(1)}$ according to our previously defined perturbation
series.

To obtain the Fourier equation we consider (75) wherein we substitute with
(73). The shear velocities arise from the derivatives of the moments,

$\displaystyle\partial_{t}^{(1)}q_{\delta}^{\textrm{(eq)}}+\partial_{\varepsilon}^{(1)}R_{\delta\varepsilon}^{\textrm{(eq)}}=2\left(p-\rho\varsigma^{2}\right)\partial_{\zeta}^{(1)}u_{\zeta}u_{\delta}$
$\displaystyle+2pu_{\varepsilon}\left(\partial_{\varepsilon}^{(1)}u_{\delta}+\partial_{\delta}^{(1)}u_{\varepsilon}\right)+2p\partial_{\delta}^{(1)}h,$
(84)

where $h=e+p/\rho$ is the enthalpy. As noted by Reyhanian  _et al._ [2, 50]
the resulting heat flux would be $\bm{q}=-\mu\nabla h$ and consequently the
$g_{i}$-kinetic equation needs to be augmented with an energy correction term
$M_{0}=2\partial_{\alpha}\bigl{(}-\mu\partial_{\alpha}h+k\partial_{\alpha}T\bigr{)}$
to recover the correct Fourier law. For real gasses
$\partial_{\alpha}h=c_{p}\partial_{\alpha}T+v\bigl{(}1-\beta\bigr{)}\partial_{\alpha}p$
with $\beta=\rho\bigl{(}\partial v/\partial T\bigr{)}_{p}$ being the thermal-
expansion coefficient. After applying the product rule to the shear-velocity
terms, the substituted result can be recombined with the first-order
perturbation contribution (68) confirming the conservation of the Fourier
equation,

$\displaystyle\rho
c_{v}D_{t}T=-\bigl{(}\tau_{\delta\varepsilon}+\widetilde{\tau}_{\delta\varepsilon}\bigr{)}\partial_{\delta}u_{\varepsilon}-T\biggl{(}\frac{\partial
p}{\partial
T}\biggr{)}_{v}\partial_{\delta}u_{\delta}-\partial_{\delta}q_{\delta},$ (85)

where the normal and fluctuating shear stresses organize as expected in the
same formulation of the tensors (77,77).

In summary, we recover the compressible, fluctuating Navier-Stokes-Fourier
hierarchy of equations,

$\displaystyle D_{t}\rho$ $\displaystyle=-\rho\nabla\cdot\bm{u},$ (86)
$\displaystyle\rho D_{t}\bm{u}$ $\displaystyle=-\nabla
p-\nabla\cdot(\bm{\tau}+\widetilde{\bm{\tau}}),$ (87) $\displaystyle\rho
c_{v}D_{t}T$
$\displaystyle=-(\bm{\tau}+\widetilde{\bm{\tau}}):\nabla\bm{u}-T\biggl{(}\frac{\partial
p}{\partial T}\biggr{)}_{v}\nabla\cdot\bm{u}-\nabla\cdot\bm{q},$ (88)

in which the Fourier heat flux is defined as $\bm{q}=-k\nabla T$. To enable
simulations we superimpose the presented viscosities into effective shear and
volume viscosities,

$\displaystyle\mu$
$\displaystyle=\mu^{\textrm{(B)}}+\mu^{\textrm{(FPB)}}_{+}+\mu^{\textrm{(FPB)}}_{-}$
$\displaystyle=\left[\left(\frac{1}{\omega^{\textrm{(B)}}}-1\right)+\left(\frac{\omega^{\textrm{(FP)}}}{2\omega^{\textrm{(B)}}}\right)-\left(\frac{1-\omega^{\textrm{(B)}}/2}{\omega^{\textrm{(FP)}}}\right)\right]p\delta
t,$ (89) $\displaystyle\zeta$
$\displaystyle=\zeta^{\textrm{(B)}}+\zeta^{\textrm{(FPB)}}_{+}+\zeta^{\textrm{(FPB)}}_{-}$
$\displaystyle=\left[\left(\frac{1}{\omega^{\textrm{(B)}}}-1\right)+\left(\frac{\omega^{\textrm{(FP)}}}{2\omega^{\textrm{(B)}}}\right)-\left(\frac{1-\omega^{\textrm{(B)}}/2}{\omega^{\textrm{(FP)}}}\right)\right]$
$\displaystyle\qquad\times\left(\frac{D+2}{D}-\frac{\rho\varsigma^{2}}{p}\right)p\delta
t,$ (90)

where known initial values of $\mu,\zeta$ can be defined and subsequently used
to tune the time-step size $\delta t$, and the kinetic relaxation frequencies
$\omega^{\textrm{(B)}},\omega^{\textrm{(FP)}}(\gamma)$, which in turn sets the
thermal noise intensity governed by the Langevin equation. This result along
with the kinetic equations (13) constitute the foundation of our numerical
method.

##  IV Discussion

As the macroscopic observables $\bm{\lambda}(\bm{x},t)=\\{\rho,\bm{u},E\\}$
inherently are driven by the Langevin forces imposed via our LFPBE, it is
important to convince ourselves that the reciprocal behaviour between the
perturbations instigated by these Langevin forces and the response in
$\bm{\lambda}$ evolves in a thermodynamically consistent manner, especially
when the perturbations bring the thermodynamic state far away from
equilibrium. A systematic approach for assessing this is provided by the
celebrated Onsager relations [59], which are used to set the thermal noise
strength while obeying the fluctuation-dissipation theorem (FDT) and ensuring
that the simulated populations in the LFPBEs meet the detailed balance
conditions [46]. Not only can they be applied to the Markovian Ornstein-
Uhlenbeck process that constitutes the FPE, they can be used to observe the
FDT applied to macroscopic observables [60]. This is particularly important in
our derived framework as we have heuristically added BGK relaxation into our
kinetic equations (13) resulting in the viscous coefficients (89, III.2.3)
which do not yet provide a physics-informed procedure for setting the
fluctuation intensity and dissipation, independently of the BGK collisions.
Setting the noise intensity should additionally be done in tandem with mixture
rules for the BGK collisions which account for hard-sphere collisions in a
liquid/vapour mixture, relative to the FP interactions that are more important
in liquids. From the Langevin equation, we stated that FDT for the noise
intensity for the Weiner process $W(t)$ in the molecular velocities evolves
only in time according to the ensemble average in (2). However, as the
hydrodynamic dissipation attributed with the fluctuating stresses
$\widetilde{\bm{\tau}}$ (and possibly a corresponding heat flux
$\widetilde{\bm{q}}$ yet to be derived) are indirectly correlated with the
fluctuating fields
$\widetilde{\bm{\lambda}}(\bm{x},t)=\bigl{\\{}\widetilde{\rho},\widetilde{\bm{u}},\widetilde{E}\bigr{\\}}$
corresponding to the Weiner process, the fluctuating variables would need to
obey tailored fluctuation-dissipation balances. In addition, $\mu,\zeta$
should similarly obey a fluctuation-dissipation balance as they are functions
of $\widetilde{p},\widetilde{\rho}$ and also have equilibrium values
$\bigl{\\{}\mu^{\textrm{(eq)}},\zeta^{\textrm{(eq)}}\bigr{\\}}\bigl{(}p^{\textrm{(eq)}},\rho^{\textrm{(eq)}}\bigr{)}$.
Examples of such balances could be inspired by the work of Gallo [15].

The assumed white noise in our Weiner process underlying our Langevin and FPE
equations represents an initial prototype for a simulation framework of
thermal noise. The independence of the white noise spectrum is an idealization
and may not apply accurately to all fluids, and especially not impure fluids.
By exemplification, previous works on Brownian motion in water showed evidence
for yellow noise scaling with the square-root of the frequency in its weak-
noise spectrum in proximity to solid walls [61], effectively rendering the
thermal noise non-Markovian, i.e. the acceleration of particles is inherently
dependent on its past motion introducing a memory effect. Moreover, the
colored-noise amplitude depended strongly on the distance to the wall. As
such, extensions of our framework to impure, multi-component cavitating and
boiling fluids may need to revisit the frequency spectrum of the thermodynamic
noise.

We furthermore append a commentary on the implications of thermal fluctuations
on nucleation phenomena. Although thermal fluctuations are irrelevant in
disordered, homogeneous fluid domains, thermal noise is very important in
proximity to critical points, such as a spinodal where fluctuations can
perturb the fluid into a new stable thermodynamic state [55]. Prospectively,
assessing the fluctuating behaviour around these critical points should be
done rigorously on the basis of variables that exhibit sensitivity to the
noise strength and correlate with the local metastable state. To that end,
order-parameters $\psi(\bm{x},t)$, such as that evolved by the Cahn-Hillard-
Cook (CHC) equation
$\partial_{t}\psi=\nabla\cdot\left[\mathcal{D}\nabla\bigl{(}\delta\mathcal{F}(\psi)/\delta\psi\bigr{)}+\tilde{\theta}(\bm{x},t)\right]$
[62], which can give an indication of the mutual interaction and similarity
between non-homogenizing regions of a phase-transitioning fluid. The CHC
equation is informed by a diffusion coefficient $\mathcal{D}$, the Helmholtz
potential $\mathcal{F}=E(\psi)-TS(\psi)$, and lastly the noise term
$\tilde{\theta}$ that obeys the FDT [14, 41] responsible for the correct
equilibration towards the Maxwell-Boltzmann distribution. Usually, this order
term contributes an additional surface-tension force
$-\psi\nabla\mu_{\mathcal{F}}$ to the RHS of (87), where
$\mu_{\mathcal{F}}(\bm{x},t)=\delta\mathcal{F}(\psi)/\delta\psi$ is the
chemical potential. Additionally, the order parameter contributes a noisy-
momentum flux $\nabla\cdot\widetilde{\bm{\tau}}^{(\psi)}$ that also obeys the
FDT. Evolving an equation such as the CHC may provide critical insight into
metastable effects interacting with local thermodynamics and surface tension
forces at interfaces. We have not considered surface-tension forces in our
current single-component implementation of the FPBE, but it naturally becomes
relevant in multi-component extensions where multiple species separate at
anomalous rates. To that end, there already exist work on $N$-component
formulations of the FPE [63].

The entropy generation associated with irreversibilities from thermal
fluctuations can be exactly achieved by further mathematical analysis of the
FPE [33]. Ultimately, this contribution to the system entropy may be non-
trivial in phase-transitioning fluids in which perturbations from noise and
metastability, and therefore the characteristic relaxation rates for the phase
transitions, are fundamentally dictated by the entropy-extremum principle
[14]. This entropy generation in turn predicts the Helmholtz and Gibbs free-
energies $G(S)$ such that there is a coupling between the noise,
thermodynamics, and disorder in phase transitions. The free energy further
enables predicting the degree of metastability under the saturation curve
[15].

##  V Conclusion

To enable the simulation of the diffusive effects from thermal noise of
Langevin particles in thermal flows, we present a lattice-Fokker-Planck-
Boltzmann model of the phase-space continuous Fokker-Planck-Boltzmann equation
based on Hermite quadrature and the conventional $DmQn$ lattice models. As
reported in the literature, the lattice-Fokker-Planck equation usually suffers
from numerical instability at lower Mach numbers than the Bhatnagar-Gross-
Krook operator. Thus, to mitigate sources of instability in the Hermite series
in highly compressible and thermal flow scenarios, we have employed the
kinetic theory with discrete velocities
$\bm{v}_{i}=\sqrt{\theta}\bm{c}_{i}+\bm{u}$ from Particles-on-Demand [47],
which for the BGK operator has reported to exhibit unconditional instability.
This is facilitated by the adaptivity of these velocities to the local
thermohydrodynamics captured in the reduced pressure $\theta$ and the
continuum velocity $\bm{u}$. Moreover, this adaptivity dovetails nicely with
the introduction of the thermal noise and associated non-equilibrium
thermodynamic contributions of the Fokker-Planck as the fluctuations are
directly incorporated into the discrete velocity space.

By Chapman-Enskog multiscale analysis we show that the two resulting lattice
kinetic equations (13) with lattice-FP operators (38, 38) by proxy solve the
Navier-Stokes-Fourier equations, and crucially that the dissipative effects of
thermal noise induce viscous stresses $\widetilde{\bm{\tau}}$ where the
shearing components are damped by the dynamic viscosity $\mu$ (89) and the
diverging terms by the volume viscosity $\zeta$ (III.2.3). The viscosities are
tuned by independently setting the relaxation frequencies
$\omega^{\textrm{(FP)}}$ and $\omega^{\textrm{(B)}}$ of Fokker-Planck
interactions and BGK collisions, respectively, where $\omega^{\textrm{(FP)}}$
directly and uniquely defines the thermal-noise spectrum via the friction
factor $\gamma$ in the Langevin equation (1). Consequently, the lattice
equations in tandem with the relaxation frequencies enable numerical
simulations of metastable fluids with phase transitions induced by thermal
white noise.

###### Acknowledgements.

We wish to thank C.S. Wang (University of British Columbia, Okanagan) for his
thorough introduction to the lattice-Boltzmann method, and A. Garcia and M.
Sharifi (University of British Columbia, Okanagan) for fruitful discussions.
K.J.P. sincerely thanks T. Colonius (California Institute of Technology) for
advisory support during his stay at the Institute, that kick-started this
research, J.R. Chreim (California Institute of Technology) for numerous,
inspiring discussions on the current research. Moreover, this research would
not have been possible without the financial support of the Natural Sciences
and Engineering Research Council of Canada (NSERC)—under Development Grant No.
519885, the American-Scandinavian Foundation, Hedorf's Foundation, and Mitacs
Globalink (IT27768). K.J.P. conceptualized this research, reviewed the
literature, conducted the mathematical derivations, and wrote the manuscript
draft. J.R.B. provided review and editing of the written material and advisory
support for the project.

##  Appendix A Eigenfunctions of the Fokker-Planck and acceleration operators

In this section, we aim to derive the explicit eigenfunctions of the
continuous FP and acceleration operators. A crucial step in discretizing the
operators is to recast the entire continuous operators in convenient,
functional forms that do not include any gradient terms. Reiterating, we seek
two Hermite series that approximates each of the continuous operators,

$\displaystyle\Bigl{\\{}\Omega^{(\eta)}_{i},\Omega^{\textrm{(FP)}}_{i}\Bigr{\\}}=w_{i}\sum_{l=0}^{K}\frac{1}{v_{T}^{2l}l!}\Bigl{\\{}\Omega^{(\eta)}_{\underline{\alpha}},\Omega^{\textrm{(FP)}}_{\underline{\alpha}}\Bigr{\\}}\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v}_{i}),$
(91a) where the Hermite coefficients integrate the operators acting on the
population $f(\bm{v};\bm{x},t)$, in phase space,
$\displaystyle\Bigl{\\{}\Omega^{(\eta)}_{\underline{\alpha}},\Omega^{\textrm{(FP)}}_{\underline{\alpha}}\Bigr{\\}}(l)=\int
d\bm{v}~{}\Bigl{\\{}\widehat{\Omega}^{(\eta)},\widehat{\Omega}^{\textrm{(FP)}}\Bigr{\\}}[f]\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v}),$
(91b)

In the following derivations we treat the hat-operators as actions on $[f]$,
but omitting the population itself. $f$ in phase-space also needs to be
discretized with its Hermite coefficients
$\mathcal{F}_{\underline{\alpha}}^{(l)}$,

$\displaystyle
f(\bm{v};\bm{x},t)=w(\bm{v})\sum_{l=0}^{K}\frac{1}{v_{T}^{2l}l!}\mathcal{F}_{\underline{\alpha}}^{(l)}(\bm{x},t)\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v}).$
(92)

As it turns out, one can exploit the eigenvalue property of the operators on
the products
$\widehat{\Omega}^{\textrm{(FP)}}\bigl{[}w(\bm{v})\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\bigr{]}$
identified in the Hermite series of both operators, to eliminate the
derivative terms and arrive at an explicit formulation. This eigenvalue
property is derived from the Gaussian,

$\displaystyle w(\bm{v})$ $\displaystyle=\frac{1}{(2\pi
v_{T}^{2})^{\nicefrac{{D}}{{2}}}}\exp\left(-\frac{\bigl{(}\bm{v}-\bm{u}\bigr{)}^{2}}{2v_{T}^{2}}\right),$
(93)

that is the common, null-space of the acceleration, FP, and BGK operators
[32]. The velocity-derivative of this results in the eigenvector
$-\sqrt{\theta}c_{\alpha}/v_{T}^{2}$ as,

$\displaystyle\partial_{v_{\alpha}}w(\bm{v})=-\frac{\sqrt{\theta}c_{\alpha}}{v_{T}^{2}}w(\bm{v}),$
(94)

where $\theta=p/(\rho RT_{L})$ is the reduced pressure [2], $c_{\alpha}$ the
conventional peculiar velocity for any $DmQn$ lattice model, and
$v_{T}^{2}=\sqrt{RT}^{2}$ the squared mean thermal speed. Considering
$\Omega^{\textrm{(FP)}}$ excluding $f$, and omitting the implied notation
$w(\bm{v})$ we can expand the product,

$\displaystyle\widehat{\Omega}^{\textrm{(FP)}}\Bigl{[}w(\bm{v})\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\Bigr{]}=\gamma\partial_{v_{\beta}}\bigl{(}v_{\beta}+v_{T}^{2}\partial_{v_{\beta}}\bigr{)}\Bigl{[}w\mathcal{H}_{\underline{\alpha}}^{(l)}\Bigr{]}$
$\displaystyle\quad=\gamma\Bigl{[}\bigl{(}w\mathcal{H}_{\underline{\alpha}}^{(l)}+v_{\beta}\mathcal{H}_{\underline{\alpha}}^{(l)}\partial_{v_{\beta}}w+v_{\beta}w\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}\bigr{)}$
$\displaystyle\quad+\bigl{(}v_{T}^{2}\partial_{v_{\beta}}\partial_{v_{\beta}}w\mathcal{H}_{\underline{\alpha}}^{(l)}\bigr{)}\Bigr{]}.$
(95)

Using the eigenproperty (94) in the first parenthesis and rewriting the
second, we obtain,

$\displaystyle\widehat{\Omega}^{\textrm{(FP)}}\Bigl{[}w(\bm{v})\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\Bigr{]}=\gamma\Biggl{[}\Biggl{(}w\mathcal{H}_{\underline{\alpha}}^{(l)}-v_{\beta}\mathcal{H}_{\underline{\alpha}}^{(l)}\frac{\sqrt{\theta}c_{\beta}}{v_{T}^{2}}w$
$\displaystyle\quad+v_{\beta}w\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}\Biggr{)}+v_{T}^{2}\partial_{v_{\beta}}\left(w\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}+\mathcal{H}_{\underline{\alpha}}^{(l)}\partial_{v_{\beta}}w\right)\Biggr{]}.$
(96)

For now we cannot manipulate the first parenthesis' terms further and we turn
our attention to the second where we can use the eigenproperty as well as the
chain-rule,

$\displaystyle(\textrm{A6})=\gamma\Biggl{[}\left(w\mathcal{H}_{\underline{\alpha}}^{(l)}-v_{\beta}\mathcal{H}_{\underline{\alpha}}^{(l)}\frac{\sqrt{\theta}c_{\beta}}{v_{T}^{2}}w+v_{\beta}w\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}\right)$
$\displaystyle\quad+v_{T}^{2}\left(\partial_{v_{\beta}}\bigl{\\{}w\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}\bigr{\\}}-\frac{1}{v_{T}^{2}}\partial_{v_{\beta}}\bigl{\\{}\mathcal{H}_{\underline{\alpha}}^{(l)}\sqrt{\theta}c_{\beta}w\bigr{\\}}\right)\Biggr{]}$
$\displaystyle\quad=\gamma\Biggl{[}\left(w\mathcal{H}_{\underline{\alpha}}^{(l)}-v_{\beta}\mathcal{H}_{\underline{\alpha}}^{(l)}\frac{\sqrt{\theta}c_{\beta}}{v_{T}^{2}}w+v_{\beta}w\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}\right)$
$\displaystyle\quad+v_{T}^{2}\Biggl{(}\bigl{\\{}w\partial_{v_{\beta}}\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}$
$\displaystyle\quad+\underbrace{\bigl{(}\partial_{v_{\beta}}w\bigr{)}}_{-\frac{\sqrt{\theta}c_{\beta}}{v_{T}^{2}}w}\bigl{(}\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}\bigr{)}\bigr{\\}}-\frac{1}{v_{T}^{2}}\partial_{v_{\beta}}\bigl{\\{}\mathcal{H}_{\underline{\alpha}}^{(l)}\sqrt{\theta}c_{\beta}w\bigr{\\}}\Biggr{)}\Biggr{]},$
(97)

which yields the expanded result,

$\displaystyle(\textrm{A7})=\gamma\Biggl{[}\left(w\mathcal{H}_{\underline{\alpha}}^{(l)}-v_{\beta}\mathcal{H}_{\underline{\alpha}}^{(l)}\frac{\sqrt{\theta}c_{\beta}}{v_{T}^{2}}w+v_{\beta}w\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}\right)$
$\displaystyle+\bigl{(}v_{T}^{2}w\partial_{v_{\beta}}\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}-\sqrt{\theta}c_{\beta}w\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}\bigr{)}-\bigl{(}\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}\sqrt{\theta}c_{\beta}w\bigr{)}\Biggr{]}.$
(98)

Now we only need to expand the term in the third parenthesis that cf. the
product rule yields,

$\displaystyle(\textrm{A8})=\gamma\Biggl{[}w\mathcal{H}_{\underline{\alpha}}^{(l)}-v_{\beta}\mathcal{H}_{\underline{\alpha}}^{(l)}\frac{\sqrt{\theta}c_{\beta}}{v_{T}^{2}}w+v_{\beta}w\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}$
$\displaystyle\quad+v_{T}^{2}w\partial_{v_{\beta}}\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}-\sqrt{\theta}c_{\beta}w\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}$
$\displaystyle\quad-\bigl{(}\sqrt{\theta}c_{\beta}w\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}+\mathcal{H}_{\underline{\alpha}}^{(l)}\sqrt{\theta}c_{\beta}\cancelto{~{}\eqref{eq:eigenvector}}{\partial_{v_{\beta}}w}\bigr{)}\Biggr{]},$
(99)

wherein $\sqrt{\theta},c_{\beta}$ are independent of $v_{\beta}$ and thus drop
out with the partial derivative, such that we are left with,

$\displaystyle(\textrm{A9})=\gamma\Biggl{[}w\mathcal{H}_{\underline{\alpha}}^{(l)}-v_{\beta}\mathcal{H}_{\underline{\alpha}}^{(l)}\frac{\sqrt{\theta}c_{\beta}}{v_{T}^{2}}w+v_{\beta}w\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}$
$\displaystyle\quad+v_{T}^{2}w\partial_{v_{\beta}}\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}-\sqrt{\theta}c_{\beta}w\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}$
$\displaystyle\quad-\left(\sqrt{\theta}c_{\beta}w\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}-\mathcal{H}_{\underline{\alpha}}^{(l)}\frac{\theta
c_{\beta}c_{\beta}}{v_{T}^{2}}w\right)\Biggr{]}.$ (100)

In there we substitute in $v_{\beta}=\sqrt{\theta}c_{\beta}+u_{\beta}$, such
that a number of terms cancel out and we find the fully expanded result,

$\displaystyle\widehat{\Omega}^{\textrm{(FP)}}\Bigl{[}w(\bm{v})\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\Bigr{]}=\gamma
w\Biggl{[}1+u_{\beta}\left(\partial_{v_{\beta}}-\frac{\sqrt{\theta}c_{\beta}}{v_{T}^{2}}\right)$
$\displaystyle\quad+v_{T}^{2}\partial_{v_{\beta}}\partial_{v_{\beta}}-\sqrt{\theta}c_{\beta}\partial_{v_{\beta}}\Biggr{]}\mathcal{H}_{\underline{\alpha}}^{(l)}.$
(101)

To transform the partial-derivative terms, we can exploit the Hermite relation
[1],

$\displaystyle
v_{T}^{2}\partial_{v_{\beta}}\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}=\bigl{(}v_{\beta}\partial_{v_{\beta}}-l\bigr{)}\mathcal{H}_{\underline{\alpha}}^{(l)},$
(102)

so that our previous result can be rewritten to,

$\displaystyle(\textrm{A12})$ $\displaystyle=\gamma
w\left[1-u_{\beta}\frac{\sqrt{\theta}c_{\beta}}{v_{T}^{2}}-l+2u_{\beta}\partial_{v_{\beta}}\right]\mathcal{H}_{\underline{\alpha}}^{(l)},$
(103)

where we exploited that $v_{\beta}-\sqrt{\theta}c_{\beta}=u_{\beta}$. Finally
the remaining partial derivative can be transformed using the recurrence
relation [1, 54],

$\displaystyle\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}=\frac{1}{v_{T}^{2}}\Bigl{[}v_{\beta}\mathcal{H}_{\underline{\alpha}}^{(l)}-\mathcal{H}_{\beta\underline{\alpha}}^{(l+1)}\Bigr{]},$
(104)

that enables us to write,

$\displaystyle(\textrm{A14})$ $\displaystyle=\gamma
w\left[1-u_{\beta}\frac{\sqrt{\theta}c_{\beta}}{v_{T}^{2}}-l\right]\mathcal{H}_{\underline{\alpha}}^{(l)}$
$\displaystyle+\gamma
w\frac{2u_{\beta}}{v_{T}^{2}}\Bigl{[}v_{\beta}\mathcal{H}_{\underline{\alpha}}^{(l)}-\mathcal{H}_{\beta\underline{\alpha}}^{(l+1)}\Bigr{]}.$
(105)

Thereto, we can again exploit that
$v_{\beta}=\sqrt{\theta}c_{\beta}+u_{\beta}$ and factorize the
$v_{\beta}\mathcal{H}_{\underline{\alpha}}^{(l)}$ term into the first
parenthesis yielding,

$\displaystyle(\textrm{A16})$ $\displaystyle=\gamma
w\left[1+u_{\beta}\frac{v_{\beta}+u_{\beta}}{v_{T}^{2}}-l\right]\mathcal{H}_{\underline{\alpha}}^{(l)}-\gamma
w\frac{2u_{\beta}}{v_{T}^{2}}\Bigl{[}\mathcal{H}_{\beta\underline{\alpha}}^{(l+1)}\Bigr{]},$
(106)

where we further exploited that
$2v_{\beta}-\sqrt{\theta}c_{\beta}=v_{\beta}+u_{\beta}$. Thus, we find the
explicit relation that can be used in the $(l\leq K)$–Hermite series expansion
for $\Omega^{\textrm{(FP)}}_{i}$ in (91a):

$\displaystyle\widehat{\Omega}^{\textrm{(FP)}}\Bigl{[}w(\bm{v})\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\Bigr{]}=$
$\displaystyle\uline{\gamma
w\left[\left(1+u_{\beta}\frac{v_{\beta}+u_{\beta}}{v_{T}^{2}}-l\right)\mathcal{H}_{\underline{\alpha}}^{(l)}-\frac{2u_{\beta}}{v_{T}^{2}}\mathcal{H}_{\beta\underline{\alpha}}^{(l+1)}\right]}.$
(107)

A similar and conveniently easier analysis can be carried out for the
acceleration-operator product
$\widehat{\Omega}^{(\eta)}\Bigl{[}w(\bm{v})\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\Bigr{]}$.
Firstly, an expansion of this exploiting the linearity of the composite
partial derivative renders,

$\displaystyle\widehat{\Omega}^{(\eta)}\Bigl{[}w(\bm{v})\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\Bigr{]}$
$\displaystyle=-\eta_{\beta}\partial_{v_{\beta}}\Bigl{[}w\mathcal{H}_{\underline{\alpha}}^{(l)}\Bigr{]}$
$\displaystyle=-\eta_{\beta}\Bigl{[}w\partial_{v_{\beta}}\mathcal{H}_{\underline{\alpha}}^{(l)}+\mathcal{H}_{\underline{\alpha}}^{(l)}\partial_{v_{\beta}}w\Bigr{]}.$
(108)

Using the eigenproperty (94) and the Hermite recurrence relation (104) this
can be rewritten to,

$\displaystyle(\textrm{A16})$
$\displaystyle=-\eta_{\beta}w\left[\frac{\Bigl{(}v_{\beta}-\sqrt{\theta}c_{\beta}\Bigr{)}\mathcal{H}_{\underline{\alpha}}^{(l)}}{v_{T}^{2}}-\frac{\mathcal{H}_{\beta\underline{\alpha}}^{(l+1)}}{v_{T}^{2}}\right].$
(109)

Thus, the explicit eigenfunction we were seeking for the acceleration operator
reads,

$\displaystyle\widehat{\Omega}^{(\eta)}\Bigl{[}w(\bm{v})\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\Bigr{]}=\uline{\frac{\eta_{\beta}w}{v_{T}^{2}}\Bigl{[}\mathcal{H}_{\beta\underline{\alpha}}^{(l+1)}-u_{\beta}\mathcal{H}_{\underline{\alpha}}^{(l)}\Bigr{]}},$
(110)

that we analogously can use in the expansion of $\Omega^{(\eta)}_{i}$ in
(91a). In App. B we use these to compute the Hermite coefficients for the
series expansions.

##  Appendix B Hermite coefficients for the operators

Now that we have computed the operator eigenfunctions
$\bigl{\\{}\widehat{\Omega}^{\textrm{(FP)}};\widehat{\Omega}^{(\eta)}\bigr{\\}}\Bigl{[}w(\bm{v})\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\Bigr{]}$
we can establish the Hermite coefficients
$\Omega^{\textrm{(FP)}}_{\underline{\alpha}},\Omega^{(\eta)}_{\underline{\alpha}}$
by evaluating the integrals (91b). For the FP operator coefficient as a
composite function of $f(\bm{v};\bm{x},t)$, we need to utilize the Hermite
series (22) for that population, such that the integral can be rewritten to
account for $\mathcal{F}_{\underline{\alpha}}^{(l)}$,

$\displaystyle\Omega^{\textrm{(FP)}}_{\underline{\gamma}}$ $\displaystyle=\int
d\bm{v}~{}\mathcal{H}_{\underline{\gamma}}^{(m)}\Bigl{[}\widehat{\Omega}^{\textrm{(FP)}}\circ
f\Bigr{]}$
$\displaystyle=\sum_{l=0}^{K}\frac{\mathcal{F}_{\underline{\alpha}}^{(l)}}{v_{T}^{2l}l!}\int
d\bm{v}~{}\mathcal{H}_{\underline{\gamma}}^{(m)}\times\widehat{\Omega}^{\textrm{(FP)}}\Bigl{[}w(\bm{v})\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\Bigr{]},$
(111)

where we also relocated $w\mathcal{H}_{\underline{\alpha}}^{(l)}$ such that
the previous result in App. A conveniently can be identified and exploited
therein. Substituting in the eigenfunction we can expand the integral to,

$\displaystyle\int
d\bm{v}~{}\mathcal{H}_{\underline{\gamma}}^{(m)}\Bigl{[}\widehat{\Omega}^{\textrm{(FP)}}\circ
f\Bigr{]}=\sum_{l=0}^{K}\frac{\mathcal{F}_{\underline{\alpha}}^{(l)}}{v_{T}^{2l}l!}\int
d\bm{v}~{}\mathcal{H}_{\underline{\gamma}}^{(m)}$ $\displaystyle\times\gamma
w\left[\left(1+u_{\beta}\frac{v_{\beta}+u_{\beta}}{v_{T}^{2}}-l\right)\mathcal{H}_{\underline{\alpha}}^{(l)}-\frac{2u_{\beta}}{v_{T}^{2}}\mathcal{H}_{\beta\underline{\alpha}}^{(l+1)}\right].$
(112)

In order to compute the integral, we can exploit the Hermite orthonormality
relation [1, 54] which reads,

$\displaystyle\int
d\bm{v}~{}w(\bm{v})\frac{1}{v_{T}^{l+m}}\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\mathcal{H}_{\underline{\beta}}^{(m)}(\bm{v})=\delta_{lm}\delta_{\underline{\alpha}\underline{\beta}}^{(l)},$
(113)

where $\delta_{\underline{\alpha}\underline{\beta}}^{(l)}$ is the sum of
products of $l$ Kronecker $\bm{\delta}$'s where each of them associates with
one index in $\underline{\alpha}$ and one in $\underline{\beta}$ such that
there are $l!$ terms in $\delta_{\underline{\alpha}\underline{\beta}}^{(l)}$
[54]. $\delta_{lm}$ simply denotes that only tensors of identical rank $l=m$
contribute to the sum, and all other combinations $l\neq m$ evaluates to zero
for the integral. As we will have to apply (113) twice in (B2), once one each
of the
$\mathcal{H}_{\underline{\alpha}}^{(l)},\mathcal{H}_{\beta\underline{\alpha}}^{(l+1)}$
terms, we treat them separately as follows, and combine their results after.
The $l$-rank term becomes,

$\displaystyle\sum_{m=l=0}^{K}\frac{\mathcal{F}_{\underline{\alpha}}^{(l)}}{l!}\times\gamma\left[1+u_{\beta}\frac{v_{\beta}+u_{\beta}}{v_{T}^{2}}-l\right]\int
d\bm{v}~{}w\frac{1}{v_{T}^{2l}}\mathcal{H}_{\underline{\gamma}}^{(m)}\mathcal{H}_{\underline{\alpha}}^{(l)}$
$\displaystyle\quad=\frac{\mathcal{F}_{\underline{\alpha}}^{(l)}\gamma}{l!}\left[1+u_{\beta}\frac{v_{\beta}+u_{\beta}}{v_{T}^{2}}-l\right]\delta_{\underline{\gamma}\underline{\alpha}}^{(l)},$
(114)

where we can observe that $1/v_{T}^{l+m}=1/v_{T}^{2l}$, as $m=l$, and that
only $\mathcal{F}_{\underline{\alpha}}^{(l)}$ associates with the indices in
$\delta_{\underline{\gamma}\underline{\alpha}}^{(l)}$. Thus, we get the
discrete sum,

$\displaystyle\textrm{(B4)}=\gamma\left[1+u_{\beta}\frac{v_{\beta}+u_{\beta}}{v_{T}^{2}}-l\right]\mathcal{F}^{(l)}_{\gamma_{1},\dots,\gamma_{l}},$
(115)

where the $l!$ sums of $\mathcal{F}_{\underline{\alpha}}^{(l)}$ cancelled out
with the $1/l!$ prefix in (B4). Now we repeat the analysis for the second part
of (B2) with
$-2u_{\beta}\mathcal{H}_{\beta\underline{\alpha}}^{(l+1)}/v_{T}^{2}$. In this
case, we note that only the $l=m-1$ ranks of
$\mathcal{H}_{\beta\underline{\alpha}}^{(l+1)}$ yield a non-zero contribution
to the integral. Consequently, we seek to compute the integral,

$\displaystyle\sum_{m=0}^{K}\frac{\mathcal{F}_{\underline{\alpha}}^{(l)}}{v_{T}^{2l}l!}\int
d\bm{v}~{}\mathcal{H}_{\underline{\gamma}}^{(m)}\times\gamma
w\left[-\frac{2u_{\beta}}{v_{T}^{2}}\mathcal{H}_{\beta\underline{\alpha}}^{(l+1)}\right]$
$\displaystyle=-\sum_{m=0}^{K}\frac{2\gamma
u_{\beta}\mathcal{F}_{\underline{\alpha}}^{(l)}}{v_{T}^{2l+2}l!}\int
d\bm{v}~{}w\mathcal{H}_{\underline{\gamma}}^{(m)}\mathcal{H}_{\beta\underline{\alpha}}^{(l+1)}$
(116) $\displaystyle=-\sum_{\begin{subarray}{c}l=m-1\\\
=0\end{subarray}}^{K}\frac{2\gamma
u_{\beta}\mathcal{F}_{\underline{\alpha}}^{(m-1)}}{(m-1)!}\int
d\bm{v}~{}w\frac{1}{v_{T}^{2(m-1)+2}}\mathcal{H}_{\underline{\gamma}}^{(m)}\mathcal{H}_{\beta\underline{\alpha}}^{(m)}.$
(117)

Therein, we note that
$\mathcal{H}_{\underline{\gamma}}^{(m)},\mathcal{H}_{\beta\underline{\alpha}}^{(m)}$
are of the same ranks, and that indeed $v_{T}^{-[2(m-1)+2]}=1/v_{T}^{l+m}$,
such that we again can apply the orthonormality relation (113). Upon
substitution we get the result,

$\displaystyle\textrm{(B11)}=-\frac{2\gamma}{(m-1)!}u_{\beta}\mathcal{F}_{\underline{\alpha}}^{(m-1)}\delta_{\underline{\gamma},\beta\underline{\alpha}}^{(m-1)},$
(118)

where both $u_{\beta},\mathcal{F}_{\underline{\alpha}}^{(m-1)}$ associate with
the indices in $\delta_{\underline{\gamma},\beta\underline{\alpha}}^{(m-1)}$.
Ultimately, we compute the discrete sum as,

$\displaystyle\textrm{(B12)}=-2\gamma\Bigl{[}u_{\gamma_{1}}\mathcal{F}_{\gamma_{2},\dots,\gamma_{m}}^{(m-1)}+\cdots+u_{\gamma_{m}}\mathcal{F}_{\gamma_{1},\dots,\gamma_{m-1}}^{(m-1)}\Bigr{]}.$
(119)

Combining the two results we get the coefficients,

$\displaystyle\Omega^{\textrm{(FP)}}_{\underline{\gamma}}(l)=\uline{\gamma\left[1+u_{\beta}\frac{v_{\beta}+u_{\beta}}{v_{T}^{2}}-l\right]\mathcal{F}_{\underline{\gamma}}^{(l)}}$
$\displaystyle\quad\uline{-2\gamma\Bigl{[}u_{\gamma_{1}}\mathcal{F}_{\gamma_{2},\dots,\gamma_{m}}^{(m-1)}+\cdots+u_{\gamma_{m}}\mathcal{F}_{\gamma_{1},\dots,\gamma_{m-1}}^{(m-1)}\Bigr{]}}.$
(120)

We repeat the procedure for the acceleration operator now, to which end we
seek to compute the coefficients,

$\displaystyle\Omega^{(\eta)}_{\underline{\gamma}}$ $\displaystyle=\int
d\bm{v}~{}\mathcal{H}_{\underline{\gamma}}^{(m)}\Bigl{[}\widehat{\Omega}^{(\eta)}\circ
f\Bigr{]}$
$\displaystyle=\sum_{l=0}^{K}\frac{\mathcal{F}_{\underline{\alpha}}^{(l)}}{v_{T}^{2l}l!}\int
d\bm{v}~{}\mathcal{H}_{\underline{\gamma}}^{(m)}\times\widehat{\Omega}^{(\eta)}\Bigl{[}w(\bm{v})\mathcal{H}_{\underline{\alpha}}^{(l)}(\bm{v})\Bigr{]}$
(121)
$\displaystyle=\sum_{l=0}^{K}\frac{\mathcal{F}_{\underline{\alpha}}^{(l)}}{v_{T}^{2l}l!}\int
d\bm{v}~{}\mathcal{H}_{\underline{\gamma}}^{(m)}\times\frac{\eta_{\beta}w}{v_{T}^{2}}\Bigl{[}\mathcal{H}_{\beta\underline{\alpha}}^{(l+1)}-u_{\beta}\mathcal{H}_{\underline{\alpha}}^{(l)}\Bigr{]}$
(122)
$\displaystyle=\sum_{l=0}^{K}\frac{\mathcal{F}_{\underline{\alpha}}^{(l)}\eta_{\beta}}{l!}\int
d\bm{v}~{}w\frac{1}{v_{T}^{2l+2}}\mathcal{H}_{\underline{\gamma}}^{(m)}\Bigl{[}\mathcal{H}_{\beta\underline{\alpha}}^{(l+1)}-u_{\beta}\mathcal{H}_{\underline{\alpha}}^{(l)}\Bigr{]}$
(123)
$\displaystyle=\sum_{l=0}^{K}\frac{\mathcal{F}_{\underline{\alpha}}^{(l)}\eta_{\beta}}{l!}\Biggl{[}\int
d\bm{v}~{}w\frac{1}{v_{T}^{2l+2}}\mathcal{H}_{\underline{\gamma}}^{(m)}\mathcal{H}_{\beta\underline{\alpha}}^{(l+1)}$
$\displaystyle\quad-\frac{u_{\beta}}{v_{T}^{2}}\int
d\bm{v}~{}w\frac{1}{v_{T}^{2l}}\mathcal{H}_{\underline{\gamma}}^{(m)}\mathcal{H}_{\underline{\alpha}}^{(l)}\Biggr{]},$
(124)

where we in the last line group the terms in a convenient for aligned with the
orthonormality relations. Again, letting $l\leftarrow m-1$ and then applying
(113), the first integral becomes:

$\displaystyle\sum_{l=0}^{K}\frac{\mathcal{F}_{\underline{\alpha}}^{(l)}\eta_{\beta}}{l!}\left[\int
d\bm{v}~{}w\frac{1}{v_{T}^{2l+2}}\mathcal{H}_{\underline{\gamma}}^{(m)}\mathcal{H}_{\beta\underline{\alpha}}^{(l+1)}\right]$
$\displaystyle\quad=\sum_{\begin{subarray}{c}l=m-1\\\
=0\end{subarray}}^{K}\frac{\mathcal{F}_{\underline{\alpha}}^{(m-1)}\eta_{\beta}}{(m-1)!}\left[\int
d\bm{v}~{}w\frac{1}{v_{T}^{2(m-1)+2}}\mathcal{H}_{\underline{\gamma}}^{(m)}\mathcal{H}_{\beta\underline{\alpha}}^{(m)}\right]$
$\displaystyle\quad=\frac{1}{(m-1)!}\mathcal{F}_{\underline{\alpha}}^{(m-1)}\eta_{\beta}\delta_{\underline{\gamma},\beta\underline{\alpha}}^{(m-1)}$
$\displaystyle\quad=\Bigl{[}\eta_{\gamma_{1}}\mathcal{F}_{\gamma_{2},\dots,\gamma_{m}}^{(m-1)}+\cdots+\eta_{\gamma_{m}}\mathcal{F}_{\gamma_{1},\dots,\gamma_{m-1}}^{(m-1)}\Bigr{]}.$
(125)

The second integral in (B18) can be directly manipulated with (113) such that
we obtain the following:

$\displaystyle\sum_{l=0}^{K}\frac{\mathcal{F}_{\underline{\alpha}}^{(l)}\eta_{\beta}}{l!}\left[-\frac{u_{\beta}}{v_{T}^{2}}\int
d\bm{v}~{}w\frac{1}{v_{T}^{2l}}\mathcal{H}_{\underline{\gamma}}^{(m)}\mathcal{H}_{\underline{\alpha}}^{(l)}\right]$
$\displaystyle\quad=-\sum_{m=l=0}^{K}\frac{\eta_{\beta}u_{\beta}\mathcal{F}_{\underline{\alpha}}^{(l)}}{v_{T}^{2}l!}\left[\int
d\bm{v}~{}w\frac{1}{v_{T}^{2l}}\mathcal{H}_{\underline{\gamma}}^{(m)}\mathcal{H}_{\underline{\alpha}}^{(l)}\right]$
$\displaystyle\quad=-\frac{\eta_{\beta}u_{\beta}}{v_{T}^{2}l!}\mathcal{F}_{\underline{\alpha}}^{(l)}\delta_{\underline{\gamma},\underline{\alpha}}^{(l)}$
$\displaystyle\quad=-\frac{\eta_{\beta}u_{\beta}}{v_{T}^{2}}\Bigl{[}\mathcal{F}^{(l)}_{\gamma_{1},\dots,\gamma_{l}}\Bigr{]}.$
(126)

Combining the two integrals yields the sought after coefficients,

$\displaystyle\Omega^{(\eta)}_{\underline{\gamma}}(l)=$
$\displaystyle\uline{\Bigl{[}\eta_{\gamma_{1}}\mathcal{F}_{\gamma_{2},\dots,\gamma_{m}}^{(m-1)}+\cdots+\eta_{\gamma_{m}}\mathcal{F}_{\gamma_{1},\dots,\gamma_{m-1}}^{(m-1)}\Bigr{]}-\frac{\eta_{\beta}u_{\beta}}{v_{T}^{2}}\Bigl{[}\mathcal{F}^{(l)}_{\gamma_{1},\dots,\gamma_{l}}\Bigr{]}},$
(127)

that we together with (B) can use to compute the full Hermite series.
Furthermore, instead of using the $f$ distribution we can also use the $g$
distributions given that we use $\mathcal{G}_{\underline{\alpha}}^{(l)}$
instead of $\mathcal{F}_{\underline{\alpha}}^{(l)}$ in the coefficients.

## References

  * Moroni  _et al._ [2006a] D. Moroni, B. Rotenberg, J.-P. Hansen, S. Succi, and S. Melchionna, Solving the fokker-planck kinetic equation on a lattice, Physical Review E 73, 066707 (2006a).