Continuous-time Monte Carlo methods for quantum impurity models

Quantum impurity models describe an atom or molecule embedded in a host material with which it can exchange electrons. They are basic to nanoscience as representations of quantum dots and molecular conductors and play an increasingly important role in the theory of “correlated electron” materials as auxiliary problems whose solution gives the “dynamical mean-field” approximation to the self-energy and local correlation functions. These applications require a method of solution which provides access to both high and low energy scales and is effective for wide classes of physically realistic models. The continuous-time quantum Monte Carlo algorithms reviewed in this article meet this challenge. Derivations and descriptions of the algorithms are presented in enough detail to allow other workers to write their own implementations, discuss the strengths and weaknesses of the methods, summarize the problems to which the new methods have been successfully applied, and outline prospects for future applications.

Figure 1

Diagrammatic representation of configurations $\mathbf{x}=\{(k;{\tau }_1,\dots ,{\tau }_k)\}\in \mathcal{C}$ showing examples with orders $k=0$, 1, 2, 3 and vertices (represented by dots) at times ${\tau }_1,\dots ,{\tau }_3$.

Figure 2

An insertion update (top to bottom) inserting a vertex at time ${\tau }_3$ and the corresponding removal update (bottom to top), removing the vertex at time ${\tau }_3$.

Figure 3

Continuous-time quantum Monte Carlo flow diagram.

Figure 4

Depiction of a third-order term in the weak-coupling expansion. Upper panel: Hubbard interaction vertices denoted by circles. Each $U{n}_{\uparrow }(\tau )n_{\downarrow }(\tau )$ vertex has four operators. Lower panel: One possible contraction of the interaction vertices.

Figure 5

Local updates for the CT-INT algorithm. (a) Starting configuration; (b) insertion of a vertex; (c) removal of a vertex; and (d) shift of a vertex in imaginary time.

Figure 6

Pictorial representation of configurations $\{k;(s_j,{\tau }_j)\}\in \mathcal{C}$ which are sampled by the CT-AUX algorithm. Diagrams for orders 0, 1, 2, and 3. In this algorithm, an auxiliary spin $s_j$ (represented here by the vertices and the direction of the arrows) needs to be sampled in addition to the imaginary-time location ${\tau }_j$ of a vertex.

Figure 7

An insertion update and its corresponding removal update within the CT-AUX algorithm.

Figure 8

Segment configuration of a $k=6$ order term in hybridization expansion of the single-orbital Anderson model. Upper line: spin-up orbital; lower line: spin-down orbital; heavy line: orbital occupied; light line: orbital empty. For each orbital, the length of the black line (occupied orbitals) determines the chemical potential contribution to the weight factor (95). Shaded areas: Regions where both up and down orbitals are filled, so the impurity is doubly occupied. The length of the shaded area enters into an overall weighting factor for the potential energy (Hubbard $U$).

Figure 9

A typical term in the expansion (93): Three “flavors” of fermionic creation and annihilation operators (denoted by filled and empty diamonds, squares, and circles) are placed at times between 0 and $\beta$. In the general case, orbital occupation is not conserved by the local Hamiltonian, so two operators of the same type may follow each other.

Figure 10

An insertion update and its corresponding removal update within the hybridization algorithm.

Figure 11

Updates of the hybridization algorithm as described in the text: (a) original configuration; (b) removal of a segment; (c) shift of an end point of a segment; (d) insertion of an antisegment; (e) removal of an antisegment; (f) removal of another antisegment such that the remaining segment “wraps” around $\beta$.

Figure 12

Hybridization algorithm: Green’s function configuration. A typical configuration for a Green’s function, created by taking the partition function configuration of order $k=3$ and identifying the creation operator at ${\tau }_3^s$ and the annihilation operator at ${\tau }_1^e$ as the Green’s function operators to obtain a Green’s function configuration corresponding to a partition function configuration at one order lower.

Figure 13

Illustration of an insertion of a segment. The diagrams represent the configurations of $\{{\tau }_i\}$ and $\{{\alpha }_i\}$.

Figure 14

Illustration of an order $n=4$ diagram for the Holstein-Hubbard model. Empty (full) circles and squares represent $V^{†}$ ($V$) hybridization events. Dashed lines indicate interactions $K(\tau )$ connecting all pairs of hybridization events. Adapted from 219.

Figure 15

Distribution of perturbation orders in converged single-site DMFT solutions of the half-filled Holstein-Hubbard model with a semicircular density of states with bandwidth $4t$, phonon frequency ${\omega }_0=0.2t$, inverse temperature $\beta t=400$, and values of electron-electron ($U$) and electron-phonon interaction strength ($\lambda$) indicated. In both the insulating ($U/t=6$) and metallic ($U/t=4$) phases, the distribution shifts little as $\lambda$ is increased except near the phase boundary to the bipolaronic phase ($\lambda /t=0.6$, $U/t=4$).

Figure 16

Single-particle spectral function of the one-dimensional Holstein model, computed as a function of frequency using CDMFT cluster dynamical mean-field methods on an $L_c=12$ cluster at filling $n=1/4$ at high temperature (left panel) and low temperature (right panel). The spectra reveal a temperature-dependent line broadening and the appearance at low $T$ of a near-Fermi-level structure associated with the development of intersite correlations. From 11.

Figure 17

Illustration of the Keldysh contour for a CT-AUX study of the Anderson model with interaction quench (top panel) and voltage quench (bottom panel). In an interaction quench starting from $U=0$, the imaginary-time branch of the contour is shifted to $t=-\infty$ and need not be explicitly considered in a weak-coupling Monte Carlo simulation. The arrows represent auxiliary Ising spin variables. The top panel shows a Monte Carlo configuration corresponding to perturbation order $n_{+}=2$, $n_{-}=2$, and the bottom panel shows a configuration corresponding to $n_{+}=3$, $n_{-}=2$, $n_{\beta }=2$. From 220.

Figure 18

Upper panel: Bethe lattice, single-site DMFT, scaling of matrix size with temperature at $U/t=4$ for the Hirsch-Fye, CT-INT, and CT-HYB algorithms. For the Hirsch-Fye algorithm, the resolution $N=\beta U$ has been chosen as a compromise between reasonable accuracy and acceptable speed, while the average matrix size is plotted for the continuous-time solvers. Lower panel: Scaling of matrix size with $U/t$ for fixed $\beta t=30$. The solutions for $U\le 4.5$ are metallic, while those for $U\ge 5.0$ are insulating. The much smaller matrix size in the relevant region of strong interactions is the reason for the higher efficiency of the hybridization expansion method. From 65.

Figure 19

Self-energy $\mathrm{Im}\Sigma (i{\omega }_0)/{\omega }_0$ at the lowest Matsubara frequency ${\omega }_0=\pi T$ as a function of $\beta$ for $U/t=4.0$. The Hirsch-Fye results exhibit large discretization errors, while the continuous-time methods CT-INT and CT-HYB agree within error bars. CT-HYB is particularly suitable for measuring quantities which depend on low-frequency components, such as the quasiparticle weight. From 65.

Figure 20

Sketch of the histogram of the expansion order for flat-histogram sampling. $x$ axis: expansion order $k$; $y$ axis: histogram $h(k)$. Left panel: no flat-histogram sampling. Middle panel: flat-histogram sampling up to one-half of the maximum order. Right panel: flat-histogram sampling up to the maximum contributing order. From 60.

Figure 21

Entropy per site of the half-filled Hubbard model with $t=1$, computed for a periodic $4\times 4$ cluster at different temperatures. From 123.

Figure 22

Sketch of results of applying rotation and block diagonalization operations to the local Hamiltonian. The Hamiltonian in the occupation number basis (shown in the left panel) is sparse but not blocked. A first permutation operation builds blocks according to the occupation number and spin of the local Hamiltonian, leading to a Hamiltonian (shown in the middle panel) which is nearly block diagonal with dense blocks. A second (rotation) matrix further reduces block size by considering rotational and translational invariance of the impurity Hamiltonian, leading to the sparse block structure shown in the right panel.

Figure 23

Sketch of one of the optimizations in the general representation: Four symmetry sectors are drawn, for which $S_z$ and $N$ are different. After the trace of $d_{\uparrow }^{†}$ and $d_{\downarrow }$ is taken only one of the symmetry sectors still contributes. In the implementation, we first identify which symmetry sectors contribute, and then compute the matrix product and trace only for these sectors. Additional symmetries vastly simplify the computation.

Figure 24

Top panel: Binning algorithm—binning of the $k$ operators into $\sqrt{⟨k⟩}$ bins, each having approximately $\sqrt{k}$ elements, reduces the effort of computing the trace after inserting or removing an operator to $O(\sqrt{⟨k⟩})$. Bottom panel: binary tree for the tree algorithm, $O(\log ⟨k⟩)$. From 60.

Figure 25

Metal-insulator phase diagram of paramagnetic two-dimensional Hubbard model in the single-site DMFT approximation, plotted against normalized interaction strength $U_r=(U-U_{\mathrm{MIT}})/U_{\mathrm{MIT}}$ with $U_{\mathrm{MIT}}=9.35t$ for this model. The transition is first order, with the coexistence region indicated by shading. The dashed line indicates the bad metal and bad insulator crossover determined from the condition in which the imaginary part of the self-energy at the few lowest Matsubara frequencies is flat at the crossover value of $U$. From 161.

Figure 26

Top panel: Kinetic energy obtained using the indicated impurity solvers plotted as a function of temperature for the Hubbard model with a semicircular density of states and bandwidth $4t$ and interactions indicated. For $U=4t$ the model is in a strongly renormalized metallic phase, for $U=4.95t$ a low $T$ metal to higher $T$ insulator transition occurs, visible as a jump in kinetic energy at $(T/t{)}^2\approx 0.0007$. From 222. Bottom panel: Doping per spin, $0.5-n$, as a function of chemical potential for $\beta t=400$ and indicated values of $U/t$. At this temperature, the transition at half filling ($\mu ={\mu }_1=U/2$) occurs at $U_c(T)\approx 5.65$. For $U>U_c$ the $n=1$ state is insulating and shifting the chemical potential induces a first-order metal-insulator transition visible as a discontinuity in $n(\mu )$. From 216.

Figure 27

Comparison of the momentum dependence of the self-energy of the two-dimensional Hubbard model with parameters $U/t=4$, $\mu /t=3.1$, and $T/t=0.4$ calculated at the Matsubara frequency ${\omega }_0=\xi =\pi /\beta$, calculated for along the cut $(0,0)-(\pi ,0)-(\pi ,\pi )-(0,0)$ in the first Brillouin zone using a numerically exact diagrammatic Monte Carlo procedure and using CT-AUX simulations of single site and 4-, 8-, 16-, and 32-site DCA DMFT approximations. From 111.

Figure 28

Metal-insulator phase diagram of the paramagnetic phase of the two-dimensional Hubbard model in the plane of temperature $T/t$ and interaction $U/t$ measured relative to the critical end-point value $U_{\mathrm{MIT}}=6.05t$ in the 4-site CDMFT cluster approximation. Band parameters are identical to those used in Fig. 25. Inset: Pie-chart histogram of occupancy probability of the two insulating states at low and high temperatures. From 161.

Figure 29

Solid line: On-site spectral function computed for different momentum sectors by maximum-entropy analytical continuation of QMC data for $U=6t$ and doping $x=0$. Dashed line: Spectral function in the $P=(0,\pi ),(\pi ,0)$-momentum sector. Dotted and dash-dotted lines: $P=(0,\pi ),(\pi ,0)$ and local spectral functions obtained by performing the DCA momentum averages of the standard spin density wave (SDW) mean-field expressions for the Green’s function, with gap $\Delta =1.3t$. From 67.

Figure 30

Left panel: Doping dependence of $P=(0,\pi ),(\pi ,0)$-sector density of states obtained by an analytical continuation of quantum Monte Carlo data obtained from the DCA approximation at $U=5.2t$, temperature $T=t/60$, and dopings $x=0.04$ (solid line), $x=0.08$ (dashed line), and $x=0.15$ (dash-dotted line). The dotted line denotes the noninteracting density of states. Right panel: Evolution of the occupation probabilities with doping at $U=5.2t$ and temperature $T=t/30$. From 67.

Figure 31

Left panel: Brillouin zone partitioning associated with the eight-site cluster DCA approximation with definition of the four inequivalent momentum sectors $A$, $B$, $C$, and $D$. The noninteracting Fermi surface for $t^{\prime }=-0.15t$ and density $n=1$ is indicated by the line. Right panel: Sketch of the paramagnetic state DCA phase diagram of the Hubbard model, calculated for the cluster shown in the left panel at half filling, as a function of interaction strength $U$ and next-nearest neighbor hopping $t^{\prime }$. A Fermi-liquid metal phase (left, shaded area), a sector selective intermediate phase (middle, shaded area) in which the sectors labeled as C are gapped but those labeled as B remain gapless, and a fully gapped insulating phase (right, shaded area) are shown. From 61.

Figure 32

$\beta G(\beta /2)$ calculated in DCA approximation to eight-site cluster for sectors $B$ (full symbols) and $C$ (empty symbols), at $U/t=7$ and $t^{\prime }/t=-0.15$. The strong temperature dependence in the sector $C$ curves arises from the van Hove divergence in the density of states. The crossing points indicate the onset of gapping in the sectors. From 61.

Figure 33

Spectral function $A_{\omega =0,k}$ at the Fermi level calculated for the Hubbard model at $t=0.25$, $t^{\prime }=-0.075$, $U=4.0$, $\beta =80$, and doping $x=0.14$ using the lowest-order momentum-dependent diagram in the dual-fermion method, with analytical continuation performed by polynomial extrapolation from Matsubara frequencies. An anisotropic destruction of the Fermi surface in the pseudogap regime is clearly visible. From 177.

Figure 34

Momentum-integrated spectral function (DOS) of the half-filled single-band Hubbard model calculated using the dual-fermion method. Left panel: Metallic (thin lines) and insulating (heavy line) spectral functions within the coexistence region of the Mott transition $U/t=6.25$, $T/t=0.08$. Right panel: DMFT (dashed line), lowest-order dual-fermion correction (thin line), and ladder dual-fermion (thick line) DOS at $U/t=4$, $T/t=0.19$. The ladder dual-fermion result exhibits a pseduogap of antiferromagnetic origin. From 69.

Figure 35

Phase diagram of the three-band model with semicircular density of states at $\beta t=50$ and $J=U/6$ in the plane of particle density $n$ and interaction strength $U$. The vertical lines indicate the Mott insulating phases at integral values of $n$. The magnetic state is labeled by $P$ (paramagnetic), $F$ (ferromagnetic), and $A$ (two sublattice antiferromagnetic), while the labels $O(N)$ denote the 3 classes of orbital ordering discussed by 31. The heavy dashed line gives the boundary of the non-Fermi-liquid frozen-moment phase discovered by 214. Inset: Hartree-Fock phase diagram for magnetic phases of the same model. Magnetic phase boundaries are indicated by solid lines and orbital ordering boundaries by dashed lines. $OO$ and $OS$ stand for the orbitally ordered and orbitally symmetric phases, respectively. All transitions are second order except the FM-AFM transition and the orbital ordering transitions at $U\gtrsim 12t$ and small $n$. From 31.

Figure 36

Imaginary time dependence of the spin-spin correlation function $⟨S_z(0)S_z(\tau )⟩$ (positive correlation function, full symbols) and orbital correlation function $⟨n_1(0)n_2(\tau )⟩$ (negative correlation function, open symbols) calculated for a three-band model with $U=8t$ and $J=U/6$ using CT-HYB at carrier concentrations $n$ indicated. The convergence of the spin correlations to a value different from zero while the orbital correlation convergence to zero indicates spin but not orbital freezing in the model. From 214.

Figure 37

Phase diagram calculated using CT-J methods with band theory input for the material family $\mathrm{Ce}{X}_2{\mathrm{Si}}_2$. The abscissa $t$ is defined as the Kondo coupling measured relative to the critical Kondo coupling for the $T=0$ magnetic transition. From 138.

Figure 38

The single-particle excitation spectrum $A(\kappa ,\omega )$ for $J=0.3$ and $n_c=0.9$ at (a) $T=0.25$ and (b) $T=0.0025$. The slanted lines represent the noninteracting spectrum $\omega =\kappa -\mu$ which is realized for $J=0$. From 156.

Figure 39

False-color plot of electron spectral function in frequency $\omega$ and scaled momentum $\kappa$ plane for Kondo lattice model with antiferromagnetic coupling $J=0.2$ in (a) the paramagnetic phase at $T=0.035$ and (b) the antiferromagnetic phase at $T=0.010$. From 86.

Figure 40

False-color plot of electron spectral function in frequency $\omega$ and scaled momentum $\kappa$ plane for Kondo lattice model with ferromagnetic interaction $J=-0.2$ in (a) the paramagnetic phase at $T=0.035$ and (b) the antiferromagnetic phase at $T=0.010$. From 86.

Figure 41

Full (all bands) spectral function for the $dpp$ Hamiltonian description of LaOFeAs. Black line: LDA; gray line: $\mathrm{LDA}+\mathrm{DMFT}$ (computed retaining only density-density interactions). From 5.

Figure 42

Histogram of occupation probabilities for each $3d$ atomic state in the DMFT calculation for ${\mathrm{BaFe}}_2{\mathrm{As}}_2$ at $T=150\mathrm{K}$. The states are sorted by total $d$ occupancy and within each manifold of fixed occupancy by energy. From 115.

Figure 43

Momentum-resolved spectral function $A(\mathbf{k},\omega )$ calculated for ${\mathrm{BaFe}}_2{\mathrm{As}}_2$. Inset: angle resolved photoemission spectroscopy (ARPES) intensity. From 22 and 115.

Figure 44

Imaginary-time density-density correlation function $W$ of two-level, two-lead model evaluated using CT-HYB at the midpoint of the imaginary-time interval and normalized to value at $T=0.01t$, as a function of interaction strength with level energies tuned to be equal. Weak $T$ and $\tau$ dependence is seen in the non-Fermi-liquid phase ($U=0.4$, 0.6) and strong $T$ and $\tau$ dependence in the Fermi-liquid phase ($U=0$, 0.2). For details see 211.

Figure 45

Comparison of CT-INT expansion of Eq. (221) with exact diagonalization for an isolated Co atom with $U=1\mathrm{eV}$, $J=0.4\mathrm{eV}$, at inverse temperature $\beta =2{\mathrm{eV}}^{-1}$, for five electrons (main panel) and eight electrons (inset). Solid lines and crosses: full Hamiltonian; dashed lines and plus symbols: model specified by diagonal terms only. From 59.

Figure 46

Left panel: Total DOS of the $3d$ orbital of the Co atom embedded in the Cu matrix. Model parameters: $U=4\mathrm{eV}$, $J=0.7\mathrm{eV}$, and $\beta =10{\mathrm{eV}}^{-1}$ for the 5-orbital impurity with seven electrons. Right panel: Total DOS of the $3d$ orbital of the Co atom embedded in the bulk of Cu, into the 1st layer and Co adatom on the Cu(111) surface. Model parameters: $U=4\mathrm{eV}$, $J=0.7\mathrm{eV}$, and $\beta =10{\mathrm{eV}}^{-1}$ for the 5-orbital impurity with seven electrons. From 59.

Figure 47

Points: Total current $I$ computed using nonequilibrium CT-HYB methods as a function of the bias voltage for a half-filled, spinless, phonon-coupled quantum dot [Eq. (133) with $U=0$] with voltage bias applied symmetrically, ${\epsilon }_d=0$, ${\mu }_{L,R}=\pm V/2$ at temperature $T=\Gamma /5$ at electron-phonon coupling strengths $\lambda$ and oscillator frequencies ${\omega }_0$ indicated. Lines: results of an approximate analytical calculation (49). From 150.

Figure 48

Left panel: Time evolution of the current for different voltage biases ($U/\Gamma =6$, $T=0$). Right panel: Current-voltage characteristics obtained from the long-time limit of the calculated currents. The symbols show Monte Carlo data for $U/\Gamma =4$, 6, 8, and 10, while the lines correspond to fourth-order perturbation theory. In the initial state, the current is given by the steady-state current through the noninteracting dot. At time $t=0$, the interaction is turned on. From 220.

Figure 49

Momentum distribution $n({\epsilon }_k,t)$ for quenches from $U/W=0$ to $U/W=0.75$ (left panel) and $U/W=1.25$ (right panel), for the Hubbard model with semicircular density of states and bandwidth $W$. $t$ is measured in units of the quarter bandwidth $W/4$. From 43.