MnLargeSymbols’164 MnLargeSymbols’171

# Dissipation in Hydrodynamics from Micro- to Macroscale:
Wisdom from Boltzmann and Stochastic Thermodynamics

Danilo Forastiere<EMAIL_ADDRESS>Complex Systems and Statistical
Mechanics, Department of Physics and Materials Science, University of
Luxembourg, L-1511 Luxembourg City, Luxembourg  Department of Physics and
Astronomy “Galileo Galilei”, University of Padova, Via F. Marzolo, 8, I-35131
Padova, Italy  INFN, Sezione di Padova, via Marzolo 8, I-35131 Padova, Italy
Francesco Avanzini<EMAIL_ADDRESS>Complex Systems and
Statistical Mechanics, Department of Physics and Materials Science, University
of Luxembourg, L-1511 Luxembourg City, Luxembourg  Department of Chemical
Sciences, University of Padova, Via F. Marzolo, 1, I-35131 Padova, Italy
Massimiliano Esposito<EMAIL_ADDRESS>Complex Systems and
Statistical Mechanics, Department of Physics and Materials Science, University
of Luxembourg, L-1511 Luxembourg City, Luxembourg

###### Abstract

We show that macroscopic irreversible thermodynamics for viscous fluids can be
derived from exact information-theoretic thermodynamic identities valid at the
microscale. Entropy production in particular is a measure of the loss of many-
particle correlations in the same way in which it measures the loss of system-
reservoirs correlations in stochastic thermodynamics (ST). More specifically,
we first show that boundary conditions at the macroscopic level define a
natural decomposition of the entropy production rate (EPR) in terms of
thermodynamic forces multiplying their conjugated currents, as well as a
change in suitable nonequilibrium potential that acts as a Lyapunov function
in absence of forces. At the microscale, for isolated Hamiltonian systems, we
identify the exact identities at the origin of these dissipative
contributions. Indeed, the molecular chaos hypothesis, which gives rise to the
Boltzmann equation at the mesoscale, leads to a positive rate of loss of many-
particle correlations which we identify with the Boltzmann entropy production
rate. By generalizing the Boltzmann equation to account for boundaries with
nonuniform temperature and nonzero velocity, and resorting to the
Chapman–Enskog expansion, we reproduce the macroscopic theory we started from.
We also show that an equivalent mesoscale description in terms of a linearized
Boltzmann equation obeys local detailed balance (LDB) and thus reproduces a ST
theory. Our work unambiguously demonstrates the information-theoretical origin
of thermodynamic notions of entropy and dissipation in macroscale irreversible
thermodynamics.

##  I Introduction

Hydrodynamics is a theory of crucial importance across science. It provides
tools to deal with systems spanning many orders of magnitude in terms of
spatial scale, and thus it is fundamental for our understanding of phenomena
ranging from chemical physics and biology to planetary science and
astrophysics de Groot and Mazur (1962); De Zarate and Sengers (2006); Phillips
_et al._ (2012); Peixoto and Oort (1992); Chandrasekhar (2013). One of its
simplest formulations describes an isotropic fluid in terms of five dynamical
fields, total energy (kinetic, potential and internal), the three components
of momentum, and mass. These are conserved at the microscopic level and thus
give rise to balance equations at the macroscopic level. These equations are
however not closed from a dynamical standpoint. In irreversible
thermodynamics, a complete thermodynamic description of a fluid assumes the
existence of a local entropy which is solely a function of internal energy,
volume and particle number, as in equilibrium. This so-called local
equilibrium assumption defines temperature, pressure and chemical potential
fields as partial derivatives of the entropy field and links them to the total
energy, momentum and mass density fields. A second assumption, which defines
linear irreversible thermodynamics, is that the fluxes featured in the balance
equations (heat flux and pressure tensor) can be linearized in terms of the
gradients of temperature and velocity. This linear regime assumption is
essential to make the balance equations a closed dynamical system de Groot and
Mazur (1962); Gallavotti (2002).

The microscopic underpinnings of irreversible thermodynamics have been the
object of intense scrutiny for more than a century and a half Chapman and
Cowling (1990); de Groot and Mazur (1962); Lanford (1975); Cercignani (1988);
Spohn (2012); Gaspard (1998, 2022). The main problems, which remain still open
today, are related to understanding the emergence of a macroscopic
irreversible behavior, from a microscopic reversible many-body dynamics. The
first problem concerns how to go from a microscopic dynamics of the gas
(formulated in terms of the many-body probability distribution of all the
fluid particles) which is hamiltonian and thus reversible, to a mesoscopic,
irreversible closed evolution equation for the single-particle probability
distribution (for example, via the Boltzmann equation Chapman and Cowling
(1990); Cercignani (1988); Gaspard (2022)). The second problem is the
derivation of an irreversible thermodynamics at the macroscopic level from the
mesoscopic level (for example, using the Chapman–Enskog expansion Chapman and
Cowling (1990); Prigogine (1949); Van Kampen (1987)). Irreversibility (_i.e._
, a nonnegative EPR) in this context is understood to arise from the molecular
chaos hypothesis which neglects correlations amongst particles building after
molecular collisions Spohn (2012); Gaspard (1998).

More recently, significant progress has been achieved in so-called stochastic
thermodynamics Jarzynski (1997); Seifert (2012); Pigolotti and Peliti (pear),
which offers a different approach to deal with the mesoscopic level of
description. Here the typical setup consists of a system in contact with one
or multiple reservoirs. At the microscopic level of description the full
probability distribution contains both the system and reservoirs degrees of
freedom. Exact identities corresponding to the first and second law of
thermodynamics have been derived when considering a class of initial
conditions where the system is prepared in an arbitrary state but uncorrelated
from the reservoirs which are at equilibrium. The entropy production is non-
negative and takes the forms of a relative entropy between the exact total
probability distribution at a given time and the exact system probability
distribution times the equilibrium reservoir distribution. It thus measures
the information lost when both the system-reservoirs correlations and the
displacement from equilibrium of the reservoirs are neglected Esposito  _et
al._ (2010); Ptaszyński and Esposito (2019). Various assumptions need to be
made in order to derive a closed dynamics for the system probability
distribution only (obtained by tracing out the reservoirs) Breuer  _et al._
(2002). The most common approach relies on the Born–Markov approximation
(assuming fast reservoirs and weak system-reservoir coupling) which plays a
role similar to the molecular chaos hypothesis in the context of fluids
McLennan Jr (1959). It also allows to express all the thermodynamics
quantities entering the first and second law in terms of the system
probability distribution only and to prove that the EPR is non-negative Soret
_et al._ (2022). The system dynamics only perceives the linear response
properties of the reservoirs via transition rates — or drift and diffusion
coefficients in a continuous description — which as a result satisfy the so-
called local detailed balance property Esposito  _et al._ (2009); Esposito
(2012); Falasco and Esposito (2021); Maes (2021). This latter condition is
known to be essential to build a consistent nonequilibrium thermodynamic
description for a system described by a stochastic dynamics Kurchan (1998);
Crooks (1999); Lebowitz and Spohn (1999); Seifert (2005); Rao and Esposito
(2018). These approaches hold for a large class of mesoscopic models such as
Markov jump processes, overdamped and underdamped Fokker–Planck equations
Seifert (2005); Esposito and Van den Broeck (2010a); Van den Broeck and
Esposito (2010); Esposito and Van den Broeck (2010b); Cavina and Esposito
(2022).

Furthermore, the formulation of ST for Markov jump processes enables a
decomposition of the EPR which provides a constructive identification of the
nonequilibrium potential and of the global thermodynamic forces driving the
system out of equilibrium Rao and Esposito (2018). Close to equilibrium, it
reproduces the classical results from linear irreversible thermodynamics, such
as Onsager reciprocal relations and minimum entropy production principle
Lebowitz and Spohn (1999); Andrieux and Gaspard (2004); Forastiere  _et al._
(2022). This decomposition of the EPR remains valid for macroscopic system
when assuming local equilibrium but _without_ resorting to the linear regime
assumption. This has been done for chemical reaction networks with diffusion
(but neglecting inertial or viscous effects) where deterministic field
theories ensue Falasco  _et al._ (2018); Avanzini  _et al._ (2019). This
construction reveals substantial differences between the global and the local
thermodynamics, because all the global conservation laws can be broken at the
local level Avanzini  _et al._ (2019). These macroscopic thermodynamic
theories with local equilibrium can be explicitly constructed from an
underlying stochastic thermodynamics by scaling up particle numbers (see
respectively Refs. Lazarescu  _et al._ (2019); Herpich  _et al._ (2020);
Freitas  _et al._ (2021) for chemical reaction networks, Potts models and
electronic circuits, and Refs. Freitas and Esposito (2022); Falasco and
Esposito (2022) for the general arguments). These approaches correspond to
overdamped dynamics where inertial effects are neglected, but generalization
to underdamped dynamics seems possible and could be used to derive
hydrodynamics.

In this paper, we adopt a different approach. We revisit the derivation of
hydrodynamics from the microscopic to the mesoscopic and macroscopic levels
and establish the corresponding thermodynamics at each of these levels, in the
light of the recent progresses in ST. In the first part of the paper we focus
on macroscopic hydrodynamics, where we decompose the EPR in terms of the
constraints imposed at the boundaries (e.g. a nonuniform temperature profile).
This allows us to identify i) the global thermodynamic forces that drive the
system out of equilibrium, and ii) the proper thermodynamic potential. When
the thermodynamic forces are switched off and no time-dependent protocol is
externally applied, the thermodynamic potential is maximized in the relaxation
process to equilibrium, for which it acts as a Lyapunov function of the
nonlinear dynamics. This shows how thermodynamic consistency (_i.e._ , the
requirement that the heat flux satisfies the second law) translates into a
dynamical result for a fully nonlinear theory. In the second part of the
paper, we focus on microscopic and mesoscopic descriptions behind macroscopic
hydrodynamics to reveal the information-theoretic origin of the thermodynamic
potential and EPR derived in the first part. To do so, we revisit the standard
derivation of hydrodynamics which starts from the microscale (hamiltonian
dynamics for particles interacting via pairwise potentials) and goes to the
macroscale via the mesoscale (Boltzmann equation) employing the Chapman–Enskog
expansion. For an isolated system, when considering an initially distribution
factorized in terms of identical one-particle distributions, we show that the
entropy production (identified with the change in Shannon entropy of the one
particle distribution) is always non-negative. Furthermore, in conceptual
analogy with the system-reservoir setup, it takes the form of a relative
entropy between the exact many-body probability distribution and the product
of the one-particle distributions. Entropy production is thus a direct
measures the information lost when neglecting inter-particle correlations.
When resorting to the molecular chaos hypothesis – used to derive a mesoscopic
description in terms of the Boltzmann equation which ensures relaxation of the
one-particle distribution to equilibrium – the EPR becomes non-negative.
Further resorting to the Chapman–Enskog expansion de Groot and Mazur (1962);
Kardar (2007) and the relaxation-time approximation Lundstrom (2002) used to
derive macroscopic hydrodynamics, we identify the parts of the EPR which
survive and give rise to the entropy production of irreversible
thermodynamics. We then move to non-isolated systems which we model at the
level of the Boltzmann equation by describing the collisions at nonisothermal
and possibly moving boundary in a thermodynamically consistent way inspired by
Refs. Bergmann and Lebowitz (1955); Lebowitz and Bergmann (1957); Cercignani
(1988); Gaspard (2022). We show that a first and second law ensue and that in
absence of thermodynamics forces, the nonequilibrium Massieu potential is
maximized by the relaxation dynamics towards equilibrium. Finally, we derive a
linear Boltzmann equation, similar to a master equation in the momentum space.
This mesoscale theory (alternative to the nonlinear Boltzmann equation) can be
seen as a ST underlying macroscopic hydrodynamics in the regime where the
Chapman–Enskog approach is justified.

The main outcome of our work is to unambiguously show how microscopic and
mesoscopic formulations of nonequilibrium thermodynamics in terms information-
theoretic quantities, are consistent with the macroscopic thermodynamics of
hydrodynamics, for both isolated and non-isolated systems.

### Outline

The paper is organized as follows. In Section II, we consider the dynamics of
continuous systems and their thermodynamics, write the balance equations for
energy and entropy, and then identify the proper thermodynamic potentials and
forces for isolated and non-isolated systems. In Section III, we derive
thermodynamic inequalities for the temporal change in information-theoretic
entropy and nonequilibrium Massieu potential at the microscopic scale
following hamiltonian dynamics by using a coarse-grained approach based on the
one-particle distribution function. In Section IV, we move to the mesoscopic
scale and revisit the thermodynamic aspects of the microscopic derivation of
the hydrodynamics of isolated systems using the Boltzmann equation. We use the
Chapman–Enskog expansion to obtain the local equilibrium condition as the
leading approximation in terms of the slowly varying variables describing
macroscopic systems. We furthermore estimate the error on the EPR during the
initial stages of the relaxation to the local equilibrium. In Section V, we
generalize the Boltzmann equation to describe non-isolated systems in a
thermodynamically consistent fashion. By writing explicitly the energy and
entropy balance equations, we show how the nonequilibrium Massieu potential
corresponds to the macroscopic thermodynamic potential identified in Section
II. Finally, in Section VI, we formulate a linear version of the generalized
Boltzmann that is consistent with both ST and the macroscopic hydrodynamics.
Future developments are discussed in Section VII.

##  II Macroscopic Scale: Navier-Stokes Equations

We revisit here the standard hydrodynamics of viscous fluids de Groot and
Mazur (1962); Kardar (2007). We start by obtaining the local formulation of
thermodynamics and then systematically show how to determine the thermodynamic
potential as well as the global thermodynamic forces keeping non-isolated
fluids out of equilibrium.

###  II.1 Dynamics

The dynamics of fluids is given by the Navier-Stokes equations de Groot and
Mazur (1962), written in lagrangian form using the material derivative
$D_{t}\equiv\partial_{t}+\bm{v}\cdot\nabla_{\bm{x}}$ as

$\displaystyle D_{t}\rho$ $\displaystyle=-\rho\nabla_{\bm{x}}\cdot\bm{v}\,,$
(1a) $\displaystyle\rho D_{t}\bm{v}$ $\displaystyle=-\nabla_{\bm{x}}\cdot
P+\rho\bm{F}\,,$ (1b)

for the spatio-temporal fields density of mass $\rho(\bm{x};t)$, velocity
$\bm{v}(\bm{x};t)$, external force (per unit mass) $\bm{F}(\bm{x};t)$, and
pressure tensor $P(\bm{x};t)$ (with $(\nabla_{\bm{x}}\cdot
P)_{j}=\sum_{i}\partial_{x_{i}}P_{ij}$, and $i$ and $j$ labeling the three
spatial components). The pressure tensor $P$ describes the stress exerted by
the fluid element on its surrounding and, in the case of isotropic,
compressible fluids, can be written in terms of the scalar pressure
$p(\bm{x};t)$ and the (symmetric) viscous tensor $\Pi(\bm{x};t)$. In turn, for
fluids sufficiently close to thermodynamic equilibrium $\Pi(\bm{x};t)$ is a
function of the derivatives of the velocity Landau and Lifshitz (1987):

$\displaystyle
P_{ij}(\bm{x};t)=p(\bm{x};t)\delta_{ij}-\Pi_{ij}(\bm{x};t)=p(\bm{x};t)\delta_{ij}-\left(\delta_{ij}\lambda\nabla_{\bm{x}}\cdot\bm{v}(\bm{x};t)+\nu\rho(\bm{x};t)\left(\partial_{x_{j}}v_{i}(\bm{x};t)+\partial_{x_{i}}v_{j}(\bm{x};t)\right)\right)\,.$
(2)

This expression of $P$ is obtained under a series of assumptions, valid in the
linear regime. First, $P$ is assumed to be invariant under global Galilean
transformations, _i.e._ ,
$\bm{v}({\bm{x}};t)\mapsto\bm{v}({\bm{x}};t)+\bm{V}$, implying that it depends
on the velocity only through its derivatives. Second, $P$ is assumed to depend
only on the first derivatives $\partial_{x_{i}}v_{j}$, which are supposed to
be small enough to contribute only linearly to $P$. Third, the pressure tensor
$P$ is assumed to be symmetric, which is justified when the molecules
composing the fluid are spherical or when the fluid density is very low de
Groot and Mazur (1962). Fourth, assuming isotropy, only two different scalar
coefficients, ${\nu}$ and $\lambda$, are needed to characterize viscous
effects in the fluid (out of six possible ones). The coefficient ${\nu}$,
accounting for the shear, is called _kinematic_ viscosity (while $\rho{\nu}$
is called _dynamic_ viscosity) and the coefficient $\lambda$, accounting for
dissipative compression or expansion, is called _second viscosity_. When it is
not specified otherwise, we will not rely on the specific expression of
$\Pi_{ij}$ in (2) in the rest of the article.

The mass conservation (1a) can equivalently be written in eulerian form as

$\displaystyle\partial_{t}\rho$
$\displaystyle=-\nabla_{\bm{x}}\cdot(\rho\bm{v})\,.$ (3)

Since we will often switch between the lagrangian and eulerian description, we
rely on this useful identity valid for any scalar quantity $o$

$\displaystyle\rho D_{t}o=\partial_{t}(\rho o)+\nabla_{\bm{x}}\cdot(\rho
o\bm{v})\,,$ (4)

which is derived using the definition of $D_{t}$ and Eq. (3).

Note that throughout this manuscript, we use a compact notation in which any
spatial and temporal dependence (of $\rho$ and $\bm{v}$, for instance) is left
implicit if no ambiguity can arise.

###  II.2 Local Energy Balance

The specific total energy, per unit mass,

$\displaystyle e=k+\phi+u\,,$ (5)

is given by three different contributions. The first,
$k(\bm{x};t)=\frac{1}{2}(v(\bm{x};t))^{2}$, is the specific kinetic energy
associated to the motion of the fluid element. The second, $\phi(\bm{x};t)$,
is the specific potential energy of the fluid element which can always be
written as the sum of two terms:
$\phi(\bm{x};t)=\phi_{\text{int}}(\bm{x};t)+\phi_{\text{w}}(\bm{x};t)$. The
term $\phi_{\text{int}}(\bm{x};t)$ accounts for the potential energy in the
bulk of the fluid. An example of this contribution is the gravitational
potential energy satisfying
$-\nabla_{\bm{x}}\phi_{\text{g}}(\bm{x})=-g\bm{u}=\bm{F}_{\text{g}}$ with
$\bm{F}_{\text{g}}$ the gravitational force (per unit mass), $g$ the
corresponding acceleration, and $\bm{u}$ the unit vector in the vertical
direction. The term $\phi_{\text{w}}(\bm{x};t)$ accounts instead for the
confinement generated by the (moving) boundaries. We assume that all possible
forces acting on a fluid element derive from a suitable potential energy, and
therefore $\bm{F}=-\nabla\phi$ in the following. The last contribution in Eq.
(5), $u(\bm{x};t)$, is the specific internal energy associated to all the
other degrees of freedom.

The balance equation for the total energy is $\rho
D_{t}e+\nabla_{\bm{x}}\cdot\bm{J}_{e}=0$, as long as the potential energy is
time-independent, while it becomes

$\displaystyle\rho
D_{t}e+\nabla_{\bm{x}}\cdot\bm{J}_{e}=\rho\partial_{t}\phi\,,$ (6)

when the potential energy is time-dependent. Here, $\bm{J}_{e}(\bm{x};t)$ is
the flux of total energy which can be specified by considering how each single
contribution to the total energy in Eq. (5) changes in time using the Navier-
Stokes equations (1).

First, we derive the balance equation for $k(\bm{x};t)$ from the momentum
conservation equation (1b) and Eq. (2):

$\displaystyle\rho D_{t}k+\sum_{i}\partial_{x_{i}}\left(\phi\rho
v_{i}+pv_{i}-\sum_{j}\Pi_{ij}v_{j}\right)=\phi\sum_{i}\partial_{x_{i}}(\rho
v_{i})+\sum_{i}p\partial_{x_{i}}v_{i}-\sum_{ij}\Pi_{ij}\partial_{x_{i}}v_{j}\,,$
(7)

where a term in the form of a divergence of fluxes has been isolated on the
left-hand side. On the right-hand side, source terms appear, even though at
this stage it is still not possible to tell apart dissipative contributions
from those that modify the potential energy of the system. Moving to the
eulerian frame using the identity (4), Eq. (7) becomes

$\displaystyle\partial_{t}(\rho
k)+\sum_{i}\partial_{x_{i}}\left(\rho(k+\phi)v_{i}+pv_{i}-\sum_{j}{\Pi_{ij}v_{j}}\right)=\phi\sum_{i}\partial_{x_{i}}(\rho
v_{i})+\sum_{i}p\partial_{x_{i}}v_{i}-\sum_{ij}\Pi_{ij}\partial_{x_{i}}v_{j}\,.$
(8)

Second, we consider the balance for the potential energy. By using Eq. (3), we
find

$\displaystyle\partial_{t}(\rho\phi)=\rho\partial_{t}\phi-\phi\nabla_{\bm{x}}\cdot(\rho\bm{v})=\rho\partial_{t}\phi-\phi\sum_{i}\partial_{x_{i}}(\rho
v_{i})\,.$ (9)

By adding Eqs. (8) and (9), we obtain the mechanical energy balance in the
eulerian frame:

$\displaystyle\partial_{t}(\rho(k+\phi))+\sum_{i}\partial_{x_{i}}\left(\rho(k+\phi)v_{i}+pv_{i}-\sum_{j}\Pi_{ij}v_{j}\right)=\rho\partial_{t}\phi+\sum_{i}p\partial_{x_{i}}v_{i}-\sum_{ij}\Pi_{ij}\partial_{x_{i}}v_{j}\,.$
(10)

Notice that the first term on the right-hand side of Eq. (8), representing the
nondissipative interconversion rate of kinetic energy into potential due to
compression and advection with the velocity field $\bm{v}$, is entirely
compensated by the corresponding term in the potential energy balance equation
(9) and, therefore, does not enter the balance equation for the mechanical
energy (10). In lagrangian form (obtained using (4)), Eq. (10) reads

$\displaystyle\rho
D_{t}(k+\phi)+\sum_{i}\partial_{x_{i}}\left(pv_{i}-\sum_{j}\Pi_{ij}v_{j}\right)=\rho\partial_{t}\phi+\sum_{i}p\partial_{x_{i}}v_{i}-\sum_{ij}\Pi_{ij}\partial_{x_{i}}v_{j}\,.$
(11)

Having derived the balance equation for the total mechanical energy, we use
Eq. (6) to obtain the balance equation for the internal energy $u(\bm{x};t)$:

$\displaystyle\rho D_{t}u=\rho\partial_{t}\phi-\rho
D_{t}(k+\phi)-\nabla_{\bm{x}}\cdot\bm{J}_{e}\,.$ (12)

Notice that the first term on the right-hand side, due to a time-dependent
change in potential energy, does not directly contribute to the balance
equation of the internal energy because it is compensated by the corresponding
term in the mechanical energy balance (11). Nevertheless, it contributes to
the internal energy since it enters the Navier-Stokes equations (via Eq. (1b))
and thus affects the state of the system.

The expression of total energy flux $\bm{J}_{e}$ can now be given as the sum
of the fluxes on the left-hand side of Eq. (11) and the, for the moment
unspecified, heat flux $\bm{J}_{q}$

$\displaystyle\bm{J}_{e}=p\bm{v}-{\Pi\bm{v}}+\bm{J}_{q}\,.$ (13)

Indeed, this leads to a (local) balance equation for the internal energy,
obtained using Eq. (11) and (13) in Eq. (12),

$\displaystyle\rho
D_{t}u+\nabla_{\bm{x}}\cdot\bm{J}_{q}=-p\nabla_{\bm{x}}\cdot\bm{v}+\operatorname{Tr}\Pi\nabla_{\bm{x}}\bm{v}\,,$
(14)

where the divergence only acts on the heat flux, while a source term appears
on the right-hand side. Here, we define
$\operatorname{Tr}\Pi\nabla_{\bm{x}}\bm{v}\equiv\sum_{ji}\Pi_{ij}\partial_{x_{j}}v_{i}$.
In eulerian form, Eq. (14) becomes

$\displaystyle\partial_{t}(\rho u)+\nabla_{\bm{x}}\cdot(\rho
u\bm{v}+\bm{J}_{q})=-p\nabla_{\bm{x}}\cdot\bm{v}+\operatorname{Tr}\Pi\nabla_{\bm{x}}\bm{v}\,.$
(15)

Equation (14) and (15) constitute two equivalent (local) formulations of the
first law of thermodynamics for fluids. Indeed, once they are integrated over
a volume element, one can recognize that the internal energy changes because
of the exchanged heat at the boundary, or the mechanical expansion of the
system, (accounted on the left-hand side of Eq. (15)) or the viscous and other
dissipative effects that take place inside the volume (accounted on the right-
hand side of Eq. (15)).

Finally, by substituting now the expression of the total energy flux (13) in
the balance equation (6) for the total energy, we obtain

$\displaystyle\rho
D_{t}e=\rho\partial_{t}\phi-\nabla_{\bm{x}}\cdot\bm{J}_{q}-\nabla_{\bm{x}}\cdot\left(p\bm{v}-\Pi\bm{v}\right)\,,$
(16)

which will be essential to prove the existence of a global Lyapunov function
in systems with a well defined equilibrium state (§ II.5).

###  II.3 Local Equilibrium and Local Entropy Balance

The second balance equation that the system must satisfy is the entropy
balance. To derive it in the macroscopic theory, we must resort to the
following assumption: each fluid element $\omega$ is locally equilibrated with
the surrounding, after a very short time scale corresponding to microscopic
collisions. This local equilibrium has well-defined values of local
temperature $T$, pressure $p$, and chemical potential $\mu$, and the fluid
element is thus characterized by the extensive quantity entropy $S_{\omega}$,
internal energy $U_{\omega}$, volume $V_{\omega}$, and number of particles
$N_{\omega}$ which can all vary in time. Furthermore, the particle density
inside the volume is uniform. Note that these conditions are satisfied if the
fluid is not extremely far from equilibrium. Thus, the fundamental relation of
equilibrium thermodynamics holds for the fluid element and can be written in
the form

$\displaystyle T\mathop{\mathrm{\mathstrut d}}\nolimits\\!_{t}S_{\omega}$
$\displaystyle=\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}U_{\omega}+p\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}V_{\omega}-M\mu\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}N_{\omega}\,,$ (17)

also called local equilibrium condition. Note that here $M$ is the molecular
mass and we used the chemical potential per unit mass $\mu=(\partial
S_{\omega}/\partial N_{\omega})/M$.

To formulate Eq. (17) in the same way as the balance equation (14) for the
internal energy, we need to introduce the specific entropy $s$, the specific
internal energy $u$, and the specific volume $\rho^{-1}$ (that coincides with
the inverse of the local mass density). This, together with Reynolds’ theorem
(see appendix A) and the identity (4), leads to

$\displaystyle\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}S_{\omega}=\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\int_{\omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}\,\rho s=\int_{\omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}\,\left(\partial_{t}(\rho
s)+\nabla_{\bm{x}}\cdot\left(\rho
s\bm{v}\right)\right)=\int_{\omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}\,\rho D_{t}s\,,$ (18a)
$\displaystyle\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}U_{\omega}=\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\int_{\omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}\,\rho u=\int_{\omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}\,\left(\partial_{t}(\rho
u)+\nabla_{\bm{x}}\cdot\left(\rho
u\bm{v}\right)\right)=\int_{\omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}\,\rho D_{t}u\,,$ (18b)
$\displaystyle\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}V_{\omega}=\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\int_{\omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}\,\rho\rho^{-1}=\int_{\omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}\,(\partial_{t}(\rho\rho^{-1})+\nabla_{\bm{x}}\cdot(\rho\rho^{-1}\bm{v}))=\int_{\omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}\,\rho D_{t}\rho^{-1}\,,$ (18c)
$\displaystyle\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}N_{\omega}=\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\int_{\omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}\,\rho
M^{-1}=\int_{\omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}\,(\partial_{t}(\rho M^{-1})+\nabla_{\bm{x}}\cdot(\rho
M^{-1}\bm{v}))=\int_{\omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}\,\rho D_{t}M^{-1}=0\,,$ (18d)

where Eq. (18d) means that particles do not leave or enter the fluid element,
which moves according to the velocity field $\bm{v}$. The local equilibrium
condition (17) in terms of specific quantities becomes

$\displaystyle\int_{\omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}\rho\left(TD_{t}s-D_{t}u-pD_{t}\rho^{-1}\right)$
$\displaystyle=0\,.$ (19)

Since all variables are assumed to be uniform inside the fluid element
$\omega$, we can also write the specific variables as
$s=S_{\omega}/(MN_{\omega})$, $u=U_{\omega}/(MN_{\omega})$,
$\rho^{-1}=V_{\omega}/(MN_{\omega})$, with $MN_{\omega}$ being the total mass
of the fluid element. Equation (19) thus becomes

$\displaystyle TD_{t}s=D_{t}u+pD_{t}\rho^{-1}\,.$ (20)

The extensivity of the entropy, _i.e._ , $S_{\omega}(\lambda
U_{\omega},\lambda V_{\omega},{\lambda N_{\omega}})=\lambda
S_{\omega}(U_{\omega},V_{\omega},{N_{\omega}})$ for all positive values of
$\lambda$, can be combined with Eq. (20) to provide a constraint on the
intensive fields $T$, $p$, $\mu$, which are therefore not independent from
each other. Indeed, by taking the derivative of $S_{\omega}(\lambda
U_{\omega},\lambda V_{\omega},{\lambda N_{\omega}})$ in $\lambda$ evaluated at
$\lambda=1$ and dividing everything by the total mass $MN_{\omega}$, we get
the Euler relation

$\displaystyle s=\frac{u}{T}+\frac{p}{T}\rho^{-1}-\frac{\mu}{T}\,.$ (21)

Then, by subtracting Eq. (20) from the material derivative of $s$ in Eq. (21),
we obtain the Gibbs-Duhem equation, that constrains the derivatives of
intensive quantities

$\displaystyle
uD_{t}\left(\frac{1}{T}\right)+\rho^{-1}D_{t}\left(\frac{p}{T}\right)-D_{t}\left(\frac{\mu}{T}\right)=0\,.$
(22)

Finally, the (local) entropy balance is obtained by using the first law (14)
and the continuity equation (1a) in the local equilibrium (20):

$\displaystyle T\rho
D_{t}s=-\nabla_{\bm{x}}\cdot\bm{J}_{q}+\operatorname{Tr}\Pi\nabla_{\bm{x}}\bm{v}\,,$
(23)

or, equivalently,

$\displaystyle\rho
D_{t}s+\nabla_{\bm{x}}\cdot\left(\frac{\bm{J}_{q}}{T}\right)=\dot{\sigma}_{\text{macro}}\,.$
(24)

In Eq. (24) we identify the entropy flux $\bm{J}_{s}\equiv\bm{J}_{q}/T$, and
the (macroscopic) entropy production rate per unit volume

$\displaystyle\dot{\sigma}_{\text{macro}}$
$\displaystyle=\bm{J}_{q}\cdot\nabla_{\bm{x}}\left(\frac{1}{T}\right)+\frac{1}{T}\operatorname{Tr}\Pi\nabla_{\bm{x}}\bm{v}$
(25a) $\displaystyle=\dot{\sigma}_{\text{hc}}+\dot{\sigma}_{\text{vf}}\,.$
(25b)

The entropy production rate is given by the sum of two different contributions
corresponding to different local dissipative mechanisms: heat conduction
($\dot{\sigma}_{\text{hc}}$) and viscous friction
($\dot{\sigma}_{\text{vf}}$). In the macroscopic theory, the assumption that
the two dissipative contributions in Eq. (25b) are separately positive
constitutes a formulation of the second law of thermodynamics. In the
microscopic and mesoscopic theory, one is instead interested in showing the
positivity of the EPR combining hypothesis concerning the dynamics and the
statistical behavior of a large number of particles (see Section III, IV, and
V).

In the linear regime, $\bm{J}_{q}$ and $\Pi$ are linear functions of
$\nabla_{\bm{x}}T$ and $\nabla_{\bm{x}}\bm{v}$. This implies that Eqs. (1a),
(1b) and (14) form a closed dynamical system for the spatio-temporal fields
$(\rho,u,\bm{v})$ when the intensive fields $p$, $T$, $\mu$ are also expressed
as functions of $(\rho,u,\bm{v})$ from the knowledge of the entropy appearing
in (21), which specifies all the needed equations of state of the system.

However, global thermodynamics provides information on the dynamics also
beyond the linear regime, _e.g._ , the existence of a Lyapunov function for
systems in contact with an equilibrium environment. In the following, we will
illustrate this point and we will also explain how different conditions
imposed at the boundary (see in Fig. 1) determine different thermodynamics
potentials, and give rise to different global thermodynamic forces.

###  II.4 Entropy Balance and Thermodynamic Potential for Isolated Systems

We formulate now the global entropy balance for an isolated fluid contained in
a domain $\Omega$ and determine its thermodynamic potential. We label
$\partial\Omega$ and $\bm{n}$ the boundary and the outward pointing unit
vector, respectively. On this boundary, the heat flow and the velocity field
must satisfy

$\displaystyle\left.\bm{J}_{q}({\bm{x}};t)\cdot\bm{n}\right|_{\partial\Omega}$
$\displaystyle=0\,,$ (26a)
$\displaystyle\left.\bm{v}({\bm{x}};t)\cdot\bm{n}\right|_{\partial\Omega}$
$\displaystyle=0\,,$ (26b)

namely, heat flux and velocity are orthogonal to the boundary which physically
means that no energy or momentum can be exchanged with the boundaries. Note
that Eq. (26b) is trivially satisfied if one imposes a no-slip condition at
the boundary, _i.e._ ,
$\left.\bm{v}({\bm{x}};t)\right|_{\partial\Omega}=\bm{0}$. We now take the
time derivative of the global entropy $S\equiv\int_{\Omega}\rho
s\mathop{\mathrm{\mathstrut d}}\nolimits\\!{\bm{x}}$ using the local entropy
balance (24), Eq. (4) and the boundary conditions in Eq. (26):

$\displaystyle\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}S=\int_{\Omega}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}\,\partial_{t}(\rho s)$
$\displaystyle=-\int_{\partial\Omega}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{n}\cdot\left(\frac{\bm{J}_{q}}{T}+\rho
s\bm{v}\right)+\int_{\Omega}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}\,\dot{\sigma}_{\text{macro}}=\int_{\Omega}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}\,\dot{\sigma}_{\text{macro}}\equiv\dot{\Sigma}_{\text{macro}}\geq
0\,,$ (27)

where we introduced the global entropy production rate
$\dot{\Sigma}_{\text{macro}}$. The inequality in Eq. (27), imposed by the
positivity of $\dot{\sigma}_{\text{macro}}$, implies that the entropy is the
thermodynamic potential that is maximized when an isolated system reaches
equilibrium Callen (1948). Prigogine introduced the notion of _evolution
criterion_ for an inequality of the type in Eq. (27), where the time
derivative of a thermodynamic potential has a definite sign Glansdorff  _et
al._ (1974). If such a potential has an upper bound, then it plays the role of
the (negative of the) Lyapunov function for the system. In the next subsection
II.5, we obtain analogous relations for non-isolated systems.

###  II.5 Entropy Balance and Thermodynamic Potential for Non-Isolated Systems

Figure 1:  Illustration of a fluid in contact with regions of the environment
on which different kinds of boundary conditions are imposed. The boundary
surface is partitioned in 6 different regions corresponding to different types
of boundary conditions. Regions 1-3 correspond to movable pistons coupling the
system to reservoirs kept at equilibrium at a fixed temperature and pressure,
described by Eq. (46). Regions 4-5 are characterized by an externally fixed
velocity on the boundary, which results in shear (4) or volume changes if a
movable wall is involved (5). The thermodynamics of a system with this kind of
boundary conditions is given in Eq. (40). In region 6 the wall is not moving
and a heat flux can be maintained using external thermal devices, _e.g._
resistors coupled to an electric circuit. Notice that the external temperature
should be defined everywhere on the boundary to allow for equilibration, which
occurs when all the regions of the boundary have the same temperature.

In non-isolated systems, the existence of a Lyapunov function satisfying the
equivalent of Eq. (27) (an evolution criterion in Glansdorff and Prigogine
terminology Glansdorff and Prigogine (1971)) depends on the boundary
conditions. We now show how the thermodynamic potential which is maximized at
equilibrium is modified when different experimentally relevant boundary
conditions are imposed.

####  II.5.1 Closed, rigid systems.

Consider a closed but non-isolated system in a rigid domain $\Omega$ with
constant volume $V$, delimited by a boundary $\partial\Omega$ satisfying

$\displaystyle\left.T({\bm{x}};t)\right|_{\partial\Omega}$
$\displaystyle=T_{B}\,,$ (28a)
$\displaystyle\left.\bm{v}({\bm{x}};t)\right|_{\partial\Omega}$
$\displaystyle=\bm{0}\,,$ (28b)

namely, with fixed temperature and a no-slip condition. The system can
exchange heat with the environment via a nonvanishing heat flux $\bm{J}_{q}$
at the boundary. If the potential $\phi$ as well as the temperature $T_{B}$
are constant in time, the system will relax to the equililbrium compatible
with the external temperature. However, the global entropy
$S=\int_{\Omega}\rho s\mathop{\mathrm{\mathstrut d}}\nolimits\\!{\bm{x}}$ is
not the thermodynamic potential anymore. Indeed, by using the local entropy
balance equation (24), Eq. (4) and the boundary conditions in Eq. (28), we
find

$\displaystyle\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}S=-\int_{\partial\Omega}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{n}\cdot\left(\frac{\bm{J}_{q}}{T}+\rho
s\bm{v}\right)+\int_{\Omega}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}\,\dot{\sigma}_{\text{macro}}=\frac{{\dot{Q}}}{T_{B}}+\dot{\Sigma}_{\text{macro}}$
(29)

where ${\dot{Q}}\equiv-\int_{\partial\Omega}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{n}\cdot\bm{J}_{q}$ is the exchanged heat at the boundary
(positive when it flows from the environment to the system) which is related
to the global total energy $E\equiv\int_{\Omega}\rho
e\mathop{\mathrm{\mathstrut d}}\nolimits\\!{\bm{x}}$ via

$\displaystyle\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}E=-\int_{\partial\Omega}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{n}\cdot\left(\bm{J}_{q}+p\bm{v}-\Pi\bm{v}+\rho
e\bm{v}\right)={\dot{Q}}\,,$ (30)

using Eqs. (16), (4), and (28). In the macroscopic theory, the second law of
thermodynamics corresponds to the assumption that $\Sigma_{\text{macro}}\geq
0$ (which will be discussed from the microscopic viewpoint starting in Section
III. Assuming the validity of the second law, the thermodynamic potential that
is maximized at equilibrium can be written as

$\displaystyle\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\left(S-\frac{E}{T_{B}}\right)=\dot{\Sigma}_{\text{macro}}\geq
0\,.$ (31)

This means that for a closed, rigid system in contact with an external
thermostat, the thermodynamic potential is the difference between the global
entropy and the global total energy divided by the temperature of the
boundary.

####  II.5.2 Closed system under pressure control.

Closed systems contained in either a flexible membrane with negligible surface
tension or in a vessel with a movable, weightless piston in contact with a gas
in equilibrium at a given temperature and pressure represent idealized
examples where the volume can change in time $V(t)$. Mathematically, the
temperature on the boundary is externally fixed according to

$\displaystyle\left.T({\bm{x}};t)\right|_{\partial\Omega(t)}$
$\displaystyle=T_{B}\,\,.$ (32)

On the other hand, the pressure $p({\bm{x}};t)$ of the fluid and that of the
environment $p_{B}({\bm{x}})$ do not have to coincide on the boundary. This
has two implications. First, the imbalance between the two pressures generates
the external force in the Navier-Stokes equation (1b) on the boundary. Second,
the boundary moves with a velocity $\bm{v}_{B}(\bm{x})$ that coincides with
the velocity field $\bm{v}(\bm{x})$ of the fluid on the boundary

$\displaystyle\left.(\bm{v}(\bm{x})-\bm{v}_{B}(\bm{x}))\right|_{\partial\Omega(t)}=\bm{0}\,,$
(33)

because of the no-slip condition. Hence, the volume changes in time according
to

$\displaystyle\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}V=\int_{\partial\Omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{n}\cdot\bm{v}_{B}(\bm{x})\,\,,$ (34)

as prescribed by Reynolds’ theorem, and $\phi$ changes in time according to

$\int_{\Omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{x}\,\rho(\bm{x})\partial_{t}\phi(\bm{x})=\int_{\partial\Omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{n}\cdot(P(\bm{x})-p_{B}(\bm{x})\mathbb{1})\bm{v}_{B}(\bm{x})$
(35)

where $\mathbb{1}$ is the identity matrix and $(P-p_{B}\mathbb{1})$ is
(proportional to) the force on the boundary. Equation (35) is derived by
considering that
$\partial_{t}\phi(\bm{x})=\partial_{t}\phi_{\text{w}}(\bm{x}-\bm{x}_{B})$ and
$\partial_{t}\phi_{\text{w}}(\bm{x}-\bm{x}_{B})$ changes in time because of
the motion of the boundary $\bm{x}_{B}$ at velocity $\bm{v}_{B}$:
$\partial_{t}\phi_{\text{w}}(\bm{x}-\bm{x}_{B})=-\nabla_{\bm{x}}\phi_{\text{w}}(\bm{x})\cdot\bm{v}_{B}(\bm{x})$.

By using Eqs. (16), (24) and Eq. (35), we can write the balance equations for
the total energy $E=\int_{\Omega}\rho e\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}$, and entropy $S=\int_{\Omega}\rho
s\mathop{\mathrm{\mathstrut d}}\nolimits\\!{\bm{x}}$ which read

$\displaystyle\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}E=\int_{\Omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}\,\rho
D_{t}e=\dot{Q}-\int_{\partial\Omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{n}\cdot p_{B}\bm{v}_{B}\,,$ (36a)
$\displaystyle\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}S=\int_{\Omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}\,\rho
D_{t}s=\frac{\dot{Q}}{T_{B}}+\dot{\Sigma}_{\text{macro}}\,,$ (36b)

respectively. By combining now Eq. (36a) with Eq. (36b), we obtain the
following expression for the global entropy production rate

$\displaystyle\dot{\Sigma}_{\text{macro}}$
$\displaystyle=\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\left(S-\frac{E}{T_{B}}\right)-\frac{1}{T_{B}}\int_{\partial\Omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{n}\cdot p_{B}\bm{v}_{B}\geq 0\,.$ (37)

If the pressure on the boundary is constant, namely,
$\int_{\partial\Omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{n}\cdot(p_{B}\mathbb{1})\bm{v}_{B}=p_{B}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}V$ (from Eq. (34)), Eq. (37) becomes

$\displaystyle\dot{\Sigma}_{\text{macro}}$
$\displaystyle=\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\left(S-\frac{E+p_{B}V}{T_{B}}\right)\geq 0\,,$ (38)

which allows us to identify the proper thermodynamic potential of a closed
systems in contact with an external thermostat and exposed to a constant
pressure as the difference between the global entropy and the sum of the
global energy plus the $p_{B}V$ divided by the temperature of the boundary.

For a reversible protocol, the entropy production rate vanishes,
$\dot{\Sigma}_{\text{macro}}=0$, and Eq. (37) identifies the maximum
mechanical work which can be extracted by letting the system relax to
equilibrium as Landau and Lifshitz (2013)

$\displaystyle W^{\text{max}}_{\text{mech}}\equiv
p_{B}(V_{\text{eq}}-V(0))=E(0)-E_{\text{eq}}-T_{B}(S_{\text{eq}}-S(0))\,.$
(39)

This maximum available work or available free energy is also known as _exergy_
in engineering Bejan (2016); Kotas (2012).

####  II.5.3 Closed system under velocity control.

The volume of closed systems can also change when the velocity profile on the
boundary $\bm{v}_{B}$ is externally controlled according to Eq. (34). If we
consider again a uniform temperature on the boundary as in Eq. (32), we can
combine the balance equations for the total energy $E=\int_{\Omega}\rho
e\mathop{\mathrm{\mathstrut d}}\nolimits\\!{\bm{x}}$ and and entropy
$S=\int_{\Omega}\rho s\mathop{\mathrm{\mathstrut d}}\nolimits\\!{\bm{x}}$ in
Eqs. (16), and (24), respectively, to express the global entropy production as

$\dot{\Sigma}_{\text{macro}}=\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\left(S-\frac{E}{T_{B}}\right)+\frac{1}{T_{B}}\underbrace{\int_{\Omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{x}\,(\rho\partial_{t}\phi-\nabla_{\bm{x}}\cdot
P\bm{v})}_{=\dot{W}_{\text{mech}}}\geq 0\,,\\\ $ (40)

where we identify $\dot{W}_{\text{mech}}$ as the mechanical work (per unit
time). Here, as well as in Subs. II.5.2, the potential energy density is time-
dependent becuase of the moving boundary.

If the velocity profile at the boundary is nonuniform, the system remains out
of equilibrium and undergoes either expansion, or compression, or shearing.
From a thermodynamic standpoint, this results from the environment
continuously performing mechanical work. Furthermore, if this mechanical
action is constant in time and the system is far from a turbulent regime
(_i.e._ when the nonlinear terms in (1b) are negligible), the system will
approach a nonequilibrium steady state. If the velocity profile at the
boundary is uniform, the mechanical work in Eq. (40) vanishes (as a result of
Eq. (40) being invariant under global Galilean transformations, see § II.1)
and the system relaxes towards equilibrium.

####  II.5.4 Closed system with multiple heat baths and time-dependent
driving.

We now consider a closed non-isolated system contained in a domain $\Omega$ of
volume $V$ whose boundary $\partial\Omega$ is composed of $b=1,\dots,B$
distinct regions $\partial\Omega_{{b}}$ called reservoirs, or baths, i.e.,
$\partial\Omega=\cup_{b}\partial\Omega_{{b}}$. The temperature of each region
$T_{{b}}(t)$ and the (uniform) pressure each region is exposed to $p_{b}(t)$
are externally controlled and may depend on time following arbitrary
protocols. On each region $b$, the fluid temperature is assumed to be
equilibrated with the environment, namely,

$\displaystyle\left.T({\bm{x}};t)\right|_{\partial\Omega_{{b}}}$
$\displaystyle=T_{{b}}(t)\,,$ (41)

while its pressure differs in general from the external pressure, i.e.,
$P(x;t)|_{\partial\Omega_{{b}}}\neq p_{b}(t)$. As we already discussed in
Subs. II.5.2, the mismatch of pressure at each region implies that i) the
boundaries move according to Eq. (33) (which now holds at the level of each
single region $b$), ii) the volume changes according to Eq. (34) (which can be
split into the contribution of each single region), and iii) the potential
energy density changes in time according to Eq. (35) (which now holds at the
level of each single region $b$).

From a thermodynamic standpoint, by using Eqs. (16), (24) and Eq. (35)
together with $\partial\Omega=\cup_{b}\partial\Omega_{{b}}$, the balance
equations for the total energy $E=\int_{\Omega}\rho
e\mathop{\mathrm{\mathstrut d}}\nolimits\\!{\bm{x}}$, and entropy
$S=\int_{\Omega}\rho s\mathop{\mathrm{\mathstrut d}}\nolimits\\!{\bm{x}}$ read

$\displaystyle\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}E=\sum_{b=1}^{B}\dot{Q}_{b}-\sum_{b=1}^{B}p_{b}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}V_{b}$ (42a) $\displaystyle\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}S=\sum_{b=1}^{B}\frac{\dot{Q}_{b}}{T_{b}}+\dot{\Sigma}_{\text{macro}}\,,$
(42b)

respectively, with
$\dot{Q}_{{b}}\equiv-\int_{\partial\Omega_{{{b}}}}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{n}\cdot\bm{J}_{q}$ and $\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}V_{b}\equiv\int_{\partial\Omega_{b}(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{n}\cdot\bm{v}_{b}$.

Without loss of generality, we choose the reservoir $b=B$ as the reference for
the equilibrium temperature and pressure. Thus, by combining now Eq. (42a)
with Eq. (42b), the global entropy production rate can be written as

$\dot{\Sigma}_{\text{macro}}=\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}S-\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}E}{T_{B}}+\sum_{b=1}^{B}\dot{Q}_{b}\bigg{(}\frac{1}{T_{B}}-\frac{1}{T_{b}}\bigg{)}-\frac{1}{T_{B}}\sum_{b=1}^{B}p_{b}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}V_{b}\geq 0\,,$ (43)

or, by summing and subtracting $(p_{B}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}V)/T_{B}$ (with $\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}V=\sum_{b}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}V_{b}$) together with collecting all the time derivatives
in a single contribution, as

$\dot{\Sigma}_{\text{macro}}=\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\bigg{(}S-\frac{E+p_{B}V}{T_{B}}\bigg{)}+\sum_{b=1}^{B}\dot{Q}_{b}\bigg{(}\frac{1}{T_{B}}-\frac{1}{T_{b}}\bigg{)}+\frac{1}{T_{B}}\sum_{b=1}^{B}(p_{B}-p_{b})\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}V_{b}+(E+p_{B}V)\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\frac{1}{T_{B}}+\frac{V}{T_{B}}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}p_{B}\geq 0\,.$ (44)

Equation (44) is the general balance equation for the global EPR in terms of
(from the right to the left) i) time-dependent driving of the reference
pressure $p_{B}$ and temperature $T_{B}$, ii) global thermodynamic forces
(_i.e._ , pressure and temperature differences between the boundaries), and
iii) the time derivative of a thermodynamic potential,

$Y(t)\equiv S-\frac{E+p_{B}(t)V(t)}{T_{B}(t)}\,.$ (45)

If the system was coupled to a single reservoir with no time-dependent driving
of $T_{B}$ and $p_{B}$, then the thermodynamic potential $Y(t)$ would be
maximized.

Note that, by using the thermodynamic potential (45), Equation (44) reads

$\dot{\Sigma}_{\text{macro}}=\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}Y+\sum_{b=1}^{B}\dot{Q}_{b}\bigg{(}\frac{1}{T_{B}}-\frac{1}{T_{b}}\bigg{)}+\frac{1}{T_{B}}\sum_{b=1}^{B}(p_{B}-p_{b})\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}V_{b}-\frac{\partial
Y}{\partial(1/T_{B})}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\frac{1}{T_{B}}-\frac{\partial Y}{\partial
p_{B}}\mathop{\mathrm{\mathstrut d}}\nolimits\\!_{t}p_{B}\geq 0\,.$ (46)

The systematic derivation of the decomposition of the EPR in terms of the
global thermodynamic forces, the time-dependent driving and the derivative of
a thermodynamic potential is a significant result of the present work. Indeed,
as it is well known, the classic formulation of nonequilibrium thermodynamics
Prigogine (1949); de Groot and Mazur (1962) is unable to differentiate systems
that relax toward equilibrium from those who are actively maintained away from
it by just looking at the expression of the EPR, in which only local
thermodynamic forces appear Gallavotti (2002). The systematic procedure to
obtain the decomposition of the EPR that we outlined here for the case of a
single component fluid, but which can be easily extended to more complex
situations, is strictly related to the one used in engineering thermodynamics
Bejan (2016); Kotas (2012) in the context of the so-called “exergy method” of
analysis. However, the latter is introduced in a phenomenological way that may
disregard some contributions, like those accounting for time-dependent driving
on the boundaries. Equation (46) may thus provide a useful point of view to
interpret new numerical and laboratory experiments in fluid mechanics Yang
_et al._ (2020); Urban  _et al._ (2022) in which time-dependent driving at the
boundaries translates into nontrivial effects on macroscopic heat fluxes
affecting the efficiency of heat transport.

Notice that an alternative formulation of Eq. (46) can be obtained by
representing each reservoir $b=1,\dots,B$ as an ideal gas at equilibrium at
temperature $T_{b}$. Indeed, this allows us to quantify the energy exchanges
between the system and the $b$ reservoir according to the equilibrium first
law:

$\displaystyle\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}E_{b}=\dot{Q}_{b}-p_{b}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}V_{b}\,,$ (47)

where $-\mathop{\mathrm{\mathstrut d}}\nolimits\\!_{t}E_{b}$ is the energy
variation of the $b$ reservoir. Thus, by using Eq. (47) in Eq. (46) and
treating $1/T_{B}$ and $p_{B}/T_{B}$ as independent intensive variables of the
thermodynamic potential $Y(t)$ in Eq. (45), Eq. (46) becomes

$\small\dot{\Sigma}_{\text{macro}}=\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}Y+\sum_{b=1}^{B}\left[\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}{E}_{{{b}}}\left(\frac{1}{T_{B}}-\frac{1}{T_{{b}}}\right)+\left(\frac{p_{B}}{T_{B}}-\frac{p_{b}}{T_{b}}\right)\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}V_{b}\right]-\frac{\partial
Y}{\partial(1/T_{B})}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\frac{1}{T_{B}}-\frac{\partial
Y}{\partial(p_{B}/T_{B})}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\left(\frac{p_{B}}{T_{B}}\right)\geq 0\,.$ (48)

Eq. (48) gives the rationale for the expressions of the intensive fields used
in stochastic thermodynamics, see _e.g._ Rao and Esposito (2018).

###  II.6 Linear regime and minimum entropy production principle

We conclude the discussion on the macroscopic aspects of the thermodynamics of
non-isolated systems by showing that systems close to equilibrium satisfy a
minimum entropy production principle in a form similar to the one obtained in
the framework of stochastic thermodynamics Forastiere  _et al._ (2022). In
turn, the latter is akin to the adiabatic-nonadiabatic decomposition Esposito
and Van den Broeck (2010b) specialized in the linear regime. For related work
in the macroscopic framework see also ref. Maes and Netočnỳ (2015). To obtain
the aformentioned decomposition of the EPR, we write the macroscopic entropy
production rate $\dot{\Sigma}_{\text{macro}}$ obtained by integrating Eq.
(25a) in the linear regime, namely, when
$\bm{J}_{\text{q}}=\kappa\nabla_{\bm{x}}(1/T)$ and $\Pi$ is evaluated at the
equilibrium density $\rho_{\text{eq}}$ and temperature $T_{\text{eq}}$:

$\displaystyle\dot{\Sigma}_{\text{macro}}\simeq\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}\begin{pmatrix}\partial_{x_{i}}\left(1/T\right)\\\
\epsilon_{ij}\end{pmatrix}\underbrace{\begin{pmatrix}\kappa\,\delta_{ij}&0\\\
0&L_{(ij),(i^{\prime}j^{\prime})}\end{pmatrix}}_{\equiv\mathbb{O}}\begin{pmatrix}\partial_{x_{j}}\left(1/T\right)\\\
\epsilon_{i^{\prime}j^{\prime}}\end{pmatrix}\,.$ (49)

Here, we use a notation such that the vectors’ entries are the three
components of $\partial_{x_{i}}\left(1/T\right)$ and the nine components of
the strain tensor
$\epsilon_{ij}\equiv{\frac{1}{2}}(\partial_{x_{i}}v_{j}+\partial_{x_{j}}v_{i})$.
Furthermore, we introduce the $9\times 9$ positive-definite matrix whose
entries are
$L_{(ij),(i^{\prime}j^{\prime})}\equiv(\lambda\delta_{ij}\delta_{i^{\prime}j^{\prime}}+2\rho_{\text{eq}}\nu\delta_{ii^{\prime}}\delta_{jj^{\prime}})/{T_{\text{eq}}}$
which together with the $3\times 3$ positive-definite matrix with entries
$\kappa\delta_{ij}$ defines the positive-definite linear response matrix
$\mathbb{O}$. In Appendix B, we prove that Eq. (49) can be alternatively
written in terms of the steady state entropy production rate
$\dot{\Sigma}^{\text{ss}}_{\text{macro}}$ as

$\displaystyle\dot{\Sigma}_{\text{macro}}$
$\displaystyle\simeq\dot{\Sigma}^{\text{ss}}_{\text{macro}}+\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!{\bm{x}}\begin{pmatrix}\partial_{x_{i}}\left(1/T-1/T^{\text{ss}}\right)\\\
\epsilon_{ij}-\epsilon_{ij}^{\text{ss}}\end{pmatrix}\begin{pmatrix}\kappa\,\delta_{ij}&0\\\
0&L_{(ij),(i^{\prime}j^{\prime})}\end{pmatrix}\begin{pmatrix}\partial_{x_{j}}\left(1/T-1/T^{\text{ss}}\right)\\\
\epsilon_{i^{\prime}j^{\prime}}-\epsilon_{i^{\prime}j^{\prime}}^{\text{ss}}\end{pmatrix}\,.$
(50)

Thus, since the response matrix $\mathbb{O}$ is positive-definite, the state
that minimizes the total entropy production rate coincides with the steady-
state solution. In Appendix B we also prove that the EPR is a Lyapunov
function of the dynamics close to equilibrium (without resorting to further
assumptions). We remark that $\dot{\Sigma}^{\text{ss}}_{\text{macro}}$ can be
computed in terms of quantities only defined at the boundaries using for
example Eq. (46) (for systems without time-dependent driving)

$\displaystyle\dot{\Sigma}_{\text{macro}}^{\text{ss}}=\sum_{b=1}^{B-1}{\dot{Q}}_{{b}}^{\text{ss}}\left(\frac{1}{T_{B}^{\text{ss}}}-\frac{1}{T_{{b}}^{{\text{ss}}}}\right)\,.$
(51)

We further remark that the formulation (50) of the minimum entropy production
theorem is different from the statement proved in Ref. Glansdorff and
Prigogine (1971) and similar in spirit, but slightly stronger, than the one
contained in Ref. de Groot and Mazur (1962).

Equation (50) is a remarkable property of close to equilibrium steady states,
which has a direct analogous in stochastic thermodynamics, and it can be
interpreted as a simple geometric formulation of Prigogine’s minimum entropy
production principle Glansdorff and Prigogine (1971). However, the positivity
of the quadratic form in the local thermodynamic forces appearing on the
right-hand side of Eq. (50) is only ensured by the equilibrium properties of
the response coefficients, and thus cannot be generalized to situations in
which the steady state is not close to equilibrium.

##  III Microscopic Scale: Hamiltonian Dynamics

In this section, we derive thermodynamic identities (and the related
inequalities) for isolated systems made of particles obeying Hamiltonian
dynamics. We use the symbol $\bm{x}_{n}$ (resp. $\bm{p}_{n}$) for the spatial
coordinates (resp. momentum) of the $n$-th microscopic particle, while the
symbol ${\bm{x}}$ (resp. $\bm{p}$) is reserved to indicate a generic value of
position (resp. momentum).

###  III.1 Dynamics

We start from a system of $N$ identical particles whose dynamics follows the
Hamilton equations

$\displaystyle\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\bm{x}_{n}=\nabla_{\bm{p}_{n}}H\,,\qquad\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\bm{p}_{n}=-\nabla_{\bm{x}_{n}}H\,,$ (52)

for $n=1,\dots,N$, with hamiltonian $H=H(\Gamma_{N};t)$ and
$\Gamma_{N}=\\{\bm{x}_{1},\dots,\bm{x}_{N};\bm{p}_{1},\dots,\bm{p}_{N}\\}$.
The evolution of the corresponding phase-space density $P_{N}(\Gamma_{N};t)$
follows the Liouville equation Gaspard (2022)

$\displaystyle\partial_{t}P_{N}=\sum_{n=1}^{N}\left(\nabla_{\bm{x}_{n}}H\cdot\nabla_{\bm{p}_{n}}P_{N}-\nabla_{\bm{p}_{n}}H\cdot\nabla_{\bm{x}_{n}}P_{N}\right)=\\{H,P_{N}\\}\,,$
(53)

where we introduced the Poisson brackets
$\\{A,B\\}=\sum_{n}\nabla_{\bm{x}_{n}}A\cdot\nabla_{\bm{p}_{n}}B-\nabla_{\bm{p}_{n}}A\cdot\nabla_{\bm{x}_{n}}B$.
By extending the definition of material derivative to the phase space, _i.e._
,
$D_{t}=\partial_{t}+\sum_{n}(\nabla_{\bm{p}_{n}}H\cdot\nabla_{\bm{x}_{n}}-\nabla_{\bm{x}_{n}}H\cdot\nabla_{\bm{p}_{n}})$,
and using the Hamilton equations (52), the Liouville equation becomes

$\displaystyle
D_{t}P_{N}=\partial_{t}P_{N}+\sum_{n}\nabla_{\bm{p}_{n}}H\cdot\nabla_{\bm{x}_{n}}P_{N}-\sum_{n}\nabla_{\bm{x}_{n}}H\cdot\nabla_{\bm{p}_{n}}P_{N}=0\,.$
(54)

A coarse-grained description of the system in terms of a one-particle
probability density $P_{1}(\bm{x},\bm{p};t)$ is obtained by marginalization
over the $N-1$ positions and momenta of the other particles. Since the $N$
particles are identical and thus there are $N$ identical ways to choose one
particle, it is convenient to define the _one-particle distribution function_

$\displaystyle f_{1}(\gamma_{1};t)=N\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\Gamma_{N-1}P_{N}(\Gamma_{N};t)\,,$ (55)

which is clearly normalized to $N$. Here and in the following the volume
element is $\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\Gamma_{l}=\prod_{n=N-l+1}^{N}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\gamma_{l}$, where $\gamma_{n}=(\bm{x}_{n},\bm{p}_{n})$ and
$\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\gamma_{n}=(2\pi\hbar)^{-3}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{x}_{n}\mathop{\mathrm{\mathstrut d}}\nolimits\\!\bm{p}_{n}$
are the vector of canonical coordinates of the $n$-th particle and its
associated volume element, respectively.

The evolution equation of the one-particle distribution function is obtained
from marginalizing the Liouville equation (53). When $H$ can be written as
$H=\sum_{n}\frac{1}{2M}\bm{p}_{n}^{2}+\sum_{n}(M\phi(\bm{x}_{n})+\sum_{m>n}\psi(\left\lvert{\bm{x}_{n}-\bm{x}_{m}}\right\rvert))$,
with $\phi$ the external single particle potential per unit mass and $\psi$
the interaction potential between pairs of particles, the time evolution of
$f_{1}$ is given by the first equation of the
Bogoliubov–Born–Green–Kirkwood–Yvon (BBGKY) hierarchy Kardar (2007):

$\displaystyle\partial_{t}f_{1}(\bm{x},\bm{p})$
$\displaystyle=N\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\Gamma_{N-1}\left(M\nabla_{\bm{x}}\phi(\bm{x})\cdot\nabla_{\bm{p}}P_{N}-\bm{p}\cdot\nabla_{\bm{x}}P_{N}+\sum_{n=2}^{N}\nabla_{\bm{x}_{n}}\psi(\left\lvert{\bm{x}-\bm{x}_{n}}\right\rvert)\cdot\nabla_{\bm{p}_{n}}P_{N}\right)$
(56a) $\displaystyle=\\{H_{1},f_{1}\\}+\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\gamma_{2}\,\nabla_{\bm{x}}\psi(\left\lvert{\bm{x}-\bm{x}_{2}}\right\rvert)\cdot\nabla_{\bm{p}}f_{2}(\bm{x},\bm{p};\bm{x}_{2},\bm{p}_{2})\,.$
(56b)

In going from Eq. (56a) to Eq. (56b) we used that $P_{N}$ vanishes at the
boundary of the integration domain. Equation (56b) shows that the evolution of
$f_{1}$ splits into i) the free evolution generated by a single particle
hamiltonian $H_{1}(\gamma)=\frac{1}{2M}\bm{p}^{2}+M\phi(\bm{x})$ and ii) a
collision integral depending on the two-particles distribution function

$\displaystyle
f_{2}(\bm{x},\bm{p};\bm{x}_{2},\bm{p}_{2})=N(N-1)\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\Gamma_{N-2}P_{N}(\Gamma_{N})\,.$ (57)

This procedure can be iterated to obtain the BBGKY hierarchy of equations
specifying the evolution of the $l$-particles distribution function
$f_{l}(\gamma_{1},\dots,\gamma_{l})=\frac{N!}{(N-l)!}\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\Gamma_{N-l}\,P_{N}(\Gamma_{N})$ in terms of the $l$-particles
Hamiltonian $H_{l}(\Gamma_{l})=\sum_{n=1}^{l}(H_{1}(\gamma_{n})+\sum_{n<m\leq
l}\psi(\left\lvert{\bm{x}_{n}-\bm{x}_{m}}\right\rvert))$ and the corresponding
collision term Kardar (2007). If the hierarchy is not truncated, all the
information of the microscopic dynamics is retained and the system evolves
according to the original hamiltonian dynamics. Nevertheless, the one-particle
distribution function, the corresponding Boltzmann entropy, and the
nonequilibrium free energy satisfy exact _thermodynamic identities_ , as we
show in § III.2.

###  III.2 Thermodynamic identities

We now derive exact thermodynamic identities for the hamiltonian dynamics, and
use them to discuss the emergence of macroscopic thermodynamics by means of
coarse graining and dynamical chaoticity.

We start by introducing two definitions of entropy. The Shannon entropy

$\displaystyle S_{\text{tot}}(t)\equiv-
k_{\text{B}}\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\Gamma_{N}\,P_{N}(\Gamma_{N};t)\ln P_{N}(\Gamma_{N};t)\,$ (58)

is defined using the probability density $P_{N}(\Gamma_{N};t)$ over the whole
phase-space, while the Boltzmann entropy de Groot and Mazur (1962)

$\displaystyle S_{\text{B}}(t)\equiv-
k_{\text{B}}\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\gamma_{1}\,f_{1}(\gamma_{1};t)(\ln f_{1}(\gamma_{1};t)-1)\,$
(59)

is defined using the one-particle distribution function (55). On the one hand,
the Shannon entropy (58) is not a thermodynamic entropy: it has not the same
properties of the macroscopic entropy introduced by the local equilibrium
(17). For instance, $S_{\text{tot}}$ is a constant of motion, _i.e._ ,
$\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}S_{\text{tot}}=\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\Gamma_{N}\\{P_{N},H\\}=0$, because the Liouville equation (53)
preserves the volume in phase space. On the other hand, the Boltzmann entropy
(59) can be related to the macroscopic entropy under certain conditions that
we now examine.

To do so, we first express the Shannon entropy (58) in terms of the Boltzmann
entropy (59):

$\displaystyle S_{\text{tot}}(t)$
$\displaystyle=S_{\text{B}}(t)-k_{\text{B}}D_{KL}\left(P_{N}(t)\left|\prod_{n=1}^{N}P_{1}(t)\right.\right)+k_{\text{B}}N(\ln
N-1)\,,$ (60a)

where we used i) multiplication and division by the logarithm of the
probability distribution
$\prod_{n=1}^{N}P_{1}(t)\equiv\prod_{n=1}^{N}P_{1}(\gamma_{n};t)$ that the
system would have if there were no correlations between particles, ii) the
fact that particles are identical, and iii) the relative entropy (or
Kullback–Leibler divergence) $D_{KL}(p|q)$ between two probability
distributions $p$ and $q$

$\displaystyle D_{KL}(p|q)=\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!x\,p(x)\ln\frac{p(x)}{q(x)}\geq 0\,.$ (61)

Equation (60a), together with $\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}S_{\text{tot}}=0$, implies that

$\displaystyle
S_{\text{tot}}(t)=S_{\text{tot}}(0)=S_{\text{B}}(0)+k_{\text{B}}N(\ln N-1)\,,$
(62)

if the system is initialized in a state of completely uncorrelated particles,
_i.e._ , $P_{N}(\Gamma_{N};0)=\prod_{n=1}^{N}P_{1}(\gamma_{n};0)$. Hence, by
combining Eqs. (62) and (60a), we obtain that the change of Boltzmann entropy
is given by the relative entropy between the $N$-particle description, in
terms of the probability density $P_{N}$, and a factorized description of the
$N$-particles, by means of products of one-particle distribution $P_{1}$
(_i.e._ the mutual information with respect to the factorized state Cover
(1999)):

$\displaystyle
S_{\text{B}}(t)-S_{\text{B}}(0)=k_{\text{B}}D_{KL}\left(P_{N}\left|\prod_{n}^{N}P_{1}\right.\right)\geq
0\,.$ (63)

The thermodynamic identity (63) represents a (weaker) microscopic analog of
the second law of thermodynamics predicting the increase of Boltzmann entropy
in isolated systems (like for the macroscopic entropy in Eq. (27)). However,
Eq. (63) does not imply the monotone increase (_i.e._ , a H-theorem) of
$S_{\text{B}}$ since recurrences may occur Esposito  _et al._ (2010), unlike
Eq. (27) for the macroscopic entropy. This becomes evident when expressing the
variation $\Delta S_{\text{B}}(t)$ during a very short time increment $\Delta
t$ via Eq. (60a) together with $\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}S_{\text{tot}}=0$:

$\displaystyle\Delta t\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}S_{\text{B}}(t)$ $\displaystyle\approx
S_{\text{B}}(t+\Delta
t)-S_{\text{B}}(t)={k_{\text{B}}}D_{KL}\left(P_{N}(t+\Delta
t)\left|\prod_{n}P_{1}(t+\Delta
t)\right.\right)-{k_{\text{B}}}D_{KL}\left(P_{N}(t)\left|\prod_{n}P_{1}(t)\right.\right)\,,$
(64)

where we did not use the any hypothesis on the initial state. Indeed,
according to the thermodynamic identity in Eq. (64), a positive rate of change
of the Boltzmann entropy, namely, $\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}S_{\text{B}}(t)\geq 0$ at all times, is granted if the
contribution coming from the second relative entropy on the right-hand side of
Eq. (64) is always negligible compared to the first one, in agreement with the
intuition gained via different approaches Lanford (1975); Spohn (2012). This
can happen when the stronger requirement of molecular chaos is satisfied,
namely, that the second relative entropy on the right-hand side of Eq. (64)
vanishes at all times. We can thus conclude that the coarse-grained
description alone, in terms of the one-particle distribution function, is not
enough to guarantee irreversible thermodynamic behavior and should be
complemented with hypotheses on the dynamics, e.g., molecular chaos.

We now introduce the microscopic nonequilibrium Massieu potential Esposito
_et al._ (2010); Rao and Esposito (2018) as

$\displaystyle\Phi(t)$ $\displaystyle\equiv
S_{\text{B}}(t)-k_{\text{B}}\beta\left\llangle H_{1}\right\rrangle(t)\,,$ (65)

where the average $\left\llangle...\right\rrangle$ stands for
$\left\llangle...\right\rrangle=\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\gamma_{1}f_{1}(\gamma_{1})(...)$, in such a way that it gives
an extensive contribution, and $\beta=(k_{\text{B}}T)^{-1}$ with $T$ an
arbitrary temperature. Its equilibrium counterpart $\Phi^{\text{eq}}$ is
obtained via Eq. (65) by using the equilibrium one particle distribution
which, in the limit of impulsive interactions (_i.e._ , short-range
interaction not contributing to the total energy on average), reads

$\displaystyle
f_{1}^{\text{eq}}=NP_{1}^{\text{eq}}=\exp\\{-\beta({H_{1}}-{\overline{\mu}}_{\text{eq}})\\}\,,$
(66)

with ${\overline{\mu}}_{\text{eq}}$ being the equilibrium chemical potential
which gives the correct normalization. Note that at equilibrium
$P_{N}^{\text{eq}}(\Gamma_{N};\beta)=\prod_{n}P_{1}^{\text{eq}}(\gamma_{n};\beta)$.

In the following, we show that i) $\Phi(t)$ is upper bounded by its
equilibrium value $\Phi^{\text{eq}}$ and ii) $\Phi(t)-\Phi(0)\geq 0$ when the
initial state satisfies $P_{N}(0)=\prod_{n=0}^{N}P_{1}(\gamma_{n};0)$. First,
by using Eq. (66) in Eq. (65) and summing and subtracting $f_{1}\ln
f^{\text{eq}}$, we obtain the following thermodynamic identity for the
nonequilibrium Massieu potential:

$\displaystyle\begin{split}\Phi(t)-\Phi^{\text{eq}}&=-k_{\text{B}}\hat{D}_{KL}(f_{1}(t)||f_{1}^{\text{eq}})\leq
0\,,\end{split}$ (67)

with

$\displaystyle\hat{D}_{KL}(f|g)=\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!x\left(f(x)\ln\frac{f(x)}{g(x)}+g(x)-f(x)\right)\geq 0$ (68)

being the generalization of the relative entropy for non-normalized, positive
functions $f$ and $g$ (also called Shear Lyapunov function in the context of
chemical reaction networks Shear (1967); Higgins (1968)).

Second, by summing and subtracting $k_{\text{B}}P_{N}\sum_{n}\ln
P_{1}^{\text{eq}}$, Eq. (63) can be rewritten as 111 Note that
$D_{KL}(P_{N}|P_{N}^{\text{eq}})\geq\hat{D}_{KL}(f_{1}|f_{1}^{\text{eq}})$
always holds (i.e., independently of the factorization of the initial state)
and implies that the one-particle distribution function $f_{1}$ is always
closer to the corresponding equilibrium than the $N$-particles density $P_{N}$

$\displaystyle S_{\text{B}}(t)-S_{\text{B}}(0)$
$\displaystyle=k_{\text{B}}D_{KL}\left(P_{N}\left|P_{N}^{\text{eq}}\right.\right)-k_{\text{B}}\hat{D}_{KL}\left(f_{1}\left|f_{1}^{\text{eq}}\right.\right)\geq
0\,,$ (69a)

which, together with Eq. (67) and the conservation of the Shannon entropy (58)
implying that $D_{KL}(P_{N}|P_{N}^{\text{eq}})$ is conserved too, leads to the
conservation of the single-particle energy:

$\displaystyle
0=D_{KL}\left(P_{N}(t)\left|P_{N}^{\text{eq}}\right.\right)-D_{KL}\left(P_{N}(0)\left|P_{N}^{\text{eq}}\right.\right)$
$\displaystyle=\beta(\left\llangle H_{1}\right\rrangle{(t)}-\left\llangle
H_{1}\right\rrangle{(0)})\,.$ (70a)

Hence,

$\displaystyle\Phi(t)-\Phi(0)=S_{\text{B}}(t)-S_{\text{B}}(0)\geq 0\,,$ (71)

when the system is initialized in a state of completely uncorrelated
particles, _i.e._ , $P_{N}(\Gamma_{N};0)=\prod_{n=1}^{N}P_{1}(\gamma_{n};0)$
and Eq. (63) holds.

We notice that similar thermodynamic identities have been already observed in
different contexts. A thermodynamic identity similar to Eq. (63) has been
already obtained using integral fluctuation relations Kawai  _et al._ (2007);
Gaspard (2022). The difference is that we compare the $N$-particles
distribution with a (completely decorrelated) product of one-particle
distributions, while they introduce the backward hamiltonian evolution which
enters the integral fluctuation relations. The thermodynamic identities (63),
(64), (67) are also similar to results obtained for systems coupled to
equilibrium baths Esposito  _et al._ (2010), but here there is no distinction
between degrees of freedom of the system and of the baths. The identities
presented in this section rely on coarse-graining the many degrees of freedom
in favor of a description in terms of the one particle distribution and the
factorization of the initial distribution. They are exact results which are
consistent with the thermodynamic intuition, but do not lead to a full
irreversible thermodynamics with a positive EPR. In fact, the microscopic
dynamics may not be chaotic enough to obtain thermodynamic behavior. A fully
consistent thermodynamic description can be obtained at the mesoscopic level
by using the assumption of molecular chaos as we now proceed to do in Section
IV.

##  IV Mesoscopic Scale: Boltzmann equation

The standard approach to introduce irreversibility in the evolution of the
one-particle distribution function relies of on the crucial assumption of
_molecular chaos_ Kardar (2007), which means that particles are assumed to be
uncorrelated before they collide Spohn (2012). The resulting mesoscopic
evolution is the Boltzmann equation. We now review the dynamics and
thermodynamics of the Boltzmann equation both in general and after applying
the Chapman–Enskog expansion to obtain the hydrodynamic equations.

###  IV.1 Dynamics

The Boltzmann equation

$\displaystyle\partial_{t}f_{1}+\\{f_{1},H_{1}\\}=C(f_{1},f_{1})\,$ (72)

is a closed dynamical equation for the one-particle distribution function
$f_{1}$ resulting from truncating the BBGKY hierarchy de Groot and Mazur
(1962); Kardar (2007); Gaspard (2022) at the first equation (56b) and
approximating the collision integral involving the two-particles distribution
function $f_{2}$ in terms of the product of two one-particle distribution
functions times a scattering rate $W$:

$\displaystyle C(f_{1},f_{1})=\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}_{2}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}_{1}^{\prime}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}_{2}^{\prime}\,\left(W(\bm{p},\bm{p}_{2}|\bm{p}_{1}^{\prime},\bm{p}_{2}^{\prime})f_{1}(\bm{x},\bm{p}_{1}^{\prime})f_{1}(\bm{x},\bm{p}_{2}^{\prime})-W(\bm{p}_{1}^{\prime},\bm{p}_{2}^{\prime}|\bm{p},\bm{p}_{2})f_{1}(\bm{x},\bm{p})f_{1}(\bm{x},\bm{p}_{2})\right)\,,$
(73)

Because of the underlying Hamiltonian dynamics Van Kampen (1992), the
scattering rate $W$ must satisfy the time-reversal symmetry

$\displaystyle
W(\bm{p}_{1},\bm{p}_{2}|\bm{p}_{1}^{\prime},\bm{p}_{2}^{\prime})=W(-\bm{p}_{1}^{\prime},-\bm{p}_{2}^{\prime}|-\bm{p}_{1},-\bm{p}_{2})\,,$
(74)

and be invariant under a generic rotation which implies invariance under
inversion of the coordinate reference frame
$(\bm{x},\bm{p})\mapsto(-\bm{x},-\bm{p})$

$\displaystyle
W(\bm{p}_{1},\bm{p}_{2}|\bm{p}_{1}^{\prime},\bm{p}_{2}^{\prime})=W(-\bm{p}_{1},-\bm{p}_{2}|-\bm{p}_{1}^{\prime},-\bm{p}_{2}^{\prime})\,.$
(75)

Combining Eqs. (74) and (75) leads to the microscopic reversibility condition:

$\displaystyle
W(\bm{p}_{1},\bm{p}_{2}|\bm{p}_{1}^{\prime},\bm{p}_{2}^{\prime})=W(\bm{p}_{1}^{\prime},\bm{p}_{2}^{\prime}|\bm{p}_{1},\bm{p}_{2})\,.$
(76)

In this framework, the so-called summational invariant functions
$A(\bm{x},\bm{p})$ represents quantities satisfing the conservation law

$A(\bm{x},\bm{p}_{1})+A(\bm{x},\bm{p}_{2})-A(\bm{x},\bm{p}_{1}^{\prime})-A(\bm{x},\bm{p}_{2}^{\prime})=0\,,$
(77)

when two particles collide in $\bm{x}$ and $\bm{p}_{1}$, $\bm{p}_{2}$,
$\bm{p}_{1}^{\prime}$, and $\bm{p}_{2}^{\prime}$ satisfy the equation of
motion, _i.e._ ,
$W(\bm{p}_{1}^{\prime},\bm{p}_{2}^{\prime}|\bm{p}_{1},\bm{p}_{2})\neq 0$.
Their integral over the momentum is zero when multiplied by the collision
integral de Groot and Mazur (1962):

$\displaystyle\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}\,A(\bm{x},\bm{p})C(f_{1},f_{1})=0\,.$ (78)

This can be proved by multiplying the conservation law (77) by the integrand
of $C(f_{1},f_{1})$, integrating over all the momenta and using the symmetry
(76) and indistinguishability of particles to obtain four times the integral
in Eq. (78). Functions of the position only $A(\bm{x})$ and constant functions
are always summational invariants. Moreover, since the scattering rate $W$
derives from an underlying hamiltonian dynamics, also the kinetic energy and
the momentum are summational invariants:
$(\bm{p}_{1})^{2}+(\bm{p}_{2})^{2}-(\bm{p}_{1}^{\prime})^{2}-(\bm{p}_{2}^{\prime})^{2}=0$
and $\bm{p}_{1}+\bm{p}_{2}-\bm{p}_{1}^{\prime}-\bm{p}_{2}^{\prime}=0$. Note
that Eq. (78) implies that $\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}\,C(f_{1},f_{1})=0$ and, therefore, that the Boltzmann
equation (72) conserves the total number of particles.

To describe the behavior of the system in the configuration space, we
introduce the (mesoscopic) spatio-temporal fields mass density and velocity as

$\displaystyle{\left\langle\rho\right\rangle}(\bm{x})$
$\displaystyle=M\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\,f_{1}(\bm{x},\bm{p})\,,$ (79a)
$\displaystyle{\left\langle\bm{v}\right\rangle}(\bm{x})$
$\displaystyle=({\left\langle\rho\right\rangle}(\bm{x}))^{-1}\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\,\bm{p}\,f_{1}(\bm{x},\bm{p})\,.$
(79b)

Their balance equations are obtained by directly using the Boltzmann equation
(72) in the time derivative of Eqs. (79), and the material derivative
$D_{t}=\partial_{t}+{\left\langle\bm{v}\right\rangle}\cdot\nabla_{\bm{x}}$.
The balance equation for mass density reads

$\displaystyle
D_{t}{\left\langle\rho\right\rangle}=-{\left\langle\rho\right\rangle}\nabla_{\bm{x}}\cdot{\left\langle\bm{v}\right\rangle}\,,$
(80)

while the balance equation for the average velocity reads

$\displaystyle{\left\langle\rho\right\rangle}D_{t}{\left\langle\bm{v}\right\rangle}=-\nabla_{\bm{x}}\cdot{\left\langle
P\right\rangle}+{\left\langle\rho\right\rangle}\bm{F}\,,$ (81)

where the pressure tensor is defined as

$\displaystyle{\left\langle
P\right\rangle}(\bm{x})=\frac{1}{M}\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\,(\bm{p}-M{\left\langle\bm{v}\right\rangle}(\bm{x}))(\bm{p}-M{\left\langle\bm{v}\right\rangle}(\bm{x}))f_{1}(\bm{x},\bm{p})\,,$
(82)

with
$(\bm{p}-M{\left\langle\bm{v}\right\rangle})(\bm{p}-M{\left\langle\bm{v}\right\rangle})$
being the matrix whose entries are $({p_{i}}-M{{\left\langle
v_{i}\right\rangle}})({p_{j}}-M{{\left\langle v_{j}\right\rangle}})$ ($i$ and
$j$ labeling here the three spatial components) and
$\bm{F}=-\nabla_{\bm{x}}\phi$. Note that these mesoscopically-derived balance
equations (80) and (81) correspond exactly to the macroscopic Navier-Stokes
equations (1).

###  IV.2 Thermodynamics

We now introduce the mesoscopic definitions, as well as their respective
balance equations, of the specific thermodynamic quantities examined at the
macroscopic level in Section II. The (specific) internal energy is given by

$\displaystyle{\left\langle u\right\rangle}(\bm{x})$
$\displaystyle=\frac{({\left\langle\rho\right\rangle}(\bm{x}))^{-1}}{2M}\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}(\bm{p}-M{\left\langle\bm{v}\right\rangle}(\bm{x}))^{2}f_{1}(\bm{x},\bm{p})\,,$
(83)

and its balance equation, obtained using again the Boltzmann equation (72),
reads

$\displaystyle{\left\langle\rho\right\rangle}D_{t}{\left\langle
u\right\rangle}$
$\displaystyle=-\nabla_{\bm{x}}\cdot{\left\langle\bm{J}_{q}\right\rangle}-\operatorname{Tr}{\left\langle
P\right\rangle}\nabla_{\bm{x}}{\left\langle\bm{v}\right\rangle}\,,$ (84)

where the heat flux is given by

$\displaystyle{\left\langle\bm{J}_{q}\right\rangle}(\bm{x})=\frac{1}{2M^{2}}\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\,(\bm{p}-M{\left\langle\bm{v}\right\rangle}(\bm{x}))^{2}(\bm{p}-M{\left\langle\bm{v}\right\rangle}(\bm{x}))f_{1}(\bm{x},\bm{p})\,.$
(85)

The specific entropy per unit mass defined, for consistency with Eq. (59), as

$\displaystyle{\left\langle s_{\text{B}}\right\rangle}(\bm{x})\equiv-
k_{\text{B}}({\left\langle\rho\right\rangle}(\bm{x}))^{-1}\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\,f_{1}(\bm{x},\bm{p})(\ln
f_{1}(\bm{x},\bm{p})-1)\,,$ (86)

satisfies the balance equation

$\displaystyle{\left\langle\rho\right\rangle}D_{t}{\left\langle
s_{\text{B}}\right\rangle}=-\nabla_{\bm{x}}\cdot{\left\langle\bm{J}_{s}\right\rangle}+{\left\langle\dot{\sigma}_{\text{B}}\right\rangle}\,.$
(87)

Here, the entropy flux is defined as

$\displaystyle{\left\langle\bm{J}_{s}\right\rangle}(\bm{x})=-k_{\text{B}}\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\,\left(\frac{\bm{p}}{M}-{\left\langle\bm{v}\right\rangle}(\bm{x})\right)f_{1}(\bm{x},\bm{p})\ln
f_{1}(\bm{x},\bm{p})\;.$ (88)

It is worth noting that, unlike the corresponding macroscopic fluxes in Eq.
(24), the mesoscopic entropy flux (88) is not proportional to the mesoscopic
heat flux (85). This difference emerges from the absence, in general, of the
local equilibrium condition. We will verify in § IV.4 that if $f_{1}$
satisfies the local equilibrium, the proportionality between
${\left\langle\bm{J}_{s}\right\rangle}$ and
${\left\langle\bm{J}_{q}\right\rangle}$ is recovered. Omitting the spatial
dependence for shortness and with the notation $\mathop{\mathrm{\mathstrut
d}}\nolimits^{4}\\!p=(2\pi\hbar)^{-3}\prod_{i=1,2}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}_{i}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}_{i}^{\prime}$, the entropy production rate
${\left\langle\dot{\sigma}_{\text{B}}\right\rangle}$ reads

$\displaystyle{\left\langle\dot{\sigma}_{\text{B}}\right\rangle}(\bm{x})=\frac{k_{\text{B}}}{4}\int\mathop{\mathrm{\mathstrut
d}}\nolimits^{4}\\!p\,W(\bm{p}_{1},\bm{p}_{2}|\bm{p}_{1}^{\prime},\bm{p}_{2}^{\prime})(f_{1}(\bm{p}_{1}^{\prime})f_{1}(\bm{p}_{2}^{\prime})-f_{1}(\bm{p}_{1})f_{1}(\bm{p}_{2}))\ln\frac{W(\bm{p}_{1},\bm{p}_{2}|\bm{p}_{1}^{\prime},\bm{p}_{2}^{\prime})f_{1}(\bm{p}_{1}^{\prime})f_{1}(\bm{p}_{2}^{\prime})}{W(\bm{p}_{1}^{\prime},\bm{p}_{2}^{\prime}|\bm{p}_{1},\bm{p}_{2})f_{1}(\bm{p}_{1})f_{1}(\bm{p}_{2})}\geq
0\,,$ (89)

where the inequality follows from microscopic reversibility (76). Note that
the balance equation (87) implies that for an isolated system

$\displaystyle\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}{{S_{\text{B}}}}=\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{x}{\left\langle\rho\right\rangle}(\bm{x}){\left\langle
s_{\text{B}}\right\rangle}(\bm{x})=\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{x}\,{\left\langle\dot{\sigma}_{\text{B}}\right\rangle}(\bm{x})\geq
0\,,$ (90)

namely, ${{S_{\text{B}}}}$ defined in Eq. (59) monotonously increases in time
according to the Boltzmann equation, in contrast to what we obtained in the
microscopic theory (see Eq. (63)). This is an alternative statement of
Boltzmann’s H-theorem de Groot and Mazur (1962).

The theory developed here is the mesoscopic analogue of the macroscopic
thermodynamics discussed in Section II, but they are still not equivalent: the
crucial local equilibrium condition (20), or equivalently, the Euler relation
(21), has not been derived yet. To recover them, we need to introduce
perturbative solutions of the Boltzmann equation as we now do in the next
subsection IV.3.

###  IV.3 Linking the Mesoscopic and Macroscopic Description via the
Chapman–Enskog Expansion

Many perturbative solutions of the Boltzmann equation have been derived de
Groot and Mazur (1962); Cercignani (1988). Since we are interested in the
dynamics of the spatio-temporal fields accounted by the balance equations
(80), (81), (84), (87), we apply the Chapman–Enskog scheme, _i.e._ , a type of
multiple-scale perturbation theory Van Kampen (1987).

####  IV.3.1 Zeroth-order of the Chapman–Enskog expansion.

We start by identifying two (time and spatial) scales of the Boltzmann
dynamics. We introduce a macroscopic length scale $L$ on which the external
potential $\phi$ varies, and the associated macroscopic timescale
$\tau=L/v_{m}$, with $v_{m}$ the typical _molecular_ velocity. The second
timescale corresponds to the mean free time between collisions $\tau_{f}$ and
is related to the mean free path $l=v_{m}\tau_{f}$. The Knudsen number is
defined as the ratio between these two scales,

$\displaystyle\varepsilon=l/L=\tau_{f}/\tau\,,$ (91)

and plays the role of the small parameter in the Chapman–Enskog (perturbative)
scheme. Indeed, by using these length and time scales to define adimensional
variables (marked by the primes),

$\displaystyle t=\tau
t^{\prime}\,,\qquad\bm{x}=L\bm{x}^{\prime}\,,\qquad\bm{p}=Mv_{m}\bm{p}^{\prime}\,,\qquad\bm{F}=\frac{L}{\tau^{2}}\bm{F}^{\prime}\,,\qquad
C=\frac{C^{\prime}}{\tau_{f}^{-1}}\,,\qquad$ (92)

the Boltzmann equation (72) becomes (removing the primes for compactness)

$\displaystyle\partial_{t}f_{1}+\,\bm{p}\cdot\nabla_{\bm{x}}f_{1}+\bm{F}\cdot\nabla_{\bm{p}}f_{1}=\varepsilon^{-1}C(f_{1},f_{1})\,.$
(93)

We are interested in the limit $\varepsilon\to 0$, in which the collision term
dominates the evolution of the one-particle distribution.

We thus consider the expansion

$\displaystyle f_{1}=f^{(0)}+\varepsilon
f^{(1)}+\varepsilon^{2}f^{(2)}\dots\,,$ (94)

and substitute it in Eq. (93). The zeroth-order solution is given by any one-
particle distribution function for which the collision integral of the
Boltzmann equation vanishes:

$\displaystyle C(f^{(0)},f^{(0)})$ $\displaystyle=0\,.$ (95)

To identify $f^{(0)}$, we notice that the one-particle distribution function
such that entropy production rate in Eq. (89) vanishes must also satisfy Eq.
(95). Thus, we look for solutions of
${\left\langle\dot{\sigma}_{\text{B}}\right\rangle}=0$ which can be rewritten
in the form (78) for the quantity $\ln f^{(0)}$ by using the symmetry in Eq.
(76) and the indistinguishability of the particles. This implies that $\ln
f^{(0)}$, as well as $f^{(0)}$ because of the monotonicity of the logarithm,
must be a function of summational invariants, like mass, momentum and kinetic
energy. We can thus verify that the local Maxwell–Boltzmann distribution

$f^{\text{MB}}(\bm{x},\bm{p})=\exp\left\\{-\frac{1}{T(\bm{x})}\left(\frac{\left(\bm{p}-{\bm{v}}(\bm{x})\right)^{2}}{2}-\mu(\bm{x})\right)\right\\}\,,$
(96)

with, in general, $\bm{x}$-dependent and $t$-dependent temperature
$T(\bm{x})$, velocity ${\bm{v}}(\bm{x})$, chemical potential of the ideal gas

$\mu(\bm{x})=T(\bm{x})\ln\left(z(\bm{x})\right)-\frac{3}{2}T(\bm{x})\ln\frac{T(\bm{x})}{2\pi\hbar^{2}}\,,$
(97)

and density $z(\bm{x})$ (satisfying $\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{x}z(\bm{x})=N$), solves Eq. (95), namely,
$f^{(0)}=f^{\text{MB}}$. Here all the quantity in Eqs. (96) and (97) are
written using the scaling defined in Eq. (92) (in particular, energy,
temperature and $\hbar$ are expressed in ${Mv_{m}^{2}}$,
${Mv_{m}^{2}}/{k_{\text{B}}}$, and $MLv_{m}$ units, respectively). At this
stage, the temperature $T(\bm{x})$, the velocity ${\bm{v}}(\bm{x})$ and the
density $z(\bm{x})$ in Eq. (96) are arbitrary parameters of the
Maxwell–Boltzmann distribution, but they will be determined in the following
using the balance equations. For convenience, we recast the expansion in Eq.
(94) truncated at the first order using the Maxwell–Boltzmann distribution as

$\displaystyle
f_{1}(\bm{x},\bm{p})=f^{\text{MB}}(\bm{x},\bm{p})(1+\varepsilon\chi(\bm{x},\bm{p}))+O(\varepsilon^{2})\,,$
(98)

namely, defining $\chi(\bm{x},\bm{p})\equiv
f^{(1)}(\bm{x},\bm{p})/f^{\text{MB}}(\bm{x},\bm{p})$.

Now we examine how the parameters of the Maxwell–Boltzmann distribution (96),
that is the zeroth-order solution of (93), are related to the the (mesoscopic)
spatio-temporal fields mass density ${\left\langle\rho\right\rangle}(\bm{x})$
and velocity ${\left\langle\bm{v}\right\rangle}(\bm{x})$ introduced in § IV.1.
By using Eq. (96) in Eqs. (79a) and (79b), we obtain

$\displaystyle{\left\langle\rho\right\rangle}(\bm{x})$
$\displaystyle=\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\,f^{\text{MB}}(\bm{x},\bm{p},0)=z(\bm{x})\,,$
(99a)
$\displaystyle{\left\langle\rho\right\rangle}(\bm{x}){\left\langle\bm{v}\right\rangle}(\bm{x})$
$\displaystyle=\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\,\bm{p}f^{\text{MB}}(\bm{x},\bm{p},0)=z(\bm{x}){\bm{v}}(\bm{x})\,,$
(99b)

and, consequently, the density $z(\bm{x})$ and velocity ${\bm{v}}(\bm{x})$ are
determined by the balance equations (80) and (81) which can be rewritten now
as

$\displaystyle D_{t}z$ $\displaystyle=-z\nabla\cdot{\bm{v}}\,,$ (100a)
$\displaystyle zD_{t}{\bm{v}}$ $\displaystyle=-\nabla p+z\bm{F}\,,$ (100b)

by solving the integral in Eq. (82) using Eq. (96) which leads to a diagonal
pressure tensor,

${\left\langle
P_{ij}\right\rangle}(\bm{x})=p(\bm{x})\delta_{ij}=z(\bm{x})T(\bm{x})\delta_{ij}\,.$
(101)

On the one hand, the zeroth-order Chapman–Enskog expansion recovers the
macroscopic mass balance equation (namely, Eq. (100a) coincides with Eq.
(1a)). On the other hand, it describes only fluids without viscous friction
since a diagonal pressure tensor enters the zeroth order momentum balance
equation (100b) unlike in Eq. (1b).

Second, we focus on the thermodynamic quantities introduced in § IV.2, as well
as their balance equations. By using Maxwell–Boltzmann distribution (96), the
(specific) internal energy (83) becomes

${\left\langle\rho\right\rangle}(\bm{x}){\left\langle
u\right\rangle}(\bm{x})=\frac{1}{2}\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\,(\bm{p}-{{\bm{v}}})^{2}f^{\text{MB}}(\bm{x},\bm{p},0)=\frac{3}{2}z(\bm{x})T(\bm{x})\,,$
(102)

while its balance equation (84) reads

$\displaystyle D_{t}T=-\frac{2}{3}T\nabla_{\bm{x}}\cdot{\bm{v}}\,,$ (103)

since the heat flux vanishes, _i.e._ ,
${\left\langle\bm{J}_{q}\right\rangle}(\bm{x})=0$. The balance equation (103)
implies that the temperature $T(\bm{x})$ changes only due to the transport due
to the velocity field. Analogously, the balance equation (87) for the local
entropy becomes

$\displaystyle D_{t}{\left\langle s_{\text{B}}\right\rangle}=0\,,$ (104)

since the Maxwell–Boltzmann distribution leads to a vanishing entropy flux
${\left\langle\bm{J}_{s}\right\rangle}(\bm{x})=0$ and entropy producation rate
${\left\langle\dot{\sigma}_{\text{B}}\right\rangle}=0$. This physically means
that, at the zeroth-order of Chapman–Enskog expansion, the evolution of a
fluid is isoentropic along the flow lines.

Finally, by using the using the Maxwell–Boltzmann distribution (96) in Eq.
(86) together with Eq. (102) and $p(\bm{x})=z(\bm{x})T(\bm{x})$, we can write
the specific entropy as

$T(\bm{x}){\left\langle s_{\text{B}}\right\rangle}(\bm{x})={\left\langle
u\right\rangle}(\bm{x})+\frac{p(\bm{x})}{z(\bm{x})}-\mu(\bm{x})\,,$ (105)

which represents the mesoscopic formulation of the Euler relation (21). Note
that, even if the zeroth-order of Chapman–Enskog expansion is consistent with
the Euler relation (21), or equivalently with the local equilibrium condition
(20), this consistency is a trivial one. In fact, the dissipative effects
characterizing the macroscopic thermodynamics and responsible for the
variation of the local entropy are not recovered (see Eq. (104)). We thus move
to the next order of the Chapman–Enskog expansion.

####  IV.3.2 First-order of the Chapman–Enskog expansion.

The first-order equation of the Chapman–Enskog expansion of the (adimensional)
Boltzmann equation (93) can be written as

$\displaystyle\begin{split}{D_{t}}f^{\text{MB}}&+(\bm{p}-{\bm{v}})\cdot\nabla_{\bm{x}}f^{\text{MB}}+\bm{F}\cdot\nabla_{\bm{p}}f^{\text{MB}}=\\\
&=\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}_{2}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}_{1}^{\prime}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}_{2}^{\prime}\,f^{\text{MB}}(\bm{x},\bm{p})f^{\text{MB}}(\bm{x},\bm{p}_{2})W(\bm{p}_{1}^{\prime},\bm{p}_{2}^{\prime}|\bm{p},\bm{p}_{2})\left(\chi(\bm{p}_{1}^{\prime})+\chi(\bm{p}_{2}^{\prime})-\chi(\bm{p})-\chi(\bm{p}_{2})\right)\,,\end{split}$
(106)

by using Eq. (98), the material derivative
$D_{t}\equiv\partial_{t}+{\bm{v}}\cdot\nabla_{\bm{x}}$, and the conservation
of the kinetic energy and momentum (implying,
$(\bm{p}_{1}-{\bm{v}})^{2}+(\bm{p}_{2}-{\bm{v}})^{2}=(\bm{p}_{1}^{\prime}-{\bm{v}})^{2}+(\bm{p}_{2}{{}^{\prime}}-{\bm{v}})^{2}$)
to collect the Maxwell–Boltzmann distributions in the collision integral. The
first order correction $\chi$ must not contribute to the value of the spatio-
temporal fields mass ${\left\langle\rho\right\rangle}(\bm{x})$, momentum
${\left\langle\bm{v}\right\rangle}(\bm{x})$, and internal energy
${\left\langle u\right\rangle}(\bm{x})$ since they are already recovered from
the the parameters $z(\bm{x})$, ${\bm{v}}(\bm{x})$, and $T(\bm{x})$ of the
zeroth-order solution (96) according to Eqs. (99a), (99b), and (102),
respectively. This implies that

$\displaystyle\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\,\chi(\bm{x},\bm{p})f^{\text{MB}}(\bm{x},\bm{p})=0\,,$
(107a) $\displaystyle\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\,\bm{p}\,\chi(\bm{x},\bm{p})f^{\text{MB}}(\bm{x},\bm{p})=0\,,$
(107b) $\displaystyle\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\,(\bm{p}-{{\bm{v}}})^{2}\chi(\bm{x},\bm{p})f^{\text{MB}}(\bm{x},\bm{p})=0\,.$
(107c)

However, the correction $\chi$ should affect the balance equations of
momentum, internal energy and entropy via the non-diagonal elements of the
pressure tensor and the heat flux and, furthermore, lead to a nonvanishing
entropy flux and entropy production rate.

To determine $\chi$, we start by expressing the terms on left-hand side of Eq.
(106) using Eqs. (100a), (100b), and (103). We obtain

$\displaystyle D_{t}\ln f^{\text{MB}}$
$\displaystyle=\left(-\frac{2}{3}\nabla_{\bm{x}}\cdot{\bm{v}}\right)\frac{(\bm{p}-{{\bm{v}}})^{2}}{{2}T}+\frac{\bm{p}-{{\bm{v}}}}{T}\cdot\left(-\nabla_{\bm{x}}T-\frac{T}{z}\nabla_{\bm{x}}z+\bm{F}\right)\,,$
(108a) $\displaystyle\nabla_{\bm{x}}\ln f^{\text{MB}}$
$\displaystyle=\nabla_{\bm{x}}T\left(\frac{(\bm{p}-{{\bm{v}}})^{2}}{2T^{2}}-\frac{3}{2T}\right)+\frac{\sum_{i,j}({p}_{i}-{{{v}}_{j}})\partial_{x_{j}}{{v}}_{i}}{T}+\frac{1}{z}\nabla_{\bm{x}}z\,,$
(108b) $\displaystyle\nabla_{\bm{p}}\ln f^{\text{MB}}$
$\displaystyle=-\left(\bm{p}-{{\bm{v}}}\right)\,,$ (108c)

with $i$ and $j$ labeling the three spatial components. By using the strain
tensor
$\epsilon_{ij}=\frac{1}{2}(\partial_{x_{i}}{{{v}}}_{j}+\partial_{x_{j}}{{{v}}}_{i})$,
the left-hand side of Eq. (106) becomes

$\displaystyle\begin{split}&D_{t}f^{\text{MB}}+(\bm{p}-{{\bm{v}}})\cdot\nabla_{\bm{x}}f^{\text{MB}}+\bm{F}\cdot\nabla_{\bm{p}}f^{\text{MB}}=\\\
&=f^{\text{MB}}\left\\{\frac{1}{T}\sum_{ij}\left((p_{i}-{{{v}}}_{i})(p_{j}-{{{v}}}_{j})-\frac{{1}}{3}(\bm{p}-{{\bm{v}}})^{2}\delta_{ij}\right)\epsilon_{ij}+\frac{1}{T}\left(\frac{(\bm{p}-{{\bm{v}}})^{2}}{2T}-\frac{5}{2}\right)(\bm{p}-{{\bm{v}}})\cdot\nabla_{\bm{x}}T\right\\}\,.\end{split}$
(109)

Comparing this expression with Eq. (106) suggests that $\chi$ is a linear
function of the temperature $\nabla_{\bm{x}}\bm{T}$ and velocity
$\\{\partial_{x_{j}}{{{v}}}_{i}\\}$ gradients de Groot and Mazur (1962). To
determine its expression, one could introduce an expansion in appropriate
orthogonal polynomials Chapman and Cowling (1990). However, for the sake of
simplicity, we will employ the _relaxation time approximation_ or _single
collision time approximation_ Kardar (2007); Lundstrom (2002).

This approximation corresponds to use the properties of the Maxwell–Boltzmann
distribution $f^{\text{MB}}$ to obtain the following estimate of the collision
operator on the right-hand side of Eq. (106),

$\displaystyle\small\begin{split}C(f,f)&\simeq\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}_{2}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}_{1}^{\prime}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}_{2}^{\prime}\,W(\bm{p},\bm{p}_{2}|\bm{p}_{1}^{\prime},\bm{p}_{2}^{\prime})f^{\text{MB}}(\bm{p}_{1}^{\prime})f^{\text{MB}}(\bm{p}_{2}^{\prime})\left(\chi(\bm{p}_{1}^{\prime})+\chi(\bm{p}_{2}^{\prime})\right)+\\\
&\quad-f^{\text{MB}}(\bm{p})\chi(\bm{p})\underbrace{\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}_{2}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}_{1}^{\prime}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}_{2}^{\prime}\,f^{\text{MB}}(\bm{p}_{2})W(\bm{p},\bm{p}_{2}|\bm{p}_{1}^{\prime},\bm{p}_{2}^{\prime})\left(1+\frac{\chi(\bm{p}_{2})}{\chi(\bm{p})}\right)}_{1/\tau_{R}}\approx-\tau_{R}^{-1}f^{\text{MB}}(\bm{p})\chi(\bm{p})\,.\end{split}$
(110)

Here, we assumed that before a collision takes place the distribution is the
equilibrium one, _i.e._ ,
$\chi(\bm{p}^{\prime}_{1})=\chi(\bm{p}^{\prime}_{2})\approx 0$. Furthermore,
we assumed that the factor
$f^{\text{MB}}(\bm{p}_{2})W(\bm{p},\bm{p}_{2}|\bm{p}_{1}^{\prime},\bm{p}_{2}^{\prime})$
decays fast enough as a function of the difference $\bm{p}_{2}-\bm{p}$ that
only the zeroth-order contribution in the expansion around $\bm{p}$, namely
$\chi(\bm{p}_{2})\simeq\chi(\bm{p})$, can be retained. Therefore, the
resulting scattering timescale $\tau_{R}$ does not depend on $\chi$ or on the
gradients of temperature or velocity.

Using the relaxation time approximation (110) together with Eq. (109) in Eq.
(106), we conclude that

$\displaystyle\chi(\bm{p})$
$\displaystyle=-\frac{\tau_{R}}{T}\left\\{\sum_{ij}\left((p_{i}-{{{v}}}_{i})(p_{j}-{{{v}}}_{j})-\frac{{1}}{3}(\bm{p}-{{\bm{v}}})^{2}\delta_{ij}\right)\epsilon_{ij}+\left(\frac{(\bm{p}-{{\bm{v}}})^{2}}{2T}-\frac{5}{2}\right)(\bm{p}-{{\bm{v}}})\cdot\nabla_{\bm{x}}T\right\\}\,.$
(111)

The relaxation time approximation is only appropriate for modeling near-
equilibrium systems Lundstrom (2002). A more rigorous treatment leads to the
emergence of two different scattering timescales associated to viscous forces,
accounted by the first term on the right-hand side of Eq. (111), and thermal
diffusion, accounted by the second term on the right-hand side of Eq. (111).
This discrepancy becomes important in determining the quantitative values of
the pressure tensor (82) and the heat flux (85), but it is not fundamental for
our qualitative discussion of thermodynamics, where we will further simplify
the expressions assuming $\tau_{R}=1$ in the following.

We notice that $\chi$ given in Eq. (111) satisfies the constraints in Eqs.
(107a), (107b), and (107c) as it can be verified using the properties of the
Gaussian distribution:

$\left\langle(p_{i}-{{v}}_{i})(p_{j}-{{v}}_{j})\right\rangle_{\text{MB}}=z(\bm{x})T(\bm{x})\delta_{ij}\,,$
(112a)
$\left\langle(p_{i}-{{v}}_{i})(p_{j}-{{v}}_{j})(p_{k}-{{v}}_{k})(p_{l}-{{v}}_{l})\right\rangle_{\text{MB}}=z(\bm{x})(T(\bm{x}))^{2}(\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})\,,$
(112b)

We have now all the elements to assess the consistency of the first order
solution (98) (with $\chi$ given in Eq. (111)) with the macroscopic theory, as
originally done by Prigogine Prigogine (1949).

First, we consider the balance equations for the (mesoscopic) spatio-temporal
fields mass ${\left\langle\rho\right\rangle}(\bm{x})$, momentum
${\left\langle\bm{v}\right\rangle}(\bm{x})$, and internal energy
${\left\langle u\right\rangle}(\bm{x})$ given in Eqs. (80), (81), and (84).
These equations have the same form of the macroscopic balance equations (1a),
(1b) and (14) using the zeroth-order solution. Therefore, we only need to use
the correction $\chi$ to evaluate the changes in the pressure tensor
${\left\langle P\right\rangle}$ in Eq. (82) and the heat flux
${\left\langle\bm{J}_{q}\right\rangle}$ in Eq. (85) which become equal to

$\displaystyle{\left\langle P_{ij}\right\rangle}(\bm{x})$
$\displaystyle=p(\bm{x})\delta_{ij}+{\left\langle\Pi_{ij}\right\rangle}(\bm{x})=\left(p+\frac{2}{3}\varepsilon
zT\sum_{k}\epsilon_{kk}\right)\delta_{ij}-2\varepsilon zT\epsilon_{ij}\,,$
(113a) $\displaystyle{\left\langle\bm{J}_{q}\right\rangle}(\bm{x})$
$\displaystyle=-\frac{5}{2}\varepsilon zT\nabla_{\bm{x}}T\,.$ (113b)

Compared to the zeroth-order expression (101) we obtain the additional viscous
part of the tensor pressure as well as a nonzero thermal conductivity.
Equations. (113a) and (113b) show a linear relation between
${\left\langle\Pi\right\rangle}$ and ${\left\langle\bm{J}_{q}\right\rangle}$,
on the one hand, and their conjugate thermodynamic forces
$\\{\epsilon_{ij}\\}$ and $\nabla_{\bm{x}}T$, on the other hand. This means
that the first-order of the Chapman–Enskog expansion describes fluids in the
linear regime, namely, close to equilibrium.

Second, the local thermodynamic equilibrium is left unaffected by the
inclusion of the first Chapman–Enskog correction: substituting the the first
order solution (98) (with $\chi$ in Eq. (111)) in Eq. (86), we obtain again
Eq. (105), as a consequence of the constraints in Eqs. (107a), (107b), and
(107c). On the other hand, the first order solution leads to a nonvanishing
entropy production rate which reads

$\displaystyle{\left\langle\dot{\sigma}_{\text{B}}\right\rangle}=\frac{\varepsilon}{4}\int\mathop{\mathrm{\mathstrut
d}}\nolimits^{4}\\!p\,W(\bm{p}_{1}^{\prime},\bm{p}_{2}^{\prime}|\bm{p}_{1},\bm{p}_{2})f^{\text{MB}}(\bm{p}_{1})f^{\text{MB}}(\bm{p}_{2})\left(\chi(\bm{p}_{1}^{\prime})+\chi(\bm{p}_{2}^{\prime})-\chi(\bm{p}_{1})-\chi(\bm{p}_{2})\right)^{2}=\dot{\sigma}_{\text{macro}}+O(\varepsilon^{2}),$
(114)

where the $\dot{\sigma}_{\text{macro}}$ is the one computed in the linear
regime of irreversible thermodynamics, _i.e._ , after linearization of
${\left\langle\Pi\right\rangle}$ and ${\left\langle\bm{J}_{q}\right\rangle}$
in terms of $(\epsilon_{ij})$ and $\nabla_{\bm{x}}T$ as given in Eq. (49).

Third, we verify that in the presence of local equilibrium the entropy and
heat flux are proportional to each other. By substitutiong the Chapman–Enskog
solution (98) in Eq. (88) and using the expression of the Maxwell–Boltzmann
and the property (107b), we find

${\left\langle\bm{J}_{s}\right\rangle}=\underbrace{-\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\left(\bm{p}-\left\langle\bm{v}\right\rangle\right)f^{\text{MB}}(\ln
f^{\text{MB}}+\varepsilon\chi)}_{=0}-\varepsilon\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\left(\bm{p}-\left\langle\bm{v}\right\rangle\right)f^{\text{MB}}\chi\ln
f^{\text{MB}}+O(\varepsilon^{2})=\frac{1}{T}{\left\langle\bm{J}_{q}\right\rangle}+O(\varepsilon^{2})\,.$
(115)

The Chapman–Enskog procedure can in principle be continued to obtain higher
order corrections in the Knudsen number. However, the resulting Burnett and
super-Burnett equations can display unphysical oscillations or other
artificial behaviors near the boundary Cercignani (1988); Struchtrup and
Torrilhon (2007). Furthermore, by considering the higher perturbative orders
one looses the crucial property of local equilibrium Prigogine (1949); de
Groot and Mazur (1962). Therefore, the higher orders of approximation have
little significance for thermodynamics and will not be considered in the
following.

###  IV.4 Thermodynamics of the Transient Relaxation to Local Equilibrium

The Chapman–Enskog solution (98) describes only a restricted class of
solutions of the Boltzmann equation (93) close to the equilibrium for which
the entropy production rate (89) coincides with the macroscopic one in the
linear regime according to Eq. (114). However, the Chapman–Enskog perturbation
series does not describe the relaxation of an arbitrary initial condition
towards the distribution (98). Even if a complete theory of this relaxation
process for arbitrary initial conditions involves the solution of the
nonlinear Boltzmann equation, and it is therefore unattainable, we employ here
a perturbative treatment following Grad’s seminal contribution to this problem
Grad (1963); Cercignani (1988). In this way, we are able to identify the
contribution of this relaxation process to the EPR in Eq. (89).

We consider an initial condition $f_{1}(0)$ for the one-particle distribution
close to the near-equilibrium solution (98). Then, we introduce two separate
timescales. One timescale is associated with the (slow) dynamics of the
Chapman–Enskog solution (98). The other one, $\tau={\varepsilon^{-1}}t$, is
related to a fast relaxation process, and is relevant only over times
comparable to the Knudsen number (91). To account for both timescales, we add
a correction term to the Chapman–Enskog solution (98) and write

$\displaystyle f_{1}(t,\tau)=f_{1}^{\text{CE}}(t)+\varepsilon
f_{1}^{\text{rel}}(\tau)\,,$ (116)

with the condition that $f_{1}^{\text{rel}}(\tau)$ decays to zero
exponentially fast for $\tau\gg 1$. By substituting Eq. (116) in the Boltzmann
equation (93), we obtain

$\displaystyle\begin{split}&\varepsilon[(\partial_{t}+\bm{p}\cdot\nabla_{\bm{x}}+\bm{F}\cdot\nabla_{\bm{p}})f_{1}^{\text{CE}}(t)]-C(f_{1}^{\text{CE}},f_{1}^{\text{CE}})=\\\
&\,=-\varepsilon^{2}[(\partial_{t}+\bm{p}\cdot\nabla_{\bm{x}}+\bm{F}\cdot\nabla_{\bm{p}})f_{1}^{\text{rel}}(\tau)]+\varepsilon\left(C(f_{1}^{\text{rel}},f_{1}^{\text{CE}})+C(f_{1}^{\text{CE}},f_{1}^{\text{rel}})\right)+\varepsilon^{2}C(f_{1}^{\text{rel}},f_{1}^{\text{rel}})\,.\end{split}$
(117)

By introducing $\partial_{\tau}=\varepsilon\partial_{t}$, keeping only terms
linear in $\varepsilon$, and using the fact that $f_{1}^{\text{CE}}$ makes the
left-hand side of Eq. (117) vanish, being a solution of Eq. (93) up to order
$O(\varepsilon)$, we obtain

$\displaystyle\partial_{\tau}f_{1}^{\text{rel}}$
$\displaystyle=L\left(\frac{f_{1}^{\text{rel}}}{f^{\text{MB}}}\right)\,.$
(118)

Here we introduced the linearized collision operator

$\displaystyle L(g)=\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}_{1}^{\prime}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}_{2}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}_{2}^{\prime}\,W(\bm{p},\bm{p}_{2}|\bm{p}_{1}^{\prime},\bm{p}_{2}^{\prime})f^{\text{MB}}(\bm{p}_{1}^{\prime})f^{\text{MB}}(\bm{p}_{2}^{\prime})\left(g(\bm{p}_{1}^{\prime})+g(\bm{p}_{2}^{\prime})-g(\bm{p})-g(\bm{p}_{2})\right)\,,$
(119)

with $f^{\text{MB}}$ the Maxwell–Boltzmann distribution (96) satisfying
$f^{\text{MB}}(\bm{p}_{1}^{\prime})f^{\text{MB}}(\bm{p}_{2}^{\prime})=f^{\text{MB}}(\bm{p}_{1})f^{\text{MB}}(\bm{p}_{2})$
because of the conservation of the kinetic energy. The dynamics of
$f^{\text{rel}}_{1}(\tau)$ is obtained by solving Eq. (118), and we refer to
the original article Grad (1963) for the quantitative theory. Here, we focus
on the entropy production (89) for distributions of the type in Eq. (116). The
result in the dominant order in $\varepsilon$ is given by

$\dot{\sigma}_{B}\simeq\dot{\sigma}_{\text{macro}}+\dot{\sigma}_{\text{rel}}+\dot{\sigma}_{\text{interf}}\,,$
(120)

where $\dot{\sigma}_{\text{macro}}$ is given in Eq. (114),

$\displaystyle\dot{\sigma}_{\text{rel}}=\frac{\varepsilon}{4}\int\mathop{\mathrm{\mathstrut
d}}\nolimits^{4}\\!p\,W(\bm{p}_{1}^{\prime},\bm{p}_{2}^{\prime}|\bm{p}_{1},\bm{p}_{2})f^{\text{MB}}(\bm{p}_{1})f^{\text{MB}}(\bm{p}_{2})\left(\frac{f^{\text{rel}}_{1}(\bm{p}_{1}^{\prime})}{f^{\text{MB}}(\bm{p}_{1}^{\prime})}+\frac{f^{\text{rel}}_{1}(\bm{p}_{2}^{\prime})}{f^{\text{MB}}(\bm{p}_{2}^{\prime})}-\frac{f^{\text{rel}}_{1}(\bm{p}_{1})}{f^{\text{MB}}(\bm{p}_{1})}-\frac{f^{\text{rel}}_{1}(\bm{p}_{2})}{f^{\text{MB}}(\bm{p}_{2})}\right)^{2}\,,$
(121)

is the entropy production rate of the relaxation to the Chapman–Enskog
distribution (the local equilibrium distribution) and becomes negligible for
long times $\tau\gg 1$. It is important to note, however, that for short times
$\tau\lesssim 1$, the two contributions are of the same order in
$\varepsilon$. The third term in Eq. (120)

$\small\dot{\sigma}_{\text{interf}}=\frac{\varepsilon}{2}\int\mathop{\mathrm{\mathstrut
d}}\nolimits^{4}\\!p\,Wf^{\text{MB}}(\bm{p}_{1})f^{\text{MB}}(\bm{p}_{2})\left(\chi(\bm{p}_{1}^{\prime})+\chi(\bm{p}_{2}^{\prime})-\chi(\bm{p}_{1})-\chi(\bm{p}_{2}\right)\left(\frac{f^{\text{rel}}_{1}(\bm{p}_{1}^{\prime})}{f^{\text{MB}}(\bm{p}_{1}^{\prime})}+\frac{f^{\text{rel}}_{1}(\bm{p}_{2}^{\prime})}{f^{\text{MB}}(\bm{p}_{2}^{\prime})}-\frac{f^{\text{rel}}_{1}(\bm{p}_{1})}{f^{\text{MB}}(\bm{p}_{1})}-\frac{f^{\text{rel}}_{1}(\bm{p}_{2})}{f^{\text{MB}}(\bm{p}_{2})}\right)\,,$
(122)

does not have, in general, a definite sign. This means that while the total
EPR (121) is positive definite, it is not necessarily greater that the one due
to the macroscopic processes.

##  V Mesoscopic Scale: generalized Boltzmann equation with external
reservoirs

The microscopic and mesoscopic theory discussed in Section III and Section IV,
respectively, are derived starting from an isolated system setup. Therefore,
they do not explicitly account for what happens at the system’s boundary which
is crucial in open systems, as it can lead to thermodynamic inconsistencies
when arbitrary conditions are fixed at the boundary Gallavotti (2002);
Cercignani (1988); Gaspard (2022). By taking inspiration from both old and
recent works in nonequilibrium thermodynamics and stochastic thermodynamics
Bergmann and Lebowitz (1955); Lebowitz and Bergmann (1957); Van den Broeck and
Toral (2015); Horowitz and Esposito (2016), we develop here a model that
includes the boundary directly in the Boltzmann equation and examine the
corresponding thermodynamics. Note that throughout this section we will use
the adimensional variables introduced in Eq. (92).

###  V.1 Dynamics

The boundary $\partial\Omega$ is split into different reservoirs, or baths,
$\partial\Omega=\cup_{{b}}\partial\Omega_{{b}}$ like in § II.5.4. We assume
that each reservoir $b$ is composed of non-interacting particles always in
equilibrium at the temperature $T_{b}$ and its one-particle distribution
function is given by

$\displaystyle{f_{1}^{{b}}}(\bm{X},\bm{P})=\frac{N_{b}}{Z(T_{b})}\exp\left\\{-\frac{1}{T_{{b}}}\left(\frac{(\bm{P}-\bm{v}_{b})^{2}}{2M_{b}}+\mathcal{U}_{b}(\bm{X})\right)\right\\}\,.$
(123)

Here, $Z(T_{b})$ is the partition function, $N_{b}$ is the total number of
particles, $\bm{v}_{b}$ is the (externally controlled) velocity of the
reservoir, $\mathcal{U}_{b}(\bm{X})$ is a (arbitrary) potential that confines
the reservoir’s particles to $\partial\Omega_{{b}}$, _i.e._ ,
$N_{b}=\int_{\partial\Omega_{{b}}}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{X}\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{P}\,f_{1}^{{b}}(\bm{X},\bm{P})$. For simplicity, we further
assume that the particles of the reservoir have the same mass of those of the
system: $M_{b}=M=1$.

In the limit of dilute fluids, collisions are rare and can only take place
either between pairs of particles of the system or between one particle of the
system and one particle of a single bath. Thus, we modify the Boltzmann
equation (93) by adding collision terms corresponding to each reservoir:

$\displaystyle\partial_{t}f_{1}+\bm{p}\cdot\nabla_{{\bm{x}}}f_{1}+\bm{F}\cdot\nabla_{\bm{p}}f_{1}=\varepsilon^{-1}C(f_{1},f_{1})+{\varepsilon^{-1}}\sum_{b=1}^{B}C_{{b}}(f_{1})\,,$
(124)

where $C(f,f)$ is given in Eq. (73), while the linear collision term due to
the reservoirs reads

$\displaystyle C_{{b}}(f)=$ $\displaystyle\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}^{\prime}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{P}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{P}^{\prime}\,(W_{{b}}(\bm{p},\bm{P}|\bm{p}^{\prime},\bm{P}^{\prime})f(\bm{p}^{\prime})f_{1}^{{b}}(\bm{P}^{\prime})-W_{{b}}(\bm{p}^{\prime},\bm{P}^{\prime}|\bm{p},\bm{P})f(\bm{p})f_{1}^{{b}}(\bm{P}))\,.$
(125)

with the scattering rate $W_{{b}}$ satisfying Eq. (76) like $W$. Note that Eq.
(124) generalizes previous formulations of linear Boltzmann equations
discussed in the framework of stochastic thermodynamics Van den Broeck and
Toral (2015); Horowitz and Esposito (2016) by including the nonlinear
collision integral $C(f,f)$.

###  V.2 Thermodynamics

####  V.2.1 First and Second Law.

To derive the thermodynamics of the generalized Boltzmann equation (124), we
start by considering the first law of thermodynamics. To maintain the
presentation simple, without loosing essential elements in the thermodynamics,
we will consider only the case in which all the components of the boundary
move with a uniform velocity $\bm{v}_{B}$. The first step is to derive the
balance equations associated to the different microscopic quantities conserved
by the collision operator $C$. The total kinetic and potential energy, the
kinetic energy of the (uniform) motion of the domain $\Omega$ with velocity
$\bm{v}_{B}$, and the momentum satisfy

$\displaystyle\frac{1}{2}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\left\llangle\bm{p}^{2}\right\rrangle$
$\displaystyle=\left\llangle\bm{p}\cdot\bm{F}\right\rrangle+\varepsilon^{-1}{\sum_{b}}\int_{\Omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{x}\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\,\frac{\bm{p}^{2}}{2}\,C_{b}(f_{1})\,,$
(126a) $\displaystyle\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\left\llangle\phi\right\rrangle$
$\displaystyle=-\left\llangle\bm{p}\cdot\bm{F}\right\rrangle\,,$ (126b)
$\displaystyle\frac{1}{2}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\left\llangle\bm{v}_{\text{B}}^{2}\right\rrangle$
$\displaystyle=\frac{1}{2}\bm{v}^{2}_{B}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}N=0\,,$ (126c) $\displaystyle\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\left\llangle\bm{p}\right\rrangle$
$\displaystyle=\left\llangle\bm{F}\right\rrangle+\varepsilon^{-1}{\sum_{b}}\int_{\Omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{x}\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\,\bm{p}\,C_{b}(f_{1})\,,$ (126d)

respectively, where we used the property of summational invariance (78) for
the collision operator $C$, $\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}\,C_{b}(f,f_{b}){=0}$, and $f_{1}=0$ on the boundaries of
the integration domain. Eqs. (126a)-(126d) allow us to write the balance for
the (global) internal energy
$U=\left\llangle(\bm{p}-\bm{v}_{B})^{2}\right\rrangle/2$ in terms of the power
$\dot{W}$ and heat flow at the boundaries $\dot{Q}_{b}$ as

$\displaystyle\mathop{\mathrm{\mathstrut d}}\nolimits\\!_{t}U$
$\displaystyle=\frac{1}{2}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\left\llangle(\bm{p}-\bm{v}_{B})^{2}\right\rrangle$ (127a)
$\displaystyle=\left\llangle(\bm{p}-\bm{v}_{B})\cdot\bm{F}\right\rrangle+\varepsilon^{-1}\sum_{b}\int_{\Omega{(t)}}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{x}\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\left(\frac{\bm{p}^{2}}{2}-\bm{v}_{\text{B}}\cdot\bm{p}\right)C_{b}(f_{1})$
(127b) $\displaystyle=\underbrace{-\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\left\llangle\phi\right\rrangle-\left\llangle\bm{v}_{B}\cdot\bm{F}\right\rrangle}_{\equiv\dot{W}}+{\sum_{b}}\underbrace{\varepsilon^{-1}\int_{\Omega{(t)}}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{x}\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\left(\frac{\bm{p}^{2}}{2}-\bm{v}_{\text{B}}\cdot\bm{p}\right)C_{b}(f_{1})}_{\equiv\dot{Q}_{b}}\,.$
(127c)

The work rate $\dot{W}$ identified here is the sum of an exact time
derivative, which gives the rate of change the potential energy due to the
changes in the particle distribution inside the volume, and which vanishes in
steady conditions, plus a macroscopic term
$\left\llangle\bm{v}_{B}\cdot\bm{F}\right\rrangle$ which is the power
externally provided to maintain the velocity $\bm{v}_{B}$ at the boundary.
Before turning to the entropy balance, we first recognize that

$\displaystyle
f_{b}^{\text{eq}}(\bm{x},\bm{p})=\exp\left\\{-\frac{1}{T_{b}}\left(\frac{(\bm{p}-\bm{v}_{B})^{2}}{2}-\mu_{\text{eq}}(\bm{x})\right)\right\\}\,,$
(128)

is the equilibrium solution of Eq. (124) if all the baths have the same
temperature $T_{b}=T_{B}$ $\forall{b}$. Indeed, it satisfies

$D_{t}{f_{b}^{\text{eq}}}(\bm{x},\bm{p})=0\,,$ (129)

where $D_{t}=(\partial_{t}+\bm{v}_{B}\cdot\nabla_{\bm{x}})$, if the chemical
potential $\mu^{\text{eq}}(\bm{x})$ satisfies
$\nabla_{\bm{x}}\mu^{\text{eq}}(\bm{x})=\bm{F}$ or, equivalently,

$\mu_{\text{eq}}(\bm{x})=-\phi(\bm{x})+{\overline{\mu}}_{\text{eq}}\,,$ (130)

with ${\overline{\mu}}_{\text{eq}}$ a constant parameter. This can be easily
proven by direct substitution and recalling that the collision terms vanish
because $f_{b}^{\text{eq}}$ is a function of summational invariants. This
physically means that the equilibrium distribution has the temperature imposed
by the reservoir and the average velocity of the fluid must coincide with the
one externally imposed by the moving boundary of the reservoir. Using the
expression of the equilibrium distribution (128) and
$\int\mathop{\mathrm{\mathstrut d}}\nolimits\\!\bm{p}\,C_{b}(f_{1})=0$, we can
write the heat flux as
$\dot{Q}_{b}=-\varepsilon^{-1}T_{b}\int_{\Omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{x}\int\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}\,\ln(f_{b}^{\text{eq}})\,C_{b}(f_{1})/(2\pi\hbar)^{3}$.
Therefore, we can write the global entropy balance as

$\displaystyle\small\mathop{\mathrm{\mathstrut d}}\nolimits\\!_{t}S$
$\displaystyle=-\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\int_{\Omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{x}\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}(f_{1}\ln f_{1}-f_{1})$ (131a)
$\displaystyle=-\varepsilon^{-1}\int_{\Omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{x}\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\left(C(f_{1},f_{1})\ln
f_{1}+{\sum_{b}}C_{b}(f_{1})\ln\frac{f_{1}}{f^{\text{eq}}_{b}}+{\sum_{b}}C_{b}(f_{1})\ln
f^{\text{eq}}_{b}\right)$ (131b)
$\displaystyle=\dot{\Sigma}_{\text{int}}+{\sum_{b}}\left(\dot{\Sigma}_{b}+\frac{\dot{Q}_{b}}{T_{b}}\right)\,.$
(131c)

The second equality is again obtained after integrating by parts, summing and
subtracting $D_{t}f_{1}\ln f_{b}^{\text{eq}}$ with $f_{b}^{\text{eq}}$ in Eq.
(128), and using the summational invariance (78). In the last step we
introduced the entropy production rates inside the system

$\displaystyle\dot{\Sigma}_{\text{int}}=$
$\displaystyle-\varepsilon^{-1}\int_{\Omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{x}\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\,C(f_{1},f_{1})\ln f_{1}\,,$ (132)

and the one at the boundary

$\displaystyle\dot{\Sigma}_{b}=$
$\displaystyle-\varepsilon^{-1}\int_{\Omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{x}\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\,C_{{b}}(f_{1})\ln\frac{f_{1}}{f_{b}^{\text{eq}}}\,,$
(133)

which are nonnegative as we now show. By using indistinguishability of the
particles and the symmetry (76), we find that Eq. (132) can be rewritten as

$\displaystyle\dot{\Sigma}_{\text{int}}=$
$\displaystyle\frac{\varepsilon^{-1}}{4}\int_{\Omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{x}\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}_{1}^{\prime}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}_{2}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}_{2}^{\prime}}{(2\pi\hbar)^{3}}\,W(\bm{p},\bm{p}_{2}|\bm{p}_{1}^{\prime},\bm{p}_{2}^{\prime})(f_{1}(\bm{p}_{1}^{\prime})f_{1}(\bm{p}_{2}^{\prime})-f_{1}(\bm{p})f_{1}(\bm{p}_{2}))\ln\frac{f_{1}(\bm{p}_{1}^{\prime})f_{1}(\bm{p}_{2}^{\prime})}{f_{1}(\bm{p})f_{1}(\bm{p}_{2})}\geq
0\,.$ (134)

This approach cannot be applied to Eq. (133) since the particles of the system
and the particles of the reservoirs are not indistinguishable. However, by
using the explicit expression of $f_{1}^{b}$ and $f_{b}^{\text{eq}}$ in Eqs.
(123) and (128), respectively, together with the conservation of the kinetic
energy, _i.e._ ,
$(\bm{p}-\bm{v}_{B})^{2}-(\bm{p}^{\prime}-\bm{v}_{B})^{2}=(\bm{P}^{\prime}-\bm{v}_{B})^{2}-(\bm{P}-\bm{v}_{B})^{2}$,
we obtain

$\displaystyle\dot{\Sigma}_{b}=$
$\displaystyle\frac{\varepsilon^{-1}}{2}\int_{\Omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{x}\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{P}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}^{\prime}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{P}^{\prime}}{(2\pi\hbar)^{3}}\,W_{{b}}(\bm{p},\bm{P}|\bm{p}^{\prime},\bm{P}^{\prime})(f_{1}(\bm{p}^{\prime})f_{1}^{{b}}(\bm{P}^{\prime})-f_{1}(\bm{p})f_{1}^{{b}}(\bm{P}))\ln\frac{f_{1}(\bm{p}^{\prime})f_{b}^{\text{eq}}(\bm{p})}{f_{1}(\bm{p})f_{b}^{\text{eq}}(\bm{p}^{\prime})}\,,$
(135a) $\displaystyle=$
$\displaystyle\frac{\varepsilon^{-1}}{2}\int_{\Omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{x}\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{P}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}^{\prime}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{P}^{\prime}}{(2\pi\hbar)^{3}}\,W_{{b}}(\bm{p},\bm{P}|\bm{p}^{\prime},\bm{P}^{\prime})(f_{1}(\bm{p}^{\prime})f_{1}^{{b}}(\bm{P}^{\prime})-f_{1}(\bm{p})f_{1}^{{b}}(\bm{P}))\ln\frac{f_{1}(\bm{p}^{\prime})f_{1}^{{b}}(\bm{P}^{\prime})}{f_{1}(\bm{p})f_{1}^{{b}}(\bm{P})}\geq
0\,.$ (135b)

####  V.2.2  Global thermodynamic potential.

We now show that, when the system is i) closed in a domain $\Omega(t)$ with
the boundary moving at constant and uniform velocity ${\bm{v}_{B}}$, and ii)
in contact with a single external reservoir (identified by $b=B$) with
temperature $T_{B}$, the dynamical equation (124) satisfies the analogous of a
H-theorem (_i.e._ , it admits a Lyapunov function) ruling the relaxation to
equilibrium.

To do so, we use the nonequilibrium Massieu potential $\Phi$ defined in Eq.
(65) (§ III.2). In Eq. (67), we proved that for a closed system $\Phi$ is
upper bounded by its equilibrium value $\Phi^{\text{eq}}$. This result
specializes here by using $f_{B}^{\text{eq}}$ given in Eq. (128) (with
${\overline{\mu}}_{\text{eq}}$ introduced in Eq. (130)), and the Hamiltonian
$H_{1}=(\bm{p}-\bm{v}_{B})^{2}/{2}+{\phi(\bm{x})}$ accounting for the moving
boundary. This nonequilibrium Massieu potential $\Phi$ plays the role of a
Lyapunov function. Indeed, for every initial $f_{1}$,

$\displaystyle\mathop{\mathrm{\mathstrut d}}\nolimits\\!_{t}\Phi(t)$
$\displaystyle=-\int_{\Omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{x}\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}D_{t}\left(f_{1}\ln\frac{f_{1}}{f_{B}^{\text{eq}}}-f_{1}+f_{B}^{\text{eq}}\right)$
(136a)
$\displaystyle=-\varepsilon^{-1}\int_{\Omega(t)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{x}\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\left(C(f_{1},f_{1})+C_{B}(f_{1},f_{1}^{B})\right)\ln\frac{f_{1}}{f_{B}^{\text{eq}}}$
(136b) $\displaystyle=\dot{\Sigma}_{\text{int}}+\dot{\Sigma}_{B}\geq 0\,,$
(136c)

where we used Eqs. (67), (167), and (124); $\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\Phi^{\text{eq}}=0$ since $D_{t}f_{B}^{\text{eq}}=0$;
$\int\mathrm{d}\bm{p}\,C(f_{1},f_{1})\ln f_{B}^{\text{eq}}=0$ since
$f_{B}^{\text{eq}}$ is a function of summational invariants; $f_{1}=0$ on the
boundaries of the integration domain; the definitions (132) and (133) in the
last step.

Finally, by using the definition of the nonequilibrium Massieu potential in
Eq. (65) (with the Hamiltonian
$H_{1}=(\bm{p}-\bm{v}_{B})^{2}/{2}+{\phi(\bm{x})}$) together with Eqs. (127c)
and (131c), we compute $\mathop{\mathrm{\mathstrut d}}\nolimits\\!_{t}\Phi$
for the general case of multiple reservoirs $b=0,...,B$. By choosing the
$B$-th reservoir as a reference and isolating the total EPR, we obtain

$\dot{\Sigma}_{\text{tot}}\equiv\dot{\Sigma}_{\text{int}}+\sum_{{b}=1}^{B-1}\dot{\Sigma}_{b}=\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\Phi-\frac{1}{T_{B}}\underbrace{\left\llangle\bm{v}_{B}\cdot\bm{F}\right\rrangle}_{\dot{W}_{\text{macro}}}+\sum_{b}\dot{Q}_{b}\left(\frac{1}{T_{B}}-\frac{1}{T_{b}}\right)\,,$
(137)

that is the mesoscopic analogue of Eq. (46). Indeed, it features the time
derivative of a thermodynamic potential $\mathop{\mathrm{\mathstrut
d}}\nolimits\\!_{t}\Phi$; a work term $\dot{W}_{\text{macro}}$ done against
the conservative force $\bm{F}=-\nabla_{\bm{x}}\phi$ via a protocol externally
imposed by the velocity of the boundary; a _flux-force_ contribution due to
the presence of heat fluxes $\dot{Q}_{b}$ exchanged with the environment and
triggered by the presence of global differences in temperature. The
nonconservative work of shearing and deformation in Eq. (46) vanishes here
because we assumed $\bm{v}_{B}$ to be uniform on the boundary. The same
happens to the terms in Eq. (46) that account for explicit time-dependence of
the externally controlled parameters. Therefore, we can identify the
nonequilibrium Massieu potential defined at the mesoscopic scale with the one
featuring in Eq. (46), concluding $\Phi=Y$, if the one particle distribution
$f_{1}$ satisfies local equilibrium.

####  V.2.3 Local thermodynamic potential.

We now identify the mesoscopic thermodynamic potential $\left\langle
y\right\rangle(\bm{x})$, upper bounded by its equilibrium value $\left\langle
y\right\rangle_{\text{eq}}$, corresponding to the macroscopic thermodynamic
potential $Y$ given in Eq. (45). We start by observing that

$\small 0\leq\int\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\bm{p}}{(2\pi\hbar)^{3}}\left(f_{1}(\bm{x},\bm{p})\ln\frac{f_{1}(\bm{x},\bm{p})}{f_{{B}}^{\text{eq}}(\bm{x},\bm{p})}+f_{B}^{\text{eq}}(\bm{x},\bm{p})-f_{1}(\bm{x},\bm{p})\right)=-{\left\langle\rho\right\rangle}{\left\langle
s_{\text{B}}\right\rangle}+{\left\langle\rho\right\rangle}_{\text{eq}}+\frac{{\left\langle\rho\right\rangle}}{T_{B}}\left(\frac{({\left\langle\bm{v}\right\rangle}-\bm{v}_{B})^{2}}{2}+\phi+{\left\langle
u\right\rangle}{-}{\overline{\mu}}_{\text{eq}}\right)\,,$ (138)

by using Eqs. (79), (83), (86), (128), (130), and
$(\bm{p}-\bm{v}_{B})^{2}=(\bm{p}-{\left\langle\bm{v}\right\rangle})^{2}+2\bm{p}\cdot({\left\langle\bm{v}\right\rangle}-\bm{v}_{B})+{\left\langle\bm{v}\right\rangle}^{2}+(\bm{v}_{B})^{2}$.
We then identify the spatio-temporal field representing the total energy in
the comoving frame of the system as

${\left\langle
e\right\rangle}\equiv\frac{({\left\langle\bm{v}\right\rangle}-\bm{v}_{B})^{2}}{2}+\phi+{\left\langle
u\right\rangle}\,,$ (139)

and recall that pressure tensor satisfies Eq. (101) when $f_{1}$ is a
Maxwell–Boltzmann distribution like $f_{B}^{\text{eq}}$ in Eq. (128) and
therefore

${\left\langle\rho\right\rangle}_{\text{eq}}=\frac{p_{B}}{T_{B}}\,,$ (140)

where $p_{B}$ is the externally imposed pressure that the system reaches at
equilibrium. This allows us to write the mescoscopic thermodynamic potential
as (minus) the right-hand side of Eq. (138) divided by the local density
${\left\langle\rho\right\rangle}$:

$\left\langle y\right\rangle={\left\langle
s_{\text{B}}\right\rangle}-\frac{{\left\langle
e\right\rangle}+p_{B}{\left\langle\rho\right\rangle}^{-1}}{T_{B}}{+}\frac{{\overline{\mu}}_{\text{eq}}}{T_{B}}\,,$
(141)

with ${\left\langle\rho\right\rangle}^{-1}$ being the local specific volume as
in § II.3. Note that $\left\langle y\right\rangle_{\text{eq}}=0$ because of
the local equilibrium (105) and $\left\langle y\right\rangle-\left\langle
y\right\rangle_{\text{eq}}\leq 0$ because of Eq. (138). Furthermore,
$\left\langle y\right\rangle$ is the mesoscopic formulation of the macroscopic
$Y$ if $f_{1}$ satisfies local equilibrium, like the Chapman–Enskog solution
(98), that turns ${\left\langle s_{\text{B}}\right\rangle}$ into the specific
entropy. The only difference between the local thermodynamic potential (141)
and the macroscopic one (45) due to the term
${\overline{\mu}}_{\text{eq}}/T_{B}$ emerges because the system is locally
open as already observed in Ref. Avanzini  _et al._ (2019).

##  VI Stochastic thermodynamics for viscous fluids

By building on Eq. (124), we are in the position to develop a linear theory
for one-particle distribution function, endowed with the property of local
detailed balance line in stochastic thermodynamics. This theory further admits
a Chapman–Enskog solution identical to the one obtained in Subs. IV.3 and