Diagrammatic Monte Carlo for dual fermions

We introduce a numerical algorithm to stochastically sample the dual fermion perturbation series around the dynamical mean field theory, generating all topologies of two-particle interaction vertices. We show results in the weak and strong coupling regime of the half-filled Hubbard model in two dimensions, illustrating that the method converges quickly where dynamical mean field theory is a good approximation, and show that corrections are large in the strong correlation regime at intermediate interaction. The fast convergence of dual corrections to dynamical mean field results illustrates the power of the approach and opens a practical avenue towards the systematic inclusion of nonlocal correlations in correlated materials simulations. An analysis of the frequency scale shows that only low-frequency propagators contribute substantially to the diagrams, putting the inclusion of higher order vertices within reach.

Figure 1

Left panel: Diagram notation for a dual fermion vertex. Two incoming lines with momenta/energies $\nu +\mathrm{\Omega }$ and ${\nu }^{\prime }$ and spin $\sigma$ and ${\sigma }^{\prime }$ scatter into two outgoing lines with momenta ${\nu }^{\prime }+\mathrm{\Omega }$ and $\nu$. A momentum $\mathrm{\Omega }$ is transferred from right to left. Right panel: lowest order diagram. For clarity we omit the arrow inside the vertex.

Figure 2

Auxiliary “worm” diagrams. Left panel: a worm tail is attached to the spin-$\sigma$ line and carries away momentum $\delta$. Right panel: a worm head, attached to the spin ${\sigma }^{\prime }$ line, carries in momentum $\delta$ (right). Note that momentum is conserved at each vertex and head or tails can be attached to either line.

Figure 3

Two complementary ways of inserting an auxiliary “worm” diagram line carrying momentum $\delta$ into a second order Luttinger-Ward diagram.

Figure 4

“Worm move” updates: starting from a worm configuration as in the left panel of Fig. 3, the worm head is (top left) moved across a vertex, changing its momentum transfer from $\mathrm{\Omega }$ to $\mathrm{\Omega }-\delta$, or (top right) its tail is moved across a vertex, changing its momentum from $\mathrm{\Omega }$ to $\mathrm{\Omega }+\delta$. Bottom panel: A worm head is moved along a propagator line, changing its momentum from ${\nu }^{\prime }$ to ${\nu }^{\prime }-\delta$.

Figure 5

Updates changing diagram topology: A worm connecting two diagrams reconnects their outgoing vertex lines.

Figure 6

Updates changing diagram expansion order: As read from left to right: update changing a worm space configuration with $n$ vertices into a Luttinger-Ward diagram configuration with $n+1$ vertices, by inserting a new vertex with momentum transfer $\delta$. As read from right to left: update changing a Luttinger-Ward diagram configuration with $n+1$ vertices into a worm configuration with $n$ vertices by removing a vertex.

Figure 7

Example of a fermion loop identified by randomly selecting a propagator and following one of its fermion lines until it reaches the original vertex. The “spin flip” update then proceeds by flipping the spin of all propagators and vertices along this loop.

Figure 8

Contribution to the real part of the Luttinger-Ward functional of the dual fermion expansion as a function of diagram expansion order. Hubbard model, half filling, $T/t=0.5$, and $U/t=4.0,8.0,10.0$.

Figure 9

Probability that a dual diagram at order 2, 3, 4, or 5 contains propagators with frequency $i{\omega }_n$, normalized to the probability of a diagram at the same order containing a propagator at the lowest Matsubara frequency. Data for $T/t=0.5$ and $U/t=8.0$.

Figure 10

Probability that a dual diagram at order 2, 3, 4, or 5 contains propagators with frequency $i{\omega }_n$, normalized to the probability of a diagram at the same order containing a propagator at the lowest Matsubara frequency. Data for $U/t=4.0$ (upper row) and $U/t=10.0$ (lower row), for $T/t=0.5$ (left panel) and $T/t=0.2$ and $T/t=0.26$ (right panel).

Figure 11

Real (upper row) and imaginary (lower row) part of the ladder dual self-energy (green lines) and real part of the self-energy from all diagrams sampled up to fifth order (red lines) plotted along the main axes of the Brillouin zone at the lowest Matsubara frequency, discretized on a $16\times 16$ grid of momenta and obtained at the coupling strengths $U/t=4$ (left column), $U/t=8$ (middle column), and $U/t=10$ (right column) and at $T/t=0.5$.

Figure 12

Real part of the lattice self-energy plotted along the main axes of the Brillouin zone at the lowest Matsubara frequency, discretized on a $16\times 16$ grid of momenta and obtained for the weak coupling strength $U/t=4$ at $T/t=0.5$, and evaluated up to second order (triangles, green lines), third order (diamonds, blue lines), and fourth order (circles, red lines). Black squares: Comparison to the piecewise constant real self-energy obtained from a dynamical cluster approximation calculation of the same system obtained on a 64-site cluster.