Improved Multi-Variable Variational Monte Carlo Method Examined by High-Precision Calculations of One-Dimensional Hubbard Model

We revisit the accuracy of the variational Monte Carlo (VMC) method by taking an example of ground state properties for the one-dimensional Hubbard model. We start from the variational wave functions with the Gutzwiller and long-range Jastrow factor introduced by Capello et al. [Phys. Rev. B 72, 085121 (2005)] and further improve it by considering several quantum-number projections and a generalized one-body wave function. We find that the quantum spin projection and total momentum projection greatly improve the accuracy of the ground state energy within 0.5% error, for both small and large systems at half filling. Besides, the momentum distribution function n(k) at quarter filling calculated up to 196 sites allows us direct estimate of the critical exponents of the charge correlations from the power-law behavior of n(k) near the Fermi wave vector. Estimated critical exponents well reproduce those predicted by the Tomonaga-Luttinger theory.


Introduction
Strongly correlated electron systems have attracted much interest over the past decades. Although rigorous ground states are obtained analytically for some simple models, numerical methods offer powerful tools for understanding more realistic systems. There are several numerical methods for the purpose of calculating low-energy properties of the strongly correlated systems. At finite temperatures, the auxiliary-field quantum Monte Carlo (AFQMC) method [1,2,3] gives exact ground states within the statistical errors. However, for the geometrically frustrated systems, the errors become exponentially large, known as the negative sign problem. To overcome this difficulty, the Gaussian-basis Monte Carlo (GBMC) method [4,5,6] has been proposed as an alternative to the AFQMC method, although the tractable range of parameters are still under debate. The dynamical mean-field theory (DMFT) [7,8] can also be applied to strongly correlated electron systems. However, when we wish to take into account sufficient spatial correlations to represent fluctuations around broken symmetry states, we need larger clusters beyond the present computational accessibility. At zero temperature, the exact diagonalization (ED) method gives exact ground states. This method, however, is not able to handle large system sizes. The density matrix renormalization group (DMRG) method [9] gives nearly exact ground states for relatively larger systems, though it is mainly applied to onedimensional systems, or two-dimensional systems with cylindrical boundary conditions. The unbiased path-integral renormalization group (PIRG) method [10,11,12] has been applied to several strongly correlated electron systems, with essentially accurate ground states, even in the presence of geometrical frustrations. However, strong correlation regime requires a set of many basis functions for accurate estimates, sometimes beyond the computational accessibility. Among them, although a bias inevitably exists, the VMC method [13] is applicable to large systems in the presence of strong correlations and geometrical frustrations.
Accuracy of the calculated properties obtained by the VMC method relies on the choice of the trial wave functions. In most cases, the trial wave functions |ψ are given by Slater determinants |D supplemented by correlation factors P as |ψ = P |D . One of the famous example of this trial function is a "projected Bardeen-Cooper-Schrieffer (BCS)" type wave function |ψ = P ∞ G |BCS , where P ∞ G = i (1 − n i↑ n i↓ ) denotes a Gutzwiller factor, which prohibits double occupancy of the sites, and |BCS = k (u k + v k c † k↑ c † −k↓ ) |0 denotes a BCS state, which is the ground state of the BCS mean-field Hamiltonian. This projected BCS wave function is equivalent to the so-called "resonating valence bond (RVB)" wave function, which is the superposition of all dimer coverings [14,15]. The projected BCS and RVB wave functions provide highly accurate descriptions in two-dimensional strongly correlated electron systems [16,17].
In principle, the trial wave functions based on Slater determinants become accurate for higher dimensions since the mean-field-like treatment becomes valid above the upper critical dimension. Although it seems difficult to capture the low-energy properties at lower spatial dimensions by using such a trial wave function, long-range correlation factors may restore the accuracy of the wave functions [18,19]. Capello et al. have calculated the ground state of the one-dimensional Hubbard model by using the projected BCS type wave functions. They have shown that the long-range Jastrow factor is essential to describe the low-energy properties. They have also reproduced correct power-law behaviors of the spin and charge correlation functions.
Although Capello et al.'s scheme appears to substantially reproduce the low-energy properties of the one-dimensional model, the errors in the energy become larger even for small system sizes, typically more than 2% for L = 18 for intermediate interaction strength [19]. This is mainly because they have used the Slater determinants with a limited number of variational parameters. The Slater determinants in their wave function are just a Fermi sea of the tight binding model, or short-range BCS wave functions with nonzero amplitudes at most up to third-nearest-neighbors. In addition, in practical numerical calculations, since the accessible system sizes are limited, it is important to prepare the wave functions preserve the correct quantum numbers within its accessible range of sizes.
In this paper, to obtain more accurate wave functions of the low-dimensional systems for the smaller as well as for larger system sizes, we revisit the variational calculations of the onedimensional Hubbard model, as a reference system. By considering the multi-variable Slater determinants and projections which restore the SU(2) spin rotational and spatially translational symmetries, we obtain more accurate ground states of the model than the previous studies. We also calculate the momentum distribution function n(k) up to 196 sites, the first high-precision variational calculation of n(k) to our knowledge, and directly estimate the critical exponents of the charge correlations. We compare them with the analysis by Kawakami et al. [20,21] as well as Ogata et al. [22]. Estimated critical exponents well reproduce those obtained by the Tomonaga-Luttinger theory.

Model and method
We consider the one-dimensional Hubbard model. The Hamiltonian is given by where c † iσ (c iσ ), n iσ , t and U denote the creation (annihilation) operator, the number operators, the nearest-neighbor hopping interaction, and the on-site Coulomb interaction, respectively.
Exact ground states of the one-dimensional Hubbard model are obtained by Lieb and Wu [23]. At half filling, the ground state is a metal for U = 0, and is an insulator for nonzero positive U . Away from half filling, the ground state is a metal for any U ≥ 0. Unlike the normal Fermi liquid, the momentum distribution of this metal does not show a finite jump at the Fermi wave vector. Instead, it shows power-law decays near the Fermi wave vector. This type of the ground state, which appears in one-dimensional electron systems, is well known as the Tomonaga-Luttinger liquid [24,25].
To obtain accurate variational description of the ground state of the model, we use a multivariable variational Monte Carlo (mVMC) method [17]. We use the trial wave functions given by where |φ pair , P G , P J , P dh , L S=0 , and L K=0 denote the one-body part, the Gutzwiller factor, the Jastrow factor, the doublon-holon factor, the spin quantum-number projection on to S = 0 space, and the total momentum projection on to K = 0 space, respectively, as we detail later. We calculate the ground states of the model under the periodic boundary conditions at half filling (n = 1) for N s = 4n + 2 sites, and at quarter filling (n = 0.5) for N s = 8n + 4 sites, so that the systems obey the closed shell conditions. For the one-body part |φ pair , we use a generalized pairing wave function defined as where f ij , N s , and N e denote the variational parameters, the number of sites, and the number of electrons, respectively. When all f ij 's are optimized simultaneously and independently, since the total number of f ij is N 2 s , the computation time for the optimization is highly demanding. To reduce the computation time to O(N s ), we assume f ij to have a sublattice structure, namely where σ(i) denotes a sublattice index at site i. Since optimized ground states by using f ij with reflectional symmetries, namely f ij = f ji , have higher energies, we do not impose reflectional symmetries on f ij . Before calculating the ground states of larger system sizes, we examine the accuracy of the wave functions using several sublattice structures of f ij for smaller system sizes less than about 20 sites. We find that optimized energies for two-sublattice structure and N s -sublattice (fully optimized) structure are nearly the same. This suggests that the ground states are well-described by the wave functions with the two-sublattice structure for the one-dimensional Hubbard model. In the following, we show results for mainly two-sublattice structure, and one-sublattice translationally invariant structure.
To include correlation effects beyond the one-body approximation, we consider the correlation factors and where g, v ij , and α ml denote variational parameters of the Gutzwiller, Jastrow, and the doublonholon correlation factors, respectively. We take into account long-range part of the Jastrow factor, which is indispensable for reproducing the low-energy properties of the one-dimensional Hubbard model [18,19], and the doublon-holon factor. A variable ξ iml in Eq. (6) is a many-body operator, which gives one if a doublon (holon) exists at the i-th site and m holons (doublons) exist at l-th nearest neighbors, and gives zero otherwise. We assume the Jastrow factor to have translational symmetry, namely v ij = v(r ij ), where r ij = |R i − R j | is the distance between the sites i and j. We note that the Gutzwiller factor and the doublon-holon factor have the translational symmetry by definition.
In addition to the one-body part and the correlation factors, we introduce a projection L S=0 , which restores the SU(2) spin rotational symmetry, and a projection L K=0 , which restores the translational symmetry. The spin quantum-number and total momentum projections are switched on at all optimization steps. Since the correlation factors have the SU(2) and translational symmetries, the quantum-number projections commute with the correlation factors.
To optimize a large number of these variational parameters, we use the stochastic reconfiguration (SR) method developed by Sorella [26]. We first optimize variational parameters f ij for typically 1000 SR steps with fixed correlation factors. After confirming the energy convergence, we simultaneously optimize all the variational parameters (g, v ij , α ml , and f ij ) for more than 4000 SR steps. The black crosses denote the exact energies obtained by the ED method for finite sizes and analytically at the thermodynamic limit. The purple down-pointing triangles and red squares denote calculated energies by the mVMC method with the translationally invariant and two-sublattice f ij trial wave functions, respectively. The former is obtained without the quantum-number projections and the latter is with the spin and momentum projections to the total singlet and zero total momentum. Two and 0.5 percent errors measured from the exact results are plotted by dash-dotted purple and dashed red curves, respectively. (b) Coulomb interaction U/t dependence of errors in calculated energies at half filling for N s = 18. The red squares denote our mVMC results obtained by using the two-sublattice f ij trial wave functions. The errors in the energy are much smaller than their results. The green triangles and blue circles denote the errors in the energy obtained by Capello et al. by using short-range BCS wave functions with the third-nearest-neighbor amplitudes [19]. The former is obtained with the assumed Jastrow factor and the latter is with the optimized Jastrow factor.

Results
We obtain the ground-state energy of the one-dimensional Hubbard model at half filling, as shown in Fig. 1. For the case of the two-sublattice f ij , the errors in the energy are within approximately 0.5%, much smaller than the previous calculations by Capello et al. [19] The extended one-body part and quantum projections improve the accuracy of the ground state energy, even for the small system sizes. We show results of the two-sublattice f ij condition, hereafter. In the gapless one-dimensional systems, spin and charge correlations show power-law decays [27,28]. The asymptotic expressions of the spin correlation functions are given as where A, k F , and K ρ denote a constant, the Fermi wave vector, and a critical exponent, respectively. The Fermi wave vector k F is written in terms of the filling n through the equation k F = nπ/2. The critical exponent K ρ is a function of the Coulomb interaction U and the filling n. At half filling, the charge degrees of freedom are gapped out with the exponential decay of correlations in distance, leaving only the gapless spin correlation in the form of Eq. (7) formally with K ρ = 0 for nonzero repulsive U . Namely, the spin correlation functions decay as C(r) ∼ A(−1) r √ ln r/r [27,28,29]. Away from half filling, the critical exponent satisfies 0 < K ρ < 1 for nonzero repulsive U . Thus the spin correlation functions decay as Eq. (7). To include finite size effects under the periodic boundary conditions, we compare the calculated correlation functions with C(r) + C(N s − r) instead of the original correlation function C(r) because of the reflection effects from the boundary. As shown in Fig. 2, calculated spin correlation functions of the Hubbard model reproduce √ ln r/r decay at half filling, and √ ln r/r 1+Kρ and 1/r 2 decays at quarter filling with the critical exponent K ρ ≈ 0.7 at U/t = 4, respectively. This critical exponent is nearly the same as the exact value, K ρ = 0.711, predicted by the Tomonaga-Luttinger theory [18,28,20,21].
In our calculations, the one-body part of the trial wave function includes long-range part in f ij and v ij . Although there should be no long-range order for the one-dimensional Hubbard model, our trial wave function could fallaciously stabilize the insulator with the antiferromagnetic long-range order. To examine whether our calculations correctly reproduce the absence of the Figure 3. Peak values of spin structure factors of Hubbard model at U/t = 4 and U/t = 10 for each system size, (a) at half filling and (b) at quarter filling. Statistical errors are smaller than the symbol sizes. At half filling, the spin structure factor shows a peak at q peak = π, while at quarter filling, it shows a peak at q peak = π/2. The filled symbols denote results obtained by the ED method for the model under the periodic boundary conditions, while the open symbols denote the mVMC results. The red squares denote results at U/t = 4, and the green circles denote results at U/t = 10. The squared magnetizations, S(q peak )/N s converge to zero in the thermodynamic limit.
long-range order in the model, we calculate the peak values of the spin structure factors for each system size. For the insulator, size dependencies of the peak values of the spin structure factor are expected to follow the scaling S(q peak ) = 1 3 dr e iq peak r S(0) · S(r) ∼ where Λ denotes the cutoff, while for the metal, since 1/r 1+Kρ ≫ 1/r 2 for r ≫ 1, we obtain By using these equations, we extrapolate the squared magnetization, the peak values of the spin structure factor S(q peak ) divided by the system size N s , in the thermodynamic limit. As shown in Fig. 3, the squared magnetizations, S(q peak )/N s converge to zero in the thermodynamic limit. These results suggest that our trial wave functions with the long-range f ij and symmetryrestoring projections L are able to reproduce ground states with no long-range magnetic order with power-law decay of correlations. For the one-dimensional Hubbard model, the momentum distribution shows a power-law singularity with a critical exponent α near the Fermi wave vector k F , contrary to the ordinary Fermi liquid, which shows a jump at the Fermi wave vector [20,21]. The critical exponent α is written in terms of K ρ through the equation α = (K ρ + 1/K ρ − 2)/4 [30,31,28]. For the infinite Coulomb interactions, Ogata et al. have calculated the momentum distribution n(k) of the metal ground states by using the Bethe-ansatz wave functions, and have estimated the critical exponents α from Eq. (12) [22]. We calculate the momentum distribution of the Hubbard model at quarter filling up to 196 sites. As shown in Fig. 4(a), calculated momentum distributions show powerlaw behavior near the wave vector k F . Another singularity also appears at k = 3k F , as expected in the previous studies [22]. Figure 5. Estimated critical exponent α from relation |n(k) − n(k F )| ∝ |k − k F | α near k F . Statistical errors are smaller than the symbol sizes. The red squares denote results at U/t = 4, and the green circles denote results at U/t = 10. These values are nearly the same as the exact results, α = 0.029 (U/t = 4) and α = 0.069 (U/t = 10), predicted by the Tomonaga-Luttinger theory [18,28,20,21].
To estimate the critical exponent α, we plot ln |n(k) − n(k F )| as a function of |k − k F | for each system size, as shown in Fig. 4(b), and extrapolate the slopes from data points by using Eq. (12). Since large finite-size effects appear near |k − k F |/π 0.1, we use data that satisfy |k − k F |/π ∈ [0.1, 0.4] for fitting. As shown in Fig. 5, as the system size N s goes to infinity, the estimated critical exponent α(N s ) gradually converges to the exact value, α = 0.029 at U/t = 4 and α = 0.069 at U/t = 10, predicted by the Tomonaga-Luttinger theory [18,28,20,21].

Summary
To improve the VMC method, we revisit the variational description of the ground states of the one-dimensional Hubbard model, by considering the multi-variable Slater determinants as well as the long-range correlation factors. We find that the quantum spin projection and total momentum projection greatly improve the accuracy of the ground state energy with less than 0.5% error, irrespective of the system size. By calculating the momentum distribution function n(k) of the model at quarter filling up to 196 sites, we have shown that our method is able to reproduce power-law decays at the momenta k F and 3k F . By directly estimating the critical exponents of the charge correlations from the power-law behavior of n(k) near k F , we have also shown that estimated critical exponents well reproduce those predicted by the Tomonaga-Luttinger theory.