# Multiscale heat transport with inertia and thermal vortices

Michal Pavelka
Mathematical Institute, Faculty of Mathematics, Charles University,
Sokolovská 83, 18675 Prague, Czech Republic Liliana Restuccia
Department of Mathematical and Computer Sciences, Physical Sciences and Earth
Sciences University of Messina,
Viale F. Stagno d’Alcontres, Salita Sperone 31, 98166, Messina, Italy David
Jou
Department of Physics, Universitat Autonoma de Barcelona
08193 Bellaterra, Catalonia, Spain

###### Abstract

In this paper, we present a Hamiltonian and thermodynamic theory of heat
transport on various levels of description. Transport of heat is formulated
within kinetic theory of polarized phonons, kinetic theory of unpolarized
phonons, hydrodynamics of polarized phonons, and hydrodynamics of unpolarized
phonons. These various levels of description are linked by Poisson reductions,
where no linearizations are made. Consequently, we obtain a new phonon
hydrodynamics that contains convective terms dependent on vorticity of the
heat flux, which are missing in the standard theories of phonon hydrodynamics.
Moreover, the equations are hyperbolic and Galilean invariant, unlike current
theories for beyond-Fourier heat transport. The vorticity-dependent terms
violate the alignment of the heat flux with the temperature gradient even in
the stationary state, which is expressed by a Fourier-Crocco equation. The new
terms also cause that temperature plays in heat transport a similar role as
pressure in aerodynamics.

###### Contents

  1. 1 Introduction
    1. 1.1 Introduction to GENERIC
  2. 2 Kinetic theory of phonons
    1. 2.1 Kinetic theory of phonons with polarization
    2. 2.2 Overall kinetic theory of phonons
  3. 3 Overall hydrodynamics of phonons
    1. 3.1 Reversible evolution
    2. 3.2 Zero-th order irreversible evolution
    3. 3.3 First-order irreversible evolution
    4. 3.4 Galilean invariance, inertia, and thermal vortices
    5. 3.5 Reduction to Fourier heat transport
    6. 3.6 A hypothetical dependence of energy on vorticity
  4. 4 Conclusion
  5. A Maximum entropy principle for energy and entropy
  6. B Simplification to the Cattaneo equation
  7. C Galilean invariance of phonon hydrodynamics
  8. D Reversibility, irreversibility, and Onsager-Casimir reciprocal relations
  9. E Details of the Chapman-Enskog approximation

##  1 Introduction

This manuscript contains a new theory of phonon hydrodynamics with corrected
behavior of thermal vortices. Phonon hydrodynamics is of particular importance
for instance in nanoscale heat transport in semiconductors, thermal diodes, or
two-dimensional materials [32, 87, 81].

Heat transport in solids is represented by transport of energy by lattice
vibrations [49, 73]. Since the lattice vibrations can be seen as
quasiparticles obeying Hamilton canonical equations (phonons), evolution
equation of their distribution function obeys kinetic theory [72, 63], which
is also a Hamiltonian system [44, 24]. Here, we start with kinetic theory of
phonons with three polarizations (two transverse and one longitudinal). This
description is then reduced to hydrodynamics with momenta and entropies of the
particular polarizations, and to the overall kinetic theory of phonons
(disregarding polarizations). Both these levels of description are then
further reduced to a new theory of phonon hydrodynamics, where state variables
are the overall momentum density and entropy density. Diagram 1 summarizes the
various levels of description and connections between them.

|
---|---
$\textstyle{f_{T+},f_{T-},f_{L}\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{\mathbf{m}_{T+},\mathbf{m}_{T-},\mathbf{m}_{T_{L}},s_{T+},s_{T-},s_{T_{L}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{f\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{s,\mathbf{m}}$
Figure 1:  Distribution functions of positively and negatively transversely
and longitudinally polarized phonons, $f_{T+},f_{T-},f_{L}$, can be summed to
the overall distribution function $f$. They can be also projected to the
momentum and entropy densities of the particular polarizations,
$\mathbf{m}_{T+}$, $\mathbf{m}_{T-}$, $\mathbf{m}_{T_{L}}$, $s_{T+}$,
$s_{T-}$, and $s_{T_{L}}$. These densities can be then summed to the overall
momentum density and entropy density of phonons, $\mathbf{m}$ and $s$.
Finally, the overall distribution function $f$ can be projected to
$\mathbf{m}$ and $s$. The projection of the Poisson bracket of kinetic theory
to the hydrodynamic fields then leads to a new theory of phonon hydrodynamics,
similarly as in [30].

A particular feature of our phonon hydrodynamics is that the evolution
equation for the conjugate entropy flux $\mathbf{w}=\mathbf{m}/s$ (the ratio
of phonon momentum and entropy, see Equation (27)) contains terms depending on
the thermal vorticity, $\boldsymbol{\omega}=\nabla\times\mathbf{w}$. In
contrast with other approaches, our theory consistently describes the motion
of thermal vorticity because we neither linearise the reversible part of the
equations nor are we restricted to one-dimensional systems [32, 35, 34, 50,
16]. The terms that depend on $\boldsymbol{\omega}$ make the equations
Galilean invariant, unlike the usual Maxwell-Cattaneo-Vernotte or Guyer-
Krumhansl equations [11, 35]. Although evidence for the presence of thermal
vorticity in heat transport is provided by kinetic theory [82, 91, 75],
standard phonon hydrodynamics can not correctly propagate such vortices
because it lack convective terms. Moreover, the missing convective terms are
the reason why the standard theories of phonon hydrodynamica violate Galilean
invariance even when all collision terms are neglected. The convective terms
are lost in the standard linearized Chapman-Enskog reduction, while our method
(based on reduction the underlying Poisson brackets) does not linearize the
reversible evolution and recovers the convective terms. In other words, phonon
hydrodynamics presented in this manuscript corrects the standard phonon
hydrodynamics so that heat vortices are transported in agreement with the
kinetic theory.

Moreover, the new terms cause that the heat flux no longer needs to be aligned
with the gradient of temperature even in the stationary state. Instead, the
stationary state is described by a Fourier-Crocco equation (49), where the
heat flux, thermal vorticity, and temperature gradient can form a triplet of
mutually perpendicular vectors in the limiting case of negligible dissipation.
Our equations also represent a system of quasilinear first-order hyperbolic
equations, so it describes hyperbolic heat conduction [74, 23, 76, 17] (if we
neglect the viscous-like dissipative terms). In the one-dimensional
approximation and after linearization, our theory simplifies to the standard
Maxwell-Cattaneo-Vernotte equation [11]. Finally, when we go beyond the zero-
th Chapman-Enskog approximation in the irreversible evolution of phonon
momentum density, we obtain higher-order (Laplacian) terms that correspond to
the dissipation in the Guyer-Krumhansl equation. Therefore, our phonon
hydrodynamics can be seen as a non-linear generalization of the (linear)
Guyer-Krumhansl equation [35]. Dynamics of thermal vorticity based on
solutions of the phonon kinetic theory was already observed in [91, 75], but
here we also formulate the corresponding phonon hydrodynamics.

A theory of hyperbolic heat transport was developed in the context of Extended
Rational Thermodynamics (RET) [77, 16, 63], including comparison with the
heat-pulse experiment [57]. In the present manuscript, we show the Hamiltonian
structure of the theory without linearization, which keeps also the non-linear
vorticity-dependent terms, that are important for inertial effects, for
effects of thermal vorticity, and for Galilean invariance.

Another description of hyperbolic heat transport was developed within the
Internal Variables Theory (IVT) [50, 5]. This theory is not constructed by
reduction from a more microscopic description (in contrast with kinetic
theory), but by extending the Fourier heat equation by including additional
internal variables and parameters, which remain unspecified. IVT shows a
better agreement with the experimental data than the kinetic theory for the
same number of state variables, but it also includes more fitting parameters.
Comparison of IVT with the GENERIC framework (used in this paper) can be found
in [84].

Apart from RET as IVT, there are other theories of hyperbolic heat transport
coming from non-equilibrium thermodynamics (NET). NET was actually founded by
Boltzmann [20], Onsager [64, 65], and others [58, 14, 18]. Subsequent theories
of NET, as Rational (Extended) Thermodynamics [88, 63], or Extended
Irreversible Thermodynamics [45, 46], bring the possibility to derive
evolution equations for additional state variables. There is plenty of
successful applications of these theories, ranging from non-Newtonian fluids
and solids to theories of charged mixtures [56]. However, it is often
difficult to fully determine the reversible terms that are not visible in the
formula for entropy production, which is why we use a different framework,
that builds upon Hamiltonian mechanics and has already been applied to
radiation hydrodynamics [29]. Here, we rely on the General Equation for Non-
Equilibrium Reversible-Irreversible Coupling (GENERIC) [31, 67, 66, 69], which
combines Hamiltonian continuum mechanics and with gradient dynamics.

The plan of this paper follows. Section 2 recalls kinetic theory of polarized
phonons, in particular the evolution of the distribution function of the
phonon quasiparticles. Section 2.2 then simplifies the description to the
overall distribution of phonons, disregarding the polarizations, and it adds
irreversible evolution in the form of the Callaway model (collisions of
phonons between themselves and with impurities of the lattice). Section 3
represents the main part of this paper as it contains the phonon hydrodynamics
that emerges by reduction of the kinetic theory. The reversible part of the
hydrodynamic theory is expressed by Equations (26b) or (30a) while their full
form (in the zero-th Chapman-Enskog approximation) is in Equations (36).
Equations (44) then show phonon hydrodynamics within the first Chapman-Enskog
approximation, where higher-order dissipative terms (Laplacian-like) appear.
The Fourier-Crocco equation (49)relates thermal vorticity, temperature
gradient, and heat flux, which leads to a similarity between temperature in
heat transfer and pressure in aerodynamics (Figure 3(b)). In order to obtain
the enhanced precision, we rely on geometric mechanics and Poisson reductions,
which represent a cornerstone of the GENERIC framework, summarized in the
following Section.

###  1.1 Introduction to GENERIC

The General Equation for Non-Equilibrium Reversible-Irreversible Coupling
(GENERIC) framework combines (reversible) Hamiltonian evolution with
(irreversible) gradient dynamics. The framework was established in [31, 67]
and further developments can be found in monographs [66, 69].111It is
equivalent with the metriplectic systems if the dissipation potential is
considered quadratic [62]. However, in cases with non-quadratic dissipation
potentials, both frameworks are different [28]. The GENERIC framework has also
a precursor in the so called one-generator formalism [6] and in the nonlinear
Onsager-Casimir (NOC) equations [26]. An advantage of GENERIC is that it
determines the reversible part of the evolution of given state variables
because it better exploits the results of geometrical mechanics. In the rest
of this Section, we shall recall the GENERIC framework.

Let $\mathbf{x}$ be a set of state variables and let $A(\mathbf{x})$ be a
functional on the space (or manifold) where $\mathbf{x}$ plays the role of
coordinates. Evolution of that functional is then given within GENERIC in two
equivalent forms,

$\dot{x}^{i}=L^{ij}\frac{\delta E}{\delta x^{j}}+\frac{\delta\Xi}{\delta
x^{*}_{i}}\Big{|}_{\mathbf{x}^{*}=\frac{\delta
S}{\delta\mathbf{x}}}\quad\text{or}\quad\dot{A}=\\{A,E\\}+\left\langle\frac{\delta
A}{\delta x^{i}},\frac{\delta\Xi}{\delta
x^{*}_{i}}\right\rangle\big{|}_{\mathbf{x}^{*}=\frac{\delta
S}{\delta\mathbf{x}}},$ (1)

where the latter is a sort of weak formulation. The evolution consists of two
parts, a reversible one and an irreversible one. Note that
$\delta\bullet/\delta\bullet$ is the functional derivative, see [66, 69, 6],
$\langle\bullet,\bullet\rangle$ is a duality (also represented by contraction
over repeated indexes), $E$ is the total energy of the isolated system under
consideration, $\Xi(\mathbf{x},\mathbf{x}^{*})$ is a dissipation potential,
and $S$ is the total entropy. We denote by $\mathbf{x}^{*}$ the dual variables
to the state variables $\mathbf{x}$, which are also called conjugate entropic
variables, eventually identified with derivatives of the total entropy $S$.

If the dissipation potential is quadratic,

$\Xi=\frac{1}{2}\left\langle
x^{*}_{i},M^{ij}(\mathbf{x}),x^{*}_{j}\right\rangle,$ (2)

the GENERIC evolution (1) can be rewritten as

$\partial_{t}x^{i}=L^{ij}E_{x^{j}}+M^{ij}S_{x^{j}},$ (3)

where the symmetric operator $M^{ij}$ (called a dissipative matrix) works as
in the Ginzburg-Landau theory [21, 41].

The reversible part is generated by a Poisson bivector $\mathbf{L}$ or a
Poisson bracket $\\{\bullet,\bullet\\}$, which are related through

$\\{A,B\\}=\left\langle dA,L^{ij},dB\right\rangle=\frac{\delta A}{\delta
x^{i}}L^{ij}\frac{\delta B}{\delta x^{j}}\quad\text{and}\quad
L^{ij}=\\{x^{i},x^{j}\\}.$ (4)

The bivector is a twice contravariant antisymmetric tensor field on the space
of the state variables, which is Lie-dragged by the Hamiltonian evolution
[19]. The latter property then ensures validity of the Jacobi identity,

$\\{A,\\{B,C\\}\\}+\\{B,\\{C,A\\}\\}+\\{C,\\{A,B\\}\\}=0\quad\forall A,B,C$
(5)

while the antisymmetry of $\mathbf{L}$ gives antisymmetry of the bracket,

$\\{A,B\\}=-\\{B,A\\}\quad\forall A,B.$ (6)

A particular implication of the antisymmetry is conservation of energy,
$\dot{E}=\\{E,E\\}=0$, by the Hamiltonian mechanics. Finally, an equivalent
way to express the reversible part of the evolution is

$(\partial_{t}\mathbf{x})_{\mathrm{rev}}=\\{\mathbf{x},E\\},$ (7)

as can be verified by the first relation of (4). Once we have recalled the
GENERIC framework, we are now in position to formulate kinetic theory of
phonons within that framework, which is the subject of the following Section.

##  2 Kinetic theory of phonons

Phonons are quasiparticles obtained by diagonalization of the Hamiltonian of a
crystal [52], and their motion is described by Hamilton canonical equations
[49], constructed from the canonical Poisson bracket and from the energy of
phonons. The purpose of this Section is to formulate Hamiltonian kinetic
theory that governs the distribution functions of the phonons. Actually,
phonons have three possible polarizations, two transverse ($T+$ and $T-$) and
one longitudinal. Therefore, there are three distribution functions
$f_{T+}(t,\mathbf{r},\mathbf{p})$, $f_{T-}(t,\mathbf{r},\mathbf{p})$, and
$f_{L}(t,\mathbf{r},\mathbf{p})$.

###  2.1 Kinetic theory of phonons with polarization

Similarly as the evolution of the distribution function of classical particles
is Hamiltonian, and thus generated by a Poisson bracket and energy, evolution
of the phononic distribution functions is also Hamiltonian, and generated by
the sum of the Boltzmann Poisson brackets

$\\{A,B\\}^{(f_{T+,T-,L})}=\sum_{\alpha\in\\{T+,T-,L\\}}\int d\mathbf{r}\int
d\mathbf{p}f_{\alpha}\cdot\left(\frac{\partial
A_{f_{\alpha}}}{\partial\mathbf{r}}\cdot\frac{\partial
B_{f_{\alpha}}}{\partial\mathbf{p}}-\frac{\partial
B_{f_{\alpha}}}{\partial\mathbf{r}}\cdot\frac{\partial
A_{f_{\alpha}}}{\partial\mathbf{p}}\right),$ (8)

where $A$ and $B$ are two arbitrary functionals of the distribution functions
[25, 69]. Symbols $A_{f_{\alpha}}$ stands for the functional derivative of $A$
with respect to $f_{\alpha}$. The reversible evolution equations implied by
this Poisson bracket are (using Equation (7))

$(\partial_{t}f_{\alpha})_{\mathrm{rev}}=-\frac{\partial
f_{\alpha}}{\partial\mathbf{r}}\cdot\frac{\partial
E_{f_{\alpha}}}{\partial\mathbf{p}}+\frac{\partial
f_{\alpha}}{\partial\mathbf{p}}\cdot\frac{\partial
E_{f_{\alpha}}}{\partial\mathbf{r}}\quad\forall\alpha\in\\{T+,T-,L\\}.$ (9)

What is the energy $E(f_{T+},f_{T-},f_{L})$? In general, we need the whole
band structure of the crystal to find the precise formula for energy [49].
However, if we only restrict to low temperatures, to the first Brillouin zone,
and acoustic phonons, we can approximate the energy as

$E=\int d\mathbf{r}\int d\mathbf{p}c|\mathbf{p}|(f_{T+}+f_{T-}+f_{L}),$ (10)

where $c$ is the Debye velocity and $|\mathbf{p}|$ is the Eulerian norm of the
momentum [15].

The evolution equations obtained from a Poisson bracket are reversible both in
the sense of the time-reversal transformation [68] and in the sense that they
preserve the bosonic entropy [51],

$S^{(f_{T+,T-,L})}=\sum_{\alpha}S^{\text{kinetic},\alpha}(f_{\alpha})\\\
\quad\text{with}\quad
S^{\text{kinetic},\alpha}(f_{\alpha})=-\frac{k_{B}}{h^{3}}\int d\mathbf{r}\int
d\mathbf{p}\left(h^{3}f_{\alpha}\ln(h^{3}f_{\alpha})-(1+h^{3}f_{\alpha})\ln(1+h^{3}f_{\alpha})\right).$
(11)

Indeed, any integral of a function of the distribution functions is a Casimir
of the Poisson bracket, that is $\\{S,B\\}^{(f_{T+,T-,L})}=0$ $\forall B$,
which can be used also to stability analysis of the Hamiltonian equations
[40].

Instead of the distribution functions, which contain a large amount of
information, we can reduce the description to the hydrodynamic fields of the
distribution functions, namely the momentum densities and entropy densities

$\displaystyle\mathbf{m}_{\alpha}$ $\displaystyle=\int
d\mathbf{p}f_{\alpha}\mathbf{p}$ (12a) $\displaystyle s_{\alpha}$
$\displaystyle=\int d\mathbf{p}\eta_{\alpha}(f_{\alpha}),$ (12b) where
$\eta_{\alpha}$ stands for a real-valued function, for instance the integrands
in entropy (11).

The Poisson bracket for the hydrodynamic fields is obtained by plugging
functionals that depend only on the hydrodynamic fields,
$A(\mathbf{m}_{\alpha}(f_{\alpha}),s(f_{\alpha}))$, to the Poisson bracket
(8). The resulting Poisson bracket reads

$\\{A,B\\}^{(\mathbf{m}_{\alpha},s_{\alpha})}=\sum_{\alpha}\int
d\mathbf{r}s_{\alpha}\left(\partial_{i}A_{s_{\alpha}}B_{m_{\alpha
i}}-\partial_{i}B_{s_{\alpha}}A_{m_{\alpha i}}\right)+\sum_{\alpha}\int
d\mathbf{r}m_{\alpha i}\left(\partial_{j}A_{m_{\alpha i}}B_{m_{\alpha
j}}-\partial_{j}B_{m_{\alpha i}}A_{m_{\alpha j}}\right),$ (13)

and the implied reversible evolution equations are

$\displaystyle\partial_{t}s_{\alpha}$
$\displaystyle=-\partial_{j}(s_{\alpha}E_{m_{\alpha j}})$ (14a)
$\displaystyle\partial_{t}m_{\alpha i}$ $\displaystyle=-\partial_{j}(m_{\alpha
i}E_{m_{\alpha j}})-s_{\alpha}\partial_{i}E_{s_{\alpha}}-m_{\alpha
j}\partial_{i}E_{m_{\alpha j}}.$ (14b)

The dependence of energy on the hydrodynamic fields can be determined by
inverting the formula for entropy of the fields (66), which was obtained by
the principle of maximization of entropy (MaxEnt) in Appendix A.

Although now we could add a dissipation potential in order to add also an
irreversible evolution of the phonons, we shall first reduce the description
to a less detailed one, parametrized only by the overall distribution function
$f=\sum_{\alpha}f_{\alpha}$, which is the subject of the following Section.

###  2.2 Overall kinetic theory of phonons

When the various polarizations of phonons are not important, we may reduce the
kinetic theory of polarized phonons to the overall kinetic theory that
describes evolution of the overall distribution function

$f=f_{T+}+f_{T-}+f_{L}.$ (15)

The Poisson bracket for the overall distribution function $f$ follows from the
Poisson bracket for state variables $f_{T+}$, $f_{T-}$, and $f_{L}$ (Equation
(8)). When we plug functionals that depend only on $f(f_{T+},f_{T-},f_{L})$ to
the Poisson bracket, we obtain bracket

$\\{A,B\\}^{(f)}=\int d\mathbf{r}\int d\mathbf{p}f\cdot\left(\frac{\partial
A_{f}}{\partial\mathbf{r}}\frac{\partial
B_{f}}{\partial\mathbf{p}}-\frac{\partial
B_{f}}{\partial\mathbf{r}}\frac{\partial A_{f}}{\partial\mathbf{p}}\right),$
(16)

which is the standard Poisson bracket for the one-particle distribution
functions (also called Boltzmann bracket), see for instance [27, 69]. The
reversible evolution equation implied by this Poisson bracket is

$(\partial_{t}f)_{\mathrm{rev}}=-\frac{\partial
f}{\partial\mathbf{r}}\cdot\frac{\partial
E_{f}}{\partial\mathbf{p}}+\frac{\partial
f}{\partial\mathbf{p}}\cdot\frac{\partial E_{f}}{\partial\mathbf{r}}.$ (17)

In order to write the evolution equation in a closed form, we need the energy
$E(f)$. In our simplified situations (restriction to the first Brillouin zone,
acoustic phonons with small $\mathbf{k}$ vectors), the energy is

$E=\int d\mathbf{r}\int d\mathbf{p}c|\mathbf{p}|f,$ (18)

where $c$ is the Debye speed [49]. Once the Debye velocity $c$ is determined,
the reversible part of the evolution equation for $f$ is fully specified. In
particular, the constant $c$ can be obtained from measurement of the second
sound speed. For instance, for a NaF crystal, Equation (48) relates the Debye
velocity with the speed of second sound, theoretical limit of which at $T=0$
is $c=2020.9\,m/s$ [37]. The heat capacity of low-temperature NaF is
$c_{V}=234N_{A}k_{B}(T/\theta_{D})^{3}$ for temperature between $0.6K$ and
$15K$; the Debye temperature is found to be $\theta_{D}=(466\pm 5)K$ [38]. The
density is approximated by $\rho\approx 2.851g/cm^{3}$ [42], and the molar
mass is $M\approx 42g/mol$.

However, the evolution of $f$ actually contains two parts, a reversible
(Hamiltonian) part, which is given by the Poisson bracket and energy, and an
irreversible part, which is given by a dissipation potential and entropy. The
entropy can be derived from the entropy of phonons with polarization, Equation
(11). By maximization of the entropy subject to the constraint (15), we obtain
the entropy

$S^{(f)}=-\frac{3k_{B}}{h^{3}}\int d\mathbf{r}\int
d\mathbf{p}\left(h^{3}\frac{f}{3}\ln\left(h^{3}\frac{f}{3}\right)-\left(1+h^{3}\frac{f}{3}\right)\ln\left(1+h^{3}\frac{f}{3}\right)\right).$
(19)

The irreversible evolution of the distribution function consists of collisions
between phonons (normal process, conserving both energy and momentum) and
collisions with impurities of the crystal (resistive process, conserves
energy, but not momentum) [63]. These two processes are described by the so
called Callaway model [9]

$(\partial_{t}f)_{\mathrm{irr}}=-\frac{1}{\tau_{R}}(f-f_{R})-\frac{1}{\tau_{N}}(f-f_{N}),$
(20)

where $f_{R}$ and $f_{N}$ are the MaxEnt estimates of the distribution
function obtained by maximization of entropy (19) subject to the constraints
given by the overall hydrodynamic momentum and energy (57), or by the
constraint that only the hydrodynamic energy is known, respectively,

$\displaystyle f_{N}$
$\displaystyle=\frac{3}{h^{3}}\frac{1}{e^{\mathbf{m}^{*}\cdot\mathbf{p}/k_{B}}e^{e^{*}c|\mathbf{p}|/k_{B}}-1}$
(21a) $\displaystyle f_{R}$
$\displaystyle=\frac{3}{h^{3}}\frac{1}{e^{e^{*}c|\mathbf{p}|/k_{B}}-1}.$ (21b)
Here, $\mathbf{m}^{*}$ and $e^{*}$ are the derivatives of the hydrodynamic
entropy (33) with respect to the momentum density and energy density,
respectively, $\mathbf{m}^{*}=\frac{\partial s^{H}}{\partial\mathbf{m}}$ and
$e^{*}=\frac{\partial s^{H}}{\partial e}$.

A feature of distribution functions $f_{N}$ and $f_{R}$ is that the Callaway
irreversible evolution conserves momentum and energy (in the normal processes)
or only the energy (in the resistive processes),

$\displaystyle\int d\mathbf{p}p_{i}(f-f_{N})=0,\int
d\mathbf{p}|\mathbf{p}|(f-f_{N})=0,\int d\mathbf{p}|\mathbf{p}|(f-f_{R})=0.$
(22)

The normal and resistive relaxation times, $\tau_{N}$ and $\tau_{R}$, have
typically values around $\tau_{R}\approx 3\cdot 10^{-6}s$ and $\tau_{N}\approx
2\cdot 10^{-7}s$ [16].

The final evolution equation is the sum of the reversible part (17) and
irreversible part (20),

$\partial_{t}f=-\frac{\partial f}{\partial\mathbf{r}}\cdot\frac{\partial
E_{f}}{\partial\mathbf{p}}+\frac{\partial
f}{\partial\mathbf{p}}\cdot\frac{\partial
E_{f}}{\partial\mathbf{r}}-\frac{1}{\tau_{R}}(f-f_{R})-\frac{1}{\tau_{N}}(f-f_{N}).$
(23)

The next Section contains a further reduction of this kinetic theory to phonon
hydrodynamics.

##  3 Overall hydrodynamics of phonons

Let us now further reduce the Boltzmann Poisson bracket, that expresses
kinematics of kinetic theory, to the Poisson bracket expressing kinematics of
the hydrodynamic fields of momentum density and entropy density,

$\displaystyle m_{i}(t,\mathbf{r})$ $\displaystyle=\int
d\mathbf{p}p_{i}f(t,\mathbf{r},\mathbf{p})$ (24a) $\displaystyle
s(t,\mathbf{r})$ $\displaystyle=\int
d\mathbf{p}\eta(f(t,\mathbf{r},\mathbf{p})).$ (24b) Here, $\eta(f)$ is the
entropy density of phonons $s^{H}$ (Equation (33)), but in general it can be
an arbitrary function of $f$.

###  3.1 Reversible evolution

Plugging arbitrary functionals $A$ and $B$ dependent only on momentum and
entropy densities, $A(\mathbf{m},s)$ and $B(\mathbf{m},s)$, into bracket (16),
we obtain a part of the Poisson brackets of fluid mechanics,

$\\{A,B\\}^{(\mathbf{m},s)}=\int
d\mathbf{r}s\left(\partial_{i}A_{s}B_{m_{i}}-\partial_{i}B_{s}A_{m_{i}}\right)+\int
d\mathbf{r}m_{i}\left(\partial_{j}A_{m_{i}}B_{m_{j}}-\partial_{j}B_{m_{i}}A_{m_{j}}\right),$
(25)

which is actually the same as the part of the Poisson bracket for fluid
mechanics of normal particles (without mass density), see [30, 69] for details
of the calculation. The reversible evolution implied by this bracket reads

$\displaystyle\partial_{t}s$ $\displaystyle=-\partial_{j}(sE_{m_{j}})$ (26a)
$\displaystyle\partial_{t}m_{i}$
$\displaystyle=-\partial_{j}(m_{i}E_{m_{j}})-s\partial_{i}E_{s}-m_{j}\partial_{i}E_{m_{i}}.$
(26b)

Similarly as in fluid mechanics, the momentum equation implied by this Poisson
bracket contains a non-linear convective term, which expresses inertial
effects of the motion of phonons and which is missing in the usual models of
hyperbolic heat conduction as the linear Guyer-Krumhansl model [35, 32].
Although the heat flux is proportional to the momentum density (Equation
(45)), the convective terms can not be written as
$\mathbf{q}\cdot\nabla\mathbf{q}$, since the derivative $E_{\mathbf{m}}$ also
depends on the field of energy density (see Equation (46)).

It is advantageous to transform the momentum density to the field

$\mathbf{w}=\mathbf{m}/s$ (27)

because then the subsequent reduction to the Fourier heat conduction becomes
simpler, see Section 3.5. Derivatives of functionals of $\mathbf{m}$ and $s$
transform as

$F_{\mathbf{m}}\rightarrow\frac{1}{s}F_{\mathbf{w}}\quad\mbox{and}\quad
F_{s}\rightarrow F_{s}-\frac{1}{s}F_{\mathbf{w}}\cdot\mathbf{w},$ (28)

and when we plug functionals dependent only on $(s,\mathbf{w})$ into the
Poisson bracket for $(\mathbf{m},s)$, we obtain a Poisson bracket for $s$ and
$\mathbf{w}$,

$\\{A,B\\}^{(Cattaneo)}=\int
d\mathbf{r}\left(\partial_{i}A_{s}B_{w_{i}}-\partial_{i}B_{s}A_{w_{i}}\right)+\int
d\mathbf{r}\frac{1}{s}(\partial_{i}w_{j}-\partial_{j}w_{i})A_{w_{i}}B_{w_{j}}.$
(29)

We refer to this Poisson bracket as the Cattaneo Poisson bracket [69] because
the resulting model simplifies to the usual Cattaneo equations when we
consider only one-dimensional processes in the linear approximation, see
Appendix B. The reversible part of the evolution equations, stemming from this
bracket, is obtained by formula (7),

$\displaystyle(\partial_{t}s)_{\mathrm{rev}}$
$\displaystyle=\\{s,E\\}^{(Cattaneo)}=-\partial_{i}E_{w_{i}}$ (30a)
$\displaystyle(\partial_{t}w_{i})_{\mathrm{rev}}$
$\displaystyle=\\{w_{i},E\\}^{(Cattaneo)}=-\partial_{i}E_{s}+\frac{1}{s}(\partial_{i}w_{j}-\partial_{j}w_{i})E_{w_{j}}.$
(30b) The latter equation can be also rewritten as
$(\partial_{t}\mathbf{w})_{\mathrm{rev}}=-\nabla
E_{s}+\frac{1}{s}E_{\mathbf{w}}\times\boldsymbol{\omega},$ (30c) where
$\boldsymbol{\omega}=\nabla\times\mathbf{w}$ is the vorticity of the field
$\mathbf{w}$. From the evolution equation for entropy, we can see that the
entropy flux is given by $E_{\mathbf{w}}$, and the field $\mathbf{w}$ can be
thus called the conjugate entropy flux.

What is the heat flux implied by Equations (30c)? If we define the heat flux
as the flux of energy, we have to first find the evolution equation for the
total energy density $e(s,\mathbf{w})$. And since we know evolution equations
for both $s$ and $\mathbf{w}$, we obtain the balance of total energy by chain
rule,

$\partial_{t}e=E_{s}\partial_{t}s+E_{w_{i}}\partial_{t}w_{i}=-\partial_{i}\left(E_{s}E_{w_{i}}\right),$
(31)

which means that the heat flux (or flux of energy) is

$\mathbf{q}=E_{s}E_{\mathbf{w}}=-T^{2}S_{\mathbf{w}}.$ (32)

In other words, the heat flux is equal to the product of temperature
($T=E_{s}$) and entropy flux ($E_{\mathbf{w}}$), which is a consequence of the
model, not an assumption. Another (equivalent) relation between the heat flux
and momentum density is derived later (Equation (45)).

What is the energy $E(s,\mathbf{w})$, that is necessary to write the evolution
equations in a closed form? It can again be derived by maximization of the
entropy for the distribution function (19) subject to the constraints given by
the knowledge of the hydrodynamic fields (Equations (24)). Appendix A contains
details of the calculation, resulting in entropy density

$s^{H}=\frac{4}{3}\left(\frac{6\sigma}{ec}\right)^{1/4}e\frac{3^{3/4}}{2^{7/4}}(3-\chi)^{1/2}(1-\chi)^{1/4},$
(33)

where $\sigma=\frac{2\pi^{5}k_{B}^{4}}{15h^{3}c^{2}}$ is the Stefan-Boltzmann
constant and
$\chi=\frac{5}{3}-\frac{4}{3}\sqrt{1-\frac{3}{4}\frac{c^{2}\mathbf{m}^{2}}{e^{2}}}$
is the Eddington factor [3]. This entropy is the same as in [63, 53]. In the
following Section, the reversible evolution is equipped with an irreversible
counterpart.

###  3.2 Zero-th order irreversible evolution

By integration of the Callaway model (Equations (20)) with respect to
momentum, we obtain the following zero-th approximation of the irreversible
evolution of the $s$ and $\mathbf{m}$ fields,

$\displaystyle(\partial_{t}\mathbf{m})_{irr}$
$\displaystyle=-\frac{1}{\tau_{R}}\mathbf{m}$ (34a)
$\displaystyle(\partial_{t}s)_{irr}$
$\displaystyle=\frac{1}{\tau_{R}E_{s}}\mathbf{m}\cdot E_{\mathbf{m}},$ (34b)
where the latter equation is obtained from the former by the requirement of
energy conservation, $(\partial_{t}e)_{irr}=0$.

These equations can be transformed into the $s,\mathbf{w}$ variables as
follows,

$\displaystyle(\partial_{t}s)_{irr}$
$\displaystyle=\frac{1}{\tau_{R}}\frac{\mathbf{w}\cdot
E_{\mathbf{w}}}{E_{s}-\frac{1}{s}\mathbf{w}\cdot E_{\mathbf{w}}},$ (35a)
$\displaystyle(\partial_{t}\mathbf{w})_{irr}$
$\displaystyle=-\frac{1}{\tau_{R}}\mathbf{w}-\frac{\mathbf{w}^{2}}{\tau_{R}s}\frac{E_{\mathbf{w}}}{E_{s}-\frac{1}{s}E_{\mathbf{w}}\cdot\mathbf{w}}.$
(35b)

The final evolution equations for the hydrodynamic fields $s$ and $\mathbf{w}$
are the sum of the reversible (Hamiltonian) evolution and irreversible
evolution,

$\displaystyle\partial_{t}s$
$\displaystyle=-\partial_{i}E_{w_{i}}+\frac{1}{\tau_{R}}\frac{\mathbf{w}\cdot
E_{\mathbf{w}}}{E_{s}-\frac{1}{s}\mathbf{w}\cdot E_{\mathbf{w}}},$ (36a)
$\displaystyle\partial_{t}w_{i}$
$\displaystyle=-\partial_{i}E_{s}+\frac{1}{s}(\partial_{i}w_{j}-\partial_{j}w_{i})E_{w_{j}}-\frac{1}{\tau_{R}}\mathbf{w}-\frac{\mathbf{w}^{2}}{\tau_{R}s}\frac{E_{\mathbf{w}}}{E_{s}-\frac{1}{s}E_{\mathbf{w}}\cdot\mathbf{w}},$
(36b) which is a nonlinear generalization of the Maxwell-Cattaneo-Vernotte
(MCV) equation, see Appendix B.

Equations (36) are a system of quasilinear first-order partial differential
equations, and they are also symmetric hyperbolic because they can be
rewritten as a system of conservation laws in the $s,\mathbf{m}$ variables,
which is suitable for the Godunov-Boillat theorem [22, 8, 77, 71]. The
dissipative part (terms that contain $\tau_{R}$) can be seen as the zero-th
order approximation in the Chapman-Enskog expansion of Equation (23). In the
following Section, we go for the first Chapman-Enskog approximation, which
brings viscous-like terms to the evolution equations.

###  3.3 First-order irreversible evolution

The first-order Chapman-Enskog approximation of the Boltzmann equation for
rarefied gases leads to the Navier-Stokes equations, which contain the viscous
terms, on top of the hyperbolic Euler equations [14, 12]. Similarly, the BGK
equation, which is actually a form of the here used Callaway model, gives the
Navier-Stokes equation as the first Chapman-Enskog approximation [7, 78]. Our
goal is now to obtain the viscous-like terms that appear in the momentum
equation in the first Chapman-Enskog approximation of the Callaway model also
in the phonon hydrodynamics, similarly as in [33] but in our state variables.

The Chapman-Enskog expansion starts with the projection of the equation for
$f$ (Equation (23)) to the momentum density $\mathbf{m}$,

$\partial_{t}m_{i}=-\frac{\partial}{\partial r^{j}}\left(\int
d\mathbf{p}p_{i}\frac{cp_{j}}{|\mathbf{p}|}f\right)-\frac{1}{\tau_{R}}m_{i}.$
(37a) For simplicity, we now disregard the upper/lower positions of the
indexes, which means that the metric tensor is implicitly present in the
calculations. In order to provide an approximation for the unknown integral on
the right hand side, $Q_{ij}=\int d\mathbf{p}p_{i}p_{j}/|\mathbf{p}|f$, the
evolution of that integral is obtained again from the evolution equation for
$f$, $\partial_{t}Q_{ij}=-\frac{\partial}{\partial r^{k}}\left(\int
d\mathbf{p}c\frac{p_{i}p_{j}p_{k}}{|\mathbf{p}|}f\right)+\frac{1}{\tau_{N}}\int
d\mathbf{p}\frac{p_{i}p_{j}}{|\mathbf{p}|}(f_{N}-f)+\frac{1}{\tau_{R}}\int
d\mathbf{p}\frac{p_{i}p_{j}}{|\mathbf{p}|}(f_{R}-f).$ (37b)

Appendix E contains approximate values of the integrals on the right hand side
of the latter equation, with which the evolution equation for $Q_{ij}$ becomes

$\partial_{t}Q_{ij}\approx\frac{\partial}{\partial
r^{k}}\left(\frac{16\pi^{5}k_{B}^{4}}{75h^{3}(e^{*})^{5}c^{4}}(m^{*}_{i}\delta^{k}_{j}+(m^{*})^{k}\delta_{ji}+m^{*}_{j}\delta_{i}^{k})\right)-\left(\frac{1}{\tau_{R}}+\frac{1}{\tau_{N}}\right)\left(Q_{ij}-\frac{4\pi^{5}k_{B}^{4}}{15h^{3}(e^{*})^{4}c^{4}}\delta_{ij}\right).$
(38)

In order to solve approximately (38), we first cast the equation into a
dimensionless form. The underlying evolution equation for the distribution
function (23) has the dimensionless form

$\frac{\partial\hat{f}}{\partial\hat{t}}=-\mathrm{St}\frac{\hat{\mathbf{p}}}{|\hat{\mathbf{p}}|}\cdot\hat{\nabla}\hat{f}+\frac{1}{\mathrm{De}_{N}}(\hat{f}_{N}-\hat{f})+\frac{1}{\mathrm{De}_{R}}(\hat{f}_{R}-\hat{f}),$
(39)

where the quantities with hats are dimensionless,
$\mathrm{St}=\frac{c\bar{t}}{\bar{r}}$ is the Strouhal number, $\bar{t}$ is a
typical time value, $\bar{r}$ is a typical spatial dimension, and
$\mathrm{De}_{N}=\tau_{N}/\bar{t}$ and $\mathrm{De}_{R}=\tau_{R}/\bar{t}$ are
the Deborah numbers corresponding to the two relaxation times. Because the
phonon gas first reaches the quasi-equilibrium manifold (defined by the
momentum density and energy density) before it reaches the equilibrium
(defined only by the energy), the normal relaxation time is assumed to be much
smaller than the resistive one, $\tau_{N}\ll\tau_{R}$, which means that
$\mathrm{De}_{N}\ll\mathrm{De}_{R}$. Therefore, the second moment $Q_{ij}$ can
be expanded in powers of the normal Deborah number,

$Q_{ij}=Q^{(0)}_{ij}+\mathrm{De}_{N}Q^{(1)}_{ij}+\mathcal{O}(\mathrm{De}_{N}^{2}).$
(40)

Plugging this expansion into the equation for $Q_{ij}$ (Equation (38)), the
$\mathrm{De}_{N}^{-1}$ order gives

$\mathbf{Q}^{(0)}=\frac{4\pi^{5}k_{B}^{4}}{15h^{3}(e^{*})^{4}c^{4}},$ (41a)
and the order $\mathrm{De}_{N}^{0}$ gives
$\partial_{t}Q^{(0)}_{ij}=\partial_{k}\left(\frac{16\pi^{5}k_{B}^{4}}{75h^{3}(e^{*})^{5}c^{4}}(m^{*}_{i}\delta^{k}_{j}+(m^{*})^{k}\delta_{ji}+m^{*}_{j}\delta_{i}^{k})\right)-{\bar{t}}^{-1}\left(1+\frac{\mathrm{De}^{N}}{\mathrm{De}^{R}}\right)Q^{(1)}_{ij}.$
(41b)

These two orders then combine to

$Q_{ij}=\frac{4\pi^{5}k_{b}^{4}}{15h^{3}(e^{*})^{4}c^{4}}+\tau\partial_{k}\left(\frac{16\pi^{5}k_{B}^{4}}{75h^{3}(e^{*})^{5}c^{4}}(m^{*}_{i}\delta^{k}_{j}+(m^{*})^{k}\delta_{ji}+m^{*}_{j}\delta_{i}^{k})\right)-\tau\frac{4\pi^{5}k_{b}^{4}}{15h^{3}c^{4}}\partial_{t}(e^{*})^{-4}\delta_{ij}$
(42)

where $\tau^{-1}=\tau_{R}^{-1}+\tau_{N}^{-1}$ is the effective relaxation
time.

The approximate $\mathbf{Q}$ from Equation (42) can be now plugged into the
momentum equation (37a). However, since we have already obtained the
reversible part of the latter equation precisely, using the Poisson reduction,
we use the approximate $\mathbf{Q}$ only in the irreversible part of the
equation,

$(\partial_{t}m_{i})_{\mathrm{irr}}=-\frac{1}{\tau_{R}}m_{i}-\partial_{j}\left(\tau\partial_{k}\left(\frac{16\pi^{5}k_{B}^{4}}{75h^{3}(e^{*})^{5}c^{3}}(m^{*}_{i}\delta^{k}_{j}+(m^{*})^{k}\delta_{ji}+m^{*}_{j}\delta_{i}^{k})\right)-\tau\frac{4\pi^{5}k_{b}^{4}}{15h^{3}c^{3}}\partial_{t}(e^{*})^{-4}\delta_{ij}\right).$
(43) The reversible part of the momentum equation is contained in Equation
(26b), and the overall evolution of the momentum density is the sum of th
reversible and irreversible parts, $\displaystyle\partial_{t}m_{i}=$
$\displaystyle\partial_{j}\left(m_{i}\frac{m^{*j}}{e^{*}}\right)-s(\mathbf{m},e)\partial_{i}\left(\frac{1}{e^{*}}\right)+m_{j}\partial_{i}\frac{m^{*j}}{e^{*}}$
(44a)
$\displaystyle-\frac{1}{\tau_{R}}m_{i}-\partial_{j}\left(\tau\partial_{k}\left(\frac{16\pi^{5}k_{B}^{4}}{75h^{3}(e^{*})^{5}c^{3}}(m^{*}_{i}\delta^{k}_{j}+(m^{*})^{k}\delta_{ji}+m^{*}_{j}\delta_{i}^{k})\right)-\tau\frac{4\pi^{5}k_{b}^{4}}{15h^{3}c^{3}}\partial_{t}(e^{*})^{-4}\delta_{ij}\right),$
where derivatives of energy were replaced with the corresponding derivatives
of entropy, $E_{\mathbf{m}}=-S_{\mathbf{m}}/S_{e}$, $E_{s}=(S_{e})^{-1}$. The
right hand side of the momentum equation contains the conjugate variables
$\mathbf{m}^{*}=\frac{\partial S^{H}}{\partial\mathbf{m}}$ and
$e^{*}=\frac{\partial S^{H}}{\partial e}$, but we know the entropy (Equation
(33)), which makes the right hand side explicit in $e$ and $\mathbf{m}$.

Finally, the momentum equation needs to be supplemented with the energy
equation,

$\partial_{t}e=\int
d\mathbf{p}c|\mathbf{p}|\partial_{t}f=-\partial_{i}(c^{2}m_{i}),$ (44b)

which also shows that the heat flux is proportional to the momentum density,

$\mathbf{q}=c^{2}\mathbf{m}.$ (45)

Equations (44) govern evolution of momentum density and energy density of
phonons in the first Chapman-Enskog approximation.

In order to better see the viscous-like character of the irreversible terms in
the momentum equation, let us use the approximate values of the Lagrange
conjugate variables (64), which give that
$\mathbf{m}^{*}/e^{*}\approx-\frac{3c^{2}}{4e}\mathbf{m}$ and
$\frac{16\pi^{5}k_{B}^{4}}{75h^{3}(e^{*})^{4}c^{3}}=\frac{8\sigma}{5c(e^{*})^{4}}\approx\frac{4}{15}e$.
Moreover, when $\mathbf{m}^{2}\approx 0$, it holds that $2\sigma
c^{-1}(e^{*})^{-4}\approx e/3$, and the momentum equation then simplifies to

$\displaystyle\partial_{t}m_{i}=$
$\displaystyle-\partial_{i}p-\partial_{j}\left(\frac{3c^{2}}{4e}m_{i}m_{j}\right)$
(46)
$\displaystyle-\frac{1}{\tau_{R}}m_{i}+\frac{\tau}{5}\partial_{j}\left(\partial_{j}m_{i}+\partial_{i}m_{j}-\frac{2}{3}\nabla\cdot\mathbf{m}\delta_{ij}\right)$

with pressure

$p=-e+\frac{4}{3}e(2\sqrt{y}-1)+\frac{3c^{2}}{2e}\frac{\mathbf{m}^{2}}{1+\sqrt{y}}=\frac{e}{3}-\frac{c^{2}\mathbf{m}^{2}}{4e}+\mathcal{O}(\mathbf{m}^{4})\approx\frac{e}{3}-\frac{c^{2}\mathbf{m}^{2}}{4e}.$
(47)

Equations (44b) and (46) are a partial linearization of the full equations for
$\mathbf{m}$ and $e$. However, they still represent a more precise version of
the Guyer-Krumhansl equation, since they still contain the convective terms in
the phonon momentum density. The precision is enhanced due to the Poisson
reduction, which makes it possible to obtain the reversible part of the
equation without linearizations. The viscous terms in the momentum equation
are similar to those found in [33], except for the terms with
$\nabla\cdot\mathbf{m}$ (here, the viscous part of the Cauchy stress tensor is
traceless). Equation (46) is similar to the evolution equation for heat flux
in [13], but it still differs by the presence of the convective term.

Equations (44b) and (46) in one dimension can be linearized around a
stationary homogeneous state, $e\approx e_{0}+\delta e$ and
$\mathbf{m}=\delta\mathbf{m}$, which leads to a wave equation (neglecting the
dissipative terms),

$\frac{\partial^{2}e}{\partial
t^{2}}=\frac{c^{2}}{3}\frac{\partial^{2}e}{\partial x^{2}}.$ (48)

Therefore, the speed of the second sound can be related with the Debye
velocity through $c_{\text{(2nd)}}=c/\sqrt{3}$, as reported in [32].

We would like to solve these equations numerically in the future, for which we
shall need boundary conditions. Although boundary conditions are out of scope
of the present paper, let us briefly mention at least some results from the
literature. Apart from no heat flux across an insulated wall
$\mathbf{q}\cdot\mathbf{n}=0$, one could consider zero heat flux also along
the insulating wall (non-slip condition, $\mathbf{q}=0$), or a non-zero value
(slip condition) analogous to the no-slip condition of velocity for rarefied
gases [2, 80]. These conditions are especially relevant in narrow channels
[79, 55]. Boundary conditions for the heat flux have also been explored within
IVT [83], and a derivation of boundary conditions from the BGK equation for
phonons can be found in [33].

Let us now discuss some features of Equations (36) (obtained by the Poisson
reduction and the zero-th Chapman-Enskog approximation).

###  3.4 Galilean invariance, inertia, and thermal vortices

Equations (36) have several features that are not present in the usual models
for hyperbolic heat conduction, Galilean invariance, inertial effects, and
effects of thermal vorticity. First, their reversible part (the Hamiltonian
part) is Galilean invariant, as shown in Appendix C. Galilean invariance of
the irreversible part of the entropy equation is achieved if we count with the
possible motion of the crystal (with velocity $\mathbf{V}$) by taking energy
$e=e(s,\mathbf{w})+s\mathbf{V}\cdot\mathbf{w}$.

Second, the evolution equation for $\mathbf{w}$ contains inertial terms, which
are proportional to $E_{\mathbf{w}}$. In other words, if we abruptly change
the temperature gradient, it takes some time to change the direction of
$\mathbf{w}$.

The zero-th Chapman-Enskog approximate evolution of $\mathbf{w}$ (36) also
expresses how the thermal vorticity
$\boldsymbol{\omega}=\nabla\times\mathbf{w}$ affects the evolution of
$\mathbf{w}$ itself. In the stationary state, the evolution equation for
$\mathbf{w}$ yields an analogy of the Crocco theorem [4]

$\nabla
T=\frac{1}{s}E_{\mathbf{w}}\times\boldsymbol{\omega}-\frac{1}{\tau_{R}}\mathbf{w},$
(49)

here called the Fourier-Crocco equation.222Here, we neglect the second
irreversible term in the equation for $\mathbf{w}$ in (36) because it is of
higher order in $\mathbf{w}$. Figure 2 shows this equation in the case of two
dimensions. In particular, for the limit of weak dissipation,
$\tau_{R}\rightarrow\infty$, the heat flux is perpendicular to the temperature
gradient. On the other hand, in the limit of strong dissipation,
$\tau_{R}\rightarrow 0$, the heat flux aligns with the temperature gradient,
as expected from the Fourier theory of heat conduction.

Figure 2:  The Fourier-Crocco equation for in the presence of thermal
vorticity $\boldsymbol{\omega}$. The conjugate entropy flux $\mathbf{w}$,
which is proportional to the heat flux, is not parallel with the gradient of
temperature.

How can be thermal vorticity $\boldsymbol{\omega}$ created? Consider for
instance a body with a gradient of thermal conductivity in the $x$-direction.
Thermal conductivity is proportional to $\tau_{R}^{-1}$ (see the following
Section). When we cast the curl operator on the evolution equation for
$\mathbf{w}$ (36b), we obtain that

$\partial_{t}\boldsymbol{\omega}|_{\boldsymbol{\omega}=0}=-\nabla\left(\frac{1}{\tau_{R}}\right)\times\mathbf{w},$
(50)

even in the case where $\boldsymbol{\omega}=0$. Note that now we have to allow
for $\tau_{R}$ not to be a constant. This means that thermal vorticity
$\boldsymbol{\omega}$ is created when the gradient of $\tau^{-1}_{R}$ is not
aligned with $\mathbf{w}$. For instance, $\mathbf{w}$ is created in the
$y$-direction (by applying a temperature gradient to the body), some amount of
thermal vorticity is created. We expect this behavior in thermal diodes [47,
90]. The thermal vorticity can be also created at the walls (by boundary
conditions) and at sharp corners as in fluid mechanics [4].

What effects can be expected in the presence of thermal vorticity? Since the
Fourier-Crocco equation (49) is similar to the classical Crocco equation from
fluid mechanics, we can find inspiration in aerodynamics. Figure 3(a) shows
the flow of phonons around an insulated cylinder. Since the cylinder is
insulated, there is no heat flux across its boundary, so that
$\mathbf{w}\cdot\mathbf{n}=0$ at the boundary. In the stationary state, when
the Fourier-Crocco equation applies, we have, moreover, that
$\nabla\cdot\mathbf{m}=0$, which brings us close to dynamics of incompressible
fluids. From the Fourier-Crocco equation, we can expect that the temperature
behind the cylinder is lower than the temperature in front of the cylinder.
The analogy in classical fluid mechanics would be that the hydrodynamic drag
force (or lower pressure behind the cylinder). In other words, an insulated
obstacle in a flow of phonons should create a region of lower temperature
behind the obstacle. This is, however, still to be verified experimentally.

(a) Flow of phonons past an insulated cylinder create a colder spot behind the
cylinder.

(b) Flow past an airfoil-like insulated body creates a temperature gradient
(hotter below).

Similarly, Figure 3(b) shows that flow of phonons along an insulated airfoil
should create a region or lower temperature above the airfoil and a region of
higher temperature below the airfoil. No such effects could be obtained
without the non-linear vorticity-dependent terms, that are usually neglected
in standard theories of hyperbolic heat conduction [35].

###  3.5 Reduction to Fourier heat transport

When we asymptotically expand the conjugate entropy flux
$\mathbf{w}=\mathbf{w}^{(0)}+\mathrm{De}_{R}\mathbf{w}^{(1)}+\mathcal{O}(\mathrm{De}_{R}^{2})$,
then from the evolution equation for $\mathbf{w}$ we obtain that
$\mathbf{w}^{(0)}=0$, $\mathbf{w}^{(1)}=-\nabla T$, and
$\mathbf{w}\approx-\tau_{R}\nabla T$. The evolution equation for entropy then
becomes

$\partial_{t}s=\nabla\cdot\left(\tau_{R}\kappa\nabla
T\right)+\frac{\tau_{R}\kappa}{T}(\nabla T)^{2},$ (51)

where $\kappa=\frac{3c^{2}s^{2}}{4e}$ is the coefficient in approximated
energy $e\approx e^{(0)}(s)+\frac{1}{2}\kappa\mathbf{w}^{2}$. In terms of
energy density, this equation reads

$\partial_{t}e=-\nabla\cdot(-\tau_{R}\kappa T\nabla T),$ (52)

which means that the heat flux, or flux of energy, is
$\mathbf{q}=-\tau_{R}T\kappa\nabla T$, that represents Fourier law with
thermal conductivity $\lambda=\tau_{R}T\kappa$. Fourier heat conduction can be
thus seen as the first-order approximation of the asymptotic expansion of
Equations (36) in the limit of strong dissipation.

###  3.6 A hypothetical dependence of energy on vorticity

In this Section, we add to the energy a hypothetical dependence on the thermal
vorticity $\boldsymbol{\omega}$, and replace the energy density derived by
maximization of the entropy of bosons by

$\bar{E}=\int d\mathbf{r}\bar{e}(s,\mathbf{w},\boldsymbol{\omega}),$ (53)

where is $\bar{e}(s,\mathbf{w},\boldsymbol{\omega})$ is a function of entropy,
the conjugate entropy flux, and thermal vorticity. In case of
$\boldsymbol{\omega}=0$, $\bar{e}$ should turn to the standard formula (33),
but we leave it undetermined because it is not clear whether such
(hypothetical) dependence is possible from the physical point of view. Phonons
express collective vibrations of atoms of a lattice, which means that
$\boldsymbol{\omega}$ encodes a sort of angular momentum of the vibrations. A
similar energy density can be seen for instance in superfluids, where
dependence of energy on vorticity of the superfluid velocity is due to the
presence of topological defects (quantum vortices) [39, 61, 89, 48]. We are
not aware, however, of any experimental observation of similar vortices in the
case of phonons.

Anyway, let us close the paper with exposition of the consequences of the
hypothetical dependence of energy density on thermal vorticity. The functional
derivative of energy with respect to the $\mathbf{w}$ field turns to

$\frac{\delta\bar{E}}{\delta\mathbf{w}}=\frac{\partial
e}{\partial\mathbf{w}}+\nabla\times\frac{\partial\bar{e}}{\partial\boldsymbol{\omega}}.$
(54)

The evolution equation for $\mathbf{w}$ (Equation (36b)) and the equation for
curl of $\mathbf{w}$ become

$\displaystyle\partial_{t}\mathbf{w}$ $\displaystyle=-\nabla
T+\frac{1}{s}\frac{\partial\bar{e}}{\partial\mathbf{w}}\times\boldsymbol{\omega}+\frac{1}{s}\left(\nabla\frac{\partial\bar{e}}{\partial\mathbf{w}}\cdot\boldsymbol{\omega}-\omega\cdot\nabla\frac{\partial\bar{e}}{\partial\mathbf{w}}\right)-\frac{1}{\tau_{R}}\mathbf{w}$
(55a) $\displaystyle\partial_{t}\boldsymbol{\omega}$
$\displaystyle=\nabla\cdot\left(\frac{1}{s}\left(\frac{\partial\bar{e}}{\partial\mathbf{w}}\otimes\boldsymbol{\omega}-\boldsymbol{\omega}\otimes\frac{\partial\bar{e}}{\partial\mathbf{w}}\right)\right)+\nabla\times\left(\frac{1}{s}\left(\nabla\frac{\partial\bar{e}}{\partial\mathbf{w}}\cdot\boldsymbol{\omega}-\omega\cdot\nabla\frac{\partial\bar{e}}{\partial\mathbf{w}}\right)\right)+\mathbf{w}\times\nabla\left(\frac{1}{\tau_{R}}\right)-\frac{1}{\tau_{R}}\boldsymbol{\omega}.$
(55b)

The equation for thermal vorticity has thus two fluxes (terms under the
$\nabla$ operator), a source term (with gradient of $\tau_{R}$), and a sink
term (proportional to $\boldsymbol{\omega}$). In particular, if
$\boldsymbol{\omega}$ is zero initially and if there is no gradient in
$\tau^{-1}_{R}$, then $\boldsymbol{\omega}$ remains zero. On the other hand,
if there is a gradient in $\tau_{R}$, as in case of a thermal diode [47, 90],
then some amount of $\boldsymbol{\omega}$ is being created until a balance
with the sink term is reached.

Heat flow also exhibits hydrodynamics behaviour in superfluid helium-4; for
instance, in the laminar regime in counterflow along a cylindrical tube or
along tubes or rectangular or square cross section, the behaviour of the heat
flow is analogous to that of viscous fluids [59, 60, 61]. In the two-fluid
model [86, 1, 48, 39] this may be easily interpreted, as the heat flow is
transported by the viscous component of the fluid [85]. Exploration of a
possible von Kármán vortex street when an insulating cylinder perpendicular to
the heat flux is introduced in a parallelepiped channel could be of interest.
However, in superfluid helium the circulation of the vortices is quantized
(due to the requirement that the $\nabla\times\mathbf{w}=0$ almost
everywhere), which makes it different from the phonon hydrodynamics explored
in the present paper.

##  4 Conclusion

In this paper, we begin with the evolution equation of distribution functions
of phonons with two transversal and one longitudinal polarizations (actually
kinetic theory of mixtures for $f_{T+}$, $f_{T-}$, and $f_{L}$). The energy is
simplified to the energy of acoustic phonons in the first Brillouin zone, and
the entropy is obtained as the entropy of an ideal bosonic gas in kinetic
theory.

This description is then reduced to hydrodynamics with momentum densities and
entropy densities of the particular polarizations ($\mathbf{m}_{\alpha}$ and
$s_{\alpha}$ for $\alpha\in\\{T+,T-,L\\}$). The underlying Poisson bracket is
obtained by reduction from the kinetic Poisson bracket and the entropy is
derived from the entropy in kinetic theory by means of MaxEnt. The kinetic
theory of polarized phonons is also reduced to the overall kinetic theory of
phonons with distribution function $f=\sum_{\alpha}f_{\alpha}$ (disregarding
polarizations), where the Poisson bracket is again obtained by reduction from
the kinetic Poisson bracket, and entropy is obtained again by MaxEnt.

Finally, these levels of description are further reduced to phonon
hydrodynamics described by state variables
$\mathbf{m}=\sum_{\alpha}\mathbf{m}_{\alpha}$ and $s=\sum_{\alpha}s_{\alpha}$
(momentum density and entropy density). The Poisson bracket and entropy are
again derived by reduction from their more detailed counterparts. Apart from
the zero-th Chapman-Enskog approximation, which leads to hyperbolic equations,
we also carry out the first Chapman-Enskog approximation, which adds higher-
order (Laplacian) terms to the evolution equation for phonon momentum density.
It turns out more practical to reformulate the hydrodynamic equations on this
level of description to variables $\mathbf{w}=\mathbf{m}/s$ and $s$, since
then we can better see the relation with standard theories of heat transport.
The obtained evolution equations are Galilean invariant, unlike usual theories
of heat transport, and they depend on the thermal vorticity,
$\boldsymbol{\omega}=\nabla\times\mathbf{w}$.

In particular, even in the stationary case, the heat flux turns out to be not
aligned with the gradient of temperature, becoming proportional to it only in
the limit of strong dissipation. In the one-dimensional case and after
linearization, the theory simplifies to the Maxwell-Cattaneo-Vernotte
equation. The thermal vorticity is created if gradients in thermal
conductivity are present that are not parallel with the heat flux. In the
presence of thermal vorticity, we can expect various phenomena similar to
aerodynamics, like a Fourier-Crocco equation (49), drop of temperature
(instead of pressure) behind and obstacle, or temperature gradient
perpendicular to the flow around an airfoil (instead of aerodynamic lift).

In future, we would like to investigate the resulting equations numerically in
order to quantitatively observe the here proposed phenomena (for instance the
aerodynamic analogies). We would also like to explore instabilities of the
equations. It would also be interesting to apply the present formalism to the
examples illustrated by [91].

## Acknowledgment

MP is grateful to Alberto Montefusco for asking the questions how exactly the
dynamics of the $s$ and $\mathbf{w}$ fields correspond to the classical
Cattaneo equation. MP is also grateful to Miroslav Grmela for the countless
discussions of the GENERIC framework, to Václav Klika for discussing the
Onsager-Casimir reciprocal relations, and to Ilya Peshkov for discussing
hyperbolic PDE’s. MP was supported by project No. UNCE/SCI/023 of Charles
University Research program and by the Czech Science Foundation (project no.
20-22092S).

##  Appendix A Maximum entropy principle for energy and entropy

Entropy of phonons, which are bosons, is expressed in terms of their
distribution function [51]

$S^{\text{kinetic}}(f)=-\frac{k_{B}}{h^{3}}\int d\mathbf{r}\int
d\mathbf{p}\left(h^{3}f\ln(h^{3}f)-(1+h^{3}f)\ln(1+h^{3}f)\right),$ (56)

where $h$ is the Planck constant. The goal is now to maximize the entropy with
the constraints given by the knowledge the energy density and momentum density
fields. This procedure is equivalent to that in [63], with the difference that
here we obtain also the entropy for an intermediate level of description where
the various polarizations of phonons have their own hydrodynamic entropies.

Indeed, because there are two transverse and one longitudinal polarizations of
phonons, we actually have three distribution functions $f_{T+}$, $f_{T-}$, and
$f_{L}$, to which correspond three entropies, given by formula (56) for each
of the distribution functions. The total entropy is then the sum of the
entropies for each polarization, Equation (11). In order to obtain analytical
results, we neglect phonon-phonon interactions, as well as phonon-crystal
interactions. We shall be constrained to the first Brillouin zone, where all
the polarizations have the same formula for energy. Moreover, we stay
restricted only to the case of acoustic phonons. This will give us the
hydrodynamic entropy of the mixture of the three types of acoustic phonons.
Eventually, we shall apply MaxEnt once more to obtain the overall hydrodynamic
entropy of all acoustic phonons, which turns out to be the same as is [63].

The projection to the fields of momentum density and energy density is

$\displaystyle\mathbf{m}(\mathbf{r})$ $\displaystyle=\int
d\mathbf{p}\,\mathbf{p}f$ (57a) $\displaystyle e(\mathbf{r})$
$\displaystyle=\int d\mathbf{p}c|\mathbf{p}|f,$ (57b)

where $c$ is the speed of sound and $|\mathbf{p}|$ is the Euclidean norm of
the momentum, $\mathbf{p}=\hbar\mathbf{k}$. By maximization of entropy (11)
subject to the constraints given by fields $e$ and $\mathbf{m}$,

$\frac{\delta}{\delta f}\left(-S^{\text{kinetic}}+\int
d\mathbf{r}\mathbf{m}^{*}(\mathbf{r})\cdot\int d\mathbf{p}\mathbf{p}f+\int
d\mathbf{r}e^{*}(\mathbf{r})\int d\mathbf{p}c|\mathbf{p}|f\right)=0,$ (58)

we obtain the MaxEnt estimate of the distribution function

$f_{\text{MaxEnt}}(\mathbf{m},e)=\frac{1}{h^{3}}\frac{1}{e^{\mathbf{m}^{*}\cdot
p/k_{B}}e^{e^{*}c|\mathbf{p}|/k_{B}}-1},$ (59)

where $\mathbf{m}^{*}$ and $e^{*}$ are the corresponding Lagrange multipliers.
The hydrodynamic formula for entropy of phonons is then given by
$S^{\text{kinetic}}(f_{\text{MaxEnt}})$, but in order to express it in terms
of state variables $\mathbf{m}$ and $e$, we have to first express the Lagrange
multipliers in terms of the state variables.

By integration of equations (57) with the MaxEnt estimate of the distribution
function (59), we obtain that

$\displaystyle m_{i}$
$\displaystyle=-\frac{16}{45}\frac{k_{B}^{4}\pi^{5}}{h^{3}(ce^{*})^{5}}\frac{m^{*}_{i}}{(1-\gamma^{2})^{3}}$
(60a) $\displaystyle e$
$\displaystyle=\frac{4}{45}\frac{\pi^{5}k_{B}^{4}}{(e^{*})^{4}c^{3}h^{3}}\frac{\gamma^{2}+3}{(1-\gamma^{2})^{3}},$
(60b)

where $\gamma=\frac{|\mathbf{m}^{*}|}{e^{*}c}$. Note that is is not trivial to
carry out the integration and several tricks are necessary. Apart from
standard techniques like swapping integrals and series or transformation to
the spherical coordinates, it is useful to rotate the coordinates so that
$\mathbf{m}^{*}$ becomes $(0,0,|\mathbf{m}^{*}|)$. Equations (60) have to be
solved so that the multipliers, $\mathbf{m}^{*}$ and $e^{*}$, are expressed in
terms of the state variables, $\mathbf{m}$ and $e$. In particular, the
equation imply that

$\frac{\beta^{2}}{16}(\gamma^{4}+6\gamma^{2}+9)=\gamma^{2},$ (61)

where $\beta=\frac{|\mathbf{m}|c}{e}$, which gives that

$\gamma^{2}=3\frac{1-\sqrt{y}}{1+\sqrt{y}}$ (62)

with $y=1-\frac{3}{4}\beta^{2}$. Note that this quadratic equation has two
roots, one of which is has no physical meaning because it would lead to a
singularity at $\beta=0$ (equilibrium, where $\mathbf{m}=0$). Since $\beta$ is
a function of the state variables, the Lagrange multipliers are now expressed
in terms of the state variables, $\mathbf{m}$ and $e$, as

$\displaystyle(e^{*})^{4}$
$\displaystyle=\frac{3\sigma}{2ec}\frac{(1+\sqrt{y})^{2}}{(2\sqrt{y}-1)^{3}}$
(63a) $\displaystyle\mathbf{m}^{*}$
$\displaystyle=-3\frac{c^{3}(3\sigma)^{1/4}}{(2ec)^{5/4}}\frac{1}{\sqrt{1+\sqrt{y}}}\frac{1}{(2\sqrt{y}-1)^{3/4}}\mathbf{m}.$
(63b)

Near the equilibrium, where
$y=1-\frac{3}{4}\frac{\mathbf{m}^{2}c^{2}}{e^{2}}\approx 1$, the multipliers
are approximated by

$e^{*}\approx\left(\frac{6\sigma}{ec}\right)^{1/4}\qquad\text{ and
}\qquad\mathbf{m}^{*}\approx-3\frac{c^{3}(3\sigma)^{1/4}}{\sqrt{2}(2ec)^{5/4}}\mathbf{m}.$
(64)

The hydrodynamic entropy is then obtained by plugging
$f_{\text{MaxEnt}}(\mathbf{m},e)$ to the kinetic formula for entropy
$S^{\text{kinetic}}$, which after some calculations results in

$S^{\text{kinetic},\alpha}(f_{\text{MaxEnt}})(\mathbf{m}_{\alpha},e_{\alpha})=\int
d\mathbf{r}\frac{4}{3}\left(\frac{2\sigma}{e_{\alpha}c}\right)^{1/4}e_{\alpha}\frac{3^{3/4}}{2^{7/4}}(3-\chi_{\alpha})^{1/2}(1-\chi_{\alpha})^{1/4},$
(65)

$\alpha\in\\{T+,T-,L\\}$, $\sigma=\frac{2\pi^{5}k_{B}^{4}}{15h^{3}c^{2}}$ is
the Stefan-Boltzmann constant, and
$\chi_{\alpha}=\frac{5}{3}-\frac{4}{3}\sqrt{1-\frac{3}{4}\frac{c^{2}\mathbf{m}_{\alpha}^{2}}{e_{\alpha}^{2}}}$
is the Eddington factor. Because we actually have three polarizations, $T+$,
$T-$, and $L$, we also have three entropies $S^{\text{kinetic},\alpha}$, while
the total entropy is their sum,

$S^{\text{HP}}=\sum_{\alpha}S^{\text{kinetic},\alpha}(f_{\text{MaxEnt}})(\mathbf{m}_{\alpha},e_{\alpha}),$
(66)

referred to as the hydrodynamic (H) entropy of polarized (P) phonons. In other
words, we have obtained the entropy for the level of description with state
variables $\mathbf{m}_{T+}$, $\mathbf{m}_{T-}$, $\mathbf{m}_{L}$, $e_{T+}$,
$e_{T-}$, and $e_{L}$.

The level of description where the various polarizations of acoustic phonons
have their own momentum and energy allows for description of processes like
change of polarization of the phonons. However, when such processes are not
important, one may opt out for a simpler level of description, where only the
overall momentum and energy,

$\displaystyle\mathbf{m}$
$\displaystyle=\mathbf{m}_{T+}+\mathbf{m}_{T-}+\mathbf{m}_{L}$ (67a)
$\displaystyle e$ $\displaystyle=e_{T+}+e_{T-}+e_{L},$ (67b)

are known. What is the entropy on such reduced level of description? The
answer is given again by maximization of the total entropy (66) subject to the
constraints given by projection (67), that is by

$\frac{\delta
S^{\text{kinetic},\alpha}}{\delta\mathbf{m}_{\alpha}}\Big{|}_{\mathbf{m}^{\text{MaxEnt}}_{\alpha},e^{\text{MaxEnt}}_{\alpha}}=\mathbf{m}^{*}\quad\text{and}\quad\frac{\delta
S^{\text{kinetic},\alpha}}{\delta
e_{\alpha}}\Big{|}_{\mathbf{m}^{\text{MaxEnt}}_{\alpha},e^{\text{MaxEnt}}_{\alpha}}=e^{*},\quad\forall\alpha.$
(68)

Since all the formulas for $S^{\text{kinetic},\alpha}$ share the same form,
their solutions are the same and the MaxEnt estimates of $\mathbf{m}_{\alpha}$
and $e_{\alpha}$ are the same for all polarizations $\alpha$. Finally, the
overall hydrodynamic entropy is obtained by plugging the MaxEnt estimates to
the formula for the total entropy (66),

$S^{\text{H}}(\mathbf{m},e)=S^{\text{HP}}(\mathbf{m}^{\text{MaxEnt}}_{\alpha},e^{\text{MaxEnt}}_{\alpha})=\int
d\mathbf{r}\frac{4}{3}\left(\frac{6\sigma}{ec}\right)^{1/4}e\frac{3^{3/4}}{2^{7/4}}(3-\chi)^{1/2}(1-\chi)^{1/4},$
(69)

where
$\chi=\frac{5}{3}-\frac{4}{3}\sqrt{1-\frac{3}{4}\frac{c^{2}\mathbf{m}^{2}}{e^{2}}}$
is the overall Eddington factor.

In the equilibrium, we obtain the equilibrium volumetric entropy density

$s^{\text{E}}=\frac{4}{3}\left(\frac{6\sigma}{c}\right)^{1/4}e^{3/4},$ (70)

which means that $e=\frac{6\sigma}{c}T^{4}$ with temperature
$T=\left(\frac{\partial s^{\text{E}}}{\partial e}\right)^{-1}$.

##  Appendix B Simplification to the Cattaneo equation

How does the dynamics for $s$ and $w$ (36) correspond to the classical
Cattaneo equation [11]? Because Cattaneo equation is a linear equation for
temperature in one dimension. Therefore, we shall be restricted to the one-
dimensional case, $\mathbf{w}=(w,0,0)$, we shall transform the
$(s,\mathbf{w})$ variables to $(T,\mathbf{w})$, and finally we will neglect
all nonlinear terms.

The derivative of energy $e$ transform to the derivatives of the free energy
$f=e-Ts$ as

$\left(\frac{\partial e}{\partial
s}\right)_{\mathbf{w}}=T,\quad\left(\frac{\partial
e}{\partial\mathbf{w}}\right)_{s}=\left(\frac{\partial
f}{\partial\mathbf{w}}\right)_{T}.$ (71)

Derivatives of the free energy (in the $(T,w)$ variables) will be denoted by
subscripts. Using the chain rule, Equations (36) become

$\displaystyle\left(\frac{\partial s}{\partial
T}\right)_{\mathbf{w}}\partial_{t}T+\left(\frac{\partial s}{\partial
w}\right)_{T}\cdot\partial_{t}w=-\partial_{x}f_{w}+\frac{1}{\tau_{R}T}f_{w}^{2}$
(72a) $\displaystyle\partial_{t}w=-\partial_{x}T-\frac{1}{\tau_{R}}f_{w}.$
(72b)

Let us assume, for simplicity, that we are restricted to a narrow temperature
range around a temperature $T_{0}$ and low values of $w$, so that the free
energy is quadratic as well

$f(T,w)\approx
f_{0}-s_{0}(T-T_{0})-\frac{c_{V}}{2T_{0}}(T-T_{0})^{2}+\frac{1}{2}\kappa
w^{2},$ (73)

where $f_{0}=f(T_{0},0)$, $s_{0}=-f_{T}$ at $T_{0}$ and $w=0$, and
$c_{V}=T\frac{\partial s}{\partial T}$ is the heat capacity. Then we have, in
particular, that $f_{TT}=-c_{V}/T_{0}=-\beta$, $\beta$ being a positive
constant, $f_{wT}=0$, and $f_{ww}=\kappa$. Note that the non-negative
coefficient $\kappa$ is a measure of the non-Fourier behavior. Equations (72)
can be then rewritten as

$\displaystyle\beta\partial_{t}T$
$\displaystyle=-\partial_{x}f_{w}+\frac{1}{\tau_{R}T}f_{w}^{2}=-\kappa\partial_{x}w+\frac{\kappa^{2}}{\tau_{R}T}w^{2}$
(74a) $\displaystyle\partial_{t}w$
$\displaystyle=-\partial_{x}T-\frac{1}{\tau_{R}}f_{w}=-\partial_{x}T-\frac{\kappa}{\tau_{R}}w.$
(74b)

Partial time-derivative of Equation (74a) then yields, assuming that
$\tau_{R}=const$,

$\displaystyle\beta\partial_{t}\partial_{t}T$
$\displaystyle=-\kappa\partial_{x}\left(-\partial_{x}T-\frac{\kappa}{\tau_{R}}w\right)+\frac{\kappa^{2}}{\tau_{R}T}\left(-\frac{\partial_{t}T}{T}w^{2}+2w\partial_{t}w\right)$
$\displaystyle=\kappa\partial_{x}\partial_{x}T-\frac{\kappa\beta}{\tau_{R}}\partial_{t}T+\mathcal{O}(\kappa^{2}\cdot
w).$ (75)

Assume that the terms of order $\kappa^{2}w$ (not necessarily the terms or
order $\kappa^{2}\partial_{x}w$) are smaller than the remaining terms, we
obtain the classical Maxwell-Cattaneo-Vernotte telegrapher’s equation,

$\beta\partial_{t}\partial_{t}T+\frac{\kappa\beta}{\tau_{R}}\partial_{t}T=\kappa\partial_{x}\partial_{x}T.$
(76)

In summary, after the transformation to variables $T$ and $\mathbf{w}$, we
obtain in the one-dimensional case and after linearization the classical
Maxwell-Cttaneo-Vernotte equation. However, Equations (36) are more general.

##  Appendix C Galilean invariance of phonon hydrodynamics

Consider an observer $\mathcal{O}$ in an inertial laboratory frame, who
measures time and space positions in $t$ and $\mathbf{r}$. Consider, moreover,
another observer $\mathcal{O}^{\prime}$ passing by the laboratory with a
constant velocity $\mathbf{V}$. Observer $\mathcal{O}^{\prime}$ measures time
and space positions $t$ and $\mathbf{r}^{\prime}=\mathbf{r}-\mathbf{V}t$.

An object with mass $M$ at position $\mathbf{r}^{\prime}$ is at position
$\mathbf{r}=\mathbf{r}^{\prime}+\mathbf{V}t$ in the laboratory frame. By
differentiation with respect to time, the velocities become related by
$\mathbf{v}=\mathbf{v}^{\prime}+\mathbf{V}$. Momentum of the object in the
laboratory frame, $\mathbf{p}=M\mathbf{v}$, relates to the momentum in the
$\mathcal{O}^{\prime}$ frame as $\mathbf{p}=\mathbf{p}^{\prime}+M\mathbf{V}$.
In particular, if the object has no mass, as a phonon, then its momentum does
not change, that is $\mathbf{m}=\mathbf{m}^{\prime}$. The $\mathbf{w}$ field
is equal to $\mathbf{m}/s$, and entropy itself is Galilean invariant
($s=s^{\prime}$), which means that the $\mathbf{w}$ field does not change
under Galilean transformations, $\mathbf{w}=\mathbf{w}^{\prime}$.

Energy of the object transforms as
$e=\frac{1}{2}M\mathbf{v}^{2}=e^{\prime}+\mathbf{p}^{\prime}\cdot\mathbf{V}+\frac{1}{2}M\mathbf{V}^{2}$,
which is a general formula for Galilean transformations of energy [54]. In
particular, if the object has no mass (as a phonon), then the formula
simplifies to $e^{\prime}=e-\mathbf{m}\cdot\mathbf{V}$, where $e$ stands for
the volumetric energy density and $\mathbf{m}$ for the volumetric momentum
density of phonons. In terms of the $\mathbf{w}$ field, the energy density
transforms as $e^{\prime}=e-s\mathbf{w}\cdot\mathbf{V}$, which means that
derivatives of energy transform as

$\frac{\partial e}{\partial s}=\frac{\partial e^{\prime}}{\partial
s^{\prime}}+\mathbf{w}^{\prime}\cdot\mathbf{V}\quad\text{and}\quad\frac{\partial
e}{\partial\mathbf{w}}=\frac{\partial
e^{\prime}}{\partial\mathbf{w}^{\prime}}+s\mathbf{V}.$ (77)

How do the time-and-space derivatives transform? For any function
$s(t,\mathbf{r})$ it holds that

$\displaystyle\left(\frac{\partial s}{\partial t}\right)_{\mathbf{r}}$
$\displaystyle=\left(\frac{\partial s}{\partial
t}\right)_{\mathbf{r}^{\prime}}-\left(\frac{\partial s}{\partial r^{\prime
j}}\right)_{t}\left(\frac{\partial r^{\prime j}}{\partial
t}\right)_{\mathbf{r}^{\prime}}=\left(\frac{\partial s}{\partial
t}\right)_{\mathbf{r}^{\prime}}-V^{j}\left(\frac{\partial s}{\partial
r^{\prime j}}\right)_{t}$ (78a) $\displaystyle\left(\frac{\partial
s}{\partial\mathbf{r}^{\prime}}\right)_{t}$
$\displaystyle=\left(\frac{\partial
s}{\partial\mathbf{r}^{\prime}}\right)_{t}\left(\frac{\partial\mathbf{r}^{\prime}}{\partial\mathbf{r}}\right)_{t}=\left(\frac{\partial
s}{\partial\mathbf{r}^{\prime}}\right)_{t}.$ (78b)

The reversible part of the evolution equations for the fields $s$ and
$\mathbf{w}$ then transform to

$\left(\frac{\partial s}{\partial
t}\right)_{\mathbf{r}^{\prime}}\cancel{-V^{j}\frac{\partial s}{\partial
r^{\prime j}}}=-\frac{\partial}{\partial r^{\prime i}}\left(\frac{\partial
e^{\prime}}{\partial w^{\prime}_{i}}+V^{i}s\right)=-\frac{\partial}{\partial
r^{\prime i}}\left(\frac{\partial e^{\prime}}{\partial
w^{\prime}_{i}}\right)\cancel{-V^{i}\frac{\partial s}{\partial r^{\prime i}}}$
(79)

and

$\left(\frac{\partial w^{\prime}_{i}}{\partial
t}\right)_{\mathbf{r}^{\prime}}\cancel{-V^{j}\frac{\partial
w^{\prime}_{i}}{\partial r^{\prime j}}}=-\frac{\partial}{\partial r^{\prime
i}}\left(\frac{\partial e^{\prime}}{\partial
s^{\prime}}+\cancel{V^{i}w^{\prime}_{i}}\right)+\frac{1}{s^{\prime}}\left(\frac{\partial
w^{\prime}_{j}}{\partial r^{\prime i}}-\frac{\partial w^{\prime}_{i}}{\partial
r^{\prime j}}\right)\left(\frac{\partial e^{\prime}}{\partial
w^{\prime}_{j}}+\cancel{V^{j}s}\right),$ (80)

which means that the reversible part of the equations is Galilean invariant.
Note, however, that without the last term on the right hand side, the
evolution equation for $\mathbf{w}$ would not be Galilean invariant. The
convective terms that we get by the Poisson reduction thus make the equations
Galilean invariant.

##  Appendix D Reversibility, irreversibility, and Onsager-Casimir reciprocal
relations

When talking about reversible and irreversible evolution, we shall first
define what we mean by (ir)reversibility. We use the definition based on the
time-reversal transformation (TRT) [58, 10]. TRT inverts velocities of all
particles and the magnetic field, and a quantity is called even with respect
to TRT (parity $\mathcal{P}=1$) if it is not affected by that inversion. On
the other hand, a quantity is odd (parity $\mathcal{P}=-1$) if its sign flips
under TRT. Note that there are also quantities that are neither even nor odd,
as for instance distribution functions, but the concept of parity still works
in a geometrical sense [68]. Most state variables, however, have a definite
parity, for instance energy, entropy, and density are even, whereas momentum
and velocity are odd.

Also, parity of a quantity can change. For instance, parity of the heat flux,
or the conjugate entropy flux $\mathbf{w}$, is odd if $\mathbf{w}$ is among
the state variables. On the other hand, after the reduction to the Fourier
heat conduction (Section 3.5), where only the energy density is the state
variable and $\mathbf{w}$ is proportional to the temperature gradient, parity
of $\mathbf{w}$ becomes even (the same as parity of the temperature gradient).

Consider now an evolution equation written in such a form that the time-
derivative of a state variable is on the left hand side. The right hand side
of the evolution equation is then interpreted as a component of a vector field
$X^{i}$. The part of the right hand side that transforms under TRT in the same
way as the left hand side is called reversible, while the part that flips the
sign is irreversible. In order to generate reversible dynamics, the Poisson
bivector must fulfill that

$TRT(L^{ij})=-\mathcal{P}(x^{i})\mathcal{P}(x^{j})L^{ij}.$ (81)

Since energy is even, this property indeed ensures that the Hamiltonian part
of the GENERIC evolution transforms as the left hand side of the evolution
equations and is thus reversible in the above sense. Property (81) is
fulfilled by the canonical Poisson bracket, from which the reversibility is
inherited by the Boltzmann Poisson bracket and all brackets obtained by its
projection (brackets within this manuscript) [68].

Entropy is required to be a Casimir of the Poisson bracket, that is
$\\{S,A\\}=0$ for all $A$, and the Hamiltonian evolution then does not produce
entropy (it can only redistribute it). The dissipative part is usually
required not to produce or destroy energy, so that $\langle
E_{\mathbf{x}},\Xi_{\mathbf{x}^{*}}\rangle|_{\mathbf{x}^{*}=S_{\mathbf{x}}}=0$.
These degeneracy requirements ensure that energy is conserved by the
irreversible part while entropy being conserved by the reversible part.

The dissipation potential is then typically required to be convex in order to
produce entropy, $\dot{S}=\langle
S_{\mathbf{x}},\Xi_{\mathbf{x}^{*}}\rangle|_{\mathbf{x}^{*}=S_{\mathbf{x}}}\geq
0$, see e.g. [27]. Note, however, that convexity is just a sufficient
condition of the positivity of entropy production, but not necessary [43].

The dissipation potential is supposed to be even with respect to TRT, which
means that it generates irreversible evolution. A part of the nonlinear
Onsager-Casimir reciprocal relations (OCRR) can be seen as the symmetry of the
Hessian of the dissipation potential [36]. However, to be closer to the
historical formulation of OCRR, let us assume, for a moment, that the
dissipation potential is quadratic as in Equation (2).

The requirement that $\Xi$ be even with respect to TRT means that

$TRT(M^{ij})=\mathcal{P}(x^{i})\mathcal{P}(x^{j})M^{ij}$ (82)

and the gradient part of the GENERIC evolution then indeed generates
irreversible evolution in the above sense. Operator $\mathbf{M}$ is often
referred to as the dissipative matrix, it is symmetric, and it can be obtained
as the second functional derivative of the dissipation potential,

$M^{ij}=\frac{\delta^{2}\Xi}{\delta x^{*}_{i}\delta x^{*}_{j}}.$ (83)

We are now in position to recall a generalized version of OCRR [66, 69].
GENERIC evolution can be rewritten in the form

$\partial_{t}x^{i}=\underbrace{\left(T_{0}L^{ij}-M^{ij}\right)}_{=K^{ij}}\frac{\partial\Phi}{\partial
x^{i}}$ (84)

with thermodynamic potential $\Phi=-S+\frac{1}{T_{0}}E$ ($T_{0}$ a constant).
The operator $\mathbf{K}$ in front of the derivatives of the potential $\Phi$
provides coupling between the state variables and it plays the role of the
Onsager-Casimir matrix of phenomenological coefficients. Let us now explicitly
discuss particular cases of the possible couplings. Due to properties (81) and
(82), we obtain the OCRR [66]:

  * •

Variables with the same parities, $\mathcal{P}(x^{i})=\mathcal{P}(x^{j})$, are
coupled through an operator symmetric with respect to the simultaneous TRT and
transposition.

  * •

Variables with opposite parities, $\mathcal{P}(x^{i})=-\mathcal{P}(x^{j})$,
are coupled through an operator skew-symmetric with respect to the
simultaneous TRT and transposition.

These properties are satisfied by both the reversible and irreversible parts
of GENERIC, in particular by the Hamiltonian mechanics, as can be checked
directly. Since GENERIC can be seen as a consequence of Hamiltonian mechanics
[70], OCRR can be seen as a consequence of Hamiltonian mechanics as well.

Let us now illustrate TRT on phonon hydrodynamics. The Poisson bivector for
the $(s,\mathbf{w})$ variables is

$\displaystyle L^{ss}$
$\displaystyle=\\{s(\mathbf{r}^{\prime}),s(\mathbf{r}^{\prime\prime})\\}=0$
(85a) $\displaystyle L^{sw_{k}}$
$\displaystyle=\\{s(\mathbf{r}^{\prime}),w_{k}(\mathbf{r}^{\prime\prime})\\}=\int
d\mathbf{r}\partial_{i}\delta(\mathbf{r}^{\prime}-\mathbf{r})\delta_{ki}\delta(\mathbf{r}^{\prime\prime}-\mathbf{r})=\partial^{\prime\prime}_{k}\delta(\mathbf{r}^{\prime}-\mathbf{r}^{\prime\prime})$
(85b) $\displaystyle L^{w_{k}s}$
$\displaystyle=\\{w_{k}(\mathbf{r}^{\prime}),s(\mathbf{r}^{\prime\prime}),\\}=-\int
d\mathbf{r}\partial_{i}\delta(\mathbf{r}^{\prime\prime}-\mathbf{r})\delta_{ki}\delta(\mathbf{r}^{\prime}-\mathbf{r})=-\partial^{\prime}_{k}\delta(\mathbf{r}^{\prime\prime}-\mathbf{r}^{\prime})$
(85c) $\displaystyle L^{w_{k}w_{l}}$
$\displaystyle=\\{w_{k}(\mathbf{r}^{\prime}),w_{l}(\mathbf{r}^{\prime\prime}),\\}=\int
d\mathbf{r}\frac{1}{s(\mathbf{r})}(\partial_{i}w_{j}-\partial_{j}w_{i})|_{\mathbf{r}}\delta_{ki}\delta(\mathbf{r}^{\prime}-\mathbf{r})\delta_{lj}\delta(\mathbf{r}^{\prime\prime}-\mathbf{r})$
$\displaystyle=\frac{1}{s(\mathbf{r}^{\prime})}(\partial^{\prime}_{k}w_{l}-\partial^{\prime}_{l}w_{k})|_{\mathbf{r}^{\prime}}\delta(\mathbf{r}^{\prime\prime}-\mathbf{r}^{\prime}).$
(85d)

We can check that the evolution generated by the bivector is reversible, since
it satisfies condition (81). Coupling between variables of the same parity,
for instance $L^{\mathbf{w}\mathbf{w}}$, is symmetric with respect to the
simultaneous transposition and time-reversal. Coupling between variables with
opposite parities, here $L^{s\mathbf{w}}$ and $L^{\mathbf{w}s}$, is skew-
symmetric in the same sense. Hamiltonian mechanics thus provides both the
symmetric and skew-symmetric coupling.

The irreversible evolution in Equations (36) can be considered as gradient
dynamics with dissipation potential

$\Xi=\frac{1}{2}\int d\mathbf{r}\xi(\mathbf{w}^{*})^{2},$ (86)

($xi$ being a constant) and the dissipative matrix can be obtained as the
second derivative of the dissipation potential,

$\mathbf{M}=\xi\delta(\mathbf{r}^{\prime}-\mathbf{r}^{\prime\prime})\begin{pmatrix}\left(\frac{\partial
s}{\partial\mathbf{w}}\right)^{2}&\frac{\partial s}{\partial\mathbf{w}}\\\
\frac{\partial s}{\partial\mathbf{w}}&\mathbf{1}\end{pmatrix}.$ (87)

Both $\mathbf{L}$ and $\mathbf{M}$ operators form the Onsager-Casimir operator
(84) and satisfy OCRR in the above sense. That is, the variables with the same
parity ($s$ and $s$, $\mathbf{w}$ and $\mathbf{w}$) are coupled through
operators that are symmetric with respect to the simultaneous transposition
(adjoint operator) and time reversal. Variables with opposite parities ($s$
and $\mathbf{w}$) are coupled through operators that are antisymmetric in the
same sense. Note that the Hamiltonian coupling contributes also to the
symmetric coupling ($\mathbf{w}$-$\mathbf{w}$) while the dissipative matrix
also contributes to the antisymmetric coupling ($s$-$\mathbf{w}$). In [66],
the composition of TRT and transposition were called bare and dressed
symmetries.

##  Appendix E Details of the Chapman-Enskog approximation

This Section contains details of the calculations needed in the first Chapman-
Enskog approximation from Section 3.3. In particular, it contains calculations
of the following integrals. The equilibrium value of the $Q_{ij}$ moment is
obtained by

$\displaystyle\int
d\mathbf{p}\frac{p_{i}p_{j}}{|\mathbf{p}|}\frac{1}{e^{e^{*}c|\mathbf{p}|/k_{B}}-1}=\frac{4\pi^{5}k_{B}^{4}}{45(e^{*})^{4}c^{4}}\delta_{ij},$
(88)

where which is calculated by observing that the integral vanishes if $i\neq
j$, by transformation to the spherical coordinates, and by expansion of the
fraction into a geometric series.

Integral

$\displaystyle\int
d\mathbf{p}\frac{p_{i}p_{j}}{|\mathbf{p}|}\frac{1}{e^{e^{*}c|\mathbf{p}|/k_{B}}e^{\mathbf{m}^{*}\cdot\mathbf{p}/k_{B}}-1}\approx\frac{4\pi^{5}k_{B}^{4}}{45(e^{*})^{4}c^{4}}\delta_{ij}$
(89)

is approximated by the Taylor expansion
$\frac{1}{ye^{x}-1}\approx\frac{1}{y-1}+\frac{yx}{(y-1)^{2}}$, where the zero-
th order term represents the equilibrium distribution function, $f_{R}$, and
the first-order term disappears after the integration.

Integral

$\displaystyle\int
d\mathbf{p}\frac{p_{i}p_{j}p_{k}}{|\mathbf{p}|}\frac{1}{e^{e^{*}c|\mathbf{p}|/k_{B}}e^{(\mathbf{m}^{*})^{l}p_{l}/k_{B}}-1}\approx-\frac{16\pi^{5}k_{B}^{4}}{225(e^{*}c)^{5}}\left(m^{*}_{i}\delta_{kj}+m^{*}_{k}\delta_{ji}+m^{*}_{j}\delta_{ik}\right)$
(90)

also starts with the above Taylor expansion, where the zero-th order term does
not contribute to the overall integral. The first-order term vanishes if an
odd power of $p_{i}$, $p_{j}$, $p_{k}$, or $p_{l}$ is present, which means
that the momenta must come either in two distinct pairs or as the fourth power
of a particular component of the momentum. Calculation in the spherical
coordinates then leads to the above result.

## References

  * [1] The theory of superfluidity of helium II.  In D. TER HAAR, editor, Collected Papers of L.D. Landau, chapter 46, pages 301 – 330. Pergamon, 1965.
  * [2] F. X. Alvarez, D. Jou, and A. Sellitto.  Phonon hydrodynamics and phonon-boundary scattering in nanosystems.  J. Appl. Phys, 105(014317), 2009.
  * [3] A. M. Anile, S. Pennisi, and M. Sammertino.  A thermodynamical approach to eddington factors.  J. Math. Phys., 32(544), 1991.
  * [4] C. Batchelor and G. Batchelor.  An Introduction to Fluid Dynamics.  Cambridge Mathematical Library. Cambridge University Press, 2000.
  * [5] A. Berezovski and P. Ván.  Internal Variables in Thermoelasticity.  Solid Mechanics and Its Applications. Springer International Publishing, 2017.
  * [6] A. Beris and B. Edwards.  Thermodynamics of Flowing Systems.  Oxford Univ. Press, Oxford, UK, 1994.
  * [7] P. L. Bhatnagar, E. P. Gross, and M. Krook.  A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems.  Phys. Rev., 94:511–525, May 1954.
  * [8] G. Boillat.  Sur l’éxistence et la recherche d’équations de conservation supplémentaires pour les systémes hyperboliques.  C. R. Acad. Sc. Paris, 278 A:909, 1974.
  * [9] J. Callaway.  Model for lattice thermal conductivity at low temperatures.  Phys. Rev., 113, 1959.
  * [10] H. B. G. Casimir.  On Onsager’s principle of microscopic reversibility.  Rev. Mod. Phys., 17:343–350, Apr 1945.
  * [11] C. Cattaneo.  Sulla conduzione del calore.  Atti del Seminario Matematico e Fisico della Universita di Modena, 3:83–101, 1948.
  * [12] S. Chapman, T. Cowling, D. Burnett, and C. Cercignani.  The Mathematical Theory of Non-uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases.  Cambridge Mathematical Library. Cambridge University Press, 1990.
  * [13] V. A. Cimmelli, A. Sellitto, and D. Jou.  Nonlinear evolution and stability of the heat flow in nanosystems: Beyond linear phonon hydrodynamics.  Phys. Rev. B, 82(184302), Nov 2010.
  * [14] S. R. de Groot and P. Mazur.  Non-equilibrium Thermodynamics.  Dover Publications, New York, 1984.
  * [15] R. Debye.  Zur theorie tier spezifischen wärmen.  Ann. Phys., 39(4):789, 1912.
  * [16] W. Dreyer and H. Struchtrup.  Heat pulse experiments revisited.  Continuum Mechanics and Thermodynamics, 5:3–50, 1993.
  * [17] M. Dumbser, I. Peshkov, E. Romenski, and O. Zanotti.  High order ADER schemes for a unified first order hyperbolic formulation of continuum mechanics: Viscous heat-conducting fluids and elastic solids.  Journal of Computational Physics, 314:824–862, 2016.
  * [18] C. Eckart.  The thermodynamics of irreversible processes I-III.  Phys. Rev., 58:267, 1940.
  * [19] M. Fecko.  Differential Geometry and Lie Groups for Physicists.  Cambridge University Press, 2006.
  * [20] L. Gesamtausgabe.  Ludwig Boltzmann Gesamtausgabe - Collected Works.  1983\.
  * [21] V. Ginzburg and L. Landau.  On the theory of superconductivity.  Zhur. Eksp. Theor. Fiz., 20:1064–1082, 1950.
  * [22] S. Godunov.  Symmetric form of the magnetohydrodynamic equation.  Chislennye Metody Mekhaniki Sploshnoi Sredy, 3(1):26–34, 1972.
  * [23] S. Godunov, T. Mikhailova, and E. Romenskii.  Systems of thermodynamically coordinated laws of conservation invariant under rotations.  Siberian Mathematical Journal, 37(4):690–705, 1996.
  * [24] M. Grmela.  Particle and bracket formulations of kinetic equations.  Contemporary Mathematics, 28:125–132, 1984.
  * [25] M. Grmela.  Bracket formulation of diffusion-convection equations.  Physica D, 21:179–212, 1986.
  * [26] M. Grmela.  Thermodynamics of driven systems.  Phys. Rev. E, 48:919–930, Aug 1993.
  * [27] M. Grmela.  Multiscale equilibrium and nonequilibrium thermodynamics in chemical engineering.  Advances in Chemical Engineering, 39:76–128, 2010.
  * [28] M. Grmela.  GENERIC guide to the multiscale dynamics and thermodynamics.  J. Phys. Commun., 2(032001), 2018.
  * [29] M. Grmela, D. Jou, and J. Casas-Vázquez.  Radiation hydrodynamics and thermodynamics.  23(3):203–225, 1998.
  * [30] M. Grmela, G. Lebon, and C. Dubois.  Multiscale thermodynamics and mechanics of heat.  Physical Review E, 83(061134), 2011.
  * [31] M. Grmela and H. C. Öttinger.  Dynamics and thermodynamics of complex fluids. I. Development of a general formalism.  Phys. Rev. E, 56:6620–6632, Dec. 1997.
  * [32] Y. Guo and M. Wang.  Phonon hydrodynamics and its applications in nanoscale heat transport.  Physics Reports, 595:1–44, 2015.  Phonon hydrodynamics and its applications in nanoscale heat transport.
  * [33] Y. Guo and M. Wang.  Phonon hydrodynamics for nanoscale heat transport at ordinary temperatures.  Phys. Rev. B, 97(035421), Jan. 2018.
  * [34] Y. Guo, Z. Zhang, M. Nomura, S. Volz, and M. Wang.  Phonon vortex dynamics in graphene ribbon by solving boltzmann transport equation with ab initio scattering rates.  International Journal of Heat and Mass Transfer, 169(120981), 2021\.
  * [35] R. A. Guyer and J. A. Krumhansl.  Solution of the linearized phonon Boltzmann equation.  Physical review, 148(2), 1966.
  * [36] I. Gyarmati.  Non-equilibrium thermodynamics: Field theory and variational principles.  Engineering Science Library. Springer, 1970.
  * [37] R. J. Hardy and S. S. Jaswal.  Velocity of second sound in naf.  Phys. Rev. B, 3:4385–4387, Jun 1971.
  * [38] J. Harrison, G. Lombardo, and P. Peressini.  Specific heat of naf.  Journal of Physics and Chemistry of Solids, 29(3):557–559, 1968\.
  * [39] D. D. Holm.  Introduction to HVBK Dynamics, pages 114–130.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2001.
  * [40] D. D. Holm, J. E. Marsden, T. Ratiu, and A. Weinstein.  Nonlinear stability of fluid and plasma equilibria.  Physics Reports, 123(1 & 2):1–116, 1985.
  * [41] M. Hütter and B. Svendsen.  Quasi-linear versus potential-based formulations of force–flux relations and the generic for irreversible processes: comparisons and examples.  Continuum Mechanics and Thermodynamics, 25(6):803–816, 2013.
  * [42] B. JA, C. JG, and W. GK.  Heat capacities and grüneisen parameters for NaF and NaBr below 25 k.  Australian Journal of Physics, 32:463–468, 1979.
  * [43] A. Janečka and M. Pavelka.  Non-convex dissipation potentials in multiscale non-equilibrium thermodynamics.  Continuum Mechanics and Thermodynamics, 30(4):917–941, 2018.
  * [44] M. J.E.  A group theoretical approach to the equations of plasma physics.  Canad. Math. Bull., 25:129–142, 1982.
  * [45] D. Jou, J. Casas-Vázquez, and G. Lebon.  Extended Irreversible Thermodynamics.  Springer-Verlag, New York, 4th edition, 2010.
  * [46] D. Jou and L. Restuccia.  Nonlinear heat transport in superlattices with mobile defects.  Entropy, 21(12), 2019.
  * [47] D. Jou and L. Restuccia.  Non-linear heat transport effects in systems with defects.  Journal of Non-Equilibrium Thermodynamics, 47(2):179–186, 2022\.
  * [48] I. Khalatnikov.  An Introduction To The Theory Of Superfluidity.  Avalon Publishing, 1989.
  * [49] C. Kittel, H. Charles Kittel, K. Charles, H. Kroemer, and K. Herbert.  Thermal Physics.  W. H. Freeman, 1980.
  * [50] R. Kovács and P. Ván.  Generalized heat conduction in heat pulse experiments.  International Journal of Heat and Mass Transfer, 83:613 – 620, 2015\.
  * [51] L. Landau and E. Lifschitz.  Statistical physics.  Number pt. 1 in Course of theoretical physics. Pergamon Press, 1969.
  * [52] L. Landau and E. Lifshitz.  Quantum Mechanics: Non-relativistic Theory.  Butterworth Heinemann. Butterworth-Heinemann, 1977.
  * [53] W. Larecki and S. Piekarski.  Phonon gas hydrodynamics based on the maximum entropy principle and the extended field theory of a rigid conductor of heat.  Arch. Mech., 43(2-3):163, 1992.
  * [54] M. Le Bellac and J. M. Lévy-Leblond.  Galilean electromagnetism.  Il Nuovo Cimento B (1971-1996), 14(2):217– 234, 1973.
  * [55] G. Lebon, D. Jou, and P. C. Dauby.  Beyond the fourier heat conduction law and the thermal non-slip condition.  Phys Lett. A, 376:2842–2846, 2012.
  * [56] G. Lebon, D. Jou, and J. Vázquez.  Understanding Non-Equilibrium Thermodynamics: Foundations, Applications, Frontiers.  SpringerLink: Springer e-Books. Springer London, Limited, 2008.
  * [57] T. McNelly, S. Rogers, D. Channin, R. Rollefson, W. Goubau, G. Schmidt, J. Krumhansl, and R. Pohl.  Heat pulses in NaF: onset of second sound.  Phys. Rev. Lett., 24:100–102, 1970.
  * [58] J. Meixner and H. Reik.  Thermodynamik der Irreversible Prozesse, in Handbuch der Physik, volume 3/II.  Springer, Berlin Heidelberg New York, 1959.
  * [59] M. Mongioví.  Superfluidity and entropy conservation in extended thermodynamics.  J. Non-Equilib. Thermodyn., 16:225–240, 1991.
  * [60] M. S. Mongiovì and L. Restuccia.  Hydrodynamical model of anisotropic, polarized turbulent superfluids. i: constraints for the fluxes.  Zeitschrift für angewandte Mathematik und Physik, 69(1):2, Nov 2017.
  * [61] M. S. Mongiovì, D. Jou, and M. Sciacca.  Non-equilibrium thermodynamics, heat transport and thermal waves in laminar and turbulent superfluid helium.  Physics Reports, 726:1–71, 2018.
  * [62] P. J. Morrison.  Bracket formulation for irreversible classical fields.  Physics Letters A, 100(8):423–427, 1984.
  * [63] I. Müller and T. Ruggeri.  Rational Extended Thermodynamics.  Springer, New York, 2nd edition, 1998.
  * [64] L. Onsager.  Reciprocal relations in irreversible processes. I.  Phys. Rev., 37:405–426, Feb. 1931.
  * [65] L. Onsager.  Reciprocal relations in irreversible processes. II.  Phys. Rev., 38:2265–2279, Dec 1931.
  * [66] H. Öttinger.  Beyond Equilibrium Thermodynamics.  Wiley, 2005.
  * [67] H. C. Öttinger and M. Grmela.  Dynamics and thermodynamics of complex fluids. II. Illustrations of a general formalism.  Phys. Rev. E, 56:6633–6655, Dec. 1997.
  * [68] M. Pavelka, V. Klika, and M. Grmela.  Time reversal in nonequilibrium thermodynamics.  Phys. Rev. E, 90(062131), 2014.
  * [69] M. Pavelka, V. Klika, and M. Grmela.  Multiscale Thermo-Dynamics.  de Gruyter (Berlin), 2018.
  * [70] M. Pavelka, V. Klika, and M. Grmela.  Generalization of the dynamical lack-of-fit reduction.  Journal of Statistical Physics, 181(1):19–52, 2020.
  * [71] M. Pavelka, I. Peshkov, and V. Klika.  On Hamiltonian continuum mechanics.  Physica D: Nonlinear phenomena, 408(132510), 2020.
  * [72] R. Peierls.  Zur kinetischen theorie der wärmeleitung in kristallen.  Ann. Phys., 3:1055, 1929.
  * [73] R. Peierls.  Quantum Theory of Solids.  Oxford University Press, 1955.
  * [74] I. Peshkov, M. Pavelka, E. Romenski, and M. Grmela.  Continuum mechanics and thermodynamics in the Hamilton and the Godunov-type formulations.  Continuum Mechanics and Thermodynamics, 30(6):1343–1378, 2018.
  * [75] M. Raya-Moreno, J. Carrete, and X. Cartoixà.  Hydrodynamic signatures in thermal transport in devices based on two-dimensional materials: An ab initio study.  Phys. Rev. B, 106:014308, Jul 2022.
  * [76] E. I. Romensky.  Hyperbolic systems of thermodynamically compatible conservation laws in continuum mechanics.  Math. Comput. Modelling, 28(10):115–130, 1998.
  * [77] T. Ruggeri and A. Strumia.  Main field and convex covariant density for quasi-linear hyperbolic systems.  Ann. Inst. H. Poincaré, 34:65, 1981.
  * [78] L. Saint-Raymond.  A mathematical PDE perspective on the Chapman–Enskog expansion.  Bulletin (New Series) of the American Mathematical Society, 51(2):247–275, 2014.
  * [79] A. Sellitto, F. X. Alvarez, and D. Jou.  Phonon-wall interactions and frequency-dependent termal conductivity in nanowires.  J. Appl. Phys, 109(064317), 2011.
  * [80] A. Sellitto, I. Carlomagno, and D. Jou.  Two-dimensional phonon hydrodynamics in narrow strips.  Proc. Roy. Soc. A, 471(20150376), 2015.
  * [81] L. Sendra, A. Beardo, P. Torres, J. Bafaluy, F. X. Alvarez, and J. Camacho.  Derivation of a hydrodynamic heat equation from the phonon Boltzmann equation for general semiconductors.  Phys. Rev. B, 103:L140301, Apr 2021.
  * [82] M.-Y. Shang, C. Zhang, Z. Guo, and J.-T. Lü.  Heat vortex in hydrodynamic phonon transport of two-dimensional materials.  Scientific Reports, 10(1):8272, 2020.
  * [83] M. Szücs, R. Kovács, and S. Simić.  Open mathematical aspects of continuum thermodynamics: Hyperbolicity, boundaries and nonlinearities.  Symmetry, 12(9), 2020.
  * [84] M. Szücs, M. Pavelka, R. Kovács, T. Fülöp, P. Ván, and M. Grmela.  A case study of non-fourier heat conduction using internal variables and generic.  Journal of Non-Equilibrium Thermodynamics, 47(1):31–60, 2021.
  * [85] M. Sýkora, M. Pavelka, M. La Mantia, D. Jou, and M. Grmela.  On the relations between large-scale models of superfluid helium-4.  Physics of Fluids, 33(12):127124, 2021.
  * [86] L. Tisza.  Transport phenomena in helium II.  Nature, 141(3577):913–913, 1938.
  * [87] P. Torres, A. Ziabari, A. Torelló, J. Bafaluy, J. Camacho, X. Cartoixà, A. Shakouri, and F. X. Alvarez.  Emergence of hydrodynamic heat transport in semiconductors at the nanoscale.  Phys. Rev. Materials, 2:076001, Jul 2018.
  * [88] C. Truesdell.  Rational thermodynamics.  Springer, New York, NY USA, 2nd ed. edition, 1984.
  * [89] S. Van Sciver.  Helium Cryogenics.  International Cryogenics Monograph Series. Springer New York, 2012.
  * [90] M. Wong, C. Tso, T. Ho, and H. Lee.  A review of state of the art thermal diodes and their potential applications.  International Journal of Heat and Mass Transfer, 164:120607, 2021\.
  * [91] C. Zhang, S. Chen, and Z. Guo.  Heat vortices of ballistic and hydrodynamic phonon transport in two-dimensional materials.  International Journal of Heat and Mass Transfer, 176:121282, 2021\.