Stability analysis for linear heat conduction with memory kernels described by Gamma functions

This paper analyzes heat equation with memory in the case of kernels that are linear combinations of Gamma distributions. In this case, it is possible to rewrite the non-local equation as a local system of partial differential equations of hyperbolic type. Stability is studied in details by analyzing the corresponding dispersion relation, providing sufficient stability condition for the general case and sharp instability thresholds in the case of linear combination of the first three Gamma functions.


Introduction
The present paper is based on three keywords: diffusion, memory and stability. The starting ingredient is a scalar quantity u satisfying where x, t ∈ R and v is the flux. Equation (1) has to be completed with an additional relation for the unknown v. The classical Fourier law (for heat conduction) prescribes that the flux is proportional to the space derivative of the unknown u, (2) v(x, t) = −α ∂ x u(x, t).
The couple (1)- (2) gives raise to the well-known heat equation, which, being parabolic, has a number of appealing features (well-posedness, smoothing effects...) At the same time, it also exhibits some deficiencies, the most striking being the infinite speed of propagation of disturbances. To mend such a flaw (and for other motivations which are not discussed here) alternative laws for the flux v has been proposed.
The most famous, suggested by Cattaneo in [1], assumes that the equilibrium between flux and gradient, istantaneous in (2), is in fact described by a dynamical relation of relaxation type, namely (3) τ ∂ t v(x, t) + v(x, t) + α ∂ x u(x, t) = 0 where τ > 0 is an additional parameter describing the relaxation timescale. The couple (1)- (3) gives raise to a damped wave equation for the unkwnon u. Hyperbolicity of the model permits to trespass many of the flaws of the heat equation. The Cattaneo's law (3) can be rephrased into a more general framework in which the flux is not given by the istantaneous value of the gradient of u, but it is rather an average of its history. Precisely, equation (3) is equivalent to the assumption that the flux v satisfies the convolution relation (4) v( with a kernel g of exponential type, i.e. g(t) = α e −t/τ /τ . As stated in [9], it would be a miracle if for some real conductor the relaxation kernel could be rigorously represented by an exponential kernel with a single time of relaxation, as is required by Cattaneo's model. Infact, the exponential choice can be considered as a stratagem which permits to trade the non-local relation (4) with the local relation (3), as a consequence of the property (5) τ dg dt + g = 0, τ g(0) = α.
The conservation equation (1) coupled with the memory law (4) gives raise to the non-local equation which is meaningful for a class of memory kernel g containing the pure exponentials as a special case. The Fourier law (2) can be recovered assuming that the kernel g is a Dirac delta distribution concentrated at t = 0. Considering relation (6) for modelling heat conduction has been proposed many years ago in [8], and it has received an increasing interest in recent years 2 (among others, see [2,4,7] and references therein). Usually, it is assumed that the function g is a function with values in [0, ∞) and such that (7) α : The value of α describes the intensity of the overall memory effect.
A key question is: under which assumption on the memory kernel g, the dynamics is stable, viz. there is no exponential growth of perturbations of constant states? In the two basic examples, Fourier (2) and Cattaneo (3) laws, stability holds as a consequence of a significant localization of the kernel g at t = 0. How is the situation for memory kernels exhibiting a wider spreading in the past history? Paraphrasing [6], does having a lasting memory cause the appearance of nontrivial patterns? Interesting answers has been given in [2,4], based on a semigroup theory approach. The spirit is to determine properties on the function g which guarantees stability of the associated semigroup. Over-simplifying (and suggesting to the interested reader to look at the originals), the first article provides a necessary condition for stability, namely g (t + s) ≥ C e −δt g (s) for t ≥ 0, s > 0 for some C ≥ 1, δ > 0, which turns to be equivalent to g + C g ≤ 0 in (0, ∞) for some C > 0. The second shows its sufficiency in the class of convex kernels g.
Here, the philosophy is different: rather than considering general memory kernels, we analyze in details a specific class, extending the pure exponentials, and still permitting a translation of the non-local term (4) in a set of local PDEs, so that the complete model turns to be equivalent to a linear hyperbolic system. For such class of systems, we analyze the stability problem by means of a frequency analysis.
The form of the memory kernels is suggested by the fact that variations of the property (5) are satisfied by a more general class of functions: the Gamma distributions. These form a two-parameter family defined as follows: given a positive real number τ and an integer j ≥ 1, the Gamma distribution with shape j and scale τ is A sketch of the graph of the Gamma functions for the first values of j (and τ = 1) is depicted in Figure 1. Shape j equal to 1 corresponds to the exponential choice. Such functions are also defined for j real, but this case is not considered here for a reason that will be clear in a few lines. In what follows, the functions g(·; j, τ ) will often be denoted by g j , omitting the dependence on τ .
The choice of memory kernel given by Gamma distributions is common in the analysis of delay differential equations, but, apparently more rare for PDEs. Indeed, apart for the fact that different values j correspond to different shape of the kernel and, in particular, different dependencies on the past history in the diffusive dynamics, there is no specific physical motivation for such a choice.
The main advantage is merely technical and it is encoded in what H. Smith calls the linear chain trick (see [12]). It consists in the fact that for such kernels it is possible to get rid of the memory term and to transform the non-local equation into a finite set of partial differential equations, exactly as in the case of the exponential memory term relation (6) is equivalent to the hyperbolic system (1)-(3).
Presentation of the main results. To enter into the details, given a natural number k, let g be the vector with components (g 1 , g 2 , . . . , g k+1 ) where g j are defined in (8), by I k the k × k identity matrix and by N k+1 the (k + 1) × (k + 1) matrix with 1 on the subdiagonal Then, there holds where e 1 = (1, 0, . . . , 0), to be considered as an extension of (5). As a consequence, for any kernel g given by linear combinations of Gamma distributions, it is possible to convert the convolution relation (4) into a set of differential equations for some auxiliary functions, whose linear combination, furnishes the flux v given in (1). Precisely, given τ > 0, k ≥ 1 integer and θ = (θ 1 , . . . , θ k+1 ) ∈ R k+1 , if the memory kernel g is given by then the non-local equation (6) is equivalent to a system for (u, v) ∈ R k+2 System (10), which is a generalization of the Cattaneo diffusion equation (recovered for k = 0), turns to be (non strictly) hyperbolic whenever θ 1 > 0 (see Section 2). The choice (9) for the kernel g allows the presence of a single time-scale τ which is from now on normalized to 1. The compresence of multiple time scales is meaningful, but it is not considered here. The goal of this paper is to analyze in details the stability/instability properties of (10) for different ranges of the parameter θ, by means of a detailed description of the dispersion relation (11) p(λ, µ) := det λu µθ Any couple (λ, µ) such that (24) holds corresponds to a solution of (10) with the form (u, v) = (U, V)e λt+µx . Then, the stability character of the model (10) is encoded in the sign of the real part of the elements of the sets Λ(ξ) for ξ ∈ R, where Such Fourier analysis can be translated into the description of the L 2 −decay in time of solutions with data in L 2 ∩ L 1 following a standard procedure, based on the use of inverse Fourier transform. General results concerning the behavior of the dispersion relation for small and large frequencies are encoded in the following statement.
Theorem 1.1. Let the memory kernel g in (4) be of the form θ·g for some θ ∈ R k+1 with θ 1 > 0 and θ j ≥ 0 for any j.
I. [small frequencies] There exist r, δ > 0 such that for |ξ| < r the set Λ(ξ) ∩ {λ : Reλ ≥ −δ} is composed by a single branch λ = λ 0 (ξ) and there holds and all of the roots of the polynomial have strictly negative real parts.
If both small and large frequencies stability holds, in order to estabilish the complete dynamical property of the system, it is necessary to analyze the behavior of the dispersion relation for intermediate values of the frequency. Specifically, it is relevant to analyze the existence/non existence of couples (λ, µ) with the form (iζ, iξ) such that (11) is satisfied. In case of non-existence of such couples, collecting with the properties that holds in case of small and large frequencies stability, one can state that there exists c 0 > 0 such that The validity of such estimate guarantees that the Shizuta-Kawashima condition is satisfied and thus the model exhibits dissipation (see [11] and descendants). This paper furnishes two additional results in this direction. The first one gives sufficient condition for stability for general values of k.
Theorem 1.2. Let the memory kernel g in (4) be of the form θ·g for some θ ∈ R k+1 with θ 1 > 0 and θ j ≥ 0 for any j. If the coefficients θ j are such that then there exists c 0 > 0 such that (14) holds.
Some of the kernels g for which the hypotheses of Theorem 1.2 hold change convexity for t ∈ (0, ∞). Indeed, for k = 1 there hold so that, when θ 1 < 2θ 2 , g is negative in (0, 2 − θ 1 /θ 2 ) and positive otherwise. In particular, for θ 2 < θ 1 < 2θ 2 , the kernel g is non-convex and nevertheless determine a stable dynamics. Similar properties hold for general k and θ j small for j > 2. Let us observe that memory kernel derived by homegenization as in [3], while exhibiting decreasing monotonicity, appear to be non-convex.
The second result is specific to the case k = 2 and furnishes a complete characterization for the values of θ guaranteeing stability.
Then, assuming θ 1 > 0 and θ 2 , θ 3 ≥ 0, estimate (14) holds if and only if (17) θ 2 < θ 1 and either 3θ 3 < 2θ 1 , The conditions on the parameters θ 1 , θ 2 , θ 3 define three regions corresponding to three different behaviors. Setting η 2 = θ 2 /θ 1 and η 3 = θ 3 /θ 1 , the region of high frequency stability is η 2 < 1. In this strip, the set where condition (14) holds is where the second set describes the interior of an ellipsis (see Figure 2). Theorem 1.3 indicates the existence of a parameter region where the instability is of "Turing type", in the sense that large and small frequencies are damped by the equation while a window of intermediate frequencies exhibit exponential growth in time. In this respect, this paper intends to contribute to the strand of the analysis of instabilities driven by the presence of memory (see [10,5] and descendants).
The organization of the paper is the following. Section 2 starts by discussing the linear chain trick for Gamma distributions. Then, it concerns with the hyperbolicity of the system (10) in term of the coefficients θ j . In Section 3, we analyze the dispersion relation of the system dividing the small and large frequencies regimes, giving the complete proofs of Theorems 1.1 and Finally, in Section 4, we analyze intermediate frequencies giving the proof of Theorem 1.2 and then concentrating on the case k = 2, providing the proof of Theorem 1.3.

Linear chain trick and hyperbolicity
For the reader's convenience, following closely [12], this Section starts with the presentation of the linear chain trick and its use to reduce the diffusion equation (1) with the flux v defined in (4) to a system of partial differential equations, in the case of memory kernels given by a linear combination of Gamma functions (8).
Fixed the scale τ (normalized to 1, without loss of generality) and given T > 0, let us consider the linear functionals defined for bounded functions φ ∈ C((−∞, T ]; R).
Proof. For j ≥ 2, integrating by parts, we infer Thus, since I 1 = 1, there holds I j = 1 for any j. Therefore we infer Moreover, each g j solves an appropriate Cauchy problem, viz.
Then, differentiating the definition ψ 1 := T 1 ψ 0 we obtain Similarly, differentiating ψ j := T j ψ 0 for j ≥ 2, we infer To conclude the proof, it has to be shown that the problem has at most one solution, viz. the one just exhibited. Assuming that ψ 1 j and ψ 2 j both solve the prescribed linear system, the difference Ψ j := ψ 1 j − ψ 2 j is a bounded solution to Then, Ψ 1 is clearly zero and, by induction, also Ψ j is zero for any j ≥ 2.
Thanks to Proposition 2.1, the non-local problem can be reduced to a system of partial differential equations when the memory kernel g has the form (9) for some natural k and θ ∈ R k+1 . Indeed, there holds so that the problem is equivalent to the system System (19) can be rewritten in vectorial form as where W = (u, ψ 1 , . . . , ψ k+1 ) and (21) where θ * = (θ 2 , . . . , θ k+1 ) and The hyperbolicity character of the system is recognized analyzing the principal term in (20) and, specifically, the spectral properties of the matrix A −1 0 A 1 .
Lemma 2.2. The matrix A −1 0 A 1 has real eigenvalues if and only if θ 1 ≥ 0. Moreover, the matrix A −1 0 A 1 i. if θ 1 = 0, has a single eigenvalue 0 with algebraic multiplicity equal to n + 2 and geometric multiplicity equal to n; ii. if θ 1 > 0, has three eigenvalues: ± √ θ 1 (each with multiplicity equal to 1) and 0 (with algebraic and geometric multiplicity equal to k).
Proof. The eigenvalue problem is (A 1 − λA 0 )W = 0. Therefore, since the statement on the eigenvalues and the algebraic multiplicity holds. Moreover, for λ = 0, the eigenvalue problem reduces to the conditions which define a vector space of dimension k in R k+2 .
As a consequence of Lemma 2.2, it is reasonable to concentrate the attention on the case θ 1 > 0, since only in such a case the principal part of the system is diagonalizable and thus the Cauchy problem is well-posed.
A diagonalizing matrix can be determined by computing the k + 2 independent eigenvectors. Since the eigenvalue problem reads as −λu + where δ ij is the Kronecker symbol. Such eigenvectors, considered as column vectors, define a change of coordinates C such that, setting Z = C −1 W , the principal part of the system (20) is diagonalized

Small and large frequencies analysis
Applying Fourier-Laplace transform, system (20) is converted into the system  is satisfied. Such relation can be read in different ways. Here, following the standard approach dictated by the Fourier transform, the variable µ will be considered as a purely imaginary number, i.e. µ = iξ with ξ ∈ R, and the attention is directed toward the corresponding values of λ ∈ C such that (23) is satisfied. Thanks to the very special structure of the matrices defining the polynomial p, it is possible to determine an explicit formula for it. Let A 0 , A 1 and B as in (21)-(22). Then, there holds Proof. The formula can be proved by induction on k. Indeed, for k = 0 the relation holds, as can be seen by direct computation. Then, making explicit the dependence on n in the dispersion relation denoting the polynomial by p k , there holds Then, assuming (24) for k − 1, we infer The proof is complete.
The analysis of the relation p(λ, µ) = 0 gives a complete information on the stability properties in the model. First of all, the behavior for small values of µ, corresponding to the large time-behavior, in case of stability, is very simple. Indeed, setting µ = 0, there holds p(λ, 0) = λ(λ + 1) k+1 .
Thus, there are k +1 modes that are dissipated with exponential rate (corresponding to the zero −1 with order k + 1) and there is a single mode with a weaker decay.
In order to explore in details the slow-mode, relative to the value λ = 0, the expansion λ = Aµ + Bµ 2 + o(µ 2 ) can be plugged into the dispersion relation, obtaining Thus, in a neighborhood of (λ, µ) ≈ (0, 0), the dispersion relation defines the curve This completes the proof of Theorem 1.1, part I. Next, we investigate the behavior of the couples (λ, µ) satisfying the dispersion relation, in the regime µ → ∞. Keeping track of the highest powers of λ and µ gives Thus, as µ → ∞, setting κ := √ θ 1 , the dispersion relation defines a set separated into three branches, given by the asymptotics which corresponds to the eigenvalues of the principal part of the system (20). The stability properties of such branches are encoded in the subsequent term in the asymptotic expansions for λ as a function of µ. Setting λ = ±κ µ + C + o(1) with C to be determined, and plugging into the dispersion relation, one obtains Similarly, choosing λ = C + o(1) and inserting in the dispersion relation, one infers −µ 2 k j=0 θ k+1−j (C + 1) j + o(µ 2 ) = 0 so that the coefficient C satisfies the algebraic relation k j=0 θ k+1−j (C + 1) j = 0.
Such equality can be rewritten as Exchanging the order of the sums, it follows that C is root of the polynomial This completes the proof of Theorem 1.1, part II.
Rephrasing, necessary and sufficient conditions for stability in the high frequency regime are θ 2 < θ 1 and the request that the polynomial Q is a (strict) Hurwitz polynomial, i.e. a polynomial whose roots have strictly negative real parts. Verification of the latter condition can be done by using the Routh-Hurwitz stability condition. Let us consider some explicative examples, when θ 1 + · · · + θ k+1 = 1. We list the explicit form of the polynomial Q and the corresponding Routh-Hurwitz array for the first integers and for specific choices of the parameters θ j : for k = 2 and general θ j , In particular, large frequency stability holds for k = 2 (any case), k = 3 and θ 1 > 1/9, k = 4 and θ 1 > 1/5. For higher values of k, the necessary and sufficient conditions appear difficult to be translated in a simple condition that is directly readable from the values of the coefficients θ k . Nevertheless, for θ 1 sufficiently large with respect to the other coefficients θ j , j = 2, . . . , k + 1, the stability condition is satisfied.
Proof. The polynomial Q can be rewritten as In the case θ j = 0 for any j = 1, the polynomial reduces to (1 + x) k and, thus, it has a single root x = −1 with algebraic multiplicity equal to k. Changing the values of θ j , the location of the roots changes smoothly. Let θ j be such that there exists y ∈ R for which Q(iy) = 0. Then, there holds Then, we infer the estimate In particular, as soon as the condition (25) is satisfied, no roots of the polynomial Q may intersect the imaginary axis.

Intermediate frequency analysis
Next, we pass to the analysis of the dispersion relation (23) (with p explicitly given in (24)) in the intermediate frequency regime, under the assumption that the system is stable in the large frequency regime. In particular, we assume θ 2 < θ 1 .
For λ = −1, the relation can be rewritten as Let (iζ, iξ) with ξ = 0 be such that the above equality is satisfied; then there hold so that, given ζ satisfying the first relation, ξ is such that if and only if the term on the righthand side is positive. Both the real part of (1−iζ) j and the imaginary part of (1 − iζ) j /ζ contain only even power of ζ; hence the above relations can be rewritten in term of the single variable s = ζ 2 .
Proof of Theorem 1.2. We claim that if (15) holds, then the first relation in (26) has no solutions. Indeed, the first condition in (26) can be rewritten as under the condition (15), the relation (28) cannot be satisfied.
Combinations of the first three Gamma distributions. Condition (15), stated in Theorem 1.2, is sufficient, but not necessary. Here, we restrict the attention to the case k = 2, that is we concentrate on kernels g of the form (29) g(t) = θ 1 + θ 2 t + 1 2 θ 3 t 2 e −t . with the aim of determining the sharp threshold separating stability and instability. Indeed, in this situation, it is possible to determine a complete characterization of the choices of the coefficients θ 1 , θ 2 , θ 3 giving raise to stable dynamics, as encoded in Theorem 1.3.
Proof of Theorem 1.3. Setting s = ζ 2 , the first equality in (26) becomes Assuming θ 2 < θ 1 (necessary condition for high-frequency stability), equation (30) has positive roots if and only if which can be equivalently rewritten as If these conditions are satisfied, for s satisfying (30), the expression (27) for ξ 2 reduces to having used the relation Thus, the unknown ξ is explictly given by where s = ζ 2 is a positive root of (30).
A graphical representation of the different stability regimes relative to the values of parameter θ 1 , θ 2 , θ 3 is given in Figure 2, where the regions lie in the plane (η 2 , η 3 ) with η j := θ j /θ 1 .
In the analysis of diffusion equations with memory, a typical requirement on the kernel g concerns with its monotonicity and convexity properties. Thus, we continue the study of the case (29) determining the parameter range in which such assumptions are satisfied. Let us set Our aim is to compare S with the regions M and C of monotonicity and mononotonicity + convexity, respectively.

C ⊂ S ⊂ M.
A graphic representation of the regions S, C and M is provided in Figure 3.
Thus, for kernels of the form (29), stability is guaranteed by convexity. This is coherent with the result in [4], taking into account that the function µ corresponds to −g . Indeed, such result is relative to convex kernels (µ non-increasing) and apply to the present case since the "flatness rate" (see details in the above reference) for kernel of the form (29) is always zero.
Also, still in the case of kernel of the form (29), non-monotonicity implies instability. Differently, we can state that monotonicity is a necessary condition for stability for kernels (29), consistently with Theorem 1.5 in [2] (proved for a general class of memory kernels).
Finally, the analysis exhibits a large region of decreasing non-convex kernels divided into two parts by the stability/instability transition line. In particular, there exist (monotone) kernels g which are stable even if not convex.