Benchmark study of the two-dimensional Hubbard model with auxiliary-field quantum Monte Carlo method

Ground-state properties of the Hubbard model on a two-dimensional square lattice are studied by the auxiliary-field quantum Monte Carlo method. Accurate results for energy, double occupancy, effective hopping, magnetization, and momentum distribution are calculated for interaction strengths of $U/t$ from 2 to 8, for a range of densities including half-filling and $n=0.3,0.5,0.6,0.75$, and $0.875$. At half-filling, the results are numerically exact. Away from half-filling, the constrained path Monte Carlo method is employed to control the sign problem. Our results are obtained with several advances in the computational algorithm, which are described in detail. We discuss the advantages of generalized Hartree-Fock trial wave functions and its connection to pairing wave functions, as well as the interplay with different forms of Hubbard-Stratonovich decompositions. We study the use of different twist angle sets when applying the twist averaged boundary conditions. We propose the use of quasirandom sequences, which improves the convergence to the thermodynamic limit over pseudorandom and other sequences. With it and a careful finite size scaling analysis, we are able to obtain accurate values of ground-state properties in the thermodynamic limit. Detailed results for finite-sized systems up to $16\times 16$ are also provided for benchmark purposes.

Figure 1

The convergence rate of the ground-state energy of the noninteracting Hubbard model computed using TABC with twist angles from QR sequence, PR sequence, and uniform grid. The system is $4\times 4$ at half-filling. The vertical axis shows the absolute value of the relative error with respect to the exact value, which is $-16/{\pi }^2$.

Figure 2

(a) The error bar of the ground-state energy computed from TABC vs the number of twist angles used with twists generated by QR sequence, PR sequence, and the uniform grid. The system is $4\times 4$ with 2 up and 2 down electrons, and $U=8$. The ground-state energies are obtained by the ED method. Note the $log\text{−}log$ scale of the plot. (b) Similar results for an $8\times 8$ lattice at $n=0.5$, with $U=4$, for QR and PR twist angles. The energies are computed by the CPMC method.

Figure 3

Ground-state energy at half-filling calculated using different boundary conditions. PBC, PBC-APBC and TABC data are represented by black triangular, blue star, and red dot, respectively. A fit of the TABC data is also shown, with solid red line. (a), (b), (c), and (d) correspond to results for $U=2,4,6,8$. In the insets of each panel, a zoom of the TABC results and the fit are shown for large supercell sizes. The cyan dot in each inset represents the TDL value and combined statistical and twist error bars and the uncertainty from the fit.

Figure 4

Double occupancy at half-filling calculated using different boundary conditions. Symbols and setup are similar to Fig. 3.

Figure 5

The dependence of the effective hopping, $t_{\mathrm{eff}}/t$, on the interaction strength $U$ at half-filling. The inset shows the corresponding potential energy in units of the noninteracting kinetic energy.

Figure 6

Spin correlation function in the ground state at half-filling. System sizes ranging from $4\times 4$ to $16\times 16$ are shown, under PBC, with $U=4$. The horizontal axis is the relative distance, $\sqrt{x^2+y^2}$. The top panel shows the spin correlation function, while the bottom panel shows the staggered correlation. The dashed horizontal line in (b) shows the final TDL value obtained from the fit.

Figure 7

Magnetization computed with TABC at half-filling for (a) $U=4$ and (b) $U=8$. For each choice of $d$, the result of a fit using the form in Eq. (21) is also plotted. The cyan dot represents the final TDL value and the estimated error bar.

Figure 8

Ground-state energy and double occupancy vs supercell size at $n=0.5$ for $U=4$ and 8. TABC is used. (a) and (b) correspond to $U=4$. (c) and (d) correspond to $U=8$. The solid lines are from a fit using $E_0/L^2\left(D\right)\sim e_0\left(D_0\right)+a/L^3$.

Figure 9

Momentum distribution at $n=0.5$ for (a) $U=4$ and (b) 8. The horizontal axis is the noninteracting energy for the given momentum normalized by the noninteracting Fermi energy of the corresponding twist.