# Methodology for determining the electronic thermal conductivity of metals via direct non-equilibrium ab initio molecular dynamics

Many physical properties of metals can be understood in terms of the free electron model kittel (), as proven by the Wiedemann-Franz law William (). According to this model, electronic thermal conductivity () can be inferred from the Boltzmann transport equation (BTE). However, the BTE does not perform well for some complex metals, such as Cu. Moreover, the BTE cannot clearly describe the origin of the thermal energy carried by electrons or how this energy is transported in metals. The charge distribution of conduction electrons in metals is known to reflect the electrostatic potential (EP) of the ion cores kittel (). Based on this premise, we develop a new methodology for evaluating by combining the free electron model and non-equilibrium ab initio molecular dynamics (NEAIMD) simulations. We demonstrate that the kinetic energy of thermally excited electrons originates from the energy of the spatial electrostatic potential oscillation (EPO), which is induced by the thermal motion of ion cores. This method directly predicts the of pure metals with a high degree of accuracy.

is one of the most important physical properties of metals. The analytical solution of based on the BTE and free electron model can be expressed as

(1) |

where is the concentration of free electrons, is the electron mass, is the Boltzmann constant, is the system temperature and is the collision time of free electrons, which is mainly determined by electron-electron, electron-hole and electron-phonon scattering. In principle, we can obtain an approximate value for from Matthiessen’s rule. However, describing every scattering process involved in the heat transfer by electrons of solid metals is too complicated. Recently, several methods have been reported for the evaluation of the electronic thermal conductivities of liquid-phase metals within the framework of density functional theory (DFT), such as ab initio molecular dynamics (AIMD) using the Kubo-Greenwood equation Collins (); Pozzo (); Nico () and DFT plus dynamical mean-field theory (DFT+DMFT) Zhang ().

In this paper, we develop a new methodology to describe the electronic heat-transport process in solid metals without explicitly addressing detailed scattering processes. From the second law of thermodynamics, we know that heat transfer in solids is driven by the temperature gradient . Typically, temperature describes the thermal motion of atoms, and the vibrations of ions can lead to spatial EPO, as can be easily deduced from the mathematical expression for the total Hamiltonian of system. Local variation of the EP can drive the collective oscillation of free electrons, and those free electrons near the Fermi surface can be excited above the Fermi surface and obtain additional thermal kinetic energy with respect to . These are called thermally excited electrons. Figs. 1(a, b) show two cartoons describing how the thermally excited electrons move in the vibrational lattice and the local EP field. Higher temperatures, which induce larger and faster ionic vibrations, lead to stronger EPO. Thus, the thermally excited electrons in high-temperature regions have more kinetic energy than those in low-temperature regions. Once a stable distribution of the thermal kinetic energy of thermally excited electrons is established along the direction of , then the heat flux carried by thermally excited electrons and can be calculated.

To prove this conjecture and quantify , we performed NEAIMD simulations Stephen (); kaviany () by modifying the Vienna Ab initio Simulation Package (VASP) vasp01 (); vasp02 (). The atomic heat flux was realized using the Müller-Plathe algorithm Muller-Plathe (), in which the kinetic energies of the atoms in the heat source and heat sink are exchanged (Supplementary Information (SI) Sec. 1). With sufficient simulation time, we can establish a stable temperature gradient in metals. Figs. 1(c, d) present the Cu model and the corresponding temperature profile, respectively. Simultaneously, we can calculate the spatial distribution and the dynamical evolution of the EP, which is expressed as

(2) |

where the test charge is the norm , and represents the ion position. Fig. 1(e) shows the theoretical results of the static distribution of the EP for a perfect Cu lattice. In the rest of this paper, we 1) illustrate the relationship between the spatial EPO and lattice vibrations, 2) demonstrate that the EPO provides additional kinetic energy to thermally excited electrons, and 3) show how to predict within our theoretical framework.

To demonstrate the relationship between EPO and lattice vibrations, we analyse the data from our AIMD simulations using the power spectral density (PSD) method Dennis (); Storch (). For a stationary signal , the PSD is defined as

(3) |

where is the autocorrelation function of Dennis (); Storch (), and denotes the expectation value. Here, we consider four signals from an AIMD simulation: atomic displacement , atomic velocity , EP displacement , and velocity of EPO (VEPO) ; these are used to calculate their respective spectral densities , Xiaoliang (), , (SI Sec. 2). and reflect the frequency-dependent lattice vibrations at a specific . Analogously, and provide information regarding the EPO with respect to frequency. We show results for Al from a 10-ps equilibrium AIMD run at 100.90 K (Figs. 2(a, b)) and a 70-ps NEAIMD simulation at 299.46 K (Figs. 2(c, d)). Fig. 2(a) clearly shows that the locations of the density peaks of and are consistent, demonstrating that the EPO is mainly caused by the lattice vibration of ion cores. Fig. 2(b) confirms this relationship. Similar results are shown in Figs. 2(c, d): for most of the frequency ranges, the peaks of and , and are consistent with each other. However, for some specific frequencies in Figs. 2(c, d) ( THz), some discrepancies exist in the peaks’ magnitudes. This phenomenon may be attributable to the heat flux applied in the NEAIMD simulations. Nevertheless, Figs. 2() provide unambiguous evidence that the EPO is directly induced by the lattice vibration of ions in metals.

To understand the dynamical evolution of spatial EPO intuitively, we present the representative case of Cu, calculated using NEAIMD at 298.49 K, in Fig. 2(e). Little variation occurs in the local electronic field between neighbouring atom layers, and the directions of these local fields continually change with time. The variations of these local fields will drive the collective vibration of free electrons, as theoretically illustrated in Fig. 2(f). We see that only free electrons near the Fermi surface can be thermally excited. Because the direction of the local field continually changes with time, the vectors of the local momentum of the thermally excited electrons should also continually change with time. Therefore, for a sufficiently long statistical time average, no net electric current should arise during the thermal transport process of metals. This is consistent with the traditional free electron model kittel ().

To prove that the EPO provides additional kinetic energy to thermally excited electrons in metals, we run a 100-ps equilibrium AIMD for Al at 329.40 K and Li at 283.97 K (both with a 222 conventional-cell and 32 total atoms). When , the total energy of the free electron system can be written as kittel (); William ()(SI Sec. 3)

where is the total energy of the free electron system at , is the thermally excited energy of the free electron system obtained from the outside environment when , is the total number of free electrons, and is the Fermi energy at . We also calculate the energy provided by EPO using

where is the average effective EPO amplitude. For Al, and . For Li, and . From these results, we have

(4) |

We conclude that lattice vibrations cause EPO in metals, and simultaneously, EPO drives the collective vibration of free electrons. The energy of these collective vibrations provides additional kinetic energy to the thermally excited electrons. This is the core concept underlying this methodology.

Within this theoretical framework, higher temperatures strengthen the spatial EPO. To prove this relationship, we perform direct FFT of the relative displacement of EP and VEPO . and were used to calculate and in Figs. 2(). describes the strength of the EPO in space, whereas reflects how fast the oscillation changes. Figs. 2(g, h) show the frequency-dependent FFT amplitudes of and , respectively. Clearly, the EPO is stronger and faster at higher temperatures.

In Fig. 2(i), we present the positive and negative integrations of the total in the same atom layer with simulation time, which can be written as , where is the index of the atom in the layer and is the total number of atoms per layer. The four quantities in Fig. 2(i) show perfect linear behaviour over time. From the absolute values, we have:

which is consistent with the evidence shown in Figs. 2(g, h). Notably, in the same temperature region, the positive and negative accumulations of are almost the same. In other words, , and thus, there is no net electric field gradient along the heat flux direction for a sufficiently long statistical time. This result confirms the physical picture illustrated in Fig. 2(f). Fig. 2(j) presents the distribution of the average effective amplitude of EPO in each atom layer along the heat flux direction, , where is the index of the atom layers. Moreover, the amplitude distribution of EPO explains how the thermal kinetic energy of thermally excited electrons is divided in space. We calculate using the RMS method John ():

where is the total number of simulation time steps, is the value of atom in a specific layer at time step , and is the average value of . Then, we define the heat flux of electrons according to the kinetic energy of thermally excited electrons between two adjacent atom layers. Because of the isotropy of the free electron model (SI Sec. 4), we take half of the difference of the thermal kinetic energy of thermally excited electrons between the two layers as

(5) | |||||

where is the cross-sectional area, is the total simulation time, is the number of free electrons per atom layer, and is the gradient of the average effective amplitude value of EPO by linear fitting of with the atom layer index number shown in Fig. 2(j). Here, a nonlinear phenomenon exists in the effective EPO amplitude distribution along the heat flux direction in some metals, such as Al, Be, and Mg. According to a case study of Be, we find that the nonlinear effect of can be reduced by increasing the system size (SI Sec. 4.2). Because of the nonlinear effect, when we calculate the of Al, Be and Mg, we fit the linear part only. For Cu and Li, the distributions exhibit perfect linear behaviour along the heat current direction. Thus, we can calculate based on Fourier’s law:

(6) |

Combining Eqs. (5, 6), we obtain the expression for :

(7) |

where is obtained by linear fitting the temperature profile with the representative case shown in Fig. 1(d).

Within this framework, we studied the of five metals (Li, Be, Mg, Al, and Cu) near room temperature. Additionally, by integrating the Müller-Plathe Muller-Plathe () atomic kinetic energy flux, as shown in Fig. 3(b), we predict the lattice (phonon) thermal conductivities of the metals (). Because of finite size effects, our NEAIMD results underestimate , especially for Cu and Mg (SI Sec. 5). By summing and from the NEAIMD simulations, we obtain the total thermal conductivities of the metals, as presented in Fig. 3(c). The results demonstrate that the thermal conductivities of metals slowly decrease with temperature near room temperature, which is consistent with traditional theory and experimental data. The error estimates in Fig. 3(c) are calculated from the expression for and error propagation theory Ku (). They mainly stem from the calculation of the gradient of and . Here, we note that because the statistical temperature fluctuation Landau () of each atom layer is large (because of the small number of atoms per layer), the conventional error estimate of will be quite large. However, NEAIMD consistently yields a stable temperature profile after a sufficiently long simulation time. Thus, we adopt the error in the linear fitting for . We also note that the aforementioned nonlinear phenomenon of the gradient of can also lead to large error bars. The details of the error-bar analysis can be found in SI Sec. 6.

In Fig. 3(d), our results show that indeed dominates the thermal transport process in metals. To compare our results with those of the traditional BTE method, we also utilize the BoltzTraP software Madsen () (based on electron energy band theory) to calculate . In obtaining , we use the constant relaxation time approximation Madsen (); Chen (). To avoid finite size effects in the calculation of , we also evaluate from the BTE method with interatomic force constants obtained from ab initio calculations broido (); chengang (), as implemented in the ShengBTE package Wu (). Then, we obtain the total thermal conductivities of the metals via the BTE method by summing from BoltzTraP and from ShengBTE. Our NEAIMD method, the traditional BTE method and experimental data are compared in Fig. 3(e). The results demonstrate that our method is superior to the traditional BTE method in predicting the electronic thermal conductivities of metals, especially for Be and Cu at room temperature. Moreover, we observe an interesting phenomenon when calculating the spectral density of VEPO () in Fig. 2(b, d). We perform exponential decay fitting of the autocorrelation function of the VEPO using the formula . Surprisingly, the exponential autocorrelation time of VEPO at room temperature is on the same approximate order of magnitude as the theoretical collision time of the free electrons Madsen (); Chen (); Campillo (). The results for Cu and Al are shown in Fig. 3(a). We also examine other metals (Be, Li, and Mg) and obtain similar results (SI Sec. 8). Therefore, we anticipate that some physical mechanisms must drive this phenomenon, i.e., it is not a coincidence.

In summary, we have developed a new methodology based on the concept of EPO to predict the electronic thermal conductivities of metals via direct NEAIMD simulation. Without explicitly addressing any complicated scattering processes of free electrons, our NEAIMD-EPO method provides better predictions of the electronic thermal conductivities of pure metals than the traditional BTE method near room temperature. We expect that this methodology will be helpful and useful for understanding and studying the heat-transfer problems of metal systems in the future. Further extension to cope with some presently challenging problems in materials, such as electron-phonon coupling, is also foreseen.

## Acknowledgments

M.H. gratefully acknowledges Prof. David G. Cahill (University of Illinois at Urbana-Champaign) for valuable comments on the manuscript. S.Y.Y. gratefully thanks Dr. Yang Han (RWTH Aachen University) and Dr. Tao Ouyang (Xiangtan University) for their helpful and fruitful discussions, and Dr. Shi-Ju Ran (ICFO, Spain) for his help in plotting Figs. 1(b, e). X.Z. greatly acknowledges Dr. Zhiwei Cui (Northwestern University) for providing the input script and MEAM potential files for the EMD simulations of Li. The authors gratefully acknowledge the computing time granted by the John von Neumann Institute for Computing (NIC) and provided on the supercomputer JURECA at Jülich Supercomputing Centre (JSC) (Project ID: JHPC25). The classical NEMD and EMD simulations were performed with computing resources granted by the Jülich Aachen Research Alliance-High Performance Computing (JARA-HPC) from RWTH Aachen University under project No. jara0135. S.S. was supported by NERC grant number NE/K006290/1.

## Author Contributions

M.H. and S.Y.Y conceived the project. M.H. supervised the project and E.D.N. co-supervised the project. S.Y.Y. developed the methodology and theoretical formula. S.Y.Y designed and executed all the AIMD simulations and calculations of the electronic thermal transport properties. G.Q., X.Z. and S.Y.Y. post-processed and analyzed the simulation data. X.Z. and S.Y.Y. performed the frequency domain analysis. S.S. provided the original code of NEAIMD. X.Z. and S.Y.Y improved the NEAIMD code. S.Y.Y. calculated the thermal conductivity of metals with the BTE method. S.Y.Y, X.Z., E.D.N., S.S and M.H. co-wrote the manuscript. All authors participated in the discussions and reviewed and revised the manuscript.

## Additional Information

Competing financial interests: The authors declare no competing financial interests.

## References

- (1) Kittel, C. Introduction to Solid States Physics (8th Ed.). Ch.6, (John.Wiley Sons, Inc., USA, 2004).
- (2) Jones, W. March, N. H. Theoretical Solid State Physics. (Courier Dover Publications, 1985).
- (3) Desjarlais, M. P., Kress, J. D. Collins, L. A. Electrical conductivity for warm, dense aluminum plasmas and liquids. Phys. Rev. E 66, 025401(R) (2002).
- (4) Pozzo, M., Davies, C., Gubbins, D. Alfè, D. Thermal and electrical conductivity of iron at Earth’s core conditions. Nature 485, 355-358 (2012).
- (5) de Koker, N., Steinle-Neumann, G. Vlek, V. Electrical resistivity and thermal conductivity of liquid Fe alloys at high P and T, and heat flux in Earth’s core. Proc. Natl. Acad. Sci. USA 109, 4070-4073 (2012).
- (6) Zhang, P., Cohen, R. E. Haule, K. Effects of electron correlations on transport properties of iron at Earth’s core conditions. Nature 517, 605-607 (2015).
- (7) Stackhouse, S., Stixrude, L. Karki, B. B. Thermal Conductivity of Periclase (MgO) from First Principles. Phys. Rev. Lett. 104, 208501 (2010).
- (8) Kim, H., Kim, M. H. Kaviany, M. Lattice thermal conductivity of using ab-initio and classical molecular dynamics. J. Appl. Phys. 115, 123510 (2014).
- (9) Kresse, G. Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B. 54, 11169-11186 (1996).
- (10) Kresse, G. Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comput. Mater. Sci. 6, 15-50 (1996).
- (11) Müller-Plathe, F. A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity. J. Chem. Phys. 106, 6082-6085 (1997).
- (12) Dennis, W. R. Echo Signal Processing. (Springer, 2003).
- (13) Storch, H. V. Zwiers, F. W. Statistical Analysis in Climate Research. (Cambridge University Press, 2001).
- (14) Zhang, X. Jiang J. Thermal conductivity of zeolitic imidazolate framework-8: a molecular simulation study. J. Phys. Chem. C. 117, 18441-18447 (2013).
- (15) Taylor, J. R. An Introduction to Error Analysis: the Study of Uncertainties in Physical Measurements (2nd Ed.). Ch.4, (University Science Books, 1997).
- (16) Ku, H. H. Notes on the use of propagation of error formulas. Journal of Research of the National Bureau of Standards - C. Engineering and Instrumentation, 70C, No. 4 (1966).
- (17) Landau, L. D. Lifshitz, E. M. Statistical Physics, Part 1 (3rd Ed.). Ch.2, (Pergamon Press, 1980).
- (18) Madsen, G. K. H. Singh, D. J. BoltzTraP. A code for calculating band-structure dependent quantities. Comput. Phys. Commun. 175, 67-71 (2006).
- (19) Chen, M. X. Podloucky, R. Electronic thermal conductivity as derived by density function theory. Phys. Rev. B 88, 045134 (2013).
- (20) Broido, D. A., Malorny, M., Birner, G., Mingo, N. Stewart, D. A. Intrinsic lattice thermal conductivity of semiconductors from first principles. Appl. Phys. Lett. 91, 231922 (2007).
- (21) Lee, S., Esfarjani, K., Mendoza, J., Dresselhaus, S. Chen, G. Lattice thermal conductivity of Bi, Sb, and Bi-Sb alloy from first principles. Phys. Rev. B. 89, 085206 (2014).
- (22) Li, W., Carrete, J., Katcho, N. A. Mingo, N. ShengBTE: A solver of the Boltzmann transport equation for phonons. Comput. Phys. Commun. 185, 1747-1758 (2014).
- (23) Campillo, I., Pitarke, J. M., Rubio, A., Zarate, E., Echenique, P. M. Inelastic lifetimes of hot electrons in real metals. Phys. Rev. Lett. 83, 2230-2233 (1999).

See pages ,,1,,2,,3,,4,,5,,6,,7,,8,,9,,10,,11,,12,,13,,14,,15,,16,,17,,18,,19, of SI.pdf