A Quantum Monte Carlo algorithm for out-of-equilibrium Green's functions at long times
Corentin Bertrand Olivier Parcollet Antoine Maillard Xavier Waintal




















    Univ. Grenoble Alpes, CEA, INAC-PHELIQS, GT F-38000 Grenoble, FranceCenter for Computational Quantum Physics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USAInstitut de Physique Théorique (IPhT), CEA, CNRS, UMR 3681, 91191 Gif-sur-Yvette, FranceLaboratoire de Physique de l’Ecole Normale Supérieure, UMR 8550, 24 Rue Lhomond, 75005 Paris, FranceUniv. Grenoble Alpes, CEA, INAC-PHELIQS, GT F-38000 Grenoble, France
    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.
    A Quantum Monte Carlo algorithm for out-of-equilibrium Green's functions at long times
        Xavier Waintal
        December 30, 2023
    ======================================================================================
    
    
    
    
    § 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 scale<cit.>, light-induced
    superconductivity <cit.> or metal-insulator transition driven by
    electric field <cit.>.
    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 <cit.> (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)<cit.>
    faces entanglement growth at long times.
    Early attempts of real time quantum Monte Carlo <cit.>
    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 <cit.>
    and the so-called “diagrammatic" Quantum Monte Carlo <cit.> (QMC).
    Diagrammatic QMC<cit.>
    computes the perturbative expansion of physical quantities in power of the interaction U,
    using an importance sampling Monte Carlo.
    In Ref. Profumo2015, 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 time<cit.> and in imaginary time<cit.>, which leads to much higher performance.
    
    
    In this article, we generalize the algorithm presented in Ref. Profumo2015 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 paper<cit.>, 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 <ref>, we introduce the necessary formalism and derive
    the basic equations that will be used to formulate the method. Section
    <ref> discusses our sampling strategy for the QMC algorithm, as well as
    relations with previous work.  Section <ref> shows our numerical
    data and the detailed comparison with our benchmarks in the weak coupling regime.
    
    
    
    
    
    
    § 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.
    
    
    
    
     §.§ Notations and main expansion formula
    
    We consider a generic class of system given by a time-dependent Hamiltonian of the form,
    
        H(t)    =   H_0(t)+U  H_ int(t) 
        H_0(t)    =   ∑_xyH^0_xy(t)c^†_xc_y
        H_ int(t)    =   ∑_xy V_xy(t) (c_x^†c_x - α_x) (c_y^†c_y - α_y)
    
    where H_0(t) is a quadratic unperturbed Hamiltonian. 
    The operators c^†_x (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.
    H_ 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. Profumo2015) 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 α, which has been
    introduced in previous works <cit.>.
    In particular, we have shown in Ref. Profumo2015 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. Gaury2014.
    
    
    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'∈{0,1} which provides the Green's function with a 2× 2 structure,
    
    G^aa'_xx'(t,t')=
    (
    [  G^T_xx'(t,t')  G^<_xx'(t,t');  G^>_xx'(t,t') G^T̅_xx'(t,t') ])
    
    where G^T_xx'(t,t'), G^<_xx'(t,t'), G^>_xx'(t,t') and G^T̅_xx'(t,t') 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'_xx'(t,t'). 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 ≡ (x, t, a).
    Using Keldysh points, we can rewrite the above definitions using the standard conventions of Keldysh formalism,
    
        G_xx'^aa'(t,t') 
        ≡ -i < T_ cc(x,t,a) c^†(x',t',a') >
    
    where the creation (c^†(X')) and annihilation operators
    (c(X) or c(x,t,a)), here in Heisenberg representation, are
    ordered using the contour time ordering operator T_ c. T_ 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
    
        G(X, X') ≡ G[(x,t,a),(x',t',a')]≡ G_xx'^aa'(t,t')
    
    and we will also note the δ function on the Keldysh contour as 
    
        δ_c[((x,t,a),(x',t',a') ] ≡δ(t-t') δ_aa'δ_xx'.
    
    
    Using the above notations (with ħ = 1), one can derive the usual
    expansion in power of U using Wick's theorem. We first assume
    α_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 α_x  0.
    We obtain<cit.>:
    
        G_xx'^aa'(t,t') = ∑_n≥ 0i^n U^n/n!∫∏_k=1^n du_k
        ∑_{x_k, y_k}∑_{a_k}×
        (
        ∏_k=1^n 
        (-1)^a_k V_x_k y_k(u_k)
        )
        (x,t,a), U_1, …, U_2n(x',t',a'), U_1, …, U_2n
    
    which forms the starting point of this work. Here, we have noted for 1≤ k ≤ n
    
        U_2k-1   = (x_k, u_k, a_k) 
        
            U_2k   = (y_k, u_k, a_k)
    
    and introduced the 
     Wick matrix:
    
        A_1, …, A_mB_1, …, B_m_ij≡[ g(A_1, B_1)             g(A_1, B_m);           ⋮           ⋱           ⋮; g(A_m, B_1)             g(A_m, B_m) ]_ij
    
    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,
    
        A_1, …, A_mB_1, …, B_m≡A_1, …, A_mB_1, …, B_m_ij
    
    In equation (<ref>), we start at t<0 with a non-interacting state and switch on the interaction for t≥ 0.
    Hence the time integrals in
    Eq. <ref> run over the segment [0, t_M], where t_M = max(t,
    t'). 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 <ref>, we will a fixed large value of t_M. 
    We emphasize that the complexity of the algorithm does  not grow with t_M.
    Eq. (<ref>) 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. (<ref>) 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. (<ref>), i.e.
        g(U_2k,U_2k)    =g^<_y_ky_k(u_k,u_k)
        
                g(U_2k-1,U_2k-1)   =g^<_x_kx_k(u_k,u_k)
         
                g(U_2k-1,U_2k)   =g^>_x_ky_k(u_k,u_k)
        
                g(U_2k,U_2k-1)   =g^<_y_kx_k(u_k,u_k).
    
    
    To proceed to the general case α_x  0, one only needs to shift the diagonal terms of the Wick matrix using the following replacement rules,
    as shown in Ref.Profumo2015, Appendix A:
    
        g(U_2k,U_2k)    → g(U_2k,U_2k) - iα_y_k
        
                g(U_2k-1,U_2k-1)    → g(U_2k-1,U_2k-1) - iα_x_k
    
    As a result, all derivations in this paper can be done by first ignoring
    α_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. (<ref>) can be extended to the
    calculations of arbitrary Green's functions. 
    The rule for doing so is as follows: whenever one introduces a product
    -ic(Y)c^†(Y') under the time ordering operator in
    Eq. (<ref>), one must add the corresponding Keldysh points in
    the Wick determinant of Eq. (<ref>):
    
        X, U_1, …, U_2nX', U_1, …, U_2n→X, Y, U_1, …, U_2nX', Y', U_1, …, U_2n
    
    
    If Y and Y' share the same time and orbital index y, we have the possibility to introduce terms of the form
    -i[c(Y)c^†(Y') + α_y] in the definition of the Green's
    function. In that case, one must replace g(Y,Y')→ g(Y,Y') - iα_y
     in the diagonal of the Wick matrix.
    Again, to improve readability, we will keep this replacement implicit
    in Wick matrices, but explicit otherwise.
    
    
    
    
     §.§ 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. (<ref>) leaves the Wick determinant unchanged up to a factor (-1). 
    An important property
    of the formalism, as already shown in Ref. Profumo2015, is that for n>0:
    
        ∑_{a_k}
        (-1)^∑_k=1^n a_kU_1, …, U_2nU_1, …, U_2n = 0
    
    for any times u_1, …, u_n and orbital indices x_1, …, x_n and y_1, …, 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
    
        A_1,…, A_mB_1,…,B_m = 
        ∑_k=1^m (-1)^k+1g(A_k,B_1) A_1, …A_k…, A_mB_1, B_2, …, B_m
    
    for the expansion along the first column and
    
        A_1,…, A_mB_1,…,B_m= 
        ∑_k=1^m (-1)^k+1g(A_1,B_k) A_1, ……, A_mB_1, B_k, …, B_m
    
    for the expansion along the first row. The notation A_k (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,
    
        (-1)^i+jA_1…A_i, ……, A_mB_1,…B_j, …, B_m =  
        A_1,…, A_mB_1,…,B_m^-1_jiA_1,…, A_mB_1,…,B_m
    
    
    
     §.§ Definition of the kernel K for the one-body Green's function
    
    
    In Ref. Profumo2015, a QMC scheme was defined directly on Eq. (<ref>) so that
    a single QMC run could provide the value of G_xx'^aa'(t,t') for a single pair of times t and t'.
    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. (<ref>) on Eq. (<ref>), we
    obtain
    
        G_xx'^aa'(t,t') = g_xx'^aa'(t,t') + ∑_n≥ 1i^n U^n/n!∫∏_k=1^n du_k
        ∑_{x_k, y_k}∑_{a_k}(
        ∏_k=1^n 
        (-1)^a_k V_x_k y_k(u_k)
        )
        ×
        (
           ∑_p=1^2n
        (-1)^p 
        g[ (x,t,a), U_p  ]
         U_1, …, U_2n (x',t',a'), U_1, …, U_p,…, U_2n
        +
        g[  (x,t,a), (x',t',a')  ]
         U_1, …, U_2n U_1, …, U_2n)
    
    The last term of the sum vanishes for n>0 due to Eq. (<ref>). Factorizing the g from the sum, we arrive at
    
        G_xx'^aa'(t,t') = g_xx'^aa'(t,t')  +  
        ∫ du ∑_b, y (-1)^b
        g_xy^ab(t,u) K_yx'^ba'(u,t')
    
    where the kernel K_yx'^ba'(u,t')=K(Y,X') with Y=(y,u,b) is defined by 
    
    
        K(Y,X') ≡ (-1)^b
        ∑_n≥ 1i^n U^n/n!∫∏_k=1^n du_k
        ∑_{x_k, y_k}∑_{a_k}(
        ∏_k=1^n (-1)^a_k V_x_k y_k(u_k)
        )
        ∑_p=1^2n
        (-1)^pδ_c(Y, U_p) 
         U_1, …, U_2n X', U_1, …, U_p,…, U_2n
    
    Equations (<ref>) and (<ref>) will be the basis of one of the method developed
    in this article. Eq. (<ref>) will provide the mean to get a full t-curve in a single calculation and
    Eq. (<ref>) to relate the corresponding kernel to the Green's function G, the target of the calculation.
    
    A symmetric kernel K̅ may be derived following the exact same route but now expanding the
    Wick determinant along the first  column using  Eq. (<ref>). We find
    
    
        G_xx'^aa'(t, t') = g_xx'^aa'(t, t') +  
        ∫ du ∑_b, y (-1)^b
        K̅_xy^ab(t,u) g_yx'^ba'(u,t')
    
    where the kernel K̅ is defined by 
    
    
        K̅(X,Y) ≡ (-1)^b
        ∑_n≥ 1i^n U^n/n!∫∏_k=1^n du_k
        ∑_{x_k, y_k}∑_{a_k}(
        ∏_k=1^n (-1)^a_k V_x_k y_k(u_k)
        )
        ∑_p=1^2n
        (-1)^pδ_c(Y,U_p) 
        X, U_1, …, U_p,…, U_2nU_1, …, U_2n
    
    
    
     §.§ 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, 
    
        F_xx'z^aa'(t,t') 
        ≡ (-i)^2 < T_ cc(x,t,a) c^†(x',t',a')[c^†(z,t',a')c(z,t',a') -α_z] >
    
    
    In the next paragraph, we shall prove that F is essentially equal to 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)<cit.>
    as well as in imaginary time QMC methods as an  improved estimator<cit.>.
    
    The expansion of F reads, 
    
        F_xx'z^aa'(t,t') =
        - ∑_n≥ 0i^nU^n/n!∫∏_k=1^n du_k
        ∑_{x_k, y_k}∑_{a_k}(
        ∏_k=1^n 
        (-1)^a_k
        V_x_k y_k(u_k)
        )
        (x,t,a), (z,t',a'), U_1,…, U_2n (x',t',a'), (z,t',a'), U_1, …, U_2n
    
    To obtain the kernel of F, we expand the determinant along its first row using Eq. (<ref>),
    
        F_xx'z^aa'(t,t') = - ∑_n≥ 0i^n U^n/n!∫∏_k=1^n du_k
        ∑_{x_k, y_k}∑_{a_k}(
        ∏_k=1^n 
        (-1)^a_k V_x_k y_k(u_k)
        )
        ×
        (
           ∑_p=1^2n
        (-1)^p+1 
        g( (x,t,a), U_p  )
         (z,t',a'), U_1, …, U_2n (x',t',a'), (z,t',a'), U_1, …, U_p,…, U_2n
        + 
         .
        
         .
        g( (x,t,a), (x',t',a')  )
         (z,t',a'),U_1, …, U_2n (z,t',a'),U_1, …, U_2n
        -
        g( (x,t,a), (z,t',a')  )
         (z,t',a'), U_1, …, U_2n (x',t',a'), U_1, …, U_2n)
    
    Identifying the two last terms with the corresponding expansion of G, we arrive at,
    
        F_xx'z^aa'(t,t') = -g_xx'^aa'(t,t') [G_zz^<(t',t') - iα_z] + g_xz^aa'(t,t') G_zx'^<(t',t') -
        ∫ du ∑_b, y (-1)^b
        g_xy^ab(t,u) L_yx'z^ba'(u,t')
    
    where the kernel L is defined by 
    
    
        L_yx'z^ba'(u,t') ≡ (-1)^b
        ∑_n≥ 1i^n U^n/n!∫∏_k=1^n du_k
        ∑_{x_k, y_k}∑_{a_k}(
        ∏_k=1^n (-1)^a_k V_x_k y_k(u_k)
        )
        ×
        ∑_p=1^2n
        (-1)^p+1δ_c((y,u,b), U_p) 
         (z,t',a'), U_1, …, U_2n (x',t',a'), (z,t',a'), U_1, …, U_p,…, U_2n
    
    
    
     §.§ Relation between F, 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 K̅, Eq. (<ref>).
    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 ↔ U_2k-1 and U_2 ↔ U_2k in the sums 
    and integration. We arrive at,
    
    
        K̅_xy^ab(t,u) = (-1)^b
        ∑_n≥ 1i^n U^n/n!∫ d u_1 ∑_x_1, y_1∑_a_1 (-1)^a_1 V_x_1, y_1(u_1)
        ∫∏_k=2^n du_k
        ∑_{x_k, y_k} k≥ 2∑_{a_k} k≥ 2(
        ∏_k=2^n (-1)^a_k V_x_k y_k(u_k)
        )
        ×
        
        2n  δ_c(U_1, (y,u,b))
        (x,t,a), U_2, U_3,…, U_2n U_1, U_2, U_3, …,  U_2n
    
    We can now perform explicitly the integral over u_1 where the delta function yields, for u ∈ [0, t_M] (K̅ is zero otherwise):
    
    
        K̅_xy^ab(t,u)    =  2iU
        ∑_n≥ 1i^n-1U^n-1/(n-1)!∑_y_1  V_y, y_1(u)
        ∫∏_k=2^n du_k
        ∑_{x_k, y_k}{a_k}k ≥ 2(
        ∏_k=2^n 
        (-1)^a_k
        V_x_k y_k(u_k)
        )
        (x,t,a), (y_1,u,b), U_3, …, U_2n (y,u,b), (y_1,u,b), U_3, …, U_2n
        K̅_xy^ab(t,u)    =  2iU
        ∑_z V_y, z(u)
        ∑_n≥ 0i^nU^n/n!∫∏_k=1^n du_k
        ∑_{x_k, y_k}{a_k}(
        ∏_k=1^n 
        (-1)^a_k
        V_x_k y_k(u_k)
        )
        (x,t,a), (z,u,b), U_1,…, U_2n (y,u,b), (z,u,b), U_1, …, U_2n
        K̅_xy^ab(t,u)    = - 2iU ∑_z V_yz(u) F^ab_xyz(t,u)
    
    This shows the kernel K̅ is no more than a sum of 2-particle Green's
    functions. The relation between K̅ and G
    in Eq. (<ref>) can then be transformed into:
    
        G_xx'^aa'(t, t') = g_xx'^aa'(t, t') -
        2 i U  ∫ du ∑_b, y (-1)^b
        ∑_z  V_yz(u) F_xyz^ab(t,u) g_yx'^ba'(u,t')
    
    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 <ref> 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.
    
    
    
    
     §.§ 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, K̅ or L.
    
    The retarded and advanced Green's functions relate to the lesser and greater Green's functions as follows:
    
        G^R_xx'(t,t') = θ(t-t') ( G^>_xx'(t,t') - G^<_xx'(t,t') ) 
           G^A_xx'(t,t') = θ(t'-t) ( G^<_xx'(t,t') - G^>_xx'(t,t') )
    
    where θ is the Heaviside function.
    From the definitions of the time-ordered and time-antiordered Green's functions, these can also be written as:
    
        G^R_xx'(t,t') = G^a0_xx'(t,t') - G^a1_xx'(t,t') 
           G^A_xx'(t,t') = G^0a_xx'(t,t') - G^1a_xx'(t,t')
    
    where a can be any Keldysh index. These are also valid for the non-interacting g.
    As K̅ is a sum of Green's function, one may define in the same way a
    retarded version of K̅, denoted K̅^R. We can see from
    Eq. (<ref>) and the definition of F in
    Eq. (<ref>) that K̅^R follows the same properties:
    
        K̅^R_xx'(t,t') = K̅^a0_xx'(t,t') - K̅^a1_xx'(t,t')
    
    
    We now show a simple relation between G^R, g^R and K̅^R. Plugging Eq. (<ref>) into Eq. (<ref>), one gets:
    
        G^R_xx'(t,t') = g^R_xx'(t,t')
            + ∫ du ∑_y (
                K̅^00_xy(t,u) [ g^00_yx'(u,t') - g^01_yx'(u,t') ]
                - K̅^01_xy(t,u) [ g^10_yx'(u,t') - g^11_yx'(u,t') ]
            )
    
    This simplifies into:
    
        G^R_xx'(t,t') = g^R_xx'(t,t')
            + ∫ du ∑_y
            K̅^R_xy(t,u) g^R_yx'(u,t')
    
    
    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. (<ref>) and Eq. (<ref>).
    As all of them are sums of Green's functions, one may show that they all verify similar properties as in Eq. (<ref>) and (<ref>).
    Then from Eq. (<ref>) follows:
    
        G^A_xx'(t,t') = g^A_xx'(t,t')
            + ∫ du ∑_y
            g^A_xy(t,u) K^A_yx'(u,t')
    
    and from Eq. (<ref>) follows:
    
        F^A_xx'z(t,t') = -g^A_xx'(t,t')[G^<_zz(t',t') - iα_z]
            + g^A_xz(t,t')G^<_zx'(t',t)
            - ∫ du ∑_y g^A_xy(t,u) L^A_yx'z(u,t')
    
    
    
    § 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. Profumo2015 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.
    
    
    
    
     §.§ Sampling of the kernel K
    
    Let us first discuss the calculation of G using the kernel K, using Eq. (<ref>) and (<ref>).
    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. Profumo2015. The kernel takes the form, 
    
        K(Y,X')    = (-1)^b
        ∑_n≥ 1∫∏_k=1^n du_k
        ∑_{x_k, y_k}∑_p=1^2n∑_a_p
        (-1)^a_pδ_c(Y, U_p) 
        W^n_p(X',{U_k}, a_p)
        
        
        W^n_p(X',{U_k}, a_p)
           ≡i^n U^n/n!(∏_k=1^n  V_x_k y_k(u_k) )
        ∑_{a_k} k≠ p( ∏_k≠ p (-1)^a_k) 
        (-1)^p
         U_1, …, U_2n (x',t',a'), U_1, …, U_p,…, U_2n
    
    We define a configuration C as 
    
      *  The order n.
       
      *  A set of times {u_k ∈ [0, t_M]} for 1≤ k ≤ n.
       
      *  Two sets of indices {x_k} and {y_k} for 1≤ k ≤ 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. (<ref>), t_M can be chosen to be any value
    larger than t and t' of the target Green's functions G(t, t').
    
    
    The kernel takes the form,
    
        K(Y,X') = (-1)^b
        ∑_ C (
        ∑_p=1^2n∑_a_p
           (-1)^a_pδ_c(Y, U_p) 
           W^n_p(X', {U_k}, a_p)  
        )
    
    where the sum over the configurations C is a compact notation for the sum and integrals of Eq.(<ref>).
    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( C) that will provide the (unnormalized) probability for the configuration C
    to be visited by the QMC algorithm. We define this weight as,
    
    
        W( C) = ∑_p=1^2n∑_a=0,1 | W^n_p(X',{U_k}, a)  |
    
    Noting Z_ QMC≡∑_ C W( C), the kernel can be rewritten as,
    
        K(Y,X') = (-1)^b Z_ QMC⟨∑_p=1^2n∑_a_p
           (-1)^a_pW^n_p(X', {U_k}, a_p)W( C)δ_c((y,u,b), U_p) 
        ⟩
    
    where the average is taken over the distribution W( C)/Z_ QMC.
    Z_ QMC is an effective partition function associated to the QMC algorithm.
    Note that by construction, the weight W( C) controls the
    measurement,  i.e.|W^n_p(X',{U_k}, a)| ≤ W( C). This is an
    essential property for a reweighting technique since it guarantees that the
    weight W^n_p(X',{U_k}, a)/W( 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. Shi2016).
    
    To sample W( C) and evaluate Z_ QMC, we use the continuous time
    Monte-Carlo technique that was discussed in details in Ref. Profumo2015.
    We use moves that change the order n by ± 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 C provides 2n values of Y=(y,u,b) which are
    recorded by binning with the weight of Eq. (<ref>).
    
    The partial weights W_p^n possess an essential  clusterization property,
    which generalizes the one discovered in Ref. Profumo2015: if one or
    several times u_k goes to infinity ( i.e. is far from
    t'), then all W_p^n goes to 0. In other words the integrand is localized
    around t'. A detailed proof is provided in Appendix <ref>. The point X'= (x',t',a') 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,
    
        W^n_p(X',{U_k}, a_p)
           =
        -i^n U^n/n!(∏_k=1^n  V_x_k y_k(u_k))
        ∑_{a_k} k≠ p( ∏_k≠ p (-1)^a_k) 
         X, U_1, …, U_2n X', U_1 ,…, U_2n_p1^-1× X, U_1, …, U_2n X', U_1 ,…, U_2n
    
    This last form is very convenient since a single Wick matrix (and its inverse) need to be stored and monitored during the calculation.
    
    
    
    
     §.§ Sampling of the kernel L
    
    
    Following the same route for L as was done for K, we can write:
    
        L_yx'z^ba'(u,t')    = (-1)^b
        ∑_n≥ 1∫∏_k=1^n du_k
        ∑_{x_k, y_k}∑_p=1^2n∑_a_p (-1)^a_pδ_c((y,u,b), U_p)
        W_p+2^n(X', {U_k}, a_p, z) 
        
        W_p+2^n(X', {U_k}, a_p, z)
           ≡i^n U^n/n!(∏_k=1^n  V_x_k y_k(u_k) )
        ∑_{a_k} k≠ p( ∏_k≠ p (-1)^a_k)
        (-1)^p+1 (z,t',a'), U_1, …, U_2n (x',t',a'), (z,t',a'), U_1, …, U_p,…, U_2n
    
    thus defining W_p^n for p = 3, 4, …, 2n+2. We also define W_1^n and W_2^n in the following way:
    
        W_1^n(X', {U_k}, z)
           ≡
        -i^n U^n/n!∑_{a_k}( ∏_k=1^n (-1)^a_k
        V_x_k y_k(u_k)  )
         (z,t',a'), U_1, …, U_2n(z,t',a'), U_1, …, U_2n
        
        W_2^n(X', {U_k}, z)
           ≡i^n U^n/n!∑_{a_k}( ∏_k=1^n (-1)^a_k
        V_x_k y_k(u_k)  )
         (z,t',a'), U_1, …, U_2n(x',t',a'), U_1, …, U_2n
    
    These two extra values are necessary to compute G^<_zz(t',t') and G^<_zx'(t',t'), which are needed to obtain F, as can be seen in Eq. (<ref>).
    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. (<ref>), the determinant within any W_p^n (for p = 1, …, 2n+2) can be replaced by:
    
        (-1)^p+1 X, (z,t',a'), U_1, …, U_2n X', (z,t',a'), U_1 ,…, U_2n_p1^-1× X, (z,t',a'), U_1, …, U_2n X', (z,t',a'), U_1 ,…, U_2n
    
    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 C is now:
    
        W( C) = | W^n_1(X',{U_k}, z) | 
            + | W^n_2(X',{U_k}, z) | 
            + ∑_p=1^2n∑_a=0,1| W^n_p+2(X',{U_k}, a, z) |
    
    We define again Z_ QMC≡∑_ C W( C) (which however has a different value than in the previous section).
    Finally, L can be written as:
    
        L_yx'z^ba'(u,t') = (-1)^b Z_ QMC⟨∑_p=1^2n∑_a_p
           (-1)^a_pW^n_p+2(X', {U_k}, a_p, z)W( C)δ_c((y,u,b), U_p)
        ⟩
    
    and, from Eq. (<ref>), we get that the values of G^<_zz(t',t') and G^<_zx'(t',t') (needed to compute F) are:
    
        G^<_zz(t',t') =   
            - Z_ QMC< W_1^n(X', {U_k}, z)W( C)>
            
        
            G^<_zx'(t',t') =   
            Z_ QMC< W_2^n(X', {U_k}, z)W( C)>
    
    The Monte-Carlo algorithm used to evaluate these averages is the same as in the previous section, except for the weight W( C).
    
    
    
    
    
     §.§ A discussion of the Werner et al. approach<cit.>
    
    
    In this paragraph, we discuss the relation of this work with a preceding 
    work <cit.> that also implements an expansion in powers of U within the Keldysh formalism. 
    Although both results are consistent, Ref. Werner2010 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. Werner2010 corresponds to a maximum time of
    t_M≈ 5/Γ between the switching of the interaction and the
    measurement of the observable while we found that t_M≈ 20/Γ is
    needed to enter the stationary result at order n=8 as shown in
    Fig. <ref>. 
    A direct consequence of this issue is that Ref. Werner2010 cannot access
    the small bias regime where the Kondo effect is present: since the Kondo
    temperature T_K is typically much smaller than Γ, long simulation
    times t_M ≫ 1/T_K are needed to capture the Kondo physics. 
    Second, the technique of Ref. Werner2010 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. Werner2010 is that some results could be
    obtained in regimes where the “sign" of the Monte-Carlo calculation was very
    small ∼ 10^-3 (see for instance Fig.5 of Ref. Werner2009). 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. (<ref>) which provides the expansion for the Keldysh Green's function G_xx'^aa'(t,t').
    A similar formula<cit.> provides the sum of vacuum diagrams, sometimes called the Keldysh “partition function" Z,
    
        Z ≡∑_n≥ 0i^n U^n/n!∫∏_k=1^n du_k
        ∑_{x_k, y_k}(
        ∏_k=1^n  V_x_k y_k(u_k)
        )
        ×
        ∑_{a_k}∏_k=1^n (-1)^a_k U_1, …, U_2nU_1, …, U_2n
    
    We have Z=1 in the Keldysh formalism, reflecting the unitarity of  quantum mechanics.
    We see that (<ref>) 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. (<ref>) (without the sum over Keldysh indices):
    
        Z({U_i}) ≡i^n U^n/n!(
        ∏_k=1^n (-1)^a_k V_x_k y_k(u_k)
        )
         U_1, …, U_2nU_1, …, U_2n
    
    Ref. Werner2009 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', {U_i})
    the integrand of Eq. (<ref>) and Z_ QMC the
    integral of |Z({U_i})|:
    
        Z_ QMC≡∑_n≥ 0∫∏_k=1^n du_k
        ∑_{x_k, y_k}∑_{a_k}
        |Z({U_i})|
    
    The weight of the QMC is |Z({U_i})|/Z_ QMC. We have 
    
        G_xx'^aa'(t,t') = ⟨ G(X,X',{U_i})/|Z({U_i})|⟩/⟨Z({U_i})/|Z({U_i})|⟩
    
    where the average is the Monte-Carlo average.
    The denominator of Eq. (<ref>) is the QMC sign mentioned above. From
    Z=1, we find that this QMC sign is simply given by 
    ⟨ Z({U_i})/|Z({U_i})|⟩ = 1/Z_ QMC.
    
    
    
    Now, we note that the probability to visit order 0 is also |Z(∅)|/Z_ QMC = 1/Z_ QMC. 
    Therefore the average sign in the denominator of Eq. (<ref>) 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. Werner2009.
    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. Profumo2015.
    A genuine sign problem can however
    results from the numerator of Eq. (<ref>).
    
    
    
    
    § 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 package<cit.>.
    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 <cit.>.
    We stress that the QMC technique
    is not restricted to impurity models and also applies to lattice models such as the Hubbard model.  
    
    
    
    
     §.§ Definition of the model
    
    In the Anderson impurity model, the impurity is described by the operators d^†_σ (d_σ) that create (destroy) 
    an electron on the impurity with spin σ. 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,
    
        g^R_σσ'(ω) = 1/ω - ϵ_d + iΓδ_σσ'
    
    where the parameter Γ sets the width of the resonance in the absence of
    electron-electron interactions and the on-site energy ϵ_d sets the
    resonance with respect to the Fermi level. Eq. (<ref>) fully
    defines the model at equilibrium. The model is made non trivial through the
    interacting terms that reads,
    
        H_ int
        = Uθ(t)(n_↑-1/2)(n_↓-1/2).
    
    where n_σ = d^†_σd_σ is the impurity electronic density of spin σ
    and the Heaviside function θ(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_ BT =
    10^-4Γ, 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' only and can be studied in the frequency domain. 
    
    G^R_σσ'(t - t') = δ_σσ'∑_n=0^+∞ G_n(t-t') U^n
    
    from which one can obtain the corresponding quantity in the frequency domain by fast Fourier transform,
    
    G^R_σσ'(ω) = δ_σσ'∑_n=0^+∞ G_n(ω) U^n
    
    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) 
    
    A(ω) = -1/π Im[G^R(ω)].
    
    and the retarded self energy Σ^R(ω),
    
        G^R(ω) = 1/ω - ϵ_d + iΓ - Σ^R(ω)
    
    
    
     §.§ Numerical Results order by order
    
    We now present the numerical data obtained by sampling the kernel K. The left panels of Fig. <ref> 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
    (50 000 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 ∫_0^t K^R(u)du is however
    noiseless as can be seen from  the example shown in
    Fig. <ref> 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(ω) is relatively noisy at high frequency. Last, we obtain
    G^R_n(ω)= g^R(ω) K^R_n(ω) for n≥ 1 (from Eq. (<ref>)) as shown in the right panel of
    Fig. <ref>. The factor g^R(ω), 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. Gull2008
    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(ω) 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(ω), 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. Profumo2015.
    
    Using the definition Eq. (<ref>) of the Self energy, we can obtain a recursive expression for Σ^R_n in term of the Green's function expansion:
    
        Σ^R_n(ω) = [g^R(ω)]^-2 G^R_n(ω) - ∑_k=1^n-1Σ^R_k(ω) G^R_n-k(ω) g^R(ω)^-1
    
    for n>1 with Σ^R_1(ω) = [g^R(ω)]^-2 G^R_1(ω). The corresponding data is shown in Fig. <ref> where we plot the coefficients
    Σ_n^R(ω) 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 Σ^R_n(ω) 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 Yamada<cit.>. The result at order 2 is compared with the result of Yamada in the left panel 
    of  Fig. <ref> 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,
    
        Σ^R(U,ω) =Γ∑_n,m i^m+1 s_n,m(ω/Γ)^m  (U/Γ)^n
    
    Table <ref> 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(ω) ∼ g(ω)^2 L_n-1(ω) so 
    that the high frequency noise is expected to behave better with L than with K (the factor g(ω)^2 effectively suppresses the high frequency).
    Fig. <ref> shows a comparison of the errors obtained on Σ^R_4(ω) 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. 
    
    
    
    
     §.§ 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,
    
        Σ^R(U,ω) ≈∑_n=1^N Σ^R_n(ω) U^n
    
    We find that the series has a convergence radius U_c ≈ 6Γ at the particle symmetry point ϵ_d=0 while this
    convergence radius decreases down to U_c ≈ 4Γ in the asymmetric case ϵ_d=Γ. 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. <ref> for the symmetric
    case (upper two panels, ϵ_d= 0 and U= 5Γ) and asymmetric case (lower two panels, ϵ_d= Γ and U= 3Γ). 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. <ref> 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 paper<cit.>, 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. <ref> that the Kondo peak starts to form
    around ω = 0, its width is significantly narrower than without
    interaction and the premisses of the side peaks at ± 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.
    
    
    
    
    § 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. Profumo2015, 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=∞.
    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 <cit.>.
    
    
    
    
    § 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.
    
    
    
    
    § PROOF OF THE CLUSTERIZATION PROPERTY OF THE KERNEL K
    
    
    In this appendix, we extend the proof of the clusterization property of Ref. Profumo2015
    for the Kernel K, K̅ and L.
    We want to show that, if some of the times u_i are sent to infinity in the integral in Eq. (<ref>),
    (<ref>), (<ref>) and (<ref>), 
    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. Profumo2015 for
    Eq. (<ref>) and examine the sum over the Keldysh indices:
    
        S ≡∑_{a_k}∏_k=1^n(-1)^ a_k(x,t,a), U_1, …, U_2n(x',t',a'), U_1, …, U_2n
        =
        ∑_{a_k}∏_k=1^n(-1)^ a_k|
        g(X, X')    g(X, U_1)       g(X, U_2n) 
        g(U_1, X')    g(U_1, U_1)       g(U_1, U_2n) 
        ⋮   ⋮   ⋱   ⋮
        g(U_2n, X')    g(U_2n, U_1)       g(U_2n, U_2n)
        |
    
    If some u_i are sent to infinity, we can relabel them u_p+1, …, u_n.
    Since g vanishes at large time (due to the presence of the bath), 
    the determinants in the sum become diagonal by block
    
        S ≈∑_{a_k}∏_k=1^n(-1)^ a_k|
        
            g(X, X')            g(X, U_1)               g(X, U_2p)           0                            0 
        
            g(U_1, X')          g(U_1, U_1)             g(U_1, U_2p)         0                            0
        ⋮   ⋮   ⋱   ⋮   ⋮   ⋱   ⋮
        
            g(U_2p, X')       g(U_2p, U_1)          g(U_2p, U_2p)      0                            0 
        
            0                   0                       0                      g(U_2p+1, U_2p+1)        g(U_2p+1, U_2n) 
        ⋮   ⋮   ⋱   ⋮   ⋮   ⋱   ⋮
        
            0                   0                       0                      g(U_2n, U_2p+1 )         g(U_2n, U_2n) 
        |
    
    The upper-left determinant does not depend on a_p+1, …, a_n, so we can apply Eq. (<ref>) to the bottom-right determinant and 
    the sum S vanishes.
    
    Let us now turn to the kernel K defined in Eq. (<ref>).
    The situation is slightly different.
    First, with a simple relabelling, we can restrict ourselves to the case p=1 in Eq. (<ref>).
    Let us first split the U into two subsets. 
    
    
        S    =
        ∑_{a_k}∏_k=1^n(-1)^ a_k|
        
            g(U_1, X')          g(U_1, U_2)             g(U_1, U_2p)         g(U_1, U_2p+1)             g(U_1, U_2n)  
        ⋮   ⋮   ⋱   ⋮   ⋮   ⋱   ⋮
        
            g(U_2p, X')       g(U_2p, U_2)          g(U_2p, U_2p)      g(U_2p, U_2p+1)          g(U_2p, U_2n) 
        
            g(U_2p+1, X')     g(U_2p+1, U_2)        g(U_2p+1, U_2p)    g(U_2p+1, U_2p+1)        g(U_2p+1, U_2n) 
        ⋮   ⋮   ⋱   ⋮   ⋮   ⋱   ⋮
        
            g(U_2n, X')       g(U_2n, U_2)          g(U_2n, U_2p)      g(U_2n, U_2p+1 )         g(U_2n, U_2n) 
        |
    
    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, …, u_n go to infinity. 
       
      * 
       If u_1 goes to infinity, we can relabel the indices so that u_1, …, u_p go to infinity.
    
    In both cases, the upper-right part of the matrix vanishes and we get a block-trigonal determinant  
    
        S    ≈∑_{a_k}∏_k=1^n(-1)^ a_k|
        
            g(U_1, X')          g(U_1, U_2)             g(U_1, U_2p)         0                            0 
        ⋮   ⋮   ⋱   ⋮   ⋮   ⋱   ⋮
        
            g(U_2p, X')       g(U_2p, U_2)          g(U_2p, U_2p)      0                            0
        
            g(U_2p+1, X')     g(U_2p+1, U_2)        g(U_2p+1, U_2p)     g(U_2p+1, U_2p+1)        g(U_2p+1, U_2n) 
        ⋮   ⋮   ⋱   ⋮   ⋮   ⋱   ⋮
        
            g(U_2n, X')       g(U_2n, U_2)          g(U_2n, U_2p)      g(U_2n, U_2p+1 )         g(U_2n, U_2n) 
        |
        
           = 
        (
        ∑_a_1,…,a_p∏_k=1^p(-1)^ a_kU_1, …, U_2pX', U_2, …, U_2p)
        ×(
        ∑_a_p+1,…,a_n∏_k=p+1^n(-1)^ a_kU_2p+1, …, U_2nU_2p+1, …, U_2n)
    
    since the first determinant does not depend on a_p+1, …, a_n.
    The second term cancels 
    because of (<ref>).
    
    
    
    
    
    
    
        
        
        
        
        
        
        
    
    
    
    
    
    
        
        
        
        
        
        
        
    
    
    
    
    
    § 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. (<ref>), we follow
    the same steps as in Section <ref>.
    We first use the fact (due to determinant symmetry) that all terms of the sum
    over p have the same contribution:
    
        L_yx'z^ba'(u,t') = (-1)^b
        ∑_n≥ 1i^n U^n/n!∫∏_k=1^n du_k
        ∑_{x_k, y_k}∑_{a_k}(
        ∏_k=1^n (-1)^a_k V_x_k y_k(u_k)
        )
        ×
        
        2n δ_c((y,u,b), U_1)
         (z,t',a'), U_1, U_2, …, U_2n (x',t',a'), (z,t',a'), U_2, …, U_2n
    
    Then we sum out the Dirac delta:
    
        L_yx'z^ba'(u,t') = 2iU
        ∑_z' V_yz'(u)
        ∑_n≥ 0i^n U^n/n!∫∏_k=1^n du_k
        ∑_{x_k, y_k}∑_{a_k}(
        ∏_k=1^n (-1)^a_k V_x_k y_k(u_k)
        )
        ×
         (z,t',a'), (y,u,b), (z',u,b), U_1, …, U_2n (x',t',a'), (z,t',a'), (z',u,b), U_1, …, U_2n
    
    The pattern of a 3-particle Green's function can be recognized:
    
        L_yx'z^ba'(u,t') = 2iU
        ∑_z' V_yz'(u)
        E_yx'zz'^ba'(u, t')
    
    where E is defined as:
    
        E_yx'zz'^ba'(u, t') ≡ (-i)^3 < T_ cc(y,u,b)c^†(x',t',a') 
            [ c^†(z,t',a')c(z,t',a') - α_z ] 
            [c^†(z',u,b)c(z',u,b) - α_z'] 
            >
    
    apsrev
