Table Of Content

Symmetry-projected Wave Functions in Quantum Monte Carlo Calculations Hao Shi1, Carlos A. Jiménez-Hoyos2, R. Rodr´ıguez-Guzmán2,3, Gustavo E. Scuseria2,3 and Shiwei Zhang1 1 Department of Physics, The College of William and Mary, Williamsburg, Virginia 23187 2 Department of Chemistry, Rice University, Houston, Texas 77005, USA 3 Department of Physics and Astronomy, Rice University, Houston, Texas 77005, USA We consider symmetry-projected Hartree–Fock trial wave functions in constrained-path Monte Carlo (CPMC) calculations. Previous CPMC calculations have mostly employed Hartree–Fock (HF)trialwavefunctions,restricted orunrestricted. Thesymmetry-projectedHFapproach results in a hierarchy of wave functions with increasing quality: the more symmetries that are broken 4 and restored in a self-consistent manner, the higher the quality of the trial wave function. This 1 hierarchyis approximately maintained in CPMC calculations: theaccuracy in theenergy increases 0 andthestatistical variancedecreases whenfurthersymmetriesarebrokenandrestored. Significant 2 improvement is achieved in CPMC with the best symmetry-projected trial wave functions over those from simple HF. We analyze and quantify the behavior using the two-dimensional repulsive n Hubbardmodelasanexample. Inthesign-problem-freeregion,whereCPMCcanbemadeexactbut a J aconstraintisdeliberately imposedhere,spin-projectedwavefunctionsremovetheconstraintbias. Away from half-filling, spatial symmetry restoration in addition to that of the spin lead to highly 1 accurate results from CPMC. Since the computational cost of symmetry-projected HF trial wave 3 functionsin CPMC can bemade toscale algebraically with system size, thisprovidesa potentially ] general approach for accurate calculations in many-fermion systems. l e PACSnumbers: 71.10.Fd,02.70.Ss,05.30.Fk,71.27.+a - r t s . I. INTRODUCTION where the local interactions lead to auxiliary-fields that t a are real, and by a phase problem in realistic electronic m systems, where the Coloumb interaction leads to com- Approximations able to account for the most relevant - plex fields.14 Inorderto removethe exponentialincrease d correlations in low-dimensional many-fermion systems of the statistical noise with system size (or inverse tem- n represent a cornerstone in condensed matter physics. In perature), we control the sign/phase problem approxi- o particular,verychallengingphenomenasuchasthehigh- mately. The quantum Monte Carlo (QMC) process is c Tc superconductivity,1 the colossal magnetic resistance2 [ formulated as open-ended random walks in the space of as well as the superconductivity in the iron-based com- Slater determinants. A constraint, formally exact, is im- 1 pounds,3,4 just to mention three of the most notable ex- posed on the paths of the random walk according to a v amples,requireustobetterunderstandthenatureofthe 8 electron-electroncorrelationsandtheirimpactonthere- signorphaseconditionofthe walker. This approachhas 1 beenreferredtoastheconstrainedpathMonteCarlo15,16 sulting properties of the considered quantum systems. 0 (CPMC)methodinthecaseof“only”asignproblemand Withinthis context,the Hubbardmodelis regardedasa 0 thephaselessAFQMC17forthegeneralcasewithaphase paradigm since, in spite of its essential simplicity, it dis- . 2 problem. In practice, the constraint is implemented by plays several challenging properties associated with the 0 the overlap of the sampled Slater determinants with a relevantphysicsinstronglycorrelatedmany-electronsys- 4 trialwavefunction. Despite its approximatenature, this tems. From the theoretical point of view, Hubbard-like 1 approach is very accurate even with a simple mean-field : models also represent valuable benchmark systems for v trialwavefunctionas ithas beenshownbefore ina wide testingdifferentstate-of-the-artapproximationsthatcan i range of applications.14,18–20 X subsequently be extended to more realistic situations. r Nowadays there are several such approximations al- Inquantummechanicalsystems,thesymmetriesofthe a ready used to study both the one (1D) and two- Hamiltonianarecrucialincharacterizingpropertiessuch dimensional (2D) Hubbard models with variable degrees as the excitation spectrum, the basic fingerprint of the of success. For sufficiently small lattices, exact diago- system.21–27 Wehavestudied28 howtheuseofsymmetry nalization is possible. However, due to its exponential influences the CPMC method and AFQMC calculations cost, such an exact diagonalization becomes impossible ingeneral,bothintheformoftheHubbard-Stratonovich beyondagivensystemsize. Onecanthenresorttoother (HS) transformation and in the trial wave function. It approximations such as the variational Monte Carlo5–7 was shown that imposition of symmetry properties re- (VMC),theCoupledCluster8 (CC)andtheDensityMa- sulted in a significant increase in the accuracy and effi- trix Renormalization Group9–11 (DMRG) methods. ciency in the QMC. Choosing a symmetry-adapted HS Auxiliary-field quantum Monte Carlo12–14 (AFQMC) transformation often does not change the cost of the is one of the most popular methods to extract collective calculation significantly. On the other hand, symmetry- properties of quantum many-body systems. However, adapting a trial wave function may result in expansions it is limited by a sign problem in Hubbard-like models, whose length depend on system size. In our prior work, 2 we obtained symmetry-adapted wave functions by using II. METHOD acompleteactivespace(CAS)expansioninthesubspace of the open shell. The problem with CAS expansions is In this work, we use as benchmark system the two- that the size of the wave function scales exponentially dimensional repulsive Hubbard model, which is written with the dimension of the active space. in second-quantized form as33 An alternative to obtain highly correlated wave func- L L tions which, at the same time, respect the original sym- Hˆ =Kˆ +Vˆ =−t c†iσcjσ+U ni↑ni↓, (1) metries of the considered many-fermion problem is pro- hXi,jiσ Xi videdbysymmetryrestorationviaprojectiontechniques. Very recently, a hierarchy of variational symmetry- where the first term represents kinetic energy from elec- projected approximations, inherited from nuclear struc- tron hopping (t > 0) and the second is the repulsive ture physics,21 has been considered to describe the on-site interaction (U > 0). The operators c† and c iσ iσ ground and excited states of the 1D (Ref. 29) and 2D create and annihilate an electron with spin direction σ (Ref. 27)Hubbardmodels. Theyhavealsobeensuccess- on the i-th lattice site. Note that we have introduced a fully incorporated into the description of molecular sys- vector notation for the site index, i.e., i = (i ,i ). The x y tems within the Projected Quasiparticle Theory (PQT) operators nˆ = cˆ† cˆ are the local number operators. framework.30,31 A symmetry-projected HF wave func- iσ iσ iσ The notation i,j is used to denote that only hopping tion considers a Slater determinant φ¯ which deliber- h i between nearest-neighbor sites is allowed. We assume | i ately breaks several symmetries of the original Hamilto- periodic boundary conditions (PBC) in both the x and nian. By the superposition of the (degenerate) Gold- y directions. Throughout this paper, energies will be re- stone manifold Rˆ φ¯ , where Rˆ represents a symmetry ported in units of t and we set t=1. | i operation, one can recover the desired symmetry in a fullyself-consistentvariation-after-projection(VAP)procedure. In addition, the intrinsic Slater determinant φ¯ | i A. Symmetry-Projected Wave function resulting from the symmetry-projected VAP procedure containsdefectsthatprovideanintuitivephysicalpicture regarding the basic units of quantum fluctuations in the In this section, we briefly discuss the form of the considered systems.29 Recently symmetry-restored trial symmetry-projected HF states used as trial wave func- wave functions have been used in CPMC in the study of tions. For more details of the symmetry-projected ap- nuclear shell models.32 proach, we refer the reader to our previous works (Refs. 29 and 27). WeworkwithHF-typeSlaterdeterminantswhichpre- It is known that symmetry-projection out of a Slater servethenumberofelectronsinthesystem. Anarbitrary determinantresultsinawavefunctionthatisnotsizeex- tensive,31 i.e., the correlation energy recovered by such HF-type Slater determinant can in generalbreak several symmetries of the original Hamiltonian. Typical exam- anansatzdoesnotscalelinearlywithsystemsize. Thisis ples are the rotational(in spin space) and spatial (space also true for the CAS expansions we have consideredbe- group for periodic systems) symmetries. In the present fore. Symmetry-projected wave functions can, neverthe- study, we have considered two different kinds of sym- less, account for some of the strong correlations (due to metry broken states: (unrestricted) UHF configurations the degeneraciesin the single-particlespectrum). In this which are eigenstates of S and therefore have a definite paper, we examine the use of symmetry-projected wave z N and N , as well as (generalized) GHF configurations functions as trial states within the CPMC framework. ↑ ↓ thatbreakS andthereforecanonlybe characterizedby We stress that, even though a single symmetry-broken z an overall N number of electrons. Slaterdeterminantisusedinthe trialwavefunction, the To restorethe spin quantumnumbers,we use the spin symmetry-projected wave functions are, via the projec- projection operator tionoperators,multi-determinantinnature. Inaddition, the symmetry-projected state can access highly-excited 2S+1 configurations. In this way, they bring a very sophisti- PˆS ′ = dΩ S∗′(Ω)Rˆ(Ω), (2) cated trial ansatz which can be expected to improve the ΣΣ 8π2 Z DΣΣ quality of the sign/phase constraint in CPMC. Indeed, ourresultsshowthattheuseofsuchtrialwavefunctions where Rˆ(Ω)=e−iαSˆze−iβSˆye−iγSˆz is the rotation opera- greatly improves the CPMC results. The combined ap- torinspinspace,thelabelΩ=(α,β,γ)standsfortheset proach demonstrates dramatically better behavior that of Euler angles, and S ′(Ω) are Wigner functions. We DΣΣ successfully removes the size-extensivity problem of pro- notethat,ifUHFwavefunctionsareused,thenumerical jected HF. Often the full symmetry-projected trial wave effort can be alleviated significantly in the evaluation of function eliminates the systematic error from the con- matrix elements, as both the integrals over α and γ can straint, and the CPMC calculations approach the exact be carried out analytically. answer, with a reduced variance. To recover the spatial symmetries we use the generic 3 projection operator • SG (space group) is used for wave functions in which the full space group of the 2D square lattice Pˆk ′ = l Γk∗ ′(g)Rˆ(g), (3) has been broken and restored. mm h mm Xg • S is used to indicate that S projection has been z z done (out of a GHF-type determinant). where Γk is an irreducible representation (irrep), which can be found by standard methods, and Rˆ(g) represents • S (spin) is used to indicate that full spin projec- thecorrespondingsymmetryoperationsintheconsidered tion has been carried out. If the determinant is of latticesparametrizedintermsofthelabelg. Inaddition, GHF-type(noncollinear),itwillinvolvethetriaxial l is the dimension of the irrep and h is the order of the integration of Eq.(2). group in Eq.(3). We note that, for the periodic Hub- bard 2D system, the irreducible representations associ- For example, SG,S-UHF means that the trial wave ated with the spatial symmetry can be characterized by functionwaspreparedbybreakingandrestoringspinand the linear momentum (k) inside the Brillouin zone and, spacegroupfromaUHF-typedeterminant,whileLM,Sz- for certain high-symmetry momenta, by additional pari- GHF indicates that linear momentum and Sz have been ties under x, y, and x-y reflections. In what follows, we broken and restored out of a GHF-type determinant. In willexplicitlyprovidethelinearmomentumoftherecov- bothcases,the determinantis variationallyoptimizedto eredirrepand,whereappropriate,wewillfurtherprovide minimize the energy in the presence of the symmetry the paritiesofthe recoveredstate(if the fullspacegroup projection operators, and before introducing the Monte projection is carried out). Carlo procedure described in the next section. We then superpose the Goldstone manifolds due to spin and/or spatial symmetries via the following ansatz B. CPMC |ΨΘKi= fKΘ′PˆKΘK′|φ¯i, (4) We briefly summarize the ground-state CPMC XK′ method15,16 below to facilitate the ensuing discussions. ThereaderisreferredtoRef.14andreferencesthereinfor wherefΘ′ arevariationalparameters,and φ¯ istherefer- further details. Aside from the symmetry-projectedtrial K | i enceHF-typesingleSlaterdeterminant. Here,theindices wave function, there is a new algorithmic element in the Θ and K denote quantum numbers associated with spin present study, namely we have generalized the AFQMC and/or spatial symmetries. Let us stress that, through approach to have GHF-type of random walkers, as dis- the actionofthe projectionoperatorPˆΘ , the multide- cussed in Sec. IIIB. This will allow the application of ′ KK terminantal character of the state characterized by the CPMC and AFQMC in generalto Hamiltonians that do quantum numbers Θ and K′ is recoveredand written in notconserveSˆ ,forexampleonethatincludesspin-orbit z terms of the quantum numbers Θ and K. The linear coupling. combination introduced in Eq.(4) guarantees the inde- All ground-state AFQMC methods are based on the pendence of ΨΘ with respect to rotations of the Slater projection determinant|φ¯K.i The wave |fuinctions ΨΘ of Eq.(4) are precisely the Ψ lim e−β(Hˆ−ET) Ψ , (5) ones used in the presen|t sKtuidy as trial states within the | 0i∝β→∞ | Ti CPMC scheme. The wave functions, Eq.(4), are deter- where E is a trial ground-state energy (an initial guess mined byresortingto the Ritz-variationof the projected T which is then improved in the calculation) and Ψ is energywithrespecttoboththe mixingcoefficients(fKΘ′) a guess of the ground state wave function, a tria|l wTiave and the HF-type determinant φ¯ . We stress that the function,whichistypicallytakentobeasingleSlaterde- | i wave functions are fully optimized in the presence of the terminantoralinearcombinationofSlaterdeterminants. projection operations, i.e., a VAP approach. In our pre- Ψ Ψ = 0 needs to be satisfied in order to project to 0 T vious work,29 we have discussed how the optimized de- h | i 6 the ground state. In a numerical method, the limit can terminant φ¯ developsdefects thatcanbe interpretedas be obtained iteratively by | i the basic units of quantum fluctuations in the system. We refer the reader to Ref. 27 for a discussion of how Ψ(n+1) =e−∆τHˆ Ψ(n) , (6) the optimization procedure is carried out in practice. | i | i Before concluding this section, let us introduce the where Ψ(0) = Ψ and ∆τ is a small positive number. T notation used in this paper for the symmetry-projected The| Troitter-|Suziuki breakup34,35 and a so-called wave functions of Eq.(4): Hubbard-Stratonovich (HS) transformation are used to break the propagatorinto many one-body propagators14 • LM (linear momentum) is used to denote wave functions in which the linear momentum has been e−∆τHˆ = dxp(x)Bˆ(x), (7) broken and restored. Z 4 wherep(x)isaprobabilitydensityfunction,forexample, with the ground state becomes zero14 auniformdistributionofoneIsingfieldpersite(x = 1 i with i = 1,2, ,L.). The one-body propagator Bˆ(±x) Ψ0 φ =0, (9) h | i ··· in Eq. (7) has a special form, namely, a product of the since they will contribute zero at any future projections. exponential of one-body operators Inpractice,inplaceofthe exactgroundstatewavefunc- Bˆ(x)=e−∆τKˆ/2evˆ(x)e−∆τKˆ/2 (8) tion, a trial wave function ψT is used instead to imple- | i ment the constraint above, wherevˆ(x)isaone-bodyoperatorwhosematrixelements are simple functions of the HS fields and of magnitude ΨT φ =0, (10) h | i (√∆τ). O whichwillintroduceasystematicbiasintheresults. The By applying each one-body propagatorto a Slater de- importance sampling transformationimposes this condi- terminantwavefunction,wewillgenerateanotherSlater determinant.36 Themany-dimensionalintegral/sumover tion naturally, terminating randomwalk paths that lead theauxiliaryIsing-fields(foreachcomponentofxandat to a zero-overlap with ΨT . As mentioned, the system- | i atic error of the constraint tends to be very small even each imaginary-time iteration n) is performed by Monte withsimplemean-fieldtrialwavefunctions. Wewillshow Carlo sampling the fields. The resulting linear combi- below that, with the symmetry-projectedwave function, nation of Slater determinants at each iteration n gives a stochastic representation of the wave function Ψ(n) . theconstrained-pathbiasisfurtherreducedandbecomes | i negligible in most of the systems studied. After convergence, all the sampled determinants from n n can be used collectively to represent the ground eq ≥ state wave function. The determinants have to be sta- III. RESULTS bilized periodically with, for example, a modified Gram- Schmidt orthogonalizationprocedure.37 CPMCusesimportancesamplingtosteerthesampling A. Half-filling toward more “likely” auxiliary-fieldconfigurations. This is achieved by a similarity transformation with a trial We first discuss the use of symmetry-projected trial wave function ψ φ to modify the probability density wave functions for the half-filled Hubbard model. Since T h | i function p(x) in Eq. (7). The importance sampling does thereis nosignprobleminthis case,CPMCcalculations notaffecttheaveragevaluesofthecomputedobservables, canbe made exactby redefining the importantsampling only the variance. The better the trial wave function to have a nonzero minimum.15,16 If we ignore this and ψ , the smaller the statistical error for a fixed amount apply the usualimportancesampling with Ψ literally, T T | i | i of Monte Carlo samples. Since the CPMC process is a the random walks can be constrained to a part of Slater branchingrandomwalk,populationcontrolhastobeap- determinant space because the paths are terminated by pliedperiodically.15,38 Asamanifestationofthevariance theconditioninEq.(10). TheCPMCresultswilldisplay reduction, the extent of weight fluctuation (branching) a constraint bias. In this section, we will use the half- willbereducedwithabetter ψ . Theuseofsymmetry- filled case in the artificial way with the constraint as a T | i projectedtrialwavefunctions leadsto areductionofthe benchmark system to study the effect of the symmetry statistical variance, as will be illustrated below. projectedtrialwavefunction. Becausetheexactresultis The sign problem14 is the leading difficulty in QMC accessible,thisstudyallowsustosystematicallyexamine simulations of many-fermion systems. The problem oc- the effect of the different symmetries in Ψ . T | i cursbecausetheprojectionisingeneralsymmetricabout Figure1isanillustrationinthe4 4Hubbardmodelat × φ and φ . In the random walks, there is no mecha- half-filling (N =N =8) with U =4. Here, the ground ↑ ↓ | i −| i nism to distinguish random walkers of opposite overall statecorrespondstoaspinsingletwithmomentum(0,0) signs. In the course of the projection, if we switched the andevenparity under allreflections. As we see fromthe sign of each random walker (e.g.,by exchanging two or- bottompanel,CPMCwithaUHFtrialwavefunctionhas bitalsin φ ,or φ , inanUHF typesimulation,ortwo abiasintheground-stateenergy,withthecalculateden- ↑ ↓ | i | i spin-orbitals in a GHF-type simulation), there would be ergy higher than the exact result ( 13.62185) by about − no noticeable change and the sampled population would 1%. Asmoresymmetriesareincludedintheprojected ∼ giveanoverlapofoppositesignwithanytrialwavefunc- UHFbasedwavefunction,thebiasisseentodecrease,as tion. Inthehalf-fillingcaseoftheHamiltonianinEq.(1), well as the corresponding Monte Carlo variance, shown itturnsoutthatsymmetryinthepropagatorandin Ψ in the top panel. The only exception is in the variance T | i can make the sign of the overlap remain non-negative, from the UHF trial wave function, which falls slightly which eliminates the sign problem. In general, however, smallerthanthose fromLMUHF andSGUHF trialwave the overlap with any trial wave function (or with Ψ ) functions. The reduction of energy bias and energy vari- 0 | i will be zero when averagedover imaginary time (or over ance indicates that the quality of the trial wave func- the entire population at any large n), leading to infinite tion is improving. This is born out by their variational Monte Carlo variance. The sign problem can be con- energies, which decrease monotonically following the se- trolled,exactly,byeliminatingthewalkerswhoseoverlap quence from left to right (not shown). 5 WealsoseefromFig.1thattheCPMCenergybecomes indistinguishable with the exact one as long as the spin symmetry is imposed. This is consistent with the find- ings from our previous studies of symmetry-preserving trial wave functions which were obtained from a CAS approach by diagonalizing a truncated active space.28 ((aa)) -4 Statistical Variance -5 -6 -7 ar) -8 V -9 g( -10 Lo -11 -12 -13 -14 -15 UHF LMUHF SGUHF S-UHF LM,S-UHF SG,S-UHF (b) 1.2 Relative Error % 1 e nc 0.8 FIG.2: (Coloronline)Size-dependenceofthecomputedener- e er 0.6 giesinthehalf-filledHubbardmodelatU =4. Thebehaviors Diff 0.4 of two different forms of trial wave functions are compared. y g 0.2 Thevariationalenergiesofthetrialwavefunctionsareshown er En 0 in the upper curves. The dashed lines are polynomial fits to -0.2 thedata. CPMCresultsusingthetwotrialwavefunctionsas UHF LMUHF SGUHF S-UHF LM,S-UHF SG,S-UHF (artificial)constraintsareshowninthelowercurves,together withexact free-projectionresults. (Notethedifferentvertical FIG. 1: (Color online) Performance of different symmetry- scales between the upper and lower curves.) Our best esti- projected trial wave functions in CPMC. Panel (a) plots the mate of the ground-state energy in the thermodynamic limit statistical uncertainty of the computed ground-state energy is given by the horizontal lines. Its statistical uncertainty is in a semi-log plot (arbitrary scale). Panel (b) plots the per- indicated bythe width of theline in theinset. centage error, (ECP−E0)/|E0|, between the CPMC energy andtheexactenergy. Thesystemisthe4×4Hubbardmodel, with N↑ = N↓ = 8 and U = 4. The number of walkers was up to 14 14 lattices with twist-average boundary con- 1000, an equilibrium phase of βeq = 4 was discarded, and × ∆τ =0.01 was used in the runs. ditions, which gives much smaller finite-size effects than PBC. A polynomial fit39 with a leading term of 1/L3/2 We next study the behavior of the simulation as a was performed on the data (not shown), and our best function of system size. The results are summarized in estimate of the energy per site for the half-filled Hub- Fig. 2, which shows the calculated energy per site from bard model at U = 4 as L is E0/L = 0.8600(1). → ∞ − various approaches versus the inverse linear dimension The CPMC/UHF results are seen to convergeto a value of the square lattice. Two types of trial wave functions slightly higher than the exact energy in the thermody- arecompared,thestandardUHFandthespin-symmetry namic limit. (onto a singlet) projected S-UHF. As mentioned ear- CPMC using the symmetry-projected S-UHF trial lier,symmetry-projectedHFwavefunctionsarenotsize- wave function is seen to agree with exact free projec- extensive. The variational energy (per site) of S-UHF tion at all finite lattice sizes studied. This extends the is considerably lower than that of UHF for smaller lat- conclusion from Fig. 1 that trial wave functions with tice sizes, but gradually approaches the UHF value as spin symmetry have negligible bias to much larger su- the system size is increased. As the fits indicate, the percells. These data indicate that, despite the lack of two give the same result in the thermodynamic limit, size-extensivityintheS-UHFtrialwavefunction,CPMC E /L 0.797. seems to restore a consistent behavior as a function of UHF ∼− The CPMC results with an artificial constraint using system size. UHFareshownasCPMC/UHF.Exactresultsfromfree- projection (i.e., by “lifting” the constraint)28 shown as FPMC/UHF are given for systems up to 16 16. The B. Away from half-filling × CPMC/UHFbiasisvisibleatalllatticesizes. Thelargest discrepancy is seen between the two at the smallest lat- We next move away from half-filling to examine the ticesize(rightmostpoint),correspondingtoa1%errorin effect of the symmetry-projected HF trial wave function Fig.1. Convergenceofthe ground-stateenergywithsys- on the sign problem. We study the system close to half- tem size is not monotonic under periodic boundary con- filling,theregionwhichhasthemostseveresignproblem. ditions,asseeninthebottomsetofcurvesandtheinset. IntermediatevaluesoftheinteractionstrengthU arecho- Toaccuratelydeterminetheexactground-stateenergyin sen, comparable to the band width, where the model is the thermodynamic limit, we have done calculations on potentiallymostrelevanttostronglycorrelatedmaterials 6 such as high-T superconductors. tems haveerrorsof 0.02%andthe maximumerroris c ∼− 0.07% at U = 8. Note, that the ground-state energy in ((aa)) 0.8 CPMC is calculated with the so-called mixed estimate, error 0.7 444xxx444---676uuu676ddd---448UUU which is not variational.40 e 0.6 4x4-7u7d-8U v ati 0.5 nal rel 00..34 0.0015 SSGG,,SS--UGHHFF o ati 0.2 ari 0.1 0.001 v 0 UHF S-UHF LM,S-UHF SG,S-UHF or 0.0005 (b) 0.035 err e ve error 0 00.0..002235 relativ 0 ati 0.015 -0.0005 el c r 0.01 qm 0.005 -0.001 p 0 c -0.005 -0.0015 UHF S-UHF LM,S-UHF SG,S-UHF U=4 4x4 6u 6d U=8 U=4 4x4 7u 7d U=8 FIG. 3: (Color online) Variational energies (upper panel) FIG. 4: (Color online) Illustration of GHF-based symmetry- of the UHF and three UHF-based symmetry-projected wave projected trial wave functions in CPMC. The relative errors functions, and the corresponding CPMC energies (lower in thecomputed correlation energy from CPMC using GHF- panel) using these as trial wave functions. The energies are typewavefunctions(SG,S-GHF)arecompared with thecor- shown in terms of relative errors in the correlation energy respondingUHF-type(SG,S-UHF).Thelatteraretakenfrom (see text). Note the different scales between the two pan- the rightmost group in Fig. 3. The same computational pa- els. Foursystemsarestudied,eachrepresentedbyadifferent rameters are used for the two typesof calculations. color/pattern of vertical bar. The CPMC results were obtained with a population size of 1,000 and 10 independent blocks of size β ∼1 each. The statistical error bars are indi- We have also investigated the effect of GHF-based cated. Trotter error is negligible. symmetry-projectedwavefunctions. IntheusualCPMC calculations discussed thus far, the randomwalkershave Figure 3 shows the results in 4 4 lattices where ex- a form φ↑ φ↓ , where the - and -spindeterminants act diagonalizations are available×for comparison. Two haveN↑| anid⊗N|↓ oirbitals,each↑specifie↓dbyLamplitudes, filling values, N = N = 6 and 7 are each studied at as in UHF. The one-body propagators in Eq. (8) decou- ↑ ↓ U = 4 and 8. The ground state in all cases is a spin ple into - and -spin components which operate on φ↑ ↑ ↓ | i singlet with momentum (0,0) transforming as the B1 ir- and |φ↓i, respectively. In the GHF-type, a walker |φi rep of the C4v group, save for the 12-electron system at contains (N↑ +N↓) spin-orbitals, each of which is spec- U =8, where the ground state switches to the A irrep. ified by 2L amplitudes that evolve stochastically. The 1 ThevariationalenergiesfortheUHFwavefunctions and one-body propagators are given by 2N 2N matrices. × three symmetry-projected HF wave functions are shown In the present study, the form of the initial/trial wave in the top panel, while the corresponding CPMC ener- function dictates which of the two approaches is used. gies are shown in the bottom panel. Here the energies Since there is more freedom in wave functions of are shown in terms of relative errors in the correlation symmetry-projected GHF form,29 lower variational en- energy,definedas(E E )/(E E ),whereE isthe ergycanbe obtained. Note thatlowerenergiesmayonly 0 RHF 0 − − calculatedenergy(eithervariationalorCPMC),E isthe be obtainedwhen the symmetry-projectedstates areop- 0 exact result, and E is the restricted HF (equivalent timized in the presence of the projection operators (a RHF to the Fermi gas wave function) energy. The variational VAP approach). This is seen from Fig. 4, which com- energy improves as more symmetries are broken and re- pares CPMC results from GHF-based and UHF-based stored in the trial wave function, as expected. The cor- trialwavefunctionswiththesamesymmetry-projections. respondingCPMCusingUHFtrialwavefunctionsyields A significant further reduction is seen in the CPMC re- relative errors in the computed correlation energy of a sults with GHF-based trial wave functions, leading to few percent,while CPMCwithsymmetry-projectedtrial essentially exact energies. In Table I, we summarize all wavefunctionsingeneralhassignificantlysmallererrors. theCPMCground-stateenergieswithbothtypesoftrial CPMC/SG,S-UHF is the most accurate for the four sys- wave functions. tems studied here, with relative errors of about 0.1% in The improvements in the energy from CPMC us- the correlation energy. In terms of percentage errors of ing symmetry-projected HF wave functions as the con- the total energy (as plotted in Fig. 1), the U = 4 sys- straint show that these wave functions give better de- 7 TABLE I: CPMC ground-state energies with UHF and various symmetry-projected trial wave functions. Two fillings (first column) and U values are studied on a 4×4 lattice Hubbard model. The exact results in the last column are from exact diagonalization. The statistical error bars in theCPMC results are on thelast digit and shown in parentheses. System U UHF S-UHF LM,S-UHF SG,S-UHF SG,Sz-GHF SG,S-GHF exact 6↑ 6↓ 4 -17.703(6) -17.727(3) -17.7428(9) -17.7327(8) -17.7316(4) -17.7301(1) -17.7296 6↑ 6↓ 8 -14.73(4) -14.63(1) -14.784(3) -14.914(3) -14.907(2) -14.920(1) -14.925 7↑ 7↓ 4 -15.688(4) -15.758(3) -15.753(2) -15.7482(5) -15.7500(8) -15.7455(2) -15.7446 7↑ 7↓ 8 -11.77(1) -11.628(9) -11.844(7) -11.872(3) -11.847(3) -11.8689(9) -11.8688 scriptions of the ground-state properties than standard suchtrialwavefunctionsareused. WeseethattheSUHF mean-field. One can then expect that such wave func- trial wave function leads to little improvement, despite tions will also improve the calculations of other ob- havingimprovedtheCPMCenergy. Asmoresymmetries servables. We study this in Fig. 5 with the exam- are included in the projected HF wave functions, the re- ple of momentum distribution. In CPMC calculations sults improve. The largeststepoccurswhenspace-group of the expectation of observables which do not com- is added. This seems reasonablegiventhe open-shellna- mute with the Hamiltonian, the back-propagation tech- ture of this system. The result with the full SG,S-GHF nique15,41 is almost always used to obtain a “pure” es- trial wave function is essentially indistinguishable from timator Ψ e−βBPHOˆ Ψ / Ψ e−βBPH Ψ ,asopposed exact diagonalization. T 0 T 0 h | | i h | | i to the mixed estimate (correspondingto βBP =0)which In Table II several CPMC results are shown using is exact for the energy but biased for a general operator symmetry-projected HF trial wave functions for larger Oˆ. For reasonable choices of the back-propagation time systems away from half-filling. Accurate energies from βBP, the technique gives exact estimators of any Oˆ ex- constraint release calculations28 with CASSCF wave h i cept for constrained-path errors. Here, for illustrative function are available in these systems, which we use for purposes, we use the mixed-estimate to calculate n(k), benchmark and are shown in the last column. The UHF which helps to magnify the effect of ΨT . trial wave functions are very good at these fillings and h | alloftendifficult tosurpass,28 soitis significantthatthe 1 ED symmetry-projectedtrialwavefunctionsexceptthemin- 0.9 UHF imal set (spin symmetry only) perform better. In fact S-UHF LM,S-UHF thebestresultsareconsistentwiththeconstraint-release 0.8 SG,S-UHF results, with much smaller statistical error bars. ution 0.7 SG,S-GHF Computational scaling with system size is an impor- b Distri 0.6 n(k)-nED(k) taanCtPMcoCnsicdaelrcautliaotnionfo,rthaenyprimmaanryy-foebrmjeicotns imnveotlhvoindg. thIne m 0.5 0.3 u trialwavefunctionthatneedto be calculatedrepeatedly nt 0.4 0.2 me 0.1 are the overlap hΨT|φi and the one-body Green func- o 0.3 tion Ψ c†c φ / Ψ φ . For a single-determinant Ψ M 0.2 0 suchhasTU|HkF,l|thiechomTp|uitations14scalelikeregularm|eaTni- -0.1 field calculations. For a Ψ in the form of Eq. (4), a 0.1 | Ti -0.2 straightforwardway is to expand it as a linear combina- 0 0.5 1 1.5 2 2.5 3 0 tion of Slater determinants after the reference determi- 0 0.5 1 1.5 2 2.5 3 |k| nant and variational parameters have been determined. This is the approachwe have taken in the present paper FIG. 5: (Color online) Momentum distribution n(k) versus as a proof-of-concept study of the effectiveness of such |k|. Results from CPMC using different trial wave functions wave functions. However, the number of Slater deter- arecomparedwitheachotherandwithexactdiagonalization minants in suchan expansiongrowsrapidly with system (ED). Theinset shows theerror from ED results. The mixed size. Forexample,inthe8 8calculationinTableII,the estimator isusedinsteadofback-propagation,tomagnifythe × effect of the trial wave function. The system is 4×4 with SG,S-UHF trial wave function involved 14,336 determi- 7↑,7↓ and U =8. The same run parameters are used as in nants. An alternative is to use Eq. (4), and incorporate Fig. 3. thesymmetryprojectionintheQMC,computingtheob- jects with numerical integrations over the projection op- As seen from Fig. 5, the calculated momentum distri- erators.27 The computational cost is then proportional bution is incorrect around the Fermi energy when the to that for a simple HF trial wave function, with a pro- UHF trial wave function is used. This is consistent with portionality constant which depends on the symmetries the fact that the UHF is a very poor wave function, and involvedbut only grows modestly with system size (e.g., back-propagationmustbeemployedforobservableswhen linearlyfortranslation). This is a veryappealing feature 8 TABLEII:CPMC ground-stateenergy of theHubbardmodel in squarelattices. Symmetry-projectedHFtrial wavefunctions are used;theground-state in all cases corresponds to a singlet transforming as theB1 irrep of theC4v group with momentum (0,0). Results using the UHF trial wave function are listed for comparison. The last column is release constraint28 results (using the CASSCF trial wave functions), which are the best estimate of the exact ground-state energy. The statistical error bars in theQMC results are on thelast digit and shown in parentheses. SYSTEM U UHF S-UHF LM,S-UHF SG,S-UHF SG,Sz-GHF exact 6×6, 12↑12↓ 4 -42.621(6) -42.605(3) -42.633(3) -42.670(2) -42.679(1) -42.669(8) 6×6, 12↑12↓ 8 -36.68(7) -36.20(2) -36.94(1) -37.188(8) -37.383(4) -37.41(6) 8×8, 22↑22↓ 4 -75.886(8) -75.875(7) -75.900(6) -75.904(3) -75.893(9) of the symmetry-projected HF wave functions as Ψ , functions with a computational cost which scales mod- T | i in contrastwith, for instance,the CASSCF-type ofwave estly with system size. The development in this paper functions, and offers the promise of scalable calculations is thus potentially a scalable and systematically improv- with our approach to reach large system sizes. able quantum Monte Carlo approach for extended systems. Moreover, one can further increase the quality of the symmetry-projected trial wave function by doing IV. CONCLUSION particle number symmetry breaking or by incorporating multi-component expansions such as those considered in We have studied symmetry-projected HF wave func- ourpreviouswork.29 We finallynotethatthecurrentde- tions as trial wave functions for CPMC. The symme- velopments in CPMC have paved the way for accurate tries are restored by projecting from an UHF or GHF many-body calculations of Hamiltonians with spin-orbit type wave function, minimizing the variational energy coupling. aftertheprojection. Thisgivesahierarchyofwavefunc- Acknowledgments tions with increasing quality. The CPMC results using these as trial wave functions generally become increas- ingly more accurate, and the Monte Carlo variance in- WethankS.Chiesa,F.Ma,W.Purwanto,andE.Wal- creasingly smaller for the same amount of samples. It ter for helpful discussions. This research was sup- is shown that CPMC largely restores size-extensivity in portedby DOE (GrantNo.de-sc0008627andDE-FG02- symmetry-projectedHFwavefunctions. Athalf-fillingin 11ER16257) and NSF (Grant No. DMR-1006217). The the Hubbard model, projected wave functions with spin work at Rice University was supported by DOE, Of- symmetryprovideanexcellentdescriptionofthesystems. fice of Basic Energy Sciences, Heavy Element Chemistry Awayfromhalf-filling,theprojectedwavefunctionswith program (Grant No. DEFG02-04ER15523) and DOE- spacegroupsymmetryandspinsymmetryleadtohighly CMCSN (Grant No. de-sc0006650). We also acknowl- accurateresults,oftenwiththeconstrained-pathsystem- edge an INCITE award for computing using resources aticerrorlessthan0.1%inthecorrelation energy. Wave- of the Oak Ridge Leadership Computing Facility at the functions projected from a GHF-type reference tend to Oak Ridge National Laboratory, which is supported by perform even better compared with the ones projected the Office of Science of the U.S. Department of Energy from UHF-types. under Contract No. DE-AC05-00OR22725 and CPD at As discussed above, the symmetry-projected wave CollegeofWilliamandMary. G.E.S.isaWelchFounda- function can be implemented in CPMC as trial wave tion Chair (C-0036). 1 J.G.BednorzandK.A.Mu¨ller,Z.Phys.B64,189(1986). 10 U. Schollwöck, Rev.Mod. Phys.77, 259 (2005). 2 E. Dagotto, Science 309, 257 (2005). 11 U. Schollwöck, Ann.Phys. 326, 96 (2011). 3 M. H.Y.Kamihara, T. Watanabe andH.Hosono, J.Am. 12 R.Blankenbecler, D.J.Scalapino, and R.L. Sugar,Phys. Chem. Soc. 130, 3296 (2008). Rev.D 24, 2278 (1981). 4 G. R. Stewart, Rev.Mod. Phys.83, 1539 (2011). 13 G. Sugiyama and S. Koonin, Ann.Phys. 168, 1 (1986). 5 W. L. McMillan, Phys.Rev. 138, A442 (1965). 14 S.Zhang,chapterintheOpenAccessbook,EmergentPhe- 6 D. Ceperley, G. V. Chester, and M. H. Kalos, Phys. Rev. nomena in Correlated Matter: Modeling and Simulation B 16, 3081 (1977). Vol 3, Eds. E. Pavarini and E. Koch and U. Schollwöck, 7 C. J. Umrigar, K. G. Wilson, and J. W. Wilkins, Phys. http://www.cond-mat.de/events/correl13 (2013). Rev.Lett. 60, 1719 (1988). 15 S. Zhang, J. Carlson, and J. E. Gubernatis, Phys. Rev. B 8 R. J. Bartlett and M. Musial , Rev. Mod. Phys. 79, 291 55, 7464 (1997). (pages 62) (2007). 16 C.-C. Chang and S. Zhang, Phys. Rev. B 78, 165101 9 S. R.White, Phys. Rev.Lett. 69, 2863 (1992). (2008). 9 17 S. Zhang and H. Krakauer, Phys. Rev. Lett. 90, 136401 30 G. E. Scuseria, C. A. Jiménez-Hoyos, T. M. Henderson, (2003). K. Samanta, and J. K. Ellis, J. Chem. Phys. 135, 124108 18 W. Purwanto, H. Krakauer, and S. Zhang, Phys. Rev. B (2011). 80, 214116 (2009). 31 C. A. Jiménez-Hoyos, T. M. Henderson, T. Tsuchimochi, 19 C.-C. Chang and S. Zhang, Phys. Rev. Lett. 104, 116402 and G. E. Scuseria, J. Chem. Phys.136, 164109 (2012). (2010). 32 J. Bonnard and O. Juillet, Phys. Rev. Lett. 111, 012502 20 W. Purwanto, H. Krakauer, Y. Virgus, and S. Zhang, J. (2013). Chem. Phys. 135, 164105 (2011). 33 J.Hubbard,Proc.R.Soc.LondonSer.A276,238(1963). 21 P.RingandP.Schuck,TheNuclearMany-BodyProblem, 34 H. F. Trotter, Proc. Amer. Math. Soc. 10, 545 (1959). Springer-Verlag, New York (1980). 35 M. Suzuki,Comm. Math. Phys. 51, 183 (1976). 22 T. Mizusaki and M. Imada, Phys. Rev. B 69, 125110 36 D. R. Hamann and S. B. Fahy, Phys. Rev. B 41, 11352 (2004). (1990). 23 F. F. Assaad, P. Werner, P. Corboz, E. Gull, and 37 S.R.White,D.J.Scalapino,R.L.Sugar,E.Y.Loh,J.E. M. Troyer,Phys. Rev.B 72, 224518 (2005). Gubernatis, and R. T. Scalettar, Phys. Rev. B 40, 506 24 K. W. Schmid, T. Dahm, J. Margueron, and H. Mu¨ther, (1989). Phys. Rev.B 72, 085116 (2005). 38 M. Calandra Buonaura and S. Sorella, Phys. Rev. B 57, 25 Q. Jie, Phys.Rev.E 77, 026705 (2008). 11446 (1998). 26 D. Tahara and M. Imada, J. Phys. Soc. Jpn. 77, 114701 39 D. A.Huse, Phys.Rev.B 37, 2380 (1988). (2008). 40 J.Carlson,J.E.Gubernatis,G.Ortiz,andS.Zhang,Phys. 27 R. Rodr´ıguez-Guzmań, K. W. Schmid, C. A. Jiménez- Rev.B 59, 12788 (1999). Hoyos,andG.E.Scuseria,Phys.Rev.B85,245130(2012). 41 W. Purwanto and S. Zhang, Phys. Rev. E 70, 056702 28 H. Shiand S. Zhang, Phys.Rev. B 88, 125132 (2013). (2004). 29 R.Rodr´ıguez-Guzmań, C.A.Jiménez-Hoyos,R.Schutski, and G. E. Scuseria, Phys.Rev.B 87, 235129 (2013).