A Quantum Monte Carlo algorithm for out-of-equilibrium Green's functions at long times

TL;DR

This paper introduces a quantum Monte Carlo algorithm capable of computing out-of-equilibrium Green's functions at long times, including the steady state, by efficiently handling perturbative expansions up to high order.

Contribution

It extends a previous Monte Carlo method to reach the infinite time limit and can be applied to both equilibrium and out-of-equilibrium systems, with automatic diagram cancellations.

Findings

Successfully computed Green's functions up to order U^{10} for the Anderson impurity model.

Achieved stable long-time and steady-state calculations due to diagram cancellations.

Validated results through benchmarking with NRG and analytical solutions.

Abstract

We present a quantum Monte-Carlo algorithm for computing the perturbative expansion in power of the coupling constant $U$ of the out-of-equilibrium Green's functions of interacting Hamiltonians of fermions. The algorithm extends the one presented in Phys. Rev. B 91 245154 (2015), and inherits its main property: it can reach the infinite time (steady state) limit since the computational cost to compute order $U^n$ is uniform versus time; the computing time increases as $2^n$. The algorithm is based on the Schwinger-Keldysh formalism and can be used for both equilibrium and out-of-equilibrium calculations. It is stable at both small and long real times including in the stationary regime, because of its automatic cancellation of the disconnected Feynman diagrams. We apply this technique to the Anderson quantum impurity model in the quantum dot geometry to obtain the Green's function and self-energy expansion up to order $U^{10}$ at very low temperature. We benchmark our results at weak and intermediate coupling with high precision Numerical Renormalization Group (NRG) computations as well as analytical results.

Tables

Table 1. Table 1: First coefficients s n , m subscript 𝑠 𝑛 𝑚 s_{n,m} (real) of the self-energy Taylor series Σ ​ ( U , ω ) / Γ = ∑ n , m i m + 1 ​ s n , m ​ ( U / Γ ) n ​ ( ω / Γ ) m Σ 𝑈 𝜔 Γ subscript 𝑛 𝑚 superscript 𝑖 𝑚 1 subscript 𝑠 𝑛 𝑚 superscript 𝑈 Γ 𝑛 superscript 𝜔 Γ 𝑚 \Sigma(U,\omega)/\Gamma=\sum_{n,m}i^{m+1}s_{n,m}(U/\Gamma)^{n}(\omega/\Gamma)^{m} on the equilibrium symmetric model ϵ d = 0 subscript italic-ϵ 𝑑 0 \epsilon_{d}=0 . Coefficients in powers of ω 𝜔 \omega have been obtained by fitting the bare data by a polynomial. We find a good agreement with analytical calculations from Ref.  Yamada, 1975 (2nd and 5th columns), as well as with Bethe ansatz exact calculations from Ref.  Horvatić and Zlatić, 1985 (3rd, 6th and 8th columns).

$m$	$n=2$ (QMC)	$n=2$ (Yamada)	$n=2$ (Bethe)	$n=4$ (QMC)	$n=4$ (Yamada)	$n=4$ (Bethe)	$n=6$ (QMC)	$n=6$ (Bethe)
0	$0\pm 1\times {10}^{-5}$	0		$0\pm 2\times {10}^{-6}$	0		$0\pm 1\times {10}^{-7}$
1	5.39(4)$\times {10}^{-2}$	5.3964$\times {10}^{-2}$	5.3964$\times {10}^{-2}$	5.7(0)$\times {10}^{-4}$	5.6771$\times {10}^{-4}$	5.6482$\times {10}^{-4}$	2.(1)$\times {10}^{-6}$	2.5119$\times {10}^{-6}$
2	5.03(6)$\times {10}^{-2}$	5.0660$\times {10}^{-2}$		1.9(9)$\times {10}^{-3}$	2.0079$\times {10}^{-3}$		3.(1)$\times {10}^{-5}$
3	3.67(5)$\times {10}^{-2}$			4.3(4)$\times {10}^{-3}$			1.(5)$\times {10}^{-4}$
4	2.17(2)$\times {10}^{-2}$			6.2(4)$\times {10}^{-3}$			4.7(1)$\times {10}^{-4}$

Equations

$$\hat{\text{H}}(t)$$

$${\hat{\text{H}}}_0(t)$$

$${\hat{\text{H}}}_{\mathrm{int}}(t)$$

$$G^{aa^{\prime}}_{xx^{\prime}}(t,t^{\prime})=\left(\begin{array}[]{cc}G^{T}_{xx^{\prime}}(t,t^{\prime})&G^{<}_{xx^{\prime}}(t,t^{\prime})\\ G^{>}_{xx^{\prime}}(t,t^{\prime})&G^{\bar{T}}_{xx^{\prime}}(t,t^{\prime})\end{array}\right)$$

$$G_{x{x}^{\prime }}^{a{a}^{\prime }}(t,t^{\prime })\equiv -i⟨T_{\mathrm{c}}\hat{\text{c}}(x,t,a){\hat{\text{c}}}^{†}(x^{\prime },t^{\prime },a^{\prime })⟩$$

$$G(X,X^{\prime })\equiv G\left[(x,t,a),(x^{\prime },t^{\prime },a^{\prime })\right]\equiv G_{x{x}^{\prime }}^{a{a}^{\prime }}(t,t^{\prime })$$

$${\delta }_c\left[((x,t,a),(x^{\prime },t^{\prime },a^{\prime })\right]\equiv \delta (t-t^{\prime }){\delta }_{a{a}^{\prime }}{\delta }_{x{x}^{\prime }}.$$

$$G_{xx^{\prime}}^{aa^{\prime}}(t,t^{\prime})=\sum_{n\geq 0}\frac{i^{n}U^{n}}{n!}\int\prod_{k=1}^{n}du_{k}\sum_{\{x_{k},y_{k}\}}\sum_{\{a_{k}\}}\times\\ \left(\prod_{k=1}^{n}(-1)^{a_{k}}V_{x_{k}y_{k}}(u_{k})\right)\left\llbracket\begin{matrix}(x,t,a),U_{1},\ldots,U_{2n}\\ (x^{\prime},t^{\prime},a^{\prime}),U_{1},\ldots,U_{2n}\end{matrix}\right\rrbracket$$

$$U_{2k-1}$$

$$U_{2k}$$

$$\left\llbracket\begin{matrix}A_{1},\ldots,A_{m}\\ B_{1},\ldots,B_{m}\end{matrix}\right\rrbracket_{ij}\equiv\begin{pmatrix}{g}(A_{1},B_{1})&\ldots&{g}(A_{1},B_{m})\\ \vdots&\ddots&\vdots\\ {g}(A_{m},B_{1})&\ldots&{g}(A_{m},B_{m})\end{pmatrix}_{ij}$$

$$\left\llbracket\begin{matrix}A_{1},\ldots,A_{m}\\ B_{1},\ldots,B_{m}\end{matrix}\right\rrbracket\equiv\det\left\llbracket\begin{matrix}A_{1},\ldots,A_{m}\\ B_{1},\ldots,B_{m}\end{matrix}\right\rrbracket_{ij}$$

$$g(U_{2k},U_{2k})$$

$$g(U_{2k-1},U_{2k-1})$$

$$g(U_{2k-1},U_{2k})$$

$$g(U_{2k},U_{2k-1})$$

$$g(U_{2k},U_{2k})$$

$$g(U_{2k-1},U_{2k-1})$$

$$\left\llbracket\begin{matrix}X,U_{1},\ldots,U_{2n}\\ X^{\prime},U_{1},\ldots,U_{2n}\end{matrix}\right\rrbracket\rightarrow\left\llbracket\begin{matrix}X,Y,U_{1},\ldots,U_{2n}\\ X^{\prime},Y^{\prime},U_{1},\ldots,U_{2n}\end{matrix}\right\rrbracket$$

$$\sum_{\{a_{k}\}}(-1)^{\sum_{k=1}^{n}a_{k}}\left\llbracket\begin{matrix}U_{1},\ldots,U_{2n}\\ U_{1},\ldots,U_{2n}\end{matrix}\right\rrbracket=0$$

$$\left\llbracket\begin{matrix}A_{1},\ldots,A_{m}\\ B_{1},\ldots,B_{m}\end{matrix}\right\rrbracket=\sum_{k=1}^{m}(-1)^{k+1}{g}(A_{k},B_{1})\left\llbracket\begin{matrix}A_{1},\ldots\cancel{A_{k}}\ldots,A_{m}\\ \cancel{B_{1}},B_{2},\ldots,B_{m}\end{matrix}\right\rrbracket$$

$$\left\llbracket\begin{matrix}A_{1},\ldots,A_{m}\\ B_{1},\ldots,B_{m}\end{matrix}\right\rrbracket=\sum_{k=1}^{m}(-1)^{k+1}{g}(A_{1},B_{k})\left\llbracket\begin{matrix}\cancel{A_{1}},\ldots\ldots,A_{m}\\ B_{1},\cancel{B_{k}},\ldots,B_{m}\end{matrix}\right\rrbracket$$

$$\displaystyle(-1)^{i+j}\left\llbracket\begin{matrix}A_{1}\ldots\cancel{A_{i}},\ldots\ldots,A_{m}\\ B_{1},\ldots\cancel{B_{j}},\ldots,B_{m}\end{matrix}\right\rrbracket=$$

$$\displaystyle\left\llbracket\begin{matrix}A_{1},\ldots,A_{m}\\ B_{1},\ldots,B_{m}\end{matrix}\right\rrbracket^{-1}_{ji}\left\llbracket\begin{matrix}A_{1},\ldots,A_{m}\\ B_{1},\ldots,B_{m}\end{matrix}\right\rrbracket$$

$$G_{xx^{\prime}}^{aa^{\prime}}(t,t^{\prime})=g_{xx^{\prime}}^{aa^{\prime}}(t,t^{\prime})+\sum_{n\geq 1}\frac{i^{n}U^{n}}{n!}\int\prod_{k=1}^{n}du_{k}\sum_{\{x_{k},y_{k}\}}\sum_{\{a_{k}\}}\left(\prod_{k=1}^{n}(-1)^{a_{k}}V_{x_{k}y_{k}}(u_{k})\right)\times\\ \left(\sum_{p=1}^{2n}(-1)^{p}g\left[(x,t,a),U_{p}\right]\left\llbracket\begin{matrix}U_{1},\ldots,U_{2n}\\ (x^{\prime},t^{\prime},a^{\prime}),U_{1},\ldots,\cancel{U_{p}},\ldots,U_{2n}\end{matrix}\right\rrbracket+g\left[(x,t,a),(x^{\prime},t^{\prime},a^{\prime})\right]\left\llbracket\begin{matrix}U_{1},\ldots,U_{2n}\\ U_{1},\ldots,U_{2n}\end{matrix}\right\rrbracket\right)$$

$$G_{x{x}^{\prime }}^{a{a}^{\prime }}(t,t^{\prime })=g_{x{x}^{\prime }}^{a{a}^{\prime }}(t,t^{\prime })+\int du\sum _{b,y}(-1{)}^b{g}_{xy}^{ab}(t,u)K_{y{x}^{\prime }}^{b{a}^{\prime }}(u,t^{\prime })$$

$$K(Y,X^{\prime})\equiv(-1)^{b}\sum_{n\geq 1}\frac{i^{n}U^{n}}{n!}\int\prod_{k=1}^{n}du_{k}\sum_{\{x_{k},y_{k}\}}\sum_{\{a_{k}\}}\left(\prod_{k=1}^{n}(-1)^{a_{k}}V_{x_{k}y_{k}}(u_{k})\right)\sum_{p=1}^{2n}(-1)^{p}\delta_{c}(Y,U_{p})\left\llbracket\begin{matrix}U_{1},\ldots,U_{2n}\\ X^{\prime},U_{1},\ldots,\cancel{U_{p}},\ldots,U_{2n}\end{matrix}\right\rrbracket$$

$$G_{x{x}^{\prime }}^{a{a}^{\prime }}(t,t^{\prime })=g_{x{x}^{\prime }}^{a{a}^{\prime }}(t,t^{\prime })+\int du\sum _{b,y}(-1{)}^b{\bar{K}}_{xy}^{ab}(t,u)g_{y{x}^{\prime }}^{b{a}^{\prime }}(u,t^{\prime })$$

$$\bar{K}(X,Y)\equiv(-1)^{b}\sum_{n\geq 1}\frac{i^{n}U^{n}}{n!}\int\prod_{k=1}^{n}du_{k}\sum_{\{x_{k},y_{k}\}}\sum_{\{a_{k}\}}\left(\prod_{k=1}^{n}(-1)^{a_{k}}V_{x_{k}y_{k}}(u_{k})\right)\sum_{p=1}^{2n}(-1)^{p}\delta_{c}(Y,U_{p})\left\llbracket\begin{matrix}X,U_{1},\ldots,\cancel{U_{p}},\ldots,U_{2n}\\ U_{1},\ldots,U_{2n}\end{matrix}\right\rrbracket$$

$$F_{x{x}^{\prime }z}^{a{a}^{\prime }}(t,t^{\prime })\equiv (-i{)}^2⟨T_{\mathrm{c}}\hat{\text{c}}(x,t,a){\hat{\text{c}}}^{†}(x^{\prime },t^{\prime },a^{\prime })\left[{\hat{\text{c}}}^{†}(z,t^{\prime },a^{\prime })\hat{\text{c}}(z,t^{\prime },a^{\prime })-{\alpha }_z\right]⟩$$

$$F_{xx^{\prime}z}^{aa^{\prime}}(t,t^{\prime})=-\sum_{n\geq 0}\frac{i^{n}U^{n}}{n!}\int\prod_{k=1}^{n}du_{k}\sum_{\{x_{k},y_{k}\}}\sum_{\{a_{k}\}}\left(\prod_{k=1}^{n}(-1)^{a_{k}}V_{x_{k}y_{k}}(u_{k})\right)\left\llbracket\begin{matrix}(x,t,a),(z,t^{\prime},a^{\prime}),U_{1},\ldots,U_{2n}\\ (x^{\prime},t^{\prime},a^{\prime}),(z,t^{\prime},a^{\prime}),U_{1},\ldots,U_{2n}\end{matrix}\right\rrbracket$$

$$F_{xx^{\prime}z}^{aa^{\prime}}(t,t^{\prime})=-\sum_{n\geq 0}\frac{i^{n}U^{n}}{n!}\int\prod_{k=1}^{n}du_{k}\sum_{\{x_{k},y_{k}\}}\sum_{\{a_{k}\}}\left(\prod_{k=1}^{n}(-1)^{a_{k}}V_{x_{k}y_{k}}(u_{k})\right)\times\\ \left(\sum_{p=1}^{2n}(-1)^{p+1}g\Bigl{(}(x,t,a),U_{p}\Bigr{)}\left\llbracket\begin{matrix}(z,t^{\prime},a^{\prime}),U_{1},\ldots,U_{2n}\\ (x^{\prime},t^{\prime},a^{\prime}),(z,t^{\prime},a^{\prime}),U_{1},\ldots,\cancel{U_{p}},\ldots,U_{2n}\end{matrix}\right\rrbracket+\right.\\ \left.g\Bigl{(}(x,t,a),(x^{\prime},t^{\prime},a^{\prime})\Bigr{)}\left\llbracket\begin{matrix}(z,t^{\prime},a^{\prime}),U_{1},\ldots,U_{2n}\\ (z,t^{\prime},a^{\prime}),U_{1},\ldots,U_{2n}\end{matrix}\right\rrbracket-g\Bigl{(}(x,t,a),(z,t^{\prime},a^{\prime})\Bigr{)}\left\llbracket\begin{matrix}(z,t^{\prime},a^{\prime}),U_{1},\ldots,U_{2n}\\ (x^{\prime},t^{\prime},a^{\prime}),U_{1},\ldots,U_{2n}\end{matrix}\right\rrbracket\right)$$

$$F_{x{x}^{\prime }z}^{a{a}^{\prime }}(t,t^{\prime })=-g_{x{x}^{\prime }}^{a{a}^{\prime }}(t,t^{\prime })[G_{zz}^{<}(t^{\prime },t^{\prime })-i{\alpha }_z]+g_{xz}^{a{a}^{\prime }}(t,t^{\prime })G_{z{x}^{\prime }}^{<}(t^{\prime },t^{\prime })-\int du\sum _{b,y}(-1{)}^b{g}_{xy}^{ab}(t,u)L_{y{x}^{\prime }z}^{b{a}^{\prime }}(u,t^{\prime })$$

$$L_{yx^{\prime}z}^{ba^{\prime}}(u,t^{\prime})\equiv(-1)^{b}\sum_{n\geq 1}\frac{i^{n}U^{n}}{n!}\int\prod_{k=1}^{n}du_{k}\sum_{\{x_{k},y_{k}\}}\sum_{\{a_{k}\}}\left(\prod_{k=1}^{n}(-1)^{a_{k}}V_{x_{k}y_{k}}(u_{k})\right)\times\\ \sum_{p=1}^{2n}(-1)^{p+1}\delta_{c}\Bigl{(}(y,u,b),U_{p}\Bigr{)}\left\llbracket\begin{matrix}(z,t^{\prime},a^{\prime}),U_{1},\ldots,U_{2n}\\ (x^{\prime},t^{\prime},a^{\prime}),(z,t^{\prime},a^{\prime}),U_{1},\ldots,\cancel{U_{p}},\ldots,U_{2n}\end{matrix}\right\rrbracket$$

Full text

I Introduction

The study of the out-of-equilibrium regime of strongly correlated many-body quantum problems is a subject of growing interest in theoretical condensed matter physics, in particular due to a rapid progress in experiments with e.g. the ability to control light-matter interaction on ultra-fast time scaleFörst et al. (2011), light-induced superconductivity Fausti et al. (2011); Nicoletti et al. (2014); Casandruc et al. (2015); Nicoletti and Cavalleri (2016); Nicoletti et al. (2018) or metal-insulator transition driven by electric field Nakamura et al. (2013). The development of high precision and controlled computational methods for non-equilibrium models in strongly correlated regimes is therefore very important. Even within an approximated framework such as Dynamical Mean Field Theory Georges et al. (1996); Kotliar et al. (2006); Aoki et al. (2014) (DMFT), which reduces bulk lattice problem to the solution of a self-consistent quantum impurity model, efficient numerically exact real time out-of-equilibrium quantum impurity solver algorithms are still lacking.

The long time steady state of non-equilibrium strongly interacting quantum systems is specially difficult to access within high precision numerical methods. Until recently, most approaches were severely limited in reaching long times, e.g. the density matrix renormalization group (DMRG)White (1992, 1993); Schollwöck (2005) faces entanglement growth at long times. Early attempts of real time quantum Monte Carlo Mühlbacher and Rabani (2008); Werner et al. (2009, 2010); Schiró and Fabrizio (2009); Schiró (2010) also experienced an exponential sign problem at long time and large interaction. Within Monte-Carlo methods, two main routes are currently explored to resolve this issue: the inchworm algorithm Cohen et al. (2014a, b, 2015); Chen et al. (2017a, b) and the so-called “diagrammatic” Quantum Monte Carlo Profumo et al. (2015) (QMC). Diagrammatic QMCProkof’ev and Svistunov (1998); Mishchenko et al. (2000); Van Houcke et al. (2010); Prokof’ev and Svistunov (2007, 2008); Gull et al. (2010); Kozik et al. (2010); Pollet (2012); Van Houcke et al. (2012); Kulagin et al. (2013a, b); Gukelberger et al. (2014); Deng et al. (2015); Huang et al. (2016); Rossi et al. (2018); Van Houcke et al. (2019) computes the perturbative expansion of physical quantities in power of the interaction $U$, using an importance sampling Monte Carlo. In Ref. Profumo et al., 2015, some of us have shown that, when properly generalized to the Schwinger-Keldysh formalism, this approach yields the perturbative expansion in the steady state, i.e. at infinite time. We showed that, by regrouping the Feynman diagrams into determinants and summing explicitly on the Keldysh indices of the times of the vertices of the expansion, we eliminate the vacuum diagrams and obtain a clusterisation property allowing us to take the long time limit. The sum over the Keldysh indices implies a minimal cost of $O(2^{n})$ to compute the order $n$, but uniformly in time, at any temperature. We refer to this class of algorithms as “diagrammatic” for historical reasons, as their first versions (in imaginary time) were using an explicit Markov chain in the space of Feynman diagrams. However, modern versions of the algorithms regroup diagrams with only an exponential number of determinants (instead of sampling the $n!$ diagrams), eliminating the vacuum diagrams, both in real timeProfumo et al. (2015) and in imaginary timeRossi (2017); Moutenet et al. (2018), which leads to much higher performance.

In this article, we generalize the algorithm presented in Ref. Profumo et al., 2015 to the calculations of Green’s functions. Indeed, in its initial formulation it only permits the calculation of physical observables at equal time such as the density or current, and the full Green’s function requires the computation of each time one by one, which is not technically viable. Here, we show how to use a kernel technique in order to obtain efficiently the perturbative expansion of the Green’s function and the self-energy, as a function of time or frequency. Computing the Green’s function is an important extension of the technique. First, it is the first step towards building a DMFT real time non-equilibrium impurity solver. Second, even in the simple context of a quantum dot, each computation provides much more information than the original algorithm (from which only a single number, e.g. the current, could be computed).

The central issue of the “diagrammatic” QMC family is to properly sum the perturbative expansion away from the weak coupling regime, specially given the fact that one has access to a limited number of orders (about $10$ in the present case) due to the exponential cost with the order $n$. Some of us will address this issue in a separate paperBertrand et al. (2019), using the building blocks introduced here. In this paper, we present the formalism and first benchmark our approach in the weak coupling regime.

This paper is organized as follows: in section II, we introduce the necessary formalism and derive the basic equations that will be used to formulate the method. Section III discusses our sampling strategy for the QMC algorithm, as well as relations with previous work. Section IV shows our numerical data and the detailed comparison with our benchmarks in the weak coupling regime.

II Wick determinant formalism

This section is devoted to the derivation of the main formula needed to set up the QMC technique. We introduce a systematic formalism that uses what we call “Wick determinants”. Although the formalism is nothing but the usual diagrammatic expansion (in Keldysh space), its Wick determinant formulation provides a route for deriving standard results (such as equation of motions) in a self contained manner that does not require to introduce Feynman diagrams. We find this approach useful for discussing numerical algorithms.

II.1 Notations and main expansion formula

We consider a generic class of system given by a time-dependent Hamiltonian of the form,

[TABLE]

where $\hat{\mathrm{{\textbf{H}}}}_{0}(t)$ is a quadratic unperturbed Hamiltonian. The operators $\hat{\mathrm{{\textbf{c}}}}^{\dagger}_{x}$ ($\hat{\mathrm{{\textbf{c}}}}_{x}$) are the usual creation (destruction) operators on site $x$. $x$ and $y$ index all discrete degrees of freedom such as sites, orbitals, spin and/or electron-hole (Nambu) degrees of freedom and will simply be called orbital indices. $\hat{\mathrm{{\textbf{H}}}}_{\rm int}(t)$ is the, possibly time-dependent, electron-electron interaction perturbation which is assumed to be switched on at $t=0$. Without loss of generality we assume $V_{xy}=V_{yx}$. We emphasize that the method described in this paper is not restricted to this form of interaction (as shown in Ref. Profumo et al., 2015) and can be generalized straightforwardly to arbitrary interactions. However, to improve readability, we will restrict hereafter our presentation to the case of density-density interactions. We also add the quadratic shift $\alpha$, which has been introduced in previous works Rubtsov and Lichtenstein (2004); Profumo et al. (2015); Wu et al. (2017). In particular, we have shown in Ref. Profumo et al., 2015 that, in the context of real time QMC, it can strongly affect the radius of convergence of the perturbative series. The non-interacting Hamiltonian is assumed to be already solved, i.e. one has calculated all the one-particle non-interacting Green’s functions. Such calculations can be done even out-of-equilibrium using e.g. the formalism discussed in Ref. Gaury et al., 2014.

Our starting point is a formula for the systematic expansion of interacting Green’s function in powers of the parameter $U$. We use the Schwinger-Keldysh formalism to produce this expansion. The Green’s functions acquire additional Keldysh indices $a,a^{\prime}\in\{0,1\}$ which provides the Green’s function with a $2\times 2$ structure,

[TABLE]

where $G^{T}_{xx^{\prime}}(t,t^{\prime})$, $G^{<}_{xx^{\prime}}(t,t^{\prime})$, $G^{>}_{xx^{\prime}}(t,t^{\prime})$ and $G^{\bar{T}}_{xx^{\prime}}(t,t^{\prime})$ are respectively the standard time ordered, lesser, greater and anti-time ordered Green’s functions. We use a similar definition for the (known) non-interacting Green’s function $g^{aa^{\prime}}_{xx^{\prime}}(t,t^{\prime})$. We introduce the Keldysh points $X$ that gather an orbital index $x$, a time $t$ and a Keldysh index $a$ to form the tuple $X\equiv(x,t,a)$. Using Keldysh points, we can rewrite the above definitions using the standard conventions of Keldysh formalism,

[TABLE]

where the creation ($\hat{\mathrm{{\textbf{c}}}}^{\dagger}(X^{\prime})$) and annihilation operators ($\hat{\mathrm{{\textbf{c}}}}(X)$ or $\hat{\mathrm{{\textbf{c}}}}(x,t,a)$), here in Heisenberg representation, are ordered using the contour time ordering operator $T_{\rm c}$. $T_{\rm c}$ orders first by Keldysh index ($a$) before ordering by increasing time within the forward branch ($a=0$), and by decreasing time within the backward branch ($a=1$) eventually multiplying the result by the usual fermionic $(-1)$ factor whenever an odd number of permutations have been performed. In several places, we will use the alternative notation for the Green’s function

[TABLE]

and we will also note the $\delta$ function on the Keldysh contour as

[TABLE]

Using the above notations (with $\hbar=1$), one can derive the usual expansion in power of $U$ using Wick’s theorem. We first assume $\alpha_{x}=0$ at all orbital indices $x$. We will explain at the end of this paragraph how to extend this formula to the general case $\alpha_{x}\neq 0$. We obtainProfumo et al. (2015):

[TABLE]

which forms the starting point of this work. Here, we have noted for $1\leq k\leq n$

[TABLE]

and introduced the Wick matrix:

[TABLE]

where $A_{i}$ and $B_{j}$ are two sets of $m$ points on the Keldysh contour. We refer to the determinant of the Wick matrix as the Wick determinant. For notation simplicity, the determinant is assumed in the absence of indices,

[TABLE]

In equation (8), we start at $t<0$ with a non-interacting state and switch on the interaction for $t\geq 0$. Hence the time integrals in Eq. 8 run over the segment $[0,t_{M}]$, where $t_{M}=\max(t,t^{\prime})$. The lower boundary simply arises from $V_{xy}(u<0)=0$. The upper boundary can be extended to any value larger than $t_{M}$ without changing the integral’s value (standard property of the Keldysh formalism). For the practical applications shown in section III, we will a fixed large value of $t_{M}$. We emphasize that the complexity of the algorithm does not grow with $t_{M}$. Eq. (8) is formally very appealing: it reduces the problem of calculating the contributions of the expansion to a combination of linear algebra and quadrature.

The definition Eq. (10) contains an ambiguity that needs to be clarified: the ordering at equal times of terms like $g(U_{2k},U_{2k})$ is ill defined. For these terms, one must keep the same ordering of the creation and destruction operators as in the original definition of the interacting Hamiltonian Eq. (1), i.e.

[TABLE]

To proceed to the general case $\alpha_{x}\neq 0$, one only needs to shift the diagonal terms of the Wick matrix using the following replacement rules, as shown in Ref.Profumo et al., 2015, Appendix A:

[TABLE]

As a result, all derivations in this paper can be done by first ignoring $\alpha_{x}$, then replacing the non-interacting Green’s functions with these rules. For this reason and for readability, we will keep these replacements implicit in Wick matrices, but explicit otherwise.

Finally, Eq. (8) can be extended to the calculations of arbitrary Green’s functions. The rule for doing so is as follows: whenever one introduces a product $-i\hat{\mathrm{{\textbf{c}}}}(Y)\hat{\mathrm{{\textbf{c}}}}^{\dagger}(Y^{\prime})$ under the time ordering operator in Eq. (5), one must add the corresponding Keldysh points in the Wick determinant of Eq. (8):

[TABLE]

If $Y$ and $Y^{\prime}$ share the same time and orbital index $y$, we have the possibility to introduce terms of the form $-i[\hat{\mathrm{{\textbf{c}}}}(Y)\hat{\mathrm{{\textbf{c}}}}^{\dagger}(Y^{\prime})+\alpha_{y}]$ in the definition of the Green’s function. In that case, one must replace $g(Y,Y^{\prime})\rightarrow g(Y,Y^{\prime})-i\alpha_{y}$ in the diagonal of the Wick matrix. Again, to improve readability, we will keep this replacement implicit in Wick matrices, but explicit otherwise.

II.2 A few properties of Wick determinants

Wick determinants have the general properties of determinants: exchanging two Keldysh points on either the first or the second line of the left hand side of Eq. (10) leaves the Wick determinant unchanged up to a factor $(-1)$. An important property of the formalism, as already shown in Ref. Profumo et al., 2015, is that for $n>0$:

[TABLE]

for any times $u_{1},\ldots,u_{n}$ and orbital indices $x_{1},\ldots,x_{n}$ and $y_{1},\ldots,y_{n}$. This relation expresses the fact that vacuum diagrams are automatically cancelled by the sum over the Keldysh indices, even before any integration over time. It is proven in the usual way in the Keldysh formalism, by considering the largest $u_{k}$ time, say $k=p$. From the properties of the bare Green’s functions, one can show that the elements of the Wick matrix, hence the determinant, are in fact all independent of $a_{p}$ (reflecting the fact that the largest time can be on the upper or the lower part of the contour). Therefore the sum over $a_{p}$ cancels the sum.

We will use the usual expansion of a determinant along one row or one column in terms of the cofactor matrix. It takes the form

[TABLE]

for the expansion along the first column and

[TABLE]

for the expansion along the first row. The notation $\cancel{A_{k}}$ ($\cancel{B_{k}}$) stands for the fact that the corresponding row or column must be removed from the Wick matrix.

Last, we will also make a systematic use of the fact that the cofactor matrix is directly related to the inverse of the matrix,

[TABLE]

II.3 Definition of the kernel $K$ for the one-body Green’s function

In Ref. Profumo et al., 2015, a QMC scheme was defined directly on Eq. (8) so that a single QMC run could provide the value of $G_{xx^{\prime}}^{aa^{\prime}}(t,t^{\prime})$ for a single pair of times $t$ and $t^{\prime}$. In order to extend the technique and obtain a full curve (as a function of $t$) in a single run, a different form must be used. Performing the expansion of Eq. (17) on Eq. (8), we obtain

[TABLE]

The last term of the sum vanishes for $n>0$ due to Eq. (15). Factorizing the $g$ from the sum, we arrive at

[TABLE]

where the kernel $K_{yx^{\prime}}^{ba^{\prime}}(u,t^{\prime})=K(Y,X^{\prime})$ with $Y=(y,u,b)$ is defined by

[TABLE]

Equations (20) and (21) will be the basis of one of the method developed in this article. Eq. (21) will provide the mean to get a full $t$-curve in a single calculation and Eq. (20) to relate the corresponding kernel to the Green’s function $G$, the target of the calculation.

A symmetric kernel $\bar{K}$ may be derived following the exact same route but now expanding the Wick determinant along the first column using Eq. (16). We find

[TABLE]

where the kernel $\bar{K}$ is defined by

[TABLE]

II.4 Definition of the kernel $L$ of the $F$ Green’s function

Let us define a new Green’s function with 4 operators, the $F$ function. As we shall see, the $F$ Green’s function can also be represented in term of a kernel so that we will be able to design a direct QMC method to calculate it. Its interest stems from the fact that it can be used to reconstruct $G$ while the corresponding QMC technique will be more precise at high frequency. It is defined as,

[TABLE]

In the next paragraph, we shall prove that $F$ is essentially equal to $\bar{K}$ (up to interacting matrix elements). The function $F$ is known to provide a better estimate of the Green’s function. It has been used in the context of the Numerical Renormalization Group (NRG)Bulla et al. (1998) as well as in imaginary time QMC methods as an improved estimator Hafermann et al. (2012).

The expansion of $F$ reads,

[TABLE]

To obtain the kernel of $F$, we expand the determinant along its first row using Eq. (17),

[TABLE]

Identifying the two last terms with the corresponding expansion of $G$, we arrive at,

[TABLE]

where the kernel $L$ is defined by

[TABLE]

II.5 Relation between $F$, $\bar{K}$ and $G$: equations of motion

Here, we show that the expressions for the different kernels can be formally integrated to provide connections between the different kernels and Green’s functions. We will arrive at expressions that are essentially what can be obtained directly using equations of motions.

We start with the expression of $\bar{K}$, Eq. (23). The first step is to realize that the sum over the $2n$ determinants provide identical contributions to the kernel. Indeed, odd $p=2k-1$ and even $p=2k$ values of $p$ provide identical contributions due to the symmetry of $V_{xy}$. Similarly, odd values $p=2k-1$ have the same contribution as $p=1$ as can be shown by using the symmetry properties of the determinant and exchanging the role of $U_{1}\leftrightarrow U_{2k-1}$ and $U_{2}\leftrightarrow U_{2k}$ in the sums and integration. We arrive at,

[TABLE]

We can now perform explicitly the integral over $u_{1}$ where the delta function yields, for $u\in[0,t_{M}]$ ($\bar{K}$ is zero otherwise):

[TABLE]

This shows the kernel $\bar{K}$ is no more than a sum of 2-particle Green’s functions. The relation between $\bar{K}$ and $G$ in Eq. (22) can then be transformed into:

[TABLE]

which can be used to reconstruct $G$ from the knowledge of $F$. This relation is the well known equation of motion for $G$. It also shows that $F$ is essentially the convolution of $G$ with the self-energy.

As a side note, we show in Appendix B that the kernel $L$ can also be expressed in terms of Green’s functions by following the same formalism, in accordance with the equation of motion for $F$.

II.6 Retarded and advanced kernels

As the retarded (or advanced) Green’s functions directly give the spectral functions, they are of particular interest. At equilibrium, they contain all information which can be obtained from the Keldysh Green’s function. We show here simple relations to compute them from the kernels $K$, $\bar{K}$ or $L$.

The retarded and advanced Green’s functions relate to the lesser and greater Green’s functions as follows:

[TABLE]

where $\theta$ is the Heaviside function. From the definitions of the time-ordered and time-antiordered Green’s functions, these can also be written as:

[TABLE]

where $a$ can be any Keldysh index. These are also valid for the non-interacting $g$. As $\bar{K}$ is a sum of Green’s function, one may define in the same way a retarded version of $\bar{K}$, denoted $\bar{K}^{R}$. We can see from Eq. (32) and the definition of $F$ in Eq. (24) that $\bar{K}^{R}$ follows the same properties:

[TABLE]

We now show a simple relation between $G^{R}$, $g^{R}$ and $\bar{K}^{R}$. Plugging Eq. (22) into Eq. (36), one gets:

[TABLE]

This simplifies into:

[TABLE]

Similar relations may be derived with $K$, $F$ and $L$. In fact, for any function from $K$, $F$ and $L$, which all depends on two times and two Keldysh indices, we choose to define a retarded and advanced function in the same way as Eq. (34) and Eq. (35). As all of them are sums of Green’s functions, one may show that they all verify similar properties as in Eq. (36) and (37). Then from Eq. (20) follows:

[TABLE]

and from Eq. (27) follows:

[TABLE]

III Quantum Monte Carlo technique

We now turn to the stochastic algorithms that will be used for the actual evaluations of the multi-dimensional integrals that define the expansion of the Green’s function. These algorithms are direct extensions of the algorithm of Ref. Profumo et al., 2015 and inherit of most of its properties. The main novelty lies in using kernels which permits the calculation of the full time dependency of the Green’s function in a single QMC run.

III.1 Sampling of the kernel $K$

Let us first discuss the calculation of $G$ using the kernel $K$, using Eq. (20) and (21). We rewrite by explicitly separating the sum over Keldysh indices (which will be summed explicitly) and the sum over space and integral over time (which will be sampled using Monte-Carlo). This separation was shown to be extremely important in Ref. Profumo et al., 2015. The kernel takes the form,

[TABLE]

We define a configuration $\cal C$ as

• •

The order $n$.

• •

A set of times $\{u_{k}\in[0,t_{M}]\}$ for $1\leq k\leq n$.

• •

Two sets of indices $\{x_{k}\}$ and $\{y_{k}\}$ for $1\leq k\leq n$.

and the sum over all configuration as the integral over the $u_{k}$ and the sum over the $x_{k},y_{k}$. For practical purpose, the time integrals run over a finite interval $[0,t_{M}]$. In accordance with the remark following Eq. (8), $t_{M}$ can be chosen to be any value larger than $t$ and $t^{\prime}$ of the target Green’s functions $G(t,t^{\prime})$.

The kernel takes the form,

[TABLE]

where the sum over the configurations $\cal C$ is a compact notation for the sum and integrals of Eq.(43). We observe that a single configuration provides values of $K$ for $2n$ different points $Y$ through the delta function in the preceding expression. To sample the sum over configurations, we need to define a positive function $W(\cal C)$ that will provide the (unnormalized) probability for the configuration $\cal C$ to be visited by the QMC algorithm. We define this weight as,

[TABLE]

Noting $Z_{\rm QMC}\equiv\sum_{\cal C}W(\cal C)$, the kernel can be rewritten as,

[TABLE]

where the average is taken over the distribution $W({\cal C})/Z_{\rm QMC}$. $Z_{\rm QMC}$ is an effective partition function associated to the QMC algorithm. Note that by construction, the weight $W({\cal C})$ controls the measurement, i.e. $|W^{n}_{p}(X^{\prime},\{U_{k}\},a)|\leq W({\cal C})$. This is an essential property for a reweighting technique since it guarantees that the weight $W^{n}_{p}(X^{\prime},\{U_{k}\},a)/W({\cal C})$ does not diverge (which can produce an infinite variance, for an example of this effect in the context of determinantal Monte-Carlo see e.g. Ref. Shi and Zhang, 2016).

To sample $W({\cal C})$ and evaluate $Z_{\rm QMC}$, we use the continuous time Monte-Carlo technique that was discussed in details in Ref. Profumo et al., 2015. We use moves that change the order $n$ by $\pm 1$ so that all orders (up to a maximum one) can be calculated in a single run. The algorithm has very good ergodicity properties since the order $n=0$ configuration is visited regularly. Each configuration $\cal C$ provides $2n$ values of $Y=(y,u,b)$ which are recorded by binning with the weight of Eq. (47).

The partial weights $W_{p}^{n}$ possess an essential clusterization property, which generalizes the one discovered in Ref. Profumo et al., 2015: if one or several times $u_{k}$ goes to infinity (i.e. is far from $t^{\prime}$), then all $W_{p}^{n}$ goes to 0. In other words the integrand is localized around $t^{\prime}$. A detailed proof is provided in Appendix  A. The point $X^{\prime}=(x^{\prime},t^{\prime},a^{\prime})$ is kept fixed during the calculation to anchor the integral around this point. An important consequence of the clusterization property is that the computational cost of the algorithm is uniform in $t_{M}$. Indeed, as one increases $t_{M}$, one simply adds regions of the configuration space that have a vanishingly small weight, hence do not contribute to the integral and do not get sampled.

Last, in order to calculate the factors $W_{p}^{n}$, we use the fact that they are made of the cofactors of the original matrix, hence can be rewritten as,

[TABLE]

This last form is very convenient since a single Wick matrix (and its inverse) need to be stored and monitored during the calculation.

III.2 Sampling of the kernel $L$

Following the same route for $L$ as was done for $K$, we can write:

[TABLE]

thus defining $W_{p}^{n}$ for $p=3,4,\ldots,2n+2$. We also define $W_{1}^{n}$ and $W_{2}^{n}$ in the following way:

[TABLE]

These two extra values are necessary to compute $G^{<}_{zz}(t^{\prime},t^{\prime})$ and $G^{<}_{zx^{\prime}}(t^{\prime},t^{\prime})$, which are needed to obtain $F$, as can be seen in Eq. (27). Moreover, they do not require extra computation time, as they are a direct by-product of the computation of the $W_{p}^{n}$ for $p>2$. Indeed, in the same spirit as Eq. (48), the determinant within any $W_{p}^{n}$ (for $p=1,\ldots,2n+2$) can be replaced by:

[TABLE]

Again, a single Wick matrix is needed to get contributions to all $W_{p}^{n}$ (given a set of Keldysh indices), which is very convenient.

Configurations are defined in the same way as in the previous section, and the weight of a configuration $\cal C$ is now:

[TABLE]

We define again $Z_{\rm QMC}\equiv\sum_{\cal C}W(\cal C)$ (which however has a different value than in the previous section). Finally, $L$ can be written as:

[TABLE]

and, from Eq. (8), we get that the values of $G^{<}_{zz}(t^{\prime},t^{\prime})$ and $G^{<}_{zx^{\prime}}(t^{\prime},t^{\prime})$ (needed to compute $F$) are:

[TABLE]

The Monte-Carlo algorithm used to evaluate these averages is the same as in the previous section, except for the weight $W({\cal C})$.

III.3 A discussion of the Werner et al. approachWerner et al. (2010)

In this paragraph, we discuss the relation of this work with a preceding work Werner et al. (2009, 2010) that also implements an expansion in powers of $U$ within the Keldysh formalism. Although both results are consistent, Ref. Werner et al., 2010 has two important limitations which are not present in the method presented here. First, it suffers from a very large sign problem that increases drastically with time, while we do not experience a sign problem. Typical data shown in Ref. Werner et al., 2010 corresponds to a maximum time of $t_{M}\approx 5/\Gamma$ between the switching of the interaction and the measurement of the observable while we found that $t_{M}\approx 20/\Gamma$ is needed to enter the stationary result at order $n=8$ as shown in Fig. 1. A direct consequence of this issue is that Ref. Werner et al., 2010 cannot access the small bias regime where the Kondo effect is present: since the Kondo temperature $T_{K}$ is typically much smaller than $\Gamma$, long simulation times $t_{M}\gg 1/T_{K}$ are needed to capture the Kondo physics. Second, the technique of Ref. Werner et al., 2010 suffered from a very large sign problem outside of the electron-hole symmetry point so that only this point could be studied.

An interesting aspect of Ref. Werner et al., 2010 is that some results could be obtained in regimes where the “sign” of the Monte-Carlo calculation was very small $\sim 10^{-3}$ (see for instance Fig.5 of Ref. Werner et al., 2009). Such a small sign is usually associated with very large error bars that prevents practical calculations to be performed. In the rest of this paragraph, we make a simple technical remark that explains the origin of this behaviour.

The main expansion formula used in the present work is Eq. (8) which provides the expansion for the Keldysh Green’s function $G_{xx^{\prime}}^{aa^{\prime}}(t,t^{\prime})$. A similar formulaProfumo et al. (2015) provides the sum of vacuum diagrams, sometimes called the Keldysh “partition function” $Z$,

[TABLE]

We have $Z=1$ in the Keldysh formalism, reflecting the unitarity of quantum mechanics. We see that (15) indeed implies $Z=1$, and that the cancellation of vacuum diagrams is due to the sum over Keldysh indices. The integrand is identically zero.

Let’s note $Z(\{U_{i}\})$ the integrand of Eq. (58) (without the sum over Keldysh indices):

[TABLE]

Ref. Werner et al., 2009 Monte-Carlo samples the absolute value of this integrand $|Z(\{U_{i}\})|$ (the authors actually used auxiliary Ising variables but that does not impact the present argument). We also note $G(X,X^{\prime},\{U_{i}\})$ the integrand of Eq. (8) and $Z_{\rm QMC}$ the integral of $|Z(\{U_{i}\})|$:

[TABLE]

The weight of the QMC is $|Z(\{U_{i}\})|/Z_{\rm QMC}$. We have

[TABLE]

where the average is the Monte-Carlo average. The denominator of Eq. (61) is the QMC sign mentioned above. From $Z=1$, we find that this QMC sign is simply given by $\left\langle Z(\{U_{i}\})/|Z(\{U_{i}\})|\right\rangle=1/Z_{\rm QMC}$.

Now, we note that the probability to visit order 0 is also $|Z(\varnothing)|/Z_{\rm QMC}=1/Z_{\rm QMC}$. Therefore the average sign in the denominator of Eq. (61) is the probability to visit a configuration at zeroth order. This probability decreases when $U$ is increased, or at long time, when higher orders are sampled which explains why the QMC sign was observed to decrease drastically in Ref. Werner et al., 2009. However, this QMC sign is always positive and could a priori be computed very efficiently as an integral of a positive function using e.g. the technique in Ref. Profumo et al., 2015. A genuine sign problem can however results from the numerator of Eq. (61).

IV Application to the Anderson impurity model

We now turn to the illustration of our new techniques with calculations done on the Anderson impurity model. The implementation of our technique was based on the TRIQS packageParcollet et al. (2015). We only present here results that showcase the technique and differ the exploration of the physics of the model, in particular the Kondo physics out of equilibrium, to the companion article of the present work Bertrand et al. (2019). We stress that the QMC technique is not restricted to impurity models and also applies to lattice models such as the Hubbard model.

IV.1 Definition of the model

In the Anderson impurity model, the impurity is described by the operators $\hat{\mathrm{{\textbf{d}}}}^{\dagger}_{\sigma}$ ($\hat{\mathrm{{\textbf{d}}}}_{\sigma}$) that create (destroy) an electron on the impurity with spin $\sigma$. The impurity is connected to two non-interacting electrodes via a tunneling Hamiltonian. Instead of providing explicitly this tunneling Hamiltonian and the Hamiltonian of the electrodes, it is simpler to write directly the non-interacting Green’s function of the impurity. We work with its wide band form which is appropriate for the low energy physics of a regular impurity model. The retarded Green’s function reads in the frequency domain,

[TABLE]

where the parameter $\Gamma$ sets the width of the resonance in the absence of electron-electron interactions and the on-site energy $\epsilon_{d}$ sets the resonance with respect to the Fermi level. Eq. (62) fully defines the model at equilibrium. The model is made non trivial through the interacting terms that reads,

[TABLE]

where $\hat{\mathrm{{\textbf{n}}}}_{\sigma}=\hat{\mathrm{{\textbf{d}}}}^{\dagger}_{\sigma}\hat{\mathrm{{\textbf{d}}}}_{\sigma}$ is the impurity electronic density of spin $\sigma$ and the Heaviside function $\theta(t)$ represents the fact that the interaction is switched on at $t=0$. The calculations are performed up to large times $t_{M}$ so that the system has relaxed to its stationary regime which corresponds to the interacting system at the bath temperature. All calculations are performed at very low temperature $k_{\rm B}T=10^{-4}\Gamma$, although the method is suitable for finite temperature as well.

The main output of our calculations is the expansion for the interacting retarded Green’s function. Restricting ourselves to the stationary limit, it is a function of $t-t^{\prime}$ only and can be studied in the frequency domain.

[TABLE]

from which one can obtain the corresponding quantity in the frequency domain by fast Fourier transform,

[TABLE]

Our technique typically provides the first $N=10$ terms of this expansion as we show next. Last, we define the spectral function (or interacting local density of state)

[TABLE]

and the retarded self energy $\Sigma^{R}(\omega)$,

[TABLE]

IV.2 Numerical Results order by order

We now present the numerical data obtained by sampling the kernel $K$. The left panels of Fig. 1 shows an example of the bare data for the retarded kernel $K^{R}(t)$ as they come out of the calculation for order $U^{2}$, $U^{4}$, $U^{6}$ and $U^{8}$ (top to bottom). Note that the noise in these data is mostly apparent, it corresponds to noise at very high frequency. This apparent noise reflects the fact that we have binned the curve $K^{R}(t)$ into a very fine grid ($50000$ grid points in this calculation). An even finer grid would show even more apparent noise (since there would be even less Monte-Carlo points per grid point). The corresponding cumulative function $\int_{0}^{t}K^{R}(u)du$ is however noiseless as can be seen from the example shown in Fig. 2 for $n=8$.

The next step is to make a Fast Fourier Transform of $K^{R}_{n}(t)$ (not shown). The resulting $K^{R}_{n}(\omega)$ is relatively noisy at high frequency. Last, we obtain $G^{R}_{n}(\omega)=g^{R}(\omega)K^{R}_{n}(\omega)$ for $n\geq 1$ (from Eq. (41)) as shown in the right panel of Fig. 1. The factor $g^{R}(\omega)$, which decays at high frequency, very efficiently suppresses the high frequency noise of the kernel data. The same noise-reduction mechanism has been reported in e.g. Ref. Gull et al., 2008 in the context of auxiliary-field Monte-Carlo. We emphasize one aspect of these data which is rather remarkable: even though the eight order contribution $G^{R}_{8}(\omega)$ is the result of an eight dimensional integral and is more than four orders of magnitude smaller than the second order contribution $G^{R}_{2}(\omega)$, it can be obtained with high precision (the error bars are of the order of the thickness of the lines here). This is due to the recursive way these integrals are calculated as discussed in Ref. Profumo et al., 2015.

Using the definition Eq. (67) of the Self energy, we can obtain a recursive expression for $\Sigma^{R}_{n}$ in term of the Green’s function expansion:

[TABLE]

for $n>1$ with $\Sigma^{R}_{1}(\omega)=[g^{R}(\omega)]^{-2}G^{R}_{1}(\omega)$. The corresponding data is shown in Fig. 3 where we plot the coefficients $\Sigma_{n}^{R}(\omega)$ for $n=2,4,6,8$ and $10$. The error bars increases with the order $n$ which we attribute to the fact that, since the self energy only contains one-particle irreducible diagrams, it is the subject of many cancellations of terms. Indeed, one finds that the decay of $\Sigma^{R}_{n}(\omega)$ with $n$ is rather rapid with seven orders of magnitude between the first and the tenth order.

Our first benchmark uses a reference calculation made by YamadaYamada (1975). The result at order $2$ is compared with the result of Yamada in the left panel of Fig. 3 and found to be in excellent agreement. In his seminal work Yamada also provided analytical calculations at order $2$ and $4$ in the form of a low frequency expansion for the particle-hole symmetric impurity,

[TABLE]

Table 1 shows the results of Yamada ($m=1,2$ and $n=2,4$ ) as well as ours (obtained by fitting our numerical data at low frequency). We find a good quantitative agreement with Yamada results. Yamada also provided numerical results at $n=4$ which are almost featureless and in very poor agreement with our data.

Our second method uses the kernel $L$ in order to calculate the Green’s function expansion. The bare data is very similar to the one obtained with the kernel $K$ method. By construction, the reconstruction of $G$ with $L$ involves $G_{n}(\omega)\sim g(\omega)^{2}L_{n-1}(\omega)$ so that the high frequency noise is expected to behave better with $L$ than with $K$ (the factor $g(\omega)^{2}$ effectively suppresses the high frequency). Fig. 4 shows a comparison of the errors obtained on $\Sigma^{R}_{4}(\omega)$ using the two methods. We find that the error using the $L$ method is essentially frequency independent while the error with the $K$ method depends strongly on frequency. In most cases the $L$ method is preferred but at small frequency, we have observed that the $K$ method can provide smaller error bars.

IV.3 Numerical results for the spectral function

Once the Green’s function or self energy has been obtained up to a certain order, the last task is to extract the physics information from this expansion. The most naive approach is to compute the truncated series up to a certain maximum order $N$,

[TABLE]

We find that the series has a convergence radius $U_{c}\approx 6\Gamma$ at the particle symmetry point $\epsilon_{d}=0$ while this convergence radius decreases down to $U_{c}\approx 4\Gamma$ in the asymmetric case $\epsilon_{d}=\Gamma$. These convergence radii fix the maximum strength of $U$ that one can study using the naive truncated series approach.

The data for the self energy (second and fourth panels) and corresponding spectral functions (first and third panels) are shown in Fig. 5 for the symmetric case (upper two panels, $\epsilon_{d}=0$ and $U=5\Gamma$) and asymmetric case (lower two panels, $\epsilon_{d}=\Gamma$ and $U=3\Gamma$). For these values of interaction, the error in our calculation is dominated by the finite truncation of the series (negligible error due to the statistical Monte-Carlo sampling) and is of the order of the line width. Fig. 5 also shows the NRG results that we use to benchmark our calculations and that are in excellent agreement with our data. The NRG calculation is the same as in the companion paperBertrand et al. (2019), where it is described in details. Note that in order to obtain this agreement, the precision of the NRG calculations had to be pushed much further than what is typically done in the field indicating that the QMC method is very competitive, in particular at large frequencies.

Qualitatively, the strength of interaction that could be reached using the truncated series corresponds to the onset of the Kondo effect: one observes in the upper panel of Fig. 5 that the Kondo peak starts to form around $\omega=0$, its width is significantly narrower than without interaction and the premisses of the side peaks at $\pm U/2$ can be seen. In order to observe well established Kondo physics, one must therefore go beyond the convergence radius wall. This is in fact rather natural, the convergence radius corresponds to poles or singularities in the complex $U$ plane which themselves correspond to the characteristic energy scales of the system. Getting past this “convergence radius wall” is crucial and is the subject of the companion article to the present manuscript.

V Conclusion

In this article, we have presented a quantum Monte-Carlo algorithm that allows one to calculate the out of equilibrium Green’s functions of an interacting system, order by order in powers of the interaction coupling constant $U$. We applied it to the Anderson model in the quantum dot geometry and obtained up to $10$ orders of the Green’s function and self-energy. A detailed benchmark was also presented against NRG computations, after a simple summation of the series at weak coupling. Our results were obtained at almost zero temperature, but we found that the method works equally well at finite temperature or out-of-equilibrium. It works equally well for transient response to an interaction quench or at long time where a stationary regime has been reached.

The method presented here has the great advantage to produce the perturbative expansion at infinite time, *i.e. * directly in the steady state. Its complexity is uniform in time: it does not grow at long time, contrary to other QMC approaches. The drawback, like any “diagrammatic” QMC, is that we have just produced the perturbative series of e.g. the Green’s function and the self-energy. At weak coupling, we can simply sum it, as shown earlier in the benchmark. However, at intermediate coupling, simply summing the series with partial sums will fail. Most quantities have a finite radius $R$ of convergence in $U$: for $U>R$, the series diverges. In Ref. Profumo et al., 2015, we showed how to use well-known conformal transformation resummation technique to solve this problem and obtain density of particle on the dot vs $U$ up to $U=\infty$. How to generalize this idea to make it work for real frequency Green’s functions, and also to control the amplification of stochastic noise due to such resummation will be addressed in a separate publication Bertrand et al. (2019).

VI Acknowledgement

We thank Laura Messio for discussions at the early stage of this work. We thank Serge Florens for interesting discussions and for providing NRG calculations. The Flatiron Institute is a division of the Simons Foundation. We acknowledge financial support from the graphene Flagship (ANR FLagera GRANSPORT), the French-US ANR PYRE and the French-Japan ANR QCONTROL.

Appendix A Proof of the clusterization property of the kernel $K$

In this appendix, we extend the proof of the clusterization property of Ref. Profumo et al., 2015 for the Kernel $K$, $\bar{K}$ and $L$. We want to show that, if some of the times $u_{i}$ are sent to infinity in the integral in Eq. (8), (21), (23) and (28), the sum under the integral vanishes (while each determinant taken individually does not). We will not try to prove here the stronger property that the integrals do indeed converge but we observe it empirically in the numerical computations.

Let us restart from the clusterization proof of Ref. Profumo et al., 2015 for Eq. (8) and examine the sum over the Keldysh indices:

[TABLE]

If some $u_{i}$ are sent to infinity, we can relabel them $u_{p+1},\dots,u_{n}$. Since $g$ vanishes at large time (due to the presence of the bath), the determinants in the sum become diagonal by block

[TABLE]

The upper-left determinant does not depend on $a_{p+1},\dots,a_{n}$, so we can apply Eq. (15) to the bottom-right determinant and the sum $S$ vanishes.

Let us now turn to the kernel $K$ defined in Eq. (21). The situation is slightly different. First, with a simple relabelling, we can restrict ourselves to the case $p=1$ in Eq. (21). Let us first split the $U$ into two subsets.

[TABLE]

Some $u_{i}$ go to infinity. We distinguish two cases.

If $u_{1}$ does not go to infinity, we can relabel the indices so that $u_{p+1},\dots,u_{n}$ go to infinity. 2. 2.

If $u_{1}$ goes to infinity, we can relabel the indices so that $u_{1},\dots,u_{p}$ go to infinity.

In both cases, the upper-right part of the matrix vanishes and we get a block-trigonal determinant

[TABLE]

since the first determinant does not depend on $a_{p+1},\dots,a_{n}$. The second term cancels because of (15).

Appendix B Expression of the kernel $L$ as a sum of Green’s functions

We show here that the kernel $L$ can be expressed in terms of Green’s functions. Starting from its definition Eq. (28), we follow the same steps as in Section II.5. We first use the fact (due to determinant symmetry) that all terms of the sum over $p$ have the same contribution:

[TABLE]

Then we sum out the Dirac delta:

[TABLE]

The pattern of a 3-particle Green’s function can be recognized:

[TABLE]

where $E$ is defined as:

[TABLE]