PHYSICAL REVIEW E 89, 032102 (2014)

Thermal conductivity of molecular chains with asymmetric potentials of pair interactions Alexander V. Savin1,* and Yuriy A. Kosevich1,2,† 1

Semenov Institute of Chemical Physics, Russian Academy of Sciences, Moscow 119991, Russia 2 Laboratoire d’Energ´etique Mol´eculaire et Macroscopique, CNRS UPR 288, Ecole Centrale Paris, Grande Voie des Vignes, 92295 Chˆatenay-Malabry, France (Received 30 September 2013; published 4 March 2014) We provide molecular-dynamics simulation of heat transport in one-dimensional molecular chains with different interparticle pair potentials. We show that the thermal conductivity is finite in the thermodynamic limit in chains with the potentials that allow for bond dissociation. The Lennard-Jones, Morse, and Coulomb potentials are such potentials. The convergence of the thermal conductivity is provided by phonon scattering on the locally strongly stretched loose interatomic bonds at low temperature and by the many-particle scattering at high temperature. On the other hand, chains with a confining pair potential, which does not allow for bond dissociation, possess anomalous thermal conductivity, diverging with the chain length. We emphasize that chains with a symmetric or asymmetric Fermi-Pasta-Ulam potential or with combined potentials, containing a parabolic and/or a quartic confining potential, all exhibit anomalous heat transport. DOI: 10.1103/PhysRevE.89.032102

PACS number(s): 44.10.+i, 05.45.−a, 05.60.−k, 05.70.Ln

I. INTRODUCTION

The heat conductivity of low-dimensional systems has stimulated intensive studies [1]. The main problem is to derive from the first principles, on the atomic level, the thermal conductivity (TC) κ as the coefficient in Fourier law between the heat flux and temperature gradient J = −κ∇T . In the one-dimensional (1D) case the Fourier law for the chain of N molecules is reduced to J = κ(T+ − T− )/(N − 1)a,

(1)

where T+,− is a temperature of the left or right chain end, respectively, and a is a chain period. For a small temperature difference T = (T+ − T− )  T = (T+ + T− )/2, Eq. (1) allows us to find the dependence of thermal conductivity on the length L = (N − 1)a and temperature T of the chain: κ(N,T ) = J (N − 1)a/T δT , δ = (T+ − T− )/T  1. (2) The Fourier law (1) is fulfilled if the following finite limit exists: κ(T ¯ ) = lim κ(N,T ). N→∞

The chain has a finite TC in this case, whereas it has an anomalous (diverging with the chain length) TC in the case of κ → ∞ for N → ∞. To date, there are numerous works devoted to the numerical modeling of heat transport in 1D lattices. Anomalies of heat transport in 1D nonlinear systems have been well known since the time of the famous work of Fermi, Pasta, and Ulam [2]. In integrable systems (e.g., a harmonic chain, Toda chain, and chain of rigid disks) the heat flux J does not depend on the chain length L, therefore the thermal conductivity diverges: κ(L) ∼ L for L → ∞. Here the heat transport is carried out by noninteracting quasiparticles, so the energy is not dissipated during the heat transport. Nonintegrability

* †

[email protected] [email protected]; [email protected]

1539-3755/2014/89(3)/032102(12)

of the system is a necessary but insufficient condition for the existence of the normal (finite) thermal conductivity. On the examples of a chain with a symmetric Fermi-Pasta-Ulam (FPU) potential [3–5], a disordered harmonic chain [6–8], a diatomic 1D gas of colliding particles [9–11], and the diatomic Toda chain [12], it was shown that some nonintegrable systems can also have anomalous (diverging with the system size) thermal conductivity. Here the thermal conductivity increases as a power function of the length: κ ∼ Lα for L → ∞, with 0 < α < 1. On the other hand, the chain with the on-site potential can have finite thermal conductivity. The convergence of the thermal conductivity with the system size was shown for the Frenkel-Kontorova chain [13,14], for the chain with the sine-Gordon on-site potential [15], for the chain with the φ 4 on-site potential [16,17], and for the chain of hard disks with the substrate potential [18]. The essential feature of all these models is the presence of the external potential, which models the interaction of the chain with the substrate. These systems do not possess translational invariance and the total momentum is not conserved therein. It was suggested in Ref. [12] that the presence of an external potential plays a key role for the convergence of the thermal conductivity in the system. An anomalous thermal conductivity for all isolated one-dimensional lattices was conjectured, where the absence of an external potential leads to the conservation of the system’s total momentum. This hypothesis was refuted in Refs. [19,20], in which it was shown that the isolated chain of coupled rotators (a chain with a periodic interparticle potential) has a normal (finite) thermal conductivity. Anomalous thermal conductivity of isolated chains is connected to a weak scattering of long-wavelength phonons, which have long mean free paths. A mechanism of effective scattering of long-wavelength phonons should always exist for the occurrence of finite thermal conductivity. For the rotator chain [19,20], such a mechanism is the scattering of phonons by rotobreathers. For the chain with the attached side atomic complexes (phonon leads) [21] and for the nanoribbon with dynamically rough edges [22], it is the destructiveinterference scattering of long-wave phonons by the local

032102-1

©2014 American Physical Society

ALEXANDER V. SAVIN AND YURIY A. KOSEVICH

PHYSICAL REVIEW E 89, 032102 (2014)

atomic oscillators, embedded in the phonon continuum of the system. In recent works [23,24], the thermal conductivity of an isolated chain with an asymmetric pair interparticle potential was studied. It is claimed therein that with a certain degree of interaction potential asymmetry, the thermal conductivity measured in nonequilibrium conditions converges in the thermodynamical limit. The authors of Refs. [23,24] attribute the convergence of the thermal conductivity to the nonuniform thermal expansion of the asymmetric chain. Without a doubt, the thermal expansion of the chain itself cannot be the reason for the convergence of the thermal conductivity, because such expansion is also present in the chain with the asymmetric Toda potential, which is an integrable system. The modeling of heat transport in Ref. [23] was performed with the help of nonequilibrium molecular-dynamics (MD) simulation with the Nos´e-Hoover heat baths. The deterministic Nos´e-Hoover thermostat is not intended to simulate heat transfer. Its usage leads to the modeling of the dynamic system, which consists of two attractors connected by a onedimensional chain [25] and may lead to incorrect results [26]. The effective mechanism of the scattering of long-wave acoustic phonons, which provides the convergence of the thermal conductivity, was also not found in Refs. [23,24]. All this makes it necessary to conduct a more detailed modeling of the heat transport in 1D systems. The aim of this work is to simulate the heat transport in molecular chains with asymmetric interparticle pair potentials in the framework of (i) nonequilibrium molecular dynamics using a Langevin thermostat and (ii) equilibrium dynamics using the Green-Kubo approach. It will be shown that the chains with unbounded asymmetric interparticle pair potentials can have finite thermal conductivity. For instance, chains with the asymmetric pair potentials that allow for bond dissociation (similar to the Lennard-Jones and Morse potentials) possess finite thermal conductivity in the thermodynamic limit. Chains with purely repulsive interaction potentials, such as the chain with a Coulomb interaction between nearest neighbors, also possess finite thermal conductivity. The convergence of thermal conductivity in such chains is caused by strong (anomalous) Rayleigh scattering of long-wave phonons by the fluctuations of local bond stretching at low temperature and by many-body scattering at high temperature. On the other hand, the presence of the interparticle pair potential that limits such fluctuations and prevents the transformation of the 1D lattice into a 1D gas of colliding particles results in anomalous (diverging with the system size) thermal conductivity of the 1D chain. Hence chains with the FPU potential, either symmetric or asymmetric, and chains with any combined pair potential, which includes parabolic and/or quartic confining potentials, have infinite thermal conductivity in the thermodynamic limit. The reason is that such chains do not allow for bond dissociation or rupture. From our studies we can conclude that the thermal transport of a one-dimensional chain will be anomalous if the asymmetric interparticle pair potential growth is not slower than the square of the relative distance. II. MODEL

We consider a molecular chain that consists of N units. In dimensionless form, the Hamiltonian of the chain can be

written as H =

N  1 n=1

2

u˙ 2n +

N−1 

V (un+1 − un − r0 ),

(3)

n=1

where un is a dimensionless coordinate of nth node, the dot denotes differentiation with respect to the dimensionless time t, V (r) is a dimensionless interaction potential between nearest neighbors normalized by V (0) = 0, V  (0) = 0, and V  (0) = 1, and r0 is an equilibrium bond length. From now on we will assume that r0 = 1. To simulate the heat transfer in the chain, we will use the stochastic Langevin thermostat. We consider a finite chain of N+ + N + N− nodes. If the interaction potential V (r) does not suppose the possibility of the rupture of the bond, we take the chain with free ends; if the potential admits the rupture, we consider a chain with fixed ends. We put the N+ right boundary nodes in the Langevin thermostat with temperature T+ and N− left nodes in the Langevin thermostat with temperature T− . The corresponding system of equations of motion of the chain is given by u¨ n = −∂H /∂un − γ u˙ n + ξn+ , n  N+ , u¨ n = −∂H /∂un , n = N+ + 1, . . . ,N+ + N, u¨ n = −∂H /∂un − γ u˙ n +

ξn− ,

(4)

n  N+ + N + 1,

where γ = 0.1 is a coefficient of particle velocity relaxation and ξn± denote random white-noise forces, which simulate the interaction with the Langevin thermostat, normalized by the conditions ξn± (t) = 0,

ξn+ (t1 )ξk− (t2 ) = 0,

ξn± (t1 )ξk± (t2 ) = 2γ T± δnk δ(t2 − t1 ). The system of the equations of motion (4) is integrated numerically. After the onset of thermal equilibrium of the chain with the thermostats, a stationary heat flow is established along the chain, following the temperature distribution Tn = u˙ 2n t , and the local heat flux is Jn = jn t ,

jn = − 12 (u˙ n + u˙ n+1 )rn V  (rn − r0 ) + u˙ n hn ,

where hn = 12 u˙ 2n + 12 [V (rn−1 − r0 ) + V (rn − r0 )], rn = un+1 − un (see [1]). Here the angular brackets denote time averaging  1 t f (t)t = lim f (τ )dτ. t→∞ t 0 The following values were used in the numerical simulation: T± = (1 ± 0.1)T , γ = 0.1, N± = 40, and N = 20,40,80, . . . ,20 480. This method of thermalization solves the problem of the thermal boundary resistance. The distributions of the local energy flux Jn and temperature profile Tn along the chain are shown in Fig. 1. In the steady-state regime, the heat flux through each site at the central part of the chain should remain the same, i.e., Jn ≡ J , N+ < n  N+ + N . This property can be employed as a criterion for the accuracy of numerical modeling and can also be used to determine the characteristic time to achieve the steady-state regime and for the calculation of Jn and Tn . Figure 1 shows that the flux is constant along

032102-2

THERMAL CONDUCTIVITY OF MOLECULAR CHAINS WITH . . .

PHYSICAL REVIEW E 89, 032102 (2014)

−5

x 10

V(r)

0.3

(a)

1

0.2

2

Jn

12

0.1

6 0 −0.5

0 40

120

240

360

480

600

680

0.055

0

0.5

r

1

1.5

FIG. 2. (Color online) Form of the symmetric [α = 0 and β = 1 (curve 1)] and the asymmetric [α = 2 and β = 1 (curve 2)] FPU potential (7).

(b)

Tn

we average the found value over 104 independent realizations of the initial thermalization of the chain. Here the problem of the convergence of the heat conductivity is reduced to the problem of the rate of decay of the correlation function c(t) for t → ∞. The chain has a finite conductivity if the rate of decay is sufficient for the convergence of the integral (6) and the chain has the anomalous thermal conductivity otherwise.

0.05

0.045 40

120

240

360

n

480

600

680

FIG. 1. (Color online) Distributions of (a) local heat flux Jn and (b) local temperature Tn in the FPU chain (for the parameters α = 2 and β = 1) for N± = 40, N = 640, T+ = 0.055, and T− = 0.045. The straight line shows the linear temperature gradient that is used in obtaining thermal conductivity κ(N ).

the central part of the chain, suggesting thus that the required regime has been achieved. At the central part of the chain, we observe an almost linear gradient of the temperature distribution, so we can define the thermal conductivity as ¯ N+ +1 − TN+ +N ), κ(N ) = J (N − 1)a/(T

(5)

where a¯ is the average value of the bond length at T = (T+ + T− )/2 (because of the thermal expansion at T > 0, the average bond length exceeds the equilibrium bond length at T = 0: a¯  r0 = 1). Thermal conductivity can also be found with the use of the Green-Kubo formula [27]  t 1 κc = lim lim c(τ )dτ, (6) t→∞ N→∞ N T 2 0 where c(t) = Js (τ )Js (τ − t)τ is the heat flux (auto)correlationfunction and the total heat flux in the chain is Js (t) = n jn (t). To find the correlation function c(t), we consider a finite cyclic chain with N = 104 units, fully immersed in the Langevin thermostat with temperature T . After the establishment of the thermal equilibrium with the thermostat, we disable the interaction of the chain with the thermostat and model further the dynamics of an isolated thermalized chain. To increase the accuracy of the correlation function measurement,

III. FERMI-PASTA-ULAM POTENTIAL

We consider here a chain with the α-β FPU pair potential V (r) = r 2 /2 − αr 3 /3 + βr 4 /4.

(7)

For nonzero α, the FPU potential is an asymmetric one (see Fig. 2). For definiteness, we take the following values of the parameters in the potential (7): α = 2 and β = 1. The FPU potential does not allow for molecular chain rupture, therefore we consider the chain with free ends. The asymmetry of the interatomic pair potential should lead to thermal expansion of the chain. The condition of the free ends of the chain allows for the thermal expansion. The condition of the fixed ends of the chain prevents it from thermal expansion, which effectively weakens the asymmetry of the interatomic pair potential. It should lead to an increase of the thermal conductivity of the chain. Therefore, we will also study, for comparison, the thermal conductivity of a chain with fixed ends. The dependence of the thermal conductivity κ on the length of the internal part N is shown in Fig. 3. As can be seen from this figure, the asymmetry of the FPU potential does not result in the convergence of the thermal conductivity; it only decreases the divergence rate. Thus, at temperature T = 0.05 the thermal conductivity increases with length as N 0.40 in the chain with the symmetric β FPU potential (for α = 0), while it increases as N 0.16 in the chain with the asymmetric α-β FPU potential (with α = 2). This divergence reduction is caused by the fact that in the chain with the asymmetric FPU pair potential a new channel of phonon scattering opens: It is the scattering of phonons by locally strongly stretched loose bonds. The condition of the fixed ends prevents thermal expansion of the chain and effectively reduces the asymmetry of the interparticle pair potential. The number of locally strongly

032102-3

ALEXANDER V. SAVIN AND YURIY A. KOSEVICH

PHYSICAL REVIEW E 89, 032102 (2014)

0

10

6.5

2 1

t−0.84

6

−0.75

t −1

4

N0.40

5

t

c(t)/c(0)

ln(κ)

−0.85

10

5.5

t−0.58

3

−2

4.5

10

4

0.07

N

4

−3

5

3.5

10

0

1

10

2.5

3

2

0.16

N 2

2

10

N

3

10

4

3

4

10

10

FIG. 4. (Color online) Power-law decay of the correlation function c(t) for the chain with an asymmetric FPU potential (7) [α = 2, β = 1, and temperature T = 0.05 (curve 1)] and for the chain with a hyperbolic potential (11) at T = 1: δ = 0 and b = 0 (curve 2); δ = 1 and b = 0 (curve 3); and δ = 0, b = 1, and c = 4 (curve 4). The straight lines give the power-law dependences t −0.58 , t −0.75 , t −0.84 , and t −0.85 .

10

FIG. 3. (Color online) Dependence of the natural logarithm of the thermal conductivity κ on the length N of the internal part of the chain with free ends and asymmetric FPU potential (7) with α = 2, β = 1, and temperature T = 0.025 (curve 1), T = 0.05 (curve 2), and T = 0.1 (curve 3) and for the chain with a symmetric FPU potential (α = 0 and β = 1) at T = 0.05 (curve 4). Curve 5 gives the dependence for a chain with an asymmetric FPU potential with fixed ends and T = 0.05. Straight solid lines show the power-law dependences N 0.07 , N 0.16 , and N 0.40 .

stretched bonds is reduced and the scattering of phonons by such bonds decreases correspondingly. This leads to a significant increase of thermal conductivity (see curves 2 and 5 in Fig. 3). Therefore, the effects associated with the asymmetry of the interparticle pair potential are most pronounced in chains with free ends. Analysis of the behavior of the correlation function c(t) for t → ∞ confirms the divergence of the thermal conductivity in the FPU chain. At t → ∞ the correlation function c(t) decays as a power function t −δ with the exponent δ < 1 (see Fig. 4, curve 1). Therefore, the integral in the Green-Kubo formula (6) diverges.

with the parameter σ = 2−1/6 and binding energy ε = 1/72, and for the Morse pair potential: (9) V (r) = ε[exp(−βr) − 1]2 , √ with the parameter β = 1/ 2ε = 6. For these potentials, one has V (0) = 0, V  (0) = 0, V  (0) = 1, and limr→+∞ V (r) = ε. The characteristic form of these two potentials is shown in Fig. 5. Since the Lennard-Jones (8) and Morse (9) pair potentials have finite binding energies, the chains with these potentials and free ends are not stable with respect to thermal fluctuations. After some time, such chains will necessarily break and heat flux over them will stop. Therefore, the thermal conductivity can be modeled for the chains with such potentials with fixed ends only. For this purpose, the equations of motion (4) should

0.02

1

0.01

IV. LENNARD-JONES AND MORSE POTENTIALS

2

0

We use the following form for the Lennard-Jones pair potential: V (r) = 4ε{[σ/(1 + r)]6 − 1/2}2 ,

10

V(r)

1

3

t

2

10

(8)

0

0.5

r

1

1.5

FIG. 5. (Color online) Form of the Morse (9) and the LennardJones (8) potentials (curves 1 and 2, respectively). The straight solid line sets the binding energy ε = 1/72.

032102-4

THERMAL CONDUCTIVITY OF MOLECULAR CHAINS WITH . . .

PHYSICAL REVIEW E 89, 032102 (2014)

4

10

5

κ

(a) 4 3

10

3

1

6

2

FIG. 7. (Color online) Time dependence of the distribution of particles in the chain with the Lennard-Jones potential with length N = 1000. The density of the chain is d = 1/1.1 and temperature T = 0.002. Each particle is represented as an interval of unit length with the center at xn . Vertical white areas show the loose locally stretched bonds and the time of life of such fluctuations correspond to the height of the area.

2

10

2

3

10

4

10

(b)

10

4

5

3

κ

10

2

10

3 7

1

10

2 1 2

10

3

10

N

4

10

FIG. 6. (Color online) Dependence of the thermal conductivity κ on the length N of the internal part of the chain with the LennardJones pair potential for chain densities (a) d = 1 and (b) d = 2/3. Dependences are shown at T = 0.002, 0.005, 0.02, and 0.2 (curves 1, 2, 3, and 4, respectively). Straight line 5 shows the power law N 0.89 . Curves 6 and 7 give the dependences for the chain with the Morse pair potential at T = 0.002 and 0.02, respectively. The horizontal dashed lines give the values of the thermal conductivity obtained using the Green-Kubo formula.

be imposed by the following boundary conditions: u1 ≡ 0,

uN+ +N+N− ≡ (N+ + N + N− − 1)a,

(10)

where the parameter a  r0 characterizes the density of the chain d = r0 /a. In these chains strong density fluctuations can occur, so we will take into account the interaction of the nth particle with the left edge thermostat only if its coordinate satisfies un (t) < N+ a and consider the interaction of the particle with the right thermostat only if its coordinate satisfies un (t) > (N+ + N )a. The dependence of the thermal conductivity of the chain with the Lennard-Jones pair potential on its length is shown in Fig. 6. We can conclude from the results of the numerical simulation that the chain with the density of d = 2/3 has a finite thermal conduction at low temperatures. The sequence κ(N ) converges at T = 0.002 and 0.005. The trend of convergence is observed at high temperatures as well, but the convergence is slower and is not fully realized for the

used lengths N  10 240. Equilibrium modeling confirms the convergence of thermal conductivity. Here the correlation function c(t) decays exponentially and thus the integral in the Green-Kubo formula always converges. Both methods give the same limiting values of thermal conductivity [see Fig. 6(b)]. For low temperatures, when the thermal transport is provided by phonons, the convergence of the thermal conductivity is caused by anomalous Rayleigh phonon scattering at the fluctuations of the local bond stretching (local force constant decrease). For the particle density d < 1, the loose locally stretched bonds are present in the chain and their time of life increases with a temperature decrease. As one can see in Fig. 7, there are three or four fluctuations of strong bond stretching with time of life tb > 103 in the chain with N = 1000 particles at temperature T = 0.002 and particle density d = 1/1.1. In such a system the phonons will have a finite mean free path, from one loose locally stretched bond to another. The phonon thermal conductivity shows the tendency of the convergence in the thermodynamic limit in the chain with particle density d = 1. The convergence is slower in such a chain than in the chain with density d < 1. In the chain with the nominal density d = 1, the probability of fluctuation of strong bond stretching is lower than in a dilute chain with density d < 1. The temperature dependence of the thermal conductivity κ is shown in Fig. 8. At nominal density d = 1, thermal conductivity is a nonmonotonic function of temperature. Thermal conductivity increases at both T → 0 and T → ∞. This is due to the fact that at (classically) low temperatures the chain behaves almost like a harmonic chain, which has an infinite thermal conductivity, while at high temperatures it behaves like a gas of rigid disks, which is also characterized by anomalous thermal transport. The situation changes, however, in the diluted chain with density d = 2/3. Here the thermal conductivity increases monotonically with temperature. For T → ∞, the thermal conductivity goes to infinity because the system behaves as a gas of rigid disks; however, the thermal conductivity goes to zero for T → 0 in such a system. It is related to the property of the dilute chain with density d < 1 that such a chain has fluctuations of strong local bond stretching, whose time of life

032102-5

ALEXANDER V. SAVIN AND YURIY A. KOSEVICH

PHYSICAL REVIEW E 89, 032102 (2014)

1.5 1200

1

2

4

V(r)

κ

1

800

1

5

0.5

400 0

2 0 0

3 0.01

0.02

T

−0.5

0.03

FIG. 8. (Color online) Temperature dependence of the thermal conductivity κ of the Lennard-Jones chain with length N = 2560 for chain densities d = 1 (curve 1) and d = 2/3 (curve 2).

diverges for T → 0. Phonon transport becomes impossible in this case and therefore the thermal conductivity goes to zero for T → 0 in such a system. Note that the modeling of the thermal conductivity of Ar nanowires performed in Ref. [28] shows that the quasi-threedimensional rod made from the particles interacting through the Lennard-Jones potential has a finite (normal) thermal conductivity. This shows that the convergence of the thermal conductivity, with the system length, is faster in a quasi-threedimensional rod than in the one-dimensional chain. In the chain with the Morse pair potential, thermal conductivity shows the same behavior as in the chain of the Lennard-Jones potential (see Fig. 6). Therefore, all the above results obtained for the Lennard-Jones chain are also valid for the Morse chain.

−1 −2

1

r

2

contrast to the symmetric FPU potential, which is characterized by hard anharmonicity, the hyperbolic potential (11) has a soft anharmonicity. Our modeling shows that thermal transport in the chain with the potential (11) is anomalous. For the symmetric potential (with the parameters δ = 0 and b = 0), at temperature T = 1 the thermal conductivity κ(N ) grows as N 0.30 (see Fig. 10, curve 2) and the correlation function c(t) decays as t −0.84 (see Fig. 4, curve 2). The asymmetry of the

4 5

0.17

N

3

0.22

2

where the parameter δ  0 characterizes the asymmetry of the potential. For r → +∞ and r → −∞, the hyperbolic potential (11) has the asymptotes r − 1 and −(1 + δ)(r + 1), respectively. Thus, for δ = 0 or 1 and b = 0, the potential (11) has the form of a symmetric or asymmetric hyperbole (see Fig. 9, curves 1 and 2). The following potential also has the form of the asymmetric hyperbole:

The potential (12) describes the interparticle interaction with the hard core. It is defined only for the relative displacements r > −1. For r → −1, the interaction energy tends to infinity [energy grows as 1/2(r + 1)] and for r → ∞ the potential increases as a linear function (r − 1)/2 (see Fig. 9, curve 4). In

N

4

(11)

(12)

0.28

ln(κ)

Here we consider the interaction potentials that have the form of the hyperbola with linear asymptotes at r → ±∞:   1 b V (r) = 1 + δ(1 − tanh r) ( 1 + r 2 − 1) − , 2 cosh2 (cr)

r > −1.

0

FIG. 9. (Color online) View of hyperbolic potential (11) with δ = 0 and b = 0 (curve 1); δ = 1 and b = 0 (curve 2); δ = 0, b = 1, and c = 4 (curve 3) and of the asymmetric potential with a hard core (12) (curve 4). Straight lines show the potential asymptotes. The dashed curve 5 shows the Toda potential (13) with the parameter b = 2.5.

V. HYPERBOLIC POTENTIALS

V (r) = r 2 /2(r + 1),

−1

3

N

0.30

N

2

1

2

10

N

3

10

4

10

FIG. 10. (Color online) Dependence of the thermal conductivity κ on the length N of the central part of the two-mass Toda chain (masses m1 = 1,m2 = 2, . . .) at T = 0.2 (curve 1) and of the chain with the hyperbolic potential (11) at T = 1 with the parameters δ = 0 and b = 0 (curve 2); δ = 1 and b = 0 (curve 3); and δ = 0, b = 1, and c = 4 (curve 4). The straight dashed lines give the power-law dependences N 0.22 , N 0.30 , N 0.17 , and N 0.28 .

032102-6

THERMAL CONDUCTIVITY OF MOLECULAR CHAINS WITH . . .

PHYSICAL REVIEW E 89, 032102 (2014) 0

10

800

5 600 −1

c(t)/c(0)

κ

10

1 400

2

4 3

200

−2

10

6 2

10

N

3

10

FIG. 11. (Color online) Dependence of the coefficient of thermal conductivity κ on the length N of the central part of the chain with the hard-core potential (12) at T = 0.05, 0.1, 0.2, and 0.5 (curves 1, 2, 3, and 4, respectively) and the chain with the Toda potential (13) at T = 0.2 taking into account the interaction of the first neighbors (curve 5) and second neighbors (curve 6). The straight lines give the values of the thermal conductivity obtained using the Green-Kubo formula.

potential reduces the rate of divergence. For the asymmetric potential (with the parameters δ = 1 and b = 0), the thermal conductivity κ(N ) grows as N 0.17 (see Fig. 10, curve 3) and the correlation function c(t) decays as t −0.85 (see Fig. 4, curve 3). For the symmetric potential with an additional potential well (with the parameters δ = 0, b = 1, and c = 4), the thermal conductivity κ(N ) grows as N 0.28 (see Fig. 10, curve 4) and the correlation function c(t) decays as t −0.58 (see Fig. 4, curve 4). In Fig. 11 we show the dependence of the thermal conductivity κ on the length of the central part N of the chain with asymmetric potential with a hard core (12). As one can see, at low temperature T = 0.05 thermal conductivity first increases as κ ∝ ln(N ) and then the growth rate starts to slow down. At high temperatures, the thermal conductivity is convergent. This is also confirmed by the behavior of the correlation functions. As t → ∞, the correlation function decays exponentially, c(t) ∝ exp(−λt) (see Fig. 12). Therefore, the Green-Kubo formula (6) gives a finite value of thermal conductivity. Moreover, both approaches (equilibrium and nonequilibrium MD modeling) give the same limiting values for N → ∞ (see Fig. 11), which confirms the correctness of the numerical simulation. VI. TODA CHAIN

For comparison, we also consider the Toda pair potential V (r) = b−2 [exp(−br) + br − 1],

0

4

10

500

1000

t

1500

2000

FIG. 12. (Color online) Exponential decrease of the correlation function c(t) for the chain with an asymmetric hyperbolic potential (12) at T = 0.2.

The chain with the Toda pair potential (13) (the Toda 1D lattice) is a completely integrable system. The MD simulation shows that the values of the boundary temperatures TN+ +1 and TN+ +N and the value of the heat flux J do not depend on the length of the chain’s central part N (there is no energy dissipation because the heat transport is performed by nonlinear excitations of the Toda chain that do not interact with each other). Thus, in view of (5), the thermal conductivity κ(N ) grows as N (see Fig. 11). The Toda chain ceases to be completely integrable if its atoms have different masses. The simplest case is the two-mass chain when the odd sites have the dimensionless mass m1 = 1 and the even sites have mass m2 = m > 1. Modeling of the thermal conductivity of such a chain was carried out in Ref. [12], where it was shown that for m = 2 the thermal conductivity κ diverges as N 0.35 . Our MD simulations show a slower divergence: κ(N ) ∼ N 0.22 for N → ∞ (see Fig. 10, curve 1). Therefore, the isotopic disorder in the chain does not lead to the convergence of thermal conductivity. The convergence of the thermal conductivity in a one-dimensional lattice can be provided only by the nonlinearity of the (pair) interatomic interaction. The Toda chain also ceases to be a completely integrable system if one takes into account the interaction between the next-to-nearest neighbors or, in other words, one uses the following Hamiltonian instead of the Hamiltonian in Eq. (3): N 2 N−k  1 2  u˙ n + H = V (un+k − un − r0 ). (14) 2 n=1 k=1 n=1 In this case the local thermal flux takes the form

(13)

where the parameter b > 0. To be definite, we take the value of the parameter b = 2.5, at which the Toda potential (13) is well approximated by the asymmetric potential (12) at low temperatures (see Fig. 9, curves 4 and 5). 032102-7

1 [(u˙ n + u˙ n+k )rn,k V  (rn,k − r0 )] + u˙ n hn , 2 k=1 2

jn = −

1 2 1 u˙ + [V (rn−k,k − r0 ) + V (rn,k − r0 )], 2 n 2 k=1 2

hn =

ALEXANDER V. SAVIN AND YURIY A. KOSEVICH

PHYSICAL REVIEW E 89, 032102 (2014)

where rn,k = un+k − un . Numerical modeling of heat transport shows that the Toda chain that takes into account the interactions of the next-to-nearest neighbors has a finite thermal conductivity and the growth of κ(N ) saturates for N → ∞ (see Fig. 11, curve 6).

500

400

κ

2

VII. PURELY REPULSIVE POTENTIALS

300

5

Now we consider the thermal conductivity of the chains with the following purely repulsive potentials: V (ρ) = 1/ρ,

(16)

V (ρ) = exp[−b(ρ − 1) ],

(17)

V (ρ) = exp[−b(ρ − 1)],

(18)

V (ρ) = exp[−2b(ρ − 1)] + 2 exp[−b(ρ − 1)],

(19)

V (ρ) = exp[−bρ]/ρ,

(20)

ρ  1, (21) ρ < 1,

where the parameters b > 0 and α > 0. In order to obtain a stable system, we fix the length of the chain and consider the chain of N sites with fixed ends: u1 (t) ≡ 0,

uN (t) ≡ (N − 1)a.

(22)

Then the ground state of the chain with a fixed density d = 1/a will be a homogeneous chain with period a, u0n = (n − 1)a, with n = 1,2, . . . ,N, in which the pair interparticle potential is given by one of the repulsive potentials (15)–(21). In a chain with a fixed density, the repulsive Coulomb potential (15) can be replaced by the asymmetrical hyperbolic potential V (ρ) =

ρ 1 2 + 2 − = r 2 /a 2 (r + a), ρ a a

4

N

2

where ρ is a distance between the interacting particles and the parameter b > 0. The following short-range potential is also purely repulsive:

V (ρ) = exp{−b[(ρ − 1) + α/(ρ − 1)2 ]} for

100

10

3

4

10

10

FIG. 13. (Color online) Thermal conductivity κ versus length N of the internal part of the chain with a repulsive Coulomb potential (15) [for chain period a = 1 (curve 1) and period a = 2 (curve 2)], short-range potential (21) [for the period a = 3 and parameters b = 2.5 and α = 0.05 (curve 3)], and repulsive potential (17) [for the period a = 1 and parameter b = 2.5 (curve 4)] and for a chain with Yukawa potential (20) [for the period a = 1 and parameter b = 1 (curve 5)]. The temperature of the chain is T = 1. The horizontal straight lines give the values of the thermal conductivity obtained using the Green-Kubo formula.

In order to understand the mechanism that provides the convergence of the thermal conductivity we consider the chain with a repulsive short-range potential (21). For definiteness we take the parameters b = 2.5 and α = 0.05. The form of the potential is shown in Fig. 14, curve 1. The repulsion of the particles takes place only when the distance between them becomes ρ < 1. The collision of the two particles with such an interaction will be elastic, but will occur for a finite period of time [see Fig. 15(a)]. If all the interactions between particles can be reduced to the two-body

(23) 2

where the relative displacement r = ρ − a is introduced. Replacement of the Coulomb potential (15) by the asymmetric hyperbolic potential (23) does not change the equations of motion of the chain with fixed ends (22). Therefore, one should expect that the chain with a repulsive Coulomb potential has finite thermal conductivity. Direct modeling of heat transfer shows that thermal conductivity of the chain with a fixed density converges with an increase of the length (see Fig. 13). Analysis of the behavior of the correlation function also confirms the convergence of thermal conductivity: The correlation function c(t) tends to zero exponentially for t → ∞. Both methods of nonequilibrium and equilibrium MD modeling give the same limiting value, for N → ∞, of the thermal conductivity.

V(r)

3

V (ρ) = 0 for

3

200

(15)

V (ρ) = 1/ρ 12 ,

1

1 1

2 3 0 −1

−0.5

0

0.5

r

1

1.5

2

FIG. 14. (Color online) Form of the short-range potential (21) (curve 1, for the parameters b = 2.5 and α = 0.05), the vibro-impact potential (30) (curve 2), and its smooth approximating potential (30) (curve 3), where r = ρ − 1 is the relative displacement.

032102-8

THERMAL CONDUCTIVITY OF MOLECULAR CHAINS WITH . . .

PHYSICAL REVIEW E 89, 032102 (2014)

5

10 8

v =−0.9 3

7

v =0.0 2

(a) 2

v =2.0 1

v3=2.0

5

v =0.0

4

v1=−0.9

κ

6

x

2

4

10

3

1 2 1 0 0

1

2

3

4 2

7

v3=−0.9

6

v =0.0

5

v1=2.0

10

(b) v =1.61

2

N

4

10

5

10

FIG. 16. (Color online) Thermal conductivity κ versus length N of the internal part of the chain with repulsive potential (16) (curve 1, for temperature T = 1), and with potential (19) (curve 2, for b = 6 and T = 10). The chain density is d = 1 (with chain period a = 1/d = 1). The horizontal straight lines give the values of the thermal conductivity, which are obtained using the Green-Kubo formula.

3

4

v2=0.76

x

3

10

3 2 1

v1=−1.27

0 −1 0

1

2

t

3

4

FIG. 15. (Color online) (a) Three consecutive pair collisions and (b) one three-body collision of three particles in a chain with a shortrange potential (21). The lines show the particle position x versus time t. Gray bars show the diameters of the particles |x(t) − x| < 0.5 (particles interact only when their diameters intersect). Particle velocities v1 , v2 , and v3 before and after all the interactions are shown.

collisions only, the momentum would be transmitted along the chain without scattering and the chain would have infinite thermal conductivity. That is the case of the chain consisting of hard disks performing instantaneous elastic collisions. If the collision does not occur instantaneously, but takes a finite interval of time, it makes possible the many-body collisions, for example, the three-particle collision (collision of one particle with a pair of interacting particles). An example of the three-particle collision is shown in Fig. 15(b). As one can see, in the three-particle collision there is no complete transfer of momentum from one extreme particle to another; some of the total energy remains at the central (intermediate) particle. Therefore, the many-body collisions in a one-dimensional chain should lead to scattering of energy and finite thermal conductivity. Modeling of the heat transfer along the chain of particles with the short-range pair interaction potential (21) shows

convergence of the thermal conductivity (see Fig. 13). Analysis of the behavior of the correlation function also confirms the convergence of the thermal conductivity. Both methods of nonequilibrium and equilibrium MD modeling give the same limiting value, for N → ∞, of thermal conductivity. From this we can conclude that the convergence of thermal conductivity in a chain with the short-range repulsive potential (21) is provided by the many-particle collisions. Note that if the collision of the particles occurs instantly, as for absolutely hard disks, the probability of the three-body collisions is equal to zero. Therefore, the chain of such particles will have an infinite (diverging with the system size) thermal conductivity [10]. A chain with a repulsive potential (16) also has a finite thermal conductivity. This potential leads to a more rigid collision than does the Coulomb potential. With this potential, the time of pair collision is shorter, and hence the effects caused by the many-particle collisions should be weaker, and the thermal conductivity must converge more slowly with the length increase. The correlation function c(t) decays exponentially. At temperature T = 1, the Green-Kubo formula gives the finite value of thermal conductivity κc = 110 000. The dependence of the thermal conductivity κ on the chain length N is shown in Fig. 16. As one can see, the convergence of the thermal conductivity should be expected at lengths N ∼ 250 000. The chain with a repulsive potential (17) also has a finite thermal conductivity. Here the correlation function c(t) decays exponentially and the Green-Kubo formula gives the limiting value of κ(N ) for N → ∞, which coincides with the one given by nonequilibrium simulations (see Fig. 13).

032102-9

ALEXANDER V. SAVIN AND YURIY A. KOSEVICH

PHYSICAL REVIEW E 89, 032102 (2014)

In a chain with a fixed density, the pure exponential repulsive potential (18) can be replaced by the Toda potential

70

V (r) = e−br + βe−ba r − e−ba [1 + ba] 60

(24)

where r = ρ − a is the relative displacement. Replacement of the exponential potential (18) by the Toda potential (24) does not lead to a change in the equations of motion in the chain with fixed ends (22). Therefore, the chain with an exponential repulsive potential is a completely integrable system and has an infinite (diverging with the system size) thermal conductivity. In the direct simulation of the heat transfer, the values of the boundary temperatures TN+ +1 and TN+ +N and the value of the heat flux J do not depend on the length of the chain between the thermostats N (the same temperature is set throughout the whole central part of the chain). The repulsive potential, which is a sum of two exponential functions (19), cannot be replaced by the Toda potential (24). The presence of the second exponential function in the potential makes the molecular chain a system that is not exactly integrable. Our MD simulations show that the chain with such a repulsive potential has a finite thermal conductivity and the correlation function decays exponentially. For b = 6 and temperature T = 10, the Green-Kubo formula gives the value of the thermal conductivity κc = 77 800. The dependence of the thermal conductivity κ on chain length N is shown in Fig. 16. As one can see, the convergence of the thermal conductivity should be expected at chain length N ∼ 105 . For a chain with the repulsive potential (20) (the Coulomb potential with screening), MD simulation shows that the chain has a finite thermal conductivity and the correlation function decays exponentially. For b = 1 and temperature T = 1, the Green-Kubo formula gives the value of the thermal conductivity κc = 330. The dependence of the thermal conductivity κ on the chain length N is shown in Fig. 13. Both methods of nonequilibrium and equilibrium MD modeling give the same limiting value, for N → ∞, of the thermal conductivity.

3

50

1

κ

= e−ba [exp(−br) + br − 1],

40

2

30

20

2

10

3

4

10

N

10

FIG. 17. (Color online) Dependence of the coefficient of thermal conductivity κ on the length N of the central part of the chain with the combined potentials (25) (curve 1, for b = 2.5 and T = 1), (26) (curve 2, for T = 0.5), and (28) (curve 3, for T = 0.25).

The sum of the parabolic potential and the hyperbolic potential with a hard core (12) has the form of V (r) = 14 r 2 /(r + 1) + 14 r 2 .

(26)

Modeling of the heat transfer in a chain with the combined potential (26) shows that at T = 0.5 the thermal conductivity κ(N ) grows as ln(N ) (see Fig. 17). Such growth is consistent with the time dependence of the correlation function, which decays slower than t −1 (see Fig. 18).

0

10

VIII. COMBINED ASYMMETRIC POTENTIALS

V (r) = 12 b−2 [exp(−br) + br − 1] + 14 r 2 ,

2

1

−1

10

c(t)/c(0)

Now we consider the potentials that are the sum of the previously considered asymmetric potentials with the symmetric unbounded potential. The additional unbounded potential leads to the stabilization of the interatomic bonds. Below, with the use of an additional parabolic and/or quartic (fourth order) anharmonic potential, we study how this stabilization affects the thermal conductivity of the chain. The sum of the Toda and the parabolic potentials has the following form:

t−1

3

−2

10

(25)

where the parameter b > 0 and r = ρ − 1 is the relative displacement. Simulations show that for b = 2.5 and T = 1, the thermal conductivity of the chain κ(N ) grows as ln(N ) (see Fig. 17). Such an increase is consistent with the observed behavior of the correlation function: The latter decays slower than t −1 (see Fig. 18). Thus the chain with the combined potential (25) has an infinite thermal conductivity.

−3

10

1

10

2

10

t

3

10

FIG. 18. (Color online) Power-law decay of the correlation function c(t) for the chain with the combined potentials (25) (curve 1, for b = 2.5 and T = 1), (26) (curve 2, for T = 0.5), and (28) (curve 3, for T = 0.25). The straight line corresponds to the power law t −1 .

032102-10

THERMAL CONDUCTIVITY OF MOLECULAR CHAINS WITH . . .

The sum of the parabolic and Morse potentials has the form of (27) V (r) = 12 ε[exp(−βr) − 1]2 + 14 r 2 , √ where the parameter β = 1/ 2ε and ε = 1/72. Moleculardynamics simulation of the heat transfer at temperatures T = 0.05 and 0.5 shows that the correlation function c(t) for the chain with the pair potential (27) decays as a power function of time proportional to t −0.9 at t → ∞. Because of this, the thermal conductivity of such chain also diverges in the limit of N → ∞. The sum of the quartic and the hyperbolic potential with a hard core (12) has the form V (r) = 12 r 2 /(r + 1) + 14 r 4 .

(28)

Modeling of the heat transfer in a chain with the combined potential (28) shows that at T = 0.25 the thermal conductivity κ(N) grows as ln(N ) (see Fig. 17). Such growth is consistent with the time dependence of the correlation function, which decays slower than t −1 (see Fig. 18). The nonsmooth vibroimpact potential V (r) = 12 r 2

for

V (r) = ∞ for

r > −1, r  −1

(29)

can be approximated by a smooth combined potential V (r) = 12 r 2 + 0.001/(r + 1)2

(30)

(see Fig. 14). Molecular-dynamics modeling of the heat transfer in the chain with the combined asymmetric pair potential (30) shows that at T = 1 the correlation function c(t) decays as a power function, as t −0.9 for t → ∞. This allows us to conclude that the chain with the vibroimpact interaction potential (30) has an infinite thermal conductivity. Therefore, the chains with the combined pair interaction potentials, containing parabolic and/or quartic confining potential, are always characterized by anomalous heat transport. The form of the potential (30), its parabolicity for small relative displacements, allows one to conclude that the anomalous thermal transport is related here to the long mean free path of small-amplitude phonons. Analyzing all the simulation results, we can conclude that the asymmetry of the pair potential only does not guarantee the convergence of the thermal conductivity. The chains with the asymmetric FPU potential (7) and with the asymmetric hyperbolic potential (11) have infinite thermal conductivity. The asymptotic behavior of the interaction potential in both approaching and separating the particles is very important for the convergence of the thermal conductivity in the thermodynamic limit. If the interaction potential has a hard core (interaction energy between two particles tends to infinity when approaching their centers) and the interaction grows no faster than the distance between the particles in separating them, the chain with such a potential will have a finite thermal conductivity. The single-well potentials (8) and (12) and purely repulsive potentials (15) and (16) belong to such a class of pair potentials. The Morse potential (9) and the short-range potential (21) do not have a hard core, but the chains with these potentials also have finite thermal conductivity. The

PHYSICAL REVIEW E 89, 032102 (2014)

finite binding energy of the pair potential can result in local bond stretching that strongly scatters phonons. On the other hand, the rapid growth of interparticle interaction energy in approaching their centers can result in energy scattering in many-particle collisions. The first and second scenarios are responsible for the convergence of the thermal conductivity at low and high temperatures, respectively. IX. CONCLUSION

Our molecular-dynamics simulations show that the onedimensional chains with unbounded asymmetric pair potentials can have finite thermal conductivity. The numerical simulation using both nonequilibrium (direct simulation of the heat transfer) and equilibrium (computation of the heat current-current correlation function) molecular dynamics clearly shows that the thermal conductivity converges in the chains with the Lennard-Jones (8) and Morse (9) potentials and in the chain with the hyperbolic potential (12). The asymmetric potentials, which lead to the convergence of the thermal conductivity, are characterized by the presence of either a final binding energy, where one branch of the potential is a limited function [the Lennard-Jones, Morse, or short-range repulsive potential (21)], or a hard-core potential [potentials (12), (15), and (16)]. In chains with fixed ends, purely repulsive potentials (15)– (17) can be reduced to single-well asymmetric potentials with an asymptotic linear branch. At high temperatures, a chain with such a potential may be considered as a one-dimensional gas of particles interacting through their collisions. An example of a chain with a short-range repulsive potential (21) shows that the finite conductivity of the gas is related to the energy dissipation during many-particle collisions. Therefore, the chains with such repulsive interparticle potentials also have a finite conductivity. The exception is the potential with exponential repulsion (18), which can be reduced to the Toda potential (24). The Toda chain is a completely integrable system in which the dynamics can always be described by means of noninteracting elementary excitations and therefore it has an infinite thermal conductivity. We emphasize that the Toda chain that takes into account long-range interactions becomes a system that is not completely integrable with finite thermal conductivity. In the combined potentials (25)–(28) and (30), the parabolic or quartic confining potential, which stabilizes the bonds, is present. The local strong bond stretching is not possible also in the chain with the FPU potential (7) and the chain with the vibroimpact potential (30). There is no anomalous scattering of phonons at the local stretching at low temperatures and there is no one-dimensional gas of particles, interacting through their collisions at high temperatures. Therefore, the thermal conductivity is divergent in the thermodynamic limit. However, the divergence can be very slow and can appear in large lengths only. That is why it is stated in Refs. [23,24] that the thermal conductivity can converge in the chains with such asymmetric pair potentials. More detailed modeling allows us to unambiguously conclude that the thermal conductivity diverges in such chains. Therefore, the chain with an asymmetric FPU potential always has an infinite thermal conductivity.

032102-11

ALEXANDER V. SAVIN AND YURIY A. KOSEVICH

PHYSICAL REVIEW E 89, 032102 (2014)

Note that the Morse potentials are commonly used in molecular dynamics to describe the stiff valence bonds and the Lennard-Jones and Coulomb potentials are used to describe the soft nonvalence bonds. Deformation of the valence and torsion angles is described by the bounded angular periodic potentials. We have shown that chains with such potentials have finite thermal conductivities. Therefore, it is expected that the quasi-one-dimensional molecular systems with such potentials must have finite thermal conductivities, but the convergence can occur at very large lengths. So far these lengths have not been reached in MD simulations of carbon nanotubes [29,30] and nanoribbons [22,31], which implies that the nanotubes and nanoribbons have infinite thermal conductivity. On the other hand, the convergence of the thermal conductivity was obtained for softer quasi-one-dimensional molecular structures, such as the double helix of DNA [32]. Thus, our numerical simulations show that chains with unbounded asymmetric potentials, which allow for strong local bond stretching (Lennard-Jones, Morse, and Coulomb potentials), are characterized by a finite thermal conductivity. The convergence of the thermal conductivity is due to scattering of phonons by locally strongly stretched bonds at low temperature and is due to many-particle collisions at high temperature. On the other hand, if the pair interactions limit strong local fluctuations and do not allow for strong local bond stretching, the chain will have an anomalous (diverging with the system

size) thermal conductivity. The thermal conductivity diverges in the thermodynamic limit in a chain with an asymmetric FPU potential and in a chain with any combined potential, which has a parabolic and/or a quartic confining potential as a component. We can conclude from our simulations that the thermal conductivity of the chain will diverge if the energy of the pair particle interaction does not grow slower than the square of the interparticle distance. Note added. Recently, we became aware of the two works [33,34] in which similar conclusions on anomalous thermal conductivity in the Fermi-Pasta-Ulam chain with asymmetric interparticle potential are reported.

[1] S. Lepri, R. Livi, and A. Politi, Phys. Rep. 377, 1 (2003). [2] E. Fermi, J. Pasta, and S. Ulam, Los Alamos National Laboratory Report No. LA-1940, 1955 (unpublished). [3] S. Lepri, R. Livi, and A. Politi, Phys. Rev. Lett. 78, 1896 (1997). [4] S. Lepri, R. Livi, and A. Politi, Physica D 119, 140 (1998). [5] S. Lepri, R. Livi, and A. Politi, Europhys. Lett. 43, 271 (1998). [6] R. Rubin and W. Greer, J. Math. Phys. 12, 1686 (1971). [7] A. Casher and J. L. Lebowitz, J. Math. Phys. 12, 1701 (1971). [8] A. Dhar, Phys. Rev. Lett. 86, 5882 (2001). [9] A. Dhar, Phys. Rev. Lett. 86, 3554 (2001). [10] A. V. Savin, G. P. Tsironis, and A. V. Zolotaryuk, Phys. Rev. Lett. 88, 154301 (2002). [11] P. Grassberger, W. Nadler, and L. Yang, Phys. Rev. Lett. 89, 180601 (2002). [12] T. Hatano, Phys. Rev. E 59, R1 (1999). [13] B. Hu, B. Li, and H. Zhao, Phys. Rev. E 57, 2992 (1998). [14] A. V. Savin and O. V. Gendelman, Phys. Rev. E 67, 041205 (2003). [15] G. P. Tsironis, A. R. Bishop, A. V. Savin, and A. V. Zolotaryuk, Phys. Rev. E 60, 6610 (1999). [16] B. Hu, B. Li, and H. Zhao, Phys. Rev. E 61, 3828 (2000). [17] K. Aoki and D. Kuznezov, Phys. Lett. A 265, 250 (2000). [18] O. V. Gendelman and A. V. Savin, Phys. Rev. Lett. 92, 074301 (2004). [19] C. Giardina, R. Livi, A. Politi, and M. Vassalli, Phys. Rev. Lett. 84, 2144 (2000).

[20] O. V. Gendelman and A. V. Savin, Phys. Rev. Lett. 84, 2381 (2000). [21] Yu. A. Kosevich, Phys. Usp. 51, 848 (2008). [22] Yu. A. Kosevich and A. V. Savin, Europhys. Lett. 88, 14002 (2009). [23] Y. Zhong, Y. Zhang, J. Wang, and H. Zhao, Phys. Rev. E 85, 060102(R) (2012). [24] S. Chen, Y. Zhang, J. Wang, and H. Zhao, arXiv:1204.5933. [25] A. Fillipov, B. Hu, B. Li, and A. Zeltser, J. Phys. A: Math. Gen. 31, 7719 (1998). [26] F. Legoll, M. Luskin, and R. Moeckel, Nonlinearity 22, 1673 (2009). [27] R. Kubo, M. Toda, and N. Hashitsume, in Statistical Physics II, edited by P. Fulde, Springer Series in Solid-State Sciences Vol. 31 (Springer, Berlin, 1991). [28] Y. Chen, D. Li, J. Yang, Y. Wu, J. R. Lukes, and A. Majumdar, Physica B 349, 270 (2004). [29] A. V. Savin, B. Hu, and Y. S. Kivshar, Phys. Rev. B 80, 195423 (2009). [30] A. V. Savin, Yu. A. Kosevich, and A. Cantarero, Phys. Rev. B 86, 064305 (2012). [31] A. V. Savin, Y. S. Kivshar, and B. Hu, Phys. Rev. B 82, 195422 (2010). [32] A. V. Savin, M. A. Mazo, I. P. Kikot, L. I. Manevitch, and A. V. Onufriev, Phys. Rev. B 83, 245406 (2011). [33] L. Wang, B. Hu, and B. Li, Phys. Rev. E 88, 052112 (2013). [34] S. G. Das, A. Dhar, and O. Narayan, arXiv:1308.5475.

ACKNOWLEDGMENTS

Y.A.K. is grateful to J. F. R. Archilla, A. Cantarero, Y. Chalopin, Yu. Gaididei, L. I. Manevitch, M. Velarde, and S. Volz for useful discussions. A.V.S. thanks the Joint Supercomputer Center of the Russian Academy of Sciences for the use of computer facilities. Y.A.K. thanks the NANOTHERM project of the European Commission for financial support through Grant Agreement No. 318117 and the Spanish Ministry of Economy and Competitivity for financial support through Grant No. CSD2010-0044. Y.A.K. acknowledges the Ecole Centrale Paris for hospitality.

032102-12