Exact Negative Solutions for Guyer – Krumhansl Type Equation and the Maximum Principle Violation

Heat propagation in the Guyer–Krumhansl model is studied. The exact analytical solutions for the one-dimensional Guyer–Krumhansl equation are obtained. The operational formalism is employed. Some examples of initial functions are considered, modeling various initial heat pulses and distributions. The effect of the ballistic heat transfer in an over–diffusive regime is elucidated. The behavior of the solutions in such a regime is explored. The maximum principle and its violation for the obtained solutions are discussed in the framework of heat conduction. Examples of negative solutions for the Guyer–Krumhansl equation are demonstrated.


Introduction
Heat conduction is one of the most common physical processes encountered in everyday life.Under normal conditions in macroscopic samples of relatively homogeneous materials the Fourier law describes well heat transport, linearly relating the temperature change with the heat flux [1].At low temperatures <15 K in some crystals and in highly inhomogeneous media even at normal conditions the Fourier law does not hold.It was first noted by Onsager [2], that Fourier law could be just an approximation, which ignores the time needed for the heat flux acceleration once the temperature has changed.In this way one of the main flaws of the Fourier law was described, i.e., an immediate rise of the heat flux in the whole body at the same time.Beyond Fourier law the most significant physical phenomena is, perhaps, the Second Sound [3], first found in heat pulse propagation experiments in crystals [4][5][6][7].The description of this phenomenon was given by Cattaneo [8,9], who proposed the proper equation τ ∂ 2  t + ∂ t T = D T ∇ 2 for phonon heat waves, where D T is the heat conductivity and τ is the relaxation time.The velocity of the heat waves propagation in matter in Cattaneo's model is finite: υ T = √ D T /τ and therefore the heat flux changes at finite rate, following the temperature change at the boundary.The respective delay τ, in heat flux change after the temperature gradient was applied, is the property of the matter and is associated with phonon-phonon interactions.The Cattaneo model was not flawless (see, for example, [10][11][12][13][14]) and its predictions did not exactly agree with experiments, giving wrong values for the second sound speed √ D T /τ in Bi and Na F crystals at low temperatures.The model of Cattaneo was succeeded by the model, developed by Guyer and Krumhansl (GK) [15], which in one dimension reduces to the following equation for the temperature [16]: where κ = 0 for GK equation; the linear term κ = 0 presumes the heat sources in Equation (1).We do not dwell on the underlying physical mechanisms for the GK model, but just note that additionally to the phonon-phonon interactions in the Cattaneo model, the GK model includes one more mechanism for the heat transfer, acting at the scales L, comparable with the phonon free path l [17,18], and describing the phenomena at the boundary of the media, which was treated neither by Fourier [1] nor by Kasimir [19].The abovementioned ballistic conditions in the framework of the phonon theory mean that the system scale is of the order of the mean free path of phonons l ≈ L and they are observed in low-dimensional systems, such as ultrathin film, in nanowires, in graphene ribbons and in carbon nanotubes, in silicon threads, etc. [20][21][22][23][24].Recently it has been shown that the GK equation also applies for the description of heat conduction in highly inhomogeneous and porous matter at other than very low temperatures [25][26][27][28][29][30][31][32], despite the fact there is neither an evident phonon mechanism, nor ballistic conditions present.Modeling of the GK phenomena and heat transfer so far has been conducted mostly by numerical methods [33][34][35][36] and considered in a more general context [37].In what follows we propose pure analytical exact solutions for GK equations, which, while being particular, explain some of the bizarre features of non-Fourier heat conduction and even go under zero temperature point in some cases.

Heat Conduction Operator
With the use of advanced mathematical methods, such as the operational approach to differential exponential operators, combined with Laplace and other integral transforms as well as special functions for differential equations (DE) solution [38][39][40][41] one can obtain particular analytical solutions for a broad range of DE, exactly accounting for all the parameters, describing relevant physical systems.
Theoretical aspects of the heat operator Ŝ: were explored by Srivastava [42].This operator provides the solution of the Fourier Equation: by means of Gauss transforms with the initial condition F(x, 0) = f (x) (see [38,41]): For the initial Gaussian distribution f (x) = exp(−x 2 ) we readily obtain the result of the action of the heat operator on it as follows: and for the Dirac δ-function we have The heat operator Ŝ has applications beyond the heat transport; the analogy can be traced to the quantum-mechanical equation process of an electric charge diffusion under a potential barrier, which in the presence of electrostatic field is described as follows: This in fact means the imaginary time scale in Schrödinger equation [43].The solution of Equation ( 7) is obtained by the action of the heat operator Ŝ = e αt∂ 2 x on the initial function F(x, 0) = f (x) and consequent action of the operator Θ f (x) ≡ e η(∂ x +α) f (x) = e ηα f (x + η): where Θ = e αβt 2 ∂ x , Φ(x, t; α, β) = αβ 2 t 3 /3 + β t x.In quantum mechanics the above solution has the meaning of the probability amplitude for the charge to appear in the point x at the moment t, moving from the point x = 0 at the moment t = 0. Same regards the concentration of the particles, diffusing in the electric field under the potential barrier, provided their interaction with each other can be neglected.
For an arbitrary function f the solution is given by the convolution of this function with the kernel χ as follows: where the kernel χ represents the solution of Equation ( 9) for the initial δ-function F(x, 0) = δ(x): The above described differential operator Ŝ plays the key role in the exact solutions for heat conduction also beyond Fourier law and for other DE [44][45][46][47][48].

Exact Operational Solution for Guyer-Krumhansl Equation
The one-dimensional Guyer-Krumhansl equation for the temperature, Equation (1), can be conveniently written in the following form: where τ = 1/ε, µ = κ/ε, k T = α/ε is the Fourier heat conductivity, k b = δ/ε is the ballistic type heat conductivity; the term ballistic will be used for k b for definiteness even when phonon transport does not apply.GK equation is in fact the particular case of the following more general type of equation: where D(x) is the differential operator, acting on x coordinate.Applying Laplace transforms, we obtain the following exact solution in the form of the integral, containing several operational exponents, consequently transforming the initial function f (x): provided the integral converges.For ε = const and for the operator D(x) = α∂ 2 x + κ we have the well-known one-dimensional telegraphers equation (TE): which solution was proposed, for example, in [49].The particular, bounded and vanishing at infinite time F(x, 0) = f (x), F(x, ∞) = 0 solution for Equation ( 14) can be obtained from Equation ( 13) by substituting proper operator D(x) in it.The more general solution for Equation ( 14) reads: where D(x) = α∂ 2 x + κ and C 1,2 (x) are determined from the initial conditions.For GK equation, Equations ( 1) and (11), ε(x) = ε − δ ∂ 2 /∂x 2 and for F(x, 0) = f (x), F(x, ∞) = 0 a particular operational solution reads as follows [50]: The second branch of the solution, satisfying the initial condition F(x, 0) = f (x), can be obtained upon simultaneous substitution t → −t, ε → −ε, δ → −δ .In what follows we will study heat conduction in the framework of GK model and employ the above sketched formalism with the exact solutions.
Above we have mentioned that the initial δ-function plays fundamental role for the heat operator and allows obtaining exact solutions for equations, containing differential operator Ŝ.In the context of heat conduction the initial distribution F(x, 0) = δ(x) models the point-like laser pulse heating.This technique is well established for heat conduction experiments [51], thus exact solutions for heat equations with F(x, 0) = δ(x) have particular value.
The exact particular solution for the telegrapher's Equation ( 14) with initial Dirac δ-function was obtained by the operational method and explored in [52].For GK Equations ( 1) and (11) the solution for f (x) = δ(x) was obtained in the form of the integral in [50]: where: and C 1,2 = const, determined from the initial condition and both signs ± are possible in the integral.
For the initial monomial the action of the heat operator yields Hermite polynomials of two variables.Proper solution of the TE was explored in [44,45,52,53].These studies constituted the base for the subsequent development of the theory of relevant relativistic heat polynomials [54].For the GK type Equation (1) with the serial initial function f (x) = ∑ n x n , the solution also takes the form of the series where C 1,2 are determined from initial conditions and H n (x, y) are the Hermite polynomials of two variables [55,56] The simplest example of f (x) = x 2 yields: Entropy 2017, 19, 440 5 of 13 where the constants C 1,2 are determined from initial conditions.The exact solution for the exponential monomial f (x) = x n e γx can be obtained with the help of the following identity: For n ∈ Integers, γ ∈ Reals we obtain the following integral: where The result of the integration in Equation ( 23) with account for Equation (20) reads in terms of elementary functions, but it is cumbersome.For F(x, 0) = e −x x, we obtain: where q = 4(κ For the following function: , which can be viewed as an initial isolated spatial smooth heat pulse, the GK equation has the following solution: where the constants C 1,2 can be determined from proper boundary or initial conditions.

Violation of Maximum Principle and Negative Solutions in GK Equation
Considering GK equation, particular attention should be paid to the values of the coefficients in it to recognize the main contributing term and the relevant physical background of the heat transfer in each case.Depending on the values of the constants α, ε, δ one or the other heat transport mechanism contributes more or less.Note, also, that for α/δ = ε the GK equation solution also satisfies the Fourier heat equation.Due to the fact that GK equation, which was originally derived in the framework of phonon theory, appears to describe not only systems at low temperature and reduced dimension, where the ballistic conditions l ≈ L can be realized, but it also describes heat condition in macroscopic samples at room temperature, where we cannot talk about phonons, we have to clarify that the regime, in which the condition δ > α/ε holds should rather be called over diffusive, following [27].In what follows we will demonstrate that the behavior of some solutions for GK equation becomes unusual from thermodynamic point of view, and we will develop this topic, earlier touched on in [57].We do not intend to consider here complicated cases since the relevant calculations are cumbersome, but we will just illustrate the possibility of bizarre solutions, considering the simple case of the solution for F(x, 0) = e −|x| x 2 : One branch of the solution for α = ε = δ = 1, κ = −0.5 is shown in Figure 1a.In this case the contribution of the Fourier, Cattaneo and ballistic terms in Equation ( 11 (26) One branch of the solution for α = ε = δ = 1, κ = −0.5 is shown in Figure 1a.In this case the contribution of the Fourier, Cattaneo and ballistic terms in Equation ( 11 1), where its solution assumes negative values, is seen in Figure 2. The initial derivative is negative.
In this case we have slightly negative initial values of the function derivative.The behaviour of the solution completely differs from that of Fourier solution.Note, that the solution in Figure 2a assumes negative values.Moreover, observe the second branch of the solution in Figure 2b, which diverges at large times.Now, consider over diffusive case α = ε = 1, δ = 10, κ = −0.5, δ α ε  ; the behavior of the solution is shown Figure 3a.Observe the rising solution due to the positive initial derivative (see Figure 3b with consequent fade.The solution is positive, but the maximum principle is violated since the maximum is inside the domain.1), where its solution assumes negative values, is seen in Figure 2. The initial derivative is negative.
In this case we have slightly negative initial values of the function derivative.The behaviour of the solution completely differs from that of Fourier solution.Note, that the solution in Figure 2a assumes negative values.Moreover, observe the second branch of the solution in Figure 2b, which diverges at large times.Now, consider over diffusive case α = ε = 1, δ = 10, κ = −0.5, δ α/ε; the behavior of the solution is shown Figure 3a.Observe the rising solution due to the positive initial derivative (see Figure 3b with consequent fade.The solution is positive, but the maximum principle is violated since the maximum is inside the domain.Now consider over diffusive solution in the case, when all transport contributions are of the same order, α = ε = δ = 10, κ = 1, (see Figure 4).
As seen in Figure 4 the solution first rapidly decreases due to negative sign of the initial derivative (see Figure 4b), assumes negative values and then flips back to positive and eventually relaxes to zero, as seen in Figure 4a.This behavior also does not comply with the maximum principle and non-negativity of the solution is not preserved.
Thus, we have demonstrated that over diffusive solution of GK Equation in some cases with strong negative or positive initial derivative can reach maximum or minimum inside the domain 2 and 4).Moreover, the solution (26) for ( ) into negative region!This occurs in over-diffusive regime for various contributions of the heat conducting mechanisms of Fourier, Cattaneo and Guyer-Krumhansl.For all comparable with each other contributions it is shown in Figure 2.For the dominating ballistic term see Figure 3 and for relatively small Cattaneo wave relaxation time see Figure 4.  Now consider over diffusive solution in the case, when all transport contributions are of the same order, α = ε = δ = 10, κ = 1, (see Figure 4).
As seen in Figure 4 the solution first rapidly decreases due to negative sign of the initial derivative (see Figure 4b), assumes negative values and then flips back to positive and eventually relaxes to zero, as seen in Figure 4a.This behavior also does not comply with the maximum principle and non-negativity of the solution is not preserved.
Thus, we have demonstrated that over diffusive solution of GK Equation in some cases with strong negative or positive initial derivative can reach maximum or minimum inside the domain 2 and 4).Moreover, the solution (26) for ( ) into negative region!This occurs in over-diffusive regime for various contributions of the heat conducting mechanisms of Fourier, Cattaneo and Guyer-Krumhansl.For all comparable with each other contributions it is shown in Figure 2.For the dominating ballistic term see Figure 3 and for relatively small Cattaneo wave relaxation time see Figure 4. Now consider over diffusive solution in the case, when all transport contributions are of the same order, α = ε = δ = 10, κ = 1, (see Figure 4).
As seen in Figure 4 the solution first rapidly decreases due to negative sign of the initial derivative (see Figure 4b), assumes negative values and then flips back to positive and eventually relaxes to zero, as seen in Figure 4a.This behavior also does not comply with the maximum principle and non-negativity of the solution is not preserved.
Thus, we have demonstrated that over diffusive solution of GK Equation in some cases with strong negative or positive initial derivative can reach maximum or minimum inside the domain t ⊂ [0, ∞], x ⊂ [−L, L] (see Figures 2 and 4).Moreover, the solution (26) for F(x, t) can evolve from initially positive values f (x) > 0 into negative region!This occurs in over-diffusive regime for various contributions of the heat conducting mechanisms of Fourier, Cattaneo and Guyer-Krumhansl.For all comparable with each other contributions it is shown in Figure 2.For the dominating ballistic term see Figure 3 and for relatively small Cattaneo wave relaxation time see Figure 4.

Solution for Guyer-Krumhansl Temperature Distribution in Thin Films
Let us now consider the GK Equation, proposed for the description of the ballistic component of the temperature distribution in thin films [17]: where ( ) t x F , is the ballistic quasi-temperature and Kn is Knudsen number, accounting for the low dimensional system.Equation ( 27) is the GK type Equation (1) with  19) and ( 23) yield two branches of the solution.For Kn = 0.1 they are slightly different from each other.It is difficult to distinguish them from each other in figures, so one of the solutions is presented in Figure 5a, the difference between the solutions is shown in Figure 5b.

Solution for Guyer-Krumhansl Temperature Distribution in Thin Films
Let us now consider the GK Equation, proposed for the description of the ballistic component of the temperature distribution in thin films [17]: where F(x, t) is the ballistic quasi-temperature and Kn is Knudsen number, accounting for the low dimensional system.Equation ( 27) is the GK type Equation ( 1) with α = 10 3 Kn 2 b , ε = 2, δ = 3Kn 2 b , κ = −1.Not entering into the details of the physical problem of heat conduction in thin films, we just solve Equation ( 27) to obtain the temperature distribution T b (x, t) with initial smooth spatial heat wave-like distribution: f (x) = x 2 e −x .Equations ( 19) and ( 23) yield two branches of the solution.For Kn = 0.1 they are slightly different from each other.It is difficult to distinguish them from each other in figures, so one of the solutions is presented in Figure 5a, the difference between the solutions is shown in Figure 5b.

Solution for Guyer-Krumhansl Temperature Distribution in Thin Films
Let us now consider the GK Equation, proposed for the description of the ballistic component of the temperature distribution in thin films [17]: where ( ) t x F , is the ballistic quasi-temperature and Kn is Knudsen number, accounting for the low dimensional system.Equation ( 27) is the GK type Equation (1) with  19) and ( 23) yield two branches of the solution.For Kn = 0.1 they are slightly different from each other.It is difficult to distinguish them from each other in figures, so one of the solutions is presented in Figure 5a, the difference between the solutions is shown in Figure 5b.For Kn = 1 the two branches of the solution behave in a different from each other way with time passing.In particular, the second branch of the solution strongly senses the initial conditions.The examples are seen in Figures 6-8.Note that both negative or strongly positive solutions, violating maximum principle, may occur for Kn = 1.However, due to the fact that the solution at infinite time should eventually fade, this second branch of the solution, diverging at t → ∞ , seems unphysical.
Entropy 2017, 19, 440 10 of 14 For Kn = 1 the two branches of the solution behave in a different from each other way with time passing.In particular, the second branch of the solution strongly senses the initial conditions.The examples are seen in Figures 6-8.Note that both negative or strongly positive solutions, violating maximum principle, may occur for Kn = 1.However, due to the fact that the solution at infinite time should eventually fade, this second branch of the solution, diverging at , seems unphysical.For Kn = 1 the two branches of the solution behave in a different from each other way with time passing.In particular, the second branch of the solution strongly senses the initial conditions.The examples are seen in Figures 6-8.Note that both negative or strongly positive solutions, violating maximum principle, may occur for Kn = 1.However, due to the fact that the solution at infinite time should eventually fade, this second branch of the solution, diverging at , seems unphysical.
(a) (b) For Kn = 1 the two branches of the solution behave in a different from each other way with time passing.In particular, the second branch of the solution strongly senses the initial conditions.The examples are seen in Figures 6-8.Note that both negative or strongly positive solutions, violating maximum principle, may occur for Kn = 1.However, due to the fact that the solution at infinite time should eventually fade, this second branch of the solution, diverging at , seems unphysical.
(a) (b) The intermediate case Kn = 0.3 is shown in Figure 9.The second branch of the solution, presented in Figure 9b, exhibits waving behaviour, but it still fades out at large times.The first branch of the solution, shown in Figure 5a for Kn = 0.1, Figure 6a for Kn = 1 and in Figure 9a for Kn = 0. Eventually, consider the propagation of the initial Dirac δ-function.In the case of Kn = 0.3 the simplest particular solution (17) of GK Equation ( 1 10a.The detail of the solution for t = 0.0001 is shown in Figure 10b.
The intermediate case Kn = 0.3 is shown in Figure 9.The second branch of the solution, presented in Figure 9b, exhibits waving behaviour, but it still fades out at large times.The first branch of the solution, shown in Figure 5a for Kn = 0.1, Figure 6a for Kn = 1 and in Figure 9a for Kn = 0. We observe in Figure 10b that the solution has small negative values around the base of the evolved δ-function.It is better seen for short times; we have verified it for the range of Knudsen numbers Kn ⊂ [0.01, 10].

Conclusions
With the help of the operational approach, combined with integral transforms and involving extended forms of Hermite polynomials, we have obtained exact analytical solutions for Guyer-Krumhansl type equations.The evolution of several initial functions has been studied, such as polynomial, exponential-polynomial and Dirac δ-function.The solutions were studied for the effect of Fourier heat diffusion, Cattaneo heat waves and of the ballistic heat transport.In particular, the example of evolution of the function f (x) = ∑ n,γ c n e −γ x n , its particular case F(x, 0) = x 2 e −|x| was explored.The cases when Cattaneo waves or Fourier diffusion or ballistic effects dominate were explored.A small Guyer-Krumhansl term ~δ enhances thermal conduction and relaxation of the solution, counteracting the Cattaneo term ~ε; the latter describes the heat wave, which propagates at finite speed and preserves the pulse shape.In over diffusive regime, δ > α/ε, heat transport may differ totally from Fourier.The solution may intensively grow or decrease, depending on the initial conditions, and oscillating solutions may occur.The maximum of the solution can be achieved inside the domain and neither in the initial moment, nor at the boundaries; it was demonstrated for α = ε = 1, δ = 10.Moreover, in Guyer-Krumhansl equation, in particular, in over diffusive regime δ > α/ε the solution may temporally assume negative values, evolving from initially positive values to zero through negative region.The minimum can be reached inside the domain and the positivity is not preserved.In diffusive regime, δ > α/ε, for example, for α = ε = δ = 10, the solution around the initial maximum point suddenly drops below zero, then rises above zero and then relaxes to zero.Again the maximum principle and the non-negativity are not preserved.Expectably, if the Fourier heat diffusion term dominates in the equation, the solution also exhibits Fourier-like behavior and the ballistic conditions as well as the Cattaneo heat waves do not have significant effect of the solution.If all heat transport mechanisms contribute equally, minimum of the solution can nevertheless be inside the domain and the solution can turn negative, then fading to zero.Eventually, if δ = α/ε, the Fourier solution is realized even is δ α, δ ε.The example of the GK heat conduction in a thin film though demonstrates converging solutions for physical boundary conditions, which behave in accordance with the Second Law of thermodynamics; the maximum principle holds for non-positive initial derivative.The second branch of the solution may, however, diverge at large times, which seems unphysical.The considered examples demonstrate that in particular in over diffusive regime, dependently on the initial conditions, for strongly positive or strongly negative initial derivatives, the solution does not satisfy the maximum principle, it may become negative and may oscillate.
In forthcoming publications we will address the extension of the obtained solutions to some two-and three-dimensional cases.We will also compare and analyze the behavior of the front heat waves of the solutions to Cattaneo and to Guyer-Krumhansl Equations and include more complicated boundaries into consideration.
) are of comparable with each other order: τ = k T = k b = 1, µ = +0.1.Despite that we observe in Figure 1 the solution, developing from the initial function F(x, 0) = e −|x| x 2 , which reaches its maximum inside the domain and then fades out, thus violating maximum principle Entropy 2017, 19, 440 7 of

.
Not entering into the details of the physical problem of heat conduction in thin films, we just solve Equation (27) to obtain the temperature distribution ( ) t x T b , with initial smooth spatial heat wave-like distribution: ( )

.
Not entering into the details of the physical problem of heat conduction in thin films, we just solve Equation (27) to obtain the temperature distribution ( ) t x T b , with initial smooth spatial heat wave-like distribution:

Figure 5 .
Figure 5.Initial heat distribution f (x) = x 2 e −x evolution in one-dimensional film with Knudsen number Kn = 0.1.(a) One solution; (b) The difference between two solutions.