Cancellation of vacuum diagrams and long-time limit in out-of-equilibrium diagrammatic Quantum Monte Carlo

TL;DR

This paper reformulates a real-time diagrammatic Quantum Monte Carlo method in the Larkin-Ovchinnikov basis, enabling efficient long-time limit calculations for out-of-equilibrium steady states, while analyzing the impact on the sign problem.

Contribution

It introduces a new formulation of the Quantum Monte Carlo algorithm that simplifies vacuum diagram cancellation and addresses long-time limits in out-of-equilibrium systems.

Findings

Vacuum diagrams vanish directly in the new formulation.

The algorithm efficiently targets interaction times near measurement.

The sign problem persists despite vacuum diagram cancellation.

Abstract

We express the recently introduced real-time diagrammatic Quantum Monte Carlo, Phys. Rev. B 91, 245154 (2015), in the Larkin-Ovchinnikov basis in Keldysh space. Based on a perturbation expansion in the local interaction $U$, the special form of the interaction vertex allows to write diagrammatic rules in which vacuum Feynman diagrams directly vanish. This reproduces the main property of the previous algorithm, without the cost of the exponential sum over Keldysh indices. In an importance sampling procedure, this implies that only interaction times in the vicinity of the measurement time contribute. Such an algorithm can then directly address the long-time limit needed in the study of steady states in out-of-equilibrium systems. We then implement and discuss different variants of Monte Carlo algorithms in the Larkin-Ovchinnikov basis. A sign problem reappears, showing that the cancellation of vacuum diagrams has no direct impact on it.

Equations

$${\hat{G}}_{\sigma }(t,t^{\prime })\equiv -i⟨T_{\mathcal{C}}{c}_{\sigma }(t)c_{\sigma }^{†}(t^{\prime })⟩,$$

$$G_{\sigma }^{<}(t,t^{\prime })$$

$$G_{\sigma }^{>}(t,t^{\prime })$$

$$G_{\sigma }^{++}(t,t^{\prime })$$

$$G_{\sigma }^{--}(t,t^{\prime })$$

$$\begin{gathered}\fmfgraph*(40,20)\fmfset{arrow_{l}en}{3mm}\fmfleft{v1}\fmfright{v2}\fmf{plain_{a}rrow,label.dist=0,label=\vspace{.5cm}$\sigma$}{v1,v2}\fmflabel{$t', \beta$}{v1}\fmflabel{$t, \alpha$}{v2}\end{gathered}\hskip 25.6073pt=i\left(\hat{g}_{\sigma}\right)_{\alpha\beta}(t,t^{\prime}).$$

$$\begin{gathered}\fmfgraph*(40,30)\fmfset{arrow_{l}en}{3mm}\fmfleft{i1,i2}\fmfright{o1,o2}\fmf{plain_{a}rrow}{i1,v}\fmf{plain_{a}rrow}{v,o2}\fmf{plain_{a}rrow}{o1,v}\fmf{plain_{a}rrow}{v,i2}\fmfdot{v}\fmflabel{$\alpha$,$\uparrow$}{i1}\fmflabel{$\alpha$,$\downarrow$}{o1}\fmflabel{$\alpha$,$\uparrow$}{i2}\fmflabel{$\alpha$,$\downarrow$}{o2}\fmflabel{$\{t, \alpha\}$}{v}\end{gathered}\hskip 14.22636pt=-i\alpha U.$$

$$-iU\left(m_{+}\otimes m_{+}-m_{-}\otimes m_{-}\right),$$

$$\begin{gathered}\fmfgraph*(40,20)\fmfset{arrow_{l}en}{3mm}\fmfleft{i1}\fmfright{o1}\fmf{plain_{a}rrow,tension=2,label.dist=0,label=\vspace{.3cm}$\sigma$}{i1,v1}\fmf{plain_{a}rrow,tension=2,label.dist=0,label=\vspace{-.6cm}$\{t,, \alpha\}$}{v1,o1}\fmffreeze\fmf{phantom,label.dist=0,label=\vspace{.6cm}$h$}{v1,o1}\fmflabel{$\alpha$}{i1}\fmflabel{$\alpha$}{o1}\fmfv{decor.shape=pentagram,decor.filled=full,decor.size=2thick}{v1}\end{gathered}\hskip 14.22636pt=-i\alpha h.$$

$$\begin{gathered}\fmfgraph*(50,20)\fmfset{arrow_{l}en}{3mm}\fmfleft{i1,i2}\fmfright{o1,o2}\fmf{plain_{a}rrow,tension=2,label.dist=0,label=\vspace{.3cm}$\sigma$}{i1,v1}\fmf{plain_{a}rrow,tension=2}{v1,o1}\fmf{phantom}{i2,v2,o2}\fmflabel{$\{t, \alpha\}$}{v1}\fmflabel{$\alpha$}{i1}\fmflabel{$\alpha$}{o1}\fmfdot{v1}\fmffreeze\fmf{plain_{a}rrow,right,tension=2}{v1,v2,v1}\fmf{phantom,label.dist=0,label=$\bar{\sigma}$}{v1,v2}\end{gathered}\hskip 14.22636pt=\alpha Ug^{\alpha\alpha}_{\bar{\sigma}}(t,t).$$

$$g_{\sigma }^{\alpha \alpha }(t,t)=g_{\sigma }^{<}(t,t).$$

$$\begin{gathered}\fmfgraph*(50,20)\fmfset{arrow_{l}en}{3mm}\fmfleft{i1,i2}\fmfright{o1,o2}\fmf{plain_{a}rrow,tension=2,label.dist=0,label=\vspace{.3cm}$\sigma$}{i1,v1}\fmf{plain_{a}rrow,tension=2}{v1,o1}\fmf{phantom}{i2,v2,o2}\fmflabel{$\{t, \alpha\}$}{v1}\fmflabel{$\alpha$}{i1}\fmflabel{$\alpha$}{o1}\fmfdot{v1}\fmffreeze\fmf{plain_{a}rrow,right,tension=2}{v1,v2,v1}\fmf{phantom,label.dist=0,label=$\bar{\sigma}$}{v1,v2}\end{gathered}\hskip 14.22636pt=\hskip 19.91684pt\begin{gathered}\fmfgraph*(40,20)\fmfset{arrow_{l}en}{3mm}\fmfleft{i1}\fmfright{o1}\fmf{plain_{a}rrow,tension=2,label.dist=0,label=\vspace{-.4cm}$\sigma$}{i1,v1}\fmf{plain_{a}rrow,tension=2,label.dist=0,label=\vspace{-.6cm}$\{t,,\alpha\}$}{v1,o1}\fmffreeze\fmf{phantom,label.dist=0,label=\vspace{.6cm}$iUg_{\bar{\sigma}}^<(t,,t)$}{v1,o1}\fmflabel{$\alpha$}{i1}\fmflabel{$\alpha$}{o1}\fmfv{decor.shape=pentagram,decor.filled=full,decor.size=2thick}{v1}\end{gathered}$$

$${\mathcal{M}}_n$$

$$={\int }_0^{t_{\text{max}}}\mathrm{d}{t}_1\dots \mathrm{d}{t}_n\sum _{{\alpha }_1\dots {\alpha }_n}{\mathcal{M}}_n^{\pm }(\{t_i,{\alpha }_i{\}}_{1\le i\le n}),$$

$$\sum _{{\alpha }_1\dots {\alpha }_n}{Z}_n^{\pm }(\{t_i,{\alpha }_i{\}}_{1\le i\le n})=0,$$

$$\begin{array}{rl}{Z}_n^{\pm }(\{t_i,{\alpha }_i{\}}_{1\le i\le n})& =(-i{\alpha }_1)\dots (-i{\alpha }_n)\times i^n{i}^n\\ & \times \prod _{\sigma }\det {\left[{\left({\hat{g}}_{\sigma }\right)}_{{\alpha }_i{\alpha }_j}(t_i,t_j)\right]}_{1\le i,j\le n}.\end{array}$$

$$\begin{gathered}\fmfgraph*(50,20)\fmfset{arrow_{l}en}{3mm}\fmfleft{i1}\fmfright{o1}\fmf{dashes_{a}rrow,tension=2,label.dist=0,label=\vspace{.4cm}\hskip 2.84544pt$0$}{i1,v1}\fmf{dashes_{a}rrow,tension=2,label.dist=0,label=\vspace{-.7cm}$t_{\text{max}}$}{v1,o1}\fmffreeze\fmf{phantom,label.dist=0,label=\vspace{.45cm}$1$}{v1,o1}\fmfv{decor.shape=circle,decor.filled=shaded,decor.size=4thick}{v1}\end{gathered}$$

$$\begin{array}{rl}{d}_n=& {\int }_0^{t_{\text{max}}}\mathrm{d}{t}_1\dots \mathrm{d}{t}_n\sum _{{\alpha }_1\dots {\alpha }_n}(-i{\alpha }_1)\dots (-i{\alpha }_n)\\ & \times \frac{i^{n+1}{i}^n}{i^2}\prod _{\sigma }\det {\mathcal{D}}_{\sigma }^{\pm }(\{t_i,{\alpha }_i{\}}_{1\le i\le n}),\end{array}$$

$${\mathcal{D}}_{\uparrow }^{\pm }(\{t_i,{\alpha }_i{\}}_{1\le i\le n})=\left(\begin{array}{cc}{\left[({\hat{g}}_{\uparrow }{)}_{{\alpha }_i{\alpha }_j}(t_i,t_j)\right]}_{1\le i,j\le n}& \begin{array}{c}({\hat{g}}_{\uparrow }{)}_{{\alpha }_{1}1}(t_1,t_{\text{max}})\\ ⋮\\ ({\hat{g}}_{\uparrow }{)}_{{\alpha }_{n}1}(t_n,t_{\text{max}})\end{array}\\ \begin{array}{ccc}({\hat{g}}_{\uparrow }{)}_{0{\alpha }_1}(t_{\text{max}},t_1)& \dots & ({\hat{g}}_{\uparrow }{)}_{0{\alpha }_n}(t_{\text{max}},t_n)\end{array}& ({\hat{g}}_{\uparrow }{)}_{01}(t_{\text{max}},t_{\text{max}})\end{array}\right),$$

$${\mathcal{D}}_{\downarrow }^{\pm }(\{t_i,{\alpha }_i{\}}_{1\le i\le n})={\left[({\hat{g}}_{\downarrow }{)}_{{\alpha }_i{\alpha }_j}(t_i,t_j)\right]}_{1\le i,j\le n}.$$

$${\hat{G}}_{\sigma }^{\text{LO}}(t,t^{\prime })\equiv L^{†}{\tau }_3{\hat{G}}_{\sigma }(t,t^{\prime })L,$$

$$R_{\sigma }(t,t^{\prime })$$

$$A_{\sigma }(t,t^{\prime })$$

$$K_{\sigma }(t,t^{\prime })$$

$$L^{†}{\tau }_3{m}_{+}L=\frac{1}{2}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)\equiv \frac{1}{2}{\tau }_{\uparrow },$$

$$L^{†}{\tau }_3{m}_{-}L=\frac{1}{2}\left(\begin{array}{cc}-1& 1\\ 1& -1\end{array}\right)\equiv \frac{1}{2}{\tau }_{\downarrow }.$$

$$-\frac{iU}{4}\left({\tau }_{\uparrow }\otimes {\tau }_{\uparrow }-{\tau }_{\downarrow }\otimes {\tau }_{\downarrow }\right)=-\frac{iU}{2}\left(1\otimes {\tau }_{\downarrow }+{\tau }_{\uparrow }\otimes 1\right),$$

$${\left({\hat{\tilde{g}}}_{\sigma }\right)}_{l{l}^{\prime }}$$

$${\left({\hat{\tilde{g}}}_{\sigma }\right)}_{l2}$$

$${\left({\hat{\tilde{g}}}_{\sigma }\right)}_{2l}$$

$${\left({\hat{\tilde{g}}}_{\sigma }\right)}_{22}$$

$${\hat{\tilde{g}}}_{\sigma }=\left(\begin{array}{ccc}{r}_{\sigma }& k_{\sigma }& r_{\sigma }+\sigma k_{\sigma }\\ 0& a_{\sigma }& \sigma a_{\sigma }\\ \sigma r_{\sigma }& \sigma k_{\sigma }+a_{\sigma }& \sigma [r_{\sigma }+a_{\sigma }]+k_{\sigma }\end{array}\right).$$

$$L^{\sigma }=$$

$$\begin{gathered}\fmfgraph*(40,20)\fmfset{arrow_{l}en}{3mm}\fmfleft{v1}\fmfright{v2}\fmf{plain_{a}rrow,label.dist=0,label=\vspace{.5cm}$\sigma$}{v1,v2}\fmflabel{$t', L'^\sigma$}{v1}\fmflabel{$t, L^\sigma$}{v2}\end{gathered}\hskip 25.6073pt=i\left(\hat{\tilde{g}}_{\sigma}\right)_{L^{\sigma}L^{\prime\sigma}}(t,t^{\prime}).$$

$$\begin{split}\begin{gathered}\fmfgraph*(40,30)\fmfset{arrow_{l}en}{3mm}\fmfleft{i1,i2}\fmfright{o1,o2}\fmf{plain_{a}rrow}{i1,v}\fmf{plain_{a}rrow}{v,o2}\fmf{plain_{a}rrow}{o1,v}\fmf{plain_{a}rrow}{v,i2}\fmfdot{v}\fmflabel{$\{t, i_\tau, l\}$}{v}\fmflabel{$L^\uparrow, \uparrow$}{i1}\fmflabel{$L^\downarrow, \downarrow$}{o1}\fmflabel{$L^\uparrow, \uparrow$}{i2}\fmflabel{$L^\downarrow, \downarrow$}{o2}\end{gathered}\hskip 22.76228pt=-\frac{iU}{2}&\left(\delta_{i_{\tau}1}\delta_{L^{\uparrow}2}\delta_{L^{\downarrow}l}\right.\\ &\left.\hskip 8.5359pt+\delta_{i_{\tau}-1}\delta_{L^{\downarrow}2}\delta_{L^{\uparrow}l}\right).\end{split}$$

Full text

I Introduction

The development of high-precision and controlled computational methods for non-equilibrium models in strongly-correlated regimes is a subject of growing interest in theoretical condensed-matter physics. Recent years have seen significant experimental progress with quantum transport through mesoscopic systems Potok et al. (2007), metal-insulator transitions driven by an electric field Nakamura et al. (2013) or light-induced superconductivity Fausti et al. (2011); Nicoletti et al. (2014); Casandruc et al. (2015); Nicoletti and Cavalleri (2016); Nicoletti et al. (2018).

Powerful tools have been designed for the study of quantum systems at equilibrium. Notably, the combination of dynamical mean-field theory Georges et al. (1996); Kotliar et al. (2006); Aoki et al. (2014) and state-of-the-art continuous-time Quantum Monte Carlo (QMC) algorithms such as CT-INT Rubtsov and Lichtenstein (2004); Rubtsov et al. (2005), CT-AUX Gull et al. (2008), or CT-HYB Werner et al. (2006); Werner and Millis (2006) have allowed for great advances. When considering out-of-equilibrium systems, however, early attempts to construct similar perturbation-expansion-based real-time QMC algorithms encountered an exponential sign problem that prevented them from reaching long times and large interactions Mühlbacher and Rabani (2008); Werner et al. (2009, 2010); Schiró and Fabrizio (2009); Schiró (2010). Other approaches such as the density matrix renormalization group (DMRG) White (1992, 1993); Schollwöck (2005) also struggle in the long-time limit due to entanglement growth. There is therefore still a great need for high-precision numerical methods that would be able to access the non-equilibrium steady states of strongly-interacting quantum systems.

Current efforts to build real-time quantum Monte Carlo methods mainly explore two routes: the inchworm algorithm Cohen et al. (2014a, b, 2015); Chen et al. (2017a, b); Antipov et al. (2017); Boag et al. (2018) and the so-called “diagrammatic” QMC Profumo et al. (2015); Bertrand et al. (2019a, b) which is the subject of this article. Using an expansion of physical quantities in powers of the interaction $U$, this algorithm has been shown to directly address the infinite-time steady states. The name “diagrammatic” refers to its imaginary-time counterparts that were historically constructing a Markov chain in the space of Feynman diagrams Prokof’ev and Svistunov (1998, 2007); Houcke et al. (2010); Bourovski et al. (2004).

First introduced in Ref. Profumo et al., 2015, the real-time diagrammatic QMC algorithm stochastically samples physical quantities using an importance sampling. At a given perturbation order $n$, its key idea is to regroup a factorial number of Feynman diagrams in a sum over Keldysh indices of $2^{n}$ determinants. This exponential sum has been shown to cancel vacuum diagrams, a property also used in recent diagrammatic QMC methods in imaginary-time Rossi (2017); Moutenet et al. (2018); Simkovic IV. and Kozik (2017). As a direct consequence, the Monte Carlo sampling only involves interaction times in a neighborhood around the measurement time $t_{\text{max}}$: we talk about the clusterization of times. The computation of the Monte Carlo weight is exponential in the perturbation order but uniform in time, at any temperature. The algorithm can therefore address long, even infinite, times in the computation of contributions to the perturbation theory. This method was recently generalized to compute the Green’s function and tested in quantum impurity models Bertrand et al. (2019a, b). The current form of the algorithm is able to compute the Kondo resonance at low temperature in the strongly-correlated Kondo regime.

Coefficients of the expansion being written in terms of high-dimensional integrals of the sum of determinants, its exponential scaling limits our capability to compute high orders with great precision (we typically are limited to 10 of them). Even though non-perturbative information and Bayesian techniques can overcome noise amplification occurring in the resummation of the series Bertrand et al. (2019b), this can prevent the algorithm to reach very large $U$.

In this article, we show that we can obtain the cancellation of diagrams and the long-time clusterization property without summing an exponential number of terms. Using the Larkin-Ovchinnikov (LO) basis in Keldysh space, we rewrite the integrand as a sum of $4^{n}$ determinants, but we show that diagrammatic rules in this basis are such that every diagram has the clusterization property. In other words, the elimination of vacuum diagrams is directly achieved in the diagrammatics without the need of an exponential sum. We then implement and compare two Monte-Carlo algorithms based on this mathematical property. Both sample single determinants at a polynomial cost, but then one measures in the LO basis (LO algorithm) while the other measures in the original basis (mixed algorithm). We obtain that a simple implementation of the real-time diagrammatic QMC in the Larkin-Ovchinnikov basis leads to a severe sign problem, which is reduced in the mixed algorithm. This shows that the main effect of the exponential sum of determinants, beyond the cancellation of vacuum disconnected diagrams, is to reduce the sign problem of this class of algorithms.

This article is organized as follows. First, we present in Sec. II the usual Keldysh formalism in the $\{\pm\}$ basis, briefly summarize the diagrammatic rules and then derive the cancellation of vacuum diagrams and the clusterization of the density when summing over Keldysh indices. We follow the same structure in Sec. III where we introduce the Larkin-Ovchinnikov basis, showing that all vacuum diagrams are equal to zero, so that density contributions directly clusterize around the measurement time. We then detail in Sec. IV the Monte Carlo implementation of the original algorithm presented in Ref. Profumo et al., 2015 ($\pm$ algorithm) and two algorithms based on the Larkin-Ovchinnikov formalism (LO and mixed algorithms). In Sec. V we compute the density of an impurity level coupled to a bath, present the results of all three algorithms and explain the origin of the observed error bars. We finally conclude in Sec. VI.

II Keldysh formalism

We work in the Keldysh formalism Schwinger (1961); Keldysh (1964); Rammer and Smith (1986); Kamenev and Levchenko (2009). In this framework, operators act on the Keldysh contour $\mathcal{C}$ consisting of a forward branch, from an initial time $t_{0}$ (that we take equal to [math] in the following) to a given time $t_{\text{max}}$, and a backward branch, from $t_{\text{max}}$ to $t_{0}$. The system is initially prepared at equilibrium without interactions. A Keldysh point $k$ on $\mathcal{C}$ is defined as a pair $k\equiv(t,\alpha)$ with a time $t\in[0,t_{\text{max}}]$ and a Keldysh index $\alpha\in\{\pm\}$ indicating which branch is to be considered. The + (resp. -) index denotes the forward (resp. backward) branch, as depicted below.

{fmffile}

keldysh_contour {fmfgraph*}(200,10) \fmfsetarrow_len3mm \fmflefti1,i2,i3 \fmfrighto1,o2,o3 \fmfphantomi2,v1 \fmfphantom,label.dist=0,tension=15.v1,o2 \fmfphantomv1,v2 \fmfphantom,label.dist=0,label= $t_{\text{max}}$v2,o2 \fmfplain_arrow, label.dist=0, label= **-**o1,i1 \fmfplain_arrow, label.dist=0, label=             **+**i3,o3 \fmfphantom,tag=1, label.dist=0, label= [math]i1,i3 \fmfipathp[] \fmfisetp1vpath1(__i1,__i3) \fmfiplainsubpath (length(p1)/3, 2length(p1)/3) of p1 \fmffreeze\fmfplaini2,v2 \fmfphantom_arrowv1,o2 \fmfplain,righto1,o3

Note that both branches are along the real axis and are displaced only for graphical purposes. In the following, Greek letters refer to $\pm$ indices unless otherwise stated. We define a contour operator $T_{\mathcal{C}}$ that follows the arrows on the above picture: $T_{\mathcal{C}}$ coincides with the usual time-ordering operator T on the $+$ branch, with the anti-time ordered operator Ť on the $-$ branch, and considers all Keldysh points on the backward branch to be later than points on the forward branch.

The formalism we develop in this section is valid for any general model described by a noninteracting Green’s function $g$ and a density-density interaction. However, for the sake of simplicity, we consider interacting electrons on a single energy level. The operator $c_{\sigma}$ (resp. $c_{\sigma}^{\dagger}$) destroys (resp. creates) an electron with spin $\sigma=\uparrow,\downarrow$. The interaction term, turned on at $t=0$, is given by the interaction vertex $Un_{\uparrow}n_{\downarrow}$, where $n_{\sigma}\equiv c_{\sigma}^{\dagger}c_{\sigma}$ is the density operator.

We define the time-dependent Green’s function

[TABLE]

where $c_{\sigma}^{(\dagger)}(t)$ is the Heisenberg representation of $c_{\sigma}^{(\dagger)}$ and the average is taken with respect to the initial noninteracting state. The Green’s function takes the form of a $2\times 2$ matrix in the $\{\pm\}$ basis: $\hat{G}_{\sigma}=\begin{pmatrix}G_{\sigma}^{++}&G_{\sigma}^{<}\\ G_{\sigma}^{>}&G_{\sigma}^{--}\end{pmatrix}$, where

[TABLE]

Throughout the article, noninteracting Green’s functions will be denoted by lower case letters, interacting ones by upper case letters, and a $~{}\hat{}~{}$ denotes a matrix.

II.1 Diagrammatic rules

In this article, we construct perturbation series in the interaction $U$ for physical observables of interest. Computing contributions at different perturbation orders relies on the evaluation of Feynman diagrams obeying rules that we briefly summarize.

A straight line represents a noninteracting Green’s function {fmffile}keldysh_g0

[TABLE]

Because the interaction has the form $Un_{\uparrow}n_{\downarrow}$, an interaction vertex is characterized by a single Keldysh point $\{t,\alpha\}$, and the indices of the four legs all have to be equal to the Keldysh index $\alpha$

{fmffile}keldysh_vertex

[TABLE]

Hence, for every interaction time $t$, there are two possible vertices. The sum of the different $\{\pm\}$ configurations can be written in the $\mathcal{H}_{\uparrow}\otimes\mathcal{H}_{\downarrow}$ space, in the form

[TABLE]

where $m_{+}=\begin{pmatrix}1&0\\ 0&0\end{pmatrix}$ and $m_{-}=\begin{pmatrix}0&0\\ 0&1\end{pmatrix}$ are matrices in the $\{\pm\}$ basis, and $\mathcal{H}_{\sigma}$ is the Hilbert space for spin $\sigma$. Furthermore, an interaction of the form $hc_{\sigma}^{\dagger}c_{\sigma}$ in the Hamiltonian would give rise to 2-leg vertices of the form {fmffile}keldysh_2_leg

[TABLE]

These do not appear directly in the diagrammatics but will be formally useful when deriving the expression of the fermionic bubble. The sum over Keldysh indices reads $-ih(m_{+}-m_{-})$ in both $\mathcal{H}_{\uparrow}$ and $\mathcal{H}_{\downarrow}$ spaces.

With the expression of the 4-leg interaction vertex, the following fermionic bubble reads {fmffile}keldysh_closed_loop

[TABLE]

Because of the form of the interaction term, we have

[TABLE]

Hence, the above diagram reduces to $\alpha Ug_{\bar{\sigma}}^{<}(t,t)$, which can be formulated as a 2-leg vertex with a $iUg_{\bar{\sigma}}^{<}(t,t)$ field {fmffile}keldysh_closed_loop2

[TABLE]

If $\mathcal{M}$ is the quantity we want to compute (later on the density), its perturbation expansion is given by $\mathcal{M}=\sum_{n}\mathcal{M}_{n}U^{n}$. Because of the form of the interaction vertex, we have

[TABLE]

where $\mathcal{M}_{n}^{\pm}(\{t_{i},\alpha_{i}\}_{1\leq i\leq n})$ can be expressed as a product of determinants, their precise form depending on the measured quantity. Throughout this paper, the $\pm$ superscript will denote quantities expressed in the $\{\pm\}$ basis. Moreover times integrated over are always considered ordered.

II.2 Cancellation of vacuum diagrams when summing over Keldysh indices

Due to the forward-backward nature of the contour $\mathcal{C}$, the partition function $Z$ is exactly equal to 1 in the real-time Keldysh formalism. Expressing $Z$ as a series in $U$ ($Z=\sum_{n}Z_{n}U^{n}$), this property implies that all $Z_{n}$ are vanishing for $n\geq 1$. Because of the form of Eq. (10), this cancellation involves both the integral over times and the sum over Keldysh indices. However, it was proven by Profumo and co-workers in Ref. Profumo et al., 2015 that only the latter is needed. For all $n\geq 1$, $\{t_{1},\dots,t_{n}\}\in[0,t_{\text{max}}]^{n}$,

[TABLE]

where

[TABLE]

Each $(-i\alpha_{k})$ comes from Eq. (4), and the two $i^{n}$ factors from the fact that a straight line actually represents an $i\hat{g}$ (Eq. (3)).

For every configuration of times $\{t_{1},\dots,t_{n}\}$, vacuum diagrams therefore cancel when performing the explicit $2^{n}$ sum over Keldysh indices. Recent developments in imaginary-time diagrammatic QMC also achieved, through an iterative procedure, the cancellation of vacuum (and, later on, non one-particle irreducible) diagrams at every Monte Carlo step at an exponential cost in the perturbation order. Rossi (2017); Moutenet et al. (2018); Simkovic IV. and Kozik (2017)

II.3 Density computation and clusterization

In the following, we compute the density $d$ of electrons with spin $\uparrow$ on the impurity level at the end point of the Keldysh contour, $d\equiv\langle n_{\uparrow}(t_{\text{max}})\rangle$.

In the $\{\pm\}$ basis, let us note that $d=(\hat{G}_{\uparrow})_{01}(t_{\text{max}},t_{\text{max}})/i$. Hence we can represent the measurement vertex as a “special” vertex bearing time $t_{\text{max}}$, such that the ingoing and outgoing Keldysh indices are 0 and 1: {fmffile}keldysh_measurement

[TABLE]

Note that the surrounding lines are dashed because they should bear a $\hat{g}$ propagator (instead of an $i\hat{g}$ one as in the rest of the formalism). The order-$n$ contribution to $d$ reads

[TABLE]

where

[TABLE]

and

[TABLE]

Using the cancellation of vacuum diagrams when summing over Keldysh indices, we reproduce in Appendix A the argument of Ref. Profumo et al., 2015 showing that the computation of $d_{n}$ only involves the sampling of interaction times close to $t_{\text{max}}$. As a direct consequence, Monte Carlo algorithms implementing this sum in the calculation of the weight can address any measurement time $t_{\text{max}}$, when earlier methods were limited to short-term measurements Mühlbacher and Rabani (2008); Werner et al. (2009, 2010); Schiró and Fabrizio (2009); Schiró (2010). We talk about the clusterization of interaction times in the computation of the density.

III Larkin-Ovchinnikov Formalism

Starting from the expression of the Green’s function in the $\{\pm\}$ basis, we define its counterpart in the LO basis, $\hat{G}^{\text{LO}}$, through the following transformation Keldysh (1964); Larkin and Ovchinnikov (1986)

[TABLE]

where $L=\frac{1}{\sqrt{2}}\begin{pmatrix}1&1\\ -1&1\end{pmatrix}$ and $\tau_{3}=\begin{pmatrix}1&0\\ 0&-1\end{pmatrix}$. The Green’s function now takes the $2\times 2$ form $\hat{G}_{\sigma}^{\text{LO}}=\begin{pmatrix}R_{\sigma}&K_{\sigma}\\ 0&A_{\sigma}\end{pmatrix}$, where $R$, $K$, and $A$ are respectively the retarded, Keldysh and advanced Green’s functions defined as

[TABLE]

In this basis, the Keldysh index $\alpha\in\{\pm\}$ is replaced by an LO index 0 or 1. In the following, $l$ will always denote such an index unless otherwise stated.

III.1 Diagrammatic rules

To expose the diagrammatic rules in this formalism, let us first determine from Eq. (5) the form of the 4-leg interaction vertex in the LO basis. The $m_{+}$ and $m_{-}$ matrices transform as

[TABLE]

Hence the sum of different LO contributions can be written

[TABLE]

where $\mathbf{1}$ is the $2\times 2$ identity matrix. Note that this is consistent with the symmetric form $-\frac{iU}{2}(\mathbf{1}\otimes\tau+\tau\otimes\mathbf{1})$ noted in Ref. Biroli and Parcollet, 2002, where $\tau=\begin{pmatrix}0&1\\ 1&0\end{pmatrix}$. The rhs form of Eq. (21) is the one we will retain in the rest of this article. We will show in Sec. III.2 and III.3 that the identity part of the vertex is essential in the proof of the cancellation of vacuum diagrams and the clusterization of times in the computation of observables.

The key point of this expression of the vertex is that we can reduce the number of indices involved in the diagrammatics using the fact that $\tau_{\uparrow}$ and $\tau_{\downarrow}$ are rank-1 matrices: $\tau_{\uparrow}=v_{\uparrow}v_{\uparrow}^{\top}$ with $v_{\uparrow}=\begin{pmatrix}1\\ 1\end{pmatrix}$ and $\tau_{\downarrow}=v_{\downarrow}(-v_{\downarrow}^{\top})$ with $v_{\downarrow}=\begin{pmatrix}1\\ -1\end{pmatrix}$. We can therefore absorb the $\tau_{\sigma}$ part of the vertex in a redifinition of the noninteracting propagator (see below).

An LO vertex can then be characterized by a tuple $\{t,i_{\tau},l\}$, where $t\in[0,t_{\text{max}}]$, $i_{\tau}\in\{-1,1\}$ and $l\in\{0,1\}$. $i_{\tau}=1$ (resp. -1) indicates that the $\uparrow$ (resp. $\downarrow$) spin is carrying the $\tau_{\uparrow}$ (resp. $\tau_{\downarrow}$) side, and $l$ is the LO index entering the identity-part of the vertex. We store the information about both the bare propagator $\hat{g}_{\sigma}^{\mathrm{LO}}=\begin{pmatrix}r_{\sigma}&k_{\sigma}\\ 0&a_{\sigma}\end{pmatrix}$ and the nature of the vertices it is connected to in the form of a $3\times 3$ matrix $\hat{\tilde{g}}_{\sigma}$. The two first indices corresponds to a connection to the identity (with $l=$ 0 or 1), and the third one to the connection to a $\tau_{\sigma}$:

[TABLE]

with the convention that $\sigma=\uparrow$ should be understood as $+1$ and $\sigma=\downarrow$ as $-1$.

We obtain

[TABLE]

To simplify upcoming equations, we express the indices of $\hat{\tilde{g}}_{\uparrow}$ and $\hat{\tilde{g}}_{\downarrow}$ at a vertex $\{t,i_{\tau},l\}$ in the form of two composite indices $L^{\uparrow}$ and $L^{\downarrow}$:

[TABLE]

Note that this $3\times 3$ form of the Green’s function comes from the absorption of the $\tau_{\sigma}$ part of the vertex and has nothing to do with the Baym-Kadanoff $L$-shaped contour used in thermal real-time computations.

With this notation, a straight line represents a noninteracting (modified) Green’s function {fmffile}lo_g0

[TABLE]

As discussed previously, the interaction vertex, proportional to the identity in the $\{\pm\}$ basis, is now proportional to $\frac{\mathbf{1}\otimes\tau_{\downarrow}+\tau_{\uparrow}\otimes\mathbf{1}}{2}$ in the $\mathcal{H}_{\uparrow}\otimes\mathcal{H}_{\downarrow}$ space

{fmffile}lo_vertex

[TABLE]

As $m_{+}-m_{-}$ transforms into the $2\times 2$ identity matrix in the LO basis, a 2-leg vertex is simply characterized by an interaction time $t$ and an LO index $l$. A term $hc_{\sigma}^{\dagger}c_{\sigma}$ in the Hamiltonian would therefore give rise to the following vertex {fmffile}lo_2_leg

[TABLE]

With this expression of the interaction vertex, the following fermionic bubble     {fmffile}lo_closed_loop {fmfgraph*}(50,20) \fmfsetarrow_len3mm \fmflefti1,i2 \fmfrighto1,o2 \fmfplain_arrow,tension=2,label.dist=0, label= $\sigma$i1,v1 \fmfplain_arrow,tension=2v1,o1 \fmfphantomi2,v2,o2 \fmflabel$\{t,i_{\tau},l\}$v1 \fmflabel$L^{\sigma}$i1 \fmflabel$L^{\sigma}$o1 \fmfdotv1 \fmffreeze\fmfplain_arrow,right,tension=2v1,v2,v1 \fmfphantom,label.dist=0, label=$\bar{\sigma}$v1,v2

    evaluates to

[TABLE]

For the equal-time limit of the retarded, Keldysh and advanced Green’s function, we choose a convention which ensures the consistency between the $\{\pm\}$ and LO basis. We consider

[TABLE]

and we show that this is consistent with Eq. (8). Using Eq. (29), the above fermionic bubble reduces to $Ug_{\bar{\sigma}}^{<}(t,t)\delta_{L^{\sigma}l}$. It can be rewritten as a 2-leg vertex with a $iUg_{\bar{\sigma}}^{<}(t,t)$ field {fmffile}closed_vertex_lo

[TABLE]

This equation is, up to a change of basis, the same as Eq. (9). The choice of equal-time limit described in Eq. (29) is therefore consistent with the $\{\pm\}$ basis formalism.

The order-$n$ contribution to the quantity $\mathcal{M}$ we want to measure in an expansion in $U$ is similar to Eq. (10), but has to take account of the new form of the vertex

[TABLE]

where the $\mathcal{M}_{n}^{\text{LO}}(\{t_{i},i_{\tau_{i}},l_{i}\}_{1\leq i\leq n})$ can once again be expressed as a product of determinants, their precise form depending on the computed quantity.

This formalism leads to $4^{n}$ LO configurations for a given set of $n$ interaction times, to be compared with $2^{n}$ possible configurations in the $\pm$ basis. However, we show in the next section that vacuum diagrams now directly cancel in this formalism, without the actual need to perform an explicit sum over all configurations.

III.2 Cancellation of vacuum diagrams

In this section, we show the main result of this article: contributions to the partition function are directly equal to zero in the LO basis. For all $n\geq 1$, $\{t_{1},\dots,t_{n}\}\in[0,t_{\text{max}}]^{n}$, $\{i_{\tau_{1}},\dots,i_{\tau_{n}}\}\in\{-1,1\}^{n}$, $\{l_{1},\dots,l_{n}\}\in\{0,1\}^{n}$,

[TABLE]

where the contributions to the partition function are

[TABLE]

Each $-\frac{i}{2}$ comes from Eq. (26) and the two $i^{n}$ factors from the fact that a straight line actually represents an $i\hat{\tilde{g}}$ (Eq. (25)).

Let us consider an order $n\geq 1$ diagram contributing to $Z$. The interaction times are denoted $t_{1},\dots,t_{n}$. We introduce $\hat{t}=\mathrm{max}_{i}t_{i}$ and $\hat{i}$ such that $t_{\hat{i}}=\hat{t}$. We label $\sigma$ the spin on the identity side of the $(\mathbf{1}\otimes\tau_{\downarrow}+\tau_{\uparrow}\otimes\mathbf{1})/2$ interaction vertex at $\hat{t}$, and $l$ the corresponding LO index. In the consider, we consider the diagrammatic line following spin $\sigma$. If $\hat{t}$ is surrounded by no other interaction vertex, the diagram is then proportional to

[TABLE]

In the case where $\hat{t}$ is surrounded by at least one other interaction vertex, we label its surrounding interaction times (that can be equal) $t_{i}$ and $t_{j}$, $i,j\neq\hat{i}$, with corresponding composite indices $L_{i}^{\sigma}$, $L_{j}^{\sigma}$. We then obtain

[TABLE]

and

[TABLE]

The full diagram is then proportional to $\delta_{l0}\delta_{l1}=0$. Hence every diagram contributing to $Z$ in the LO basis is exactly equal to 0. This formalism directly cancels vacuum diagrams.

Finally, we note that this proof relies only on having the identity on one side of the interaction vertex, and not on the explicit contraction with $\tau_{\uparrow},\tau_{\downarrow}$. Had we kept the diagrammatics with $\hat{g}$ lines instead of $\hat{\tilde{g}}$ ones, we would also obtain the cancellation of vacuum diagrams.

III.3 Density computation and clusterization

In order to understand how to write the density of $\uparrow$ electrons on the energy level in the LO basis, we use the following property of the Keldysh formalism: the average value of an operator does not depend on the branch of $\mathcal{C}$ where it is computed. Considering $d$ on the $+$ branch of the contour, the computation of the density can be understood as the action of the $m_{+}$ matrix in the $\{\pm\}$ basis, which transforms in the $\frac{1}{2}\tau_{\uparrow}$ matrix in the LO basis according to Eq. (20a). Hence we can represent the measurement vertex as a “special” interaction vertex at time $t_{\text{max}}$ with $i_{\tau}=1$: {fmffile}lo_measurement

[TABLE]

As previously, surrounding lines are dashed because they bear a $\hat{\tilde{g}}$ (and not an $i\hat{\tilde{g}}$). Hence the order-$n$ contribution to $d$ reads

[TABLE]

The $\mathcal{D}_{\sigma}^{\text{LO}}$ matrices are defined as

[TABLE]

and

[TABLE]

Before considering the clusterization of interaction times, we note that half of the contributions to the density vanish. Let us consider a given set $\{t_{i},i_{\tau_{i}},l_{i}\}_{1\leq i\leq n}$ of LO vertices at order $n$, and let us label $\hat{t}=\max_{i}t_{i}$ and $\hat{i}$ such that $t_{\hat{i}}=\hat{t}$. If $i_{\tau_{\hat{i}}}=1$, then the $\downarrow$ spin is carrying the identity side of the vertex. As we measure the density on the $\uparrow$ spin, the argument used in the cancellation vacuum diagrams (see III.2) applies again and $\mathcal{D}_{\downarrow}^{\text{LO}}(\{t_{i},i_{\tau_{i}},l_{i}\}_{1\leq i\leq n})$ is the $n\times n$ null matrix. If $i_{\tau_{\hat{i}}}=-1$, the contribution does not vanish. Hence, when computing the density, at every order $n$ and for every set of $n$ interaction times, $4^{n}/2$ LO configurations (out of $4^{n}$) are exactly zero.

The clusterization of interaction times around $t_{\text{max}}$ in the calculation of the density is then a direct consequence of the cancellation of vacuum diagrams and is very similar to the proof in the $\{\pm\}$ basis (now without the exponential sum). Let $n$ be a given perturbation order, and $t_{1}<t_{2}<\dots<t_{n}$ $n$ interaction times. Let us assume that the first $j$ times are located far away from the measurement time $t_{\text{max}}$, and that the last $n-j$ times are located in the vicinity of $t_{\text{max}}$. We can formally consider

[TABLE]

Because the Green’s function is a local quantity in time, this means that for all $t\in\{t_{1},\dots,t_{j}\}$, $t^{\prime}\in\{t_{j+1},\dots,t_{n};t_{\text{max}}\}$

[TABLE]

We therefore have

[TABLE]

with

[TABLE]

and $B_{\uparrow}$ is the $\left[(\hat{\tilde{g}})_{L_{i}^{\uparrow}L_{i^{\prime}}^{\uparrow}}(t_{i},t_{i^{\prime}})\right]_{j+1\leq i,i^{\prime}\leq n}$ matrix where a last line and column corresponding to $t_{\text{max}}$ are added, similar to Eq. (39). However, $\prod_{\sigma}A_{\sigma}$ is a contribution to $Z$ at order $j$, and it vanishes according to (32). Therefore $\prod_{\sigma}\det D_{\sigma}\simeq 0$, and this proves the clusterization of times around $t_{\text{max}}$ in the computation of the density.

In the next section, we present different algorithms to stochastically sample Eqs (15) and (38).

IV Monte Carlo implementation

In this section, we describe how to compute the density $d$ introduced above using quantum Monte Carlo (MC) techniques. We present three different algorithms to compute this quantity, one using the $\pm$ algorithm presented in Ref. Profumo et al., 2015 and the other two based on the LO formalism presented above.

IV.1 Monte Carlo algorithms

We first describe how to stochastically generate MC configurations to sample the order-$n$ contribution, $d_{n}$, as expressed in Eqs (15) and (38).

The $\pm$ algorithm works directly on the Keldysh contour. A configuration $\mathfrak{c}$ is determined by a given perturbation order $n$ and a set of $n$ interaction times (and not Keldysh points): $\mathfrak{c}=\{n;t_{1},\dots,t_{n}\}$. The contribution to $d_{n}$ of a given configuration is

[TABLE]

In the Monte Carlo, configurations are sampled stochastically according to their weight, which we choose to be $|w_{\mathfrak{c}}^{\pm}|$. We then have

[TABLE]

Note that it was shown in Ref. Profumo et al., 2015 that $w_{\mathfrak{c}}^{\pm}\in\mathbb{R}$.

In the LO algorithm, a configuration $\mathfrak{c}$ is now determined by a given perturbation order $n$ and a set of $n$ interaction LO vertices: $\mathfrak{c}=\{n;y_{1},\dots,y_{n}\}$, where $y_{i}=\{t_{i},i_{\tau_{i}},l_{i}\}$. Because the density is a real quantity, the contributions to $d_{n}$ of a configuration $\mathfrak{c}$ can be written as

[TABLE]

If $|w_{\mathfrak{c}}^{\text{LO}}|$ is the statistical weight of $\mathfrak{c}$ in the Monte Carlo process, then

[TABLE]

The third algorithm that we study is a mixed algorithm that samples the configurations according to their LO weight $|w_{\mathfrak{c}}^{\text{LO}}|$ but computes $d_{n}$ in the original $\{\pm\}$ basis, from the contributions $w_{\mathfrak{c}}^{\pm}$ at the sampled times. A configuration $\mathfrak{c}$ is then determined by a given perturbation order $n$ and a set of $n$ interaction LO vertices: $\mathfrak{c}=\{n;y_{1},\dots,y_{n}\}$ and the MC weight is $|w_{\mathfrak{c}}^{\text{mixed}}|=|w_{\mathfrak{c}}^{\text{LO}}|$, so that

[TABLE]

where $\mathcal{N}$ is the number of non-zero LO configurations. When computing the density, $\mathcal{N}=4^{n}/2$ at order $n$ (see Section III.3).

In all three techniques, we use a standard Metropolis algorithm Metropolis et al. (1953) to generate Markov chains distributed according to the weights $|w_{\mathfrak{c}}|$. Starting from a given configuration $\mathfrak{c}$, a new configuration $\mathfrak{c}^{\prime}$ is proposed according to one of the following two Monte Carlo updates:

Remove a randomly chosen interaction time (for the $\pm$ algorithm) or interaction LO vertex (for the LO and mixed algorithms) from $\mathfrak{c}$. 2. 2.

Add a new interaction time (for the $\pm$ algorithm) or an interaction LO vertex (for the LO and mixed algorithms). In all three techniques, because of the clusterization of times around $t_{\text{max}}$, we choose the new interaction time according to a Cauchy law (see below). We randomly choose the $i_{\tau}$ and $l$ indices.

The new configuration $\mathfrak{c}^{\prime}$ is accepted or rejected with the usual Metropolis ratio

[TABLE]

where $T_{\mathfrak{c}\mathfrak{c}^{\prime}}$ is the probability to propose $\mathfrak{c}^{\prime}$ after $\mathfrak{c}$.

IV.2 Proposition of times

We have shown previously that times clusterize around $t_{\text{max}}$. It is therefore more efficient to propose times located around it compared to uniformly distributed between [math] and $t_{\text{max}}$. We consider a Cauchy law determined by two parameters $t_{0}$ and $a$

[TABLE]

$C$ is a normalization factor such that the integral of $\rho$ between [math] and $t_{\text{max}}$ gives 1, defined as $C=a\left[C_{2}-C_{1}\right]$, where $C_{1}=\arctan\left(-\frac{t_{0}}{a}\right)$ and $C_{2}=\arctan\left(\frac{t_{\text{max}}-t_{0}}{a}\right)$.

To obtain a new time that follows this probability law, one can perform these three steps:

Choose a random number $u$ uniformly distributed between 0 and 1. 2. 2.

Construct

[TABLE]

uniformly distributed between $\frac{1}{2}+\frac{1}{\pi}C_{1}$ and $\frac{1}{2}+\frac{1}{\pi}C_{2}$. 3. 3.

Compute

[TABLE]

distributed between [math] and $t_{\text{max}}$ according to $\rho$.

The parameters $t_{0}$ and $a$ are then fitted to the 1D projection of times visited by the Monte Carlo, accumulated during the first part of the computation.

IV.3 Redefinition of noninteracting propagators

As shown in previous worksRubtsov et al. (2005); Wu et al. (2017); Profumo et al. (2015); Rossi et al. (2016), there is some freedom in the choice of the noninteracting propagator used to construct the perturbation expansion, since the interaction can be redefined as

[TABLE]

Note that in this subsection $\alpha$ does not denote a Keldysh index but a scalar, in order to be consistent with the existing literature. In particular, it was shown that $\alpha$ can strongly modify the radius of convergence of the perturbation series Profumo et al. (2015); Wu et al. (2017). This redefinition of the interaction term in Eq. (54) is taken into account by subtracting $\alpha$ on the diagonal of the determinants as explained and proved in Ref. Profumo et al., 2015. The second term in Eq. (54) acts as a shift in the chemical potential and can be absorbed in a redefinition of the noninteracting propagators.

Let us first consider the LO basis. This shift acts a diagonal term in the self-energy and hence in

[TABLE]

$\alpha$ therefore modifies $r_{\sigma}$ and $a_{\sigma}$ into

[TABLE]

As $k_{\sigma}/|r_{\sigma}|^{2}$ is not impacted by the shift, the modified Keldysh Green’s function is then

[TABLE]

From these expressions, we can then deduce the modified Green’s functions in the $\{\pm\}$ basis through a change of basis tranformation.

IV.4 Normalization procedure

All Monte Carlo algorithms presented above compute the order-$n$ contribution to the density $d$, however the MC results need to be normalized. Hence we restrict our calculation to two consecutive orders, $n$ and $n+1$, and a time or vertex can be added (resp. removed) only if the current configuration $\mathfrak{c}$ is at order $n$ (resp $n+1$). We measure both the density ($d_{n}$ and $d_{n+1}$) and a normalization factor ($\eta_{n}$ and $\eta_{n+1}$). In all algorithms, the normalization factor is chosen to be the sum of the absolute value of the contributions to the density:

[TABLE]

where the proportionality constant is the same as in the calculation of $d_{n}$. If $\tilde{d}_{n}$ and $\tilde{\eta}_{n}$ are the unrenormalized sums of the contributions accumulated in the Monte Carlo procedure, then the normalized values for $d_{n}$ and $\eta_{n}$ are obtained as

[TABLE]

and $\eta_{n}$ is then used to normalize the following simulation between orders $n+1$ and $n+2$. The lowest order is computed analytically to close the equations.

V Results

V.1 Density

In this section, we present actual computations of the density according to the algorithms described in the previous section and compare their efficiency. In the following, we consider an energy level $\epsilon_{d}$ coupled to a bath described by a semi-circular density of states of bandwidth $4t$. The Green’s function describing this bath is defined on the complex plane as Georges et al. (1996)

[TABLE]

The noninteracting retarded Green’s function of the impurity level is

[TABLE]

where $\gamma$ is a coupling term between the energy level and the bath. The Keldysh Green’s function is then deduced using the fluctuation-dissipation theorem

[TABLE]

In the following, $t=1$ is our energy unit. We consider $\beta t=100$, $\gamma^{2}=0.04t^{2}$, $\epsilon_{d}=-0.36t$. Electrons on the impurity experience a local Coulomb interaction $U=1.2t$. We choose the $\alpha$ shift to be $\alpha=0.3$ (see Sec. IV.3), such that $U\alpha=-\epsilon_{d}$. The bath being particle-hole symmetric, this creates a shifted retarded Green’s function $\bar{r}(\omega)$ that is itself particle-hole symmetric (see Eq. (56)). However, we have checked that this particular choice of $\alpha$ does not influence our conclusions. We provide in Appendix B a table benchmarking the LO and mixed algorithms against the original $\pm$ algorithm. This shows in particular that the LO and mixed algorithms yield correct results and that we can indeed reach long times in the LO algorithm without an exponential sum of determinants.

Our main result is shown on Figure 1 where we compare the relative error bar in the density computation as a function of the perturbation order. Blue dots denote the $\pm$ algorithm, orange stars the LO algorithm, and green dots the mixed algorithm. The order-9 relative error is not shown for the LO algorithm as it exceeds 1 and is therefore meaningless. In all three cases, dotted lines are guides to the eye. The computational time is 240 CPU*hours for each order.

We see that all three relative error bars increase with perturbation order. This can either come from the increasing difficulty of computing the series coefficients, or an error propagation coming from the normalization factor $\eta$. We plot in Appendix C the relative error bar on $\eta$, which is much smaller than the final relative error on the density, showing that the latter mainly comes from the increasing difficulty to compute higher order coefficients. Moreover, the LO relative error bars very quickly become much larger than the $\pm$ ones, their difference nearly reaching two orders of magnitude at order $8$. The mixed algorithm is found to perform better than the LO algorithm but its error bars slowly grow larger than the $\pm$ ones. This is surprising, as one could have expected to at least gain the decorrelation time over the algorithm of Ref. Profumo et al., 2015. We discuss the origin of the error bars in both algorithms in the next section.

V.2 The return of the sign problem

In this section, we discuss the origin of the large variance in the computation of the density in the LO algorithm in terms of a sign problem in the Monte Carlo sampling and we show how this impacts the error bars of the mixed algorithm.

On the upper panel of Figure 2, we plot as blue dots the non-zero LO weights for two different time configurations, sorted according to their absolute value. The left and right panel correspond to two different time configurations (Cf caption). In both cases, the red line indicates the full sum over all LO indices, normalized to 1 (which coincides with the $\pm$ weights). The lower panel shows the partial sum, from left to right, of the LO weights plotted above. The last point, equal to 1 by construction, is emphasized as a red dot. As roughly half of the weights are positive and half negative, we see that the sum of the LO weights over the indices at fixed time configuration is characterized by a massive cancellation. This is the origin of the large error bar in the Monte-Carlo, i.e. another manifestation of the sign problem. Furthermore, the partial sum shows that there is no clear feature or cutoff from which one could extract the value of the full sum.

Let us now turn to the mixed algorithm. On both the left and right panels of Figure 2, the sum over all LO indices, which coincides with the $\pm$ weight, is normalized to 1. However, on the left panel, the weights of the different LO configurations are small compared to the final result, reaching at most 20% of it. On the right panel, those same weights are much bigger, reaching up to 1700% of the full sum. Hence the Monte Carlo implemented in the LO basis does not sample the same time configurations as the algorithm in the $\{\pm\}$ basis. This is illustrated in Figure 3 where the histograms of the times visited by the Monte Carlo, projected in one dimension, are plotted for both the $\pm$ algorithm (blue line) and LO one (orange line). First, we observe the clusterization of times proved at the beginning of this article: interaction times contributing to the density tend to be in the vicinity of $t_{\text{max}}$. Then, we see that some times located far away from the measurement but still contributing significantly to the $\pm$ algorithm are almost never visited in the LO algorithm. On the other hand, times close to $t_{\text{max}}$ are more sampled in the latter. As times visited by the mixed algorithm coincide with the LO ones, this explains the difference in error bars between the mixed and $\pm$ algorithms observed in Figure 1.

VI Conclusion

In conclusion, the explicit sum over the Keldysh indices of the original $\pm$ algorithm of Ref. Profumo et al., 2015 has two functions: i) it allows to reach the very long times due to the clusterization of the integrand caused by the cancellation of vacuum diagrams; *ii) * it strongly reduces the error bar by performing a massive cancellation of terms. In this article, we have shown that one can obtain the first properties for each determinant using the Larkin-Ovchinnikov basis, hence without the exponentially large sum of determinants. A direct implementation of the algorithm in the LO basis indeed reaches the steady state, but also has an error bar growing quickly with the order $n$ due to a sign problem. An interesting possibility would be the existence of an optimum between the LO and original $\pm$ algorithms, using partial groupings of terms in the LO basis with less than $2^{n}$ terms that would reduce the sign problem and yields a better scaling than the original algorithm in the $\{\pm\}$ basis. Work is in progress in this direction.

Acknowledgements.

We are grateful to Xavier Waintal and Antoine Georges for useful discussions. This work was partly supported by the European Research Council grant ERC-278472-MottMetals (PS, OP). The Flatiron Institute is a division of the Simons Foundation. Part of this work was performed using HPC resources from GENCI (Grant No. A0050510609).

Appendix A Clusterization of the density in the $\{\pm\}$ basis

We reproduce here the argument of Ref. Profumo et al., 2015 showing that the cancellation of vacuum diagrams when summing over Keldysh indices implies the clusterization of interaction times near $t_{\text{max}}$.

Let $n$ be a given perturbation order, and $t_{1}<t_{2}<\dots<t_{n}$ $n$ interaction times. Let’s assume that the first $j$ times are located far away from the measurement time $t_{\text{max}}$, and that the last $n-j$ times are located in the vicinity of $t_{\text{max}}$. We can formally consider

[TABLE]

Because the Green’s function is a local quantity, this means that for all $t\in\{t_{1},\dots,t_{j}\}$, $t^{\prime}\in\{t_{j+1},\dots,t_{n};t_{\text{max}}\}$

[TABLE]

We therefore have

[TABLE]

with

[TABLE]

and $B_{\uparrow}$ is the $\left[(\hat{g}_{\uparrow})_{\alpha_{i}\alpha^{\prime}_{i}}(t_{i},t_{i^{\prime}})\right]_{j+1\leq i,i^{\prime}\leq n}$ matrix where a last line and column corresponding to $t_{\text{max}}$ are added, similar to Eq. (16). However, $\sum_{\alpha_{1}\dots\alpha_{j}}\alpha_{1}\dots\alpha_{j}\prod_{\sigma}A_{\sigma}$ is a contribution to $Z$ at order $j$, and it vanishes according to (12). Therefore $\sum_{\alpha_{1}\dots\alpha_{n}}\alpha_{1}\dots\alpha_{n}\prod_{\sigma}\det D^{\pm}_{\sigma}\simeq 0$, and this proves the clusterization of times in around $t_{\text{max}}$ in the computation of the density.

Appendix B Benchmark

The table below benchmarks the contributions to the density between the $\pm$, LO, and mixed algorithms. We take $t=1$ as our energy unit, and parameters are $\beta t=100$, $\gamma^{2}=0.04t^{2}$, $\epsilon_{d}=-0.36t$, $U=1.2t$, $\alpha=0.3$. Computation effort is 240 CPU*hours for each perturbation order.

[TABLE]

Appendix C Origin of error bar

We have seen in Sec. IV.4 that the contributions to the density have to be normalized by a factor $\eta$, see Eq. (60). To verify that the error bars on the density are not due to this normalization factor, we plot its relative error bars on Figure 4. Blue dots denote the $\pm$ algorithm, orange stars the LO algorithm, and green dots the mixed algorithm. Comparing it to Figure 1, we see that the relative error bars on $\eta$ are much smaller than the ones on the density.