Spectral Functions from Auxiliary-Field Quantum Monte Carlo without Analytic Continuation: The Extended Koopmans' Theorem Approach
Joonho Lee, Fionn D. Malone, Miguel A. Morales, David R. Reichman
SSpectral Functions from Auxiliary-Field Quantum Monte Carlo without AnalyticContinuation: The Extended Koopmans’ Theorem Approach
Joonho Lee, ∗ Fionn D. Malone, † Miguel A. Morales, ‡ and David R. Reichman § Department of Chemistry, Columbia University, New York, NY 10027, USA. Quantum Simulations Group, Lawrence Livermore National Laboratory, 7000 East Avenue, Livermore, CA 94551 USA.
We explore the extended Koopmans’ theorem (EKT) within the phaseless auxiliary-field quantumMonte Carlo (AFQMC) method. The EKT allows for the direct calculation of electron addition andremoval spectral functions using reduced density matrices of the N -particle system, and avoids theneed for analytic continuation. The lowest level of EKT with AFQMC, called EKT1-AFQMC, isbenchmarked using small molecules, 14-electron and 54-electron uniform electron gas supercells,and diamond at the Γ-point. Via comparison with numerically exact results (when possible) andcoupled-cluster methods, we find that EKT1-AFQMC can reproduce the qualitative features ofspectral functions for Koopmans-like charge excitations with errors in peak locations of less than0.25 eV in a finite basis. We also note the numerical difficulties that arise in the EKT1-AFQMCeigenvalue problem, especially when back-propagated quantities are very noisy. We show how asystematic higher order EKT approach can correct errors in EKT1-based theories with respect tothe satellite region of the spectral function. Our work will be of use for the study of low-energycharge excitations and spectral functions in correlated molecules and solids where AFQMC can bereliably performed. I. INTRODUCTION
The dynamical response to external perturbation isone of the most powerful means of experimentally prob-ing molecules and materials. Examples include angle-resolved photoemission spectroscopy, electron energyloss spectroscopy, and inelastic neutron scattering, eachof which encodes the excitation spectrum of a many-bodysystem. The theoretical description of such experimentscan be modeled (in the linear response regime) by con-sidering many-body Green’s functions. For example,differential cross sections from direct and inverse pho-toemission experiments can be related to the (retarded)single-particle Green’s function. In a general sense, theseobservables are connected to spectral functions describ-ing electron removal and addition via the single-particleGreen’s function.
Given the above facts, the theoretical descriptionof dynamical response properties have been dominatedby Green’s function-based approaches mainly due tothe direct access to the spectral function that theyafford.
Among Green’s function methods, the G W approach is perhaps the most widely used approxima-tion to Hedin’s equations for the description of quasi-particle spectra in solids. This approach has providednumerous valuable insights for the interpretation of ex-perimental data.
Despite these successes, G W hasseveral deficiencies, such as a notable dependence on theinput Green’s function G and the screened Coulomboperator W , a poor description of satellite regionof the spectral function, and the absence of certainconserving properties. Attempts have been made to ad-dress these deficiencies via (partially) self-consistent GWapproaches, as well as with the incorporation of ver-tex corrections including via the use of cumulant-based approaches.
There has also been a sizable effort to construct spec-tral functions based on wavefunction methods. Theseapproaches include the use of matrix product states(MPS), algebraic diagrammatic constructions, and coupled-cluster Green’s function (or equation-of-motion coupled-cluster) methods.
These and relatedapproaches have distinct strengths and weaknesses interms of both cost and accuracy and continue to be ac-tively pursued.Another useful path to the description of spectral in-formation is based on projector quantum Monte Carlo(PQMC) approaches.
PQMC methods provide ahighly accurate means to simulate the ground state prop-erties of correlated solids. Unlike the aforementionedwavefunction-based approaches, PQMC methods do notprovide direct access to real-time and real-frequencyGreen’s functions. This is a direct consequence of theimaginary-time propagation at the heart of all PQMCapproaches. A popular way around this hurdle is to firstobtain the imaginary-time Green’s function and then per-form analytical continuation to obtain the real-frequencyGreen’s function.
Unfortunately, analytic continu-ation is numerically ill-conditioned, and the methodsto perform analytic continuation such as the maximumentropy method can exhibit difficulties in resolv-ing sharp features in the real-frequency spectral func-tion even if high quality imaginary-time Green’s func-tions are used as input.
Therefore, it is highly desir-able to develop an alternative means to obtain spectralfunctions which can work with PQMC methods withoutsacrificing its ground state accuracy. There have beenapproaches based on diffusion Monte Carlo (DMC) and the Krylov-projected full-configuration QMC (KP-FCIQMC) where one samples a low-energy Hamilto-nian matrix and solves an eigenvalue problem to obtaina low-energy spectrum. Similar to its ground state coun-terpart, the excited states from a DMC-based approach a r X i v : . [ phy s i c s . c h e m - ph ] J a n would be biased due to the fixed node error. On theother hand, KP-FCIQMC is numerically exact but scalesexponentially in system size analogously to its groundstate counterpart, FCIQMC. Therefore, the scope ofKP-FCIQMC has been limited to small systems.
Fur-thermore, because one is trying to obtain excited stateinformation from imaginary-time propagation in theseapproaches, the higher-lying excited states are exponen-tially harder to obtain. This makes it challenging forthese approaches to estimate high energy spectral infor-mation.The approach that we will examine in this work iscalled the extended Koopmans’ theorem (EKT).
The EKT generalizes the Koopmans’ theorem in Hartree-Fock (HF) theory for arbitrary many-body wavefunc-tions. Its working ingredients are reduced density ma-trices (RDMs) for an N -particle system and it pro-duces approximate ( N -1)-particle and ( N +1)-particle wavefunctions even without the N -particle ground statewavefunction. Due to this desirable feature, EKTmethods have been widely used as a means to ob-tain spectral information for approaches for which onehas access neither to real-time Green’s functions norwavefunctions. Examples include direct RDM-basedmethods, density matrix functional theory, natu-ral orbital functional methods, and second-orderGreen’s function methods. EKT has also been ex-plored with wavefunction methods such as configura-tion interaction methods, Møller-Plesset perturbationmethods, and coupled-cluster methods. It is also apromising way for any QMC method to compute excitedstate and spectral information if the necessary RDMscan be constructed. The EKT has also been used toobtain the quasi-particle band structure of silicon, ion-ization potentials and electron affinities of atoms, andthe Fermi velocity of graphene using VMC. Lastly, theEKT has been combined with DMC to study similarsystems. A PQMC approach that can be readily combinedwith the EKT is the phaseless auxiliary-field QMC (ph-AFQMC) method. ph-AFQMC has emerged as aflexible, accurate and scalable many-body method. Itimposes an approximate gauge boundary condition (i.e.,the phaseless constraint) on the imaginary-time evolu-tion of Slater determinant walkers, completely remov-ing the Fermionic phase problem. While the resultingenergy is neither exact nor a variational upper boundto the exact ground state energy, many benchmarkstudies have demonstrated the accuracy of ph-AFQMCand its related variants. Furthermore, with recentadvances in local energy evaluation techniques in ph-AFQMC, the cost for obtaining each statisticalsample scales cubically with system size, which rendersit less expensive than many other many-body methods.With the advent of the back-propagation (BP) methodin ph-AFQMC, with some additional effort onecan compute pure estimators for any operator, includ-ing those that do not commute with the Hamiltonian. Therefore, one can compute the relevant input for theEKT directly from ph-AFQMC using the BP algorithm.This is the direction we persue in this work.This paper is organized as follows. We first present thegeneral framework of the EKT, its most common formEKT-1, and its extension, EKT-3. We then discuss howto obtain the relevant input for EKT-1 using BP and ph-AFQMC. We assess the accuracy of EKT-1-AFQMC ona variety of small molecules and the uniform electron gasmodel. We further show the qualitative failure of EKT-1for the core spectra of the 14-electron uniform electrongas model and illustrate the drastic improvement uponthis result from using EKT-3 on the same model. We alsoapply EKT1-AFQMC to the 54-electron uniform electrongas model and diamond at the Γ-point. We conclude andsummarize our most important findings.
II. THEORYA. The Extended Koopmans’ Theorem
In order to compute quasi-particle gaps and spectralfunctions, one must compute ionization potential andelectron attachment energies along with the associatedwavefunctions (or at least squared amplitudes for spec-tral weights). While we focus in this work on electronremoval processes, we keep our presentation of theorygeneral so that it is also applicable to electron additionprocesses.In the EKT approach, we consider wavefunctions | Ψ N ± ν (cid:105) = ˆ O ± ν | Ψ N (cid:105) , (1)where the electron addition operator ˆ O + ν isˆ O + ν = (cid:88) p ( c + ) νp ˆ a † p (2)for 1-particle excitations, and the electron removal oper-ator ˆ O − ν is ˆ O − ν = (cid:88) p ( c − ) νp ˆ a p (3)for 1-hole excitations. We obtain the linear coefficients c ± by minimizing the following variational energy expres-sion:∆ E ± ν = E ( N ± ν − E ( N )0 = (cid:104) Ψ N | ( ˆ O ± ν ) † [ ˆ H , ˆ O ± ν ] | Ψ N (cid:105)(cid:104) Ψ N | ( ˆ O ± ν ) † ˆ O ± ν | Ψ N (cid:105) , (4)where we have defined E ( N ± ν = (cid:104) Ψ N | ( ˆ O ± ν ) † ˆ H ˆ O ± ν | Ψ N (cid:105)(cid:104) Ψ N | ( ˆ O ± ν ) † ˆ O ± ν | Ψ N (cid:105) , (5)and assumed ˆ H| Ψ N (cid:105) = E ( N )0 | Ψ N (cid:105) . (6)We refer this approach to as EKT-1. The excitation levelsin Eq. (3) and Eq. (2) can be systematically increased toachieve a greater accuracy in principle at the expense ofgreater computational costs. The next level of theorywould incorporate 2h1p and 2p1h excitations instead ofEq. (3) and Eq. (2), respectively:ˆ O + ν = (cid:88) pqr ( c + ) νpqr ˆ a r ˆ a † q ˆ a † p (7)and ˆ O − ν = (cid:88) pqr ( c − ) νpqr ˆ a † r ˆ a q ˆ a p (8)These operators include EKT-1 excitations because when r = q we recover the 1h and 1p excitations, as in Eq. (3)and Eq. (2). We refer this higher level of theory to asEKT-3.
1. EKT-1
We consider the following Lagrangian for 1h and 1pexcitations L [ c ν ] = (cid:104) Ψ N | ( ˆ O ± ν ) † [ ˆ H , ˆ O ± ν ] | Ψ N (cid:105) − (cid:15) ν ± (( c ν ) † S ± c ν − , (9)where S ± is a pertinent metric matrix for normalizationand (cid:15) ν ± is a Lagrange multiplier. We note that (cid:15) ν ± = ± ∆ E ± ν . (10)The normalization of | Ψ N ± ν (cid:105) is ensured by the constraintin Eq. (9). The stationary condition of Eq. (9) withrespect to ( c ν ± ) † then leads to a generalized eigenvalueequation, F ± c ν ± = (cid:15) ν ± S ± c ν ± (11)where the generalized Fock matrix is defined as (assum-ing that | Ψ N (cid:105) is normalized)( F − ) pq = (cid:104) Ψ N | ˆ a † p [ ˆ H , ˆ a q ] | Ψ N (cid:105) , (12)and ( F + ) pq = (cid:104) Ψ N | ˆ a p [ ˆ H , ˆ a † q ] | Ψ N (cid:105) , (13)and the corresponding metric matrix S ± is S − = P , (14)and S + = I − P T . (15)Here, P is the one-body reduced density matrix (1-RDM), P pq = (cid:104) Ψ N | ˆ a † p ˆ a q | Ψ N (cid:105) . (16) The electron attachment and ionization potential simplyfollow (cid:15) + = − EA and (cid:15) − = IP (assuming ν correspondsto the lowest energy state). Then the quasiparticle gap isgiven as ∆ E qp = (cid:15) + + (cid:15) − . We note that these Fock matri-ces are not Hermitian unless | Ψ N (cid:105) is an exact eigenstateof ˆ H .To provide more detailed expressions, let us define ageneric ab-initio Hamiltonian,ˆ H = ˆ H + ˆ H , (17)with ˆ H = (cid:88) pq h pq a † p ˆ a q , (18)ˆ H = 12 (cid:88) pqrs (cid:104) pq | rs (cid:105) ˆ a † p ˆ a † q ˆ a s ˆ a r , (19)where h pq is the one-body Hamiltonian matrix elementand (cid:104) pq | rs (cid:105) is the two-electron integral tensor in Dirac no-tation. Substituting Eq. (17) into Eq. (12) and Eq. (13),it can be shown that F ± can be evaluated with the 1-RDM and the two-body RDM (2-RDM):( F − ) pq = − (cid:88) q h qr ( S − ) pr + 12 (cid:88) trs (cid:104) tq || rs (cid:105) Γ rspt , (20)and ( F + ) pq = (cid:88) q h qr ( S + ) pr + 12 (cid:88) trs (cid:104) rt || qs (cid:105) Γ sprt + (cid:88) rs ( S − ) rs (cid:104) pr || qs (cid:105) , (21)where the 2-RDM Γ isΓ rspt = (cid:104) Ψ N | ˆ a † p ˆ a † t ˆ a s ˆ a r | Ψ N (cid:105) , (22)and the antisymmetrized two-electron integral tensor isdefined as (cid:104) pq || rs (cid:105) = (cid:104) pq | rs (cid:105) − (cid:104) pq | sr (cid:105) . (23)
2. EKT-3
For many solid state systems, including 1h or 1p ex-citations only is not sufficient as the restriction to suchexcitations would not be capable of describing satellitepeaks.
A straightforward way to obtain satellite peaksin addition to dominant quasiparticle peaks is to includehigher-order excitations in the EKT ansatz. Thus, EKT-3 is the next level in this hierarchy that can be attempted.Although EKT-3 has been mentioned in literature and approximately implemented (neglecting the oppo-site spin term) before, to the best of our knowledgethis work presents the first complete implementation ofEKT-3 along with numerical results.The corresponding generalized Fock operator for theIP problem reads( F − ) pqr,stu = (cid:104) Ψ N | ˆ a † p ˆ a † q ˆ a r [ ˆ H, ˆ a † u ˆ a t ˆ a s ] | Ψ N (cid:105) . (24)Using the SecondQuantizationAlgebra package devel-oped by Neuscamman and others, we derived acomplete spin-orbital equation of the generalized Fockoperator:( F − ) ijk,lmn = − h kn Γ mlij − (cid:88) a ( h la Γ amij δ kn + h ma δ kn Γ alij + h la Γ amkijn − h ma Γ alkijn + h na Γ mlkija )+ 12 (cid:88) ab ( (cid:104) lm || ab (cid:105) δ kn Γ baij − (cid:104) kl || ab (cid:105) Γ bamijn + (cid:104) km || ab (cid:105) Γ balijn − (cid:104) ka || nb (cid:105) Γ bmlija − (cid:104) lm || ab (cid:105) Γ bakijn )+ 12 (cid:88) abc ( (cid:104) ma || bc (cid:105) δ kn Γ cblija − (cid:104) la || bc (cid:105) δ kn Γ cbmija − (cid:104) la || bc (cid:105) Γ cbmkijna + (cid:104) ma || bc (cid:105) Γ cblkijna − (cid:104) na || bc (cid:105) Γ amlkijbc ) , (25)where the three-body RDM (3-RDM) isΓ npqijk = (cid:104) Ψ N | ˆ a † i ˆ a † j ˆ a † k ˆ a q ˆ a p ˆ a n | Ψ N (cid:105) , (26)and the four-body RDM (4-RDM) isΓ mnpqijkl = (cid:104) Ψ N | ˆ a † i ˆ a † j ˆ a † k ˆ a † l ˆ a q ˆ a p ˆ a n ˆ a m | Ψ N (cid:105) . (27)The pertinent metric, S , for this generalized eigenvalueproblem is S pqr,stu = (cid:104) Ψ N | ˆ a † p ˆ a † q ˆ a r ˆ a † u ˆ a t ˆ a s | Ψ N (cid:105) = δ ur Γ stpq − Γ strpqu . (28)The storage requirement of the 4-RDM scales as O ( N )and it becomes prohibitively expensive for more than 16orbitals. To circumvent this problem, we approximatethe 4-RDM via a cumulant expansion. The cumulantapproximation to the 4-RDM has been used in multi-reference perturbation theory and configuration interac-tion methods previously. In essence, the 4-RDMis approximately constructed from four classes of terms:(1) 1-RDM × × × × × × × Practical implementations may be achieved using spin-orbital expressions where we consider two spin-blocks in c − (( ααα ) and ( αββ )) for removing an α electron:ˆ O − ν = (cid:88) p α q α r α ( c − ) νp α q α r α ˆ a † r α ˆ a q α ˆ a p α + (cid:88) p α q β r β ( c − ) νp α q β r β ˆ a † r β ˆ a q β ˆ a p α . (29)Consequently, this leads to four distinct spin-blocks for F : ( αααααα ), ( αββααα ), ( ααααββ ), and ( αββαββ ). B. Spectral Functions from EKT
We write the retarded single-particle Green’s functionin a finite basis as iG Rpq ( t, t (cid:48) ) = θ ( t − t (cid:48) ) (cid:104) Ψ N |{ ˆ a p ( t ) , ˆ a † q ( t (cid:48) ) }| Ψ N (cid:105) , (30)where θ ( t ) is the Heaviside step function. Assuming theˆ H is time independent, we can write the Green’s functionin the frequency domain as G Rpq ( ω + iη ) = (cid:104) Ψ N | ˆ a p ω − ( ˆ H − E ( N )0 ) + iη ˆ a † q | Ψ N (cid:105) + (cid:104) Ψ N | ˆ a † q ω − ( ˆ H − E ( N )0 ) + iη ˆ a p | Ψ N (cid:105) , (31)where η is a small positive constant and E N is theground-state energy of N -particle system.The EKT approach offers a systematically improvableway to approximate the evaluation of Eq. (31). This isbecause one can form projection operators on the sub-space of EKT excitations,ˆ P ± = (cid:88) µν | Ψ N ± µ (cid:105) ( S ± µν ) − (cid:104) Ψ N ± ν | (32)where | Ψ N ± ν (cid:105) are approximate wavefunctions obtainedvia the EKT as defined in Eq. (1) and S ± is a metricin the pertinent space. Using Eq. (32), we obtain anapproximate G R , G Rpq ( ω + iη ) (cid:39) (cid:104) Ψ N | ˆ a p ˆ P + ω − ( ˆ H − E ( N )0 ) + iη ˆ P + ˆ a † q | Ψ N (cid:105) + (cid:104) Ψ N | ˆ a † q ˆ P − ω − ( ˆ H − E ( N )0 ) + iη ˆ P − ˆ a p | Ψ N (cid:105) (33)This approximation can be systematically improved ashigher-order excitations are included in Eq. (3) andEq. (2). It is exact when | Ψ N ± ν (cid:105) spans the entire( N ± P ± ˆ H ˆ P ± | Ψ N ± ν (cid:105) = E N ± ν | Ψ N ± ν (cid:105) (from Eq. (5)), we obtain the (approximate) Lehmannrepresentation of the Green’s function G Rpq ( ω ) = (cid:88) ν (cid:104) Ψ N | ˆ a p | ˜Ψ N +1 ν (cid:105)(cid:104) ˜Ψ N +1 ν | ˆ a † q | Ψ N (cid:105) ω − ( E ( N +1) ν − E ( N )0 ) + iη + (cid:88) ν (cid:104) Ψ N | ˆ a † q | ˜Ψ N − ν (cid:105)(cid:104) ˜Ψ N − ν | ˆ a p | Ψ N (cid:105) ω + ( E ( N − ν − E ( N )0 ) + iη , (34)where we orthogonalized the eigenvectors by | ˜Ψ N ± ν (cid:105) = (cid:88) µ | Ψ N ± µ (cid:105) ( S ± ) − / µν . (35)Details concerning the numerical implementation of thisorthogonalization procedure are given in Appendix A.A spectral function in a finite basis set can be com-puted from A pq ( ω ) = − π lim η → + Im (cid:2) G Rpq ( ω + iη ) (cid:3) = A >pq ( ω ) + A While the ph-AFQMC formalism has been presentedbefore in detail, we review the essence of the algorithmto provide a self-contained description. The imaginarypropagation is given as | Ψ (cid:105) ∝ lim τ →∞ exp (cid:16) − τ ˆ H (cid:17) | Φ (cid:105) = lim τ →∞ | Ψ( τ ) (cid:105) , (44)where τ is the imaginary time, | Ψ (cid:105) is the exact groundstate of a Hamiltonian ˆ H and | Φ (cid:105) is an initial start-ing wavefunction with a non-zero overlap with | Ψ (cid:105) . Weassume no special structure in the underlying Hamilto-nian and work with the generic ab-initio Hamiltonians ofEq. (17).In ph-AFQMC, this imaginary-time propagation isstochastically implemented. One discretizes the imag-inary time τ with a time step of ∆ τ such thatfor N time steps we have τ = N ∆ τ . Using theTrotter approximation and the Hubbard-Stratonovichtransformation, a single time step many-bodypropagator can be written in integral form,exp( − ∆ τ ˆ H ) = (cid:90) d N α x p ( x ) ˆ B (∆ τ, x ) , (45)where p ( x ) is the standard normal distribution, x is avector of N α auxiliary fields and ˆ B is defined asˆ B (∆ τ, x ) = e − ∆ τ ˆ H e −√ ∆ τ x · ˆ v e − ∆ τ ˆ H + O (∆ τ ) , (46)where ˆ v is defined fromˆ H = − N α (cid:88) α ˆ v α . (47)The computation of the integral in Eq. (45) is carried outvia Monte Carlo sampling where each walker samples aninstance of x .The global wavefunction is, with importance sampling,represented as a linear combination of walker wavefunc-tions: | Ψ( τ ) (cid:105) = (cid:88) i w i ( τ ) | ψ i ( τ ) (cid:105)(cid:104) Ψ T | ψ i ( τ ) (cid:105) , (48)where w i is the weight of the i -th walker, | ψ i ( τ ) (cid:105) is thesingle Slater determinant of the i -th walker, and | Ψ T (cid:105) is the trial wavefunction. At each time step, each walkersamples a set of x , forms ˆ B (∆ τ, x ), and updates its wave-function by applying ˆ B (∆ τ, x ) to it. Practical implemen-tations employ the so-called “optimal” force bias whichshifts the Gaussian distribution, ¯x i (∆ τ, τ ) = −√ ∆ τ (cid:104) Ψ T | ˆ v (cid:48) | ψ i ( τ ) (cid:105)(cid:104) Ψ T | ψ i ( τ ) (cid:105) . (49)With the optimal force bias, a single time step propaga-tion can be summarized with two equations w i ( τ + ∆ τ ) = I ph ( x i , ¯x i , τ, ∆ τ ) × w i ( τ ) , (50) | ψ i ( τ + ∆ τ ) (cid:105) = ˆ B (∆ τ, x i − ¯x i ) | ψ i ( τ ) (cid:105) , (51)where the phaseless importance function in hybrid formis defined as I ph ( x i , ¯x i , τ, ∆ τ ) = | I ( x i , ¯x i , τ, ∆ τ ) |× max(0 , cos( θ i ( τ ))) , (52)with I ( x i , ¯x i , τ, ∆ τ ) = S i ( τ, ∆ τ ) e x i · ¯x i − ¯x i · ¯x i / , (53)and S i ( τ, ∆ τ ) = (cid:104) Ψ T | ˆ B (∆ τ, x i − ¯x i ) | ψ i ( τ ) (cid:105)(cid:104) Ψ T | ψ i ( τ ) . (cid:105) , (54)With this specific walker update instruction, all walkerweights in Eq. (48) remain real and positive and therebyit completely eliminates the fermionic phases problem.In ph-AFQMC, the simplest way to evaluate the ex-pectation value of an operator ˆ O is using the mixed esti-mator (cid:104) ˆ O (cid:105) mixed := (cid:104) Ψ T | ˆ O | Ψ( τ ) (cid:105)(cid:104) Ψ T | Ψ( τ ) (cid:105) (55)= (cid:80) i w i ( τ ) (cid:104) Ψ T | ˆ O | ψ i ( τ ) (cid:105)(cid:104) Ψ T | ψ i ( τ ) (cid:105) (cid:80) i w i ( τ ) . (56)The mixed estimator is an unbiased estimator only foroperators that commute with the Hamiltonian. For op-erators that do not commute with the Hamiltonian, themixed estimator can introduce significant biases due to the approximate trial wavefunctions that can be prac-tically used. To overcome this do we use the back-propagation algorithm and write (cid:104) ˆ O (cid:105) ≈ lim κ →∞ (cid:104) Ψ T | e − κ ˆ H ˆ O | Ψ( τ ) (cid:105)(cid:104) Ψ T | e − κ ˆ H | Ψ( τ ) (cid:105) (57)= lim κ →∞ (cid:80) i w i ( τ + κ ) (cid:104) ψ i ( κ ) | ˆ O | ψ i ( τ ) (cid:105)(cid:104) ψ i ( κ ) | ψ i ( τ ) (cid:105) (cid:80) i w i ( τ + κ ) . (58)To summarize, we propagate | Ψ (cid:105) until κ + τ , storingthe walker wavefunction at time τ . We can then splitthe propagation into κ back-propagation and τ forward-propagation as in Eq. (58). The back propagated wave-function is constructed by applying a walker’s propa-gators to the trial wavefunction from the κ portion ofthe path. Practically the convergence of the expecta-tion value has to be monitored with respect to the backpropagation time κ . It should be emphasized that in ph-AFQMC the walker wavefunction is a single determinantwavefunction.It was found in Ref. 115 that the standard back-propagation algorithm described in Eq. (58) can yieldpoor results in ph-AFQMC when applied to ab-initio systems. The authors devised a number of additionalsteps to reduce the phaseless error, the most accurate ofwhich was to partially restore the phase and cosine fac-tors along the back propagation portion of the path. Inthis work we restore phases along the back propagatedpath as well as along the forward direction. Practicallythis amounts to storing the phases and cosine factors be-tween [ τ − κ, τ + κ ] and multiplying these by the weightsappearing in Eq. (58). This additional restoration ofpaths along the forward direction was not described inRef. 115 but was used in practice and we found itnecessary to obtain more accurate results for the systemsstudied here.In EKT1-AFQMC, we directly sample F ± using theback-propagated estimator form. This boils down to theevaluation of the 2-RDM appearing in Eq. (22) using theback propagated 1-RDM via Wick’s theorem:Γ rspt = P pr P ts − P ps P tr . (59)With these ingredients, we can evaluate F ± by contract-ing the one- and two-body matrix elements with the backpropagated 1- and 2-RDMs.While efficient implementations are not the focus ofour efforts in this work, we mention the computationalcost of producing one back-propagated sample of F ± and S ± . A sample of S ± has the same overhead as comput-ing a back-propagated 1-RDM sample which scales as O ( N M ) where M is the number of orbitals and N isthe number of occupied orbitals. The cost for producinga sample of F ± is more involved and depends on the in-tegral factorization that one chooses to use. Using themost common integral factorization, i.e., the Choleskyfactorization, it can be shown that the cost scales as O ( M X ) where X is the number of Cholesky vectors(also note that the 2-RDM is never explicitly formed).With tensor hypercontraction, the cost can bebrought down to overall cubic. If one were to just imple-ment a matrix-vector product for iterative eigensolvers,the Cholesky factorization can achieve cubic scaling permatrix-vector product as well. The cost is increased inEKT3-AFQMC where each matrix-vector product sam-ple costs O ( M ). It is potentially possible to reduce thiscost further by also factorizing EKT amplitudes in a THCformat, as is done in Ref. 128.We leave the exploration of EKT3-AFQMC for the fu-ture study and focus on EKT1-AFQMC in this work. D. Uniform Electron Gas Aside from small molecular benchmarks, we also studythe spectral properties of the uniform electron gas (UEG)model. The UEG model is usually defined in the plane-wave basis set, which gives the one-body operatorˆ H = (cid:88) K | K | a † K a K (60)and the electron-electron interaction operator is (in aspin-orbital basis)ˆ H = 12Ω (cid:88) K (cid:54) = , K , K π | K | a † K + K a † K − K a K a K , (61)where K here is a planewave vector and Ω is the vol-ume of the unit cell. In addition to ˆ H and ˆ H , thereis a constant term that arises due to a finite-size effect.Specifically, the Madelung energy E M should be includedto account for self-interactions associated with the Ewaldsum under periodic boundary conditions via E M = N ξ, (62)with ξ = − × . × (cid:18) π (cid:19) / N − / r − s , (63)where N is the number of electrons in the unit cell and r s is the Wigner-Seitz radius. We define the UEG Hamil-tonian as a sum of these three terms,ˆ H UEG = ˆ H + ˆ H + E M . (64)The Madelung constant can be either included in theHamiltonian as written in Eq. (64) or it can be includedas an a posteriori correction to the simulation done with-out it. When the latter choice is made, the spectral func-tions have to be shifted accordingly in order to comparethe results obtained from the former approach. The cor-responding shift can be derived from a shift in the polesIP = E ( N − − E ( N ) → − ξ , (65)EA = E ( N + 1) − E ( N ) → ξ . (66) Regardless of whether we included the Madelung con-stant in the Hamiltonian, there is an additional correc-tion of − ξ for both IP and EA to remove spurious imageinteractions coming from the excess charge created. Therefore, overall electron removal poles are shifted by − ξ and electron addition poles remain the same.For many molecular quantum chemistry methods, of-ten the two-electron integral tensor is assumed to be 8-fold symmetric. Practical implementations utilize thissymmetry to simplify equations as well. As such, with-out the 8-fold symmetry, molecular quantum chemistrymethods would not produce correct answers even thoughthe UEG Hamiltonian only contains real-valued matrixelements. This complication of the UEG Hamiltonianbecomes more obvious once we write it in the form ofEq. (19) with (cid:104) pq | rs (cid:105) = 1Ω 4 π | K p − K r | δ K p − K r , K q − K s × (1 − δ K p , K r ) . (67)The permutation between p and r or between q and s alters the value of the integral tensor because of the Kro-necker delta term. This is a direct consequence of us-ing a planewave basis which is complex, unlike the usualGaussian orbitals. To circumvent any complications dueto this, we perform a unitary transformation that rotatesthe planewave basis into a real-valued basis. Namely, forgiven K and − K (assuming K (cid:54) = ), we use U = 1 √ (cid:18) − i i (cid:19) . (68)We apply this transformation to every pair of K and − K in the two-electron integral tensor in Eq. (67). Theresulting transformed integral tensor now recovers thefull 8-fold symmetry. One can also transform observablessuch as spectral functions back to the original basis using U † when necessary. III. COMPUTATIONAL DETAILS All quantum chemistry calculations are performedwith PySCF which include mean-field (HF) calcula-tions, coupled-cluster with singles and doubles (CCSD)and CCSD with perturbative triples (CCSD(T)). All one-and two-electron integrals needed for ph-AFQMC werealso generated with PySCF . ph-AFQMC calculations weremostly performed with QMCPACK and PAUXY wasused to crosscheck some results.We use a selected configuration interaction (CI)method called heat-bath CI (HCI) to produce nu-merically exact IPs within a basis whenever possible.Furthermore, we compute the HCI IPs within EKT1 (i.e.,EKT1-HCI) using the 1- and 2-RDM from variationalHCI wavefunctions along with Eq. (20). Since there canbe an inherent bias of EKT1 itself, we provide EKT1-HCIas an “exact” result for IPs within the limits of the EKT1approach. This can be used to quantify the phaselessand back-propagation errors in EKT1-AFQMC. In HCI,there is a single tunable parameter, (cid:15) that controls thevariational energy which is used to select determinantsto be included in the variational expansion. We also usethe 3-RDM of the variational wavefunction to computethe 4-RDM via the cumulant construction. These3- and 4-RDMs are further used to construct the EKT3Fock matrix in Eq. (25). This approach is referred toas EKT3-HCI. The eigenvalues and eigenvectors of theEKT3 Fock matrix (with its pertinent metric) can thenbe used to produce IPs and spectral functions of EKT3.We tuned (cid:15) to be such that the resulting second-orderEpstein-Nesbet perturbation energy is no greater than 1m E h for every system except the UEG model. In theUEG model, we observed a PT2 correction of 3 m E h at r s = 4. This was found to be sufficient to produceaccurate EKT IPs for systems studied here. All calcu-lations are performed with a locally modified version of Dice . We used a timestep of 0.01 au and the pair branch pop-ulation control method used the hybrid propagationscheme for all ph-AFQMC simulations. For the smallatoms and molecules and 14 electron UEG examples, weused 2880 walkers while for the 54 electron UEG we used1152 walkers. All calculations used restricted Hartree–Fock (RHF) trial wavefunctions except for the chargedspecies where we instead used unrestricted Hartree–Fockwavefunctions (UHF). The exception to this was CH where we used the same RHF orbitals for the chargedspecies as we found that the UHF solution broke spa-tial symmetry and led to a large phaseless constraintbias. All AFQMC results are performed with the phase-less approximation so we simply refer ph-AFQMC to asAFQMC in the following sections.We adapted the standard dynamical Lanczosalgorithm to obtain spectral functions of Eq. (64).Even by exploiting symmetry and using a distributedsparse Hamiltonian, dynamical Lanczos results couldonly be obtained for the smallest UEG system of14 electrons in 19 plane waves, corresponding to an N -electron Hilbert space size of 2 . × determinants.In the dynamical Lanczos algorithm, one first obtainsthe N -particle ground state, | Ψ N (cid:105) iterating withinthe Lanczos Krylov subspace. Ultimately, our goal isto compute Eq. (36) which then requires another runof the Lanczos algorithm. For an electron removalproblem, we pick an orbital index i and generate aninitial vector in the ( N − | f (cid:105) =ˆ a i | Ψ N (cid:105) / (cid:104) Ψ N | ˆ a † i ˆ a i | Ψ N (cid:105) . Each Lanczos iteration thengenerates coefficients, { a k } and { b k } , for the followingcontinued fraction expression: A ii ( ω, η ) = − π Im (cid:104) Ψ N | ˆ a † i ˆ a i | Ψ N (cid:105) z − a − b z − a − b z − a ··· , (69)where z = E ( N )0 − ω + iη with some spectral broadeningconstant η . We take a total of 50 Lanczos iterations to Experiment ∆HCIHe 24.59 -0.23Be 9.32 -0.03Ne 21.56 -0.13FH 16.19 -0.12N O 12.62 -0.07TABLE I. Experimental first ionization potentials (eV) anddeviation (eV) of the numerically exact ∆HCI from these re-sults for chemical species considered. ∆HCI employed theaug-cc-pVDZ basis. generate the continued fraction coefficients and this wasenough to converge the low-energy spectrum within theenergy scale that is relevant in this work. IV. RESULTS AND DISCUSSION While the EKT approach is valid for both electron re-moval and electron addition, for numerical results wefocus on electron removal processes (IP energies andelectron removal spectral functions) for simplicity. Webenchmark the IP energies from the proposed EKT1-AFQMC approach over several small chemical systemsand the UEG model. Furthermore, we also show promis-ing improvements over EKT1 using EKT3-HCI whensatellite peaks are important. We use a “∆ method”to denote a scheme where we run the pertinent methodfor both N - and ( N − A. Small chemical systems in the aug-cc-pVDZbasis In this section, we study seven small chemical systems(He, Be, Ne, FH, N , CH , and H O) that have well-documented experimental IPs. We use the nomencla-ture FH for the hydrogen fluoride molecule to distinguishit from the abbreviation for Hartree-Fock (HF). All ge-ometries were taken from ref. 141. We used a relativelysmall basis set, namely aug-cc-pVDZ, to obtain goodstatistics in back-propagated estimators. We have usedmore than 3000 back-propagated estimator samples inall cases considered here, each of which requires a back-propagation time of greater than 4 a.u. This results in atotal propagation time longer than 12000 a.u., which isunusually long for standard AFQMC calculations. Theuse of this basis set also allows for a direct comparisonbetween AFQMC and numerically exact HCI within thisbasis set.The goal of this numerical section is to quantify thethree sources of error in addition to the basis set incom-pleteness error in EKT1-AFQMC based on simple exam-ples where exact simulations are possible. These three Ground state IPHe 0.00 0.00Be 0.01 -0.01Ne -0.03 0.07FH -0.03 0.07N -0.04 0.06CH -0.03 0.06H O -0.03 0.04TABLE II. Error (eV) in the AFQMC N -electron systemground state energy and in the first ionization potential withrespect to the corresponding HCI results. The statistical errorbar of AFQMC is less than 0.01 eV and therefore we do notpresent them here. sources of error are:1. Phaseless constraint errors. As mentioned, thephaseless constraint is necessary to remove thephase problem that arises in the imaginary-timepropagation. However, due to this constraint, theresulting ground state energies and properties (e.g.,RDMs) are biased.2. Back-propagation errors. The back-propagation al-gorithm incurs additional errors. This was notedand studied in detail in ref. 115. For instance, inref. 115, it was shown that for neon the phaselesserror with a simple trial wavefunction is negligible(below 1 m E h ) but the error in the one-body en-ergy from the back-propagated 1-RDM was about5 m E h .3. EKT1 errors. While systematically improvablewith higher-order excitations, EKT1 is not an exactapproach to quasiparticle spectra unless all ordersof excitations are included. Nonetheless, for thefirst IP, it has been numerically and analyticallysuggested that EKT1 approaches the exact IP inthe basis set limit if the exact 1- and 2-RDMs areused. Beyond the first IP, we will showthat EKT1 qualitatively fails to capture satellitepeaks that arise in the case of the core spectrum ofthe UEG model.In Table I, we present numerically exact first IPs ofmolecules within this basis set using ∆HCI. The basis setincompleteness error can be as large as 0.3 eV in thesemolecules and therefore we will only compare AFQMCresults to these numerically exact results in the same ba-sis set as opposed to comparing to the experimental data.We do not expect the qualitative conclusions of our studyto change with larger basis sets.Next, we assess the phaseless bias in the N -electronsystem ground state energy as well as the error in the∆AFQMC IP energies compared to the ∆HCI IPs. InTable II, we present numerical data that detail the phase-less bias in these quantities. The ground state energyerror is less than 0.04 eV which is in the neighborhood ofthe usual standard of accuracy, 1 m E h . Unfortunately, EKT1-HCI KT EKT1-AFQMCHe 24.36 0.60 0.02Be 9.29 -0.87 0.06Ne 21.48 1.73 -0.04FH 16.13 1.58 0.01N O 12.60 1.27 -0.04TABLE III. Error (eV) in the first IP obtained from EKT1-AFQMC and the Koopmans’ theorem (KT) relative to EKT1-HCI. The statistical error bar of EKT1-AFQMC cannot beestimated without bias (see main text for discussion). in many cases ( N − N -electron systems, albeit withan opposite sign of the error. Thus, AFQMC does notbenefit from a cancellation of errors for the IP energy,which results in IP errors that are larger than those forthe ground state energy. The largest IP error we find isaround 0.07 eV.We present the EKT1 results in Table III. We referreaders to Appendix A where a detailed description ofthe EKT1 calculation is given. Theoretically “exact”EKT1 results can be obtained by using exact RDMs fromHCI. To quantify the back-propagation error of AFQMC(given that the phaseless error is very small for these sys-tems), we shall compare EKT1-AFQMC to EKT1-HCI.We also computed simple Koopmans’ theorem (KT) IPsusing HF and report these results in Table III. The er-ror of EKT1-AFQMC is small for most chemical species,but it becomes as large as 0.19 eV for CH . Even thoughthe phaseless bias in the ground state energy was foundto be very small (0.04 eV or less), EKT1-AFQMC errorsfrom using the back-propagated 1-RDM and the EKTFock matrix can be five times larger. Therefore, we at-tribute this error mainly to the back-propagation error.Nonetheless, the comparison against the simpler KT sug-gests that EKT1-AFQMC can readily recover the corre-lation contribution to quasiparticle energies with errorsless than 0.2 eV in these examples.We note that EKT1-AFQMC IPs do not have any sta-tistical error bars. This is not because these numbers aredeterministic but arises because EKT eigenvalues are notunbiased estimators. This was also observed in a similarapproach called Krylov-projected full configuration inter-action QMC. Interested readers are referred to Ref. 70 fordetails. In essence, eigenvalues of a noisy matrix whereeach element is normally distributed are not normallydistributed. Therefore, statistical error bars are difficultto estimate and are not simply associated with variancesof Gaussian distributions. Given the small statistical er-ror bars (on the order of 3 × − or less) on the diagonalelements of the EKT1 Fock matrix and 1-RDM, we ex-pect that these results are reproducible up to 0.01 eV ifone follows exactly the same numerical protocol given inAppendix A.Finally, we discuss the inherent error of EKT1 by com-0 EKT1-HCI EKT1-AFQMC EOM-IP-CCSDHe 0.00 0.02 0.00Be 0.00 0.06 0.00Ne 0.05 0.02 -0.28FH 0.06 0.07 -0.22N -0.15 0.04 -0.02H O 0.05 0.01 -0.15TABLE IV. Error (eV) in the first IP obtained from EKT1-HCI, EKT1-AFQMC, and EOM-IP-CCSD relative to ∆HCI(given in Table I) in the aug-cc-pVDZ basis. The statisti-cal error bar of EKT1-AFQMC cannot be estimated withoutbiases (see main text for discussion). paring EKT1-HCI and EKT1-AFQMC against ∆HCI asin Table IV. We also compare a more widely used ap-proach called equation-of-motion coupled cluster ioniza-tion potential with singles and doubles (EOM-IP-CCSD)to gauge the magnitude of EKT1 errors. Since this isa benchmark on the first IP of molecules, EKT1-HCIis expected to be quite accurate. This expectation isdue to the general belief that EKT1 with exact RDMsyields exact first IP. Given this, the small errorsof EKT1-HCI in Table IV are not surprising. The onlyexception is CH with an error of -0.15 eV, which webelieve will still approach the exact IP in the completebasis set limit. EKT1-AFQMC appears to be as good asEKT1-HCI, with an outlier for the case of N . The errorof EKT1-AFQMC relative to EKT1-HCI is comparableto the error of EKT1-HCI relative to ∆HCI in this ba-sis set. EOM-IP-CCSD generally does not work as wellas the EKT1 approaches with a maximum error of -0.28eV on the neon atom. While these finite basis set com-parisons are informative, we emphasize that more faircomparisons should be conducted in the complete basisset limit and we hope to carry these out in the future.Regardless, the EKT1-AFQMC results are encouraging.The main motivation for performing EKT1 withinAFQMC was to obtain spectral functions. Poles alonecan be obtained using the ground state AFQMC al-gorithm (i.e., ∆AFQMC) by imposing a proper con-straint for excited state descriptions. While the choiceof proper trials is a challenge for this purpose, suchan approach avoids the complications due to back-propagation. However, spectral weights cannot be ob-tained from ∆AFQMC. We show EKT1-AFQMC spec-tral functions for FH and N in Fig. 1. Based on the re-sults shown in Table III, EKT1-AFQMC is in good agree-ment with EKT1-HCI for FH, but not for N . There-fore, comparing these two cases is useful for understand-ing how back-propagation errors are reflected in spectralfunctions. In the case of FH, we do not see any visi-ble differences between EKT1-HCI and EKT1-AFQMCon the plotted energy scale. However, for N we canclearly see some deviation between EKT1-AFQMC andEKT1-HCI. Nonetheless, the main features of the spec-tral function are reproduced, namely three large quasi- particle peaks with the middle peak being the largest.We note that there are peaks with very small spectralweights in EKT1-HCI close to -30 eV, -27 eV, and -5 eVand that these features these do not appear in EKT1-AFQMC. We attribute this to a relatively large lineardependency cutoff (10 − ) needed in EKT1-AFQMC tostabilize the generalized eigenvalue problem as explainedin Appendix A. In both molecules, there are insignificantdifferences between the EKT1 and HF spectra in termsof the peak heights, locations and the number of peakswith a brodening parameter of η = 0 . B. The uniform electron gas (UEG) model AFQMC has emerged as a unique tool for simulatingcorrelated solids. A model solid that de-scribes the basic physics of metallic systems is the UEGmodel. The accuracy and scope of AFQMC in study-ing the UEG model has been well documented at zerotemperature and finite temperature. Motivated bythese studies, we investigate the spectral properties ofthe UEG model within the EKT approaches (EKT1 andEKT3). 1. 14-electrons/19-planewaves The first example that we consider is a relatively smallUEG supercell with only 19 planewaves. It is far from thebasis set limit as well as from the thermodynamic limit.However, it is small enough for one to produce unbiasedEKT results using HCI and numerically exact dynamicalLanczos results. We note that the spectral function ofthis benchmark UEG model at r s = 4 was first presentedin Ref. 50. We produced around 3000 back-propagationsamples, which yielded the largest statistical error in theFock matrix and 1-RDM on the order of 5 × − .We consider two Wigner-Seitz radii, r s = 0 . r s = 4. Based on our previous benchmark study ofAFQMC on this system, we expect that the phaselesserror in the ground state at r s = 0 . r s = 4. Com-pared to the numerically exact energies, it was foundthat the constraint bias in AFQMC is only -0.0118(6)eV at r s = 0 . r s = 4. Given thissmall ground state bias, we expect the EKT1-AFQMCapproach would be as accurate as EKT1-HCI if goodstatistics and accuracy in the back-propagated estima-tors can be achieved.In Fig. 2, we present spectral functions of this modelat r s = 0 . K = [0 , , However, such features are very simple in asmall supercell like this one. As can be seen in Fig. 2(a),1 30 20 10 0Energy (eV)0.0000.0050.0100.0150.0200.0250.0300.0350.040 D e n s it y o f s t a t e s ( e V ) (a) FH EKT1-HCIEKT1-AFQMCHF 30 20 10 0Energy (eV)0.0000.0050.0100.0150.0200.0250.0300.035 (b) N FIG. 1. Electron removal spectral functions from various methods (EKT1-HCI, EKT1-AFQMC, and HF) of (a) FH and(b) N in the aug-cc-pVDZ basis set. Note that in (a) EKT1-AFQMC is right on top of EKT1-HCI on the plotted scale. Abroadening parameter η = 0 . 120 100 80 60 40 20Energy (eV)0.00.10.20.30.4 D e n s it y o f s t a t e s ( e V ) (a) K = [0, 0, 0] LanczosEKT1-HCIEKT1-AFQMCHFEKT3-HCI 75 80 85 90 95 100Energy (eV)0.00.10.20.30.4 (b) K = [0, 0, 2 / L ] FIG. 2. Electron removal spectral functions of the 14-electron in 19-planewave UEG model from various methods at r s = 0 . K = [0 , , 0] and (b) K = [0 , , π/L ]. In (a), note that EKT1-AFQMC, EKT1-HCI, and HF are right on top of each other.In (b), EKT1-AFQMC and EKT1-HCI are right on top of each other. In both (a) an (b), Lanczos and EKT3-HCI are righton top of each other. A broadening parameter of 0.2 eV was used for all plots. we only have two peaks from the Lanczos method. Nei-ther of these peaks corresponds to a single Koopmans-likeexcitation. Namely, they cannot be found from a simplesingle ionization process that either HF (i.e., Koopmanstheorem) or EKT1 describes well. As a result of this,HF, EKT1-AFQMC, and EKT-HCI all fail to capture the peak split and only yield a single peak. The cor-relation effect is very marginal in the sense that nearlyno improvement was observed with the EKT1 methodscompared to HF.A significant improvement can be made by incorporat-ing higher-order terms such as 2h1p excitations. In other2 D e n s it y o f s t a t e s ( e V ) (a) K = [0, 0, 0] LanczosEKT1-HCIEKT1-AFQMCHFEKT3-HCI (b) K = [0, 0, 2 / L ] FIG. 3. Electron removal spectral functions of the 14-electron in 19-planewave UEG model from various methods at r s = 4 . K = [0 , , 0] and (b) K = [0 , , π/L ]. In (a), note that EKT1-AFQMC, EKT1-HCI, and HF are right on top of each other.In (b), EKT1-AFQMC and EKT1-HCI are nearly on top of each other. A broadening parameter of 0.2 eV was used for allplots. words, core ionization satellite states (up to leading or-der) require excitations such as( (cid:88) K (cid:48) (cid:54) = C K (cid:48) ˆ a † K (cid:48) ˆ a K (cid:48) )ˆ a K = (70)where C K (cid:48) are the excitation amplitudes. All of these ex-citations are included in EKT3. EKT3-HCI can nearlycompletely reproduce the exact spectral function despitethe use of the cumulant approximation for the 4-RDM.The cumulant approximation error is small especially inweakly correlated cases such as r s = 0 . K = [0 , , π/L ] corresponds to the top of thevalence band corresponding to the first IP of the UEGmodel. There is no satellite peak visible at this momen-tum and the single peak found from the EKT1 methodsis reasonable. The peak location of HF is displaced byabout +0.9 eV from the correct location while the EKT1methods yield a peak shifted by about -0.6 eV. We notethat EKT1-AFQMC is practically indistinguishable fromEKT1-HCI for both momenta which indicates a smallphaseless bias, a small back-propagation error, and goodstatistics for the estimators.A similar conclusion can be drawn for r s = 4 as shownin Fig. 3. For the core excitation spectrum in Fig. 3(a),we observe the same split peak structure observed at r s = 0 . 5. We see another smaller peak emerging on theleft shoulder of the peak near -12 eV. While EKT3-HCIis no longer exact, it reproduces most of the features in the exact spectral function including the emergence ofthe third peak. While EKT1-HCI and EKT1-AFQMCagree well, there is no visible improvement over HF. TheEKT1 methods all yield a single peak which is qualita-tively wrong. The valence excitation structure illustratedin Fig. 3(b) is relatively featureless, but there are smallpeaks emerging in the high energy (more negative) re-gion of the spectrum. EKT3-HCI shows good agreementwith Lanczos for the main quasiparticle peak and alsoproduces satellite features. There is a visible improve-ment of EKT1 approaches (with a deviation of the peakenergy of approximately -0.25 eV) compared to HF (withan approximate deviation of +0.34 eV). A slight devia-tion of EKT1-AFQMC from EKT1-HCI is observed, butthe difference in the main quasiparticle peak location isonly about 0.01 eV.Overall, in this small benchmark study, the EKT1 ap-proaches provide some improvement over HF for valenceexcitations and qualitatively fails for the core region.The agreement between the EKT1 valence peaks and theLanczos peaks is not perfect, with an error of about -0.6eV for r s = 0 . r s = 4. However, we em-phasize again that we expect EKT1 to become exact forthe first IP as the complete basis set limit is approached.Finally, EKT1-AFQMC is able to reproduce EKT1-HCInearly perfectly even for r s = 4, where the phaseless errorin the ground state energy is about 0.185(2) eV.3 D e n s it y o f s t a t e s ( e V ) EKT1-AFQMCHFAFQMCEOM-IP-CCSD FIG. 4. Electron removal spectral functions of the 54-electron in 287-planewave UEG model at r s = 2 from EKT1-AFQMC and HF. The first IPs from ∆AFQMC and EOM-IP-CCSD are shown for comparisons. A broadening parameterof 0.2 eV was used for all plots. 2. 54-electrons/287-planewaves Next, we consider a larger UEG supercell (54 elec-trons in 287 planewaves) where obtaining many back-propagation samples is difficult. We study r s = 2,where AFQMC can be reliably extrapolated to the basisset limit. We produced 600 back-propagation sampleswith a back propagation time of 8 a.u. This amounts toa total of 4800 a.u. propagation time, which is a longpropagation for this system size. Our approach yielded amaximum statistical error in 1-RDM and Fock matricesof 4 × − . While this error is not small, the proceduredescribed in Appendix A was enough to stabilize the fi-nal results. Unlike previous cases, we take the uppertriangular part of the Fock matrix and explicitly sym-metrize the Fock matrix. A linear dependency cutoff of10 − was used in EKT1-AFQMC. It is difficult to gen-erate highly accurate 1- and 2-RDMs from HCI for thissystem size, so for this system we do not have an exactbenchmark reference to compare to our EKT1-AFQMCresults. Similarly, EKT3-HCI is also intractable for thissystem size. Instead, we have performed EOM-IP-CCSDand ∆AFQMC to compare with and to gauge the mag-nitude of the errors of EKT1-AFQMC.In Fig. 4, EKT1-AFQMC and HF spectral functionsare shown. As expected, EKT1-AFQMC does not showany satellite peaks at all and EKT1-AFQMC only in-troduces a shift to the HF spectrum. Correlation ef-fects in the peak height do not appear large, but thepeak location changes by about 1 eV going from HF toEKT1-AFQMC. We also produced ∆AFQMC and EOM-IP-CCSD for comparison. The ∆AFQMC IP is 2.18(1)eV, EOM-IP-CCSD yields 1.91 eV, and EKT1-AFQMCgives 2.51 eV in this basis set. The deviation of EKT1-AFQMC from ∆AFQMC is 0.33(1) eV whereas EOM-IP-CCSD deviates by 0.27(1) eV. These deviations are sim-ilar in magnitude but with opposite signs. The accuracyof EOM-IP-CCSD is unclear because the ground state en- 20 15 10 5 0 5 10 15Energy (eV)0.000.020.040.06 D e n s it y o f s t a t e s ( e V ) EKT1-HCIEKT1-AFQMCHFEOM-IP-CCSD FIG. 5. EKT1-HCI and EKT1-AFQMC electron removalspectral functions of diamond at the Γ-point. A broadeningparameter η = 0 . ergy found from CCSD is higher than that of AFQMC by0.0291(1) eV per electron. In the basis set limit, AFQMCwas found to be as accurate as state-of-the-art diffusionMonte Carlo for the ground state energy, differing by0.0088(9) eV per electron or less. This suggests thatthe ground state correlation energy of CCSD may be onthe order of 0.03 eV per electron. How much of this erroris propagated to the EOM-IP-CCSD calculation remainsunclear. In the complete basis set limit, we believe thatthe first IP from EKT1-AFQMC will become more accu-rate and closer to that of ∆AFQMC. A more completecomparison should be conducted in this limit, but suchcalculations are very difficult due to the need to procuremany back-propagation samples in EKT1-AFQMC. C. Towards ab-initio solids: diamond at the Γ -point With recent advances in open-source software suchas PySCF , performing calculations on ab-initio solidsis relatively straightforward. While an implementa-tion of AFQMC with k -points has been previouslypresented, we only present a Γ-point result in thiswork. This is mainly because our current EKT1 im-plementation does not explicitly consider k -points. Wechose to study diamond because it is one of the simplestsolids just with two carbon atoms in its unit cell. Weused the GTH-PADE pseudopotential and the GTH-DZVP basis set. This system is overall as small as the smaller sys-tems considered in this work. Therefore, we could obtainover 6000 back-propagation samples with a 4 a.u. back-propagation time. Even with these many samples, thelargest error bar in the Fock and 1-RDM matrices was6 × − . Due to this large statistical error in the matrixelements, we used a linear dependency cutoff of 0.01 andsymmetrized the Fock matrix by taking only the uppertriangular part of it (see Appendix A for details). Wealso included the shift by the Madelung constant in thespectral functions.4As shown in Fig. 5, EKT1-AFQMC successfully repro-duces EKT1-HCI. A small quasiparticle peak near -5 eVis not accurately captured by EKT1-AFQMC largely dueto statistical noise. Small peaks are difficult to resolvein EKT1-AFQMC without reducing the statistical errorson each matrix element further. Nonetheless, two otherlarger quasiparticle peaks are well represented. Bothpeaks are well reproduced within 0.25 eV from EKT1-HCI. The improvement over HF is about 1 eV or so andthe first IP from EOM-IP-CCSD is within 0.02 eV of theEKT1-HCI result. It will be interesting to revisit theassessment of EKT1-AFQMC in the basis set and ther-modynamic limits via direct comparisons to experiments. V. CONCLUSIONS In this work, we have explored the extended Koop-mans’ theorem (EKT) approach to the computationof spectral functions via phaseless auxiliary-field quan-tum Monte Carlo (AFQMC). Previous attempts inAFQMC to obtaining spectral functions have resortedto analytic continuation which has well-documenteddrawbacks. The EKT approach is attractive be-cause it requires neither an explicit representation of theground state wavefunction nor analytic continuation tocompute spectral functions. Instead, its only inputs areN-particle reduced density matrices (N-RDMs) which canbe computed in AFQMC via the back-propagation algo-rithm. The motivation of our work was thus to use theEKT approach with the aim of avoiding numerical prob-lems arising in analytic continuation for the accurate as-sessment of real-frequency spectral information. Whilemany studies have so far focused on the simplest level ofthe EKT, the EKT approach is systematically improv-able with increasing order of excitations: 1h, 1p2h, etc.for electron removal and 1p, 1h2p, etc. for electron addi-tion. We presented the implementation of EKT1 (1h or1p) and EKT3 (1p2h or 1h2p). For EKT3, we proposedthe use of a cumulant approximation to the 4-RDM toavoid the steep storage requirements.We produced preliminary results using EKT1 withinAFQMC (EKT1-AFQMC) for small molecular systems,uniform electron gas (UEG) 14-electron and 54-electronsupercells, and diamond at the Γ-point. We focused onstudying the first ionization potential (IP) and electronremoval spectral functions of these systems. By compar-ing numerically exact EKT1 results based on heat-bathconfiguration interaction (i.e., EKT1-HCI), we showedthat despite statistical noise, EKT1-AFQMC can cap-ture most qualitative features of EKT1-HCI. We providea more detailed summary on our findings as follows:1. In small molecular benchmarks within the aug-cc-pVDZ basis, we found the maximum deviation ofEKT1-AFQMC from EKT1-HCI in the first IPto be 0.19 eV. These molecules have quite smallphaseless biases in the ground state energy ( ≤ . 04 eV) so we attributed additional biases to back-propagation. Electron removal spectral functionsfrom EKT1-AFQMC look qualitatively similar tothat of EKT1-HCI even in the least accurate case(N ).2. For the 14-electron UEG supercell (19-planewave)benchmark, we observed a qualitative failure ofEKT1 due to its inability to describe satellite statesat K = . We showed that EKT3 (within HCI)significantly improves this. Despite these failuresof EKT1, we found EKT1-AFQMC to have peaklocations that are nearly identical (within 0.01 eV)to EKT1-HCI for both r s = 0 . r s = 4. Giventhe noticeable phaseless bias at r s = 4, this resultis quite encouraging. Lastly, for the valence regionof the electron removal spectral function, we ob-served reasonable accuracy of EKT1 compared tothe exact spectral function. The location of thefirst IP was off by 0.4 eV for r s = 0 . r s = 4 . 0, which we expect to improve in largerbases.3. For the 54-electron UEG supercell (257-planewave)benchmark, we could not obtain EKT1-HCI duecomputational expense. Therefore, we attemptedto assess the accuracy of EKT1-AFQMC by com-paring the first IP of EKT1-AFQMC with thatof equation-of-motion IP coupled-cluster with sin-gles and doubles (EOM-IP-CCSD) and ∆-AFQMC.However, all three methods differ from each otherby more than 0.25 eV and a more thorough bench-mark in the basis set limit is highly desirable.4. For diamond at the Γ-point, EKT1-AFQMC pro-duced a qualitatively correct electron removal spec-tral function which agrees well with EKT1-HCI.However, EKT1-AFQMC peak locations were offby 0.2 eV from those of EKT1-HCI. We also notedthat EKT1-HCI first IP agrees with that of EOM-IP-CCSD within 0.01 eV.While a more extensive benchmark study is highly de-sirable, we cautiously conclude that EKT1-AFQMC isuseful for charge excitations that are heavily dominatedby Koopmans-like excitations. EKT1-AFQMC errors inpeak locations can be as large as 0.25 eV compared toEKT1-HCI, but the line shapes of EKT1-AFQMC closelyfollow those of EKT1-HCI in all systems considered inthis work.The greatest challenge of EKT1-AFQMC is currentlythe statistical inefficiency in obtaining relevant back-propagated quantities with error bars small enough toenable the construction of stable EKT1-AFQMC spectralfunctions. Future work must first be dedicated to improv-ing the statistical efficiency of back-propagation. Fur-thermore, better back-propagation algorithms are neededto reduce the back-propagation error further. A prac-tical implementation of EKT3-AFQMC using an itera-tive eigenvalue solver will be an interesting topic to ex-5plore in the future. Several interesting extensions areimmediately possible. First, extending the EKT frame-work to neutral excitations is relatively straightfor-ward and could be interesting to explore. Next, theextension of the EKT framework for finite-temperaturecoupled electron-phonon problems would provide a wayto compute temperature-dependent vibronic spectra di-rectly from AFQMC. We also leave the compari-son of these EKT-based spectral functions to analyticallycontinued spectral functions for a future study. VI. ACKNOWLEDGEMENTS We thank Sandeep Sharma for providing access to aversion of Dice with the 3-RDM and 4-RDM capabil-ity which was used in testing and producing EKT3-HCIresults presented in this work. We also thank GarnetChan for discussion about finite-size corrections to ion-ization energies. DRR acknowledges support of NSFCHE- 1954791. The work of FDM and MAM was per-formed under the auspices of the U.S. Department ofEnergy (DOE) by LLNL under Contract No. DE-AC52-07NA27344 and was supported by the U.S. DOE, Officeof Science, Basic Energy Sciences, Materials Sciences andEngineering Division, as part of the Computational Ma-terials Sciences Program and Center for Predictive Sim-ulation of Functional Materials (CPSFM). Computingsupport for this work came from the LLNL InstitutionalComputing Grand Challenge program. Appendix A: Numerical Details of EKT Determining the eigenvalues and eigenvectors ofEq. (11) from noisy QMC density matrices is non-trivial.We first diagonalise the metric matrices S ± , S ± = U Λ ± U † , (A1)and discard eigenvalues below a given threshold. In thiswork, unless noted otherwise, in the main text we used10 − for all EKT1-AFQMC results, which is on the orderof the largest statistical error in the EKT1 Fock matrix.We next construct the transformation matrix X ± = UΛ − / ± . (A2) This procedure is often referred to as canonical orthogo-nalization in quantum chemistry. Then, we solve˜ F ± ˜ c ν ± = (cid:15) ν ± ˜ c ν ± , (A3)where ˜ F ± = X ±† F ± X ± . (A4)Finally the eigenvectors in the original basis can be de-termined from ( c ν ± ) p = (cid:88) I ( X ± ) pI ( ˜ c ν ± ) I . (A5)Following Kent et al. , we explicitly zero out all matrixelements whose magnitude is smaller than two times thecorresponding statistical error bar. We explicitly sym-metrize RDMs, but leave the Fock matrix asymmetricas required for approximate wavefunctions. However, forthe more difficult problems considered in this work, suchas the 54-electron UEG electron and diamond, we foundthat symmetrizing the Fock matrix is useful so we chooseto symmetrize the Fock matrix in such cases (by takingonly the upper triangular part of the Fock matrix). Thesesteps improved the numerical stability of the eigenvalueproblem.We also note that there are generic numerical issuesarising in EKT even without any statistical sampling er-ror. This was observed in both EKT1-HCI and EKT3-HCI, where spurious solutions with large negative IPsappear. These states stem from the fact that the metricmatrix in EKT problems are generally low-rank, as theyare related to RDMs. For instance, in the EKT1 for-mulation, the metric we diagonalize for the IP problemis a 1-RDM whose rank is not so much larger than thenumber of electrons in the system. These spurious statescan be removed with larger cutoffs while often affectingpeak locations of quasiparticle states. Interestingly, thesespurious states all carry negligible spectral weights anddo not appear in spectra. Motivated by this observation,the most satisfying solution we found was to use as asmall threshold as possible and to measure the overlapbetween Koopmans states and eigenvectors to identifyquasiparticle-like eigenvectors. This was enough to iden-tify physical IP excitations that are quasi-particle-like.The same principle is applicable to EKT3-HCI and alsothe EA calculations.Spectral functions are plotted by approximating the δ -function in the spectral function expression as aLorentzian function δ ( ω ) (cid:39) π (cid:18) ηω + η (cid:19) (A6)for some small constant η . To ensure reproducibility, wespecified the value of η in all relevant figures. ∗ [email protected] † fi[email protected] ‡ [email protected] § [email protected] D. Lu, I. M. Vishik, M. Yi, Y. Chen, R. G. Moore, andZ.-X. Shen, Annu. Rev. Condens. Matter Phys. , 129(2012). R. F. Egerton, Rep. Prog. Phys. , 016502 (2008). C. Andreani, D. Colognesi, J. Mayers, G. Reiter, andR. Senesi, Adv. Phys. , 377 (2005). A. L. Fetter and J. D. Walecka, Quantum Theory of Many-Particle Systems (Dover Books on Physics) (Dover Pub-lications, 2003). G. Onida, L. Reining, and A. Rubio, Rev. Mod. Phys. , 601 (2002). L. Cederbaum and W. Domcke, Adv. Chem. Phys , 205(1977). L. Hedin, J. Michiels, and J. Inglesfield, Phys. Rev. B , 15565 (1998). I. Duchemin and X. Blase, J. Chem. Theory Comput. ,1742 (2020). M. S. Hybertsen and S. G. Louie, Phys. Rev. B , 5390(1986). L. Hedin, Phys. Rev. , A796 (1965). L. Reining, Wiley Interdiscip. Rev. Comput. Mol. Sci. ,e1344 (2018). P. Rinke, A. Qteish, J. Neugebauer, C. Freysoldt, andM. Scheffler, New J. Phys. , 126 (2005). F. Fuchs, J. Furthm¨uller, F. Bechstedt, M. Shishkin, andG. Kresse, Phys. Rev. B , 115109 (2007). N. Marom, F. Caruso, X. Ren, O. T. Hofmann,T. K¨orzd¨orfer, J. R. Chelikowsky, A. Rubio, M. Scheffler,and P. Rinke, Phys. Rev. B , 245127 (2012). F. Bruneval, J. Chem. Phys. , 194107 (2012). D. C. Langreth, Phys. Rev. B , 471 (1970). F. Aryasetiawan, L. Hedin, and K. Karlsson, Phys. Rev.Lett. , 2268 (1996). M. Guzzo, G. Lani, F. Sottile, P. Romaniello, M. Gatti,J. J. Kas, J. J. Rehr, M. G. Silly, F. Sirotti, and L. Rein-ing, Phys. Rev. Lett. , 166401 (2011). J. Lischner, D. Vigil-Fowler, and S. G. Louie, Phys. Rev.Lett. , 146801 (2013). A. Stan, N. E. Dahlen, and R. Van Leeuwen, J. Chem.Phys. , 114105 (2009). U. von Barth and B. Holm, Phys. Rev. B , 8411 (1996). B. Holm and U. von Barth, Phys. Rev. B , 2108 (1998). S. V. Faleev, M. Van Schilfgaarde, and T. Kotani, Phys.Rev. Lett. , 126406 (2004). M. van Schilfgaarde, T. Kotani, and S. Faleev, Phys. Rev.Lett. , 226402 (2006). M. Shishkin and G. Kresse, Phys. Rev. B , 235102(2007). C. Rostgaard, K. W. Jacobsen, and K. S. Thygesen, Phys.Rev. B , 085103 (2010). P. Koval, D. Foerster, and D. S´anchez-Portal, Phys. Rev.B , 155417 (2014). P. Bobbert and W. Van Haeringen, Phys. Rev. B ,10326 (1994). A. Schindlmayr and R. W. Godby, Phys. Rev. Lett. ,1702 (1998). P. Romaniello, S. Guyot, and L. Reining, J. Chem. Phys. , 154111 (2009). W. Chen and A. Pasquarello, Phys. Rev. B , 041115(2015). E. Maggio and G. Kresse, J. Chem. Theory Comput. ,4765 (2017). B. Holm and F. Aryasetiawan, Phys. Rev. B , 12825(1997). J. J. Kas, J. J. Rehr, and L. Reining, Phys. Rev. B ,085112 (2014). J. Lischner, D. Vigil-Fowler, and S. G. Louie, Phys. Rev.B , 125430 (2014). F. Caruso and F. Giustino, Eur. Phys. J. B (2016). M. Z. Mayers, M. S. Hybertsen, and D. R. Reichman,Phys. Rev. B , 081109 (2016). E. Jeckelmann, Phys. Rev. B , 045114 (2002). U. Schollw¨ock and S. R. White, AIP Conf. Proc. , 155(2006). J. Schirmer, Phys. Rev. A , 2395 (1982). J. Schirmer, L. S. Cederbaum, and O. Walter, Phys. Rev.A , 1237 (1983). A. Tarantelli and L. S. Cederbaum, Phys. Rev. A , 1656(1989). J. Wenzel, M. Wormit, and A. Dreuw, J. Chem. TheoryComput. , 4583 (2014). A. Dreuw and M. Wormit, WIREs Comput. Mol. Sci. ,82 (2015). A. Yu. Sokolov, J. Chem. Phys. , 204113 (2018). H. J. Monkhorst, Int. J. Quantum Chem. , 421 (1977). M. Nooijen and J. G. Snijders, Int. J. Quantum Chem. , 55 (1992). M. Nooijen and J. G. Snijders, Int. J. Quantum Chem. , 15 (1993). B. Peng and K. Kowalski, Phys. Rev. A , 062512 (2016). J. McClain, J. Lischner, T. Watson, D. A. Matthews,E. Ronca, S. G. Louie, T. C. Berkelbach, and G. K.-L.Chan, Phys. Rev. B , 235139 (2016). J. McClain, Q. Sun, G. K.-L. Chan, and T. C. Berkel-bach, J. Chem. Theory Comput. , 1209 (2017). Y. Furukawa, T. Kosugi, H. Nishi, and Y.-i. Matsushita,J. Chem. Phys. , 204109 (2018). B. Peng and K. Kowalski, Mol. Phys. , 561 (2018). A. Shee and D. Zgid, J. Chem. Theory Comput. , 6010(2019). R. Blankenbecler and R. L. Sugar, Phys. Rev. D , 1304(1983). F. Becca and S. Sorella, Quantum Monte Carlo Ap-proaches for Correlated Systems (Cambridge UniversityPress, Cambridge, England, UK, 2017). W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Ra-jagopal, Rev. Mod. Phys. , 33 (2001). R. N. Silver, D. S. Sivia, and J. E. Gubernatis, Phys.Rev. B , 2380 (1990). J. E. Gubernatis, M. Jarrell, R. N. Silver, and D. S. Sivia,Phys. Rev. B , 6011 (1991). M. Jarrell and J. E. Gubernatis, Phys. Rep. , 133(1996). M. Motta, D. E. Galli, S. Moroni, and E. Vitali, J. Chem.Phys. , 024107 (2014). M. Motta, D. E. Galli, S. Moroni, and E. Vitali, J. Chem.Phys. , 164108 (2015). J. Otsuki, M. Ohzeki, H. Shinaoka, and K. Yoshimi, Phys.Rev. E , 061302 (2017). G. Bertaina, D. E. Galli, and E. Vitali, Adv. Phys-X ,302 (2017). D. R. Reichman and E. Rabani, J. Chem. Phys. ,054502 (2009). O. Goulko, A. S. Mishchenko, L. Pollet, N. Prokof’ev,and B. Svistunov, Phys. Rev. B , 014102 (2017). T. Dornheim, S. Groth, J. Vorberger, and M. Bonitz, Phys. Rev. Lett. , 255001 (2018). D. M. Ceperley and B. Bernu, J. Chem. Phys. , 6316(1988). N. S. Blunt, A. Alavi, and G. H. Booth, Phys. Rev. Lett. , 050603 (2015). N. S. Blunt, A. Alavi, and G. H. Booth, Phys. Rev. B , 085118 (2018). G. H. Booth, A. J. W. Thom, and A. Alavi, J. Chem.Phys. , 054106 (2009). O. W. Day, D. W. Smith, and R. C. Morrison, J. Chem.Phys. , 115 (1975). D. W. Smith and O. W. Day, J. Chem. Phys. , 113(1975). B. T. Pickup, Chem. Phys. Lett. , 422 (1975). R. C. Morrison, O. W. Day, and D. W. Smith, Int. J.Quantum Chem. , 229 (1975). J. C. Ellenbogen, O. W. Day, D. W. Smith, and R. C.Morrison, J. Chem. Phys. , 4795 (1977). R. C. Morrison and G. Liu, J. Comp. Chem. , 1004(1992). R. C. Morrison, J. Chem. Phys. , 3718 (1992). D. Sundholm and J. Olsen, J. Chem. Phys. , 3999(1993). J. Cioslowski, P. Piskorz, and G. Liu, J. Chem. Phys. , 6804 (1997). P. Kent, R. Q. Hood, M. Towler, R. Needs, and G. Ra-jagopal, Phys. Rev. B - Condensed Matter and MaterialsPhysics , 15293 (1998). J. Olsen and D. Sundholm, Chem. Phys. Lett. , 282(1998). K. Pernal and J. Cioslowski, J. Chem. Phys. , 4359(2001). J. D. Farnum and D. A. Mazziotti, Chem. Phys. Lett. , 90 (2004). K. Pernal and J. Cioslowski, Chem. Phys. Lett. , 71(2005). M. Ernzerhof, J. Chem. Theory Comput. , 793 (2009). D. Vanfleteren, D. Van Neck, P. W. Ayers, R. C. Mor-rison, and P. Bultinck, J. Chem. Phys. (2009),10.1063/1.3130044. M. Piris, J. M. Matxain, X. Lopez, and J. M. Ugalde, J.Chem. Phys. , 174116 (2012). M. Piris, J. M. Matxain, X. Lopez, and J. M. Ugalde,in (Springer, Berlin, Ger-many, 2013) pp. 5–15. U. Bozkaya, J. Chem. Phys. (2013),10.1063/1.4825041. A. R. Welden, J. J. Phillips, and D. Zgid, , 1 (2015),arXiv:1505.05575. U. Bozkaya and A. ¨Unal, J. Phys. Chem. A , 4375(2018). Y. Pavlyukh, Phys. Rev. A , 1 (2018). Y. Pavlyukh, Phys. Status Solidi B Basic Res. , 1(2019). H. Zheng, “First principles quantum Monte Carlo study ofcorrelated electronic systems,” (2016), [Online; accessed12. Jan. 2021]. S. Zhang, J. Carlson, and J. E. Gubernatis, Phys. Rev.Lett. , 3652 (1995). S. Zhang, J. Carlson, and J. E. Gubernatis, Phys. Rev.B , 7464 (1997). S. Zhang and H. Krakauer, Phys. Rev. Lett. , 136401 (2003). J. Carlson, J. Gubernatis, G. Ortiz, and S. Zhang, Phys.Rev. B , 12788 (1999). J. P. LeBlanc, A. E. Antipov, F. Becca, I. W. Bulik,G. K.-L. Chan, C.-M. Chung, Y. Deng, M. Ferrero, T. M.Henderson, C. A. Jim´enez-Hoyos, et al. , Phys. Rev. X ,041041 (2015). B.-X. Zheng, C.-M. Chung, P. Corboz, G. Ehlers, M.-P.Qin, R. M. Noack, H. Shi, S. R. White, S. Zhang, andG. K.-L. Chan, Science , 1155 (2017). M. Motta, D. M. Ceperley, G. K.-L. Chan, J. A. Gomez,E. Gull, S. Guo, C. A. Jim´enez-Hoyos, T. N. Lan, J. Li,F. Ma, et al. , Phys. Rev. X , 031059 (2017). S. Zhang, F. D. Malone, and M. A. Morales, J. Chem.Phys. , 164102 (2018). M. Motta, C. Genovese, F. Ma, Z.-H. Cui, R. Sawaya,G. K.-L. Chan, N. Chepiga, P. Helms, C. Jim´enez-Hoyos,A. J. Millis, et al. , Phys. Rev. X , 031058 (2020). J. Lee, F. D. Malone, and M. A. Morales, J. Chem. Phys. , 064122 (2019). J. Lee, F. D. Malone, and D. R. Reichman, J. Chem.Phys. , 126101 (2020). K. T. Williams, Y. Yao, J. Li, L. Chen, H. Shi, M. Motta,C. Niu, U. Ray, S. Guo, R. J. Anderson, et al. , Phys. Rev.X , 011041 (2020). F. D. Malone, S. Zhang, and M. A. Morales, J. Chem.Theory Comput. , 4286 (2020). J. Lee, F. D. Malone, and M. A. Morales, J. Chem. The-ory Comput. , 3019 (2020). M. Qin, C.-M. Chung, H. Shi, E. Vitali, C. Hubig,U. Schollw¨ock, S. R. White, S. Zhang, et al. , Phys. Rev.X , 031016 (2020). F. D. Malone, A. Benali, M. A. Morales, M. Caffarel,P. R. C. Kent, and L. Shulenburger, Phys. Rev. B ,161104 (2020). F. D. Malone, S. Zhang, and M. A. Morales, J. Chem.Theory. Comput. , 256 (2019). J. Lee and D. R. Reichman, J. Chem. Phys. , 044131(2020). W. Purwanto and S. Zhang, Phys. Rev. E , 056702(2004). M. Motta and S. Zhang, J. Chem. Theory Comput. ,5367 (2017). L. Hedin, Phys. Scr. , 477 (1980). See https://github.com/msaitow/SecondQuantizationAlgebrafor details on how to obtain the source code. E. Neuscamman, T. Yanai, and G. K.-L. Chan, J. Chem.Phys. , 124102 (2009). M. Saitow, Y. Kurashige, and T. Yanai, J. Chem. Phys. , 044118 (2013). D. Zgid, D. Ghosh, E. Neuscamman, and G. K.-L. Chan,J. Chem. Phys. , 194107 (2009). M. Motta and S. Zhang, WIREs Comput. Mol. Sci. ,e1364 (2018). J. Hubbard, Phys. Rev. Lett. , 77 (1959). J. E. Hirsch, Phys. Rev. B , 4059 (1983). Mario Motta, private communication. S. Chen, M. Motta, F. Ma, and S. Zhang, arXiv preprintarXiv:2011.08335 (2020). E. G. Hohenstein, R. M. Parrish, and T. J. Mart´ınez, J.Chem. Phys. , 044103 (2012). R. M. Parrish, E. G. Hohenstein, T. J. Mart´ınez, andC. D. Sherrill, J. Chem. Phys. , 224106 (2012). E. G. Hohenstein, R. M. Parrish, C. D. Sherrill, and T. J. Mart´ınez, J. Chem. Phys. , 221101 (2012). J. Lee, L. Lin, and M. Head-Gordon, J. Chem. TheoryComput. , 243 (2019). T. Schoof, S. Groth, J. Vorberger, and M. Bonitz, Phys.Rev. Lett. , 130402 (2015). Y. Yang, V. Gorelov, C. Pierleoni, D. M. Ceperley, andM. Holzmann, Phys. Rev. B , 085115 (2020). Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth,S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova,S. Sharma, S. Wouters, and G. K.-L. Chan, WIREs Com-put. Mol. Sci. , e1340 (2018). P. R. C. Kent, A. Annaberdiyev, A. Benali, M. C. Ben-nett, E. J. Landinez Borda, P. Doak, H. Hao, K. D. Jor-dan, J. T. Krogel, I. Kyl¨anp¨a¨a, J. Lee, Y. Luo, F. D. Mal-one, C. A. Melton, L. Mitas, M. A. Morales, E. Neuscam-man, F. A. Reboredo, B. Rubenstein, K. Saritas, S. Upad-hyay, G. Wang, S. Zhang, and L. Zhao, J. Chem. Phys. , 174105 (2020). See https://github.com/pauxy-qmc/pauxy for details onhow to obtain the source code. A. A. Holmes, N. M. Tubman, and C. J. Umrigar, J.Chem. Theory Comput. , 3674 (2016). S. Sharma, A. A. Holmes, G. Jeanmairet, A. Alavi, andC. J. Umrigar, J. Chem. Theory Comput. , 1595 (2017). J. E. T. Smith, B. Mussard, A. A. Holmes, andS. Sharma, J. Chem. Theory Comput. , 5468 (2017). L. K. Wagner, M. Bajdich, and L. Mitas, J. Comput.Phys. , 3390 (2009). R. Haydock, V. Heine, and M. J. Kelly, J. Phys. C: SolidState Phys. , 2591 (1975). E. Dagotto, Rev. Mod. Phys. , 763 (1994). M. J. van Setten, F. Caruso, S. Sharifzadeh, X. Ren,M. Scheffler, F. Liu, J. Lischner, L. Lin, J. R. Deslippe,S. G. Louie, C. Yang, F. Weigend, J. B. Neaton, F. Evers,and P. Rinke, J. Chem. Theory Comput. , 5665 (2015). T. H. Dunning, J. Chem. Phys. , 1007 (1989). J. Lee, M. A. Morales, and F. D. Malone, arXiv (2020),2012.12228. M. Motta, S. Zhang, and G. K.-L. Chan, Phys. Rev. B , 045127 (2019). S. Goedecker, M. Teter, and J. Hutter, Phys. Rev. B ,1703 (1996). J. VandeVondele and J. Hutter, J. Chem. Phys. ,114105 (2007). E. Vitali, H. Shi, M. Qin, and S. Zhang, Phys. Rev. B , 085140 (2016).148