{'fraction_non_alphanumeric': 0.046991672856640215, 'fraction_numerical': 0.01954728488212987, 'mean_word_length': 4.242877638860422, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Benchmark study of an auxiliary-field quantum Monte Carlo technique for the Hubbard model with shifted-discrete Hubbard-Stratonovich transformations

Within the ground-state auxiliary-field quantum Monte Carlo technique, we introduce discrete Hubbard-Stratonovich transformations (HSTs) that are suitable also for spatially inhomogeneous trial functions. The discrete auxiliary fields introduced here are coupled to local spin or charge operators fluctuating around their Hartree-Fock values. The formalism can be considered as a generalization of the discrete HSTs by Hirsch [J. E. Hirsch, Phys. Rev. B 28, 4059 (1983)] or a compactification of the shifted-contour auxiliary-field Monte Carlo formalism by Rom et al. [N. Rom et al., Chem. Phys. Lett. 270, 382 (1997)]. An improvement of the acceptance ratio is found for a real auxiliary field, while an improvement of the average sign is found for a pure-imaginary auxiliary field. Efficiencies of the different HSTs are tested in the single-band Hubbard model at and away from half filling by studying the staggered magnetization and energy expectation values, respectively.


I. INTRODUCTION
The numerical solution of the Hubbard model with strong correlations is one of the most challenging issues in the theory of strongly correlated electron systems [1,2]. Attempts to determine the ground state are often based on iterative techniques based on a repeated application of a short imaginary time propagator, or by using the simple power method and more advanced Krylov-subspace techniques, as, for instance the Lanczos algorithm, where a Hamiltonian operator is repeatedly applied to a properly chosen trial state. In both cases the ground-state component of the trial state is filtered out after several iterations.
Among these projection techniques, the auxiliary-field quantum Monte Carlo (AFQMC) [3][4][5][6] is one of the most powerful schemes, as it allows us to study, for example, the ground-state properties of the Hubbard model with several thousands electrons and lattice sites, when the negative-sign problem is absent [7][8][9]. In the ground-state AFQMC, even if the Hamiltonian is the same, there exists some arbitrariness in choosing the trial wave function and the type of auxiliary fields (e.g., real, complex, continuous, or discrete). Experience has shown that an appropriate choice of these ingredients may significantly improve the efficiency of the Monte Carlo simulations [10].
It has been demonstrated [2,11] that a Slater-determinant obtained from an unrestricted Hartree-Fock (UHF) approximation [12] provides a good trial wave function for the doped Hubbard model in the constrained-path AFQMC [13]. Recently, for a particular parameter set at doping δ = 1/8 and electron-electron repulsion U/t = 8, the ground state of the Hubbard model on the square lattice has been predicted to exhibit a vertical stripe order [2], where the stripe states with periods λ = 5, 6, 7 and 8 in units of the lattice constant are nearly degenerate, while a spatially homogeneous d-wave superconducting state should have, according to their study, a higher energy. Recent variational Monte Carlo (VMC) calculations [14][15][16] have also shown that various vertical-stripe orders with different periods appear depending on the doping and the hopping parameter. In most of the calculations in Ref. [2], the symmetry of finite-size clusters is broken due to the use of UHF trial wave functions or by applying pinning magnetic fields, and the results are extrapolated to the thermodynamic limit. The success of utilizing symmetry-broken wave functions is rather surprising, because symmetry breakings do not occur in the exact ground state of finite-size systems. The similar issue is known as the symmetry dilemma in first-principles calculations for molecules [17,18]. Recently, it has been shown that the quality of the trial wave function can be improved by restoring the symmetries that are once broken by UHF or mean-field treatments [19][20][21]. However, in the present work, we do not enter into the issue on symmetry breakings of trial wave functions, and rather focus on the arbitrariness of the auxiliary field to improve the efficiency of AFQMC simulations with such symmetry-broken trial wave functions.
The way of transforming a quartic interaction term into a quadratic one via the Hubbard-Stratonovich transformation (HST) [22] is not unique and affects the efficiency of simulations [23][24][25][26][27][28]. Recently the popularity of this technique is substantially increased, because it has been realized that, with continuous auxiliary fields, one can treat interaction terms beyond the on-site Hubbard interaction, up to the complete treatment of the long-range Coulomb interaction [29][30][31][32][33], or of the long-range electron-phonon interaction [34], and even both of them on the same footing [35], without being vexed by the sign problem in a certain parameter region on bipartitle lattices. Interestingly, such a parameter region coincides with the one where rigorous statements on the ground state of an extended Hubbard-Holstein model are available [36,37]. It is also noteworthy that, even when the sign problem cannot be eliminated completely, continuous auxiliary fields with a proper shift [38][39][40] can improve the efficiency of simulations compared to the one without the shift. A similar idea has been employed also in the AFQMC [41,42] within the constrainedpath approximation [13].
In this paper, we introduce shifted-discrete HSTs, where auxiliary fields are coupled to the fluctuation of local spin or charge. The method is applied for AFQMC simulations of the Hubbard model on the square lattice. It is shown that the shifted-discrete HSTs can improve the efficiency of the AFQMC simulations. Moreover we present results on the magnetic order parameter as a function of U/t with high statistical accuracy, that represents important benchmark, useful also for comparison with experiments. The rest of the paper is organized as follows. In Sec. II, the Hubbard model is defined and the AFQMC method is described. In Sec. III, the shifted-discrete HSTs are introduced. In Sec. IV, numerical results of the AFQMC simulations for the Hubbard model are presented. Section V is devoted to conclusions and discussions.

II. MODEL AND METHOD
We consider the Hubbard model whose Hamiltonian is defined byĤ =K +V, wherê c † iσ (ĉ iσ ) creates (annihilates) a fermion with site index i and spin index σ(=↑, ↓),n iσ =ĉ † iσĉ iσ , t is the hopping parameter between the nearest-neighbor sites on the square lattice and U > 0 is the on-site electron-electron repulsion. We consider the Hubbard model on N = L × L-site clusters. Boundary conditions will be specified for each calculation in Sec. IV. The lattice constant is set to be unity.
In the AFQMC, the expectation value of an operatorÔ is calculated as where β is the projection time and |Ψ T is a trial wave function. If β is infinitely large, one can obtain the ground-state expectation value as long as |Ψ T has a finite overlap with the ground state [43]. If β is finite, the results depend on the trial wave function (see for example Ref. [44]). If β = 0, Eq. (3) reduces to the expectation value ofÔ with respect to the trial wave function. At finite dopings, |Ψ T is obtained by solving the eigenvalue problem of the following UHF Hamiltonian self consistently: where U eff is an arbitrary parameter and the expectation value · · · 0 in Eq. (4) is defined in Eq. (3) with β = 0. A fine tuning of U eff can improve the quality of the trial wave function [11]. We set U eff /t = 2.5 which has turned out to provide a good trial wave function for the doped cases studied here, in the sense that the energy expectation value decreases quickly with increasing β. By adding a small bias in the initial condition for the self consistent UHF loop to pin the direction of the stripe, |Ψ T shows a vertical stripe order with period λ = 8 around δ = 1/8 doping on the 16 × 16 cluster. At half filling, |Ψ T is obtained as a ground state of noninteracting electrons on the square lattice under a staggered magnetic field along the spin-quantized axis (z direction): where (−1) i = 1(−1) if the site i belongs to the A (B) sublattice and ∆ AF can be chosen arbitrarily. The value of ∆ AF will be specified with the numerical results in Sec. IV. By using the second-order Suzuki-Trotter decomposition [45,46], the imaginary-time propagator can be expressed as where the projection time β is discretized into N τ time slices and ∆ τ = β/N τ . For the doped cases, we set ∆ τ t = 0.05 so that the discretization error is within statistical errors. For the half filled case, we perform extrapolations of ∆ τ → 0 to eliminate the discretization error, which becomes non negligible for large U/t as compared to statistical or extrapolation errors for the results shown in Sec. IV B. An HST is applied to e −∆ τV and the summation over the auxiliary fields is performed by the Monte Carlo method with the importance sampling, where a proposed auxiliary-field configuration is accepted or rejected according to the Metropolis algorithm. In the next Section, we introduce shifted-discrete HSTs for e −∆ τV .

III. SHIFTED-DISCRETE HUBBARD-STRATONOVICH TRANSFORMATIONS
In this Section we derive shifted-discrete HSTs which couple the auxiliary field to the local spin fluctuation in Sec. III A and to the local charge fluctuation in Sec. III B. Although the two HSTs can be formulated almost in parallel, we provide both of them separately for completeness.

A. Auxiliary field coupled to spin fluctuation
The Hubbard interaction in Eq. (2) can be written aŝ wherem i is an arbitrary number. Then e −∆ τ V can be written as Let us consider the first exponential factor in the right-hand side of Eq. (8). For each site i, we consider the following HST: where s i = ±1 is the discrete auxiliary field, and the undetermined four parameters α i , m i ,m i and C i are related through the following three equations (see Appendix A for derivation): Therefore, if say m i is given, α i ,m i , and C i are determined from Eqs. (10)- (12). Finally we obtain Note that in general m i m i and C i 's are irrelevant for results of simulations because they cancel out from the numerator and the denominator in Eq. (3). If m i = 0, the HST reduces to the one introduced by Hirsch [23]. However, the arbitrariness of m i can be utilized to improve the efficiency of AFQMC simulations as shown in Sec. IV.
In the right-hand side of Eq. (13), the auxiliary field α i s i is shifted by ∆ τ Um i as compared to the case of m i =m i = 0. To obtain more physical intuitions for m i , we rewrite the exponent of the right-hand side of Eq. (13) as In the first term, the auxiliary field α i s i is coupled to the fluctuation of the local magnetization (n i↑ −n i↓ − m i ), while the shift of the local magnetization by −m i in the first term is compensated by the spatially inhomogeneous magnetic field ∆ τ Um i in the second term.
We set the parameter m i as the local magnetization in the trial wave function This m i can be easily calculated and is expected to stabilize the simulation by keeping the first term in Eq. (14) "small" during the imaginary-time evolution. For a given m i , α i can be determined from Eq (10),m i from Eq. (11), and C i from Eq. (12). The solution α i of Eq. (10) can be found by the Newton method with an initial guess α i,initial = cosh −1 e ∆ τ U/2 , for example.

B. Auxiliary field coupled to charge fluctuation
In this subsection, α i and C i will be re-defined. The Hubbard interaction in Eq. (2) can be written aŝ whereñ i is an arbitrary number. Then e −∆ τ V can be written as Let us consider the first exponential factor in the right-hand side of Eq. (17). For each site i, we consider the following HST: where s i = ±1 is the discrete auxiliary field, and the undetermined four parameters α i , n i ,ñ i and C i are related through the following three equations (see Appendix A for derivation): Therefore, if say n i is given, α i ,ñ i , and C i are determined from Eqs. (19)- (21). Finally we obtain Note that in general n i ñ i and C i 's are irrelevant for results of simulations because they cancel out between the numerator and the denominator in Eq. (3). If n i = 1, the HST reduces to the one introduced by Hirsch [23]. However, the arbitrariness of n i can be utilized to improve the efficiency of AFQMC simulations as shown in Sec. IV.
In the right-hand side of Eq. (22), the auxiliary field iα i s i is shifted by ∆ τ U(1 −ñ i ) as compared to the case of n i =ñ i = 1. To obtain more physical intuitions for n i , we rewrite the exponent of the right-hand side of Eq. (22) as In the first term, the auxiliary field iα i s i is coupled to the fluctuation of the local density (n i↑ +n i↓ −n i ), while the shift of the local density by −(1 − n i ) in the first term is compensated by the spatially inhomogeneous chemical potential ∆ τ U(1 −ñ i ) in the second term.
We set the parameter n i as the local charge density in the trial wave function This n i can be easily calculated and is expected to stabilize the simulation by keeping the first term in Eq. (23) "small" during the imaginary-time evolution. For a given n i , α i can be determined from Eq (19),ñ i from Eq. (20), and C i from Eq. (21). The solution α i of Eq. (19) can be found by the Newton method with an initial guess α i,initial = cos −1 e −∆ τ U/2 , for example.

A. Finite dopings
At finite dopings, the sign problem occurs [47,48]. In the presence of the sign problem, the projection time β cannot be taken as large as that for the half filled case because the average sign (of the statistical weight) decreases exponentially in β [48], otherwise the number of statistical samplings has be increased exponentially to keep the statistical error small. We set the maximum β at which the average sign is ∼ 0.05. It will be shown that even in the presence of the sign problem, the AFQMC can still provide a good upper bound of the groundstate energy. Figure 1 shows the energy per site E(β) = Ĥ β /N, the average sign, and the acceptance rate as a function of β at U/t = 8 for the 16 × 16 cluster with 224 electrons, corresponding to δ = 1/8 = 0.125. Note that since E(β) is a decreasing function of β and its slope dE(β)/dβ is proportional to the energy variance [43]. The energies calculated by different HSTs coincide within the statistical errors, though their reachable β is different, as shown in Fig. 1(a). In Fig. 1(b), the acceptance rate of the real auxiliary field with the shift (HST spin with shift) is increased the one from without the shift (HST spin without shift). The reason can be attributed to that since the first term of Eq. (14) with a relevant m i is expected to be "smaller" than that with m i = 0, the factor e ±2α i (n i↑ −n i↓ −m i ) is closer to unity, so that the fluctuation of the norm of the determinant ratio is stabilized. On the other hand, the shift of the real auxiliary field does not affect the average sign significantly because the shift does not affect the sign of the determinant ratio, as can be seen in Fig. 1(c). The situation is different for the pure-imaginary auxiliary fields. Without the shift (HST charge without shift), the average sign diminishes significantly, even at βt = 0.1. By introducing the shift (HST charge with shift), the average sign is improved significantly. The reason can be attributed to that since the first term of Eq. (23) with a relevant n i is expected to be "smaller" than that with n i = 1, the factor e ±2iα i (n i↑ +n i↓ −n i ) is closer to unity, so that the fluctuation of the phase of the determinant ratio is stabilized. However, the average sign is still quite smaller than that with the real auxiliary fields. Although the acceptance rate is higher than the real auxiliary fields, the pure-imaginary fields may not be practical in the presence of the sign problem.
To show the usefulness of the AFQMC with a short imaginary-time propagation, we make a comparison with the state-of-the-art variational wave functions for the Hubbard model [14,15]. To this end, we move to the smaller doping with the larger U/t, where the more severe sign problem is expected. Figure 2 shows the energy per and the average sign as a function of β at U/t = 10 for the 16 × 16 cluster with 228 electrons, corresponding to δ = 0.109375. Here, only the shifted real auxiliary field is employed because it turned out to be the most efficient, as shown in Fig. 1 for U/t = 8 and δ = 0.125. We use periodic-(antiperiodic-) boundary condition in the x (y) direction to compare directly with the reference VMC results [14,15]. Notice that our AFQMC energy, computed at finite projection time β when the average sign is sufficiently large, respects the Ritz's variational principle [see Eq. (3)], because it corresponds to the variational expectation value ofĤ over the state e − β 2Ĥ |Ψ T / Ψ T |e −βĤ |Ψ T 1 2 . This is a useful property of an approximate technique that is not always satisfied, as for instance for the constrained-path AFQMC [13,41]. At βt = 0.7, where the average sign remains ∼ 0.99, the AFQMC energy is already lower than the VMC energy without variance extrapolation. At βt = 1.1, the AFQMC energy almost coincides with the VMC variance-extrapolated one, while the slope dE(β)/dβ is still finite, indicating that the AFQMC energy variance is nonzero [see Eq. (25)]. At βt = 1.5, the AFQMC energy is E/t = −0.6552 (4), which is lower than  [14,15]. FTTN stands for fat-tree tensor network and Var. ext. for variance extrapolation. the variance-extrapolated VMC energy E/t = 0.6538(9) [14] which may be compatible with our number within two standard deviations. This result suggest that the ground-state AFQMC method remains very useful for providing upper bound values of the ground-state energy even in the presence of the negative-sign problem.

B. Half filling
At half filling, the sign problem is absent. Therefore the AFQMC can provide exact results which often serve as a reference benchmark for other numerical techniques. An excellent agreement in the ground-state energies of the twodimensional Hubbard model between the AFQMC and other many-body techniques has been reported in Ref. [1]. Moreover, within the AFQMC, the staggered magnetization m, i.e., the order parameter at half filling, can be estimated accurately by using the twist-averaged boundary condition for small U/t, e.g., U 4 [49,50]. However, for large U/t, AFQMC simulations still face a difficulty of large fluctuations of the magnetization, which often lead to a relatively large error bar in m [1,49]. The same difficulty arises also in finite-temperature determinant QMC simulations [51,52]. In previous works, in order to overcome the difficulty, a pinning-field method has been proposed with a clear improvement for the determination of m in the thermodynamic limit [53,54]. In the following, we report an accurate estimate m especially for large U/t by making use of a symmetry-broken trial wave function. Figure 3 shows the staggered magnetization along the z di- rection as a function of the Monte Carlo sweep with different HSTs. The calculations are done for U/t = 8, βt = 24, and ∆ τ t = 0.1 on the L = 16 cluster with periodic-boundary conditions. We use ∆ AF /t = 0.001 to give a finite staggered magnetization in the trial wave function. This small value of ∆ AF is effective to pin a sizable value of the finite size order parameter m(0) because the single-particle states at U/t = 0 have a large degeneracy (∝ L) at the Fermi level,and are therefore strongly renormalized upon an arbitrary small ∆ AF . Since m(0) is finite [see Eq. (5)], m(β) remains finite even for finite L. Note that, at half filling, the HST in the charge channel with shift is equivalent to the one without shift because n i = 1. A very large equilibration time of ∼ 5000 Monte Carlo sweeps is found for m with the standard real HST coupled to the on-site electron spins. In this case, our shifted HST improves the equilibration time, allowing also higher acceptance rate (not shown) as in the doped cases, but the improvement is not really important. Amazingly, m is equilibrated almost immediately for the complex HST coupled to the on-site electron charges. This result implies that this pure-imaginary auxiliary field, which was first introduced by Hirsch [23], is very useful to estimate m at half filling for large U/t. Here we emphasize also that, not only the correlation time is highly reduced with this technique, but also fluctuations, thanks to this pinning strategy in the trial wave function, do not show any problem of large fluctuations, even at very large U/t and values. Figure 4 shows the finite-size scaling of m for U/t = 10. The cluster sizes used are L = 6, 8, 10, 12, 14, 16, 18, 20, 22, and 24. Here, the projection time β is chosen proportional to L, i.e., βt = αL, with α = 0.5, 1 and 1.5. These βt values are an order of magnitude smaller than those used with the pinning-field method [53,54], because in our approach we can reach the thermodynamic limit consistently without unnecessarily large values of β. Indeed, the extrapolated values at 1/β = 1/L = 0 are consistent for all values of α, which validates our approach. Our best estimate is obtained from the βt = 1.5L set of data, yielding m = 0.3046(1) in the ∆ → 0 limit, where the number in the parentheses indicates the extrapolation error in the last digit. Calculations are done for ∆ τ t = 0.2, 0.1, and 0.05 and the extrapolations to ∆ τ → 0 are obtained by a linear fit in (∆ τ t) 2 , determined by the leastsquares method. The ground-state expectation value m in the thermodynamic limit is obtained by extrapolating the results to L → ∞. In this case we fit the data in the range 6 L 24 with quadratic polynomials in 1/L. As it can be seen in Fig. 4, the time-discretization error is not negligible for U/t = 10. The extrapolated value is certainly smaller than the one in the Heisenberg model [55][56][57][58], where the latest Monte Carlo estimate is m = 0.30743(1) [59,60].
In Table I and Fig. 5 we show the values for m in the thermodynamic limit for U/t = 2, 4, 6, 8, 10 and 12 and compare them with the ones available in the literature [35,49,50,61]. The main outcome of this work is the estimated value of m for U/t 8, which is usually the accepted value for cuprates. Here, our error bar at U/t = 8 is two orders of magnitude smaller than the previous AFQMC estimate [1,49]. Thanks to this high statistical accuracy, our results clearly show that m increases monotonically in U/t. This is consistent with a strong-coupling expansion around the Heisenberg limit [62]. Here, finite-size scaling analyses are performed as follows. For U/t 6, the scheme of finite-size scaling analyses is the same as that for U/t = 10 which has been described before. For U/t = 4 (U/t = 2), cluster sizes up to L = 32 (L = 50) with twist-averaged boundary conditions [49,50,[63][64][65][66] are used because the finite-size effects are more important than those we have found at larger U/t values. A much larger value of ∆ AF /t = 10 is used for U/t 4 because the twists remove the degeneracy of the single-particle states at U/t = 0, as discussed before. All the results are obtained in the ∆ τ → 0 limit using data at ∆ τ t = 0.2, 0.1, and 0.05 for U/t 4 and ∆ τ t = 0.25, 0.2, and 0.1 for U/t = 2, respectively.

V. CONCLUSIONS AND DISCUSSIONS
In this work we have shown that, within the ground-state AFQMC technique, the choice of the trial function and the one for the auxiliary field are extremely important. In particular we have improved the efficiency of the method, by introducing shifted-discrete HSTs, that are useful for performing the imaginary-time evolution of symmetry-broken trial wave functions. The formalism can be considered as a generalization of the discrete HSTs in Ref. [23] or a compactification of the shifted-contour auxiliary-field Monte Carlo formalism in Ref. [38,39] specialized to the on-site Hubbard interaction.
Properly chosen auxiliary fields can improve the efficiency of AFQMC simulations. The shifted real auxiliary fields can improve the acceptance ratio, while the shifted pure- imaginary auxiliary fields can improve the average sign. The reason is that the shift in the real auxiliary field can stabilize the fluctuations of the norm of the determinant ratio, while the shift in the pure-imaginary auxiliary field can stabilize the fluctuations of the phase of the determinant ratio. However, even after the improvement, the average sign with the pureimaginary auxiliary field remains worse than the one obtained with real auxiliary field for the doped cases. Therefore, in the presence of the sign problem, the real auxiliary field is still recommended for achieving longer imaginary-time propagations. On the other hand, at half filling with large U/t, the pure-imaginary auxiliary field is dramatically more efficient than the real-auxiliary fields for evaluating the staggered magnetization m.
Finally, we remark on the d-wave superconducting order which has not been considered in the present study. It is noteworthy that an early study on a t-t -J model [68] has shown that a stripe state with spatially oscillating d-wave superconductivity is favored around 1/8 hole doping. Considering such an inhomogeneous superconductivity in a trial wave function might be of interest for a possible improvement of AFQMC simulations for doped Hubbard models with large U/t.