Static and dynamic properties of the spinless Falicov-Kimball model K.W. Beckera, S. Sykoraa, and V. Zlatica,b 7 0 aInstitut fu¨r Theoretische Physik, Technische 0 2 Universit¨at Dresden, D-01062 Dresden, Germany n a bInstitute of Physics, Bijenicka c. 46, P.O.B 304, 10000 Zagreb, Croatia J 8 (Dated: February 6, 2008) 1 ] Abstract l e - Thespinless Falicov-Kimball model is studied by the use of a recently developed projector-based r t s . renormalization method (PRM) for many-particle Hamiltonians. The method is used to evaluate t a m static and dynamic quantities of the one-dimensional model at half-filling. To these belong the - d quasiparticle excitation energy ε˜ and the momentum distribution n of the conduction electrons k k n o and spatial correlation functions of the localized electrons. One of the most remarkable results is c [ the appearance of a gap in ε˜k at the Fermi level of the order of the Coulomb repulsion U, which is 1 v accompanied by a smooth behavior for n . The density of states for the conduction electrons and k 3 3 the one-particle spectral functions for the localized electrons are also discussed. In both quantities 4 1 a gap opens with increasing U. 0 7 0 PACS numbers: 71.10.Fd, 71.27.+a,75.30.Mb / t a m - d n o c : v i X r a 1 I. INTRODUCTION The Falicov-Kimball model is a widely used model to study the properties of interacting fermions1 and has also been used to describe the valence change transition in YbInCu and 4 similar compounds2. (For a recent review see Ref. 3). The model considers a lattice of localized f-sites, which are either empty or singly-occupied site, and of conduction electron sites, which are delocalized by the nearest-neighbor hopping. There are n f-electrons and f n conduction electrons interacting by an on-site Coulomb interaction U. The two types of c electrons share the common chemical potential µ, which is always adjusted so as to keep n +n = 1. The local f-charge is a constant of motion but thermal fluctuations can change f c the average f-occupation by transferring electrons or holes from the conduction band to f-states. The relative occupation of f-states and conduction states is determined by the competition between the entropy of the f- and the conduction states, the excitation energy of the f-states, the kinetic energy of band states, and the interaction energy. The general solution of the Falicov-Kimball model is not known except in infinite dimension3,4,5. In low dimensions, the zero-temperature phase diagram is highly non-trivial and has attracted a lot of theoretical attention (see (3) for the list of references). In one dimension the thermodynamic properties of the model are more or less understood, whereas dynamical properties have been much less studied. Kennedy and Lieb6,7 proved that at low enough temperatures the half-filled Falicov-Kimball model for dimensions d 2 possesses ≥ long-range order, i.e., the ions form a checkerboard pattern, the same as in the ground state. In the ordered phase, the lattice can be divided into two inter-penetrating sublattices A and B in such a way that all nearest neighbors of a site from sublattice A belong to sublattice B and vice versa. This result holds for any value of the Coulomb interaction U. For one dimension longe range order exists only at temperature T = 0. In what follows, we first describe the Hamiltonian and the projector-based renormaliza- tion method (PRM) which was recently introduced (Sect. I). In Sect. II the method is then applied to the one-dimensional Falicov-Kimball model and renormalization equations arede- rived for the parameters of the Hamiltonian. The results for static and dynamic quantities are discussed in Sect. III and IV. To simplify the calculation, we consider only the case with translation symmetry. This excludes the discussion of physical properties which depend on lattice sites i, i.e. we do 2 not discuss the ordered phase which emerges at low temperatures. However, as discussed in Sect. V, we can still evaluate physical quantities which are given by averages over all lattice sites. Note that the PRM approach canalso be extended to situations with long rangeorder. This was shown for instance in Ref. (8), where the quantum phase transition for the spinless Holstein model was discussed as function of the electron-phonon coupling. A. Model The spinless Falicov-Kimball model is described by the Hamiltonian = t c†c +ε n +U n nc, (1) H i,j i j f i i i <i,j> i i X X X where, c† and c are the creation and annihilation operators for conduction electrons at i i site i, t are the hopping matrix elements, ε is the energy level of the localized electrons, i,j f and n = f†f and nc = c†c are the occupation number operators of the localized and i i i i i i conduction electrons at site i. The summation runs over the N sites of a periodic lattice, andU istheinteractionstrengthofthelocalCoulombrepulsionatsitesibetweenconduction and localized electrons. Note that n at each site i commutes with the Hamiltonian. Thus, i the f-electron occupation number is a good quantum number, taking only two values 1 or 0, according to whether or not site i is occupied by the localized f-electron. B. Projector-based renormalization method (PRM) Let us first discuss a recently introduced many-particle method, the projector-based renormalization method (PRM)9, which provides an approximate solution to many-particle problems defined by the Hamiltonian , H = + . (2) 0 1 H H H (0) Here, is the unperturbed part with a known eigenvalue spectrum, n = E n , 0 0 n H H | i | i and is the perturbation which does not commute with and has off-diagonal matrix 1 0 H H elements, n m = 0, but no diagonal elements. In a usual perturbative approach one 1 h |H | i 6 evaluates the corrections to the eigenstates and eigenvalues of due to . The n-th 0 1 H H order correction requires the evaluation of all matrix elements up to n ( )n m , which are 1 h | H | i difficult to calculate. 3 An alternative insight is obtained by making a unitary transformation to a new basis which has no transition matrix elements with energy differences larger than some chosen cutoff λ < Λ. Here Λ is the largest energy difference between any two eigenstates of 0 H which are connected by . This generates a new Hamiltonian 1 H = eXλ e−Xλ, (3) λ H H which can be written as a sum of two terms, = + , (4) λ 0,λ 1,λ H H H such that m n = 0 for Eλ Eλ > λ, where Eλ and n are the new eigenvalues h |H1,λ| i | n − m| n | i and eigenstates of . is chosen in such a way that it has no diagonal elements with 0,λ 1,λ H H respect to . The generator of the unitary transformation is denoted by X which has 0,λ λ H to be anti-hermitian X† = X . Note that the elimination of high energy transitions may λ − λ also generate new interaction terms which have operator structures different from that of . However, they do not connect states with an energy separation larger than λ. 1 H In the PRM method the elimination of transition matrix elements is carried out by defining a generalized projection operator P which removes from any operator those λ A parts which give rise to ’forbidden’ transitions, i.e., P = n m n m , (5) λ A | ih |h |A| i m,n X Eλ−Eλ ≤λ n m | | i.e. we retain in (5) only the states m and n with Eλ Eλ λ. The orthogonal | i | i | n − m| ≤ complement of P is Q = 1 P . Here only dyads n m contribute with E E > λ. λ λ − λ | ih | | n− m| The generator X in (3) is determined by the condition λ Q = 0, (6) λ λ H as can be seen by expanding the exponentials in (3). This gives ∞ 1 Q + X n = 0, (7) λ 1 λ H n! H ( ) n=1 X where X is a superoperator defined by X = [X , ], X 2 = [X ,[X , ]], etc. In λ λ λ λ λ λ H H H H lowest order, Q = 0 reduces to Q L X(1) = 0, where L = [ , ] defines λHλ λ H1 − 0 λ 0A H0 A n o 4 the unperturbed Liouville operator L , which commutes with Q . Obviously, the condition 0 λ Q = 0 is satisfied by the expression X(1) = L−1Q , i.e., to lowest order, X(1) can λHλ λ 0 λH1 λ be obtained from the decomposition of into the eigenmodes of L . The higher order 1 0 H correction terms to X follow systematically from the higher order commutators of (7) (for λ details see, e.g., Ref. 9). We see from Eq. (7) that all transitions in with energy transfers between the original 1 H cutoff Λ and the new cutoff λ are eliminated in one step. However it is more convenient to perform the elimination procedure step-wise, such that each step reduces the cutoff energy λ by a small amount ∆λ. Thus, the first step removes all the transitions which involve energy transfers between (the original cutoff) Λ and Λ ∆λ, the subsequent steps remove − all transitions larger than Λ 2∆λ, Λ 3∆λ, and so on. The unitary transformation for − − the step from an intermediate cutoff λ to the new cutoff λ ∆λ reads (in analogy to (3)) − = eXλ,∆λ e−Xλ,∆λ (8) λ−∆λ λ H H where = + and = + .The generator X is now fixed λ 0,λ 1,λ λ−∆λ 0,λ−∆λ 1,λ−∆λ λ,∆λ H H H H H H by the condition Q = 0. (9) λ−∆λ λ−∆λ H whichspecifiesthat containsnomatrixelementswhichconnecteigenstatesof λ−∆λ 0,λ−∆λ H H with energy differences larger than λ ∆λ. − Let us assume that the operator structure of is invariant with respect to further λ H unitary transformations. Then, we can use Eq. (7) to derive difference equations for the λ-dependence of the coupling constants of the Hamiltonian. These equations connect the parameters of the Hamiltonian with cutoff λ to those with cutoff λ ∆λ. By using a finite − number of steps we can proceed to λ 0 and obtain a set of nonlinear equations for the → renormalized parameters. The solution determines the final, fully renormalized Hamiltonian = , which depends on the initial parameter values of the original model at λ→0 0,λ→0 H H cutoff Λ. In this limit all transitions due to the interaction with nonzero transition energies have been eliminated so that identically vanishes, i.e. = 0. 1,λ→0 1,λ→0 H H The underlying idea of the PRM, namely the elimination of the interaction terms by the unitary transformations, has been used before in the literature. For instance, in Ref.10 a unitary transformation has been employed to eliminate off-diagonal matrix elements of 5 the Coulomb interaction in order to study ground state properties of the Hubbard model. Similarly, in the Schrieffer-Wolff transformation11, which maps the Anderson to the Kondo model, matrix elements connecting different charge configurations of magnetic ions are elim- inated also by use of a unitary transformation. Note however that in the PRM approach unitarytransformationsareperformedinsmallstepsincontrasttoearlier applications. Note also that the approach removes high energy transitions but does not decimate the Hilbert space. This is different to the poor man’s scaling12 which removes high energy states. It should also be mentioned that the PRM resembles the previous flow equation method13 by Wegner and the similarity renormalization approach14 by Glazek and Wilson which can be considered as continuous versions of the PRM method. II. RENORMALIZATION OF THE FALICOV-KIMBALL MODEL In order to apply the PRM we express the Falicov-Kimball interaction in terms of the fluctuations with respect to the thermal averages of nc and n , and write the diagonal and i i off-diagonal part of as H = ε +U nc n + ε +U n c†c UN n nc , (10) H0 f i i k i k k − i i i k X(cid:0) (cid:10) (cid:11)(cid:1) X(cid:0) (cid:10) (cid:11)(cid:1) (cid:10) (cid:11) (cid:10) (cid:11) and U H1 = N e−i(k−k′)Ri δni δ(c†kck′), (11) ikk′ X where δnc = nc nc , δn = n n , denotes the thermal averaging with the i i − i i i − i ··· full Hamiltonian ,(cid:10)an(cid:11)d c† = (1/√N)(cid:10) (cid:11)eik(cid:10)Ric†(cid:11)is the Fourier transform of c†. Of course, H k i i i due to the assumed translation invariaPnce hnii is equivalent to (1/N) ihnii, where N is the number of lattice sites. The sum in (11) is restricted to k = k′,Psince the diagonal 6 part of the Coulomb repulsion is included in . From translation invariance it follows 0 H δ(c†kck′) = c†kck′ − c†kck δk,k′. In principle, the average values in (10) can be defined with an arbitrary ensem(cid:10)ble s(cid:11)ince Eq. (1) and Eqs. (10), (11) are equivalent. It has turned out that the best choice is to evaluate the averages with respect to the original full Hamiltonian = , which can be done by the procedure explained below. Λ H H The perturbation causes transitions between the eigenstates of by creating 1 0 H H electron-hole pairs in the conduction band, whereas the number of f-electrons is conserved 6 at each site. Using L0 δni δ(c†kck′) = (εk −εk′) δni δ(c†kck′), H1 can be decomposed into a sum of eigenmodes of the unperturbed Liouville operator. The lowest order solution of the generator of the unitary transformation is X(1) = L−1Q , which is given by λ 0 λH1 U ei(k−k′)Ri Xλ(1) = N ikk′ εk −εk′ Θ(|εk −εk′|−λ) δni c†kck′. (12) X Here Θ( εk εk′ λ) is the projection operator which removes from 1 all transitions | − | − H with energy transfers larger than λ. Note that the same operator form also appears in the non-perturbative calculation below. A. Renormalized Hamiltonian The unitary transformation leads to the renormalized Hamiltonian = eXλ e−Xλ, λ H H which can be written as = + , where gives rise to the transitions between λ 0,λ 1,λ 1,λ H H H H the eigenstates of but has no diagonal elements. The choice of is not unique, and 0,λ 0,λ H H each particular case requires a physical intuition. But in any case it has to be such that the thermal averages which appear in the procedure can be evaluated. For the Falicov-Kimball model, the choice which emphasizes the weak coupling limit, is = ε δn + g δn δn + ε c†c +E , (13) H0,λ f,λ i ij,λ i j k,λ k k λ i i6=j k X X X U H1,λ = N Θ(λ−|εk,λ−εk′,λ|) ei(k−k′)Ri δni δ(c†kck′) , (14) ikk′ X (cid:0) (cid:1) where ε , g , and ε are the renormalized parameters which are calculated below. Note f,λ ij,λ k,λ that the interaction parameter U will not change in the renormalization procedure below so that a λ-dependence of U is not considered from the beginning. The second term in (13) is a new density-density interaction between the localized electrons which is generated during the renormalization procedure. Additional higher order operator terms, which are also generated by the renormalization, are neglected. The original model (λ = Λ) is defined by the initial condition, according to (10), (11) ε = ε +U nc ε = ε +U n (15) f,Λ f i k,Λ k i g = 0 (cid:10) (cid:11) E = N n (cid:10)(ε(cid:11)+U nc ). ij,Λ Λ − i f i (cid:10) (cid:11) (cid:10) (cid:11) 7 Once we have made the ansatz (13) and (14) for the operator structure of and , 0,λ 1,λ H H we reduce the cutoff to λ ∆λ to find the new effective Hamiltonian − = eXλ,∆λ e−Xλ,∆λ = + (16) λ−∆λ λ 0,λ−∆λ 1,λ−∆λ H H H H with renormalized parameters. The new Hamiltonian is such that does not give 1,λ−∆λ H rise to transitions with energy transfers larger than (λ ∆λ) (with respect to the new − unperturbed part ). should have the same operator structure as , however 0,λ−∆λ λ−∆λ λ H H H with renormalized parameters = ε δn + g δn δn + ε c†c +E (17) Hλ−∆λ f,(λ−∆λ) i ij,(λ−∆λ) i j k,(λ−∆λ) k k (λ−∆λ) i i6=j k X X X U + N Θ(λ−∆λ−|εk,λ−∆λ−εk′,λ−∆λ|) ei(k−k′)Ri δni δ(c†kck′) . ikk′ X (cid:0) (cid:1) B. Unitary transformation For the explicit evaluation of the unitary transformation (16) we need the generator X . We make the following ansatz λ,∆λ 1 Xλ,∆λ = N Aλk,,k∆′λ Θλkk,∆′λ e−i(k−k′)Ri δni δ(c†kck′). (18) k,k′,i X Here Θλ,∆λ is a product of Θ-functions kk′ Θλkk,∆′λ = Θ(λ−|εk,λ−εk′,λ|)×Θ |εk,(λ−∆λ) −εk′,(λ−∆λ)|−(λ−∆λ) , (19) (cid:2) (cid:3) which project all the operators in X ’on the shell’ between λ and λ ∆λ. More precisely, λ,∆λ − the energy differences εk,λ εk′,λ and εk,(λ−∆λ) εk′,(λ−∆λ) in (19) refer to the two different | − | | − | Hamiltonians and . Therefore, the two Θ-functions in (19) take into account that λ λ−∆λ H H λ possesses only transition elements with energy transfer εk,λ εk′,λ < λ, whereas in H | − | λ−∆λ no transitions with εk,(λ−∆λ) εk′,(λ−∆λ) > λ ∆λ are allowed. The coefficients H | − | − Aλ,∆λ aredetermined below, usingtheconditionQ = 0. Again,theformofX k,k′ λ−∆λHλ−∆λ λ,∆λ is suggested by its lowest order expression according to (12) (compare (9)). The unitarity of the transformation requires X = X† , so that Aλ,∆λ has to be antisymmetric with λ,∆λ − λ,∆λ k,k′ respecttotheinterchangeofk andk′. Expandingtheexponentialsin(16)thetransformation can be written as, ∞ 1 = eXλ,∆λ e−Xλ,∆λ = + Xn (20) Hλ−∆λ Hλ Hλ n! λ,∆λHλ n=1 X 8 where X = [X , ], X2 = [X ,[X , ]], etc. The commutators are λ,∆λHλ λ,∆λ Hλ λ,∆λHλ λ,∆λ λ,∆λ Hλ evaluated in Appendix A. The basic approximation of the PRM approach is to replace some of the operators generated by X by ensemble averages, so as to keep the structure of λ,∆λ λ H invariant during the renormalization procedure. This additional factorization enables us to resum the series in (20). C. Renormalization of the coupling constants By comparing the respective operators in (17) with those obtained from (20) one finds renormalization equations for the parameters of . First, from the prefactors of δn one λ i H obtains (1 2 n ) (1 cosvλ ) εf,(λ−∆λ) = εf,λ + −2N i Θλk,,k∆′λ (εk′,λ −εk,λ) C−(k kk′k)′ c†kck (21) (cid:10) (cid:11) kk′ ρ − X (cid:10) (cid:11) U(1 2 n ) sinvλ + − i Θλ,∆λ kk′ c†c . N√N k,k′ C (k k′) k k (cid:10) (cid:11) Xkk′ ρ − (cid:10) (cid:11) p Similarly, by comparing the coefficients of c†c one finds k k 1 εk,(λ−∆λ) = εk,λ+ 2 Θλk,,k∆′λ (εk′,λ −εk,λ) (1−cosvkλk′) (22) k′ X U + Θλ,∆λ C (k k′) sinvλ . √N k,k′ ρ − kk′ Xk′ q The renormalization equation for the interaction parameter g between f-electrons is best ij,λ expressed in terms of its Fourier transform g , q,λ 1 c†c c† c g = g + k k − k+q k+q Θλ,∆λ (ε ε ) (1 cosvλ ) q,(λ−∆λ) q,λ 4 C (q) k,k+q k+q,λ− k,λ − k,k+q (cid:10) (cid:11) ρ(cid:10) (cid:11) k X † † U c c c c + k k − k+q k+q Θλ,∆λ sinvλ (23) 2√N C (q) k,k+q k,k+q k (cid:10) (cid:11) (cid:10)ρ (cid:11) X 1 p ( ) − N ··· q′ X where ( ) denotes the second and third term on the r.h.s. with the wave vector q replaced ··· by q′. In this way, the exact sum rule g = 0 is fulfilled which guarantees that only q q,λ sites with i = j contribute to gij. FinallPy, also the renormalization equation for the energy 6 9 shift E can be obtained which will not be given explicitly. The other new quantities in λ Eqs. (21) to (23) are defined as follows C (k k′) 1 vλ = 2 ρ − Aλ,∆λ, C (q) = eiq(Ri−Rj) δn δn . (24) kk′ N k,k′ ρ N i j r ij X (cid:10) (cid:11) Finally we have to determine the coefficients Aλ,∆λ of the unitary transformation. They k,k′ follow from the condition Q = 0, which removes from all the operators λ−∆λ λ−∆λ λ−∆λ H H giving rise to the transitions with energy transfers larger than λ ∆λ. One finds − 1 N C (k k′) U Θλ,∆λ Aλ,∆λ = Θλ,∆λ arctan 2 ρ − . (25) kk′ k,k′ kk′ 2 sCρ(k k′) r N εk,λ εk′,λ − (cid:16) − (cid:17) Note that the Aλ,∆λ depend on the parameters at cutoff λ as well as from the reduced k,k′ cutoff λ ∆λ. The presence of Θλ,∆λ on both sides of Eq.(25) means that Aλ,∆λ has to be − kk′ k,k′ specified only for those values of k and k′ which are ’on the shell’ defined by the condition Θλ,∆λ = 1. ’Off-the-shell’ coefficients can take any value and, in what follows, we use kk′ Aλ,∆λ = 0 for Θλ,∆λ = 0. Note that expression (25) includes U to all orders. However, in k,k′ kk′ the thermodynamic limit N , the coefficients Aλ,∆λ become linear in U since all higher → ∞ k,k′ order terms vanish for N . In the actual numerical evaluation of the renormalization → ∞ equations on a lattice of finite size N the excitation energies (εk′,λ εk,λ) may become very − small for λ 0 so that the expansion of Aλ,∆λ to linear order in U breaks down. In this → k,k′ case, the full expression (25) has to be taken. By help of (25) one finds for the parameters vλ of (24) kk′ sinvλ = 2U Cρ(k −k′) sign (εk,λ−εk′,λ) kk′ r N (εk,λ εk′,λ)2 +4U2Cρ(k k′)/N − − cosvλ = |εk,λp−εk′,λ| . (26) kk′ (εk,λ εk′,λ)2 +4U2Cρ(k k′)/N − − p Equations (21) to (23) represent the final renormalization equations for the parameters of theHamiltonianastheenergycutoffisreducedfromλtoλ ∆λ. Theoverallrenormalization − starts from the original cutoff λ = Λ, where the initial parameters are given by (15). For the parameterssinvλ , cosvλ , andC (k k′), which containexpectation values, anappropriate kk′ kk′ ρ − choice is taken at cutoff Λ. Using these quantities we obtain from the renormalization 10