Exact factorization-based density functional theory of electron-phonon systems

TL;DR

This paper extends density functional theory to include electron-phonon interactions using an exact factorization approach, revealing the importance of Berry curvature and providing a functional approximation that captures nonadiabatic effects on band structures.

Contribution

It introduces a novel exact factorization-based DFT framework for electron-phonon systems, incorporating Berry curvature effects and developing an orbital-dependent functional for nonadiabatic interactions.

Findings

Berry curvature influences phonon modes beyond potential energy surfaces.

The functional reproduces nonadiabatic band structure renormalization.

The approach generalizes DFT to more accurately model electron-phonon coupling.

Abstract

Density functional theory is generalized to incorporate electron-phonon coupling. A Kohn-Sham equation yielding the electronic density $n_U(\mathbf{r})$, a conditional probability density depending parametrically on the phonon normal mode amplitudes $U=\{U_{\mathbf{q}\lambda}\}$, is coupled to the nuclear Schr\"odinger equation of the exact factorization method. The phonon modes are defined from the harmonic expansion of the nuclear Schr\"odinger equation. A nonzero Berry curvature on nuclear configuration space affects the phonon modes, showing that the potential energy surface alone is generally not sufficient to define the phonons. An orbital-dependent functional approximation for the non\-adiabatic exchange-correlation energy reproduces the leading-order nonadiabatic electron-phonon-induced band structure renormalization in the Fr\"ohlich model.

Equations

$$\Psi (r,R)={\Phi }_R(r)\chi (R),$$

$$n_R(\mathbf{r})$$

$$=⟨{\Phi }_R|{\hat{\psi }}^{†}(\mathbf{r})\hat{\psi }(\mathbf{r})|{\Phi }_R⟩,$$

$$\Psi (r,R)=\sum _{J=0}^{\infty }{\Phi }_{JR}^{BO}(r){\chi }_J^{BO}(R),$$

$${\hat{H}}^{BO}$$

$$-\sum _{i,\mu }\frac{Z_{\mu }{e}^2}{4\pi {ϵ}_0|{\mathbf{r}}_i-{\mathbf{R}}_{\mu }|}.$$

$$n_R(\mathbf{r})=\frac{{\sum }_{J=0}^{\infty }|{\chi }_J^{BO}(R){|}^2⟨{\Phi }_{JR}^{BO}|{\hat{\psi }}^{†}(\mathbf{r})\hat{\psi }(\mathbf{r})|{\Phi }_{JR}^{BO}⟩}{{\sum }_{J=0}^{\infty }|{\chi }_J^{BO}(R){|}^2}.$$

$$n_R^{BO}(\mathbf{r})=⟨{\Phi }_{0R}^{BO}|{\hat{\psi }}^{†}(\mathbf{r})\hat{\psi }(\mathbf{r})|{\Phi }_{0R}^{BO}⟩.$$

$$\displaystyle E=T_{n,marg}+\int|\chi(R)|^{2}\big{[}\mathcal{E}^{BO}(R)+\mathcal{E}_{geo}(R)\big{]}dR{,}$$

$$T_{n,marg}$$

$${\mathcal{E}}^{BO}(R)$$

$${\mathcal{E}}_{geo}(R)$$

$$\displaystyle\big{[}\hat{H}^{BO}+\hat{U}_{en}\big{]}|\Phi_{R}\rangle$$

$$\displaystyle\left[\sum_{\mu}\frac{\big{(}\mathbf{P}_{\mu}+\mathbf{A}_{\mu}\big{)}^{2}}{2M_{\mu}}+\mathcal{E}(R)\right]\chi(R)$$

$$\left[\frac{{\mathbf{p}}^2}{2m}+v_{en}(\mathbf{r},R)+v_{hxc}^{EF}(\mathbf{r},R)\right]{\psi }_{iR}(\mathbf{r})={ϵ}_{iR}{\psi }_{iR}(\mathbf{r}),$$

$${\mathbf{R}}_{l\kappa }^{(0)}={\mathbf{R}}_l^{(0)}+{\mathbf{\tau }}_{\kappa },$$

$$\Psi (r,U)={\Phi }_U(r)\chi (U).$$

$${\mathbf{u}}_{l\kappa }$$

$${\hat{U}}_{\mathbf{q}\lambda }$$

$$\left[\frac{{\mathbf{p}}^2}{2m}+v_{en}(\mathbf{r},u)+v_{hxc}^{EF}(\mathbf{r},u)\right]{\psi }_{n\mathbf{k}u}(\mathbf{r})={ϵ}_{n\mathbf{k}u}{\psi }_{n\mathbf{k}u}(\mathbf{r}),$$

$$v_s^{(0)}(\mathbf{r})={v_{en}(\mathbf{r},u)|}_{{\mathbf{u}}_{l\kappa }=0}+{v_{hxc}^{EF}(\mathbf{r},u)|}_{{\mathbf{u}}_{l\kappa }=0}$$

$$\displaystyle\left[\sum_{l\kappa}\frac{\big{(}\mathbf{p}_{l\kappa}+\mathbf{A}_{l\kappa}(u)\big{)}^{2}}{2M_{\kappa}}+\mathcal{E}(u)\right]\chi(u)$$

$${\mathcal{E}}^{BO}(u)$$

$$\displaystyle\mathcal{E}_{geo}(u)=\sum_{l\kappa}\frac{\hbar^{2}\big{<}\nabla_{\mathbf{u}_{l\kappa}}\Phi_{u}\big{|}\big{(}1-|\Phi_{u}\rangle\langle\Phi_{u}|\big{)}\big{|}\nabla_{\mathbf{u}_{l\kappa}}\Phi_{u}\big{>}}{2M_{\kappa}}{,}$$

$$\mathcal{E}(u)$$

$$A_{l\kappa \alpha }(u)$$

$$\displaystyle C_{l\kappa\alpha l^{\prime}\kappa^{\prime}\alpha^{\prime}}=\frac{\partial^{2}\mathcal{E}}{\partial u_{l\kappa\alpha}\partial u_{l^{\prime}\kappa^{\prime}\alpha^{\prime}}}\bigg{|}_{u=0}$$

$${\hat{H}}_{nucl}^{(2)}$$

$$+\frac{1}{2}\sum _{l\kappa \alpha l^{\prime }{\kappa }^{\prime }{\alpha }^{\prime }}{C}_{l\kappa \alpha l^{\prime }{\kappa }^{\prime }{\alpha }^{\prime }}{\hat{u}}_{l\kappa \alpha }{\hat{u}}_{l^{\prime }{\kappa }^{\prime }{\alpha }^{\prime }}.$$

$${\hat{u}}_{l\kappa \alpha }$$

$${\hat{u}}_{\mathbf{q}\kappa \alpha }$$

$${\hat{p}}_{l\kappa \alpha }$$

$${\hat{p}}_{\mathbf{q}\kappa \alpha }$$

$${\hat{H}}_{nucl}^{(2)}$$

$$+\frac{M_0}{2}\sum _{\mathbf{q}}\sum _{\kappa \alpha {\kappa }^{\prime }{\alpha }^{\prime }}{\hat{u}}_{-\mathbf{q}\kappa \alpha }{D}_{\kappa \alpha {\kappa }^{\prime }{\alpha }^{\prime }}(\mathbf{q}){\hat{u}}_{\mathbf{q}{\kappa }^{\prime }{\alpha }^{\prime }},$$

$$A_{\kappa \alpha {\kappa }^{\prime }{\alpha }^{\prime }}^{(1)}(\mathbf{q})$$

Full text

I Introduction

The standard picture of interacting electrons and phonons in solids is a product of the Born-Oppenheimer (BO) approximation.Born and Oppenheimer (1927); Born and Huang (1954); Ziman (1960) State-of-the-art first-principles calculations of electron-phonon-coupling effects start from a density functional theory (DFT) calculation for the equilibrium crystal structure.Hohenberg and Kohn (1964); Kohn and Sham (1965) The resulting “clamped nuclei” electronic band structure depends on the BO approximation. The normal modes of vibration (phonons) and first- and second-order electron-phonon coupling matrix elements are calculated from the response of the BO potential energy surface and Kohn-Sham (KS) orbitals to small displacements in the atomic positions. The electronic band structure, phonon modes and electron-phonon coupling terms define a BO reference state that contains enough information to evaluate several observables, such as the electron-phonon coupling constant $\lambda$ and transition temperature $T_{c}$ in conventional superconductors such as MgB2,Kortus et al. (2001); Kong et al. (2001); Liu et al. (2001); Bohnen et al. (2001); Choi et al. (2002a, b); Golubov et al. (2002); Choi et al. (2003); Mitrović (2004); Choi et al. (2006); Cappelluti (2006); Eiguren and Ambrosch-Draxl (2008); Choi et al. (2009a); Calandra et al. (2010); De la Peña-Seaman et al. (2010); Margine and Giustino (2013); Aperis et al. (2015) electronic band structure renormalization, Allen and Heine (1976); Allen and Cardona (1981a); *allen1981erratum; Allen and Cardona (1983); King-Smith et al. (1989); Eiguren et al. (2003); Park et al. (2007); Giustino et al. (2008); Park et al. (2009); Eiguren et al. (2009); Giustino et al. (2010); Cannuccia and Marini (2011); Gonze et al. (2011); Antonius et al. (2015); Poncé et al. (2015); *ponce2017erratum; Antonius and Louie (2016); Nery and Allen (2016); Monserrat and Vanderbilt (2016); Allen and Nery (2017) and electronic mass enhancement and specific heat.Allen (1972); Golubov et al. (2002); Choi et al. (2003); Lin et al. (2008); Delaire et al. (2008); Choi et al. (2009a, b); Subedi and Singh (2009); Diakhate et al. (2011); Wei et al. (2013); Tütüncü and Srivastava (2013, 2015); Zocco et al. (2015); Li et al. (2015); Wiendlocha et al. (2016); Brown et al. (2016); Zheng and Margine (2017) Nevertheless, there is growing interest in capturing nonadiabatic electron-phonon effectsEngelsberg and Schrieffer (1963); Jarlborg (1992); Falter et al. (1995); Maksimov and Shulga (1996); Kulić (2000); Ferrari (2007); Basko et al. (2009); Dean et al. (2010); De Fillipis et al. (2010); Klimin et al. (2016); Ponosov and Streltsov (2016, 2017) by ab initio approaches that go beyond this BO reference state.van Leeuwen (2004); Lazzeri and Mauri (2006); Bock et al. (2006); Pisana et al. (2007); Piscanec et al. (2007); Calandra et al. (2007); Caudal et al. (2007); Saitta et al. (2008); Calandra et al. (2010); Gonze et al. (2011); Cannuccia and Marini (2012); Marini et al. (2015); Poncé et al. (2015); Antonius et al. (2015); M. d’Astuto et al. (2016); Gali et al. (2016); Nery and Allen (2016); Allen and Nery (2017); Giustino (2017); Long and Prezhdo (2017); Zhou et al. (2017); Caruso et al. (2017); Nery et al. (2018); Marini and Pavlyukh (2018); Novko (2018); Caruso et al. (2018)

Reliance on the BO approximation complicates subsequent many-body calculations. Since the electronic Hamiltonian is already included in the adiabatic potential energy surface from which the reference BO phonons are calculated, it is not straightforward to rigorously divide the original electron-nuclear Hamiltonian into electronic $\hat{H}_{elec}$, phononic $\hat{H}_{phon}$, and electron-phonon coupling $\hat{H}_{elec-phon}$ terms,van Leeuwen (2004) as typically done in setting up many-body perturbation theory. It is therefore difficult to avoid double counting electronic interactions, and the many-body formalisms that have been proposed Hedin and Lundqvist (1969); van Leeuwen (2004); Marini et al. (2015); Antonius et al. (2015) are still more complicated than the widely-used BO-based approach outlined above, although work in this direction is ongoing.Giustino (2017); Marini and Pavlyukh (2018); Karlsson and van Leeuwen (2018)

To avoid double-counting issues, it would be desirable to be able to calculate electronic and phononic observables within a formally-exact DFT-like framework. This is possible in multicomponent DFT,Kreibich and Gross (2001) where the functionals depend on both the electronic density $n(\mathbf{r})$ in the body-fixed frame and the $N_{n}$-body nuclear density $\Gamma(\mathbf{R}_{1},\mathbf{R}_{2},\ldots,\mathbf{R}_{N_{n}})$, provided these densities can be realized in a noninteracting system with appropriate scalar potentials (noninteracting $v$ representablility). A series of works Lüders et al. (2005); Marques et al. (2005); Floris et al. (2005, 2007) on superconducting DFT Oliveira et al. (1988) have also been formulated to include $\Gamma(\mathbf{R}_{1},\mathbf{R}_{2},\ldots,\mathbf{R}_{N_{n}})$. However, it has proven difficult to approximate the $\Gamma$-dependence of the exchange-correlation potentials in both multicomponent DFT and superconducting DFT. Additionally, at temperature $T=0$, there does not exist an auxiliary noninteracting system capable of reproducing the density and anomalous density in superconducting DFT.Schmidt et al. (2019)

In this paper, we focus on normal-state properties and show that a recent generalization of density functional theory Requist and Gross (2016); Li et al. (2018) based on the exact factorization (EF) of the electron-nuclear wavefunction into electronic and nuclear factors Hunter (1975); Gidopoulos and Gross (2014); Abedi et al. (2010) offers a promising alternative for calculating electronic and phononic observables. In contrast to multicomponent DFT’s, the basic variable is a conditional electronic density $n_{R}(\mathbf{r})$, a function which encodes the electronic density for each different set of nuclear coordinates $R=(\mathbf{R}_{1},\mathbf{R}_{2},\ldots)$. Working with $n_{R}(\mathbf{r})$ instead of the body-fixed-frame density $n(\mathbf{r})$ makes the exchange-correlation functionals in this theory closer to those of standard BO-based DFT.

Since EF-based DFT preserves the density-functional description of electronic structure that has made DFT so successful for solids, one can hope to obtain accurate approximations by building on the functionals of standard DFT. Following Ref. Li et al., 2018, we consider an approximation strategy that consists in adding a nonadiabatic correction term to a standard DFT functional, such as a local density approximation (LDA) Kohn and Sham (1965); von Barth and Hedin (1972) or a generalized gradient approximation (GGA).Perdew et al. (1996) Analytical calculations for the Fröhlich model prove that this approximation achieves the correct leading-order electron-phonon-coupling induced band structure renormalization, including the velocity renormalization near the Fermi energy.

Section II presents the general formalism of exact factorization-based density functional theory and its application to electron-phonon systems; Sec. III defines phonons; Sec. IV introduces our functional approximation; and Sec. V applies the theory to the Fröhlich model. Conclusions and an outlook on future developments are given in Sec. VI.

II Exact factorization DFT

II.1 Electron-nuclear DFT

The exact factorization method Hunter (1975); Gidopoulos and Gross (2014); Abedi et al. (2010) expresses the full electron-nuclear wavefunction as

[TABLE]

where $r=\{\mathbf{r}_{i}\}$ denotes the set of electronic coordinates and $R=\{\mathbf{R}_{\mu}\}$ denotes the set of nuclear coordinates. The key variable in exact factorization-based DFTRequist and Gross (2016); Li et al. (2018) is the conditional electronic density

[TABLE]

where $P(\mathbf{r},R)$ is the joint probability to find an electron at position $\mathbf{r}$ and the nuclei at positions $R$ and $P(R)$ is the marginal probability of finding the nuclei at $R$ regardless of where the electrons are.

The electronic density $n(\mathbf{r})$ in a standard DFT calculation, which we hereafter denote as $n_{R}^{BO}(\mathbf{r})$, is also a conditional density depending parametrically on $R$. To see what beyond-BO contributions $n_{R}^{BO}(\mathbf{r})$ is missing, consider the Born-Huang expansion Born and Huang (1954)

[TABLE]

where $\Phi_{JR}^{BO}(r)$ is the $J^{th}$ eigenstate of the BO Hamiltonian

[TABLE]

The exact conditional electronic density can be written in terms of the Born-Huang expansion as

[TABLE]

A standard DFT calculation gives only the single term

[TABLE]

In terms of the nuclear wavefunction $\chi(R)$ and conditional electronic wavefunction $\Phi_{R}(r)$, the total energy of the electron-nuclear system can be expressed as

[TABLE]

where

[TABLE]

and $\mathbf{A}_{\mu}\equiv\hbar\mathrm{Im}\langle\Phi_{R}|\nabla_{\mathbf{R}_{\mu}}\Phi_{R}\rangle$. Making the energy stationary with respect to variations of $\chi(R)$ and $\Phi_{R}(r)$ subject to the partial normalization condition $\int|\Phi_{R}(r)|^{2}dr=1$ for all $R$ leads to the following equations:Gidopoulos and Gross (2014)

[TABLE]

where $\hat{U}_{en}$ is a complicated operator that depends nonlinearly on $\chi$ and $\Phi_{R}$. The nuclear equation has the form of a conventional Schrödinger equation with an exact potential energy surface $\mathcal{E}=\mathcal{E}^{BO}+\mathcal{E}_{geo}$ and an exact induced vector potential $\mathbf{A}_{\mu}$.

Exact factorization-based DFTRequist and Gross (2016); Li et al. (2018) seeks to bypass the many-body electronic equation, Eq. (9), using in its place the conditional KS equation

[TABLE]

where $v_{en}(\mathbf{r},R)=-\sum_{\mu}Z_{\mu}e^{2}/4\pi\epsilon_{0}|\mathbf{r}-\mathbf{R}_{\mu}|$ and $v_{hxc}^{EF}(\mathbf{r})$ is a nonadiabatic Hartree-exchange-correlation potential.

II.2 Electron-phonon DFT

We now consider a stable crystal with a well-defined equilibrium lattice structure. Adopting notations similar to those in Refs. Kwok, 1967; Maradudin and Vosko, 1968; Giustino, 2017, we specify the equilibrium position of nucleus $\kappa$ in primitive cell $l=(l_{1},l_{2},l_{3})$ as

[TABLE]

where $\mathbf{R}_{l}^{(0)}=l_{1}\mathbf{a}_{1}+l_{2}\mathbf{a}_{2}+l_{3}\mathbf{a}_{3}$ is the position of primitive cell $l$ and $\bm{\tau}_{\kappa}$ is the position of nucleus $\kappa$ within the primitive cell; $\mathbf{a}_{i}$ are the primitive lattice vectors. The displacement of a nucleus from its equilibrium position is defined to be $\mathbf{u}_{l\kappa}=\mathbf{R}_{l\kappa}-\mathbf{R}_{l\kappa}^{(0)}$, and we denote the set of nuclear displacements as $u=\{\mathbf{u}_{l\kappa}\}$.

For electron-phonon systems, it is convenient to adopt Born-von Kármán boundary conditions and work with the phonon normal mode coordinates $U=\{U_{\mathbf{q}\lambda}\}$. Therefore, we introduce the factorization

[TABLE]

The phonon normal mode coordinates $U$ and their relationship to $u$ will be derived in the following section. There are important differences with respect to standard DFT, where the relationship between $u$ and $U$ is

[TABLE]

Here, $\mathbf{e}_{\kappa}(\mathbf{q}\lambda)$ is the polarization vector of the phonon normal mode, $M_{0}$ is an arbitrary reference mass, e.g. the proton mass, and $N$ is the number of primitive cells under Born-von Kármán boundary conditions. Throughout the paper, it is to be understood that the $\mathbf{q}=0$ acoustic modes are excluded from sums over the phonon quasimomentum.

The conditional KS equation in Eq. (11) becomes

[TABLE]

where $v_{en}(\mathbf{r},u)=-\sum_{l\kappa}Z_{\kappa}e^{2}/4\pi\epsilon_{0}|\mathbf{r}-\mathbf{R}_{l\kappa}^{(0)}-\mathbf{u}_{l\kappa}|$. If we set the displacements $\mathbf{u}_{l\kappa}$ to zero, then the KS potential

[TABLE]

has lattice translational symmetry. As in standard DFT, this allows us to label the KS orbitals with a band index $n$ and wavevector $\mathbf{k}$. In terms of the displacement coordinates $u$, Eq. (10) becomes

[TABLE]

where $\mathbf{p}_{l\kappa}=-i\hbar\nabla_{\mathbf{u}_{l\kappa}}$ and $\mathbf{A}_{l\kappa}=\hbar\mathrm{Im}\langle\Phi_{u}|\nabla_{\mathbf{u}_{l\kappa}}\Phi_{u}\rangle$. Equations (15) and (17) are the fundamental equations of EF-based DFT for electron-phonon systems. The exact potential energy surface $\mathcal{E}=\mathcal{E}^{BO}+\mathcal{E}_{geo}$ comprises a BO-like term

[TABLE]

and a geometric termRequist and Gross (2016); Requist et al. (2016)

[TABLE]

which is similar to a term that can be derived in the BO approximation.Berry (1989); Berry and Lim (1990); Berry and Robbins (1993) $\mathcal{E}_{geo}$ is a geometric quantity that can be written as the contraction of a Riemannian metric tensor and an inverse mass tensor.

As a consequence of imposing Born-von Kármán boundary conditions, the exact potential energy surface, induced vector potential, and total energy in Eq. (17) acquire a parametric dependence on $\mathbf{R}_{l}^{(0)}$, i.e., on the lattice vectors $\mathbf{a}_{i}$. The equilibrium values of $\mathbf{a}_{i}$ can be obtained by minimizing the total energy at the end of the calculation.

Although separating off the center-of-mass motion, as we did in writing Eq. (13), modifies the electronic and nuclear kinetic energy operators,Sutcliffe (2000) the exact factorization scheme can still be straightforwardly applied to the resulting Schrödinger equation (see the supplemental material of Ref. Requist and Gross, 2016). To keep our focus on the essential differences between the present theory and standard DFT calculations of electron-phonon systems, we neglect these modifications and, moreover, we restrict our attention to nonpolar solids.

The induced vector potential $\mathbf{A}_{l\kappa}$ in Eq. (17) is said to be trivial if there exists a gauge choice such that $\mathbf{A}_{l\kappa}=0$. This is not always the case.Requist et al. (2016, 2017) In the following section, we show that the induced magnetic field (the curl of $\mathbf{A}_{l\kappa}$) affects the phonons.

III Exact phonons

Phonons are usually calculated in the BO approximation. The nuclear Schrödinger equation (17) affords us a way of defining “exact phonons.”

We start by expanding $\mathcal{E}(u)$ and $\mathbf{A}_{l\kappa}(u)$ as

[TABLE]

where $\alpha=(x,y,z)$, $A_{l\kappa\alpha,l^{\prime}\kappa^{\prime}\alpha^{\prime}}^{(1)}=\partial A_{l\kappa\alpha}/\partial u_{l^{\prime}\kappa^{\prime}\alpha^{\prime}}|_{u=0}$, and

[TABLE]

is the force constant matrix. We have assumed that the equilibrium coordinates $\mathbf{R}_{l\kappa}^{(0)}$ coincide with the coordinates $\mathbf{R}_{l\kappa}^{(min)}$ that minimize $\mathcal{E}(R)$. This may not always be the case, particularly if the phononic wavefunction is delocalized on a strongly anharmonic potential energy surface. The constant term $A_{l\kappa\alpha}^{(0)}$ in the expansion of $A_{l\kappa\alpha}$ can be removed by a gauge transformation $\Phi_{u}=\tilde{\Phi}_{u}e^{-(i/\hbar)A^{(0)}_{l\kappa\alpha}u_{l\kappa\alpha}}$ and is therefore inconsequential. Thus, to second order, the Hamiltonian in Eq. (17) is

[TABLE]

Within the BO approximation, this form of Hamiltonian has been considered previously.Holz (1972); Zhang et al. (2010) To find the eigenstates of $\hat{H}_{nucl}^{(2)}$, we first define the Fourier transformations

[TABLE]

and

[TABLE]

In the $\mathbf{q}$ representation, the Hamiltonian becomes

[TABLE]

where

[TABLE]

and the operators $\hat{u}_{\mathbf{q}\kappa\alpha}$ and $\hat{p}_{\mathbf{q}\kappa\alpha}$ satisfy the commutation relations

[TABLE]

To diagonalize the Hamiltonian in Eq. (25), we first apply the canonical transformation

[TABLE]

The inverse transformation is

[TABLE]

$\hat{U}_{\mathbf{q}\lambda}$ and $\hat{P}_{\mathbf{q}\lambda}$ should satisfy the commutation relations

[TABLE]

The transformation in Eqs. (28) and (29) will preserve the commutation relations if the polarization vectors $e_{\kappa\alpha}(\mathbf{q}\lambda)$ and $d_{\kappa\alpha}(\mathbf{q}\lambda)$ satisfy the orthonormality conditions

[TABLE]

and

[TABLE]

The transformation can be summarized with the help of a matrix $V$ as

[TABLE]

where $[\hat{u}_{\mathbf{q}}]$, $[\hat{p}_{\mathbf{q}}]$, $[\hat{U}_{\mathbf{q}}]$, $[\hat{P}_{\mathbf{q}}]$, $[e(\mathbf{q})]$ and $[d(\mathbf{q})]$ are columns and matrices indexed by $\kappa\alpha$ and $\lambda$. Then, the orthonormality constraints can be succinctly expressed as

[TABLE]

Applying the above transformation to Eq. (25) and requiring the coefficients of the $\hat{U}_{-\mathbf{q}\lambda}\hat{P}_{\mathbf{q}\lambda}$ and $\hat{P}_{-\mathbf{q}\lambda}\hat{U}_{\mathbf{q}\lambda}$ terms to vanish leads to the following eigenvalue equations:

[TABLE]

where $A^{(1)}$ and $D$ are the matrices in Eq. (26), $I$ is the identity matrix, and $1/\mathcal{M}$ and $\mathcal{M}\Omega^{2}$ are eigenvalues; $\mathbf{q}\lambda$ has been suppressed. The two types of eigenvectors in Eq. (59) are

[TABLE]

The identities $D_{\kappa\alpha\kappa^{\prime}\alpha^{\prime}}(-\mathbf{q})=D_{\kappa\alpha\kappa^{\prime}\alpha^{\prime}}^{*}(\mathbf{q})$ and $A_{\kappa\alpha\kappa^{\prime}\alpha^{\prime}}^{(1)}(-\mathbf{q})=A_{\kappa\alpha\kappa^{\prime}\alpha^{\prime}}^{(1)*}(\mathbf{q})$, which follow from the definitions in Eq. (26), together with $e_{\kappa\alpha}(-\mathbf{q})=e_{\kappa\alpha}^{*}(\mathbf{q})$ and $d_{\kappa\alpha}(-\mathbf{q})=d_{\kappa\alpha}^{*}(\mathbf{q})$, imply that if $v_{P}(\mathbf{q}\lambda)$ is an eigenvector with eigenvalue $1/\mathcal{M}_{\mathbf{q}\lambda}$, then $v_{P}(-\mathbf{q}\lambda)$ is also an eigenvector with the same eigenvalue, i.e. $\mathcal{M}_{\mathbf{q}\lambda}=\mathcal{M}_{-\mathbf{q}\lambda}$. Similar considerations for $v_{U}(\mathbf{q}\lambda)$ imply $\Omega_{\mathbf{q}\lambda}=\Omega_{-\mathbf{q}\lambda}$. Equation (59) replaces the standard eigenvalue equation $De=\omega^{2}e$ defining the phonons in terms of the dynamical matrix $D$. In Sec. IV.2, we discuss the differences between phonon calculations in our theory and standard calculations in density functional perturbation theory.Baroni et al. (2001); Savrasov and Savrasov (1996); Gonze and Lee (1997)

After these preliminaries, the transformed Hamiltonian can be expressed as

[TABLE]

where $\mathcal{M}_{\mathbf{q}\lambda}$ enters as a $\mathbf{q}\lambda$-dependent effective mass. Finally, in terms of the creation and annihilation operators

[TABLE]

the harmonic phonon Hamiltonian becomes

[TABLE]

Since this is bilinear in $a_{\mathbf{q}\lambda}^{{\dagger}}$ and $a_{\mathbf{q}\lambda}$, there are no phonon interactions at this order. It was for the purpose of obtaining this result that the expansion of $A_{l\kappa\alpha}$ was terminated at the first order; the second-order terms would have generated terms in Eq. (75) that are cubic and quartic in $a_{\mathbf{q}\lambda}^{{\dagger}}$ and $a_{\mathbf{q}\lambda}$.

The occurrence of a nonvanishing Berry curvature $B_{l\kappa\alpha,l^{\prime}\kappa^{\prime}\alpha^{\prime}}=A^{(1)}_{l^{\prime}\kappa^{\prime}\alpha^{\prime}l\kappa\alpha}-A^{(1)}_{l\kappa\alpha,l^{\prime}\kappa^{\prime}\alpha^{\prime}}$ on nuclear configuration space implies time-reversal symmetry breaking; this occurs naturally if the electronic state breaks time-reversal symmetry, e.g. in magnetic or (anomalous) quantum Hall systems. Our analysis is similar to that of Ref. Holz, 1972, where external rather than induced magnetic fields were considered.

In the special case $A^{(1)}=0$, $\mathcal{M}_{\mathbf{q}\lambda}\rightarrow M_{0}$, $\Omega_{\mathbf{q}\lambda}^{2}\rightarrow\omega_{\mathbf{q}\lambda}^{2}$, and Eq. (59) reduces to a single eigenvalue equation

[TABLE]

and the transformation in Eqs. (28) and (29) reduces to Eq. (14). Hence, in this case Eq. (73) recovers the standard Hamiltonian

[TABLE]

It will be helpful to write the explicit harmonic ground state wavefunction in the $U$ representation. Using the orthogonality of the phonon modes, it is

[TABLE]

where $L_{\mathbf{q}\lambda}$ is the amplitude of zero-point motion

[TABLE]

Phonon-phonon interactions arise from the anharmonicity of the potential energy surface and higher-order terms in the expansion of $A_{l\kappa\alpha}$, e.g. at third-order

[TABLE]

The full Schrödinger equation for the phonons then has the form

[TABLE]

and, in practice, the phonon-phonon interactions must be truncated at some order.

IV Nonadiabatic Hartree-exchange-correlation functional

From now on, we assume that the induced vector potential $\mathbf{A}_{l\kappa}$ in Eq. (17) is trivial; this assumption can be relaxed. Using the transformation [Eq. (14)] from nuclear displacements $u=\{\mathbf{u}_{l\kappa}\}$ to phonon normal mode coordinates $U=\{U_{\mathbf{q}\lambda}\}$, Eqs. (15) and (17) become

[TABLE]

where $\hat{P}_{\mathbf{q}\lambda}=-i\hbar\partial/\partial U_{-\mathbf{q}\lambda}$, $\mathcal{E}(U)=\mathcal{E}^{BO}(U)+\mathcal{E}_{geo}(U)$,

[TABLE]

and we recall that $\mathbf{q}=0$ acoustic modes are omitted from all sums.

The total energy of the electron-phonon system is

[TABLE]

One of the advantages of unifying electrons and phonons in a DFT framework is that a single density functional approximation for $\mathcal{E}(U)$ determines, on equal footing, all of the potentials in Eqs. (81) and (82).

As in standard DFT, the conditional electronic density is obtained from the occupied orbitals according to

[TABLE]

where $f_{n\mathbf{k}U}$ is a $U$-dependent occupation number. From now on, we suppress the subscript $U$ on the occupation numbers and orbitals.

Reference Requist and Gross, 2016 introduced an exact factorization-based DFT in which the energy is expressed as a variational functional of $(n_{R},\mathbf{j}_{pR},A_{\mu},\mathcal{T}_{\mu\nu},\chi)$, where $\mathbf{j}_{pR}$ is the conditional electronic paramagnetic current density, $A_{\mu}$ is the induced vector potential and $\mathcal{T}_{\mu\nu}$ is the quantum geometric tensor.Berry (1989); Provost and Vallee (1980) Reference Li et al., 2018 showed that the energy can also be expressed as a functional of $(n_{R},|\chi|^{2})$. We take a similar approach here and interpret the energy of an electron-phonon system as a functional of $(n_{U},|\chi|^{2})$. Equations (81) and (82) are then coupled through the functional dependence of the potentials: $\mathcal{E}$ depends on the density $n_{U}$, while $v_{hxc}^{EF}$ depends on $|\chi|^{2}$ (and $n_{U}$).

IV.1 Approximation strategy

Our strategy for approximating $\mathcal{E}[n_{U}]$ is the following. First, noting that $\mathcal{E}^{BO}[n_{U}]$ can be written as in DFT as

[TABLE]

we approximate $E_{hxc}[n_{U}]$ by a standard semilocal BO-based DFT functional $E_{hxc}^{BO}[n_{U}]$ such as a GGA. Second, in the nonadiabatic term $\mathcal{E}_{geo}$ in Eq. (83), we approximate the correlated electronic wavefunction $|\Phi_{U}\rangle$ by the Slater determinant of occupied KS orbitals. This defines an orbital-dependent functional $\mathcal{E}_{geo}[\psi_{n\mathbf{k}}]$, which is only implicitly a functional of $n_{U}$. The essential feature of this approximation is that nuclear mass-dependent, nonadiabatic effects are described by a simple additive correction to an existing DFT functional.

With this approximation for $\mathcal{E}_{geo}$, we have

[TABLE]

Via a chain rule for orbital-dependent functionals,Kümmel and Kronik (2008) e.g.

[TABLE]

the above approximations yield the scalar nonadiabatic Hartree-exchange-correlation potential in Eq. (81), i.e.

[TABLE]

where $v_{hxc}^{BO}(\mathbf{r},U)=\delta E_{hxc}^{BO}/\delta n_{U}(\mathbf{r})$ is the standard DFT potential and $v_{geo}(\mathbf{r},U)$ is a nonadiabatic correction.

Our analytical calculations for the Fröhlich model in Sec. V suggest that using a nonlocal (orbital-dependent) exchange-correlation potential is a more natural way to incorporate electron-phonon coupling. Since $\mathcal{E}^{BO}[n_{U}]$ can be converted into an orbital-dependent functional by substituting $n_{U}=\sum_{n\mathbf{k}}f_{n\mathbf{k}}|\psi_{n\mathbf{k}}|^{2}$, our approximation $\mathcal{E}_{geo}[\psi_{n\mathbf{k}}]$ implies an approximation for $\mathcal{E}[\psi_{n\mathbf{k}}]=\mathcal{E}^{BO}[\psi_{n\mathbf{k}}]+\mathcal{E}_{geo}[\psi_{n\mathbf{k}}]$ and, in turn, the total energy in Eq. (84). The stationary conditions with respect to $\psi_{n\mathbf{k}}^{*}(\mathbf{r})$ lead to a nonadiabatic and nonlocal generalized KS potential of the form:

[TABLE]

where $\hat{v}_{en}(U)$ and $\hat{v}_{hxc}^{\rm BO}(U)$ are the usual local potentials $v_{en}(\mathbf{r},U)$ and $v_{hxc}^{\rm BO}(\mathbf{r},U)$, and $\hat{v}_{geo}(U)$ is defined by its matrix elements

[TABLE]

The first term, hereafter denoted as $\langle\psi_{m\mathbf{k}+\mathbf{q}}|\hat{v}_{geo}^{(1)}|\psi_{n\mathbf{k}}\rangle$, is first order in $g$. The second term is second order. For the nuclear wavefunction in Eq. (78), we have

[TABLE]

so that to leading order we can write

[TABLE]

The factor $\langle\psi_{m\mathbf{k}+\mathbf{q}}|\partial\psi_{n\mathbf{k}}/\partial U_{\mathbf{q}\lambda}\rangle$ would in practice be determined self-consistently during the solution of Eqs. (81) and (82). Similarly, from the stationary condition with respect to variations of $f_{n\mathbf{k}}$, we obtain a second-order diagonal contribution

[TABLE]

Equations (81) and (82) together with $E_{hxc}^{BO}[n_{U}]$ and Eqs. (87), (93) and (94) completely determine the exact factorization DFT equations to second order in $g$. In Sec. V, we apply these equations to the Fröhlich model and demonstrate that they exactly recover the leading-order nonadiabatic electron-phonon coupling effects.

The operator $\hat{v}_{s}(U)$ in Eq. (90) is not a scalar multiplicative potential and therefore takes us outside a strict KS framework. As a result, the single-particle orbitals $\psi_{n\mathbf{k}}(\mathbf{r})$ will not generally equal the KS orbitals.Kümmel and Kronik (2008)

IV.2 Evaluating the force constant matrix

Given any approximation for $\mathcal{E}(u)$, we can evaluate the force constant matrix $C_{l\kappa\alpha l^{\prime}\kappa^{\prime}\alpha^{\prime}}=\partial^{2}\mathcal{E}/\partial u_{l\kappa\alpha}\partial u_{l^{\prime}\kappa^{\prime}\alpha^{\prime}}$. By virtue of the transformation in Eq. (14), our approximation $\mathcal{E}(U)=T_{s}[n_{U}]+\int V_{en}(\mathbf{r},U)n_{U}(\mathbf{r})d\mathbf{r}+V_{nn}(U)+E_{hxc}^{BO}[n_{U}]+\mathcal{E}_{geo}[\psi_{n\mathbf{k}U}]$ implies an approximation for $\mathcal{E}(u)$. Since the $u$ dependence enters both explicitly, through $V_{en}$ and $V_{nn}$, and implicitly, through the functional dependence on $n_{u}$ and $\psi_{n\mathbf{k}u}$, we obtain

[TABLE]

The last term generally leads to many terms involving the chain rule, e.g.

[TABLE]

The first-order density response $\partial n_{u}(\mathbf{r})/\partial u_{l\kappa\alpha}$ plays an essential role in Eq. (95), just as it does in density functional perturbation theory (DFPT)Baroni et al. (2001); Savrasov and Savrasov (1996); Gonze and Lee (1997) [c.f. Eq. (10) in Ref. Baroni et al., 2001]. There are a few important distinctions between phonon calculations in EF-based DFPT and standard BO-based DFPT. First, we cannot use the Hellmann-Feynman theorem, since $\mathbf{u}_{l\kappa}$ are not merely parameters in the exact electronic Schrödinger equation [the operator $\hat{U}_{en}$ in Eq. (9) contains the gradient $\nabla_{\mathbf{R}_{l\kappa}}=\nabla_{\mathbf{u}_{l\kappa}}$]. As a result, the Hessian of $\mathcal{E}$ also depends on the second-order density response $\partial^{2}n_{u}(\mathbf{r})/\partial u_{l\kappa\alpha}\partial u_{l^{\prime}\kappa^{\prime}\alpha^{\prime}}$. Second, our theory includes an induced vector potential $\mathbf{A}_{l\kappa}$ in the exact nuclear Schrödinger equation, which, if nontrivial, affects the phonon modes, showing that the force constant matrix alone is generally not sufficient to define the exact phonons. Through the self-consistent solution of Eqs. (81) and (82), we achieve a nonadiabatic extension of standard DFPT. Only marginally more computational time and resources are needed for an EF-based DFPT calculation than for a standard DFPT calculation.

V Fröhlich model

Here we consider an application of the above theory to the Fröhlich model with Hamiltonian

[TABLE]

where the electron-phonon interaction is

[TABLE]

We further simplify this to a single free-electron-like band and a single phonon mode in one dimension, i.e.

[TABLE]

The electronic states are denoted as $|\psi_{k}\rangle=c_{k}^{{\dagger}}|\mathrm{elec\;vac}\rangle$ and a general multiphonon state as

[TABLE]

where $n_{q}$ stands for the number of phonons in mode $q$.

V.1 First-order conditional density

We will apply perturbation theory for weak electron-phonon coupling $g_{k+q,k}$. The ground state for $g_{k+q,k}=0$ will be denoted as

[TABLE]

where $|\rm{FS}\rangle$ is the electronic Fermi sea and $|\chi_{0}\rangle$ is the vibrational ground state in Eq. (78). The excited states to which $|\Psi_{0}\rangle$ couples under $\hat{H}_{1}$ will be denoted as

[TABLE]

The first-order contribution to the wave function is

[TABLE]

In the $U$-coordinate representation, the wavefunction is

[TABLE]

Now we use the exact factorization method to derive the conditional electronic density to first order in $g$. The nuclear wave function is

[TABLE]

where the inner product is on the electronic Hilbert space only. From the conditional electronic wave function

[TABLE]

we obtain the zeroth-order and first-order contribution to the conditional electronic density

[TABLE]

$n_{U}^{(1)}(r)$ encodes how the conditional density is perturbed by the electron-phonon interaction.

V.2 Geometric correction

As a preliminary step, we expand $|\chi(U)|^{2}$ as

[TABLE]

The series has only even contributions. Choosing a gauge in which $\chi$ is real, we can write $\chi(U)=\chi^{(0)}(U)+\chi^{(2)}(U)+\mathcal{O}(g^{4})$ with

[TABLE]

where

[TABLE]

The conditional electronic wavefunction to third order is

[TABLE]

From Eq. (111), we find that the only contribution to $\mathcal{E}_{geo}$ through second order comes from the term

[TABLE]

V.3 Generalized Kohn-Sham system

The generalized KS potential

[TABLE]

depends parametrically on $U_{q}$, e.g. in the potential

[TABLE]

the atomic coordinates $\{R_{\mu}\}$ are implicit functions of the phonon amplitudes $U=\{U_{q}\}$.

We now apply perturbation theory to the KS system to see how it reproduces the results of the previous sections, particularly Eq. (107). The KS potential is expanded as

[TABLE]

where the superscript denotes the order in powers of $U_{q}$ (or, equivalently, in powers of $g$), $\hat{v}_{en}^{(0)}$ is the potential at the equilibrium atomic coordinates $R_{\mu}^{(0)}$, and $\hat{v}_{hxc}^{(n)}=0$ for all $n$. The unperturbed potential leads to the zeroth-order KS orbitals $\psi_{k}(r)$ through solution of the unperturbed KS equation

[TABLE]

The perturbations $\hat{v}_{en}^{(1)}$, $\hat{v}_{geo}^{(1)}$ and $\hat{v}_{geo}^{(2)}$ are defined by their matrix elements [cf. Eqs. (93) and (94)]

[TABLE]

To proceed, we need to recognize that the right-hand side of Eq. (118) is itself dependent on $\hat{v}_{geo}^{(1)}$. From perturbation theory applied to the KS equation, we obtain

[TABLE]

Equations (118) and (120) lead to the following differential equation for $G\equiv\langle\psi_{k+q}|\hat{v}_{geo}^{(1)}|\psi_{k}\rangle$:

[TABLE]

We choose the particular solution

[TABLE]

Using the first-order KS orbitals

[TABLE]

we immediately recover the result in Eq. (107), namely

[TABLE]

Thus, the generalized KS system with our nonadiabatic functional approximation reproduces the exact linear response density. Remarkably, the nonadiabatic potential $\hat{v}_{geo}$ has the effect of inserting $\hbar\omega_{q}$ into the denominator, thus recovering the expected nonadiabatic correction.

To determine $\langle\psi_{k}|\hat{v}_{geo}^{(2)}|\psi_{k}\rangle$ in Eq. (119), we first use perturbation theory to show that

[TABLE]

which implies

[TABLE]

Finally, is easy to show that

[TABLE]

reproduces Eq. (112).

V.4 Electronic band structure renormalization

The first-order correction to the KS eigenvalues

[TABLE]

vanishes since $\hat{v}_{en}^{(1)}$ and $\hat{v}_{geo}^{(1)}$ are off-diagonal.

The electronic velocity renormalization (the “wiggle”) at the Fermi energy appears in the second-order correction

[TABLE]

To obtain the observable perturbation we average over $U$ using $|\chi(U)|^{2}$ as a weighting function. The final result

[TABLE]

agrees with real part of the Fan-Migdal self-energyFan (1951) $\mathrm{Re}\Sigma_{nnk}^{FM}(\epsilon_{k})$ at $T=0$ and therefore encodes the correct electronic velocity renormalization.

A second-order electron-phonon interaction, called the Debye-Waller term, provides another contribution to electronic band structure renormalization.Antončik (1955); Allen and Heine (1976); Allen and Cardona (1981a) We do not consider it here, as it is not present in the Fröhlich model, although in real materials its contribution can be of the same order as the Fan-Migdal contribution.

VI Conclusions

Exact factorization-based DFT has been applied to interacting electrons and phonons in solids. The equations to be solved are (i) a generalized KS equation with a nonadiabatic Hartree-exchange-correlation potential that depends on the nuclear wavefunction $\chi$ and (ii) a nuclear Schrödinger equation with a beyond-BO potential energy surface and induced vector potential. Exact phonons are defined from the harmonic expansion of the nuclear Schrödinger equation without additional approximations.

We have proposed an approximation strategy in which nonadiabatic contributions to the KS potential and nuclear PES appear as simple additive corrections. For the Fröhlich model, the self-consistent solution of (i) and (ii) within our approximation recovers the exact electron-phonon-induced first-order density response and second-order electronic band structure renormalization. This suggests that we can obtain good results for electron-phonon effects in real materials by adding these nonadiabatic corrections to existing DFT functionals such as the LDA and GGA.

Subjects for future work are the formulation of a finite temperature theory and an investigation of the simultaneous effects of electron-electron and electron-phonon interactions, which, in principle, can be described exactly through the nonadiabatic Hartree-exchange-correlation potential $v_{hxc}^{EF}$. Lastly, the formalism introduced here provides an efficient methodology for predicting the effect of lattice degrees of freedom on geometric and topological properties of electronic Bloch states, such as the macroscopic polarization and topological invariants.

Acknowledgements.

R. R. thanks P. B. Allen for comments on the manuscript and R. van Leeuwen for discussions. C. R. P. thanks Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET) for partial financial support, grant PIP 2014-2016, and ANCyT under grant PICT 2016-1087.