Avoiding ergodicity problems in lattice discretizations of the Hubbard model

The Hubbard model arises naturally when electron-electron interactions are added to the tight-binding descriptions of many condensed matter systems. For instance, the two-dimensional Hubbard model on the honeycomb lattice is central to the ab initio description of the electronic structure of carbon nanomaterials, such as graphene. Such low-dimensional Hubbard models are advantageously studied with Markov chain Monte Carlo methods, such as hybrid Monte Carlo (HMC). HMC is the standard algorithm of the lattice gauge theory community, as it is well suited to theories of dynamical fermions. As HMC performs continuous, global updates of the lattice degrees of freedom, it provides superior scaling with system size relative to local updating methods. A potential drawback of HMC is its susceptibility to ergodicity problems due to so-called exceptional configurations, for which the fermion operator cannot be inverted. Recently, ergodicity problems were found in some formulations of HMC simulations of the Hubbard model. Here, we address this issue directly and clarify under what conditions ergodicity is maintained or violated in HMC simulations of the Hubbard model. We study different lattice formulations of the fermion operator and provide explicit, representative calculations for small systems, often comparing to exact results. We show that a fermion operator can be found which is both computationally convenient and free of ergodicity problems.

Figure 1

In the four panels, we show as lighter red points the spectrum of a fermion matrix with a particular $\alpha$ and discretization, as indicated, for a four-site honeycomb lattice with $N_t=96$. The darker, black points are the spectrum of the noninteracting case ($\varphi =0$), which depends only on which discretization is used. We show light gray rings to guide the eye for the eigenvalue rings as described by (33). There are four rings because there are four spatial sites; each ring has $N_t=96$ points. In all four cases, we use the same configuration $\varphi$ which was randomly sampled from a Gaussian with width $\sqrt{U\beta /N_t}$, where $U\beta =60$.

Figure 2

Eigenvalues $\lambda$ and $f$ as in (47) for $\varphi =\tau {\varphi }^{-}$, where $detM\left[{\varphi }^{-}\right]<0$ in all cases. The fictional time $\tau$ runs from 0 (purple), so that $\varphi$ is the noninteracting case and the determinant is positive, to 1 (yellow) so that the determinant is negative; intermediate values of $\tau$ are shown in the color bar on the right. The eigenvalue trajectories are described in the text. The eigenvalue locations were computed on an $N_t=20$ four-site honeycomb lattice with $\kappa \beta =4$. See Ref. [52] for animations of these eigenvalue dances.

Figure 3

Contours of constant probability, (69) (left, $\alpha =0$) and (71) (right, $\alpha =1$) with $U\beta =18$ and $N_t=2$ for the one-site problem, where the subscript on the Hubbard-Stratonovich fields $\varphi$ indicates the time slice. The dark contours represent high probability density, which decreases as the contours get lighter. The $\mathrm{\Phi }$ axis runs from bottom left to top right. (Right) The diagonal (black) lines correspond to locations where the probability weight (71) is exactly zero, creating barriers to HMC; their regular spacing is a consequence of the periodicity of the determinant, discussed in Sec. 2d3. We show every tenth configuration of 10 000 configurations (generated with a very precise integrator) as points. It is evident that the HMC algorithm was trapped in a region of high probability density in the $\alpha =0$ case by the wide separation with low probability density and in the $\alpha =1$ case by the lines of zero weight.

Figure 4

Histograms of $\mathrm{\Phi }$ and correlators for the one-site problem with very fine integration. Monte Carlo evolution was started in different lobes and the trajectories are stuck in those lobes, shown by the different colors. All ensembles consist of 50 k thermalized configurations. Single-particle correlators to the right are color coded to match the histograms in the left column. Points where $C\left(\tau \right)\le 0$ are not shown because of the logarithmic scale. The black lines show the exact results. The histograms are all normalized but, because HMC was trapped, the numerically generated histograms are much taller and thinner than the exact results.

Figure 5

Contours of constant probability weight (79) using the exponential discretization (left column) and the diagonal discretization (right column) with $\alpha =0$ (top row) and $\alpha =1$ (bottom row) for the two-site problem with $U\beta =18,\kappa \beta =1$, and $N_t=1$. Contour colors range from orange (light) for large to blue (dark) for small values in arbitrary units. In the bottom panels, the black diagonal squiggly lines and regularly spaced black dots correspond to locations where (79) is exactly zero; their regular spacing is a consequence of the periodicity of the determinant, discussed in Sec. 2d3. Crossing a black line changes the sign of the fermion determinant. The colored points show 10 k HMC trajectories. For $\alpha =0$, evolution was started in different modes of the field space, shown in different colors. For $\alpha =1$, evolution was started at the origin with a small random deviation. The outer histograms show the marginal distributions in ${\varphi }_1$ and ${\varphi }_2$ where the black lines are the exact results obtained from (79).

Figure 6

Analogous contours to Fig. 5 but using $\kappa =0.1$ to mimic approaching the continuum limit. The exponential and diagonal discretizations now have very similar probability contours. For the $\alpha =1$ case, the black lines of zero weight of the exponential discretization pinch and nearly close, while the zero-weight configurations remain isolated in the diagonal discretization.

Figure 7

HMC history of ${\mathrm{\Phi }}_x={\sum }_t{\varphi }_{xt}$ for the two-site problem with 100 k trajectories with $U\beta =60,\kappa \beta =6$, and $N_t=40$. Only every 100th trajectory of each ensemble is shown. Several HMC streams for each case were produced by starting in different modes in configuration space as indicated by color. Only for $\alpha =1$, diagonal there is no difference, where the run is started, so only one stream is shown. The black lines in the marginal distributions for ${\mathrm{\Phi }}_1$ and ${\mathrm{\Phi }}_2$ are the exact one-site distributions which are recovered in the strong-coupling limit. The contours show the product of the one-site distributions and are not exact results for this case.

Figure 8

Correlators from Eq. (80) color coded to match the ensembles shown in Fig. 7. Points where $C\left(\tau \right)\le 0$ are not shown because of the logarithmic scale. The black lines show exact results from direct diagonalization of the two-site Hubbard model. For clarity, the figure shows only $C_{+},C_{-}$ behaves in the same way.

Figure 9

Dependence of correlators on initial HMC configuration. The red correlators are calculated at $a\mathit{k}=(\pi /2,\pi /2)$ which for this case means that $C_{+}=C_{-}$. The green and blue correlators (lower, off center minima) are calculated at $a\mathit{k}=(0,0)$. For $\alpha =0$ the initial configuration was chosen such that $detM\left[\varphi \right]<0$. For $\alpha =1$, exponential discretization the starting configuration satisfies $0>e^{-i\mathrm{\Phi }/2}detM\left[\varphi \right]\in \mathbb{R}$. For $\alpha =1$, diagonal no particular starting criterion was set. All plots show data from 2400 trajectories. The starting criteria stayed fulfilled for all of those configurations.

Figure 10

Histograms of $detM$ for the $\alpha =0$ basis (top row) and diagonal $\alpha =1$ case (bottom right), and $e^{-i\mathrm{\Phi }/2}detM\left[\varphi \right]$ for the exponential $\alpha =1$ case (bottom left) for the ensembles used for Fig. 9. Except for the diagonal $\alpha =1$ case, all quantities are real. In the diagonal $\alpha =1$ case, $detM$ is complex and the domain of the histogram is the complex plane.

Figure 11

Same as Fig. 4 but with coarser integration with an acceptance rate of about $60\%$. For $\alpha =0$, Monte Carlo evolution was started in the right hand lobe and remained stuck there. Both ensembles consist of 50k thermalized configurations. The purple correlator is $C_{\mathbf{▾}}\left(\varphi \right)=(C_{•}\left(\varphi \right)+C_{•}(-\varphi ))/2$, the result of the symmetrization discussed in Sec. 4b.

Figure 12

HMC history of ${\mathrm{\Phi }}_x$ for the two-site problem with 100 k trajectories and a coarse integrator. Only every 100th trajectory of each ensemble is shown. The black lines in the marginal distributions for ${\mathrm{\Phi }}_1$ and ${\mathrm{\Phi }}_2$ are the exact one-site distributions which are recovered in the strong-coupling limit. The contours show the product of the one-site distributions and are not exact results for this case.

Figure 13

Correlators from Eq. (80) for the ensembles shown in Fig. 12. Points where $C\left(\tau \right)\le 0$ are not shown because of the logarithmic scale. The black lines show exact results from direct diagonalization of the two-site Hubbard model. For clarity, the figure shows only $C_{+},C_{-}$ behaves in the same way.

Figure 14

HMC history of ${\mathrm{\Phi }}_x={\sum }_t{\varphi }_{xt}$ for the two-site problem with 100 k trajectories and a sign flip every 100 trajectories according to (83). Only every 100th trajectory of each ensemble is shown. The black lines in the marginal distributions for ${\mathrm{\Phi }}_1$ and ${\mathrm{\Phi }}_2$ are the exact one-site distributions which are recovered in the strong-coupling limit. The contours show the product of the one-site distributions and are not exact results for this case.

Figure 15

Correlators from Eq. (80) from the ensembles with sign flips shown in Fig. 14. Points where $C\left(\tau \right)\le 0$ are not shown because of the logarithmic scale. The black lines show exact results from direct diagonalization of the one-site Hubbard model. For clarity, the figure shows only $C_{+},C_{-}$ behaves in the same way.

Figure 16

HMC trajectories and correlator for the exponential $\alpha =1$ case with the same parameters as Figs. 7 and 8, but large jumps as in Eqs. (88) and (89) every 100 trajectories. Only every 100 th trajectory is shown in the left panel and the ensemble consists of 100k trajectories in total.

Figure 17

Comparison of correlators on a $4\times 4$ square lattice. The data points for BSS are taken from Ref. [40]. The exponential $\alpha =0$ HMC run has negative $detM\left[\varphi \right]$ for all configurations. Started from a positive determinant this case produces results consistent with BSS and $\alpha =1$.

Figure 18

A variety of correlation functions (shown in different colors) of local bilinear operators as a function of Euclidean time, for the exponential and diagonal discretizations, projected to the uniform and staggered irreps on both spatial indices. Exact results are shown as thin solid lines, while numerical results are shown as markerless error bars, and the respective equal-time order parameters are indicated in text, though the volume factors are omitted for clarity. In the left panel, the dark and light results correspond to different ensembles, as described in the text. For additional details and discussion, see the text.