# Strong coupling quantum impurity solver on the real and imaginary axis

Kristjan Haule  Center for Materials Theory, Department of Physics and
Astronomy, Rutgers University, Piscataway, NJ 08854, USA

###### Abstract

The diagramatic Monte Carlo method has so far been primarily used in
connection with the weak coupling expansion. Here we show that the strong
coupling expansion offers a significant advantage: it can be efficiently
implemented on both the real and the imaginary axis at finite temperature.
Using the example of a quantum impurity solver for the Dynamical Mean Field
Theory (DMFT) problem, we illustrate rapid convergence with respect to the
expansion order. We derive a closed-form expression for the Feynman diagrams
of arbitrary order on the real axis. Employing these Feynman rules, we
implement the bold hybridization-expansion quantum Monte Carlo (BHQMC)
impurity solver and compare its performance to state-of-the-art results from
Numerical Renormalization Group calculations of the Mott transition within
DMFT applied to the Hubbard model. We demonstrate its power in providing a
very precise frequency dependent scattering rate at finite temperature,
enabling accurate spectroscopy calculations and delivering benchmark results
for transport within DMFT.

Introduction. Numerous quantum many body problems can be solved by expanding
the action in powers of certain parameter and summing the resulting terms to
high orders in perturbation theory. This has been successfully demonstrated
through the diagrammatic Monte Carlo method [1, 2] in various physical
systems, such as the problem of the uniform electron gas [3, 4], the unitary
gas [2], the polaron problem [5], the quantum impurity problem [6, 7, 8, 9],
and even the Hubbard model in some parameter regime [1, 10, 11, 12]. Most of
these expansions are carried out from so-called weak coupling starting point,
where the unperturbed action is quadratic, and the application of the Wick’s
theorem is straightforward. In contrast, the strong coupling expansion is less
commonly attempted due to the absence of Wick’s theorem and the resultant
complexity of the diagrammatic rules [13, 14]. Nevertheless, these rules are
well known [14, 15] and have been carried out for many quantum problems. In
particular the hybridization expansion continuous time quantum Monte Carlo
(CTQMC) method [7, 8] is an example of a very successful strong coupling
expansion to very high orders, which has proven to be the most robust
numerically exact quantum impurity solver. This solver is implemented in the
imaginary time formalism, which requires numerical analytic continuation
method to obtain the observables on the real-frequency axis. As the
continuation method is mathematically ill-posed and numerically unstable,
accurate spectra on the real axis are rarely available. To overcome this
issue, the expansion on the real axis was recently developed for the weak
coupling approach [16, 17, 18, 19], which does not require analytic
continuation and produces highly accurate spectra directly on the real axis.
However, this method exponentially increases computational cost with expansion
order, as it generates at least $2^{n}$ terms at expansion order $n$ from a
single term in imaginary axis.

Using example of the quantum impurity problem, we will here show that in
contrast to the weak coupling expansion, the strong coupling expansion on the
real axis can be carried out with similar efficiently as on the imaginary
axis. We will derive a simple closed form expressions for the Feynman
diagrams, which closely resemble the expansion in the Matsubara formalism.

We carry out the bold strong coupling expansion on both the imaginary and the
real axis, and demonstrate that highly accurate real axis spectra can be
obtained using this method, comparable to established state of the art
Numerical Renormalization group (NRG) results. This method is similar to the
conventional hybridization CTQMC method [7, 8], with the difference that the
atomic propagators are here boldified and self-consistently determined. The
computational effort is comparable to the bare expansion, but is expected to
have much less severe minus sign problem in cluster DMFT applications as the
expansion order is much smaller. More importantly, this method delivers both
the imaginary axis and the real axis observables. The former is best suited to
calculate the total charge density, total energy and forces in realistic
materials [20, 21, 22], while the latter enables precise calculations of
response functions for theoretical spectroscopy.

Action. The strong coupling expansion starts from the sum of the atomic action
$S_{atm}$, which is solved exactly, and the perturbative part, which is
expanded in power series [15]. For the case of quantum impurity model, the
latter is the hybridization $S_{\Delta}$ of the form

$S_{\Delta}=\int_{0}^{\beta}\int_{0}^{\beta}d\tau
d\tau^{\prime}\sum_{\alpha,\beta}\psi^{\dagger}_{\alpha}(\tau)\Delta_{\alpha,\beta}(\tau-\tau^{\prime})\psi_{\beta}(\tau^{\prime}),$
(1)

and $\Delta_{\alpha,\beta}(\tau-\tau^{\prime})$ represents the matrix of the
hybridization function. Here, $\alpha$ and $\beta$ are orbital and spin
indices, $\psi$ represents Grassmann variables associated with the electron
annihilation operator.

We denote the energy levels of the atomic Hamiltonian ($H_{atm}$) as $E_{m}$
and its eigenvectors as $\ket{m}$. The Feynman perturbative expansion in
powers of $S_{\Delta}$ can be carried out conveniently by rewriting the atomic
action in terms of the quadratic action
$S_{atm}=\int_{0}^{\beta}d\tau\sum_{m}a_{m}^{\dagger}(\tau)\left(\partial_{\tau}+E_{m}+\lambda\right)a_{m}(\tau)$,
where the atomic eigenstate $\ket{m}$ is created by the pseudo-particle
operator from a new vacuum state, denoted as $\ket{m}\equiv
a^{\dagger}_{m}\ket{0}$. [15, 23] The atomic state with even (odd) number of
electrons is pseudo boson (fermion) [15]. The variable $\lambda$ serves as a
chemical potential for pseudo-particles and is used later to enforce
projection onto the physical Hilbert space. The entire action can now be
rewritten in terms of these pseudo-particles, and the hybridization part takes
the form of a retarded interaction:

$\displaystyle
S_{\Delta}=\sum_{\alpha\beta,mnm^{\prime}n^{\prime}}\braket{m}{\psi_{\alpha}^{\dagger}}{n}\braket{m^{\prime}}{\psi_{\beta}}{n^{\prime}}\times$
(2) $\displaystyle\int_{0}^{\beta}\int_{0}^{\beta}d\tau
d\tau^{\prime}a_{m}^{\dagger}(\tau)a_{n}(\tau)\Delta_{\alpha,\beta}(\tau-\tau^{\prime})a_{m^{\prime}}^{\dagger}(\tau^{\prime})a_{n^{\prime}}(\tau^{\prime}).$

In the quantum impurity model it is a local interaction, while in the more
general lattice model it is non-local retarded interaction. The total action
$S=S_{atm}+S_{\Delta}$ conserves the pseudo-particle number on each given
site, but allows any number of pseudo-particles in the system. This pseudo
representation thus expands the Hilbert space, incorporating both the physical
and unphysical part. The latter hence needs to be projected out. The advantage
of this representation is that Wick’s theorem in the grand canonical space is
valid, hence the expansion in terms of Feynman diagrams is straightforward.

Projection. Because the atomic eigenstates satisfy the completeness relation
$\sum_{m}\ket{m}\bra{m}=1$, the physical part of the Hilbert space allows for
the existence of exactly one pseudo-particle per correlated site in the system
at any given time. This is because the completness relation in pseudo-particle
representation is $\sum_{m}a^{\dagger}_{m}a_{m}=1$. If we define pseudo-charge
$Q=\sum_{m}a_{m}^{\dagger}a_{m}$, the projection than requires $Q=1$. The
expanded grand canonical space permits any number of pseudo-particles, and the
extraction of the $Q=1$ part is accomplished through the chemical potential
for pseudo-particles, $\lambda$, introduced earlier, which must approach
infinity. [15, 24, 25] It’s worth noting that this chemical potential is
simply added to the action (or the Hamiltonian) as $H\rightarrow H+\lambda Q$,
effectively splitting the pseudo-particle spectra corresponding to different
occupations $Q$ by the value $\lambda$. $Q$ is constant in time as $[H,Q]=0$.
If $G_{p}(\omega)$ represents the pseudo Green’s function on the real axis,
the spectra corresponding to charge $Q$ appears at $G_{p}(\omega+Q\lambda)$,
where $\omega$ is small (at most of the order of bandwidth), and $\lambda$
tends towards infinity. The pseudo-particle spectra and the impurity spectra
vanish at $Q=0$ (because this corresponds to the absence of atomic degrees of
freedom). Therefore, the first nonzero spectra appears at $\omega+\lambda$,
corresponding to the physical $Q=1$ spectra. All higher charging states are
removed by setting $\lambda\rightarrow\infty$. The projection onto the
physical space (when calculation is performed on the real axis) is thus
achieved by shifting the frequency integrations of all pseudo-Green’s
functions such that they contain a single $\lambda$ in their argument, i.e.,
$G_{p}(\omega+\lambda)$, while the hybridization functions $\Delta(\omega)$
should not be shifted. More details on projection can be found in Sect.1 of
Suppl. [25].

Figure 1:  a) A third order Feynman diagram with frequency loops in color
coding. b) Probability for perturbation order in bare and bold expansion of
action Eq. 2 at $U=8t$, $T=0.04t$ and DMFT solution of the 2D square lattice
Hubbard model. c) Contributions to the spectral functions from different
orders for the same parameters. d) The imaginary part of the self-energy of
the 2D square lattice and Bethe lattice solved by DMFT and $U=9t$, $T=0.04t$.
Black thin line stands for the linear function of frequency.

Expansion gets simplified by projection. While enforcing the constraint
through projection ($Q=1$) may initially seem like a tedious task, it turns
out that it greatly simplifies the perturbative expansion. It limits the
topology of Feynman diagrams to those with a single loop of pseudo-particles
per correlated site with strict time ordering, which limits quantum impurity
diagrams to those with a single backbone (example in Fig. 1a). While this
simplification has been well known, its impact on the analytic continuation of
generic Feynman diagrams was not appreciated. Specifically, a single diagram
at order $n$ contains $n$ hybridization lines, and its expression in grand-
canonical space along the imaginary axis can be readily written down following
standard Feynman rules. However, in the process of the analytic continuation
to the real axis, the number of terms normally balloons to at least $2^{n}$,
as each Matsubara sum necessitates integration across all branch cuts of
propagators in the loop (at least 2).

We will demonstrate that the projection in this formulation of the problem
results in an enormous simplification, where only a single term (or at most
$2n$ for physical spectral functions) survives the projection. This implies
that the evaluation of Feynman diagrams on both the real and imaginary axis
becomes equally straightforward within a Monte Carlo random walk framework.
Both random samplings are equally efficient.

To demonstrate the existence of a concise expression for each Feynman diagram
on the real axis, we will perform a generic Matsubara sum for a given
frequency loop and establish that all but one integration over the branch cut
vanish due to the projection. A generic closed loop diagram (see example in
Fig.1a) has $n+1$ frequency loops, and the self-energy contribution, which is
a functional derivative, requires $n$ Matsubara summations. We select
Matsubara frequencies on hybridization propagators to be simple fermionic
Matsubara frequencies, either $\nu_{m}$ or $-\nu_{m}$, where
$\nu_{m}=(2m+1)\pi T$, with $m$ as an integer and $T$ as the temperature (see
color coded loops in Fig.1a). To evaluate any pseudo self-energy, the
summation over $n$-fermionic $\nu_{m}$, associated with the hybridizations,
should be performed. Assume that we already performed summation over several
$i\nu^{\alpha}_{m}$, which were replaced by integral over real variables
$y_{\alpha}$, and we still need to carry out the summation over
$i\nu^{\beta}_{l}$, while $i\nu_{m}$ Matsubara loop is the focus of this step.
We have

$\displaystyle T\sum_{i\nu_{m}}\Delta(\pm
i\nu_{m})\prod_{j=1}^{N}G_{j}(i\nu_{m}+i\omega+\sum_{\alpha}y_{\alpha}a_{j\alpha}+\sum_{\beta}i\nu_{l}^{\beta}a_{j\beta})$

Here, $N$ represents the number of pseudo Green’s functions in the considered
loop $\nu_{m}$, $a_{j\alpha}$ and $a_{j\beta}$ can take values of either $1$
or $-1$. Here $y_{\alpha}$ appear as arguments in
$\Delta_{\alpha}(y_{\alpha})$ already summed over.

The Matsubara summation entails contour integral over the complex plain
encompassing all Matsubara frequencies while avoiding all branch cuts of the
integrand, including the branch-cut of $\Delta$, which appears on the real
axis, and all other branch cuts of the $G_{j}$, located off the real axis
[25]. This process generates $N+1$ terms in this step. However, we will show
that only the branch cut of $\Delta$ survives the projection, resulting in the
expression:

$\pm\int dyA_{c}(\pm
y)f(y)\prod_{j=1}^{N}G_{j}(y+i\omega+\sum_{\alpha}y_{\alpha}a_{j\alpha}+\sum_{\beta}i\nu_{l}^{\beta}a_{j\beta}).$
(3)

Here, $A_{c}=-\textrm{Im}(\Delta)/\pi$ represents the spectral function of the
hybridization. It’s worth noting that this term remains finite and unchanged
as the limit $\lambda\rightarrow\infty$ is taken, because the common frequency
in all $G_{j}$ is $\omega$, which is shifted by $\lambda$ to infinity through
projection, and allows all arguments of the hybridization function and the
Fermi functions to remain finite. Next we demonstrate that the other terms
generated by this Matsubara sum vanish. Let’s consider the branch-cut
associated with function $G_{j}$. Initially, we replace this specific $G_{j}$
with its spectral representation: $G_{j}(z)=\int
dx_{j}A_{j}(x_{j})/(z-x_{j})$, where $A_{j}$ is the pseudo-particle spectral
function. We then evaluate its contour integral, yielding:

$\int
dx_{j}A_{j}(x_{j})f(x_{j}-\sum_{\alpha}y_{\alpha}a_{j\alpha})R(x_{j},i\omega,y_{\alpha},i\nu_{l}^{\beta})$
(4)

Here, $R(\cdots)$ contains other terms generated by the sum, which we omit for
brevity, but note they are finite [25]. The crucial observation is that the
projection requires shifting $x_{j}$ by $\lambda$ because it appears as the
argument of the pseudo-spectral function $A_{j}$, while all $y_{\alpha}$ are
unshifted, as they are arguments of the hybridization functions. As a result
we have $A(x_{j}+\lambda)f(x_{j}+\lambda-O(1))$. Consequently, the Fermi
function $f(x_{j}+\lambda-\cdots)$ evaluates to $e^{-\beta\lambda-\beta
x_{j}}$, and vanishes after the projection. Thus, we have shown that each
summation over Matsubara frequency associated with the hybridization function
generates a single term on the real axis, easily obtained using the Feynman
rules, specified in [25]. The pseudo self-energies do not require summation
over the frequency $i\omega$, and as a result, the final expression always
consists of a single term. On the other hand, the impurity Green’s function
necessitates summation over $i\omega$ and generates $2n$ terms after
projection. This is negligible compared to over $2^{n}$ terms in the weak
coupling expansion.

Algorithm. We have implemented a diagrammatic Monte Carlo (MC) algorithm for
sampling arbitrary high orders in the expansion of the action Eq. 2. This
algorithm samples skeleton diagrams, meaning that the pseudo-propagators are
fully dressed with self-energy, which is obtained through the functional
derivative of the Luttinger-Ward functional. The electron Green’s function is
computed by taking the functional derivative with respect to the hybridization
function [25]. In our numerical sampling, a MC configuration consists of a set
of up to four Luttinger-Ward functionals (an example is provided in Fig. 1a).
These functionals are differentiated on the fly to obtain all physical
observables. Sampling only skeleton diagrams presents a considerable challenge
because the space of skeletons is not well-connected, and simple moves do not
create an ergodic random walk. To address this challenge, we have employed the
following algorithm: we utilize rejection-free MC sampling technique [26],
where we keep up to four diagrams within a single MC configuration to improve
ergodicity [26]. Our MC steps involve changing time/frequency, adding or
removing a single hybridization, or exchanging any two vertices where
hybridization starts or ends. It’s important to note that most diagram-
changing steps generate non-skeleton diagrams, which we skip over until a
skeleton diagram is found. This approach ensures that the space of diagrams is
simply connected, much like in conventional CTQMC. However, in this algorithm,
we need to compute the volume of the space of skeleton diagrams, which we
achieve by keeping track of each unique diagram visited during the random
walk. This allows us to properly normalize the results by knowing the volume
of the reduced phase space of skeleton diagrams.

In Fig. 1b, we compare the probability for expansion order of the bare
expansion (conventional CTQMC [7, 8]), and the bold expansions, i.e., BHQMC.
As illustrated, the bold expansion converges significantly faster compared to
the bare expansion. The average order $n$ of the latter is at
$|\mathrm{Tr}(\Delta G)|/T$, which is the kinetic energy divided by
temperature, as we are sampling the partition function $Z$. On the other hand,
the average $n$ of the bold algorithm corresponds to the logarithm of the same
expression, since we are sampling $\log(Z)$ rather than $Z$. Therefore the
actual efficiency of the bold algorithm is better than bare expansion for the
systems tested here, even though hybridizations can not be grouped into
determinants here.

In Fig. 1c, we present the contribution of each order to the electron Green’s
function, and it’s evident that higher-order contributions fall off rapidly
and are more concentrated near zero frequency. In Fig. 1d, we display the
electron self-energy near zero frequency to demonstrate that the Fermi liquid
regime is attained at this temperature, with a scattering rate that is
quadratic at low frequencies. In our work, we compute the electron self-energy
through the two-particle response function using the Bulla trick [27] (see
[25]). This method facilitates rapid convergence of the electron self-energy
with perturbation order. We also compare the Bethe lattice self-energy with
that of the 2D square lattice to illustrate that the former exhibits a very
high coherence scale with a quadratic scattering rate over a wide frequency
range. However, the 2D square lattice features a Van-Hove singularity at half-
filling, which results in a self-energy that is quasi-linear in frequency.
This phenomenon has been discussed in previous works on the same [28] and
related model [29].

Figure 2: The DMFT spectral function for the 2D square lattice Hubbard model.
Comparison between published NRG results from Ref. [28] with thin lines and
current method with bold lines. We note that Mott transition at $T=0.04t$ is
around $U=10t$ [30], while at zero temperature is $11.6t$ [28].

In Fig. 2, we show the spectral function across the DMFT Mott transition for
the 2D square lattice. The bold lines represent the results by BHQMC, while
the thin lines correspond to the NRG results, reproduced from Ref. [28].
Remarkably, we observe excellent agreement across the entire frequency range
and for all interaction strengths, except at $U=11.2t$ where NRG still has the
narrow peak while in BHQMC the quasiparticle peak is absent. This disagreement
arises because the NRG results are obtained at zero temperature and are below
the critical $U$ ($U_{c}(T=0)=11.6t$ [28]), wheres the BHQMC results are
obtained at $T=0.04t$ ($U_{c}(T=0.04t)\approx 10t$ [30]), and are above the
corresponding critical $U$. It is quite interesting to see that the Hubbard
bands around the Mott transition ($U=11.2t$) are essentially the same at
finite and zero temperature. The agreement for the rest of the parameters is
excellent. It’s also worth noting that the NRG seems to slightly over-broaden
features at large $U$ near the edges of the Hubbard bands, whereas the BHQMC
seems to provide a more abrupt and fine details at the edges of Hubbard bands.

Figure 3:  a) The resistivity of the 2D Hubbard model as measured by the cold
atom experiment [31]. (red symbols) compared to finite temperature Lanczos
method (blue thick line) and the DMFT method using exact diagonalization
(green thick line). The green line with dots correspond to BHQMC. Here
$U=7.5t$ and doping $n_{d}=0.825$.

The computation of transport quantities within DMFT is particularly
challenging as it demands highly accurate self-energy on the real axis. Recent
advancements in cold atom experiments have allowed the measurement of
resistivity in the 2D Hubbard model at high temperatures, and we reproduce
these results from Ref. [31] in Fig. 4a (red symbols). In the same work a
comparison was drawn with theoretical predictions of the finite-temperature
Lanczos method (blue thick line) and DMFT (green thick line). In this study,
we recalculated the DMFT curve (green dots) and compared it with earlier
estimates derived from a finite-size exact diagonalization solver. Notably,
the latter method is sensitive to the broadening of spectral peaks. In Fig.
4b, we present the corresponding spectral function and its temperature
evolution. Our precise values of resistivity align well with the earlier
results at low temperatures but at high temperatures are somewhat larger than
previously obtained approximate results. The current understanding behind the
discrepancy between experiment and DMFT is that the current vertex corrections
in 2D are significant, even at elevated temperatures [32, 33].

Conclusions. The diagramatic Monte Carlo method has been so far mainly used in
connection with the weak coupling expansion. Here we showed that the strong
coupling expansion has an important advantage: it is rapidly converging when
used in connection with DMFT as the hybridization is a small parameter, and
can be very efficiently implemented on the real axis as the closed form
expression for Feynman diagrams exists. By implementing the bold expansion
solver on both the real and the imaginary axis we demonstrated that the finite
frequency spectroscopies can be calculated very accurately within the DMFT at
finite temperature. The comparison with NRG results illustrates the
exceptional accuracy of BHQMC. The extension in several directions is now
possible: a) the multi-orbital imurity solver, as needed for combination of
DFT and DMFT applied to realistic materials; b) the cluster-DMFT solver for
which conventional CTQMC shows a strong fermionic sign problem, such as the
p-d model of cuprates; c) the steady state non-equilibrium problem on the real
axis.

Acknowledgements. We acknowledge support of NSF DMR-2233892 and NSF
OAC-2311557.

## References

  * Kozik  _et al._ [2010] E. Kozik, K. V. Houcke, E. Gull, L. Pollet, N. Prokof’ev, B. Svistunov, and M. Troyer, Europhysics Letters 90, 10004 (2010).
  * Van Houcke  _et al._ [2012] K. Van Houcke, F. Werner, E. Kozik, N. Prokof’ev, B. Svistunov, M. J. H. Ku, A. T. Sommer, L. W. Cheuk, A. Schirotzek, and M. W. Zwierlein, Nature Physics 8, 366 (2012).
  * Chen and Haule [2019] K. Chen and K. Haule, Nature Communications 10, 3725 (2019).
  * Haule and Chen [2022] K. Haule and K. Chen, Scientific Reports 12, 2294 (2022).
  * Mishchenko  _et al._ [2014] A. S. Mishchenko, N. Nagaosa, and N. Prokof’ev, Phys. Rev. Lett. 113, 166402 (2014).
  * Rubtsov  _et al._ [2005] A. N. Rubtsov, V. V. Savkin, and A. I. Lichtenstein, Phys. Rev. B 72, 035122 (2005).
  * Werner  _et al._ [2006] P. Werner, A. Comanac, L. de’ Medici, M. Troyer, and A. J. Millis, Phys. Rev. Lett. 97, 076405 (2006).
  * Haule [2007] K. Haule, Phys. Rev. B 75, 155113 (2007).
  * Gull  _et al._ [2008] E. Gull, P. Werner, O. Parcollet, and M. Troyer, Europhysics Letters 82, 57003 (2008).
  * Rossi [2017] R. Rossi, Phys. Rev. Lett. 119, 045701 (2017).
  * Wu  _et al._ [2017] W. Wu, M. Ferrero, A. Georges, and E. Kozik, Phys. Rev. B 96, 041105(R) (2017).
  * Schäfer  _et al._ [2021] T. Schäfer, N. Wentzell, F. Šimkovic, Y.-Y. He, C. Hille, M. Klett, C. J. Eckhardt, B. Arzhang, V. Harkov, F. m. c.-M. Le Régent, A. Kirsch, Y. Wang, A. J. Kim, E. Kozik, E. A. Stepanov, A. Kauch, S. Andergassen, P. Hansmann, D. Rohe, Y. M. Vilk, J. P. F. LeBlanc, S. Zhang, A.-M. S. Tremblay, M. Ferrero, O. Parcollet, and A. Georges, Phys. Rev. X 11, 011058 (2021).
  * Keiter and Kimball [1971] H. Keiter and J. Kimball, International Journal of Magnetism 1, 233 (1971).
  * Grewe and Keiter [1981] N. Grewe and H. Keiter, Phys. Rev. B 24, 4420 (1981).
  * Coleman [1984] P. Coleman, Phys. Rev. B 29, 3035 (1984).
  * Taheridehkordi  _et al._ [2019] A. Taheridehkordi, S. H. Curnoe, and J. P. F. LeBlanc, Phys. Rev. B 99, 035120 (2019).
  * LeBlanc  _et al._ [2022] J. P. F. LeBlanc, K. Chen, K. Haule, N. V. Prokof’ev, and I. S. Tupitsyn, Phys. Rev. Lett. 129, 246401 (2022).
  * Vucicevic and Ferrero [2020] J. Vucicevic and M. Ferrero, Phys. Rev. B 101, 075113 (2020).
  * Vučičević  _et al._ [2021] J. Vučičević, P. Stipsić, and M. Ferrero, Phys. Rev. Res. 3, 023082 (2021).
  * Kotliar  _et al._ [2006] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Rev. Mod. Phys. 78, 865 (2006).
  * Haule and Pascut [2016] K. Haule and G. L. Pascut, Phys. Rev. B 94, 195146 (2016).
  * Haule and Birol [2015] K. Haule and T. Birol, Phys. Rev. Lett. 115, 256402 (2015).
  * Haule  _et al._ [2010] K. Haule, C.-H. Yee, and K. Kim, Physical Review B 81, 195107 (2010).
  * Abrikosov [1965] A. A. Abrikosov, Physics Physique Fizika 2, 5 (1965).
  * [25] “Supplemental material,” http://www.example.com/supplementary_material.
  * Šimkovic and Rossi [2021] F. Šimkovic and R. Rossi, arXiv preprint arXiv:2102.05613 (2021).
  * Bulla  _et al._ [1998] R. Bulla, A. C. Hewson, and T. Pruschke, Journal of Physics: Condensed Matter 10, 8365 (1998).
  * Žitko  _et al._ [2009] R. Žitko, J. Bonča, and T. Pruschke, Phys. Rev. B 80, 245112 (2009).
  * Zang  _et al._ [2022] J. Zang, J. Wang, J. Cano, A. Georges, and A. J. Millis, Phys. Rev. X 12, 021064 (2022).
  * Park  _et al._ [2008] H. Park, K. Haule, and G. Kotliar, Phys. Rev. Lett. 101, 186403 (2008).
  * Brown  _et al._ [2019] P. T. Brown, D. Mitra, E. Guardado-Sanchez, R. Nourafkan, A. Reymbaut, C.-D. Hebert, S. Bergeron, A.-M. S. Tremblay, J. Kokalj, D. A. Huse, P. Schauss, and W. S. Bakr, Science 363, 379 (2019), https://www.science.org/doi/pdf/10.1126/science.aat4134 .
  * Vučičević  _et al._ [2019] J. Vučičević, J. Kokalj, R. Žitko, N. Wentzell, D. Tanasković, and J. Mravlje, Phys. Rev. Lett. 123, 036601 (2019).
  * Vranić  _et al._ [2020] A. Vranić, J. Vučičević, J. Kokalj, J. Skolimowski, R. Žitko, J. Mravlje, and D. Tanasković, Phys. Rev. B 102, 115142 (2020).
  * Tupitsyn  _et al._ [2021] I. S. Tupitsyn, A. M. Tsvelik, R. M. Konik, and N. V. Prokof’ev, Phys. Rev. Lett. 127, 026403 (2021).

##  I Supplementary information

###  I.1 Spectral function

Figure 4:  The temperature evolution of the single-particle spectral function
for $U=7.5t$ and doping $n_{d}=0.825$, parameters from Fig.3.

The resistivity of the 2D Hubbard model was shown Fig.3 for parameters
$U=7.5t$ and doping $n_{d}=0.825$, and reproduced from Ref. [31]. The
temperature evolution of the corresponding single-particle spectral function
for the same parameters is shown here. It demonstrates a gradual reduction in
the pseudo-particle peak, contrasting with the typical sudden collapse
observed in cases of analytic continuation from imaginary axis data.

###  I.2 Projection to the physical Hilbert space

By introducing the pseudo-particle creation operator $a_{m}^{\dagger}$ such
that each generates an atomic eigen-state $\ket{m}$, i.e.,
$a_{m}^{\dagger}\ket{0}=\ket{m}$, the Hilbert space is substantially expanded.
In the physical part of the Hilbert space, only one atomic state is allowed,
represented by the constraint
$Q=\sum_{m}a_{m}^{\dagger}a_{m}=\sum_{m}\ket{m}\bra{m}=1$. To eliminate the
unphysical portion of the Hilbert space, an exact projection to the $Q=1$
Hilbert subspace is necessary. This is achieved by introducing a chemical
potential for pseudo-particles, i.e., $H\rightarrow H+\lambda Q$. It is then
realized that at large $\lambda$, the excitations corresponding to different
Hilbert spaces with different $Q$ will be separated by approximately
$\lambda$. The physical excitations of $Q=1$ will, therefore, appear at a
frequency of $\omega+\lambda$, where $\omega$ is of the order of the
bandwidth, and $\lambda$ tends to infinity. Consequently, all frequencies in
pseudo-particle propagators need to be shifted by $\lambda$ as
$G_{p}(\omega+\lambda)$.

In addition to shifting the frequency variable in all pseudo-particle
quantities, it is essential to ensure that physical observables are computed
by taking the proper trace over the physical Hilbert space. Before the
projection, we calculate physical observables in the grand-canonical ensemble,
where any value of $Q$ is allowed, and is represented as:

$\braket{A}=\frac{\mathrm{Tr}(e^{-\beta(H+\lambda
Q)}A)}{\mathrm{Tr}(e^{-\beta(H+\lambda Q)})},$ (5)

As $Q$ is a conserved quantity, we can separately perform the trace over
states with $Q=0$, $Q=1$, and so on. These traces are denoted as
$\mathrm{Tr}_{Q=0}$, $\mathrm{Tr}_{Q=1}$, and so forth. Consequently, we
obtain:

$\braket{A}=\frac{\mathrm{Tr}_{Q=0}(e^{-\beta
H}A)+e^{-\beta\lambda}\mathrm{Tr}_{Q=1}(e^{-\beta
H}A)+O(e^{-2\beta\lambda})}{\mathrm{Tr}_{Q=0}(e^{-\beta
H})+O(e^{-\beta\lambda})}$ (6)

What we are interested in is the trace over the physical Hilbert space,
denoted in this notation as

$\braket{A}_{Q=1}=\frac{\mathrm{Tr}_{Q=1}(e^{-\beta
H}A)}{\mathrm{Tr}_{Q=1}(e^{-\beta H})}.$

Impurity quantities vanish in the absence of atomic states; therefore, if the
observable is the impurity Green’s function or the impurity free energy,
$\mathrm{Tr}_{Q=0}(e^{-\beta H}A)$ vanishes, hence the first term in the
numerator of Eq. (6) vanishes. If we consider as the observable charge $Q$, we
obtain:

$\braket{Q}=\frac{e^{-\beta\lambda}\mathrm{Tr}_{Q=1}(e^{-\beta
H})+O(e^{-2\beta\lambda})}{\mathrm{Tr}_{Q=0}(e^{-\beta
H})+O(e^{-\beta\lambda})}.$ (7)

The average over the physical Hilbert space of an observable, denoted as
$\braket{A}_{Q=1}$, can thus be expressed as the ratio of these two
quantities:

$\braket{A}_{Q=1}=\lim_{\lambda\rightarrow\infty}\frac{\braket{A}}{\braket{Q}}.$
(8)

For example, the electron single-particle Green’s function ${\cal G}$ must be
computed as

${\cal G}=\frac{1}{Q}\frac{\delta\Phi}{\delta\Delta}.$ (9)

Further details are provided in the next chapter.

####  I.2.1 Numerical treatment of projected quantities

In the next paragraph, we provide some details on the numerical implementation
of the projection, which requires special care due to the somewhat different
properties of pseudo-particle Green’s functions compared to regular electron
Green’s functions.

We begin with the imaginary axis quantities. The pseudo-particle Green’s
functions in imaginary time before projection are expressed by standard
formulas:

$\displaystyle
G_{p}(\tau)=\frac{1}{\beta}\sum_{i\omega}e^{-i\omega\tau}G_{p}(i\omega)=\int\frac{dx}{\pi}f(-x)e^{-x\tau}G^{\prime\prime}_{p}(x+i\delta).$
(10)

This expression is valid for fermionic pseudo-particles, but after projection,
the statistics does not matter. Additionally, it is valid only for positive
$\tau$. For negative $\tau$, the standard expression is:
$G_{p}(\tau)=-\int\frac{dx}{\pi}f(x)e^{-x\tau}G^{\prime\prime}_{p}(x+i\delta)$,
since this is well behaved at $\tau<0$, while the former expression is well
behaved for $\tau>0$. The projection requires shifting the variable $x$ by
$\lambda$ to the physical $Q=1$ part of the spectra, as discussed above. For
positive times $\tau>0$, we have:

$\displaystyle
G_{p}(\tau)=\lim_{\lambda\rightarrow\infty}\int\frac{dx}{\pi}f(-x-\lambda)e^{-x\tau-\lambda\tau}G^{\prime\prime}_{p}(x+\lambda+i\delta)=e^{-\lambda\tau}\int\frac{dx}{\pi}e^{-x\tau}G^{\prime\prime}_{p}(x+\lambda+i\delta)$
(11)

which is numerically challenging for large $\lambda$ but is finite at small
$\tau$. We note that the same variable shift reveals that for negative times
$G_{p}(\tau<0)\approx-e^{-\lambda\beta}\int\frac{dx}{\pi}e^{-x(\tau+\beta)}G_{p}^{\prime\prime}(x+\lambda+i\delta)$
indicating that the Green’s function vanishes for any time $\tau<0$. This
means that pseudo-particles do not propagate back in time, and consequently
the times on the backbone are time ordered. We also note in passing that
choosing bose statistics for pseudo-particle would require one to replace
$f(-x)$ in Eq.10 by $-n(-x)$, and when the limit $\lambda\rightarrow\infty$ is
taken they both give unity, hence statistics of pseudo-particles is
irrelevant.

To make numerics stable, we define numerically more appropriate pseudo-
particles quantities

$\displaystyle\widetilde{G}_{p}(\tau)\equiv e^{\lambda\tau}G_{p}(\tau).$ (12)

and we will recast all equations into the form that contain only
$\widetilde{G}_{p}(\tau)$, which from Eq. (11), is expressed as:

$\displaystyle\widetilde{G}_{p}(\tau)\equiv\int\frac{dx}{\pi}e^{-x\tau}G^{\prime\prime}_{p}(x+\lambda+i\delta)$
(13)

We will later prove order by order that $\widetilde{G}_{p}(\tau)$ is
numerically stable quantity.

On the real axis we also encounter a numerical ill-posed problem because at
negative values of $x$, the function $G_{p}^{\prime\prime}(x+\lambda+i\delta)$
becomes very small, of the order $e^{-|x|\beta}$. To compute physical
observables (or, in this case, $\widetilde{G}_{p}(\tau)$), we must multiply it
by a very large number $e^{|x|\tau}$.

To avoid this numerical instability, we devised two different strategies. One
commonly used approach in the literature is to employ
$G_{p}^{\prime\prime}(x+\lambda+i\delta)/f(-x)$. In other words, we rewrite
all equations in terms of:

$\widetilde{G}_{p}(x)\equiv G_{p}^{\prime\prime}(x+\lambda+i\delta)/f(-x).$
(14)

The second strategy, which we implement in Monte Carlo (MC) sampling, involves
storing $\log(G_{p}^{\prime\prime}(x+\lambda+i\delta))$ in addition to
$G_{p}(x+\lambda+i\delta)$. This is beneficial because, even though
$G_{p}^{\prime\prime}$ at negative $x$ is of the order $e^{-|x|\beta}$, its
logarithm is of the order $-|x|\beta$, which is easy to manage numerically.

We note in passing that $\widetilde{G}_{p}(\tau)$ and $\widetilde{G}_{p}(x)$
are interconnected. Their respective definitions in Eqs. (13) and (14)
establish the relationship:

$\displaystyle\widetilde{G}_{p}(\tau)=\int\frac{dx}{\pi}e^{-x\tau}f(-x)\widetilde{G}_{p}(x).$
(15)

This equation represents the analytic continuation kernel for fermionic
quantities, indicating that $\widetilde{G}_{p}(x)$ and
$\widetilde{G}_{p}(\tau)$ share the same spectral density of a fermionic-like
quantity in imaginary time or real frequency. Although the imaginary part of
the pseudo-particle spectra behaves as a fermionic quantity, this resemblance
is a result of our chosen numerical treatment. Pseudo-particles before
projection obey either fermionic or bosonic commutation relations (depending
on the number of electrons in the atomic state). After the projection, their
statistic does not matter, as both the bosonic and the fermionic pseudo-
particles have the same type of response functions. Additionally, it’s worth
noting that while the imaginary part of the pseudo Green’s functions, when
divided by the Fermi function, exhibits femion-like properties, the real part
does not. Namely, the real part of the pseudo Green’s functions is Kramers-
Kronig related to $G_{p}^{\prime\prime}(x+\lambda+i\delta)$:

$\displaystyle G_{p}^{\prime}(x)=-\frac{1}{\pi}\int
dy\frac{\widetilde{G}_{p}(y)f(-y)}{x-y}$ (16)

Finally we discuss the Dyson equation for pseudo-particles. It directly
follows from the effective action that the Dyson equation on the real axis is
given by:

$\displaystyle G_{m}(\omega)=\frac{1}{\omega-
E_{m}-\lambda-\Sigma_{m}(\omega)}.$ (17)

Here, $E_{m}$ represents the atomic energy levels, and $\Sigma_{m}$ is the
corresponding pseudo-particle self-energy. The projection requires shifting
the frequency $\omega\rightarrow\omega+\lambda$, resulting in:

$\displaystyle G_{m}(\omega+\lambda)=\frac{1}{\omega-
E_{m}-\Sigma_{m}(\omega+\lambda)}.$ (18)

This form is numerically easy to solve, with the real part computed directly.
The imaginary part, on the other hand, can be conveniently computed using:

$\displaystyle\widetilde{G}_{m}(\omega)=G^{\prime\prime}_{m}(\omega+\lambda)/f(-\omega)=\frac{\widetilde{\Sigma}_{m}(\omega)}{(\omega-
E_{m}-\Sigma^{\prime}(\omega+\lambda))^{2}+(\Sigma^{\prime\prime}(\omega+\lambda))^{2}}$
(19)

Solving the Dyson equation on the imaginary axis is more challenging. We first
rewrite the Dyson Eq. 17 as:

$G=G^{0}+G^{0}\Sigma G.$ (20)

Here, the product should be understood as the integral over time, and
$G^{0}_{m}(\tau)=e^{-(E_{m}+\lambda)\tau}$. It’s important to note that pseudo
Green’s functions do not allow negative times, as proven above. Hence, we
have:

$\displaystyle
G_{m}(\tau)=e^{-(E_{m}+\lambda)\tau}+\int_{0}^{\tau}d\tau_{2}\int_{0}^{\tau_{2}}d\tau_{1}e^{-(E_{m}+\lambda)(\tau-\tau_{2})}\Sigma_{m}(\tau_{2}-\tau_{1})G_{m}(\tau_{1})$
(21)

This expression is numerically challenging to compute when $\lambda$ is
large.However, we introduced above more convenient functions $\widetilde{G}$,
in terms of which we can write

$\displaystyle\widetilde{G}_{m}(\tau)=e^{-E_{m}\tau}+\int_{0}^{\tau}d\tau_{2}\int_{0}^{\tau_{2}}d\tau_{1}e^{-E_{m}(\tau-\tau_{2})}\widetilde{\Sigma}_{m}(\tau_{2}-\tau_{1})\widetilde{G}_{m}(\tau_{1}).$
(22)

Here we used $G_{m}(\tau)=e^{-\lambda\tau}\widetilde{G}_{m}(\tau)$ and
$\Sigma_{m}(\tau)=e^{-\lambda\tau}\widetilde{\Sigma}_{m}(\tau)$. Eq. (22) can
be solved for $\widetilde{G}_{m}$ by matrix inversion on a discrete mesh or by
iteratively evaluating the integral.

###  I.3 Method to calculate the electron self-energy

The electron self-energy is easier to compute from the two particle response
function than from the Dyson equation, as first discovered by Bulla in the
context of NRG [27]. This is particularly crucial here, as the convergence
with perturbation order is significantly faster when the self-energy is
computed from the two-particle response function. In this case, even the first
order (non-crossing approximation method) yields the correct high-frequency
behavior. Of course, at high order, both the Dyson equation and the two-
particle response function can be used to compute the self-energy, and by
comparing the two results, we can also test the convergence.

For generic impurity problem the connection between the electron self-energy
$\Sigma(\omega)$ and two particle response function is given by:

$\displaystyle({{\cal
G}}(i\omega)\Sigma(i\omega))_{\beta\alpha}=-\sum_{ijk}\int_{0}^{\beta}d\tau
e^{i\omega\tau}\braket{T_{\tau}\psi_{\beta}(\tau)\psi^{\dagger}_{i}(0)\psi^{\dagger}_{j}(0)\psi_{k}(0^{-})}\frac{1}{2}(U_{ijk\alpha}-U_{jik\alpha}).$
(23)

Here ${\cal G}$ is the electron single-particle Green’s function (not to be
confused with pseudo-particle Green’s function $G_{p}$), and the interaction
$\hat{U}$ has the form:
$\hat{U}=\sum_{ijk\alpha}U_{ijk\alpha}\psi_{i}^{\dagger}\psi_{j}^{\dagger}\psi_{k}\psi_{\alpha}$.
$\Sigma$ here is the electron self-energy. This expression can be obtained
either diagrammatically or using the equation of motion. Before demonstrating
how to compute such a quantity within the current pseudo-particle formulation,
we first want to show how the single-particle Green’s function ${\cal G}$ is
computed, and how it relates to pseudo-particles.

Let’s start by revisiting how to compute the electron single-particle Green’s
function ${\cal G}$ for a generic impurity problem. The impurity partition
function is given by:

$\displaystyle Z=\int{\cal
D}[\psi,\psi^{+}]e^{-S_{atm}-\int_{0}^{\beta}\int_{0}^{\beta}d\tau
d\tau^{\prime}\sum_{\alpha,\beta}\psi^{\dagger}_{\alpha}(\tau)\Delta_{\alpha,\beta}(\tau-\tau^{\prime})\psi_{\beta}(\tau^{\prime})},$
(24)

hence we have

$\displaystyle-\frac{\delta\log
Z}{\delta\Delta_{\alpha\beta}(\tau-\tau^{\prime})}=\frac{1}{Z}\int{\cal
D}[\psi,\psi^{+}]e^{-S}\psi^{\dagger}_{\alpha}(\tau)\psi_{\beta}(\tau^{\prime}),$
(25)

which is the expression for the electron Green’s function. Therefore we have

$\displaystyle{{\cal G}}_{\beta\alpha}(\tau^{\prime}-\tau)=-\frac{\delta\log
Z}{\delta\Delta_{\alpha\beta}(\tau-\tau^{\prime})}.$ (26)

Eq. 26 is used in the conventional CTQMC (bare expansion). In the bold
version, BHQMC, we dress the pseudo-particles with their self-energy, and
hence the quantity to sample is not $Z$ but rather the Luttinger-Ward
functional $\Phi_{\Delta}[\\{G_{p}\\}]$ in terms of which all self-energies
are obtained as $\Sigma_{p}=\frac{\delta\Phi_{\Delta}[\\{G_{p}\\}]}{\delta
G_{p}}$. This is related to $\log Z$ (before the projection) by the standard
Klein functional:

$\displaystyle-\log(Z)=\beta\mathrm{Tr}\log(G_{p})-\beta\mathrm{Tr}(({G_{p}^{0}}^{-1}-{G_{p}}^{-1})G_{p})+\Phi_{\Delta}[\\{G_{p}\\}]$
(27)

Here $G_{p}$ and $G_{p}^{0}$ are the pseudo-particle Green’s functions. The
Klein functional 27 is stationary with respect to $G_{p}$, hence
$\delta\log(Z)/\delta G_{p}=0$. The derivative with respect to other
quantities, like $\Delta$ is thus

$\frac{\delta\log(Z)}{\delta\Delta}=\frac{\delta
G_{p}}{\delta\Delta}\frac{\delta\log(Z)}{\delta
G_{p}}+\left(\frac{\partial\log(Z)}{\partial\Delta}\right)_{G_{p}}=-\left(\frac{\partial\Phi_{\Delta}}{\partial\Delta}\right)_{G_{p}}$
(28)

The first term vanishes because $\delta\log(Z)/\delta G_{p}=0$, and since
$G_{p}^{0}$ does not depend on $\Delta$, we just need to take the derivative
of $\Phi$ functional.

Notice that $\Phi_{\Delta}[\\{G_{p}\\}]$ depends on $\Delta$ in a similar way
as $\Phi_{U}[\\{{\cal G}\\}]$ depends on interaction $U$ in the in weak
coupling theory: the theory is not dressed in this channel. The key difference
with the weak coupling is that in the weak coupling the interaction can be
dressed with the many-body fluctuations (we can use dressed $W$ instead of $U$
to develop a theory), while here any dressing of $\Delta$ by the Dyson
equation vanishes exactly because of the projection. This is because
projection requires a single loop of pseudo-particles.

Because the theory is not dressed in $\Delta$, it is obvious from Eqs. 28 and
26 that

${{\cal
G}}_{\beta\alpha}(\tau^{\prime}-\tau)=\frac{\partial\Phi_{\Delta}(\\{G_{p}\\})}{\partial\Delta_{\alpha\beta}(\tau-\tau^{\prime})},$
(29)

where we cut any hybridization lines, but do not cut $G_{p}$ (treat $G_{p}$ as
a constant).

However, this form is valid in the grand-canonical ensemble (before the
projection), and as discussed in chapter Projection to the physical Hilbert
space, we need to divide this value by $\braket{Q}$ to account for the exact
projection. Therefore,

${{\cal
G}}_{\beta\alpha}(\tau^{\prime}-\tau)=\frac{1}{\braket{Q}}\frac{\partial\Phi_{\Delta}(\\{G_{p}\\})}{\partial\Delta_{\alpha\beta}(\tau-\tau^{\prime})}.$
(30)

In the MC sampling we cut each hybridization line, and we obtain $n$
contribution to the singe-particle Green’s function, where $n$ is the
perturbation order. Fig. 5 shows a generic contribution to $\Phi_{\Delta}$.

Figure 5:  The sketch of the Luttinger-Ward functional
$\Phi_{\Delta}[\\{G_{p}\\}]$ where we concentrate on a single hybridization
line $\Delta$.

In the pseudo-particle formulation, the single-particle electron Green’s
function takes the form of the two particle correlation function:

$\displaystyle{\cal
G}_{\beta\alpha}(\tau^{\prime}-\tau)=-\sum_{\\{m\\}}\braket{m_{4}}{\psi_{\beta}}{m_{3}}\braket{m_{2}}{\psi^{\dagger}_{\alpha}}{m_{1}}\braket{T_{\tau}a^{\dagger}_{m_{4}}(\tau^{\prime})a_{m_{3}}(\tau^{\prime})a^{\dagger}_{m_{2}}(\tau)a_{m_{1}}(\tau)}$
(31)

This can be easily derived by noting that the action of the annihilation
operator on the atomic states is

$\displaystyle\psi_{\beta}=\sum_{\\{m\\}}\ket{m_{4}}\braket{m_{4}}{\psi_{\beta}}{m_{3}}\bra{m_{3}}\rightarrow\sum_{\\{m\\}}a_{m_{4}}^{\dagger}\braket{m_{4}}{\psi_{\beta}}{m_{3}}a_{m_{3}}$
(32)

and similarly

$\displaystyle\psi^{\dagger}_{\alpha}=\sum_{\\{m\\}}\ket{m_{2}}\braket{m_{2}}{\psi^{\dagger}_{\alpha}}{m_{1}}\bra{m_{1}}\rightarrow\sum_{\\{m\\}}a_{m_{2}}^{\dagger}\braket{m_{2}}{\psi^{\dagger}_{\alpha}}{m_{1}}a_{m_{1}}$
(33)

We then recognize that each contribution to ${\cal G}$ can be assigned a set
of atomic states $m_{1},m_{2},m_{3},m_{4}$, in which $m_{1}$ and $m_{2}$
appear just before and after the action of $\psi^{\dagger}_{\alpha}$ operator,
and $m_{3}$ and $m_{4}$ just before and after $\psi_{\beta}$ operator, as
sketched in Fig. 5. We will denote such contribution with
$\left(\frac{\partial\Phi}{\partial\Delta_{\alpha\beta}(\tau-\tau^{\prime})}\right)_{m_{1},m_{2}\in\tau;m_{3},m_{4}\in\tau^{\prime}}$
so that

${\cal
G}_{\beta\alpha}(\tau^{\prime}-\tau)=\sum_{\\{m\\}}\left(\frac{\partial\Phi}{\partial\Delta_{\alpha\beta}(\tau-\tau^{\prime})}\right)_{m_{1},m_{2}\in\tau;m_{3},m_{4}\in\tau^{\prime}}.$
(34)

This means that

$\left(\frac{\partial\Phi}{\partial\Delta_{\alpha\beta}(\tau-\tau^{\prime})}\right)_{m_{1},m_{2}\in\tau;m_{3},m_{4}\in\tau^{\prime}}=-\braket{m_{4}}{\psi_{\beta}}{m_{3}}\braket{m_{2}}{\psi^{\dagger}_{\alpha}}{m_{1}}\braket{T_{\tau}a^{\dagger}_{m_{4}}(\tau^{\prime})a_{m_{3}}(\tau^{\prime})a^{\dagger}_{m_{2}}(\tau)a_{m_{1}}(\tau)}$
(35)

Now we are in a position to rewrite Eq. 23 in pseudo-particle formulation. We
first define the imaginary axis equivalent ${\cal
G}(i\omega)\Sigma(i\omega)=\int_{0}^{\beta}d\tau e^{i\omega\tau}({\cal
G}\cdot\Sigma)(\tau)$ and than the imaginary time quantity is

$\displaystyle({\cal
G}\cdot\Sigma)_{\beta\alpha}(\tau^{\prime}-\tau)=-\sum_{\\{m\\}}\braket{m_{4}}{\psi_{\beta}}{m_{3}}F^{\alpha}_{m_{2},m_{1}}\braket{T_{\tau}a^{\dagger}_{m_{4}}(\tau^{\prime})a_{m_{3}}(\tau^{\prime})a^{\dagger}_{m_{2}}(\tau)a_{m_{1}}(\tau)}$
$\displaystyle
F^{\alpha}_{m_{2},m_{1}}=\sum_{ijk}\frac{1}{2}(U_{ijk\alpha}-U_{jik\alpha})\braket{m_{2}}{\psi^{\dagger}_{i}\psi^{\dagger}_{j}\psi_{k}}{m_{1}}.$
(36)

This can be derived by noting that

$\displaystyle\psi^{\dagger}_{i}\psi^{\dagger}_{j}\psi_{k}=\sum_{\\{m\\}}\ket{m_{2}}\braket{m_{2}}{\psi^{\dagger}_{i}\psi^{\dagger}_{j}\psi_{k}}{m_{1}}\bra{m_{1}}\rightarrow\sum_{\\{m\\}}a_{m_{2}}^{\dagger}\braket{m_{2}}{\psi^{\dagger}_{i}\psi^{\dagger}_{j}\psi_{k}}{m_{1}}a_{m_{1}}$
(37)

We than see that ${\cal G}\cdot\Sigma$ in Eq. 36 requires the same correlation
function as ${\cal G}$ in Eq. 31, except that the prefactors are different. We
have

$\displaystyle({\cal
G}\cdot\Sigma)_{\beta\alpha}(\tau^{\prime}-\tau)=\sum_{\\{m\\}}\frac{F^{\alpha}_{m_{2},m_{1}}}{\braket{m_{2}}{\psi^{\dagger}_{\alpha}}{m_{1}}}\left(\frac{\partial\Phi}{\partial\Delta_{\alpha\beta}(\tau-\tau^{\prime})}\right)_{m_{1},m_{2}\in\tau;m_{3},m_{4}\in\tau^{\prime}}.$
(38)

These matrix elements
$F^{\alpha}_{m_{2},m_{1}}/\braket{m_{2}}{\psi^{\dagger}_{\alpha}}{m_{1}}$ can
be precomputed, and hence the contributions to $({\cal G}\cdot\Sigma)$ is
easily computed along with contributions to ${\cal G}$ with no extra cost.
Once the sampling is concluded, we perform the Fourier transform of $({\cal
G}\cdot\Sigma)(\tau^{\prime}-\tau)$ to get ${\cal G}(i\omega)\Sigma(i\omega)$.
We compute ${\cal G}(i\omega)$ in the same process, hence $\Sigma(i\omega)$
can be obtained by dividing the two quantities.

On the real axis the equation has a very similar form, namely,

$\displaystyle({\cal G}\cdot\Sigma)_{\beta\alpha}(\omega)$ $\displaystyle=$
$\displaystyle\sum_{\\{m\\}}\frac{F^{\alpha}_{m_{2},m_{1}}}{\braket{m_{2}}{\psi^{\dagger}_{\alpha}}{m_{1}}}\widetilde{\left(\frac{\partial\Phi}{\partial\Delta_{\alpha\beta}(\omega)}\right)}_{m_{1},m_{2}\in\tau;m_{3},m_{4}\in\tau^{\prime}}$
(39) $\displaystyle{\cal G}_{\beta\alpha}(\omega)$ $\displaystyle=$
$\displaystyle\sum_{\\{m\\}}\widetilde{\left(\frac{\partial\Phi}{\partial\Delta_{\alpha\beta}(\omega)}\right)}_{m_{1},m_{2}\in\tau;m_{3},m_{4}\in\tau^{\prime}}$
(40)

however
$\widetilde{\left(\frac{\partial\Phi}{\partial\Delta_{\alpha\beta}(\omega)}\right)}$
is not just the functional derivative, but needs to remove also corresponding
fermi function from the result, as discussed in chapter Feynman rules on
imaginary axis.

For completness, we also state the equations for the two particle response
functions, i.e., susceptibilities. We start by taking the second derivative of
the partition function Eq. 24:

$\displaystyle\chi_{\alpha\beta\gamma\delta}(\tau_{1},\tau_{2},\tau_{3},\tau_{4})\equiv-\braket{T_{\tau}\psi^{\dagger}_{\alpha}(\tau_{1})\psi_{\beta}(\tau_{2})\psi^{\dagger}_{\gamma}(\tau_{3})\psi_{\delta}(\tau_{4})}=-\frac{\delta^{2}\log
Z}{\delta\Delta_{\alpha\beta}(\tau_{1}-\tau_{2})\delta\Delta_{\gamma\delta}(\tau_{3}-\tau_{4})}.$
(41)

This is what is used in the bare expansion (conventional CTQMC). But here we
have self-consistent propagators and we want to express them in terms of
$\Phi_{\Delta}$ and $G_{p}$. We realize that the Dyson equation is always
satisfied, therefore, according to Eq. 27 we can take the derivative of
$\Phi_{\Delta}$ functional only

$\displaystyle\chi_{\alpha\beta\gamma\delta}(\tau_{1},\tau_{2},\tau_{3},\tau_{4})$
$\displaystyle=\frac{\delta}{\delta\Delta_{\alpha\beta}(\tau_{1}-\tau_{2})}\left(\frac{\partial\Phi_{\Delta}}{\partial\Delta_{\gamma\delta}(\tau_{3}-\tau_{4})}\right)_{G}$
$\displaystyle=\left(\frac{\partial^{2}\Phi_{\Delta}}{\partial\Delta_{\alpha\beta}(\tau_{1}-\tau_{2})\partial\Delta_{\gamma\delta}(\tau_{3}-\tau_{4})}\right)_{G}+\sum_{i}\left(\frac{\partial^{2}\Phi_{\Delta}}{\partial\Delta_{\gamma\delta}(\tau_{3}-\tau_{4})\partial
G_{i}(\tau_{i+1}-\tau_{i})}\right)\frac{\delta
G_{i}(\tau_{i+1}-\tau_{i})}{\delta\Delta_{\alpha\beta}(\tau_{1}-\tau_{2})}$

Here the last equation shows that we might cut two hybridizations, or, just
one hybridization and one pseudo-particle propagator. The latter requires the
mixed derivative, which plays the central role in this formalism. The quantity
that appears the last in Eq. I.3 can be expressed as follows:

$\displaystyle\frac{\delta
G_{i}(\tau_{i+1}-\tau_{i})}{\delta\Delta_{\alpha\beta}(\tau_{1}-\tau_{2})}=\int_{\tau_{i}}^{\tau_{i+1}}d\tau_{2}\int_{\tau_{i}}^{\tau_{2}^{\prime}}d\tau_{1}G_{i}(\tau_{i+1}-\tau_{2}^{\prime})\frac{\partial\Sigma_{i}(\tau_{2}^{\prime}-\tau_{1}^{\prime})}{\partial\Delta_{\alpha\beta}(\tau_{1}-\tau_{2})}G_{i}(\tau_{1}^{\prime}-\tau_{i}).$
(43)

This follows from taking the derivative of $GG^{-1}=1$ and noticing that
$\delta G^{-1}/\delta\Delta=-\delta\Sigma/\delta\Delta$, because
$G^{-1}(\omega)=\omega-E-\Sigma$. In grand canonical ensemble one also
generates higer order terms in expansion, but those will require more than one
pseudo-loop, and hence vanish after projection. Now we realize that the self-
energy $\Sigma=\partial\Phi_{\Delta}/\partial G$, hence its derivative with
respect to $\Delta$ is again the mixed derivative, which we introduced before:

$\displaystyle\frac{\delta
G_{i}(\tau_{i+1}-\tau_{i})}{\delta\Delta_{\alpha\beta}(\tau_{1}-\tau_{2})}=\int_{\tau_{i}}^{\tau_{i+1}}d\tau_{2}\int_{\tau_{i}}^{\tau_{2}^{\prime}}d\tau_{1}G_{i}(\tau_{i+1}-\tau_{2}^{\prime})\frac{\partial^{2}\Phi_{\Delta}}{\partial\Delta_{\alpha\beta}(\tau_{1}-\tau_{2})\partial
G_{i}(\tau_{2}^{\prime}-\tau_{1}^{\prime})}G_{i}(\tau_{1}^{\prime}-\tau_{i}).$
(44)

We note that here $\Phi_{\Delta}$ needs to contain all Feynman diagrams, not
just the diagram that we currently concentrate on. To evaluate all such
contributions, we can sample two quantities:

$\displaystyle\chi^{(1)}(\tau_{1},\tau_{2},\tau_{3},\tau_{4})=\left(\frac{\partial^{2}\Phi_{\Delta}}{\partial\Delta_{\alpha\beta}(\tau_{1}-\tau_{2})\partial\Delta_{\gamma\delta}(\tau_{3}-\tau_{4})}\right)_{G}$
(45)
$\displaystyle\Gamma_{\gamma\delta,p}(\tau_{p+1},\tau_{p},\tau_{3},\tau_{4})=\left(\frac{\partial^{2}\Phi_{\Delta}}{\partial\Delta_{\gamma\delta}(\tau_{3}-\tau_{4})\partial
G_{p}(\tau_{p+1}-\tau_{p})}\right)$ (46)

And with postprocessing, we can compute $\chi=\chi^{(1)}+\chi^{(2)}$, where
the second part $\chi^{(2)}$ is:

$\displaystyle\chi^{(2)}_{\alpha\beta\gamma\delta}(\tau_{1},\tau_{2},\tau_{3},\tau_{4})=\sum_{p}\int_{0}^{\beta}d\tau_{4}\int_{0}^{\tau_{4}}d\tau_{3}\int_{0}^{\tau_{3}}d\tau_{2}\int_{0}^{\tau_{2}}d\tau_{1}\Gamma_{\alpha\beta,p}(\tau_{2}^{\prime},\tau_{1}^{\prime},\tau_{1},\tau_{2})G_{p}(\beta+\tau_{1}^{\prime}-\tau_{4}^{\prime})\Gamma_{\gamma\delta,p}(\tau_{4}^{\prime},\tau_{3}^{\prime},\tau_{3},\tau_{4})G_{p}(\tau_{3}^{\prime}-\tau_{2}^{\prime})$

Notice that $\partial^{2}\Phi_{\Delta}$ appears twice in $\chi^{(2)}$ term.
Also notice that at the lowest order, within Non-crossing approximation, the
first order $\chi^{(1)}$ vanishes, but the second part $\chi^{(2)}$ is finite.

###  I.4 Feynman rules on imaginary axis

Figure 6:  The sketch of the Luttinger-Ward functional
$\Phi_{\Delta}[\\{G_{p}\\}]$.

In imaginary time, the Feynman rules are standard. However, the pseudo-
particles $G_{p}(\tau)$ vanish for negative times. Consequently, the time
variables on the backbone increase (appear sorted) and are distributed in the
interval $[0,\beta)$. The hybridization functions have a fermionic nature,
allowing negative times for hybridization functions, and they exhibit
antiperiodic behavior in imaginary time. To make discussion easier, we show in
Fig. 6 an example of a typical third-order Feynman diagram for the Luttinger-
Ward functional. Any self-energy contribution is obtained by cutting the
respective propagator. Without loos of generality, we can choose the first
time to be $0$, hence the expression for this particular diagram is

$\displaystyle\Phi_{\Delta}[\\{G_{p}\\}]=\prod_{i}\int
d\tau_{i}G_{0}(\tau_{1})G_{1}(\tau_{2}-\tau_{1})G_{2}(\tau_{3}-\tau_{2})G_{3}(\tau_{4}-\tau_{3})G_{4}(\tau_{5}-\tau_{4})G_{5}(\beta-\tau_{5})\Delta_{1}(\tau_{2})\Delta_{3}(\tau_{5}-\tau_{3})\Delta_{2}(\tau_{4}-\tau_{1}),$

where the integral spans the space
$0\leq\tau_{1}\leq\tau_{2}\cdots\leq\tau_{5}\leq\beta$.

The overall sign could be computed from the number of fermionic loops,
following conventional Feynman rules. However, this process is tedious because
choosing what constitutes a fermionic loop in these diagrams is not unique.
Alternatively, it could be computed by the same method as in conventional
CTQMC. There are two contributions to the sign: i) the parity of permutation
of hybridizations, and ii) the matrix elements of the $\psi^{\dagger}$ and
$\psi$ operators. The parity of permutation for hybridization lines is easy to
compute by inspection. In the case of Fig. 6 the permutation is
$(0,2)(1,4)(3,5)$. Here, the numbers in the bracket denotes the start and the
end of hybridization line. Since this is an even permutation of the sorted
list $(0,1)(2,3)(4,5)$, the sign contribution (i) is positive. The matrix
elements of $\psi$ operators in Fig. 6 are
$\braket{m_{5}}{\psi_{5}^{\dagger}}{m_{4}}\braket{m_{4}}{\psi_{4}^{\dagger}}{m_{3}}\braket{m_{3}}{\psi_{3}}{m_{2}}\braket{m_{2}}{\psi_{2}^{\dagger}}{m_{1}}\braket{m_{1}}{\psi_{1}}{m_{0}}\braket{m_{0}}{\psi_{0}}{m_{5}}$.
In a single orbital model, only four states of the atom are possible: the
empty and doubly occupied states, as well as $\ket{\uparrow}$ and
$\ket{\downarrow}$. If we choose hybridizations $\Delta_{1}$ and $\Delta_{3}$
to correspond to $\Delta_{\downarrow}$ and $\Delta_{2}$ to $\Delta_{\uparrow}$
spin, we must have atomic states $\ket{1}$ and $\ket{3}$ to be the empty
impurity state, $\ket{5}$ the doubly occupied, and $\ket{0}$ and $\ket{4}$ to
be $\ket{\uparrow}$, while $\ket{2}$ is $\ket{\downarrow}$. The product of
matrix elements in this case is also $+1$, resulting in an overall sign of
$+1$.

To compute the self-energy for the pseudo-particles involved in this
functional we need to take the functional derivative in a standard way:

$\Sigma_{m}(\beta-\tau)=\frac{\delta\Phi_{\Delta}[\\{G_{p}\\}]}{\delta
G_{m}(\tau)}$ (47)

It’s important to note that only positive times are allowed, hence
$\Sigma_{m}(\beta-\tau)$ appears instead of $\Sigma_{m}(-\tau)$.

Next we discuss how to use numerically stable quantities to evaluate all self-
energies. First, all $G_{p}(\tau)$ should be substituted by
$G_{p}(\tau)=e^{-\lambda\tau}\widetilde{G}_{p}(\tau)$, where
$\widetilde{G}_{p}$ are numerically stable. It is easy to see that this
substitution always gives an overal factor of $e^{-\beta\lambda}$ to
$\Phi_{\Delta}[\\{\widetilde{G}\\}]$, while otherwise the form of
$\Phi_{\Delta}$ remains the same. We will hence define
$\widetilde{\Phi}_{\Delta}\equiv\Phi_{\Delta}/e^{-\beta\lambda}$ so that
$\lambda$ disappears in expression for
$\widetilde{\Phi}_{\Delta}[\\{\widetilde{G}\\}]$. The pseudo self-energies and
the pseudo-green’s functions also need to be substituted by
$\Sigma_{p}(\tau)=e^{-\lambda\tau}\widetilde{\Sigma}_{p}(\tau)$. We recognize
that in Eq. (47) these exponential factors precisely cancel, and hence we can
use

$\widetilde{\Sigma}_{m}(\beta-\tau)=\frac{\delta\widetilde{\Phi}_{\Delta}[\\{\widetilde{G}_{p}\\}]}{\delta\widetilde{G}_{m}(\tau)}$
(48)

where we removed $\lambda$ from all equations and the projected quantities are
now numerically stable.

Finally, the contribution to the electron single particle Green’s function
${\cal G}$, derived above, is

${\cal
G}_{\beta\alpha}(-\tau)=\frac{1}{\braket{Q}}\frac{\partial\Phi_{\Delta}[\\{G_{p}\\}]}{\partial\Delta_{\alpha\beta}(\tau)},$
(49)

which also requires some adjustments for numeric stability. As discussed
earlier $\Phi_{\Delta}$ carries a factor of $e^{-\beta\lambda}$ due to the
projection to the $Q=1$ subspace. This factor is canceled by $\braket{Q}$,
which is also of the order of $e^{-\beta\lambda}$. More precisely, the
occupation of any pseudo-particle is
$\braket{a_{p}^{\dagger}a_{p}}=G_{p}(\beta^{-})$, hence

$\braket{Q}=e^{-\beta\lambda}\sum_{p}\widetilde{G}_{p}(\beta).$ (50)

We can thus define $\braket{\widetilde{Q}}=\braket{Q}/e^{-\beta\lambda}$, and
in terms of these quantities, we can write

${\cal
G}_{\beta\alpha}(-\tau)=\frac{1}{\braket{\widetilde{Q}}}\frac{\partial\widetilde{\Phi}_{\Delta}[\\{\widetilde{G}_{p}\\}]}{\partial\Delta_{\alpha\beta}(\tau)}$
(51)

where

$\displaystyle\widetilde{Q}=\sum_{p}\widetilde{G}_{p}(\beta).$ (52)

We thus managed to rewrite all equations in terms of the numerically stable
quantities. For the above example we have

$\displaystyle\widetilde{\Phi}_{\Delta}[\\{\widetilde{G}_{p}\\}]=\prod_{i}\int
d\tau_{i}\widetilde{G}_{0}(\tau_{1})\widetilde{G}_{1}(\tau_{2}-\tau_{1})\widetilde{G}_{2}(\tau_{3}-\tau_{2})\widetilde{G}_{3}(\tau_{4}-\tau_{3})\widetilde{G}_{4}(\tau_{5}-\tau_{4})\widetilde{G}_{5}(\beta-\tau_{5})\Delta_{1}(\tau_{2})\Delta_{3}(\tau_{5}-\tau_{3})\Delta_{2}(\tau_{4}-\tau_{1})$

and for example

$\displaystyle\widetilde{\Sigma}_{2}(\beta-\tau_{3}+\tau_{2})=\prod_{i\neq
2,3}\int
d\tau_{i}\widetilde{G}_{0}(\tau_{1})\widetilde{G}_{1}(\tau_{2}-\tau_{1})\widetilde{G}_{3}(\tau_{4}-\tau_{3})\widetilde{G}_{4}(\tau_{5}-\tau_{4})\widetilde{G}_{5}(\beta-\tau_{5})\Delta_{1}(\tau_{2})\Delta_{3}(\tau_{5}-\tau_{3})\Delta_{2}(\tau_{4}-\tau_{1})$
$\displaystyle{\cal
G}_{3}(\tau_{3}-\tau_{5})=\frac{1}{\braket{\widetilde{Q}}}\prod_{i\neq
3,5}\int
d\tau_{i}\widetilde{G}_{0}(\tau_{1})\widetilde{G}_{1}(\tau_{2}-\tau_{1})\widetilde{G}_{2}(\tau_{3}-\tau_{2})\widetilde{G}_{3}(\tau_{4}-\tau_{3})\widetilde{G}_{4}(\tau_{5}-\tau_{4})\widetilde{G}_{5}(\beta-\tau_{5})\Delta_{1}(\tau_{2})\Delta_{2}(\tau_{4}-\tau_{1})$

Notice that all quantities in these expressions are numerically stable, and
Feynman rules are the same as for the original quantities without tilde.

###  I.5 Feynman rules on the real axis

From the expression on the imaginary axis for Feynman diagrams, one can always
find their counterparts on the real axis by replacing the Matsubara summations
with integrals over the complex plane. However, each Matsubara summation in
this step leads to several terms, as many as the number of propagators in the
loop. For a diagram of order $n$, where $n$ successive summations are
necessary, and each loop involves at least $2$ propagators, this results in
generating at least $2^{n}$ terms on the real axis. Unfortunately, this often
makes Monte Carlo (MC) sampling on the real axis exponentially more expensive
than on the imaginary axis.

As explained in the main text, the strong coupling expansion method offers an
important simplification: it generates precisely one term for any Matsubara
sum over the hybridization frequency $i\nu$, and exactly $2n$ terms are
generated when evaluating the sum over $i\omega$. To compute the pseudo
particle self-energy we need to evaluate $n$ Matsubara sums over the
hybridization frequencies $i\nu$, while $i\omega$ becomes the external
frequency. As a result, we generate a single term. On the other hand, the
electron Green’s function requires $n-1$ Matsubara sums over hybridization
frequencies and also the sum over $i\omega$. The former gives a single term,
while the latter generates $2n$ terms. Hence, we end up with $2n$ terms, which
we will derive below. In contrast, in the weak coupling expansion, at least
$2^{n}$ terms are generated.

We will first derive expressions for the pseudo particle self-energies. In
this case we need to consider only the integrals over the branch-cuts of
hybridization functions, which enormously simplifies the Matsubara summations.
In the main part of the text, we explained that the Matsubara sum over
hybridization frequency $i\nu_{n}$ results in a simple integral over the real
variable $y$:

$\displaystyle T\sum_{i\nu_{n}}\Delta(\pm
i\nu_{n})\prod_{j=1}^{N}G_{j}(i\nu_{n}+i\omega+\sum_{\alpha}y_{\alpha}a_{\alpha}^{j}+\sum_{\beta}i\nu_{m}^{\beta}a_{\beta}^{j})=$
$\displaystyle\pm\int dyA_{c}(\pm
y)f(y)\prod_{j=1}^{N}G_{j}(y+i\omega+\sum_{\alpha}y_{\alpha}a_{\alpha}^{j}+\sum_{\beta}i\nu_{m}^{\beta}a_{\beta}^{j}).$
(53)

This expression is derived by integration over the complex plane of the
following integral:

$\displaystyle T\sum_{i\nu_{n}}\Delta(\pm
i\nu_{n})\prod_{j=1}^{N}G_{j}(i\nu_{n}+i\omega+\sum_{\alpha}y_{\alpha}a_{\alpha}^{j}+\sum_{\beta}i\nu_{m}^{\beta}a_{\beta}^{j})=$
$\displaystyle-\oint\frac{dz}{2\pi i}\Delta(\pm
z)f(z)\prod_{j=1}^{N}G_{j}(z+i\omega+\sum_{\alpha}y_{\alpha}a_{\alpha}^{j}+\sum_{\beta}i\nu_{m}^{\beta}a_{\beta}^{j}).$
(54)

where the integral must encircle all Matsubara points but has to avoid all
branch-cuts of functions $\Delta$ and $G_{j}$. As explained in the main part
of the text, the branch-cut on the real axis (at $z=y\in\mathbb{R}$) gives Eq.
53, while we argued that all other branch cuts of $G_{j}$ give vanishing
contributions after projection. This is because integral around the branch-cut
of $G_{j}$ for $j=p$ requires the substitution of $z\rightarrow y\pm
i\delta-i\omega-\sum_{\alpha}y_{\alpha}a_{\alpha}^{p}-\sum_{\beta}i\nu_{m}^{\beta}a_{\beta}^{p}$,
and we get

$\int
dyA_{p}(y)f(y-\sum_{\alpha}a_{\alpha}^{p}y_{\alpha}-\cdots)\Delta(\pm(y-i\omega-\sum_{\alpha}y_{\alpha}a_{\alpha}^{p}-\sum_{\beta}i\nu_{m}^{\beta}a_{\beta}^{p}))\prod_{j=1,j\neq
p}^{N}G_{j}(y+\sum_{\alpha}y_{\alpha}(a_{\alpha}^{j}-a_{\alpha}^{p})+\sum_{\beta}i\nu_{m}^{\beta}(a_{\beta}^{j}-a_{\beta}^{p}))$

Here $\cdots$ stand for a sum of Matsubara frequencies, which can result in a
bosonic or fermionic frequency. If bosonic, we can ignore it in the Fermi
function, as $f(x+i\Omega)=f(x)$ ($i\Omega=2n\pi T$). If fermionic, it gives
$f(x+i\Omega)=-n(x)$. But we also apply projection, which requires us to shift
the argument of $A_{p}(y)$ to $A_{p}(y+\lambda)$ (with $\lambda$ large), so
$f(y+\cdots)$ will also shift to
$f(y+\lambda+\cdots)=e^{-\beta(y+\lambda+\cdots)}$. Hence, irrespective of the
nature of the Matsubara sum denoted by $\cdots$, we get the same exponent
$e^{-\beta(y+\lambda+\cdots)}$, with additional minus sign in the prefactor
when the Matsubara sum is fermionic. After projection we hence have

$\displaystyle\pm\int
dyA_{p}(y+\lambda)e^{-\beta(y+\lambda-\sum_{\alpha}a_{\alpha}^{p}y_{\alpha})}\Delta(\pm(y+\lambda-i\omega-\sum_{\alpha}y_{\alpha}a_{\alpha}^{p}-\sum_{\beta}i\nu_{m}^{\beta}a_{\beta}^{p}))\times$
$\displaystyle\times\prod_{j=1,j\neq
p}^{N}G_{j}(y+\lambda+\sum_{\alpha}y_{\alpha}(a_{\alpha}^{j}-a_{\alpha}^{p})+\sum_{\beta}i\nu_{m}^{\beta}(a_{\beta}^{j}-a_{\beta}^{p}))$
(55)

The crucial observation is that all terms in the integral are non-divergent,
while the exponent contains a factor $e^{-\beta\lambda}$. Since $\lambda$ is
set to infinity, the entire term vanishes after projection.

In Eq. 55 there are also terms like $e^{\beta a_{\alpha}^{p}y_{\alpha}}$ which
become large for large positive $y_{\alpha}$. However, a more careful
derivation shows that for each such term, we also have corresponding Fermi
function (from the already performed Matsubara summation), and the combined
$f(a_{\alpha}^{p}y_{\alpha})e^{\beta a_{\alpha}^{p}y_{\alpha}}$ is well-
behaved.

We want to point out that if the argument of $G_{j}$ has the opposite sign,
$G_{j}(-i\nu_{n}+\cdots)$ then we have to choose $-f(-z)$ instead of $f(z)$ to
generate an equally simple result with a single non-vanishing term. In other
words:

$\displaystyle T\sum_{i\nu_{n}}\Delta(\mp
i\nu_{n})\prod_{j=1}^{N}G_{j}(-i\nu_{n}+i\omega+\sum_{\alpha}y_{\alpha}a_{\alpha}^{j}+\sum_{\beta}i\nu_{m}^{\beta}a_{\beta}^{j})=$
$\displaystyle\pm\int dyA_{c}(\mp
y)f(-y)\prod_{j=1}^{N}G_{j}(-y+i\omega+\sum_{\alpha}y_{\alpha}a_{\alpha}^{j}+\sum_{\beta}i\nu_{m}^{\beta}a_{\beta}^{j}).$
(56)

Of course Eq. 56 and 53 are equivalent because $i\nu_{n}$ go over both the
positive and negative frequencies, and the integral over $y$ is unrestricted,
allowing a simple substitution, $y\rightarrow-y$. However, notice that the
sign in the fermi function is tied to the sign in the pseudo Green’s function:
$f(\pm y)G_{j}(\pm y+\cdots)$.

Of course, normally, we have the freedom to choose between $f(z)$ and
$-f(-z)$, which have the same residues. We just need to be careful to make the
resulting expression converge for large $z$. However, here we select the sign
in $f(\pm z)$ such that only a single term survives the projection. For the
sake of clarity, lets assume that $G_{j}$ has a positive argument
$G_{j}(i\nu_{n}+\cdots)$, like in Eq. 53. Still, we would select $-f(-z)$
instead of $f(z)$, which has the same residues. In this case, several terms
would survive the projection, needing to be combined together. Of course, they
would give an identical result, but the derivation in the latter case is more
tedious.

In summary, we have just shown that the sign of the fermi function argument,
$f(y)$, in the final results is always equal to the sign in the pseudo-Green’s
function argument, namely, $f(\pm y)G_{j}(\pm y+...)$ in Eq. (53). We now
recognize that all Matsubara sums can be carried out by this simple rule.
Hence, the analytic continuation of a Feyman diagram is really simple: we
replace $\Delta(i\nu)$ with $A_{c}(y)f(\pm y)$, where the sign in $f$ argument
needs to coinside with the sign of $y$ in pseudo Green’s functions.

To make the rules very clear, we provide a concrete form for the Luttinger-
Ward functional corresponding to a generic diagram on the imaginary axis:

$\displaystyle\Phi\propto
T\sum_{i\omega}\left(\prod_{\alpha=1}^{n}T\sum_{i\nu_{\alpha}}\Delta_{\alpha}(a_{\alpha}i\nu_{\alpha})\right)\prod_{j=1}^{2n}G_{j}(i\omega+\sum_{\beta}t_{j,\beta}i\nu_{\beta}).$
(57)

Here $t_{j,\beta}$ is a set of numbers, which can only take values of either
+1,-1 or 0, and $a_{\alpha}$ is either +1 or -1. The self-energy contribution
to pseudo-particle $G_{p}$ is a functional derivative, i.e.,

$\displaystyle\Sigma_{p}(i\omega)=\frac{\delta\Phi}{\delta G_{p}(i\omega)}$
(58)

and becomes

$\displaystyle\Sigma_{p}(i\omega)\propto\left(\prod_{\alpha=1}^{n}T\sum_{i\nu_{\alpha}}\Delta_{\alpha}(a_{\alpha}i\nu_{\alpha})\right)\prod_{j=1;j\neq
p}^{2n}G_{j}(i\omega+\sum_{\beta}t_{j,\beta}^{(p)}i\nu_{\beta}).$ (59)

Here we used a short notation $t_{j,\beta}-t_{p,\beta}\equiv
t_{j,\beta}^{(p)}$. As we will demonstrate below $t_{j,\beta}^{(p)}$ also
takes the values +1,-1, or 0. Next we slightly manipulate the expression, so
that we can use Eqs. 56 or 53 on it:

$\displaystyle\Sigma_{p}(i\omega)\propto\left(\prod_{\alpha=1}^{n}T\sum_{i\nu_{\alpha}}\Delta_{\alpha}(a_{\alpha}i\nu_{\alpha})\right)\prod_{j=1;j\neq
p}^{2n}G_{j}(t^{p}_{j,\alpha}i\nu_{\alpha}+i\omega+\sum_{\beta\neq\alpha}t^{p}_{j,\beta}i\nu_{\beta}).$
(60)

Before we carry out the summations, we note that due to projection, only
diagrams with a single back-bone survive. As a result, for any $j=1...2n$, the
variables $t^{(p)}_{j,\alpha}$ must have the same sign. Specifically, at fixed
$\alpha$ and $p$, we either have $t^{(p)}_{j,\alpha}\geq 0$ or
$t^{(p)}_{j,\alpha}\leq 0$ for any $j$. This restriction arises because, for a
given loop $\alpha$, we can only have either clockwise
($t^{(p)}_{j,\alpha}\geq 0$) or anti-clockwise ($t^{(p)}_{j,\alpha}\leq 0$)
orientation of the loop through pseudo Green’s functions. To make further
derivation easier, we break $t_{j,\alpha}$ into two parts
$t^{(p)}_{j,\alpha}=c^{(p)}_{j,\alpha}b^{(p)}_{\alpha}$, where
$c^{(p)}_{j,\alpha}$ can only take values of 0 or 1, and $b^{(p)}_{\alpha}$
can only be +1 or -1. Therefore, $b^{(p)}_{\alpha}$ carries the sign and
indicates whether the particular $\alpha$ loop is oriented clockwise or
counter-clockwise, while $c^{(p)}_{j,\alpha}=|t^{(p)}_{j,\alpha}|$ is nonzero
only for those pseudo-particles $j$ that are in a particular frequency loop
$\alpha$. We can now apply Eqs. 56 or 53 to obtain

$\displaystyle\Sigma_{p}(\omega+\lambda)\propto\left(\prod_{\alpha=1}^{n}a_{\alpha}b^{(p)}_{\alpha}\right)\left(\prod_{\alpha=1}^{n}\int
dy_{\alpha}A^{c}_{\alpha}(a_{\alpha}y_{\alpha})f(b^{(p)}_{\alpha}y_{\alpha})\right)\prod_{j=1;j\neq
p}^{2n}G_{j}(\omega+\lambda+i\delta+\sum_{\beta}t^{p}_{j,\beta}y_{\beta}).$
(61)

In this equation we also changed $i\omega$ to $\omega+\lambda+i\delta$, which
accounts for analytic continuation and projection. We now notice that appart
from overal sign $\prod_{\alpha=1}^{n}a_{\alpha}b^{(p)}_{\alpha}$, which we
here determine in alternative way through topology and direct matrix element
calculation, the two expressions on real axis (Eq. 61) and imaginary axis
(Eq.59) have nearly identical form. The rule to obtain the real axis
expression from the imaginary axis expression is simple:

  * •

replace $i\nu_{\alpha}$ with real frequency $y_{\alpha}$ and replace the
Matsubara sum over $i\nu_{\alpha}$ with integral over $y_{\alpha}$.

  * •

replace $i\omega$ with $\omega+\lambda+i\delta$

  * •

replace hybridization function $\Delta_{\alpha}$ with its spectral function
$A^{c}_{\alpha}$.

  * •

each hybridization function should be accompanied by the fermi function
$f(b^{(p)}_{\alpha}y_{\alpha})$, in which the sign is determined by the sign
of the arguments in pseudo Green’s functions
$b_{\alpha}=\textrm{sign}(t^{(p)}_{,\alpha})$.

To be concrete, we will next discuss the particular diagram in Fig. 6. When we
calculate the pseudo self-energy $\Sigma_{2}(\omega)$ the arguments in the
pseudo Green’s functions are given in table $t^{(2)}_{j,\alpha}$ in Eq. 105
and are also denoted in the figure. Each row represents an argument for one of
the functions. For example, $G_{0}$ appears as $G_{0}(\omega+y_{1}+y_{2})$
once the Matsubara frequencies $i\nu_{1}$ and $i\nu_{2}$ are replaced by
integration over $y_{1}$ and $y_{2}$. However, when we calculate the self-
energy $\Sigma_{0}(\omega)$, we need to shift the arguments by
$\omega\rightarrow\omega-y_{1}-y_{2}$, which results in arguments displayed in
the table $t^{(0)}_{j,\alpha}$. We notice that they can be obtained by the
above stated formula
$t^{(0)}_{j,\alpha}=t^{(i)}_{j,\alpha}-t^{(i)}_{0,\alpha}$. The rest of the
self-energies require arguments in the tables $t^{(1)}\cdots t^{(5)}$. Now
that we know the argument in all $G_{j}$’s, we can determine the sign of the
fermi functions. We realize that at fixed representation $i$ and frequency
loop $\alpha$ all $t^{(i)}_{j,\alpha}$ have the same sign for any $j$, i.e.,
all pseudo-Green’s functions $G_{j}$ contain the same sign for argument
$y_{\alpha}$. This sign was before denoted by
$b^{(i)}_{\alpha}=\textrm{sign}(t^{(i)}_{,\alpha})$, and is displayed in Eq.
112. These numbers $b^{(i)}_{\alpha}$ give the sign of the arguments in the
fermi functions Eq. 61

$\displaystyle\left(\begin{array}[]{c|rrr}t^{(2)}&y_{1}&y_{2}&y_{3}\\\
\hline\cr G_{0}&1&1&0\\\ G_{1}&1&0&0\\\ G_{2}&0&0&0\\\ G_{3}&0&0&1\\\
G_{4}&0&1&1\\\
G_{5}&0&1&0\end{array}\right)\qquad\left(\begin{array}[]{c|rrr}t^{(0)}&y_{1}&y_{2}&y_{3}\\\
\hline\cr G_{0}&0&0&0\\\ G_{1}&0&-1&0\\\ G_{2}&-1&-1&0\\\ G_{3}&-1&-1&1\\\
G_{4}&-1&0&1\\\
G_{5}&-1&0&0\end{array}\right)\qquad\left(\begin{array}[]{c|rrr}t^{(1)}&y_{1}&y_{2}&y_{3}\\\
\hline\cr G_{0}&0&1&0\\\ G_{1}&0&0&0\\\ G_{2}&-1&0&0\\\ G_{3}&-1&0&1\\\
G_{4}&-1&1&1\\\ G_{5}&-1&1&0\end{array}\right)$ (83)
$\displaystyle\left(\begin{array}[]{c|rrr}t^{(3)}&y_{1}&y_{2}&y_{3}\\\
\hline\cr G_{0}&1&1&-1\\\ G_{1}&1&0&-1\\\ G_{2}&0&0&-1\\\ G_{3}&0&0&0\\\
G_{4}&0&1&0\\\
G_{5}&0&1&-1\end{array}\right)\qquad\left(\begin{array}[]{c|rrr}t^{(4)}&y_{1}&y_{2}&y_{3}\\\
\hline\cr G_{0}&1&0&-1\\\ G_{1}&1&-1&-1\\\ G_{2}&0&-1&-1\\\ G_{3}&0&-1&0\\\
G_{4}&0&0&0\\\
G_{5}&0&0&-1\end{array}\right)\qquad\left(\begin{array}[]{c|rrr}t^{(5)}&y_{1}&y_{2}&y_{3}\\\
\hline\cr G_{0}&1&0&0\\\ G_{1}&1&-1&0\\\ G_{2}&0&-1&0\\\ G_{3}&0&-1&1\\\
G_{4}&0&0&1\\\ G_{5}&0&0&0\end{array}\right)$ (105)

$a=[-1,1,-1]$ $\displaystyle\begin{array}[]{lrrr}b^{(2)}=[&1,&1,&1]\\\
b^{(0)}=[&-1,&-1,&1]\\\ b^{(1)}=[&-1,&1,&1]\\\ b^{(3)}=[&1,&1,&-1]\\\
b^{(4)}=[&1,&-1,&-1]\\\ b^{(5)}=[&1,&-1,&1]\end{array}$ (112)

For completness we also show the form of the pseudo self-energies for the
example in Fig. 6. We first construct the modified generating functional
$\overline{\Phi}$, which contains only hybridizations functions $A^{c}(y)$ and
pseudo-functions $G_{j}$. The functional form of $\overline{\Phi}$ is
identical to its Matsubara imaginary axis form, except that $i\nu_{l}$ are
replaced by real variables $y_{l}$ and $i\omega$ by $\omega$. For the example
in Fig. 6 the expression is

$\displaystyle\overline{\Phi}=\int
A^{c}_{1}(-y_{1})A^{c}_{2}(y_{2})A^{c}_{3}(-y_{3})G_{0}(\omega+y_{1}+y_{2})G_{1}(\omega+y_{1})G_{2}(\omega)G_{3}(\omega+y_{3})G_{4}(\omega+y_{2}+y_{3})G_{5}(\omega+y_{2})$
(113)

Next we add $n$ fermi functions, one for each hybridization frequency $y_{i}$,
using the above derived rule. The sign of the argument is determined from the
tables $b^{(i)}$ and is identical to the sign of that frequency in pseudo
Green’s function $f(\pm y_{i})G_{p}(\pm y_{i}+\cdots)$. We note that the
tables $t^{(i)}_{j,\alpha}$ are easily generated
($t^{(i)}_{j,\alpha}=t_{j,\alpha}-t_{i,\alpha}$) for each Feynman diagram once
topology is known and frequency loops have been chosen. For the above example
in Fig. 6, we have

$\displaystyle\Sigma_{2}(\omega)=\frac{\delta\overline{\Phi}}{\delta
G_{2}(\omega)}f(y_{1})f(y_{2})f(y_{3})$ (114)
$\displaystyle\Sigma_{0}(\omega+y_{1}+y_{2})=\frac{\delta\overline{\Phi}}{\delta
G_{0}(\omega+y_{1}+y_{2})}f(-y_{1})f(-y_{2})f(y_{3})$ (115)
$\displaystyle\Sigma_{1}(\omega+y_{1})=\frac{\delta\overline{\Phi}}{\delta
G_{1}(\omega+y_{1})}f(-y_{1})f(y_{2})f(y_{3})$ (116)
$\displaystyle\Sigma_{3}(\omega+y_{3})=\frac{\delta\overline{\Phi}}{\delta
G_{3}(\omega+y_{3})}f(y_{1})f(y_{2})f(-y_{3})$ (117)
$\displaystyle\Sigma_{4}(\omega+y_{2}+y_{3})=\frac{\delta\overline{\Phi}}{\delta
G_{4}(\omega+y_{2}+y_{3})}f(y_{1})f(-y_{2})f(-y_{3})$ (118)
$\displaystyle\Sigma_{5}(\omega+y_{2})=\frac{\delta\overline{\Phi}}{\delta
G_{5}(\omega+y_{2})}f(y_{1})f(-y_{2})f(y_{3})$ (119)

Here we emphasized that contribution to all pseudo self-energies can be
computed from the same form of the functional $\overline{\Phi}$ and we do not
actually need to shift variables in $\overline{\Phi}$. This is numerically
much more efficient than shifting variables. The shifting of variables was
used above only to determine the sign of the fermi function variables $f(\pm
y_{\alpha})$.

Finally, we are going to derive the expression for the electron Green’s
function ${\cal G}(\omega)$. On imaginary axis this is obtained by cutting one
of the hybridization propagators, i.e.,

$\displaystyle{\cal
G}(i\nu_{\alpha})=\frac{\partial\Phi_{\Delta}[\\{G_{p}\\}]}{\partial\Delta_{\alpha}(i\nu_{\alpha})}$
(120)

As the external variable is now one of the hybridization frequencies, the
summation over that particular $i\nu_{\alpha}$ is omitted, while the summation
over $i\omega$ must be carried out instead.

The end result on the real axis has quite simple structure, but the derivation
is somewhat lengthly. Above we defined $2n$ matrices $t^{(i)}$, which all
represent equivalent ways of choosing frequency arguments in pseudo Green’s
functions $G_{p}$, and $t^{(i)}$ stands for representation in which $G_{i}$
contains only $i\omega$, i.e., $G_{i}(i\omega)$. Sometimes we omit this
superscript $(i)$, when any representation could be chosen, and we did not yet
determine which one is being chosen.

Next we want to derive some simple relations between above introduced
quantities. We already stated above how to obtain representation $t^{(p)}$
from any other representation (say $t^{(l)}$), namely
$t^{(p)}_{j,\alpha}=t^{(l)}_{j,\alpha}-t^{(l)}_{p,\alpha}$. Consequently
$t^{(p)}_{p,\alpha}=0$ for any $\alpha$, as required in representation $p$
that $G_{p}(i\omega)$ has no hybridization frequency $y_{\alpha}$ in the
argument. Second we want to derive $b^{(p)}$ in representation $p$ from
knowing $b^{(l)}$ and $c^{(l)}$ in $l$ representation. The connection is
simply $b^{(p)}_{\alpha}=b^{(l)}_{\alpha}(-1)^{c^{(l)}_{p,\alpha}}$. Note that
this comes simply from the algebraic solution of eliminating hybridization
frequencies $y_{\alpha}$ in a chosen pseudo-particle. Also note that because
$t^{(p)}_{p,\alpha}={c^{(p)}_{p,\alpha}}=0$ for any $\alpha$ and $p$, the
above equation just leads to identity for the case $l=p$.

In Eq. 57 we stated the Luttinger functional expression on the imaginary axis
for a generic Feynman diagram. Using Eq. 120, we obtain

$\displaystyle{\cal G}_{\alpha}(a_{\alpha}i\nu_{\alpha})\propto
T\sum_{i\omega}\prod_{\beta=1,\beta\neq\alpha}^{n}T\sum_{i\nu_{\beta}}\Delta_{\beta}(a_{\beta}i\nu_{\beta})\prod_{j=1}^{2n}G_{j}(i\omega+\sum_{\beta}t_{j,\beta}i\nu_{\beta})$
(121)

Here $a_{\alpha}=\pm 1$ is the relative orientation of the loop and the
hybridization propagator.

Now we start to carry out Matsubara summations. First we perform $n-1$
summations over the hybridization frequencies $i\nu_{\beta}$, which were
explained above, and give a single term on the real axis, which takes the form

$\displaystyle{\cal
G}_{\alpha}(a_{\alpha}i\nu_{\alpha})\propto\prod_{\beta=1,\beta\neq\alpha}^{n}a_{\beta}b_{\beta}\int
dy_{\beta}A^{c}_{\beta}(a_{\beta}y_{\beta})f(b_{\beta}y_{\beta})\;T\sum_{i\omega}\prod_{j=1}^{2n}G_{j}(i\omega+\sum_{\beta\neq\alpha}t_{j,\beta}y_{\beta}+t_{j,\alpha}i\nu_{\alpha}).$
(122)

Next we sum over $i\omega$, which can be either fermionic or bosonic. We first
rewrite all pseudo Green’s functions in their spectral representation
$G_{j}(z)=\int dx_{j}A(x_{j})/(z-x_{j})$, which gives

$\displaystyle{\cal
G}_{\alpha}(a_{\alpha}i\nu_{\alpha})\propto\prod_{\beta=1,\beta\neq\alpha}^{n}a_{\beta}b_{\beta}\int
dy_{\beta}A^{c}_{\beta}(a_{\beta}y_{\beta})f(b_{\beta}y_{\beta})\;T\sum_{i\omega}\prod_{j=1}^{2n}\int
dx_{j}A_{j}(x_{j})\frac{1}{i\omega+\sum_{\beta\neq\alpha}t_{j,\beta}y_{\beta}+t_{j,\alpha}i\nu_{\alpha}-x_{j}}$
(123)

Now we recognize that for this summation we can use the generalized residue
formula, which can for example be found in Ref. [34],

$\displaystyle
T\sum_{i\omega}\prod_{j=1}^{N}\frac{1}{i\omega+z_{j}}=\sum_{l=1}^{N}f(-z_{l})\prod_{j=1,j\neq
l}^{N}\frac{1}{z_{j}-z_{l}}$ (124)

Here $f(z)$ is the fermi function $f(z)=1/(\exp(\beta z)+1)$ or minus bose
function $f(z)=-1/(\exp(\beta z)-1)$ if $i\omega$ is fermionic or bosonic,
respectively. In our case
$z_{j}=\sum_{\beta\neq\alpha}t_{j,\beta}y_{\beta}+t_{j,\alpha}i\nu_{\alpha}-x_{j}$,
hence the result is

$\displaystyle{\cal
G}_{\alpha}(a_{\alpha}i\nu_{\alpha})\propto\prod_{\beta=1,\beta\neq\alpha}^{n}a_{\beta}b_{\beta}\int
dy_{\beta}A^{c}_{\beta}(a_{\beta}y_{\beta})f(b_{\beta}y_{\beta})$
$\displaystyle\sum_{l=1}^{2n}\int
dx_{l}A_{l}(x_{l})f(x_{l}-\sum_{\beta\neq\alpha}t_{l,\beta}y_{\beta}-t_{l,\alpha}i\nu_{\alpha})\prod_{j=1,j\neq
l}^{2n}\int
dx_{j}\frac{A_{j}(x_{j})}{\sum_{\beta\neq\alpha}(t_{j,\beta}-t_{l,\beta})y_{\beta}+(t_{j,\alpha}-t_{l,\alpha})i\nu_{\alpha}+x_{l}-x_{j}}.$
(125)

Now we rewrite the spectral representation of $G_{j}$ back into their closed
form for $j\neq l$ terms, to obtain

$\displaystyle{\cal
G}_{\alpha}(a_{\alpha}i\nu_{\alpha})\propto\prod_{\beta=1,\beta\neq\alpha}^{n}a_{\beta}b_{\beta}\int
dy_{\beta}A^{c}_{\beta}(a_{\beta}y_{\beta})f(b_{\beta}y_{\beta})$
$\displaystyle\sum_{l=1}^{2n}\int
dx_{l}A_{l}(x_{l})f(x_{l}-\sum_{\beta\neq\alpha}t_{l,\beta}y_{\beta}-t_{l,\alpha}i\nu_{\alpha})\prod_{j=1,j\neq
l}^{2n}G_{j}(x_{l}+t^{(l)}_{j,\alpha}i\nu_{\alpha}+\sum_{\beta\neq\alpha}t^{(l)}_{j,\beta}y_{\beta}),$
(126)

where we used $t^{(l)}_{j,\alpha}=t_{j,\alpha}-t_{l,\alpha}$. We introduced
before the split of $t_{l,\alpha}$ into $b$ and $c$, i.e.,
$t_{l,\alpha}=b_{\alpha}c_{l,\alpha}$, where $c_{l,\alpha}=|t_{l,\alpha}|$,
hence $e^{-\beta i\nu_{\alpha}t_{l,\alpha}}=(-1)^{c_{l,\alpha}}$, because
external Matsubara frequency $i\nu_{\alpha}$ is fermionic. Because of
projection, we also need to shift $x_{l}$ frequency to $x_{l}+\lambda$, hence
the $f$ function becomes

$f(x_{l}-\sum_{\beta\neq\alpha}t_{l,\beta}y_{\beta}-t_{l,\alpha}i\nu_{\alpha})\rightarrow\pm(-1)^{c_{l,\alpha}}\exp\left(-\beta(x_{l}+\lambda-\sum_{\beta\neq\alpha}t_{l,\beta}y_{\beta})\right).$

Here $+1$ and $-1$ corresponds to $i\omega$ being fermionic or bosonic. We
will therefore introduce notation $(-1)^{\omega_{boson}}$. After the frequency
$x_{l}$ is shifted by $\lambda$ to account for projection, we obtain

$\displaystyle{\cal
G}_{\alpha}(a_{\alpha}i\nu_{\alpha})\propto(-1)^{\omega_{boson}}\prod_{\beta=1,\beta\neq\alpha}^{n}a_{\beta}b_{\beta}\int
dy_{\beta}A^{c}_{\beta}(a_{\beta}y_{\beta})f(b_{\beta}y_{\beta})$
$\displaystyle\sum_{l=1}^{2n}e^{\beta\sum_{\beta\neq\alpha}t_{l,\beta}y_{\beta}}(-1)^{c_{l,\alpha}}\int
dx_{l}A_{l}(x_{l}+\lambda)e^{-\beta(x_{l}+\lambda)}\prod_{j=1,j\neq
l}^{2n}G_{j}(x_{l}+\lambda+t^{(l)}_{j,\alpha}i\nu_{\alpha}+\sum_{\beta\neq\alpha}t^{(l)}_{j,\beta}y_{\beta})$
(127)

We next combine the term $e^{\beta\sum_{\beta\neq\alpha}t_{l,\beta}y_{\beta}}$
with the fermi functions under the product
$\prod_{\beta\neq\alpha}f(b_{\beta}y_{\beta})$. We defined above that
$t_{l,\beta}=c_{l,\beta}b_{\beta}$ and $c_{l,\beta}$ is either 0 or 1. We use
the fact that $f(x)e^{\beta x}=f(-x)$ to write

$f(b_{\beta}y_{\beta})e^{\beta
y_{\beta}t_{l,\beta}}=f(b_{\beta}y_{\beta}(-1)^{c_{l,\beta}}).$

This now gives

$\displaystyle{\cal
G}_{\alpha}(a_{\alpha}i\nu_{\alpha})\propto\sum_{l=1}^{2n}(-1)^{c_{l,\alpha}+\omega_{boson}}\prod_{\beta=1,\beta\neq\alpha}^{n}a_{\beta}b_{\beta}\int
dy_{\beta}A^{c}_{\beta}(a_{\beta}y_{\beta})f(b_{\beta}y_{\beta}(-1)^{c_{l,\beta}})$
$\displaystyle\int dxA_{l}(x+\lambda)e^{-\beta(x+\lambda)}\prod_{j=1,j\neq
l}^{2n}G_{j}(x+\lambda+t^{(l)}_{j,\alpha}i\nu_{\alpha}+\sum_{\beta\neq\alpha}t^{(l)}_{j,\beta}y_{\beta})$
(128)

Next we use the fact $b_{\alpha}(-1)^{c_{l,\alpha}}=b^{(l)}_{\alpha}$, and
$b_{\alpha}^{2}=1$ and $a_{\alpha}^{2}=1$, therefore the sum over $l$ in the
above equation is the sum over $2n$ representations $t^{(l)}$. We have

$\displaystyle{\cal
G}_{\alpha}(a_{\alpha}i\nu_{\alpha})\propto\left((-1)^{\omega_{boson}}\prod_{\beta=1}^{n}a_{\beta}b_{\beta}\right)\sum_{l=1}^{2n}a_{\alpha}b^{(l)}_{\alpha}\prod_{\beta=1,\beta\neq\alpha}^{n}\int
dy_{\beta}A^{c}_{\beta}(a_{\beta}y_{\beta})f(b^{(l)}_{\beta}y_{\beta})$
$\displaystyle\int dxA_{l}(x+\lambda)e^{-\beta(x+\lambda)}\prod_{j=1,j\neq
l}^{2n}G_{j}(x+\lambda+t^{(l)}_{j,\alpha}i\nu_{\alpha}+\sum_{\beta\neq\alpha}t^{(l)}_{j,\beta}y_{\beta})$
(129)

Next we perform analytic continuation and let $i\nu_{\alpha}\rightarrow
y_{\alpha}+i\delta$. We will generate $G_{j}(\cdots)$,
$\textrm{Re}(G_{j}(\cdots))$, or $G^{*}_{j}(\cdots)$ if the coefficient in
front of $i\nu$ is 1, 0, or -1. We will therefore define

$\displaystyle G_{j}^{(0)}(x)\equiv\textrm{Re}(G_{j}(x))$ (130) $\displaystyle
G_{j}^{(1)}(x)\equiv G_{j}(x)$ (131) $\displaystyle G_{j}^{(-1)}(x)\equiv
G_{j}^{*}(x)$ (132) $\displaystyle\widetilde{A}_{j}(x)=A_{j}(x)e^{-\beta x}$
(133)

and in terms of these pseudo Green’s function, we have

$\displaystyle{\cal
G}_{\alpha}(a_{\alpha}y_{\alpha})\propto\left((-1)^{\omega_{boson}}\prod_{\beta=1}^{n}a_{\beta}b_{\beta}\right)\sum_{l=1}^{2n}a_{\alpha}b^{(l)}_{\alpha}\prod_{\beta=1,\beta\neq\alpha}^{n}\int
dy_{\beta}A^{c}_{\beta}(a_{\beta}y_{\beta})f(b^{(l)}_{\beta}y_{\beta})$
$\displaystyle\int dx\widetilde{A}_{l}(x+\lambda)\prod_{j=1,j\neq
l}^{2n}G^{(a_{\alpha}t^{(l)}_{j,\alpha})}_{j}(x+\lambda+\sum_{\beta}t^{(l)}_{j,\beta}y_{\beta})$
(134)

The first bracket gives just an overal sign, which needs to be combined with
other signs (like number of fermionic loops, and perturbation order) that we
ignored at the beginning, because we determine them in alternative way, i.e.,
in the same way as in conventional CTQMC. We will use notation

$\displaystyle(-1)^{sg}\equiv(-1)^{\omega_{boson}}\prod_{\beta=1}^{n}a_{\beta}b_{\beta}.$
(135)

It can be shown that this sign is independent of the choice of the loops (the
choice of representation) and is a property of the topology of diagram.

Finally, we got the following expression for the electron Green’s function

$\displaystyle{\cal
G}_{\alpha}(a_{\alpha}y_{\alpha})\propto(-1)^{sg}\sum_{l=1}^{2n}a_{\alpha}b^{(l)}_{\alpha}\prod_{\beta=1,\beta\neq\alpha}^{n}\int
dy_{\beta}A^{c}_{\beta}(a_{\beta}y_{\beta})f(b^{(l)}_{\beta}y_{\beta})$
$\displaystyle\int dx\widetilde{A}_{l}(x+\lambda)\prod_{j=1,j\neq
l}^{2n}G^{(a_{\alpha}t^{(l)}_{j,\alpha})}_{j}(x+\lambda+\sum_{\beta}t^{(l)}_{j,\beta}y_{\beta})$
(136)

While this is a closed expression that could be coded, it is numerically more
convenient to compute all Green’s function in the $l$ sum using the same set
of arguments. To get that, we need to transform from $(l)$ representation back
to the original representation in which the diagram was initially constructed.
We will therefore shift the frequency $x$ to
$x=\omega+\sum_{\beta}t_{l,\beta}$ to obtain

$\displaystyle{\cal G}_{\alpha}(a_{\alpha}y_{\alpha})\propto$
$\displaystyle(-1)^{sg}\sum_{l=1}^{2n}a_{\alpha}b^{(l)}_{\alpha}\int
d\omega\prod_{\beta=1,\beta\neq\alpha}^{n}\int
dy_{\beta}A^{c}_{\beta}(a_{\beta}y_{\beta})f(b^{(l)}_{\beta}y_{\beta})$ (137)
$\displaystyle\widetilde{A}_{l}(\omega+\lambda+\sum_{\beta}t_{l,\beta})\prod_{j=1,j\neq
l}^{2n}G^{(a_{\alpha}t^{(l)}_{j,\alpha})}_{j}(\omega+\lambda+\sum_{\beta}t_{j,\beta}y_{\beta})$

so that all pseudo Green’s functions can be evaluated in any representation in
which the diagram is initially expressed. This is the final expression for the
real axis contribution to the Green’s function and is the central result of
this paper. We notice that all the arguments of the propagators are now
written in an arbitrary representation. We still have the sum over $l$, which
runs over $2n$ representations, but we only use representation $t^{(l)}$ to
determine the sign of the fermi function $f(b^{(l)}_{\beta}y_{\beta})$ and to
figure out whether pseudo Green’s function needs to be conjugated or we need
to take the real part $G^{(a_{\alpha}t^{(l)}_{j,\alpha})}_{j}$. Furthermore,
$a_{\alpha}t^{(l)}_{j,\alpha}$ when not zero, has the same sign for all $j$,
and its value is $a_{\alpha}b_{\alpha}^{(l)}$, because
$t^{(l)}_{j,\alpha}=b_{\alpha}^{(l)}c^{(l)}_{j,\alpha}$ where
$c^{(l)}_{j,\alpha}$ can only be 0 or 1. As a consequence, the sign of each
term in the sum ($a_{\alpha}b_{\alpha}^{(l)}$) is easy to determine: If the
$l$ term has any pseudo-Green’s function propagator conjugated, then the sign
is negative as all pseudo-Green’s functions in this term are either conjugated
or we need to take their real part. On the other hand, if any pseudo-Green’s
function in the $l$ term is non conjugated, then none is conjugated, and the
sign is positive.

Finally we state the Feynman rules for the real axis electron Green’s function
calculation If we compare this real axis expression Eq. 137 with the imaginary
axis expression Eq. 121 we notice that the following substitutions need to be
made:

  * •

hybridization function has to be evaluated on the real axis substituting
$\Delta_{\beta}(a_{\beta}i\nu_{\beta})\rightarrow
A^{c}_{\beta}(a_{\beta}y_{\beta})=-\frac{1}{\pi}\textrm{Im}(\Delta_{\beta}(a_{\beta}y_{\beta}))$

  * •

the sums over hybridization function Matsubara sums $i\nu_{\beta}$ are
replaced by integrals over real variables $\int dy_{\beta}$, and the sum over
$i\omega$ is replaced by the integral over $\omega$.

  * •

We need to sum over $2n$ terms in which one of the pseudo-propagators
$G_{l}(i\omega+\cdots)$ is replaced by its spectral function
$\widetilde{A}_{l}(\omega+\lambda+\cdots)$ while the rest $2n-1$ pseudo
Green’s functions are substituted with the real part, the retarder, or the
advanced Green’s function on the real axis. This is determined by the argument
$a_{\alpha}t^{(l)}_{j,\alpha}$, which can be 0, 1, or -1, i.e.,
$G^{(a_{\alpha}t^{(l)}_{j,\alpha})}_{j}(\omega+\lambda+\cdots)$. The overal
sign $a_{\alpha}b^{(l)}_{\alpha}$ of the term is positive when retarder
$G_{j}$ appear, and negative when advanced $G_{j}$ appear. Note that mixed
term with both advanced and retarded $G_{j}$ is not possible.

  * •

Finally, each term requires different set of $n-1$ fermi functions, which are
given by $b^{(l)}_{\beta}$ in representation $l$, i.e.,
$\prod_{\beta\neq\alpha}f(b^{(l)}_{\beta}y_{\beta})$. The same set of fermi
functions are required to compute pseudo self-energy $\Sigma_{l}$ on the real
axis.

We next apply these rules to our example in Fig. 6. We start with the modified
generating functional $\overline{\Phi}$, defined in Eq. 113. We than construct
2n=6 terms in which exactly one $G$ is replaced by its spectral function
$\widetilde{A}$ and some are replaced by their real part. For example, when we
compute ${\cal G}_{1}(-y_{1})$, and we concentrate on $(l)=(0)$, we have
$a_{1}=-1$ and $t^{(0)}_{:,1}=(0,0,-1,-1,-1,-1)$, hence
$a_{1}t^{(0)}_{:,1}=(0,0,1,1,1,1)$. The zeroth-component $G_{0}$ is replace by
the spectral function $\widetilde{A}_{0}$, while the component 1 vanishes, and
hence requires $G_{1}$ to be replaced by the real part $\textrm{Re}G_{1}$. The
rest of the Green’s functions should be retarded. The sign of the term is
positive. In the second term $(l)=1$ we have the same
$a_{1}t^{(1)}_{:,1}=(0,0,1,1,1,1)$, hence only for $G_{0}$ we take the real
part, while $G_{1}$ is replaced by $\widetilde{A}_{1}$. The next four terms
have the same $a_{1}t^{(2)}_{:,1}=(-1,-1,0,0,0,0)$, which requires complex
conjugated $G_{0}$ and $G_{1}$, while the rest of $G$’s are replaced by their
real part, except for one, which requires the spectral function. The explicit
form is:

$\displaystyle{\cal G}_{1}(-y_{1})\braket{Q}=$
$\displaystyle\left(\frac{\delta\overline{\Phi}}{\delta
A^{c}_{1}(-y_{1})}\right)\left(f(-y_{2})f(y_{3})\frac{\widetilde{A}_{0}(\omega+y_{1}+y_{2})G^{\prime}_{1}(\omega+y_{1})}{G_{0}(\omega+y_{1}+y_{2})G_{1}(\omega+y_{1})}+f(y_{2})f(y_{3})\frac{G^{\prime}_{0}(\omega+y_{1}+y_{2})\widetilde{A}_{1}(\omega+y_{1})}{G_{0}(\omega+y_{1}+y_{2})G_{1}(\omega+y_{1})}\right)$
(138) $\displaystyle-$ $\displaystyle\left(\frac{\delta\overline{\Phi}}{\delta
A^{c}_{1}(-y_{1})}\right)^{*}\left(f(y_{2})f(y_{3})\frac{\widetilde{A}_{2}(\omega)G^{\prime}_{3}(\omega+y_{3})G^{\prime}_{4}(\omega+y_{2}+y_{3})G^{\prime}_{5}(\omega+y_{2})}{G^{*}_{2}(\omega)G^{*}_{3}(\omega+y_{3})G^{*}_{4}(\omega+y_{2}+y_{3})G^{*}_{5}(\omega+y_{2})}\right)$
$\displaystyle-$ $\displaystyle\left(\frac{\delta\overline{\Phi}}{\delta
A^{c}_{1}(-y_{1})}\right)^{*}\left(f(y_{2})f(-y_{3})\frac{G^{\prime}_{2}(\omega)\widetilde{A}_{3}(\omega+y_{3})G^{\prime}_{4}(\omega+y_{2}+y_{3})G^{\prime}_{5}(\omega+y_{2})}{G^{*}_{2}(\omega)G^{*}_{3}(\omega+y_{3})G^{*}_{4}(\omega+y_{2}+y_{3})G^{*}_{5}(\omega+y_{2})}\right)$
$\displaystyle-$ $\displaystyle\left(\frac{\delta\overline{\Phi}}{\delta
A^{c}_{1}(-y_{1})}\right)^{*}\left(f(-y_{2})f(-y_{3})\frac{G^{\prime}_{2}(\omega)G^{\prime}_{3}(\omega+y_{3})\widetilde{A}_{4}(\omega+y_{2}+y_{3})G^{\prime}_{5}(\omega+y_{2})}{G^{*}_{2}(\omega)G^{*}_{3}(\omega+y_{3})G^{*}_{4}(\omega+y_{2}+y_{3})G^{*}_{5}(\omega+y_{2})}\right)$
$\displaystyle-$ $\displaystyle\left(\frac{\delta\overline{\Phi}}{\delta
A^{c}_{1}(-y_{1})}\right)^{*}\left(f(-y_{2})f(y_{3})\frac{G^{\prime}_{2}(\omega)G^{\prime}_{3}(\omega+y_{3})G^{\prime}_{4}(\omega+y_{2}+y_{3})\widetilde{A}_{5}(\omega+y_{2})}{G^{*}_{2}(\omega)G^{*}_{3}(\omega+y_{3})G^{*}_{4}(\omega+y_{2}+y_{3})G^{*}_{5}(\omega+y_{2})}\right)$

The fermi functions are constructed from quantities $b^{(l)}_{\alpha}$ listed
above. For the second contribution to the electron Green’s function
$\alpha=2$, we again need to generate 6 terms, which are:

$\displaystyle{\cal G}_{2}(y_{2})\braket{Q}=$
$\displaystyle\left(\frac{\delta\overline{\Phi}}{\delta
A^{c}_{2}(y_{2})}\right)\left(f(-y_{1})f(y_{3})\frac{\widetilde{A}_{1}(\omega+y_{1})G^{\prime}_{2}(\omega)G^{\prime}_{3}(\omega+y_{2})}{G_{1}(\omega+y_{1})G_{2}(\omega)G_{3}(\omega+y_{2})}+f(y_{1})f(y_{3})\frac{G^{\prime}_{1}(\omega+y_{1})\widetilde{A}_{2}(\omega)G^{\prime}_{3}(\omega+y_{2})}{G_{1}(\omega+y_{1})G_{2}(\omega)G_{3}(\omega+y_{2})}\right)$
(139) $\displaystyle+$ $\displaystyle\left(\frac{\delta\overline{\Phi}}{\delta
A^{c}_{2}(y_{2})}\right)\left(f(y_{1})f(-y_{3})\frac{G^{\prime}_{1}(\omega+y_{1})G^{\prime}_{2}(\omega)\widetilde{A}_{3}(\omega+y_{2})}{G_{1}(\omega+y_{1})G_{2}(\omega)G_{3}(\omega+y_{2})}\right)$
$\displaystyle-$ $\displaystyle\left(\frac{\delta\overline{\Phi}}{\delta
A^{c}_{2}(y_{2})}\right)^{*}\left(f(-y_{1})f(y_{3})\frac{\widetilde{A}_{0}(\omega+y_{1}+y_{2})G^{\prime}_{4}(\omega+y_{2}+y_{3})G^{\prime}_{5}(\omega+y_{2})}{G^{*}_{0}(\omega+y_{1}+y_{2})G^{*}_{4}(\omega+y_{2}+y_{3})G^{*}_{5}(\omega+y_{2})}\right)$
$\displaystyle-$ $\displaystyle\left(\frac{\delta\overline{\Phi}}{\delta
A^{c}_{2}(y_{2})}\right)^{*}\left(f(y_{1})f(-y_{3})\frac{G^{\prime}_{0}(\omega+y_{1}+y_{2})\widetilde{A}_{4}(\omega+y_{2}+y_{3})G^{\prime}_{5}(\omega+y_{2})}{G^{*}_{0}(\omega+y_{1}+y_{2})G^{*}_{4}(\omega+y_{2}+y_{3})G^{*}_{5}(\omega+y_{2})}\right)$
$\displaystyle-$ $\displaystyle\left(\frac{\delta\overline{\Phi}}{\delta
A^{c}_{2}(y_{2})}\right)^{*}\left(f(y_{1})f(y_{3})\frac{G^{\prime}_{0}(\omega+y_{1}+y_{2})G^{\prime}_{4}(\omega+y_{2}+y_{3})\widetilde{A}_{5}(\omega+y_{2})}{G^{*}_{0}(\omega+y_{1}+y_{2})G^{*}_{4}(\omega+y_{2}+y_{3})G^{*}_{5}(\omega+y_{2})}\right)$