Quantum fluctuations beyond the Gutzwiller approximation
aa r X i v : . [ c ond - m a t . s t r- e l ] J a n Quantum fluctuations beyond the Gutzwiller approximation
Michele Fabrizio International School for Advanced Studies (SISSA), Via Bonomea 265, I-34136 Trieste, Italy (Dated: September 10, 2018)We present a simple scheme to evaluate linear response functions including quantum fluctuationcorrections on top of the Gutzwiller approximation. The method is derived for a generic multi-band lattice Hamiltonian without any assumption about the dynamics of the variational correlationparameters that define the Gutzwiller wavefunction, and which thus behave as genuine dynamicaldegrees of freedom that add on those of the variational uncorrelated Slater determinant.We apply the method to the standard half-filled single-band Hubbard model. We are able to recoverknown results, but, as by-product, we also obtain few novel ones. In particular, we show thatquantum fluctuations can reproduce almost quantitatively the behaviour of the uniform magneticsusceptibility uncovered by dynamical mean field theory, which, though enhanced by correlations,is found to be smooth across the paramagnetic Mott transition. By contrast, the simple Gutzwillerapproximation predicts that susceptibility to diverge at the transition.
PACS numbers: 71.10.-w,71.30.+h,71.10.Fd
I. INTRODUCTION
The Gutzwiller approximation is likely the simplesttool to deal with strong correlations in lattice models ofinteracting electrons. It consists in a recipe for approx-imate analytical expressions of expectation values in aclass of wavefunctions, named Gutzwiller wavefunctions,of the form | Ψ i = Y i P ( i ) | Ψ i , (1)where | Ψ i is a variational Slater determinant, and P ( i )a linear operator that acts on the local Hilbert space atsite i and depends on a set of variational parameters.Curiously, the Gutzwiller approximation often pro-vides physically more sound results than a directevaluation of expectation values in wavefunctions likeEq. (1). For instance, the numerical optimisation ona finite-dimensional lattice of a variational Gutzwillerwavefunction for a single-band half-filled Hubbard modelnever stabilises a genuine Mott insulating phase , i.e.an insulator that does not break any symmetry, whichintuitively is to be expected beyond a critical strengthof the on-site repulsion. By contrast, the Gutzwillerapproximation is instead able to describe such a genuineMott transition . The explanation of this strangeoutcome relies on the following observations. The first isthat, in order to describe a genuine Mott insulator, oneneeds to add to the Gutzwiller wavefunction, Eq. (1),long range density-density Jastrow factors . However,the effect of such Jastrow factors disappears in latticeswith coordination number z → ∞ , therefore, only inthat limit, wavefunctions like Eq. (1) can faithfullydescribe Mott insulators. Moreover, right in that limitof z → ∞ , the Gutzwiller approximation provides theexact expression of expectation values . Therefore theGutzwiller approximation should better be regarded asa recipe to evaluate approximate expectation values inGutzwiller-Jastrow wavefunctions, which becomes exact when the coordination number tends to infinity, ratherthan in Gutzwiller-only wavefunctions. In other words,the Gutzwiller approximation applied on a lattice withfinite z is just the variational counterpart of dynamicalmean field theory (DMFT) applied on that same lattice.Recently, several attempts to include the Gutzwiller ap-proximation inside DFT electronic structure codes havebeen performed with quite encouraging outcomes . Inthis perspective, it might be useful to have at disposal asimple and flexible method to calculate linear responsefunctions within the Gutzwiller approximation, inview of an extension of the so-called linear responseTDDFT to the case when DFT is combined withthe Gutzwiller approximation.There are already several works dealing with linearresponse in the Gutzwiller approximation, most of whichlimited to the single-band Hubbard model . Exten-sions to multi-band models have been attempted ,though under an assumption about the dynamics ofthe variational parameters that determine the linearoperators P ( i ) in Eq. (1).Here we shall instead present a very simple and generalmethod to evaluate linear response functions withinthe Gutzwiller approximation without any preliminaryassumption. The method is essentially an extensionof the time-dependent Gutzwiller approximation ofRef. 31 to a generic multi-band Hamiltonian, where thedynamics of the linear operators P ( i ) and of the Slaterdeterminant | Ψ i , see Eq. (1), are treated on equalfooting. Linearisation of the equations of motion aroundthe stationary solution, which is the equilibrium state,thus allows calculating linear response functions.We note that the results of the Gutzwiller approxi-mation at equilibrium coincide with the saddle pointsolution of the slave-boson theory in the path-integralformulation , which, in multi-band models, correspondsto the so-called rotationally invariant slave bosonformalism (RISB) . Our present results in the linearresponse regime can therefore be considered equivalentto the quantum fluctuations corrections above theRISB saddle-point solution. We preferred here toderive such corrections to the action directly from thetime-dependent Gutzwiller approximation rather thanfrom the RISB theory, since the former is at least a wellcontrolled variational scheme in lattices with infinitecoordination number. However, both the notations aswell as the language we shall use are actually closelyrelated to RISB theory.The paper is organised as follows. In Sec. II we brieflypresent the time-dependent Gutzwiller approximation,with some additional technical details postponed to theAppendix. In Sec. III we linearise the equations of motionaround the stationary solution and derive an effective ac-tion for the fluctuations in the harmonic approximation.In Sec. IV we apply the method to the single-band half-filled Hubbard model, which allows a comparison withalready existing results. Section V is devoted to conclud-ing remarks. II. THE GUTZWILLER APPROXIMATION INBRIEF
Besides the original works where M. Gutzwiller in-troduced a novel class of variational wavefunctions aswell as an approximate scheme to compute expectationvalues, after him called Gutzwiller wavefunctions andapproximation, and the subsequent demonstration thatsuch an approximation becomes exact in the limit ofinfinite-coordination lattices , there are by now manyarticles where the Gutzwiller approximation is describedin detail. Here we shall follow Ref. 34 and use its samenotations.The time-dependent Gutzwiller wavefunction is definedthrough | Ψ( t ) i = Y i P ( i, t ) | Ψ ( t ) i , (2)which is the analogous of Eq. (1) where now | Ψ ( t ) i is a time-dependent variational Slater determinant, and P ( i, t ) linear operators on the local Hilbert space thatdepend on time-dependent variational parameters. Forsake of simplicity, we shall not include in our analysisBCS wavefunctions nor operators P ( i, t ) that are chargenon-conserving. The extension to those cases is simple,though notations get more involved.Suppose that the Hamiltonian is written in terms offermionic operators c iα and c † iα , α = 1 , . . . , M , thatcorrespond to annihilating or creating a fermion at site i in a chosen basis of Wannier functions φ i α ( x , t ), where α indicates both spin and orbital indices. Let us imaginea U (2 M ) unitary transformation W ( i, t ) = exp (cid:18) i X αβ K αβ ( i, t ) c † iα c iβ (cid:19) , (3) with K αβ ( i, t ) = K βα ( i, t ) ∗ , which maps c i α into a newbasis set d i α of single particle operators d i α = W ( i, t ) † c i α W ( i, t ) = X β U αβ ( i, t ) c i β . (4)Evidently, if we consider the gauge transformation P ( i, t ) → P ( i, t ) W ( i, t ) † , (5) | Ψ ( t ) i → Y i W ( i, t ) | Ψ ( t ) i , (6)the Gutzwiller wavefunction | Ψ( t ) i in (2) stays invariantand the transformed | Ψ ( t ) i remains a Slater determi-nant. Such gauge invariance, analogous to that of theRISB theory , repeatedly appears in the calculationsthat follow.The most general P ( i, t ) can be written as P ( i, t ) = X n ¯ m λ n ¯ m ( i, t ) | n ; i ih ¯ m ; i | , (7)where n and ¯ m can be chosen to belong to the local basisof Fock states built with the operators c i α . Alternatively,one can use a mixed-basis representation where n labelsFock states in the original basis c i α , and ¯ m Fock states ina different basis , e.g. the basis of the operators d iα inEq. (4), which is also used to built the Slater determinant | Ψ ( t ) i . We define the uncorrelated local probability dis-tribution ˆ P ( i, t ), which is positive definite, by its matrixelements P n ¯ m ( i, t ) = (cid:10) Ψ ( t ) (cid:12)(cid:12)(cid:12) | ¯ m ; i ih ¯ n ; i | (cid:12)(cid:12)(cid:12) Ψ ( t ) (cid:11) , (8)as well as the Gutzwiller variational matrixˆΦ( i, t ) ≡ ˆ λ ( i, t ) q ˆ P ( i, t ) , (9)with matrix elements Φ n ¯ m ( i, t ). Expectation values oflocal and non-local operators in the Gutzwiller wavefunc-tion (2) can be calculated explicitly in infinite coordina-tion lattices if one imposes the following two constraintsat any time :Tr (cid:16) ˆΦ( i, t ) † ˆΦ( i, t ) (cid:17) = 1 , (10)Tr (cid:16) ˆΦ( i, t ) † ˆΦ( i, t ) ˆ c † i α ˆ c i β (cid:17) ≡ n αβ ( i, t )= h Ψ ( t ) | c † i α c i β | Ψ ( t ) i , (11)where the fermionic operators within the spur must beregarded as their matrix representation in the local Fockspace. The second constraint Eq. (11) plays the role of a gauge-fixing condition, exactly as in the RISB model .Another important ingredient is the wavefunction renor-malisation matrix ˆ R ( i, t ) with elements R αβ ( i, t ), definedby solving the set of equations h Ψ ( t ) | c † iγ P ( i, t ) † c iα P ( i, t ) | Ψ ( t ) i = X β n γβ ( i, t ) R αβ ( i, t ) , (12)where the left hand side can be straightforwardly evalu-ated by the Wick’s theorem. As shown in the AppendixA, the solution of the above equation readsˆ R ( i, t ) = ˆ Q ( i, t ) ˆ S ( i, t ) , (13)where ˆ Q ( i, t ) has matrix elements Q αβ ( i, t ) = Tr (cid:16) ˆΦ( i, t ) † ˆ c iα ˆΦ( i, t ) ˆ c † iβ (cid:17) , (14)and the hermitian matrix ˆ S ( i, t ) is defined through4 ˆ S ( i, t ) − = 1 − ˆ∆( i, t ) (15)where the matrix elements of ˆ∆( i, t ) are∆ αβ ( i, t ) = Tr (cid:18) ˆΦ( i, t ) † ˆΦ( i, t ) h ˆ c iα , ˆ c † iβ i(cid:19) . (16)The meaning of ˆ R ( i, t ) is that the action of the anni-hilation operator c iα on the Gutzwiller wavefunction isequivalent to the action of the operator P ( i, t ) † c i P ( i, t ) → ˆ R ( i, t ) c i (17)on the Slater determinant | Ψ ( t ) i , where c i is a spinorwith components c iα . One can readily show that underthe gauge transformation Eq. (5),ˆ R ( i, t ) → ˆ R ( i, t ) W = ˆ R ( i, t ) ˆ U ( i, t ) † , (18)where ˆ U ( i, t ) has the matrix elements U αβ ( i, t ) of Eq. (4),so that Eq. (17) transforms into W ( i, t ) P ( i, t ) † c i P ( i, t ) W ( i, t ) † → ˆ R ( i, t ) W d i . Since we have complete freedom in choosing W ( i, t ), aconvenient choice is the unitary transformation that diag-onalises the local single-particle density matrix, in whichcase the operators d iα are associated to the natural or-bitals and satisfyTr (cid:16) ˆΦ( i, t ) † ˆΦ( i, t ) ˆ d † i α ˆ d i β (cid:17) = δ αβ n α ( i, t ) , (19)while the matrix elements of ˆ R ( i, t ) W acquire the simpleexpression R αβ ( i, t ) W = Tr (cid:16) ˆΦ( i, t ) † ˆ c iα ˆΦ( i, t ) ˆ d † iβ (cid:17)r n β ( i, t ) (cid:16) − n β ( i, t ) (cid:17) . (20)The matrix ˆΦ( i, t ) is in this case conveniently definedin the mixed-basis representation, where n in Φ n ¯ m ( i, t )refers to a Fock state in the original basis, and ¯ m to aFock state in the natural one. Such a mixed-basis repre-sentation is useful since,throughout all calculations, onedoes not actually need to know what the natural basisis in terms of the original one . Such a nice property islinked to the gauge-invariance, equations (5) and (6), ofthe theory . A. The model
We shall assume the generic Hamiltonian H = X i = j c † i ˆ t ij c † j + X i H i , (21)where H i includes all on-site terms. If the constraintsEq. (10) and Eq. (11) are satisfied at any time t , then, ininfinite coordination lattices, it holds that E ( t ) = h Ψ( t ) | H | Ψ( t ) i = h Ψ ( t ) | H ∗ ( t ) | Ψ ( t ) i + X i Tr (cid:16) ˆΦ( i, t ) † ˆ H i ˆΦ( i, t ) (cid:17) ≡ E ∗ ( t ) + X i Tr (cid:16) ˆΦ( i, t ) † ˆ H i ˆΦ( i, t ) (cid:17) , (22)where H ∗ ( t ) = X i = j c † i ˆ R ( i, t ) † ˆ t ij ˆ R ( j, t ) c i , (23)may be interpreted as the Hamiltonian of the quasipar-ticles. Evidently, all expectation values can be straight-forwardly evaluated since the uncorrelated wavefunction | Ψ ( t ) i allows using Wick’s theorem. B. The action
In the time-domain the variational principle corre-sponds to searching for the saddle point of the action S = Z dt (cid:20) i h Ψ( t ) | ˙Ψ( t ) i − E ( t ) (cid:21) ≡ Z dt ( i X i Tr ˆΦ( i, t ) † ∂ ˆΦ( i, t ) ∂t ! + i h Ψ ( t ) | ˙Ψ ( t ) i − E ( t ) ) , (24)where the equivalence holds on provision that the con-straints (10) and (11) are fulfilled at any time. The sad-dle point equations are readily obtained: i ∂ ˆΦ( i, t ) ∂t = ˆ H i ˆΦ( i, t ) + ∂E ∗ ( t ) ∂ ˆΦ( i, t ) † , (25) i | ˙Ψ ( t ) i = H ∗ ( t ) | Ψ ( t ) i , (26)where ∂E ∗ ( t ) ∂ ˆΦ( i, t ) † = D Ψ ( t ) (cid:12)(cid:12)(cid:12)(cid:12) ∂ H ∗ ( t ) ∂ ˆΦ( i, t ) † (cid:12)(cid:12)(cid:12)(cid:12) Ψ ( t ) E ≡ ˆ T ( i, t ) ˆΦ( i, t ) . (27)ˆ T ( i, t ) is a tensor with components T nm ; n ′ m ′ ( i, t ), whichis still functional of the matrices ˆΦ and ˆΦ † at site i aswell as at all sites connected to i by the hopping. Onecan show that this tensor is hermitean, ˆ T ( i, t ) = ˆ T ( i, t ) † ,which implies that the normalisation Eq. (10) is con-served by the time evolution. C. Fate of the constraint
Concerning the second constraint, Eq. (11), we nowprove that, if it is satisfied at the initial time, it willremain so at the saddle point solutions of Eq. (25) andEq. (26). Suppose we have indeed found the saddle pointˆΦ( i, t ) and | Ψ ( t ) i . By definition, any small variationwith respect to that solution must lead to a vanishingvariation of the action. Let us consider the infinitesimalgauge transformationˆΦ( i, t ) + δ ˆΦ( i, t ) = ˆΦ( i, t ) (cid:16) − i ˆ K ( i, t ) (cid:17) , | Ψ ( t ) i + | δ Ψ ( t ) i = i X i K ( i, t ) ! | Ψ ( t ) i , where the operator K ( i, t ) = X αβ K αβ ( i, t ) c † i α c i β , (28)has infinitesimal matrix elements K αβ ( i, t ) = K βα ( i, t ) ∗ ,and ˆ K ( i, t ) is its matrix representation in the Fock space.We already mentioned that the energy E ( t ) is gauge in-variant so that the variation of the action, δ S = S W − S ,simply reads δ S = Z dt ( i X i Tr (cid:18) ˆΦ( i, t ) † ˆΦ( i, t ) ∂ ˆ W ( i, t ) † ∂t ˆ W ( i, t ) (cid:19) + i X i h Ψ ( t ) | W ( i, t ) † ˙ W ( i, t ) | Ψ ( t ) i ) ≃ X i Z dt ( Tr (cid:18) ˆΦ( i, t ) † ˆΦ( i, t ) ˙ K ( i, t ) (cid:19) − h Ψ ( t ) | ˙ K ( i, t ) | Ψ ( t ) i ) = X i X αβ Z dt ˙ K αβ ( i, t ) ( Tr (cid:18) ˆΦ( i, t ) † ˆΦ( i, t ) ˆ c † i α ˆ c † i β (cid:19) − h Ψ ( t ) | c † i α c i β | Ψ ( t ) i ) = − X i X αβ Z dt K αβ ( i, t ) ∂∂t ( Tr (cid:18) ˆΦ( i, t ) † ˆΦ( i, t ) ˆ c † i α ˆ c † i β (cid:19) − h Ψ ( t ) | c † i α c i β | Ψ ( t ) i ) . Since ˆΦ( i, t ) and | Ψ ( t ) i are solutions of the saddle pointequations, it follows that δS must strictly vanish forany choice of the infinitesimally small matrix elements K αβ ( t ), which implies ∂∂t ( Tr (cid:18) ˆΦ( i, t ) † ˆΦ( i, t ) ˆ c † i α ˆ c † i β (cid:19)) − h Ψ ( t ) | c † i α c i β | Ψ ( t ) i ) = 0 , thus just the desired result. It actually means thatthe term in parenthesis is conserved in the evolution.Therefore, if it is initially vanishing, it will remain so atany time, which thus implies that the constraint Eq. (11)is fulfilled during the whole time evolution. D. Stationary problem
At equilibrium one needs to find the minimum of theenergy with the two constraints Eqs. (10) and (11), whichcan be enforced e.g. by Lagrange multipliers, leading tothe set of equationsΛ( i ) ˆΦ( i ) = (cid:16) ˆ H i + ˆ T ( i ) (cid:17) ˆΦ( i ) + X αβ µ αβ ( i ) ˆΦ( i ) ˆ d † iα ˆ d iβ , (29) E ∗ | Ψ i = (cid:16) H ∗ − X i µ αβ ( i ) d † i α d i β (cid:17) | Ψ i , (30)where Λ( i ) enforces Eq. (10), and the hermitean matrixˆ µ ( i ) with components µ αβ ( i ) enforces Eq. (11). In what-ever follows we shall assume to work in a mixed-basisrepresentation where the operators d iα are associated tothe natural orbitals, so that we must also ensure thatTr (cid:16) ˆΦ † ( i ) ˆΦ( i ) ˆ d † iα ˆ d iβ (cid:17) = h Ψ | d † iα d iβ | Ψ i = δ αβ n α ( i ) . The quasiparticle Hamiltonian in the natural basis, in-cluding explicitly the Lagrange multipliers, is therefore H ∗ → X i = j d † i ˆ R ( i ) † ˆ t ij ˆ R ( j ) d i − X i d † i ˆ µ ( i ) d i , (31)with ˆ R defined in Eq. (20). Working in the mixed-basisrepresentation with the natural orbitals considerably sim-plifies all calculations.Recalling that ˆ T ( i ) is still functional of ˆΦ, Eq. (29) lookslike a stationary non-linear Schrœdinger equation .One can for instance solve it as in any Hartree-Fock cal-culation. Namely, one can find the eigenstates and eigen-values of Eq. (29) assuming ˆ T ( i ) fixed, and impose that,when ˆ T ( i ) is calculated substituting the actual expres-sion of the lowest energy solution ˆΦ ( i ), the two valuescoincide. The Lagrange multiplier ˆ µ is fixed by imposingEq. (11) and Eq. (19). In this way one finally gets theself-consistent ˆ T ( i ), which we shall hereafter denote asˆ T (0) ( i ) ≡ ˆ T h ˆΦ , ˆΦ † i . (32)Once the latter is known, as well as the value of ˆ µ , onecan also solve (29) for all eigenvectors, ˆΦ n ( i ) and corre-sponding eigenvalues E n ( i ), with E ( i ) = Λ( i ). We shalldenote H ∗ , ˆ R , ˆ Q , ˆ n and ˆ S calculated with ˆΦ as H (0) ∗ ,ˆ R (0) , ˆ Q (0) , ˆ n (0) and ˆ S (0) , respectively, with the latter twomatrices diagonal in the natural basis, n (0) αβ = δ αβ n (0) α , (33) S (0) αβ = δ αβ S (0) α = δ αβ (cid:16) n (0) α (cid:16) − n (0) α (cid:17)(cid:17) − / . (34)We conclude by noting that the saddle point Hamil-tonian Eq. (31) with the inclusion of the Lagrange mul-tipliers is not anymore invariant under the most general U (2 M ) gauge transformation, but only under a subgroup G with generators ˆ T a that commute with ˆ µ . This is com-mon in theories where the gauge invariance implementsconstraints about physical states. In the natural basisrepresentation, µ i,αβ = δ αβ µ iα is diagonal, so that thematrix elements of ˆ T a must satisfy T ai,αβ (cid:0) µ iα − µ iβ (cid:1) = 0 , (35)whose solution is straightforward. For any non-degenerate α , i.e. such that µ iα = µ iβ , ∀ β = α , weassociate the generators T αi,γβ = δ αβ δ γβ of U (1) abeliangroups. On the contrary, for any set of α i , i = 1 , . . . , k ,such that µ iα i = µ iα j = µ iβ , ∀ β = α , . . . , α k , we canassociate generators of a U ( k ) Lie algebra. III. FLUCTUATIONS ABOVE THE SADDLEPOINT SOLUTION
Our goal is to determine the action of the fluctuationsbeyond the saddle point within the harmonic approxima-tion. To that purpose we assume thatˆΦ( i, t ) = e − iE t X n φ n ( i, t ) ˆΦ n ( i ) ˆ W ( i, t ) † , (36)where φ n ( i, t ) for n > φ ( i, t ) = 1 − X n> | φ n ( i, t ) | . (37)In addition, the Slater determinant is defined through | Ψ ( t ) i → e − iE ∗ t W ( t ) | Ψ ( t ) i , (38) where | Ψ ( t ) i is properly normalised and includes thezeroth order | Ψ (0)0 i , solution of the saddle point, as well asa fluctuation correction | δ Ψ ( t ) i . The unitary operator W ( i, t ) = exp (cid:18) − i t d † i ˆ µ ( i ) d i (cid:19) , (39)where ˆ µ ( i ) is the equilibrium Lagrange multiplier, andˆ W ( i, t ) is the matrix representation of W ( i, t ).Through the above definitions, the action becomes S = Z dt ( i X i X n> φ n ( i, t ) ∗ ˙ φ n ( i, t ) + i h Ψ ( t ) | ˙Ψ ( t ) i− X i X nm φ n ( i, t ) ∗ V nm ( i ) φ m ( i, t )+ E + E ∗ − E ∗ ( t ) ) , (40)where E ∗ ( t ) = h Ψ ( t ) | H ∗ ( t ) | Ψ ( t ) i , being now H ∗ ( t ) = X i = j d † i ˆ R ( i, t ) † ˆ t ij ˆ R ( j, t ) d i − X i d † i ˆ µ ( i ) d i , (41)and V nm ( i ) = Tr (cid:16) ˆΦ n ( i ) † ˆ H i ˆΦ m ( i ) (cid:17) (42)+Tr (cid:16) ˆΦ n ( i ) † ˆΦ m ( i ) ˆ d † i ˆ µ ( i ) ˆ d i (cid:17) . We expand H ∗ ( t ) up to second order in the fluctu-ations. The zeroth order is just H (0) ∗ . Since the sta-tionary solution is the saddle point of the action, theexpectation value of the first order expansion H (1) ( t )over the saddle point Slater determinant | Ψ (0)0 i can-cels with the first order expansion of the local energy P i P nm φ n ( i, t ) ∗ V nm ( i ) φ m ( i, t ). Therefore H (1) ( t ) con-tributes to E ∗ ( t ) with a second order term that, by linearresponse theory, reads δ E ∗ ( t ) = h δ Ψ ( t ) | H (1) ( t ) | Ψ (0)0 i + c.c. = − i Z t dτ (cid:10) h H (1) ( t ) , H (1) ( τ ) i (cid:11) , (43)where, hereafter, h . . . i will denote average over | Ψ (0)0 i ,and the operators in Eq. (43) have an additional time de-pendence since are evolved with the saddle point Hamil-tonian H (0) ∗ . The explicit expression of H (1) ( t ) is H (1) ( t ) = X i = j " d † j ˆ R (0) ( j ) † ˆ t ji ˆ R (1) ( i, t ) d i + H.c. , (44)where ˆ R (0) ( i ) is the stationary value, while the explicitexpression of the first order Taylor expansion ˆ R (1) ( i, t )is given in Appendix A 1, see Eq. (A17).There are several second order terms upon expanding H ∗ ( t ), which we shall consider separately. The first issimply H (2)1 ( t ) = X i = j d † i ˆ R (1) ( i, t ) † ˆ t ij ˆ R (1) ( j, t ) d j , (45)whose expectation value over | Ψ (0)0 i is an additional sec-ond order contribution δ E ∗ ( t ) = (cid:10) H (2)1 ( t ) (cid:11) , (46)which, together with δ E ∗ ( t ) in Eq. (43), endow the ac-tion with spatial correlations among the φ n ( i, t )’s at dif-ferent sites.The next second order corrections to H ∗ ( t ) derive fromthe second order expansion of ˆ R ( i, t )ˆ R (2) ( i, t ) = ˆ R (2)1 ( i, t ) + ˆ R (2)2 ( i, t ) , (47)where we distinguish two different contributions, seeequations (A19) and (A20) in Appendix A 1. The rea-son of this distinction is that δ E ∗ ( t ) = X i X nm φ n ( i, t ) ∗ V nm ( i ) φ m ( i, t )+ h X i = j (cid:16) d † j ˆ R (0) ( j ) ˆ t ji ˆ R (2)1 ( i, t ) d i + H.c. (cid:17) i = X n> (cid:16) E n − E (cid:17) φ n ( i, t ) ∗ φ n ( i, t ) , (48)reproduces the bare excitation energy of the fluctuations.The last contribution to the energy of the fluctuations istherefore δ E ∗ ( t ) = h X i = j (cid:16) d † j ˆ R (0) ( j )ˆ t ji ˆ R (2)2 ( i, t ) d i + H.c. (cid:17) i . (49)If we define new variables x n ( i, t ) = 1 √ (cid:16) φ n ( i, t ) + φ n ( i, t ) ∗ (cid:17) , (50) p n ( i, t ) = − i √ (cid:16) φ n ( i, t ) − φ n ( i, t ) ∗ (cid:17) , (51)and the quadratic potential U (cid:0) t, { x, p } (cid:1) = δ E ∗ ( t ) + δ E ∗ ( t ) + δ E ∗ ( t ) , (52)which has a retarded component δ E ∗ ( t ), see Eq. (43),the action of the fluctuations reads, upon defining ω n = E n − E , δ S = Z dt ( X i X n> (cid:20) p n ( i, t ) ˙ x n ( i, t ) (53) − ω n (cid:16) x n ( i, t ) + p n ( i, t ) (cid:17)(cid:21) − U (cid:0) t, { x, p } (cid:1)) , which is just the action of coupled harmonic oscillators. δS in Eq. (53) can be for instance used to evaluate the fluctuation corrections to linear response functions of lo-cal operators. For any local observable ˆ O ( i ), let us definethe matrix element O n ( i ) ≡ Tr (cid:16) ˆΦ n ( i ) † ˆ O ( i ) ˆΦ ( i ) (cid:17) . (54)Suppose we add a perturbation that couples to the localdensity matrix δ H ( t ) = X i c † i ˆ V ( i, t ) c i , (55)where the matrix ˆ V ( i, t ) with elements V αβ ( i, t ) repre-sents the external field. Without loss of generality we canassume that the expectation value of δ H ( t ) in Eq. (55)vanishes at the stationary solution. Since by assumptionthe external field is first order, the perturbation adds asecond order correction to the action (53) that is V ( t ) = X i X n " φ n ( i, t ) ∗ Tr (cid:16) ˆΦ n ( i ) † ˆ c † i ˆ V ( i, t ) ˆ c i ˆΦ ( i ) (cid:17) + φ n ( i, t ) Tr (cid:16) ˆΦ ( i ) † ˆ c † i ˆ V ( i, t ) ˆ c i ˆΦ n ( i ) (cid:17) (56) ≡ X i X n (cid:16) φ n ( i, t ) ∗ V n ( i, t ) + φ n ( i, t ) V n ( i, t ) ∗ (cid:17) = √ X i n (cid:18) ℜ eV n ( i, t ) x n ( i, t ) + ℑ mV n ( i, t ) p n ( i, t ) (cid:19) . In the presence of V ( t ) the action transforms into that offorced harmonic oscillators, whose solution allows calcu-lating the expectation value of any local operator ˆ O ( i ),see Eq. (54), O ( i, t ) = Tr (cid:16) ˆΦ( i, t ) † ˆ O ( i ) ˆΦ( i, t ) (cid:17) ≃ √ X n (cid:18) ℜ eO n ( i ) x n ( i, t ) + ℑ mO n ( i ) p n ( i, t ) (cid:19) , at linear order in the external field. A. Residual gauge invariance and would-beGoldstone modes
As we mentioned, the action Eq. (40), with the time de-pendent quasiparticle Hamiltonian defined in Eq. (41), isinvariant under a subgroup G of the initial U (2 M ) gaugesymmetry. This implies the existence of massless modeswith singular propagators that diverge as 1 /ω at low fre-quency, which are the would-be Goldstone modes relatedto the fact that the saddle-point ˆΦ ( i ) is not invariantunder G . Let us consider for instance a U (1) subgroupof G related to the non-degenerate state α in the naturalbasis. The associated adjoint charge is n α ( i, t ) ≃ X n> (cid:16) φ n ( i, t ) ∗ Tr (cid:16) ˆΦ n ( i ) † ˆΦ ( i ) ˆ d † iα ˆ d iα (cid:17) + c.c. (cid:17) , and its conjugate variable is readily found to be ϕ α ( i, t ) ≃ i n (0) α X n> (cid:18) φ n ( i, t ) ∗ Tr (cid:16) ˆΦ n ( i ) † ˆΦ ( i ) ˆ d † iα ˆ d iα (cid:17) − φ n ( i, t ) Tr (cid:16) ˆΦ ( i ) † ˆΦ n ( i ) ˆ d † iα ˆ d iα (cid:17)(cid:19) . The role of ϕ α ( i, t ) is just to enforce the constraintEq. (11), i.e. n α ( i, t ) = h Ψ ( t ) | c † iα c iα | Ψ ( t ) i ≡ h c † iα c iα i t . Indeed we can always perform a gauge transformation onthe fermions c iα → e − iϕ α ( i,t ) c iα , which makes ϕ α ( i, t ) to disappear from the energy leavingjust the time derivative term in the action, δ S = − Z dt ˙ ϕ α ( i, t ) (cid:16) n α ( i, t ) − h c † iα c iα i t (cid:17) . The condition of vanishing derivative with respect to ϕ α ( i, t ) is therefore just the condition that the constraintis conserved.It follows that we can always drop from the action allterms that contain the variables conjugate to the adjointcharges associated with the gauge symmetry G , on pro-vision that, wherever n α ( i, t ) appears, we replace it with h c † iα c iα i t .However, the above procedure does not involve all the co-efficients φ n ( i, t ); some of their linear combinations areuntouched by gauge-fixing and remain genuine indepen-dent dynamical degrees of freedom . This fact, ratherthan being a limitation, it endows the theory with a richerdynamics. IV. APPLICATION TO THE HALF-FILLEDHUBBARD MODEL
We now apply the above formalism to the simple caseof a single band Hubbard model at half-filling, whereall calculations can be worked out analytically andwhich also allows for a direct comparison with previousworks . We will show that we can indeedrecover known results, but also find few novel ones.The Hamiltonian is in this case H = − t √ z X
00 ˆΦ s (cid:19) , (58)where the charge component, i.e. the matrix elements inthe subspace (cid:0) | i , | i (cid:1) , isˆΦ c = (cid:18) φ c + φ c φ c − φ c (cid:19) = φ c σ + φ c σ , (59)with σ the 2 × σ i , i = 1 , . . . , (cid:0) |↑i , |↓i (cid:1) , is insteadˆΦ s = X i =0 φ si σ i = φ s σ + φ s · σ , (60)which allows a full spin- SU (2) invariant analysis .Normalisation implies that1 = (cid:12)(cid:12) φ c (cid:12)(cid:12) + (cid:12)(cid:12) φ c (cid:12)(cid:12) + (cid:12)(cid:12) φ s (cid:12)(cid:12) + φ ∗ s · φ s . One can readily verify that the matrix ˆ Q with compo-nents Q σσ ′ = Tr (cid:16) ˆΦ † c σ ˆΦ c † σ ′ (cid:17) , (61)can be written as ˆ Q = Q σ + Q · σ , (62)where2 Q = (cid:16) φ ∗ c φ s + φ ∗ s φ c (cid:17) + (cid:16) φ ∗ c φ s − φ ∗ s φ c (cid:17) , (63)2 Q i = (cid:16) φ ∗ c φ si − φ ∗ si φ c (cid:17) + (cid:16) φ ∗ c φ si + φ ∗ si φ c (cid:17) , (64)with i = 1 , . . . ,
3. Seemingly,ˆ∆ ≡ (cid:16) φ ∗ c φ c + φ ∗ c φ c (cid:17) σ − (cid:16) φ ∗ s φ s + φ s φ ∗ s + i φ ∗ s ∧ φ s (cid:17) · σ ≡ ∆ σ + ∆ · σ . (65) A. Stationary solution
As common when discussing the Mott transition in thesingle band Hubbard model, we shall be interested in thestationary solution within the paramagnetic sector, i.e.neglecting spontaneous breakdown of spin SU (2) sym-metry. Such solution at half-filling is characterised by asite independentˆΦ = 1 √ φ (0) c σ φ (0) s σ ! , with 1 = (cid:12)(cid:12) φ (0) c (cid:12)(cid:12) + (cid:12)(cid:12) φ (0) s (cid:12)(cid:12) . Under this assumptionˆ R ( i ) = (cid:16) φ (0) c ∗ φ (0) s + φ (0) s ∗ φ (0) c (cid:17) σ = R (0) σ , ∀ i , (66)so that the quasiparticle Hamiltonian is just a tight-binding model with renormalised hopping, i.e. H (0) ∗ = − t √ z R (0)2 X
00 cos θ σ ! , (71)where cos θ = min (1 , U/U c ), and has eigenvalue E = − U θ = − Max (
U, U c )4 . (72)We can now find all other eigenvalues and eigenvectors.The highest energy one isˆΦ = 1 √ cos θ σ − sin θ σ ! , (73)with eigenvalue E = − E . (74)This eigenstate actually corresponds to the high energyHubbard bands.The lowest excited eigenstate is threefold degenerate ( i =1 , ,
3) ˆΦ i = 1 √ (cid:18) σ i (cid:19) , (75)with eigenvalue E = − U , (76)and describes spin fluctuations. We note that above theBrinkmann-Rice transition, U > U c , this magnetic statebecomes degenerate with the ground state. In what fol-lows we shall anyway expand always around ˆΦ , and, toavoid problems, we will mostly consider the metal phaseat U ≤ U c .Finally, the last eigenstate isˆΦ = 1 √ (cid:18) σ
00 0 (cid:19) , (77)with eigenvalue E = + U , (78)and describes instead charge fluctuations. This modebecomes degenerate with ˆΦ above the transition. B. Action of the fluctuations
Following section III we writeˆΦ( i, t ) = φ ( i, t ) ˆΦ + X i =1 φ i ( i, t ) ˆΦ i + X n =2 φ n ( i, t ) ˆΦ n , (79)with φ ( i, t ) fixed by normalisation. Through equations(62), (63) and (64) we find thatˆ R (1) ( i, t ) = sin θ (cid:16) φ ( i, t ) − φ ( i, t ) ∗ (cid:17) · σ − cos θ (cid:16) φ ( i, t ) − φ ( i, t ) ∗ (cid:17) + cos θ (cid:16) φ ( i, t ) + φ ( i, t ) ∗ (cid:17) ≡ i √ θ p ( i, t ) · σ − i √ θ p ( i, t )+ √ θ x ( i, t ) , (80)where we have introduced the conjugate variables asso-ciated with φ n and φ ∗ n . Eq. (44) reads explicitly H (1) ∗ = X i ( √ (cid:18) θ (cid:19) − ∇ · J s ( i ) · p ( i, t ) −√ (cid:18) θ (cid:19) − ∇ · J c ( i ) p ( i, t )+2 √ θ h ∗ ( i ) x ( i, t ) ) , (81)where ∇ is the lattice divergence, J s ( i ) and J c ( i ) thespin and charge currents, respectively, defined throughthe continuity equations i ∂∂t (cid:16) c † i σ c i (cid:17) = h c † i σ c i , H (0) ∗ i ≡ − i ∇ · J c ( i ) , (82) i ∂∂t (cid:16) c † i σ c i (cid:17) = h c † i σ c i , H (0) ∗ i ≡ − i ∇ · J s ( i ) . (83)and finally h ∗ ( i ) the Hamiltonian density h ∗ ( i ) = − t √ z R (0)2 X j n.n. i (cid:16) c † i c j + H.c. (cid:17) . (84)Therefore δ E ∗ ( t ) defined in Eq. (43) becomes, due toparticle-hole and spin SU (2) symmetry δ E ∗ ( t ) = X i,j Z dτ (
11 + cos θ χ ∇ J ∇ J ( i − j, t − τ ) p ( i, t ) · p ( j, τ ) + 11 − cos θ χ ∇ J ∇ J ( i − j, t − τ ) p ( i, t ) p ( j, τ )+ 8 cot θ χ h ∗ h ∗ ( i − j, t − τ ) x ( i, t ) x ( j, τ ) ) , (85)where χ ∇ J ∇ J is the linear response function of ∇ J withthe Hamiltonian H (0) ∗ , which is actually the same forcharge and spin currents, and χ h ∗ h ∗ the response functionof h ∗ . We observe that, because of charge and spin con-tinuity equations, in Fourier space the following equiva-lence holds2 T sin θ (cid:16) γ − γ q (cid:17) + χ ∇ J ∇ J ( q , ω ) = ω χ ( q , ω ) , (86)where χ ( q , ω ) is the density-density response function,which is the same both in the charge and spin channels, and by definition γ q = 2 z d X i =1 cos q i ∈ [ − , +1] . (87)Without going into further details, we find that thefollowing expressions for the remaining contributions δ E ∗ ( t ) in Eq. (46), and δ E ∗ ( t ) in Eq. (49): δ E ∗ ( t ) = − T z X
As we mentioned, the Hubbard bands may be associ-ated with the excited state ˆΦ , hence with the operators0 x and p . Their equations of motion in Fourier spaceare − iω x ( q , ω ) = ω p ( q , ω ) , (90) − iω p ( q , ω ) = − h ω − T γ q + 8 cot θ χ h ∗ h ∗ ( q , ω ) i x ( q , ω ) . (91)Within the metal phase, U < U c , ω = E − E = 4 T , sothat, upon defining cos θ = U/U c ≡ u , and noting that,for small | q | , χ hh ( q , ω ) = O ( q ), the eigenmode energyis solution of the equation ω q = 4 T (cid:20) T (cid:0) − u (cid:1) + 4 T u (cid:0) γ − γ q (cid:1) +8 u − u χ h ∗ h ∗ ( q , ω q ) (cid:21) (92) ≃ T "(cid:0) − u (cid:1) + u (cid:0) γ − γ q (cid:1) , thus describes an optical mode that softens at the metalinsulator transition, ω = 4 T √ − u → u →
1. We observe that the continuum of quasiparticle-quasihole excitations extends up to an energy of order T (cid:0) − u (cid:1) , so that, upon approaching the transition, ω q must detach from the continuum and become agenuine coherent excitation.This coherent mode actually corresponds to the spin-wave excitations of the Ising field within the Z slave-spinrepresentation of the Hubbard model . This is notsurprising since, as shown in Ref. 27, the Gutzwillerwavefunction is just the mean-field variational state ofthe Z slave-spin theory. At the mean-field level, theMott transition in this representation translates intothe order-disorder transition of a quantum Ising model.Therefore the mode x seems to be the real fingerprintof the Mott transition. D. Dynamical charge susceptibility
We assume to perturb the system in the metal phase, u ≤
1, by an external potential that couples to the chargedeviation from half-filling, namely δ H ( t ) = X i v ( i, t ) (cid:0) n i − (cid:1) ≃ −√ θ X i v ( i, t ) x ( i, t ) . (93)Since ω = E − E = 2 T (cid:0) u (cid:1) and by means ofEq. (86), we find in the presence of the field the followingequations of motion for the conjugate variables x and p − iω x ( q , ω ) = ω − u χ ( q , ω ) p ( q , ω ) , − iω p ( q , ω ) = √ θ v ( q , ω ) − T (cid:0) u (cid:1) u (cid:0) − u (cid:1) x ( q , ω ) , from which it follows that the dynamical charge suscep-tibility is χ c ( q , ω ) = (1 − u ) χ ( q , ω )(1 − u ) − T (cid:0) u (cid:1) u (cid:0) − u (cid:1) χ ( q , ω ) ≡ χ ( q , ω )1 + Γ c χ ( q , ω ) , (94)where it is evident the analogy with conventional RPA,though with a renormalised coupling constantΓ c = − U u − u (cid:16) − u (cid:17) < . (95)We note that χ ( q → , ω = 0) = −N ∗ , where N ∗ = N − u , (96)is the quasiparticle density of states (DOS) at the chem-ical potential, as opposed to the bare DOS N , anddiverges approaching the Mott transition. Therefore,through Eq. (94), the charge compressibility is readilyobtained κ = N ∗ − Γ c N ∗ ≡ N ∗ F S , and defines the Landau F S parameter F S = −N ∗ Γ c . (97)Since approaching the transition, u → F S ∼ (1 − u ) − diverges faster than N ∗ ∼ (1 − u ) − , we find that thecharge compressibility correctly vanishes at the MIT. Theexpression of F S coincides with that originally obtainedby Vollhardt .In the opposite limit of small | q | with respect to fre-quency, χ ( q , ω ) ≃ T (cid:0) − u (cid:1)(cid:0) γ − γ q (cid:1) ω , which, inserted into Eq. (94), allows calculating the polesof the dynamical charge susceptibility, which are ω c q = 4 T (1 + u ) u (2 − u ) (cid:0) γ − γ q (cid:1) . (98)This acoustic mode is above the quasiparticle-quasiholecontinuum and actually corresponds to the Landau’s zerosound. Once again this result is compatible with Voll-hardt’s description of the correlated metal within theGutzwiller approximation in the framework of Landau-Fermi liquid theory . Indeed the zero sound velocity hasthe expected Landau’s expression, once one realises that1in a lattice with infinite coordination F S = 0 and it isunrelated to the enhancement of the effective mass.We conclude highlighting that the velocity of the zerosound stays constant approaching the Mott transition.In particular, for ω ≫ T (cid:0) − u (cid:1) (cid:0) γ − γ q (cid:1) , the dynam-ical charge susceptibility can be written as χ c ( q → , ω ) = 2 T (1 − u ) (cid:0) γ − γ q (cid:1) ω − ω c q , (99)hence the pole at the zero sound has vanishing weight asthe transition u → ( q , ω )of p ( q , ω )Π ( q , ω ) = − ω (1 − u )Γ c c χ ( q , ω ) , is singular at ω = 0, although this singularity doesnot appear in the physical response function, which isproportional to the propagator of the conjugate vari-able x ( q , ω ). Indeed, p ( q , ω ) is one of the would-beGoldstone modes that we mentioned in section III A.The action of the single-band Hubbard model is U (2) = U (1) × SU (2) gauge invariant, and p ( q , ω ) is just thewould-be Goldstone mode associated with the abelian U (1), whereas we shall see that p ( q , ω ) are instead thoseassociated with SU (2). In fact, the RPA form of thecharge susceptibility could be very easily obtained by thegauge-fixing prescription of section III A. If we drop allterms that contain p ( i, t ) and replace −√ θ x ( i, t ) → h n i − i t , we get an effective Hamiltonian of the quasiparticles, ne-glecting for convenience all other variables but x ( i, t ), H ∗ ( t ) = H (0) ∗ + X i v ∗ ( i, t ) (cid:0) n i − (cid:1) , where v ∗ ( i, t ) = v ( i, t ) − Γ c h n i − i t , (100)which readily leads to Eq. (94). E. Dynamical spin susceptibility
In order to study the spin response, we imagine to addan external field that couples to the spin density, e.g. toits z component, namely δ H ( t ) = − X i B ( i, t ) (cid:0) n i ↑ − n i ↓ (cid:1) = −√ θ X i B ( i, t ) x , ( i, t ) . (101) In the metal phase ω = E − E = 2 T (cid:0) − u (cid:1) , and re-peating all calculations done for the charge susceptibility,we finally obtain the dynamical spin susceptibility χ s ( q , ω ) = χ ( q , ω )1 + Γ s χ ( q , ω ) , (102)where Γ s = U − u u (cid:16) u (cid:17) > . (103)The above expression reproduces the small u Stoner’senhancement of the magnetic susceptibility. In addi-tion it satisfies the relationship Γ s ( U ) = Γ c ( − U ) validat particle-hole symmetry . Since Γ s ∼ (1 − u ) van-ishes linearly approaching the transition, the Landau’sparameter F A = −N ∗ Γ s < , (104)is constant for u →
1, which implies that the uniformstatic spin susceptibility diverges at the MIT. This resultagrees with previous ones also obtained within theGutzwiller approximation, but contrasts DMFT, whichinstead finds a finite uniform spin susceptibility at thetransition.Such negative outcome critically depends from the factthat the effective interaction Γ s , Eq. (103), vanishes atthe transition. We are going to show that beyond theharmonic approximation this cancellation does not occuranymore.We note that p a ( i, t ), a = 1 , . . . ,
3, are now the Gold-stone modes associated with SU (2) gauge invariance, andtheir propagatorsΠ a ( q , ω ) = − ω (1 + u ) Γ s s χ ( q , ω ) , diverge at ω = 0. We can, as in section IV D, drop p ( i, t )from the action and replace √ θ x ( i, t ) → h c † i σ c i i t , whose effect could be absorbed into an effective magneticfield B ∗ a ( i, t ) = δ a B ( i, t ) − Γ s h c † i σ a c i i t , (105)that straightforwardly leads to Eq. (102). F. Beyond RPA in the x mode We observe that all the above results in the metal phasecorrespond to expanding the action at second order inthe fluctuations but treating the linear coupling betweenthe latter and the fermions just within RPA, i.e. not ac-counting for exchange processes. While this procedure issomehow forced by gauge invariance for what it concerns2charge and spin modes, see the ending parts of sectionsIV D and IV E, it is not really compulsory for the x ( i, t )mode that describes the Hubbard bands. We can there-fore take a first step forward when dealing with x ( i, t ) inthe direction of the so called RPA+Exchange. Accordingto Eq. (81), promoting x and p to quantum conjugatevariables, after defining t ∗ = t sin θ and X ( i ) = 1 + √ θ x ( i ) , the Hamiltonian reads H ∗ = − t ∗ √ z X
In this paper we have presented a quite simplemethod to calculate linear response functions within theGutzwiller approximation, including in a consistent wayquantum fluctuations in the harmonic approximation.The calculation is straightforward and just requires alittle more effort than the equilibrium one. In fact,besides the variational matrix ˆΦ that minimises theenergy at equilibrium, and which can be regarded asthe lowest energy eigenstate of a local Hamiltonian ,see Eq. (29), one also needs all excited eigenstates andeigenvalues. In a model that involves M correlatedorbitals in each unit cell, this local Hamiltonian isdefined in a Hilbert space of dimension (cid:0) M M (cid:1) , and canbe conveniently recast into the problem of an impuritywith M orbitals hybridised to a single bath site with thesame number of orbitals, the coupled system being athalf-filling .As a check we have applied the method to the single-band Hubbard model at half-filling and recovered allknown results . As a by-product, we alsoshowed how to cure one flaw of the Gutzwiller approxima-tion, i.e. the divergence of the uniform magnetic suscep-tibility approaching the Mott transition from the metalside.3 ACKNOWLEDGMENTS
This work has been supported by the European Unionunder H2020 Framework Programs, ERC AdvancedGrant No. 692670 “FIRSTORM”.
Appendix A: The wavefunction renormalisationmatrix ˆ R ( i ) At equilibrium and in the natural basis, the constraintEq. (11) readsTr (cid:16) ˆΦ ( i ) † ˆΦ ( i ) ˆ d † iα ˆ d iβ (cid:17) = Tr (cid:16) ˆ P (0)0 ( i ) ˆ d † iα ˆ d iβ (cid:17) = δ αβ n (0) α ( i ) , where ˆ P ( i ) is the local probability distribution of theSlater determinant. Hereafter we shall drop for simplicitythe site index i .We can always write ˆ P (0)0 as the Boltzmann distributionof a non-interacting Hamiltonian H = X α ǫ α n α , where f (cid:0) ǫ α (cid:1) = n (0) α is the Fermi distribution function. IfˆΦ is varied, also the probability distribution must varyin such a way as to preserve the constraint. This changewill generally correspond to H → H + δH . Since H must still be a one body Hamiltonian it followsthat d α ( τ ) = e τH d α e − τH = (cid:16) e − ˆ H τ d (cid:17) α = X β U βα ( τ ) d β , where ˆ H is the matrix representation of H in the single-particle basis, so that d α ( τ ) remains a combination ofcreation operators. Since ˆ U ( τ ) ˆ U ( τ ) = ˆ U ( τ + τ ), ittrivially holds that ˆ U ( τ ) ˆ U ( − τ ) = 1 andˆ U ( β/
2) ˆ U ( β/
2) = ˆ U ( β ) . (A1)The local probability distributionˆ P = e − β ˆ H Tr (cid:16) e − β ˆ H (cid:17) , so thatTr (cid:16) ˆ P ˆ d β ( β ) ˆ d † α (cid:17) = Tr (cid:16) ˆ P ˆ d † α ˆ d β (cid:17) ≡ n αβ = X γ U γβ ( β ) Tr (cid:16) ˆ P d γ d † α (cid:17) = X γ U γβ ( β ) (cid:16) δ αγ − n αγ (cid:17) = U αβ ( β ) − X γ n αγ U γβ ( β ) , namelyˆ U ( β ) = (cid:16) − ˆ n (cid:17) − ˆ n = − (cid:16) − ˆ n (cid:17) − , (A2)which relates ˆ U ( β ) to ˆ n . It also follows thatˆ U ( − β ) = (cid:16) − ˆ n (cid:17) ˆ n − = ˆ n − (cid:16) − ˆ n (cid:17) = ˆ n − − . (A3)The renormalisation coefficients R is obtained by solvingfor any α and γ Tr q ˆ P ˆΦ † ˆ c † α ˆΦ 1 q ˆ P ˆ d γ ! = X β Tr (cid:16) ˆΦ † ˆΦ ˆ d † β ˆ d γ (cid:17) R ∗ αβ , (A4)where 1 q ˆ P ˆ d γ q ˆ P = e β ˆ H/ ˆ d γ e − β ˆ H/ = ˆ d γ ( β/
2) = X β U βγ ( β/
2) ˆ d δ . Therefore, once we define Q ∗ αβ ≡ Tr (cid:16) ˆΦ † ˆ c † α ˆΦ ˆ d β (cid:17) , then Eq. (A4) is equivalent to X β Q ∗ αβ U βγ (cid:0) β/ (cid:1) = X β R ∗ αβ n βγ , or, in matrix form, and observing that ˆ n = ˆ U ( β ) − ˆ n ˆ U ( β ),ˆ Q ∗ ˆ U ( β/
2) = ˆ R ∗ ˆ n = ˆ R ∗ (cid:16) ˆ U ( β ) − ˆ n ˆ U ( β ) (cid:17) = ˆ R ∗ ˆ U ( β ) − ˆ Q ∗ ˆ U (3 β/ , so that, multiplying both sides on the right by ˆ U ( − β ) wefinally getˆ R ∗ = ˆ Q ∗ (cid:16) ˆ U ( β/
2) + ˆ U ( − β/ (cid:17) = ˆ Q ∗ (cid:18)q ˆ U ( β ) + q ˆ U ( − β ) (cid:19) = ˆ Q ∗ (cid:18)r ˆ n − ˆ n + r − ˆ n ˆ n (cid:19) = ˆ Q ∗ (cid:18) q ˆ n (cid:0) − ˆ n (cid:1) (cid:19) − . We denote asˆ S ∗ = (cid:18) q ˆ n (cid:0) − ˆ n (cid:1) (cid:19) − = ˆ S T , S = ˆ S † , so thatˆ R ∗ = ˆ Q ∗ ˆ S ∗ −→ ˆ R † = ˆ S † ˆ Q † = ˆ S ˆ Q † , namely the desired resultˆ R = ˆ Q ˆ S . (A5)One can rewrite4 ˆ S − = 4ˆ n T (cid:16) − ˆ n T (cid:17) = 1 − (cid:16) − n T (cid:17) ≡ − ˆ∆ , where the matrix elements of ˆ∆ are∆ αβ = δ αβ − (cid:16) ˆΦ † ˆΦ ˆ d † β ˆ d α (cid:17) = Tr (cid:18) ˆΦ † ˆΦ h ˆ d α , ˆ d † β i(cid:19) . (A6)At equilibrium∆ (0) αβ = δ αβ (cid:16) − n (0) α (cid:17) , (A7) S (0) αβ = δ αβ / r n (0) α (cid:16) − n (0) α (cid:17) ≡ δ αβ S (0) α , (A8)are diagonal, which allow an explicit evaluation of matrixderivatives. It follows that the equilibrium renormalisa-tion matrix has elements R (0) αβ = Tr (cid:16) ˆΦ † ˆ c α ˆΦ ˆ d † β (cid:17) S (0) β ≡ Q (0) αβ S (0) β . (A9)
1. Derivatives of ˆ R We writeˆΦ = X n φ n ˆΦ n , ˆΦ † = X n φ ∗ n ˆΦ † n , where ˆΦ n is a basis set,Tr (cid:16) ˆΦ † n ˆΦ m (cid:17) = δ nm , with ˆΦ the equilibrium solution. By inspection we re-alise that ∂R αβ ∂ ˆΦ † = ˆΓ αβ h ˆΦ , ˆΦ † i ˆΦ , where the tensor ˆΓ αβ h ˆΦ , ˆΦ † i is still functional of ˆΦ andˆΦ † . Therefore ∂R αβ ∂φ ∗ n = Tr (cid:18) ˆΦ † n ˆΓ αβ h ˆΦ , ˆΦ † i ˆΦ (cid:19) . The equilibrium value is obtained by setting φ n = δ n .In particular, exploiting the fact that ˆ S is diagonal atequilibrium, the first order derivatives evaluated at equi-librium read explicitly ∂R αβ ∂φ ∗ n = ∂Q αβ ∂φ ∗ n S (0) β + X γ Q (0) αγ S (0) γ F γβ ∂ ∆ γβ ∂φ ∗ n , (A10) ∂R αβ ∂φ n = ∂Q αβ ∂φ n S (0) β + X γ Q (0) αγ S (0) γ F γβ ∂ ∆ γβ ∂φ n , (A11)while the second derivative, still calculated at equilib-rium, is ∂ R αβ ∂φ ∗ n ∂φ m = Tr (cid:18) ˆΦ † n ˆΓ αβ h ˆΦ , ˆΦ † i ˆΦ m (cid:19) +Tr ˆΦ † n ∂ ˆΓ αβ h ˆΦ , ˆΦ † i ∂φ m (cid:12)(cid:12) ˆΦ ! (A12) ≡ (cid:18) ∂ R αβ ∂φ ∗ n ∂φ m (cid:19) + (cid:18) ∂ R αβ ∂φ ∗ n ∂φ m (cid:19) , where (cid:18) ∂ R αβ ∂φ ∗ n ∂φ m (cid:19) = X γ " ∂ Q αβ ∂φ ∗ n ∂φ m S (0) β + Q (0) αγ F γβ ∂ ∆ γβ ∂φ ∗ n ∂φ m , (A13) (cid:18) ∂ R αβ ∂φ ∗ n ∂φ m (cid:19) = X γ " ∂Q αγ ∂φ ∗ n F γβ ∂ ∆ γβ ∂φ m + ∂Q αγ ∂φ m F γβ ∂ ∆ γβ ∂φ ∗ n + Q (0) αγ (cid:18) ∂ S γβ ∂φ ∗ n ∂φ m (cid:19) . (A14)The terms that appear in the above equations are ∂Q αβ ∂φ ∗ n = Tr (cid:16) ˆΦ † n ˆ c α ˆΦ ˆ d † β (cid:17) , ∂Q αβ ∂φ n = Tr (cid:16) ˆΦ † ˆ c α ˆΦ n ˆ d † β (cid:17) ,∂ ∆ αβ ∂φ ∗ n = Tr (cid:18) ˆΦ † n ˆΦ h ˆ d α , ˆ d † β i(cid:19) , ∂ ∆ αβ ∂φ n = Tr (cid:18) ˆΦ † ˆΦ n h ˆ d α , ˆ d † β i(cid:19) ,F αβ = 12 (cid:16) S (0) α S (0) β (cid:17) S (0) α + S (0) β (cid:16) − n (0) α − n (0) β (cid:17) ,∂ Q αβ ∂φ ∗ n ∂φ m = Tr (cid:16) ˆΦ † n ˆ c α ˆΦ m ˆ d † β (cid:17) , ∂ ∆ αβ ( i ) ∂φ ∗ n ∂φ m = Tr (cid:18) ˆΦ † n ˆΦ m h ˆ d α , ˆ d † β i(cid:19) , and, lastly, (cid:18) ∂ S αβ ∂φ ∗ n ∂φ m (cid:19) = (cid:16) S (0) α S (0) β (cid:17) S (0) α + S (0) β X γ " ∂ ∆ αγ ∂φ ∗ n ∂ ∆ γβ ∂φ m + ∂ ∆ αγ ∂φ m ∂ ∆ γβ ∂φ ∗ n
14 + F αγ F γβ S (0) α S (0) γ + S (0) γ S (0) β + S (0) β S (0) α (cid:16) S (0) α S (0) γ S (0) β (cid:17) . In addition ∂ R αβ ∂φ ∗ n ∂φ ∗ m = (cid:18) ∂ R αβ ∂φ ∗ n ∂φ ∗ m (cid:19) , (A15) ∂ R αβ ∂φ n ∂φ m = (cid:18) ∂ R αβ ∂φ n ∂φ m (cid:19) , (A16)where the right hand sides are obtained straightforwardlythrough Eq. (A14). The above derivatives calculated atthe equilibrium solution allow calculating the Taylor ex-pansion of ˆ R . In particular, through equations (A10) and(A11), the first order expansion isˆ R (1) = X n (cid:20) φ ∗ n ∂R αβ ∂φ ∗ n + φ n ∂R αβ ∂φ n (cid:21) , (A17) while the second order expansion mentioned in Eq. (47),is ˆ R (2) = ˆ R (2)1 + ˆ R (2)1 , (A18)where, explicitly,ˆ R (2)1 = X nm φ ∗ n φ m ∂ ˆ R∂φ ∗ n ∂φ m ! , (A19)andˆ R (2)2 = 12 X nm " φ ∗ n φ m ∂ ˆ R∂φ ∗ n ∂φ m ! + φ ∗ n φ ∗ m ∂ ˆ R∂φ ∗ n ∂φ ∗ m ! + φ n φ m ∂ ˆ R∂φ n ∂φ m ! . (A20) M. C. Gutzwiller, Phys. Rev. Lett. , 159 (1963), URL http://link.aps.org/doi/10.1103/PhysRevLett.10.159 . M. C. Gutzwiller, Phys. Rev. , A1726 (1965). H. Yokoyama and H. Shiba, Journal of the Phys-ical Society of Japan , 1490 (1987), URL http://dx.doi.org/10.1143/JPSJ.56.1490 . M. Capello, F. Becca, M. Fabrizio, S. Sorella, andE. Tosatti, Phys. Rev. Lett. , 026406 (2005), URL http://link.aps.org/doi/10.1103/PhysRevLett.94.026406 . W. F. Brinkman and T. M. Rice,Phys. Rev. B , 4302 (1970), URL http://link.aps.org/doi/10.1103/PhysRevB.2.4302 . W. Metzner and D. Vollhardt, Phys.Rev. Lett. , 324 (1989), URL http://link.aps.org/doi/10.1103/PhysRevLett.62.324 . J. B¨unemann, W. Weber, and F. Geb-hard, Phys. Rev. B , 6896 (1998), URL http://link.aps.org/doi/10.1103/PhysRevB.57.6896 . A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg,Rev. Mod. Phys. , 13 (1996). X. Deng, X. Dai, and Z. Fang, EPL (Eu-rophysics Letters) , 37008 (2008), URL http://stacks.iop.org/0295-5075/83/i=3/a=37008 . K. M. Ho, J. Schmalian, and C. Z. Wang,Phys. Rev. B , 073101 (2008), URL http://link.aps.org/doi/10.1103/PhysRevB.77.073101 . X. Y. Deng, L. Wang, X. Dai, and Z. Fang,Phys. Rev. B , 075114 (2009), URL http://link.aps.org/doi/10.1103/PhysRevB.79.075114 . G. T. Wang, Y. Qian, G. Xu, X. Dai, andZ. Fang, Phys. Rev. Lett. , 047002 (2010), URL http://link.aps.org/doi/10.1103/PhysRevLett.104.047002 . M.-F. Tian, X. Deng, Z. Fang, and X. Dai,Phys. Rev. B , 205124 (2011), URL http://link.aps.org/doi/10.1103/PhysRevB.84.205124 . Y. X. Yao, C. Z. Wang, and K. M. Ho, International Jour-nal of Quantum Chemistry , 240 (2012), ISSN 1097-461X, URL http://dx.doi.org/10.1002/qua.23238 . T. Schickling, F. Gebhard, J. B¨unemann,L. Boeri, O. K. Andersen, and W. Weber,Phys. Rev. Lett. , 036406 (2012), URL http://link.aps.org/doi/10.1103/PhysRevLett.108.036406 . N. Lanat`a, Y.-X. Yao, C.-Z. Wang, K.-M.Ho, J. Schmalian, K. Haule, and G. Kotliar,Phys. Rev. Lett. , 196801 (2013), URL http://link.aps.org/doi/10.1103/PhysRevLett.111.196801 . N. Lanat`a, Y.-X. Yao, C.-Z. Wang, K.-M. Ho, andG. Kotliar, Phys. Rev. B , 161104 (2014), URL http://link.aps.org/doi/10.1103/PhysRevB.90.161104 . G. Borghi, M. Fabrizio, and E. Tosatti,Phys. Rev. B , 125102 (2014), URL http://link.aps.org/doi/10.1103/PhysRevB.90.125102 . M.-F. Tian, H.-F. Song, H.-F. Liu, C. Wang, Z. Fang,and X. Dai, Phys. Rev. B , 125148 (2015), URL http://link.aps.org/doi/10.1103/PhysRevB.91.125148 . N. Lanat`a, Y. Yao, C.-Z. Wang, K.-M. Ho, andG. Kotliar, Phys. Rev. X , 011008 (2015), URL http://link.aps.org/doi/10.1103/PhysRevX.5.011008 . E. Runge and E. K. U. Gross, Phys.Rev. Lett. , 997 (1984), URL http://link.aps.org/doi/10.1103/PhysRevLett.52.997 . E. K. U. Gross and M. Dobson, J. F.and Peter-silka,
Density functional theory of time-dependent phe-nomena (Springer Berlin Heidelberg, Berlin, Heidel-berg, 1996), pp. 81–172, ISBN 978-3-540-49946-6, URL http://dx.doi.org/10.1007/BFb0016643 . D. Vollhardt, Rev. Mod. Phys. , 99 (1984), URL http://link.aps.org/doi/10.1103/RevModPhys.56.99 . G. Seibold and J. Lorenzana, Phys. Rev. Lett. , 2605(2001). G. Seibold, F. Becca, and J. Lorenzana,Phys. Rev. B , 085108 (2003), URL http://link.aps.org/doi/10.1103/PhysRevB.67.085108 . G. Seibold, F. Becca, P. Rubin, and J. Loren-zana, Phys. Rev. B , 155113 (2004), URL http://link.aps.org/doi/10.1103/PhysRevB.69.155113 . M. Schir´o and M. Fabrizio, Phys. Rev. B , 165105(2011). J. Bnemann, M. Capone, J. Lorenzana, and G. Sei-bold, New Journal of Physics , 053050 (2013), URL http://stacks.iop.org/1367-2630/15/i=5/a=053050 . E. v. Oelsen, G. Seibold, and J. B¨unemann,Phys. Rev. Lett. , 076402 (2011), URL http://link.aps.org/doi/10.1103/PhysRevLett.107.076402 . E. von Oelsen, G. Seibold, and J. Bnemann,New Journal of Physics , 113031 (2011), URL http://stacks.iop.org/1367-2630/13/i=11/a=113031 . M. Schir´o and M. Fabrizio, Phys. Rev. Lett , 076401(2010). G. Kotliar and A. E. Ruckenstein,Phys. Rev. Lett. , 1362 (1986), URL http://link.aps.org/doi/10.1103/PhysRevLett.57.1362 . F. Lechermann, A. Georges, G. Kotliar, and O. Parcollet,Phys. Rev. B , 155102 (2007). M. Fabrizio, in
New Materials for Thermoelectric Ap-plications: Theory and Experiment , edited by V. Zlaticand A. Hewson (Springer Netherlands, 2013), NATOScience for Peace and Security Series B: Physics andBiophysics, pp. 247–273, ISBN 978-94-007-4983-2, URL http://dx.doi.org/10.1007/978-94-007-4984-9_16 . M. Fabrizio, Phys. Rev. B , 165110 (2007), URL http://link.aps.org/doi/10.1103/PhysRevB.76.165110 . N. Lanat`a, P. Barone, and M. Fabrizio,Phys. Rev. B , 155127 (2008), URL http://link.aps.org/doi/10.1103/PhysRevB.78.155127 . N. Lanat`a, H. U. R. Strand, X. Dai, and B. Hells-ing, Phys. Rev. B , 035133 (2012), URL http://link.aps.org/doi/10.1103/PhysRevB.85.035133 . T. Jolicoeur and J. C. Le Guillou,Phys. Rev. B , 2403 (1991), URL http://link.aps.org/doi/10.1103/PhysRevB.44.2403 . R. Raimondi and C. Castellani, Phys.Rev. B , 11453 (1993), URL http://link.aps.org/doi/10.1103/PhysRevB.48.11453 . M. Sandri, M. Schir´o, and M. Fabrizio,Phys. Rev. B , 075122 (2012), URL http://link.aps.org/doi/10.1103/PhysRevB.86.075122 . T. Li, P. W¨olfle, and P. J. Hirschfeld,Phys. Rev. B , 6817 (1989), URL http://link.aps.org/doi/10.1103/PhysRevB.40.6817 . S. D. Huber and A. R¨uegg, Phys. Rev. Lett. , 065301(2009). A. R¨uegg, S. D. Huber, and M. Sigrist, Phys. Rev. B81