# Continuous-time auxiliary field Monte Carlo for quantum impurity models

## Abstract

We present a continuous-time Monte Carlo method for quantum impurity models, which combines a weak-coupling expansion with an auxiliary-field decomposition. The method is considerably more efficient than Hirsch-Fye and free of time discretization errors, and is particularly useful as impurity solver in large cluster dynamical mean field theory (DMFT) calculations.

###### pacs:

71.10.Fd###### pacs:

02.70.Ss###### pacs:

71.27.+a###### pacs:

71.30.+hLattice fermion models (Hubbard model, etc.) Quantum Monte Carlo methods Strongly correlated electron systems; heavy fermions Metal-insulator transitions and other electronic transitions

## 1 Introduction

The development of efficient numerical methods for solving quantum impurity models has been driven in recent years by the success of dynamical mean field theory (DMFT) [1, 2, 3] and its extensions. DMFT is an approximate framework for the study of fermionic lattice models, which replaces the lattice by a quantum impurity embedded in a self-consistent bath. Both cluster-extensions of DMFT[4, 5, 3, 6, 7] and realistic electronic structure calculations, which combine DMFT with band structure methods [3], involve multi-site or multi-orbital impurity models (e.g. for - and -electron systems), whose solution is computationally expensive and in practice the bottleneck of the calculations. In order to facilitate progress in this field, it is therefore important to develop fast and accurate impurity solvers.

Until recently, the Hirsch-Fye auxiliary field method[8] has been the only Quantum Monte Carlo impurity solver used in DMFT. It suffers from one major drawback: it requires a discretization of the imaginary time interval into a large number of time slices and therefore the calculation of determinants of matrices (for models with sites and on-site interactions only), which is computationally expensive. Furthermore, this discretization introduces a systematic error which needs to be dealt with (in principle) through tedious extrapolations .

Important progress was achieved recently with the development of continuous-time impurity solvers, which are based on the stochastic sampling of a diagrammatic expansion of the partition function. These methods do not suffer from time discretization errors and allow the simulation of models with more general interactions. The first continuous-time impurity solver was proposed by Rubtsov et al. [9], who expanded the partition function in the interaction terms and used Wick’s theorem. Another powerful and flexible diagrammatic solver for small impurity problems, based on a diagrammatic expansion in the impurity-bath hybridization, has been proposed in Refs. [10, 11, 12]. Since this method perturbs around an exactly solved atomic limit, it is particularly efficient at moderate and strong interactions [13]. While the sign problem in this algorithm is less severe than in Hirsch-Fye or in the weak-coupling continuous-time method, the computational effort scales exponentially with the number of sites and orbitals, making it difficult or impossible to solve clusters with eight or more sites. As a consequence, Hirsch-Fye is currently still considered the method of choice for large cluster DMFT computations.

In this paper, we present a new continuous-time impurity solver which combines a weak-coupling expansion with an auxiliary field decomposition, and which was inspired by the work of Rombouts et al. [14] for lattice models. Our method is formally similar to the Hirsch-Fye algorithm, but as a weak-coupling solver performs comparable to Rubtsov’s method.

## 2 Method

We present the algorithm for the single-impurity model corresponding to the DMFT solution of the one-band Hubbard model. The extension to multi-band or cluster models with density-density coupling is straightforward and will be briefly discussed at the end of this section.

The partition function for the impurity model can be written as a path integral over Grassman variables and , , with effective action

(1) |

Here, and is related to the “conventional” non-interacting Green’s function of Ref. [1] by the expression , which means that the chemical potential is shifted by and for .

In order to closely follow the standard derivation of the Hirsch-Fye algorithm (see e.g. Ref. [1]), we switch to the Hamiltonian formulation

(2a) | |||||

(2b) | |||||

(2c) |

where is the Gaussian term containing both the impurity () and the bath () degrees of freedom. Following Rombouts et al. [14] we introduce a constant , express the partition function in an interaction representation,

(3) |

and expand the time ordered exponential in powers of (dropping the irrelevant factor ):

(4) |

We then decouple the interaction terms as follows [14]:

(5a) | ||||

(5b) |

Expressions (5a) and (5b) are valid for arbitrary (complex) parameters . If , is real and the expansion parameter is positive. After the decoupling, the partition function is of the form

(6a) | |||

(6b) |

with for and . is very similar to the expression for the partition function in the Hirsch-Fye algorithm after the Trotter approximation, see for example Eq. (117) of Ref. [1], except that the time arguments of the auxiliary spins are not regularly spaced on . Indeed, one can straightforwardly generalize the calculation in Ref. [1] to rewrite () as:

(7) | |||||

(8) | |||||

(9) |

with the notations , and for , .

While we tried to emphasize in our derivation the similarities to the Hirsch-Fye algorithm, let us note at this point also the essential differences between Hirsch-Fye and our continuous-time auxiliary field method (CT-AUX): i) CT-AUX is based on a weak-coupling expansion, not a Suzuki-Trotter decomposition of the partition function; ii) the auxiliary fields in CT-AUX originate from Rombout’s decoupling formula (5a). In particular, CT-AUX does not require any time discretization. The number and position of auxiliary spins on the imaginary time interval is arbitrary and changes constantly during the simulation.

The formulae above are easily generalized for cluster and multiorbital DMFT problems with density-density interactions by performing a similar expansion for all the interaction terms. For clusters of size with local density-density interaction (relevant e.g. for cluster DMFT approximations of the Hubbard model), expression (5b) for remains unchanged and other formulas, like Eq. (7), can be generalized straightforwardly by replacing the time and spin indices by (time, site) and (spin, site) multi-indices respectively. The multiplicative factor dropped from the partition function is in this case.

### 2.1 Sampling procedure

Our algorithm samples time ordered configurations consisting of spins at times with weight

(10) |

For ergodicity it is sufficient to insert/remove spins with random orientation at random times.

The detailed balance condition can be implemented as follows. Assuming that we pick a random time in the interval and a random direction for this new spin ()), and propose to remove it with probability , we get

(11) |

The matrices are stored and manipulated using fast update formulas analogous to those of Refs. [9, 10]. When inserting a spin we add a new row and column to . Following the notation of Ref. [15], we define the blocks (omitting the index until the end of this section)

(12) |

where , , denote , , matrices, respectively, which contain the contribution of the added spin. The determinant ratio needed for the acceptance/rejection probability is then given by

(13) |

As we store computing the acceptance/rejection probability of an insertion move involves a matrix-vector multiplication followed by an inner product, i.e. an operation. If a move is accepted, a rank one update is performed to compute the new matrix out of , and :

(14a) | ||||

(14b) | ||||

(14c) | ||||

(14d) |

### 2.2 Measurement of the Green’s function

The main observable of interest in the simulations is the Green’s function . First, let us note from (6) that one can add two additional “non-interacting” spins at any fixed times and (we denote with a tilde the corresponding matrices of size ). is then given by an expression similar to Eqs. (6), with an insertion of and at the corresponding times. We can again use the standard Hirsch-Fye formula for the discretized Green function (Eq. (118) of Ref. [1]) to obtain

(15) |

with Since , a simple block calculation leads to

(16) |

In order to compute the Green’s function, one can not just accumulate its values at the discrete times of the auxiliary spins, since the are correlated. Rather, the Green’s function is accumulated using Eqs. (2.2) and (2.2):

(17) | |||

(18) | |||

(19) |

where we have used translational invariance, set , and denoted the Monte Carlo average with angular brackets (our convention is for ). Hence, we measure only the quantity , which we bin into fine bins. After the simulation is completed, the Green’s function is constructed using Eq. (17).

Note that the Dyson equation

(20) |

implies that this procedure amounts to accumulating . Besides the higher efficiency with respect to the direct accumulation of the Green’s function, an important advantage of such a measurement is the reduction in high-frequency noise by the multiplication with (see also Ref. [16] for similar ideas in the NRG-DMFT context).

### 2.3 Four point functions

Four point correlation functions can also be computed in a similar way as in Hirsch-Fye using the fact that for a fixed auxiliary spin configuration the problem is Gaussian and Wick’s theorem can therefore be used together with Eq. (2.2). Thus the problem reduces to the accumulation of the determinant of a matrix

(21) |

with defined in Eq. (19). If only a few correlation functions are measured, Eq. (21) is best evaluated directly during the simulation. If many or all correlation functions have to be measured at time points and the size of is comparatively small, it is advantageous to accumulate only and and reconstruct the correlation function at the end of the computation. Indeed, while binning the latter expression is in memory, it is only computationally (using time translation invariance).

### 2.4 Role of the expansion parameter

The average perturbation order is related to the parameter , potential energy and filling by

(22) |

This expression is obtained by applying the operator to both in its original form (3) and to (6), including the factor dropped after Eq. (3) (see also Ref. [14]). In the case of the weak-coupling algorithm [9], , where is the small parameter which must be introduced to reduce the sign problem. Hence, the perturbation order in the continuous-time auxiliary field method grows linearly with (see inset of Fig. 1) and .

Figure 1 shows the perturbation orders for the two methods as a function of . For these small values of and , the perturbation orders are essentially identical. Both weak-coupling methods scale roughly linearly with , with a kink visible at the Mott critical value. It also follows from Eq. (22) that the perturbation order is essentially linear in the inverse temperature .

Similar to the weak-coupling expansion parameter [9], the parameter can be freely adjusted. While a larger yields a larger expansion order, it also reduces the value of (see Eq. (5)). This makes it easier to flip auxiliary spins. Therefore the auxiliary spins have less tendency to polarize for larger . In practice, however, -values of order turned out to be adequate. Although we found that the sign problem improves slightly with larger , this small gain is more than compensated by the increase in computational cost at larger values of .

## 3 Comparison with other QMC methods

To compare to other methods we implemented single-site, as well as 4 and 8 site cluster calculations in the dynamical cluster approximation (DCA)[4, 17], and expect similar results for other cluster schemes such as cellular dynamical mean field theory (CDMFT) [7]. The upper panel of Fig. 2 shows a typical real space cluster Green’s function for a 4-site DCA calculation. The CT-AUX results are identical to the other QMC results, showing the accuracy of the new approach. The lower panel of Fig. 2 shows Green’s functions for an 8-site cluster () obtained from a converged in 8 CPU hours on a 1.6 GHz Opteron 244. Symbols show the result for Hirsch-Fye with 80 time slices, and the lines indicate the result from CT-AUX, measured at 500 time points.

As a continuous time method CT-AUX has a definite advantage over Hirsch-Fye QMC, since it removes the necessity of careful extrapolations to the continuous time limit. However, in order to be a useful replacement for Hirsch-Fye in practice, CT-AUX has to satisfy two requirements: i) The average expansion order , which determines the complexity of the calculation (), has to be smaller than the number of times slices required in Hirsch-Fye (close to the continuous limit, where extrapolation is meaningful); ii) The sign problem must not be worse than in previous algorithms.

### 3.1 Expansion order

First we compare the expansion order to Hirsch-Fye for a high temperature 2x2 DCA calculation in Fig. 3. In Hirsch-Fye the number of time slices was fixed a priori using the optimistic criterion which corresponds to , where is the size of the time slices – just barely in the region of validity of the Trotter approximation underlying the Hirsch-Fye method. CT-AUX with its roughly two times lower average perturbation order is much more efficient since both algorithms scale like the cube of the matrix size. In the 8-site cluster simulation of Fig. 2, the Hirsch-Fye algorithm with 80 time slices had to update matrices of size 640, while CT-AUX merely had to operate on matrices of average size 136. This means that CT-AUX allows to reach substantially lower temperatures, even with modest computational resources.

To make the comparison with Hirsch-Fye more precise and get rid of the arbitrariness of the choice of the number of times slices in Hirsch-Fye we have reproduced in Fig. 4 the self-energy calculation presented in Fig. 15 of Ref. [1], where the same single-site calculation was performed with Hirsch-Fye (32, 64 and 128 time slices) and with exact diagonalization (ED). Even with 128 times slices, the Hirsch-Fye results still have substantial discretization errors while CT-AUX produces the exact result (comparable to the ED result) with only spins. This shows that CT-AUX is indeed much more efficient than the Hirsch-Fye method: not only does it compute the numerically exact result directly, but it does so using significantly less auxiliary spins. This is due to the fact that CT-AUX, like Rubtsov’s method, is based on a weak-coupling expansion (see Fig. 1 and Ref. [13] for a comparison of the weak-coupling method with Hirsch-Fye).

### 3.2 Sign problem

For single site impurity models, there is no sign problem since the proof of Ref. [18] can be extended to CT-AUX. For cluster calculations, as and is increased, the sign becomes smaller than one. However, for the unfrustrated plaquette at temperatures down to we did not observe a significant sign problem (). In order to produce a severe sign problem at high temperature, we frustrated our plaquette with a hopping along the diagonal. For the almost triangular case the Hirsch-Fye method, weak-coupling method and our solver exhibit a sign problem that becomes increasingly severe as the interaction strength is increased or the temperature lowered. For , the average sign is less than , as seen in Fig. 5, making it difficult to access temperatures below . Remarkably, the average signs in CT-AUX, Hirsch-Fye and the weak-coupling algorithm are almost identical.

Since one of the likely applications of the CT-AUX method is the solution of large DMFT clusters (not accessible to the hybridization expansion solver), we performed a similar study for an 8-site cluster, with a similar conclusion. The sign problem at a reference point on the eight site Betts cluster ( , )[19] turned out to be the same in Hirsch-Fye as in our new algorithm (Fig. 5).

## 4 Conclusion

We have presented a continuous time impurity solver based on a weak-coupling expansion of the partition function and an auxiliary field decomposition of the interaction terms. The algorithm relies on fast local updates of auxiliary Ising spin variables, whose number and position are not fixed. As a continuous time solver, our method does not suffer from the deficiencies of the Hirsch-Fye algorithm and its variants [13]. In particular, it does not require multiple runs for several discretizations of the imaginary time interval and subsequent extrapolations. Moreover, it requires fewer auxiliary field variables than a standard Hirsch-Fye calculation. This translates into substantial gains in computational efficiency, especially for large clusters and at low temperature.

For all regions of parameter space considered the sign problem is approximately the same in the weak coupling, Hirsch-Fye and CT-AUX algorithms. Further investigation is however needed to determine if this is the case in all regions of parameter space and for all cluster geometries. We expect the new solver to be particularly useful in the simulation of large clusters, and to completely replace the Hirsch-Fye algorithm.

###### Acknowledgements.

The calculations have been performed on the Hreidar beowulf cluster at ETH Zürich, using the ALPS-library,[20] and our cluster simulations employed a DCA self-consistency loop implemented by S. Fuchs. PW acknowledges support from NSF-DMR-040135.### References

- \NameGeorges A., Kotliar G., Krauth W. Rozenberg M. J. \REVIEWRev. Mod. Phys. 68199613.
- \NameGeorges A. Krauth W. \REVIEWPhys. Rev. Lett. 6919921240.
- \NameKotliar G., Savrasov S. Y., Haule K. \etal \REVIEWRev. Mod. Phys. 782006865.
- \NameHettler M. H., Tahvildar-Zadeh A. N., Jarrell M. \etal \REVIEWPhys. Rev. B 581998R7475.
- \NameMaier T., Jarrell M., Pruschke T. Hettler M. H. \REVIEWRev. Mod. Phys. 7720051027.
- \NameLichtenstein A. I. Katsnelson M. I. \REVIEWPhys. Rev. B 6220009283.
- \NameKotliar G., Savrasov S. Y., Pálsson G. Biroli G. \REVIEWPhys. Rev. Lett. 872001186401.
- \NameHirsch J. E. Fye R. M. \REVIEWPhys. Rev. Lett. 5619862521.
- \NameRubtsov A. N., Savkin V. V. Lichtenstein A. I. \REVIEWPhys. Rev. B 722005035122.
- \NameWerner P., Comanac A., de’ Medici L. \etal \REVIEWPhys. Rev. Lett. 972006076405.
- \NameWerner P. Millis A. J. \REVIEWPhys. Rev. B 742006155107.
- \NameHaule K. \REVIEWPhys. Rev. B 752007155113.
- \NameGull E., Werner P., Millis A. Troyer M. \REVIEWPhys. Rev. B 762007235123.
- \NameRombouts S. M. A., Heyde K. Jachowicz N. \REVIEWPhys. Rev. Lett. 8219994155.
- \NameVetterling W. T., Flannery B. P., Press W. H. Teukolski S. A. \BookNumerical Recipes in FORTRAN - The Art of Scientific Computing - Second Edition (University Press, Cambridge) 1992.
- \NameBulla R., Hewson A. C. Pruschke T. \REVIEWJournal of Physics: Condensed Matter 1019988365.
- \NameHettler M. H., Mukherjee M., Jarrell M. Krishnamurthy H. R. \REVIEWPhys. Rev. B 61200012739.
- \NameYoo J., Chandrasekharan S., Kaul R. K. \etal \REVIEWJournal of Physics A Mathematical General 38200510307.
- \NameMaier T. private communication (2007).
- \NameAlbuquerque A., Alet F., Corboz P. \etal \REVIEWJournal of Magnetism and Magnetic Materials 31020071187.