Adiabatic theory of the polaron spectral function

An analytic theory for the spectral function for electrons coupled with phonons is formulated in the adiabatic limit. In the case when the chemical potential is large and negative $\mu \to -\infty$ the ground state does not have the adiabatic deformation and the spectral function is defined by the standard perturbation theory. In this limit we present the spectral function calculated up to fourth order in the electron-phonon coupling constant which satisfies the sum rules up to the 5th moment. In the case when the chemical potential is pinned at the polaron binding energy the spectral function is defined by the ground state with a nonzero adiabatic deformation. We calculate the spectral function with the finite polaron density in the adiabatic limit. We also demonstrate how the sum rules for higher moments may be evaluated in the adiabatic limit. Contrary to the case of zero polaron density the spectral function with the finite polaron concentration has some contributions which are characteristic for polarons.


INTRODUCTION
The properties of different types of polarons were well studied a long time ago. A number of review articles (see for example [1][2][3][4][5][6][7]) and textbooks [8][9][10][11][12][13] is available, which describe the spectroscopic, thermodynamic, kinetic, and other physical properties of polarons. Usually, the polaron theory is used in order to describe electric transport in low mobility crystalline or organic semiconductors (see review article of I. G. Austin, N. F. Mott [1] and references therein or more recent papers on organic semiconductors [14]). It was also successfully used in order to describe equilibrium and photo-induced mid-infrared optical absorption spectra of high-T c superconductors at low doping [15][16][17][18]. The description is based on the well developed theory of polaron optical absorption [19][20][21][22][23]. The small Jahn-Teller polarons [24] are also found in colossal magnetoresistance manganites [25,26].
Note that until recently most investigations were related to the ground state properties of polaronic systems, such as polaron binding energy, effective mass, transport properties as well as mid-infrared absorption. The development of angle resolved photoemission spectroscopy (ARPES) demands the analysis of the whole spectral function of polaron systems [41]. This problem was successfully addressed by a number of different numerical techniques [27-29, 39, 42-47]. Nevertheless, most of these techniques are restricted to studying the so called Holstein model [48] of molecular crystals, where both the electron-phonon coupling constant and the phonon frequency is momentum independent. The spatial dispersion of the electron-phonon interaction leads to a substantial increase in the number of nonzero matrix elements in the Hamiltonian and poses substantial restrictions for this type of calculation. For the same reason, calculations are usually performed in the non-realistic limit for crystalline materials when the phonon frequency ω 0 is of the order of the hopping integral t for electrons [28-31, 42, 43, 47]. Moreover, the calculations are usually restricted by a 1D finite system. The system size is usually restricted by N ∼ 10 [47] leading to poorly controllable finite size effects. On the other hand, it is well known that the polaron formation crucially depends on the dimension of the system [13,49].
The investigation of the polaron spectral function was stimulated when the exact sum rules for dilute electronphonon systems were formulated by Kornilovitch [50] (see also Ref. [39]). It was shown that the spectral function proposed in Ref. [41] satisfies only the zero order sum rule. On the other hand, the numerical calculations in Ref. [47] obey the sum rules derived for the limit of zero polaron density n = 0. In Ref. [39] an uncontrollable approximation of the polaron Greens function which neglects all momentum dependence of the self-energy was proposed.
Here we present an adiabatic theory for the spectral function of the Holstein model. The theory is based on equations, formulated for the case of the Holstein model [48] in Ref. [49]. These equations are similar to that, derived by Pekar for the polaron in polar crystals [8]. The equations have two different sets of solutions [49]. Near the trivial solution in the limit of zero polaron density n = 0 (chemical potential µ → −∞) we derive an equation for the self-energy. Vertex in this equation is found from an equation that accurately takes into account threshold effects.
This theory corresponds to the summation of infinite series of diagrams where each phonon line has not more than three crossings. The theory takes into account exactly all diagrams up to 6th order. As a result, the spectral function obeys the sum rules up to the 7th order and has accuracy better than 3% for λ ≤ 3 in the adiabatic limit. Here λ is dimensionless electron-phonon coupling constant. It is interesting to note that the polaron contribution to the spectral function is at least exponentially small ∼ exp (−const M/m) where m is the effective mass of the electron and M is the ionic mass. As it was mentioned in Ref. [39] these small terms are absent in the sum rules for the spectral function at zero polaron density, indicating that contributions from the polaron state to the spectral function are negligible. For the case when the chemical potential is pinned to the polaron binding energy and small but finite concentration of polarons n = 0, the ground state has nonzero adiabatic deformation. In this limit, we derive the spectral function which obeys the exact known sum rules for finite polaron density.
The paper is organized as follows. In the next section, we briefly discuss some important details of the polaron theory, and then in the section "Results" we discuss the polaron spectral function and the sum rules in the limit of zero and finite polaron density.

Electron-phonon interaction Hamiltonian
The Hamiltonian of interacting electrons and phonons has the form: where c k and b q are the electron and phonon annihilation operators with momentum k and q respectively, ǫ(k) and ω q are the electron band energy and the phonon frequency, and N is the number of sites in the lattice. We also assume here that = 1. In the following, we consider dispersionless phonons with ω q = ω 0 . The dimensionless matrix element of the electron-phonon interaction γ(q) has different q dependence for different types of crystals. In ionic crystals (the Fröhlich model [51]) with strong dispersion of the dielectric permittivity in the long wavelength limit |γ(q)| 2 = 4πe 2 /2κω 0 a 3 q 2 , where e is the elementary charge, κ −1 = ε −1 ∞ − ε −1 0 , ε 0 , ε ∞ are the static and the high frequency dielectric constants, and a is the lattice constant. In the case of molecular crystals the Holstein model is valid and γ(q) = g is momentum independent.
For the vast majority of crystalline solids, the adiabatic approximation is valid [52]. It means that the ratio of the electron band mass to the ionic mass is a small parameter. Indeed, except for compounds with heavy fermions, the effective mass of electrons or holes is of the order of the free electron mass. Ion masses are at least 1700 times larger. In organic solids, the applicability of the adiabatic approximation may be questionable. Indeed, organic materials which contain large molecules may have relatively small hopping integrals t ∼ 0.1 eV because the molecular orbitals are spread over the whole molecule. On the other hand, intramolecular vibrations may be as hard as 0.1-0.2 eV because of the presence of light atoms in the molecule.
Let us estimate the matrix elements in Hamiltonian Eq.(1) in the Fröhlich model [51]. The interaction between ions is dominated by the Coulomb attraction V (r) ∼ e 2 /ε ∞ r. The phonon frequency is determined by the second derivative of V (r) at r = a, ω 0 ∼ e/ε It is easy to see that the ratio e 2 ma ε∞ ≈1. It means that everywhere in the Brillouin zone except the close vicinity of q = 0 point we have the following hierarchy t : γ(q)ω 0 : ω 0 ∼ 1 : ( m M 1/4 : ( m M 1/2 . This hierarchy is also valid in the case of the Holstein model. The electron-phonon matrix element in the Hamiltonian is always smaller than the kinetic energy of electrons. Nevertheless, from the second and the third terms in Hamiltonian (1) we can construct adiabatic (i.e. ion mass independent) energy which is called the polaron binding energy or the polaron shift. This energy does not depend on ionic mass and may be as large as t or even larger. If following the works of Eliashberg [53] we define the Eliashberg function and determine the dimensionless electron-phonon coupling constant it is determined exactly by Eq. (2) λ = E p /zt (z is the number of nearest neighbors in the lattice).
One can apply perturbation theory to Hamiltonian (1). In the continuum approximation of the Fröhlich model when t → ∞ a → 0 keeping the effective mass 1/ta 2 constant the perturbation theory is the expansion in the dimensionless parameter α = e 2 κ m 2ω0 [51,54]. In the lattice Holstein model the self energy represents the expansion in the dimensionless parameter g 2 ω 2 0 /t 2 ∝ m/M [55]. Some diagrammatic expansions for the polaron self energy were used in the past in order to identify the instability associated with the polaron formation [28,56]. Nevertheless, within the standard perturbation theory the clear instability was not found.
Translation symmetry broken ground state equations in the adiabatic limit M → ∞ The adiabatic theory was initially formulated for the Holstein model Eq.(1) with γ(k) = g in Ref. [48]. Later the theory was reformulated using the field theory in the vicinity of nonzero classical solutions [57] in Ref. [49]. This theory allows also us to analyze the nonadiabatic corrections to the adiabatic solutions. In Ref. [49] the equations which describe the saddle points of the adiabatic potential were derived. The central equation in the adiabatic theory is the Schrödinger equation for an electron moving in the external potential of the lattice deformation. In the discrete lattice case (the tight binding approximation) it has the form [49]: Here ψ n is the electronic wave function on the site n, ϕ n is the deformation at the site n, l describes the quantum numbers of the problem, and the summation over m is taken over the nearest neighbours. The important assumption of the adiabatic approximation is that the deformation field is very slow and we assume that ϕ n is time independent, ∂ϕ/∂t = 0 in Eq.(3). Therefore, ω 0 ∝ M −1/2 → 0 and g 2 ∝ M 1/2 → ∞ but the polaron shift g 2 ω 0 = E p is finite. The equation for ϕ n has the form [49]: After substitution of Eq.(4) to Eq.(3) we obtain: Similar equations may be formulated in the case of Fröhlich model [51]. In the continuum case, it was done by Pekar [8] and later for the discrete lattice (see for example Refs. [58,59]). The first equation represents the Schrödinger equation for an electron in the potential φ n generated by the displaced ions: The equation for the potential φ n reads: Note that the left hand side of Eq.(7) represents the discrete version of the Laplacian △. If we solve Eq.(7) and substitute the solution back in to the Schrödinger equation we obtain exactly the equation derived by Pekar [8]. Note, that Eqs. (3,6) define the energy spectrum of the electron in the deformation ϕ n or polarization φ n fields. The total energy must include the positive energy of the deformation and the polarization itself. In the case of the Holstein model the deformation energy is defined as E def = ω 0 n |ϕ n | 2 /2. Equation (5) has two types of solutions. The first solution breaks the translation invariance and corresponds to the self-trapped state. This is the solution that has nonzero deformation when the electron occupies the ground state level. The properties of the self-trapped solutions depend strongly on the system dimensionality. In the 1D case, the self-trapped solution exists at any value of the polaron shift E p . Therefore, the polaron is always stable in 1D. In 2D and 3D, the self-trapped solution exists only when the electron-phonon coupling constant λ > λ c , where λ c is of the order of 1 and depends on the system dimensionality. The self-trapping may occur with or without barrier formation depending on the value of the coupling constant λ and the system dimensionality [13,49,60].
The first correction to the adiabatic solution describes the renormalization of the local phonon mode. The electron self-trapping causes an increase of the local density near the polaron center and leads to the shift of the local vibrational frequency. In the strong coupling limit E p ≫ t the renormalized mode has the frequency [49]: here z is the number of nearest neighbours. The second nonadiabatic correction describes the tunneling of self-trapped polarons. The tunneling splitting was first derived by Holstein in Ref. [48]. A more comprehensive formula was derived in Ref. [28] (see also Ref. [61]). In the adiabatic limit the polaron tunneling is exponentially suppressed, The second type of solution preserves translation invariance and represents itinerant band states. The deformation around the electron is absent ϕ n = 0 therefore all nonadiabatic corrections to this solution may be described in terms of ordinary perturbation theory and are described by the standard diagram technique, taking into account that m/M is small.

The spectral function.
As it follows from the previous section the eigenvalues and eigenstates of Hamiltonian (1) in the adiabatic limit m/M → 0 are determined by the two sets of states. The first one is determined by a nonzero lattice deformation in the ground state and represents all eigenstates of the Hamiltonian in the presence of this deformation. This type of eigenstate corresponds to a isystem with one polaron in the ground state of the grand-canonical Hamiltonian. The second set of eigenstates is the eigenstates of the Hamiltonian without any lattice deformation. These states describe the empty system and electrons may appear in the system only due to an external perturbation like photoexcitation or due to an injection. The parameter which controls the density of carriers is the chemical potential µ. The spectral function is defined as [62,63]: here G R (k, ω), G A (k, ω) is the retarded and advanced Green's functions of an electron. Using the Lehmann representation the spectral function may be written as [63][64][65]: where |i represents the eigenstates of the grand-canonical Hamiltonian k c k is the particle number operator and µ is the chemical potential. Z = n exp (−βE n ) is the grand-canonical statistical sum, β = 1/k B T (k B is the Boltzmann constant and T is temperature). The zero temperature limit of this formula has a direct physical meaning: here |0 and E 0 are the ground state and the ground state energy of the grand-canonical Hamiltonian H ′ . The first term in Eq.(11) describes the inverse photo-emission spectrum, when the number of electrons in the ground state is increased by 1. Contrary, the second term describes the direct photo-emission spectrum, when the number of electrons is reduced by 1.
Very often only the first term in Eq. (11) is calculated [28,29,39,47] assuming that the chemical potential is large and negative (µ < 0). Indeed, in that case, |0 is the phonon vacuum without any electrons, and therefore the second term in Eq.(11) is equal to zero. Calculating eigenstates and eigenvalues of Hamiltonian Eq.(1) with 1 electron it is easy to evaluate Eq. (11). On the other hand, the calculation of the spectral function at arbitrary µ requires diagonalization of many-body grand-canonical Hamiltonian H ′ with an arbitrary number of electrons, which is quite difficult task even for the exact diagonalization calculations [28,29,47].
Note that the spectral density, calculated with the large negative chemical potential is not sensitive to the polaron states, which are described by the solutions of Eq.(5) with a nonzero lattice deformation in the ground state of the grand-canonical Hamiltonian. Indeed, the overlap of the wave functions with and without the polaron deformation is exponentially small ∝ exp (−const M/m). That is the reason why the sum rules for the spectral density do not contain any adiabatic and nonadiabatic contributions, related to the polaron formation. To find polaronic features in this spectral function it is necessary to perform calculations with exponential accuracy in the adiabatic parameter m/M .
In the next section, we present the calculations of the spectral function in the limit of large and negative chemical potential. Then we assume that chemical potential is fixed in close vicinity to the polaron level in the way that the ground state has exactly one polaron. Here we consider spinless electrons and assume, that bipolaron formation is suppressed.

RESULTS.
The spectral function for µ → −∞ and the sum rules.
In this section, we construct the spectral function for the Holstein model [48] in 1D, where the polaron represents the excited state of the grand canonical Hamiltonian and the spectral function is given by the first term in Eq. (11).
Here we neglect all exponentially small terms, i.e. the overlap between the wave functions with zero and nonzero lattice deformation is neglected. Therefore, only terms with ϕ n = 0 give nonzero contributions to the spectrum, given by Eq. (3). This theory corresponds to the so-called sudden approximation. Note that it is exactly the reason why the sum rules do not have any contribution associated with the polaron formation as mentioned in Ref. [39]. Solution of Eq.(3) with ϕ n = 0 in 1D case gives the spectrum ǫ(k) = −2t cos (ka), where k is the electron momentum. Therefore, the spectral function in the adiabatic limit is given by the formula: This spectral function satisfies the sum rules M n = ∞ −∞ ω n A(k, ω)dω for n = 0, 1, derived by Kornilovich for the case of zero polaron density n → 0 [50]. Here we do not write the chemical potential explicitly and absorb it to the definition of ω. Note that all higher moments of the spectral function contain explicitly the terms proportional to the electron-phonon coupling constant and are proportional to some powers of (m/M ). To demonstrate that, we write the third sum rule M 3 = ǫ(k) 3 + 2ǫ(k)g 2 ω 2 0 + g 2 ω 3 0 . Taking into account that E p = g 2 ω 0 is independent of the ion mass M we conclude that the second term in the expression for M 3 is proportional to (m/M ) 1/2 and the third term ∝ (m/M ). Higher order moments M n have the mass independent term ǫ(k) n and the sum of m/M terms of powers l where 1 ≤ l ≤ n − 1. Therefore, if we neglect all nonadiabatic terms Eq.(12) satisfies all sum rules.
We can demonstrate that the spectral function of the Holstein model is indeed represented by the single δ-function in the adiabatic limit by plotting in Fig. 1 the spectral function of the two-site Hamiltonian A(k = 0, ω) (See for example Ref. [27]). This figure demonstrates that when ω 0 /t → 0 the spectral density at k = 0 is represented by a single peak centered at −t. At k = π/a the spectral function is converging to a single peak at ω = t. The width of this peak is decreasing to 0 when ω 0 /t → 0. Since M 0 = 1, this corresponds to the definition of the δ-function.
In the case when the ground state does not have polarons the solution of Eq.(4) is trivial ϕ n = 0. The nonadiabatic correction to the spectral function may be calculated by the standard perturbation theory. The dimensionless parameter, which determines the perturbation series is gω 0 /2t = λω 0 /2t ∝ (m/M ) 1/4 . Therefore, the perturbation theory may be considered as an expansion in series over the adiabatic parameter (m/M ) 1/2 because the expansion contains only even powers of the coupling constant. In order to calculate the spectral function with the given accuracy δ it is necessary to calculate all irreducible diagrams for the self-energy up to the order n which satisfies the inequality (gω 0 /2t) n = (λω 0 /2t) n/2 < δ. If we consider realistic parameters for crystalline materials ω 0 /2t ≤ 0.1 for required accuracy better than δ = 10% for λ ≤ 3 we have to take into account fourth order diagrams for the self-energy. Indeed, λ 2 (ω 0 /2t) 2 ≤ 0.09 < 0.1.
Nevertheless, in the 1D Holstein model, the self-energy diverges near the threshold of the phonon emission. This divergence becomes even more pronounced in higher orders of the perturbation theory. It is easy to see by the comparison of the contribution of the second-order and the fourth-order diagrams for the self-energy. The second order contribution diverges near the threshold as |x| −1/2 while the fourth order contribution diverges as |x| −3/2 , here x = 1 − (ǫ − ω 0 )/2t describes the deviation from the threshold. It demonstrates that in order to describe the behaviour of the self-energy correctly the summation of all diagrams which contain the same divergencies should be performed [66,67]. Following the procedure [66,67] we formulate Dyson's equation for the electron self-energy Σ(ǫ, k) (Fig.(2) where the vertex part Γ(ǫ, k, q) satisfies the equation schematically represented in Fig.(3). The kernel K(ǫ, k, q, x) of the integral equation for the vertex part Γ(ǫ, k, q) is represented by the square, which contains all irreducible diagrams with one incoming and one outcoming electron lines and one incoming and one outcoming phonon lines. The kernel may be evaluated by perturbation theory, schematically represented in Fig.(4). If we sum all diagrams contributing to this kernel we present the exact solution of the problem with ϕ n = 0. The exponentially small terms corresponding to the overlap of the states with and without deformation cannot be evaluated in this procedure because of the nonanalytic nature of these terms. The diagrams shown in Fig.(4) allow evaluating vertex (Fig.(3) which is exact in sixth order of perturbation expansion and correctly describes the threshold effects. Finally, the solution of Dyson's equation (Fig.(2)) with the vertex defined by Fig.(3) and with the kernel defined by diagrams, plotted in Fig.(4) represents the summation of all diagrams, where the number of crossings in each phonon line is less than four. This theory is accurate up to the sixth order in the dimensionless parameter (λω 0 /2t) 1/2 ∝ (m/M ) 1/4 and correctly describes the behaviour of the self-energy near the phonon emission threshold. Therefore, the accuracy of our calculations for λ < 3 and ω 0 ≤ 0.2t is better than 3%.
After integration over energies Dyson's equation represented in Fig.(2) and the equation for the vertex part represented in Fig.(3) have the form: and The kernel K(ǫ, k, q, x) calculated up to the sixth order in perturbation theory (Fig.(4)) is represented by the equation: Here the electron Green's function is defined as: As a result, the problem of calculation of the spectral function is reduced to a solution of two integral equations for the self-energy Eq.(13) and for the vertex part Eq. (14) with the kernel, defined by Eq. (15). Here we present the numerical solution to these equations. We also propose an accurate approximation for the vertex part and compare the numerical solution with the approximate solution.
In order to solve Eqs. (13,14) the energy was defined in 3001 points on the interval −6t < ǫ < 6t. The momentum was defined in 101 points of the Brillouin zone −π < k < π with the step ∆h = 2π/100. The integration was performed by the trapezoidal rule with the accuracy ∼ (∆h) 2 better than 1%. As a staring point of the iteration procedure, we introduce the vertex and the self-energy which are averaged over momenta and represent the solution of Eqs. (13,14) averaged over momenta Here z n is the smallest (|z n | ≤ 1) root of the quadratic equation z 2 + ǫ n /tz + 1 = 0, and ǫ n = ǫ − nω 0 − Σ(ǫ − nω 0 ). The self-energy Σ(ǫ) is defined from Eq.(13) with approximate vertex Γ av (ǫ). Note, that the Green's function Eq. (16) with this self-energy satisfies the sum rules up to the seventh order. The approximate self-energy and vertex part are used as a starting point to solve Eqs. (13,14) iteratively. The iteration procedure looks as follows. On every step the self-energy is calculated from Eq.(13) with the vertex from the previous step. Then the new self-energy is used to obtain the new vertex from Eq. (14). Dyson's equation Eq. (13) was solved by iterations. On the other hand, the standard routine for the solution of the system of the linear equation from the NAG library was used in every step of the iteration. Usually, only a few iterations are necessary to obtain solutions for the vertex part and the self-energy Eqs. (13,14). The results of calculations are presented in Fig(5a).
To construct an approximate solution for the vertex Γ(ǫ, k, q) we notice first, that the main contribution to the integral in Eq. (14) comes from the vicinity of x = k point. Note, that k − x = 0 point represents the Van Hove singularity in Eq. (14), because at this point group velocity is equal to 0. Therefore, we may take out the vertex part Γ(ǫ, k, k) from the integral. As a result, we obtain: where F (ǫ, k, q) = π −π K(ǫ, k, q, x)G(ǫ − ω 0 , k − x) dx 2π Rewriting this equation at q = k and then solving equation for Γ(ǫ, k, k), we obtain approximation for the vertex: Γ app (ǫ, k, q) = gω 0 + gω 0 F (ǫ, k, q)/(1 − F (ǫ, k, k)). Integrals in formula for F (ǫ, k, q) may be calculated analytically (See Appendix A). Substituting this vertex with analytic formulae for F (ǫ, k, q) Eqs. (29,30,31,32) to Dyson equation Eq.(13) we calculate spectral function, presented in Fig.(5b). Comparison of spectral functions presented in Fig.(5a) and Fig.(5b) demonstrates that indeed the main contribution to the vertex Γ(ǫ, k, q) is coming from the vicinity of the Van Hove singularity in Eq. (14). Therefore, there is very good agreement between numerical solution of Eqs. (13,14,15) and approximate results for the spectral function. Note, that both spectral functions satisfy the sum rules up to the seventh moment (Fig(6)). Fig.(5) clearly demonstrates that the lowest energy band has minimum at ǫ = ǫ(k = 0) ∝ (E 2 p ω 2 0 /t) 1/3 ∼ (m/M ) 1/3 . With increasing k the energy disperses relatively quick up to the energy of the order of ǫ(k) = ǫ(k = 0) + ω 0 at k ≈ √ 2mω 0 and then remains unchanged with further increase of k up to the Brillouin zone edge. This effect is well known in literature [66] and called the threshold phenomenon. Note that the threshold phenomena were first discussed in connection with the excitation spectrum of the superfluid He by Pitaevskii [67]. In the case of the 1D Holstein model, the anomaly in the self-energy belongs to type "c" as discussed by Levinson and Rashba [66] and more recently by Goodvin and Beciu [46]. The self-energy shows a strong divergence −(−2t + ω 0 − ǫ) −1/2 when ǫ → −2t + ω 0 − 0. This first branch of the spectrum does not have any damping, because any inelastic scattering of electrons is forbidden by the conservation of energy. The anomalies similar to the threshold effect are also present at higher energies because the higher order diagrams have additional divergence at −2t + nω 0 , n = 1, 2, 3, .... These anomalies are much less pronounced because of the finite imaginary part of the self-energy at these energies. And finally, there is an increase of A(k, ω) near ka = π and ω > 2t. This peak is less pronounced than the peak at the bottom of the band, because it has finite damping. Now let us compare the present theory, which takes into account all diagrams up to the sixth order and represents the partial summation of infinite perturbation series with the momentum average approximation, proposed in Ref. ([39]). This theory is accurate in the second-order approximation. The fourth-order contribution in this theory is approximated by the self-energy contribution: which is equal to an averaged over k fourth order diagram with phonon lines crossings (22), and I(z) = 2/(z − z −1 ). This expression should be compared with the exact forth order contribution Σ 4 (k, ǫ) = Σ N C 4 (k, ǫ) + Σ C 4 (k, ǫ), where  here z 1 and z 2 (|z 1,2 | < 1) roots of quadratic equations, defined after Eq. (17) and are calculated with Σ B (ω) = −iδ.
Comparing these results we conclude that momentum average approximation is rather poor because it provides the self-energy which is functionally different from the perturbation theory. Therefore, momentum average approximation is accurate only in the second order in coupling constant. There were few attempts to improve this approximation [68]. The first step of corrections leads to the correct expression for fourth order noncrossing diagram Σ N C 4 (k, ω) Eq.(21) but fails to take into account k-dependence of the diagram with crossing of phonon lines Σ C 4 (k, ǫ). The next and the last level of corrections leads to some infinite system of inhomogeneous equations, which was solved by truncation. This procedure is not well justified because the coefficient in this equation does not fall with the increasing of the system size. Nevertheless, the solution of this system expanding in powers of the coupling constant may be performed analytically and we recover the correct expression for the fourth-order diagrams. Higher order diagrams require the next level of corrections which were not discussed in Ref. [68]. Therefore, we conclude that momentum average approximation with the corrections of levels 1 and 2 is accurate up to the fourth order in perturbation theory. In Fig.(7) the spectral function is plotted for three different ka = 0, π/2, π. In the same graphs, the spectral functions calculated within momentum average approximations with the corrections of level 2 [68] are presented. There is quite good agreement between spectral functions calculated within these approaches at large momenta (ka = π/2, π). At ka = 0 the incoherent part of the spectrum is also quite similar in both cases. Nevertheless, the present theory shows much sharper peaks. These spectral functions satisfy the sum rules up to 7th moment (Fig(6)).

The spectral function in the polaron state
In order to describe the spectral function of the system where the ground state of the grand canonical Hamiltonian H ′ corresponds to a finite polaron density, we have to tune the chemical potential close to the single polaron level. The requirement is that the lowest energy state of the grand canonical Hamiltonian H ′ with at least one polaron should be lower than the lowest energy level of the Hamiltonian without polarons. In order to prevent bipolaron formation, we consider the spinless fermions and assume that bipolarons represent an excited state on the grand canonical Hamiltonian when the chemical potential is pinned near the polaron level.
The spectral function Eq. (11) has two terms. The first term describes the creation of an additional carrier in the ground state of the grand canonical Hamiltonian H ′ with one polaron which is proportional to 1 − n, where n = 1/N is the polaron density. More important is the second term which describes the emission of the electron from the ground state which is proportional to n = 1/N . Importantly, exactly this term is measured in direct photoemission spectroscopy. Note that the second term carries information about polaron formation. This was pointed out in Ref. [50] where the exact formula for the first moment of the spectral function was derived. This first moment has the contribution, which is proportional to the adiabatic polaron shift 2nE p . Therefore, this term may be derived from the adiabatic theory without involving nonadiabatic corrections.
In the adiabatic approximation, the polaron state is the self-trapped state localized in the translation symmetry broken deformation field ϕ n . This state is degenerate because the polaron energy does not depend on the polaron position and is described by the solution of Eq.(5) ψ 0 n with the lowest energy E 0 . Within the sudden approximation the second term in Eq.(11) is proportional to the square of the Fourier transform of the F k = n exp (ikn)ψ 0 n / √ N . Therefore, Eq.(11) is given by a single δ function: Since the polaron wave function is localized F k ∝ 1/N 1/2 the spectral function A − (k, ω) is proportional to polaron density N −1 (one polaron per N sites). The chemical potential µ is pinned to the polaron energy µ = E 0 +E p n |ψ 0 n | 4 [49] which is the sum of the lowest eigenvalue of Eq(5) and the deformation energy caused by the polaron. In the strong coupling limit µ = −E p − t 2 /E p therefore the peak of the spectral function A(k, ω) is at energy ω = −E p + t 2 /E p . Note, that the spectral function is equal to 0 at the chemical potential (ω = 0). This is because the overlap of the wave functions with and without polaron deformation is exponentially small and the annihilation of the polaron together with the deformation is forbidden due to the Franck-Condon principle. The part of the spectral function associated with the inverse photo-emission is given by the first part of Eq.(11) and may be written as where F l k = n exp (ikn)ψ l n / √ N is the Fourier transform of the wave function which represents the l-th eigenstate of Eq.(3) where ϕ n is determined by Eq.(4) with ψ l n = ψ 0 n . Here the sum does not include l = 0 because the absorption of electron directly to the polaronic state is proportional to the overlap of the lattice wave functions with zero and nonzero deformation which is exponentially small ∝ exp (− M/m). As in the photoemission case this process is forbidden due to the Franck-Condon principle. Note that the spectrum has a clear pseudogap, which corresponds to |E 0 | − 2t ≈ 2E p − 2t. The physics of this pseudogap is simple. The probability to remove polaron from the Fermi level or the probability to add polaron to the Fermi level is exponentially suppressed and these processes are too weak to observe them experimentally. Therefore, the main features of the spectral function are shifted from the Fermi energy because of the Franck-Condon principle.
In Fig.(8) the spectral function of the polaron is plotted for the 1D system with 100 sites with periodic boundary conditions in the adiabatic limit for the polaron binding energy E p = 2.5t Fig.(8a) and for E p = 0.5t Fig.(8b). The main spectral features repeat the free electron spectrum except that the presence of polaron deformation breaks the translation invariance. Therefore, the linewidth has finite broadening due to the scattering of a free electron on the polaron deformation. As expected, this spectral density satisfies the zero M 0 = 1 and the first M 1 = −2tcos(ka)−2nE p sum rules, a is the lattice constant and n = 1/N is the polaron density.
In the limit when the polaron size is larger than the lattice constant Eqs. (3,4,5) have analytic solutions for both the localized and itinerant part of the spectrum (See Appendix B). Using these solutions the matrix elements F k and F l k in Eqs. (23,24) may be evaluated analytically: Here a is the lattice constant, N is the number of sites in the system, κ = l − k, where l represents the solutions of Eq.(38) k = 2πn/N a where n = 0, 1, 2, ., N − 1. These two expressions substituted to Eqs. (23,24) with E 0 = −2t − E 2 p /4t, E l = −2t + tl 2 a 2 , and the chemical potential µ = −2t − E 2 p /12t provides the analytic description of the spectral function which is accurate up to the nonadiabatic corrections which are small as m/M . Note that this analytic expression satisfies the exact sum rules (i.e. M 0 and M 1 ) [50] when the polaron density is finite.
The spectral function A(k, ω) calculated on the basis of Eqs. (24,23) using Eqs. (25,26) for two different momenta is plotted in Fig.9. As expected the results are very similar to that plotted in Fig.8. There is a weak spectral intensity at E 0 − µ. This intensity is proportional to the polaron density and is suppressed quickly when momentum moves away from the k = 0 point. Main intensity is centered at ω = E 2 p /12t + tk 2 a 2 . Because the exact wave functions Eq.(37) are not plane waves there is a natural broadening of the spectral line due to the scattering of an electron on the polaron deformation. Therefore, the line, corresponding to an almost free electron has finite width in energy. This broadening is very important because it compensates the contribution from the polaron state at negative energies in the sum rules. Therefore, the exact sum rules M 0 = 1 and M 1 = −2nE p + tk 2 a 2 + E 2 p /12t [50] are satisfied. Note that in Ref. ([50]) the sum rule M 2 for the spectral function was derived (See Eq.(15) in [50]). This moment contains averages of ϕ n and ϕ 2 n , which cannot be evaluated exactly. Nevertheless, within the adiabatic approximation, these averages may be easily evaluated. Indeed, Eq.(4), defines the deformation field at site n. Therefore, integration of this equation with ψ 0 (x) defined by Eq.(34) leads to the result, obtained in [50] ϕ = − √ 2g which corresponds to the case of one polaron per N sites. This immediately leads to the correct equation for the momentum M 1 . Similarly, ϕ 2 = 2g 2 n |ψ 0 n | 4 . Calculation of the sum leads to the following result ϕ 2 = g 2 E p /3t in the weak coupling case E p << t and ϕ 2 = g 2 (1 − t 2 /E 2 p ). Therefore, the second moment of the spectral function has the form: In Fig.10 the calculated sum rule M 2 is plotted as a function of momentum k in comparison with expression Eq. (27). The calculations are performed for the tight binding model (Eq.(3,4)) and for analytic differential equation (35). In both cases the results perfectly match the analytic formula Eq. (27). CONCLUSION We formulate the analytic theory for the spectral function of the electron-phonon system in the adiabatic limit. In the limit of the dilute polaron system where the polaron density n → 0 and when the chemical potential is large and negative µ → −∞ the polaron deformation is absent and the spectrum of electrons coupled to phonons is defined by the standard perturbation theory. Using the diagram technique with an accurate account of threshold effects we were able to formulate an integral equation for the vertex, which takes into account all diagrams up to the sixth order. Then we solve Dyson's equation for the self-energy with the renormalized vertex and calculate the polaron spectral function in 1D with an accuracy better than 3% for λ ≤ 3 in the adiabatic limit. The moments of the spectral function satisfy the exact sum rules up to the 7th moment.
For the system with a finite polaron density when the chemical potential is pinned at the polaron level the ground state has nonzero adiabatic deformation due to the presence of polarons. In this case, the spectral function is calculated in the adiabatic limit without any non-adiabatic corrections. We also show how the adiabatic terms, proportional to the polaron density n, may be evaluated for higher-order moments. The spectral function shows weak spectral intensity at E 0 − µ, which is proportional to n. At the Fermi level in the adiabatic limit the spectral density is absent. Continuous spectrum starts at −µ − 2t > 0. Therefore, the spectrum has a gap, which is proportional to the polaron binding energy E 0 . Contrary to the case of n = 0 the spectral function with a finite polaron density has some adiabatic contributions which are characteristic of the polaron formation.
Eq.(5) in 1D in the weak coupling limit t >> E p when the polaron radius is much larger than the lattice constant has an analytic solution. In that limit Eq.(5) may be rewritten in a differential form. Indeed, expanding ψ n+m in a Taylor series up to the second order in the lattice constant Eq.(5) has the form: − ta 2 d 2 ψ(x) dx 2 − 2E p |ψ(x)| 2 ψ(x) = (E + 2t)ψ(x). (33) It is easy to check that the function is the solution of Eq.(33) corresponding to the bound state energy E 0 = −2t − E 2 p /4t. Therefore, Eq.(3) with deformation field ϕ 0 (x) from Eq.(4) may be rewritten as: Rescaling the spatial variable E p x/2at → x this equation is reduced to the well known equation [69]: This equation is integrable and has analytic eigenfunctions. Therefore, the excited states of Eq.(35) have the following form: with the excited state energies E k = −2t + tk 2 a 2 . These wave functions satisfy the periodic boundary conditions when where L = N a. This equation has N − 1 roots for corresponding to N − 1 itinerant states, because one state splits to localised polaron level E 0 . These values of k together with ψ 0 (x) (Eq.(34)) define the complete set of eigenfunctions of Eq. (35). The Fourier transform of these eigenfunctions defines the matrix elements Eqs. (25,26) which determine the polaron spectral function.