Diagrammatic Monte Carlo procedure for the spin-charge transformed Hubbard model

Using a dual representation of lattice fermion models that is based on spin-charge transformation and fermionization of the original description, I derive an algorithm for diagrammatic Monte Carlo simulation of strongly correlated systems. This scheme allows eliminating large expansion parameters, as well as large corrections to the density matrix that generally prevent diagrammatic methods from being efficient in this regime. As an example, I compute the filling factor for the Hubbard model at infinite on-site repulsion and compare the results to controllable data obtained from numerical linked-cluster expansion. I find excellent agreement between the two methods, as well as rapid convergence of the diagrammatic series. I also report results for the momentum distribution and kinetic energy of the electrons.

Figure 1

The elementary diagrammatic objects that emerge from the model (17): (a) Terms originating from hopping electrons and (b) projection terms used to enforce the restricted Hilbert space. The line types are described in (c). This work uses a description based on imaginary time $\tau$ and momentum $\overline{k}$, although it is possible to work with any combination of real-space/momentum and imaginary-time/Matsubara frequencies. Thus, all lines are functions of $({\tau }_2-{\tau }_1,\overline{k})$.

Figure 2

Basic updates: “N” is an artificial configuration that exists outside the diagrammatic space and is used for normalization. It is connected to the set of first-order diagrams through leaving or entering the normalization sector [(a) $\iff$ (b)]. Adding or removing a vertex [(b) $\iff$ (c)] changes the diagram order by one and either creates or destroys a worm. Commuting operators [(c) $\iff$ (d)] changes the topo- logy while preserving the diagram order. The worm absorbs the momentum difference $k_w=k_2-k_1$ between the commuted particle lines so that momentum is preserved at every vertex. Moving one of the worm ends [(d) $\iff$ (e)] along a particle line changes its momentum by $k_w$. Adding or removing a worm [(e) $\iff$ (f)] connects configurations where a worm closes on itself with a worm-free state. Finally, we allow translating a vertex or pair of vertices in imaginary time (not depicted).

Figure 3

Schematic view of the measurement process. A measuring line is chosen stochastically. In (a) a spinful fermion (marked by an X) is selected, so the measuring line is given by $G_{\sigma }(\tau ,\overline{k})$. The tagged line is removed in (b), and the diagram is interpreted as a contribution to the self-energy ${\mathrm{\Sigma }}_{\sigma }(-\tau ,\overline{k})$. In (c), a hopping line is selected, and the interpretation in (d) is a contribution to the polarization operator ${\mathrm{\Pi }}_t(-\tau ,\overline{k})$. If a measurement is attempted while in the normalization sector, then a number is added to a special normalization bin (denoted “N” in Fig. 2). The properly normalized self-energy or polarization operator is obtained by dividing by this bin.

Figure 4

Projected Green's function for the holes ${\tilde{G}}_{h,\sigma }(\tau ,\overline{k})$ [Eq. (28)] at the high-symmetry points $\overline{k}=\{0,0\}$ and $\overline{k}=\{\pi ,\pi \}$. The model parameters are $\mu =2,U=\infty$, and the temperatures are shown in the legend. The expansion order is $O=3$ in all cases.

Figure 5

Filling factor [Eq. (31)] as a function of temperature (in units of hopping) for the 2D Hubbard model when $U=\infty$ and $\mu /t=2$. The colored circles show results obtained through SCT-BDMC, with expansion orders ranging from $O=2$ to $O=4$. The solid and dashed lines show results obtained by Khatami et al. from numerical linked-cluster expansion (NLCE) [42]. Data are given both for bare expansions to orders 10 and 11 and after Wynn resummation. At $T/t=2$ there is little discrepancy between techniques or expansion orders, indicating rapid convergence with either method. With decreasing temperature, bare NLCE and the lowest-order expansions start to diverge. Bare NLCE remains reliable down to $T/t\approx \sqrt{1/2}$, and these results are reproduced by SCT-BDMC within error bars. When $T/t<\sqrt{1/2}$, a comparison can be made only to data obtained via Wynn resummation. In this region diagrammatic expansion to orders 3 and 4 start to diverge, indicating that higher-order terms are required to see convergence. Still, expansion to fourth order gives results in good agreement with the resummed data.

Figure 6

Electron density as a function of momentum $n\left(k\right)$ obtained from Eq. (30) at temperatures $1/8\le T/t\le 2$ and with $U=\infty ,\mu =2t$, plotted along the high-symmetry axes shown in the inset. The data correspond to the highest expansion order shown in Fig. 5 at every temperature. The error estimates were obtained by conducting multiple simulations and taking $\mathrm{\Delta }n=\text{max}\left[\right|n_i\left(\overline{k}\right)-n_j\left(\overline{k}\right)\left|\right]$ for any $\overline{k}$, where $n_i,n_j$ are obtained from independent simulations.

Figure 7

(a) Kinetic energy as a function of temperature at expansion orders $2\ge O\ge 4$ for the case $\mu =2t,U=\infty$ (the same as in Fig. 5). The two highest orders give results in good agreement down to $T=1/2$, i.e., slightly lower than the filling factor. At lower temperatures the correction from higher-order diagrams is substantial. (b) Normalizing by the carrier density we obtain the kinetic energy per hole, which can be compared to the corresponding figure for an ideal ferromagnet, shown by the dashed line [see also Eq. (34)]. The carrier density (number of holes per lattice site) is given by $1-\rho$ [see Eq. (31)].