rarefied (low-pressure) gases, which needs to generalize Newton’s law of
fluids, too [266, 267, 268, 269].

Their kinetic interpretation helps understand their common property: a
significant part of the energy carriers (such as phonons or actual molecules)
can propagate without collisions, i.e., the mean free path becomes comparable
to the characteristic length of the object. It does not mean that all the
energy carriers propagate that way, but the ballistic effects notably
contribute to the overall energy transport. In the cases mentioned above, the
observations of ballistic propagation are performed differently as it is not
the temperature history only for which the ballistic effects manifest or the
temperature itself is not directly measurable, contrary to the room-
temperature heat pulse experiments.

###  5.1. Modeling of low-temperature heat conduction.

Let us recall Figure 1, which shows the temperature history recorded in a
remarkably intriguing heat pulse experiment, called NaF experiments, performed
by McNelly et al. [1]. Interestingly, all propagation modes are present in a
single experiment: the longitudinal and transversal ballistic propagation
modes (which waves are the fastest), the second sound, and diffusion. That
remarkable outcome makes this experiment a benchmark problem in this field.

####  5.1.1. Phonon models.

Starting with a kinetic background, we note that, unfortunately, the thermo-
mechanical GK system (61)-(63) is not tested; thus, it is unknown whether this
coupled model is suitable to model these experiments. On the contrary, the
momentum series expansion of the Boltzmann equation, Eq. (50), recalled here,

$\frac{\partial u_{\langle m\rangle}}{\partial
t}+\frac{m^{2}}{4m^{2}-1}c\frac{\partial u_{\langle m-1\rangle}}{\partial
x}+c\frac{\partial u_{\langle m+1\rangle}}{\partial
x}=\left\\{\begin{array}[]{ll}\displaystyle 0&\ m=0\\\
-\frac{1}{\tau_{R}}u_{\langle 1\rangle}&\ m=1\\\
-\left(\frac{1}{\tau_{R}}+\frac{1}{\tau_{N}}\right)u_{\langle m\rangle}&\
2\leq m\leq M\end{array}\right.,$

performs well with minor discrepancies, which mostly originate from the
properties of that approach. First, one has to utilize at least $M=30$
momentum equations to obtain a good approximation for the speed of sound
(ballistic speed of phonons), that is, for a general three-dimensional problem
a 30${}^{\textrm{th}}$-order tensor would appear in the model. Moreover, in
principle, one would need infinitely many momentum equations for the exact
value [32]. Despite these shortcomings, the phonon hydrodynamic equations can
be reasonably simplified, and the moments up to $M=3$ provide an acceptable
approximation [3], which reads

$\displaystyle\dot{e}+c^{2}p_{x}$ $\displaystyle=0,$ (87)
$\displaystyle\dot{p}+\frac{1}{3}e_{x}+N_{x}$
$\displaystyle=-\frac{1}{\tau_{R}}p,$ (88)
$\displaystyle\dot{N}+\frac{4}{15}c^{2}p_{x}$
$\displaystyle=-\left(\frac{1}{\tau_{R}}+\frac{1}{\tau_{N}}\right)N,$ (89)

where

$e=hcu,\quad p_{i}=hu_{i},\quad N_{\langle ij\rangle}=hcu_{\langle
ij\rangle},$ (90)

are the energy density, momentum density, and the deviatoric part of the
pressure tensor, respectively, have been found for the one-dimensional case in
Eq. (50).

Staying with the phonon hydrodynamic background, we have to mention here the
work of Ma [270, 271] as well, who developed a so-called hybrid phonon gas
model based on the papers of Rogers [272] and Landau [273]. That work was
motivated to include both the longitudinal and transversal ballistic modes
simultaneously with keeping a one-dimensional model instead of a proper three-
dimensional description. To achieve this goal, the internal energy $e$ is
divided into two parts, $e=e_{l}+e_{t}$, together with the heat flux
$q=q_{l}+q_{t}$, each of them being connected to the corresponding
longitudinal and transversal ballistic mode. Furthermore, following Rogers, it
is supposed that the heat flux $\mathbf{q}$ is proportional to the velocity of
phonon gas, hence $\mathbf{q}=e\mathbf{v}$, and the classical equation of
motion, the Navier-Stokes equation,

$\displaystyle\rho\partial_{t}\mathbf{v}+\rho(\mathbf{v}\cdot\nabla)\mathbf{v}=-\nabla
p+\eta\Delta\mathbf{v}+\left(\xi+\frac{1}{3}\eta\right)\nabla\nabla\cdot\mathbf{v},$
(91)

is modified accordingly and reduced to one dimension. Here, $p$ denotes the
scalar pressure, $\eta$ and $\xi$ are for the shear and bulk viscosity.
According to Rogers, in such a low-temperature situation, the shear viscosity
tends to zero when the normal processes dominate the propagation.
Consequently, bulk viscosity plays an essential role in the damping mechanism.
For that model, Landau’s complex viscosity equation is adopted,

$\displaystyle\xi=\frac{2\tau e_{t}}{3(1-i\omega\tau)},$ (92)

in which $\tau^{-1}=\tau_{R}^{-1}+\tau_{N}^{-1}$, and $\omega$ stands for the
angular frequency of heat waves, and obtains the evolution equation

$\displaystyle\partial_{t}q+\frac{1}{\tau_{R}}q=-\frac{1}{3}\partial_{x}e+\frac{2\tau}{3(1-i\omega\tau)}\partial_{xx}q.$
(93)

Eventually, Eq. (93) stands for both $q_{t}$ and $q_{l}$, with separate
balances for the internal energies $e_{l}$ and $e_{t}$, and the relation
$e=e_{l}+e_{t}$ realizes the coupling to obtain the temperature evolution.
This appears to be a particular extension of a two-temperature model in which
two constitutive equations of GK-type are coupled, thus it is neither a pure
phonon nor a continuum model.

####  5.1.2. Continuum models.

In the analogy of the MCV and GK equations, it is also possible to achieve
analogy with the phonon hydrodynamic approach (87)-(89), in the same way as
the GK equation is derived [236], alternatively, in a GENERIC framework [168].
However, as Eq. (89) suggests, we need to include a second-order tensor in the
state space as well, viz.,
$s(e,\mathbf{q},\mathbf{Q})=s_{\textrm{le}}(e)-m_{1}/2\mathbf{q}^{2}-m_{2}/2\mathbf{Q}:\mathbf{Q}$.
Concerning the entropy flux, Eq. (65) remains suitable, and the current
multiplier $\mathbf{B}$ is helpful to realize the coupling between
$\mathbf{q}$ and $\mathbf{Q}$. It becomes more transparent in the
corresponding entropy production,

$\displaystyle\sigma_{s}=\mathbf{q}\cdot\left(-\rho
m_{1}\partial_{t}\mathbf{q}+\nabla\frac{1}{T}\mathbf{I}+\nabla\cdot\mathbf{B}\right)-\rho
m_{2}\mathbf{Q}:\partial_{t}\mathbf{Q}+\mathbf{B}:\nabla\mathbf{q},$ (94)

resulting in the Onsagerian relations

$\displaystyle-\rho
m_{1}\partial_{t}\mathbf{q}+\nabla\frac{1}{T}+\nabla\cdot\mathbf{B}$
$\displaystyle=l\mathbf{q},$ (95) $\displaystyle-\rho
m_{2}\partial_{t}\mathbf{Q}$
$\displaystyle=L_{11}\mathbf{Q}+L_{12}\nabla\mathbf{q},$ (96)
$\displaystyle\mathbf{B}=L_{21}\mathbf{Q}+L_{22}\nabla\mathbf{q},$ (97)

with $l\geq 0$, $L_{11}\geq 0$, $L_{22}\geq 0$, and
$L_{11}L_{22}-L_{12}L_{21}\geq 0$, expressing the positive semi-definiteness
of the inequality (94). For the compatibility, we take $L_{22}=0$, thus
identifying $\mathbf{B}=L_{21}\mathbf{Q}$. In that approach, $\mathbf{Q}$ is
no longer directly proportional with $\nabla\mathbf{q}$ on the contrary to the
continuum GK model, Eq. (96) describes its time evolution, The removal of the
non-equilibrium contribution of $\mathbf{Q}$ from the entropy density
naturally relaxes $\mathbf{Q}$ to $\nabla\mathbf{q}$. That picture is also
helpful in the interpretation of the current multiplier. Consequently, Eq.
(74) is extended with the time derivative of $\mathbf{Q}$, with a relaxation
time denoted by $\tau_{Q}$. Reducing the system (95)-(97) to the one spatial
dimension leads to

$\displaystyle\tau_{q}\partial_{t}q+q$
$\displaystyle=-\lambda\partial_{x}T+l\partial_{x}Q,$ (98)
$\displaystyle\tau_{Q}\partial_{t}Q+Q$ $\displaystyle=l\partial_{x}q,$ (99)

where $\tau_{q}$ and $\tau_{Q}$ are the corresponding relaxation times. Eqs.
(98)-(99) found as a hyperbolic approximation of a more general parabolic set
of equations [236]. Its complete isotropic representation is presented in
[274], showing the detailed relationships between the coefficients for various
situations. Here, we emphasize that these relaxation times are independent of
each other, and the compatibility with the phonon model (87)-(89) is formal
but still comparable [236, 275]. One can apply the kinetic interpretation,
which is a matter of choice, but not the only option. When one attempts to
model the NaF experiments, for instance, then the coefficient $l$ can be
adjusted to obtain the exact speed of sound of ballistic mode, together with
the relaxation times; overall, having the same number of parameters to be fit
[276, 277]. Figure 9 summarizes the modeling capabilities on experimental data
of different models. That difference becomes more crucial for rarefied (real)
gases, and what flexibility appears to be an advantage here will become a
disadvantage there. Therefore, no such ‘universal’ model can be suggested for
any heat conduction problems, but there are alternatives from which one can
choose the best one for a particular task.

Evaluation of the NaF experiments with both the phonon and continuum models
showed the temperature dependence of relaxation times, besides the thermal
conductivity [277]. Interestingly, as mentioned in the case of the MCV
equation, the direct implementation of such nonlinearities would open new
perspectives toward modeling low-temperature phenomena.

Figure 9. McNelly’s NaF experiments are evaluated with (A) phonon [32, 277],
(B) continuum [277], and (C) hybrid phonon models [271]. Part (B) also shows
the originally recorded data for comparison.

It also stands as further motivation to investigate thermo-mechanical models,
sharing that perspective with Frischmuth and Cimmelli [135, 136, 278]. Here,
we restrict ourselves to the small deformation regime. They attempted to
improve the modeling capabilities of the semi-empirical temperature approach
by introducing mechanics through thermal expansion. While the mechanical
subsystem remains completely classical, they modify the thermal part by
exchanging the heat flux $\mathbf{q}$ to $\mathbf{g}=\nabla\beta$, postulated
as a new state variable, for which an MCV-like equation is supposed to be
valid such as

$\displaystyle C\partial_{t}T$ $\displaystyle=\lambda\nabla\cdot\mathbf{g}+r,$
(100) $\displaystyle\tau\partial_{t}\mathbf{g}$ $\displaystyle=\nabla
T-\mathbf{g},$ (101)

where $r$ stands for a source term, including the mechanical contribution
$r\sim T_{0}\nabla\cdot\mathbf{v}$, and $C$ is specific heat capacity. We note
that the energy balance (100) is also arbitrarily modified, its thermodynamic
compatibility is unclear. However, the modification was necessary to match the
units in Eq. (101), otherwise, the thermal conductivity disappears, and
eliminating $\mathbf{g}$ still can recover the $T$-representation of the MCV
equation in the form

$\displaystyle\tau\partial_{tt}T+\partial_{t}T=\alpha\Delta
T+\frac{r}{C}+\tau\frac{\dot{r}}{C},$ (102)

but Eq. (102) cannot reflect the underlying assumptions and model details.
More importantly, it cannot justify the thermodynamic compatibility of the
system (100)-(101). Eq. (101) must be supplemented with the dynamic equation
$\dot{\beta}=f(T,\beta)$, adding further flexibility to the model. Sadly, they
did not manage to compare its solutions directly to experimental data, but in
a test solution, both the ballistic and second sound modes are present.

####  5.1.3. Notes on the thermo-mechanical couplings.

Despite the previous model’s questionable points, it is reasonable to couple
thermal expansion to the MCV equation [191]. Keeping thermodynamic
compatibility in mind, we start with the internal energy $e$ (in 1D) by
recalling Eq. (41)

$\displaystyle
e=cT+\frac{E}{2\rho}\varepsilon^{2}+\frac{E\chi}{\rho}T_{0}\varepsilon,$ (103)

where, we recall that, $\chi$ stands for the thermal expansion coefficient,
and $E$ is Young’s modulus. We note that for a low-temperature situation, this
could be valid for minimal temperature excitation for a narrow temperature
interval due to the emerging nonlinearities. This is avoided here.
Furthermore, in a general three-dimensional treatment, separating the
spherical and deviatoric terms of the strain tensor $\boldsymbol{\varepsilon}$
would be expedient. The mechanical contribution in the energy balance appears
as a $\sigma\partial_{t}\varepsilon$ source term

$\displaystyle\rho\partial_{t}e+\partial_{x}q=\sigma\partial_{t}\varepsilon,$
(104)

in which the stress $\sigma=E\varepsilon-E\chi(T-T_{0})$, and that energy
balance is coupled to the system

$\displaystyle\rho\partial_{t}v-\partial_{x}\sigma$ $\displaystyle=0,$ (105)
$\displaystyle\partial_{t}\varepsilon$ $\displaystyle=\partial_{x}v,$ (106)
$\displaystyle\tau\partial_{t}q+q$ $\displaystyle=-\lambda\partial_{x}T.$
(107)

In summary, this is a classical thermo-mechanical model including isotropic
Hooke’s elasticity for stress, and exchanging Fourier’s law to the MCV
equation, referred to as MCV-TE model in [191]. It is worth investigating the
temperature representation for a simplified ($\tau=0$, i.e., Fourier), but a
three-dimensional system, derived in [279],

$\displaystyle\frac{1}{c_{||}^{2}}\partial_{tt}(\gamma_{1}\partial_{t}T-\lambda\Delta
T)=\Delta\left[\left(\rho
c-\frac{6E^{\textrm{dev}}E^{\textrm{sph}}\chi^{2}T_{0}}{E^{\textrm{sph}}+2E^{\textrm{dev}}}\right)\partial_{t}T-\lambda\Delta
T\right],$ (108)

with parameters

$\displaystyle
c^{2}_{||}=\frac{E^{\textrm{sph}}+2E^{\textrm{dev}}}{3\rho},\quad\gamma_{1}=\rho
c-3E^{\textrm{sph}}\chi^{2}T_{0},\quad\gamma_{2}=\rho
c-\frac{6E^{\textrm{dev}}E^{\textrm{sph}}\chi^{2}T_{0}}{E^{\textrm{sph}}+2E^{\textrm{dev}}},\quad
E^{\textrm{sph}}=3K,\quad E^{\textrm{dev}}=2G$ (109)

in which $c_{||}$ is the longitudinal elastic wave propagation speed, $K$ and
$G$ are the bulk and shear moduli. Interestingly, the classical thermal
diffusivity $\alpha=\lambda/(\rho c)$ is modified by thermal expansion as
$\hat{\alpha}=\alpha(\rho c/\gamma_{2})$; if the thermal expansion coefficient
$\chi=0$, $\gamma_{1}=\gamma_{2}=\rho c$ follows, and the correction factor
$(\rho c/\gamma_{2})=1$ holds. In other words, neglecting thermal expansion
effects in the evaluation of the experimental data detunes the thermal
diffusivity, e.g., by 6% for aluminum.

Considering the NaF experiments, the coefficients $\rho(T)$ and $\chi(T)$ are
unknown, and therefore it seems reasonable to find $\rho$ and $\chi$ for each
reference temperature, but this has not been done so far. Similarly, for Eqs.
(98)-(99), one must apply the measured speeds of first and second sounds as
constraints for the fitted parameters, and that again reduces the number of
free parameters. The analytical and numerical solutions of ballistic equations
are analogous to the previous models, and the staggered spatial discretization
is more advantageous from a numerical point of view as it eases the
implementation of boundary conditions.

In fact, the MCV-TE model (105)-(107) is not a completely novel idea. Lord and
Shulman [280] were the first who proposed to couple the MCV equation with
thermal expansion. At that time, much before the work of Coleman et al. [125],
they formulated the anisotropic MCV model as

$\displaystyle\mathbf{A}\cdot\dot{\mathbf{q}}+a\dot{\mathbf{q}}+\mathbf{q}=b\nabla
T+\mathbf{B}\cdot\nabla T,$ (110)

where $a,b,\mathbf{A}$, and $\mathbf{B}$ are all considered as material
parameters so that the requirements for such model provided later by [125] are
missing (see Eqs. (38) and (39)). Furthermore, these parameters are not
necessarily positive, e.g., Fourier’s law is recovered with $b=-\lambda$,
($\lambda>0$). However, only the simplest isotropic form is used later with
constant scalar coefficients ($\mathbf{A}=\mathbf{0}$, $a=\tau$, $b=-\lambda$,
$\mathbf{B}=\mathbf{0}$). In the state space, they used the small strain
tensor $\boldsymbol{\varepsilon}$ and $T$, so that
$e=e(T,\boldsymbol{\varepsilon})$, and $s=s(T,\boldsymbol{\varepsilon})$. The
energy balance is formulated with the help of the Helmholtz free energy
$\varphi(T,\boldsymbol{\varepsilon})=e(T,\boldsymbol{\varepsilon})-Ts(T,\boldsymbol{\varepsilon})$,

$\displaystyle\nabla\cdot\mathbf{q}=\rho
T\left(\frac{\partial^{2}\varphi}{\partial
T^{2}}\dot{T}+\frac{\partial^{2}\varphi}{\partial\boldsymbol{\varepsilon}\partial
T}\dot{\boldsymbol{\varepsilon}}\right).$ (111)

Eliminating $\mathbf{q}$ using Eqs. (24) and (111) leads to a complex
nonlinear model in which terms $\dot{T}^{2}$,
$\dot{\boldsymbol{\varepsilon}}\dot{T}$, and
$\dot{\boldsymbol{\varepsilon}}^{2}$ are neglected, and yields

$\displaystyle\lambda\Delta T=-\rho T\left[\frac{\partial^{2}\varphi}{\partial
T^{2}}(\dot{T}+\tau\ddot{T})+\frac{\partial^{2}\varphi}{\partial{\boldsymbol{\varepsilon}}\partial
T}(\dot{\boldsymbol{\varepsilon}}+\tau\ddot{\boldsymbol{\varepsilon}})\right],$
(112)

that can be solved together with the equation of motion and with a particular
free energy function, resulting in a more complicated structure for the time
evolution than Eq. (108). Dhaliwal and Sherief utilize the same concept [281]
with anisotropic thermal conductivity, though with scalar relaxation time, for
which uniqueness is proved for homogeneous boundary conditions. That property
is inherited to (112), too. Compared to the MCV-TE model, the one-dimensional
versions can be identical; thus, the Lord-Shulman theory is helpful in
modeling ballistic heat conduction experiments. Since the heat flux evolution
equation is arbitrarily chosen, further, more general thermoelastic models can
be derived on the same basis.

The next class of thermoelastic models is called temperature-rate dependent
equations, where the time derivative of the temperature $\dot{T}$ also
influences the stress and heat flux evolutions. This is developed by Green and
Lindsay [282]. Considering the material to be anisotropic with constant
coefficients, the heat flux, and stress are

$\displaystyle\mathbf{q}=-T_{0}(\mathbf{b}\dot{T}+\mathbf{K}\nabla
T),\quad\boldsymbol{\sigma}=\mathbf{C}^{(4)}\boldsymbol{\varepsilon}+\mathbf{A}(T-T_{0}+a\dot{T}),$
(113)

where $a,\mathbf{b}$ and $\mathbf{A}$ are material parameters, $a$ might be
related to the relaxation time, but basically, these parameters have an
undiscovered physical interpretation.

The concept of elastic heat flow vector field by Ignaczak [283] is also worth
mentioning. This is based on the assumption that an additional vector field
$\boldsymbol{\beta}$ influences the resulting heat flux (similarly to an
internal variable),

$\displaystyle\mathbf{q}=-\lambda\nabla
T+\boldsymbol{\beta},\quad\dot{\boldsymbol{\beta}}=-a\frac{\nabla T}{T}$ (114)

with $a>0$ parameter. Consequently, both the internal energy and entropy
depend on $\boldsymbol{\beta}$, having importance in modeling low-temperature
wave propagation. This inevitably results in a nonlinear evolution for the
temperature field, specifically developed to find thermal soliton-like
solutions. In a one-dimensional case, two characteristic wave speeds are
propagating in the same direction, therefore describing different modes of
heat conduction. However, that theory is not yet directly tested on ballistic
experimental data.

For a more thorough review and for further details with mathematical proofs,
let us refer to [284, 285, 286, 287], and [288]. Furthermore, we also want to
highlight the book of Berezovski and Ván [42] for a broader discussion on the
role of internal variables and thermodynamics concerning thermoelasticity, the
corresponding numerical techniques are presented, too. In regard to
thermomechanics, it is worth mentioning that from an effective point of view,
although on significantly longer time scales, rheology can also stand as a
factor for non-Fourier heat conduction [192, 289]. For a profound mathematical
treatment of even more complex situations related to thermopiezoelecricity,
requiring anisotropy and further couplings with the electric field, let us
refer to the recent works of Tibullo et al. [290, 291]. These topics still
have immense research potential. From a theoretical point of view, the
development (or improvement) of the thermodynamic background is crucial, as it
is helpful to consistently introduce the necessary couplings, revealing the
functional dependencies and, in parallel, ensuring compatibility with the
second law. From a practical point of view, such advanced models have various
application possibilities, and their presence additionally motivates the
improvement of analytical and numerical solution techniques, and therefore the
referred mathematical analysis is inevitable.

###  5.2. Micro and nanoscale heat conduction.

While the Knudsen number ($\textrm{Kn}=l/L$) serves as a natural connection
between the macroscopic low-temperature and microscopic room-temperature heat
conduction problems, there are some essential differences. First, the
experimental background, e.g., the realization of experiments and how to
measure the temperature, can greatly differ. In this respect, we refer to the
paper of Goodson and Ju [292] for further details. Second, the low-temperature
conditions can enhance couplings and nonlinearities, which are not necessarily
present in a room-temperature environment despite the nanometer length scale.
In contrast, the size dependence of thermal conductivity stands as a central
question since even the steady-state heat transport is influenced by ballistic
contributions [263]. In other words, that means $\lambda=\lambda(\textrm{Kn})$
in which $\lambda\rightarrow\lambda_{0}$ takes the macroscopic bulk limit when
$\textrm{Kn}\ll 1$. However, the other direction ($\textrm{Kn}\sim 1$) is
still unclear, and there are multiple approaches in this respect. First,
according to [293], thermal conductivity cannot be defined when only ballistic
propagation is present in the material, and thus radiative transfer occurs
between the boundaries instead of conduction, standing as an alternative
approach. Consequently, it is worth noting that the following models are valid
only when the diffusive mode is also present, which provides a lower bound for
the length scale. Second, these ballistic boundary scattering contributions
add further thermal resistance to the structure and thus reduce the thermal
conductivity [294, 295]. From a practical point of view, the superlattices,
thin layers, and thermoelectric devices stand as outstanding examples [296,
297]. From an experimental point of view, the so-called time-domain
thermoreflectance (TDTR) method is a widely used standard procedure to detect
the thermal conductivity of such material structure as the reflected optical
signal is sensitive to any temperature change [298, 299]. For a detailed
experimental overview, we refer to the papers of Jiang et al. [300] and Saha
et al. [301, 302]. Since thin layers (and therefore superlattices) are used as
a periodic structure in a microelectronic device, these naturally show
anisotropic properties with strong temperature dependence [303, 304]. Here,
the size dependence is meant to be related to the cross-plane thermal
conductivity, depending on the number of layers (or period thickness) [302,
10, 305]. Figure 10/(A and B) present an example of thin films and
superlattices, and (C) demonstrates the appearance of ballistic effects when
the mean free path changes.

Figure 10. A) Magnified image of a diamond film applied to enhance heat
removal [302]. B) TEM micrograph of a superlattice of Ti0.7W0.3N and
Al0.72Sc0.28N layers [306]. C) The change of effective thermal conductivity of
nanostructures with respect to the ratio of the nanoparticle size ($r$) and
mean free path ($l$) [307].

####  5.2.1. Size-dependent thermal conductivity.

A continuum model does not build the governing equations employing detailed
transport mechanisms, we again begin with the kinetic theory background. It
would be more accurate to use ‘structure-dependent’ thermal conductivity
instead of size dependence since it is not only the size that matters. For
instance, the interface roughness, the material type, and the temperature
conditions influence the effective thermal conductivity. Also, it is apparent
now that Fourier’s law using Debye’s definition for thermal conductivity (8)
cannot be valid at such scales, and one must utilize a more general treatment.
This is more visible from the general definition of the bulk thermal
conductivity from kinetic theory [308], called the dispersion model,

$\displaystyle\lambda=\frac{1}{3}\sum_{p}\int
c(\omega)_{p}v(\omega)_{p}l(\omega)_{p}\textrm{d}\omega,$ (115)

in which the summation denotes different phonon branches, such as optical and
acoustic, and the latter can be further separated into transverse and
longitudinal modes. The integration considers the whole spectrum for each
branch. Furthermore, $v(\omega)$ considers each branch’s characteristic group
velocity, and $c(\omega)$ stands for the frequency-dependent specific heat
capacity; hence, the internal energy is also $\omega$-dependent. Eq. (115) can
be simplified with the so-called gray medium approximation for which the mean
free path $l$ becomes independent of $\omega$. For instance, consider a
superlattice with two components, one can estimate its bulk thermal
conductivity by knowing the component properties

$\displaystyle\lambda=\sum_{i=1}^{2}\chi_{i}\lambda_{i}\big{(}1-1.5p(1-\alpha_{j}/\alpha_{i})A_{1}/\xi_{i}-1.5(1-p)A_{2}/\xi_{i}\big{)},\quad
p\approx\exp(-16\pi^{3}\delta^{2}/\omega^{2}),$ (116)

where $\chi$, $\xi$, and $\delta$ are the relative layer thickness and film
thickness normalized to the mean free path and mean interface surface
roughness, respectively [308]. Here, $p$ stands for the interface scattering
parameter suggested by Ziman [309], ranging from $0$ (diffuse) to $1$
(specular), see Figure 11/A for its influence. For a more flexible treatment,
one can consider $p$ independent of the frequency $\omega$ and used as a
fitting parameter as $p$ can significantly influence how the thermal
conductivity decreases with the layer thickness. Additionally, $A_{1}$ and
$A_{2}$ express the direction and spectral-dependent transmissivity, and their
exact expression can be found in [308].

More recent studies by Saha et al. [302, 310] provides two methods to estimate
the interface thermal resistance. Their starting point is that the overall
thermal resistance of a superlattice is expressed as the sum of each
component, i.e., $R_{\textrm{total}}=R_{\textrm{layer 1}}+R_{\textrm{layer
2}}+R_{\textrm{interface}}$, yielding

$\displaystyle\frac{NR_{\textrm{interface}}}{L}=\frac{1}{\lambda_{\textrm{total}}}-\frac{1}{2}\left(\frac{1}{\lambda_{1}}+\frac{1}{\lambda_{2}}\right),$
(117)

for which the two components have the same thickness, and $L$ denotes the
total thickness of the superlattice in which there are $N$ number of
interfaces. In their first scenario, they suppose that $\lambda_{1}$ and
$\lambda_{2}$ thermal conductivities are the same as the measured $240$ nm
thin film, which implies a diffusive phonon scattering mechanism for
thicknesses larger than the mean free path. The other scenario assumes the
ballistic propagation of phonons; therefore, the thermal conductivities
$\lambda_{1}$ and $\lambda_{2}$ do not contribute to the overall resistance,
which could be reasonable for layer thicknesses being much smaller than the
mean free path. Hence, after measuring the effective total thermal
conductivity $\lambda_{\textrm{total}}$ for multiple situations, one can
estimate how the interfaces contribute to the overall thermal behavior in a
superlattice.

The above-discussed results are further extended with the work of Alvarez et
al. [311, 312], in which the compatibility with kinetic theory is considered
less rigorously but obtained multiple helpful $\lambda(\textrm{Kn})$
expressions depending on the boundary contributions of phonons and
investigated their thermodynamic compatibility in the framework of EIT. They
focus only on longitudinal flow in nanowires with radius $R$. Their starting
point is the GK equation (60) under stationary conditions, viz. there is no
change in time ($\nabla\cdot\mathbf{q}=0$). Moreover, $\mathbf{q}$ is
considered to be much smaller than $l^{2}\Delta\mathbf{q}$ (but
$\mathbf{q}\not\equiv 0$, otherwise the solution would be trivial), thus
solving

$\displaystyle\Delta\mathbf{q}=\frac{\lambda_{0}}{l^{2}}\nabla T.$ (118)

They also assume that the heat flux $\mathbf{q}$ can be divided into a bulk
$\mathbf{q}_{b}$ and a wall contribution part $\mathbf{q}_{w}$, that is,
$\mathbf{q}=\mathbf{q}_{b}+\mathbf{q}_{w}$. They solve Eq. (118) first with
$q_{b}=0$ boundary condition at both ends ($r=\pm R$) in a one-dimensional
setting (along the radius $r$), and supposing that a slip boundary condition
is valid for the wall contribution,

$\displaystyle q_{w}=Cl\left(\frac{\partial q_{b}}{\partial r}\right)_{r=R},$
(119)

the constant $C$ expresses the diffusive and specular boundary scattering
(similarly to the parameter $p$ before). Let us realize that $C=0$ is
excluded, otherwise, the system could not be in a stationary state and would
result in $\mathbf{q}=0$ and $\nabla T=0$. Then, the effective, size-dependent
thermal conductivity is found as

$\displaystyle\lambda(\textrm{Kn})=\frac{1}{\pi R^{2}}\frac{L}{\delta
T}Q_{\textrm{tot}}=\frac{\lambda_{0}}{8\textrm{Kn}^{2}}(1+4C\textrm{Kn}),\quad
Q_{\textrm{tot}}=\int\displaylimits_{0}^{R}2\pi\big{(}q_{b}(r)+q_{w}(r)\big{)}\textrm{d}r,$
(120)

where $\delta T$ expresses the temperature difference, and, more
interestingly, this $\lambda(\textrm{Kn})$ formula cannot provide the bulk
limit $\lambda_{0}$ for very small Knudsen numbers. They test the
thermodynamic compatibility of Eq. (120) by substituting it into the classical
entropy production (5) with local equilibrium assumption for which
$\sigma_{s}\geq 0$ is satisfied. They continue their analysis by extending the
boundary contribution with a second-order term as

$\displaystyle q_{w}=Cl\left(\frac{\partial q_{b}}{\partial
r}\right)_{r=R}-\gamma l^{2}\left(\frac{\partial^{2}q_{b}}{\partial
r^{2}}\right)_{r=R}$ (121)

for which the parameter $\gamma>0$ can be chosen to recover various slip
boundary conditions. Repeating the previous approach,

$\displaystyle\lambda(\textrm{Kn})=\frac{\lambda_{0}}{8\textrm{Kn}^{2}}(1+4C\textrm{Kn}-4\gamma\textrm{Kn}^{2}),$
(122)

is found and investigated later with more detail by Vázquez and Márkus [313,
314, 315]. Here, the situation seems more complicated since, in the bulk
limit, it would lead to negative thermal conductivity, thus, there exists a
critical Kn, which acts as an upper bound. When first testing its
thermodynamic compatibility, using Eq. (5) (the local equilibrium version)
again, they found that such second-order extension for the slip boundary
condition is not thermodynamically compatible, however, exchanging it to a
generalized one with extended state space (for which both $\mathbf{q}$ and
$\mathbf{Q}\sim\nabla\mathbf{q}$ are included as before), then it restores the
compatibility.

At this point, we want to place some necessary criticism regarding the
thermodynamic treatment to be clear and avoid misunderstandings. First, while
it could be reasonable under some special conditions to neglect $\mathbf{q}$,
such a term can significantly modify the solutions of the equation (118) and
put the thermodynamic compatibility into a different view. Second, utilizing a
local equilibrium entropy production for the (simplified) GK equation is not
physically reasonable since the GK equation itself is derived from a weakly
nonlocal framework. Third, it seems unnatural that these
$\lambda(\textrm{Kn})$ expressions have various validity limits, basically
depending on an arbitrary choice of boundary conditions (e.g., how to choose
$\gamma$), not explicitly on the rarefiedness of the medium what an intuition
would suggest. However, despite these flaws, Eq. (122) still can be helpful in
modeling experiments [316], but one has to keep in mind the restrictions
together with all the related assumptions.

Alvarez and Jou [317, 318] also investigated a different approach, exploiting
the hierarchical structure of the EIT evolution equations, having the same
characteristic hierarchical structure as the system (50), with one notable
difference: in EIT, there is a freedom to employ the coefficients from kinetic
theory. Briefly, the EIT system reads,

$\displaystyle
a_{1}\partial_{t}\mathbf{q}^{(1)}-\nabla\frac{1}{T}-b_{1}\nabla\cdot\mathbf{q}^{(2)}$
$\displaystyle=-c_{1}\mathbf{q}^{(1)}$ (123) $\displaystyle
a_{n}\partial_{t}\mathbf{q}^{(n)}-b_{n-1}\nabla\mathbf{q}^{(n-1)}-b_{n}\nabla\cdot\mathbf{q}^{(n+1)}$
$\displaystyle=-c_{n}\mathbf{q}^{(n)},$ (124)

where the coefficients $a_{i}$, $b_{i}$, and $c_{i}$ are now phenomenological,
and the upper indices $(n)$ denote the tensorial order. This is also
characteristic of the internal variable approach [42, 196]. Nevertheless, EIT
requires a hierarchical structure similar to the momentum series expansion to
have a hyperbolic system, but this is not a requirement in the internal
variable approach. Interestingly, when a current multiplier is supposed to
contribute to the entropy flux, it always leads to a parabolic set of
evolution equations. The system (123)-(124) is closed by truncation, i.e., for
the $n$${}^{\textrm{th}}$ variable, the $(n+1)$${}^{\textrm{th}}$ flux is
considered to be zero. This is the most straightforward closure, the framework
of RET is more sophisticated, requiring Galilean covariance and involving the
maximum entropy principle to close the hierarchy [33]. We refer to [319] for a
more comprehensive comparative study.

After taking the Fourier transform of (123)-(124) (interestingly, without
using the energy balance, thus this system is neither mathematically nor
physically closed), they express the thermal conductivity from the resulting
dispersion relation and substitute it into the Fourier law,

$\displaystyle\mathbf{q}(\omega,k)=-ik\lambda(\omega,k)\hat{T}(\omega,k),$
(125)

with $i$ is the imaginary unit, and $\lambda(\omega,k)$ is found to be a
continued-fraction, depending on which level one truncates the hierarchy,

$\displaystyle\lambda(\omega,k)=\frac{\lambda_{0}}{1+i\omega\tau_{1}+\frac{k^{2}l^{2}_{1}}{1+i\omega\tau_{2}+\frac{k^{2}l_{2}^{2}}{1+i\omega\tau_{3}+\dots}}}.$
(126)

Here, the relaxation times $\tau_{i}$ are formed from $a_{i}/c_{i}$, and the
intrinsic length scales $l_{i}$ originate from $b_{i}/(c_{i}c_{i+1})$. While
the bulk term $\lambda_{0}$ can depend on temperature, all the other
coefficients must be constant. The characteristic length of the system $L$ is
introduced through the wave number $k=2\pi/L$; thus, the expression (126)
differs from the usual dispersion relations in this respect. For a stationary
case, they assume $\omega=0$, and the continuation requires further knowledge
on the relaxation times $\tau_{i}$ and the intrinsic length scales $l_{i}$.
Figure 11/B compares its outcome for numerous experimental results [318]. For
further possibilities, we refer to [317, 318], and we add that such size-
dependent thermal conductivity is applied within the transient theory of EIT
[311]. These aspects are not restricted to superlattices and also prevail
concerning nanotubes [320, 8, 321, 9].

####  5.2.2. Short note on nanofluids.

We close this part by clarifying the role of nanofluids as it might be
misleading that, e.g., in [322], nanofluids are mentioned concerning nanoscale
heat conduction. Nanofluids are “ordinary” fluids (such as water) mixed with a
small portion of nanoparticles in order to enhance the fluid thermal
properties and to improve the heat transfer capabilities [92, 323], e.g., in
solar panels [324, 325]. Hence, indeed, the presence of nanoparticles induces
nanoscale phenomena, but usually, these are effectively modeled in a
macroscale approach, for instance, by determining the effective thermal
conductivity for a particular composition, including particle size effects as
well [326]. Nevertheless, the determination of the effective thermal
conductivity of a nanofluid stands as a challenging and still open question,
and there are no general approaches as $\lambda$ might depend on the particle
shape, volume and mass fractions, and sedimentation attributes, see for
instance [327, 328, 329, 330]. Despite the seeming analogy, non-Fourier
effects are neither expected nor observed in nanofluids on macroscale. We
refer to [331] for a detailed overview of nanofluids and their technological
background.

Figure 11. Comparing model predictions to experimental data. A) Phonon model
with adjustable interface scattering parameter [295]. B) EIT continuum model
utilizing Eq. (126) with assuming $\omega=0$ and all $l_{i}$ mean free paths
are equal [318].

####  5.2.3. Transient heat conduction.

Similarly to the steady-state situation, the central question remains about
how the ballistic and diffusive phonon propagation modes contribute to the
overall heat transport. As a particular extremum of these modes, Majumdar et
al. [293, 332] suppose only ballistic modes, and in such case, phonons become
analogous to photons from many aspects, and thus the heat flux becomes

$\displaystyle\mathbf{q}=\sigma_{\textrm{phonon}}(T_{1}^{4}-T_{2}^{4}),$ (127)

between two parallel plates (e.g., in thin films), with
$\sigma_{\textrm{phonon}}$ being the phonon Stefan-Boltzmann constant
depending on the particular phonon modes and their corresponding propagation
speeds. It must be emphasized that the relation (127) is valid only in the
Casimir limit [333] and when the temperature is much below the Debye
temperature [334] (i.e., only the acoustic modes are present). Otherwise, the
heat transfer is not purely ballistic, the phonon-phonon scattering becomes
significant. Due to the analogy with photons, it is suitable to introduce the
phonon intensity $I_{\omega}$ as

$\displaystyle
I_{\omega}(\mathbf{r},\mathbf{v},t)=\frac{\hbar\omega}{4\pi}D(\omega)|\mathbf{v}|f(\mathbf{r},\mathbf{v},t),$
(128)

where $\mathbf{v}$ is the group velocity of heat carriers and $D(\omega)$
describes the phonon density of states per unit volume. In a one-dimensional
case, the Boltzmann transport equation is written as

$\displaystyle\frac{\partial f_{\omega}}{\partial t}+v_{x}\frac{\partial
f_{\omega}}{\partial
x}=\frac{f_{\omega}^{0}-f_{\omega}}{\tau}\Rightarrow\big{(}\textrm{Eq.~{}}\eqref{maj2}\big{)}\Rightarrow\frac{1}{v}\frac{\partial
I_{\omega}}{\partial t}+\mu\frac{\partial I_{\omega}}{\partial
x}=\frac{I_{\omega}^{0}(T(x))-I_{\omega}}{v\tau(\omega,T)},$ (129)

with $v_{x}=\mu v$ with $-1<\mu<1$ being a geometrical factor ($\mu>0$ means
forward, and $\mu<0$ means backward directions), and also, the relaxation time
approximation is applied with the equilibrium intensity corresponding to a
black body at a temperature below the Debye temperature, and following the
Bose-Einstein statistics. For gray bodies, as before, the relaxation time
becomes independent of the frequency $\omega$, and the equilibrium intensity
distribution can be approximated with

$\displaystyle
I_{\omega}^{0}(T)=\frac{1}{2}\int\displaylimits^{1}_{-1}I_{\omega}\textrm{d}\mu,\quad\textrm{and}\quad\sigma_{\textrm{phonon}}T^{4}=\pi\int\displaylimits_{0}^{\omega_{D}}I_{\omega}^{0}(T)\textrm{d}\omega$
(130)

with $\omega_{D}$ being the cut-off frequency related to the Debye
temperature. Equations (129) and (130) together form the so-called ‘Equations
of Phonon Radiative Transfer’ (EPRT) model [332], valid below the Debye
temperature.

Chen’s approach [263], however, explicitly separates the ballistic and
diffusive contributions as
$I_{\omega}(\mathbf{r},t)=I_{\omega,b}(\mathbf{r},t)+I_{\omega,d}(\mathbf{r},t)$,
consequently, $\mathbf{q}=\mathbf{q}_{b}+\mathbf{q}_{d}$ and
$e=e_{b}+e_{d}=CT$, but it is not a two-temperature model since the
temperature $T$ is associated to the internal energy $e$ only, not separately
to $e_{b}$ and $e_{d}$, and $C$ is the specific heat capacity, not assigned to
each propagation mode either. The ballistic contribution
$I_{\omega,b}(\mathbf{r},t)$ is calculated using Eq. (129) with one essential
difference, for ballistic phonons, there is no associated equilibrium state.
For the diffusive part, $I_{\omega,d}(\mathbf{r},t)$, however, the equilibrium
$I_{\omega}^{0}(T)$ remains meaningful, and the time evolution is again given
in the form of Eq. (129). It significantly eases the solution procedure to
apply the diffusion approximation for thermal radiation, being valid only for
optically thick media for which the radiative heat flux is obtained in the
form of Fourier’s law. The model is further simplified by considering the
spherical harmonic expansion of $I_{\omega,d}(\mathbf{r},t)$ as
$I_{\omega,d}(\mathbf{r},t)=J_{\omega,0}(\mathbf{r},t)+\mathbf{n}\cdot\mathbf{J}_{\omega,1}(\mathbf{r},t)$,
which breaks down the original integral-partial differential evolution
equation into a set of coupled partial differential equations, similar to the
momentum expansion. However, here the orthogonality of spherical harmonics can
be exploited. The coupled system reads

$\displaystyle\frac{4\pi}{|\mathbf{v}|}\frac{\partial J_{\omega,0}}{\partial
t}+\frac{4\pi}{3}\nabla\cdot\mathbf{J}_{\omega,1}$
$\displaystyle=-\frac{4\pi}{|\mathbf{v}|\tau_{\omega}}(J_{\omega,0}-I_{\omega}^{0}),$
(131)
$\displaystyle\frac{1}{|\mathbf{v}|}\frac{\partial\mathbf{J}_{\omega,1}}{\partial
t}+\nabla J_{\omega,0}$
$\displaystyle=-\frac{1}{|\mathbf{v}|\tau_{\omega}}\mathbf{J}_{\omega,1},$
(132)

where Eq. (131) is of balance type for the (diffusive) internal energy, and
Eq. (132) recalls the MCV-type constitutive equation, indeed,

$\displaystyle
u_{d}=\int\frac{4\pi}{|\mathbf{v}|}J_{\omega,0}\textrm{d}\omega,\quad
u_{b}=\int\int\frac{1}{|\mathbf{v}|}I_{\omega,b}\textrm{d}\hat{\boldsymbol{\Omega}}\textrm{d}\omega,\quad\mathbf{q}_{d}=\frac{4\pi}{3}\int\mathbf{J}_{\omega,1}\textrm{d}\omega,\quad\mathbf{q}_{b}=\int\int
I_{\omega,b}\cos\theta\textrm{d}\hat{\boldsymbol{\Omega}}\textrm{d}\omega,$
(133)

where the integration of $\textrm{d}\hat{\boldsymbol{\Omega}}$ expresses the
summation over the entire solid angle. Overall, these equations yield

$\displaystyle\frac{\partial e_{d}}{\partial t}+\frac{\partial e_{b}}{\partial
t}$ $\displaystyle=-\nabla\cdot\mathbf{q}_{d}-\nabla\cdot\mathbf{q}_{b},$
(134) $\displaystyle\tau\frac{\partial e_{b}}{\partial
t}+\nabla\cdot\mathbf{q}_{b}$ $\displaystyle=-e_{b},$ (135)
$\displaystyle\tau\frac{\partial\mathbf{q}_{d}}{\partial t}+\mathbf{q}_{d}$
$\displaystyle=-\frac{\lambda_{d}}{C}\nabla e_{d},$ (136)

and the ballistic heat flux contribution $\mathbf{q}_{b}$ can directly be
expressed using Eqs. (129) and (133). This model is called the ‘ballistic-
diffusive’ equation. The balance of $e_{d}$ can be obtained by subtracting
Eqs. (135) from (134), this is only a matter of choice, but only one of them
can be used to avoid the over-determination of the model.

However, according to Lebon et al., one could approximate $\mathbf{q}_{b}$
with an MCV-type equation, too, with different relaxation time and thermal
conductivity, and thus the compatibility with EIT would be easily accessible
[322, 335]. Therefore, this is the moment when we also introduce their EIT
approach. This is purely a macroscopic approach for which it remains arbitrary
whether to imply any assumptions from kinetic theory. The EIT model also
applies the separation of ballistic and diffusive contributions such as
$e=e_{b}+e_{d}$ and $\mathbf{q}=\mathbf{q}_{b}+\mathbf{q}_{d}$, and renders a
sort of ‘quasi-temperature’ $T_{b}$ and $T_{d}$ to the internal energies with
heat capacities $c_{b}$ and $c_{d}$, choosing $c_{d}=c_{b}$, therefore
$T=T_{b}+T_{d}$. We must keep in mind the intensive attribute of $T$, thus
$T_{b}$ and $T_{d}$ are not real thermodynamics quantities, but both field
variables are driven by their balance equations,

$\displaystyle\frac{\partial e_{i}}{\partial
t}+\nabla\cdot\mathbf{q}_{i}=r_{j},\quad\textrm{with}\quad i,j=\\{b,\ d,|\
i\neq j\\}$ (137)

and the total internal energy $e$ fulfills the first law of thermodynamics
with a source term $r=r_{b}+r_{d}$, and using Chen’s model,
$r_{b}=-e_{b}/\tau_{b}$, and $r_{d}$ is constrained through the given $r$.
Then, based on the earlier results from kinetic theory, it is supposed that
the diffusive part of phonons can properly be modeled with the MCV equation,
and the ballistic part, however, is assumed to be driven by the GK equation,

$\displaystyle\tau_{d}\frac{\partial\mathbf{q}_{d}}{\partial
t}+\mathbf{q}_{d}$ $\displaystyle=-\lambda_{d}\nabla T_{d},$ (138)
$\displaystyle\tau_{b}\frac{\partial\mathbf{q}_{b}}{\partial
t}+\mathbf{q}_{b}$ $\displaystyle=-\lambda_{b}\nabla
T_{b}+l_{b}^{2}(\Delta\mathbf{q}_{b}+2\nabla\nabla\cdot\mathbf{q}_{b}),$ (139)

with the coefficients supposed to follow the kinetic theory,

$\displaystyle\lambda_{d}=\frac{1}{3}c_{d}v_{d}^{2}\tau_{d},\quad\lambda_{b}=\frac{1}{3}c_{b}v_{b}^{2}\tau_{b},\quad
v_{d}=\frac{l_{d}}{\tau_{d}},\quad v_{b}=\frac{l_{b}}{\tau_{b}}.$ (140)

At this point, we must note some further comments about the system
(138)-(140). First, sadly, this system is not derived in the same way as
before, i.e., defining the elements of the state space and how the entropy
flux constitutes the heat fluxes. It might be necessary since the Onsagerian
solution of the entropy production would naturally offer a coupling between
$\mathbf{q}_{b}$ and $\mathbf{q}_{d}$, and these parts are omitted in
(138)-(140), although could have outstanding importance even in more
straightforward situations. Second, although [322] states that the system
(138)-(140) provides a macroscopic description, it is not entirely true since
that approach is only applicable for nanosystems considering the size
dependence of $\lambda$, together with the kinetic assumptions. Despite these
hidden attributes, the EIT model is a valuable approach [153], however, with
further undiscovered potential, which would be helpful to recognize the limits
of continuum approaches and discover the modeling capabilities. For a more
detailed overview of the microscale heat transfer and further coupled heat
equations, let us refer to [336, 11].

###  5.3. Comments on the quantum concepts

So far, we reviewed one branch of quantum transport focusing on the phonon-
phonon interaction without introducing (and exploiting) the concepts of wave
functions, quantum statistics, and density matrices. Moreover, the phonon
transport models have counterparts from continuum frameworks, recalling the
exact structure of evolution equations. This does not hold for electron
transport, and even properly formulating Fourier’s law based on quantum
mechanics has been an old problem, especially when the heat flux needs further
generalizations. While the utilization of phonons is also a valid quantum
approach that has been discussed so far, a more general framework should
include, e.g., the electrons and electron-phonon interactions. However, these
particles possess a notably different quantum nature, thus influencing the
corresponding set of wave functions and their symmetry attributes, summarized
in [178]. The introduction of probability distributions is inevitable. The
existence of fluctuations stands as further challenging properties of quantum
systems and might result in negative entropy production [337].

Let us imagine a stationary problem with two heat reservoirs characterized by
the temperatures $T_{1}$ and $T_{2}\neq T_{1}$, and connected via a conductor.
The electron transport simultaneously induces charge ($I_{c}$) and heat
currents ($I_{h}$) between the reservoirs. Instead of the temperature
difference, the resulting currents are proportional to the difference between
the distribution functions ($\vartheta_{1,2}$), following [338],

$\displaystyle I_{c}$
$\displaystyle=\tilde{q}\sum_{p}\int\displaylimits_{0}^{\infty}\frac{1}{2\pi}\mathbf{v}_{p}(\mathbf{k})(\vartheta_{1}(\mathbf{k})-\vartheta_{2}(\mathbf{k}))\mathcal{T}_{p}(\mathbf{k})\textrm{d}\mathbf{k},\quad$
$\displaystyle\Rightarrow\quad$ $\displaystyle
I_{c}=\frac{\tilde{q}}{h}\sum_{p}\int\displaylimits_{\epsilon(0)}^{\infty}(\vartheta_{1}(\epsilon)-\vartheta_{2}(\epsilon))\mathcal{T}_{p}(\epsilon)\textrm{d}\epsilon$
(141) $\displaystyle I_{h}$
$\displaystyle=\sum_{p}\int\displaylimits_{0}^{\infty}\frac{1}{2\pi}\epsilon_{p}(\mathbf{k})\mathbf{v}_{p}(\mathbf{k})(\vartheta_{1}(\mathbf{k})-\vartheta_{2}(\mathbf{k}))\mathcal{T}_{p}(\mathbf{k})\textrm{d}\mathbf{k},\quad$
$\displaystyle\Rightarrow\quad$ $\displaystyle
I_{h}=\frac{1}{h}\sum_{p}\int\displaylimits_{\epsilon(0)}^{\infty}\epsilon(\vartheta_{1}(\epsilon)-\vartheta_{2}(\epsilon))\mathcal{T}_{p}(\epsilon)\textrm{d}\epsilon$
(142)

where $\tilde{q}$ is the particle charge, and
$\epsilon=\epsilon_{p}(\mathbf{k}=\mathbf{0})$ indicates the energy of the
given particle as a function of the wave number $\mathbf{k}$, for all
independent modes $p$. Furthermore, $\mathcal{T}_{p}$ stands for the particle
transmission probability, for ballistic transport
$\mathcal{T}_{p}(\mathbf{k})\equiv 1$. Eqs. (141) and (142) are applicable for
both bosons and fermions with their corresponding distribution function and
charge, thus, the resulting electron transport can also be characterized.
Assuming a single ideal ballistic channel between the reservoirs leads to the
observation that the thermal conductivity becomes quantized,

$\displaystyle\lambda_{0}=\frac{\pi^{2}k_{B}^{2}}{3h}T,$ (143)

similarly to the electric conductance, being the same for both bosons and
fermions and first experimentally observed by Scwab et al. [339]. For further
experimental and theoretical details, let us refer to [340, 338]. Although
such an approach provides a quantitative estimate for the energy (and charge)
transfer through a given conductor (not necessarily in ballistic mode), being
helpful in quantum electronics design, it is limited in a sense to offer a
field equation to model the transients of the studied domain.

The paper of Greiner et al. [341] provides the closest formalism that we have
used so far for a one-dimensional conductor (valid for, e.g., a nanowire) of
length $\ell$, considering electrons to be heat carriers. The solution for the
entropy inequality is provided in the form of

$\displaystyle\mathbf{j}_{\mu}=\sum_{\nu}L_{\mu\nu}\mathbf{X}_{\nu},\quad\mu,\nu=1,2$
(144)

in which $L_{\mu\nu}$ are the Onsagerian coefficients, $\mathbf{j}_{\mu}$
represent the electrical current and heat flux densities, and
$\mathbf{X}_{\nu}$ express the corresponding gradients. Exploiting the
fluctuation-dissipation theorem, the linear response is found using a
correlation function $C_{\mu\nu}(t)$. In that particular case, it reads

$\displaystyle
C_{\mu\nu}(t)=\frac{1}{\pi}\left(\frac{e\hbar}{mL}\right)^{2}\left(\frac{\hbar^{2}}{2me}\right)^{\mu+\nu-2}e^{-t/\tau_{c}}\int\displaylimits_{-L/2}^{L/2}\textrm{d}z\int\displaylimits_{k^{-}}^{k^{+}}k^{2(\mu+\nu-1)}f_{0}(1-f_{0})\textrm{d}k,\quad
k^{\pm}=\frac{m}{\hbar t}\left(\pm\frac{L}{2}-z\right),$ (145)

where $e$ is the electron charge, $m$ is the carrier’s effective mass, $\hbar$
stands for the reduced Planck’s constant, $k$ is the wave vector, and $f_{0}$
is the Fermi distribution. The exponential term $e^{-t/\tau_{c}}$ primarily
determines the transport regime through the average scattering time
(scattering on the lattice imperfections). That is, the transport is called
ballistic if $e^{-t/\tau_{c}}=1$, implying that $\tau_{c}\rightarrow\infty$,
and also providing a continuous transition towards the diffusive transport
when $C_{\mu\nu}(t)\sim e^{-t/\tau_{c}}$. The system’s response is determined
on the spectral properties of $C_{\mu\nu}(t)$; therefore, they apply a
Fourier-Laplace transform on Eq. (145), finding $S_{\mu\nu}(\omega)$, hence
the conduction coefficients $L_{\mu\nu}$ and the thermal conductivity
$\lambda$ are

$\displaystyle
L_{\mu\nu}(\omega)=\frac{LS_{\mu\nu}(\omega)}{k_{\textrm{B}}T}\quad\Rightarrow\quad\lambda(\omega)=\frac{L_{11}(\omega)L_{22}(\omega)-L_{12}(\omega)L_{21}(\omega)}{L_{11}(\omega)T}.$
(146)

That approach naturally includes the size dependence in $\lambda$. Moreover,
it also returns the universal thermal conductance (reciprocal of the thermal
resistance) of $\lambda_{0}=\pi^{2}k_{\textrm{B}}^{2}T/(3h)$ being valid for
perfect (or ideal) channels. That ideal thermal conductance can be transferred
to non-ideal situations [342].

Interestingly, Michel et al. [343, 344] simulated heat conduction in a one-
dimensional chain, consisting of $N$ identical and weakly coupled subsystems,
analogously to the Fermi-Pasta-Ulam problem [345]. The system is described by
the Hamiltonian

$\displaystyle
H=\sum_{\mu=1}^{N}H_{\textrm{loc}}(\mu)+g\sum_{\mu=1}^{N-1}V(\mu,\mu+1),$
(147)

in which the weakness of the coupling is expressed by $g\ll\Delta E$, i.e.,
the interaction $V(\mu,\mu+1)$ between each subsystem is re-scaled with a
parameter $g$, being much smaller than the local energy gap $\Delta E$ between
the subsystems. Each subsystem has $n$ numbers of equally distributed levels
plus their ground level. Additionally, the local Hamiltonian describes the
local energy using $\Psi(t)H_{\textrm{loc}}\Psi(t)=\Delta EP_{\mu}$, where
$\Psi(t)$ is a wave function and $P_{\mu}$ is the probability of finding the
$\mu^{\textrm{th}}$ subsystem in an excited state. With introducing a time
evolution operator $D(\hat{\tau})$ in a way
$\Psi(t+\hat{\tau})=D(\hat{\tau})\Psi(t)$ (the notation $\hat{\tau}$ is
introduced for the transition time to distinguish it from the relaxation time
$\tau$). Consequently, the $P_{\mu}$ probabilities inherit that time
dependence and thus form a dynamics, a set of ordinary differential equations
for all $P_{\mu}$ for small values of $\hat{\tau}$. It will reduce to the
Fourier heat equation only if the conditions

$\displaystyle 2g\frac{n}{\delta\epsilon}\geq 1,\quad
g^{2}\frac{n}{\delta\epsilon^{2}}\ll 1$ (148)

are satisfied, and $\delta\epsilon\ll\Delta E$ is the bandwidth.

Manzano et al. [346] follow an analogous procedure, imaging a $N$-element
chain of a two-level quantum system with the Hamiltonian

$\displaystyle
H=\sum_{j=1}^{N}\frac{\hbar\omega}{2}\sigma_{j}^{z}+\sum_{j=1}^{N-1}\hbar
g(\sigma_{j}^{+}\sigma_{j+1}^{-}+\sigma_{j}^{-}\sigma_{j+1}^{+}),$ (149)

where $\sigma_{j}^{z}$, $\sigma_{j}^{+}$, and $\sigma_{j}^{-}$ are Pauli’s
$z$, raising and lowering operator, and $g$ is again a coupling constant. The
master equation of Lindblad form describes the time evolution of the quantum
system,

$\displaystyle\dot{\rho}=-\frac{i}{\hbar}[H,\rho]+\mathcal{L}_{1}\rho+\mathcal{L}_{2},$
(150)

in which $\rho$ denotes the density matrix, and $\mathcal{L}_{1,2}$ models the
interactions with the bosonic heat reservoirs. The resulting heat current can
be found using $\textrm{Tr}(H\dot{\rho})$. However, one must add a so-called
dephasing term to Eq. (150) to obtain a heat flux compatible with Fourier’s
law. The dephasing term introduces a coherence decay into the system with

$\displaystyle\mathcal{L}_{\textrm{deph}}\rho=\gamma\sum_{j=1}^{N}\left(\sigma_{j}^{+}\sigma_{j}^{-}\rho\sigma_{j}^{+}\sigma_{j}^{-}-\frac{1}{2}\\{\sigma_{j}^{+}\sigma_{j}^{-},\rho\\}\right),$
(151)

with an adjustable parameter $\gamma$ to investigate the sensitivity of the
resulting heat current to dephasing. Although the two-level chain itself
strictly restricts the heat capacity, and thus the dependence on the
reservoir’s temperature, a proper size-dependent behavior can be modeled
[346].

These approaches can reconstruct the diffusive energy transport, but the non-
diffusive modes are missing, and one can define multiple Hamiltonians to
create quantum chains with diffusive behavior and proper scaling. From that
point of view, the work of Mazza et al. [347] could be a notable progress.
They investigate the dynamics of a so-called layered single-band Hubbard
model, considering the sub-picosecond time scale for electronic interactions.
The topmost layer is excited with an ultra-fast light pulse, creating hot
electrons. The interaction between the cold and hot electrons introduces
different time scales, which are implemented into the time evolution of the
density matrix. The results include the presence of all propagation modes. The
heat flux shows a strong correlation with the density of hot electrons. During
ballistic propagation, the temperature gradient is negligible. It starts to be
significant in the hydrodynamic regime, and lastly, in the diffusive regime,
the temperature gradient and heat flux are synchronized. These relationships
between $\nabla T$ and $q$ are often disregarded and would be essential to
study in the future for different systems, whether these are more general or
not. For further details, see Figure 12 and [347].

Figure 12. The time evolution of heat flux and temperature gradient,
presenting the characteristics for each propagation mode [347].

We find it insightful to mention the recent work of Ván [348], about
investigating the relations between capillary fluids, superfluids, and quantum
mechanics. Interestingly, these seemingly different models (with different
physical meanings) are connected due to a strong common property called
classical holographic property, meaning that the pressure of perfect fluids
($\mathbf{P}$) is connected to a scalar potential ($\phi$) via the relation

$\displaystyle\nabla\cdot\mathbf{P}=\rho\nabla\phi.$ (152)

This happens to be a thermodynamic consequence under zero dissipation and thus
can be helpful to find a quantized form of the heat conduction field equation
directly without assuming particular chains and interactions among them.

####  5.3.1. Thermal rectification

Let us imagine a situation in which we have a one-dimensional conductor
between two thermal reservoirs. In a steady-state, $Q_{1\rightarrow 2}$
current occurs. Now let us exchange the temperature of reservoirs, i.e.,
reversing the direction of the temperature gradient. The resulting
$Q_{2\rightarrow 1}$ current is not necessarily the same if certain conditions
are satisfied. This thermal rectification can be characterized in various
ways, such as the ratio of these heat currents. Although that phenomenon does
not demand either the presence of quantum effects or particular length, time,
or energy scales, it is still strongly connected to quantum systems due to its
practical applications. Thermal rectification can be exploited in electronic
systems as thermal diodes (on the analogy of electronic diodes).

Asymmetry stands as the most critical requirement to achieve such a thermal
rectification effect. A straightforward macroscopic example is a functionally
graded material in which the material properties depend on space, consequently
also possessing different thermal conductivity. Moreover, their temperature
dependence is also different, hence nonlinearities can induce asymmetry as
well [349]. The additional, macroscopically feasible cause is the presence of
distinct interface attributes. The different surface and mechanical properties
can also induce thermal rectification, which is called thermal strain or
warping [350].

However, such an effect can be realized in practice in particular
heterogeneous structures (recall the aim of thermal metamaterials), for which
the material structure with heterogeneities is manufactured in an intended
way. For very insightful examples, we refer to the paper of Zhao et al. [351],
in which paper they collected and reviewed a vast number of experimental works
with discussing the possible causes. Interestingly, shape is also a factor,
Yang et al. [352] designed and then experimentally tested several trapezoid-
shaped microelectrodes.

From a quantum mechanical point of view, thermal rectification can be realized
with a graded system suited to a proper long-range interaction [353, 354]. The
corresponding Hamiltonian reads

$\displaystyle
H=\sum_{j=1}^{N}\left(\frac{p_{j}^{2}}{2m_{j}}+\frac{u_{j}^{4}}{4}\right)+\sum_{j,k}\frac{(u_{j}-u_{k})^{2}}{2(1+|j-k|^{\beta})},$
(153)

where $p_{j}$ denotes the momentum of the j${}^{\textrm{th}}$ particle with
mass $m_{j}$, and $u_{j}$ is its displacement. The graded system means that
the mass distribution of the particles satisfies $m_{j}>m_{j-1}$,
nevertheless, at the ends $m_{0}=m_{1}$ and $m_{N}=m_{N-1}$. The last term
represents the interaction potential, where $\beta$ characterizes the
polynomial decay. Both the mass asymmetry and the long-range interaction are
in favor of the rectification phenomenon and remain significant even for long
chains $N>200$. We stress that asymmetry remains crucial, e.g., a chain
consisting of equal masses can still present rectification if an additional,
space-dependent external field can act on the system [355]. Let us recall that
Eq. (80) is motivated by the same reasons, and a particular nonlocal (on the
analogy of long-range interaction) is implemented into a continuum framework,
leading to further modifications of the Guyer-Krumhansl equation [242]. We
refer to [354, 356, 357] for further quantum mechanical aspects and insight.
Thermal transistors are constructed analogously, in that regard, we only refer
to [358, 359] as that topic points beyond the limits of the present paper.

###  5.4. Rarefied gases.

It is worth starting this topic again with kinetic theory as there is a clear
analogy between rarefied phonon gas and rarefied real gaseous materials, the
Knudsen number is similarly large in both cases, but the carrier is different.
The modeling of that difference over the energy carriers distinguishes between
the models based on kinetic theory. For instance, it does matter whether the
gas is monatomic or polyatomic and hence how the source term in the Boltzmann
equation is formulated. Such details are missing from the continuum models,
but they appear again to be more flexible, especially when nonlinearities
enter the picture. From an experimental point of view, the rarefied state is
achieved by decreasing the system’s pressure while keeping its volume
constant. There should be a transition from the classical Navier-Stokes-
Fourier (NSF) equations to the more advanced, generalized systems by having
the mass density dependence of the transport coefficients such as viscosities,
thermal conductivity, and, therefore, the new coefficients (relaxation times
and various coupling parameters). Such transition is observed as the change in
the speed of sound, starting from the dense (normal) state to the rarefied
region. However, modeling that transition is not straightforward, and thus
various models can be derived on kinetic and continuum bases. To be more
precise, Figure 13 presents a typical experimental arrangement and data
recorded for the change of speed of sound as a function of frequency over
pressure ($\omega/p$), which scaling property will be discussed soon. The
classical NSF system cannot explain these observations, and these experiments
act as a benchmark for extended theories.

Figure 13. Rarefied gases from an experimental point of view. A) Typical
measurement arrangement [360]. B) Change in the speed of sound with respect to
mass density [361]. C) Mass density dependence of shear viscosity towards the
dense states [362, 363]. D) Mass density dependence of shear viscosity towards
the rarefied states [364] (different from part C)).

####  5.4.1. Kinetic approach.

The starting point is a suitable approximation for the collision integral in
the Boltzmann transport equation (49), such as the Bhatnagar-Gross-Krook (BGK)
or ellipsoidal-statistical (ES)-BGK model [365], in the form of
$S=-\nu(f-f_{\textrm{eq}})$, where $\nu$ is the collision frequency ($\sim
1/\tau$), and $f_{\textrm{eq}}$ expresses the corresponding equilibrium
distribution for a monatomic gas. These models differ in how the collision
frequency and the equilibrium distribution are utilized (e.g., on what
quantities they depend on), what molecule type with what potential is included
[366, 267]. Moreover, it is also a matter of choice at which level we want to
approximate the Boltzmann equation together with its collision integral.
Consequently, here we restrict ourselves to the conceptual presentation only
due to the large number of modeling possibilities. Since the transport
coefficients are given as an outcome, these approximations also restrict their
ratio (like the Prandtl number) and hence how realistic the model is for a
particular medium. For instance, for the BGK model, the Prandtl number is 1,
which could be quite inaccurate, but in parallel, it offers a simpler approach
than the ES-BGK approximation. However, the latter model allows for adjustment
of the Prandtl number through the equilibrium distribution and also consists
of the BGK model as a particular choice.

Both the kinetic and continuum approaches require the basic conservation laws
for mass, momentum, and energy. In the kinetic theory, these balances are
found through the collision invariants $\psi=\psi(\mathbf{r},\mathbf{v},t)$,
such as

$\displaystyle\psi=1:\ \rho=m\int
f_{\textrm{eq}}\textrm{d}\mathbf{c},\quad\psi=\mathbf{c}:\
\rho\mathbf{v}=m\int\mathbf{c}f_{\textrm{eq}}\textrm{d}\mathbf{c},\quad\psi=\frac{1}{2}c^{2}:\
\rho e=\frac{3}{2}\rho
T=\frac{m}{2}\int(\mathbf{c}-\mathbf{v})^{2}f_{\textrm{eq}}\textrm{d}\mathbf{c},$
(154)

for which $\mathbf{v}$ stands as the barycentric velocity and $\mathbf{c}$ is
the phase velocity, and the corresponding balances are

mass: $\displaystyle\dot{\rho}+\rho\nabla\cdot\mathbf{v}$ $\displaystyle=0,$
(155) momentum: $\displaystyle\rho\dot{\mathbf{v}}+\mathbf{P}\cdot\nabla$
$\displaystyle=\rho\mathbf{f},$ (156) total energy:
$\displaystyle\rho\left(e+\frac{1}{2}v^{2}\right)\accentset{\mbox{.}}{\vphantom{\left(e+\frac{1}{2}v^{2}\right)}}+(\mathbf{v}\cdot\mathbf{P}+\mathbf{q})\cdot\nabla$
$\displaystyle=\rho\mathbf{f}\cdot\mathbf{v},$ (157) internal energy:
$\displaystyle\rho\dot{e}+\nabla\cdot\mathbf{q}$
$\displaystyle=-\mathbf{P}:\nabla\mathbf{v},$ (158)

with $\mathbf{P}$ being the pressure tensor with a suitable decomposition of
$\mathbf{P}=p\mathbf{I}+\boldsymbol{\Pi}$ in which $p$ is the hydrostatic
pressure (not identical with the spherical part), and
$\boldsymbol{\Pi}=\boldsymbol{\Pi}^{\textrm{dev}}+\Pi\ \mathbf{I}$ is the
viscous pressure, where $\boldsymbol{\Pi}^{\textrm{dev}}$ is the traceless
deviatoric part, and $\Pi$ is called dynamic pressure. Furthermore,
$\mathbf{f}$ is a force density (such as gravitation). These balances are the
starting point in continuum theories. In a kinetic approach, these follow from
the Boltzmann equation. Furthermore, any approximations of the collision
integral $S$ must satisfy these balances as well as the $h$-theorem. Here, we
note that similarly to the balance equations (155)-(158), the entropy balance
can also be derived with a collision invariant $\psi=-k_{\textrm{B}}\ln(f/y)$
with $k_{\textrm{B}}$ being the Boltzmann constant, resulting in the kinetic
counterpart of Eq. (1), using the same notations as Eq. (1),

$\displaystyle s=-k_{\textrm{B}}\int
f\ln\frac{f}{y}\textrm{d}\mathbf{c},\quad\mathbf{J}_{s}=-k_{\textrm{B}}\int\mathbf{c}f\ln\frac{f}{y}\textrm{d}\mathbf{c},\quad\sigma_{s}=-k_{\textrm{B}}\int
Sf\ln\frac{f}{y}\textrm{d}\mathbf{c}.$ (159)

Again, we emphasize that it becomes a matter of choice on what level we want
to approximate the Boltzmann equation, viz., how to close these set of
balances. The calculations and the resulting models can greatly vary in this
respect, and therefore we refer to these works for a more detailed overview
[366, 267].

Conceptually, the aim is to find a suitable and manageable approximation of
the Boltzmann equation, with a given set of potentials, and molecule types,
for a given state space, and the calculability strongly depends on the chosen
setting. Each procedure results in a unique distribution function $f$, which
is used to close the hierarchical system with constitutive equations and
construct the corresponding entropy function. We emphasize again that due to
the numerous possibilities and technical details, we keep focusing on the
conceptual questions.

Let us begin with the Chapman-Enskog (CE) expansion for which the main idea is
to expand the phase density $f$ into a series $f=f^{(0)}+\varepsilon
f^{(1)}+\varepsilon^{2}f^{(2)}+\dots$, where $\varepsilon$ stands for a small
parameter. In the CE expansion, this is identified with the Knudsen number,
however, this is not necessary and could depend on the particular molecular
model. Furthermore, any expansion must satisfy the compatibility conditions,
i.e., when $f_{\textrm{eq}}$ substituted with $f^{(0)}$ in Eq. (154), hence
$f^{(0)}$ is found as the local Maxwellian distribution; and these are
supplemented with

$\displaystyle 0=m\int
f^{(\alpha)}\textrm{d}\mathbf{c},\quad\mathbf{0}=m\int\mathbf{c}f^{(\alpha)}\textrm{d}\mathbf{c},\quad
0=\frac{m}{2}\int(\mathbf{c}-\mathbf{v})^{2}f^{(\alpha)}\textrm{d}\mathbf{c},\quad\alpha\geq
1.$ (160)

The CE expansion includes five fields ($\rho$, $\mathbf{v}$ and $T$), and
closes the balances with

$\displaystyle q^{(\alpha)}_{i}=\frac{m}{2}\int
c^{2}c_{i}f^{(\alpha)}\textrm{d}\mathbf{c},\quad\sigma^{(\alpha)}_{ij}=m\int
c_{\langle i}c_{j\rangle}f^{(\alpha)}\textrm{d}\mathbf{c},$ (161)

which can be used to construct the series of heat flux $q_{i}=\varepsilon
q_{i}^{(1)}+\varepsilon^{2}q_{i}^{(2)}+\dots$, and pressure
$\sigma_{ij}=\varepsilon\sigma_{ij}^{(1)}+\varepsilon^{2}\sigma_{ij}^{(2)}+\dots$
up to an arbitrary order, and $\sigma_{ij}$ stands for the symmetric part of
$\boldsymbol{\Pi}^{\textrm{dev}}$, and thus it implies that the dynamic
pressure $\Pi$ is assumed to be zero together with the antisymmetric part of
$\boldsymbol{\Pi}^{\textrm{dev}}$. The notations are following Eq. (50). Since
$f^{(0)}$ is given as a local Maxwellian, hence $P_{ij}^{(0)}=0$, and
$q_{i}^{(0)}=0$, and the equations of Euler fluid is given. In the first
order, the classical NSF equations are found with

$\displaystyle q_{i}^{(1)}=-\frac{5}{2}\frac{p}{\nu}\frac{\partial T}{\partial
x_{i}},\quad\sigma_{ij}^{(1)}=-\frac{2}{1-b}\frac{p}{\nu}\frac{\partial
v_{\langle i}}{\partial
x_{j\rangle}},\quad\lambda=\frac{5}{2}\frac{p}{\nu},\quad\mu=\frac{1}{1-b}\frac{p}{\nu},$
(162)

where the parameter $-1/2\leq b<1$ is used to adjust the Prandtl number (Pr
$\sim\mu/\lambda$) in the ES-BGK approximation (where $\mu$ is shear
viscosity). Higher-order expansions result in the Burnett (second-order) and
super-Burnett (third-order) equations. However, these are found to be unstable
under particular conditions but might be used only for steady-state problems
[367, 368]. Further orders are extremely cumbersome to determine, and other
methods are more suitable. This highlights an issue that the sequential
approximations and series expansions are not necessarily convergent and might
result in nonphysical solutions. Despite these drawbacks, the kinetic approach
could be undoubtedly valuable when the transport coefficients must be
determined and measurements are unavailable, especially for the coupled
thermo-diffusion problems [369, 370], or heat transfer coefficient in a
rarefied medium.

Interestingly, one might find the thermal conductivity and viscosity to depend
only on the temperature, not on the mass density [371, 364, 372, 373]. The
literature can be divisive in that sense as some argue that these quantities
remain constant for a given temperature, but other measurements show its
opposite [374]. Furthermore, the mentioned $\omega/p$ scaling is strongly
related to this property. Such scaling means that by calculating the
dispersion relation of the transport equations, one finds the phase velocity
as a function of $\omega/p$ only, however, that is true only for an ideal gas
with constant transport properties, and the mass density dependence would
violate that scaling. This also holds for higher-order approximations [375,
376], and it has a crucial role in experimental modeling. To be clear, in
numerous experiments [361, 360, 377], only the temperature is kept constant,
and either the $\omega$ or $p$ are varied, sometimes only the pressure is
changed in a considerable interval; thus only the mass density dependence
could cause the deviation from the NSF equations. We may add that such
theoretical result leads to nonzero transport coefficients in the zero density
limit [33] contrary to the experiments [378], we refer to Fig. 13 for the
related experimental background. Kinetic theory is relatively inflexible from
this perspective compared to a continuum approach. Besides, numerous
experiments show the mass density dependence of viscosity, for instance, for
both high and low-pressure intervals [379, 362], which require further
corrections, e.g., as a function of the Knudsen number [372, 380, 381].

Furthermore, it turned out that the CE expansion can lead to negative phase
densities [267]. This property is shared with the so-called Grad’s method
[382, 383], for which the set of variables is extended with the heat flux
$\mathbf{q}$ and deviatoric pressure $\mathbf{\sigma}$ (13-moment equations),
or even further, higher tensorial order moments can be included (26 or more
moment equations). Similarly to the previous theories, such extension of state
space leads to the appearance of time derivatives of these quantities. The
constitutive equations are found in a balance form, generalizing the classical
NSF equations in a particular way. The resulting system of equations respects
the symmetric hyperbolic structure, together with the balances (155)-(158),
the 13-moment equations reads

$\displaystyle\dot{q}_{i}+\frac{5}{2}\rho T\frac{\partial T}{\partial
x_{i}}+\frac{5}{2}\sigma_{ik}\frac{\partial T}{\partial
x_{k}}-\sigma_{ik}T\frac{\partial\ln\rho}{\partial
x_{k}}-\frac{\sigma_{ik}}{\rho}\frac{\partial\sigma_{kl}}{\partial
x_{l}}+T\frac{\partial\sigma_{ik}}{\partial
x_{k}}+\frac{7}{5}q_{i}\frac{\partial v_{k}}{\partial
x_{k}}+\frac{7}{5}q_{k}\frac{\partial v_{i}}{\partial
x_{k}}+\frac{2}{5}q_{k}\frac{\partial v_{k}}{\partial x_{i}}$
$\displaystyle=-\frac{2}{3}\frac{p}{\mu}q_{i},$ (163)
$\displaystyle\dot{\sigma}_{ij}+\frac{4}{5}\frac{\partial q_{\langle
i}}{\partial x_{j\rangle}}+2\sigma_{k\langle i}\frac{\partial
v_{j\rangle}}{\partial x_{k}}+\sigma_{ij}\frac{\partial v_{k}}{\partial
x_{k}}+2\rho T\frac{\partial v_{\langle i}}{\partial x_{j\rangle}}$
$\displaystyle=-\frac{p}{\mu}\sigma_{ij},$ (164)

being valid only for Maxwell molecules (for which the potential $V(r)\sim
r^{-4}$, and that represents a borderline between hard and soft potentials,
too [83]). It is also claimed that the convergence with moments is rather
slow, nonetheless, it still needs to model fast phenomena or large gradients
(such as shocks). It means that the Guyer-Krumhansl-type equations are
excluded, and non-local terms (such as $\Delta\mathbf{q}$) are not allowed.
Furthermore, a critical Mach number exists, depending on the number of moments
for which sub-shocks and discontinuities occur. These can be contradictory
with experiments and molecular-level simulations, and their appearance depends
on the number of moments [267, 384].

Since the phase density can be negative, even with small amplitudes, it is
impossible to construct the corresponding entropy function, which is feasible
only for linearized equations. Moreover, besides the second law of
thermodynamics, even hyperbolicity can be broken [385], and thus this approach
is accompanied by additional conditions on the hyperbolicity radius around
equilibrium [32]. This shortcoming can partially be solved using the so-called
Maximum Entropy (abbreviated as ’MaxEnt’ [386, 387]) closure, being valid only
for ideal gases [388]. Here, we seek an appropriate $f$, which maximizes the
entropy density from Eq. (159), and the resulting phase density is

$\displaystyle
f_{\textrm{max}}=y\exp\big{(}-\sum\displaylimits_{A}\Lambda_{A}\phi_{A}\big{)},\quad
A=\\{1,2,\dots,N\\},$ (165)

where $\Lambda_{A}$ are the Lagrange multipliers for the $N$-moment system
with moments $F_{A}$, and $\phi_{A}$ are even-order polynomials of
$\mathbf{c}$, otherwise, the resulting system will not be thermodynamically
compatible (hence unstable). That requirement excludes the 13-moment equations
but supports the 14-moment system. Additionally, above the 10-moment equations
(13-moment equations without $\mathbf{q}$), this closure becomes feasible only
numerically. The symmetry
$F_{A}(\Lambda_{A})\longleftrightarrow\Lambda_{A}(F_{A})$ shows the symmetric
property of the hyperbolic system, and it becomes visible only when the system
is represented in with the variables $\Lambda_{A}$; therefore the variables
$\Lambda_{A}$ are together called main field [267].

Another approach to obtain a more reliable set of equations is a so-called
regularization procedure, analogous to the current multipliers seen earlier
regarding Eq. (65). In other words, one introduces additional relaxed
variables (i.e., without time dependence). Here, they appear as auxiliary
(relaxed) moments with a much faster time scale than the main field. For the
detailed derivation, see [267]. The resulting model becomes parabolic (such as
the GK equation). In the example of 13-moment equations, the regularization
yields,

$\displaystyle\dot{q}_{i}+\frac{5}{2}\rho T\frac{\partial T}{\partial
x_{i}}+\frac{5}{2}\sigma_{ik}\frac{\partial T}{\partial
x_{k}}-\sigma_{ik}T\frac{\partial\ln\rho}{\partial
x_{k}}-\frac{\sigma_{ik}}{\rho}\frac{\partial\sigma_{kl}}{\partial
x_{l}}+T\frac{\partial\sigma_{ik}}{\partial x_{k}}+$
$\displaystyle+\frac{7}{5}q_{i}\frac{\partial v_{k}}{\partial
x_{k}}+\frac{7}{5}q_{k}\frac{\partial v_{i}}{\partial
x_{k}}+\frac{2}{5}q_{k}\frac{\partial v_{k}}{\partial x_{i}}+$
$\displaystyle\frac{1}{2}\frac{\partial R_{ik}}{\partial
x_{k}}+\frac{1}{6}\frac{\partial\Delta}{\partial x_{i}}+m_{ikl}\frac{\partial
v_{k}}{\partial x_{l}}=-\frac{2}{3}\frac{p}{\mu}q_{i},$ (166)
$\displaystyle\dot{\sigma}_{ij}+\frac{4}{5}\frac{\partial q_{\langle
i}}{\partial x_{\rangle j}}+2\sigma_{k\langle i}\frac{\partial
v_{j\rangle}}{\partial x_{k}}+\sigma_{ij}\frac{\partial v_{k}}{\partial
x_{k}}+$ $\displaystyle 2\rho T\frac{\partial v_{\langle i}}{\partial
x_{j\rangle}}+\frac{\partial m_{ijk}}{\partial
x_{k}}=-\frac{p}{\mu}\sigma_{ij},$ (167) $\displaystyle
m_{ijk}=\frac{4}{3}\frac{\sigma_{\langle ij}q_{k\rangle}}{p}-2\mu
T\frac{\partial}{\partial x_{\langle
i}}\left(\frac{\sigma_{jk\rangle}}{p}\right),\quad
R_{ij}=\frac{20}{7}\frac{\sigma_{k\langle i}\sigma_{j\rangle
k}}{\rho}+\frac{64}{25}\frac{q_{\langle i}q_{j\rangle}}{p}-\frac{24}{5}\mu
T\frac{\partial}{\partial x_{\langle i}}\left(\frac{q_{j\rangle}}{p}\right),$
$\displaystyle\Delta=5\frac{\sigma_{kl}\sigma_{kl}}{\rho}+\frac{56}{5}\frac{q_{k}q_{k}}{p}-12\mu
T\frac{\partial}{\partial x_{k}}\left(\frac{q_{k}}{p}\right),$ (168)

for which we have kept the original notations for the additional moments
following for the sake of comparability and valid only for Maxwell molecules.
The set of equations (166)-(168) are denoted with R13, and taking the
additional moments ($m_{ijk}$, $R_{ij}$, and $\Delta$) to be zero, it reduces
to Eqs. (163)-(164). These relaxed variables appear as a combination of the
lower-order moments contrary to the current multipliers, and the coefficients
strongly depend on the molecule model. For a more exhaustive systematic
overview of the kinetic models, we refer to Struchtrup [267, 389, 367] and
Cercignani [366]. Furthermore, it is also possible to derive the 13-moment
equations in the framework of GENERIC, providing further insights into the
proper thermodynamic formulation and, in parallel, avoiding the original
shortcomings (e.g., loss of hyperbolicity) [390]. An alternative approach to
stabilize the momentum equations and find a stable version of the Burnett
equation is presented in [391].

For monatomic gases, Rational Extended Thermodynamics inherits the
hierarchical structure of moment expansion, therefore, it becomes possible to
obtain identical evolution equations. However, the difference occurs in the
closure as RET does not attempt to rely on any particular distribution
function or its approximations but is based on Galilean covariance and entropy
principle with Liu’s procedure. Concerning Grad’s 13-moment model, Eqs.
(163)-(164), these approaches can result in a completely identical outcome
[33].

The situation changes for polyatomic gases due to the additional degrees of
freedom for which the existing single-hierarchy ($F$-series) is replaced with
a double-hierarchy ($F$ and $G$-series),

$\displaystyle\partial_{t}F+\partial_{k}F_{k}$ $\displaystyle=0,$
$\displaystyle\partial_{t}F_{i}+\partial_{k}F_{ik}$ $\displaystyle=0,$
$\displaystyle\partial_{t}F_{ij}+\partial_{k}F_{ijk}$ $\displaystyle=S_{ij},$
$\displaystyle\quad\partial_{t}G_{ll}+\partial_{k}G_{llk}$ $\displaystyle=0,$
$\displaystyle\quad\partial_{t}G_{lli}+\partial_{k}G_{llik}$
$\displaystyle=Q_{lli},$ (169)

where the densities can be expressed in terms of physical variables as
follows:

$\displaystyle F$ $\displaystyle=\rho,$ $\displaystyle F_{i}$
$\displaystyle=\rho v_{i},$ $\displaystyle F_{ij}$ $\displaystyle=\rho
v_{i}v_{j}+(p+\Pi)\delta_{ij}+\Pi_{\langle ij\rangle},$ $\displaystyle\quad
G_{ll}$ $\displaystyle=\rho v^{2}+2\rho e,$ $\displaystyle\quad G_{lli}$
$\displaystyle=\rho v^{2}v_{i}+2(\rho e+p+\Pi)v_{i}+2\Pi_{\langle
ki\rangle}v_{k}+2q_{i},$ (170)

thus the $G$-hierarchy becomes independent implying that $G_{ll}\neq F_{ll}$,
viz., $2\rho e\neq 3p$, moreover, these elements of the $G$-hierarchy
represent Galilean transformation rules for energy and heat flux [71].
Consequently, that structure also enables the treatment of the dynamic
pressure $\Pi$ as a new field variable. We note that the size of the
$G$-hierarchy must always be smaller by $1$ than the $F$-hierarchy to respect
Galilean covariance. It has particular importance in modeling rarefied gases,
especially when $\omega~{}\sim 1/\tau$, hence $\Pi$ has its evolution equation
together with a relaxation time. Its simplest version is called Meixner’s
theory [392]. In RET, Meixner’s theory is called ET6 ($6$-field) theory, in
which the Euler fluid model is extended with a relaxation equation in the
rarefied gas limit, it reads

$\displaystyle\dot{\rho}+\rho\frac{\partial v_{k}}{\partial x_{k}}$
$\displaystyle=0,$ (171) $\displaystyle\rho\dot{v}_{i}+\frac{\partial
p}{\partial x_{i}}+\frac{\partial\Pi}{\partial x_{i}}$ $\displaystyle=0,$
(172)
$\displaystyle\dot{T}+\frac{2m}{Dk_{\textrm{B}}\rho}\Big{(}p+\Pi\Big{)}\frac{\partial
v_{k}}{\partial x_{k}}$ $\displaystyle=0,$ (173)
$\displaystyle\tau_{\Pi}\dot{\Pi}+\Pi+\eta\frac{\partial v_{k}}{\partial
x_{k}}+\frac{5D-6}{3D}\Pi\tau_{\Pi}\frac{\partial v_{k}}{\partial x_{k}}$
$\displaystyle=0,$ (174)

where the relaxation time $\tau_{\Pi}$ and the bulk viscosity $\eta$ are given
as a function of degrees of freedom $D$,

$\displaystyle\tau_{\Pi}=\frac{2pT(D-3)}{3D\zeta},\quad\eta=\frac{2(D-3)}{3D}p\tau_{\Pi},\quad
p=\frac{k_{\textrm{B}}}{m}\rho T,\quad
e=\frac{D}{2}\frac{k_{\textrm{B}}}{m}T,\quad\zeta=\zeta(\rho,T)>0\ \
\textrm{for any}\ \ \\{\rho,T\\}.$ (175)

It is visible that Meixner’s theory reduces to the Euler fluid for monatomic
gases as $D=3\Rightarrow\tau_{\Pi}=0,\eta=0\Rightarrow\Pi=0$ [393, 376]. The
system can only be obtained from kinetic theory if the phase space of $f$ is
extended with a variable denoted by $I$ and that is related to (but not
identical to) the internal degrees of freedom [40], thus
$f=f(t,\mathbf{x},\mathbf{c},I)$, together with a non-negative measure
$\varphi(I)\textrm{d}I$, and $\varphi(I)=I^{n}$, $n=(D-5)/2$, called state
density. It yields

$\displaystyle 2\rho
e=\int\displaylimits_{-\infty}^{\infty}\int\displaylimits_{0}^{\infty}(mC^{2}+2I)f(t,\mathbf{x},\mathbf{c},I)\varphi(I)\textrm{d}I\textrm{d}\mathbf{c},\quad
3(p+\Pi)=\int\displaylimits_{-\infty}^{\infty}\int\displaylimits_{0}^{\infty}mC^{2}f(t,\mathbf{x},\mathbf{c},I)\varphi(I)\textrm{d}I\textrm{d}\mathbf{c}.$
(176)

Sadly, it is proved that the Eqs. (171)-(174) cannot adequately model rarefied
gases, thus further extensions are necessary. Since the $13$-moment equation
proved to have stability and convergence issues, it is reasonable to consider
the ET14 theory (14-momentum equations) as the next candidate. One has to
consider the $F$ and $G$-hierarchies plotted in Eqs. (5.4.1) and (5.4.1), and
thus the balances are

$\displaystyle\dot{\rho}+\rho\partial_{k}v_{k}$ $\displaystyle=0,$ (177)
$\displaystyle\rho\dot{v}_{i}+\partial_{j}\left[(p+\Pi)\delta_{ij}+\Pi_{\langle
ij\rangle}\right]$ $\displaystyle=0_{i},$ (178)
$\displaystyle\rho\dot{e}+\left[(p+\Pi)\delta_{ij}+\Pi_{\langle
ij\rangle}\right]\partial_{i}v_{j}+\partial_{i}q_{i}$ $\displaystyle=0,$ (179)

and their closure results

$\displaystyle\dot{\Pi}+\left(\frac{2\hat{c}_{v}-3}{3\hat{c}_{v}}p+\frac{5\hat{c}_{v}-3}{3\hat{c}_{v}}\Pi\right)\partial_{k}v_{k}+\frac{2\hat{c}_{v}-3}{3\hat{c}_{v}}\Pi_{\langle
ik\rangle}\partial_{\langle
i}v_{k\rangle}-\frac{5}{3}\frac{1}{(1+\hat{c}_{v})^{2}}\frac{\textrm{d}\hat{c}_{v}}{\textrm{d}T}q_{k}\partial_{k}T+\frac{2\hat{c}_{v}-3}{3\hat{c}_{v}(1+\hat{c}_{v})}\partial_{k}q_{k}$
$\displaystyle=-\frac{1}{\tau_{\Pi}}\Pi,$ (180)
$\displaystyle\dot{\Pi}_{\langle ij\rangle}+\Pi_{\langle
ij\rangle}\partial_{k}v_{k}+2\partial_{k}v_{\langle i}\Pi_{\langle j\rangle
k\rangle}+2(p+\Pi)\partial_{\langle
i}v_{j\rangle}-\frac{2}{(1+\hat{c}_{v})^{2}}\frac{\textrm{d}\hat{c}_{v}}{\textrm{d}T}\partial_{k}Tq_{\langle
i}\delta_{j\rangle k}+\frac{2}{1+\hat{c}_{v}}\partial_{\langle j}q_{i\rangle}$
$\displaystyle=-\frac{1}{\tau_{S}}\Pi_{\langle ij\rangle},$ (181)

for the viscous and dynamic pressure with $\hat{c}_{v}=c_{v}m/k_{\textrm{B}}$,
and for the heat flux,

$\displaystyle{\dot{q}_{i}+\frac{2+\hat{c}_{v}}{1+\hat{c}_{v}}q_{i}\partial_{k}v_{k}+\frac{1}{1+\hat{c}_{v}}q_{k}\partial_{i}v_{k}+\frac{2+\hat{c}_{v}}{1+\hat{c}_{v}}q_{k}\partial_{k}v_{i}}+\frac{k_{\mathrm{B}}}{m}\left\\{(1+\hat{c}_{v})p\delta_{ki}+(2+\hat{c}_{v})(\Pi\delta_{ki}+\Pi_{\langle
ki\rangle})\right\\}\partial_{k}T$
$\displaystyle{\quad-\frac{k_{\mathrm{B}}}{m}T\partial_{i}p+\frac{1}{\rho}\left\\{(p-\Pi)\delta_{ki}-\Pi_{\langle
ki\rangle}\right\\}\partial_{l}\left\\{(p+\Pi)\delta_{kl}+\Pi_{\langle
kl\rangle}\right\\}=-\frac{1}{\tau_{q}}q_{i}.}$ (182)

The linearized form of Eqs. (180)-(5.4.1) could have significant practical
relevance, as it turned out from the evaluations of acoustic experimental data
[394]. ET14 has its counterpart from a continuum point of view [275], however,
their compatibility is not yet fully discovered in the nonlinear regime.
Overall, the linearized-ET14 equations are

$\displaystyle\dot{\Pi}+\frac{2\hat{c}_{v}^{(0)}-3}{3\hat{c}_{v}^{(0)}}p^{(0)}\partial_{k}v_{k}+\frac{2\hat{c}_{v}^{(0)}-3}{3\hat{c}_{v}^{(0)}(1+\hat{c}_{v}^{(0)})}\partial_{k}q_{k}$
$\displaystyle=-\frac{1}{\tau_{\Pi}}\Pi,$ (183)
$\displaystyle\dot{\Pi}_{\langle ij\rangle}+2p^{(0)}\partial_{\langle
j}v_{i\rangle}+\frac{2}{1+\hat{c}_{v}^{(0)}}\partial_{\langle j}q_{i\rangle}$
$\displaystyle=-\frac{1}{\tau_{S}}\Pi_{\langle ij\rangle},$ (184)
$\displaystyle\dot{q}_{i}-\frac{k_{\mathrm{B}}}{m}T^{(0)}\partial_{i}p+\frac{k_{\mathrm{B}}}{m}(1+\hat{c}_{v}^{(0)})p^{(0)}\partial_{i}T+\frac{p^{(0)}}{\rho^{(0)}}\partial_{k}\left\\{(p+\Pi)\delta_{ik}\Pi_{\langle
ik\rangle}\right\\}$ $\displaystyle=-\frac{1}{\tau_{q}}q_{i},$ (185)

for which the linearization is performed around a set of reference values
$(\rho^{(0)},v_{i}^{(0)},T^{(0)},\Pi=0,\Pi_{\langle
ij\rangle}=0_{ij},q_{i}=0_{i})=\mathrm{const}$. Eqs. (183)-(185) can exactly
be reproduced with both EIT and NET-IV approaches [275]. The closure of the
nonlinear model (177)-(5.4.1) yields the following (approximate) entropy
density and entropy flux,

$\displaystyle s(e,\rho,q_{i},\Pi,\Pi_{\langle ij\rangle})$
$\displaystyle=s_{\textrm{eq}}(e,\rho)-\frac{3\hat{c}_{v}}{2(2\hat{c}_{v}-3)pT}\Pi^{2}-\frac{1}{4\rho
T}\Pi_{\langle ij\rangle}\Pi_{\langle
ij\rangle}-\frac{\rho}{2p^{2}T(1+\hat{c}_{v})}q_{i}q_{i}+\mathcal{O}(3),$
(186) $\displaystyle J_{k}$
$\displaystyle=\frac{1}{T}q_{k}-\frac{1}{pT(1+\hat{c}_{v})}\Pi
q_{k}-\frac{1}{pT(1+\hat{c}_{v})}q_{i}\Pi_{\langle
ik\rangle}+\mathcal{O}(3),\quad 2\hat{c}_{v}-3>0,$ (187)

in which the third-order non-equilibrium terms are neglected. Interestingly,
while these entropy relations are found in the closure procedure as an
outcome, EIT and NET-IV thermodynamic approaches use them as a starting point,
constraining the evolution equations from the beginning. Hence Eqs.
(186)-(187) serve as a strong connection between various thermodynamic
approaches and provide insight into the compatibility conditions.
Additionally, from Eq. (187), the appearance of a current multiplier, the
viscous pressure is clear, however, it must be separated into spherical and
deviatoric parts in order to implement the proper coupling to the heat flux.

Further refining on the ET14 model is possible by splitting the variable $I$
into rotational $I^{\textrm{r}}$ and vibrational modes $I^{\textrm{v}}$, and
that leads to a triple ($F$, $G$ and $H$) hierarchies, using
$f=f(t,\mathbf{x},\mathbf{c},I^{\textrm{r}},I^{\textrm{v}})$ with the
corresponding state densities. It does not influence the derivation procedure,
and the resulting equations again fulfill the second law of thermodynamics and
Galilean covariance. For further details, we refer to the works of Ruggeri et
al. [375, 70].

####  5.4.2. Continuum models.

Although, in the light of knowing the relations (186)-(187), it would be
straightforward to achieve compatibility, at least in the linear regime, this
is not the primary aim for a continuum approach. Kinetic theory and RET models
are specific, valid only for ideal gases, and implement various molecule
properties. On the contrary, the continuum approach can be more adaptable,
facilitate the adaptation of measured state-dependent transport coefficients,
and extend the region of validity towards the dense states. Moreover, in the
analogy of the effective modeling we have seen concerning heterogeneous
materials, the generalized NSF equations could be helpful in modeling two-
phase flows for which the droplets in the gas flow behave like molecules in
the rarefied flow, but these analogies are not yet experimentally tested.
Thus, compatibility is a secondary but still important attribute in comparing
and benchmarking different approaches, as it is a reasonable expectation that
all (acceptable) theories should reproduce the modeling of certain phenomena.
In that case, the compatibility is tested on the rarefied gas experiments
mentioned earlier [374, 275, 395].

The EIT approach for rarefied gas modeling results in a model more similar to
the R13 equations than the system of ET14 since the dynamic pressure $\Pi$ is
not among the state variables, but the current multipliers in the entropy flux
are different [396, 397, 398],

$\displaystyle\mathbf{J}_{s}=\frac{1}{T}\mathbf{q}+\alpha_{1}\boldsymbol{\Pi}^{\textrm{dev}}\cdot\mathbf{q}+\alpha_{2}\nabla\mathbf{q}\cdot\mathbf{q}+\alpha_{3}\nabla\boldsymbol{\Pi}^{\textrm{dev}}:\boldsymbol{\Pi}^{\textrm{dev}},$
(188)

and the non-local $\nabla\mathbf{q}$ and
$\nabla\boldsymbol{\Pi}^{\textrm{dev}}$ terms make the resulting evolution
equations to be parabolic. Together with the $s=s(e,\rho,q_{i},\Pi_{\langle
ij\rangle})$ entropy density, the evolution equations are

$\displaystyle\tau_{q}\dot{\mathbf{q}}+\mathbf{q}$
$\displaystyle=-\lambda\nabla
T-\lambda_{2}\boldsymbol{\Pi}^{\textrm{dev}}\cdot\nabla-\lambda_{3}\Delta\mathbf{q},$
(189)
$\displaystyle\tau_{S}\dot{\boldsymbol{\Pi}}^{\textrm{dev}}+\boldsymbol{\Pi}^{\textrm{dev}}$
$\displaystyle=-\mu(\nabla\mathbf{v})^{\textrm{dev}}-\mu_{2}(\nabla\mathbf{q})^{\textrm{dev}}-\mu_{3}\Delta\boldsymbol{\Pi}^{\textrm{dev}},$
(190)

with the coefficients

$\displaystyle\tau_{q}=\frac{3}{2}\tau_{S},\quad\lambda=\frac{15}{4}p\tau_{S},\quad\lambda_{2}=\frac{3}{2}RT\tau_{S},\quad\mu=2p\tau_{S},\quad\mu_{2}=\frac{4}{5}\tau_{S},\quad\mu_{3}=-RT\tau_{S}^{2}$
(191)

inheriting the coefficients from the R13-moment equations (166)-(167), except
$\lambda_{3}$, which could also be found from the Eq. (168), and the nonlinear
terms are not present here. Overall, both $\lambda_{3}$ and $\tau_{S}$ are to
be fitted to experimental data. We note that such correspondence between the
coefficients remains arbitrary, but for a rarefied gas model, it seems
reasonable and advantageous to adopt the kinetic relations to reduce the
number of free parameters from $6$ to $2$ (or $1$, if suitable).

Figure 14. Model benchmarking for two different problems. First row: shock
structure studies, in which the NSF (A), Grad13 (B), and Burnett (C) equations
are compared to direct-simulation Monte Carlo data [399], $\lambda_{0}$
denotes the mean free path. Second row: velocity measurements compared to NSF,
ET14 (D) [394], and NET-IV (E) theories [374]. Additionally, in part E), the
original $\omega/p$ scaling is exchanged for mass density as the frequency is
constant throughout the measurement.

The situation becomes slightly more complicated for the NET-IV since the
corresponding continuum model aims to achieve the ET14 level [275].
Consequently, the dynamic pressure is also present in the state space,
$s=s(e,\rho,q_{i},\Pi,\Pi_{\langle ij\rangle})$, and
$s=s_{\textrm{eq}}(e,\rho)-\frac{m_{1}}{2}q_{i}q_{i}-\frac{m_{2}}{2}\Pi_{\langle
ij\rangle}\Pi_{\langle ij\rangle}-\frac{m_{3}}{6}\Pi_{ii}\Pi_{jj}$, where
$m_{1},m_{2},m_{3}>0$ are positive coefficients and, similarly to the previous
cases, they could also be a positive function of the state variables.
Regarding the entropy flux, adopting (187) is also possible, leading to
further nonlinearities. The potential of these latter (nonlinear) cases are
not yet discovered; hence we apply only constant $m_{i}$ ($i=1,2,3$)
coefficients with $J_{i}=(b_{\langle ij\rangle}+b_{kk}\delta_{ij}/3)q_{j}$,
for which the second law provides the necessary closure to find $b_{\langle
ij\rangle}$ and $b_{kk}$, still keeping the possibility to be compatible with
(187). Hence the Onsagerian relations are

$\displaystyle-\rho
m_{1}\dot{q}_{i}+\frac{1}{3}\partial_{i}b_{kk}+\partial_{j}b_{\langle
ji\rangle}$ $\displaystyle=$ $\displaystyle nq_{i},$ (192)
$\displaystyle-\frac{1}{T}\partial_{\langle i}v_{j\rangle}-\rho
m_{2}\dot{\Pi}_{\langle ij\rangle}$ $\displaystyle=$ $\displaystyle
l_{11}\Pi_{\langle ij\rangle}+l_{12}\partial_{\langle i}q_{j\rangle},$ (193)
$\displaystyle b_{\langle ij\rangle}$ $\displaystyle=$ $\displaystyle
l_{21}\Pi_{\langle ij\rangle}+l_{22}\partial_{\langle i}q_{j\rangle},$ (194)
$\displaystyle-\frac{1}{T}\partial_{j}v_{j}-\rho m_{3}\dot{\Pi}_{ii}$
$\displaystyle=$ $\displaystyle
k_{11}\frac{\Pi_{ii}}{3}+k_{12}\partial_{i}q_{i},$ (195) $\displaystyle
b_{kk}-\frac{1}{T}$ $\displaystyle=$ $\displaystyle
k_{21}\frac{\Pi_{ii}}{3}+k_{22}\partial_{i}q_{i}.$ (196)

where the coefficients must fulfill

$\displaystyle n\geq 0,\quad l_{11}l_{22}-(l_{12}+l_{21})^{2}/4\geq 0,\quad
k_{11}k_{22}-(k_{12}+k_{21})^{2}/4\geq 0.$ (197)

to preserve the positive definiteness of the entropy production. Compatibility
can be achieved by making the system (192)-(196) to be hyperbolic with
implying $l_{22}=k_{22}=0$ to eliminate the nonlocal terms. Therefore, to
satisfy (197), $l_{12}=-l_{21}$ and $k_{12}=-k_{21}$ holds, and yields

$\displaystyle\tau_{q}\dot{q}_{i}+q_{i}+\lambda\partial_{i}T-\alpha_{21}\partial_{i}\Pi_{kk}-\beta_{21}\partial_{i}\Pi_{\langle
ij\rangle}$ $\displaystyle=0,$ (198) $\displaystyle\tau_{S}\dot{\Pi}_{\langle
ij\rangle}+\Pi_{\langle ij\rangle}+\mu\partial_{\langle
i}v_{j\rangle}+\beta_{12}\partial_{\langle i}q_{j\rangle}$ $\displaystyle=0,$
(199)
$\displaystyle\tau_{\Pi}\dot{\Pi}_{ii}+\Pi_{ii}+\eta\partial_{i}v_{i}+\alpha_{12}\partial_{i}q_{i}$
$\displaystyle=0,$ (200)

for which the coefficients can either be adopted from the ET14 equations or to
be fitted, and they constitute the following Onsagerian coefficients:

$\displaystyle\tau_{q}=\frac{\rho m_{1}}{n},\quad\tau_{S}=\frac{\rho
m_{2}}{l_{11}},\quad\tau_{\Pi}=\frac{3\rho
m_{3}}{k_{11}},\quad\lambda=\frac{1}{3nT^{2}},\quad\mu=\frac{1}{Tl_{11}},\quad\eta=\frac{3}{Tk_{11}},$
(201)
$\displaystyle\alpha_{12}=\frac{3k_{12}}{k_{11}}\quad\alpha_{21}=\frac{k_{21}}{9n},\quad\beta_{12}=\frac{l_{12}}{l_{11}},\quad\beta_{21}=\frac{l_{21}}{n}.$
(202)

We want to emphasize again that the NET-IV model (198)-(200) are valid only
for constant coefficients, and the nonlinear options are not yet discovered.
It is straightforward to accept that the NET-IV model offers much freedom,
which freedom is not necessary to exploit completely. It remains a matter of
choice how the particular modeling task requires, primarily when a measured
state dependence of the transport coefficients must be implemented. The
compatibility between these approaches is proved and extensively discussed in
[275]. Figure 14 presents two benchmarking situations concerning shock
structures [399, 400, 401] and speed of sound measurements [374, 368],
comparing both the kinetic and continuum models to experimental data and
detailed Monte Carlo simulations.

##  6\. Dual-phase-lag concept

Tzou’s original idea was to unify the various heat equations and find their
common root on a constitutive level [43]. Tzou observed the resemblance
between the $T$-representation of various heat equations, including two-
temperature models, and proposed the DPL concept as

$\displaystyle\mathbf{q}(\mathbf{x},t+\tau_{q})=-\lambda\nabla
T(\mathbf{x},t+\tau_{T}),,$ (203)

assuming a delay between the heat flux and temperature gradient, contradicting
the time homogeneity principle, and that Eq. (203) can approximate numerous
heat equations, but not in the same way. While the first-order Taylor series
expansion of the left hand side yields the MCV constitutive equation, the
structure of the GK equation is recovered solely in its $T$-representation.
Eq. (203) reduces to Fourier’s law for $\tau_{q}=\tau_{T}=0$ but also
reproduces Fourier’s solutions for any $\tau_{q}=\tau_{T}$, likewise to the
Fourier resonance condition [104]. Although Tzou mentions the possible
microstructural effects behind non-Fourier heat conduction, it is also claimed
that the time lags $\tau_{q}$ and $\tau_{T}$ are effective parameters,
collectively modeling the corresponding microstructural phenomena on a
macroscopic level [43]. Therefore, in the absence of any (thermodynamic or
kinetic) theory, it is unattainable to attach direct interpretation to these
coefficients [16].

###  6.1. Jeffreys heat equation.

There are two, partially valid approximations of Eq. (203), the MCV equation
($\\{1,0\\}$-type DPL, i.e., first order in $\mathbf{q}$, and $\tau_{T}=0$);
the second one is the Jeffreys equation ($\\{1,1\\}$-type) [402, 2]. The
expression ‘partially’ is reasonable since Eq. (203) has no such thermodynamic
background, which has been developed for these models, and thus there are
consequences, but first, let us present the Jeffreys equation.

Its derivation is more common in NET-IV [40], compatible with GENERIC [168],
and cannot be derived within the framework of RET. The corresponding state
space is $s=s(e,\mathbf{y})$, where $\mathbf{y}$ ($\neq\mathbf{q}$) is an
extensive vectorial state variable, and thus
$s(e,\mathbf{y})=s_{\textrm{eq}}(e)-\frac{m}{2}\mathbf{y}\cdot\mathbf{y}$
($m>0$), and the classical entropy flux is applied
($\mathbf{J}_{s}=\mathbf{q}/T$). The resulting Onsagerian relations are

$\displaystyle\mathbf{q}$
$\displaystyle=l_{11}\nabla\frac{1}{T}-l_{12}m\mathbf{y},$ (204)
$\displaystyle\rho\partial_{t}\mathbf{y}$
$\displaystyle=l_{21}\nabla\frac{1}{T}-l_{22}m\mathbf{y},$ (205)

with the requirements on $l_{ij}$,

$\displaystyle l_{11}\geq 0,\quad l_{22}\geq 0,\quad
l_{11}l_{22}-l_{12}l_{21}\geq 0.$ (206)

After eliminating $\mathbf{y}$, the Jeffreys-type equation is formed,

$\displaystyle\tau\frac{\partial\mathbf{q}}{\partial
t}+\mathbf{q}=\lambda_{1}\nabla\frac{1}{T}+\lambda_{2}\frac{\partial}{\partial
t}\nabla\frac{1}{T},$ (207)

with the coefficients

$\displaystyle\tau=\frac{\rho}{l_{22}m},\quad\lambda_{1}=\frac{\textrm{det}l_{ij}}{l_{22}},\quad\lambda_{2}=\frac{\rho
l_{11}}{l_{22}m}=l_{11}\tau,$ (208)

and remarkably, the $\\{0,1\\}$-type model is excluded in virtue of (208), and
the MCV equation ($\\{1,0\\}$-type) is recovered when $l_{11}=0$. Furthermore,
the Fourier resonance condition requires that

$\displaystyle\lambda_{1}=\frac{\lambda_{2}}{\tau}\Rightarrow\frac{\textrm{det}l_{ij}}{l_{22}}=l_{11}\Rightarrow\frac{l_{12}l_{21}}{l_{22}}=0,$
(209)

that is, the coupling between $\mathbf{y}$ and $\mathbf{q}$ vanishes, and
$\mathbf{y}$ relaxes towards an asymptotically stable equilibrium point.
Finally, we note that Rukolaine found nonphysical solutions for the Eq. (207)
[403, 404]. However, in that test problem, a particular source term includes a
Heaviside step function in time, which could fall from the physically
admissible function space for $T$ and $\mathbf{q}$. Moreover, apparently, the
initial conditions ($T=$const., and $\partial_{t}T=0$) are not
thermodynamically compatible for a space-dependent source term, and thus
further investigations and experimental testing might be necessary before
excluding the Jeffreys heat equation from the physically admissible models.
Nevertheless, the DPL concept introduces further disadvantages for higher-
order expansions of (203).

The $T$-representation of Eq. (207) formally resembles to the GK equation in
the linear regime, i.e., when $2\lambda_{2}T^{-3}\partial_{t}T\nabla T$ is
omitted from the last term, and all coefficients are constants, thus

$\displaystyle\tau\partial_{tt}T+\partial_{t}T=\alpha_{1}\Delta
T+\alpha_{2}\partial_{t}\Delta T+\frac{Q_{v}}{\rho c_{v}}+\frac{\tau}{\rho
c_{v}}\partial_{t}Q_{v},\quad\alpha_{1,2}=\frac{\hat{\lambda}_{1,2}}{\rho
c_{v}},\quad\hat{\lambda}_{1,2}=\frac{\lambda_{1,2}}{T^{2}}.$ (210)

The $q$-representation also bears a similar structure,

$\displaystyle\tau\partial_{tt}\mathbf{q}+\partial_{t}\mathbf{q}=\alpha_{1}\nabla\nabla\cdot\mathbf{q}+\alpha_{2}\nabla\nabla\cdot\partial_{t}\mathbf{q}-\alpha_{1}\nabla
Q_{v}-\alpha_{2}\nabla\partial_{t}Q_{v}.$ (211)

These forms are valid only for linear models with constant coefficients,
otherwise it is not suggested to eliminate any field variables. For higher-
order approximations of Eq. (203), the model structure is preserved, only
higher-order time derivatives enter either the $T$ or the $q$-representations.

###  6.2. Stability conditions.

Although finding a constitutive equation capable of connecting or unifying
many other heat equations is interesting, it also has limitations and
shortcomings. First, the Taylor series expansion on both sides remains
arbitrary, and its convergence is not proven. Even the meaning of the
relaxation times can change with the order of expansion, however, every
expansion approximates the same model (203). Due to the lack of thermodynamic
background, its compatibility with the second law remains questionable and
could change with the expansion order. Further aspects emerge as follows.

  * •

Nonlinearities: according to the second law of thermodynamics, the
coefficients are not independent. Here, no functional connections are
established, therefore, it is not recommended to apply the model for, e.g.,
$T$-dependent coefficients.

  * •

Anisotropy: the tensorial properties of each physical quantity become crucial
for anisotropic materials, e.g., the scalar relaxation times become a second-
order tensor, together with the thermal conductivity. It could allow further
couplings and result in a complex time evolution equation for non-Fourier heat
equations. Starting with Eq. (203) restricts the modeling capabilities on the
isotropic behavior.

  * •

Coupled problems: thermodynamics provides a consistent approach in deriving
coupled equations such as thermo-diffusion [29] and thermo-mechanics [285].
Due to the missing background of Eq. (203), the DPL approach could only be
helpful for pure heat conduction problems.

  * •

Time shift paradox: sadly, Eq. (203) directly contradicts the homogeneity of
time; only the difference between the two relaxation times should be essential
[405]. This is in agreement with Rukolaine [403], Fabrizio and Franchi [406],
and it turned out that $\tau_{T}<\tau_{q}$ must be satisfied to have a well-
posed model. Otherwise, stability and ill-posedness issues arise.

Because of these shortcomings, the region of validity of the DPL equations is
firmly limited. For instance, many papers present experimental data
evaluations using the DPL model without checking the stability conditions. The
positivity of relaxation times is not sufficient. The $\\{2,1\\}$-type DPL,
i.e., it is second-order for $\mathbf{q}$, and first-order for $\nabla T$,

$\displaystyle\frac{\tau_{q}^{2}}{2}\frac{\partial^{2}\mathbf{q}}{\partial
t^{2}}+\tau_{q}\frac{\partial\mathbf{q}}{\partial t}+\mathbf{q}=-\lambda\nabla
T-\lambda\tau_{T}\frac{\partial}{\partial t}\nabla T,$ (212)

has stable solution if [407],

$\displaystyle\frac{\tau_{T}}{\tau_{q}}\geq\frac{1}{2}.$ (213)

For a $\\{2,2\\}$-type DPL, the stability condition (213) modifies to [408,
409],

$\displaystyle 2+\sqrt{3}>\frac{\tau_{T}}{\tau_{q}}>2-\sqrt{3},$ (214)

in other words, the stability conditions are different for any approximation
of Eq. (203), and one must be cautious when choosing a seemingly suitable
approximation. The over-diffusive region ($\tau_{T}>\tau_{q}$) is preferred by
(213) but begins to be strongly restrained. For further mathematical analysis,
we refer to the study of Shen and Zhang [410].

Now let us turn our attention to the asymptotic behavior of the DPL model with
a series expansion up to arbitrary orders of $n$ and $m$,

$\displaystyle\mathbf{q}+\frac{\tau_{q}}{1!}\partial_{t}\mathbf{q}+\frac{\tau_{q}^{2}}{2!}\partial_{t}^{2}\mathbf{q}+\dots+\frac{\tau_{q}^{n}}{n!}\partial_{t}^{n}\mathbf{q}=-\lambda\left(\nabla
T+\frac{\tau_{T}}{1!}\partial_{t}\nabla
T+\frac{\tau_{T}^{2}}{2!}\partial_{t}^{2}\nabla
T+\dots+\frac{\tau_{T}^{m}}{m!}\partial_{t}^{m}\nabla T\right).$ (215)

It turned out by investigating the characteristic equations of the resulting
operators that if $n$ (or $m$) is $\geq 5$, then there will be at least one
pair of complex conjugate roots with positive real parts, i.e., instability
persists without any conditions [411]. In other words, the Taylor series
expansion of such constitutive equations is not convergent.

Figure 15. Roots of the characteristic equations. Left: all roots have
negative real parts. Right: after $n\geq 5$, a root with a positive real part
appears and remains (red squares highlight them), always resulting in
instability [411].

The further merit of [411] emerges by investigating all situations up to the
order of {4,4}-type DPL from a thermodynamic compatibility point of view,
exploiting the fading memory concept of Gurtin and Pipkin [173]. The $n,m\geq
5$ condition does not mean that all combinations of $n,m<5$ are
thermodynamically compatible. This is studied in view of the Clausius-Duhem
inequality formulated for cycles with period $p=2\pi/\omega$,

$\displaystyle\int\displaylimits_{0}^{p}\mathbf{q}(t)\cdot\nabla
T(t)\textrm{d}t\leq 0,$ (216)

with $\nabla
T(t-s)=\mathbf{f}\cos\big{(}\omega(t-s)\big{)}+\mathbf{q}\sin\big{(}\omega(t-s)\big{)}$,
and $\omega>0$, $\mathbf{f}^{2}+\mathbf{g}^{2}>0$.

We have already shown that the {1,0} and {1,1}-types are thermodynamically
compatible. Furthermore, the {0,1}-type, which can be rewritten as

$\displaystyle\nabla
T(\mathbf{x},t)=-\frac{1}{\tau_{T}}\int\displaylimits_{0}^{\infty}e^{-s/\tau_{T}}\lambda\mathbf{q}(\mathbf{x},t-s)\textrm{d}s,$
(217)

is still thermodynamically compatible for all $\tau_{T}>0$ in the light of
(216). However, the {2,0}-type model cannot satisfy (216),

$\displaystyle\int\displaylimits_{0}^{p}\mathbf{q}(t)\cdot\nabla
T(t)\textrm{d}t=-\frac{2\pi\lambda(\mathbf{f}^{2}+\mathbf{g}^{2})}{\omega(\tau_{q}^{4}\omega^{4}+4)}(2-\tau_{q}^{2}\omega^{2})$
(218)

cannot preserve the sign for all $\omega>0$, hence thermodynamically
incompatible, similarly to the {0,2}-type. Previous stability conditions can
be achieved similarly [411]. For higher-order models, the stability conditions
become more and more complicated. For demonstration, we show only the
condition for the {3,3}-type, and for the others, we refer to [411] for
convenience. The stability condition of the {3,3}-type model is

$\displaystyle\sqrt{\left(b_{1}^{2}+3a_{1}c_{1}\right)^{3}}\leq\frac{1}{2}\left(27a_{1}^{2}d_{1}-2b_{1}^{3}-9a_{1}b_{1}c_{1}\right),$
(219)

with the coefficients as a function of $\tau_{T}/\tau_{q}$,

$\displaystyle a_{1}=0.677637\frac{\tau_{T}^{3}}{\tau_{q}^{3}};\quad
b_{1}=\frac{\tau_{T}}{\tau_{q}}\left(4.06574\frac{\tau_{T}^{2}}{\tau_{q}^{2}}-6.09877\frac{\tau_{T}}{\tau_{q}}+4.06582\right);$
$\displaystyle
c_{1}=12.1972\left(\frac{\tau_{T}}{\tau_{q}}-1\right)^{2}+0.0003;\quad
d_{1}=24.3944.$ (220)

So the conditions summed up in Table 2 are required to achieve thermodynamic
compatibility. For well-posedness, proper initial and boundary conditions are
also needed, and in that respect, we want to recall the difficulties mentioned
concerning the MCV equation in Sec. 3.8. For a more thorough review of the DPL
concept and its applications, we refer to [222, 225, 118, 412, 413, 414], and
in regard to the well-posedness and stability properties, we refer to the
papers of Quintanilla [415, 416, 417].

Table 2. Summary of the possible DPL models, based on [411], and ✗ means that the resulting model is unconditionally incompatible with thermodynamics. Model type | Compatibility | Condition | Model type | Compatibility
---|---|---|---|---
{0,0} |  | - | {2,0} | ✗
{1,0} |  | $\tau_{q}>0$ | {0,2} | ✗
{0,1} |  | $\tau_{T}>0$ | {3,0} | ✗
{1,1} |  |  $\tau_{q}>0$, $\tau_{T}>0$ | {0,3} | ✗
{2,1} |  | $\tau_{T}\geq\tau_{q}/2$ | {4,0} | ✗
{1,2} |  | $\tau_{q}\geq\tau_{T}/2$ | {0,4} | ✗
{2,2} |  | $2+\sqrt{3}>\frac{\tau_{T}}{\tau_{q}}>2-\sqrt{3}$ | {3,1} | ✗
{3,2} |  | $1.4902\tau_{q}\geq\tau_{T}>0.28441\tau_{q}$ | {1,3} | ✗
{2,3} |  | $1.4902\tau_{T}\geq\tau_{q}>0.28441\tau_{T}$ | {0,0} | ✗
{3,3} |  | See Eqs. (219)-(220). | {4,1} | ✗
{4,3} |  | See [411]. | {1,4} | ✗
{3,4} |  | See [411]. | {4,2} | ✗
{4,4} |  | See [411]. | {2,4} | ✗

###  6.3. Generalized DPL model.

Practically, the DPL concept is the only non-Fourier model applied in
biological heat conduction problems, therefore, it is inevitable to summarize
the essential aspects here briefly. It is also crucial to emphasize that for a
macroscopic heterogeneous system in a room temperature environment, such a
model should only be interpreted as an effective model, the time lags are due
to the interaction of multiple coupled phenomena. Therefore only over-
diffusive solutions are expected, thermal waves (under-diffusive) do not exist
under such conditions, and that is often misinterpreted in the literature so
that $\tau_{T}>\tau_{q}$ should be satisfied for physically sound solutions,
e.g., for copper, $\tau_{T}/\tau_{q}=163$ [284].

Let us recall the two-temperature model of Roetzel and Xuan [114],

$\displaystyle\varepsilon\rho_{b}c_{b}\left(\partial_{t}T_{b}+\mathbf{v}_{b}\cdot\nabla
T_{b}\right)$ $\displaystyle=\nabla\cdot(\lambda_{b}\nabla
T_{b})+h(T_{t}-T_{b}),$ (221)
$\displaystyle(1-\varepsilon)\rho_{t}c_{t}\partial_{t}T_{t}$
$\displaystyle=\nabla\cdot(\lambda_{t}\nabla
T_{t})-h(T_{t}-T_{b})+(1-\varepsilon)q_{\textrm{met}},$ (222)

as it serves as the basis for the analysis of Zhang [222]. After eliminating
the blood temperature from Eqs. (221) and (222), with exploiting that there is
an additional relaxation process between the tissue and blood following [418],

$\displaystyle\frac{\varepsilon\rho_{b}c_{b}}{h}\partial_{t}T_{b}+T_{b}=T_{t},$
(223)

they obtain an evolution equation for the tissue,

$\displaystyle\tau_{q}\partial_{tt}T_{t}+\partial_{t}T_{t}+\frac{\varepsilon\rho_{b}c_{b}}{(\rho
c)_{\textrm{eff}}}\mathbf{v}_{b}\cdot\nabla
T_{t}=\alpha_{\textrm{eff}}\Big{(}\Delta T_{t}+\tau_{T}\partial_{t}\Delta
T_{t}\Big{)}+\left(\frac{1-\varepsilon}{(\rho
c)_{\textrm{eff}}}+\frac{\varepsilon\rho_{b}c_{b}}{h(\rho
c)_{\textrm{eff}}}\partial_{t}\right)q_{\textrm{met}}$ (224)

in which metabolic and convective heat transfer sources are included, and the
acceleration of the blood flow is omitted. This is called the generalized DPL
model. The effective quantities are

$\displaystyle(\rho
c)_{\textrm{eff}}=\varepsilon\rho_{b}c_{b}+(1-\varepsilon)\rho_{t}c_{t},\quad\lambda_{\textrm{eff}}=\varepsilon\lambda_{b}+(1-\varepsilon)\lambda_{t},\quad\alpha_{\textrm{eff}}=\frac{\lambda_{\textrm{eff}}}{(\rho
c)_{\textrm{eff}}},$ (225)

for which let us recall that $\varepsilon$ is the known porosity of a given
tissue, and the indices refer to the blood and tissue properties. Here, the
notable result is that in such a way, $\tau_{q}$ and $\tau_{T}$ are found
(analogously to the Eq. (84)) as

$\displaystyle\tau_{q}=\frac{\varepsilon(1-\varepsilon)\rho_{b}c_{b}\rho_{t}c_{t}}{h(\rho
c)_{\textrm{eff}}},\quad\textrm{and}\quad\tau_{T}=\frac{\varepsilon(1-\varepsilon)\rho_{b}c_{b}\lambda_{t}}{h\lambda_{\textrm{eff}}},$
(226)

which also strengthens the observation about functionally connected
coefficients, but – for instance, in [419] – it is not treated carefully, and
therefore the physical consequences are omitted, [420] is a positive
counterexample. According to [48],

$\displaystyle-\varepsilon\rho_{b}c_{b}\mathbf{v}_{b}\cdot\nabla T_{t}\approx
h(T_{b}-T_{t}),$ (227)

therefore Eq. (224) can be further simplified, and one does not need to know
the blood velocity field. Zhang tested [222] the relaxation times of Eq. (226)
for various sets of realistic parameters and found that $\tau_{T}>\tau_{q}$ is
satisfied for all cases, also confirmed by [421]. That is a promising result
and proves that the DPL concept can be helpful when the equations are handled
properly. The model, sadly, is not directly solved and tested on experimental
data.

However, $\tau_{T}>\tau_{q}$ is not always the case. Hooshmand et al. [223]
implemented Eq. (224) (without Eq. (226)) to predict the transient temperature
history under laser irradiation (the heating is included in a source term) and
compared to experimental data. Significant deviation occurred between the
observations and predictions. Besides, $\tau_{q}=16$ s $>\tau_{T}=0.05\dots
16$ s varied, but $\tau_{T}$ always kept smaller than $\tau_{q}$, similarly as
[224]. On the contrary, [225] found much better agreement with measurement
data with parameters $\tau_{T}>\tau_{q}$. If one would continue the Taylor
series expansion of that {1,1}-type model to {2,1} or {2,2}-types, then either
the problem becomes ill-posed, or the solution becomes unstable. It is of
crucial importance to keep clear the physically sound domains of coefficients.

##  7\. Green-Naghdi models

The unique concept of Green and Naghdi [422] is based on three pillars. First,
they introduce a so-called thermal displacement field $\alpha$ as

$\displaystyle\alpha=\int\displaylimits_{t_{0}}^{t}T(\mathbf{r},s)\textrm{d}s+\alpha_{0},\quad\alpha_{0}=\alpha(t_{0}),$
(228)

resembling the fading memory concept of Gurtin and Pipkin, and $\alpha$ must
be at least two times continuously differentiable both in space and time, and
$\dot{\alpha}=T$ holds. Second, they modify the heat-entropy flux relation,
i.e.,

$\displaystyle\mathbf{J}_{s}=\frac{1}{\theta}\mathbf{q},\quad\theta=\theta(T,\alpha;\Omega),\quad\Omega=\\{a,b,d_{1}\\}\in\mathbb{R}^{+},$
(229)

in which $\theta$ is a positive function and monotonous in $T$, and might
include some positive constants as well [202] ($c$ is preserved to denote the
specific heat). Third, the state space can be constructed in various ways,
three of them is included in the original study of Green and Naghdi,

$\displaystyle\mathbb{S}_{\textrm{I}}=(T,\nabla
T),\quad\mathbb{S}_{\textrm{II}}=(\alpha,T,\nabla\alpha),\quad\mathbb{S}_{\textrm{III}}=(\alpha,T,\nabla\alpha,\nabla
T),$ (230)

for which the subscript distinguishes the model type. The type-I includes
$\theta=T$, and thus leads to the classical Fourier equation. The type-II is
more interesting, as it extends Fourier’s law in a particular way,

$\displaystyle\theta=a+b\cdot
T,\quad\psi=c(\theta-\theta\ln\theta)+\frac{k}{2}\nabla\alpha\cdot\nabla\alpha,\quad\mathbf{J}_{s}=-\frac{\rho}{\frac{\partial\theta}{\partial
T}}\frac{\partial\psi}{\partial\nabla\alpha}=-\frac{\rho k}{b}\nabla\alpha,$
(231)

where $k>0$ is also a thermal conductivity, and thus the heat flux is obtained
through $\mathbf{q}=\theta\mathbf{J}_{s}=-(a+b\cdot T)\frac{\rho
k}{b}\nabla\alpha$ utilizing the particular form of Helmholtz free energy
$\psi$. When the state space is extended with mechanical variables such as the
right Cauchy-Green tensor [423, 424], then that approach is viable to derive a
thermo-mechanical model in which the coupling is realized through $\alpha$.
The nonlinear term $T\nabla\alpha$ is then omitted, and consequently, the
type-II heat equation reads

$\displaystyle\mathbf{q}=-\frac{a\rho
k}{b}\nabla\alpha,\quad\ddot{T}=\frac{ak}{cb^{2}}\Delta T,$ (232)

appearing as a wave equation for temperature, propagating without dissipation
with the characteristic speed of $\sqrt{ak/(cb^{2})}$. We emphasize that this
occurs as a special case but is not forbidden in this framework and,
therefore, cannot be compatible with any previous continuum or kinetic
approaches. Accordingly, the internal energy density is also modified,

$\displaystyle e=c\theta+\frac{k}{2}\nabla\alpha\cdot\nabla\alpha.$ (233)

The setting of $\mathbb{S}_{\textrm{III}}$ includes

$\displaystyle\theta=a+b\cdot
T+d_{1}\alpha,\quad\psi=\frac{k}{2}\nabla\alpha\cdot\nabla\alpha-\frac{d_{2}}{2}\alpha^{2}-b_{2}\alpha
T-\frac{b_{3}}{2}T^{2},\quad\mathbf{J}_{s}=-\frac{\rho}{\frac{\partial\theta}{\partial
T}}\frac{\partial\psi}{\partial\nabla\alpha}=-\frac{\rho k}{b}\nabla\alpha$
(234)

with the restrictions

$\displaystyle kd_{1}\geq 0,\quad b_{2}b-b_{3}d_{1}\geq 0,$ (235)

resulting in a heat equation

$\displaystyle\mathbf{q}=-\kappa_{1}\nabla\alpha-\kappa_{2}\nabla
T,\quad\frac{\rho
a}{b}(b_{2}\dot{\alpha}+b_{3}\ddot{\alpha})=\kappa_{1}\Delta\alpha+\kappa_{2}\Delta\dot{\alpha},\quad\textrm{or}\quad
b_{2}\dot{T}+b_{3}\ddot{T}=\frac{b}{\rho a}(\kappa_{1}\Delta
T+\kappa_{2}\Delta\dot{T}),$ (236)

with recalling $\dot{\alpha}=T$, and the transport coefficients, $\kappa_{1}$
and $\kappa_{2}$, are positive. Remarkably, Eq. (236) resembles the GK and
Jeffreys equations in its $T$-representation, however, with a completely
different constitutive background. These models are found to be well-posed
[425]. Regarding its finite element solution methods and further discussion,
we refer to [203, 426, 427]. Theoretically, the type-III model is applicable
for second sound modeling and, with the mechanical coupling, could be valid in
simulating ballistic heat conduction. In [203], a brief comparison with second
sound experiments is made, but a thorough investigation of ballistic
propagation has not been performed.

####  7.0.1. Triple-phase-lag concept

The idea comes from Choudhuri [428], combining the Green-Naghdi and DPL
models, particularly Eqs. (203) and (236), assuming that three time lags are
present,

$\displaystyle\mathbf{q}(\mathbf{x},t+\tau_{q})=-\kappa_{1}\nabla\alpha(\mathbf{x},t+\tau_{\alpha})-\kappa_{2}\nabla
T(\mathbf{x},t+\tau_{T}),$ (237)

with $0\leq\tau_{\alpha}<\tau_{T}<\tau_{q}$, i.e., immediately restricting the
model validity for wave-like propagation. Again, analogously to the DPL model,
it also becomes arbitrary at which order we choose to approximate Eq. (237) by
a Taylor series expansion, therefore, the same stability issues might occur
[429]. In more detail, following [429], let us assume a first-order expansion
in all three time lags, hence the $T$-representation of Eq. (237) is

$\displaystyle\tau_{q}\rho c_{v}\partial_{t}^{3}T+\rho
c_{v}\partial_{t}^{2}T=\kappa_{1}\Delta
T+\hat{\tau}_{\alpha}\partial_{t}\Delta
T+\kappa_{2}\tau_{T}\partial_{t}^{2}\Delta T$ (238)

with $\hat{\tau}_{\alpha}=\kappa_{1}\tau_{\alpha}+\kappa_{2}$. If
$\hat{\tau}_{\alpha}-\kappa_{1}\tau_{q}<0$, unstable solutions emerge for a
homogeneous $T=0$ boundary condition. In fact, the sufficient requirement is

$\displaystyle\phi_{1}>\frac{\rho
c_{v}(\tau_{q}\kappa_{1}-\hat{\tau}_{\alpha})}{\kappa_{2}\tau_{T}},$ (239)

where $\phi_{1}$ denotes the smallest eigenvalue of the Laplacian $-\Delta$
for the given homogeneous boundary condition. Consequently, the original
assumption of $0\leq\tau_{\alpha}<\tau_{T}<\tau_{q}$ should be replaced with
Eq. (239) for a first order approximation. Similarly to the DPL model, the
situation changes with the order of expansion. If one takes the second-order
term as well in $\tau_{q}$, then a fourth-order term appears as

$\displaystyle\frac{\tau_{q}^{2}}{2}\rho c_{v}\partial_{t}^{4}T+\tau_{q}\rho
c_{v}\partial_{t}^{3}T+\rho c_{v}\partial_{t}^{2}T=\kappa_{1}\Delta
T+\hat{\tau}_{\alpha}\partial_{t}\Delta
T+\kappa_{2}\tau_{T}\partial_{t}^{2}\Delta T,$ (240)

for which the solution will be always unstable if
$\tau_{q}\hat{\tau}_{\alpha}>2\kappa_{2}\tau_{T}$. If
$\tau_{q}\hat{\tau}_{\alpha}<2\kappa_{2}\tau_{T}$, but
$\hat{\tau}_{\alpha}>\kappa_{1}\tau_{q}$, then unstable solutions may also
exist. In order to ensure stability, the necessary condition is

$\displaystyle\kappa_{1}\tau_{q}<\hat{\tau}_{\alpha}<\frac{2\kappa_{2}\tau_{T}}{\tau_{q}},$
(241)

so that the stability condition is far from being trivial, and thus the
triple-phase-lag model suffers from the same weaknesses as the DPL approach.
For $\kappa_{1}=0$, the requirements are reduced to Eq. (213) [429]. For
further details on stability, we refer to [430].

To achieve a mechanical coupling, Choudhuri considered the first-order
expansion of (237), and a thermal expansion source term in the energy balance
(similarly to (103)), thus Eq. (238) is extended as well,

$\displaystyle\tau_{q}\rho c_{v}\dddot{T}+\rho c_{v}\ddot{T}+\gamma
T_{0}\ddot{D}+\tau_{q}\gamma T_{0}\dddot{D}=\kappa_{1}\Delta
T+\hat{\tau}_{\alpha}\Delta\dot{T}+\kappa_{2}\tau_{T}\Delta\ddot{T}+\dot{Q}_{v}+\tau_{q}\ddot{Q}_{v},$
(242)

where the partial derivatives are exchanged with material time derivative, and
$D$ denotes the dilatation. Eq. (242) becomes a mathematically and physically
closed system together with the equation of motion, following the formalism of
[428], it reads

$\displaystyle\rho\ddot{\mathbf{u}}=\mu\Delta\mathbf{u}+(\lambda+\mu)\nabla\nabla\cdot\mathbf{u}-\gamma\nabla
T+\rho\mathbf{f},$ (243)

in which $\lambda$ and $\mu$ are the Lame coefficients,
$\gamma=(3\lambda+2\mu)\chi$ (recalling that $\chi$ is the thermal expansion
coefficient), and $\mathbf{u}$ is the displacement field on which a volumetric
force $\mathbf{f}$ might act as well.

For further details on well-posedness and stability, let us refer to the
papers [431, 432, 429]. We note that in [429, 432], further modifications are
studied based on the work of Gurtin et al. [433, 434], namely, the triple-
phase-lag model (237) is combined with a special two-temperature approach. The
model is based on the assumption that there are two distinct temperatures, one
conductive ($\theta$) and one thermodynamic ($T$), connected through

$\displaystyle T=\theta-a\Delta\theta,\quad a\geq 0,$ (244)

and $T$ is in the energy balance (4), and $\theta$ substitutes $T$ in the
constitutive equation (237). Moreover, analogously, the thermal displacement
$\alpha$ is also exchanged with $v$ in (237), and

$\displaystyle\alpha=v-a\Delta v,\quad a\geq 0$ (245)

assumed to be valid. So that Eqs. (244) and (245) are added to the model,
resulting in a much more complex approach, and due to the Laplacian terms,
these add further higher-order derivatives to the constitutive equation (237).
For the corresponding stability and well-posedness conditions, we refer to
[432] in which various sub-cases are studied.

##  8\. Outlook on further concepts

While in the previous Sections, we systematically went through the models with
thermodynamic origins, here we would like to look towards two further concepts
having no deeper relationship to and direct compatibility with the second law
of thermodynamics. These concepts are the thermomass and fractional
derivatives, both falling outside the systematic structure of evolution
equations, however, it does not exclude the possibility that the resulting
model cannot be analogous with any of the previous ones. In parallel, we also
emphasize that if a model fits into the systematic generalization of Fourier’s
law, it does not guarantee its physical admissibility (such as the DPL
concept), but undoubtedly eases its interpretation, finding its solution and
applications since the resulting model can inherit (at least partially) the
existing methods. Neither concepts of thermomass nor fractional derivatives
can inherit any previous properties due to their initial characteristic
assumptions from which the evolution equations are obtained, and that
attribute distinguishes them from the previous models. Nevertheless, if one
accepts their theoretical limitations, these approaches could still provide
model equations with acceptable outcomes for specific heat conduction
problems.

###  8.1. Thermomass concept.

Based on the work of Nie et al. [435], and Guo et al. [44, 436], the starting
point originates from Einstein’s mass-energy relation $E=mc^{2}$, which is
assumed to be valid for the internal energy, too. In other words, they assume
a thermomass expressed as

$\displaystyle
m_{h}=\frac{E_{h}}{c^{2}},\quad\Rightarrow\quad\rho_{h}=\frac{\rho
c_{v}T}{c^{2}},\quad\mathbf{u}_{h}=\frac{\mathbf{q}}{\rho c_{v}T},$ (246)

hence the thermomass gas is characterized with its density $\rho_{h}$ and the
drift velocity $\mathbf{u}_{h}$. Analogously with phonons, the thermomass gas
is constituted of thermons. Eq. (246) identifies the mass balance of the
thermomass gas directly with the balance of internal energy,

$\displaystyle\frac{\partial\rho_{h}}{\partial
t}+\nabla\cdot(\mathbf{u}_{h}\rho_{h})$ $\displaystyle=Q_{m},$ (247)
$\displaystyle\frac{\partial(\rho_{h}\mathbf{u}_{h})}{\partial
t}+\nabla\cdot(\rho_{h}\mathbf{u}_{h}\mathbf{u}_{h})$ $\displaystyle=-\nabla
P_{h}+\mathbf{f}_{h},$ (248) $\displaystyle\frac{\partial}{\partial
t}\left(\frac{1}{2}\rho_{h}\mathbf{u}_{h}^{2}\right)+\nabla\left(\mathbf{u}_{h}\cdot\frac{1}{2}\rho_{h}\mathbf{u}_{h}^{2}\right)$
$\displaystyle=-\nabla
P_{h}\cdot\mathbf{u}_{h}+\mathbf{f}_{h}\cdot\mathbf{u}_{h}.$ (249)

The pressure-like quantity, $P_{h}$, appearing in the balances of momentum
(248) and kinetic energy (249), is defined as

$\displaystyle P_{h}=\frac{\gamma_{h}\rho}{c^{2}}(c_{v}T)^{2},$ (250)

with $\gamma_{h}$ being the so-called Grüneisen constant, which could be used
to take the thermo-mechanical effects into account as $\gamma_{h}$ depends on
the thermal expansion coefficient and bulk modulus. However, due to the
initial assumption (246), this concept is restricted only to rigid heat
conductors, no coupling is possible, and the mass of the thermomass gas
becomes a constitutive quantity. Moreover, any volumetric heat source results
in a $Q_{m}$ source term of the mass balance (247), i.e., the thermomass
cannot be a conserved quantity here.

####  8.1.1. Flux limiters and heat choking

Fourier’s law can be recovered with the assumption that the so-called
resistance $\mathbf{f}_{h}$ is proportional to the drift velocity as
$\mathbf{f}_{h}=-\mu\mathbf{u}_{h}$, and the entire left hand side of the
momentum balance (248) is zero, therefore

$\displaystyle\nabla
P_{h}=\mathbf{f}_{h},\quad\Rightarrow\quad\nabla\left(\frac{\gamma_{h}\rho}{c^{2}}(c_{v}T)^{2}\right)=-\mu\frac{\mathbf{q}}{\rho