Nonequilibrium Generalised Langevin Equation for the calculation of heat transport properties in model 1D atomic chains coupled to two 3D thermal baths
We use a Generalised Langevin Equation (GLE) scheme to study the thermal transport of low dimensional systems. In this approach, the central classical region is connected to two realistic thermal baths kept at two different temperatures [H. Ness et al., Phys. Rev. B 93, 174303 (2016)]. We consider model Al systems, i.e. one-dimensional atomic chains connected to three-dimensional baths. The thermal transport properties are studied as a function of the chain length and the temperature difference between the baths. We calculate the transport properties both in the linear response regime and in the non-linear regime. Two different laws are obtained for the linear conductance versus the length of the chains. For large temperatures ( K) and temperature differences ( K), the chains, with atoms, present a diffusive transport regime with the presence of a temperature gradient across the system. For lower temperatures( K) and temperature differences ( K), a regime similar to the ballistic regime is observed. Such a ballistic-like regime is also obtained for shorter chains (). Our detailed analysis suggests that the behaviour at higher temperatures and temperature differences is mainly due to anharmonic effects within the long chains.
pacs:05.10.Gg, 05.70.Ln, 02.70.-c, 63.70.+h
Understanding the physical properties of low-dimensional thermal transport is an active field of research since such properties are often counter intuitive and present intriguing features AtleeJackson:1968 (); Nakawa:1968 (); Lepri:2003 (); Dahr:2006 (); Dhar:2008 (); Segal:2016 (). A recent review on the subject can be found in [Lepri:2016, ]. For instance, some one-dimensional (1D) models violate the well known Fourier law of heat transport, while some others do not.
It has been known, since the seminal work of Rieder et al. [Rieder:1967, ], that in 1D homogeneous harmonic systems (also referred to as integrable systems), the thermal conductivity diverges in the thermodynamic limit. No temperature gradient is formed in the bulk of the system, since the dominating energy “carriers” are not scattered and propagate ballistically. A large variety of harmonic (integrable) classical Rieder:1967 (); Rich:1975 (); Spohn:1977 (); Casati:1979 (); Hu:2000 (); Segal:2003 (); Segal:2008 (); Kannan:2012 (); Saaskilahti:2012 (); Landi:2013 () and quantum Spohn:1977 (); Zurcher:1990 (); Saito:1996 (); Saito:2000 (); Dahr:2003 (); Segal:2003 (); Gaul:2007 (); Asadian:2013 () systems have been studied using analytical and/or numerical approaches. All these studies show that there is no temperature gradient inside the system (except for small regions in the vicinity of the contacts between the central system and the baths). One usually obtains a constant-temperature profile Rieder:1967 (); Rich:1975 (); Spohn:1977 (); Casati:1979 (); Hu:2000 (); Kannan:2012 (); Landi:2013 (); Zurcher:1990 (); Saito:1996 (); Saito:2000 (); Dahr:2003 (); Gaul:2007 (); Asadian:2013 () in the central system around the averaged temperature .
On the other hand, in classical or quantum non-integrable systems, a temperature gradient is formed inside the system. The temperature gradient is uniform, and the heat conductivity is finite. The transport is said to be diffusive and these systems obey Fourier’s law. In order to obtain a diffusive transport regime, one has to introduce any form of anharmonic effects in the system Jackson:1968 (); Bolsterli:1970 (); Nakazawa:1970 (); Spohn:1977 (); Eckmann:1999 (); Hatano:1999 (); Hu:2000 (); Zhang:2002 (); Pereira:2004 (); Mai:2006 (); Bricmont:2007 (); Gaul:2007 (); Segal:2003 (); Segal:2008 (); Hu:2010 (); Giberti:2011 (); Pereira:2011 (); Saaskilahti:2012 (); Shah:2013 (); Landi:2013 (). Diffusive transport is also obtained when extra local stochastic processes Bolsterli:1970 (); Rich:1975 (); Spohn:1977 (); Zurcher:1990 (); Pereira:2004 (); Bernardin:2005 (); Kannan:2012 (); Landi:2013 () or extra collision processes Lepri:2009 (); Bernardin:2012 () are introduced. A vibrational mode coupling in classical systems Wang:2004 () is also responsible for diffusive transport.
The introduction of configurational defects Jackson:1968 (); Rubin:1971 (); Tsironis:1999 () or disorder Rubin:1971 (); Rich:1975 (); Kipnis:1982 (); Zhang:2002 (); Kannan:2012 () in harmonic systems can also lead to diffusive transport, with the build up of a temperature gradient across the system. Defects and disorder introduce some form of localisation of the vibration modes Nakazawa:1970 (); Stoneham:1975 () which do not favour ballistic transport Nakazawa:1970 ().
In the previously mentioned studies, the system is either ballistic or diffusive. However it is very important to understand if a crossover between the two regimes can be obtained as it has been observed experimentally in a wide range of low-dimensional systems (from carbon nanotubes Kim:2001 (), graphene nanoribbons Bae:2013 (), nanowires Hsiao:2013 () and polymer nanofibres Shen:2010 ()). From the theoretical point of view, the ballistic to diffusive crossover has been studied phenomenologically by introducing and tuning a viscosity-like coefficient (in a simple Langevin-like dynamics) on each atom in the system Landi:2013 (); Roy:2008 (), by modifying the interatomic potential Segal:2003 (), or by introducing a simple description of phonon-phonon interaction Yamamoto:2009 (). A single scheme that can describe both regimes is central and crucially needed to understand the origin of the crossover. Recently, a unified microscopy formalism to study the ballistic-to-diffusive crossover was provided in Bagchi:2014 () on the basis of a scaling ansatz for model systems.
In this manuscript, we present an application of our recently developed Generalised Langevin Equation (GLE) method to the study of the heat transport properties of low-dimensional systems. Within this single GLE scheme, we simulate the dissipative dynamics of model Al systems. Our main objective is to determine if different transport regimes can be established depending on the chain’s length and temperatures involved.
We consider 1D atomic chains connected to two realistic thermal 3D baths kept at temperatures and . By realistic it is meant that (1) the bath are described at the atomic level, with a proper 3D structure, (2) the coupling between the system and the baths obtained in our GLE approach is not simply given by some arbitrary constant(s), (3) the inter-atomic interaction is given by a N-body type interaction designed in materials science and not by a model pairwise potential. Such a N-body potential gives a realistic description of the covalent bonding between atoms in metallic systems which is not well described by pair-wise potentials. Therefore, we use the Embedded Atom Method (EAM) Daw:1984 (); NISTweb (); Mishin:1999 () to describe the interatomic potential between the Al atoms in the 1D chains and between the chains and the two thermal baths. The only “free parameters” that can be changed are the length of the 1D chains and the temperatures and of the baths. We calculate the transport properties of the 1D atomic chains in both the linear reponse regime and at full nonequilibrium conditions. We study the temperature and length dependence of the linear heat conductance and of the nonequilibrium heat current. We find different transport regimes, which appear to be more characteristic of either the diffusive or ballistic regime, depending on the temperatures and temperature differences and on the length of the chains. We analyse the results in terms of anharmonic versus harmonic effects in the effective interatomic potential.
Ii Method and systems
To study the thermal transport of 1D atomic chains, we used our recently developed GLE scheme Kantorovich:2008 (); Kantorovich:2008b (). The method and its implementation in the molecular dynamics package LAMMPS Plimpton:1995 () have been well documented in recent papers Stella:2014 (); Ness:2015 (). The extension of the original GLE scheme to situations where the central system is interacting with several independent baths has been recently given in Ness:2016 (). We call it the GLE-2B scheme. We stress that this method, at least in principle, enables one to obtain an exact solution of the problem of heat transport in the classical case. The method is based on a physically realistic picture of infinite leads kept at equilibrium (at corresponding temperatures) and an arbitrary central region which interacts and exchanges heat with the leads. The only approximations made are that the leads are treated as harmonic baths and the interaction between the leads and the central system is linear in terms of displacements of the atoms in the leads; the central system is treated without any approximations and can be arbitrarily anharmonic.
In a recent paper [Ness:2016, ], we have discussed the advantages of the GLE-2B approach over other more conventional thermostatting approaches such as Nose-Hoover or simple Langevin thermostats. Here, we use our GLE-2B approach as it does not rely on the use of adjustable parameter(s) to describe the relaxation processes in the baths. In Ref. [Ness:2016, ], we have shown that depending upon the value used for such parameters, one can obtain completely different physical results. In the present work, we want to study exclusively the influence of the length of the system and the baths’ temperatures on the thermal transport properties.
We now briefly recall the main “ingredients” of the method. We consider a central system (the 1D atomic chain) with a general Hamiltonian dynamics for the positions and momenta degrees of freedom (DOF) associated with atom (mass ) and Cartesian coordinate . The central system is connected to two () harmonic baths, with DOF and , via coupling . The coupling is linear with respect to the atomic displacement () of the atom (), in bath , around its equilibrium position. is the mass of atom . The coupling force between the central system and the bath is arbitrary with respect to the central system DOF .
By solving Newton’s equations of motion for all the DOFs, and integrating out the baths’ DOFs, one obtains the following “embedded” dissipative dynamics for the DOF of the central system Kantorovich:2008 (); Ness:2016 ():
The total force acting on atom of the central system (in direction ) arises not only from the interaction between atoms in the central system, but also from the interaction between these atoms and the atoms of the harmonic baths which are kept fixed at their equilibrium positions. The force also contains a polaronic-like term which reflects the fact that, due to the coupling between the central system and the baths, the harmonic oscillators of the baths are displaced Kantorovich:2008 (); Stella:2014 (); Ness:2015 (); Ness:2016 ().
Eq. (1) also contains a generalised memory kernel which depends on both times and separately (i.e. not on their difference) and also on the spatial coordinates (DOFs) of all atoms of the central system. The memory kernel is given by
where is the derivative of the coupling force and the polarization matrix is related to the harmonic dynamical matrix of the bath [Kantorovich:2008, ; Stella:2014, ; Ness:2015, ; Ness:2016, ]. In the energy representation, the polarization matrix represents the vibrational density of states of the corresponding bath [Kantorovich:2008, ; Stella:2014, ; Ness:2015, ; Ness:2016, ]. Finally, the terms contain all the information about the initial conditions of the bath DOFs.
By considering these initial conditions as random processes, the corresponding stochastic forces can be described by a multi-dimensional Gaussian stochastic process with correlation functions Kantorovich:2008 (); Stella:2014 ()
and consequently Eq.(1) corresponds now to a generalised Langevin equation for the central system DOFs, with a non-Markovian memory kernel and coloured noise.
Eq. (1) can be efficiently solved numerically by introducing a set of auxiliary DOFs (aDOF) (the superscript is used to count the aDOFs) Stella:2014 (); Kupferman:2004 (); Bao:2004 (); Luczka:2005 (); Ceriotti:2009 (); Ferrario:1979 (); Marchesoni:1983 () and by mapping the polarization matrix onto a specific analytical form Stella:2014 ()
Each pair of aDOF is associated with the corresponding mapping coefficients . Then it can be shown that solving Eq. (1) is equivalent to solving an extended set of Langevin equations, for the central system DOFs and the aDOFs, as a multivariate Markovian process, where all the DOF are independent Wiener stochastic processes with (white noise) correlation functions Stella:2014 (); Ness:2015 (); Ness:2016 ().
The systems we consider are 1D atomic chains (of length ) connected to 3D thermal baths as shown in Figure 1. The electronic transport properties of similar Al nanowires have been studied some decades ago Wan:1997 (); Taraschi:1998 (); Taylor:2001 (). It is now interesting to study their thermal transport properties using our method. We have taken the Embedded Atom Method Daw:1984 () to model the metallic Al system. The tabulated interatomic potential, provided by the NIST Interatomic Potential Repository Project NISTweb (), is a typical non-pairwise potential.
It is important to note that, when using realistic interatomic interaction potential, the model goes beyond the ball-and-spring (with nearest-neighbour interaction) toy models. In Fig. 1 we show that more than two end-atoms of the 1D chains are coupled directly to the thermal baths. In the present case, the first (last) 3 atoms of the 1D chains are coupled directly to the left (right) thermal bath. Furthermore the coupling of these atoms to the baths is given by the matrix elements which are not constants throughout the GLE dynamics. These matrix elements are explicitly dependent on the positions of the atoms in the 1D chains and change accordingly along the GLE dynamics.
Finally, the thermal baths have their own specific spectral signatures, given by the polarisation matrix (obtained from the Fourier transform of ). The frequency dependence of is not trivial and depends on the very nature (atomic configuration, chemical nature, etc …) of the corresponding bath described at the atomic level. Since in our approach we are dealing with a realistic description of the baths, it is not straight forward to determine if our baths are ohmic (sub- or super-ohmic) as is usually done in other simulations of less realistic models.
We focus our study on the steady state only and consider chains of length atoms. To obtain interesting low-dimensional transport properties, the dynamics of the atoms in the chains is constrained to the 1D motion along the chain axis (-direction) by imposing the condition for every atom in the chain.
We calculate the heat current from the time derivative of the lead Hamiltonian, which results in the sum, over the chain atoms, of the products of their velocity and the corresponding force Lepri:2003 (); Segal:2003 () (see Appendix A for further details).
The GLE-2B simulations are typically performed for 300 ps (150000 timesteps with = 2 fs). The steady state is reached after around 100 ps. We calculate a time average of the heat current over the time range between 200 to 300 ps where the system has reached the steady state (Appendix A). Furthermore, for information, all the results obtained for the heat current flowing across the chains and for all set of temperatures and are given in Figure 9 of Appendix B.
From the slope of at small temperature difference Segal:2003 (), we extract the value of the linear thermal conductance . Figure 2 shows the linear thermal conductance versus the length for different temperatures . We can observe two different regimes for the behaviour of versus : one for short chains (typically ), and the other for longer chains (typically ).
For the longer chains, the length dependence of follows quite well the typical power law as shown in Fig. 3. We have calculated the the function for different values of , the results for are shown in Fig. 3 for the different temperatures corresponding to Fig. 2. Qualitatively speaking, the function is more flat for the longer chains (), i.e. follows a power law in with (diffusive Fourier law). For the shorter chains (), the function clearly presents a stronger dependence on (i.e. is not as flat as for the longer chains) implying that follows a length dependence law different from .
The dependence of the heat current versus length for different sets of baths’ temperatures is shown in Figure 4. For large temperature differences or , the current decreases with the chain length which is a characteristic of a diffusive system. For such large temperatures and temperature differences, one should keep in mind that the system is well beyond the linear response regime, and full non-linear and nonequilibrium effects are present. Determining the true nature of such non-linear effects requires further investigations which are beyond the scope of the present paper. For small temperature differences , the heat current appears almost independent of the chain length (such a behaviour would characterise a ballistic transport regime). The inset in Figure 4 also clearly shows, on a logarithmic scale, that does not saturate for the longer chains, especially for larger temperatures and temperature differences (the blue curves in Fig. 4).
Before we can establish the nature of the transport regime, we now concentrate on the temperature profiles across the chains. As we mentioned in the Introduction, the diffusive transport regime is also associated with the establishment of a linear temperature profile across the system. The temperature profiles for different chain lengths and different temperature differences are shown in Figure 5. One can see the presence of a temperature gradient across the longer chains. The gradient becomes more pronounced for larger temperature differences. [Note that, for the systems considered here (one-dimensional chains connected to 3D baths), most of the temperature drop occurs at the contact between the chain and the bath (i.e. large thermal contact resistance) in contrast to systems with a larger (constant and finite) cross section].
For shorter chains, apart from the case of a very large temperature difference (), the temperature profile across the chain is always flat with the temperature given by , see Fig. 5. As we mentioned in the Introduction, this behaviour is a typical characteristic of the ballistic thermal transport regime. Therefore we can conclude that, by changing the temperature differences and/or the length of the system, we can obtain two different transport regimes.
The short chains appear to behave like a harmonic system with ballistic transport properties and no temperature gradient. For the longer chains and at high temperature differences, the dependence of and vs , and the presence of a temperature gradient across the chain, indicate a more diffusive transport regime. This could be a signature either of an anharmonicity in the systems or of the existence of localised vibration modes. Note that for small temperatures and small temperature differences, all chains appear to behave as ballistic harmonic systems.
It is important now to understand what kind of physics is behind the two transport regimes. For that we check first how the eigenmodes of vibration of the 1D chains change when increases. The eigenmodes of vibration of the chains are obtained from the diagonalisation of the dynamical matrix of the chains connected to the baths, as shown on Figure 1. Figure 6 shows the eigenvalues for four different chain lengths. One observes more (nearly) degenerate modes in the longer chains and a larger number of long wavelength modes in the longer chains, as expected. Furthermore, we do not see any eigenvalues outside the energy spectrum (i.e. or eV) that could correspond to localised modes (i.e. bound states). The presence of such more localised modes (in the longer chains) would lead to the breakdown of the ballistic properties.
Another possible mechanism is related to anharmonic effects in the interatomic potential. In order to understand such effects, we first consider the work of Segal et al. Segal:2003 (). In that work, a simpler bath model and coupling to the bath were used in comparison with our GLE-2B. However, in their model calculations, the authors were able to modify the interatomic potential used for their description of the interaction in the 1D chains. It was found that for purely harmonic chains, the heat current is roughly independent of the chain length (a harmonic chain is an integrable model and presents ballistic properties no matter what the chain length is). When anharmonic effects are introduced in the interatomic potential, becomes length dependent and decreases when chain length increases (see Fig. 10 in Ref. [Segal:2003, ]). Such results are fully compatible with our calculations of shown in Figure 4.
In order to quantify the presence of the anharmonic effects, we consider the averaged atomic displacement . The averaged displacement is obtained from where is the relative displacement of atom in the chain of length . The latter is calculated from where is the time averaged position of atom in the time window , i.e. with .
Figure 7 shows the evolution of versus time when the system has reached the steady state. The results show that, for a given chain length, the average variance increases when the temperature increases, as expected (top panels in Fig. 7). However, for a given couple of temperatures , the average variance is more pronounced in longer chains than in the shorter ones (see the bottom panel in Fig. 7). This suggests that for short chains, the atoms are more confined around the bottom of the potential energy wells, while for longer chains, the atoms have more “room” to move and are able to sample the anharmonic part of the potential energy well.
Fig. 8 shows typical configurations of the atoms in the long chain with . Clearly, during the GLE run, some atoms get closer to each other (with distances smaller than the average interatomic distance), forming some of kind of clusters. The distance between “clusters” is also larger than the average interatomic distance. Such a feature is also characteristic of the presence of long wavelength modes, which are more numerous in longer chains than in the shorter ones.
Therefore, we can conclude that, for longer chains and higher temperatures, the atomic motions sample a larger range of distances and a considerable part of the anharmonic EAM potential. This leads to the expected diffusive transport regime and to the presence of a temperature gradient across the long chains. However this behaviour is less apparent for low temperatures (as shown for example for and in Fig. 4).
For shorter chains, the motion of the atoms is more “restrained” and most probably their motion samples only the harmonic part of the potential. This leads to a more harmonic-like system with a more ballistic-like transport regime with no temperature gradient across the short chains (except in the regime of very large temperature difference ). Note that such a behaviour for the atomic motion could also be “accentuated” by the fact that the three first (last) atoms of the chains are also directly interacting with the atomic baths.
We have used our recently developed Generalised Langevin Equation with 2 baths (GLE-2B) method Ness:2016 () to study the thermal transport properties of 1D atomic chains coupled to two realistic 3D thermal baths kept at their own temperature. The results presented in this paper should be understood as a proof of principle of the robustness and efficiency of the numerical GLE-2B methodology that we have recently developed.
We have found that two different laws are obtained for the linear conductance versus the length of the chains. Furthermore, for large temperatures and temperature differences, the chains present a diffusive transport regime with the presence of a temperature gradient across the system. In such a regime, nonequilibrium effects are present and require an in-depth analysis. For lower temperatures and temperature differences, a regime reminiscent to the ballistic regime is observed. In short chains, except for the largest temperature differences considered, the temperature profile does not present any gradient, a characteristic of a ballistic transport property. Our detailed analysis suggests that the increase in anharmonic effects at higher temperatures/temperature differences is mainly responsible for the diffusive transport regime in the longer chains.
Acknowledgements.We acknowledge financial support from the UK EPSRC, under Grant No. EP/J019259/1.
Appendix A Calculation of the heat current
For the systems considered here, the Hamiltonian of the central region is given by,
where and are the momentum and position of the DOFs of the central region. The Hamiltonian of the harmonic bath is given by
with and being the the momentum and position of the DOFs of the bath , and is the harmonic potential energy of the bath. The coupling between the bath and the central region is given by
where the coupling force between the central system and the bath is arbitrary with respect to the central system DOF.
The EOM for the DOF in the central region are given by
Now, we define the heat current flowing between the central region and the bath as follows:
This definition arises from the local continuity equation (between the Hamiltonian density and the corresponding flux Lepri:2003 ()) integrated over the volume encompassing the bath DOFs. The volume integration of the time derivative of the Hamiltonian density gives the RHS of Eq. (8). The volume integration of the divergence of the flux is transformed into a surface integral of the flux over the interface between the bath and the central region. The latter surface integral (with the surface normal vector pointing from the bath towards the central region) provides the total heat current flowing through the interface between the bath and the central region.
Using elementary calculus ( and from the definition of and , one finds that
By using the EOM Eq. (6) of the bath DOF, the above equation reduces to
Note that, by definition, is the spatial derivative of the force between the DOFs and , and are small displacements of the bath DOFs around their equilibrium positions. Therefore the quantity is a force and its sum can be seen as the total force acting on the DOF due to its coupling to the bath . Consequently the hear current can be expressed as the sum of the products of velocity times force:
For our 1D system, Eq. (10) simply reduces to .
In the steady-state, the current is supposed to be the same (up to thermal fluctuations) in between any pair of atoms in the 1D chain. In practice, we calculate the heat current between each pair of atoms contained in the 1D chains. To make further connections with previous studies Lepri:2003 (); Segal:2003 (), we use the following notations and . The heat current, between the pair , is given by , and we average over all pairs of atoms of the chain.
We perform a further average over time, in the time range typically 200 to 300 ps, where the system has reached the steady-state to get the steady state heat current .
Note that, when the force on atom in the chain is given by a (short-ranged) pairwise potential, it can be decomposed into two contributions from the nearest-neighbours, i.e. with the nearest-neighbours forces . Our definition of the heat current becomes then equivalent (after a few manipulation of the indices ) to the more commonly used expression for the current usually found in the literature Lepri:2003 (); Segal:2003 ().
Appendix B Heat current
For information, we present in Figure 9 most of the results we have obtained for the heat current flowing across all chains (length ) and for all sets of temperatures and . From these results, it can be seen that the heat current increases with increasing temperature differences (as expected). Furthermore, decreases when the length of the chain increases; such a behaviour is more pronounced for large temperature differences .
Furthermore, we can estimate an error on the numerical values calculated for the heat current . We recall that the latter is calculated as an average of the heat current between the pairs in the chain:
The local heat current is obtained from the time average of over the time period with .
From the different simulations, we estimate an absolute error of [eV/ps] for the local current. As the heat current is an average over the different pairs, we assume that the corresponding standard error of the mean behaves as [eV/ps].
Within this error margin, we can safely assume that the heat current calculated for and and for and (shown in Fig. 4 by the black curves with circles and triangles respectivelly) is nearly independent of the chain length.
- (1) E. Atlee Jackson, J. R. Pasta and J. F. Waters, J. Comp. Phys. 2, 207 (1968).
- (2) H. Nakawa, Supp. Prog. Theor. Phys. 45, 231 (1970).
- (3) S. Lepri, R. Livi and A. Politi, Phys. Rep. 377, 1 (2003).
- (4) A. Dhar and D. Roy, J. Stat. Phys. 125, 805 (2006).
- (5) A. Dhar, Adv. Phys. 57 457 (2008).
- (6) D. Segal and B. K. Agarwalla, Annu. Rev. Phys. Chem. 67, 185 (2016).
- (7) S. Lepri (Edt), Thermal Transport in Low Dimensions - From Statistical Physics to Nanoscale Heat Transfer, Notes in Physics, Volume 921 (Springer International Publishing, Switzerland, 2016).
- (8) Z. Rieder, J. L. Lebowitz and E. Lieb, J. Math. Phys. 8, 1073 (1967).
- (9) M. Rich and W. M. Visscher, Phys. Rev. B 11, 2164 (1975).
- (10) H. Spohn and J. L. Lebowitz, Commun. Math. Phys. 54, 97 (1977).
- (11) G. Casati, Nuovo Cimento 52, 257 (1979).
- (12) B. Hu and B. Li and H. Zhao, Phys. Rev. E 61, 3828 (2000).
- (13) D. Segal, A. Nitzan and P. Hänggi, J. Chem. Phys. 119, 6840 (2003).
- (14) D. Segal, J. Chem. Phys. 128, 224710 (2008).
- (15) V. Kannan, A. Dhar and J. L. Lebowitz, Phys. Rev. E 85, 041118 (2012).
- (16) K. Sääskilahti, J. Oksanen, R. P. Linna and J. Tulkki, Phys. Rev. E 86, 031107 (2012).
- (17) G. T. Landi and M. J. de Oliveira, Phys. Rev. E 87, 052126 (2013).
- (18) U. Zürcher and P. Talkner, Phys. Rev. A 42, 3278 (1990).
- (19) K. Saito, S. Takesue and S. Miyashita, Phys. Rev. E 54, 2404 (1996).
- (20) K. Saito, S. Takesue and S. Miyashita, Phys. Rev. E 61, 2397 (2000).
- (21) A. Dahr and B. Sriram Shastry, Phys. Rev. B 67, 195405 (2003).
- (22) C. Gaul and H. Büttner, Phys. Rev. E 76, 011111 (2007).
- (23) A. Asadian, D. Manzano, M. Tiersch, and H. J. Briegel, Phys. Rev. E 87, 012109 (2013).
- (24) E. A. Jackson. J. R. Pasta and J. F. Waters, J. Comput. Phys. 2, 207 (1968).
- (25) R. J. Rubin and W. L. Greer, J. Math. Phys. 12, 1686 (1971).
- (26) M. Bolsterli, M. Rich, and W. M. Visscher, Phys. Rev. A 1, 1086 (1970)
- (27) H. Nakazawa, Prog. Theor. Phys. Suppl. 45, 231 (1970).
- (28) A. M. Stoneham, Theory of Defects in Solids: Electronic Structure of Defects in Insulators and Semiconductors (Oxford University Press, Oxford, 1975). Chapter 11.
- (29) J.-P. Eckmann, C.-A. Pillet and L. Rey-Bellet, Commun. Math. Phys. 201, 657 (1999).
- (30) T. Hatano, Phys. Rev. E 59, R1 (1999).
- (31) Y. Zhang and H. Zhao, Phys. Rev. E 66, 026106 (2002).
- (32) E. Pereira and R. Falcao, Phys. Rev. E 70, 046105 (2004).
- (33) T. Mai and O. Narayan, Phys. Rev. E 73, 061202 (2006).
- (34) J. Bricmont and A. Kupiainen, Commun. Math. Phys. 274, 555 (2007).
- (35) T. Hu and Y. Tang, J.Phys.Soc.Jpn. 79, 064601 (2010).
- (36) C. Giberti and L. Rondoni, Phys. Rev. E 83, 041115 (2011).
- (37) E. Pereira, Physica A 390, 4131 (2011).
- (38) T.N. Shah and P.N. Gajjar, Commun. Theor. Phys. 59, 361 (2013).
- (39) G.P. Tsironis, A. R. Bishop, A. V. Savin and A. V. Zolotaryuk, Phys. Rev. E 60, 6610 (1999).
- (40) C. Kipnis, C. Marchioro and E. Presutti, J. Stat. Phys. 27, 65 (1982).
- (41) C. Bernardin and S. Olla, J. Stat. Phys. 121, 271 (2005).
- (42) S. Lepri, C. Mejía-Monasterio and A. Politi, J. Math. A: Math. Theor. 42, 025010 (2009).
- (43) C. Bernardin, V. Kannan, J.L. Lebowitz and J. Lukkarinen, J. Stat. Phys. 146, 800 (2012).
- (44) E.B. Davis, J. Stat. Phys. 18, 161 (1978).
- (45) J.-S. Wang and B. Li, Phys. Rev. E 70, 021204 (2004).
- (46) P. Kim, L. Shi, A. Majumdar, and P. L. McEuen, Phys. Rev. Lett. 87, 215502 (2001).
- (47) M.-H. Bae, Z. Li, Z. Aksamija, P. N. Martin, F. Xiong, Z.-Y. Ong, I. Knezevic and E. Pop, Nat. Commun. 4, 1734 (2013).
- (48) T.-K. Hsiao, H.-K. Chang, S.-C. Liou, M.-W. Chu, S.-C. Lee and C.-W. Chang, Nat. Nanotechnol. 8, 534 (2013).
- (49) S. Shen, A. Henry, J. Tong, R. Zheng and G. Chen, Nat. Nanotechnol. 5, 251 (2010).
- (50) D. Roy, Phys. Rev. E 77, 062102 (2008).
- (51) T. Yamamoto, S. Konabe, J. Shiomi and S. Maruyama, Appl. Phys. Express 2, 095003 (2009).
- (52) D. Bagchi and P. K. Mohanty, J. Stat. Mech. 11, P11025 (2014).
- (53) M. S. Daw and M. I. Baskes, Phys. Rev. B 29, 6443 (1984).
- (54) We have use the data file Al99.eam.alloy from the NIST Interatomic Potential Repository Project (http://www.ctcms.nist.gov/potentials/) as provided by Ref. [Mishin:1999, ].
- (55) Y. Mishin, D. Farkas, M.J. Mehl, and D.A. Papaconstantopoulos, Phys. Rev. B 59, 3393 (1999).
- (56) L. Kantorovich, Phys. Rev. B 78 094304 (2008).
- (57) L. Kantorovich and N. Rompotis, Phys. Rev. B 78 094305 (2008).
- (58) S. J. Plimpton, J. Comput. Phys. 117, 1 (1995)
- (59) L. Stella, C. D. Lorenz and L. Kantorovich, Phys. Rev. B 89, 134303 (2014).
- (60) H. Ness, L. Stella, C. D. Lorenz and L. Kantorovich, Phys. Rev. B 91, 014301 (2015).
- (61) H. Ness, A. Genina, L. Stella, C. D. Lorenz and L. Kantorovich, Phys. Rev. B 93, 174303 (2016).
- (62) R. Kupferman, J. Stat. Phys. 114, 291 (2004).
- (63) J.-D. Bao, J. Stat. Phys. 114, 503 (2004).
- (64) J. Łuczka, Chaos 15, 026107 (2005).
- (65) M. Ceriotti, G. Bussi and M. Parrinello, Phys. Rev. Lett. 102, 020601 (2009).
- (66) M. Ferrario and P. Grigolini, J. Math. Phys. 20, 2567 (1979).
- (67) F. Marchesoni and P. Grigolini, J. Chem. Phys. 78, 6287 (1983).
- (68) C. C. Wan, J.-L. Mozos, J. Wang and H. Guo, Phys. Rev. B 55, R13393(R) (1997).
- (69) G. Taraschi, J.-L. Mozos, C. C. Wan, H. Guo and J. Wang Phys. Rev. B 58, 13138 (1998).
- (70) J. Taylor, H. Guo and J. Wang, Phys. Rev. B 63, 245407 (2001).