Nonequilibrium Generalised Langevin Equation for the calculation of heat transport properties in model 1D atomic chains coupled to two 3D thermal baths
NNonequilibrium Generalised Langevin Equation for the calculation of heat transportproperties in model 1D atomic chains coupled to two 3D thermal baths
H. Ness, a) L. Stella, C.D. Lorenz, and L. Kantorovich Department of Physics, Faculty of Natural and Mathematical Sciences,King’s College London, Strand, London WC2R 2LS, UK Atomistic Simulation Centre, School of Mathematics and Physics,Queen’s University Belfast, University Road, Belfast BT7 1NN, Northern Ireland,UK
We use a Generalised Langevin Equation (GLE) scheme to study the thermal trans-port of low dimensional systems. In this approach, the central classical region isconnected to two realistic thermal baths kept at two different temperatures [H. Ness et al. , Phys. Rev. B , 174303 (2016)]. We consider model Al systems, i.e. one-dimensional atomic chains connected to three-dimensional baths. The thermal trans-port properties are studied as a function of the chain length N and the temperaturedifference ∆ T between the baths. We calculate the transport properties both in thelinear response regime and in the non-linear regime. Two different laws are obtainedfor the linear conductance versus the length of the chains. For large temperatures( T (cid:38)
500 K) and temperature differences (∆ T (cid:38)
500 K), the chains, with
N > T (cid:46)
500 K) and temperature dif-ferences (∆ T (cid:46)
400 K), a regime similar to the ballistic regime is observed. Such aballistic-like regime is also obtained for shorter chains ( N ≤ a) Electronic mail: [email protected] a r X i v : . [ c ond - m a t . m e s - h a ll ] A p r . INTRODUCTION Understanding the physical properties of low-dimensional thermal transport is an activefield of research since such properties are often counter intuitive and present intriguingfeatures . A recent review on the subject can be found in [7]. For instance, some one-dimensional (1D) models violate the well known Fourier law of heat transport, while someothers do not.It has been known, since the seminal work of Rieder et al. [8], that in 1D homogeneousharmonic systems (also referred to as integrable systems), the thermal conductivity divergesin 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. Alarge variety of harmonic (integrable) classical and quantum systems have beenstudied using analytical and/or numerical approaches. All these studies show that thereis no temperature gradient inside the system (except for small regions in the vicinity ofthe contacts between the central system and the baths). One usually obtains a constant-temperature profile in the central system around the averaged temperature T av =( T L + T R ) / . Diffusive transport is also obtained when extralocal stochastic processes or extra collision processes are introduced. Avibrational mode coupling in classical systems is also responsible for diffusive transport.The introduction of configurational defects or disorder in harmonic sys-tems can also lead to diffusive transport, with the build up of a temperature gradientacross the system. Defects and disorder introduce some form of localisation of the vibrationmodes which do not favour ballistic transport .In the previously mentioned studies, the system is either ballistic or diffusive. However itis very important to understand if a crossover between the two regimes can be obtained as ithas been observed experimentally in a wide range of low-dimensional systems (from carbonnanotubes , graphene nanoribbons , nanowires and polymer nanofibres ). From the the-2retical point of view, the ballistic to diffusive crossover has been studied phenomenologicallyby introducing and tuning a viscosity-like coefficient (in a simple Langevin-like dynamics) oneach atom in the system , by modifying the interatomic potential , or by introducing asimple description of phonon-phonon interaction . A single scheme that can describe bothregimes 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 on the basis of a scaling ansatz for model systems.In this manuscript, we present an application of our recently developed GeneralisedLangevin Equation (GLE) method to the study of the heat transport properties of low-dimensional systems. Within this single GLE scheme, we simulate the dissipative dynamicsof model Al systems. Our main objective is to determine if different transport regimes canbe established depending on the chain’s length and temperatures involved.We consider 1D atomic chains connected to two realistic thermal 3D baths kept at tem-peratures T L and T R . By realistic it is meant that (1) the bath are described at theatomic level, with a proper 3D structure, (2) the coupling between the system and the bathsobtained in our GLE approach is not simply given by some arbitrary constant(s), (3) theinter-atomic interaction is given by a N-body type interaction designed in materials scienceand not by a model pairwise potential. Such a N-body potential gives a realistic descriptionof the covalent bonding between atoms in metallic systems which is not well described bypair-wise potentials. Therefore, we use the Embedded Atom Method (EAM) to describethe interatomic potential between the Al atoms in the 1D chains and between the chains andthe two thermal baths. The only “free parameters” that can be changed are the length N of the 1D chains and the temperatures T L and T R of the baths. We calculate the transportproperties of the 1D atomic chains in both the linear reponse regime and at full nonequi-librium conditions. We study the temperature and length dependence of the linear heatconductance and of the nonequilibrium heat current. We find different transport regimes,which appear to be more characteristic of either the diffusive or ballistic regime, dependingon the temperatures and temperature differences and on the length of the chains. We anal-yse the results in terms of anharmonic versus harmonic effects in the effective interatomicpotential. 3 I. METHOD AND SYSTEMS
To study the thermal transport of 1D atomic chains, we used our recently developedGLE scheme . The method and its implementation in the molecular dynamics packageLAMMPS have been well documented in recent papers . The extension of the originalGLE scheme to situations where the central system is interacting with several independentbaths has been recently given in . We call it the GLE-2B scheme. We stress that thismethod, at least in principle, enables one to obtain an exact solution of the problem ofheat transport in the classical case. The method is based on a physically realistic picture ofinfinite leads kept at equilibrium (at corresponding temperatures) and an arbitrary centralregion which interacts and exchanges heat with the leads. The only approximations madeare that the leads are treated as harmonic baths and the interaction between the leads andthe central system is linear in terms of displacements of the atoms in the leads; the centralsystem is treated without any approximations and can be arbitrarily anharmonic.In a recent paper [61], we have discussed the advantages of the GLE-2B approach overother more conventional thermostatting approaches such as Nose-Hoover or simple Langevinthermostats. Here, we use our GLE-2B approach as it does not rely on the use of adjustableparameter(s) to describe the relaxation processes in the baths. In Ref. [61], we have shownthat depending upon the value used for such parameters, one can obtain completely differentphysical results. In the present work, we want to study exclusively the influence of the lengthof the system and the baths’ temperatures on the thermal transport properties.We now briefly recall the main “ingredients” of the method. We consider a centralsystem (the 1D atomic chain) with a general Hamiltonian dynamics for the positions r iα andmomenta p iα degrees of freedom (DOF) associated with atom i (mass m i ) and Cartesiancoordinate α = x, y, z . The central system is connected to two ( ν = 1 ,
2) harmonic baths,with DOF u b ν and ˙ u b ν , via coupling (cid:80) b ν µ l ν f b ν ( r ) u b ν . The coupling is linear with respect tothe atomic displacement u b ν ( b ν ≡ l ν γ ) of the atom l ν ( γ = x, y, z ), in bath ν , around itsequilibrium position. µ l ν is the mass of atom l ν . The coupling force µ l ν f b ν ( r ) between thecentral system and the bath ν is arbitrary with respect to the central system DOF r iα .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 central4ystem : m i ¨ r iα = F iα − (cid:90) t −∞ d t (cid:48) (cid:88) ν,i (cid:48) α (cid:48) K ( ν ) iα,i (cid:48) α (cid:48) ( t, t (cid:48) ; r ) ˙ r i (cid:48) α (cid:48) ( t (cid:48) )+ (cid:88) ν η ( ν ) iα ( t ; r ) (1)The total force F iα acting on atom i of the central system (in direction α ) arises notonly from the interaction between atoms in the central system, but also from the interactionbetween these atoms and the atoms of the harmonic baths which are kept fixed at theirequilibrium positions. The force F iα also contains a polaronic-like term which reflects thefact that, due to the coupling between the central system and the baths, the harmonicoscillators of the baths are displaced .Eq. (1) also contains a generalised memory kernel which depends on both times t and t (cid:48) separately (i.e. not on their difference) and also on the spatial coordinates (DOFs) r = ( r iα )of all atoms of the central system. The memory kernel is given by K ( ν ) iα,i (cid:48) α (cid:48) ( t, t (cid:48) ; r ) = (cid:88) b ν ,b (cid:48) ν √ µ l ν µ l (cid:48) ν g iα,b ν ( r ( t )) × Π b ν ,b (cid:48) ν ( t − t (cid:48) ) g i (cid:48) α (cid:48) ,b (cid:48) ν ( r ( t (cid:48) )) , (2)where g iα,b ν is the derivative of the coupling force g iα,b ν = ∂f b ν ( r ) /∂r iα and the polarizationmatrix Π b ν ,b (cid:48) ν ( t − t (cid:48) ) is related to the harmonic dynamical matrix of the bath ν [56, 59–61].In the energy representation, the polarization matrix Π b ν ,b (cid:48) ν ( ω ) represents the vibrationaldensity of states of the corresponding bath [56, 59–61]. Finally, the terms η ( ν ) iα contain allthe information about the initial conditions of the bath ν DOFs.By considering these initial conditions as random processes, the corresponding stochas-tic forces η ( ν ) iα can be described by a multi-dimensional Gaussian stochastic process withcorrelation functions (cid:104) η ( ν ) iα ( t ; r ) (cid:105) = 0 (cid:104) η ( ν ) iα ( t ; r ) η ( ν (cid:48) ) i (cid:48) α (cid:48) ( t (cid:48) ; r ) (cid:105) = δ νν (cid:48) k B T ν K ( ν ) iα,i (cid:48) α (cid:48) ( t, t (cid:48) ; r )and consequently Eq.(1) corresponds now to a generalised Langevin equation for the centralsystem 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) { s ( k ν ) ν, ( t ) , s ( k ν ) ν, ( t ) } (the superscript k ν is used to count the aDOFs) and by5 IG. 1. Schematic representation of the different one-dimensional Al chains (length N =7 , , ,
27 from top to bottom) connected to two three-dimensional L and R baths (each made of30 Al atoms). The baths shown here are in fact what is called the reduced bath regions . Theseregions contain atoms b ν which interact with atoms iα of the central region via non-zero matrixelements g iα,b ν . The EAM is used for the interatomic interaction. The three first (last) atoms inthe chains (shown in a different color in the top left corner) are the atoms interacting directly withthe left (right) bath respectively. mapping the polarization matrix onto a specific analytical form Π b ν ,b (cid:48) ν ( t − t (cid:48) ) → (cid:88) k ν A ( k ν ) b ν b (cid:48) ν e −| t − t (cid:48) | /τ kν cos( ω k ν | t − t (cid:48) | ) . Each pair of aDOF { s ( k ν ) ν, , s ( k ν ) ν, } is associated with the corresponding mapping coefficients { τ k ν , ω k ν , c ( k ν ) b ν } . Then it can be shown that solving Eq. (1) is equivalent to solving anextended set of Langevin equations, for the central system DOFs and the aDOFs, as a mul-tivariate Markovian process, where all the DOF are independent Wiener stochastic processeswith (white noise) correlation functions .The systems we consider are 1D atomic chains (of length N ) connected to 3D thermalbaths as shown in Figure 1. The electronic transport properties of similar Al nanowires havebeen studied some decades ago . It is now interesting to study their thermal transportproperties using our method. We have taken the Embedded Atom Method to model themetallic Al system. The tabulated interatomic potential, provided by the NIST Interatomic6otential Repository Project , is a typical non-pairwise potential.It is important to note that, when using realistic interatomic interaction potential, themodel goes beyond the ball-and-spring (with nearest-neighbour interaction) toy models. InFig. 1 we show that more than two end-atoms of the 1D chains are coupled directly to thethermal baths. In the present case, the first (last) 3 atoms of the 1D chains are coupleddirectly to the left (right) thermal bath. Furthermore the coupling of these atoms to thebaths is given by the matrix elements g iα,b ν which are not constants throughout the GLEdynamics. These matrix elements are explicitly dependent on the positions of the atoms inthe 1D chains and change accordingly along the GLE dynamics.Finally, the thermal baths have their own specific spectral signatures, given by the polari-sation matrix Π b ν ,b (cid:48) ν ( ω ) (obtained from the Fourier transform of Π b ν ,b (cid:48) ν ( t − t (cid:48) )). The frequencydependence of Π b ν ,b (cid:48) ν ( ω ) is not trivial and depends on the very nature (atomic configuration,chemical nature, etc ...) of the corresponding bath described at the atomic level. Sincein our approach we are dealing with a realistic description of the baths, it is not straightforward to determine if our baths are ohmic (sub- or super-ohmic) as is usually done in othersimulations of less realistic models. III. RESULTS
We focus our study on the steady state only and consider chains of length N =7 , , , , ,
27 atoms. To obtain interesting low-dimensional transport properties, thedynamics of the atoms in the chains is constrained to the 1D motion along the chain axis( x -direction) by imposing the condition p iy = p iz = 0 for every atom i in the chain.We calculate the heat current J heat ( N, T L , T R ) from the time derivative of the lead Hamil-tonian, which results in the sum, over the chain atoms, of the products of their velocity andthe corresponding force (see Appendix A for further details).The GLE-2B simulations are typically performed for 300 ps (150000 timesteps with ∆ t = 2 fs). The steady state is reached after around 100 ps. We calculate a time average of theheat current J heat ( N, T L , T R ) over the time range between 200 to 300 ps where the systemhas reached the steady state (Appendix A). Furthermore, for information, all the resultsobtained for the heat current J heat flowing across the chains and for all set of temperatures T L and T R are given in Figure 9 of Appendix B.7 K li n [ e V / p s ] / [ K ] x - T = 300T = 400T = 500
FIG. 2. Linear thermal conductance K lin for the 1D chains of different length N and for differenttemperatures T . (Note that, in the linear regime, ∆ T → T L → T R ≡ T ). It appears that thelength dependence of K lin is not the same for the shorter chains (i.e. N <
16) and for the longestchains (
N >
10 15 20 25N50100 f( N ) f(N) = K lin (N) x N α α = 1.0 FIG. 3. The function f ( N ) = K lin ( N, T ) N α is plotted for α =1.0 for the different temperatures.It appears that the function f ( N ) is more constant for the longer chains than for the shorter chains.Hence, for the longer chains, we have K lin ( N ) ∝ /N while for the shorter chains K lin ( N ) does notfollow the same length dependence. Note that the color of the lines correspond to the temperaturesgiven in Fig. 2. The units of f ( N ) are the same as the unit of K lin in Fig. 2, i.e. × − [eV/ps K]. From the slope of J heat ( N, T, ∆ T ) at small temperature difference , we extract the valueof the linear thermal conductance K lin ( N, T ) = lim ∆ T → J heat ( N, T, ∆ T ) / ∆ T . Figure 2shows the linear thermal conductance K lin versus the length N for different temperatures T L → T R . We can observe two different regimes for the behaviour of J heat versus N : one forshort chains (typically N <
N > J h ea t [ e V / p s ] ∆ T = 100 ( T
L,R = 400, 300 ) ∆ T = 100 ( T
L,R = 500, 400 ) ∆ T = 500 ( T
L,R = 800, 300 ) ∆ T = 500 ( T
L,R = 1000, 500 ) ∆ T = 1200 ( T
L,R = 1500, 300 ) ∆ T = 1200 ( T
L,R = 1700, 500 )
10 15 20 25N-8-7-6-5-4 l n J h ea t FIG. 4. Heat current J heat versus length N of the chains for different sets of bath temperatures T L , T R . For small temperature differences ∆ T (cid:28) T L , T R , the heat current is almost independentof the chain length (i.e. ballistic-like regime). For larger temperature differences ∆ T > T L or T R , J heat decreases with the chain lengths. This behaviour is more typical for the diffusive transportregime. The inset shows that J heat , on a logarithmic scale, decreases continuously with increasing N and does not saturate for longer chains. For the longer chains, the length dependence of K lin follows quite well the typical 1 /N power law as shown in Fig. 3. We have calculated the the function f ( N ) = K lin ( N ) N α for different values of α ∼
1, the results for α = 1 are shown in Fig. 3 for the differenttemperatures corresponding to Fig. 2. Qualitatively speaking, the function f ( N ) is moreflat for the longer chains ( N > K lin follows a power law in 1 /N α with α ∼ N ≤ f ( N ) clearly presentsa stronger dependence on N (i.e. is not as flat as for the longer chains) implying that K lin follows a length dependence law different from 1 /N .The dependence of the heat current J heat versus length N for different sets of baths’ tem-peratures T L , T R is shown in Figure 4. For large temperature differences ∆ T > T L or T R , thecurrent J heat 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 thesystem is well beyond the linear response regime, and full non-linear and nonequilibriumeffects are present. Determining the true nature of such non-linear effects requires furtherinvestigations which are beyond the scope of the present paper. For small temperaturedifferences ∆ T (cid:28) T L , T R , the heat current appears almost independent of the chain length9such a behaviour would characterise a ballistic transport regime). The inset in Figure 4also clearly shows, on a logarithmic scale, that J heat 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 thetemperature profiles across the chains. As we mentioned in the Introduction, the diffusivetransport regime is also associated with the establishment of a linear temperature profileacross the system. The temperature profiles for different chain lengths and different temper-ature differences are shown in Figure 5. One can see the presence of a temperature gradientacross the longer chains. The gradient becomes more pronounced for larger temperaturedifferences. [Note that, for the systems considered here (one-dimensional chains connectedto 3D baths), most of the temperature drop occurs at the contact between the chain and thebath (i.e. large thermal contact resistance) in contrast to systems with a larger (constantand finite) cross section].For shorter chains, apart from the case of a very large temperature difference ( T L =1700 , T R = 500), the temperature profile across the chain is always flat with the temperaturegiven by T av = ( T L + T R ) /
2, see Fig. 5. As we mentioned in the Introduction, this behaviour isa typical characteristic of the ballistic thermal transport regime. Therefore we can concludethat, by changing the temperature differences and/or the length of the system, we can obtaintwo different transport regimes.The short chains appear to behave like a harmonic system with ballistic transport proper-ties and no temperature gradient. For the longer chains and at high temperature differences,the dependence of K lin ( N ) and J heat ( N ) vs N , and the presence of a temperature gradientacross the chain, indicate a more diffusive transport regime. This could be a signature eitherof an anharmonicity in the systems or of the existence of localised vibration modes. Notethat for small temperatures and small temperature differences, all chains appear to behaveas ballistic harmonic systems.It is important now to understand what kind of physics is behind the two transportregimes. For that we check first how the eigenmodes of vibration of the 1D chains changewhen N increases. The eigenmodes of vibration of the chains are obtained from the diagonal-isation 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 the10 T [ K ] N = 7N = 15N = 23N = 27 ∆ T = 200 0 5 10 15 20 25 304005006007008009001000 ∆ T = 2000 5 10 15 20 25 30N2004006008001000120014001600 T [ K ] ∆ T = 1200 0 5 10 15 20 25 30N40060080010001200140016001800 ∆ T = 1200
FIG. 5. Temperature profiles of 1D chains of different length N = 7,15,23 and 27 for differentbath temperatures T L (filled triangle up) and T R (filled triangle down). For the shortest chain( N = 7), the temperature profile is always flat (except for the large temperature difference when T L = 1700 and T R = 500) and is characteristic of the ballistic transport. For the longest chain N = 27, the temperature profile in the chain always presents a gradient even for T L = 500 and T R = 300, characteristic of the diffusive transport regime (Fourier law). An intermediate behaviouris obtained for the other chains N = 15 , longer chains, as expected. Furthermore, we do not see any eigenvalues outside the energyspectrum (i.e. ω (cid:46) ω (cid:38)
24 eV) that could correspond to localised modes (i.e. boundstates). The presence of such more localised modes (in the longer chains) would lead to thebreakdown 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. . 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 interatomicpotential used for their description of the interaction in the 1D chains. It was found thatfor 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 whatthe chain length is). When anharmonic effects are introduced in the interatomic potential, J heat becomes length dependent and decreases when chain length increases (see Fig. 10 inRef. [13]). Such results are fully compatible with our calculations of J heat ( N ) shown in11 ω n [eV] N = 7N = 15N = 23N = 27 FIG. 6. Eigenvalues (in eV) of the vibrational modes in the chains with N = 7 , , ,
27 (bottomto top lines with circle). More (nearly) degenerate modes exist in the longer chains. However the“accumulation” of a lot of long wave-length modes (around ω →
0) is not yet achieved for N = 27.There are no obvious bound states outside the energy spectrum ( ω (cid:46) ω (cid:38)
24 eV) that wouldcorrespond to more localised vibrational modes.
Figure 4.In order to quantify the presence of the anharmonic effects, we consider the averagedatomic displacement ∆ x av . The averaged displacement ∆ x av ( t ) is obtained from ∆ x av ( t ) =[ (cid:80) Ni =1 (∆ x i ( t )) /N ] / where ∆ x i ( t ) is the relative displacement of atom i in the chain oflength N . The latter is calculated from ∆ x i ( t ) = x i ( t ) − (cid:104) x i (cid:105) where (cid:104) x i (cid:105) is the time averagedposition of atom i in the time window [ t start , t stop ], i.e. (cid:104) x i (cid:105) = (cid:80) Mm =0 x i ( t start + m ∆ t ) / ( M + 1)with t stop = t start + M ∆ t .Figure 7 shows the evolution of ∆ x av versus time when the system has reached the steadystate. The results show that, for a given chain length, the average variance ∆ x av increaseswhen the temperature increases, as expected (top panels in Fig. 7). However, for a givencouple of temperatures T L,R , the average variance is more pronounced in longer chains thanin the shorter ones (see the bottom panel in Fig. 7). This suggests that for short chains, theatoms are more confined around the bottom of the potential energy wells, while for longerchains, the atoms have more “room” to move and are able to sample the anharmonic partof the potential energy well.Fig. 8 shows typical configurations of the atoms in the long chain with N = 27. Clearly,during the GLE run, some atoms get closer to each other (with distances smaller than the12
00 225 250 275 300time [ps]00.20.40.60.8 ∆ x a v [ Å ] T L = 500, T R = 300T L = 1700, T R = 500N = 7200 210 220 230 240 250 260 270 280 290 300time [ps]00.511.52 ∆ x a v [ Å ] N = 7N = 15N = 27200 225 250 275 300time [ps] 00.40.81.21.6 ∆ x a v [ Å ] T L = 500, T R = 300T L = 1700, T R = 500N = 15T L = 500, T R = 300 FIG. 7. Averaged displacement ∆ x av versus time. The value of ∆ x av increases when the tempera-ture increases (see the two top panels). However such an increase is bigger when one increases thelength of the chains. Hence in short chains, the atoms are more confined around the bottom of thepotential energy wells, while for longer chains, the atoms have more “space” to move and are ableto sample the anharmonic part of the potential well. average interatomic distance), forming some of kind of clusters. The distance between “clus-ters” is also larger than the average interatomic distance. Such a feature is also characteristicof the presence of long wavelength modes, which are more numerous in longer chains thanin the shorter ones.Therefore, we can conclude that, for longer chains and higher temperatures, the atomicmotions sample a larger range of distances and a considerable part of the anharmonic EAMpotential. This leads to the expected diffusive transport regime and to the presence of atemperature gradient across the long chains. However this behaviour is less apparent forlow temperatures (as shown for example for T L = 400 and T R = 300 in Fig. 4).For shorter chains, the motion of the atoms is more “restrained” and most probably theirmotion samples only the harmonic part of the potential. This leads to a more harmonic-likesystem with a more ballistic-like transport regime with no temperature gradient across theshort chains (except in the regime of very large temperature difference ∆ T ). Note that sucha behaviour for the atomic motion could also be “accentuated” by the fact that the threefirst (last) atoms of the chains are also directly interacting with the atomic baths.13 IG. 8. Typical atomic configurations of the long chain with N = 27 during a GLE run with T L = 800 and T R = 300. The configuration in the top panel corresponds to the initial conditions,with equally spaced atoms. During the GLE run, a form of clustering appears where the distancesbetween atoms in (between) the clusters is smaller (larger) than the average interatomic distance,hence sampling the anharmonic parts of the potential. IV. CONCLUSION
We have used our recently developed Generalised Langevin Equation with 2 baths (GLE-2B) method to study the thermal transport properties of 1D atomic chains coupled totwo realistic 3D thermal baths kept at their own temperature. The results presented in thispaper should be understood as a proof of principle of the robustness and efficiency of thenumerical GLE-2B methodology that we have recently developed.We have found that two different laws are obtained for the linear conductance versus thelength of the chains. Furthermore, for large temperatures and temperature differences, thechains present a diffusive transport regime with the presence of a temperature gradient acrossthe system. In such a regime, nonequilibrium effects are present and require an in-depthanalysis. For lower temperatures and temperature differences, a regime reminiscent to theballistic regime is observed. In short chains, except for the largest temperature differencesconsidered, the temperature profile does not present any gradient, a characteristic of aballistic transport property. Our detailed analysis suggests that the increase in anharmoniceffects at higher temperatures/temperature differences is mainly responsible for the diffusivetransport regime in the longer chains. ACKNOWLEDGMENTS
We acknowledge financial support from the UK EPSRC, under Grant No. EP/J019259/1.14 ppendix A: Calculation of the heat current
For the systems considered here, the Hamiltonian of the central region is given by, H C = (cid:88) iα m i p iα + V ( { r iα } ) , (A1)where p iα and r iα are the momentum and position of the DOFs of the central region. TheHamiltonian H ν of the harmonic bath ν = 1 , H ν = (cid:88) b ν µ l ν p b ν + V harm( ν ) ( { u b ν } ) (A2)with p b ν and u b ν being the the momentum and position of the DOFs of the bath ν , and V harm( ν ) is the harmonic potential energy of the bath. The coupling between the bath ν = 1 , H νC = (cid:88) b ν µ l ν f b ν ( { r iα } ) u b ν , (A3)where the coupling force f b ν ( { r iα } ) between the central system and the bath ν is arbitrarywith respect to the central system DOF.In Refs. [56, 59, and 61], it was shown that the equation of motion (EOM) of the DOFsin the bath are given by µ l ν ¨ u b ν = − ∂V harm( ν ) ∂u b ν − µ l ν f b ν ( { r iα } ) . (A4)The EOM for the DOF in the central region are given by m i ¨ r iα = − ∂V ( r ) ∂r iα − (cid:88) ν =1 (cid:88) b ν µ l ν g iα,b ν u b ν (A5)where g iα,b ν = ∂f b ν /∂r iα .Now, we define the heat current J ν flowing between the central region and the bath ν = 1as follows: J ν = ddt ( H ν + H νC ) . (A6)This definition arises from the local continuity equation (between the Hamiltonian densityand the corresponding flux ) integrated over the volume encompassing the bath ν DOFs.The volume integration of the time derivative of the Hamiltonian density gives the RHS of15q. (A6). The volume integration of the divergence of the flux is transformed into a surfaceintegral of the flux over the interface between the bath ν and the central region. The lattersurface integral (with the surface normal vector pointing from the bath towards the centralregion) provides the total heat current J ν flowing through the interface between the bath ν and the central region.Using elementary calculus ( dH ( p, q ) /dt = ˙ p dH/dp + ˙ q dH/dq ) and from the definitionof H ν and H νC , one finds that J ν = (cid:88) b ν ˙ u b ν (cid:32) µ l ν ¨ u b ν + ∂V harm( ν ) ∂u b ν (cid:33) + µ l ν f b ν ˙ u b ν + µ l ν u b ν (cid:88) iα g iα,b ν ˙ r iα By using the EOM Eq. (A4) of the bath DOF, the above equation reduces to J ν = (cid:88) iα ˙ r iα (cid:32)(cid:88) b ν µ l ν g iα,b ν u b ν (cid:33) . (A7)Note that, by definition, µ l ν g iα,b ν is the spatial derivative of the force between the DOFs iα and b ν , and u b ν are small displacements of the bath DOFs around their equilibriumpositions. Therefore the quantity µ l ν g iα,b ν u b ν is a force and its sum F ( ν ) iα = (cid:80) b ν µ l ν g iα,b ν u b ν can be seen as the total force acting on the DOF iα due to its coupling to the bath ν .Consequently the hear current can be expressed as the sum of the products of velocity timesforce: J heat = (cid:88) iα ˙ r iα F ( ν ) iα (A8)For our 1D system, Eq. (A8) simply reduces to J heat = (cid:80) i ˙ x i F ( ν ) ix .In the steady-state, the current is supposed to be the same (up to thermal fluctuations) inbetween any pair of atoms in the 1D chain. In practice, we calculate the heat current betweeneach pair ( n, n + 1) of atoms contained in the 1D chains. To make further connections withprevious studies , we use the following notations ˙ x n ≡ v n and F ( ν ) nx ≡ f n . The heat current,between the pair ( n, n + 1), is given by j n,n +1 = ( v n f n + v n +1 f n +1 ) /
2, and we average j n,n +1 over all pairs of atoms of the chain.We perform a further average over time, in the time range typically 200 to 300 ps, wherethe system has reached the steady-state to get the steady state heat current J heat .16ote that, when the force f n on atom n in the chain is given by a (short-ranged) pairwisepotential, it can be decomposed into two contributions from the nearest-neighbours, i.e. f n = F n − ,n + F n,n +1 with the nearest-neighbours forces F n,m = F m,n . Our definition ofthe heat current (cid:80) n v n f n becomes then equivalent (after a few manipulation of the indices n − , n, n + 1) to the more commonly used expression for the current (cid:80) n ( v n +1 + v n ) F n +1 ,n usually found in the literature . Appendix B: Heat current
For information, we present in Figure 9 most of the results we have obtained for the heatcurrent J heat flowing across all chains (length N = 7 , , , , ,
27) and for all sets oftemperatures T L and T R . From these results, it can be seen that the heat current increaseswith increasing temperature differences ∆ T = T L − T R (as expected). Furthermore, J heat decreases when the length of the chain increases; such a behaviour is more pronounced forlarge temperature differences ∆ T .Furthermore, we can estimate an error on the numerical values calculated for the heatcurrent J heat . We recall that the latter is calculated as an average of the heat current betweenthe pairs ( n, n + 1) in the chain: J heat = 1 N pair N − (cid:88) n =1 (cid:104) j n,n +1 (cid:105) τ = 1 N pair N − (cid:88) n =1 ( (cid:104) v n f n (cid:105) τ + (cid:104) v n +1 f n +1 (cid:105) τ ) / . (B1)The local heat current (cid:104) v n f n (cid:105) τ is obtained from the time average of v n ( t ) f n ( t ) over the timeperiod t ∈ [ t start , t end ] with t end − t start = τ .From the different simulations, we estimate an absolute error of ∆( j n,n +1 ) ∼ .
008 [eV/ps]for the local current. As the heat current J heat is an average over the different pairs, we as-sume that the corresponding standard error of the mean behaves as ∆( J heat ) ∼ . / (cid:112) N pair [eV/ps].Within this error margin, we can safely assume that the heat current calculated for T L = 400 and T R = 300 and for T L = 500 and T R = 400 (shown in Fig. 4 by the black curveswith circles and triangles respectivelly) is nearly independent of the chain length.17
500 1000 1500 ∆ T = T L - T R J h ea t [ e V / p s ] N = 7, T R = 300N = 7, T R = 400N = 7, T R = 500N = 7, T R = 600N = 11, T R = 200N = 11, T R = 300N = 11, T R = 400N = 11, T R = 500N = 11, T R = 600 ∆ T = T L - T R J h ea t [ e V / p s ] N = 15, T R = 300N = 15, T R = 400N = 15, T R = 500N = 19, T R = 300N = 19, T R = 400N = 19, T R = 500 ∆ T = T L - T R J h ea t [ e V / p s ] N = 23, T R = 300N = 23, T R = 400N = 23, T R = 500N = 27, T R = 300N = 27, T R = 400N = 27, T R = 500 FIG. 9. Compilation of most of the results for the heat current J heat through the system, 1D chainsof different length N = 7 , , , ,
27 and for different baths’ temperatures T L and T R . Here thetemperature difference is ∆ T = T L − T R . REFERENCES E. Atlee Jackson, J. R. Pasta and J. F. Waters, J. Comp. Phys. , 207 (1968). H. Nakawa, Supp. Prog. Theor. Phys. , 231 (1970). S. Lepri, R. Livi and A. Politi, Phys. Rep. , 1 (2003). A. Dhar and D. Roy, J. Stat. Phys. , 805 (2006).18
A. Dhar, Adv. Phys.
457 (2008). D. Segal and B. K. Agarwalla, Annu. Rev. Phys. Chem. , 185 (2016). S. Lepri (Edt),
Thermal Transport in Low Dimensions - From Statistical Physics toNanoscale Heat Transfer , Notes in Physics, Volume 921 (Springer International Publishing,Switzerland, 2016). Z. Rieder, J. L. Lebowitz and E. Lieb, J. Math. Phys. , 1073 (1967). M. Rich and W. M. Visscher, Phys. Rev. B , 2164 (1975). H. Spohn and J. L. Lebowitz, Commun. Math. Phys. , 97 (1977). G. Casati, Nuovo Cimento , 257 (1979). B. Hu and B. Li and H. Zhao, Phys. Rev. E , 3828 (2000). D. Segal, A. Nitzan and P. H¨anggi, J. Chem. Phys. , 6840 (2003). D. Segal, J. Chem. Phys. , 224710 (2008). V. Kannan, A. Dhar and J. L. Lebowitz, Phys. Rev. E , 041118 (2012). K. S¨a¨askilahti, J. Oksanen, R. P. Linna and J. Tulkki, Phys. Rev. E , 031107 (2012). G. T. Landi and M. J. de Oliveira, Phys. Rev. E , 052126 (2013). U. Z¨urcher and P. Talkner, Phys. Rev. A , 3278 (1990). K. Saito, S. Takesue and S. Miyashita, Phys. Rev. E , 2404 (1996). K. Saito, S. Takesue and S. Miyashita, Phys. Rev. E , 2397 (2000). A. Dahr and B. Sriram Shastry, Phys. Rev. B , 195405 (2003). C. Gaul and H. B¨uttner, Phys. Rev. E , 011111 (2007). A. Asadian, D. Manzano, M. Tiersch, and H. J. Briegel, Phys. Rev. E , 012109 (2013). E. A. Jackson. J. R. Pasta and J. F. Waters, J. Comput. Phys. , 207 (1968). R. J. Rubin and W. L. Greer, J. Math. Phys. , 1686 (1971). M. Bolsterli, M. Rich, and W. M. Visscher, Phys. Rev. A , 1086 (1970) H. Nakazawa, Prog. Theor. Phys. Suppl. , 231 (1970). A. M. Stoneham,
Theory of Defects in Solids: Electronic Structure of Defects in Insulatorsand Semiconductors (Oxford University Press, Oxford, 1975). Chapter 11. J.-P. Eckmann, C.-A. Pillet and L. Rey-Bellet, Commun. Math. Phys. , 657 (1999). T. Hatano, Phys. Rev. E , R1 (1999). Y. Zhang and H. Zhao, Phys. Rev. E , 026106 (2002). E. Pereira and R. Falcao, Phys. Rev. E , 046105 (2004). T. Mai and O. Narayan, Phys. Rev. E , 061202 (2006).19 J. Bricmont and A. Kupiainen, Commun. Math. Phys. , 555 (2007). T. Hu and Y. Tang, J.Phys.Soc.Jpn. , 064601 (2010). C. Giberti and L. Rondoni, Phys. Rev. E , 041115 (2011). E. Pereira, Physica A , 4131 (2011). T.N. Shah and P.N. Gajjar, Commun. Theor. Phys. , 361 (2013). G.P. Tsironis, A. R. Bishop, A. V. Savin and A. V. Zolotaryuk, Phys. Rev. E , 6610(1999). C. Kipnis, C. Marchioro and E. Presutti, J. Stat. Phys. , 65 (1982). C. Bernardin and S. Olla, J. Stat. Phys. , 271 (2005). S. Lepri, C. Mej´ıa-Monasterio and A. Politi, J. Math. A: Math. Theor. , 025010 (2009). C. Bernardin, V. Kannan, J.L. Lebowitz and J. Lukkarinen, J. Stat. Phys. , 800 (2012). E.B. Davis, J. Stat. Phys. , 161 (1978). J.-S. Wang and B. Li, Phys. Rev. E , 021204 (2004). P. Kim, L. Shi, A. Majumdar, and P. L. McEuen, Phys. Rev. Lett. , 215502 (2001). M.-H. Bae, Z. Li, Z. Aksamija, P. N. Martin, F. Xiong, Z.-Y. Ong, I. Knezevic and E.Pop, Nat. Commun. , 1734 (2013). T.-K. Hsiao, H.-K. Chang, S.-C. Liou, M.-W. Chu, S.-C. Lee and C.-W. Chang, Nat.Nanotechnol. , 534 (2013). S. Shen, A. Henry, J. Tong, R. Zheng and G. Chen, Nat. Nanotechnol. , 251 (2010). D. Roy, Phys. Rev. E , 062102 (2008). T. Yamamoto, S. Konabe, J. Shiomi and S. Maruyama, Appl. Phys. Express , 095003(2009). D. Bagchi and P. K. Mohanty, J. Stat. Mech. , P11025 (2014). M. S. Daw and M. I. Baskes, Phys. Rev. B , 6443 (1984). Y. Mishin, D. Farkas, M.J. Mehl, and D.A. Papaconstantopoulos, Phys. Rev. B 59, 3393(1999). L. Kantorovich, Phys. Rev. B L. Kantorovich and N. Rompotis, Phys. Rev. B S. J. Plimpton, J. Comput. Phys. , 1 (1995) L. Stella, C. D. Lorenz and L. Kantorovich, Phys. Rev. B , 134303 (2014).20 H. Ness, L. Stella, C. D. Lorenz and L. Kantorovich, Phys. Rev. B , 014301 (2015). H. Ness, A. Genina, L. Stella, C. D. Lorenz and L. Kantorovich, Phys. Rev. B , 174303(2016). R. Kupferman, J. Stat. Phys. , 291 (2004). J.-D. Bao, J. Stat. Phys. , 503 (2004). J. (cid:32)Luczka, Chaos , 026107 (2005). M. Ceriotti, G. Bussi and M. Parrinello, Phys. Rev. Lett. , 020601 (2009). M. Ferrario and P. Grigolini, J. Math. Phys. , 2567 (1979). F. Marchesoni and P. Grigolini, J. Chem. Phys. , 6287 (1983). C. C. Wan, J.-L. Mozos, J. Wang and H. Guo, Phys. Rev. B , R13393(R) (1997). G. Taraschi, J.-L. Mozos, C. C. Wan, H. Guo and J. Wang Phys. Rev. B , 13138 (1998). J. Taylor, H. Guo and J. Wang, Phys. Rev. B63