Reconstructing non-equilibrium regimes of quantum many-body systems from the analytical structure of perturbative expansions

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.


I. 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 lightmatter interaction on ultra-fast time scale 1 , light-induced superconductivity [2][3][4][5][6] or metal-insulator transition driven by electric field 7 , proposed e.g. to build artificial neurons 8 . 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 9,10 and overscreened 11,12 (multichannel) Kondo effects, characterized by non-Fermi liquid fixed points. Other notable examples of new quantum states induced by interactions are Luttinger liquids 13 that take place at edges in the fractional quantum Hall regime, or the "0.7 anomaly" [14][15][16] 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 17 .
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 18,19 .
Early approaches used a range of approximate techniques including 4th order perturbation theory 20 , equation of motion techniques 21 and the Non Crossing Approximation (NCA) 22 . State of the art techniques include the time-dependent Numerical Renormalization Group (NRG) and the density matrix renormalization group (DMRG) 19,[23][24][25][26][27][28] . Early attempts of real time quantum Monte Carlo [29][30][31][32][33] 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 [34][35][36][37][38] and the Schwinger-Keldysh diagrammatic Quantum Monte Carlo 39 (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 [40][41][42] (DMFT) which reduce bulk lattice problem to the solution of a self-consistent quantum impurity model.
A straightforward approach to study out-ofequilibrium 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 43 . Such a perturbative expansion is a central tool in quantum mechanics and quantum arXiv:1903.11646v2 [cond-mat.str-el] 29 Mar 2019 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" [44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59] 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 39 , and in the imaginary time Matsubara formalism [60][61][62][63] .
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. 1). 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. 1). We also introduce a Bayesian technique to take into account the amplification of the Monte-Carlo noise in the resummation process using some simple nonperturbative 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. 1). 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. 39 . The algorithm is discussed in details in a companion paper 64 , its imlementation is based on the TRIQS library 65 . 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 II introduces our notations for the single impurity Anderson model. Section III develops the resummation technique and illustrates it on the Kondo temperature. Section IV performs a benchmark of the method against NRG for the equilibrium dynamics. Section V 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 VI concludes this article and presents perspectives for our conformal approach to the perturbative expansions of strongly interacting quantum systems.

II. 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 66 . The Anderson impurity Hamiltonian reads: 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. 1. Here d is the on-site energy of the impurity (relative to the particle-hole symmetric point),n σ =ĉ † 0,σĉ 0,σ is the impurity density of spin σ electrons.ĉ † i,σ andĉ 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 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 nonequilibrium Keldysh formalism 67 , the non-interacting lesser and upper Green's functions are given by: 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), 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 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 68 , The retarded self energy Σ R (ω) is defined from the interacting Green's function by: Physical quantities have systematic expansion in power of U from which we obtain the corresponding quantity in the frequency domain by Fourier transform, First non zero orders of the self-energy series Σ R (U, ω) in powers of U for the equilibrium particle-hole symmetric Anderson model (real part in blue, imaginary part in red). This series has been computed with a real-time diagrammatic quantum Monte-Carlo method detailed in a companion article 64 . The statistical error is shown as shaded areas. Due to particle-hole symmetry, odd orders are zeros.
We obtain the functions G R n (ω) (typically up to n = 10) using the QMC algorithm of Ref. 39 and 64. The expansion of the self-energy is obtained from the G R n (ω) using a formal series expansion order by order of the Dyson equation (8). As an illustration, Fig. 2 shows the self-energy series, up to order 10, for the equilibrium particle-hole symmetric model as obtained from diagrammatic QMC 64 . 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.

III. 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. 69. Various techniques have been used in physics problems including Padé approximants 70 , Lindelöf extrapolation 52,71 or Cesàro-Riesz technique 48 . 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 40 , or e.g. Trilex 72,73 . 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 n=0 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. 74 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.
A. 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 75 . In a previous work 39 , 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. 1. 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. 1). 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 76 , 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 seriesF (W ) ≡ pF p W p defined by the com-positionF (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 seriesF (W 0 ) at the point of interest W 0 = W (U 0 ). Indeed, by construction W 0 ∈ D C and, sinceF (W ) is holomorphic in D C , D C is included in the convergence disk of the series F . Hence the seriesF converges at W 0 .
The result is independent of the choice of the domain C but the speed of convergence of the series forF (W 0 ) versus N is not, since it is determined by the relative position of W 0 compared to the radius of convergence RF ofF , i.e. η C ≡ |W 0 /RF |. 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 seriesF at W 0 , hence a clear failure of the method rather than a wrong result. Conversely the study of the convergence radius of thē 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 ofF (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. 39 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. 77 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. 77 for a detailed analysis. For a given U 0 and domain C, we now have to sum the transformed seriesF 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, N p=0 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 toF 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. 6.
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 sys- Poles of TK (U ) identified from the zeros of the 1/TK (U ) function. These are found by looking for the zeros of its truncated series. Here they are shown in the U/Γ complex plane with truncation at order 6 (red squares), 8 (blue points) and 10 (black stars). The black circle corresponds to |U | = RT K where RT K is the radius of convergence of the series of TK . The stable points ±i5Γ correspond to true nonperturbative poles of TK (U ). tematic 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 n=0 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 and the errors of the F 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. 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 distribu- and define the posterior distribution as, . (13) In practice, one proceeds as follows: (i) one generates many series according to P prior (F = f ) (this is straightforward since all the coefficients are assumed to be Gaussian and independent). (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.

B. 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: 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.

Singularities in the complex U plane
The dashed blue lines of Fig. 3 shows the truncated series of T K = N n=0 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. IV A for details on the used NRG implementation).
The inset of Fig. 3 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. 3. 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 78 ).
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. 4 for N = 6 (red squares), N = 8 (blue circles) and N = 10 (stars). One pair of zeros U ≈ ±i5Γ 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. 4 also shows the circle |U | = R T K extracted from the analysis of the T K (U ) series done in the inset of Fig. 3. We find that the two poles ±i5Γ 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 and the "parabola" map which is defined as where p is an adjustable complex parameter. Fig. 5 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 se-riesF p in the W variable, as explained above. The inset of Fig. 3 showsF n W n 0 (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 ofF (W ) whose radius is found to be RF ≈ 2. The final result T K (U ) using both transforms is shown in Fig. 3 and agrees quantitatively with the reference NRG results.
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. 5, 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 A. 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. 6 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. III A 3, our final result is obtained by averaging the different traces (thin red lines) with the weight given by Eq. (13). The right panel of Fig. 6 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. IV A for details).

C. Application to the Green's function and self energy
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. 2 where we plot the coefficients Σ R n (ω) 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 Σ R n (ω) is explained in the companion paper to this article 64 .
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. 7 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. 7 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. 7 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. 8. 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. 7. 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 The apparent convergence radius decreases with time. For small values of t, we can observe that the series coefficients decrease faster than exponentially, which indicates an infinite convergence radius. The thick dashed line shows the corresponding fit with (tΓ/2) n /n!. For large enough t, the series converges toward the steady state limit.
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).
First, at finite t, the radius of convergence of this series is infinite. Indeed, the average is given by where the integral goes along the Keldysh contour 0 → t → 0, the operators are taken in the interaction representation, T c is the usual Keldysh contour ordering operator andĤ int (u) is the interacting part of the Hamiltonian. At finite time t, the above expression is "just" the exponential of a (admittedly large) matrix times U , hence has an infinite radius of convergence, with an expected scaling of the coefficients controlled by O( c n n! ) with c some constant.
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: This behaviour was already noted in Fig. 14 of Ref. 39. It is also illustrated on Fig. 9, which shows various orders n of the expansion of the current through the dot versus 0.0 0.5 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 −1 K , see e.g. Ref. 79. 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 ).

IV. 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.

A. NRG implementation
The Numerical Renormalization Group (NRG) 80 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 81 was used to reduce finite size effects due to the NRG truncation. Second, symmetries of the problem were heavily exploited 82 , 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 Σ(ω) 83 , according to its exact representation as the ratio of two retarded correlation functions in the frequency domain: where G R (t) = −iθ(t) {d σ (0), d † σ (t)} is the usual single particle retarded Green's function in the time domain, and } 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.
B. Comparison to NRG in equilibrium Fig. 10 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. 8). 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. 11 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 31 . 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 64 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. 12 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.

V. 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.
A. Splitting of the spectral function Fig. 13 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 22 and perturbative 20 calculations that the Kondo resonance simultaneously broadens and get split into two peaks. Previous results on the spectral function 35 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. 14 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 [84][85][86] . 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. 13).

B. I-V transport characteristics
The left panel in Fig. 15 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 B T K the conductance experiences an extra suppression by the interaction (Coulomb blockade). We find a very good match with a previous calculation from Ref. 31. The present technique allows one to lift the main limitations that Ref. 30 and 31 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. 31. As a consequence, we can reach the low bias regime, which was not accessible in Ref. 31. Another important point is that the method is not limited to the symmetry point as we now demonstrate.
The right panel in Fig. 15 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. 12). As one increases further U , the Kondo width T K shrinks and the current decreases as Coulomb blockade starts to set in.

C. 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 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. 16. 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 [87][88][89][90] . 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.

VI. 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 chief 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. This methodology may have implications for a large class of other problems within or beyond condensed matter physics. It includes e.g. building a real time DMFT quantum impurity solver, or directly addressing lattice problems such as the Hubbard model. Work in these directions is in progress.
It was noticed recently 77 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. 77 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.