Self-Consistent Fluctuation Theory for Strongly Correlated Electron Systems\abst
A self-consistent theory for two-particle fluctuations with renormalized irreducible vertices is proposed. Using the Parquet formalism, we construct the fully antisymmetric full vertex in terms of the two-particle fluctuations in the charge, the spin and the particle-particle channels on an equal footing to satisfy the Pauli principle. The fluctuations are determined self-consistently, which are reflected into the one-particle self-energy via the Schwinger-Dyson equation. We demonstrate the application of the present theory to the impurity Anderson model and the Hubbard model on a square lattice mainly for the particle-hole symmetric case. In both models the vertex renormalization in the spin channel eliminates magnetic instabilities of mean-field theory to ensure the Mermin-Wagner theorem. The present theory gives the same critical exponents of the self-consistent renormalization theory in the quantum critical region. \kwordtwo-particle self-consistency, parquet equation, self-consistent renormalization theory, impurity Anderson model, Hubbard model, Mermin-Wagner theorem, quantum critical point
A variety of phenomena such as a metal-insulator transition[1, 2, 3, 4, 5, 6] and an anisotropic superconductivity[7, 8, 9, 10] emerges by strong electron-electron correlations. The fundamental models like the Hubbard and the Anderson lattice models have been investigated extensively using various analytical and numerical techniques. It has revealed that microscopic detail of the non-interacting band structures is capable to yield various type of metallic states and superconducting pairings in spite of simplicity of the models. Insulating phases with or without a magnetic ordering are described by the models as well.
The attempts to describe effects of strong correlations in metallic phases often begin with the Hartree-Fock (HF) and the random-phase approximation (RPA). Then, the low-order perturbation theory and the fluctuation exchange (FLEX) approximation[11, 12, 13] that includes systematically the collective fluctuations improve the treatment of correlations. They are mainly based on the approximation for the one-particle self-energy, and the latter follows the Baym and Kadanoff’s conserving arguments[14, 15] within the one-particle level. So-called the dynamical mean-field theory (DMFT) taking account of local fluctuations exactly[17, 18, 19] is also classified into the one-particle based theories. The theories mentioned above have achieved considerable success in some aspects, however, at the same time they have significant inconsistency in the two-particle level such as a violation of the Mermin-Wagner theorem[20, 21]. This is inherent from the one-particle based theories and/or the lack of (spatial) vertex corrections[22, 23]. These theories also violate the Pauli principle that is the essential ingredient behind the Fermi-liquid (FL) theory[24, 25, 26, 27]. We also note that an appropriate description for a magnetic instability especially near a quantum critical point (QCP) facilitates a better understanding of anisotropic superconductivity mostly appeared near the QCP[7, 8].
Meanwhile, the theories based on the two-particle self-consistency have been developed. The successful phenomenological method is the self-consistent renormalization (SCR) theory by Moriya and co-workers[8, 28, 29, 30, 8]. The SCR is thermodynamically consistent due to the self-consistent treatment of spin fluctuations in the mode-mode coupling term, and it gives the non-trivial exponents near the QCP. Since the SCR theory relies on the phenomenological expression for the spin fluctuation that is valid in the FL regime, a range of the applicability of the SCR is a delicate issue especially in heavy-fermion systems where a crossover to a local-moment regime occurs at elevated temperatures. Note that in the quantum critical region, the renormalization-group treatments give the same exponents of the SCR[32, 33]. The similar idea of the SCR has been used in the microscopic theory referred as the two-particle self-consistent (TPSC)[34, 35], where the charge and the spin fluctuations in the longitudinal channel are determined by enforcing the local sum rules to satisfy the Pauli principle. The double occupancy in the theory is expressed empirically in terms of the renormalized irreducible vertex used to describe the spin fluctuation, which introduces the self-consistent relation of the fluctuations. Once the two-particle fluctuations are obtained, the single-particle self-energy is computed by the way similar to the paramagnon theories.
In order to go beyond the FLEX, Bickers and the collaborators have argued the full consistency among one and two-particle levels using the parquet formalism[13, 35, 36, 37]. Note that the parquet formalism with the unrenormalized bare irreducible vertices is equivalent to the FLEX up to the two-particle level. They have approximated the irreducible vertices as the contact-type pseudo potentials at the representative point in the momentum-frequency space in order to solve the parquet equation practically. The quantitative agreements have been demonstrated by the comparison with the quantum monte carlo (QMC) data for the intermediate repulsions. Nevertheless, there exists a methodological difficulty in choosing the representative point in the 4-vector space for general models, which is usually unknown a priori. For the impurity Anderson model Janiš and the co-worker have found the semi-analytic solution for the simplified parquet equation in the critical regime of the spin fluctuations[39, 40]. They have reproduced the correct Kondo scale in the strong-coupling limit.
These findings strongly encourage to construct a reasonable microscopic theory for strongly correlated electron systems using the parquet formalism with appropriate modifications to the previous attempts. The parquet formalism is also suitable for the purpose that the two-particle self-consistency should be enforced prior to the single-particle one in view of the loop-expansion argument of the renormalization-group theory. Furthermore, the parquet equation is to be thermodynamically consistent, since it can be derived by functional derivatives of the generating free-energy functional with respect to the two-particle vertex functions. It is not a conserving approximation in the sense of Baym and Kadanoff[13, 35]. Needless to say, conserving approximations do not necessarily mean quantitatively reliable approximations. Indeed, the FLEX shows the considerable deviations from the finite-size QMC results even in the weak-coupling region[13, 34]. Moreover, one should notice that any expression for the self-energy following the Baym-Kadanoff scheme is inconsistent with the exact self-energy formula of the Schwinger-Dyson equation[13, 34, 41]. The both relations could be satisfied simultaneously, only when the exact expression for the irreducible vertices would be found.
In this paper, we propose yet another self-consistent fluctuation (SCF) theory based on the parquet formalism, which will be explained in the next section. The formulation of the SCF theory will be given in §3, where the fully antisymmetric full vertex is constructed in terms of the two-particle fluctuations in all channels to ensure the Pauli principle. Regarding the Pauli principle, the present theory is contrast to the TPSC, which satisfies it implicitly through the local sum rules. These fluctuations are determined self-consistently, and they are used to improve the one-particle self-energy in the same spirit of the TPSC. Note that the numerical cost of the present theory is much less than that of the FLEX. In §4, we demonstrate the applications to the impurity Anderson model and the Hubbard model on a square lattice as a benchmark of the present theory. In the last section, the relation to the SCR theory in the QCP region is discussed. The present theory gives the same critical exponents and it can describe a crossover to non-universal region away from the QCP without introducing any phenomenological cutoffs as in the SCR. We also argue the possible improvements of the present theory, and we summarize the paper. In Appendix we summarize the derivation of the parquet equation for the single-band model in the case of a nonlocal general two-body interaction conserving the -component of the spins.
2 Parquet Formalism
In this section, we introduce briefly the parquet formalism and set up the notational convention through it. A detailed pedagogical review for the parquet formalism can be found in the literature.
2.1 Full vertex and crossing symmetry
Let us first define an effective interaction with full vertex in the particle-particle (PP) and the particle-hole (PH) channels in the form,
where the number in the parenthesis represents a set of indices of the band (orbital), the spin, the spatial coordinate and the imaginary time. The superscript “p” indicates a quantity concerning the PP channel. The repeated indices that do not appear in the left-hand side are assumed to be summed hereafter. We define the fourier and inverse fourier transformations between the spatial coordinates and the four-vector momenta , , as
where , are the fermionic Matsubara frequencies, while is the bosonic one. We have assigned the four momenta at each legs of the vertex as in the PH channel and in the PP channel. The assignment of the variables for the full vertices are schematically shown in Fig. 1.
The antisymmetric nature of the fermion fields gives the following relations for the full vertices,
where we have introduced sets of coordinates/momenta as
for notational convenience. These relations are called the crossing symmetry, which is the direct consequence of the Pauli principle. Note that frequently used theories such as the HF-RPA and FLEX violate the crossing symmetry.
2.2 Bethe-Salpeter equation
The full vertex can be decomposed into multiple scatterings of the PP or the PH pair, respectively, in terms of the irreducible vertices, and . Namely, so-called the Bethe-Salpeter (BS) equation in the PH channel is given by
where is the lowest-order two-particle green’s function. , is irreducible with respect to a vertical cut of two green’s function lines. Here the “matrix product” is introduced as
The identity matrix is given by .
Similarly, the BS equation in the PP channel is given by
where . The factor is to avoid double counting of twisted diagrams. Exchanging the indices and in the BS equation, we can show that also satisfies the crossing symmetry (6), i.e., . Note that the irreducible vertex in one channel is reducible in another in general. The Feynman diagrams for the BS equations are given in Fig. 2.
Using the full vertices, we express two-particle green’s functions in the PP and the PH channels as
The ordinary two-particle correlation functions are obtained by setting and , i.e., and , where with being the number of lattice sites.
2.3 Parquet equation
Although the equations given in the previous subsections are formally exact, it is necessary to introduce physically motivated approximations to make them practically solvable. However, a simple approximation to the irreducible vertices often results in full vertices violating the crossing symmetry. In order to obtain a solution that maintains the crossing symmetry, it is adequate to use a set of exact relations between the irreducible vertices. This is the Parquet equation[35, 36], which is give by
Here we have introduced the fully irreducible vertex, , in all channels, which satisfies the crossing symmetry like (6). A choice of is ambiguous. In the simplest case, an anti-symmetrized bare interaction is used for , and it is often called as “basic” parquet equation.
Using the BS equations, we reexpress the full vertices as
where the crossing symmetry is apparent. It should be emphasized that resultant full vertices always satisfy the crossing symmetry, irrespective of the explicit expressions of and . The Feynman diagrams of the parquet equations are shown in Fig. 3.
The BS and the parquet equations constitute a set of the self-consistent integral equations among the full and the irreducible vertices for an arbitrary given one-particle green’s function .
2.4 Self energy in terms of full vertex
Once we obtain the full vertex, the self-energy is calculated through the Schwinger-Dyson equation as shown schematically in Fig. 4,
This equation together with the BS and the parquet equations provide a complete set of the self-consistent equations for correlated fermion systems.
2.5 Case of Hubbard-type interaction with SU(2) symmetry
In the presence of SU(2) symmetry in spin space, the two-particle quantities can be classified by the total spin of the PH and the PP pairs, and the one-particle quantities are independent of spin indices (See, Appendix in detail). Namely, the charge (even-parity) and the spin (odd-parity) vertices belong to and PH (PP) channels, respectively. For later convenience, we quote the relevant expressions in the case of the Hubbard-type bare interaction, with SU(2) symmetry.
where and .
Two-particle correlation function
where and . Note that the (pair) density operators with are used in the definition of the correlation functions, i.e.,
Basic parquet equation ()
Here is used instead of in order to keep track of the fully irreducible vertex.
where is the electron density per spin.
With this preliminary, we will discuss a practical approximation of the parquet formalism in the next section.
3 Self-Consistent Fluctuation Approximation
The parquet formalism is a reliable nontrivial renormalization scheme, where two-particle self-consistency and the crossing symmetry are respected. However, a practical solution is limited to very small systems since the central nonlinear integral equations with three 4-vector variables are too large to handle even at modern numerical computation facilities. Therefore, a practical simplification preserving essential ingredient of the parquet formalism is necessary.
In order to simplify the parquet formalism, we introduce the contact-type effective irreducible vertices in the charge, the spin, the even-parity and the odd-parity channels,
which act as the renormalized vertices on an average, and they will be determined later without losing two-particle self-consistency and the crossing symmetry. The minus sign and the factor 2 are only for the notational simplicity of the following equations. Note that the effective vertex in the odd-parity channel identically vanishes due to the crossing symmetry, . With this approximation, the formal solution of the BS equation (29) takes the form,
where we have introduced two “bosonic” variables, and . The Pauli principle can be confirmed directly by these expressions.
Substituting these vertices into the second definition of the two-particle green’s functions (33), we obtain
where the vertex corrections are defined as
Here we have defined the integrals corresponding to the Maki-Thompson processes as
Since we attempt to find a solution within the contact-type vertex approximation, the above equations make sense with an appropriate averaging procedure with respect to . We adopt a simple average over , and we finally obtain the self-consistent condition for the effective vertices,
Note that the average of , and appeared in the vertex corrections (62) reduces to the simple average,
By solving (62), (68) and (70) self-consistently for given , we obtain the renormalized vertices, , and , and simultaneously the two-particle correlation functions and the fully antisymmetric vertices via (55).
The self-energy in this approximation is schematically shown in Fig. 7. Note that when we neglect the vertex corrections, , the present formalism reduces to the FLEX.
Let us now comment on the self-consistent procedure between two-particle and one-particle levels. In the two-particle self-consistency, the one-particle propagator is an external input. It is well known that the one-particle self-consistency without using the consistent frequency-dependent irreducible vertices always results in a universal high-frequency behavior and the Hubbard satellite bands are completely smeared out[34, 39]. This is because the high-frequency asymptotic behavior sets in, not above the bandwidth but above an incorrect energy scale, , which is extensively discussed by Vilk and Tremblay, and Janiš and Augustinský.
In order to balance the use of the contact-type irreducible vertices, and , we use the HF propagator as an input for the two-particle self-consistency, and the one-particle self-energy is updated only once. Specifically, we first determine the HF propagator by adjusting the chemical potential to match the given electron density per spin, . Using , we determine , and by solving the two-particle self-consistent equations. Then, the self-energy is updated via the Schwinger-Dyson equation, and the renormalized chemical potential, , is determined by with .
Empirically, the contribution from the PP fluctuation seems to be slightly overestimated, and it yields unphysical behaviors in one-particle spectrum. By neglecting the last term in (71), the SCF works quite well as will be shown in the next section.
4 Application to Fundamental Models
In this section, we consider two fundamental models of many-body fermion systems, i.e., the impurity Anderson model and the Hubbard model on a square lattice, as applications of the SCF approximation. In the latter, the first Brillouin zone is discretized by meshes, and Matsubara frequencies are used in both models. To obtain the real-frequency spectrum, we adopt the Padé interpolation to carry out an analytic continuation.
4.1 Impurity Anderson model
The impurity Anderson model is given by
where the level is measured from , and . We use the practically infinite bandwidth of the conduction electron, so that the Hartree green’s function is given by
where . The hybridization strength is , which is taken as a unit of energy.
Let us consider the symmetric case (), in which and are pure imaginary, and , and hold.
The dependence of the effective irreducible vertices at are shown in Fig. 8. The charge (even-parity) vertex is roughly the order of , while the spin vertex is strongly renormalized as increases. In the HF-RPA theory, the Stoner condition is given by , while the spin vertex in the present theory never reaches the magnetic instability condition as it should in the essentially zero-dimensional system. Note that gives a measure of the Kondo temperature, i.e., .
The dependences of and at are shown in Fig. 9. Both the charge and the spin vertices have very weak temperature dependence, and the dependence of the susceptibilities thus mainly comes from that of .
The dependence of the static charge and spin susceptibilities, and at are shown in Fig. 10. The charge fluctuation is suppressed below , while the saturation of the spin susceptibility occurs below reflecting the singlet formation of the Kondo effect. The exact expression of the Kondo temperature in Kondo regime[45, 46] () is given by
For the intermediate coupling, , and the present theory slightly overestimates the Kondo temperature.
The dynamical spin susceptibility is shown in Fig. 11. As the temperature decreases, the intensity of transfers to the lower-energy region. The low-energy limit of is proportional to as the FL theory.
Let us now turn to the one-particle density of states, . The dependence of at is shown in Fig. 12. The central Kondo quasiparticle peak is well developed and its width becomes narrower as increases. The position of the Hubbard satellite peaks is slightly smaller than the values of the atomic limit, . The overall features of the DOS are well reproduced and is quantitatively similar to the best know results by the numerical renormalization-group method[46, 47].
The dependence of at is shown in Fig. 13. The Kondo peak begins to develop below where the static spin susceptibility also begins to saturate. Note that the zero-temperature limit of is fixed as due to the Friedel sum rule.
The quantitative success of the density of states has a common origin of the second-order perturbation theory in the symmetric case[49, 50, 51, 52]. Namely, the dominant contribution of the fluctuations in (71) is the second-order term (the first term in the curly bracket). By the same reason, the overall features of in the asymmetric case (not shown here) are similar to those of the second-order perturbation theory[51, 52]. Nevertheless, the low-order perturbation theory is inapplicable in higher-dimensional systems in principle, where the phase transitions are possible to occur. The present approach can describe potential phase transitions in higher dimensions, and at the same time the artificial instabilities of the mean-field theory can be successfully eliminated by the self-consistent treatment of the fluctuations. The concrete example of the two-dimensional Hubbard model will be shown in the next subsection.
4.2 Hubbard model on square lattice
Next, we consider the Hubbard model with the nearest-neighbor hopping,
where , with being the lattice constant. The number operator at site with the spin is given by . In the following we use . In the particle-hole symmetric case at the half filling, we have the relations, , , and , where and .
Let us discuss the symmetric case at the half filling, . Figures 14 and 15 show the and dependences of and . The overall behaviors are similar to those of the impurity Anderson model, e.g., is of the order of , is renormalized to eliminate the magnetic instability, and both of and have weak dependences.
For in Fig. 15, approaches to unity, but it never exceeds the value of the magnetic instability toward . In order to elucidate the Mermin-Wagner theorem, the dependence of the inverse staggered susceptibility is shown in Fig. 16. follows the Curie-Weiss behavior in higher temperatures, and it begins to deviate below . For later convenience, we define the fictitious Néel temperature, by extrapolating the Cruie-Weiss dependence. corresponds to “mean-field” Néel temperature, and it is regarded as the true critical temperature in the presence of a weak three dimensionality. The inverse susceptibility is then approaching to zero with the exponential dependence. The dashed lines in the figure show the fitting by or with or without the quantum correction[31, 34]. The relation between the present theory and the SCR in the quantum critical region will be discussed in the later section.
The intensity map of for and along the high-symmetry line of the Brillouin zone is shown in Fig. 17, where , , and . Since the spin susceptibility is of the RPA-type with the renormalized vertex , the support of the continuum excitations is given by the possible combinations of the particle-hole pair excitations. Note that the critically enhanced collective paramagnon excitation is prominent at lower temperatures.
The imaginary part of the dynamical staggered susceptibility for is shown in Fig 18. It shows the strong dependence to enhance the low-energy excitations as decreases. As will be shown below, the critical AF spin fluctuations result in the pseudogap behavior in the one-particle excitation, which should give a feedback to the low-energy excitation of as well. Since the present theory lacks the self-consistency between the one-particle and the two-particle levels, no suppressions in the low-energy excitations appears as a feedback. The full self-consistent treatment remains an open and very challenging problem.
Let us now turn to the one-particle properties. The dependence of the DOS for is shown in Fig. 19. As decreases the central quasiparticle peak develops, and simultaneously the pseudogap in the quasiparticle peak becomes deeper at . Note that the development of the pseudogap begins at higher temperature where the staggered spin susceptibility deviates from the Curie-Weiss behavior. The critical AF spin fluctuations have little influence on the incoherent part, and the transfer of the weight mainly takes place within the quasiparticle peak. This tendency is in contrast to that of the DMFT with nonlocal corrections where considerable amount of the weight transfers to the incoherent Hubbard satellite peaks from the quasiparticle peak[19, 23, 53].
The spectral intensity map of along the high-symmetry line of the Brillouin zone for at is shown in Fig. 20. Since there is a perfect nesting at half filling and the magnetic Brillouin zone coincides with the Fermi surface, the pseudogap behaviors appear in the coherent quasiparticle bands everywhere on the Fermi surface (see, the X-M’ line). The overall structure of is very similar to that obtained by the DMFT with nonlocal corrections[23, 53].
Now, let us consider the hole doping and the associated superconducting transitions. In the present theory, the even-pair and the odd-pair irreducible vertices are given by
Provided the irreducible vertices, we determine the superconducting transition temperature by the linearized BS equation,
where signals the superconducting instability of the normal state, and its eigenvector provides the symmetry of the gap function. Here, and are given by (49). For the approximate vertices (76), we easily solve the linearized BS equation by using the power method with the fast fourier transformation (FFT) algorithm. On the Hubbard model, we only consider the -type (-type) order parameter, which is the favorable symmetry mediated by the AF spin fluctuations.
Away from the half filling, the “mean-field” decreases and it vanishes at as shown in Fig. 22. Solving the linearized BS equation, the maximum eigenvalue does not reach unity at the lowest temperature (not shown) for the doping range, . As compared with the FLEX, which is often used to discuss the superconducting instability, the present theory uses the much smaller coupling constant to the spin fluctuations, renormalized from . It is the main cause to suppress in the present work. However, the constant approximation for the irreducible vertices may underestimate the coupling constant to the spin fluctuations. If we considered the momentum and/or the frequency dependence, would have a minimum at and it would approach to away from . Since the strongest fluctuation tends to break the Cooper pairs, the electrons make use of a broad range of fluctuations around the maximum point of where the coupling constant would be larger than its averaged value. In order to mimic this effect, we merely replace with (). The simple enhancement of of course needs to be justified by more elaborate calculations with -dependent vertices. Nevertheless, it is useful to grasp the doping dependence of . For instance, setting (), we have the dependence of the eigenvalue as shown in Fig. 21. In this case, the eigenvalues reach unity with the non-monotonous doping dependence. Note that is suppressed strongly when the system goes into the pseudogap region. The superconducting phase boundary is indicated by the closed triangle in Fig. 22. In order to elucidate the retardation effect in the self-energy, we also calculate using the bare propagator, i.e., , in the linearized BS kernel (76) with no enhancement factor, . In contrast to the decreasing behavior of to the pseudogap region, we obtain the monotonously increasing as increases (the open triangle in Fig. 22). Therefore, we conclude that the suppression of is mainly due to the retardation effect of the self-energy.
5 Discussions and Summary
5.1 Relation to the SCR theory in the vicinity of the QCP
We discuss the relation between the SCR and the present theory near the QCP. The core of the SCR theory is the self-consistent equation,
where the dynamical spin susceptibility is approximated as
in the critical region, with being the dynamical exponent. The is a magnetic ordering vector. The bracket means the -average, , and the subscript represents that the significant dependence is included in the average. The coefficient is related to the mode-mode coupling constant in the 4th-order free-energy functional. is a square of the inverse correlation length, and it is regarded as a measure of the distance from the magnetic QCP. The peculiar dependence of determined by (77) governs the critical exponents of various quantities in the QCP region. Setting at , we obtain the condition for to give the QCP,
In the present theory, if we have the single peak in at , we have the relation,
where the weak dependence in can be neglected in the QCP region as long as additional singularities such as the van Hove singularity and the perfect nesting are not involved. The self-consistent equation (68) for in the critical region may read
where is constant contributions from the charge and the PP fluctuations, and is an constant. In the leading order in , it reduces to
which is the same form as (77) with different coefficients. The formal equivalence of the self-consistent equations thus leads to the same exponent of , and other relevant exponents in the QCP region.
Note that in general the basic parquet equations lead to the same critical exponents of a self-consistent expansion for the corresponding order-parameter field theories where is the number of order-parameter components. In contrast to a theory that handles a mode divergence in a single channel, the present theory naturally takes account of interference effects among different modes in the presence of the plural divergences in different channels or multiple characteristic wave vector .
5.2 Possible improvements for the present theory
Let us consider possible improvements for the effective irreducible vertices. The extension of -dependent vertex, is formally straightforward. Namely, the self-consistent equations (68) now become
where the average has been replaced by the average. Although a formal change is very little, it increases a practical computational effort considerably in the evaluation of the vertex correction terms (62) where the simple FFT algorithm cannot be applied.
One may consider that the extension of the full