Orthogonal Polynomial Representation of Imaginary-Time Green’s Functions

# Orthogonal Polynomial Representation of Imaginary-Time Green’s Functions

Lewin Boehnke I. Institut für Theoretische Physik, Universität Hamburg, D-20355 Hamburg, Germany    Hartmut Hafermann Centre de Physique Théorique, Ecole Polytechnique, CNRS, 91128 Palaiseau Cedex, France    Michel Ferrero Centre de Physique Théorique, Ecole Polytechnique, CNRS, 91128 Palaiseau Cedex, France    Frank Lechermann I. Institut für Theoretische Physik, Universität Hamburg, D-20355 Hamburg, Germany    Olivier Parcollet Institut de Physique Théorique (IPhT), CEA, CNRS, URA 2306, 91191 Gif-sur-Yvette, France
###### Abstract

We study the expansion of single-particle and two-particle imaginary-time Matsubara Green’s functions of quantum impurity models in the basis of Legendre orthogonal polynomials. We discuss various applications within the dynamical mean-field theory (DMFT) framework. The method provides a more compact representation of the Green’s functions than standard Matsubara frequencies and therefore significantly reduces the memory-storage size of these quantities. Moreover, it can be used as an efficient noise filter for various physical quantities within the continuous-time quantum Monte Carlo impurity solvers recently developed for DMFT and its extensions. In particular, we show how to use it for the computation of energies in the context of realistic DMFT calculations in combination with the local density approximation to the density functional theory (LDA+DMFT) and for the calculation of lattice susceptibilities from the local irreducible vertex function.

###### pacs:
71.27.+a, 71.10.Fd

In recent years, significant progress has been made in the study of strongly-correlated fermionic quantum systems with the development of methods combining systematic analytical approximations and modern numerical algorithms. The Dynamical Mean-Field Theory (DMFT) (for a review see Ref. Georges96, ) and its various extensions. MaierReview (); toschi (); dualfermion1 (); clusterdualfermion (); ldfa () serve as successful examples for this theoretical advance. On the technical side, important progress was made in the solution of quantum impurity problems, i.e. local quantum systems coupled to a bath (self-consistently determined in the DMFT formalism). In particular, a new generation of continuous-time quantum Monte Carlo (CTQMC) impurity solvers Rubtsov05 (); Werner06 (); WernerPRBLong (); CTAUX () has emerged that provide unprecedented efficiency and accuracy (for a recent review, see Ref. ctqmcrmp, ).

In practice, several important technical issues still remain. Firstly, while the original DMFT formalism is expressed in terms of single-particle quantities (Green’s function and self-energy), two-particle quantities play a central role in the formulation of some DMFT extensions (e.g. dual-fermions dualfermion1 (); dualfermion2 (); gangli (); brener (), Dtoschi ()) and in susceptibility and transport computations in DMFT itself. They typically depend on three independent times or frequencies, and spatial indices. Therefore, they are quite large objects that are hard to store, manipulate and analyze, even with modern computing capabilities. Developing more compact representations of these objects and using them to solve, e.g., the Bethe-Salpeter equations is therefore an important challenge.

A natural route is to use an orthogonal polynomial representation of the imaginary-time dependence of these objects. While the application of orthogonal polynomials has had productive use in other approaches to correlated electrons, KPM (); CheMPS () in this paper we show how to use Legendre polynomials to represent various imaginary-time Green’s functions in a more compact way and show their usefulness in some concrete calculations.

A second aspect is that modern CTQMC impurity solvers still have limitations. One well-known problem is the high-frequency noise observed in the Green’s function and the self-energy (see e.g. Fig. 6 of Ref. Gull07, ). Even though this is in general of little concern for the DMFT self-consistency itself, it can become problematic when computing the energy, since the precision depends crucially on the high-frequency expansion coefficients of the Green’s function and self-energy. An important field of application involves realistic models of strongly correlated materials through the combination with the local density approximation (LDA+DMFT). KotliarRMP06 () In this paper, we show that physical quantities such as the Green’s function, kinetic energy, and even the coefficients of the high-frequency expansion of the Green’s function can be measured directly in the Legendre representation within CTQMC and that the basis truncation acts as a very efficient noise filter: the statistical noise is mostly carried by high-order Legendre coefficients, while the physical properties are determined by the low-order coefficients.

This paper is organized as follows: Section I is devoted to single-particle Green’s functions. More precisely, in Sec. I.A, we introduce the Legendre representation of the single-particle Green’s function and how it appears in the CTQMC context; we then illustrate the method on the imaginary-time (I.B) and imaginary-frequency (I.C) Green’s function of a standard DMFT computation; in I.D, we discuss the use of the Legendre representation to compute the energy in a realistic computation for SrVO. Section II is devoted to two-particle Green’s functions: We first present the expansion in Sec. II.A and illustrate it on an explicit DMFT computation of the antiferromagnetic susceptibility in Sec. II.B, followed by the example of a calculation of the dynamical wave-vector resolved magnetic susceptibility. Additional information can be found in the appendixes. Appendix A gives some properties of the Legendre polynomials relevant for this work. Appendix B discusses the rapid decay of the Legendre coefficients of the single-particle Green’s function. Appendix C first derives the accumulation formulas for the single-particle and two-particle Green’s functions in the hybridization expansion CTQMC (CT-HYB) algorithm Werner06 () (while these formulas have been given before, Werner06 (); Gull07 () the proof presented here aims to explain their resemblance to a Wick’s theorem). We then give the explicit formulas in the Legendre basis. For completeness, we provide an accumulation formula for the continuous- time interaction expansion (CT-INT)Rubtsov05 () and auxiliary field (CT-AUX)CTAUX () algorithms in Appendix D. Finally, in Appendix E, we derive the expression for the matrix that relates the coefficients of the Green’s function in the Legendre representation to its Matsubara frequency representation.

## I Single-particle Green’s function

### i.1 Legendre representation

We consider the single-particle imaginary-time Green’s function defined on the interval , where is the inverse temperature. Expanding in terms of Legendre polynomials defined on the interval , we have

 G(τ)= ∑l≥0√2l+1βPl(x(τ))Gl, (1) Gl= √2l+1∫β0dτPl(x(τ))G(τ). (2)

where and denote the coefficients of in the Legendre basis. The most important properties of the Legendre polynomials are summarized in Appendix A.

We note that a priori different orthogonal polynomial bases (e.g., Chebyshev instead of Legendre polynomials) may be used, and many of the conclusions in this paper would remain valid. The advantage of the Legendre polynomials is that the transformation between the Legendre representation and the Matsubara representation can be written in terms of a unitary matrix, since Legendre polynomials are orthogonal with respect to a scalar product that does not involve a weight function (see below and Appendix E). In this paper, therefore we restrict our discussion to the Legendre polynomials.

On general grounds, one can expect the Legendre representation of to be much more compact than the standard Matsubara representation: in order to perform a Fourier series expansion in terms of Matsubara frequencies, has to be anti-periodized for all , while the full information is already contained in the interval . As a result, the Green’s function contains discontinuities in that result in a slow decay at large frequencies (typically ). On the other hand, expanding , which is a smooth function of on the interval , in terms of Legendre polynomials yields coefficients that decay faster than the inverse of any power of (as shown in Appendix B). As a result, the information about a Green’s function can be saved in a very small storage volume. As we will show in Sec. II, this is particularly relevant when dealing with more complex objects such as the two-particle Green’s function, which depends on three frequencies.

CTQMC algorithms usually measure the Green’s function in one of the two following ways: (i) using a very fine grid for the interval or (ii) measuring the Fourier transform of the Green’s function on a finite set of Matsubara frequencies. Rubtsov05 (); Haule07 () We show in Appendix C explicitely for the CT-HYB Werner06 (); ctqmcrmp () algorithm, that one can also directly measure the coefficients during the Monte Carlo process (we expect our conclusions to hold for any continuous-time Monte Carlo algorithm).

As an illustration, we will focus on the Green’s function obtained by DMFT for the Hubbard model at half-filling described by the Hamiltonian

 H=−t∑⟨ij⟩σc†iσcjσ+U∑ini↑ni↓, (3)

where creates (annihilates) an electron with spin on the site of a Bethe lattice bethe (); Georges96 () and on the sum denotes nearest neighbors. In the following, quantities will be expressed in units of the hopping and we set the on-site Coulomb repulsion to and use the temperature . We solve the DMFT equations using the TRIQSTRIQS () toolkit and its implementation of the CT-HYB Werner06 (); ctqmcrmp () algorithm. In Fig. 1, we show the coefficients that we obtain. Note that coefficients for odd must be zero due to particle-hole symmetry. Indeed, the coefficients in our data for odd ’s all take on very small value, compatible with a vanishing value within their error bars. The even coefficients instead show a very fast decay, as discussed above. For , all coefficients eventually take values of the order of the statistical error bar.

Let us now discuss the specific issue of the statistical Monte Carlo noise. We observe that the high-order Legendre coefficients have a larger relative noise than small coefficients. On general grounds, we expect the coefficients of the exact Green’s function to continue to decrease faster than any power of to zero (cf. Appendix B). Hence, physical quantities computed from are likely to have a very weak dependence on the for large . A good approximation then is to truncate the expansion in Legendre polynomials at an order and set for . The choice for has to be such that the quantity of interest is accurately represented. On the other hand, if is too large, we would start to include coefficients that have increasingly large error bars compared to their value and this would eventually pollute the calculation. A systematic method is therefore to examine the physical quantity as a function of the cutoff . We expect that it will first reach a plateau where it is well converged. The existence of a plateau means that the contribution of higher-order coefficients is indeed negligible. For larger , the statistical noise in the will destabilize this plateau whose size will increase with the precision of the CTQMC computation. The existence of such a plateau provides a controlled way to determine the adequate value of . In the remaining paragraphs of this section, we will illustrate this phenomenon on different physical quantities by studying their dependence on .

### i.2 Imaginary-time Green’s function

It is instructive to analyze the effect of on the reconstructed imaginary-time Green’s function (using Eq. (1)). In Fig. 2, we show the evolution of at , , , and with the cutoff. It is apparent that these values very rapidly converge as a function of . We observe a well-defined and extended plateau. As the cutoff grows bigger, noise reappears in because of the comparatively large error bars in higher-order ’s.

In Fig. 3, the Green’s function is reconstructed on the full interval and compared to a direct measurement on a 1500-bin mesh. For , where the individual values of have not yet converged to their plateau (see Fig. 2), the resulting Green’s function is smooth but not compatible with the scattered direct measurements. For and , is smooth and nicely interpolates the scattered data. Moreover is virtually identical for both values of . This is expected because both of these values lie on the plateau. When is very large, i.e., of the order of the number of imaginary-time bins, the noise in eventually reappears and begins to resemble that of the direct measurement. We emphasize that all measurements have been performed within the same calculation and hence contain identical statistics. Hence the information in both measurements is identical up to the error committed by truncating the basis.

It is clear from this analysis, that the truncation of the Legendre basis acts as a noise filter. We note that no information is lost by the truncation: the high-order coefficients correspond to information on very fine details of the Green’s function, which cannot be resolved within a Monte Carlo calculation, as is obvious from the noisy .

### i.3 Matsubara Green’s function and high-frequency expansion

It is common to use the Fourier transform of to manipulate Green’s functions. This representation is, for example, convenient to compute the self-energy from Dyson’s equation or to compute correlation energies. In terms of , we can obtain the Matsubara Green’s function with

 G(iνn) =∑l≥0Gl√2l+1β∫β0dτeiνnτPl(x(τ)) =∑l≥0TnlGl. (4)

where the unitary transformation is shown in Appendix E to be

 Tnl=(−1)nil+1√2l+1jl((2n+1)π2) (5)

with denoting the spherical Bessel functions. Note that is independent of .

In Fig. 4, we display the Matsubara Green’s function as measured directly on the Matsubara axis and as computed from Eq. (4) with a fixed cutoff . The direct measurement of has been done within the same Monte Carlo simulation as the one used to compute the discussed above. It is clear from the plot that the truncation to has filtered the high-frequency noise, and that for large the Matsubara Green’s function has a smooth power-law decay. Let us emphasize here that the Matsubara Green’s function is obtained in an unbiased manner that does not involve any model-guided Fourier transform (see also Ref. gunnarsson, ).

We will now show that the coefficients that control this power-law decay can also be accurately computed. Let us consider the high-frequency expansion of

 G(iνn)=c1iνn+c2(iνn)2+c3(iνn)3+… (6)

Using the known high-frequency expansion of (cf. Appendix E),

 Tnl=t(1)liνnβ+t(2)l(iνnβ)2+t(3)l(iνnβ)3+…, (7)

one can directly relate the and the . Indeed, from (4), (6), and (7), it follows that

 cp=1βp∑l≥0t(p)lGl. (8)

The general expression of the coefficients is shown in (69). For the first three moments, we have the following expressions

 c1 =−∑l≥0,even2√2l+1βGl (9a) c2 =+∑l>0,odd2√2l+1β2Gll(l+1) (9b) c3 =−∑l≥0,even√2l+1β3Gl(l+2)(l+1)l(l−1). (9c)

Since , with the fast decay of the discussed above, we can expect a stable convergence of the as a function of . Note, however, that when increases, the coefficients grow, so we expect to need more and more Legendre coefficients to compute the series in practice.

The convergence of the moments is illustrated in Fig. 5. For the model we consider, the first moments are explicitly given by

 c1 =1c2=0 c3 =5c4=0.

We see that and smoothly converge to a plateau. For the higher moment , a larger number of Legendre coefficients is required. A plateau is reached but is (depending on the accuracy of the data accumulated in the QMC simulation) quickly destabilized when gets bigger and noisy are included in the calculation. This clearly shows that has to be chosen carefully to get sensitive . For larger cutoff the error in the moments grows rapidly. This shows that a large error on the high-frequency moments is committed when measuring in a basis in which it is not possible to filter the noise, i.e., the conventional imaginary-time or Matsubara representation.

Note that it is easy to incorporate a priori information on the moments . For example, in the model we consider above, we have . From (9a), we see that this is a linear constraint on the coefficients, which we can therefore enforce by projecting the Legendre coefficients onto the ()-dimensional hyperplane defined by the constraint (9a). A correction to impose, e.g., a particular is straightforwardly found to be

 Gl→Gl+(βc1−lmax∑l′=0t(1)l′Gl′)t(1)l∑l|t(1)l|2. (10)

This is easily generalized to other constraints.

### i.4 Energy

The accurate determination of the high-frequency coefficients is of central importance, since many quantities are computed from sums over all Matsubara frequencies involving . Because slowly decreases as to leading order, these sums are usually computed from the actual data up to a given Matsubara frequency and the remaining frequencies are summed up analytically from the knowledge of the . Thus, an incorrect determination of the leads to significant numerical errors. This is a particularly delicate issue when is measured directly on the Matsubara axis. In this case one usually needs to fit the noisy high-frequency data to infer the high-frequency moments. As discussed above, such a procedure is not required when using Legendre coefficients and the can be computed in a controlled manner. In the following, we illustrate this point in an actual energy calculation.

Based on an LDA+DMFT calculation for the compound SrVOamadon2008 (); aichhorn2009 () we compute the kinetic energy and the correlation energy ( denotes the number of lattice sites) resulting from the implementation and parameters of Ref. aichhorn2009, . These terms are contributions to the LDA+DMFT total energy Amadon06 () which depend explicitly on the results of the DMFT impurity solver.

The results are shown in Fig. 6. Here the parameter , against which these quantities are plotted, represents the number of Legendre coefficients used throughout the LDA+DMFT self-consistency. It is also the number of coefficients used to evaluate from the lattice Green’s function . Note that has been accumulated directly within the CTQMC simulation.

In agreement with an analysis of the convergence with respect to the number of Legendre coefficients similar to the ones shown in Figs. 2, 5 for an individual DMFT iteration, we find a plateau for both energies at . While the energy can be accurately computed within a single DMFT iteration, the error here mainly stems from the fluctuations between successive DMFT iterations. The plateau remains up to the largest values of . However, as gets larger, so do the error bars, due to the feedback of noise from the largest Legendre coefficients. Note that the error bars on the correlation energy, computed directly within the CTQMC algorithm, are of the same order of magnitude as those on the kinetic energy. The existence of a plateau implies that for a well-chosen cutoff , the energy can be computed in a controlled manner. We want to emphasize that such an approach is simpler and better controlled than delicate fitting procedures of high-frequency tails of the Green’s function on the Matsubara axis.

## Ii Two-particle Green’s function

### ii.1 Legendre representation for two-particle Green’s functions

The use of Legendre polynomials proves very useful when dealing with two-particle Green’s functions. We will show that it brings about improvements both from the perspective of storage size and convergence as a function of the truncation. The object one mainly deals with is the generalized susceptibility,

 ˜χσσ′ (τ12,τ34,τ14)=˜χσσ′(τ1−τ2,τ3−τ4,τ1−τ4)= ⟨Tc†σ(τ1)cσ(τ2)c†σ′(τ3)cσ′(τ4)⟩ −⟨Tc†σ(τ1)cσ(τ2)⟩⟨Tc†σ′(τ3)cσ′(τ4)⟩. (11)

Let us emphasize that is a function of three independent time-differences only. With the particular choice made above, is -antiperiodic in and , while it is -periodic in . Consequently, its Fourier transform is a function of two fermionic frequencies , , and one bosonic frequency .

We introduce a representation of in terms of the coefficients such that

 ˜χ(τ12,τ34,τ14)=∑l,l′≥0∑m∈Z√2l+1√2l′+1β3(−1)l′+1 Pl(x(τ12))Pl′(x(τ34))eiωmτ14˜χll′(iωm). (12)

In this mixed basis representation, the and dependence of is expanded in terms of Legendre polynomials, while the dependence is described through Fourier modes . The motivation behind this choice is that many equations involving generalized susceptibilities (like the Bethe-Salpeter equation) are diagonal in . The inverse of (12) reads

 ˜χll′(iωm)= ∭dτ12dτ34dτ14√2l+1√2l′+1(−1)l′+1 Pl(x(τ12))Pl′(x(τ34))e−iωmτ14˜χ(τ12,τ34,τ14). (13)

We show in Appendix C how the Legendre expansion coefficients of the one- and two-particle Green’s function (hence of ) can be measured directly within CT-HYB. With the above definition, the Fourier transform is easily found with

 ˜χ(iνn,iνn′,iωm)=∑l,l≥0Tnl˜χll′(iωm)T∗n′l′. (14)

was already defined in Eq. (4). Using the additional unitarity property of in Eq. (14) one can in general easily rewrite equations involving the Fourier coefficients in sole terms of the .

In the DMFT framework, the lattice susceptibility is obtained from Georges96 (); MaierReview ()

 [˜χlatt––––––––––]−1(iωm,q)= [˜χloc––––––––]−1(iωm) − [˜χ0loc––––––––]−1(iωm)+[˜χ0latt––––––––––]−1(iωm,q), (15)

where the double underline emphasizes that this is to be thought of as a matrix equation for the coefficients expressed either in in the Fourier representation or in in the mixed Legendre-Fourier representation. The bare susceptibilities are given by

 ˜χ0loc(iνn,iνn′,iωm)=−Gloc(iνn+iωm)Gloc(iνn)δn,n′, ˜χ0latt(iνn,iνn′,iωm,q) =−∑kGlattk+q(iνn+iωm)Glattk(iνn)δn,n′, (16)

where

 Glattk(iνn)=[iνn+μ−ϵk−Σloc(iνn)]−1, (17)

and , are the Green’s function and self-energy of the local DMFT impurity problem, respectively. The equivalent susceptibilities in the mixed Legendre-Fourier representation are simply obtained as the inverse of Eq. (14)

 ˜χ0ll′(iωm,q)=∑n,n′∈ZT∗nl˜χ0(iνn,iνn′,iωm,q)Tn′l′, (18)

where the high-frequency behavior of and can easily be considered in the frequency sums. Evaluation of lattice susceptibilities from can also directly be propagated to the mixed Legendre-Fourier representation, abolishing altogether the need to transform back to Fourier representation

 χ(i ωm,q)=1β2∑nn′∈Z˜χlatt(iνn,iνn′,iωm,q) =1β2∑ll′≥0(−1)l+l′√2l+1√2l′+1˜χlatt,ll′(iωm,q). (19)

Note that can be written as the sum of a free two-particle propagation Eq. (16) (bubble part) and a connected rest (vertex part). These two terms can be separately summed in Eq. (19).

The present mixed basis representation has been successfully used in a recent investigation of static finite-temperature lattice charge and magnetic susceptibilities for the NaCoO system at intermediate-to-larger doping Boehnke10 () A first example for the dynamical, i.e., finite-frequency, case will be discussed in Sec. II.3.

### ii.2 Antiferromagnetic susceptibility of the three-dimensional Hubbard model

In order to benchmark our approach, we investigate the antiferromagnetic susceptibility of the half-filled Hubbard model (3) on a cubic lattice within the DMFT framework. All quantities are again expressed in units of the hopping and with and . This temperature is sufficiently close to the DMFT Néel temperature to yield a dominant vertex part, while still having a non-negligible bubble contribution.

We compute the susceptibility of the DMFT impurity problem using the CT-HYB algorithm. In Fig. 7 we compare the mixed Legendre-Fourier coefficients to the Fourier coefficients . For clarity, we focus on the first bosonic frequency . We observe that the have a very fast decay except in the direction. This contrasts with the behavior of which exhibits slower decay in the three major directions , and .

The generalized susceptibility in -differences (11) has discontinuities along the planes and as well as non-analyticities (kinks) for and . These planes induce corresponding slow decay in the Fourier representation (14).gunnarsson () When it comes to the mixed Legendre-Fourier representation (13) however, the planes and are on the border of the imaginary-time region being expanded in this basis, which renders the coefficients insensitive toward these.

Computing lattice susceptibilities from Eq. (15), it is necessarily required to truncate the matrices. This leads to difficulties when computing the susceptibility from the Fourier coefficients . As we can see from Fig. 7, the Fourier coefficients have a slow decay along three directions. The inversion of is delicate because many coefficients are involved even for large . One needs to use a very large cutoff to obtain a precise result. Alternatively, one can try to separate the high- and low-frequency parts of the equation and replace the susceptibilities with their asymptotic form at high frequency (see Ref. kunes, ). While is is effectively possible to treat larger matrices, it is still required to impose a cutoff on the high-frequency part for the numerical computations.

In the mixed Legendre-Fourier representation, the situation is different. Only the coefficients along the diagonal decay slowly. In the inversion of the matrix, the elements on the diagonal for large are essentially recomputed from themselves. One can expect that there will be a lot less mixing and thus a much faster convergence as a function of the truncation.

In Fig. 8, we display the vertex part of the generalized lattice susceptibility obtained from Eq. (15) in both representations. In both cases, we see that the diagonal part quickly becomes very small. In other words, the diagonal of the lattice susceptibility is essentially given by the bubble part . However, while essentially all the information is condensed close to in the mixed Legendre-Fourier representation, the Fourier coefficients still have a slow decay along the directions given by and . From this figure one can speculate that a quantity computed from the Legendre-Fourier coefficients will converge rapidly as a function of a cutoff . However, we need to make sure that the coefficients close to are not affected much by the truncation.

In order to assess the validity of these speculations we compute the static antiferromagnetic () susceptibility as a function of the cutoff in both representations. It is obtained from Eq. (19) using the magnetic susceptibility .

Since the diagonal of the lattice susceptibility is essentially given by the bubble (see Fig. 8), the sums above are performed in two steps. The vertex part shown in Fig. 8 is summed up to the chosen cutoff, while the bubble part is summed over all frequencies with the knowledge of its high-frequency behavior. The result is shown in Fig. 9. It reveals a major benefit of the Legendre representation: the susceptibility converges much faster as a function of the cutoff. The static susceptibility is essentially converged at . This corroborates the idea that the small- part of is only weakly dependent on the further diagonal elements of .

### ii.3 Dynamical susceptibility of the two-dimensional Hubbard model

As a final benchmark, we demonstrate that our method is not restricted to the static case. To this end, we show the momentum resolved dynamical magnetic susceptibility for a DMFT calculation for the half-filled two-dimensional (2D) square lattice Hubbard model in Fig. 10. We have chosen an on-site interaction and temperature , which is slightly above the DMFT Néel temperature. The susceptibility was computed from the Legendre representation according to Eq. (19) using Legendre coefficients, which was sufficient for all bosonic frequencies. In general, for higher bosonic frequencies more Legendre coefficients are needed to represent the vertex part of the generalized magnetic lattice susceptibility. However, no additional structure appears in the high , region. We then analytically continued the data using Padé approximants. pade () The figure shows the typical magnon spectrum Preuss97 (); Hochkeppel08 (); ldfa () reminiscent of a spin wave in this paramagnetic state with strongly enhanced weight at the antiferromagnetic wave vector due to the proximity of the mean-field antiferromagnetic instability.

## Iii Conclusion

In this paper, we have studied the representation of imaginary-time Green’s functions in terms of a Legendre orthogonal polynomial basis. We have shown that CTQMC can directly accumulate the Green’s function in this basis. This representation has several advantages over the standard Matsubara frequency representation: (i) It is much more compact, i.e., coefficients decay much faster; this is particularly interesting for storing and manipulating the two-particle Green’s functions. Moreover, two-particle response functions can be computed directly in the Legendre representation, without the need to transform back to the Matsubara representation. In particular, the matrix manipulations required for the solution of the Bethe-Salpeter equations can be performed in this basis. We have shown that this greatly enhances the accuracy of the calculations, since in contrast to the Matsubara representation the error due to the truncation of the matrices becomes negligible. (ii) The Monte Carlo noise is mainly concentrated in the higher Legendre coefficients, the contribution of which is usually very small; this allows us to develop a systematic method to filter out noise in physical quantities and to obtain more accurate values for, e.g., the correlation energy in LDA+DMFT computations.

###### Acknowledgements.
L.B. and F.L. thank A.I. Lichtenstein for helpful discussions. M.F. and O.P. thank M. Aichhorn for helpful discussions and providing the parameters and data from Ref. aichhorn2009, . O.P. thanks J.M. Normand for pointing out Ref. BatemanBook, . We thank P. Werner for a careful reading of the manuscript. Calculations were performed with the TRIQSTRIQS () hier project using HPC resources from The North-German Supercomputing Alliance (HLRN) and from GENCI-CCRT (Grant No. 2011-t2011056112). TRIQS uses some libraries of the ALPS ALPS20 () project. This work was supported by the DFG Research Unit FOR 1346.

## Appendix A Some Properties of the Legendre Polynomials

In this appendix, we summarize for convenience some basic properties of the Legendre polynomials. Further references can be found in Refs. GradshteynBook, ; WhittakerWatsonBook, ; BatemanBook, . We use the standardized polynomial defined on through the recursive relation

 (l+1)Pl+1(x) =(2l+1)xPl(x)−lPl−1(x) (20) P0(x) =1,P1(x)=x (21)

are orthogonal and their normalization is given by

 ∫1−1dxPk(x)Pl(x)=22l+1δkl (22)

The are bounded on the segment by WhittakerWatsonBook ()

 |Pl(x)|≤1 (23)

with the special points

 Pl(±1)=(±1)l. (24)

The primitive of that vanishes at is (cf. Ref. BatemanBook, , Vol II, section 10.10)

 ∫x−1dyPl(y)=Pl+1(x)−Pl−1(x)2l+1,l≥1 (25)

By orthogonality or (24), it also vanishes at . The Fourier transform of the Legendre polynomial restricted to the segment is given by formula 7.243.5 of Ref. GradshteynBook,

 ∫1−1eiaxPl(x)dx =il√2πaJl+12(a) =2iljl(a), (26)

where denotes the Bessel function and denotes the spherical Bessel functions.

## Appendix B Fast Decay of the Legendre Coefficients

Let us consider a function smooth on the segment (i.e. to be precise , indefinitely differentiable), and antiperiodic, like a Green’s function. In this appendix, we show that its Legendre coefficients decay faster than any power law contrary to its standard Fourier expansion coefficients which decay as power laws determined by the discontinuities of the function and its derivatives.

Let us start by reminding the asymptotics of the standard Fourier expansion coefficients on fermionic Matsubara frequencies. These coefficients are given by

 ^g(iνn) =∫β0dτ g(τ)eiνnτ (27) =g(τ)eiνnτ∣∣β0iνn−∫β0dτ g′(τ)eiνnτiνn (28)

The coefficients vanish for , and applying the same result to , one obtains

 ^g(iνn)=−g(β−)+g(0+)iνn+O(1ν2n) (29)

Let us now turn to the Legendre expansion. Using the same rescaling as before, we can consider for simplicity a function smooth on . We can proceed in a similar way using the primitive of the Legendre polynomial (which is also given by a simple formula, (25)). For , we have

 fl√2l+1= ∫1−1dx f(x)Pl(x) = f(x)(∫x−1dy Pl(y))∣∣∣1−1− ∫1−1dx f′(x)(∫x−1dy Pl(y)) = −∫1−1dx f′(x)Pl+1(x)−Pl−1(x)2l+1 (30)

The crucial difference with the Fourier case is that, for , the boundary terms always cancel, whatever the function due to the orthogonality property of the polynomials (it can also be checked directly from (24)). So we are left with just the integral term. Since the Legendre coefficients of vanish at large (by applying the previous formula to ), we get instead of (29)

 fl√2l+1=o(1l) (31)

In both cases, the reasoning can be reproduced recursively, by further differentiating the function, as long as no singularity are encountered. In the Fourier case, it produces the well-known high-frequency expansion in terms of the discontinuity of the function and its derivatives. In the Legendre case, we find that the coefficients are as soon as is times differentiable. Hence if the function is smooth on , the coefficients decays asymptotically faster than any power law.

The only point that remains to be checked is that indeed is smooth on . It is clear from its spectral representation

 G(τ)=−∫∞−∞dνe−τν1+e−βνA(ν) (32)

if we admit that the spectral function has compact support, by differentiating under the integral.

Finally, while this simple result of “fast decay” is enough for our purposes in this paper, it is possible to get much more refined statements on the asymptotics of the Legendre coefficients of the function , in particular when it has some analyticity properties. For a detailed discussion of these issues, and in particular of the conditions needed to get the generic exponential decay of the coefficients, we refer to Ref. BoydBook, .

## Appendix C Direct Accumulation of the Legendre Coefficients for the Ct-Hyb Algorithm

In this appendix, we describe how to compute directly the Legendre expansion of the one-particle and the two-particle Green’s function.

### c.1 The accumulation formulas in CT-HYB

For completeness, let us first recall the accumulation formula for the one-particle and the two-particle Green’s functions in the CT-HYB algorithm Werner06 (); Haule07 (); WernerPRBLong (); ctqmcrmp (), which sums the perturbation theory in the hybridization function on the Matsubara axis. While these formulas have appeared previously in the literature, this simple functional derivation emphasizes the “Wick”-like form of the high-order correlation function.

The partition function of the impurity model reads

 Z=∫Dc†Dcexp(−Seff) (33)

where the effective action has the form

 Seff=−∬β0dτdτ′∑A,Bc†A(τ)G−10,AB(τ,τ′)cB(τ′) +∫β0dτHint({c†A(τ),cA(τ)}) (34) G−10AB(iνn)=(iνn+μ)δAB−h0AB−ΔAB(iνn), (35)

To simplify the notations, we use here a generic index . In the case where there are symmetries, like the spin symmetry in the standard DMFT problem, the Green’s functions are block diagonal. For example, the generic index can be , where is an orbital or site index, and spin index is the block index.

The partition function is expanded in powers of the hybridization as

 Z=∑n≥0∫n∏i=1dτidτ′i∑λi,λ′iw(n,{λj,λ′j,τj,τ′j}) (36) w(n,{λj,λ′j,τj,τ′j})≡1n!2det1≤i,j≤n[Δλi,λ′j(τi−τ′j)]× Tr(Te−βHlocn∏i=1c†λi(τi)cλ′i(τ′i)), (37)

where is time ordering and is the local Hamiltonian Werner06 (); Haule07 (); WernerPRBLong (); ctqmcrmp (). are the weights of the Quantum Monte Carlo Markov chain. Introducing the short notation for the QMC configuration, the partition function and the average of any function over the configuration space (denoted by angular bracket in this section) are given by

 Z =∑Cw(C) (38) ⟨f(C)⟩ =1Z∑Cw(C)f(C) (39)

The one-particle and two-particle Green’s functions are obtained as functional derivatives of with respect to the hybridization function, as

 GAB(τ1,τ2) =−1Z∂Z∂ΔBA(τ2,τ1) (40a) G(4)ABCD(τ1,τ2,τ3,τ4) =1Z∂2Z∂ΔBA(τ2,τ1)∂ΔDC(τ4,τ3) (40b)

In order to use the expansion of , we need to compute the derivative of a determinant with respect to its elements. Let us consider a general matrix , its inverse and use Grassman integral representation

 det¯Δ=∫∏idηid¯ηie∑ij¯ηi¯Δijηj (41)

Using the Wick theorem, we have

 ∂det¯Δ∂¯Δba =∫∏idηid¯ηi(¯ηbηa)e∑ij¯ηi¯Δijηj =det¯Δ×¯Mab (42a) ∂2det¯Δ∂¯Δba∂¯Δdc =∫∏idηid¯ηi(¯ηbηa¯ηdηc)e∑ij¯ηi¯Δijηj =det¯Δ(¯Mab¯Mcd−¯Mad¯Mcb) (42b)

Let us now apply (42) by introducing for each configuration the matrix of size given by

 ^Δ(C)ij≡Δλi,λ′j(τi−τ′j) (43)

and its inverse . We obtain

 ∂w(C)∂ΔBA(τ2,τ1) =w(C)det^Δ(C)n∑α,β=1∂det^Δ(C)∂^Δ(C)βα∂^Δ(C)βα∂ΔBA(τ2,τ1) =w(C)n∑α,β=1MCαβ∂^Δ(C)βα∂ΔBA(τ2,τ1) (44a) and ∂2w(C)∂ΔBA(τ2,τ1)∂ΔDC(τ4,τ3)=w(C)det^Δ(C)×n∑αβγδ=1∂2det^Δ(C)∂^Δ(C)βα∂^Δ(C)δγ∂^Δ(C)βα∂ΔBA(τ2,τ1)∂^Δ(C)δγ∂ΔDC(τ4,τ3)