Heat conductivity in the \beta-FPU lattice. Solitons and breathers as energy carriers

Heat conductivity in the -FPU lattice.
Solitons and breathers as energy carriers

T.Yu. Astakhova    V.N. Likhachev    G.A. Vinogradov Emanuel Institute of Biochemical Physics RAS,
ul. Kosygina 4, Moscow 119334, Russian Federation
Abstract

This paper consists of two parts. The first part proposes a new methodological framework within which the heat conductivity in 1D lattices can be studied. The total process of heat conductivity is decomposed into two contributions where the first one is the equilibrium process at equal temperatures of both lattice ends and the second – non-equilibrium process with the temperature of one end and zero temperature of the other. This approach allows to isolate and analyze the heat transfer in explicit form. The heat conductivity in the limit is reduced to the heat conductivity of harmonic lattice with stochastic rigidities determined by the equilibrium process at temperature . A threshold temperature is found which separates two regimes: small perturbations exponentially decay at and tend to constant value at . The threshold temperature scales with the lattice size and in the thermodynamic limit. Some unusual properties of heat conductivity can be exhibited on nanoscales at low temperatures. The thermodynamics of the -FPU lattice can be adequately approximated by the harmonic lattice with temperature renormalized coefficients of rigidity. The second part testifies in the favor of the soliton and breather contribution to the heat conductivity in contrast to conclusions made in [N. Li, B. Li, S. Flach, PRL 105 (2010) 054102]. In the long-wavelength continuum limit the discrete -FPU lattice is reduced to the modified Korteweg – de Vries equation. This equation has soliton and breather solutions. Numerical simulations demonstrate their high stability. New method for the visualization of moving solitons and breathers is suggested. An accurate expression for the dependence of the sound velocity on temperature is also obtained. Our results support the conjecture on the solitons and breathers contribution to the heat conductivity. The fraction of total heat flux transferred by solitons and breathers merits additional analysis.

Keywords: heat conductivity, -FPU lattice, soliton, breather

The corresponding author: gvin@deom.chph.ras.ru.

I Introduction

The problem of heat conductivity in low dimensional systems attracts much attention in last decades (see review Lep03 ()) 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 conductivity was initiated by the celebrated preprint of E. Fermi, J. Pasta, and S. Ulam Fer55 (), 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 Cha05 () and Lecture Notes in Physics Lec08 () devoted to the 50th anniversary of the FPU preprint).

The dynamical properties of nonlinear systems in microcanonical ensemble (total energy const) were thoroughly analyzed. It allows to investigate the dynamics and to get exact results (soliton Kru64 (); Zab65 (); Dod82 () and breather Cam96 (); Sie88 (); Dau93 (); Aub94 (); Mac94 () 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” Por09 () .

About ten decades ago P. Debye argued that the nonlinearity can be responsible for the finite value of heat conductivity in insulating materials Deb14 (). But modern analysis shows that it is not always the case. There are many examples where the coefficient of heat conductivity diverges with the increasing 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., Lep03 (); Cas05 (); Lic08 () ). Moreover, some other systems, – two- Bar07 (); Lip00 (); Sai10 () and three-dimensional lattices Shi08 (), polyethylene chain Hen09 (), carbon nanotubes Yao05 (); Yu_05 (); Mar02 (); Min05 (); Cao04 () have analogous property – diverging heat conductivity with the increasing size of the system.

There were some conjectures explaining the anomalous heat conductivity. 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 Lep05 (). 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 or/and breathers starts to play more significant role, and M. Toda was the first who suggested the possibility of heat transport by solitons Tod79 ().

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 Fri99 (); Fri02 (); Fri04a (); Fri04b () 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 Hof08 () that these solutions are asymptotically stable. Thus most acceptable point of view on the origin of anomalous heat conductivity in nonlinear lattices is as follows: phonons are responsible for heat conductivity at low temperatures, and at high temperatures – solitons Li_05 (); Vil02 ().

Rather confusing experimental and numerical results are demonstrated in literature about the dependence of heat conductivity on different parameters, – model under consideration, types of boundary conditions, used thermostat and temperature. For instance, temperature dependence of heat conductivity in carbon nanotubes decreases as at K Mar03 (); experimentally is found Yu_05 () that also decreases with the growth of temperature. Different temperature dependencies vs. were found in 1D nonlinear lattices. For -FPU lattice: at and at Aok01 () what is usually observed in insulating crystals. For the interparticle harmonic potentials and on-site potentials (e.g. Klein-Gordon chains) , i.e. heat conductivity decreases with the growth of temperature Aok00 (). But there exists firm theoretical background Nar02 () that the exponent in the dependence is the universal constant in momentum-conserving systems.

The calculation of heat conductivity at small temperature gradients is an additional problem. 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 the heat conductivity is considered when the temperature gradient is small. The explicit contribution to the heat conductivity is extracted by the decomposing of the total process into two parts. The first one is the equilibrium process at temperatures of both lattice ends, and the second – non-equilibrium, when one lattice end has temperature and the other – zero temperature. Namely the latter process is responsible for the heat conductivity. This method allows to find the threshold temperature . And though scales with the lattice length as , unusual dynamics can be revealed on nanoscale when both and are small. The low temperature thermodynamics can be adequately described in terms of harmonic lattice with temperature renormalized rigidity coefficients.

Section III is independent of the previous one and is aimed at elucidating the role of solitons and breathers in the heat conductivity. Solitons and breathers are found as the solutions of the modified Korteweg – de Vries equation. This equation is obtained in the continuum limit from the discrete -FPU lattice. Soliton and breather solution are checked in numerical simulation and demonstrate very high stability.

Necessary details of the derivation of accurate soliton and breather solutions in the continuum limit are given in Appendix.

Ii Heat conductivity in the -FPU lattice

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

(1)

(usually the dimensionless potential will be used below, that is ). 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 Lep03 (). The usage of the Langevin forces with the noise terms and friction forces acting on the left and right oscillators is the common practice ( is also put for brevity). are independent Wiener processes with zero mean and the correlator . is the temperature difference. The generalized Langevin dynamics with a memory kernel and colored noises is also suggested Wan07 () to correctly account for the effect of the heat baths.

The following set of stochastic differential equations (SDEs)

(2)

is usually solved to find the heat flux . And the local heat flux (power transmitted from th to th oscillator) is Zha02 ()

(3)

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

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

If then the process of heat conductivity can be formally decomposed into two contributions: the first one – equilibrium process with equal temperatures of both lattice ends; and the 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).

Figure 1: Schematic representation of the total process as sum of equilibrium and non-equilibrium processes.

Namely the second process is responsible for the heat transport taking place on(?) the background of equilibrium process. Once this approach is utilized then the noise terms in (2), owing to their independence, are and for the left and right lattice ends, correspondingly Superscripts ‘0’ and ‘1’ refer to equilibrium and nonequilibrium processes. The total process can be represented as the sum

(4)

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

(5)
(6)

and the sum of equations (5) and (6) is 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 replaced by coordinates of the equilibrium process . Expression in the square brackets in (6) is the difference of forces acting on the th oscillator from the total process and equilibrium process . It is worth mentioning 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 conductivity in the random time-independent potentials was analyzed in Joh08 ().

Equation (5) describes the system embedded in the heat reservoir at temperature of both lattice ends. And is the stationary equilibrium process described by the canonical Gibbs distribution. Process is responsible for the heat transport and the Wiener process on the left lattice end defines temperature . Right lattice end has zero temperature (only the friction force acts on this oscillator). An expression for the local heat flux is

(7)

and the equilibrium process does not transfer energy: , where stands for the time average. It is essential that the heat flux (7) is the small difference of large values from processes and . This is the reason why MD simulation gives large fluctuation when and is not low. It can be shown that the time of computation increases if the standard error is fixed. The comparison of two approaches (solving of standard SDEs (2) and (5)-(6)) is shown in Fig. 2 and results coincide with very good accuracy.

Figure 2: Temperature dependence of heat conductivity for the lattice of oscillators. Filled circles: solution of standard SDEs (2); empty circles: SDEs (5)-(6). Averaging over MD trajectories time units (t.u.) each. , . Triangle up at is the exact value in the harmonic approximation ().

Some results in this paper are obtained for the number of oscillators . It may appear that this value is too small. For instance, “standard” simulations require up to particles and integration steps plus ensemble averaging Lep03 (). But our results are aimed at establishing some new issues where number of particles is unessential. Lattices with larger number of oscillators were tested when necessary.

Relative displacements and for processes and are shown in Fig. 3. These values characterize the energy fluxes in the lattice. Energy fluxes to the left and to the right are equal on average for process as it is the equilibrium process without energy transfer. But the energy flux is directed mainly to the right for process (right panel). Low temperature is chosen for the better illustration.

Figure 3: Spatiotemporal evolutions of relative displacements and for the equilibrium (left panel) and nonequilibrium (right panel) processes. , , .

The dependence of heat conductivity on the oscillators number is shown in Fig. 4 at two value of temperature . Results coincide with very good accuracy. Inharmonicity becomes negligible at low temperature and heat conductivity at (cirles in Fig. 4) coincides with the heat conductivity of the harmonic lattice (dashed line) with good accuracy. The analytical solution of the heat conductivity for the harmonic lattice is given in Rie67 ().

Figure 4: Coefficient of heat conductivity for the -FPU lattice for oscillators. Squares: , circles: . Filled symbols – results obtained by the solution of standard SDEs (2), empty symbols – SDEs (5)-(6). Averaging over 200 MD trajectories t.u. . Dashed line – harmonic approximation. Filled symbols are practically fully covered by empty symbols and are invisible.

There should be solved twice as large SDEs (5)-(6) in suggested approach as that in standard scheme (2). But this approach has some undoubted merits as discussed below.

ii.2 Heat conductivity at small temperature gradients

One of the goals of the present paper is the computation of heat conductivity at very small temperature gradients. With this in mind an expression for the heat flux is analyzed in more details. The expression for the local heat flux (7) can be rewritten as

(8)

where – velocity of th oscillator in process . Taking in mind that , the expression (8) can be transformed to

(9)

Processes and have different ranges of specific energies. Noise terms , which provide temperature , are of the order (as ). And one can expect that process has the same order because equations (6) become linear in the limit when . Expression in curly brackets in (9) is the polynomial of the third degree in the square root of temperature difference . Taking into account that the velocity is also of the order , expression (9) is the polynomial of the forth degree in . But the coefficient of heat conductivity is determined by the relation . Then terms of the third and forth orders can be neglected at . Then (9) is simplified to

(10)

An expression for the potential energy, corresponding to process , can be derived analogously. This energy is the difference of potential energies and again, using coordinates and , and preserving only terms quadratic in , one can get the potential energy for process in the form

(11)

where are time-dependent random coefficients of rigidities determined by process .

It is illuminating to note that in the limit the problem of heat conductivity in the -FPU lattice is reduced to the harmonic lattice with random coefficients. Corresponding SDEs have noise terms with friction forces on the left oscillator and zero temperature (only viscous forces) on the right oscillator:

(12)

and are defined in (11). If the 1D lattice with an arbitrary interaction potential is analyzed then the corresponding equations are the same with rigidities ) where is the potential energy. SDE for the system with arbitrary neighbor radius of interaction can be written in the general form as

(13)

where – matrix of second derivatives of potential energy depending on , and is the number of neighbors.

ii.3 Unusual dynamics of process at high temperatures

The process can be characterized by some time average correlators. The correlator was analyzed for the better understanding of the dynamics of process . One can expect that as discussed above. Two temperatures of the background process were tested: and (hereafter subindex ‘–’ is omitted, that is ). Results are shown in Fig. 5.

As expected, the correlator linearly depends on : at . But the case is quite different at : reaches the stationary value in the limit . It means that there exists some undamped stationary process at high temperatures of the background process even in the limit . These results also imply an existence of a threshold temperature separating two regimes – damped at low temperatures and undamped at high temperatures.

Figure 5: Dependence of correlator on the temperature difference . Filled circles: , empty circles: . Asymptotic value at (coefficient of linear regression ). Averaging over 100 MD trajectories t.u. each. . The range of : .

ii.4 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 defines different damping rates. Bearing this in mind, it is convenient to excite some auxiliary process over the background process and to analyze it.

Coordinates and velocities of process get random increments . The particular choice of initial conditions does not influence the final results. The total dynamics is the sum of two processes .

Stochastic differential equations for the process are

(14)

and only viscous forces acts at the extreme left and right oscillators. and are potential energies with coordinates and , correspondingly. Stochastic dynamics (14) is implicitly ruled out by the temperature of process .

To find the threshold temperature we initially consider the case of small temperature when process is damped out. Its damping is determined by the viscous friction of left and right oscillators in (14). Gradually increasing the temperature its threshold value can be found when process becomes undamped. The damping of mean squared displacement of the first oscillator was calculated. Process exponentially decays and depends on (see Fig. 6a).

Figure 6: a) Exponential damping of process at different temperatures: (circles), (squares), (triangles up), (triangles down). Solid lines – linear regressions. Averaging time t.u. 20 trajectories were used to estimate the standard error. b) Damping coefficient as the function of temperature of process . Damping stops at . .

The temperature dependence of coefficient is shown in Fig. 6b and when .

Next method to find the threshold temperature is moving “from up to down”, going from higher to lower temperatures. At high temperatures there exists the stationary process arising from random forces (see (14)). And process decreases in a sense that all quadratic mean values tend to zero as temperatures approaches from above. When the temperature reaches its threshold value, process is totally damped (see Fig. 7). The found threshold temperature is .

Figure 7: Stationary values at . Time of averaging t.u. The temperature dependence is approximated by the function (solid line). .

Strange behavior of process is basically explained by time-dependent random forces (see (6)) rather then random Langevin forces . And the plateau for the correlator at is determined exclusively by the background process .

This conjecture can be additionally supported. Let us consider the case of low temperatures when process is “weak”. Then the rigidity coefficients in (11) are close to unity. The lattice where actual rigidity coefficients (11) are substituted by the mean values taken from the equilibrium Gibbs distribution: and is considered as an example. This harmonic model is exactly solvable and results are shown in Fig. 8

Figure 8: Temperature dependence of the mean squared displacement . Circles – MD simulation of SDEs (12); dashed line – model of mean rigidities in the harmonic approximation. Averaging over 20 trajectories t.u. each.

One can see that process damps out in the model with constant rigidity in contrast to the case when actual values (11) are used. And the growth of process , when temperature increases, is mainly governed by fluctuations rather then the increasing rigidities.

Additional evidence of threshold phenomena is the one-dimensional analogue of the Mathieu equation

(15)

Different types of solutions depend on the parameter and initial conditions. There exists such critical value that the solution is the superposition of periodic functions at , and the solution diverges at .

We consider an equation for the harmonic oscillator for one variable with friction force

(16)

where – stochastic rigidity. This equation is similar to equation (12) for process if and is the stochastic process generated by the dynamics of harmonic oscillator with noise term and friction force at temperature :

(17)

and spectral property . The substitution excludes damping and (16) becomes

(18)

what is similar to the Mathieu’s equation (15).

Eqs. (17)-(18) have rich family of solutions depending on initial conditions, temperature and the sequence . The solution is nearly harmonic function at very low temperatures. The superposition of harmonic functions is the solution at higher temperatures. At last there exists such threshold temperature that the solution diverges and is the product of harmonic functions by . The solution is the product of some stochastic process by at much higher temperatures, and . There are many interesting intermediate solutions, and this problem merits more attention. Considered examples show that an existence of threshold phenomena is not exceptional and can occur in different dynamical systems.

The threshold temperature was found for the lattice length . Larger lattice lengths were considered and the dependence of on the lattice length is shown in Fig. 9. Fitting gives dependence .

Figure 9: Dependence of vs. lattice length in log-log coordinates. Solid line is the fitting .

ii.5 Time-resolved dynamics of process

Dynamics of process at high temperatures was analyzed above in terms of time-average correlators. And time-resolved behavior of process at two temperatures of the background process is shown in Fig. 10. As above, was calculated.

Figure 10: Dependence of on time at different temperatures of process . Left panel: ; right panel: . , integration step . Time average are shown in horizontal solid lines.

One can see that behaves highly irregular. And numbers and heights of observed peaks increase with the growth of temperature until the process becomes chaotic at high . At temperatures process consists of individual rare peaks what can point to the possibility that the energy can be transmitted by impulses.

Iii Sound, solitons and breathers in the -FPU lattice

In accordance with the conjecture on the soliton contribution to the heat conductivity Li_05 (); Vil02 (), an attempt was made to shed some light on this problem. The spatiotemporal correlator was calculated, where is the relative displacement of neighboring particles in time instant . Solitons, being highly correlated displacements of particles, can leave a trace in correlation function. Time shift was fixed and spatial correlation were calculated. Result is shown in Fig. 11.

Figure 11: Correlator as a function of the lattice coordinate . -FPU lattice, , . Arrow points to

Two peaks in the correlation function, shifted by , are visible. The velocity of their propagation is .

iii.1 Sound velocity in the -FPU lattice

The temperature dependence of the sound velocity in the -FPU lattice was discovered about a decade ago Lep98 (). The sound velocity was estimated , where – parameter of renormalized frequencies depending on the temperature. Asymptotic value of the sound velocity in the high temperature limit was derived recently in Li_10 (). If this formula apply to then what differs from found from correlation functions above.

Below we derive more accurate expression for the sound velocity at low temperatures. It was shown Ala01 () that in nonlinear systems there exists a spectrum of frequencies which are proportional to the harmonic ones, according to a well defined law. Then the -FPU potential can be rewritten as

(19)

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

(20)

As a result the lattice becomes harmonic and it is necessary to find . It can be done in terms of mean field approximation (MFA). Mean value of potential energy is

(21)

where is the mean of . According to the virial theorem, mean values of potential and kinetic energies are equal in the harmonic lattice, that is . But the identity holds for 1D systems. Then

(22)

The condition of self consistency of the MFA is

(23)

and it follows that and substitution of this expression into (23) gives the self-consistent equation for :

(24)

with the solution

(25)

Eq. (25) defines the rigidity coefficient for the harmonic lattice with the renormalized spectrum depending on temperature . Thereby the temperature renormalized sound velocity

(26)

and for what coincides with good accuracy with found from correlation functions. The high temperature asymptotic of the sound velocity . The temperature dependence of sound velocity for temperatures in the range is shown in Fig. 13 and very good agreement between analytical and “experimental” results is observed.

Figure 12: The temperature dependence of sound velocity : solid line – dependence (26); empty circles – MD simulation.

If the mean field approximation is applied to the lattices with cubic nonlinearity then no dependence of the sound velocity on temperature is expected. Really, and as .

The obtained results point to the fact that the nonlinearity can be also ignored in describing the thermodynamics of the -FPU lattice, at least at , and the harmonic lattice with renormalized rigidity coefficients (25) is an adequate model. This conjecture was checked at different temperatures by the comparison of the total energies computed in MD simulation of the -FPU lattice and its harmonic model with renormalized rigidity coefficients. Very good coincidence of both energies testifies this hypothesis.

There exists an accurate expression for the specific potential energy (mean potential energy of one oscillator) Lik09 ()

(27)

derived from the thermodynamics of the -FPU lattice; – modified Bessel functions. Specific potential energy (27) and computed in harmonic approximation with rigidity coefficients (25) coincidence with good accuracy.

iii.2 Solitons and breathers in the -FPU lattice

The discrete -FPU lattice can be reduced to the modified Korteweg – de Vries (mKdV) equation in the continuum long-wavelength approximation (see Appendix). The mKdV equation has solutions in the form of solitons and moving breathers Lam80 ().

Solitons of compression and elongation has the form

(28)

where plus/minus signs stand for elongation/compression solitons; – single parameter which simultaneously determines amplitude, width and velocity of soliton; – local deformation of the lattice; – soliton coordinate at time .

The two-parameteric breather solution is

(29)

where

(30)

and are free parameters (see Appendix for more details). It is necessary to make transformation from continuous variables to discrete variables and in an attempt to use soliton and breather solutions on the discrete -FPU lattice.

If soliton (28) and breather (29) solutions exist in the continuum limit, then the question arises: whether these moving localized excitations can be observed on discrete lattice? The answer is ‘yes’ and below the visualization method is suggested.

Recall that the visualization method for standing discrete breathers is well known Bik99 (); Aub06 (): boundaries with friction forces absorb thermal noise (phonons) and standing breathers can be easily seen. Other method is necessary to visualize the moving excitations. Let we have the thermolized lattice with oscillators at temperature . If “cold” lattice (with zero velocities and displacements) is switched to the thermalized lattice then solitons and breathers should “run out” to the cold lattice and could be observed. Results are shown in Fig. 13 and solitons and breather are immediately seen.

Figure 13: Solitons and breather running out of the initially thermalized lattice. A – breather, B – soliton of compression, C – pair of antisolitons (solitons of elongation). Arrow at shows the border separating initially thermalized and “cold” parts of the lattice. Initial temperature of the left lattice part (), .

Solitons and breather (shown in inserts to Fig. 13) move faster then the sound front. At the first glance it is inconsistent with the relation between velocities of sound and solitons. The sound velocity at is what is larger then the maximal soliton velocity (see Appendix). But the temperature of expanding thermal excitations gradually decreases and the sound velocity also decreases according to (26). And there comes a point when solitons, which have constant velocity, keep ahead the sound front.

There are good grounds for believing that solitons and breathers do exist in the -FPU lattice. Very likely that the soliton contribution to the heat conductivity increases with the growth of temperature. Really, the temperature dependence of the soliton density obeys the relation for the Toda lattice at low temperatures Mar91 (). Conceivably the growth of the solitons density with temperature might be an inherent property of nonlinear systems.

Our results on the soliton contribution to the heat conductivity are inconsistent with previous publication Li_10 () where the energy carriers are effective phonons rather than solitons.

The possibility of energy transfer by solitons was conjectured three decades ago Tod79 (). Less studied is the possibility of energy transfer by breathers. One suggested mechanism is the Targeted Energy Transfer Kop01 (); Man04 () 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 Aub06 ().

Iv Conclusions

In conclusion we briefly summarize our results. A new method for the calculation of the heat conductivity is suggested. This is done by the decomposing of the total dynamics into two parts: equilibrium process at equal temperatures of both lattice ends, and nonequilibrium process at temperature of one end and zero temperature of the other. This approach allows to extract and analyze the heat conductivity in an explicit form.

The primary goal of the paper was to develop a method which would allow to decrease the computational time at small temperature gradients when fluctuations of the heat flux are usually too large. It was supposed that at small temperature gradients, when the harmonic approximation is valid and an expression for the heat flux has the form (10), an analytical averaging over random Langevin noise terms can be done. This approach is very efficient for the calculation of quadratic in terms – the gain was thousand-fold. But formulaes are very complex for the linear terms and an efficient algorithm for their realization was not found yet. By this expedient the objective has not been met in full: the gain in computational time is obvious for small lattices, but decreases as the lattice length increases. Nevertheless we suppose that the further analysis of process can be useful as it is responsible for the energy transfer.

The threshold phenomena are familiar in microcanonical ensembles Lic08 (). There exists two values of specific energy . One separates dynamical regime and weak chaos, and higher separates weak and strong chaos. It may be inferred that an existence of threshold phenomena is also a common occurrence in canonical ensembles. Really, a threshold temperature was found. The threshold temperature separates the different behavior of process : process damps out at and reaches the stationary value at . It may be conceived that the soliton and breather contributions to the heat conductivity increases with the growth of temperature if . Solitons and breathers can emerge from either thermal fluctuations or higher order phonon interactions. Additional experiments for nanosized systems at low temperatures when can reveal some new features omitted in the present work.

The modified Korteweg – de Vries equation is derived in the continuum approximation for the -FPU lattice. mKdV has solutions in the form of compression/elongation solitons and breathers. The stability of these quasi-particles was checked in numerical experiments. Both types of excitations were directly visualized.

On the other hand, it was found that the non-linear -FPU lattice can be reduced to the harmonic lattice with the temperature renormalized frequency spectrum. This reduction allows to reproduce adequately the heat conductivity and thermodynamics of the parent lattice. These two, mutually contradictory, properties of the -FPU lattice, – an inherent existence of solitons and breathers, and its “harmonic” behavior, seem to be very strange. Further analysis is necessary to solve this dilemma.

It is likely that some fraction of total heat conductivity is conditioned by solitons and moving breathers. But their contribution to the energy transfer deserves further investigation.

The -FPU lattice is unique in the sense that in the continuum limit it has stable solutions in the form of solitons of compression and elongation. If the number of compression and elongation solitons is equal on average, then no macroscopic changes in the lattice length appear and an additional energy of deformation is negligible. It is an additional energetic factor favoring the solitons existence.

The case is quite different for unsymmetrical potentials which in the lowest order of the Taylor expansions have the form . -FPU, Toda, Morse, Lennard-Jones potentials are the examples. All these “cubic” potentials can be easily reduced to the ordinary KdV equation with the solution in the form of soliton of compression. And the large number of solitons is highly energetically unfavorable due to macroscopical compression of the lattice length.

Appendix A The modifief Korteweg – de Vries equation
for the -FPU lattice. Solitons and breathers.

An approximate solution for the soliton of compression in the -FPU lattice was obtained about two decades ago Wat93 (). Analogous soliton solution with the profile was used recently Li_10 () and this solution is the function of a single parameter – sound velocity .

Below we derive an equation for the continuum analogue of the -FPU lattice with more accurate soliton and breather solutions. The -FPU potential has general form

(31)

where is the relative displacement of neighboring oscillators. The corresponding equations of motion are

(32)

Starting from (32), the continuum approximation can be derived supposing small deviations from equilibrium. The series expansion in terms of up to the forth order is:

(33)

Substitution of this expansion into (32) gives the continuum equation:

(34)

Below we follow the well known reductive perturbation method (RPM) Sas81 (); Leb08 () to get necessary equation in partial derivatives. The receipt consists in introducing new variables

(35)

Next the hierarchy of the expansions in terms of small parameter should be used. If (35) is substituted in (34) and terms of the order and higher are neglected, then and . is the sound velocity in the harmonic approximation. Transformation (35) means that the new coordinate system moves with velocity relative to the old coordinate system .

Equation with the accuracy of the order is

(36)

and one can see that (36) reminds the well known modified KdV equation. Additional variables substitution and should be done to get the exact form of the mKdV equation:

(37)

and are spatial and time variables.

Equation (37) has two types of solutions Lam80 (). The soliton solution is

(38)

Plus/minus signs are related to elongation/contraction solitons, correspondingly. Soliton (38) is the one-parametric solution: free parameter defines simultaneously amplitude (), half-width () and velocity (); defines the soliton coordinate at .

Returning back to initial coordinates , the soliton solution is

(39)

where – coordinate, – time, – soliton coordinate at . The soliton velocity is and is “supersonic” relative to the sound velocity in the harmonic approximation . If increases then soliton has larger amplitude, becomes more narrow and its velocity increases. The soliton solution for the -FPU lattice can be written if discrete variables in (39) are used: , , .

Parameter can be arbitrary large in the continuum limit (39). But the lattice discreteness imposes limitations on the soliton width: solitons with the half-width less then become unstable. That is to say, soliton amplitude and velocity also have upper limit: ,    and the free parameter for discrete -FPU lattice .

Two-parametric breather solution is

(40)

where

(41)

are free parameters; – initial phases; the group and phase velocities are and , correspondingly. Returning back to coordinate system , the breather solution takes the form

(42)

where