Improved Estimators for the Self-Energy and Vertex Function in Hybridization Expansion Continuous-Time Quantum Monte Carlo Simulations

Improved Estimators for the Self-Energy and Vertex Function in Hybridization Expansion Continuous-Time Quantum Monte Carlo Simulations

Hartmut Hafermann Centre de Physique Théorique, École Polytechnique, CNRS, 91128 Palaiseau Cedex, France    Kelly R. Patton Department of Physics and Astronomy, Louisiana State University, Baton Rouge, Louisiana 70803    Philipp Werner Theoretische Physik, ETH Zurich, 8093 Zürich, Switzerland
July 3, 2019

We propose efficient measurement procedures for the self-energy and vertex function of the Anderson impurity model within the hybridization expansion continuous-time quantum Monte Carlo algorithm. The method is based on the measurement of higher-order correlation functions related to the quantities being sought through the equation of motion, a technique previously introduced in the NRG context. For the case of interactions of density-density type, the additional correlators can be obtained at essentially no additional computational cost. In combination with a recently introduced method for filtering the Monte Carlo noise using a representation in terms of orthogonal polynomials, we obtain data with unprecedented accuracy. This leads to an enhanced stability in analytical continuations of the self-energy or in two-particle based theories such as the dual fermion approach. As an illustration of the method we reexamine the previously reported spin-freezing and high-spin to low-spin transitions in a two-orbital model with density-density interactions. In both cases, the vertex function undergoes significant changes, which suggests significant corrections to the dynamical mean-field solutions in dual fermion calculations.


I Introduction

Continuous-time quantum Monte Carlo solvers (for a recent review, see Ref. Gull et al., 2011) have become an important tool for the study of the Anderson impurity model (AIM) and its multi-orbital generalizations, due to their accuracy, efficiency and ability to treat arbitrary interaction terms. While the AIM plays a fundamental role in various areas of condensed matter physics, the continuous-time solvers have been developed and primarily applied in the context of dynamical mean-field theory (DMFT),Georges et al. (1996) which maps correlated lattice models to an appropriately defined quantum impurity model.

Impurity models with multiple orbitals are important for two different reasons: (i) The combination of density functional theory and DMFT Kotliar et al. (2006) provides a formalism for the calculation and prediction of properties of strongly correlated materials. The description of materials with open d- or f-shells requires the solution of impurity models with up to five or seven correlated orbitals, respectively. (ii) In the context of cluster extensions of DMFT,Maier et al. (2005) the lattice is mapped onto a cluster of impurities in order to account for spatial correlations, and it is desirable to solve clusters with as many sites as possible, in order to be able to infer reliable predictions for the infinite system. For cluster DMFT calculations of simple models, the interaction expansion or weak-coupling algorithm, which is based on an expansion of the partition function in the interaction (henceforth referred to as CT-INTRubtsov et al. (2005)) and the related continuous-time auxiliary field algorithm (CT-AUXGull et al. (2008)) have proven useful.Werner et al. (2009); Gull et al. (2009) Here the number of interaction terms and hence the perturbation order grow linearly with the number of cluster sites. In the context of ab initio calculations of correlated materials, however, the number of interaction terms grows at least as the square of the number of orbitals and hence the hybridization expansion algorithm Werner et al. (2006); Werner and Millis (2006); Haule (2007) (abbreviated CT-HYB) is the method of choice.

For the latter, the calculation of the one-electron self-energy has proven problematic. The calculation from Dyson’s equation, where it is evaluated as the difference between the inverses of two Green’s functions, is highly susceptible to noise in the numerical data. In contrast to CT-INT, the Green’s function in CT-HYB is not measured as a correction to a known function (the noninteracting Green’s function ), which leads to large noise in the intermediate to high-frequency region. Except for the low-frequency region, better statistics is needed for the calculation of the self-energy in CT-HYB to obtain results of comparable accuracy as in CT-INT.Gull et al. (2007)

Similar problems arise in the calculation of the impurity vertex function. While the vertex allows one to access response functions of a system, interest in this quantity has recently grown in particular due to the advent of novel diagrammatic extensions of the dynamical mean-field theory.Kusunose (2006); Toschi et al. (2007); Slezak et al. (2009); Rubtsov et al. (2008); Hafermann et al. (2009) In these approaches, spatial correlations beyond DMFT are included through two-particle field theories, which involve the two-particle irreducible (in the dynamical vertex approximation Toschi et al. (2007)) or reducible (in the dual fermion approach Rubtsov et al. (2008); Hafermann et al. (2009)) vertex function. Thus far these approaches mainly relied on exact diagonalization (ED) or CT-INT for the calculation of the impurity vertex functions. Only few applications employing the CT-HYB algorithm for measuring the two-particle function have been reported because of the aforementioned problems. The measurement of frequency dependent two-particle functions is particularly challenging for multiorbital models. For example, in a recent letter by Park et al., this problem was circumvented by neglecting the frequency dependence of the irreducible vertex function and hence setting it to a constant.Park et al. (2011)

On the other hand, the CT-HYB approach with its low perturbation order and favorable sign statisticsGull et al. (2007) appears to be the most suitable, in principle, for measuring vertex functions in realistic multiorbital models.

In this paper, we propose an efficient method for calculating the self-energy and vertex function which eliminates the noise problems. Through the equation of motion, the self-energy and vertex function are related to higher-order correlation functions, which can be measured in CT-HYB. In combination with a recently developed method to eliminate the Monte Carlo noise through a representation in terms of orthogonal polynomials,Boehnke et al. (2011) the present approach yields considerably more accurate results compared to the standard measurements.

Ii Method

ii.1 Model

Figure 1: (Color online) Typical segment configuration in the hybridization expansion continuous-time quantum Monte Carlo simulation of the one-orbital AIM. Each segment marks an imaginary-time interval in which an electron of given spin resides on the impurity. The overlap between the up- and down-spin segments (blue or gray shaded area) determines the interaction energy.

We consider the single-site, multiorbital AIM with interactions of density-density type, described by the Hamiltonian


where latin subscripts denote a ‘flavor’ index, i.e. a combined index for spin- and orbital degrees of freedom. The operators create (destroy) an electron with flavor on the impurity site, while creates (destroys) a conduction band electron in band with momentum ( is the corresponding dispersion). The impurity levels are denoted by and the matrix contains the various interaction parameters for the interactions of density-density type. Electrons with band-flavor index and momentum are allowed to couple to electrons with any other flavor index as described by the hybridization matrix . Integrating out the noninteracting conduction band electrons gives rise to the hybridization function


Here we use an implementation based on the segment picture of the hybridization expansion algorithm.Werner et al. (2006) The segment picture applies whenever operators of a given flavor appear in alternating order, which is always the case for an interaction of density-density type.

In this case the trace over the sequence of impurity creation and annihilation operators (which defines the Monte Carlo configuration) does not have to be computed explicitly. In order to evaluate the ratio of traces one instead only needs to compute the length of the segments for the different flavors and their overlaps. This yields a very efficient algorithm for multiorbital problems with density-density interaction, which — as long as the determinant calculation dominates the computational effort — scales linearly in the number of flavors.Gull et al. (2011) A possible segment configuration for the one-orbital AIM is illustrated in Fig. 1. The segments represent imaginary-time intervals in which an electron of given flavor (spin) resides on the impurity. Overlapping segments correspond to time intervals in which the impurity is doubly occupied. As we will see below, the additional higher-order correlation functions which arise in the expressions for the self-energy and vertex function in the improved estimator can be measured at essentially no additional computational cost in this segment representation. The generalization to the case of interactions of non-density-density type, such as spin-flip and pair hopping terms, is in principle straightforward. The measurement of the required correlation functions in the hybridization expansion algorithm, however, becomes more involved.

In the following subsection, we first introduce and illustrate the procedure for the self-energy and consider the generalization to the vertex function in a second separate subsection.

ii.2 Self-energy

Figure 2: (Color online) The ways of connecting hybridization lines in a configuration with 3 segments for a timeline of a given flavor. The sign of each diagram is indicated. For non-diagonal hybridization, the lines also connect segments with different flavors.

The asymptotic tail of the self-energy can be obtained from a high-frequency expansion, which requires measurements of single- and two-particle equal-time correlators (see e.g. Ref. Wang et al., 2011 and references therein). This eliminates the noise at high frequencies. However, the choice of the cutoff value cannot easily be automatized, and as we will show (see, e.g., the inset of Fig. 6), the noise problem of the standard measurement is significant even at frequencies for which the self-energy clearly has not yet reached its asymptotic behavior.

The method presented here yields accurate results over the entire frequency range and considerably reduces the noise at intermediate to high frequencies. An alternative approach to reduce the Monte Carlo noise is to measure observables in an orthogonal polynomial representation.Boehnke et al. (2011) This corresponds to an advantageous change of basis, which yields a compact representation of observables and allows one to filter the Monte Carlo noise, without any loss of information. We note that such a procedure does not reduce the required statistics and hence for optimal results the two approaches should be combined.

The computation of the self-energy from higher-order correlation functions of the impurity model has been introduced previously in the context of numerical renormalization group (NRG) Bulla et al. (1998) calculations. Here we provide the expression for the multiorbital case and demonstrate the usefulness of this scheme for the CT-HYB algorithm. An expression for the self-energy in terms of a two-particle correlator is obtained from the equation of motion for the Green’s function. One may write the equation of motion in terms of a derivative with respect to either the first or second of its arguments. For reasons outlined in the appendix, we choose the second option. Here we only present a brief outline of the derivation. Details can be found in Appendix A.

The time derivative of the Green’s function


is given by


where is the usual time-ordering operator. Its equation of motion follows from the one for the operator taken in the Heisenberg representation:


We introduce the following correlation functions:


together with their Fourier transforms and . With an application in the context of DMFT in mind, we furthermore introduce the noninteracting Green’s function of the impurity model,


where is the hybridization function. Evaluating the commutator in Eq. (5) with the Hamiltonian (1) yields


so that Eq. (4) in Fourier space can be written


The function in turn can be eliminated by expressing it in terms of the impurity Green’s function through its equation of motion. In Fourier space, it reads


Inserting this into (10) and using the expression for the hybridization function (2), the hybridization terms cancel. Comparing the resulting expression with Dyson’s equation (see Fig. 3),


we see that the impurity self-energy can be expressed in the following form:


which for a single orbital model reduces to the result given in Ref. Bulla et al., 1998. This equation relates the self-energy to two measurable quantities, the impurity Green’s function and an additional correlation function .

Figure 3: Diagrammatical representation of the Green’s function in terms the self energy , Eq. (12). Thick lines denote fully dressed propagators, thin lines denote bare propagators.

In order to show how correlation functions are measured in the CT-HYB, we recall that in this algorithm, one samples over configurations whose weight is given by a determinant of hybridization functions times a trace over a sequence of operators. In the segment picture, a configuration is visualized by a collection of segments on the timeline from to for each flavor ( is the inverse temperature). A typical configuration for a single-orbital model with Hubbard interaction is depicted in Fig. 1. The start of a segment (open circles) is associated with an impurity creation operator, while an impurity annihilation operator is associated with the segment endpoint (closed circles). A segment hence marks the imaginary time interval in which the impurity is occupied by an electron of a given flavor. The segments are connected by hybridization lines in all possible ways, which is the interpretation of the determinant (see Fig. 2). For the case of non-diagonal hybridization considered here, the segments are also connected by hybridization lines linking different flavors (not shown).

Figure 4: (Color online) Self-energy as a function of Matsubara frequencies for the half-filled Hubbard model on the Bethe lattice (bandwidth ) as calculated from Dyson’s equation to illustrate the noise problem. The parameters are the same as in Ref. Gull et al., 2007, and . The results were obtained by measuring the Green’s function directly on Matsubara frequencies and in imaginary time using 1500 bins and subsequent Fourier transform for the latter. The high-frequency tail is shown for comparison.

A contribution to the Green’s function of a particular Monte Carlo configuration is obtained by cutting all hybridization lines connected to a given pair of a creation and an annihilation operator, leaving a configuration with two unconnected operators. This corresponds to removing row and column from the matrix of hybridization functions . Denoting the resulting matrix as , the contribution of the particular configuration is essentially the ratio between the matrices after and before removing the hybridization functions, i.e. . This ratio is precisely the -element of the inverse of the matrix of hybridization functions, denoted by . Hence the Green’s function defined on the interval from 0 to is measured according toWerner et al. (2006); Werner and Millis (2006)




The only difference in the measurement for the function is the presence of the additional density operator. Hence the measurement formula reads


In the segment picture, the matrix element (one or zero) of this operator is simply determined by examining whether or not a segment is present in the timeline for flavor at time . Hence this quantity can be measured at essentially no additional computational cost. Note that this function according to (13) only contributes if is different from the index (and ) and therefore is never evaluated at the position of the creator of flavor at time .

Figure 5: (Color online) CT-HYB data for the self-energy computed from Eq. (13) for the same model and parameters and measured within the same simulation as the results in Fig. 4. The Green’s function and the two-particle correlator have been measured on the Matsubara axis. The line with closed symbols (blue) shows the result from an independent DMFT calculation using the interaction expansion algorithm (CT-INT), where the self energy has been obtained via Dyson’s equation with the Green’s function measured on the Matsubara axis. The noise problem at intermediate and high frequencies does not exist in the CT-INT and is resolved in CT-HYB using the improved estimator. The inset shows a blowup of the intermediate frequency range.

It is convenient to accumulate


directly instead of the individual quantities . This is then analogous to the measurement of times in the CT-AUX algorithm.Gull et al. (2008) Note that these correlators can be measured in any basis by appropriately transforming the measurement rules, e.g. by taking the Fourier transform to measure directly the Matsubara coefficients.

If the interaction is of non-density-density type, the segment picture has to be abandoned. In this case, the equation of motion generates additional terms, which require the measurement of correlation functions of the form (see appendix C). In addition to computing the ratio of determinants, the measurement of such functions requires the evaluation of the ratio of traces over the operators in the particular configuration with and without the operators inserted at the corresponding time.

Figure 6: (Color online) Same as in Fig. 5, albeit with the Green’s function and the two-particle correlator measured in the Legendre polynomial basis and transformed back to Matsubara frequencies. This allows to efficiently filter the Monte Carlo noise. The inset compares the result at intermediate frequencies to that of Fig. 4 obtained from the Fourier transformed imaginary time measurement (data are measured within the same simulation).

The major benefit of rewriting the self-energy in the form (13) instead of computing it from Dyson’s equation


is that it is expressed in terms of a ratio of two measured quantities instead of a difference of two functions. Forming the difference between an exactly known and an approximate quantity is susceptible to numerical errors, since, as already noted in Ref. Bulla et al., 1998, for a ratio only the relative error propagates, while in a difference, the absolute error (here the absolute error of ) propagates. In DMFT, is computed from the self-energy of the previous iteration and is not known exactly. Error cancellation however cannot be expected since and are computed in two independent Monte Carlo runs. Since decays as it is clear that computing the self-energy according to Eq. (18) will lead to large errors in particular at intermediate to large frequencies.

In the following we show results illustrating the advantages of the improved measurement for the CT-HYB. We concentrate on the DMFT solution of the Hubbard model on the Bethe lattice with bandwidth . The parameters are the same as in Ref. Gull et al., 2007. In order to ensure comparability, we have performed the various measurements within the same simulation, i.e. all quantities have been averaged over the identical sequence of Monte Carlo configurations (including the results of Fig. 4).

Figure 4 shows the self-energy determined in the standard way, i.e. by measuring the Green’s function and calculating the self-energy using Dyson’s equation. The Monte Carlo noise is clearly visible for intermediate to large frequencies, as anticipated. It is important to note that there is also considerable noise in the region where the self-energy has clearly not yet reached its asymptotic behavior (cf. also inset of Fig. 6 below). While the results are similar for the frequency and imaginary time measurement at small to intermediate frequencies, the Monte Carlo error steadily grows with the frequency for the measurement in Matsubara frequencies. The noise sets in much below the Nyquist frequency (of approximately ). In the imaginary-time measurement the noise is suppressed above the Nyquist frequency and the result smoothly approaches the tail, which is enforced in the calculation of the Fourier transform.

Figure 7: (Color online) Comparison of Monte Carlo data for the imaginary part of the self-energy (obtained with the same method as in Fig. 6) with the exact result for a correlated site coupled to a single bath level. Some points of the Monte Carlo result are skipped to improve visualization. The result computed from Dyson’s equation (with the Green’s function measured on the Matsubara axis within the same simulation) illustrates the substantial improvement in accuracy.

Figure 5 clearly shows the advantage of the improved measurement: The error is considerably reduced over the entire frequency range. The number of measurements required for a given accuracy at intermediate and large frequencies is hence greatly reduced. In the same figure we compare the improved measurement of CT-HYB with a result obtained from a comparable run of the CT-INT using Dyson’s equation, with the Green’s function also measured on the Matsubara axis. We note that this is meant as an illustration of the convergence properties at high frequency. Here we do not attempt the intricate task of separating the various effects that influence the performance and accuracy of these complementary algorithms, which furthermore scale very differently with the number of flavors.Gull et al. (2011) Comparable refers to the fact that for both runs we have used similar simulation times that yield converged results in practice. A detailed performance comparison of the algorithms was presented in Ref. Gull et al., 2007.

The result illustrates that the noise problem does not exist in the CT-INT algorithm. The reason is that in CT-INT, the Green’s function is measured as an correction to the bare Green’s function . As a consequence, the evaluation of the self-energy is considerably more stable (see below). We find no advantage of calculating the self-energy according to Eq. (13) in CT-INT, at least for the model considered here.

We note that in principle one may evaluate the self-energy without explicitly measuring the Green’s function . Combining Eqs. (12) and (13), one sees that the Green’s function may be written in the form


where the matrix is defined as


The self-energy is obtained by multiplying (17) by the inverse of (19). We find that this procedure yields results of similar but somewhat worse quality than the ones shown in Fig. 5, so that it does not reduce the effective computation time for given accuracy.

In Fig. 6 we show results obtained by combining the improved estimator with an efficient method to suppress the Monte Carlo noise. The latter is based on a representation in terms of Legendre orthogonal polynomials. The transformation to the Legendre basis is exact, as in the Fourier case. Both and have been measured directly in the Legendre basis and transformed back exactly to Matsubara frequencies after the simulation. Appropriately choosing the cutoff in this basis allows to eliminate the Monte Carlo noise without loosing physical information (for details see Ref. Boehnke et al., 2011). This method yields very accurate, noise-free results over the entire frequency range and captures the high-frequency tail correctly. We expect that Monte Carlo data measured in this way will allow a more stable analytical continuation.

In order to show that the proposed method not only yields smooth, but indeed correct results, we compare the Monte Carlo data to an exact result. A suitable test case is a correlated impurity site coupled to a single bath level at for which the Hamiltonian (1) reduces to


and which is trivially solved by exact diagonalization. The hybridization function for this case is . Although the hybridization function has no structure in this case, this is not a trivial problem for the Monte Carlo approach, since all diagrams in principle have to be sampled to obtain the exact solution. The result at half-filling is shown in Fig. 7 for and temperature . To improve the visualization, we have skipped some points from the calculated result. The curves lie on top of each other. That this is indeed not trivial can be seen from the result computed from Dyson’s equation using a Green’s function measurement on the Matsubara axis, which has again been accumulated within the same simulation.

ii.3 Vertex function

Figure 8: Diagrammatical representation of the two-particle Green’s function in terms of its connected part and its disconnected part and definition of the vertex function , Eqs. (23)-(25). The lines with arrows denote fully dressed propagators.
Figure 9: (Color online) Real part of the spin-up-up component of the connected part of the two-particle Green’s function for a correlated impurity coupled to a single bath-level for for fixed and bosonic frequencies (circles), (triangles) and (squares). The parameters are otherwise the same as in Fig. 7. Solid lines show the exact diagonalization data and open symbols the results from the improved Monte Carlo measurement (measured on the Matsubara axis). The usual Matsubara axis measurement based on Eqs. (22)-(25) is shown by crosses. The two results are hardly distinguishable on this scale.

To set the stage for the discussion of the vertex function, we define the two-particle Green’s function as


Its Fourier transform is . With the disconnected part of the two-particle Green’s function,


we obtain the connected part (cf. Fig. 8):


in terms of which the two-particle vertex function is defined as

Figure 10: (Color online) Real part of the spin up-up component of the reducible vertex function computed from the connected part shown in Fig. 9. The standard computation (Eqs. (22)-(25), crosses connected by dashed lines) is clearly numerically unstable for large values of , while the improved estimators yield the proper asymptotic behavior. The improved estimator computed from Eqs. (25), (26) (closed symbols) and from Eqs. (25), (31) (open symbols) yield similar results. They also yield more accurate results even before the vertex approaches its asymptotic behavior and for small frequencies (clearly visible for ).

The connected part of the two-particle Green’s function is the difference between a two-particle quantity and a product of two single-particle Green’s functions. The calculation of this difference is susceptible to numerical errors. Since single- and two-particle quantities have vastly different statistical errors, one cannot expect these to cancel. The error in the vertex function itself is further amplified by multiplying the connected part with inverse Green’s functions. This leads to deviations in particular when is small, i.e. for large frequencies or in an insulator, where extrapolates to zero for .

Hence we seek to express the connected part in a similar fashion as was done for the self-energy using the equation of motion. As shown in appendix B, the connected part of the two-particle Green’s function can be written in the alternative form


A priori it is not clear that computing the connected part according to this expression has any advantage compared to Eq. (24), since it still has the form of a difference of products of correlation functions. However we find that it indeed yields more accurate results.

In addition to the single- and two-particle Green’s function this expression involves the correlation functions and the Fourier transform of


In order to see how this function is measured, recall that in CT-HYB the two-particle Green’s function with its arguments varying in the interval from 0 to is measured according toGull et al. (2011)


where is defined as before (15) and . Due to time translation invariance, the two-particle Green’s function needs to be measured as a function of three independent time differences only. These have been chosen such that is antiperiodic in and , while it is periodic in . When taking the Fourier transform, the time differences and are associated with fermionic Matsubara frequencies , , while is associated with a bosonic frequency . Note that the relation to the representation with four frequencies is such that


The measurement formula for the correlator reads


in analogy to . The correlator can thus be measured at essentially no additional computational cost together with the two-particle Green’s function.

Figure 11: (Color online) Spin-up-down component of the real part of the reducible vertex function for a correlated impurity coupled to a single bath-level. The parameters are the same as in Fig. 7. Lines show the exact diagonalization result and open symbols show the improved Monte Carlo measurement. Black symbols correspond to the usual Matsubara axis measurement based on Eqs. (22)-(25) (note that the measurement in the Legendre basis leads to systematic deviations rather than the irregular noise observed in the Matsubara measurement).

Single-particle quantities can be measured directly in imaginary time with essentially arbitrarily fine resolution, since the number of imaginary time bins does not influence the performance of the algorithm. For two-particle quantities, which depend on three independent time differences, this is not the case, due to memory restrictions. Hence we measure these quantities directly in the Matsubara or the Legendre basis. Instead of measuring the correlator for all flavors , we accumulate in order to save memory. However, this still requires one to measure, in addition to the two-particle Green’s function, an object which is of the same size. In order to avoid this, one can replace on the right hand side of Eq. (26) by the sum of its connected and disconnected parts, which leads to


where is defined as in (20). Note that here the disconnected part instead of the full two-particle Green’s function appears on the right hand side. In the following we assess the quality of the results obtained from the different measurement formulas for the vertex function.

As a test, we again consider a correlated impurity hybridized to a single bath level. This allows us to compare Monte Carlo data to exact results. Figure 9 shows the real part of the spin-up-up component of the connected part of the two-particle Green’s function, , as a function of the fermionic Matsubara frequency for fixed (, ). Circles, squares and triangles correspond to the bosonic frequencies , , and , respectively. The improved Monte Carlo measurements (open and closed symbols), as well as the usual Matsubara-axis measurements (crosses) agree very well with the exact diagonalization result (lines). The problems caused by the stochastic noise are not evident at the level of the two-particle Green’s function, but only at the level of the vertex.

The real part of the spin-up-up component of the reducible vertex is shown in Fig. 10. Here, due to the multiplication with inverse Green’s functions, the noise in the data obtained from the standard Matsubara measurement grows considerably at large frequencies, while the improved measurement reproduces the correct high-frequency behavior. Both improved estimators, Eqs. (25), (26) and Eqs. (25), (31), yield results of comparable accuracy. In the case of larger bosonic frequencies, where the connected part becomes small, the improved accuracy of the new measurement procedure is evident even at the lowest .

We found that using the vertex function obtained via the improved estimators enhances the stability of dual fermion calculations. This is particularly the case for calculations in the insulating phase, where the improved estimators appear to be significantly more accurate and where it is otherwise difficult to obtain sufficient statistics due to the low perturbation order.

Figure 11 shows the real part of the spin-up-down component of the reducible vertex at low temperature and for the same parameters as in Fig. 7 (, ), measured in the Legendre basis. Again, one sees that at large bosonic , deviations between the exact result (lines) and the standard Matsubara measurement (black symbols) appear already at small , while the improved estimators (symbols) yield more accurate data.

Iii Application

iii.1 Two-orbital model

As an application of the methods described above, we calculate the self-energy and vertex function for a two-orbital model within the single-site DMFT. We consider the Hubbard model on the Bethe lattice with semi-elliptical density of states with bandwidth . First, we re-examine the spin-freezing transition reported in Ref. Werner et al., 2008 for a three-orbital model, and show that qualitatively the same physics is found in the two-orbital model, with the interaction restricted to density-density terms. We will work with the Hamiltonian (1), with the interaction part explicitly given by


where is the orbital index, denotes spin, and are the intra- and inter-orbital Coulomb interaction parameters, is the Hund’s rule coupling coefficient and .

iii.2 Spin-freezing transition

Figure 12: (Color online) Phasediagram in the space of density and interaction of the two-orbital Hubbard model on the Bethe lattice with and . The thick black lines mark regions of Mott insulating behavior. Circles, squares and stars, respectively, mark the parameters for which the calculations in Figs. 15, 16 and 17 have been performed.

The phasediagram in the plane of filling and interaction is shown in Fig. 12 for and temperature . It reproduces the qualitative features reported in Ref. Werner et al., 2008 for the three-orbital model: For small values of the density and interaction, we find a Fermi liquid phase, while for larger values, the model exhibits a frozen moment phase. The frozen moment phase is characterized by a spin-spin correlation function that approaches a non-zero constant at large times, as shown in Fig. 13. In this phase, the spin-spin correlation function at (which we refer to as ), is expected to be independent of temperature. In a Fermi liquid at low temperature, on the other hand, the spin-spin correlation function behaves as for times sufficiently far from or , respectively, so that . In Fig. 14 we show the ratio , as a function of filling, which clearly confirms this behavior. The ratio changes from the value expected in the Fermi liquid phase to expected for frozen moments, passing through near the spin-freezing transition, indicating a -linear behavior. Hence, in the vicinity of the transition, the spin-spin correlation function decays unusually slowly, , consistent with the findings in Ref. Werner et al., 2008.

Figure 13: (Color online) Imaginary-time dependence of the spin-spin correlation function for and densities and (from top to bottom). For small densities, in the Fermi liquid regime, the spin-spin correlation at is proportional to and hence approaches zero as , while for high densities the spin correlation function approaches a constant, indicating the presence of frozen moments.

We note that a similar phasediagram for a two-orbital model, which plots the “quasi-particle weight” in the space of filling and interaction strength, has recently been reported in Ref. de’ Medici et al., 2011. Our spin-freezing transition line seems to resemble the contour-lines for fixed in Ref. de’ Medici et al., 2011, although at the temperature of our calculation, the strong deviations from Fermi-liquid behavior near the transition mean that in this regime a quasi-particle weight cannot be properly defined. It is an interesting open question whether and how the properties of a (strongly renormalized) Fermi-liquid are recovered as .

Figure 14: (Color online) Ratio of the spin-spin correlation function for temperatures and as a function of density. Vertical lines mark the approximate location of the transition, where the self-energy is consistent with a behavior with . On the Fermi liquid side, obviously , while in the frozen moment phase.

The development of frozen moments is accompanied by a simultaneous change in the low-frequency behavior of the self-energy. In Fig. 15 we plot the imaginary part of the self-energy as a function of Matsubara frequencies for and the densities indicated by circles in the phasediagram in Fig. 12. For small densities, the imaginary part of the self-energy extrapolates to zero as expected for a Fermi liquid. In the presence of frozen moments, however, the electrons are expected to be scattered by these moments, resulting in a non-zero value . This behavior is evident from the imaginary time data. We find that in the transition region the low-frequency behavior of the self-energy is well described by a power-law . The phase boundary in Fig. 12 corresponds to a “square-root” frequency dependence, i. e. . Note that the horizontal axis is , so that corresponds to a linear increase at low frequencies (black line in Fig. 15).

Figure 15: (Color online) Imaginary part of the self-energy for various fillings for the two-orbital model at , obtained using the improved estimator measured in the Legendre basis. The positions of the corresponding parameters are marked by circles in the phasediagram in Fig. 12. The solid line is a fit proportional to , with .
Figure 16: (Color online) Self-energy on the real axis obtained from the Matsubara axis data by analytical continuation using Padé approximants. Thick lines (red, gray in print) show the result obtained from the improved estimator measured in the Legendre basis. The dashed and dotted lines show results obtained from the improved estimator and Dyson’s equation, respectively, measured in Matsubara basis. The positions of the parameters in the phasediagram are marked by squares in Fig. 12. In the Fermi liquid phase (left panel) the data obtained from the improved estimators measured in the Legendre basis is well described by the Fermi liquid form over a wide energy range (thin solid lines), while the data obtained from the two other methods exhibits spurious features.

It is instructive to compute the self-energy on the real axis. Analytical continuation is a delicate issue, in particular for Monte Carlo data, due to the presence of statistical noise. Given the high quality of the Matsubara-axis data of Fig. 15, which have been obtained using the combination of improved estimator and Legendre noise filter, we nevertheless attempt to perform an analytical continuation using Padé approximants.Vidberg and Serene (1977) The results are shown by thick red (gray) lines in Fig. 16. The three values of the density correspond to the positions marked by squares in the phasediagram in Fig. 12. We first observe that for , the real axis self-energy has the expected Fermi liquid form and is well described by the expression over a wide energy range (thin solid line). For comparison we have also plotted the results obtained from the self-energy measured in Matsubara representation, using Dyson’s equation (dotted lines) and the improved estimator (dashed). In particular the self-energy obtained from Dyson’s equation shows spurious features and a departure from the expected Fermi liquid behavior already for energies close to the Fermi level.

By construction, the Padé analytical continuation yields the most accurate result for small . As a consistency check, we performed a least-squares fit of the Padé data to the Fermi liquid form in the energy window (thin solid line), finding and , which corresponds to a weakly renormalized Fermi liquid with . The value of agrees remarkably well with the value obtained from a second-order polynomial extrapolation of the self-energy, for the improved estimator measured in Legendre ( for the self-energy obtained from Dyson’s equation). The deviation is less than 0.5%, while the fit of the Padé analytical continuation of the self-energy obtained from Dyson’s equation already disagrees by more than 7%. We note that the Monte Carlo error on the self-energy is of the order of . A summary of the parameters extracted from the Padé analytical continuation for the various measurement procedures is given in Table 1.

Leg. + impr. est. impr. est. Dyson
0.02 0
0.1 12
0.25 32
0.5 64
Table 1: Coefficient obtained from a least-squares fit of the Padé data in the energy window using #df degrees of freedom. The values should be compared with the result obtained from a second-order polynomial extrapolation of the imaginary part of the Matsubara self-energy, which yields .
Figure 17: (Color online) Real and imaginary parts of the reducible vertex function for the two-orbital Hubbard model for and across the spin-freezing transition. The left panels shows intra- and the right panels interorbital components of the spin (closed symbols) and charge (open symbols) channels. Circles (red), triangles (green) and squares (blue) correspond to densities and , respectively. The position of these parameters in the phasediagram is marked by stars in Fig. 12 (in the corresponding color). The y-axis range in the upper left panel has been restricted to improve visualization.

The self-energy for (middle panel) shows incipient non-Fermi liquid behavior. Both the real and imaginary parts develop a peak around the Fermi level, qualitatively consistent with the appearance of a “square-root” non-analyticity near the spin-freezing transition. While the results from the improved estimators agree fairly well, the result obtained using Dyson’s equation deviates rather strongly. Finally, for the self-energy clearly exhibits non-Fermi liquid behavior, with even more pronounced peaks around the Fermi level. The non-zero value of the imaginary part at zero frequency is again consistent with the value obtained from a polynomial extrapolation of the Matsubara self-energy. As a final illustration for this model, we calculate the reducible vertex function according to Eqs. (25)-(26). Figure 17 illustrates the evolution of the vertex for fixed and bosonic Matsubara frequency across the spin freezing transition, both for intra-orbital (left panels) and inter-orbital (right panels) components. The vertex is essentially featureless on the Fermi liquid side (, red circles), but develops structure as the spin-freezing transition is approached (, green triangles), and notably in the the frozen moment phase (, blue squares). In particular the spin-channel (closed symbols) is strongly enhanced in the frozen moment phase.

One can anticipate that these features will induce significant changes when the vertex is used to calculate the momentum dependent self-energy in diagrammatic extensions of dynamical mean-field theory.

iii.3 High-spin to low-spin transition

Figure 18: (Color online) Phasediagram in the space of crystal field splitting and interaction strength of the half-filled two-orbital Hubbard model on the Bethe lattice with bandwidth , and . The parameters for which the vertex function in Fig. 20 has been calculated are marked by stars.
Figure 19: (Color online) Total moment , moment of orbital 1, and density of orbital 1, as a funciton of for various values of the crystal field splitting . From left to right (and top to bottom) the curves correspond to the values .
Figure 20: (Color online) Real and imaginary parts of the reducible vertex function for the two-orbital Hubbard model for and at for the points indicated by stars (in corresponding color) in the phase diagram of the low-spin to high-spin transition in Fig. 18. The left panel shows intraorbital (for orbital 2) and the right panel interorbital components of the spin (closed symbols) and charge (open symbols) channels. Circles (red), triangles (green) and squares (blue) correspond to and . The corresponding points in the phasediagram are marked by stars in Fig. 18 (in corresponding color).

As a second application, we study the same model at half-filling and additionally consider a crystal-field splitting term , defined as


in Eq. (1). This model has recently been considered as a minimal model for the physics of LaCoO.Kuneš and Křápek (2011) The condition for half-filling is found by taking the chemical potential to be half the sum over the elements of any row or column of the interaction matrix, which gives . For this model, a high-spin to low-spin, and associated orbital polarization transition have been reported.Werner and Millis (2007) To quantify the effect of our density-density approximation, we compute the phasediagram in the space of crystal field splitting and interaction strength. Comparison of the result shown in Fig. 18 with the phasediagram for the rotationally invariant interaction (including spin-flip and pair-hopping terms) reported in Ref. Werner and Millis, 2007 shows that the differences are rather marginal and that all qualitative features are reproduced by the density-density approximation. In particular we find a high-spin Mott insulating phase at large interaction and small crystal field splitting, and a low-spin orbitally polarized insulator at small interaction and large crystal field splitting. These two qualitatively distinct insulating phases are separated by a metallic phase with an end-point near .

The nature of the two insulating phases becomes evident if one plots the spin expectation value (total moment or moment of orbital 1 ) and the occupancy of orbital 1 () as a function of interaction strength. Several such traces corresponding to fixed values are shown in Fig. 19. In the orbitally polarized insulator, the occupancy of orbital 1 (which is shifted up in energy) is very low. The filling continuously increases with through the metallic phase and eventually jumps to at the transition into the high-spin Mott insulating phase. A similar behavior is seen in the traces for .

We finally show in Fig. 20 the evolution of the reducible vertex function in the metallic phase, as one moves from the phase boundary to the orbitally polarized insulator to the phase boundary with the high-spin Mott insulator. The three data points, corresponding to and , and are marked in the phasediagram (Fig. 18) by star symbols. While rather featureless near the phase boundary to the orbitally polarized insulator, the vertex develops structure as the moment increases (cf. the arrows in the upper panel of Fig. 19).

Iv Conclusions and outlook

We presented efficient measurement procedures for the self-energy and vertex function within the hybridization expansion CTQMC approach. The improved estimators are based on higher-order correlation functions, which can be particularly easily measured in single-site (multi-orbital) models with density-density interactions, but also for generic models. In combination with a recently proposed noise-filtering scheme (Legendre representation), we were able to completely eliminate the noise problems at intermediate to high frequencies, which have plagued the “standard” evaluation of the self-energy and vertex function using the hybridization expansion method. With the noise-problem solved, one can fully exploit the performance advantages of CT-HYB in the strong-correlation regime and in the case of (single-site) multiorbital models. The accurate measurement of two-particle Green’s functions and vertex functions should enable dual fermionRubtsov et al. (2008); Hafermann et al. (2009) or dynamical vertex approximationToschi et al. (2007) calculations of multi-band systems, and thus (in combination with band-structure input) the ab-initio simulation of transition metal compounds which capture the effect of non-local correlations. While we have presented results for two-orbital models, the computational effort for the measurement of the vertex function scales as the square of the number of orbitals (for diagonal hybridization) and thus simulations are feasible even for five-orbital models on a small computer cluster.

We have further demonstrated the efficiency of the improved measurements with self-energy and vertex data for a two-orbital model with density-density interactions. In particular, we have revisited the phenomenon of the “spin freezing” transition, which was originally discovered in a three-orbital model, but which manifests itself in a very analogous manner also in the two-orbital case. In combination with the results by Ishida and Liebsch Ishida and Liebsch (2010) for five-orbital models, this establishes the spin-freezing phenomenon as a generic feature of multi-orbital models with large Hund’s rule coupling. This phenomenon leads to strong non-Fermi liquid effects in certain regions of the paramagnetic metallic phase.

We would like to thank L. Boehnke, M. Eckstein, M. Ferrero, J. Mravlje, and O. Parcollet for valuable discussions. The simulations have been performed on the Brutus cluster at ETH Zürich, using an implementation based on the ALPS-libraries.Bauer et al. (2011)

Appendix A Self-energy

Following Ref. Bulla et al., 1998, we derive the expression for the self-energy for the Hamiltonian (1). The single-particle Green’s function is defined as


In order to arrive at the equation of motion, for reasons outlined in appendix B, we take the derivative with respect to its second argument,


The equation of motion for the creation operator is


and the commutator with the Hamiltonian (1) evaluates to


Note that the commutator in (37) involving the interaction term generates two terms with different order of the operators and . These can be commuted since the interaction matrix for identical flavors is zero: . Inserting this expression into (35), we obtain


where we have used the above definition of Green’s function and introduced the functions


Using the definition of the Fourier transform


(38) can be written in Fourier space as


Note that here no summation over repeated indices is implied (unless explicitly indicated). In DMFT it is customary to define the bare Green’s function of the impurity model as


where is the hybridization function. Subtracting on both sides of (41), we can rewrite its left hand side as