Validity of local thermal equilibrium in anomalous heat diffusion

Local thermal equilibrium (LTE) is a general presumption of many theoretical analyses in nonequilibrium statistical physics. It describes a situation that although the system is not in global thermal equilibrium, each small portion of the system may still be described approximately by the laws of thermal equilibrium. The validity of LTE has however seldom been investigated carefully. Here, by studying the ensemble velocity distribution and its spatial correlation, we present strong evidences for the lack of LTE in anomalous heat diffusion processes in one dimensional harmonic lattices and the Fermi–Pasta–Ulam-β lattices. In particular, clear nonzero excess kurtosis and long range correlations have been observed, with values scaling linearly with the initial temperature difference. Therefore near thermal equilibrium is not sufficient to achieve LTE. Some existing studies that are based on the existence of LTE should be revisited. Using the same methods, we also show that LTE is still valid in the ϕ4 lattice, in which heat diffuses following Fourier’s law.


Introduction
Temperature is the pivotal concept of thermodynamics and statistical physics, which predominantly affects almost all fields of natural science, ranging from material properties, phase transition, chemical reaction to weather, climate, as well as our life process. However, From a statistical physics point of view, the definition of temperature is not applicable to most of the systems and environment that we encounter in both our daily life and scientific experiments. Since the definition has a solid basis only for systems in thermal equilibrium states [1] -following the axiomatic zeroth law of thermodynamics-in which all properties reaches spatial and temporal uniform [2].
For a canonical ensemble of classical Hamiltonian systems, H(p i , q i ), temperature is the unique external quantity characterizing the phase space distribution ρ∝e −H/T , knowing as the canonical distribution (we choose a temperature scale so that k B =1). This in turn leads to the generalized equipartition theorem [2], i.e. T H c x = á ¶ ñ x , where ξ is any canonical coordinate of the system (p i or q i ); and c á ñ · denotes the canonical ensemble average under the measure ρ. The choice of ξ=p i for any i then gives a formula for the temperature commonly referred as the kinetic temperature. It has long been established that the only possible steady state of a system in contact with heat baths at the same temperature is a thermal equilibrium state depicted by the canonical distribution, leading to a well defined temperature for the system. With differing heat bath temperatures, however, the system is nonequilibrium in the long run and a canonical distribution is no longer achievable, prohibiting a rigorous definition of temperature for the system [3]. Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
To deal with such nonequilibrium systems, the concept of local thermal equilibrium (LTE) is introduced, with the hope that each small portion of the system can be treated as if it is in a locally thermal equilibrium state with a local temperature (LT), T(x, t), varying slowly in space and/or in time-so that the existing theories describing global thermal equilibrium can still be applied approximately to each small region [4]. Based on it, other slowly varying thermal properties of the system can then be defined locally as well; and various transport phenomena induce by this thermal nonequilibrium can be studied, such as the diffusion of heat. Microscopically, LTE also implies that the generalized equipartition is approximately satisfied and LT can still be obtained using the kinetic temperature definition T p i i 2 = á ñwith á ñ · denoting the average over the nonequilibrium ensemble under study.
Regardless of its ubiquitous applications in studying nonequilibrium phenomena, the validity of LTE has rarely been investigated from a fundamental basis. Empirically, its validity in three-dimensional bulk materials, as well as the validity of subsequent Fourier's law, is supported by the repeated experimental [5,6] and numerical observations [7,8]. In the past two decades, the discovery of breakdown of Fourier's law in one-dimensional (1D) Fermi-Pasta-Ulam (FPU) lattice [9] has inspired numerous studies and debates on anomalous heat transport in low-dimensional lattices, both numerically and theoretically [10,11]. Those studies were typically setup in nonequilibrium steady states maintained by heat baths at different temperatures. For homogeneous systems, smooth local kinetic temperature profiles were generally observed and the existence of LTE was not questioned. An counterexample is the FPU lattice with alternating masses-when attached to heat baths with different temperatures, the LT profile exhibits oscillation, manifesting the lack of LTE [12].
In this paper, we study numerically the validity of LTE in heat diffusion in 1D homogeneous systems without heat baths, starting from well prepared nonequilibrium states with LTE and undergoing relaxation towards equilibrium. These systems are thus in nonequilibrium transient states governed by the Liouvillian dynamics. Unlike the studies in [12] that the systems are in stationary states, the systems that we shall study are in nonstationary states. By introducing varies measures, we demonstrate that the initial LTE condition is not capable of ensuring LTE afterwards in anomalous heat diffusion under general settings, even for systems with homogeneity and arbitrarily close to thermal equilibrium. For normal heat diffusion obeying Fourier's law, however, LTE is plausible.

Model and setup
We study 1D homogeneous lattices with general Hamiltonian where the masses are set to unity without loss of generality, V and U depicts the interparticle interaction and onsite potential, respectively. We focus more on the FPU-β lattice with V x x x 1 2 2 4 4 = + b ( ) and U(x)=0. Very recently, a quantitative dynamical characterization of the degree of equipartition of a given microscopic state has been worked out for this model with however different parameters [13]. For comparison, we also study a harmonic lattice with V x x = + ( ) . These systems do not exhibit thermal expansion and thus macroscopic work is absent between adjacent parts, which make us more concentrated in the subject under study.
We simulate heat diffusion in these lattices in two stages following [14]. Specifically, in the first stage we prepare the system in a nonequilibrium steady state by attaching a Langevin heat bath at temperature T i to each of the particles i. The temperatures follow a Gaussian profile: Throughout the study, we fix the background temperature T 0 =1 and σ=4. Such a temperature is far beyond the equipartition threshold of the FPU-β lattice, which would then prohibit the generation of ordered orbits in stage two (see below) and the system will evolve towards thermal equilibrium autonomously [15,16].
In order to make sure that the lattice is in its linear response region during the diffusion process, we apply a small value of ΔT=0.1, unless otherwise stated. Such an 10% relative temperature difference is smaller than most of the values applied in existing heat conduction and heat diffusion simulations, in which LTE was commonly presumed.
In this stage, the particles follow Langevin dynamics: The parameter γ is tuned to make the kinetic temperature distribution of the lattice as close to the target distribution equation (3) as possible. The resulting initial heat current is generally also maximized. Our choice is γ=3, which is typically more or less the best. After a long enough warm-up period, the system eventually reaches a nonequilibrium steady state. Then we define time t=0 and start the second stage by removing all the heat baths so that heat starts to diffuse. The system evolves autonomously towards equilibrium according to Liouvillian dynamics, i.e. equation (4) with the second bracket removed. We define the local energy and local heat current as , where E E i c 0 = á ñ stands for the average energy density in thermal equilibrium at the background temperature T 0 =1 [14].

The kurtosis of velocity distribution
The immediate quantity that needs to be checked for a valid LTE state is the distribution of the single particle velocities. If a part of a system is in a LTE state, then the velocity distribution should be Gaussian [12,17]. A direct check of the single-particle velocity distributions reveals that they are still in bell shape. The symmetry of the velocity distributions is guaranteed by the symmetry of the system, leading to a vanishing skewness. A standard criteria to measure the deviation of a symmetric distribution from Gaussian is thus the kurtosis κ: Comparing to Gaussian distribution, a positive κ implies a narrower central peak and two fatter tails; and vice versa. Kurtosis profiles in the FPU-β lattice for β=1 at different time t are plotted in figure 1(a). At time t=0, the kurtosis of the particles are uniformly zero, implying a potential LTE state maintained by the Langevin heat baths with close temperatures. Nonvanishing kurtosis at t=0 can be observed only for very large ΔT-a phenomenon not surprising at all since the large variation of temperature along the lattice is likely to destroy the precondition of LTE.
During the diffusion process with t>0, however, the kurtosis quickly changes from zero to nonzero and forms two side peaks that spread from the center. The pictures are partly similar to those for the local excess energy distributions ΔE i that are plotted in figure 1(b). This can be understood from the theory of nonlinear fluctuating hydrodynamics which predicts Kardar-Parisi-Zhang [18] like sound modes and Levy heat modes, see [19] and references therein for detail. However, some distinctions between ΔE i and κ are also clearly observed. Firstly, the central part of the kurtosis profile forms a deep well at the beginning and then diminishes rapidly; in contrast with the slowly decaying central peak of excess energy distribution. In other words, while energy is still concentrated in the center and gradually spread towards the two ends of the lattice, the velocity of the particles in the central part approaches Gaussian much faster. Secondly and more interestingly, the two side peaks of the kurtosis profile moves faster than those of the ΔE i . The speed of the heat propagation in a lattice is determined by the speed of the heat carriers, possible candidates include solitons and renormalized phonons for the FPU-β lattices. The peak of the energy distribution should move at the sound speed in the lattices. Here we see that the information of an upcoming heat diffusion is actually propagated supersonically and outruns the heat itself. The change in the velocity distribution is thus an herald of heat diffusion which can be detected before the heat wave is reached. As a side effect, we also observe that when ΔE i (t) reaches its local maximums, the kurtosis κ crosses zero and changes sign.
To rule out the possibility that the nonzero kurtosis profile is manifested only because the initial temperature difference is too large so that we are actually in the far-from-equilibrium region, we study the role of the initial ΔT on the velocity distribution. we simulate heat diffusion for various ΔT and plot the rescaled kurtosis, κ/ΔT, in figure 1(c). It is clearly seen that the curves for ΔT=0.1, 0.3 and 1.0 overlap with each other, indicating that κ/ΔT approaches a constant distribution in the small ΔT limit. In other words, as ΔT approaches zero, κ and ΔT are infinitesimals of the same order and therefore its deviation from Gaussian distribution is still observable in the near-equilibrium limit, just as the presence of heat current under the influence of temperature gradient as described by Fourier's law.
However, the non-zero kurtosis distribution can be clearly seen only in the large β cases when the lattice is highly nonlinear. In the low nonlinearity case, e.g. when β decreases to as low as 0.01, the kurtosis is hardly distinguishable from the background statistical fluctuations, see figure 1(d). This is because when the nonlinearity strength is lower than an equipartition threshold, the FPU-β lattice behaves like a harmonic one [15,20]; and for harmonic lattice it can be analytically shown that the single-particle velocity distribution remains Gaussian throughout the diffusion process, yielding an exactly zero kurtosis.

Velocity correlation
Since the kurtosis checking loses its effectiveness in the weak nonlinear cases, we now consider another criterion for LTE, the long-range correlation of particles' velocities. It is not surprising that due to the existence of some conserved quantities the long-range correlation in a system at a nonequilibrium state can be observed. However, if each small part of the system is in LTE then the correlation should be either vanished or short-ranged. So that locally thermal equilibrated regions can be well defined. Therefore, in the next step we check the spatial distribution of velocity correlation in the lattice.
In figure 2, we plot the heatmap of simultaneous two-particle velocity correlation v t v t i j á ñ ( ) ( ) for all the particle pairs in a lattice of length N=201. Since the number of unique pairs n=N(N+1)/2 increases rapidly with N, it is very challenging to further increase the size N with the limited computation power. á ñ ( ) ( ) emerge and move along the anti-diagonal direction. The most noticeable ones are the two dots in yellow and red. This implies that particles located symmetrically with respect to the center are also inclined to be correlated. It is because the changes of v i and v −i at the same time t are due to the fluctuation at the center of the lattice at the same early time.
In figure 3(a), we plot the velocity correlation v t v t i i á ñ -( ) ( ) at various time t. For t=0, the correlation displays a δ function shape, which confirms again that no long range velocity correlation exists. This indicates an LTE description is adoptable for the system at t=0 when the system is in contact with the heat baths and rests in a non-equilibrium steady state. For t>0 however, two side peaks appear and move out with a basically constant speed, the sound speed.
In figure 3(b) we plot the correlation at the same time t=120 for various values of β. The side peaks can be clearly observed for all cases, even for the β=0.01 case where the nonzero kurtosis is merely observable, see figure 1(d). Again it can be confirmed that such a long range velocity correlation is not a far-from-equilibrium phenomenon, either. In ( ) ( ) is linear with ΔT even in the small ΔT limit. For comparison, we also plot the correlation between each particle and the central particle, v t v t i 0 á ñ ( ) ( ) , in figure 3(d). We see that in all cases v t v t i 0 á ñ ( ) ( ) are localized around the center. It is interesting to see that although the central particle is not correlated with a far-away particle with large i | |, two farther-away particles can be correlated.
Is this long range correlation ubiquitous for the heat diffusion in 1D lattices? The answer is no. In the f 4 lattice displaying normal heat diffusion and obeying the Fourier's law, we observe a distinct yet expected behavior, as shown in figure 4. There we see that nonzero correlation only exists in a very narrow strip surrounding the diagonal line. Therefore, the velocity correlation is short-ranged and a locally well-defined thermal equilibrium state is not prohibited.
So far we have presented that in the heat diffusion process in the FPU-β lattices, LTE cannot be perfectly reached. Therefore, an local thermodynamic temperature cannot be well defined. If one still insists on defining the LT in terms of the kinetic temperature, i.e. equation (1), something paradoxical appears. In - º ---( )at various times, respectively. Not surprisingly, for any t>0, j i (t)>0 for i>0 and j i (t)<0 for i<0. Namely, heat always flows from the center towards the two ends. −∇T i should also be always positive and negative in the right and left half of the lattice, respectively. However, we see in figure 5(b) that this is not always true: at some places the local heat current j i flows against the opposite LT gradient −∇T i . Such a phenomenon is not allowed by the second law of thermodynamics.

Summary and discussion
We have investigated the validity of the LTE presumption in heat diffusion processes in a few 1D lattices. It was observed that the velocity distribution in a nonlinear lattice, e.g. the FPU-β lattice systematically departs from Gaussian. The kurtosis distribution displays a rapidly decaying central well as well as two slowly decaying side peaks spreading from the center. Interestingly, the speed of the peaks is even higher than the sound speed in the lattice. This indicates that the effect from the central heat package reaches a far-away particle before the heat reaches. By detecting this effect we would be able to sense in advance the upcoming heat wave which transports at the sound speed. Obvious practical value can be expected. Such a phenomenon does not vanish in the nearequilibrium limit because the rescaled value κ/Δt approaches a constant distribution in the small ΔT limit.
Besides kurtosis, velocity correlation between particles is another measure that we used to check the validity of LTE. For the FPU-β lattice, it was observed that particles at long distances can be correlated in the heat diffusion process. The profile of v v i i á ñ -, i.e. the velocity correlation between one particle and its symmetric particle as a function of i, looks alike to the profile of the energy distribution. Similar to the kurtosis, the rescaled correlation does not vanish but approaches a constant distribution in the small ΔT limit, which demonstrates that this is also a phenomenon observable for near equilibrium situations. We have also confirmed that the above findings are universal for a series of 1D nonlinear lattices. Some results are presented in the appendix.
Due to the breakdown of LTE, it would be natural to expect that the local version of Fourier's law j(x, t)=−κ(x, t) ∇ T(x, t), with T the kinetic temperature, is not applicable to the transient states in the anomalous diffusion process. This conjecture was confirmed by our intensive numerical simulations. It was observed that the instantaneous local heat current is not proportional to the instantaneous local kinetic temperature gradient; though the kinetic temperature profile still remains smooth. In some regions the directions of these two quantities are even opposite, i.e. the heat flows against the LT drop, which violates the local version of second thermodynamics law. This fact indicates that although a smooth kinetic temperature profile is commonly regarded as a good indicator for LTE in a stationary state [12], it loses its effectiveness in nonstationary anomalous heat diffusion process.
On the other side, when energy diffuses normally as in the f 4 lattice, all the abnormal phenomena studies above disappear and thus LTE is still a good approximation. Being the prototypes of anomalous and normal heat diffusion/conduction in 1D nonlinear lattices respectively, we expect that the above conclusions for the FPU-β lattice and the f 4 lattice to be generally valid in 1D lattices.