Efficient variational approach to the impurity problem and its application to the dynamical mean-field theory
aa r X i v : . [ c ond - m a t . s t r- e l ] J u l An efficient variational approach to the impurity problem and its application to thedynamical mean field theory
Chungwei Lin and Alexander A. Demkov ∗ Department of Physics, University of Texas at Austin (Dated: July 27, 2018)Within the framework of exact diagonalization (ED), we compute the ground state of Andersonimpurity problem using the variational approach based on the configuration interaction (CI) expan-sion. We demonstrate that an accurate ground state can be obtained by iteratively diagonalizing amatrix with the dimension that is less than 10% of the full Hamiltonian. The efficiency of the CIexpansion for different problems is analyzed. By way of example, we apply this method to the single-site dynamical mean field theory using ED as the impurity solver. Specifically, to demonstrate theusefulness of this approach, we solve the attractive Hubbard model in the grand-canonical ensemble,where the s-wave superconducting solution is explicitly obtained.
PACS numbers: 71.10.Fd, 71.30.+h, 71.20.b, 71.15.m
I. INTRODUCTION
The electron-electron interactions are of fundamental importance in wide classes of systems and account for manyfascinating phenomena, such as the metal insulator transition [1], superconductivity [2], and magnetic phase [3, 4]. Inprinciple, for a system whose orbitals near the Fermi energy are spatially localized, such as the transition metal oxideswith occupied d orbitals [5], the electron-electron interaction is important and cannot be neglected. When takingthe inter-particle interaction into account, the problem becomes “strongly correlated” and is very difficult to solve,because an exponentially large number of states couples to one another and a direct determination of the eigenstatesbecomes computationally impossible. There are several well-developed approaches to the strongly correlated problem.Almost all of them, in one way or another, involve identifying a reduced space of the workable dimension that includes(or at least overlaps strongly with) the true ground state. For example, the Hartree-Fock approximation and densityfunctional theory [6, 7] map the correlated problem onto a non-interacting one, with the dimension depending linearlyon the system size. The renormalization approach keeps only low-excited states [8, 9] in momentum space, or onlythe most relevant states [10, 11] in real space. Using a variational wave function can greatly reduce the dimension,but this approach is system specific and requires insight into the physics of the problem. One of the most successfulcases is Laughlin’s variational wave function for the two-dimensional electron gas under a strong magnetic field [12].If only the ground state properties of the system are of interest, the recently proposed density matrix embeddingtheory requires solving only a two-site problem [13].The dynamical mean field theory (DMFT) [14–18] has emerged as a powerful tool to study the strongly correlatedproblem. The key step of DMFT is to map the lattice problem onto a “simpler” impurity model [19], where a finite set ∗ E-mail: [email protected] of interacting impurity orbitals is coupled to a large number of noninteracting bath orbitals. We will demonstrate inSection III that the impurity problem is simpler in the sense that it requires a much smaller number of determinantalstates to capture the ground state than what is needed in the lattice model. An accurate numerical solver for correlatedimpurity models is crucial for DMFT calculations, and several exact and numerically tractable approaches have beendeveloped. The frequently used solvers include the continuous-time quantum Monte Carlo [18], renormalization [9, 20],and exact diagonalization (ED) [21, 22].Recently, Zgid et al . have applied the truncated configuration interaction (CI) expansion, a technique developed inquantum chemistry [23, 24], to the Anderson impurity model [25] and DMFT [26, 27]. The CI solver can be viewedas an approximation to the full ED, as it finds the ground state in the truncated Hilbert space (HS) based on energyconsiderations. Zgid et al . demonstrate that only a small number of determinantal states is needed to obtain theaccurate ground state of the Anderson impurity model. This finding makes the CI solver a strong candidate forstudying more complicated problems such as multi-orbital models, cluster DMFT [27], or even the non-equilibriumproblem [28]. Very recently, the CI method has been applied to the density matrix renormalization group method[29].The CI algorithm and its application to DMFT have been discussed in detail in Ref. [26, 27]. In this paper,we further analyze the efficiency of the CI method based on the variational principle. In particular, we formulatean efficient and general procedure, based on the successive orbital transformation, to accurately express the groundstate in a small number of determinantal states, from which the Green’s function can be calculated with minimumcomputational resources. The same procedure determines whether the CI truncation is valid or not.The rest of the paper is organized as follows. In Section II we briefly introduce the CI method, and describe severalCI truncation schemes. We define the correlation entropy to quantify the degree of correlation. In Section III wesolve an impurity model and an 8-site Hubbard model to demonstrate under which circumstances the CI solver canbe efficient. A general strategy to compute the ground state and Green’s function is formulated. In Section IV wesolve the DMFT equation for the Hubbard model with attraction in the grand-canonical ensemble to describe thesuperconducting state. In the appendices we provide some key steps of successive orbital transformation, and detailsof solving the impurity problem in the grand-canonical ensemble.
II. CONFIGURATION INTERACTION EXPANSIONA. Overview and basis classification
Using the ED to compute the ground state of a given Hamiltonian requires two basic ingredients: first, all basisvectors of the HS (or its invariant subspace) have to be labeled, and second, the matrix elements between any two basisvectors have to be computed. Once these are available, the ground state can be computed by directly diagonalizing theHamiltonian matrix, or by a power method such as Lanczos [2] or Davidson algorithms [30]. The CI solver [23, 24, 27]is a non-perturbative approximation to solve a many-body problem based on the ED formalism. In essence, it is atruncation of HS based on the energy considerations.The many-body fermionic HS is constructed from the antisymmetrized single-particle orbitals (determinantalstates), which are concisely expressed in the second quantization notation (Fock states). For n e fermions occu- FIG. 1: (Color online) (a) Basis states for 4 fermions occupying 8 spin-orbitals. Picking up one determinantal state as thereference, all basis states can be classified by the number of particle-hole pairs with respect to the reference state. (b) Sixdeterminantal states constructed from CAS(2,4). Two secondary and active orbitals are always empty and filled respectively.(c) Some allowed and truncated states in CAS(2,4)CIS or equivalently RAS(2,-1;2,1). The inactive orbitals can have at mostone hole, while secondary orbitals can at most have one particle, therefore the last two states are not included. pying n o spin-orbitals, the HS dimension D H is C n o n e ≡ n o ! n e !( n o − n e )! . As an example, we consider 4 fermions occupying8 spin-orbitals, with each orbital specified by a creation operator d † j [ j =1,2,..., n o ( n o = 8)]. The 4-particle Hilbertspace contains D H = C = 70 states, with each basis state labeled by | i i [ e.g. , | i i = d † d † d † d † | i ] with i = 1 , , ..., | GS i = D H (=70) X i =1 c i | i i . (1)To avoid the confusion in terminology, throughout this paper the “vectors” or “states” are reserved for the many-body wave functions; while “orbitals” are reserved for the single-particle wave functions (specified by { d † i } ) that are usedto build the many-body states.By choosing one determinantal state as a reference, all basis states can be classified by the number of particle-hole(p-h) pairs with respect to the reference state [31]. In our example shown in Fig. 1 (a), there are 16 [( C ) ] basisstates with one and three p-h pairs, 36 [( C ) ] with two p-h pairs, and 1 [( C ) ] with zero and four p-h pairs. Thetruncated CI expansion only keeps the basis states with a small number of p-h pairs. When keeping only the referencestate, the CI method is equivalent to the Hartree-Fock approximation. Keeping states up to 1, 2, 3, and 4 p-h pairsare referred to as CIS, CISD, CISDT, and CISDTQ expansions (S: single, D: double, T: triple, Q: quadruple), andthe corresponding truncated spaces are referred to as the CI subspaces [23, 24]. The non-reference states (stateswith non-zero p-h pairs) are sometimes referred to as “excitation levels” in the Hartree-Fock sense, but they are notactually the excitations of the Hamiltonian and we avoid using this term here.The efficiency of the CI truncation depends crucially on the choice of the single-particle orbitals. As a trivialexample to illustrate its importance, we note that the many-body ground state of a non-interacting Hamiltonian canbe described by a single determinantal state – if choosing the orbitals diagonalizing the single-particle Hamiltonian,all non-reference states can be discarded; otherwise more determinantal states are needed. B. Multi-reference configurations, orbital classification and truncation schemes
When several determinantal states simultaneously have a significant weight in the ground state, the CI expansionfrom a single reference may become inadequate, and a generalization to multi-reference configurations is needed.To deal with this situation, the orbitals are classified according to their occupation: they are (1) inactive (alwaysoccupied); (2) secondary (always empty); or (3) active (partially filled). Building the CI space based on these threeclasses of orbitals is referred to as a “complete active space” (CAS) approximation [23]. The dimension of CAS isdetermined by the number of electrons in active orbitals only. The notation CAS( m, n A ) indicates filling m electronsin n A active orbitals. For our example, six determinantal states built from CAS(2,4) (2 electrons occupying 4 activeorbitals) are shown in Fig. 1 (b). The CAS is designed to unbiasedly capture contributions from determinantal statesmade of inactive and active orbitals, which serve as the multi-reference states.One can further expand the CI space by including states of p-h pairs with respect to multi-reference states definedby CAS. Following the notation in Ref. [27], the expanded CI spaces are denoted by CAS( m, n A )CIS [single p-hpair upon CAS( m, n A )], CAS( m, n A )CISD [single and double p-h pairs upon CAS( m, n A )],.. etc. The generalizationis equivalent to the “restricted active space” (RAS) approximation [23]. In RAS, the inactive orbitals are mostlyoccupied but are allowed to have small number of holes; the secondary orbitals are mostly empty but are allowed tohave a small number of particles; the active orbitals have no constraints. We use the notation RAS( n I , − k ; n S , l ) toindicate allowing maximum k holes (we use minus sign to indicate the holes) in n I inactive orbitals, and maximum l particles in n S secondary orbitals. Typically RAS( n I ,-1; n S ,1) = CAS( m, n A )CIS, with n I and n S are respectivelythe number of inactive and secondary orbitals.Note that in CAS( m, n A ), we specify the the number of active orbitals whereas in RAS( n I , − k ; n S , l ) we specifyinactive and secondary orbitals, and n I + n S + n A = n o . The single-reference expansion means no active orbitals. InTable I we summarize CI truncation schemes based on electron distributions over different orbital classes.TABLE I: Electron distributions over different orbital classes for CI truncation schemes.Scheme inactive secondary activeCIS(D,T,Q) maximum 1 (2,3,4) hole maximum 1 (2,3,4) electron –CAS( m, n A ) 0 hole 0 electron no constraintRAS( n I , − k ; n S , l ) maximum k hole maximum l electron no constraintIn terms of the variational principle, different truncation schemes correspond to different subspaces, and one finds thebest single-particle orbitals in the restricted space. C. Natural orbitals, degree of correlation, and iteration to ground state
If the many-body ground state | GS i is known, one can determine the “natural orbital” basis which leads to themost rapid convergence of the CI calculation. Natural orbitals are eigen basis of the one-particle density matrix D [ { d } ] = h GS | d † i d j | GS i , (2)with operators { d } specifying the orbitals and i, j = 1 , , ...n o . The eigenvalues of D give the natural occupancies.When only one determinantal state is needed (non-interacting Hamiltonian), the natural occupancies are either 1 or0. Any natural occupancy deviating from 0 or 1 indicates the need of more than one determinantal state to properlydescribe the ground state.The degree of correlation can be thus characterized by the minimum number of determinantal states needed todescribe the ground state, which is directly related to the natural occupancies. Similar to the entanglement entropy[32, 33], we can define the “correlation entropy” as S = n o X i =1 S i = − n o X i =1 p i log p i , (3)to measure the degree of correlation. Here p i is the occupation of i th natural orbital . This definition ensures thatwhen p i = 0 or 1, the resulting entropy is zero. Therefore the larger the correlation entropy, the more determinantalstates are needed, and the more correlated the ground state is.Since the ground state is not known until the calculation is performed, in practice the natural orbitals are determinedin the following way. Starting from an orbital basis { d (0) } (can be arbitrary), the ground state is determined undersome CI truncation scheme, from which the single-particle density matrix ( D (0) ) and the natural orbitals ( { d (1) } which diagonalize D (0) ) are determined. This procedure can be iterated as follows: { d (0) } ⇒ | GS (0) i ⇒ D (0) = h GS (0) | ( d (0) i ) † d (0) j | GS (0) i , d (1) i = X j U (1) i,j d (0) j { d (1) } ⇒ | GS (1) i ⇒ D (1) = h GS (1) | ( d (1) i ) † d (1) j | GS (1) i , d (2) i = X j U (2) i,j d (1) j ........ { d ( p ) } ⇒ | GS ( p ) i ⇒ D ( p ) = h GS ( p ) | ( d ( p ) i ) † d ( p ) j | GS ( p ) i , (4)where the { d ( i +1) } ( d ( i +1) i = P j U ( i +1) i,j d ( i ) j ) orbitals diagonalize D ( i ) . Some key steps of solving the problem duringthis successive orbital transformation is given in Appendix A.Eq. (4) is useful in computing the ground state based on the following observations. First, we observe that theiteration converges (see Section III B and D), and the converged ground state corresponds to a variational solution(local energy minimum), subject to the given truncation scheme. This result can be understood as follows. Thenatural orbital basis is the single-particle basis which requires the least number of determinantal states to describethe full | GS i , therefore the CI method performs best in this basis. Indeed, we find that during the iterations definedby Eq. (4), more and more natural orbital occupations approach 0 or 1 (justifying the use of the secondary andinactive orbitals), and the ground state energy is lowered (see Table III and V).The second observation is that when the CI schemes allow all orbitals to mix during iterations, the convergedsolution is the global minimum subject to the CI scheme. This is an empirical statement and numerical evidence willbe provided in the next section. This observation implies that a CAS scheme is not suitable for finding the globalminimum because during iterations, secondary, active, and inactive orbitals do not mix [23]. The CIS(D,T) and RASschemes, however, have no problem finding the global minimum. Therefore, if the ground state falls into the subspacedefined by CIS(D,T) or RAS schemes, one can find the full ground state in the restricted subspace [34]. D. Limitations of the CI truncation
We conclude this section by discussing the limitations of the CI truncation. When none of the natural occupationsapproach 0 or 1 (within 10 − ), the CI method does not have any meaningful advantage, because the exact groundstate involves too many determinantal states, no matter which orbital basis one uses. This limit intrinsically originatesfrom the strong correlation of the ground state (or the Hamiltonian), where an expansion based on the determinantalbasis is not a good starting point. Examples will be provided in the next section. III. APPLICATION TO A FINITE-SIZE MODELA. Overview
In this section we provide two examples demonstrating the efficiency of the CI method. The efficiency, for ourpurpose, is quantified by the number of states needed during the calculation. The first example is the Andersonimpurity model, where the on-site repulsion is present only at one (impurity) site. In this case, the CI methodperforms well – including only a tiny fraction of the total HS gives the ground state with high accuracy. The secondexample is a one-dimensional Hubbard model where on-site repulsion is present at all sites. In this case the CI methodperforms poorly – essentially all states in the Hilbert space are needed to get the full ground state.We consider 8 electrons and 8 spatial orbitals (16 spin-orbitals) where the full ED calculation can be performed withease. The CI schemes are formulated in spin-orbitals, i.e. same spatial orbitals of opposite spins have different labelswhen specifying inactive ( n I ), active ( n A ), and secondary ( n S ) orbitals. In the current case n A + n I + n S = n o = 16.Using spin-orbitals facilitates describing the superconducting state, which does not conserve the total number ofparticles (see Appendix B). The ground state is searched in the M z = 0 ( M z is the z -component of total spin)invariant subspace of dimension 4900 [( C ) ], and is obtained using Davidson algorithm [30]. B. Anderson impurity model
The Hamiltonian of the Anderson impurity model is H = ε d X σ c † ,σ c ,σ + U n , ↑ n , ↓ + X p =2 X σ t p [ c † ,σ c p,σ + h.c. ] + X p =2 X σ ε p c † p,σ c p,σ (5)We fix the parameters ε p = ± , ± / , ± / , t p = 1 / √
8, and vary U with ε d = − U/ U = 20 , , , ,
1. The degenerate occupancies are due to spin symmetry.TABLE II: Anderson impurity model. Natural orbital occupancies for U = 20 , , , , µ = U/ n i − = n i ( i = 1 , , ..., U n , n n , n n , n n , n n , n n , n n , n n , n
20 1 0.999999 0.999566 0.582742 0.417258 0.000434141 8.28594e-07 3.08453e-1010 1 0.999999 0.999197 0.67482 0.32518 0.000802745 1.45276e-06 1.03073e-095 1 0.999999 0.999331 0.81756 0.18244 0.000669318 1.25351e-06 1.06296e-092 1 1 0.999826 0.954373 0.0456267 0.000174072 3.46593e-07 2.57253e-101 1 1 0.999954 0.987583 0.0124174 4.64471e-05 9.20178e-08 8.14356e-11From Table II we see that for a wide range of U values, out of 16 spin-orbitals, only eight orbitals – from 5 to 12 – areactive (partially occupied), and one in principle only needs determinantal states composed of 8 orbitals to accuratelydescribe the ground state. The intermediate U = 10 needs the largest number of determinantal states, we thereforeapply several CI schemes for U = 10 to Eq. (4). The CI schemes, their dimensions, and the ground state energy ateach iteration are given in Table III. The starting orbital basis is chosen as the local basis which is very different fromthe natural basis. We have tried several starting orbital basis, and they all converge to the same ground state energy.TABLE III: Calculated ground state energy using iterated orbital basis for different CI schemes.These results are for U = 10 and ε d = −
5. The number in parenthesis indicates the dimensionof the CI subspace. “-” indicates the same energy as the previous iteration.iteration CISD (361) CISDT (1545) RAS(4,-1; 4,1) (523) RAS(6,-1; 6,1) (118) Exact (4900)0 -13.2455 -13.2465 -13.1178 -13.1167 -13.24651 -13.2464 - -13.2465 -13.2463 -2 - - - -13.2464 -3 - - - - -We can point out some general features. First, the ground state energy decreases during iterations, indicating thatit converges to at least a local minimum in the subspace. Moreover, the converged ground state energy is very closeto the exact one, indicating the solution is the global minimum in the subspace (the statements in Section II.C).Depending on the required accuracy, all four CI schemes give good ground state energy with the error smaller than10 − . Typically, RAS schemes perform better than a single-reference CI. Once the natural orbital basis is obtained,we can perform the calculation again, using the CAS scheme based on the natural orbital occupation. We foundthat performing CAS(4,8), which contains only 35 determinantal states ( C plus M z = 0 constraint), after obtainingnatural orbitals using any schemes listed in Table III gives the practically exact ground state energy of -13.2465. Thismeans that one in principle only needs to keep 118 states [RAS(6,-1; 6,1) in Table III] out of 4900 states during theentire iteration to get the ground state of high accuracy, although in practice one needs to check the convergenceagainst other CI schemes keeping more states. Based on the natural occupancy analysis (Table II), 35 is the minimum number of determinantal basis states needed. C. Strategy to compute the ground state and Green’s functions
The impurity Green’s function is defined as G imp ( iω n ) = h GS | c [ iω n − ( H − E GS )] − c † | GS i + h GS | c † [ iω n + ( H − E GS )] − c | GS i . (6)with H defined in Eq. (5) and E GS the ground state energy. To compute G imp , one needs to specify the HS of n e ± n e ± • Use RAS to obtain the natural orbital basis. To save time, one starts by using severely truncated scheme, thenchecks the convergence by comparing the ground state energy obtained from the less-truncated CI schemes. • Once the natural orbital basis is determined, perform a CAS calculation [according to the natural occupancies]to further reduce the dimension of the subspace without losing accuracy. The reason is that the corresponding n e ± • With the ground state in the CAS scheme [35], compute the Green’s function using the Lanczos algorithm[2, 36].
D. Hubbard model on a ring
We test this strategy by considering a Hubbard model on an 8-site ring described by the following Hamiltonian: H = − µ X i,σ c † i,σ c i,σ + U X i,σ n i, ↑ n i, ↓ − t X i =1 ,σ [ c † ,σ c i +1 ,σ + h.c. ] , (7)with the periodic boundary condition c ≡ c . We take t = 1 for different values of U at half filling ( µ = U/ U = 20 , , , U = 20 , , , µ = U/ n i − = n i ( i = 1 , , ..., U n , n n , n n , n n , n n , n n , n n , n n , n
20 0.636333 0.59957 0.59957 0.5 0.5 0.40043 0.40043 0.36366710 0.748239 0.695164 0.695164 0.5 0.5 0.304836 0.304836 0.2517614 0.914041 0.882618 0.882616 0.500001 0.499999 0.117384 0.117382 0.08595912 0.974145 0.963764 0.963764 0.500002 0.499998 0.0362365 0.0362365 0.0258546From Table IV, we see that the natural occupancies are not close to either 0 or 1, even for U = 2. From thediscussion in Section II.C and D, we conclude that the CI method is not suitable for this situation. For U = 2 case,we perform CISDTQ (=CI[4 p-h]), CI[5 p-h], CI[6 p-h] calculations (the number in square brackets indicates themaximum number of p-h pairs), and find that the efficiency of the CI method is poor – to converge to 10 − in energy,the required number of determinantal states is of the order of the original HS. The results are summarized in TableV. One notices that the iteration still finds the global minimum in the restricted space.TABLE V: Calculated ground state energy using iterated orbital basis for different CI schemes.These results are for U = 2 and µ = U/
2. The number in parenthesis indicates the dimensionof the CI space. “-” indicates the same energy as the previous iteration.iteration CISDTQ (3355) CI[5 p-h] (4539) CI[6 p-h] (4867) Exact (4900)0 -14.1328 -14.5148 -14.5681 -14.56821 -14.5581 -14.5642 -14.5682 -2 -14.5625 - - -3 - - - -As was mentioned in the introduction, the DMFT maps a lattice problem onto an impurity one. By comparing twoexamples in Section III, the lattice model is intrinsically more complicated because its ground state involves manymore determinantal states. This provides a quantitative explanation why the impurity model is more numericallytractable than the lattice model. Within DMFT, our two model calculations suggest that the main advantage of CImethod over the full ED is the ability to include more bath orbitals.
IV. APPLICATION TO DYNAMICAL MEAN FIELD THEORYA. Overview
We now apply the CI method introduced in Section II to the DMFT, using ED as the impurity solver. TheDMFT involves self-consistently determining the parameters of an impurity problem, and CI method enters as anapproximation to the ED solver. The impurity Green’s function, which is the most resource-consuming part of DMFTequations (in terms of time and computer memory), is computed using the strategy formulated in Section III C. Wechoose to solve the Hubbard model with attraction on the infinite-dimension Bathe lattice (bandwidth 4 t ), under thes-wave superconducting self-consistent condition. This model has been studied in Refs [37–39]. The main complicationarising from the superconducting state is that the total particle number is not conserved anymore, and one has toformulate the problem in the grand-canonical ensemble including particle number fluctuations. The grand-canonicalensemble enlarges the HS dimension and complicates the calculation, which we choose as a platform to test theefficiency of the CI method.We solve the problem at half and quarter fillings. We are able to consider up to 9 bath sites (totally 10 sites includingthe impurity) using moderate computational resources (a single-core laptop with 2GB memory). The convergencecan be accelerated by using severely truncated CI schemes in the first few DMFT self-consistency iterations, as theseiterations improve the bath parameters at a relatively low computational cost. We found that the results for 5 bathsites are essentially identical to those with 7 and 9 bath sites, with the exception of the spectral functions [39].For spectral functions, however, the gap (or more precisely the reduction of the density of states) around the Fermienergy is robust for 5, 7, and 9 bath sites. As for the ED solver, we use the fictitious temperature of 0.1 t , andLevenberg-Marquardt algorithm [40] to determine the bath parameters. In the next subsection we will formulate theself-consistent equation and show the results. Details of solving the impurity model in the grand-canonical ensembleare given in the Appendix B.0 B. Self-consistent equation
The self-consistent equation for the s-wave superconducting state has been derived in Refs [14, 21]. To describe thesuperconducting state, one uses the Nambu Green’s function which breaks the U(1) symmetry. By defining a Nambuspinor as Ψ † k = ( c † k , ↑ , c − k , ↓ ), the lattice Green function (a 2 × G ( k , τ ) = − T h Ψ k ( τ )Ψ † k (0) i = − T h c k , ↑ ( τ ) c †− k , ↓ ( τ ) (cid:16) c † k , ↑ (0) c − k , ↓ (0) (cid:17) i = − T h c k , ↑ ( τ ) c † k , ↑ (0) i − T h c k , ↑ ( τ ) c − k , ↓ (0) i− T h c †− k , ↓ ( τ ) c † k , ↑ (0) i − T h c †− k , ↓ ( τ ) c − k , ↓ (0) i ≡ G ( k , τ ) F ( k , τ ) F ( − k , − τ ) ∗ − G ( − k , − τ ) = G ( k , τ ) F ( k , τ ) F ( k , τ ) ∗ − G ( − k , − τ ) (8)where h ... i represents the ground state expectation value. The matrix structure of ˆ G ( k , τ ) is obtained by using time-translational symmetry and anti-communication relation. For example, two off-diagonal terms are related (neglectingthe spin and momentum subscripts): h c ( τ ) c (0) i ≡ f ( τ ) ⇒ h c † ( τ ) c † (0) i = h [ c (0) c ( τ )] † i = h [ c ( − τ ) c (0)] † i = f ∗ ( − τ ) . For the third line in Eq. (8), we further use the property that for a singlet pairing field, F ( − k , − τ ) = F ( k , τ ) [fortriplet pairing (p-wave), one instead obtains F ( − k , − τ ) = − F ( k , τ )] [14]. Under the single-site DMFT approximation,where self energy has no momentum dependence, the lattice Green’s function at Matsubara frequencies isˆ G ( k , iω n ) − = iω n − ( ε k − µ ) 00 iω n + ( ε − k − µ ) − ˆΣ( iω n ) , (9)with ˆΣ( iω n ) = Σ( iω n ) S ( iω n ) S ( iω n ) − Σ( iω n ) ∗ computed from the impurity problem. The local Green’s functionˆ G loc ( iω n ) = 1 N X k ˆ G ( k , iω n ) = G ( iω n ) F ( iω n ) F ( iω n ) − G ( − iω n ) , (10)is required to be the impurity Green’s function by the self-consistent condition.The pairing order parameter is defined as φ = T P n F ( iω n ). For any negative U , | φ | is non-zero [39]. When | U | issmall, the order parameter is exponentially small and cannot be resolved in our calculation. The calculated spectralfunction at quarter filling, using 5 bath sites for U = − , − , − | U | becomes larger. This corresponds to the formationof the superconducting gap and its quasi-particle excitations. Using the ED solver, however, we cannot observe thegap of too small amplitudes due to the finite bath discretization. In the inset we show the pairing order parameteras a function of U . Fig. 2 (b) shows the spectral function at half-filling for U = − FIG. 2: (Color online) (a) The spectral function of the spin-up electron for U = − , − , − | U | increases. (Inset) The corresponding order parameter φ as a function of U . (b) Thespectral function of the spin-up electron for U = − at the imaginary frequencies, as well as the order parameter, are well converged (the relative change | ∆ φ | / | φ | in theorder parameters is smaller than 10 − ). From Fig. 2 (b), we see that the spectral functions around the Fermi level,especially the gap structure, do not change as the bath size increases. Finally, we mention (not shown) that the sameself-consistent equations can apply to the positive U regime. There the pairing order parameter vanishes, and knownresults such as the metal-insulator transition are recovered. V. CONCLUSIONS
In this paper, we analyze the efficiency of the configuration interaction method as an approximation to the fullexact diagonalization. Not surprisingly, the CI method performs well if the exact ground state can be described by asmall number of determinantal states; but does not hold any meaningful advantage (compared to full ED) when theground state is composed of many determinantal states. In terms of the variational principle, the validity of the CItruncation depends intrinsically on the overlap between the exact ground state and the CI truncated subspace. Theperformance can be estimated by examining the natural occupancies computed by a truncated CI scheme (not thefull ED). If there are natural orbitals with occupations close to 1 and 0, the CI method is useful, otherwise it is not.We illustrate this point with two examples, the Anderson impurity model and the Hubbard model. When the CItruncation is valid, we demonstrate that by iteratively searching for the natural orbital basis under the CI schemes2that allow all orbitals to mix (such as CISD or RAS, but not
CAS), a very accurate ground state can be obtained.This procedure can be viewed as searching for the ground state in the variational subspace. Once the natural orbitalsare determined, a CAS can be performed to express the ground state in the minimum number of determinantal stateswithout losing accuracy. The impurity Green’s function can then be efficiently computed from the CAS ground state.As an illustration, we apply this method to the single-site dynamical mean field theory using ED as the impuritysolver. We solve the Hubbard model with attraction in the grand-canonical ensemble ( i.e. including the particlenumber fluctuation) for the s-wave superconducting state. All expected results, such as non-zero pairing and gapopening, are reproduced. Moreover, we are able to check the convergence with respect to the bath size (up to 9 bathsites) using only moderate computational resources (a single-core laptop of 2GB memory).Our two model calculations, Anderson impurity and ring Hubbard model, suggest that CI method is intrinsicallyefficient for the impurity model (only a few orbitals are correlated), but not so for the lattice model (all orbitalsare correlated). Based on this observation, we think the CI method is most useful, in terms of studying correlatedsystems, when combining with DMFT. As a DMFT impurity solver, the main advantage of CI method over the fullED is its ability to include more bath orbitals. This is consistent with our superconducting calculation and the clustercalculation done in Ref. [27]. The CI solver, by itself, is also efficient when the impurity problem is far away from thehalf filling.
Acknowledgements
C.L. thanks Hung The Dang, Ara Go, Andrew Millis, Takashi Oka, Hsiang-Hsuan Hung, and Andreas R¨uegg forhelpful discussions. We thank Richard Hatch, Hosung Seo, and Alex Slepko for insightful comments. Support forthis work was provided through Scientific Discovery through Advanced Computing (SciDAC) program funded by U.S.Department of Energy, Office of Science, Advanced Scientific Computing Research and Basic Energy Sciences underaward number DESC0008877. [1] M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. , 1039 (1998), URL http://link.aps.org/doi/10.1103/RevModPhys.70.1039 .[2] E. Dagotto, Rev. Mod. Phys. , 763 (1994).[3] H. v. L¨ohneysen, A. Rosch, M. Vojta, and P. W¨olfle, Rev. Mod. Phys. , 1015 (2007), URL http://link.aps.org/doi/10.1103/RevModPhys.79.1015 .[4] S. Bl¨ugel, Phys. Rev. Lett. , 851 (1992), URL http://link.aps.org/doi/10.1103/PhysRevLett.68.851 .[5] J. Hubbard, Proc. R. Soc. Lond. A , 238 (1963).[6] P. Hohenberg and W. Kohn, Phys. Rev. , B864 (1964), URL http://link.aps.org/doi/10.1103/PhysRev.136.B864 .[7] W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965), URL http://link.aps.org/doi/10.1103/PhysRev.140.A1133 .[8] K. G. Wilson, Rev. Mod. Phys. , 773 (1975), URL http://link.aps.org/doi/10.1103/RevModPhys.47.773 .[9] R. Bulla, T. A. Costi, and T. Pruschke, Rev. Mod. Phys. , 395 (2008), URL http://link.aps.org/doi/10.1103/RevModPhys.80.395 .[10] S. R. White, Phys. Rev. B , 10345 (1993), URL http://link.aps.org/doi/10.1103/PhysRevB.48.10345 . [11] U. Schollw¨ock, Rev. Mod. Phys. , 259 (2005), URL http://link.aps.org/doi/10.1103/RevModPhys.77.259 .[12] R. B. Laughlin, Phys. Rev. Lett. , 1395 (1983), URL http://link.aps.org/doi/10.1103/PhysRevLett.50.1395 .[13] G. Knizia and G. K.-L. Chan, Phys. Rev. Lett. , 186404 (2012), URL http://link.aps.org/doi/10.1103/PhysRevLett.109.186404 .[14] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. , 13 (1996).[15] T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler, Rev. Mod. Phys. , 1027 (2005).[16] G. Kotliar and D. Vollhardt, Physics Today (2004).[17] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Rev. Mod. Phys. , 865(2006).[18] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Rev. Mod. Phys. , 349 (2011), URL http://link.aps.org/doi/10.1103/RevModPhys.83.349 .[19] P. W. Anderson, Phys. Rev. , 41 (1961), URL http://link.aps.org/doi/10.1103/PhysRev.124.41 .[20] D. J. Garcia, K. Hallberg, and M. J. Rozenberg, Phys. Rev. Lett. , 246403 (2004), URL http://link.aps.org/doi/10.1103/PhysRevLett.93.246403 .[21] M. Caffarel and W. Krauth, Phys. Rev. Lett. , 1545 (1994), URL http://link.aps.org/doi/10.1103/PhysRevLett.72.1545 .[22] A. Liebsch and I. Ishida, Journal of Physics: Condensed Matter (2012).[23] T. Helgaker, P. Jørgensen, and J. Olsen, Molecular Electronic-Structure Theory (Job Wiley & Sons, LTD., 2000).[24] C. D. Sherrill and H. F. Schaefer III, in The Configuration Interaction Method: Advances in Highly Correlated Approaches,Advances in Quantum Chemistry, Vol. 34, (edited by P. O. Lowdin, J. R. Sabin, M. C. Zerner, and E. Brandas, Academic,New York, pp. 143-269, 1999).[25] O. Gunnarsson and K. Sch¨onhammer, Phys. Rev. B , 4315 (1983), URL http://link.aps.org/doi/10.1103/PhysRevB.28.4315 .[26] D. Zgid and G. K.-L. Chan, J. Chem. Phys. , 094115 (2011).[27] D. Zgid, E. Gull, and G. K.-L. Chan, Phys. Rev. B , 165128 (2012).[28] E. Arrigoni, M. Knap, and W. von der Linden, Phys. Rev. Lett. , 086403 (2013), URL http://link.aps.org/doi/10.1103/PhysRevLett.110.086403 .[29] Y. Ma and H. Ma, arXiv:1303.0616 (2013).[30] E. R. Davidson, Journal of Computational Physics , 87 (1975), ISSN 0021-9991.[31] Note that when mentioning states of non-zero particle-hole pairs, a reference state is already (implicitly) defined.[32] I. Peschel, Brazilian J. of Phys. , 267 (2012).[33] K. Byczuk, J. Kuneˇs, W. Hofstetter, and D. Vollhardt, Phys. Rev. Lett. , 087004 (2012), URL http://link.aps.org/doi/10.1103/PhysRevLett.108.087004 .[34] R. S. Grev and Schaefer, J. Chem. Phys. , 6850 (1992).[35] If the CAS space is inadequate because no natural occupancies approach close enough 1 and 0, one then uses the RASspace of the smallest dimension.[36] G. Grosso and G. P. Parravicini, Solid State Physics (Acedemic Press, 2000).[37] M. Keller, W. Metzner, and U. Schollw¨ock, Phys. Rev. Lett. , 4612 (2001), URL http://link.aps.org/doi/10.1103/PhysRevLett.86.4612 .[38] M. Capone, C. Castellani, and M. Grilli, Phys. Rev. Lett. , 126403 (2002), URL http://link.aps.org/doi/10.1103/PhysRevLett.88.126403 .[39] A. Toschi, P. Barone, M. Capone, and C. Castellani, New J. Phys. , 7 (2005), URL http://stacks.iop.org/1367-2630/7/i=1/a=007 .[40] We use minpack package. [41] P. B. Allen and B. Mitrovic, Solid State Physics , 1 (1983). Appendix A: Successive orbital transformations
We assume the Hamiltonian is specified in the basis { d (0) } as H ( { d (0) } ). Two sets of basis orbitals, { d (0) } and { d ( p ) } , are related by an unitary matrix d ( p ) i = X j ˜ U ( p ) i,j d (0) j . (A1)Once ˜ U ( p ) is known, one can express H ( { d (0) } ) in terms of H ( { d ( p ) } ), and re-diagonalize the problem in the { d ( p ) } basis. As shown in Eq. (4), ˜ U ( p ) can be obtained by˜ U ( p ) = U ( p ) ...U (2) U (1) ≡ Π pi =1 U ( i ) = U ( p ) Π p − i =1 U ( i ) = U ( p ) ˜ U ( p − . (A2)Note that in our notation, U ( p ) connects { d ( p − } and { d ( p ) } basis sets, whereas ˜ U ( p ) connects { d (0) } and { d ( p ) } basissets. For the first iteration, ˜ U (0) can be arbitrary. At the iteration using { d ( p ) } basis set [see Eq. (4)], ˜ U ( p − isneeded. After getting the CI ground state | GS ( p ) i , we diagonalize the single particle density matrix D ( p ) and get U ( p ) , from which ˜ U ( p ) = U ( p ) ˜ U ( p − is obtained for the next iteration using the { d ( p +1) } basis set.We would like to point out an important step concerning the evaluation of the single-particle density matrix in thetruncated CI space [ D ( p ) in Eq. (4)]. Since the CI ground state | GS ( p ) i does not include all n e -particle Hilbert space,the state ( d ( p ) i ) † d ( p ) j | GS ( p ) i generally leads to components outside the truncated CI space. This numerical error canlead to a non-Hermitian D ( p ) and therefore non-unitary U ( i ) . To avoid this error, we construct the ( n e − n e -particle CI space. For example, if the n e -particle CI groundstate is approximated in RAS( n I , − k ; n S , l ), then we construct RAS( n I , − k − n S , l ), which contains all componentsof d ( p ) j | GS ( p ) i . The D ( p ) is obtained by computing the inner-product of d ( p ) i | GS ( p ) i and d ( p ) j | GS ( p ) i .In the Green’s function calculation, we need to construct both RAS( n I , − k − n S , l ) [an ( n e − d ( p ) i | GS ( p ) i components] and RAS( n I , − k ; n S , l + 1) [an ( n e + 1)-particle space containing all ( d ( p ) i ) † | GS ( p ) i components], to completely describe hole and particle excitations. Once they are constructed, the corresponding holeand particle Krylov spaces respectively spanned by {| ψ ′ i = c | GS ( p ) i , | ψ ′ i = H | ψ ′ i , ..., | ψ ′ i +1 i = H | ψ ′ i i , ... } , {| ψ i = c † | GS ( p ) i , | ψ i = H | ψ i , ..., | ψ i +1 i = H | ψ i i , ... } (A3)are iteratively generated within the ( n e ±
1) subspace, in which the Green’s function is computed using the Lanczosalgorithm. Certainly the smaller n e -particle CI space leads to the smaller corresponding ( n e ± Appendix B: Grand canonical ensemble and Bogoliubov particles
Here we formulate the impurity model using the grand canonical ensemble that can describe the superconductingstate. To describe the superconducting state, the impurity model is generalized to [compare to Eq. (5)] H And,SC = ε d X σ c † ,σ c ,σ + N X p =1 t sc,p ( c † p, ↑ c † p, ↓ + h.c. ) + U n , ↑ n , ↓ + N X p =2 ,σ t p [ c † ,σ c p,σ + h.c. ] + N X p =2 ,σ ε p c † p,σ c p,σ , (B1)which breaks the U(1) symmetry via the term t sc,p ( c † p, ↑ c † p, ↓ + h.c. ). We have assumed site 1 to be the impurity site, and N to be the total number of sites. Using the standard procedure [41], one performs the particle-hole transformationon one of the spin species and rewrites Eq. (B1) as H And,SC = ε d ( c † , ↑ c , ↑ − c , ↓ c † , ↓ ) + N X p =1 t sc,p ( c † p, ↑ c † p, ↓ + h.c. ) + U c † , ↑ c , ↑ c † , ↓ c , ↓ + N X p =2 t p [ c † , ↑ c p, ↑ − c , ↓ c † p, ↓ + h.c. ] + N X p =2 ε p ( c † p, ↑ c p, ↑ − c p, ↓ c † p, ↓ ) + K, (B2)with the constant K = P Np =1 ε p . When computing the Green’s function, the constant K does not play any role.Using a Nambu spinor, the quadratic part of H And,SC can be expressed as (cid:16) c † , ↑ c , ↓ ... c † N, ↑ c N, ↓ (cid:17) H quad c , ↑ c † , ↓ ... c N, ↑ c † N, ↓ = (cid:16) γ † ... γ † N (cid:17) U T H quad U γ γ ... γ N . (B3)Here U is a 2 N × N unitary matrix, and one can choose it freely to simplify the calculation. Note that γ † i , whichis a linear combination of c †↑ and c ↓ , obeys the same anti-commutation relation of fermions, i.e. { γ † i , γ j } = δ ij and { γ i , γ j } = 0, and can be used to specify determinantal basis states. We refer to γ † i as a creation operator ofa Bogoliubov particle. Note that the interaction conserves the number Bogoliubov particles, and we consider theHilbert space of N Bogoliubov particles, which is the largest invariant subspace of 2 N orbitals.As an example, we take N = 4 (3 bath sites) which generates 8 (2 ×
4, 2 for spin) Bogoliubov orbitals. The space of 4Bogoliubov particles has the dimension of C = 70, and the basis state is expressed as γ † i γ † j γ † k γ † l | i . Note the vacuumstate is a state of zero Bogoliubov particles, not a state of zero “real” particles. For a given impurity Hamiltonian, oneactually performs the particle-hole transformation on spin down electrons. In terms of original particles, it impliesthe proper vacuum state fills all spin-down electrons, i.e. | i = Π N (=4) i =1 c † i, ↓ | no particle i = c † , ↓ c † , ↓ c † , ↓ c † , ↓ | no particle i . (B4)A state in subspace specified by N (= 4) Bogoliubov particles is γ † i γ † j γ † k γ † l h Π N (=4) i =1 c † i, ↓ | no particle i i . Since γ † ∼ c †↑ + c ↓ ,when expressing in the original { c i } basis, it contains subspace of 0, 2, 4, 6, 8 original particles whose representatives6are | i = | no particle i [1] | i = c † , ↑ c † , ↓ | no particle i [16] | i = c † , ↑ c † , ↑ c † , ↓ c † , ↓ | no particle i [36] | i = c † , ↑ c † , ↑ c † , ↑ c † , ↓ c † , ↓ c † , ↓ | no particle i [16] | i = c † , ↑ c † , ↑ c † , ↑ c † , ↑ c † , ↓ c † , ↓ c † , ↓ c † , ↓ | no particle ii