Automatic Order Detection and Restoration Through Systematically Improvable Variational Wave Functions

Variational wave function ansatze are an invaluable tool to study the properties of strongly correlated systems. We propose such a wave function, based on the theory of auxiliary fields and combining aspects of auxiliary-field quantum Monte Carlo and modern variational optimization techniques including automatic differentiation. The resulting ansatz, consisting of several slices of optimized projectors, is highly expressive and systematically improvable. We benchmark this form on the two-dimensional Hubbard model, using both cylindrical and large, fully periodic supercells. The computed ground-state energies are competitive with the best variational results. Moreover, the optimized wave functions predict the correct ground-state order with near full symmetry restoration (i.e. translation invariance) despite initial states with incorrect orders. The ansatz can become a tool for local order prediction, leading to a new paradigm for variational studies of bulk systems. It can also be viewed as an approach to produce accurate and systematically improvable wave functions in a convenient form of non-orthogonal Slater determinants (e.g., for quantum chemistry) at polynomial computational cost.


I. INTRODUCTION
Understanding the ground state properties of strongly correlated systems has been an important goal in the study of electronic systems.There has been much success in utilizing variational methods, such as variational Monte Carlo (VMC) [1,2] and density matrix renormalization group (DMRG) [3][4][5], however each has drawbacks.Tensor network methods, while having incredible variational power, often have challenges of converging large width cylinders or periodic boundary conditions with finite resources [6].In order to understand these large bulk systems, quantum Monte Carlo (QMC) methods [7,8] are invaluable, however interesting regimes can be plagued with a sign problem, which renders many variational calculations exponentially difficult.Mitigation of this problem, for example using constrained path auxiliary field QMC (CP-AFQMC) [9], can provide accurate local energies but require more sophisticated methods, such as back-propagation, to measure observables that do not commute with the Hamiltonian.
In VMC, a parameterized wave function ansatz |ψ⟩ is optimized against the variational energy expectation ⟨E⟩ = ⟨ψ| Ĥ|ψ⟩ / ⟨ψ|ψ⟩.For many common ansatze, this evaluation takes place by stochastic sampling via Monte Carlo in an occupation basis.The key to this approach is the interplay between expressivity of the ansatz, and the ability to optimize it to a global minimum.Often the limiting factor has been that the commonly used forms for fermionic systems lack sufficient expressivity to capture different quantum phases under the same wave function ansatz, thereby making it difficult to detect orders in an unbiased manner.
Recently there has been a number of fermionic wave functions parameterizations introduced that are inspired * rlevy@flatironinstitute.org from machine learning algorithms [10][11][12][13].These are generally systematically improvable and have been shown to obtain very close agreement with the true ground state in small to medium-sized systems.However, in order to optimize vast numbers of parameters, the state may need to be constrained to certain physics, e.g.superconductivity or the wavelength of a charge density wave [12,14].Furthermore, a wave function ansatz with many parameters can make it more challenging to connect with the underlying physics.In strongly correlated systems with competing orders being separated by small energy differences, the variational energy can cease being an effective signal for detection of different phases.An ansatz with a smaller number of parameters and a simpler, more physically-inspired construction can be advantageous.
We present a variational wave functions ansatz based on the theory of auxiliary fields, which aims to have O(N 2 ) or less optimizable variational parameters and yields a simple form of the wave function that is straightforward to systematically improve.The ansatz builds on a form originally introduced in Ref. [15] by Sorella, as a modified version of auxiliary field quantum Monte Carlo (AFQMC), and given the name variational AFQMC, or VAFQMC.We introduce a further generalized variational form, which defines a framework with which the expressivity can be easily expanded.An extension was recently presented to treat molecules in Ref. [16], with a different form than what we consider.We benchmark our variational wave function in the Hubbard model against existing VMC results, as well as DMRG and two projector QMC methods: AFQMC [9,17] and fixed-node diffusion Monte Carlo [14,18].
The rest of the paper is organized as follows.In section II, we introduce our variational formulation and describe the implementation details.Then in section III we benchmark against state of the art VMC, projector QMC, and tensor network methods, showing competitive variational energies and predictive power of local observ- ables.Finally in section IV we conclude and discuss the outlook of VAFQMC in conjunction with other methods.

II. METHODS
We consider a variational form modeled after the projector method [19] such that the ground state |ψ 0 ⟩ of a Hamiltonian H can be obtained via for some trial state |ψ T ⟩ which has non-zero overlap with the ground state.Traditionally, as in AFQMC, a stochastic random walk is formulated to project the trial wave function into the ground state via a Monte Carlo routine.A cartoon representation of this process is depicted in Fig. 1a.Instead, by fixing the number of projection steps and optimizing the projectors as a wave function ansatz, we can produce an integral representation of the variational state shown in Fig. 1b.Our formalism will apply to all Hamiltonians containing general one-body and two-body terms.For concreteness we will use the Hubbard model [20] in this work to both describe our approach and perform benchmark tests.The Hubbard model is given by where ⟨ij⟩ denotes nearest-neighbor hopping, U/t is the strength of the on-site interaction.We restrict the model to 2D such that i = (i x , i y ), and consider lattices of size N x × N y .Our focus is on describing ground-state order and approaching the 2D limit.For this we study both cylindrical geometry, with open boundary condition along x and periodic boundary condition (PBC) along y, and square lattices with fully PBCs along both direction.The cylindrical cells allow us to perform systematic benchmarks on the description of spin and charge density waves, the primary order in the ground state of this model.The fully periodic supercells allow us to test the quality of the wave function -including symmetry restoration -and scaling as we approach the 2D bulk limit.
To apply the projection operator exp(−τ Ĥ) to the trial wave function, in AFQMC the propogator is split up into a kinetic and potential part exp(−τ Ĥ) ≈ exp(−τ T ) exp(−τ V ).This step incurs an error of O(τ ) and must be corrected, but this is not of concern here due to the variational nature of our wave function as dis-cussed below.The kinetic term can be applied to the trial wave function as a transformation of the Slater determinant, via Thouless theorem [21].For the potential term, a Hubbard-Stratonovich (HS) decomposition can be applied [22,23], for example the discrete version where cosh λ i = exp(U τ ).The form of our variational wave function can be thought of as a sum over auxiliary fields: where T l σ = t lσ ij c † σi c σj and constant terms are ignored.This can be recast into an effective Hamiltonian form as where Ĥl is determined by optimizable parameters T l σ , τ l , as well as the form of the HS, all of which can be different from that of the actual Hamiltonian H.In other words, we consider variational wave functions of the form where B l (x l ) is a one-body propagator of the form B l (x l ; α) = exp[ ij σσ ′ h l iσ;jσ ′ (x l ij ; {α})c † i,σ c j,σ ′ ], with i denoting lattice sites (or basis functions), |x⟩ is the resulting application of B l (x) |ϕ 0 ⟩, and σ (σ ′ ) denoting spins.The effective one-body "hamiltonian" in the exponent in B l contains optimization parameters, denoted by α.For the Hubbard model example we consider here, B l (x l ) is of the simpler, product form in Eq. ( 4), with α consisting of α = {λ i , τ l σ , t lσ ij }.Here we have separated h l in the one-body auxiliary-field-independent term into an overall scaling factor τ l and the hopping terms T l iσ;jσ , so that no overall scaling factor is involved in the optimization of the latter.These parameters are unconstrained: λ i is repeated for each projector slice l but can be non-uniform over sites i; τ l σ is now distinct from λ i and is allowed to be positive or negative; and T l is no longer constrained to be Hermitian.In practice the values of these parameters are optimized directly.Most of these parameters are in T l σ , with O(N ℓ • N 2 ) parameters.This general variational ansatz, which is restricted to one-body terms, can also be related to shadow wave functions where a kernel of both physical and auxiliary degrees of freedom is parameterized and auxiliary variables are integrated out [24,25].Shadow wave functions however, have historically been used in a particle configuration basis unlike Eq. ( 6).Unlike Ref. [15], our variational form has no fixed projection time . This allows the ansatz to 'grow' as needed in projection time, to potentially reach the ground state with few effective projector steps N ℓ .In fact, completely decoupling τ l from λ i allows for not only cooling, i.e. β l+1 = β l + τ , but intermediate heating i.e. β l+1 = β l − τ .We show below that the optimized results indeed takes advantage of this structure.
Given a variational form of Eq. ( 4) or generally Eq. ( 6), we can evaluate the energy via a double integral where each x (or x ′ ) is a field configuration that produces a Slater determinant |x⟩ The probability distribution | ⟨x ′ |x⟩ | can be sampled using Markov chain Monte Carlo (MCMC) through conventional Metropolis algorithm methods [8] where each state is a configuration of fields (x, x ′ ).This procedure is distinctly different from VMC however, and is instead closer to the evaluation of a N ℓ -slice path integral in imaginary time.Note that because ⟨x ′ |x⟩ is not necessarily positive semi-definite a sign problem can occur.However, for the problems we consider, the number of slices, N l , is modest and the average ⟨S⟩ remains positive as further discussed below.Energy optimization is carried out through gradient descent.The derivation of the derivative is provided in Ref. [15] which we repeat here for convenience.The derivative of of parameter α is given by where O = ∂ ln(Sρ) ∂α .Both ∂E L ∂α and O can be easily evaluated using automatic differentiation.By sampling ρ and accumulating the values in eq.( 10), optimization of the wave function can be performed.This optimization procedure of VAFQMC is therefore dominated by both matrix multiplication (the projection of |ϕ 0 ⟩) and determinant calculations for the local energy which implies VAFQMC scales as ≈ O(N ℓ • N 3 ).We note that standard practices in AFQMC, for example, compact decomposition of long-range interaction (or Jastrow), the use of force bias [26] or other update schemes, stabilization, low-rank decomposition, fast computations of local energy etc [27], can all be applied to this method.

A. Implementation Details
We consider the standard Hubbard model with only nearest neighbor hopping.Two types of calculations are considered, using two different boundary conditions: fully periodic (PBC) and cylindrical boundary conditions (CBC).When using cylindrical boundaries, we apply a weak edge pinning field H pin = i∈edge (−1) ix+iy h Ŝz i where Ŝz i = (n ↑i − n↓i )/2 to induce a local antiferromagnetic order, where h controls the strength of pinning and is chosen as h = 1/2 unless otherwise specified.When pinning is present, energies are reported including the contribution from the additional field, consistent with other publications [28].
Our implementation utilizes Jax [29] to perform automatic differentiation and sampling, alongside the Optax package [30] for optimization.The matrix exponential is computed approximately using a k = 3 product formula [31].Optimization generally uses around 180 optimization steps with varying and decaying step size (see appendix); gradients are accumulated with 189,000 samples and parameters are updated using a random step with the computed sign [11].
This work focuses on the case of N ℓ = 4, which gives substantial improvement over e.g.N ℓ = 2 or a single slice.While any given calculation can always be reoptimized with additional slices, by fixing N ℓ we can observe the limits of this form.
We can additionally add a term to optimize the sign similar to Ref. [16], however in the applications presented here the average sign is almost always very large ⟨S⟩ ≥ 0.93.For U = 8, the sign can get as low as ⟨S⟩ ≈ 0.74 but this is sufficiently large to resolve observables and derivatives without needing any further constraint in the optimization.We renormalize (i.e.'clip') both the sign term derivative and the energy derivative term to be on the same order.An example of the optimization is shown in Fig. 2.
After optimization, parameters are selected based on the lowest observed average energy.Using this parameter set, observables such as the local energy and electron occupation are measured with approximately 10 6 samples with a thermalization of 10 3 sweeps.In order to properly treat the 1/⟨S⟩ term, error bars are reported from a jackknife estimate [8].

A. Energetic Benchmarks
To test the variational ansatz, we first benchmark against doped large square Hubbard models with PBC as shown in Table I.Some systems at intermediate doping have a reference value from AFQMC [17] which is often effectively exact; we compare to AFQMC the following systems: 6 × 6 with n = 2/3, 8 × 8 with n = 0.6875, and 10 × 10 with n = 0.8 at U/t = 4.The VAFQMC ansatz is able to obtain an energy error that is O(10 −3 ) with an energy variance per site (⟨H 2 ⟩ − ⟨H⟩ 2 )/N of less than 1, which improves over traditional variational approaches in such systems.
For systems without exact energies, we compare with several state of the art methods, projector and also variational QMC methods.Note that of the projector methods, AFQMC and fixed-node projector diffusion Monte Carlo (DMC) with a VMC optimized trial wave function, only fixed-node is variational.Then we compare with VMC energy and energy variance.Additionally note that the fixed-node DMC energies of Ref. [14] are higher than those of Ref. [18]; we report the lowest energy available.
Overall, the variational energy of the VAFQC ansatz is competitive or improves upon those from state-of-the-art VMC under the standard formalism of Slater determinants (or antisymmetrized germinal powers), backflow, and Jastrow factors.For the 8 × 8 periodic system, at U = 4, VAFQMC is lower in energy across densities and, at U = 8, E/N from VAFQMC is approximately 0.004t lower at n = 0.78125 and 0.01t higher at n = 0.875.More impressively, at n = 0.78125 we obtain nearly the same energies for U = 4 and U = 8 (with U = 8 having worse performance) to a very sophisticated VMC wave function, enforcing symmetries such as SU(2) and k = 0 momentum projection, augmented by additional Lanczos steps.Likewise when comparing to fixed-node DMC  [17,18], fixed-node DMC [18], VMC [14], and VAFQMC (this work), on a PBC square Hubbard model with density n and interaction strength U (in units of t).Energies with a * are effectively exact from release constraint calculations of AFQMC [17], and energies with † include SU(2) and momentum projection augmented with a Lanczos step [14].Additionally shown is the variance σ 2 = ⟨H 2 ⟩ − ⟨H⟩ 2 of the energy, which can be computed more conveniently in variational methods and which vanishes (σ 2 = 0) for a true eigenstate of the Hamiltonian.VAFQMC results are optimized with N ℓ = 4 projector slices.
energies across system sizes and densities, we find similar or slightly lower energies for U = 4 and competitive energies for U = 8.This is particularly favorable as VAFQMC can conveniently compute both local and nonlocal observables and correlation functions, which can be challenging or inaccessible within fixed-node DMC and which require back-propagation [9] in AFQMC.We also note for the e.g. 10 × 10 U = 4 point, the local magnetization ⟨S z ⟩ is < 10 −2 on each site, recovering an unpolarized spin solution despite the initial parameters breaking spin symmetry explicitly.This result comes at a small cost however, and lower energy results can be obtained when starting from a UHF trial wave function (e.g.Fig. 1c).Using a UHF initial state |ϕ 0 ⟩ with U U HF = 0.5 (the same as used by AFQMC for its trial wave function), the local magnetization is ⟨S z i ⟩ ∼ O(10 −2 ) but about 0.05t lower in energy.This is not unique to VAFQMC but rather to variational wave functions generally, where local minima of symmetry broken states can be hard to optimize away.In VAFQMC this manifests as requiring more projector slices N ℓ to modify a strong local magnetization or density in the trial state, compared to a spin-restricted initial ansatz, and thus we suggest utilizing the restricted form to minimize total parameters.

B. Local Order Prediction and Automatic Symmetry Restoration
As a proxy for phase diagram exploration, we use VAFQMC to explore the local order of cylindrical Hubbard models with antiferromagnetic edge pinning.This model was used to study the stripe behavior of the ground state [28,38], In Ref. [28], the stability of the stripe state was determined by comparing the staggered spin and hole densities at different interaction strength U and electron density n as a function of system size.We perform calculations similar to those, but first with small enough system sizes that DMRG can provide quasi-exact reference results.We vary the electron density n, or equivalently hole doping of δ ≡ 1 − n, and consider cylinders of width r/δ × N y for some integer r.
With different interaction strengths and system sizes, we find remarkable agreement in the magnetic order between VAFQMC, DMRG, and AFQMC in the top of Fig. 3.The variational ansatz for VAFQMC is completely unconstrained by symmetries, yet can nearly recover both translational symmetry in the y direction and symmetry across the middle of the cylinder.The corresponding hole densities are shown in bottom of Fig. 3.We again see almost complete recovery of the symmetry of the system.There is small disagreement between DMRG and AFQMC as to the hole densities at large U .As DMRG is expected to be nearly exact in these width-4 systems, the discrepancy is likely an indication of small residual bias from the constraint in AFQMC [37].Remarkably our VAFQMC results are nearly indistinguishable from the DMRG results.This extends to other dopings across similar values of U as well (see appendix), and illustrates the predictive power of the VAFQMC method.At larger values of U , the variational energy from VAFQMC is 1-2% higher than DMRG and AFQMC.These results are a reminder that the total energy, while very important for a variational ansatz, should not be unduly emphasized, and order parameters and correlation functions often provide a more stringent measure of the predictive power of an ansatz.
We next move beyond width-4 cylinders and consider wider systems, which are important [28,34] to access the thermodynamic limit.In Fig. 4 a comparison of staggered S z and hole density between AFQMC and VAFQMC is shown for a 12 × 8 cylinder with n = 5/6 density and U = 2, 4, 6, 8. Exceptionally, VAFQMC continues to obtain the correct order and stripe wavelength, with reasonable quantitative agreement with AFQMC.Combined the results for larger periodic cells discussed in the previous section, these results show the promise of the VAFQMC ansatz as a predictive approach for ex- tended two-dimensional systems and beyond.
As mentioned, the projection times are unconstrained as full variational parameters in our approach.Unlike the ansatz in Ref. [15], a number of solutions have a negative τ l value within one of the projector slices, as shown in the left of Fig. 5, similar to that of higher order Trotter expansions [39,40].As both T σl ij and τ l can change, we renormalize the τ l values by the average value of nearestneighbor hopping terms ⟨T σl ⟨ij⟩ ⟩ which is almost always O(−t).By only considering τ l values, we can see that the effect of each time slice is not equivalent over a given U .As U increases, the value of τ l tends to fluctuate between positive and negative.We believe this is to effectively 'raise the temperature' and broaden the trial wave function in order to (in the next slice) project into the correct manifold.For width-4 systems, |τ l | values are around ≈ 0.13 while width-8 systems they are slightly larger at 0.15 on average, but both width cylinders have slices l where τ l σ ≈ 0.25 for example.Width-4 systems also tend to have larger negative τ l σ , with min l,σ τ l σ ∼ −0.2 for width-4 cylinders and min l,σ τ l σ ∼ −0.07 for width-8.We emphasize that while some values of |τ l | are quite small, the interaction term remains O(10 −1 ) and thus the ansatz is still not reduced to fewer slices.
To examine the role that the negative τ values play, in the right of Fig. 5 we plot the ⟨E l ⟩ for the optimized four slice VAFQMC ansatz truncated to l (l ≤ 4) time slices.(This is to be distinguished from a stand-alone l-slice or N ℓ = l wave function.In the latter the parameters would be optimized for l slices, while here they are optimized for the entire 4-slice wave function.)We subtract the energy of that of the |ϕ 0 ⟩ state (denoted RHF ) and rescale by 1/N .We see that for the first (l = 1) and last (l = 4) slice the energy always decreases.At intermediate values of l there may be energy values that are significantly worse than the previous time slice, unlike the near monotonic decrease in a QMC projector method.For example consider the 12 × 4 system at U = 4, the energy per site at an intermediate point l = 2 → 3 jumps by ∆E/N ∼ +0.17, then lowers in the final slice (l = 3 → 4) by ∆E/N ∼ −0.2.The large value of |τ 4 | at the end of the ansatz appears to be linked to this substantial decrease in energy, to compensate for the heating and cooling of the path.
For lower values of U , often a given slice may have a very small value of τ l , suggesting that an ansatz with fewer N ℓ may suffice to capture the dominate physics.This is further suggested by the intermediate values of the energy ⟨E l ⟩, which only occupies a range of ∼ 0.04/site.We conjecture that the flexibility in {τ l } can improve the size extensivity of the ansatz.The ansatz can both increase the overall effective β, similar to Ref. [15], to produce a size extensive ansatz, or may be something more complex within the path.The new ansatz based on an effective Hamiltonian given in Eq. ( 5) is the key for allowing better size consistency with a fixed N ℓ .More systematic studies on this will be very valuable.

IV. DISCUSSION AND SUMMARY
We presented a variational approach for many-electron systems with an expressive ansatz of O(N ℓ • N 2 ) parameters.Although inspired by the idea of an approximate realization of the imaginary-time projection in AFQMC, the ansatz is much more general.The approach shares the idea of variational optimization of a small number of projections with AFQMC-like time slices as in Ref. [15], and we have retained the name VAFQMC from there.However, our approach takes a different philosophy of viewing the ansatz as a variational effective Hamiltonian rather than a variation projection using the physical Hamiltonian.Instead of retaining the form of the Hamiltonian Ĥ in the projection (e.g., keeping T and the basic form of the HS transformation dictated by V ) and creating different orders via |ϕ 0 ⟩, we view the ansatz as seeking the most effective projection (to the true ground state) via variational optimization.
We have emphasized the underlying theory of optimizing a set of effective one-body Hamiltonians or actions defined by auxiliary fields.This wave function, even in our still rudimentary implementation, is already competitive with state of the art VMC and even fixed-node Green's function Monte Carlo.We have demonstrated this on large square Hubbard models, where the approach can systematically obtain the ground state energy to within or better than 1-2% of the exact results.Furthermore, it nearly restores symmetry automatically via its optimization process and provides accurate physical observables which allow predictive resolution of magnetic and charge correlations and orders.Broken e.g.SU (2) symmetry from pinning fields can be additionally restored by hand and used in the original model.
Historically, many VMC wave function ansatze were built around a specific mean-field state (for example in Ref. [15]), or an explicit form that targets a particular ⟨ij⟩ ⟩.The electron density is n = 5/6.Right: Energy per site (relative to that from |ψRHF ⟩ = |ϕ0⟩) at a given time slice l (in units of t) for the systems shown on the left.Note that the energy at each l is evaluated with the final parameters fixed at the optimized values for l = 4. Marker size is proportional to the renormalized τ l = max(τ l ↑ , τ l ↓ ) with red denoting positive (equivalent to cooling projection) and purple denoting negative (equivalent to heating projection).Error bars are smaller than the marker size.type of order.Under this approach, optimization is constrained to the fixed state type, and prediction of the actual order is made by comparing the relative energetic behaviors of two or more such fixed types of states.Despite its considerable success, this approach lacks true predictive power in many situations, either because the order is not known, or because the different types being compared are not well balanced in the ansatz (e.g., one with more variational flexibility than the other).With the onset of neural quantum states (e.g.Refs.[11,12,41]), wave functions are parameterized with partially or fully unconstrained variational freedom to improve expressibility.Likewise, the aim of our VAFQMC approach is to have an ansatz which is similarly more flexible and more expressive such that local order can instead be predicted from optimization.As our examples have shown, with large optimization freedom VAFQMC can function as a predictive ansatz in addition to providing good variational energies in a challenging system operating at realistic parameters and under "real-life" conditions.It is also notable that this is achieved without the use of neural networks or tensor networks, bringing a new class of unconstrained systemically improvable variational wave functions.
Optimization of an expressive ansatz alone is some-times still not enough to resolve the hardest problems.For models in which there are many competitive states, e.g.stripes or intertwined orders, the application of pinning fields has become a very useful and essential tool to elucidate the underlying order.The success with which this paradigm has been used [28,34,35,38,42], and this benchmark study of VAFQMC suggests that this is a potential paradigm to better solve problems inaccessible to 1D tensor networks.This should prove to be an invaluable method in the toolbox of computational many-body physics.
Despite the 'inexpensive' ansatz presented here with fixed N ℓ , we emphasize that VAFQMC can still be systematically improved.One can continue to add more projector slices to potentially improve the energy, but after some point computational costs will increase (and further care must be taken to properly stabilize the 'long time' projection of the trial state).At N ℓ = 4 this was functionally unnecessary for the systems we considered.Neural networks can be added within the ansatz, for example as in Ref. [16], to further improve expressibility, particularly in the strongly correlated regime where correlations between auxiliary fields may be useful.Additional considerations may be needed to improve scalability, however.Likewise while the optimization proce-dure was able to regain translation invariance and minimal local spin symmetry breaking, this should not be relied upon when scaling.Applying symmetry projection [43,44], either within |ϕ 0 ⟩ or for each projected |x⟩ state, is a compelling future direction.
Given the promise the ansatz has already shown at this very early stage of development, we believe the method presents many exciting and interesting new possibilities.In addition to the directions discussed above, there are a number of issues to be addressed or better understood.For example, optimizing to machine precision is currently intractable due to an infinite variance problem.Because the local energies do not have a zero variance principle, sampling can be dominated by rare events and thus are unable to reduce the statistical noise with more samples.We have checked our results on small systems using the methods of Refs.[45] (see appendix), however there may be a better methods to be explored in future work.lower in energy as U gets larger while VAFQMC results trend higher.The largest discrepancies occur at U = 8 where VAFQMC is approximately 1 − 1.5% from DMRG while AFQMC is < 1% lower (relatively).
We show further local observables comparisons, of av-erage staggered S z and hole density, at n = 5/6 for 12×4 and 24 × 4 system sizes in Fig. A7.The VAFQMC results of spin and hole densities mirror those at in Fig. 3 are seen to be in excellent agreement with those from DMRG.While staggered spin densities agree well be-

FIG. 1
FIG.1.a) Cartoon of standard AFQMC where random walks are taken with many small steps of size τ , following a pre-defined projector, to reach an imaginary time of β. b) cartoon of VAFQMC random walk, where a number of potentially large steps are taken, following an effective projector that is variationally optimized, to reach the equivalent imaginary time of β.Each dot represents a projector slice and is a function of optimizable parameters α. c) Staggered magnetization measurement for the 3rd column of a 10 × 10 U = 4 Hubbard model with n = 0.8 density.The first pane shows the initial unrestricted Hartree-Fock (UHF) solution, and subsequent panes show VAFQMC slowly projecting out the incorrect initial order to create a spatially uniform ground state.d) Staggered magnetization for a 20 × 4 Hubbard model at U = 4 and n = 0.9 density, with magnetic pinning fields applied at the edges.Starting from a qualitatively incorrect RHF initial state and optimizing the projection can uncover the correct ground state local order.

FIG. 2 .
FIG. 2. Variational energy ⟨ Ĥ⟩ and variance per site (⟨ Ĥ2 ⟩ − ⟨ Ĥ⟩ 2 )/N (inset) during optimization of a 4×4 Hubbard model at U = 6 and n = 0.625 using VAFQMC with N ℓ = 4.The ground state from exact diagonalization is shown in a black dashed line.Fluctuations are due to the limited sampling during optimization.

FIG. 5 .
FIG.5.Left: Optimized values for τ l maximized over spin (τ l = max(τ l ↑ , τ l ↓ )) for the Hubbard model in cylindrical geometry with edge pinning, shown for four values of U and three system sizes.Values have been renormalized by ⟨T σl ⟨ij⟩ ⟩.The electron density is n = 5/6.Right: Energy per site (relative to that from |ψRHF ⟩ = |ϕ0⟩) at a given time slice l (in units of t) for the systems shown on the left.Note that the energy at each l is evaluated with the final parameters fixed at the optimized values for l = 4. Marker size is proportional to the renormalized τ l = max(τ l ↑ , τ l ↓ ) with red denoting positive (equivalent to cooling projection) and purple denoting negative (equivalent to heating projection).Error bars are smaller than the marker size.