Transition from Diffusion to Wave Propagation in Fractional Jeffreys-Type Heat Conduction Equation

The heat conduction equation with a fractional Jeffreys-type constitutive law is studied. Depending on the value of a characteristic parameter, two fundamentally different types of behavior are established: diffusion regime and propagation regime. In the first case, the considered equation is a generalized diffusion equation, while in the second it is a generalized wave equation. The corresponding memory kernels are expressed in both cases in terms of Mittag–Leffler functions. Explicit representations for the one-dimensional fundamental solution and the mean squared displacement are provided and analyzed analytically and numerically. The one-dimensional fundamental solution is shown to be a spatial probability density function evolving in time, which is unimodal in the diffusion regime and bimodal in the propagation regime. The multi-dimensional fundamental solutions are probability densities only in the diffusion case, while in the propagation case they can have negative values. In addition, two different types of subordination principles are formulated for the two regimes. The Bernstein functions technique is extensively employed in the theoretical proofs.


Introduction and Model Formulation
The need to go beyond the Fourier heat conduction equation is experimentally proved under different conditions since decades. There are various situations where the classical Fourier constitutive law for heat conduction is not applicable (heat flux is not directly proportional to temperature gradient), such as heat conduction in heterogeneous materials or at low temperatures [1][2][3].
Several generalizations of the Fourier heat conduction equation using integer-order derivatives are reviewed in [4]. Among them is the Jeffreys-type heat conduction equation where q is the heat flux vector, τ q and τ T are generalized relaxation times, k is the thermal conductivity, D α t denotes the Riemann-Liouville fractional time-derivative of order α ∈ (0, 1), and ∇ denotes the gradient operator, acting with respect to the space variables. Combining Equation (2) with the energy balance equation where C denotes the heat capacity, one finds the following fractional Jeffreys-type heat conduction equation where D = k/C is the thermal diffusivity. For notational simplicity we set in the present work D = 1. We assume arbitrary initial temperature distribution, and zero initial heat flux q(x, 0) = 0.
The one-dimensional Cauchy problem for Equation (4) is solved analytically in ( [23] Chapter 7) by means of the Laplace and Fourier transforms in the time and space domains, respectively, and some numerical examples are given. In [24] Equation (4) is studied in the context of unidirectional flows of viscoelastic fluids. A fractional Jeffreys-type equation of a more general form is proposed and studied in [6] to address the shortcomings in the DPL model reported in [12] when τ q > τ T . Equation (4) can be represented as a generalized diffusion-wave equation with memory kernel. The generalized time-fractional diffusion equation with memory kernel is studied, e.g., in [25][26][27][28]. Generalized time-fractional wave equations with memory kernel are reviewed in [29,30]. Properties of multi-dimensional fundamental solutions with the emphasis on their positivity, are studied in the context of time-fractional and space-time fractional diffusion-wave equations in [31][32][33]. The Bernsten functions technique plays a crucial role in the definition of the classes of generalized diffusion-wave equations, as well as in the proof of positivity of their fundamental solutions.
Motivated by the aforementioned developments, in the present work we revisit the fractional Jeffreys-type heat Equation (4) with the main emphasis on the differences in behavior in the two cases: τ q < τ T and τ q > τ T . Two fundamentally different types of behavior are established: diffusion regime for τ q < τ T and propagation regime for τ q > τ T . Based on the theory of Bernstein functions, we recast Equation (4) into a generalized diffusion equation in the first case and a generalized wave equation in the second. Explicit representation for the one-dimensional fundamental solution is provided and used for numerical computation and visualization. The one-dimensional fundamental solution is shown to be a spatial probability density function evolving in time, which is unimodal in the diffusion regime and bimodal in the propagation regime. The multi-dimensional fundamental solutions are probability densities only in the diffusion case, while in the propagation case they can have negative values. In addition, two different types of subordination principles are formulated for the two regimes.
This paper is organized as follows. Section 2 contains preliminaries: definitions and basic properties of fractional derivatives, Mittag-Leffler functions, Bernstein functions and related classes of functions. In Section 3, Equation (4) is rewritten as a generalized diffusion-wave equation. The fundamental solution of the one-dimensional Cauchy problem and the mean squared displacement are obtained in Section 4. Section 5 contains numerical examples and plots. Subordination principles are formulated in Section 6 and multi-dimensional fundamental solutions are discussed briefly. Section 7 contains concluding remarks.

Preliminaries
The Riemann-Liouville fractional derivative D α t of order α ∈ (0, 1) is defined as [34,35] The Caputo fractional derivative D α t is defined by D α t u(t) = D α t (u(t) − u(0)), α ∈ (0, 1). The Laplace transform of a function u(t), t ∈ R + , is denoted by In this work we consider only functions, which Laplace transform exists for s > 0. For a function of several variables we denote by u(x, s) or L{u}(x, s) the Laplace transform of u(x, t) with respect to t.
The Laplace transform pair related to the Riemann-Liouville fractional derivative is provided u is continuous function satisfying u(0) < ∞ (see, e.g., [35], Chapter 1, Equation (1.29)). The Fourier transform of a function v(x), x ∈ R n , is given by The Fourier transform pair corresponding to the Laplace operator ∆ of a function v(x), x ∈ R n , such that lim |x|→∞ v(x) = 0, is (see, e.g., [36] Chapter 15) The two-parameter Mittag-Leffler function is defined by the series For β = 1 this is the one-parameter Mittag-Leffler function, E α (z) = E α,1 (z). The Mittag-Leffler function admits the asymptotic expansion (see, e.g., [35], Equation (E.30)) and satisfies the Laplace transform identity (see, e.g., [35], Equation (E.53)) The following relation, which can be easily verified directly from the definition (9), is often useful: Bernstein functions and related special classes of functions play a crucial role in the present work. The necessary definitions and basic properties are summarized below. For a detailed study on these classes of functions we refer to [37], see also ([38] Chapter 4).
The characterization of the class CMF is given by the Bernstein's theorem, which states that a function is completely monotone if and only if it can be represented as the Laplace transform of a non-negative measure (non-negative function or generalized function).
The class of Stieltjes functions (S F ) consists of all functions defined on R + which have the representation (see [25]) where a, b ≥ 0, ψ ∈ CMF and the Laplace transform of ψ exists for any λ > 0. A non-negative function ϕ ∈ C ∞ (R + ) is said to be a Bernstein function (BF ) if ϕ (λ) ∈ CMF ; ϕ(t) is said to be a complete Bernstein functions (CBF ) if and only if ϕ(λ)/λ ∈ SF .
Any of the four sets of functions defined above, CMF , SF , BF and CBF is a convex cone, i.e., if ϕ 1 and ϕ 2 belong to a specific set, then a 1 ϕ 1 + a 2 ϕ 2 belongs to the same set for all a 1 , a 2 ≥ 0.
For proofs and further details on these special classes of functions we refer to [37].

Representations as Generalized Diffusion-Wave Equations
We are concerned with the fractional Jeffreys' equation where τ q , τ T ≥ 0, 0 < α ≤ 1, as our first goal is to recast it into a generalized diffusion-wave equation with a memory kernel. For this purpose we apply Laplace transform to (15) and by using (7) obtain Here, and in what follows, the commutativity relation L (∆u) = ∆ (Lu) is used. For sufficiently well behaved functions u this property can be easily justified by the theorem for differentiation under the integral sign.
The next proposition shows that the two cases τ q < τ T and τ q > τ T should be considered separately.
The following assertions hold true for the function Proof. Let τ q < τ T and use the representation Since s α ∈ CBF by property (F), then s α + 1/τ T ∈ CBF and therefore, (D) implies (s α + 1/τ T ) −1 ∈ SF . Therefore, f (s) is a Stieltjes function for τ q < τ T as a linear combination with positive coefficients of two Stieltjes functions. The case τ q > τ T follows by interchanging τ q and τ T and taking into account property (D).

Generalized Diffusion Equation (τ q < τ T )
A generalized diffusion equation has the form [25][26][27] where ξ(t) ∈ L 1 loc (R + ) is a non-negative function, such that its Laplace transform ξ(s) ∈ SF . Applying inverse Laplace transform to (16) implies the generalized diffusion Equation (19) with kernel ξ(t), such that Proposition 1 states that ξ(s) is a Stieltjes function for τ q < τ T . Taking into account (18) and (11) we get from (20) the explicit form of the kernel ξ(t) where δ(t) is the Dirac delta function. Therefore, in the considered in this subsection case 0 < τ q /τ T < 1 the function ξ(t) is non-negative. In this way we proved that the required conditions on the kernel ξ(t) in Equation (19) are satisfied. Plugging the expression (21) for the kernel into the diffusion Equation (19) gives the following representation of the fractional Jeffreys' equation in the diffusion case

Generalized Wave Equation (τ q > τ T )
In the case τ q > τ T we are looking for a representation of Equation (15) as a generalized wave equation of the form [29,30] where η(t) ∈ L 1 loc (R + ) is a non-negative function, such that η(s) ∈ SF . In view of the balance Equation (3) the assumption (5) yields the following initial condition for the temperature distribution which implies With the help of this identity we rewrite Equation (16) as Therefore, the fractional Jeffreys' Equation (15) is equivalent to the generalized wave Equation (23) with kernel η(t), such that Applying inverse Laplace transform in (26) gives by the use of (11) Since τ q /τ T > 1 the kernel η(t) is a non-negative function. Moreover, η(t) ∈ CMF by property (G) and therefore η(s) ∈ SF . Hence, the requirements on the kernel η(t) are satisfied.
Alternatively, in the case τ q > τ T we can rewrite Equation (15) as the following integro-differential equation where µ(s) = 1+τ T s α 1+τ q s α . This follows by rewriting (16) in the form For the kernel µ(t) we get by applying (11) and Equation (28) admits the form Let us note that in this form the fractional Jeffreys' Equation (15) is studied in [24] in the context of viscoelastic flows of fractional Jeffreys' fluids.

Fundamental Solution
In this section we are concerned with the solution of the Cauchy problem for the fractional Jeffreys' heat conduction equation Problem (32), (33) and (34) is conveniently treated using Laplace transform with respect to the temporal variable and Fourier transform with respect to the spatial variables. By applying Laplace and Fourier transforms to Equation (32) and taking into account initial conditions (33), the boundary condition (34), and identities (7) and (8) we derive the solution in Fourier-Laplace domain where g(s) denotes the characteristic function Therefore, the solution of the Cauchy problem (32), (33) and (34) is given by the integral where G n (x, t) is the fundamental solution (Green function), defined in Fourier-Laplace domain as The behavior of the fundamental solution is studied next, with the emphasis on the one-dimensional case n = 1. The cases n = 2 and n = 3 are briefly discussed.
Let us note that in the special case τ q = τ T (i.e., g(s) = s) Equation (32) is the classical diffusion equation. The fundamental solution in this special case is the Gaussian function For the study of positivity of the fundamental solution in the general case we use Bernstein's theorem, stating that positivity of a function is equivalent with complete monotonicity of its Laplace transform. For this purpose, relevant properties of the characteristic function g(s) are derived next. Proposition 2. Let 0 < α ≤ 1. For any τ q , τ T ≥ 0 the function g(s) is a complete Bernstein function. If moreover 0 ≤ τ q < τ T then g(s) is a complete Bernstein function.
Proof. We use the relation g(s) = s f (s), where f (s) is defined in (17).
In the case τ q < τ T Proposition 1 implies f (s) ∈ SF . Therefore, g(s) ∈ CBF , according to the definition of the class CBF . Since square root of a complete Bernstein function is again a complete Bernstein function (by property (E) and the fact that 1 ∈ CBF ), this also gives g(s) ∈ CBF .
In the case τ q > τ T Proposition 1 states that f (s) ∈ CBF . Therefore, g(s) is a product of two complete Bernstein functions (s and f (s)) and property (E) implies g(s) ∈ CBF .
Let us note that g(s) / ∈ CBF for τ q > τ T . In fact, g(s) is not even a Bernstein function in this case. Indeed, for small s we have the asymptotic expansion g(s) ≈ s + (τ q − τ T )s α+1 and thus g (s) ≈ (τ q − τ T )(α + 1)αs α−1 > 0 for τ q > τ T , hence g(s) / ∈ BF . We also point out some general relations to the kernels ξ and η. First, in the diffusion regime, g(s) = s ξ(s) ∈ CBF if and only if ξ(s) ∈ SF . For the propagation regime we note that g(s) = s 2 η(s). Then the property η(s) ∈ SF implies that g(s) is a product of two complete Bernstein functions (s and s η(s)), therefore g(s) ∈ CBF .

One-Dimensional Solution
Consider the Cauchy problem (32), (33) and (34) in the one-dimensional case. The one-dimensional fundamental solution G 1 (x, t) is given in Fourier-Laplace domain by (37) with n = 1. By the Fourier inversion and using the well-known formula we get the Laplace transform of the fundamental solution of the one-dimensional problem Theorem 1. The fundamental solution G 1 (x, t) is a spatial probability density function evolving in time.
Proof. For the proof we use representation (39). First, cccording to Proposition 2 g(s) ∈ CBF ⊂ BF . Then property (B) yields exp −|x| g(s) ∈ CMF with respect to s > 0 for any x ∈ R considered as a parameter. On the other hand, g(s) ∈ CBF implies g(s)/s ∈ SF ⊂ CMF . The desired result, G 1 (x, s) ∈ CMF , follows by noting that the class CMF is closed under point-wise multiplication (property (A)). Therefore, by Bernstein's theorem, G 1 (x, t) ≥ 0. Further, (39) yields and, applying inverse Laplace transform we obtain The theorem is proved.
In the presented plots in the next section we see that G 1 (x, t) is a unimodal probability density function (p.d.f.) in the diffusion regime (τ q < τ T ), while in the propagation regime (τ q > τ T ) it is a bimodal p.d.f. The numerical results in the next section are based on an explicit integral representation of the fundamental solution, which we derive next.
To find explicit integral expression for the solution G 1 (x, t) we apply Bromwich integral inversion formula to (39), which yields By the Cauchy's theorem, the integration on the Bromwich path (γ − i∞, γ + i∞) can be replaced by integration on the contour D ∪ D 0 , where Indeed, the integrals on the contours {s = σ ± iR, σ ∈ (0, γ)} vanish for R → ∞ due to the following asymptotic expression Moreover, since lim s→0 s g(s) s exp st − |x| g(s) = 0, it follows that the integral on the semi-circular contour D 0 also vanishes. Therefore, integration on the contour D yields after letting ε → 0 To express the imaginary part under the integral sign in terms of elementary real functions we apply the formula for real and imaginary parts of the square root of a complex number and obtain after some standard manipulations the following result.

Theorem 2.
The fundamental solution G 1 (x, t) of the one-dimensional Cauchy problem admits the following integral representation for x ∈ R\{0}, t > 0: where the functions K ± (r) are defined by We note that the convergence of the integral in (41) is guaranteed by the following properties of the functions K ± (r): K ± (r) > 0, K ± (r) ∼ r (1−α)/2 as r → +∞ and K ± (r) ∼ r (1+α)/2 as r → 0.
The explicit integral representation (41) is used in the next section for numerical computation and visualization of the fundamental solution G 1 (x, t) for different values of the parameters. Let us note that another integral representation for the fundamental solution G 1 (x, t) can be found in ([23] Section 7.2.1). A comparison of some numerical examples confirms that the two representations are equivalent, see Figure 1a in the next section and Figure 7.14 in [23].

Mean Squared Displacement
Next, we study the temporal behavior of the mean squared displacement (MSD) in the one-dimensional case. Representation (39) implies for the MSD in Laplace domain Calculation of the integral yields where g(s) is the characteristic function (36). By the use of Laplace transform pair (11) we invert (42) and get two equivalent expressions in terms of Mittag-Leffler functions Both expressions are valid for all τ q , τ T ≥ 0. However we give the two different forms, since (43) seems more natural for the diffusion regime (τ q < τ T ) and (44) for the propagation regime (τ q > τ T ). Let us mote that the Mittag-Leffler functions in the MSD representations are positive functions, due to property (G) and the relation From the definition (9) of the Mittag-Leffler functions and their asymptotics (10) we derive the following asymptotic behavior for the MSD for short and long times The established asymptotic expansions show linear asymptotic behavior for short and long times. We also observe that the dominating term in the gradient of the MSD is 2τ T /τ q for t → 0 versus 2 for t → ∞. Therefore, in the diffusion regime (τ q < τ T ) the MSD increases faster near the origin than for large times. The opposite behavior is observed in the propagation regime (τ q > τ T ): the MSD increases slower near the origin than for large times. Let us note that qualitatively comparable asymptotic behavior of the MSD is observed in the fractional diffusion-wave equation, where MSD∼ t β with β ∈ (0, 1) in the diffusion regime and β ∈ (1, 2) in the propagation regime.

Numerical Examples
The integral representation (41) is used for numerical computation and visualization of the fundamental solution G 1 (x, t). The results for different values of the parameters covering the two regimes: diffusion and propagation, are given in this section. For the numerical computation of the improper integral in (41) the MATLAB subroutine "integral" is used. In all figures the fundamental solution G 1 (x, t) is plotted versus x, for x > 0. The part for x < 0 is symmetric, G 1 (−x, t) = G 1 (x, t). Figure 1 shows the evolution in time of the p.d.f. G 1 (x, t), starting from a delta function δ(x) at t = 0. The solution is plotted for five different time instances. In the diffusion regime (a) the maximum remains at t = 0, i.e., the p.d.f. is unimodal. In the propagation regime (b) the maximum moves away from the origin, the p.d.f. is bimodal.
In Figure 2 the solution is plotted for different values of the fractional parameter α ∈ (0, 1]. For α → 1 the solution approaches that of the Jeffreys' heat conduction Equation (1). For α → 0 the fractional derivatives become identity operators and (32) approaches the classical diffusion equation with the one-dimensional Gaussian as fundamental solution (see the plots for α = 0.05, which are qualitatively close to a Gaussian function).
In Figure 3 several plots are given with one and the same ratio of the relaxation times, in (a) It is seen that different plots correspond to one and the same ratio, that is the fundamental solution depends not only on the ratio of the relaxation times, but on their specific values. In all figures we observe behavior, typical for a diffusion process for τ q < τ T : the fundamental solution is monotonically decreasing in x for x > 0, representing a unimodal p.d.f. For τ q > τ T the behavior is typical for a propagation process, with a maximum moving away from the origin. In this respect there is a strong analogy with the fractional diffusion-wave equation with Caputo time-derivative of order β ∈ (1, 2) with the two corresponding regimes: diffusion (0 < β < 1) and propagation (1 < β < 2), c.f. ( [35] Figure 6.1), [31].

Subordination Principles and Multi-Dimensional Fundamental Solutions
In this section, we give one more standard representation of the fractional Jeffreys-type heat conduction Equation (15)-as a Volterra integral equation-and use some facts from the theory in the monograph [38]. As in the beginning of this work we consider again Equation (15) in general form without specifying the domain of the Laplace operator (finite or infinite spatial domain and type of boundary conditions). By the use of (7) Equation (15) reads in Laplace domain where T 0 (x) = T(x, 0). Taking the inverse Laplace transform in (45) yields the Volterra integral equation with kernel κ(t) satisfying κ(s) = 1/g(s), where the function g(s) is defined in (36). The explicit representation for the kernel is deduced by the use of the Laplace transform pair (11) as Let us note that κ(t) > 0. Indeed, according to (12), the first derivative of κ(t) is Therefore, for τ q < τ T the function κ(t) is decreasing from κ(0) = τ T /τ q > 1 to κ(+∞) = 1 and for τ q > τ T the function κ(t) is increasing from κ(0) = τ T /τ q < 1 to κ(+∞) = 1.
Let us look briefly at propagation speed of a disturbance in model (46). From general theory, see ( [38] Chapter 5), the velocity of propagation of the head of a disturbance is Therefore, c = ∞, except in the special case τ T = 0, τ q = 0, α = 1, in which we recognize the classical Cattaneo's equation. Therefore, with the exception of this specific case, a disturbance propagates with infinite speed, and therefore the Equation (15) is parabolic.
Let us point out that the symbol ∆ in (46) can represent different realizations of the Laplace operator, such that it is a closed and densely defined operator in the chosen space of functions and the wave equation W tt = ∆W is well posed. For details on abstract Volterra integral equations we refer [38]. Thanks to the properties of the function g(s) given in Proposition 2 above, we can apply Theorems 4.2 and 4.3 in [38] and formulate the following subordination principles. For completeness we briefly present the main idea and steps of the proof.
where the function W(x, t) is the solution of the classical wave equation and the kernel ϕ 1 (t, τ) is defined via the Laplace transform In the diffusion case τ q < τ T the following stronger subordination relation holds true where U(x, τ) is the solution of the classical diffusion equation and the kernel ϕ 2 (t, τ) is defined via the Laplace transform The functions ϕ j (t, τ), j = 1, 2, are probability densities in τ, that is Proof. By applying Laplace transform in (46) it follows where the function g(s) is defined in (36) Then application of the Laplace transform in (48) yields by the use of (49) and (54) Consider now the function T(x, t) defined by (50). The solution U(x, t) of the diffusion equation obeys in Laplace domain Therefore, application of the Laplace transform in (50) yields by the use of (51) and (56) with It remains to check properties (52). The proof is analogous to that of Theorem 1. To prove non-negativity, it suffices to establish that ϕ j (s, τ) ∈ CMF , s > 0. The nonnegativity of ϕ 1 (t, τ) is already proved in Theorem 1, ϕ 2 (t, τ) is treated in an analogous way, taking into account that g(s) ∈ CBF for τ q < τ T .
Further, from the definitions (49) and (51) of the kernels ϕ j (t, τ), j = 1, 2 we obtain by the use and the Bromwich inversion formula This implies Therefore, the kernels are normalized and properties (52) are established.
Theorem 3 provides a useful tool for the study of multi-dimensional problems for Equation (15) or problems on finite spatial domain by deriving the solution and relevant properties from known results concerning the corresponding diffusion and wave equations. An example of application of Theorem 3 is given next.
Proof. According to Theorem 3 subordination identity (50) holds true. Therefore where G * n is the Gaussian function given in (38), for which it is well known that it is a spatial p.d.f. Therefore, the non-negativity of G * n imply directly the non-negativity of G n and where for the last identity we have used (52). This proves the theorem.
To prove strictly that G 3 (x, t) can have negative values it suffices to show that G 3 (x, s) / ∈ CMF with respect to s (Bernstein's theorem). From (58) we deduce the asymptotic expansion which is not a completely monotone function with respect to s provided τ q > τ T . Indeed, in this case its first derivative ∂ ∂s { G 3 (x, s)} → +∞ as s → 0 and, thus, G 3 (x, s) as a function of s does not satisfy inequalities (13).

Concluding Remarks
The heat conduction equation with a fractional Jeffreys-type constitutive law is studied. By employing the Bernstein functions technique, two fundamentally different types of behavior are established, depending on the value of the characteristic parameter τ q /τ T : diffusion regime for τ q /τ T < 1 and propagation regime for τ q /τ T > 1. The one-dimensional fundamental solution is shown to be a spatial probability density function evolving in time, which is unimodal in the diffusion regime and bimodal in the propagation regime. The multidimensional fundamental solutions are probability densities only in the diffusion case, while in the propagation case they may have negative values.
The fractional Jeffreys-type heat conduction equation in the two different regimes is represented as a generalized diffusion, respectively wave, equation. The memory kernel in both cases is expressed in terms of Mittag-Leffler functions. The kernel is singular in the diffusion regime if α = 1 and regular in the propagation regime. This shows that one and the same model may support different types of kernels.
Subordination principles are formulated, which are useful for the study of multi-dimensional problems on finite or infinite spatial region.
All results hold true also for the case of classical Jeffreys-type heat conduction equation (α = 1). The established properties indicate a strong analogy between the fractional Jeffreys-type equation and the fractional diffusion-wave equation with the Caputo fractional time-derivative D β t u(x, t) = ∆u(x, t), 0 < β < 2, with its two different regimes: diffusion (0 < β < 1) and propagation (1 < β < 2). The technique developed in this work can be extended to the case of generalized diffusion-wave equation with a general memory kernel.