Symmetric Improved Estimators for Continuous-time Quantum Monte Carlo

TL;DR

This paper introduces symmetric improved estimators for continuous-time quantum Monte Carlo simulations of the multi-orbital Anderson impurity model, significantly reducing noise at high frequencies and enhancing computational accuracy.

Contribution

It derives new equations of motion for Green's functions using symmetric differentiation, enabling improved estimators that mitigate high-frequency noise in Monte Carlo simulations.

Findings

Reduced noise at large Matsubara frequencies

Enhanced accuracy of Green's function measurements

Applicability to multi-orbital Anderson impurity models

Abstract

We derive equations of motion for Green's functions of the multi-orbital Anderson impurity model by differentiating symmetrically with respect to all time arguments. The resulting equations relate the one- and two-particle Green's function to correlators of up to six particles at four times. As an application we consider continuous-time quantum Monte Carlo simulations in the hybridization expansion, which hitherto suffered from notoriously high noise levels at large Matsubara frequencies. Employing the derived symmetric improved estimators overcomes this problem.

Tables

Table 1. Table 1: Terms of the symmetric improved estimators. We use the abbreviations ( U ​ n ) a ​ b = ∑ j ​ k U [ a ​ j ] ​ [ b ​ k ] ​ d j † ​ d k subscript 𝑈 𝑛 𝑎 𝑏 subscript 𝑗 𝑘 subscript 𝑈 delimited-[] 𝑎 𝑗 delimited-[] 𝑏 𝑘 subscript superscript 𝑑 † 𝑗 subscript 𝑑 𝑘 (Un)_{ab}=\sum_{jk}U_{[aj][bk]}d^{\dagger}_{j}d_{k} , ( U ​ d † ​ d † ) a ​ b = ∑ j ​ k U j ​ k ​ [ a ​ b ] ​ d j † ​ d k † subscript 𝑈 superscript 𝑑 † superscript 𝑑 † 𝑎 𝑏 subscript 𝑗 𝑘 subscript 𝑈 𝑗 𝑘 delimited-[] 𝑎 𝑏 subscript superscript 𝑑 † 𝑗 subscript superscript 𝑑 † 𝑘 (Ud^{\dagger}d^{\dagger})_{ab}=\sum_{jk}U_{jk[ab]}d^{\dagger}_{j}d^{\dagger}_{k} and ( U ​ d ​ d ) a ​ b = ∑ j ​ k U [ a ​ b ] ​ j ​ k ​ d k ​ d j subscript 𝑈 𝑑 𝑑 𝑎 𝑏 subscript 𝑗 𝑘 subscript 𝑈 delimited-[] 𝑎 𝑏 𝑗 𝑘 subscript 𝑑 𝑘 subscript 𝑑 𝑗 (Udd)_{ab}=\sum_{jk}U_{[ab]jk}d_{k}d_{j} . Frequencies ν ( i ) subscript 𝜈 𝑖 \nu_{(i)} and ω 𝜔 \omega are fermionic and bosonic Matsubara frequencies, respectively.

Name	Formula
${\xi }_{ab}^{\nu }$	$-{\int }_0^{\beta }𝑑{\tau }_1𝑑{\tau }_2⟨T_{\tau }{q}_a\left({\tau }_1\right)d_b^{†}\left({\tau }_2\right)⟩e^{i\nu \left({\tau }_1-{\tau }_2\right)}/\beta$
${\vartheta }_{ab}^{\nu }$	$-{\int }_0^{\beta }𝑑{\tau }_1𝑑{\tau }_2⟨T_{\tau }{q}_a\left({\tau }_1\right)q_b^{†}\left({\tau }_2\right)⟩e^{i\nu \left({\tau }_1-{\tau }_2\right)}/\beta$
${\varphi }_{abcd}^{\omega }$	$\int 𝑑{\tau }_1𝑑{\tau }_2⟨T_{\tau }{\left(Un\right)}_{ab}\left({\tau }_1\right){\left(Un\right)}_{cd}\left({\tau }_2\right)⟩e^{i\omega \left({\tau }_1-{\tau }_2\right)}/\beta$
${\psi }_{abcd}^{\omega }$	$\int 𝑑{\tau }_1𝑑{\tau }_2⟨T_{\tau }{\left(Udd\right)}_{ac}\left({\tau }_1\right){\left(U{d}^{†}{d}^{†}\right)}_{bd}\left({\tau }_2\right)⟩e^{i\omega \left({\tau }_1-{\tau }_2\right)}/\beta$
$f_{abcd}^{\nu \omega }$	$\int 𝑑{\tau }_1𝑑{\tau }_2𝑑{\tau }_3⟨T_{\tau }{q}_a\left({\tau }_1\right)q_b^{†}\left({\tau }_2\right){\left(Un\right)}_{cd}\left({\tau }_3\right)⟩e^{i\nu \left({\tau }_1-{\tau }_2\right)+i\omega \left({\tau }_2-{\tau }_3\right)}/\beta$
$g_{abcd}^{\nu {\nu }^{\prime }}$	$\int 𝑑{\tau }_1𝑑{\tau }_2𝑑{\tau }_3⟨T_{\tau }{q}_a\left({\tau }_1\right)q_c\left({\tau }_2\right){\left(U{d}^{†}{d}^{†}\right)}_{bd}\left({\tau }_3\right)⟩e^{i\nu \left({\tau }_1-{\tau }_3\right)+i{\nu }^{\prime }\left({\tau }_2-{\tau }_3\right)}/\beta$
$h_{abcd}^{{\nu }_1{\nu }_2{\nu }_3{\nu }_4}$	$\int 𝑑{\tau }_1𝑑{\tau }_2𝑑{\tau }_3𝑑{\tau }_4⟨T_{\tau }{q}_a\left({\tau }_1\right)q_b^{†}\left({\tau }_2\right)q_c\left({\tau }_3\right)q_d^{†}\left({\tau }_4\right)⟩e^{i\left({\nu }_1{\tau }_1-{\nu }_2{\tau }_2+{\nu }_3{\tau }_3-{\nu }_4{\tau }_4\right)}/\beta$

Equations

$$H_{\mathrm{AIM}}$$

$${\Delta }_a^{\nu }=\sum _K\frac{V_K^{aa}(V_K^{aa}{)}^{\ast }}{i\nu -{\epsilon }_{Ka}},$$

$${\mathcal{G}}_a^{\nu }=\frac{1}{i\nu -{\tilde{\epsilon }}_a-{\Delta }_a^{\nu }}.$$

$$G_a({\tau }_1,{\tau }_2)=-⟨T_{\tau }{d}_a^{}({\tau }_1)d_a^{†}({\tau }_2)⟩,$$

$$G_{abcd}({\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4)=⟨T_{\tau }{d}_a^{}({\tau }_1)d_b^{†}({\tau }_2)d_c^{}({\tau }_3)d_d^{†}({\tau }_4)⟩.$$

$$G_a^{\nu }=\frac{1}{\beta }{\int }_0^{\beta }d{\tau }_{1}d{\tau }_2{e}^{i\nu ({\tau }_1-{\tau }_2)}{G}_a({\tau }_1,{\tau }_2)$$

$$\begin{gathered} G_{abcd}^{{\nu }_1{\nu }_2{\nu }_3{\nu }_4}=\frac{1}{{\beta }^2}{\int }_0^{\beta }d{\tau }_{1}d{\tau }_{2}d{\tau }_{3}d{\tau }_4 \\ {e}^{i({\nu }_1{\tau }_1-{\nu }_2{\tau }_2+{\nu }_3{\tau }_3-{\nu }_4{\tau }_4)}{G}_{abcd}({\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4) \end{gathered}$$

$${\nu }_1+{\nu }_3={\nu }_2+{\nu }_4.$$

$$G_{a}^{\nu}=\mathcal{G}_{a}^{\nu}\big{(}1+\xi_{a}^{\nu}\big{)},$$

$${\xi }_{ab}^{\nu }={\xi }_a^{\nu }{\delta }_{ab}=\frac{1}{\beta }{\int }_0^{\beta }d{\tau }_{1}d{\tau }_2\underset{\equiv {\xi }_a({\tau }_1,{\tau }_2)}{\underbrace{⟨-T_{\tau }{q}_a^{}({\tau }_1)d_b^{†}({\tau }_2)⟩}}{e}^{i\nu ({\tau }_1-{\tau }_2)}$$

$$q_a$$

$$q_a^{†}$$

$$G_{a}^{\nu}=\mathcal{G}_{a}^{\nu}\bigg{(}1+\mathcal{G}_{a}^{\nu}\Big{(}2\sum_{j}U_{[aj][aj]}n_{j}+\vartheta_{a}^{\nu}\Big{)}\bigg{)},$$

$$\frac{1}{2}\left(U_{ijkl}-U_{jikl}\right)=:U_{[ij]kl}$$

$$\frac{1}{2}\left(U_{ijkl}-U_{ijlk}\right)=:U_{ij[kl]}.$$

$${\vartheta }_{ab}^{\nu }={\vartheta }_a^{\nu }{\delta }_{ab}=-\frac{1}{\beta }{\int }_0^{\beta }d{\tau }_{1}d{\tau }_2⟨T_{\tau }{q}_a^{}({\tau }_1)q_b^{†}({\tau }_2)⟩e^{i\nu ({\tau }_1-{\tau }_2)}.$$

$${\Sigma }_a^{\nu }=\frac{2{\sum }_j{U}_{[aj][aj]}{n}_j+{\vartheta }_a^{\nu }}{1+{\mathcal{G}}_a^{\nu }(2{\sum }_j{U}_{[aj][aj]}{n}_j+{\vartheta }_a^{\nu })}.$$

$$\mathrm{var}[\Sigma ]={\left|\frac{\partial \Sigma }{\partial \vartheta }\right|}^2\mathrm{var}[\vartheta ].$$

$$G^{\nu_{1}\nu_{2}\nu_{3}\nu_{4}}_{abcd}=\mathcal{G}^{\nu_{1}}_{a}\Bigg{(}-R_{1,abcd}^{\nu_{1}\nu_{2}\nu_{3}\nu_{4}}+\mathcal{G}^{\nu_{2}}_{b}\bigg{(}R_{2,abcd}^{\nu_{1}\nu_{2}\nu_{3}\nu_{4}}+\mathcal{G}^{\nu_{3}}_{c}\Big{(}-R_{3,abcd}^{\nu_{1}\nu_{2}\nu_{3}\nu_{4}}+\mathcal{G}^{\nu_{4}}_{d}\big{(}R_{4,abcd}^{\nu_{1}\nu_{2}\nu_{3}\nu_{4}}+h_{abcd}^{\nu_{1}\nu_{2}\nu_{3}\nu_{4}}/\beta\big{)}\Big{)}\bigg{)}\Bigg{)}\mbox{ where }$$

$$R_1^{{\nu }_1{\nu }_2{\nu }_3{\nu }_4}$$

$$R_2^{{\nu }_1{\nu }_2{\nu }_3{\nu }_4}$$

$$R_3^{{\nu }_1{\nu }_2{\nu }_3{\nu }_4}$$

$$R_4^{{\nu }_1{\nu }_2{\nu }_3{\nu }_4}$$

$$G_{abcd}^{{\nu }_1{\nu }_2{\nu }_3{\nu }_4}$$

$${\mathfrak{F}}_{abcd}^{{\nu }_1{\nu }_2{\nu }_3{\nu }_4}$$

$$\displaystyle+2U_{[ac][bd]}+2\sum_{i}\Big{(}U_{[ic][bd]}\xi_{ai}^{\nu_{1}}+U_{[ac][id]}\xi_{bi}^{\nu_{2}}+U_{[ai][bd]}\xi_{ci}^{\nu_{3}}+U_{[ac][bi]}\xi_{di}^{\nu_{4}}\Big{)}$$

$$\underset{(i)}{\underbrace{-4{\varphi }_{abcd}^{{\nu }_1-{\nu }_2}+2f_{abcd}^{{\nu }_1({\nu }_1-{\nu }_2)}+2f_{cdab}^{{\nu }_3({\nu }_2-{\nu }_1)}}}\underset{(ii)}{\underbrace{+4{\varphi }_{adcb}^{{\nu }_1-{\nu }_4}-2f_{adcb}^{{\nu }_1({\nu }_1-{\nu }_4)}-2f_{cbad}^{{\nu }_3({\nu }_4-{\nu }_1)}}}\underset{(iii)}{\underbrace{-{\psi }_{abcd}^{{\nu }_1+{\nu }_3}-g_{abcd}^{{\nu }_1{\nu }_3}-g_{dcba}^{{\nu }_4{\nu }_2}}}+h_{abcd}^{{\nu }_1{\nu }_2{\nu }_3{\nu }_4}.$$

$$G_{abcd}^{{\nu }_1{\nu }_2{\nu }_3{\nu }_4}=({\delta }_{12}-{\delta }_{14})G_a^{{\nu }_1}{G}_c^{{\nu }_3}-\frac{1}{\beta }{G}_a^{{\nu }_1}{G}_b^{{\nu }_2}{G}_c^{{\nu }_3}{G}_d^{{\nu }_4}{F}_{abcd}^{{\nu }_1{\nu }_2{\nu }_3{\nu }_4}.$$

$$F_{abcd}^{{\nu }_1{\nu }_2{\nu }_3{\nu }_4}=\frac{{\mathcal{G}}_a^{{\nu }_1}{\mathcal{G}}_b^{{\nu }_2}{\mathcal{G}}_c^{{\nu }_3}{\mathcal{G}}_d^{{\nu }_4}}{G_a^{{\nu }_1}{G}_b^{{\nu }_2}{G}_c^{{\nu }_3}{G}_d^{{\nu }_4}}{\mathfrak{F}}_{abcd}^{{\nu }_1{\nu }_2{\nu }_3{\nu }_4}.$$

$$G_{abcd}^{{\nu }_1{\nu }_2}=\frac{1}{\beta }$$

$$\times ⟨T_{\tau }{d}_a^{}({\tau }_1)d_b^{†}({\tau }_2)[d_c^{}{d}_d^{†}]({\tau }_3)⟩,$$

$$G_{abcd}^{{\nu }_1{\nu }_2}=$$

$$-{\mathcal{G}}_a^{{\nu }_1}{\mathcal{G}}_b^{{\nu }_2}[\beta {\delta }_{12}{n}_c(2(Un{)}_a+{\vartheta }_a^{{\nu }_1})$$

$$-{\delta }_{14}{\xi }_a^{{\nu }_1}{\xi }_b^{{\nu }_2}-2{\hat{\varphi }}_{abcd}^{{\nu }_1-{\nu }_2}+{\hat{f}}_{abcd}^{{\nu }_1({\nu }_1-{\nu }_2)}],$$

$${\hat{\varphi }}_{abcd}^{\omega }=\frac{1}{\beta }$$

$$⟨T_{\tau }\sum _{jk}{U}_{[aj][bk]}[d_j^{†}{d}_k]({\tau }_1)[d_d^{†}{d}_c]({\tau }_2)⟩$$

Full text

I Introduction

The Hubbard model (HM)[1] and the Anderson impurity model (AIM)[2], which are related through the dynamical mean-field theory (DMFT)[3], are two of the basic models for strongly correlated electrons. Thus, numerical and analytic solutions of these models over a wide range of parameters are of great interest in condensed matter physics. As of today, the continuous-time quantum Monte Carlo (CT-QMC)[4] method is the workhorse for obtaining numerical solutions in terms of one- and two-particle Green’s functions.

CT-QMC algorithms are based on a stochastic sampling of the diagrammatic series expansion of either the partition function or, directly, the thermal expectation value of some operators, which is also known as worm sampling.[5, 6] For the AIM, one distinguishes between weak coupling expansions in the interaction CT-INT[7] (and related CT-AUX[8]) and strong coupling expansions in the hybridization CT-HYB.[9] While traditionally the former are employed for single-orbital impurity model calculations and impurity clusters, the latter are primarily used for multi-orbital impurity models. The reason for this is that in weak coupling, the exponential scaling of the sign problem makes multi-orbital calculations with non-density-density interactions difficult, while in strong coupling, the exponential scaling of the local Hamiltonian dimensions forbids medium- to large scale clusters.

Conceptionally one might expect CT-INT and CT-HYB to behave similarly apart from their differences in scaling with interaction and hybridization strength, respectively. The estimators for one- and two-particle Green’s functions are instead considerably different in weak- and strong coupling approaches. While Green’s function estimators in CT-INT are formulated as corrections to the non-interacting Green’s function[10], this is not the case for CT-HYB. This results in poor asymptotic behavior of derived quantities, which in turn triggered a series of developments in the CT-HYB community attempting to remedy this problem. These developments include orthogonal polynomial representations as effective low-pass filters[11], moment expansions in the one-particle self-energy[12, 13], asymptotic expressions for the two-particle vertex functions[14, 15] and approaches based on the equation of motion (EOM), often referred to as improved estimators[16, 17, 18]. Similar to CT-INT, the improved estimators of CT-HYB calculate the Green’s function as a correction to the non-interacting Green’s function.

In this paper we introduce and explore the concept of symmetric improved estimators. That is, we derive EOMs by differentiating with respect to all time arguments. The derived equations relate the one- and two-particle Green’s function to correlation functions of higher order in the number of creation and annihilation operators. These EOMs may prove useful in various contexts. We employ them for CT-HYB estimators of the self-energy, two-particle (four-leg) and three-leg vertex. We give theoretical arguments showing that the symmetric improved estimators have a strongly reduced noise level at large Matsubara frequencies $\nu$, and even a different scaling with $\nu$, which is confirmed in actual CT-HYB calculations.

Section II introduces the AIM Hamiltonian and our notation for the one- and two-particle Green’s function. Section III sketches the derivations for the symmetric improved estimators on the one- and two-particle level. We note that the expressions are based on hierarchies of EOMs and are also useful outside the context of CT-QMC.[19, 20] We further discuss the improved Monte Carlo error scaling of Green’s functions, self-energies and vertex functions at large Matsubara frequencies. Section IV discusses the implementation of the symmetric improved estimators in CT-HYB. An implementation for density-density interactions in segment CT-HYB[21] is straight-forward, an implementation for general interactions requires worm sampling. We discuss drawbacks for a classical worm sampling implementation and propose methods to increase the sampling efficiency. In Section V we show results for self-energies as well as three- and four-leg vertex functions and validate them by comparing to results from exact diagonalization (ED). Finally, Section VI summarizes our work. We furthermore provide an appendix that contains the derivation of the presented formulas in more detail.

II Definitions and model

The AIM consists of an impurity site that is coupled to a bath and its Hamiltonian reads

[TABLE]

Here, $d_{i}$ ($d_{i}^{\dagger}$) is the annihilation (creation) operator of an electron with spin-orbital flavor $i$ on the impurity; $c_{Ki}$ ($c_{Ki}^{\dagger}$) is the annihilation (creation) operator of an electron with impurity flavor $i$ in the non-interacting bath and $K$ subsumes the remaining bath degrees of freedom (e. g. the momentum $\mathbf{k}$). The impurity is described by a local one-particle potential $\tilde{\varepsilon}_{i}$ (e.g. the crystal field), the interaction matrix $U_{ijkl}$, the bath dispersion $\varepsilon_{Ki}$, and the hybridization strength $V_{K}^{ij}$. Since the bath degrees of freedom appear in Eq. (II) at most quadratically, they can be formally integrated out yielding a one-body term in the impurity operators. This contains the hybridization function $\Delta_{ab}(\tau_{1},\tau_{2})$ encoding the entire influence of the bath. In the following we assume a diagonal hybridization $\Delta_{ab}(\tau_{1},\tau_{2})\equiv\Delta_{a}(\tau_{1},\tau_{2})\delta_{ab}$ resulting in a diagonal one-particle Green’s function $G^{\tau_{1}\tau_{2}}_{ab}=G^{\tau_{1}\tau_{2}}_{a}\delta_{ab}$. But the same concept of symmetric improved estimator equations can be extended straight-forwardly to non-diagonal hybridizations.

The hybridization function is most conveniently written in Matsubara frequencies as

[TABLE]

and it relates to the non-interacting Green’s function of the impurity as

[TABLE]

We define the interacting one-particle Green’s function of the AIM as

[TABLE]

where $d_{a}(\tau)$ ($d_{a}^{\dagger}(\tau)$) are now the annihilation (creation) operators for electrons of flavor $a$ at (imaginary) time $\tau$. Furthermore, $T_{\tau}$ is the imaginary-time ordering operator, and $\langle\ldots\rangle=(\text{Tr}e^{-\beta H}\ldots)/Z$ the thermal expectation value at temperature $T$ ($\beta=1/T$), $Z$ is the partition function. The Green’s function is related to the density as $n_{a}\equiv\langle d^{\dagger}_{a}d^{\phantom{{\dagger}}}_{a}\rangle=1+G_{a}(0+,0)$.

Analogously, we define the two-particle Green’s function of the AIM as

[TABLE]

The Fourier transforms to fermionic Matsubara frequencies[22] $\nu=(2l+1)\pi/\beta$ at integer numbers $l$ are given by

[TABLE]

and

[TABLE]

for the one- and two-particle Green’s function, respectively.

The time translation symmetry of $H_{\mathrm{AIM}}$ results in time translation symmetry of Green’s functions in imaginary time. This is equivalent to energy conservation, which for the two-particle Green’s function reads

[TABLE]

For the two-particle Green’s function it is more common to assume a mixed bosonic-fermionic frequency representation with two fermionic- and one bosonic frequency. However, the choice of these frequencies is ambiguous and therefore the derivations of this paper are done in four fermionic frequencies. With the definitions above, the reader can easily adapt the results to their favorite convention. (Sometimes we use generic bosonic frequencies $\omega$ which are however replaced by $\nu_{1}\ldots\nu_{4}$ in the final expression.)

III Symmetric improved estimators

III.1 One-particle estimator

It is well-known that the Heisenberg equation of motion for a one-particle Green’s function is an equation that connects the one- and two-particle Green’s functions. This has been exploited in the CT-HYB[16, 17, 18] and numerical renormalization group (NRG)[23] algorithms. It leads to the so-called “improved estimator” equation

[TABLE]

where

[TABLE]

is a two-particle Green’s function and $q$ contracts three operators at equal time,111 One may also refer to $\xi$ as $G\Sigma$. i.e.,

[TABLE]

The explicit derivation of this equation of motion can be found in Appendix A and is based on the derivative with respect to the first time argument of $G_{a}$.

In this paper, we now express $\xi_{a}(\tau_{1},\tau_{2})$ by using the equation of motion again. This time, we apply it to the second time argument of $\xi_{a}(\tau_{1},\tau_{2})$, insert it into Eq. (9), and finally arrive at (for the detailed calculation see Appendix A)

[TABLE]

which is also shown by Feynman diagrams in Fig. 1. Here, we make use of the anti-symmetric $U$-matrix Eq. (39) in Appendix A

[TABLE]

and the following three-particle Green’s function with only two distinct time arguments is employed

[TABLE]

Error propagation.

The declared goal of improved estimators is to reduce high-frequency noise in the quantity of interest. Among one-particle quantities, noise is most prominent in the self-energy. Using Eq. (13) and the Dyson equation, we get

[TABLE]

The variance of the self-energy can hence approximately be computed by error propagation as (neglecting the error in the density $n_{i}$):

[TABLE]

Since the derivative goes to 1 in the limit of high frequencies, we conclude that the high-frequency noise amplitude is indeed identical to the noise amplitude of the estimator $\vartheta_{a}^{\nu}$, which should be independent of $\nu$ for large $\nu$. This result for the symmetric improved estimator can be compared to the behavior of the conventional improved estimator for which a discussion of the error propagation can be found elsewhere.[25] The main result is that the error of the conventional improved estimator grows linearly with $\nu$.

III.2 Two-particle estimator

For the symmetric improved estimator of the two-particle Green’s function, we obtain (again see the Appendix B for the derivation):

[TABLE]

[TABLE]

All estimators occurring on the right-hand side of Eqs. (20)-(23) are defined in Table 1. In order to obtain a more symmetric form, we further made use of the relations $\delta_{\nu_{1}\nu_{2}}\equiv\delta_{\nu_{1}\nu_{2}}\delta_{\nu_{3}\nu_{4}}$ and $\delta_{\nu_{1}\nu_{4}}\equiv\delta_{\nu_{2}\nu_{3}}\delta_{\nu_{3}\nu_{4}}$, which are valid due to energy conservation. Additionally we employed $\delta_{12}\equiv\delta_{ab}\delta_{cd}\delta_{\nu_{1}\nu_{2}}\delta_{\nu_{3}\nu_{4}}$ and $\delta_{14}\equiv\delta_{ad}\delta_{bc}\delta_{\nu_{1}\nu_{4}}\delta_{\nu_{2}\nu_{3}}$ to make the expressions shorter.

Inserting the $R_{i}$ terms from Eqs. (20)-(23) into Eq. (19) and regrouping the terms leads to the following expression for the symmetric improved two-particle Green’s function:

[TABLE]

(A pictorial representation of this formula is found in Fig. 2.) The quantities entering here are the antisymmetrized Coulomb interaction $U$, the one-particle density $n_{i}$ and the correlators $\vartheta$ and $\xi$ which already appeared for the one-particle Green’s function in Section III.1. In order to shorten the expression, we singled out the disconnected part $(\delta_{12}-\delta_{14})G_{a}^{\nu_{1}}G_{c}^{\nu_{3}}$ of the two-particle Green’s function by applying Eq. (13). However, there remains another term proportional to $(\delta_{12}-\delta_{14})$ in $\mathfrak{F}$. This is rooted in the expansion via equations of motion, which always leads to expressions involving non-interacting Green’s functions $\cal G$. The 10 terms in the last line of Eq. (III.2) are genuinely related to the two-particle Green’s function. Further, the frequency structure of the terms $(i)$-$(iii)$ resembles the contributions from the particle-hole, transversal particle-hole and particle-particle channel [26]. The last term $h_{abcd}^{\nu_{1}\nu_{2}\nu_{3}\nu_{4}}$ must hence include fully irreducible contributions.

Two-particle vertex.

As for two-particle quantities, one is often interested in vertex functions instead of the Green’s function itself, e.g. when calculating susceptibilities [27, 28] or for diagrammatic extensions of DMFT [29, 30, 26]. The full vertex $F$ is related to the two-particle Green’s function by

[TABLE]

This bears a certain similarity to Eq. (III.2), and it becomes apparent, that the full vertex is given by

[TABLE]

Error propagation.

Assuming that quantities measured in CT-QMC have approximately the same noise amplitude over the whole frequency range, we can conclude that this holds for $\mathfrak{F}$ as well. We furthermore conclude from Eq. (26) that the noise amplitude is rescaled by $\mathcal{G}^{\nu}/G^{\nu}$ in every frequency variable. Since this ratio goes to 1 in the limit of high frequency, we may finally conclude that the noise amplitude of the vertex $F$ in the high-frequency region is identical to the noise amplitude of $\mathfrak{F}$ and thus directly proportional to the error of the Monte Carlo simulation. We hence expect (and will confirm this later) a constant noise level for large Matsubara frequencies.

In contrast, for conventional CT-QMC calculations in the hybridization expansion a strong increase of noise in $F$ with increasing Matsubara frequencies is observed. This yields a too noisy vertex at high frequencies so that approaches to circumvent the calculation of $F$ at large frequencies have been developed, as e.g. replacing the vertex by its high frequency asymptotics [15]. This high noise level of conventional CT-HYB calculations can be understood from Eq. (25). If we assume a constant noise level of the two-particle Green’s function on the left hand side, extracting $F$ by dividing through four Green’s functions increases the error four times by a factor $\sim\nu_{i}$ for large $\nu_{i}$.

III.3 Three-leg vertex

Let us further define the particle-hole three-leg Green’s function as

[TABLE]

where $\nu_{1}$ and $\nu_{2}$ are fermionic Matsubara frequencies. Differentiation with respect to $\tau_{1}$ and $\tau_{2}$ leads to the symmetric improved version

[TABLE]

with the auxiliary definitions

[TABLE]

and

[TABLE]

Since the $R_{2}$ term of the symmetric improved estimator of the two-particle Green’s function essentially contains a threeleg Green’s function, the derivation of the threeleg symmetric improved estimator is completely analogous to the one shown in the last subsection of Appendix B.

For diagrammatic extensions of DMFT, the associated three-leg vertex functions are often particularly interesting, see e.g. Refs. [31, 32, 33, 34, 35]. These are related to the three-leg Green’s function by subtraction of disconnected parts and division by two one-particle Green’s functions. Representative for the variety of definitions, we show the so-called Kernel-2 function in the particle-hole channel (for a definition see [36, 14, 15]), which can be obtained through the following symmetric improved estimator

[TABLE]

IV Implementation

As discussed before, the expressions for e.g. the one- or two-particle Green’s functions derived using the EOM (Eq. (13) and Eq. (III.2)) can be employed in a CT-HYB simulation to obtain results with asymptotically smaller error for high frequencies. In the simulation, we need to get the QMC estimates for the values of each individual term contributing to the symmetric improved estimators (cf. Tab. 1), i.e. of correlation functions consisting of up to 12 operators with up to 4 different imaginary time or 3 different Matsubara frequency arguments with several components corresponding to possible combinations of the discrete quantum numbers of the operators. We accomplish this by performing worm sampling in our CT-HYB w2dynamics[37] program package.

A full introduction to the CT-HYB algorithm is given in Ref. 4, but for the sake of understanding let us briefly recall some main aspects. The starting point is the AIM Hamiltonian Eq. (II). After integrating out the bath degrees of freedom one can evaluate the thermodynamic partition function $Z=\operatorname{Tr}\exp(-\beta H)$ for the impurity by summing over all impurity field configurations. The CT-HYB choice of expanding the exponential in the hybridization turns it into a series of “local” traces. These consist of pairs of impurity operators that evolve in imaginary time according to the local part of the Hamiltonian. The hybridization with the bath is described by a determinant of a matrix that contains, order by order, the hybridization function $\Delta$ where the impurity operators in the local trace are.

Each combination of expansion order, imaginary times, orbitals and spins of the local operators describes one point in the space of partition function configurations $\mathcal{C}_{Z}$. The quantity of interest, e.g. a Green’s function $G_{ab}=\operatorname{Tr}T_{\tau}\exp(-\beta H)d_{a}(\tau_{a})d^{\dagger}_{b}(\tau_{b})/Z$, can be obtained from $\mathcal{C}_{Z}$ either by manipulating a $\mathcal{C}_{Z}$ configuration accordingly in the measurement step or by directly sampling $\mathcal{C}_{G_{ab}}$ configurations. These are like the $\mathcal{C}_{Z}$ configurations, but explicitly contain (in this case) two additional operators (“worm”) that appear in the definition of $G_{ab}$. The worm algorithm consists in sampling both $\mathcal{C}_{G_{ab}}$ and $\mathcal{C}_{Z}$ (for normalization) in one simulation. More details on worm sampling can be found in Ref. 6.

This kind of worm sampling is employed for the calculation of each component $k_{a_{1}\dots a_{n}}$ of each correlation function $k$ needed for the symmetric improved estimator (see Table 1). When we perform a sampling run in the extended configuration space $\mathcal{C}_{Z}\oplus\mathcal{C}_{k_{a_{1}\dots a_{n}}}$ which includes partition function configurations and worm configurations for the specific correlation function, the measurement procedure itself trivially consists of counting samples.

As we have seen, the weight in the worm spaces differs from the weight of a similar $Z$ space configuration in the value of the local trace, but we also explicitly add a suitably chosen weight factor $\eta_{k_{a_{1}\dots a_{n}}}$ to balance the number of steps spent in the current worm space and the partition function space. It is not necessary to perform separate runs per component and quantity, but this also makes it simpler in practice to select the appropriate amount of measurements for the desired target error of the final result (to which all components of all quantities may contribute differently).

In both configuration spaces, we allow all moves that only change operators connected to hybridization events, i.e. in our case pair insertions, pair removals and some global moves. Additionally, worm insertion and removal steps must be employed to change between the two subspaces $\mathcal{C}_{Z}$ and $\mathcal{C}_{k_{a_{1}\ldots a_{n}}}$. For estimators with density-like parts only, such as $\phi_{aabb}^{\omega}$, this should even in practice be enough to ensure ergodic sampling. For other estimators, possible quantum number violations (i.e. configurations with two sequential operators that are zero because of commutation relations) and changes in the energies of states occurring in the local time evolution strongly suppress insertions with large time differences between “compensating” operators (cf. Ref. 6, 38). If the worm operators’ positions could only be changed in the worm insertion step, this would lead to problematically bad statistics for large distances (towards $\Delta\tau=\beta/2$).

Therefore, we introduce further moves that shift or replace some of the estimator’s worm operators analogous to the worm replacement moves[6]. In Ref. 6, these moves transfer the “worm status” from a worm operator to a hybridization operator, i.e. they change which operators are connected with hybridization events (and accordingly only change the bath weight). Since, contrary to Ref. 6, our estimators contain several operators at equal times, the procedure needs to be slightly modified as compared to a simple replacement: After we select one impurity operator connected to a hybridization event at $\tau_{h}$ and one worm operator at $\tau_{w}$ for a replacement, we not only “exchange” them (changing only the bath part of the weight), but also move any other worm operators at $\tau_{w}$ to $\tau_{h}$ to reconstruct the same equal-time object at another position (cf. Fig. 3). This may cause a lower acceptance rate compared to simpler replacements, but especially if the other worm operators are density-like (which may e.g. be the case when performed on a $q$), these moves are reasonably effective.

V Validation

The best way to validate CT-QMC results is benchmarking against exact diagonalization (ED). To this end CT-QMC results were calculated in w2dynamics[37] and, at the same time, all estimators were calculated exactly by evaluating their Lehmann representation for a small Hamiltonian, i.e. an impurity model with one orbital and a discrete bath. The ED results were also used to confirm the validity of Eq. (III.2).

Specifically, we choose a bath with one energy level $\varepsilon_{\mathrm{bath}}=0.5$ and hybridization amplitude $V=0.3$. The chemical potential is set to $\mu=-0.1$, the inverse temperature is chosen to be $\beta=10$ and the local interaction $U=2$. At these parameters, every spin-orbital is on average occupied by 0.307 electrons.

V.1 Self-energy

The quantities $G^{\nu}$, $\xi^{\nu}$ and $\theta^{\nu}$ were evaluated by performing 1.44$\times 10^{9}$ QMC measurements on their respective estimators. Subsequently the one-particle Green’s function and the self-energy were calculated from Eq. (9) and Eq. (13). Fig. 4 shows a comparison of the improved and the symmetric improved self-energy as well as the one obtained directly from the CT-HYB Green’s function as calculated by worm sampling without improved estimators. We note that both “improved” ways to calculate the self-energy suffer from lower precision in the low-frequency regime. However, except for the first few Matsubara frequencies, where precise quantities can be obtained by conventional $Z$-sampling of the Green’s function, the symmetric improved one-particle estimator yields considerably better results. As anticipated in Section III.1, we also observe a better scaling of the error at high Matsubara frequencies.

V.2 Vertex functions

Vertex functions are related to the two-particle Green’s function. In order to assemble the symmetric improved form by Eq. (III.2), one needs to sample all seven occurring estimators. The full vertex $F^{\nu\nu^{\prime}\omega}\equiv F^{\nu(\nu-\omega)(\nu^{\prime}-\omega)\nu^{\prime}}$ can be obtained from the connected part of the two-particle Green’s function by “amputation” of its legs, i. e. division by a product of four Green’s functions, cf. Eq. (25). In the high-frequency case this leads to massive noise amplification, if the two-particle Green’s function is directly computed in QMC. However, as discussed in Sec. III.2, this is healed by symmetric improved estimators. In Fig. 5 we show slices through $F^{\nu\nu^{\prime}\omega}_{\uparrow\uparrow\downarrow\downarrow}$ at two fixed bosonic frequencies $\omega$. In analogy to the one-particle estimator, we get precise results over the whole frequency range, and in particular, there is no increase of noise at high Matsubara frequencies.

V.3 Three-leg vertex

Eq. (III.3) allows us to compute the kernel-2 function $K^{(2),ph,\nu(\nu-\omega)}$ from QMC-estimators that were also used for the full two-particle vertex. In order to judge the improvement introduced by Eq. (III.3), we compute $K^{(2),ph,\nu(\nu-\omega)}_{\uparrow\uparrow\downarrow\downarrow}$ not only in this new way, but also in the way of Ref. 15 from the three-leg Green’s function Eq. (III.3) measured in QMC. In Fig. 6 we compare the result obtained by symmetric improved estimators to the exact result and to the result of the conventional calculation. Notably also here the increase of noise at high Matsubara frequencies is absent.

VI Conclusion

We derived equations of motion for the one- and two-particle Green’s function differentiating symmetrically with respect to all time arguments. With these symmetric improved estimators, we found a way to compute self-energy and vertex functions in CT-HYB without suffering from noise that strongly increases at high Matsubara frequencies as in conventional CT-HYB calculations. In fact, our results rapidly converge towards the exact results, with the exception of the lowest few Matsubara frequencies, where our estimators do not lead to an improvement. For very weak hybridization the self-energy and vertex on the lowest few Matsubara frequencies can actually be calculated with higher accuracy if no improved estimators are used. We hence conclude that in some cases it will be best to combine conventional and improved estimators, using the former for small frequencies and the latter at large frequencies.

Acknowledgements.

We thank Markus Wallerberger, Andreas Hausoel, Tin Ribic, Oleg Janson, and Dominique Geffroy for valuable discussions. This work has been supported by the Vienna Scientific Cluster (VSC) Research Center funded by the Austrian Federal Ministry of Science, Research and Economy (bmwfw), the Deutsche Forschungs Gemeinschaft (DFG) through research unit FOR 1346, the Austrian Science Fund (FWF) through project P 30997-N32, and the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC grant agreement n. 306447 (AbinitioD$\Gamma$A). J.K. was partially supported by the Leibniz society through the Leibniz competition. A.K. and G.S. have been supported by the DFG (through SFB 1170 “ToCoTronics”). G. S. further acknowledges financial support from the DFG through the Würzburg-Dresden Cluster of Excellence on Complexity and Topology in Quantum Matter – ct.qmat (EXC 2147, project-id 39085490). Some of the computational results were obtained using the Vienna Scientific Cluster (VSC). The authors also gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (www.lrz.de).

Appendix A One-particle Symmetric Improved Estimator

The first-order equations (i.e. improved estimators) were already derived in Ref. 16. We nevertheless give a detailed derivation to introduce notation and concepts necessary to derive higher-order estimators (i.e. symmetric improved estimators) in the following. This further sets the stage for the two-particle symmetric improved estimators.

The time derivative of the one-particle Green’s function follows as:

[TABLE]

The commutator between the local Hamiltonian and the impurity-annihilation operator on the right hand side can be calculated as:

[TABLE]

with hybridization $V_{K}^{aa}$ and the anti-symmetrizations of the interaction matrix:

[TABLE]

It is convenient to introduce new operators for the contractions of the interaction matrix with three local operators:

[TABLE]

With these equal-time three-particle operators, the one-particle Green’s function follows as:

[TABLE]

which is the equation of motion for the impurity Green’s function of the AIM.

The mixed impurity-bath Green’s function can be further calculated by applying the above procedure onto the bath operator once again:

[TABLE]

Here, the commutator between the local Hamiltonian and the bath-annihilation operator can be calculated as:

[TABLE]

Thus we can relate the time-derivative of the mixed impurity-bath Green’s function to the impurity Green’s function:

[TABLE]

In order to insert this into Eq. (A), we have to Fourier-transform the above expression with respect to $\tau_{1}$:

[TABLE]

Now the expression on the left-hand side of Eq. (A) can be simplified by integration by parts and we have

[TABLE]

Applying the same Fourier transform and integration by parts also to Eq. (A) allows us to insert the above result and we get

[TABLE]

Now the terms can be regrouped to express the Green’s function:

[TABLE]

This expression can be Fourier-transformed with respect to $\tau_{2}$ by applying $\int_{0}^{\beta}\mathrm{d}\tau_{2}\mathrm{exp}(-i\nu_{2}\tau_{2})$ to both sides of the equation. Considering the definitions of the Green’s function and the function $\xi$, one arrives at Eq. (9).

In order to obtain the symmetric improved estimator, the equal-time two-particle Green’s function is differentiated with respect to the impurity creation operator at time $\tau_{2}$, such that:

[TABLE]

The anti-commutator $\{q_{a},d^{\dagger}_{a}\}$ for diagonal hybridization functions follows as:

[TABLE]

The commutator between the local Hamiltonian and the creation operator is calculated in analogy to Eq.(A) and follows as:

[TABLE]

Thus:

[TABLE]

The mixed bath-impurity expectation value is calculated in analogy to Eqs.(A)-(47) and follows as:

[TABLE]

This allows one to express the mixed bath-impurity expectation value as an impurity expectation value. Again the equation can be made algebraic by Fourier transforming it, but this time with respect to $\tau_{2}$:

[TABLE]

Re-inserting into the (Fourier-transformed) Eq.(56) gives:

[TABLE]

Rearranging gives:

[TABLE]

Together with the definitions of Table 1 we obtain the one-particle symmetric improved estimator Eq. (13):

[TABLE]

Appendix B Two-particle Symmetric Improved Estimator

In the following we derive the two-particle symmetric improved estimators. Again the conventional improved estimators were already derived in Ref. 16. The procedure in deriving the higher order (up to fourth-order) equations is in principle equivalent to the one-particle symmetric improved estimators. Nevertheless, the equations are more involved due to the necessity of considering multiple hierarchies of equations of motions. Repeating derivations (such as the explicit calculation of mixed impurity-bath expectation values) are omitted.

First Order

Applying the time derivative onto the first annihilation operator of the two-particle Green’s function gives:

[TABLE]

Again, Fourier transformation with respect to the first time argument converts the differential equation to an algebraic one, and we have

[TABLE]

$R_{1}$ describes the equal-time contribution due to the derivation of the time-ordering operator:

[TABLE]

At this stage one can already observe that it is possible to apply further (imaginary) time derivatives, or equivalently, hierarchies of equations of motion onto the $R_{1}$ term. The corresponding time-derivatives of the one-particle Green’s function are given by the one-particle symmetric improved estimators and where already calculated explicitly in Appendix A.

Second Order

Applying the time derivative onto the first creation operator of $S_{1}$ gives:

[TABLE]

Fourier-transforming this also with respect to $\tau_{1}$ and plugging it into Eq. (64) yields:

[TABLE]

$R_{2}$ describes the equal-time contribution due to the derivation of the time-ordering operator of $S_{1}$:

[TABLE]

Time derivatives of $R_{2}$ with respect to $\tau_{3}$ and $\tau_{4}$ are shown below.

Third Order

Applying the time derivative onto the remaining annihilation operator of $S_{2}$ gives:

[TABLE]

Fourier-transforming this expression for $S_{2}$ now also with respect to $\tau_{1}$ and $\tau_{2}$, and plugging it into Eq. (70) yields:

[TABLE]

$R_{3}$ describes the equal-time contribution due to the derivation of the time-ordering operator of $S_{2}$:

[TABLE]

Time derivatives of $R_{3}$ with respect to $\tau_{4}$ are shown below.

Fourth Order

Applying the time derivative onto the remaining creation operator of $S_{3}$ gives:

[TABLE]

Fourier-transforming $S_{3}$ with respect to the remaining time arguments and inserting it into Eq. (76) yields:

[TABLE]

where $S_{4,abcd}^{\nu_{1},\nu_{2},\nu_{3},\nu_{4}}=h_{abcd}^{\nu_{1},\nu_{2},\nu_{3},\nu_{4}}/\beta$, as defined in Table 1. $R_{4}$ describes the equal-time contribution due to the derivation of the time-ordering operator of $S_{3}$:

[TABLE]

Unlike the previous terms $R_{1},R_{2},R_{3}$, the above $R_{4}$ cannot be derived any further.

Third and fourth order of $R_{2}$ and $R_{3}$

Looking at equation Eq. (82) one notices that in fact only $R_{4}$ and $S_{4}$ are multiplied by a product of four non-interacting Green’s functions. $R_{1}$ consists of one-particle Green’s functions, which can be expressed by their symmetric improved estimators. However, $R_{2}$ and $R_{3}$ are multiplied only by products of 2 and 3 non-interacting Green’s functions, respectively. This means that we need to express them by their equations of motion recursively.

First, let us take a closer look at the first term of $R_{2}$, Eq. (72). It is very similar to the three-leg Green’s function of Eq. (III.3), the only difference being in the operator ordering and the additional $U$-matrix. The following steps may therefore also serve as a derivation of the threeleg improved estimator.

Taking the time-derivative with respect to $\tau_{3}$ generically yields

[TABLE]

Analogously as before, we express the time derivative by the Heisenberg equation of motion and Fourier-transform the equation with respect to $\tau_{3}$ to make it algebraic. We thus obtain

[TABLE]

The rest-term $R_{21}$ that originates from the derivative of the time-ordering operator, is

[TABLE]

We want to emphasize that due to the bosonic operator $d^{\dagger}d$, a commutator occurs in the first line. After evaluation of the commutator and anti-commutator, we get

[TABLE]

In order to obtain the final form, we also have to apply the equation of motion with respect to $\tau_{4}$. In $R_{21}$ this concerns only the first term. This is, however, just a one-particle Green’s function and we can insert the improved-estimator formula Eq. (9) and we thus get

[TABLE]

Whereas the rest-term $R_{21}$ contains only the one-particle improved estimator and the occupation number, we can differentiate the last term in equation Eq. (B) by its last time argument $\tau_{4}$. The same procedure as above yields now

[TABLE]

The rest-term $R_{22}$ again requires some precaution, since it contains both a commutator and an anti-commutator:

[TABLE]

After evaluation of the commutator and the anti-commutator, we obtain

[TABLE]

We are now able to write down the full expression for $R_{2}$ by combining Eq. (72) with Eq. (B) and Eq. (B). Since in the end we need it in Matsubara frequencies, it is of advantage to perform a Fourier transform with respect to all time arguments already here. Furthermore this allows us to make the equation more compact by using the definitions of Table 1 and the relation $\langle T_{\tau}q_{j}(\tau)d_{k}^{\dagger}\rangle=\langle T_{\tau}d_{k}(\tau)q_{j}^{\dagger}\rangle$, such that we finally arrive at Eq. (21).

For $R_{3}$ we only have to apply the equation of motion with respect to $\tau_{4}$ to the first two terms of Eq. (78). Again we perform a Fourier transform with respect to all time arguments and compactify the expression by the definitions of Table 1. Thus, we obtain Eq. (22).

Appendix C Lehmann representations for two-frequency objects

For objects depending on a single fermionic/bosonic frequency, the standard Lehmann representations for fermionic/bosonic Green’s functions can be used. For objects depending on two fermionic and one bosonic Green’s function, we use the form published in [29]. Additionally, we need it for

[TABLE]

Here, $F_{i}$ ($B$) are fermionic (bosonic) operators, and $\nu$ ($\omega$) are fermionic (bosonic) frequencies. Inserting the eigenbasis of the Hamiltonian and evaluating the integrals analytically, we obtain

[TABLE]

Furthermore, we need

[TABLE]

with two fermionic frequencies $\nu_{1}$, $\nu_{2}$. The result is

[TABLE]