Theory of the electron-ion temperature relaxation rate spanning the hot solid metals and plasma phases

TL;DR

This paper develops a comprehensive quantum-mechanical theory for electron-ion energy exchange rates applicable across hot solids and plasmas, unifying and extending existing models with broad applicability.

Contribution

It derives a general expression for the electron-ion coupling factor that includes quantum, thermal, disorder, and correlation effects, unifying various limiting models.

Findings

Reduces to the standard electron-phonon coupling in hot solids.

Extends the electron-phonon model beyond harmonic approximation.

Recovers the Spitzer formula and Fermi golden rule in plasma regimes.

Abstract

We present a theory for the rate of energy exchange between electrons and ions -- also known as the electron-ion coupling factor -- in physical systems ranging from hot solid metals to plasmas, including liquid metals and warm dense matter. The paper provides the theoretical foundations of a recent work [J. Simoni and J. Daligault, Phys. Rev. Lett. {\bf 122}, 205001 (2019)], where first-principles quantum molecular dynamics calculations based on this theory were presented for representative materials and conditions. We first derive a general expression for the electron-ion coupling factor that includes self-consistently the quantum mechanical and statistical nature of electrons, the thermal and disorder effects, and the correlations between particles. We show that our theory reduces to well-known models in limiting cases. In particular, we show that it simplifies to the standard electron-phonon coupling formula in the limit of hot solids with lattice and electronic temperatures much greater than the Debye temperature, and that it extends the electron-phonon coupling formula beyond the harmonic phonon approximation. For plasmas, we show that the theory readily reduces to well-know Spitzer formula in the hot plasma limit, to the Fermi golden rule formula in the limit of weak electron-ion interactions, and to other models proposed to go beyond the latter approximation. We explain that the electron-ion coupling is particularly well adapted to averaged atom models, which offer an effective way to include non-ideal interaction effects to the standard models and at a much reduced computational cost in comparison to first-principles quantum molecular dynamics simulations.

Equations

$$v_{\mathrm{C}}(\mathbf{r})=\frac{e^2}{|\mathbf{r}|}$$

$$H$$

$$H_{\mathrm{e}}(\mathbf{R})$$

$$\frac{d{P}_n}{dt}$$

$$M{\ddot{R}}_{\alpha }$$

$$W_{nm}=2\pi \hslash |{\mathbf{d}}_{nm}\cdot \mathbf{V}{|}^2{e}^{-\frac{(E_m-E_n)|{\mathbf{d}}_{nm}{|}^2}{2M|{\mathbf{d}}_{nm}\cdot \mathbf{V}{|}^2}}\delta (E_n-E_m),$$

$$F_{\alpha }^{\mathrm{BO}}(t)=\sum _n{P}_n(t)f_{nn}^{\alpha }(t),$$

$${\gamma }_{\alpha \beta }^{[\mathbf{R}]}=-\frac{\pi \hslash }{M}\sum _{n\ne m}\frac{P_n-P_m}{E_n-E_m}{f}_{nm}^{\alpha }{f}_{mn}^{\beta }\delta (E_n-E_m),$$

$$\ll {\xi }_{\alpha }^{[\mathbf{R}]}(t)\gg$$

$$\ll {\xi }_{\alpha }^{[\mathbf{R}]}(t){\xi }_{\beta }^{[\mathbf{R}]}(t^{\prime })\gg$$

$$B_{\alpha \beta }^{[\mathbf{R}]}=\pi \hslash \sum _{n\ne m}(P_n+P_m)f_{nm}^{\alpha }{f}_{mn}^{\beta }\delta (E_n-E_m).$$

$$c_{\mathrm{i}}^0\frac{d{T}_{\mathrm{i}}}{dt}$$

$$c_{\mathrm{e}}\frac{d{T}_{\mathrm{e}}}{dt}$$

$$g(T_{\mathrm{e}},T_{\mathrm{i}})=3k_{\mathrm{B}}n_{\mathrm{i}}\bigg{\langle}\frac{1}{3N}\sum_{\alpha=1}^{3N}\gamma_{\alpha\alpha}^{[\mathbf{R}]}(T_{e},T_{i})\bigg{\rangle}\,.$$

$$\displaystyle f({\bf R},{\bf V},t)=\left\langle\delta\Big{(}{\bf R}-{\bf R}(t)\Big{)}\delta\Big{(}{\bf V}-\dot{{\bf R}}(t)\Big{)}\right\rangle\,,$$

$$\frac{\partial f}{\partial t}+\sum _{\alpha =1}^{3N}{V}_{\alpha }\frac{\partial f}{\partial R_{\alpha }}+\sum _{\alpha =1}^{3N}\frac{F_{\alpha }^{BO}}{M}\frac{\partial f}{\partial V_{\alpha }}$$

$$\frac{d{K}_{ion}}{dt}$$

$$\frac{1}{V}\frac{d}{dt}⟨E_{\mathrm{elec}}⟩=-g(T_e-T_i),$$

$$c_e\frac{d{T}_e}{dt}=-g(T_e-T_i).$$

$$P_n=\frac{e^{-E_n/k_B{T}_e}}{\mathcal{Z}},$$

$${\gamma }_{\alpha \beta }^{[\mathbf{R}]}=\frac{1}{2M{k}_B{T}_e}{\int }_{-\infty }^{\infty }{K}_{\alpha \beta }^{[\mathbf{R}]}(t)dt=\frac{1}{2M{k}_B{T}_e}{K}_{\alpha \beta }^{[\mathbf{R}]}(\omega =0)$$

$$K_{\alpha \beta }^{[\mathbf{R}]}(t)=k_B{T}_e{\int }_0^{1/k_B{T}_e}d\lambda {⟨e^{\lambda {\hat{H}}_{\mathrm{e}}^{[\mathbf{R}]}}\delta {\hat{f}}_{\beta }{e}^{-\lambda {\hat{H}}_{\mathrm{e}}^{[\mathbf{R}]}}\delta {\hat{f}}_{\alpha }(t)⟩}_{\mathrm{e}},$$

$${\gamma }_{\alpha \beta }^{[\mathbf{R}]}=\frac{1}{2M{k}_B{T}_e}\mathrm{Re}{\int }_0^{\infty }dt⟨\delta {\hat{f}}_{\alpha }(t)\delta {\hat{f}}_{\beta }(0){⟩}_{\mathrm{e}}.$$

$${\gamma }_{\alpha \beta }^{[\mathbf{R}]}=\frac{1}{M{k}_B{T}_e}$$

$$f_{\alpha =Ix}(\mathbf{r})=-{\nabla }_{R_{\alpha =Ix}}{v}_{ie}(\mathbf{r}-{\mathbf{R}}_I),$$

$${\gamma }_{\alpha \beta }^{[\mathbf{R}]}=-\frac{1}{M}$$

$$g(T_e,T_i)$$

$$=$$

$$\sum _{\mathrm{I},\mathrm{J}}{\gamma }_{\mathrm{I}x,\mathrm{J}y}^{[\mathbf{R}]}=0\text{for all }x,y,$$

$$\sum _{\alpha ,\beta }{\gamma }_{\alpha \beta }^{[\mathbf{R}]}=0.$$

$$\delta n_e(\mathbf{r},\omega )$$

Full text

I Introduction

Nonequilibrium states of matter where the constituent electrons and ions are separately strongly driven out of equilibrium are routinely created in the laboratory. Such conditions typically occur when a material is subjected to an impulsive perturbation such as caused by an intense femtosecond laser pulse Gamaly2011 , the irradiation by swift neutrons and charged projectiles Race2010 , or a strong shock wave Celliers1992 . Understanding the energy exchanges between and among the electrons and the ions that follow the excitation and drive the system towards equilibirum is a long standing problem in condensed matter Kaganov1957 ; Fujimoto1984 ; Schoenlein1987 ; Elsayed_et_al_1987 ; Allen1987 ; Fann1992 ; Groeneveld1995 ; Hohlfeld2000 ; DelFatti2000 ; Bonn2000 ; Rethfeld2002 ; Lin_et_al_2008 ; Mueller2013 ; JiZhang2016 ; Henighan2016 ; Waldecker2016 ; Sadasivam2017 ; Maldonado2017 ; Lu1 and plasma physics Spitzer1962 ; Ramazashvili1962 ; Kihara1963 ; Brysk1975 ; Hazak2001 ; DharmaWardana_2001 ; Gericke2002 ; Daligault_Mozyrsky_2007 ; Daligault_Mozyrsky_2008 ; Glosli2008 ; DimonteDaligault2008 ; BrownSingleton2009 ; DaligaultDimonte2009 ; Faussurier2010 ; Vorberger2012 ; Baalrud2014 ; Ng2012 ; Leguay2013 ; Cho_et_al_2015 ; Jourdain_et_al_2018 ; Ogitsu2018 ; SimoniDaligaultPRL2019 . Because of the small electron to ion mass ratio, it is customary to distinguish two distinct time scales, namely a short time scale that characterizes the fast internal thermalization of each particle species, and a longer time scale that characterizes the slower equilibration of the electron and ion temperatures. The latter, which is the subject of this work, is governed by the strength of the electron-ion coupling factor.

For solids and weakly coupled plasmas, the central mechanisms that govern the energy exchanges between electrons and ions have been known for a long time. In the case of solids, the energy exchanges are well described as resulting from the interactions between electrons and phonons Allen1987 ; Patterson2010 . Yet, the modeling of energy exchanges between electrons and phonons remains an active subject of research driven by increasingly accurate measurements and numerical simulationsMueller2013 ; Waldecker2016 ; Henighan2016 ; Maldonado2017 ; Lu1 . In the case of weakly coupled plasmas, the energy exchanges are well described as resulting from individual binary collisions between charged particles screened by the surrouding plasma Spitzer1962 . Different methods have been proposed and validated with simulations to self-consistently take into account the effects of the plasma on binary collisions Glosli2008 ; BrownSingleton2009 ; DimonteDaligault2008 ; DaligaultDimonte2009 ; Baalrud2014 .

For systems at the confluence of solids and plasmas, however, different models Lin_et_al_2008 ; Brysk1975 ; Hazak2001 ; DharmaWardana_2001 ; Daligault_Mozyrsky_2008 ; Faussurier2010 ; Vorberger2012 have been proposed that offer diverging predictions even for simple materials (see table I in SimoniDaligaultPRL2019 ). This intermediate regime, which is characterized by the coexistence and interplay of significant quantum, thermal, disorder and strong Coulomb interaction effects, challenges the standard simplifying approximations of either ordinary condensed-matter physics (e.g., band structure, phonons, etc.) or plasma physics (classical statistics, binary collisions, etc.) T1 . The last decade has seen remarkable progress in our ability to form and interrogate in the laboratory materials under conditions at the confluence of solids and plasmas Fortov2016 . These experiments typically produce transient, non-equilibrium conditions and measurements may be misleading if recorded while the plasma species are still out of equilibrium. Remarkably, the electron-ion energy relaxation rate is now accessible to experimental measurements thanks to the diagnostic capabilities offered by the new generation of x-ray light sources Leguay2013 ; Cho_et_al_2015 ; Jourdain_et_al_2018 ; Ogitsu2018 .

Motivated by these recent developments, we present a theory for the electron-ion coupling factor, the rate of energy exchange between electrons and ions, that applies to physical systems ranging from hot solid metals to plasmas, including liquid metals and warm dense matter. In a recent paper SimoniDaligaultPRL2019 , we have presented first-principles calculations based on this theory for representative materials of various electronic complexity and over a range of conditions, but the present theory was only briefly outlined without justification. The purpose of this paper is to provide a detail exposition of the theory, its properties and its relation to previous models. We focus on the theoretical aspects only, the numerical algorithms used to obtain the results shown in SimoniDaligaultPRL2019 will be presented in another publication. We do not present calculations on specific physical systems in this paper and refer the reader to SimoniDaligaultPRL2019 and future publications for such applications.

The paper is organized as follows.

In Sec. II, we derive the general expression for the electron-ion coupling factor $g$. The derivation relies on a recent work presented in Ref. Daligault_Mozyrsky_2018 , where it was shown, under mild assumptions suitable for the physical systems of interest here, that the coupled dynamics of ions and electrons can be replaced by a simpler, effective classical-quantum dynamics, in which the atomic motions are governed by a stochastic Langevin-like equation and the electron dynamics is described by a master equation for the populations of the electronic states. By assuming that the electron and ion systems can be characterized by separate temperatures, the theory implies that the two temperatures satisfy coupled rate equations and yields an explicit expression for the relaxation rate – the electron-ion coupling factor – in terms of the friction coefficients felt by individual ions due to their non-adiabatic interactions with the electrons. Each friction coefficient satisfies a Kubo relation given by the time integral of the autocorrelation function of the interaction force of an ion with the electrons. Several equivalent expressions are given for $g$.

In Sec. III, we recast the general result of Sec. II in terms of quantities that highlight the many-body screening and correlation effects, and that are more easily amenable to theoretical analysis, practical approximations and numerical evaluations. We consider three different reformulations, each being based on a different representation of the many-body electron response function in terms of reference response functions, namely the free-particle, the proper and the Kohn-Sham response functions commonly used in condensed matter and plasma physics. The formulation in terms of the Kohn-Sham response function is at the basis of the first-principles molecular dynamics simulations presented in SimoniDaligaultPRL2019 , which will be discussed at length elsewhere.

In Sec. IV, we show that our theory reduces to well-known models in limiting cases and suggests improved practical models. We show that it reduces to the traditional Spitzer formula in the hot plasma limit Spitzer1962 , to the Fermi golden rule formula in the limit of weak electron-ion interactions Hazak2001 , and to the model of Daligault and Dimonte in their attempt to include non-linear electron-ion effects DaligaultDimonte2009 . We then explain that the electron-ion coupling is particularly well adapted to averaged atom models, which have proved to be accurate and computationally much more expedient than quantum molecular dynamics simulations. We then demonstrate that our theory also applies to hot solids, namely to solid metals with lattice and electronic temperatures much greater than the Debye temperature. The theory reduces to the standard electron-phonon coupling formula Allen1987 in the appropriate limit and extends the latter by including ionic motions beyond the harmonic approximation. We finally relate our theory to a simple model due to Wang et al. Wang1994 , which has served as a reference in recent works on the temperature relaxation rates in hots solids and warm dense matter Lin_et_al_2008 ; Cho_et_al_2015 .

For clarity, the technical details are presented in the appendices. Throughout the paper, $\hbar$ is the reduced Planck constant, $k_{B}$ is the Boltzmann constant, and $m_{u}$ is the atomic mass unit, and $e^{2}=q_{e}^{2}/4\pi\epsilon_{0}$, where $q_{e}$ is the elementary charge and $\epsilon_{0}$ is the vacuum permittivity. $\mathrm{Re}z$ and $\mathrm{Im}z$ denote the real and imaginar parts of a complex number $z$. Throughout the paper, $t$ and $\omega$ denote the time and frequency variables and, with no risk of confusion, $f(\omega)=\int_{-\infty}^{\infty}{dte^{i\omega t}f(t)}$ denotes the Fourier transform of the function $f(t)$. Finally,

[TABLE]

denotes the Coulomb potential energy.

II Derivation and general expressions of the electron-ion coupling factor

II.1 Definitions and assumptions

We consider a material containing one atomic species enclosed in a three-dimensional cubic box of volume $V=L^{3}$. The material is described as a two-component system comprised of ions (mass $m_{i}=Am_{u}$, number density $n_{i}=N/V$, charge $Ze$) and of electrons (mass $m_{e}$, density $n_{e}=Zn_{i}$), where each ion consists of an atomic nucleus and its most tightly bound, unresponsive core electrons. Throughout the paper, the calculations are performed by imposing three-dimensional periodic (Born-von Karman) boundary conditions along each direction of the box; physical results are then obtained in the thermodynamic limit where both $N$ and $V$ tend to infinity in such a way that $n_{i}$ remains constant.

The total Hamiltonian of the system is

[TABLE]

with the electron Hamiltonian

[TABLE]

where $v_{ii}$ and $v_{ie}$ are the ion-ion and ion-electron interaction energies. For simplicity of exposition, the electron-ion interaction is described by a local pseudopotential $v_{ie}(r)$; in practice, the formalism allows to deal with more elaborate descriptions, e.g., using plane-augmented wave pseudopotentials as will be discussed elsewhere.

Below, $R_{\alpha=Ix}$ denotes the position of ion $I$ along the $x$-direction and ${\bf R}=\{R_{\alpha}\}_{\alpha=1,\dots,3N}$ denotes the set of all ionic positions.

II.2 Effective dynamics of ions and electrons

For the metallic systems of interest here, the dynamics governed by the complete quantum Hamiltonian (2) can be replaced by a simpler, effective classical-quantum dynamics by making use of the naturally small electron-ion mass ratio and the existence of a manifold of infinitesimally separated electronic excitations. More precisely, we make the following three assumptions:

    (i) The dynamics of each ion can be described by that of the center ${\bf R}_{i}(t)$ of its narrowly localized wavepacket. This is justified here, since the thermal de Broglie wavelength $\Lambda=\hbar\sqrt{2\pi/m_{i}k_{B}T_{i}}(\simeq 0.3/\sqrt{AT_{i}[eV]}$ Bohr) of ions is generally much smaller than the spatial variations of forces acting on them due to their large mass and the relatively high temperatures.

    (ii) The typical ion velocities are small compared to the typical electronic velocities. For instance, for the two-temperature systems considered later, we assume $T_{i}/m_{i}\!\ll\!T_{F}/m_{e}$ or $T_{i}/m_{i}\!\ll\!T_{e}/m_{e}$ in the degenerate $T_{e}/T_{F}\!\ll\!1$ or non-degenerate limit $T_{e}/T_{F}\!\gg\!1$, respectively, where $T_{F}=\frac{\hbar^{2}}{2m_{e}k_{B}}(3\pi^{2}n_{e})^{\frac{2}{3}}$ ($\simeq 1.69\,\left(n_{e}[{\rm cm}^{-3}]/10^{22}\right)^{\frac{2}{3}}$ eV) is the electronic Fermi temperature. This condition is generally respected due to the natural smallness of $m_{e}/m_{i}$, and is challenged only if $T_{i}\gg T_{e}$.

    (iii) We assume that there is a quasi-continuum of electronic states, as it is the case for the metallic systems of interest here.

Under these conditions, the electron and ion dynamics can be described by the coupled set of equations Daligault_Mozyrsky_2018

[TABLE]

The electron dynamics is described by the master equation (3a) for the populations $P_{n}(t)$ of the adiabatic electronic states $\ket{n(\mathbf{R}(t))}$ defined by $\hat{H}_{\mathrm{e}}(\mathbf{R}(t))\ket{n(\mathbf{R}(t))}=E_{n}(\mathbf{R}(t))\ket{n(\mathbf{R}(t))}$. The transition rates between different electronic states are given by

[TABLE]

where $\mathbf{V}=\dot{\mathbf{R}}$ represents the full set of atomic velocities and the non adiabatic couplings $\mathbf{d}_{nm}=\matrixelement{n}{\nabla_{\mathbf{R}}}{m}$.

Each ionic position follows a stochastic Langevin-like equation (3b), where

[TABLE]

is the adiabatic Born-Oppenheimer force, which includes the interactions between ions and with the instantaneous electrostatic potential of electrons, where $f_{nn}^{\alpha}(t)=\matrixelement{n}{-\nabla_{R_{\alpha}}\hat{H}_{\mathrm{e}}({\bf R}(t))}{n}$ represents the force exerted on the degree of freedom $R_{\alpha}$ by the state $\ket{n(\{\mathbf{R}(t)\})}$. The other terms describe the effect of non-adiabatic transitions between closely spaced electronic states induced by the atomic motions and electronic excitations. These terms, which are not accounted for in current quantum molecular dynamics simulations, are responsible for the constant, non-reversible, energy exchanges between electron and ions. Like the buffeting of light liquid particles on a heavy Brownian particle, the non-adiabatic effects produce the friction forces $-M\gamma_{\alpha\beta}^{[\mathbf{R}]}\dot{R}_{\beta}$, where

[TABLE]

Here, the out-of-diagonal force matrix elements are defined as $f_{nm}^{\alpha}=\matrixelement{n}{-\nabla_{R_{\alpha}}\hat{H}_{\mathrm{e}}({\bf R}(t))}{m}$. The symbol $[\mathbf{R}]$ is used to indicate that the quantity depends on the instantaneous atomic configuration of the system. However, in order to avoid cluttering the mathematical expressions, we do not always indicate the dependence on $[\mathbf{R}]$.

The second term is a $\delta$-correlated Gaussian random fluctuating force, $\xi_{\alpha}^{[\mathbf{R}]}(t)$, that satisfies the following two conditions

[TABLE]

where $\ll\ldots\gg$ denotes an average over the Gaussian noise, and

[TABLE]

II.3 Application to a two-temperature plasma

II.3.1 General result

For our present purpose, we assume throughout the rest of the paper that the material can be described as an isolated, homogeneous, two-temperature system characterized at all times $t$ by the temperatures $T_{e}(t)$ and $T_{i}(t)$ of the electronic (e) and ionic (i) subsystems. (We could in principle relax the homogeneity condition and add an external energy source, e.g. a laser, but this is not necessary here.) The instantaneous electronic populations are then $P_{n}=e^{-E_{n}/k_{B}T_{e}}/\mathcal{Z}$, with the canonical partition function $\mathcal{Z}=\mathrm{Tr}e^{-\hat{H}_{\mathrm{e}}^{[\mathbf{R}]}/k_{B}T_{e}}=\sum_{n}{e^{-E_{n}/k_{B}T_{e}}}$.

The temporal evolution of the temperatures can be obtained by applying the evolution equations (3b)-(3a) to the ensemble averaged kinetic energy of the ions $K_{\mathrm{ion}}(t)=\left\langle\frac{1}{2}M\dot{\bf R}^{2}(t)\right\rangle$ and to the internal energy of electrons E_{\mathrm{elec}}(t)=\Big{\langle}\sum_{n}{P_{n}(t)E_{n}(t)}\Big{\rangle}. As shown below, this yields the rate equations

[TABLE]

where $c_{\mathrm{i}}^{0}=V^{-1}\partial[3Nk_{B}T_{i}/2]/\partial T_{i}=3n_{i}k_{B}/2$ is the kinetic component of the heat capacity of the ions, c_{\mathrm{e}}=V^{-1}\partial\big{\langle}\sum_{n}{P_{n}E_{n}}\big{\rangle}/\partial T_{e} is the electronic heat per unit volume, and $g$ is the electron-ion coupling factor of interest here given by

[TABLE]

We show in Sec. II.4 that the friction coefficients $\gamma_{\alpha\alpha}^{[\mathbf{R}]}$ and the resulting coupling factor $g$ can be written in the form of standard Kubo relations like the ordinary electronic and ionic transport coefficients. The physical content of Eq.(10) is discussed in the following sections. Before, we present a derivation of Eqs.(9) and (10).

II.3.2 Proof of Eqs.(9) and (10).

Equation (9a) for the ionic temperature can be readily obtained by recalling the close relationship between the Langevin equation and the Fokker-Planck-Kramers equation vanKampenbook . The swarm of trajectories generated by Eq.(3b) can be described by the probability distribution function

[TABLE]

with ${\bf R},{\bf V}\in\mathbb{R}^{3N}$, and $\langle\dots\rangle$ represents the double average $\langle\ll\dots\gg\rangle$ over the noise and over the initial distribution function. The distribution function $f$ satisfies the Fokker-Planck-Kramers equation

[TABLE]

Remark that this equation, which governs the full distribution functions of the ions, should not be confused with the celebrated Fokker-Planck equation for the single-particle distribution function that is widely used in ordinary plasma physics ThomsonHubbard1960 . From this evolution equation, we find the time evolution of the kinetic energy of ions $K_{ion}(t)=\iint{\frac{1}{2}M{\bf V}^{2}f({\bf R},{\bf V},t)d{\bf R}d{\bf V}}$,

[TABLE]

By assuming a Maxwellian velocity distribution at temperature $T_{i}(t)$, i.e. $f({\bf R},{\bf V},t)\propto e^{-M{\bf V}^{2}/2k_{B}T_{i}(t)}$, the kinetic energy $K_{ion}(t)=3Nk_{B}T_{i}(t)/2$ and Eq.(II.3.2) simplifies

[TABLE]

Since we also assume a thermal distribution of adiabatic states at temperature $T_{e}$, $\gamma_{\alpha\alpha}$ and $B_{\alpha\alpha}$ satisfy the relation $B_{\alpha\alpha}=2Mk_{B}T_{e}\gamma_{\alpha\alpha}$ (compare Eqs.(5) and (8), see details in Daligault_Mozyrsky_2018 ), and Eq.(II.3.2) reduces to the desired result, Eq.(9a).

The equation (9b) for the electronic temperature is obtained by combining Eq.(II.3.2) with the conservation equation $\frac{d}{dt}\left\langle\frac{1}{2}M{\dot{\bf R}}^{2}+V_{ii}+\sum_{n}{P_{n}E_{n}}\right\rangle=0$ and with the property $d\langle V_{ii}\rangle/dt=0$ that is easily shown using the Fokker-Planck-Kramers equation. A proof of the energy conservation equation can be found in Ref. Daligault_Mozyrsky_2018 . This yields the desired result for the rate of change of the electronic energy $E_{\mathrm{elec}}$,

[TABLE]

or, with $c_{\mathrm{e}}\equiv V^{-1}\partial\langle E_{\mathrm{elec}}\rangle/\partial T_{e}$,

[TABLE]

II.4 Kubo relations for the friction coefficients and the electron-ion coupling. Sum rule.

Under the two temperature assumption considered here, the electronic populations are

[TABLE]

and the general expression (5) for the friction coefficients can be effectively written compactly in the form of ordinary Kubo relations, i.e. as time integrals of correlation functions. The different expressions below result from well-known relations between thermal correlation functions, response function and their Lehmann representations Kubo_book ; their defintions are recalled for convenience in appendix A. Below, \langle\ldots\rangle_{\mathrm{e}}=\mathrm{Tr}\bigg{[}\frac{e^{-\hat{H}_{\mathrm{e}}/k_{B}T_{e}}}{\mathcal{Z}}\ldots\bigg{]} indicates the thermal average over the electronic subsystem.

i)

With Eq.(17), the expression (5) can be written compactly as (see Eq.(LABEL:K_lehmann))

[TABLE]

in terms of the Kubo correlation function

[TABLE]

where $\hat{f}_{\alpha=Ix}(t)=-e^{i\hat{H}_{\mathrm{e}}t/\hbar}\,\partial_{R_{\alpha}}\hat{H}_{\mathrm{e}}\,e^{-i\hat{H}_{\mathrm{e}}t/\hbar}$ is the force operator at time $t$ between ion $I$ and the electronic subsystem along the $x$-direction and $\delta\hat{f}_{\alpha}(t)$ indicates the same operator deprived of its diagonal matrix elements, i.e. $(\delta f_{\alpha})_{nm}=f_{nm}^{\alpha}(1-\delta_{n,m})$.

ii)

Using the property (84), the relation (18) can be written in terms of the electron-ion force-force correlation function $\langle\,\delta\!\hat{f}_{\alpha}(t)\delta\!\hat{f}_{\beta}(0)\,\rangle_{\mathrm{e}}$,

[TABLE]

iii)

Equation (20) can be expressed in terms of the symmetric electronic density correlation function $S_{\mathrm{ee}}^{[\mathbf{R}]}(\mathbf{r}_{1},\mathbf{r}_{2},t)=\frac{1}{2}\langle\,\delta\hat{n}_{\mathrm{e}}(\mathbf{r}_{1},t)\delta\hat{n}_{\mathrm{e}}(\mathbf{r}_{2},0)+\delta\hat{n}_{\mathrm{e}}(\mathbf{r}_{2},0)\delta\hat{n}_{\mathrm{e}}(\mathbf{r}_{1},t)\,\rangle_{\mathrm{e}}$ as follows

[TABLE]

where

[TABLE]

is the force along the $x$-direction between the ion $I$ and an electron located at ${\bf r}$. Equation (LABEL:gammaalphabetaSee) is easily obtained using $\partial_{{\bf R}_{I}}\hat{H}_{\mathrm{e}}=\int_{V}{d{\bf r}\,\partial_{{\bf R}_{I}}\,v_{ie}(\mathbf{r}-{\bf R}_{I})\hat{n}_{\mathrm{e}}(\mathbf{r})}$ in Eq.(20).

iv)

The expression (18) can also be written in terms of the electron response function (a.k.a. susceptibility), $\chi_{\mathrm{ee}}^{[\mathbf{R}]}(\mathbf{r}_{1},\mathbf{r}_{2},t)=-\frac{i}{\hbar}\theta(t)\langle\,[\delta\hat{n}_{\mathrm{e}}(\mathbf{r}_{1},t),\delta\hat{n}_{\mathrm{e}}(\mathbf{r}_{2},0)]\,\rangle_{\mathrm{e}}$ of the electronic subsystem in the frozen ionic configuration ${\bf R}$, as follows

[TABLE]

This is easily found using Eq.(LABEL:gammaalphabetaSee) and the fluctuation-dissipation relation (85) between $\chi_{\mathrm{ee}}^{[\mathbf{R}]}$ and $S_{\mathrm{ee}}^{[\mathbf{R}]}$. We shall mainly rely on the expressions (20) and (LABEL:Eq:5) in the reminder of the paper.

v)

With the help of the previous expressions, the electron-ion coupling factor (10) writes as

[TABLE]

vi)

The electron-ion coupling (10) equals the trace of the matrix $\tensor{\gamma}^{[\mathbf{R}]}=\left\{\gamma^{[\mathbf{R}]}_{\alpha,\beta}\right\}$ of friction coefficients. Other combinations of matrix elements satisfy remarkable properties. Most remarkably,

[TABLE]

and, therefore,

[TABLE]

for the sum of all matrix elements. These sum rules, which are physically related to the conservation of momentum, are proved in appendix C.

III Electronic screening, exchange and correlation effects

Thus far we have given general expressions for the friction coefficients and the electron-ion coupling factor in which the electrons are not described individually but are described as a single entity. For instance, the electronic states $|n({\bf R})\rangle$ in Eq.(5) are many-body states, and the response function $\chi_{\mathrm{ee}}^{[\mathbf{R}]}$ in Eq.(LABEL:Eq:5) is the full density-density response function of the electronic subsystem. In this section, we recast these results in terms of quantities that instead emphasize the individual character of electrons. The many-body screening and correlation effects are displayed more distinctly in terms of dielectric functions and local field corrections. To accomplish this, we express the full density-density response function $\chi_{\mathrm{ee}}^{[\mathbf{R}]}$ in terms of the response of a reference system of independent particles. We shall consider three different reference response functions that are often used GiulianiVignale2005 , namely the proper response function, the non-interacting response function, and the Kohn-Sham response function. As we shall show in Sec. IV, previous models for $g$ are easily recovered from these reformulations.

For pedagogical clarity, we follow the same line of presentation in each case. To this end, we recall that the response function $\chi_{\mathrm{ee}}^{[\mathbf{R}]}$ gives the change in the ground-state electronic density $n_{e}^{(0)}$ through

[TABLE]

when the electron subsystem in the frozen ionic configuration ${\bf R}$ is perturbed by a weak time-dependent scalar potential $\delta v_{ext}({\bf r},t)$ GiulianiVignale2005 .

III.1 Relation to the proper response function

The proper density-density response function $\tilde{\chi}^{[\mathbf{R}]}$ allows to write Eq. (27) as GiulianiVignale2005

[TABLE]

in terms of the screened potential

[TABLE]

This is the potential experienced by a test particle (i.e., a fictitious particle that does not disturb the system in which it is embedded) in the electron gas due to both the external potential and the Coulomb field created by the density pertubation $\delta n_{e}$ induced by $\delta v_{ext}$. This potential does not account for the correlation that exists between a given electron of the electron gas and the other electrons. By definition, these correlation effects are implicitely incorporated in the proper response function $\tilde{\chi}^{[\mathbf{R}]}$.

The proper response function $\tilde{\chi}^{[\mathbf{R}]}$ is related to the full response function $\chi_{\rm ee}^{[\mathbf{R}]}$ through the integral (Dyson) equation,

[TABLE]

Using this relation in Eq.(LABEL:Eq:5), it is straigthforward to write the friction coefficient $\gamma_{\alpha\beta}^{[\mathbf{R}]}$ as (see appendix B.1)

[TABLE]

where the $f_{\alpha=Ix}^{\mathrm{L(R)}}$ represents the force of interaction along the $x$-direction of ion $I$ and a test charge,

[TABLE]

in terms of the inverse of the left (L) and right (R) dielectric functions

[TABLE]

Note that the definition of the left (L) and right (R) dielectric functions is needed at this level of generality since the system is embedded in the inhomogenous background of the ionic configuration ${\bf R}$ (for homogeneous systems, $\epsilon_{\mathrm{L}}=\epsilon_{\mathrm{R}}$ (see Sec. IV)).

With Eq.(LABEL:Eq:5b), the electron-ion coupling factor becomes

[TABLE]

The expressions (LABEL:Eq:5b) and (34) make a good starting point of further theoretical analysis because the proper response function – or irreducible reponse function – lends itself well to advanced perturbative methods in order to systematically include the correlation effects beyond the mean-field approximation GiulianiVignale2005 ; Mahan2000 ; BruusFlensberg2004 . This goes beyond the scope of the present work and, instead, we shall now consider another representation where these important many-body effects appear even more visibly in terms of local field corrections.

III.2 Relation to the ideal gas response function

It is also common to express the deviation $\delta n_{e}$ as the induced density of a noninteracting (free) electron gas such as GiulianiVignale2005

[TABLE]

where $\chi_{0}^{[\mathbf{R}]}$ is the density-density response function of the inhomogeneous, free-electron gas in the static ionic configuration ${\bf R}$, and $\delta v_{\mathrm{eff}}$ is the effective potential

[TABLE]

In contrast with the previous section, here, it is the free electron response function that does not include the correlations existing between a given electron and the other electrons of the gas. These correlations are included through the last term of the effective potential by means of the so-called local field correction $G_{ee}^{[\mathbf{R}]}(\mathbf{r}_{1},\mathbf{r}_{2},\omega)$ (we extend standard definitions for the homogeneous electron gas GiulianiVignale2005 ; Kugler1975 to the non-homogeous case).

The response functions $\chi_{\mathrm{ee}}^{[\mathbf{R}]}$ and $\chi_{0}^{[\mathbf{R}]}$ are related through the integral equation

[TABLE]

with the interaction kernel

[TABLE]

As shown in appendix B.2, with this relation, the friction coefficient (LABEL:Eq:5) can be written as the sum of two terms,

[TABLE]

The first term writes as

[TABLE]

the spatial derivative of the effective screened potentials are given respectively by

[TABLE]

where the electronic screening effect is described by the inverse of the (electron-test charge) dielectric functions

[TABLE]

By comparison with Eqs.(24), we see that $\bar{\gamma}_{\alpha\beta}^{[\mathbf{R}]}$ can also be expressed as follows

[TABLE]

in terms of the corrrelation function $\left\langle\,\delta{\bar{f}}_{\alpha}^{L}(t)\delta{\bar{f}}_{\alpha}^{R}(0)\,\right\rangle_{\mathrm{e}}$ of the electron-ion force screened by the electron-test charge dielectric functions $\bar{\varepsilon}_{\mathrm{L,R}}$, where ${\bar{f}}_{\alpha}^{L,R}(t)=e^{i\hat{h}_{\mathrm{e}}^{0}t/\hbar}\bar{f}_{\alpha}^{\,\mathrm{L,R}}e^{-i\hat{h}_{\mathrm{e}}^{0}t/\hbar}$ and $\hat{h}_{e}^{0}=\frac{\hat{{\bf p}}^{2}}{2m_{e}}+\sum_{I=1}^{N}{v_{ie}(\hat{\bf r}-{\bf R}_{I})}$ is the Hamiltonian of a free electron moving in the background of the ionic potential.

The second term in Eq.(39) writes as

[TABLE]

where $n_{\alpha}^{\prime}(\mathbf{r})=\int_{V}d{{\bf r}^{\prime}}f_{\alpha}(\mathbf{r^{\prime}})\chi_{\mathrm{ee}}^{[\mathbf{R}]}(\mathbf{r^{\prime}},\mathbf{r},\omega=0)$.

We expect that for many systems the second term will be a small correction to the first term. In practice, the frequency dependence of the local field correction is poorly known and the so-called static approximation $G_{ee}^{[\mathbf{R}]}(\mathbf{r}_{1},\mathbf{r}_{2},\omega)=G_{ee}^{[\mathbf{R}]}(\mathbf{r}_{1},\mathbf{r}_{2},\omega=0)$ is commonly used; in this so-called static local-field correction approximation, the last term cancels out and $\gamma_{\alpha\beta}^{[\mathbf{R}]}=\bar{\gamma}_{\alpha\beta}^{[\mathbf{R}]}$.

With Eq.(39), the electron-ion coupling factor becomes

[TABLE]

III.3 Relation to the Kohn-Sham response function

The previous expressions in terms of the proper response function or in terms of the local field corrections are useful for theoretical analysis using many-body techniques and for the developments of practical models (see, e.g., Sec. IV). Another approach consists in expressing Eq.(24) in terms of quantities directly accessible to ab-initio simulations. The formulation in terms of the Kohn-Sham response function described here is at the basis of the first-principles molecular dynamics simulations presented in SimoniDaligaultPRL2019 . We restrict ourselves to the presentation of the exact expressions; the practical algortihms used to numerically evaluate them will be discussed at length elsewhere SimoniDaligaultQMDpaper .

The most important among these methods is Density Functional Theory (DFT) Ulrich2012 , in which the density deviation Eq.(27) can be written as

[TABLE]

where the effective interaction

[TABLE]

is the sum of the external perturbation $\delta v_{\rm ext}$ and of the linearized Kohn-Sham potential. The latter consists of the sum of the linearized Hartree potential and of the linearized exchange-correlation potential, where $f_{\mathrm{xc}}^{[\mathbf{R}]}[n_{e}](\mathbf{r},\mathbf{r^{\prime}},\omega)$ is the so-called exchange-correlation kernel Ulrich2012 . The Kohn-Sham density-density response function is given by Ulrich2012

[TABLE]

where $p_{i}^{\rm eq}=1/\left[1+e^{-(\mu(T_{e})-\epsilon_{i})/k_{B}T_{e}}\right]$ is the Fermi-Dirac function at temperature $T_{e}$,and $\innerproduct{\bf r}{i}=\phi_{i}(\mathbf{r})$ are the Kohn-Sham wave functions, obtained from the solution of the following set of equations

[TABLE]

where

[TABLE]

is the Kohn-Sham potential, $v_{\mathrm{H}}[n_{\rm e}](\mathbf{r})=e^{2}\int_{V}d{{\bf r}^{\prime}}\frac{n_{\mathrm{e}}(\mathbf{r^{\prime}})}{|\mathbf{r}-\mathbf{r^{\prime}}|}$ is the Hartree potential, $v_{\mathrm{xc}}[n_{\rm e}](\mathbf{r})$ is the exchange-correlation potential, and

[TABLE]

is the ground state electron density of the physical system. The present approach differs from the one of the previous section in the fact that, unlike $\chi_{0}^{[\mathbf{R}]}$ that describes the response of a free electron gas to the potential $\delta v_{\rm{eff}}^{[{\bf R}]}({\bf r})$, the Kohn-Sham response function is the response of a system of free particles moving under the effect of the self-consistent Kohn-Sham potential $\delta v_{\mathrm{KS}}({\bf r})$, Eq.(47).

The response functions $\chi_{\mathrm{ee}}^{[\mathbf{R}]}$ and $\chi_{\mathrm{KS}}^{[\mathbf{R}]}$ are related through the integral equation Ulrich2012

[TABLE]

where the kernel $K_{\mathrm{KS}}^{[\mathbf{R}]}(\mathbf{r},\mathbf{r^{\prime}},\omega)$ is now given by the sum of the Hartree and exchange-correlation contributions

[TABLE]

Using Eq. (LABEL:DysonchichiKS) into Eq. (LABEL:Eq:5), each friction coefficient becomes the sum of two contributions,

[TABLE]

The first term is

[TABLE]

where the spatial derivative of the effective screened potentials are given respectively by

[TABLE]

and the screening effect of the electronic subsystem is described by the inverse of the dielectric functions

[TABLE]

By comparison with Eqs.(24), we see that $\tilde{\gamma}_{\alpha\beta}^{[\mathbf{R}]}$ can also be expressed as follows

[TABLE]

in terms of the corrrelation function $\left\langle\,\delta{\tilde{f}}_{\alpha}^{L}(t)\delta{\tilde{f}}_{\alpha}^{R}(0)\,\right\rangle_{\mathrm{e}}$ of the KS particle-ion force screened by the KS particle -test charge dielectric function $\tilde{\varepsilon}_{\mathrm{L,R}}$, where ${\tilde{f}}_{\alpha}^{L,R}(t)=e^{i\hat{h}_{\mathrm{KS}}t/\hbar}\tilde{f}_{\alpha}^{\,\mathrm{L,R}}e^{-i\hat{h}_{\mathrm{KS}}t/\hbar}$, and $\hat{h}_{KS}=\frac{\hat{{\bf p}}^{2}}{2m_{e}}+v_{KS}(\hat{\bf r})$ is the KS Hamiltonian.

The second term in the expression (51), instead, represents a quantum many body correction including all the memory effects

[TABLE]

where $n_{\alpha}^{\prime}(\mathbf{r})=\int_{V}d{{\bf r}^{\prime}}f_{\alpha}(\mathbf{r^{\prime}})\chi_{\mathrm{ee}}^{[\mathbf{R}]}(\mathbf{r^{\prime}},\mathbf{r},\omega=0)$. In practice SimoniDaligaultQMDpaper , the frequency dependence of the exchange-correlation kernel remains poorly known Nazarov2009 ; Thiele2014 ; we expect it to be a small correction to the first term for many systems.

We note that the sum rules (25)-(26) discussed in Sec. II.4 become (see Appendix B.2)

[TABLE]

for all x,y. In practice, these sum rules provide a useful test of the accuracy or the consistency of the numerical calculations.

By using Eq.(51), the electron-ion coupling factor becomes

[TABLE]

We find important to remark that Eq.(51) is an exact rewriting of the Kubo formula (20). It does not correspond to the Kubo-Greenwood approximation that is widely used to calculate other transport properties such as the electrical and thermal conductivities. In the Kubo-Greenwood (KG) approximation, the Kohn-Sham wave functions are assumed to represent the single-particle excitations of the electronic system, i.e. the many-body electronic states are approximated by Slater determinant of KS orbitals. Here, this approximation amounts to set $\chi_{\mathrm{ee}}^{[\mathbf{R}]}$ in Eq.(LABEL:Eq:5) equal to $\chi_{\mathrm{KS}}^{[\mathbf{R}]}$, yielding

[TABLE]

which should be compared with Eq.(51). Here $\hat{f}_{\alpha}^{KG}(t)=-e^{i\hat{h}_{\mathrm{KS}}t/\hbar}\partial_{R_{\alpha}}\hat{V}_{ie}e^{-i\hat{h}_{\mathrm{KS}}t/\hbar}$ is the bare electron-ion force propagated by the Kohn-Sham Hamiltonian. Unlike the force $\tilde{f}_{\alpha}$ in Eq.(53), $\hat{f}_{\alpha}^{KG}$ fails to include the effects of the electronic screening on the electron-ion interactions. As a consequence, the KG approximation (57) can be shown to diverge logarithmically at large interparticle distances $|{\bf r}_{1}-{\bf r}_{2}|$. This is fully analogous to the well-known infrared divergence that occurs in the basic calculations of scattering cross sections in plasma physics when screening effects are neglected Spitzer1962 .

IV Relation to other models

The electron-ion coupling formula (24) includes self-consistently the quantum mechanical and statistical nature of electrons, the thermal effects and the correlations between all particles. In this section we show that our theory reduces to well-known models developed in plasma and condensed matter physics when one or several of these effects are treated approximately or neglected.

We begin in Sec. IV.1 and IV.2 with discussing the reduction to the Spitzer formula in the hot plasma limit Spitzer1962 , to the Fermi golden rule formula in the limit of weak electron-ion interactions Hazak2001 , and to the model developed by Daligault and Dimonte that includes important electron-ion correlation effects DaligaultDimonte2009 . We then explain that the problem of the calculation of the electron-ion coupling is particularly well adapted to average atom model calculations. In Sec. IV.3, we show that our theory also applies to hot solids, namely to solid metals with lattice and electronic temperatures much greater than the Debye temperature. When the ionic motions are described within the phonon approximation, our theory reproduces the standard electron-phonon coupling formula Allen1987 in the high temperature limit; the full formula includes effects beyond the harmonic approximation. Finally, in Sec. IV.4, we relate our theory to a simplified model of the electon-phonon coupling factor due to Wang et al. Wang1994 .

In the following, $v_{ie}(k)=\int{\!d{\bf r}v_{ie}(r)e^{-i{\bf k}\cdot{\bf r}}}$ and $v_{C}(k)=4\pi e^{2}/k^{2}$ are the spatial Fourier transforms of the electron-ion and Coulomb interaction potentials.

IV.1 Reduction to the Spitzer model and the Fermi golden rule formula

In the general formula (24), the effects of the interactions between the ionic and electronic subsystems are included non-perturbatively. Simpler models are obtained if one assumes that the interaction $v_{ie}$ is weak. Then, to lowest order in $v_{ie}$, the response function $\chi_{ee}^{[{\bf R}]}$ can be approximated by the density-density response function $\chi_{jel}$ of the homogeneous electron gas (a.k.a. jellium) at temperature $T_{e}$, i.e.

[TABLE]

With this approximation, the electron-ion coupling formula (24b) becomes

[TABLE]

In the thermodynamic limit, this expression is conveniently rewritten such as

[TABLE]

where $\chi_{jel}(k,\omega)=\int{\!d{\bf r}\chi_{jel}(r,\omega)e^{-i{\bf k}\cdot{\bf r}}}$. Similarly, the formulas (34) and (44) become

[TABLE]

and

[TABLE]

where $\tilde{\chi}$, $\chi_{0}$ and $G_{ee}$ are the proper response function, the non-interacting (a.k.a. Lindhard) response function and the local field correction of the jellium model GiulianiVignale2005 . Equations (59), (LABEL:gchi0approx1) and (LABEL:gchi0approx2) correspond to the so-called ‘Fermi golden rule formula’ for the electron-ion coupling derived in Ref. Hazak2001 by first calculating the energy exchanges between the electron and ion subsytems within the framework of linear response theory and then by taking the small ion to electron velocity ratio into account.

The relation of Eq.(59) to the celebrated Spitzer formula was discussed elsewhere (e.g., see DaligaultDimonte2009 ), and we only briefly recall the result for completeness. When electron-electron correlation effects are neglected, $\tilde{\chi}=\chi_{0}$ and $G_{ee}=0$, and

[TABLE]

The familar plasma physics results are recovered using the Coulomb interaction $v_{ie}(r)=-Ze^{2}/r$ in Eq.(LABEL:g_towards_plasmas), yielding

[TABLE]

in terms of the Coulomb logarithm

[TABLE]

with the screening length $\lambda_{sc}=1/\sqrt{4\pi e^{2}\chi_{0}(k,0)}$ and $f(k)$ is the Fermi-Dirac function (see Eq.(90)). In the classical limit $\hbar\to 0$, $f(k)=1$ and $\lambda_{sc}=\sqrt{k_{B}T_{e}/4\pi n_{e}e^{2}}$ is the Debye-Hückel screening length, and Eq.(63) becomes the celebrated Spitzer formula, which logarithmically diverges at large $k$. In the non-degenerate limit, Eq.(64) becomes $\ln\Lambda=\int_{0}^{\infty}{\frac{dk}{k}\frac{k^{4}}{(k^{2}+\lambda_{De}^{-2})^{2}}e^{-\lambda_{e}^{2}k^{2}/8}}$ is convergent, where $\lambda_{e}=\hbar/\sqrt{m_{e}k_{B}T_{e}}$ is the electronic thermal de Broglie wavelength.

IV.2 Going beyond the weak electron-ion interaction approximation using average atom models

In the previous section, we have discussed the general expression (24) in the limit of weak electron-ion interactions. Here, we discuss approaches to go beyond this approximation. To this end, we remark that the formula (10) can be written as

[TABLE]

where

[TABLE]

can be regarded as the averaged friction coefficient felt by any ion in the system along any direction of motion.

A natural approximation for $\Gamma$ consists in identifying it with the friction felt by a single impurity embedded in an electron gas (jellium) or by a slow projectile (i.e., whose velocity is much smaller than the electronic velocities),

[TABLE]

where ${\bf f}$ is the force between the impurity and the electrons. This problem has been extensively studied in the past (see, e.g. dAgliano1975 ) and, for completeness, we recall results useful to the present work. When the electron-impurity interaction is treated within the framework of linear response,

[TABLE]

When used in Eq.(65), we, not surprisingly, retrieve the Fermi golden rule formula (59) discussed above.

As suggested in Ref. Dufty1996 , the so-called ‘disconnected approximation’ can be used to extend Eq.(68) beyond the weak electron-ion interaction approximation. This approximation, which orginates from works in the classical kinetic theory of strongly coupled plasmas Gould1977 ; Boercker1982 , neglects the effect of the slow impurity on the electron dynamics, but accounts for the average distortion of the electronic density around the impurity. It amounts to replacing in Eq.(68) the term $v_{ie}(k)^{2}$ by $v_{ie}(k)^{2}[1-G_{ie}(k)]$, where $G_{ie}$ is the electron-ion local field correction $G_{ie}(k)$, yielding

[TABLE]

This corresponds to the model derived by Daligault and Dimonte in DaligaultDimonte2009 using a very different method. We refer to Ref. DaligaultDimonte2009 for a detailed discussion of Eq.(LABEL:g_DaligaultDimonte).

Let assume for the moment that electron interactions can be neglected ($V_{ee}\equiv 0$). Electrons then move independently in the potential of the ionic impurity and the Kubo relation (67) can then be expressed in terms of the basic scattering properties of the electron-ion potential $v_{ie}$ (see e.g. dAgliano1975 ), which gives

[TABLE]

where $\epsilon_{k}=\hbar^{2}k^{2}/2m_{e}$, $n_{\rm FD}(\epsilon)=1/(1+e^{-(\mu-\epsilon)/k_{B}T_{e}})$ is the Fermi-Dirac distribution, and $\sigma_{tr}$ is the cross section for binary collisions

[TABLE]

where $\delta_{l}(k)$ is the phase shift of the $l^{\rm th}$ partial wave at momentum $\hbar k$ calculated for the spherically symmetric $v_{ie}(r)$ (see appendix E). This formula is applicable to any temperature $T_{e}$. In particular, at $T_{e}=0$,

[TABLE]

and, at high $T_{e}$, the result (IV.2) agrees with that obatined from the classical Boltzmann-Lorentz kinetic theory FerzigerKaper1972

[TABLE]

with

[TABLE]

The result (IV.2) can be effectively used in conjunction with an average atom model to include the effects of electron and ion interactions that affect the electron-ion cross section in a plasma. Average atom models have been a quite popular approximate method to model both the equation-of-state and transport properties of dense ionized matter Faussurier2010 ; Faussurier2016 ; StarrettSaumon2013 ; Starrett2017 . They have proved to be accurate enough to be useful while being computationally much more expedient. An average atom model assumes that the physical system is spherically symmetric about a central nucleus and one calculates with finite-temperature density functional theory the electronic structure of the central ions and of the surrounding conduction electrons. As we have seen above the electron-ion coupling factor is related to the averaged friction coefficient felt by any ion in the system along any direction of motion; its calculation is thus particularly well adapted for this transport property. Many variations exist that differ in the description of the surrounding plasma, e.g. via boundary conditions or by coupling the model with the theory of fluids. Electrons are treated as independent particles subject to the Kohn-Sham potential $v_{KS}(r)$ given by the sum of the central nuclear potential, the spherically averaged distribution of surrounding ions, and the self-consistent Hartree and exchange correlation potentials. The expressions (IV.2) and (71) still apply but should be computed using the phase shifts $\delta_{l}(k)$ corresponding to the effective potential $v_{KS}(r)$. Results for the friction coefficient felt by a charge immersed in a homogeneous electron gas were presented in Ref. Nazarov2005 .

IV.3 Relation to the electron-phonon coupling factor

In hte case of solids, the energy exchanges between electrons and ions are generally described in terms of interactions between electrons and phonons. In the phonon approximation, the total Hamiltonian is given by $\hat{H}=\hat{H}_{\rm e}({\bf R}_{0})+\hat{H}_{\rm ph}+\hat{H}_{\rm e-ph}$, where $\hat{H}_{\rm e}({\bf R}_{0})$ is the electronic Hamiltonian in the potential of the Bravais lattice (${\bf R}^{0}$ denote the equilibrium lattice positions), $\hat{H}_{\rm ph}$ is the Hamiltonian of phonons, and $\hat{H}_{\rm e-ph}$ is the electron-phonon Hamiltonian (see Appendix D). Here we show that, for hot solids, the standard electron-phonon coupling factor can be readily derived from the electron-ion coupling formula (24). For simplicity, we use the assumption of thermal phonons characterized by a single temperature $T_{i}$, although recent investigations suggest that it could lead to marked disagreement with experimental observations Waldecker2016 ; Henighan2016 ; Maldonado2017 ; Lu1 .

Most of the works in condensed matter physics rely on a formula for the electron-ion coupling derived by Allen Allen1987 . Here, we prefer to work with a generalization of Allen’s formula that, unlike the latter, does not approximate from the outset the electrons with the Bloch states. This generalized formula naturally reduces to Allen’s formula in the appropriate limit as shown in appendix D. As shown in appendix D.1, the rate of change of the total electron energy due to the absorption and emission of phonon excitations can be written as

[TABLE]

with $n_{B}(x)=1/(e^{x}-1)$. Here $\chi_{{\bf G},{\bf G^{\prime}}}^{\rm ph}({\bf q},\omega)$ is the Fourier transformed density-density response function of the ion subsystem calculated in the phonon approximation (see appendix D.5), and $\chi_{{\bf G},{\bf G^{\prime}}}^{[{\bf R}_{0}]}({\bf q},\omega)$ is the Fourier transformed density-density response function $\chi_{ee}^{[{\bf R}_{0}]}$ of the electron system described by the Hamiltonian $\hat{H}_{e}$, i.e. of the system of interacting electrons in the background of the perfect lattice potential,

[TABLE]

As shown in appendix D.4, Eq.(IV.3) reduces to Allen’s work (see Eq.(6) in Allen1987 ) in the limit of Bloch electrons.

Equation (IV.3) is valid for any temperatures $T_{e,i}$. The rate of energy change (IV.3) can be simplified at high temperatures, i.e. for electronic and ion temperatures greater than the Debye temperature $\Theta_{D}$ ($k_{B}\Theta_{D}=\hbar\omega_{D}$ with the Debye frequency $\omega_{D}=\underset{{\bf q},\lambda}{\rm max}\{\omega_{{\bf q}\lambda}\}$ (typically $0.01-0.04$ eV AshcroftMerminbook )). In that limit, as shown in appendix D.2, Eq.(IV.3) simplifies to

[TABLE]

where the electron-phonon coupling factor $g_{e-ph}$ is given by

[TABLE]

Remarkably, this expression is readily obtained from our general formula (24) for the electron-ion coupling by substituting in Eq. (24b) the expression (LABEL:chi_eeph) for the electron-electron density response function (see appendix D.3 for the details). Since Eq.(24b) is an approximation of the density response function of electrons in a solid, the formula (24) extends the electron-phonon coupling formula (78). In particular, it goes beyond the harmonic phonon approximation as the ionic configurations ${\bf R}$ in Eq.(24) include those thermally sampled that are not described by the small harmonic lattice vibrations. Such anharmonic motions become increasingly important as one approaches melting conditions. We refer the reader to SimoniDaligaultPRL2019 for preliminary first-principle calculations of the electron-phonon coupling based on Eq.(24).

IV.4 Relation to the Lin et al. model

We finally relate our approach to the following model due to Wang et al. Wang1994 ; Lin_et_al_2008 ,

[TABLE]

obtained as a simplification valid at high temperatures of Allen’s electron-phonon coupling fomula Allen1987 . Here $g(\epsilon)$ is the electron density of states (DOS), which is computable with DFT, and $G_{0}^{\rm e-ph}=\pi\hbar k_{B}\lambda\langle\omega^{2}\rangle g(\epsilon_{F})$, where $\epsilon_{F}=k_{B}T_{F}$ is the Fermi energy, $\langle\omega^{2}\rangle$ is the second moment of the phonon spectrum, and $\lambda$ is the electron-phonon mass enhancement factor. In previous works, the prefactor $G_{0}^{\rm e-ph}$ was either set to match an experimental measurement at low electronic temperature Lin_et_al_2008 , or was calculated ab-initio Waldecker2016 ; JiZhang2016 . Although derived for crystalline solids, the model (79) was used in recent works on warm dense matter systems Leguay2013 ; Cho_et_al_2015 ; Jourdain_et_al_2018 ; Ogitsu2018 . Remarkably, an expression similar to Eq.(79) also results from Eq.(56) if one neglects the second term and if one assumes that the matrix elements between the Kohn-Sham states of the force operator $\delta\tilde{f}^{L,R}$ depend weakly on the energies and spatial directions, which yields

[TABLE]

as shown in appendix F. Here $G_{0}^{ei}=|\delta\tilde{f}|^{2}g^{[{\bf R}]}(\epsilon_{F})^{2}$, where $g^{[{\bf R}]}(\epsilon)$ is the density of states of the Kohn-Sham system in the frozen ionic configuration ${\bf R}$, and $\delta\tilde{f}$ is the characteristic matrix element. The formulas (79) and (80) highlight the interplay between the DOS and the distribution of electronic states, which, as shown by Lin et al. Lin_et_al_2008 , results in a strong dependence on the chemical composition and often on sharp variations with $T_{e}$. In SimoniDaligaultPRL2019 , we have compared our results to predictions based on (79) reported by others and on Eq.(80) with $G_{0}^{ei}$ set to reproduce the value of $G_{ei}$ at the lowest $T_{e}$ considered. We find that the simplified models (79) and (80) tend to overestimate the dependence on $T_{e}$ or predict variations at odds with the full calculation.

Acknowledgements.

This work was supported by the US Department of Energy through the Los Alamos National Laboratory through the LDRD Grant No. 20170490ER and the Center of Non-Linear Studies (CNLS). Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001).

Appendix A Miscellaneous properties of correlation and response functions

A detailed exposition of the definitions and properties recalled below can be found in standard textbooks, e.g., Kubo_book ; GiulianiVignale2005

A.1 Quantum correlation functions

In quantum statistical mechanics, there are several ways to measure the temporal correlations between two observables $\hat{A}$ and $\hat{B}$ at thermal equilibrium. These include: the canonical Kubo relation

[TABLE]

the symmetrized correlation function,

[TABLE]

and the unsymmetrized correlation function,

[TABLE]

where in this appendix $\langle\dots\rangle={\rm Tr}\left(e^{-\beta\hat{H}}...\right)/{\rm Tre^{-\beta\hat{H}}}$ indicates a canonical thermal average. In the classical limit, the three definitions are equivalent.

Lehmann representation. By expanding over the eigenspectrum of the Hamiltonian $\hat{H}$, $\hat{H}|n\rangle=E_{n}|n\rangle$, the Fourier transforms of the time correlation functions can be written as

[TABLE]

and

[TABLE]

where $P_{n}^{eq}=e^{-E_{n}/k_{B}T}/\sum_{m}{e^{-E_{m}/k_{B}T}}$ is the thermal population of state $n$. We recall the relation,

[TABLE]

A.2 Density correlation and response function

We recall the well-known (fluctuation-dissipation) relation

[TABLE]

between the symmetric density-density correlation function $S(\mathbf{r}_{1},\mathbf{r}_{2},t)=\frac{1}{2}\langle\,\delta\hat{n}_{\mathrm{e}}(\mathbf{r}_{1},t)\delta\hat{n}_{\mathrm{e}}(\mathbf{r}_{2},0)+\delta\hat{n}_{\mathrm{e}}(\mathbf{r}_{2},0)\delta\hat{n}_{\mathrm{e}}(\mathbf{r}_{1},t)\,\rangle$ and the density-density response function $\chi(\mathbf{r}_{1},\mathbf{r}_{2},t)=-\frac{i}{\hbar}\theta(t)\langle\,[\delta\hat{n}_{\mathrm{e}}(\mathbf{r}_{1},t),\delta\hat{n}_{\mathrm{e}}(\mathbf{r}_{2},0)]\,\rangle$.

A.3 Lehmann representations of the density response functions

The Fourier transform of the density-density response function is

[TABLE]

For a system of independent particles,

[TABLE]

where $n_{FD}(\epsilon)=1/\left[1+e^{-(\mu(T)-\epsilon)/k_{B}T}\right]$ is the Fermi-Dirac population, and $\mu(T)$ is the chemical potential.

A.4 Lindhard response function

The density-density response function $\chi^{0}(k,\omega)$ of a non-interacting electron gas at temeprature $T$ is given by GiulianiVignale2005

[TABLE]

where $\epsilon({\bf p})={\bf p}^{2}/2m$ is the energy of a particle of momentum ${\bf p}$. Equation (88) implies DaligaultDimonte2009

[TABLE]

with

[TABLE]

In the classical limit ($\hbar\to 0$),

[TABLE]

Appendix B Details on the derivation of the relations (LABEL:Eq:5b), (39) and (51).

Here we drop the superscript $[{\bf R}]$ indicating the dependence on the instantaneous ionic configuration, and the integral relations between the response functions are written in operator notations. The identity operator is denoted by $\rm I$, i.e. $I({\bf r},{\bf r}^{\prime})=\delta({\bf r}-{\bf r}^{\prime})$.

B.1 Relation (LABEL:Eq:5b) to the proper response

Using obvious notations, the Dyson equation (LABEL:Dysonproper) can be written as

[TABLE]

where we introduced the left and right dielectric functions

[TABLE]

with inverses

[TABLE]

The relation ${\varepsilon}_{\mathrm{L}}(\omega)*{\varepsilon}_{\mathrm{L}}^{-1}(\omega)={\rm I}$ implies

[TABLE]

Therefore, from Eq.(B.1),

[TABLE]

By using the last expression into Eq. (LABEL:Eq:5), we obtain the desired Eq.(LABEL:Eq:5b).

B.2 Relations (39) and (51) to the Lindhard and Kohn-Sham responses

Because of the close similarity between the Dyson equations (LABEL:chieechi0) and (LABEL:DysonchichiKS) satisfied by the free-electron and Kohn-Sham response, the derivations of Eqs.(39) and (51) are analogous. We here consider the case involving the free electron response function. the Dyson equation (39) can be written as

[TABLE]

with the frequency-dependent kernel

[TABLE]

and the left and right dielectric functions

[TABLE]

with inverses

[TABLE]

Following the steps used to derive Eq.(93), we now obtain

[TABLE]

where the additional term results from the dependence of the kernel on the frequency. Similarly, we find

[TABLE]

By introducing the last expressions into Eq. (LABEL:Eq:5), we readily obtain the desired relations (39) and (51).

B.3 Homogeneous limit of the response and dielectric functions

We give properties satisfied by the reponse function $\chi_{ee}$ and related quantities in the limit of a homogeneous electron gas. Similar relations are satisfied by $\tilde{\chi}$, $\chi_{0}$ and $\chi_{KS}$.

In the limit of a homogeneous electron gas,

[TABLE]

and the left and right dielectric functions are equal,

[TABLE]

The spatial Fourier transform, generally defined as

[TABLE]

satisfies

[TABLE]

The inverse dielectric function satisfies

[TABLE]

Finally, the integral equations in Sec. B.1 and B.2 become algebraic equation, e.g.,

[TABLE]

with $K({\bf k},\omega)=v_{C}({\bf k})\left[1-G_{ee}({\bf k},\omega)\right]$ and $K_{KS}({\bf k},\omega)=v_{C}({\bf k})+f_{xc}({\bf k},\omega)$

Appendix C Sum rules

In this appendix we provide the proof for a set of sum rules, the $[{\bf R}]$ superscript is dropped in order to simplify the notation. Below,

[TABLE]

i)

The force matrix elements and the total linear momentum matrix elements satisfy

[TABLE]

*Proof:

*By starting from the definition

[TABLE]

we may rewrite the term on the right hand side as follows

[TABLE]

where $\hat{{\bf P}}=\sum_{i=1}^{N_{e}}{\hat{\mathbf{p}}_{i}}$ is the linear momentum operator of the many body system. In deriving Eq.(105), we used the relation $\left[\sum_{i}\hat{{\bf p}}_{i},\hat{V}_{ee}\right]=0$ that results from the symmetry of the Coulomb interaction. By substituting the previous result into Eq. (104) we easily obtain the final result (103).

ii)

The density response function satisfies the relations:

[TABLE]

[TABLE]

*Proof:

*Here we limit ourself to prove the first expression, for the second one the procedure is completely analogous. By using the Lehmann representation for the electron-electron susceptibility (86)

[TABLE]

where we have used (103) for the force matrix elements, that proves Eq. (106).

iii)

The Kohn-Sham matrix elements satisfy

[TABLE]

*Proof:

*From the definition of the force matrix elements

[TABLE]

we rewrite the term on the right hand side as follows

[TABLE]

where $\hat{H}_{\rm KS}({\bf R})=-\frac{\hbar^{2}\nabla^{2}}{2m}+\hat{V}_{\rm KS}({\bf r})$ is the Hamiltonian of the Kohn-Sham system, such that $\hat{H}_{\rm KS}({\bf R})\ket{n}=\epsilon_{n}\ket{n}$. From this result (108) then follows immediately.

iv)

The Kohn-Sham density response function satisfies the relations:

[TABLE]

[TABLE]

*Proof:

*The procedure is analogous to the one used in $ii)$, with the only difference that now we need to use (108) instead of (103)

[TABLE]

that proves (109), while for (110) the proof is identical.

v)

The Kohn-Sham dielectric functions satisfy the relations:

[TABLE]

[TABLE]

*Proof:

*By using the definition of the left dielectric function

[TABLE]

into the left hand side of (111) we obtain, by using Eq. (106), (we omit the frequency dependence in the derivation)

[TABLE]

that proves (111), while (112) may be obtained in the same way by using the definition of the right dielectric function, $\tilde{\varepsilon}_{\rm R}$.

vi)

The friction coefficients satisfy the sum rule

[TABLE]

*Proof:

*From the definition of the many body friction frictions, we can write

[TABLE]

while the frequency derivative of the electron-electron susceptibility is

[TABLE]

by using (104) the combination of the previous two expressions leads to

[TABLE]

with $P_{nm}^{\rm eq}=P_{n}^{\rm eq}-P_{m}^{\rm eq}$ and $E_{nm}=E_{n}-E_{m}$, proving the sum rule (113).

vii)

An analogous result is valid also for the Kohn-Sham friction tensor

[TABLE]

*Proof:

*From the definition of the Kohn-Sham tensor (LABEL:Eq:8) and by using the sum rules (111) and (112) for the gradient of the external potential it is easy to write

[TABLE]

that finally proves (114). As a consequence of (114) and (111) we also have

[TABLE]

completing the set of sum rules we seek to prove.

Appendix D The electron-phonon coupling formula

D.1 Derivation of the electon-phonon coupling formula (IV.3)

We treat the general case of any monatomic Bravais lattice, whose ionic equilibrium positions are denoted by ${\bf R}^{0}=\{{\bf R}_{I}^{0}\}$.

The Hamiltonian of an electron gas interacting with a periodic lattice of ions oscillating around their equilibrium positions is $\hat{H}=\hat{H}_{\rm e}+\hat{H}_{\rm ph}+\hat{H}_{\rm e-ph}$, where BruusFlensberg2004

[TABLE]

is the Hamiltonian of the electron gas interacting with the ions in their equilibrium positions ${\bf R}^{0}$,

[TABLE]

is the Hamiltonian of the phonons, and

[TABLE]

is the interaction between electrons and phonons with the phonon coupling

[TABLE]

Here $\hat{b}_{{\bf q}\lambda},\hat{b}_{{\bf q}\lambda}^{\dagger}$ are the annihilation and creation operator of a phonon of frequency $\omega_{{\bf q},\lambda}$, $\bm{\epsilon}_{{\bf q}\lambda}$ are the polarization vectors, $\hat{\rho}_{\bf k}$ is the Fourier transform of the electron density, $\sum_{\lambda}$ is the sum over polarizabilities, $\sum_{{\bf q}}$ means ${\bf q}$ in first Brillouin zone, $\sum_{{\bf G}}$ means ${\bf G}$ in reciprocal lattice, k is a Brillouin zone’s vector, q is localized in the first Brillouin zone and G is a reciprocal space’s vector BruusFlensberg2004 . The interpretation of the previous expression is straightforward, the electron in fact can be scattered from any initial state $\ket{{\bf k}}$ to a final state $\ket{{\bf k}+{\bf G}+{\bf q}}$ either by absorbing a phonon in the state $\ket{{\bf q},\lambda}$ or by emitting a phonon in the state $\ket{-{\bf q},\lambda}$.

We shall apply the Fermi golden rule to calculate the rate of change of the total electron energy

[TABLE]

where $W_{\rm abs}({\bf q},\lambda)$ is the rate of absorption and $W_{\rm em}({\bf q},\lambda)$ is the rate of emission of a phonon of energy $\hbar\omega_{{\bf q}\lambda}$ by the electronic states $\ket{m}$ defined by $\hat{H}_{\rm e}({\bf R}^{0})\ket{m}=E_{m}\ket{m}$. We calculate these rates to lowest order of pertrubation theory by applying the Fermi golden rule.

[TABLE]

The propability per unit time of transition between state $\ket{m^{\prime}}\otimes|\dots n_{{\bf q}\lambda}\dots\rangle$ and state $\ket{m}\otimes|\dots(n_{{\bf q}\lambda}-1)\dots\rangle$ is

[TABLE]

By averaging over a thermal distribution of electronic states at temperature $T_{e}$ and of phonon states at temperature $T_{i}$, we obtain the rate of phonon absorption

[TABLE]

where $N_{{\bf q}\lambda}=1/(e^{\hbar\omega_{{\bf q},\lambda}/k_{B}T_{i}}-1)$ is the Bose population of the phonon mode (${\bf q},\lambda)$ at temperature $T_{i}$ and $P_{m}=e^{-E_{m}/k_{B}T_{e}}/{\cal{Z}}$ is the thermal population of the electronic state $|m\rangle$. In the second line, we have introduced the (non-symmetrical) electron density correlation function (see also appendix A)

[TABLE]

Similarly, the emission calculation gives

[TABLE]

where in the second line we used the detailed balance property GiulianiVignale2005 .

The rate (120) of energy exchange between electrons and phonons becomes

[TABLE]

In deriving Eq.(125), we used the simple relation

[TABLE]

and the fluctuation-dissipation relation

[TABLE]

between the correlation function and the density-density response function of electrons described by the Hamiltonian $\hat{H}_{\rm e}({\bf R}^{0})$ (see Eq.(LABEL:chi_eeph)). In going from equation (125) to the desired result (125), we used the expression for the density-density response of ions in the phonon approximation derived below in Sec. D.5,

[TABLE]

in terms of the phonon spectral function $\mathcal{A}_{\lambda}({\bf q},\omega)=\delta(\hbar\omega-\hbar\omega_{{\bf q}\lambda})-\delta(\hbar\omega+\hbar\omega_{{\bf q}\lambda})$.

Equation (125) was derived by treating the electrons as a many-body system, i.e. the states $|m\rangle$ are the many-body eigenstates of the Hamiltonian $\hat{H}_{e}({\bf R}^{0})$, Eq.(116). The considerations of Sec. III, where the many-body properties are expressed in terms of single-particle properties, can be straightforwardly adapted to effectively deal with electrons in a crystalline solid. For instance, instead of using the free electron response function as in Sec. III.2, the response function $\chi_{{\bf GG^{\prime}}}^{[{\bf R}_{0}]}({\bf q},\omega)$ can be expressed in terms of the response function of non-interacting electrons immersed in the perfect ion lattice ${\bf R}^{0}$ described by the single particle Hamiltonian,

[TABLE]

whose eigenstates are the so-called Bloch electron states.

D.2 High temperature limit $T_{i,e}\gg\Theta_{D}$.

We show that in the hot solid limit, $T_{i,e}\gg\Theta_{D}$, the relation (IV.3) simplifies to Eq.(77). The quantity ${\rm Im}\chi_{{\bf G},{\bf G^{\prime}}}^{\rm ph}({\bf q},\omega)$, which is simply related to the phonon spectrum (LABEL:Imchiphphspectrum), is non-zero only for frequency $|\omega|$ smaller than the Debye frequency $\omega_{D}$. For $k_{B}T_{i,e}\gg\hbar\omega_{D}$ and $0\leq\omega\leq\omega_{D}$, we have

[TABLE]

and, as a consequence of the small electron to ion mass ratio,

[TABLE]

Using these two approximations in Eq.(IV.3), we obtain

[TABLE]

The expression (LABEL:Imchiphphspectrum) implies the relation

[TABLE]

which, when introduced in Eq.(LABEL:dEedintermediatehighTlimit), implies the desired results (77) and (78).

D.3 Derivation of $g_{e-ph}$, Eq.(78), from the general formula (24)

By introducing the expression (LABEL:chi_eeph) in $g_{ei}$, Eq.(24b), the average over ions disappear (it is set to ${\bf R}^{0}$) and we obtain

[TABLE]

D.4 Reduction to Allen’s formula

In his paper, Allen describes the electrons in terms of Block states $\psi_{k}$, eigenstates of the Hamiltonian (127),

[TABLE]

where $k=(n,{\bf k},\sigma)$ is short for the Bloch electron quantum number (${\bf k}$ is in the first Brillouin zone, $n\in\mathbb{N}$ is the band index and $\sigma$ denotes the spin).

[TABLE]

The density response function of Bloch electrons satisfies

[TABLE]

with $\rho_{k,k^{\prime}}({\bf K})=\int_{V}{d{\bf r}u_{k^{\prime}}^{*}({\bf r})u_{k}({\bf r})e^{-i({\bf k}^{\prime}-{\bf k}+{\bf K})\cdot{\bf r}}}\delta_{\sigma,\sigma^{\prime}}$

[TABLE]

[TABLE]

where in the last equation we used the energy conservation described by the delta function in Eq.(D.4). where we have introduced the scattering amplitude probability

[TABLE]

Equation (D.4) corresponds to the starting point of Allen’s deriation, see equation (6) in Ref. Allen1987 .

D.5 The ionic density response function in the phonon approximation

Here we derive an expression for the density-density response function of ions

[TABLE]

in a solid at temperature $T_{i}$ in the phonon approximation, where $n_{i}({\bf k})=\sum_{I}e^{-i{\bf k}\cdot{\bf R}_{I}}$ is the ion density and $\langle\dots\rangle_{i}$ denotes the thermal average at temperature $T_{i}$, and $\delta\hat{n}_{i}=\hat{n}_{i}-\langle\hat{n}_{i}\rangle_{i}$. In the phonon approximation, ${\bf R}_{I}(t)={\bf R}_{I}^{0}+{\bf u}_{I}(t)$, where the harmonic displacement ${\bf u}_{I}$ of ion $I$ around its equilibrium position ${\bf R}_{I}^{0}$ is given by

[TABLE]

where ${\bf q}$ is in the first Brillouin zone. To lowest order in the displacements,

[TABLE]

This implies

[TABLE]

where $D_{\rm R}^{\lambda}({\bf q},t)=-\frac{i}{\hbar}\theta(t)\langle[\hat{A}_{{\bf q}\lambda}(t),\hat{A}_{{\bf q}\lambda}^{\dagger}(0)]\rangle$ with $\hat{A}_{{\bf q}\lambda}=\hat{b}_{{\bf q}\lambda}(t)+\hat{b}_{-{\bf q}\lambda}^{\dagger}$ is the retarded phonon Green’s function BruusFlensberg2004 . Using the definition of the spectral function $\mathcal{A}_{\lambda}({\bf q},\omega)=-{\rm Im}D_{\rm R}^{\lambda}({\bf q},\omega)/\pi$, we find the relation

[TABLE]

used in Sec. D.1.

Appendix E Quick remarks about Eq.(IV.2).

The passage from Eq.(67) to Eq.(IV.2) in the limit of non-interacting electrons is non-trivial but can be found in several papers; e.g., see Sec. Section III-c of Ref. dAgliano1975 . For completeness, we recall the main steps, which rely on standard results of scattering theory. The Kubo relation is developed as follows

[TABLE]

over the basis of the so-called scattering states defined as

[TABLE]

Here $|{\bf k}\sigma\rangle$ is a plane wave of momentum $\hbar{\bf k}$, energy $\epsilon_{k}=(\hbar{\bf k})^{2}/2m_{e}$ and spin $\sigma$; $\hat{t}^{\pm}=\hat{t}(\epsilon_{k}\pm 0^{+})$ and $\hat{G}^{\pm}=\hat{G}(\epsilon_{k}\pm 0^{+})$ with the t-matrix $\hat{t}(z)$ and resolvent operator $\hat{G}(z)=\left[z-\frac{\hat{p}^{2}}{2m_{e}}\right]$ satisfy the Lippman-Schwinger equation

[TABLE]

We then use the property

[TABLE]

which results from $\hat{F}_{x}=\frac{i}{\hbar}\left[\hat{p}_{x},\hat{v}_{ie}\right]$ and of properties of Eq.(134). We obtain

[TABLE]

Equation (IV.2) is then obtained from the well-known representation of the matrix elements $\langle{{\bf k}^{\prime}}\sigma^{\prime}|\hat{t}^{+}|{\bf k}\sigma\rangle$ in terms of the phase shifts $\delta_{l}(k)$ for the spherically symmetric potential $v_{ie}(r)$,

[TABLE]

with $\cos\Omega={{\bf k}^{\prime}}\cdot{\bf k}/k^{\prime}k$.

The $T_{e}=0$ limit (72) is obtained using $dn_{FD}(\epsilon)/d\epsilon\to\delta(\epsilon-\epsilon_{F})$ with $\epsilon_{F}=\hbar^{2}k_{F}^{2}/2m_{e}$ is the Fermi energy.

The non-degenerate limit (73) is obtained using $n_{FD}(\epsilon)[1-n_{FD}(\epsilon)]\sim e^{(\mu-\epsilon)/k_{B}T_{e}}$ and $e^{\mu/k_{B}T_{e}}=n_{e}\frac{(2\pi\hbar)^{3}}{2}\left(2\pi m_{e}k_{B}T_{e}\right)^{2/3}$.

Appendix F Derivation of Eq.(80).

First, we neglect the correction term $\delta\tilde{\gamma}_{\alpha\beta}^{[{\bf R}]}$ in Eq.(56) and expand the expression (53) for $\tilde{\gamma}_{\alpha\beta}^{[{\bf R}]}$ over the Kohn-Sham states as follows

[TABLE]

where the matrix elements f_{Ix}^{nm}=\left\langle n\big{|}\hat{f}_{Ix}^{(sc)}\big{|}m\right\rangle and $\hat{f}_{Ix}^{(sc)}$ is the effective force along the $x$-direction between ion $I$ and a Kohn-Sham electron screened by other electrons. Assuming that the matrix elements depend weakly on the energies, they can be factorized outside the sum in Eq.(135), and we obtain

[TABLE]

where $g^{[{\bf R}]}(\epsilon)=\sum_{n}\delta(\epsilon-\epsilon_{n})$ is the density of states of the Kohn-Sham system in the frozen ionic configuration ${\bf R}$.