Padé resummation of many-body perturbation theories

In a typical scenario the diagrammatic many-body perturbation theory generates asymptotic series. Despite non-convergence, the asymptotic expansions are useful when truncated to a finite number of terms. This is the reason for the popularity of leading-order methods such as the GW approximation in condensed matter, molecular and atomic physics. Appropriate truncation order required for the accurate description of strongly correlated materials is, however, not known a priori. Here an efficient method based on the Padé approximation is introduced for the regularization of perturbative series allowing to perform higher-order self-consistent calculations and to make quantitative predictions on the convergence of many-body perturbation theories. The theory is extended towards excited states where the Wick theorem is not directly applicable. Focusing on the plasmon-assisted photoemission from graphene, we treat diagrammatically electrons coupled to the excited state plasmons and predict new spectral features that can be observed in the time-resolved measurements.

Introduction of the Green's function methods to electronic structure calculations is the most prominent achievement of the field-theoretic methods 1-3 on par with the density functional theory having immediate technological applications 4,5 . Even in the lowest (beyond the mean field) order one obtains significant improvements of e. g. the band gap through the correlation shifts (Δ). Including higher-order diagrams (vertex corrections) is numerically demanding and so far the truncation of perturbative expansions has been done in ad hoc manner. Having a tool to systematically perform higher-order self-consistent (sc) calculations would allow to make calculations for some representative systems and extrapolate these results in order to make quantitative statements on the convergence and accuracy of perturbative expansions for specific cases. However, there are fundamental obstacles on the way that arise from dealing with diverging series as the following consideration illustrates. . Nonetheless, a sensible spectral function, = π A z g z ( ) Im ( ) 1 , in the domain of interest can be reconstructed by using the Padé approximation. The procedure is outlined at Fig. 1 where the original function g model (z), the series expansion (1) and the Padé reconstruction are shown. The Padé approximation (PA) allows to obtain very accurate values also in the domain where the series (1) is divergent. The method works so well here because it is known in advance that GF consists of one pole only and this fact is used for the reconstruction: according to the exact form of g model (z) we use the [0/1] approximant (the Padé approximation has form of rational function denoted as [M/N]) with M + N + 1 coefficients, M, N are the orders of the numerator and denominator, respectively 6 ). For realistic calculations we do not have this knowledge and have to rely on some additional assumptions about the analytic structure of the Green's function. As an illustration let us consider the electron-boson Hamiltonian -a model which is ubiquitous in condensed matter physics.

Padé approximation. Let
where c is the creation operator of the deep hole with energy , b † is the bosonic creation operator of the plasmon with the energy Ω. The generalization to the case of multiple fermion kinds as in the mixed-valence impurity model or the plasmon dispersion is possible and will be commented on after the presentation of the diagrammatic solution. The Hamiltonian (3) is quite versatile and is applicable to other scenarios such as resonant-tunneling through a single level coupled to wide-band phonons 16 . Remarkably, also the two particle GF can be found analytically 17 ; the model can be solved at finite temperatures, and some its non-equilibrium properties have been studied 18,19 . Let us consider the following Green's function where |ψ〉 is the exact ground state of the no-hole system with = . It can be diagrammatically found by writing the cumulant expansion for the Green's function () in terms of non-interacting GF g (0) (t) and the cumulant function C(t). Observing that only a single diagram contributes to the cumulant function (in a more general scenario it fulfills an integral equation 20,21 ) one writes The corresponding exact spectral function is compared with the zeroth order and the self-consistent GW (sc-GW) approximations in Fig. 2. The latter is computed using the first diagram in Fig. 3(b); illuminating discussion of the GW approximation in electron-electron vs. electron-boson cases can be found in ref. 18. The results using the Padé approximation (PA). Parameters are as follows: Δ = 1.1,  = 1, η = 0.5. The series expansion (1) is restricted at n max = 10, and PA is applied at the point z = 6. Notice that the original (magenta) and reconstructed (white) densities of states are practically indistinguishable.
Scientific RepoRts | 7: 504 | DOI:10.1038/s41598-017-00355-w are plotted for the strongly correlated regime (γ > Ω) and can be characterized as follows: (i) The quasiparticle (QP) peak is shifted by the energy ω ∆ = γ Ω 2 compared to the noninteracting case; (ii) this main peak is followed by the ladder of plasmonic satellites; (iii) the self-consistent GW method predicts the satellites. However, the position of even the main peak is wrong. This inaccuracy is the main motivation for performing higher-order diagrammatic calculations. Diagrammatic properties. Because the ground state is a no-hole state, ψ † c vanishes and, hence, the time-ordered Green's function only consists of the hole propagator, i. e., can be expressed solely in terms of the lesser GF component, and the corresponding non-interacting GF takes a form: . This fact simplifies the diagrams considerably: (i) in the expansion of the Green's function (g) and the self-energy (Σ) all intermediate points are time-ordered ( Fig. 3(a)); (ii) diagrams containing loops necessarily yield a zero contribution (this is not the case for nonequilibrium states where renormalization of bosonic propagators by fermionic loops needs to be additionally considered 18,22 ). These properties allow to write the self-energy (SE) for this model in analytic form. Because there is no spatial degrees of freedom, the problem is similar to that of the Feynman diagrams enumeration which can be solved by collapsing the space-time variables to one point (the zero-dimensional model [23][24][25][26]. Let ) be an n th-order self-energy term corresponding to a particular diagram, which will be denoted as v. We will prove below that the corresponding expression in the frequency representation is given by the product: . Bosonic propagators are denoted as wavylines. In the frequency representation the given SE diagram yields ω γ ω ω ω 2 ) ( ) . (b) Four lowest orders of the SE diagrammatic expansion of the electron-boson model. Notice that two diagrams of the third-order containing loops are not shown because they are equal to zero. The fourth-order self-energy is given in terms of chord diagrams with color-coding. Only one representative for each class is shown. Due to the absence of loops, an isomorphism between the Feynman diagrams and the chord diagrams can be established.
Scientific RepoRts | 7: 504 | DOI:10.1038/s41598-017-00355-w where ν k i n ( , ) is the integer number of absorbed plasmons in each fermionic line. Let us position 2n − 1 vertical lines such that they cut each fermionic line ( Fig. 3(a)). Then ν k i n ( , ) is computed as a number of bosonic lines crossing i th vertical line. Equation (5) can be derived by using the nonequilibrium Green's function (NEGF) formalism. Let a vertical line separate times lying on the forward and backward branches of the Keldysh contour in an expression for the lesser self-energy (Σ < ). Consider, for instance, a third vertical line at Fig. 3 is the lesser bosonic propagator. Performing three frequency integrals (over y 1 , y 2 , y 3 ) a contribution proportional to g < (ω + 3Ω) is obtained. Similar considerations can be repeated for each vertical line and fermionic propagator yielding in total 2n − 1 terms for each n th-order SE diagram . Now, since f i (ω) are non-singular the generic expression for the time-ordered self-energy (5) is obtained.
Equation (5) serves as the starting point for numerics; complexity goes into the generation of Feynman diagrams and the determination of the coefficients ) . Together with the Padé approximation this is the second important ingredient of our approach. The coefficients are computed purely algebraically by introducing an external time-dependent potential φ(1) (for brevity time variables are denoted as ≡ t i i ) and using the variational derivative technique 27 as in the derivation of Hedin's equations 28 . As was shown above the bosonic propagator in the present model does not renormalize (loops give zero contribution), i. e., = , leading to a simpler set of equations: where Γ(12, 3) is the vertex function, Σ(12) is the electron self-energy, and V(3) is the external plus the induced field in the system. All these quantities are functionally dependent on the external field φ(3) and on the full electron propagator g (12). The set of equations (6, 7 and 8) can now be iterated starting from leading to the diagrams shown at Fig. 3(b). The chord diagram 29,30 representation is natural in this case because according to the analysis above the fermionic loops yield zero contribution. In order to further facilitate the interpretation of the graphs in frequency space we use color coding for the coefficients ν k i n ( , ) entering the GF arguments. The graphs were generated by our symbolic algorithm in mathematica computer algebra system. Conversion from the time to the frequency domain is likewise performed using a symbolic algorithm. The self-energy which is accurate to the sixth order comprises 1, 1, 4, 27, 248, and 2830 diagrams of the first to sixth orders, respectively, has the following algebraic representation (see Supplementary Information for higher order terms): . Our explicit form for the self-energy dictates that the singularities of Σ(ω) should be located exactly at the GF poles. It is physically wrong as it is well known that the SE poles lie between the poles of the corresponding exact Green's function 31 . These two facts can be reconciled by noticing that already starting with the second order term γ ω ω ω the self-energy contains higher-order poles. As in our toy model (Fig. 1) they are responsible for the energy shift. However, the convergence is mathematically more intricate. In Fig. 4 the exact Green's function is compared with its reconstruction using SE of different orders (non-self-consistent calculation with Σ[g] being a functional of exact g). With increasing the perturbative order, the higher-order poles accumulate on the edges of lagoons (which enclose true simple poles) where the Green's function is not properly represented (vanishing real and imaginary parts). Outside the lagoons, the perturbative series converge rather rapidly as can be seen from the contour plots of tanh |g(ω)|. As in the toy model, PA can be used to obtain the Green's function in the whole complex plane. This will be illustrated now by sc calculations for the electron-boson model at equilibrium and at zero temperature.

Results
Self-consistent calculations at equilibrium and at T = 0. Assume that in the course of a self-consistent calculation an approximate GF (g (i) (ω)) has been obtained. Using the diagrammatic expansion (9) we evaluate the self-energy at a chosen frequency point. The point ω * should belong to the domain of convergence. In order to obtain the self-energy in the vicinity of the Green's function poles where the series diverges (note the unphysical multiple poles in the complex ω-plane on Fig. 5(b)), we perform the Padé approximation and use the new self-energy in the Dyson equation 1 and typically converge within some tens of cycles. Convergence is improved by using PA of variable order: on the first iteration cycle the non-interacting GF is used as an input leading to relatively simple self-energy that can regularized using PA of low order ([0/1]). In the course of sc calculations GF develops more satellites which require higher order PA (typically Scientific RepoRts | 7: 504 | DOI:10.1038/s41598-017-00355-w [11/12]) in order to accurately represent the self-energy. The quality of the resulting spectral function (cf. Fig. 5(b,c)) strongly depends on the order of perturbative expansions and on the electron-plasmon interaction strength γ. For the weakly correlated regime (γ . Ω  0 2 , we will show below that such a coupling strength is typical for monolayer graphene) already the GW approximation faithfully reproduces the exact spectral function, but it ceases to be valid in the correlated regime as demonstrated in Fig. 2. The energy of the QP peak is the major discrepancy. For γ = 0.65 (this value would be typical for the valence electrons in Al, Cu, Au metals), the third-order treatment substantially improves its position and strength, Fig. 5(a). Yet, the first satellite, which has a rather large contribution to the density of states at this value of γ (notice the logarithmic scale), represents a considerably more complicated feature. It can only be captured with a self-energy that is accurate to the 6th order (thin dark red line). However, even 3111 diagrams are not sufficient to reproduce the second-order satellite! There are known examples of GF calculations performed with even much larger number of diagrams, such as in diagrammatic quantum Monte Carlo 32 study of the Fröhlich polaron 33 , Anderson 34 and Hubbard 35 models. The major distinction of our approach is that it operates with skeletonic (in terms of dressed propagators from sc calculations) expansions. SE expansions in terms of bare propagators are reviewed by Cini and D' Andrea 12 . They can be represented in terms of continued fractions.  Estimates for realistic systems. Being able to perform higher order self-consistent diagrammatic calculations with high accuracy allows us to make predictions about the convergence of many-body perturbation theory (MBPT) for real-world systems, which is the main goal of this work. To this end, all relevant parameters for the mapping onto the electron-boson model need to be determined. Since  barely shifts the electronic spectrum, there is only one relevant dimensionless parameter a = γ/Ω; in its terms the correlation shift is given by Δ = a 2 Ω. The a-parameter for the hole states can be obtained by comparing the spectral strength of SE in the vicinity of  − Ω and  − 2Ω poles and be expressed in terms of the corresponding residues of the first and second order lesser self-energies (for the particle sector, similar expressions can be written in terms of the greater self-energies): (1) The second form in terms of the zeroth spectral moments of the lesser self-energies is preferable for realistic systems, where due to the momentum dispersions of electrons and plasmons the SE singularities are partially smeared out. As a paradigmatic system we consider here the homogeneous electron gas (HEG) 3 . It is believed to capture main electronic properties of simple metals, is a prototypic system for deriving approximations for the exchange-correlation functional 36 and has been widely studied using a variety of approaches [37][38][39][40][41][42] . Comparing with the electron-boson model, HEG is a considerably more complicated system: in addition to the excitations of multiple plasmon branches 43,44 , the electron-electron interaction is also accompanied by the excitation of particle-hole (p − h) pairs. These effects are more pronounced for the states close to the Fermi sphere (k = k F ). At the band bottom (the momentum k = 0) the phase-space for the excitation of p − h pairs is reduced and the mapping onto the electron-boson model is more justified. Therefore, this case will be considered here selecting (in conformance with the electron-boson model) only the SE diagrams describing the scattering processes with generation of plasmons. Such a selection is possible using the methods developed in refs 45 and 46. In particular it means that in the second order (in the screened Coulomb interaction) only has to be considered (pluses and minuses here denote the position of time-arguments on the Keldysh contour: −/+ for forward/backward branches). After analytically performing the frequency integrations over ω and over the internal frequencies we are left with the following momentum integrals: where α π = (4/9 ) 1/3 and r s is the Wigner-Seitz radius (in atomic units). The six-dimensional integration in M (2) can be somewhat simplified in spherical coordinates and reduced so to a 4d integral amenable to the Monte-Carlo approach introduced in ref. 47 and applied in refs 45, 46 and 48 to more complicated problems. Nonetheless it is instructive to introduce an approximation with the goal of estimating the a-parameter analytically.
We focus on the limit of large densities, α → r 0 s , where random phase approximation holds. In this case the critical momentum likewise approaches zero and we set ω ω ω ≈ t y y ( ) / ( ) pl pl pl 2 in the first integral leading to α ω ≈ M r q 2 s pl c (1) . In the second integral ω −  y y y ( ) Our numerical calculations (Fig. 6) confirm that despite a number of approximations the analytical estimate holds well in the weakly correlated limit. In the opposite limit the phase-space restriction due to n F (y 3 ) and the plasmonic dispersion gain significance reducing the electron-plasmon coupling as compared to the estimate (13).
Having the electron-plasmon interaction strength at our disposal, we can make some concrete predictions about the perturbative order required to accurately describe certain features in the electronic spectrum. As such, let us consider here the position of qp-peak and its first plasmon satellite. The horizontal lines in Fig. 6 denote the values of γ/Ω at which the sc theory yields these features with 10% accuracy. For instance, already the second-order self-energy is sufficient to describe both features in Al, whereas for Na 4th-order diagrams would be required to correctly compute the energetic position of the first plasmonic satellite! In contrast to these simple metals, the 2d systems are much more diverse in terms of their electron density parameter r s . By changing the substrate or by doping the system, it is possible to tune r s from the weakly correlated limit, such as in the case of monolayer graphene (for MLG r s ≈ 2.2/ M , M  is the background dielectric constant), to a strongly correlated regime in, e. g., 2DEG in GaAs where for the carrier concentration of n ≈ 10 9 cm 2 the density r s ≈ 13 has been reported 45 ( = r n 4/ s , = − n n/10 cm 10 2 ). At even higher densities a transition to the Wigner crystal phase takes place 46 . Our analysis (Fig. 6) thus justifies the use of GW approximation for the monolayer graphene 47,48 . Similar conclusion holds for fullerenes (r s ≈ 1.0), which are molecules of graphene wrapped by the introduction of pentagons on the hexagonal lattice, endorsing the use of GW for these systems 49,50 .
As can be seen from Fig. 6, with increasing r s the skeletonic diagrammatic SE expansion quickly becomes impractical. Therefore, different resummation methods such as sc parquet approximation (see Supplementary Information) need to be used.

Self-consistent calculations for excited states
So far results at zero temperature and zero bosonic occupation number have been presented. Equilibrium finite temperature scenario seems to be an obvious extension because the Wick theorem still holds 51 . Quite unexpectedly, calculations indicate that sc approach can only be realized at the lowest order, i. e., at the level of GW approximation. Higher-order self-energies cannot be reγularized with the help of PA. To understand this behavior, it is instructive to analyze the exact electron Green's function for a finite boson occupation number . A solution in terms of a continued fraction is known due to Cini 52 , but it also represent no difficulty to generalize the cumulant function (4) to this case. Using the standard finite temperature expression where as above = γ Ω a , and ( ) n k is the binomial coefficient, and the averaged boson occupation number is n b . The poles are now situated on both sides of the qp-peak. In addition to the poles associated with the excitation of bosons, there are poles associated with the energy absorption from the thermally excited bosons. Thus, real ω * points cannot be used for the Padé approximation because they are inevitably situated in the proximity of poles. In fact, by plotting an approximate GF in the complex plane (Fig. 7) we see that the situation cannot be cured even by shifting ω * away from the real axis: the lagoon that encompasses regions of non-convergence is extended along real and imaginary axes of the complex plane posing problems for the regularization.
In thermal equilibrium at finite β (the inverse temperature) the mixed bosonic state is represented by the density matrix ρ β . It is of interest to find the electron Green's function for pure bosonic states corresponding to a given n b . These states are particularly relevant for the state-of-the-art ultrafast experiments 53,54 where the interaction with laser pulses cannot be assumed to follow the adiabatic path 22, 55, 56 . Excitation of confined plasmons in graphene by impinging free electrons 57 can be considered as a paradigmatic system. The analysis of García de Abajo suggests that 100 eV electrons interacting with a monolayer graphene excite on average one plasmon per electron. This is a sizeable effect that can be detected with standard spectroscopic methods as discussed below. To describe photoemission from the system excited by impinging electrons, we compute the single-particle GF for n b > 0 using the Feynman disentangling of operators (see Supplemental Information). Expansions of g [nb] (ω) in terms of the shifted n b = 0 propagators ω ≡ g g( ) [0] provide an exact solution for the particular case of electron-boson model (3). However, perturbative sc calculations using some approximations for the electron SE have additional advantage that they can be generalized to more complicated scenarios, e. g. include dispersion and multiple electronic and bosonic bands. Let us recall that perturbative expansions of correlators in MBPT (including the self-energies) are generated by expanding the contour evolu- in powers of the interaction and expressing averages of the operator products in terms of products of simple propagators (the Wick theorem). For the ground state (generalizations to arbitrary initial states are also possible 51 ) of a system of fermions the time-ordered product of any number of field operators splits up into the sum of the products of normal products of pairs, the averages of normal products being equal to zero. For bosonic systems with n b > 0 a straightforward generalization of Eq. (9) would not work since it relies upon the Wick theorem, which cannot be formulated here because averages of the normal product of bosonic operators are non-zero 58 . Nonetheless a method to generate excited states self-energies can be devised.
The method can be illustrated by considering the computation of bosonic averages, , over the states with fixed particle number, where square brackets are used to distinguish the state with fixed boson number from the thermal state. As a particular example, n b = 1 will be computed in accordance with the scenario of plasmons in graphene excited by means of impinging electrons. Recalling that bosonic displacement operator is given by one can write: n n 1 2 where  is the contour ordering operator, and for the computation of the correlator on the right hand side we first (where  is the order relation with respect to  ) and subsequently evaluate the limit τ → 0. Now, on the right hand side we have a n b = 0 correlator which can be computed using the standard Wick theorem. The procedure is also suitable for the computation of the electron self-energy because the conditions of validity of the Wick theorem for fermionic degrees of freedom are not affected by the choice of a reference bosonic state. By computing corresponding bosonic correlators we first arrive at the SE being a functional of bare GFs [1] ( 0) , see Supplemental Information for the explicit form), and by iterating the Dyson equation further express the self-energy in terms of the full propagators Σ = Σ g [ ] [1] [ 1] . This leads to a generalization of Eq. (9) for the n b = 1 pure bosonic state: g g g g g g g g g g g g g g g g g g gg g g g g  In contrast to the n b = 0 case, SE in terms of the dressed propagators given by the equation above possesses less economical series than Σ g [ ] r (0) . This can be understood by considering the lowest order Σ [1] which is identical for the nonequilibrium and thermal states (latter can be inferred from the first term of Equation (9) Figure 7. Electron Green's function g β (ω) resolved in complex plane for Ω = 1, γ = 0.65 and n b = 1/2 computed using the forth order SE. Color coding as in Fig. 4. Exact spectral function is superimposed.
the single-shot calculation yields a single peak above the qp state, however, multiple satellites above the qp state are created if the procedure is iterated to self-consistency. All of them but one quickly diminish when higher-order SE terms are included. The cancellation is achieved owing to extra terms containing g −1 (i > 1); they do not appear (0) . Corresponding sc results are obtained using the Padé approximation, Fig. 8(a). The extra satellite at ε + + Ω γ Ω 2 , is a marked feature of the spectrum that can be observed in time-resolved photoemission. A possible scenario for such an experiment is depicted in Fig. 8(b): as in the proposal of García de Abajo 57 impinging free electrons excite the bosonic subsystem -plasmons. The same effect can, in principle, be achieved with ordinary laser pulses, however, for confined systems, direct optical excitation of plasmons is less efficient. On the second step, optical or UV pulse detects the change in the electronic density of states provided the delay between the electron pump and the photon probe does not exceed the plasmon relaxation time. Thus, the proposed experiment is capable of directly measuring the electron-plasmon interaction strength and, if resolved in time, yielding the plasmon relaxation time.

Discussion
There is more than a computational complexity which prevents the applications of the MBPT beyond the leading order. The resulting asymptotic series lead to Green's functions with incorrect physical properties such as non-positive densities and higher-order poles already at the second order [59][60][61] . Besides the interaction strength, the domain of convergence strongly depends on the microscopic details of the model: continuous space vs. lattice formulation 62 and also on the temperature as discussed below. For various statistical and many-body models PA has been used to extend perturbative expansions beyond their domain of convergence [63][64][65][66] . The same mathematical approach is used here in a different context, to regularize the electron SE. Using NEGF formalism the self-energy of the electron-boson model is derived in an explicit form for the ground n b = 0 (Equation 9) and excited n b = 1 (Equation 16) bosonic states and a connection of its diagrammatic expansion to a certain class of chord diagrams is demonstrated. With the help of the developed symbolic algorithm higher-order self-consistent calculations were performed, Figs 5 and 8. More complicated models can be treated starting from the same idea and applying the momentum average approximation 67 , which was proved to be accurate for the description of dressed particles in the Holstein polaron model 68 .
To achieve the main goal of this work and make quantitative predictions on the convergence of MBPT for realistic systems, we map the homogeneous electron gas onto the studied model. It allows to express our finding in terms of the density parameter r s (Fig. 6). We find that the number of self-energy terms needed to be included in the sc calculations grows rapidly with r s : while for r s = 1 (typical for Carbon materials) the first order accurately describes the first plasmon satellite, even the third order SE is not sufficient at metallic densities (r s = 4).
Finally, extension of the formalism towards electron-boson systems in an excited bosonic state represents an interesting conceptual problem and is relevant for the description of ultrafast spectroscopic experiments. For thermal states at finite temperature, there is a clear mathematical argument precluding self-consistent calculations. However, we demonstrate, the calculations can be performed for pure states with finite boson occupation number (n b > 1). A method is derived, allowing to obtain perturbative expansions even though the Wick theorem is not applicable in this case. For n b = 1, our calculations show additional plasmonic satellite above the quasiparticle peak. Even at rather low electronic density, such as found in graphene, the modification of the spectral function compared to the n b = 0 state is substantial and can be experimentally observed. However, for confined systems 69 , direct optical excitation of plasmons is inefficient 70,71 . Therefore, we argue, by exciting the system using a beam of free electrons, the appearance of extra satellites in the electronic density of states can be probed by the time-resolved photoemission 22,72,73 . This setup allows to directly measure the electron-plasmon interaction strength and the plasmon relaxation time.