Transport and Helfand moments in the Lennard-Jones fluid. II. Thermal conductivity S. Viscardy, J. Servantie, and P. Gaspard Center for Nonlinear Phenomena and Complex Systems, Universit´e Libre de Bruxelles, 7 Campus Plaine, Code Postal 231, B-1050 Brussels, Belgium 0 0 The thermal conductivity is calculated with the Helfand-moment method in the Lennard-Jones 2 fluid near thetriple point. The Helfand moment of thermal conductivity is here derived for molec- ular dynamics with periodic boundary conditions. Thermal conductivity is given by a generalized n a EinsteinrelationwiththisHelfandmoment. Wecomputethermalconductivitybythisnewmethod J andcompareitwithourownvaluesobtainedbythestandardGreen-Kubomethod. Theagreement is excellent. 1 1 PACSnumbers: 02.70.Ns;05.60.-k;05.20.Dd ] h c I. INTRODUCTION Successful results have already been obtained for sim- e ple systems: diffusion in the Lorentz gas [24] and vis- m cosity in a two-particle system [25]. Furthermore, the - This paper completes the analysis carried out in the hydrodynamic-mode method, aiming to construct at the t companion paper [1] by the study of thermal conduc- a microscopiclevelthehydrodynamicmodeswhicharethe t tivity with the Helfand-moment method [2]. Since the solutions ofthe diffusion or Navier-Stokesequations also s seventies,severalworkshavebeencarriedoutforthecal- . depend on the possibility to define Helfand moments in t a culation of the transport coefficients. Shear viscosity is periodicsystems. Thisconstructionhasalreadybeencar- m the most studied transport property, while thermal con- riedoutfordiffusion inthe Lorentzgas[26]. Bydefining ductivity hasbeenlessstudied. Thefirstcomputationof - theHelfandmoment,wehopefully projecttoextendthis d thermal conductivity goes back to the work of Alder et construction to the other transport processes, in partic- n al. forhard-spheresystems[3]. Soft-spherepotentialsys- ular thermal conductivity. The paper is organized as o temswereconsideredforthefirsttimeafewyearslaterby follows. In Section II, the theoretical results of Helfand c Levesqueetal. [4]usingthewell-knownGreen-Kubofor- [ [2] are briefly outlined. In Section III, we present our mula [5, 6, 7, 8]. Nevertheless, since the autocorrelation Helfand-moment method for the calculation of thermal 1 function of the microscopic flux decreases slowly, espe- conductivity. Section IV gives the results of the molecu- v cially near the triple point, most of the studies devoted 4 lardynamicssimulations. Wecompareourresultsbythe to the computation of thermal conductivity were per- 5 Heldand-momentandGreen-Kubotechniquesandalsoto formedbynonequilibriummethods. Thisisusuallydone 2 the litterature. Finally, conclusions aredrawninSection 1 by fixing the temperature gradient and measuring the V. 0 heatflux [9, 10, 11, 12, 13, 14], exceptin Ref. [15] where 7 the opposite is done, i.e., the heat flux is fixed and the 0 temperature gradient is measured. However, some equi- II. THEORETICAL BACKGROUND / t librium molecular dynamics studies using the standard a Green-kuboformulawerecarriedout[4,12,14,16]. More m In the fifties, the Boltzmann and Enskog theories of recently,assuggestedinthenineties[17,18],thegeneral- transportprocesseswerecompletedbyanalternativeap- - izedEinsteinrelationexpressingthermalconductivityas d proachdevelopedby Green[5, 6], Kubo [7]and Mori[8]. n the variance of the time integral of the microscopic flux This new approach relates the transport coefficients to o was computed [19]. The purpose of the present paper is thetimeautocorrelationfunctionsoftheircorresponding c to extend the Helfand-moment method of the compan- microscopicflux. Inparticular,the thermalconductivity : v ion paper [1] from viscosity to thermal conductivity. As coefficient κ can be expressed as i an application of our method, we calculate the thermal X conductivity of the Lennard-Jonesfluid at a phase point κ= lim 1 ∞ J(κ)(0)J(κ)(t) dt. (1) ar near the triple point. We show that the values obtained N,V→∞kBT2V Z0 D E bytheHelfand-momentmethodandtheGreen-Kuboval- where J(κ) is the microscopic flux associated with ther- uesareinexcellentagreement. Suchamethodbringsnot mal conductivity. This theory plays an important role only some interest as an alternative equilibrium molecu- in numerical simulations. On the other hand, Helfand lar dynamics method, but also in the context of recent aimed to express the transport coefficients in terms of theoriesin nonequilibriumstatisticalmechanics. Indeed, generalized Einstein relations [2] and showed that ther- the escape-rate formalism propose clear and direct rela- mal conductivity is given by tionships between the transport coefficients and quanti- tiescharacterizingthemicroscopicchaossuchastheLya- 1 2 κ= lim G(κ)(t)−G(κ)(0) , (2) punov exponents and fractal dimensions [20, 21, 22, 23]. N,V,t→∞2kBT2Vt(cid:28)h i (cid:29) 2 where G(κ)(t) is the corresponding Helfand moment de- Following a way similar to the case of viscosity, we con- fined as the centroid clude that the Helfand moment which should be used in periodic systems is given by N G(κ)(t)= x (E −hE i) , (3) a a a G(κ)(t) = x (E −hE i) a=1 a a a X a X of the energies of the particles − ∆x(s)(E −hE i)θ(t−t ) a a a s p2 1 Xa Xs Ea ≡ 2ma + 2 u(rab). (4) − 1 tdτ L pa+pb ·F(r ). (9) bX(6=a) 4 a6=bZ0 b|ax m ab X It can be proved that the Green-Kubo formula (1) and The derivation is given in Appendix A. As for shear vis- theHelfandequation(2)areequivalentifthemicroscopic cosity,thefirstextratermisduetothejumpsofparticles flux is related to the Helfand moment by from a boundary to the opposite one. The second extra term takes into account of interactions of particles with dG(κ)(t) J(κ)(t)= . (5) the image of other particles. By adding the two extra dt terms to the original Helfand moment, we recover the relation G= dtJ. Thanks to Eq. (2), the thermal conductivity coefficient is manifestly a non-negative quantity, as required by the R second law of thermodynamics. IV. NUMERICAL RESULTS We carried out molecular dynamics simulation to ap- III. HELFAND-MOMENT METHOD ply our Helfand-moment method to the calculation of thermal conductivity in a fluid with the standard 6-12 Since the transport coefficients are bulk properties, Lennard-Jones potential u(r) = 4ǫ (σ/r)12−(σ/r)6 . they can be calculated with equilibrium molecular dy- All calculations we performed are done with the cut- namics with periodic boundary conditions. This dynam- off r = 2.5σ. In the reduced un(cid:2)its with the tim(cid:3)e c ics is ruled by Newton’s equations modified in order to t∗ = t ǫ/(mσ2) and the space r∗ = r/σ, the reduced take into account the periodicity: thermal conductivity is given by p ddrta = pma + ∆r(as) δ(t−ts), κ∗ =κ σ2 m. (10) Xs kBr ǫ dp ∂u(r) a = F(r ); F(r)=− , (6) The details of the simulations are given in the preceding dt ab ∂r bX(6=a) paper [1]. To show the validity of the Helfand-moment method, we first compare with the results of the Green- where ∆r(s) is the jump of the particle a at time t with Kubo formula. We depict in Fig. 1 the time derivative a s k∆r(as)k=Lequalto thelengthLofthe simulationbox. o(3f)thaendmtehaen-tsiqmuearinetdeigsrpallacoefmtheentauoftotchoerHreelalftainondfmunomcteionnt The force F(r) is here assumed of finite range r < L/2. of the microscopic flux. The calculation is carried out The relative position appearing in the force is defined as for N =1372atoms at the phase point near of the triple pointwithareducedtemperatureofT∗ =k T/ǫ=0.722 r =r −r −L (7) B ab a b a|b and a density n∗ =nσ3 =0.8442. As seen in Fig. 1, the two methods are in perfect agreement. We estimated wherethetime-dependentcelltranslationvectorL [18] a|b thermal conductivity by a linear fit on the mean-square ischosenfortheminimum-imageconventionkr k<L/2 ab displacement of the Helfand moment. The fit is done in to be satisfied. For more details, see Section III of the region between 2 and 8 time units to guarantee that the companion paper [1]. We propose here a Helfand- the linear regime is reached. We calculated the thermal momentmethodforthermalconductivityinsystemswith conductivity for increasing system sizes from N = 108 periodic boundary conditions. As for shear viscosity [1], atoms to N = 1372. We give the results in Table I and Helfand’s originalexpression(3) must be modified in or- depict them in Fig. 2 as a function of the inverse N−1 der to take into account the periodicity of the system by of the system size. The linear extrapolation gives the the addition of extra terms: followingestimate ofthermalconductivity for aninfinite N system, G(κ)(t)= x (E −hE i)+I(t). (8) a a a κ∗ =6.990±0.030. (11) a=1 X 3 8 is based on generalized Einstein relations (2) with the Helfand moment defined to take into account periodic 7 boundary conditions in the molecular dynamics. This allows us to calculate the thermal conductivity coeffi- 6 cient in terms of the variance of the Helfand moment (9). WeshowthattheoriginalHelfandmoment(3)must * 5 k be modified by adding two extra terms. After the addi- ctivity 4 tbieotnwoeefnthtehsee emxitcrraostceormpisc,flwuexreJcoavnedr tthheerceolartrieosnpoJnd=inGg˙ u ond 3 Helfand moment G. This is the generalization of the c case of diffusion for which the velocity v (the micro- x scopic flux) is relatedto the position x (the Helfand mo- 2 ment) in the same way v = x˙. Consequently, we have x 1 Green-Kubo showed that it is possible to calculate the thermal con- Helfand ductivity coefficient by the mean-square displacement of 0 the Helfand moment. Indeed, our molecular dynamics 0 0.5 1 1.5 2 2.5 3 3.5 4 simulationsshowperfectagreementbetweentheHelfand- time moment and Green-Kubo methods. We think that the FIG. 1: Thermal conductivity at the phase point T∗ =0.722 new expressions for the Helfand moment given in the andn∗ =0.8442forN =1372. Theplainlineisthederivative presentandthecompanionpapers[1]solvethe problems ofthemean-squaredisplacement oftheHelfand moment and on the use of generalized Einstein relations in periodic thecirclestheintegralofthemicroscopicfluxautocorrelation systems. In addition to the viscosities [1] and thermal function. conductivity, it is possible to use a similar method for the electric conductivity σ in periodic systems. Indeed, by a similar derivation to the one found in Appendix A, We see in Table I that the values of the present work one can show that the electric conductivity can be writ- compare very well with the values of previous studies. ten as: 7.2 1 2 σ = lim G(σ)(t)−hG(σ)(t)i , (12) N,V,t→∞2kBTVt(cid:28)h i (cid:29) 7.1 wherethemodifiedHelfandmomentforperiodicsystems is expressed as: 7 * k y 6.9 G(σ)(t)= eZ x − ∆x(s)θ(t−t ) , (13) ctivit Xa a" a Xs a s # du 6.8 on with the electric charge eZa of the particles and c G(σ)(0)=0. Finally,wenoticethattheHelfand-moment 6.7 method plays an important role in the escape-rate for- malism [20, 21] and the hydrodynamic-mode method 6.6 [26]developedinnonequibriumstatisticalmechanics[22]. They establish relationships between microscopic and 6.5 macroscopiclevelsinthecontextoftheunderstandingof 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 1/N the origin of the irreversibility of chemical-physical phe- nomena. ThederivationoftheHelfand-momentmethods FIG. 2: Thermal conductivity at the phase point T∗ =0.722 give one the possibility to confirm numerically the theo- and n∗ = 0.8442 as a function of the inverse of the num- retical predictions of these formalisms. ber N of atoms. The circles are the results of the numerical simulations and the dashed line thelinear extrapolation. Acknowledgments We thank K. Meier for useful discussions. This re- V. CONCLUSIONS search is financially supported by the “Communaut´e franc¸aise de Belgique” (contract “Actions de Recherche In this paper, a new method is proposed for the cal- Concert´ees” No. 04/09-312) and the National Fund culation of the thermal conductivity in soft-sphere sys- for Scientific Research (F. N. R. S. Belgium, contract tems. ThetechniquewecalltheHelfand-moment method F. R. F. C. No. 2.4577.04). 4 APPENDIX A: DERIVATION OF THE HELFAND + ∆x(s)(E −hE i)δ(t−t ) MOMENT FOR THE THERMAL a a a s a s CONDUCTIVITY IN PERIODIC SYSTEMS XX 1 p +p + x a b ·F(r ) ab ab 4 m The time derivative of the modified Helfand moment a6=b X (8) is given by 1 p +p + L a b ·F(r ) 4 b|ax m ab dG(κ)(t) dxa Xa6=b = (E −hE i) dt dt a a dI(t) a + . (A4) X dt dE dI(t) a + x + (A1) a dt dt a X Onthe other hand,itis wellknownthat the microscopic where the time derivative of E is flux J(κ) for thermal conductivity is given by a dE p dp 1 ∂u p −p a a a ab a b = · + · dt m dt 2bX(6=a) ∂rab m J(κ)(t) = pmax (Ea−hEai) = 1 pa+pb ·F(rab). (A2) + X1a x pa+pb ·F(r ). (A5) 2 m ab ab 4 m bX(6=a) a6=b X We notice that there is here no jump in position to con- sider because drab concerns a relative position r = ra −rb −Lb|a wdhtich satisfies the minimum imageacbon- Since dG(dκt)(t) = J(κ)(t) by definition, we thus obtain vention within the range of the force. Symmetrizing the that: second term in Eq. (A1), we find xaddEta = 41 xapam+pb ·F(rab) dId(tt) = − ∆x(as)(Ea−hEai)δ(t−ts) Xa Xa6=b 1 Xa Xs p +p + 1 xbpa+pb ·F(rba) − 4 Lb|ax am b ·F(rab). (A6) 4 m a6=b a6=b X X 1 p +p = x a b ·F(r ) ab ab 4 m Finally, the quantity to be added to the usual Helfand a6=b X moment (3) is: 1 p +p + L a b ·F(r ) (A3) 4 b|ax m ab a6=b X I(t) = − ∆x(s)(E −hE i)θ(t−t ) where F(rab) = −F(rba) according to Newton’s third a a a s law andbecause xab =xa−xb−Lb|ax. Substituting this Xa Xs Nexepwrteossni’osneqinuaEtqio.ns(Ao1f)m, owthioenre(6dd)xt,awies ggeivtetnhebyexmproedsisfiioedn − 41 a6=bZ0tdτ Lb|axpam+pb ·F(rab). (A7) X dG(κ)(t) p ax = (E −hE i) a a dt m a X [1] S. Viscardy, J. Servantie, and P. Gaspard, companion [9] W.T.Ashurst,paperpresentedatthe13thInternational paperto appear (2007). Conference on Thermal Conductivity, Lake Orzak, Mis- [2] E. Helfand, Phys. Rev. 119, 1 (1960). souri,Nov5-7,1973,inAdvances inThermal Conductiv- [3] B.J.Alder,D.M.Gass,andT.E.Wainwright,J.Chem. ity (UniversityofMissouri-Rolla, Rolla, Missouri, 1976), Phys.53, 3813 (1970). p. 89. [4] D. Levesque, L. Verlet, and J. Ku¨rkijarvi, Phys. Rev. A [10] D. Evans, Phys.Lett. A 91, 457 (1982). 7, 1690 (1973). [11] G. Ciccotti and A. Tenenbaum, J. Stat. Phys. 23, 767 [5] M. S. Green, J. Chem. Phys. 19, 1036 (1951). (1980). [6] M. S. Green, Phys.Rev. 119, 829 (1960). [12] C. Massobrio and G. Ciccotti, Phys. Rev. A 30, 3191 [7] R.Kubo, J. Phys. Soc. Jpn. 12, 570 (1957). (1984). [8] H.Mori, Phys. Rev. 112, 1829 (1958). [13] D. Evans, Phys.Rev. A 34, 1449 (1986). 5 [14] D.M. Heyes,Phys. Rev.B 37, 5677 (1988). (1995). [15] F. Mller-Plathe, J. Chem. Phys. 106, 6082 (1997). [22] P. Gaspard, Chaos, Scattering and Statistical Mechanics [16] C. Hoheisel, Comput. Phys.Rep. 12, 29 (1990). (Cambridge University Press, Cambridge, 1998). [17] M. P. Allen, in Computer Simulation in Chemical [23] J. R. Dorfman, An Introduction to Chaos in Nonequilib- Physics, edited by M. P. Allen and D. J. Tildesley riumStatistical Mechanics (CambridgeUniversityPress, (Kluwer, Amsterdam,1993), pp.49–92. Cambridge, 1999). [18] J.M.Haile,MolecularDynamicsSimulation(JohnWiley [24] P.Gaspard andF.Baras, Phys.Rev.E51, 5332(1995). & Sons, New York,1997). [25] S. Viscardy and P. Gaspard, Phys. Rev. E 68, 041205 [19] K. Meier, Ph.D. thesis, Department of Mechanical En- (2003). gineering, University of the Federal Armed Forces Ham- [26] P. Gaspard, Phys.Rev. E 53, 4379 (1996). burg(2002). [27] G.V.Paolini,G.Ciccotti, andC.Massobrio, Phys.Rev. [20] J. R. Dorfman and P. Gaspard, Phys. Rev. E 51, 28 A 34, 1355 (1986). (1995). [21] P. Gaspard and J. R. Dorfman, Phys. Rev. E 52, 3525 6 Authors year method r∗ T∗ N κ∗ ∆κ∗ cut Levesqueet al. [4] 1973 GK N.C. 0.715 256 7.067 0.416 0.722 864 7.136 0.277 Evans[10] 1982 EM 2.5 0.722 108 6.54 N. C. Massobrio and Ciccotti [12] 1984 GK 3.5 0.721 256 6.880 0.485 DM 3.5 0.721 256 6.873 0.229 Paolini et al. [27] 1986 NEM 2.5 0.718-0.721 averagea) 6.776 0.166 Heyes[14] 1988 GK N.C. 0.72 108 6.7 0.34 (5%) 256 6.9 0.35 (5%) 500 6.5 0.33 (5%) This work 2007 HM (MD) 2.5 0.722 108 6.663 0.10 256 6.867 0.12 500 6.932 0.10 864 6.948 0.07 1372 6.946 0.12 ∞ 6.990 0.061 TABLEI:Resultsfoundintheliteratureforthethermalcon- ductivity in the Lennard-Jones fluid near the triple point. The reduced density is n∗ = 0.8442 except for Heyes [14] (n∗ = 0.848). Abbreviations: DM, nonequilibrium results obtained by the differential method. EM, Evans method [10]. GK, Green-Kubo formula with equilibrium molecu- lar dynamics. NEM, results by nonequilibrium method ob- tained by Paolini, Ciccotti and Massobrio [27]. HM, the present Helfand-moment method with molecular dynamics (MD). N.C. means thevaluehas not been communicated. a) Thermalconductivityobtainedbyaveragingthevaluesfor N =108,256,500,864.