Compact numerical solutions to the two-dimensional repulsive Hubbard model obtained via non-unitary similarity transformations

Similarity transformation of the Hubbard Hamiltonian using a Gutzwiller correlator leads to a non-Hermitian effective Hamiltonian, which can be expressed exactly in momentum-space representation, and contains three-body interactions. We apply this methodology to study the two-dimensional Hubbard model with repulsive interactions near half-filling in the intermediate interaction strength regime ($U/t=4$). We show that at optimal or near optimal strength of the Gutzwiller correlator, the similarity transformed Hamiltonian has extremely compact right eigenvectors, which can be sampled to high accuracy using the Full Configuration Interaction Quantum Monte Carlo (FCIQMC) method, and its initiator approximation. Near-optimal correlators can be obtained using a simple projective equation, thus obviating the need for a numerical optimisation of the correlator. The FCIQMC method, as a projective technique, is well-suited for such non-Hermitian problems, and its stochastic nature can handle the 3-body interactions exactly without undue increase in computational cost. The highly compact nature of the right eigenvectors means that the initiator approximation in FCIQMC is not severe, and that large lattices can be simulated, well beyond the reach of the method applied to the original Hubbard Hamiltonian. Results are provided in lattice sizes upto 50 sites and compared to auxiliary-field QMC. New benchmark results are provided in the off half-filling regime, with no severe sign-problem being encountered. In addition, we show that methodology can be used to calculate excited states of the Hubbard model and lay the groundwork for the calculation of observables other than the energy.


I. INTRODUCTION
The Fermionic two-dimensional Hubbard model [1][2][3] with repulsive interactions is a minimal model of itinerant strongly correlated electrons that is believed to exhibit extraordinarily rich physical behaviour.Especially in the past thirty years, it has been intensively studied as a model to understand the physics of high-temperature superconductivity observed in layered cuprates 4 .Its phase diagram as a function of temperature, interaction strength and filling includes antiferromagnetism, Mott metal-insulator transition, unconventional superconductivity 5 with d-wave pairing off halffilling, striped phases, a pseudo gap regime, charge and spin density waves 6 .Confronted with such a plethora of physical phenomena, accurate numerical results are indispensable in resolving various competing theoretical scenarios.
Thermodynamic limit extrapolations have been carried out with the aim of assessing the accuracy of the methodologies in various regimes of interaction, filling factor and temperature [27][28][29] .On the other hand each of these methods incur systematic errors which are extremely difficult to quantify and there is an urgent need to develop methods in which convergence behaviour can be quantified internally.
In this paper, rather than attempting a direct numerical attack on the 2D Hubbard Hamiltonian with a given technique, we ask if there is an alternative exact reformulation of the problem, the solution of which is easier to approximate than that of the original problem.If this is the case (and this is obviously highly desirable), it should be demonstrable within the framework of a given technique, without reference to any other method.The physical basis for any observed simplification should be transparent.Such an approach turns out to be possible, at least for intermediate interactions strengths based, on a Gutzwiller non-unitary similarity transformation of the Hubbard Hamiltonian.
The Gutzwiller Ansatz 30,31 and Gutzwiller approximation [32][33][34][35][36] are intensively studied methods to solve the Hubbard model.The parameter of the Ansatz is usually optimized to minimize the energy expectation value by variational Monte Carlo schemes 37,38 based on a single Fermi-sea reference state.It has been long realized that the simple Gutzwiller Ansatz misses important correlations [39][40][41] , especially in the strong interaction regime of the Hubbard model.More general, Jastrow-like 42 , correlators, including density-density 43 and holon-doublon 44,45 , have been proposed to capture more physical features within the Ansatz.In addition the Fermi-sea reference function have been extended to HF spin-density waves 46 and BCS 47 wavefunctions 8,48-53 including anti-ferromagnetic 11,54 and charge order 55 .
An alternative strategy is to use a Gutzwiller correlator to perform a non-unitary similarity transformation of the Hubbard Hamiltonian, whose solution can be well approximated using a Slater determinant.Such an approach is reminiscent of the quantum chemical transcorrelated method of Boys and Handy 56,57 , as well as Hirschfelder 58 , in which a non-Hermitian Hamiltonian is derived on the basis of a Jastrow factorisation of the wavefunction.This idea was applied by Tsuneyuki 59 to the Hubbard model by minimizing the variance of the energy based on projection on the HF determinant.Scuseria and coworkers 60 and Chan et.al. 61ve recently generalized to general two-body correlators and more sophisticated reference states, where the correlator optimization was not performed in a stochastic VMC manner, but in the spirit of coupled-cluster theory, by projecting the transformed Hamiltonian in the important subspace spanned by the correlators.These methods have in common that they are based on a single reference optimization of the correlation parameters and thus the energy obtained is on a mean-field level.We instead would like to fully solve the similarity transformed Hamiltonian in a complete momentum space basis.We will use a single reference optimization, based on projection 60,61 , to generate a similarity-transformed Hamiltonian (non-Hermitian with 3-body interactions), whose ground-state solution (right-eigenvector) will be using the projective FCIQMC 62 method.
The remainder of the paper is organized as follows: In II we recap the derivation of the Gutzwiller similarity transformed Hubbard Hamiltonian and the projective solution based on the restricted Hartree-Fock determinant.We also present analytic and exact diagonalization results, to illustrate the influence of the transformation on the energies and eigenvectors.In III we recap the basics of the FCIQMC method and necessary adaptations for its application to the similarity transformed Hubbard Hamiltonian in a momentum-space basis, named similarity transformed FCIQMC(ST-FCIQMC).In IV we benchmark the ST-FCIQMC method for the exact diagonalizable 18-site Hubbard model, present ground-and excited-state energies.We observe an increased compactness of the right eigenvector of the non-Hermitian transformed Hamil-tonian.We also compare the results obtained with our method for non-trivial 36-and 50-site lattices, at and off half-filling with interaction strengths up to U/t = 4.In V we conclude our findings and explain future applications for observables other than the energy and correct calculation of left and right excited state eigenvectors.

II. THE SIMILARITY TRANSFORMED HAMILTONIAN
We would like to solve for the ground-state energy of the two-dimensional, single-band Hubbard model [1][2][3] with the Hamiltonian in a real-space basis i,σ being the fermionic annihilation(creation) operator for site i and spin σ, n l,σ the number operator, t the nearest neighbor hopping amplitude and U ≥ 0 the on-site Coulomb repulsion.We employ a Gutzwiller-type Ansatz 30,33,63 for the ground-state wavefunction where D is the sum of all double occupancies in |Φ , which are repressed with 0 ≤ g ≤ 1 → −∞ < J ≤ 0.
In the Gutzwiller Ansatz, |Φ is usually chosen to be a single-particle product wavefunction 30,64 , |Φ 0 , such as the Fermi-sea solution of the non-interacting U = 0 system, or other similar forms such as unrestricted Hartree-Fock spin-density waves 46 , or superconducting BCS wavefunctions 48 .The parameter J is usually optimized via Variational Monte Carlo(VMC) 36 , minimizing the expectation value In this work, however, |Φ is taken to be a full CI expansion in terms of Slater determinants with which we aim to solve an equivalent exact eigenvalue equation H denotes a similarity-transformed Hamiltonian.Eq.( 6) is obtained by substituting eq.( 2) as an eigenfunction Ansatz into eq.( 1) and multiplying with e −τ from the left, and due to [n i,σ , τ ] = 0.
The similarity transformation of eq.( 1) moves the complexity of the correlated Ansatz for the wavefunction |Ψ into the Hamiltonian, without changing its spectrum.It is a non-unitary transformation, and the resulting Hamiltonian is not Hermitian.Such similarity transformations have been introduced in quantum chemistry 56,58,65 in the context of a Slater -Jastrow Ansatz, were it is known as the transcorrelated -method of Boys and Handy.It was first applied to the Hubbard model by Tsuneyuki 59 .The transcorrelated method has been quite widely applied in combination with explicitly correlated methods in quantum chemistry [66][67][68][69] , with approximations being employed to terminate the commutator series arising from the evaluation of e −τ Ĥ e τ 70,71 .The explicit similarity transformation of the Hubbard Hamiltonian(1) with a Gutzwiller(2) 59,72 or more general correlator 60,61 , which can be obtained without approximations due to a terminating commutator series, has been solved on a mean-field level 59 .In the present work, we will not restrict ourselves to a mean-field solution, but solve for the exact ground state of H with the FCIQMC method 62,73 .
A. Recap of the derivation of H Tsuneyuki 59 and Scuseria et al. 60 have provided a derivation of the similarity transformed Hubbard Hamiltonian, based on the Gutzwiller and more general two-body correlators, respectively.
Their derivations result in a Hamiltonian expressed in real-space.Here we go one step further and obtain an exact momentum space representation of the similarity transformed Hamiltonian, which is advantageous in the numerical study of the intermediate correlation regime.In this representation, the total momentum is an exact quantum number, resulting in a block diagonalised Hamiltonian.This is computationally useful in projective schemes, especially where there are near-degeneracies in the exact spectrum close to the ground-state, which can lead to very long projection times and be problematic to resolve.Additionally, it turns out that even in the intermediate strength regime, the ground-state right eigenvector is dominated by a single Fermi determinant for the half-filled system.This is in stark contrast with the ground-state eigenvector of the original Hubbard Hamiltonian, which is highly multi-configurational in this regime.
As seen in eq.( 7) we need to compute the following transformation which can be done by introducing a formal variable x and performing a Taylor expansion (cf. the Baker-Campbell-Hausdorff expansion).The derivatives of ( 8) can be calculated as With this closed form (9) the Taylor expansion can be summed up as F (1) = a † i,σ a j,σ e J(n j,σ −n i,σ ) and eq.( 6) takes the final form of 59,60,72 Due to the idempotency of the (Fermionic) number operators, n 2 i,σ = n i,σ , we have for m ≥ 1: With eq.( 11) the exponential factor in eq.( 10) can be calculated as With eq.( 12) we can write the final similarity transformed Hamiltonian as Formulated in a real-space basis the additional factor in eq.( 13) is simply a nearest-neighbor density dependent renormalisation of the hopping amplitude.For large interaction U/t ≫ 1, as already pointed out by Fulde et al. 39 , the simple Ansatz (2) shows the incorrect asymptotic energy behav- better suited choice of basis, due to a dominant Fermi-sea determinant and thus a single reference character of the ground-state wavefunction.Thus we transform eq.( 13) with with M being the size of the system and c k,σ the annihilation(creation) operator of a state with momentum k and spin σ into a momentum space representation.The terms of eq.( 13) become <ij>,σ <ij>,σ <ij>,σ with ǫ k being the dispersion relation of the lattice.The original Hubbard Hamiltonian in k-space is while the similarity transformed Hamiltonian in k-space is a function of the correlation parameter Comparing to the original Hubbard Hamiltonian in k-space (18) to determine the optimal J var .
Projective methods such as the Power method 75 , or a stochastic variant such as FCIQMC 73 , can converge the right/left eigenvectors by multiple application of a suitable propagator, without recourse to a variational optimisation procedure, and this is the technique we shall use here.Since the matrix elements of ( 19) can be calculated analytically and on-the-fly, the additional cost of the 3-body term is no hindrance in our calculations and we do not need to apply additional approximations, unlike other explicitly correlated approaches 76,77 .While complicating the calculation of observables other than the energy, due to the need to have both the left and right eigenvector of the now non-Hermitian Hamiltonian (19), the difference between the left and right eigenvectors actually proves to be beneficial for the sampling of the ground-state wavefunction in the FCIQMC method.
This will be numerically demonstrated below in II C. As a side note, the use of more elaborate correlators, like density-density 43 or holon-doublon 44,45,78 is no hindrance in the real-space formulation of the Hubbard model and is currently being investigated 79 , but in the momentum-space basis would lead to even higher order interactions and have not been further explored.

B. Analytic results for the Hubbard model
As a starting point we optimize the strength of the correlation factor, controlled by the single parameter J from the Ansatz (2), by projecting the single determinant eigenvalue equation H(J) − E |Φ HF = 0 to the single basis of the correlation factor 60,61,72 where . . .c denotes a cumulant expression, where only linked diagrams are evaluated.HF denotes the state with all orbitals with |k| ≤ k F being doubly occupied and k F being the Fermi surface.
Eq.( 22) is similar to a Coupled-Cluster equation.We simply report the results here (further information on the solution of eq.( 22) can be found in Appendix A).For an infinite system at half-filling, and only considering the two-body contribution of eq.( 22), we can express the optimal J which fulfils eq.( 22), and the corresponding total energy per site, as (see B) The results of eq.(22-24) compared with AFQMC results  to the solution of eq.( 22) for the actual finite lattices.At half-filling AFQMC does not suffer from a sign problem 81 and is numerically exact.One can see that the results obtained from eq. (22-24)   capture most of the correlation energy for low U/t as expected.For larger U/t, due to the missing correlation between empty and doubly occupied sites in the Ansatz (2), the energies progressively deteriorate compared to the reference results.The optimal value J opt is also displayed in Table I.We use these values of J obtained by solving eq.( 22) as a starting point for our FCIQMC calculations to capture the remaining missing correlation energy.
To compare this most basic combination of a on-site Gutzwiller correlator and a single restricted Hartree-Fock determinant as a reference in eq.( 22) we show in Table II the percentage of the energy obtained with this method to more elaborate correlators and reference states, for different system sizes M, number of electrons n el and interaction strengths U/t.E (S)U GST in II denotes a on-site Gutzwiller correlator with a (symmetry-projected) unrestricted Hartree-Fock reference state 72 .At half filling and U/t ≤ 4 we can capture more than 80% of the energy obtained with a more elaborate reference determinant.Off half-filling the recovered energy is above 90% up to U/t = 4.For a more dilute filling of n = 0.8, for M = 100 and U/t = 2, the energies agree to better than 99%.Although the absolute error in energy increases off half-filling, as already mentioned in Ref. [60 and 72] the relative error actually decreases 37,82,83 , as can be seen in the comparison with the AFQMC reference results 27,60,80,84 , E ref in Table II.As expected, for U/t > 4 the results from eq.( 22) drastically deteriorate compared to E (S)U GST .
E R/U J in Table II refer to energies obtained with restricted/unrestricted Hartree-Fock reference states with a general two-body correlator 60 , which includes all possible density-density correlations in addition to the on-site Gutzwiller factor.The comparison with E RJ shows that, as already found in Ref. [60], the Gutzwiller factor is by far the most important term in a general two-body correlator for low to intermediate values of U/t ≤ 4, with an agreement of over 98% with E RJ .Off half-filling, as can be seen in the N = 36, n el = 24 and U/t = 8 case, the relative error remains small even for large interaction.The comparison with the available AFQMC reference results 27,60,80,84 , E ref , shows that the solution of eq.( 22) with a on-site Gutzwiller correlator and a restricted Hartree-Fock reference, retrieves above 80% of the energy for U/t ≤ 4.This gives us confidence that the optimal J obtained by this method is appropriate in the context of the Gutzwiller similarity transformed Hamiltonian, which we further solve with the FCIQMC method.
C. Exact diagonalization study.

Due to the non Hermitian nature of H (19) the left and right ground-state eigenvectors |Φ
L/R 0 differ and depending on the strength of the correlation parameter J they can have a very different form.
The most important characteristic for the projective FCIQMC method is the compactness of the sampled wavefunction.As a measure of this compactness we chose the L 2 norm of the exact |Φ L/R 0 contained in the leading HF-determinant |Φ HF and double excitations |Φ ab ij = c † a c † b c i c j |Φ HF (spinindex omitted) thereof, i.e. the sum over the squares of the coefficients of these determinants: As a simple example, in the top panel of figure 1 we show L 2 (0,2) for the 1D half-filled 6-site Hubbard model with periodic boundary conditions at U/t = 4 and k = 0, as a function of the correlation parameter J. J = 0 corresponds to the original Hamiltonian (18).In the bottom panel of fig. 1 the Hartree-Fock energy E HF and results of minimizing the variance of the energy E var by Tsuneyuki 59 , E HF (J opt ) with J opt obtained by solving eq.( 22) and Variational Monte Carlo(VMC) results 85 E V M C are shown.Due to the fact that H is not Hermitian any more, and hence not bounded by below, E HF can drop below the exact ground-state energy E ex , also displayed in fig. 1, so following Tsuneyuki 59 we termed the energy axis "pseudo energy".There is a huge increase in the L 2 (0,2) norm of the |Φ R 0 until an optimal value of J max , close to the J opt obtained by solving eq.( 22), see Tab.VI, where L 2 (0,2) ≈ 1, followed by subsequent decrease.The result obtained by minimizing the energy variance E var is higher in energy and farther from J max than J opt .And, although E V M C is closer to E ex , the optimized correlation parameter obtained by VMC is also farther from J max than J opt .At the same time L 2 (0,2) of |Φ L 0 shows a monotonic decrease with increasing −J.This shows that the amount of relevant information contained within the HF determinant and double excitations thereof can be drastically increased in the right eigenvector, whilst decreased in the left one.For the calculation of the energy, where only the right eigenvector is necessary, a more efficient sampling with the stochastic FCIQMC method should be possible.

III. THE FCIQMC METHOD
The FCIQMC 62,73 method is a projector Monte Carlo method, based on the integrated imaginarytime Schrödinger equation where t is an imaginary-time parameter and |Ψ(t = 0) is an arbitrary initial wave-function with non-zero overlap with |Ψ 0 .One obtains the ground-state energy and wave-function by repeatedly 0.7 0.8 0.9 applying a first-order difference approximated projector of (25) to the initial state Eq.( 27) governs the dynamics of a population of N w signed walkers, which stochastically sample the ground-state wave-function |Ψ 0 .Since the number of states, N d , grows combinatorially with system size, only a stochastic "snap-shot" of |Ψ 0 is stored every iteration, where only states occupied by at least one walker are retained.The diagonal term of eq.( 27), 1 − ∆t (H ii − E S ), increases or decreases the number of walkers on state i.The shift energy E S (t) is dynamically adapted after the chosen number of walkers N w is reached to keep it constant over time.The off-diagonal term, −∆tH ij , creates new walkers from an occupied determinant i to a connected state j.The sum is is sampled stochastically by only performing one of these "spawning" events with a probability and the sign of the new walker is: −sign(H ij ).At the end of each iteration, walkers with opposite sign on the same determinant, which is a reflection of the fermionic sign problem, are removed from the simulation.For sufficiently many walkers the sign problem can be controlled for many systems.In the intermediate to high interaction regime of the Hubbard model, the number of necessary walkers is proportional to the Hilbert space size, making this "original" FCIQMC method impractical.The initiator approximation i-FCIQMC 62 overcomes this exponential bottleneck at the cost of introducing an initiator bias.It does so by allowing only walkers on determinants above a certain population threshold n init to spawn onto empty determinants (thereby dynamically truncating the Hamiltonian matrix elements between low-population determinants and empty ones).This is the source of the initiator error, which can be systematically reduced by increasing the walker population.Nevertheless, convergence can be slow, especially if the ground state wavefunction is highly spread out over the Hilbert space, as is often the case for strongly correlated systems.On the other hand, convergence can be rapidly obtained if the ground-state eigenvector is relatively compact, and does not require any prior knowledge of this fact, nor of the nature of the compactness.
In fact, it is precisely for this reason that the similarity transformations can be of use in the i-FCIQMC method.
In addition to the shift energy E S (t), a projected energy with |D ref being the most occupied determinant in a simulation, is an estimate for the ground-state energy, if |Ψ(t) ≈ |Ψ 0 .An improved estimate (with a smaller variance) can also be obtained by projection onto a multi-determinant trial wave-function Φ trial |, where Φ trial | is obtained as the eigenvector of a small sub-space diagonalised similarity-transformed Hamiltonian.This is particularly useful in open shell problems, where there are several dominant determinants in the ground-state wave-function, and as a result E trial (t) can exhibit notably smaller fluctuations than E P (t).

A. The ST-FCIQMC approach
In variational approaches the lack of a lower bound of the energy due to the non-Hermiticity of the similarity transformed Hamiltonian poses a severe problem.As a projective technique, however, the FCIQMC method has no inherent problem sampling the right ground-state eigenvector, obtaining the corresponding eigenvalue by repetitive application of the projector (26).Additionally, the increased compactness of |Φ R 0 observed in section II C, due to the suppressed double occupations via the Gutzwiller Ansatz, tremendously benefits the sampling dynamics of i-FCIQMC.On the other hand, the implementation of the additional 3-body term in (19) necessitate major technical changes to the FCIQMC algorithm.We changed the NECI 87 code to enable triple excitations.
Due to momentum conservation and the specific spin-relations (σσσ) of the involved orbitals and efficiently analytically calculable 3-body integrals of ( 19), these could be implemented without a major decrease of the performance of the algorithm.In fact the contractions of the 3-body term in (19), namely k = 0 ⊻ k ′ = 0 ⊻ k = k ′ ⊻ q + k ′ = s lead to an O(M) additional cost of the 2-body matrix element, which have the largest detrimental effect on the performance.(There is an O(M 2 ) scaling for the diagonal matrix elements, coming from the k = k ′ = 0 contraction, but this has a negligible overall effect, since we store this quantity for each occupied determinant and is thus not computed often).The additional cost for 2-body integrals is similar to the calculation of 1-body integrals in conventional ab-initio quantum chemistry calculations and unavoidably hampers the performance, but is manageable.Surprisingly, the actual performance improves with increasing strength of the correlation parameter J, even though the three-body interactions are increasing in magnitude.This is due to the following fact: the performance of the FCIQMC method depends heavily on the "worst-case" |H ij |/p(i|j) ratio, where p(i|j) is the probability to spawn a new walker on determinant |D i from |D j and |H ij | is the absolute value of the corresponding matrix element D i | Ĥ|D j .The time-step ∆t of the FCIQMC simulation is on-the-fly adapted to ensure the "worstcase" product remains close to unity, ∆t|H ij |/p(i|j) ≈ 1.An optimal sampling would be achieved, if for every pair (i, j) : p(i|j) ≈ |H ij | and thus ∆t ≈ min(1, E −1 W ). Since H is not Hermitian, the off-diagonal matrix elements are not uniform, as in the original Hamiltonian (18).We therefore need to ensure an efficient sampling by a more sophisticated choice of p(i|j).Hubbard model with U/t = 4.While there is a 7-fold increase of the time per iteration of the mixed scheme compared to the original uniform, the time-step is almost a third larger and the accepted rate of spawning events a third higher.n abort indicates those spawning attempts which originally are proposed in the uniform scheme, but are finally rejected, due to zero matrix elements or are Fermi blocked.This quantity is also decreased by more than a half in the mixed method compared to the uniform original scheme.n accept indicates the number of accepted proposed spawning events and is directly related to the p spawn (28).The choice of the excitation generator is therefore not straightforward and depends on the interaction strength and J: the uniform scheme performs better than expected at small U/t, whilst the mixed scheme performs better at large U/t.

IV. RESULTS
We assessed the performance of initiator ST-FCIQMC (i-ST-FCIQMC) for different Hubbard lattices, as a function of the Gutzwiller correlation factor J. As a starting guide for J, we use J opt obtained by solving eq.( 22) for the specific lattice size M, number of electrons n el and interaction strength U/t, and calculate the ground-state and excited states energies with i-ST-FCIQMC.In particular, we were interested in the rate of convergence of the energy with respect walker number, or in other words, how quickly the initiator error disappeared with increasing walker number.The optimal values of J for each studied system can be found in table VI in the appendix A. All energies are given per site and in units of the hopping parameter t and the lines in the figures 3 to 7 are guides to the eye.

A. 18-site Hubbard model
We first study the 18-site Hubbard model on a square lattice with tilted boundary conditions (see fig. 8), which can be exactly diagonalised: at half-filling and zero total momentum (k = 0) it has a Hilbert space of ∼ 10 8 determinants.All the exact reference results were obtained by a Lanczos diagonalization 88 .For this system ST-FCIQMC could be run either in "full" mode or with the initiator approximation, i-ST-FCIQMC.This enables us to assess to two separate questions, namely the performance of i-ST-FCIQMC with regards to initiator error on the one hand, and compactness of the wave-functions resulting from the similarity transformation (without the complicating effects of the initiator approximation), on the other.to be more efficient in the low U/t regime 39 .Nevertheless, the results shown in fig. 3  (0,2) and J opt is the result of eq.( 22).
since Ĥ † = Ĥ and τ † = τ .And Similar to the exact results for the 6-site model in fig. 1, the right eigenvector shows a huge compactification compared to the original J = 0 result, going from 0.65 to over 0.9.The "optimal" value of J = J max = −0.57444831,where L 2 (0,2) is maximal, is close to the analytical obtained J opt = −0.5234470,indicating that we can simply use J opt without further numerical optimisation of J, and still be close to optimal conditions.Fig. 5 shows the L 2 norm contained in each excitation level relative to the HF determinant for the half-filled, 18-site, U/t = 4 Hubbard model for different values of J.For J = −1/2 there is a huge increase in the L 2 norm of the HF determinant, indicated by an excitation level of 0, while it drops of very quickly for higher excitation levels and remains one order of magnitude lower than the J = 0 result above an excitation level of 5.
Our second analysis on the compactness of |Φ R 0 consisted of running truncated CI 89 calculations, analogous to the CISD, CISDTQ, etc. methods of quantum chemistry.Here we truncate the Hilbert space by only allowing excitation up to a chosen value n trunc relative to the HF determinant.Fig. 6 shows the error of the energy per site as a function of n trunc for different J. from the n trunc = 8 simulations for each value of J, which do not differ much from n trunc = 5 to n trunc = 8 for each simulation.Already at sextuple excitations we are well within error bars of the exact result for J = −1/2, with an error that is two orders of magnitude smaller than the J = 0 result.

Off half-filling 14 e − in 18-sites
We have also investigated the applicability of the i-ST-FCIQMC method to the off-half-filling case, and also to excited states calculations.To this end we calculated the ground, first and second excited states of the 14 e − in 18-sites, U/t = 4, k = 0 system.Such a system can be prepared by removing 4 electrons (2 α and 2 β spins) from the corners of the Fermi-sea determinant, and using this as a starting point for an i-ST-FCIQMC simulation.Excited states are obtained by running multiple independent runs in parallel and applying a Gram-Schmidt orthogonalization to a chosen number of excited states 90 with Pi (t) being the orthogonal projector = 0 for i = j (i and j indicate the excited states), the sampled excited states will in general not be identical to the exact right eigenvectors of H. On the other hand, since the spectrum of H does not change due to the similarity transformation ( 6), the shift energy E S in (33), dynamically adapted to keep the walker number constant, remains a proper estimate for the excited states energy.This interesting fact is developed further in appendix C. Additionally, if the excited state belongs to a different spatial or total-spin symmetry sector the overlap to the ground-state is zero, so our excited state approach within the FCIQMC formalism, via orthogonalisation, correctly samples these orthogonal excited states.Fig. 7 shows the energy per site error of the ground-, first and second excited state of the 14 e − in 18-site, U/t = 4, k = 0 system, compared to exact Lanczos reference results 88 for different values of J versus the walker number N w , obtained via the shift energy E S,i .All states show a similar behaviour of the energy per site error.The closer J gets to the optimal value of J opt = −0.557941for U/t = 4, which is determined for E 0 , one observes that more than an order of magnitude fewer walkers are necessary to achieve the same accuracy as the J = 0 case.This is true for all the states considered.For E 1 , the energy difference of the N w = 10 7 and J = −1/2 calculation is = −0.958466, the relative error is in fact smaller off half-filling.As already mentioned above and shown in table I and II, the projective solution based on the restricted Hartree-Fock determinant ( 22) also yields smaller relative errors off half-filling.These results give us confidence to also apply the i-ST-FCIQMC method to systems off half-filling and for excited states energy calculations.

Symmetry Analysis
As mentioned above, the set of right eigenvectors of a non-Hermitian operator are in general not orthogonal, except when they belong to different irreducible representations and/or total spin symmetry sectors.Here we investigate the interesting influence of the similarity transformation on the symmetry properties of the truncated low-energy subspace of the 14e − in 18-sites system with total k = 0.There are 8 important low energy determinants with the 5 lowest energy k points double occupied and 4 e − distributed among the 4 degenerate orbitals of the corner of the square ) and k 4 = (1, 1) to preserve the total k = 0 symmetry.This is illustrated in fig.8, where red indicates the doubly occupied k-points and green the singly solutions for different values of J.For a large enough correlation parameter J the ground-state of the low-energy subspace resembles the correct symmetry structure as for the full solution.
occupied ones.The point group of the square lattice is D 4h .There are 2 closed shell determinants in this set, with opposite k-points doubly occupied and 6 open-shell determinants with all 4 corners of the BZ singly occupied.Without a correlation parameter all these 8 determinants are degenerate in energy, while with J = 0 this degeneracy is lifted.To study the low energy properties of this system we diagonalized H in this sub-space.Table IV shows the results.We found that with J = 0 the ground state of this subspace has a different spatial and spin symmetry, 5 B 1g , than the ground state of the full system, which belongs to 1 A 1g .At approximately J ≈ −0.71 there is a crossover and the subspace ground state changes to 1 A 1g symmetry.The first excited state in the subspace is then the 5 B 1g , which is also the symmetry of the first excited state of the full system and the 2nd excited state is of 1 B 2g symmetry, the same as 2nd excited state of the not truncated system.Therefore the similarity transformation not only ensures a more compact form of the ground-and excited state wavefunctions, but also correctly orders the states obtained from subspace diagonalisations.The implication is that, in the off half-filling Hubbard model , the structure of ground state has very important contributions arising from high-lying determinants, so much so that they are necessary to get a qualitatively correct ground-state wavefunction (i.e. one with the correct symmetry and spin).With the similarity transformed Hamiltonian, however, this is not the case.Even small sub-space diagonalisations yield a ground-state wavefunction with the same symmetry and spin as the exact one.In other words, the similarity transformation effectively downfolds information from higher lying regions of the Hilbert space to modify the matrix elements between the low-lying determinants.Since the structure of the ground-state eigenvector already has the correct symmetry (and therefore signs) in small subspaces, the rate of convergence of the solution with respect to the addition of further determinants is much more rapid.We believe this is a crucial property which leads to the observed greatly improved convergence rate of the i-ST-FCIQMC method in the off half-filling regime.respectively.The 36-site lattice is studied with periodic and mixed, periodic along the x-axis and antiperiodic boundary conditions along the y-axis.The red points in the 18-site lattice indicate the doubly occupied states and the green points the singly occupied states in the sub-space study in IV A 2.

B. Results for the 36-and 50-site Hubbard model
To put the i-ST-FCIQMC method to a stern test, we applied it to two much larger systems, namely 36-site and 50-site lattices, which are well beyond the capabilities of exact diagonalisation.
In the case of the 36-site (6 × 6) lattice, we considered two boundary conditions, namely fully periodic (PBC) and a mixed periodic-anti-periodic (along the x-and y-axes respectively), the latter being used in some studies to avoid degeneracy of the non-interacting solution 82 .We considered two fillings, namely half-filling and 24e − , at U/t = 2 and U/t = 4.The optimal J opt was determined by solving eq.( 22) and is listed in table VI in the appendix A. For the 6 × 6 by lattice we compared our results to AFQMC calculations 80 , which are numerically exact at half-filling 81 .The results are shown in table V.While the original i-FCIQMC method shows a large error even at walker numbers up to N w = 5 • 10 8 the i-ST-FCIQMC method agrees with the AFQMC reference to within one σ error bars in all but one case (PBC U/t = 4 half-filled), where the agreement is within 2σ.Even in that case the energies agree to better than 99.8%.The small discrepancy could be due to this system being strongly open-shell, making equilibration more challenging.
The 50-site Hubbard lattice corresponds to a 5 √ 5 × 5 √ 5 tilted square, which has been widely investigated using the AFQMC method.We considered half-filling and various off half-filling, n el = 26, 42, 44, 46 and 48 cases for U/t = 1, 2, 3 and 4 and calculated the ground-state energy.The optimal J are listed in table VI in the appendix A. This system size, especially with increasing U/t and off half-filling, was previously unreachable with the FCIQMC method.We compare our halffilling results to AFQMC [91][92][93] reference values, which do not have a sign problem at half-filling 81 .The remaining sources of error are extrapolation to zero temperature and finite steps, both of which are expected to be very small.Off half-filling, exact AFQMC results are not available, and we compare against constrained-path AFQMC(CP-AFQMC) 94,95 and linearized-AFQMC(L-AFQMC) 96 .
Table V shows the results for various fillings and U/t values the reference calculations, the original i-FCIQMC and the i-ST-FCIQMC method.We converged our results for this system size up to a walker number of N w = 10 9 .We can see that the original i-FCIQMC method performs well for the weakly correlated half-filled U/t = 1 system, but fails to reproduce the reference results at U/t = 2 for this system size, and the discrepancy worsens with increasing interaction.The i-ST-FCIQMC method on the other hand agrees within error bars with the reported reference calculation up to U/t = 3 at half-filling.Similar to the half-filled 36-site lattice, the i-ST-FCIQMC results are slightly below the AFQMC reference results at U/t = 4, which could be a finite temperature effect of the AFQMC reference results.
We investigated the half-filled 50-site U/t = 4 system further by performing the convergence of a truncated CI expansions, similar to the 18-site lattice.The results are shown in Fig. 9.The convergence with excitation level truncation shows that convergence occurs from above, and at 6fold excitations we are converged to statistical accuracy to the fully unconstrained simulation.The energy at 6-fold truncation is indeed slightly below the AFQMC result, although the discrepancy is small (approximately 0.1%).It is intriguing that the CI expansion of the 50-site lattice is converged at 6-fold excitations, which is the same as observed for the 18-site lattice.This suggests that linear solutions to the similarity-transformed Hamiltonian may be size-consistent to a greater degree than similar truncations to the original untransformed Hamiltonian.This question however is left for a future study.
For U/t = 4 off half-filling the i-ST-FCIQMC energies are consistently slightly above the reference AFQMC results.However the approximations in the off half-filling AFQMC calculations can lead to energies below the exact ones.For example 84 , CP-AFQMC on a 4×4 lattice with 14 e − and U/t = 4 gives an energy of −0.9863(1) compared to an exact energy of −0.9840, i.e. an overshoot 0.2%.
Similar overshoots are observed at other fillings.In the off half-filling regime in the 50-site system at U/t = 4, CP-AFQMC overshoots our i-ST-FCIQMC results by similar amounts.Therefore our results are in line with ED results for smaller lattices, and thus represent a new set of benchmarks for the off half-filling 50-site Hubbard model.for various interaction strengths U/t, number of electrons n el and periodic (PBC) and mixed (anti-)periodic boundary conditions along the (y-)x-axis, obtained with the initiator FCIQMC and the i-ST-FCIQMC method compared with available (CP-)AFQMC and linearised-AFQMC reference results 80,84,[96][97][98] .The differences to the AFQMC reference energies are displayed as ∆E.The correlation parameter J was chosen close to the optimal J opt obtained by solving eq.( 22) listed in table VI of App.A for the specific U/t value.An initiator threshold of n init = 2.0 was chosen and convergence of the energy up to a walker number of N w = 10 We have used a projective solution based on the restricted Hartree-Fock determinant to obtain an optimized Gutzwiller correlation parameter.For low to intermediate interaction strength, this method generally recovers over 80% of the ground-state energy.Based on this mean-field solution we derived a similarity transformed Hubbard Hamiltonian, generated by the Gutzwiller Ansatz, in a momentum-space basis.We solved for the exact ground-and excited states energy of this non-Hermitian operator with the FCIQMC method.We have shown that the right eigenvector Similar to the 18-site system at half-filling (Fig. 6.), the energies are well-converged at 6-fold excitations.
of the non-Hermitian Hamiltonian has a dramatically more compact form, due to suppression of energetically unfavorable double occupancies, via the Gutzwiller Ansatz.This increased compactness of the right eigenvectors allowed us to solve the Hubbard model for system sizes, which were previously unreachable with the i-FCIQMC method.We benchmarked our results with highly accurate AFQMC reference results and find extremely good agreement at and off half-filling up to interaction strengths of U/t = 4.We hope this combination of a similarity transformation based on a correlated Ansatz for the ground-state wavefunction and subsequent beyond mean-field solution with FCIQMC can aid the ongoing search for the phase diagram of the two-dimensional Hubbard model in the thermodynamic limit.
An important extension of the present work will be to compute observables other than the energy.To As already observed in IV A, applying H with −J yields the left eigenvector |Φ L = e τ |Ψ .To perform this in the FCIQMC we only need to run two independent simulations in parallel, as is aleady done in replica-sampling of reduced density matrices 99 , where the two runs use an opposite sign of the correlation parameter J. Observables, Ô, which commute with the chosen Gutzwiller correlator [ τ , Ô ] = 0, such as the double occupancy n ↑ n ↓ , can be calculated by the 2-body RDM obtained with the left and right eigenvector with normalized Φ L | Φ R = 1 and p, q, r and s denoting spin-orbital labels in the momentum space.Non-commuting observable, [ τ , Ô ], have to be similarity transformed Ō = e −τ Ô e τ and might require higher order density matrices.
Simultaneous calculation of the left eigenvectors |Φ i L also allows us to obtain the correct excited state wave-functions, in addition to the already-correct excited state energy via the shift energy E S,i mentioned in IV A 1 and App.C, in the following manner: For m excited states we run 2m simulations in parallel, where every odd numbered calculation solves for a right eigenstate |Φ i R , which is orthogonalized against all |Φ j L with E j < E i .And vice versa, every even numbered simulation solves for a left eigenvector |Φ i L , orthogonalised against each |Φ j R with E j < E i .In this shoelace-manner m left and right excited state eigenvectors are obtained based on the bi-orthogonal property of left and right eigenvectors of non-Hermitian Operators Φ i L Φ j R = 0 for i = j.Results on observables other than the energy and correct left and right eigenvectors of excited states will be reported in future work.
To perform accurate thermodynamic-limit extrapolations, we also need to reduce the finite size errors of the kinetic term in 7.This can be done by twist averaged boundary conditions   With the coordinate rotation (B5), the integrand of T 1 can be factorized as and T 1 can also be found as a function of |k x | and |k y | In a similar way T 2 can be calculated as The exchange part of the three body contribution in (22) to the correlation energy can be calculated as (using here again the rotation (B5) for p) The final results are and the summations can also be calculate as integrals (e J − e −J = 0, (B19) which, for small U/t, can be approximated as At half-filling Hartree-Fock energy of the original Hubbard Hamiltonian (18), with k = 0 in the two-body term, in the thermodynamic limit (TDL).The additional contributions arising due to the similarity transformation p,q,k,σ ǫ p+k n p,σ n q+k,σ n q,s HF (B23) can be estimated, with cosh(J − 1) ≈ J 2 for small J and eq.(B13), as Hence, the energy per site in the TDL for an un-polarized system at half-filling is given by  33) as a valid estimator for the excited state energies.In fig. 10 the difference to the exact energy, obtained by the projected e p and shift e s energy estimator, for the first 10 states of the 1D 6 e − in 6 site, periodic, U/t = 4, k = 0 Hubbard model with a correlation parameter J = −0.1 are shown.Also shown is the difference of the sum of the overlap of the i-th excited states to all lower lying states j with E j < E i , for the exact right eigenvectors obtained by exact diagonalization and the sampled eigenvectors within FCIQMC But as the i-th excited state is only orthogonalised to all the lower lying excited states to converge to the next higher energy governed by the dynamics (26) and the spectrum of the Hamiltonian (1), which is unchanged by the similarity transformation ( 6), the shift energy remains a good energy estimator.This can clearly be seen in fig. 10, as the shift energy is a good estimate of all the targeted eigenstates.
The only exception in fig.To see why the shift energy is a valid estimate for the exact excited states energy, let's look at the where |Ψ i is the i-th right eigenvector of Ĥ.We now want to show that there exists a vector |Φ i , which is a eigenvector of the composite operator Pi Ĥ with the same eigenvalue where Pi is the Gram-Schmidt projector (C3) and |Φ 0 = |Ψ 0 , which creates an orthonormal basis out of the linear-independent, but not necessarily orthonormal set {|Ψ i }.We assume all states to be normalized.Multiplying eq.(C4) with Pi from the left, we obtain And we assume |Φ i to be the desired eigenvector of Pi Ĥ.To show that we plug (C6) into eq.(C5) with b ij = Φ j | Ψ i .We can express |Φ j in eq.(C7) and all subsequent appearances of |Φ k with k < i as |Φ k = Pk |Ψ k until we reach |Ψ 0 .So the remaining thing to show is that Pi |Ψ j = 0 for i > j.
For i > j is easy to show since {|Φ j } is a orthonormal basis.We prove Pi |Ψ j = 0, j < i by induction.For Let's assume Pi |Ψ j = 0 for i < j, performing the induction step i → i + 1 yields |Ψ j |Φ i = 0 (C10) where we used the Hermiticity P † i = Pi and idempotency P 2 i = Pi of the projection operator.With Pi |Ψ j = 0 eq.(C7) gives the desired Pi Ĥ |Φ i = E i |Φ i . (C12) And this eigenvector |Φ i of the composite operator Pi Ĥ is the stationary vector we sample in FCIQMC.Since it has the same eigenvalue E i , we obtain the correct excited state energy estimate from the shift energy E S i in the propagator (C2).Since the same argument holds for the long-time limit of the projection Qi (E S i ) This |Φ i is sampled by the walkers in a FCIQMC simulation and the shift energy E S i (t) is adapted to keep the walker population fixed.The projected energy is in general not a good energy estimate, since Where we can estimate the overlap Φ j | Ψ i from the orthogonalisation procedure.
Actually for the correct projected energy one needs to calculate Unfortunately the numerator of eq.(C24) takes the following form To calculate Φ j | Ĥ|Φ i we would need the transition (reduced) density matrices (t-(R)DM) between all states j < i.And for the similarity transformed momentum-space Hubbard Hamiltonian even up to the 3-body t-RDM.So we have to rely on the shift energy to yield the correct excited state energy in the ST-FCIQMC method our apply the mentioned shoelace technique in V.

FIG. 1 .
FIG. 1. Top: L 2 norm within the HF determinant and double excitations, L 2 (0,2) , for a half-filled 6-site Hubbard chain with periodic boundary conditions at U/t = 4 and k = 0 for the left |Φ L 0 and right |Φ R 0 ground-state eigenvector of H as a function of −J.Bottom: The Hartree-Fock energy E HF (J) = Φ HF | H(J)|Φ HF as a function of −J.The dotted line indicates the exact ground-state energy E ex for this system, and since H is not Hermitian E HF can drop below the exact energy.Also indicated are the results of minimizing the variance of the energy of the similarity transformed Hamiltonian E var by Tsuneyuki 59 , the result of solving eq.(22) E HF (J opt ) and the result from a VMC optimization E V M C 85 .The vertical dotted line indicates J max where L 2 (0,2) of |Φ R 0 is maximal.All energies are in units of t and the two panel share the x-axis.

FIG. 2 .
FIG. 2. Histogram of | Hij |/p ij ratios for the half-filled, 50-site, U/t = 4 Hubbard model with periodic boundary conditions for uniform, weighted and mixed generation probabilities.

Fig. 3 2 10 3 10 4 10 5 10 6 5 FIG. 3 .
Fig.3shows the error (on a double-logarithmic scale) of the energy per site in the initiator calculation, as a function of walker number.The left panel shows results for the U/t = 2 system.As one can see there is a steep decrease in the error and even with only 10 4 walkers, for a correlation parameter of J = −1/4 (close to the J opt ) the error is below 10 −4 .At 10 6 walkers it is well below 10 −6 , almost two orders of magnitude lower than the original (i.e.J = 0) Hamiltonian at this value of N w .This also confirms the assumption that the chosen Ansatz for the correlation function (2) is particularly useful in the low U/t regime.Results for an intermediate strength, U/t = 4, are shown in the right panel of fig.3.Compared to the U/t = 2, more walkers are needed to achieve a similar level of accuracy.The two sources for this behaviour are: Firstly, i-FCIQMC calculations on the momentum-space Hubbard model are expected to become more difficult with increasing interaction strength U/t, due to the enhanced multi-configurational character of the ground-state wave-function.Secondly, the chosen correlation Ansatz (2) is proven

FIG. 4 .
FIG.4.L 2 norm captured within the HF determinant and double excitations, L 2 (0,2) , for the half-filled 18-site Hubbard model at U/t = 4 for the left and right ground-state eigenvector of the non-Hermitian similarity transformed Hamiltonian(19) as a function of −J.The results were obtained by a non-initiator ST-FCIQMC calculation.|Φ L was sampled by running the simulation with positive J, which corresponds to conjugating H. J max indicates the position of the maximum of L 2 (0,2) and J opt is the result of eq.(22).

1 FIG. 5 .
FIG. 5.The L 2 norm contained in each excitation level relative to the HF determinant, indicated by excitation level 0, for the half-filled, 18-site, U/t = 4, k = 0 Hubbard model for different values of J.The inset shows the tail of the same data on a logarithmic scale.

5 FIG. 6 .
FIG. 6. Error of the energy per site versus the excitation level truncation n trunc in the half-filled, 18-site, U/t = 4, k = 0 Hubbard model for different values of J.The inset shows the absolute error on a logarithmic scale.The dashed lines in the inset indicate the statistical error for the n trunc = 8 results for each value of J.

10 3 10 4 10 5 10 6 10 7 10 E 1 10 3 10 4 10 5 10 6 10 7 (c) E 2 J 5 FIG. 7 .
FIG. 7. Energy per site error compared to exact Lanczos results 88 for the 14 e − in 18-site, U/t = 4, k = 0 Hubbard model for the (a) ground-, (b) first and (c) second excited state as a function of walker number N w .All three panels share the x and y axes.The dashed lines indicate the statistical errors of the N w = 10 7 for each value of J.

FIG. 9 .
FIG.9.The energy per site versus the excitation level truncation n trunc in the half-filled, 50-site, U/t = 4, k = 0 Hubbard model.AFQMC reference98 and non-truncated i-ST-FCIQMC results are also shown.

a
anti-periodic BC along y-axis b open-shell k = 0 reference is possible for non-Hermitian operators, and is the case for states 3, 4 and 5 shown in fig.10, indicated by a large value of ∆O i , since within FCIQMC the incorrect Φ tried to be enforced.The partially incorrect wave-function form is additionally marked by an large error in the projected energy e p compared to the exact result.
FIG. 10.Error of the first 10 eigenstate energies obtained by the projected energy e p and shift energy e s for the 6 e − in 6 site 1D periodic Hubbard model at U/t = 4 and k = 0 compared to exact diagonalization results.The horizontal dashed lines indicate the averaged statistical errors.The green pluses show the difference of the exact overlaps to the overlaps obtained within FCIQMC, see eq. (C1).A correlation parameter of J = −0.1,initiator threshold of n init = 1.2 and a maximum walker number of N w = 10 5 were used.(a) Exact overlap is ill-defined due to degeneracy of states 6 and 7.

i = 1
we have P1 |Ψ 0 = |Ψ 0 − Ψ 0 | Ψ 0 |Ψ 0 = 0. (C9) i − j<i b ij D I | Ĥ|Φ j c I,i − j<i d I,j b ij = E i c I,i − j<i b ij d I,j E P j c I,i − j<i d I,j b ij (C16) with c I,i = D I | Ψ i , b ij = Φ j | Ψ i , d I,j = D I | Φ j (C17)and |D I being the reference determinant of state i.With eq.(C16) and knowledge of the exact eigenfunctions {|Ψ i } the excited state energy could be calculated asE i = E P i c I,i − j<i d I,j b ij + j<i b ij d I,j E P j c −1 I,i .(C18)Forstates where D I | Φ i ≈ c I,i and b ij ≈ 0 the projected energy remains a good estimator for the exact E i .But especially in cases where the exact right eigenvectors are not orthogonal to all lower lying ones, as demonstrated in fig.10, the projected energy should not be trusted.Another correction for the projected energy would beE P i = D I | Ĥ|Φ i D I | Φ i (C19) D I | Φ i E P i = D I | Ĥ Pi |Ψ i = D I | Ψ I E i − j<i Φ j | Ψ i D I | Ĥ|Φ j = D I |Φ j E P j (C20)

TABLE I .
Ground-state energy per site for the half-filled 16 × 16 square lattice with periodic (PBC) and

TABLE II .
27,60,80,84 the total energy obtained with the Gutzwiller Ansatz(2) based on the Hartree-Fock determinant(22)compared with a Gutzwiller correlator with a unrestricted Hartree-Fock reference E U GST and subsequent symmetry projection E SU GST72and a general two-body correlator with a Hartree-Fock reference E RJ and unrestricted Hartree-Fock reference E U J 60 and numerically exact AFQMC reference results27,60,80,84for different number of sites M , number of electrons n el and interaction strengths U/t.M n el U/t %E U GST %E SU GST %E RJ %E U J %E ref the spectral width of Ĥ.If the energy shift E S = E 0 , convergence to a non-diverging and non-zero solution can occur.In practice the shift is dynamically adapted to keep the walker number, explained below, constant, which corresponds to keeping the L 1 norm of the sampled wave-function constant.If the sampled wave-function is a stationary solution to the projector, adapting E S (t) to keep the L 1 norm constant guarantees E S (t) → E 0 .|Ψ(t) is expanded in an orthonormal basis of N d Slater determinants, |Ψ(t) = N d i c i (t) |D i and the working equation for the coefficients c i is

TABLE III
. Ratio of the time-step ∆t, time per iteration ∆t, aborted n abort and accepted excitations n accept of the different excitation generation probabilities compared to the J = 0 uniform reference for the half-filled, 50 site, U/t = 4 Hubbard model with periodic boundary conditions.

TABLE IV .
Irreducible representations and spin symmetry of the ground-state E 0 and first two excited states E 1 , E 2 of the U/t = 4, k = 0, 14 e − in the 18-site Hubbard model for the full and sub-space(subsp.)

TABLE V .
Zero temperature, k = 0 ground-state energy results for the 36-site and 50-site Hubbard model 9was checked.
80,[100][101][102], which is readily applicable for the similarity transformed Hamiltonian in FCIQMC, and will be reported in future work.

TABLE VI .
J opt obtained by solving eq.(22) for the specific lattice sizes, fillings and U/t values used in this paper.J ex is the value, which sets set the J-dependent Hartree-Fock energy, E HF J , to the exact energy, if available, or to the AFQMC reference energies for larger systems. ) Hermitianoperator, since they are sampled by orthogonalizing the n-th excited state to all lower energy states m < n.But it turns out that we are still able to use the dynamically adapted shift energy E S i of eq. ( 90Appendix C: Excited states As discussed in IV A 1 the right eigenvectors of a non-Hermitian operator H |Φ R i = E i |Φ R i are in general not orthogonal to each other.And hence the way excited states are obtained withthe FCIQMC method90should in general not be applicable to excited states of a non- 90, which could be misleading, is state number 7, which appears to have a large error in ∆O i , but the projected energy is still a good estimator for the energy.This comes from the fact that state 6 and 7 are actually degenerate and thus the exact eigenvectors |Φ ex Lapack 105 are an arbitrary linear combination and could be chosen to be both orthogonal to the states i < 6.The n-th excited state in FCIQMC is obtained90by|Φ n (t + δt) = Pn (t + δt) 1 − δt Ĥ − E S Schmidtprojector, which removes all contributions of lower lying states |Φ m and thus orthogonalises |Φ n to each state with E m < E n .For the set of right eigenvector of a non-Hermitian Hamiltonian the assumption of them being orthogonal to each other does not hold in general.So this method of obtaining the excited states of H should in principle not work.But the results above indicate, that the shift energy still provides a correct energy estimate.
C13) with stationary |Ψ i for E S i = E i .There is an eigenvector |Φ i of the composite operator Pi Qi (E S i ) |Φ i = |Φ i (C14) for E S i = E i with |Φ i = P i |Ψ i , since Pi Qi (E S i ) |Ψ i = Pi |Ψ i .(C15)