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

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

Rosario E. V. Profumo Univ. Grenoble Alpes, INAC-SPSMS, F-38000 Grenoble, France CEA, INAC-SPSMS, F-38000 Grenoble, France    Christoph Groth Univ. Grenoble Alpes, INAC-SPSMS, F-38000 Grenoble, France CEA, INAC-SPSMS, F-38000 Grenoble, France    Laura Messio IPhT, CEA, CNRS, URA 2306, 91191 Gif-sur-Yvette, France LPTMC, UMR 7600 CNRS, Université Pierre et Marie Curie, 75252 Paris, France    Olivier Parcollet IPhT, CEA, CNRS, URA 2306, 91191 Gif-sur-Yvette, France    Xavier Waintal Univ. Grenoble Alpes, INAC-SPSMS, F-38000 Grenoble, France CEA, INAC-SPSMS, F-38000 Grenoble, France
July 3, 2019

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 systemsLee and Ramakrishnan (1985), the Fermi- edge singularity in a quantum dotMatveev and Larkin (1992), a Kondo impurity embedded in an electronic interferometerTakada et al. (2014) and possibly the 0.7 anomaly in a quantum point contactThomas et al. (1996). While for a few situations, e.g. zero-dimensional (Kondo effects) and one-dimensional (Luttinger liquids) systems there exist exact analytical and numerical techniquesSchollwöck (2005); Bulla et al. (2008), 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.Gukelberger et al. (2015) 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 studiedProkof’ev and Svistunov (2008); Van Houcke et al. (2008); Kozik et al. (2010); Van Houcke et al. (2012), 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 systemsRubtsov et al. (2005); Werner et al. (2006); Gull et al. (2008, 2011), 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 Mühlbacher and Rabani (2008); Schiró and Fabrizio (2009); Schiró (2010); Werner et al. (2009, 2010). 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. Werner et al., 2009, 2010 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 Cohen et al. (2014a, b) 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 packageGroth et al. (2014) to tackle electron-electron interactions, or of the Triqs package Parcollet et al. (2014) 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 Rubtsov et al. (2005)), 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 I, 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 II 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 III 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 IV 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 V describes our QMC algorithm. In Section VI 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 VII, 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.

I Summary 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-timeWerner et al. (2010); Schiró and Fabrizio (2009); Schiró (2010) 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 problemShi and Zhang (2013); Rom et al. (1997). For instance in Ref. Werner et al., 2010 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.

Ii Models and basic formalism

Figure 1: 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 .

ii.1 Models

We consider a general time-dependent model for a confined nanoelectronic system connected to metallic electrodes, following the approach of Ref. Meir and Wingreen, 1992. A sketch of a generic system is given in Fig. 1. 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 Groth et al. (2014) and time-dependent properties Gaury et al. (2014). 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. (3) 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. (3). 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,




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. Rubtsov et al., 2005, 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 .

ii.2 Interaction 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 stepRammer and Smith (1986) 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 . The expansion in powers of can now be performed,


Of particular interest to us are one-particle observables (say current or electronic density) which can be directly expressed in terms of the lesser Green’s function at equal times:


Note at this stage that the following derivation is presented for the one particle correlator, but can be straightforwardly generalized to higher correlators.

To proceed, one evaluates the average of the (large) products of creation/destruction operators using Wick theorem, in a form of the determinant of a matrix,


where is given by


and the zeroth order term is . Eq. (14) 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. (14) 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 IV 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




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.

Iii The 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 extensivelyGaury et al. (2014) and we briefly summarize the main aspects here. Note that in Ref. Gaury et al., 2014, 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.

iii.1 General 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 Gaury et al. (2014),


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 ( sitesGaury and Waintal (2014)). The bound states contribution was not considered in Ref. Gaury et al., 2014 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 Groth et al. (2014) 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.

iii.2 Bound 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. (23) 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. (26) 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.

Iv Analysis 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. (14) 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. 2 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. (28) is convergent, but only that its coefficients are well defined.

Figure 2: as a function of for , , and . The calculations are performed using a direct evaluations of the integrals in Eq. (14) 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. 3 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. 3 should be contrasted with Fig. 4 which shows the same integrand but now summed over the four Keldysh indices. The integrand shown in Fig. 4 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. 5 shows the same as Fig. 3 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. 4 for . In the next section, we will use these observations to design a better sampling strategy for the Monte-Carlo method.

Figure 3: 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 .
Figure 4: Same parameters as in Fig. 3 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 .
Figure 5: Same parameters as in Fig. 3 but the integrand now uses the matrix instead of , i.e. is associated to the partition function instead of an observable. The colorplot represents

V Quantum 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).

v.1 Sampling 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. 4 and Fig. 5 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. (13)):


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. (32) 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




which relates the observable () to the output of the computations (). As always, Monte-Carlo computations only provide ratios between quantities, see Eq. (36). Here, we use our knowledge of the non-interacting observable (hence of ) to obtain the and finally the recursively. Eq. (36) 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.

v.2 Moves and detailed balance

Before we can actually perform calculations, we are left with a last task: designing a random walk that actually samples Eq. (30) 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 Prokof’ev et al. (1996); Rubtsov et al. (2005); Werner et al. (2006); Gull et al. (2008), 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.

v.3 Remarks

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.

v.3.1 Comparison with the sampling of the partition function

Our sampling strategy should be contrasted with the usual approach, used for instance in Ref. Werner et al., 2010, 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. (42) 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. (44) 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. (21) 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. 5.

(2). It is not guaranteed that the most probable configurations sampled by Eq. (44) 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. 3, Fig. 4 and Fig. 5.

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

v.3.2 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. (21) 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.

v.3.3 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 point