Causality in non-fourier heat conduction

We present a study of the causal response of non-Fourier heat conduction by introducing a dispersive generalized thermal susceptibility and show it obeys Kramers-Kronig relations. In particular, we discuss the Cattaneo-Vernotte and the dual-phase lag models. The constitutive equations are written as a linear response theory and determine the conditions for which the dual-phase lag model does not satisfy Kramers-Kronig relations. An alternative model is presented to describe non-Fourier transport using the similarity with the causal response of viscoelastic materials to avoid this unphysical issue.


Introduction
The Fourier law of heat conduction relates the temperature gradient ∇T and the heat flux q  , through the linear relation q r t T r t , , , where the thermal conductivity κ defines the response of the system. As written, this relation implies that there is an immediate heat flux when the temperature gradient is applied [1], or that the heat flux has a zero-time response. This unphysical effect, common in classical physics [2], was corrected by introducing a time delay or phase lag between the temperature and heat flux by Cattaneo and Vernotte (CV) [3][4][5]. The consequence of having a time lag is that the partial differential equation for the temperature field is no longer diffusive but becomes hyperbolic (telegraphers equation) and accepts wave-like solutions, known as thermal waves [6-8, 8, 9]. The CV equation assumes that once the temperature gradient is established, the heat flux occurs with a time delay or time lag, τ q . However, systems like metallic films or dielectrics, where electronphonon interactions and phonon scattering influence the microscopic behavior of heat transport and fasttransient processes [10], are not adequately described by the CV equation [11,12]. An alternative to CV is the dual-phase lag model (DPL). Besides the delay time τ q an additional delay time, τ T , is introduced in the temperature gradient [13][14][15][16], this is q r t T r t , , The non-Fourier behavior has been experimentally observed [15,17] in organic tissue such as pork muscle and blood [18,19], in porous systems such as metals, foams [20,21] and sands [22]. The wave-like behavior has been considered in many systems, exploiting similarities with other wave phenomena. In a periodic system, the appearance of a band structure in the dispersion relation is possible. For example, in temperonic crystals [23] a graphene layer is illuminated by a modulated heat source producing a periodic change in the temperature and a band structure. This band-structure is also obtained with a superlattice made of a regular array of materials for both finite periodic (thermal Bragg mirrors) and semi-infinite systems [24,25]. Another case is the use of thermal waves for cloaking employing core-shell spheres. By selecting the thermal conductivity of the core and the surrounding shell, the scattering cross-section of the thermal wave is reduced, obtaining a cloaking effect. This is achieved by applying Mie's theory of scattering to thermal waves [26][27][28].
Although the topic has been studied for several years, the possible applications at the meso and nanoscale [29]. Besides the Cattaneo-Vernotte and the dual-phase lag models, other models can be derived from a linear irreversible thermodynamic framework, such as the Guyer-Krumhansl equation, the Jeffreys and the Green-Nahdi type equations [30]. All of these falls in the general category of heat transport equations with a lagging behavior [31][32][33]. The DPL model has also been tested experimentally by Tang [34] that reported the temperature profile as a function of time in biological tissue. Within the experimental error, the DPL model provides a better description than the Fourier law. However, the DPL model can present some unphysical effects [35] that lead to negative temperature. Even considering higher orders in the delay time τ T , DPL model shows unphysical results depending on the initial conditions [36]. From the point of view of a strict thermodynamics approach, the DPL offers consistency with the second law provided that the two relaxation times τ T and τ q are positive [37]. Further studies have determined the possible values that the delay times τ p and τ q have to satisfy for the DPL equation to be consistent with the second law of thermodynamics and have stable solutions [38][39][40]. Despite the shortfalls of the DPL, it has been verified experimentally [11]. However, other systems require a three-phase lag model that has also been proposed to describe the temperature in biological systems where due to blood profusion, there is a displacement of the temperature gradient [41].
The non-Fourier behavior is an important topic in applications such as photo-thermal therapy, where it is important to understand the temperature profile around nanoparticles that are heated by an external source [32,[42][43][44]. In micro and nanofluidics also a non-Fourier behavior is observed [45,46]. At these size scales, the role of the dimension of the materials and the transport properties become important [47].
In non-Fourier problems, a missing discussion is a causal response between the heat flux and the temperature gradient as described by the Kramers-Kronig relations (KK). Causality given by the KK relations is known as strict causality, meaning that a system can not respond before being excited. This is not to be confused with causality in special relativity, which states that no signal can propagate faster than the speed of light [48,49]. In frequency space, there is a complex response or transfer function whose real and imaginary parts are not independent but related through the KK relations [50]. The most common use of these relations is in optical spectroscopy [51] where the refractive index and absorption coefficient are related. Thus, any measurement of the refraction index or attenuation can be verified for consistency if the KK relations are satisfied [52,53]. Another advantage of the causal relations is that it provides criteria to find permissible values for the electric permittivity and magnetic permeabilities [54].
The Kramers-Kronig relations extend to other linear systems. The real and imaginary parts of the complex compliance of viscoelastic materials also satisfy causal relations. The real part represents the in-phase elastic response and the imaginary part the out-of-phase response related to the viscous energy loss in the material [55,56]. For the case of acoustic waves, causality implies a relation between the reciprocal of the wave group velocity and the attenuation [57,58].
The issue of the temporal response in thermal processes was discussed by Vázquez et al [59] analyzing the causality principle and the connection between different mechanisms contributing to the thermal transport for hyperbolic and parabolic equations. To give closure to the issues surrounding the causality of non-Fourier heat conduction, in this paper, we discuss the Kramers-Kronig relations for the Cattaneo-Vernotte equation and the dual-phase model that are the most common approaches to describe causal heat conduction.

The heat conduction as linear response theory
The problem of heat conduction can be cast as a linear response equation where the stimulus is the temperature gradient T r t , ( )   and the causal response is the heat flux q r t , ( )   needed to restore thermal equilibrium. This can be written as the integral relation: The specific form of H t t ( ) -¢ depends on the constitutive equations between the heat flux and the temperature gradient [59]. For the Fourier law of heat conduction, we have where κ is the thermal conductivity. Evaluating equation (1) we have The memory function H F describes the instantaneous response of the system to the thermal perturbation. The Fourier law is an example of a linear response that does not obey causality.
By introducing a delay time, τ q , between the heat flux and temperature gradient, the response function [1,59] is (1), H CV is the solution of the Cattaneo-Vernotte (CV) equation given by The additional term q t q  t ¶ ¶ depicts the thermal inertia needed to restore thermal equilibrium. In the limit τ q → 0 we recover the delta function characteristic of Fourier's law. Otherwise, we observe that system's memory decays exponentially depending on the relaxation time value. The CV model removes the unphysical result of an infinite propagation speed of thermal signals, and instead, it predicts a speed of heat propagation given by v The other case is known as the dual-phase lag model [13][14][15][16]61], it describes heat conduction assuming two thermal relaxation times, τ q and τ T . In this model, causality is also imposed in the energy conservation equation. From this point of view, the heat flux can be a response of the system to restore the thermal equilibrium after a gradient of temperature is imposed, or it can be the cause of the thermal disturbance [12]. The way of taking into account τ T is by modifying the CV equation as follows Considering small values of τ T we expand the right side of equation (6): If τ q > τ T , the heat flux is the response system as happens with the CV model, while when τ q < τ T , the heat flux acts as the excitation.
The memory function of dual-phase lag for τ q > τ T is The memory function when τ q < τ T is obtained exchanging, in equation (1), τ q by τ T , and the heat flux by the temperature gradient. The dual-phase lag model is also known as the Jeffreys equation [1,62].
For each one of the models, we introduce the time variable where we use the memory function or kernel as with the step function defined as Θ(η) = 1 for η > 0 and zero for η < 0. The introduction of this function allows us extend the lower integration limit. Thus, for the three heat conduction models (Fourier, CV, and DPL), we get the following kernels: where the heat flux is the response system, and , we apply it on equation (9), and considering the convolution theorem we have as the complex generalized thermal susceptibility [60].
When we solve either the CV equation or the DPL equation, there are two choices for the frequency and wavevector: (i) to assume that the frequency takes complex values while the wave vector is real [23,61], in general, this is useful for infinite and semi-infinite systems; (ii) to assume a real frequency and complex wavevector, generally used for finite systems. We adopt the second case since it is a more common approach in wave propagation phenomena, particularly problems involving boundaries [24,43].
Taking the Fourier transform of the kernels ( )  h we obtain for Fourier's law and for the DPL model The last equation, with the change τ q by τ T and κ by 1/κ, corresponds to the dual-phase lag model where the temperature gradient is the system response; this is In figure 1, we plot the real and imaginary parts of the thermal susceptibility for the CV model (solid lines) as a function of ωτ q and for the DPL model (dashed lines) with the ratio τ T /τ q = 0.5.
The CV response (equation (14)) is square-integrable [63], this is where C is a constant. This is equivalent to requiring that On the other hand, the dual-phase lag model is not square integrable but satisfies the weaker condition Taking the limit in equation (15), we obtain C T q 2 2 ( ) k t t = .
In the usual notation, this limit is where is redefined the real part as˜  ¢ -¥ , while the imaginary part stays as.
According to Titchmarsh's theorem [64] a square-integrable dispersive function must be causal. A consequence of this theorem is that the real and imaginary parts of the response function obey Kramers-Kronig relations. The KK relations describe how the imaginary () and real (¢) parts of the generalized susceptibility are related if the kernel of the model is causal and analytic [58]. The KK relations are given by where P denotes the principal value of the integral. For the CV equation 0  = ¥ , and for the dual phase-lag model The proof of the Kramers-Kronig relations is well established, see, for example [58,65,66].
If we consider the real part of the thermal susceptibility for the DPL model, the CV case is recovered for τ T = 0. We want to verify the validity of the equation (21), thus the improper integral must be equal to. It is computed using the Cauchy integral theorem for the simple poles on the upper halfplane [65] , i.e., z 1 = ω and z 2 = i/τ q , and the respective residues Res f z , 1 1 The contour of integration is the upper-half plane as indicated in figure 2; crosses indicate the two poles.
Thus, the imaginary part is and has be proven the second KK relation (equation (21)). Similarly, its real part can be derived from the imaginary part of the generalized thermal susceptibility. The CV susceptibility also satisfies the KK relations by taking, in the previous equation, the limit τ T → 0 and 0  = ¥ . As further proof of causality and the KK relations, it is common to use the Cole-Cole (Argand) diagram representation [67]. This diagram is a parametric plot of˜( ) for different values of the frequency. In figure 3, we show the Cole-Cole diagram for the CV model and for the DPL model for different values of the ratio r = τ T /τ q . The Cattaneo case corresponds to r = 0; we observe that the curve intercepts the real axis ( k ¢ ) at 0 and 1. We identify the last values as the static case 1 0  k = , and the arrows indicate the two crossings. The DPL case is characterized by 0 < r < 1, for these curves, we observe that the first crossing with the horizontal axis is at k ¥ and the second one is always at 0  k. The limiting case, ¥ , for each r is indicated by a cross. The highest value of corresponds to ωτ q = 1 given by 2 The two delay times are not independent, and r is bounded by [40] r 2 3 2 3 -< < + . These bounds are needed to have stability in the DPL and to avoid unphysical effects. The case r = 0 corresponds to the Cattaneo equation. All the curves for r ≠ 0 satisfy Kramers-Kronig relations, including the one corresponding to r = 0.1 that does not fall within the desired bounds. This shows that satisfying KK relations is a necessary condition but is not sufficient to warrant that the chosen parameters for the DPL equation are physically acceptable.
At the value of τ q = τ T , the curves become a point and k ¢ = , the static weight of the thermal susceptibility, while the imaginary part is zero. We observe that for this case, the dual-phase lag model does not satisfy the KK relations; in other words, it breaks causality. Figure 2. Contour of integration, where the limit R → ∞ has be taken. Crosses indicate the location of the poles. The CV and DPL (for τ q > τ T ) models have the same poles. In case τ T > τ q , the pole along the y axis will be located at ζ = i/τ T .
for U 0 the initial condition [12]. Equation (27) is the diffusion equation with a source term that decays exponentially. Thus, for τ T = τ q , we obtain a case of diffusion heat conduction, which does not satisfy the Kramers-Kronig relations.
Additional comments on the dual-phase lag model Linear response functions similar to those described here are standard in other fields, such as viscoelastic materials. These materials can deform after a force is applied but tend to restore their initial shape. Phenomenologically, their behavior is modeled using arrays of springs that describe the elastic behavior and dashpots for the viscous part [68].
The most straightforward representation of an elastic body consists of a spring with elasticity constant E; by the Hooke law, the relation between the stress σ and the strain ε is σ = Eε, where the spring has the function to restore the mechanical equilibrium. This case is analog to Fourier's law of heat conduction, where we associate the stress to the temperature gradient, this is σ → ∇T, the strain with the heat flux, q  e  and 1/E with the thermal conductivity, as shown in figure 4(a). The viscous behavior is added through dashpots which have the effect of opposing deformation. An example is the Kelvin-Voigt model, which consists of a spring and a dashpot in parallel; both components feel the same deformation while the resulting stress of each element is added. The spring tends to restore the equilibrium while the dashpot presents a resistance or dissipation; this last is measured by the viscosity η. For simplicity, we obviated the tensorial nature of the stress and the strain. Even when the applied stress is removed, the tendency for the strain to continue deforming the body has its counterpart in heat conduction with the so-called thermal inertia introduced in the Cattaneo model; we show a depiction in figure 4(b). Finally, we have Zener's model of viscoelastic materials that includes two relaxation mechanisms; the stress dashpot is relaxed and provides a deformation in the whole system represented by a Based on these similarities and the fact that the DPL model does not satisfy KK relations only for τ T = τ q , we propose that equation (7) can be modified assuming that the thermal conductivity κ for the static case is different in the high-frequency limit, say κ T . This will modify the thermal susceptibility of equation (15) as The physical interpretation for this proposal is that the heat conduction has a different relaxation mechanism when the delay τ T is included, as happens with the inclusion of a spring in series in the Zener configuration. At this point, this mechanism is not certain, but this modified version of DPL has a well-behaved dispersion behavior even for τ q = τ T and also the CV equation is recovered when τ q → 0. The KK relations provide consistency of a physical model; this is the basis of our proposal for the modified DPL model. There is a caveat to this approach. Modifying dpl ( )  w introducing a static conductivity κ and a high-frequency one, κ T , can be done if there is consistency with thermodynamics and will require an analysis of whether the possible values of the relaxation times τ q and τ T are still valid as done in [38][39][40]. We emphasize that by analogy we follow the definition of Dragoman [69] that there is a resemblance between the equations although the phenomena are different. These similarities have been proven useful in the past, as is the case of quantum and electromagnetic analogies leading to the concept of photonic crystals [70].
As mentioned in the introduction, there are other time-lag model equations, such as the Guyer-Krumhansl (GK) equation [30] that involve second-spatial derivatives of the heat flux. In a one-dimensional case, it reads The additional term γ 2 is related to processes that preserve phonon momentum [62]. Because of the spacederivative of the heat flux, the corresponding memory function will have a spatial dependence. This is, the complex thermal susceptibility will be of the form k , ( )   w , where k  is a wave vector that has to be defined. It means we will have a nonlocal response function. However, the formalism presented in this paper will not be applicable, and a different approach for analyzing the Kramers-Kronig relations will have to be developed. To our knowledge, the nonlocal KK relations have only been studied in the case of optical response [71]. Since the CV equation is embedded into the GK equation, and it preserves the same time-dependence form, it will be expected that the GK model should also satisfy the KK relations.

Conclusions
The Cattaneo-Vernotte equation and the dual-phase model are examples of non-Fourier heat conduction. In this work, we showed that the CV and DPL models could be written as a linear response equation with a frequency-dependent generalized complex thermal susceptibility,. From this complex susceptibility, we derived the Kramers-Kronig relations that relate the real and imaginary parts needed to satisfy strict causality. The DPL model presents a problem when the time delay of the heat flux τ q and the temperature gradient τ T are equal. It does not satisfy the KK relations, and it has been shown [38][39][40] that it leads to stability issues. To recover causality for τ q = τ T , a modified version of the DPL model is introduced following the empirical models for viscoelastic materials.
There are many of non-Fourier models, and none shows a definite answer of which one is the most suitable. In the case of a spatial dependent susceptibility, usually named nonlocal response, the appropriate causal relations are to be derived. Unlike Cattaneo's model, which is hyperbolic, the DPL and GK equations are parabolic. However, by including additional terms in the Taylor expansion of the delay time, there is a 'hyperbolic' version of these models [62], whose response kernels and corresponding complex thermal susceptibilities are yet to be determined.
Although the experimental evidence must be the decisive factor for rejecting any model, it is mandatory to verify that any model satisfies the basic requirements of physical consistency given by the Kramers-Kronig relations.