Reconstructing non-equilibrium regimes of quantum many-body systems from the analytical structure of perturbative expansions
Corentin Bertrand Serge Florens Olivier Parcollet Xavier Waintal




















    Univ. Grenoble Alpes, CEA, INAC-PHELIQS, GT  F-38000 Grenoble, FranceUniv. Grenoble Alpes, CNRS, Institut Néel, 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, FranceUniv. Grenoble Alpes, CEA, INAC-PHELIQS, GT F-38000 Grenoble, France
    We propose a systematic approach to the non-equilibrium dynamics of strongly
    interacting many-body quantum systems, building upon the standard perturbative expansion 
    in the Coulomb interaction. High order series are derived from the Keldysh version of 
    determinantal diagrammatic Quantum Monte Carlo, and the reconstruction beyond the weak 
    coupling regime of physical quantities is obtained by considering them as analytic 
    functions of a complex-valued interaction U. Our advances rely on two crucial 
    ingredients: i) a conformal change of variable, based on the approximate 
    location of the singularities of these functions in the complex U-plane;
    ii) a Bayesian inference technique, that takes into account additional
    known non-perturbative relations, in order to control the amplification of noise 
    occurring at large U. This general methodology is applied to the strongly
    correlated Anderson quantum impurity model, and is thoroughly tested both in- and 
    out-of-equilibrium. In the situation of a finite voltage bias, our method is able 
    to extend previous studies, by bridging with the regime of unitary conductance, and
    by dealing with energy offsets from particle-hole symmetry. We also confirm the 
    existence of a voltage splitting of the impurity density of states, and find that 
    it is tied to a non-trivial behavior of the non-equilibrium distribution function.
    Beyond impurity problems, our approach could be directly applied to Hubbard-like 
    models, as well as other types of expansions.
    Reconstructing non-equilibrium regimes of quantum many-body systems from
    the analytical structure of perturbative expansions
        Xavier Waintal
        
    ============================================================================================================================
    
    
    
    § INTRODUCTION
    
    
    The study of the out-of-equilibrium regime of strongly correlated many-body quantum problems is a major challenge in theoretical condensed matter physics.
    Its interest has grown rapidly in the past few years with new
    experiments, 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.>, proposed e.g. to build artificial neurons<cit.>.  
    These experiments raise the question whether the combination of strong correlation effects and out of equilibrium
    regimes could lead to genuinely new physics and phases of matter that do not have an
    equilibrium counterpart. 
    Quantum nanoelectronics also provide many examples of such systems.
    A classic example is the spin-1/2 Kondo effect occurring in a quantum dot, but recent experiments have 
    also managed to study in great detail underscreened <cit.> and
    overscreened <cit.> (multi-channel) Kondo effects, characterized by non-Fermi liquid 
    fixed points.
    Other notable examples of new quantum states induced by interactions are
    Luttinger liquids<cit.> that take place 
    at edges in the fractional quantum Hall regime, or the “0.7
    anomaly"<cit.> occurring in a simple
    quantum point contact geometry.
    Last, solid state based quantum computers such as spin qubits devices are nothing but out-of-equilibrium quantum 
    many-body systems (few sites Hubbard like models, possibly connected to electrodes) that bring new questions into 
    the scope of correlated systems<cit.>.   
    
    It is worth noting that even the simplest of these out-of-equilibrium problems,
    the single impurity Anderson model, is still the subject of active
    research<cit.>. 
    Early approaches used a range
    of approximate techniques including 4th order perturbation
    theory<cit.>, equation of motion techniques<cit.> and
    the Non Crossing Approximation (NCA)<cit.>.  
    State of the art techniques include the time-dependent Numerical Renormalization Group (NRG) and the
    density matrix renormalization group (DMRG)<cit.>. 
    
    Early attempts of real time quantum Monte Carlo <cit.>
    have 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 Schwinger-Keldysh diagrammatic Quantum Monte Carlo <cit.> (QMC).
    The later, which we use in this paper, reaches the infinite time steady state limit 
    and has a complexity which does not grow with time.
    The development of controlled computational methods is critical for the development of the theory
    in this field. Beyond its direct application to impurities and quantum dot
    physics, the Anderson model is of direct interest for quantum embedding methods
    such as Dynamical Mean Field Theory <cit.> (DMFT)
    which reduce bulk lattice problem to the solution of a self-consistent quantum
    impurity model.
    
    A straightforward approach to study out-of-equilibrium many-body quantum problem
    is to compute the systematic perturbative expansion of some physical quantity
    F in power of the electron-electron interaction U: F(U) ≡∑_n=0^∞ F_n U^n. In practice, F may depend on time (or frequency)
    as well as voltage-bias, temperature, etc.
    The coefficients F_n are given by the out-of-equilibrium Schwinger-Keldysh version of the Feynman diagrams<cit.>.
    Such a perturbative expansion is a central tool in quantum mechanics and quantum field theory.  
    In weak coupling theories, a few orders are sufficient to explain many physical phenomena, even
    quantitatively, as e.g. in Quantum Electrodynamics (QED). 
    However, at intermediate or strong coupling, this approach faces two main challenges: 
     (i) the computation of the coefficients for n large enough and 
     (ii) the reconstruction of the physical quantities as a function of U from a finite number of coefficients.
    
    Using the standard Wick theorem, an explicit expression of F_n to order n
    can be written as n-dimensional integrals.
    While the computation of F_n can hardly been achieved analytically beyond a few orders, 
    Quantum Monte-Carlo (QMC) algorithms known as “diagrammatic
    Monte-Carlo"<cit.> are able to 
    compute a finite number of these coefficients F_n for a general class of
    quantum many-body problems, in practice up to 8 or 15 depending on the model and the physical quantity.
    The first generation of these algorithms
    explicitly sampled the Feynman diagrams one by one with a complex Markov chain,
    moving from one diagram to another.  A second generation of algorithms handles the diagrams
    collectively using combinations of determinants to cancel disconnected diagrams
    in physical quantities.
    This was achieved in the real time Schwinger-Keldysh formalism<cit.>,
    and in the imaginary time Matsubara formalism<cit.>.
    
    The resummation of the series is a non-trivial mathematical task outside of the weak
    coupling regime, even with a perfect knowledge of the coefficients F_n.  
    The issue comes from the finite radius of convergence of the series.  When U is
    larger than this radius, the truncated series to the first N-th terms does not converge with N
    and some resummation technique must be used to compute F(U).
    Moreover, there are two additional difficulties associated with numerical methods: 
     i) only a finite number of coefficients F_n can be computed since
    the computation cost is exponential in n and  ii) the F_n are
    only known with a finite precision, typically of a few digits in QMC.
    
    
    
    In this paper, we approach this problem from the angle of complex analysis.
    Indeed, the divergence of the series originates from the singularity structure 
    of the function F(U) in the complex plane U (lower left panel in
    Fig. <ref>). We discuss how to locate the singularities closest to 0,
    and how to construct an analytic change of variable to resum the series beyond weak coupling (lower 
    right panel in Fig. <ref>). We also introduce a Bayesian technique to take into account the 
    amplification of the Monte-Carlo noise in the resummation process using some simple non-perturbative 
    additional information on the model.
    
    While our approach is quite general, we will focus here on the non-equilibrium Anderson quantum impurity model
    in the quantum dot configuration (upper panel in Fig. <ref>).
    Our starting point is an expansion of the Green's function in power of the
    Hubbard interaction U, using an extension of the algorithm of Ref.<cit.>.
    The algorithm is discussed in details in a companion paper <cit.>, 
    its implementation is based on the TRIQS library<cit.>.
    This algorithm provides a numerically
    exact computation of the perturbative series of
    physical quantities in power of the interaction U, at a cost which is  uniform in time but exponential with the expansion order.
    Hence it allows to compute in a transient regime as well as directly in a long time steady state, a regime in which most competing methods have severe limitations. 
    
    This paper is organized as follows.
    Section <ref> introduces our notations for the single impurity
    Anderson model. Section <ref> develops the resummation technique
    and illustrates it on the Kondo temperature. Section <ref> performs a
    benchmark of the method against NRG for the equilibrium dynamics. Section <ref> 
    presents new results in the non-equilibrium regime, including the voltage-split spectral 
    function, extended-range current-voltage characteristics, and a non-trivial dot distribution 
    function.
    Section <ref> concludes this article and presents perspectives for
    our conformal approach to the perturbative expansions of strongly interacting
    quantum systems.
    
    
    
    
    
    § THE ANDERSON IMPURITY MODEL
    
    
    
    In this paper, we focus on the single impurity Anderson model both at
    and out-of equilibrium.  While originally formulated to describe the effect of
    magnetic impurities in metals, this model is widely used in theoretical
    condensed matter, both as a simple model for quantum dots in mesoscopic
    physics and as a building block of “quantum embedding"  approximations like
    DMFT and its generalizations.  At the core of the
    Anderson model lies Kondo physics. The repulsive interaction on the quantum dot leads
    to an effective antiferromagnetic interaction between the
    electronic reservoirs and the spin of the (unique) electron trapped in the
    quantum dot in the local moment regime. This interaction leads to the formation of the Kondo resonance, a
    thin peak in the local density of state pinned at the Fermi
    energy<cit.>. 
    The Anderson impurity Hamiltonian reads:
    
        H   =   ∑_i=-∞^+∞∑_σγ_i
        c^†_i,σc_i+1,σ+h.c.
        +ϵ_d (n_↑+n_↓)
        
               +Uθ(t)(n_↑-1/2)(n_↓-1/2).
    
    It connects an impurity on site 0 to two semi-infinite electrodes i<0 and
    i>0.  The model corresponds to a single level artificial atom as sketched in the upper panel of 
    Fig. <ref>. Here ϵ_d is the on-site energy
    of the impurity (relative to the particle-hole symmetric point), 
    n_σ = c^†_0,σc_0,
       σ is the impurity density of spin σ electrons.
    c^†_i,σ and c_i,σ are the creation and
    annihilation operators for electrons on site i with spin σ.
    We use ħ = e = 1.
    θ(t) is the Heaviside function: We switch the interaction on at time t=0.  
    Typical calculations will be performed for large
    times so that the system has relaxed to its stationary regime.  The hopping
    parameters are given by γ_i = 1 except for γ_0=γ_-1=
    γ which connect the impurity to the electrodes.  The calculations can be
    performed for arbitrary values of γ. However, since we are not
    interested in the large energy physics of the electrodes, we suppose that
    γ≪ 1, i.e. that the tunneling rate from the impurity to the
    electrodes is energy independent Γ = 2 πγ^2 ρ_F where
    ρ_F is the density of states of the electron reservoirs at the Fermi level.
    The non-interacting retarded Green's function of the free impurity
    is given by
    
        g^R(ω) = 1/ω - ϵ_d + iΓ.
    
     The two electrodes have a chemical potential symmetric with respect to zero ± V_ b/2 which 
     corresponds to a bias voltage V_ b. 
     They share the same temperature that we take very low T=10^-4Γ. Within the standard non-equilibrium Keldysh formalism <cit.>, the non-interacting
    lesser and upper Green's functions are given by:
    
        g^<(ω)    =   iΓ[ n_F ( ω+V_ b/2 ) + n_F (ω-V_ b/2) ]
        /(ω - ϵ_d)^2 + Γ^2,
        
        
        g^>(ω)    =   iΓ[  n_F (ω+V_ b/2 ) + n_F (ω-V_ b/2) - 2]
        /(ω - ϵ_d)^2 + Γ^2,
    
    where n_F(ω)=1/(e^ω/T + 1) is the Fermi function.
    g^>(ω) and g^<(ω) are the starting point for the expansion in power of U that will be performed with real-time diagrammatic quantum Monte-Carlo.
    
    
    The quantities of interest in this article are the interacting Green's functions (denoted with capital letters),
    
        G^R(t,t')    =    -iθ(t-t') <{c_0↑(t), c_0↑^†(t')}>, 
        
            G^<(t,t')    =    i <c^†_0↑(t') c_0↑(t)>, 
        
            G^>(t,t')    =    -i <c_0↑(t) c_0↑^†(t')>,
    
    where the operators have been
    written in Heisenberg representation. Since we will restrict ourselves to the
    stationary limit, these functions are a function of t-t' only and can be
    studied in the frequency domain. Of particular interest is the spectral
    function (or interacting local density of state) given by
    
    A(ω) = -1/π Im[G^R(ω)].
    
    The equilibrium spectral function displays the important features of Kondo physics: 
    a sharp Kondo resonance at the Fermi level, and satellite peaks around
    ω=± U/2 in the case of particle-hole symmetry.
    
    
    
    The out of equilibrium spectral function can be used for the computation of the current-voltage characteristic 
    using the Wingreen-Meir formula<cit.>,
    
    I = Γ/2∫  A(ω) [ n_F (ω+V_
    b/2) - n_F (ω-V_ b/2) ] dω.
    
    The retarded self energy Σ^R(ω) is defined from the interacting
    Green's function by:
    
        G^R(ω) = 1/ω - ϵ_d + iΓ - Σ^R(ω).
    
    Physical quantities have systematic expansion in power of U
    G^R(t - t') = ∑_n=0^+∞ G^R_n(t-t') U^n,
    
    from which we obtain the corresponding quantity in the frequency domain by Fourier transform,
    
    G^R(ω) = ∑_n=0^+∞ G^R_n(ω) U^n.
    
    We obtain the functions G^R_n(ω) (typically up to n=10) using the QMC algorithm
    of Ref. Bertrand2019a, Profumo2015. The expansion of the self-energy 
    Σ^R(ω) = ∑_n=0^+∞Σ^R_n(ω) U^n
    
    is obtained from the G^R_n(ω)
    using a formal series expansion order by order of the Dyson equation (<ref>). 
    As an illustration, Fig. <ref> shows the self-energy series, up to order 10, for the equilibrium particle-hole symmetric model as obtained from diagrammatic QMC<cit.>. 
    These series are the starting point of this paper, 
    which is devoted to the resummation of the perturbative expansion for the Green's function and the self-energy beyond weak coupling.
    
    
    
    
    
    § THE PERTURBATIVE SERIES BEYOND THE WEAK COUPLING REGIME
    
    
    Diagrammatic Quantum Monte Carlo yields the first orders of the perturbation
    expansion of physical quantities, with some error bars.
    In weak coupling, we can directly sum this series and obtain the physical quantities
    with a few orders. Beyond weak coupling however, we face a more complex problem.
    For a given physical quantity F, we want to evaluate F(U) from 
    the first N (typically N∼ 10) coefficients F_0, F_1,
    F_2…F_N of a series F(U) ≡∑_n=0^∞ F_n U^n.
    In the following, F will stand for the width of the Kondo peak, the Green's function G or the self-energy Σ of the impurity.
    In the latter cases, the coefficients are functions of the frequencies, G_n(ω) and Σ_n(ω).
    We also want to know, for a given physical quantity F and interaction U, how many orders N_0
    are needed to obtain F(U) at a given precision.
    Since the cost of the diagrammatic QMC approach is exponential in N_0, 
    the answer to this question gives the ultimate limit of the method.
    
    The mathematical problem of series resummation is a quite old topic, e.g. 
    Ref. Hardy1949. Various techniques have been used in physics problems including Padé approximants
    <cit.>, Lindelöf extrapolation<cit.> or Cesàro-Riesz technique<cit.>.
    In diagrammatic QMC, this is typically a post-processing step: the Monte-Carlo
    produces the values of the various orders of the expansion, and 
    one then attempts to sum the series to obtain the final result.
    However, the situation is quite different if we want to use
    such technique to solve quantum impurity models in the context of 
    the quantum embedding methods like DMFT<cit.>, 
    or e.g. Trilex<cit.>.
    Indeed, in such cases, the method require multiple solutions of impurity model
    to solve their self-consistency loop.
    Therefore, it is necessary to develop more robust methods to sum the perturbative 
    series for impurity systems, which could be automatized.
    
    In the cases considered in this paper (quantum
    impurity models), and in general for lattice models at finite temperature (such
    as the Hubbard model), the series for F is expected to have a non-zero radius of
    convergence R_F. Note that R_F not only depends on the chosen physical
    quantity F, but may also depend on frequency, voltage, temperature, etc.
    R_F separates the  weak coupling regime (|U|<R_F) from the  strong
    coupling regime (|U|>R_F).  At weak coupling, the truncated
    series ∑_n=0^N F_n U^n provides an accurate estimate of F(U) and is
    controlled exponentially with the number of coefficients N (like a geometric
    series since F_n∼ (1/R_F)^n).  At strong coupling however, this truncated
    series diverges.
    Note that in some problems like e.g. the unitary fermionic gas, the series has a 
    zero radius of convergence at zero temperature,
    see e.g. Ref. Rossi2018 for a recent example with diagrammatic QMC. 
    We will not consider these cases in this paper, as they require other techniques as the ones presented here, 
    e.g. Borel summation techniques. 
    
    In this paper, we consider the series summation problem with the angle
    of reconstructing the function F(U) in the complex U plane.
    The divergence of the series is due to the presence of singularities
    in the complex U plane, starting on the circle |U|=R_F.
    The question is to reconstruct F beyond the radius of convergence.
    
    
    
    
     §.§ General theory
    
    
      §.§.§ Conformal transformation
    
    
    Conformal transformations can be used to deform
    the complex plane and bring the point to be computed back into the convergence disk of a transformed series.
    This technique was used a long time ago e.g. in statistical physics <cit.>. 
    In a previous work<cit.>, some of us have shown that a simple conformal Euler transform allows 
    to compute the density on the impurity up to U= ∞, at very low temperature, from the first 12 coefficients of the series.
    However, this Euler transform is not always successful in resumming other quantities like the Green's function 
    and the self-energy, and needs to be generalized.
    
    Suppose that we aim at evaluating F(U) at U= U_0 with U_0 real, positive
    and U_0>R_F.
    First, we assume a  separation property,  i.e. that we can find a simply connected domain delimited by a curve C 
    containing 0 and U_0 but no singularities of the function F, as
    illustrated in the lower left panel of Fig. <ref>. The singularities of the function F(U) will be located outside the domain C.
    We then proceed as follows:
    
      *  First, according to the Riemann mapping theorem, we can construct a biholomorphic change of variable W(U) such that  i) W(0) =
    0,  ii)  it maps the interior of C into a disk D_ C centered at 0 in the W plane
    (see the lower right panel in Fig. <ref>).
    In practice, we seek C to separate the singularities from the half straight line of real positive U.
    In the following, we will use two simple transformations, but in general we could use a Schwarz–Christoffel map if C is a polygon <cit.>, 
    composed with a Möbius transformation of the disk to enforce  i) .
    
    
      *  Second, we form the series for the reciprocal function U(W) of W(U) which is defined term by term by the equation U(W(U))=U.
    We then construct the series F̅(W) ≡∑_p F̅_p W^p  defined by the composition
    F̅(W) = F(U(W)).
    Since W(0)= 0,
    the first N terms of F(W) can be computed from the first N terms of F(U).
    
    
      *  We evaluate the series F̅(W_0) at the point of interest W_0= W(U_0). 
    Indeed, by construction W_0∈ D_ C and, since F̅(W) is holomorphic in  D_ C, D_ C is included in the convergence disk 
    of the series F̅. Hence the series F̅ converges at W_0.
    
    
    The result is independent of the choice of the domain C but 
    the speed of convergence of the series for F̅(W_0) versus N is not, since it is determined by the relative position
    of W_0 compared to the radius of convergence R_F̅ of F̅, i.e.η_ C≡ |W_0/R_F̅|.
    Therefore, there are ways to optimize the domain C.
    For example, we can not simply take a narrow domain close to the real axis, for the convergence in W would
    be really slow: we need to have U_0 and 0 as “far" as possible from the curve C
    (the precise meaning of “far" being given by η_ C).
    For each domain C satisfying the separation property, there is a minimum number of orders 
    N_ C needed to obtain the result at a given precision ϵ.
    There is therefore an optimal domain, which minimize N_ C to N_opt = min_ C N_ C.
    This is the absolute minimum of orders needed to sum the series, and therefore 
    determine  in fine the complexity of the diagrammatic QMC algorithm.
    Our next goal will be to approach such optimum.
    
    
    Note that a failure of the separation assumption, i.e. the choice of a domain
    containing singularities, may result simply in the divergence of the series F̅ at W_0,
    hence a clear failure of the method rather than a wrong result.
    Conversely the study of the convergence radius of the F̅(W)
    series provides direct information on the singularity free regions of the U plane. Indeed, the region of the U plane that maps towards the inside of the convergence radius of F̅(W) are singularity/branch cut free. 
    Hence, using several conformal transforms, one may perform a step by step construction of the domain
    C. 
    Another note is that, as a consistency check, one can also check the stability of the final result upon small deformations of the domain (or the W(U) function), 
    as was discussed in details in Ref. Profumo2015 for the Euler transform.
    
    The existence of the domain C and the transformation W(U) has a
    direct consequence on the algorithmic complexity of diagrammatic Quantum Monte-Carlo.
    It was shown in Ref. Rossi2017a that, 
    for values of U inside the convergence radius, connected diagrammatic quantum Monte-Carlo techniques
    provide a systematic route for calculating
    the many-body quantum problem in a computational time that only increases 
    polynomially  with the requested precision.
    The result also applies to the Keldysh diagrammatic QMC.
    For completeness, the core of the argument is as follows: 
    inside the radius of convergence R, the precision of a calculation
    ϵ increases exponentially with the number of orders N used ϵ∼ (U/R)^N. Hence, although the computational time C increases
    exponentially with N,  C∼ a^N, 
    the overall computational time scales as C ∼ (1/ϵ)^log a/log (R/U), i.e.
    polynomially, see Ref. Rossi2017a for a detailed analysis.
    For a given U_0 and domain C, we now have to sum the transformed
    series F̅ inside the radius of convergence. Hence the same argument also apply for this series, 
    and therefore we conclude that,  even outside the disk of convergence, 
    we expect the algorithm to have a polynomial complexity as a function of the precision.
    Let us emphasize however that this result is largely academic, since in practice the power law can be large.
    Moreover, as we will discuss, for some physical quantities
    the transformation to W can lead to a dramatic increase of the noise
    which induces a large computation time for a given precision.
    
    
    
    
    
      §.§.§ Location of singularities in the complex U plane
    
    
    In order to choose C properly, we need to have some information on the location
    of the singularities in the U plane. In this paper, we use the following technique to approximately
    locate the poles of F(U) in the complex plane. 
    
    
      *  We form an inverse of F of the form  K(U)=1/(F(U)+a) as a formal series (i.e. order by order).
    	 a is a constant that we choose at our convenience. In order for the series K(U) to exist, we must have F_0 + a ≠ 0.
       
      *  We estimate the radii of convergence R_F (resp. R_K) of F (resp. K), 
             by plotting |F_n| and |K_n| versus n, and fitting the asymptote |F_n|∼ (1/R_F)^n
      *  In most of situations we found R_F ≠ R_K. If not, we used a different a so as to obtain R_F ≠ R_K. Without loss of generality, let us assume that R_K is the largest.
             We use the truncated polynomial of the series,  ∑_p=0^N K_p U^p
    	 to compute K(U) within its disk of convergence and therefore locate its zeros, which are the poles of F.
             They will appear as the accumulation of the zeros of the polynomials at large enough N.
    	 If R_F > R_K, we simply reverse the roles of the series and reconstruct K(U).
    
    
    This technique has a quite large degree of generality, but also limitations.  It
    assumes for example that the leading singularities in F are poles and that the
    radius of convergence of F and K are different. 
    Also it does not give us indications of poles that would be far from the origin but close to the real axis.
    However, in practice, we will see
    below that for the quantities and the physical problem considered in this paper (Green's function and self-energy in
    real frequency, and Kondo temperature), this technique is sufficient.
    Finally, once F(U) has been re-summed, it can be used to locate its zeros, hence for the resummation of K(U)
    which provides another consistency check of the method.
    
    
    
    
      §.§.§ Controlling the noise amplification using non-perturbative information and Bayesian inference
    
    
    The transformation from F_n to F̅_p is a linear one (with a lower triangular matrix), for a given transformation W(U).
    Depending on the eigenvalues of the corresponding matrix, the Monte-Carlo error bar in F_n may be strongly amplified by the transformation.
    As a result, the method may become unusable at strong coupling, as will be
    illustrated below on Fig. <ref>.
    
    However, if we add some  non-perturbative information, such as the fact that the Kondo
    temperature vanishes at infinite U, or a sum rule, we can
    construct a Bayesian inference technique that may be used
    to decrease the statistical uncertainty.
    Bayesian inference provides a systematic and rigorous way to incorporate this information into
    the results and improve their accuracy.
    In the rest of this paragraph, we describe the general theory for this technique.
    We will illustrate it in the following section.
    
    Let us consider a series F(U)=∑_n=0^N F_n U^n where the F_n are known with a finite precision.
    We note F = { F_0, F_1, … F_N} the corresponding (vectorial) random
    variable.
    We calculate the mean values ⟨ F_n⟩ and the corresponding errors δ_n within the
    quantum Monte-Carlo technique. We assume that the coefficients F_n
    are given by independent Gaussian variables. This forms the “prior"
    distribution P_ prior(F=f) in the absence of additional information. 
    
        P_ prior(F=f) = ∏_i=0^N 1/√(2πδ_n) e^-(f_n-⟨ F_n⟩)^2/2δ_n^2
    
    Let us note the additional information X. X is a random variable that can be
    directly calculated from the series, X=g(F) but whose actual value is also
    known very precisely by other means. In the example below, X will be the
    value of F(U) at large U. Bayesian inference amounts to replacing the prior
    distribution with the posterior distribution P(F=f| X=x_0) that incorporates
    the knowledge of the actual value of X (we note P(A|B) the conditional
    probability of event A knowing event B). The value of X is often known
    exactly. However, due to the presence of truncation errors, its value cannot be enforced exactly,
    and we suppose that it is known
    with a small error ε. Eventually, we take the limit
    ε→ 0. Hence, we assign to X a Gaussian probability
    distribution
    P_X(X=x) = 1/(ε√(2π)) exp[-(x-x_0)^2/(2ε^2)] and define the posterior distribution as,
    
        P_ posterior(F=f) ≡∫ dx P(F=f| X=x) P_X(X=x).
    
    Using Bayes formula P(F=f| X=x) = P(X=x| F=f)P_ prior(F=f)/P_ prior(X=x) and the deterministic relation
    P(X=x| F=f)=δ[x-g(f)], one arrives at,
    
        P_ posterior(F=f) = P_X(X = g(f))  P_ prior(F=f)/P_
        prior(X=g(f)).
    
    In practice, one proceeds as follows: (i) one generates many series according
    to P_ prior(F=f). We emphasize that these series result from a single
    QMC run, hence are trivially generated (independent Gaussian numbers). Bayesian
    inference implies no significant computational overhead (ii) One construct a
    histogram of the values of X
    to obtain P_ prior(X=g(f)). (iii) Each series is given a weight P_X(X =
    g(f))/P_ prior(X=g(f)) which is used to calculate other observables such as
    the value of F(U) at different values of U. In
    practice the results are insensitive to the choice of ε as long as
    it is chosen large enough so that a finite fraction of the sample contributes
    to the final statistics.  
    
    
    
    
    
    
     §.§ Illustration with the Kondo temperature
    
    
    Let us first apply the method described above to the Kondo temperature T_K (which will be F in this section).
    T_K corresponds roughly to the width of the low energy Kondo peak,
    and is defined more specifically in this paper as the dimensionful Fermi liquid quasi-particle weight extracted from the 
    retarded self-energy at low energy:
    
        T_K(U) ≡2Γ/1 - .∂_ω ReΣ^R(U, ω)|_ω=0.
    
    Our first goal is to illustrate how the method actually works, and benchmark it against
    the calculation of the same quantity from the Numerical Renormalization Group (NRG) technique
    and Bethe ansatz<cit.>.
    
    
    
      §.§.§ Singularities in the complex U plane
    
    
    The dashed blue lines of Fig. <ref> shows the truncated series of
    T_K = ∑_n=0^N F_n U^n for various orders N≤ 10.
    These truncated series diverge around R_T_K≈ 5Γ which is the convergence
    radius of the series for these parameters. Increasing the value of N helps to obtain a reliable
    value of T_K closer to R_T_K. However, as expected, even with a very large number of terms,
    the bare series cannot be summed near or above R_T_K. Anticipating the final results, the plain red 
    line corresponds to the results after resummation which matches very well what was obtained with our 
    benchmark NRG calculation (see Sec. <ref> for details on the used NRG
    implementation).
    
    The inset of Fig. <ref> shows the value of |F_n U^n| (blue
    circles) as a function of n for U/Γ=9 which lies above the
    convergence radius of the series. The log-linear plot shows an exponential
    increase of |F_n U^n| ∼ (U/R_T_K)^n with n which we use to extract the
    convergence radius of the series. Note that for other series, it can happen
    that |F_n| oscillates with n. Whenever F_n changes sign, it becomes
    close to zero which provides deviations from the clear exponential behaviour
    shown in the inset of Fig. <ref>. 
    Hence, to obtain convergence radii which are robust to these outliers, we used a robust
    regression method on the log|F_n| versus n data (we compute the regression
    slope as the median of all slopes between pairs of data points, this is known
    in statistics as the Theil-Sen estimator<cit.>).
    
    
    We now compute the first 10 terms of the series of 1/T_K(U). This series
    has a radius of convergence of the order of 10Γ. 
    We look for the zeros, in the complex plane, of the series 1/T_K(U) truncated at order N.
    Since the truncated series is a polynomial, it has
    (generically) N zeros, which are shown in Fig. <ref> for N=6
    (red squares), N=8 (blue circles) and N=10 (stars). One pair of zeros U
    ≈± i 5Γ is converged for all the truncations, hence corresponds
    to a true zero of 1/T_K(U), i.e. to a pole of T_K(U).
    Fig. <ref> also shows the circle |U|=R_T_K extracted from the
    analysis of the T_K(U) series done in the inset of
    Fig. <ref>. We find that the two poles ± i 5Γ do
    indeed lie right on this circle. 
    
    
    
      §.§.§ Conformal transformation
    
    
    Let us now turn to the conformal transformation W(U), which maps the
    two poles ± i5Γ  away and brings the values of interest U > 0 (real)
    closer to zero. 
    
    We illustrate the technique with two maps:
    the Euler map defined by 
    
        W = U/U-p,
    
    and the “parabola" map which is defined as
    
        W = -  tan^2 ( π/2√(U/p)),
    
    where p is an adjustable complex parameter.
    
    Fig. <ref> shows the various regions (different colors) in
    the U plane that are mapped onto concentric circles of the W plane.  0 is
    mapped onto 0 and p onto ∞ in both transforms. The Euler map (left
    column) maps one half of the plane into the unit disk and the other half into
    the outside of the unit disk (separated by a black line).  The parabola 
    transform (right column) maps the inside of a parabola (black line) into the
    unit disk and the outside of the parabola into the outside of the unit disk.
    In the case where there are no singularities on the positive half plane Re[U]>0,
    the Euler transform should be preferred since real values of
    U>0 are typically mapped closer to U=0 than with the parabola transform
    (compare the size of the blue region of the parabola and Euler case for
    instance).  However, the parabola map is more agnostic about the positions of
    the singularities and will work even if there are singularities on the positive
    half plane Re[U]>0 as long as they lie outside the parabola.
    
    We now perform the resummation of T_K(U). The series contains only 
    even power of U due to particle-hole symmetry, so that it can be considered as
    a function of U^2. The two poles U=± i5Γ correspond to a
    single one U^2 = -25Γ^2. In the U^2 plane, the pole being on the negative real axis,
    the Euler maps works very effectively. The resummation can also be performed with the parabola transform.
    
    
    Once the conformal map is selected, we form the series F̅_p in the W variable, as explained above.
    The inset of Fig. <ref> shows F̅_n W_0^n (red squares) as a function of n for
    W_0=0.7=W(U_0=9Γ), using the Euler map with p=-35Γ^2 (the parabola yields similar results with p=-15Γ^2).
    As expected, U_0 is way beyond the radius of convergence in the original variable U, 
    while W_0 lies within the disk of convergence of F̅(W) whose radius is found to be R_F̅≈ 2.
    The final result T_K(U) using the Euler transforms is shown in
    Fig. <ref>. The parabola transform (not shown) is undistinguishable from
    the Euler at this scale.
    
    In this work, singularities were never found near the real positive axis, so
    that all U > 0 can be reached using the conformal transforms of
    Fig. <ref>, given that enough orders of the series are
    known. However, one may very well build a conformal transform to reach a regime
    beyond a singularity by considering a concave contour C, as it is shown in Appendix <ref>. This may become
    interesting if a phase transition occurs when interaction is increased.
    
    
    
      §.§.§ Noise reduction with Bayesian inference
    
    
    Let us now apply the Bayesian inference technique described above to the
    computation of T_K(U).
    In the left panel of Fig. <ref> we have
    re-sampled the series for the Kondo temperature, i.e. we have generated many
    series (typically 10^3 to 10^5 samples).
    For each sample we perform the
    conformal transformation and plot the result for the Kondo temperature as a
    function of U (thin red lines). While we find that all results agree for U
    ≤ 6Γ, the bundle of curves start to diverge for larger values of U. In the
    middle panel, we plot (black thin line) the corresponding histogram of the
    values obtained for T_K(U=∞), which is P_ prior(T_K=g(f)). 
    
    We use the non-perturbative relation lim_U→∞ T_K(U) = 0.
    Hence we want to “post-select" the configuration of F_n which 
    give a vanishing Kondo temperature at large U, at precision ϵ.
    Following the  procedure described in Sec. <ref>, our final
    result is obtained by averaging the different traces (thin red lines) with the
    weight given by Eq. (<ref>). The right panel of
    Fig. <ref> shows the result for three different
    values of U as a function of ε which confirms that the results
    are insensitive to the actual value of ε. We find a very good
    agreement with the results obtained from NRG even at large values of U, 
    noting that NRG spectra have typical relative error bars of a few percents 
    (see Sec. <ref> for details).
    
    
    
    
      §.§.§ Benchmark with the Bethe Ansatz exact solution
    
    
    The series expansion for 1/T_K(U) has been calculated explicitly and exactly
    using the Bethe Ansatz technique by Horvatic and Zlatic<cit.>.
    Ref.Horvatic1985 provides an iterative formula for calculating the
    coefficients of the expansions and shows that the corresponding series has an
    infinite radius of convergence.  This provide another independent benchmark of
    the calculation of T_K(U) as well as of the method itself. We checked that
    the 10 first coefficients of this series agree with the one that we computed
    with QMC.
    
    Fig. <ref> shows our final result together with the NRG result
    (black circles) and the Bethe ansatz results. At this scale, the agreement is
    perfect. Using the exact series for 1/T_K(U) (truncated to around 50
    coefficients), we studied its zeros which are the poles of T_K(U). We find
    that they are situated on the imaginary axis. The poles closest to the origin
    are U/Γ≈± 4.89059579i in agreement with our findings, see Fig.
    <ref>. The next poles are U/Γ≈± 13.79i, 21.77i,
    29.89i, 37.87i and 45.9i but are too far to be accessible with only the
    first ten coefficients. 
    
    
    
    The right panel of Fig.<ref> provides a detailed benchmark of our
    results versus both NRG and the exact Bethe Ansatz solution.
    
    
    
    
    
    
    
    We find that the QMC results for T_K are slightly more accurate than NRG, because
    the extraction of T_K from the NRG self-energy (see Eq. (<ref>)) contains
    inherent broadening errors. The agreement between all three methods is nevertheless
    excellent. In addition, we can extract from the Bethe Ansatz the exact QMC
    error, and this error matches the measured 1 sigma statistical error bars.
    
    
    
    
     §.§ Equilibrium dynamical correlation functions
    
    
    Let us now apply our method to the 
    Green's function and self-energy as a function of the real frequency ω.
    
    
    
      §.§.§ Singularities in the long time (stationary) limit
    
    
    Let us now turn to the full Green's function G^R(ω,U) and self-energy
    Σ^R(ω,U). 
    An example of our bare data is shown in Fig. <ref> where we plot the coefficients
    Σ_n^R(ω) obtained from real time diagrammatic quantum Monte-Carlo for n=2,4,6,8 and 10.
    The description of the method used to calculate these coefficients Σ_n^R(ω) is explained in the companion paper to this
    article<cit.>.
    
    We focus on the quantity
    Σ^R(ω)-iΓ and denote its inverse
    F_ω(U)=1/(Σ^R(ω)-iΓ). The retarded Green's function can
    be recovered from F_ω(U) using G^R(ω) = 1/(ω
    -F_ω(U)^-1) (using ω -Σ^R(ω) +iΓ turns out to be less 
    convenient especially at high frequency).
    
    
    Fig. <ref> shows the convergence radius of F_ω(U) as
    a function of frequency, extracted from a study of the exponential decay of the
    corresponding series with n. We have also performed a systematic study of the
    zeros of Σ^R(ω)-iΓ in order to localize the poles of
    F_ω(U). We find one pair of poles at each frequency. The results are
    shown in the inset of Fig. <ref> for a set of frequencies
    from ω=0 to ω=10Γ in the complex plane for U^2.  The
    absolute value of the poles of F_ω(U) is also plotted in the main frame
    of Fig. <ref> as a function of frequency (circles of varying
    colors from blue to red). We observe a perfect match with our estimation of the
    convergence radius reflecting the fact that these poles are responsible for the
    divergence of the series. It is important to note here that working in the
    real frequency domain is very helpful: we found a single pole per
    frequency (at least for the range of interactions that we could study). 
    Hence, we expect that performing the resummation in real time or imaginary frequencies
    could be more complex, since all these poles would be involved simultaneously.
    
    
    
    
    The results for three frequencies (ω/Γ = 1, 2 and 6) are given in 
    Fig. <ref>. We show the convergence of the imaginary
    part of the self-energy using two different resummed series: F_ω(U)
    (green symbols) and 1/F_ω(U) (purple symbols).
    The former has been resummed with an Euler transform with a frequency
    dependant p set close to the poles shown in Fig. <ref>.
    The latter, for which our method did not detect poles, has been resummed with
    the parabola transform (in the U plane) with p=-4.5Γ. Again, Bayesian
    inference has been used to enforce lim_U→∞ G(U, ω) = 0 for all ω≠ 0.
    For comparison, we also include the NRG results (which are
    very accurate at small frequency and possibly less accurate at large frequency). 
     The slight difference between the purple and
    green curves is due to the truncation error.  We find
    that the series which has (initially) the largest convergence radius is less
    sensitive to truncation error or statistical noise than the other. 
    We attribute the small discrepancy
    between the QMC results and NRG at large frequency to a lack of convergence of
    the latter. These results are obtained for a rather strong interaction
    U=9Γ. At smaller interaction the QMC and NRG results become
    undistinguishable. At larger interactions, the QMC results become increasingly
    inaccurate due to truncation errors.   
    
    
    
      §.§.§ The long time limit
    
    
    In the Keldysh formalism, the interactions are switched on at an initial time (0),
    and one follows the evolution of the system with time t.
    We assume here that the system relaxes to a steady state at long time.
    Let us consider the average of an operator Ô as a function of time, and 
    its expansion ⟨Ô(t) ⟩ = ∑_n O_n(t)U^n (the extension of the following arguments 
    to Green's function is straightforward).
    
    At finite time t, the radius of convergence of this series is infinite, as shown in 
    Appendix <ref>.
    Each order in the perturbation expansion O_n(t) relaxes with t to a long time limit, 
    but the time t_relax(n) it takes to reach this limit can increase with n.
    The long time and large n limit do not commute in general:
    
        lim_n→∞lim_t→∞O_n(t) ≠lim_t→∞lim_n→∞ O_n(t).
    
    This behaviour was already noted in Fig. 14 of Ref. Profumo2015.
    It is also illustrated on Fig. <ref>, which shows 
    various orders n of the expansion of the current through the dot versus n, 
    for different times.
    We observe that at small times the orders I_n decreases faster than exponentially with n, 
    consistent with the bound mentioned above. The coefficients converge
    to the steady state limit at long time.
    
    
    At finite time t, since the series converges, it is sufficient to have enough orders.
    In the steady state, as explained above, we have a minimal order N_0 needed
    to compute the quantity at a given precision. One should then simply compute at a time
    t > t_relax(N_0).
    
    In the Anderson model, some quantities like the spectral function are known to relax
    on a long time scale t_K ∼ T_K^-1, see e.g. Ref. Nordlander1999.
    The previous remarks explain how the algorithm deals with this long time.
    For a given U, we need N_0(U) orders, hence to compute at a time larger than t_relax(N_0(U)).
    The larger U is, the longer this time becomes. However, it is still finite at fixed U, 
    and since our calculation of the perturbative expansion is  uniform in time, 
    it is not an issue (the computation effort does not grow with time).
    However, the existence of the Kondo time indicates that the number of orders necessary to compute e.g. the low frequency spectral function
    at a given U increases with U (otherwise the relaxation time of the physical quantity would be bounded at large U).
    
    
    
    
    
    § BENCHMARK OF THE DYNAMICS IN EQUILIBRIUM
    
    
    We now benchmark our results in the case of equilibrium, testing 
    various regimes of the Anderson impurity model.
    Let us first describe the high-precision NRG computations that were performed.
    
    
    
     §.§ NRG implementation
    
    
    The Numerical Renormalization Group (NRG) <cit.>
    was used to benchmark our QMC calculations in equilibrium,
    and to test the reliability of the series extrapolation method
    for spectral functions at various values of U and ϵ_d.
    In order to obtain precise NRG data for the spectral function 
    of the Anderson impurity model, the computations were performed 
    using several improvements over the simplest implementations of
    the NRG. First, the full density matrix formulation of 
    NRG <cit.> was used to reduce finite size effects
    due to the NRG truncation. Second, symmetries of the problem were heavily 
    exploited <cit.>, allowing to reduce significantly the Hilbert 
    space dimension of various multiplets. In the particle-hole symmetric case, 
    the full SU(2)_charge⊗SU(2)_spin
    symmetry was used, while the charge sector was reduced to U(1)_charge 
    away from particle-hole symmetry. Third, the impurity Green's function
    was extracted from a direct computation of the d-level
    self-energy Σ(ω) <cit.>, according to its exact 
    representation as the ratio of two retarded correlation functions in the
    frequency domain:
    
        Σ(ω) = U F^R(ω)/G^R(ω),
    
    where G^R(t)=-iθ(t)⟨{d^†_σ(0),d_σ^†(t)}⟩ is the
    usual single particle retarded Green's function in the time domain, and
    F^R(t)=-iθ(t)⟨{d^†_σ(0) 
    d_-σ^†(0) d^†_-σ(0),
    d_σ^†(t) }⟩ is a composite fermionic correlation function.
    In practice, Im[G^R(ω)] and Im[F^R(ω)] are
    computed from the Källén-Lehmann representation using the broadened NRG 
    spectra, and the real parts of both G^R(ω) and F^R(ω) are obtained 
    via a Kramers-Kronig relation.
    Finally, the truncation parameters of the NRG simulations were taken 
    to model as closely as possible a continuous density of states for
    the electronic bath. Although the use of the logarithmic Wilson
    discretization grid, ω_n = D Λ^-n, is inherent to 
    the practical success of NRG, we found that values of Λ
    as low Λ=1.4 could be managed in practice within the NRG, 
    taking a very large number N_kept=3200 of kept multiplets.
    Up to N_iter=120 NRG iterations were used, so that the
    effective temperature can be considered to be practically zero.
    With such small value of Λ, the broadening parameter b 
    of the discrete NRG spectra could be decreased down to b=0.2, without
    z-averaging, which further enhanced the spectral resolution of 
    the Hubbard satellites in the spectral function.
    
    
    
    
    
    
     §.§ Comparison to NRG in equilibrium
    
    
    Fig. <ref> shows the spectral function as well as the imaginary and real part of the
    self energy for the symmetric Anderson impurity in the strong correlation
    regime U=9Γ (same data as the purple curve of Fig. <ref>). 
    The spectral function shows a clear Kondo peak and the two satellites at
    ω≃± 4.5Γ=± U/2 in good agreement with the NRG data. For this calculation, 
    a simple second order calculation of the self-energy already provides a reasonably good result 
    (thin black line), due to near cancellations in higher order diagrams in
    the peculiar case of particle-hole symmetry. 
    
    Fig. <ref> shows the same plot in the asymmetric case
    ϵ_d = 1. This case is more complex because the resonance at U=0 is
    offset with respect to the Fermi level, hence to the position of the Kondo peak. We note
    that previous real time QMC techniques suffered from a strong sign problem and
    could not access the asymmetric regime <cit.>.  We also stress that the
    second order approximation is now very different from the correct result.
    The comparison to the NRG data is still excellent.
    
    
    
    Another advantage of the techniques described in this article and its
    companion article <cit.> is that a single QMC run provides the full
    dependence in both ω and U, which is very time consuming in the NRG. 
    This is illustrated in
    Fig. <ref> where the color map shows the spectral function as a
    function of ω and U. One can clearly observe the formation of the
    Kondo peak (which gets thinner as one increases U and shifts toward
    ω=0 in the asymmetric case) as well as the Hubbard bands at ω=±
    U/2. Note that the results are perfectly well behaved (qualitatively correct)
    up to very large U (even above U=12Γ shown in
    the plot) but  become quantitatively inaccurate at too large values of U.
    Improving them would require the use of higher perturbation orders.
    
    
    
    
    
    § OUT OF EQUILIBRIUM RESULTS
    
    
    We finally turn to the out-of-equilibrium regime, and present some 
    accurate computation of current-voltage characteristics, as well as 
    novel predictions for dynamical observables in presence of a finite
    bias voltage.
    
    
    
     §.§ Splitting of the spectral function
    
    
    Fig. <ref> shows
    the spectral function of the symmetric impurity in the presence of various bias
    voltages from V_ b=0 to 4Γ. The results were obtained using
    the parabolic map on the series of Σ_ω(U^2) - iΓ (with an
    optimized frequency dependent parameter  p/Γ^2∈[-25, -200]). Upon
    increasing the bias voltage, we find as expected from NCA<cit.>
    and perturbative<cit.> calculations that the Kondo resonance
    simultaneously broadens and get split into two peaks.  Previous results on
    the spectral function<cit.> were based on the bold diagrammatic
    approach and were calculated at relatively high temperature (T = Γ/3)
    while using a third terminal for computing the spectral function.
    
    
    
    Most of the results of this paper have been obtained at very low temperature.
    We emphasize however that increasing the temperature makes the calculations easier:
    indeed at finite temperature, the non-interacting Green's functions decrease exponentially 
    as e^-t/T instead of the algebric decay at zero temperature. It follows that the support of the integrals to be calculated is smaller, hence the convergence of the calculation faster.
    We show a calculation at finite temperature in Fig. <ref> where
    we have computed the spectral density of the symmetric impurity at temperature T = Γ / 50
    under a bias voltage V_ b = 0.6Γ and V_ b = 1.5Γ. A
    single Monte-Carlo run allows us to observe the splitting of
    the Kondo resonance as U is increased (upper panel). The result is
    quantitatively accurate up to U ≈ 8Γ (lower panel) but remains
    qualitatively meaningful at higher interaction (upper panel).
    
    The fate of the Kondo resonance out-of-equilibrium,
    in presence of a bias voltage, can be understood qualitatively from the
    interplay of two phenomena. On the one hand, the bias voltage induces a
    splitting of the Fermi energies of the two reservoirs, hence one expects a
    corresponding splitting of the Kondo resonance. On the other hand, 
    the voltage, like the temperature, increases the energy and phase space 
    for the spin fluctuations, leading eventually to the disappearance of the
    Kondo resonance<cit.>. The
    competition between both effects leads to the appearance of the splitting
    only above a finite voltage threshold (about V_b≃Γ in the plot 
    of Fig. <ref>).
    
    
    
    
     §.§ I-V transport characteristics
    
    
    Fig. <ref> shows the results obtained for the I-V characteristics in
    the symmetric case ϵ_d=0. The resummation has been done for the series of 1/I(U^2)
    using a parabolic transform with p=-40Γ^2. At small bias, we recover a
    perfect transmission I=(e^2/h) V_ b due to the unitary Kondo
    resonance, while for eV_ b> k_BT_K the conductance experiences an extra
    suppression by the interaction (Coulomb blockade). We find a very good match with a previous calculation
    from Ref. Werner2010.  The present technique
    allows one to lift the main limitations that Ref. Werner2009,Werner2010
    was facing: we can now access long times (here we have used ∼ 20/Γ but
    it could be increased further if necessary) to be compared with maximum times
    of the order of ∼ 3-5/Γ in Ref. Werner2010. As a consequence,
    we can reach the low bias regime, which was not accessible in
    Ref. Werner2010. Another important point is that the method is not
    limited to the symmetry point as we now demonstrate.
    
    Fig. <ref> shows the I-V characteristics 
    for an  asymmetric model with ϵ_d/Γ=1. 
    The results have been obtained from the
    resummation of 1/I(U) with a parabolic transform (p=-6Γ) and no
    Bayesian inference. The I-V characteristics is particularly interesting
    because, due to the asymmetry, the non-interacting low bias transmission is
    modified by interactions and one must first build up the Kondo resonance to approach 
    I≃ (e^2/h) V_ b
    (note that the unitary limit is strictly exact only at ϵ_d=0,
    and the conductance is slightly lower than e^2/h otherwise in the Kondo
    regime).
    This behavior leads to a non monotonous current versus U: as one
    increases U, the current first increases until the Kondo resonance is fully built
     (see the bottom panel of Fig. <ref>).
    As one increases further U, the Kondo width T_K shrinks and
    the current decreases as Coulomb blockade starts to set in.
    
    
    
    
     §.§ Biased distribution function
    
    
    Finally, we discuss the out-of-equilibrium distribution function of the
    impurity, i.e. its energy-dependent probability of occupation.  We define the
    distribution function n(ω) as 
    
    n(ω) = G^<(ω)/2π i A(ω),
     
    so that at equilibrium n(ω) is
    simply the Fermi function n_F(ω). 
    Without interaction, the distribution function amounts (at zero temperature) to a double step 
    function
    n(ω)_U=0 = [n_F(ω-V_ b/2)+n_F(ω+V_ b/2)]/2. 
    We want to investigate the behaviour of n(ω) as U increases, a
    question that was not addressed in previous literature to the best of our
    knowledge.
    
    The results are shown in Fig. <ref>. In this
    particular case, the series are fully alternated which means that the
    singularity lies on the negative real axis. We could sum the series using an
    Euler transform (p=-8Γ^2) up to U=+∞. 
    We find that the function n(ω) is not thermal,  i.e. it can not be fitted 
    by a Fermi function n_F with an effective temperature.
    In particular, it still exhibits discontinuities at the position of the lead Fermi surfaces, 
    which we expect to be rounded at finite temperature. Interestingly, these
    discontinuities are comparable to the equilibrium quasiparticle weight for
    U=4Γ, do not seem to vanish in the limit U=∞. Also very
    striking is the quasi-linear behavior of n(ω) that is observed for 
    -V_b/2<ω<V_b/2.
    
    Experiments that measure the non-equilibrium distribution function quantity typically use a 
    third (for instance superconducting) terminal weakly coupled to the system 
    <cit.>. 
    To the best of our knowledge, this quantity has not been measured in quantum
    dots, and we hope that the present prediction may stimulate some experimental activity.
    
    
    
    
    § CONCLUSION: THE FALL OF THE CONVERGENCE WALL
    
    
    We have presented a systematic computation of the perturbative expansion of the
    Anderson impurity model in and out of equilibrium in power of the interaction
    strength U. The main advantage of our Keldysh expansion approach is its
    ability to calculate directly in the long time steady state regime.
    Using our approach, we were able to obtain improved or novel results regarding the
    non-equilibrium dynamics of strongly interacting quantum dots.
    
    The main contribution of this article lies in the systematic construction of a
    set of conformal transformations that provide a practical route for a
    mathematically controlled resummation of series.  We have shown how to use
    analytic conformal transform guided by an approximated location of the
    singularities of the physical quantities in the U complex plane.  We also
    presented a Bayesian method to control the strong amplification of statistical
    noise during this procedure, using some simple non-perturbative information.
    The combination of singularity location, conformal transform crafting and
    Bayesian inference provides a robust and generic resummation methodology. 
    
    It was noticed recently<cit.> that for values of U inside
    the convergence radius, connected diagrammatic quantum Monte-Carlo techniques
    provide a systematic route for calculating the many-body quantum problem in a
    computational time that only increases  polynomially  with the requested
    precision. We argue that the argument of Ref. Rossi2017a can be directly
    extended to systems where the separation hypothesis holds (switching from
    working with the series in U to the series in W). We conclude that, in
    general, systems where the separation hypothesis hold can be computed with a
    computing time that increases polynomially with the requested precision.
    
    The approach presented here may have implications for a large class of other
    problems within or beyond condensed matter physics. In particular, a possible
    extension is to build a real time (equilibrium or non-equilibrium) quantum
    impurity solver for DMFT or its extensions, or directly addressing lattice
    problems such as the Hubbard model.  At its core, it consists in techniques to
    efficiently compute the bare perturbation series and to sum it. 
    Its limitations remain to be explored.  They could come from a resurgence of the sign problem,
    which would manifest itself in a very oscillatory nature of the integrals for
    expansion coefficients, making them hard to evaluate, or from a difficulty to
    sum the perturbative series, in particular for systems with a phase transition,
    or a non-Fermi liquid fixed point at low temperature.
    In order to address these questions, the technique needs to be applied
    to more complex models. Work is in progress in this direction.
    
    
    
    
    The Flatiron Institute is a division of the Simons Foundation.
    We acknowledge useful discussions with Laura Messio, Volker Meden and Christophe Mora.
    We thank our anonymous referee for pointing out Ref.<cit.>.
    We acknowledge financial support from the graphene Flagship (ANR FLagera GRANSPORT),
    the French-US ANR PYRE and the French-Japan ANR QCONTROL.
    
    
    § A TOY MODEL FUNCTION WITH A SINGULARITY ON THE REAL AXIS
    
    
    We present here on  Fig. <ref> 
    a toy model for the resummation of a function 
    f(U)= 1/ln(i(1-U)+1) that has a pole on the real axis at U=1 as well as a branch
    cut on the curve U=1-i (1+x) with x∈ [0,∞]. The aim of this toy
    model is to show that even though f(U) has a singularity on the real axis
    (and hence it will be difficult to calculate close to this singularity), it is
    possible to calculate the function beyond the singularity using a conformal
    transformation. We use the conformal map W = χ(U) - χ(0)/χ(U)
    - χ(0)^*, with χ(U) = i√((U-1)/p) - i that maps the inside of a
    parabola into the outside of the unit disk (see the upper left and right panels
    of Fig. <ref>).  The lower panel of
    Fig. <ref> shows the corresponding resummed series using
    N=10, 20 and 30 terms in the expansion of f(U).  Although we cannot
    calculate close to U=1, we find that with as little as N=20 terms in the
    expansion of f(U), we can recover an accurate description of f(U) for
    U>1.2 from an expansion around U=0.
    
    
    
    § CONVERGENCE OF THE PERTURBATION SERIES AT FINITE TIME
    
    
    In this appendix, we show that  at finite time t, the radius of convergence of the perturbation series for an operator O is infinite, 
    for a system with an interaction on a  finite number of sites and an infinite bath. 
    Indeed, the average is given by 
    
        ⟨Ô(t) ⟩∝⟨ T_c e^ -i U∫ du H_ int(u)Ô (t)
        ⟩,
    
    where the integral goes along the forward-backward Keldysh contour 0 → t → 0, 
    the operators are taken in the interaction representation, T_c is the
    usual Keldysh contour ordering operator and H_ int(u) is the
    interacting part of the Hamiltonian.
    
    
    
    
    More precisely, each of the 2^n terms of the expansion of the exponential has the form,
    
        U^n/n!∫_[0,t]^n du_1...du_n ⟨Ô (t) C(u_1,… u_n)⟩.
    
    where C is a product of c, c^†, and unitary time evolution operators.
    The terms ⟨Ô (t) C(u_1,… u_n)⟩ are amplitudes of probability for quantum processes and are therefore bounded.
    Explicitly, 
    
        ⟨Ô (t) C(u_1,… u_n)⟩   =
           Tr( e^-β H_0/Z_0Ô (t) C(u_1,… u_n) ) 
           =  ∑_ψ⟨ψ | e^-β H_0/Z_0Ô (t) C(u_1,… u_n) | ψ⟩
    
        | ⟨Ô (t) C(u_1,… u_n)⟩ |
              ≤∑_ψe^-β H_0/Z_0Ô (t) C(u_1,… u_n)  ψ
           ≤∑_ψe^-β H_0/Z_0Ô (t) C(u_1,… u_n)ψ
           ≤Ô,
    
    where v is the norm for a vector and the induced norm for an operator.
    We note that the norm is not modified by the unitary evolution e^iH_0 u Ae^-iH_0 u = A for any operator A, 
    and for the canonical operators c = 1, as can be checked in the Fock basis,  independently of the size of the bath.
    Since the norm is sub-multiplicative, we obtain the last inequality.
    
    Therefore, the term of order n in the expansion of (<ref>) is controlled by 
    a bound O(2 U t L)^n/n!, so the series
    has an infinite radius of convergence. Note that this argument is valid because
    the electron-electron interaction is present on a finite number of sites only.
    It would not apply directly to e.g. the Hubbard model in the thermodynamic limit.
    
    
    
    apsrev
