Operational Approach and Solutions of Hyperbolic Heat Conduction Equations

The differential equations (DE) are of paramount importance both in pure mathematics and in physics since they describe a very broad range of physical processes. Rapid development of computer methods and machine calculations in the 21st century facilitated equation solving. A good description of the major numerical methods is given, for example, in [1,2,3,4,5,6]. Here they allow numerical modelling of complicated physical processes [7,8,9,10,11,12,13,14,15,16,17,18,19], including multidimensional heat transfer in rectangles and cylinders [20,21,22,23]. However, proper understanding of the solutions and of the obtained results can be best done when they are obtained in analytical form. Analytical studies generally are more suitable for the analysis of the undergoing physical processes rather than numerical models, the latter giving precise description of the performance of the devices and of the specific studied cases. Analytical solutions are highly appreciated, but only a few types of DE allow explicit, if any, exact analytical solutions. Recently some fractional ordinary DE and partial differential equations (PDE) were analyzed and analytically solved in [24,25,26,27,28,29,30,31,32,33,34,35,36,37]. They benefit from the use of special functions [38,39,40]. The mathematical instruments, used to solve DE, generally range from a variety of integral transforms [41,42] to expansion in a series of generalized orthogonal polynomials [43] with many variables and indices [44,45,46], which arise naturally in studies of physical problems, such as the radiation and dynamics of beams of charges [47,48,49,50,51,52,53,54], heat and mass transfer [55,56,57,58,59], etc. Moreover, exponential operators and matrices are currently used also for description of such nature fundamentals as neutrino and quarks in theoretical [60,61,62,63,64,65] and in experimental [66,67,68] frameworks. The method of inverse differential and exponential operators has multiple applications for treating the above mentioned problems and related processes; some examples of DE solution by the inverse derivative method with regard to the heat equation, the diffusion equation, and their extensions, involving the Laguerre derivative, were given in [46,69,70,71,72,73]. Orthogonal polynomials can be defined in forms through operational relations [74], although we will also use their series presentations.

In what follows we treat the problem of heat conduction by the operational method in the framework of classical thermodynamics. We obtain, compare, and explore exact solutions for relevant DE in the framework of heat conduction models for Fourier [75], Cattaneo [76], and Guyer-Krumhansl [77] heat laws.

Fourier’s law of heat propagation imposes a linear relation between the temperature gradient and the heat flux. This is one of the most popular laws in continuum physics and it is in excellent agreement with everyday life and with more than 90% experiments. Recently the Fourier heat equation with a linear term was studied operationally as a special case of the Schrödinger equation in imaginary time [70,73]. The solution of such DE with the initial condition $F(x,0)=f(x)$ is given by the Gauss transform: where $\mathrm{\Phi }(x,t;\alpha ,\beta )=\frac{1}{3}\alpha {\beta }^2{t}^3+\beta tx$, $\widehat{\mathrm{\Theta }}=e^{\alpha \beta t^2{\partial }_x}$, $\widehat{S}=e^{\alpha t{\partial }_x^2}$. The heat diffusion operator was thoroughly explored by Srivastava in [78]. The exponential differential operator $e^{{\partial }_x^2}$ (3) reduces to the exponential differential operator $e^{{\partial }_x}$ upon the application of the following integral presentation [41]: where $\widehat{p}=\sqrt{t}\tilde{D}$ in our case. Thus, the above formula reads as follows:

The action of the operator of translation $\exp (\eta \tilde{D})$ for $\tilde{D}=D+\alpha$ produces a shift

The solution (2) consists of the action of the evolution operator on the initial condition $F(x,0)=f(x)$, which is transformed by $\widehat{\overline{S}}$ and $\widehat{\overline{\mathrm{\Theta }}}$. Interestingly, for the Airy initial condition $f(x)=Ai\left(\frac{x}{C}\right)=\frac{1}{\pi }{\displaystyle \underset{0}{\overset{\infty }{\int }}\cos \left(\frac{1}{3}{\zeta }^3+\frac{x}{C}\zeta \right)d\zeta }$ the operational method readily yields the fading oscillating solution without any spread [71], while the solution for the Gauss initial function demonstrates the spread of the packet.

In direct analogy with (1) the following two dimensional DE with the initial condition $F\left(x,y,0\right)=f\left(x,y\right)$ has the solution, which reads: and it includes the phase $\mathrm{\Psi }=(\alpha b^2+\gamma c^2)t^3/3+t(bx+cy)$. The solution consists of the two-dimensional Gauss transform due to the action of the heat diffusion operators (3) ${\widehat{S}}_x{\widehat{S}}_y$, and of the consequent shift along both coordinates, executed by the translation operators ${\widehat{\mathrm{\Theta }}}_x=e^{t^2\left(\alpha b/2\right){\partial }_x}$, ${\widehat{\mathrm{\Theta }}}_y=e^{t^2\left(\gamma c/2\right){\partial }_y}$.

In the context of the operational method the orthogonal polynomials are useful and their operational definitions are necessary. Indeed, the action of the heat diffusion operator on the monomial $x^n$ yields $e^{a{\partial }_x^2}{x}^n=H_n\left(x,a\right)$ according to the operational definition of the Hermite polynomials of two variables $H_n(x,y)$ [74] and we obtain the solution of DE (1) $F\left(x\right)\propto H_n\left(x+ab,a\right)$, where $a=\alpha t,b=\beta t$ as follows:

It is easy to demonstrate that for the initial condition $f\left(x\right)=x^k{e}^{\delta x}$ the following operational rule applies: so that we obtain $\widehat{S}f\left(x\right)=e^{\delta \left(x+\delta a\right)}{H}_k\left(x+2\delta a,a\right)=f\left(x,t\right)$, $a=\alpha t$. The consequent action of the translation operator $\widehat{\mathrm{\Theta }}$ yields the shift along the x argument and results in where ${\Delta }_1=\delta \left(x+\delta \alpha t+\alpha \beta t^2\right)$.

Now we can easily solve the following extended heat equation: which can be considered as the generalization of Equation (1) upon the substitution of the derivative ${\partial }_x\to {\partial }_x+\delta /\alpha$. To address Equation (13) we note that if $G\left(x,t\right)$ is the solution of Equation $\psi ({\partial }_x)G(x,t)=\widehat{D}\left(t\right)G(x,t)$ with the initial function $g(x)=G(x,0)=e^{\lambda x}F(x,0)=e^{\lambda x}f(x)$, then $F(x,t)=e^{-\lambda x}G(x,t)$ is the solution of the equation $\psi ({\partial }_x+\lambda )F(x,t)=\widehat{D}\left(t\right)F(x,t)$ with the initial function $f\left(x\right)=F(x,0)$. Distinguishing the perfect square of the operator ${\partial }_x+\delta /\alpha$, we come to the solution of Equation (13) in the following form: where $G\left(x,t\right)$ satisfies Equation (1) for $G$ with the initial condition $g\left(x\right)=G\left(x,0\right)=e^{\delta x/\alpha }f\left(x\right)$. Let us set $\alpha =1$ without substantial loss of generality and choose the monomial $f\left(x\right)=x^k$ initial condition for (13). Then $g\left(x\right)=x^k{e}^{\delta x}$ is the initial condition for the Equation ${\partial }_{t}G=({\partial }_x^2+\beta x)G$. Its solution is given by (12) upon the substitution of $F\to G$ and we end up with the following solution of the extended heat Equation (13) where ${\Delta }_2\text{}=t\gamma +t^2\delta \beta$ is the additional phase. This solution reduces to (10) in the proper limiting case.

In heat conduction experiments the relaxation of the instant point-like heat surge is the standard technique, approved for thermal diffusivity measurements [79]. In practical terms it is executed with the flash method—the standard engineering procedure for measuring thermal diffusivity with intense and ultra-short laser heat pulses. The latter can be modelled by the initial δ-function space distribution. Let us consider this example now: $f\left(x,y\right)=\delta \left(x,y\right)$. With the help of Equation (5) we immediately obtain the action of the $\widehat{S}$ operator on δ-function as follows:

In two dimensions Equation (16) becomes $f\left(x,y,t\right)={\widehat{S}}_x{\widehat{S}}_y\delta \left(x,y\right)=e^{-\left(x^2/\alpha +y^2/\gamma \right)/4t}/4\pi t\alpha \gamma$. Then the operators ${\widehat{\mathrm{\Theta }}}_{x,y}$ induce the space shift and yield the particular solution (8) for the initial function $f\left(x,y\right)=\delta \left(x,y\right)$, which evolves in the Gaussian as follows:

One-dimensional solution for the initial flash condition $f\left(x\right)=\delta \left(x\right)$ simply reads $F(x,t)=e^{\mathrm{\Psi }}\frac{e^{-\frac{{\left(x+ab{t}^2\right)}^2}{4t\alpha }}}{2\sqrt{\pi t\alpha }}$, and unsurprisingly resembles the solution $F(x,t)=\frac{e^{\mathrm{\Phi }(x,t)}}{\sqrt{1+4t}}{e}^{-\frac{{(x+\beta t^2)}^2}{1+4t}}$ [71] for the initial Gaussian function $f(x)=e^{-x^2}$.

Consider the Equation (13) for $\beta =\delta =0$, $\gamma \ne 0$. This case represents the Fourier heat propagation, including heat exchange with the environment. Consider the flash initial condition $f\left(x\right)=\delta \left(x\right)$ for it. The solution is shown in Figure 1 and it evidences significant spread of the initial function already for the times $t\in [{10}^{-3},{10}^{-2}]$.

The plot in Figure 1 is compiled for α = 1, γ = 1, which represents the Fourier heat solution with some heat exchange with the environment. At the moment t = 0.001 we see the Gaussian as the result of the evolution of the δ-function; at the moment t = 0.01 this Gaussian has faded and the contribution of the heat exchange term γ is noticeable. The relaxed solution slowly grows due to non-zero heat exchange with $\gamma \ne 0$.

The Fourier heat law [75] has some shortcomings, noted by L. Onsager in 1931, who said [80] that it “contradicts the principle of microscopic reversibility, but this contradiction is removed if we recognize that it is only an approximate description of the process of conduction, neglecting the time needed for acceleration of the heat flow”. The Fourier law does not properly describe heat conduction at low temperature <25 K in dielectric crystals and in systems with reduced dimensions. Moreover, it has some unphysical properties, such as lack of inertia. In other words, if an instant temperature perturbation is applied at a point in the solid, it will be felt everywhere instantaneously. This contradicts the phenomenon of the second sound, when the temperature perturbation propagates like a wave with damping. To overcome these problems Cattaneo [76] proposed a time-dependent relaxational model, which yielded the following equation: $(\tau {\partial }_t^2+{\partial }_t)T=D_T{\nabla }^{2}T$, where D_T is the heat conductivity and the relaxation time τ in heat conduction is extremely small (τ ≈ 10⁻¹³ s) at room temperature. This equation models the phenomenon of the second sound, first observed in liquid helium [81]. With the development of the second sound theoretical background [77], it was then later detected also in solid crystals [82,83,84,85]. In the relevant tests heat flash technology [79] was used for the sensitive measurement of the thermal diffusivity.

The application of the operational method for solution of the second order PDE was given in [72]. The following PDE with initial conditions: where $\widehat{D}\left(x\right)$ and $\widehat{\epsilon }\left(x\right)$ are differential operators, acting over the coordinate, dependently on the specific form of the operators $\widehat{\epsilon }\left(x\right)$ and $\widehat{D}\left(x\right)$, describes a variety of physical processes in thermodynamics, electrodynamics etc. The particular fading with time solution of Equation (18) formally reads as follows [72]:

The other solution of Equation (18) has a positive sign in the exponential: $F_2\propto e^t$, and does not satisfy the requirement $F\left(x,\infty \right)=0$, which otherwise can be formulated as $F\left(x,\infty \right)<\infty$, i.e., the solution converges at infinite time. Both of these assumptions are reasonable for physical applications. We perform the Laplace transforms [41] in (19): and, provided the integral converges, we obtain for Equation (18) the following vanishing at infinite time solution: whose particular form depends on the differential operators $\widehat{D}\left(x\right)$, $\widehat{\epsilon }\left(x\right)$ and on the initial condition $f\left(x\right)$. In what follows, we consider some examples.

Let us first of all consider $\widehat{\epsilon }\left(x\right)=\epsilon =\mathrm{const}$ and $\widehat{D}\left(x\right)=\alpha {\partial }_x^2+\kappa$, keeping the linear term κ in the r.h.s. of Equation (18). Then we end with the telegraph equation, also known as a hyperbolic heat conduction equation (HHE):

Its fading at $t\to \infty$ solution reads as follows:

The action of the heat diffusion operator $\hat{S}$ can be accomplished with the help of the identity (4), resulting in where $\zeta =4u\alpha$. Consider an example of the initial function $f\left(x\right)=e^{\gamma x}{x}^n$, which, is useful for the description of heat pulses of a variety of shapes, custom modelled by the sum ${\displaystyle {\sum }_{n,\gamma }{e}^{\gamma x}{x}^n},\gamma <0$. With the help of the operational identity (11) we obtain the exact form of the solution:

We omit the bulky result of the above integration for arbitrary $n\in \mathrm{Integers},\gamma \in \mathrm{Reals}$, which can be computed with account for (9). In the particular case of given n and γ, for the example for $n=-\gamma =1$, we obtain

For $n=2$, $\gamma =-1$, $\kappa =0$ we have the following solution:

The solution for the initial function $f\left(x\right)={\left|x\right|}^3{e}^{-\left|x\right|}$ is too bulky to be presented in its analytical form; it describes the damped wave propagation as shown in Figure 2.

By choosing $\alpha =\epsilon =5$ in (22), i.e., $D_T=\alpha /\epsilon =1$, $\tau =1/\epsilon =0.2$, the influence of the second time derivative ${\partial }_t^2$ in the Cattaneo equation is reduced, but the heat conductivity $D_T=\alpha /\epsilon =1$ remains unchanged, as compared with the case $\alpha =\epsilon =1$. The relaxation of the solution occurs much earlier for $\alpha =\epsilon =5$, than for $\alpha =\epsilon =1$, as follows from the comparison of the plot in Figure 2 with that in Figure 3; the reason is that the diffusive heat transfer in this case prevails over the wave heat transfer.

We now consider the example of initial δ-function $f\left(x\right)=\delta \left(x\right)$ for Equation (22). The action of the operator $\widehat{S}=e^{-4\alpha \xi {\partial }_x^2}$ in our solution (23) can be accomplished with the help of Equation (16), where $\nu =-4\alpha \xi$, as follows: $\widehat{S}\delta \left(x\right)=\frac{e^{\raisebox{1ex}{$x^2$}\!\left/ \!\raisebox{-1ex}{$16\alpha \xi $}\right.}}{4\sqrt{\pi \left|\alpha \xi \right|}}$. Now the integration in (23) can be completed for the times $t>x/\sqrt{\alpha }$, where the solution exists:

Physical meaning of the condition $\upsilon =t^2-x^2/\alpha >0$, i.e., $t>x/\sqrt{\alpha }$ is that the time, needed for the initial δ(x) function to reach the point x in space, is exactly $t_0=x/\sqrt{\alpha }$ and before that the heat surge is just not felt in the point x at all. The value $\sqrt{\alpha }$ is the velocity of the surge propagation, finite in the Cattaneo model, contrary to the infinite signal propagation speed in the Fourier law. The solution of Equation (22) for $\epsilon =10,\kappa =1,\alpha =100$ for the initial function $f\left(x\right)=\delta \left(x\right)$ is shown in Figure 4 for the moments of time t = 0.01 (light blue line), t = 0.1 (lilac line), t = 0.3 (pink line), t = 0.6 (blue line) and t = 1.1 (yellow–green line). Differently from the Fourier law, the initial flash reaches the point x at the moment $t_0=x/\sqrt{\alpha }$ and there is still no non-trivial solution in the point x before this moment of time. The initial δ-function is not shown in Figure 4 because of its infinite amplitude.

Higher values of α and ε reduce the contribution of the second order time-derivative ${\partial }_t^2$ in the Cattaneo Equation respectively to the role of other terms in (22). The damping of the solution in this case occurs sooner. We omit the proper plot for conciseness. Note that in the case of negative κ, i.e., when positive κ is in the l.h.s. of the hyperbolic Equation (22), the integral in (28) converges for ${\epsilon }^2>4\kappa$, $\kappa <0$. Interestingly, despite the fact that the solution (28) was obtained for $\omega >0\mathrm{and}\upsilon >0$, it holds true for any $\omega \ne 0\mathrm{and}\upsilon \ne 0$. For $\omega <0\mathrm{and}\upsilon <0$ we obtain positive values for $F(x,t)$, while for $\omega <0\mathrm{and}\upsilon >0$ as well as for $\omega >0\mathrm{and}\upsilon <0$ we obtain complex values for the resulting function $F(x,t)$. For $\omega \to 0,\upsilon \ne 0$ we obtain ${F\left(x,t\right)|}_{\omega \to 0}\text{}=e^{-\frac{\epsilon t}{2}}\frac{t}{2\pi {\alpha }^{1/2}}\frac{2}{\left|\upsilon \right|}$. For $\omega \ne 0,\upsilon \to 0$ the solution (28) diverges.

Heat propagation in three spatial dimensions in the Cattaneo model is governed by the direct three-dimensional generalization of (22) for κ = 0: ${\nabla }^{2}T=\frac{1}{D_T}\frac{\partial T}{\partial t}+\frac{\tau }{D_T}\frac{{\partial }^{2}T}{\partial t^2}$, where $D_T$ is the heat conductivity. The ratio $C=\sqrt{D_T/\tau }$ has the dimension of velocity and it stands for the speed of the heat wave propagation in the medium, τ is the material parameter, describing the time, needed for the initiation of a heat flow after a temperature gradient was imposed at the boundary of the domain. Thus, it determines the time lag for the appearance or disappearance of the heat flow after the temperature gradient is imposed or removed. This relaxation time is associated in the framework of the Cattaneo model with the linkage time of phonon-phonon collision and it measures the thermal inertia of the medium through the time, needed for the heat flow to fade in or out. The role of the constant term in the HHE becomes clearer in the context of the signal propagation in long telephone lines. We will consider it in what follows. At the moment we will keep for generality the constant term in the three-dimensional telegraph equation, thus writing it as follows: where $D_T=\alpha /\epsilon =k/c_p\rho$ is the thermal diffusivity, k is the thermal conductivity, ρ is the mass density, c_p is the specific heat capacity, $\tau =1/\epsilon$ is the relaxation time, often related to the speed of the second sound C in media: $\tau =D_T/C^2$. The operational solution of (29) includes the action of the heat diffusion operators for each coordinate on the initial function:

The action of ${\widehat{S}}_x{\widehat{S}}_y{\widehat{S}}_{z}f\left(x,y,z\right)$, where ${\widehat{S}}_i=e^{-4\alpha \xi {\partial }_i^2}f\left(i\right),i=x,y,z$ is easy to obtain: the result for each coordinate heat operator action is given by the inner integral in (24) and then the integration over $d\xi$ can be performed if the integral converges.

For experimental measurement of the thermal conductivity the flash technique is commonly approved. For this reason we choose the initial function $f\left(x,y,z\right)=\delta \left(x,y,z\right)$, modelling an intense instant point-like volume heating. The application of the heat operators ${\widehat{S}}_i=e^{-4\alpha \xi {\partial }_{x_i}^2}$ gives ${\widehat{S}}_x{\widehat{S}}_y{\widehat{S}}_z\delta \left(x,y,z\right)=e^{\frac{x^2+y^2+z^2}{16\alpha \xi }}/16{(\pi \left|\alpha \xi \right|)}^{3/2}$ yields the following solution for HHE (29) with $f\left(x,y,z\right)=\delta \left(x,y,z\right)$:

The above integral converges for $t>\sqrt{\left(x^2+y^2+z^2\right)/\alpha }$ and we obtain the following simple analytical expression, involving modified Bessel functions K₂:

Note, that the above solution also holds for $\omega \ne 0\mathrm{and}\upsilon \ne 0$. As well as in the one-dimensional case (28) for $\omega <0,\upsilon <0$, that is for the times, inferior $t_0=\sqrt{\left(x^2+y^2+z^2\right)/\alpha }$, when the heat wave reaches point x, we obtain positive values for the solution $F(x,t)$ for negative κ, such that $4\kappa <-{\epsilon }^2$. For $\omega <0,\upsilon >0$ as well as for $\omega >0,\upsilon <0$ we obtain complex solution $F(x,t)$. In the limiting case of $\omega \to 0,\upsilon \ne 0$ we obtain the finite limit for the solution: ${F\left(x,y,z,t\right)|}_{\omega \to 0}\text{}=e^{-\frac{\epsilon t}{2}}\frac{4t}{{\pi }^2{\alpha }^{3/2}}\frac{1}{{\upsilon }^2}$, but for $\omega \ne 0,\upsilon \to 0$ the solution diverges.

Despite that the Cattaneo heat Equation (22) for $\kappa =0$ is reversible on the time scale of the thermal relaxation time τ, and despite that it predicts a finite value for the heat wave velocity $C=\sqrt{D_T/\tau }$, its quantitative value disagrees with the experimental data at high frequencies and low temperatures. However, the HHE is widely used in radio engineering.

In the context of signal propagation in cable lines without radiation loss the HHE (22) is perfectly usable for the description of the voltage and current 86. For this reason we choose the initial harmonic function $f\left(x\right)=e^{inx}$. The action of the operator $\widehat{S}$ with account for (5) yields $\widehat{S}{e}^{inx}=e^{\nu {\partial }_x^2}{e}^{inx}=e^{inx-n^2\nu }$, $\nu =-4\alpha \xi$. Integration in (24) in turn yields the following solution for the telegraph Equation (22) for $F\left(x,0\right)=f\left(x\right)=e^{inx}$:

The harmonic solution (33) is not spreading in space, but only fading in time. Its physical meaning is clarified if we attribute the constants $\alpha ,\epsilon ,\kappa$ in HHE (22) the values according to the schematic diagram of the electric circuit in Figure 5, where the electric line with finite resistance, inductivity, capacitance, and leakage is shown.

In terms of the resistance R_L, inductivity L, capacitance C and leakage resistance R_C the voltage $u\left(x,t\right)$ along the transmission line, shown in Figure 5, is described by the one-dimensional telegraph equation [86]:

In our notations $\alpha =c^2=1/LC,\epsilon =\left(a+b\right)=R_C/C+R_L/L,\kappa =-ab=-R_C{R}_L/LC<0$, $V={\left(a-b\right)}^2-{\left(2cn\right)}^2$ and ${\epsilon }^2+4\kappa =a^2+b^2-2ab={\left(a-b\right)}^2={\left(R_C/C-R_L/L\right)}^2\ge 0$. The solutions (28) and (32) for the initial δ-surge are valid and converge also for $R_c/C=R_L/L$ provided $t\ne t_0=\sqrt{\left(x^2+y^2+z^2\right)/\alpha }$.

For ${\left(R_C/C-R_L/L\right)}^2\text{}\ge 4n^2/LC$ the periodic solution $u\left(x,t\right)=e^{inx-\frac{t}{2}\left(\epsilon +\sqrt{V}\right)},V\ge 0$ (33) is realized; it fades without the space shift (see examples in Figure 6). In the case of ${\left(R_C/C-R_L/L\right)}^2\text{}<4n^2/LC$ the voltage behavior in the circuit has space shift and its time fading depends exclusively on $\epsilon$ (see examples in Figure 7):

The essential difference between the plots in Figure 6 and Figure 7 is in the time dependent spatial phase shift, which is seen in Figure 7 for $n=4,5$ and is absent in Figure 6 for $n=2,3$.

Thus the first three harmonics for R_C = 7, C = L = R_L = 1 develop in time without the space shift with decent fading (see examples for n = 2, 3 in Figure 6); for them the solution (33) is realized. The harmonics with $n\ge 4$ are governed by the solution (35), they have a time-dependent space shift (see examples for n = 4, 5 in Figure 7). With increase of n the relaxation time of the harmonics also increases, but from a certain harmonic, for which the space shift appears, the relaxation time stabilizes and remains equal for all higher harmonics. In the context of heat conduction it means that the heat conduction for higher harmonics is lower than the heat conduction for low harmonics. It may also occur that, dependent on the values of the parameters in the HHE all harmonics have time-dependent space shift, i.e., the solution (35) is realized for $n\ge 1$. This frequency dependent heat conductivity is typical for HHE and it is absent in the Fourier law.

Concluding the study of the hyperbolic heat Equation (22) we note that it can be easily modified by adding non-commuting with ${\partial }_x$ terms and mixed derivatives over time and coordinate. Such a modified hyperbolic equation can be solved with the help of the above developed operational technique.

The Cattaneo hyperbolic equation gives a qualitatively correct description of the second sound through the heat wave propagation at finite velocity, but numerical results do not match the experiment. The predicted values for the speed of the heat wave in matter disagree with the data on ultrasonic wave propagation in dilute gases. Neither the heat pulse propagation at very low temperature is described correctly in non-metallic very pure crystals of Bi and Na F. In response to this disagreement, further generalizations of the Fourier law were developed. Among them the Guyer-Krumhansl (GK) equation [87] is distinguished for simplicity and coherence with observed data. In addition to the heat waves the GK model takes into consideration the so-called ballistic transport, which can be observed when the mean free path of a particle significantly exceeds the dimension of the medium in which it travels. The mean free path increases at low temperatures For example, for electrons in a medium with negligible electrical resistance, their motion is altered by collisions with the walls. Ballistic conduction applies also to phonons. It is typically observed in low dimensional structures, such as very thin films, silicon nanowires, carbon nanotubes, graphene etc. (see, for example, [88,89,90,91,92]). When the bulk phonon mean free path l is comparable to the structure size L, neither the Casimir phonon theory [93] nor Fourier diffusion [75] describe the heat transfer well, which is affected by boundary as well as internal scattering. The common rule that determines which type of transport dominates is the following: when l << L, the heat diffusion prevails, but for l >> L, or when the temperature gradient becomes large, the ballistic transport cannot be ignored. Recently the conditions for the transition region between the ballistic heat transport and the diffusion were described in [92]. Whether the heat transfer is ballistic or diffusive is also important for experimentalists, because the account for the ballistic phonon transport requires the knowledge of the boundary quality, while the wave and diffusive transport constants are intrinsic material properties [94,95,96]. The three dimensional Guyer-Krumhansl (GK) heat law was derived from the solution of the linearized Boltzmann equation for a phonon field in dielectric crystals at low temperature. It is good for the description of phonon gases in low temperature samples and even for some heat propagation processes on the microscopic [97] and macroscopic [98,99,100,101,102] level at room temperature [103,104,105]. The relaxation of the laser flash was shown to follow the GK law [106], which was reconsidered in the framework of weakly non-local thermodynamics in [107]. Due to the wide interest in the GK equation and its relatively good agreement with an array of experiments over a broad temperature range, conducted in various materials, we will consider the GK equation solution for propagation of spatial heat waves and flash pulses. To this end we will use the operational technique.

Let us choose the operators $\widehat{\epsilon }=\epsilon -\delta {\partial }_x^2$ and $\widehat{D}=\alpha {\partial }_x^2+\kappa$ in (18), which results in the following GK type equation with mixed derivative:

The one-dimensional GK equation is essentially Equation (36) with κ = 0; we will keep $\kappa \ne 0$ for the sake of generality. Its role is not indifferent and it will be discussed in what follows. Introducing the operators ${\widehat{\Delta }}_1={\partial }_t-(\alpha /\epsilon ){\partial }_x^2$ and ${\widehat{\Delta }}_2={\partial }_t-\delta {\partial }_x^2$, we can rewrite (36) with positive coefficients $\alpha ,\epsilon ,\delta >0$ in the following form:

Note that equations ${\widehat{\Delta }}_{1,2}{F}_{F_1,F_2}\left(x,t\right)=0$ represent Fourier law. If $\delta =\alpha /\epsilon$ and $\kappa =0$, we have ${\Delta }_1={\Delta }_2$, and, provided the Fourier Equation ${\widehat{\Delta }}_{1,2}{F}_F\left(x,t\right)=0$ is fulfilled, we have the particular case of GK equation, which is Equation (37) with $\kappa =0$, also satisfied by this solution $F\left(x,t\right)=F_F\left(x,t\right)$. The form of the heat conduction equation Equation (36) with $\kappa \ne 0$ arises in the study of heat propagation in a one-dimensional thin film with account for the phonon transport [108]. From (21), where $\widehat{\epsilon }=\epsilon -\delta {\partial }_x^2$, $\widehat{D}=\alpha {\partial }_x^2+\kappa$, we write the solution of Equation (36) as follows:

By means of the integral presentation $e^{-\xi {\delta }^2{\partial }_x^4}f(x)={\displaystyle {\int }_{-\infty }^{\infty }{e}^{-{\zeta }^2\text{}+2i\zeta \delta \sqrt{\xi }{\partial }_x^2}f(x)d\zeta }/\sqrt{\pi }$ (see (5)) we obtain after grouping the 2nd order derivative terms in the exponential ${\displaystyle {\int }_{-\infty }^{\infty }{e}^{-{\zeta }^2\text{}+\left(t\delta /2-4\xi \alpha +2\xi \epsilon \delta +2i\zeta \sqrt{\xi }\delta \right){\partial }_x^2}f(x)d\zeta }$ the particular bounded solution for the modified hyperbolic heat conduction Equation (36) with initial function $F\left(x,0\right)=f\left(x\right)$ as follows: where $\widehat{S}=e^{\nu {\partial }_x^2}$, $\nu =a+ib\zeta$, $a=t\delta /2-4\xi \alpha +2\xi \epsilon \delta ,b=2\sqrt{\xi }\delta$. Note, that the obtained integral solution is valid also if the constant term in the r.h.s. of the Equation (36) is negative: $\kappa <0$ and $\epsilon >2\sqrt{\left|\kappa \right|}$, which insures the integral convergence. Equations with negative sign of κ have sense and they arise for ballistic heat propagation in thin films [108].

The application of the above solution is direct: it was demonstrated that heat propagation in thin films obeys the extended form of GK Equation (36). In particular, the following equation for thin films was proposed [108]: where $F\left(x,t\right)$ is the ballistic component of the dimensionless energy (or quasi-temperature) and Kn is the Knudsen number, describing the molecular and the boundary effects. Equation (40) is a GK type Equation (36) with the following coefficients: $\alpha =\frac{10}{3}K{n}_b^2$, $\epsilon =2$, $\delta =3K{n}_b^2$, $\kappa =-1$.

In what follows, we will consider some examples of heat pulse propagation.

Consider the example of the initial flash $F\left(x,0\right)=\delta \left(x\right)$, corresponding to an intense and instant laser point heating at $x=0$ at the moment of time $t=0$. From Equation (16) for $\nu =a+ib\zeta ,a=t\delta /2-4\xi \alpha +2\xi \epsilon \delta ,b=2\sqrt{\xi }\delta$ we obtain $\widehat{S}\delta \left(x\right)=\frac{e^{-\raisebox{1ex}{$x^2$}\!\left/ \!\raisebox{-1ex}{$4\left(a+ib\zeta \right)$}\right.}}{2\sqrt{\pi \left|a+ib\zeta \right|}}$. Then from (39) the analytical solution of GK type Equation (36) with $f\left(x\right)=\delta \left(x\right)$ follows: where we denote the special function

Function $\mathrm{\Phi }\left(x;a,b\right)$ reduces for some particular values of a and b to the hypergeometric function ${}_0{F}_2\left(\left\{{\beta }_1,{\beta }_2\right\};z\right)$, this last being the particular case of ${}_p{F}_q\left(\left\{{\alpha }_1\mathrm{...}{\alpha }_p\right\};\left\{{\beta }_1\mathrm{...}{\beta }_q\right\};z\right)$. Moreover, $\mathrm{\Phi }(0;0,b)=\frac{1}{\sqrt{\left|b\right|}}\mathrm{\Gamma }\left(\frac{1}{4}\right)$ and ${\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\frac{e^{-{\zeta }^2-\frac{x^2}{4b\zeta }}d\zeta }{\sqrt{\left|\zeta \right|}}}=\left(\begin{array}{l}\mathrm{\Gamma }\left(\frac{1}{4}\right){\text{}}_0{F}_2\left(;\left\{\frac{1}{2},\frac{3}{4}\right\};-\frac{x^4}{64b^2}\right)-\\ \sqrt{\frac{\pi }{b}}\left|x\right|{}_0{F}_2\left(;\left\{\frac{3}{4},\frac{5}{4}\right\};-\frac{x^4}{64b^2}\right)\end{array}\right)\mathrm{sgn}b$. The numerical calculation of the double integral (41), taken with care around the point a = 0, corresponding to the time $t=4\xi \left(2\frac{\alpha }{\delta }-\epsilon \right)$, yields real values. The example of the solution (41) of the GK type Equation (36) with $\alpha =\delta =\epsilon =\kappa =1$ for the initial flash $f\left(x\right)=\delta \left(x\right)$ is shown in Figure 8.

From the comparison of Figure 8 with Figure 1 we conclude that the solution of the Fourier heat Equation (1) (see Figure 1) is more pronounced at short times, it has a higher peak, while the solution of GK Equation (36) spreads faster (see Figure 8) when the additional terms ε = δ = 1 are comparable with the others order.

In the context of heat conduction in thin films we solved the respective GK equation for $Kn=0.2$ and for $Kn=1$. The results for the ballistic dimensionless energy change after a laser pulse heating of a thin film, modelled by the GK type Equation (40), are shown in Figure 9 for Knudsen number Kn = 1 and in Figure 10 for Knudsen number Kn = 0.2. Their comparison confirms much faster heat propagation in a thin film with Kn = 1, where the ballistic transport mechanism contributes noticeably, as compared with the film with small value of Kn = 0.2.

Interesting features of the heat conduction, governed by the Guyer-Krumhansl equation, arise from consideration of the evolution of the power–exponential initial function $F\left(x,0\right)=e^{\gamma x}{x}^n$. Consideration of such a function is useful not only because it allows better approximation of experimental data than the usual polynomial $x^n$ and $f\left(x\right)=e^{\gamma x}{x}^n$, γx < 0, but also as itself it represents a pulse, and practically any solitary space wave or surge is easy to approximate by its sum $f\left(x\right)=e^{\gamma x}{x}^n$. The operational method allows the exact analytical solutions for these functions to be obtained. Indeed, we make use of the operational identity (11) to obtain for $\widehat{S}=e^{\nu {\partial }_x^2}$, $\nu =t\delta /2-4\xi \alpha +2\xi \epsilon \delta +i2\sqrt{\xi }\delta \zeta$, $\widehat{S}{e}^{\gamma x}{x}^n=e^{\gamma x+{\gamma }^2\nu }{H}_n\left(x+2\gamma \nu ,\nu \right)$. The solution for arbitrary values of $n\in \mathrm{Integers},\gamma \in \mathrm{Reals}$ can be obtained upon the following integration: where $\nu =a+ib\zeta$, $a=t\delta /2-4\xi \alpha +2\xi \epsilon \delta$, $b=2\sqrt{\xi }\delta$, $d={\epsilon }^2+4\kappa +2{\gamma }^2(2\alpha -\epsilon \delta )$. The integration can be done in elementary functions if we account for the series presentation $H_n(x+2\gamma \nu ,\nu )=n!{\displaystyle \sum _{r=0}^{[n/2]}\frac{{(x+2\gamma \nu )}^{n-2r}{\nu }^r}{\left(n-2r\right)!r!}}$ $={(-i)}^n{(a+ib\zeta )}^{\frac{n}{2}}\text{}{H}_n\left(\frac{i(x+2\gamma (a+ib\zeta ))}{2\sqrt{a+ib\zeta }}\right)$. We omit the final result for conciseness. For example, in the simplest case of $n=-\gamma =1$, i.e., $f\left(x\right)=e^{-x}x$ we obtain the following solution of the GK equation: