A Comparison of efficiencies of different methods in computation of heat conductance

# Heat conductance in nonlinear lattices at small temperature gradients

## Abstract

This paper proposes a new methodological framework within which the heat conductance in 1D lattices can be studied. The total process of heat conductance is separated into two parts where the first one is the equilibrium process at equal temperatures of both ends and the second one – non-equilibrium with the temperature of one end and zero temperature of the other. This approach allows significant decrease of computational time at . The threshold temperature is found which scales with the lattice size and by convention separates two mechanisms of heat conductance: phonon mechanism dominates at and the soliton contribution increases with temperature at . Solitons and breathers are directly visualized in numerical experiments. The problem of heat conductance in non-linear lattices in the limit can be reduced to the heat conductance of harmonic lattice with time-dependent stochastic rigidities determined by the equilibrium process at temperature . The detailed analysis is done for the -FPU lattice though main results are valid for one-dimensional lattices with arbitrary potentials.

Keywords : FPU lattice, nonlinearity, heat conductance, solitons

PACS numbers: 44.05.+e; 05.45.Yv

## I Introduction

The problem of heat conductance in low dimensional systems attracts much attention in last decades (see review (1)) and is motivated by the discovery of quasi-one-dimensional (nanotubes, nanowires, etc.) and two-dimensional (graphen, graphan, etc.) systems.

The modern theory of heat conductance was initiated by the celebrated preprint of E. Fermi, J. Pasta and S. Ulam (2), though the primary aim was “of establishing, experimentally, the rate of approaching to the equipartition of energy among the various degrees of freedom”. Subsequent investigations demonstrated wide area of consequences in many physical and mathematical phenomena (see reviews in special issues of journals CHAOS (3) and Lecture Notes in Physics (4) devoted to the 50th anniversary of the FPU preprint).

The dynamical properties of nonlinear systems in microcanonical ensemble (total energy const) were thoroughly analyzed in most papers. It allows to investigate the dynamics and to get exact results (soliton (5); (6); (7) and breather (8); (9); (10); (11); (12) solutions), to analyze regular and stochastic regimes and to find the corresponding thresholds. The FPU preprint also initiated the investigations in the field of “experimental mathematics” (13) .

About ten decades ago P. Debye argued that the nonlinearity can be responsible for the finite value of heat conductance in insulating materials (14). But modern analysis shows that it is not always the case. There are many examples where the coefficient of heat conductance diverges with the increasing of the system size as where , and in the thermodynamic limit (). Most of momentum conserving one-dimensional nonlinear lattices with various types of nearest-neighbor interactions have this unusual property (see, e.g., (1); (15); (16) ). Moreover, some other systems, – two- (17); (18); (19); (20) and three-dimensional lattices (21); (22), polyethylene chains (23), carbon nanotubes (24); (25); (26); (27); (28) have analogous property – diverging heat conductance with the increasing size of the system.

There were some conjectures explaining the anomalous heat conductance. Generally speaking, whenever the equilibrium dynamics of a lattice can be decomposed into that of independent “modes” or quasi-particles, the system is expected to behave as an ideal thermal conductor (29). Thereby, the existence of stable nonlinear excitations is expected to yield ballistic rather than diffusive transport. At low temperatures normal modes are phonons. At higher temperatures noninteracting “gas” of solitons starts to play more significant role and M. Toda was the first, who suggested the possibility of heat transport by solitons (30).

Though analytical expressions for solitons can be derived only for few continuum models described by partial differential equations, Friesecke and Pego in a series of recent papers (31); (32); (33); (34) made a detailed study of the existence and stability of solitary wave solutions on discrete lattices with the Hamiltonian , where . It has been proven that the systems with this Hamiltonian and with the following generic properties of nearest-neighbor interactions: has a family of solitary wave solutions which in the small amplitude, long-wavelength limit have a profile close to that of the KdV soliton. It was also shown (35) that these solutions are asymptotically stable. Thus most acceptable point of view on the origin of anomalous heat conductance in nonlinear lattices is as follows: phonons are responsible for heat conductance at low temperatures, and at high temperatures – solitons (36); (37).

A set of generic properties were found in a series of papers in investigation of dynamics of nonlinear lattices, starting from the celebrated preprint of FPU (2). And one is an existence of stochasticity thresholds. The weak stochasticity threshold is characterized by a specific energy below the which the trajectory in the phase space is almost regular (with near zero Lyapunov exponents) and only small part of normal modes is excited (it is just the case observed and analyzed by FPU). The strong stochasticity threshold corresponds to the value of above which energy equipartition between normal modes is established, and Lyapunov exponents are positive (38); (39); (40); (41); (42); (43); (16); (44); (45).

A major part of results was obtained using microcanonical ensemble for isolated systems. Physically more justified is the usage of canonical ensemble where temperature is kept (on average) constant by some or other type of heat baths. If the constant temperature is maintained by the Langevin sources (random forces with viscous friction), then from the Fokker-Planck equation the equilibrium Gibbs distribution immediately follows.

If one starts calculations from arbitrary initial conditions in canonical ensemble then some time is necessary to achieve the state of thermodynamic equilibrium. And this stage is not a trivial one (46). Firstly and surprisingly, kinetic and potential energy can relax to equilibrium with different rates; secondly, obeying the Maxwell velocity distribution function is not the sufficient condition of achievement the equilibrium. And the critical stage of achievement the equilibrium (energy equipartition between normal modes) is the excitation of the most longwave normal mode. Characteristic times of achievement the equilibrium can cover very wide range. For instance, there are well localized excitations in the harmonic lattice with random masses or, equivalently, random interparticle potentials (Anderson localization (47); (48)) where (49). And this phenomenon is explained by very weak interaction of localized excitations, if they are centered near the lattice center, with the heat reservoir located at the lattice ends.

If to return back to the problem of heat conductance, then one meets rather confusing experimental and numerical results, e.g. exponent in the dependence depends on the model under consideration, types of boundary conditions, used thermostat (Langevin or Nóse-Hoover (50); (51)), and also on temperature. For instance, temperature dependence of heat conductance in carbon nanotubes decreases as at K (52); experimentally is found (25) that also decreases with the growth of temperature. Different temperature dependencies vs. were found in 1D nonlinear lattices. For -FPU lattice: at and at (53) what is usually observed in insulating crystals. For the interparticle harmonic potentials and on-site potentials (e.g. Klein-Gordon chains) , i.e. heat conductance decreases with the growth of temperature (54). One more problem is the calculation of heat conductance at small temperature gradients. Usually these calculations are very time consuming because of great fluctuations of heat current and statistical averaging over large number of MD trajectories is necessary.

The paper organized as follows: in Section II we introduce new method for the calculation of the heat conductance which significantly decreases the computation time and diminishes the standard error. The method is based on the separation of the total process of heat conductance into two contributions: equilibrium and non-equilibrium and the latter one is responsible for the energy transfer. In the next Section we found that some quadratic mean values do not exhibit the expected tendency to reach zero values as the temperature difference . The threshold temperature , separating two regimes, – damped and undamped, is revealed. And the dependencies of on temperature and lattice length are found. Some modifications of the calculation of heat conductance in the limit are introduced in Section V. Direct evidences of the solitons contribution to the heat conductance are given in the next Section. -FPU lattice is considered as an example.

## Ii Heat conductance in the β-FPU lattice

We consider the one-dimensional lattice of oscillators with the interaction of nearest neighbors

 U=∑iu(yi),yi=xi−xi−1 (1)

and the -FPU potential (usually we put ).

Nonequilibrium conditions are necessary for the heat transport simulation. The most abundant method is the placement of the lattice into the heat bath with different temperatures of left and right ends (). Different types of heat reservoirs are thoroughly analyzed in (1). We utilize the Langevin forces acting on the left and right oscillators. are independent Wiener processes with zero mean and . is the temperature difference. The generalized Langevin dynamics with a memory kernel and colored noises is also suggested (55) to correctly account for the effect of the heat baths.

The following set of stochastic differential equations (SDEs)

 ¨xi=−∂U∂xi+δi1F++δiNF− (2)

are to be solved to find the heat flux . Then from the Fourier low the coefficient of heat conductance is

 ϰ=NJ/ΔT, (3)

and the problem is to find the heat current . The local heat flux (from th to th oscillator) is defined (56) by

 Ji→i+1=⟨Fi→i+1˙xi+1⟩;Fi→i+1≡−U′(xi+1−xi), (4)

where is a shorthand notation for the force exerted by the th on the th oscillator and is the time averaged. The total heat flux can be found as the mean value .

### ii.1 Equilibrium and non-equilibrium contributions to the heat conductance

If then the process of heat conductance can be formally separated into two parts: the first one – equilibrium process with equal temperatures of both lattice ends; and second – nonequilibrium process with temperature of the left lattice end and zero temperature of the right end (see Fig. 1) (by ‘process’ we hereafter assume for brevity the solution of the corresponding SDEs).

Namely the second process defines the heat transport realized against the background of the equilibrium process. Once we utilize this approach then the Langevin forces in (2) can be written as and for the left and right lattice ends, correspondingly; superscripts ‘0’ and ‘1’ refer to equilibrium and nonequilibrium processes. Then the total dynamical process can be represented as the sum of two processes

 x(t)=x0(t)+x1(t), (5)

where is the equilibrium (Gibbs’s) process at temperature , and – nonequilibrium, responsible for the energy transport, process. Then the Langevin dynamics is

 ¨x0i=−∂U0∂xi+δi1(ξ0−γ˙x01)+δiN(ξ0−γ˙x0N), (6)
 ¨x1i=−[∂U∂xi−∂U0∂xi]+δi1(ξ1−γ˙x11)+δiN(−γ˙x1N), (7)

and the sum of equations (6) and (7) is virtually identical to the parent equation (2). Random values and obey the identities and ; is the total energy (1) where the arguments of the total process are substituted to the coordinates of the equilibrium process . Expression in the square brackets in (7) is the difference of forces acting on the th particle from the total process and equilibrium process . It is significant to note that this force is the random value, and the process (heat transport) is realized in the lattice with time-dependent random potentials. The problem of heat conductance in the random time-independent potentials was analyzed in (57)

Equation (6) describes the system embedded in the heat reservoir at temperature . And is the stationary equilibrium process described by the canonical Gibbs distribution (equilibrium thermodynamics of the -FPU lattice in the canonical ensemble was considered in (46)).

Process is responsible for the heat transport and the Wiener’s process on the left lattice end defines small temperature . Right lattice end has zero temperature. An expression for the local heat flux is

 Ji→i+1=⟨Fi→i+1(x)˙xi+1−Fi→i+1(x0)˙x0i+1⟩, (8)

and the equilibrium process does not transfer energy: .

One of the goals of the present paper is the calculation of heat conductance at small temperature gradients. Usually these calculations are realized by solving SDEs (2) and are very time consuming because of great fluctuations of heat current (below we show that the time of computation increases if the accuracy of calculations is predetermined).

The comparison of two approaches (solving of standard SDEs (2) and (6)-(7)) is shown in Fig. 2 and results coincide with very good accuracy.

Note, that most of results in this paper are presented for the number of oscillators in the lattice. It may appear that this value is too small. For instance, the best estimate so far required simulations of up particles and integration steps plus ensemble averaging (1). But our results are aimed at founding some basic issues where number of particles is less essential. Lattices with larger number of oscillators were tested where necessary.

The dependence of heat conductance on the particles number is shown in Fig. 3 at two value of temperature . Inharmonicity becomes negligible in the limit and the analytical solution of the heat conductance for the harmonic lattice is given in (58).

There should be solved twice as large SDEs (6)-(7) in suggested approach as that in standard scheme (2), and this the price which is paid for the facility with using small temperature gradients. As one would expect, the accuracy of the suggested approach is higher (provided that all computational terms and conditions are identical). The comparison of accuracies is given in Appendix A

## Iii Strange behavior of process x1(t) at high temperatures

Usually the temperature difference at the lattice ends is an appropriate choice. Then the Fourier law (at fixed ) is valid with good accuracy. Actually, the corrections to the heat current are of the order as the current is the odd function of the temperature difference, and this ensures the reasonable accuracy of the linear approximation.

Now we concentrate our efforts on the elucidating the heat conductance dependence via temperature of the background process . Langevin forces , which provide temperature , are of the order (as ). And one can expect that process should have the same order because equation (7) becomes linear in the limit when . Thus any quadratic mean values should be of the order .

Two temperatures of the background process were tested: and (from here we omit subindex ‘–’ for brevity). Mean value was analyzed as an example and results are shown in Fig. 4 (fully identical properties have all quadratic values (correlators) of the types ). As one expects, the quadratic form linearly depends on : at . But the case is quite different at : mean value tends to a stationary value in the limit . It means that there exists some undamped stationary process at high temperatures even in the limit . These results also can imply an existence of a threshold temperature separating two regimes – damped at low temperatures and undamped at high temperatures.

## Iv The threshold temperature

Any process damps out at low temperatures and flattens out to a stationary value at higher temperatures even in the limit , and the temperature of process determines the different damping rates. And an illustrative process was analyzed to determine an existence of a threshold temperature and its value (’tilde’ marks the process at to avoid confusions).

Process can be exited in some or other manner. Usually and get random values in such a way that . The particular choice of initial conditions does not influence the final results.

Stochastic differential equation for the process are

 ¨˜x1i=−[∂U∂xi−∂U0∂xi]−δi1˙˜x11−δiN˙˜x1N;(γ=1), (9)

with random forces determined by the difference of processes and , and damping at the extreme left and right oscillators; and are potential energies with coordinates and , correspondingly. Stochastic dynamics (9) is implicitly ruled out by the temperature of process .

### iv.1 Two methods to find Tthr

We consider the case of small temperature when process is damped out. The damping is determined by the viscous friction of left and right oscillators in (9). Gradually increasing the temperature we find its threshold value when process becomes undamped.

The damping of mean squared displacement of the first oscillator was calculated. It was found that this process exponentially decays and depends on (see Fig. 5a). One can see that the damping stops in the range .

The dependence of coefficient on the temperature of process is shown in Fig. 5b and at .

Now we find the threshold temperature going “from up to down”, going from higher temperatures. At high temperatures there exists the stationary process outcoming from random forces (see (7) – expression in square brackets). Process decreases in the sense that all quadratic mean values tend to zero as temperatures approaches . When the threshold temperature reaches it threshold value, process disappears (see Fig. 6). The found threshold temperature is .

### iv.2 Time-resolved dynamics of process x1(t)

To elucidate the reasons of strange dynamics of process at high temperatures we analyzed it more thoroughly. As above, was calculated but without averaging over time. Results are shown in Fig. 7 for three temperatures of process .

One can see that behaves highly irregular. And numbers and heights of observed peaks increases with the growth of temperature until it becomes chaotic at high . The mean values at different temperatures averaged over time interval t.u. increase with temperature. Mean values , shown in horizontal solid lines, are nothing else than the stationary values calculated above. It was specially checked out that the dynamics observed in Fig. 7 is not due to numerical artifacts.

### iv.3 Heat conductance at small temperature gradients

Our main concern is the computation of heat conductance at small temperature gradients. With this in mind we analyze an expression for the heat current in more details. And this analysis can also shed some light upon the problem why process behaves in such strange manner. Remind an expression for force in the heat current: is the force acting on the th oscillator from left to right, and the derivative of the potential energy between oscillators is taken with respect to . Then the expression for the local heat current (8) can be rewritten as

 Ji→i+1=⟨[Fi→i+1(x)−Fi→i+1(x0)]˙x1i+1⟩=⟨[−(xi+1−xi)−(xi+1−xi)3+(x0i+1−x0i)+(x0i+1−x0i)3]˙x1i+1⟩, (10)

where – velocity of th oscillator in process . Difference of forces (in square brackets in the second line) is the polynomial of the third degree in the square root of temperature difference (process is of the order of as discussed above), and taking into account that the velocity is also of the order of , it is the polynomial of the forth degree in . But the coefficient of heat conductance is determined by the relation therefor terms of the third and forth orders can be neglected at small values of . Then (10) is simplified to

 Ji→i+1=⟨−(x1i+1−x1i)[1+3(x0i+1−x0i)2]˙x1i+1⟩. (11)

It is significant that the total heat current in (10) at large temperature and small temperature gradient is the difference of finite terms from processes and , vanishing in the limit . And it is the reason why direct MD simulation is highly inefficient in this case and gives very large fluctuations. But, as will be shown below, there exists an efficient method to overcome this difficulty.

The behavior of process is explained by the fact that it is determined not only by random Langevin forces , but also (and more essentially) by time-dependent random forces (see (7)). The plateau for the correlator equal to at is determined exclusively by random forces from stationary process (an illustrative example of one variable is considered in Appendix B). Thus, the dynamical process becomes the stationary one, determined by the background process at high temperatures.

## V Threshold temperature in the limit ΔT→0

In this section process is considered at an arbitrary temperature and in the limit . Remind that process is completely suppressed at . To realize the limiting transition in (7) it is convenient to divide both sides by . Then new coordinates are . It is also convenient to introduce normalized to unity random force . Then the linear equation for can be obtained as quadratic and cubic forces can be neglected (see (10)). The corresponding equation for can be derived if to rearrange one term in potential energy (1) keeping in mind that . Then .

Transforming variables to and retaining terms quadratic in , one can get the potential energy in the form

 u=12∑igi(t)(yi−yi−1)2,gi(t)=1+3[x0i(t)−x0i−1(t)]2, (12)

where are time-dependent random coefficients of rigidity determined by the dynamical process . It is illuminating to note that the problem of heat conductance can be reduced to the quadratic potential energy in the limit . Corresponding SDEs have Langevin source with unit temperature at the left oscillator and zero temperature at the right oscillator:

 ¨yi=−gi(yi−yi−1)+gi+1(yi+1−yi)+δi1(θ−˙y1)−δiN˙yN. (13)

It should be also noted that if the 1D lattice with an arbitrary interaction (Morse, Toda, LJ, etc) is analyzed then the corresponding equation will be the same, and random rigidities are ) where is some or other type of potential energy. Equation for arbitrary systems (with arbitrary neighbor radius of interaction) can be also written in the general form as

 ¨yi=−M∑j=1Λ0ijyj+δi1(θ−˙y1)−δiN˙yN, (14)

where – matrix of second derivatives of potential energy depending on , and is the number of neighbors. Equation (14) is valid for arbitrary systems.

Equations (14) define the stationary random process only if temperature . As temperature approaches the value , quadratic mean values diverge. It is shown in Fig. 8. And this is the third method to find .

There was analyzed the case of low temperatures when process is “weak”. Then rigidity coefficients are close to unity (see (12)). And as an example we consider the lattice where actual rigidity coefficients (12) are substituted by the mean value taken from the equilibrium Gibbs distribution and . This harmonic model is exactly solvable and results are shown in Fig. 9

One can see that process is damped out in the model with constant rigidity in contrast to the case when actual values (12) are used. And one can conclude that the growth of process , when temperature increases, is determined by an increase of fluctuations but not only by the increase of rigidities.

Process diverges at high temperatures. And it gives one more possibility, the fourth one, to find the threshold temperature. To attain this end the equilibrium process at temperature is established. Then process is excited in some or other way (its initial conditions do not influence the final results). And the evolution of the is analyzed. One can see (Fig. 10) that the process exponentially damps out at and exponentially grows at .

It should be stressed out that the method just described differs from the previous one (Fig. 6 and discussion). The nonlinear case was considered there and its stationarity was conditioned by nonlinear terms in forces which are absent in the harmonic approximation.

Four methods give the threshold temperature . This temperature was found for a fixed lattice length . Larger lattice lengths were considered and the dependence of on the lattice length is shown in Fig. 11. Approximate fitting gives dependence .

It means that the majority of usually studied lattices are in the state when their temperatures are much higher then the threshold temperature (e.g. if then ).

## Vi Sound velocity and solitons in β-FPU lattice

Heat conductance is observed in both regimes, – higher and below the threshold temperature. And an attempt was undertaken to find the soliton contribution to the heat conductance at . With this in mind, the correlator was analyzed ( is the displacement of th oscillator from equilibrium at time instant ). In numerical simulations we fixed the time shift t.u. and calculated the corresponding correlator (). Results are shown in Fig. 12.

The correlator has peaks at the coordinate shifts . It allows to calculate the velocity of excitation propagation and . This velocity is higher then the sound velocity calculated in the harmonic approximation at .

Initially these peaks were attributed to solitons. But more thorough analysis shows that this concepts is not valid. Let we have the -FPU potential ( and ). In (59) it was shown that there exists a spectrum of frequencies which are proportional to the harmonic ones, according to a well defined law. Therefor the -FPU potential can be represented as

 u(y)=(1+12y2)12y2 (15)

and an expression in brackets can be replaced by an effective harmonic rigidity

 u(y)=keff12y2 (16)

The problem is to find . It can be done in terms of a mean field approximation (MFA). Mean value of potential energy is

 ⟨up(y)⟩=keff12⟨y2⟩, (17)

where is the mean value of .

In the harmonic approximation (at not too high temperatures) mean values of potential and kinetic energies are equal . In canonical ensemble the identity is valid for 1D systems. Then

 keff12⟨y2⟩=T2 (18)

The self consistency of the MFA is (expression in brackets in (15) = )

 (1+12⟨y2⟩)=keff (19)

From (18) it follows that and substitution of into (19) gives the self-consistent equation for

 1+T/(2keff)=keff (20)

with the solution

 keff=12+√14+T2 (21)

Thereby the the effective (“nonlinear”) sound velocity

 veff=√keff=√12+√14+T2,  (m=1) (22)

is higher then the velocity found in the harmonic approximation. Eq. 22 gives for what coincides with the value found from correlation functions. The dependence of effective sound velocity versus temperature is shown in Fig. 13.

Note that the MFA is valid up to very high temperatures , while this approach originally is well suited only for low temperatures, and the effective sound velocity exceeds its harmonic value (at ) by .

Next we try to find direct evidences on the solitons participation in energy transfer. It was done in the following manner. Initially lattice of oscillators was thermalized for some time to reach the thermodynamic equilibrium. Then the “cold” lattice (with zero velocities and displacements) with 1000 oscillators was switched to the right end of the lattice. Solitons, if they exist in the initial lattice, should “run out” to the cold lattice. The same is valid for the moving breathers. (Note that in the continuum approximation the mKdV equation corresponds to the discrete -FPU potential. And one can find analytical expressions for solitons of compression, antisolitons of elogation and different types of breathers in (60) ). We waited some time till excitations run out of the lattice to its cold part where they can be observed. Results are shown in Fig. 14

Analogous approach to visualize breathers in 2D lattice with three different on-site potentials was utilized in (61) where initially thermalized lattice was cooled from the borders and breathers were detected after thermal noise was deleted through damping boundaries.

The possibility of energy transfer due to solitons was conjectured three decades ago (30). Less studied is the possibility of energy transfer by breathers. One suggested mechanism is the Targeted Energy Transfer (62); (63) when an efficient energy transfer can occur under a precise condition of nonlinear resonance between discrete breathers. Various aspects and possible applications of energy transfer by breathers are considered in (64).

## Vii conclusions

In conclusion we briefly summarize our main results. A new method is developed which allows considerable decreasing of the computation time in calculations of the heat conductance at low temperature gradients (temperature of the left lattice end and – of the right end and ). This success was achieved by the separation of the total process of heat conductance into two parts: an equilibrium process at equal temperatures of both lattice ends and non-equilibrium process , responsible for the energy transport, which occurs at temperature of one end and zero temperature of other end. The equilibrium (background) process strongly influences the transport properties: there exists the threshold temperature above which some undamped characteristics are observed; more precisely, correlators of the types do not tend to zero at , as expected, but have certain nonzero values; at “normal” dependence is observed, i.e. these correlators have zero values when . The reason of two distinct behaviors is not due to the temperature of the background process but sooner to the temperature fluctuations. An illustrative example of one variable is briefly analyzed where the threshold temperature is also found. The model of one variable has a rich family of solutions depending on the parameters and serves to be investigated in more thoroughly.

The threshold temperature was found by few methods and scales with the lattice size . All practically interesting systems lies above . The threshold temperature is not sharply pronounced and arbitrarily separates two mechanisms of the heat conduction: the phonon mechanism prevails at , and at the soliton contribution starts to play more significant role with the increase of temperature. Highly probable that the temperature fluctuations are responsible for the solitons generation.

Analytical solutions for solitons and breathers are known for the -FPU lattice and these excitations were directly observed in numerical experiments. Our findings on the soliton contribution to the heat conductance are in accordance with the general scenario of heat conductance: phonons gave main contribution to the heat conductance at low temperatures and solitons more and more dominate when temperature increases.

We found no relations between the well known weak and strong stochasticity thresholds and the threshold temperature in the present paper: they have different energy ranges and different dependencies on . Additional difference is due to different statistical ensembles used: traditionally stochasticity thresholds are found in microcanonical ensemble, but critical temperature is observed in canonical ensemble where energy equipartition is realized at any temperature. Statistical properties do coincide in the thermodynamical limit for -canonical and canonical ensembles, but the dynamical properties can differ.

## Appendix A Comparison of efficiencies of different methods in computation of heat conductance

The main problem in the calculation of heat conductance is to find the heat current with an appropriate standard error. And we consider below the effectiveness of the separation of the total process into the sum of two: in the sense of computer time expenditure (see (2) and (6)-(7) ), where is the equilibrium background process and process is responsible for the heat transport.

The comparative efficiency of two approaches (‘old’ – solving of SDEs (2) and ‘new’ – two SDEs (6)-(7) ) to the calculation of the heat flux can be estimated as the relation of their standard errors at equal conditions of computation: . The result is shown in Fig. 15. One can see that the standard error is systematically less in the suggested approach as compared to the usually utilized.

More impressive is the behavior of efficiency at decreasing of the temperature gradient (see Fig. 16) and at equal parameters of numerical simulations. Here it should be emphasized that standard errors increase with the diminishing of in wide range and the efficiency Eff remains as before. The standard error increases as .

## Appendix B An example of one variable

We consider an example of an equation with one variable for the harmonic oscillator with damping

 ¨x=−k(t)x−γ˙x. (23)

where is the stochastic rigidity. This equation is the illustrative analogue of multi-variable SDEs equations (7) for the process . The variable substitution eliminates the damping and (23) can be reduced to

 ¨X=−k(t)X. (24)

Potential energy (12) has the form in the case of one variable, where and – stochastic process generated by the background process . And it is reasonable to choose the random rigidity in (24) in the form

 k(t)=1+εz2(t), (25)

where is free parameter and is the stationary random process describing the dynamic of the harmonic oscillator influenced by Langevin source with temperature :

 ¨z=−z+ξ−γ˙z (26)

and

We compare the solution (23)–(25) with the solution of well known deterministic Mathieu equation

 ¨y=−[1+gcos2(t)]y. (27)

Different types of solutions of the Mathieu equation depend on the parameter and initial conditions. There exists such that the solution is the sum of periodic functions at , and the solution is the superposition of periodic functions multiplied by the exponentially increasing and decreasing function at .

Equations (24)–(26) also have a rich family of solutions depending on initial conditions and parameter values. As an illustrative example we consider the following set of parameters: .

Below we demonstrate only the qualitative behavior of process depending on the temperature of stochastic process (26) and results are shown in Figs. 17(a–c). One can see different regimes as increases. And there exists some critical temperature above which the process diverges . It should be noted that the overall scenario strongly depends on the choice of initial conditions () and the particular sequence of random Langevin forces in (26). At larger times process becomes more complex. The full analysis of system (24)–(26) is not our primary goal, but these equations serve more intensive attention.

### References

1. S. Lepri, R. Livi, A. Politi, Phys. Reports 377, 1 (2003).
2. E. Fermi, J. Pasta, S. Ulam, Document LA–1940 (May 1955); in Collected Papers of E. Fermi (University of Chicago Press, Chicago, 1965), Vol. 2, p. 78.
3. CHAOS 15(1) (2005) (Focus Issue: The Fermi-Pasta-Ulam Problem - The First 50 Years, Ed. by D.K. Campbell, Ph. Rosenau, and G.M. Zaslavsky).
4. G. Gallavotti (Ed.), The Fermi-Pasta-Ulam Problem: A Status Report, Lect. Notes Phys. 728 (Springer, Berlin Heidelberg 2008).
5. M.D. Kruskal, N.J. Zabusky, J. Math. Phys. 5, 231 (1964)
6. N.J. Zabusky, M.D. Kruskal, Phys. Rev. Lett. 15, 240 (1965).
7. R.K. Dodd, J.C. Eilbeck, J.D. Gibbon, H.C. Morris, Solitons and Nonlinear Wave Equations (Academic Press, New York, 1982).
8. D.K. Campbell, M. Peyrard, Physica D 18, 47 (1986).
9. A. J. Sievers, S. Takeno, Phys. Rev. Lett. 61, 970 (1988).
10. T. Dauxois, M. Peyrard, Phys. Rev. Lett. 70, 3935 (1993).
11. S. Aubry, Physica D 71, 196 (1994).
12. R.S. MacKay, S. Aubry, Nonlinearity 7, 1623 (1994).
13. M.A. Porter, N.J. Zabusky, B. Hu, D.K. Campbell, American Scientist 97(3), 214, (2009).
14. P. Debye, Vorträge über die Kinetische Theorie der Wärme (Teubner, 1914).
15. G. Casati, B. Li, arXiv:cond-mat/0502546.
16. A.J. Lichtenberg, R. Livi, M. Pettini, S. Ruffo, Lect. Notes Phys. 728, 21 (2008).
17. X.Yi, J.A.D. Wattis, H. Susanto, L.J. Cummings, J. Phys. A 42, 355207 (2009).
18. I.A. Butt, J.A.D. Wattis, J. Phys. A 40, 1239 (2006).
19. I.A. Butt, J.A.D. Wattis, J. Phys. A 39, 4955 (2006).
20. S. Flach, K. Kladko, C.R. Willis, Phys. Rev. 50, 2293 (1994).
21. S. Flach, K. Klado, R.S. MacKay, Phys. Rev. Lett. 78, 1207 (1997).
22. H. Shiba, N. Ito, J. Phys. Soc. Jap. 77, 054006 (2008).
23. A. Henry, G. Chen, Phys. Rev. B 79, 144305 (2009).
24. Z. Yao, J.-S. Wang, B. Li, G.-R. Liu, Phys. Rev B 71, 085417 (2005).
25. C. Yu, L. Shi, Z. Yao, D. Li, A. Majumdar, Nano Lett. 5, 1842 (2005).
26. S. Maruyama, Physica B 323, 193 (2002).
27. N. Mingo, and D.A. Broido, Nano Lett. 5, 1221 (2005)
28. J.X. Cao, X.H. Yan, Y. Xiao, J.W. Ding, Phys. Rev. B 69, 073407 (2004).
29. S. Lepri, R. Livi, A. Politi, CHAOS 15, 015118 (2005).
30. M. Toda, Phys. Scr. 20, 424 (1979).
31. G. Friesecke, R.L. Pego, Nonlinearity 12, 1601 (1999).
32. G. Friesecke, R.L. Pego, Nonlinearity 15, 1343 (2002).
33. G. Friesecke, R.L. Pego. Nonlinearity 17, 207 (2004).
34. G. Friesecke, R.L. Pego. Nonlinearity, 17, 229 (2004).
35. A. Hoffman, C.E. Wayne, J. Gynamics and Diff. Equations 21, 343 (2009).
36. B. Li, J. Wang, L. Wang G. Zhang, CHAOS 15, 015121 (2005).
37. H.J. Viljoen, L.L. Lauderback D. Sornette, Phys. Rev. E 65, 026609 (2002).
38. P. Bocchieri, A. Scotti, B. Bearzi, A. Loinger, Phys. Rev. A 2, 2013 (1970).
39. M. Pettini, M. Landolfi, Phys. Rev. A 41, 768 (1990).
40. L. Casetti, R. Livi, M. Pettini, Phys. Rev. Lett. 74, 375 (1995).
41. L. Casetti, M. Cerruti-Sola, M. Pettini, E.G.D. Cohen, Phys. Rev. E 55, 6566 (1997).
42. L. Galgani and A. Giorgilli, J. Math. Sci. 128, 2761 (2005).
43. G. Gallavotti, Lect. Notes Phys. 728, 1 (2008).
44. M. Pettini, L. Casetti, M. Cerruti-Sola, R. Franzosi, E.G.D. Cohen, CHAOS 15, 015106 (2005).
45. G.M. Zaslavsky, CHAOS 15, 015103 (2005).
46. V.N. Likhachev, T.Yu. Astakhova, W. Ebeling, M.G. Velarde, and G.A. Vinogradov, Eur. Phys. J. B 72, 247 (2009).
47. P.W. Anderson, Phys.Rev. 109, 1492 (1958).
48. P.W. Anderson, Phys.Rev. 124, 41 (1961).
49. V.N. Likhachev, G.A. Vinogradov, T.Yu. Astakhova, A.E. Yakovenko, Phys. Rev. E 72, 016701 (2005).
50. S. Nóse, J. Chem. Phys. 81, 511 (1984).
51. W.G. Hoover, Phys. Rev. A 31, 1695 (1985).
52. S. Maruyama, Microscale Thermophysical Engineering, 7, 41 (2003).
53. K. Aoki, D. Kusnezov, Phys. Rev. Lett. 86, 4029 (2001).
54. K. Aoki, D. Kusnezov, Phys. Lett. A 265, 250 (2000).
55. J.-S. Wang, Phys. Rev. Lett. 99, 160601 (2007).
56. Y. Zhang, H. Zhao, Phys. Rev.E 66, 026106 (2002).
57. M. Johansson, G. Kopidakis, S. Lepri, S. Aubry, Europh. Lett. 86, 10009 (2009).
58. Z. Rieder, J.L. Lebowitz, E.Lieb, J. Math. Phys. 8, 1073 (1967).
59. C. Alabiso, M. Casartelli, J. Phys. A 34, 1223 (2001).
60. G. Lamb Elements of soliton theory (John Willey & Sons, New York Chichester Brisbane Toronto 1980).
61. A. Bikaki, N.K. Voulgarakis, S. Aubry, G.P. Tsironis, Phys. Rev. E 59, 1234 (1999).
62. G. Kopidakis, S. Aubry, G.P. Tsironis, Phys. Rev. Lett 87, 165501 (2001).
63. P. Maniadis, G. Kopidakis, S. Aubry, Physica D 188, 153 (2004).
64. S. Aubry, Physica D 216, 1 (2006).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters