{'fraction_non_alphanumeric': 0.03792209636757448, 'fraction_numerical': 0.011131626602279586, 'mean_word_length': 4.6136137133745265, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Eliashberg Theory: a short review

Eliashberg theory is a theory of superconductivity that describes the role of phonons in providing the attractive interaction between two electrons. Phonon dynamics are taken into account, thus giving rise to retardation effects that impact the electrons, in the form of a frequency-dependent electron self-energy. In the superconducting state, this means that the order parameter, generally considered to be a static quantity in the Bardeen-Cooper-Schrieffer (BCS) theory, also becomes frequency-dependent. Here we review the finite temperature formulation of Eliashberg theory, both on the imaginary and real frequency axis, and briefly display some examples of the consequences of a dynamical, as opposed to static, interaction. Along the way we point out where further work is required, concerning the validity of some of the assumptions used.


I. INTRODUCTION
Superconductivity is a remarkable phenomenon, not least because it represents a manifestation of the quantum world on a macroscopic scale. It is spectacularly demonstrated with levitating train sets, 1 and indeed this property and many others of superconductors are slowly being utilized in everyday applications. 2 However, the established practice of incorporating superconductors into the real world should not be taken as an indication that "the last nail in the coffin [of superconductivity]" 3 has been achieved. On the contrary, in the intervening halfcentury since this quote was written, many new superconductors have been discovered, and we have reached a point where it is clear that a deep lack of understanding 4 of superconductivity currently exists. Reference [5] compiles a series of articles reviewing the various "families" or classes of superconductors, where one can readily see common and different characteristics. At the moment many of these classes require a class-specific mechanism for superconductivity, a clearly untenable situation, in my opinion. Further classes have been discovered or expanded upon since, such as nickelates, 7 and the hydrides under pressure, 6 for example.
Our "deep lack of understanding" should not be taken to indicate that theoretical contributions have not been forthcoming. In fact there have been remarkable contributions to key theoretical ideas in physics that stem from research in superconductivity, starting with London theory 8 and Ginzburg-Landau 9 theory, through to BCS theory. 10 When Gor'kov 11 recast the BCS theory of superconductivity in the language of Green functions, then the stage was set for Eliashberg 12,13 to formulate the theory that bears his name. It is fitting that we honour the lasting impact of his work with this brief review, on the occasion of his 90th birthday, and the 60th anniversary of the publication of two papers that paved the path for considerable future quantitative work in superconductivity. Based on an index I am fond of using for famous people, his name appears in titles of papers 248 times, and in abstracts and keywords 1439 times. 16 Before proceeding further, we wish to make some remarks about the nature of this review. It will necessarily repeat material from previous reviews, which we catalogue as follows. Scalapino 17 and McMillan and Rowell 18 perhaps gave one of the first comprehensive reviews of both calculations and experiments that provide remarkable evidence for the validity of Eliashberg theory for various superconductors. These reviews were provided in the comprehensive monogram by Parks; 19 the reader should refer to this monogram and the references therein, as we cannot possibly properly reference all the primary literature sources before "Parks", as this would consume too many pages here. The author list in Parks is a who's who of experts in superconductivity, with two notable exceptions, John Bardeen and Gerasim (Sima) Eliashberg.
A subsequent very influential review was that of Allen and Mitrović, 20 where mostly superconducting T c was discussed. These authors highlighted the expediency of doing many calculations on the imaginary frequency axis, a possibility first noted in Ref. [21] and utilized to great advantage in subsequent years. [22][23][24][25] A few years later Rainer wrote a "state-of-the-union" address 26 on first principles calculations of superconducting T c in which a challenge was issued to both band structure and many-body theorists. For the former, the missing ingredient was a complete (italics are mine) calculation of the electron-phonon coupling. These were first calculated in the 1960's (e.g. Ref. [27]) but have experienced vast improvement over the past 50 years, through the adoption and improvement of Density Functional Theory methods, plus the increased computational ability achieved in the intervening decades. Excellent summaries of this progress is provided in Refs. [28,29,30], where two alternative procedures are described. The first follows the original route of determining the electronphonon interaction and including this as input to the Eliashberg equations, while the second aims to treat both the electron-phonon-induced electron-electron interaction and the direct Coulomb interaction on an equal footing. The result is advertised to meet Rainer's challenge and calculate T c and other superconducting properties without any experimental input. It is noteworthy that in the second formulation (see also Refs. [31] and [32]) the equations are BCS-like and do not depend on frequency, but only momentum. In the usual formula-arXiv:1911.05065v1 [cond-mat.supr-con] 12 Nov 2019 tion, the calculation of µ * , the effective direct electronelectron Coulomb interaction, is often simply assigned a (small) numerical value, and therefore is treated phenomenologically as a fitting parameter. The challenge to many-body theorists remains, as more superconductors have been discovered that seem to extend beyond the weak coupling regime, and likely require descriptions beyond BCS and Eliashberg theory.
In 1990 the review by Carbotte 33 provided a comprehensive update for a number of thermodynamic properties of various superconductors known at the time, including the high temperature cuprate materials. Partly for this reason his review is titled "Properties of bosonexchange superconductors," since there was a feeling at the time (and still is in parts of the community) that the Eliashberg framework might apply to these superconductors, but with exchange of a boson other than the phonon. Quite a few years later we wrote a review jointly, 34 focusing on "electron-phonon superconductivity." This review summarized known properties and extended results to dynamical properties such as the optical conductivity, building on earlier work in Ref. [35] and mini-reviews in Ref. [36]. More recently Ummarino has published a minireview with some generalizations to multi-band and the iron pnictide superconductors. 37 The other remark we should make is that while Eliashberg theory has been extremely successful, we will also point out the limitations that exist. Indeed, these were recognized right from the beginning, with both Eliashberg 12 and Migdal 38 emphasizing that limitations exist on the value of the dimensionless coupling parameter, λ, due to the expected phonon softening that would occur as λ increases. They claimed an upper limit of λ ≈ 1, which then significantly restricts the domain of validity of the theory. Constraints on the parameters would be a constant theme over the ensuing years. In 1968 McMillan 39 gave more quantitative arguments for a maximum T c , based on the expected relationship between the coupling strength and the phonon frequency. This was reinforced by Cohen and Anderson 40 and has been discussed critically a number of times since. 41,42 Alexandrov 43 has also raised objections, based on polaron collapse, a topic we will revisit later.
Some of the early history regarding the origins of the electron-phonon interaction was provided in Ref. [34] and will be omitted here. By the early to mid 1950's Fröhlich 44 and Bardeen and Pines 45 had established that the effective Hamiltonian for the electron-phonon interaction had a potential interaction of the form 46 (1) where k TF is the Thomas-Fermi wave vector, and ω(q) is the dressed phonon frequency. This part of the Hamiltonian represents the pairing interaction between two electrons with wave vectors k and k in the First Brillouin Zone (FBZ) and energies k and k . The interaction Hamiltonian written in this form is often said to have "the phonons integrated out." It was on the basis of this Hamiltonian that Bardeen, Cooper and Schrieffer (BCS) 10 formulated a model Hamiltonian with an attractive (negative in sign) interaction for electron energies near the Fermi energy, F .

II. THE ELIASHBERG EQUATIONS
Eliashberg, 12 following what Migdal 38 had calculated in the normal state, did not "integrate out the phonons," but instead adopted the Hamiltonian where c kσ (c † kσ ) is the annihilation (creation) operator for an electron with spin σ and wave vector k, and a q (a † q ) is the annihilation (creation) operator for a phonon with wave vector q. The electron-phonon coupling function, g(k, k ) is generally a function of both wave vectors (and not just their difference), and in principle is calculable with the Density Functional Theory Methods mentioned earlier. Very often models are adopted based on simple (e.g. tight-binding) considerations.
Eliashberg then applied the apparatus of field theory to formulate a pairing theory that accounts for the dynamics of the interaction, i.e. for retardation effects. A sketch of the derivation, taken from Rickayzen 47 (see also Ref. [34]), is provided in the Appendix. This is my favourite derivation, as it does not rely on a formalism (e.g. the Nambu formalism) whose validity requires an act of faith (or, you simply work through everything anyways, to ensure that the formalism "works"). Following Eliashberg 13 with more modern notation, 17,20 the "normal" self energy Σ(k, iω m ) is separated out into even and odd (in Matsubara frequency) parts, so that two new functions, Z and χ are defined: The equations that emerge are along with the equation for the order parameter: These are supplemented with the electron number equation, which determines the chemical potential, µ: Written in this way both Z and χ are even functions of iω m (and, as we've assumed from the beginning, they are also even functions of k). With electron-phonon pairing the anomalous self energy, which determines the anomalous pairing amplitude φ(k, iω m ), is also an even function of Matsubara frequency. A generalization of this latter result, giving rise to so-called Berezinskii 48 "oddfrequency" pairing, is beyond the scope of this review. A survey of Berezinskii pairing is given in Ref. [49].
Other symbols in Eqs. (4)(5)(6)(7) are as follows.The number of lattice sites is given by N , the parameter β ≡ 1/(k B T ), where k B is the Boltzmann constant and T is the temperature, µ is the chemical potential, and g( F ) is the electronic density of states at the Fermi level in the band. These equations are generally valid for multi-band systems, and then the labels k and k are to be understood to include band indices. However, we shall proceed for simplicity with the assumption of a single band, with single particle energy k . Because we are assuming finite temperature right from the start, the equations are written on the imaginary frequency axis, and are functions of the Fermion Matsubara frequencies, iω m ≡ πk B T (2m − 1), with m an integer. Similarly the Boson Matsubara frequencies are given by iν n ≡ 2πk B T n, where n is an integer. Finally, we have also included a direct Coulomb repulsion in the form of V k,k , which in principle represents the full (albeit screened) Coulomb interaction between two electrons.
The key ingredient of Eliashberg theory (as opposed to BCS theory) is the presence of the electron-phonon propagator, contained in with α 2 kk F (ν) the spectral function of the phonon Green function. This function is sometimes written as α 2 kk (ν)F (ν) to emphasize that the coupling part (α 2 ) can have significant frequency dependence. This spectral function is often called the Eliashberg function. Equation (8) has been used as a "bosonic glue" to generalize the application of the Eliashberg/BCS formalism to beyond that of phonon exchange. Very often the boson is a collective mode of the very degrees of freedom that are superconducting, i.e. the conduction electrons. Examples include spin fluctuations or plasmons, but this work is on more questionable footing. 50 A significant anisotropy may exist, specifically through the nature of the coupling in the Eliashberg function. Since the important physical attribute of the Eliashberg formalism beyond BCS is retardation, and therefore in the frequency domain, we will nonetheless neglect anisotropy in what follows. 51 More nuanced arguments for the wave vector dependence of the electronphonon coupling are provided in Ref. [20], connected to the energy scale hierarchy F >> ν phonon >> T c , where ν phonon is a typical phonon energy scale (note that we have adopted the standard practice of droppingh and k B , and therefore we refer to temperatures and phonon frequencies as energies). Indeed very often the neglect of anisotropy was justified by the study of so-called "dirty" superconductors, where the presence of impurities served to self-average over anisotropies. We will also drop the wave vector dependence in the direct Coulomb repulsion, although this step is less justified. It means that the direct Coulomb repulsion is represented by a single parameter, which we will call U , since this is what we would obtain by reducing the long-range Coulomb repulsion with an on-site Hubbard interaction with strength given by U . This is one of the weak points of the Eliashberg description of superconductivity -an inadequate description of correlations due to Coulomb interactions. In what follows we will focus more attention on the retardation effects, since this is the part that Eliashberg theory is best designed to handle properly.
Once we drop the wave vector dependence in the coupling function, all quantities (Z, χ and φ) become independent of wave vector. The integration over the First Brillouin Zone can then be performed, although here a series of approximations are utilized. The result can lead to confusion, so we provide some detail here. First, once it is determined that the unknown functions in Eqs. (4)(5)(6)(7) do not depend on wave vector, we can replace the sum over wave vectors in the first Brillouin zone with an integration over the electronic density of states, where g( ) is the single electron density of states and min and max are the minimum and maximum energies of the electronic band. Since typically the energy scales are such that T c << ν phonon << F , the variation of the electronic density of states away from the Fermi level is of little importance, so the approximation g( ) ≈ g( F ) is used. Exceptions to this case have been discussed previously, and one is referred to Ref. [34] for references. The remaining integration is now elementary and yields a combination of inverse tangent and logarithmic functions. 52 However, the "standard" practice is to extend these integrations over energy to ±∞ (and similarly adopt particle-hole symmetry so min ≡ − max , along with µ = 0, and now max → ∞) with the philosophy that the remainder of the integrand ensures that these additional contributions are negligible. This is in fact not true and the term proportional to V k,k in Eq. (6) (now replaced by U as described above) will result in a Matsubara sum that diverges if this procedure is carried out without thought. So in fact the bandwidth parameters are required in the integration for this term and the resulting inverse tangent function results in a soft cutoff at ω m ≈ max . This is most often replaced with a hard (step-function) cutoff. Because of the assumptions about particle-hole symmetry the function χ(k, iω m ) in Eq. (5) is identically zero, and Eq. (7) becomes meaningless (so the occupation is no longer considered an input parameter).
The more highly simplified equations that result are where W is the bandwidth (i.e. max = − min ≡ W/2 and the band has been centred around zero for convenience), and θ(x) is the usual step function, and u ≡ U g( F ). Very often a different order parameter is favoured over φ(ω m ), defined as ∆(ω m ) ≡ φ(ω m )/Z(ω m ). Then Eqs. (10,11) are written as These are the standard Eliashberg equations, written on the imaginary axis. These components of the self energy (refer back to Eq. (3)) are all real functions. One further 'simplification' is usually made before computations are performed. The Matsubara summation is in principle infinite; in practice the summation with the electronphonon kernel converges with relatively few terms, corresponding to a frequency cutoff ω c << W/2. Yet in the term with the Coulomb interaction one is required to carry out a summation over many more terms, corresponding to a frequency cutoff of W/2. Inspection of Eqs (12,13) shows that over this range (ω c < ω m < W/2), the gap function is a constant, i.e. ∆(ω m ) ≈ ∆ ∞ , and Z(ω m ) ≈ 1. Making use of this allows one to sum this part analytically, with the result that U changes to an and we have approximated digamma functions with their asymptotic logarithmic form, since it is assumed that both W/2 >> T c and ω c >> T c . Now the equations are where u * (ω c ) ≡ g( F )U * (ω c ). One should note that U * (ω c ) < U , physically corresponding to the fact that retardation effects allow two electrons to exchange a phonon with one another while not being at the same place at the same time. This means they do not feel the full direct Coulomb interaction with one another. Thus far we have written the Eliashberg equations as functions of imaginary frequency. As we will see in the next subsection one can solve these equations as they are, to determine many thermodynamic quantities of interest, in particular T c . However, later we will extend these equations to the upper half-plane, and in particular just above the real axis. This is required for the evaluation of dynamic quantities like the tunneling density of states and the optical conductivity. 17,34,36 In anticipation of these results we note here that we use functions Z(z) and φ(z) [and therefore ∆(z)] with the following properties 53 as a function of complex frequency z A. Results on the imaginary axis: Tc To compute actual results for T c , along with the gap function ∆(ω m ) and the renormalization function Z(ω m ), we need to specify α 2 F (ν) (now assumed to be isotropic) and u * (ω c ). The latter quantity is very difficult to compute, and the former is more tractable through Density Functional Theory. Historically it has been "measured" through tunnelling measurements. 18 We use quo-tation marks around the word "measured" because in fact the current is measured while the spectral function is extracted through an inversion process that requires theoretical input through the Eliashberg equations themselves. 18 We will simply adopt a model spectral function given by that is, a Lorentzian line shape cut off in such a way that the function goes smoothly to zero in the positive frequency domain. This Lorentzian has a centroid given by ν 0 and a half-width given by . The cutoff frequency parameter ν c makes the Lorentzian go to zero at frequency ν 0 + ν c and frequency ν 0 − ν c . For concreteness we will use a variety of values of ν 0 with = 0, or ≈ ν 0 /10. The first choice results in a δ-function spectrum with weight such that the mass enhancement parameter, λ, defined by is simply given by λ 0 . As increases λ decreases from λ 0 ; however in what follows we will adjust λ 0 to keep λ constant. 54 Since the main focus of Eliashberg theory is the effect of retardation, we will often set u * (ω c ) = 0, but we will nonetheless note how this quantity affects the gap function and T c . Superconducting T c is determined by linearizing the gap equations, Eqs. (15,16) so that they become The latter of these two equations is an eigenvalue equation and can be solved as such. We use a power method that iterates the eigenvalue and eigenvector simultaneously by requiring that the gap function at the lowest Matsubara frequency, 55 ∆(iω 1 ), remain at unity. This procedure tends to converge very quickly for stronger coupling, but requires more care for weaker coupling. 56 We begin with a standard plot of T c vs λ in Fig. 1, for a simple δ-function spectral function with frequency as shown. There are scaling relations for T c with typical phonon frequency, but we choose to show the results explicitly in real units to make the result clear. The The points marked with asterisks are again for the same broadened Lorentzian centred at ν0 = 100 meV, but now with u * (ωc = 1 eV. Note that νc is used as a practical cutoff for the phonon spectrum whereas ωc is used for the Matsubara cutoff for the direct Coulomb repulsion. See the discussion in the text for how plausible these parameters might or might not be. possibility of determining an expression for T c analytically has been discussed extensively in the literature 20 and will not be done here. The trends are clear; higher T c comes from higher values of λ and from higher values of the characteristic phonon frequency, which in the present case is provided by ν 0 . The width of the spectrum plays a minor role, and the Coulomb repulsion suppresses T c , as indicated by the marked points. It is well known that Eliashberg theory predicts that T c increases with both frequency and coupling strength as T c ≈ ν 0 √ λ in the asymptotic limit. 24,57 B. Validity of the Theory A perhaps more important question is the validity of the parameters used in the calculation. We have shown results up to a value of λ = 1. Are higher values allowed? In particular, is it possible for a material to have a sizeable value of electron-phonon coupling while maintaining a large phonon frequency? As discussed in the introduction, this question has been the subject of previous investigation, [40][41][42] although with only qualitative conclusions. The discovery of (very) high temperature superconductivity in the hydrides 58,59 under intense pressure has spurred a reassessment of this type of analysis, since a part of the community believes that these superconductors are electron-phonon driven. The main evidence has been an observed isotope shift. 58 Moreover, the prediction of superconductivity in some of these compounds through density functional theory calculations 60 adds plausibility to this explanation. However, very high characteristic phonon frequencies (60 -120 meV) and rather large electron-phonon coupling values (λ ≈ 2 or more) are required. The latter is well outside the range considered reasonable, especially given that the characteristic phonon energy remains so high. Moreover, the superconductivity literature has unfortunately lapsed into simply accepting as "standard" or "conventional" a value for the Coulomb pseudopotential u * = 0.1, and, especially given the high value of phonon frequency, the anticipated reduction of the Coulomb interaction through retardation will be much lower than previously thought, and the value of the direct Coulomb repulsion is undoubtedly higher when the phonon frequency is so high.
An additional direction of addressing this question comes from microscopic calculations involving Quantum Monte Carlo (QMC) and Exact Diagonalization (ED) techniques, utilizing specific microscopic models. These methods provide controlled approximations and are therefore suitable for benchmarking more approximate theories like Eliashberg theory. By far the most work in this direction has been done on the Holstein model. 61 The Holstein model retains only the short-range (on-site) interaction between local (Einstein) oscillators and the electron charge density. Because it is a very local model it is more amenable to the exact or controlled methods developed over the past 40 years, and therefore is a "favourite" for understanding the electron-phonon interaction, much like the Hubbard model 62 is heavily used for the study of electron-electron interactions. A short historical account of this activity is provided in the Appendix of Ref. [34].
Briefly, early Quantum Monte Carlo studies in one dimension 63 and two dimensions 64,65 established that charge-density-wave (CDW) correlations dominate at half-filling and close to half-filling. The critical question is whether, sufficiently away from half-filling, where the susceptibility for superconductivity is stronger than that for CDW formation, do the "remnant" CDW correlations enhance or suppress suppress superconducting T c ? In Ref. [66] the present author, based on a comparison of QMC and Migdal-Eliashberg calculations on (very!) small lattices, argued that CDW fluctuations actually suppress superconducting T c . In the so-called renormalized Migdal-Eliashberg calculations a phonon self-energy was included; in this manner CDW fluctuations impacted In the bottom figure the renormalized (solid curve) and unrenormalized (dashed curves) are plotted for a larger system. The renormalized calculations stop at an electron density close to n ≈ 0.9 because a CDW instability occurs there. The unrenormalized calculations carry on to half-filling, because they are oblivious to the CDW instability (and fluctuations). This result indicates that CDW fluctuations at densities less than n = 0.9 suppress pairing, and presumably Tc, even though λ eff → ∞. Reproduced from Ref. [66]. the phonon propagator, resulting in softer phonons and an enhanced coupling constant. Comparisons with the QMC results served to benchmark the Eliashberg-like calculations. Figure 2, reproduced from Ref. [66], illustrates that the so-called renormalized calculations agree with the QMC results. These calculations (solid curves) include phonon self-energy effects which are essentially the CDW fluctuations. 65 In contrast, the unrenormalized calculations (dashed curves) are the standard Migdal-Eliashberg calculations that omit phonon self-energy effects. Figure 2(a) illustrates that the renormalized calculations are more accurate, and Fig. 2(b) shows that including CDW fluctuations suppresses the pairing susceptibility, χ SP . We understand these results to indicate that in the vicinity of a CDW instability, while the effective coupling constant (λ eff in Ref. [65 and 66]) increases, superconducting T c actually decreases.
More recently similar calculations have been performed 67 and other methodologies have been employed. 68 In the latter reference the role of retardation in reducing the direct Coulomb interaction was also addressed; while they found qualitative agreement with the standard arguments, quantitative agreement was lacking, especially for the expected large values of direct Coulomb repulsion. An older calculation with just two electrons, 69 based on Exact Diagonalization studies, also found qualitative support for a retardation-related reduction in the direct Coulomb repulsion. It is worth mentioning that finding this insensitivity to increased Coulomb repulsion, known as the "pseudopotential effect," 70,71 has been looked for in QMC studies, but these have mostly been unsuccessful. They may still be there; part of the problem is that QMC results become more difficult as the electronic and phonon energy scales differ from one another by a significant amount. Moreover, it may be that if larger lattices and more realistic phonon frequencies (i.e. significantly less than the electron hopping parameter, t) are used, the result illustrated in Fig. 2 could change qualitatively.
An additional concern has been raised about the electron-phonon coupling becoming too strong -that of polaron collapse. 43 Exact studies in the thermodynamic limit 72 have established that a single electron, interacting with Einstein oscillators through the Holstein model, acquires an additional mass which is modest for λ < ≈ 1, but becomes quickly (though smoothly!) very, very large beyond this point. This is true independent of dimension, 73,74 and is especially acute when ν 0 << W , 75 where W is the electronic bandwidth. So we have the intriguing situation where the standard Migdal approximation (upon which Eliashberg theory is based) utilizes a Fermi sea of electrons strongly coupled to phonons, while close examination of just one of these constituent quasiparticles (polarons) shows that it acquires a tremendous effective mass. Thus, the properties of the single polaron, out of which a Fermi sea is constructed, appear to be incompatible with the properties of the electrons in that Fermi sea.
For example, in Migdal theory, the effective mass for electrons near the Fermi energy is m * /m e ≈ 1 + λ, even for λ ≈ 2 or more, whereas a single electron with this coupling would have an effective mass many orders of magnitude higher. It is important to note that in the quantum treatment polarons never "self-localize," essentially because of Bloch's Theorem. However, with such large effective mass ratios, any impurities (including surfaces), would readily act as localization sites.
There are perhaps a few scenarios to work one's way out of this dilemma. First, as we have already mentioned, perhaps the Holstein model itself is pathological. For this reason it is important to study other models. Unfortunately other models are more difficult to work with, but thus far the conclusions arrived at with the Holstein model seem to hold for these other models as well. For example, in Ref. [76] a variational approach was used with the Fröhlich Hamiltonian in the continuum with acoustic phonons, and in Refs. [77 and 78] the Barisić-Labbé-Friedel/Su-Schrieffer-Heeger (BLF/SSH) model was examined with perturbation theory and the adiabatic approximation. In either case more definitive results as achieved with the Holstein model are still lacking, although all indications are that these models have strong polaronic tendencies as well.
A second scenario is that as one assembles a Fermi sea of polarons, they somehow become increasingly undressed, presumably due to some argument involving Pauli blocking. There are no calculations that we are aware of, however, that provide a demonstration of this. 79 Part of the reason may be psychological; the Migdal approximation is more often called the Migdal Theorem, and so one may be inclined to take it for granted that this is what will happen when we assemble a Fermi sea -the "theorem" will be fulfilled. However, in my opinion this is more a belief than an established argument, as the Migdal approximation does not foresee or account for polaron physics.
This discussion has been a long digression concerning the domain of applicability of Eliashberg theory, and clearly a lot more investigation is required on this question. For now we return to the properties of the solutions to the Eliashberg equations.

C. Results on the imaginary axis: in the superconducting state
Returning to Eqs. (15,16), or their linearized counterparts, Eqs. (22,23), once T c is determined then the gap function can be determined both at T c and below T c . The gap function is a generalized (frequency-dependent) order parameter. It will grow continuously from zero at T c to its full value at zero temperature, but it depends on frequency. In Fig. 3(a) we show the gap function as a function of Matsubara frequency for a variety of temperatures, for λ = 1 and ν 0 = 10 meV. Note that ωm. Now the results look very similar to one another, which makes convergence from one temperature to the next relatively easy. Also shown is the weak coupling expectation 56 at Tc, ∆(iωm)/∆(iω1) = ω 2 E /(ω 2 E + ω 2 m ), indicated with a solid red curve. This result clearly does not resemble the data, since we the numerical results are for λ = 1. However, the slightly modified result, ∆(iωm)/∆(iω1) = ω 2 E /(ω 2 E + ω 2 m /(1 + λ) 2 ), shown as a dashed black curve, is a fairly good fit. these functions are defined on a discrete set of points (the Fermion Matsubara frequencies) that become more closely spaced as the temperature is lowered. For temperatures close to T c the gap function diminishes gradually to zero at all frequencies, while at the lowest temperature the gap function attains a maximum. In Fig. 3(b) we plot the normalized values, ∆(iω m )/∆(iω 1 ) vs. Matsubara frequency, and it is clear that they differ from one another by very little. Returning to Fig. 3(a), the lowest frequency function value can be thought of as an order parameter. In Fig. 4 we show the lowest frequency gap function value, ∆(iω 1 ), now normalized to the value at T = 0 vs. reduced temperature, T /T c . These are shown as blue squares, for about 9 temperatures. Also shown is the BCS weak coupling result, given as a red curve. One can see that the differences are small. We have also a plotted many more points (green asterisks) for a weaker coupling, λ = 0.3 (same phonon frequency), which fall exactly on the BCS curve. The main point is that deviations from the weak coupling BCS result are minor. For much stronger coupling than given here deviations are similarly very small, and experiment confirms this to be the case. 80 Many measurable properties of the superconducting state can be calculated from the imaginary axis solutions to the gap function. The renormalization function, Z(iω m ), is also required but this does not change by very much in the superconducting state. Examples of measurable properties include all thermodynamic quantities like the specific heat, and various critical fields. Systematic changes with coupling strength, as measured by λ, or alternatively the ratio of the critical temperature to a particular phonon frequency moment, T c /ω ln , have been reviewed elsewhere, 33,34 and will not be repeated here.

D. Results on the real axis
For any dynamical property (tunneling current, optical response, dynamical penetration depth, etc.), the relevant Green function (and therefore self-energy) is required as a function of real frequency. More precisely, for the retarded Green function it is needed at ω + iδ, i.e. infinitesimally above the real axis. In the original literature 12,13,71,81,82 the spectral representation was introduced to replace Matsubara sums with real frequency integrals, and these equations were then solved, either analytically (with approximations) or numerically. This was a difficult task (especially with the computers available at that time), and eventually the procedure was adopted that first required a solution on the imaginary axis and then analytic continuation (iω m → ω + iδ) through some approximate process. For this type of analytic function, the method of Padé approximants was used, 83 although the degree of precision needed for the gap function on the imaginary axis was very stringent (10 −12 for relative errors) in order to achieve accurate results on the real axis.
An appreciation for the information imbedded in imaginary axis solutions can be attained by considering the simple example of g(iω m ) = sech(ω m /ν 0 ), a very smooth function without structure, and monotonically decreasing with frequency on the positive imaginary axis. The analytic continuation can be easily done analytically; it is g(ω + iδ) = sec[(ω + iδ)/ν 0 ]. This function, in contrast to its imaginary axis counterpart, is riddled with divergences and discontinuities. And yet, in principle at least, this information is embedded in the (smooth) results on the imaginary axis. In practice, the information is contained in the 10 th significant digit and beyond.
An alternative, numerically exact procedure was developed in the late 1980's. 84 Here we simply write down the resulting expressions, the derivation of which is available in Ref. [84]. They are and and of course ∆(ω + iδ) ≡ φ(ω + iδ)/Z(ω + iδ). Here f (ω) ≡ 1/(exp(βω) + 1) and N (ν) ≡ 1/(exp(βν) − 1) are the Fermi and Bose functions respectively. Note that in cases where the square-root is complex, the branch with positive imaginary part is to be chosen. The reason for this can be traced back to Eq. (6) [or Eq. (4)] where the integration over k (with the assumptions made there) requires that the pole (given the same square-root that appears here) be above the real axis. It can easily be verified that substituting ω + iδ → iω n instantly recovers the imaginary axis equations (all the Fermi and Bose factors cancel to give zero contributions beyond the initial terms that require Matsubara summations). Clearly the inverse is not true -replacing the Matsubara frequency iω m where it appears in Eqs. (15,16) produces the first lines in Eqs. (24,25) involving Matsubara sums, but leaves out the remaining two lines in each case. The strategy for solving these equations is straightforward; the imaginary axis equations [Eqs. (15,16) ] are first solved self-consistently. These are then used in Eqs. (24,25) and these equations are iterated to convergence. The presence of the first lines in these equations provides a "driving term" that makes the iteration process quite rapid. For example, performing the entire operation (solution of imaginary axis and real axis equations) for a given temperature takes about a tenth of a second on a laptop.
Moreover, T = 0 is a special case, as is clear from these equations. In fact, this was recognized a long time ago, 85 where they established the following low frequency behaviour at T = 0, where c and d are constants. In contrast, the behaviour at any non-zero temperature is For conventional parameter choices this distinction has very little consequence, as the differences are barely discernible. For example, the expression of the imaginary part of Z(ω + iδ) in the normal state is given by for a δ-function spectrum [ → 0 in Eq. (20)] with strength λ and central frequency ν 0 . For ν 0 = 10 meV and λ = 1 then Fig. 1 indicates T c ≈ 10 K, and at T /T c = 0.1, the exponent in the Bose and Fermi functions makes this part ≈ 10 −50 . Thus, true to Eq. (27) the limiting behaviour is ∝ ν 0 /ω. However, to see this requires ω/ν 0 < 10 −45 , making it unobservable, and indistinguishable from the T = 0 case.
While this expression is for the normal state, the corre-sponding one of the superconducting state is even more severe, due to the development of the superconducting order parameter, which gives rise to a gap in the spectrum. Again, given the first two lines in Eq. (27) there is technically no gap, but in practice for reasons like we have just indicated, the practical results more closely follow the behaviour indicated in Eq. (26). Where the finite temperature result becomes pronounced and noticeably different than the zero temperature behaviour is in the strong coupling limit, 86 but in this case the parameters are not realistic and undoubtedly beyond the domain of validity of the theory. Returning to T = 0, for cases like the present where the phonon spectrum has a gap there is a special simplification. Basically, no iteration is required -the low frequency gap and renormalization functions come entirely from the first lines of Eqs. (24,25), and these can be constructed explicitly from the imaginary axis results. For a δ-function phonon spectrum, however, one has to be careful to convert the Matsubara summation to an actual integral, as a discontinuity will occur at the phonon frequency (at non-zero temperature this discontinuity is broadened into a gradual drop). Once the low frequency gap and renormalization functions are so constructed, higher frequency values require the Matsubara sum and real axis values of these functions at lower frequencies only. Eventually the entire functional forms are so constructed, without the need for iteration. As we will see, low temperature results converge quite rapidly to the zero temperature result, so this non-iterative method can be used as an alternative, for the lowest temperatures. Nonetheless, we will proceed with fully converged (iterative) finite temperature results, since these require such few iterations anyways.
To show real axis results, we utilize a phonon spectrum as in Eq. (20) with ν 0 = 10 meV, ν c = 8 meV, and = 1 eV, and λ = 1. Results with a δ-function spectrum tend to have a series of singularities, that are anyways artifacts of the singular spectrum, so we prefer to show results corresponding to this model spectrum. A series of such plots was also shown in Ref. [34] for the tunnelingderived Pb spectrum, with a much larger value of λ, and many more such results have been shown in the literature, often using this same method. 87 In Fig. 5 we show (a) the real part of the gap function, (b) the imaginary part of the gap function, (c) the real part of the renormalization function, (d) the imaginary part of the renormalization function, and (e) the tunneling density of states, which is measurable in single-particle tunneling experiments. The first observation, difficult to make with just these results, is that an image of the α 2 F (ν) spectrum is contained in both the real and imaginary parts of the gap function. Here it is the peak structure clearly evident in (a) centred around 10 meV for the highest temperature shown. As the temperature decreases this peak shifts to higher frequency, roughly by an amount equal to the value of the gap function at low frequency (about 2 meV in the present case). Experimentation with different spectral functions makes this observation more selfevident. See, for example, the distinctive spectrum for Pb in Fig. 4.35(a) of Ref. [34]. Both functions in (a) and (b) go to zero as the critical temperature is approached from below. As discussed earlier, they both go to zero at zero frequency at all temperatures shown, according to Eqs. (27), although one cannot see this on the scale shown. Even for the highest temperatures shown this behaviour can barely be seen, but becomes evident when one expands the low frequency scale. For the lowest temperatures shown even expanding the frequency scale by a few orders of magnitude is not enough to reveal the low-frequency behaviour indicated by Eqs. (27). For this reason one cannot use ∆(ω + iδ) as an order parameter at any frequency; either one has to revert to φ(ω + iδ) at zero frequency, or one can use the imaginary axis result for ∆(iω m ), as we did in Fig. 4. Also note that these functions approach zero at high frequency. If a Coulomb repulsion is included then the real part of ∆(ω + iδ) approaches a negative constant at high frequency.
In contrast the real and imaginary parts of Z(ω + iδ) plotted in (c) and (d) have changed very little in the superconducting state, and remain non-zero at the superconducting critical temperature, as indicated by the black curve. An image of α 2 F (ν) is present in this function as well, particularly in the imaginary part (see also Fig. 4.35 (c) and (d) in Ref. [34]). Finally, the tunneling density of states is shown in Fig. 5(e), and reveals a "gap" that opens from zero at T = T c rather quickly and then saturates to a low temperature value as indicated. In fact a plot of this "gap" vs. temperature would follow the result displayed in Fig. 4 very closely. However, as first pointed out by Karakozov et al. 85 there is no "gap" (hence the parentheses) and in fact this is evident in Fig. 5(e), where there is a noticeable rounding of the curves at frequencies below the sharp peak at almost all temperatures shown. The peak is a remnant of the square-root singularity known from BCS theory, which is evident from Eq. (29) if a constant gap function is used, ∆(ω + iδ) = ∆ 0 . In fact Eliashberg theory predicts smearing of this singularity simply due to the presence of imaginary components of all the functions involved in Eqs. (27) for all finite temperatures. It is also worth pointing out that the BCS limit of Eliashberg theory is not achieved by setting the gap function to a constant, ∆(ω + iδ) → ∆ 0 , but in fact the gap function is a decaying function of frequency in this limit. 85 This frequency dependence and its implications for the weak coupling limit has been further explored in Refs. [56, 88, and 89].
Finally, though not so evident in Fig. 5(e), there are "ripples" in the density of states beyond the "gap" region, caused by the coupling of electrons to phonons. The presence of these ripples in experiments (see, in particular, Refs. [18, 90, and 91]) is perhaps the strongest evidence of the validity of Eliashberg theory. In fact the most intense scrutiny has been superconducting Pb, where the electron-phonon coupling is particularly strong, λ ≈ 1.5, with a value well beyond the expected domain of validity. These experiments, coupled with an inversion technique that use the Eliashberg equations themselves, result in a consistent description of the superconducting state for Pb and other so-called "strong coupling" superconductors.
Many systematic deviations from BCS theory have been characterized, for example the gap ratio, 2∆ 0 /(k B T c ), 92 the specific heat jump, and many other dimensionless ratios 93,94 These have all been reviewed in Ref. [33], and show very systematic behaviour as a func-tion of the strong coupling index, T c /ω ln . On the other hand, when systematics are examined with purely experimental parameters, the picture is not so clear. 95

III. SUMMARY AND OUTLOOK
I have provided just a sketch of what we consider the essence of Eliashberg theory -retardation effects, in the context of a single featureless band. The generalization of these types of calculations to more complicated scenarios is well documented in a number of places, and have not been reviewed here. These include order parameter anisotropy, multi-band superconductivity, Berezinskii "odd-frequency" pairing, sharply varying electronic density of state, impurity effects, and so on. These additional complications are increasingly taken into account to understand new classes of compounds that exhibit superconductivity, such as the hydrides, MgB 2 , and the pnictides. In some cases, these additional effects have been invoked to explain higher critical temperatures as well, but for the most part they are motivated by matching theory to experiment.
In its bare form, Fig. 1 presents the possibilities for T c provided by Eliashberg theory. The conscious decision was made to extend the domain of coupling strength to unity only and not beyond, because there are reasons to believe that going beyond this regime is not viable. At the same time, large values of the characteristic phonon frequency have been used, and this is why the plot extends to beyond ≈ 50 K for the vertical axis, T c . Are these values of frequency, together with large values of λ ≈ 1 viable? Probably not, but given these sorts of parameter values, this is what Eliashberg in its standard form predicts.
I have also tried to touch on aspects of the theory where more critical scrutiny is possible, by comparing results to those obtained with microscopic models, such as the Holstein model. We believe there are significant difficulties that arise when these comparisons are made. One reaction is to dismiss such comparisons, as the Holstein model (or the Hubbard model, for that matter) may be regarded as "toy models," possibly fraught with pathologies. However, if the Holstein model is lacking in some way, it is important to know why, and what other aspect of the electron-phonon interaction (wave-vector dependence?) is essential to the success of Eliashberg theory. For example, if the super-high-T c of the hydrides is confirmed to originate in the electron-phonon interaction, then clearly one or more missing gaps in our understanding of how this happens needs to be filled.
Moreover, as presented here, Eliashberg theory focusses on the superconducting instability, and does not consider other, possibly competing, or potentially enhancing, instabilities. This possibility has come up as more and more phase diagrams of families of materials exhibit a nearby antiferromagnetic or charge-density-wave instability, as a function of some tuning parameter (dop-ing, pressure, etc.). It would be desirable to generalize Eliashberg theory to be more "self-regulating," and have the theory itself indicate when a competing instability is limiting superconducting T c , for example.
The other aspect that goes hand in hand with the electron-phonon coupling is the direct Coulomb interaction. We cannot claim to understand superconductivity to the point of having predictive power until we understand the role of Coulomb correlations, and their detrimental (or perhaps favourable?) effects on pairing. A key advancement has to come in further understanding the role that competitive tendencies or instabilities play in superconductivity. Many of the modern-day methods (Dynamical Mean Field Theory, for example) seek to address the question of competing interactions. Studies with Quantum Monte Carlo methods, like the ones mentioned here, will also aid in furthering our understanding of interacting electrons, and similar studies with the nowaccessible much larger lattice sizes would be welcome.

ACKNOWLEDGMENTS
This work was supported in part by the Natural Sciences and Engineering Research Council of Canada (NSERC). I would like to thank the many students and postdoctoral fellows in my own group who have contributed to the work described here. I also want to thank in particular Jules Carbotte, who first taught me about Eliashberg theory, and Ewald Schachinger, who was instrumental in teaching me about programming the Eliashberg imaginary axis equations. I also want to thank Sasha Alexandrov, who over the years continued to question the validity of MIgdal-Eliashberg theory. In the same way, Jorge Hirsch, with whom I have worked a great deal on other matters, has been instrumental in discussions concerning the validity of the work reviewed here. I also want to thank him for critical comments concerning parts of this review. Finally, Jules sadly passed away earlier this year (April 5, 2019), and this review is dedicated to his memory. He was a wonderful man, and I feel extremely fortunate to have first entered the physics world under his guidance.

APPENDIX: Derivation of Eliashberg Theory
In this Appendix, we will first outline a derivation of Eliashberg theory, based on a weak coupling approach. Our primary source for this derivation is Ref. [47]. Migdal theory of the normal state follows by simply dropping the anomalous amplitudes in what follows.
If we know the many-body wave function of system, we can calculate the expectation value for any observable. However, usually this is something we do not know, and instead we calculate multi-electron Green functions, which are themselves related to observables. The Green functions are necessarily almost always approximate, and those computed in Eliashberg theory are no exception. In fact, Eliashberg theory is essentially a mean-field theory, though because of the inherent frequency dependence in the self energy, it is in many ways a precursor to Dynamical Mean Field Theory. 96 We begin with the definition of the one-electron Green function, defined in momentum space, 97 as a function of imaginary time, τ , where k is the momentum and σ is the spin. The angular brackets denote a thermodynamic average. We can Fourier-expand this Green function in imaginary frequency: The frequencies iω m are the Fermion Matsubara frequencies, given by iω m = iπT (2m − 1), m = 0, ±1, ±2, ..., where T is the temperature. Because the c's are Fermion operators, the Matsubara frequencies are odd multiples of iπT. The imaginary time τ takes on values from 0 to β ≡ 1/(k B T ). A similar definition holds for the phonon Green function, where A q (τ ) ≡ a q (τ ) + a † −q (τ ).
The Fourier transform is similar to that given in (31) except that the Matsubara frequencies are iν n ≡ iπT 2n, n = 0, ±1, ±2, ... and occur at even multiples of iπT . These are the Boson Matsubara frequencies.
To derive the Eliashberg equations, we follow Ref. [47], and use the equation-of-motion method. The starting point is the time derivative of Eq. (30), where we have put τ = 0, without loss of generality. We use the Hamiltonian (2), and assume, for the Coulomb interaction, the simple Hubbard model, H Coul = U i n i↑ n i↓ . Including the Coulomb repulsion, the result is repeated here, where the various symbols have already been defined in the text. We consider only the Green function with σ =↑; the commutator in Eq. (34) is straightforward and we obtain Various higher order propagators now appear; to determine them another equation of motion can be written, which would, in turn, generate even higher order propagators, and this eventually leads to an infinite set of equations with hierarchical structure. This infinite series is normally truncated at some point by the process of decoupling, which is simply an approximation procedure. For example, in (36) the Coulomb term is normally decoupled at this point and becomes The case of the electron-phonon term is more difficult; in this case we define a 'hybrid' electron/phonon Green function, G 2 (k, k , τ, τ 1 ) ≡< T τ A k−k (τ )c k ↑ (τ 1 )c † k↑ (0) >, (38) and write out an equation of motion for it. We simply get ∂ ∂τ G 2 (k, k , τ, τ 1 ) = −ω k−k < T τ P k−k (τ )c k ↑ (τ 1 )c † k↑ (0) >, (39) where P q (τ ) = a q (τ ) − a −q (τ ). Taking a second derivative yields At this point we do not simply decouple the last line of Eq. (40). We first need to take the phonon propagator into account, and the standard procedure is to use the "non-interacting" phonon propagator. The adjective "non-interacting" is in quotes because part of the philosophy of proceeding in this way was a desire to not compute corrections to the phonon propagator, because the information going into this part of the calculations (e.g. the phonon spectral function) was going to come from experiment. Coming from experiment means that nature "had already done the calculation," and we did not want to double count. Clearly, if the purpose of this calculation is to compare to Quantum Monte Carlo calculations where this is not the case, then something different should be done, and this is what motivated the renormalized Migdal-Eliashberg calculations of Refs. [65 and 66]. For now, we proceed with the standard Eliashberg calculations. The equation of motion for the non-interacting phonon propagator is ∂ 2 ∂τ 2 − ω 2 q D(q, τ − τ ) = 2ω q δ(τ − τ ).
Inserting this expression into Eq. (40) then yields where now τ 1 has been set equal to τ as is required in (36). It is important that this be done only after applying Eq. (41). The result can be substituted into Eq. (36), and then Fourier transformed (from imaginary time to imaginary frequency). Before doing this however, we recall that Gor'kov 11 realized the important role of the so-called Gor'kov anomalous amplitude, in the Wick decomposition 97 of the various two-particle Green functions encountered above. We therefore have to account for these in addition to the pairing of fermion operators used in Eq. (37). The anomalous amplitudes are defined to be F (k, τ ) ≡ − < T τ c k↑ (τ )c −k↓ (0) > andF (k, τ ) ≡ − < T τ c † −k↓ (τ )c † k↑ (0) > .
Now we need to repeat the same steps as above with F andF as we did with G. Skipping the intermediate steps, the result is an equation analogous to Eq. (36) and similarly for the function F . This leads to the need for another 'hybrid' electron/phonon anomalous Green function, and, following the same process as for the regular Green function, we find where again τ 1 has been set equal to τ after applying Eq. (41). The Fourier definitions of the anomalous Green function are the same as Eq. (31): and similarly for F . In frequency space one finds that two self energies naturally arise, One can show that g k k = g * kk ; normally one expects a similar relation with negative wave vectors, and we assume it in what follows. These equations are then written self-consistently and lead to Eqs. (4-7) once Eq. (3) is used.