Coupled-cluster impurity solvers for dynamical mean-field theory
Tianyu Zhu, Carlos A. Jimenez-Hoyos, James McClain, Timothy C. Berkelbach, Garnet Kin-Lic Chan
CCoupled-cluster impurity solvers for dynamical mean-field theory
Tianyu Zhu, Carlos A. Jiménez-Hoyos, James McClain, Timothy C. Berkelbach,
3, 4, ∗ and Garnet Kin-Lic Chan † Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena CA 91125 Department of Chemistry, Wesleyan University, Middletown CT 06457 Department of Chemistry, Columbia University, New York NY 10027 Center for Computational Quantum Physics, Flatiron Institute, New York NY 10010
We describe the use of coupled-cluster theory as an impurity solver in dynamical mean-field theory (DMFT)and its cluster extensions. We present numerical results at the level of coupled-cluster theory with single anddouble excitations (CCSD) for the density of states and self-energies of cluster impurity problems in the one- andtwo-dimensional Hubbard models. Comparison to exact diagonalization shows that CCSD produces accuratedensity of states and self-energies at a variety of values of 𝑈 ∕ 𝑡 and filling fractions. However, the low costallows for the use of many bath sites, which we define by a discretization of the hybridization directly on the realfrequency axis. We observe convergence of dynamical quantities using approximately 30 bath sites per impuritysite, with our largest 4-site cluster DMFT calculation using 120 bath sites. We suggest coupled cluster impuritysolvers will be attractive in ab initio formulations of dynamical mean-field theory. I. INTRODUCTION
Dynamical mean-field theory (DMFT) and its clus-ter extensions (such as cluster dynamical mean-field the-ory (CDMFT) and the dynamical cluster approximation(DCA) ) approximate the single-particle Green’s functionof an interacting quantum lattice Hamiltonian using the self-energy of a self-consistent impurity model. Computing theimpurity self-energy and Green’s function is thus the mainnumerical task, and falls to the so-called quantum impuritysolver, the focus of this work.In DMFT impurity models, the impurity sites retain the fullinteraction, while the rest of the lattice is replaced by a self-consistent hybridization 𝚫 ( 𝜔 ) . Impurity solvers can be dividedinto two classes based on how they treat this hybridization. Inthe first class (diagrammatic), which includes methods suchas diagrammatic Monte Carlo (diagMC) and some forms ofcontinuous-time quantum Monte Carlo (CT-QMC), 𝚫 ( 𝜔 ) is directly included in the evaluation of the diagrams. In thesecond class (bath-based), which includes methods such as ex-act diagonalization (ED), configuration interaction (CI), numerical and density matrix renormalization group methods(NRG, DMRG), and the variational Gutzwiller ansatz, the hybridization is first unfolded into a discrete bath. The im-purity and bath sites then together define an impurity Hamilto-nian, from which the Green’s function and self-energy can becomputed from a finite system simulation. At zero tempera-ture, this usually involves computing the impurity ground-statewavefunction, and the impurity Green’s function as a correla-tion function. In this work, we will explore an impurity solverbased on coupled-cluster (CC) theory, which falls into thissecond class of solvers.The bottleneck in all bath-based methods is the number ofsites in the impurity Hamiltonian. Even for a small numberof impurity sites, one requires several bath sites per impuritysite to adequately discretize the hybridization. The exponen-tial complexity of exact diagonalization limits calculations toabout 16 sites in total, and thus to only small impurity clus-ters, typically with no more than 4 impurity sites and 2 bathsites per impurity site. There are two common strategies to ameliorate this bottleneck. The first is to restrict the Hilbertspace in which the diagonalization is performed. This was ex-plored by Zgid and Chan with truncated CI, which defines asystematic selection of a reduced set of Slater determinants inwhich to solve the quantum impurity problem. Bravyi pro-vided a formal justification for this truncation, showing thatfor fixed impurity size, the Green’s function can be convergedby a linear combination of Gaussian states with a cost that ispolynomial in the accuracy. CI solvers significantly extend thesize of systems that can be treated by diagonalization meth-ods particularly with respect to the number of bath sites. A second strategy is to parametrize the impurity wavefunctionthrough a non-linear ansatz. NRG, DMRG, and the variationalGutzwiller approximation adopt this latter approach.The CC solvers we explore here also correspond to a non-linear ansatz for the impurity wavefunction. CC theory hassimilar strengths to configuration interaction, for example itcan be used at zero-temperature as well as finite-temperature and with arbitrary interactions and couplings, but it addressesseveral formal deficiencies of the truncated CI wavefunction.Most importantly, the CC parametrization is size-extensive. This means that (much like tensor networks) one can repre-sent product states of disjoint interacting clusters, with a num-ber of parameters linear in the number of clusters, instead ofthe exponential number of parameters in exact diagonaliza-tion. When interactions are not too strong, this makes theCC parametrization exponentially more compact than the CIparametrization for a given accuracy. In ab initio quantumchemistry calculations, CC theory is generally used in placeof CI, except when there are simultaneously strong interac-tions and a large number of degenerate sites.In this work, we study the performance of truncated CC the-ory as an impurity solver in cluster dynamical mean-field the-ory calculations. We will use the 1D and 2D Hubbard models,canonical models of correlated materials, as our test systems.We focus on the lowest order truncation in the CC theory ofthe Green’s function, namely, equation-of-motion CC theorywith single and double excitations (EOM-CCSD), a level oftheory that is widely implemented and available in quantumchemistry packages. To understand the quality of our low-order CC truncation, we will compare to dynamical quanti- a r X i v : . [ c ond - m a t . s t r- e l ] A ug ties obtained from an exact diagonalization solver with smallbath discretizations. We will further exploit the low cost ofthe CC solver to carry out calculations with a large number ofbath sites. Overall, we seek to shed light on the regimes inwhich low-order CC methods are a promising class of impu-rity solvers for dynamical mean-field theories. After this workwas submitted, Shee and Zgid have presented related work that also explores the use of CC impurity solvers in Green’sfunction embedding methods. II. THEORY RECAPITULATIONA. Cluster Dynamical Mean-Field Theory
Cluster dynamical mean-field theory (CDMFT) has beenextensively reviewed and we present only a minimal de-scription sufficient for the numerical considerations in ourwork. We first consider a general translationally invariant lat-tice Hamiltonian with one- and two-particle interactions ̂𝐻 = ∑ 𝑝𝑞 ℎ 𝑝𝑞 𝑎 † 𝑝 𝑎 𝑞 + 12 ∑ 𝑝𝑞𝑟𝑠 𝑉 𝑝𝑞𝑟𝑠 𝑎 † 𝑝 𝑎 † 𝑞 𝑎 𝑟 𝑎 𝑠 (1)where 𝑝, 𝑞, 𝑟, 𝑠 label lattice sites (including spin) and 𝑎 (†) arefermion annihilation (creation) operators. We then consider acomputational unit cell with 𝑁 cluster sites, that tiles thelattice through a set of translation vectors T . Taking k as acorresponding reciprocal space vector taken from an evenlyspaced mesh of 𝑁 𝑘 points in the first Brillouin zone of , thereciprocal space Hamiltonian becomes ̂𝐻 = ∑ 𝑝𝑞 ∈ ∑ k ̃ℎ 𝑝𝑞 ( k ) ̃𝑎 † 𝑝 k ̃𝑎 𝑞 k + 12 ∑ 𝑝𝑞𝑟𝑠 ∈ ∑ k 𝑝 k 𝑞 k 𝑟 ̃𝑉 𝑝𝑞𝑟𝑠 ( k 𝑝 , k 𝑞 , k 𝑟 ) ̃𝑎 † 𝑝 k 𝑝 ̃𝑎 † 𝑞 k 𝑞 ̃𝑎 𝑟 k 𝑟 ̃𝑎 𝑠 k 𝑝 + k 𝑞 − k 𝑟 (2)with ̃𝑎 (†) 𝑝 k = √ 𝑁 𝑘 ∑ T 𝑎 (†) 𝑝 𝑒 𝑖 k ⋅ T where ̃ℎ and ̃𝑉 are the ma-trix elements of ℎ and 𝑉 in the symmetry-adapted basis. Thenon-interacting and interacting Green’s functions, 𝐠 ( k , 𝜔 ) and 𝐆 ( k , 𝜔 ) , are block diagonal in reciprocal space, and are relatedby the block-diagonal lattice self-energy 𝚺 ( k , 𝜔 ) via the Dysonequation 𝐆 ( k , 𝜔 ) = [ ( 𝜔 + 𝜇 ) − ̃𝐡 ( k ) − 𝚺 ( k , 𝜔 ) ] −1 . (3)The cellular Green’s function is related to the reciprocal spaceGreen’s function by 𝐆 ( 𝜔 ) = 𝑁 −1 𝑘 ∑ k 𝐆 ( k , 𝜔 ) . (4)The key quantity to approximate is the lattice self-energythat contains the effects of interactions. In CDMFT, the lat-tice self-energy is taken to be equal to the self-energy of animpurity with 𝑁 sites, i.e. 𝚺 ( k , 𝜔 ) = 𝚺 imp ( 𝜔 ) (5) The impurity model is characterized by a cellular hybridiza-tion 𝚫 ( 𝜔 ) that describes the delocalization effects betweenthe cell and the lattice. Defining the cellular non-interactingHamiltonian ̂ℎ as ̂ℎ = ∑ 𝑝𝑞 ∈ ℎ 𝑝𝑞 𝑎 † 𝑝 𝑎 𝑞 , (6)the hybridization follows as 𝚫 ( 𝜔 ) = ( 𝜔 + 𝜇 ) − 𝐡 − 𝚺 imp ( 𝜔 ) − 𝐆 −1 ( 𝜔 ) (7)The impurity Green’s function 𝐆 imp ( 𝜔 ) is formally definedfrom the zero-temperature generating functional 𝑊 [ 𝐉 ] 𝑊 = ∫ ∫ 𝐜 ̄ 𝐜 𝑒 𝑖𝑆 ( 𝐉 ) (8) 𝑆 = ∫ ∫ 𝑑𝑡𝑑𝑡 ′ [ ̄ 𝐜 𝑇 ( 𝑡 )[( 𝑖𝜕 𝑡 − 𝐡 ) 𝛿 ( 𝑡 − 𝑡 ′ ) − 𝚫 ( 𝑡, 𝑡 ′ )+ 𝐉 ( 𝑡, 𝑡 ′ )] 𝐜 ( 𝑡 ′ ) ] + 𝑉 [ 𝐜 , ̄ 𝐜 ] (9)where 𝐜 , ̄ 𝐜 are vectors of 𝑁 Grassmann variables, 𝑉 [ 𝐜 , ̄ 𝐜 ] is the interaction contribution to the action 𝑆 , 𝐆 imp ( 𝑡, 𝑡 ′ ) = 𝛿𝑊 ∕ 𝛿 𝐉 ( 𝑡, 𝑡 ′ ) , and 𝐆 imp ( 𝜔 ) = 𝜋 ∫ 𝑑𝑡 𝐆 imp (0 , 𝑡 ) 𝑒 𝑖𝜔𝑡 .From the impurity Green’s function, 𝚺 imp ( 𝜔 ) follows as 𝚺 imp ( 𝜔 ) = [ ( 𝜔 + 𝜇 ) − 𝐡 − 𝚫 ( 𝜔 ) ] − 𝐆 −1imp ( 𝜔 ) . (10)Using 𝚺 imp ( 𝜔 ) as the lattice self-energy in (5) leads to a newlattice Green’s function and hybridization, and thus a new im-purity Green’s function. The self-consistency in CDMFT isthen achieved when the cellular Green’s function and impurityGreen’s function agree, 𝐆 imp ( 𝜔 ) = 𝐆 ( 𝜔 ) (11)Note that after self-consistency, if the primitive cell of the lat-tice is smaller than the computational cell, the cellular Green’sfunction may break translational invariance. There are vari-ety of related formulations, that restore the translationalinvariance (including the dynamical cluster approximation )but we will not use them here.In the bath-based CDMFT, the hybridization function is rep-resented by a set of fictitious bath sites and couplings. For-mally, we consider each element of 𝚫 ( 𝜔 ) as the Hilbert trans-form 𝚫 ( 𝜔 ) = ∫ 𝑑𝜀 𝐉 ( 𝜀 ) 𝜔 − 𝜀 (12)with the spectral density 𝐉 ( 𝜔 ) = − 1 𝜋 Im 𝚫 ( 𝜔 + 𝑖𝜂 ) . (13)A bath parametrization can be considered a discrete repre-sentation of the above integral, where 𝜀 and 𝐉 ( 𝜀 ) take on adiscrete set of values given by the bath energies and cou-plings. One common choice is to approximate 𝚫 ( 𝜔 ) alongthe imaginary axis where it is smooth and can be more eas-ily fit to a bath representation by numerical optimization ofa cost function. This is beneficial for exact diagonalization(and related) solvers which can only handle a very small num-ber of bath sites but can lead to a loss of accuracy whenanalytically continuing to the real axis for spectral computa-tion.
Another choice is to approximate 𝚫 ( 𝜔 ) along thereal axis directly, which is commonly done with NRG andDMRG solvers. Here we view the Hilbert transform asa quadrature along the real axis 𝚫 ( 𝜔 ) = 𝑁 𝜔 ∑ 𝑛 =1 𝑤 𝑛 𝐉 ( 𝜀 𝑛 ) 𝜔 − 𝜀 𝑛 (14)where 𝑤 𝑛 are quadrature weights for 𝑁 𝜔 integration gridpoints. Then, we can define an approximate hybridization ofthe form ̃ Δ 𝑝𝑞 ( 𝜔 ) = 𝑁 𝜔 ∑ 𝑛 =1 𝑁 ∑ 𝑘 =1 𝑉 ( 𝑛 ) 𝑝,𝑘 𝑉 ( 𝑛 ) 𝑞,𝑘 𝜔 − 𝜀 𝑛 (15)where 𝐕 ( 𝑛 ) = [ 𝑤 𝑛 𝐉 ( 𝜀 𝑛 )] . This leads to the discrete impurityHamiltonian, ̂𝐻 imp = ̂𝐻 + 𝑁 𝜔 ∑ 𝑛 =1 𝑁 ∑ 𝑘 =1 [ 𝜀 𝑛 𝑎 † 𝑛𝑘 𝑎 𝑛𝑘 + ∑ 𝑝 ( 𝑉 ( 𝑛 ) 𝑝,𝑘 𝑎 † 𝑝 𝑎 𝑛𝑘 + H . c . )] (16)where ̂𝐻 is the cellular Hamiltonian with interactions andthe creation and annihilation operators indexed by 𝑛𝑘 in thelast two terms act on the fictitious bath space. From the impu-rity Hamiltonian, the Green’s functions can then be defined ascorrelation functions. For example, the addition and removalparts of the impurity Green’s function 𝐆 imp ( 𝜔 ) are given by 𝐺 + 𝑝𝑞 ( 𝜔 ) = ⟨ Ψ | 𝑎 𝑝 [ 𝜔 + 𝜇 − ( ̂𝐻 imp − 𝐸 ) + 𝑖𝜂 ] −1 𝑎 † 𝑞 | Ψ ⟩ (17a) 𝐺 − 𝑝𝑞 ( 𝜔 ) = ⟨ Ψ | 𝑎 † 𝑞 [ 𝜔 + 𝜇 − ( 𝐸 − ̂𝐻 imp ) + 𝑖𝜂 ] −1 𝑎 𝑝 | Ψ ⟩ (17b)The goal of the impurity solver in the bath-based represen-tation is thus to approximate the impurity model ground-statewavefunction | Ψ ⟩ and the impurity Green’s function 𝐆 imp ( 𝜔 ) via the expressions in (17). We next describe how this is donein CC theory. B. Coupled-cluster theory of the ground-state and Green’sfunction
We now describe the basics of CC theory. Detailed discus-sions of ground-state CC and the CC Green’s function canbe found in Refs. 40–45. The CC ansatz is an exponentialparametrization of the wavefunction, | Ψ ⟩ = 𝑒 ̂𝑇 | Φ ⟩ (18) where | Φ ⟩ is a reference determinant, and ̂𝑇 is the clusterexcitation operator ̂𝑇 = ∑ 𝑖𝑎 𝑡 𝑎𝑖 𝑎 † 𝑎 𝑎 𝑖 + 14 ∑ 𝑖𝑗𝑎𝑏 𝑡 𝑎𝑏𝑖𝑗 𝑎 † 𝑎 𝑎 † 𝑏 𝑎 𝑖 𝑎 𝑗 + …= ̂𝑇 + ̂𝑇 + … (19)where indices 𝑖, 𝑗, … and 𝑎, 𝑏, … , label particle (p) and hole (h)orbitals in the reference determinant, and 𝑡 𝑎𝑖 , 𝑡 𝑎𝑏𝑖𝑗 , … are 1p1h,2p2h, etc. cluster amplitudes. Truncating the CC amplitudesat 1p1h, 2p2h, etc. gives the coupled-cluster singles (CCS),coupled-cluster singles and doubles (CCSD), etc. approxima-tions. CCS constitutes a mean-field ansatz. Thus we use theCCSD ansatz, the lowest truncation that includes correlations,in this work.Given a Hamiltonian ̂𝐻 , the reference determinant | Φ ⟩ isoften chosen to be the mean-field (Hartree-Fock) ground-stateof ̂𝐻 (although including 𝑒 ̂𝑇 in the cluster operator renders theapproximation somewhat insensitive to the choice of determi-nant, as an arbitrary determinant satisfies | Φ ′ ⟩ = 𝑒 ̂𝑇 | Φ ⟩ ).The amplitudes are chosen to satisfy the projected CC equa-tions; for CCSD, these are 𝐸 = ⟨ Φ | 𝑒 − ̂𝑇 ̂𝐻𝑒 ̂𝑇 | Φ ⟩ ⟨ Φ 𝑎𝑖 | 𝑒 − ̂𝑇 ̂𝐻𝑒 ̂𝑇 | Φ ⟩ ⟨ Φ 𝑎𝑏𝑖𝑗 | 𝑒 − ̂𝑇 ̂𝐻𝑒 ̂𝑇 | Φ ⟩ (20)where ⟨ Φ 𝑎𝑖 | = ⟨ Φ | 𝑎 𝑖 𝑎 † 𝑎 … and 𝐸 is the (non-variational) CCapproximation to the ground-state energy. The operator 𝑒 ̂𝑇 de-fines a similarity transformation. Thus ̄𝐻 = 𝑒 − ̂𝑇 ̂𝐻𝑒 ̂𝑇 is an ef-fective Hamiltonian whose mean-field energy is 𝐸 , or equiv-alently, the CC wavefunction is a product state of similaritytransformed quasi-particles, defined by the quasi-particle op-erators ̄𝑎 † 𝑖 = 𝑒 − ̂𝑇 𝑎 † 𝑖 𝑒 ̂𝑇 .The advantages of truncated CC versus truncated CI derivefrom the exponentiation of the cluster operator. A CI wave-function can be written in similar notation as | Ψ ⟩ = ̂𝐶 | Φ ⟩ ̂𝐶 = 1 + ∑ 𝑖𝑎 𝑐 𝑎𝑖 𝑎 † 𝑎 𝑎 𝑖 + 14 ∑ 𝑖𝑗𝑎𝑏 𝑐 𝑎𝑏𝑖𝑗 𝑎 † 𝑎 𝑎 † 𝑏 𝑎 𝑖 𝑎 𝑗 + …= 1 + ̂𝐶 + ̂𝐶 + … (21)The truncated ̂𝐶 can be seen as a linearization of a truncated 𝑒 ̂𝑇 , or equivalently, the terms in ̂𝑇 are the cumulants of ̂𝐶 ,rewriting the latter as products of cluster excitations, e.g. ̂𝐶 = ̂𝑇 + ̂𝑇 . The extensivity of the truncated CC approximationalso derives from the exponential structure, since for a system 𝐴𝐵 consisting of two separated clusters of sites 𝐴 and 𝐵 , wecan write 𝑒 ̂𝑇 𝐴𝐵 | Φ 𝐴𝐵 ⟩ = 𝑒 ̂𝑇 𝐴 | Φ 𝐴 ⟩ 𝑒 ̂𝑇 𝐵 | Φ 𝐵 ⟩ .If the interactions between particle-hole excitations are nottoo strong, we expect higher-order cluster operators to becomesmall. In this situation, the truncated CC approximation pro-vides a significant improvement over the truncated CI ansatz.Low-order CC truncations are accurate if the state of interestcan be described in terms of low-order products of fluctua-tions around the chosen reference mean-field state | Φ ⟩ . Inthis work, we will consider two possible CCSD formulationsbased on different reference mean-field states corresponding torestricted Hartree-Fock (RHF) and unrestricted Hartree-Fock(UHF). In restricted CCSD (RCCSD), the reference mean-field state is required to be a singlet ( 𝑆 = 0 ) eigenfunction of ̂𝑆 , and the cluster operators ̂𝑇 , ̂𝑇 each commute with ̂𝑆 . Inunrestricted CCSD (UCCSD), the reference mean-field stateis only required to be an eigenfunction of ̂𝑆 𝑧 , and the clus-ter operators commute with ̂𝑆 𝑧 . Consequently, spin symmetrycan be broken in UCCSD, for example, to yield an antiferro-magnetic state. RCCSD is expected to be most accurate whenthe interacting state can be described by fluctuations around aparamagnetic state, while UCCSD is more appropriate to de-scribe fluctuations around an (anti)ferromagnetic state.From the CC ground state, we then compute the Green’sfunction as a correlation function. Using the cluster operator,similarity transformed Hamiltonian, and quasiparticle opera-tors, we can rewrite the exact expressions for the Green’s func-tion in (17) in the equation-of-motion coupled-cluster (EOM-CC) form 𝐺 + 𝑝𝑞 ( 𝜔 ) = ∑ 𝑚𝑛 ⟨ Λ | ̄𝑎 𝑝 | Φ 𝑚 ⟩⟨ Φ 𝑛 | ̄𝑎 † 𝑞 | Φ ⟩ × ⟨ Φ 𝑚 | [ 𝜔 + 𝜇 − ( ̄𝐻 imp − 𝐸 ) + 𝑖𝜂 ] −1 | Φ 𝑛 ⟩ (22a) 𝐺 − 𝑝𝑞 ( 𝜔 ) = ∑ 𝑚𝑛 ⟨ Λ | ̄𝑎 † 𝑞 | Φ 𝑚 ⟩⟨ Φ 𝑛 | ̄𝑎 𝑝 | Φ ⟩ × ⟨ Φ 𝑚 | [ 𝜔 + 𝜇 − ( 𝐸 − ̄𝐻 imp ) + 𝑖𝜂 ] −1 | Φ 𝑛 ⟩ (22b)where ⟨ Λ | is the left eigenstate of ̄𝐻 , and ∑ 𝑚 | Φ 𝑚 ⟩⟨ Φ 𝑚 | =1 , where | Φ 𝑚 ⟩ is a determinant. Defining response vectors | 𝑅 ± 𝑝 ( 𝜔 ) ⟩ , | 𝑅 + 𝑞 ( 𝜔 ) ⟩ = ̂𝑃 + [ 𝜔 + 𝜇 − ( ̄𝐻 imp − 𝐸 ) + 𝑖𝜂 ] −1 ̂𝑃 + ̄𝑎 † 𝑞 | Φ ⟩ (23a) | 𝑅 − 𝑝 ( 𝜔 ) ⟩ = ̂𝑃 − [ 𝜔 + 𝜇 − ( 𝐸 − ̄𝐻 imp ) + 𝑖𝜂 ] −1 ̂𝑃 − ̄𝑎 𝑝 | Φ ⟩ , (23b)where ̂𝑃 ± is the projector onto 1p, 2p1h, ... states or 1h, 2h1p,... states, allows the CC Green’s functions to be efficiently com-puted as 𝐺 + 𝑝𝑞 ( 𝜔 ) = ⟨ Φ | ̄𝑎 𝑝 | 𝑅 + 𝑞 ( 𝜔 ) ⟩ (24a) 𝐺 − 𝑝𝑞 ( 𝜔 ) = ⟨ Φ | ̄𝑎 † 𝑞 | 𝑅 − 𝑝 ( 𝜔 ) ⟩ . (24b)A truncated EOM-CC approximation to the Green’s func-tion contains two truncations, one of the ̂𝑇 operator that de-fines the ground-state, and another of the resolution of theidentity determinants in the Lehmann sum. In this work wewill use the EOM-CCSD approximation, where ̂𝑇 correspondsto the ground-state CCSD truncation, and where | Φ 𝑚 ⟩ , | Φ 𝑛 ⟩ are restricted to 1h, 2h1p ( 𝐆 − ) and 1p, 2p1h ( 𝐆 + ) excitations out of the reference determinant. Note that ̄𝑎 𝑝 | Φ ⟩ = ( 𝑎 𝑝 + [ 𝑎 𝑝 , ̂𝑇 ]) | Φ ⟩ (25a) ̄𝑎 † 𝑝 | Φ ⟩ = ( 𝑎 † 𝑝 + [ 𝑎 † 𝑝 , ̂𝑇 ]) | Φ ⟩ (25b)thus the truncation of ̂𝑇 at the CCSD level implies that ̄𝑎 𝑝 | Φ ⟩ and ̄𝑎 † 𝑝 | Φ ⟩ are expressible in terms of 1h, 2h1p and 1p, 2p1hspaces, respectively.From the Green’s function one can compute all one-particleexpectation values (such as the particle number and the single-particle density matrix) as well as the total energy, from theMigdal formula 𝐸 = − 𝜋 Im ∫ 𝜇 −∞ 𝑑𝜔 Tr( 𝜔 + 𝐡 ) 𝐆 ( 𝜔 ) . How-ever, EOM-CCSD is not a conserving approximation, thusthe energy computed using the Migdal formula is differentfrom the ground-state CC energy 𝐸 . Nonetheless, the single-particle density matrix 𝛾 𝑝𝑞 = (2 𝜋𝑖 ) −1 ∫ 𝜇 −∞ ( 𝐺 − 𝑝𝑞 ( 𝜔 ) + 𝐺 + 𝑝𝑞 ( 𝜔 )) is correctly normalized and equal to its definition as an en-ergy derivative of the CC energy functional. Furthermore, theEOM-CCSD Green’s function is not strictly causal. One im-plication of this is that the impurity self-energy calculated viaEq. (10) may not have an imaginary part that is negative def-inite, similar to the behavior observed in adaptively truncatedCI solvers. A futher discussion of this point and partial so-lution is described in the Appendix.The computational cost of ground-state CCSD for a general(e.g. quantum chemistry) two-particle interaction is ( 𝑁 ) where 𝑁 is the number of sites (more specifically, the costis ( 𝑜 𝑣 + 𝑜 𝑣 ) where 𝑜 is the number of electrons and 𝑣 = 𝑁 − 𝑜 is the number of unoccupied states). There are somespecial considerations, however, when using CC to determinethe ground-state of the CDMFT impurity Hamiltonian. Theimpurity Hamiltonian only has two-particle interactions on theimpurity sites when working in the site basis. However, the CCequations are usually implemented in the mean-field molecu-lar orbital basis, for which the impurity model has two-particleinteractions over all orbitals. Although we do not take advan-tage of it here, in principle, using the locality of the interactionin the site basis can significantly lower the computational costof CCSD, particularly if we consider the scaling of the costwhen increasing the number of bath sites while the numberof impurity sites is kept fixed. Another consideration is thatthe CC theory presented here is a zero-temperature pure statetheory, with a fixed particle number. This is in contrast to thegrand canonical formulation of (C)DMFT in Sec. II A. Thismeans that it is necessary to search over all particle numbers(and spin sectors) in the ground-state calculations of the impu-rity to find the quantum numbers of the (lowest energy) groundstate, as done with other zero-temperature impurity solvers. The cost of EOM-CCSD (without accounting for local-ity in the interaction) is ( 𝑁 ) per response vector and thus ( 𝑁 𝜔 𝑁 𝑁 ) for all elements of the frequency-dependent im-purity Green’s function. Computing the response vector canbe done either by solving a system of linear equations or bycomputing | 𝑅 ± 𝑝 ( 𝑡 ) ⟩ in the time domain followed by Fouriertransformation (similarly as done in td-DMRG solvers ).We have found that a generalized minimum residual (GMRES)solver or a simplified two-parameter generalized conjugateresidual method with inner orthogonalization and outer trun-cation (GCROT( 𝑚, 𝑘 )) works well for the linear equationsand converges in (10) iterations when the response vector ata nearby frequency is used to initialize the solution for the re-sponse vector at a new frequency and the matrix diagonal isused as a preconditioner.Once the impurity Green’s function is computed, it maythen be used in all the expressions in Sec. II B. We note thatthe paramagnetic formulation of CDMFT is incompatible withsymmetry breaking that may occur in the UCCSD impuritysolver. To study a paramagnetic phase using the UCCSDsolver, we spin-average the impurity Green’s function, 𝐆 imp = ∑ 𝜎 = ↑ , ↓ 𝐺 𝜎𝜎 imp . C. Algorithm
For completeness, we outline the full computational proce-dure for our CDMFT algorithm at fixed chemical potential us-ing a CC solver. The loop is initialized with the Hartree-Fockhybridization.1. Discretize hybridization via (29).2. Solve impurity problem with Hartree-Fock at fixed 𝜇 (electron number can change).3. Calculate impurity Green’s function with CC via (23)and (24) using GMRES or GCROT( 𝑚, 𝑘 ).4. Calculate hybridization 𝚫 ( 𝜔 ) via (7) and check for con-vergence. If not converged, update hybridization (usingDIIS ) and return to 1 .If results are desired at a fixed occupancy, then the occupancycan be calculated via 𝑛 = Tr γ imp ∕ 𝑁 using the ground-stateCC solution. If the occupancy does not equal the target occu-pancy, then the chemical potential is updated and the CDMFTloop is repeated.We have implemented the above algorithm using the HF,CCSD, and EOM-CCSD routines from the PySCF quantumchemistry package. III. APPLICATIONS TO THE 1D AND 2D HUBBARDMODELS
The Hubbard model is defined by the lattice Hamiltonian ̂𝐻 = − 𝑡 ∑ ⟨ 𝑝𝑞 ⟩ ,𝜎 𝑎 † 𝑝𝜎 𝑎 𝑞𝜎 + 𝑈 ∑ 𝑝 𝑛 𝑝 ↑ 𝑛 𝑝 ↓ (26)where 𝑛 𝑝𝜎 = 𝑎 † 𝑝𝜎 𝑎 𝑝𝜎 . It is the canonical model for DMFTstudies as there are no nonlocal interactions. A. Implementation in the 1D and 2D Hubbard models
We will consider 1 and 2 site clusters for the 1D Hub-bard model and 1 and 4 site (2 ×
2) clusters for the 2D Hub-bard model. The 2 site and 4 site clusters contain additional point group symmetry which allows us to define a symmetry-adapted impurity orbital basis, associated with operators 𝑎 †Γ + = 1 √ 𝑎 †1 + 𝑎 †2 ) (27a) 𝑎 †Γ − = 1 √ 𝑎 †1 − 𝑎 †2 ) (27b)for two sites, and 𝑎 †Γ = 12 ( 𝑎 †1 + 𝑎 †2 + 𝑎 †3 + 𝑎 †4 ) (28a) 𝑎 †Γ = 12 ( 𝑎 †1 + 𝑎 †2 − 𝑎 †3 + 𝑎 †4 ) (28b) 𝑎 †Γ = 12 ( 𝑎 †1 − 𝑎 †2 − 𝑎 †3 + 𝑎 †4 ) (28c) 𝑎 †Γ = 12 ( 𝑎 †1 − 𝑎 †2 + 𝑎 †3 − 𝑎 †4 ) (28d)for four sites. The impurity Green’s function, self-energy,and hybridization are diagonal in the symmetry-adapted or-bital basis. The diagonal hybridization leads to a simpler formof the bath discretization (15), which can now be written as ̃ Δ Γ ( 𝜔 ) = 𝑁 𝜔 ∑ 𝑛 =1 𝑁 ∑ 𝑘 =1 | 𝑉 Γ ,𝑛𝑘 | 𝜔 − 𝜀 𝑛 . (29)The one-electron couplings in the site basis are 𝑉 𝑝,𝑛𝑘 = ∑ Γ 𝑈 𝑝, Γ 𝑉 Γ ,𝑛𝑘 where 𝐔 is the frequency-independent matrix ofsymmetry-adapted eigenvectors associated with the change ofbasis.The bath energies and weights used to discretize the hy-bridization are chosen according to Gauss-Legendre quadra-ture on the interval [−7 𝑡 + 𝑈 ∕2 , +7 𝑡 + 𝑈 ∕2] for the 1D Hub-bard model and [−9 𝑡 + 𝑈 ∕2 , +9 𝑡 + 𝑈 ∕2] for the 2D Hubbardmodel. In almost all calculations, we achieve convergence us-ing a small imaginary broadening 𝜂 ∕ 𝑡 = 0 . . In some strong-coupling cases where convergence was especially challenging,we use a slightly large value 𝜂 ∕ 𝑡 = 0 . . This latter value resultsin a weaker hybridization and also facilitates the solution ofthe CC linear response equations (23), similar to previous ob-servations with DMRG solvers. At convergence, quantitieson the real frequency axis are plotted with a larger broadening 𝜂 ∕ 𝑡 = 0 . , for visual clarity only. B. 1D Hubbard model
We first consider the 1D Hubbard model. We present re-sults at 𝑈 ∕ 𝑡 = 2 and 𝑈 ∕ 𝑡 = 6 (weak and strong couplingrelative to the single-particle bandwidth of 𝑡 ). We showthe impurity density of states (DOS) on the real frequencyaxis 𝜌 ( 𝜔 ) = −( 𝜋𝑁 ) −1 TrIm 𝐆 ( 𝜔 ) and the imaginary partof the impurity self-energy on the imaginary frequency axis ImΣ ( 𝑖𝜔 𝑛 ) . For visual clarity, the DOS is plotted with abroadening of 𝜂 ∕ 𝑡 = 0 . .We first assess the accuracy of CCSD compared to exact di-agonalization. In Fig. 1, we show the DOS and self-energy U / t = U / t = U / t = U / t = DO S ρ ( ω ) Frequency ( ω − U / / t CCSD (1,9)Exact (1,9)0.000.050.100.150.200.250.30 -6 -4 -2 0 2 4 6 Frequency ( ω − U / / t I m Σ ( i ω n ) / t Imag. frequency ω n / t CCSD (2,8)Exact (2,8)-0.25-0.20-0.15-0.10-0.050.00 0 2 4 6 8 10 Imag. frequency ω n / t CCSD (1,9)Exact (1,9)CCSD (2,8)Exact (2,8)-4.0-3.0-2.0-1.00.0 0 2 4 6 8 10
FIG. 1. Single-site DMFT and two-site CDMFT results for the 1DHubbard model, comparing the use of the CC and exact diagonal-ization (ED) impurity solvers. Numbers in parentheses indicate thenumber of impurity sites and the total number of bath sites. Resultsare shown at half-filling with weak interactions ( 𝑈 ∕ 𝑡 = 2 , left) andstrong interactions ( 𝑈 ∕ 𝑡 = 6 , right). The chemical potential is fixedat 𝜇 = 𝑈 ∕2 and the DOS is plotted with a broadening of 𝜂 = 0 . 𝑡 forclarity. from single-site DMFT with 9 bath sites and from two-siteCDMFT with 4 bath sites per impurity site (8 in total). Bothcorrespond to an impurity problem with 10 sites which is read-ily accessible with exact diagonalization. At both values of 𝑈 ∕ 𝑡 , the RCCSD and ED plots are indistinguishable. Note thatat small 𝑈 ∕ 𝑡 , UHF does not break spin symmetry and UCCSDand RCCSD give identical results. At large 𝑈 ∕ 𝑡 , UHF stronglybreaks spin symmetry, however the mean-field AFM order isreduced by UCCSD such that the final results are almost indis-tinguishable from those of RCCSD. Therefore, all results forthe 1D Hubbard model are presented for RCCSD only.We next assess convergence of the DOS with respect to thenumber of bath sites, which requires impurity problem sizesbeyond the reach of ED. In Fig. 2, we show the single-siteDMFT DOS and self-energy computed using 9, 19, and 29bath sites and the RCCSD solver. The largest impurity prob-lem involves 30 particles in 30 orbitals. The plots are qualita-tively converged with 19 bath sites and converged to the eyewith 29 bath sites. Consistent with previous studies, wefind that the single-site DMFT produces a Kondo-like reso-nance in the DOS at large 𝑈 , as compared with exact Mottinsulating behavior at all 𝑈 . We note that this Kondo-like be-havior is a result of using a single impurity site in DMFT andis not due to the CCSD approximation.Based on the converged results for single-site DMFT, wenext use two-site CDMFT with at least 30 bath sites per impu-rity, i.e. 60 bath sites in total. In Fig. 3, we present the DOSas a function of occupancy. At half-filling ( 𝑛 = 1 , 𝜇 = 𝑈 ∕2 ), U / t = U / t = U / t = U / t = DO S ρ ( ω ) Frequency ( ω − U / / t ω − U / / t I m Σ ( i ω n ) / t Imag. frequency ω n / tN ω = N ω = N ω = ω n / tN ω = N ω = N ω = FIG. 2. Single-site DMFT results for the 1D Hubbard at half-fillingwith weak interactions ( 𝑈 ∕ 𝑡 = 2 , left) and strong interactions ( 𝑈 ∕ 𝑡 =6 , right), showing convergence with respect to the number of bathsites 𝑁 𝜔 . The chemical potential is fixed at 𝜇 = 𝑈 ∕2 and the DOS isplotted with a broadening of 𝜂 ∕ 𝑡 = 0 . for clarity. we see a clear Mott gap proportional to 𝑈 ∕ 𝑡 . We access otherfilling fractions by changing the chemical potential. At small 𝑈 ∕ 𝑡 , only minor changes are seen in the DOS for moderatechanges in the occupancy. At larger 𝑈 ∕ 𝑡 , there is significantredistribution of the spectral weight towards lower energy, cre-ating a metallic DOS around the chemical potential, which isindicated by a vertical line in the DOS plots. In the bottompanel of Fig. 3, we show the occupancy as a function of thechemical potential for the two values of 𝑈 studied. The dis-crete nature of the bath allows access to only a discrete set ofoccupations, and so this latter data was obtained using largerbath sizes. The appearance of a Mott plateau and suppressedcompressibility is clearly seen for 𝑈 ∕ 𝑡 = 6 . Due to the dis-cretized bath, there are artificial plateaus in the regions wherethe Bethe ansatz (BA) varies monotonically, even with 40 bathsites per impurity site. By increasing the bath size to 80 and60 bath sites per impurity site for 𝑈 ∕ 𝑡 = 2 and , we showthat such artificial behavior can be removed, and the occu-pancy converges to the BA solution. This result may suggestthat when discretizing the hybridization along real axis, staticquantities (e.g. occupancy) converge more slowly with respectto the number of bath orbitals as compared to dynamic quan-tities (such as the spectral function) although a more carefulcomparison with imaginary axis techniques is required. C. 2D Hubbard model
We next study the 2D Hubbard model at half-filling. Weperform calculations at 𝑈 ∕ 𝑡 = 2 and 𝑈 ∕ 𝑡 = 8 . This is aboveand below the CDMFT paramagnetic Mott transition U / t = U / t = DO S ρ ( ω ) Frequency ( ω − U / / tn = . n = . n = . ω − U / / tn = . n = . n = . S it e o cc up a ti on n Chemical potential ( µ − U / / tU / t = N ω = U / t = N ω = U / t = N ω = U / t = N ω = FIG. 3. Two-site CDMFT results as a function of doping. The DOSobtained using 𝑁 𝜔 = 30 (60 bath sites in total) is shown in the top twopanels, where the chemical potential is indicated by a vertical line.The site occupancy as a function of the chemical potential using 𝑁 𝜔 =40 (open symbols) and 𝑁 𝜔 = 60 , (filled symbols) is shown inthe bottom panel, along with the numerically exact result from Betheansatz and the mean-field result. around 𝑈 ∕ 𝑡 ≈ 6 . In Fig. 4, we show results for single-site DMFT and
CDMFT. As before, we fix the numberof bath sites to be about 30 per impurity site, such that ourlargest calculations have 120 particles in 120 orbitals.At 𝑈 ∕ 𝑡 = 2 DMFT and CDMFT give very similar DOSand self-energies; i.e. the effect of the impurity size is small.Even if we allow for symmetry breaking in the cluster,we find it to be very weak, thus (paramagnetic) RCCSD and(weakly antiferromagnetic) UCCSD give very similar results.In both cases, the gap is smaller than the bath discretizationand thus indistinguishable from a metal. At 𝑈 ∕ 𝑡 = 8 , para-magnetic single-site DMFT with RCCSD or UCCSD givesidentical results, showing a Kondo-like feature centered be-tween the lower and upper Hubbard bands. Paramagnetic 2 ×2 CDMFT with RCCSD fails to converge in the DMFT loop.The first few iterations of the DMFT cycle are physically rea-sonable and show the opening of a paramagnetic Mott gap.However, eventually an impurity problem is constructed forwhich the iterative solution of the ground-state RCCSD equa-tions fails to converge, because the associated CC amplitudesbecome large. If we instead carry out an antiferromagneticCDMFT calculation allowing for symmetry breaking in theUCCSD solver, the DMFT cycle converges smoothly. As ex-pected, the antiferromagnetic DOS exhibits a gap, with a qual-itatively reasonable size of about 𝑡 . U / t = U / t = U / t = U / t = DO S ρ ( ω ) Frequency ( ω − U / / tn imp =
1, RCCSD n imp =
4, RCCSD n imp =
4, UCCSD0.000.050.100.150.200.250.30 -8 -6 -4 -2 0 2 4 6 8 Frequency ( ω − U / / tn imp =
1, RCCSD n imp =
4, UCCSD0.000.050.100.150.200.250.30 -8 -6 -4 -2 0 2 4 6 8 I m Σ ( i ω n ) / t Imag. frequency ω n / t -0.15-0.10-0.050.00 0 5 10 15 20 Imag. frequency ω n / t -3.0-2.5-2.0-1.5-1.0-0.50.0 0 5 10 15 20 FIG. 4. Single-site DMFT and four-site ( ) CDMFT results forthe 2D Hubbard at half-filling with weak interactions ( 𝑈 ∕ 𝑡 = 2 , left)and strong interactions ( 𝑈 ∕ 𝑡 = 8 , right). We used 30 bath sites forsingle-site DMFT and 120 bath sites for the four-site CDMFT. Thechemical potential is fixed at 𝜇 = 𝑈 ∕2 and the DOS is plotted with abroadening of 𝜂 ∕ 𝑡 = 0 . for clarity. IV. CONCLUSIONS
We have demonstrated the promise of coupled-cluster (CC)theory as an impurity solver in single-site dynamical mean-field theory and multi-site cluster dynamical mean-field the-ory. In particular, the polynomial cost of truncated CC allowsthe use of many bath sites per impurity site, which enables afaithful discretization of the hybridization directly on the realfrequency axis. In our studies of the Hubbard model, we findthat the Green’s functions and self-energies are well convergedusing approximately 30 bath sites per impurity site.Despite the strongly correlated nature of the lattice problem,the low-order CCSD truncation provides very accurate resultsfor spectral functions, self-energies, and occupation numbersas an impurity solver in DMFT, when there is a modest num-ber of impurities. This is consistent with the impurity problembeing less strongly correlated than the lattice problem, even fora dense discretization of the hybridization. Consequently, weexpect that CC impurity solvers will find most use with clus-ters of a moderate size in real-space. The low cost of the CCimpurity solvers also makes them very promising for applica-tions to real materials with many electrons and many orbitalsper impurity site in ab initio quantum chemical formulationsof DMFT. In particular, the embedding approach to CC re-sponse functions in solids is a promising alternative to full pe-riodic CC calculations, and work along these lines is cur-rently in progress.
ACKNOWLEDGMENTS
TZ and GKC were supported by the US Department of En-ergy, Office of Science, via grant number SC19390. GKCis also supported by the Simons Foundation, via the Many-Electron Collaboration, and via the Simons Investigator Pro-gram. The Flatiron Institute is a division of the Simons Foun-dation.
Appendix: Causality in coupled-cluster theory
As discussed in Sec. II B, EOM-CCSD is not a conservingapproximation and thus the resulting Green’s functions andself-energies need not be causal. The lack of strict causality isnot seen in any of the results we showed above. However, it ispossible to observe small violations of causality if we carry outcalculations with a very small broadening (e.g. 𝜂 = 0 . 𝑡 ).An example of non-causal behaviour that can occur in this set-ting is shown in Fig. A.1. As defined in Eq. (10), 𝚺 imp ( 𝜔 ) should always have a negative imaginary part if it is computedfrom a causal impurity Green’s function. However, we seethat the CCSD impurity self-energy develops a positive imag-inary part exactly around the bath frequencies. One way tounderstand this is that truncated CC theory does not includeall possible diagrams involving interactions and hopping to thebath (i.e. the hybridization). Thus the hybridization contribu-tion to 𝐆 −1 imp in Eq. (10) and 𝚫 ( 𝜔 ) do not precisely cancel togive a causal self-energy. -0.010-0.0050.0000.0050.010 0 1 2 3 4 5 U / t = , η/ t = . I m Σ ( ω ) Frequency ( ω − U / / t Σ imp ( ω ) Σ imp + bath ( ω ) ω bath -0.010-0.0050.0000.0050.010 0 1 2 3 4 5 FIG. A.1. Impurity self-energy on the real frequency axis from atwo-site CDMFT calculation for the 1D Hubbard at half-filling with 𝑈 ∕ 𝑡 = 2 . RCCSD is used as the impurity solver and two bath sitesare coupled to each impurity site. A small broadening of 𝜂 = 0 . 𝑡 is used to show the non-causal behavior. This non-causality only appears when the broadening isvery small and does not affect the results we presented. Nev-ertheless, we observe that there is a simple procedure that re-moves this issue in practice. Instead of computing 𝚺 imp ( 𝜔 ) inEq. (10), one can first compute the self-energy for the wholeimpurity plus bath system 𝚺 imp+bath ( 𝜔 ) = 𝐆 −10 , imp+bath ( 𝜔 ) − 𝐆 −1imp+bath ( 𝜔 ) . (A.1)where 𝐆 −10 , imp+bath ( 𝜔 ) is the non-interacting Green’s functionof the impurity plus bath system. The impurity self-energyis then defined as the impurity block of 𝚺 imp+bath ( 𝜔 ) . (In anexact solver, this is the only non-zero part of 𝚺 imp+bath ( 𝜔 ) ).As shown in Fig. A.1, the non-causal behavior disappears ifthe impurity self-energy is defined in this way, at the cost ofmore computation. While we do not claim that this proceduremakes the self-energy strictly causal under all circumstances,we find that it removes even the small degree of non-causalbehaviour observed at very small broadenings in this work. ∗ [email protected] † [email protected] A. Georges and G. Kotliar, Phys. Rev. B , 6479 (1992). A. Georges, G. Kotliar, W. Krauth, and M. Rozenberg, Rev. Mod.Phys. , 13 (1996). A. I. Lichtenstein and M. I. Katsnelson, Phys. Rev. B - Con-dens. Matter Mater. Phys. , 9283 (2000), arXiv:9911320v1[arXiv:cond-mat]. G. Kotliar, S. Y. Savrasov, G. Pálsson, and G. Biroli, Phys. Rev. Lett. , 186401 (2001), arXiv:0010328 [cond-mat]. M. Hettler, A. Tahvildar-Zadeh, M. Jarrell, T. Pruschke, andH. Krishnamurthy, Phys. Rev. B , 7475 (1998). M. Hettler, M. Mukherjee, M. Jarrell, and H. Krishnamurthy,Phys. Rev. B , 12739 (2000). E. Kozik, K. Van Houcke, E. Gull, L. Pollet, N. Prokof’ev, B. Svis-tunov, and M. Troyer, EPL (Europhysics Lett. , 10004 (2010). P. Werner, A. Comanac, L. de’ Medici, M. Troyer, and A. J. Millis,Phys. Rev. Lett. , 076405 (2006). E. Gull, P. Werner, O. Parcollet, and M. Troyer, EPL (EurophysicsLett. , 57003 (2008). E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer,and P. Werner, Rev. Mod. Phys. , 349 (2011). M. Caffarel and W. Krauth, Phys. Rev. Lett. , 1545 (1994). M. Capone, L. de’ Medici, and A. Georges, Phys. Rev. B ,245116 (2007). D. Zgid and G. K.-L. Chan, J. Chem. Phys. , 094115 (2011). K. G. Wilson, Rev. Mod. Phys. , 773 (1975). S. R. White, Phys. Rev. Lett. , 2863 (1992). R. Bulla, T. A. Costi, and T. Pruschke, Rev. Mod. Phys. , 395(2008). M. C. Gutzwiller, Phys. Rev. Lett. , 159 (1963). J. Čížek, J. Chem. Phys. , 4256 (1966). S. Bravyi and D. Gosset, Commun. Math. Phys. , 451 (2017). D. Zgid, E. Gull, and G. K. L. Chan, Phys. Rev. B - Condens.Matter Mater. Phys. , 1 (2012), arXiv:1203.1914. C. Lin and A. A. Demkov, Phys. Rev. B , 035123 (2013). Y. Lu, M. Höppner, O. Gunnarsson, and M. Haverkort, Phys. Rev.B , 085102 (2014). A. Go and A. J. Millis, Phys. Rev. B , 1 (2017),arXiv:1703.04928. A. F. White and G. K.-L. Chan, J. Chem. Theory Comput. , 5690(2018). R. J. Bartlett, Annu. Rev. Phys. Chem. , 359 (1981). R. J. Bartlett and M. Musiał, Rev. Mod. Phys. , 291 (2007). A. Shee and D. Zgid, arXiv:1906.04079 (2019),arXiv:1906.04079. T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler, Rev. Mod.Phys. , 1027 (2005). G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcol-let, and C. A. Marianetti, Rev. Mod. Phys. , 865 (2006). K. Held, Adv. Phys. , 829 (2007). G. Biroli, O. Parcollet, and G. Kotliar, Phys. Rev. B , 205108(2004). M. Capone, M. Civelli, S. S. Kancharla, C. Castellani, andG. Kotliar, Phys. Rev. B , 195105 (2004). I. de Vega, U. Schollwöck, and F. A. Wolf, Phys. Rev. B ,155126 (2015). A. Liebsch and H. Ishida, J. Phys. Condens. Matter , 053201(2012). E. Koch, G. Sangiovanni, and O. Gunnarsson, Phys. Rev. B ,115102 (2008). F. A. Wolf, A. Go, I. P. McCulloch, A. J. Millis, and U. Scholl-wöck, Phys. Rev. X , 041032 (2015). R. Peters, Phys. Rev. B , 075139 (2011). M. Ganahl, P. Thunström, F. Verstraete, K. Held, and H. G. Evertz,Phys. Rev. B , 045144 (2014). M. Ganahl, M. Aichhorn, H. G. Evertz, P. Thunström, K. Held,and F. Verstraete, Phys. Rev. B , 155132 (2015). I. Shavitt and R. J. Bartlett,
Many-Body Methods in Chemistry and Physics (Cambridge University Press, Cambridge, 2009). G. E. Scuseria, C. L. Janssen, and H. F. Schaefer, J. Chem. Phys. , 7382 (1988). M. Nooijen and J. G. Snijders, Int. J. Quantum Chem. , 55(1992). M. Nooijen and J. G. Snijders, Int. J. Quantum Chem. , 15(1993). 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). K. Bhaskaran-Nair, K. Kowalski, and W. A. Shelton, J. Chem.Phys. , 144101 (2016). J. F. Stanton and R. J. Bartlett, J. Chem. Phys. , 7029 (1993). H. J. Monkhorst, Int. J. Quantum Chem. , 421 (2009). A. I. Krylov, Annu. Rev. Phys. Chem. , 433 (2008). V. M. Galitskii and A. B. Migdal, Sov. Phys. JETP (1958). M. F. Lange and T. C. Berkelbach, J. Chem. Theory Comput. ,4224 (2018). S. R. White and A. E. Feiguin, Phys. Rev. Lett. , 076401 (2004). E. Ronca, Z. Li, C. A. Jimenez-Hoyos, and G. K.-L. Chan, J.Chem. Theory Comput. , 5560 (2017). Y. Saad and M. H. Schultz, SIAM J. Sci. Stat. Comput. , 856(1986). E. de Sturler, SIAM J. Numer. Anal. , 864 (1999). P. Pulay, Chem. Phys. Lett. , 393 (1980). 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, Wiley Interdiscip. Rev. Comput. Mol. Sci. ,e1340 (2018). J. Hubbard, Proc. R. Soc. London. Ser. A. Math. Phys. Sci. ,238 (1963). A. Liebsch, H. Ishida, and J. Merino, Phys. Rev. B - Condens.Matter Mater. Phys. , 1 (2008), arXiv:0810.4837. A. Liebsch and N. H. Tong, Phys. Rev. B - Condens. Matter Mater.Phys. , 1 (2009). M. Rozenberg, G. Kotliar, and H. Kajueter, Phys. Rev. B - Con-dens. Matter Mater. Phys. , 8452 (1996), arXiv:9509182 [cond-mat]. C. J. Bolech, S. S. Kancharla, and G. Kotliar, Phys. Rev. B ,075110 (2003). E. H. Lieb and F. Y. Wu, Phys. Rev. Lett. , 1445 (1968). H. Park, K. Haule, and G. Kotliar, Phys. Rev. Lett. , 186403(2008). G. Sordi, K. Haule, and A.-M. S. Tremblay, Phys. Rev. Lett. ,226402 (2010). J. McClain, Q. Sun, G. K.-L. Chan, and T. C. Berkelbach, J. Chem.Theory Comput. , 1209 (2017). T. Gruber, K. Liao, T. Tsatsoulis, F. Hummel, and A. Grüneis,Phys. Rev. X8