A f-sum rule for the conductivity

# Optical and DC conductivity of the two-dimensional Hubbard model in the pseudogap regime and across the antiferromagnetic quantum critical point, including vertex corrections

## Abstract

The conductivity of the two-dimensional Hubbard model is particularly relevant for high-temperature superconductors. Vertex corrections are expected to be important because of strongly momentum dependent self-energies. To attack this problem, one must also take into account the Mermin-Wagner theorem, the Pauli principle and crucial sum rules in order to reach non-perturbative regimes. Here, we use the Two-Particle Self-Consistent approach that satisfies these constraints. This approach is reliable from weak to intermediate coupling. A functional derivative approach ensures that vertex corrections are included in a way that satisfies the f sum-rule. The two types of vertex corrections that we find are the antiferromagnetic analogs of the Maki-Thompson and Aslamasov-Larkin contributions of superconducting fluctuations to the conductivity but, contrary to the latter, they include non-perturbative effects. The resulting analytical expressions must be evaluated numerically. The calculations are impossible unless a number of advanced numerical algorithms are used. These algorithms make extensive use of fast Fourier transforms, cubic splines and asymptotic forms. A maximum entropy approach is specially developed for analytical continuation of our results. These algorithms are explained in detail in appendices. The numerical results are for nearest neighbor hoppings. In the pseudogap regime induced by two-dimensional antiferromagnetic fluctuations, the effect of vertex corrections is dramatic. Without vertex corrections the resistivity increases as we enter the pseudogap regime. Adding vertex corrections leads to a drop in resistivity, as observed in some high temperature superconductors. At high temperature, the resistivity saturates at the Ioffe-Regel limit. At the quantum critical point and beyond, the resistivity displays both linear and quadratic temperature dependence and there is a correlation between the linear term and the superconducting transition temperature. A hump is observed in the mid-infrared range of the optical conductivity in the presence of antiferromagnetic fluctuations.

PACS number

## I Introduction

The calculation of transport quantities in strongly correlated electron systems is particularly challenging, but is a necessary step to make contact with a wide class of experiments. Even for the simplest model, namely the single-band Hubbard model, this is a formidable task. Taking up the challenge is all the more important for the two-dimensional case, that acts as the minimal model for the high temperature cuprate superconductors,(1) layered organic superconductors, (2) and a number of other materials.

Even in cases where one has a good handle on the single-particle Green’s function, the difficulty of calculating transport in the 2D Hubbard model stems from the fact that one cannot neglect the effect of vertex corrections when strong momentum-dependent correlations are present. Those vertex corrections are the analog of the self-energy, but for the two-particle response functions. When vertex corrections are not included, conservation laws can be violated and results inaccurate. In the case of small finite systems tractable by exact diagonalization or quantum Monte Carlo calculations (QMC), the correlation function is directly evaluated and vertex corrections are not an issue. However, those results are more relevant for finite frequency conductivity and strong coupling, where correlations are mainly local.(3); (4); (5); (6); (7); (8)

Consider for example the electrical conductivity for the two-dimensional Hubbard model. For the infinite system, optical and DC conductivity calculations have been performed without vertex corrections using Dynamical Mean-Field Theory (DMFT) (9) and Cellular-DMFT (CDMFT) (10). Those calculations have also been done with the composite operator method (11) but vertex corrections cannot all be taken into account. For the model, the strong coupling limit of the Hubbard model, a number of approaches have been used, in particular the extended dynamical cluster approximation (12); (13), but vertex corrections (14); (15) have been neglected. However, recent optical conductivity calculations for the Hubbard model with the dynamical cluster approximation (DCA) took vertex corrections into account(16); (17). The effects were found to be important only at high frequency. Despite this recent advance, the calculation of vertex corrections with DCA or CDMFT, considered the best available ones at strong coupling, (18); (19); (20) is still an open problem.

At weak coupling, the Boltzmann equation offers a tractable approach that satisfies conservation laws when one includes scattering-in terms. For linear response, its variational formulation has been used for example to compute the effect of spin fluctuations within the self-consistent renormalized approach (21) and also to investigate the resistivity near the quantum critical point in the clean (22) and disordered cases.(23); (24) The drawback of this approach is that it assumes the existence of quasiparticles and that this assumption is not valid in two dimensions, especially near the pseudogap regime and quantum critical point. Green’s function approaches that do not assume quasiparticles are preferable. Hence, some resistivity calculations without vertex corrections were done with the T-matrix approximation (25) and with the fluctuation-exchange (FLEX) approximation(26). Other FLEX calculations take into account some vertex corrections due to spin and charge fluctuations (27); (28). In these works, only the antiferromagnetic Maki-Thompson (MK) diagrams were included. The antiferromagnetic analogs of the Aslamasov-Larkin (AL) diagrams were neglected, claiming that the latter are negligible. We will see that in the DC case this is true only in the presence of particle-hole symmetry. A review on those calculations is given in Ref. (29). In Ref. (30) that takes into account superconducting fluctuations, AL diagrams are taken into account for superconducting fluctuations but they are neglected for the spin fluctuations. A recent calculation using field-theory methods of the conductivity at the quantum critical point (31) includes all vertex corrections but only at . The renormalized classical regime where a pseudogap appears has also been considered but focussing only on the hot spots (32) or neglecting the AL contribution (33). There are also analytical results for the conductivity with vertex corrections using Fermi liquid theory (34); (35).

Despite all these results, the electrical resistivity at weak to intermediate coupling has not been reliably computed for all dopings and temperatures because particle-hole symmetry, where the AL term vanishes, does not generally hold and because there are regimes where Fermi liquid theory is no-longer valid. Fermi-liquid theory breaks down in the pseudogap regime, near the antiferromagnetic quantum critical point and in the Ioffe-Regel limit. Hence, in this paper, we extend the Two-Particle Self-Consistent (TPSC) approach(36); (37); (38); (39) to include the effect of vertex corrections in the calculation of the resistivity and optical conductivity of the one-band, square lattice, nearest-neighbor two-dimensional Hubbard model for weak to intermediate coupling. This regime corresponds to values of the interaction strength below the critical value for the Mott transition. We present numerical results as examples and discuss possible links with experiments on cuprates. In particular, we consider the origin of the mid-infrared hump in the electron-doped materials, the Ioffe-Regel limit, insulating behavior in the pseudogap regime and the link between linear resistivity, quantum critical behavior and superconductivity.

The TPSC approach has the following strengths that make it a good choice for the present purposes. In two dimensions, the Mermin-Wagner theorem(40); (41) prevents the occurrence of antiferromagnetic long-range order at finite temperature. Not many theories can handle that constraint. Because long-range order is prohibited, there is a wide range of temperatures where there are huge antiferromagnetic [or spin-density wave (SDW)] fluctuations in the paramagnetic state. It is in this regime that a fluctuation-induced pseudogap can appear.(42); (37); (43); (44) The standard way to treat fluctuations in many-body theory, the Random Phase Approximation (RPA), leads instead to long-range order and misses this effect. The RPA also violates the Pauli principle in an important way.(36) The FLEX(45); (46) approximation and the self-consistent renormalized theory of Moriya-Lonzarich (47); (48); (49) satisfy the Mermin-Wagner theorem but they do not satisfy the Pauli principle and consistency between one and two-particle quantities. Strengths and weaknesses of these approaches are discussed further in Refs. (37); (38). Weak coupling renormalization group approaches become uncontrolled when the antiferromagnetic fluctuations begin to diverge (50); (51); (52); (53). Other approaches include the effective spin-Hamiltonian approach (54). The TPSC approach does not assume a Migdal theorem for spin fluctuations and Kanamori-Brueckner renormalization of the bare interaction is included without adjustable parameter. The conditions for the appearance of a pseudogap induced by antiferromagnetic fluctuations have been verified experimentally in electron-doped cuprates.(55) In addition to the above theoretical considerations, the TPSC approach has been extensively benchmarked against Quantum Monte Carlo calculations in regimes where the latter is available.(36); (37); (56); (43); (57); (58); (39) The agreement between the one-particle spectral function of the approach and QMC calculations in the pseudogap regime is remarkable(43). The TPSC approach however fails when temperature is too far below that where the pseudogap appears.

In the TPSC method, the calculation proceeds in two steps. At the first step, spin-spin and charge-charge correlation functions are obtained with irreducible vertices that are determined self-consistently. That is the origin of the name of the approach. At the second level of approximation, a non-trivial self-energy that is consistent with the spin and charge fluctuations and that can explain fluctuation-induced momentum-dependent pseudogaps is then calculated. The charge-charge correlation function at the first level of approximation satisfies the f-sum rule with the Green’s function at the same level, but it misses lifetime effects necessary to obtain non-trivial conductivity. What is needed is a calculation of the current-current correlation function that includes Green’s functions dressed at the second level of approximation with the corresponding irreducible vertices. What has been missing up to now is an expression for the corresponding irreducible vertices. Following Baym and Kadanoff,(59); (60) here we use a functional derivative approach to obtain irreducible vertices that satisfy conservation laws. We check that the f-sum rule is then satisfied with the Green’s function obtained at the second-level of approximation. For the conductivity, we show that, not only the Maki-Thompson-type contributions from spin-density wave (SDW) fluctuations, but also the Aslamasov-Larkin contributions are important. The latter have a dramatic effect in the pseudogap regime.

The paper is structured as follow. The next section contains the details of the methodology and is divided into five subsections: II.1, the model; II.2, the conductivity in linear response theory; II.3, a derivation of the TPSC approach for the spin and charge response functions and the one-particle self-energy; II.4, the conductivity in the TPSC approach; and finally II.5, a description of the numerical algorithms that were used to calculate the expression given in subsection II.4. Section III presents the results for the system considered, followed by a discussion and a conclusion in sections IV and V, respectively. Also, some useful derivations, such as those for the conductivity, the f-sum rule, Ward identities, are given in appendices along with details of algorithms that are of more general applicability, such as calculating response functions with Fast-Fourier transforms and cubic splines and analytical continuation of numerical data.

## Ii Methodology

We first define the model, recall the conductivity formula and introduce the TPSC method in the functional derivative formalism. This approach allows us, in the fourth subsection, to derive the conductivity formula including vertex corrections. The last subsection briefly describes the numerical algorithms that we implemented. Those are detailed in appendices. In this section, we use units where , and , being the lattice parameter. In section III, we use both those units and physical units to allow comparison with typical experimental values.

### ii.1 Model

Our model is the two-dimensional Hubbard Hamiltonian in the presence of an electromagnetic field that is treated classically. With the usual Peierls substitution, we have

 H(t)=−∑ijσtijc†iσcjσe−i∫jidrij⋅A(r,t)+U∑ini↑ni↓, (1)

where , are the hopping matrix elements between the sites of the lattice, destroys a particle with spin at site and creates a particle at site , is the vector potential, is the on-site repulsion energy and is the spin particle number operator at site . Note that it is perfectly general to use only the vector potential to represent the field, since a scalar potential can always be removed using the proper gauge transformation. The form (1) for the Hamiltonian is justified by gauge invariance. Further discussion of the Peierls substitution may be found in Ref. (61).

### ii.2 Conductivity in linear response theory

This derivation of the linear response result will allow us to set the notation. To obtain the conductivity, we first need the expression for the current operator. In the direction, for example, it is given by

 jSx(r,t)=−δHδAx(r,t), (2)

where the superscript indicates that is in Shrödinger representation, despite its dependance on . If we apply this definition to Eq. (1), keep terms up to linear order in the vector potential, assuming that , with varying slowly on the scale of a lattice spacing so that

 ∫jidxijAx(r,t)≈xij2(Ax(ri,t)+Ax(rj,t)), (3)

where is the component of the vector , we obtain,

 jSx(rl,t)=i2∑δσδxtδ(c†lσcl−δ,σ+c†l+δ,σcl,σ)−12Ax(rl,t)∑δσδ2xtδ(c†lσcl−δ,σ+c†l+δ,σclσ), (4)

where is the projection along of the vector between neighbors. is the corresponding hopping matrix element. For a uniform electric field we can take the vector potential independent of position, which means that we need only the component of the current

 Extra open brace or missing close brace (5)

where is the dispersion relation and the number of lattice sites or wave vectors in the Brillouin zone. For the following we need to define the paramagnetic current,

 jpx=−1N∑kσ∂ϵk∂kxc†kσck,σ (6)

and the diamagnetic current,

 jdx(t)=−Ax(t)1N∑kσ∂2ϵk∂k2xc†kσck,σ, (7)

so that .

According to linear response theory, the frequency dependent current in response to the field is

 Missing or unrecognized delimiter for \right (8)

where is the Fourier transform of

 jx(t)=U†(t,−∞)jSx(t)U(t,−∞), (9)

where is the time evolution operator for the Hamiltonian , Eq.(1),

 Missing or unrecognized delimiter for \left (10)

is the current-current correlation function and the notation stands for the interaction representation of the operator , namely, , where is the Hamiltonian Eq.(1) with . Also, means an equilibrium average, namely, in the system .

Since , the equilibrium average of the current operator in the interaction representation is given by the diamagnetic term

 ⟨^jx(t)⟩=⟨^jdx(t)⟩=⟨jdx(t)⟩=−Ax(t)1N∑kσ∂2ϵk∂k2x⟨c†kσck,σ⟩. (11)

Defining

 ⟨kx⟩=−1N∑k∂2ϵk∂k2x⟨nk⟩, (12)

where , we have

 ⟨^jx(t)⟩=⟨kx⟩Ax(t) (13)

or, in frequency,

 ⟨^jx(ω)⟩=⟨kx⟩Ax(ω), (14)

and Eq.(8) therefore becomes

 ⟨jx(ω)⟩=[⟨kx⟩+χjxjx(ω)]Ax(ω). (15)

Note that, in the case where the Hamiltonian has only nearest-neighbor hoppings, is proportional to the local kinetic energy, (8) but in general it cannot be regarded as such as is clear from Eq.(12).

To find the conductivity, it suffices to relate the electric field to the vector potential through

 Ex(ω)=i(ω+iη)Ax(ω), (16)

where is positive and infinitesimal. Thus, the current is related to the electric field by

 ⟨jx(ω)⟩=⟨kx⟩+χjxjx(ω)i(ω+iη)Ex(ω) (17)

and finally, since by definition, the expression for the optical conductivity in linear response theory is

 σxx(ω)=⟨kx⟩+χjxjx(ω)i(ω+iη). (18)

In this work, we calculate only the real part of , the dissipative part, and the expression that we use in practice is

 Reσxx(ω)=χ′′jxjx(ω)ω, (19)

which is derived in Appendix A from Eq.(101) to Eq.(103). If one is interested in the imaginary part, it can be obtained using Kramers-Krönig relations.

### ii.3 Two-Particle Self-Consistent approach

In the TPSC approach, one- and two-particle Green’s functions for the Hubbard model are calculated in a non-perturbative way. The approach enforces conservation laws, key sum rules and the Pauli principle. It was shown, from benchmarks with quantum Monte Carlo results (43); (37); (36); (56); (62); (58), to be accurate within a few percent for interaction strengths up to about . We will derive the TPSC approach below, but the reader can also resort to Refs. (37), (38),(58) and (39) for a more detailed discussion of the approach itself and a comparison with other approaches.

In this subsection, we present the key equations for the theory. The following subsection contains details of the derivation. We use the short-hand notation for space and imaginary time coordinates and for reciprocal space and Matsubara frequency coordinates.

The spin and charge response functions must be computed first from the expressions,

 χsp(q)=1N∫dτeiωnτ⟨TτSz(q,τ)Sz(−q)⟩=χ0(q)1−Usp2χ0(q) (20)

and

 Missing or unrecognized delimiter for \right (21)

where is the imaginary time-ordering operator and is the Lindhard function, given by

 χ0(q)=−2TN∑kG(1)(k+q)G(1)(k). (22)

Here, corresponds numerically to a non-interacting Green’s function because the initial approximation for the self-energy, that is given below in Eq.(49), is constant and the chemical potential is adjusted to obtain the correct filling. This point will become clear in the next subsection. Note that we assume here that the system is paramagnetic so that the spin index will often be omitted and the sum over it replaced by a factor of two, as in Eq.(22). In Eq.(20) and Eq.(21), the parameters and are the spin and charge irreducible vertices, respectively. First, is defined by

 Usp=U⟨n↑(1)n↓(1)⟩⟨n↑(1)⟩⟨n↓(1)⟩ (23)

(this definition is used for hole doping, electron doping will be discussed below), so that it can be determined from the fluctuation-dissipation theorem

 Missing or unrecognized delimiter for \left (24)

where we have used the Pauli principle . Note that all quantities on right-hand side of Eq.(24) are local ones. Then, once double occupancy is known, can be determined from

 TN∑qχch(q)=⟨n⟩+2⟨n↑n↓⟩−⟨n⟩2. (25)

We also call Eq.(24) and Eq.(25) the local spin and local charge sum rules.

Finally, with the response functions (20) and (21) we can obtain the one-particle self-energy,

 Σ(2)σ(k)=Un−σ+U8TN∑q[3Uspχsp(q)+Uchχch(q)]G(1)σ(k+q) (26)

whose form will be explained below. Expressions (20), (21) and (26) are the basic TPSC equations.

#### First step: spin and charge susceptibilities

First let us derive expressions (20) and (21) for the spin and charge susceptibilities. In the following, we use Einstein’s convention for the sums (or integrals), namely that an index appearing twice or more in an expression is summed over lattice sites and integrated over imaginary time. An overbar helps clarify which indices are involved. The approach follows the Martin-Schwinger techniques (63) described in Kadanoff and Baym’s book (60).

It is convenient to introduce a “source field” that couples to one-particle excitations in the system. It allows us to easily obtain correlation functions from functional derivatives. The source field can be set to zero at the end of the calculation. The source-field dependent Green’s function in the grand canonical ensemble is then given by

 Gσ(1,2;{ϕ})=−Tr[e−βKTτe−c†¯σ(¯1)ϕ¯σ(¯1,¯2)c¯σ(¯2)cσ(1)c†σ(2)]Tr[e−βKTτe−c†¯σ(¯1)ϕ¯σ(¯1,¯2)c¯σ(¯2)]=−⟨Tτcσ(1)c†σ(2)⟩ϕ (27)

where , being the chemical potential and the number operator. We have also used the notation which means that the average is taken with the source field turned on. Response functions can be obtained from functional derivatives of with respect to since

 δGσ(1,2;{ϕ})δϕσ′(3,4)=Gσ′(4,3;{ϕ})Gσ(1,2;{ϕ})−⟨Tτc†σ′(3)cσ′(4)c†σ(2)cσ(1)⟩ϕ. (28)

We also have

 Missing or unrecognized delimiter for \right (29)

where we have used spin rotation symmetry to obtain the last line from the previous one. For the charge response function, the corresponding expression is

 Missing or unrecognized delimiter for \right (30)

According to Eq. (28), the last two results Eqs. (29) and (30) can be written as

 χch/sp(1,2)=−2(δG↑(1,1+;{ϕ})δϕ↑(2+,2)±δG↑(1,1+;{ϕ})δϕ↓(2+,2))∣∣∣{ϕ}=0, (31)

where the expressions with the plus and minus sign correspond respectively to the charge and spin response functions. Here, , where is positive and infinitesimal. For the remainder of this section, we implicitly assume that derivatives with respect to are evaluated at .

To obtain integral equations for the response functions, one begins with

 Gσ(1,¯1)G−1σ(¯1,2)=δ(1−2). (32)

Taking the functional derivative of this equation with respect to , taking , multiplying on the right by and summing over , we obtain

 δGσ(1,2)δϕσ′(3,4)=−Gσ(1,¯1)δG−1σ(¯1,¯2)δϕσ′(3,4)Gσ(¯2,2). (33)

On the other hand, Dyson’s equation in the presence of the field reads

 G−1σ(1,2)=G(0)−1σ(1,2)−ϕσ(1,2)−Σσ(1,2), (34)

where is the non interacting Green’s function and is the self-energy, so that

 δG−1σ(1,2)δϕσ′(3,4)=−δ(1−3)δ(2−4)δσσ′−δΣσ(1,2)δϕσ′(3,4), (35)

and therefore, from (33),

 δGσ(1,2)δϕσ′(3,4)=Gσ(1,3)Gσ(4,2)δσσ′+Gσ(1,¯1)δΣσ(¯1,¯2)δϕσ′(3,4)Gσ(¯2,2). (36)

Following Luttinger and Ward(64), is a functional of and , we find, applying the chain rule, that

 δGσ(1,2)δϕσ′(3,4)=Gσ(1,3)Gσ(4,2)δσσ′+Gσ(1,¯1)δΣσ(¯1,¯2)δG¯σ(¯3,¯4)δG¯σ(¯3,¯4)δϕσ′(3,4)Gσ(¯2,2). (37)

This is the analog of the Bethe-Salpeter equation for the particle-hole channel. Defining the particle-hole irreducible vertex

 Γσσ′(1,2;3,4)=δΣσ(1,2)δGσ′(3,4), (38)

 δGσ(1,2)δϕσ′(3,4)=Gσ(1,3)Gσ(4,2)δσσ′+Gσ(1,¯1)Γσ¯σ(¯1,¯2;¯3,¯4)δG¯σ(¯3,¯4)δϕσ′(3,4)Gσ(¯2,2). (39)

The spin and charge response functions in Eq. (31) are special cases of the more general functions

 χ±(1,2;3,4)=−2(δG↑(1,2;{ϕ})δϕ↑(3,4)±δG↑(1,2;{ϕ})δϕ↓(3,4)). (40)

It is straightforward to show, from the Bethe-Salpeter equation Eq. (39) and spin-rotational invariance, that

 χ±(1,2;3,4)=−2G(1,3)G(4,2)±Γch/sp(¯1,¯2;¯3,¯4)G(1,¯1)G(¯2,2)χ±(¯3,¯4|3,4) (41)

where

 Γch/sp(¯1,¯2;¯3,¯4)=δΣ↑(¯1,¯2)δG↓(¯3,¯4)±δΣ↑(¯1,¯2)δG↑(¯3,¯4) (42)

and .

Up to now, all the results are exact. The first step in the TPSC method is to obtain a first approximation for the self-energy that will be used to obtain irreducible vertices. First, let us rewrite Dyson’s equation Eq. (34) with zero source field in the form,

 G(0)−1σ(1,¯3)Gσ(¯3,2)=δ(1−2)+Σσ(1,¯3)Gσ(¯3,2). (43)

Then, note that the equation of motion for the Green’s function reads

 ∑l[(−∂∂τ+μ)δil+til]Gσ(l−j,τ)=δ(τ)δij−U⟨Tτni~σ(τ)ciσ(τ)c†jσ⟩. (44)

or, in compact form,

 G(0)−1σ(1,¯3)Gσ(¯3,2)=δ(1−2)−U⟨Tτn~σ(1)cσ(1)c†σ(2)⟩, (45)

where . By comparing Eqs.(43) and (45) one concludes that

 Σσ(1,¯3)Gσ(¯3,2)=−U⟨Tτn~σ(1)cσ(1)c†σ(2)⟩. (46)

Now, comes our first approximation. We replace the last expression Eq. (46) by

 Σσ(1,¯1)Gσ(¯1,2)≈Ug↑↓(1)G~σ(1,1+)Gσ(1,2), (47)

where, for the hole-doped case,

 g↑↓(1)=⟨n↑(1)n↓(1)⟩⟨n↑(1)⟩⟨n↓(1)⟩. (48)

For the electron-doped case, one replaces in this expression by . Thus, when the lattice is bipartite, particle-hole symetry of the phase diagram is preserved. It also gives a better agreement with quantum Monte Carlo results when the lattice is not bipartite. Those two different approximation are equivalent to assume that the approximation (47) is made for holes when hole density is smaller than particle density and for particles otherwise. Now, substituting the definition (48) in Eq. (47), it is clear that one recovers the exact result for the self-energy Eq. (46) in the special case . The approximation (65); (36) is justified if it is correct to assume that the four-point correlation function factorizes only when points and are different.

Within this approximation, the self-energy can be obtained by multiplying Eq. (47) by from the right and summing over , to obtain

 Σ(1)σ(1,2)=Ug↑↓(1)G~σ(1,1+)δ(1−2), (49)

which is the self-energy ansatz that is used to define the Green’s function in the Lindhard function (22), from which are defined the TPSC susceptibilities (20) and (21). Note that once the source field is turned off and translational invariance is restored, becomes independent of position and its Fourier transform is simply a constant that can be absorbed in the chemical potential in , so that this Green’s function is in practice a non-interacting one.

Given the self-energy, as well as the result

 δg↑↓(1)δGσ(2,3)=δg↑↓(1)δG−σ(2,3), (50)

valid in the paramagnetic phase, and the definition of the spin vertex, Eq.(42), we obtain the spin vertex

 Γsp(1,2;3,4)=Uspδ(1−3)δ(1+−4)δ(1−−2) (51)

where . It is possible to obtain an analytical expression for to compute the charge vertex from the definition (42) and the ansatz (49), but to this day the most successful approach has been to approximate this functional derivative as local, i.e. proportional to , which leads to

 Γch(1,2;3,4)=Uchδ(1−3)δ(1+−4)δ(1−−2). (52)

Substituting in the Bethe-Salpeter equation for the susceptibilities, Eq.(41), we find the corresponding TPSC expressions for the spin (20) and charge (21) susceptibilities. They suffice to determine also .

#### Second step: improved self-energy

Collective modes are generally less influenced by details of the single-particle properties than the other way around since collective modes depend more strongly on general principles like conservation laws. We thus wish now to obtain an improved approximation for the self-energy that takes advantage of the fact that we have found accurate approximations for the low-frequency spin and charge fluctuations. We begin from the general definition of the self-energy Eq.(46) obtained from Dyson’s equation.

We start with the longitudinal channel ( diagonal in spin indices) and use the corresponding expression for the correlation function in terms of response function Eq.(28). In that case, the right-hand side of the general definition of the self-energy Eq.(46) may be written as

 Σσ(1,¯1)Gσ(¯1,2)=−U[δGσ(1,2)δϕ−σ(1+,1)−G−σ(1,1+)Gσ(1,2)]. (53)

The last term is the Hartree-Fock contribution, which gives the exact result for the left-hand side in the limit .(37) The term is thus a contribution to lower frequencies and it comes from the spin and charge fluctuations. Right-multiplying the last equation by , replacing the lower energy part by its general expression in terms of irreducible vertices, Eq.(37), and taking and on the right-hand side to be the first level approximations and , respectively, we find

 Σ(2)σ(1,2)=UG(1)−σ(1,1+)δ(1−2)−UG(1)σ(1,¯3)δΣ(1)σ(¯3,2)δG(1)¯σ(¯4,¯5)δG(1)¯σ(¯4,¯5)δϕ−σ(1+,1). (54)

If we expand the sum over spins on the right-hand side to express the irreducible vertices in terms of their spin and charge versions Eqs.(42) we find, after using the TPSC vertices, Eqs.(51,52),

 Σl(2)σ(k)=Un−σ+U4TN∑q[Uspχsp(q)+Uchχch(q)]G(1)σ(k+q). (55)

There is, however, an ambiguity in obtaining the self-energy formula. Indeed, we can obtain an expression for the self-energy by using a transverse source field in Eq.(27). By taking functional derivatives of with respect to this , we first obtain the transverse spin correlation functions and . Then, using a derivation analogous to that for the longitudinal case, (38) we find

 Σt(2)σ(k)=Un−σ+U2TN∑kUspχsp(q)G(1)−σ(k+q). (56)

During the derivation, was used, taking spin rotational invariance into account.

The two previous expressions for the self-energy clearly show that our approximations for the fully reducible vertex does not preserve crossing symmetry, that is the symmetry under the exchange of two particles or two holes. To improve our approximation and restore crossing symmetry,(43) we average the two expressions (55) and (56), which gives the final result Eq.(26) that we use in the rest of this paper. It turns out that this “symmetric” expression of the self energy gives a better agreement with quantum Monte Carlo results.(43) In addition, one verifies numerically that the exact sum rule (Appendix A in Ref. (37))

 −∫dω′πImΣRσ(k,ω′)=U2n−σ(1−n−σ), (57)

determining the high-frequency behavior of , is satisfied to a higher degree of accuracy with the symmetrized self-energy expression Eq.(26).

#### Internal accuracy checks

Our final expression for the self-energy Eq.(26), is in principle an improvement over the constant self-energy entering the calculation of the susceptibilities. But clearly the approach is not one-particle self-consistent. All the Green’s functions entering the right-hand side of Eq.(26) are evaluated with , which has a constant self-energy. The advantage is that all quantities, including the vertices, on the right-hand side of Eq.(26) are computed at the same level of approximation. In fact, one can miss some important physics if this is not the case (37).

Apart from comparisons with quantum Monte Carlo (QMC) calculations, we can check the accuracy in other ways. For example, the f-sum rule,

 ∫dωπωχ′′ch(q,ω)=1N∑kσ(ϵk+q+ϵk−q−2ϵk)⟨nkσ⟩, (58)

is exactly satisfied at the first level of approximation (i.e. with on the right-hand side) and the charge susceptibility obtained with . Suppose that on the right-hand side of that equation, one uses obtained from instead of the Fermi function. In cases where the agreement with QMC calculations is good, one should find that the right-hand side does not change by more than a few percent.

When we are in the Fermi liquid regime, another way to verify the accuracy of the approach is to verify if the Fermi surface obtained from satisfies Luttinger’s theorem very closely.

Finally, the consistency relation between one- and two-particle quantities ( and ), Eq. (46), should be satisfied exactly in the Hubbard model. In the special case where , this relation can be written as

 12Tr(ΣG)=limτ→0−TN∑ke−iknτΣσ(k)Gσ(k)=U⟨n↑n↓⟩. (59)

In standard many-body books (66), this expression is encountered in the calculation of the free energy through a coupling-constant integration. In the TPSC approach, it is not difficult to show (Appendix B of Ref. (37)) that the following equation

 12Tr(Σ(2)G(1))=U⟨n↑n↓⟩ (60)

is satisfied exactly with the self-consistent obtained from the sum rule (24). An internal accuracy check consists in verifying by how much differs from Again, in regimes where we have agreement with Quantum Monte Carlo calculations, the difference is only a few percent.

The above relation between and gives us another way to justify our expression for , Eq.(26). Suppose one starts from (55) to obtain a self-energy that contains only the longitudinal spin fluctuations and the charge fluctuations, as was done in the first papers on the TPSC approach(36). One finds that the spin part and the charge part each contribute an amount to the consistency relation Eq.(60). Similarly, if we work only in the transverse spin channel (43); (38), we find that each of the two transverse spin components also contributes to . Hence, averaging the two expressions also preserves rotational invariance since each spin component contributes equally to Eq.(60). Note that Eq.(26) for is different from so-called Berk-Schrieffer type expressions (67) that contains only bare vertices and do not satisfy (Appendix E in Ref.  (37)) the consistency condition between one- and two-particle properties, Eq.(59).

### ii.4 Conductivity in the Two-Particle Self-Consistent approach

To compute Re  from the Kubo formula Eq.(19), we have to obtain the current-current correlation function given by Eq. (10), which contains the component of the paramagnetic current. From the general expression for the current Eq.(4), we have

 jpx=∑ljpx(rl)=i∑δσδxtδ∑lc†lσcl−δ,σ. (61)

Since, the actual calculation is done in Matsubara frequency instead of the real-frequency definition of , Eq. (10), we use

 Missing or unrecognized delimiter for \right (62)

where and , with integer, is a bosonic Matsubara frequency. Substituting (61) in this expression gives

 ⟨Tτ^jpx(τ)^jpx(τ′)⟩=−∑ilδ1δ2σ1σ2δx1δx2tδ1tδ2×⟨Tτc†iσ(τ)ci−δ1,σ(τ)c†lσ(τ′)cl−δ2,σ(τ′)⟩. (63)

If we substitute the functional derivative expression for the correlation function, Eq.(28), in this equation, we obtain

 Missing or unrecognized delimiter for \right (64)

where have used the notation and . Note that the first term on the right-hand side of Eq. (28) does not contribute to the sum since . Now, summing over spin indices and using spin rotational invariance, we obtain

 Missing or unrecognized delimiter for \right (65)

where we have used the definition Eq. (40) for the generalized susceptibility .

The most general way of thinking about the last result is that it comes from a functional derivative of the current with respect to an applied vector potential representing the electric field. The remaining of the derivation, that we give in the rest of this subsection, is based on the idea that we should evaluate this functional derivative in a systematic way to obtain a conserving approximation, as shown by Baym (59). We will keep working with the source field but it is useful to remember that within simple prefactors, it is equivalent to working with the vector potential as the source field.

The susceptibility is defined in terms of a functional derivative in Eq. (40). That functional derivative leads to a Bethe-Salpeter like equation (36) that contains two different types of terms

 χ+(1,2|3,4)=−2G(1,3)G(4,2)−2G(1,¯1)(δΣσ(¯1,¯2)δϕσ(3,4)+δΣσ(¯1,¯2)δϕ−σ(3,4))G(¯2,2). (66)

The product comes from the explicit dependence of on the source field while the last one comes from the implicit dependence of the self-energy on that source field. The product is what leads to the ”bubble” in standard conductivity calculations. The other term is the vertex correction.

Up to now, everything is exact in the present subsection. To work within the TPSC approach, it suffices to use everywhere the results obtained at the second level of approximation, namely and since at the first level of approximation the conductivity is infinite. Our explicit formula for Eq. (26) is a functional of . We can thus write

 δΣ(2)σ(¯1,¯2)δϕσ′(3,4)=δΣ(2)σ(¯1,¯2)δG(1)¯σ(¯3,¯4)δG(1)¯σ(¯3,¯4)δϕσ′(3,4). (67)

Expanding the sum over spin, using the chain rule, spin rotation invariance and the definition of Eq. (40) we obtain

 δΣ(2)σ(¯1,¯2)δϕσ(3,4)+δΣ(2)σ(¯1,¯2)δϕ−σ(3,4)=−12⎛⎝δΣ(2)σ(¯1,¯2)δG(1)σ(¯3,¯4)+δΣ(2)σ(¯1,¯2)δG(1)−σ(¯3,¯4)⎞⎠χ(1)+(¯3,¯4|3,4) (68)

where is computed with . Substituting in the exact expression for Eq. (66), we find

 χ(2)+(1,2|3,4)=−2G(2)(1,3)G(2)(4,2)+G(2)(1,¯1)G(2)(¯2,2)⎛⎝δΣ(2)σ(¯1,¯2)δG(1)σ(¯3,¯4)+δΣ(2)σ(¯1,¯2)δG(1)−σ(¯3,¯4)⎞⎠×χ(1)+(¯3,¯4|3,4). (69)

If had been a functional of we would have had an infinite series to sum. Instead, here the series ends as we now show. First note that, from Eq.(41) and the approximation (52) for the charge vertex,

 χ(1)+(¯3,¯4|3,4)=−2G(1)(¯3,3)G(1)(4,¯4)+G(1)(¯3,¯5)G(1)(¯6,¯4)Γ