Self-consistently renormalized quasiparticles under the electron-phonon interaction

Asier Eiguren,1 Claudia Ambrosch-Draxl,2 and Pedro M. Echenique1,3 1Donostia International Physics Center (DIPC), Paseo Manuel de Lardizabal 4, 20018 Donostia-San Sebastián, Spain 2Chair of Atomistic Modelling and Design of Materials, University of Leoben, Franz-Josef-Straße 18, A-8700 Leoben, Austria 3Departamento de Física de Materiales, Facultad de Ciencias Químicas, UPV/EHU, Apartado 1072, 20080 San Sebastián, Basque Country, Spain and Centro Mixto CSIC-UPV/EHU, Apartado 1072, 20080 San Sebastián, Basque Country, Spain Received 5 March 2009; published 2 June 2009


I. INTRODUCTION
2][3] Physically, a quasiparticle state is viewed as the original noninteracting state, but dressed-or surrounded-by an "interaction cloud," 4 which is responsible for the renormalization.The actual velocities and lifetimes of quasiparticles have an important impact on thermodynamical and transport properties of real materials. 5,6This is the reason why understanding the renormalization effects caused by electron-phonon interaction has been subject of intense experimental and theoretical activities during the last years.The decisive impulse came from the experimental side, with the advent of high-resolution angle-resolved photoemission spectroscopy ͑ARPES͒. 7,8Recent developments have led to a very high resolution, i.e., in the meV range, which, in turn, brings out effects that pose a challenge for many-body theories.In this line, surfaces and Shockley-type surface states represent privileged systems for probing many-body theories, since surfaces are directly accessible by ARPES.There are several outstanding surface systems where the strong influence of the electron-phonon interaction on quasiparticle dispersion and lifetimes has been clearly demonstrated, [9][10][11][12][13] and some review papers [14][15][16][17][18] covering the subject.
The many-body Green's-function theory 1,19 represents the most appropriate framework for describing and even defining quasiparticles.In this theoretical framework, an ideal quasiparticle appears as an excitation with a perfect Lorentzianshaped spectral function A͑͒ = ⌫ / ͓͑ − ⑀ qp ͒ 2 + ⌫ 2 ͔.In this simplest case, the energy-time uncertainty principle and the energy width ⌫ allow definition of the decay time =1/ ⌫.The position of the peak, ⑀ qp , provides the quasiparticle en-ergy.The quasiparticle picture is then well justified as long as the actual spectral function ͑measured or calculated͒ shows an approximate Lorentzian profile.If only one such peak is present, one can say, in turn, that a singlequasiparticle picture is valid.However, it has been found that electron-phonon interaction introduces more complex structures than simple Lorentzian functions, and that several peaks can be observed, as shown theoretically 6,20,21 as well as experimentally in many surface systems 9,11,12,17,18,22 and in high-T c materials, 8,[23][24][25][26] among others.In all the above examples, the electron spectral functions show ͑at least͒ two peaks for electron states close to typical phonon energies, and make it clear that the single-quasiparticle picture breaks down. 9This was first pointed out by Engelsberg and Schrieffer 20 ͑ES͒ in a pioneering theoretical work investigating simple ideal systems.In particular, ES studied the electron-phonon interaction in the Einstein and Debye models at T = 0 and demonstrated that considering the Dyson equation 1,6,19 in the full complex plane gives rise to two quasiparticle solutions, one behaving similarly to a polaron and the other one behaving as a damped dispersing electron state.As far as the authors know, these ideas have only been considered as a conceptual tool, but without application in real material systems, until very recently. 21he state-of-the-art of theoretical framework for obtaining quasiparticle energies and lifetimes consists of calculating the appropriate Eliashberg function and the subsequent computation of the second-order perturbative electron self-energy ⌺. 10,[27][28][29] The electron lifetimes are directly obtained from the imaginary part of the self-energy, ⌺ I .The real part ⌺ R is commonly considered to be uniquely responsible for the band distortion or renormalization.This is a simplified picture because the effects of real and imaginary parts appear separated, in contrast to the conclusions of the work by ES, who first showed that both should be considered selfconsistently.The need for a self-consistent treatment is understood taking into account that the real part of the electron self-energy is not only a function of real energies, but a func-tion of both the real and imaginary components.As will be shown, the separation of the original complex Dyson equation into imaginary and real components can be achieved only by ignoring the effect of the imaginary energies on the electron self-energy.This would be well justified for electrons with negligible lifetime broadening ͑in comparison to the real part of the self-energy͒, but this is not generally the case, particularly in the electron-phonon problem.The only general method for estimating renormalization corrections is the so-called quasiparticle expansion approach ͑to be introduced later͒, which assumes a weak imaginary part of the self-energy, and demonstrates a very limited range of validity.
In this paper we consider the electron Dyson equation in the full complex plane, similarly as ES did, but including finite temperatures and combining it with ab initio electron structure theory.This allows us to find quasiparticle poles in a realistic and self-consistent way, without neglecting the imaginary part of the self-energy.In particular, we show that the full ͑complex͒ Dyson equation shows several quasiparticle solutions.Thus, the formalism provides a broader quasiparticle interpretation, since several quasiparticles are found in a natural way, without any ad hoc assumption.In this line, we investigate under which circumstances the spectral functions are well represented not within a single-quasiparticle, but rather a multiple-quasiparticle picture, checking if the multiple peaks observed in the spectral functions can be traced back to calculated multiple quasiparticles.
We consider three example systems in order to illustrate the theory.These are the simple Einstein and Debye models, which we treat quite in detail, bulk MgB 2 , and the 1 ϫ 1 hydrogen-covered ͑and deuterium-covered͒ W͑110͒ surface.

II. THEORY
In standard many-body perturbation theory, [1][2][3][4]19 the selfenergy ⌺ k ͑rЈ , r , E͒ and the electron Green's function G k ͑rЈ , r , E͒ are the two key quantities. G  ͑rЈ , r , E͒ describes the propagation of an electron with momentum k from rЈ to r within an interacting many-body system, 1 and ⌺ k ͑rЈ , r , E͒ encodes the information about the scattering processes.Denoting by G 0 the noninteracting Green's function, the geometric sum of the perturbative series ¯formally leads to the well-known Dyson equation, expressing the exact relation between G and ⌺, G͑r,rЈ,E͒ = ͓G 0 ͑rЈ,r,E͒ −1 − ⌺͑r,rЈ,E͔͒ −1  ͑r , rЈ , E͒ were an energy-independent function, the system would become linear, and one would obtain a unique solution for each electron wave vector k.Thus it is the energy dependence of ⌺ ͑the nonlinearity͒ that enables the existence of several quasiparticle states for the same k.
At the moment, we have considered E as a purely real variable.However, it is well known that quasiparticles with mere real energy ͑infinitely long living͒ do not exist in the thermodynamic limit, 3 and thus the corresponding poles must have finite imaginary components.The review by Farid in Ref. 3 treats this subject in detail.In order to obtain the complex poles of the Dyson equation, we must extend the theory into the complex plane.In other words, ⌺ and G need to be considered not only as complex functions, but also as functions of complex energy.However, simply replacing E with its complex counterpart z does not work; i.e., it leads to a Dyson equation without solution.This is easily shown as follows.Replacing E with z in G and ⌺, one obtains that these functions follow the reflection property Thereby G͑z͒ is the analytic continuation of G͑E͒ from the real axis into the entire upper half plane, 1 where G is analytic everywhere for Im͑z͒ Ͼ 0 as implied by causality. 30This is the reason why the branch in the upper ͑lower͒ half plane is commonly referred to as the ͑un͒physical Riemann sheet.This leads us to the contradiction, however, that Eq. ͑3͒ implies poles in both half planes, because if z qp is a pole of G͑z͒, so should ͑z qp ͒ ‫ء‬ be, and causality implies the existence of poles only in the lower half plane ͑if any͒.Thus, one concludes that substituting G͑z͒ and ⌺͑z͒ in the Dyson equation does in general not provide any solution, except for the real axis in some exceptional ͑academic͒ cases. 3From Eq. ͑3͒, it is also clear that ⌺͑z͒ and G͑z͒ are discontinuous across the real axis if The solution to this apparent problem is well known 1,3 and consists of considering the analytic continuation of G from the upper into the lower half complex plane.Thus, in order to find a meaningful Dyson equation, one must proceed as follows.Replace E with z in G for the upper half plane, and analytically continue ͑by whatever method͒ the soobtained G͑z͒ from the upper to the lower half plane.In this way we obtain a Green's function which is analytic for Im͑z͒ Ͼ 0, but at the same time has all its possible poles located in the lower half plane.As an example, one could calculate G͑z͒ for Im͑z͒ Ն 0 by means of standard secondorder perturbation theory, and then consider a Padé approximation for Im͑z͒ Ͻ 0. In the following, we will refer by G ˜͑z͒ the so-obtained analytically continued function, in contrast to G͑z͒ which exhibits reflection symmetry, i.e., G͑z͒ ‫ء‬ = G͑z ‫ء‬ ͒.
The true quasiparticle equation is then Det͉G ˜k 0 ͑r,rЈ,z͒ −1 − ⌺ ˜k͑r,rЈ,z͉͒ = 0, ͑4͒ with G ˜͑z͒ and ⌺ ˜͑z͒ referring to the functions obtained by the above described procedure.Writing Eq. ͑2͒ in the basis of the unperturbed set of electron eigenstates i,k 0 , we obtain where ⌺ ˜i,j,k ͑z͒ϵ͗ j,k 0 ͑rЈ͉͒⌺ ˜k͑r , rЈ , z͉͒ i,k 0 ͑r͒͘.Equation ͑5͒ demonstrates that neglecting the nondiagonal components, we obtain a simpler Dyson equation with the same form as found for homogeneous systems,

͑6͒
In principle, the quasiparticle wave functions ⌽ i,k qp ͑r͒ are also available from Eq. ͑4͒.Calculating ⌽ k qp ͑r͒ in real space would directly visualize the "dressing" of quasiparticles, but as far as the authors know, such calculations have not been performed so far.
In the following, we will only consider the diagonal form of the Dyson equation.The real and imaginary components of Eq. ͑6͒ constitute two coupled nonlinear equations The complex solutions z qp ͑k , n͒ = E R ͑k , n͒ + iE I ͑k , n͒ represent the energies ͓E R ͑k , n͔͒ and lifetimes ͓1 / = E I ͑k , n͔͒ of the quasiparticle states, with the index n labeling them.In principle, Eq. ͑7͒ must be solved self-consistently, and the main problem consists of obtaining the analytically continued self-energy ⌺ ˜͑z͒.Staying on the real axis, one can, of course, avoid any analytic continuation.Hence, as already mentioned, it is customary to simply neglect the effect of the imaginary part of the energies in the argument of the selfenergy: The upper index ͑1͒ indicates that we consider this as a first approximation.We will refer to these equations also as first renormalization.The advantage of Eqs.͑8͒ and ͑9͒ is that they appear to be decoupled, and that the imaginary part of the self-energy, E I 1 , directly provides an approximate inverse lifetime.Unfortunately, there is no justification for the above approximation, since the basic assumption is to ignore E I inside the self-energy, while it is well known that in the electron-phonon problem, real and imaginary parts of ⌺ are generally of the same order.In any case, Eq. ͑8͒ gives the correct quasiparticle dispersion close to the Fermi level and T = 0 because E I → 0 in this limit.At T = 0 and just at the Fermi level, where Eqs.͑8͒ and ͑9͒ are still valid, the effective mass of quasiparticles is enhanced by m = m 0 / ͑1+͒, where m 0 is the bare electron mass, and is a dimensionless parameter. 6reflects the strength of the electron-phonon coupling and is defined as In the low-energy, low-temperature range, the quasiparticle dispersions are renormalized compared to the unperturbed values ⑀ i,k , according to the first renormalization, i.e., E R ͑1͒ = ⑀ i,k / ͑1+͒.In Sec.II A we will make a step beyond the above approximation ͓Eqs.͑8͒ and ͑9͔͒, by including a first correction arising from the imaginary part of the poles.

A. Quasiparticle expansion
The so-called quasiparticle expansion is the most popular approximation for improving over the results obtained by Eqs.͑8͒ and ͑9͒, which completely ignore the effect of the imaginary part of the poles in the argument of the selfenergy.Opposite to the procedure in which one obtains the analytic continuation ⌺ ˜exactly for the lower half complex plane, here the idea is to treat the continuation only approximately.This is done by means of a first-order Taylor expansion of ⌺͑z͒ around z 0 = E R ͑1͒ + i0 + , i.e., by considering ⌺ ˜͑z͒ Ӎ ⌺͑z 0 ͒ + ͑z − z 0 ͒⌺Ј͑z 0 ͒ in the Dyson equation.In this way, E R ͑1͒ , given by Eq. ͑8͒, is considered as a first estimate, and the improved but still approximate Dyson equation reads

͑11͒
The inferred pole z k ͑2͒,qp is straightforwardly equated from the above equation, and its real and imaginary parts are Thereby we have made use of Eqs.͑8͒ and ͑9͒.The denominators in Eqs.͑12͒ and ͑13͒ represent the imaginary and the real parts of the quasiparticle renormalization factor Thus the above equations may be rewritten as Equation ͑15͒ shows that E R 1 ͑k͒ is shifted due to the imaginary parts of both ⌺ and Z qp .The imaginary part of the poles is renormalized by the real part of Z qp .
We will refer to this step as a second renormalization, representing the first correction due to the imaginary part of the poles in the energy argument of the self-energy.Unfortunately, the range of validity of the above equations is really limited, mainly because of the poor quality of the analytic continuation obtained by means of a first-order Taylor expansion.Note, for instance, that the quasiparticle renormalization factor Z qp diverges for ⌺Ј ϳ 1, a limit easily reached in electron-phonon interaction.Therefore the importance of Eq. ͑15͒ is mainly conceptual, as it shows that the real and imaginary parts of the quasiparticle poles are strongly coupled, and that a simple expansion cannot cope with the problem.
In the following sections, we will deal with the solution of the full Dyson equation instead of the simplified Eq. ͑8͒ or ͑15͒, and the energies without index ͑E R and E I ͒ will indicate the complete solution.

B. Quasiparticle spectral function
It is convenient to introduce the spectral function before analyzing the complex Dyson equation.The imaginary part of the electron Green's function gives the electron spectral function as a function of energy,

͑17͒
Let us suppose for a moment that the analytically continued Green's function has only one pole located at z k qp = E R ͑k͒ + iE I ͑k͒ in the lower half of the complex plane.Then the first-order Laurent expansion of G ˜͑z͒ around z qp is given by with the renormalization factor Z k qp ͓Eq.͑14͔͒ evaluated at z k qp .Note that Z k qp is complex in general, and the imaginary part of the expansion gives rise to two terms, where we have defined

͑20͒
If k qp = 0, i.e., if Z k qp is a real number, the ideal quasiparticle spectral function is recovered: The real part of the renormalization factor reflects the total spectral weight of the quasiparticle state, A finite imaginary part of Z k qp ͑ k qp 0͒ plays the role of distorting the ideal Lorentzian shape but keeping the total weight constant ͓see the approximate Eq. ͑15͔͒.If Z k qp were purely imaginary it would not give any net contribution to the spectral function.This is schematically shown in Fig. 1, taking a pole located at z qp =2−i as an example.The left panels refer to the case where Z qp is purely real, while the ones to the right correspond to the case where = / 6.In the bottom and top panels we show the model spectral functions A͑z͒ =−Im͓Z qp / ͑z − z qp ͔͒ / for a part of the complex plane and for the real axis, respectively.The contribution of a pole to G ˜͑z͒ behaves like a dipole for A͑z͒ but rotated by an angle with respect to the y axis.A pole with a purely real Z qp produces a perfectly Lorentzian-shaped spectral function ͑left panel͒, while in case of a finite imaginary component of Z qp ͑right panel͒, i.e., qp 0, the spectra shows a distorted shape, with the peak slightly displaced from the real part of the pole.We will see that even in the simplest Einstein model the quasiparticle renormalization factors acquire finite imaginary components, with the consequence that the spectral peaks do not exactly correspond to the real parts of the quasiparticle poles.As mentioned, nothing excludes the possibility of finding several quasiparticle poles when considering a Dyson equation in the entire complex plane.In this case, one would straightforwardly generalize Eq. ͑19͒, for several quasiparticle poles, and a multiple-quasiparticle picture could be defined as the case in which the expansion approximately coincides with the actual spectral function

C. Analytic continuation and the Einstein model
The Einstein model is the simplest and most fundamental problem in the context of electron-phonon interaction.It was first studied by Engelsberg and Schrieffer at zero temperature, 20 where the analytic continuation was easily obtained due to the special properties of ⌺ ˜͑z͒ at T =0.
The Einstein model describes a linearly dispersing electron state, interacting with a single phonon mode with energy 0 .This is a highly idealized case, but as will be demonstrated below it is very useful even when treating realistic systems with ab initio quality.For this system, the secondorder self-energy in integral form 6 is As will be shown in Appendix A, the analytic continuation of ⌺ Ei ͑E , 0 ͒ from the upper to the lower half complex plane can be expressed in terms of the digamma function ⌿.

͑24͒
Thus, ⌺ Ei ͑z͒ and ⌺ ˜Ei ͑z͒ are identical for the entire upper half plane including the real axis, but for the lower half, ⌺ Ei ͑z͒ is simply the reflection of the upper half ⌺ Ei ͑z ‫ء‬ ͒ = ⌺ Ei ͑z͒ ‫ء‬ , and hence it is a discontinuous function across the real axis.In contrast, ⌺ ˜Ei is analytic everywhere except at isolated simple poles at negative Matsubara frequencies ͑plus Ϯ 0 ͒, z j p = Ϯ 0 − i͑2j +1͒k B T, with j =0, ... ,ϱ.In the following, we set k B = 1 for simplicity.Equation ͑24͒ is the required continuation to be considered in the complex Dyson equation, i.e., Eq. ͑6͒.As mentioned, an alternative could be a Padé approximation ͑or a Taylor expansion͒, but we have found that using Eq.͑24͒ offers various practical advantages.The Dyson equation for complex poles is then and the solutions must be found numerically.A simple "brute force" method would consist of calculating ⌺ ˜Ei ͑z , 0 ͒ for a sufficiently fine grid spanning the lower half of the complex plane and to check whether ͉z − ⑀ k − ⌺ ˜Ei ͑z , 0 ͉͒ Ͻ ␦ for arbitrary small ␦.Instead, we have considered the complex version of the Newton method, i.e., iterating the map until convergence.Note that we need an initial guess for the poles z 0 to start the iterative process.If one were choosing z 0 = E R ͑1͒ ͑k͒, the first iteration of Eq. ͑26͒ would be equivalent to the quasiparticle expansion or the second renormalization.Hence it seems convenient at this point to refer to the above approach as the self-consistent renormalization, in contrast to the first or second renormalization as already described.
In order to find all possible solutions of Eq. ͑25͒ ͑for the same k͒, one must consider several initial values z j=0 , distributed over the relevant part of the complex plane, say, 0 Ͻ Re͑z 0 ͒ Շ 4 0 and 0 Ͼ Im͑z 0 ͒ տ −4 0 .In practice, the basin of attraction for each solution of the map ͓Eq.͑26͔͒ is generally a fractal, and a relatively dense sampling is need.We typically consider about 10 3 of such z 0 points.
ES studied the Einstein model at zero temperature.In this temperature limit, the continued electron self-energy has the following symmetry properties: In other words, for real electron energies smaller than the phonon frequency 0 , the continued self-energy exhibits the reflection property ⌺ ˜Ei ͑z ‫ء‬ ͒ = ⌺ ˜Ei ͑z͒ ‫ء‬ , and for ͉E R ͉ Ն 0 one must consider the jump due to the imaginary part.ES followed these properties in order to obtain the continuation of ⌺, but, unfortunately, this procedure is only valid precisely for T = 0.The advantage of considering Eq. ͑24͒ instead is that it is analytic for all temperatures and contains also the correct T = 0 limit.The left panels of Fig. 2 show the quasiparticle dispersions, i.e., the dependence of the real part of the poles, E R , on the unperturbed energy ⑀ k 0 for the three different temperatures, T =0 ͑top panel͒, T = 0 / 10 ͑middle panel͒, and T = 0 ͑bottom panel͒.The solid ͑black͒ lines represent the solution of the complex Dyson equation ͓Eq.͑25͔͒, while the thin ͑gray͒ line is the solution of approximate equation considering only real quantities ͓Eq.͑8͔͒.The renormalization factor Z qp , evaluated at the self-consistently calculated poles, is depicted in the right panels as a function of E R .The solid thick lines ͑red͒ represent the real part of Z qp , while the thin ͑green͒ line corresponds to the imaginary part Im͑Z qp ͒.In Fig. 3 the imaginary part of the poles, E I , is shown as a function of E R for T =0 ͑top panel͒, T = 0 / 10 ͑middle panel͒, and T = 0 ͑bottom panel͒.The solid ͑black͒ line is the solution of the complex Dyson equation ͓Eq.͑25͔͒, while the dashed ͑red͒ line is the imaginary part of the bare selfenergy, E I ͑1͒ ͓see Eq. ͑9͔͒.In Fig. 2, the quasiparticle dispersions at zero temperature ͑top panel͒ are similar to the ones calculated by ES, 20 except for the case ͉⑀ k 0 ͉ Ͻ 0 , where ES found two quasiparticle states, while we obtain only one.However, the expansion of Eq. ͑24͒ around z = 0 and T =0, ⌺ ˜Ei ͑z͒ϳz, shows indeed that only one quasiparticle state exists for E R → 0, namely, As in the work by ES, for ⑀ k 0 Ͼ 0 we also find two solutions of the Dyson equation.The one at E R Շ 0 behaves similarly to a polaron state, practically without dispersion.This electron state is localized due to the heavy phonon cloud originating from the virtual-phonon processes, which have maximum amplitude at E R ϳ 0 .At zero temperature, the lifetime of this state is infinite ͑E I =0͒, explaining why it coincides with the approximate solution ͓Eq.͑8͒; Fig. 2, thin gray line͔.The other state at E R Ͼ 0 is not restricted to virtualphonon processes only, but is energetically allowed to emit also real phonons.Hence it acquires a finite lifetime.This is appreciated in Fig. 3.At T = 0, the self-consistent imaginary part of the self-energy, E I , ͑solid line͒ follows the same qualitative behavior as the bare one ͑dashed line͒, but is significantly enhanced for electron energies just above the phonon frequency ͑E R տ 0 ͒.This means that the electron state for E R Ͼ 0 is damped at an even higher rate than predicted by the bare self-energy ͉⌺ I ͉ ͓Eq.͑9͔͒, where the former is basically equivalent to Fermi's golden rule.The enhancement of E I with respect to ⌺ I is due to band renormalization, i.e., due to the deformation of the quasiparticle bands by electron-phonon interaction.This behavior was anticipated by the quasiparticle expansion approach in a qualitative way, where the renormalization effect was partially accounted for by considering a Taylor expansion E I ͑2͒ = E I ͑1͒ Re͓Z qp ͑E R ͑1͒ ͔͒.However, Z qp calculated using Eq.͑8͒ ͑not shown͒ diverges for certain energies E ͑2͒ , showing an unphysical behavior.Note that this problem is solved when considering the selfconsistent solution of Eq. ͑25͒.At T = 0 / 10, the polaron kind of state remains practically unchanged compared to the T = 0 result, the main difference being that several additional states appear close to E R Ӎ 0 .The middle panel of Fig. 2 ͑right͒ demonstrates that the spectral weight of most of these states is negligible since Re͑Z qp ͒ is relatively small.However, the states with ͉E I ͉ Շ 2 0 are still physically meaningful.For example, the damped electron-state shearing part of the E R Ӎ 0 range shows a reasonable behavior with the renormalization factor developing continuously from zero to unity for high energies.Since we have considered a coupling constant given by =−‫ץ‬⌺ R ͑E͒ / ‫ץ‬E ͉ E=0 = 1, the real part of the renormalization factor at the Fermi level ͑E R =0͒ is Re͑Z qp ͒Ӎ1 / ͑1+͒ Ӎ 1 / 2 at E = 0, for both T = 0 and T = 0 / 10.Note that Re͑Z qp ͒ vanishes continuously when going to E R → 0 .At higher temperatures, such as T = 0 , the same qualitative behavior is observed, but Re͑Z qp ͒ is significantly weakened for E R Ͻ 0 .For all temperatures and E R ӷ 0 the damped electron dominates the full spectral weight since Re͑Z qp ͒ → 1 for high energies.
Figures 4 and 5 show the electron spectral functions for ⑀ k = 2 at temperatures T = 0 and T = 0 / 10, respectively.The bottom panels represent the analytically continued spectral function A k ͑z͒ =−Im͓G ˜k͑z͔͒ / for ⑀ k =2 0 .The top panels demonstrate A k ͑E R ͒ in the real axis.At T = 0, two poles are clearly visible ͑bottom panels͒, one located at the real axis at E R ϳ 0.8 0 and the other one at ͑E R ϳ 1.8 0 , E I ϳ −2 0 ͒.The first one corresponds to the polaron state, while the other one is a damped electron.At finite temperatures, these two states are again discernible, but several additional states show up with smaller spectral weight.

D. Momentum-dependent self-energy in real materials
Besides the electron and phonon band structures, a key ingredient for calculating electron-phonon interaction related properties is the Eliashberg function 6 This function is basically the phonon density of states weighted by the electron-phonon matrix elements, and gives the probability of an electron-phonon scattering event transferring energy at T =0.i and j label the different electron bands, k and q are the electron and phonon wave vectors, and q, and g q, i,j stand for the phonon energy and the electron-phonon matrix elements ͑related to the phonon mode ͒, respectively.
The Eliashberg function corresponding to the Einstein model with phonon energy 0 and matrix element g, is basically a Dirac delta function.Thus Eq. ͑31͒ may be reinterpreted as the superposition of effective Einstein modes with energies and ␣ 2 F i,k ͑͒ playing the role of the interaction strength.In this way, the total second-order selfenergy for an electron with band index i and momentum k may be written as a sum of contributions of effective Einstein modes.
⌺ ˜i,k ͑z͒ is analytic across the real axis because so is ⌺ ˜Ei ͑z , ͒.
Hence one can use Eq.͑33͒ in the complex Dyson equation to describe any arbitrary system.

III. IMPLEMENTATION
Equation ͑33͒ gives the analytic continuation of the retarded electron self-energy from the upper into the lower half complex plane, i.e., from the physical into the unphysical Riemann sheet.Since analyticity implies that real and imaginary parts of the self-energy fulfill the Cauchy-Riemann relations separately, these functions are harmonic: As a consequence, the problem of obtaining the analytic continuation is equivalent to solving the Laplace equation for ⌺ ˜R and ⌺ ˜I.The boundary conditions would be given by the function values at the real axis together with the normal derivatives ͑with respect to the real axis͒.In addition, it must be imposed that these functions are bounded in the upper half plane due to causality.These conditions read This is a Cauchy boundary-value problem and, in principle, ⌺ ˜R͑E , ⌫͒ and ⌺ ˜I͑E , ⌫͒ could be obtained if ⌺ ˜R and ⌺ ˜I were known in the real axis.Thus, Eq. ͑33͒ formally generalizes the solution given by the Poisson kernel 31 for the upper half plane.The practical difficulty of numerically evaluating Eq. ͑33͒ resides in that the pole contributions of ⌺ ˜Ei ͑z , ͒ at z p =−i͑2n +1͒T Ϯ , with n =0,1,..., induce a difference between the integral of the analytical continuation and the analytical continuation of the integral. 19This problem would be unimportant if the integral of Eq. ͑33͒ were known analytically, which does not occur except for some ideal systems such as the Debye model.In any case, an accurate numerical approximation can be achieved by considering the analytic integral of the piecewise cubic ͑or higher-order͒ interpolation of the ␣ 2 F͑͒ function.An equivalent procedure is to twice, or more times, integrate by parts Eq.͑33͒ in such a way that the continuity of the second derivative of the Eliashberg function is imposed explicitly: The subindex ͑n͒ in Eq. ͑38͒ denotes the nth integral of the Einstein self-energy, i.e.,

⌺ ˜͑n+1͒
Ei ͑z,͒ ϵ ͵ 0 dЈ⌺ ˜͑n͒ Ei ͑z,Ј͒.͑39͒ Opposite to ⌺ ˜Ei ͑z , ͒, it is verified that ⌺ ˜͑n͒ Ei ͑z , ͒ for n Ն 2 are bounded for z p =−i͑2j +1͒T Ϯ , with j =0,1,....In this way, the numerical evaluation of ⌺ ˜i,k ͑z͒ using Eq.͑38͒ becomes plausible.The explicit formulas for the functions ⌺ ˜͑n͒ Ei are given in Appendix C. Appendix B provides the technical details regarding the numerical treatment of Eq. ͑38͒ in terms of ⌺ ˜͑n͒ Ei .In the next subsection, where we will discuss the Debye model, we show that if the Eliashberg function ͑or its first or second derivative͒ contains a discontinuity, it must, however, be handled explicitly.

Debye model
The Debye model describes the interaction of an electron state with acoustic phonons in a solid.The corresponding Eliashberg function is defined to be simply proportional to the phonon density of states, is the Heaviside function.
Proceeding to integrate Eq. ͑33͒ by parts, the first derivative of ␣ D 2 F͑͒ is given by In this way, and repeating once more the same integration procedure one finally obtains Equation ͑43͒ gives the temperature-dependent electron self-energy for the Debye model compactly, and represents the analytic continuation from the upper into the lower half complex plane.ES published a result ͑Eq.5.5 of Ref. 20͒ for the analytically continued ⌺ for the Debye model at T =0, but a careful check shows that it does not fulfill the harmonicity requirement.In our formalism, the T = 0 limit is reached by considering the first term of the expansion of Eq. ͑43͒ around T =0, ͑44͒ Thus, Eq. ͑44͒ revises the result of ES for this case.
For the sake of completeness, we have solved the temperature-dependent complex Dyson equation.The conclusions are very similar to the ones found for the Einstein model; hence the results are not described in detail here.In any case, the left panels of Fig. 6 show the quasiparticle dispersions and the renormalization factors for T =0 ͑top panel͒, T = D / 10 ͑middle panel͒, and T = D / 2 ͑bottom panel͒.The right panels show the complex renormalization factors for the same temperatures.
Figure 7 shows the imaginary part of the poles, E I , as a function of E R for T =0 ͑top͒, T = D / 10 ͑middle panel͒, and T = D / 2 ͑bottom panel͒.The solid ͑black͒ line is the solution of the complex Dyson equation, Eq. ͑25͒, while the dashed ͑red͒ line is the imaginary part of the bare selfenergy, E I ͑1͒ ͓see Eq. ͑9͔͒.These figures demonstrate the existence of a polaronic state also in the Debye model, and that the lifetimes are enhanced for quasiparticles below D -approximately by a factor of 2-while suffering a reduction for states above D .

IV. APPLICATIONS A. 1 Ã 1 H/W(110) and D/W(110) surfaces
Probably, the strongest experimental evidence for quasiparticle band splitting comes from the hydrogen-and deuterium-covered 1 ϫ 1 W͑110͒ surface ͓1 ϫ 1 H/W͑110͒ and 1 ϫ 1 D/W͑110͒, respectively͔.For one of the surface states the angle-resolved photoemission ͑ARPES͒ experiment 12 demonstrated the unambiguous appearance of at least two features in the ARPES spectra, where only one peak should be expected in the absence of electron-phonon interaction.The hydrogen and deuterium coverages naturally have the same unperturbed electron structure, but different vibrational frequencies and amplitudes.In this way, this system provided a textbook example of electron-phonon interaction.
Another interesting feature of this surface is that the detected surface states, labeled commonly as S 1 and S 2 , are spin split and circularly spin polarized with respect to the high-symmetry point S ¯.This behavior could be traced back to spin-orbit interaction, and relativistic ab initio calculations 21,32 have shown very good agreement with the ARPES Fermi-surface and band dispersions. 33The spin polarization introduces other interesting effects such as suppressing or enhancing the electron-phonon matrix elements depending on the angle between the initial and final electron spins. 34For the spin-polarized ARPES and relativistic electronic structure calculations, we refer to Refs. 12 and 33 as well as Refs.21 and 32.In this paper, we include all this information only implicitly, i.e., built into the Eliashberg function.Thus, here we concentrate on the quasiparticle spectra, without worrying about the precise role of the spin.For more details on the calculated vibrational structure, we point to Refs.21  Some results for the 1 ϫ 1 H/W͑110͒ surface were published recently by two of the authors, 21 where we resolved the angular dependence of the electron self-energy and the quasiparticle spectra.Below we show the results for both the hydrogen and deuterium coverages, including additional complementary details.As in Ref.21, we focus on the electron-phonon interaction for the surface state S 1 , for which ample measurements exist. 12The Fermi surface of this system is constituted by two spin-split surface states S 1 and S 2 , as described in Refs.21 and 32.These states originate from a band which is spin degenerate in absence of spinorbit coupling.Within a window of ϳ150 meV, corresponding to typical phonon energies in this system, these bands are well approximated by a linear dispersion.
The vibrational spectra presented in Fig. 8 for the two systems were calculated within the linear-response approach. 39The lower phonon bands ͑Ͻ30 meV͒ are almost not modified going from H to D coverage, because these modes basically originate from the W͑110͒ substrate.In contrast, the three vibrational modes higher in energy, which are related to vibrations of H and D, exhibit an almost perfect scaling by a factor of ͱ 2. These hydrogen-derived ͑deuterium-derived͒ modes are commonly labeled 36 as the wagging mode ͓ϳ90͑60͒ meV͔, and the asymmetric ͓ϳ110͑75͒ meV͔ and symmetric ͓ϳ160͑115͒ meV͔ stretch modes.The connection between the ab initio calculations and self-consistent treatment of the quasiparticle poles is made through the first-principles Eliashberg functions, as described in Sec.II D. In Fig. 9 we show the corresponding calculated Eliashberg functions ␣ 2 F k ͑͒, of the S 1 state in ⌫ ¯S ¯direction for the hydrogen ͑solid black line͒ and deuterium ͑dotted red line͒ coverages.According to the phonon band structure, the lower-energy range of the Eliashberg function is dominated by the contributions from the tungsten modes, while the three structures at higher energies are directly connected to the hydrogen ͑deuterium͒ displacements.To obtain the Eliashberg functions, we used the method presented in Ref. 40 considering a 40ϫ 40ϫ 1 grid for both the electron and phonon wave vectors, where the electron-phonon matrix elements g q,k, = ͚ , Ј ͗ Ј ,q+k ͉␦V q, ͉ ,k ͘␦ , Ј were calculated considering the full spinor states.
For both systems, we solved Eq. ͑6͒ in ⌫ ¯S ¯direction for the S 1 surface electron.The top and bottom panels of tively low temperatures, since for the hydrogen coverage, the maximum phonon energy is m ϳ 160 meV;.Thus, T 1 ϳ 0.02 m and T 2 ϳ 0.08 m .The dashed ͑orange͒ lines indicate the calculated dispersions at T 1 using the approximation given by Eq. ͑8͒.
The lowest-energy splitting at ϳ20 meV roughly corresponds to the main peak of ␣ 2 F͑͒ observed in the lowenergy range in Fig. 10, and is common to both surfaces.As in the Einstein model, a band splitting reflects a boundary between two different interaction regimes.For example, the electron above ϳ20 meV is energetically allowed to emit all tungsten-related phonon modes ͑ Շ 20 meV͒ but below this energy, the virtual excitations become the protagonists.As a consequence, the effective mass of electrons at E Շ 20 meV gets considerably enhanced and reduced for larger energies.This picture is similarly true for the rest of the main structures of the Eliashberg function, and is the origin of the predicted multiple band splittings.The next phonon modes are located at ϳ90 meV for hydrogen, and at ϳ65 meV for deuterium on top of the tungsten surface.For electron states between ϳ20 meV and those energies ͓ϳ90 meV for 1 ϫ 1 H/W͑110͒ and ϳ65 meV for 1 ϫ 1 D/W͑110͔͒, there is an interplay between the real emission processes of the lowest-energy phonon modes and the virtual emission processes of the highest-energy phonon modes.A minor additional splitting is also predicted at ϳ160 meV for 1 ϫ 1 H/W͑110͒ and at ϳ120 meV for 1 ϫ 1 D/W͑110͒.From these results it is clear that the origin of the band splitting observed by Rotenberg and co-workers 12 is a consequence of two effects, i.e., the full dressing of the electron states and the appearance of multiple quasiparticles.In contrast, the approximate solution ͓Eq.͑8͒; Fig. 10, dashed line͔ provides only one quasiparticle state for all wave vectors.
The imaginary part of the quasiparticle poles is shown in the right panels of Fig. 10.The gray and black lines ͑solid and dashed͒ correspond to the calculations for T 1 and T 2 , respectively.The solid lines indicate the fully renormalized inverse lifetimes ͓Eq.͑6͔͒, while the dashed lines stand for the bare imaginary part of the self-energy, given by the approximate expression ͓Eq.͑9͔͒.We should keep in mind that the former approximation is equivalent to the well-known Fermi's golden rule.As an example, at T 1 the inverse lifetime of the lowest-energy quasiparticle state in the left panel of Fig. 10 monotonically increases until a value of ͉E I ͉ ϳ 15 meV is achieved at E R ϳ 20 meV.But there are some energy ranges for which E I decreases for increasing energy.Note, for instance, that the quasiparticles at very low energies live longer than predicted by the bare imaginary selfenergy E I ͑1͒ .For electron energies close to the main peaks of ␣ 2 F͑͒, a similar trend is found.The electron lifetimes are considerably enhanced just below the most important peaks of the Eliashberg function and reduced above.For the damped electron state above all phonon energies, we find that ͉E I ͉ ͑solid line͒ approaches ͉E I ͑1͒ ͉ asymptotically, but maintaining ͉E I ͉ Ͼ⌺ I ͑E R ͒ as in the Einstein model.
In Fig. 11 we show the calculated renormalization factor Z qp as a function of the real part of the quasiparticle poles.The top ͑bottom͒ panel corresponds to the hydrogen ͑deute-rium͒ coverage.The thick ͑red͒ line indicates the real part of Z qp , while the thin ͑green͒ line is the imaginary part.This figure is better understood in comparison with Fig. 10, since the relative spectral weight of each quasiparticle band in Fig. 10 is given by the real part of Z qp .In the case of the hydrogen ͑deuterium͒ coverage, we detect four ͑three͒ main physically meaningful quasiparticle bands with a non-negligible spectral weight.As mentioned, the imaginary part of Z qp indicates the distortion of the quasiparticle contributions to the spectral functions.Thus the calculated poles are not in exact correspondence with the peaks in the spectral functions. 21ote, however, that ͉Im͑Z qp ͉͒ is generally less than ϳ0.2, and we should expect a close correspondence between the real part of the poles and the spectral features.
The top and bottom panels of Fig. 12 depict the calculated spectral functions at T 1 =40 K ͑thick black line͒ and T 2 = 150 K ͑thin gray line͒, respectively.The left panels show the full calculated spectral functions A͑E R ͒ =−Im͓G͑E R ͔͒ / , while the right panels display the same quantity, but calculated out of the quasiparticle contributions alone, as given by Eq. ͑22͒.An ideal quasiparticle state would give one perfect Lorentzian contribution to the spectral function, but Fig. 12 demonstrates that at least two main features exist for several ranges of the electron wave vector.It is remarkable that the function A QP ͑E R ͒ ͑right͒, being constructed from the contributions of the QP states alone, reproduces so well the full spectra ͑left͒.Thus this figure demonstrates that even when a single-quasiparticle picture breaks down, the quasiparticle picture itself could still be valid, but with the extension to several quasiparticles.

B. Bulk MgB 2
A prototypical class of materials where electron-phonon coupling plays an important role are that of natural superconductors with phonon-mediated Cooper pairing.Such a material with a particularly high critical temperature of 40 K is MgB 2 .Hence in this subsection, we have applied our method to bulk MgB 2 , predicting similar behavior as found above for the surface system.The left panel of Fig. 13 shows the Fermi surface.It is composed of two concentric "cylindroids" in the inner part of the Brillouin zone.The Eliashberg function depicted in the right panel corresponds to the outermost cylindroid along the high-symmetry direction AL, as indicated in the figure.We refer to Ref. 40 for more details regarding the computational details as well as the electron and vibrational structures.
The left panels of Fig. 14 show the calculated quasiparticle dispersions at T 1 =0 ͑top panel͒ and at T 2 = 8 meV ϳ m / 10 ͑bottom panel͒.In these panels, the real part of the poles is depicted as a function of the unperturbed electron energy ͑red line͒, and the black envelope indicates the real part of the renormalization factor, Re͑Z qp ͒.Thus the "width" of each band highlights the relative importance of the state regarding its contribution to the spectral function.The gray line represents the solution of the approximate Dyson equation considering only real quantities.Note that the selfconsistent calculation predicts a splitting of the quasiparticle bands at ⑀ k 0 ϳ 100 meV, while in the solution of Eq. ͑8͒ one must go to ⑀ k 0 ϳ 150 meV to find a strong deviation from the bare dispersion.Again, this is due to the imaginary part of the poles that are self-consistently accounted for in Eq. ͑6͒ but not in Eq. ͑8͒.At T 1 = 0 the fully renormalized bands are qualitatively similar to the ones found in the Einstein model, and the splitting occurs approximately at the maximum of the phonon energies.At T 2 = 8 meV ͑lower right panel͒ the polaron and the damped electron bands are detected again, but even more quasiparticle bands appear at typical phonon energies.Most of these bands have a relatively weak quasiparticle renormalization factor although it is not negligible for some of the bands.In the right panels, we show the imaginary part of the poles ͑x axis͒ as a function of the real part ͑y axis͒, again at

V. CONCLUSIONS
In conclusion, we have investigated the Dyson equation of the electron-phonon problem, considering the effect of the imaginary part of the poles self-consistently.The approach proposed by Engelsberg and Schrieffer 20 several decades ago has been extended in such a way that calculations for real materials are possible, and the effect of the temperature can be studied.We compare this theory with the most popular approximations for the Dyson equation in which the imaginary parts of the poles are neglected or considered partially through a quasiparticle expansion.The formalism demonstrates that the former two approximations can be understood in a more general context considering the self-consistent solution of the Dyson equation.In order to illustrate the presented theory, we have investigated the Einstein model, and the Debye model, as well as the spin-polarized 1 ϫ 1 H͑D͒/ W͑110͒ surface and bulk MgB 2 at several temperatures.

ACKNOWLEDGMENTS
A.E. appreciates encouraging discussions with Eugene Sherman.This work was supported by Gipuzkoako Foru Aldundia and Donostia International Physics Center ͑DIPC͒ as well as by the Austrian Science Fund, Project No. P16227.

APPENDIX A: RELATION BETWEEN THE SELF-ENERGY AND THE POLYGAMMA FUNCTIONS
In this appendix we transform the second-order perturbative expression for the electron self-energy into an equivalent form in terms of the polygamma functions ͑see

͑A1͒
If we replace E with the complex-valued quantity z in the integrand of Eq. ͑A1͒, one generally obtains a discontinuous function across the real axis.The objective is to find the analytic continuation from the upper into lower half of the complex plane, i.e., a complex function I ˜͑z , ͒ which is the same as I͑z , ͒ for the upper half of the complex plane but analytic across the real axis.At the moment, we assume that E is above the real axis but infinitesimally close to it, E ϵ E + i0 + .First, we split the integral into two pieces, I͑E,͒ = I 1 ͑E,͒ + I 2 ͑E,͒, ͑A2͒ where As the expansion of the digamma function reads we obtain The reflection property of the digamma function allows us to further simplify the expression and finally we obtain 2T ͪ.

͑A11͒
We now have evaluated Eq. ͑A1͒ for E with an infinitesimal positive imaginary component.As opposite to Eq. ͑A1͒, substituting E with complex z in Eq. ͑A11͒ gives now a function which is analytical across the real axis.This is so because the integral has been already accomplished analytically, ͑A12͒ Thus, Eq. ͑A12͒ is the required analytic continuation of Eq. ͑A1͒ from the upper to the lower half of the complex plane.

FIG. 1 .
FIG. 1. ͑Color online͒ Schematic representations of the analytically continued Green's function G ˜͑z͒ and the corresponding spectral function.We consider the two examples where Z qp is purely real ͑left panels͒, i.e., qp = 0, and where qp = / 6 ͑right panels͒.The quasiparticle pole is located at z qp =2−i in both cases.The two lower panels show a contour plot of the spectral function A i,k L ͑z͒ ϳ Im G ˜͑z͒, with the dashed lines indicating the position of the real axis.The upper panels exhibit the corresponding spectral functions evaluated at the real axis.

FIG. 2 .
FIG. 2. ͑Color online͒ Left: QP band structure for the Einstein model for the coupling parameter = ͉g͉ 2 / 0 = 1 at temperatures T Ӎ 0 ͑top panel͒, T = 0 / 10 ͑middle panel͒, and T = 0 / 10 ͑bottom panel͒.The solid thick ͑black͒ lines represent the real part of the poles, E R , obtained from the complex Dyson equation ͓Eq.͑6͔͒.The solid thin ͑gray͒ line is the solution of the approximate Dyson equation considering only real quantities ͓Eq.͑8͔͒.Right: the quasiparticle renormalization factor Z qp as a function of E R .The solid thick ͑red͒ line shows the real part of Z qp , while the thin ͑green͒ line is the imaginary part Im͑Z qp ͒.

FIG. 4 .
FIG. 4. ͑Color online͒ Representations of the analytically continued spectral function A k ͑z͒ =−Im͓G ˜k͑z͔͒ / for ⑀ k =2 0 and T = 0.The bottom panel shows A k ͑z͒ in the complex plane, and the top panel demonstrates A k ͑E R ͒ in the real axis.

FIG. 6 . 2 D 2 ͑
FIG. 6. ͑Color online͒ Left: QP band structures for the Debye model ͓␣ 2 F͑͒ = 2 D 2 ͑ D − ͔͒ for the coupling parameter =1 at temperatures T Ӎ 0 ͑top panel͒, T = D / 10 ͑middle panel͒, and T = D / 2 ͑bottom panel͒.The solid thick ͑black͒ lines represent the real part of the poles, E R , obtained from the complex Dyson equation ͓Eq.͑6͔͒.The solid thin ͑gray͒ line is the solution of the approximate Dyson equation considering only real quantities ͓Eq.͑8͔͒.Right: the quasiparticle renormalization factor Z qp as a function of E R .The solid thick line ͑red͒ shows the real part of Z qp , while the thin ͑green͒ line is the imaginary part Im͑Z qp ͒.

FIG. 8 .
FIG. 7. ͑Color online͒ Imaginary part of the quasiparticle poles, E I , as a function of the real part E R for different temperatures, i.e., T =0 ͑top panel͒, T = D / 10 ͑middle panel͒, and T = D / 2 ͑bottom panel͒.The thick solid lines ͑black͒ are the self-consistent solutions obtained by means of Eq. ͑6͒, while the dashed ͑red͒ line is the imaginary part of the self-energy E I ͑1͒ = ⌺ I ͑E R + i0 0 ͒ ͓see Eq. ͑9͔͒.

FIG. 10 .FIG. 11 .
FIG. 10. ͑Color online͒ Calculated quasiparticle bands and inverse lifetimes for the hydrogen-͑top panels͒ and deuteriumcovered ͑bottom panels͒ W͑110͒ surfaces.Left: fully renormalized electron bands ͑solid lines͒, calculated from the complex Dyson equation ͓Eq.͑6͔͒ at T 1 = 3.5 meVӍ 40 K ͑gray line͒ and T 2 = 12.93 meVӍ 150 K ͑black line͒.The thin dashed ͑orange͒ line is the solution of the approximate Dyson equation at T 1 considering only real quantities ͓Eq.͑8͔͒.Right: imaginary part of the QP energies ͉E I ͉ =1/ ͑black line͒ versus the real energy of the QP poles, E R .The imaginary part of the bare electron self-energy ͑dashed line͒ is shown for comparison.
FIG. 12. Spectral functions A k ͑͒ at T =40 K ͑thick black line͒ and T = 150 K ͑thin gray line͒ of 1 ϫ 1 H/W͑110͒ ͑top panels͒ and 1 ϫ 1 D/W͑110͒ ͑bottom panels͒ for different momenta compared to their counterparts in the multiple-quasiparticle approximation A k QP ͑͒ as given in Eq. ͑22͒.

FIG. 14 .
FIG. 14. ͑Color online͒ Calculated quasiparticle bands and inverse lifetimes of MgB 2 in direction AL ͑ 2 band͒.Left: fully renormalized electron bands, calculated from the complex Dyson equation ͓Eq.͑6͔͒ at T =8 meVϳ 92 K ϳ m / 10.The black envelope indicates the real part of the renormalization factor and the thin ͑red͒ line inside stands for the calculated quasiparticle pole.The thin ͑gray͒ line is the solution of the approximate Dyson equation considering only real quantities ͓Eq.͑8͔͒ at T = 0. Right: the fully renormalized imaginary part of the quasiparticle energies is indicated by the solid ͑black͒ line, while the dashed ͑red͒ line stands for the bare imaginary part.
0 and ͉g͉ as the electron-phonon matrix element.Each of the two terms inside the large parenthesis contains two different scattering processes.One of them describes the scattering from an initial state with energy E into another state with energy EЈ through the emission of a single phonon with energy 0 .It is weighted by the statistical factor ͓1− f͑EЈ͔͓͒1+n͑ 0 ͔͒.The other process is of truly manybody nature taking into account how the occupied states at EЈ are prevented to scatter by absorption into the ͑now oc-cupied͒ state at E, because of the Pauli principle.This gives a factor of f͑EЈ͒n͑ 0 ͒ and the summation of both gives n͑ 0 ͒ +1− f͑EЈ͒.The second term is similarly obtained by considering the absorption processes of state E and the forbidden emission processes for the occupied states with energy EЈ, leading to the weight ͓1− f͑EЈ͔͒n͑ 0