Improved estimators for the self-energy and vertex function in hybridization-expansion continuous-time quantum Monte Carlo simulations

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 numerical renormalization-group 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.

Figure 1

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 areas) determines the interaction energy.

Figure 2

Diagrammatical representation of the Green's function $G_{ab}$ in terms of the self-energy $\Sigma$, Eq. (12). Thick lines denote fully dressed propagators; thin lines denote bare propagators.

Figure 3

The $3!$ ways of connecting hybridization lines in a configuration with three segments for a timeline of a given flavor. The sign of each diagram is indicated. For nondiagonal hybridization, the lines also connect segments with different flavors.

Figure 4

Self-energy as a function of Matsubara frequencies for the half-filled Hubbard model on the Bethe lattice (bandwidth $4t$) as calculated from Dyson's equation to illustrate the noise problem. The parameters are the same as in Ref. 12, $U/t=4$ and $T/t=1/45$. 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 ${\lim }_{\nu \to \infty }\Sigma \left(i\nu \right)=U^2\langle n\rangle (1-\langle n\rangle )/\left(i\nu \right)$ is shown for comparison.

Figure 5

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.

Figure 6

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 us 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).

Figure 7

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 8

Diagrammatical representation of the two-particle Green's function ${\chi }_{abcd}$ in terms of its connected part ${\chi }_{abcd}^{\text{con.}}$ and its disconnected part ${\chi }_{abcd}^0$ and definition of the vertex function $\gamma$,Eqs. (23, 24, 25). The lines with arrows denote fully dressed propagators.

Figure 9

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 $T/t=1/10$ for fixed ${\nu }^{\prime }$ and bosonic frequencies $\omega =0$ (circles), $4\pi T$ (triangles), and $20\pi T$ (squares). The parameters are otherwise the same as in Fig. 7. Solid lines show the exact diagonalization data and open symbols show the results from the improved Monte Carlo measurement (measured on the Matsubara axis). The usual Matsubara axis measurement based on Eqs. (22, 23, 24, 25) is shown by crosses. The two results are hardly distinguishable on this scale.

Figure 10

Real part of the spin-up-up component of the reducible vertex function ${\gamma }_{\nu {\nu }^{\prime }\omega }$ computed from the connected part shown in Fig. 9. The standard computation [Eqs. (22, 23, 24, 25), crosses connected by dashed lines] is clearly numerically unstable for large values of ${\nu }_n$, while the improved estimators yield the proper asymptotic behavior. The improved estimators computed from Eqs. (25) and (26) (closed symbols) and from Eqs. (25) and (31) (open symbols) yield similar results. They also yield more accurate results even before the vertex approaches its asymptotic behavior and for small frequencies ${\nu }_n$ (clearly visible for $\omega =20\pi T$).

Figure 11

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, 23, 24, 25) (note that the measurement in the Legendre basis leads to systematic deviations rather than the irregular noise observed in the Matsubara measurement).

Figure 12

Phase diagram in the space of density $\langle n\rangle$ and interaction $U$ of the two-orbital Hubbard model on the Bethe lattice with $J/U=1/6$ and $T/t=0.02$. 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, 17 have been performed.

Figure 13

Imaginary time dependence of the spin-spin-correlation function for $U/t=8$ and densities $\langle n\rangle =1.90,1.72,1.53,1.29,1.10$, and $0.72$ (from top to bottom). For small densities, in the Fermi-liquid regime, the spin-spin-correlation at $\tau =1/\left(2T\right)$ is proportional to $T$ and hence approaches zero as $T\to 0$, while for high densities the spin-correlation function approaches a constant, indicating the presence of frozen moments.

Figure 14

Ratio of the spin-spin-correlation function $C_{\text{SS}}(\tau =1/\left(2T\right))$ for temperatures $T=0.02t$ and $0.01t$ as a function of density. Vertical lines mark the approximate location of the transition, where the self-energy is consistent with a ${({\nu }_n/t)}^{\alpha }$ behavior with $\alpha \sim 0.5$. On the Fermi-liquid side, obviously $C_{\text{SS}}(\tau =1/\left(2T\right))\sim T^2$, while $C_{\text{SS}}(\tau =1/\left(2T\right))=\mathrm{const}.$ in the frozen moment phase.

Figure 15

Imaginary part of the self-energy for various fillings for the two-orbital model at $U/t=8$, obtained using the improved estimator measured in the Legendre basis. The positions of the corresponding parameters are marked by circles in the phase diagram in Fig. 12. The solid line is a fit proportional to ${({\nu }_n/t)}^{\alpha }$, with $\alpha =0.49$.

Figure 16

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 the Matsubara basis. The positions of the parameters in the phase diagram 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 are well described by the Fermi-liquid form $\Sigma \left(\omega \right)\sim (1-1/Z)\omega -i\mu {\omega }^2$ over a wide energy range (thin solid lines), while the data obtained from the two other methods exhibit spurious features.

Figure 17

Real and imaginary parts of the reducible vertex function ${\gamma }_{\nu {\nu }^{\prime }\omega }$ for the two-orbital Hubbard model for ${\nu }^{\prime }=-\pi T$ and $\omega =2\pi T$ across the spin freezing transition. The left panels show intraorbital 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 $\langle n\rangle =0.72,1.10$, and $1.72$, respectively. The position of these parameters in the phase diagram 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.

Figure 18

Phase diagram in the space of crystal-field splitting ${\Delta }_{\text{cf}}$ and interaction strength $U$ of the half-filled two-orbital Hubbard model on the Bethe lattice with bandwidth $4t$, $J/U=1/4$ and $T/t=0.02$. The parameters for which the vertex function in Fig. 20 has been calculated are marked by stars.

Figure 19

Total moment $\langle S_z^2\rangle$, moment of orbital 1, $\langle S_{1z}^2\rangle$ and density of orbital 1, $\langle n_1\rangle$ as a function of $U$ for various values of the crystal-field splitting ${\Delta }_{\text{cf}}$. From left to right (and top to bottom) the curves correspond to the values ${\Delta }_{\text{cf}}=0,0.5,1,1.5,2,2.5,3,3.5$.

Figure 20

Real and imaginary parts of the reducible vertex function ${\gamma }_{\nu ,{\nu }^{\prime },\omega }$ for the two-orbital Hubbard model for ${\nu }^{\prime }=-\pi T$ and $\omega =2\pi T$ at ${\Delta }_{cf}/t=3.0$ 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 $U/t=6.0,7.0$, and $8.0$.