# Investigating the whirling heat current density in the Guyer–Krumhansl
equation

Mátyás Szücs  Department of Energy Engineering, Faculty of Mechanical
Engineering, Budapest University of Technology and Economics, Műegyetem rkp.
3., H-1111 Budapest, Hungary  Department of Theoretical Physics, HUN-REN
Wigner Research Centre for Physics, Institute for Particle and Nuclear
Physics, Konkoly Thege Miklós St. 29–33, H-1121 Budapest, Hungary  Montavid
Thermodynamic Research Group, Budapest, Hungary  Carmelo Filippo Munafo
Department of Energy Engineering, Faculty of Mechanical Engineering, Budapest
University of Technology and Economics, Műegyetem rkp. 3., H-1111 Budapest,
Hungary  Department of Mathematical and Computer Sciences, Physical Sciences
and Earth Sciences, University of Messina, Viale F. Stagno d’Alcontres 31,
98166, Messina, Italy  Róbert Kovács  Department of Energy Engineering,
Faculty of Mechanical Engineering, Budapest University of Technology and
Economics, Műegyetem rkp. 3., H-1111 Budapest, Hungary  Department of
Theoretical Physics, HUN-REN Wigner Research Centre for Physics, Institute for
Particle and Nuclear Physics, Konkoly Thege Miklós St. 29–33, H-1121 Budapest,
Hungary  Montavid Thermodynamic Research Group, Budapest, Hungary

###### Abstract

Among the numerous heat conduction models, the Guyer–Krumhansl equation has a
special role. Besides its various application possibilities in nanotechnology,
cryotechnology, and even in case of modeling heterogeneous materials, it poses
additional mathematical challenges compared to the Fourier or Cattaneo (a.k.a.
Maxwell–Cattaneo–Vernotte) equations. Furthermore, the Guyer–Krumhansl
equation is the first heat conduction model, which includes the curl of the
heat flux density in the evolution equation. In the present paper, we place
our focus on the consequences of the existence of such whirling heat current
density by solving the two-dimensional Guyer–Krumhansl equation with a space
and time-dependent heat pulse boundary condition. The discretization poses
further challenges in regard to the boundary condition for which we propose a
particular extrapolation method. Furthermore, with the help of the Helmholtz
decomposition, we show the analogy with the linearized acoustics of Newtonian
fluids, which reveals how the heat flux density plays the role of the velocity
field. Our solutions also reveal an unexpected temperature evolution caused by
the whirling heat flux density, namely, the temperature can locally be
decreased for a short time in a case when the curl of the heat flux density
dominates the heat conduction process.

##  1 Introduction

Several heat conduction models and approaches have been developed and tested
in previous decades, such as the Cattaneo [1], Guyer–Krumhansl (GK) [2], two-
temperature [3, 4, 5], and Jeffreys equations [6]. From an engineering point
of view, the GK and Jeffreys equations could be viable alternatives to the
Fourier equation [7]. In some particular situation, also the two-temperature
model can be useful, however, with strict restrictions on the composition,
time scales and measurements [8]. Both models contain Fourier’s law as a
particular case, furthermore, they can model second sound (similarly to the
Cattaneo equation) and provide a useful, effective description for
heterogeneous materials [7] using a continuum approach. Therefore, the
continuum background of GK equation hold broad potential for more and more
advanced practical applications, also including heat foam-based heat
exchangers and thermal storage technologies [9, 10, 11].

The GK equation, however, is a notably more complicated model than the
Jeffreys equation, especially in two or three spatial dimensions. The reason
lies in their background. The Jeffreys equation doubles Fourier’s law and does
not introduce second-order tensors along with further complications
(couplings, isotropic representation, boundary conditions, etc.). On the
contrary, the GK equation originates from phonon hydrodynamics [2, 12, 13],
and therefore, it is a special fluid model, thus allowing the curl of the heat
flux. This phenomenon is particularly interesting for superfluids [14, 15,
16], and also makes us revisit our expectations about how it can appear in the
temperature history, even for solids in their effective description. In our
work, we place the emphasis on the transient evolution which highlights
particular temperature decrease effects occurring locally due to the heat
current vorticity in agreement with the findings of [17]. This is apparent
only in a transient setting, the time evolution of the heat current vorticity
induces the local temperature decrease relative to the initial state.
Therefore, these observations remain hidden in the work of Beardo et al. [18]
since they applied the GK equation on a stationary problem.

The phonon hydrodynamic description is mainly characteristic of Rational
Extended Thermodynamics, constraining the coefficients and possible
nonlinearities (state dependence) [12]. However, it is possible to derive the
GK equation in a continuum framework using the internal variable approach. We
wish to introduce an additional interesting aspect compared to the work of
Fülöp and Ván [19] in regard to the constitutive properties of the internal
variable [20]. The continuum background offers more flexibility but keeps the
structure of the GK equation. Consequently, as the model becomes free from the
limitations of phonon hydrodynamics, the coefficients are merely restricted by
the second law of thermodynamics. Thus, the GK equation can be used in the
effective description of heterogeneous materials [7] and still can inherit the
coefficients from phonon hydrodynamics. However, we wish to emphasize that we
do not want to delve into the general details of two-dimensional materials and
the related phonon hydrodynamic properties. In this respect, we want to refer
to the recent work of Shang et al. [21] in which they also conduct a similar
investigation as ours, but they slightly modify the GK equation based on
particular phonon hydrodynamic assumptions being characteristic at nanoscale
for two-dimensional materials. Since phonon hydrodynamics strictly restrict
the applicable parameters, the effects of the whirling heat current density
remain hidden, however, the flexibility given the present continuum approach
allows to find further insights.

In the present paper, we show the continuum derivation of the GK equation
together with its complete isotropic representation – which, to the best of
our knowledge, was first derived by Ván [22] –, highlighting the differences
compared to phonon hydrodynamics. Longitudinal and transversal modes are also
discussed, highlighting the role of the rotational part of the heat flux field
and providing further insights into the structure of the GK equation.
Additionally, we solve a two-dimensional problem with space and time-dependent
heat pulse boundary condition, using a staggered scheme. Our particular
interest is to show the consequences of space-dependent boundaries and how
they affect their numerical realization. We evaluate the findings in relation
to the curl of the heat flux, the necessary boundary conditions, and how such
a complex model can relate to the Fourier equation.

##  2 The Guyer–Krumhansl equation within the frame of internal variables

Considering a rigid, homogeneous, and isotropic heat conductor at rest w.r.t.
a given reference frame, the balance of internal energy with neglected
mechanical interactions and volumetric heat sources reads as

$\displaystyle\varrho\partial_{t}e=-\nabla\cdot\mathbf{q},$ (1)

where $\varrho$, $e$ and $\mathbf{q}$ are the (mass) density – which is
constant due to the assumption of rigidity –, the (mass) specific internal
energy and the heat current density, respectively, $\partial_{t}$ is the
partial time derivative and $\nabla$ is the nabla operator. Let $\bm{\varphi}$
denote a physical quantity with arbitrary tensorial order, i.e., its
components w.r.t. the corresponding basis built from Cartesian base vectors
$\mathbf{e}_{j},\ j=1,2,3$ are denoted as $\varphi_{i_{1}\dots i_{N}},\
i_{1},\dots,i_{N}=1,2,3$. In accordance with the usual notation and usage of
the nabla operator in continuum mechanics, right gradient, right divergence
and right curl are defined as

$\displaystyle\operatorname{grad}_{\rm R}\bm{\varphi}$
$\displaystyle:=\bm{\varphi}\otimes\nabla=\partial_{j}\varphi_{i_{1}\dots
i_{N}}\mathbf{e}_{i_{1}}\otimes\dots\otimes\mathbf{e}_{i_{N}}\otimes\mathbf{e}_{j},$
(2) $\displaystyle\operatorname{div}_{\rm R}\bm{\varphi}$
$\displaystyle:=\bm{\varphi}\cdot\nabla=\partial_{j}\varphi_{i_{1}\dots
i_{N-1}j}\mathbf{e}_{i_{1}}\otimes\dots\otimes\mathbf{e}_{i_{N-1}},$ (3)
$\displaystyle\operatorname{curl}_{\rm R}\bm{\varphi}$
$\displaystyle:=\bm{\varphi}\times\nabla=\epsilon_{i_{N}jk}\partial_{j}\varphi_{i_{1}\dots
i_{N}}\mathbf{e}_{i_{1}}\otimes\dots\otimes\mathbf{e}_{i_{N-1}}\otimes\mathbf{e}_{k},$
(4)

where Einstein summation is applied and $\epsilon_{ijk}$ is the Levi-Civita
permutation symbol. Especially, for an arbitrary scalar field $a$

$\displaystyle a\otimes\nabla$ $\displaystyle=\nabla\otimes a=:\nabla a$ (5)

and, divergence and rotation is not interpreted. For an arbitrary vector field
$\mathbf{v}$

$\displaystyle\mathbf{v}\otimes\nabla$
$\displaystyle=\left(\nabla\otimes\mathbf{v}\right)^{\rm T},$ (6)
$\displaystyle\mathbf{v}\cdot\nabla$ $\displaystyle=\nabla\cdot\mathbf{v},$
(7) $\displaystyle\mathbf{v}\times\nabla$
$\displaystyle=-\nabla\times\mathbf{v},$ (8)

where T denotes the transpose of a second order tensor.

Assuming that the thermodynamical state space is spanned by the specific
internal energy and an internal variable denoted by $\bm{\zeta}$, then the
Gibbs relation for the specific entropy function $s$ reads as

$\displaystyle{\rm d}s=\frac{1}{T}{\rm d}e-\mathbf{z}\cdot{\rm d}\bm{\zeta},$
(9)

i.e., the partial derivatives of entropy w.r.t. variables $e$ and $\bm{\zeta}$
are

$\displaystyle\frac{\partial s}{\partial e}$ $\displaystyle=\frac{1}{T},$
$\displaystyle\frac{\partial s}{\partial\bm{\zeta}}$
$\displaystyle=-\mathbf{z},$ (10)

where $T$ denotes the temperature and $\mathbf{z}$ is the entropy conjugate of
the internal variable. Let us assume that entropy current density $\mathbf{J}$
is also generalized through a second order tensor, called Nyíri multiplier
[23] and denoted by $\mathbf{C}$, i.e.,

$\displaystyle\mathbf{J}=\left(\frac{1}{T}\mathbf{1}+\mathbf{C}\right)\mathbf{q},$
(11)

where $\mathbf{1}$ is the identity tensor. Therefore, entropy production rate
density $\Sigma_{S}$ is calculated through the balance of entropy

$\displaystyle 0\leq\Sigma_{S}$
$\displaystyle=\varrho\partial_{t}s+\nabla\cdot\mathbf{J}\stackrel{{\scriptstyle(\ref{eq:Gibbs})-(\ref{eq:Nyiri})}}{{=}}\frac{1}{T}\varrho\partial_{t}e-\mathbf{z}\cdot\varrho\partial_{t}\bm{\zeta}+\left[\left(\frac{1}{T}\mathbf{1}+\mathbf{C}\right)\cdot\nabla\right]\cdot\mathbf{q}+\left(\frac{1}{T}\mathbf{1}+\mathbf{C}\right):\left(\mathbf{q}\otimes\nabla\right)$
$\displaystyle\stackrel{{\scriptstyle(\ref{eq:bal-e})}}{{=}}-\frac{1}{T}\nabla\cdot\mathbf{q}-\mathbf{z}\cdot\varrho\partial_{t}\bm{\zeta}+\left[\left(\frac{1}{T}\mathbf{1}+\mathbf{C}\right)\cdot\nabla\right]\cdot\mathbf{q}+\frac{1}{T}\nabla\cdot\mathbf{q}+\mathbf{C}:\left(\mathbf{q}\otimes\nabla\right)$
$\displaystyle=-\mathbf{z}\cdot\varrho\partial_{t}\bm{\zeta}+\left[\left(\frac{1}{T}\mathbf{1}+\mathbf{C}\right)\cdot\nabla\right]\cdot\mathbf{q}+\mathbf{C}:\left(\mathbf{q}\otimes\nabla\right).$
(12)

According to classical irreversible thermodynamics [24, 25], entropy
production rate density is a quadratic expression of thermodynamical fluxes
and forces, among which functional relations are prescibed, i.e., fluxes
depend on the forces. However, to differentiate between fluxes and forces in
(12) is complicated, these correspondences are rather formal. Nevertheless,
the proven Onsagerian method works, we choose and group the quantities
according to Table 1.

To be determined (‘fluxes’) | $\varrho\partial_{t}\bm{\zeta}$ | $\mathbf{q}$ | $\mathbf{C}$
---|---|---|---
In relation with (‘forces’) | $\mathbf{z}$ | $\left(\frac{1}{T}\mathbf{1}+\mathbf{C}\right)\cdot\nabla$ | $\mathbf{q}\otimes\nabla$
Table 1: Thermodynamic ‘fluxes’ and ‘forces’.

Assuming isotropic material, due to Curie’s principle the different tensorial
orders and characters do not couple, hence positive semi-definiteness of (12)
is ensured via the linear Onsagerian equations

$\displaystyle\varrho\partial_{t}\bm{\zeta}$
$\displaystyle=-l_{11}\mathbf{z}+l_{12}\left(\frac{1}{T}\mathbf{1}+\mathbf{C}\right)\cdot\nabla,$
(13) $\displaystyle\mathbf{q}$
$\displaystyle=-l_{21}\mathbf{z}+l_{22}\left(\frac{1}{T}\mathbf{1}+\mathbf{C}\right)\cdot\nabla,$
(14) $\displaystyle\mathbf{C}$
$\displaystyle=\mathbb{L}:\left(\mathbf{q}\otimes\nabla\right)$ (15)

with the criteria on the scalar coefficients

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

and the fourth order constitutive tensor $\mathbb{L}$ with appropriate
conditions. Since isotropy of the material is assumed (15) can be given in the
isotropic decomposition as

$\displaystyle\mathbf{C}$ $\displaystyle=\frac{L^{\rm s}-L^{\rm
d}}{3}\left(\nabla\cdot\mathbf{q}\right)\mathbf{1}+\frac{L^{\rm
d}}{2}\left(\mathbf{q}\otimes\nabla+\nabla\otimes\mathbf{q}\right)+\frac{L^{\rm
A}}{2}\bm{\epsilon}:\left(\mathbf{q}\times\nabla\right),$ (17)

with the Levi–Civita tensor $\bm{\epsilon}$ and the scalar coefficients

$\displaystyle L^{\rm s}$ $\displaystyle\geq 0,$ $\displaystyle L^{\rm d}$
$\displaystyle\geq 0,$ $\displaystyle L^{\rm A}$ $\displaystyle\geq 0,$ (18)

which map the spherical, symmetric traceless (i.e., deviatoric), and
antisymmetric parts of the gradient of heat current density, respectively. In
general, coefficients $l_{11}$, $l_{12}$, $l_{21}$, $l_{22}$, $L^{\rm s}$,
$L^{\rm d}$ and $L^{\rm A}$ are state-dependent parameters, however, for
simplicity, we treat them as constants. Equations (13) and (14) can be
reformulated as

$\displaystyle\varrho\partial_{t}\bm{\zeta}$
$\displaystyle=-l_{11}\mathbf{z}+l^{\rm
S}\left(\frac{1}{T}\mathbf{1}+\mathbf{C}\right)\cdot\nabla+l^{\rm
A}\left(\frac{1}{T}\mathbf{1}+\mathbf{C}\right)\cdot\nabla,$ (19)
$\displaystyle\mathbf{q}$ $\displaystyle=-l^{\rm S}\mathbf{z}+l^{\rm
A}\mathbf{z}+l_{22}\left(\frac{1}{T}\mathbf{1}+\mathbf{C}\right)\cdot\nabla$
(20)

with $l^{\rm S}=\frac{1}{2}\left(l_{12}+l_{21}\right)$ and $l^{\rm
A}=\frac{1}{2}\left(l_{12}-l_{21}\right)$, then entropy production rate
density can be given as

$\displaystyle
0\leq\Sigma_{S}=\begin{pmatrix}-\mathbf{z}&\frac{1}{T}\mathbf{1}+\mathbf{C}\end{pmatrix}\cdot\begin{pmatrix}l_{11}&l^{\rm
S}\\\ l^{\rm S}&l_{22}\end{pmatrix}\begin{pmatrix}-\mathbf{z}\\\
\frac{1}{T}\mathbf{1}+\mathbf{C}\end{pmatrix}+\mathbf{C}:\left(\mathbf{q}\otimes\nabla\right),$

therefore, the terms in (19) and (20) with coefficients $l^{\rm A}$ do not
increase entropy. Let us now assume that $l^{\rm A}=-1$, $l^{\rm S}=0$ and
$l_{22}=0$, then one obtains for (19) and (20)

$\displaystyle\varrho\partial_{t}\bm{\zeta}$
$\displaystyle=-\left(\frac{1}{T}\mathbf{1}+\mathbf{C}\right)\cdot\nabla-
l_{11}\mathbf{z},$ (21) $\displaystyle\mathbf{q}$ $\displaystyle=-\mathbf{z}.$
(22)

According to the latter equation, the entropy-conjugated internal variable
$-\mathbf{z}$ is identified with the heat current density. Furthermore,
assuming the linear equation of state

$\displaystyle\mathbf{z}=m\bm{\zeta}$ (23)

with a constant $m$ – which is positive to ensure concavity of specific
entropy –, the internal variable can be eliminated, hence (21) and (22)
reduces to the single equation

$\displaystyle-\frac{\varrho}{m}\partial_{t}\mathbf{q}$
$\displaystyle=-\left(\frac{1}{T}\mathbf{1}+\mathbf{C}\right)\cdot\nabla+l_{11}\mathbf{q}.$
(24)

Finally, replacing (17) into (24), we obtain the Guyer–Krumhansl equation

$\displaystyle\tau\partial_{t}\mathbf{q}+\mathbf{q}=-\lambda\nabla
T+\eta_{1}\Delta\mathbf{q}+\eta_{2}\nabla\left(\nabla\cdot\mathbf{q}\right)$
(25)

with the coefficients defined as

$\displaystyle\lambda$ $\displaystyle:=\frac{1}{l_{11}T^{2}}\geq 0,$
$\displaystyle\tau$ $\displaystyle:=\frac{\varrho}{ml_{11}}\geq 0,$
$\displaystyle\eta_{1}$ $\displaystyle:=\frac{L^{\rm d}+L^{\rm
A}}{2l_{11}}\geq 0,$ $\displaystyle\eta_{2}$ $\displaystyle:=\frac{2L^{\rm
s}+L^{\rm d}-3L^{\rm A}}{6l_{11}}.$ (26)

The balance of internal energy (1), the Guyer–Krumhansl equation together with
the (thermostatic caloric) equation of state

$\displaystyle e=cT,$ (27)

where $c$ is the specific heat capacity of the material, forming a closed
system of equations, which, with appropriate initial and boundary conditions,
can be solved.

We must observe that $\eta_{1}$ and $\eta_{2}$ are linearly independent
coefficients and are merely restricted by the second law of thermodynamics.
Therefore, they are positive semidefinite and free from the restrictions given
by the phonon hydrodynamic background, i.e., $\eta_{1}$ is not identical to
the square of the mean free path, and the ratio of $\eta_{2}$ to $\eta_{1}$ is
not necessarily 2 either. The present continuum background provides a notably
more flexible adjustment for the coefficients, either by means of an
experiment such as [7] or inheriting the phonon hydrodynamic approach [12,
13].

Furthermore, we wish to call attention to the coefficients (26) in which
relations their functional relationships become apparent. It is clear that if
$\lambda=\lambda(T)$ holds (e.g, a linear or exponential one), then
$l_{11}=l_{11}(T)$ can be immediately given and that $T$-dependence is
inherited in all the other coefficients. The different parameters can
optionally adjust the necessary $T$-dependence. Additionally, if such
nonlinearities are required, then equations (17), (19) and (20) [or instead of
these latter ones (21) and (22)] must be the starting point in order to take
into account the additional contributions correctly. For instance, due to an
assumed temperature dependence in the coefficients $L^{\rm s}$, $L^{\rm d}$
and $L^{\rm A}$ a direct multiplicative coupling between temperature gradient
and gradient of heat current density emerges.

Comparing our derivation to the work of Fülöp and Ván [19], we have introduced
two notable modifications here.

First, we distinguished between the internal variable and its entropy
conjugate, revealing the need for an additional equation of state. In other
words, the particular form of specific entropy

$\displaystyle
s(e,\bm{\zeta})=\hat{s}(e)-\frac{m}{2}\bm{\zeta}\cdot\bm{\zeta}$ (28)

with $m>0$ constant immediately implicitly imposes equation (23), hence, in
this linear case with constant coefficients the two derivations are
equivalent. When the equation of state (23) is assumed to be linear, then the
internal variable itself can be identified with the heat current density, as
done, for example, in [26]. Futhermore, distinction between the internal
variable and its conjugate enables the compatibility with the GENERIC (General
Equation for the Non-Equilibrium Reversible–Irreversible Coupling) framework
[27]. Finally we note that if equation of state (23) is beyond a linear
relationship, then elimination of the internal variable and its conjugate can
be complicated or impossible, therefore, in this case deeper knowledge (e.g,
microscopic or mesoscopic interpretations) about these variables are required.

Second, we interpret the current multiplier as a separate, relaxed state
variable in accordance with [27], hence a direct interpretation of Müller’s
$\mathbf{K}$-vector as $\mathbf{K}=\mathbf{C}\mathbf{q}$ is obtained [28, 29].

##  3 Longitudinal and transversal heat propagation

Via the vector Laplacian identity

$\displaystyle\Delta\mathbf{q}=\nabla\left(\nabla\cdot\mathbf{q}\right)-\nabla\times\nabla\times\mathbf{q}$
(29)

the Guyer–Krumhansl constitutive equation (25) can also be given as

$\displaystyle\tau\partial_{t}\mathbf{q}+\mathbf{q}=-\lambda\nabla
T+\left(\eta_{1}+\eta_{2}\right)\nabla\left(\nabla\cdot\mathbf{q}\right)-\eta_{1}\nabla\times\nabla\times\mathbf{q},$
(30)

therefore, theoretically, in a GK-type heat conductor transversal heat
propagation can also be observed. Via the Helmholtz decomposition, the heat
current density can be uniquely – up to a space-independent but arbitrary
time-dependent vector function – given as a sum of irrotational (i.e., curl-
free, $\mathbf{q}^{*}$) and solenoidal (i.e., divergence-free,
$\mathbf{q}^{\circ}$) vector fields, i.e.,

$\displaystyle\mathbf{q}$ $\displaystyle=\mathbf{q}^{*}+\mathbf{q}^{\circ},$
with $\displaystyle\nabla\cdot\mathbf{q}$
$\displaystyle=\nabla\cdot\mathbf{q}^{*}$ and
$\displaystyle\nabla\times\mathbf{q}$
$\displaystyle=\nabla\times\mathbf{q}^{\circ}.$ (31)

Curl-free and divergence-free components of the vector field are usually
referred to as longitudinal and transversal components, respectively. Since
$\nabla T$ is curl-free, the governing equations [together with (27)] of a GK-
type heat conductor can be reformulated as

$\displaystyle\varrho c\partial_{t}T$
$\displaystyle=-\nabla\cdot\mathbf{q}^{*},$ (32)
$\displaystyle\tau\partial_{t}\mathbf{q}^{*}+\mathbf{q}^{*}$
$\displaystyle=-\lambda\nabla
T+\left(\eta_{1}+\eta_{2}\right)\Delta\mathbf{q}^{*},$ (33)
$\displaystyle\tau\partial_{t}\mathbf{q}^{\circ}+\mathbf{q}^{\circ}$
$\displaystyle=\eta_{1}\Delta\mathbf{q}^{\circ},$ (34)

which decomposition reveals the time evolutions of the longitudinal and
transversal components of heat current density. Here, the evolution of
$\mathbf{q}^{\circ}$ is decoupled from evolutions of $\mathbf{q}^{*}$ and $T$,
thus $\mathbf{q}^{\circ}$ can only be introduced through a spatially-dependent
boundary condition or by a particular initial condition. A non-homogeneous
temperature field alone can not induce non-zero $\mathbf{q}^{\circ}$. Since
boundary conditions on the heat current density defines $\mathbf{q}$ itself,
it is inevitable to introduce disturbances into both parts of $\mathbf{q}$,
but unique separation of these on the boundary can be complicated.

Later on, we refer on $\eta_{1}+\eta_{2}$ as the longitudinal GK coefficient
and on $\eta_{1}$ as the transversal GK coefficient. Based on our previous 1
spatial dimensional studies of the GK equation [32], the longitudinal GK
coefficient (denoted usually in 1 spatial dimension with $\kappa^{2}$) has to
be positive semi-definite, hence the seemingly indefinite $\eta_{2}$ parameter
in (26) is constraint by

$\displaystyle\eta_{2}\geq-\eta_{1}.$ (35)

An additional interesting property is related to the Fourier resonance
condition, i.e., the temperature history given by the GK equation is identical
with the Fourier’s one. Namely, replacing the gradient of (32) into the
partial time derivative of (33) and taking the advantage of commutation of
$\nabla$ and $\partial_{t}$, one obtains

$\displaystyle\tau\partial_{t}\left(\partial_{t}\mathbf{q}^{*}-\frac{\eta_{1}+\eta_{2}}{\tau}\Delta\mathbf{q}^{*}\right)+\left(\partial_{t}\mathbf{q}^{*}-\frac{\lambda}{\varrho
c}\Delta\mathbf{q}^{*}\right)$ $\displaystyle=0,$ (36)

which is the sum of a Fourier heat conduction equation and the partial time
derivative of a slightly modified Fourier heat conduction equation. If the
additional time scale becomes identical given by the thermal diffusivity,
i.e.,

$\displaystyle\frac{\eta_{1}+\eta_{2}}{\tau}=\frac{\lambda}{\varrho c},$ (37)

Fourier resonance occurs. The resonance condition (37) is the same as obtained
in the one-dimensional case [32], however, in the three-dimensional setting,
Fourier resonance appears only in the longitudinal direction, the transversal
contribution may distort this behaviour. When

$\displaystyle\frac{\left(\eta_{1}+\eta_{2}\right)/\tau}{\lambda/(\varrho
c)}>1,$ (38)

then over-diffusive solutions are obtained, while in the under-damped opposite
case attenuated wave-like propagation of the temperature field is observable.
From an engineering point of view, Fourier resonance also seems to be a
natural requirement, since the GK equation also consists of the Fourier
equation as a particular case, and thus it is highly advantageous in practice
if the GK equation can reproduce the simpler Fourier solution by the
particular alignment of coefficients, without any modifications on the model.

Let us refer here on the governing equations of linearized acoustics of
Newtonian fluids, which reads – via applying Helmholtz decomposition on the
velocity field $\mathbf{v}$ –

$\displaystyle\partial_{t}\rho$
$\displaystyle=-\bar{\rho}\nabla\cdot\mathbf{v}^{*},$ (39)
$\displaystyle\bar{\rho}\partial_{t}\mathbf{v}^{*}$
$\displaystyle=-\bar{a}_{s}^{2}\nabla\rho+\left(\bar{\eta}_{\rm
Vol}+\frac{4}{3}\bar{\eta}_{\rm Sh}\right)\Delta\mathbf{v}^{*},$ (40)
$\displaystyle\bar{\rho}\partial_{t}\mathbf{v}^{\circ}$
$\displaystyle=\bar{\eta}_{\rm Sh}\Delta\mathbf{v}^{\circ}$ (41)

with density difference $\rho$ measured from density $\bar{\rho}$ of the
unperturbed state, isentropic speed of sound $\bar{a}_{s}$, volume and shear
viscosities $\bar{\eta}_{\rm Vol}$ and $\bar{\eta}_{\rm Sh}$, respectively
(cf. equations (3.100) and (3.109) in [30]). Via the curl of velocity the
vorticity $\bm{\omega}=\nabla\times\mathbf{v}=\nabla\times\mathbf{v}^{\circ}$
is defined, hence (41) can also be written on $\bm{\omega}$ instead of
$\mathbf{v}^{\circ}$. Therefore, Helmholtz decomposition highlights that in
linear approximation of acoustics the evolution of $\bm{\omega}$ is decoupled
from $\rho$ and $\mathbf{v}^{*}$, hence resulting that sound is a transversal
wave [cf. (39) and (40)] and vorticity can not be introduced through the
acoustic fields $\rho$ and $\mathbf{v}^{*}$ (nor through the pressure field).
We are dealing with something very similar in case of GK equations (32), (33)
and (34).

##  4 A staggered grid finite difference method demonstrated on the heat pulse
experiment in three spatial dimensions

We aim to numerically model the heat pulse experiment in which a single short
pulse thermally excites the sample. We take into account both the spatial and
time dependence of the pulse, serving as an outstanding example to demonstrate
the role of boundary conditions due to the appearance of in-plane derivatives.

In order to ease the discretization of the second-order derivative of
$\mathbf{q}$, it is advantageous to introduce the gradient of heat current
density as an auxiliary variable

$\displaystyle\mathbf{Q}:=\mathbf{q}\otimes\nabla,$ (42)

hence, a system of first-order equations has to be solved. Additionally, that
$\mathbf{Q}$ feels natural to introduce in the case of the GK equation as it
strongly resembles the current multiplier $\mathbf{C}$. After that
reformulation, we obtain

$\displaystyle\tau\partial_{t}\mathbf{q}+\mathbf{q}=-\lambda\nabla
T+\eta_{1}\mathbf{Q}\cdot\nabla+\eta_{2}\nabla\operatorname{tr}\mathbf{Q}$
(43)

for which it is crucial to correctly discretize $\mathbf{Q}$ on the boundary.

We wish to solve the GK equation in Cartesian coordinate system for a
rectangular domain with the size of $X\times Y\times Z$. We assume that the
heat pulse excites the entire $Z$ direction uniformly, therefore we can reduce
the problem to two-dimensional in $X$ and $Y$. Furthermore, we consider
symmetry at $y=0$, so that $-\frac{Y}{2}\leq y\leq\frac{Y}{2}$, and thus we
deal only half of the rectangle. Consequently, the boundary conditions are

$\displaystyle q_{x}(t,x=0,y)$ $\displaystyle=\begin{cases}Q_{\rm
P,Z}\frac{1}{\tau_{\rm P}Y_{\rm P}}\left[1-\cos\left(2\pi\frac{t}{\tau_{\rm
P}}\right)\right]\left[1+\cos\left(2\pi\frac{y}{Y_{\rm
P}}\right)\right]&\text{if\ }0\leq t\leq\tau_{\rm P}\text{\ and\ }0\leq
y\leq\frac{Y_{\rm P}}{2}\leq\frac{Y}{2}\\\ 0&\text{otherwise}\end{cases},$
(44) $\displaystyle q_{x}(t,x=X,y)$ $\displaystyle=0,$ (45) $\displaystyle
q_{y}(t,x,y=0)$ $\displaystyle=0,$ (46) $\displaystyle
q_{y}(t,x,y=+\frac{Y}{2})$ $\displaystyle=0.$ (47)

where $Q_{\rm P,Z}$ is the $Z$-length specific amount of heat introduced
during the heat pulse, measured in $\rm\frac{J}{m}$. The initial condition
describes equilibrium with homogeneous temperature distribution
$T(t=0,x,y)=T_{0}$, hence $\mathbf{q}(t=0,x,y)=\mathbf{0}$, and
$\mathbf{Q}(t,x,y)=\mathbf{0}$.

We wish to transform the GK equation to a non-dimensional one using the
following characteristic scales. We choose $X$ for the length scale and
$\frac{X^{2}}{\alpha}$ for the time scale using the thermal diffusivity
$\alpha:=\frac{\lambda}{\varrho c}$ (this leads to the usual Fourier number).
Therefore, we obtain the following non-dimensional variables

$\displaystyle\hat{t}$ $\displaystyle:=\frac{t}{\frac{X^{2}}{\alpha}},$
$\displaystyle\hat{x}$ $\displaystyle:=\frac{x}{X},$ $\displaystyle\hat{y}$
$\displaystyle:=\frac{y}{X},$ (48)

and derivatives

$\displaystyle\partial_{\hat{t}}$
$\displaystyle=\frac{X^{2}}{\alpha}\partial_{t},$
$\displaystyle\partial_{\hat{x}}$ $\displaystyle=X\partial_{x},$
$\displaystyle\partial_{\hat{y}}$ $\displaystyle=X\partial_{y}.$ (49)

The non-dimensional fields read

$\displaystyle\hat{q}_{i}:=\frac{q_{i}}{\frac{\alpha Q_{\rm
P,Z}}{X^{3}R_{Y}}},\quad\hat{Q}_{ij}=\frac{Q_{ij}}{\frac{\alpha Q_{\rm
P,Z}}{X^{2}R_{Y}}},\quad\hat{T}:=\frac{T-T_{0}}{T_{\rm
max}-T_{0}}=\frac{\varrho cX^{2}R_{Y}\left(T-T_{0}\right)}{Q_{\rm P,Z}}$ (50)

in which $i,j=x,y$, $R_{Y}=\frac{Y}{X}$ and $T_{\rm max}$ is calculated via
integrating (1) [together with (27)] on the whole sample volume in time from
the initial homogeneous temperature state $T_{0}$ to the final homogeneous
temperature state $T_{\rm max}$.

Let us summarize the complete system of non-dimensional equations,

$\displaystyle\partial_{\hat{t}}\hat{T}$
$\displaystyle=-\left(\partial_{\hat{x}}\hat{q}_{x}+\partial_{\hat{y}}\hat{q}_{y}\right),$
(51) $\displaystyle\hat{\tau}\partial_{\hat{t}}\hat{q}_{x}+\hat{q}_{x}$
$\displaystyle=-\partial_{\hat{x}}\hat{T}+\left(\hat{\eta}_{1}+\hat{\eta}_{2}\right)\partial_{\hat{x}}\hat{Q}_{xx}+\hat{\eta}_{1}\partial_{\hat{y}}\hat{Q}_{xy}+\hat{\eta}_{2}\partial_{\hat{x}}\hat{Q}_{yy},$
(52) $\displaystyle\hat{\tau}\partial_{\hat{t}}\hat{q}_{y}+\hat{q}_{y}$
$\displaystyle=-\partial_{\hat{y}}\hat{T}+\left(\hat{\eta}_{1}+\hat{\eta}_{2}\right)\partial_{\hat{y}}\hat{Q}_{yy}+\hat{\eta}_{1}\partial_{\hat{x}}\hat{Q}_{yx}+\hat{\eta}_{2}\partial_{\hat{y}}\hat{Q}_{xx},$
(53) $\displaystyle\hat{Q}_{xx}$
$\displaystyle=\partial_{\hat{x}}\hat{q}_{x},$ (54)
$\displaystyle\hat{Q}_{xy}$ $\displaystyle=\partial_{\hat{y}}\hat{q}_{x},$
(55) $\displaystyle\hat{Q}_{yx}$
$\displaystyle=\partial_{\hat{x}}\hat{q}_{y},$ (56)
$\displaystyle\hat{Q}_{yy}$ $\displaystyle=\partial_{\hat{y}}\hat{q}_{y}$ (57)

where

$\displaystyle\hat{\tau}$ $\displaystyle=\frac{\tau}{\frac{X^{2}}{\alpha}},$
$\displaystyle\hat{\eta}_{1}$ $\displaystyle=\frac{\eta_{1}}{X^{2}},$
$\displaystyle\hat{\eta}_{2}$ $\displaystyle=\frac{\eta_{2}}{X^{2}},$ (58)

and the non-dimensional heat pulse boundary conditions reads

$\displaystyle\hat{q}_{x}(\hat{t},\hat{x}=0,\hat{y})$
$\displaystyle=\begin{cases}\frac{R_{Y}}{\hat{\tau}_{\rm P}R_{Y,\rm
P}}\left[1-\cos\left(2\pi\frac{\hat{t}}{\hat{\tau}_{\rm
P}}\right)\right]\left[1+\cos\left(2\pi\frac{\hat{y}}{R_{Y,\rm
P}}\right)\right]&\text{if\ }0\leq\hat{t}\leq\hat{\tau}_{\rm P}\text{\ and\
}0\leq\hat{y}\leq R_{Y,\rm P}\leq R_{Y}\\\ 0&\text{otherwise}\end{cases},$
(59)

where $R_{Y,\rm P}=\frac{Y_{\rm P}}{X}$ stands.

Spatial discretization is realized via a staggered scheme [31], which is
depicted in Figure 1, and the structure of governing equations (51)–(57)
restricts how each field can be represented on the discrete lattice with
directional equidistant grid points with distance $\Delta\hat{x}$ and
$\Delta\hat{y}$, hence $\hat{x}_{m}=m\Delta\hat{x},\ m=1,\dots,M$,
$\hat{y}_{n}=n\Delta\hat{y},\ n=1,\dots,N$, where $M=\frac{1}{\Delta\hat{x}}$
and $N=\frac{R_{Y}}{2\Delta\hat{y}}$. The investigated time interval is also
discretized through equidistance time steps $\Delta\hat{t}$, i.e.,
$\hat{t}^{j}=j\Delta\hat{t},\ j=1,\dots,J$. Therefore, the approximated value
of a function $f$ in the discrete time and space coordinates
$\left(\hat{t}^{j},\hat{x}_{m},\hat{y}_{n}\right)$ is denoted by
$f^{j}_{m,n}$. Temperature, as a state variable characterizing homogeneously
one discrete cell, is placed in the middle of the cell, while heat current
density characterizing fluxes through the boundaries of the cell, therefore,
the corresponding normal components of the heat current density are placed on
the boundaries of the cell in line with the discrete temperature values.
Discretization of $\mathbf{Q}$ follows directly from the discrete values of
heat current density and equations (54)–(57), consequently, its diagonal
elements are also in the middle of the cell, but its off-diagonal elements are
placed in the corners of the cell. Since we apply $\mathbf{q}$-boundaries,
only complete cells are used to discretize the entire spatial domain.

Figure 1: Discretization of the 2D GK equation, showing the staggered grid,
and how $\mathbf{Q}$ points are allocated and extrapolated from the bulk.
Additionally, the lower left figure demonstrated the extrapolation we use to
determine the unknown off-diagonal elements.

Furthermore, let us note that according to Figure 1, one needs to prescribe
$Q_{xy}$ and $Q_{yx}$ on the boundaries, in accordance with the
$\mathbf{q}$-boundaries. On each side, one of these off-diagonal quantities
can be determined analytically and represented on the discrete lattice. For
instance, for a given $q_{x}(t,x=0,y)$, $Q_{xy}$ can be determined
immediately, however, $Q_{yx}$ must be extrapolated from the bulk nodes. This
procedure holds for all four boundaries. For the extrapolation, we use
quadratic Lagrange polynomials in order to preserve the sign given by three
bulk points next to each boundary. This extrapolation is schematically
demonstrated in Figure 1 for one setting in which the $m^{th}$ value of
$Q_{yx}$ is calculated based on the bulk point of $m+1$, $m+2$, and $m+3$. We
want to emphasize that any direct definition of the $\mathbf{Q}$ on the
boundaries can significantly distort the physical content of the solution, and
most probably, a different problem is solved than expected in that case. With
our method utilizing Lagrange polynomials, we can avoid the definition of any
additional, unnecessary boundary conditions or introduction of virtual nodes.

For the time derivatives, we choose the simplest forward time stepping method,
the explicit Euler, which is eligible for our aim to investigate the solutions
of the two-dimensional GK equation. Summarizing the numerical scheme built in
accordance of Figure 1:

$\displaystyle\hat{T}_{m+\nicefrac{{1}}{{2}},n+\nicefrac{{1}}{{2}}}^{j+1}$
$\displaystyle=T_{m+\nicefrac{{1}}{{2}},n+\nicefrac{{1}}{{2}}}^{j}+\Delta\hat{t}\left(\frac{\left(\hat{q}_{x}\right)^{j}_{m,n+\nicefrac{{1}}{{2}}}-\left(\hat{q}_{x}\right)^{j}_{m+1,n+\nicefrac{{1}}{{2}}}}{\Delta\hat{x}}+\frac{\left(\hat{q}_{y}\right)^{j}_{m+\nicefrac{{1}}{{2}},n}-\left(\hat{q}_{y}\right)^{j}_{m+\nicefrac{{1}}{{2}},n+1}}{\Delta\hat{y}}\right),$
(60) $\displaystyle\left(\hat{q}_{x}\right)^{j+1}_{m,n+\nicefrac{{1}}{{2}}}$
$\displaystyle=\left(1-\frac{\Delta\hat{t}}{\hat{\tau}}\right)\left(\hat{q}_{x}\right)^{j}_{m,n+\nicefrac{{1}}{{2}}}+\frac{\Delta\hat{t}}{\hat{\tau}}\Bigg{[}\frac{\hat{T}^{j}_{m-\nicefrac{{1}}{{2}},n+\nicefrac{{1}}{{2}}}-\hat{T}^{j}_{m+\nicefrac{{1}}{{2}},n+\nicefrac{{1}}{{2}}}}{\Delta\hat{x}}\mathop{+}$
$\displaystyle\hskip
21.52771pt+\left(\hat{\eta}_{1}+\hat{\eta}_{2}\right)\frac{\left(\hat{Q}_{xx}\right)^{j}_{m+\nicefrac{{1}}{{2}},n+\nicefrac{{1}}{{2}}}-\left(\hat{Q}_{xx}\right)^{j}_{m-\nicefrac{{1}}{{2}},n+\nicefrac{{1}}{{2}}}}{\Delta\hat{x}}+\hat{\eta}_{1}\frac{\left(\hat{Q}_{xy}\right)^{j}_{m,n+1}-\left(\hat{Q}_{xy}\right)^{j}_{m,n}}{\Delta\hat{y}}\mathop{+}$
$\displaystyle\hskip
21.52771pt+\hat{\eta}_{2}\frac{\left(\hat{Q}_{yy}\right)^{j}_{m+\nicefrac{{1}}{{2}},n+\nicefrac{{1}}{{2}}}-\left(\hat{Q}_{yy}\right)^{j}_{m-\nicefrac{{1}}{{2}},n+\nicefrac{{1}}{{2}}}}{\Delta\hat{x}}\Bigg{]},$
(61) $\displaystyle\left(\hat{q}_{y}\right)^{j+1}_{m+\nicefrac{{1}}{{2}},n}$
$\displaystyle=\left(1-\frac{\Delta\hat{t}}{\hat{\tau}}\right)\left(\hat{q}_{y}\right)^{j}_{m+\nicefrac{{1}}{{2}},n}+\frac{\Delta\hat{t}}{\hat{\tau}}\Bigg{[}\frac{\hat{T}^{j}_{m+\nicefrac{{1}}{{2}},n-\nicefrac{{1}}{{2}}}-\hat{T}^{j}_{m+\nicefrac{{1}}{{2}},n+\nicefrac{{1}}{{2}}}}{\Delta\hat{x}}\mathop{+}$
$\displaystyle\hskip
21.52771pt+\left(\hat{\eta}_{1}+\hat{\eta}_{2}\right)\frac{\left(\hat{Q}_{yy}\right)^{j}_{m+\nicefrac{{1}}{{2}},n+\nicefrac{{1}}{{2}}}-\left(\hat{Q}_{yy}\right)^{j}_{m+\nicefrac{{1}}{{2}},n-\nicefrac{{1}}{{2}}}}{\Delta\hat{y}}+\hat{\eta}_{1}\frac{\left(\hat{Q}_{yx}\right)^{j}_{m+1,n}-\left(\hat{Q}_{xy}\right)^{j}_{m,n}}{\Delta\hat{x}}\mathop{+}$
$\displaystyle\hskip
21.52771pt+\hat{\eta}_{2}\frac{\left(\hat{Q}_{xx}\right)^{j}_{m+\nicefrac{{1}}{{2}},n+\nicefrac{{1}}{{2}}}-\left(\hat{Q}_{xx}\right)^{j}_{m+\nicefrac{{1}}{{2}},n-\nicefrac{{1}}{{2}}}}{\Delta\hat{y}}\Bigg{]},$
(62)
$\displaystyle\left(\hat{Q}_{xx}\right)^{j}_{m+\nicefrac{{1}}{{2}},n+\nicefrac{{1}}{{2}}}$
$\displaystyle=\frac{\left(\hat{q}_{x}\right)^{j}_{m+1,n+\nicefrac{{1}}{{2}}}-\left(\hat{q}_{x}\right)^{j}_{m,n+\nicefrac{{1}}{{2}}}}{\Delta\hat{x}},$
(63) $\displaystyle\left(\hat{Q}_{xy}\right)^{j}_{m,n}$
$\displaystyle=\frac{\left(\hat{q}_{x}\right)^{j}_{m,n+\nicefrac{{1}}{{2}}}-\left(\hat{q}_{x}\right)^{j}_{m,n-\nicefrac{{1}}{{2}}}}{\Delta\hat{y}},$
(64) $\displaystyle\left(\hat{Q}_{yx}\right)^{j}_{m,n}$
$\displaystyle=\frac{\left(\hat{q}_{y}\right)^{j}_{m+\nicefrac{{1}}{{2}},n}-\left(\hat{q}_{y}\right)^{j}_{m-\nicefrac{{1}}{{2}},n}}{\Delta\hat{x}},$
(65)
$\displaystyle\left(\hat{Q}_{yy}\right)^{j}_{m+\nicefrac{{1}}{{2}},n+\nicefrac{{1}}{{2}}}$
$\displaystyle=\frac{\left(\hat{q}_{y}\right)^{j}_{m+\nicefrac{{1}}{{2}},n+1}-\left(\hat{q}_{y}\right)^{j}_{m+\nicefrac{{1}}{{2}},n}}{\Delta\hat{y}}.$
(66)

##  5 Numerical results

###  5.1 2D Fourier solutions

Here, let us start with the Fourier heat equation and present the plots we use
to visualize the temperature and heat flux density history in time and space.
We wish to emphasize that we have a different code for the Fourier equation in
which we solve only the Fourier equation and not the simplification of the GK
equation, hence there are no difficulties with the boundary conditions in this
case, i.e., no $\mathbf{Q}$ is used. Figure 2 shows the temperature histories
at specified spatial points. The first column shows the front face, and each
column increases the spatial step by $1/4$, and the last one is related to the
rear face of the sample. Additionally, the first row shows the top side, and
the last one presents the temperature history on the symmetry axis. We use the
same plot concept for the GK equation as well. Furthermore, it is insightful
to check the vector plot for the heat flux field in which we can observe the
effect of space-dependent excitation and compare it immediately with the
corresponding temperature distribution, see Figure 3. Also, Figure 3 shows
that it is straightforward to reproduce the well-known 1D solutions by
applying a spatially homogeneous boundary condition. The curl of the heat flux
field is identically zero. For all the subsequent calculations, we fix
$R_{Y,\rm P}=0.4$ with fixed spatial resolution
($\Delta\hat{x}=\Delta\hat{y}=0.02$), and $\hat{\tau}_{\rm P}=0.01$.

Figure 2: 2D Fourier solution in space and time, also indicating the
corresponding boundary conditions. Figure 3: 2D vector plot of the heat flux
field for the Fourier equation at $\hat{t}=0.1$, using the same color bar for
both figures. Left: spatially inhomogeneous heat pulse boundary condition.
Right: spatially homogeneous heat pulse boundary condition.

###  5.2 2D GK solutions

#### Demonstrating the Fourier resonance

First, we start by showing that the GK equation reproduces the Fourier
solutions when $\hat{\eta}_{1}=0$, and $\hat{\eta}_{2}=\hat{\tau}=0.05$,
Figure 4 shows the difference between the temperature fields. The observed
errors are practically zero. This solution also supports our handling of
$\mathbf{Q}$ on the boundary and the extrapolation method as it reproduces the
Fourier solution in a particular parameter setting. In this parameter setting,
the $\nabla\times\nabla\times\mathbf{q}$ term vanishes in the Eq. (30), thus
the Laplacian of $\mathbf{q}$ remains. Notable differences emerge when we keep
$\hat{\tau}=0.05$, but $\hat{\eta}_{1}=\hat{\eta}_{2}=0.025$, i.e., despite
that $\hat{\eta}_{1}+\hat{\eta}_{2}=\hat{\tau}$, the resonance condition is
violated and the solution significantly differs from the Fourier equation, see
Figure 5 for details. Let us turn our attention to the more interesting and
intriguing solutions when $\hat{\eta}_{1}\neq 0$, and thus the curl of the
heat flux field becomes meaningful.

Figure 4: The difference between the Fourier and GK temperature histories when
$\hat{\eta}_{1}=0$, and $\hat{\eta}_{2}=\hat{\tau}=0.05$. Figure 5: The
difference between the Fourier and GK temperature histories when
$\hat{\eta}_{1}=\hat{\eta}_{2}=0.025$, and $\hat{\tau}=0.05$.

#### Vorticity-free solution

We call the solution to be over-diffusive if
$\hat{\eta}_{1}+\hat{\eta}_{2}>\hat{\tau}$. Let us start with the case when we
keep $\hat{\eta}_{1}=0$ in order to avoid the effects of the nonzero curl of
the heat flux field. It is worth comparing the characteristics of the
temperature history in the middle to the rear side. In the middle, the two
distinct time scales are apparent on the contrary to the rear side (see Figure
6). In a heat pulse experiment using heterogeneous materials, such two
distinct time scales are also visible, and that numerical solution reflects
the size dependence of the observation [32]. This effect might not be apparent
for thicker samples, and it also depends on the material properties. Figure 7
shows the 2D vector plot of the heat flux in which the over-diffusive behavior
remains hidden; the temperature contours (isothermal lines) are slightly
distorted compared to the Fourier case. Furthermore, since $\eta=0$, the curl
of the heat flux field is expected to be zero, and this is reflected by Figure
7, too.

Figure 6: The GK temperature history when $\hat{\eta}_{1}=0$,
$\hat{\eta}_{2}=0.1$, and $\hat{\tau}=0.05$. Figure 7: Left: 2D vector plot of
the heat flux field when $\hat{\eta}_{1}=0$, $\hat{\eta}_{2}=0.1$, and
$\hat{\tau}=0.05$ at $\hat{t}=0.1$. Right: the curl of the heat flux field is
practically zero in accordance with the parameter settings.

#### Solutions with strong vorticity

We turn our attention to the reverse case in which we keep $\hat{\eta}_{2}=0$,
and now let us investigate the effects of the parameter $\hat{\eta}_{1}$, let
it be $\hat{\eta}_{1}=0.075$. Figure 8 presents the temperature history for
the given spatial points. Comparing it to the previous case, we can observe a
notably different behavior. First, near the front face, the temperature
decreases due to the significant curl effects. We emphasize that this is not
identical to reaching a negative absolute temperature. We wish to emphasize
that the apparent negative temperature is relative to the initial temperature
on contrary to the observations of Zhukovsky [33]. However, the temperature
field can exhibit unusual evolution in such a particular parameter setting
since the GK equation is based on a hydrodynamic analogy, and the curl of the
heat flux field can naturally appear. Here, we particularly strengthened this
effect to make it easily observable. This is only meaningful in a two- or
three-dimensional setting. This temperature-decreasing effect disappears soon
and is not observable for any other spatial domains. This is also depicted in
Figure 9. Furthermore, contrary to the previous situations, the curl of the
heat flux field becomes significantly larger (Figure 10).

Figure 8: The GK temperature history when $\hat{\eta}_{1}=0.075$,
$\hat{\eta}_{2}=0$, and $\hat{\tau}=0.05$. Figure 9: Left: 2D vector plot of
the heat flux field when $\hat{\eta}_{1}=0.075$, $\hat{\eta}_{2}=0$, and
$\hat{\tau}=0.05$ at $\hat{t}=0.1$. Right: the contour plot of the temperature
field, highlighting the temperature-decrease effect next to the heat pulse at
$\hat{t}=0.01$. Figure 10: Left: the curl of the heat flux field when
$\hat{\eta}_{1}=0.075$, $\hat{\eta}_{2}=0$, and $\hat{\tau}=0.05$ right after
the heat pulse at $\hat{t}=0.01$. Right: the curl vanishes with time and
becomes practically zero (at $\hat{t}=1$).

#### Devating from the phonon hydrodynamic ratio

We wish to recall that the ratio $\hat{\eta}_{2}/\hat{\eta}_{1}=2$ is fixed in
a phonon hydrodynamic approach, but this does not necessarily hold in a
continuum framework. In order to make this difference apparent, we provide
solutions with respect to $\hat{\eta}_{2}/\hat{\eta}_{1}$. Figure 11 shows the
rear side temperature history for three situations in which
$\hat{\eta}_{2}/\hat{\eta}_{1}=\\{1.5,2,2.5\\}$ with fixed
$\hat{\eta}_{1}=0.05$. This does not show any remarkable properties compared
to the one-dimensional case, for which only the effect of
$\hat{\eta}_{1}+\hat{\eta}_{2}$ is observable, and the increase of the over-
diffusion makes the temperature signal propagation to be faster. However, if
we consider the front-side temperature history in the middle (Figure 12), then
it makes more visible how the ratio of $\hat{\eta}_{2}/\hat{\eta}_{1}$
modifies the solution. Decreasing $\hat{\eta}_{2}/\hat{\eta}_{1}$ amplifies
the temperature-decrease effect near the heat pulse since $\hat{\eta}_{1}$ –
the rotational part – becomes more dominant.

Figure 11: The GK temperature history on the rear side ($\hat{x}=1$,
$\hat{y}=0.25$) when $\hat{\eta}_{1}=0.05$, $\hat{\tau}=0.05$, and
$\hat{\eta}_{2}=\\{0.075,0.1,0.125\\}$, demonstrating the deviation from the
phonon hydrodynamic ratio $\hat{\eta}_{2}/\hat{\eta}_{1}=2$. Figure 12: The GK
temperature history on the front side ($\hat{x}=0$, $\hat{y}=0.25$) when
$\hat{\eta}_{1}=0.05$, $\hat{\tau}=0.05$, and
$\hat{\eta}_{2}=\\{0.1,0.05,0.025\\}$.

#### How does the curl of heat current density behave on the boundary?

Previously, we introduced the auxiliary field $\mathbf{Q}$ in order to ease
the discretization of the second-order spatial derivatives and make it easier
to realize the boundary conditions properly. Since the diagonal elements of
$\mathbf{Q}$ are inside the spatial domain, they are not directly related to
the boundary conditions. However, the off-diagonals are not independent of the
$\mathbf{q}$-boundary. We introduced an extrapolation from the bulk in order
to avoid the definition of incompatible boundary data for the unknown off-
diagonals, and thus avoiding the introduction of any artificial distortion.
Now let us depict the difference between the $Q_{xy}$ and $Q_{yx}$ components
in two cases. In fact, this difference is the only component of the curl of
the in-plane heat current density. In the first one, $\hat{\eta}_{1}=0$ is
considered, hence, the rotational term is zero
($\hat{\eta}_{2}=\hat{\tau}=0.05$). In the second case,
$\hat{\eta}_{1}=\hat{\eta}_{2}=0.05$, and notable differences are expected.
Figure 13 presents their difference, highlighting that $\hat{\eta}_{1}$ indeed
introduces significant changes in the evolution of off-diagonals of
$\mathbf{Q}$, especially near the boundaries, but also notably affecting the
bulk behavior. Figure 14 shows the time evolution in agreement with Eq. (34),
presenting an exponential decay in time.

Figure 13: Left: $\left(\mathbf{q}\times\nabla\right)_{z}$ without the
rotational term at $\hat{t}=0.1$. Right:
$-\left(\mathbf{q}\times\nabla\right)_{z}$ with the rotational term at
$\hat{t}=0.1$. The sign is changed in order to have colors correctly
emphasizing the differences in magnitudes. Figure 14: The time evolution of
$\left(\mathbf{q}\times\nabla\right)_{z}$ including the rotational term on the
front face, indicating that the vorticities are the strongest at the heat
pulse, and quickly attenuate after $\hat{y}=0.2$. Furthermore, shortly after
the heat pulse, the exponential decay is apparent in agreement with Eq. (34).

##  6 Discussion and summary

In the present paper, we revisited the continuum thermodynamic background of
the Guyer–Krumhansl equation. With distinguishing between the internal
variable and its entropy-conjugate, we revealed the need for a further
constitutive relationship, which connect these two together. This additional
possibility remains hidden when a particular extension for the entropy density
is assumed, especially when the internal variable is immediately identified
with the heat flux density. Although we do not investigate its outcomes in
greater detail here, we note that it poses further potential and could have
non-trivial consequences.

Furthermore, we applied a staggered discretization approach to numerically
solve a two-dimensional setting including a space and time-dependent heat flux
boundary condition. The discretization is eased by introducing a second-order
tensory, $\mathbf{Q}$ as an auxiliary quantity, and $\mathbf{Q}$ is also
helpful in the proper realization of boundary conditions. We discussed that
$\mathbf{Q}$ is not independent of the given $\mathbf{q}$ boundary condition,
however, not all components of $\mathbf{Q}$ can be found immediately. In order
to avoid artificial and unnecessary assumptions, we proposed to use a
quadratic Lagrangian extrapolation based on the bulk points to update the
unknown $\mathbf{Q}$ components on the boundary. That approach successfully
reproduced the solutions of the Fourier equation applying the resonance
condition ($\hat{\eta}_{1}=0$, i.e., with vanishing rotational term).

Additionally, the continuum background of the GK equation allowed to adjust
the parameters in a range, which is beyond the validity of phonon
hydrodynamics. In a case when the rotational terms dominate the heat flux
density time evolution, we could observe that the temperature can be
significantly decreased, even under the initial temperature. However, we want
to strongly emphasize that this phenomenon is not identical to obtain negative
temperature, and occurs only locally for a brief period. This is a
characteristic outcome of the whirling heat current density.

Our analysis also revealed that when the GK equation is used as an effective
description of heterogeneous materials, $\hat{\eta}_{1}=0$ is a necessary
choice in order to keep the possibility for Fourier resonance, and avoid
solutions being characteristic only for the hydrodynamic domain, e.g, for low-
temperature or nanoscale problems. However, the continuum background always
offers the possibility to inherit the coefficients from phonon hydrodynamics,
and therefore our numerical approach and findings can be useful even in these
situations.

##  7 Acknowledgement

The authors express their gratitude to Péter Ván for his useful suggestions.
The research of C.F.M. has been carried out under the auspices of GNFM
(National Group of Mathematical-Physics) of INdAM (National Institute of
Advanced Mathematics), through the grant ‘Progetto Giovani’ CUPE53C22001930001
for financial support. Project no. TKP-6-6/PALY-2021 has been implemented with
the support provided by the Ministry of Culture and Innovation of Hungary from
the National Research, Development and Innovation Fund, financed under the
TKP2021-NVA funding scheme. The research was funded by the Sustainable
Development and Technologies National Programme of the Hungarian Academy of
Sciences (FFT NP FTA). This work was partially supported in part by the
Hungarian Scientific Research Fund under Grant agreement FK 134277.

## References

  * [1] C. Cattaneo.  Sur une forme de lequation de la chaleur eliminant le paradoxe dune propagation instantanee.  Comptes Rendus Hebdomadaires Des Seances De L’Academie Des Sciences, 247(4):431–433, 1958.
  * [2] R. A. Guyer and J. A. Krumhansl.  Thermal Conductivity, Second Sound, and Phonon Hydrodynamic Phenomena in Nonmetallic Crystals.  Physical Review, 148:778–788, 1966.
  * [3] S. L. Sobolev.  Heat conduction equation for systems with an inhomogeneous internal structure.  Journal of Engineering Physics and Thermophysics, 66(4):436–440, 1994.
  * [4] S. L. Sobolev.  Local non-equilibrium transport models.  Physics-Uspekhi, 40(10):1043–1053, 1997.
  * [5] S. L. Sobolev.  Nonlocal two-temperature model: Application to heat transport in metals irradiated by ultrashort laser pulses.  International Journal of Heat and Mass Transfer, 94:138–144, 2016\.
  * [6] D. D. Joseph and L. Preziosi.  Heat waves.  Reviews of Modern Physics, 61(1):41–86, 1989.
  * [7] A. Fehér and R. Kovács.  On the dynamic thermal conductivity and diffusivity observed in heat pulse experiments.  Journal of Non-Equilibrium Thermodynamics, 49:161–170, 2024.
  * [8] R. Kovács, A. Fehér, and S. Sobolev.  On the two-temperature description of heterogeneous materials.  International Journal of Heat and Mass Transfer, 194:123021, 2022\.
  * [9] X. Chen, X. Li, X. Xia, C. Sun, and R. Liu.  Thermal storage analysis of a foam-filled PCM heat exchanger subjected to fluctuating flow conditions.  Energy, 216:119259, 2021.
  * [10] A. NematpourKeshteli, M. Iasiello, G. Langella, and N. Bianco.  Enhancing PCMs thermal conductivity: A comparison among porous metal foams, nanoparticles and finned surfaces in triplex tube heat exchangers.  Applied Thermal Engineering, 212:118623, 2022.
  * [11] S. Zhang, L. Pu, S. Mancin, Z. Ma, and L. Xu.  Experimental study on heat transfer characteristics of metal foam/paraffin composite PCMs in large cavities: Effects of material types and heating configurations.  Applied Energy, 325:119790, 2022.
  * [12] I. Müller and T. Ruggeri.  Rational Extended Thermodynamics.  Springer, 1998.
  * [13] W. Dreyer and H. Struchtrup.  Heat pulse experiments revisited.  Continuum Mechanics and Thermodynamics, 5:3–50, 1993.
  * [14] M. S. Mongiovi, D. Jou, and M. Sciacca.  Non-equilibrium thermodynamics, heat transport and thermal waves in laminar and turbulent superfluid helium.  Physics Reports, 2018.
  * [15] M. Sỳkora, M. Pavelka, M. La Mantia, D. Jou, and M. Grmela.  On the relations between large-scale models of superfluid helium-4.  Physics of Fluids, 33(12), 2021.
  * [16] D. Jou, M. S. Mongiovì, and M. Sciacca.  Hydrodynamic equations of anisotropic, polarized and inhomogeneous superfluid vortex tangles.  Physica D: Nonlinear Phenomena, 240(3):249–258, 2011.
  * [17] M. Sỳkora, M. Pavelka, L. Restuccia, and D. Jou.  Multiscale heat transport with inertia and thermal vortices.  Physica Scripta, 98(10):105234, 2023.
  * [18] J. Tur-Prats, M. Gutiérrez-Pérez, J. Bafaluy, J. Camacho, F. X. Alvarez, and A. Beardo.  Microscopic origin of heat vorticity in quasi-ballistic phonon transport.  International Journal of Heat and Mass Transfer, 226:125464, 2024.
  * [19] P. Ván and T. Fülöp.  Universality in heat conduction theory – weakly nonlocal thermodynamics.  Annalen der Physik (Berlin), 524(8):470–478, 2012.
  * [20] A. Berezovski and P. Ván.  Internal Variables in Thermoelasticity.  Springer, 2017.
  * [21] M.-Y. Shang, C. Zhang, Z. Guo, and J.-T. Lü.  Heat vortex in hydrodynamic phonon transport of two-dimensional materials.  Scientific Reports, 10(1):8272, 2020.
  * [22] P. Ván.  Weakly nonlocal irreversible thermodynamics – the Guyer-Krumhansl and the Cahn-Hilliard equations.  Physics Letters A, 290: 88–92, 2001.
  * [23] B. Nyíri.  On the entropy current.  Journal of Non-Equilibrium Thermodynamics, 16(2): 179–186, 1991.
  * [24] S. R. de Groot and P. Mazur.  Non-equilibrium thermodynamics.  Dover Publications, Amsterdam, 1962.
  * [25] I. Gyarmati.  Non-Equilibrium Thermodynamics—Field Theory and Variational Principles.,  Springer-Verlag, Berlin, Heidelberg, New York, 1970.
  * [26] R. Kovács and P. Ván.  Generalized heat conduction in heat pulse experiments.  International Journal of Heat and Mass Transfer, 83:613–620, 2015.
  * [27] M. Szücs, M. Pavelka, R. Kovács, T. Fülöp, P. Ván, and M. Grmela.  A case study of non-Fourier heat condction using internal variables and GENERIC.  Journal of Non-Equilibrium Thermodynamics, 47:31–60, 2022.
  * [28] I. Müller.  A Thermodynamic Theory of Mixtures of Fluids.  Archive for Rational Mechanics and Analysis, 28:1–39, 1968.
  * [29] I. Müller.  Die Kältefunktion, eine universelle Funktion in der Thermodynamik viskoser wärmeleitender Flüssigkeiten.  Archive for Rational Mechanics and Analysis, 40:1–36, 1971.
  * [30] T. D. Rossing (ed.).  Springer Handbook of Acoustics.  Spinger New York, NY, 2007.
  * [31] T. Fülöp, R. Kovács, M. Szücs, and M. Fawaier.  Thermodynamical extension of a symplectic numerical scheme with half space and time shifts demonstrated on rheological waves in solids.  Entropy, 22:155, 2020.
  * [32] A. Fehér, N. Lukács, L. Somlai, T. Fodor, M. Szücs, T. Fülöp, P. Ván, and R. Kovács.  Size effects and beyond-Fourier heat conduction in room-temperature experiments.  Journal of Non-Equilibrium Thermodynamics, 46:403–411, 2021.
  * [33] K. Zhukovsky.  Violation of the maximum principle and negative solutions for pulse propagation in Guyer–Krumhansl model.  International Journal of Heat and Mass Transfer, 98:523–529, 2016.