---
abstract: 'The order-by-order renormalization of the self-consistent mean-field potential in many-body perturbation theory for normal Fermi systems is investigated in detail. Building on previous work mainly by Balian and de Dominicis, as a key result we derive a thermodynamic perturbation series that manifests the consistency of the adiabatic zero-temperature formalism with perturbative statistical mechanics—for both isotropic and anisotropic systems—and satisfies at each order and for all temperatures the thermodynamic relations associated with Fermi-liquid theory. These properties are proved to all orders.'
author:
- Corbinian Wellenhofer
bibliography:
- 'refs.bib'
title: 'Zero-Temperature Limit and Statistical Quasiparticles in Many-Body Perturbation Theory'
---

Introduction
============

Many-body perturbation theory (MBPT) represents the elementary framework for calculations aimed at the properties of nonrelativistic many-fermion systems at zero and finite temperature. In general, for Fermi systems the correct ground-state is not a normal state but involves Cooper pairs [@PhysRev.150.202; @PhysRevLett.15.524; @doi:10.1142/S0217979292001249; @RevModPhys.66.129; @Salmhofer:1999uq]. However, pairing effects can often be neglected for approximative calculations of thermodynamic properties close to zero temperature. For such calculations there are two formalisms: first, there is grand-canonical perturbation theory, and second, the zero-temperature formalism based on the adiabatic continuation of the ground state [@PhysRev.84.350; @Goldstone:1957zz; @Nozbook; @Runge; @negele; @Fetter; @Abrikosov]. In their time-dependent (i.e., in frequency space) formulations, these two formalisms give matching results if all quantities are derived from the exact Green’s functions, i.e., from the self-consistently renormalized propagators [@Luttinger:1960ua; @Fetter; @Abrikosov; @Feldman1999]. The renormalization of MBPT in frequency space can be generalized to vertex functions [@dommar2; @dommar1; @Hausmann; @Rossi:2015lda; @PhysRevB.93.161102; @DICKHOFF2004377; @VanHoucke:2011ux], and is essential to obtain a fully consistent framework for calculating transport properties [@Baym:1961zz; @Baym:1962sx; @stefanuc].

Nevertheless, the use of bare propagators has the benefit that in that case the time integrals can be performed analytically. With bare propagators, MBPT in its most basic form corresponds to a perturbative expansion in terms of the interaction Hamiltonian $V$ about the noninteracting system with Hamiltonian $H_0$, where $H=H_0+V$ is the full Hamiltonian. First-order self-energy effects can be included to all orders in bare MBPT by expanding instead about a reference Hamiltonian $H_\text{ref}=H_0+U_1$, where $U_1$ includes the first-order contribution to the (frequency-space) self-energy $\Sigma_{1,{{\textbf}{k}}}$ as a self-consistent single-particle potential (mean field). The renormalization of $H_\text{ref}$ in terms of $U_1$ has the effect that all two-particle reducible diagrams with first-order pieces (single-vertex loops) are canceled. At second order the self-energy becomes frequency dependent and complex, so the equivalence between the propagator renormalization in frequency space and the renormalization of the mean-field part of $H_\text{ref}$ in bare MBPT is restricted to the Hartree-Fock level.

Zero-temperature MBPT calculations with bare propagators and a Hartree-Fock reference Hamiltonian ${H_\text{ref}=H_0+U_1}$ are common in quantum chemistry and nuclear physics. With a Hartree-Fock reference Hamiltonian (or, with ${H_\text{ref}=H_0}$), however, the adiabatic zero-temperature formalism is inconsistent with the zero-temperature limit (${T\rightarrow 0}$) of grand-canonical MBPT. The (main) fault however lies not with zero-temperature MBPT, but with the grand-canonical perturbation series: in the bare grand-canonical formalism (with $H_\text{ref}\in\{H_0,H_0+U_1\}$) there is a mismatch in the Fermi-Dirac distribution functions caused by using the *reference* spectrum $\varepsilon_{{\textbf}{k}}$ together with the *true* chemical potential $\mu$, and in general this leads to deficient results [@Fritsch:2002hp; @PhysRevC.89.064009; @Wellenhofer:2017qla]. The adiabatic formalism on the other hand uses the reference chemical potential, i.e., the reference Fermi energy ${\varepsilon}_{\text{F}}$. Related to this is the presence of additional contributions from two-particle reducible diagrams, the so-called anomalous contributions, in the grand-canonical formalism.

This issue is usually dealt with by modifying the grand-canonical perturbation series for the free energy in terms of an expansion about the chemical potential $\mu_\text{ref}\xrightarrow{T\rightarrow 0}{\varepsilon}_{\text{F}}$ of the reference system [@Kohn:1960zz; @brout2] (see also Sec. \[sec42\]). This expansion introduces additional anomalous contributions, and for isotropic systems these can be seen to cancel the old ones for ${T\rightarrow 0}$ [@Luttinger:1960ua]. Thus, the modified perturbation series for the free energy $F(T,\mu_\text{ref})$ reproduces the adiabatic series in the isotropic case. For anisotropic systems, however, the anomalous contributions persist at ${T=0}$ (for ${H_\text{ref}=H_0+U_1}$, at fourth order and beyond). Negele and Orland [@negele] interpret this feature as follows: there is nothing fundamentally wrong with the bare zero-temperature formalism, but for anisotropic systems the adiabatic continuation must be based on a better reference Hamiltonian $H_\text{ref}$. Since the convergence rate[^1] of MBPT depends on the choice of $H_\text{ref}$, this issue is relevant also for finite-temperature calculations, and for isotropic systems.

Recently, Holt and Kaiser [@PhysRevC.95.034326] have shown that including the real part of the bare second-order contribution to the (on-shell) self-energy, $\text{Re}\,[\Sigma_{2,{{\textbf}{k}}}(\varepsilon_{{\textbf}{k}})]$, as the second-order contribution to the self-consistent mean field has a significant effect in perturbative nuclear matter calculations with modern two- and three-nucleon potentials (see, e.g., Refs. [@RevModPhys.81.1773; @MACHLEIDT20111; @Bogner:2009bt]). However, a formal clarification for the renormalization of $H_\text{ref}$ in terms of $\text{Re}\,[\Sigma_{2,{{\textbf}{k}}}(\varepsilon_{{\textbf}{k}})]$ was not included in Ref. [@PhysRevC.95.034326]. In particular, from the discussion of Ref. [@PhysRevC.95.034326] it is not clear whether the use of this second-order mean field should be considered an improvement or not, compared to calculations with a Hartree-Fock mean field.[^2]

A general scheme where the reference Hamiltonian is renormalized at each order in grand-canonical MBPT was introduced by Balian, Bloch, and de Dominicis [@Balian1961529] (see also Refs. [@Balian1961529b; @Balianph; @BlochImperFerm; @boer; @1964mbpdedom]). This scheme however leads to a mean field whose functional form is given by $U[n_{\textbf}{k},T]$, where $n_{\textbf}{k}(T,\mu)$ is the Fermi-Dirac distribution and the explicit temperature dependence involves factors $\operatorname{\operatorname{e}}^{\pm({\varepsilon}_{\textbf}{k}-\mu)/T}$. Because of the $\operatorname{\operatorname{e}}^{\pm({\varepsilon}_{\textbf}{k}-\mu)/T}$ factors, the resulting perturbation series is well-behaved only at sufficiently large temperatures, and its ${T\rightarrow 0}$ limit does not exist.[^3]

A different renormalization scheme was outlined by Balian and de Dominicis (BdD) in Refs. [@statquasi3; @statquasi1] (see also Refs. [@1964mbpdedom; @boer]). At second order, this scheme leads to the mean field employed by Holt and Kaiser [@PhysRevC.95.034326]. The outline given in Refs. [@statquasi3; @statquasi1] indicates the following results:

1.  The functional form of the mean field is to all orders given by $U[n_{\textbf}{k}]$, i.e., there is no explicit temperature dependence (apart from the one given by the Fermi-Dirac distributions), so the ${T\rightarrow 0}$ limit exists.

2.  The zero-temperature limit of the renormalized grand-canonical perturbation series for the free energy $F(T,\mu)$ reproduces the (correspondingly renormalized) adiabatic series for the ground-state energy $E^{(0)}({\varepsilon}_{\text{F}})$ to all orders; i.e., the reference spectrum ${\varepsilon}_{\textbf}{k}$ has been adjusted to the true chemical potential $\mu$, with ${\varepsilon}_{\text{F}}=\mu$ at ${T=0}$.

3.  One obtains at each perturbative order and for all temperatures the thermodynamic relations associated with Fermi-liquid theory [@Landau]. This result corresponds to the notion statistical quasiparticles [@PhysRevA.7.304; @pethick2; @Fermiliq].

The most intricate part in establishing these results is as follows. For ${T\neq 0}$, there are no energy denominator poles in the (proper) expressions for the perturbative contributions to the grand-canonical potential. The BdD renormalization scheme however introduces such poles, and therefore a regularization procedure is required to apply the scheme. So far, this issue has been studied in more detail only for the case of impurity systems [@BALIAN1971229; @LUTTINGER19731; @keitermorandi].

Motivated by this situation, in the present paper we revisit the order-by-order renormalization of the reference Hamiltonian in bare MBPT.[^4] First, in Sec. \[sec2\] we give a short review of grand-canonical perturbation theory with bare propagators and introduce the various order-by-order renormalizations of the reference Hamiltonian. We also discuss how dynamical quasiparticles arise in (frequency-space) MBPT, and show that their energies are distinguished from the ones of the statistical quasiparticles associated with result (2). In Sec. \[sec3\] we discuss the regularization procedure for the BdD renormalization scheme, and analyze the resulting expressions for the second- and third-order contributions to the grand-canonical potential and the BdD mean field. In Sec. \[sec4\] we prove to all orders that the BdD renormalized perturbation series satisfies the Fermi-liquid relations (2) and, as a consequence, manifests the consistency of the adiabatic zero-temperature formalism (1). The paper is concluded in Sec. \[summary\]. In Appendix \[app1\], we derive explicitly the renormalized contribution from two-particle reducible diagrams at fourth order. In Appendix \[app2\], we discuss in more detail the various forms of the self-energy, derive various expressions for the mean occupation numbers, and examine the functional relations between the grand-canonical potential and the (various forms of the) self-energy in bare MBPT.

Grand-Canonical Perturbation Theory {#sec2}
===================================

Setup {#sec21}
-----

We consider a homogeneous but not necessarily isotropic system of nonrelativistic fermions in thermodynamic equilibrium. The Hamiltonian is given by ${H= H_0+V}$, where $V$ is a two-body operator representing pair interactions. Multi-fermion interactions do not raise any new formal or conceptual issues, and are therefore neglected. For notational simplicity and without loss of generality we assume a single species of spinless fermions. If there is no external potential, then $H_0$ is the kinetic energy operator. We now introduce an additional one-body operator $U$, and write $$\begin{aligned}
\label{Eq0}
H= \underbrace{(H_0+U)}_{H_\text{ref}}+(V- U).\end{aligned}$$ The operator $U$ represents a mean field, i.e., an effective one-body potential which allows to define a solvable reference system that includes the effects of pair interactions in the system to a certain degree. For a homogeneous system the mean field should preserve translational invariance, so the eigenstates $\ket{\psi_{{\textbf}{k}}}$ of the momentum operator are eigenstates of $H_\text{ref}=H_0+U$, i.e., $$\begin{aligned}
H_{\text{ref}} \ket{\psi_{{\textbf}{k}}} = {\varepsilon}_{{\textbf}{k}}\ket{\psi_{{\textbf}{k}}}.\end{aligned}$$ Because the mean field is supposed to include interaction effects self-consistently, the single-particle energies ${\varepsilon}_{{\textbf}{k}}$ are determined by the self-consistent equation $$\begin{aligned}
{\varepsilon}_{{\textbf}{k}} = {\varepsilon}_{0,{\textbf}{k}} + U_{{\textbf}{k}}[{\varepsilon}_{{\textbf}{k}}],\end{aligned}$$ where ${\varepsilon}_{0,{\textbf}{k}}=\braket{\psi_{{\textbf}{k}}|H_0|\psi_{{\textbf}{k}}}$ and $U_{{\textbf}{k}}[{\varepsilon}_{{\textbf}{k}}]=\braket{\psi_{{\textbf}{k}}|U|\psi_{{\textbf}{k}}}$.[^5] The occupation number representation of the reference Hamiltonian $H_{\text{ref}}=H_0+ U$ is then given by $$\begin{aligned}
\mathcal{H}_{\text{ref}} &= \sum_{{\textbf}{k}} \braket{\psi_{{\textbf}{k}}|H_{\text{ref}}|\psi_{{\textbf}{k}} }  a^\dagger_{{\textbf}{k}} a_{{\textbf}{k}},\end{aligned}$$ where $a^\dagger_{\textbf}{k}$ and $a_{\textbf}{k}$ are creation and annihiliation operators with respect to momentum eigenstates. If not indicated explicitly otherwise, we assume the thermodynamic limit where $\sum_{\textbf}{k}\rightarrow\int \! d^3 k/(2\pi)^3$.[^6] The occupation number representation of the perturbation Hamiltonian ${H-H_\text {ref}=V-U}$ is given by $$\begin{aligned}
\label{VHamilton}
\mathcal{V}&= \frac{1}{2!}\sum_{{\textbf}{k}_1,{\textbf}{k}_2,{\textbf}{k}_3,{\textbf}{k}_4} 
\braket{\psi_{{\textbf}{k}_1}\psi_{{\textbf}{k}_2}| V|\psi_{{\textbf}{k}_3}\psi_{{\textbf}{k}_4} } a^\dagger_{{\textbf}{k}_1} a^\dagger_{{\textbf}{k}_2} a_{{\textbf}{k}_4} a_{{\textbf}{k}_3}
{\nonumber \\}& \quad
-\sum_{{\textbf}{k}}
\braket{\psi_{{\textbf}{k}}|  U|\psi_{{\textbf}{k}} } a^\dagger_{{\textbf}{k}}  a_{{\textbf}{k}},\end{aligned}$$ where momentum conservation is implied, i.e., ${\textbf}{k}_1+{\textbf}{k}_2$$={\textbf}{k}_3+{\textbf}{k}_4$. We assume that the potential $V$ is sufficiently regular(ized) such that no ultraviolet [@Hammer:2000xg; @Wellenhofer:2018dwh] or infrared [@PhysRev.106.364] divergences appear in perturbation theory.[^7] Further, we require that $V$ has a form (e.g., finite-rangedinteractions) for which the thermodynamic limit exists; see, e.g., Refs. [@haagbook; @ruelle; @lieb].

Perturbation series and diagrammatic analysis {#sec22}
---------------------------------------------

### Grand-canonical perturbation series

For truncation order $N$, the perturbation series for the grand-canonical potential $\Omega(T,\mu)$ is given by $$\begin{aligned}
 \label{MBPTgc0}
\Omega(T,\mu) = \Omega_\text{ref}(T,\mu) 
+\Omega_{U}(T,\mu) 
+  \sum_{n=1}^{N} \Omega_n(T,\mu),\end{aligned}$$ where $$\begin{aligned}
\Omega_\text{ref}(T,\mu) &= T\sum_{{\textbf}{k}}\ln(\bar n_{\textbf}{k}),
{\nonumber \\}\Omega_{U}(T,\mu)&= -
\sum_{{\textbf}{k}} U_{{\textbf}{k}} n_{\textbf}{k}.\end{aligned}$$ Here, $\bar n_{\textbf}{k}=1-n_{\textbf}{k}$, with $n_{\textbf}{k}=[1+\operatorname{\operatorname{e}}^{\beta({\varepsilon}_{\textbf}{k}-\mu}]^{-1}$ the Fermi-Dirac distribution function, and $\beta=1/T$. From the grand-canonical version of Wick’s theorem one obtains the following formula [@Fetter; @BDDnuclphys7] for $\Omega_{n}(T,\mu)$: $$\begin{aligned}
 \label{OmegaT}  
\Omega_n^{\text{direct}[P]}=
-\frac{1}{\beta}\frac{(-1)^{n}}{n!} 
\int\limits_{0}^\beta \!
d \tau_n \cdots d \tau_1 \;
\Braket{  \mathcal{T}\big[
\mathcal{V}(\tau_n) \cdots
\mathcal{V}(\tau_1) \big]
}_{L},\end{aligned}$$ where $\mathcal{T}$ is the time-ordering operator and $\mathcal{V}(\tau)=\operatorname{\operatorname{e}}^{\mathcal{H}_\text{ref} \tau}\mathcal{V}\operatorname{\operatorname{e}}^{-\mathcal{H}_\text{ref} \tau}$ is the interaction picture representation (in imaginary time) of the perturbation operator $\mathcal{V}$ given by Eq. .

### Classification of diagrams

The various ways the Wick contractions in the unperturbed ensemble average $\braket{\ldots}$ can be performed can be represented by Hugenholtz diagrams, i.e., diagrams composed of $V$ and $-U$ vertices,[^8] and directed lines attached to vertices at both ends. Left-pointing lines are called holes and correspond to factors $n_{\textbf}{k}$, right-pointing lines are called particles and have factors $\bar n_{\textbf}{k}$. In the case of two-particle reducible diagrams, momentum conservation implies that there are two or more lines with identical three-momenta. We refer to these lines as articulation lines. The diagrammatic parts connected via articulation lines are referred to as pieces. Two-particle irreducible diagrams have only $V$ vertices. Two-particle reducible diagrams where at least one set of lines with identical three-momenta includes both holes and particles are called anomalous, with the indicative lines referred to as anomalous articulation lines. All other (two-particle reducible or irreducible) diagrams are called normal. The parts of anomalous diagrams connected via anomalous articulation lines are called normal pieces.[^9] In general, normal two-particle reducible diagrams transform into anomalous diagrams under vertex permutations, see Figs. \[fig3red\], \[figx\] and \[fig4\].

In Eq. , the subscript $L$ means that only linked diagrams are taken into account. By virtue of the time integration and the time-ordering operator, in Eq.  there is no distinction between diagrams connected via vertex permutations; in particular, there is no distinction between normal and anomalous two-particle reducible diagrams. The distinction between the different diagrams in the permutation invariant sets of diagrams is however relevant for the time-independent formulas discussed below.

### Time-independent formulas

From Eq. , Bloch and de Dominicis [@BDDnuclphys7] (see also Refs. [@1964mbpdedom0; @boer; @keitermorandi]) have derived several time-independent formulas for $\Omega_{n}(T,\mu)$. One of them, here referred to as the , is given by $$\begin{aligned}
 \label{direct}
\Omega_{n}^\text{direct}&=\frac{1}{\beta} \frac{(-1)^{n}}{2\pi {\text{i}}} 
\oint_{C}  dz \frac{\operatorname{\operatorname{e}}^{-\beta z}}{z^2}
\Braket{ 
\mathcal{V} \frac{1}{{D}_n-z} \cdots
\mathcal{V} \frac{1}{{D}_1-z}
\mathcal{V}
}_{\!L},\end{aligned}$$ where the contour $C$ encloses all the poles ${z=0,{D}_1,\ldots,{D}_n}$, with ${D}_{\nu\in\{1,\ldots,n\}}$ the energy denominators for the respective diagrams. Furthermore, in Eq. , it is implied that the contributions from all poles are summed before the momentum integration, i.e., the $z$ integral is performed inside the momentum integrals. This has the consequence that the integrands of the momentum integrals have no poles (for ${T\neq 0}$, see below) from vanishing energy denominators. The expressions obtained from the direct formula deviate from the ones obtained from the time-dependent formula Eq. , but—as evident from the derivation of direct formula [@BDDnuclphys7]—the sum of the direct expressions obtained for a set of diagrams that is closed under vertex permutations is equivalent (but not identical) to the expression obtained from Eq. .

From the cyclic property of the trace, another time-independent formula can be derived [@BDDnuclphys7], here referred to as the cyclic formula, i.e., $$\begin{aligned}
 \label{cyclic}
\Omega_{n}^\text{cyclic}&=\frac{1}{n} \frac{(-1)^{n+1}}{2\pi {\text{i}}} 
\oint_{C}  dz \frac{\operatorname{\operatorname{e}}^{-\beta z}}{z}
\Braket{ 
\mathcal{V} \frac{1}{{D}_n-z} \cdots
\mathcal{V} \frac{1}{{D}_1-z}
\mathcal{V}
}_{\!L},\end{aligned}$$ where again it is implied that the $z$ integral is performed inside the momentum integrals; again, this has the consequence that the integrands have no poles (for ${T\neq 0}$). The direct and the cyclic formula give equivalent (but not identical) expressions only for the sums of diagrams connected via cyclic vertex permutations, and the cyclic expressions for the individual diagrams in these cyclic groups are equivalent.

Finally, from the analysis of the contributions from the different poles in Eq.  one can *formally* write down a reduced form of the cyclic formula [@BDDnuclphys7], here referred to as the reduced formula, i.e., $$\begin{aligned}
 \label{reduced}
\Omega_{n}^\text{reduced}&=  \frac{(-1)^{n+1}}{\mathcal{O}}
\underset{z=0}{\text{Res}}  \frac{\operatorname{\operatorname{e}}^{-\beta z}}{z}
\Braket{ 
\mathcal{V} \frac{1}{{D}_n-z} \cdots
\mathcal{V} \frac{1}{{D}_1-z}
\mathcal{V}
}_{\!L},\end{aligned}$$ where $\mathcal{O}$ is the order of the pole at $z=0$. The reduced expressions for normal diagrams are identical to the usual expressions of zero-temperature MBPT, except that the step functions are replaced by Fermi-Dirac distributions. As a consequence, while at ${T=0}$ the energy denominator poles in these expressions are at the integration boundary, for ${T\neq 0}$ they are in the interior. This entails that the reduced expressions for individual diagrams are not well-defined for ${T\neq 0}$.

Last, we note that each of the time-independent formulas can be applied also to unlinked diagrams (the only change being the omission of the subscript $L$); this will become relevant in Sec. \[sec4\].

### Classification of perturbative contributions

Anomalous diagrams give no contribution in zero-temperature MBPT. However, the contributions from anomalous diagrams in grand-canonical MBPT do *not* vanish for ${T\rightarrow 0}$ (in the thermodynamic limit). The reduced integrands (which are well-defined at ${T=0}$) for diagrams with identically vanishing energy denominators[^10] have terms of the form $$\begin{aligned}
 \label{anomT0lim}
\frac{\partial^\nu n_{\textbf}{k}}{\partial \mu^\nu} \xrightarrow{T\rightarrow 0}
\delta^{(\nu)}(\mu-{\varepsilon}_{\textbf}{k}),\end{aligned}$$ e.g., $\beta  n_{\textbf}{k} \bar n_{\textbf}{k}=\partial n_{\textbf}{k}/\partial \mu  \xrightarrow{T\rightarrow 0}
\delta(\mu-{\varepsilon}_{\textbf}{k})$. Contributions with such terms are called anomalous contributions. There are also contributions that vanish for ${T\rightarrow 0}$, e.g., $$\begin{aligned}
n_{\textbf}{k} \bar n_{\textbf}{k}=T\frac{\partial n_{\textbf}{k}}{\partial \mu} \xrightarrow{T\rightarrow 0}
0.\end{aligned}$$ Such pseudoanomalous contributions can be associated also with *normal* two-particle reducible diagrams via the relation $$\begin{aligned}
\label{ndouble}
\bar n_{\textbf}{k}=1-n_{\textbf}{k},\end{aligned}$$ i.e., $$\begin{aligned}
\label{ndouble1}
n_{\textbf}{k} n_{\textbf}{k} &= n_{\textbf}{k} -n_{\textbf}{k} \bar n_{\textbf}{k}, 
\\ \label{ndouble2}
\bar n_{\textbf}{k} \bar n_{\textbf}{k} &= \bar n_{\textbf}{k} -n_{\textbf}{k} \bar n_{\textbf}{k},
\\ \label{ndouble3}
n_{\textbf}{k} n_{\textbf}{k} n_{\textbf}{k} &= n_{\textbf}{k} -2 n_{\textbf}{k} \bar n_{\textbf}{k}+n_{\textbf}{k} \bar n_{\textbf}{k} \bar n_{\textbf}{k},
\\ \label{ndouble4}
\bar n_{\textbf}{k} \bar n_{\textbf}{k} \bar n_{\textbf}{k} &= \bar n_{\textbf}{k} -2 n_{\textbf}{k} \bar n_{\textbf}{k}+n_{\textbf}{k} n_{\textbf}{k} \bar n_{\textbf}{k},\end{aligned}$$ etc.[^11] Contributions which are not anomalous or pseudoanomalous are referred to as normal contributions. Following loosely Balian, Bloch, and de Dominicis [@Balian1961529], we refer to the application of Eq.  according to Eqs. –, etc. as disentanglement, denoted symbolically by $\div$.

For the ${T\rightarrow 0}$ limit, the energy denominator exponentials present in the direct and cyclic formula all have to be evaluated via $$\begin{aligned}
 \label{Dexp}
\bar n_{\textbf}{k} \operatorname{\operatorname{e}}^{-\beta ({\varepsilon}_{\textbf}{k}-\mu)} = n_{\textbf}{k}.\end{aligned}$$ The simple relations given by Eqs. (\[ndouble\]) and (\[Dexp\]) play a crucial role in many of the issues and results discussed in the present paper.

Discrete spectrum inconsistency and anomalous contributions {#sec23a}
-----------------------------------------------------------

Apart from being essential for practical many-body calculations, the thermodynamic limit is in fact essential for the thermodynamic consistency of the grand-canonical perturbation series at low $T$, in particular for ${T\rightarrow 0}$, in the general case (see below). For a finite system with a discrete spectrum at ${T=0}$ one has either $\mu\in\{{\varepsilon}_{\textbf}{k}\}$ or $\mu\not\in\{{\varepsilon}_{\textbf}{k}\}$. Both cases are inconsistent.

In the first case the ${T\rightarrow 0}$ limit is singular, because in that case the anomalous contributions diverge. In addition, for discrete systems and $\mu\in\{{\varepsilon}_{\textbf}{k}\}$ the ${T\rightarrow 0}$ limit is singular due to energy denominator singularities.[^12]

In the second case the anomalous contributions vanish for ${T\rightarrow 0}$. From $F(T,\mu)=\Omega(T,\mu)+\mu\varrho(T,\mu)$ and the fact that all contributions to $\varrho(T,\mu)=-\partial\Omega(T,\mu)/\partial\mu$ except the ones from $\Omega_\text{ref}(T,\mu)$ are anomalous, for $\mu\not\in\{{\varepsilon}_{\textbf}{k}\}$ we obtain $$\begin{aligned}
 \label{gcfinite}
F(T,\mu)\xrightarrow{T\rightarrow 0} E^{(0)}(\mu), 
\;\;\;\;\;\;\;\;\;\;
\varrho(T,\mu)\xrightarrow{T\rightarrow 0}\sum_{{\textbf}{k}}\theta(\mu-{\varepsilon}_{\textbf}{k}),\end{aligned}$$ where $E^{(0)}({\varepsilon}_{\text{F}})$ corresponds to the adiabatic series. As noted by Kohn and Luttinger [@Kohn:1960zz], the two parts of Eq.  are inconsistent with each other. A possible definition of the chemical potential at ${T=0}$ in the finite case is $$\begin{aligned}
\mu(T=0,\varrho)=\frac{E^{(0)}(\varrho+1)+E^{(0)}(\varrho)}{2}.\end{aligned}$$ The second part of Eq.  however is equivalent to $$\begin{aligned}
\mu(T=0,\varrho)=\frac{E_\text{ref}^{(0)}(\varrho+1)+E_\text{ref}^{(0)}(\varrho)}{2} \equiv \mu_\text{ref}(T=0,\varrho),\end{aligned}$$ which contradicts the previous equation. For a given particle number the true chemical potential deviates from the chemical potential of the reference system, and Eq.  would imply that they are equal at $T=0$. Thus, in the discrete case the ${T\rightarrow 0}$ limit of the grand-canonical perturbation series is inconsistent also for $\mu\not\in\{{\varepsilon}_{\textbf}{k}\}$.

The same inconsistency can arise in the thermodynamic limit if the reference spectrum has a gap $\Delta$ and $\mu\not\in\{{\varepsilon}_{\textbf}{k}\}$. The ${T\rightarrow 0}$ limit is still smooth in the discrete (gapped) case, so the inconsistency is still present for nonzero $T$, although it is washed out at sufficiently high $T$. Qualitatively, in the discrete case the inconsistency is relevant if the spectrum does not resolve the anomalous terms $\partial^\nu n_{\textbf}{k}/\partial \mu^\nu$. As discussed in Sec. \[sec23\], contributions with such terms can be seen to account for the mismatch generated by using the reference spectrum together with the true chemical potential. If the anomalous terms are not sufficiently resolved the information about this mismatch gets lost and one approaches the paradoxical result that $\mu(T,\varrho)= \mu_\text{ref}(T,\varrho)$.[^13]

There are two ways the discrete (gapped) spectrum inconsistency for $\mu\not\in\{{\varepsilon}_{\textbf}{k}\}$ can be partially resolved, i.e.,

1.  by using the reference chemical potential instead of the true one,

2.  by choosing a mean-field that leads to $\mu({T=0},\varrho)= \mu_\text{ref}({T=0},\varrho)$.

Case (i) corresponds to the modified perturbation series $F(T,\mu_\text{ref})$. The partial resolution of the discrete spectrum inconsistency in that case is as follows:

1.  $F(T,\mu_\text{ref})$ involves additional anomalous contributions, and the failure to resolve the old ones is balanced (for anisotropic systems, only partially) by not resolving the new ones.

In the gapped case $F(T,\mu_\text{ref})$ reproduces the adiabatic series in the ${T\rightarrow 0}$ limit (if $\mu_\text{ref}\not\in\{{\varepsilon}_{\textbf}{k}\}$). In the gapless case the adiabatic series is reproduced only for isotropic systems; in that case the old and new anomalous contributions cancel for ${T\rightarrow 0}$. Thus, there is still a remainder of the discrete spectrum inconsistency. The information about anisotropy encoded in the anomalous contributions is not resolved in the discrete case at low $T$. In particular, for $F(T,\mu_\text{ref})$ the thermodynamic limit (and ${\Delta\rightarrow 0}$ limit, respectively) and the ${T\rightarrow 0}$ limit are noncommuting limits in the anisotropic case.

Regarding case (ii), there are three mean-field renormalization schemes that lead to $\mu(T,\varrho)= \mu_\text{ref}(T,\varrho)$, and accordingly, $F(T,\mu)=F(T,\mu_\text{ref})$: the direct, the cyclic, and the BdD scheme; see Sec. \[sec23\]. The anomalous diagrams are removed in each scheme, but in the direct and cyclic schemes there are still anomalous contributions. Hence, for the direct and cyclic schemes there is no discrete spectrum inconsistency despite anomalous contributions. However, these schemes are well-behaved only at high $T$ where the inconsistency ceases to be relevant. In particular, the ${T\rightarrow 0}$ limit does not exist for the direct and cyclic scheme.

The ${T\rightarrow 0}$ limit exists for the BdD scheme, but for ${N>2}$ this scheme exists only in the thermodynamic limit. The commutativity of the ${T\rightarrow 0}$ and ${\Delta \rightarrow 0}$ limits is fully restored in the BdD scheme, irrespective of isotropy. The anomalous contributions can be removed and the result $\mu({T=0},\varrho)= \mu_\text{ref}({T=0},\varrho)$ can be achieved also for finite systems, via the mean-field renormalization scheme specified by Eq.  below.[^14] For ${N\leq 2}$ the Eq.  scheme converges to the BdD scheme, but for ${N>2}$ it becomes ill-defined (singular, for ${N\geq 4}$) in the thermodynamic limit. Altogether, we have:

1.  The result $\mu({T=0},\varrho)= \mu_\text{ref}({T=0},\varrho)$ can be achieved for finite systems and in the thermodynamic limit, irrespective of isotropy, but for ${N>2}$ these two cases are not smoothly connected.

The commutativity of limits can however be fully restored for $F(T,\mu_\text{ref})$:

1.  For $F(T,\mu_\text{ref})$ together with the mean-field renormalization scheme given by Eq. \[Ufermi\] below the limit ${T\rightarrow 0}$ commutes with both the thermodynamic limit and the ${\Delta \rightarrow 0}$ limit, irrespective of isotropy. This is because for $F(T,\mu_\text{ref})$, at ${T=0}$ the Eq. \[Ufermi\] scheme removes the anomalous contributions.

Case (ii) and case (iii) both lead to the adiabatic formalism, irrespective of isotropy. There are however still anomalous contributions at finite $T$ in case (iii), and the reference chemical potential is identified with the true chemical potential only in case (ii).[^15]

Mean-field renormalization schemes {#sec23}
----------------------------------

The usual choices for the mean-field potential are ${U=0}$ (free reference spectrum) or ${U=U_1}$ (Hartree-Fock spectrum). In general, one expects that the choice ${U=U_1}$ leads to an improved perturbation series, compared to ${U=0}$. For first-order MBPT this can be seen from the fact that ${U=U_1(T,\mu)}$ and ${U=U_1(\varepsilon_{\text{F}})}$, respectively, are stationary points of the right-hand sides of the inequalities $$\begin{aligned}
\Omega(T,\mu) &\leq \Omega_\text{ref}(T,\mu) 
+\Omega_{U}(T,\mu) 
+ \Omega_1(T,\mu),
\\
E^{(0)}(\varepsilon_\text{F}) &\leq E^{(0)}_\text{ref}(\varepsilon_\text{F}) 
+E^{(0)}_{U}(\varepsilon_\text{F})
+ E^{(0)}_1(\varepsilon_\text{F}),\end{aligned}$$ for grand-canonical and adiabatic MBPT, respectively. For truncation orders $N> 1$, however, no similar formal argument is available for ${U=U_1}$ representing the best choice.

For both ${U=0}$ and ${U=U_1}$, in the thermodynamic limit the grand-canonical perturbation series does not reproduce the adiabatic one for $T\rightarrow 0$. The adiabatic series is also not reproduced in the discrete case, and in that case the grand-canonical series is inconsistent, in general (in particular, for ${U\in\{\,0,U_1\}}$); see Sec. \[sec23a\]).

It is now important to note that, at least for ${U=0}$, in general bare grand-canonical MBPT leads to deficient results also in the thermodynamic limit. This is particularly evident for a system with a first-order phase transition: for ${U=0}$ it is impossible to obtain the nonconvex single-phase constrained free energy from $\Omega(T,\mu)$, since $\Omega(T,\mu)$ is necessarily a single-valued function of $\mu$ for ${U=0}$; see also Refs. [@Fritsch:2002hp; @PhysRevC.89.064009; @Wellenhofer:2017qla].

This deficiency can be repaired by modifying the expression for $F(T,\mu)$ in terms of a (truncated) formal expansion about the chemical potential $\mu_\text{ref}\xrightarrow{T\rightarrow 0}{\varepsilon}_{\text{F}}$ of the reference system; see Sec. \[sec42\] for details. This expansion introduces additional contributions, and the structure of these contributions is very similar to anomalous diagrams. For isotropic systems it can be seen that the anomalous parts of these additional contributions cancel the old ones for ${T\rightarrow 0}$, leading to $$\begin{aligned}
F(T,\mu_\text{ref})\xrightarrow{T\rightarrow 0}E^{(0)}(\varepsilon_{\text{F}})\end{aligned}$$ in the isotropic case.[^16]

For $\Omega(T,\mu)$ with ${U=U_1(T,\mu)}$ one has $\varrho(T,\mu)=\sum_{\textbf}{k} n_{\textbf}{k}$ at truncation order ${N=1}$. Thus, for ${N=1}$ (but not for ${N>1}$) it is $$\begin{aligned}
F(T,\mu)=F(T,\mu_\text{ref}),\end{aligned}$$ with $\mu(T,\varrho)=\mu_\text{ref}(T,\varrho)$, where $F(T,\mu_\text{ref})$ now corresponds to the modified series with ${U=U_1(T,\mu_\text{ref})}$. For $\Omega(T,\mu)$ the change from ${U=0}$ to ${U=U_1}$ removes all anomalous (and normal) diagrams with single-vertex loops. For $F(T,\mu_\text{ref})$ both the reference spectrum and the reference chemical potential get renormalized, and both the anomalous diagrams and the additional ones with single-vertex loops are removed.

Now, these features make evident that there is a deficiency in the grand-canonical series with ${U=0}$ irrespective of the presence of a first-order phase transition: there is a mismatch in the Fermi-Dirac distribution functions generated by using the spectrum of $H_0$ together with the true chemical potential, leading to decreased perturbative convergence, compared to $F(T,\mu_\text{ref})$ with the same setup.[^17] One may interpret the anomalous contributions as a symptom of this mismatch. In that sense, the “expanding away” of the mismatch, i.e., the construction of $F(T,\mu_\text{ref})$, corresponds to a symptomatic treatment that provides as remedy additional anomalous contributions that counteract the old ones.[^18] The mismatch can however be ameliorated (cured, for $N=1$, in the $U=U_1$ case) by improving the quality of the reference Hamiltonian: the change from ${U=0}$ to ${U=U_1}$ removes the main symptom (and the corresponding remedy, in the modified case), the anomalous diagrams with single-vertex loops.

Altogether, this suggests that one can expect that the convergence behavior of $\Omega(T,\mu)$ is inferior to the one of $F(T,\mu_\text{ref})$ also for ${U=U_1}$. Moreover, one can suspect that both $\Omega(T,\mu)$ and $F(T,\mu_\text{ref})$ may be further improved by using a mean field beyond Hartree-Fock. In the best case, the additional mean-field contributions should remove all the remaining anomalous diagrams (and additional diagrams, for the modified series), i.e., the ones with higher-order pieces, and lead to $F(T,\mu)=F(T,\mu_\text{ref})$ for truncation orders ${N>1}$.

In the following, we introduce three different renormalization schemes where the mean field receives additional contributions for each $N$, i.e, $$\begin{aligned}
 \label{mbptren1}
U^{\aleph,(\ast\ast),\div}=U_1+\sum_{n=2}^N U^{\aleph,(\ast\ast),\div}_n.\end{aligned}$$ Here, $\aleph$ refers to one of the three time-independent formulas (direct, cyclic, or reduced), $\div$ to the disentanglement, and $\ast\ast$ to the regularization of energy denominators required to make the reduced formula well-defined. Equation  is understood to imply a reordering of the perturbation contributions such that a given order $n \in\{1,\ldots, N\}$ involves only diagrams for which $$\begin{aligned}
\label{Ncounting}
\mathscr{N}(V)+\mathscr{N}(U_1)+\sum_{m=2}^N m \mathscr{N}(U_m) =n \in\{1,\ldots, N\} ,\end{aligned}$$ where $\mathscr{N}(V)$ is the number of $V$ vertices, and $\mathscr{N}(U_1)$ and $\mathscr{N}(U_n)$ the number of $-U_1$ and $-U^{\aleph,(\ast\ast),\div}_n$ vertices, respectively. This can be implemented by writing Eq.  as $$\begin{aligned}
H = \underbrace{(U_0+U^{\aleph,(\ast\ast),\div})}_{H_\text{ref}} + \lambda V - \lambda U_1-\sum_{n=2}^N\lambda^n  U^{\aleph,(\ast\ast),\div}_n\end{aligned}$$ and ordering the perturbation series with respect to powers of $\lambda$ (which is at the end set to $\lambda=1$).

For each truncation order $N$, the three schemes constitute three (different, for ${N>1}$) stationary points of MBPT. Related to this, in each of the three schemes the (direct, cyclic and reduced, respectively) contributions from anomalous diagrams are removed, and in each scheme the relation between the particle number and the chemical potential matches the adiabatic relation, i.e., $$\begin{aligned}
U=U^{\aleph,(\ast\ast),\div}:\;\;\;\;\; 
\varrho(T,\mu)=
-\frac{\partial \Omega(T,\mu)}{\partial\mu}=\sum_{{\textbf}{k}} n_{\textbf}{k},\end{aligned}$$ so $F(T,\mu)=F(T,\mu_\text{ref})$ holds in each scheme. The ${T\rightarrow 0}$ limit exists however only for the case where $U=U^{\text{reduced},\ast\ast,\div}$. In that case, the grand-canonical formalism and zero-temperature MBPT are consistent with each other for both isotropic and anisotropic systems, given that the adiabatic continuation is based on $H_\text{ref}=H_0+U^{\text{reduced},\ast\ast,\div}$.

### Scheme by Balian, Bloch, and de Dominicis (direct scheme)

In the renormalization scheme by Balian, Bloch, and de Dominicis [@Balian1961529], the mean-field potential is, for truncation order $N$, defined as $$\begin{aligned}
 \label{Udirect}
U^{\text{direct},\div}_{{\textbf}{k}} = 
\sum_{n=1}^{N} 
U^{\text{direct},\div}_{n,{\textbf}{k}}
&=
\sum_{n=1}^{N} \frac{\delta 
\Omega_{n,\text{normal}}^{\text{direct},\div}}{\delta n_{\textbf}{k}}
= 
\frac{\delta 
\mathcal{D}^{\text{direct},\div}}{\delta n_{\textbf}{k}},\end{aligned}$$ i.e., only the direct contributions from normal diagrams are included, and $\div$ (disentanglement) means that for each set of (normal) articulation lines with identical three-momenta only one (hole or particle) distribution function appears \[i.e., only the first term of the right-hand sides of Eqs. –, etc., is included\]. The ${n=1}$ contribution to the mean field corresponds (as in the other schemes) to the usual Hartree-Fock single-particle potential, i.e., $$\begin{aligned}
 \label{UHF}
U_{1,{\textbf}{k}} =
\sum_{{\textbf}{k}'}
\braket{\psi_{{\textbf}{k}}\psi_{{\textbf}{k}'}| V|\psi_{{\textbf}{k}} \psi_{{\textbf}{k}'}}
n_{{\textbf}{k}'},\end{aligned}$$ where antisymmetrization is implied. For the higher-order contributions, the functional derivative $\delta/\delta n_{\textbf}{k}$ has to be evaluated *without* applying Eq. , i.e., the energy denominator exponentials have to be kept in the form that results from the contour integral. Otherwise, the functional derivative would be ill-defined (due to the emergence of poles). For ${n=2}$, one finds $$\begin{aligned}
 \label{U2dir}
U^{\text{direct},(\div)}_{2,{\textbf}{k}} &= 
\frac{1}{2}
\sum_{{\textbf}{k}_2, {\textbf}{k}_3, {\textbf}{k}_4} \!
|\braket{\psi_{{\textbf}{k}}\psi_{{\textbf}{k}_2}| V|\psi_{{\textbf}{k}_3} \psi_{{\textbf}{k}_4}}|^2 
{\nonumber \\}& \quad \times
\Big[ n_{{\textbf}{k}_2} \bar n_{{\textbf}{k}_3} \bar n_{{\textbf}{k}_4}\mathcal{F}^\text{direct}(D) 
- n_{{\textbf}{k}_3} n_{{\textbf}{k}_4} \bar n_{{\textbf}{k}_2} \mathcal{F}^\text{direct}(-D) \Big] ,\end{aligned}$$ where $$\begin{aligned}
\label{F2T0lim}
\mathcal{F}^\text{direct}(D)
&= \frac{1-\beta D-\operatorname{\operatorname{e}}^{-\beta D}}{\beta D^2} \xrightarrow{D\rightarrow 0} -\frac{\beta}{2},\end{aligned}$$ with $D={\varepsilon}_{{\textbf}{k}_3}+{\varepsilon}_{{\textbf}{k}_4} - {\varepsilon}_{{\textbf}{k}_2} - {\varepsilon}_{{\textbf}{k}}$. The ${T\rightarrow 0}$ limit of Eq.  is singular, due to the energy denominator exponential in $\mathcal{F}^\text{direct}(D)$: the functional derivative has removed one distribution function in the integrand, inhibiting the complete elimination of the energy denominator exponential via Eq. . Hence, the renormalization scheme of Balian, Bloch, and de Dominicis is of interest only for systems which are sufficiently close to the classical limit.[^19]

The direct contributions from anomalous diagrams composed of two normal pieces that are not (but may involve) $-U$ vertices have the factorized form $$\begin{aligned}
 \label{directfactorized}
\Omega_{n_1+n_2,\text{anomalous}}^{\text{direct},\div}
&= -\frac{\beta}{2} \sum_{\textbf}{k} 
U^{\text{direct},\div}_{n_1,{\textbf}{k}} n_{\textbf}{k} \bar n_{\textbf}{k} \,
U^{\text{direct},\div}_{n_2,{\textbf}{k}}
(2-\delta_{n_1,n_2})
,\end{aligned}$$ and similar for anomalous diagrams with several normal (non $-U$) pieces; see Sec. \[sec42\]. Given that for normal diagrams with $-U$ vertices the functional derivative in Eq.  acts only on the diagrammatic lines,[^20] Eq.  implies that the direct contributions from these diagrams are all canceled by the contributions from the corresponding diagrams with $-U$ pieces. The resulting perturbation series is then given by $$\begin{aligned}
 \label{MBPTrenorm1a}
\Omega(T,\mu) &= \Omega_\text{ref}(T,\mu) + 
\Omega_{U}(T,\mu) +
\mathcal{D}^{\text{direct},\div}(T,\mu).\end{aligned}$$ Using ${\varepsilon}_{{\textbf}{k}}={\varepsilon}_{0,{\textbf}{k}}+U^{\text{direct},\div}_{{\textbf}{k}}$, Eq.  can be written in the equivalent form $$\begin{aligned}
 \label{MBPTrenorm1}
\Omega[n_{\textbf}{k},T] &=T 
\sum_{\textbf}{k} \big( n_{\textbf}{k}\ln n_{\textbf}{k} +  
\bar n_{\textbf}{k}\ln \bar n_{\textbf}{k}\big)
+\sum_{\textbf}{k} \left({\varepsilon}_{0,{\textbf}{k}} 
- \mu \right) n_{\textbf}{k} 
{\nonumber \\}&  \quad + \mathcal{D}^{\text{direct},\div}[n_{\textbf}{k},T],\end{aligned}$$ which, using Eqs.  and , can be seen to be stationary under variations of the distribution functions, $\delta \Omega[n_{\textbf}{k},T]/\delta n_{\textbf}{k} =0$. From this one readily obtains the following expressions for the fermion number $\varrho$, the entropy $S$, and the internal energy $E$: $$\begin{aligned}
 \label{StatQP1a}
\varrho &= 
\sum_{\textbf}{k} n_{\textbf}{k},
\\ \label{StatQP2a}
S &=-\sum_{\textbf}{k} \big( n_{\textbf}{k}\ln n_{\textbf}{k} +  
\bar n_{\textbf}{k}\ln \bar n_{\textbf}{k}\big)
- \frac{\partial 
\mathcal{D}^{\text{direct},\div}[n_{\textbf}{k},T]}{\partial T} , 
\\
E &= \sum_{\textbf}{k} {\varepsilon}_{0,{\textbf}{k}} n_{\textbf}{k} +  
\mathcal{D}^{\text{direct},\div}[n_{\textbf}{k},T] 
- T\frac{\partial \mathcal{D}^{\text{direct},\div}[n_{\textbf}{k},T]}{\partial T}.\end{aligned}$$ The variation of the internal energy $\delta E[n_{\textbf}{k},T]/\delta  n_{\textbf}{k}$ is given by $$\begin{aligned}
 \label{StatQP3a}
\frac{\delta E}{\delta n_{\textbf}{k}} &=  {\varepsilon}_{{\textbf}{k}} 
-T\frac{\partial U^{\text{direct},\div}_{{\textbf}{k}}[n_{\textbf}{k},T]}{\partial T}.\end{aligned}$$ The relations given by Eqs. , and match those of Fermi-liquid theory [@Landau; @Landau2; @Landau3], except for the terms due to the explicit temperature dependence of $\mathcal{D}^{\text{direct},\div}[n_{\textbf}{k},T]$ and $U^{\text{direct},\div}_{\textbf}{k}[n_{\textbf}{k},T]$.

### Cyclic scheme

There is a straightforward variant of the scheme by Balian, Bloch, and de Dominicis: the cyclic scheme, with mean-field potential $$\begin{aligned}
 \label{Ucyclic}
U^{\text{cyclic},\div}_{{\textbf}{k}} = 
\sum_{n=1}^{N} 
U^{\text{cyclic},\div}_{n,{\textbf}{k}}
&=
\sum_{n=1}^{N} \frac{\delta 
\Omega_{n,\text{normal}}^{\text{cyclic},\div}}{\delta n_{\textbf}{k}}
= 
\frac{\delta 
\mathcal{D}^{\text{cyclic},\div}}{\delta n_{\textbf}{k}}.\end{aligned}$$ At second order one has $$\begin{aligned}
 \label{U2cyc}
U^{\text{cyclic},(\div)}_{2,{\textbf}{k}} &= 
-\frac{1}{4}
\sum_{{\textbf}{k}_2, {\textbf}{k}_3, {\textbf}{k}_4} \!
|\braket{\psi_{{\textbf}{k}}\psi_{{\textbf}{k}_2}| V|\psi_{{\textbf}{k}_3} \psi_{{\textbf}{k}_4}}|^2 
{\nonumber \\}& \quad \times
\Big[ n_{{\textbf}{k}_2} \bar n_{{\textbf}{k}_3} \bar n_{{\textbf}{k}_4} \mathcal{F}^\text{cyclic}(D)
-n_{{\textbf}{k}_3} n_{{\textbf}{k}_4} \bar n_{{\textbf}{k}_2} \mathcal{F}^\text{cyclic}(-D) \Big],\end{aligned}$$ where $$\begin{aligned}
\label{F2cycT0lim}
\mathcal{F}^\text{cyclic}(D)
&= \frac{1-\operatorname{\operatorname{e}}^{-\beta D}}{ D} \xrightarrow{D\rightarrow 0} \beta.\end{aligned}$$ In the cyclic scheme, the perturbation series and thermodynamic relations have the same structure as in the direct scheme. In particular, the same factorization property holds (see Sec. \[sec42\]), and again the zero-temperature limit does not exist \[as evident from Eq. \]. The direct scheme is, however, distinguished from the cyclic scheme in terms of it leading to the identification of the Fermi-Dirac distribution functions with the exact mean occupation numbers [@Balian1961529; @Balian1961529b; @boer] (see also Appendix \[app22\]) and (in the classical limit) the virial expansion [@Balian1961529b]. This indicates that, for calculations close to the classical limit, the direct scheme is preferable to the cyclic scheme.

### Reduced scheme(s)

In the renormalization scheme outlined by Balian and de Dominicis (BdD) [@statquasi3; @statquasi1], the term $\mathcal{D}^{\text{direct},\div}[n_{\textbf}{k},T]$ (or, $\mathcal{D}^{\text{cyclic},\div}[n_{\textbf}{k},T]$) is replaced by a term $\mathcal{D}^\text{BdD}[n_{\textbf}{k}]$ that has no explicit temperature dependence in addition to the one given by the functional dependence on $n_{\textbf}{k}(T,\mu)$, and satisfies $$\begin{aligned}
 \label{T0BdD}
\mathcal{D}^\text{BdD}(T,\mu) \xrightarrow{T\rightarrow 0} \sum_{n=1}^N E^{(0)}_n({\varepsilon}_{\text{F}}),\end{aligned}$$ where, by Eq. , $\mu\xrightarrow{T\rightarrow 0}{\varepsilon}_{\text{F}}$, and $E^{(0)}_n({\varepsilon}_{\text{F}})$ corresponds to the sum of all contributions of order $n$ in zero-temperature MBPT. This implies consistency with the adiabatic zero-temperature formalism irrespective of isotropy. The BdD mean field is given by $$\begin{aligned}
 \label{UBdDfirst}
U^\text{BdD}_{\textbf}{k}[n_{\textbf}{k}]=\frac{\delta\mathcal{D}^\text{BdD}[n_{\textbf}{k}]}{\delta n_{\textbf}{k}}.\end{aligned}$$ Since $\mathcal{D}^\text{BdD}[n_{\textbf}{k}]$ is supposed to have no explicit temperature dependence, it must be constructed by eliminating all energy denominator exponentials via Eq. . But then the functional derivative will lead to poles. To make the functional derivative well-defined, the energy denominators have to be regularized.

Now, as first recognized by Balian and de Dominicis [@BALIAN1960502] as well as Horwitz, Brout and Englert [@brout1], for a finite system with a discrete spectrum the following renormalized perturbation series can be constructed $$\begin{aligned}
 \label{MBPTrenorm2a}
\Omega(T,\mu) &= \Omega_\text{ref}(T,\mu) + 
\Omega_{U}(T,\mu) +
\mathcal{D}^{\text{reduced,}\ast,\div}(T,\mu),\end{aligned}$$ with mean field $$\begin{aligned}
  \label{Uredast}
U^{\text{reduced,}\ast\ast,\div}_{\textbf}{k}=\frac{\delta \mathcal{D}^{\text{reduced,}\ast,\div}}{\delta n_{\textbf}{k}},\end{aligned}$$ where $$\begin{aligned}
\label{MBPTrenorm2b}
\mathcal{D}^{\text{reduced,}\ast,\div}(T,\mu) 
=
\sum_{n=1}^{N} \frac{\delta 
\Omega_{n,\text{normal}}^{\text{reduced,}\ast,\div}}{\delta n_{\textbf}{k}}\xrightarrow{T\rightarrow 0} \sum_{n=1}^N E^{(0)}_n({\varepsilon}_{\text{F}}),\end{aligned}$$ with $\mu\xrightarrow{T\rightarrow 0}{\varepsilon}_{\text{F}}$. Here, ${\ast}$ means that the energy denominator poles are excluded in the discrete state sums (which makes the reduced formula well-defined, for a finite system). Equation  entails another factorization property, i.e., (see Sec. \[sec43\]) $$\begin{aligned}
\label{reducedfactorized1}
\Omega_{n_1+n_2,\text{anomalous}}^{\text{reduced,}\ast,\div}
&=  -\frac{\beta}{2}   \sum_{\textbf}{k} 
U^{\text{reduced,}\ast,\div}_{n_1,{\textbf}{k}} n_{\textbf}{k} \bar n_{\textbf}{k} 
U^{\text{reduced,}\ast,\div}_{n_2,{\textbf}{k}}
{\nonumber \\}&\quad \times
(2-\delta_{n_1,n_2})
.\end{aligned}$$ In Eq. , $\div$ implies that the pseudoanomalous terms from the reduced expressions for normal two-particle reducible diagrams with the same pieces are added (to the reduced expressions for the corresponding anomalous diagrams).

Equations  and lead to the Fermi-liquid relations for $\varrho$, $S$, and $\delta E/\delta n_{\textbf}{k}$. The validity of the ${\ast}$ prescription for finite systems is however somewhat questionable, since it disregards the contributions from the energy denominator poles present in the cyclic and direct case for ${T\neq 0}$ \[see Eqs.  and \].[^21] In the thermodynamic limit the contributions from energy denominator poles have measure zero. However, the thermodynamic limit of $\mathcal{D}^{\text{reduced,}\ast,\div}$ is singular at ${T\neq 0}$, due to terms with energy denominator poles of even degree.[^22] In addition, there are terms with several (odd) energy denominator poles for which the thermodynamic limit is not well-defined, as evident from the Poincaré-Bertrand transformation formula Eq. ; this implies that in the thermodynamic limit $U^{\text{reduced,}\ast,\div}_{\textbf}{k}$ is ill-defined also at ${T=0}$.

All in all, Eqs. \[MBPTrenorm2a\], , and indicate that the BdD renormalization scheme should correspond to $$\begin{aligned}
 \label{bddassume}
\mathcal{D}^\text{BdD}[n_{\textbf}{k}] 
=
\mathcal{D}^{\text{reduced,}\ast\ast,\div}[n_{\textbf}{k}],\end{aligned}$$ where ${\ast\ast}$ refers to the energy denominator regularization for infinite systems.

Statistical versus dynamical quasiparticles {#sec24}
-------------------------------------------

The statistical quasiparticles associated with the BdD renormalization scheme are distinguished from the dynamical quasiparticles [@Noz1; @Noz2; @Benfatto2006] associated with the asymptotic stability of the low-lying excited states. In the following, we examine how dynamical quasiparticles arise in grand-canonical MBPT, and compare their energies to the ones of the statistical quasiparticles (i.e., the single-particle energies in the BdD scheme). More details on the (various forms of the) self-energy are given in Appendix \[app2\]. In particular, in Appendix \[app22\] we show that (only) in the direct scheme the exact mean occupation numbers $f_{{\textbf}{k}}(T,\mu)$ are identified with the Fermi-Dirac distributions. Note that since the ${T\rightarrow 0}$ limit does not exist for the direct scheme, this result is consistent with the discontinuity of $f_{{\textbf}{k}}(T,\mu)$ at ${T=0}$. The consistency of $f_{{\textbf}{k}}({T\neq 0},\mu)=n_{{\textbf}{k}}({T\neq 0},\mu)$ with the results discussed below is examined in Appendix \[app22\].

### Dynamical quasiparticles without mean field

In MBPT (for normal systems), dynamical quasiparticles arise as follows. The perturbative contributions $\Sigma_{n,{\textbf}{k}}(z,T,\mu)$ to the frequency-space self-energy $\Sigma_{{\textbf}{k}}(z,T,\mu)$ are given by a specific analytic continuation (see Appendix \[app21\]) of the perturbative contributions to the Matsubara self-energy $\Xi_{{\textbf}{k}}(z_l,T,\mu)$, where $$\begin{aligned}
 \label{Matsubarafreq}
z_l=\frac{{\text{i}}(2l+1)\pi}{\beta}+\mu\end{aligned}$$ are the Matsubara frequencies, with $l\in\mathbb{Z}$. For example, in bare MBPT (with ${U=0}$) the two-particle irreducible second-order contribution to $\Xi_{{\textbf}{k}}(z_l,T,\mu)$ is given by \[see Eq. \] $$\begin{aligned}
 \label{sigma2}
\Xi_{2,{\textbf}{k}}(z_l,T,\mu) &= 
\frac{1}{2}
\sum_{{\textbf}{k}_2, {\textbf}{k}_3, {\textbf}{k}_4} \!
|\braket{\psi_{{\textbf}{k}}\psi_{{\textbf}{k}_2}| V|\psi_{{\textbf}{k}_3} \psi_{{\textbf}{k}_4}}|^2  n_{{\textbf}{k}_2} \bar n_{{\textbf}{k}_3} \bar n_{{\textbf}{k}_4}
{\nonumber \\}& \quad \times
\frac{e^{-\beta({\varepsilon}_{{\textbf}{k}_3}+{\varepsilon}_{{\textbf}{k}_4} - {\varepsilon}_{{\textbf}{k}_2} - z_l)}-1}{{\varepsilon}_{{\textbf}{k}_3}+{\varepsilon}_{{\textbf}{k}_4} - {\varepsilon}_{{\textbf}{k}_2} - z_l}.\end{aligned}$$ From this, the expression for $\Sigma_{2,{\textbf}{k}}(z,T,\mu)$ is obtained by *first* substituting $\operatorname{\operatorname{e}}^{\beta(z_l-\mu)}=-1$ and *then* performing the analytic continuation. Using Eq. , one gets $$\begin{aligned}
 \label{sigma2b}
\Sigma_{2,{\textbf}{k}}(z,T,\mu) &= 
-\frac{1}{2}
\sum_{{\textbf}{k}_2, {\textbf}{k}_3, {\textbf}{k}_4} \!
|\braket{\psi_{{\textbf}{k}}\psi_{{\textbf}{k}_2}| V|\psi_{{\textbf}{k}_3} \psi_{{\textbf}{k}_4}}|^2  
{\nonumber \\}& \quad \times
\frac{
n_{{\textbf}{k}_2} \bar n_{{\textbf}{k}_3} \bar n_{{\textbf}{k}_4} +
n_{{\textbf}{k}_3} n_{{\textbf}{k}_4} \bar n_{{\textbf}{k}_2}}
{{\varepsilon}_{{\textbf}{k}_3}+{\varepsilon}_{{\textbf}{k}_4} - {\varepsilon}_{{\textbf}{k}_2} - z}.\end{aligned}$$ As evident from the second-order contribution, setting $z = \omega\pm{\text{i}}\eta$, with $\omega$ real and $\eta$ infinitesimal, leads to the general relation [@kadanoffbaym] $$\begin{aligned}
 \label{dynquasi1a}
\Sigma_{\textbf}{k}(\omega\pm{\text{i}}\eta,T,\mu) = \mathcal{S}_{\textbf}{k}(\omega,T,\mu) \mp {\text{i}}\mathcal{J}_{\textbf}{k}(\omega,T,\mu),\end{aligned}$$ where $\mathcal{S}_{\textbf}{k}$ and $\mathcal{J}_{\textbf}{k}$ are real, and ${\mathcal{J}_{\textbf}{k}\geq 0}$ (see Appendix \[app21\]). From the property that at ${T=0}$ the energy denominators in the expressions for the perturbative contributions to the self-energy, $\Sigma_{n,{\textbf}{k}}(z,T,\mu)$, vanish only for $z\rightarrow \mu$, Luttinger [@PhysRev.121.942] showed that $$\begin{aligned}
 \label{dynquasi1}
\mathcal{J}_{\textbf}{k}(\omega,0,\mu) \xrightarrow{\omega\rightarrow \mu} C_{\textbf}{k}(\mu)\, (\omega-\mu)^2,\end{aligned}$$ with ${C_{\textbf}{k}(\mu)\geq 0}$. Crucial for our discussion (i.e., in particular for the next paragraph), this result holds not only if $\Sigma_{\textbf}{k}$ is calculated using self-consistent propagators but also if $\Sigma_{\textbf}{k}$ is calculated using bare propagators.

In Ref. [@PhysRev.119.1153], Luttinger showed that Eq.  implies a discontinuity at ${T=0}$ and ${{\textbf}{k}={\textbf}{k}_{\text{F}}}$ of the exact mean occupation numbers $f_{{\textbf}{k}}(T,\mu)$ of the momentum eigenstates $\ket{\psi_{{\textbf}{k}}}$, i.e., [@Luttinger:1960ua; @Parry] $$\begin{aligned}
f_{{\textbf}{k}}(T,\mu)=\braket{\!\braket{a_{\textbf}{k}^\dagger a_{\textbf}{k}}\!} 
= 
\int \limits_{-\infty}^\infty \!\frac{d\omega}{2\pi} \frac{1}{1+\operatorname{\operatorname{e}}^{\beta(\omega-\mu)}}\mathcal{A}_{\textbf}{k}(\omega,T,\mu),\end{aligned}$$ where $\braket{\!\braket{\ldots}\!}$ denotes the true ensemble average, and the spectral function $\mathcal{A}_{\textbf}{k}(\omega,T,\mu)$ is given by [@kadanoffbaym] (see also Appendix \[app21\]) $$\begin{aligned}
 \label{spectral}
\mathcal{A}_{\textbf}{k}(\omega,T,\mu)=\frac{2 \mathcal{J}_{\textbf}{k}(\omega,T,\mu)}{\left[\omega-{\varepsilon}_{0,{\textbf}{k}}-\mathcal{S}_{\textbf}{k}(\omega,T,\mu)\right]^2+\left[\mathcal{J}_{\textbf}{k}(\omega,T,\mu)\right]^2}.\end{aligned}$$ The (true) Fermi momentum ${\textbf}{k}_{\text{F}}$, defined in terms of the discontinuity of $f_{{\textbf}{k}}(0,\mu)$, is determined by [@PhysRev.119.1153] $$\begin{aligned}
 \label{dynquasi2}
\mu={\varepsilon}_{0,{\textbf}{k}_{\text{F}}} + \mathcal{S}_{{\textbf}{k}_{\text{F}}}(\mu,0,\mu).\end{aligned}$$ The lifetime of a single-mode excitation with momentum ${\textbf}{k}$ of the ground state is determined by the width of the spectral function at ${T=0}$ [@Fetter; @kadanoffbaym]. From Eqs.  and , the width vanishes (i.e., the excitation becomes stable against decay into collective modes) for ${\omega\rightarrow \mu}$ and ${{\textbf}{k}\rightarrow{\textbf}{k}_{\text{F}}}$. The energies $\mathcal{E}_{\textbf}{k}$ of the dynamical quasiparticles are therefore determined by $$\begin{aligned}
 \label{dynquasi2b}
\mathcal{E}_{\textbf}{k}(\mu)={\varepsilon}_{0,{\textbf}{k}} + \mathcal{S}_{{\textbf}{k}}(\mathcal{E}_{\textbf}{k},0,\mu),\end{aligned}$$ where ${\mathcal{E}_{\textbf}{k}\approx \mu}$ and ${{\textbf}{k}\approx{\textbf}{k}_{\text{F}}}$ (low-lying excitations).[^23]

### Dynamical quasiparticles with mean field

The distinction between the energies of statistical and dynamical quasiparticles can now be made explicit, in a specific sense. For bare MBPT with mean field $U_{\textbf}{k}(T,\mu)$ the self-energy is given by $$\begin{aligned}
\label{dynquasi3b}
\Sigma_{{\textbf}{k}}(z,T,\mu)
&=
-U_{\textbf}{k}(T,\mu) +\Sigma'_{{\textbf}{k}}(z,T,\mu).\end{aligned}$$ Here, the first term corresponds to the contribution from the self-energy diagram composed of a single $-U$ vertex. Since bare propagators are used, $\Sigma'_{{\textbf}{k}}(z,T,\mu)$ involves not only one- and two-particle irreducible but also two-particle reducible self-energy diagrams (including diagrams with $-U$ vertices); see, e.g., Ref. [@PLATTER2003250].[^24] It can be seen that \[see Eq. \] $$\begin{aligned}
\label{dynquasi4}
\Sigma'_{n,{\textbf}{k}}(z) =  
\left[\frac{\delta \Omega_n^{\text{reduced}}[n_{\textbf}{k}]}{\delta n_{\textbf}{k}}\bigg|_{{\varepsilon}_{\textbf}{k}= z}
\right]_{{\textbf}{k} \notin \{\text{articulation lines}\}},\end{aligned}$$ (with $\text{Im}[z]\neq 0$). Instead of Eq.  we have $$\begin{aligned}
\Sigma_{{\textbf}{k}}(z,T,\mu)
&=
-U_{\textbf}{k}(T,\mu) +\mathcal{S}'_{{\textbf}{k}}(z,T,\mu)+{\text{i}}\mathcal{J}'_{{\textbf}{k}}(z,T,\mu),\end{aligned}$$ with $$\begin{aligned}
 \label{dynquasi3d}
\mathcal{J}'_{\textbf}{k}(\omega,0,\mu) = C'_{\textbf}{k} (\omega-\mu)^2.\end{aligned}$$ The spectral function is now given by $$\begin{aligned}
 \label{spectralSJ}
\mathcal{A}_{\textbf}{k}(\omega,T,\mu)=\frac{2 \mathcal{J}'_{\textbf}{k}
(\omega,T,\mu)}{\Big[\omega-{\varepsilon}_{{\textbf}{k}}-\text{Re}\left[\Sigma_{{\textbf}{k}}(z,T,\mu)\right]\Big]^2+\left[\mathcal{J}'_{\textbf}{k}(\omega,T,\mu)\right]^2}.\end{aligned}$$ Using ${\varepsilon}_{{\textbf}{k}}={\varepsilon}_{0,{\textbf}{k}}+U_{{\textbf}{k}}(T,\mu)$, this becomes $$\begin{aligned}
\mathcal{A}_{\textbf}{k}(\omega,T,\mu)=\frac{2 \mathcal{J}'_{\textbf}{k}(\omega,T,\mu)}{\left[\omega-{\varepsilon}_{0,{\textbf}{k}}-\mathcal{S}'_{\textbf}{k}(\omega,T,\mu)\right]^2+\left[\mathcal{J}'_{\textbf}{k}(\omega,T,\mu)\right]^2},\end{aligned}$$ so the (true) Fermi-momentum ${\textbf}{k}_{\text{F}}$ is determined by $$\begin{aligned}
 \label{kfdef}
\mu={\varepsilon}_{0,{\textbf}{k}_{\text{F}}} + \mathcal{S}'_{{\textbf}{k}_{\text{F}}}(\mu,0,\mu),\end{aligned}$$ and the dynamical quasiparticle energies $\mathcal{E}_{\textbf}{k}$ are given by $$\begin{aligned}
 \label{quasiU}
\mathcal{E}_{\textbf}{k}(\mu)={\varepsilon}_{0,{\textbf}{k}} + \mathcal{S}'_{{\textbf}{k}}(\mathcal{E}_{\textbf}{k},0,\mu),\end{aligned}$$ where ${\mathcal{E}_{\textbf}{k}\approx \mu}$ and ${{\textbf}{k}\approx{\textbf}{k}_{\text{F}}}$. It is ${\mathcal{E}_{\textbf}{k} = {\varepsilon}_{\textbf}{k}}$ for ${N\leq 2}$ within the BdD renormalization scheme, but from Eq.  as well as Eqs.  and it is clear that this correspondence breaks down for truncation orders ${N>2}$. To have ${\mathcal{E}_{\textbf}{k} = {\varepsilon}_{\textbf}{k}}$ for ${N>2}$ the mean field must satisfy $$\begin{aligned}
\label{luttward}
U_{\textbf}{k}(0,\mu) = \mathcal{S}'_{{\textbf}{k}}(\mu,0,\mu),\end{aligned}$$ but then no statistical quasiparticle relations are obtained. In particular, formally extending Eq.  to momenta ${\textbf}{k}\in[0,{\textbf}{k}_{\text{F}}]$, the mean-field renormalization specified by Eq.  leads to $$\begin{aligned}
\sum_{\textbf}{k} \theta(\mu-\mathcal{E}_{\textbf}{k})
=
\sum_{\textbf}{k} \theta(\mu-{\varepsilon}_{\textbf}{k}),\end{aligned}$$ but the relation (i.e., Luttinger’s theorem [@PhysRev.121.942; @stefanuc; @Baym:1962sx]) $$\begin{aligned}
\sum_{\textbf}{k} \theta(\mu-\mathcal{E}_{\textbf}{k})=\varrho(T=0,\mu)\end{aligned}$$ is satisfied only for truncation orders $N\leq 2$.

Regularization of energy denominators {#sec3}
=====================================

An energy denominator regularization scheme is a procedure that allows to evaluate the contributions associated with the various parts ${\mathcal{F}_\alpha=f_\alpha/D}$ of the energy denominator terms ${\mathcal{F}=\sum_\alpha\mathcal{F}_\alpha}$ separately \[cf., e.g., Eq. \]. The (formal) splitting of the $\mathcal{F}$’s into parts $\mathcal{F}_\alpha$ introduces poles, so the essence of any regularization scheme must be a change in the way the contributions near the zeros of the denominators $D=\prod_\nu D_\nu^{n_\nu}$ of these terms are evaluated (in particular for the case where some $n_\nu$ are even). This change must be such that, for a fixed mean field, the same results are obtained as from the original unregularized expressions for the $\mathcal{F}$’s (e.g., the expressions obtained from the direct or cyclic formula).

For the second-order normal contribution the regularization is (essentially) unique and corresponds to evaluating the two parts of Eq.  separately via principal value integrals. For the higher-order contributions, the regularization scheme introduced here starts by adding infinitesimal imaginary parts to the individual energy denominators $D_\nu$, i.e., $\prod_\nu D_\nu^{n_\nu}\rightarrow \prod_\nu (D_\nu+i\eta_\nu)^{n_\nu}$. The regularization then corresponds to evaluating the various parts with energy denominator terms $\mathcal{F}_{\alpha,[\{\eta_\nu\}]}=f_\alpha/[\prod_\nu (D_\nu+i\eta_\nu)]^{n_\nu}$ via the Sokhotski-Plemelj-Fox formula. That this is a valid procedure can be seen from the fact that (after adding infinitesimal imaginary parts) the Sokhotski-Plemelj-Fox formula can be applied (formally) also to the unsplit expressions with energy denominator terms $\mathcal{F}_{[\{\eta_\nu\}]}=\sum_\alpha\mathcal{F}_{\alpha,[\{\eta_\nu\}]}$, and after its application the splitting corresponds again (i.e., as in the second-order case) to a separation into principal value integrals, by virtue of Eq.  below.

The crucial point of this particular regularization scheme is that it allows to separate the normal, anomalous, and pseudoanomalous contributions (at finite $T$) such that these contributions have a form that matches the (regularized) disentangled reduced formula. This feature is essential for the cancellation of the pseudoanomalous contributions and the factorization of the anomalous contributions, and these properties lead to the thermodynamic Fermi-liquid relations via the BdD scheme. In other terms, the Fermi-liquid relations uniquely determine the regularization of the energy denominators.[^25]

In Sec. \[sec31\] we introduce the formal approach to the energy denominator regularization for the BdD scheme.[^26] The numerical evaluation of the resulting expressions is discussed in Sec. \[sec32\].

Formal regularization {#sec31}
---------------------

From the cyclic expressions, the regularized (${\ast\ast}$) disentangled ($\div$) reduced expressions are obtained by performing the following steps:

1.  add infinitesimal imaginary parts $\eta_\nu$ to the energy denominators $D_\nu$ (where $\eta_1\neq \eta_2\neq \ldots$),

2.  eliminate the energy denominator exponentials via Eq. ,

3.  apply Eq. .

Here, the first step is part of $\ast\ast$, the second step is part of the reduction, and the third step is associated with $\div$. Then

1.  for two-particle reducible diagrams, average over the signs $\text{sgn}(\eta_\nu)$ of the imaginary parts,

2.  split the integrals such that the various parts of the cyclic energy denominator terms are integrated separately, then suitably relabel indices in some integrals, and finally recombine the integrals that lead to normal, pseudoanomalous and anomalous contributions,

3.  observe that the pseudoanomalous contributions vanish (this is proved to all orders in Sec. \[sec4\]),

4.  observe that the anomalous contributions factorize (this is proved to all orders in Sec. \[sec4\]),

where the first step is part of $\ast\ast$, and the second, third and fourth steps are associated with $\div$ and reduction. To show how these rules arise we now regularize, disentangle, and reduce the expressions for the contributions from the normal second-order diagram and from selected third-order diagrams.

![The normal second-order diagram. It is invariant under vertex permutations.[]{data-label="fig2normal"}](fig1.pdf){width="12.00000%"}

The cyclic expression for the normal second-order diagram shown in Fig. \[fig2normal\] is given by $$\begin{aligned}
\label{Omega2}
\Omega_{2,\text{normal}}^\text{cyclic}= -\frac{1}{8} \sum_{ijab}
\zeta^{ijab} n_{ij}\bar n_{ab} \frac{1-\operatorname{\operatorname{e}}^{-\beta D_{ab,ij}}}{D_{ab,ij}},\end{aligned}$$ where $\zeta^{ijab}=V^{ij,ab}V^{ab,ij}$, with $V^{ij,ab}=
\braket{\psi_{{\textbf}{k}_i}\psi_{{\textbf}{k}_j}|V|\psi_{{\textbf}{k}_a}\psi_{{\textbf}{k}_b} }$. Moreover, $\sum_i=\int d^3 k_i/(2\pi)^3$, $n_{ij}=n_{{\textbf}{k}_i}n_{{\textbf}{k}_j}$ and $\bar n_{ij}=(1-n_{{\textbf}{k}_a})(1-n_{{\textbf}{k}_b})$, and $D_{ab,ij}={\varepsilon}_{{\textbf}{k}_a}+{\varepsilon}_{{\textbf}{k}_b}-{\varepsilon}_{{\textbf}{k}_i}-{\varepsilon}_{{\textbf}{k}_j}$.

In Eq. , the term $(1-\operatorname{\operatorname{e}}^{-\beta D_{ab,ij}})/D_{ab,ij}$ is regular for ${D_{ab,ij}=0}$. To evaluate the two parts of the numerator of this term separately, we add an infinitesimal imaginary term ${\text{i}}\eta$ to the energy denominator. This leads to $$\begin{aligned}
\label{Omega2b}
\Omega_{2,\text{normal}}^\text{cyclic}&= -\frac{1}{8} \sum_{ijab}
\zeta^{ijab} n_{ij} \bar n_{ab} \frac{1-\operatorname{\operatorname{e}}^{-\beta D_{ab,ij}}}{D_{ab,ij}+{\text{i}}\eta}.
{\nonumber \\}&=
-\frac{1}{8} \sum_{ijab}
\zeta^{ijab} n_{ij} \bar n_{ab} \frac{1}{D_{ab,ij}+{\text{i}}\eta}
{\nonumber \\}& \quad +\frac{1}{8} \sum_{ijab}
\zeta^{ijab} n_{ab} \bar n_{ij} \frac{1}{D_{ab,ij}-{\text{i}}\eta}.\end{aligned}$$ where we have applied Eq.  to eliminate the energy denominator exponential in the second part. Relabeling indices $(i,j)\leftrightarrow (a,b)$ and recombining the two terms leads to $$\begin{aligned}
\label{Omega2c}
\Omega_{2,\text{normal}}^\text{cyclic}&=
-\frac{1}{8} \sum_{ijab}
\zeta^{ijab} n_{ij} \bar n_{ab} \left[\frac{1}{D_{ab,ij}+{\text{i}}\eta}+\frac{1}{D_{ab,ij}-{\text{i}}\eta}\right]
{\nonumber \\}& \equiv \Omega_{2,\text{normal}}^{\text{reduced,}\ast\ast,(\div)}=\Omega_{2,\text{normal}}^\text{BdD}.\end{aligned}$$ From this, one obtains for the second-order contribution to the BdD mean field the expression $$\begin{aligned}
\label{Omega2d}
U^\text{BdD}_{2,i} = U_{2,i}^{\text{reduced,}\ast\ast,(\div)} &=
-\frac{1}{4} \sum_{jab}
\zeta^{ijab} \big( n_{j} \bar n_{ab}+n_{ab} \bar n_{j}\big)
{\nonumber \\}& \quad \times
\left[\frac{1}{D_{ab,ij}+{\text{i}}\eta}+\frac{1}{D_{ab,ij}-{\text{i}}\eta}\right].\end{aligned}$$ Note that the expressions for $\Omega_{2,\text{normal}}^\text{BdD}$ and $U^\text{BdD}_{2,i}$ are real. Given that the integration variables include $D_{ab,ij}$ \[or an equivalent variable, see Eq. \], this can be seen explicitly from the Sokhotski-Plemelj theorem $$\begin{aligned}
 \label{plemelj}
\frac{1}{x+{\text{i}}\eta}= \frac{P}{x}-{\text{i}}\pi\,\text{sgn}(\eta)\,\delta(x),\end{aligned}$$ where $P$ refers to the Cauchy principal value. For actual numerical calculations it is however more practical not to use $D_{ab,ij}$ as an integration variable, and then the application of the Sokhotski-Plemelj theorem requires further attention. This issue is discussed in Sec. \[sec32\].

It will be useful now to examine how Eq.  can be derived from the direct formula. The direct expression is given by $$\begin{aligned}
 \label{Omega2dira}
\Omega_{2,\text{normal}}^\text{direct}= \frac{1}{4} \sum_{ijab}
\zeta^{ijab} n_{ij}\bar n_{ab} \frac{1-\operatorname{\operatorname{e}}^{-\beta D}-\beta D}{\beta D^2},\end{aligned}$$ where $D=D_{ab,ij}$. Adding an imaginary part to the energy denominator we have $$\begin{aligned}
 \label{Omega2dirb}
\Omega_{2,\text{normal}}^\text{direct}= \frac{1}{4} \sum_{ijab}
\zeta^{ijab} n_{ij}\bar n_{ab} \frac{1-\operatorname{\operatorname{e}}^{-\beta D}-\beta D}{\beta (D+{\text{i}}\eta)^2}.\end{aligned}$$ Here, the integral can be evaluated in terms of the Sokhotski-Plemelj-Fox formula [@fox] $$\begin{aligned}
\label{plemelj2}
\frac{1}{(x+{\text{i}}\eta)^n}= \frac{P}{x^n} + {\text{i}}\pi(-1)^{n} \,\text{sgn}(\eta)\,\delta^{(n-1)}(x),\end{aligned}$$ where now $P$ denotes the Hadamard finite part [@hadamard] (see also Refs. [@MONEGATO2009425; @galap; @dispersions]), i.e., $$\begin{aligned}
\label{hadamard}
\int \!\! dx\, \varphi(x)\frac{P}{x^{n+1}}\equiv
\frac{1}{n!}
\lim_{y\rightarrow 0}
\frac{\partial^{n}}{\partial y^{n}}\!\!
\int \!\! dx\, \varphi(x)\frac{P}{x-y},\end{aligned}$$ and $\delta^{(n-1)}(x)=\partial \delta(x)/\partial x^n$. Note that this prescription satisfies $x^k/(x+{\text{i}}\eta)^n=1/(x+{\text{i}}\eta)^{n-k}$. Since $\partial(1-\operatorname{\operatorname{e}}^{-\beta D}-\beta D)/\partial D=0$ for $D=0$, evaluating Eq.  with the Sokhotski-Plemelj-Fox formula gives the same result as Eq. . This equivalence is maintained if the three parts of the $1-\operatorname{\operatorname{e}}^{-\beta D}-\beta D$ are integrated separately (and evaluated with the Sokhotski-Plemelj-Fox formula). That is, applying *first* the Sokhotski-Plemelj-Fox formula and *then* Eq.  and the relabeling of indices we find $$\begin{aligned}
 \label{Omega2dirc}
\Omega_{2,\text{normal}}^\text{direct} &= -\frac{1}{4} \sum_{ijab}
\zeta^{ijab} n_{ij}\bar n_{ab} \, D\frac{P}{D^2}
=-\frac{1}{4} \sum_{ijab}
\zeta^{ijab} n_{ij}\bar n_{ab} \, \frac{P}{D}
{\nonumber \\}& \equiv \Omega_{2,\text{normal}}^{\text{reduced,}\ast\ast,(\div)}=\Omega_{2,\text{normal}}^\text{BdD}.\end{aligned}$$ It is now important to note that applying Eq.  and relabeling indices in the second part (which implies $D\rightarrow -D$) *before* applying the Sokhotski-Plemelj-Fox formula would lead to incorrect results, i.e., this procedure would leave the real part invariant but produce a finite imaginary part. This is because $$\begin{aligned}
 \label{Omega2wronga}
f(D)\frac{\operatorname{\operatorname{e}}^{-\beta D}}{(D+{\text{i}}\eta)^2}
&=f(D)\operatorname{\operatorname{e}}^{-\beta D}\frac{P}{D^2}+{\text{i}}\pi \,\text{sgn}(\eta)\,\delta(D)\,\beta f(0)
{\nonumber \\}& \quad
-{\text{i}}\pi \,\text{sgn}(\eta)\,\delta(D)\,\frac{\partial f(D)}{\partial D},\end{aligned}$$ whereas $$\begin{aligned}
\frac{f(-D)}{(-D+{\text{i}}\eta)^2}
=f(-D)\frac{P}{(-D)^2}+{\text{i}}\pi \,\text{sgn}(\eta)\,\delta(D)\,f(0).\end{aligned}$$ In general, for ${n>1}$ it is $$\begin{aligned}
f(D)\frac{\operatorname{\operatorname{e}}^{-\beta D}}{(D+{\text{i}}\eta)^n} \neq \frac{f(-D)}{(-D+{\text{i}}\eta)^n}.\end{aligned}$$ However, note that $f(-D)\frac{P}{(-D)^n}=f(D)\frac{P}{D^n}$ since $D$ is integrated in the whole real domain, and therefore $$\begin{aligned}
\text{Re}\left[f(D)\frac{\operatorname{\operatorname{e}}^{-\beta D}}{(D+{\text{i}}\eta)^n}\right] \equiv \text{Re}\left[\frac{f(-D)}{(-D+{\text{i}}\eta)^n}\right]\end{aligned}$$ for the considered $f(D)$. Hence, applying Eq.  and relabeling indices without first applying the Sokhotski-Plemelj-Fox formula becomes valid if we average over the sign of $\eta$, i.e., $$\begin{aligned}
\frac{1}{2}\sum_{\text{sgn}(\eta)}f(D)\frac{\operatorname{\operatorname{e}}^{-\beta D}}{(D+{\text{i}}\eta)^n} = \frac{1}{2}\sum_{\text{sgn}(\eta)} \frac{f(-D)}{(-D+{\text{i}}\eta)^n}.\end{aligned}$$ Note that the average has to be taken for all three parts of Eq. , otherwise imaginary parts would remain.

![The third-order two-particle irreducible diagrams. Each diagram is invariant under cyclic vertex permutations. The first (pp) and second (hh) diagram transform into each other under noncyclic permutations, and the third (ph) diagram is permutation invariant.[]{data-label="fig3ppph"}](fig2.pdf){width="45.00000%"}

The cyclic expressions for the third-order two-particle irreducible diagrams shown in Fig. \[fig3ppph\] are given by $$\begin{aligned}
 \label{Omega3pp}
\Omega_{3,\text{pp}}^\text{cyclic}&= \frac{1}{24} \sum_{ijabcd}
\zeta_{\text{pp}}^{ijabcd} n_{ij}\bar n_{abcd} 
\mathcal{F}_\text{pp}^\text{cyclic},
\\ \label{Omega3hh}
\Omega_{3,\text{hh}}^\text{cyclic}&= \frac{1}{24} \sum_{ijklab}
\zeta_{\text{hh}}^{ijklcd} n_{ijkl}\bar n_{ab} 
\mathcal{F}_\text{hh}^\text{cyclic},
\\ \label{Omega3ph}
\Omega_{3,\text{ph}}^\text{cyclic}&= \frac{1}{3} \sum_{ijkabc}
\zeta_{\text{ph}}^{ijkabc} n_{ijk}\bar n_{abc} 
\mathcal{F}_\text{ph}^\text{cyclic},\end{aligned}$$

where $\zeta_{\text{pp}}^{ijabcd}=V^{ij,ab}V^{ab,cd}V^{cd,ij}$, $\zeta_{\text{hh}}^{ijabcd}=V^{ij,ab}V^{kl,ij}V^{ab,kl}$, and $\zeta_{\text{ph}}^{ijabcd}=V^{ij,ab}V^{kb,ic}V^{ac,jk}$. The energy denominator terms are given by $$\begin{aligned}
\mathcal{F}_\text{pp,hh,ph}^\text{cyclic}
=
\left[\frac{1}{D_1 D_2}
+\frac{\operatorname{\operatorname{e}}^{-\beta D_1}}{D_1 (D_1 - D_2)}
-\frac{\operatorname{\operatorname{e}}^{-\beta D_2}}{D_2 (D_1 - D_2)}\right],\end{aligned}$$ with $D_1=D_{ab,ij}$ and $D_2=D_{cd,ij}$ for the pp diagram, $D_1=D_{ab,ij}$ and $D_2=D_{ab,kl}$ for the hh diagram, and $D_1=D_{ab,ij}$ and $D_2=D_{ac,jk}$ for the ph diagram. In each case, substituting $D_1\rightarrow D_1+{\text{i}}\eta_1$ and $D_2\rightarrow D_2+{\text{i}}\eta_2$, with $\eta_1\neq\eta_2$, splitting the integrals, eliminating the energy denominator exponentials and relabeling indices leads to $$\begin{aligned}
\mathcal{F}_\text{pp,hh,ph}^{\text{cyclic},\ast\ast} 
&=
\left[\frac{1}{(D_1+{\text{i}}\eta_1) (D_2+{\text{i}}\eta_2)}
+\frac{1}{(D_1-{\text{i}}\eta_1) (D_2+{\text{i}}\eta_2)}\right.
{\nonumber \\}&\quad  \left.
+\frac{1}{(D_1-{\text{i}}\eta_1) (D_2-{\text{i}}\eta_2)}\right] \equiv \mathcal{F}_\text{pp,hh,ph}^{\text{reduced},\ast\ast} ,\end{aligned}$$ which is real. Substituting this for $\mathcal{F}_\text{pp,hh,ph}^\text{cyclic}$ in Eqs. , , and  and performing the functional derivative one obtains the third-order contribution to $U^\text{BdD}$.

![The six third-order two-particle reducible diagrams composed of one second-order and one first-order piece. Articulation lines are shown as dashed lines. The shaded blobs represent vertices with loops (first-order pieces). In each row, the diagram on the left is a normal diagram, and the other two are anomalous. The diagrams in each row transform into each other under cyclic vertex permutations. The set of all six diagrams is closed under general vertex permutations.[]{data-label="fig3red"}](fig3.pdf){width="48.00000%"}

The normal third-order two-particle reducible diagrams are shown in Fig. \[fig3red\]. Also shown are the cyclically related anomalous diagrams. The cyclic expression for the sum of these diagrams is given by $$\begin{aligned}
\label{Omega3_21}
\Omega_{3,(21)}^\text{cyclic}&= -\frac{1}{4} \sum_{ijab}
\zeta^{ijab} n_{ij}\bar n_{ab} \mathcal{F}_{(21)}^\text{cyclic}
\big(n_i U_{1,i}-\bar n_a U_{1,a}\big),\end{aligned}$$ where $$\begin{aligned}
\mathcal{F}_{(21)}^\text{cyclic}
= \frac{1-\operatorname{\operatorname{e}}^{-\beta D}-\beta D\operatorname{\operatorname{e}}^{-\beta D}}{D^2}
,\end{aligned}$$ with $D=D_{ab,ij}$. In Hartree-Fock MBPT, the contribution from these diagrams is (as is well known) canceled by the corresponding diagrams where the first-order pieces are replaced by $-U_{1}$ vertices. Nevertheless, it will be still be useful to regularize these contributions. We will then find that, if $U_1$ were left out, the anomalous part of these diagrams can still be canceled via $U_2^\text{BdD}$.[^27]

Substituting $D\rightarrow D+{\text{i}}\eta$ and applying Eqs.  and and the relabeling $(i,j)\leftrightarrow (a,b)$, we can separate $\Omega_{3,(21)}^\text{cyclic}$ into the three contributions $$\begin{aligned}
\label{Omega3_21a}
\Omega_{3,(21),\text{normal}}^{\text{reduced,}\circ\circ,\div}&= -\frac{1}{4} \sum_{ijab}
\zeta^{ij,ab} n_{ij}\bar n_{ab} \frac{1}{(D+{\text{i}}\eta)^2}
\big(U_{1,i}-U_{1,a}\big),
\\ \label{Omega3_21b}
\Omega_{3,(21),\text{anom.}}^{\text{reduced,}\circ\circ,\div}&= \frac{\beta}{4} \sum_{ijab}
\zeta^{ij,ab} n_{ij}\bar n_{ab} \frac{1}{D+{\text{i}}\eta}
\big(\bar n_iU_{1,i}- n_aU_{1,a}\big),
\\ \label{Omega3_21c}
\Omega_{3,(21),\text{pseudo-a.}}^{\text{reduced,}\circ\circ,\div}&= \frac{1}{4} \sum_{ijab}
\zeta^{ij,ab} n_{ij}\bar n_{ab} \left[\frac{1}{(D+{\text{i}}\eta)^2}-\frac{1}{(D-{\text{i}}\eta)^2}\right]
{\nonumber \\}&\quad \times
\big(n_iU_{1,i}-\bar n_aU_{1,a}\big).\end{aligned}$$ Here, ${\circ\circ}$ refers to an incomplete (in fact, incorrect) regularization: none of the three contributions given by Eq. , and is real, and (more severely) also their sum is not real. As explained below Eq. , the reason for this deficiency is that we have applied Eq.  and relabeled indices without applying the Sokhotski-Plemelj-Fox formula first. To repair this we have to average over the signs of the imaginary parts, which leads to $$\begin{aligned}
 \label{T0canc}
\Omega_{3,(21),\text{normal}}^{\text{reduced,}\ast\ast,\div}&= -\frac{1}{8} \sum_{ijab}
\zeta^{ij,ab} n_{ij}\bar n_{ab} 
{\nonumber \\}&\quad 
\times \left[\frac{1}{(D+{\text{i}}\eta)^2}+\frac{1}{(D-{\text{i}}\eta)^2}\right]
 \big(U_{1,i}- U_{1,a}\big),
\\
\Omega_{3,(21),\text{anom.}}^{\text{reduced,}\ast\ast,\div}&= \frac{\beta}{8} \sum_{ijab}
\zeta^{ij,ab} n_{ij}\bar n_{ab} 
\left[\frac{1}{(D+{\text{i}}\eta)}+\frac{1}{(D-{\text{i}}\eta)}\right],
{\nonumber \\}&\quad 
\times (\bar n_iU_{1,i}- n_aU_{1,a}),
\\ \label{321pseudovanish}
\Omega_{3,(21),\text{pseudo-a.}}^{\text{reduced,}\ast\ast,\div}&= 0.\end{aligned}$$ The pseudoanomalous contribution has vanished: this feature, which is essential to obtain the Fermi-liquid relations at $T\neq 0$ (but not ${T=0}$), holds to all orders (see Sec. \[sec4\]). Note that the vanishing of the pseudoanomalous contributions holds only if all vertex permutations are included, i.e., it holds not separately for cyclically closed sets (in the present case, the two rows in Fig. \[fig3red\]).

The anomalous contribution has the factorized form given by Eq.  (with $\ast\ast$ instead of $\ast$), i.e., $$\begin{aligned}
\Omega_{3,(21),\text{anom.}}^{\text{reduced,}\ast\ast,\div}&= -\beta \,  U_{2,i}^{\text{reduced,}\ast\ast,(\div)} n_i \bar n_i
\, U_{1,i}.\end{aligned}$$ Thus, the anomalous contribution from the diagrams of Fig. \[fig3red\] gets canceled by the contribution from the diagram shown in Fig. \[fig4UU\] where one piece is a first-order diagram and the other one either a $-U_1$ vertex or a $-U_{2}^\text{BdD}$ vertex. The same cancellation occurs between the rotated diagram and the one with two mean-field vertices, and similar for the case where both $U_1$ and $U_{2}^\text{BdD}$ are included.

![The anomalous diagram composed of two pieces of the mean-field or single-vertex loop type.[]{data-label="fig4UU"}](fig4.pdf){width="14.00000%"}

Integration variables {#sec32}
---------------------

We now discuss how the formulas derived in Sec. \[sec31\] can be evaluated in numerical calculations. Nonvanishing contributions with poles of even degree appear first at fourth order in the BdD renormalization scheme. These have to be evaluated in terms of the Hadamard finite part, which obviously represents a major difficulty in the numerical application of the BdD scheme at high orders. We leave out the discussion of methods to evaluate the Hadamard finite part numerically, and defer numerical applications of the BdD scheme (and the other schemes) to future research.

For an isotropic system and MBPT without a mean-field potential (${U_{\textbf}{k}=0}$) where ${\varepsilon}_{\textbf}{k}={\varepsilon}_{0,{\textbf}{k}}={\textbf}{k}^2/(2M)$, using as integration variables relative momenta ${\textbf}{p}=({\textbf}{k}_i-{\textbf}{k}_j)/2$ and ${\textbf}{A}=({\textbf}{k}_a-{\textbf}{k}_b)/2$ as well as the average momentum ${\textbf}{K}=({\textbf}{k}_i+{\textbf}{k}_j)/2=({\textbf}{k}_a+{\textbf}{k}_b)/2$, one obtains from Eq.  the following expression for the second-order normal contribution: $$\begin{aligned}
\label{Omega2num1}
\Omega_{2,\text{normal}}^\text{BdD}&=
-2 M \sum_{{\textbf}{K},{\textbf}{p},{\textbf}{A}}
\zeta^{ijab} n_{ij} \bar n_{ab} \frac{P}{A^2-p^2}.\end{aligned}$$ The functional derivative of this expression with respect to $n_{{\textbf}{k}_i}$ is given by $$\begin{aligned}
\label{U2num1}
U_{2,{\textbf}{k}_i}^\text{BdD}[{\varepsilon}_{0,{\textbf}{k}_i}]=
-4 M \sum_{{\textbf}{p},{\textbf}{A}}
\zeta^{ijab} \big(n_{j} \bar n_{ab}+n_{ab} \bar n_{j}\big) \frac{P}{A^2-p^2}.\end{aligned}$$ For truncation order ${N=2}$, the single-particle energies in the BdD scheme are obtained from the self-consistent equation $$\begin{aligned}
{\varepsilon}_{\textbf}{k}=\frac{{\textbf}{k}^2}{2M}+ U_{1,{\textbf}{k}}+ U_{2,{\textbf}{k}}^\text{BdD}[{\varepsilon}_{\textbf}{k}],\end{aligned}$$ where one may use for $U_{2,{\textbf}{k}}^\text{BdD}[{\varepsilon}_{\textbf}{k}]$ the expression obtained by substituting in Eq.  the term $D_{ab,ij}={\varepsilon}_{{\textbf}{k}_a}+{\varepsilon}_{{\textbf}{k}_b}-{\varepsilon}_{{\textbf}{k}_i}-{\varepsilon}_{{\textbf}{k}_j}$ for $M/(A^2-p^2)$ if this substitution does not introduce additional poles; otherwise one must go back to the expression with infinitesimal imaginary parts, Eq. . This issue can be seen also in the ${U_{\textbf}{k}=0}$ case if ${\textbf}{k}_a$, ${\textbf}{k}_i$ and ${\textbf}{k}_j$ are used as integration variables to evaluate Eq. . Considering a one-dimensional system for simplicity, we have $$\begin{aligned}
\label{Omega2num2a}
\Omega_{2,\text{normal}}^\text{BdD}&=
-\frac{ M}{4}\!\! \sum_{k_a,k_i,k_j}\!\! \zeta^{ijab} n_{ij} \bar n_{ab} \,
\left[\frac{1}{\kappa+{\text{i}}\eta}
+\frac{1}{\kappa-{\text{i}}\eta} \right],\end{aligned}$$ with $k_b=k_i+k_j-k_a$ and $\sum_k=\int dk/(2\pi)$. Moreover, $\kappa=(k_a-k_i)(k_a-k_j)$, i.e., now there are two poles. To bring Eq.  into a form where the Sokhotski-Plemelj theorem can be applied, we note that $$\begin{aligned}
&(k_a-k_i+ {\text{i}}\eta)(k_a-k_j+{\text{i}}\eta)
{\nonumber \\}&\quad \quad=
(\kappa +{\text{i}}\eta) \, \theta(2k_a-k_i-k_j)
+(\kappa - {\text{i}}\eta) \, \theta(k_i+k_j-2k_a),\end{aligned}$$ so $$\begin{aligned}
\sum_{\text{sgn}(\eta)}\frac{1}{\kappa+{\text{i}}\eta}
=
\sum_{\text{sgn}(\eta)}
\frac{1}{(k_a-k_i+{\text{i}}\eta)(k_a-k_j+{\text{i}}\eta)}.\end{aligned}$$ The Sokhotski-Plemelj theorem can now be applied (assuming that $k_a$ is integrated after $k_i$ or $k_j$), which leads to $$\begin{aligned}
\label{Omega2num2b}
\Omega_{2,\text{normal}}^\text{BdD}&=
-\frac{ M}{2} \sum_{k_a,k_i,k_j} \zeta^{ijab} n_{ij} \bar n_{ab} 
{\nonumber \\}&\quad \times
\bigg[\frac{P}{k_a-k_i}\frac{P}{k_a-k_j}
+
\pi^2 \delta(k_a-k_i)\delta(k_a-k_j) \bigg],\end{aligned}$$ where the integration order is fixed. Changing the integration order such that $k_a$ is integrated first would lead to an incorrect result, as evident from the Poincaré-Bertrand transformation formula [@hardy; @poincare; @bertrand; @Muskheli; @dispersions] $$\begin{aligned}
 \label{bertrand}
&\int \!\! dx \!\int \!\! dy 
\,\,\varphi(x,y)
\frac{P}{x-y}\frac{P}{x-z} 
{\nonumber \\}&\quad = 
\int \!\! dy \!\int \!\! dx 
\,\,\varphi(x,y)
\frac{P}{x-y}\frac{P}{x-z} 
+
\pi^2 \varphi(z,z).\end{aligned}$$ Since it has only one pole, the expression given by Eq.  is however preferable compared to the one where ${\textbf}{k}_a$, ${\textbf}{k}_i$ and ${\textbf}{k}_j$ are used as integration variables.

At third order the issue manifested by the Poincaré-Bertrand transformation formula becomes unavoidable. For an isotropic system and ${U_{\textbf}{k}=0}$, using relative and average momenta as integration variables one obtains for $\Omega_{3,\text{pp}}^\text{BdD}$ the expression (see also Refs. [@Kondo; @YosidaHirosi]) $$\begin{aligned}
\label{Omega3ppnum}
\Omega_{3,\text{pp}}^\text{BdD}&=
\frac{M^2}{3} \!
\sum_{{\textbf}{K},{\textbf}{p},{\textbf}{A},{\textbf}{B}}
\zeta^{ijabcd}_\text{pp} n_{ij} \bar n_{abcd} 
{\nonumber \\}& \quad \times
\left[3\frac{P}{A^2-p^2}\frac{P}{B^2-p^2}+\pi^2 \frac{\delta(A-p)\delta(B-p)}{(A+p)(B+p)}\right],\end{aligned}$$ where ${\textbf}{p}=({\textbf}{k}_i-{\textbf}{k}_j)/2$, ${\textbf}{A}=({\textbf}{k}_a-{\textbf}{k}_b)/2$, ${\textbf}{B}=({\textbf}{k}_c-{\textbf}{k}_d)/2$ and ${\textbf}{K}=({\textbf}{k}_i+{\textbf}{k}_j)/2$. In Eq. , the integration order is such that $p$ is integrated after $A$ or $B$.[^28] The expression for $\Omega_{3,\text{hh}}^\text{BdD}$ is similar to Eq. . For $\Omega_{3,\text{ph}}^\text{BdD}$, however, using relative and average momenta as integration variables leads to $$\begin{aligned}
\label{Omega3phnum}
\Omega_{3,\text{ph}}^\text{BdD}&=
\frac{8 M^2}{3} 
\sum_{{\textbf}{K},{\textbf}{p},{\textbf}{A},{\textbf}{Y}}
\zeta^{ijkabc}_\text{ph} n_{ijk} \bar n_{abc}
{\nonumber \\}& \quad \times
\left[\mathcal{F}_{[\eta_1,\eta_2]}+\mathcal{F}_{[-\eta_1,\eta_2]}+\mathcal{F}_{[-\eta_1,-\eta_2]} \right],\end{aligned}$$ where ${\textbf}{p}=({\textbf}{k}_i-{\textbf}{k}_j)/2$, ${\textbf}{A}=({\textbf}{k}_a-{\textbf}{k}_b)/2$, ${\textbf}{Y}=({\textbf}{k}_a-{\textbf}{k}_c)/2$ and ${\textbf}{K}=({\textbf}{k}_i+{\textbf}{k}_j)/2$, and $$\begin{aligned}
\mathcal{F}_{[\eta_1,\eta_2]}=
\frac{1}{\big[A^2-p^2+{\text{i}}\eta_1\big]\, \big[({\textbf}{p}-{\textbf}{A})\cdot({\textbf}{A}-2{\textbf}{Y}+{\textbf}{p})+{\text{i}}\eta_2\big]}.\end{aligned}$$ From here one would have to proceed similar to the steps that lead from Eq.  to Eq. .

Factorization to All Orders {#sec4}
===========================

Here, we prove to all orders that the BdD renormalization scheme implies the thermodynamic relations associated with Fermi-liquid theory and (consequently) leads to a perturbation series that manifests the concistency of the adiabatic zero-temperature formalism, for both isotropic and anisotropic systems.

First, in Sec. \[sec41\], we examine more closely how the linked-cluster theorem manifests itself. Second, in Sec. \[sec42\] we systematize the disentanglement ($\div$) of the grand-canonical perturbation series. These two steps provide the basis for Sec. \[sec43\], where we prove to all orders the reduced factorization property for finite systems, Eq. . In Sec. \[sec44\] we then infer that the reduced factorization property holds also for the BdD renormalization scheme. This implies the Fermi-liquid relations and the consistency of the adiabatic formalism. Finally, in Sec. \[sec45\] we point out that the BdD renormalization scheme maintains the cancellation of the divergencies (at ${T=0}$) from energy denominator poles and discuss the minimal renormalization requirement for the consistency of the adiabatic formalism with the modified perturbation series for the free energy, $F(T,\mu_\text{ref})$, in the anisotropic case.

Linked-cluster theorem {#sec41}
----------------------

Letting the truncation order (formally) go to infinity, the sum of all perturbative contributions to $\Omega(T,\mu)$ can be written as $$\begin{aligned}
 \label{OmegaY}
\Delta \Omega = \sum_{n=1}^\infty \Omega_n
=-\frac{1}{\beta}\ln \left[1-\beta\sum_{n=1}^\infty \Upsilon_n \right],\end{aligned}$$ where $\Upsilon_n$ denotes the contribution of order $n$ from both linked and unlinked diagrams. We refer to the various linked parts of an unlinked diagram as subdiagrams. Further, we denote the contribution—evaluated via a given *time-independent* ($\aleph$) formula (i.e., direct, cyclic, or reduced with $\ast$ or $\ast\ast$)—to $\Upsilon_n$ from a diagram composed of $K=k\sum_{i=0}^{k} \alpha_i$ linked parts involving $k$ different subdiagram species $\Gamma_1\neq \Gamma_2 \neq \ldots \neq \Gamma_k$ where each $\Gamma_i$ appears $\alpha_i$ times in the complete diagram, by $\Upsilon^\aleph_{[\Gamma_{1}^{\alpha_1}\cdots\Gamma_{k}^{\alpha_k}]_n}$. In this notation, Eq.  reads $$\begin{aligned}
 \label{redfac1}
\Delta\Omega =
-\frac{1}{\beta} \ln \Big[1-\beta  \sum_{n=1}^\infty
\sum_{[\Gamma_{1}^{\alpha_1}\cdots\Gamma_{k}^{\alpha_k}]_n} 
\sum_O
\Upsilon^\aleph_{[\Gamma_{1}^{\alpha_1}\cdots\Gamma_{k}^{\alpha_k}]_n} \Big],\end{aligned}$$ where $\sum_{[\Gamma_{1}^{\alpha_1}\cdots\Gamma_{k}^{\alpha_k}]_n}$ is the sum over all possible (i.e., those consistent with order $n$) combinations of subdiagrams (including repetitions), and $\sum_O$ denotes the sum over all distinguishable vertex permutations of the unlinked diagram that leave the subdiagrams invariant. This is illustrated in Fig. \[figx\]. We write $$\begin{aligned}
 \label{Gammasum}
\sum_{n=1}^\infty   \sum_{[\Gamma_{1}^{\alpha_1}\cdots\Gamma_{k}^{\alpha_k}]_n}
\sum_O
\Upsilon^\aleph_{[\Gamma_{1}^{\alpha_1}\cdots\Gamma_{k}^{\alpha_k}]_n}
=
\sum_{\Gamma_{1}^{\alpha_1}\cdots\Gamma_{k}^{\alpha_k}}
\sum_O
\Upsilon^\aleph_{\Gamma_{1}^{\alpha_1}\cdots\Gamma_{k}^{\alpha_k}}.\end{aligned}$$ It is $$\begin{aligned}
 \label{redfac1b}
\sum_{\Gamma_{1}^{\alpha_1}\cdots\Gamma_{k}^{\alpha_k}} 
\sum_O
\Upsilon^\aleph_{\Gamma_{1}^{\alpha_1}\cdots\Gamma_{k}^{\alpha_k}} 
=
\sum_{\widetilde{\Gamma}_{1}^{\alpha_1}\cdots\widetilde{\Gamma}_{k}^{\alpha_k}} 
\sum_P
\Upsilon^\aleph_{\widetilde{\Gamma}_{1}^{\alpha_1}\cdots\widetilde{\Gamma}_{k}^{\alpha_k}} ,\end{aligned}$$ where $\sum_P$ denotes the sum over all distinguishable vertex orderings, and $\sum_{\widetilde{\Gamma}_{1}^{\alpha_1},
\ldots,\widetilde{\Gamma}_{k}^{\alpha_k}}$ sums over all combinations of subdiagrams where in the underlying set of linked diagrams $\{\Gamma_i\}$ only one (arbitrary) element is included for each set of diagrams that is closed under vertex permutations. For example, among the first two diagrams of Fig. \[fig3ppph\] only one is included, and only one of the six diagrams of Fig. \[fig3red\].

![Vertex permutations for an unlinked diagram with two linked parts (subdiagrams). If the diagram in the first row represents the original vertex ordering, then the second and third diagram correspond to (nonoverlapping and overlapping, respectively) orderings $\in O$, and the fourth diagram to an ordering $\in P/O$.[]{data-label="figx"}](fig5.pdf){width="25.00000%"}

The generalization of Eq.  for $\Upsilon_n$ is given by $$\begin{aligned}
 \label{OmegadirectT}  
\Upsilon_n^{\text{direct}[P]}=
-\frac{1}{\beta}\frac{(-1)^{n}}{n!} 
\int\limits_{0}^\beta \!
d \tau_n \cdots d \tau_1 \;
\Braket{  \mathcal{T}\big[
\mathcal{V}(\tau_n) \cdots
\mathcal{V}(\tau_1) \big]
}.\end{aligned}$$ We denote the expressions obtained from Eq.  for the contribution from a given permutation invariant set of (linked or unlinked) diagrams by $\Upsilon^{\text{direct}[P]}_{\widetilde{\Gamma}_{1}^{\alpha_1}\cdots
\widetilde{\Gamma}_{k}^{\alpha_k}}$. As noted in Sec. \[sec22\], these expressions are equivalent to the summed expressions obtained from any of the time-independent formulas (direct, cyclic, or reduced with $\ast$ or $\ast\ast$), i.e., $$\begin{aligned}
 \label{sumuppi}
\Upsilon^{\text{direct}[P]}_{\widetilde{\Gamma}_{1}^{\alpha_1}\cdots
\widetilde{\Gamma}_{k}^{\alpha_k}}
&=
\sum_P\Upsilon^\aleph_{\widetilde{\Gamma}_{1}^{\alpha_1}\cdots
\widetilde{\Gamma}_{k}^{\alpha_k}}.\end{aligned}$$ Now, the number of ways the $n$ perturbation operators in Eq.  can be partitioned into the subgroups specified by $\Upsilon^{\text{direct}[P]}_{\widetilde{\Gamma}_{1}^{\alpha_1}\cdots
\widetilde{\Gamma}_{k}^{\alpha_k}}$ is given by [@Abrikosov; @Fetter] $$\begin{aligned}
\label{partitioningcount}
\frac{1}{\alpha_1!
\cdots 
\alpha_k!}\frac{n!}{(n_1!)^{\alpha_1}
\cdots (n_k!)^{\alpha_k}},\end{aligned}$$ where $n_i$ are the orders of the respective subdiagrams. From Eq. , this leads to $$\begin{aligned}
 \label{redfac2}
\sum_P\Upsilon^\aleph_{\widetilde{\Gamma}_{1}^{\alpha_1}\cdots
\widetilde{\Gamma}_{k}^{\alpha_k}}
&= 
-\frac{1}{\beta} \prod_{i=1}^k
\frac{\left(-\beta \,
\Upsilon^{\text{direct}[P]}_{\widetilde{\Gamma}_{i}}\right)^{\alpha_i}}{\alpha_i!}
{\nonumber \\}&= 
-\frac{1}{\beta} \prod_{i=1}^k
\frac{\left(-\beta \,\sum_{P_i}
\Upsilon^\aleph_{\widetilde{\Gamma}_{i}}\right)^{\alpha_i}}{\alpha_i!},\end{aligned}$$ where in the second step we have applied Eq. . In Sec. \[sec43\] we will see that Eq.  implies the (direct, cyclic, and reduced) factorization properties for anomalous diagrams.

It is now straightforward to verify by explicit comparison with Eq.  that $$\begin{aligned}
 \label{redfac3}
 \sum_{\widetilde{\Gamma}_{1}^{\alpha_1}\cdots\widetilde{\Gamma}_{k}^{\alpha_k}} 
\sum_P
\left(-\beta\Upsilon^\aleph_{\widetilde{\Gamma}_{1}^{\alpha_1}\cdots\widetilde{\Gamma}_{k}^{\alpha_k}}\right)
&= 
\sum_{\alpha=1}^\infty \frac{1}{\alpha!}
\bigg[  -\beta 
\sum_{\Gamma}
\Upsilon^\aleph_{\Gamma}
\bigg]^\alpha 
{\nonumber \\}&=
-1+\exp\Big( -\beta \sum_{\Gamma} \Upsilon^\aleph_{\Gamma}\Big).\end{aligned}$$ Applying this to Eq.  leads to $$\begin{aligned}
\Delta \Omega &=
\sum_{\Gamma} \Upsilon^\aleph_{\Gamma},\end{aligned}$$ which constitutes the linked-cluster theorem. Now, expanding the logarithm in Eq.  we find $$\begin{aligned}
\label{Insert2}
\Delta \Omega&=
\sum_{n=1}^\infty  \, \Bigg[ \sum_{\nu,\{n_i\},\{k_i\}}
\beta^{k_1+\ldots+k_\nu-1}
\binom{k_1+\ldots+k_\nu}{k_1,\ldots,k_\nu}
\frac{(\Upsilon_{n_1})^{k_1}\cdots (\Upsilon_{n_\nu})^{k_\nu}}{k_1 +\ldots +k_\nu}\Bigg],\end{aligned}$$ where the inner sum is subject to the constraints $\sum_{i=1}^\nu n_i k_i=n$, $1\leq n_1<n_2<\ldots< n_\nu$, $k_i\geq 1$, and $\nu\geq 1$. The linked-cluster theorem implies that, if evaluated in terms of the usual Wick contraction formalism [@Fetter], in Eq.  the contributions with $\nu=1$ and $k_1=1$ from unlinked diagrams are all canceled by the contributions with $\nu>1$ or $k_1>1$. The individual expressions from these canceling terms are not size extensive, i.e., in the thermodynamic limit they diverge with higher powers of the confining volume.

Disentanglement {#sec42}
---------------

Here, we first introduce the cumulant formalism,[^29] which allows systematizing the disentanglement ($\div$). Then, we show that this formalism provides a new representation and evaluation method for the contributions associated (in the usual Wick contraction formalism) with anomalous diagrams and the subleading parts of Eqs. –, etc. Finally, we construct and discuss the modified perturbation series for the free energy $F(T,\mu_\text{ref})$.

### Cumulant formalism

We define $\mathcal{C}_{i_1 \ldots i_n}$ as the unperturbed ensemble average of a fully-contracted (indicted by paired indices) but not necessarily linked sequence of creation and annihilation operators, i.e., $$\begin{aligned}
\label{chi1def}
\mathcal{C}_{i_1 \ldots i_n}=\braket{a_{i_1}^\dagger a_{i_1} \cdots a_{i_n}^\dagger a_{i_n}},\end{aligned}$$ where some of the index tuples may be identical (articulation lines). In Eq. , all contractions are of the hole type. For the case where there are also particles we introduce the notation $$\begin{aligned}
\label{chi1defpart}
\mathcal{C}_{i_1 \cdots i_n}^{a_1 \cdots a_m}=\braket{a_{i_1}^\dagger a_{i_1} \cdots a_{i_n}^\dagger a_{i_n}
a_{a_1} a_{a_1}^\dagger \cdots a_{a_m} b_{a_m}^\dagger }.\end{aligned}$$ This can be expressed in terms of functional derivatives of the unperturbed partition function $\mathcal{Y}_{\!\text{ref}}=\operatorname{\operatorname{Tr}}[\operatorname{\operatorname{e}}^{-\beta(\mathcal{H}_\text{ref}-\mu \mathcal{N})}]$, i.e.,[^30] $$\begin{aligned}
 \label{chi1defpart2}
\mathcal{C}_{i_1 \cdots i_n}^{a_1 \cdots a_m}&=\frac{1}{\mathcal{Y}_{\!\text{ref}}} 
\frac{\delta}{\delta[-\beta {\varepsilon}_{i_1}]} \cdots 
\frac{\delta}{\delta[-\beta {\varepsilon}_{i_n}]}
{\nonumber \\}&\quad\times
\bigg(1-\frac{\delta}{\delta[-\beta {\varepsilon}_{a_1}]}\bigg)
\cdots 
\bigg(1-\frac{\delta}{\delta[-\beta {\varepsilon}_{a_m}]}\bigg)\;\mathcal{Y}_{\!\text{ref}}.\end{aligned}$$ This shows that the upper indices can be lowered iteratively, i.e., $\mathcal{C}_{i_1 \cdots i_n}^{a_1 \cdots a_m}=\mathcal{C}_{i_1 \cdots i_n}^{a_1 \cdots a_{m-1}}-\mathcal{C}_{i_1 \cdots i_n a_m}^{a_1 \cdots a_{m-1}}$, which leads to $$\begin{aligned}
 \label{chi1defpart3}
\mathcal{C}_{i_1 \cdots i_n}^{a_1 \cdots a_m}=
\sum_{\mathcal{P}\subset\{1,\ldots,m\}} (-1)^{|P|}
\mathcal{C}_{i_1 \cdots i_n \{a_k\}_{k \in P}  }.\end{aligned}$$ The cumulants $\mathcal{K}_{i_1 \ldots i_n}$ are defined by $$\begin{aligned}
\label{Maverage}
\mathcal{K}_{i_{1}\ldots i_{n}}= \frac{\delta^n \ln \mathcal{Y}_{\!\text{ref}}}{\delta[-\beta {\varepsilon}_{i_1}] \cdots \delta[-\beta {\varepsilon}_{i_n}]}.\end{aligned}$$ The relation between the $\mathcal{K}$’s and the $\mathcal{C}$’s is given by [@sokal] $$\begin{aligned}
  \label{generalG}
\mathcal{C}_{i_1 \cdots i_n}=\sum_{\substack{\mathcal{P}\in\,\text{partitions}\\\text{of}\,\{1,\ldots,n\}}} \prod_{I\in \mathcal{P}}
\mathcal{K}_{  \{i_k\}_{k\in I}  }.\end{aligned}$$ These formulas provide an alternative way (compared to the Wick contraction formalism) to evaluate the various contributions $\Upsilon_{[\Gamma_{1}^{\alpha_1}\cdots\Gamma_{k}^{\alpha_k}]_n}$ in Eq. .

### Simply connected unlinked diagrams

For linked diagrams without articulation lines (i.e., two-particle irreducible diagrams) the contributions from higher cumulants have measure zero for infinite systems. For such diagrams, the (sums of the) contributions from higher cumulants vanish also in the finite case, via exchange antisymmetry. This is clear, since these (nonextensive) contributions are absent in the Wick contraction formalism. For instance, for the first-order diagram $\mathcal{C}_{ij}=\mathcal{K}_{i}\mathcal{K}_{j}+\delta_{ij}\mathcal{K}_{ii}$ the part $\delta_{ij}\mathcal{K}_{ii}$ gives no contribution (by antisymmetry). Overall, it is $$\begin{aligned}
\label{highK1}
i_1\neq i_2\neq  \ldots \neq  i_n:\;\;\;
\mathcal{C}_{i_1\cdots i_n}=\prod_{\nu=1}^n \mathcal{K}_{i_\nu}\end{aligned}$$ for linked diagrams. This means that for linked two-particle irreducible diagrams the cumulant formalism leads to the same expressions as the Wick contraction formalism. For two-particle reducible diagrams, however, there are additional size extensive contributions from higher cumulants corresponding to articulation lines with identical three-momenta. This has the effect that for each set of normal articulation lines with identical three-momenta there is only a single distribution function, i.e., $$\begin{aligned}
 \label{Gnormal}
\mathcal{C}_{i_1 \cdots i_n j \cdots j}^{a_1 \cdots a_m}&=\mathcal{C}_{i_1 \cdots i_n j}^{a_1 \cdots a_m},
\\ \label{Gnormal2}
\mathcal{C}_{i_1 \cdots i_n}^{a_1 \cdots a_m  b \cdots b}&=\mathcal{C}_{i_1 \cdots i_n }^{a_1 \cdots a_m b},\end{aligned}$$ see Ref. [@Wellenhofer:2017qla]. For example, $\mathcal{C}_{i_1\ldots i_n jj}=\mathcal{C}_{i_1\ldots i_n}(\mathcal{K}_{j}\mathcal{K}_{j}+\mathcal{K}_{jj})=\mathcal{C}_{i_1\ldots i_n j}$, since $\mathcal{K}_{j}\mathcal{K}_{j}+\mathcal{K}_{jj}=n_jn_j+n_j\bar n_j=n_j$. Equations  and together with Eq.  imply that the contributions from anomalous diagrams are zero: $$\begin{aligned}
\label{anomzero}
\mathcal{C}_{i_1 \cdots i_n j \cdots j}^{a_1 \cdots a_m j \cdots j}=0.\end{aligned}$$ The contributions from anomalous diagrams and from the subleading parts of Eqs. –, etc. , now arise instead from unlinked diagrams (with normal subdiagrams). That is, for unlinked diagrams composed of $N$ subdiagrams the contributions from higher cumulants connecting $N$ lines with distinct three-momenta are size extensive. For example, for the case of three first-order subdiagrams with indices $(i,j)$, $(k,l)$, and $(m,n)$ one has the contributions $$\begin{aligned}
\delta_{ik}\delta_{jm}
\mathcal{K}_{i i}\mathcal{K}_{j j} \mathcal{K}_{l} \mathcal{K}_{n}, \ldots, \;\; 
\delta_{ik}\delta_{im}
\mathcal{K}_{i i i} \mathcal{K}_{j} \mathcal{K}_{l} \mathcal{K}_{n},\ldots,\end{aligned}$$ where $\mathcal{K}_{i i i}=n_i\bar n_i \bar n_i-n_i n_i \bar n_i$ and the ellipses represent terms with other index combinations. By virtue of the linked-cluster theorem, the size extensive contributions from unlinked diagrams where not all higher-cumulant indices correspond to different subdiagrams cancel against the corresponding terms with $\nu>1$ or $k_1>1$ in Eq. . The remaining size extensive contributions from unlinked diagrams are exactly those where the different (normal) subdiagrams are simply connected via higher cumulants. This provides a new representation for the contributions associated (in the Wick contraction formalism) with anomalous diagrams and the contributions not included in Eqs.  and .

### Modified thermodynamic perturbation series

There are two methods for the construction of the modified perturbation series for the free energy $F(T,\mu_\text{ref})$. The first, introduced by Kohn and Luttinger [@Kohn:1960zz], is based on grand-canonical MBPT; it constructs $F(T,\mu_\text{ref})$ in terms of a truncated formal expansion[^31] of $F(T,\mu)$ about $\mu_\text{ref}$, see Refs. [@Kohn:1960zz; @PhysRevC.89.064009; @Wellenhofer:2017qla] for details. The second method, due to Brout and Englert [@brout2], starts from the canonical ensemble. In canonical perturbation theory [@PhysRev.115.1374; @Parry], Eq.  is replaced by $$\begin{aligned}
\mathscr{C}_{i_1 \ldots i_n}=\braket{a_{i_1}^\dagger a_{i_1} \cdots a_{i_n}^\dagger a_{i_n}}_{\!\varrho},\end{aligned}$$ where $\braket{\ldots}_{\!\varrho}$ denotes the unperturbed canonical ensemble average which involves only Fock states $\ket{\Psi_{\!\varrho}}$ with fixed $\varrho=\Braket{\Psi_{\!\varrho}|\mathcal{N}|\Psi_{\!\varrho}}$. From this, we proceed analogously to the grand-canonical case, with $\mathcal{Y}_{\!\text{ref}}$ replaced by the unperturbed canonical partition function $\mathcal{Z}_\text{ref}=\sum_{\Psi_{\!\varrho}}\braket{\Psi_{\!\varrho}|\operatorname{\operatorname{e}}^{-\beta \mathcal{H}_\text{ref}}|\Psi_{\!\varrho}}$, i.e., the cumulants are now given by $$\begin{aligned}
\mathscr{K}_{i_{1}\ldots i_{n}}= \frac{\delta^n \ln \mathcal{Z}_\text{ref}}{\delta[-\beta {\varepsilon}_{i_1}] \cdots \delta[-\beta {\varepsilon}_{i_n}]}.\end{aligned}$$ The decisive new step is now to evaluate the cumulants not directly (which would be practically impossible) but using the Legendre transformation $$\begin{aligned}
 \label{legendre}
\ln \mathcal{Z}_\text{ref}(T,\varrho)&=\ln \mathcal{Y}_{\!\text{ref}}(T,\mu_\text{ref})- \mu_\text{ref} \frac{\partial \ln \mathcal{Y}_{\!\text{ref}}(T,\mu_\text{ref})}{\partial \mu_\text{ref}},\end{aligned}$$ where $\mu_\text{ref}$ is the chemical potential of an unperturbed grand-canonical system with the same mean fermion number as the fully interacting canonical system, i.e., $$\begin{aligned}
\varrho=-\frac{1}{\beta} \frac{\partial \ln \mathcal{Y}_{\!\text{ref}}(T,\mu_\text{ref})}{\partial \mu_\text{ref}}
=-\frac{\partial \Omega_\text{ref}(T,\mu_\text{ref})}{\partial \mu_\text{ref}}
=\sum_{\textbf}{k} \tilde{n}_{\textbf}{k},\end{aligned}$$ where $\tilde{n}_{\textbf}{k}$ denotes the Fermi-Dirac distribution with $\mu_\text{ref}$ as the chemical potential. With $\varrho$ being fixed, $\varrho=\sum_{\textbf}{k} \tilde{n}_{\textbf}{k}$ determines $\mu_\text{ref}$ as a functional of the spectrum $\varepsilon_{\textbf}{k}$. From this and Eq. , the expression for $\mathscr{K}_{i}$ is given by $$\begin{aligned}
\mathscr{K}_{i}
&=
\frac{\delta
\ln \mathcal{Y}_{\!\text{ref}}}{\delta[-\beta \varepsilon_i]}
-\frac{1}{\beta}\frac{\partial 
\ln \mathcal{Y}_{\!\text{ref}}}{\partial \mu_\text{ref}}
\Bigg(\frac{\delta \mu_\text{ref}}{\delta \varepsilon_i}\Bigg)_{\!\varrho}
-\varrho\,\Bigg( \frac{\delta\mu_\text{ref}}{\partial\varepsilon_i}\Bigg)_{\!\varrho}
=
\tilde{n}_{i}. \end{aligned}$$ The higher $\mathscr{K}$’s can then be determined iteratively, i.e., $$\begin{aligned}
\mathscr{K}_{i_1 i_2}
&=\Bigg(\frac{\delta \mathscr{K}_{i_1}}{\delta[-\beta \varepsilon_{i_2}]}\Bigg)_{\!\varrho}
{\nonumber \\}&=
\delta_{i_1i_2}
\frac{\partial \mathscr{K}_{i_1}}{\partial[-\beta \varepsilon_{i_1}]}
-
\frac{\partial \mathscr{K}_{i_1}}{\partial \mu_\text{ref}}\,
\Bigg[\frac{\partial \varrho}{\partial\mu_\text{ref}}\Bigg]^{-1}\!\!
\frac{\delta \varrho}{\delta[-\beta \varepsilon_{i_2}]}
{\nonumber \\}&=
\delta_{i_1i_2} \tilde{n}_{i_1} (1-\tilde{n}_{i_1})
-\frac{\tilde{n}_{i_1} (1-\tilde{n}_{i_1})\tilde{n}_{i_2} (1-\tilde{n}_{i_2})}{\sum_{i} \tilde{n}_{i} (1-\tilde{n}_{i})}.\end{aligned}$$ For $\mathscr{K}_{i_1 i_2 i_3}$ and beyond there is also a contribution where the energy derivative acts on $[\partial \varrho/\partial\mu_\text{ref}]^{-1}=\sum_{i} \tilde{n}_{i} (1-\tilde{n}_{i})$, i.e.,[^32] $$\begin{aligned}
 \label{Kscr3}
\big[\mathscr{K}_{i_1 i_2 i_3}\big]_{i_1\neq i_2 \neq i_3}
&=
-
\frac{\partial \mathscr{K}_{i_1i_2}}{\partial \mu_\text{ref}}
\Bigg[\frac{\delta \varrho}{\delta\mu_\text{ref}}\Bigg]^{-1}
\frac{\delta \varrho}{\delta[-\beta \varepsilon_{i_3}]}
{\nonumber \\}&\quad 
-
\frac{\partial \mathscr{K}_{i_1}}{\beta \mu_\text{ref}}
\frac{\delta \varrho}{\delta[-\beta \varepsilon_{i_2}]}
\frac{\partial}{\partial[-\beta\varepsilon_{i_3}]} \Bigg[\frac{\partial \varrho}{\partial\mu_\text{ref}}\Bigg]^{-1}.\end{aligned}$$ One can show that $\big[\mathscr{K}_{i_1 \ldots i_n}\big]_{i_a \neq i_b \,\forall a,b\in[1,n]}=\mathcal{O}(1/\varrho^{n-1})$, see Ref. [@Wellenhofer:2017qla], so the size extensive contributions from unlinked diagrams are again given by simply connected diagrams. For isotropic systems the anomalous parts of these contributions cancel at each order in the zero-temperature limit,[^33] thus $$\begin{aligned}
\text{isotropy:}\;\;
F(T,\mu_\text{ref}) \xrightarrow{T\rightarrow 0} E^{(0)}({\varepsilon}_{\text{F}}),\end{aligned}$$ with $\mu_\text{ref} \xrightarrow{T\rightarrow 0} {\varepsilon}_{\text{F}}$. By construction, within each of the order-by-order renormalization schemes (direct, cyclic, BdD), at each order the modified perturbation series $F(T,\mu_\text{ref})$ matches the grand-canonical perturbation series for the free energy $F(T,\mu)=\Omega(T,\mu)-\mu \,\partial \Omega(T,\mu)/\partial \mu$. The zero-temperature limit exists however only for the BdD scheme (see Sec. \[sec2\]).

Factorization theorem(s) {#sec43}
------------------------

Using the direct formula, the cyclic formula, or the reduced formula for finite systems ($\ast$) and applying the cumulant formalism to Eq.  leads to $$\begin{aligned}
 \label{redfac2cumulant}
\sum_{P/A}\Upsilon^{\text{direct},\div}_{[\widetilde{\Pi}_{1}^{\alpha_1}\cdots
\widetilde{\Pi}_{k}^{\alpha_k}]_n}
&= 
-\frac{1}{\beta} \prod_{i=1}^k
\frac{\left(-\beta \,\sum_{P_i/A_i}
\Upsilon^{\text{direct},\div}_{\widetilde{\Pi}_{i}}\right)^{\alpha_i}}{\alpha_i!},
\\
\label{redfac2cumulantcyc}
\sum_{P/A}\Upsilon^{\text{cyclic},\div}_{[\widetilde{\Pi}_{1}^{\alpha_1}\cdots
\widetilde{\Pi}_{k}^{\alpha_k}]_n}
&= 
-\frac{1}{\beta} \prod_{i=1}^k
\frac{\left(-\beta \,\sum_{P_i/A_i}
\Upsilon^{\text{cyclic},\div}_{\widetilde{\Pi}_{i}}\right)^{\alpha_i}}{\alpha_i!},\end{aligned}$$ $$\begin{aligned}
\label{redfac7}
\sum_{P/A} \Upsilon^{\text{reduced,}\ast,(\div)}_{[\widetilde{\Pi}_{1}^{\alpha_1}\cdots
\widetilde{\Pi}_{k}^{\alpha_k}]_n}
&= 
-\frac{1}{\beta} \prod_{i=1}^k
\frac{\left(-\beta \,\sum_{P_i/A_i}
\Upsilon^{\text{reduced,}\ast,\div}_{\widetilde{\Pi}_{i}}\right)^{\alpha_i}}{\alpha_i!},\end{aligned}$$ where the $\widetilde{\Pi}_i$ are all normal diagrams, and $P/A$ excludes those permutations that lead to anomalous diagrams. The combinatorics (and sign factors) of the higher-cumulant connections matches the combinatorics of the functional derivatives that generate the mean-field contributions from the perturbative contributions to the grand-canonical potential. Hence, Eqs.  and prove the direct and the cyclic version of the factorization property given by Eq.  and its cyclic analog, and Eq.  proves the reduced factorization property for finite systems $(\ast)$, Eq. . Note that Eq.  implies that in the reduced finite case the pseudoanomalous contributions vanish at each order.

The reduced version of the factorization theorem can also be proved as follows. For a given unlinked diagram where none of the linked parts are overlapping (see Fig. \[figx\]), the reduced formula has the form $$\begin{aligned}
\Upsilon^{\text{reduced,}\ast,(\div)}_{[{\Pi}_{1}^{\alpha_1}\cdots
{\Pi}_{k}^{\alpha_k}]_n} 
&\sim  
\underset{z=0}{\text{Res}}  \frac{\operatorname{\operatorname{e}}^{-\beta z}}{z} \frac{1}{(-z)^K}
\prod_i\frac{1}{{D}_i-z} 
{\nonumber \\}& \sim 
\beta^{K-1}
\prod_i\frac{1}{{D}_i} + \text{extra terms},\end{aligned}$$ where the extra terms are proportional to $\beta^{K-n}$, with $n\in\{2,\ldots,K\}$. The reduced expressions for unlinked diagrams with overlapping linked parts are composed entirely of such extra terms. These extra terms are incompatible with the linked-cluster theorem: they do not match the temperature dependence of (the disentangled reduced expressions) for the corresponding contributions with $\nu>1$ or $k_1>1$ in Eq. . The extra terms must therefore cancel each other at each order in the sum $\sum_{P/A}$.[^34] Thus, symbolically we have $$\begin{aligned}
 \label{redfac7b}
\sum_{P/A}
\Upsilon^{\text{reduced,}\ast,(\div)}_{[\widetilde{\Pi}_{1}^{\alpha_1}\cdots
\widetilde{\Pi}_{k}^{\alpha_k}]_n}
& \sim 
\beta^{K-1} \sum_{P/A}\left[
\prod_i\frac{1}{{D}_i} \right],\end{aligned}$$ which is equivalent to Eq. .

Statistical quasiparticles {#sec44}
--------------------------

The energy denominator regularization maintains the linked-cluster theorem. From the proof of the (reduced) factorization theorem it can be inferred that this suffices to establish that $$\begin{aligned}
 \label{redfac7BdD}
\sum_{P/A} \Upsilon^{\text{reduced,}\ast\ast,\div}_{[\widetilde{\Pi}_{1}^{\alpha_1}\cdots
\widetilde{\Pi}_{k}^{\alpha_k}]_n}
&= 
-\frac{1}{\beta} \prod_{i=1}^k
\frac{\left(-\beta \,\sum_{P_i/A_i}
\Upsilon^{\text{reduced,}\ast\ast,\div}_{\widetilde{\Pi}_{i}}\right)^{\alpha_i}}{\alpha_i!},\end{aligned}$$ which (by virtue of the cumulant formalism) implies the BdD factorization property $$\begin{aligned}
\label{reducedfactorizedBdD}
\Omega_{n_1+n_2,\text{anomalous}}^{\text{reduced,}\ast\ast,\div}
&= -\frac{\beta}{2} \sum_{\textbf}{k} 
U^{\text{reduced,}\ast\ast,\div}_{n_1,{\textbf}{k}} n_{\textbf}{k} \bar n_{\textbf}{k} \,
U^{\text{reduced,}\ast\ast,\div}_{n_2,{\textbf}{k}}
{\nonumber \\}& \quad \times 
(2-\delta_{n_1,n_2})
,\end{aligned}$$ and similar \[i.e., as specified by Eq. \] for anomalous contributions with several pieces (subdiagrams, in the cumulant formalism).

It is now clear how the cancellation between the contributions from simply connected diagrams composed of $V$ vertices and those where also $-U$ vertices are present works. For a given simply connected diagram, only the subdiagrams with *single* higher-cumulant connections can be replaced by $-U$ vertices, so at truncation orders $2N+1$ and $2N+2$ all anomalous contributions are removed if the mean field includes all contributions $U^{\aleph,\div}_{n,{\textbf}{k}}$ with $n\leq N$. However, this does *not* imply consistency with the adiabatic formalism for $U^{\aleph,\div}_{n,{\textbf}{k}}=U^{\text{reduced},\ast\ast,\div}_{n,{\textbf}{k}}$ (irrespective of isotropy), since the relation between chemical potential $\mu$ and the fermion number $\varrho$ does not match the adiabatic relation $\varrho=\sum_{\textbf}{k} \theta({\varepsilon}_{\text{F}}-{\varepsilon}_{\textbf}{k})$. For the consistency of the grand-canonical and the adiabatic formalism, the BdD mean field must include all contributions up to the truncation order; only then one preserves the thermodynamic relations of the pure mean-field theory (where ${H=H_0+U}$, with $U\equiv U[n_{\textbf}{k}]$), i.e., the Fermi-liquid relations $$\begin{aligned}
 \label{StatQP1aFERMILIQ}
\varrho &= 
\sum_{\textbf}{k} n_{\textbf}{k},
\\ \label{StatQP2aFERMILIQ}
S &=-\sum_{\textbf}{k} \big( n_{\textbf}{k}\ln n_{\textbf}{k} +  
\bar n_{\textbf}{k}\ln \bar n_{\textbf}{k}\big) , 
\\  \label{StatQP3aFERMILIQ}
\frac{\delta E}{\delta n_{\textbf}{k}} &=  {\varepsilon}_{{\textbf}{k}} .\end{aligned}$$ These relations are valid for all temperatures.

Zero-temperature limit {#sec45}
----------------------

At zero temperature, the energy denominator poles are at the boundary of the integration region,[^35] which implies that the contributions from two-particle reducible diagrams with several identical energy denominators diverge [@Wellenhofer:2018dwh; @Feldman1996]. For MBPT with ${U=0}$ or ${U=U_1}$ one finds that the divergent contributions cancel each other at each order.[^36] This cancellation is maintained in the BdD renormalization scheme: the cancellation occurs separately for normal contributions, and for the sum of the matching contributions the Sokhotski-Plemelj-Fox formula is consistent with the ${T\rightarrow 0}$ limit.

Notably, the energy denominator regularization is not required to construct a thermodynamic perturbation series that is consistent with the adiabatic formalism in the anisotropic case: at ${T=0}$, the BdD factorization theorem takes the form $$\begin{aligned}
\label{reducedfactorizedT0}
T=0:\;\;\Omega_{n_1+n_2,\text{anomalous}}
&= -\frac{1}{2} \sum_{\textbf}{k} 
U^{\text{reduced},\ast\ast,\div}_{n_1,{\textbf}{k}} 
\delta({\varepsilon}_{\textbf}{k}-\mu) 
{\nonumber \\}&\quad \times U^{\text{reduced},\ast\ast,\div}_{n_2,{\textbf}{k}}
(2-\delta_{n,m}),\end{aligned}$$ and similar for anomalous diagrams with several pieces. (At ${T=0}$, the symbols $\ast\ast$ and $\div$ (and the specification of $\aleph$ to reduced) are not needed for the separation of normal and anomalous contributions to the grand-canonical potential.) Thus, as recognized by Feldman *et al.* [@Feldman1999], to cancel the anomalous contributions to $\Omega(T=0,\mu)$ the following mean field is sufficient (for truncation orders below $2N+2$): $$\begin{aligned}
 \label{Ufermi}
U^{L_{\text{F}}}_{{\textbf}{k}} = U_{1,{\textbf}{k}} + \sum_{n=2}^N L_{\text{F}}\left[ U^{\text{reduced,}\ast\ast,\div}_{n,{\textbf}{k}}(T=0,\mu)\right],\end{aligned}$$ where $L_{\text{F}}$ satisfies $L_{\text{F}}[g({\textbf}{k})] = g({\textbf}{k})$ for ${\varepsilon}_{\textbf}{k}=\mu$ and is smoothed off away from ${\varepsilon}_{\textbf}{k}=\mu$. There are still anomalous contributions to the particle number, so the adiabatic series is not reproduced. The renormalization given by Eq.  (with $\mu$ replaced by $\mu_\text{ref}$) is however sufficient for the consistency of the adiabatic formalism with the modified perturbation series $F(T,\mu_\text{ref})$ in the anisotropic case.

Conclusion {#summary}
==========

In the present paper, we have, substantiating the outline by Balian and de Dominicis (BdD) [@statquasi3; @statquasi1],[^37] derived a thermodynamic perturbation series for infinite Fermi systems that (1) is consistent with the adiabatic zero-temperature formalism for both isotropic and anisotropic systems and (2) satisfies at each order and for all temperatures the thermodynamic relations associated with Fermi-liquid theory. This result arises, essentially, as a corollary of the linked-cluster theorem. The proof of (2) \[which implies (1)\] given here relies, apart from the earlier analysis of the disentanglement ($\div$) conducted by Balian, Bloch, and de Dominicis [@Balian1961529] and the outline provided by Balian and de Dominicis, on the application of the cumulant formalism (as a systematic method to perform $\div$) introduced by Brout and Englert [@brout2; @PhysRev.115.824]. The statistical quasiparticles associated with the thermodynamic Fermi-liquid relations are distinguished from the dynamical quasiparticles associated with the asymptotic stability of the low-lying excited states; in particular, the energies of dynamical and statistical quasiparticles are different.

In the perturbation series derived in the present paper the reference Hamiltonian is renormalized at each order in terms of additional contributions to the self-consistent mean-field potential. Conceptually, such an order-by-order renormalization is appealing: At each new order, not only is new information about interaction effects included, but this information automatically improves the reference point. Nevertheless, the relevance of this perturbation series depends on its convergence rate compared to the modified perturbation series for the free energy $F(T,\mu_\text{ref})$ with a fixed reference Hamiltonian; e.g., Hartree-Fock, or the (modified) second-order BdD mean field. In addition to the complete removal of anomalous contributions, the higher-order mean-field contributions lead also to partial cancellations of normal two-particle reducible contributions. This suggests that the convergence rate may indeed improve by renormalizing the mean field at each order.[^38] Apart from the question of convergence, beyond second order the practicality of the BdD renormalization scheme is impeded by the increasingly complicated regularization procedure required for its numerical application.

An alternative renormalization scheme, the direct scheme, was introduced by Balian, Bloch, and de Dominicis [@Balian1961529] (and rederived in the present paper, together with yet another scheme, the cyclic one). The thermodynamic relations resulting from the direct scheme however deviate from the Fermi-liquid relations. More severely, for the direct (and the cyclic) scheme the zero-temperature limit does not exist. The direct scheme may however still be useful for numerical calculations close to the classical limit. In particular, the corresponding perturbation series reproduces the virial expansion in the classical limit [@Balian1961529b]. The BdD renormalization scheme is thus mainly targeted at calculations not too far from the degenerate limit, in particular perturbative nuclear matter calculations (see, e.g., Refs. [@BOGNER200559; @PhysRevC.82.014314; @PhysRevC.83.031301; @PhysRevLett.110.032504; @PhysRevC.88.025802; @PhysRevC.89.025806; @PhysRevC.93.054314; @PhysRevC.94.054307; @Drischler:2017wtt; @PhysRevC.87.014322; @PhysRevC.89.044321; @PhysRevC.91.054311; @Sammarruca:2018bqh]. Notably, the statistical quasiparticle relations may be useful for the application of the Sommerfeld expansion [@prakash] and to connect with phenomenological parametrizations [@Yasin:2018ckc].

In conclusion, future research in the many-fermion problem will investigate the perturbation series derived in the present paper.[^39]\
I thank A. Carbone, C. Drischler, K. Hebeler, J. W. Holt, F. Hummel, N. Kaiser, R. Lang, M. Prakash, S. Reddy, A. Schwenk and W. Weise for useful discussions. Moreover, I thank the referees for helpful comments. Finally, I thank for their warm hospitality the group T39 (TU München), the INT (Seattle), the CEA (Saclay) and the ECT\* (Trento), where parts of this work have been presented. This work is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Projektnummer 279384907– SFB 1245 as well as the DFG and NSFC through the CRC 110 “Symmetries and the Emergence of Structure in QCD”.

Two-particle reducible diagrams at fourth order {#app1}
===============================================

Here we derive explicitly the regularized ($\ast\ast$) disentangled ($\div$) reduced expressions for the contributions from two-particle reducible diagrams at fourth order. Diagrams with single-vertex loops are canceled by the ones with $-U_1$ vertices; the remaining diagrams with $V$ vertices only are shown in Fig. \[fig4\]. One can choose indices such that for each diagram the matrix elements are given by $$\begin{aligned}
\zeta=
V^{ij,ab}V^{ik,cd}V^{cd,ki}V^{ab,ij},\end{aligned}$$ and the energy denominators corresponding to the two second-order pieces are given by $D_1=D_{ab,ij}$ and $D_2=D_{cd,ik}$. The cyclic expression for the sum of the diagrams in each row $\nu \in\{1,2,3,4\}$ can then be written as $$\begin{aligned}
\Omega^\text{cyclic}_{4,\nu}
=
\xi_{\nu} \sum_{ijkabcd}
\zeta \;
\mathcal{N}_{\nu} \, 
\mathcal{F}_{\nu}^\text{cyclic},\end{aligned}$$

where $\xi_{1,2,3,4}=(-1/4,-1/4,1/8,1/8)$. For the chosen indices the $\mathcal{N}_\nu$ are fixed as $\mathcal{N}_{1}=\mathcal{N}_{2}= n_{iijk} \bar{n}_{abcd}$,$\mathcal{N}_{3}= n_{iijcd} \bar{n}_{abk}$, and $\mathcal{N}_{4}= n_{abk} \bar{n}_{iijcd}$. Finally, from Eq. , the $\mathcal{F}^\text{cyclic}_\nu$ are given by $$\begin{aligned}
\mathcal{F}_{1}^\text{cyclic}
&=
\frac{1}{D_1^2 D_{1+2}}
-\frac{\operatorname{\operatorname{e}}^{-\beta (D_{1+2})}}{D_2^2 D_{1+2}}
+\frac{\operatorname{\operatorname{e}}^{-\beta D_1}D_{1-2}}{D_1^2 D_2^2}
-\beta\frac{\operatorname{\operatorname{e}}^{-\beta D_1}}{D_1 D_2},
\\
\mathcal{F}_{2}^\text{cyclic}
&=
\frac{1}{D_1D_2 D_{1+2}}
-\frac{\operatorname{\operatorname{e}}^{-\beta (D_{1+2})}}{D_1D_2 D_{1+2}}
+\frac{\operatorname{\operatorname{e}}^{-\beta D_1}}{D_1 D_2 D_{1-2}}
{\nonumber \\}&\quad
-\frac{\operatorname{\operatorname{e}}^{-\beta D_2}}{D_1 D_2 D_{1-2}},
\\
\mathcal{F}_{3}^\text{cyclic}
&=
\frac{1}{D_1^2 D_{1-2}}
-\frac{\operatorname{\operatorname{e}}^{-\beta (D_{1-2})}}{D_2^2 D_{1-2}}
+\frac{\operatorname{\operatorname{e}}^{-\beta D_1}D_{1+2}}{D_1^2 D_2^2}
+\beta\frac{\operatorname{\operatorname{e}}^{-\beta D_1}}{D_1 D_2},
\\
\mathcal{F}_{4}^\text{cyclic}
&=
-\frac{1}{D_1^2 D_{1-2}}
+\frac{\operatorname{\operatorname{e}}^{\beta (D_{1-2})}}{D_2^2 D_{1-2}}
-\frac{\operatorname{\operatorname{e}}^{\beta D_1}D_{1+2}}{D_1^2 D_2^2}
+\beta\frac{\operatorname{\operatorname{e}}^{\beta D_1}}{D_1 D_2},\end{aligned}$$ where $D_{1\pm 2}=D_1 \pm D_2$. Although their individual parts have poles, the $\mathcal{F}^\text{cyclic}_{\nu}$ are regular for any zero of $D_1 D_2 D_{1+2} D_{1-2}$. To separate the various parts, we add infinitesimal imaginary parts to the energy denominators, i.e., $$\begin{aligned}
{D}_{1} \rightarrow \,& {D}_{1,\eta_{1}}={D}_{1} +{\text{i}}\eta_{1},
\\
{D}_{2} \rightarrow \,& {D}_{2,\eta_{2}}={D}_{2} +{\text{i}}\eta_{2},\end{aligned}$$ where $|\eta_1| \neq |\eta_2|$, since otherwise ${D}_{1+2,\eta_{1+2}}$ or ${D}_{1-2,\eta_{1-2}}$ has zeros. Averaging over the signs of the imaginary parts and applying Eqs. \[ndouble\] and , we can reorganize the sum of the 12 diagrams according to $$\begin{aligned}
\sum_{\nu=1}^4\Omega^\text{cyclic}_{4,\nu}
=
\Omega_{4,\text{normal}}^{\text{reduced,}\ast\ast,\div} + 
\Omega_{4,\text{anom.}}^{\text{reduced,}\ast\ast,\div} + 
\Omega_{4,\text{pseudo-a.}}^{\text{reduced,}\ast\ast,\div},\end{aligned}$$

where $$\begin{aligned}
 \label{4normal}
\Omega_{4,\text{normal}}^{\text{reduced,}\ast\ast,\div} = & 
\sum_{\alpha=1}^4 
\bigg[
\frac{1}{8} \sum_{ijkabcd} \!\!
\zeta \; 
\mathcal{N}^\text{normal}_{\alpha}\;
\mathcal{F}^{\text{reduced},\ast\ast}_{\alpha,\text{normal}} 
\bigg]
,
\\ \label{4anom}
\Omega_{4,\text{anom.}}^{\text{reduced,}\ast\ast,\div} = & 
\sum_{\alpha=1,3,4} 
\bigg[
\frac{\beta}{8} \sum_{ijkabcd} \!\!
\zeta \; 
\mathcal{N}^\text{anom.}_{\alpha} \;
\mathcal{F}^{\text{reduced},\ast\ast}_{\alpha,\text{anom.}} 
\bigg]
, 
\\ \label{4pseud}
\Omega_{4,\text{pseudo-a.}}^{\text{reduced,}\ast\ast,\div} = & 
\sum_{\alpha=1,3,4}  
\bigg[
\frac{1}{8} \sum_{ijkabcd} \!\!
\zeta \; 
\mathcal{N}^\text{anom.}_{\alpha}\;
\mathcal{R}^{\text{reduced},\ast\ast}_{\alpha}
\bigg]
,\end{aligned}$$ with $$\begin{aligned}
\mathcal{F}^{\text{reduced},\ast\ast}_{\alpha,\text{normal}}
=
\!\!\!\!\!  \sum_{\text{sgn}(\eta_{1}),\text{sgn}(\eta_{2})} \!\!
\mathcal{F}^{\text{reduced},\ast\ast}_{\alpha,\text{normal},[\eta_1,\eta_2]}
,\end{aligned}$$ and similar for the anomalous and pseudoanomalous contributions. The correspondence $\alpha \cong \nu$ holds only for the anomalous contributions, and the normal ones with $\alpha=3,4$. For the normal contributions with $\alpha=1,2$, we combine the (disentangled) contributions from the first two ($\alpha=1$) and the third two ($\alpha=2$) diagrams in the first two rows. Regarding the pseudoanomalous contributions, each $\alpha$ corresponds to several $\nu$’s, by virtue of the application of Eq. . In the anomalous contribution $$\begin{aligned}
 \label{Nanom1}
\mathcal{N}^\text{anom.}_1&= n_{iabk} \bar{n}_{ijcd},\\ \label{Nanom2}
\mathcal{N}^\text{anom.}_3&= n_{iabcd} \bar{n}_{ijk},\\ \label{Nanom3}
\mathcal{N}^\text{anom.}_4&= n_{ijk} \bar{n}_{iabcd},\end{aligned}$$

![The 12 fourth-order two-particle reducible diagrams composed of two second-order normal pieces. In each of the four rows (1, 2, 3, 4), the first (and third) diagram is a normal diagram, the other anomalous. The diagrams in each row transform into each other under cyclic vertex permutations. The set of all 12 diagrams is closed under general vertex permutations. []{data-label="fig4"}](fig6.pdf){width="48.00000%"}

and $$\begin{aligned}
\mathcal{F}^{\text{reduced},\ast\ast}_{1,\text{anom.}}
&=\!\!\!\!\!\!
\sum_{\text{sgn}(\eta_{1}),\text{sgn}(\eta_{2})} \!\!
\left[-\frac{2}{{D}_{1,\eta_1}{D}_{2,\eta_2}}\right]
,
\\
\mathcal{F}^{\text{reduced},\ast\ast}_{3,\text{anom.}}
&=\!\!\!\!\!\!
\sum_{\text{sgn}(\eta_{1}),\text{sgn}(\eta_{2})} \!\!\frac{1}{{D}_{1,\eta_1}{D}_{2,\eta_2}},
\\
\mathcal{F}^{\text{reduced},\ast\ast}_{4,\text{anom.}}
&=\!\!\!\!\!\!
\sum_{\text{sgn}(\eta_{1}),\text{sgn}(\eta_{2})} \!\!\frac{1}{{D}_{1,\eta_1}{D}_{2,\eta_2}}.\end{aligned}$$ Suitably relabeling indices, we obtain the BdD factorization property $$\begin{aligned}
 \label{4thorderfactorized}
\Omega_{4,\text{anom.}}^{\text{reduced,}\ast\ast,\div}
=-\frac{\beta}{2}\sum_i  U_{2,i}^{\text{reduced,}\ast\ast,(\div)}   n_i \bar n_i \, U_{2,i}^{\text{reduced,}\ast\ast,(\div)},\end{aligned}$$ with $U_{2,i}^{\text{reduced,}\ast\ast,(\div)}$ given by Eq. . Relabeling indices according to Eqs. , , and , the energy denominators in the pseudoanomalous contribution are $$\begin{aligned}
\mathcal{R}^{\text{reduced},\ast\ast}_{1,[\eta_1,\eta_2]}
&=
-\frac{2}{{D}_{1,\eta_1}({D}_{2,\eta_2})^2} 
+\frac{2}{({D}_{1,\eta_1})^2{D}_{2,\eta_2}}
{\nonumber \\}& \quad
-\frac{2}{{D}_{1,\eta_1}{D}_{2,\eta_2}{D}_{1-2,\eta_{1-2}}}
-\frac{2}{{D}_{1,\eta_2}{D}_{2,\eta_1}{D}_{1-2,\eta_{2-1}}}
{\nonumber \\}& \quad
+\frac{1}{({D}_{2,\eta_1})^2{D}_{1-2,\eta_{2-1}}}
+\frac{1}{({D}_{2,\eta_2})^2{D}_{1-2,\eta_{1-2}}}
{\nonumber \\}& \quad
+\frac{1}{({D}_{1,\eta_1})^2{D}_{1-2,\eta_{1-2}}}
+\frac{1}{({D}_{1,\eta_2})^2{D}_{1-2,\eta_{2-1}}}
,
\\
\mathcal{R}^{\text{reduced},\ast\ast}_{3,[\eta_1,\eta_2]}
&=
-\frac{2}{({D}_{2,\eta_2})^2{D}_{1+2,\eta_{1+2}}}
-\frac{2}{{D}_{1,\eta_1}{D}_{2,\eta_2}{D}_{1+2,\eta_{1+2}}}
{\nonumber \\}& \quad
+\frac{1}{{D}_{1,\eta_1}({D}_{2,\eta_2})^2} 
+\frac{1}{({D}_{1,\eta_1})^2{D}_{2,\eta_2}}
,
\\
\mathcal{R}^{\text{reduced},\ast\ast}_{4,[\eta_1,\eta_2]}
&=
\frac{2}{({D}_{2,\eta_2})^2{D}_{1+2,\eta_{1+2}}}
+\frac{2}{{D}_{1,\eta_1}{D}_{2,\eta_2}{D}_{1+2,\eta_{1+2}}}
{\nonumber \\}& \quad
-\frac{1}{{D}_{1,\eta_1}({D}_{2,\eta_2})^2} 
-\frac{1}{({D}_{1,\eta_1})^2{D}_{2,\eta_2}}.\end{aligned}$$ In these expressions, the parts with three different denominators require special attention: the formal application of the Sokhotski-Plemelj-Fox formula assumes that each energy denominator is used as an explicit integration variable, but this is not possible for terms with denominators of the form ${D}_{1,\eta_1}{D}_{2,\eta_2}{D}_{1\pm 2,\eta_{1\pm 2}}$. To evaluate these terms, we use the relations $$\begin{aligned}
 \label{combineD}
\frac{1}{({D}_{1,\eta_1})^2{D}_{1\pm 2,\eta_{1\pm 2}}}\pm
\frac{1}{{D}_{1,\eta_1}{D}_{2,\eta_2}{D}_{1\pm 2,\eta_{1\pm 2}}}
&=
\pm \frac{1}{({D}_{1,\eta_1})^2{D}_{2,\eta_{2}}},
\\ \label{combineD2}
\frac{1}{({D}_{2,\eta_2})^2{D}_{1\pm 2,\eta_{1\pm 2}}}\pm
\frac{1}{{D}_{1,\eta_1}{D}_{2,\eta_2}{D}_{1\pm 2,\eta_{1\pm 2}}}
&=
\pm\frac{1}{{D}_{1,\eta_1}({D}_{2,\eta_{2}})^2}.\end{aligned}$$ This leads to

$$\begin{aligned}
\mathcal{R}^{\text{reduced},\ast\ast}_{1}
&=0
,
\\ \label{4pseudoav2}
\mathcal{R}^{\text{reduced},\ast\ast}_{3}
&=\!\!\!
\sum_{\text{sgn}(\eta_{1}),\text{sgn}(\eta_{2})} 
\left[
-\frac{1}{{D}_{1,\eta_1}({D}_{2,\eta_2})^2} 
+\frac{1}{({D}_{1,\eta_1})^2{D}_{2,\eta_2}}
\right]
,
\\ \label{4pseudoav3}
\mathcal{R}^{\text{reduced},\ast\ast}_{4}
&=\!\!\!
\sum_{\text{sgn}(\eta_{1}),\text{sgn}(\eta_{2})}
\left[
\frac{1}{{D}_{1,\eta_1}({D}_{2,\eta_2})^2} 
-\frac{1}{({D}_{1,\eta_1})^2{D}_{2,\eta_2}}
\right]
.\end{aligned}$$

One sees that $\mathcal{R}^{\text{reduced},\ast\ast}_{3}$ and $\mathcal{R}^{\text{reduced},\ast\ast}_{4}$ are antisymmetric under $D_1 \leftrightarrow D_2$. In each case, the remaining part of the integrand is symmetric under $D_1 \leftrightarrow D_2$. Thus, the pseudoanomalous contribution is zero: $$\begin{aligned}
\Omega_{4,\text{pseudo-a.}}^{\text{reduced,}\ast\ast,\div}
=0.\end{aligned}$$ Finally, in the normal contribution $$\begin{aligned}
 \label{Nnorm1}
\mathcal{N}^\text{normal}_1&= n_{ijk} \bar{n}_{abcd},\\ \label{Nnorm2}
\mathcal{N}^\text{normal}_2&= n_{abcd} \bar{n}_{ijk},\\ \label{Nnorm3}
\mathcal{N}^\text{normal}_3&= n_{ijcd} \bar{n}_{abk},\\ \label{Nnorm4}
\mathcal{N}^\text{normal}_4&= n_{abk} \bar{n}_{ijcd},\end{aligned}$$ and

$$\begin{aligned}
 \label{Enorm1}
\mathcal{F}^{\text{reduced},\ast\ast}_{1,\text{normal},[\eta_1,\eta_2]}
&=
-\frac{2}{({D}_{1,\eta_1})^2{D}_{1+2,\eta_{1+2}}}
{\nonumber \\}& \quad
-\frac{2}{{D}_{1,\eta_1}{D}_{2,\eta_2}{D}_{1+2,\eta_{1+2}}}
,
\\ \label{Enorm2}
\mathcal{F}^{\text{reduced},\ast\ast}_{2,\text{normal},[\eta_1,\eta_2]}
&=
\frac{2}{({D}_{2,\eta_2})^2{D}_{1+2,\eta_{1+2}}}
+\frac{2}{{D}_{1,\eta_1}{D}_{2,\eta_2}{D}_{1+2,\eta_{1+2}}}
,
\\ \label{Enorm3}
\mathcal{F}^{\text{reduced},\ast\ast}_{3,\text{normal},[\eta_1,\eta_2]}
&=
\frac{1}{({D}_{1,\eta_1})^2{D}_{1-2,\eta_{1-2}}}
+\frac{1}{({D}_{1,\eta_2})^2{D}_{1-2,\eta_{2-1}}}
,
\\ \label{Enorm4}
\mathcal{F}^{\text{reduced},\ast\ast}_{4,\text{normal},[\eta_1,\eta_2]}
&=-\frac{1}{({D}_{1,\eta_1})^2{D}_{1-2,\eta_{1-2}}}
-\frac{1}{({D}_{1,\eta_2})^2{D}_{1-2,\eta_{2-1}}}
,\end{aligned}$$

where we have suitably relabeled indices. Applying Eqs.  and , the averaged expressions are given by $$\begin{aligned}
 \label{Enorm1av}
\mathcal{F}^{\text{reduced},\ast\ast}_{1,\text{normal}}
&=\!\!\!
\sum_{\text{sgn}(\eta_{1}),\text{sgn}(\eta_{2})}
\left[-\frac{2}{({D}_{1,\eta_1})^2{D}_{2,\eta_{2}}}\right]
,
\\ \label{Enorm2av}
\mathcal{F}^{\text{reduced},\ast\ast}_{2,\text{normal}}
&=\!\!\!
\sum_{\text{sgn}(\eta_{1}),\text{sgn}(\eta_{2})}
\frac{2}{{D}_{1,\eta_{1}}({D}_{2,\eta_2})^2}
,
\\ \label{Enorm3av}
\mathcal{F}^{\text{reduced},\ast\ast}_{3,\text{normal}}
&=\!\!\!
\sum_{\text{sgn}(\eta_{1}),\text{sgn}(\eta_{2})}
\frac{1}{({D}_{1,\eta_1})^2{D}_{1-2,\eta_{1-2}}}
,
\\ \label{Enorm4av}
\mathcal{F}^{\text{reduced},\ast\ast}_{4,\text{normal}}
&=\!\!\!
\sum_{\text{sgn}(\eta_{1}),\text{sgn}(\eta_{2})} 
\left[-\frac{1}{({D}_{1,\eta_1})^2{D}_{1-2,\eta_{1-2}}}\right]
.\end{aligned}$$ In addition to the contribution from the 12 diagrams shown in Fig. \[fig4\], in the BdD scheme the two-particle reducible contribution at fourth order involves the six diagrams of Fig. \[fig3red\] with the first-order subdiagrams replaced by $-U_{2}^\text{BdD}$ vertices, and the diagram composed of two $-U_{2}^\text{BdD}$ vertices (Fig. \[fig4UU\]). The anomalous contributions from these 19 diagrams cancel each other (as a consequence of Eq. ). Notably, there is also a partial analytic cancellation between the contributions from Eqs.  and and the normal contribution from the third-order diagrams with one $-U_{2}^\text{BdD}$ vertex (Fig. \[fig3red\]), see Refs. [@Becker:1971asg; @Jones:1970hwh; @Wellenhofer:2017qla]. Such partial cancellations can be found also at higher orders for the normal contribution from certain (normal) two-particle reducible diagrams, i.e., for those where cutting the articulation lines and closing them such that an unlinked diagram (with two linked parts) is generated leaves the number of holes invariant.[^40]

Finally, the two-particle reducible fourth-order contribution to $U^\text{BdD}$ is given by the functional derivative of the regularized disentangled reduced normal contributions from the diagrams of Fig. \[fig4\] and the ones of Fig. \[fig3red\] with the first-order subdiagrams replaced by $-U_{2}^\text{BdD}$ vertices.[^41]

Self-energy, mass function, mean field and all that {#app2}
===================================================

Here, we discuss the various forms of the self-energy, and their relation to the grand-canonical potential, the mean occupation numbers, and the (various forms of the) mean field.

Analytic continuation(s) of the Matsubara self-energy {#app20}
-----------------------------------------------------

Although it is defined in terms of the self-consistent Dyson equation, the proper Matsubara self-energy $\Xi_{\textbf}{k}(z_l)$ can also be calculated using bare propagators; in that case, also two-particle reducible self-energy diagrams contribute to $\Xi_{\textbf}{k}(z_l)$; see, e.g., Ref. [@PLATTER2003250]. From $\Xi_{\textbf}{k}(z_l)$ the frequency-space self-energy $\Sigma_{\textbf}{k}(z)$ is obtained as the analytic continuation of $\Xi_{\textbf}{k}(z_l)$ that has the following properties:[^42]

1.  $\Sigma_{\textbf}{k}(z)$ is analytic in both the upper and lower half plane, vanishes at infinity, and has a branch cut along the real axis where $\text{Im}[\Sigma_{\textbf}{k}(z)]$ changes sign, with ${\text{Im}[\Sigma_{\textbf}{k}(z)]\lessgtr 0}$ for ${\text{Im}[z] \gtrless 0}$.

With these properties, $\Sigma_{\textbf}{k}(z)$ leads to the spectral representation of the mean occupation number, see Eq.  below.

Now, as shown below, in bare MBPT another analytic continuation of the Matsubara self-energy $\Xi_{\textbf}{k}(z_l)$ can be defined, here referred to as the mass function $\mathcal{M}_{\textbf}{k}(z)$. It has the following properties:

1.  $\mathcal{M}_{\textbf}{k}(z)$ is entire and real on the real axis for ${T\neq 0}$.

2.  It vanishes at infinity, except for $1/\text{Re}[z]=0^{+}$ where it has an essential singularity.

3.  For $\text{Re}[z]>\mu$ it has an essential singularity at ${T=0}$.

If $z={\varepsilon}_{\textbf}{k}+\mathcal{M}_{\textbf}{k}(z)$ has no solutions off the real axis, then one can obtain from $\mathcal{M}_{\textbf}{k}(z)$ another simple expression for the mean occupation numbers: the mass function representation $f_{\textbf}{k}=n(\mathscr{E}_{\textbf}{k})$, see Eq. . The ${T\rightarrow 0}$ limit of this representation is singular for $\mathscr{E}_{\textbf}{k}>\mu$, and gives $f_{\textbf}{k}=n_{\textbf}{k}$ for $\mathscr{E}_{\textbf}{k}<\mu$. The issue whether $z={\varepsilon}_{\textbf}{k}+\mathcal{M}_{\textbf}{k}(z)$ may have nonreal solutions is discussed further below. We did not see an argument that guarantees that nonreal solutions exist.

The functional forms of the bare perturbative contributions to $\Xi_{\textbf}{k}(z_l)$, $\Sigma_{\textbf}{k}(z)$, and $\mathcal{M}_{\textbf}{k}(z)$ are related to ones of the different time-independent formulas ($\aleph$) for the perturbative contributions to the grand-canonical potential, $\Omega_n^\aleph$. This is examined in part \[app23\] of this Appendix.

Mean occupation numbers from Dyson equation {#app21}
-------------------------------------------

Here, we first derive the mass function representation for the mean occupation number, Eq. . Then, we derive the spectral representation Eq. . Only the spectral representation can be derived also from the real-time propagator.[^43] Last, we examine the relation between the collision self-energy $\Sigma^{\text{coll}}_{\textbf}{k}(\omega)$ and the frequency-space self-energy $\Sigma_{\textbf}{k}(z)$ at ${T=0}$.

### Mass function

The imaginary-time propagator is given by $$\begin{aligned}
\mathscr{G}_{\textbf}{k}(\tau-\tau')=-\Braket{\! \Braket{ \mathcal{T}\left[a_{\textbf}{k}(\tau) a_{\textbf}{k}^\dagger(\tau') \right] } \!},\end{aligned}$$ where $a_{\textbf}{k}(\tau)=a_{\textbf}{k}\operatorname{\operatorname{e}}^{-{\varepsilon}_{\textbf}{k} \tau}$ and $a^\dagger_{\textbf}{k}(\tau)=a^\dagger_{\textbf}{k}\operatorname{\operatorname{e}}^{{\varepsilon}_{\textbf}{k} \tau}$, and $\braket{\!\braket{\ldots}\!}$ is the true ensemble average. Its Fourier series is $$\begin{aligned}
 \label{Gfourierser}
\mathscr{G}_{\textbf}{k}(\tau)=\frac{1}{\beta} 
\lim_{l_\text{max} \rightarrow \infty}
\sum_{l\in \mathcal{L}(l_\text{max}) } \mathscr{G}_{\textbf}{k}(z_l) \operatorname{\operatorname{e}}^{-z_l\tau},\end{aligned}$$ where $\mathcal{L}(l_\text{max})=\{-l_\text{max},\ldots,l_\text{max}\}$, and $z_l$ are the Matsubara frequencies (see Eq. ). The Fourier coefficients are given by $$\begin{aligned}
\label{Gfourierserinv}
\mathscr{G}_{\textbf}{k}(z_l)=\int \limits_{0}^\beta d\tau  \,\mathscr{G}_{\textbf}{k}(\tau)\operatorname{\operatorname{e}}^{z_l\tau}.\end{aligned}$$ The Dyson equation in Fourier (Matsubara) space is given by $$\begin{aligned}
\mathscr{G}_{\textbf}{k}(z_l) 
&=
g_{\textbf}{k}(z_l) + g_{\textbf}{k}(z_l)\, \Xi_{\textbf}{k}(z_l)\, \mathscr{G}_{\textbf}{k}(z_l),\end{aligned}$$ where $\Xi_{\textbf}{k}(z_l)$ is the Matsubara self-energy and $g_{{\textbf}{k}}(\omega_l)$ is the unperturbed propagator in Matsubara space, i.e., $$\begin{aligned}
 \label{g0fourier}
g_{{\textbf}{k}}(z_l)
=\frac{1}{z_l-{\varepsilon}_{\textbf}{k}}.\end{aligned}$$ Iterating the Dyson equation and summing the resulting geometric series leads to $$\begin{aligned}
 \label{Dysonsummed}
\mathscr{G}_{\textbf}{k}(z_l) 
&=
\frac{1}{z_l-{\varepsilon}_{\textbf}{k}-\Xi_{\textbf}{k}(z_l)}.\end{aligned}$$ Inserting this into the Fourier series Eq.  and replacing the discrete frequency sums by a contour integral leads to $$\begin{aligned}
 \label{Gim0C0}
\mathscr{G}_{\textbf}{k}(\tau)
&= \!
\oint\limits_{C_0[l_\text{max}]} \!\!\!
\frac{d z}{2\pi {\text{i}}} \, 
\operatorname{\operatorname{e}}^{-z\tau}n(z)\,
\frac{1}{z-{\varepsilon}_{\textbf}{k}-\mathcal{M}_{\textbf}{k}(z)}  ,\end{aligned}$$ where $l_\text{max}\rightarrow \infty$ is implied, and $n(z)=[1+\operatorname{\operatorname{e}}^{\beta(z-\mu}]^{-1}$. The contour $C_0[l_\text{max}]$ encloses all the Matsubara poles $z=z_{l\in\mathcal{L}(l_\text{max})}$ but not the pole at $z={\varepsilon}_{\textbf}{k}+\mathcal{M}_{\textbf}{k}(z)$, see Fig. \[figc1\]. Note that by construction $C_0[l_\text{max}]$ crosses the real axis. Thus, for Eq.  to be equivalent to Eq. , the mass function $\mathcal{M}_{\textbf}{k}(z)$ must be an analytic continuation of the Matsubara self-energy $\Xi_{\textbf}{k}(z_l)$ that is analytic on the real axis (and near the Matsubara poles). This is easy to get: for the second-order two-particle irreducible contribution to $\mathcal{M}_{\textbf}{k}(z)$ we obtain from Eq.  the expression $$\begin{aligned}
 \label{sigma2bMass}
\mathcal{M}_{2,{\textbf}{k}}(z) &= 
\frac{1}{2}
\sum_{{\textbf}{k}_2, {\textbf}{k}_3, {\textbf}{k}_4} \!
|\braket{\psi_{{\textbf}{k}}\psi_{{\textbf}{k}_2}| V|\psi_{{\textbf}{k}_3} \psi_{{\textbf}{k}_4}}|^2  n_{{\textbf}{k}_2} \bar n_{{\textbf}{k}_3} \bar n_{{\textbf}{k}_4}
{\nonumber \\}& \quad \times
\frac{e^{-\beta({\varepsilon}_{{\textbf}{k}_3}+{\varepsilon}_{{\textbf}{k}_4} - {\varepsilon}_{{\textbf}{k}_2} - z)}-1}{{\varepsilon}_{{\textbf}{k}_3}+{\varepsilon}_{{\textbf}{k}_4} - {\varepsilon}_{{\textbf}{k}_2} - z},\end{aligned}$$ i.e., in contrast to Eq. , we do *not* substitute$\operatorname{\operatorname{e}}^{\beta(z_l-\mu)}=-1$ before performing the analytic continuation, and similar for higher-order contributions. Since with this prescription there are no poles, $\mathcal{M}_{\textbf}{k}(z)$ is entire and real on the real axis, and regular everywhere except for $\text{Re}[z]\rightarrow \infty$. However, the ${T\rightarrow 0}$ limit of $\mathcal{M}_{\textbf}{k}(z)$ is singular for $\text{Re}[z]>\mu$, due to terms $\operatorname{\operatorname{e}}^{\beta (z-\mu)}$ as in Eq. .

![Contours $C_0[l_\text{max}]$ (left panel) and $C_1$ (right panel).[]{data-label="figc1"}](C1.pdf){width="45.00000%"}

In Eq. , for $\tau<0$ the term $\operatorname{\operatorname{e}}^{-z\tau} n_{\textbf}{k}(z)$ is regular at infinity, and the term $[z-{\varepsilon}_{\textbf}{k}-\mathcal{M}_{\textbf}{k}(z)]^{-1}$ vanishes at infinity. Hence, if we assume that $z={\varepsilon}_{\textbf}{k}+\mathcal{M}_{\textbf}{k}(z)$ has no solutions off the real axis, for $\tau<0$ the contour $C_0[l_\text{max}]$ can be deformed into the contour $C_1$ (see Fig. \[figc1\]) that encloses the pole on the real axis at $$\begin{aligned}
 \label{massSC}
\mathscr{E}_{\textbf}{k}={\varepsilon}_{\textbf}{k}+\mathcal{M}_{\textbf}{k}(\mathscr{E}_{\textbf}{k}).\end{aligned}$$ For nonreal $z=x+iy$, $z={\varepsilon}_{\textbf}{k}+\mathcal{M}_{\textbf}{k}(z)$ is equivalent to the two coupled nonlinear equations $x={\varepsilon}_{\textbf}{k}+\text{Re}[\mathcal{M}_{\textbf}{k}(x+iy)]$ and $y=\text{Im}[\mathcal{M}_{\textbf}{k}(x+iy)]$. At second order, this is given by $$\begin{aligned}
 \label{Mcoup1}
x &= {\varepsilon}_{\textbf}{k} - U_{{\textbf}{k}} + U_{1,{\textbf}{k}} + \mathcal{M}^{\ddagger}_{2,{\textbf}{k}}
{\nonumber \\}& \quad  +
\frac{1}{2}
\sum_{{\textbf}{k}_2, {\textbf}{k}_3, {\textbf}{k}_4} \!
|\braket{\psi_{{\textbf}{k}}\psi_{{\textbf}{k}_2}| V|\psi_{{\textbf}{k}_3} \psi_{{\textbf}{k}_4}}|^2 
\frac{n_{{\textbf}{k}_2} \bar n_{{\textbf}{k}_3} \bar n_{{\textbf}{k}_4}}{\left[D(x)\right]^2+y^2}
{\nonumber \\}& \quad \times
\{ D(x)[\cos(\beta y)e^{-\beta D(x)}-1] -y \sin(\beta y) e^{-\beta D(x)} \}
,
\\ \label{Mcoup2}
y &= 
\frac{1}{2}
\sum_{{\textbf}{k}_2, {\textbf}{k}_3, {\textbf}{k}_4} \!
|\braket{\psi_{{\textbf}{k}}\psi_{{\textbf}{k}_2}| V|\psi_{{\textbf}{k}_3} \psi_{{\textbf}{k}_4}}|^2 
\frac{n_{{\textbf}{k}_2} \bar n_{{\textbf}{k}_3} \bar n_{{\textbf}{k}_4}}{\left[D(x)\right]^2+y^2}
{\nonumber \\}& \quad \times
\{ y[ \cos(\beta y)e^{-\beta D(x)}-1] +D(x) \sin(\beta y) e^{-\beta D(x)} \},\end{aligned}$$ with $D(x)={\varepsilon}_{{\textbf}{k}_3}+{\varepsilon}_{{\textbf}{k}_4} - {\varepsilon}_{{\textbf}{k}_2} - x$. Here, $-U_{{\textbf}{k}}$ corresponds to the diagram composed of a single $-U$ vertex, $U_{1,{\textbf}{k}}$ to the one with a single $V$ vertex, and $\mathcal{M}^{\ddagger}_{2,{\textbf}{k}}$ denotes the second-order two-particle reducible contribution to $\mathcal{M}_{2,{\textbf}{k}}(z)$. We did not see an argument that guarantees that Eqs.  and have solutions for ${y\neq 0}$.

Assuming that there is only the pole given by Eq. , we get $$\begin{aligned}
 \label{fMass0}
\mathscr{G}_{\textbf}{k}(\tau<0)
&=n\left(\mathscr{E}_{\textbf}{k}\right)\operatorname{\operatorname{e}}^{-\mathscr{E}_{\textbf}{k}\tau},\end{aligned}$$ and the expression for the mean occupation numbers is given by $$\begin{aligned}
 \label{fMass}
f_{\textbf}{k}
&=
\mathscr{G}_{\textbf}{k}(0^{-})
=n\left(\mathscr{E}_{\textbf}{k}\right),\end{aligned}$$ i.e., the exact mean occupation numbers are given by the Fermi-Dirac distribution with the reference spectrum renormalized in terms of the on-shell mass function $\mathcal{M}_{\textbf}{k}(\mathscr{E}_{\textbf}{k})$ defined via Eq.  and the analytic continuation of the Matsubara self-energy $\Xi_{\textbf}{k}(z_l)$ that is real analytic on the real axis. As discussed, the ${T\rightarrow 0}$ limit of $\mathcal{M}_{\textbf}{k}(\mathscr{E}_{\textbf}{k})$ is ill-behaved.

### Frequency-space self-energy

The real-time propagator is[^44] $$\begin{aligned}
{\text{i}}G_{\textbf}{k}(t-t')= \Braket{\! \Braket{ \mathcal{T}\left[a_{\textbf}{k}(t) a_{\textbf}{k}^\dagger(t') \right] } \!},\end{aligned}$$ with $a_{\textbf}{k}(t)=a_{\textbf}{k}\operatorname{\operatorname{e}}^{-{\text{i}}{\varepsilon}_{\textbf}{k} t}$ and $a^\dagger_{\textbf}{k}(t)=a^\dagger_{\textbf}{k}\operatorname{\operatorname{e}}^{{\text{i}}{\varepsilon}_{\textbf}{k}t}$. It can be decomposed as $$\begin{aligned}
{\text{i}}G_{\textbf}{k}(t-t')= 
\theta(t-t') \underbrace{\braket{\! \braket{a_{\textbf}{k}(t) a_{\textbf}{k}^\dagger(t')  }\!} }_{ {\text{i}}G^{>}_{\textbf}{k}(t-t')}
-
\theta(t'-t)\underbrace{\braket{\!  \braket{a_{\textbf}{k}^\dagger(t') a_{\textbf}{k}(t)  }\!} }_{-{\text{i}}G^{<}_{\textbf}{k}(t-t')},\end{aligned}$$ where we have defined the correlation functions ${\text{i}}G^{>}_{\textbf}{k}(t-t')$ and ${\text{i}}G^{<}_{\textbf}{k}(t-t')$. The Fourier transforms of the real-time propagator and the correlation functions are given by $$\begin{aligned}
G_{\textbf}{k}(\omega)&=\int \limits_{-\infty}^\infty \!\! dt \, {\text{i}}G_{\textbf}{k}(t)\operatorname{\operatorname{e}}^{{\text{i}}\omega t},
\\
G^{>}_{\textbf}{k}(\omega)&=\int \limits_{-\infty}^\infty \!\! dt \, {\text{i}}G_{\textbf}{k}^{>}(t)\operatorname{\operatorname{e}}^{{\text{i}}\omega t},
\\
G^{<}_{\textbf}{k}(\omega)&=\int \limits_{-\infty}^\infty \!\! dt \, \left(-{\text{i}}G_{\textbf}{k}^{<}(t)\operatorname{\operatorname{e}}^{{\text{i}}\omega t}\right),\end{aligned}$$ with inverse transforms $$\begin{aligned}
 \label{Gfouriertransform}
{\text{i}}G_{\textbf}{k}(t)=\int \limits_{-\infty}^\infty \!\! \frac{d \omega}{2 \pi} \, G_{\textbf}{k}(\omega)  \operatorname{\operatorname{e}}^{-{\text{i}}\omega t},
\\
{\text{i}}G^{>}_{\textbf}{k}(t)=\int \limits_{-\infty}^\infty \!\! \frac{d \omega}{2 \pi} \, G^{>}_{\textbf}{k}(\omega)  \operatorname{\operatorname{e}}^{-{\text{i}}\omega t},
\\
-{\text{i}}G^{<}_{\textbf}{k}(t)=\int \limits_{-\infty}^\infty \!\! \frac{d \omega}{2 \pi} \, G^{<}_{\textbf}{k}(\omega)  \operatorname{\operatorname{e}}^{-{\text{i}}\omega t}.\end{aligned}$$ The Fourier transforms of the correlation functions satisfy the KMS relation [@kadanoffbaym; @ThesisRios] (see also Refs. [@PhysRev.115.1342; @Haag1967; @doi:10.1143/JPSJ.12.570]) $$\begin{aligned}
G^{<}_{\textbf}{k}(\omega) = \operatorname{\operatorname{e}}^{-\beta(\omega-\mu)} G^{>}_{\textbf}{k}(\omega).\end{aligned}$$ From this relation it follows that we can write $$\begin{aligned}
G^{>}_{\textbf}{k}(\omega)
&=
\bar n(\omega) \mathcal{A}_{\textbf}{k}(\omega) ,
\\
G^{<}_{\textbf}{k}(\omega)
&=
n(\omega) \mathcal{A}_{\textbf}{k}(\omega),\end{aligned}$$ with $$\begin{aligned}
\mathcal{A}_{\textbf}{k}(\omega) 
&=
G^{>}_{\textbf}{k}(\omega)+G^{<}_{\textbf}{k}(\omega),\end{aligned}$$ and $\bar n(\omega)=1-n(\omega)$. From the Lehmann representations of $G^{>}_{\textbf}{k}(t)$ and $G^{<}_{\textbf}{k}(t)$ it can be seen that the spectral function is semipositive, i.e., $$\begin{aligned}
 \label{Ageq0}
\mathcal{A}_{\textbf}{k}(\omega) & \geq 0 ,\end{aligned}$$ and satisfies the sum rule $$\begin{aligned}
 \label{sumrule}
\int \limits_{-\infty}^\infty \!\! \frac{d \omega}{2\pi}\mathcal{A}_{\textbf}{k}(\omega) & =1,\end{aligned}$$ see e.g., Refs. [@Fetter; @ThesisRios]. Consider now the function $\Gamma_{\textbf}{k}(z)$ defined by $$\begin{aligned}
 \label{Gammadef}
\Gamma_{\textbf}{k}(z)
&=
\int \limits_{-\infty}^\infty \!\! \frac{d \omega}{2\pi}
\frac{\mathcal{A}_{\textbf}{k}(\omega)}{z-\omega}.\end{aligned}$$ From the Lehmann representation of the imaginary-time propagator $\mathscr{G}_{\textbf}{k}(\tau)$ it can be seen that [@Fetter] $$\begin{aligned}
\mathscr{G}_{\textbf}{k}(z_l) 
&=
\int \limits_{-\infty}^\infty \!\! \frac{d \omega}{2\pi}
\frac{\mathcal{A}_{\textbf}{k}(\omega)}{z_l-\omega}.\end{aligned}$$ From the sum rule for $\mathcal{A}_{\textbf}{k}(\omega)$, Eq. , it then follows that $\Gamma_{\textbf}{k}(z)$ corresponds to the (unique [@baym]) analytic continuation of $\mathscr{G}_{\textbf}{k}(z_l)$ that satisfies $\Gamma_{\textbf}{k}(z)\sim z^{-1}$ for $|z|\rightarrow \pm\infty$. From Eq. , this can be obtained via[^45] $$\begin{aligned}
\label{Sigmadef}
\Gamma_{\textbf}{k}(z)
&=
\frac{1}{z-{\varepsilon}_{\textbf}{k}-\Sigma_{\textbf}{k}(z)},\end{aligned}$$ where the frequency-space self-energy $\Sigma_{\textbf}{k}(z)$ is defined as the analytic continuation of the Matsubara self-energy $\Xi_{\textbf}{k}(z_l)$ that satisfies $\Sigma_{\textbf}{k}(z)\rightarrow 0$ for $|z|\rightarrow \pm\infty$. In bare MBPT, this is given by the prescription noted before Eq.  in Sec. \[sec24\], i.e., one first substitutes $\operatorname{\operatorname{e}}^{\beta(\omega_l-\mu)}=-1$ and then performs the analytic continuation. For convenience, we give again the irreducible part of the bare second-order contribution to $\Sigma_{\textbf}{k}(z)$, i.e., $$\begin{aligned}
 \label{sigma2again}
\Sigma_{2,{\textbf}{k}}(z) &= 
-\frac{1}{2}
\sum_{{\textbf}{k}_2, {\textbf}{k}_3, {\textbf}{k}_4} \!
|\braket{\psi_{{\textbf}{k}}\psi_{{\textbf}{k}_2}| V|\psi_{{\textbf}{k}_3} \psi_{{\textbf}{k}_4}}|^2  
\frac{
n_{{\textbf}{k}_2} \bar n_{{\textbf}{k}_3} \bar n_{{\textbf}{k}_4} +
n_{{\textbf}{k}_3} n_{{\textbf}{k}_4} \bar n_{{\textbf}{k}_2}}
{{\varepsilon}_{{\textbf}{k}_3}+{\varepsilon}_{{\textbf}{k}_4} - {\varepsilon}_{{\textbf}{k}_2} - z}.\end{aligned}$$ From Eqs.  and , we obtain for the spectral function the expression $$\begin{aligned}
 \label{spectralfundef}
\mathcal{A}_{\textbf}{k}(\omega)
&={\text{i}}\Big[\Gamma_{\textbf}{k}(\omega+{\text{i}}\eta) - \Gamma_{\textbf}{k}(\omega-{\text{i}}\eta) \Big].
{\nonumber \\}&={\text{i}}\left[\frac{1}{\omega-{\varepsilon}_{\textbf}{k}-\Sigma_{\textbf}{k}(\omega+{\text{i}}\eta)+{\text{i}}\eta} - \text{c.c.} \right],\end{aligned}$$ where c.c. denotes the complex conjugate. Note that inserting Eq.  into Eq.  leads to the Breit-Wigner form of the spectral function, Eq. . The relation between $\mathcal{A}_{\textbf}{k}(\omega)$ and the Fourier transform of the real-time propagator $G_{\textbf}{k}(\omega)$ is obtained as follows: $$\begin{aligned}
 \label{Gomega}
G_{\textbf}{k}(\omega) 
&=
\int \limits_{-\infty}^\infty \!\! dt \,\operatorname{\operatorname{e}}^{{\text{i}}\omega t} \Big[\theta(t) {\text{i}}G^{>}_{\textbf}{k}(t)+\theta(-t) {\text{i}}G^{<}_{\textbf}{k}(t) \Big]
{\nonumber \\}&=
-\int \limits_{-\infty}^\infty \!\! dt \,\operatorname{\operatorname{e}}^{{\text{i}}\omega t} 
\left[
\int \limits_{-\infty}^\infty \!\! \frac{d \xi}{2 \pi{\text{i}}} \,\frac{\operatorname{\operatorname{e}}^{-{\text{i}}\xi t} }{\xi+{\text{i}}\eta} {\text{i}}G^{>}_{\textbf}{k}(t)
\right.
 \left.
-
\int \limits_{-\infty}^\infty \!\! \frac{d \xi}{2 \pi{\text{i}}} \,\frac{\operatorname{\operatorname{e}}^{-{\text{i}}\xi t} }{\xi-{\text{i}}\eta}  {\text{i}}G^{<}_{\textbf}{k}(t)
\right]
{\nonumber \\}&= 
-\int \limits_{-\infty}^\infty \!\! \frac{d \xi}{2\pi {\text{i}}} 
\left[ 
\frac{G^{>}_{\textbf}{k}(\omega-\xi)}{\xi+{\text{i}}\eta}+\frac{G^{<}_{\textbf}{k}(\omega-\xi)}{\xi-{\text{i}}\eta}
\right]
{\nonumber \\}&= 
-\int \limits_{-\infty}^\infty \!\! \frac{d \xi}{2\pi{\text{i}}} 
\left[ 
\frac{\bar n(\xi) \mathcal{A}_{\textbf}{k}(\xi)}{\omega-\xi+{\text{i}}\eta} +\frac{n(\xi) \mathcal{A}_{\textbf}{k}(\xi)}{\omega-\xi-{\text{i}}\eta} 
\right],\end{aligned}$$ where we have used the relation [@Fetter] $$\begin{aligned}
\theta(\pm t) = \mp\int\limits_{-\infty}^\infty \!\!  \frac{d \xi}{2 \pi{\text{i}}} \, \frac{\operatorname{\operatorname{e}}^{-{\text{i}}\xi t} }{\xi\pm {\text{i}}\eta}.\end{aligned}$$ From Eq.  we then have $$\begin{aligned}
{\text{i}}G_{\textbf}{k}(t)= 
-\int \limits_{-\infty}^\infty \!\!  
\frac{d \omega}{2\pi}
\operatorname{\operatorname{e}}^{-{\text{i}}\omega t}\!\!
\int \limits_{-\infty}^\infty \!\! \frac{d \xi}{2\pi{\text{i}}} 
\left[ 
\frac{\bar n(\xi) \mathcal{A}_{\textbf}{k}(\xi)}{\omega-\xi+{\text{i}}\eta} +\frac{n(\xi) \mathcal{A}_{\textbf}{k}(\xi)}{\omega-\xi-{\text{i}}\eta} 
\right].\end{aligned}$$ For $t<0$ we can close the $\omega$ integral in the upper half plane. Interchanging the integration order, we then get $$\begin{aligned}
 \label{Gtless0}
{\text{i}}G_{\textbf}{k}(t<0)=
-\int \limits_{-\infty}^\infty  \!\! \frac{d \xi}{2\pi} \,
 \operatorname{\operatorname{e}}^{-{\text{i}}\xi t} n(\xi)  \mathcal{A}_{\textbf}{k}(\xi).\end{aligned}$$ Thus, the expression for the exact mean occupation numbers is $$\begin{aligned}
 \label{fspectralrep}
f_{\textbf}{k}
&=
-{\text{i}}G_{\textbf}{k}(0^{-})
=\int \limits_{-\infty}^\infty \!\! \frac{d \xi}{2\pi} \,
n(\xi)  \mathcal{A}_{\textbf}{k}(\xi).\end{aligned}$$ Here, in contrast to Eq. , the ${T\rightarrow 0}$ limit is well-behaved, and its analysis reveals that $f_{\textbf}{k}(T=0,\mu)$ has a discontinuity at ${\textbf}{k}={\textbf}{k}_{\text{F}}$, see Ref. [@PhysRev.119.1153] and Sec. \[sec24\].

![Contours $C^{\pm}_0$ (left panel) and $C_2$ (right panel).[]{data-label="figc2"}](C2.pdf){width="45.00000%"}

The result given by Eq.  can also be obtained directly from the Fourier expansion of the imaginary-time propagator, Eq. . That is, taking first the limit ${l_\text{max}\rightarrow\infty}$ and then performing the analytic continuation of $\Xi_{\textbf}{k}(z_l)$ to $\Sigma_{\textbf}{k}(z)$ we get $$\begin{aligned}
\mathscr{G}_{\textbf}{k}(\tau)
&=
\oint\limits_{C^{\pm}_0} \frac{d z}{2\pi {\text{i}}} \, 
\operatorname{\operatorname{e}}^{-z\tau} n(z) 
\frac{1}{z-{\varepsilon}_{\textbf}{k}-\Sigma_{\textbf}{k}(z)} ,\end{aligned}$$ with $C^{\pm}_0=C^{+}_0+C^{-}_0$, where $C^{+}_0$ encloses the Matsubara poles in the upper half plane without crossing the real axis, and $C^{-}_0$ the poles in the lower half plane. Since $\Sigma_{\textbf}{k}(z)$ is analytic in the two half planes and vanishes at complex infinity, and $\operatorname{\operatorname{e}}^{-z\tau} n_{\textbf}{k}(z)$ is regular at infinity for $\tau<0$, for $\tau<0$ these two contours can be deformed into the contour $C_2$ that encloses the real axis, see Fig. \[figc2\], i.e., $$\begin{aligned}
 \label{Gim0C0again}
\mathscr{G}_{\textbf}{k}(\tau<0)
&= \!\!
\int\limits_{-\infty}^\infty \!\! \frac{d \xi}{2\pi {\text{i}}}  \,
\operatorname{\operatorname{e}}^{-\xi\tau} n(\xi) 
\underbrace{\left[\frac{1}{\xi-{\varepsilon}_{\textbf}{k}-\Sigma_{\textbf}{k}(\xi+{\text{i}}\eta)+{\text{i}}\eta}
-\text{c.c.}\right]}_{-{\text{i}}\mathcal{A}_{\textbf}{k}(\xi)},\end{aligned}$$ which is just the Wick rotation of Eq. .

### Collision self-energy at zero temperature

The self-energy corresponding to the real-time propagator, here referred to as the collision self-energy $\Sigma_{\textbf}{k}^\text{coll}(\omega)$, can be defined by [@Fetter][^46] $$\begin{aligned}
 \label{Ecoll}
G_{\textbf}{k}(\omega) = {\text{i}}\frac{1}{\omega-{\varepsilon}_{\textbf}{k}- \Sigma_{\textbf}{k}^\text{coll}(\omega) }.\end{aligned}$$ In the following, we examine how at ${T=0}$ the collision self-energy relates to the frequency-space self-energy $\Sigma_{\textbf}{k}(z)$. For this, using the Sokhotski-Plemelj theorem we rewrite Eq.  as $$\begin{aligned}
G_{\textbf}{k}(\omega) = 
\int \limits_{-\infty}^\infty \!\! \frac{d \xi}{2\pi{\text{i}}}  \mathcal{A}_{\textbf}{k}(\xi) \frac{P}{\omega-\xi} 
-
\frac{\bar n(\omega) \mathcal{A}_{\textbf}{k}(\omega)}{2}
+
\frac{ n(\omega) \mathcal{A}_{\textbf}{k}(\omega)}{2}
.\end{aligned}$$ From $n(\omega)\xrightarrow{T\rightarrow 0}\theta(\mu-\omega)$, at zero temperature we have $$\begin{aligned}
G_{\textbf}{k}(\omega) = 
\theta(\omega-\mu)\,G_{\textbf}{k}^\text{R}(\omega)
+
\theta(\mu-\omega)\,G_{\textbf}{k}^\text{A}(\omega)
,\end{aligned}$$ with the Fourier transforms of the retarded and advanced propagators given by $$\begin{aligned}
G_{\textbf}{k}^R(\omega)&= 
\int \limits_{-\infty}^\infty \!\! \frac{d \xi}{2\pi{\text{i}}}  \mathcal{A}_{\textbf}{k}(\xi) \frac{P}{\omega-\xi} 
-
\frac{\mathcal{A}_{\textbf}{k}(\omega)}{2},
\\
G_{\textbf}{k}^A(\omega)&= 
\int \limits_{-\infty}^\infty \!\! \frac{d \xi}{2\pi{\text{i}}}  \mathcal{A}_{\textbf}{k}(\xi) \frac{P}{\omega-\xi} 
+
\frac{\mathcal{A}_{\textbf}{k}(\omega)}{2}.\end{aligned}$$ Comparing with Eq.  we see that $$\begin{aligned}
G_{\textbf}{k}^\text{R}(\omega) = 
{\text{i}}\Gamma_{\textbf}{k}(\omega+ {\text{i}}\eta),
\;\;\;\;\;\;\;\;
G_{\textbf}{k}^\text{A}(\omega) = 
{\text{i}}\Gamma_{\textbf}{k}(\omega- {\text{i}}\eta).\end{aligned}$$ Thus, from Eq.  we can at ${T=0}$ make the identification $$\begin{aligned}
\Sigma_{\textbf}{k}^\text{coll}(\omega)=
\theta(\omega-\mu)\,\Sigma_{\textbf}{k}(\omega+ {\text{i}}\eta)
+
\theta(\mu-\omega)\,\Sigma_{\textbf}{k}(\omega- {\text{i}}\eta).\end{aligned}$$ From Eq. , i.e., $$\begin{aligned}
 \label{dynquasi1aagain}
\Sigma_{\textbf}{k}(\omega\pm{\text{i}}\eta) = \mathcal{S}_{\textbf}{k}(\omega) \mp {\text{i}}\mathcal{J}_{\textbf}{k}(\omega),\end{aligned}$$ we have (at ${T=0}$) $$\begin{aligned}
\Sigma_{\textbf}{k}^\text{coll}(\omega) &= 
\theta(\omega-\mu)\,\big[\mathcal{S}_{\textbf}{k}(\omega)-{\text{i}}\mathcal{J}_{\textbf}{k}(\omega) \big]
{\nonumber \\}& \quad 
+\theta(\mu-\omega)\,\big[\mathcal{S}_{\textbf}{k}(\omega)+{\text{i}}\mathcal{J}_{\textbf}{k}(\omega) \big].\end{aligned}$$ In particular, (as discussed in Sec. \[sec24\]), at zero temperature it is ${\mathcal{J}_{\textbf}{k}(\omega) \xrightarrow{\omega \rightarrow \mu} C_{\textbf}{k}(\mu) (\omega-\mu)^2}$, with $C_{\textbf}{k}(\mu)\geq 0$, so (at ${T=0}$) $$\begin{aligned}
\text{Im}\big[\Sigma_{\textbf}{k}^\text{coll}(\omega)\big] &= 
-\theta(\omega-\mu)\,\mathcal{J}_{\textbf}{k}(\omega)
+\theta(\mu-\omega)\,\mathcal{J}_{\textbf}{k}(\omega) 
{\nonumber \\}&\xrightarrow{\omega \rightarrow \mu}
-C_{\textbf}{k}(\mu)\, (\omega-\mu)|\omega-\mu|.\end{aligned}$$ Finally, for the on-shell collision self-energy this leads to $$\begin{aligned}
\text{Im}\big[\Sigma_{\textbf}{k}^\text{coll}({\varepsilon}_{\textbf}{k})\big] 
&\xrightarrow{{\varepsilon}_{\textbf}{k} \rightarrow \mu}
-C_{\textbf}{k}(\mu)\, ({\varepsilon}_{\textbf}{k}-\mu)|{\varepsilon}_{\textbf}{k}-\mu|,\end{aligned}$$ which is the property quoted in Refs.[@Holt13a; @Kaiser:2001ra].[^47]

Mean occupation numbers from direct mean-field renormalization {#app22}
--------------------------------------------------------------

In the direct renormalization scheme the exact mean occupation numbers are identified with the Fermi-Dirac distributions (i.e., with the mean occupation numbers in the reference system), i.e., $$\begin{aligned}
 \label{fisn}
\text{direct scheme:}\;\;\;\;
f_{\textbf}{k}=n_{\textbf}{k}.\end{aligned}$$ From this one may conclude that in the direct scheme the mass function is zero, $\mathcal{M}_{\textbf}{k}(z)=0$, and the spectral function is given by the unperturbed one, $\mathcal{A}_{\textbf}{k}(\omega)=2\pi \delta(\omega-{\varepsilon}_{\textbf}{k})$. More generally, one may conclude that the Matsubara self-energy is zero, $\Xi_{\textbf}{k}(z_l)=0$.

However, these conclusions come with two caveats:

1.  The cancellations that lead to Eq.  are not available in Matsubara space. That is, the result $\Xi_{\textbf}{k}(z_l)=0$ is obtained only from the Fourier expansion of the direct expression for the propagator, Eq. . If one instead Fourier expands the (unperturbed) propagators (cf. Appendix \[app23\]) in the time-integral representation, Eq. , then one obtains the usual result, i.e., $\Xi_{\textbf}{k}(z_l)\neq 0$, also in the direct scheme.

2.  The (proper) Matsubara self-energy $\Xi_{\textbf}{k}(z_l)$ is defined in terms of the Dyson equation, Eq. . The Dyson equation is inconsistent with a perturbative truncation order. In contrast, Eq.  relies on a finite truncation order $N$.

Caveat (1) implies that the ${T\neq 0}$ part[^48] of the general results of Sec. \[sec24\] can be obtained also in the direct scheme, and caveat (2) makes evident that these results and Eq.  do not contradict each other; they correspond to different partial summations of a divergent asymptotic series.

### Proof of Eq. (\[fisn\])

The perturbation series for the imaginary-time propagator $\mathscr{G}_{{\textbf}{k}}(\tau-\tau')$ is given by $$\begin{aligned}
\mathscr{G}_{n,{\textbf}{k}}(\tau-\tau')
&=
g_{{\textbf}{k}}(\tau-\tau')
+\sum_{n=1}^N \mathscr{G}_{n,{\textbf}{k}}(\tau-\tau').\end{aligned}$$ Here, the unperturbed propagator $g_{{\textbf}{k}}(\tau-\tau')$ is given by $$\begin{aligned}
 \label{smallg}
g_{{\textbf}{k}}(\tau-\tau')
&= - \Braket{ \mathcal{T}\left[a_{\textbf}{k}(\tau) a_{\textbf}{k}^\dagger(\tau') \right] }
{\nonumber \\}&= \theta(\tau-\tau')\, n_{\textbf}{k} \operatorname{\operatorname{e}}^{{\varepsilon}_{\textbf}{k} (\tau-\tau')} 
-
\theta(\tau'-\tau)\, \bar n_{\textbf}{k} \operatorname{\operatorname{e}}^{{\varepsilon}_{\textbf}{k} (\tau-\tau')},\end{aligned}$$ and its Matsubara coefficients are given by Eq. . The perturbative contributions $\mathscr{G}_{n,{\textbf}{k}}(\tau-\tau')$ are given by the expression [@Fetter] $$\begin{aligned}
 \label{Gmbpt}
\mathscr{G}_{n,{\textbf}{k}}(\tau-\tau')
&=
\frac{(-1)^{n+1}}{n!}
\int \limits_{0}^{\beta} d\tau_n\cdots d\tau_1 
{\nonumber \\}&\quad \times
\Braket{  \mathcal{T}\big[ a_{\textbf}{k}(\tau) \,
\mathcal{V}(\tau_n) \cdots
\mathcal{V}(\tau_1) \,
a_{\textbf}{k}^\dagger(\tau') \big]
}_{L}
{\nonumber \\}&\equiv\mathscr{G}^{\text{direct}(P)}_{n,{\textbf}{k}}(\tau-\tau').\end{aligned}$$ This can be written as $$\begin{aligned}
\mathscr{G}_{n,{\textbf}{k}}(\tau-\tau')
&=
(-1)^{n+1}
\int \limits_{0}^{\beta} d\tau_{n} \int \limits_{0}^{\tau_n} d\tau_{n-1}\cdots \int \limits_{0}^{\tau_2} d\tau_1 
{\nonumber \\}&\quad \times
\Braket{  \mathcal{T}\big[ a_{\textbf}{k}(\tau) \,
\mathcal{V}(\tau_{n}) \cdots
\mathcal{V}(\tau_1) \,
a_{\textbf}{k}^\dagger(\tau') \big]
}_{\!L}
{\nonumber \\}&\equiv\mathscr{G}^\text{direct}_{n,{\textbf}{k}}(\tau-\tau')
.\end{aligned}$$ From here, we can follows the steps that lead Bloch and de Dominicis [@BDDnuclphys7] to the direct formula[^49] for the perturbative contributions to the grand-canonical potential (see Sec. \[sec22\]). Because $\mathscr{G}_{n,{\textbf}{k}}(\tau)$ is antiperiodic with period $\beta$ we can, without loss of generality, set $\tau<0$ and $\tau'<0$. For $\tau-\tau'<0$, this leads to $$\begin{aligned}
 \label{Gdirect}
\mathscr{G}^\text{direct}_{n,{\textbf}{k}}(\tau<0)&=
\operatorname{\operatorname{e}}^{{\varepsilon}_{\textbf}{k}\tau}
\frac{(-1)^{n}}{2\pi {\text{i}}} 
\oint_{C}  dz \frac{\operatorname{\operatorname{e}}^{-\beta z}}{z^2}
{\nonumber \\}&\quad \times
\Braket{ 
\mathcal{V} \frac{1}{{D}_n-z} \cdots
\mathcal{V} \frac{1}{{D}_1-z}
\mathcal{V} \,a_{\textbf}{k}^\dagger\, a_{\textbf}{k}
}_{\!L}.\end{aligned}$$ For truncation order $N$, the contributions to $\mathscr{G}_{{\textbf}{k}}(\tau)$ are given by all linked (one-particle irreducible and reducible) propagator diagrams that satisfy Eq. . Applying the cumulant formalism, the contributions to Eq.  are given by normal propagator diagrams[^50] with normal Hugenholtz diagrams attached via higher-cumulant connections, plus diagrams composed of multiple normal propagator diagrams simply-connected via higher-cumulant connections attached to normal Hugenholtz diagrams. With the mean field given by $$\begin{aligned}
U_{\textbf}{k}=\sum_{n=1}^N U_{n,{\textbf}{k}}^{\text{direct},\div},\end{aligned}$$ the contributions with higher-cumulant connections are removed. *Furthermore*, because propagator diagrams involve all possible orderings of the vertices, (an analog of) the direct factorization theorem applies also for the remaining contributions; e.g., for a one-particle irreducible propagator diagram with non $-U$ self-energy part (i.e., at least one $V$ vertex is involved) we have $$\begin{aligned}
 \label{Gdirectfactorized1}
\mathscr{G}_{n,{\textbf}{k}}^{\div}(\tau)&= g_{{\textbf}{k}}(\tau) \;
U^{\text{direct,}\div}_{n,{\textbf}{k}},\end{aligned}$$ for a one-particle reducible diagram with two non $-U$ self-energy parts we have $$\begin{aligned}
 \label{Gdirectfactorized2}
\mathscr{G}_{n_1+n_2,{\textbf}{k}}^{\div}(\tau)&= g_{{\textbf}{k}}(\tau) \;
U^{\text{direct,}\div}_{n_1,{\textbf}{k}} \;  U^{\text{direct,}\div}_{n_2,{\textbf}{k}},\end{aligned}$$ etc. Hence, in the direct scheme these contributions are canceled by the diagrams where the self-energy parts are replaced by $-U^{\text{direct},\div}$ vertices. Thus,[^51] $$\begin{aligned}
\text{direct scheme:}\;\;\;\;
\mathscr{G}_{{\textbf}{k}}(\tau)=
g_{{\textbf}{k}}(\tau),\end{aligned}$$ and Eq. (\[fisn\]) is proved.

Self-energy, mass function, and grand-canonical potential {#app23}
---------------------------------------------------------

The (proper) Matsubara self-energy $\Xi_{\textbf}{k}(z_l)$ can be calculated using self-consistent propagators or using bare propagators (or, anything in between). In the bare case, also two-particle reducible self-energy diagrams contribute to $\Xi_{\textbf}{k}(z_l)$; see, e.g., Ref. [@PLATTER2003250].

Below, we first explain how the bare perturbative contributions to the improper Matsubara self-energy $\Xi^\star_{\textbf}{k}(z_l)$ can be obtained. From this, the bare contributions to $\Xi_{\textbf}{k}(z_l)$ are obtained via the restriction to one-particle irreducible diagrams.

Second, we derive the functional relations between the bare perturbative contributions to the (various forms of the) improper self-energy and the grand-canonical potential.[^52] In particular, we find the simple relation for the proper frequency-space self-energy $\Sigma_{\textbf}{k}(z)$ given by Eq. .

### Matsubara self-energy

The improper Matsubara self-energy $\Xi^\star_{{\textbf}{k}}(z_l)$ is defined by [@Luttinger:1960ua] $$\begin{aligned}
 \label{Xistar0}
\mathscr{G}_{{\textbf}{k}}(z_l)
&=
g_{{\textbf}{k}}(z_l)
+
g_{{\textbf}{k}}(z_l)\,
\Xi^\star_{{\textbf}{k}}(z_l)\,
g_{{\textbf}{k}}(z_l),\end{aligned}$$ i.e., the perturbative contributions to $\Xi^\star_{{\textbf}{k}}(z_l)$ are defined by $$\begin{aligned}
 \label{Xistar1}
\mathscr{G}_{n,{\textbf}{k}}(z_l)
&=
g_{{\textbf}{k}}(z_l)\,
\Xi^\star_{n,{\textbf}{k}}(z_l)\,
g_{{\textbf}{k}}(z_l).\end{aligned}$$ For example, from Eq.  the second-order irreducible contribution to $\mathscr{G}_{{\textbf}{k}}(\tau)$ is given by $$\begin{aligned}
\mathscr{G}_{2,{\textbf}{k}}(\tau)
&=
-\frac{1}{2} \sum_{{\textbf}{k}_2,{\textbf}{k}_3,{\textbf}{k}_4}
|\braket{\psi_{{\textbf}{k}}\psi_{{\textbf}{k}_2}| V|\psi_{{\textbf}{k}_3}\psi_{{\textbf}{k}_4}}|^2 
\int \limits_{0}^{\beta} \! d\tau_1 \int \limits_{0}^{\beta} \!d\tau_2
{\nonumber \\}&\quad \times
g_{\textbf}{k}(\tau-\tau_1)g_{\textbf}{k}(\tau_2-0)
g_{{\textbf}{k}_2}(\tau_{21})g_{{\textbf}{k}_3}(\tau_{12}) g_{{\textbf}{k}_5}(\tau_{12}),\end{aligned}$$ with $\tau_{ij}=\tau_i-\tau_j$. Inserting the Fourier series of the unperturbed propagators $g_{{\textbf}{k}}(\tau)=\beta^{-1}\sum_l g_{{\textbf}{k}}(z_l) e^{-z_l \tau}$ we obtain the expression $$\begin{aligned}
  \label{Xistar2}
\mathscr{G}_{2,{\textbf}{k}}(\tau) 
&=
-\frac{1}{2} 
\sum_{{\textbf}{k}_2,{\textbf}{k}_3,{\textbf}{k}_4} \!\!
|\braket{\psi_{{\textbf}{k}}\psi_{{\textbf}{k}_2}| V|\psi_{{\textbf}{k}_3}\psi_{{\textbf}{k}_4}}|^2 \sum_{l,l_2,l_3,l_4}
{\nonumber \\}&\quad \times
\frac{1}{\beta^4} \operatorname{\operatorname{e}}^{-z_{l}\,\tau} 
\int \limits_{0}^{\beta} \! d\tau_1
\operatorname{\operatorname{e}}^{-(z_{l_3}+z_{l_4}-z_{l_2}-z_{l})\,\tau_1} 
{\nonumber \\}&\quad \times
\left[g_{{\textbf}{k}}(z_{l}) \right]^2\,g_{{\textbf}{k}_2}(z_{l_2})g_{{\textbf}{k}_3}(z_{l_3}) g_{{\textbf}{k}_4}(z_{l_4}),\end{aligned}$$ where we have eliminated the $\tau_2$ integral and one Matsubara sum via the relation $$\begin{aligned}
 \label{kroneck}
\frac{1}{\beta}\int \limits_{0}^{\beta} \! d\tau_2 \operatorname{\operatorname{e}}^{\pm(z_{l_3}+z_{l_4}-z_{l_2}-z_{l'})\,\tau_2}=\delta_{l_3+l_4,l_2+l'}.\end{aligned}$$ From $\mathscr{G}_{2,{\textbf}{k}}(\tau)=\beta^{-1}\sum_l \mathscr{G}_{2,{\textbf}{k}}(z_l) e^{-z_l \tau}$ and Eq.  we then find that $$\begin{aligned}
 \label{Xistar3}
\Xi_{2,{\textbf}{k}}[g_{\textbf}{k}(z_l),z_l]
&=
-\frac{1}{2} 
\sum_{{\textbf}{k}_2,{\textbf}{k}_3,{\textbf}{k}_4} \!\!
|\braket{\psi_{{\textbf}{k}}\psi_{{\textbf}{k}_2}| V|\psi_{{\textbf}{k}_3}\psi_{{\textbf}{k}_4}}|^2 \;
\sum_{l_2,l_3,l_4}
{\nonumber \\}&\quad \times
\frac{1}{\beta^3}\int \limits_{0}^\beta \! d\tau \operatorname{\operatorname{e}}^{-(z_{l_3}+z_{l_4}-z_{l_2}-z_{l})\,\tau} 
{\nonumber \\}&\quad \times
g_{{\textbf}{k}_2}(z_{l_2})g_{{\textbf}{k}_3}(z_{l_3}) g_{{\textbf}{k}_4}(z_{l_3}),\end{aligned}$$ i.e., $$\begin{aligned}
 \label{Xistar4}
\Xi_{2,{\textbf}{k}}[g_{\textbf}{k}(\tau),z_l]
&=
-\frac{1}{2} 
\sum_{{\textbf}{k}_2,{\textbf}{k}_3,{\textbf}{k}_4} \!\!
|\braket{\psi_{{\textbf}{k}}\psi_{{\textbf}{k}_2}| V|\psi_{{\textbf}{k}_3}\psi_{{\textbf}{k}_4}}|^2 \;
{\nonumber \\}&\quad \times
\int \limits_{0}^\beta \! d\tau \operatorname{\operatorname{e}}^{z_l\,\tau} 
g_{{\textbf}{k}_2}(\tau)g_{{\textbf}{k}_3}(\tau) g_{{\textbf}{k}_4}(\tau).\end{aligned}$$ Since $\tau>0$ in the time integral, from Eq.  we have $$\begin{aligned}
 \label{Xi2calc}
\Xi_{2,{\textbf}{k}}[n_{\textbf}{k},z_l]
&=
-\frac{1}{2} 
\sum_{{\textbf}{k}_2,{\textbf}{k}_3,{\textbf}{k}_4}
|\braket{\psi_{{\textbf}{k}}\psi_{{\textbf}{k}_2}| V|\psi_{{\textbf}{k}_3}\psi_{{\textbf}{k}_4}}|^2 
{\nonumber \\}&\quad \times 
\int \limits_{0}^{\beta} \! d\tau \operatorname{\operatorname{e}}^{-({\varepsilon}_{{\textbf}{k}_3}+{\varepsilon}_{{\textbf}{k}_4}-{\varepsilon}_{{\textbf}{k}_2}-z_l)\,\tau} n_{{\textbf}{k}_2}\bar n_{{\textbf}{k}_3}\bar n_{{\textbf}{k}_4},\end{aligned}$$ and carrying out the time integral we get Eq. .

### Functional relations

The functional relations between the perturbative contributions to the improper Matsubara self-energy and the grand-canonical potential are given by (see, e.g., Ref. [@Luttinger:1960ua]) $$\begin{aligned}
 \label{funcrel1a}
\Omega^{\aleph}_{n}[g_{{\textbf}{k}}(z_l)] &= \frac{1}{2n\beta} \sum_{\textbf}{k}  \sum_{l}
g_{{\textbf}{k}}(z_l)\,
\Xi^\star_{n,{\textbf}{k}}[g_{{\textbf}{k}}(z_l),z_l],
\\ \label{funcrel1b}
\Xi^\star_{n,{\textbf}{k}}[g_{{\textbf}{k}}(z_l),z_l]
&=
\beta\frac{\delta \Omega^{\aleph}_{n}[g_{{\textbf}{k}}(z_l)]}{\delta [g_{{\textbf}{k}}(z_l)]},\end{aligned}$$ and similar for $\Omega^{\aleph}_{n}[g_{{\textbf}{k}}(\tau)]$ and $\Xi^\star_{n,{\textbf}{k}}[g_{{\textbf}{k}}(\tau),\tau]$. The question is, what does $\aleph$ correspond to?

To find this out, we first evaluate the expression obtained from Eq.  for the second-order normal contribution, i.e., $$\begin{aligned}
 \label{O2calc1}
\Omega^{\text{direct}(P)}_{2,\text{normal}}[g_{{\textbf}{k}}(\tau)]
&=
\frac{1}{8 \beta} 
\sum_{{\textbf}{k}_1,{\textbf}{k}_2,{\textbf}{k}_3,{\textbf}{k}_4} \!\!
|\braket{\psi_{{\textbf}{k}}\psi_{{\textbf}{k}_2}| V|\psi_{{\textbf}{k}_3}\psi_{{\textbf}{k}_4}}|^2
W^{\text{direct}(P)}_{{\textbf}{k}_1,{\textbf}{k}_2,{\textbf}{k}_3,{\textbf}{k}_4},\end{aligned}$$ where $$\begin{aligned}
 \label{Wcalc1}
W^{\text{direct}(P)}_{{\textbf}{k}_1,{\textbf}{k}_2,{\textbf}{k}_3,{\textbf}{k}_4}
&=
\int \limits^{\beta}_0 \!\!d\tau_1 \int \limits^{\beta}_0 \!\! d\tau_2 \,
g_{{\textbf}{k}_1}(\tau_{21})g_{{\textbf}{k}_2}(\tau_{21})g_{{\textbf}{k}_3}(\tau_{12}) g_{{\textbf}{k}_4}(\tau_{12})
{\nonumber \\}&=
\int \limits^{\beta}_0 \!\!d\tau_1 \!\!\!\! \int \limits^{\beta-\tau_1}_{-\tau_1} \!\! d\tau' \,
g_{{\textbf}{k}_1}(\tau')g_{{\textbf}{k}_2}(\tau')g_{{\textbf}{k}_3}(-\tau') g_{{\textbf}{k}_4}(-\tau')
{\nonumber \\}&=
\int \limits^{\beta}_0 \!\!d\tau_1 \!\!\ \int \limits^{0}_{-\tau_1} \!\! d\tau' \,
\bar n_{{\textbf}{k}_1}\bar n_{{\textbf}{k}_2}n_{{\textbf}{k}_3} n_{{\textbf}{k}_4}
\operatorname{\operatorname{e}}^{-D\tau'}
{\nonumber \\}& \quad +
\int \limits^{\beta}_0 \!\!d\tau_1 \!\! \int \limits^{\beta-\tau_1}_{0} \!\!\!\! d\tau' \,
n_{{\textbf}{k}_1}n_{{\textbf}{k}_2}\bar n_{{\textbf}{k}_3} \bar n_{{\textbf}{k}_4}
\operatorname{\operatorname{e}}^{-D\tau'}
{\nonumber \\}&=
\bar n_{{\textbf}{k}_1}\bar n_{{\textbf}{k}_2}n_{{\textbf}{k}_3} n_{{\textbf}{k}_4} \frac{\beta D-1+\operatorname{\operatorname{e}}^{\beta D}}{D^2}
{\nonumber \\}& \quad -
n_{{\textbf}{k}_1}n_{{\textbf}{k}_2}\bar n_{{\textbf}{k}_3} \bar n_{{\textbf}{k}_4} \frac{\beta D-1+\operatorname{\operatorname{e}}^{-\beta D}}{D^2},\end{aligned}$$ with $D={\varepsilon}_{{\textbf}{k}_3}+{\varepsilon}_{{\textbf}{k}_4} - {\varepsilon}_{{\textbf}{k}_2} - {\varepsilon}_{{\textbf}{k}_1}$.[^53] Now, we can evaluate $W_{{\textbf}{k}_1,{\textbf}{k}_2,{\textbf}{k}_3,{\textbf}{k}_4}$ also by inserting in the first expression in Eq.  the Fourier expansion of the unperturbed propagators. This leads to $$\begin{aligned}
 \label{Wcalc2}
W^{\aleph}_{{\textbf}{k}_1,{\textbf}{k}_2,{\textbf}{k}_3,{\textbf}{k}_4}
&=
\frac{1}{\beta^3} 
\sum_{l_1,l_2,l_3,l_4} \delta_{l_3+l_4,l_2+l_1}
g_{{\textbf}{k}_1}(z_{l_1})
g_{{\textbf}{k}_2}(z_{l_2})g_{{\textbf}{k}_3}(z_{l_3}) g_{{\textbf}{k}_4}(z_{l_3}).\end{aligned}$$ Note that this is the expression we get by substituting $\Xi^\star_{2,{\textbf}{k}}[g_{{\textbf}{k}}(z_l),z_l]$ into Eq. . Using Eq.  we find $$\begin{aligned}
 \label{Wcalc3}
W^{\aleph}_{{\textbf}{k}_1,{\textbf}{k}_2,{\textbf}{k}_3,{\textbf}{k}_4}
&=
\beta
\int \limits_{0}^{\beta} \! d\tau \operatorname{\operatorname{e}}^{-({\varepsilon}_{{\textbf}{k}_3}+{\varepsilon}_{{\textbf}{k}_4}-{\varepsilon}_{{\textbf}{k}_2}-{\varepsilon}_{{\textbf}{k}_1})\,\tau} 
n_{{\textbf}{k}_1} n_{{\textbf}{k}_2}\bar n_{{\textbf}{k}_3}\bar n_{{\textbf}{k}_4}
{\nonumber \\}&= \beta n_{{\textbf}{k}_1} n_{{\textbf}{k}_2}\bar n_{{\textbf}{k}_3}\bar n_{{\textbf}{k}_4} \frac{\operatorname{\operatorname{e}}^{-\beta D}-1}{D},\end{aligned}$$ which corresponds to the *cyclic* expression, Eq. . However, we could have easily evaluated Eq.  such that the expression given by Eq.  would be obtained (i.e., by reversing the step that lead to Eq. ). Thus, the $\aleph$ in Eq.  depends on how the Matsubara sums are carried out. The identification of $\aleph$ with cyclic can however be fixed (formally) by substituting $\Xi^\star_{n,{\textbf}{k}}[n_{{\textbf}{k}},z_l]$ for $\Xi^\star_{n,{\textbf}{k}}[g_{{\textbf}{k}}(z_l),z_l]$, i.e., $$\begin{aligned}
\label{funcrel1c}
\Omega^\text{cyclic}_{n}[g_{{\textbf}{k}}(z_l)] &= \frac{1}{2n\beta} \sum_{\textbf}{k}  \sum_{l}
g_{{\textbf}{k}}(z_l)\,
\Xi^\star_{n,{\textbf}{k}}[n_{{\textbf}{k}},z_l].\end{aligned}$$ Note that no functional derivative relation is available for $\Xi^\star_{n,{\textbf}{k}}[n_{\textbf}{k},z_l]$.

Now, from Eqs.  and , we obtain by analytic continuation the relations $$\begin{aligned}
\label{funcrel2a}
\Omega_{n}^\aleph[g_{{\textbf}{k}}(z)] &= \frac{1}{2n} \sum_{\textbf}{k} 
\oint\limits_{C_0} 
\frac{d z}{2\pi {\text{i}}} \,
g_{{\textbf}{k}}(z)\,n_{{\textbf}{k}}(z)\,
\Xi^\star_{n,{\textbf}{k}}[g_{{\textbf}{k}}(z),z],
\\ \label{funcrel2b}
\Xi^\star_{n,{\textbf}{k}}[g_{{\textbf}{k}}(z')] 
&=
\frac{\delta \Omega_{n}^\aleph[g_{{\textbf}{k}}(z)]}{\delta[ g_{{\textbf}{k}}(z')]},\end{aligned}$$ where $C_0 \in \{C_0[l_\text{max}],C_0^{\pm}\}$, with $C_0[l_\text{max}]$ from Fig. \[figc1\] and $C_0^{\pm}$ from Fig. \[figc2\]. Note that these relations require that $\Xi^\star_{n,{\textbf}{k}}$ is represented as a functional of $g_{{\textbf}{k}}(z)$. Replacing $\Xi^\star_{n,{\textbf}{k}}[g_{{\textbf}{k}}(z),z]$ by the mass function $\mathcal{M}^\star_{n,{\textbf}{k}}[n_{\textbf}{k},z]$ leads to $$\begin{aligned}
\label{funcrel2c}
\Omega_{n}^\text{cyclic}[n_{\textbf}{k}] &= \frac{1}{2n} \sum_{\textbf}{k} 
\oint\limits_{C_0[l_\text{max}]} \!\!\!\!\!
\frac{d z}{2\pi {\text{i}}} \,
g_{{\textbf}{k}}(z)\,n_{{\textbf}{k}}(z)\,
\mathcal{M}^\star_{n,{\textbf}{k}}[n_{\textbf}{k},z].\end{aligned}$$ Because $g_{{\textbf}{k}}(z)\,n_{{\textbf}{k}}(z)\,
\mathcal{M}^\star_{n,{\textbf}{k}}[n_{\textbf}{k},z]$ vanishes at infinity we can deform the contour $C_0[l_\text{max}]$ into the contour $C_1$ from Fig. \[figc1\]. Since $\mathcal{M}^\star_{n,{\textbf}{k}}(z)$ is entire, we get only the contributions from the pole at $1/g_{{\textbf}{k}}(z)=0$, i.e., at $z={\varepsilon}_{\textbf}{k}$, so $$\begin{aligned}
\label{funcrel2d}
\Omega_{n}^\text{cyclic}[n_{\textbf}{k}] &= \frac{1}{2n} \sum_{\textbf}{k} n_{{\textbf}{k}}
\mathcal{M}^\star_{n,{\textbf}{k}}[n_{\textbf}{k},{\varepsilon}_{\textbf}{k}].\end{aligned}$$ Finally, as discussed above, from the expressions for $\mathcal{M}^\star_{n,{\textbf}{k}}[n_{\textbf}{k},z]$ the ones for the perturbative contributions to the frequency-space self-energy $\Sigma^\star_{n,{\textbf}{k}}[n_{\textbf}{k},z]$ are obtained by substituting $\operatorname{\operatorname{e}}^{\beta(z_l-\mu)}=-1$ and removing the remaining energy denominator exponentials via Eq. . From this, we find (analogous to zero-temperature MBPT [@Kaiser:2001ra; @Kaiser:2013bua]) that $$\begin{aligned}
\label{sigmastarred}
\Sigma^\star_{n,{\textbf}{k}}[n_{\textbf}{k},z] =  
\frac{\delta\Omega_n^{\text{reduced}}[n_{\textbf}{k}]}{\delta n_{\textbf}{k}}\bigg|_{{\varepsilon}_{\textbf}{k}= z},\end{aligned}$$ which implies the relation for the proper frequency-space self-energy $\Sigma_{n,{\textbf}{k}}(z)$ given by Eq. .

[^1]: In general, MBPT corresponds to divergent asymptotic series [@negele; @PhysRevLett.121.130405; @rossi; @refId0; @Marino:2019wra], so convergence rate should be understood in terms of the result at optimal truncation.

[^2]: To be precise, since Ref. [@PhysRevC.95.034326] uses the adiabatic formalism, the considered self-energy is not the frequency-space self-energy of the imaginary-time formalism, $\Sigma_{{{\textbf}{k}}}(z)$, but the collisional one $\Sigma^\text{coll}_{{\textbf}{k}}(\omega)$, which however satisfies $\text{Re}\,[\Sigma^\text{coll}_{{{\textbf}{k}}}(\omega)]=\text{Re}\,[\Sigma_{{{\textbf}{k}}}(\omega\pm {\text{i}}\eta)]$. We use the notion frequency somewhat generalized, i.e., by frequency we refer mostly to the argument $z\in \mathbb{C}$ of the (usual) analytic continuation $\Sigma_{{{\textbf}{k}}}(z)$ of the Matsubara self-energy $\Xi_{{{\textbf}{k}}}(z_l)$. See Appendix \[app21\] for details.

[^3]: Note that in Luttinger’s analysis [@PhysRev.174.263] of the scheme by Balian, Bloch, and de Dominicis it is incorrectly assumed that the mean field has the form $U[n_{\textbf}{k}]$. This is why in Ref. [@Fermiliq] Luttinger’s paper has been (incorrectly) associated with statistical quasiparticles.

[^4]: Several of the results presented here are also discussed in the authors dissertation [@Wellenhofer:2017qla], but note that some technical details have been missed and several typos appear there.

[^5]: In the Hartree-Fock case the self-consistency requirement can be evaded for isotropic systems at ${T=0}$ by replacing in the expression for $U_{{\textbf}{k}}=U_{1,k}$ the distribution functions $n_k=\theta(\mu-{\varepsilon}_k)$ by $\theta(k_{{\text{F}},\text{ref}}-k)$, where the unperturbed Fermi momentum $k_{{\text{F}},\text{ref}}$ is defined via ${{\varepsilon}_{k_{{\text{F}},\text{ref}}}=\mu}$. In that case, first-order MBPT is identical for ${U=0}$ and ${U=U_1}$ (more generally, ${U\propto U_1}$).

[^6]: In the thermodynamic limit the expressions for all size extensive quantities scale linearly with the confining volume. For notational simplicity, we neglect the scale factors. For discussions regarding our choice of basis states, see Refs. [@2016arXiv160900014C; @PhysRevA.95.062124].

[^7]: At ${T=0}$, in MBPT there are still divergences due to vanishing energy denominators, but these cancel each other at each order [@Wellenhofer:2018dwh] (see also Sec. \[sec44\]).

[^8]: The diagram composed of a single $-U$ vertex corresponds to $\Omega_{U}(T,\mu)$ and is excluded here. For the diagrammatic rules, see, e.g., Refs. [@szabo; @Runge].

[^9]: That is, normal pieces correspond to the linked normal subdiagrams of the normal unlinked diagram generated by cutting all anomalous articulation lines and closing them in each separated part.

[^10]: That is, diagrams with energy denominators involving only articulation lines with identical three-momenta. Such diagrams are anomalous. Note that one must distinguish between anomalous (normal) diagrams and anomalous (normal) contributions.

[^11]: In the case of the direct formula there are also pseudoanomalous contributions (of a different kind, i.e., terms $\sim T$) from the pole at $z=0$. Furthermore, in both the direct and the cyclic case the expressions for diagrams with several identical energy denominators involve terms $\sim T^{-\nu}$ with $\nu\geq 1$. In Hartree-Fock MBPT such diagrams appear first at sixth order, i.e., normal two-particle reducible diagrams composed of three second-order pieces. Because terms $\sim T^{-\nu}$ with $\nu\geq 1$ do not appear in the reduced formula, for the cyclic sums of diagrams these terms cancel each other in the ${T\rightarrow 0}$ limit.

[^12]: These singularities are present also in the adiabatic case for ${\varepsilon}_{\text{F}}\in\{{\varepsilon}_{\textbf}{k}\}$ (i.e., for open-shell systems).

[^13]: This inconsistency has been overlooked in Ref. [@SANTRA2017355].

[^14]: As discussed in Sec. \[sec23\], the Eq.  scheme leads to $\mu(T,\varrho)= \mu_\text{ref}(T,\varrho)$ for nonzero $T$ if the (pseudoanomalous) contributions from energy denominator poles are excluded, but it is not clear whether this is justified.

[^15]: Accordingly, for the grand-canonical series the Eq. \[Ufermi\] scheme does not remove the anomalous contributions to $\rho(T,\mu)$ at ${T=0}$ (or ${T\neq 0}$), so in that case the adiabatic series is not reproduced (in any case) and the discrete spectrum inconsistency persists.

[^16]: We note that $F(T,\mu_\text{ref})$ with $U\in\{\,0,U_1\}$ and ${N=2}$ has been employed in nuclear matter calculations in Refs. [@Tolos:2007bh; @Fritsch:2002hp; @Fiorilla:2011sr; @Holt:2013fwa; @PhysRevC.89.064009; @PhysRevC.92.015801; @PhysRevC.93.055802]. For nuclear matter calculations with self-consistent propagators, see, e.g., Refs. [@PhysRevC.88.044302; @PhysRevC.90.054322; @PhysRevC.98.025804].

[^17]: For additional details and numerical evidence, see Ref. [@Wellenhofer:2017qla].

[^18]: At second order the two types of anomalous contributions have been found to give individually very large but nearly canceling contributions in nuclear matter calculations [@Wellenhofer:2017qla; @Tolos:2007bh; @PhysRevC.89.064009].

[^19]: See also the next paragraph, and Sec. \[summary\].

[^20]: That is, the functional dependence on $n_{\textbf}{k}$ of the $-U$ vertices is *not* taken into account in Eq. .

[^21]: If the pole contributions are included for a finite system then Eq.  is valid only for ${T\rightarrow 0}$ (and the ${T\rightarrow 0}$ limit exists only for $\mu\not\in\{{\varepsilon}_{\textbf}{k}\}$). In that sense, the construction of the thermodynamic Fermi-liquid relations via MBPT depends on the thermodynamic limit.

[^22]: At ${T=0}$, these singular terms cancel each other, see Ref. [@Wellenhofer:2018dwh] and Sec. \[sec45\].

[^23]: Note that the relation ${\mathcal{E}_{{\textbf}{k}_{\text{F}}}=\partial E({T=0},\varrho)/\partial \varrho}$ (Hugenholtz-Van Hove theorem [@HUGENHOLTZ1958363; @Baym:1962sx]) is trivial if $E({T=0},\varrho)$ is derived from $\Omega(T,\mu)$.

[^24]: Note that this implies that there are diagrams with several identical energy denominators, i.e., the Hadamard finite part appears.

[^25]: A different regularization scheme can for example be set up via $\prod_\nu D_\nu^{n_\nu}\rightarrow (\prod_\nu D_\nu)^{n_\nu}+i\eta$. The parts $\mathcal{F}_{\alpha,[\eta]}$ then have a form that deviates from the reduced formula (in particular, the pseudoanomalous contributions do not cancel; see also Appendix \[app1\]), so the Fermi-liquid relations cannot be obtained in this scheme.

[^26]: Rules for the formal regularization have been presented also in Refs. [@BALIAN1971229; @LUTTINGER19731; @keitermorandi] for the case of impurity systems.

[^27]: It should be noted that, while the complete cancellation of two-particle reducible diagrams (with first-order pieces) is specific to $U_1$, including $U_2^\text{BdD}$, $U_3^\text{BdD}$, etc. does not only eliminate anomalous contributions but also partially cancels normal contributions. Note also that the reduced contributions from normal two-particle reducible diagrams with single-vertex loops can be resummed as geometric series; in zero-temperature MBPT this is equivalent to the change from ${U=0}$ to ${U=U_1}$ for isotropic systems (only).

[^28]: Notably, the same expression results if one naively introduces principal values in $\Omega_{3,\text{pp}}^\text{reduced}$ and averages over three different integration orders (where in one case $p$ is integrated before $A$ or $B$); for Eq.  this procedure would, however, lead to an incorrect result.

[^29]: In the context of MBPT for Fermi systems this method was introduced by Brout and Englert [@brout2; @PhysRev.115.824] (see also Ref. [@horwitz2]).

[^30]: The number operator is given by $\mathcal{N}=\sum_{\textbf}{k} a_{\textbf}{k}^\dagger a_{\textbf}{k}$.

[^31]: The mean field is not expanded; i.e., the expansion is performed after $U(T,\mu)$ is replaced by $U(T,\mu_\text{ref})$. This (and the truncation of the expansion) makes evident that at a given order the modified and the unmodified perturbation series lead to different results; see also Sec.\[sec23\] and Refs. [@PhysRevC.89.064009; @Wellenhofer:2017qla].

[^32]: Note that Eq. (B.12) of Ref. [@brout2] is not valid; e.g., it misses the second part of Eq. .

[^33]: This feature is expected from indirect arguments [@Luttinger:1960ua; @Wellenhofer:2017qla]. The cancellation has been shown explicitly to all orders for certain subclasses of diagrams [@Wellenhofer:2017qla], but no direct proof to all orders exists.

[^34]: This cancellation is not always purely algebraic, see Eqs.  and .

[^35]: For an interesting implication of this feature, i.e., the singularity at fourth order and ${T=0}$ of the Maclaurin expansion in terms of $x=\mu_\uparrow-\mu_\downarrow$ (or, $x=\varrho_\uparrow-\varrho_\downarrow$) for a system of spin one-half fermions with spins $\uparrow$ and $\downarrow$, see Refs. [@PhysRevC.91.065201; @PhysRevC.93.055802; @Wellenhofer:2017qla]. Note however that the statement in Refs. [@PhysRevC.93.055802; @Wellenhofer:2017qla] that the convergence radius of the expansion is still zero (instead of just very small) near (but not at) the degenerate limit appears somewhat questionable. In particular, Fig. 6of Ref. [@PhysRevC.93.055802] should be interpreted not in terms of the radius of convergence but in terms of convergence at $x=\pm 1$.

[^36]: See Ref. [@Wellenhofer:2018dwh] \[and Eq. \] for an example of this. We defer a more detailed analysis of these cancellations to a future publication.

[^37]: Other studies regarding the derivation of statistical quasiparticle relations can be found in Refs. [@nort1; @nort2; @article; @PhysRevA.1.1243; @RevModPhys.39.771; @PhysRev.153.263; @TUTTLE1966510; @1964quasip; @PhysRev.121.957; @PhysRev.124.583].

[^38]: Note also that somewhat similar methods have been applied with considerable success for (certain) finite systems [@Tichai:2018qge; @PhysRevLett.121.032501; @PhysRevLett.75.2787].

[^39]: More generally, the effect on convergence of higher-order contributions to the mean field (in the modified perturbation series for the free energy) will be investigated.

[^40]: We defer a more detailed analysis of these partial analytic cancellations to a future publication.

[^41]: Here, the functional derivative is supposed to disregard the implicit dependence on $n_{\textbf}{k}$ of $U_{2}^\text{BdD}$; see Sec. \[sec23\].

[^42]: See the second paragraph of part \[app21\] of this Appendix.

[^43]: In that sense, the mass function representation (as well as the direct representation of part \[app22\]) represents a purely statistical result, while the spectral representation corresponds to a statistical-dynamical result. Only the statistical-dynamical result has a well-behaved ${T\rightarrow 0}$ limit.

[^44]: In this paragraph we follow for the most part Kadanoff and Baym [@kadanoffbaym], Fetter and Walecka [@Fetter], and Ref. [@ThesisRios].

[^45]: Note that Eqs.  and imply that $\text{Im}[\Sigma_{\textbf}{k}(z)]\lessgtr 0$ for $\text{Im}[z] \gtrless 0$.

[^46]: Note that Fetter and Walecka omit the factor ${\text{i}}$ in the definition of $G_{\textbf}{k}(\omega)$, so no ${\text{i}}$ appears in their version of our Eq. , i.e., in Eq. (9.33) of Ref. [@Fetter].

[^47]: In the adiabatic formalism only real-time propagators appear, so it is the collisional self-energy that is calculated.

[^48]: We note again that in the direct (and cyclic) scheme the ${T\rightarrow 0}$ limit is nonexistent.

[^49]: Because of the two external lines no cyclic and reduced versions of Eq.  are available; see Ref. [@BDDnuclphys7] for details on the derivation of the cyclic formula and the reduced formula.

[^50]: We use the notion normal propagator diagrams to refer to diagrams that have no anomalous articulation lines and are either (i) one-particle irreducible propagator diagrams or (ii) one-particle reducible propagator diagrams where all cuttable propagator lines go in the same direction.

[^51]: In the cyclic and the BdD scheme only the contributions with Hugenholtz diagrams attached via higher-cumulant connections can be canceled. The remaining propagator contributions in these schemes are then given by Eqs. , , etc., with the $-U$ vertices (but not the self-energy parts) given by $-U^{\text{cyclic},\div}_n$ and $-U^{\text{BdD}}_n$, respectively, plus diagrams that have self-energy parts consisting of $-U^{\text{cyclic},\div}_n$ and $-U^{\text{BdD}}_n$ vertices (with $2\leq n\leq N$), respectively.

[^52]: For the self-consistent functional relations between the proper self-energy and the grand-canonical potential, see, e.g., Refs. [@Abrikosov; @Luttinger:1960ua; @Baym:1962sx].

[^53]: It can be seen by regularizing the energy denominators that the expression obtained from Eqs.  and is equivalent to the direct, cyclic, and regularized reduced expressions for the (permutation invariant) second-order normal diagram; see Sec. \[sec31\].
