## Abstract

We show that the constitutive relation for the thermal flux proposed by Xu & Hu (2011) admits an unconditional instability. We also highlight the difference between mathematical models containing delay and those that include relaxation effects.

## 1. Introduction

In a recent article, Xu & Hu [1] put forth what they describe as a ‘ballistic-diffusive heat conduction model’ to relate ** q**, the heat flux vector, to

*T*, the absolute temperature. After a detailed and, seemingly, rigorous derivation starting from Boltzmann's transport equation, these authors arrive at the delay-containing constitutive relation [1], eqn (4.17): 1.1where the positive constants

*τ*

_{e},

*λ*(

*L*) and

*L*denote a

*delay*(or

*lag*) time, an effective thermal conductivity, and a macroscopic length scale, respectively, and

**=(**

*r**x*,

*y*,

*z*) is the position vector. According to Xu & Hu [1], p. 1858, ‘

*L*is the size of the heat conduction medium’.

Xu & Hu [1] then proceed to eliminate ** q** between equation (1.1) and the energy conservation equation [1], eqn (4.11), the result of which (see [1], eqn (5.1)) is
1.2a partial differential equation (PDE) known as the

*delayed heat equation*(see, e.g. [2] and the references therein). In equation (1.2), the positive constants

*ϱ*and

*c*

_{p}denote the mass density and the specific heat at constant pressure, respectively, and we observe that the coefficient of the Laplacian operator carries (SI) units of m

^{2}s

^{−1}. Note that [1], eqn (4.11), and all equations in [1] stemming from it, contain

*c*

_{v}, the specific heat at constant

*volume*, instead of

*c*

_{p}, the form of the specific heat used in this article. While

*not*incorrect in the case of many, if not most, solids, the use of

*c*

_{v}in the energy conservation equation is both contrary to standard practice and unnecessarily complicates the Laplacian's coefficient (see e.g. [3], p. 9; [4], §57).

Xu & Hu [1] next turn their attention to the problem of heat flow in thin films, for which they formulate a one-dimensional (1D) initial-boundary-value problem (IBVP) of the signalling type, where in the present setting ‘1D’ implies that ** q**=(

*q*(

*x*,

*t*),0,0) and

**∇**

*T*=(∂

*T*(

*x*,

*t*)/∂

*x*,0,0). However, rather than using equation (1.2) as the temperature equation for the thin film, which they regard as a thermally conducting rigid solid, these authors use [1], eqn (5.2): 1.3a hyperbolic PDE that many refer to as the

*damped wave equation*.

^{1}

As Xu & Hu [1] point out, equation (1.3) was obtained by expanding the left-hand side of the 1D version of equation (1.2) in a Taylor series, and then neglecting all but the first two terms, an approach often followed when applying non-classical diffusion models (see, e.g. [6,7], as well as [8], §1.5 and the references therein). Apparently unbeknownst to Xu & Hu [1], however, this attempt at approximation resulted in an equation for the temperature field (i.e. equation (1.3)) that, from the mathematical standpoint, is *equivalent* to the one which would have resulted had these authors assumed the constitutive relation known as the Maxwell–Catteneo (MC) law [8–10], namely,
1.4where *k*(>0) is the thermal conductivity and *τ*(>0) is the thermal *relaxation* time.

In the next section, we show that the Taylor-series approximation used to simplify the 1D version of equation (1.2) has resulted in the masking^{2} of a ‘fatal’ flaw contained in equation (1.1).

## 2. Mathematical analysis

In this section, we reconsider Xu & Hu's thin film problem under the assumption that the temperature field is described by the 1D version of equation (1.2); first, however, we introduce the following dimensionless variables:
2.1where *T*_{0}(>0) is the amplitude of the thermal signal which shall be applied to the boundary *x*=0, *Q*_{c}=*λ*(*L*)*T*_{0}/*L* is a characteristic value of the heat flux, *κ*=*λ*(*L*)/(*ϱc*_{p}) is the coefficient of thermal diffusivity, and *L* is now taken to represent the film's thickness. So, with a word of caution regarding differences in units vis-a-vis [1], IBVP (5.3), we now seek an analytical solution to the *delayed-*initial-boundary-value problem (dIBVP)
2.2a
2.2b
2.2cHere, *H*(⋅) denotes the Heaviside unit step function; *τ*_{0}=*τ*_{e}(*κ*/*L*^{2}) is the dimensionless thermal lag time; all asterisk superscripts have been suppressed but should remain understood; and for future reference we note that, in terms of the dimensionless variables defined above, the 1D version of equation (1.2) becomes
2.3

Owing to the nature of dIBVP (2.2), we elect to follow a dual transform approach in determining its solution. To this end, we first apply the finite sine transform [12] to equation (2.2a). This yields, after making use of the boundary conditions given in equation (2.2b), the subsidiary ordinary differential equation
2.4where a bar over a quantity denotes the image of that quantity under the finite sine transform. In turn, applying the Laplace transform to equation (2.4), and then making use of the delayed initial condition (dIC) given in equation (2.2c), we find, after solving the ensuing algebraic equation, that
2.5where *s* denotes the Laplace transform parameter and a hat over a quantity denotes the image of that quantity in the Laplace transform domain.

Now, by the Laplace inversion theorem [13],^{3}
2.6where the constant *γ*(>0) is chosen so that the line *s*=*γ* lies to the right of all the singularities of the integrand. As the singularities of the integrand consist only of poles, it follows that the integral in equation (2.6) can be evaluated using the residue theorem [13]. Determining these poles is, of course, equivalent to finding the roots of the transcendental equation
2.7Clearly, this equation is satisfied either when *s*=0 or when
2.8As shown in [2], the roots of equation (2.8) can be expressed in closed form. First, for *n*≠*n**:=*π*^{−1}(*eτ*_{0})^{−1/2}:
2.9all of which are simple poles. Second, for *n*=*n**:
2.10where except for , the single pole of order two, all are simple poles. Here, *W*_{r}(⋅) denotes the *r*th branch of the Lambert *W*-function [14], with *W*_{0}(⋅) denoting the *principal* branch.

On calculating and summing up the residues at each of the poles, it is a straightforward matter to show that the integral in equation (2.6) evaluates to
2.11Here, for simplicity of presentation, we have defined
2.12wherein, post-inversion, we performed the shifting *t*↦*t*−*τ*_{0} so that the dIC given in equation (2.2c) would be satisfied.

Finally, on applying [12], theorem 27 to equation (2.11) and simplifying, the exact solution to dIBVP (2.2) is found to be 2.13the substitution of which into equation (2.3) gives the corresponding expression for the flux, namely 2.14Equations (2.13) and (2.14) should be contrasted with equations (A 2) and (A 5) in Appendix A, which correspond to the IBVP for equation (1.3).

Now, an inspection of equation (2.12) for the case *n*≠*n** reveals that, in particular, both Re[*W*_{−1}(−*n*^{2}*π*^{2}*τ*_{0})] and Re[*W*_{0}(−*n*^{2}*π*^{2}*τ*_{0})] are *positive* for every *n*>*n*^{•}, where Re[⋅] and Im[⋅] denote the real and imaginary parts, respectively, of a complex quantity and *n*^{•}:=(2*πτ*_{0})^{−1/2}, where *n*^{•}>*n**. Therefore, , as , for every *n*>*n*^{•}; i.e. the solution to dIBVP (2.2) is strictly, and unconditionally, *unstable*. Criteria that ensure each of the remaining poles in equation (2.9) have a positive real part can also be derived, highlighting the severity of the instability.

## 3. Closing remarks

(i) While the delay-induced instability reported here is similar to the one studied in [2], there are two important differences to take note of: first, in dIBVP (2.2) a thermal excitation is delivered via the boundary, versus that of [2] wherein a pre-initial temperature profile evolves under the delayed heat equation and, second, whereas the instability exhibited by dIBVP (2.2) is unavoidable, given the manner of excitation, the instability encountered in [2] can be prevented, at the cost of restricting the value of

*τ*_{0}, by switching to a (problem-)specific class of pre-initial data.(ii) While [1], eqn (5.5) is correct for the value of

*c*_{0}(dimensionless lag time) chosen therein, it does not represent the*general*solution to [1], IBVP (5.3). Hence, for completeness, in appendix A, [1], IBVP (5.3) is reformulated in terms of the present dimensionless variables, its general solution is given, and the correct expression for the heat flux through the boundary*x*=0 is calculated.(iii) Along with the aforementioned instability, equation (2.13) also permits

*T*to assume negative values, a violation of the thermodynamic requirement*T*>0; again, see [2].(iv) Xu & Hu are correct in remarking that ‘the Maxwell–Cattaneo law violates Galileo's principle of relativity’ [1], p. 1852, and citing [15] in support, but

*only*to the extent that their remark is confined to equation (1.4). In 2009, Christov [16] presented the following*frame-indifferent*re-formulation of the MC law: , where represents a particular Lie–Oldroyd derivative. This generalization of*both*[15], eqn (20) and equation (1.4) can be obtained from the latter by replacing ∂/∂*q**t*with 3.1whereis the velocity vector, and it should be clear that reduces to ∂/∂*v**t*for a medium at rest; see also Straughan [8].(v) Xu & Hu state that ‘the Maxwell–Cattaneo violates the second law of thermodynamics [15]’ [1], p. 1861; however, as a reading of [15] of this communication will establish, no such claim is made by Christov & Jordan. For more on issues relating to the MC law's compatibility with the Clausius–Duhem inequality,

^{4}see [17,19] and the references therein.(vi) There are a number of models that should be considered as well-posed alternatives for obtaining a ‘non-Fourier heat conduction model … [that] takes into account of the memory effect and avoids the infinite heat propagation speed’ [1], p. 1863. For example, a number of models of heat conduction at finite wave speed(s) that account for memory of the heat flux are summarized in [9]; see also [8], §1.6 and the references therein. More recently, ‘higher gradient’ heat conduction models have also been proposed in this context by Christov [20] and Lebon

*et al.*[21].

## Funding statements

I.C.C. was supported by the National Science Foundation (NSF) under grant no. DMS-1104047 (at Princeton University) and by the LANL/LDRD Program through a Feynman Distinguished Fellowship (at Los Alamos National Laboratory). LANL is operated by Los Alamos National Security, L.L.C. for the National Nuclear Security Administration of the US Department of Energy under Contract No. DE-AC52-06NA25396. P.M.J. was supported by ONR funding.

## Acknowledgements

The authors wish to express their deep gratitude to (the late) Prof. C. I. Christov for both sharing his insight on non-classical heat conduction and his encouragement to complete the present communication.

## Appendix A. Exact solution to the thin film ‘signalling problem’ for the damped wave equation

The MC-based counterpart to dIBVP (2.2) involves the 1D version of equation (1.3); in terms of the set of dimensionless variables introduced in equation (2.1), this *well-posed* IBVP reads
A1a
A1b
A1cwhere all asterisk superscripts have been suppressed but should remain understood. IBVP (A 1) can be solved analytically using either an integral transform approach or via an eigenfunction expansion [22]. Either way, the exact, physically admissible, series solution (e.g. [22]) is
A2where, for simplicity of presentation, we have defined
A3and Δ_{n}:=1−4*τ*_{0}*n*^{2}*π*^{2}. This solution should be contrasted with [1], eqn (5.5), which, strictly speaking, is valid *only* for *c*_{j}≠0 (i.e. *τ*_{0} such that Δ_{n}≠0 ∀*n*).

To determine the corresponding expression for the flux, we must solve the MC law, which in the case of IBVP (A 1) assumes the (dimensionless) form *τ*_{0}∂*q*/∂*t*+*q*=−∂*T*/∂*x*, for *q*. Omitting the details, it is readily established that
A4which, on carrying out the indicated integration and simplifying, becomes
A5Here, *q*_{0}:=*q*(*x*,0), which is necessarily a constant in the case of IBVP (A 1), is arbitrary because the initial condition (∂*T*/∂*t*)|_{t=0}=0 is insufficient to determine its value uniquely. Hence, the flux through the boundary *x*=0 is easily found to be
A6which, even after setting *q*_{0}=0, does not appear to be equivalent to [1], eqn (5.6).

## Appendix B. Inversion of Laplace-domain solution via the series method

If, following Zwillinger [11] and setting aside (for the moment) the fact that , we Taylor-expand the right-hand side of equation (2.5) under the assumption *n*^{2}*π*^{2}|*e*^{−sτ0}*s*^{−1}|<1, then we obtain
B1

Setting *m*=*j*+1 and inverting term-by-term using any standard table of Laplace inverses, it is readily established that
B2which we observe is equivalent to the finite sum
B3Here, ⌈⋅⌉ denotes the ceiling function and, subsequent to inversion, we performed the shifting *t*↦*t*−*τ*_{0} so that the dIC given in equation (2.2c) would be satisfied. Now applying [12], theorem 27 to equation (B3), the exact solution to dIBVP (2.2) is found to be
B4

However, as the index *n* runs to , ∃*n*_{crt} such that *n*^{2}*π*^{2}|*e*^{−sτ0}*s*^{−1}|>1 for every *n*≥*n*_{crt}, where . This clearly renders the (naive) assumption on which going from equation (2.5) to (B1) was based invalid; it also highlights what, evidently, is another consequence of the unconditional instability exhibited by the solution of dIBVP (2.2).

## Footnotes

↵1 See, e.g. [5], wherein an in-depth analysis of equation (1.3) is presented.

↵2 When considered in the broader context of delay equation theory, this should not be surprising; as Zwillinger [11], p. 235 points out, ‘… the [Taylor series based] approximations that are obtained are often unrelated to the original [delay] equation’.

↵3 An alternative approach to inverting the Laplace-domain solution, and another representation of the exact solution

*T*(*x*,*t*), are given in appendix B.↵4 While generally regarded as a consequence of the second law of thermodynamics [17], the Clausius–Duhem inequality, which can be expressed as

⋅(*q***∇***T*)≤0, has been shown to imply a number of questionable results; see [18] and the references therein.

- Received August 16, 2013.
- Accepted October 17, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.