Methodology for determining the electronic thermal conductivity of metals via direct nonequilibrium ab initio molecular dynamics

Many physical properties of metals can be understood in terms of the free electron model, as proven by the Wiedemann-Franz law. 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 of the ion cores. Based on this premise, we develop a methodology for evaluating electronic thermal conductivity of metals by combining the free electron model and nonequilibrium ab initio molecular dynamics simulations. We confirm that the kinetic energy of thermally excited electrons originates from the energy of the spatial electrostatic potential oscillation, which is induced by the thermal motion of ion cores. This method directly predicts the electronic thermal conductivity of pure metals with a high degree of accuracy, without explicitly addressing any complicated scattering processes of free electrons. Our methodology offers a route to understand the physics of heat transfer by electrons at the atomistic level. The methodology can be further extended to the study of similar electron-involved problems in materials, such as electron-phonon coupling, which is underway currently.

Figure 1

Cartoons of free electrons in a metal moving in (a) the vibrating lattice and (b) electrostatic potential field.

Figure 2

Overview of the simulation model, temperature profile, and EP field of copper. (a) Model of copper used in NEAIMD simulations and (b) the corresponding temperature profile. One unit cell length comprises two layers of atoms. We use fixed boundary conditions with the layers of fixed atoms and vacuum layers along the direction of $\mathbf{\nabla }T$. Periodic boundary conditions are adopted in the other two dimensions. (c) Theoretical EP field of a perfect copper structure (the test charge number is norm 1).

Figure 3

Overview of the relationship between EPO at ion cores and lattice vibrations. Blue and red lines represent the spectral density of (a) atom displacements $\left(S_D\right)$ and (c) the EP displacement $\left(S_U\right)$, respectively, at specific ion cores. Green and orange lines represent the spectral density of (b) atom velocities $\left(S_V\right)$ and (d) the EPO velocities $\left(S_{\mathrm{\Delta }U}\right)$, respectively, at specific ion cores. Data in (a) and (b) are from a 10-ps equilibrium AIMD simulation of Al at 100.90 K, whereas (c) and (d) present the same physical quantities from a 70-ps NEAIMD simulation of Al at 299.46 K.

Figure 4

(a) The variation of EPO in space over time (the test charge number is 1). Data shown are from a 20-ps NEAIMD simulation of Cu at 298.49 K. (b) Schematic of the whole Fermi sphere oscillation as local electric field vibration along the $\vec{z}$ direction.

Figure 5

The direct fast Fourier transform (FFT) amplitudes of (a) the displacement of EP (DEP, $U_{\mathrm{ion}}$) and (b) VEPO $\left(\mathrm{\Delta }{U}_{\mathrm{ion}}\right)$ of ion cores at different temperatures. Data shown are from a 20-ps NEAIMD simulation of Cu at 298.49 K. The high and low temperatures correspond to 374.92 K and 198.39 K, respectively.

Figure 6

(a) Integration of positive and negative VEPO $\left[{\sum }_{j=1}^{N_{al}}{\sum }_t\mathrm{\Delta }{U}_j\left(t\right)\right]$ at different temperatures. (b) The effective amplitude of EPO $[\frac{1}{N_{al}}{\sum }_{j=1}^{N_{al}}\sqrt{\frac{1}{n_{\mathrm{steps}}}{\sum }_{i=1}^{n_{\mathrm{steps}}}{(U_{ij}-{\overline{U}}_j)}^2}$, average root mean square (rms) [16] of EPO, where $N_{al}$ is the atom number per layer] in atom layers along the $\vec{z}$ direction. Data shown are from a 20-ps NEAIMD simulation of Cu at 298.49 K.

Figure 7

Integration of the atomic kinetic energy flux with time based on the NEAIMD (Müller-Plathe) simulations of Al, Li, Be, and Mg.

Figure 8

NEAIMD-EPO simulation results for metals. (a) Total thermal conductivities of metals from NEAIMD simulation (with error bars determined from the calculation of $\mathbf{\nabla }T$ and $\frac{\partial {\overline{U}}_{\mathrm{EPO}}\left(l\right)}{\partial N_l}$ along the heat-transport direction) and experimental data at 300 K. (b) Pie graphs showing the electronic and phononic contributions to the total thermal conductivity of Al, Li, Be, and Mg at 300 K.

Figure 9

Bar graph comparing the thermal conductivities of metals calculated using the NEAIMD method (with error bars), Boltzmann method, and experimental data at 300 K.

Figure 10

Autocorrelation function of VEPO from NEAIMD simulations for the case of (a) Cu and (b) Al. The exponential fitting follows the function $y=Aexp(-t/{\tau }_{\exp })$, where ${\tau }_{\exp }$ represents the exponential autocorrelation time of VEPO.