Quantum Monte-Carlo for correlated out-of-equilibrium nanoelectronics devices

Quantum Monte-Carlo for correlated out-of-equilibrium nanoelectronics devices

Abstract

We present a simple, general purpose, quantum Monte-Carlo algorithm for out-of-equilibrium interacting nanoelectronics systems. It allows one to systematically compute the expansion of any physical observable (such as current or density) in powers of the electron-electron interaction coupling constant . It is based on the out-of-equilibrium Keldysh Green’s function formalism in real-time and corresponds to evaluating all the Feynman diagrams to a given order (up to in the present work). A key idea is to explicitly sum over the Keldysh indices in order to enforce the unitarity of the time evolution. The coefficients of the expansion can easily be obtained for long time, stationary regimes, even at zero temperature. We then illustrate our approach with an application to the Anderson model, an archetype interacting mesoscopic system. We recover various results of the literature such as the spin susceptibility or the “Kondo ridge” in the current-voltage characteristics. In this case, we found the Monte-Carlo free of the sign problem even at zero temperature, in the stationary regime and in absence of particle-hole symmetry. The main limitation of the method is the lack of convergence of the expansion in for large , i.e. a mathematical property of the model rather than a limitation of the Monte-Carlo algorithm. Standard extrapolation methods of divergent series can be used to evaluate the series in the strong correlation regime.

The field of electronic correlations is largely dominated by applications to strongly correlated materials such as high- superconductors or heavy fermions. As a result, the large effort made by the community to build new numerical techniques to address correlations aims chiefly at reaching strongly correlated regimes for systems whose one-body dynamics is rather simple (the archetype of these systems being the Hubbard model). There are, however, many situations where the correlations are either small or moderate, yet their interplay with one-body dynamics might be very interesting. Examples include, for instance, the zero-bias anomaly in disordered systems[?], the Fermi- edge singularity in a quantum dot[?], a Kondo impurity embedded in an electronic interferometer[?] and possibly the 0.7 anomaly in a quantum point contact[?]. While for a few situations, e.g. zero-dimensional (Kondo effects) and one-dimensional (Luttinger liquids) systems there exist exact analytical and numerical techniques[?], the vast majority of these problems remains elusive to theoretical approaches. The aim of this article is to design a technique that could address moderate interactions for a large variety of out-of-equilibrium situations.

A natural route for dealing with electron-electron interactions is to compute the expansion of physical quantities in powers of the interaction coupling constant, denoted hereafter by . This expansion is traditionally written in terms of Feynman diagrams. One can then compute the first orders, or try various resummation strategies that have been elaborated in order to choose the relevant Feynman diagrams for a given problem.[?] From a numerical point of view, systematic expansions in powers of have also been intensively studied. In this context, various diagrammatic Monte-Carlo have been developed and studied[?], which aim at explicitly summing the series of Feynman diagrams numerically, for example for the self-energy. Concerning quantum impurity models, there has been an intense activity in the recent years in the development of new continuous (mostly imaginary) time quantum Monte-Carlo techniques, based on an expansion in (or around the strong coupling limit). These new algorithms are of huge practical value in solving the self-consistent impurity problems that arise from the dynamical mean-field theory of correlated bulk systems[?], even though they still suffer from the sign problem. They have been extended to the non-equilibrium case in a relatively straightforward way, simply adapting the Monte-Carlo method to the Keldysh formalism [?]. However, these out-of-equilibrium versions suffer from a severe dynamical sign problem, compared to their equilibrium counterparts, which has severely limited their usage in practice. In particular, they can not reach the long-time steady-state limit in several regimes of parameters. Moreover, the approach of Ref. has only been shown to work with sufficient accuracy for an Anderson impurity with particle-hole symmetry, i.e. a very special point of the phase diagram. More recently, bold diagrammatic Monte-Carlo for impurity models have also been extended to the Keldysh context and used in combination with the master equation for the density matrix [?] to reach longer time. Finally, testing these approaches in large systems, even at moderate interaction, remains also an open question.

In this paper, we first present a simple, systematic and general Monte-Carlo method to compute the first 10 to 15 coefficients of the expansion of any physical observable for a nanoelectronic system, in an out-of-equilibrium situation. The system can be a nanoscopic system connected to leads or a quantum impurity in a (possibly self-consistently determined) bath. Our method can be applied in various non-equilibrium contexts, either at short time after a quench, or in the long-time steady-state. In particular, it can easily reach the steady-state limit, even at zero temperature, as well as any intermediate time. It is also not limited to particle-hole symmetric case. The software developed can be seen as an extension of the Kwant package[?] to tackle electron-electron interactions, or of the Triqs package [?] to deal with non-equilibrium situations. Second, we discuss the issue of the summation of the perturbative series of the physical quantity, which is well-known to be a prominent topic in the quantum many-body problem. We will show that in the out-of-equilibrium Anderson model, for the parameters studied here, the current through the dot or the density on it have a finite apparent radius of convergence as a function of . We will also show that by simply modifying the quadratic part of the action (i.e. playing with the so-called term in [?]), one can significantly extend the radius of convergence, hence in practice compute for higher values of the interaction. Finally, we show that extrapolation technique for divergent series, e.g. Lindelöf method, can also significantly improve the range of applicability of our method.

In Section 1, we summarize our method and explain the main differences between our work and previous ones. This short section is mostly for QMC experts and can be skipped for people new to the field. Section 2 introduces our models and notations as well as the basic many-body perturbation expression that forms our starting point. This expression relates the interacting observables (such as current or magnetization) to the non-interacting Green’s function of the system. Section 3 discusses how to obtain the latter one, a prerequisite to any QMC scheme. While this step is relatively straightforward for simple impurity problems (the vast majority of the systems considered so far), its generalization to non-trivial geometries requires some care, or can become a computationally prohibitive task. Section 4 discusses a direct calculation of the first few orders of the interaction expansion by a brute force numerical integration. The discussion of the structure of the functions to be integrated will lead to a key insight for the Monte-Carlo. Section 5 describes our QMC algorithm. In Section 6 we use the QMC algorithm to calculate the first 10-15 terms in the expansion in powers of the interaction strength of the local charge on an Anderson impurity in equilibrium. We analyze the radius of convergence of the series in presence/absence of a mean-field term in the non-interacting Hamiltonian. In Section 7, we use the method in the ouf-of-equilibrium regime, to obtain some results associated with the Kondo effect. The article ends with a discussion and various appendices that contain some proofs or technical details.

1Summary of the approach

Let us briefly sketch our algorithm and its properties. Non-QMC experts can skip this part, since its content will be detailed and explained in the next sections. We start with a general Hamiltonian where is a non-interacting quadratic Hamiltonian of an infinite system (typically a nanoelectronic system connected to several electrodes) and contains the interacting part which is switched on at . We aim at calculating the expansion of an observable (say the current or the local occupation of an orbital) in powers of :

The are given by many-body perturbation theory in the Keldysh formalism in the form of a multi-dimensional integral of a determinant, according to Wick’s theorem:

The sum contains an -dimensional integral over internal times as well as a sum over Keldysh indices and a sum over the different interaction matrix elements. contains interaction matrix elements. is the determinant of a matrix built up with the non-interacting Green’s function of . Our algorithm works as follows:

(1) We compute directly the . In contrast, the imaginary-time or real-time[?] continuous-time algorithms sample the partition function of the problem . In the real-time Keldysh formalism however, the partition function is by construction and as we shall see, its sampling is not well suited for obtaining the . Technically, the integrand is concentrated around times which are close to while in the sampling of the are spread over the full interval . This first step ensures that our technique converges well as and that this limit can be taken order by order in .

(2) In the Keldysh formalism, one typically starts from a non interacting system, switches on the interaction, let the system evolve for a time , measures the observable and then evolves back to the non-interacting initial state. The Keldysh indices are reminiscent of the two evolutions [from 0 to () and back to 0 ()]. The time evolution is unitary. To keep this unitarity order by order, we choose to sum explicitly over the Keldysh indices. Hence our algorithm samples directly . Indeed, performing the sum over the Keldysh indices only with the Monte-Carlo Markov chain implies that unitarity is only respected on average (i.e. not for a single configuration), and in other words relies on the Monte-Carlo to perform massive cancellations. Obviously the explicit sum comes at an exponential cost: for each Monte-Carlo move one needs to calculate terms. However, we shall see that the gain in signal to noise ratio more than compensates for this additional computational cost. In the problems that we have computed so far, the dynamical sign/phase problem entirely disappears even for large times. From a diagrammatic perspective, the summation over Keldysh indices also implements the automatic cancellation of disconnected Feynman diagrams.

(3) We compute the coefficients individually rather than . Indeed, once the are obtained, we can analyze the convergence of the series for , which is a separate mathematical problem and has nothing to do with the Monte-Carlo or any other technique used to obtain the . Moreover, in continuous-time algorithms, the interaction very often takes the form of a density-density interaction and very special values of must be used for the computation not to be plagued by the sign problem[?]. For instance in Ref. the algorithm fluctuates wildly away from . Here we find that the convergence of the series for depends strongly on the value of . Note that this is a property of the perturbation series and therefore independent of the procedure used to obtain the coefficients. In practice, the dependence of the radius of convergence on the parameters can be used to access larger values of the interaction: we add (to ) and substract (to ) an on-site potential to our Hamiltonian such that the full Hamiltonian is unchanged but the series acquires a larger radius of convergence. This step allows us to tackle systems away from the particle-hole symmetry point. It is not linked to the QMC technique per se but to the choice of the initial quadratic Hamiltonian around which one performs the interaction expansion.

2Models and basic formalism

  Sketch of a typical mesoscopic system: a central interacting region (red) is connected to several (semi-infinite) non-interacting electrodes (blue) with finite temperatures T_i and chemical potentials \mu_i.
Sketch of a typical mesoscopic system: a central interacting region (red) is connected to several (semi-infinite) non-interacting electrodes (blue) with finite temperatures and chemical potentials .

2.1Models

We consider a general time-dependent model for a confined nanoelectronic system connected to metallic electrodes, following the approach of Ref. . A sketch of a generic system is given in Fig. ?. The Hamiltonian consists of a quadratic term and an electron-electron interacting term,

where the parameter controls the magnitude of the interaction. The non-interacting Hamiltonian takes the following form,

where () are the usual fermionic creation (annihilation) operators of a one-particle state on the site . The site index is general and can include different kinds of degrees of freedom: space, spin, orbitals. A crucial aspect is that the number of “sites” is infinite so that the non-interacting system has a well-defined density of states (as opposed to a sum of delta functions for a finite system) while interactions only take place in a finite region. Typically, the system will consist of a central part connected to semi-infinite periodic non-interacting leads. The dynamics of such non-interacting systems is well known and mature techniques exist to calculate both their stationary [?] and time-dependent properties [?]. The interaction Hamiltonian reads

In contrast to the non-interacting part, it is confined to a finite region. We also supposed that the interaction vanishes for negative time and is slowly or abruptly switched on at . A typical system described by Eq. is a quantum dot where electrostatic gates confine the electrons in a small, highly interacting, region while the electrodes have high electronic density, hence weak interactions.

The techniques described below are rather general and will be discussed within the framework of Eq. . The practical calculations however will be performed on the following Anderson impurity models. Model A corresponds to one interacting site, “0” inside an infinite one-dimensional chain,

where

is the level on-site energy, the (Zeeman) magnetic field and is the Heaviside function so that the interaction is switched on at . The hopping parameter is equal to unity for all sites except . We apply a bias voltage between the two (Left and Right) electrodes which are characterized by their chemical potential and and temperature . Model B is very close to model A with additional parameters , ,

One easily realizes that the two models are in fact equivalent in the stationary limit for : . However they have very different large limit at fixed (small) : model A corresponds to the degeneracy point between 0 and 1 electrons on the impurity (where Coulomb blockade is lifted) while model B corresponds to the Kondo regime. More importantly, the perturbation series in powers of of the same observable will be different between these two models, with different convergence radius for fixed . The parameters have been introduced in Ref. , to improve the sign problem in imaginary-time Quantum Monte-Carlo. An important energy scale for these models is the (non-interacting) tunneling rate from the impurity to the reservoirs. It is given by with .

2.2Interaction Expansion

Our starting point for this work is a formal expansion of the out-of-equilibrium (Keldysh) Green’s function in powers of electron-electron interactions. This is a standard step[?] which we briefly sketch to introduce our notations.

Using the interaction representation, one defines where is the evolution operator from to associated with . Introducing the Keldysh index , one defines the contour ordering for pairs : for all , if and if . The contour ordering operator acts on products of fermionic operators labeled by various “contour times” and reorder them according to the contour ordering: if and if . The non-interacting contour Green’s function is defined as

where is just , the Keldysh index serving only to define the position of the operator after contour ordering. The contour Green’s function has a matrix structure in which reads

where , , and are respectively the time ordered, lesser, greater and anti-time ordered Green’s functions. Efficient techniques to obtain these non-interacting objects for large systems will be discussed in the next section. Last, one defines the full Green’s function with definitions identical to the above except that is replaced by , the evolution operator associated to the full Hamiltonian . The fundamental expression for reads

where the integral over is taken along the Keldysh contour, i.e. increasing for and decreasing for . is equal to with the operators replaced by .

Eq. might look cumbersome at first sight, yet it is a compact expression: provided one knows how to calculate non-interacting Green’s functions (which will be taken care of in the next section), Eq. expresses the full interacting Keldysh Green’s function, hence the physical observables, in terms of integrals and sums of determinant of known quantities. All that remains is to find a suitable numerical way to perform those integrals and sums. For a local interaction present on sites, calculating all contributions to order in the stationary regime corresponds to a numerical complexity of the order of where the measurement time has to be large enough for the effect of the electron-electron interaction to be well established. This can be performed using standard integration routines for the first few orders (In Section 4 we calculate contributions for model A) but becomes quickly prohibitive for larger values of . For larger orders, a stochastic sampling of the integrals is compulsory.

To simplify the notations, we introduce the notion of configuration

and note

Introducing,

we get the following compact expression:

where the factor has dropped out due to the ordering of the . Note that in the Keldysh formalism the partition function is unity which translates into

where the matrix is identical to with the last row and column deleted. Actually, a much stronger statement can be made on : for any and configuration , one has,

The proof is straightforward and standard: one first locates the largest time in the configuration , say . When goes from to , the ordering of with respect to the other times is unchanged ( is larger than all the times on the upper part of the contour and smaller than all those on the lower part of the contour), hence the contour Green’s functions are unchanged and the matrix is also unchanged. As a result of the sign these two contributions cancel each other.

3The non-interacting Green’s function

In order to proceed with evaluating the interaction corrections to observables, the first step is an efficient way to calculate the various real-time non-interacting Green’s functions of the problem. For a small dot problem or a DMFT model, this step is easy. For larger systems, this question is more delicate, and it has been studied extensively[?] and we briefly summarize the main aspects here. Note that in Ref. , only quantum transport was of interest so that contributions coming from bound states could have been omitted. Here however, they will have to be taken into account properly.

3.1General method

Our starting point for calculating non-interacting Green’s functions is an expression that relates them to the (Scattering) wave functions in the system [?],

Here, labels the various propagating channels of the leads, the scattering state at energy (in the electrode) and the corresponding Fermi distribution function. labels a bound state of energy and wave functions . The greater Green’s function is obtained with an identical expression with the Fermi functions replaced by . Efficient techniques for calculating the scattering wave functions have been designed so that these objects can be obtained for large systems ( sites[?]). The bound states contribution was not considered in Ref. and will be discussed below. The actual calculations performed in this article were restricted to a stationary non-interacting system, where the above expression further simplifies into

Here again, the stationary wave functions are standard objects. They are in fact direct outputs of the Kwant software [?] which we use for their calculations. Once the lesser and greater Green’s functions are known, one completes the Keldysh matrix with the standard relations

To obtain those Green’s function numerically, i.e. for many values of , one needs to perform the integration over the energy many times. In practice, the stationary wave functions are calculated once using Kwant and cached. The integration itself is performed using standard numerical routines. For the single site model A or B, the above technique in its full generality can be avoided: one can simply compute the Green’s function in energy analytically and perform a numerical Fourier transform. We have checked explicitely that both techniques provide identical non-interacting Green’s functions in this special case.

3.2Bound states contribution

The presence of the electrodes in the system is very important physically: it provides the system with a relaxation mechanism. Mathematically, the integral in Eq. mixes nearby energies so that the resulting non-interacting Green’s functions decay (and oscillate) at large times. However, in presence of a large enough confining energy (far from zero parameter in model A), true bound states can appear in the system. They have energies outside of the electrode bands and therefore cannot hybridize with the plane waves of the electrodes. They satisfy the stationary Schrodinger equation for the infinite system. Upon integrating over the electrode degrees of freedom, they satisfy a simpler (yet non-linear) equation for the interacting region only:

where is the retarded self-energy due to the electrode. For a practical calculation, we do as follows: first we truncate and keep the interacting region plus a rather large (yet finite) fraction of the electrodes. We diagonalize the corresponding finite matrix and locate the eigenvalues that are outside the conducting bands of the electrodes. These eigenvalues are used as initial guess and we compute the bound states by iteratively solving Eq. until convergence. Note that there is an easy check to make sure that one uses a complete basis of the problem: one must have,

which is not verified if some bound states are forgotten. Note also that in most of this article, we focus on situations where there are no bound states in the system. This can be easily achieved by using leads which have a larger bandwidth than that of the central system, so that any bound state that could take place there hybridize with the continuum of the lead.

4Analysis of the first terms of the perturbative series

Knowing how to get the non-interacting Green’s functions, we are now ready to calculate the perturbation series. A first, rather naive, technique would consist in calculating the integrals in Eq. using a simple discretization scheme (Simpson in our case). Only the first few orders can be obtained that way, at large computational cost. Nevertheless, it is rather instructive and also serves as a check for the QMC algorithms discussed in the next section. We focus on model A and compute the local charge at various orders in ,

Fig. ? shows the resulting for . With a parallel implementation, results for can also be obtained (not shown) at important computational cost and is prohibitive. All these results will be reproduced using the quantum Monte-Carlo sampling with a tiny fraction of the computational time. Note that the stationary results obtained for do not mean that the series Eq. is convergent, but only that its coefficients are well defined.

 Q_n as a function of \epsilon_d for n=0, 1, 2 and 3. The calculations are performed using a direct evaluations of the integrals in Eq.  using the Simpson rule for t=20, \gamma=0.5 and T=0.
as a function of for , , and . The calculations are performed using a direct evaluations of the integrals in Eq. using the Simpson rule for , and .

It is very instructive to have a look at the quantity which is actually integrated to obtain the . Fig. ? shows the integrand of for the 4 values of the pair of Keldysh indices . We see that this integrand decays slowly as a function of and even more slowly as or get away from the time where the charge is measured. The sign of the integrand changes as one changes the Keldysh indices. Fig. ? should be contrasted with Fig. ? which shows the same integrand but now summed over the four Keldysh indices. The integrand shown in Fig. ? now decays fast as or gets away from . This observation can be proven and generalized for higher orders: the integrand decays to 0 when a group of is far from the time where the physical observable is measured, Cf. Appendix B. Finally, Fig. ? shows the same as Fig. ? but for the matrix associated with the partition function. Note that for , the sum on the Keldysh indices simply vanishes, so there is no analogous Figure as Fig. ? for . In the next section, we will use these observations to design a better sampling strategy for the Monte-Carlo method.

 Colorplot of the integrand of Q_2 as a function of the two times u_1 and u_2 for model A with \mu_L=\mu_R=0, \epsilon_d=0, T=0 and t=10. The four panels correspond to the 4 possible values of the two Keldysh indices a_1 and a_2. The explicit form of the integrand is f(u_1,u_2,a_1,a_2)= -\Im m (-1)^{\sum_i a_i}  \det  \mathrm{{\textbf{M}}} _{2}(u_1,u_2,a_1,a_2).
Colorplot of the integrand of as a function of the two times and for model A with , , and . The four panels correspond to the 4 possible values of the two Keldysh indices and . The explicit form of the integrand is .
 Same parameters as in Fig.  but the integrand has now been summed over Keldysh indices. The colorplot represents f(u_1,u_2)= i \sum_{a_1,a_2} (-1)^{\sum_i a_i}  \det  \mathrm{{\textbf{M}}} _{2}(u_1,u_2,a_1,a_2) (f is real). Note that the integrand is now real, positive and concentrated around u_1 = u_2 = t.
Same parameters as in Fig. but the integrand has now been summed over Keldysh indices. The colorplot represents ( is real). Note that the integrand is now real, positive and concentrated around .
 Same parameters as in Fig.  but the integrand now uses the matrix \mathrm{{\textbf{P}}} _2 instead of \mathrm{{\textbf{M}}} _2, i.e. is associated to the partition function instead of an observable. The colorplot represents f(u_1,u_2,a_1,a_2)= -\Im m (-1)^{\sum_i a_i}  \det  \mathrm{{\textbf{P}}} _{2}(u_1,u_2,a_1,a_2)
Same parameters as in Fig. but the integrand now uses the matrix instead of , i.e. is associated to the partition function instead of an observable. The colorplot represents

5Quantum Monte-Carlo

The direct method of the previous section works in principle but is limited in practice to very small orders due to its prohibitive computational cost. Stochastic methods, such as the Metropolis algorithm, can be extremely efficient at calculating integrals in high dimensions. In this section, we propose a new route to sample the interacting series by constructing a Markov process in the Fock configuration space (i.e. that not only samples the integrals themselves but also samples the various orders within one process).

5.1Sampling strategy

Our algorithm is inspired by the conclusion of the previous section. It consists in i) sampling directly the physical quantity to be computed (and not the partition function, which is anyway in the Keldysh formalism), and ii) summing explicitly over the Keldysh indices to restore unitarity (the symmetry between the two Keldysh contours) for all configurations. Indeed, it is clear from the contrast observed between Fig. ? and Fig. ? that it is a much better choice, since the integration region is in the first case well localized around the time at which the quantity is computed. Sampling would result in sampling large regions of the Fock space which are irrelevant to the actual observable.

We introduce:

which is a real number, as proven in Appendix C. We construct a Markov process that samples the following density of probability,

where the denominator ensures that the probability is normalized. We use the notation and to show explicitly that the “partition function” and interaction parameter belong to the QMC technique. In particular, the physical value of the interaction can (and will) be distinct from the one(s) used in the Monte-Carlo; also the physical partition function is in the Keldysh formalism.

Introducing , the contribution to order to the observable (see Eq. ):

the terms of the perturbative expansion are given by

where stands for the average over the probability . All is left is to construct a Markov process that samples this distribution. A key aspect of the approach is that is sampled for several (at least two) values of simultaneously: in the average , varies between and . Introducing

we find that the probability to be in the order (in practice the fraction of the Monte-Carlo spent in configurations at order ) is

The normalization of the total probability provides the partition function in terms of the , . Note that the are by definition independent of the QMC technique used to calculate them and in particular of the value of . is simply given by the non-interacting value of the observable: . The last item to introduce is the probability for the fluctuating sign in Eq. to be (and to be ). Note that is always real (Cf Appendix C) so that we only average fluctuating signs, not phases. We note so that is the average sign at a given order. and are the direct outputs of the computations. With these notations, and one finally arrives at

and

which relates the observable () to the output of the computations (). As always, Monte-Carlo computations only provide ratios between quantities, see Eq. . Here, we use our knowledge of the non-interacting observable (hence of ) to obtain the and finally the recursively. Eq. needs to be applied for all up to the maximal order needed, but its evaluation for various needs not to be done within the same QMC run.

5.2Moves and detailed balance

Before we can actually perform calculations, we are left with a last task: designing a random walk that actually samples Eq. with varying between and . This step is very similar to the construction of other continuous-time QMC. We introduce two sorts of moves: the moves where we increase by one unit by adding one vertex and the moves where one vertex is deleted. The algorithm to obtain the configuration at step from the configuration goes as follows:

(i) Move selection. We choose the move () with probability () with . When () we have () otherwise we typically choose .

(ii) Move . We select uniformly in . We select uniformly among the different terms . The overall probability to propose the move is

(iii) Move . We select the vertex to be removed uniformly between . The overall probability to propose the move is

(iv) Detailed balance. We choose the acceptance probability and so that it satisfies the detailed balance equation

Using the Metropolis algorithm this leads to

Note that, as usual in continuous-time quantum Monte-Carlo methods [?], the factor is present on both sides of the detailed balance equation and eventually drops, so that the Monte-Carlo can be performed directly in the (time) continuum.

5.3Remarks

We now have a complete practical scheme for calculating many-body perturbation to a given observable. Before we embark in concrete examples, let us make a few remarks.

Comparison with the sampling of the partition function

Our sampling strategy should be contrasted with the usual approach, used for instance in Ref. , which has its origin in the imaginary-time techniques and where the (density of) probability to be in the configuration with the Keldysh indices is given by

i.e. one samples the matrix (instead of ) and one samples the Keldysh indices (instead of summing on them exactly and explicitely). In this scheme, an observable (say the charge , possibly resolved in spin) is given by

while the partition function is given by

Constructing a Markov process that samples Eq. can be done similarly to the construction presented in the previous subsection. The observable can be estimated by taking the ratio of the above two equations.

Although this approach has shown some success, one can make the following remarks.

(1). The sign of can fluctuate strongly so that the statistical average in Eq. is very small and suffers from a very bad signal to noise ratio. This is the sign problem that plagues Quantum Monte-Carlo techniques for fermions. Eq. indicates that this problem is most probably worse in presence of the Keldysh indices where one expects wildly fluctuating signs. This is shown by the data in Fig. ?.

(2). It is not guaranteed that the most probable configurations sampled by Eq. are also the ones that contribute most to . On the contrary: the determinant of depends only on the relative positions of the with respect to the others (it is essentially a sum of terms of the form ) and not at all of . The integrals contributing to on the other hand decay when the get away from . Hence for large times, the above scheme samples values of the very far from which therefore contribute very little to the actual observable. This was shown explicitly in Fig. ?, Fig. ? and Fig. ?.

(3). The signal to noise ratio usually deteriorates rapidly with in these algorithms making it difficult to reach the stationary regime.

Role of the explicit sum over the Keldysh indices

The sampling of discussed above does not preserve the symmetry between the two parts of the Keldysh contour for a given configuration (it preserves it in average, as it should): in the Keldysh formalism, one starts from a non-interacting density matrix, switches on the interaction for some time , measures the observable, then unwinds the effect of the interaction until one is back to the original density matrix. Here however, a given configuration might have a few Keldysh indices on one branch and the rest on the other, meaning that a typical configuration does not enforce the symmetry of Keldysh indices, which reflects unitarity. From a different perspective, the corresponding expansion includes all the Feynman diagrams, including the disconnected diagrams although they have a vanishing contribution to the observables.

In our scheme in contrast, we explicitly sum over all Keldysh indices. Eq. indicates that all contributions from disconnected diagrams explicitly drop from the calculation. We calculate only connected diagrams which should be advantageous. From a technical perspective, one finds that the quantity is always real (not complex, see Appendix C) so that one averages signs instead of phases. One expects that the resulting potential sign problem is milder than the so-called phase problem which originates from averaging a random phase.

However, our scheme has an obvious drawback: one evaluation of corresponds to the calculation of determinants, so that this algorithm complexity now increases exponentially with the maximum order . We show below results for up to and it is reasonable to expect that one could calculate up to . Usual algorithms to calculate determinants have complexities which scale as . We show in Appendix D that the fast updates of the determinants (with complexities that scale as ) can easily be extended to the present case using Gray code, so that the overall complexity of our algorithm scales as . Actually, “mirror” Keldysh configurations have equal contributions (see Appendix C) so that only configurations need to be included. Overall, we shall see that the additional computational complexity is more than compensated by the important gain in signal to noise ratio.

Statistical errors and the sign problem

In practice, we calculate the moments sequentially with a separate QMC computation for each value of (typically with , or , ) so that one gets a fully controlled (optimally flat) histogram of the various orders. Starting from (known without error), we iteratively compute from different runs that use different interaction strengths ,

One can use for instance such that (this is always possible but not strictly necessary as long as remains of order unity).

Note that a naive scheme where one would try to evaluate all the values of in a single run would run into an artificial difficulty: in a single run, the histogram is usually sharply peaked around one value of and one cannot get a good statistics both for and . A small statistics in, say leads to a large error in the estimate of which further corrupts the evaluation of all , hence .

Coming back to our scheme, the error made on the estimate of is bounded by where is the number of independent points. Hence we find that the (one standard deviation) error on is bounded by

which can be controlled to arbitrary precision provided . Hence, the calculation of the , which involves only positive numbers, can be done with extremely good accuracy.

The limitation to the precision – the sign problem – takes its origin in the average of the sign contained in Eq. . Note that in contrast to other techniques, this sign is here in the numerator and is not present in the denominator, i.e. a small sign does not necessarily mean a sign problem: it can simply mean a small value of . This analysis is close to e.g. Refs . The error made on the sign is given by a Bernoulli law ( with probability , with probability ), hence is given by which is always smaller than

(the upper bound is reached when the sign becomes small). Putting everything together, implies that and we arrive at the relative error,

In this form, it seems that the smaller the sign, the larger the error, hence the sign problem. However, one must remember that while and depend on the QMC algorithm, their product does not, so that the error can be recast into

In this second form, it becomes apparent that a bad sign problem (small ) is equivalent to a bad sampling choice which leads to a large . The behaviour of the error is therefore intimately linked with the growth of with which itself depends strongly on the actual probability sampled. For instance, if one samples the sum over Keldysh indices (instead of summing explicitely over the indices as we do), one gets identical but much larger and consequently much smaller . In fact the would increase by more than a factor so that the overall method would be much less efficient than the one proposed here. The global relative error therefore contains three contributions, which reflect respectively the total computing time (), the intrinsic physics of the problem () and the quality of the choice of the sampling procedure (). From the above analysis, the colorplots shown in Fig. ?, Fig. ? and Fig. ? take a different meaning. Indeed, one can see that (the integral of the function displayed in Fig. ?) is rather small: the function decays rather quickly when get away from . If on the other hand we have chosen to sample the Keldysh indices in addition (as in the standard schemes), then would have been the sum of the integral of (the absolute value of) the different panels of Fig. ?. One immediatly realizes that the signal to noise ratio would have been much smaller.

Convergence of the interacting series

Our approach separates the calculation of the from the study of (the convergence of) the series itself. This could be used to obtain the full dependence of the observable, but more importantly, it allows one to disentangle physical aspects (for instance the convergence or lack of of the series) from technical ones (the calculation of its elements). The convergence of the series will be discussed next. In principle, one could use various resummation procedures such as Pade approximant, Lindelöf analytical continuation and/or Borel resummation in order to extrapolate the series from its first coefficients. An example of such a procedure is given in the appendices.

It is important to note already at this stage that the parameter , as introduced in model B, plays a crucial role in the algorithms sampling the partition function at equilibrium [?] as only special values of are free of the sign problem. It is also known that the typical perturbation order in those algorithms is strongly reduced by using the best value of . We shall find in the next section that the parameter has a drastic influence on the convergence of the series but not on the actual calculation of the itself.

6First results: analysis of the series convergence

 QMC results for model A at \gamma=1/2, T=0 and t=10 for \epsilon_d=0 (green squares) and \epsilon_d=-0.5 (blue circles), as a function of the order n. Red diamonds: model B with \gamma=1/2, T=0 and t=10 for \epsilon_d=-0.5 and \alpha=0.5. Top panel: Q_n, central panel: c_n and bottom panel: average sign s_n.
QMC results for model A at , and for (green squares) and (blue circles), as a function of the order . Red diamonds: model B with , and for and . Top panel: , central panel: and bottom panel: average sign .
 QMC results for model A at \gamma=1/2, T=0 and t=10. Top panel: charge Q(N, U) as a function of U, for different N and for \epsilon_d=0: N=5 (Black), N=7 (red), N=13 (blue) and N=15 (green). Bottom panel: Q(N, U) as a function of 1/N for different U: U=0.2 (Orange), U=0.4 (red), U=0.6 (black) and U=0.8 (blue).
QMC results for model A at , and . Top panel: charge as a function of , for different and for : (Black), (red), (blue) and (green). Bottom panel: as a function of for different : (Orange), (red), (black) and (blue).

6.1Bare results

As a first application, we compute the interaction corrections to the charge on the impurity in model A. Fig. ? shows the for up to as well as the corresponding and . We find that the magnitudes of the do not appear to decrease with but rather remain of order unity . This implies (Hadamard’s theorem) that the series has an apparent convergence radius of order unity (apparent because it relies on extrapolating the behaviour of the first known coefficients to large orders). We can already notice that the curve for model B () has a very different behaviour with a fastly decreasing hence ; this aspect will be discussed later in the text. Fig. ? shows the truncated series

as a function of (upper panel) and (lower panel). We find a nice convergence for but a divergence beyond, i.e. we cannot access the physics beyond () by simply summing the series. This divergence has nothing to do with the QMC technique which is just a way to calculate the - it belongs to the physics of the problem. In the following, we will discuss two ways to bypass this problem: performing the interaction expansion from a different starting point and resumming the series by moving its singularities away from the expansion point.

6.2Hartree-Fock series

Before going on with the QMC results, let us analyse a simpler, approximate, series which is obtained from the Hartree-Fock approximation (which reduces to Hartree for model A). This series will be useful in identifying a possible source for the observed apparent radius of convergence. Introducing,

the Fourier transform of the retarded non-interacting Green’s function , the non-interacting charge at equilibrium and takes the form

In the Hartree approximation, one replaces the on-site energy by its mean-field value. The fully self-consistent Hartree would be defined as . Here however, we restrict ourselves to summing the ladder of tadpole diagrams which is sufficient to illustrate our point: . The corresponding series is given by with

Note that the first two moments are the exact ones: and . On the other hand, from the above construction Eq.(Equation 15), we find that has branch cuts in the complex plane, given by

so that one expects the series in powers of to have a finite radius of convergence given by the branch closest to , i.e. equal to unity (see the inset of Fig. ?). Indeed, Fig. ? shows an example of the partial sums which diverges around . This is very reminiscent to what we have found for model A. This calculaiton illustrates in a simple and tractable approximation that a finite radius of convergence is due to the existence of singularities or branch cuts in the complex plane.

 Hartree-Fock series of model A as a function of U: Exact curve Q^{HF}(U) (thick line) and partial sums Q^{HF}(N,U) for N=10,20,30, \epsilon_d=0 and \gamma=1/2. The series has a divergence at U=1. The inset shows the analytical structure of Q^{HF}(U) in the complex plane ({\rm Re}\ U, {\rm Im}\ U): branch cut of Q^{HF}(U) (thick magenta line) and convergence radius of the Hartree-Fock series (dashed line).
Hartree-Fock series of model A as a function of : Exact curve (thick line) and partial sums for , and . The series has a divergence at . The inset shows the analytical structure of in the complex plane : branch cut of (thick magenta line) and convergence radius of the Hartree-Fock series (dashed line).

6.3Using a different non-interacting problem

In this section, we show that by playing with the parameter, one can greatly enhance the radius of convergance of the series hence access correlated regimes. Instead of starting the perturbation from the Hamiltonian, one can incorporate the mean-field treatment into the non-interacting Hamiltonian so that only the fluctuations of the interaction need to be taken into account in the perturbation. In fact, one can even push this idea further and add an arbitrary quadratic Hamiltonian to and remove it accordingly from . The idea is to start with one-body propagators as close as possible to the interacting ones, so that the role of the perturbation becomes very weak. Within our current implementation we can easily add an on-site potential to the non-interacting Hamiltonian and use in the perturbation where the new parameter is our targetted value of the interaction (see the definition of model B in Eq. ). For , one recovers the original model A for , hence the corresponding results in the long-time limit.

The parameters are exactly the same as the ones used in the interaction expansion continuous-time quantum Monte-Carlo introduced by Rubtsov [?], in equilibrium. In that algorithm, it was shown [?] that the sign problem strongly depends on the value of . It also reduces the average order of perturbation of these QMC methods Here, we will now show that the apparent radius of convergence of the interaction expansion strongly depends on these parameters, and we will use this to our advantage.

The technique is illustrated in Fig. ? first for a value of that could be reached with the initial approach () and secondly for a value that could not be reached (). We find that this approach works remarkably well: not only we can recover the former results, but we can also go to regimes that were not accessible before. As a self-consistency check, we find that the results do not depend on for . Of course a disadvantage of this approach is that one must perform a separate computation for each value of needed.

 Role of the non-interacting Hamiltonian around which is done the expansion; Q(U) at \epsilon_d=0 with an extra potential \delta\epsilon_d in the non-interacting Hamiltonian and a corresponding compensating term \alpha=\delta\epsilon_d/\bar U in the perturbation. The green line corresponds to the original series with \delta\epsilon_d=0 in its regime of convergence. All curves for \gamma=1/2, T=0 and t=10.
Role of the non-interacting Hamiltonian around which is done the expansion; at with an extra potential in the non-interacting Hamiltonian and a corresponding compensating term in the perturbation. The green line corresponds to the original series with in its regime of convergence. All curves for , and .

The faster convergence of the new series can be seen in Fig. ? where we compute the corresponding series versus (upper panel) as well as the convergence of the partial sum versus (lower panel). We find that only 2 or 3 orders are sufficient to obtain the exact result provided one uses a “starting point” close enough to the final solution. A pragmatic way to perform the calculation is therefore to optimize the value of so that the corresponding decreases as fast as possible.

 Top panel: absolute value of the charges Q_n for different values of \alpha,\delta\epsilon_d, at \epsilon_d=0, \gamma=1/2, T=0 and t=10: \alpha=\delta\epsilon_d=0 (Black), \alpha=\delta\epsilon_d=0.5 (Blue), \alpha=\delta\epsilon_d=0.6 (Green), \alpha=\delta\epsilon_d=0.7 (Red). Bottom panel: corresponding dependence of the charge Q(N,U=1) as a function of 1/N.
Top panel: absolute value of the charges for different values of ,, at , , and :