{'fraction_non_alphanumeric': 0.043589807665075676, 'fraction_numerical': 0.020578754029537447, 'mean_word_length': 4.668153918791313, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Electron-phonon interactions from first principles

This article reviews the theory of electron-phonon interactions in solids from the point of view of ab-initio calculations. While the electron-phonon interaction has been studied for almost a century, predictive non-empirical calculations have become feasible only during the past two decades. Today it is possible to calculate from first principles many materials properties related to the electron-phonon interaction, including the critical temperature of conventional superconductors, the carrier mobility in semiconductors, the temperature dependence of optical spectra in direct and indirect-gap semiconductors, the relaxation rates of photoexcited carriers, the electron mass renormalization in angle-resolved photoelectron spectra, and the non-adiabatic corrections to phonon dispersion relations. Here we review the theoretical and computational framework underlying modern electron-phonon calculations from first principles, as well as landmark investigations of the electron-phonon interaction in real materials. In the first part of the article we summarize the elementary theory of electron-phonon interactions and their calculations based on density-functional theory. In the second part we discuss a general field-theoretic formulation of the electron-phonon problem, and establish the connection with practical first-principles calculations. In the third part we review a number of recent investigations of electron-phonon interactions in the areas of vibrational spectroscopy, photoelectron spectroscopy, optical spectroscopy, transport, and superconductivity.

The interaction between fermions and bosons is one of the cornerstones of many-particle physics. It is therefore unsurprising that, despite being one of the most thoroughly studied chapters of solid state physics, the interaction between electrons and phonons in solids continues to attract unrelenting attention.
Electron-phonon interactions (EPIs) are ubiquitous in condensed matter and materials physics. For example, they underpin the temperature dependence of the electrical resistivity in metals and the carrier mobility in semiconductors, they give rise to conventional superconductivity, and contribute to optical absorption in indirectgap semiconductors. In addition, EPIs enable the thermalization of hot carriers, determine the temperature dependence of electron energy bands in solids, and distort band structures and phonon dispersion relations of metals, leading to characteristic kinks and Kohn anomalies in photoemission and Raman/neutron spectra, respectively. EPIs also play a role in the areas of spintronics and quantum information, for example by coupling lattice and spin degrees of freedom in electromagnons, or by modulating the lifetimes of electron spins in color centers.
Given the fundamental and practical importance of electron-phonon interactions, it is perhaps surprising that the majority of theoretical studies in this area still rely on semi-empirical model Hamiltonians, especially in times when ab initio calculations have become pervasive in every area of condensed matter and materials physics. The reason for this lag can be found in the complexity of electron-phonon calculations: while density functional theory (DFT) calculations of total energies and structural properties were already well established in the early 1980s (Martin, 2004), systematic ab initio calculations of EPIs had to wait for the development of density functional perturbation theory (DFPT) for lattice dynamics between the late 1980s and the mid 1990s (Baroni et al., 1987;Gonze et al., 1992;Savrasov, 1992).
Despite this delayed start, the past two decades have witnessed tremendous progress in this area, and new exciting applications are becoming accessible as firstprinciples techniques for studying EPIs catch up with more established DFT methods. These advances are driving the evolution from qualitative and descriptive theories of electron-phonon effects in model solids to quantitative and predictive theories of real materials. As the methodology for calculating EPIs from first principles is rapidly reaching maturity, it appears that the time is ripe for reviewing this vast, complex and fascinating landscape.
One of the most authoritative reviews on the theory of EPIs is the classic book by Grimvall (1981). This monumental work represents an unmissable reference for the specialist. However, as this book pre-dates the rise of ab initio computational methods based on DFT, it inevitably misses the most recent developments in this area. The present article constitutes an attempt at filling this gap by reflecting on what DFT calculations can contribute to the study of electron-phonon physics. In addition, this article is also an opportunity to establish a unified conceptual and mathematical framework in this incredibly diverse landscape, shed light on the key approximations, and identify some of the challenges and opportunities ahead.
As emphasized by the title 'Electron-phonon interactions from first principles', the aim of this article is to review the ab initio theory of EPIs and to survey modern advances in ab initio calculations of EPIs. The reader interested in the fundamentals of electron-phonon physics or in theoretical developments relating to model Hamiltonians is referred to the outstanding monographs by Ziman (1960), Grimvall (1981), Schrieffer (1983), Mahan (1993), and Alexandrov and Devreese (2010).
Among significant recent advances that are covered in this review we mention the zero-point renormalization and the temperature dependence of electronic band structures; the calculation of phonon-assisted optical absorption spectra; the electron mass renormalization and the kinks in angle-resolved photoemission spectra; the thermalization of hot carriers in semiconductors; the calculation of phonon-limited mobility; the development of efficient computational techniques for calculating EPIs; and efforts to improve the predictive power of EPI calculations by going beyond standard density functional theory.
The review is organized as follows: Sec. II provides an historical perspective on the development of theories of the EPI, from early semi-empirical approaches to modern first-principles calculations. In Sec. III we ex-amine the various components of DFT calculations of EPIs in solids, and set the formalism which will be used throughout this article. Section IV provides a synthesis of the most advanced field-theoretic approaches employed to study EPIs, and Sec. V makes the link between the most general formalism and DFT calculations for real materials. In this section the reader will find a number of expressions which are useful for practical implementations. Section VI reviews advanced computational techniques for performing calculations of EPIs efficiently and accurately, such as Wannier interpolation and Fermi surface harmonics. Here we also discuss recent progress in the study of electron-phonon couplings in polar semiconductors. In Sec. VII we discuss recent calculations of phonons beyond the adiabatic Born-Oppenheimer approximation. Section VIII reviews calculations of EPIs in the context of photoelectron spectroscopy. Section IX focuses on the optical properties of semiconductors and insulators, in particular the temperature dependence of the band structure and phonon-assisted optical processes. In Sec. X we review calculations on the effects of EPIs on carrier dynamics and transport, including carrier thermalization rates and mobilities. Section XI discusses EPI calculations in the area of phonon-mediated superconductivity. Attempts at improving the accuracy and predictive power of ab initio EPI calculations by using more sophisticated electronic structure methods are discussed in Sec. XII. Finally in Sec. XIII we highlight the most pressing challenges in the study of EPIs from first principles, and we present our conclusions. We leave to the appendices some notational remarks and more technical discussions.

II. HISTORICAL DEVELOPMENT
The notion of 'electron-phonon interactions' is as old as the quantum theory of solids. In fact in the very same work where Bloch (1929) discussed the formal solutions of the Schrödinger equation in periodic potentials, Sec. V begins with the all-telling title: "The interaction of the electrons and the elastic waves of the lattice". In this work the first quantum theory of the temperaturedependent electrical resistivity of metals was developed. It took only a few years for Bloch's 'elastic waves' to be replaced by the brand-name 'phonon' by Frenkel (1932), thus establishing a tradition that continues unaltered almost a century later (Walker and Slack, 1970).
In order to discuss the early approaches to the electronphonon problem, it is useful to state right from the start the standard form of the Hamiltonian describing a cou-pled electron-phonon system: H = nk ε nkĉ † nkĉ nk + qν ω qν (â † qνâqν + 1/2) + N − 1 2 p k,q mnν g mnν (k, q)ĉ † mk+qĉ nk (â qν +â † −qν ) In this expression the first line describes the separate electron and phonon subsystems using the usual secondquantized formalism, while the second line specifies the mutual coupling between electrons and phonons to first order in the atomic displacements (Mahan, 1993). Here ε nk is the single-particle eigenvalue of an electron with crystal momentum k in the band n, ω qν is the frequency of a lattice vibration with crystal momentum q in the branch ν, andĉ † nk /ĉ nk (â † qν /â qν ) are the associated fermionic (bosonic) creation/destruction operators. N p is the number of unit cells in the Born-von Kármán supercell (see Appendix A). The third and fourth lines of Eq. (1) describe the electron-phonon coupling Hamiltonian to second order in the atomic displacements. This contribution is rarely found in the early literature (hence the square brackets), but it plays an important role in the theory of temperature-dependent band structures (Sec. V.B.1). The matrix elements g mnν (k, q) and g DW mnνν (k, q, q ) measure the strength of the coupling between the electron and the phonon subsystems, and have physical dimensions of an energy. Here the superscript 'DW' stands for Debye-Waller, and relates to the Debye-Waller self-energy to be discussed in Sec. V.B.2. Complete details as well as a derivation of Eq. (1) will be provided in Sec. III.
The formal simplicity of Eq. (1) conceals some important difficulties that one faces when attempting to use this equation for predictive calculations. For example, the electronic Hamiltonian relies on the assumption that the system under consideration can be described in terms of well-defined quasi-particle excitations. Similarly, the phonon term is meaningful only within the harmonic and the adiabatic approximations. More importantly, Eq. (1) does not provide us with any prescription for determining the numerical parameters ε nk , ω qν , g mnν (k, q), and g DW mnνν (k, q, q ).
In a sense the history of the study of electron-phonon interactions is really the history of how to calculate the parameters entering Eq. (1) using procedures that can be at once rigorous, reliable, and practical. As it will become clear in Sec. IV, despite enormous progress in this area, some conceptual difficulties still remain.
A clear account of the theory of EPIs until the late 1950s is given by Ziman (1960). In the following we highlight only those aspects that are relevant to the subsequent discussion in this article.
Early studies of electron-phonon interactions in solids were motivated by the quest for a quantum theory of the electrical resistivity in metals (Hoddeson and Baym, 1980). The common denominator of most early approaches is that the electronic excitations in Eq. (1) were described using the free electron gas model, ε nk = 2 k 2 /2m e − ε F , m e being the electron mass and ε F the Fermi energy; the lattice vibrations were described as acoustic waves using the Debye model, ω qν = v s |q|, v s being the speed of sound in the solid. Both approximations were reasonable given that the systems of interest included almost exclusively elemental metals, and primarily monovalent alkali and noble metals (Mott and Jones, 1936). While these approximations were fairly straightforward, it was considerably more challenging to determine the EPI matrix elements g mnν (k, q) using realistic approximations.
The very first expression of the electron-phonon matrix element was derived by Bloch (1929); using contemporary notation it can be written as: Here M κ is the mass of the κ-th nucleus, and e κν (q) is the polarization of the acoustic wave corresponding to the wavevector q and mode ν. The term V 0 represents a unit-cell average of the 'effective' potential experienced by the electrons in the crystal. Equation (2) was meant to describe the scattering from an initial electronic state with wavevector k to a final state with wavevector k + q, via an acoustic phonon of wavevector q and frequency ω qν . The formula was developed for free electron metals, and neglects so-called 'umklapp' (folding) processes, i.e. scattering events whereby k goes into k + q + G with G being a reciprocal lattice vector. A derivation of Eq.
(2) is provided in Sec. III.B.5. In order to determine V 0 Bloch (1929) argued that the crystal may be described as a continuous deformable medium. Starting from this assumption he reached the conclusion that the average potential can be approximated as V 0 = 2 /(16m e a 2 0 ) (a 0 is the Bohr radius). Even though Bloch's matrix element is no longer in use, this model provides helpful insight into the nature of EPIs in monovalent metals. For example the so-called 'polarization factor' in Eq. (2), q · e κν (q), shows that (in the absence of umklapp processes) only longitudinal sound waves scatter electrons. Nordheim (1931) proposed a refinement of Bloch's model whereby the average potential V 0 in Eq. (2) is replaced by the Fourier component V κ (q) of the ionic Coulomb potential (see Sec. III.B.5). The key assumption underlying this model is that the effective potential experienced by the electrons is simply the sum of the individual bare ionic potentials of each nucleus. When a nucleus is displaced from its equilibrium position, the corresponding potential also shifts rigidly. This is the so-called the 'rigid-ion' approximation.
The main difficulty that arises with the rigid-ion model is that the Fourier transform of the Coulomb potential diverges as q −2 for q = |q| → 0; this leads to unrealistically strong EPIs. In order to circumvent this difficulty Mott and Jones (1936) proposed to truncate the ionic potential at the boundary of the Wigner-Seitz unit cell of the crystal. This choice represents the first attempt at including the electronic screening of the nuclear potential in a rudimentary form. In practice Mott and Jones (1936) calculated the Fourier transform of V κ (r) by restricting the integration over a Wigner-Seitz cell; the resulting potential is no longer singular at long wavelengths. A detailed discussion of this model can be found in (Ziman, 1960).
Despite some initial successes in the study of the electrical conductivity of metals, the descriptive power of these early models was undermined by the complete neglect of the electronic response to the ionic displacements. The first attempt at describing the effect of the electronic screening was made by Bardeen (1937). In his model the average potential V 0 in Eq. (2) is replaced by: where (q) is the Lindhard function (Mahan, 1993): (q) = 1 + (k TF /q) 2 F (q/2k F ).
Here k TF and k F are the Thomas-Fermi screening wavevector and the Fermi wavevector, respectively, and F (x) = 1/2 + (4x) −1 (1 − x 2 ) log |1 + x|/|1 − x|. A derivation of Bardeen's model is provided in Sec. III.B.5. Since (q) → (k TF /q) 2 for q → 0, the sigularity of the electronnuclei potential is removed in Bardeen's matrix element. The work of Bardeen (1937) can be considered as a precursor of modern ab initio approaches, insofar the calculation of the matrix element was carried out using a self-consistent field method within the linearized Hartree theory. This strategy is similar in spirit to modern DFPT calculations.
The key qualitative difference between the approach of Bardeen (1937) and modern techniques lies in the neglect of exchange and correlation effects in the screening. A possible route to overcome this limitation was proposed by Bardeen and Pines (1955). In this work the authors considered the role of a screened exchange interaction in the electron-phonon problem (see Appendix B of their work), however the mathematical complexity of the formalism prevented further progress along this direction. Similar efforts were undertaken by Hone (1960), and a more detailed account of the early approximations to exchange and correlation can be found in (Grimvall, 1981).
The most interesting aspect of the work by Bardeen and Pines (1955), as well as previous work along the same lines by Nakajima (1954), is that for the first time the electron-phonon problem was addressed using a fieldtheoretic approach.
One interesting feature in the theory of Bardeen and Pines is that their field-theoretic formulation naturally leads to a retarded electron-phonon vertex: the effective potential experienced by electrons upon the displacement of nuclei depends on how fast this displacement takes place. In this approach the effective potential V 0 in Eq. (2) is replaced by the dynamically screened potential: Here (q, ω) is the frequency-dependent Lindhard function (Mahan, 1993), and the effect of electronic screening is evaluated at the phonon frequency, ω = ω qν . Somewhat surprisingly, this development was not followed up in the literature on ab initio calculations of EPIs.

Semiconductors
While the investigation of electron-phonon effects was initially restricted to monovalent metals, the formal developments were soon extended to the case of more complex systems such as semiconductors. Carriers in semiconductors are typically confined within a narrow energy range near the band extrema; consequently it is expected that the dominant electron-phonon scattering mechanisms will involve long-wavelength phonons (q → 0). This concept was formalized by Bardeen and Shockley (1950) and Shockley and Bardeen (1950), laying the foundations of the 'deformation-potential' method.
In the deformation potential approach it is assumed that the atomic displacements can be described by longwavelength acoustic waves, and these can be related in turn to the elastic strain of the crystal. Using concepts from the effective mass theory, Bardeen and Shockley showed that in this approximation the potential V 0 in Eq. (2) can be replaced by: where Ω represents the volume of the unit cell, and the electron eigenvalues correspond to the valence or conduction band extrema. The derivation of this result can be found in Appendix B of ( Bardeen and Shockley, 1950). The deformation potentials E 1 were obtained empirically; for example Bardeen and Shockley determined these values for the band extrema of silicon by fitting mobility data. More complex scenarios such as anisotropic constant-energy surfaces in semiconductors were subsequently addressed by considering the effects of shear deformations (Dumke, 1956). While the concept of deformation potentials has become a classic in semiconductor physics, this method relies on a semi-empirical approach and lacks predictive power.

Ionic crystals
A class of materials that played an important role in the development of the theory of EPIs is that of ionic crystals. The qualitative difference between ionic solids and the systems discussed in Secs. II.A.1-II.A.2 is that the atomic displacements can generate long-ranged electric fields; these fields provide a new scattering channel for electrons and holes.
The theory of polar electron-phonon coupling started with the investigation of the electron mean free path in ionic crystals, in search for a theoretical model of dielectric breakdown in insulators (Fröhlich, 1937;Fröhlich and Mott, 1939). The central idea of these models is that in insulators the density of free carriers is very low, therefore it is sensible to consider a single electron interacting with the polarization field of the ionic lattice.
The Fröhlich model is similar in spirit to the contemporary work of Bardeen (1937) for metals. The main difference is that Fröhlich considered the screening arising from the dielectric polarization of an insulating crystal, while Bardeen considered the screening arising from the response of the Fermi sea. Fröhlich et al. (1950) showed that in the case of isotropic ionic crystals the effective potential V 0 appearing in Eq. (2) must be replaced by: In this expression e is the electron charge, 0 is the dielectric permittivity of vacuum, 0 and ∞ are the static and the high-frequency relative permittivities, respectively. This result is derived in Sec. VI.A.3. Using Eqs. (7) and (2) we see that when 0 > ∞ the matrix element g mnν (k, q) diverges as |q| −1 at long wavelengths. This singular behavior can lead to very strong EPIs, and provides the physical basis for the phenomenon of electron self-trapping in polarons (Pekar, 1946;Emin, 2013). The initial studies in this area were rapidly followed by more refined approaches based on field-theoretic methods (Lee et al., 1953). A comprehensive discussion of the various models can be found in the original review article by Fröhlich (1954).

B. The pseudopotential method
The approximations underpinning the models discussed in Sec. II.A become inadequate when one tries to study EPIs for elements across the periodic table. This and other limitations stimulated the development of the pseudopotential method, starting in the late 1950s with the work of Phillips and Kleinman (1959). The theory of pseudopotentials is too vast to be summarized in a few lines, and the reader is referred to Chapter 11 of (Martin, 2004) for a thorough discussion. Here we only highlight the aspects that are relevant to the calculation of EPIs.
The genesis of the pseudopotential method is linked with the question on how the valence electrons of metals could be described using the electron gas model, even though the orthogonality to the core states imposes rapid fluctuations of the valence wavefunctions near the atomic cores. In order to address this question, it is useful to go through the key steps of the orthogonalized planewaves method (Herring, 1940). In this method one considers planewaves |k + G for the wavevector k + G, and projects out the component belonging the Hilbert subspace spanned by core electrons. This is done by defining |k + G OPW = |k + G − c |φ c φ c |k + G , where the |φ c represent the core states of all atoms in the system. The functions |k + G OPW are by construction orthogonal to core states, therefore they can be used to expand the valence electron wavefunctions |ψ nk using only a few basis elements: |ψ nk = G c k (G)|k + G OPW . In the language of pseudopotential theory |ψ nk is referred to as the 'all-electron' wavefunction, while the function |ψ nk = G c k (G)|k + G is referred to as the 'pseudo' wavefunction. The all-electron and the pseudo wavefunctions are simply related as follows: Here we used a modern notation borrowed from the projector-augmented wave (PAW) method of Blöchl (1994). By construction, the pseudo-wavefunction |ψ nk does not exhibit rapid fluctuations near the atomic cores. The projector operatorT is now used to rewrite the single-particle Schrödinger equation for the all-electron wavefunction (e.g. the Kohn-Sham equations) in terms of the pseudo-wavefunctions. UsingĤ|ψ nk = ε nk |ψ nk and Eq. (8) we have: which is a generalized eigenvalue problem. By replacing the definition ofT given above one finds (Phillips and Kleinman, 1959): withV rep = c (ε nk − ε c )|φ c φ c | and ε c being the eigenvalue of a core electron. Clearly the additional potentialV rep is strongly repulsive and is localized near the atomic cores. Cohen and Heine (1961) showed that this extra potential largely cancels the attractive potential of the nuclei. This is the reason why valence electrons in metals behave almost like free electrons.
The practical consequence of these developments is that it is possible to define smooth effective 'pseudopotentials' for systematic band structure calculations, whose form factors include only a few Fourier components (Phillips, 1958;Heine and Abarenkov, 1964;Animalu and Heine, 1965;Cohen and Bergstresser, 1966).
The use of pseudopotentials in electron-phonon calculations started with the works of Sham (1961) and Sham and Ziman (1963). Sham (1961) showed that, if the pseudopotential can be described by a local function, then the electron-phonon matrix element g mnν (k, q) can be calculated by replacing the all-electron potentials and wavefunctions by the corresponding pseudopotentials and pseudo-wavefunctions. In this approach the pseudo-potentials move around rigidly with the ionic cores, therefore we are dealing effectively with an improved version of the rigid-ion approximation discussed in Sec. II.A.
The pseudopotential method was employed by Shuey (1965) in order to calculate the electron-phonon matrix elements in germanium. Shortly afterwards many calculations of electron-phonon interactions based on the pseudopotential method appeared in the literature, including work on the resistivity of metals (Carbotte and Dynes, 1967;Dynes and Carbotte, 1968;Hayman and Carbotte, 1971;Kaveh and Wiser, 1972), the electron massenhancement in metals (Ashcroft and Wilkins, 1965;Grimvall, 1969;Allen and Cohen, 1970;Allen and Lee, 1972;Allen, 1972a) the superconducting transition temperatures within the McMillan formalism (Allen et al., 1968;Allen and Cohen, 1969), the mobility of semiconductors (Ralph, 1970), and the temperature dependence of semiconductor band structures Cardona, 1981, 1983). These calculations were mostly based on phonon dispersion relations extracted from neutron scattering data, and the results were in reasonable agreement with experiment. It seems fair to say that the pseudopotential method enabled the evolution from qualitative to quantitative calculations of electron-phonon interactions.
Before proceeding we note that, although Eqs. (8) and (9) were introduced starting from the method of orthogonalized planewaves, there exists considerable freedom in the choice of the operatorT . In practiceT can be chosen so as to makeψ nk as smooth as possible, while retaining information on the all-electron wavefunctions near the ionic cores. This was achieved by the PAW method of Blöchl (1994). Broadly speaking it is also possible to re-interpret the historical development of the pseudopotential method as the evolution of the projectorT . In fact Blöchl showed how the most popular pseudopotential methods (Hamann et al., 1979;Bachelet et al., 1982;Troullier and Martins, 1991;Vanderbilt, 1990) can be derived from the PAW method under specific approximations.
C. Ab initio self-consistent field calculations Predictive calculations of EPIs became possible with the development of ab initio DFT techniques. The key advantage of DFT methods is the possibility of calculating electron band structures, phonon dispersion relations, and electron-phonon matrix elements entirely from first principles. Historically, DFT started with the works of Hohenberg and Kohn (1964) and Kohn and Sham (1965). However, its widespread use had to wait for the development of accurate parametrizations of the exchange and correlation energy of the electron gas (Hedin and Lundqvist, 1971;von Barth and Hedin, 1972;Gunnarsson et al., 1974;Ceperley and Alder, 1980;Perdew and Zunger, 1981). An introduction to DFT techniques can be found in the books by Parr and Yang (1994, advanced), Martin (2004, intermediate), and Giustino (2014, elementary).
The first calculation of electron-phonon interactions using DFT was carried out by Dacorogna et al. (1985b) using a 'frozen-phonon' approach (see Sec. III.B.3). In this work the authors computed electron bands, phonon dispersions, and electron-phonon matrix elements of Al entirely from first principles. Quoting from the original manuscript: "This calculation is ab initio since only information about the Al atom, i.e. the atomic number and atomic mass, is used as input". Dacorogna et al. calculated the so-called electron-phonon coupling strength λ qν for several phonon branches ν and momenta q throughout the Brillouin zone, as well as the phonon linewidths arising from the EPI (see Secs. VII and XI.A). The average coupling strength was found to be in good agreement with that extracted from the superconducting transition temperature. In the approach of Chang et al. (1985); Dacorogna et al. (1985a,b); and Lam et al. (1986) the electron-phonon matrix element was calculated using: with u nk and u mk+q being the Bloch-periodic components of the Kohn-Sham electron wavefunctions, ∆ qν v KS being the phonon-induced variation of the self-consistent potential experienced by the electrons, and the integral extending over one unit cell. Equation (11) will be discussed in Sec. III.B.2. The scattering potential ∆ qν v KS was calculated by explicitly taking into account the rearrangement of the electronic charge following a small displacement of the nuclei. The inclusion of the selfconsistent response of the electrons constitutes a considerable step forward beyond the rigid-ion approximation of Sec. II.B.
The next and most recent step in the evolution of electron-phonon calculations came with the development of DFPT for lattice dynamics (Baroni et al., 1987;Gonze et al., 1992;Savrasov, 1992). In contrast to frozenphonon calculations, which may require large supercells, DFPT enables the calculations of vibrational frequencies and eigenmodes at arbitrary wavevectors in the Brillouin zone. This innovation was critical in the context of electron-phonon physics, since the calculation of many physical quantities requires the evaluation of nontrivial integrals over the Brillouin zone. The first calculations of EPIs using DPFT were reported by Savrasov et al. (1994), Liu and Quong (1996), Mauri et al. (1996), andBauer et al. (1998). They calculated the electrical resistivity, thermal conductivity, mass enhancement and superconducting critical temperature of a number of elemental metals (e.g. Al, Au, Cu, Mo, Nb, Pb, Pd, Ta, V, and Te), and reported good agreement with experiment.
By the late 1990s most of the basic ingredients required for the ab initio calculation of EPIs were available; subsequent studies focused on using these techniques for calculating a variety of materials properties, and on improving the efficiency and accuracy of the methodology. The most recent advances will be reviewed in Secs. VI-XII.

III. ELECTRON-PHONON INTERACTION IN DENSITY-FUNCTIONAL THEORY
In this section we review the basic formalism underlying the calculation of EPIs using DFT, and we establish the link with the Hamiltonian in Eq. (1). We start by introducing the standard formalism for lattice vibrations, and the electron-phonon coupling Hamiltonian. Then we briefly summarize established methods of DFPT for calculating electron-phonon matrix elements. For the time being we describe electrons and phonons as separate subsystems; a rigorous theoretical framework for addressing the coupled electron-phonon system will be discussed in Sec. IV.

A. Lattice vibrations in crystals
The formalism for studying lattice dynamics in crystals is covered in many excellent textbooks such as (Born and Huang, 1954;Ziman, 1960;Kittel, 1963;Ashcroft and Mermin, 1976;Kittel, 1976). Here we introduce the notation and summarize those aspects which will be useful for subsequent discussions in this section and in Secs. IV and V.
We consider M nuclei or ions in the unit cell. The position vector and Cartesian coordinates of the nucleus κ in the primitive unit cell are denoted by τ κ and τ κα , respectively. We describe the infinitely extended solid using Born-von Kármán (BvK) boundary conditions. In this approach, periodic boundary conditions are applied to a large supercell which contains N p unit cells, identified by the direct lattice vectors R p , with p = 1, . . . , N p . The position of the nucleus κ belonging to the unit cell p is indicated by τ κp = R p + τ κ . The Bloch wavevectors q are taken to define a uniform grid of N p points in one unit cell of the reciprocal lattice, and the vectors of the reciprocal lattice are indicated by G. In Appendix A we provide additional details on the notation, and we state the Fourier transforms between direct and reciprocal lattice.
Using standard DFT techniques it is possible to calculate the total potential energy of electrons and nuclei in the BvK supercell. This quantity is denoted as U ({τ κp }), where the braces are a short-hand notation for the coordinates of all the ions. The total potential energy refers to electrons in their ground state, with the nuclei being represented as classical particles clamped at the coordinates τ κp . Every DFT software package available today provides the quantity U as a standard output.
In order to study lattice vibrations, one begins by making the harmonic approximation. Accordingly, the total potential energy is expanded to second order in the displacements ∆τ καp of the ions in the BvK supercell away from their equilibrium positions τ 0 κp : where U 0 denotes the total energy calculated for the ions in their equilibrium positions, and the derivatives are evaluated for the equilibrium structure. The second derivatives of the total energy with respect to the nuclear coordinates define the matrix of 'interatomic force constants': The Fourier transform of the interatomic force constants yields the 'dynamical matrix' (Maradudin and Vosko, 1968): where M κ is the mass of the κ-th ion. The superscript 'dm' is there to distinguish this quantity from the manybody phonon propagators D(12) and D καp,κ α p that will be introduced in Sec. IV.B. The dynamical matrix is Hermitian and therefore admits real eigenvalues, which we denote as ω 2 qν : κ α D dm κα,κ α (q) e κ α ,ν (q) = ω 2 qν e κα,ν (q).
Here the index ν runs from 1 to 3M . The column vectors e κα,ν (q) for a given ν are called the 'normal modes of vibration' or the 'polarization' of the vibration wave. The following relations can be derived from Eq. (14): These relations between normal modes carry a degree of arbitrariness in the choice of phases; here we have chosen to follow the same phase convention as Maradudin and Vosko (1968). Using Eqs. (12) and (13) the Hamiltonian for nuclei considered as quantum particles can be written as: where the ground-state energy U 0 has been omitted and the last term is the kinetic energy operator. The Hamiltonian in the above expression corresponds to the energy of an entire BvK supercell. Equation (19) relies on two approximations: (i) the harmonic approximation, which coincides with the truncation of Eq. (12) to second order in the displacements; and (ii) the Born-Oppenheimer adiabatic approximation. This latter approximation is made when one calculates the interatomic force constants with the electrons in their ground state. The scope and validity of the adiabatic approximation will be discussed in detail in Sec. V.A.1. We note incidentally that, strictly speaking, the Born-Oppenheimer approximation does not need to be invoked were one to use the generalization of DFT to multicomponent systems introduced by Kreibich and Gross (2001) and Kreibich et al. (2008).
For practical purposes it is convenient to rewrite Eq. (19) by introducing the quanta of lattice vibrations. This is accomplished by defining the standard creation (â † qν ) and destruction (â qν ) operators for each phonon of energy ω qν and polarization e κα,ν (q). This operation is not entirely trivial and is described in detail in Appendix B. The formal definition of the ladder operators is given in Eqs. (B13)-(B14). These operators obey the commutation relations [â qν ,â † q ν ] = δ νν δ qq and [â qν ,â q ν ] = [â † qν ,â † q ν ] = 0, where δ is the Kronecker symbol. From these relations we know that the quanta of the harmonic oscillations in crystals obey Bose-Einstein statistics. In Appendix B it is shown that the atomic displacements can be expressed in terms of the ladder operators as follows: with l qν being the 'zero-point' displacement amplitude: Here M 0 is an arbitrary reference mass which is introduced to ensure that l qν has the dimensions of a length and is similar in magnitude to ∆τ καp . Typically M 0 is chosen to be the proton mass. Using Eqs. (13)-(21) the nuclear Hamiltonian can be written in terms of 3M N p independent harmonic oscillators as follows: where the sum is over all wavevectors. The groundstate wavefunction of this Hamiltonian is a product of Gaussians, and all other states can be generated by acting on the ground state with the operatorsâ † qν . In the case of |q| = 0 there are three normal modes for which ω qν = 0. For these modes, which correspond to global translations of the crystal, the zero-point displacement l qν is not defined. Throughout this article it is assumed that these modes are skipped in summations containing zero-point amplitudes. A detailed derivation of Eq. (22) and a discussion of the eigenstates ofĤ p are provided in Appendix B.

B. Electron-phonon coupling Hamiltonian
Having outlined the standard formalism for addressing lattice vibrations in crystals, we now proceed to make the connection between DFT calculations and the remaining terms of Eq. (1). The electronic band structure ε nk and electron-phonon matrix elements g mnν (k, q) are almost invariably calculated by using the Kohn-Sham (KS) Hamiltonian (Hohenberg and Kohn, 1964;Kohn and Sham, 1965). A justification for these choices will be provided in Sec. V; for now we limit ourselves to outline the key elements of practical calculations.

Kohn-Sham Hamiltonian
Let us denote the Kohn-Sham eigenfunctions by ψ nk (r), and use k to indicate both the wavevector and spin. We shall restrict ourselves to systems with collinear spins. The KS eigenfunctions satisfy the equa-tionĤ KS ψ nk (r) = ε nk ψ nk (r), with the Hamiltonian given by:Ĥ Here the potential V KS is the sum of the nuclear (or ionic) contribution V en , the Hartree electronic screening V H , and the exchange and correlation potential V xc (Martin, 2004): The potentials appearing in Eq. (24) are defined as follows. The electron-nuclei potential energy is given by: where V κ (r) is the interaction between an electron and the nucleus κ located at the center of the reference frame, and T denotes a lattice vector of the BvK supercell. In the case of all-electron DFT calculations, V κ (r) is the Coulomb interaction: where Z κ is the atomic number of the nucleus κ. In the case of pseudopotential implementations V κ is a function that goes as in Eq. (26) at large |r|, but remains finite at |r| = 0. Furthermore the nuclear charge is replaced by the ionic charge, that is the difference between the nuclear charge and the number of core electrons described by the pseudopotential. In all modern pseudopotential implementations V κ (r) is nonlocal due to the separation of the angular momentum channels (Martin, 2004). However, since this nonlocality is short-ranged and is inconsequential in the following discussion, it will be ignored here in order to maintain a light notation. The Hartree term is obtained from the electron density, n(r ; {τ καp }): where the integral extends over the supercell. The exchange and correlation potential is the functional derivative of the exchange and correlation energy with respect to the electron density (Kohn and Sham, 1965): The eigenfunctions ψ nk ofĤ KS can be expressed in the Bloch form: with u nk a lattice-periodic function. The wavefunction ψ nk is taken to be normalized to one in the supercell, while the periodic part u nk (r) is normalized to one in the crystal unit cell. The electron density is n(r) = vk |ψ vk (r)| 2 , where v indicates occupied states. In order to determine ψ nk and ε nk the Kohn-Sham equations are solved self-consistently. This requires one to start from a reasonable guess for the electron density (for example a superposition of atomic electron densities), calculate the potentials in Eq. (24), and determine the solutions of the KS Hamiltonian in Eq. (23). The electron density is re-calculated using these solutions, and the cycle is repeated until convergence. In order to establish the link with Eq. (1), we can regard the KS Hamiltonian as an effective one-body operator, and make the transition to a second-quantized formalism by using the standard prescription (Merzbacher, 1998): (30) This expression is useful for performing formal manipulations in the study of coupled electron-phonon systems. However, Eq. (30) implicitly introduces the drastic approximation that the electronic system can be described in terms of sharp quasiparticle excitations. A field-theoretic approach that does not rely on any such approximation is discussed in Sec. IV.

Electron-phonon coupling Hamiltonian to first-and second-order in the atomic displacements
Within the DFT Kohn-Sham formalism, the coupling Hamiltonian appearing in the second line of Eq. (1) is obtained by expanding the Kohn-Sham effective potential in terms of the nuclear displacements ∆τ κp from their equilibrium positions τ 0 κp . The potential to first order in the displacements is: This expression can be rewritten into normal mode coordinates using Eq. (20): having defined: From the last expression we see that ∂ κα,q v KS and ∆ qν v KS are lattice-periodic functions. The transition to second quantization is performed as in Eq. (30) (Merzbacher, 1998): where the brakets indicate an integral over the supercell. After using Eqs. (29), (32)-(35), and (A1) we have: where the electron-phonon matrix element is given by: Here the subscript 'uc' indicates that the integral is carried out within one unit cell of the crystal. The coupling Hamiltonian in Eq. (37) yields the energy of an entire supercell. In the case of the three translational modes at |q| = 0 we set the matrix elements g mnν (k, q) to zero, as a consequence of the acoustic sum rule (see discussion in Sec.IX.A.1). Taken together, Eqs. (22), (30), and (37) constitute the starting point of most first-principles calculations of electron-phonon interactions. It remains to be seen how one calculates the electron-phonon matrix elements g mnν (k, q); the most common procedures are described in Sec. III.B.3.
Before proceeding, we discuss briefly the second-order coupling Hamiltonian which appears in the third and fourth lines of Eq. (1). The rationale for incorporating this extra term is that the expansion of the Kohn-Sham potential to first order in the atomic displacements, Eq. (31), is somewhat inconsistent with the choice of expanding the total potential energy in Eq. (12) to second order in the atomic displacements. This aspect was discussed by Allen and Heine (1976) and Allen (1978). In order to obtain an electron-phonon coupling Hamiltonian including terms of second-order in the displacements, we must include the second derivatives of the , and follow the same steps which led to Eq. (37). By calling the extra termĤ (2) ep we have: where (40) The variations ∆ qν are the same as in Eqs. (33)-(35).
The second-order coupling Hamiltonian in Eq. (39) is considerably more involved than its first-order counterpart; the increased complexity partly explains why in the literature this term has largely been ignored. So far the HamiltonianĤ (2) ep has only been described using an approximation based on first-order perturbation theory (Allen and Heine, 1976). In this special case, the only terms in Eq. (39) that can modify the electron excitation spectrum are those with q = −q. The corresponding energy shift is ∆ε nk = N −1 p qν g DW nnνν (k, q, −q)(2n qν +1), with n qν being the number of phonons in each mode. We will come back to this point in Sec. IX.A.1.

Calculation of electron-phonon matrix elements using density-functional perturbation theory
In this section we review how the scattering potential ∆ qν v KS appearing in Eq. (38) is calculated in firstprinciples approaches. The most intuitive approach is to evaluate the derivatives appearing in Eq. (35) by using finite atomic displacements in a supercell: In this expression b is a small displacement of the order of the zero-point amplitude (∼ 0.1 Å), and the atom κ in the unit cell p is displaced along the direction α. The first calculations of electron-phonon interactions within DFT employed a variant of this 'supercell approach' whereby all atoms are displaced according to a chosen vibrational eigenmode (Chang et al., 1985;Dacorogna et al., 1985a,b;Lam et al., 1986); this strategy is usually referred to as the 'frozen-phonon' method. One disadvantage of the frozen-phonon method is that the supercell may become impractically large when evaluating matrix elements corresponding to long-wavelength phonons. This difficulty can be circumvented by using DFPT (Baroni et al., 1987;Gonze et al., 1992;Savrasov, 1992). The main strength of DFPT is that the scattering potential ∆ qν v KS in Eq. (38) is obtained by performing calculations within a single unit cell. Since the computational workload of standard (non linear-scaling) DFT calculations scales as the cube of the number of electrons, the saving afforded by DFPT over the frozenphonon method is proportional to N 2 p , and typically corresponds to a factor > 10 3 .
In the DFPT approach of Baroni et al. (2001) one calculates the lattice-periodic scattering potential ∂ κα,q v KS defined by Eq. (35). By differentiating Eq. (24) via Eq. (35) this potential is written as: The variation of the ionic potential is obtained from Eqs. (25) and (35). The result is conveniently expressed in reciprocal space: where the convention for the Fourier transform is f (G) = Ω −1 uc dr e −iG·r f (r), and Ω is the volume of the unit cell. In order to keep the presentation as general as possible we avoid indicating explicitly the non-locality of V κ which arises in pseudopotential implementations. The adaptation of this equation and the following ones to the case of nonlocal pseudopotentials, ultrasoft pseudopotentials, and the projector-augmented wave method can be found in (Giannozzi et al., 1991), (Dal Corso et al., 1997), and (Audouze et al., 2006), respectively. The variation of the Hartree and exchange-correlation contributions to the Kohn-Sham potential is obtained from the self-consistent charge density response to the perturbation in Eq. (43). After a few manipulations using Eqs. (27) and (35) one obtains: where v C (q) = Ω −1 dr e −iq·r e 2 /4π 0 |r| is the Fourier transform of the Coulomb potential. For the exchange and correlation potential we use Eq. (28) and the Taylor expansion of a functional to find: where f xc indicates the standard exchange and correlation kernel, which is the second-order functional derivative of the exchange and correlation energy E xc with respect to the electron density (Hohenberg and Kohn, 1964): In the case of the local density approximation (LDA) to DFT the exchange and correlation kernel reduces to a local function (Parr and Yang, 1994), and Eq. (45) is more conveniently evaluated in real space. Today DFPT calculations can be performed using one of several exchange and correlation kernels. The effect of the kernel on the calculation of lattice-dynamical properties of solids has been analyzed in several works, see for example Dal Corso (2013) and He et al. (2014). The formal structure of the DFPT equations discussed in this section remains unchanged if we replace the DFT kernel in Eq. (46) by more sophisticated versions. For example both DFPT calculations based on Hubbard-corrected DFT (Floris et al., 2011) and DFPT coupled with dynamical mean-field theory (Savrasov and Kotliar, 2003) have been demonstrated. It should be noted that in Eqs. (45)-(46) we are implicitly assuming a spin-unpolarized system. The adaptation of these equations as well as the other DFPT equations to the most general case of non-collinear spin systems can be found in (Dal Corso, 2007Verstraete et al., 2008).
From Eqs. (44) and (45) we see that the evaluation of g mnν (k, q) goes through the calculation of the variation of the electron density induced by the perturbation ∂ κα,q v KS (r) e iq·r . Within DFPT such a variation is obtained by evaluating the change of the Kohn-Sham wavefunctions to first order in perturbation theory. After inspection of the perturbed Hamiltonian it becomes evident that the wavefunction change must be of the form ∂u nk,q e iq·r , with ∂u nk,q a lattice-periodic function. Using this observation the first-order variation of the Kohn-Sham equations can be written as a Sternheimer equation (Sternheimer, 1954): withĤ KS k+q = e −i(k+q)·rĤ KS e i(k+q)·r . In this equation the index v indicates an occupied state. For |q| = 0 one needs also to consider a shift of the energy eigenvalues which introduces an additional term u vk |∂ κα,0 v|u vk uc u vk on the right-hand side of Eq. (47). In practice this term is canceled by the use of the projectors described in Eq. (48) below, unless one is dealing with metallic systems. This aspect is discussed in detail by de Gironcoli (1995) and Baroni et al. (2001). The principal advantage of Eq. (47) over standard perturbation theory is that it does not involve unoccupied electronic states.
A practical problem arises when attempting to solve Eq. (47): the linear system on the left-hand side is ill-conditioned owing to small eigenvalues corresponding to ε vk ε v k+q ; furthermore in the case of accidental degeneracies, ε vk = ε v k+q , the system becomes singular. In order to make the system non-singular Giannozzi et al. (1991) noted that the variation of the electron density only involves the component of ∂u vk,q belonging to the unoccupied manifold of Kohn-Sham states. As a consequence, what is really needed is only ∂ũ vk,q = (1 −P occ k+q ) ∂u vk,q , having denoted byP occ k+q = v |u vk+q u vk+q | the projector over the occupied states with wavevector k + q. The equation for this 'trimmed' wavefunction variation is simply obtained by projecting both side of Eq. (47) onto (1 −P occ k+q ), and noting that P occ k+q andĤ KS k+q do commute: At this point it is possible to remove all small or null eigenvalues of the operator on the left-hand side by adding a term αP occ k+q to the Hamiltonian. This term has no effect on the wavefunction variation, sincê P occ k+q ∂ũ vk,q = 0 by construction. The operator is made non-singular by choosing the parameter α larger than the valence bandwidth (Baroni et al., 2001). From the wavefunction variation obtained by solving Eq. (48), it is now possible to construct the density response associated with the wavevector q: For simplicity a spin-degenerate system has been assumed (a factor of 2 is implicitly included in the sum over k), and time-reversal symmetry has been used in order to make the expression more compact (yielding the factor of 2 on the right-hand side). In practical DFPT calculations, Eq. (48) is solved using an iterative procedure which is similar to standard DFT total energy calculations. One sets the starting perturbation ∂ κα,q v KS to be equal to the electron-nuclei potential in Eq. (43). By solving Eq. (48) for each occupied state v and each wavevector k using standard linear algebra techniques, one obtains the induced density in Eq. (49). The new density is now used to construct the variations of the Hartree and exchange-correlation potentials in Eqs. (44) and (45). These induced potentials are added to the electron-nuclei potential, yielding a 'screened' perturbation ∂ κα,q v KS in Eq. (48). The cycle is repeated until the change of ∂n κα,q between two successive cycles is smaller than a set tolerance.
It can be shown that the screened perturbation ∂ κα,q v KS described in this section is also the key ingredient required for calculating the interatomic force constants in Eq. (13) (Baroni et al., 2001). As a practical consequence, every software implementation that supports DFPT calculations already contains all the information necessary for evaluating the electron-phonon matrix elements g mnν (k, q).
All the quantities introduced in this section can equivalently be calculated using an alternative, variational formulation of density-functional perturbation theory (Gonze et al., 1992;Gonze, 1995aGonze, , 1997Gonze and Lee, 1997). A thorough discussion of the connection between the Sternheimer approach and the variational approach to DFPT is provided by Gonze (1995b).
The second-order matrix elements g DW mn,νν (k, q, q ) given by Eq. (40) involve the second derivative of the Kohn-Sham potential with respect to the nuclear displacements. The evaluation of these quantities would require the solution of second-order Sternheimer equations for the second variations of the Kohn-Sham wavefunctions. The general structure of second-order Sternheimer equations can be found in Sec. IV.H of (Gonze, 1995b). Since these calculations are rather involved, most practical implementations employ an approximation whereby the Debye-Waller matrix elements are expressed in terms of products of the standard matrix elements g mnν (k, q). Such an alternative formulation was developed by Allen and Heine (1976) and Allen and Cardona (1981), and will be discussed in Sec. IX.A.1. All recent ab initio calculations of electron-phonon interactions based on DFPT employed this latter approach. 1

The dielectric approach
Besides the DFPT method described in the previous section, it is also possible to calculate the screened perturbation ∂ κα,q v KS using the so-called 'dielectric approach' (Pick et al., 1970;Quong and Klein, 1992). This latter approach did not find as widespread an application as those of Baroni et al. (1987); Gonze et al. (1992);and Savrasov (1992), but it is useful to establish a link between DFT calculations of electron-phonon matrix elements and the field-theoretic formulation to be discussed in Sec. IV.
For consistency with Sec. III.B.3, we derive the key expressions of the dielectric approach starting from DFPT. To this aim we expand the variation of the wavefunction ∂ũ vk,q using the complete set of states u nk+q (with n referring to both occupied and empty Kohn-Sham states). Then we replace this expansion inside Eq. (48), project onto an arbitrary conduction state, and insert the result in Eq. (49). After taking into account time-reversal symmetry, these steps lead to the following result: having defined: In this expression f nk and f mk+q are the occupations of each state, and the indices run over all bands. A factor of 2 for the spin degeneracy is implicitly included in the sum over k. The quantity χ 0 q in Eq. (51) is the lattice-periodic component for the wavevector q of the 'independent-electron polarizability' (Adler, 1962;Wiser, 1963;Pick et al., 1970;Quong and Klein, 1992).
For ease of notation we can write Eq. (50) in symbolic form as ∂n = χ 0 ∂v KS . Using the same symbolic notation it is also possible to formally rewrite Eqs. (42), (44), and (45) as follows: from which one obtains: having defined the dielectric matrix: The superscript 'Hxc' refers to the Hartree and exchange and correlation components of the screening. In the language of many-body perturbation theory Hxc is referred to as the 'test electron' dielectric matrix, hinting at the fact that the electron density redistribution in response to a perturbation arises both from classical electrostatics (the Hartree term v C χ 0 ) and from quantum effects (the exchange and correlation term f xc χ 0 ). If we neglect the kernel f xc in this expression, then we obtain the 'test charge' dielectric matrix, which is most commonly known as the dielectric matrix in the random-phase approximation (RPA) (Pines and Bohm, 1952): The symbolic expressions outlined here remain almost unchanged when using a reciprocal-space representation.
As an example, Eq. (55) becomes simply: Taken together Eqs. (38) and (53) show that the calculation of electron-phonon matrix elements using DFPT is equivalent to screening the bare electron-nucleus interaction using Hxc ; in this case we say that the screening is described at the 'RPA+xc' level of approximation. At this point it is worth to point out that so far we only considered the screening of static perturbations: in fact ∂v en was implicitly taken to be frequency-independent. Physically this choice corresponds to describing phonons as quasi-static perturbations, so that at each set of instantaneous atomic positions during a vibration cycle, the electrons have enough time to re-adjust and reach their ground state. This is a statement of the adiabatic approximation (Born and Oppenheimer, 1927). The importance of retardation effects in the electron-phonon problem was already recognized in the early work of Bardeen and Pines (1955), but the first ab initio calculations of these effects appeared much later (see . The formal framework required to incorporate retardation in the study of EPIs will be presented in Sec. IV.

Connection with early formulations
For completeness, we illustrate the link between electron-phonon matrix elements obtained within DFPT (Sec. III.B.3) and the early approaches of Bloch (1929) and Bardeen (1937) The Bloch matrix element can be derived as follows. We assume that the scattering potential is unscreened and corresponds to the bare pseudopotentials V κ in Eq. (43); that there is only one atom at the origin of the unit cell; and the Kohn-Sham wavefunctions can be approximated by free electrons, u nk (r) = Ω − 1 2 exp(iG n · r). In the last expression, the subscript in G n is used in order to stress the one-to-one correspondence between the reciprocal lattice vectors and the energy bands of the free electron gas in the reduced zone scheme. Using these ap-proximations in Eqs. (21), (34), (38), and (43), we find: By further neglecting umklapp processes (G m = G n ) the previous result becomes (Grimvall, 1981, Sec. 3.4): The expression obtained by Bloch (1929) and reproduced in Eq.
(2) is simply obtained by replacing V κ (q) with the effective potential V 0 . The Bardeen matrix element is more elaborate and can be derived as follows. We describe the screening of the bare ionic potential within the RPA approximation, and determine the dielectric matrix by replacing the Kohn-Sham wavefunctions by free electrons. Using u nk (r) = Ω − 1 2 exp(iG n · r) and ε nk = 2 (k + G n ) 2 /2m e − ε F in Eq. (51), the polarizability reduces to: where F is the function defined below Eq. (4). The derivation of this result requires making the transition from the first Brillouin zone to the extended zone scheme. If we use Eq. (59) inside Eq. (56), neglect the exchange and correlation kernel, and use the Fourier transform of the Coulomb potential, we find: where the Thomas-Fermi screening length is given by k TF = [4me 2 k F /(4πε 0 π 2 )] 1/2 . Equation (60) is the wellknown Lindhard dielectric matrix, and the diagonal matrix elements are the same as in Eq. (4) (see Mahan, 1993, andGiuliani andVignale, 2005, for in-depth discussions of the Lindhard function). By following the same steps that led to Eq. (58), replacing the bare ionic potential by its screened counterpart, and using Eq. (53) with instead of Hxc , we obtain: (61) Here we considered only one atom at the center of the unit cell, and we neglected umklapp processes. This is essentially the result derived by Bardeen (1937) and reproduced in Eq. (3).

IV. FIELD-THEORETIC APPROACH TO THE ELECTRON-PHONON INTERACTION
In Sec. III we discussed how the materials parameters entering the electron-phonon Hamiltonian in Eq. (1), namely ε nk , ω qν , and g mnν (k, q), can be calculated from first principles using DFT and DFPT. Today the formalism and techniques described in Sec. III constitute de facto the standard tool in quantitative studies of electronphonon interactions in solids (see Secs. VII-XII).
However, it should be noted that the DFT approach to EPIs does not rest on strong theoretical foundations. For one, the evaluation of the EPI matrix elements via Eq. (38) relies on the assumption that the interaction between electrons and nuclei is governed by the effective Kohn-Sham potential; therefore we can expect the matrix elements to be sensitive to the exchange and correlation functional (see Sec. XII). Furthermore, the very definition of phonons starting from Eq. (12) relies on the Born-Oppenheimer approximation, and one might ask whether this choice is accurate enough in metals and narrow-gap semiconductors (see Sec. VII). Finally, if one were to go beyond the Born-Oppenheimer approximation, then it would seem sensible to also incorporate retardation effects in the calculation of the EPI matrix elements.
On top of these practical points, and at a more fundamental level, we expect that the electron-phonon interaction will modify both the electronic structure and the lattice dynamics of a solid, and these modifications will in turn affect the coupling between electrons and phonons. It is therefore clear that a complete theory of interacting electrons and phonons must be self-consistent. In order to address these issues it is necessary to formulate the electron-phonon problem using a rigorous and general theory of interacting electrons and phonons in solids.
The most systematic and elegant approach is based on quantum field theory (Schwinger, 1951), and is tightly connected with the development of the GW method (Hedin, 1965). The first attempts in this direction were from Nakajima (1954), Bardeen and Pines (1955), Migdal (1958), and Engelsberg and Schrieffer (1963). However, from the point of view of the present article, these works are of limited usefulness since they were mostly developed around the homogeneous electron gas.
A completely general formulation of the problem, which seamlessly applies to metals, semiconductors, and insulators, was first provided by Baym (1961) and subsequently by Hedin and Lundqvist (1969). The formalism developed in these articles constitutes today the most complete theory of the electron-phonon problem. In fact, many aspects of this formalism are yet to be explored within the context of ab initio calculations. After these seminal works several authors contributed to clarifying various aspects of the many-body theory of the coupled electron-phonon system, including Keating (1968), Gillis (1970), Sjölander and Johnson (1965), Maksimov (1976), Vogl (1976), and more recently van Leeuwen (2004) and Marini et al. (2015). In particular, van Leeuwen focused on the issues of translational and rotational invariance of the resulting theory, while Marini et al. analyzed the connection between many-body perturbation theory approaches and DFT calculations.
Since the mathematical notation of the original articles is obsolete and rather difficult to follow, in Secs. IV.A-IV.D we cover the theory in some detail using contemporary notation. The following derivations can be found across the works of Kato et al. (1960), Baym (1961), Hedin and Lundqvist (1969), and Maksimov (1976). Here we provide a synthesis of these contributions using a unified notation, and we fill the gaps wherever it is necessary. The presentation requires some familiarity with field operators (see for example Merzbacher, 1998 for a succinct introduction).

A. Operators and distinguishability
The starting point for studying EPIs using a fieldtheoretic approach is to define the Fock space and the field operators for electrons and nuclei. In the case of electrons the choice is unambiguous, since any manybody state can be represented as a linear combination of Slater determinants constructed using a basis of singleparticle wavefunctions. In the case of nuclei the situation is slightly more ambiguous: in principle we might proceed in a very general way by choosing to focus on the nuclei as our quantum particles, as opposed to their displacements from equilibrium. In practice this choice leads to a dead end for two reasons. Firstly, the quantum statistics of nuclei would be dependent on their spin, therefore we would end up with an unwieldy mix of fermions and bosons depending on the solid. Secondly, the notion of 'indistinguishable' particles, which is central to second quantization, does not apply to nuclei in solids (at least in thermodynamic equilibrium and far from a solid-liquid phase transition). In fact, in many cases we can directly label the nuclei, for example by means of experimental probes such as scanning tunneling microscopy and electron diffraction. In order to avoid these issues, it is best to study the electron-phonon problem by considering (i) indistinguishable electrons, for which it is convenient to use second-quantized operators; (ii) distinguishable nuclei, for which it is best to use first quantization in the displacements; (iii) indistinguishable phonons, resulting from the quantization of the nuclear displacements; in this latter case the distinction between first and second quantization is irrelevant. These aspects are briefly mentioned by Baym (1961) and Maksimov (1976).
With these choices, the dynamical variables of the problem are the electronic field operatorsψ (discussed below) and the nuclear displacements from equilibrium ∆τ (discussed in Sec. IV.C). In this theory the equilibrium coordinates of the nuclei are regarded as external parameters, and are to be obtained for example from crystallography or DFT calculations. Throughout this section, we limit ourselves to consider equilibrium Green's functions at zero temperature. As a result, all expectation values will be evaluated for the electron-nuclei ground state |0 . The extension of the main results to finite temperature is presented in Sec. V. We will not specify how to obtain the ground state, since the following discussion is independent on the precise shape of this state. In order to derive expressions that are useful for first-principles calculations, at the very end the ground state will be approximated using standard DFT wavefunctions and phonons (see Sec. V).
The electronic field creation/destruction operators are denoted byψ † (x)/ψ(x), where the variable x indicates both the position r and the spin label σ. These operators obey the anti-commutation relations (Merzbacher, 1998): . The most general non-relativistic Hamiltonian for a system of coupled electrons and nuclei can be written as:Ĥ where each term will be introduced hereafter. The electron kinetic energy is: with m e being the electron mass, and the integrals dx denoting the sum over spin and the integration over space, σ dr. The electron-electron interaction is: where the electron particle density operator is given bŷ n e (r) = σψ † (x)ψ(x), and where v(r, r ) = e 2 /(4πε 0 |r− r |) is the Coulomb interaction between two particles of charge e. In Eqs. (63) and (64) the integrals are over the entire crystal. This corresponds to considering a supercell of infinite size (therefore the lattice vectors T of the supercell drop out) and a dense sampling of wavevectors q in the Brillouin zone. This choice is useful in order to maintain the formalism as light as possible. Accordingly, all sums over q are replaced using The nuclear kinetic energy operator is the same as the last term in Eq. (19). Using the same notation as in Sec. III the nucleus-nucleus interaction energy is: Here τ 0 κp denotes the classical equilibrium position of each nucleus, and the displacement operators ∆τ κp will later be expressed in terms of the ladder operators from Appendix B. The electron-nucleus interaction energy is: where the nuclear charge density operator is given by: Here the density operators are expressed in units of the electron charge, so that the expectation value of the total charge density is −e 0|n(r)|0 withn(r) =n e (r) +n n (r).
We underline the asymmetry between Eqs. (64) and (65): in the case of electrons one considers the electrostatic energy of a continuous distribution of charge, and the unphysical self-interaction is removed by the Dirac delta; whereas in the case of nuclei, the particles are distinguishable therefore one has to take into account all pairwise interactions individually.

Equation of motion and self-energy
In this section we focus on the electrons. By combining Eqs. (62)-(66) and using the anti-commutation relations for the field operators one finds the standard expression: where the nuclear potentialV n is given by: In order to study the excitation spectrum of the manybody HamiltonianĤ at equilibrium we need to determine the time-ordered one-electron Green's function (Kato et al., 1960;Fetter and Walecka, 2003). At zero temperature this function is defined as: whereT is Wick's time-ordering operator for fermions, and ensures that the times of the subsequent operators increase towards the left. The formal definition of the Wick operator is: Based on this definition we see that for t > t the Green's function in Eq. (70) corresponds to the scalar product between the initial state ψ † (x t )|0 and the final state ψ † (xt)|0 . This product is precisely the probability amplitude for finding an electron in the position x at the time t, after having introduced an electron in x at an earlier time t . In the case t < t the situation is reversed and the Green's function describes the propagation of a hole created in the system at the time t . In order to determine G(xt, x t ) we need to establish an equation of motion for the field operators. This can be done by describing the time-dependence of the operators within the Heisenberg picture: whereĤ was defined in Eq. (68). From this definition it follows immediately: By combining Eqs. (68) and (72) and using the anticommutation relations for the field operators one obtains: where the time-dependence inn(r t) is to be understood in the Heisenberg sense, as in Eq. (71). This equation of motion allows us to write the corresponding equation for the electron Green's function in Eq. (70): Here , the brakets · · · are a short-hand notation for 0| · · · |0 , and the additional term ϕ is discussed below. In order to obtain Eq. (74) we used once again the anti-commutation relations, and we noted that the derivative of the Heaviside function is a Dirac delta.
The new term ϕ(rt) which appeared in Eq. (74) is a scalar electric potential which couples to both electronic and nuclear charges. This potential has been introduced in order to perturb the system via the additional Hamil-tonianĤ 1 (t) = drn(rt)ϕ(rt). The physical idea behind this choice is to use ϕ(rt) in order to induce forced oscillations in the system. When the system resonates with the perturbation we know that the resonant frequency must correspond to a free oscillation, that is a many-body eigenmode. From a formal point of view, the potential ϕ(rt) is introduced in order to exploit Schwinger's functional derivative technique (Kato, 1960, Appendix II) and is set to zero at the end of the derivation.
One complication arising from the introduction of ϕ(rt) in Eq. (74) is that the time evolution in Eq. (71) is no longer valid, since the perturbed Hamiltonian now depends on the time variable. The way around this complication is to switch from the Heisenberg picture to the interaction picture. This change amounts to replacing the exponentials in Eq. (71) by the time-ordered (Fetter and Walecka, 2003, p. 57). Since this would lead to an overlong derivation, we prefer to leave this aspect aside and refer the reader to Aryasetiawan and Gunnarsson (1998) for a more comprehensive discussion.
In order to write Eq. (74) in a manageable form, we use the identity (Kato et al., 1960): (75) In this and the following expressions δ/δϕ(r t ) denotes the functional derivative with respect to ϕ(r t ), and should not be confused with the Dirac delta functions δ(xt, x t ). Equation (75) is proven by Kato et al. (1960, Appendix II) and by Hedin and Lundqvist (1969, Appendix B.a). After identifyingâ andb withψ(xt) and ψ † (x t ), respectively, Eq. (74) becomes: where η is a positive infinitesimal arising from the timeordering, and V tot (rt) is the total potential acting on the electronic and nuclear charges, averaged over the manybody quantum state |0 : Equation (76) was first derived by Kato et al. (1960). In order to avoid a proliferation of variables, it is common practice to replace the letters by integer numbers, using the convention (xt) or (rt) → 1, (x t ) or (r t ) → 2, (rt+η) → 1 + , and so on. Using this convention the last two equations become: and In these expressions the spin labels are implied for the Green's function and for the Dirac delta. At this point, a set of self-consistent equations for coupled electrons and phonons can be generated by eliminating the functional derivative in Eq. (78). For this purpose one first relates the total screened electrostatic potential V tot to the external potential ϕ by introducing the inverse dielectric matrix −1 as a functional derivative: The function −1 (12) is the many-body counterpart of the dielectric matrix discussed in Sec. III.B.4. The form given by Eq. (80) is the most general field-theoretic formulation for a system of interacting electrons and nuclei.
The next step is to rewrite δG/δϕ inside Eq. (78) in terms of the inverse Green's function, G −1 . By using the fact that δ d2 G(12) G −1 (23) = 0 and the rule for the functional derivative of a product (Kadanoff and Baym, 1962) one obtains: In order to eliminate any explicit reference to ϕ we can express the functional derivative on the right-hand side using the chain rule for functional differentiation (Kadanoff and Baym, 1962): It is customary to call 'vertex' the three-point quantity defined by: By combining Eqs. (78) and (80)-(83) one finds: having introduced the so-called 'electron self-energy' Σ: which in turn contains the 'screened Coulomb interaction' W , defined as: The last equality can be obtained by observing that δ n(1) /δϕ(2) = δ n(2) /δϕ(1) after Eq. (75), therefore W (12) = W (21). Now, by inverting Eq. (84) and using Eq. (83), we can express the vertex Γ in terms of Σ and G: (87) The derivation of this result is rather lengthy: it requires the use of the chain rule, in symbols δ/δV tot = (δG/δV tot )δ/δG, as well as Eq. (81) with Σ and V tot instead of G and ϕ, respectively.
Equations (84)-(87) form a nonlinear system of equations for the electron Green's function, G, the electron self-energy, Σ, the total screened Coulomb interaction, W , and the vertex, Γ. In order to close the loop it remains to specify the relation between W and the other quantities. The next section is devoted to this aspect.

The screened Coulomb interaction
An equation for the screened Coulomb interaction can be found by combining Eqs. (79), (80), and (86): By defining the 'polarization propagator' as: the previous expression takes the usual form (Hedin, 1965): This result can be combined with Eq. (86) in order to express the dielectric matrix in terms of the polarization: We now consider the special case whereby the nuclei are regarded as classical point charges clamped to their equilibrium positions. In this situation, the variation of the charge density δ n in Eq. (89) corresponds to the redistribution of the electronic charge in response to the perturbation δV tot . In order to describe this special case it is convenient to introduce a new polarization propagator, P e , associated with the electronic response only: (92) The last equality in this expression is obtained by using Eq. (81) with V tot instead of ϕ, together with Eq. (83), and by considering that the electron density is related to the Green's function via the relation: In conjunction with P e it is natural to define the Coulomb interaction screened by the electronic polarization only: as well as the associated dielectric matrix, in analogy with Eq. (91): Taken together, Eqs. (84)-(87) with W replaced by W e constitute the well-known Hedin's equations for a system of interacting electrons when the nuclei are clamped to their equilibrium positions (Hedin, 1965).
In order to go back to the most general case whereby the nuclei are not clamped to their equilibrium positions, one has to describe the re-adjustment of both electronic and nuclear charge. To this aim we combine Eq. (80), (86), (88), (92), (94), and (95). The result is: An explicit expression for δ n n /δϕ can be obtained using the following reasoning. We go into the details since this is a delicate passage. Equation (75) provides a recipe for evaluating the variation of any operator with respect to a potential ϕ(rt) which couples to the total charge density operatorn(rt) viaĤ 1 (t) = drn(rt)ϕ(rt). Therefore we can replaceâb in Eq. (75) byn n to obtain: In addition, if we introduce a second perturbation, H 2 (t) = drn n (rt)J(rt), which couples only to the nuclear charges, we can repeat the same reasoning as in Eq. (97) after replacing ϕ by J andn byn n : The comparison between Eqs. (97) and (98) yields: This can be restated by using the chain rule, δ n e /δJ = δ n e /δV tot × δV tot /δ n × δ n /δJ: The variation δ n n /δJ on the right-hand side can be expressed as in Eq. (98): and since n n − n n = 0 this can also be rewritten as: having defined: (103) This quantity is called the 'density-density correlation function' for the nuclei.
Finally, we can combine Eqs. (96), (99), (100), and (103) to obtain: where W ph is the nuclear contribution to the screened Coulomb interaction, and is given by: This important result was derived first by Hedin and Lundqvist (1969).

Nuclear contribution to the screened Coulomb interaction
In view of the forthcoming discussion, it is useful to derive a more explicit expression for the screened interaction W ph in Eq. (105). Here we follow Baym (1961) and Maksimov (1976). The Taylor expansion of the Dirac delta to second order in the displacement u reads: where u · ∇∇ · u is a short-hand notation for the secondorder derivative αα u α u α ∇ α ∇ α . The above expression derives from the Fourier representation of the Dirac delta. Using Eq. (106) inside (67) we deduce: where n 0 n (r) is the density of nuclear point charges at the classical equilibrium positions τ 0 κp . After taking into account the time evolution in the Heisenberg picture as in Eq. (71), we can replace this expansion inside Eq. (103) to obtain: On the right-hand side we introduced the 'displacementdisplacement correlation function': If we insert the last two equations in Eq. (105) we find: In this expression V κ is the bare Coulomb potential of a nucleus or its ionic pseudo-potential. At this point of the derivation, Hedin and Lundqvist introduce the approximation that the electronic dielectric matrix in Eq. (110) can be replaced by its static counterpart. This choice implies the Born-Oppenheimer adiabatic approximation. In view of maintaining the formalism as general as possible, we prefer to keep retardation effects, following the earlier works by Bardeen and Pines (1955) and Baym (1961). We will come back to this aspect in Secs. V.A.3 and V.B.2.
We stress that the sole approximation used until this point is to truncate the density operator for the nuclei to the second order in the atomic displacements. This is the standard harmonic approximation. Apart from this approximation, which is useful to express W ph in a tractable form, no other assumptions are made. Gillis (1970) proposed a generalization of the results by Baym (1961) which does not use the harmonic approximation. However, the resulting formalism is exceedingly complex, and has not been followed up.

C. Phonon Green's function
In order to complete the set of self-consistent manybody equations for the coupled electron-phonon system, it remains to specify a prescription for calculating the displacement-displacement correlation function, D καp,κ α p (tt ). This function is seldom referred to as the 'phonon Green's function', even though stricly speaking this name should be reserved for the quantity −(i/ ) Tâ qν (t)â † q ν (t ) which will be discussed in Sec. V.A.1. In the following we describe the procedure originally devised by Baym (1961), and subsequently analyzed by Keating (1968), Hedin and Lundqvist (1969), Gillis (1970), and Maksimov (1976).
The starting point is the equation of motion for the displacement operators ∆τ καp (t). In analogy with Eq. (72) we have: i ∂/∂t∆τ κp (t) = [∆τ κp (t),Ĥ]. Since we are considering the harmonic approximation and we expect the nuclei to oscillate around their equilibrium positions, it is convenient to aim for an expression resembling Newton's equation. This can be done by taking the timederivative of the equation of motion: After evaluating the commutators using Eqs. (62)-(66) and performing the derivatives with respect to the nuclear displacements by means of Eq. (106), we obtain (the steps are laborious but straightforward): Heren (κp) (r) is the total charge density of electrons and nuclei, except for the contribution of the nucleus κ in the unit cell p. In the second line ∇ indicates that the derivatives are taken with respect to the variable r . At this point, we can use the functional derivative technique as in Sec. IV.B in order to determine an expression involving the displacement-displacement correlation function from Eq. (109). Here, instead of using J(r) as in Sec. IV.B.2 for the nuclear density, it is convenient to work with the individual displacements, and introduce a third pertur-bationĤ 3 (t) = κp F κp (t) · ∆τ κp (t). The extra terms F κp (t) have the meaning of external forces acting on the nuclei. Using this perturbation, it is possible to write the displacement-displacement correlation function in a manner similar to Eq. (101): This result was derived by Baym (1961) using a finitetemperature formalism. As in the case of the electron Green's function in Sec. IV.B, it can only be obtained by working in the interaction picture, and by taking into account the explicit time-dependence of the Hamiltonian H +Ĥ 3 (t). Also in the present case, we omit these details for the sake of conciseness.
If we take the expectation value of Eq. (112) in the ground state, after having added the new force term −F κp (t), and carry out the functional derivative with respect to F κ p (t ), we obtain: where the sum over the Cartesian directions γ is implied. The derivation of this result is rather cumbersome and involves the following considerations: (i) the Dirac deltas in the first line come from the force terms −F κp added to the Hamiltonian; (ii) within the Harmonic approximation the expectation value n (κp) ∆τ κp can be replaced by n (κp) ∆τ κp (Baym, 1961); (iii) the expectation value ∆τ κp can be set to zero, because at the end of the derivation one sets |F κp | = 0 hence the expectation values of the displacements vanish. We note that Hedin and Lundqvist (1969) omitted the last line of Eq. (114) in their derivation, but this term was correctly included by Baym (1961) and Maksimov (1976). The remaining functional derivative in Eq. (114) can be expressed in terms of the nuclear charge density using the same strategy which led to Eq. (100). The result is: By inserting this result inside Eq. (114) and using the expansion in Eq. (107), we finally obtain the equation of motion for the displacement-displacement correlation function: The quantity Π κpα,κ p α (tt ) in this expression is called the 'phonon self-energy' and is given by: The derivation of Eq. (117) is nontrivial and is not found consistently in the literature; it requires converting the derivatives with respect to the position variables r, r into derivatives with respect to the nuclear coordinates; integrating by parts in order to re-arrange the derivatives with respect to r and r ; invoking the harmonic approximation; and considering that, after setting the forces F kp (t) = 0 and the field ϕ(rt) = 0 at the end, the expectation value n(rt) does not depend on time. The term in the third line of Eq. (117) is what Baym (1961) called the 'static force', since it arises from the forces experienced by the nuclei in their equilibrium configuration. In order to simplify Eq. (117) it is convenient to move from the time to the frequency domain. We use the following convention for the Fourier transform of a function Since we are considering equilibrium properties in absence of time-dependent external potentials, the time variables enter in the above quantities only as differences (Abrikosov et al., 1975); for example W e (rt, r t ) = W e (r, r , t−t ). As a consequence, Eq. (116) is rewritten as: whereas the phonon self-energy in the frequency-domain at equilibrium reads: The second line in this expression is conveniently rewritten by making use of the acoustic sum rule. This sum rule is well-known in the theory of lattice dynamics of crystals (Born and Huang, 1954), and can be generalized to the case of many-body Green's function approaches as follows (Hedin and Lundqvist, 1969): Eq. Description Expression This relation was first derived by Baym (1961) by imposing the condition that the nuclei in the crystal must remain near their equilibrium positions due to fictitious restoring forces. Physically this condition corresponds to considering a crystal which is held fixed in the laboratory reference frame. In this approach, the crystal cannot translate or rotate as a whole. Similar relations were derived by Sjölander and Johnson (1965) and Gillis (1970). If we combine Eqs. (119) and (120), perform integrations by parts, and carry out the integrations in r and r we obtain: which fulfils the sum rule in Eq. (120).
Eqs. (118) and (121) completely define the nuclear dynamics in the harmonic approximation. After obtaining the displacement-displacement correlation function D καp,κ α p (tt ) by solving this set of equations, it is possible to construct the expectation value of the nuclear density using Eqs. (107) and (109): (122) We should emphasize that, according to Eq. (121), the coupling of the nuclear displacements to the electrons is completely defined by the electronic dielectric matrix through W e . Similarly, the nuclei affect the electronic structure via the dielectric matrix which enters W ph in Eq. (110) and via the nuclear density inside V tot in Eq. (79). From these considerations it should be clear that the electronic dielectric matrix e (r, r , ω) plays an absolutely central role in the the field-theoretic approach to the electron-phonon problem.

D. Hedin-Baym equations
Apart from making use of the harmonic approximation, the set of equations given by Eqs. (79) (121), and (122) describe the coupled electron-phonon system entirely from first principles. This set of equations can be regarded as the most sophisticated description of interacting electrons and phonons available today. Since the self-consistent equations for the electrons were originally derived by Hedin (1965), and those for the nuclei were derived first by Baym (1961), we will refer to the complete set as the Hedin-Baym equations. Given the importance of these relations, we summarize them schematically in Table I. The standard Hedin's equations for interacting electrons in the potential of clamped nuclei (Hedin, 1965) are immediately recovered from the Hedin-Baym equations by setting to zero the displacement-displacement correlation function of the nuclei, D κpα,κ p α = 0.
Table I provides a closed set of self-consistent equations whose solution yields the Green's functions of a fully-interacting electron-phonon system, within the harmonic approximation. We stress that these relations are fundamentally different from diagrammatic approaches. In fact, here the coupled electron-phonon system is not addressed using Feynman-Dyson perturbation theory as it was done for example by Keating (1968). Instead, in Table I, electrons and phonons are described nonperturbatively by means of a coupled set of nonlinear equations for the exact propagators. In particular we emphasize that this approach does not require the Born-Oppenheimer adiabatic approximation, and therefore it encompasses insulators, intrinsic as well as doped semiconductors, metals, and superconductors.
Almost every property related to electron-phonon interactions in solids that can be calculated today from first principles can be derived from these equations. Examples to be discussed in Secs. V-X include the renormalization of the Fermi velocity, the band gap renormalization in semiconductors and insulators, the non-adiabatic corrections to vibrational frequencies, the Frölich interaction, and the lifetimes of electrons and phonons. The generalization of these results to the case of finite temperature should also be able to describe phonon-mediated superconductivity, although this phenomenon is best addressed by studying directly the propagation of Cooper pairs (see Sec. XI).
Baym's theory can in principle be extended to go beyond the harmonic approximation (Gillis, 1970). However, the mathematical complexity of the resulting formalism is formidable, due to the appearance of many additional terms which are neglected in the harmonic approximation.

V. FROM A MANY-BODY FORMALISM TO PRACTICAL CALCULATIONS
The Hedin-Baym equations summarized in Table I define a rigorous formalism for studying interacting electrons and phonons in metals, semiconductors, and insulators entirely from first principles. However, a direct numerical solution of these equations for real materials is currently out of reach, and approximations are needed for practical calculations. The following sections establish the connection between the Hedin-Baym equations and standard expressions which are currently in use in ab initio calculations of electron-phonon interactions.
A. Effects of the electron-phonon interaction on phonons

Phonons in the Born-Oppenheimer adiabatic approximation
The vibrational eigenmodes of the nuclei can be identified with the resonances of the displacementdisplacement correlation function D κpα,κ p α (tt ) in the frequency domain. If we denote by M the diagonal matrix having the nuclear masses M κ along its diagonal, then the formal solution of Eq. (118) can be written as: where D is the matrix with elements D κpα,κ p α . The resonant frequencies of the system correspond to the solutions of the nonlinear equations: where Ω 2 ν (ω) is an eigenvalue of M −1/2 Π(ω) M −1/2 , parametric in the variable ω.
As expected, the study of lattice vibrations within a field-theoretic framework resembles the standard eigenvalue problem reviewed in Sec. III.A. In particular, the matrix Π(ω) represents the many-body counterpart of the matrix of interatomic force constants C καp,κ α p introduced in Eq. (13). However, despite its formal simplicity, Eq. (124) conceals the full wealth of information associated with the many-body electronic screening e (r, r , ω) via Eq. (121). In fact, the phonon self-energy is generally complex and frequency-dependent. Therefore we can expect to find roots of Eq. (124) outside of the real frequency axis, as well as multiple roots for the same 'eigenmode'.
The link between Eq. (124) and phonon calculations by means of DFT is established by noting that DFT relies on the Born-Oppenheimer adiabatic approximation. In the adiabatic approximation the nuclei are considered immobile during characteristic electronic timescales. Formally, this approximation is introduced by setting ω = 0 in Eq. (121) (Keating, 1968). In practice, this assumption corresponds to stating that e (r, r , ω) can be replaced by e (r, r , 0) in the frequency range of the vibrational excitations. Obviously this is not always the case, and important exceptions will be discussed in Sec. V.A.2.
In order to see more clearly the connection with the formalism discussed in Sec. III.A, we partition the phonon self-energy into 'adiabatic' and 'non-adiabatic' contributions: with Π A = Π(ω = 0). As we will see below, the adiabatic term Π A will be taken to describe 'non-interacting' phonons, and the non-adiabatic self-energy Π NA will be used to describe the effects of electron-phonon interactions.
In the early literature it is common to find a different partitioning, whereby the non-interacting system is defined by the bare interatomic force constants, corresponding to nuclei in the absence of electrons (Grimvall, 1981). This alternative choice is not useful in modern calculations, because the resulting non-interacting phonon dispersions are very different from the fully-interacting dispersions. The present choice of using instead adiabatic phonons as the non-interacting system, is more convenient in the context of modern ab initio techniques, since calculations of adiabatic phonon spectra are routinely performed within DFPT.
In the remainder of this section we concentrate on the adiabatic term, and we defer the discussion of the nonadiabatic self-energy to Sec. V.A.2. Using Eq. (121), we can rewrite the adiabatic self-energy as follows: In this expression U nn is the nucleus-nucleus interaction energy from Eq. (65), V en is the electron-nuclei interaction from Eq. (25), and all the derivatives are taken at the equilibrium coordinates. The derivation of Eq. (126) requires the use of the identity: This identity follows from the same reasoning leading to Eq. (100), after considering an external potential which modifies the position of the nucleus κ in the cell p. Equation (126) can be recast in a familiar form by exploiting the acoustic sum rule in Eq. (120). Indeed after a few tedious but straightforward manipulations we obtain: In this form, one can see that the adiabatic selfenergy gives precisely the interatomic force constants that we would obtain using the Born-Oppenheimer approximation and the Hellman-Feynman theorem, compare Eq. (128) for example with Baroni et al. (2001, p. 517). The difference between the Π A καp,κ α p in Eq. (128) and the C καp,κ α p in Eq. (13) is that, in the former case, the electron density response to atomic displacements is governed by the exact many-body dielectric matrix e (r, r , 0) and electron density n e (r) , as shown by Eqs. (127) and (128). As a result, Π A καp,κ α p corresponds to force constants and electron density dressed by all many-body interactions of the system (both electronelectron and electron-phonon interactions). In contrast, when the force constants in Eq. (13) are calculated using DFT, the electron density response to an atomic displacement is evaluated using the RPA+xc screening, that is Hxc (r, r ) from Sec. III.B.4, and the ground-state electron density is calculated at clamped nuclei.
The use of the adiabatic approximation in the study of phonons carries the important advantage that the manybody force constants Π A form a real and symmetric matrix. This can be seen by rewriting Eq. (121) for ω = 0, and using the relation W e (r, r , ω) = W e (r , r, −ω) which follows from the property W (12) = W (21) (see Sec. IV.B). Since Π A is real and symmetric, all its eigenvalues are guaranteed to be real. In this approximation, the excitations of the lattice correspond to sharp resonances in the displacement-displacement correlation function D(ω), and it is meaningful to talk about phonons as long-lived excitations of the system. In fact these excitations are infinitely long-lived in the harmonic approximation. In practical calculations, the many-body Π A is invariably replaced by the DFT interatomic force constants, and in this case the agreement of the calculated phonon frequencies with experiment is excellent in most cases. Illustrative examples can be found among others in (Yin and Cohen, 1982;Giannozzi et al., 1991;Dal Corso et al., 1993;Lee et al., 1994;Kresse et al., 1995;Dal Corso and de Gironcoli, 2000;Bungaro et al., 2000;Karki et al., 2000;Baroni et al., 2001;Díaz-Sánchez et al., 2007;Dal Corso, 2013).
The most obvious criticism to the adiabatic approximation is that, in the case of metals, the assumption e (r, r , ω) e (r, r , 0) is inadequate. This can intuitively be understood by recalling that the dielectric function of the homogeneous electron gas diverges when ω, q → 0 (Mahan, 1993). In practical calculations, this divergence is connected with vanishing denominators in Eq. (51) for q → 0. An approximate, yet very successful strategy for overcoming this problem, is to replace the occupation numbers in Eq. (51) by smoothing functions such as the Fermi-Dirac distribution, and to describe the singular terms analytically (de Gironcoli, 1995). Most first-principles calculations of phonon dispersion relations in metals have been carried out using this strategy. Improvements upon this strategy will be discussed in Sec. V.A.2.
The adiabatic approximation leads naturally to the definition of an 'adiabatic' propagator D A (ω), which can be obtained from Eq. (123) after replacing the phonon self-energy by its static limit: Now, if we identify Π A καp,κ α p with the interatomic force constant C καp,κ α p in Eq. (13), we can obtain an explicit expression for the adiabatic phonon propagator in terms of the eigenmodes e κα,ν (q) and eigenfrequencies ω qν introduced in Sec. III.A. To this end we invert Eq. (129) using Eqs. (14)- (17), and recall that the dynamical matrix is Hermitian and obeys the relation D dm, * κα,κ α (q) = D dm κα,κ α (−q) (Maradudin and Vosko, 1968). After te-dious but straightforward steps we find: with the definition: This result suggests that, as expected, the propagator should take a simple form in the eigenmodes representation. In fact, by using the inverse transform of Eq. (131) we have: This result can alternatively be obtained starting from the ladder operators of Appendix B. In fact, after using Eqs. (20), (109), and (131) we find: (133) An explicit evaluation of the right-hand side using the Heisenberg time evolution generated by the phonon Hamiltonian in Eq. (22) yields precisely Eq. (132), with the added advantage that it is easier to keep track of the time-ordering. The result is: with η a positive real infinitesimal. This alternative approach is very common in textbooks, see for example Schrieffer (1983) and Schäfer and Wegener (2002). However, it does not carry general validity in a field-theoretic framework since it rests on the adiabatic approximation.

Phonons beyond the adiabatic approximation
In order to go beyond the adiabatic approximation, it is necessary to determine the complete propagator D(ω) in Eq. (123). Formally this can be done by combining Eqs. (123) and (125) to obtain the following Dyson-like equation: In this form it is apparent that the non-adiabatic phonon self-energy Π NA (ω) 'dresses' the non-interacting phonons obtained within the adiabatic approximation, as shown schematically in Fig. 1(a). It is convenient to rewrite the Dyson equation in such a way as to show more clearly the poles of the propagator. Using Eqs. (131)-(132) we find: where Π NA qνν and D −1 qνν are obtained using the transform of Eq. (131) and its inverse, respectively.
From Eq. (136) we see that the non-adiabatic selfenergy Π NA modifies the adiabatic phonon spectrum in four distinct ways: (i) the real part of the diagonal elements Π NA qνν shifts the adiabatic frequencies; (ii) the imaginary part introduces spectral broadening; (iii) the offdiagonal elements of Π NA qνν introduce a coupling between the adiabatic vibrational eigenmodes; (iv) the frequencydependence of Π NA qνν (ω) might lead to multiple poles for the same mode ν, thereby introducing new structures in the phonon spectrum.
Today it is relatively common to calculate phonon linewidths arising from electron-phonon interactions (Allen, 1972b;Bauer et al., 1998). Recently it has also become possible to study the frequency renormalization due to non-adiabatic effects (Calandra et al., 2010;Saitta et al., 2008).
The possibility of observing new features in vibrational spectra arising from the EPI has not been studied from first principles, but the underlying phenomenology should be similar to that of plasmon satellites in photoelectron spectra (see Sec. V.B.6). Generally speaking we expect satellites whenever e (r, r , ω) exhibits dynamical structure close to vibrational frequencies. This can happen, for example, in the case of degenerate polar semicondutors, when phonon and plasmon energies are in resonance. In these cases, phonons and plasmons can combine into 'coupled plasmon-phonon modes' (Richter, 1984), which are the electronic analogue of photon polaritons (Yu and Cardona, 2010, p. 295). This phenomenon was predicted theoretically (Varga, 1965), and subsequently confirmed by Raman measurements on GaAs (Mooradian and Wright, 1966). We speculate that it should be possible to obtain coupled plasmon-phonon modes from the frequency-dependence of the phonon selfenergy in Eq. (136); it would be interesting to perform first-principles calculations in order to shed light on these aspects.
In practical calculations the non-adiabatic corrections to the adiabatic phonon spectrum are evaluated from Eq. (136) using first-order perturbation theory, by retaining only the diagonal elements of Π NA . If we denote the complex zeros of D −1 qνν (ω) byΩ qν = Ω qν − iγ qν , in the case of non-degenerate eigenmodes Eq. (136) gives: therefore: Apart from the small γ 2 qν term in the second line, these expressions are identical to those provided by Allen (1972b) and Grimvall (1981). Since non-adiabatic corrections are usually small as compared to the adiabatic phonon frequencies, the above expressions are often simplified further by using the additional approximations |Ω qν − ω qν | ω qν and |γ qν | ω qν , leading to γ qν −Im Π NA qνν (ω qν ) and Ω qν ω qν + Re Π NA qνν (ω qν ). In these forms it becomes evident that the real part of the selfenergy shifts the adiabatic phonon frequencies, and the imaginary part is responsible for the spectral broadening of the resonances. Using these expressions in Eq. (136) and going back to the time domain, it is seen that, as a result of the EPI, phonons acquire a finite lifetime given by τ ph qν = (2γ qν ) −1 .

Expressions for the phonon self-energy used in ab initio calculations
In the literature on electron-phonon interactions, the phonon self-energy Π is almost invariably expressed in terms of an electron-phonon vertex g and the electron Green's function G as Π = |g| 2 GG in symbolic notation, see for example Grimvall (1981, p. 195). While this has become common practice also in ab initio calculations, the origin of this choice is not entirely transparent. One could derive the above expression directly from Eq. (1), using standard Green's function techniques. However, this procedure does not answer the key question on how to calculate the electron-phonon matrix elements g.
Closer inspection of the theory reveals that this is a rather nontrivial point. In fact, on the one hand, a straightforward expansion of the second-quantized Hamiltonian of Eq. (68) in terms of the nuclear coordinates leads to 'bare' electron-phonon matrix elements, g b , which contain the bare Coulomb interaction between electrons and nuclei. On the other hand, if we go back to Sec. III.B.2, we see that the electron-phonon matrix elements in DFT are 'dressed' by the self-consistent response of the electrons. The difference between bare and dressed vertex is not only quantitative, but also qualitative: for example in metals the bare vertex is long-ranged, while the screened vertex is short-ranged.
The relation between bare and dressed electron-phonon vertices and the derivation of explicit expressions for the phonon self-energy have been discussed by many authors, see for example Rickayzen (1980) and Scalapino (1969). In short the argument is that the lowest-order Feynman diagram starting and ending with a phonon line must contain precisely two bare electron-phonon vertices, as shown in Fig. 1(b). By construction this diagram corresponds to having Π = |g b | 2 GG. In order to make the transition from the bare vertex to the dressed vertex it is necessary to collect together all the proper electronic polarization diagrams around the vertex. However, these steps have been carried out only for the homogeneous electron gas (Rickayzen, 1980;Scalapino, 1969). In the following we show how the dressed electron-phonon vertex emerges from a non-perturbative analysis based on the Hedin-Baym equations. The non-adiabatic phonon self-energy Π NA introduced in Sec. V.A.2 can be written explicitly by combining Eqs. (119) and (125): (140) Using the Dyson equation for the screened Coulomb interaction, it can be seen that this expression does indeed contain electron-phonon matrix elements. In fact, by inserting Eq. (94) into Eq. (140) we find terms like vP e W e , and the electron-phonon matrix elements will arise from taking the gradients of v and W e , respectively. By working in the eigenmodes representation via Eq. (131), after lengthy manipulations this procedure yields: where we introduced electron-phonon 'coupling functions' as follows. The bare coupling g b is defined as: where V en is the potential of the nuclei from Eq. (25); in practical calculations this quantity is replaced by the usual ionic pseudo-potentials. The meaning of the variation ∆ qν is the same as in Eqs. (33)-(35). The dressed couplings g and g cc are defined as (Hedin and Lundqvist, 1969): Since the dielectric matrix is real at ω = 0, we have the simple relation g cc qν (r, 0) = g * qν (r, 0). In order to derive Eq. (141) it is best to carry out the algebra in the time domain. We emphasize that the result expressed by Eq. (141) is non-perturbative, and relies solely on the harmonic approximation.
Equation (141) is in agreement with the standard result for the homogeneous electron gas (Scalapino, 1969). The same expression was also obtained by Keating (1968) using a detailed diagrammatic analysis. Keating's diagrammatic representation of the self-energy is shown in Fig. 1(c), and can be obtained from Eq. (141) by noting that, in symbolic notation, P e = GGΓ from Eq. (92), therefore Π = g b GGΓg. For completeness we also show in Fig. 1(d) a diagrammatic representation of the dressed electron-phonon coupling function g as given by Eq. (144). This representation is obtained by observing that g = −1 g b , = 1 − vP , and P = GGΓ, therefore g = g b + vGGΓg.
In view of practical first-principles calculations it is useful to have a simplified expression for the nonadiabatic phonon self-energy in Eq. (141). To this aim we make the following approximations: (i) The vertex function in Eq. (92) is set to Γ(123) = δ(12)δ(13). This is the same approximation at the heart of the GW method (Hedin, 1965;Hybertsen and Louie, 1986;Onida et al., 2002); (ii) The fully-interacting electron Green's function G is replaced by its non-interacting counterpart, using the Kohn-Sham eigenstates/eigenvalues evaluated with the nuclei held in their equilibrium positions; (iii) The fully-interacting dielectric matrix in Eq. (144) is replaced by the RPA+xc response obtained from a DFT calculation, as discussed in Sec. III.B.2; (iv) The frequency-dependence of the screened electronphonon coupling defined in Eq. (144) is neglected: g cc qν (r, ω) g cc qν (r, 0) = g * qν (r, 0). This approximation is ubiquitous in the literature but it is never mentioned explicitly; (v) For notational simplicity, we consider a spindegenerate system with time-reversal symmetry; this simplification is easily removed.
Using these assumptions we can rewrite the component of Eq. (141) for q = q as: We note that the components of the phonon self-energy for q = q vanish due to the periodicity of the crystalline lattice. In Eq. (145) the sums run over all the Kohn-Sham states, with occupations f nk , and η is a real positive infinitesimal. In this case we indicate explicitly the factor of 2 arising from the spin degeneracy. The matrix element g mnν (k, q) is the same as in Eq. (38), and it is precisely the quantity calculated by most linearresponse codes. The matrix element g b mnν (k, q) is obtained from g mnν (k, q) by replacing the variation of the Kohn-Sham potential by the corresponding variation of the ionic (pseudo)potentials. The field-theoretic phonon self-energy given by Eq. (145) is in agreement with the expression derived by Calandra et al. (2010) starting from time-dependent density-functional perturbation theory.
The presence of both the bare electron-phonon matrix element and the screened matrix element in Eq. (145) has not been fully appreciated in the literature, and most ab initio calculations employ an approximate self-energy whereby g b is replaced by g. The replacement of the bare matrix elements by their screened counterparts in the phonon self-energy goes a long way back, and can be found already in the seminal work by Allen (1972b). As a result many investigators (including the author) calculated phonon lifetimes using the following expression (Grimvall, 1981): This is obtained from Eq. (145) by taking the imaginary part and by making the replacement g b → g. While Eq. (146) can be derived from the Fermi golden rule in a independent-particle approximation (see Albers et al., 1976, Appendix B), the choice of the electron-phonon matrix elements is somewhat arbitrary. In future calculations of the phonon self-energy it will be important to assess the importance of using the correct vertex structure, that is replacing |g mnν (k, q)| 2 by g b mnν (k, q)g * mnν (k, q) in Eq. (146).
In general, the effects of the non-adiabatic self-energy on the phonon spectrum are expected to be significant only in the case of metals and small-gap semiconductors. In fact, by combining Eqs. (138), (139), and (145) it is seen that Π NA can be large only when occupied and empty single-particle states are separated by an energy of the order of the characteristic phonon energy. In such a case, we can expect a shift of the adiabatic phonon frequencies, and a concomitant broadening of the lines. A clear illustration of these effects was provided by Maksimov and Shulga (1996), who analyzed a simplified model of a metal with linear bands near the Fermi level.
Calculations of phonon linewidths based on Eq. (146) have been reported by several authors 2 and have become commonplace in first-principles studies of electronphonon physics. On the other hand, calculations of the non-adiabatic phonon frequencies using Eq. (145) have only been reported in Piscanec et al., 2007;Caudal et al., 2007;Saitta et al., 2008;Calandra et al., 2010), using the approximation that the bare vertex g b can be replaced by the screened vertex g. Examples of such calculations will be reviewed in Sec. VII.
Equation (145) suggests several avenues worth exploring in the future: firstly, the use of the bare vertex should not pose a challenge in practical calculations, since this quantity is already being calculated in linear-response DFT codes. Testing the impact of the bare vertex on phonon linewidths and frequency renormalizations will be important. Secondly, Eq. (145) contains off-diagonal couplings, which are usually ignored. It will be interesting to check the effect of using the complete matrix self-energy. Thirdly, the dynamical structure of the selfenergy may contain interesting information, such as for example spectral satellites and coupled phonon-plasmon modes. Lastly, the move from Eq. (141) to Eq. (145) involves the approximation that the frequency-dependence of the electron-phonon matrix elements can be neglected. The validity of this approximation is uncertain, and there are no reference ab initio calculations on this. However, we note that frequency-dependent electron-phonon matrix elements have been employed systematically in theoretical models of doped semiconductors (Mahan, 1993, Sec. 6.3).
Before closing this section we note that the formalism discussed here is based on zero-temperature Green's functions. In order to extend the present results to finite temperature it is necessary to repeat all derivations using the Matsubara representation, and then perform the analytic continuation of the self-energy to the real frequency axis. Detailed derivations can be found in (Baym, 1961) and (Gillis, 1970), and will not be repeated here. Fortunately it turns out that Eq. (145) Table I and ab initio calculations of phonons. We first identified a Hermitian eigenvalue problem for the vibrational frequencies via the adiabatic approximation, and then we improved upon this description by means of a non-adiabatic selfenergy. In this section, we adopt a similar strategy in order to discuss electronic excitation energies: first we identify an approximation to the Hedin-Baym equations which does not include any electron-phonon interactions, and then we introduce an electron self-energy to incorporate such interactions.
The single most common approximation in firstprinciples electronic structure calculations is to describe nuclei as classical particles clamped in their equilibrium positions. Within this approximation the expectation value of the nuclear charge density operator in Eq. (67), n n (r) , is replaced by the first term in Eq. (107), n 0 n (r). From Eq. (122) we see that this approximation formally corresponds to setting to zero the displacementdisplacement correlation function of the nuclei. This observation suggests that, in order to unambiguously single out electron-phonon interactions in the Hedin-Baym equations, we need to define a non-interacting problem by setting D καp,κ α p = 0, and identify the electronphonon interaction with the remainder. In the following, we write an equation of motion for the electrons analogous to Eq. (84), except with the nuclei clamped in their equilibrium positions; then we use a Dyson-like equation to recover the fully-interacting electron Green's function.
Here the potential V cn tot differs from its counterpart V tot of Eq. (79) in that the total density of electrons and nuclei, n , is replaced by the density calculated at clamped nuclei, n cn : where we defined: The term Σ cn e in Eq. (147)  Equations (147)-(150) lead directly to the well-known Hedin's equations (Hedin, 1965). Hedin's equations and the associated GW method at clamped nuclei are addressed in a number of excellent reviews (Hedin and Lundqvist, 1969;Hybertsen and Louie, 1986;Aryasetiawan and Gunnarsson, 1998;Onida et al., 2002) hence they will not be discussed here.
We emphasize that Eqs. (147)-(155) are just an alternative formulation of the Hedin-Baym equations in Table I. The advantage of this formulation is that it better reflects standard practice, whereby the DFT equations and the GW quasiparticle corrections are evaluated at clamped nuclei. Equations (147)-(155) are formally exact within the harmonic approximation. A schematic representation of the Dyson equation for the electron Green's function and the decomposition of the electron self-energy are given in Fig. 2. The selfenergy contribution Σ FM in Eq. (153) is a dynamic correction to the electronic excitation energies, and is analogous to the GW self-energy at clamped nuclei. Indeed, in the same way as the correlation part of the standard GW self-energy describes the effect of the dynamic electronic polarization upon the addition of electrons or holes to the system, the term GW ph in Eq. (153) describes the effect of the dynamic polarization of the lattice.
The static term Σ DW in Eq. (154) corresponds to the difference between the self-consistent potential V tot calculated for the fully-interacting system, and the same potential evaluated with the nuclei clamped in their equilibrium positions, V cn tot . Intuitively this term corresponds to a time-independent correction to the 'crystal potential' that arises from the fuzziness of the nuclear charge density around the equilibrium nuclear positions. This term is similar to the one appearing in the study of the temperature dependence of X-ray diffraction and neutron diffraction spectra (Mermin, 1966;Ashcroft and Mermin, 1966), and is commonly referred to as the Debye-Waller (DW) term (Antoncík, 1955;Walter et al., 1970;Cardona and Thewalt, 2005). Hedin and Lundqvist (1969) did not include this term in their classic work, however this contribution was discussed by Allen and Heine (1976; see also Allen, 1978).
The last term Eq. (152), Σ dGW , is the correction to the standard Hedin self-energy arising from the fact that the fully-interacting electron Green's function and density are slightly different from those evaluated at clamped nuclei, owing to the electron-phonon interaction. The magnitude of this term corresponds to a fraction of the GW Γ quasiparticle corrections at clamped nuclei. Since Σ dGW has never been investigated so far, we will not discuss this term further.

Expressions for the electron self-energy used in ab initio calculations
A complete self-consistent solution of Eqs. (147)-(155) from first principles is not possible at present, and one has to replace the various entries of Eq. (152) by the best approximations available. In practice, one resorts to either DFT or to GW calculations; recent progress will be reviewed in Secs. VIII and IX.
Using Eqs. (110), (130)-(131), and (144) we can rewrite the Fan-Migdal self-energy as follows: ω). (156) This shows that the Fan-Migdal self-energy is, in symbolic notation, of the type Σ = g 2 DGΓ; a graphical representation of this term is given in Fig. 2(c). In order to make the last expression amenable to ab initio calculations, it is common to make the following approximations, which are similar to those introduced earlier for the phonon self-energy: (i) The vertex Γ(123) is set to δ(13)δ(23); (ii) The fully-interacting electron Green's function is replaced by the Kohn-Sham Green's function evaluated at clamped nuclei; (iii) The fully-interacting phonon propagator D qνν (ω) is replaced by the adiabatic propagator D A qνν (ω) from Eq. (134); (iv) The screened electron-phonon vertex is evaluated using the RPA+xc electronic screening from a DFT calculation; (v) The frequency dependence of the electron-phonon coupling is neglected, g qν (r, ω) g qν (r, 0).
After using these approximations in Eq. (156), we obtain the following result for the k = k matrix elements of the Fan-Migdal self-energy in the basis of Kohn-Sham states: Here f mk+q = 1 for occupied Kohn-Sham states and 0 otherwise, and the matrix element g mnν (k, q) is obtained from Eq. (38). As for the phonon self-energy in Eq. (145), also in this case the components of the electron self-energy for k = k vanish due to the periodicity of the lattice. The result in Eq. (157) is obtained by closing the contour of the frequency integration in the upper complex plane, owing to the t + 1 in the exponential of Eq. (156). The infinitesimals inside the electron and phonon propagators, which reflect the time-ordering, are crucial to obtain the correct result (Schrieffer, 1983). The spin label is omitted in Eq. (157) since this contribution to the self-energy is diagonal in the spin indices.
The result above is only valid at zero temperature. The extension to finite temperature requires going through the Matsubara representation, and then continuing the self-energy from the imaginary axis to the real axis. The procedure is described in many textbooks, see for example Sec. 3.5 of (Mahan, 1993). The result is that at finite temperature the square brackets of Eq. (157) are to be modified as follows: where f mk+q and n qν are now Fermi-Dirac and Bose-Einstein distribution functions, respectively (see footnote 3). The change of sign in the imaginary infinitesimal on the second fraction has to do with the fact that in the Matsubara formalism the analytic continuation from the imaginary frequency axis to the real axis through the upper complex plane leads to the so-called retarded selfenergy, that is a self-energy with all poles below the real axis (Abrikosov et al., 1975;Mahan, 1993).
The Debye-Waller term in Eq. (154) can also be written in a form which is convenient for practical calculations, by expanding the total density operatorn(3) to second order in the atomic displacements. Using Eqs. (79) and (109) we find: (159) In order to arrive at this result, it is necessary to make the additional approximation that the electronic field operators and the operators for the nuclear displacements are uncorrelated, that is n ∆τ καp = n ∆τ καp , and similarly for the second power of the displacements. This requirement was noticed by Gillis (1970), and is trivially satisfied if we describe phonons within the adiabatic approximation of Sec. V.A.1. Equation (159) motivates the diagrammatic representation of the Debye-Waller selfenergy shown in Fig. 2(d), whereby the phonon line begins and ends at the same time point. We note that Eq. (159) involves the variation of the screened potential V tot ; this result, which we here derived starting from Schwinger's functional derivative technique, is also obtained when starting from a perturbative diagrammatic analysis . The Debye-Waller selfenergy can be simplified further if we make use of the following approximations, in the same spirit as for the Fan-Migdal self-energy: (vi) The fully-interacting phonon propagator is replaced by the adiabatic propagator D A qνν (ω) from Eq. (134); (vii) The total many-body potential V tot of Eq. (79) is replaced by the Kohn-Sham potential V KS (r) evaluated at clamped nuclei. Strictly speaking, the Kohn-Sham effective potential includes also contributions from exchange and correlation, which in the Hedin-Baym equations are all contained in the electron self-energy. However, the present discussion holds unchanged if we add any local and frequencyindependent potential to V tot in Eq. (84), while removing the same potential from the self-energy (Keating, 1968).
Using these simplifications together with Eqs. (130) and (131), we can write the Σ DW in the basis of Kohn-Sham eigenstates as follows: where the Debye-Waller matrix element g DW is obtained from Eq. (40), and the presence of only the diagonal terms ν = ν is a result of the Kronecker delta in Eq. (132).
In going from Eq. (159) to Eq. (160) the frequency integration is performed by using Eq. (134), after closing the contour in the lower half plane. The resulting expression is diagonal in the spin indices.
The expression for the Debye-Waller term in Eq. (160) is only valid at zero temperature. In this case the extension to finite temperature is immediate since the selfenergy does not involve the electron Green's function, hence we only need to evaluate the canonical average of Eq. (133) at equal times. The result is that Eq. (160) is simply to be modified as follows: with n qν being the Bose-Einstein occupations (see footnote 3). The Debye-Waller contribution to the electron selfenergy is almost invariably ignored in the literature on metals and superconductors, but it is well-known in the theory of temperature-dependent band structures of semiconductors (Allen and Heine, 1976;Allen and Cardona, 1981;Marini, 2008;Giustino et al., 2010;Poncé et al., 2015). Neglecting Σ DW in metals is partly justified by the fact that this term is frequency-independent, therefore it is expected to be a slowly-varying function over each Fermi surface sheet. A detailed first-principles analysis of this aspect is currently lacking.

Temperature-dependence of electronic band structures
Once determined the electron self-energy as in Sec. V.B.2 it is possible to study the modification of the electronic structure induced by the EPI. To this aim it is convenient to rewrite Eq. (151) in the basis of Kohn-Sham eigenstates: Assuming that the electronic structure problem at clamped nuclei has been solved using DFT or DFT+GW calculations, the Green's function G cn can be written in terms of simple poles at the Kohn-Sham or quasiparticle eigenvalues ε nk (Hedin and Lundqvist, 1969). In this case Eq. (162) reduces to: whereε nk = ε nk ± i η with the upper/lower sign corresponding to occupied/unoccupied states. The spin indices are omitted since these self-energy contributions do not mix states with opposite spin. In order to gain insight into the effects of the electronphonon interaction, we start from the drastic approximation that Σ ep only leads to a small shift of the quasiparticle poles, from the 'non-interacting' energies ε nk to the renormalized energiesẼ nk = E nk + iΓ nk . In this approximation, the fully interacting Green's function is expressed as a sum of simple poles, given by the zeros of Eq. (163): As in the case of vibrational frequencies in Eq. (137), we are considering for simplicity non-degenerate electronic states, and making the assumption that the off-diagonal elements of the self-energy Σ ep nn k with n = n can be neglected. In more general situations the right-hand side of Eq. (163) needs to be to diagonalized, or alternatively the off-diagonal terms Σ ep nn k need to be treated perturbatively. The energies E nk obtained from Eq. (164) yield the band structure renormalized by the EPIs, to be discussed below. The imaginary part Γ nk in Eq. (165) is connected with the quasiparticle lifetimes and will be discussed in Sec. V.B.4.
Equation (164) is to be solved self-consistently for E nk and Γ nk . When Eq. (164) is used in combination with the standard approximations to the Fan-Migdal and Debye-Waller self-energies given by Eqs. (157) and (160), the result that one obtains is equivalent to describing electron-phonon couplings to second order in Brillouin-Wigner perturbation theory (Mahan, 1993). Similarly one recovers the more basic Rayleigh-Schrödinger perturbation theory by making the replacements E nk → ε nk and Γ nk → 0 in Eq. (164).
By combining Eqs. (157)-(158), (160)-(161), and (164), we obtain the temperature-dependent 'band structure renormalization' arising from the EPI: For practical calculations it is important to bear in mind that this result rests on the approximations (i)-(vii) introduced at p. 29, as well as the harmonic approximation. The theory of temperature-dependent band structures developed by Allen and Heine (1976) makes two additional approximations on top of Eq. (166): Brillouin-Wigner perturbation theory is replaced by Rayleigh-Schrödinger perturbation theory; and the phonon energies in the denominators are neglected. Using these additional approximations Eq. (166) becomes: which is referred to as the 'adiabatic Allen-Heine formula'. By setting T = 0 the Bose-Einstein factors n qν vanish and we have the so-called 'zero-point renormalization' of the energy bands, ∆E ZP nk = E nk (T = 0) − ε nk . This is the modification of the electronic energies evaluated at clamped nuclei, which arises from the zero-point fluctuations of the atoms around their equilibrium sites.
An expression that is essentially identical to Eq. (167) can also be obtained directly from Eq. (1) using secondorder Raleigh-Schrödinger perturbation theory in Fock space, following the same lines as in (Kittel, 1963, p. 134). A detailed derivation of the formalism starting from Eq. (1) was given by Chakraborty and Allen (1978).
Historically, the Allen-Heine theory (Allen and Heine, 1976) was developed by starting from a straightforward perturbation theory expansion of the electron energies in terms of the atomic displacements within the adiabatic approximation, followed by a canonical average of the displacements using Bose-Einstein statistics. It is reassuring that, after making a few well-defined approximations, a field-theoretic method leads to the same result.
Equation (167) was employed in many semi-empirical calculations. 4 More recently, this expression was used in the context of first-principles DFT calculations by Marini (2008) and Giustino et al. (2010). DFT calculations of band structure renormalization based on Eqs. (166) or (167) are becoming increasingly popular, and the latest developments will be reviewed in Sec. IX.A.1.
The nature of the band structure renormalization by electron-phonon interactions can be understood at a qualitative level by considering a drastically simplified model, consisting of a semiconductor with parabolic and nondegenerate valence and conduction bands, with the band extrema coupled to all other states by a dispersionless phonon mode of frequency ω 0 . If the Debye-Waller matrix elements are much smaller than the Fan-Migdal matrix elements, then the dominant contributions to Eq. (167) arise from denominators such as ε nk − ε nk+q ± 2 |q| 2 /2m * n , where the upper/lower sign is for the valence/conduction band, and m * n are the effective masses. As a result the temperature dependence of the band gap will take the form: where E cn g is the gap at clamped nuclei, and |∆E ZP g | is the zero-point correction. The negative sign in the last expression arises from the opposite curvatures of the valence and conduction bands. In this example the band gap decreases with temperature: this is a well known effect in semiconductor physics, and is often referred to as the 'Varshni effect' (Varshni, 1967). The first measurements of such effects were performed by Becker and Fan (1949), and stimulated the development of the first theory of temperature-dependent band gaps (Fan, 1951). A schematic illustration of this qualitative model is provided in Fig. 3(a). The redshift of the gap as a function of temperature is seen in many albeit not all semiconductors. For example, copper halides (Göbel et al., 1998) and lead halide perovskites (D'Innocenzo et al., 2014) exhibit an 'inverse Varshni' effect, that is a blueshift of the gap with temperature; in addition some chalcopyrites exhibit a non-monotonic temperature dependence of the band gap (Bhosale et al., 2012). We also point out that the qualitative model shown in Fig. 3(a) does not take into account the subtle temperature-dependence of the band gap renormalization at very low temperature. These effects were recently investigated by Allen and Nery (2017).

Carrier lifetimes
While the real part of the poles in Eq. (164) describes the energy level renormalization induced by the electronphonon coupling, the imaginary part Γ nk in Eq. (165) carries information on the spectral broadening, which will be discussed in Sec. (V.B.5), and on quasiparticle lifetimes, which we discuss below.
After transforming G nn k (ω) from Eq. (163) into the time domain it is seen that, for an electron or hole added to the system at time t in the state |nk , the probability amplitude to persist in the same state decreases as exp[Γ nk (t − t)/ ]. Using Eqs. (157) and (165) it can be seen that Γ nk < 0 for an electron added to the system and Γ nk > 0 for a hole. Therefore the average time spent by the particle in the state |nk is τ nk = /(2|Γ nk |).
A popular expression for the electron and hole lifetimes is obtained by making the replacementẼ nk → ε nk in Eq. (157), and by taking the absolute value of the imaginary part. We find: A more accurate expression is discussed after Eq. (174) in the next section. The extension of the above result to finite temperature is obtained by taking the absolute value of the imaginary part of Eq. (158): We emphasize the change of sign in the third line, resulting from the analytic continuation to the retarded self-energy. Equation (170) coincides with the expression that one would obtain by using the standard Fermi golden rule (Grimvall, 1981). The intuitive interpretation of this result is that the quasiparticle lifetime is reduced by processes of phonon emission and absorption, corresponding to the second and third lines of Eq. (170), respectively. We note that in deriving Eq. (170) we did not consider the Debye-Waller self-energy; this is because the diagonal matrix elements of Σ DW are purely real, hence they do not contribute to the quasiparticle widths (Lautenschlager et al., 1986). Ab initio calculations of carrier lifetimes using Eq. (170) were first reported by Eiguren et al. (2003aEiguren et al. ( , 2002. These applications and more recent developments will be reviewed in Sec. X.A. If we evaluate Eq. (170) for the same simplified model introduced for the temperature renormalization, and we neglect the phonon energy in the Dirac delta functions, we obtain Γ nk (T ) = Γ ZP nk [1 + 2 n( ω 0 /k B T )] where Γ ZP nk is the linewidth at T = 0. The dependence of the linewidth and the corresponding lifetime on temperature for this model are shown in Fig. 3(b). This trend is typical in semiconductors (Lautenschlager et al., 1987a,b).

Kinks and satellites
In many cases of interest, the use of Brillouin-Wigner perturbation theory as given by Eqs. (164)-(165) is not sufficient to provide an adequate description of EPIs, and it becomes necessary to go back to the complete Dyson equation, Eq. (163). Generally speaking a direct solution of the Dyson equation is important in all those cases where the electronic energy scales are comparable to phonon energies, namely in metals (including superconductors), narrow-gap semiconductors, and doped semiconductors. In order to study these systems, it is convenient to introduce an auxiliary function called the 'spectral density function', or simply spectral function.
In its simplest version the spectral function is defined as (Abrikosov et al., 1975;Mahan, 1993): where the superscript 'ret' stands for 'retarded', and simply indicates that all poles of the Green's function G nn k (ω) in the upper complex plane must be replaced by their complex conjugate. The spectral function is positive definite and carries the meaning of a 'many-body momentum-resolved density of states' (Abrikosov et al., 1975). This is precisely the function that is probed by angle-resolved photoelectron spectroscopy experiments or ARPES (Damascelli et al., 2003). Using Eq. (163) the spectral function can be rewritten as: (172) In order to obtain the correct spectral function, it is important to use the retarded self-energy. This is done by using Eq. (158) for the Fan-Migdal term, while the static Debye-Waller term remains unchanged.
It is often convenient to approximate the spectral function as a sum of quasiparticle peaks. To this aim, one performs a linear expansion of Eq. (172) around each quasiparticle energy E nk , to obtain: (173) This is a sum of Lorentzians with strength Z nk and width Z nk Im Σ ep nnk (E nk / ). Here the 'quasiparticle strength' is defined as the homonymous quantity appearing in GW calculations (Hedin and Lundqvist, 1969): The result expressed by Eq. (173) shows that, in a rigorous field-theoretic approach, the quasiparticle broadening and lifetime given by Eqs. (169) and (170) should be renormalized by Z nk and Z −1 nk , respectively, and should be evaluated using the quasiparticle energy E nk instead of ε nk . This result can also be derived from Eq. (165) by performing a Taylor expansion of the self-energy along the imaginary axis and using the Cauchy-Riemann conditions.
In order to illustrate the typical features of the spectral function, we consider a model system characterized by one parabolic conduction band. The occupied electronic states couple to all states within an energy cutoff via a dispersionless phonon mode and a constant electronphonon matrix element. A simplified version of this model was discussed by Engelsberg and Schrieffer (1963) by considering a constant density of electronic states. By evaluating the spectral function in Eq. (172) using the Fan-Migdal self-energy and neglecting the Debye-Waller term, we obtain the results shown in Fig. 4 for two sets of parameters.
In Fig. 4(a) the Fermi energy is much larger than the characteristic phonon energy. This case is representative of a metallic system with electron bands nearly linear around the Fermi level. Here the electron-phonon interaction leads to (i) a reduction of the band velocity in proximity of the Fermi level, and (ii) a broadening of the spectral function beyond the phonon energy ω 0 . A detailed analysis of these features for a slightly simpler model system, including a discussion of the analytic properties of the Green's function, can be found in the work by Engelsberg and Schrieffer.
The solid line (red/dark gray) in Fig. 4(a) shows the renormalized band structure obtained from Brillouin-Wigner perturbation theory, Eq. (164). We see that these solutions track the maxima of the spectral function A(k, ω). The renormalized bands exhibit a characteristic 'S-shape' near the Fermi level, corresponding to multiple solutions of Eq. (164) for the same wavevector k. Starting from the late 1990s such S-shaped energymomentum dispersion curves have been observed in a number of ARPES experiments, and have become known in the literature as the 'photoemission kink' (Valla et al.,FIG. 4 (Color online) Two-dimensional maps of the electron spectral function A(k, ω) for electrons coupled to a dispersionless phonon of frequency ω0. The non-interacting bands are given by ε(k) = −εF + 2 |k| 2 /2m * , and the Fermi level coincides with the top of the energy window. The matrix element is |g| 2 = ω0/NF when the electron energies differ by less than the cutoff εmax, and zero otherwise (NF is the density of states at the Fermi level). (a) Spectral function for the case εF = 10 ω0 (white on blue/black), non-interacting band structure (solid line, yellow/light gray), and fully-interacting band structure within Brillouin-Wigner perturbation theory (solid line, red/dark gray). (b) Spectral function for the case εF = 2 ω0. The model parameters are: m * = 0.1 me, ω0 = 100 meV, η = 20 meV, εc = 5 eV, For clarity the calculated spectral functions are cut off at the value 3 eV −1 and normalized. The self-energy is shifted by a constant so as to have Σ ep (0) = 0; this correction guarantees the fulfillment of Luttinger's theorem about the volume enclosed by the Fermi surface (Luttinger, 1960). 1999a). First-principles calculations of kinks were first reported by Eiguren et al. (2003a) and Giustino et al. (2008), and will be reviewed in Sec. VIII.
In Fig. 4(b) the Fermi energy is comparable to the characteristic phonon energy. This case is representative of a degenerately doped semiconductor close to a conduction band minimum. Here the electron-phonon interaction leads to two distinct spectral features: (i) a parabolic band with a heavier mass, which is well described by the Brillouin-Wigner solutions (solid line, red/dark gray), and (ii) a polaron satellite that is visible further down. In this example, it is clear that Eq. (164) is unable to describe the satellite, and that the spectral function carries qualitatively new information about the system. Polaron satellites resembling Fig. 4(b) have been observed in ARPES experiments on doped oxides (Moser et al., 2013;Chen et al., 2015;Cancellieri et al., 2016;Wang et al., 2016) and recently calculated from first principles .

Model Hamiltonians, polarons, and the cumulant expansion
At the end of this section it is worth mentioning complementary non first-principles approaches for studying the effects of EPIs on the electronic properties of solids. Model EPI Hamiltonians can be derived from Eq. (1) by choosing a priori explicit expressions for the electron band energies, the vibrational frequencies, and the coupling matrix elements. Examples of model Hamiltonians are those of Fröhlich (1954), Holstein (1959), Su, Schrieffer, and Heeger (1979), the Hubbard-Holstein model (Berger et al., 1995), the Peierls-Hubbard model (Campbell et al., 1984), the 't-J' Holstein model (Rösch and Gunnarsson, 2004), and the Su-Schrieffer-Heeger-Holstein model (Perroni et al., 2004). These models involve the tight-binding approximation, the Einstein phonon spectrum, and electron-phonon couplings to first order in the atomic displacements. Using these model Hamiltonians it is possible to go beyond the approximations introduced in Sec. V.B.2, and obtain nonperturbative solutions by means of canonical Lang-Firsov transformations, path-integral methods, exact diagonalization, variational or quantum Monte Carlo techniques (Alexandrov, 2008;Alexandrov and Devreese, 2010). These models have been used extensively to explore many aspects of polaron physics, for example the ground-state energy of polarons (weak or strong coupling), their spatial extent (large or small polarons), and transport properties (band-like or hopping-like).
Given the considerable body of literature on model EPI Hamiltonians, it is natural to ask whether one could bring ab initio calculations of EPIs to a similar level of sophistication. The main limitation of current first-principles approaches is that, given the complexity of the calculations, the electron self-energies are evaluated using the bare propagators, as in Eq. (157). As a consequence, higher-order interaction diagrams beyond the Migdal approximation (Migdal, 1958) are omitted altogether.
A promising avenue for going beyond the Migdal approximation consists of introducing higher-order diagrams via the 'cumulant expansion' approach (Hedin and Lundqvist, 1969;Langreth, 1970;Aryasetiawan et al., 1996). In the cumulant expansion method, instead of calculating the electron Green's function via a Dyson equation, one evaluates the time evolution of the Green's function by formulating the problem in the interaction picture, in symbols: (Aryasetiawan et al., 1996). The distinctive advantage of this approach is that the 'cumulant' C nnk (t) can be obtained from a low-order self-energy, for example the Fan-Migdal self-energy in Eq. (157), and the exponential 'resummation' automatically generates higherorder diagrams (Mahan, 1993, pag. 523). Detailed discussions of the cumulant expansion formalism are given by Zhou et al. (2015) and Gumhalter et al. (2016).
The cumulant method provides an interesting point of contact between ab initio and model Hamiltonian approaches. In fact, the cumulant expansion is closely related to the 'momentum average approximation' introduced by Berciu (2006) for studying the Green's function of the Holstein polaron.
The cumulant expansion has proven successful in ab initio calculations of electron-electron interactions, in particular valence band satellites in semiconductors (Kheifets et al., 2003;Guzzo et al., 2011Guzzo et al., , 2012Guzzo et al., , 2014Lischner et al., 2013;Kas et al., 2014;Caruso et al., , 2015b. In the context of EPIs, the ab initio cumulant expansion method has been applied to elemental metals by Story et al. (2014), and to the ARPES spectra of ndoped TiO 2 by Verdi et al. (2016). In the latter work the cumulant method correctly reproduced the polaron satellites observed in the experiments of Moser et al. (2013).
The study of polarons using ab initio many-body techniques is yet to begin, however a first calculation of the spectral function of Fröhlich polarons and an approximate polaron wavefunction have recently been reported by Verdi et al. (2016).

VI. EFFICIENT CALCULATIONS OF MATRIX ELEMENTS AND THEIR INTEGRALS
The study of EPIs from first principles requires evaluating Brillouin-zone integrals of functions that exhibit strong fluctuations. This requirement can be appreciated by inspecting Eqs. (145) and (157): there the denominators become large whenever the difference between two electronic eigenvalues approaches a phonon energy. As a result, while in DFT total energy calculations the Brillouin zone is typically discretized using meshes of the order of 10×10×10 points, the numerical convergence of EPI calculations requires much finer grids, sometimes with as many as 10 6 wavevectors (Giustino et al., 2007a;Poncé et al., 2015). Determining vibrational frequencies ω qν and perturbations ∆ qν v KS (r) for such a large number of wavevectors is a prohibitive task, since every calculation is roughly as expensive as one total energy minimization.
These difficulties stimulated the development of specialized numerical techniques for making calculations of EPIs affordable. In the following sections two such techniques are reviewed: electron-phonon Wannier interpolation and Fermi-surface harmonics.

Maximally-localized Wannier functions
In addition to the standard description of electrons in solids in terms of Bloch waves, as in Eq. (29), it is possible to adopt an alternative point of view whereby electrons are described as linear combinations of localized orbitals called 'Wannier functions' (Wannier, 1937). The most general relation between Wannier functions and Bloch waves can be written as follows. One considers electron bands ε nk with eigenfunctions ψ nk , where the index n is restricted to a set of bands that are separated from all other bands by finite energy gaps above and below. These bands are referred to as 'composite energy bands' (Marzari and Vanderbilt, 1997). Wannier functions are defined as: where U nmk is a unitary matrix in the indices m and n. From this definition and Eq. (A1) it follows that Wannier functions are normalized in the supercell, w mp |w m p sc = δ mp,m p . Furthermore, since u nk is lattice-periodic, Wannier functions have the property w mp (r) = w m0 (r − R p ). The inverse transformation of Eq. (175) is obtained by using the unitary character of U nmk together with Eq. (A1): The unitary matrix U nmk is completely arbitrary, therefore there exists considerable freedom in the construction of Wannier functions. For example by requiring that U nm,−k = U * nmk one can make Wannier functions realvalued. Marzari and Vanderbilt exploited this degree of freedom to construct Wannier functions that are maximally localized.
A comprehensive and up-to-date review of the theory and applications of maximally-localized Wannier functions (MLWFs) is given by Marzari et al. (2012). Here we only recall that, in order to minimize the spatial extent of a function in a periodic solid, one needs to use a modified definition of the position operator, since the standard position operator is unbounded in an infinite crystal. This procedure is now well-established and it is linked to the development of the modern theory of dielectric polarization (King-Smith and Vanderbilt, 1993;Resta, 1994). Nowadays it is possible to determine ML-WFs routinely (Mostofi et al., 2008). The original algorithm of Marzari and Vanderbilt (1997) was also extended to deal with situations where a composite set of bands cannot be identified. This happens notably in metals for electronic states near the Fermi energy. For these cases, Souza et al. (2001) developed a band 'disentanglement' procedure, which extracts a subset of composite bands out of a larger set of states.
For the purposes of the present article, the most important property of MLWFs is that they are exponentially localized in insulators, in the sense that |w m0 (r)| ∼ |r| −α exp(−h|r|) for large |r|, with α, h > 0 real parameters. This property was demonstrated in one spatial dimension by Kohn (1959a) and He and Vanderbilt (2001), and in two and three dimensions by Brouder et al. (2007), under the condition that the system exhibits time-reversal symmetry. In the case of metallic systems, no exponential localization is expected. However, the Wannier functions obtained in metals using the disentanglement procedure of Souza et al. (2001) are typically highly localized.
MLWFs are usually comparable in size to atomic orbitals, and this makes them ideally suited for Slater-Koster interpolation of band structures, as shown by Souza et al. (2001). This concept was successfully employed in a number of applications requiring accurate calculations of band velocities, effective masses, density of states, Brillouin-zone integrals, and transport coefficients (Wang et al., 2006;Yates et al., 2007;Wang et al., 2007;Pizzi et al., 2014).

Interpolation of electron-phonon matrix elements
Wannier functions were introduced in the study of EPIs by Giustino et al. (2007a,b). The starting point is the definition of the electron-phonon matrix element in the Wannier representation: 5 where the subscript 'sc' indicates that the integral is over the BvK supercell. The relation between these quantities and the standard EPI matrix elements g mnν (k, q) is found by replacing Eq. (176) inside Eq. (38), and using Eqs. (34)-(35) (Giustino et al., 2007a): where we defined u κα,qν = ( /2M κ ω qν ) 1 2 e κα,ν (q) and e κα,ν (q) are the vibrational eigenmodes of Eq. (15). The inverse relation is: with u −1 κα,qν = ( /2M κ ω qν ) − 1 2 e * κα,ν (q). The last two equations define a generalized Fourier transform of the electron-phonon matrix elements between reciprocal space and real space. In Eq. (179) we have N p and N p to indicate that the BvK supercells for electronic band structures and phonon dispersions may not coincide.
If the quantity g mnκα (R p , R p ) decays rapidly as a function of |R p | and |R p |, then only a small number of matrix elements in the Wannier representation will be sufficient to generate g mnν (k, q) anywhere in the Brillouin zone by means of Eq. (178). The dependence of the matrix elements on R p and R p can be analyzed by considering the following bound: |g mnκα (R p , R p )| ≤ sc dr |w * m0 (r)w n0 (r−R p )|× sc dr |∂V KS /∂τ κα (r−R p )|. The first term guarantees that the matrix element decays in the variable R p at least as fast as MLWFs. As a result the worst case scenario corresponds to the choice R p = 0. In this case, the matrix element |g mnκα (0, R p )| decays with the variable R p at the same rate as the screened electric dipole potential generated by the atomic displacement ∆τ κα . In non-polar semiconductors and insulators, owing to the analytical properties of the dielectric matrix (Pick et al., 1970), this potential decays at least as fast as a quadrupole, that is |R p | −3 . As a result, all matrix elements in reciprocal space are finite for q → 0 (Vogl, 1976) and hence amenable to interpolation. In the case of metals the asymptotic trend of ∂V KS /∂τ κα is dictated by Fermi-surface nesting, leading to Friedel oscillations that decay as |R p | −4 (Fetter and Walecka, 2003, pp. 175-180). These oscillations are connected to the Kohn anomalies in the phonon dispersion relations (Kohn, 1959b). In practical calculations, Friedel oscillations are usually not an issue since they are suppressed by the numerical smearing of the Fermi-Dirac occupations, and a Yukawa-type exponential decay is recovered. The case of polar materials is more subtle and will be discussed in Sec. VI.A.3. Figure 5 illustrates the spatial decay of |g mnκα (R p , R p )| as a function of R p and R p for the prototypical case of diamond.
The interpolation strategy is entirely analogous to standard techniques for generating phonon dispersion relations using the interatomic force constants (Gonze and Lee, 1997): one first determines matrix elements in the Bloch representation using DFPT on a corse grid in the Brillouin zone, as in Sec. III.B.3. Then ML-WFs are determined using the procedures of Marzari and Vanderbilt (1997) and Souza et al. (2001). This yields the rotation matrices U mnk to be used in Eq. (179). The Fourier transform to real space is performed via Eq. (179). At this point, one assumes that matrix elements outside of the Wigner-Seitz supercell defined by the coarse Brillouin-zone grid can be neglected, and uses Eq. (178) in order to obtain the matrix elements g mnν (k, q) on very fine grids. The last step requires the knowledge of the rotation matrices U mnk also on the fine grids; these matrices are obtained from the Wannier interpolation of the band structures, as described by Souza et al. (2001). The operation is computationally inexpensive and enables the calculation of millions of electron-phonon matrix elements. The procedure can now be applied routinely (Noffsinger et al., 2010;Poncé et al., 2016). Figure 6 shows the matrix elements obtained using this method, as compared to explicit DFPT calculations. FIG. 6 (Color online) Comparison between Wannierinterpolated electron-phonon matrix elements and explicit DFPT calculations, for diamond. The interpolated matrix elements were calculated starting from a coarse 4 3 Brillouinzone grid (dotted line, black), a 6 3 grid (dashed line, blue), and a 8 3 grid (solid line, red). The dots indicate explicit DFPT calculations. In this example |nk is set to the valence band top at Γ; |mk + q spans Λ3, ∆5, and Σ2 bands, and the phonon is set to the highest optical branch. Reproduced with permission from (Giustino et al., 2007a), copyright (2007) by the American Physical Society.
Wannier interpolation of electron-phonon matrix elements was successfully employed in a number of applications, ranging from metal and superconductors to semi-conductors and nanoscale systems. 6

Electron-phonon matrix elements in polar materials
In the case of polar materials, that is systems exhibiting nonzero Born effective charges (Pick et al., 1970), the interpolation scheme discussed in Sec. VI.A.2 breaks down. In fact, in these systems the dominant contribution to the potential ∂V KS /∂τ κα in Eq. (177) is a dipole, which decays as |R p | −2 . As a consequence some of the matrix elements in reciprocal space diverge as |q| −1 for q → 0, and cannot be interpolated straightforwardly from a coarse grid to a fine grid. Physically this singularity corresponds to the 'Fröhlich electron-phonon coupling' (Fröhlich, 1954).
The adaptation of the Wannier interpolation method to the case of polar materials was recently given by Sjakste et al. (2015) and by Verdi and Giustino (2015). In both works the basic idea is to separate the matrix elements into a short-range contribution, g S mnν (k, q), which is amenable to standard Wannier interpolation, and a long-range contribution, g L mnν (k, q), which is singular and is dealt with analytically. The strategy is analogous to that in use for calculating LO-TO splittings in polar materials (Gonze and Lee, 1997). The starting point is to define the long-range component of the matrix elements by considering the potential generated by the Born charges of all the atoms, when displaced according to a given vibrational eigenmode. The derivation relies on standard electrostatics and can be found in (Verdi and Giustino, 2015): In this expression Z * κ and ∞ denote the Born effective charge tensors and the electronic permittivity tensor (that is, the permittivity evaluated at clamped nuclei). This expression is the generalization of Fröhlich's model to the case of anisotropic crystalline lattices and multiple phonon modes (Fröhlich, 1954). The result can be derived alternatively using the analytical properties of the dielectric matrix (Pick et al., 1970) as discussed by Vogl (1976). The dots correspond to explicit DFPT calculations. The red dashed line is the short-range component of the matrix elements, g S . The solid curve in blue represents the matrix elements g S + g L , as obtained from the modified Wannier interpolation of Sec. VI.A.3. The interpolation was performed starting from a coarse 4×4×4 unshifted grid. Reproduced with permission from (Verdi and Giustino, 2015), copyright (2015) by the American Physical Society.
In order to perform Wannier interpolation, one subtracts Eq. (180) from the matrix elements computed on a coarse grid, interpolates the remainig short-range part, and then adds back Eq. (180) on the fine grid. This process requires the interpolation of the brakets · · · sc in the second line of Eq. (180). Verdi and Giustino showed that, for small q + G, these brakets can be interpolated via the relation ψ mk+q |e i(q+G)·r |ψ nk sc = U k+q U † k mn , where the rotation matrices U mnk are obtained as usual from the procedure of Marzari and Vanderbilt (1997) and Souza et al. (2001). 7 Figure 7 shows an example of Wannier interpolation for the prototypical polar semiconductor TiO 2 : it is seen that the singularity is correctly captured by the modified interpolation method.
At the end of this section, we mention that other interpolation schemes are equally possible (Eiguren and Ambrosch-Draxl, 2008b;Prendergast and Louie, 2009;Agapito et al., 2013;Gunst et al., 2016). For example Eiguren and Ambrosch-Draxl (2008b)  orbitals, and all the concepts discussed in this section remain valid. An interpolation scheme using local orbitals was recently demonstrated by Gunst et al. (2016).

B. Fermi surface harmonics
In the study of metallic systems, one is often interested in describing EPIs only for electronic states in the vicinity of the Fermi surface. In these cases, besides the Wannier interpolation discussed in Secs. VI.A.2-VI.A.3, it is possible to perform efficient calculations using 'Fermi-surface harmonics' (FSH). FSHs were introduced by Allen (1976) and recently revisited by Eiguren and Gurtubay (2014).
The basic idea underlying FSHs is to replace expensive three-dimensional Brillouin-zone integrals by inexpensive one-dimensional integrals in the energy variable. To this aim, Allen (1976) proposed to expand functions of the band index n and wavevector k, say A nk , in products of pairs of functions, one depending on the energy, A L (ε), and one depending on the wavevector, Φ L (k): In this expression, the Fermi-surface harmonics Φ L (k) (to be defined below) are constructed so as to obey the following orthogonality condition: where N (ε) = N −1 p nk δ(ε nk −ε) is the density of states. Using Eqs. (181)-(182) one finds: Allen showed that in the FSH representation a linear system of the kind A nk = N −1 p n k M nk,n k B n k transforms into A L (ε) = L dε N (ε )M LL (ε, ε )B L (ε ). Linear systems of this kind are common in the solution of the Boltzmann transport equation (Sec. X) and the Eliashberg equations for the superconducting gap (Sec. XI.B). If one could perform the expansion using only a few harmonics, then the transformation would be advantageous, since the integrals over the wavevectors would have been absorbed in the expansion coefficients.
In the original proposal of Allen, the harmonics Φ L (k) were defined as polynomials in the band velocities, however the completeness of the basis set was not established. In a recent work, Eiguren and Gurtubay (2014) proposed to construct these functions as eigenstates of a modified Helmholtz equation: where v k = −1 ∇ k ε nk is the band velocity for states at the Fermi surface, and ω L is the eigenvalue for the harmonic Φ L . The new definition in Eq. (184) maintains the properties of the original FSHs, and carries the added advantage that the basis set is complete. In this case the subscript L in Φ L (k) labels the eigenstates of the Helmholtz equation. Eiguren and Gurtubay demonstrated the construction of 'Helmholtz FSHs' for prototypical metals such as Cu, Li,and MgB 2 . Recent examples of the application of Fermi surface harmonics to first-principles calculations of EPIs include work on the photoemission kink of YBa 2 Cu 3 O 7 (Heid et al., 2008), and on the Seebeck coefficient of Li (Xu and Verstraete, 2014).

VII. NON-ADIABATIC VIBRATIONAL FREQUENCIES AND LINEWIDTHS
As discussed in Secs. V.A.2-V.A.3, the electron-phonon interaction can lead to a renormalization of the adiabatic vibrational frequencies and to a broadening of the spectral lines.
The first ab initio investigations of the effects of the non-adiabatic renormalization of phonon frequencies were reported by  and Pisana et al. (2007). In these works the authors concentrated on the E 2g phonon of graphene, which is found at the wavenumber ω/2πc = 1585 cm −1 at room temperature (c is the speed of light). This phonon corresponds to an in-plane C-C stretching vibration with q = 0, and has been studied extensively via Raman spectroscopy. In the graphene literature this mode is referred to as the 'Raman G band'. Figure 8 shows a comparison between calculated and measured E 2g phonon frequencies, as a function of doping, from (Pisana et al., 2007). The calculations were performed (i) within the adiabatic approximation and (ii) by including the non-adiabatic frequency renormalization using Eq. (145). 8 From Fig. 8 we see that the adiabatic theory is unable to reproduce the experimental data. On the contrary, the calculations including non-adiabatic effects nicely follow the measured Raman shift. This is a clear example of the limits of the adiabatic Born-Oppenheimer approximation and a demonstration of the importance of the phonon self-energy in Eq. (145).
The fact that the adiabatic approximation is inadequate for the E 2g phonon of doped graphene should have been expected from the discussion on p. 27. In fact, graphene is a zero-gap semiconductor, therefore electrons residing in the vicinity of the Dirac points can make 'virtual' transitions with |q| = 0 and energies comparable to that of the E 2g mode. As a result, the condition underlying the adiabatic approximation, |ε mk+q − ε nk | ω qν , does not hold in this case.
The importance of non-adiabatic effects was confirmed also in the case of metallic single-walled carbon nanotubes (Piscanec et al., 2007;Caudal et al., 2007). In these works, the authors studied the phonon dispersion relations in the vicinity of |q| = 0, and they found that the difference between adiabatic and non-adiabatic dispersions is concentrated around the zone center. This finding is consistent with earlier models of non-adiabatic effects. In fact Maksimov and Shulga (1996) showed that, for metals with linear electron bands crossing the Fermi level, Π NA is only significant for wavevectors |q| ∼ ω/v F , where ω is the phonon energy and v F the Fermi velocity. This result can be derived from Eq. (145).
In the previous examples, the non-adiabatic renormalization of the vibrational frequencies is measurable but very small, typically of the order of 1% of the corresponding adiabatic frequencies. Saitta et al. (2008) considered the question as to whether one could find materials exhibiting large non-adiabatic renormalizations, and considered several graphite intercalation compounds, namely LiC 6 , LiC 12 , KC 8 , KC 24 , RbC 8 , CaC 6 , SrC 6 , BaC 6 , as well as other metallic systems such as MgB 2 , Mg, and Ti. In order to calculate the non-adiabatic renormalization at a reduced computational cost, the real part of the phonon self-energy was approximated as follows: Re Π NA q=0,νν N F |g nnν (k, q = 0)| 2 BZ , where · · · BZ stands for the average taken over the wavevectors k in the Brillouin zone, and N F is the density of states at the Fermi level. This expression can be derived from Eq. (145) by replacing the bare matrix elements by their screened counterparts, and by neglecting the 'interband' contributions m = n in the sum. Figure 9 shows a comparison between vibrational frequencies from experiment and those calculated with or without including the nonadiabatic self-energy. It is clear that the non-adiabatic frequencies are in much better agreement with experiment than the corresponding adiabatic calculations. Furthermore, in these compounds the renormalization can reach values as large as ∼300 cm −1 . In contrast to this, the renormalization in Mg and Ti was found to be of only a few wavenumbers in cm −1 .
The case of MgB 2 proved more puzzling: here the calculated non-adiabatic frequency is 761 cm −1 , while experiments reported 600 cm −1 . In order to explain this discrepancy, Saitta et al. reasoned that a more accurate calculation would require taking into account the relaxation time of the electrons, as pointed out by Maksimov and Shulga (1996). This would act so as to partly wash out non-adiabatic effects. In the field-theoretic language of Sec. V, this observation corresponds to stating that when one approximates the dressed electron propagator G in Fig. 2(c) using the non-interacting Kohn-Sham Green's function, one should include at the very least the effects of finite electron lifetimes (due to electronelectron, electron-impurity, and electron-phonon scattering), for example via the imaginary part of Eq. (157).
The calculations discussed so far in this section addressed the non-adiabatic renormalization of zone-center phonons. The generalization to calculations of complete phonon dispersions was made by Calandra et al. (2010). In their work, Calandra et al. employed Wannier interpolation (see Sec. VI) in order to calculate the nonadiabatic phonon self-energy of Eq. (145) throughout the Brillouin zone. Figure 10 shows a comparison between the standard DFPT phonon dispersion relations of CaC 6 and the dispersions obtained after incorporating the nonadiabatic self-energy. It is seen that also in this case nonadiabatic effects are most pronounced at small q, and can be as large as 7% of the adiabatic frequency.
In their work Calandra et al. approximated the bare matrix element g b mnν (k, q) appearing in Eq. (145) by the screened matrix element g mnν (k, q). This replacement was justified by reasoning that the error is of secondorder in the induced electron density, hence it should be negligible.
The broadening of vibrational spectra arising from the electron-phonon interaction is almost invariably calculated from first principles using Eq. (146). Since the integration of the Dirac delta is computationally costly, it is common to rewrite that equation by neglecting the phonon energy in the delta function and by taking the low-temperature limit, as proposed by Allen (1972b): where ε F is the Fermi energy. Oddly enough, this is a sort of adiabatic approximation to the non-adiabatic theory. The main advantage is that this expression is positive definite, hence easier to converge numerically as compared to the complete expression in Eq. (146). The disadvantages are that the temperature dependence is lost, and that one cannot resolve fine features on the scale of a phonon energy. There exists a vast literature on first-principles calculations of phonon linewidths using Eq. (185), mostly in connection with electron-phonon su-perconductors. 9 Equation (185) is now implemented in several large software projects (Giannozzi et al., 2009;Gonze et al., 2009), and it is used routinely.
Phonon linewidths range from very small values such as ∼1 meV in Nb (Bauer et al., 1998) to large values such as ∼30 meV in MgB 2 (Shukla et al., 2003). The agreement between calculations and neutron scattering or Raman measurements is usually reasonable.
Calculations of phonon linewidths using the more accurate expression in Eq. (146) are computationally more demanding and have been reported less extensively in the literature. 10 So far we considered the effect of EPIs on the frequencies and lifetimes of vibrational excitations in solids. Another important phenomenon which modifies frequencies and lifetimes is anharmonicity. Anharmonic effects result from third and higher-order terms in the Taylor expansion of the total potential energy U in the atomic displacements (Sec. III.A). These effects can be interpreted as additional interactions that couple the oscillators of the harmonic lattice; for example, third-order anharmonic effects reduce phonon lifetimes via energy upor down-conversion processes involving three phonons.
Anharmonic effects can be described using a manybody field-theoretic formalism (Cowley, 1963), in complete analogy with the theory of EPIs discussed in Sec. IV. The calculation of anharmonic effects from first principles goes through the evaluation of thirdand fourth-order derivatives of the total potential energy U in the adiabatic approximation. Third-order coefficients are routinely computed using DFPT (Debernardi et al., 1995;Deinzer et al., 2003;Lazzeri et al., 2003;Broido et al., 2007;Bonini et al., 2007). In thoses cases where the harmonic approximation fails completely, 'selfconsistent phonon' techniques can be employed (Hooton, 1955;Koehler, 1966); recent implementations and calculations can be found in Errea et al. (2011Errea et al. ( , 2014, Hellman et al. (2011Hellman et al. ( , 2013; and Monserrat et al. (2013).

VIII. ELECTRON-PHONON INTERACTIONS IN PHOTOELECTRON SPECTROSCOPY
In Sec. V.B.5 we have seen how the electron-phonon interaction in metals can lead to band structure 'kinks' near the Fermi energy. This is illustrated by the model calculation in Fig. 4(a). The experimental investigation of these features started in the late 1990s, following the 9 See for example early frozen-phonon calculations (Dacorogna et al., 1985b;Chang and Cohen, 1986;Lam et al., 1986) and more recent DFPT calculations (Savrasov et al., 1994;Bauer et al., 1998;Shukla et al., 2003;Heid et al., 2010). Earlier calculations not based on DFT are reviewed by Grimvall (1981). development of high-resolution angle-resolved photoelectron spectroscopy (ARPES). Since in ARPES only the component of the photoelectron momentum parallel to the sample surface is conserved (Damascelli et al., 2003), complete energy vs. wavevector dispersion relations can be measured directly only for 2D or quasi-2D materials. Accordingly, the first observations of kinks were reported for the surface states of elemental metals 11 and for the CuO 2 planes of copper oxide superconductors. 12 Ab initio calculations of ARPES kinks can be performed by using the diagonal part of the Fan-Migdal self-energy (the Debye-Waller self-energy will be discussed at the end of this section). To this aim, it is common to rewrite Eqs. (157)-(158) at finite temperature using a spectral representation: Here the function α 2 F is the so-called the 'Eliashberg function' and is defined as: This quantity is positive, temperature-independent, and contains all the materials-specific parameters. Physically it is proportional to the scattering rate of an electron in the state |nk into any electronic state at the energy ε, by emitting or absorbing any one phonon of energy ε . One complication related to the Eliashberg function is that in the literature many variants of Eq. (187) can be found, each stemming from specific approximations; some of these expressions are summarized by Grimvall (1981, pp. 107-109). The first ab initio calculations of the phonon-induced electron self-energies and photoemission kinks were reported by Eiguren et al. (2003a) for the surface state of the Be(0001) surface. Since the evaluation of Eqs. (186) and (187) is computationally demanding, Eiguren et al. employed simplified expressions which involve three approximations: (i) the Eliashberg function is replaced by its isotropic average, α 2 F n (ε, ε ) = dk δ(ε nk − ε) α 2 F nk (ε, ε )/ dk δ(ε nk −ε); (ii) phonon energies are neglected next to electron energies (as in the adiabatic approximation), and (iii) particle-hole symmetry is assumed. Using these approximations the imaginary part  (0001) surface. (a) Imaginary part of the self-energy, obtained from Eq. (188). The dashed (black) line is the self-energy evaluated using the DFT/LDA bands; the solid lines (color/grayscale) correspond to the selfenergy calculated by taking into account the renormalization of the DFT band structure by the electron-phonon interaction. (b) Real part of the self-energy, using the same color/grayscale code as in (a). The inset in (a) compares the renormalized band structure (color) with the 'bare' DFT band (black dashed line). The inset of (b) shows the renormalization of the band velocity induced by the electron-phonon interaction. Reproduced with permission from (Eiguren et al., 2003a(Eiguren et al., ), copyright (2003 by the American Physical Society. of the Fan-Migdal self-energy becomes: where the average of the self-energy is defined as for the Eliashberg function. The real part of the self-energy can be found starting from the same approximations, and is given by Grimvall (1981). Figure 11 shows the self-energy of the surface state at the Be(0001) surface calculated by Eiguren et al. (2003a) using Eq. (188). The imaginary part resembles a stepfunction, with an onset around the energy threshold for phonon emission by a hole (40-80 meV in this case). At a qualitative level, this trend can be rationalized by replacing the Eliashberg function in Eq. (188) by a Dirac delta at a characteristic phonon energy ω ph . In this case, the hole self-energy becomes proportional to f [( ω + ω ph )/k B T ]. At T = 0 this is precisely a step function with onset at − ω ph . The real part of the selfenergy vanishes for |ω| ω max , with ω max being the largest phonon frequency. This can be seen in Fig. 11(b), and is a consequence of the approximation of particle-hole symmetry. Eiguren et al. also determined the renormalization of the surface state band structure arising from electron-phonon interactions, using Eq. (164); this is shown in the inset of Fig. 11(a). Overall the calculations of Eiguren et al. (2003a) showed good agreement with photoelectron spectroscopy experiments, both in the shape and magnitude of the self-energy (LaShell et al., 2000).
In the case of complex systems the validity of the approximations (i)-(iii) leading to Eq. (188) is not warranted, and a direct calculation of the Fan-Migdal selfenergy using Eqs. (186) and (187) is necessary. Calculations of the complete self-energy were reported by Park et al. (2007) for graphene, by Giustino et al. (2008) for the high-temperature superconductor La 1−x Sr x Cu 2 O 4 , and by Margine et al. (2016) for Ca-intercalated bilayer graphene. Figure 12 shows the calculated self-energy and spectral function of graphene calculated by Park et al. (2007Park et al. ( , 2009). The structure of the self-energy is similar to that shown in Fig. 11, with one important exception: Σ FM nnk (ω) does not vanish a few phonon energies away from the Fermi level, but tends instead towards a linear asymptote. Calandra and Mauri (2007) performed a combined ab initio/analytical study of the effects of the electron-phonon interaction on the electron bands of graphene and obtained very similar results. A linear asymptote in the real-part of the self-energy is a general feature of systems which do not exhibit particle-hole symmetry. For another example see the work on copper oxides by Giustino et al. (2008). Using the Fan-Migdal self-energy, it is possible to calculate the renormalization of the band velocity induced by the electron-phonon interaction. Let us denote by v nk = −1 ∇ k ε nk the DFT band velocity and V nk = −1 ∇ k E nk the band velocity after including electronphonon interactions. Using Eq. (164) with Γ nk = 0 we find that these two quantities are simply related by V nk = Z nk v nk = v nk /(1 + λ nk ), where Z nk is the quasiparticle strength of Eq. (174), and λ nk is the 'mass enhancement parameter' or 'electron-phonon coupling strength' (Grimvall, 1981): In the study of EPIs in metals, the coupling strength λ nk is of significant interest since it is related to the superconducting transition temperature of phonon-mediated superconductors (see Secs. VIII.1 and XI).
The velocity renormalization in graphene calculated using Eq. (189) is shown in Fig. 12(c) and (d), while a calculation of the complete spectral function A(k, ω) is shown in Fig. 12(e). Here the characteristic photoemission kink is visible between 150-200 meV but it is not very pronounced, since in this case λ nk ∼ 0.1 (Park et al., 2009). These results are in good agreement with measured photoelectron spectra (Bostwick et al., 2007;Zhou et al., 2007).
Incidentally, we remark that in the analysis of ARPES data it is common to extract the coupling strength λ nk directly from the ratio of the band velocities above and below the electron-phonon kink. However, this procedure is subject to a significant uncertainty, since the 'bare' velocity is not known and must be approximated by fitting the fully-interacting dispersions using ad hoc models. For example, in the vicinity of Van Hove singularities this procedure leads to a significant overestimation of the electron-phonon coupling strength λ nk (Park et al., 2008b;Bianchi et al., 2010).
In addition to photoemission kinks, recent ARPES measurements revealed the existence of polaron satellites in doped oxides, namely TiO 2 (Moser et al., 2013) and SrTiO 3 (Cancellieri et al., 2016;Chen et al., 2015;Wang et al., 2016). The phenomenology is similar to what was discussed in relation to Fig. 4(b). The first theoretical studies along this direction were reported by Story et al. (2014), who applied the cumulant expansion approach to the case of the electron-phonon self-energy; by Antonius et al. (2015) who identified satellites in the spectral functions of LiF and MgO; and by Verdi et al. (2016), who calculated the ARPES spectra of n-doped TiO 2 .
At the end of this section, it is worth coming back to the Debye-Waller self-energy. So far we only discussed the Fan-Migdal self-energy, starting from Eq. (186), and we ignored the Debye-Waller self-energy appearing in Eq. (152). This omission reflects the fact that, in the literature on electron-phonon interactions in metals, the DW term has always been disregarded. In order to rationalize this approximation, we rewrite the DW self-energy as follows, by combining Eqs. (160) and (40): with V DW (r) = Ω −1 BZ ν dq (n qν +1/2)∆ qν ∆ −qν v KS (r). The subscript 'uc' indicates that the integral is over one unit cell. From Eq. (190) we see that V DW acts like a static local potential; indeed the first calculations including DW effects were performed by directly modifying the ionic pseudopotentials (Antoncík, 1955;Walter et al., 1970). From Eq. (190) we also see that the only variation in the DW self-energy comes from the Bloch amplitudes u nk . Let us consider the limiting situation of the homogeneous electron gas. In this case |u nk (r)| 2 = 1/Ω (Sec. III.B.5), therefore Σ DW nnk is a constant, independent of k. This behavior should be contrasted with the Fan-Migdal self-energy, which exhibits significant structure near the Fermi energy, as it can be seen in Fig. 4.
In more realistic situations, such as doped semiconductors, k · p perturbation theory (Kittel, 1963) can be used to show that Σ DW nnk varies smoothly as a function of k within the same band. In contrast with this scenario, Σ DW nnk exhibits large variations across different bands. This carries important consequences for the calculation of temperature-dependent band gaps (Sec. IX.A.1).

Electron mass enhancement in metals
We now come back to the mass enhancement parameter λ nk introduced in Eq. (189), since this quantity played a central role in the development of the theory of EPIs in metals.
The notion of 'mass enhancement' becomes clear when we consider a parabolic band as in the model calculations of Fig. 4. Near the Fermi surface the non-interacting dispersions are given by ε nk = k F · (k − k F )/m * , where k F is a wavevector on the Fermi surface, and the electron velocity is v nk = k F /m * . After taking into account the EPI, the electron velocity is renormalized to V nk = v nk /(1 + λ nk ). Since the magnitude of the Fermi momentum is unchanged (see caption of Fig. 4) this renormalization can be interpreted as an effective increase of the band mass from m * to m ep = m * (1 + λ nk ). This reasoning holds for metals with parabolic bands and for doped semiconductors near band extrema, and does not take into account the Debye-Waller self-energy.
The electron mass enhancement is reflected into the increase of the heat capacity of metals at low temperature. In fact, below the Debye temperature the electronic contribution to the heat capacity dominates over the lattice contribution (Kittel, 1976). Since the heat capacity is proportional to the density of states at the Fermi level, and the density of states is inversely proportional to the band velocity, it follows that the heat capacity is directly proportional to the electron mass. This property can be used as a means to determine the strength of the electron-phonon coupling in simple metals from specific heat measurements (Grimvall, 1975).
The general theory of the effects of electron-phonon interactions on the heat capacity and other thermodynamic functions was developed by Eliashberg (1963), Prange andKadanoff (1964), andGrimvall (1969). A field-theoretic analysis of the effect of EPIs on thermodynamic functions was developed by Eliashberg (1963) starting from the identities of Luttinger and Ward (1960) in the zero-temperature limit. Eliashberg's analysis was subsequently extended to all temperatures by Grimvall (1969). Here we only quote Grimvall's result relating the electronic entropy to the Fan-Migdal self-energy of Sec. V.B.1: (191) In order to derive this relation, Grimvall started by expressing the thermodynamic potential of the coupled electron-phonon system in terms of the electron and phonon propagators and self-energies, and identified the electronic contribution by neglecting terms of order (m e /M 0 ) 1/2 as well as electron-electron interactions (Grimvall, 1969, Appendix). Below the Debye temperature, an explicit expression for the entropy in Eq. (191) can be obtained by noting that the function cosh −2 ( ω/2k B T ) is nonvanishing only for ω 2k B T ; in this range Eq. (189) yields Σ FM nnk (ω, T ) −λ nk ω, therefore the integration in Eq. (191) can be carried out explicitly. As a result, the low-temperature heat capacity can be written as: If we ignore the EPI by setting λ nk = 0, this expression reduces to the standard textbook result for the free electron gas (Kittel, 1976). At high temperature Eq. (192) is no longer valid, and one has to evaluate the integral in Eq. (191) using the complete frequency-dependent FM self-energy. The main result of this procedure is that at high temperature the electronic heat capacity is no longer renormalized by EPIs. A detailed discussion of this aspect is provided by Grimvall (1969;1981). Early examples of DFT calculations of mass enhancement parameters and comparisons with specific heat measurements in simple metals can be found in Dacorogna et al. (1985b), Savrasov et al. (1994, and Liu and Quong (1996).

IX. ELECTRON-PHONON EFFECTS IN THE OPTICAL PROPERTIES OF SEMICONDUCTORS
A. Temperature dependence of band gaps and band structures

Perturbative calculations based on the Allen-Heine theory
In Sec. V.B.3 we discussed how the electron-phonon interaction induces a 'renormalization' of the electronic energy levels, and thereby gives rise to 'temperaturedependent band structures'. These effects have been studied in detail using the Fan-Midgal and the Debye-Waller self-energies, either via the Raleigh-Schrödinger approximation to Eq. (166), or via its adiabatic counterpart given by Eq. (167).
Equation (167) was first employed in a number of calculations based on empirical pseudopotentials, following the seminal work of Allen and Heine (1976). The key references can be found on p. 31, footnote. Allen and Cardona (1981) offer a clear introduction to the basic theory, a discussion of the computational methodology, as well as an historical perspective on earlier calculations.
The evaluation of the Debye-Waller contribution to the self-energy requires the calculation of the second-order variations of the Kohn-Sham potential with respect to the ionic displacements, Eq. (40). From a computational standpoint, this is challenging because one would have to use second-order DFPT, as discussed at the end of Sec. III.B.3. In order to avoid this complication, it is common practice to recast all second-order derivatives as products of first-order derivatives. This strategy was introduced by Allen and Heine in the case of monoatomic crystals, and extended to polyatomic unit cells by Allen and Cardona. The key observation behind this approach is that one can impose translational invariance of the theory to second order in the nuclear displacements. Specifically, the variation of the Kohn-Sham eigenvalues ensuing an arbitrary displacement of the nuclei should not change if all nuclei are further displaced by the same amount. Using time-independent perturbation theory, this condition yields the following two sum rules: Here ∂ καp V KS is a short-hand notation for ∂V KS /∂τ καp and similarly for the second derivative; the primed summation indicates that eigenstates ψ mk+q such that ε nk = ε mk+q are skipped. The first sum rule is equivalent to stating that the electron-phonon matrix elements g mnν (k, q) associated with the three translational modes at |q| = 0 must vanish; this is an alternative formulation of the 'acoustic sum rule'. The second sum rule, Eq. (194), suggests to express the matrix elements of the second derivatives of the potential in terms of firstorder derivatives. However, Eq. (194) cannot be used as it stands, since it involves a sum of matrix elements on the left-hand side. In order to proceed, Allen and Heine employed the 'rigid-ion' approximation, whereby V KS is written as a sum of atom-centered contributions (see Sec. II.A.1). Under this approximation all the terms κp = κ p on the left-hand side of Eq. (194) are neglected, and an explicit expression for ψ nk |∂ 2 καp,καp V KS |ψ nk sc is obtained.
In view of practical DFT implementations, Giustino et al. (2010) used the sum rule in Eq. (194) in order to rewrite the Debye-Waller self-energy as follows: Here g DW mnν (k, q) is an 'effective' matrix element, and it is obtained from the standard DFPT matrix elements by means of inexpensive matrix multiplications: In the case of the three translational modes at |q| = 0, these definitions are replaced by g DW mnν (k, q) = 0, see the discussion at the end of p. 57. The derivation of Eqs. The first ab initio calculations using the formalism of Allen and Heine were reported by Marini (2008), who investigated the temperature dependence of the optical absorption spectrum of silicon and boron nitride. In this work Marini included excitonic effects by combining the Bethe-Salpeter formalism (Onida et al., 2002) with the Allen-Heine theory, and obtained good agreement with experiments by calculating the direct absorption peaks using DFT/LDA phonons and matrix elements (indirect optical absorption will be discussed in Sec. IX.B).
The second application of the Allen-Heine theory using DFT/LDA was reported by Giustino et al. (2010) for the case of diamond. Here the authors investigated the temperature dependence of the direct band gap of diamond, and found that the Fan-Migdal and the Debye-Waller self-energies are of comparable magnitude. The calculations captured the characteristic Varshni effect (Fig. 3), and were able to reproduce the measured redshift of the band gap in the temperature range 80-800 K. These calculations were based on the adiabatic version of the  for the Debye-Waller self-energy. The calculations by Giustino et al. confirmed the large (> 0.5 eV) zero-point renormalization of the direct gap of diamond predicted by Zollner et al. (1992) using the empirical pseudopotential method.
The unusually large zero-point correction to the electronic structure of diamond stimulated further work on this system: Marini (2011, 2012) calculated the gap renormalization in diamond by employing both the adiabatic version of the Allen-Heine theory, as well as the non-adiabatic Green's function approach, as FIG. 13 (Color online) Temperature dependence of the direct band gap of diamond calculated using the Allen-Heine theory. The upper curve shows the results obtained within DFPT at the LDA level. The lower curve was obtained via GW calculations in the frozen-phonon approach. Triangles are experimental data. The zero-point renormalization calculated by including GW quasiparticle corrections is 628 meV. Reproduced with permission from (Antonius et al., 2014), copyright (2014) by the American Physical Society. described in Sec. V.B.5. Their calculations confirmed the large zero-point renormalization, and showed that the adiabatic theory underestimates the effect to some extent. Cannuccia and Marini (2012; also analyzed the quasiparticle renormalization and the spectral function. Antonius et al. (2014) revisited the electron-phonon interaction in diamond by assessing the reliability of the rigid-ion approximation and the importance of manybody GW quasiparticle corrections to the DFT/LDA band structure. The main findings were that the rigid-ion approximation introduces a very small error in diamond, of the order of ∼10 meV, while GW quasiparticle corrections can increase the electron-phonon renormalization of the band gap by as much as ∼200 meV. The temperature dependence of the band gap of diamond calculated by Antonius et al. is shown in Fig. 13.
Further work on diamond was reported by Poncé et al. (2014b), who compared ab initio calculations based on the Allen-Heine theory with explicit frozen-phonon calculations (Sec. IX.A.2). The corrections to the rigid-ion approximation were found to be smaller than 4 meV in all cases considered. Poncé et al. (2014a,b) also reported a detailed assessment of the accuracy of the various levels of approximation in the calculation of the zero-point renormalization of energy levels, as well as a thorough comparison between the results of different first-principles implementations.
The electron-phonon renormalization of band structures was also investigated in a number of other systems. For example Kawai et al. (2014) studied zinc-blend GaN by combining the Allen-Heine theory with the Bethe-Salpeter approach. Poncé et al. (2015) investigated silicon, diamond, BN, α-AlN, and β-AlN using both the adiabatic version of the Allen-Heine theory and the non-adiabatic Green's function method of Eqs. (157)-(158). Friedrich et al. (2015) investigated the zero-point renormalization in LiNbO 3 using the adiabatic Allen-Heine theory. Villegas et al. (2016) studied the anomalous temperature dependence of the band gap of black phosphorous. Antonius et al. (2015) investigated diamond, BN, LiF, and MgO, focusing on the dynamical aspects and the spectral function (see Sec. VIII). The works by Antonius et al. (2015) and Poncé et al. (2015) were the first to report complete band structures at finite temperature.
In  the authors paid special attention to the numerical convergence of the self-energy integrals with respect to the limit η → 0 of the broadening parameter; they noted that in the case of polar crystals the adiabatic correction to the electron energies of band extrema, as given by Eq. (167), diverges in the limit of dense Brillouin-zone sampling. This behavior stems from the polar singularity in the electron-phonon matrix elements, Eq. (180). In fact, near band extrema the bands are approximately parabolic, and the integrand in the adiabatic Fan-Migdal self-energy goes as q −4 for q → 0, while the volume element goes only as dq = 4πq 2 dq. This problem can be avoided by first performing the integration over q in principal value, without neglecting phonon frequencies, and then taking the limit ω qν → 0 so as to recover the adiabatic approximation (Fan, 1951, Sec. IV); in this way the adiabatic approximation can still be employed without incurring into a singularity in the calculations. A practical strategy to correctly perform the principal value integration in first-principles calculations was recently proposed by Nery and Allen (2016); here the authors treat the singularity via an explicity analytic integration near q = 0. The complications arising in polar materials can also be avoided at once by using directly the more accurate expression in Eq. (166) based on Brillouin-Wigner perturbation theory, or even better by calculating the spectral functions as in (Antonius et al., 2015;Kawai et al., 2014). In particular, Eq. (166) shows that in more accurate approaches the infinitesimal η should be replaced by the finite physical linewidth Γ nk .
Although temperature-dependent band structures of polar materials were recently reported (Antonius et al., 2015;Kawai et al., 2014;Poncé et al., 2015), the specific role of the Fröhlich coupling discussed in Sec. VI.A.3 received only little attention so far. The only ab initio investigations which specifically addressed the role of polar couplings in this context are from Botti and Marques (2013) and Nery and Allen (2016). In order to understand the strategy of Botti and Marques, we refer to the Hedin-Baym equations in Sec. IV.B.1. Botti and Marques proposed that, instead of splitting the screened Coulomb interaction W into electronic and nuclear contributions as in Eq. (104), one could try to directly calculate the screened Coulomb interaction W including the lattice screening, as in Eqs. (85) and (86). In order to make the calculations tractable, Botti and Marques evaluated the total dielectric matrices using a simplified model based on the Lyddane-Sachs-Teller relations. The resulting formalism combines GW calculations and experimentally measured LO-TO splittings. The zero-point renormalization of the band gaps calculated by Botti and Marques for LiF, LiCl, NaCl, and MgO were all > 1 eV. This is an interesting result and deserves further investigation. We note incidentally that the Allen-Heine theory and that of Botti and Marques can both be derived from the Hedin-Baym equations. Therefore the approach of Botti and Marques should effectively correspond to calculating the Fan-Migdal self-energy by retaining only the long-range part of the polar electron-phonon matrix elements. In the case of the work by Nery and Allen (2016), the authors reported a Fröhlich contribution to the zeropoint renormalization of the band gap of GaN of 45 meV, to be compared with the total renormalization arising from all modes of 150 meV.

Non-perturbative adiabatic calculations
An alternative approach to the calculation of temperature-dependent band structures consists of avoiding perturbation theory and electron-phonon matrix elements altogether, and replacing the entire methodology discussed in Sec. IX.A.1 by straightforward finite-differences calculations. To see how this alternative strategy works we perform a Taylor expansion of the Kohn-Sham eigenvalues ε nk to second order in the atomic displacements ∆τ καp , and then average the result on a nuclear wavefunction identified by the quantum numbers {n qν }. After using Eq. (20) one obtains: where we used the formal definition ∂/∂n qν = ∆ qν ∆ −qν , and the variations ∆ qν are the same as in Eqs. (33)-(35). The nuclear wavefunctions are obtained from the groundstate in Eq. (B10) by applying the ladder operators, as discussed in Appendix B. The above expression can be generalized to finite temperature by considering a canonical average over all possible nuclear states. The result maintains the same form as in Eq. (199), except that we now have the Bose-Einstein occupations n qν (T ) (see footnote 3). Equation (199) is precisely the conceptual starting point of the Allen-Heine theory of temperaturedependent band structures, and appeared for the first time in (Allen and Cardona, 1981). If the variations ∆ qν ∆ −qν ε nk are calculated in second-order perturbation theory, one obtains precisely the formalism of Allen and Heine. It has been proposed that the coefficients ∂ε nk /∂n qν could alternatively be obtained from the derivatives of the vibrational frequencies with respect to the electronic occupations, ∂ω qν /∂f nk (Allen and Hui, 1980;King-Smith et al., 1989;Poncé et al., 2014b). A formal derivation of the link between these alternative approaches can be found in (Allen and Hui, 1980, Appendix; the authors refer to this as Brooks' theorem). Incidentally, the first ab initio calculation of temperature-dependent band gaps relied on this approach (King-Smith et al., 1989).
Most commonly, the coefficients ∂ε nk /∂n qν in Eq. (199) are evaluated using frozen-phonon supercell calculations, via the second derivative of the eigenvalue ε nk with respect to collective atomic displacements along the vibrational eigenmodes e καν (q). This approach was employed by Capaz et al. (2005) to study the temperature dependence of the band gaps in carbon nanotubes (within a tight-binding model), and by Han and Bester (2013) to obtain the zero-point renormalization and temperature dependence of the gaps of silicon and diamond quantum dots. Recent examples include works on diamond, silicon, SiC , as well as CsSnI 3 .
Frozen-phonon supercell calculations based on Eq. (199) carry the advantage that the rigid-ion approximation, which is necessary to obtain Eqs. (195)-(198), is no longer required. Therefore this approach is more accurate in principle. In practice, however, the calculations are challenging as they require large supercells, and the derivatives must be evaluated for every vibrational mode of the supercell. Several computational strategies were developed to tackle this challenge. Patrick and Giustino (2013) proposed to perform the averages leading to Eq. (199) via importance-sampling Monte Carlo integration. Monserrat (2016b) described a constrained Monte Carlo scheme which improves the variance of the Monte Carlo estimator. Recently Zacharias and Giustino (2016) showed that it is possible to perform these calculations more efficiently by replacing the stochastic sampling by a suitable choice of an 'optimum' configuration; the result becomes exact in the thermodynamic limit of large supercells. In order to reduce the computational cost associated with the use of large supercells, Lloyd- Williams and Monserrat (2015) introduced 'non-diagonal' supercells, which allow one to access phonon wavevectors belonging to a uniform grid of N p points using supercells containing only N 1/3 p unit cells.
The merit of these non-perturbative approaches is that they treat explicitly the nuclear wavefunctions, and enable exploring effects which go beyond the Allen-Heine theory. For example Monserrat et al. (2013Monserrat et al. ( , 2015 and Engel et al. (2015) were able to investigate effects beyond the harmonic approximations in several systems, such as LiH, LiD, high-pressure He, molecular crystals of CH 4 , NH 3 , H 2 O, HF, as well as Ice. In all these cases the authors found large zero-point effects on the band gaps.
Finally, we mention that the calculation of electronic properties at finite temperature via the Allen-Heine the-ory and its variants is closely related to what one would obtain using path-integral molecular dynamics simulations (Della Sala et al., 2004;Ramírez et al., 2006Ramírez et al., , 2008, or even classical molecular dynamics simulations at high enough temperatures (Franceschetti, 2007).

B. Phonon-assisted optical absorption
In addition to modifying the electron energy levels in solids, the electron-phonon interaction plays an important role in the optical properties of semiconductors and insulators, as it is responsible for phonon-assisted optical transitions. Phonon-assisted processes could be analyzed by considering the many-body electronic screening function e (12) in Eq. (95), by using the electron Green's function G dressed by the electron-phonon self-energy Σ ep as in Eq. (151). Since this would require us a lengthy detour, here we simply reproduce the standard result of secondorder time-dependent perturbation theory (Bassani and Parravicini, 1975;Ridley, 1993): In this expression α(ω) is the absorption coefficient for visible light, e is the photon polarization, v mn are the matrix elements of the electron velocity operator, and n r (ω) is the real part of the refractive index. The two denominators in the second line corresponds to indirect processes whereby a photon is absorbed and a phonon is absorbed or emitted (left), and processes whereby a phonon is absorbed or emitted, and subsequently a photon is absorbed (right). The above expression relies on the electric dipole approximation and is therefore valid for photon energies up to a few electronvolts. The theory leading to Eq. (200) was originally developed by Hall, Bardeen, and Blatt (1954). The first ab initio calculation employing Eq. (200) was reported by Noffsinger et al. (2012) for the prototypical case of silicon. The authors employed DFT for computing phonons and electron-phonon matrix elements, and the GW method for the quasiparticle band structures. The sampling of the Brillouin zone was achieved by means of the interpolation strategy described in Sec. VI. Figure 14 shows that the calculations by Noffsinger et al. are in very good agreement with experiment throughout the energy range of indirect absorption.
Further work along similar lines was reported by Kioupakis et al. (2010), who calculated the indirect optical absorption by free carriers in GaN;and Peelaers et al. (2015), who studied the indirect absorption by free carriers in the transparent conducting oxide SnO 2 . Recently, FIG. 14 (Color online) Phonon-assisted optical absorption in silicon: comparison between first-principles calculations (solid lines, orange) and experiment (circles, blue). The calculations were performed using the theory of Hall et al. (1954), as given by Eq. (200). Spectra calculated at different temperatures were shifted horizontally so as to match the experimental onsets. Reproduced with permission from (Noffsinger et al., 2012), copyright (2012) by the American Physical Society.
the ab initio theory of phonon-assisted absorption was also extended to the case of indirect Auger recombination by Kioupakis et al. (2015).
One limitation of the theory by Hall et al. is that the indirect absorption onset is independent of temperature. This is seen by noting that the Dirac delta functions in Eq. (200) contain the band structure energies at clamped nuclei. The generalization to incorporate temperaturedependent band structures as discussed in Sec. IX.A.1 is nontrivial. Patrick and Giustino (2014) and Zacharias and Giustino (2016) showed that the electron-phonon renormalization of the band structure modifies the energies of real transitions but leaves unchanged the energies of virtual transitions; in other words the Allen-Heine renormalization should be incorporated only in the Dirac delta functions and in the first denominator in Eq. (200). In order to avoid these complications at once, Zacharias et al. (2015) developed an alternative approach which relies on the 'semiclassical' approximation of Williams (1951) and Lax (1952). In this approximation, the initial states in the optical transitions are described quantummechanically, and the final states are replaced by a quasiclassical continuum. In the formulation of Zacharias et al. the imaginary part of the temperature-dependent dielectric function takes the form: where 2 (ω) denotes the imaginary part of the dielectric function at clamped nuclei, and the expectation values have the same meaning as in Eq. (199). E {nqν } is the energy of the quantum nuclear state specified by the quantum number {n qν } and Z is the canonical partition function. Zacharias et al. demonstrated that this ap-proach provides an adiabatic approximation to Eq. (200), and seamlessly includes the temperature dependence of the electronic structure within the Allen-Heine theory.
Using techniques similar to those of Sec. IX.A.2, the authors calculated the indirect optical absorption lineshape of silicon at various temperatures and obtained very good agreement with experiment. These results were recently extended to the temperature-dependent optical spectra of diamond and gallium arsenide by Zacharias and Giustino (2016).

A. Electron linewidths and lifetimes
In Sec. V.B.4 we have seen how the Fan-Migdal selfenergy can be used in order to evaluate the quasiparticle lifetimes (or equivalently linewidths) resulting from the electron-phonon interaction. The first ab initio calculations of this kind were reported by Eiguren et al. (2003aEiguren et al. ( , 2002 in the study of the decay of metal surface states, and by Sklyadneva et al. (2005Sklyadneva et al. ( , 2006 and Leonardo et al. (2007) in the study of electron lifetimes of elemental metals. Some of these calculations and the underlying approximations were reviewed in Sec. VIII. Calculations of quasiparticle linewidths were also employed to study the temperature-dependent broadening of the optical spectra in semiconductors. For example Giustino et al. (2010) and Poncé et al. (2015) investigated the broadening of the direct absorption edge of diamond and silicon, respectively. In both cases good agreement with experiment was obtained. More recently, the same approach was employed to study the broadening of photoluminescence peaks in lead-iodide perovskites (Wright et al., 2016). In this case, it was found that the standard Fermi golden rule expression, Eq. (170), significantly overestimates the experimental data. The agreement with experiment is restored by taking into account the quasiparticle renormalization Z nk ; see discussion following Eq. (174).
While these works were primarily concerned with the broadening of the spectral lines in photoemission or optical experiments, Eq. (170) can also be used to study carrier lifetimes in time-resolved experiments. The first ab initio study in this direction was reported from Sjakste et al. (2007), who investigated the thermalization of hot electrons in GaAs and the exciton lifetimes in GaP. In the case of GaAs, Sjakste et al. found thermalization rates in quantitative agreement with time-resolved luminescence and transient optical absorption measurements. Work along similar lines was also reported for the intervalley scattering times in Ge . Recently, the thermalization rates of hot electrons in GaAs were revisited by Bernardi et al. (2015). The authors employed Eq. (170) and the Wannier interpolation technique described in Sec. VI in order to finely sample the electron-phonon scattering processes near the bottom of the conduction band, see Fig. 15. Based on these calculations they were able to interpret transient absorption measurements in terms of the carrier lifetimes within each valley. Another interesting application of Eq. (170) was reported by Bernardi et al. (2014), who investigated the rate of hot carrier thermalization in silicon, within the context of photovoltaics applications.
Very recently Sangalli and Marini (2015a,b) employed the lifetimes calculated using Eq. (170) in order to study carrier dynamics in silicon in real time. Strictly speaking, these developments lie outside of the scope of equilibrium Green's functions discussed in Sec. IV, and require concepts based on non-equilibrium Green's functions (Kadanoff and Baym, 1962). However, the basic ingredients of the electron-phonon calculations remain unchanged.
In all calculations discussed in this section, the electron-phonon matrix elements were obtained within DFT. However, in order to accurately describe electronphonon scattering near band extrema Bernardi et al. (2014Bernardi et al. ( , 2015 and Wright et al. (2016) employed GW quasiparticle band structures. This is important in order to obtain accurate band effective masses, which affect the carrier lifetimes via the density of states.

B. Phonon-limited mobility
The carrier lifetimes τ nk of Eq. (170) are also useful in the calculation of electrical mobility, conductivity, and resistivity, within the context of the semiclassical model of electron dynamics in solids (Ashcroft and Mermin, 1976). In the semiclassical model, one describes the electronic response to an external perturbation by taking the fermionic occupations f nk to represent the probability density function in the phase space defined by the unperturbed band structure. The probabilities f nk are then determined using a standard Boltzmann equation. A comprehensive discussion of these methods can be found in the classic book of Ziman (1960).
Here we only touch upon the key result which is needed in ab initio calculations of electrical conductivity. In the semiclassical model, the electrical current is calculated as J = −2e Ω −1 BZ n dk v nk f nk . In the absence of external fields and thermal gradients, the occupations f nk reduce to the standard Fermi-Dirac occupations at equilibrium, f 0 nk , and the current vanishes identically. Upon introducing an electric field E, the electrons respond by adjusting their occupations. In this model it is assumed that the variation f nk − f 0 nk is so small that the electronic density is essentially the same as in the unperturbed system. The modified occupations can be calculated using the linearized Boltzmann transport equation (Ziman, 1960): The left-hand side of Eq. (202) represents the collisionless term of the Boltzmann equation, that is the change in occupations due to the particle drift under the electric field. The right-hand side represents the change of occupations resulting from electrons scattered in or out of the state |nk by phonon emission or absorption. The rates in Eq. (203) are simply derived from Fermi's golden rule (Grimvall, 1981). By solving Eq. (202) self-consistently for all f nk it is possible to calculate the current, and from there the conductivity. The connection with the carrier lifetimes τ nk of Eq. (170) is obtained within the so-called 'energy relaxation time approximation'. In this approximation the incoming electrons are neglected in Eq. (202), that is the last term (f mk+q − f 0 mk+q ) in the second line is ignored. As a result the entire right-hand side of the equation simplifies to −(f nk − f 0 nk )/τ nk . The direct solution of Eq. (202) is computationally challenging, and fully ab initio calculations were reported only very recently by Li (2015) for Si, MoS 2 , and Al, and by Fiorentini and Bonini (2016) for n-doped Si. Figure 16 shows that the mobility of n-type silicon calculated by Li is in good agreement with experiment. The theory overestimates the measured values to some extent, and this might have to do with the limitations of the DFT matrix elements (see Sec. XII). In addition to the carrier mobility, Fiorentini and Bonini employed the ab initio Boltzmann formalism to calculate thermoelectric properties, such as the Lorenz number and the Seebeck coefficient.
Ab initio calculations of the resistivity of metals are less challenging than for semiconductors, and started appearing already with the work of Bauer et al. (1998). Most calculations on metals are based on Ziman's resistivity formula, see Grimvall (1981, p. 210). An interesting recent example can be found in the work by Xu and Verstraete (2014) on the transport coefficients of lithium. We also highlight related work on phonon-limited transport in organic crystals Bobbert, 2004a, 2004b;Ortmann et al., 2009;Vukmirović et al., 2012).

XI. PHONON-MEDIATED SUPERCONDUCTORS
The last application of ab initio calculations of EPIs that we will consider is the study of phonon-mediated superconductivity (Schrieffer, 1983). This research field is so vast that any attempt at covering it in a few pages would not make justice to the subject. For this reason, it was decided to limit the discussion to those novel theoretical and methodological developments which are aiming at fully predictive calculations, namely the 'anisotropic Migdal-Eliashberg theory' (Sec. XI.B), and the 'density functional theory for superconductors' (Sec. XI.C). For completeness, in Sec. XI.A we also summarize the most popular equations employed in the study of phononmediated superconductors. All calculations described in this section were performed at the DFT level.

A. McMillan/Allen-Dynes formula
Most ab initio calculations on phonon-mediated superconductors rely on a semi-empirical expression for the critical temperature, first introduced by McMillan (1968) and then refined by Allen and Dynes (1975): .
Here T c is the superconducting critical temperature, ω log is a 'logarithmic' average of the phonon frequencies (Allen and Dynes, 1975), λ is the electron-phonon 'coupling strength', and µ * is a parameter describing the Coulomb repulsion. The functional form of Eq. (204) (McMillan, 1968;Allen and Dynes, 1975;Grimvall, 1981;Allen and Mitrovic, 1982): where N F is the density of states at the Fermi level and the matrix elements g mnν (k, q) are the same as in Eq. (38). The remaining parameter µ * (Morel and Anderson, 1962) is obtained as 1/µ * = 1/µ + log(ω p /ω ph ), where ω p is the characteristic plasma energy of the system, ω ph the largest phonon energy, and µ is the average electron-electron Coulomb repulsion across the Fermi surface. More specifically: µ = N F V nk,n k FS , where · · · FS denotes a double average over the Fermi surface, and V nk,n k = k n , −k n |W |kn, −kn , with W being the screened Coulomb interaction of Sec. IV.B.2 (Lee et al., 1995;Lee and Chang, 1996).
The coupling strength λ is related to the mass enhancement parameter λ nk discussed in Sec. VIII..1. The main difference between λ and λ nk is that the former represents an average over the Fermi surface, while the latter refers to the Fermi velocity renormalization of a specific electron band. While these quantities are related, they do not coincide and hence cannot be used interchangeably.
Equations (204)-(207) involve a number of approximations. For example, it is assumed that the superconductor is isotropic and exhibits a single superconducting gap. Furthermore, almost invariably the effective Coulomb potential µ * is treated as an adjustable parameter, on the grounds that it should be in the range µ * = 0.1-0.2. This procedure introduces a large uncertainty in the determination of T c , especially at moderate coupling strengths.

B. Anisotropic Migdal-Eliashberg theory
A first-principles approach to the calculation of the superconducting critical temperature is provided by the anisotropic Migdal-Eliashberg theory (Migdal, 1958;Eliashberg, 1960). This is a field-theoretic approach to the superconducting pairing, formulated in the language of finite-temperature Green's functions. At variance with the Hedin-Baym equations of Table I, the Migdal-Eliashberg theory is best developed within the Nambu-Gor'kov formalism (Gor'kov, 1958;Nambu, 1960), which enables describing the propagation of electron quasiparticles and of superconducting Cooper pairs on the same footing (Scalapino, 1969;Schrieffer, 1983). A comprehensive presentation of the Migdal-Eliashberg theory is provided by Allen and Mitrovic (1982). Their article served as the starting point of current first-principles implementations of the theory.
In the Migdal-Eliashberg theory, one solves the two coupled equations: where k stands for Ω −1 BZ dk . In these equations, T is the absolute temperature, Z nk (iω j ) is the quasiparticle renormalization function, and is analogous to Z nk in Eq. (189). ∆ nk (iω j ) the superconducting gap function. The functions Z nk (iω j ) and ∆ nk (iω j ) are determined along the imaginary frequency axis, at the fermion Matsubara frequencies iω j = i(2j + 1)πk B T / with j an integer. The anisotropic and frequency-dependent generalization of Eq. (206) to be used in the Migdal-Eliashberg equations is: with q = k −k. Equations (208)-(209) are to be solved self-consistently for each temperature T . The superconducting critical temperature is then obtained as the highest temperature for which a nontrivial solution is obtained, that is a solution with ∆ nk (iω j ) = 0. From the superconducting gap along the imaginary axis it is then possible to obtain the gap at real frequencies by analytic continuation (Marsiglio et al., 1988), and from there one can compute various thermodynamic functions. The first ab initio implementation of the anisotropic Migdal-Eliashberg theory was reported by Choi et al. (2002a,b) and Choi et al. (2003) in a study of the superconducting properties of MgB 2 . The authors succeeded to explain the anomalous heat capacity of MgB 2 in terms of two distinct superconducting gaps, and obtained a T c in good agreement with experiment. These calculations were later extended to MgB 2 under pressure (Choi et al., 2009a) and other hypothetical borides (Choi et al., 2009b). Margine and Giustino (2013) demonstrated an implementation of the Migdal-Eliashberg theory based on the Wannier interpolation scheme of Sec. VI, and reported applications to Pb and MgB 2 . The superconducting gap and superconducting density of states of MgB 2 calculated by Margine and Giustino are shown in Fig. 17. In all these calculations, the Coulomb repulsion was described empirically via µ * , and this partly accounts for the slight discrepancy between the calculated T c of 50 K and the experimental T c of 39 K (Nagamatsu et al., 2001). Additional calculations based on the anisotropic Migdal-Eliashberg theory include a study of doped graphene (Margine and Giustino, 2014), as well as investigations of Li-decorated monolayer graphene (Zheng and Margine, 2016) and Ca-intercalated bilayer graphene . In this latter work the authors incorporated Coulomb interactions from first principles, after calculating µ * via the screened Coulomb interaction in the random-phase approximation. The calculated T c = 7-8 K was in reasonable agreement with the experimental value of 4 K (Ichinokura et al., 2016). The Migdal-Eliashberg theory has also been extended to describe the superconducting state as a function of applied magnetic field; a complete ab initio implementation was successfully demonstrated with an application to MgB 2 (Aperis et al., 2015). Very recently Sano et al. (2016) performed ab initio Migdal-Eliashberg calculations including retardation effects on high-pressure sulfur hydrides, obtaining good agreement with experiment. Interestingly in this work the authors also checked the effect of the zero-point renormalization of the electron bands within the Allen-Heine theory, and found that it accounts for a change in T c of up to 20 K.

C. Density functional theory for superconductors
Another promising ab initio approach to the calculation of the superconducting critical temperature is the density functional theory for superconductors Marques et al., 2005). The starting point of this approach is a generalization of the Hohenberg-Kohn theorem (Hohenberg and Kohn, 1964) to a system described by three densities: the electron density in the normal state, the density of superconducting pairs, and the nuclear density. Based on this premise, Lüders et al. mapped the fully-interacting system into an equivalent Kohn-Sham system (Kohn and Sham, 1965) of non-interacting nuclei and non-interacting, yet superconducting, electrons. The resulting Kohn-Sham equations for the electrons take the form of Bogoliubov-de Gennes equations (Bogoliubov, 1958), whereby electrons are paired by an effective gap function ∆(r, r ).
In its simplest formulation, the density functional theory for superconductors determines the expectation value of the pairing field over Kohn-Sham eigenstates, ∆ nk = u nk (r)|∆(r, r )|u nk (r ) , using the following gap equation: where E 2 nk = ε 2 nk + |∆ nk | 2 . In this expression, the kernel K contains information about the phonon-mediated pairing interaction and the Coulomb repulsion between electrons, K = K ep + K ee , and Z contains information about the electron-phonon interaction. More specifically, K ep and Z are evaluated starting from the electronphonon matrix elements g mnν (k, q) and the DFT electron band structure and phonon dispersions, as in the Migdal-Eliashberg theory. K ee is approximated using the screened Coulomb interaction V nk,n k introduced below Eq. (207). Complete expressions for Z and K can be found in .
Equation (211) is reminiscent of the gap equation in the Bardeen-Cooper-Schrieffer (BCS) theory (Schrieffer, 1983), with the difference that the ab initio kernel K replaces the model interaction of the BCS theory, and the function Z introduces quasiparticle renormalization as in the Migdal-Eliashberg theory, see Eq. (208). At variance with the Migdal-Eliashberg theory, the gap function in the density functional theory for superconductors does not carry an explicity frequency dependence. Nevertheless, retardation effects are fully included through the dependence of the kernels Z and K on the electron bands and the phonon dispersions. An important advantage of this theory is that the Coulomb potential µ * is not required, since the electron-electron repulsion is seamlessly taken into account by means of the kernel K ee .
Further developments of the superconducting density functional theory include the study of non-phononic pairing mechanisms, such as plasmon-assisted superconductivity (Akashi and Arita, 2013a), and the extension to magnetic systems (Linscheid et al., 2015a,c).

XII. ELECTRON-PHONON INTERACTIONS BEYOND THE LOCAL DENSITY APPROXIMATION TO DFT
The calculations of electron-phonon interactions reviewed in Sec. VII-XI have in common the fact that most investigations used the local density approximation to DFT or a generalized gradient approximation (GGA) such as the PBE functional (Perdew et al., 1996a). Although the LDA and the GGA do represent the workhorse of electron-phonon calculations from first principles, there is growing evidence that these choices can lead to an underestimation of the electron-phonon coupling strength. At a conceptual level we can understand this point by rewriting the electron-phonon matrix element after combining Eqs. (38) and (142), (143): g mnν (k, q) = u mk+q | dr −1 e (r, r , ω)∆ qν v en (r )|u nk uc .
(212) In DFT the many-body dielectric matrix e appearing in this expression is replaced by the RPA+xc screening Hxc from Eq. (54). Given the DFT band gap problem, we expect Hxc to overestimate the screening, thereby leading to matrix elements g mnν (k, q) which are underestimated to some extent.
Several groups investigated this point on quantitative grounds. Zhang et al. (2007) studied the electronphonon coupling in a model copper oxide superconductor, CaCuO 2 . By calculating the vibrational frequencies of the half-breathing Cu-O stretching mode, the authors established that the local spin-density approximation (LSDA) yields phonons which are too soft (65.3 meV) as compared to experiment (80.1 meV). In contrast, the introduction of Hubbard corrections in a LSDA+U scheme restored agreement with experiment (80.9 meV). Since the electron-phonon matrix elements are connected to the phonon frequencies via the phonon self-energy, Eq. (145), a corresponding underestimation of the matrix elements can be expected. These results were supported by the work of Floris et al. (2011), who developed DFPT within LSDA+U , and applied their formalism to the phonon dispersions of antiferromagnetic MnO and NiO. Here the authors found that the DFT underestimates measured LO energies by as much as 15 meV in MnO, while the use of LSDA+U leads to good agreement with experiment. Related work was reported by Hong et al. (2012), who investigated the multiferroic perovskites CaMnO 3 , SrMnO 3 , BaMnO 3 , LaCrO 3 , LaFeO 3 , and the double perovskite La 2 CrFeO 6 . Here the authors calculated the variation of the vibrational frequencies between the ferromagnetic and the antiferromagnetic phases of these compounds as a function of the Hubbard U parameter, and compared DFT+U calculations with hybrid-functional calculations. Lazzeri et al. (2008) investigated the effect of quasiparticle GW corrections on the electron-phonon coupling of graphene and graphite, for the A 1 phonon at K and the E 2g phonon at Γ. They evaluated the intraband electronphonon matrix elements using a frozen-phonon approach, noting that g nnν (k, q = 0) represents precisely the shift of the Kohn-Sham energy ε nk upon displacing the atoms according to the ν-th phonon eigenmode at q = 0. Lazzeri et al. found that the matrix elements increase by almost 40% from DFT to GW . The GW values led to slopes in the phonon dispersions near K in very good agreement with inelastic X-ray scattering data (Grüneis et al., 2009b). Similar results, albeit less dramatic, were obtained by Grüneis et al. (2009a) for the potassiumintercalated graphite KC 8 .
Laflamme Janssen et al. (2010) studied the electronphonon coupling in the C 60 molecule as a model for superconducting alkali-doped fullerides. They employed the PBE0 hybrid functional (Perdew et al., 1996b) with a fraction of exact exchange α = 30%, and obtained an enhancement of the total coupling strength λ of 42% as compared to PBE. This work was followed up by Faber et al. (2011), who used the GW approximation and obtained a similar enhancement of 48%. We also point out an earlier work by Saito (2002) based on the B3LYP functional, reporting similar results. Yin et al. (2013) investigated the effects of using the GW approximation and the HSE hybrid functional (Heyd et al., 2003) on the electron-phonon coupling in the superconducting bismuthates Ba 1−x K x BiO 3 and chloronitrides β-ZrNCl, as well as MgB 2 . In the case of Ba 1−x K x BiO 3 the authors obtained a three-fold increase in the coupling strength λ from PBE to HSE. This enhancement brought the critical temperature calculated using Eq. (204) to 31 K, very close to the experimental value of 32 K. Similarly, in the case of β-ZrNCl, Yin et al. obtained a 50% increase of λ, bringing the calculated critical temperature, 18 K, close to the experimental value of 16 K. Instead, in the case of MgB 2 , they noticed only a slight increase of the electron-phonon coupling as compared to the standard LDA.
Another application of hybrid functionals to the study of EPIs was reported by Komelj and Krakauer (2015). Here the authors investigated the sensitivity of the superconducting critical temperature of the H 3 S phase of sulfur hydride to the exchange and correlation functional.
They found that the PBE0 functional enhances the critical temperature by up to 25% as compared to PBE, bringing T c from 201-217 K to 253-270 K (the spread in values is related to the choice of the parameter µ * ). Mandal et al. (2014) reported work on the superconductor FeSe based on dynamical mean-field theory (DMFT). In this case DMFT yielded a three-fold enhancement of the coupling strength for selected modes.
As already mentioned in Sec. IX.A.1, Antonius et al. (2014) performed GW calculations of the electronphonon coupling in diamond using a frozen-phonon approach. They found that quasiparticle corrections lead to a uniform enhancement of the electron-phonon matrix elements. The net effect is an increase of the zero-point renormalization of the band gap by 40% as compared to standard LDA calculations. Monserrat (2016a) confirmed this result and found a GW correction of comparable magnitude in the case of silicon. However, Monserrat also found that the GW corrections to the zero-point band gap renormalization of LiF, MgO, and TiO 2 are very small (∼5% of the PBE value), therefore at present it is not possible to draw general conclusions.
Finally, we mention that Faber et al. (2015) examined possible strategies for systematically incorporating GW corrections in electron-phonon calculations. By using diamond, graphene, and C 60 as test cases, the authors showed that a 'constant screening' approximation is able to reproduce complete GW results with an error below 10% at reduced computational cost. This approximation amounts to evaluating the variation of the Green's function G in a frozen-phonon calculation, while retaining the screened Coulomb interaction W of the unperturbed ground state.
All these recent developments point to the need of moving beyond local exchange and correlation functionals in the study of electron-phonon interactions from first principles. In the future, it will be important to devise accurate computational methods for calculating not only the intraband electron-phonon matrix elements (as in the frozen-phonon method) but also matrix elements between all states and for scattering across the entire Brillouin zone.
For the sake of completeness we emphasize that the underestimation of the EPI matrix elements by semilocal DFT functionals does not propagate in the same way into different materials properties. This is readily understood by examining two fundamental quantities, the Allen-Heine renormalization of electron bands, Eq. (167), and the adiabatic phonon frequencies, as obtained from Eqs. (126) and (127). In the former case the electronic screening enters as −2 ∞ ; in the latter case the screening contributes through a term which scales with −1/2 ∞ . As a result, in the hypothetical case of a semiconductor for which DFT underestimated the electronic permittivity by 20%, we would have an error of ∼40% in the energy renormalization, and of ∼10% in the phonon frequen-cies. This example is an oversimplification of the problem, but it shows that different properties relating to the EPI could be affected to a very different degree by the inherent limitations of DFT functionals.

XIII. CONCLUSIONS
The study of electron-phonon interactions has a long and distinguished history, but it is only during the past two decades that quantitative and predictive calculations have become possible. First-principles calculations of electron-phonon couplings are finding an unprecedented variety of applications in many areas of condensed matter and materials physics, from spectroscopy to transport, from metals to semiconductors and superconductors. In this article we discussed the standard DFT formalism for performing calculations of electron-phonon interactions, we showed how most equations can be derived from a field-theoretic framework using a few well-defined approximations, and we reviewed recent applications of the theory to many materials of current interest.
As calculation methods improve relentlessly and quantitative comparisons between theory and experiment become increasingly refined, new and more complex questions arise. Much is still left to do, both in the fundamental theory of electron-phonon interactions, and in the development of more accurate and more efficient computational methods.
For one, we are still using theories where the coupling matrix elements are calculated using the adiabatic local density approximation to DFT. The need for moving beyond standard DFT and beyond the adiabatic approximation can hardly be overemphasized. Progress is being made on the incorporation of nonlocal corrections into electron-phonon matrix elements, for example using hybrid functionals or GW techniques, but very little is known about retardation effects. It is expected that such effects may be important in the study of heavily doped oxides and semiconductors, both in their normal and superconducting states (Mahan, 1993, Sec. 6.3.A), but ab initio investigations are currently missing. This is truly uncharted territory and a systematic effort in this direction is warranted.
In this article we emphasized that it is possible to formulate a compact, unified theory of electron-phonon interactions starting from a fully ab initio field-theoretic approach. The only assumption which is absolutely crucial to the theory is the harmonic approximation. Abandoning the harmonic approximation leads to the appearance of several new terms in the equations, and the resulting formalism becomes considerably more complex than in Table I. Despite these difficulties, given the importance of anharmonic effects in many systems of current interest, extending the theory to the case of anharmonic phonons and multi-phonon interactions consti-tutes a pressing challenge. Ab initio investigations of anharmonic effects on the temperature dependence of band gaps have recently been reported (Monserrat et al., 2013;Antonius et al., 2015). Since these studies rely on nonperturbative adiabatic calculations in supercells, it would be highly desirable to establish a clear formal connection of these methods with the rigorous field-theoretic approach of Sec. IV. Along the same line, it would be important to clarify the relation between many-body approaches, adiabatic supercells calculations, and more traditional classical or path-integral molecular dynamics simulations.
The study of electron-phonon interactions has long been dominated by Fröhlich-like Hamiltonians, whereby the electron-phonon coupling is retained only to linear order in the atomic displacements. This is the case for all the model Hamiltonians mentioned in Sec. V.B.6. It is now clear that quadratic couplings, leading to the Debye-Waller contributions in the optical spectra of semiconductors, are by no means negligible and should be investigated more systematically. For example, in the current literature it is invariably assumed that Debye-Waller contributions are negligible in metals near the Fermi surface; while this is probably the case for the simplest elemental metals, what happens in the case of multiple Fermisurface sheets is far from clear, and should be tested by direct calculations.
The identification of the correct matrix elements to be calculated is not always a trivial task, as it was discussed for the case of the non-adiabatic phonon self-energy. In the future it will be important to pay attention to these aspects, especially in view of detailed comparison with experiment. For now, the issue on whether the phonon self-energy arising from EPIs should be calculated using bare or screened EPI matrix elements (Sec. VII) is to be considered an open question, and calls for further investigation.
The theory and applications reviewed in this article focused on non-magnetic systems. The rationale for this choice is that a complete many-body theory of electronphonon interactions for magnetic systems is not available yet. Recent investigations of spin-phonon couplings were conducted by assuming that the spin and the vibrational degrees of freedom can be decoupled, as in the Born-Oppenheimer approximation. Under this assumption it is possible to investigate how the spin configuration responds to a frozen phonon, or alternatively how the vibrational frequencies depend on the spin configuration (see for example Chan et al., 2007;Łażewski et al., 2010;Lee and Rabe, 2011;Cao et al., 2015). In all these cases it would be desirable to employ a more rigorous manybody theory of spin-phonon interactions. The Hedin-Baym equations discussed in Sec. IV maintain their validity in the case of spin-polarized systems, provided collinear spins are assumed. In more general situations, where it is important to consider noncollinear spins, ex-ternal magnetic perturbations, or spin-dependent interactions such as spin-orbit and Rashba-Dresselhaus couplings, it becomes necessary to generalize the equations in Table I. Although such a generalization has not been reported yet, the work of Aryasetiawan and Biermann (2008) constitutes a promising starting point. In that work the Schwinger functional derivative technique (see Sec. IV.B.1) was used to extend Hedin's equations at clamped nuclei to systems containing spin-dependent interactions. Generalizing Aryasetiawan and Biermann's work to incorporate nuclear vibrations will be important for the study of electron-phonon interactions in many systems of current interest, from multifunctional materials to topological quantum matter.
At this time it is not possible to predict how this fastmoving field will evolve over the years to come. However, the impressive progress made during the past decade gives us confidence that this interesting research area will continue to thrive, and will keep surprising us with fascinating challenges and exciting new opportunities. carries some redundancy. Indeed by combining Eqs. (18) and (B1) it is seen that: The inverse relation of Eq. (B1) is: The right-hand side is real-valued after Eqs. (18) and (B2). In preparation for the transition to a quantum description of lattice vibrations, it is useful to identify 3M N p independent normal coordinates. This can be done by partitioning the grid of q-vectors in three sets. We call A the set of vectors which are invariant under inversion, that is −q + G = q for some reciprocal lattice vector G (including |G| = 0). The center of the Brillouin zone and the centers of its faces belong to this set. The remaining vectors can be separated further in B and C, in such a way that all the vectors in C are obtained from those in B by inversion (modulo a reciprocal lattice vector). After defining z qν = x qν + iy qν , Eq. (B3) can be rewritten as: q∈A,ν e κα,ν (q) x qν + 2Re q∈B,ν e iq·Rp e κα,ν (q)(x qν + iy qν ) .(B4) The q-vectors of the set C have been grouped together with those in B by taking the real part in the second line. It can be verified that in this expression there are exactly 3M N p real coordinates, therefore we can choose the x qν for q in A and the pairs x qν , y qν for q in B as the independent variables. These variables are referred to as 'real normal coordinates' (Brüesch, 1982). Using Eqs. (12)-(18), (A1), and (B4) the nuclear Hamiltonian can be written in terms of 3M N p independent harmonic oscillators in the real normal coordinates: H p = 1 2 q∈B,ν ω qν (−∂ 2 /∂x 2 qν − ∂ 2 /∂ỹ 2 qν +x 2 qν +ỹ 2 qν ) + 1 2 q∈A,ν ω qν (−∂ 2 /∂x 2 qν +x 2 qν ), where for ease of notation we performed the scaling: x qν = x qν /2 l qν for q in A, (B6) x qν = x qν /l qν ,ỹ qν = y qν /l qν for q in B, (B7) with l qν being the zero-point displacement amplitude of Eq. (21). In the case of |q| = 0 there are three normal modes for which ω qν = 0, and the corresponding potential termsx 2 qν must be removed from Eq. (B5). The eigenstates of Eq. (B5) are found by introducing the real ladder operators for each normal coordinate (Cohen-Tannoudji et al., 1977): and similarly forâ qν,y . With these definitions Eq. (B5) becomes: H p = q∈B,ν ω qν â † qν,xâqν,x +â † qν,yâqν,y + 1 + q∈A,ν ω qν â † qν,xâqν,x + 1/2 .
The eigenstates of this Hamiltonian are products of simple harmonic oscillators (Merzbacher, 1998), and the ground state is: The eigenstates ofĤ p can be generated by applyinĝ a † qν,x andâ † qν,y to the ground state χ 0 . However this approach is not entirely satisfactory, since we cannot assign separate quantum numbers to modes with wavevectors q or −q. In order to avoid this inconvenience we observe that, for each normal mode, the first line of Eq. (B5) defines an effective isotropic two-dimensional harmonic oscillator. The degenerate eigenstates of this oscillator can be combined to form eigenstates of the angular momentum; this leads to right and left circular quanta with the same energy and definite angular momentum (Cohen-Tannoudji et al., 1977). This analogy motivates the consideration of the following linear combinations, for q in B: a + qν = (â qν,x + iâ qν,y )/ √ 2, a − qν = (â qν,x − iâ qν,y )/ √ 2.
Since bothâ qν,x andâ qν,y lower the energy of an eigenstate by the same quantum of energy ω qν , the resulting states are degenerate and their linear combinations are also eigenstates for the same eigenvalue. As a consequence we can generate all the eigenstates of the Hamil-tonianĤ p by acting on the ground state χ 0 with the creation operatorsâ +, † qν andâ −, † qν . In this reasoning the wavevectors q belong to B; if we now consider Eqs. (B2), (B11), and (B12) we see that formally we also havê a − qν =â + −qν . Therefore it is natural to associateâ − qν to phonons propagating along the direction −q. These observations suggest replacing the real ladder operators of Eq. (B8) by the complex ladder operatorŝ a + qν andâ − −qν for q in B and C, respectively. In the case of q in A we keep the real operatorsâ qν,x . These definitions can be turned into the compact expressions: a qν =â qν,x for q in A, (B13) a qν = (â qν,x + iâ qν,y )/ √ 2 for q in B, C. (B14) Using these operators the nuclear Hamiltonian of Eq. (B5) takes the well-known form given by Eq. (22). Any eigenstate ofĤ p can now be generated as qν (n qν !) − 1 2 (â † qν ) nqν χ 0 . In this form we see that it is possible to assign independently a number of phonons n qν to each wavevector q and each mode ν. Using Eqs. (B6)-(B8) and (B13)-(B14) we also have the basic identity: z qν = l qν (â qν +â † −qν ).