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

Kris van Houcke, Evgeny Kozik, Riccardo Rossi, Youjin Deng, Félix Werner

► To cite this version:

Kris van Houcke, Evgeny Kozik, Riccardo Rossi, Youjin Deng, Félix Werner. Physical and un-physical regimes of self-consistent many-body perturbation theory. SciPost Physics, 2024, pp.133. 10.21468/SciPostPhys.16.5.133. hal-03365949

HAL Id: hal-03365949

https://hal.science/hal-03365949

Submitted on 31 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Physical and unphysical regimes of self-consistent many-body perturbation theory

Kris van Houcke, Evgeny Kozik, Riccardo Rossi, Youjin Deng, Félix Werner

► To cite this version:

Kris van Houcke, Evgeny Kozik, Riccardo Rossi, Youjin Deng, Félix Werner. Physical and unphysical regimes of self-consistent many-body perturbation theory. SciPost Physics, 2024, 16 (5), pp.133. 10.21468/SciPostPhys.16.5.133. hal-03365949

HAL Id: hal-03365949

https://hal.science/hal-03365949v1

Submitted on 31 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Physical and unphysical regimes of self-consistent many-body perturbation theory

K. Van Houcke¹, E. Kozik², R. Rossi¹,³,⁴,⁵, Y. Deng⁶,⁷, and F. Werner⁸*

1 Laboratoire de Physique de l'École normale supérieure, ENS - Université PSL, CNRS, Sorbonne Université, Université Paris Cité, 75005 Paris, France

2 Physics Department, King's College, London WC2R 2LS, United Kingdom

3 Center for Computational Quantum Physics, Flatiron Institute, New York, NY 10010, USA

4 CNRS, LPTMC, Sorbonne Université, 75005 Paris, France

5 Institute of Physics, EPFL, 1015 Lausanne, Switzerland

6 National Laboratory for Physical Sciences at Microscale and Department of Modern Physics, University of Science and Technology of China, Hefei, Anhui 230026, China

7 Shanghai Research Center for Quantum Science, Shanghai 201315, China

8 Laboratoire Kastler Brossel, ENS - Université PSL, CNRS, Collège de France, Sorbonne Université, 75005 Paris, France

* werner@lkb.ens.fr

May 27, 2024

## Abstract

In the standard framework of self-consistent many-body perturbation theory, the skeleton series for the self-energy is truncated at a finite order $N$ and plugged into the Dyson equation, which is then solved for the propagator $G_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_N$ converges when $N \to \infty$ to a limit $G_\infty$, which coincides with the exact physical propagator $G_{exact}$ at small enough coupling, while $G_\infty \neq G_{exact}$ at strong coupling. This follows from the findings of [1] and an additional subtle mathematical mechanism elucidated here. Second, we demonstrate that it is possible to discriminate between the $G_\infty = G_{exact}$ and $G_\infty \neq G_{exact}$ regimes thanks to a criterion which does not require the knowledge of $G_{exact}$, as proposed in [2].

## Contents

1 Introduction 2

2 Zero space-time dimensional toy-model 3

2.1 Definitions and reminders 3

2.2 Limit of the skeleton sequence 5

2.3 Diagnosing the misleading convergence 6

3 Hubbard atom 7

4 Conclusion 10

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

K. Van Houcke¹, E. Kozik², R. Rossi¹,³,⁴,⁵, Y. Deng⁶,⁷, and F. Werner⁸*

1 Laboratoire de Physique de l'École normale supérieure, ENS - Université PSL, CNRS, Sorbonne Université, Université Paris Cité, 75005 Paris, France

2 Physics Department, King's College, London WC2R 2LS, United Kingdom

3 Center for Computational Quantum Physics, Flatiron Institute, New York, NY 10010, USA

4 CNRS, LPTMC, Sorbonne Université, 75005 Paris, France

5 Institute of Physics, EPFL, 1015 Lausanne, Switzerland

6 National Laboratory for Physical Sciences at Microscale and Department of Modern Physics, University of Science and Technology of China, Hefei, Anhui 230026, China

7 Shanghai Research Center for Quantum Science, Shanghai 201315, China

8 Laboratoire Kastler Brossel, ENS - Université PSL, CNRS, Collège de France, Sorbonne Université, 75005 Paris, France

* werner@lkb.ens.fr

May 27, 2024

## Abstract

In the standard framework of self-consistent many-body perturbation theory, the skeleton series for the self-energy is truncated at a finite order $N$ and plugged into the Dyson equation, which is then solved for the propagator $G_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_N$ converges when $N \to \infty$ to a limit $G_\infty$, which coincides with the exact physical propagator $G_{exact}$ at small enough coupling, while $G_\infty \neq G_{exact}$ at strong coupling. This follows from the findings of [1] and an additional subtle mathematical mechanism elucidated here. Second, we demonstrate that it is possible to discriminate between the $G_\infty = G_{exact}$ and $G_\infty \neq G_{exact}$ regimes thanks to a criterion which does not require the knowledge of $G_{exact}$, as proposed in [2].

## Contents

1 Introduction 2

2 Zero space-time dimensional toy-model 3

2.1 Definitions and reminders 3

2.2 Limit of the skeleton sequence 5

2.3 Diagnosing the misleading convergence 6

3 Hubbard atom 7

4 Conclusion 10

Appendix 11
# 1 Introduction

Self-consistent perturbation theory is a particularly elegant and powerful approach in quantum many-body physics [3–5]. The single-particle propagator $G$ is expressed through the Dyson equation

$$G^{-1} = G_0^{-1} - \Sigma \qquad (1)$$

in terms of the non-interacting propagator $G_0$ and the self-energy $\Sigma$, which itself is formally expressed in terms of $G$ through the skeleton series,

$$\Sigma = \Sigma_{\text{bold}}[G] \equiv \sum_{n=1}^{\infty} \Sigma_{\text{bold}}^{(n)}[G] \qquad (2)$$

where $\Sigma_{\text{bold}}^{(n)}[G]$ is the sum of all skeleton self-energy Feynman diagrams of order $n$ (these diagrams are built with bold propagator lines representing $G$, and remain connected if one cuts one or two $G$-lines).

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<fnref n="1" />

$$G_N^{-1} = G_0^{-1} - \Sigma_{\text{bold}}^{(\le N)}[G_N] \qquad (3)$$

with

$$\Sigma_{\text{bold}}^{(\le N)} := \sum_{n=1}^{N} \Sigma_{\text{bold}}^{(n)}.$$

The natural expectation is that one obtains the exact propagator by sending the truncation order to infinity: $G_N \to G_{\text{exact}}$ for $N \to \infty$.

However, as was discovered in [1], the series $\Sigma_{\text{bold}}^{(\le N)}[G_{\text{exact}}]$ can converge when $N \to \infty$ to a result which differs from the exact physical self-energy $\Sigma_{\text{exact}} = G_0^{-1} - G_{\text{exact}}^{-1}$. 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_{\text{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 [6]. Numerous works [2, 7–15] have studied various aspects of the problem found in [1], as well as the related divergences of irreducible vertices ([9, 11, 12, 14, 16–20] and Refs. therein). In particular, Ref. [8] 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 [1], as well as the related multivaluedness of the Luttinger-Ward functional also discovered in [1].

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_{\text{exact}}$ is unknown. For the toy model of [8], we find that $G_N$ converges when $N \to \infty$ to a limit $G_\infty$ which differs from $G_{\text{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. This misleading convergence of $G_N$ is the first result reported in this article. Secondly, we present data, again for the toy model of [8] and for the Hubbard atom, demonstrating that a criterion proposed in [2] enables one to discriminate between the $G_\infty \neq G_{\text{exact}}$ and $G_\infty = G_{\text{exact}}$ regimes without using the knowledge of $G_{\text{exact}}$.

<footnote n="1">¹We assume that the solution $G_N$ of (3) is unique, or at least that there is no difficulty in identifying a unique potentially physical solution (e.g., by starting from the weakly interacting limit where $G_N \to G_0$, and following the solution as a function of interaction strength).</footnote>
# 1 Introduction

Self-consistent perturbation theory is a particularly elegant and powerful approach in quantum many-body physics [3–5]. The single-particle propagator $G$ is expressed through the Dyson equation

$$G^{-1} = G_0^{-1} - \Sigma \qquad (1)$$

in terms of the non-interacting propagator $G_0$ and the self-energy $\Sigma$, which itself is formally expressed in terms of $G$ through the skeleton series,

$$\Sigma = \Sigma_{\text{bold}}[G] \equiv \sum_{n=1}^{\infty} \Sigma_{\text{bold}}^{(n)}[G] \qquad (2)$$

where $\Sigma_{\text{bold}}^{(n)}[G]$ is the sum of all skeleton self-energy Feynman diagrams of order $n$ (these diagrams are built with bold propagator lines representing $G$, and remain connected if one cuts one or two $G$-lines).

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<fnref n="1" />

$$G_N^{-1} = G_0^{-1} - \Sigma_{\text{bold}}^{(\le N)}[G_N] \qquad (3)$$

with

$$\Sigma_{\text{bold}}^{(\le N)} := \sum_{n=1}^{N} \Sigma_{\text{bold}}^{(n)}.$$

The natural expectation is that one obtains the exact propagator by sending the truncation order to infinity: $G_N \to G_{\text{exact}}$ for $N \to \infty$.

However, as was discovered in [1], the series $\Sigma_{\text{bold}}^{(\le N)}[G_{\text{exact}}]$ can converge when $N \to \infty$ to a result which differs from the exact physical self-energy $\Sigma_{\text{exact}} = G_0^{-1} - G_{\text{exact}}^{-1}$. 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_{\text{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 [6]. Numerous works [2, 7–15] have studied various aspects of the problem found in [1], as well as the related divergences of irreducible vertices ([9, 11, 12, 14, 16–20] and Refs. therein). In particular, Ref. [8] 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 [1], as well as the related multivaluedness of the Luttinger-Ward functional also discovered in [1].

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_{\text{exact}}$ is unknown. For the toy model of [8], we find that $G_N$ converges when $N \to \infty$ to a limit $G_\infty$ which differs from $G_{\text{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. This misleading convergence of $G_N$ is the first result reported in this article. Secondly, we present data, again for the toy model of [8] and for the Hubbard atom, demonstrating that a criterion proposed in [2] enables one to discriminate between the $G_\infty \neq G_{\text{exact}}$ and $G_\infty = G_{\text{exact}}$ regimes without using the knowledge of $G_{\text{exact}}$.

<footnote n="1">¹We assume that the solution $G_N$ of (3) is unique, or at least that there is no difficulty in identifying a unique potentially physical solution (e.g., by starting from the weakly interacting limit where $G_N \to G_0$, and following the solution as a function of interaction strength).</footnote>
The misleading convergence of $G_N$ reported here is a non-trivial fact. It comes from a subtle mathematical mechanism (as we will see), and does not merely follow from the misleading convergence of $\Sigma_{\text{bold}}^{(\le N)}[G_{\text{exact}}]$ discovered in [1]. Indeed, a naive reasoning would suggest that if the misleading convergence of $\Sigma_{\text{bold}}^{(\le N)}[G_{\text{exact}}]$ takes place, then $G_N$ should not converge at all.<fnref n="2" />

We restrict here to the scheme (1,2) where $G$ is the only bold element (as in, e.g., Ref. [21]). Nevertheless, 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 $\Gamma$. 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. [22–25]), while for several paradigmatic lattice models, bold diagrammatic Monte Carlo (BDMC) made it possible to reach larger $N$ and claim a small residual truncation error [26–29]. The scheme built with $G$ and $\Gamma$ is natural for contact interactions; truncation at order $N=1$ then corresponds to the self-consistent T-matrix approximation [30–32], and precise large-$N$ results were obtained by BDMC in the normal phase of the Hubbard model [33, 34] and of the unitary Fermi gas [35–37]. Other BDMC results were obtained for models of coupled electrons and phonons, where it is natural to introduce a bold phonon propagator [26, 38], and for frustrated spins [39–41]. Schemes containing three- or four-point bold vertices were also employed, to construct extensions of dynamical mean-field theory [18, 42].

# 2 Zero space-time dimensional toy-model

## 2.1 Definitions and reminders

We begin with some reminders from [8] (see [8] for the derivations). 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 [4, 43], in this simplified toy model the Grassmann fields are replaced with Grassmann numbers $\varphi_s$ and $\bar{\varphi}_s$ that do not depend on anything, apart from a spin index $s \in \{\uparrow, \downarrow\}$. The partition function, the action and the propagator are then defined by

$$ 
\begin{align*}
Z &= \int \left( \prod_s d\varphi_s \, d\bar{\varphi}_s \right) e^{-S} \\
S &= -\mu \sum_s \bar{\varphi}_s \varphi_s + U \bar{\varphi}_\uparrow \varphi_\uparrow \bar{\varphi}_\downarrow \varphi_\downarrow \\
G &= -\frac{1}{Z} \int \left( \prod_s d\varphi_s \, d\bar{\varphi}_s \right) \varphi_{s'} \bar{\varphi}_{s'} e^{-S},
\end{align*}
$$

the dimensionless parameters $\mu$ and $U$ being the analogs of chemical potential and interaction strength. Since $G$ is spin-independent, we omit its spin index. We restrict for convenience to $\mu > 0$ (changing the sign of $\mu$ essentially amounts to the change of variables $\varphi \leftrightarrow \bar{\varphi}$) and to $U < 0$ (as in [8]).

The coefficients of the skeleton series have the analytical expression

$$ \Sigma_{\text{bold}}[G] = \sum_{n=1}^{\infty} a_n G^{2n-1} U^n \quad \text{with} \quad a_n = \frac{(-1)^{n-1} (2n-2)!}{n!(n-1)!}. $$

<footnote n="2">²The naive reasoning goes as follows: If $G_N$ would converge to some $G_\infty$ for $N \to \infty$, then, from (3), one can expect $G_\infty^{-1} = G_0^{-1} - \Sigma_{\text{bold}}[G_\infty]$, and hence, assuming unicity of the solution of the Dyson equation, $G_\infty = G_{\text{exact}}$. Thus $\Sigma_{\text{bold}}[G_{\text{exact}}] = \Sigma_{\text{exact}}$, in contradiction with the misleading convergence of $\Sigma_{\text{bold}}^{(\le N)}[G_{\text{exact}}]$.</footnote>
The misleading convergence of $G_N$ reported here is a non-trivial fact. It comes from a subtle mathematical mechanism (as we will see), and does not merely follow from the misleading convergence of $\Sigma_{\text{bold}}^{(\le N)}[G_{\text{exact}}]$ discovered in [1]. Indeed, a naive reasoning would suggest that if the misleading convergence of $\Sigma_{\text{bold}}^{(\le N)}[G_{\text{exact}}]$ takes place, then $G_N$ should not converge at all.<fnref n="2" />

We restrict here to the scheme (1,2) where $G$ is the only bold element (as in, e.g., Ref. [21]). Nevertheless, 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 $\Gamma$. 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. [22–25]), while for several paradigmatic lattice models, bold diagrammatic Monte Carlo (BDMC) made it possible to reach larger $N$ and claim a small residual truncation error [26–29]. The scheme built with $G$ and $\Gamma$ is natural for contact interactions; truncation at order $N=1$ then corresponds to the self-consistent T-matrix approximation [30–32], and precise large-$N$ results were obtained by BDMC in the normal phase of the Hubbard model [33, 34] and of the unitary Fermi gas [35–37]. Other BDMC results were obtained for models of coupled electrons and phonons, where it is natural to introduce a bold phonon propagator [26, 38], and for frustrated spins [39–41]. Schemes containing three- or four-point bold vertices were also employed, to construct extensions of dynamical mean-field theory [18, 42].

# 2 Zero space-time dimensional toy-model

## 2.1 Definitions and reminders

We begin with some reminders from [8] (see [8] for the derivations). 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 [4, 43], in this simplified toy model the Grassmann fields are replaced with Grassmann numbers $\varphi_s$ and $\bar{\varphi}_s$ that do not depend on anything, apart from a spin index $s \in \{\uparrow, \downarrow\}$. The partition function, the action and the propagator are then defined by

$$ 
\begin{aligned}
Z &= \int \left( \prod_s d\varphi_s \, d\bar{\varphi}_s \right) e^{-S} \\
S &= -\mu \sum_s \bar{\varphi}_s \varphi_s + U \bar{\varphi}_\uparrow \varphi_\uparrow \bar{\varphi}_\downarrow \varphi_\downarrow \\
G &= -\frac{1}{Z} \int \left( \prod_s d\varphi_s \, d\bar{\varphi}_s \right) \varphi_{s'} \bar{\varphi}_{s'} \, e^{-S},
\end{aligned}
$$

the dimensionless parameters $\mu$ and $U$ being the analogs of chemical potential and interaction strength. Since $G$ is spin-independent, we omit its spin index. We restrict for convenience to $\mu > 0$ (changing the sign of $\mu$ essentially amounts to the change of variables $\varphi \leftrightarrow \bar{\varphi}$) and to $U < 0$ (as in [8]).

The coefficients of the skeleton series have the analytical expression

$$ \Sigma_{\text{bold}}[G] = \sum_{n=1}^{\infty} a_n G^{2n-1} U^n \quad \text{with} \quad a_n = \frac{(-1)^{n-1} (2n-2)!}{n!(n-1)!}. $$

<footnote n="2">²The naive reasoning goes as follows: If $G_N$ would converge to some $G_\infty$ for $N \to \infty$, then, from (3), one can expect $G_\infty^{-1} = G_0^{-1} - \Sigma_{\text{bold}}[G_\infty]$, and hence, assuming unicity of the solution of the Dyson equation, $G_\infty = G_{\text{exact}}$. Thus $\Sigma_{\text{bold}}[G_{\text{exact}}] = \Sigma_{\text{exact}}$, in contradiction with the misleading convergence of $\Sigma_{\text{bold}}^{(\le N)}[G_{\text{exact}}]$.</footnote>
It is convenient to work with rescaled variables, multiplying propagators with $\sqrt{|U|}$ and dividing self-energies with the same factor,

$$g := G\sqrt{|U|}, \quad \sigma := \Sigma/\sqrt{|U|}. \qquad (4)$$

The rescaled skeleton series is then given by

$$\sigma_{\text{bold}}(g) = \sum_{n=1}^{\infty} \sigma_{\text{bold}}^{(n)}(g) \quad \text{with} \quad \sigma_{\text{bold}}^{(n)}(g) = a_n(-1)^n g^{2n-1}$$

and accordingly $\sigma_{\text{bold}}^{(\le N)}(g) \equiv \sum_{n=1}^{N} \sigma_{\text{bold}}^{(n)}(g)$.

The exact self-energy and propagator are given by

$$\sigma_{\text{exact}}(g_0) = -g_0 \\ g_{\text{exact}}(g_0) = \frac{g_0}{1+g_0^2}$$

in terms of the rescaled free propagator $g_0 := \sqrt{|U|} G_0 = \sqrt{|U|}/\mu$.

If one evaluates the skeleton series at the exact $G$, one obtains the correct physical self-energy for $|U| < \mu^2$ and an incorrect result for $|U| > \mu^2$. This is directly related to the fact that the self-energy functional (which reduces to a function in this toy model) has two branches,

$$\sigma^{(\pm)}(g) = \frac{-1 \pm \sqrt{1-4g^2}}{2g} \qquad (5)$$

as represented in Fig. 1 (this corresponds to the derivative of the two branches of the Luttinger-Ward functional, restricting to spin-independent $G$ for simplicity). The physical branch is the (+) branch for $g_0 < 1$, and the (-) branch for $g_0 > 1$; i.e., $\sigma_{\text{exact}}(g_0) = \sigma^{(\text{sign}(1-g_0))}(g_{\text{exact}}(g_0))$. On the other hand, the skeleton series, evaluated at the exact physical propagator, always converges to the (+) branch; i.e., $\sigma_{\text{bold}}(g_{\text{exact}}(g_0)) = \sigma^{(+)}(g_{\text{exact}}(g_0))$ for all $g_0 > 0$.

Note that $\sigma_{\text{bold}}(g)$ is the expansion of $\sigma^{(+)}(g)$ in powers of $g$, and thus from (5) the convergence radius of the series $\sigma_{\text{bold}}(g)$ is 1/2.

<figure><img src="image_2.png" /><figcaption>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: $\sigma_{\text{bold}}(g) = \sigma^{(+)}(g)$ for $g \le 1/2$.</figcaption></figure>
It is convenient to work with rescaled variables, multiplying propagators with $\sqrt{|U|}$ and dividing self-energies with the same factor,

$$g := G\sqrt{|U|}, \quad \sigma := \Sigma/\sqrt{|U|}. \qquad (4)$$

The rescaled skeleton series is then given by

$$\sigma_{\text{bold}}(g) = \sum_{n=1}^{\infty} \sigma_{\text{bold}}^{(n)}(g) \quad \text{with} \quad \sigma_{\text{bold}}^{(n)}(g) = a_n(-1)^n g^{2n-1}$$

and accordingly $\sigma_{\text{bold}}^{(\le N)}(g) \equiv \sum_{n=1}^{N} \sigma_{\text{bold}}^{(n)}(g)$.

The exact self-energy and propagator are given by

$$\sigma_{\text{exact}}(g_0) = -g_0 \\ g_{\text{exact}}(g_0) = \frac{g_0}{1+g_0^2}$$

in terms of the rescaled free propagator $g_0 := \sqrt{|U|} G_0 = \sqrt{|U|}/\mu$.

If one evaluates the skeleton series at the exact $G$, one obtains the correct physical self-energy for $|U| < \mu^2$ and an incorrect result for $|U| > \mu^2$. This is directly related to the fact that the self-energy functional (which reduces to a function in this toy model) has two branches,

$$\sigma^{(\pm)}(g) = \frac{-1 \pm \sqrt{1-4g^2}}{2g} \qquad (5)$$

as represented in Fig. 1 (this corresponds to the derivative of the two branches of the Luttinger-Ward functional, restricting to spin-independent $G$ for simplicity). The physical branch is the (+) branch for $g_0 < 1$, and the (-) branch for $g_0 > 1$; i.e., $\sigma_{\text{exact}}(g_0) = \sigma^{(\text{sign}(1-g_0))}(g_{\text{exact}}(g_0))$. On the other hand, the skeleton series, evaluated at the exact physical propagator, always converges to the (+) branch; i.e., $\sigma_{\text{bold}}(g_{\text{exact}}(g_0)) = \sigma^{(+)}(g_{\text{exact}}(g_0))$ for all $g_0 > 0$.

Note that $\sigma_{\text{bold}}(g)$ is the expansion of $\sigma^{(+)}(g)$ in powers of $g$, and thus from (5) the convergence radius of the series $\sigma_{\text{bold}}(g)$ is $1/2$.

<figure><img src="image_2.png" /><figcaption>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: $\sigma_{\text{bold}}(g) = \sigma^{(+)}(g)$ for $g \le 1/2$.</figcaption></figure>
## 2.2 Limit of the skeleton sequence

We now go beyond Ref. [8] and study the “skeleton sequence” $G_N$ defined by Eq. (3). Rescaling variables as in (4), in particular setting $g_N := G_N \sqrt{|U|}$, Equation (3) becomes

$$ \frac{1}{g_N} = \frac{1}{g_0} - \sigma_{\text{bold}}^{(\le N)}(g_N). \qquad (6) $$

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

* for $g_0 < 1$, $g_N \underset{N \to \infty}{\to} g_{\text{exact}}(g_0)$

* for $g_0 > 1$, $g_N \underset{N \to \infty}{\to} g_\infty \neq g_{\text{exact}}(g_0)$

i.e., the skeleton sequence converges to the correct physical result below a critical coupling strength, and to an unphysical result above it.

<figure><img src="image_3.png" /><figcaption>Figure 2: Illustrative example of misleading convergence of the skeleton sequence for the toy model. The rescaled propagator $g_N$, obtained from the self-consistency equation with the skeleton series truncated at order $N$, converges for $N \to \infty$ to the limit 0.5, which differs from the exact result (dashed line). This happens when the rescaled free propagator $g_0 > 1$ (here, $g_0 = 1.5$).</figcaption></figure>

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_\infty = 1/2$ of the skeleton sequence $g_N$ is equal to the radius of convergence of the skeleton series $\sigma_{\text{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 skeleton 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

$$ \lim_{N \to \infty} \sigma_{\text{bold}}^{(\le N)}(g_N) \neq \lim_{N \to \infty} \sigma_{\text{bold}}^{(\le N)}(g_\infty) \equiv \sigma_{\text{bold}}(g_\infty). \qquad (7) $$
## 2.2 Limit of the skeleton sequence

We now go beyond Ref. [8] and study the “skeleton sequence” $G_N$ defined by Eq. (3). Rescaling variables as in (4), in particular setting $g_N := G_N \sqrt{|U|}$, Equation (3) becomes

$$ \frac{1}{g_N} = \frac{1}{g_0} - \sigma_{\text{bold}}^{(\le N)}(g_N). \qquad (6) $$

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

* for $g_0 < 1$, $g_N \underset{N \to \infty}{\to} g_{\text{exact}}(g_0)$

* for $g_0 > 1$, $g_N \underset{N \to \infty}{\to} g_\infty \neq g_{\text{exact}}(g_0)$

i.e., the skeleton sequence converges to the correct physical result below a critical coupling strength, and to an unphysical result above it.

<figure><img src="image_3.png" /><figcaption>Figure 2: Illustrative example of misleading convergence of the skeleton sequence for the toy model. The rescaled propagator $g_N$, obtained from the self-consistency equation with the skeleton series truncated at order $N$, converges for $N \to \infty$ to the limit 0.5, which differs from the exact result (dashed line). This happens when the rescaled free propagator $g_0 > 1$ (here, $g_0 = 1.5$).</figcaption></figure>

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_\infty = 1/2$ of the skeleton sequence $g_N$ is equal to the radius of convergence of the skeleton series $\sigma_{\text{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 skeleton 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

$$ \lim_{N \to \infty} \sigma_{\text{bold}}^{(\le N)}(g_N) \neq \lim_{N \to \infty} \sigma_{\text{bold}}^{(\le N)}(g_\infty) \equiv \sigma_{\text{bold}}(g_\infty). \qquad (7) $$
<figure><img src="image_4.png" /><figcaption>Figure 3: Explanation for the misleading convergence. The two branches of the self-energy $\sigma^{(\pm)}(g)$, together with the partial sums of the skeleton series $\sigma_{\text{bold}}^{(\le N)}(g)$ for different values of the truncation order $N$. Also shown is the curve corresponding to the Dyson equation, $-\sigma = 1/g - 1/g_0$. This Dyson-equation curve intersects $\sigma_{\text{bold}}^{(\le N)}(g)$ at $g = g_N$, whereas the exact propagator $g = g_{\text{exact}}$ is given by the intersection of the Dyson-equation curve with the physical branch $\sigma^{(\text{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.</figcaption></figure>

## 2.3 Diagnosing the misleading convergence

In the general case where $G_{\text{exact}}$ is unknown, when one observes numerically that $G_N$ converges to some limit, one needs a way to tell whether this limit is equal to $G_{\text{exact}}$, i.e., whether the result can be trusted. To this end, we consider

$$ \Sigma_{N, \xi} := \sum_{n=1}^{N} \Sigma_{\text{bold}}^{(n)}[G_N] \xi^n. \qquad (8) $$

Assuming that $G_N \to G_\infty$ for $N \to \infty$, the following criterion [2] is a sufficient condition for $G_\infty$ to be equal to $G_{\text{exact}}$:

$$ \left\{ \begin{array}{l} \text{There exists } \epsilon > 0 \text{ such that:} \\ \text{For any } \xi \text{ in the disc } \mathcal{D} = \{|\xi| < 1 + \epsilon\}, \Sigma_{N, \xi} \text{ converges for } N \to \infty; \\ \text{moreover, this sequence is uniformly bounded for } \xi \in \mathcal{D}. \end{array} \right. \qquad (9) $$

The derivation of this criterion is contained in [2], and its main steps are reproduced in the Appendix for convenience.

For all practical purpose, we expect the criterion (9) to be essentially equivalent to the following simpler one:

$$ \text{There exists } \xi > 1 \text{ such that } \Sigma_{N, \xi} \text{ converges for } N \to \infty. \qquad (10) $$

Indeed, (9) implies (10), and a situation where (10) would hold while (9) would not hold seems unlikely to occur. In what follows we will use the simplified criterion (10). We also introduce an extra factor $1/\xi^{N_0}$ in the definition (8) of $\Sigma_{N, \xi}$, 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 $\Sigma_{N, \xi}$ converges).
<figure><img src="image_4.png" /><figcaption>Figure 3: Explanation for the misleading convergence. The two branches of the self-energy $\sigma^{(\pm)}(g)$, together with the partial sums of the skeleton series $\sigma_{\text{bold}}^{(\le N)}(g)$ for different values of the truncation order $N$. Also shown is the curve corresponding to the Dyson equation, $-\sigma = 1/g - 1/g_0$. This Dyson-equation curve intersects $\sigma_{\text{bold}}^{(\le N)}(g)$ at $g = g_N$, whereas the exact propagator $g = g_{\text{exact}}$ is given by the intersection of the Dyson-equation curve with the physical branch $\sigma^{(\text{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.</figcaption></figure>

## 2.3 Diagnosing the misleading convergence

In the general case where $G_{\text{exact}}$ is unknown, when one observes numerically that $G_N$ converges to some limit, one needs a way to tell whether this limit is equal to $G_{\text{exact}}$, i.e., whether the result can be trusted. To this end, we consider

$$ \Sigma_{N, \xi} := \sum_{n=1}^{N} \Sigma_{\text{bold}}^{(n)} [G_N] \xi^n. \qquad (8) $$

Assuming that $G_N \to G_\infty$ for $N \to \infty$, the following criterion [2] is a sufficient condition for $G_\infty$ to be equal to $G_{\text{exact}}$:

$$ \left\{ \begin{array}{l} \text{There exists } \epsilon > 0 \text{ such that:} \\ \text{For any } \xi \text{ in the disc } \mathcal{D} = \{|\xi| < 1 + \epsilon\}, \Sigma_{N, \xi} \text{ converges for } N \to \infty; \\ \text{moreover, this sequence is uniformly bounded for } \xi \in \mathcal{D}. \end{array} \right. \qquad (9) $$

The derivation of this criterion is contained in [2], and its main steps are reproduced in the Appendix for convenience.

For all practical purpose, we expect the criterion (9) to be essentially equivalent to the following simpler one:

$$ \text{There exists } \xi > 1 \text{ such that } \Sigma_{N, \xi} \text{ converges for } N \to \infty. \qquad (10) $$

Indeed, (9) implies (10), and a situation where (10) would hold while (9) would not hold seems unlikely to occur. In what follows we will use the simplified criterion (10). We also introduce an extra factor $1/\xi^{N_0}$ in the definition (8) of $\Sigma_{N, \xi}$, 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 $\Sigma_{N, \xi}$ converges).
For the toy-model, this means that assuming $g_N \to g_\infty$ for $N \to \infty$, a sufficient condition for $g_\infty$ to be equal to the correct physical $g_{\text{exact}}(g_0)$ is that there exists $\xi > 1$ such that

$$ \sigma_{N,\xi} := \sum_{n=1}^{N} \sigma_{\text{bold}}^{(n)}(g_N) \xi^{n-1} = \sigma_{\text{bold}}^{(\le N)}(g_N \sqrt{\xi}) / \sqrt{\xi} $$

converges for $N \to \infty$. As illustrated in Fig. 4, this criterion indeed enables one to detect the misleading convergence for $g_0 > 1$, and to trust the result for $g_0 < 1$.<fnref n="3" />

<figure><img src="image_5.png" /><figcaption>Figure 4: Detecting the misleading convergence for the toy model. Introducing a finite $\xi$, the sequence becomes divergent which enables one to detect the problem (left panel), or remains convergent which enables one to trust the result (right panel).</figcaption></figure>

## 3 Hubbard atom

We turn to the single-site Hubbard model, defined by the grand-canonical Hamiltonian $-\mu \sum_s n_s + U n_\uparrow n_\downarrow$. The propagator can be expressed as a functional integral over $\beta$-antiperiodic Grassmann fields [4, 43],

$$ G_s(\tau) = - \langle \varphi_s(\tau) \bar{\varphi}_s(0) \rangle_s \equiv - \frac{\int \mathcal{D}\varphi \, \mathcal{D}\bar{\varphi} \, \varphi_s(\tau) \bar{\varphi}_s(0) e^{-S}}{\int \mathcal{D}\varphi \, \mathcal{D}\bar{\varphi} \, e^{-S}} \quad (11) $$

with the action

$$ S = \int_0^\beta d\tau \left[ -\sum_s \bar{\varphi}_s(\tau)(G_0^{-1}\varphi_s)(\tau) + U (\bar{\varphi}_\uparrow\bar{\varphi}_\downarrow\varphi_\downarrow\varphi_\uparrow)(\tau) \right] \quad (12) $$

and

$$ G_0^{-1} = \mu - \frac{d}{d\tau}. \quad (13) $$

We restrict for simplicity to the half-filled case $\mu = U/2$, which should be the most dangerous case, since it is at and around half-filling that the misleading convergence of $\Sigma_{\text{bold}}[G_{\text{exact}}]$ was discovered in [1]. We use the BDMC method [6, 35, 44, 45] to sum all skeleton diagrams

<footnote n="3">³For the toy model, the criterion is easily understood from Fig. 3. For $g_0 > 1$, $g_N \to 1/2$ for $N \to \infty$, so that for any $\xi > 1$, $\lim_{N \to \infty} g_N \sqrt{\xi}$ is strictly larger than the convergence radius $1/2$, leading to the divergence of $\sigma_{\text{bold}}^{(\le N)}(g_N \sqrt{\xi})$. On the contrary, for $g_0 < 1$, the $g_N$'s stay at a finite distance on the left of the convergence radius $1/2$.</footnote>
For the toy-model, this means that assuming $g_N \to g_\infty$ for $N \to \infty$, a sufficient condition for $g_\infty$ to be equal to the correct physical $g_{\text{exact}}(g_0)$ is that there exists $\xi > 1$ such that

$$ \sigma_{N,\xi} := \sum_{n=1}^{N} \sigma_{\text{bold}}^{(n)}(g_N) \xi^{n-1} = \sigma_{\text{bold}}^{(\le N)}(g_N \sqrt{\xi}) / \sqrt{\xi} $$

converges for $N \to \infty$. As illustrated in Fig. 4, this criterion indeed enables one to detect the misleading convergence for $g_0 > 1$, and to trust the result for $g_0 < 1$.<fnref n="3" />

<figure><img src="image_5.png" /><figcaption>Figure 4: Detecting the misleading convergence for the toy model. Introducing a finite $\xi$, the sequence becomes divergent which enables one to detect the problem (left panel), or remains convergent which enables one to trust the result (right panel).</figcaption></figure>

### 3 Hubbard atom

We turn to the single-site Hubbard model, defined by the grand-canonical Hamiltonian $-\mu \sum_s n_s + U n_\uparrow n_\downarrow$. The propagator can be expressed as a functional integral over $\beta$-antiperiodic Grassmann fields [4, 43],

$$ G_s(\tau) = - \langle \varphi_s(\tau) \bar{\varphi}_s(0) \rangle_s = - \frac{\int \mathcal{D}\varphi \, \mathcal{D}\bar{\varphi} \, \varphi_s(\tau) \bar{\varphi}_s(0) e^{-S}}{\int \mathcal{D}\varphi \, \mathcal{D}\bar{\varphi} \, e^{-S}} \quad (11) $$

with the action

$$ S = \int_0^\beta d\tau \left[ -\sum_s \bar{\varphi}_s(\tau)(G_0^{-1}\varphi_s)(\tau) + U (\bar{\varphi}_\uparrow\bar{\varphi}_\downarrow\varphi_\downarrow\varphi_\uparrow)(\tau) \right] \quad (12) $$

and

$$ G_0^{-1} = \mu - \frac{d}{d\tau}. \quad (13) $$

We restrict for simplicity to the half-filled case $\mu = U/2$, which should be the most dangerous case, since it is at and around half-filling that the misleading convergence of $\Sigma_{\text{bold}}[G_{\text{exact}}]$ was discovered in [1]. We use the BDMC method [6, 35, 44, 45] to sum all skeleton diagrams

<footnote n="3">³For the toy model, the criterion is easily understood from Fig. 3. For $g_0 > 1$, $g_N \to 1/2$ for $N \to \infty$, so that for any $\xi > 1$, $\lim_{N \to \infty} g_N \sqrt{\xi}$ is strictly larger than the convergence radius $1/2$, leading to the divergence of $\sigma_{\text{bold}}^{(\le N)}(g_N \sqrt{\xi})$. On the contrary, for $g_0 < 1$, the $g_N$'s stay at a finite distance on the left of the convergence radius $1/2$.</footnote>
and solve the self-consistency equation (3) for truncation orders $N \le 8$ (note that at half filling, $\Sigma_{\text{bold}}^{(n)} = 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 $\Sigma_{\text{bold}}^{(\le N)}[G_N] =: \Sigma_N$ can converge to an unphysical result. Let us first consider the double occupancy

$$D = \langle n_{\uparrow} n_{\downarrow} \rangle = U^{-1} \text{tr}(\Sigma G)$$

and the corresponding sequence $D_N := U^{-1} \text{tr}(\Sigma_N G_N)$. At large enough $U$, our data strongly indicate that this sequence does converge (albeit slowly) towards an unphysical result, see left panel of Fig. 5. For small enough $U$, there is a fast convergence to the correct result, see right panel of Fig. 5.

<figure><img src="image_6.png" /><figcaption>Figure 5: For the Hubbard atom at half filling, the double occupancy, as obtained from the skeleton sequence, converges to an unphysical result for large $U$ (left panel) and to the correct result for small enough $U$ (right panel) when the truncation order $N \to \infty$ (dashed line: exact result).</figcaption></figure>

The next question is whether the criterion (10) enables us to discriminate between these two situations. We therefore plot the sequence $\Sigma_{N,\xi}$ in Figs. 6 and 7. We only show the imaginary part because in the considered half-filled case, the real part of $\Sigma_N(\omega_n)$ automatically equals $U/2$; moreover we focus for simplicity on the lowest Matsubara frequency $\omega_0 = \pi/\beta$, and we choose $N_0 = 2$.

For $\xi = 1$, $\Sigma_{N,\xi}$ reduces to the original skeleton sequence $\Sigma_N$, and the behavior is similar to the double occupancy: The sequence appears to converge, albeit slowly, towards an unphysical result for $\beta U = 8$ (Fig. 6), while fast convergence to the correct physical result takes place for $\beta U = 1$ (Fig. 7). For $\xi > 1$, the sequence does not appear to converge any more for $\beta U = 8$, see Fig. 6: The data do not satisfy the criterion, indicating that the results cannot be trusted in this case. In contrast, for $\beta U = 1$, the criterion enables one to validate the results, since the sequence remains convergent at $\xi > 1$, see Fig. 7.

Regarding the choice of $\xi$, 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. More precisely, $\xi - 1$ should not be too small, so that $\xi^N$ differs significantly from 1 (and hence $\Sigma_{N,\xi}$ differ significantly from $\Sigma_{N,1} = \Sigma_N$) at the largest accessible order $N_{\text{max}}$. This necessity to work with a finite $\xi - 1$ implies that the criterion is conservative: It leads to discarding results in a region of the parameter space near but outside the misleading-convergence region.<fnref n="4" /> For $N_{\text{max}} = 8$, the choices $\xi - 1 = 0.1$ and $0.2$ are a priori large enough (since $\xi^8 \approx 2$ and 4) and Fig. 6 confirms

<footnote n="4">⁴The required $\xi - 1$ scales as $1/N_{\text{max}}$; for the toy model this leads to discarding results not only in the misleading-convergence region $g_0 > 1$, but also for $0 \le 1 - g_0 \lesssim 1/N_{\text{max}}$.</footnote>
and solve the self-consistency equation (3) for truncation orders $N \le 8$ (note that at half filling, $\Sigma_{\text{bold}}^{(n)} = 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 $\Sigma_{\text{bold}}^{(\le N)}[G_N] =: \Sigma_N$ can converge to an unphysical result. Let us first consider the double occupancy

$$D = \langle n_{\uparrow} n_{\downarrow} \rangle = U^{-1} \text{tr}(\Sigma G)$$

and the corresponding sequence $D_N := U^{-1} \text{tr}(\Sigma_N G_N)$. At large enough $U$, our data strongly indicate that this sequence does converge (albeit slowly) towards an unphysical result, see left panel of Fig. 5. For small enough $U$, there is a fast convergence to the correct result, see right panel of Fig. 5.

<figure><img src="image_6.png" /><figcaption>Figure 5: For the Hubbard atom at half filling, the double occupancy, as obtained from the skeleton sequence, converges to an unphysical result for large $U$ (left panel) and to the correct result for small enough $U$ (right panel) when the truncation order $N \to \infty$ (dashed line: exact result).</figcaption></figure>

The next question is whether the criterion (10) enables us to discriminate between these two situations. We therefore plot the sequence $\Sigma_{N,\xi}$ in Figs. 6 and 7. We only show the imaginary part because in the considered half-filled case, the real part of $\Sigma_N(\omega_n)$ automatically equals $U/2$; moreover we focus for simplicity on the lowest Matsubara frequency $\omega_0 = \pi/\beta$, and we choose $N_0 = 2$.

For $\xi = 1$, $\Sigma_{N,\xi}$ reduces to the original skeleton sequence $\Sigma_N$, and the behavior is similar to the double occupancy: The sequence appears to converge, albeit slowly, towards an unphysical result for $\beta U = 8$ (Fig. 6), while fast convergence to the correct physical result takes place for $\beta U = 1$ (Fig. 7). For $\xi > 1$, the sequence does not appear to converge any more for $\beta U = 8$, see Fig. 6: The data do not satisfy the criterion, indicating that the results cannot be trusted in this case. In contrast, for $\beta U = 1$, the criterion enables one to validate the results, since the sequence remains convergent at $\xi > 1$, see Fig. 7.

Regarding the choice of $\xi$, 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. More precisely, $\xi - 1$ should not be too small, so that $\xi^N$ differs significantly from 1 (and hence $\Sigma_{N,\xi}$ differ significantly from $\Sigma_{N,1} = \Sigma_N$) at the largest accessible order $N_{\text{max}}$. This necessity to work with a finite $\xi - 1$ implies that the criterion is conservative: It leads to discarding results in a region of the parameter space near but outside the misleading-convergence region.<fnref n="4" /> For $N_{\text{max}} = 8$, the choices $\xi - 1 = 0.1$ and $0.2$ are a priori large enough (since $\xi^8 \approx 2$ and 4) and Fig. 6 confirms

<footnote n="4">⁴The required $\xi - 1$ scales as $1/N_{\text{max}}$; for the toy model this leads to discarding results not only in the misleading-convergence region $g_0 > 1$, but also for $0 \le 1 - g_0 \lesssim 1/N_{\text{max}}$.</footnote>
that they allow to detect the divergence of $\Sigma_{N,\xi}$ in the misleading-convergence regime, while on the other hand Fig. 7 shows that $\beta U = 1$ is at a sufficient distance from the misleading-convergence region for $\Sigma_{N,\xi}$ to remain convergent for these $\xi$ values.

<figure><img src="image_7.png" /><figcaption>Figure 6: For the half-filled Hubbard atom at large coupling, the original skeleton sequence ($\xi = 1$) converges to an unphysical result. At $\xi > 1$, the sequence does not converge any more: The data do not satisfy the criterion.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: For the half-filled Hubbard atom at small enough coupling, the original skeleton sequence ($\xi = 1$) converges to the correct physical result. At $\xi > 1$, the sequence remains convergent: The criterion enables one to trust the result.</figcaption></figure>

We remark that when solving Eq. (3) by iterations, in the case where the convergence to the unphysical result for $N \to \infty$ occurs, convergence as a function of iterations at fixed $N$ only takes place if we use a damping procedure, where $G_N$ at iteration $(i+1)$ is obtained as $G_N^{(i+1)} = [G_0^{-1} - \Sigma^{(i)}]^{-1}$ with $\Sigma^{(i)}$ a weighted average of $\Sigma_{\text{bold}}^{(\le N)}[G_N^{(i)}]$ and $\Sigma^{(i-1)}$, while the fixed point is unstable for the undamped iterative procedure $\Sigma^{(i)} := \Sigma_{\text{bold}}^{(\le N)}[G_N^{(i)}]$. Such a damping procedure is commonly used in BDMC where it also reduces the statistical error [45, 46]. In the toy model, one can easily show that an increasingly strong damping is required when $N$ is increased, because for $N \to \infty$, the slope $[d\sigma_{\text{bold}}^{(\le N)}(g)/dg]_{g=g_N}$ diverges, making the undamped iterative procedure unstable. This observation could also be useful for misleading-convergence detection.
that they allow to detect the divergence of $\Sigma_{N,\xi}$ in the misleading-convergence regime, while on the other hand Fig. 7 shows that $\beta U = 1$ is at a sufficient distance from the misleading-convergence region for $\Sigma_{N,\xi}$ to remain convergent for these $\xi$ values.

<figure><img src="image_7.png" /><figcaption>Figure 6: For the half-filled Hubbard atom at large coupling, the original skeleton sequence ($\xi = 1$) converges to an unphysical result. At $\xi > 1$, the sequence does not converge any more: The data do not satisfy the criterion.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: For the half-filled Hubbard atom at small enough coupling, the original skeleton sequence ($\xi = 1$) converges to the correct physical result. At $\xi > 1$, the sequence remains convergent: The criterion enables one to trust the result.</figcaption></figure>

We remark that when solving Eq. (3) by iterations, in the case where the convergence to the unphysical result for $N \to \infty$ occurs, convergence as a function of iterations at fixed $N$ only takes place if we use a damping procedure, where $G_N$ at iteration $(i+1)$ is obtained as $G_N^{(i+1)} = [G_0^{-1} - \Sigma^{(i)}]^{-1}$ with $\Sigma^{(i)}$ a weighted average of $\Sigma_{\text{bold}}^{(\le N)}[G_N^{(i)}]$ and $\Sigma^{(i-1)}$, while the fixed point is unstable for the undamped iterative procedure $\Sigma^{(i)} := \Sigma_{\text{bold}}^{(\le N)}[G_N^{(i)}]$. Such a damping procedure is commonly used in BDMC where it also reduces the statistical error [45, 46]. In the toy model, one can easily show that an increasingly strong damping is required when $N$ is increased, because for $N \to \infty$, the slope $[d\sigma_{\text{bold}}^{(\le N)}(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 [1], 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. (11,12) are viewed as a mapping G₀ ↦ G[G₀], then there exists G₀,unphys such that Σ<sub>bold</sub>[G<sub>exact</sub>] = G₀,unphys⁻¹ − G<sub>exact</sub>⁻¹ and G[G₀,unphys] = G<sub>exact</sub> ≡ G[G₀]. As noted in [1], this G₀,unphys does not belong to the set of physical bare propagators which are of the form (13) for some value of chemical potential; therefore, by looking at G₀,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 sequence G<sub>N</sub> reported here cannot be detected in this way. Indeed, the self-consistency equation (3) is enforced with the original physical G₀.

# 4 Conclusion

We have demonstrated 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 the findings of [1] combined with an additional subtle mathematical mechanism which we have elucidated by analyzing the zero space-time dimensional model of [8]. In this problematic regime, lowest order calculations can be off by one order of magnitude, but access to higher orders enables one to detect the problem numerically through the divergence of a slightly modified sequence, whereas seeing convergence of this modified sequence enables one to rule out misleading convergence and to trust the result, as proposed in [2] 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 [47–68] (for which it was shown that misleading convergence generically does not occur [2]) and with strong-coupling expansions [69–71], self-consistent BDMC remains among the state of the art approaches. In particular, during finalization of the present manuscript, its main findings have been used 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 [72].

# Acknowledgments

We thank N. Prokof'ev, B. Svistunov and L. Reining for useful discussions and comments.

**Funding information** We acknowledge support from ERC (Critisup2, H2020 Adv-743159) (F.W.), ANR (LODIS, ANR-21-CE30-0033) (K.V.H. and F.W.), the Simons Foundation through the Simons Collaboration on the Many Electron Problem and EPSRC (EP/P003052/1) (E.K.), the National Natural Science Foundation of China (grant No. 11625522) and the Science and Technology Committee of Shanghai (grant No. 20DZ2210100) (Y.D.). The Flatiron Institute is a division of the Simons Foundation.
Finally, we comment on the link with the multivaluedness of the self-energy functional Σ[G] (i.e., of the Luttinger-Ward functional). In [1], 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. (11,12) are viewed as a mapping $G_0 \mapsto G[G_0]$, then there exists $G_{0,\text{unphys}}$ such that $\Sigma_{\text{bold}}[G_{\text{exact}}] = G_{0,\text{unphys}}^{-1} - G_{\text{exact}}^{-1}$ and $G[G_{0,\text{unphys}}] = G_{\text{exact}} \equiv G[G_0]$. As noted in [1], this $G_{0,\text{unphys}}$ does not belong to the set of physical bare propagators which are of the form (13) for some value of chemical potential; therefore, by looking at $G_{0,\text{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 sequence $G_N$ reported here cannot be detected in this way. Indeed, the self-consistency equation (3) is enforced with the original physical $G_0$.

## 4 Conclusion

We have demonstrated 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 the findings of [1] combined with an additional subtle mathematical mechanism which we have elucidated by analyzing the zero space-time dimensional model of [8]. In this problematic regime, lowest order calculations can be off by one order of magnitude, but access to higher orders enables one to detect the problem numerically through the divergence of a slightly modified sequence, whereas seeing convergence of this modified sequence enables one to rule out misleading convergence and to trust the result, as proposed in [2] 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 [47–68] (for which it was shown that misleading convergence generically does not occur [2]) and with strong-coupling expansions [69–71], self-consistent BDMC remains among the state of the art approaches. In particular, during finalization of the present manuscript, its main findings have been used 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 [72].

## Acknowledgments

We thank N. Prokof'ev, B. Svistunov and L. Reining for useful discussions and comments.

**Funding information** We acknowledge support from ERC (Critisup2, H2020 Adv-743159) (F.W.), ANR (LODIS, ANR-21-CE30-0033) (K.V.H. and F.W.), the Simons Foundation through the Simons Collaboration on the Many Electron Problem and EPSRC (EP/P003052/1) (E.K.), the National Natural Science Foundation of China (grant No. 11625522) and the Science and Technology Committee of Shanghai (grant No. 20DZ2210100) (Y.D.). The Flatiron Institute is a division of the Simons Foundation.
## Appendix: Main steps of the derivation of the criterion

For convenience, we reproduce here the main steps of the derivation of the criterion (9), see [2, 73] for more details. For definiteness, we consider the Hubbard model at finite temperature; the reasoning is also directly applicable to the Hubbard atom (by removing the position variable) and to the zero space-time dimensional toy model (by also removing the imaginary time variable). Let $\Sigma_{\infty, \xi} := \lim_{N \to \infty} \Sigma_{N, \xi}$. By making use of Morera's theorem, Cauchy's integral formula and the dominated convergence theorem, the condition (9) allows one to show the key property

$$ \Sigma_{\infty, \xi} = \sum_{n=1}^{\infty} \Sigma_{\text{bold}}^{(n)}[G_{\infty}] \xi^n, \quad \forall \xi \in \mathcal{D}. \qquad (14) $$

Setting $\xi = 1$ in (14) yields<fnref n="5" />

$$ \lim_{N \to \infty} \Sigma_{\text{bold}}^{(\le N)}[G_N] = \Sigma_{\text{bold}}[G_{\infty}]. \qquad (15) $$

Substituting (15) into (3) yields

$$ G_{\infty}^{-1} = G_{0}^{-1} - \Sigma_{\text{bold}}[G_{\infty}]. \qquad (16) $$

The next step is to consider the action

$$ S_{\text{bold}}^{(\xi)} := - \sum_{\mathbf{r},s} \int_{0}^{\beta} d\tau \left[ \bar{\varphi}_{s} \left( G_{\infty}^{-1} + \sum_{n=1}^{\infty} \Sigma_{\text{bold}}^{(n)}[G_{\infty}] \xi^{n} \right) \varphi_{s} \right](\mathbf{r}, \tau) \\ + \xi U \sum_{\mathbf{r}} \int_{0}^{\beta} d\tau (\bar{\varphi}_{\uparrow} \bar{\varphi}_{\downarrow} \varphi_{\downarrow} \varphi_{\uparrow})(\mathbf{r}, \tau) $$

and the corresponding propagator $G_{\text{bold}}^{(\xi)}(\mathbf{r}, \tau) := -\langle \varphi_s(\mathbf{r}, \tau) \bar{\varphi}_s(\mathbf{0}, 0) \rangle_{S_{\text{bold}}^{(\xi)}}$. The action $S_{\text{bold}}^{(\xi)}$ is designed in such a way that

$$ \left. \frac{\partial^n G_{\text{bold}}^{(\xi)}}{\partial \xi^n} \right|_{\xi=0} = 0, \quad \forall n \ge 1. \qquad (17) $$

Obviously, $G_{\text{bold}}^{(\xi=0)} = G_{\infty}$. On the other hand, (16) implies that $S_{\text{bold}}^{(\xi=1)}$ is equal to the physical action of the Hubbard model, so that $G_{\text{bold}}^{(\xi=1)} = G_{\text{exact}}$. Now, since $S_{\text{bold}}^{(\xi)}$ depends analytically on $\xi$ in $\mathcal{D}$, we expect (at least for fermions on a lattice at finite temperature) that one of the following alternatives holds:

(i) $G_{\text{bold}}^{(\xi)}$ depends analytically on $\xi$ in $\mathcal{D}$

(ii) $G_{\text{bold}}^{(\xi)}$ has a non-removable singularity at a point $\xi_c \in \mathcal{D}$ (e.g., a phase transition), and $G_{\text{bold}}^{(\xi)}$ is analytic in the disc $\{|\xi| < |\xi_c|\}$.

In case (ii), the convergence radius of the Taylor series of $G_{\text{bold}}^{(\xi)}$ at the origin would be $|\xi_c|$, in contradiction with (17). Hence (i) holds, and we have

$$ G_{\text{exact}} = G_{\text{bold}}^{(\xi=1)} = G_{\text{bold}}^{(\xi=0)} + \sum_{n=1}^{\infty} \left. \frac{1}{n!} \frac{\partial^n G_{\text{bold}}^{(\xi)}}{\partial \xi^n} \right|_{\xi=0} = G_{\text{bold}}^{(\xi=0)} = G_{\infty}. $$

This concludes the derivation of the equality $G_{\infty} = G_{\text{exact}}$ under the assumption (9).

<footnote n="5">⁵Equation (15) breaks down for the toy model in the problematic regime ($g_0 > 1$), as pointed out in Eq. (7).</footnote>
## Appendix: Main steps of the derivation of the criterion

For convenience, we reproduce here the main steps of the derivation of the criterion (9), see [2, 73] for more details. For definiteness, we consider the Hubbard model at finite temperature; the reasoning is also directly applicable to the Hubbard atom (by removing the position variable) and to the zero space-time dimensional toy model (by also removing the imaginary time variable). Let $\Sigma_{\infty, \xi} := \lim_{N \to \infty} \Sigma_{N, \xi}$. By making use of Morera's theorem, Cauchy's integral formula and the dominated convergence theorem, the condition (9) allows one to show the key property

$$ \Sigma_{\infty, \xi} = \sum_{n=1}^{\infty} \Sigma_{\text{bold}}^{(n)}[G_{\infty}] \xi^n, \quad \forall \xi \in \mathcal{D}. \qquad (14) $$

Setting $\xi = 1$ in (14) yields<fnref n="5" />

$$ \lim_{N \to \infty} \Sigma_{\text{bold}}^{(\le N)}[G_N] = \Sigma_{\text{bold}}[G_{\infty}]. \qquad (15) $$

Substituting (15) into (3) yields

$$ G_{\infty}^{-1} = G_{0}^{-1} - \Sigma_{\text{bold}}[G_{\infty}]. \qquad (16) $$

The next step is to consider the action

$$ S_{\text{bold}}^{(\xi)} := - \sum_{\mathbf{r},s} \int_0^{\beta} d\tau \left[ \bar{\varphi}_s \left( G_{\infty}^{-1} + \sum_{n=1}^{\infty} \Sigma_{\text{bold}}^{(n)}[G_{\infty}] \xi^n \right) \varphi_s \right](\mathbf{r}, \tau) \\ + \xi U \sum_{\mathbf{r}} \int_0^{\beta} d\tau (\bar{\varphi}_{\uparrow} \bar{\varphi}_{\downarrow} \varphi_{\downarrow} \varphi_{\uparrow})(\mathbf{r}, \tau) $$

and the corresponding propagator $G_{\text{bold}}^{(\xi)}(\mathbf{r}, \tau) := -\langle \varphi_s(\mathbf{r}, \tau) \bar{\varphi}_s(\mathbf{0}, 0) \rangle_{S_{\text{bold}}^{(\xi)}}$. The action $S_{\text{bold}}^{(\xi)}$ is designed in such a way that

$$ \left. \frac{\partial^n G_{\text{bold}}^{(\xi)}}{\partial \xi^n} \right|_{\xi=0} = 0, \quad \forall n \ge 1. \qquad (17) $$

Obviously, $G_{\text{bold}}^{(\xi=0)} = G_{\infty}$. On the other hand, (16) implies that $S_{\text{bold}}^{(\xi=1)}$ is equal to the physical action of the Hubbard model, so that $G_{\text{bold}}^{(\xi=1)} = G_{\text{exact}}$. Now, since $S_{\text{bold}}^{(\xi)}$ depends analytically on $\xi$ in $\mathcal{D}$, we expect (at least for fermions on a lattice at finite temperature) that one of the following alternatives holds:

(i) $G_{\text{bold}}^{(\xi)}$ depends analytically on $\xi$ in $\mathcal{D}$

(ii) $G_{\text{bold}}^{(\xi)}$ has a non-removable singularity at a point $\xi_c \in \mathcal{D}$ (e.g., a phase transition), and $G_{\text{bold}}^{(\xi)}$ is analytic in the disc $\{|\xi| < |\xi_c|\}$.

In case (ii), the convergence radius of the Taylor series of $G_{\text{bold}}^{(\xi)}$ at the origin would be $|\xi_c|$, in contradiction with (17). Hence (i) holds, and we have

$$ G_{\text{exact}} = G_{\text{bold}}^{(\xi=1)} = G_{\text{bold}}^{(\xi=0)} + \sum_{n=1}^{\infty} \left. \frac{1}{n!} \frac{\partial^n G_{\text{bold}}^{(\xi)}}{\partial \xi^n} \right|_{\xi=0} = G_{\text{bold}}^{(\xi=0)} = G_{\infty}. $$

This concludes the derivation of the equality $G_{\infty} = G_{\text{exact}}$ under the assumption (9).

<footnote n="5">⁵Equation (15) breaks down for the toy model in the problematic regime ($g_0 > 1$), as pointed out in Eq. (7).</footnote>
References

[1] E. Kozik, M. Ferrero and A. Georges, *Nonexistence of the Luttinger-Ward Functional and Misleading Convergence of Skeleton Diagrammatic Series for Hubbard-Like Models*, Phys. Rev. Lett. **114**, 156402 (2015), doi:10.1103/PhysRevLett.114.156402.

[2] R. Rossi, F. Werner, N. Prokof'ev and B. Svistunov, *Shifted-Action Expansion and Applicability of Dressed Diagrammatic Schemes*, Phys. Rev. B **93**, 161102(R) (2016), doi:10.1103/PhysRevB.93.161102.

[3] R. M. Martin, L. Reining and D. M. Ceperley, *Interacting Electrons: Theory and Computational Approaches*, Cambridge University Press, ISBN 9781139050807, doi:10.1017/CBO9781139050807 (2016).

[4] N. Dupuis, *Field Theory of Condensed Matter and Ultracold Gases*, World Scientific, ISBN 978-1-80061-390-4, doi:10.1142/q0409 (2023).

[5] G. Stefanucci and R. van Leeuwen, *Nonequilibrium Many-Body Theory of Quantum Systems: A Modern Introduction*, Cambridge University Press, ISBN 9781139023979, doi:10.1017/CBO9781139023979 (2013).

[6] K. Van Houcke, E. Kozik, N. Prokof'ev and B. Svistunov, *Diagrammatic Monte Carlo*, in Computer Simulation Studies in Condensed Matter Physics XXI. CSP-2008. Eds. D.P. Landau, S.P. Lewis, and H.B. Schüttler, Phys. Procedia **6**, 95 (2010), doi:10.1016/j.phpro.2010.09.034.

[7] A. Stan, P. Romaniello, S. Rigamonti, L. Reining and J. Berger, *Unphysical and physical solutions in many-body theories: from weak to strong correlation*, New J. Phys. **17**, 093045 (2015), doi:10.1088/1367-2630/17/9/093045.

[8] R. Rossi and F. Werner, *Skeleton series and multivaluedness of the self-energy functional in zero space-time dimensions*, J. Phys. A **48**, 485202 (2015), doi:10.1088/1751-8113/48/48/485202.

[9] T. Schäfer, S. Ciuchi, M. Wallerberger, P. Thunström, O. Gunnarsson, G. Sangiovanni, G. Rohringer and A. Toschi, *Nonperturbative landscape of the Mott-Hubbard transition: Multiple divergence lines around the critical endpoint*, Phys. Rev. B **94**, 235108 (2016), doi:10.1103/PhysRevB.94.235108.

[10] W. Tarantino, P. Romaniello, J. A. Berger and L. Reining, *Self-consistent Dyson equation and self-energy functionals: An analysis and illustration on the example of the Hubbard atom*, Phys. Rev. B **96**, 045124 (2017), doi:10.1103/PhysRevB.96.045124.

[11] O. Gunnarsson, G. Rohringer, T. Schäfer, G. Sangiovanni and A. Toschi, *Breakdown of Traditional Many-Body Theories for Correlated Electrons*, Phys. Rev. Lett. **119**, 056402 (2017), doi:10.1103/PhysRevLett.119.056402.

[12] J. Vučičević, N. Wentzell, M. Ferrero and O. Parcollet, *Practical consequences of the Luttinger-Ward functional multivaluedness for cluster DMFT methods*, Phys. Rev. B **97**, 125141 (2018), doi:10.1103/PhysRevB.97.125141.

[13] L. Lin and M. Lindsey, *Variational structure of Luttinger-Ward formalism and bold diagrammatic expansion for Euclidean lattice field theory*, PNAS **115**(10), 2282 (2018), doi:10.1073/pnas.1720782115.
References

[1] E. Kozik, M. Ferrero and A. Georges, *Nonexistence of the Luttinger-Ward Functional and Misleading Convergence of Skeleton Diagrammatic Series for Hubbard-Like Models*, Phys. Rev. Lett. **114**, 156402 (2015), doi:10.1103/PhysRevLett.114.156402.

[2] R. Rossi, F. Werner, N. Prokof'ev and B. Svistunov, *Shifted-Action Expansion and Applicability of Dressed Diagrammatic Schemes*, Phys. Rev. B **93**, 161102(R) (2016), doi:10.1103/PhysRevB.93.161102.

[3] R. M. Martin, L. Reining and D. M. Ceperley, *Interacting Electrons: Theory and Computational Approaches*, Cambridge University Press, ISBN 9781139050807, doi:10.1017/CBO9781139050807 (2016).

[4] N. Dupuis, *Field Theory of Condensed Matter and Ultracold Gases*, World Scientific, ISBN 978-1-80061-390-4, doi:10.1142/q0409 (2023).

[5] G. Stefanucci and R. van Leeuwen, *Nonequilibrium Many-Body Theory of Quantum Systems: A Modern Introduction*, Cambridge University Press, ISBN 9781139023979, doi:10.1017/CBO9781139023979 (2013).

[6] K. Van Houcke, E. Kozik, N. Prokof'ev and B. Svistunov, *Diagrammatic Monte Carlo*, in Computer Simulation Studies in Condensed Matter Physics XXI. CSP-2008. Eds. D.P. Landau, S.P. Lewis, and H.B. Schüttler, Phys. Procedia **6**, 95 (2010), doi:10.1016/j.phpro.2010.09.034.

[7] A. Stan, P. Romaniello, S. Rigamonti, L. Reining and J. Berger, *Unphysical and physical solutions in many-body theories: from weak to strong correlation*, New J. Phys. **17**, 093045 (2015), doi:10.1088/1367-2630/17/9/093045.

[8] R. Rossi and F. Werner, *Skeleton series and multivaluedness of the self-energy functional in zero space-time dimensions*, J. Phys. A **48**, 485202 (2015), doi:10.1088/1751-8113/48/48/485202.

[9] T. Schäfer, S. Ciuchi, M. Wallerberger, P. Thunström, O. Gunnarsson, G. Sangiovanni, G. Rohringer and A. Toschi, *Nonperturbative landscape of the Mott-Hubbard transition: Multiple divergence lines around the critical endpoint*, Phys. Rev. B **94**, 235108 (2016), doi:10.1103/PhysRevB.94.235108.

[10] W. Tarantino, P. Romaniello, J. A. Berger and L. Reining, *Self-consistent Dyson equation and self-energy functionals: An analysis and illustration on the example of the Hubbard atom*, Phys. Rev. B **96**, 045124 (2017), doi:10.1103/PhysRevB.96.045124.

[11] O. Gunnarsson, G. Rohringer, T. Schäfer, G. Sangiovanni and A. Toschi, *Breakdown of Traditional Many-Body Theories for Correlated Electrons*, Phys. Rev. Lett. **119**, 056402 (2017), doi:10.1103/PhysRevLett.119.056402.

[12] J. Vučičević, N. Wentzell, M. Ferrero and O. Parcollet, *Practical consequences of the Luttinger-Ward functional multivaluedness for cluster DMFT methods*, Phys. Rev. B **97**, 125141 (2018), doi:10.1103/PhysRevB.97.125141.

[13] L. Lin and M. Lindsey, *Variational structure of Luttinger-Ward formalism and bold diagrammatic expansion for Euclidean lattice field theory*, PNAS **115**(10), 2282 (2018), doi:10.1073/pnas.1720782115.
[14] P. Chalupa, P. Gunacker, T. Schäfer, K. Held and A. Toschi, *Divergences of the irreducible vertex functions in correlated metallic systems: Insights from the Anderson impurity model*, Phys. Rev. B **97**, 245136 (2018), doi:10.1103/PhysRevB.97.245136.

[15] A. J. Kim and V. Sacksteder, *Multivaluedness of the Luttinger-Ward functional in the fermionic and bosonic system with replicas*, Phys. Rev. B **101**, 115146 (2020), doi:10.1103/PhysRevB.101.115146.

[16] T. Schäfer, G. Rohringer, O. Gunnarsson, S. Ciuchi, G. Sangiovanni and A. Toschi, *Divergent Precursors of the Mott-Hubbard Transition at the Two-Particle Level*, Phys. Rev. Lett. **110**, 246405 (2013), doi:10.1103/PhysRevLett.110.246405.

[17] O. Gunnarsson, T. Schäfer, J. P. F. LeBlanc, J. Merino, G. Sangiovanni, G. Rohringer and A. Toschi, *Parquet decomposition calculations of the electronic self-energy*, Phys. Rev. B **93**, 245102 (2016), doi:10.1103/PhysRevB.93.245102.

[18] G. Rohringer, H. Hafermann, A. Toschi, A. A. Katanin, A. E. Antipov, M. I. Katsnelson, A. I. Lichtenstein, A. N. Rubtsov and K. Held, *Diagrammatic routes to nonlocal correlations beyond dynamical mean field theory*, Rev. Mod. Phys. **90**, 025003 (2018), doi:10.1103/RevModPhys.90.025003.

[19] M. Reitner, P. Chalupa, L. Del Re, D. Springer, S. Ciuchi, G. Sangiovanni and A. Toschi, *Attractive Effect of a Strong Electronic Repulsion: The Physics of Vertex Divergences*, Phys. Rev. Lett. **125**, 196403 (2020), doi:10.1103/PhysRevLett.125.196403.

[20] S. Adler, F. Krien, P. Chalupa-Gantner, G. Sangiovanni and A. Toschi, *Non-perturbative intertwining between spin and charge correlations: A “smoking gun” single-boson-exchange result*, SciPost Phys. **16**, 054 (2024), doi:10.21468/SciPostPhys.16.2.054.

[21] C. Naya, T. Bertolini and J. Carlström, *Stability of line-node semimetals with strong Coulomb interactions and properties of the symmetry-broken state*, Phys. Rev. Res. **5**, 013069 (2023), doi:10.1103/PhysRevResearch.5.013069.

[22] F. Aryasetiawan and O. Gunnarsson, *The GW method*, Rep. Prog. Phys. **61**(3), 237 (1998), doi:10.1088/0034-4885/61/3/002.

[23] K. T. Williams, Y. Yao, J. Li, L. Chen, H. Shi, M. Motta, C. Niu, U. Ray, S. Guo, R. J. Anderson, J. Li, L. N. Tran et al., *Direct Comparison of Many-Body Methods for Realistic Electronic Hamiltonians*, Phys. Rev. X **10**, 011041 (2020), doi:10.1103/PhysRevX.10.011041.

[24] N. E. Dahlen and R. van Leeuwen, *Self-consistent solution of the Dyson equation for atoms and molecules within a conserving approximation*, J. Chem. Phys. **122**(16), 164102 (2005), doi:10.1063/1.1884965.

[25] J. J. Phillips and D. Zgid, *The description of strong correlation within self-consistent Green's function second-order perturbation theory*, J. Chem. Phys. **140**(24), 241101 (2014), doi:10.1063/1.4884951.

[26] I. S. Tupitsyn, A. S. Mishchenko, N. Nagaosa and N. Prokof'ev, *Coulomb and electron-phonon interactions in metals*, Phys. Rev. B **94**, 155145 (2016), doi:10.1103/PhysRevB.94.155145.

[27] I. S. Tupitsyn and N. V. Prokof'ev, *Stability of Dirac Liquids with Strong Coulomb Interaction*, Phys. Rev. Lett. **118**, 026403 (2017), doi:10.1103/PhysRevLett.118.026403.
[14] P. Chalupa, P. Gunacker, T. Schäfer, K. Held and A. Toschi, *Divergences of the irreducible vertex functions in correlated metallic systems: Insights from the Anderson impurity model*, Phys. Rev. B **97**, 245136 (2018), doi:10.1103/PhysRevB.97.245136.

[15] A. J. Kim and V. Sacksteder, *Multivaluedness of the Luttinger-Ward functional in the fermionic and bosonic system with replicas*, Phys. Rev. B **101**, 115146 (2020), doi:10.1103/PhysRevB.101.115146.

[16] T. Schäfer, G. Rohringer, O. Gunnarsson, S. Ciuchi, G. Sangiovanni and A. Toschi, *Divergent Precursors of the Mott-Hubbard Transition at the Two-Particle Level*, Phys. Rev. Lett. **110**, 246405 (2013), doi:10.1103/PhysRevLett.110.246405.

[17] O. Gunnarsson, T. Schäfer, J. P. F. LeBlanc, J. Merino, G. Sangiovanni, G. Rohringer and A. Toschi, *Parquet decomposition calculations of the electronic self-energy*, Phys. Rev. B **93**, 245102 (2016), doi:10.1103/PhysRevB.93.245102.

[18] G. Rohringer, H. Hafermann, A. Toschi, A. A. Katanin, A. E. Antipov, M. I. Katsnelson, A. I. Lichtenstein, A. N. Rubtsov and K. Held, *Diagrammatic routes to nonlocal correlations beyond dynamical mean field theory*, Rev. Mod. Phys. **90**, 025003 (2018), doi:10.1103/RevModPhys.90.025003.

[19] M. Reitner, P. Chalupa, L. Del Re, D. Springer, S. Ciuchi, G. Sangiovanni and A. Toschi, *Attractive Effect of a Strong Electronic Repulsion: The Physics of Vertex Divergences*, Phys. Rev. Lett. **125**, 196403 (2020), doi:10.1103/PhysRevLett.125.196403.

[20] S. Adler, F. Krien, P. Chalupa-Gantner, G. Sangiovanni and A. Toschi, *Non-perturbative intertwining between spin and charge correlations: A “smoking gun” single-boson-exchange result*, SciPost Phys. **16**, 054 (2024), doi:10.21468/SciPostPhys.16.2.054.

[21] C. Naya, T. Bertolini and J. Carlström, *Stability of line-node semimetals with strong Coulomb interactions and properties of the symmetry-broken state*, Phys. Rev. Res. **5**, 013069 (2023), doi:10.1103/PhysRevResearch.5.013069.

[22] F. Aryasetiawan and O. Gunnarsson, *The GW method*, Rep. Prog. Phys. **61**(3), 237 (1998), doi:10.1088/0034-4885/61/3/002.

[23] K. T. Williams, Y. Yao, J. Li, L. Chen, H. Shi, M. Motta, C. Niu, U. Ray, S. Guo, R. J. Anderson, J. Li, L. N. Tran et al., *Direct Comparison of Many-Body Methods for Realistic Electronic Hamiltonians*, Phys. Rev. X **10**, 011041 (2020), doi:10.1103/PhysRevX.10.011041.

[24] N. E. Dahlen and R. van Leeuwen, *Self-consistent solution of the Dyson equation for atoms and molecules within a conserving approximation*, J. Chem. Phys. **122**(16), 164102 (2005), doi:10.1063/1.1884965.

[25] J. J. Phillips and D. Zgid, *The description of strong correlation within self-consistent Green's function second-order perturbation theory*, J. Chem. Phys. **140**(24), 241101 (2014), doi:10.1063/1.4884951.

[26] I. S. Tupitsyn, A. S. Mishchenko, N. Nagaosa and N. Prokof'ev, *Coulomb and electron-phonon interactions in metals*, Phys. Rev. B **94**, 155145 (2016), doi:10.1103/PhysRevB.94.155145.

[27] I. S. Tupitsyn and N. V. Prokof'ev, *Stability of Dirac Liquids with Strong Coulomb Interaction*, Phys. Rev. Lett. **118**, 026403 (2017), doi:10.1103/PhysRevLett.118.026403.
[28] I. S. Tupitsyn and N. V. Prokof'ev, *Phase diagram topology of the Haldane-Hubbard-Coulomb model*, Phys. Rev. B **99**, 121113(R) (2019), doi:10.1103/PhysRevB.99.121113.

[29] M. Motta, D. M. Ceperley, G. K.-L. Chan, J. A. Gomez, E. Gull, S. Guo, C. A. Jiménez-Hoyos, T. N. Lan, J. Li, F. Ma, A. J. Millis, N. V. Prokof'ev et al., *Towards the Solution of the Many-Electron Problem in Real Materials: Equation of State of the Hydrogen Chain with State-of-the-Art Many-Body Methods*, Phys. Rev. X **7**, 031059 (2017), doi:10.1103/PhysRevX.7.031059.

[30] R. Haussmann, *Crossover from BCS superconductivity to Bose-Einstein condensation: a self-consistent theory*, Z. Phys. B **91**, 291 (1993), doi:10.1007/BF01344058.

[31] R. Haussmann, *Properties of a Fermi liquid at the superfluid transition in the crossover region between BCS superconductivity and Bose-Einstein condensation*, Phys. Rev. B **49**, 12975 (1994), doi:10.1103/PhysRevB.49.12975.

[32] R. Haussmann, W. Rantner, S. Cerrito and W. Zwerger, *Thermodynamics of the BCS-BEC crossover*, Phys. Rev. A **75**, 023610 (2007), doi:10.1103/PhysRevA.75.023610.

[33] Y. Deng, E. Kozik, N. V. Prokof'ev and B. V. Svistunov, *Emergent BCS regime of the two-dimensional fermionic Hubbard model: Ground-state phase diagram*, Europhys. Lett. **110**, 57001 (2015), doi:10.1209/0295-5075/110/57001.

[34] F. Šimkovic, Y. Deng and E. Kozik, *Superfluid ground state phase diagram of the two-dimensional Hubbard model in the emergent Bardeen-Cooper-Schrieffer regime*, Phys. Rev. B **104**, L020507 (2021), doi:10.1103/PhysRevB.104.L020507.

[35] K. Van Houcke, F. Werner, E. Kozik, N. Prokof'ev, B. Svistunov, M. J. H. Ku, A. T. Sommer, L. W. Cheuk, A. Schirotzek and M. W. Zwierlein, *Feynman diagrams versus Fermi-gas Feynman emulator*, Nature Phys. **8**, 366 (2012), doi:10.1038/NPHYS2273.

[36] R. Rossi, T. Ohgoe, K. Van Houcke and F. Werner, *Resummation of diagrammatic series with zero convergence radius for strongly correlated fermions*, Phys. Rev. Lett **121**, 130405 (2018), doi:10.1103/PhysRevLett.121.130405.

[37] R. Rossi, T. Ohgoe, E. Kozik, N. Prokof'ev, B. Svistunov, K. Van Houcke and F. Werner, *Contact and momentum distribution of the unitary Fermi gas*, Phys. Rev. Lett **121**, 130406 (2018), doi:10.1103/PhysRevLett.121.130406.

[38] A. S. Mishchenko, N. Nagaosa and N. Prokof'ev, *Diagrammatic Monte Carlo Method for Many-Polaron Problems*, Phys. Rev. Lett. **113**, 166402 (2014), doi:10.1103/PhysRevLett.113.166402.

[39] S. A. Kulagin, N. Prokof'ev, O. A. Starykh, B. Svistunov and C. N. Varney, *Bold Diagrammatic Monte Carlo Method Applied to Fermionized Frustrated Spins*, Phys. Rev. Lett. **110**, 070601 (2013), doi:10.1103/PhysRevLett.110.070601.

[40] Y. Huang, K. Chen, Y. Deng, N. Prokof'ev and B. Svistunov, *Spin-Ice State of the Quantum Heisenberg Antiferromagnet on the Pyrochlore Lattice*, Phys. Rev. Lett. **116**, 177203 (2016), doi:10.1103/PhysRevLett.116.177203.

[41] T. Wang, X. Cai, K. Chen, N. V. Prokof'ev and B. V. Svistunov, *Quantum-to-classical correspondence in two-dimensional Heisenberg models*, Phys. Rev. B **101**, 035132 (2020), doi:10.1103/PhysRevB.101.035132.
[28] I. S. Tupitsyn and N. V. Prokof'ev, *Phase diagram topology of the Haldane-Hubbard-Coulomb model*, Phys. Rev. B **99**, 121113(R) (2019), doi:10.1103/PhysRevB.99.121113.

[29] M. Motta, D. M. Ceperley, G. K.-L. Chan, J. A. Gomez, E. Gull, S. Guo, C. A. Jiménez-Hoyos, T. N. Lan, J. Li, F. Ma, A. J. Millis, N. V. Prokof'ev et al., *Towards the Solution of the Many-Electron Problem in Real Materials: Equation of State of the Hydrogen Chain with State-of-the-Art Many-Body Methods*, Phys. Rev. X **7**, 031059 (2017), doi:10.1103/PhysRevX.7.031059.

[30] R. Haussmann, *Crossover from BCS superconductivity to Bose-Einstein condensation: a self-consistent theory*, Z. Phys. B **91**, 291 (1993), doi:10.1007/BF01344058.

[31] R. Haussmann, *Properties of a Fermi liquid at the superfluid transition in the crossover region between BCS superconductivity and Bose-Einstein condensation*, Phys. Rev. B **49**, 12975 (1994), doi:10.1103/PhysRevB.49.12975.

[32] R. Haussmann, W. Rantner, S. Cerrito and W. Zwerger, *Thermodynamics of the BCS-BEC crossover*, Phys. Rev. A **75**, 023610 (2007), doi:10.1103/PhysRevA.75.023610.

[33] Y. Deng, E. Kozik, N. V. Prokof'ev and B. V. Svistunov, *Emergent BCS regime of the two-dimensional fermionic Hubbard model: Ground-state phase diagram*, Europhys. Lett. **110**, 57001 (2015), doi:10.1209/0295-5075/110/57001.

[34] F. Šimkovic, Y. Deng and E. Kozik, *Superfluid ground state phase diagram of the two-dimensional Hubbard model in the emergent Bardeen-Cooper-Schrieffer regime*, Phys. Rev. B **104**, L020507 (2021), doi:10.1103/PhysRevB.104.L020507.

[35] K. Van Houcke, F. Werner, E. Kozik, N. Prokof'ev, B. Svistunov, M. J. H. Ku, A. T. Sommer, L. W. Cheuk, A. Schirotzek and M. W. Zwierlein, *Feynman diagrams versus Fermi-gas Feynman emulator*, Nature Phys. **8**, 366 (2012), doi:10.1038/NPHYS2273.

[36] R. Rossi, T. Ohgoe, K. Van Houcke and F. Werner, *Resummation of diagrammatic series with zero convergence radius for strongly correlated fermions*, Phys. Rev. Lett **121**, 130405 (2018), doi:10.1103/PhysRevLett.121.130405.

[37] R. Rossi, T. Ohgoe, E. Kozik, N. Prokof'ev, B. Svistunov, K. Van Houcke and F. Werner, *Contact and momentum distribution of the unitary Fermi gas*, Phys. Rev. Lett **121**, 130406 (2018), doi:10.1103/PhysRevLett.121.130406.

[38] A. S. Mishchenko, N. Nagaosa and N. Prokof'ev, *Diagrammatic Monte Carlo Method for Many-Polaron Problems*, Phys. Rev. Lett. **113**, 166402 (2014), doi:10.1103/PhysRevLett.113.166402.

[39] S. A. Kulagin, N. Prokof'ev, O. A. Starykh, B. Svistunov and C. N. Varney, *Bold Diagrammatic Monte Carlo Method Applied to Fermionized Frustrated Spins*, Phys. Rev. Lett. **110**, 070601 (2013), doi:10.1103/PhysRevLett.110.070601.

[40] Y. Huang, K. Chen, Y. Deng, N. Prokof'ev and B. Svistunov, *Spin-Ice State of the Quantum Heisenberg Antiferromagnet on the Pyrochlore Lattice*, Phys. Rev. Lett. **116**, 177203 (2016), doi:10.1103/PhysRevLett.116.177203.

[41] T. Wang, X. Cai, K. Chen, N. V. Prokof'ev and B. V. Svistunov, *Quantum-to-classical correspondence in two-dimensional Heisenberg models*, Phys. Rev. B **101**, 035132 (2020), doi:10.1103/PhysRevB.101.035132.
[42] T. Ayral and O. Parcollet, *Mott physics and spin fluctuations: A functional viewpoint*, Phys. Rev. B **93**, 235124 (2016), doi:10.1103/PhysRevB.93.235124.

[43] J. W. Negele and H. Orland, *Quantum Many-particle Systems*, Addison-Wesley, doi:10.1201/9780429497926 (1988).

[44] N. V. Prokof'ev and B. V. Svistunov, *Bold diagrammatic Monte Carlo: A generic sign-problem tolerant technique for polaron models and possibly interacting many-body problems*, Phys. Rev. B **77**, 125101 (2008), doi:10.1103/PhysRevB.77.125101.

[45] K. Van Houcke, F. Werner, T. Ohgoe, N. Prokof'ev and B. Svistunov, *Diagrammatic Monte Carlo algorithm for the resonant Fermi gas*, Phys. Rev. B **99**, 035140 (2019), doi:10.1103/PhysRevB.99.035140.

[46] N. Prokof'ev and B. Svistunov, *Bold Diagrammatic Monte Carlo Technique: When the Sign Problem Is Welcome*, Phys. Rev. Lett. **99**, 250201 (2007), doi:10.1103/PhysRevLett.99.250201.

[47] R. Rossi, *Determinant Diagrammatic Monte Carlo in the Thermodynamic Limit*, Phys. Rev. Lett. **119**, 045701 (2017), doi:10.1103/PhysRevLett.119.045701.

[48] F. Šimkovic and E. Kozik, *Determinant Monte Carlo for irreducible Feynman diagrams in the strongly correlated regime*, Phys. Rev. B **100**, 121102(R) (2019), doi:10.1103/PhysRevB.100.121102.

[49] A. Moutenet, W. Wu and M. Ferrero, *Determinant Monte Carlo algorithms for dynamical quantities in fermionic systems*, Phys. Rev. B **97**, 085117 (2018), doi:10.1103/PhysRevB.97.085117.

[50] F. Šimkovic, J. P. F. LeBlanc, A. J. Kim, Y. Deng, N. V. Prokof'ev, B. V. Svistunov and E. Kozik, *Extended Crossover from a Fermi Liquid to a Quasiantiferromagnet in the Half-Filled 2D Hubbard Model*, Phys. Rev. Lett. **124**, 017003 (2020), doi:10.1103/PhysRevLett.124.017003.

[51] A. J. Kim, F. Simkovic and E. Kozik, *Spin and Charge Correlations across the Metal-to-Insulator Crossover in the Half-Filled 2D Hubbard Model*, Phys. Rev. Lett. **124**, 117602 (2020), doi:10.1103/PhysRevLett.124.117602.

[52] J. Carlström, *Diagrammatic Monte Carlo procedure for the spin-charge transformed Hubbard model*, Phys. Rev. B **97**, 075119 (2018), doi:10.1103/PhysRevB.97.075119.

[53] K. Chen and K. Haule, *A combined variational and diagrammatic quantum Monte Carlo approach to the many-electron problem*, Nature Comm. **10**, 3725 (2019), doi:10.1038/s41467-019-11708-6.

[54] K. Haule and K. Chen, *Single-particle excitations in the uniform electron gas by diagrammatic Monte Carlo*, Sci. Rep. **12**(1), 2294 (2022), doi:10.1038/s41598-022-06188-6.

[55] J. Li, M. Wallerberger and E. Gull, *Diagrammatic Monte Carlo method for impurity models with general interactions and hybridizations*, Phys. Rev. Research **2**, 033211 (2020), doi:10.1103/PhysRevResearch.2.033211.

[56] F. Šimkovic, R. Rossi and M. Ferrero, *Efficient one-loop-renormalized vertex expansions with connected determinant diagrammatic Monte Carlo*, Phys. Rev. B **102**, 195122 (2020), doi:10.1103/PhysRevB.102.195122.
[42] T. Ayral and O. Parcollet, *Mott physics and spin fluctuations: A functional viewpoint*, Phys. Rev. B **93**, 235124 (2016), doi:10.1103/PhysRevB.93.235124.

[43] J. W. Negele and H. Orland, *Quantum Many-particle Systems*, Addison-Wesley, doi:10.1201/9780429497926 (1988).

[44] N. V. Prokof'ev and B. V. Svistunov, *Bold diagrammatic Monte Carlo: A generic sign-problem tolerant technique for polaron models and possibly interacting many-body problems*, Phys. Rev. B **77**, 125101 (2008), doi:10.1103/PhysRevB.77.125101.

[45] K. Van Houcke, F. Werner, T. Ohgoe, N. Prokof'ev and B. Svistunov, *Diagrammatic Monte Carlo algorithm for the resonant Fermi gas*, Phys. Rev. B **99**, 035140 (2019), doi:10.1103/PhysRevB.99.035140.

[46] N. Prokof'ev and B. Svistunov, *Bold Diagrammatic Monte Carlo Technique: When the Sign Problem Is Welcome*, Phys. Rev. Lett. **99**, 250201 (2007), doi:10.1103/PhysRevLett.99.250201.

[47] R. Rossi, *Determinant Diagrammatic Monte Carlo in the Thermodynamic Limit*, Phys. Rev. Lett. **119**, 045701 (2017), doi:10.1103/PhysRevLett.119.045701.

[48] F. Šimkovic and E. Kozik, *Determinant Monte Carlo for irreducible Feynman diagrams in the strongly correlated regime*, Phys. Rev. B **100**, 121102(R) (2019), doi:10.1103/PhysRevB.100.121102.

[49] A. Moutenet, W. Wu and M. Ferrero, *Determinant Monte Carlo algorithms for dynamical quantities in fermionic systems*, Phys. Rev. B **97**, 085117 (2018), doi:10.1103/PhysRevB.97.085117.

[50] F. Šimkovic, J. P. F. LeBlanc, A. J. Kim, Y. Deng, N. V. Prokof'ev, B. V. Svistunov and E. Kozik, *Extended Crossover from a Fermi Liquid to a Quasiantiferromagnet in the Half-Filled 2D Hubbard Model*, Phys. Rev. Lett. **124**, 017003 (2020), doi:10.1103/PhysRevLett.124.017003.

[51] A. J. Kim, F. Simkovic and E. Kozik, *Spin and Charge Correlations across the Metal-to-Insulator Crossover in the Half-Filled 2D Hubbard Model*, Phys. Rev. Lett. **124**, 117602 (2020), doi:10.1103/PhysRevLett.124.117602.

[52] J. Carlström, *Diagrammatic Monte Carlo procedure for the spin-charge transformed Hubbard model*, Phys. Rev. B **97**, 075119 (2018), doi:10.1103/PhysRevB.97.075119.

[53] K. Chen and K. Haule, *A combined variational and diagrammatic quantum Monte Carlo approach to the many-electron problem*, Nature Comm. **10**, 3725 (2019), doi:10.1038/s41467-019-11708-6.

[54] K. Haule and K. Chen, *Single-particle excitations in the uniform electron gas by diagrammatic Monte Carlo*, Sci. Rep. **12**(1), 2294 (2022), doi:10.1038/s41598-022-06188-6.

[55] J. Li, M. Wallerberger and E. Gull, *Diagrammatic Monte Carlo method for impurity models with general interactions and hybridizations*, Phys. Rev. Research **2**, 033211 (2020), doi:10.1103/PhysRevResearch.2.033211.

[56] F. Šimkovic, R. Rossi and M. Ferrero, *Efficient one-loop-renormalized vertex expansions with connected determinant diagrammatic Monte Carlo*, Phys. Rev. B **102**, 195122 (2020), doi:10.1103/PhysRevB.102.195122.
[57] R. Rossi, F. Šimkovic and M. Ferrero, Renormalized perturbation theory at large expansion orders, Europhys. Lett. **132**, 11001 (2020), doi:10.1209/0295-5075/132/11001.

[58] F. Šimkovic, R. Rossi and M. Ferrero, Two-dimensional Hubbard model at finite temperature: Weak, strong, and long correlation regimes, Phys. Rev. Res. **4**, 043201 (2022), doi:10.1103/PhysRevResearch.4.043201.

[59] C. Lenihan, A. J. Kim, F. Šimkovic and E. Kozik, Evaluating Second-Order Phase Transitions with Diagrammatic Monte Carlo: Néel Transition in the Doped Three-Dimensional Hubbard Model, Phys. Rev. Lett. **129**, 107202 (2022), doi:10.1103/PhysRevLett.129.107202.

[60] G. Spada, R. Rossi, F. Simkovic, R. Garioud, M. Ferrero, K. Van Houcke and F. Werner, High-order expansion around BCS theory, arXiv:2103.12038.

[61] F. Simkovic, R. Rossi, A. Georges and M. Ferrero, arXiv:2209.09237.

[62] A. Wietek, R. Rossi, F. Šimkovic, M. Klett, P. Hansmann, M. Ferrero, E. M. Stoudenmire, T. Schäfer and A. Georges, Mott Insulating States with Competing Orders in the Triangular Lattice Hubbard Model, Phys. Rev. X **11**, 041013 (2021), doi:10.1103/PhysRevX.11.041013.

[63] P.-C. Hou, B.-Z. Wang, K. Haule, Y. Deng and K. Chen, Exchange-correlation effect in the charge response of a warm dense electron gas, Phys. Rev. B **106**, L081126 (2022), doi:10.1103/PhysRevB.106.L081126.

[64] I. S. Tupitsyn, A. M. Tsvelik, R. M. Konik and N. V. Prokof'ev, Real-Frequency Response Functions at Finite Temperature, Phys. Rev. Lett. **127**, 026403 (2021), doi:10.1103/PhysRevLett.127.026403.

[65] J. P. F. LeBlanc, K. Chen, K. Haule, N. V. Prokof'ev and I. S. Tupitsyn, Dynamic Response of an Electron Gas: Towards the Exact Exchange-Correlation Kernel, Phys. Rev. Lett. **129**, 246401 (2022), doi:10.1103/PhysRevLett.129.246401.

[66] I. S. Tupitsyn and N. V. Prokof'ev, Landau Damping in an Electron Gas, arXiv:2311.05611.

[67] E. Kozik, Combinatorial summation of Feynman diagrams: Equation of state of the 2D SU(N) Hubbard model, arXiv:2309.13774.

[68] R. Rossi, F. Simkovic, M. Ferrero, A. Georges, A. M. Tsvelik, N. V. Prokof'ev and I. S. Tupitsyn, Interaction-enhanced nesting in Spin-Fermion and Fermi-Hubbard models, arXiv:2402.13238.

[69] J. Carlström, Strong-coupling diagrammatic Monte Carlo technique for correlated fermions and frustrated spins, Phys. Rev. B **103**, 195147 (2021), doi:10.1103/PhysRevB.103.195147.

[70] J. Carlström, In situ controllable magnetic phases in doped twisted bilayer transition metal dichalcogenides, Phys. Rev. Res. **4**, 043126 (2022), doi:10.1103/PhysRevResearch.4.043126.

[71] J. Carlström, Spectral topology and its relation to Fermi arcs in strongly correlated systems, Phys. Rev. Res. **5**, 033160 (2023), doi:10.1103/PhysRevResearch.5.033160.

[72] A. J. Kim, P. Werner and E. Kozik, Strange Metal Solution in the Diagrammatic Theory for the 2d Hubbard Model, arXiv:2012.06159.

[73] R. Rossi et al., to be published.
[57] R. Rossi, F. Šimkovic and M. Ferrero, Renormalized perturbation theory at large expansion orders, Europhys. Lett. 132, 11001 (2020), doi:10.1209/0295-5075/132/11001.

[58] F. Šimkovic, R. Rossi and M. Ferrero, Two-dimensional Hubbard model at finite temperature: Weak, strong, and long correlation regimes, Phys. Rev. Res. 4, 043201 (2022), doi:10.1103/PhysRevResearch.4.043201.

[59] C. Lenihan, A. J. Kim, F. Šimkovic and E. Kozik, Evaluating Second-Order Phase Transitions with Diagrammatic Monte Carlo: Néel Transition in the Doped Three-Dimensional Hubbard Model, Phys. Rev. Lett. 129, 107202 (2022), doi:10.1103/PhysRevLett.129.107202.

[60] G. Spada, R. Rossi, F. Simkovic, R. Garioud, M. Ferrero, K. Van Houcke and F. Werner, High-order expansion around BCS theory, arXiv:2103.12038.

[61] F. Simkovic, R. Rossi, A. Georges and M. Ferrero, arXiv:2209.09237.

[62] A. Wietek, R. Rossi, F. Šimkovic, M. Klett, P. Hansmann, M. Ferrero, E. M. Stoudenmire, T. Schäfer and A. Georges, Mott Insulating States with Competing Orders in the Triangular Lattice Hubbard Model, Phys. Rev. X 11, 041013 (2021), doi:10.1103/PhysRevX.11.041013.

[63] P.-C. Hou, B.-Z. Wang, K. Haule, Y. Deng and K. Chen, Exchange-correlation effect in the charge response of a warm dense electron gas, Phys. Rev. B 106, L081126 (2022), doi:10.1103/PhysRevB.106.L081126.

[64] I. S. Tupitsyn, A. M. Tsvelik, R. M. Konik and N. V. Prokof'ev, Real-Frequency Response Functions at Finite Temperature, Phys. Rev. Lett. 127, 026403 (2021), doi:10.1103/PhysRevLett.127.026403.

[65] J. P. F. LeBlanc, K. Chen, K. Haule, N. V. Prokof'ev and I. S. Tupitsyn, Dynamic Response of an Electron Gas: Towards the Exact Exchange-Correlation Kernel, Phys. Rev. Lett. 129, 246401 (2022), doi:10.1103/PhysRevLett.129.246401.

[66] I. S. Tupitsyn and N. V. Prokof'ev, Landau Damping in an Electron Gas, arXiv:2311.05611.

[67] E. Kozik, Combinatorial summation of Feynman diagrams: Equation of state of the 2D SU(N) Hubbard model, arXiv:2309.13774.

[68] R. Rossi, F. Simkovic, M. Ferrero, A. Georges, A. M. Tsvelik, N. V. Prokof'ev and I. S. Tupitsyn, Interaction-enhanced nesting in Spin-Fermion and Fermi-Hubbard models, arXiv:2402.13238.

[69] J. Carlström, Strong-coupling diagrammatic Monte Carlo technique for correlated fermions and frustrated spins, Phys. Rev. B 103, 195147 (2021), doi:10.1103/PhysRevB.103.195147.

[70] J. Carlström, In situ controllable magnetic phases in doped twisted bilayer transition metal dichalcogenides, Phys. Rev. Res. 4, 043126 (2022), doi:10.1103/PhysRevResearch.4.043126.

[71] J. Carlström, Spectral topology and its relation to Fermi arcs in strongly correlated systems, Phys. Rev. Res. 5, 033160 (2023), doi:10.1103/PhysRevResearch.5.033160.

[72] A. J. Kim, P. Werner and E. Kozik, Strange Metal Solution in the Diagrammatic Theory for the 2d Hubbard Model, arXiv:2012.06159.

[73] R. Rossi et al., to be published.