Numerical results on the short-range spin correlation functions in the ground state of the two-dimensional Hubbard model

TL;DR

This study uses advanced quantum Monte Carlo methods to compute short-range spin correlations in the ground state of the 2D Hubbard model, providing benchmarks for ultracold atom experiments and revealing how correlations evolve with interaction strength and density.

Contribution

It offers numerically exact and constrained path AFQMC calculations of spin correlations in the 2D Hubbard model, including finite doping, with optimized trial wave-functions and large supercells.

Findings

Nearest-neighbor spin correlation increases with U.

Next nearest neighbor correlation changes sign with density.

Results serve as benchmarks for ultracold atom experiments.

Abstract

Optical lattice experiments with ultracold fermion atoms and quantum gas microscopy have recently realized direct measurements of magnetic correlations at the site-resolved level. We calculate the short-range spin correlation functions in the ground state of the two-dimensional repulsive Hubbard model with the auxiliary-field Quantum Monte Carlo (AFQMC) method. The results are numerically exact at half filling where the fermion sign problem is absent. Away from half-filling, we employ the constrained path AFQMC approach to eliminate the exponential computational scaling from the sign problem. The constraint employs unrestricted Hartree-Fock trial wave-functions with an effective interaction strength U , which is optimized self-consistently within AFQMC. Large supercells are studied, with twist averaged boundary conditions as needed, to reach the thermodynamic limit. We find that the nearest-neighbor spin correlation always increases with the interaction strength U , contrary to the finite-temperature behavior where a maximum is reached at a finite U value. We also observe a change of sign in the next nearest neighbor spin correlation with increasing density, which is a consequence of the buildup of the long-range anti-ferromagnetic correlation. We expect the results presented in this work to serve as a benchmark as lower temperatures are reached in ultracold atom experiments.

Tables

Table 1. Table 1: Values of the spin correlation functions plotted in Fig.  2 for U = 2 , 4 , 6 , 8 , 10 𝑈 2 4 6 8 10 U=2,4,6,8,10 , and 12 12 12 , with distance 𝐫 = ( 0 , 1 ) , ( 1 , 1 ) , ( 0 , 2 ) 𝐫 0 1 1 1 0 2 \mathbf{r}=(0,1),(1,1),(0,2) , and ( 1 , 2 ) 1 2 (1,2) .

$U/t$	$𝐫=(0,1)$	$(1,1)$	$(0,2)$	$(1,2)$
$2$	$-0.0996\left(7\right)$	$0.0260\left(8\right)$	$0.0211\left(8\right)$	$-0.0294\left(7\right)$
$4$	$-0.1782\left(7\right)$	$0.091\left(3\right)$	$0.0811\left(8\right)$	$-0.0875\left(9\right)$
$6$	$-0.237\left(2\right)$	$0.137\left(4\right)$	$0.123\left(4\right)$	$-0.125\left(3\right)$
$8$	$-0.273\left(3\right)$	$0.162\left(9\right)$	$0.149\left(5\right)$	$-0.15\left(1\right)$
$10$	$-0.293\left(4\right)$	$0.180\left(5\right)$	$0.165\left(3\right)$	$-0.167\left(7\right)$
$12$	$-0.307\left(2\right)$	$0.190\left(6\right)$	$0.176\left(6\right)$	$-0.169\left(6\right)$

Equations

$$H=K+V=-\sum _{i,j,s}{t}_{ij}\left(c_{i,s}^{†}{c}_{j,s}+H.c.\right)+U\sum _i{n}_{i\uparrow }{n}_{i\downarrow },$$

$$c(\mathbf{r})=⟨{\psi }_g|\overrightarrow{S}(0)\cdot \overrightarrow{S}(\mathbf{r})|{\psi }_g⟩,$$

$$\overrightarrow{S}(\mathbf{r})=\frac{1}{2}\sum _{s{s}^{\prime }}{c}_{i,s}^{†}{\overrightarrow{\sigma }}_{s,s^{\prime }}{c}_{i,s^{\prime }},$$

$$\Psi (\dots ,{\mathbf{r}}_m+\mathbf{L},\dots )=e^{i\hat{\mathbf{L}}\cdot \mathbf{\Theta }}\Psi (\dots ,{\mathbf{r}}_m,\dots ),$$

$$|{\psi }_g⟩\propto \underset{\beta \to \infty }{\lim }{e}^{-\beta H}|{\psi }_0⟩$$

$$⟨O⟩=\frac{⟨{\psi }_0|e^{-\beta H}O{e}^{-\beta H}|{\psi }_0⟩}{⟨{\psi }_0|e^{-2\beta H}|{\psi }_0⟩}.$$

$$e^{-\beta H}=(e^{-\tau H}{)}^n=(e^{-\frac{1}{2}\tau K}{e}^{-\tau V}{e}^{-\frac{1}{2}\tau K}{)}^n+O({\tau }^2),$$

$$e^{-\tau U{n}_{\uparrow }{n}_{\downarrow }}=e^{-\tau U(n_{\uparrow }+n_{\downarrow })/2}\sum _{x=\pm 1}\frac{1}{2}{e}^{{\gamma }_{s}x(n_{\uparrow }-n_{\downarrow })},$$

$$e^{-\tau U{n}_{\uparrow }{n}_{\downarrow }}=e^{-\tau U(n_{\uparrow }+n_{\downarrow }-1)/2}\sum _{x=\pm 1}\frac{1}{2}{e}^{{\gamma }_{c}x(n_{\uparrow }+n_{\downarrow }-1)},$$

$$⟨O⟩=\frac{{\sum }_{\{X_i,X_j\}}⟨{\psi }_0|{\prod }_{i=1}^n{P}_i(X_i)O{\prod }_{j=1}^n{P}_j(X_j)|{\psi }_0⟩}{{\sum }_{\{X_i,X_j\}}⟨{\psi }_0|{\prod }_{i=1}^n{P}_i(X_i){\prod }_{j=1}^n{P}_j(X_j)|{\psi }_0⟩}$$

$$⟨O{⟩}_{\mathrm{mixed}}=\frac{{\sum }_k{w}_k⟨{\psi }_T|O|{\psi }_k⟩}{{\sum }_k{w}_k⟨{\psi }_T|{\psi }_k⟩},$$

$$c_L=c_{\infty }+a/L^3$$

$$c(0)=\frac{3}{4}(1-2⟨n_{\uparrow }{n}_{\downarrow }⟩).$$

$${\delta }_{\text{UHF}}=\frac{{\sum }_{\mathbf{r}}|c_{\text{UHF}}(\mathbf{r})-c_{\text{AFQMC}}(\mathbf{r})|}{{\sum }_{\mathbf{r}}|c_{\text{AFQMC}}(\mathbf{r})|}.$$

$${\delta }_{\text{CP}}=\frac{{\sum }_{\mathbf{\Theta }}|c_{\text{AFQMC}}^{\mathbf{\Theta }}-c_{\text{ED}}^{\mathbf{\Theta }}|}{{\sum }_{\mathbf{\Theta }}|c_{\text{ED}}^{\mathbf{\Theta }}|}.$$

Full text

I introduction

The Hubbard model Hubbard_origional is one of the most studied models in physics. The model is believed to be relevant to many correlated electron phenomena including interaction-driven metal-insulator transitions Imada_rmp_1998 , magnetism ph-sym and spin and charge density waves cdw_sdw , and most importantly high-temperature superconductivity htc . Except in one dimension lieb_wu , however, there is no analytic solution to the model. Numerical simulations have thus played an increasingly larger role in the study of the Hubbard model paper_simons .

The development in ultracold atom experiments provides another possibility for direct “simulation” of the Hubbard model rev_zoller ; rev_bloch . Recently, the two-dimensional Hubbard model was realized with ultracold atoms in optical lattices, and a flurry of activities have been reported where both local quantities and short-range (spin and charge) correlations were measured nat_519_2015 ; prl_116_175301 ; prl_116_235301 ; sci_1253 ; sci_1257 ; sci_1260 . Numerical results have usually been used to benchmark the experimental results. In addition to serving as a thermometry for ultracold atoms, computational results have been integrated in these studies to guide the experiments and interpretation.

Understanding the properties of the doped Hubbard model is challenging, because of the existence of different competing orders for the ground state which are separated by tiny energy scales stripe . High accuracy and resolution is required to distinguish and characterize the different candidate states. To date, the temperature that can be reached by ultracold atom experiments is still relatively “high” compared to the energy scales of the competing ground-state orders. In this regime, reliable numerical results have been provided by a multitude of computational methods including finite temperature determinant quantum Monte Carlo (DQMC) prl_104_066406 ; prl_106_035301 , dynamical cluster approximation (DCA) prb_88_155108 , and numerical linked-cluster expansion (NLCE) pra_84_053611 ; prl_109_205301 , all of which work at finite temperatures.

As experiments continue the development of cooling technology, we can expect lower temperatures to be reached in the near future. In fact a lower temperature of $T=0.25(2)t$ has already been reported in the most recent experiment Mazurenko_2016 . This offers the exciting prospect of determining and understanding the low-temperature and ground-state phases of the Hubbard model by the optical lattice experiments. Reliable numerical data at lower temperatures and ground state will be crucial for benchmarking experimental results and assisting interpretation and analysis. Reaching lower temperatures present significant challenges paper_simons for numerical methods, however, and numerical data will be less readily available or reliable.

In anticipation of these developments and challenges, we calculate the short-range spin correlation functions of the 2D Hubbard model at zero temperature. Applying the state-of-the-art auxiliary-field quantum Monte Carlo (AFQMC) method, we study systems with large size and take advantage of twist averaged boundary conditions to reach the thermodynamic limit. At half-filling, where the sign problem sign_1 ; sign_2 is absent, the results are numerically exact. For doped systems, we employ the constrained-path (CP) approximation zhang_prb_1997 to deal with the sign problem. In CP-AFQMC the sign problem is eliminated by a constraint on the random walks in Slater determinant space which is dependent on the trial wave-function. Previous studies have shown that this bias is usually small chia-chen_EOS ; CPMC_sym_1 . A further recent advance allows a self-consistent optimization of the trial wave-function prb_self-consistent . For each filling and interaction strength, we test a series of trial wave-functions generated from unrestricted Hartree-Fock (UHF) with different $U$ values, relying on the feedback from the CP-AFQMC calculation to self-consistently determine the optimal effective $U$ value, $U_{\rm eff}^{\rm UHF}$, for the UHF. As a result, the computational approach we use can treat sufficiently large system sizes in the ground state, and obtain highly accurate results paper_simons ; stripe

We find major differences in the behavior of the short-range magnetic correlations from what has been observed to date experimentally and computationally. At half-filling, the spin correlation always increases with the interaction strength $U$ and monotonically approaches the Heisenberg limit with $U\rightarrow\infty$. This behavior of the spin correlation is in sharp contrast with the finite-temperature situation, where the correlation reaches its maximum at a finite $U$. For the more important doped cases, we observe a change of sign in the next nearest neighboring (NNN) correlation functions which is a precursor for the onset of spin-density wave orders and the buildup of the anti-ferromagnetic correlation with increase of density. We find that the filling factor where the NNN correlations function changes sign is essentially independent of interaction strengths.

The rest of the paper is organized as follows. In Sec. II, we first define the Hubbard model and give a brief summary of the method used in this work. We will comment on the approach to the thermodynamic limit with twist averaged boundary conditions (TABC), and mention some of the latest technical advances employed in the study to improve accuracy and computational capabilities and efficiency. In Sec. III we present the results on spin correlations at half-filling as a function of $U$ values. In Sec. IV, the nearest and next nearest neighboring spin correlation functions for systems away from half-filling are presented. A short summary in Sec. V will conclude this paper.

II Model and Method

II.1 Hubbard Model

The Hubbard model is defined as

[TABLE]

where $K$ and $V$ are the kinetic and on-site interaction terms, respectively. The creation (annihilation) operator on site $i$ is $c_{i,s}^{\dagger}$ ($c_{i,s}$), with $s=\uparrow,\downarrow$ the spin of the electron, and $n_{i,s}$ is the corresponding number operator. We denote the total number of electrons with $\uparrow$- and $\downarrow$-spin by $N_{\uparrow}$ and $N_{\downarrow}$, respectively. In this work, only spin-balanced systems ($N_{\uparrow}=N_{\downarrow}$) are considered. The filling factor is defined as $n=(N_{\uparrow}+N_{\downarrow})/{N}$ with $N$ being the total number of lattice sites in the supercell. So $n=1$ means half-filling. We will typically use supercells of square lattice with size $N=L\times L$. We only consider nearest neighbor and uniform hopping in this work, $t_{ij}=t$ for each near-neighbor pair $\langle ij\rangle$, and $t$ is set as the energy unit. The strength of the repulsive interaction is given by $U/t$.

The quantity we investigate in this work is the spin correlation function,

[TABLE]

where $|\psi_{g}\rangle$ is the ground state of $H$ in Eq. (1). The spin operator at $\mathbf{r}$ is given by

[TABLE]

where $i$ is the lattice site label of the position $\mathbf{r}$, and $\overrightarrow{\sigma}$ are the Pauli matrices. Because of translational invariance, the reference site ‘$\mathbf{0}$’ can be averaged over in the periodic supercell, and the correlation function is only a function of lattice vector ${\mathbf{r}}$ and satisfies all lattice symmetries.

In order to extrapolate more reliably to the thermodynamic limit (TDL), we adopt twist averaged boundary conditions TBC . Under twist boundary conditions (TBC), an electron gains a phase when hopping across the boundaries:

[TABLE]

where $\hat{\mathbf{L}}$ is the unit vector along $\mathbf{L}$, $\mathbf{r}_{m}$ is the position of the $m$-th electron, and the twist angle $\mathbf{\Theta}=(\theta_{x},\theta_{y})$ is a two dimensional parameter, with $\theta_{x}$ ($\theta_{y}$) $\in[0,2\pi)$. This is equivalent to placing the lattice on a torus and applying a magnetic field which induces a flux of $\theta_{x}$ along the $x$-direction (and a flux of $\theta_{y}$ along the $y$-direction). In Eq. (4), the translational symmetry is explicitly broken, however it can be easily restored with an alternative gauge, i.e., adjust $t$ to $t\times e^{i\theta_{x}/L}$ along $x$ ($t\times e^{i\theta_{y}/L}$ along $y$).

To implement TABC, a set of $N_{\theta}$ twist angles is chosen and we carry out an independent calculation for each twist. The final result of a physical quantity is the averaged value from these independent calculations. The constrained path condition can be generalized straightforwardly to the case of TBC chia-chen_EOS . By imposing a random TBC, the possible degeneracy of the non-interacting energy levels is lifted by explicitly breaking the rotational symmetry of the lattice. This eliminates the so called open-shell effects and helps to reduce the constraint bias when a free-electron trial wave function is used. In most of our calculations, we use UHF trial wave functions so that this is not especially relevant. Nevertheless we have retained the same procedure.

When implementing the TABC, we use a quasi-random instead of a pseudo-random sequence to generate the twist. As shown in Ref. prb_benchmark , the TABC physical quantities converge faster (with respect to the number of twists used) with quasi-random twists than with a pseudo-random sequence. A quasi-random sequence is distributed more uniformly in the sampled space. Twists from a uniform grid have essentially the same convergence rate as a quasi-random sequence, however it can lead to open-shell situations in the many-body calculation, which is less advantageous especially when a simple trial wave function is used prb_benchmark . Additionally, a uniform grid is typically not cumulative which means that, if the targeted statistical accuracy is not reached with a given set, we would likely not be able to re-use the data, and new calculations would be necessary for each twist in the new set. A quasi-random sequence avoids these problems and combines the advantages of a pseudo-random sequence and a uniform grid.

II.2 Auxiliary-field quantum Monte Carlo method

In this section, we will briefly introduce the methods used in this work. (A more comprehensive discussion of this method can be found in Ref. lecture-notes .) Our AFQMC calculations share much of the same basic framework as standard DQMC method AFQMC and its ground-state variant koonin . By repeatedly applying the projection operator to a state $|\psi_{0}\rangle$ with non-zero overlap with the ground state $|\psi_{g}\rangle$ of the Hamiltonian $H$ in Eq. (1), we can obtain $|\psi_{g}\rangle$:

[TABLE]

and the expectation value of an operator $O$ can be represented as

[TABLE]

Through the Trotter Suzuki decomposition, the kinetic and interaction parts in the projection operator can be decoupled as:

[TABLE]

where $\beta=\tau n$. The Trotter error can be eliminated by an extrapolation of $\tau$ to [math]. We typically choose $\tau=0.01$ in this work. We have verified that the Trotter error with $\tau=0.01$ is smaller than the targeted statistical errors.

The initial state $|\psi_{0}\rangle$ is usually chosen as a Slater determinant in AFQMC. The one-body term $e^{-\frac{1}{2}\tau K}$ can be directly applied to it and the result is another Slater determinant. This is not true for the two-body term $e^{-\tau V}$. However, the two-body term can be decomposed into an integral of one-body terms through the so-called Hubbard-Stratonovich (HS) transformation. Different types of HS transformations TREMBLAY_1992 for $e^{-\tau V}$ exist in literature. The two commonly used types are the so called spin decomposition

[TABLE]

with the constant $\gamma_{s}$ determined by $\cosh(\gamma_{s})\equiv\exp({\tau}U/2)$, and the charge decomposition

[TABLE]

with $\cosh(\gamma_{c})\equiv\exp(-{\tau}U/2)$ hirsch_prb_1983 . Here $x$ is an Ising-spin-like auxiliary field. Transformations with continuous Gaussian auxiliary-fields exist which can be made CPMC_sym_1 essentially as efficient as the discrete forms. Different choices of the HS can lead to different accuracies or efficiencies, because of symmetry considerations CPMC_sym_1 ; CPMC_sym_2 or other factors Hao-inf-var .

With the HS transformation, Eq. (6) turns into

[TABLE]

where $X_{i}$ is the collection of the $N$ auxiliary fields from the HS transformation, and $P_{i}$ is the product of the kinetic term and the one-body terms from the HS transformation at time slice $i$. The multi-dimensional integrals can then be evaluated by Monte Carlo methods, e.g., with the Metropolis algorithm.

At half-filling, each individual term in the sum in the denominator of Eq. (10) is always non-negative, because of particle-hole symmetry ph-sym . So we can use it as probability density and the sign problem is absent. Our calculations in this work use mostly the path integral form outlined above, but introduce several recent algorithmic advances such as an acceleration technique hao_pra_2015 (with force bias lecture-notes ; zhang_prl_2003 in the Metropolis sampling) and control of the divergence of Monte Carlo variance Hao-inf-var .

Away from half-filling, a direct evaluation of Eq. (10) by Monte Carlo will suffer from the sign problem sign_1 ; sign_2 , since terms in the denominator of Eq. (10) can become negative for some auxiliary fields lecture-notes . The sign problem can be eliminated by the constrained path approximation. The framework within which this has been implemented in Hubbard-like model has often been referred to as the constrained path Monte Carlo (CPMC) method zhang_prb_1997 . Here we will refer to the approach as CP-AFQMC to be consistent with recent conventions. A description of the CP-AFQMC method as applied to Hubbard-like models can be found in Ref. Huy_CPC_2014 .

In CP-AFQMC, the wave function is represented as an ensemble of a set of Slater determinants which are called walkers. The evolution of wave function in the imaginary time is represented as random walks in the Slater determinant space by sampling the auxiliary fields. The paths of auxiliary-fields are constrained to ensure the overlap of any propagated walker with the trial wave-function, $|\psi_{T}\rangle$, computed at each time slice, remains non-negative. Physical quantities can be evaluated using the mixed estimate as

[TABLE]

where $|\psi_{k}\rangle$ is the kth walker, and $w_{k}$ is the corresponding weight. The mixed estimate is used to compute the energy (and other observables which commute with the Hamiltonian). For observables which do not commute with the Hamiltonian, the mixed estimate is biased, and we use back propagation to correct for this zhang_prb_1997 ; Wirawan-PRE .

As mentioned, the CP-AFQMC overcomes the sign problem and restores a low computational scaling with system size, at the cost of a systematic error which depends on the trial wave-function $|\psi_{T}\rangle$. Previous studies have shown the systematic error is small even with a free-electron or Hartree-Fork trial wave-function chia-chen_EOS . In this work, we adopt a recent advance prb_self-consistent which allows the self-consistent construction of an optimal optimal UHF trial wave-function with an effective interaction $U_{\rm eff}^{\rm UHF}$. In the self-consistent procedure, we couple the CP-AFQMC calculation with a mean-field calculation. We first carry out an ordinary CP-AFQMC calculation with a free electron or UHF trial wave-function. Then instead of solving the mean-field Hamiltonian self-consistently, we feed the local density from the CP-AFQMC calculation as the “field” in the mean-field Hamiltonian and scan the interaction $U$ to find the solution which yields a local density closest to the input density. Next we use the mean field solution with this “optimal” $U$ as the trial wave-function for the next-step CP-AFQMC calculation. The same process is repeated until the local density is converged.

The use of the self-consistent CP-AFQMC further improved the accuracy, especially for determining spin- and charge-orders prb_self-consistent . We have carried out additional benchmarks here in smaller systems specifically targeting short-range correlations. Therefore results presented here, consistent with previous experience paper_simons ; stripe , are expected to be very accurate.

III Results at half-filling

In this section, we present results for and discuss the half-filling case. Following the procedure in Ref. prb_benchmark , we take advantage of the TABC to remove errors from the finite size effect from the use of supercells. Lattice size up to $22\times 22$ are studied, which is sufficient for convergence. In Fig. 1, we show the results of the nearest neighbor (NN), i.e., ${\mathbf{r}}=(0,1)$, spin-spin correlation. It is seen that the NN correlation is already converged to within $1\%$ with a $14\times 14$ lattice. We also perform a finite-size fit of the results. Since the NN spin correlation function (in the infinite large U case) is actually given by the energy in the Heisenberg model on a square lattice, we fit the short-range spin correlation function using the result from spin-wave theory finite_size_scaling_C :

[TABLE]

where $c_{L}$ is a spin correlation function defined in Eq. (2) for system with size $L\times L$ and $c_{\infty}$ is the thermodynamic value. The quality of the fits, as shown in Fig. 1, indicate that the finite-size effects in the NN spin correlation function is indeed captured well by the form in Eq. (12).

The spin correlation functions for all ${\mathbf{r}}$ between $(0,1)$ and $(1,2)$ at thermodynamic limit are plotted in Fig. 2. Spin correlations at large distances, in the context of determining the long-range order, have been computed in Ref. prb_benchmark . From Fig. 2, we see that the spin correlation functions always increase as $U$ is increased. The NN spin correlation function approaches the Heisenberg value, $-0.334718(3)$ sandvik_1997 when $U$ approaches infinity. Also at large $U$ the square root of the correlation function in the infinite distance limit should approach the value of magnetization in the 2D Heisenberg limit ($0.3070(3)$ in Ref. sandvik_1997 ). We also list all the spin correlation function values in Table 1.

We observe that the behavior of spin correlation functions at zero temperature is different from that at finite temperatures where thermal fluctuation is present. For comparison, we include in Fig. 2 the NN spin correlation at $T=0.31$ from Ref. prl_106_035301 . While the $T=0\,$K result increases monotonically with $U$, the $T=0.31$ NN spin correlation reaches a maximum value at a finite $U$ ($\sim 8$) prl_106_035301 ; pra_84_053611 . The difference results from the competition between quantum and thermal fluctuations. At $T=0$ K, with no thermal fluctuations, the increase of $U$ always drives the system towards the Heisenberg limit where no double occupancy is allowed. At a given non-zero temperature, however, the effective Heisenberg anti-ferromagnetic coupling MacDonald_prb_1988 , $t^{2}/U$, decreases with $U$. (This implies that, in the infinite-$U$ limit, the corresponding Heisenberg model would be at infinitely high temperature.) So we expect that the value of $U$ where the NN spin correlation function reaches its maximum will increase as the temperature $T$ is decreased, and become infinitely large at the limit of $T=0$ K. This behavior should be easy to confirm experimentally as lower temperatures are reached with ultracold atoms in optical lattices.

The value of the on-site spin correlation function at half-filling can be inferred from the double occupancy results in Ref. prb_benchmark through the following:

[TABLE]

Unlike the spin correlation function, the double occupancy always decreases with increasing $U$. Correspondingly, the on-site spin correlation function always increases with interaction, both at zero and finite temperatures.

IV Results away from half-filling

IV.1 Self-consistent constraint and error quantification

In order to control the sign problem we a constraint on the paths of the importance-sampled random walks in Slater determinant space, as mentioned in Sec. II.2. A UHF trial wave function is used, which is generated with an effective interaction, $U_{\rm eff}^{\rm UHF}$, via a self-consistent iteration with the AFQMC calculation prb_self-consistent . The optimal $U_{\rm eff}^{\rm UHF}$ is the value of effective $U$ with which the corresponding UHF solution yields results closest to those from AFQMC (using the UHF as constraining trial wave-function and always performed with the physical $U$). To measure “closeness” we seek to minimize the following metric in the present work:

[TABLE]

Since we treat periodic supercells and the UHF solution breaks translational symmetry, a subtlety arises in using the UHF as a trial wave-function and also in comparing the results as in Eq. (14) above. One can think of this as applying a small pinning field as was the situation in Refs. prb_self-consistent ; stripe . In the TDL one would expect the effect on the short-range correlations to be small. This is supported by our finding below that the optimal $U_{\rm eff}^{\rm UHF}$ is slightly larger for smaller supercells but approaches the results in Ref. prb_self-consistent in large supercell sizes. It is important to note that the AFQMC results are insensitive to variations in the trial wave function that result from small changes in the value of $U_{\rm eff}^{\rm UHF}$.

In Fig. 3, we show the $8\times 8$ system with $U=6,8$ and $12$ and $n=7/8$ as an example for the search of the effective $U$ for the UHF trial wave-function. The optimal $U_{\rm eff}^{\rm UHF}$ is seen to fall between $2.5$ and $3.0$ which is close to the results ($U_{\text{eff}}\sim 2.7$ for $U=8$) in Ref. prb_self-consistent , where fully self-consistent calculations were performed. We apply this procedure to all other fillings and $U$ values and determine the corresponding optimal $U_{\rm eff}^{\rm UHF}$ in the same way. We then study systems with larger supercell sizes using the UHF trial wave-functions with the same corresponding $U_{\rm eff}^{\rm UHF}$. As correlation effects are diminished in systems with small $U$ and at very low filling factors, the optimal trial wave-functions are found to become the free electron wave-functions ($U_{\rm eff}^{\rm UHF}=0$).

We also carried out an additional benchmark of CP-AFQMC, comparing results with exact diagonalization (ED) in the worst case scenario of the $4\times 4$ lattice, with $U=12$ and $N_{\uparrow}=N_{\downarrow}=7$. Results are shown in Fig. 4. We plot the NN spin correlation functions for $7$ quasi-random twist angles, where the CP-AFQMC calculations were performed with UHF trial wave-functions generated with different $U_{\rm eff}^{\rm UHF}$. We calculate the relative absolute error of the NN spin correlation with respect to the ED values as:

[TABLE]

From the inset of Fig. 4, we see that $U_{\rm eff}^{\rm UHF}\sim 3.5$ yields a minimum $\delta_{\text{CP}}$, of $\sim 5\%$. This is consistent with the optimal $U_{\rm eff}^{\rm UHF}$ value determined by the minimization of $\delta_{\text{UHF}}$ in Fig. 3. As mentioned above, it is reasonable in a smaller system that $U_{\rm eff}^{\rm UHF}$ is slightly larger, due to a need to overcompensate for lack of dynamic correlations. It should be noted that the final TABC result will have a smaller error, because of error cancellations from different twist angles. (For the system shown in Fig. 4, it is 3%.) The results in this system provide a kind of upper bound estimate to the CP bias. Lower $U$ and other doping parameters all make the method much more accurate. Larger supercell sizes are also expected to reduce the sensitivity of the short-range correlations.

IV.2 Results

In all our calculations we study sufficiently large supercells to ensure that finite-size effects are negligible in the computed spin correlation functions. In Fig. 5, we plot the NN and NNN spin correlation functions for $U=4$ and $8$ at a filling factor of $n=0.75$. Linear fits of the results with $1/L$ are also shown. It can be seen that, in these systems, the finite size effect is smaller than the targeted statistical accuracy even at $12\times 12$.

In Fig. 6, we plot the NN and NNN spin correlation functions for interaction strengths ranging from $U=2$ to $U=12$. The negative sign of the NN spin correlation function reflects the short-range anti-ferromagnetic correlations in the Hubbard model. The dependence of the correlation on $U$ is mild at the dilute limit, which is reasonable since double occupancy is significantly reduced. As the density is increased, stronger dependence of the correlation on $U$ is seen from weak to moderate interactions. For even larger interaction strengths, the NN spin correlation functions approach the value of infinite-$U$ limit where no double occupancy is allowed.

For the NNN spin correlation function, the sign is negative from the dilute limit through $n\sim 0.8$. This is a reflection of the exchange-correlation hole which continues from on-site to NN to NNN correlations and so on. Similar to the NN spin correlation function, it only shows mild dependence on the interaction strength in the dilute limit. As the filling factor is further increased, the NNN spin correlation changes sign from negative to positive, which is a precursor of the buildup of anti-ferromagnetic order in the system. Interestingly, the crossover points in density are very close for different interaction strengths. The change of sign in the NNN spin correlation function, needless to say, does not necessarily imply long-range order in the systems. To establish the existence of long-range anti-ferromagnetic orders, a more systematic examination of the behavior of the tail of $c({\mathbf{r}})$ is necessary (see, e.g. Ref. prb_benchmark for analysis at half-filling).

Previous AFQMC calculations have shown the existence of spin-density wave states at intermediate interaction ($U\sim 4$) for density larger than $\sim 0.85$ chang_prl_2010 . For larger interactions they evolve into stripe states chang_prl_2010 ; stripe . The determination of such phases requires the study of long-range spin correlation functions chang_prl_2010 or long-range density and spin-density variations in the presence of pinning fields stripe , as well as careful removal of the influence of finite size and supercell shape chia-chen_EOS ; chang_prl_2010 ; stripe . The effect of these collective modes in the ground state is less direct in the short-range correlations computed here, although clearly they can lead to quantitative modifications (reductions) in the magnitude of NN and NNN spin correlations. The slopes of the curves in Fig. 6 show large increases in magnitude as we approach half-filling, which is consistent with and could be a manifestation of these states.

V Summary

In summary, we have calculated the short-range spin correlation functions in the two-dimensional Hubbard model at zero temperature with AFQMC. At half-filling, the results are numerically exact. The absolute values of spin correlation functions are found to always increase with $U$, which is different from the finite temperature behavior. Away from half-filling, we eliminate the sign problem with a self-consistent constraint. The systematic errors from the constraint are examined and quantified. We observe a change of sign from negative to positive in the next nearest neighbor spin correlation function as a function of doping, with a crossing point slightly below $0.8$ which shows very weak dependence on the interaction strength. The results in this paper can serve as a valuable benchmark as optical lattice experiments with ultra-cold atoms reach lower temperatures. The detailed and quantitative results on spin correlation functions provide useful information in the search for theoretical understanding and experimental realization of exotic phases of magnetic and accompanying or competing charge and possibly superconducting orders.

Acknowledgements.

We acknowledge support from NSF (DMR-1409510). MQ and SZ were also supported by the Simons Foundation. The calculations were carried out at the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575, and the computational facilities at the College of William and Mary. We gratefully acknowledge a Director’s discretionary allocation at OLCF.