Physical and unphysical regimes of self-consistent many-body perturbation theory

In the standard framework of self-consistent many-body perturbation theory, the skeleton series for the self-energy is truncated at a finite order $\mathcal{N}$ and plugged into the Dyson equation, which is then solved for the propagator $G_{\mathcal{N}}\,$. We consider two examples of fermionic models, the Hubbard atom at half filling and its zero space-time dimensional simplified version. First, we show that $G_{\mathcal{N}}\,$ converges when $\mathcal{N}\to\infty$ to a limit $G_\infty\,$, which coincides with the exact physical propagator $G_{\rm exact}\ $ at small enough coupling, while $G_\infty \neq G_{\rm exact}\ $ at strong coupling. This follows from the findings of [Kozik, Ferrero and Georges, PRL 114, 156402 (2015)] and an additional subtle mathematical mechanism elucidated here. Second, we demonstrate that it is possible to discriminate between the $G_\infty=G_{\rm exact}\ $ and $G_\infty\neq G_{\rm exact}\ $ regimes thanks to a criterion which does not require the knowledge of $G_{\rm exact}\ $, as proposed in [Rossi et al., PRB 93, 161102(R) (2016)].


Introduction
Self-consistent perturbation theory is a particularly elegant and powerful approach in quantum many-body physics [2,3]. The single-particle propagator G is expressed through the Dyson equation in terms of the non-interacting propagator G 0 and the self-energy Σ, which itself is formally expressed in terms of G through the skeleton series where Σ (n) bold [G] is the sum of all two-particle irreducible Feynman diagrams of order n (built with bold propagator lines representing G).
The standard procedure for solving Eqs. (1,2) is to truncate the skeleton series at a finite order N , and to look for the solution G N of the self-consistency equation 1 with Σ bold . The natural expectation is that one obtains the exact propagator by sending the truncation order to infinity: G N → G exact for N → ∞.
However, as was discovered in [4], the series Σ (≤N ) bold [G exact ] can converge when N → ∞ to a result which differs from the exact physical self-energy Σ exact = G −1 0 − G −1 exact . This misleading convergence phenomenon was observed for three fermionic textbook models -Hubbard atom, Anderson impurity model, and half-filled 2D Hubbard model-in a region of the parameter space (at and around half filling, at strong interaction and low temperature). G exact was computed with a numerically exact quantum Monte Carlo method, and the skeleton series was evaluated up to N = 6 or 8 by diagrammatic Monte Carlo [5]. Numerous works [1,[6][7][8][9][10][11][12][13][14] have studied various aspects of the problem found in [4], as well as the related divergences of irreducible vertices ( [8,10,11,13,[15][16][17][18] and Refs. therein). In particular, Ref. [7] introduced an exactly solvable toy model, which has the structure of a fermionic model in zero space-time dimensions, and features the misleading convergence problem of [4].
In this article, we study the consequences of this problem for the sequence G N , which is the crucial question in the most relevant cases where G exact is unknown. For the toy model of [7], we find that G N converges when N → ∞ to a limit G ∞ which differs from G exact at strong coupling. For the Hubbard atom, our numerical data strongly indicate that such misleading convergence of the sequence G N also occurs at large coupling and half filling. Moreover, we demonstrate that a criterion proposed in [1] allows to discriminate between the G ∞ = G exact and G ∞ = G exact regimes without using the knowledge of G exact .
We note that although we restrict here to the scheme (1,2) where G is the only bold element, our findings may also be relevant to other schemes containing additional bold elements, such as a bold interaction line W , or a bold pair propagator line Γ. The scheme built with G and W is natural for Coulomb interactions, and is widely used for solids and molecules with a truncation order N = 1 (the GW approximation) and sometimes with N = 2 (see, e.g., Refs. [19][20][21][22]), while for several paradigmatic lattice models, bold diagrammatic Monte Carlo (BDMC) made it possible to reach larger N and claim a small residual trunction error [23][24][25][26]. The scheme built with G and Γ is natural for contact interactions; truncation at order N = 1 then corresponds to the self-consistent T-matrix approximation [27][28][29], and precise large-N results were obtained by BDMC in the normal phase of the Hubbard model [30,31] and of the unitary Fermi gas [32][33][34]. Other BDMC results were obtained for models of coupled electrons and phonons, where it is natural to introduce a bold phonon propagator [23,35], and for frustrated spins [36][37][38]. Schemes containing three-or four-point bold vertices were also employed, to construct extensions of dynamical mean-field theory [17,39].
2 Zero space-time dimensional toy-model

Definitions and reminders
We begin with some reminders from [7]. While fermionic many-body problems can be represented by a functional integral over Grassmann fields, which depend on d space coordinates and one imaginary time coordinate [40], in this simplified toy model the Grassmann fields are replaced with Grassmann numbers ϕ s andφ s that do not depend on anything, apart from a spin index s ∈ {↑, ↓}. The partition function, the action and the propagator are then defined by the dimensionless parameters µ and U being the analogs of chemical potential and interaction strength. We restrict for convenience to µ > 0 (changing the sign of µ essentially amounts to the change of variables ϕ ↔φ) and to U < 0 (as in [7]). The coefficients of the skeleton series have the analytical expression It is convenient to work with rescaled variables, multiplying propagators with |U | and dividing self-energies with the same factor, The rescaled skeleton series is then given by bold (g) = a n (−1) n g 2n−1 and accordingly σ bold (g). The exact self-energy and propagator are given by in terms of the rescaled free propagator g 0 := |U | G 0 = |U |/µ.
If one evaluates the bold series at the exact G, one obtains the correct physical selfenergy for |U | < µ 2 and an incorrect result for |U | > µ 2 . More precisely, the self-energy functional (which reduces to a function in this toy model) has the two branches as represented on Fig. 1. The physical branch is the (+) branch for g 0 < 1, and the (−) branch for g 0 > 1; i.e., σ exact (g 0 ) = σ (sign(1−g 0 )) (g exact (g 0 )). On the other hand, the bold series, evaluated at the exact physical propagator, always converges to the (+) branch; i.e., σ bold (g exact (g 0 )) = σ (+) (g exact (g 0 )) for all g 0 > 0. Note that σ bold (g) is the expansion of σ (+) (g) in powers of g, and thus from (5) the convergence radius of the series σ bold (g) is 1/2. Figure 1: The two branches of the self-energy as a function of the full propagator, for the toy model in zero space-time dimensions. The skeleton series converges up to g = 1/2 and coincides with the (+) branch: σ bold (g) = σ (+) (g) for g ≤ 1/2.

Limit of the skeleton sequence
We will refer to the sequence G N as the skeleton sequence. Rescaling variables as in (4), in particular setting g N := G N |U |, the self-consistency equation (3) becomes This equation is readily solved for g N numerically: The solutions are roots of a polynomial of order 2N , and we observe that there is a unique real positive root, which we take to be g N (recall that the exact g is always real and positive); alternatively, we solved Eq. (6) by iterations (with a damping procedure described in the next Section), and we found convergence to this same g N . We find that i.e., the skeleton sequence converges to the correct physical result below a critical coupling strength, and to an unphysical result above it. Let us focus on the regime g 0 > 1, where the convergence to an unphysical result takes place (as demonstrated in Fig. 2). The fact that the skeleton sequence converges at all in this regime is non-trivial. The value of the unphysical limit g ∞ = 1/2 of the skeleton sequence g N is equal to the radius of convergence of the skeleton series σ bold (g). This is not a coincidence, and the reason for this self-tuning towards the convergence radius becomes clear from Fig. 3: For a large truncation order, the curve representing the truncated bold series as a function of g becomes an almost vertical line above the position of the convergence radius (g = 1/2), so that it intersects the Dyson-equation curve near this value of g. It also becomes clear that we are in an unusual situation where  bold (g) at g = g N , whereas the exact propagator g = g exact is given by the intersection of the Dyson-equation curve with the physical branch σ (sign(1−g 0 )) (g). It appears clearly that for g 0 < 1, g N converges to the exact g, while for g 0 > 1, g N always tends to 1/2, the convergence radius of the skeleton series.

Diagnosing the misleading convergence
In the general case where G exact is unknown, when one observes numerically that G N converges to some limit, one needs a way to tell whether or not this limit is equal to G exact . Assuming that G N → G ∞ for N → ∞, the following criterion [1] is a sufficient condition for G ∞ to be equal to G exact : There exists > 0 such that for any ξ in the disc D = {|ξ| < 1 + }, For all practical purpose, we expect this criterion to be essentially equivalent to the following simpler one: There exists ξ > 1 such that Σ N,ξ converges for N → ∞.
In what follows we will use this simplified criterion. We also introduce an extra factor 1/ξ N 0 in the definition (7) of Σ N,ξ , where the value of N 0 will be conveniently chosen; such an N -independent factor does not matter for the criterion (it does not change whether or not the sequence Σ N,ξ converges).
For the toy-model, this means that assuming g N → g ∞ for N → ∞, a sufficient condition for g ∞ to be equal to the correct physical g exact (g 0 ) is that there exists ξ > 1 such that bold (g N ξ) / ξ (9) converges for N → ∞. As illustrated in Fig. 4, this criterion indeed allows to detect the misleading convergence for g 0 > 1, and to trust the result for g 0 < 1.

Hubbard atom
We turn to the single-site Hubbard model, defined by the grand-canonical Hamiltonian −µ s n s + U n ↑ n ↓ . The propagator can be expressed as a functional integral over βantiperiodic Grassmann fields [40], with the action and We restrict for simplicity to the half-filled case µ = U/2, which should be the most dangerous case, since it is at and around half-filling that the misleading convergence of Σ bold [G exact ] was discovered in [4]. We use the BDMC method [5,32,41,42] to sum all skeleton diagrams and solve the self-consistency equation (3) for truncation orders N ≤ 8 (note that at half filling, Σ (n) bold = 0 for all odd n > 1). The first question is whether the skeleton sequence G N can also converge to an unphysical result, or equivalently, whether Σ The next question is whether the criterion (8) allows us to discriminate between these two situations. We therefore plot the sequence Σ N,ξ in Figs. 6 and 7. We only show the imaginary part because in the considered half-filled case, the real part of Σ N (ω n ) automatically equals U/2; moreover we focus for simplicity on the lowest Matsubara frequency ω 0 = π/β, and we choose N 0 = 2.
For ξ = 1, Σ N,ξ reduces to the original skeleton sequence Σ N , and the behavior is similar to the double occupancy: The sequence appears to converge, albeit slowly, towards an unphysical result for βU = 8 (Fig. 6), while fast convergence to the correct physical result takes place for βU = 1 (Fig. 7). For ξ > 1, the sequence does not appear to converge any more for βU = 8, see Fig. 6: The criterion correctly indicates that the results should not be trusted in this case. In contrast, for βU = 1, the criterion allows to validate the results, since the sequence remains convergent at ξ > 1, see Fig. 7. Regarding the choice of ξ, it should be neither too small in order to have an effect at the accessible orders, nor too large to avoid making the criterion too conservative; here we see that ξ = 1.1 and 1.2 are appropriate.    . Such a damping procedure is commonly used in BDMC where it also reduces the statistical error [42,43]. In the toy model, one can easily show that an increasingly strong damping is required when N is increased, because for N → ∞, the slope [dσ (≤N ) bold (g)/dg] g=g N diverges, making the undamped iterative procedure unstable. This observation could also be useful for misleading-convergence detection.
Finally, we comment on the link with the multivaluedness of the self-energy functional Σ[G] (i.e., of the Luttinger-Ward functional). In [4], the misleading convergence of the skeleton series was found to be towards an unphysical branch of the self-energy functional, in the sense that if Eqs. (10,11) are viewed as a mapping G 0 → G[G 0 ], then there exists ]. As noted in [4], this G 0,unphys does not belong to the set of physical bare propagators which are of the form (12) for some value of chemical potential; therefore, by looking at G 0,unphys , one can tell that the result is on an unphysical branch, and hence detect the misleading convergence of the skeleton series. In contrast, the misleading convergence of the skeleton sequence found here cannot be detected in this way. Indeed, the self-consistency equation (3) is enforced with the original physical G 0 .

Conclusion
We have observed that there is a regime where the solution of self-consistent many-body perturbation theory converges to an unphysical result in the limit of infinite truncation order of the skeleton series. This surprising breakdown of the standard framework results from a subtle mathematical mechanism which we have clarified by analyzing the zero space-space dimensional model. In this problematic regime, lowest order calculations can be off by one order of magnitude, but access to higher orders allows to detect the problem numerically through the divergence of a slightly modified sequence, whereas seeing convergence of this modified sequence allows to rule out misleading convergence and to trust the result, as proposed in [1] and demonstrated here for the Hubbard atom. Such a proof of principle is relevant for many-body problems in regimes where, in spite of important progress with non self-consistent frameworks [44][45][46][47][48][49][50][51][52][53][54] (for which misleading convergence generically does not occur [1]), self-consistent BDMC remains the state of the art to date. In particular, the present findings served as a basis to discriminate between physical and unphysical BDMC results for the doped two-dimensional Hubbard model at strong coupling in a non-Fermi liquid regime [55]. and Y.D. from the National Natural Science Foundation of China (grant No. 11625522) and the Science and Technology Committee of Shanghai (grant No. 20DZ2210100). The Flatiron Institute is a division of the Simons Foundation.