Operational Approach and Solutions of Hyperbolic Heat Conduction Equations

We studied physical problems related to heat transport and the corresponding differential equations, which describe a wider range of physical processes. The operational method was employed to construct particular solutions for them. Inverse differential operators and operational exponent as well as operational definitions and operational rules for generalized orthogonal polynomials were used together with integral transforms and special functions. Examples of an electric charge in a constant electric field passing under a potential barrier and of heat diffusion were compared and explored in two dimensions. Non-Fourier heat propagation models were studied and compared with each other and with Fourier heat transfer. Exact analytical solutions for the hyperbolic heat equation and for its extensions were explored. The exact analytical solution for the Guyer-Krumhansl type heat equation was derived. Using the latter, the heat surge propagation and relaxation was studied for the Guyer-Krumhansl heat transport model, for the Cattaneo and for the Fourier models. The comparison between them was drawn. Space-time propagation of a power–exponential function and of a periodic signal, obeying the Fourier law, the hyperbolic heat equation and its extended Guyer-Krumhansl form were studied by the operational technique. The role of various terms in the equations was explored and their influence on the solutions demonstrated. The accordance of the solutions with maximum principle is discussed. The application of our theoretical study for heat propagation in thin films is considered. The examples of the relaxation of the initial laser flash, the wide heat spot, and the harmonic function are considered and solved analytically.


Introduction
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 Heat Equation and Its Operational Solution
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: F(x, t) = e Φ(x,t;α,β) Θ Ŝ f (x) = e Φ(x,t;α,β) f (x + αβ t 2 , t), where Φ(x, t; α, β) = 1 3 αβ 2 t 3 + β t x, Θ = e αβt 2 ∂ x , Ŝ = e αt ∂ 2 x .The heat diffusion operator was thoroughly explored by Srivastava in [78].The exponential differential operator e ∂ 2 x (3) reduces to the exponential differential operator e ∂ x upon the application of the following integral presentation [41]: where p = √ t D in our case.Thus, the above formula reads as follows: The action of the operator of translation exp(η D) for D = D + α produces a shift e η(D+α) f (x) = e ηα f (x + η), (6) The solution (2) consists of the action of the evolution operator on the initial condition F(x, 0) = f (x), which is transformed by Ŝ and Θ.Interestingly, for the Airy initial condition ζ dζ 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 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∂ 2 x x n = H n (x, a) according to the operational definition of the Hermite polynomials of two variables H n (x, y) [74] H n (x, y) = e y ∂ 2 and we obtain the solution of DE (1) F(x) ∝ H n (x + ab, a), where a = αt, b = βt as follows: It is easy to demonstrate that for the initial condition f (x) = x k e δ x the following operational rule applies: e y∂ 2 x x k e α x = e (α x+α 2 y) H k (x + 2αy, y), (11) so that we obtain Ŝ f (x) = e δ(x+δa) H k (x + 2δa, a) = f (x, t), a = αt.The consequent action of the translation operator Θ yields the shift along the x argument and results in where 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 Distinguishing the perfect square of the operator ∂ x + δ/α, we come to the solution of Equation ( 13) in the following form: where G(x, t) satisfies Equation (1) for G with the initial condition g(x) = G(x, 0) = e δ x/α f (x).Let us set α = 1 without substantial loss of generality and choose the monomial f (x) = x k initial condition for (13).Then g(x) = x k e δ x is the initial condition for the Equation ∂ t G = (∂ 2 x + βx)G.Its solution is given by ( 12) upon the substitution of F → G and we end up with the following solution of the extended heat Equation ( 13) where ∆ 2 = tγ + t 2 δβ is the additional phase.This solution reduces to (10) in the proper limiting case.

Propagation of a Heat Surge with Fourier Heat Diffusion Type Equation
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 (x, y) = δ(x, y).With the help of Equation ( 5) we immediately obtain the action of the Ŝ operator on δ-function as follows: In two dimensions Equation ( 16) becomes f (x, y, t) = Ŝx Ŝy δ(x, y) = e −(x 2 /α+y 2 /γ)/4t /4πtαγ.Then the operators Θx,y induce the space shift and yield the particular solution (8) for the initial function f (x, y) = δ(x, y), which evolves in the Gaussian as follows: One-dimensional solution for the initial flash condition f (x) = δ(x) simply reads , and unsurprisingly resembles the solution F(x, t) = e Φ(x, t) Consider the Equation ( 13) for β = δ = 0, γ = 0.This case represents the Fourier heat propagation, including heat exchange with the environment.Consider the flash initial condition f (x) = δ(x) for it.The solution is shown in Figure 1 and it evidences significant spread of the initial function already for the times t ∈ [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 γ = 0.

Operational Solution of the Hyperbolic Heat Conduction Equation
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

Operational Solution of the Hyperbolic Heat Conduction Equation
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: , where D T is the heat conductivity and the relaxation time τ in heat conduction is extremely small (τ ≈ 10 −13 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 D(x) and ε(x) are differential operators, acting over the coordinate, dependently on the specific form of the operators ε(x) and D(x), 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 ∝ e t , and does not satisfy the requirement F(x, ∞) = 0, which otherwise can be formulated as F(x, ∞) < ∞, 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 D(x), ε(x) and on the initial condition f (x).In what follows, we consider some examples.Let us first of all consider ε(x) = ε = const and D(x) = α∂ 2 x + κ, 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 → ∞ solution reads as follows: The action of the heat diffusion operator Ŝ can be accomplished with the help of the identity (4), resulting in where ζ = 4uα.Consider an example of the initial function f (x) = e γx x n , which, is useful for the description of heat pulses of a variety of shapes, custom modelled by the sum ∑ n,γ e γx x n , γ < 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 ∈ Integers, γ ∈ Reals, which can be computed with account for (9).In the particular case of given n and γ, for the example for n = −γ = 1, we obtain For n = 2, γ = −1, κ = 0 we have the following solution: The solution for the initial function f (x) = |x| 3 e −|x| is too bulky to be presented in its analytical form; it describes the damped wave propagation as shown in Figure 2.
The solution for the initial function ( ) is too bulky to be presented in its analytical form; it describes the damped wave propagation as shown in Figure 2. , 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.By choosing α = ε = 5 in (22), i.e., D T = α/ε = 1, τ = 1/ε = 0.2, the influence of the second time derivative ∂ 2 t in the Cattaneo equation is reduced, but the heat conductivity D T = α/ε = 1 remains unchanged, as compared with the case α = ε = 1.The relaxation of the solution occurs much earlier for α = ε = 5, than for α = ε = 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. .Now the integration in ( 23) can be completed for the times , where the solution exists:

Propagation of a Heat Surge with the Hyperbolic Heat Conduction Equation
We now consider the example of initial δ-function f (x) = δ(x) for Equation (22).The action of the operator Ŝ = e −4αξ∂ 2 x in our solution ( 23) can be accomplished with the help of Equation ( 16), where . Now the integration in ( 23) can be completed for the times t > x/ √ α, where the solution exists: Physical meaning of the condition υ = t 2 − x 2 /α > 0, i.e., t > x/ √ α is that the time, needed for the initial δ(x) function to reach the point x in space, is exactly t 0 = x/ √ α and before that the heat surge is just not felt in the point x at all.The value √ α 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 ε = 10, κ = 1, α = 100 for the initial function f (x) = δ(x) 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/ √ α 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 ∂ 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 κ ε Heat propagation in three spatial dimensions in the Cattaneo model is governed by the direct three-dimensional generalization of ( 22) for κ = 0: 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  Higher values of α and ε reduce the contribution of the second order time-derivative ∂ 2 t 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 ε 2 > 4κ, κ < 0. Interestingly, despite the fact that the solution (28) was obtained for ω > 0 and υ > 0, it holds true for any ω = 0 and υ = 0.For ω < 0 and υ < 0 we obtain positive values for F(x, t), while for ω < 0 and υ > 0 as well as for ω > 0 and υ < 0 we obtain complex values for the resulting function F(x, t).
Heat propagation in three spatial dimensions in the Cattaneo model is governed by the direct three-dimensional generalization of (22) for κ = 0: , where D T is the heat conductivity.The ratio C = √ D T /τ 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 = α/ε = k/c p ρ is the thermal diffusivity, k is the thermal conductivity, ρ is the mass density, c p is the specific heat capacity, τ = 1/ε is the relaxation time, often related to the speed of the second sound C in media: τ = 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 Ŝx Ŝy Ŝz f (x, y, z), where Ŝi = e −4α ξ ∂ 2 i f (i), 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ξ 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 (x, y, z) = δ(x, y, z), modelling an intense instant point-like volume heating.The application of the heat operators Ŝi = e −4αξ∂ 2 x i gives Ŝx Ŝy Ŝz δ(x, y, z) = e x 2 +y 2 +z 2 16αξ /16(π|αξ|) 3/2 yields the following solution for HHE (29) with f (x, y, z) = δ(x, y, z): The above integral converges for t > (x 2 + y 2 + z 2 )/α and we obtain the following simple analytical expression, involving modified Bessel functions K 2 : Note, that the above solution also holds for ω = 0 and υ = 0.As well as in the one-dimensional case (28) for ω < 0, υ < 0, that is for the times, inferior t 0 = (x 2 + y 2 + z 2 )/α, when the heat wave reaches point x, we obtain positive values for the solution F(x, t) for negative κ, such that 4κ < −ε 2 .For ω < 0, υ > 0 as well as for ω > 0, υ < 0 we obtain complex solution F(x, t).In the limiting case of ω → 0, υ = 0 we obtain the finite limit for the solution: Despite that the Cattaneo heat Equation ( 22) for κ = 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 = √ D T /τ, its quantitative value disagrees with the experimental data at high frequencies and low temperatures.However, the HHE is widely used in radio engineering.

Propagation of a Harmonic Signal with the Telegraph Equation
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 (x) = e inx .The action of the operator Ŝ with account for (5) yields Ŝe inx = e ν ∂ 2 x e inx = e inx−n 2 ν , ν = −4αξ.Integration in (24) in turn yields the following solution for the telegraph Equation (22) for The harmonic solution (33) is not spreading in space, but only fading in time.Its physical meaning is clarified if we attribute the constants α, ε, κ 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.
The harmonic solution (33) is not spreading in space, but only fading in time.Its physical meaning is clarified if we attribute the constants κ ε α , , 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 our notations In terms of the resistance R L , inductivity L, capacitance C and leakage resistance R C the voltage u(x, t) along the transmission line, shown in Figure 5, is described by the one-dimensional telegraph equation [86]: In our notations The solutions (28) and (32) for the initial δ-surge are valid and converge also for is realized; it fades without the space shift (see examples in Figure 6).In the case of (R C /C − R L /L) 2 < 4n 2 /LC the voltage behavior in the circuit has space shift and its time fading depends exclusively on ε (see examples in Figure 7): The essential difference between the plots in Figures 6 and 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.    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 . 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 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.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 ≥ 1.This frequency dependent heat conductivity is typical for HHE and it is absent in the Fourier law.

Operational Solution of Guyer-Krumhansl Type Heat Equation and Heat Conduction in Thin Films
Concluding the study of the hyperbolic heat Equation ( 22) we note that it can be easily modified by adding non-commuting with ∂ 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.

Operational Solution of Guyer-Krumhansl Type Heat Equation and Heat Conduction in Thin Films
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 ε = ε − δ∂ 2 x and D = α∂ 2 x + κ 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 κ = 0 for the sake of generality.Its role is not indifferent and it will be discussed in what follows.Introducing the operators ∆1 = ∂ t − (α/ε)∂ 2 x and ∆2 = ∂ t − δ∂ 2 x , we can rewrite (36) with positive coefficients α, ε, δ > 0 in the following form: Note that equations ∆1,2 F F 1 ,F 2 (x, t) = 0 represent Fourier law.If δ = α/ε and κ = 0, we have ∆ 1 = ∆ 2 , and, provided the Fourier Equation ∆1,2 F F (x, t) = 0 is fulfilled, we have the particular case of GK equation, which is Equation (37) with κ = 0, also satisfied by this solution F(x, t) = F F (x, t).The form of the heat conduction equation Equation (36) with κ = 0 arises in the study of heat propagation in a one-dimensional thin film with account for the phonon transport [108].From (21), where ε = ε − δ∂ 2 x , D = α∂ 2 x + κ, we write the solution of Equation ( 36) as follows: By means of the integral presentation e −ξδ 2 ∂ 4 (5)) we obtain after grouping the 2nd order derivative terms in the exponential x f (x)dζ the particular bounded solution for the modified hyperbolic heat conduction Equation (36) with initial function F(x, 0) = f (x) as follows: where Note, that the obtained integral solution is valid also if the constant term in the r.h.s. of the Equation ( 36) is negative: κ < 0 and ε > 2 |κ|, 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(x, t) 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: In what follows, we will consider some examples of heat pulse propagation.

Propagation of a Heat Surge with the Guyer-Krumhansl Heat Equation
Consider the example of the initial flash F(x, 0) = δ(x), corresponding to an intense and instant laser point heating at x = 0 at the moment of time t = 0. From Equation ( 16) where we denote the special function Function Φ(x; a, b) reduces for some particular values of a and b to the hypergeometric function 0 F 2 ({β 1 , β 2 }; z), this last being the particular case of p F q α 1 ...α p ; β 1 ...β q ; z .Moreover, The numerical calculation of the double integral (41), taken with care around the point a = 0, corresponding to the time , yields real values.The example of the solution (41) of the GK type Equation ( 36) with for the initial flash ( ) ( ) is shown in 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.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.

Solution of the Guyer-Krumhansl Equation for the Exponential-Polynomial Initial Function
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(x, 0) = e γ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 (x) = e γ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 (x) = e γ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 Ŝ = e ν ∂ 2 x , ν = tδ/2 − 4ξα + 2ξεδ + i2 √ ξδζ, Ŝ e γx x n = e γ x+γ 2 ν H n (x + 2γν, ν).The solution for arbitrary values of n ∈ Integers, γ ∈ Reals can be obtained upon the following integration: where The integration can be done in elementary functions if we account for the series presentation . We omit the final result for conciseness.
For example, in the simplest case of n = −γ = 1, i.e., f (x) = e −x x we obtain the following solution of the GK equation: Let us now consider the initial smooth function F(x, 0) = (|x| + 2) 2 e −|x| .By means of the above developed operational technique we obtained the following exact analytical solution of the GK equation in elementary functions: Direct substitution of the solution satisfies the GK equation.Earlier it was demonstrated that such solutions can show bizarre behavior with local maximums and local minimums even with negative values, occurring in the middle of the domain, so that theoretically for a certain set of parameters of the GK equation the maximum principle can be violated [58].However, for applications, such as those in the thin films, modelled by the GK type Equation ( 40), there is no problem with the second law of thermodynamics violation.Indeed, the propagation of the smooth spatial heat distribution F(x, 0) = (|x| + 2) 2 e −|x| in the range of the values of the Knudsen number between Kn = 0.05 ÷ 2 changes no more than three percent.
The difference between these plots is not distinguishable visually and we present here just one example for Kn = 0.2 in Figure 11.It corresponds to the values α = 2 15 , δ = 0.12, ε = 2, κ = −1.The behavior of the solution surprisingly resembles that for α ≈ ε ≈ δ ≈ 1, which we omit for brevity, as well as that for Kn = 2, for which α = 13 1  3 , δ = 12, ε = 2, κ = −1.To contrast it, we demonstrate the behavior of the solution of GK equation (see (36), κ = 0) for relatively small contribution of the Cattaneo's term: α = ε = δ = 10, which does not depend on the sign and on the exact value of the constant term κ (we omit proper figures with κ = 0) for conciseness) and it has bizarre non-Fourier behavior, shown in Figure 12.
The initial pulse rapidly decreases and assumes negative values, then it gradually approaches zero.It maintains negative values in the whole space domain, but for the vicinity of x = 0. Around x = 0 the solution becomes positive: F(x ≈ 0, t > 0.8) > 0 and then it relaxes to zero: F(x, t → ∞) = 0 .
In this case, based on the obtained exact vanishing analytical solution of the GK equation with the initial function F(x, 0) = (|x| + 2) 2 e −|x| and F(x, ∞) = 0, we conclude the minimum of the solution over the time-space domain is in the middle of the domain (see Figure 12).The local maximum of the solution may also occur in the middle of the domain (see [58,59])., we conclude the minimum of the solution over the time-space domain is in the middle of the domain (see Figure 12).The local maximum of the solution may also occur in the middle of the domain (see [58,59]).

Harmonic Solution of Guyer-Krumhansl Equation and Temperature Distribution in Thin Films
The other good example of a given initial temperature distribution, is given by a harmonic function ( ) ( ) , which is necessary to consider approximations or expansion of the initial function into the Fourier series to fit experimental data distributions.The action of the heat conduction operator yields the following bounded at infinite times solution: ( ) ( Upon integration we end up with the explicit solution of the GK Equations ( 36) and (37) for the periodic initial function , expressed in elementary functions as follows: ( ) Interestingly, the solution (47) of GK type Equation (36) for the initial harmonic function , where (33).In other words, the solution (47) of Equation GK type (36) for is also the solution of the HHE (22) with the harmonic dependent coefficient δ ε for the first order time derivative:  36) and (48), these share identical particular solution (47).For the initial function ( ) ( ) , the solution of Equations ( 36) and ( 48) will be given by the series , , where ( ) t x F , is the solution of (47).In this sense the GK type equation with

Harmonic Solution of Guyer-Krumhansl Equation and Temperature Distribution in Thin Films
The other good example of a given initial temperature distribution, is given by a harmonic function f (x) = exp(inx), which is necessary to consider approximations or expansion of the initial function into the Fourier series to fit experimental data distributions.The action of the heat conduction operator yields the following bounded at infinite times solution: Upon integration we end up with the explicit solution of the GK Equations ( 36) and (37) for the periodic initial function f (x) = exp(inx), expressed in elementary functions as follows: Interestingly, the solution (47) of GK type Equation (36) for the initial harmonic function f (x) = exp(inx) repeats the solution of telegraph Equation (33) upon the substitution ε → ε + n 2 δ , where δ > 0, ε > 0, in (33).In other words, the solution (47) of Equation GK type (36) for f (x) = e inx is also the solution of the HHE (22) with the harmonic dependent coefficient ε + n 2 δ for the first order time derivative: Thus, for the initial harmonic function f (x) = e inx and same values of the coefficients α, ε, κ, δ in Equations ( 36) and (48), these share identical particular solution (47).For the initial function φ(x) = ∑ n c n exp(inx), the solution of Equations ( 36) and (48) will be given by the series Φ = ∑ n c n F(x, t), where F(x, t) is the solution of (47).In this sense the GK type equation with the harmonic initial function can be viewed as the telegraph equation with the harmonic dependent coefficient for ∂/∂t: (ε + n 2 δ)∂/∂t and the GK equation for f (x) = e inx is in fact the Cattaneo equation for f (x) = e inx with (ε + n 2 δ)∂ t instead of ε∂ t term.This radically changes the solution behavior.Indeed, the relaxation time for all higher harmonics e inx in the telegraph Equation (33) is the same; only low harmonics may fade out faster.On the contrary, the solution (47) of the GK type Equation (36) for the harmonic function f (x) = e inx and the solution of the equivalent telegraph equation with ε → ε + n 2 δ (48) for high harmonics show that these harmonics fade out faster than the fundamental one e ix (n = 1) and the higher the harmonic number n, the lower is its relaxation time.The examples for n = 3, 5 are shown in Figure 13.Compare Figure 13 with the solution of the telegraph Equation ( 22) (see Figures 6 and 7).High harmonics of the telegraph equation, shown in Figure 7, fade out slower than the low harmonics, shown in Figure 6.For the GK type equation or the HHE with δ ε ε 2 n + → on the contrary, high harmonics fade out faster than the fundamental and the low ones, as evidenced in Figure 13.The effective thermal conductivity in the common telegraph Equation ( 22) is constant for higher harmonics and in certain cases even for all n.On the contrary, in the GK type Equation ( 36  Compare Figure 13 with the solution of the telegraph Equation ( 22) (see Figures 6 and 7).High harmonics of the telegraph equation, shown in Figure 7, fade out slower than the low harmonics, shown in Figure 6.For the GK type equation or the HHE with ε → ε + n 2 δ on the contrary, high harmonics fade out faster than the fundamental and the low ones, as evidenced in Figure 13.The effective thermal conductivity in the common telegraph Equation ( 22) is constant for higher harmonics and in certain cases even for all n.On the contrary, in the GK type Equation (36), which is HHE with ε → ε + n 2 δ for f (x) = exp(inx) (see (48)), the thermal conductivity for the harmonic function f (x) = e inx rises with increase of the harmonic number n.
In the context of the heat conduction in thin films, for Kn = 1 we plotted the bounded solutions of Equation ( 40), for f (x) = e inx , n = 1, 3, in Figures 14 and 15 respectively.High harmonics fade out rapidly as follows from their comparison with teach other.For Kn = 0.2 we present the solutions in Figures 16 and 17 for n = 1, 3 respectively.The comparison with Figure 14, Figure 15 shows that the solutions for Kn = 0.2 relax much slower than those for Kn = 1.This means that high values of the Knudsen number provide the ballistic transport contribution to the heat conduction, and the latter increases in the GK model, in particular, for high harmonics.The proper relaxation times for the first harmonic f (x) = e ix read as follows: n=1 τ Kn=1 ≈ 0.5 and n=1 τ Kn=0.2 ≈ 2.5; for the third harmonic f (x) = e 3ix the relaxation times are n=3 τ Kn=1 ≈ 0.2 and n=3 τ Kn=0.2 ≈ 3. The ratios of proper relaxation times are n=1 τ Kn=1 / n=1 τ Kn=0.2 = 1/5 and n=3 τ Kn=1 / n=3 τ Kn=0.2 = 1/15.The heat conduction is evidently better for Kn = 1 than for Kn = 0.2.

Conclusions
We obtained exact analytical solutions for heat conduction in the Fourier, Cattaneo, and Guyer-Krumhansl heat transport models.Extended forms of these equations with additional constant terms were solved by the operational method.The obtained solutions directly satisfy the considered equations.The initial functions, modelling various types of pulses and signals were considered.The power-exponential pulse and the flash heat pulse were considered, modelling the laser heat pulse experimental technique, common for thermal conductivity measurements.
The analytical solutions were compared with each other.The finite speed of the heat propagation in the Cattaneo model was demonstrated.By reducing the effect of the second time derivative in the equation for α = ε = 10, maintaining the heat conductivity unchanged, we obtain faster fading of the solution.In this case diffusive heat conduction prevails over the wave-like propagation process.On the contrary, for α = ε = 0.1 the wave-like propagation of the initial function dominates.For α = ε = 1, both heat transport mechanisms contribute; the wave propagates at a constant speed, accompanied by fading.
The exact analytical solution of HHE in three space dimensions for initial flash f = δ(x, y, z) was obtained in terms of modified Bessel functions.The HHE was solved with the harmonic initial function f = e inx .The obtained solution fades in time exponentially and does not spread.The space phase shift of the solution depends on the sign of the linear term in HHE.
The GK type heat equation with linear term was solved operationally; the solution was written in terms of integrals and special functions.The propagation of the initial heat flash was studied and demonstrated for the model of heat conduction in thin films.Small values of Knudsen number result in slower relaxation of the initial heat surge; for Kn ≈ 1 the ballistic heat transport contributes and the heat surge fades out much faster.An exact solution for the evolution of the exponential-monomial function F(x, 0) = x n e −γx in the GK equation was obtained.This allows not only the propagation of a surge with power rise, followed by the common exponential fade, to be studied but also the result for the initial function to be obtained, which can be expanded in the series f (x) = ∑ n,γ c n e γ x n or approximated by them.The example of the initial smooth heat pulse relaxation in a thin film was considered.
We found that variation of the Knudsen number values in a reasonable range Kn ∈ [0.05, 2] practically does not influence the relaxation time of the initial smooth pulse: the values of the solution of GK type equation vary less than 3% for 0.05 < Kn < 2. Thus in thin films the conditions on the surface, expressed by the Knudsen number, are particularly important for laser heat flash, but not that important for relaxation of a wide heat spot without a sharp boundary.
We obtained exact analytical solutions in terms of the orthogonal polynomials and special functions for three different models of heat conduction: Fourier, Cattaneo, and Guyer-Krumhansl.These solutions allow analytical modelling of the space-time propagation for heat flashes and pulses, which is the established experimental technique.We can also model the relaxation of any initial function, expandable in the Fourier series and in the series of the exponential-polynomial function.This series approximation is applicable practically to any pulse.
In conclusion we stress that the operational technique, used for the solution of hyperbolic heat equations, yields exact analytical solutions.They have clear physical meaning and allow easy understanding of the role of different terms in the equations.The validity of the obtained analytical solutions was checked by their substitution in the equations, the latter were satisfied.Our study not only has theoretical interest, but it also provides practical results of experimental technique and measurements.The obtained exact analytical solutions can be used as a benchmark for numerical solutions of more sophisticated forms of the studied equations.

Figure 1 .
Figure 1.Evolution of the initial δ(x) function as the solution of the Fourier Equation for

Figure 4
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

.
Interestingly, despite the fact that the solution (28) was obtained for

1 =.
is the thermal diffusivity, k is the thermal conductivity, ρ is the mass density, cp is the specific heat capacity, ε τ is the relaxation time, often related to the speed of the second sound C in media: The operational solution of (29) includes the action of the heat diffusion operators for each coordinate on the initial function:

Figure 5 .
Figure 5. Schematic diagram of an electric cable line with leakage.In terms of the resistance RL, inductivity L, capacitance C and leakage resistance RC the voltage ( ) t x u , along the transmission line, shown in Figure 5, is described by the one-dimensional telegraph equation [86]:

Figure 5 .
Figure 5. Schematic diagram of an electric cable line with leakage.

Figure 6 .
Figure 6.Space-time distribution of

Figure 7 .
Figure 7. Space-time distribution of

2 
sgnb.The numerical calculation of the double integral(41), taken with care around the point a = 0, corresponding to the time t = 4ξ 2 α δ − ε , yields real values.The example of the solution (41) of the GK type Equation (36) with α = δ = ε = κ = 1 for the initial flash f (x) = δ(x) is shown in Figure8.

1 =
of heat conduction in thin films we solved the respective GK equation for Kn .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 Figure9for Knudsen number Kn = 1 and in Figure10for 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.
Axioms 2016, 5,28 19 of 29 and on the exact value of the constant term κ (we omit proper figures with κ ≠ 0) for conciseness) and it has bizarre non-Fourier behavior, shown in Figure12.The initial pulse rapidly decreases and assumes negative values, then it gradually approaches zero.It maintains negative values in the whole space domain, but for the vicinity of case, based on the obtained exact vanishing analytical solution of the GK equation with the initial function ( ) ( )

Figure 11 .
Figure 11.Solution of GK type equation for heat transport in thin films with Knudsen number Kn = 0.2 for the initial function ( ) ( ) x e x x F −
the harmonic initial function can be viewed as the telegraph equation with the harmonic dependent coefficient for radically changes the solution behavior.Indeed, the relaxation time for all higher harmonics inx e in the telegraph Equation (33) is the same; only low harmonics may fade out faster.On the contrary, the solution(47) of the GK type Equation(36) for the harmonic function ( ) for high harmonics show that these harmonics fade out faster than the fundamental one ix e (n = 1) and the higher the harmonic number n, the lower is its relaxation time.The examples for n = 3, 5 are shown in Figure13.

Figure 13 . 3 = 5 =
Figure 13.Solution of GK type Equation: of the harmonic number n.In the context of the heat conduction in thin films, for Kn = 1 we plotted the bounded solutions of Equation (40), for ( ) inx e x f = , n = 1, 3, in Figure14and 15 respectively.High harmonics fade out rapidly as follows from their comparison with teach other.For Kn = 0.2 we present the solutions in Figures16 and 17for n = 1, 3 respectively.The comparison with Figure14, Figure15shows that the solutions for Kn = 0.2 relax much slower than those for Kn = 1.

Figure 17 .x f 3 =.
Figure 17.The behavior of the 3rd harmonic (n = 3) of the GK type Equation (36) solution for Kn = 0.2.This means that high values of the Knudsen number provide the ballistic transport contribution to the heat conduction, and the latter increases in the GK model, in particular, for high harmonics.The proper relaxation times for the first harmonic ( ) ix e x f = read as follows: n=1τKn=1 ≈ 0.5 and

Figure 17 .x f 3 =Figure 17 .
Figure 17.The behavior of the 3rd harmonic (n = 3) of the GK type Equation (36) solution for Kn = 0.2.This means that high values of the Knudsen number provide the ballistic transport contribution to the heat conduction, and the latter increases in the GK model, in particular, for high harmonics.The proper relaxation times for the first harmonic ( ) ix e x f = read as follows: n=1τKn=1 ≈ 0.5 and