Localized heat perturbation in harmonic 1D crystals. Solutions for an equation of anomalous heat conduction

TL;DR

This paper derives exact solutions for anomalous heat propagation in 1D harmonic lattices, revealing how different initial temperature perturbations decay over time and highlighting observable peaks near the wavefront.

Contribution

It provides explicit solutions for various initial perturbations in 1D harmonic crystals, advancing understanding of anomalous heat conduction mechanisms.

Findings

Decay near wavefront proportional to 1/√t

Decay at center proportional to 1/t

Peaks remain visible near the wavefront over time

Abstract

In this work exact solutions for the equation that describes anomalous heat propagation in 1D harmonic lattices are obtained. Rectangular, triangular, and sawtooth initial perturbations of the temperature field are considered. The solution for an initially rectangular temperature profile is investigated in detail. It is shown that the decay of the solution near the wavefront is proportional to $1/ \sqrt{t}$. In the center of the perturbation zone the decay is proportional to $1/t$. Thus the solution decays slower near the wavefront, leaving clearly visible peaks that can be detected experimentally.

Equations

$$\dot{T}=\beta T^{\prime \prime },$$

$$q=-\kappa \nabla T,$$

$$\tau \ddot{T}+\dot{T}=\beta T^{\prime \prime },$$

$${\ddot{u}}_k={\omega }_e^2(u_{k-1}-2u_k+u_{k+1}),{\omega }_e\stackrel{\text{mbox}def}{=}\sqrt{\frac{C}{m}},$$

$$u_k{|}_{t=0}=0,\dot{u}{|}_{t=0}=\sigma (x){\rho }_k,$$

$$k_{B}T=m⟨\dot{u_k}{⟩}^2,$$

$$\ddot{T}+\frac{1}{t}\dot{T}=c^2{T}^{\prime \prime },$$

$$\dot{T}{|}_{t=0}=0,T{|}_{t=0}=T_0(x).$$

$$T(t,x)=\frac{1}{\pi }\underset{-t}{\overset{t}{\int }}\frac{T_0(x-c\tau )}{\sqrt{t^2-{\tau }^2}}\mathrm{d}\tau .$$

$$T_0(x)=\{\begin{array}{ll}0,& x<-l,\\ \Phi (x),& -l<x<l,\\ 0,& x>l,\end{array}$$

$$T_0(x)=A\left(\mathcal{H}\left(x+l\right)-\mathcal{H}\left(x-l\right)\right),$$

$$\mathcal{H}(x)=\{\begin{array}{ll}0,& x<0,\\ 1,& x\ge 0.\end{array}$$

$$T(t,x)=\frac{A}{\pi }\underset{-t}{\overset{t}{\int }}\frac{\mathcal{H}(x+l)}{\sqrt{t^2-{\tau }^2}}\mathrm{d}\tau -\frac{A}{\pi }\underset{-t}{\overset{t}{\int }}\frac{\mathcal{H}(x-l)}{\sqrt{t^2-{\tau }^2}}\mathrm{d}\tau .$$

$$T(x,t)\stackrel{\text{mbox}def}{=}{T}_S(x,t)=\{\begin{array}{ll}0,& x\le -ct,\\ \frac{A}{\pi }\arccos \left(\frac{x}{ct}\right)& -ct\le x\le ct,\\ A,& x\ge ct\text{par}\end{array}$$

$$t\le {\tau }_0:T(x,t)=\{\begin{array}{ll}0,& l+ct\le x,\\ \frac{A}{\pi }\arccos (\frac{x-l}{ct}),& l-ct\le x\le l+ct,\\ A,& 0\le x\le l-ct,\end{array}$$

$$t\ge {\tau }_0:T(x,t)=\{\begin{array}{ll}0,& ct+l\le x,\\ \frac{A}{\pi }\arccos (\frac{x-l}{ct}),& ct-l\le x\le ct+l,\\ \frac{A}{\pi }\left(-\arccos (\frac{x+l}{ct})+\arccos \frac{x-l}{ct}\right),& 0\le x\le ct-l,\end{array}$$

$$\dot{T}=\beta T^{\prime \prime }.$$

$$T(x,t)=\frac{1}{2}\mathrm{erf}\left(\frac{x}{\sqrt{4\beta t}}\right),$$

$$T(x,t)=\frac{1}{2}\mathrm{erf}\left(\frac{x+l}{\sqrt{4\beta t}}\right)-\frac{1}{2}\mathrm{erf}\left(\frac{x-l}{\sqrt{4\beta t}}\right).$$

$$T(t,0)=\frac{A}{\pi }\left[\pi -2\arccos \left(\frac{l}{ct}\right)\right]=2\epsilon +O({\epsilon }^3),$$

$$T(t,-l+ct)=T(t,l-ct)=\frac{A}{\pi }\left[\pi -\arccos \left(\frac{2l}{ct}-1\right)\right]=2\sqrt{\epsilon }+O\left({\epsilon }^{\frac{3}{2}}\right),$$

$$T(t,0)\stackrel{t\to \infty }{\sim }2\epsilon \sim \frac{1}{t},T(t,-l+ct)=T(t,l-ct)\stackrel{t\to \infty }{\sim }2\sqrt{\epsilon }\sim \frac{1}{\sqrt{t}}.$$

$$T_{\mathrm{env}}(x)=\frac{A}{\pi }\left[\pi -\arccos \left(\frac{2l}{x+l}-1\right)\right].$$

$$T(\xi ,t)=\frac{A}{\pi }\arccos \left[1-\frac{l}{ct}(\xi +1)\right]\stackrel{t\to \infty }{\sim }\frac{A}{\pi }\sqrt{\frac{2l}{ct}}\sqrt{\xi +1},$$

$$T(\xi ,t)=$$

$$\stackrel{t\to \infty }{\sim }$$

$$T=\frac{A}{\pi }\sqrt{\frac{2l}{ct}}F(\xi ).$$

$$T_0(x)=\{\begin{array}{ll}0,& x<0;\\ Bx,& x\ge 0,\end{array}$$

$$T(x,t)=Bx\left(\frac{1}{\pi }\arcsin x+\frac{1}{2}\right)+\frac{B}{\pi }\sqrt{t^2{c}^2-x^2}\stackrel{\text{mbox}def}{=}f(x),$$

$$T_0(x)=\{\begin{array}{ll}0,& x<-l,\\ (x+l)B,& -l\le x<0,\\ (-x+l)B,& 0\le x<l,\\ 0,& x\ge l.\end{array}$$

$$T(t,x)=T_L(t,x-l)+T_L(t,x+l)-2T_L(t,x),$$

Full text

1 Introduction

Nowadays investigation of nonlinear thermomechanical processes in low-dimensional structures attracts high interest due to the rapid development of nanoelectronical devices based on materials with microstructure [1, 2, 3, 4]. Achievements in nanotechnology allowed for an experimental proof of the wave nature and finite propagation velocity of thermal perturbations [5, 6]. This can provide a foundation for a universal theory of thermal conduction, applicable both on micro and macroscales.

The classical heat equation is a parabolic partial differential equation that describes the distribution of heat in a given region over time,

[TABLE]

where $T$ is temperature, $\beta$ is thermal diffusivity, dot $\dot{(\leavevmode\nobreak\ \leavevmode\nobreak\ )}$ denotes differentiation with respect to $t$, prime $(\leavevmode\nobreak\ \leavevmode\nobreak\ )^{\prime}$ denotes differentiation with respect to $x$. The classical heat equation is derived on the basis of the Fourier’s law [7, 8],

[TABLE]

where $\kappa$ is thermal conductivity, $q$ is heat flux, and $T$ is temperature. Practical applications show that at the macroscale Fourier’s law of heat conduction is well applicable in order to describe temperature processes. However, it predicts an infinite speed of signal propagation, which is paradoxical from a physical point of view. A study of processes at the microscale, when the characteristic length is proportional to several atomic bond lengths, requires us to use more complicated models of heat transfer, which take the finite velocity of heat propagation into account. Well observed abnormalities from the Fourier’s law happen in thermal processes occurring in one-dimensional crystalline structures [9]. Recent experimental works show length dependence of thermal conductivity of nanostructures [10]. Significant deviations from Fourier’s law were shown for $\mathrm{C}$ and $\mathrm{BN}$ nanotubes  [11]. Thermal anomalies for nanoscale structures can be used in practice for designing perspective devices such as thermal diodes [4].

The anomalous nature of heat processes for one-dimensional lattices was demonstrated analytically in paper [12], where a problem of heat flow between two heat baths was considered. A hyperbolic heat equation is one of the alternatives to describe heat processes, which take finite speed of temperature propagation [13, 14] into account,

[TABLE]

where $\tau$ is a relaxation time. However, Eqn. (3) has serious difficulties in describing heat transfer in one-dimensional crystals, since no unique relaxation time can be determined [15]. A perspective approach for description of un-steady heat processes in 1D crystals is presented in papers [16, 17, 18]. By using correlational analysis, the initial stochastic problem for individual particles is reduced to a deterministic problem for statistical characteristics of the crystal motion. Finally, a continuum equation (7) describing anomalous heat transfer in 1D harmonic lattices is obtained in paper [17]. In the current work exact solutions for this equation will be obtained. Solutions for a number of problems such as rectangular, triangular, and sawtooth initial perturbations will be obtained. Properties of the solution for rectangular initial perturbation such as decay behavior and asymptotics of the wavefront will be investigated. These results can be used for analysis of the anomalous heat transfer in more complex systems, such as 1D crystals on elastic foundation [19] and 2D-3D crystals [20]. An understanding of the anomalous heat conduction is important for analysis of the experimental results, which are to be obtained in the nearest future due to the rapid development of nanotechnologies.

2 Localized perturbations in a harmonic chain

The harmonic chain is a simple and powerful model in order to investigate anomalous heat conduction phenomena. Following on to the work [17] let us consider an infinite harmonic chain. Each particle with mass $m$ is connected to its neighbor by Hookean springs with stiffness $C$. The equation of motion of the particles reads:

[TABLE]

where $u_{k}$ is displacements of particle with index $k$. The following initial conditions are considered:

[TABLE]

where $\rho_{k}$ are independent random variables with zero expectation and unit variance; $\sigma$ is the variance of the initial particle velocity. The variance is a slowly changing function of the spatial coordinate $x=ka$, where $a$ is the initial distance between neighboring particles. Such initial conditions can be realized by ultrafast heating, for example with a laser [21]. Let us introduce the kinetic temperature $T$ as

[TABLE]

where $\langle...\rangle$ is an operator averaging over realizations, and $k_{B}$ is the Boltzmann constant. In paper [17] a continuum partial differential equation for the kinetic temperature was obtained:

[TABLE]

where $c$ is the speed of sound in a one-dimensional crystal. Eqn. (7) describes the evolution of the spatial temperature distribution in the chain. The following initial conditions for the equation [17] corresponds to stochastic initial conditions (5):

[TABLE]

The solution of the initial problem (7)–(8) can be obtained in integral form [17]:

[TABLE]

Eqn. (7) is a particular case of the Darboux equation [22]. This type of equation was investigated earlier in context with spherical averages for solutions of 2D and 3D wave equations. However, it was not investigated well in connection with the problems of heat conduction. Eqn. (7) looks similar to the hyperbolic heat equation (3), however, it has a variable coefficient. This peculiarity is due to anomalous heat transfer in a 1D chain. From the form of Eqn. (7) it seems that it has a singularity. However, it does not matter because Eqn. (7) is to be solved together with the initial conditions (8). The absence of singularity is confirmed by the general analytical solution (9) and solutions of particular problems that will be considered below.

This work is dedicated to finding exact analytical solutions of Eqn. (7) for cases when the initial thermal distribution $T_{0}(x)$ is a localized function of coordinate $x$,

[TABLE]

where $\varPhi(x)$ is an arbitrary function, and $l$ is the half width of the localized perturbation. Experimentally such an initial temperature distribution (10) can be realized by performing superfast laser heating of a localized region of the chain.

3 Rectangular perturbation

3.1 Solution

Let us consider the case when the initial temperature perturbation has a rectangular shape:

[TABLE]

where $\mathcal{H}(x)$ is the Heaviside function:

[TABLE]

$A$ is the amplitude of the temperature perturbation. After substituting formula (11) into the solution (9) we obtain:

[TABLE]

By substituting the solution for a single Heaviside initial impulse, which was obtained in [17],

[TABLE]

to (13) we obtain the solution of the given problem as a linear combination of these solutions. Solution for positive $x$:

[TABLE]

[TABLE]

where $\tau_{0}=l/c$. For negative $x$ the solution is symmetric and can be obtained by $T(x,t)=T(-x,t)$.

For comparison let us consider the same initial problem for the classical heat equation:

[TABLE]

The solution for an initial Heaviside step temperature perturbation has the form [23]

[TABLE]

where $\operatorname{erf}(x)$ is the Gaussian error function. Then solution of the initial problem (8), (11), (17) is:

[TABLE]

The time evolution plots for the solution of anomalous heat equation (7) and the Fourier equation (17) are shown in Fig. 1. Let us compare the two solutions. The Fourier solution is forming a peak at $x=0$ which decays exponentially. For the case of the anomalous heat equation the solution decays in the area near $x=0$ more rapidly than near the wavefronts forming two peaks. The peaks travel in negative and positive directions with coordinates $x=-l+ct$ and $x=l-ct$.

3.2 Decay behavior

Let us consider the decay behavior of the solution (16) at $x=0$. We perform a series expansion of the solution,

[TABLE]

where $\varepsilon=\frac{l}{ct}$ is a small parameter.

Now let us consider the decay behavior of the peaks $x=l-ct$ and $x=-l+ct$. From formula (16) it follows that:

[TABLE]

Summarizing the above:

[TABLE]

Thus, the solution decays faster in the area between wavefronts (proportional to $1/t$) rather than near the wavefront (proportional to $1/\sqrt{t}$). Thus the peaks remain strongly pronounced even for long times.

3.3 Envelope curve for the peaks

The solution (16) has two peaks. The peaks travel in positive and negative directions at speed $c$. Since the solution is symmetric let us consider only the peak with the coordinate $x=ct-l$. We shall consider the curve drawn by the peak of the solution as it travels in positive direction. By substituting $t=\frac{x+l}{c}$ into formula (21) we obtain the expression for the enveloping curve:

[TABLE]

For any $x$ we have: $T(x)\leq T_{\mathrm{env}}\left(|x|\right)$.

The enveloping curve is shown in Fig.2. The expression decays as $1/\sqrt{x}$, which corresponds to the fact that the solution decays as $1/\sqrt{t}$ near the wavefront (the wavefront travels at constant speed).

3.4 Asymptotic behavior of the wavefront

Let us consider the solution (16) near the wavefront at long times $t$. Let $\xi=\frac{-x+ct}{l}$. For $x\in[-l+ct;l+ct]$ we have:

[TABLE]

where $\cfrac{\xi+1}{ct}$ is a small parameter, used for expansion. For $x\in[l-ct;-l+ct]$ we have:

[TABLE]

The functions (24) and (25) have the following structure:

[TABLE]

The plot of the solution (16), expressions (24) and (25) with the corresponding dimensionless time parameter $t/\tau_{0}=100$ is shown in Fig. 3. The relation (26) means that the shape of the solution shrinks vertically with time, but it does not change horizontally. The asymptotic solutions (24) and (25) give peak values of $T$ for $\xi=1$, where $F(\xi)=\sqrt{2}$. Thus the solution at this point is continuous, but not smooth (the derivative of the solution has a jump).

4 Triangular perturbation

In order to obtain the solution for a triangular initial function we consider the following auxiliary problem where the initial temperature distribution is a linearly heated semispace,

[TABLE]

where $B=A/l$ is a constant of proportionality. After substituting (27) in (9) we obtain the solution for $|x|<ct$,

[TABLE]

and for $|x|>ct$ the initial temperature distribution is preserved.

Now we consider the problem for a triangular initial heat perturbation, which can be expressed by the following piecewise function:

[TABLE]

The solution for the initial temperature distribution (29) will be a linear combination of solutions for the linearly heated semispace. Denote the solution (28) by $T_{L}$. The solution for the initial distribution (29) will then be as follows:

[TABLE]

the solution is symmetric, $T(x,t)=T(-x,t)$. The part corresponding to positive $x$ has the following piecewise form:

[TABLE]

[TABLE]

[TABLE]

The plot of the solution for the triangular initial perturbation is shown in Fig. 4.

Unlike the solution for a rectangular initial perturbation, which has a wavefront with vertical tangent and infinite derivative and a break of the temperature profile at the peak, the solution for a triangular perturbation has a smooth beginning at the wavefront and smooth behavior at the peak.

5 Sawtooth perturbation

We consider an initial heat perturbation as a sawtooth spatial temperature distribution. It can be written in the following form:

[TABLE]

The initial conditions (34) can be written as linear combinations of step function and linearly heated semispace. Then the solution for a sawtooth initial perturbation (34) can be obtained from the corresponding combination of the solutions for a step initial distribution $T_{S}(x,t)$ and a linearly heated semispace initial distribution $T_{L}(x,t)$:

[TABLE]

it has the following piecewise form:

[TABLE]

[TABLE]

The plot of the solution is shown in Fig. 5. The left wavefront has a smooth beginning and an infinite derivative at the peak. On the other hand, the right wavefront has an infinite derivative and vertical tangent at the beginning, smooth behavior at the peak, and a horizontal tangent and zero derivative at the peak.

6 Conclusions

The process of heat transfer in a 1D infinite harmonic chain was investigated. Localized initial perturbations were considered. Solutions for an equation describing anomalous heat conduction (7) derived in [17] were obtained. Exact analytical solutions for rectangular, triangular, and sawtooth initial impulses were considered. It was shown that solutions for (7) unlike solutions for classical heat equation have a strongly pronounced wavefront. For the rectangular case it was shown that the decay of the solution near the wavefront is proportional to $1/\sqrt{t}$. Near zero the decay is proportional to $1/t$. Thus the solution decays slower near the wavefront, leaving clearly observable peaks. The shape of the wavefront is described by a function inversely proportional to the square root of time and has the form $T=\cfrac{A}{\pi}\,\sqrt{\cfrac{2l}{ct}}\,F(\xi)$.

The solution for a triangular initial temperature perturbations has a smooth beginning at the wavefront and a smooth behavior at the peak.

In case of a sawtooth initial perturbation we have a non-symmetrical solution. The left wavefront has a smooth beginning and an infinite derivative and vertical tangent at the peak; the right wavefront has an infinite derivative and vertical tangent at the beginning, smooth behavior at the peak, zero derivative and a horizontal tangent at the peak.

The obtained solutions demonstrate the wave behavior accompanied with power decay. This differs from the results obtained from the solutions of the classic heat equation (1) (diffusive behavior, exponential decay) and the hyperbolic heat equation (3) (wave behavior, exponential decay). Such properties of the obtained solutions can be applied for analysis of the experimental data and choosing the right model for the description of the heat processes.

Acknowledgement

The work was supported by Russian Science Foundation [Grant No. 14-11-00599].