Multi-reference symmetry-projected variational approaches for ground and excited states of the one-dimensional Hubbard model
R. Rodríguez-Guzmán, Carlos A. Jiménez-Hoyos, R. Schutski, Gustavo E. Scuseria
aa r X i v : . [ c ond - m a t . s t r- e l ] J un Multi-reference symmetry-projected variational approaches for ground and excitedstates of the one-dimensional Hubbard model
R. Rodr´ıguez-Guzm´an , , Carlos A. Jim´enez-Hoyos , R. Schutski and Gustavo E. Scuseria , Department of Chemistry, Rice University, Houston, Texas 77005, USA Department of Physics and Astronomy, Rice University, Houston, Texas 77005, USA (Dated: August 7, 2018)We present a multi-reference configuration mixing scheme for describing ground and excited states,with well defined spin and space group symmetry quantum numbers, of the one-dimensional Hubbardmodel with nearest-neighbor hopping and periodic boundary conditions. Within this scheme, eachstate is expanded in terms of non-orthogonal and variationally determined symmetry-projectedconfigurations. The results for lattices up to 30 and 50 sites compare well with the exact Lieb-Wu solutions as well as with results from other state-of-the-art approximations. In addition tospin-spin correlation functions in real space and magnetic structure factors, we present results forspectral functions and density of states computed with an ansatz whose quality can be well-controlledby the number of symmetry-projected configurations used to approximate the systems with N e and N e ± PACS numbers: 71.27.+a, 74.20.Pq, 71.10.Fd
I. INTRODUCTION.
Studies of correlations arising from electron-electroninteractions remain a central theme in condensed mat-ter physis to better understand challenging phenomenasuch as high-Tc superconductivity or colossal magneticresistance. There is a need for better theoretical modelsthat can account for relevant correlations in ground andexcites states of fermionic systems with as much sim-plicity as possible. Within this context, the repulsiveHubbard Hamiltonian has received a lot of attentionsince it is considered the generic model of strongly cor-related electron systems. Hubbard-like models have alsoreceived renewed attention in the study of cold fermionicatoms in optical lattices and the electronic properties ofgraphene. Unlike the one-dimensional (1D) Hubbard model,which is exactly solvable using the Bethe ansatz, anexact solution of the two-dimensional (2D) problem isnot known. Therefore, it is highly desirable to developapproximations that, on the one hand, can capture themain features of the exact 1D Bethe solution and, on theother hand, can be extended to higher dimensions. Forsmall lattices, one can resort to exact diagonalization us-ing the Lanczos method. For larger systems, severalother methods have been extensively used to study the1D and 2D Hubbard models as well as their strong cou-pling versions. Among such approximations, we havethe quantum Monte Carlo, the variational MonteCarlo, the density matrix renormalization group aswell as approximations based on matrix product and ten-sor network states. Both the dynamical mean field the-ory and its cluster extensions have made important contributions to our present knowledge of the Hubbardmodel. Other embedding approaches are also available. Finally, we refer the reader to the recent state-of-the-artapplications of the coupled cluster method to frustratedHubbard-like models.
Although routinely used in nuclear structure physics,especially within the Generator Coordinate Method, symmetry restoration via projection techniques has re-ceived little attention in condensed matter physics. Nev-ertheless, these techniques offer an alternative for ob-taining accurate correlated wave functions that respectthe symmetries of the considered many-fermion prob-lem. The key idea is, on the one hand, to considera mean-field trial state |Di which deliberately breaksseveral symmetries of the original Hamiltonian. On theother hand, the Goldstone manifold ˆ R |Di , where ˆ R rep-resents a symmetry operation, is degenerate and thesuperposition of such Goldstone states, can be usedto recover the desired symmetry by means of a self-consistent variation-after-projection procedure. Sucha single-reference (SR) scheme provides the optimal Ritz-variational representation of a given state by meansof only one symmetry-projected mean-field configura-tion. This kind of SR variation-after-projection scheme,has already been applied to the 1D and 2D Hubbardmodels as well as in quantum chemistry within theframework of the Projected Quasiparticle Theory. One of the main advantages of the symmetry-projectedapproximations is that they offer compact wavefunctions as well as a systematic way to improve theirquality by adopting a multi-reference (MR) approach. Inthis case, a set of symmetry-broken mean-field states |D i i is used to build Goldstone manifolds ˆ R |D i i whose super-position can be used to recover the desired symmetriesof the Hamiltonian. The key idea is then to expanda given state in terms of several symmetry-projected andvariationally determined mean-field configurations. Theresulting wave functions encode more correlations thanthe ones obtained within SR methods while still keepingwell defined symmetry quantum numbers. There are differents flavors of MR approximationsavailable in the literature.
In the present study, weadopt a MR scheme well known in nuclear structurephysics, which to the best of our knowledge has notbeen applied to lattice models. The key ingredient insuch MR scheme is the inclusion of relevant correlationsin both ground and excited states on an equal footing.As a benchmark test, we concentrate on the 1D Hubbardmodel for which exact solutions are known. In particu-lar, we consider the case of half-filled lattices. Neverthe-less, the present MR approximation can be extended tothe 2D case as well as to doped systems with arbitraryon-site interaction strengths.For a given single-electron space we resort to gener-alized HF-transformations (GHF) mixing all quantumnumbers of the single-electron basis states. The corre-sponding Slater determinants deliberately break spin andspatial symmetries of the 1D Hubbard Hamiltonian. We restore these broken symmetries with the help of pro-jection operators. The resulting MR ground state wavefunctions are obtained applying the variational principleto the projected energy.The structure of our MR ground state wave functions isformally similar to the one adopted within the Resonat-ing HF (ResHF) method, i.e., they are expanded interms of a given number of non-orthogonal symmetry-projected configurations. Nevertheless, while in the lat-ter all the underlying HF-transformations and mixing co-efficients are optimized simultaneously, in our casethe orbital optimization is performed sequentially, onlyfor the last added HF-transformation (all our mixing co-efficients are still optimized at the same time) renderingour calculations easier to handle. This is particularly rel-evant for alleviating our numerical effort if one keeps inmind that, for both ground and excited states, we usethe most general GHF-transformations and therefore afull 3D spin projection is required.Our MR scheme is also used to compute spin-spin cor-relation functions (SSCFs) in real space, magnetic struc-ture factors (MSFs) as well as dynamical properties ofthe 1D Hubbard model like spectral functions (SFs) anddensity of states (DOS).
On the other hand, onemay wonder whether there is any relevant information inthe intrinsic symmetry-broken GHF-determinants asso-ciated with our MR wave functions. As will be shownbelow, the structure of such intrinsic determinants canbe interpreted in terms of basic units of quantum fluctu-ations for the lattices considered. In addition to ground state properties, our MRframework treats excited states with well definedquantum numbers as expansions in terms of non- orthogonal symmetry-projected configurations usingchains of variation-after-projection (VAP) calculations.As a byproduct, we also obtain a (truncated) basis con-sisting of a few Gram-Schmidt orthonormalized states, which may be used to perform a final diagonalization ofthe Hamiltonian in order to account for further correla-tions in both ground and excited states.The layout of the theory part of this paper is as follows.First, we introduce the methodology of our MR VAPscheme in Sec.II. Symmetry restoration based on a sin-gle Slater determinant (i.e., SR symmetry restoration) isdescribed in Sec.II A. This section will serve to set our no-tation as well as to introduce some key elements of our 3Dspin and full space group projection techniques. Subse-quently, symmetry restoration based on several Slater de-terminants (i.e., MR symmetry restoration) is discussedin Sec.II B. In particular, the MR description of groundand excited states is presented in Secs.II B 1 and II B 2,respectively. In Sec.II C, we will briefly discuss the com-putation of the SFs and DOS within our theoreticalframework.The results presented in this paper test the perfor-mance of our approximation in a selected set of illus-trative examples. In most cases, calculations have beencarried out for on-site repulsions U = 2 t, t and 8 t takenas representatives of weak, intermediate-to-strong, andstrong correlation regimes. In Sec.III, we first considerthe ground states of half-filled lattices with up to 50sites. We compare our ground state and correlation en-ergies with the exact ones as well as with those obtainedusing other theoretical methods. We then discuss thedependence of the predicted correlation energies on thenumber of non-orthogonal symmetry-projected configu-rations used to expand our ground state wave functions,the computational performance of our scheme as well asthe structure of the intrinsic GHF-determinants result-ing from our MR VAP procedure. Next, we consider theresults of our calculations for SSCFs in real space andMSFs for half-filled lattices with up to 30 sites. Theseresults are compared with density matrix renormalizationgroup (DMRG) ones obtained with the open sourceALPS software. This comparison is very valuable asDMRG represents one of the most accurate approxima-tions in the 1D case. Subsequently, we compare the DOSprovided by our theoretical framework with the exactone, obtained with an in-house diagonalization code, ina lattice with 10 sites. Results for hole SFs are also dis-cussed for a 30-site lattice. We end Sec.III by presentingresults for excitation spectra in various lattices and dis-cussing the structure of the intrinsic GHF-determinantsresulting from our MR VAP procedure for excited states.Finally, Sec.IV is devoted to concluding remarks andwork perspectives.
TABLE I: Ground state energy of the half-filled lattices with N sites = 12 and 20, as predicted with the GHF-FED schemebased on n = 10 GHF-transformations, are compared withexact results for on-site repulsions of U = 2 t, t and 8 t . En-ergies obtained with the RHF and UHF approximations areincluded as a reference. The ratio of correlation energies κ ob-tained with the UHF and GHF-FED aproximations, is com-puted according to Eq.(27). For more details, see the maintext. N sites = 12 κ (%) N sites = 20 κ (%)U=2t RHF -8.9282 -15.2551UHF -9.3379 36.79 -15.6411 24.05GHF-FED -10.0401 99.85 -16.8565 99.79EXACT -10.0418 -16.8599U=4t RHF -2.9282 -5.2550UHF -5.6290 67.65 -9.3821 66.08GHF-FED -6.9201 99.99 -11.4954 99.92EXACT -6.9204 -11.5005U=8t RHF 9.0718 14.7450UHF -2.9532 92.26 -4.9219 92.23GHF-FED -3.9625 99.99 -6.5612 99.96EXACT -3.9626 -6.5699 II. THEORETICAL FRAMEWORK
In what follows, we describe the theoretical frameworkused in the present study. First, SR symmetry restora-tion is presented in Sec.II A. Subsequently, in Sec.II B,we consider our MR scheme to describe both ground(Sec.II B 1) and excited (Sec.II B 2) states of the 1D Hub-bard model. The computation of SFs and DOS is brieflydiscussed in Sec.II C.
A. Single-reference (SR) symmetry restoration
We consider the 1D Hubbard Hamiltonian ˆ H = − t X j,σ n ˆ c † j +1 σ ˆ c jσ + ˆ c † jσ ˆ c j +1 σ o + U X j ˆ n j ↑ ˆ n j ↓ (1)where the first term represents the nearest-neighbor hop-ping ( t >
0) and the second is the repulsive on-site inter-action (
U > operators ˆ c † jσ and ˆ c jσ create and destroy an electron with spin-projection σ = ± / σ = ↑ , ↓ ) along an arbitrary chosenquantization axis on a lattice site j = 1 , . . . , N sites . Theoperators ˆ n jσ = ˆ c † jσ ˆ c jσ are the local number operators.We assume periodic boundary conditions, i.e., the sites j and j + N sites are identical. Furthermore, we assume alattice spacing ∆ = 1. In the standard HF-approximation, the groundstate of an N e -electron system is represented by a Slaterdeterminant |Di = Q N e h =1 ˆ b + h | i in which the energeticallylowest N e single-fermion states (holes h, h ′ , . . . ) are oc-cupied while the remaining 2 N sites − N e states (particles p, p ′ , . . . ) are empty. For a set of single-fermion opera-tors ˆ c † , the HF-quasiparticle operators ˆ b † are given bythe following canonical transformation ˆ b † i = X jσ D ∗ jσ,i ˆ c † jσ (2)where D is a general 2 N sites × N sites unitary matrix,i.e., DD † = D † D = 1. In all the calculations to bediscussed below, we have used generalized HF (GHF)transformations. As it is well known, the most generalGHF-determinant |Di deliberately breaks several sym-metries of the original Hamiltonian.
Typical ex-amples are the rotational (in spin space) and spatialsymmetries. To restore the spin quantum numbers ina symmetry-broken GHF-determinant, we explicitly usethe full 3D projection operator ˆ P S ΣΣ ′ = 2 S + 18 π Z d Ω D S ∗ ΣΣ ′ (Ω) R (Ω) (3)where R (Ω) = e − iα ˆ S z e − iβ ˆ S y e − iγ ˆ S z is the rotation opera-tor in spin space, the label Ω = ( α, β, γ ) stands for theset of Euler angles and D S ΣΣ ′ (Ω) are Wigner functions. To recover the spatial symmetries, we introduce the pro-jection operatorˆ P kmm ′ = 12 N sites X g Γ kmm ′ ( g ) ˆ R ( g ) (4)where Γ kmm ′ ( g ) is the matrix representation of an irre-ducible representation, which can be found by standardmethods, and ˆ R ( g ) represents the corresponding sym-metry operations (i.e., translation by one lattice site andthe reflection x → − x ) parametrized in terms of the la-bel g . The linear momentum k = (2 π/N sites ) ξ is givenin terms of the quantum number ξ that takes the values ξ = − N sites , . . . , N sites Equivalently, itcan take all integer values between 0 and N sites −
1. For k = 0 , π an additional label b = ± x → − x . In whatfollows, we do not explictly write this label b but thereader should keep in mind that it is taken into accountwhenever needed.We introduce the shorthand notation Θ = ( S, k ) forthe set of symmetry (i.e., spin and linear momentum)
TABLE II: Ground state energy of the half-filled lattices with N sites = 30 and 50 predicted with the GHF-FED schemebased on n = 25 GHF-transformations, are compared withexact results for on-site repulsions of U = 2 t, t and 8 t . Re-sults obtained with the UHF-ResHF approximation, basedon n = 30 UHF-transformations, as well as the RHF andUHF energies are also included in the table. The ratio ofcorrelation energies κ obtained with the UHF, UHF-ResHFand the GHF-FED aproximations, is computed according toEq.(27). N sites = 30 κ (%) N sites = 50 κ (%)U=2t RHF -23.2671 -38.7039UHF -23.4792 10.02 -39.1294 12.02UHF-ResHF -25.3436 98.11 -41.9535 91.78UHF-FED -25.3508 98.45 -41.9963 92.99GHF-FED -25.3730 99.50 -42.1219 96.46EXACT -25.3835 -42.2443U=4t RHF -8.2671 -13.7039UHF -14.0732 64.75 -23.4553 65.02UHF-ResHF -17.0542 98.00 -27.9633 95.09UHF-FED -16.9420 96.75 -27.3518 91.01GHF-FED -17.1789 99.39 -27.9788 95.19EXACT -17.2335 -28.6993U=8t RHF 21.7329 36.2961UHF -7.8329 93.65 -12.3048 92.26UHF-ResHF -9.5378 99.04 -15.6422 98.59UHF-FED -9.3524 98.46 -14.8461 97.08GHF-FED -9.7612 99.75 -15.6753 98.65EXACT -9.8387 -16.3842 quantum numbers as well as K = (Σ , m ). The totalprojection operator readsˆ P Θ KK ′ ≡ ˆ P S ΣΣ ′ ˆ P kmm ′ (6)We then superpose the Goldstone manifoldˆ R (Ω) ˆ R ( g ) |Di to recover the spin and spatialsymmetries via the following SR ansatz |D ; Θ; K i = X K ′ f Θ K ′ ˆ P Θ KK ′ |Di (7)where f Θ are variational parameters. Note, that thestate Eq.(7) is already multi-determinantal via theprojection operator ˆ P Θ KK ′ . For a given symmetry Θ,the energy (independent of K) associated with the stateEq.(7) E Θ = f Θ † H Θ f Θ f Θ † N Θ f Θ (8)is given in terms of the Hamiltonian and norm H Θ KK ′ = hD| ˆ H ˆ P Θ KK ′ |DiN Θ KK ′ = hD| ˆ P Θ KK ′ |Di (9)matrices. It has to be minimized with respect to thecoefficients f Θ and the underlying GHF-transformation D . The variation with respect to the former yields thefollowing resonon-like eigenvalue equation (cid:0) H Θ − E Θ N Θ (cid:1) f Θ = 0 (10)with the constraint f Θ † N Θ f Θ = 1 ensuring the orthog-onality of the solutions. On the other hand, the unre-stricted minimization of the energy [Eq.(8)] with respectto D is carried out via the Thouless theorem. For a given symmetry Θ, we only retain the en-ergetically lowest solution of our VAP equations. Both the GHF-transformation D and the mixing coef-ficients f Θ are complex, therefore one needs to mini-mize n var = 2(2 N sites − N e ) × N e + 4 S real variables.We use a limited-memory quasi-Newton method for suchminimization. In practice, the integration over theset of Euler angles in Eq.(3) is discretized. For example,for a lattice with N sites = 30 we have used 13, 26, and 13grid points for the integrations over α , β , and γ , respec-tively. In this case, a total of 263,640 grid points are usedin the discretization of the projection operator of Eq.(6).We have afforded such a task by developing a parallel im-plementation for all the VAP schemes discussed in thispaper. B. Multi-reference (MR) symmetry restoration
For each symmetry Θ, the SR procedure described inSec.II A provides us with the optimal variational repre-sentation of the corresponding ground state via a sin-gle symmetry-projected GHF-determinant. However, asthe lattice size increases one may adopt a MR perspec-tive to keep and/or improve the quality of the wavefunctions.
The key features of our MR approach,known in nuclear structure physics as the FED VAMP (Few Determinant Variation After Mean-field Projec-tion) strategy, for the considered ground states are de-scribed in the next subsection. We use the acronymGHF-FED to refer to it in the present work. On the otherhand, our MR approach for excited states, known as EX-CITED FED VAMP, will be presented in Sec.II B 2.We will use the acronym GHF-EXC-FED to refer to itin what follows.
1. MR symmetry restoration for ground states (GHF-FED)
Our goal in this section is to obtain, through a chainof VAP calculations, non-orthogonal symmetry-projected n k ( % ) N sites =20, U =2 tN sites =20, U =4 tN sites =20, U =8 tN sites =30, U =2 tN sites =30, U =4 tN sites =30, U =8 t FIG. 1: (Color online) The ratio of correlation energies κ obtained with the GHF-FED approximation is plotted as afunction of the inverse of the number of GHF-transformationsfor the half-filled lattices with N sites = 20 and 30. Results areshown for on-site repulsions of U = 2 t, t and 8 t . For moredetails, see the main text. GHF-configurations used to build a MR expansion of agiven ground state with well defined symmetry quan-tum numbers Θ.Suppose we have generated a ground state solution | φ K i = |D ; Θ; K i [Eq.(7)]. Note that at this point,we have added the superscript 1 to explicitly indicatethat only one GHF-transformation has been used withinthe SR approximation discussed in Sec.II A. On the otherhand, the subscript 1 has been added to indicate that theground state is considered. As we will see in Sec.II B 2,this subscript will allow us to distinguish between ground(i.e., i = 1) and excited (i.e., i = 2 , , . . . , m ) states. Onthe other hand, both indices are also added to the (in-trinsic) GHF-transformation to explicitly indicate that itis variationally optimized for the state | φ K i . We thenkeep the transformation D fixed and consider the ansatz | φ K i = X K ′ X i =1 f i Θ1 K ′ ˆ P Θ KK ′ |D i i (11)which approximates the ground state (subscript 1) bymeans of two (superscript 2) non-orthogonal symmetry-projected GHF-determinants. It is obtained applying thevariational principle to the energy functional with re-spect to the last added transformation D and all thenew mixing coefficients f i Θ1 . A similar procedure canbe followed to approximate the ground state by a largernumber of non-orthogonal symmetry-projected configu-rations. Let us assume that n − D n , a new set of mixing coefficients f i Θ1 , and makes the MR GHF-FED ansatz | φ n Θ1 K i = X K ′ n X i =1 f i Θ1 K ′ ˆ P Θ KK ′ |D i i (12)which superposes the Goldstone manifoldsˆ R (Ω) ˆ R ( g ) |D i i . The corresponding energy E n Θ1 = f n Θ † H n Θ f n Θ f n Θ † N n Θ f n Θ (13)is given in terms of the Hamiltonian and norm H n Θ iK,jK ′ = hD i | ˆ H ˆ P Θ KK ′ |D j iN n Θ iK,jK ′ = hD i | ˆ P Θ KK ′ |D j i (14)kernels, which require the knowledge of the symmetry-projected matrix elements between all the GHF-determinants used in the expansion Eq.(12). Thewave function Eq.(12) is determined varying the energyEq.(13) with respect to all the new mixing coefficients f i Θ1 and the last added transformation D n . In the for-mer case, we obtain an eigenvalue equation similar toEq.(10), with the constraint f n Θ † N n Θ f n Θ = 1, whilethe unrestricted minimization with respect to D n is car-ried out via the Thouless theorem. Let us stress that theGHF-FED MR approximation Eq.(12) of a given groundstate enlarges the flexibity in our wave functions to a totalnumber of n var = 2 n (2 N sites − N e ) × N e +4 n S +2( n −
2. MR symmetry restoration for excited states(GHF-EXC-FED)
In this section, we construct non-orthogonalsymmetry-projected GHF-configurations to expanda given excited state. The orthogonalization be-tween ground and excited states is achieved via theGram-Schmidt procedure. As a byproduct, our MRGHF-EXC-FED method also yields a (truncated) basisconsisting of a few orthonormal states which may be usedto diagonalize the Hamiltonian and account for furthercorrelations in both ground and excited states.
Let us assume that we have already obtained a GHF-FED ground state | φ n i = | φ n Θ1 K i [Eq.(12)] along thelines discussed in the previous Sec.II B 1. We then lookfor the first excited state (subscript 2) with the samesymmetry Θ, approximated by a given n number of non-orthogonal symmetry-projected configurations. We startwith the ansatz | ϕ i = α | φ n i + β | φ i (15)
60 120 180 240 300 360number of proccessors102030 s peedup U =4 t N sites =50(a) n t ( s e c ) U =4 t N sites =50(b) FIG. 2: Speedup of a typical GHF-FED calculation is shownin panel (a) as a function of the number of proccessors. Thecorresponding scaling (for a fixed number of proccessors) withthe number of transformations n is presented in panel (b).Results are for the half-filled lattice with N sites = 50 and U = 4 t. where | φ K i has a form similar to Eq.(7) but written interms of the coefficients f and the GHF-determinant |D i . Both α and β can be obtained by requiring or-thonormalization with respect to the ground state thatwe already have. The state Eq.(15) is determined vary-ing the energy functional with respect to f and D .When n − | ϕ n i = α n | φ n i + β n | φ n i (16)where the state | φ n Θ2 K i has a form similar to Eq.(12) butwritten in terms of the new coefficients f i Θ2 and the GHF-transformations D i (i= 1, . . . , n ). Once again, the coef-ficients α n and β n are obtained by requiring orthonor-malization with respect to the ground state we alreadyhave. The wave function Eq.(16) is determined vary-ing the energy functional with respect to the last addedGHF-transformation D n and all the coefficients f i Θ2 .Now, we consider the most general situation in whichthe ground state | ϕ n i = | φ n i as well as a set of m − | ϕ n i , | ϕ n i , . . . , | ϕ n m − m − i , all of them with the same symmetryquantum numbers Θ, are already at our disposal. Eachof these m − m th excited state by n m non-orthogonal symmetry-projected configurations. We also need to ensure orthog-onality with respect to all the m − n m − m th excited state withsymmetry Θ. Then an approximation in terms of n m non-orthogonal symmetry-projected GHF-configurationsis obtained with the GHF-EXC-FED ansatz | ϕ n m m i = m − X i =1 ω mi | φ n i i i + τ m | φ n m m i (17)where the state | φ n m Θ mK i has a form similar to Eq.(12)but is written in terms of new coefficients f i Θ m and GHF-transformations D im (i= 1, . . . , n m ). The coefficients τ m and ω mi read τ m = h φ n m m | (cid:16) − ˆ S m − (cid:17) | φ n m m i − / ω mi = − m − X k =1 (cid:0) A − (cid:1) ik h φ n k k | φ n m m i τ m (18)in terms of the projector (i.e., ˆ S m − = ˆ S m − )ˆ S m − = m − X i,k =1 | φ n i i i (cid:0) A − (cid:1) ik h φ n k k | (19)with the overlap matrix A ik = h φ n i i | φ n k k i . The MRGHF-EXC-FED wave function Eq.(17) is determined byvarying all the coefficients f i Θ m and the last added GHF-transformation D n m m . The energy is E n m Θ m = f n m Θ † H n m Θ f n m Θ f n m Θ † N n m Θ f n m Θ (20)where the Hamiltonian H n m Θ and norm N n m Θ ker-nels account for the fact that m − f i Θ m leads to a generalized eigen-value equation similar to Eq.(10) with the constraint f n m Θ † N n m Θ f n m Θ = 1, while the unrestricted minimiza-tion with respect to the last added GHF-transformation D n m m is carried out via the Thouless theorem.The GHF-EXC-FED scheme outlined in this sectionprovides, for each set of symmetry quantum numbersΘ, a (truncated) basis of m (orthonormalized) states | ϕ n Θ1 K i , . . . , | ϕ n m Θ mK i , each of them expanded by n , . . . , n m non-orthogonal symmetry-projected GHF-determinants,
10 20 30 40 50j-160-80080160 x · (a) 10 20 30 40 50j U =4 t (b) 10 20 30 40 50j (c) 10 20 30 40 50j (d) FIG. 3: (Color online) The quantity ξ i [Eq.(28)] is plotted as a function of lattice site j for some typical symmetry-brokenGHF-determinants | D i i resulting from the GHF-FED VAP optimization for the half-filled lattice with N sites = 50 and U = 4 t .The UHF spin-density wave is plotted in red for comparison. For more details, see the main text. respectively. Finally, the diagonalization of the Hamilto-nian Eq.(1) in such a basis m X j =1 h h ϕ n i i | ˆ H | ϕ n j j i − ǫ Θ α δ ij i C Θ jα = 0 (21)provides ground and excited states | Ω Θ α i = m X j =1 C Θ jα | ϕ n j j i (22)which may account for additional correlations. Never-theless, because many of these correlations have alreadybeen accounted for in the MR expansion of each of the m basis states (as discussed above), one may expect therole of this final diagonalization to be, in general, lessimportant than in the scheme used in Ref. 36.Both the GHF-FED and GHF-EXC-FED VAP approx-imations could be extended to the use of general Hartree-Fock-Bogoliubov (HFB) transformations. This, how-ever, would require an additional particle number sym-metry restoration that increases our numerical effort byaround one order of magnitude and has hence not beenincluded in the present study.
C. Spectral functions and density of states
In this section, we briefly discuss the computation ofthe SFs and DOS within our theoretical framework. Letus assume that for an N e -electron system we have alreadyobtained, along the lines described in Sec.II B 1, a GHF-FED ground state solution | φ n Θ1 K i . For all the latticesconsidered in the present study the ground state has spin S = 0 but not necessarily linear momentum zero [i.e.,Θ = (0 , k )]. In all cases, the ground state transforms asan irrep of dimension 1. Therefore, for this specific case,we can simply write the ground state wave function as | n , N e , k i . The ground state energy will be denoted as E n k .Usually, the SFs are defined as the imaginary part ofthe time-ordered Green’s function and can be calculatedfrom the Lehmann representation. In order to computethem, we approximate the ground states of the ( N e ± ± =(1 / , k ± ). For the ( N e − R (Ω) ˆ R ( g )ˆ b h ( D i ) |D i i inthe ansatz | n T , N e − , k − i = X ihM f i Θ − hM ˆ P Θ − KM ˆ b h ( D i ) |D i i (23)where i = 1 , . . . , n T . The number n T of GHF-transformations used to expand the state Eq.(23) maybe different from the one (i.e., n ) in the GHF-FEDground state wave function. We write ˆ b h ( D i ) to ex-plicitly indicate that holes are made on different in-trinsic determinants |D i i corresponding to the lowest-energy states of the N e -electron system approximatedby a single symmetry-projected configuration along thelines described in Sec.II A. The hole index h runs as h = 1 , . . . , N e . Note, that the label b = ± N e + 1)-electron system we superpose theGoldstone (particle) manifolds ˆ R (Ω) ˆ R ( g )ˆ b † p ( D i ) |D i i andwrite | n T , N e + 1 , k + i = X ipM g i Θ + pM ˆ P Θ + KM ˆ b † p ( D i ) |D i i (24)where the index i runs again as i = 1 , . . . , n T and p = N e + 1 , . . . , N sites . The mixing coefficients f i Θ − and g i Θ + are determined by solving eigenvalue equationssimilar to Eq.(10). This yields a maximun of 2 n T × N e × d hole solutions with energies E n T k − and a maximum of j -0.4-0.3-0.2-0.100.10.20.30.40.50.60.70.80.91 S p i n C o rr e l a t i on F un c t i on U =2 t (a) DMRG
UHF N sites =30 N sites =26 N sites =22 N sites =18 N sites =14GHF-FED5 6 7 8 9 10 11 12 13 14 15 -0.04-0.0200.020.04 -0.4-0.3-0.2-0.100.10.20.30.40.50.60.70.80.91 S p i n C o rr e l a t i on F un c t i on U =4 t (b) -0.4-0.3-0.2-0.100.10.20.30.40.50.60.70.80.911.1 S p i n C o rr e l a t i on F un c t i on U =8 t (c) FIG. 4: (Color online) GHF-FED ground state spin-spin cor-relation functions in real space for half-filled lattices of differ-ent sizes (red, brown, magenta, blue and cyan curves). DMRGvalues are plotted with black triangles. Results are shown forU=2 t (a), 4 t (b) and 8 t (c). In each panel, the inset dis-plays a close-up of the long-range behavior predicted by theGHF-FED (red curve) and DMRG (black triangles) schemescompared with the one obtained within the standard UHF ap-proximation (green curve) in the case of the N sites =30 lattice.For more details, see the main text. n T × (2 N sites − N e ) × d particle solutions with ener-gies E n T k + for each irreducible representation of thespace group. The quantity d is the dimension of thecorresponding irreducible representations, i.e., d = 1 for k ± = 0 , π and d = 2 for k ± = 0 , π . The hole B ( q, ω ) andparticle A ( q, ω ) SFs are then written in their standardform B ( q, ω ) = X k − σ |h n T , N e − , k − | ˆ c qσ | n , N e , k i| × δ (cid:16) ω − E n k + E n T k − (cid:17) A ( q, ω ) = X k + σ |h n T , N e + 1 , k + | ˆ c † qσ | n , N e , k i| × δ (cid:16) ω − E n T k + + E n k (cid:17) (25)and the DOS can be computed as N ( ω ) = X q ( B ( q, ω ) + A ( q, ω )) (26)Due to the finite size of the considered lattices, boththe hole and particle SFs consist of a finite number of δ functions with different weights. Therefore, we introducean artificial Lorentzian width Γ for each state. III. DISCUSSION OF RESULTS
In this section, we discuss the results of our calcula-tions for some illustrative examples. In most cases, weconsider on-site repulsions of U = 2 t, t and 8 t represent-ing weak, intermediate-to-strong (i.e., non-interactingband width) and strong correlation regimes. First, inSec.III A, we consider the ground states of half-filled lat-tices of various sizes. We compare our ground state andcorrelation energies with the exact ones, as well as withthose obtained using other theoretical approaches. Wethen discuss the dependence of the predicted correlationenergies on the number n of non-orthogonal symmetry-projected configurations used to expand our ground statewave functions. The computational performance of ourscheme is also addressed. The structure of the intrinsicGHF-determinants resulting from our GHF-FED VAPprocedure is discussed in Sec.III B. Our results for SS-CFs in real space and MSFs, for lattices with up to30 sites, are presented in Sec.III C. They are comparedwith DMRG results. For all the considered lattices, wehave retained 1024 states in the renormalization proce-dure. In Sec.III D, we compare the DOS provided byour theoretical framework with the exact one, obtainedwith an in-house full diagonalization code in a latticewith N sites = 10. Hole SFs are also discussed in thecase of N sites = 30. Finally, in Sec.III E, we present re-sults obtained for the excitation spectra in lattices with N sites = 12 ,
14 and 20 and also discuss the structure ln N sites M agne t i c S t r u c t u r e F a c t o r ( q = p ) UHF , U =2 tGHF-FED , U =2 tUHF , U =4 tGHF-FED , U =4 tUHF , U =8 tGHF-FED , U =8 t DMRG
FIG. 5: (Color online) The GHF-FED magnetic structurefactor, evaluated at the wave vector q = π , is plotted as afunction of ln N sites for half-filled lattices of different sizes.GHF-FED results are shown for on-site repulsions of U = 2 t (red diamonds), 4 t (black diamonds), and 8 t (blue diamonds).The corresponding DMRG values are plotted with open cir-cles. A straight line has been fitted to guide the eye. Themagnetic structure factors predicted by the UHF approxima-tion for U = 2 t (continuous red curve), 4 t (continuous blackcurve), and 8 t (continuous blue curve) are also included forcomparison purposes. of the underlying symmetry-broken GHF-determinantsresulting from our GHF-EXC-FED VAP procedure forexcited states. A. Ground state and correlation energies
Let us start by considering lattices of 12 and 20 siteswith the GHF-FED scheme discussed in Sec.II B 1. Thecorresponding Θ = (0 , π ) ground states have B symme-try, i.e., they are symmetric under the reflection x → − x .In Table I, we compare the predicted ground state ener-gies with the exact ones. For completeness, we also in-clude energies provided by the standard (i.e., one trans-formation) restricted (RHF) and unrestricted (UHF) HFframeworks. Ours is a VAP approach whose quality canbe checked by studying how well it reproduces the ex-act ground state correlation energies. To this end, weconsider the ratio κ GHF − F ED = E RHF − E GHF − F ED E RHF − E EXACT ×
100 (27)between the GHF-FED and the exact correlation ener-gies. For the UHF approximation, κ UHF is obtained froma similar expression.We observe from Table I that the inclusion of n =10non-orthogonal symmetry-projected configurations with the GHF-FED approach significantly improves correla-tion energies with respect to UHF. In fact, κ GHF − F ED ≥ .
79% in all considered correlation regimes even for N sites = 20 which is out of reach with exact diagonal-ization.In the case of N sites = 14, whose Θ = (0 ,
0) groundstate has A symmetry, i.e., it is symmetric under thereflection x → − x , our calculations with n = 10transformations predict energies of − . t , − . t ,and − . t compared to the exact ones of − . t , − . t , and − . t for U = 2 t, t, and 8 t , respec-tively. This yields κ GHF − F ED values of 99 . , .
97, and99 . For half-filled latticeswith sizes comparable to the ones already mentioned, theGutzwiller method provides κ ratios around 85 , , and50%, respectively. Let us also mention that our GHF-FED energies for N sites = 12 and 14 improve upon pre-viously reported VAP values of − . t and − . t for U = 4 t . Calculations have also been carried out for N sites =16, whose Θ = (0 , π ) ground state has B symmetry.We have obtained ground state energies of − . t for U = t and − . t for U = 4 t while the exactones are − . t and − . t , respectively. Previ-ous DMRG results for this lattice, have been reported inthe literature. For all the lattices with sizes N sites ≤
18 our DMRG calculations, retaining 1024 states in therenormalization procedure, reproduce the exact Lieb-Wuground state energies (to all quoted figures) for the con-sidered U values.The ground state energies for the lattices with N sites =30 and 50 are compared in Table II with the exactones. In this case, the corresponding Θ = (0 ,
0) groundstates have A symmetry. In the same table, we alsopresent ground state energies predicted with the ResHFmethod based on n = 30 UHF-transformations (i.e.,UHF-ResHF). It is very satisfying to observe that boththe GHF-FED and the UHF-ResHF VAP schemes canaccount for κ ≥
98% in a relatively large lattice with N sites = 30. In fact, the GHF-FED scheme provides κ ≥ .
39% with 45,048 variational parameters that rep-resents a small fraction of the dimension of the restricted(i.e., accounting for all symmetries) Hilbert space in thislattice. In this case, our GHF-FED energy also improvesthe variational value − . t obtained in Ref. 35 for U = 4 t using a single symmetry-projected configuration.Note, that the ResHF method is not intrinsicallylimited to the use of UHF-transformations and, there-fore, the UHF-ResHF ground state energies shown inTable II can still be improved by, for example, adopt-ing GHF-transformations as basic building blocks. Onthe other hand, for N sites = 30 our DMRG calculationsprovide the energies − . t , − . t , and − . t for U = 2 t, t, and 8 t , respectively.Let us now comment on our results for N sites = 32whose Θ = (0 , π ) ground state has B symmetry. We0 -3 -2 -1 0 1 2 3 w -U/ ( in t units )04080120160 N ( w ) (a)U=2t -3 -2 -1 0 1 2 3 w -U/ ( in t units ) (b) n , n T =10 ED U=4t -9 -6 -3 0 3 6 9 w -U/ ( in t units ) (c)U=8t FIG. 6: (Color online) The DOS (black) for the half-filled lattice with N sites = 10 at U = 2 t, t, and 8 t is plotted in panels(a), (b), and (c), respectively, as a function of the shifted excitation energy ω − U/ t units). Results have been obtainedby approximating the N e and ( N e ± n = 10 and n T = 10 GHF-determinants. Our results are hardlydistinguishable from the DOS obtained with exact diagonalization (red). A Lorentzian folding of width Γ = 0 . t has beenused. For more details, see the main text. have used n = 25 GHF-transformations. For on-siterepulsions of U = t and 2 t , we have obtained ener-gies of − . t and − . t while the UHF-ResHFones are − . t and − . t , respectively. Theseenergies, should be compared with the exact ones of − . t and − . t as well as with the DMRGvalues − . t and − . N sites = 50. In particular, even when our descrip-tion of the ground state in this lattice is poorer thanin the N sites = 30 case since we have kept the samenumber of GHF-transformations, it is remarkable that weobtain (with 125,048 variational parameters) the values κ GHF − F ED = 96 . , . , and 98 . predicts κ values of around 87 , , and 96%.The corresponding UHF-ResHF values are also listedin Table II.In Fig.1, we have plotted the ratio κ GHF − F ED , as afunction of the inverse 1 /n of the number of transfor-mations n included in the GHF-FED ansatz, for lat-tices with N sites = 20 and 30. They increase smoothlywith the number of non-orthogonal symmetry-projectedconfigurations used to expand the wave function. FromFig.1, it is apparent (see also Tables I and II) that withincreasing lattice size, we need a larger number n ofsymmetry-projected configurations to keep and/or im-prove the quality of the GHF-FED wave functions. Forexample, comparing the N sites = 20 and 30 lattices,we see that in the former n = 10 transformations are enough to obtain κ GHF − F ED ≥ .
79% while in the lat-ter 98 . ≤ κ GHF − F ED ≤ . N sites = 50 case, n = 10 transformations leads to93 . ≤ κ GHF − F ED ≤ . , whereas with n = 25we reach the κ GHF − F ED values shown in Table II.Obviously, as in many other approaches to many-fermion systems, we are always limited to a finite num-ber of configurations in practical calculations. Never-theless, the GHF-FED scheme provides compact groundstate wave functions whose quality can be systematicallyimproved by adding new (variationally determined) non-orthogonal symmetry-projected configurations. In fact,both ours and the ResHF wave functions are noth-ing else than a discretized form of the exact coherent-state representation of a fermion state and, therefore,become exact in the limit n → ∞ . Our aim in thepresent work is not to lower the ground state energy asmuch as possible but to test to which extent our schemecan account for relevant correlations in the consideredlattices. Therefore, for the largest lattices here studied(i.e., N sites = 30 and 50), we have restricted ourselvesin practice to a maximun number n = 25 of GHF-transformations.A few words concerning the computational perfor-mance of our method are in order here. In panel (a)of Fig.2, we have plotted the speedup of a typical calcu-lation as a function of the number of proccessors. Resultsare shown for N sites = 50 and U = 4 t but similar behav-ior was also obtained for U = 2 t and 8 t . As demonstratedin the plot, the GHF-FED speedup grows linearly withthe number of processors used in the calculations. Onthe other hand, panel (b) of the same figure shows that(for a fixed number of proccessors) an efficient imple-mentation of our variational scheme scales linearly withthe number n of GHF-transformations used while theResHF scaling is quadratic. Concerning the scaling ofour method with system size, Fig.1 shows that as the1 -4 -3 -2 -1 0 w -U/ ( in t units ) p / - p / p / - p / p U =2 t (a) p / p / U =2 t (d) -4 -3 -2 -1 w -U/ ( in t units ) U =4 t (b) N sites =30 U =4 t (e) n T =5 n T =15 n T =25 -6 -5 -4 -3 w -U/ ( in t units ) U =8 t (c) U =8 t (f) FIG. 7: (Color online) The hole SFs for the half-filled lattice with N sites = 30 at U = 2 t, t, and 8 t are plotted in panels (a),(b), and (c) as functions of the shifted excitation energy ω − U/ t units). Results have been obtained by approximatingthe N e and ( N e ± n = 25 and n T = 25 GHF-determinants. The hole SFs for momenta identical to the Fermimomentum are displayed in brown color. The shapes of some selected hole SFs (i.e., k = 0 , π/ π/ N e ± n T = 5 (blue), 15 (red) and 25 (black) but the ground stateof the N e -system always with n =25 GHF-transformations, are compared in panels (d), (e), and (f). A Lorentzian folding ofwidth Γ = 0 . t has been used. For more details, see the main text. system becomes larger a larger number of transforma-tions is required to keep the quality of our wave func-tions. We cannot currently determine how the numberof transformations scales with system size as this wouldrequire to consider larger lattices than the ones studiedin this paper. B. Structure of the intrinsic determinants andbasic units of quantum fluctuations
An interesting issue is whether there is any relevantinformation in the symmetry-broken (i.e., intrinsic) de-terminants |D i i resulting from the GHF-FED VAP opti-mization. We are interested in comparing the structureof these determinants with the spin-density wave solutionobtained with the standard UHF approximation. Here,one should keep in mind that a variationally optimizedGHF-determinant has the same energy as the optimalUHF one. We have studied the quantity ξ i ( j ) = ( − ) j − hD i | S ( j ) |D i i · hD i | S (1) |D i i (28)where j = 1 , . . . , N sites is the lattice index while i =1 , . . . , n enumerates the GHF-determinants in the GHF-FED ground state (subscript 1) solution. Among the n = 25 transformations D i used for the N sites =50 lat-tice, we have selected some typical examples to plot the quantity ξ i ( j ). Results are displayed in panels (a) to (d)of Fig.3. Other determinants |D i i , not shown in the fig-ure, exhibit the same qualitative features. Similar resultsare also found for other on-site repulsions, as well as forother lattices.For the standard UHF spin-density wave solution, thequantity ξ i ( j ) has nearly constant positive values plot-ted with red lines in Fig.3. A very different behavior ap-pears in the intrinsic GHF-determinants |D i i associatedwith the GHF-FED solution. First, we observe a broadspin feature distributed all over the lattice, which is aconsequence of the richer spin textures provided by theuse of GHF-transformations and full 3D spin projection.In addition, pairs of points (black squares) appear where ξ i ( j ) changes its sign (i.e., the spin-density wave reversesit phase). These defects of the spin-density wave phaserepresent soliton-antisoliton ( S − S ) pairs in the case ofhalf-filled lattices. In particular, our analysis ofthe charge densities ρ i ( j ) = 1 − P σ hD i | ˆ n jσ |D i i revealthat they correspond to neutral S − S pairs. Let usstress that the presence of at least one S − S pair is agenuine VAP effect appearing even if we approximate agiven ground state within a SR framework, as discussedin Sec.II A.Furthermore, Fig.3 illustrates how the S − S pairsappear at different lattice locations j with varying dis-tance R S − S among the members of the pairs. Thelatter represents the breathing motion of the S − S S =0, k = p -6.9093t S =1, k =0-6.6314t S =0, k =0-6.4710t S =1, k =5 p / S =1, k = p / S =0, k = p -6.9201t S =1, k =0-6.6694t S =0, k =0-6.4987t S =1, k =5 p / S =1, k = p / S =0, k = p -6.9204t S =1, k =0-6.6701t S =0, k =0-6.4993t S =1, k =5 p / N sites =12 S =0, k =0-8.0577 S =1, k = p -7.7595t S =0, k = p -7.6290t S =1, k = p / S =1, k =6 p / S =0, k =6 p / S =0, k =0-8.0874t S =1, k = p -7.8228t S =0, k = p -7.6846t S =1, k = p / S =1, k =6 p / S =0, k =6 p / S =0, k =0-8.0883t S =1, k = p -7.8322tEXACT N sites =14 FIG. 8: The energies of some selected states obtained withinthe GHF-EXC-FED approximation for the half-filled latticeswith N sites = 12 and 14 are compared with the ones providedby single-reference VAP (SR VAP) calculations (only 3D spinand linear momentum projections) as well as with the exactones from Lanczos diagonalization. Results are shown for U = 4 t . pairs. The S − S pairs are present in all the intrinsicdeterminants |D i i associated with the GHF-FED expan-sion, which as already mentioned above, superposes theGoldstone manifolds ˆ R S (Ω) ˆ R ( g ) |D i i containing defectsin the spin-density wave. We are then left with an in-tuitive physical picture in which the soliton pairs can beregarded as basic units of quantum fluctuations in ourGHF-FED states. On the other hand, the interferencebetween S − S pairs belonging to different symmetry-broken determinants |D i i is accounted for in our calcula-tions through a resonon-like equation similar to Eq.(10).This interpretation has been suggested in previous stud-ies with the ResHF method. C. Spin-spin correlation functions and magneticstructure factors
Let us now consider the ground state spin-spin corre-lation functions (SSCFs) in real space. For a given set ofsymmetry quantum numbers Θ, they can be computedas F n Θ m ( j ) = h φ n Θ1 K | S ( j ) · S (1) | φ n Θ1 K ih φ n Θ1 K | φ n Θ1 K i (29)Note that if a wave function has good spin, as it isthe case with the GHF-FED one, the SSCFs have to bethe same for all the members of a (2S+1)-multiplet and, therefore, they cannot depend on the Σ quantum number.However, a dependence with respect to the particular row m of the space group irreducible representation that weare using in the projection still remains and is explicitlyincluded in F n Θ m ( j ).The SSCFs corresponding to the ground states for N sites = 14 , , , , and 30, approximated by n =10 , , , , and 25 GHF-transformations, respectively,are depicted in panels (a), (b) and (c) of Fig.4. In thesame figure, we have also plotted the values resultingfrom our DMRG calculations. We observe a good agree-ment between the GHF-FED and DMRG SSCFs withslight deviations for the N sites = 30 lattice at U=8t,which can be improved by increasing the number of trans-formations. In particular, both SSCFs display a rapiddecrease for j ≤
3. A similar feature has been studiedin previous works.
Regardless of the on-site inter-action, the short range part of the SSCFs runs parallelfor the lattices considered in Fig.4, pointing to convergedbehavior as a function of lattice size. Morevover, the midand long range amplitude of the SSCF for a given latticeincreases with increasing U values.In each panel of Fig.4, the inset displays a close-up ofthe long range behavior of the SSCF predicted by theGHF-FED and DMRG approximations compared withthe one obtained within the standard UHF approachfor N sites = 30. As can be observed, the amplitude ofthe UHF SSCF remains constant for j ≥ ) and having pure spin states(i.e., no spin contamination ). Our wave functions meetboth conditions.The magnetic structure factors (MSFs), evaluated atthe wave vector q = π , can be computed as S n Θ m ( π ) = 1 N sites X ij ( − ) i + j h φ n Θ1 K | S ( i ) · S ( j ) | φ n Θ1 K ih φ n Θ1 K | φ n Θ1 K i (30)and the ones corresponding to the ground states for N sites = 14 , , , , and 30 are displayed in Fig.5 asfunctions of ln N sites . The corresponding DMRG resultsare shown in the same plot. We have also included theUHF MSFs for comparison purposes. At variance withthe UHF MSFs which diverge exponentially, both theGHF-FED and DMRG results display an almost linearbehavior. A previous work has shown that the SSCFsin real space behave for a half-filled system as ≈ ( ln σ j ) /j .This implies that as a funcion of the lattice size, the MSFsshould behave as ln σ N sites . In Fig.5, we have simplyfitted a straight line using the DMRG MSFs to guidethe eye. We have not attempted to determine logarith-mic corrections as this would require larger lattices thanthose studied in the present paper.3 (a) (b) (c)(d) (e) (f) π/ π/ π/ π/ π U = 2 tU = 2 t π/ π/ π/ π/ π U = 4 tU = 4 tN sites = 14 N sites = 14 π/ π/ π/ π/ π U = 8 tU = 8 t U = 2 tU = 2 t U = 4 tU = 4 tN sites = 20 N sites = 20 S = 0 S = 1 S = 2 U = 8 tU = 8 tk ∆ E ( i n t un i t s ) FIG. 9: (Color online) Energy spectrum obtained via Eq.(21) for the half-filled lattices with N sites = 14 [panels (a), (b), and(c)] and 20 [panels (d), (e), and (f)]. For each irreducible representation of the space group, the lowest-energy and first excitedstates with the spins S = 0 (red bars),1 (blue bars), and 2 (green bars) have been plotted. Results are shown for U = 2 t, t, and8 t . The exact dispersion curves for N sites → ∞ (thin black lines) are also included. In order to guide the eye, the lowest-lyingstates with spin S = 0 ( S = 1) have been connected by long (short) dashed lines. For more details, see the main text. D. Spectral functions and density of states
In panels (a), (b), and (c) of Fig.6, we have plotted(black) the DOS N ( ω ) for N sites = 10. In the same fig-ure, we have also plotted (red) the exact DOS obtainedwith an in-house full diagonalization code. There is excel-lent agreement between ours and the exact DOS concern-ing the position and relative heights of all the prominentpeaks. Both ours and the exact DOS exhibit the particle-hole symmetry well known for half-filled systems and asplitting into the lower and upper Hubbard bands. TheHubbard gap between these bands increases with larger U . Both dynamical cluster approximation and cel-lular dynamical mean field theory studies suggest thatthis gap is preserved for any finite value of the on-siteinteraction at sufficiently low temperatures even in thethermodynamic limit, with U = 0 t being the only singu-lar point.Tendencies to spin-charge separation as well as otherrelevant shadow features inside the Brillouin zone, similarto the ones expected in the infinite- U limit of the 1DHubbard model, have been found in previous cellulardynamical mean field theory and cluster perturbationtheory studies of the spectral weigths in the case offinite on-site repulsions. In panels (a), (b), and (c) of Fig.7, we have plotted the hole SFs for the N sites = 30lattice.The first feature observed from Fig.7 is the Hubbardgap opening at the Fermi momentum k F = 7 π/
15. Thespectral weight concentrates on the prominent peaks be-longing to the spinon band. Our calculations for smallerlattices with N sites = 14 and 20 indicate that this spinonband is quite stable in terms of lattice size althoughthe relative height of its peaks decreases for increas-ing U values. The holon singularities are clearly visiblein some of the SFs shown in Fig.7 for linear momenta − π/ < k < π/
2. On the other hand, the holon bandscan also be followed for linear momenta k > π/ k < − π/
2. They are the mirror images of the oneswith opposite ω − U/ and become apparentfor U = 4 t and 8 t . However, besides the spinon band,the most relevant feature in our SFs is the very extendeddistribution of the spectral weight for linear momenta − π/ < k < π/
2. The comparison with our SFs for N sites = 14 and 20 , obtained with n = 10 and n T = 10GHF-transformations, reveals that the increase of latticesize produces more pronounced shadow features due tothe fragmentation of the spectral strength over a widerinterval of ω − U/ The shapes of some selected hole SFs are compared inpanels (d), (e), and (f) of Fig.7. As can be seen frompanel (d), n T = 15 transformations are enough to ac-count for all relevant details of the SFs shown in panel(a) for U = 2 t . On the other hand, panels (e) and (f)show that a larger number of transformations is requiredfor U = 4 t and 8 t . In particular, increasing the numberof transformations for the ( N e ± n T = 5to n T = 15 and/or 25 leads to a shift of the main peaksand redistributes the spectral strength of some of thepeaks found in the SFs for n T = 5 as a result of thesmall number of configurations used in the calculations.This explains the differences between ours and the SFsreported, for the same lattice at U = 4 t , in the previousVAP study of the 1D Hubbard model using only n = 1and n T = 5 GHF-transformations. E. Excitation spectra
In this section, we consider the low-lying excitationspectra obtained for N sites = 12 , , and 20 with theGHF-EXC-FED scheme discussed in Sec.II B 2. For eachirreducible representation of the space group, we havecomputed the lowest-energy and first excited states withspins S = 0 , , and 2. Each of these states has beenapproximated by 10 non-orthogonal symmetry-projectedGHF-determinants. A final 2 × | ϕ n =10 , Θ1 K i and first excited | ϕ n =10 , Θ2 K i states are veryweakly coupled through the Hamiltonian. Due to this,the energies corresponding to the states | Ω Θ1 K i and | Ω Θ2 K i resulting from the 2 × | ϕ n =10 , Θ1 K i and | ϕ n =10 , Θ2 K i . However, this cannot be anticipated apriori and the final diagonalization of the Hamiltonianshould always be carried out.In Fig.8, we compare the energies of some selectedstates for N sites = 12 and 14 with the ones obtainedin the previous variational study of the 1D Hubbardmodel where 3D spin and linear momentum projectionswere carried out. The exact results in Fig.8 correspondto Lanczos diagonalizations. As can be observed, our MRcalculations, where in addition to 3D spin projection thefull space group of the 1D Hubbard model is taken intoaccount, improve the energies reported in Ref. 35 forboth ground and excited states.In Fig.9, we show the low-lying spectrum obtainedvia Eq.(21), for N sites = 14 [panels (a), (b), and (c)]and 20 [panels (d), (e), and (f)] taken as representativesexamples of systems whose Θ = (0 ,
0) and Θ = (0 , π )ground states have A and B symmetries, respectively.We observe that both the lowest-lying singlet and tripletstates obtained in our calculations nicely follow the sine-like dispersion trend in the exact curve for N sites → ∞ . The anomaly observed in the GHF-EXC-FED k -dispersion for the lowest-energy singlets and triplets hasalso been found in previous studies within the ResHFframework as well as in finite versions of the exactLieb-Wu solutions. For any finite U value, the exact N sites → ∞ curves exhibit gapped excitations, excep-tion made for the k = 0 and k = π states which aredegenerate. In our calculations such a degeneracy is bro-ken due to finite size effects. However, we observe thatfor a given finite U value, the energy difference betweenthe lowest singlet and triplet states decreases with in-creasing lattice size. For example, for U = 2 t we haveobtained ∆ E s − t = 0 . t in the case N sites = 14 while∆ E s − t = 0 . t for N sites = 20. For increasing on-site respulsions, irrespective of the lattice size, an overallcompression of the spectra takes place. This is consistentwith the fact that in the limit U → ∞ all the configura-tions shown in Fig.9 should become degenerate. F. Structure of the intrinsic determinants andbasic units of quantum fluctuations in theGHF-EXC-FED wave functions
In Sec.III B, we have discussed the structure of theintrinsic determinants associated with the GHF-FEDstates. Here, we pay attention to the symmetry-brokenones used to expand the GHF-EXC-FED wave functions.To illustrate our results, we consider states belonging tothe spectrum shown in panel (e) of Fig.9. In particular,we have plotted in panels (a) to (d) of Fig.10 the quanti-ties ξ i ( j ) (black) and ξ i ( j ) (red) computed [see, Eq.(28)]with some of the n = 10 and n = 10 symmetry-brokendeterminants |D i i and |D i i used to expand the lowest-energy | ϕ n =10 , Θ1 K i and first excited | ϕ n =10 , Θ2 K i states withΘ = (1 ,
0) and A symmetry. Other determinants, notshown in the figure, exhibit the same qualitative features.Similar results are also obtained for U = 2 t and 8 t as wellas for other lattices.As can be observed, both ξ i ( j ) and ξ i ( j ) display de-fects similar to the ones already discussed for the S = 0ground states provided by the GHF-FED approximation(see, Fig.3). From this we conclude that not only theground but also the excited state wave functions providedby our MR VAP scheme superpose Goldstone manifoldsbuilt in terms of intrinsic GHF-determinants containingdefects (i.e., solitons) that can be regarded as basic unitsof quantum fluctuations. In general, the intrinsic de-terminants associated with different symmetry-projectedstates may develop local variations of the charge densityas seen (dashed blue curve) from Fig.10 where we havealso plotted the quantity ρ i ( j ) = 1 − P σ hD i | ˆ n jσ |D i i . IV. CONCLUSIONS
The accurate description of the most relevant correla-tions in the ground and low-lying excited states of a given5 j -200-160-120-80-4004080120160200 x · U=4t (a) 2 4 6 8 10 12 14 16 18 20 j U=4t(b) 2 4 6 8 10 12 14 16 18 20 j U=4t (c) 2 4 6 8 10 12 14 16 18 20 j U=4t (d)
FIG. 10: (Color online) Structure of some typical symmetry-broken GHF-determinants |D i i (black) and |D i i (red) used toexpand the lowest-energy and first excited states with spin S = 1, linear momentum k = 0 and A symmetry. The chargedensity (dashed blue) corresponding to the determinants |D i i is also included in each panel. Results are shown for the half-filledlattice with N sites = 20 at U = 4 t . For more details, see the main text. many-fermion system, with as few configurations as pos-sible, is a central problem in quantum chemistry, solidstate, and nuclear structure physics. In the present study,we have explored a VAP MR avenue for the 1D Hubbardmodel. The main acomplishments of the present workare listed below.(i) We have presented a powerful methodology of aVAP MR configuration mixing scheme, originally devisedfor the nuclear many-body problem, but not yet consid-ered to study ground and excited states, with well definedsymmetry quantum numbers, of the 1D Hubbard modelwith nearest-neighbor hopping and periodic boundaryconditions. Both ground and excited states are expandedin terms of non-orthogonal and Ritz-variationally opti-mized symmetry-projected configurations. The simplestructure of our projected states allows an efficient par-allelization of our variational scheme, which scales lin-early with the number of processors as well as with thenumber of transformations used in the calculations. Themethod also provides a (truncated) basis consisting ofa few Gram-Schmidt orthonormalized states. This ba-sis may be used to diagonalized the Hamiltonian to ac-count, in a similar fashion, for additional correlations inthe ground and excited states with well defined symmetryquantum numbers.(ii) We have shown that our MR approximation givesaccurate ground state energies and correlation energiesas compared with the exact Lieb-Wu solutions for rel-atively large half-filled lattices up to 30 and 50 sites.The comparison with other theoretical approaches alsoreveals that our scheme can be considered as a reason-able starting point for obtaining correlated ground statewave functions in the case of the 1D Hubbard model.We have computed the full low-lying spectrum for the N sites = 14 and 20 lattices. The momentum dispersionof the lowest-lying singlet and triplet states follows theexact shape predicted by the Lieb-Wu solution in thethermodynamic limit. With increasing U we also observea general compresion of the spectrum. iii) From the analysis of the structure of the intrinsicdeterminants associated with our MR ground and excitedstate wave functions, we observe that they all containdefects (i.e., solitons) that can be regarded as basic unitsof quantum fluctuations for the considered lattices.(iv) Our results for the ground state SSCFs in realspace show long range decay that is not observed in theUHF case. The MSFs computed from such correlationfunctions show a behavior approximately linear in lnN sites consistent with previous results available in theliterature.(v) Our approximation also allows to compute SFs andthe DOS. To this end, we considered ans¨atze, whose flex-ibility is determined by the numbers n and n T of HF-transformations used to expand the wave functions ofsystems with N e and ( N e ± N sites = 10 we have compared the DOS pre-dicted within our approach with the one obtained usinga full diagonalization and found an excellent agreementbetween both. For a larger lattice with N sites = 30 ourscheme provides hole SFs that agree qualitatively wellwith the ones obtained with other approximations andexhibit tendencies beyond a simple quasiparticle distri-bution.We believe that the finite size calculations discussedin the present study already show that VAP approx-imations, based on MR expansions in terms of non-orthogonal symmetry-projected HF-determinants, rep-resent useful tools that complement other existing ap-proaches to study the physics of low-dimensional corre-lated electronic systems. Within this context, the schemepresented in this work leaves ample space for furtherimprovements and research. First, the number of non-orthogonal symmetry-projected configurations used inthe corresponding MR expansions can be increased toimprove the quality of our wave functions. Second, wecould still incorporate particle number symmetry break-ing (i.e., general HFB-transformations) and restoration(i.e., particle number projection) to access even more cor-6relations. Third, our scheme can be easily extended tothe 2D case as well as to doped systems with arbitraryon-site interaction strengths. Our approximation is alsogeneral enough so as to be implemented for the molecularHamiltonian as well as for lattices like the honeycomb,the Kagome or the Shastry-Sutherland ones. The sameVAP MR scheme can also be applied to study frustratedHubbard models in the 1D and 2D cases. Finally, the MRscheme discussed in the present work could also be usedas a powerful solver in the framework of fragment-bathembedding approximations. In particular, it could re-place exact diagonalizations for fragment sizes where it isnot feasible while still providing highly correlated (frag-ment) wave functions. Obviously, a careful analysis ofthe corresponding symmetries should be carried out ineach case. Work along these avenues is in progress and will be reported elsewhere.
Acknowledgments
This work was supported by the Department of En-ergy, Office of Basic Energy Sciences, Grant No. DE-FG02-09ER16053. GES is a Welch Foundation Chair (C-0036). Some of the calculations in this work have beenperformed at the Titan computational facility, Oak RidgeNational Laboratory National Center for ComputationalSciences, under project CHM048. One of us (R.R-G.)would like to thank Prof. K. W. Schmid, Institut f¨urTheoretische Physik der Universit¨at T¨ubingen, for valu-able discussions. E. Dagotto, Rev. Mod. Phys. , 763 (1994). J. G. Bednorz and K. A. M¨uller, Z. Phys. B , 189 (1986). E. Dagotto, Science , 257 (2005). J. Hubbard, Proc. R. Soc. (London) A , 238 (1963). R. J¨ordens, N. Strohmaier, K. G¨unter, H. Moritz and T.Esslinger, Nature , 204 (2008); U. Schneider, L. Hack-erm¨uller, S. Will, Th. Best, I. Bloch, T. A. Costi, R. W.Helmes, D. Rasch and A. Rosch, Science , 1520 (2008);I. Bloch, J. Dalibard and W. Zwerger, Rev. Mod. Phys. ,885 (2008). A. H. Castro Neto, F. Guinea, N. M. R. Peres, K., S.Novosolev and A. K. Geim, Rev. Mod. Phys. , 109(2009). E. H. Lieb and F. Y. Wu, Phys. Rev. Lett. , 1445 (1968). H. Bethe, Z. Phys. , 205 (1931). G. Fano, F. Ortolani and A. Parola, Phys. Rev. B , 1048(1992). F. H. L. Essler, H. Frahm, F. G¨ohmann, A. Kl¨umper andV. E. Korepin,
The One-Dimensional Hubbard Model (Uni-versity Press, Cambridge, 2005). Quantum Monte Carlo Methods in Physics and Chemistry edited by M. P. Nightingale and C. J. Umrigar, NATOAdvanced Studies Institute, Series C: Mathematical andPhysical Sciences (Kluwer, Dordrecht, 1999), Vol. 525. H. De Raedt and W. von der Linden,
The Monte CarloMethod in Condensed Matter Physics , edited by K. Binder(Springer-Verlag, Heidelberg, 1992). E. Neuscamman, C. J. Umrigar and Garnet Kin-Lic Chan,Phys. Rev. B , 045103 (2012). S. R. White, Phys. Rev. Lett. , 2863 (1992). J. Dukelsky and S. Pittel, Rep. Prog. Phys. , 513 (2004). U. Schollw¨ock, Rev. Mod. Phys. , 259 (2005). T. Xiang, Phys. Rev. B , 10445 (1996). L. Tagliacozzo, G. Evenbly and G. Vidal, Phys. Rev. B ,235127 (2009); C. V. Kraus, N. Schuch, F. Verstraete andJ. I. Cirac, Phys. Rev. A , 052338 (2010); U. Schollw¨ock,Ann. Phys. , 96 (2010); G. K. -L. Chan and S. Sharma,Ann. Rev. Phys. Chem. , 465 (2011). D. Zgid, E. Gull and G. Chan, Phys. Rev. B , 165128(2012). A. Georges, G. Kotliar, W. Krauth and M. J. Rozenberg,Rev. Mod. Phys. , 13 (1996). T. Maier, M. Jarrell, T. Pruschke and M. H. Hettler, Rev.Mod. Phys. 77, 1027 (2005). T. D. Stanescu, M. Civelli, K. Haule and G. Kotliar, Ann.Phys. 321, 1682 (2006). S. Moukouri and M. Jarrell, Phys. Rev. Lett. , 167010(2001). C. Huscroft, M. Jarrell, Th. Maier, S. Moukouri and A. N.Tahvildarzadeh, Phys. Rev. Lett. , 139 (2001). K. Aryanpour, M. H. Hettler and M. Jarrell, Phys. Rev. B , 085101 (2003). M. Potthoff, Eur. Phys. J. B , 429 (2003). G. Knizia and G. K.-L. Chan, Phys. Rev. Lett. , 186404(2012). R. F. Bishop, P. H. Y. Li, D. J. J. Farnell, J. Richter andC. E. Campbell, Phys. Rev. B , 205122 (2012). P. H. Y. Li, R. F. Bishop, D. J. J. Farnell, and C. E.Campbell, Phys. Rev. B , 144404 (2012). R. Rodr´ıguez-Guzm´an, J. L. Egido and L. M. Robledo,Nucl. Phys. A , 201 (2002) R. Rodr´ıguez-Guzm´an, L. M. Robledo and P. Sarriguren,Phys. Rev. C , 034336 (2012). P. Ring and P. Schuck,
The Nuclear Many-Body Problem (Springer, Berlin, 1980). R.R.Rodr´ıguez-Guzm´an and K.W.Schmid, Eur. Phys. J. A , 45 (2004). J.-P. Blaizot and G. Ripka,
Quantum Theory of FiniteFermi Systems (The MIT Press, Cambridge, MA, 1985). K.W. Schmid, T. Dahm, J. Margueron and H. M¨uther,Phys. Rev. B , 085116 (2005). R. Rodr´ıguez-Guzm´an, K.W. Schmid, C. A. Jim´enez-Hoyos and G. E. Scuseria, Phys. Rev. B , 245130 (2012). G. E. Scuseria, C. A. Jim´enez-Hoyos, T. M. Henderson,K. Samanta and J. K. Ellis, J. Chem. Phys. , 124108(2011). C. A. Jim´enez-Hoyos, T. M. Henderson, T. Tsuchimochiand G. E. Scuseria, J. Chem. Phys. , 164109 (2012). K. Samanta, C. A. Jim´enez-Hoyos and G. E. Scuseria, J.Chem. Theory Comput. , 4944 (2012). K. W. Schmid, Prog. Part. Nucl. Phys. , 565 (2004). H. Fukutome, Prog. Theor. Phys. , 417 (1988); , 342(1989). N. Tomita, Phys. Rev. B , 045110 (2004). S. Yamamoto, A. Takahashi and H. Fukutome, J. Phys.Soc. Jpn. , 3433 (1991). S. Yamamoto and H. Fukutome, J. Phys. Soc. Jpn. ,3209 (1992). A. Ikawa, S. Yamamoto, and H. Fukutome, J. Phys. Soc.Jpn. , 1653 (1993). N. Tomita and S. Watanabe, Phys. Rev. Lett. , 116401(2009). T. Mizusaki and M. Imada, Phys. Rev. B , 125110(2004). T. Mizusaki and M. Imada, Phys. Rev. B , 014421(2006). N. Tomita, Phys. Rev. B , 075113 (2009). J. L. Stuber and J.Paldus,
Symmetry Breaking in the In-dependent Particle Model . Fundamental World of Quan-tum Chemistry: A Tribute Volume to the Memory of Per-Olov L¨owdin; Edited by E. J. Brandas and E. S Kryachko(Kluwer Academic Publishers: Dordrecht, The Nether-lands, 2003). A. L. Fetter and J. D. Walecka,
Quantum Theory of Many-Particle systems (McGraw-Hill, New York, 1971). B. Horovitz,
Solitons , edited by S. E. Trullinger , V. E. Za-kharov and V. L. Pokrovsky (Elsevier Science Publishers,Amsterdam, 1986). A. F. Albuquerque et al., J. Magn. Magn. Mater. 310, 1187(2007); B. Bauer et al., J. Stat. Mech. P05001 (2011). For a recent account on symmetry-projection based onnon-unitary HF-transformations the reader is referred toC. A. Jim´enez-Hoyos, R. Rodr´ıguez-Guzm´an and G. E.Scuseria, Phys. Rev. A , 052102 (2012). A.R. Edmonds,
Angular Momentum in Quantum Mechan-ics , Princeton University Press, Princeton (1957). N.W. Ashcroft and N.D. Mermin,
Solid State Physics ,Saunders College, 1976. M.C. Gutzwiller, Phys. Rev. Lett. , 159 (1963). D.C. Liu and J. Nocedal, Math. Program. B , 503(1989). W. Metzner and D. Vollhard, Phys. Rev. B , 7382(1988); F. Gebhard and D. Vollhard, Phys. Rev. B ,6911 (1988). H. Yokoyama and H. Shiba, J. Phys. Soc. Jpn. , 3582(1987). F. Satoh, M. Ozaki, T. Maruyama and N. Tomita, Phys.Rev. B , 245101 (2011). S. Nishimoto, E. Jeckelmann and F. Gebhard, Phys. Rev.B , 165114 (2002). A.M. Perlemov, Sov. Phys. Usp. , 703 (1977). V. Bach, E.H. Lieb and J.P. Solovej, J. Stat. Phys. , 3(1994). H. Shiba and P. A. Pincus, Phys. Rev. B , 1966 (1972). M. Imada, N. Furukawa and T. M. Rice, J. Phys. Soc. Jpn. , 3861 (1992). A. Go and G. S. Jeon, J. Phys.: Condens. Matter ,485602 (2009). K. Penc, K. Hallberg, F. Mila and H. Shiba, Phys. Rev.Lett. , 1390 (1996); J. Favan, S. Haas, K. Penc, F. Milaand E. Dagotto, Phys. Rev. B , R4859 (1997); A. Parolaand S. Sorella, Phys. Rev. B , R13156 (1992); M. Ogata,T. Sugiyama and H. Shiba, Phys. Rev. B , 8401 (1991);Phys. Rev. B , 2326 (1990). D. S´en´echal, D. Perez and M. Pioro-Ladri´ere, Phys. Rev.Lett. , 522 (2000). F. Woynarovich, J. Phys. C , 5293 (1983). K. Hashimoto, Int. J. Quant. Chem. , 633 (1986). C. A. Jim´enez-Hoyos, R. Rodr´ıguez-Guzm´an and G. E.Scuseria, in preparation. B. Sutherland and B. S. Shastry, J. Stat. Phys.33