format: text

text:
# Irreducible Green Functions Method and Many-Particle Interacting Systems on a Lattice ”Rivista del Nuovo Cimento” vol.25, N 1 (2002) pp.1-91

## 1 Introduction

The basic problems of field theory and statistical mechanics are much similar in many aspects, especially, when we use the method of second quantization and Green functions. In both the cases, we are dealing with systems possessing a large number of degrees of freedom (the energy spectrum is practically a continuous one) and with averages of quantum mechanical operators . In quantum field theory, we mostly consider averages over the ground state, while in statistical mechanics, we consider finite temperatures (ensemble averages) as well as ground-state averages. Great advances have been made during the last decades in statistical physics and condensed matter theory through the use of methods of quantum field theory  \- . It was widely recognized that a successful approximation for determining exited states is based on the quasi-particle concept and the Green function method. For example, the study of highly correlated electron systems has attracted much attention recently  \- , especially after discovery of copper oxide superconductors, a new class of heavy fermions , and low-dimensional compounds , . Although much work for strongly correlated systems has been performed during the last years, it is worthy to remind that the investigation of excitations in many-body systems has been one of the most important and interesting subjects for the last few decades.
The quantum field theoretical techniques have been widely applied to statistical treatment of a large number of interacting particles. Many-body calculations are often done for model many-particle systems by using a perturbation expansion. The basic procedure in many-body theory  is to find a suitable unperturbed Hamiltonian and then to take into account a small perturbation operator. This procedure that works well for weakly interacting systems needs a special reformulation for many-body systems with complex spectra and strong interaction. For many practically interesting cases (e.g. in quantum chemistry problems ), the standard schemes of perturbation expansion must be reformulated greatly  \- . Moreover, many-body systems on a lattice have their own specific features and in some important aspects differ greatly from continuous systems.
In this review that is largely pedagogical we are primarily dealing with the spectra of elementary excitations to learn about quasi-particle many-body dynamics of interacting systems on a lattice. Our analysis is based on the equation-of-motion approach, the derivation of the exact representation of the Dyson equation and construction of an approximate scheme of calculations in a self-consistent way. In this review only some topics in the field are discussed. The emphasis is on the methods rather than on a detailed comparison with the experimental results. We attempt to prove that the approach we suggest produces a more advanced physical picture of the problem of the quasi-particle many-body dynamics.
The most characteristic feature of the recent advancement in the basic research on electronic properties of solids is the development of variety of the new classes of materials with unusual properties: high-$`T_c`$ superconductors, heavy fermion compounds, complex oxides, diluted magnetic semiconductors, perovskite manganites, etc. Contrary to simple metals, where the fundamentals are very well known and the electrons can be represented so that they weakly interact with each other, in these materials, the electrons interact strongly, and moreover their spectra are complicated, i.e. have many branches. This gives rise to interesting phenomena  such as magnetism, metal-insulator transition in oxides, heavy fermions, colossal negative magnetoresistance in manganites, etc., but the understanding of what is going on is in many cases only partial.
The subject of the present paper is a microscopic many-body theory of strongly correlated electron models. A principle importance of these studies is concerned with a fundamental problem of electronic solid state theory, namely with the tendency of 3d(4d) electrons in transition metal compounds and 4f(5f) electrons in rare-earth metal compounds and alloys to exhibit both the localized and delocalized (itinerant) behaviour. Interesting electronic and magnetic properties of these substances are intimately related to this dual behaviour of electrons-.
The problem of adequate description of strongly correlated electron systems has been studied intensively during the last decade,, especially in context of the physics of magnetism, heavy fermions and high-$`T_c`$ superconductivity . The understanding of the true nature of electronic states and their quasi-particle dynamics is one of the central topics of the current experimental and theoretical studies in the field. A plenty of experimental and theoretical results show that this many-body quasi-particle dynamics is quite nontrivial. A vast amount of theoretical searches for a suitable description of strongly correlated fermion systems deal with simplified model Hamiltonians. These include, as workable patterns, the single-impurity Anderson model (SIAM) and Hubbard model. In spite of certain drawbacks, these models exhibit the key physical feature: the competition and interplay between kinetic energy (itinerant) and potential energy (localized) effects. A fully consistent theory of quasi-particle dynamics of both the models is believed to be crucially important for a deeper understanding of the true nature of electronic states in the above-mentioned class of materials. In spite of experimental and theoretical achievements, it remains still much to be understood concerning such systems ,.
Recent theoretical investigations of strongly correlated systems have brought forth a significant variety of the approaches to solve these controversial problems. There is an important aspect of the problem under consideration, namely, how to take adequately into account the lattice (quasi-localized) character of charge carriers, contrary to simplified theories of the type of a weakly interacting electron gas. To match such a trend, we need to develop a systematic theory of correlated systems, to describe, from the first principles of the condensed matter theory and statistical mechanics, the physical properties of this class of materials.
In previous papers, we set up the practical technique of the method of the irreducible Green functions (IGF)  -. This IGF method allows one to describe quasi-particle spectra with damping for systems with complex spectra and strong correlation in a very general and natural way. This scheme differs from the traditional methods of decoupling or terminating an infinite chain of the equations and permits one to construct the relevant dynamic solutions in a self-consistent way on the level of the Dyson equation without decoupling the chain of the equations of motion for the double-time temperature Green functions. The essence of our consideration of dynamic properties of many-body system with strong interaction is related closely with the field theoretical approach, and we use the advantage of the Green- function language and the Dyson equation. It is possible to say that our method emphasizes the fundamental and central role of the Dyson equation for the single-particle dynamics of many-body systems at finite temperature. This approach has been suggested as essential for various many-body systems, and we believe that it bears the real physics of interacting many-particle interacting systems , .
It is the purpose of the present paper to introduce the concepts of irreducible Green functions (or irreducible operators) and Generalized Mean Fields (GMF ) in a simple and coherent fashion to assess the validity of quasi-particle description and mean field theory. The irreducible Green function method is a reformulation of the equation-of-motion approach for the double-time thermal GFs, aimed of operating with the correct functional structure of the required solutions. In this sense, it has all advantages and shortcomings of the Green-function method in comparison, say, with the functional integration technique, that, in turn, has also its own advantages and shortcomings. The usefulness of one or another method depends on the problem we are trying to solve. For the calculation of quasi-particle spectra, the Green-function method is the best. The irreducible-Green-function method adds to this statement: ”for the calculation of the quasi-particle spectra with damping” and gives a workable recipe how to do this in a self-consistent way.
The distinction between elastic and inelastic scattering effects is a fundamental one in the physics of many-body systems, and it is also reflected in a number of other ways than in the mean-field and finite lifetimes. The present review attempts to offer a balanced view of quasi-particle interaction effects in terms of division into elastic- and inelastic-scattering characteristics. For this aim, in the present paper, we discuss the background of the IGF approach more thoroughly. To demonstrate the general analysis, we consider here the calculation of quasi-particle spectra and their damping within various types of correlated electron models to extend the applicability of the general formalism and show flexibility and practical usage of the IGF method.

## 2  Varieties of Green Functions

It is appropriate to remind the ideas underlying the Green- function method, and to discuss briefly why they are particularly useful in the study of interacting many-particle systems.
The Green functions of potential theory  were introduced to find the field which is produced by a source distribution (e.g. the electromagnetic field which is produced by current and charge distribution). The Green functions in field theory are the so-called propagators which describe the temporal development of quantized fields, in its particle aspect, as was shown by Schwinger in his seminal works  \- . The idea of the Green function method is contained in the observation that it is not necessary to attempt to calculate all the wave functions and energy levels of a system. Instead, it is more instructive to study the way in which it responds to simple perturbations, for example, by adding or removing particles, or by applying external fields.
There is a variety of Green functions  and there are Green functions for one particle, two particles…, n particles. A considerable progress in studying the spectra of elementary excitations and thermodynamic properties of many-body systems has been for most part due to the development of the temperature-dependent Green-functions methods.

### 2.1 Temperature Green Functions

The temperature dependent Green functions were introduced by Matsubara . He considered a many-particle system with the Hamiltonian

(1)
$$H=H_0+V$$

and observed a remarkable similarity that exists between the evaluation of the grand partition function of the system and the vacuum expectation of the so-called S-matrix in quantum field theory

$$Z=Tr\mathrm{exp}[(\mu NH_0)\beta ]S(\beta );S(\beta )=1_0^\beta V(\tau )S(\tau )𝑑\tau $$

where $`\beta =(kT)^1`$. In essence, Matsubara observed and exploited, to great advantage, formal similarities between the statistical operator $`\mathrm{exp}(\beta H)`$ and the quantum-mechanical time-evolution operator $`\mathrm{exp}(iHt)`$. As a result, he introduced thermal ( temperature-dependent ) Green functions which we call now the Matsubara Green functions.
We note that the thermodynamic perturbation theory has been invented by Peierls . For the free energy of a weakly interacting system he derived the expansion up to second order in perturbation:

$$F=F_0+\underset{n}{}V_{nn}\rho _n+\underset{m,n}{}\frac{|V_{nm}|^2\rho _n}{E_n^0E_m^0}\frac{\beta }{2}\underset{n}{}V_{nn}^2\rho _n+\frac{\beta }{2}\left(\underset{n}{}V_{nn}\rho _n\right)^2$$

where $`\rho _n=\mathrm{exp}[\beta (F_0E_n^0)]`$ and $`\mathrm{exp}(\beta F_0)=_n\mathrm{exp}(\beta E_n^0)`$. By using the expansion of $`S(\beta )`$ up to second order

(2)
$$S(\beta )=1_0^\beta V(\tau )𝑑\tau +_0^\beta 𝑑\tau _1_0^{\tau _1}𝑑\tau _2V(\tau _1)V(\tau _2)+\mathrm{}$$

and rearranging the terms in the expression for $`Z`$, it can be shown that the Peierls result for the thermodynamic potential $`\mathrm{\Omega }`$ can be reproduced by the Matsubara technique (for a canonical ensemble).
Though the use of Green functions is related traditionally with the perturbation theory through the use of diagram techniques, in paper  a prophetic remark has been made:

> ”… it is desirable to avoid founding the formal theory of the Green functions on the restricted basis provided by the assumption of expandability in powers of coupling constants”.

Since the most important aspect of the many-body theory is the necessity of taking properly into account the interaction between particles, that changes ( sometimes drastically ) the behaviour of non-interacting particles, this remark of Schwinger is still extremely actual and important.
Since that time, a great deal of work has been done, and many different variants of the Green functions have been proposed for studies of equilibrium and non-equilibrium properties of many-particle systems. We can mention, in particular, the methods of Martin and Schwinger  and of Kadanoff and Baym . Martin and Schwinger formulated the GF theory not in terms of conventional diagrammatic techniques, but in terms of functional-derivative techniques that reduces the many-body problem directly to the solution of a coupled set of nonlinear integral equations (see also). The approach of Kadanoff and Baym establishes general rules for obtaining approximations which preserve the conservation laws ( sometimes called conserving approximations  ). As many transport coefficients are related to conservation laws, one should take care of it when calculating the two-particle and one-particle Green functions. The random-phase approximation, that is an essential point of the whole Kadanoff-Baym method, does this and so preserves the appropriate conservation laws. It should be noted, however, that the Martin-Schwinger and Kadanoff-Baym methods in their initial form were formulated for treating the continuum models and should be adapted to study lattice models, as well.
However, as was claimed by Matsubara in his subsequent paper , the most convenient way to describe the equilibrium average of any observable or time-dependent response of a system to external disturbances is to express them in terms of a set of the double-time, or Bogoliubov-Tyablikov, Green functions.
The aim of the present paper is to suggest and justify that an approach , the irreducible Green functions (IGF) method , that is in essence a suitable reformulation of an equation-of-motion approach for the double-time temperature-dependent Green functions provides an effective and self-consistent scheme for description the many-body quasi-particle dynamics of strongly interacting many-particle systems with complex spectra. This IGF method provides some systematization of approximations and removes (at least partially) the difficulties usually encountered in the termination of the hierarchy of equations of motion for the GF.

### 2.2 Double-time Green Functions

In this Section, we briefly review the double-time temperature-dependent Green functions .
The double-time temperature-dependent Green functions were introduced by Bogoliubov and Tyablikov  and reviewed by Zubarev  and Tyablikov .
Consider a many-particle system with the time-independent Hamiltonian $`𝐇=H\mu N`$; $`\mu `$ is the chemical potential, $`N`$ is the operator of the total number of particles, and we have chosen our units so that $`\mathrm{}=1`$. Let $`A(t)`$ and $`B(t^{})`$ be some operators . The time development of these operators in the Heisenberg representation is given by:

(3)
$$A(t)=\mathrm{exp}(i𝐇t)A(0)\mathrm{exp}(i𝐇t)$$

We define three types of Green functions, the retarded, advanced, and causal Green functions:

(4) $`G^r=<<A(t),B(t^{})>>^r=i\theta (tt^{})<[A(t)B(t^{})]_\eta >,\eta =\pm 1.`$
(5) $`G^a=<<A(t),B(t^{})>>^a=i\theta (t^{}t)<[A(t)B(t^{})]_\eta >,\eta =\pm 1.`$
(6) $`G^c=<<A(t),B(t^{})>>^c=iT<A(t)B(t^{})>=`$
$`i\theta (tt^{})<A(t)B(t^{})>+\eta i\theta (t^{}t)<B(t^{})A(t)>,\eta =\pm 1.`$
where $`<\mathrm{}>`$ is the average over a grand canonical ensemble, $`\theta (t)`$ is a step function, and square brackets represent the commutator or anticommutator

(7)
$$[A,B]_\eta =AB\eta BA$$

Differentiating a Green function with respect to one of the arguments, for example, the first argument, we can obtain the equation (equation-of-motion) describing the development of this function with time

(8)
$$id/dtG^\alpha (t,t^{})=\delta (tt^{})<[A,B]_\eta >+<<[A,H](t),B(t^{})>>^\alpha ;\alpha =r,a,c$$

Since this differential equation contains an inhomogeneous term with $`\delta `$-type factors, we are dealing formally with the equation similar to the usual one for the Green function  and for this reason, we use the term the Green function. We note that the equation of motion is of the same functional form for all the three types of Green functions ( i.e. retarded, advanced, and causal ). However, the boundary conditions for $`t`$ are different for the retarded, advanced, and causal functions .
The next differentiation gives an infinite chain of coupled equations of motion

(9) $`(i)^nd^n/dt^nG(t,t^{})=`$
$`{\displaystyle \underset{k=1}{\overset{n}{}}}(i)^{nk}d^{nk}/dt^{nk}\delta (tt^{})<[[\mathrm{}[A,\underset{k1}{\underset{}{H]\mathrm{}H]}},B]_\eta >`$
$`+<<[[\mathrm{}[A,\underset{n}{\underset{}{H]\mathrm{}H]}}(t),B(t^{})>>`$
To solve the differential equation-of-motion, we should consider the Fourier time transforms of the Green functions:

(10)
$$G_{AB}(tt^{})=(2\pi )^1_{\mathrm{}}^{\mathrm{}}𝑑\omega G_{AB}(\omega )\mathrm{exp}[i\omega (tt^{})],$$

(11)
$$G_{AB}(\omega )=<<A|B>>_\omega =_{\mathrm{}}^{\mathrm{}}𝑑tG_{AB}(t)\mathrm{exp}(i\omega t),$$

By inserting (10) into (8) and (9), we obtain

(12) $`\omega G_{AB}(\omega )=<[A,B]_\eta >+<<[A,H]|B>>_\omega ;`$
(13) $`\omega ^nG_{AB}(\omega )={\displaystyle \underset{k=1}{\overset{n}{}}}\omega ^{nk}<[[\mathrm{}[A,\underset{k1}{\underset{}{H]\mathrm{}H]}},B]_\eta >`$
$`+<<[[\mathrm{}[A,\underset{n}{\underset{}{H]\mathrm{}H}}]|B>>_\omega `$
It is often convenient to differentiate of the Green function with respect to the second time $`t^{}`$ . In terms of Fourier time transforms, the corresponding equations of motion read

(14) $`\omega G_{AB}(\omega )=<[A,B]_\eta >+<<A|[B,H]>>_\omega ;`$
$`(1)^n\omega ^nG_{AB}(\omega )={\displaystyle \underset{k=1}{\overset{n}{}}}(1)^{nk}\omega ^{nk}<[A,[\mathrm{}[B,\underset{k1}{\underset{}{H]\mathrm{}H]}}]_\eta >`$
(15) $`+<<A|[\mathrm{}[B,\underset{n}{\underset{}{H]\mathrm{}H}}]>>_\omega `$
It is rather difficult problem to solve the infinite chain of coupled equations of motion (2.2) and (15). It is well established now that the usefulness of the retarded and advanced Green functions is deeply related with the dispersion relations , that provide the boundary conditions in the form of spectral representations of the Green functions.

### 2.3 Spectral Representations

The GFs are linear combinations of the time correlation functions:

(16) $`F_{AB}(tt^{})=<A(t)B(t^{})>={\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\omega \mathrm{exp}[i\omega (tt^{})]A_{AB}(\omega )`$
(17) $`F_{BA}(t^{}t)=<B(t^{})A(t)>={\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\omega \mathrm{exp}[i\omega (t^{}t)]A_{BA}(\omega )`$
Here, the Fourier transforms $`A_{AB}(\omega )`$ and $`A_{BA}(\omega )`$ are of the form

(18) $`A_{BA}(\omega )=`$
$`Q^12\pi {\displaystyle \underset{m,n}{}}\mathrm{exp}(\beta E_n)(\psi _n^{}B\psi _m)(\psi _m^{}A\psi _n)\delta (E_nE_m\omega )`$
(19) $`A_{AB}=\mathrm{exp}(\beta \omega )A_{BA}(\omega )`$
The expressions (18) and (19) are spectral representations of the corresponding time correlation functions. The quantities $`A_{AB}`$ and $`A_{BA}`$ are spectral densities or spectral weight functions.
It is convenient to define

(20) $`F_{BA}(0)=<B(t)A(t)>={\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\omega A(\omega )`$
(21) $`F_{AB}(0)=<A(t)B(t)>={\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\omega \mathrm{exp}(\beta \omega )A(\omega )`$
Then, the spectral representations of the Green functions can be expressed in the form

(22) $`G^r(\omega )=<<A|B>>_\omega ^r=`$
$`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\omega ^{}}{\omega \omega ^{}+iϵ}}[\mathrm{exp}(\beta \omega ^{})\eta ]A(\omega ^{})`$
(23) $`G^a(\omega )=<<A|B>>_\omega ^a=`$
$`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\omega ^{}}{\omega \omega ^{}iϵ}}[\mathrm{exp}(\beta \omega ^{})\eta ]A(\omega ^{})`$
The most important practical consequence of spectral representations for the retarded and advanced GFs is the so-called spectral theorem. The spectral theorem can be written as

(24) $`<B(t^{})A(t)>=`$
$`{\displaystyle \frac{1}{\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\omega \mathrm{exp}[i\omega (tt^{})][\mathrm{exp}(\beta \omega )\eta ]^1ImG_{AB}(\omega +iϵ)`$
(25) $`<A(t)B(t^{})>=`$
$`{\displaystyle \frac{1}{\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\omega \mathrm{exp}(\beta \omega )\mathrm{exp}[i\omega (tt^{})][\mathrm{exp}(\beta \omega )\eta ]^1ImG_{AB}(\omega +iϵ)`$
Expressions (24) and (25) are of fundamental importance. They directly relate the statistical averages with the Fourier transforms of the corresponding GFs. The problem of evaluating the latter is thus reduced to finding their Fourier transforms , providing the practical usefulness of the Green functions technique .

## 3 Irreducible Green Functions Method

In this Section, we discuss the main ideas of the IGF approach which allows one to describe completely quasi-particle spectra with damping in a very natural way.
We reformulated the two-time GF method ,  to the form, which is especially adjusted ,  for correlated fermion systems on a lattice and systems with complex spectra ,. A similar method was proposed in paper  for Bose systems ( anharmonic phonons and spin dynamics of Heisenberg ferromagnet ). The very important concept of the whole method is the Generalized Mean Field (GMF), as it was formulated in ref. . These GMFs have a complicated structure for the strongly correlated case and complex spectra and are not reduced to the functional of mean densities of the electrons or spins when one calculates excitation spectra at finite temperatures.

### 3.1  Outline of IGF Method

To clarify the foregoing, let us consider a retarded GF of the form

(26)
$$G^r=<<A(t),A^{}(t^{})>>=i\theta (tt^{})<[A(t)A^{}(t^{})]_\eta >,\eta =\pm 1$$

As an introduction to the concept of IGFs, let us describe the main ideas of this approach in a symbolic and simplified form. To calculate the retarded GF $`G(tt^{})`$, let us write down the equation of motion for it:

(27)
$$\omega G(\omega )=<[A,A^{}]_\eta >+<<[A,H]_{}A^{}>>_\omega $$

The essence of the method is as follows :
It is based on the notion of the ”IRREDUCIBLE” parts of GFs (or the irreducible parts of the operators, $`A`$ and $`A^{}`$, out of which the GF is constructed) in terms of which it is possible, without recourse to a truncation of the hierarchy of equations for the GFs, to write down the exact Dyson equation and to obtain an exact analytic representation for the self-energy operator. By definition, we introduce the irreducible part (ir) of the GF

(28)
$$^{(ir)}<<[A,H]_{}|A^{}>>=<<[A,H]_{}zA|A^{}>>$$

The unknown constant z is defined by the condition (or constraint)

(29)
$$<[[A,H]_{}^{(ir)},A^{}]_\eta >=0$$

which is an analogue of the orthogonality condition in the Mori formalism ( see ref.). From the condition (29) one can find:

(30)
$$z=\frac{<[[A,H]_{},A^{}]_\eta >}{<[A,A^{}]_\eta >}=\frac{M_1}{M_0}$$

Here $`M_0`$ and $`M_1`$ are the zeroth and first order moments of the spectral density. Therefore, the irreducible GFs are defined so that they cannot be reduced to the lower-order ones by any kind of decoupling. It is worth noting that the term ”irreducible” in a group theory means a representation of a symmetry operation that cannot be expressed in terms of lower dimensional representations. Irreducible (or connected ) correlation functions are known in statistical mechanics (cf.). In the diagrammatic approach, the irreducible vertices are defined as graphs that do not contain inner parts connected by the $`G^0`$-line. With the aid of the definition (30) these concepts are translated into the language of retarded and advanced GFs. This procedure extracts all relevant (for the problem under consideration) mean-field contributions and puts them into the generalized mean-field GF which is defined here as

(31)
$$G^0(\omega )=\frac{<[A,A^{}]_\eta >}{(\omega z)}$$

To calculate the IGF $`{}_{}{}^{(ir)}<<[A,H]_{}(t),A^{}(t^{})>>`$ in (27), we have to write the equation of motion for it after differentiation with respect to the second time variable $`t^{}`$. It should be noted that the trick of two-time differentiation with respect to the first time $`t`$ and second time $`t^{}`$ (in one equation of motion) was introduced for the first time by Tserkovnikov .
The condition of orthogonality (31) removes the inhomogeneous term from this equation and is a very crucial point of the whole approach. If one introduces the irreducible part for the right-hand side operator as discussed above for the “left” operator, the equation of motion (27) can be exactly rewritten in the following form

(32)
$$G=G^0+G^0PG^0$$

The scattering operator $`P`$ is given by

(33)
$$P=(M_0)^1(^{(ir)}<<[A,H]_{}|[A^{},H]_{}>>^{(ir)})(M_0)^1$$

The structure of equation (34) enables us to determine the self-energy operator $`M`$, by analogy with the diagram technique

(34)
$$P=M+MG^0P$$

From the definition (34) it follows that the self-energy operator $`M`$ is defined as a proper (in the diagrammatic language, “connected”) part of the scattering operator $`M=(P)^p`$. As a result, we obtain the exact Dyson equation for the thermodynamic double-time Green functions:

(35)
$$G=G^0+G^0MG$$

The difference between $`P`$ and $`M`$ can be regarded as two different solutions of two integral equations (32) and (35). But from the Dyson equation (35) only the full GF is seen to be expressed as a formal solution of the form

(36)
$$G=[(G^0)^1M]^1$$

Equation (38) can be regarded as an alternative form of the Dyson equation (35) and the definition of $`M`$ provided that the generalized mean-field GF $`G^0`$ is specified. On the contrary , for the scattering operator $`P`$, instead of property $`G^0G^1+G^0M=1`$, one has the property

$$(G^0)^1G^1=PG^0G^1$$

Thus, the very functional form of the formal solution (38) determines the difference between $`P`$ and $`M`$ precisely.
Thus, by introducing irreducible parts of GF (or irreducible parts of the operators, out of which the GF is constructed) the equation of motion (27) for the GF can exactly be ( but using orthogonality constraint (29)) transformed into the Dyson equation for the double-time thermal GF (35). This result is very remarkable , because the traditional form of the GF method does not include this point. Notice that all quantities thus considered are exact. Approximations can be generated not by truncating the set of coupled equations of motions but by a specific approximation of the functional form of the mass operator $`M`$ within a self-consistent scheme, expressing $`M`$ in terms of initial GF

$$MF[G]$$

Different approximations are relevant to different physical situations.
The projection operator technique  has essentially the same philosophy, but with using the constraint (29) in our approach we emphasize the fundamental and central role of the Dyson equation for the calculation of single-particle properties of many-body systems. The problem of reducing the whole hierarchy of equations involving higher-order GFs by a coupled nonlinear set of integro-differential equations connecting the single-particle GF to the self-energy operator is rather nontrivial ( cf.). A characteristic feature of these equations is that, besides the single-particle GF, they involve also higher-order GF. The irreducible counterparts of the GFs, vertex functions, etc, serve to identify correctly the self-energy as

$$M=G_0^1G^1$$

The integral form of Dyson equation (35) gives $`M`$ the physical meaning of a nonlocal and energy-dependent effective single-particle potential. This meaning can be verified for the exact self-energy through the diagrammatic expansion for the causal GF.
It is important to note that for the retarded and advanced GFs, the notion of the proper part $`M=(P)^p`$ is symbolic in nature . In a certain sense, it is possible to say that it is defined here by analogy with the irreducible many-particle $`T`$-matrix. Furthermore, by analogy with the diagrammatic technique, we can also introduce the proper part defined as a solution to the integral equation (34). These analogues allow us to understand better the formal structure of the Dyson equation for the double-time thermal GF but only in a symbolic form . However, because of the identical form of the equations for GFs for all three types ( advanced, retarded, and causal ), we can convert in each stage of calculations to causal GF and, thereby, confirm the substantiated nature of definition (34)! We therefore should speak of an analogy of the Dyson equation. Hereafter, we drop this stipulating, since it does not cause any misunderstanding. In a sense, the IGF method is a variant of the Gram-Schmidt orthogonalization procedure (see Appendix A ).
It should be emphasized that the scheme presented above gives just a general idea of the IGF method. A more exact explanation why one should not introduce the approximation already in $`P`$, instead of having to work out $`M`$, is given below when working out the application of the method to specific problems.
The general philosophy of the IGF method is in the separation and identification of elastic scattering effects and inelastic ones. This latter point is quite often underestimated, and both effects are mixed. However, as far as the right definition of quasi-particle damping is concerned, the separation of elastic and inelastic scattering processes is believed to be crucially important for many-body systems with complicated spectra and strong interaction.
From a technical point of view, the elastic GMF renormalizations can exhibit quite a nontrivial structure. To obtain this structure correctly, one should construct the full GF from the complete algebra of relevant operators and develop a special projection procedure for higher-order GFs in accordance with a given algebra. Then the natural question arises how to select the relevant set of operators $`\{A_1,A_2,\mathrm{}A_n\}`$ , describing the ”relevant degrees of freedom”. The above consideration suggests an intuitive and heuristic way to the suitable procedure as arising from an infinite chain of equations of motion (12). Let us consider the column

$$\left(\begin{array}{c}A_1\\ A_2\\ \mathrm{}\\ A_n\end{array}\right)$$

where

$$A_1=A,A_2=[A,H],A_3=[[A,H],H],\mathrm{}A_n=[[\mathrm{}[A,\underset{n}{\underset{}{H]\mathrm{}H}}]$$

Then the most general possible Green function can be expressed as a matrix

$$\widehat{G}=<<\left(\begin{array}{c}A_1\\ A_2\\ \mathrm{}\\ A_n\end{array}\right)|\left(\begin{array}{cccc}A_1^{}& A_2^{}& \mathrm{}& A_n^{}\end{array}\right)>>$$

This generalized Green function describes the one-, two- and $`n`$-particle dynamics. The equation of motion for it includes, as a particular case, the Dyson equation for single-particle Green function, the Bethe-Salpeter equation, which is the equation of motion for the two-particle Green function and which is an analogue of the Dyson equation, etc. The corresponding reduced equations should be extracted from the equation of motion for the generalized GF with the aid of the special techniques such as the projection method and similar techniques. This must be a final goal towards a real understanding of the true many-body dynamics. At this point, it is worthwhile to underline that the above discussion is a heuristic scheme only but not a straightforward recipe. The specific method of introducing the IGFs depends on the form of operators $`A_n`$, the type of the Hamiltonian, and conditions of the problem. The irreducible parts in higher-order equations and connection with Mori formalism was considered by Tserkovnikov . The incorporation of irreducible parts in higher-order equations and connection with the moment expansion was studied in ref.  ( see Appendix B ).
Here a sketchy form of the IGF method is presented. The aim to introduce the general scheme and to lay the groundwork for generalizations and specific applications is expounded in the next Sections. We demonstrate below that the IGF method is a powerful tool for describing the quasi-particle excitation spectra, allowing a deeper understanding of elastic and inelastic quasi-particle scattering effects and corresponding aspects of damping and finite lifetimes. In the present context, it provides a clear link between the equation-of-motion approach and the diagrammatic methods due to derivation of the Dyson equation (37). Moreover, due to the fact that it allows the approximate treatment of the self-energy effects on a final stage, it yields a systematic way of the construction of approximate solutions.
It is necessary to emphasize that there is an intimate connection between an adequate introduction of mean fields and internal symmetries of the Hamiltonian. To test these ideas further, in the following Sections, we analyze the mean field and generalized mean field concepts for various many-body systems on a lattice.

## 4 Many-Particle Interacting Systems on a Lattice

### 4.1 Spin Systems on a Lattice

There exists a big variety of magnetic materials. The group of magnetic insulators is of a special importance. For the group of systems considered in this Section, the physical picture can be represented by a model in which the localized magnetic moments originating from ions with incomplete shells interact through a short-range interaction. Individual spin moments form a regular lattice. The first model of a lattice spin system was constructed to describe a linear chain of projected electron spins with nearest-neighbor coupling. This was the famous Lenz-Izing model which was thought to yield a more sophisticated description of ferromagnetism than the Weiss uniform molecular field picture. However, in this model, only one spin component is significant. As a result, the system has no collective dynamics. The quantum states that are eigenstates of the relevant spin components are stationary states. The collective dynamics of magnetic systems is of great importance since it is related to the study of low-lying excitations and their interactions. This is the main aim of the present consideration. Although the Izing model was an intuitively right step forward from the uniform Weiss molecular field picture, the physical meaning of the model coupling constant remained completely unclear. The concept of the exchange coupling of spins of two or more nonsinglet atoms appeared as a result of the Heitler-London consideration of chemical bond. This theory and the Dirac analysis of the singlet-triplet splitting in the helium spectrum stimulated Heisenberg to make a next essential step. Heisenberg suggested that the exchange interaction could be the relevant mechanism responsible for ferromagnetism.

#### 4.1.1  Heisenberg Ferromagnet

The Heisenberg model of a system of spins on various lattices ( which was actually written down explicitly by van Vleck ) is termed the Heisenberg ferromagnet and establishes the origin of the coupling constant as the exchange energy. The Heisenberg ferromagnet in a magnetic field $`H`$ is described by the Hamiltonian

(37)
$$H=\underset{ij}{}J(ij)\stackrel{}{S}_i\stackrel{}{S}_jg\mu _BH\underset{i}{}S_i^z$$

The coupling coefficient $`J(ij)`$ is the measure of the exchange interaction between spins at the lattice sites $`i`$ and $`j`$ and is defined usually to have the property J(i - j = 0) = 0. This constraint means that only the inter-exchange interactions are taken into account. However, in some complicated magnetic salts, it is necessary to consider an ”effective” intra-site (see) interaction (Hund-rule-type terms). The coupling, in principle, can be of a more general type (non-Heisenberg terms). These aspects of construction of a more general Hamiltonian are very interesting, but we do not pause here to give the details.
For crystal lattices in which every ion is at the centre of symmetry, the exchange parameter has the property

$$J(ij)=J(ji)$$

We can rewrite then the Hamiltonian (37) as

(38)
$$H=\underset{ij}{}J(ij)(S_i^zS_j^z+S_i^+S_j^{})$$

Here $`S^\pm =S^x\pm iS^y`$ are the raising and lowering spin angular momentum operators. The complete set of spin commutation relations is

$`[S_i^+,S_j^{}]_{}=2S_i^z\delta _{ij};[S_i^+,S_i^{}]_+=2S(S+1)2(S_i^z)^2;`$
$`[S_i^{},S_j^z]_{}=\pm S_i^{}\delta _{ij};S_i^z=S(S+1)(S_i^z)^2S_i^{}S_i^+;`$
$`(S_i^+)^{2S+1}=0,(S_i^{})^{2S+1}=0`$
We omit the term of interaction of the spin with an external magnetic field for the brevity of notation. The statistical mechanical problem involving this Hamiltonian was not exactly solved, but many approximate solutions were obtained.
To proceed further, it is important to note that for the isotropic Heisenberg model, the total $`z`$-component of spin $`S_{tot}^z=_iS_i^z`$ is a constant of motion, i.e.

$$[H,S_{tot}^z]=0$$

There are cases when the total spin is not a constant of motion, as, for instance, for the Heisenberg model with the dipole terms added.
Let us define the eigenstate $`|\psi _0>`$ so that $`S_i^+|\psi _0>=0`$ for all lattice sites $`R_i`$. It is clear that $`|\psi _0>`$ is a state in which all the spins are fully aligned and for which $`S_i^z|\psi _0>=S|\psi _0>`$. We also have

$$J_\stackrel{}{k}=\underset{i}{}e^{(i\stackrel{}{k}\stackrel{}{R}_i)}J(i)=J_\stackrel{}{k}$$

, where the reciprocal vectors $`\stackrel{}{k}`$ are defined by cyclic boundary conditions. Then we obtain

$$H|\psi _0>=\underset{ij}{}J(ij)S^2=NS^2J_0$$

Here $`N`$ is the total number of ions in the crystal. So, for the isotropic Heisenberg ferromagnet, the ground state $`|\psi _0>`$ has an energy $`NS^2J_0`$.
The state $`|\psi _0>`$ corresponds to a total spin $`NS`$.
Let us consider now the first excited state. This state can be constructed by creating one unit of spin deviation in the system. As a result, the total spin is $`NS1`$. The state

$$|\psi _k>=\frac{1}{\sqrt{(2SN)}}\underset{j}{}e^{(i\stackrel{}{k}\stackrel{}{R}_j)}S_j^{}|\psi _0>$$

is an eigenstate of $`H`$ which corresponds to a single magnon of the energy

(39)
$$\omega _0^{(fm)}(k)=2S(J_0J_k)$$

Note that the role of translational symmetry, i.e. the regular lattice of spins, is essential, since the state $`|\psi _k>`$ is constructed from the fully aligned state by decreasing the spin at each site and summing over all spins with the phase factor $`e^{i\stackrel{}{k}\stackrel{}{R}_j}`$. It is easy to verify that

$$<\psi _k|S_{tot}^z|\psi _k>=NS1$$

The above consideration was possible because we knew the exact ground state of the Hamiltonian . There are many models where this is not the case. For example, we do not know the exact ground state of a Heisenberg ferromagnet with dipolar forces and the ground state of the Heisenberg antiferromagnet.

#### 4.1.2 Heisenberg Antiferromagnet

We now discuss the Heisenberg model of the antiferromagnet which is more complicated to analyse. The fundamental problem here is that the exact ground state is unknown. We consider, for simplicity, a two-sublattice structure in which nearest neighbour ions on opposite sublattices interact through the Heisenberg exchange. For a system of ions on two sublattices, the Hamiltonian is

(40)
$$H=J\underset{m,\delta }{}\stackrel{}{S}_m\stackrel{}{S}_{m+\delta }+J\underset{n,\delta }{}\stackrel{}{S}_n\stackrel{}{S}_{n+\delta }$$

Here the notation $`m=\stackrel{}{R}_m`$ means the position vectors of ions on one sublattice ($`a`$) and $`n`$ for the ions on the other ($`b`$). Nearest neighbor ions on different sublattices are a distance $`|\stackrel{}{\delta }|`$ apart. ( The anisotropy field $`\mu H_A(_mS_m^z_nS_n^z)`$, which is not written down explicitly, is taken to be parallel to the z axis. ) The simplest crystal structures that can be constructed from two interpenetrating identical sublattices are the body-centered and simple cubic.
The exact ground state of this Hamiltonian is not known. One can use the approximation of taking the ground state to be a classical ground state, usually called the Neel state, in which the spins of the ions on each sublattice are oppositely aligned along the z axis. However, this state is not even an eigenstate of the Hamiltonian (40). Let us remark that the total $`z`$-component of the spin commutes with the Hamiltonian (40). It would be instructive to consider here the construction of a spin wave theory for the low-lying excitations of the Heisenberg antiferromagnet in a sketchy form to clarify the foregoing.
To demonstrate the specifics of Heisenberg antiferromagnet more explicitly, it is convenient to rotate the axes of one sublattice through $`\pi `$ about the x-axis. This transformation preserves the spin operator commutation relations and therefore is canonical. Let us perform the transformation on the $`\stackrel{}{R}_n`$, or $`b`$-sublattice

$$S_n^z\stackrel{~}{S}_n^z;S_n^\pm \stackrel{~}{S}_n^{}$$

The operators $`S_m^\alpha `$ and $`\stackrel{~}{S}_n^\beta `$ commute, because they refer to different sublattices.
The transformation to the momentum representation is modified in comparison with the ferromagnet case

$$S_m^\pm =\frac{1}{N}\underset{\stackrel{}{q}}{}e^{(\pm i\stackrel{}{q}\stackrel{}{R}_m)}S_q^\pm ;\stackrel{~}{S}_m^\pm =\frac{1}{N}\underset{\stackrel{}{q}}{}e^{(i\stackrel{}{q}\stackrel{}{R}_m)}\stackrel{~}{S}_q^\pm $$

Here $`\stackrel{}{q}`$ is the reciprocal lattice vectors for one sublattice, each sublattice containing $`N`$ ions. After these transformations, the Hamiltonian (40) can be rewritten as

(41)
$$H=\frac{1}{2SN}\underset{q}{}2zJS[(S_q^{}S_q^++\stackrel{~}{S}_q^{}\stackrel{~}{S}_q^+)+\gamma _q(S_q^+\stackrel{~}{S}_q^++S_q^{}\stackrel{~}{S}_q^{})]$$

In (41), $`\gamma _q`$ is defined as $`z\gamma _q=_{m=n.n.}exp(i\stackrel{}{q}\stackrel{}{R}_m)`$, and $`z`$ is the number of nearest neighbors; the constant terms and the products of four operators are omitted. Thus the Hamiltonian of the Heisenberg antifferomagnet is more complicated than that for the ferromagnet. Because it contains two types of spin operators that are coupled together, the diagonalization of (41) has its own specificity.
To diagonalize (41), let us make a linear transformation to new operators ( Bogoliubov transformation )

(42)
$$S_q^+=u_qa_q+v_qb_q^{};\stackrel{~}{S}_q^{}=u_qb_q^{}+v_qa_q$$

with

$$[a_q,a_q^{}^{}]=\delta _{q,q^{}};[b_q,b_q^{}^{}]=\delta _{q,q^{}}$$

The transformation coefficients $`u_k`$ and $`v_k`$ are purely real. To preserve the commutation rules for the spin operators

$$[S_k^+,S_k^{}^{}]=2SN\delta _{k,k^{}}$$

, they should satisfy $`u^2(k)v^2(k)=2SN`$. The transformations from the operators $`(S_q^+,\stackrel{~}{S}_q^{})`$ to the operators $`(a_q,b_q^{})`$ give

$`[(S_q^{}S_q^++\stackrel{~}{S}_q^{}\stackrel{~}{S}_q^+)+\gamma _q(S_q^+\stackrel{~}{S}_q^++S_q^{}\stackrel{~}{S}_q^{})]=`$
$`(a_q^{}a_q+b_q^{}b_q)[(u^2(q)+v^2(q))+2u_qv_q\gamma _q]`$
$`+(a_qb_q+a_q^{}b_q^{})[(u^2(q)+v^2(q))\gamma _q+2u_qv_q]`$
(43) $`+2u_qv_q\gamma _q+2v^2(q)`$
We represented Hamiltonian (41) as a form quadratic in the Bose operators $`(a_q,b_q^{})`$ . We shall now consider the problem of diagonalization of this form. To diagonalize (41), we should require that

$$2u_qv_q+(u^2(q)+v^2(q))\gamma _q=0$$

Then we obtain

(44)
$$2u^2(q)=2SN\frac{(1+\kappa _q)}{\kappa _q};2v^2(q)=2SN\frac{(1\kappa _q)}{\kappa _q}$$

Here the following notation was introduced: $`\kappa _q=\sqrt{(1\gamma _q^2)}`$ and $`2u_qv_q=2SN\gamma _q/\kappa _q`$ After the transformation (42), we get, instead of (41),

(45)
$$H=\underset{k}{}\omega _0^{(afm)}(k)(a_q^{}a_q+b_q^{}b_q)$$

with

(46)
$$\omega _0^{(afm)}(k)=2zJS\sqrt{1\gamma _k^2}$$

Expression (45) contains two terms, each with the same energy spectrum. Thus, there are two degenerate spin wave modes, because there can be two kinds of precession of the spin about the anisotropy direction. The degeneracy is lifted by the application of an external magnetic field in the $`z`$ direction, because in this case the two sublattices become nonequivalent. These results should be kept in mind when discussing the quasi-particle many-body dynamics of the spin lattice models.

### 4.2 Correlated Electrons on a Lattice

The importance of intra-atomic correlation effects in determining the magnetic properties of transition metals and their compounds and oxides was recognized many years ago. The essential basis of studies of metallic magnetism, namely, that the dominant physical mechanism responsible for the observed magnetic properties of the transition metals and their compounds and alloys is the strong intra-atomic correlation in an otherwise tight-binding picture, is generally accepted as being most suitable. The problem of the adequate description of strongly correlated electron systems on a lattice was studied intensively during the last decade, especially in the context of metallic magnetism, heavy fermions, and high-$`T_c`$ superconductivity . The understanding of the true nature of electronic states and their quasi-particle dynamics is one of the central topics of the current experimental and theoretical efforts in the field. The source of spin magnetism in solids is, of course, the Pauli exclusion principle as manifested in the exchange interaction and higher order mechanism. Of particular interest is the fact that the Hartree-Fock or mean field theory, i.e. the theory including exchange but not correlation effects, invariably overestimates the tendency to magnetism. This fact obviously complicated the already complicated problem of magnetism in a metal with the $`d`$ band electrons which, as was mentioned above, are really neither ”local” nor ”itinerant” in a full sense.
The strongly correlated electron systems are systems in which electron correlations dominate. The theoretical studies of strongly correlated systems had as a consequence the formulation of two model Hamiltonians which play a central role in our attempts to get an insight into this complicated problem. These are the Anderson single-impurity model (SIAM)  and Hubbard model . It was only relatively recently recognized that both the models have a very complicated many-body dynamics, and their ”simplicity” manifests itself in the dynamics of two-particle scattering, as was shown via elegant Bethe-anzatz solutions.

#### 4.2.1  Hubbard Model

The model Hamiltonian usually referred to as the Hubbard Hamiltonian,

(47)
$$H=\underset{ij\sigma }{}t_{ij}a_{i\sigma }^{}a_{j\sigma }+U/2\underset{i\sigma }{}n_{i\sigma }n_{i\sigma }$$

includes the intra-atomic Coulomb repulsion $`U`$ and the one-electron hopping energy $`t_{ij}`$. The electron correlation forces electrons to localize in the atomic orbitals which are modelled here by a complete and orthogonal set of the Wannier wave functions $`[\varphi (\stackrel{}{r}\stackrel{}{R}_j)]`$. On the other hand, the kinetic energy is reduced when electrons are delocalized. The main difficulty in solving the Hubbard model correctly is the necessity of taking into account both these effects simultaneously. Thus, the Hamiltonian (47) is specified by two parameters: $`U`$ and the effective electron bandwidth

$$\mathrm{\Delta }=(N^1\underset{ij}{}|t_{ij}|^2)^{1/2}.$$

The band energy of Bloch electrons $`ϵ(\stackrel{}{k})`$ is defined as follows

$$t_{ij}=N^1\underset{\stackrel{}{k}}{}ϵ(\stackrel{}{k})\mathrm{exp}[i\stackrel{}{k}(\stackrel{}{R}_i\stackrel{}{R}_j],$$

where $`N`$ is the number of lattice sites. It is convenient to count the energy from the center of gravity of the band, i.e. $`t_{ii}=t_0=_kϵ(k)=0`$ ( sometimes it is useful to retain $`t_0`$ explicitly ).
This conceptually simple model is mathematically very complicated. The effective electron bandwidth $`\mathrm{\Delta }`$ and Coulomb intra-site integral $`U`$ determine different regimes in 3 dimensions depending on the parameter $`\gamma =\mathrm{\Delta }/U`$. In addition, the Pauli exclusion principle that does not allow two electrons of common spin to be at the same site, i.e. $`n_{i\sigma }^2=n_{i\sigma }`$, plays a crucial role, and it should be taking into account properly while making any approximations. It is usually rather a difficult task to find an interpolating solution for dynamic properties of the Hubbard model for various mean particle densities. To solve this problem with a reasonably accuracy and to describe correctly an interpolated solution from the “band” limit ($`\gamma 1`$) to the “atomic” limit ($`\gamma 0`$), one needs a more sophisticated approach than usual procedures developed for description of the interacting electron gas problem. We have evidently to improve the early Hubbard theory taking into account of variety of possible regimes for the model depending on the electron density, temperature, and values of $`\gamma `$. The single-electron GF

(48)
$$G_{ij\sigma }(\omega )=<<a_{i\sigma }|a_{j\sigma }^{}>>=N^1\underset{\stackrel{}{k}}{}G_\sigma (\stackrel{}{k},\omega )\mathrm{exp}[i\stackrel{}{k}(\stackrel{}{R}_i\stackrel{}{R}_j)],$$

calculated by Hubbard , , has the characteristic two-pole functional structure

(49)
$$G_\sigma (k,\omega )=[F_\sigma (\omega )ϵ(k)]^1$$

where

(50)
$$F_\sigma ^1(\omega )=\frac{\omega (n_\sigma ^+E_{}+n_\sigma ^{}E_+)\lambda }{(\omega E_+n_\sigma ^{}\lambda )(\omega E_{}n_\sigma ^+\lambda )n_\sigma ^+n_\sigma ^{}\lambda ^2}$$

Here $`n^+=n`$ , $`n^{}=1n`$; $`E_+=U`$, $`E_{}=0`$, and $`\lambda `$ is a certain function which depends on parameters of the Hamiltonian. In this approximation, Hubbard took account of the scattering effect of electrons with spins $`\sigma `$ by electrons with spin $`\sigma `$ which are frozen as well as the ”resonance broadening” effect due to the motion of the electrons with spin $`\sigma `$. The ”Hubbard III” decoupling procedure suffered of serious limitations. However, in spite of the limitations, this solution gave the first clue to the qualitative understanding of the property of narrow-band system like the metal-insulator transition.
If $`\lambda `$ is small ($`\lambda 0`$), then expression (50) takes the form:

$$F_\sigma ^1(\omega )\frac{n_\sigma ^{}}{\omega E_{}n_\sigma ^+\lambda }+\frac{n_\sigma ^+}{\omega E_+n_\sigma ^{}\lambda },$$

which corresponds to two shifted subbands with the gap

$$\omega _1\omega _2=(E_+E_{})+(n_\sigma ^{}n_\sigma ^+)\lambda =U+\lambda 2n_\sigma ^+.$$

If $`\lambda `$ is very big, then we obtain

$$F_\sigma ^1(\omega )\frac{\lambda }{[(\omega E_{})n_\sigma ^{}+(\omega E_+)n_\sigma ^+]\lambda }=\frac{1}{\omega (n_\sigma ^+E_+n_\sigma ^{}E_{})}.$$

The latter solution corresponds to a single band centered at the energy $`\omega n_\sigma ^+U`$. Thus, this solution explains qualitatively the appearance of a gap in the density of states when the value of the intra-atomic correlation exceeds a certain critical value, as it was first conjectured by N. Mott.
The two-pole functional structure of the single-particle GF is easy to understand within the formalism that describes the motion of electrons in binary alloys ,. If one introduces the two types of the scattering potentials $`t_\pm (\omega E_\pm )^1`$, then the two kinds of the t-matrix $`T_+`$ and $`T_{}`$ appears which satisfy the following system of equations:

$$T_+=t_++t_+G_{++}^0T_++t_+G_+^0T_{}$$

$$T_{}=t_{}+t_{}G_{}^0T_{}+t_{}G_+^0T_+,$$

where $`G_{\alpha \beta }^0`$ is the bare propagator between the sites with energies $`E_\pm `$. The solution of this system is of the following form

$`T_\pm ={\displaystyle \frac{t_\pm +t_\pm G_\pm ^0t_\pm }{(1t_+G_{++}^0)(1t_{}G_{}^0)G_+^0G_+^0t_+t_{}}}=`$
(51) $`{\displaystyle \frac{t_{}^1+G_\pm ^0}{(t_+^1G_{++}^0)(t_{}^1G_{}^0)G_+^0G_+^0}}.`$
Thus, by comparing this functional two-pole structure and the “Hubbard III” solution ,

$$\mathrm{\Sigma }_\sigma (\omega )=\omega F_\sigma (\omega )$$

, it is possible to identify the “scattering corrections” and “resonance broadening corrections” in the following way:

$$F_\sigma (\omega )=\frac{\omega (\omega U)(\omega Un_\sigma )A_\sigma (\omega )}{\omega U(1n_\sigma )A_\sigma (\omega )}$$

$$A_\sigma (\omega )=Y_\sigma (\omega )+Y_\sigma (\omega )Y_\sigma ^{}(U\omega )$$

$$Y_\sigma =F_\sigma (\omega )G_{0\sigma }^1(\omega );G_{0\sigma }(\omega )=N^1\underset{k}{}G_{k\sigma }(\omega )$$

If we put $`A_\sigma (\omega )=0`$, we immediately obtain the “Hubbard I” solution

(52)
$$G_\sigma ^{(H1)}(k,\omega )\frac{n_\sigma }{\omega Uϵ(k)n_\sigma }+\frac{1n_\sigma }{\omega ϵ(k)(1n_\sigma )}$$

Despite that this solution is exact in the atomic limit ( $`t_{ij}=0`$), the ”Hubbard I” solution has many serious drawbacks. The corresponding spectral function consists of two $`\delta `$-function peaks. The ”Hubbard III” solution includes several corrections, including scattering corrections which broadens the peaks and shift them when $`U`$ is changed.
The “alloy analogy” approximation corresponds to $`A_\sigma (\omega )Y_\sigma (\omega )`$. An interesting analysis of the ”Hubbard III” solution was performed in paper. The Hubbard sub-band structure was obtained in an analytic form in the ”Hubbard III” approximation, using the Lorentzian form for the density of states for non-interacting electrons. This resulted in an analytical form for the self-energy and the density of states for interacting electrons. Note that the “Hubbard III” self-energy operator $`\mathrm{\Sigma }_\sigma (\omega )`$ is local, i.e. does not depend on the quasi-momentum. Another drawback of this solution is a very inconvenient functional representation of elastic and inelastic scattering processes.
The conceptually new approach to the theory of very strong but finite electron correlation for the Hubbard model was proposed by Roth . She clarified microscopically the origination of the two-pole solution of the single-particle GF in the strongly correlated limit

(53) $`G_\sigma ^{(R)}(k,\omega )`$
$`{\displaystyle \frac{n_\sigma }{\omega Uϵ(k)n_\sigma W_{k\sigma }(1n_\sigma )}}+{\displaystyle \frac{1n_\sigma }{\omega ϵ(k)(1n_\sigma )n_\sigma W_{k\sigma }}}`$
We see that, in addition to a band narrowing effect, there is an energy shift $`W_{k\sigma }`$ given by

$`n_\sigma (1n_\sigma )W_{k\sigma }={\displaystyle \underset{ij}{}}t_{ij}<a_{i\sigma }^{}a_{j\sigma }(1n_{i\sigma }n_{j\sigma })>{\displaystyle \underset{ij}{}}t_{ij}\mathrm{exp}[ik(ji)]`$
(54) $`(n_\sigma ^2<n_{i\sigma }n_{j\sigma }>+<a_{j\sigma }^{}a_{i\sigma }^{}a_{j\sigma }a_{i\sigma }>+<a_{j\sigma }^{}a_{j\sigma }^{}a_{i\sigma }a_{i\sigma }>)`$
This energy shift corrects the situation with the ”Hubbard I” spectral function and recovers, in principle, the possibility of describing the ferromagnetic solution. Thus, the Roth solution gives an improved version of ”Hubbard I” two-pole solution and includes the band shift, that is most important in the case of a nearly-half-filled band. It is worth noting that this result was a very unusual fact from the point of view of the standard Fermi-liquid approach, showing that the naive one-electron approximation of band structure calculations is not valid for the description of electron correlations of lattice fermions.
It is this feature - the strong modification of single-particle states by many-body correlation effects - whose importance we wish to emphasize here.
Various attempts were made to describe the properties of the Hubbard model in both the strong and weak coupling regimes and to find a better solution ( e.g.  \-  ). Different schemes of decoupling of the equations of motion for the GFs analysed and compared in paper, when calculating the electron contribution to the cohesive energy in a narrow band system. These calculations showed importance of the correlation effects and the right scheme of approximation.
Thus, a sophisticated many-body technique is to be used for calculating the excitation spectra and other characteristics at finite temperatures. We shall show here following papers , that the IGF method permits us to improve substantially both the solutions, Hubbard and Roth, by defining the correct Generalized Mean Fields for the Hubbard model.

#### 4.2.2 Single Impurity Anderson Model (SIAM)

The Hamiltonian of SIAM can be written in the form

(55) $`H={\displaystyle \underset{k\sigma }{}}ϵ_kc_{k\sigma }^{}c_{k\sigma }+{\displaystyle \underset{\sigma }{}}E_{0\sigma }f_{0\sigma }^{}f_{0\sigma }+U/2{\displaystyle \underset{\sigma }{}}n_{0\sigma }n_{0\sigma }+`$
$`{\displaystyle \underset{k\sigma }{}}V_k(c_{k\sigma }^{}f_{0\sigma }+f_{0\sigma }^{}c_{k\sigma })`$
where $`c_{k\sigma }^{}`$ and $`f_{0\sigma }^{}`$ are, respectively, the creation operators for conduction and localized electrons; $`ϵ_k`$ is the conduction electron energy, $`E_{0\sigma }`$ is the localized electron energy level, and $`U`$ is the intra-atomic Coulomb interaction at the impurity site; $`V_k`$ represents the $`s(d)f`$ hybridization interaction term and was written in paper in the following form

(56)
$$V_k=\frac{1}{\sqrt{N}}\underset{j}{}V_f(R_j)exp(ikR_j)$$

The hybridization matrix element is

$$V_f(R_j)=\psi _k^{}(\stackrel{}{r})H^{HF}\varphi (\stackrel{}{r}\stackrel{}{R}_j)𝑑r$$

The use of Hartree-Fock term here is notable, since it justifies the initial treatment of SIAM in entirely in the H-F approximation. A number of approaches for SIAM and other correlated electronic systems was proposed, aimed at answering the Anderson question: ”…whether a real many-body theory would give answer radically different from the Hartree-Fock results?”.
Our goal is to propose a new combined many-body approach for the description of many-body quasi-particle dynamics of SIAM at finite temperatures. The interplay and competition of the kinetic energy ($`ϵ_k`$), potential energy ($`U`$), and hybridization ($`V`$) substantially influence the electronic spectrum. The renormalized electron energies are temperature-dependent, and electronic states have finite lifetimes. These effects are described most suitable by the Green functions method. The purpose of the present approach is to find the electronic quasi-particle spectrum renormalized by interactions (U- and V-terms) in a wide range of temperatures and model parameters and to calculate explicitly the damping of the electronic states.

#### 4.2.3 Periodic Anderson Model (PAM)

Let us now consider a lattice generalization of SIAM, the so-called periodic Anderson model (PAM). The basic assumption of the periodic impurity Anderson model is the presence of two well-defined subsystems, i.e. the Fermi sea of nearly free conduction electrons and the localized impurity orbitals embedded into the continuum of conduction electron states ( in rare-earth compounds, for instance, the continuum is actually a mixture of $`s`$, $`p`$, and $`d`$ states, and the localized orbitals are $`f`$ states). The simplest form of PAM

(57) $`H={\displaystyle \underset{k\sigma }{}}ϵ_kc_{k\sigma }^{}c_{k\sigma }+{\displaystyle \underset{i\sigma }{}}E_0f_{i\sigma }^{}f_{i\sigma }+U/2{\displaystyle \underset{i\sigma }{}}n_{i\sigma }n_{i\sigma }+`$
$`{\displaystyle \frac{V}{\sqrt{N}}}{\displaystyle \underset{ik\sigma }{}}(exp(ikR_i)c_{k\sigma }^{}f_{i\sigma }+exp(ikR_i)f_{i\sigma }^{}c_{k\sigma })`$
assumes a one-electron energy level $`E_0`$, hybridization interaction $`V`$, and the Coulomb interaction $`U`$ at each lattice site. Using the transformation

$$c_{k\sigma }^{}=\frac{1}{\sqrt{N}}\underset{j}{}exp(ikR_j)c_{j\sigma }^{};c_{k\sigma }=\frac{1}{\sqrt{N}}\underset{j}{}exp(ikR_j)c_{j\sigma }$$

the Hamiltonian ( 57) can be rewritten in the Wannier representation:

(58) $`H={\displaystyle \underset{ij\sigma }{}}t_{ij}c_{i\sigma }^{}c_{j\sigma }+{\displaystyle \underset{i\sigma }{}}E_0f_{i\sigma }^{}f_{i\sigma }+U/2{\displaystyle \underset{i\sigma }{}}n_{i\sigma }n_{i\sigma }+`$
$`V{\displaystyle \underset{i\sigma }{}}(c_{i\sigma }^{}f_{i\sigma }+f_{i\sigma }^{}c_{i\sigma })`$
If one retains the $`k`$-dependence of the hybridization matrix element $`V_k`$ in (58), the last term in the PAM Hamiltonian describing the hybridization interaction between the localized impurity states and extended conduction states and containing the essence of a specificicity of the Anderson model, is as follows

$$\underset{ij\sigma }{}V_{ij}(c_{i\sigma }^{}f_{i\sigma }+f_{i\sigma }^{}c_{i\sigma });V_{ij}=\frac{1}{N}\underset{k}{}V_kexp[ik(R_jR_i)]$$

The on-site hybridization $`V_{ii}`$ is equal to zero for symmetry reasons. A detailed analysis of the hybridization problem from a general point of view and in the context of PAM was made in paper. The Hamiltonian of PAM in the Bloch representation takes the form

(59) $`H={\displaystyle \underset{k\sigma }{}}ϵ_kc_{k\sigma }^{}c_{k\sigma }+{\displaystyle \underset{i\sigma }{}}E_kf_{k\sigma }^{}f_{k\sigma }+U/2{\displaystyle \underset{i\sigma }{}}n_{i\sigma }n_{i\sigma }+`$
$`{\displaystyle \underset{k\sigma }{}}V_k(c_{k\sigma }^{}f_{k\sigma }+f_{k\sigma }^{}c_{k\sigma })`$
Note that as compared to the SIAM, the PAM has its own specific features. This can lead to peculiar magnetic properties for concentrated rare-earth systems where the criterion for magnetic ordering depends on the competition between indirect RKKY-type interaction ( not included into SIAM ) and the Kondo-type singlet-site screening ( contained in SIAM ). The inclusion of inter-impurity correlations makes the problem more difficult. Since these inter-impurity effects play an essential role in physical behaviour of real systems,, it is instructive to consider the two-impurity Anderson model (TIAM) too.

#### 4.2.4 Two-Impurity Anderson Model ( TIAM )

The two-impurity Anderson model was considered by Alexander and Anderson. They put forward a theory which introduces the impurity-impurity interaction within a game of parameters. The Hamiltonian of TIAM reads

(60) $`H={\displaystyle \underset{ij\sigma }{}}t_{ij}c_{i\sigma }^{}c_{j\sigma }+{\displaystyle \underset{i=1,2\sigma }{}}E_{0i}f_{i\sigma }^{}f_{i\sigma }+U/2{\displaystyle \underset{i=1,2\sigma }{}}n_{i\sigma }n_{i\sigma }+`$
$`{\displaystyle \underset{i\sigma }{}}(V_{ki}c_{i\sigma }^{}f_{i\sigma }+V_{ik}f_{i\sigma }^{}c_{i\sigma })+{\displaystyle \underset{\sigma }{}}(V_{12}f_{1\sigma }^{}f_{2\sigma }+V_{21}f_{2\sigma }^{}f_{1\sigma })`$
where $`E_{0i}`$ are the position energies of localized states ( for simplicity, we consider identical impurities and $`s`$-type (i.e. non-degenerate) orbitals: $`E_{01}=E_{02}=E_0`$. Let us recall that the hybridization matrix element $`V_{ik}`$ was defined in (56). As for the TIAM, the situation with the right definition of the parameters $`V_{12}`$ and $`V_{ik}`$ is not very clear. The definition of $`V_{12}`$ in is the following

$$V_{12}=V_{21}^{}=\varphi _1^{}(\stackrel{}{r})H_f\varphi _2(\stackrel{}{r})𝑑r$$

(now $`H_f`$ without ”H-F” mark). The essentially local character of the Hamiltonian $`H_f`$ clearly shows that $`V_{12}`$ describes the direct coupling between nearest neighboring sites ( for a detailed discussion see where the hierarchy of the Anderson models was discussed too ).

## 5 Effective and Generalized Mean Fields

### 5.1 Molecular Field Approximation

The most common technique for studying the subject of interacting many-particle systems is to use the mean field theory. This approximation is especially popular in the theory of magnetism . Nevertheless, it was pointed  that

> ”the Weiss molecular field theory plays an enigmatic role in statistical mechanics of magnetism”.

To calculate the susceptibility and other characteristic functions of a system of localized magnetic moments, with a given interaction Hamiltonian, the approximation, termed the ”molecular field approximation” was used widely. However, it is not an easy task to give the formal unified definition what the mean field is. In a sense, the mean field is the umbrella term for a variety of theoretical methods of reducing the many-particle problem to the single-particle one. Mean field theory, that approximates the behaviour of a system by ignoring the effect of fluctuations and those spin correlations which dominate the collective properties of the ferromagnet usually provides a starting and estimating point only, for studying phase transitions. The mean field theories miss important features of the dynamics of a system. The main intention of the mean field theories, starting from the works of van der Waals and P.Weiss, is to take into account the cooperative behaviour of a large number of particles. It is well known that earlier theories of phase transitions based on the ideas of van der Waals and Weiss lead to predictions which are qualitatively at variance with results of measurements near the critical point. Other variants of simplified mean field theories such as the Hartree-Fock theory for electrons in atoms, etc lead to discrepancies of various kinds too. It is therefore natural to analyze the reasons for such drawbacks of earlier variants of the mean field theories.

### 5.2 Effective Field Theories

A number of effective field theories which are improved versions of the ”molecular field approximation” were proposed. It is the purpose of this study to stress a specificity of strongly correlated many-particle systems on a lattice contrary to continuum (uniform) systems. Although many important questions remain still unresolved, a vision of useful synthesis begins to emerge. As a workable eye-guide , the set of mean field theories ( most probably incomplete ) is shown in Table 1.

The meaning of many these entries and terms will become clearer in the forthcoming discussion and will put them in a clearer perspective. My main purpose is to elucidate ( at least in the mathematical structure ) and to give plausible arguments for the tendency, which expounded in Table 1. This tendency shows the following. The earlier concepts of molecular field were described in terms of a functional of mean magnetic moments (in magnetic terminology ) or mean particle densities ( Hartree-Fock field ). The corresponding mean-field functional $`F[<n>,<S^z>]`$ describes the uniform mean field.
Actually, the Weiss model was not based on discrete ”spins” as is well known, but the uniformity of the mean internal field was the most essential feature of the model. In the modern language, one should assume that the interaction between atomic spins $`S_i`$ and its neighbors is equivalent to a mean ( or molecular ) field, $`M_i=\chi _0[h_i^{(ext)}+h_i^{(mf)}]`$ and that the molecular field $`h_i^{(mf)}`$ is of the form $`h^{(mf)}=_iJ(R_{ji})<S_i>`$ (above $`T_c`$ ). Here $`h^{ext}`$ is an applied conjugate field, $`\chi _0`$ is the response function, and $`J(R_{ji})`$ is an interaction. In other words, the mean field approximation reduces the many-particle problem to a single-site problem in which a magnetic moment at any site can be either parallel or antiparallel to the total magnetic field composed of the applied field and the molecular field. The average interaction of $`i`$ neighbors was taken into account only, and the fluctuations were neglected. One particular example, where the mean field theory works relatively well is the homogeneous structural phase transitions; in this case the fluctuations are confined in phase space.
The next important step was made by L. Neel . He conjectured that the Weiss internal field might be either positive or negative in sign. In the latter case, he showed that below a critical temperature ( Neel temperature ) an ordered arrangement of equal numbers of oppositely directed atomic moments could be energetically favorable. This new magnetic structure was termed antiferromagnetism. It was conjectured that the two-sublattice Neel ( classical ) ground state is formed by local staggered internal mean fields.
There is a number of the ”correlated effective field” theories, that tend to repair the limitations of simplified mean field theories. The remarkable and ingenious one is the Onsager ”reaction field approximation”. He suggested that the part of the molecular field on a given dipole moment which comes from the reaction of neighboring molecules to the instantaneous orientation of the moment should not be included into the effective orienting field. This ”reaction field” simply follows the motion of the moment and thus does not favor one orientation over another. The meaning of the mean field approximation for the spin glass problem is very interesting but specific, and we will not discuss it here. A single-site molecular-field model for randomly dilute ferro- and antiferromagnets in the framework of the double-time thermal GFs was presented in paper.

### 5.3 Generalized Mean Fields

It was shown  that mean-field approximations, for example the molecular field approximation for a spin system, the Hartree-Fock approximation and the BCS-Bogoliubov approximation for an electron system are universally formulated by the Peierls-Bogoliubov-Feynman (PBF) inequality:

$`\beta ^1ln(Tre^{(\beta H)})`$
(61) $`\beta ^1ln(Tre^{(\beta H^{mf})})+{\displaystyle \frac{Tre^{(\beta H^{mf})}(HH^{mf})}{Tre^{(\beta H^{mf})}}}`$
Here $`F`$ is the free energy, and $`H^{mf}`$ is a ”trial” or a ”mean field” approximating Hamiltonian. This inequality gives the upper bound of the free energy of a many-body system. It is important to emphasize that the BCS-Bogoliubov theory of superconductivity , was formulated on the basis of a trial Hamiltonian which consists of a quadratic form of creation and annihilation operators, including ”anomalous” ( off-diagonal ) averages . The functional of the mean field ( for the superconducting single-band Hubbard model ) is of the following form :

(62)
$$\mathrm{\Sigma }_\sigma ^c=U\left(\begin{array}{cc}<a_{i\sigma }^{}a_{i\sigma }>& <a_{i\sigma }a_{i\sigma }>\\ <a_{i\sigma }^{}a_{i\sigma }^{}>& <a_{i\sigma }^{}a_{i\sigma }>\end{array}\right)$$

The ”anomalous” off-diagonal terms fix the relevant BCS-Bogoliubov vacuum and select the appropriate set of solutions.
Another remark about the BCS-Bogolubov mean-field approach is instructive. Speaking in physical terms, this theory involves a condensation correctly, in spite that such a condensation cannot be obtained by an expansion in the effective interaction between electrons. Other mean field theories, e.g. the Weiss molecular field theory and the van der Waals theory of the liquid-gas transition are much less reliable. The reason why a mean-field theory of the superconductivity in the BCS-Bogoliubov form is successful would appear to be that the main correlations in metal are governed by the extreme degeneracy of the electron gas. The correlations due to the pair condensation, although they have dramatic effects, are weak ( at least in the ordinary superconductors ) in comparison with the typical electron energies, and may be treated in an average way with a reasonable accuracy. All above remarks have relevance to ordinary low-temperature superconductors. In high-$`T_c`$ superconductors, the corresponding degeneracy temperature is much lower, and the transition temperatures are much higher. In addition, the relevant interaction responsible for the pairing and its strength are unknown. From this point of view, the high-$`T_c`$ systems are more complicated. It should be clarified what governs the scale of temperatures, i.e. critical temperature, degeneracy temperature, interaction strength or their complex combination, etc. In this way a useful insight into this extremely complicated problem would be gained.
Generalization of the molecular field approximation on the basis of the PBF inequality is possible when we know a particular solution of the model (e.g., for one-dimensional Ising model we know the exact solution in the field). One can use this solution to get a better approximation than the mean field theory. There are some other methods of improvement of the molecular field theory , . Unfortunately, these approaches are nonsystematic.
From the point of view of quantum many-body theory, the problem of adequate introduction of mean fields for system of many interacting particles can be most consistently investigated in the framework of the IGF method. A correct calculation of the quasi-particle spectra and their damping, particularly, for systems with a complicated spectrum and strong interaction  reveals, as it will be shown below, that the generalized mean fields can have very complicated structure which cannot be described by a functional of the mean-particle density.
To illustrate the actual distinction of description of the generalized mean field in the equation-of-motion method for the double-time Green functions, let us compare the two approaches, namely, that of Tyablikov  and of Callen . We shall consider the Green function $`<<S^+|S^{}>>`$ for the isotropic Heisenberg model

(63)
$$H=\frac{1}{2}\underset{ij}{}J(ij)\stackrel{}{S}_i\stackrel{}{S}_j$$

The equation of motion (12) for the spin Green function is of the form

(64) $`\omega <<S_i^+|S_j^{}>>_\omega =`$
$`2<S^z>\delta _{ij}+{\displaystyle \underset{g}{}}J(ig)<<S_i^+S_g^zS_g^+S_i^z|S_j^{}>>_\omega `$
The Tyablikov decoupling expresses the second-order GF in terms of the first (initial) GF:

(65)
$$<<S_i^+S_g^z|S_j^{}>>=<S^z><<S_i^+|S_j^{}>>$$

This approximation is an RPA-type; it does not lead to the damping of spin wave excitations (cf. (39) )

(66)
$$E(q)=\underset{g}{}J(ig)<S^z>\mathrm{exp}[i(\stackrel{}{R}_i\stackrel{}{R}_g)\stackrel{}{q}]=2<S^z>(J_0J_q)$$

The reason for this is rather transparent. This decoupling does not take into account the inelastic magnon-magnon scattering processes. In a sense, the Tyablikov approximation consists of approximating the commutation relations of spin operators to the extent of replacing the commutation relation $`[S_i^+,S_j^{}]_{}=2S_i^z\delta _{ij}`$ by $`[S_i^+,S_j^{}]_{}=2<S^z>\delta _{ij}`$ .
Callen  has proposed an improved decoupling approximation in the method of Tyablikov in the following form:

$$<<S_g^zS_f^+|B>><S^z><<S_f^+|B>>\alpha <S_g^{}S_f^+><<S_g^+|B>>$$

Here $`0\alpha 1`$. To clarify this point, it should be reminded that for spin $`1/2`$ ( the procedure was generalized by Callen to an arbitrary spin), the spin operator $`S^z`$ can be written as $`S_g^z=SS_g^{}S_g^+`$ or $`S_g^z=\frac{1}{2}(S_g^+S_g^{}S_g^{}S_g^+)`$. It is easy to show that

$$S_g^z=\alpha S+\frac{1\alpha }{2}S_g^+S_g^{}\frac{1+\alpha }{2}S_g^{}S_g^+$$

The operator $`S_g^{}S_g^+`$ represents the deviation of $`<S^z>`$ from $`S`$. In the low-temperature region, this deviation is small, and $`\alpha 1`$. Similarly, the operator $`\frac{1}{2}(S_g^+S_g^{}S_g^{}S_g^+)`$ represents the deviation of $`<S^z>`$ from 0. Thus, when $`<S^z>`$ approaches to zero, one can expect that $`\alpha 0`$. Thus, in this way, it is possible to obtain a correction to the Tyablikov decoupling with either a positive or negative sign, or no correction at all, or any intermediate value, depending on the choice of $`\alpha `$. The above Callen arguments are not rigorous , for, although the difference in the operators $`S^+S^{}`$ and $`S^{}S^+`$ is small if $`<S^z>0`$, each operator makes a contribution of the order of $`S`$, and it is each operator which is treated approximately, not the difference. There are some other drawbacks of the Callen decoupling scheme. Nevertheless, the Callen decoupling was the first conceptual attempt to introduce the interpolation decoupling procedure. Let us note that the choice of $`\alpha =0`$ over the entire temperature range is just the Tyablikov decoupling (65).
The energy spectrum for the Callen decoupling is given by

(67)
$$E(q)=2<S^z>\left((J_0J_q)+\frac{<S^z>}{NS^2}\underset{k}{}[J(k)J(kq)]N(E(k))\right)$$

Here $`N(E(k))`$ is the Bose distribution function $`N(E(k))=[\mathrm{exp}(E(k)\beta )1]^1`$. This is an implicit equation for $`N(E(k))`$, involving the unknown quantity $`<S^z>`$ . For the latter an additional equation is given . Thus, both these equations constitute a set of coupled equations which must be solved self-consistently for $`<S^z>`$.
This formulation of the Callen decoupling scheme displays explicitly the tendency of the improved description of the mean field. In a sense, it is possible to say that the Callen work dates really the idea of the generalized mean field within the equation-of-motion method for double-time GFs, however, in a semi-intuitive form. The next essential steps were made by Plakida  for the Heisenberg ferromagnet and by Kuzemsky , for the Hubbard model. As was mentioned above, the correct definition of Generalized Mean Fields depends on the condition of the problem, the strength of interaction, the choice of relevant operators, and on the symmetry requirements.

### 5.4 Symmetry Broken Solutions

In many-body interacting systems, the symmetry is important in classifying different phases and in understanding the phase transitions between them . According to Bogoliubov ( cf. refs. , ,) in each condensed phase, in addition to the normal process, there is an anomalous process (or processes) which can take place because of the long-range internal field, with a corresponding propagator. Additionally, the Goldstone theorem states that, in a system in which a continuous symmetry is broken ( i.e. a system such that the ground state is not invariant under the operations of a continuous unitary group whose generators commute with the Hamiltonian ), there exists a collective mode with frequency vanishing, as the momentum goes to zero. For many-particle systems on a lattice, this statement needs a proper adaptation. In the above form, the Goldstone theorem is true only if the condensed and normal phases have the same translational properties. When translational symmetry is also broken, the Goldstone mode appears at a zero frequency but at nonzero momentum, e.g., a crystal and a helical spin-density-wave (SDW) ordering (see for discussion-).
The anomalous propagators for an interacting many-fermion system corresponding to the ferromagnetic (FM) , antiferromagnetic (AFM), and superconducting (SC) long-range ordering are given by

(68) $`FM:G_{fm}<<a_{k\sigma };a_{k\sigma }^{}>>`$
$`AFM:G_{afm}<<a_{k+Q\sigma };a_{k+Q^{}\sigma ^{}}^{}>>`$
$`SC:G_{sc}<<a_{k\sigma };a_{k\sigma }>>`$
In the SDW case, a particle picks up a momentum $`QQ^{}`$ from scattering against the periodic structure of the spiral ( nonuniform ) internal field, and has its spin changed from $`\sigma `$ to $`\sigma ^{}`$ by the spin-aligning character of the internal field. The Long-Range-Order (LRO) parameters are:

(69) $`FM:m=1/N{\displaystyle \underset{k\sigma }{}}<a_{k\sigma }^{}a_{k\sigma }>`$
$`AFM:M_Q={\displaystyle \underset{k\sigma }{}}<a_{k\sigma }^{}a_{k+Q\sigma }>`$
$`SC:\mathrm{\Delta }={\displaystyle \underset{k}{}}<a_k^{}a_k^{}>`$
It is important to note that the long-range order parameters are functions of the internal field, which is itself a function of the order parameter. There is a more mathematical way of formulating this assertion. According to the paper , the notion ”symmetry breaking” means that the state fails to have the symmetry that the Hamiltonian has.
A true breaking of symmetry can arise only if there are infinitesimal ”source fields”. Indeed, for the rotationally and translationally invariant Hamiltonian, suitable source terms should be added:

(70) $`FM:\epsilon \mu _BH_x{\displaystyle \underset{k\sigma }{}}a_{k\sigma }^{}a_{k\sigma }`$
$`AFM:\epsilon \mu _BH{\displaystyle \underset{kQ}{}}a_{k\sigma }^{}a_{k+Q\sigma }`$
$`SC:\epsilon v{\displaystyle \underset{k}{}}(a_k^{}a_k^{}+a_ka_k)`$
where $`\epsilon 0`$ is to be taken at the end of calculations.
For example, broken symmetry solutions of the SDW type imply that the vector $`Q`$ is a measure of the inhomogeneity or breaking of translational symmetry. The Hubbard model is a very interesting tool for analyzing the symmetry broken concept. It is possible to show that antiferromagnetic state and more complicated states ( e.g. ferrimagnetic) can be made eigenfunctions of the self-consistent field equations within an ”extended” mean-field approach, assuming that the ”anomalous” averages $`<a_{i\sigma }^{}a_{i\sigma }>`$ determine the behaviour of the system on the same footing as the ”normal” density of quasi-particles $`<a_{i\sigma }^{}a_{i\sigma }>`$. It is clear, however, that these ”spin-flip” terms break the rotational symmetry of the Hubbard Hamiltonian. For the single-band Hubbard Hamiltonian, the averages $`<a_{i\sigma }^{}a_{i,\sigma }>=0`$ because of the rotational symmetry of the Hubbard model. The inclusion of ”anomalous” averages leads to the so-called ”unresricted” H-F approximation (UHFA). This type of approximation was used sometimes also for the single-band Hubbard model for calculating the density of states. For this aim, the following definition of UHFA

(71)
$$n_{i\sigma }a_{i\sigma }<n_{i\sigma }>a_{i\sigma }<a_{i\sigma }^{}a_{i\sigma }>a_{i\sigma }$$

was used. Thus, in addition to the standard H-F term, the new so-called “spin-flip” terms are retained. This example clearly shows that the structure of mean field follows from the specificity of the problem and should be defined in a proper way. So, one needs a properly defined effective Hamiltonian $`H_{\mathrm{eff}}`$. In paper  we thoroughly analyzed the proper definition of the irreducible GFs which includes the “spin-flip” terms for the case of itinerant antiferromagnetism of correlated lattice fermions. For the single-orbital Hubbard model, the definition of the ”irreducible” part should be modified in the following way:

$`{}_{}{}^{(ir)}<<a_{k+p\sigma }a_{p+q\sigma }^{}a_{q\sigma }|a_{k\sigma }^{}>>_\omega =<<a_{k+p\sigma }a_{p+q\sigma }^{}a_{q\sigma }|a_{k\sigma }^{}>>_\omega `$
(72) $`\delta _{p,0}<n_{q\sigma }>G_{k\sigma }<a_{k+p\sigma }a_{p+q\sigma }^{}><<a_{q\sigma }|a_{k\sigma }^{}>>_\omega `$
From this definition it follows that this way of introduction of the IGF broadens the initial algebra of operators and the initial set of the GFs. This means that the “actual” algebra of operators must include the spin-flip terms from the beginning, namely: $`(a_{i\sigma }`$, $`a_{i\sigma }^{}`$, $`n_{i\sigma }`$, $`a_{i\sigma }^{}a_{i\sigma })`$. The corresponding initial GF will be of the form

$$\left(\begin{array}{cc}<<a_{i\sigma }|a_{j\sigma }^{}>>& <<a_{i\sigma }|a_{j\sigma }^{}>>\\ <<a_{i\sigma }|a_{j\sigma }^{}>>& <<a_{i\sigma }|a_{j\sigma }^{}>>\end{array}\right)$$

With this definition, one introduces the so-called anomalous (off-diagonal) GFs which fix the relevant vacuum and select the proper symmetry broken solutions. In fact, this approximation was investigated earlier by Kishore and Joshi . They clearly pointed out that they assumed a system to be magnetised in the $`x`$ direction instead of the conventional $`z`$ axis.
The problem of finding the ferromagnetic and antiferromagnetic ”symmetry broken” solutions of the correlated lattice fermion models within IGF method was investigated in ref. . A unified scheme for the construction of Generalized Mean Fields ( elastic scattering corrections ) and self-energy ( inelastic scattering ) in terms of the Dyson equation was generalized in order to include the ”source fields”. The ”symmetry broken” dynamic solutions of the Hubbard model which correspond to various types of itinerant antiferromagnetism were discussed. This approach complements previous studies of microscopic theory of the Heisenberg antiferromagnet  and clarifies the concepts of Neel sublattices for localized and itinerant antiferromagnetism and ”spin-aligning fields” of correlated lattice fermions.

## 6 Quasi-Particle Many Body Dynamics

In this Section, we discuss the microscopic view of a dynamic behaviour of interacting many-body systems on a lattice. It was recognized for many years that the strong correlation in solids exist between the motions of various particles ( electrons and ions, i.e. the fermion and boson degrees of freedom ) which arise from the Coulomb forces. The most interesting objects are metals and their compounds. They are invariant under the translation group of a crystal lattice and have lattice vibrations as well as electron degrees of freedom. There are many evidences for the importance of many-body effects in these systems. Within the Landau semi-phenomenological theory it was suggested that the low-lying excited states of an interacting Fermi gas can be described in terms of a set of ”independent quasi-particles”. However, this was a phenomenological approach and did not reveal the nature of relevant interactions.

### 6.1 Green Function Picture of Quasi-Particles

An alternative way of viewing quasi-particles, more general and consistent, is through the Green function scheme of many-body theory, which we sketch below for completeness and for pedagogical reasons.
We should mention that there exist a big variety of quasi-particles in many-body systems. At sufficiently low temperatures, few quasi-particles are excited, and therefore this dilute quasi-particle gas is nearly a non-interacting gas in the sense that the quasi-particles rarely collide. The success of the quasi-particle concept in an interacting many-body system is particularly striking because of a great number of various applications. However, the range of validity of the quasi-particle approximation, especially for strongly interacting lattice systems, was not discussed properly in many cases. In systems like simple metals, quasi-particles constitute long-lived, weakly interacting excitations, since their intrinsic decay rate varies as the square of the dispersion law, thereby justifying their use as the building blocks for the low-lying excitation spectrum.
Unfortunately, there are many strongly correlated systems on a lattice for which we do not have at present the truly the first-principles proof of a similar correspondence of the low-lying excited states of noninteracting and interacting systems, adiabatic switching on of the interaction, a simple effective mass spectrum, long lifetimes of quasi-particles, etc. These specific features of strongly correlated systems are the main reason of why the usual perturbation theory starting from noninteracting states does not work properly. Many other subtle nonanalytic effects which are present even in normal systems have the similar nature . This lack of a rigorous foundation for the theory of strongly interacting systems on a lattice is not only a problem of the mathematical perfectionism, but also that of the correct physics of interacting systems.
As we mentioned earlier, to describe a quasi-particle correctly, the Green functions method is a very suitable and useful tool. What concerns us here are formal expression for the single-particle GF (38) and the corresponding quasi-particle excitation spectrum. From the equation ( 22) it is thus seen that the GF is completely determined by the spectral weight function $`A(\omega )`$. The spectral weight function reflects the microscopic structure of the system under consideration. The other term in ( 22) is a separation of the purely statistical aspects of GF. From the equation ( 18) it follows that the spectral weight function can be written formally in terms of many-particle eigenstates. Its Fourier transform origination ( 16) is then the density of states that can be reached by adding or removing a particle of a given momentum and energy.
Consider a system of interacting fermions as an example. For a noninteracting system, the spectral weight function of the single-particle GF $`G_k(\omega )=<<a_{k\sigma };a_{k\sigma }^{}>>`$ has the simple peaked structure

$$A_k(\omega )\delta (\omega ϵ_k)$$

. For an interacting system, the spectral function $`A_k(\omega )`$ has no such a simple peaked structure, but it obeys the following conditions

$$A_k(\omega )0;A_k(\omega )𝑑\omega =<[a_{k\sigma },a_{k\sigma }^{}]_+>=1$$

Thus, we can see from these expressions that for a noninteracting system, the sum rule is exhausted by a single peak. A sharply peaked spectral function for an interacting system means a long-lived single-particle-like excitation. Thus, the spectral weight function was established here as the physically significant attribute of GF. The question of how best to extract it from a microscopic theory is the main aim of the present review.
The GF for the non-interacting system is $`G_k(\omega )=(\omega ϵ_k)^1`$. For a weakly interacting Fermi system, we have $`G_k(\omega )=(\omega ϵ_kM_k(\omega ))^1`$ where $`M_k(\omega )`$ is the mass operator. Thus, for a weakly interacting system, the $`\delta `$-function for $`A_k(\omega )`$ is spread into a peak of finite width due to the mass operator. We have

$$M_k(\omega \pm iϵ)=ReM_k(\omega )ImM_k(\omega )=\mathrm{\Delta }_k(\omega )\mathrm{\Gamma }_k(\omega )$$

The single-particle GF can be written in the form

(73)
$$G_k(\omega )=\{\omega [ϵ_k+\mathrm{\Delta }_k(\omega )]\pm \mathrm{\Gamma }_k(\omega )\}^1$$

In the weakly interacting case, we can thus find the energies of quasi-particles by looking for the poles of single-particle GF (73)

$$\omega =ϵ_k+\mathrm{\Delta }_k(\omega )\pm \mathrm{\Gamma }_k(\omega )$$

. The dispersion relation of a quasi-particle

$$ϵ(k)=ϵ_k+\mathrm{\Delta }_k[ϵ(k)]\pm \mathrm{\Gamma }_k[ϵ(k)]$$

and the lifetime $`1/\mathrm{\Gamma }_k`$ then reflects the inter-particle interaction. It is easy to see the connection between the width of the spectral weight function and decay rate. We can write

(74) $`A_k(\omega )=(\mathrm{exp}(\beta \omega )+1)^1(i)[G_k(\omega +iϵ)G_k(\omega iϵ)]=`$
$`(\mathrm{exp}(\beta \omega )+1)^1{\displaystyle \frac{2\mathrm{\Gamma }_k(\omega )}{[\omega (ϵ_k+\mathrm{\Delta }_k(\omega ))]^2+\mathrm{\Gamma }_k^2(\omega )}}`$
In other words, for this case, the corresponding propagator can be written in the form

$$G_k(t)\mathrm{exp}(iϵ(k)t)\mathrm{exp}(\mathrm{\Gamma }_kt)$$

This form shows under which conditions, the time-development of an interacting system can be interpreted as the propagation of a quasi-particle with a reasonably well-defined energy and a sufficiently long lifetime. To demonstrate this, we consider the following conditions

$$\mathrm{\Delta }_k[ϵ(k)]ϵ(k);\mathrm{\Gamma }_k[ϵ(k)]ϵ(k)$$

Then we can write

(75)
$$G_k(\omega )=\frac{1}{[\omega ϵ(k)][1\frac{d\mathrm{\Delta }_k(\omega )}{d\omega }|_{\omega =ϵ(k)}]+i\mathrm{\Gamma }_k[ϵ(k)]}$$

where the renormalized energy of excitations is defined by

$$ϵ(k)=ϵ_k+\mathrm{\Delta }_k[ϵ(k)]$$

In this case, we have, instead of ( 74),

(76) $`A_k(\omega )=`$
$`[\mathrm{exp}(\beta ϵ(k))+1]^1[1{\displaystyle \frac{d\mathrm{\Delta }_k(\omega )}{d\omega }}|_{ϵ(k)}]^1{\displaystyle \frac{2\mathrm{\Gamma }(k)}{(\omega ϵ(k))^2+\mathrm{\Gamma }^2(k)}}`$
As a result, we find

(77) $`G_k(t)=<<a_{k\sigma }(t);a_{k\sigma }^{}>>=`$
$`=i\theta (t)\mathrm{exp}(iϵ(k)t)\mathrm{exp}(\mathrm{\Gamma }(k)t)[1{\displaystyle \frac{d\mathrm{\Delta }_k(\omega )}{d\omega }}|_{ϵ(k)}]^1`$
A widely known strategy to justify this line of reasoning is the perturbation theory. A detailed analysis of various successful approximations for the determination of excited states in the framework of the quasi-particle concept and the Green functions method in metals, semiconductors, and insulators was done in review paper.
There are examples of weakly interacting systems, i.g. the superconducting phase, which are not connected perturbatively with noninteracting systems. Moreover, the superconductor is a system in which the interaction between electrons qualitatively changes the spectrum of excitations. However, quasi-particles are still of use even in this case, due to the correct redefinition of the relevant generalized mean field which includes the anomalous averages (see (69)). In a strongly interacted system on a lattice with complex spectra, the concept of a quasi-particle needs a suitable adaptation and a careful examination. It is therefore useful to have the workable and efficient IGF method which, as we shall see, permits one to determine and correctly separate the elastic and inelastic scattering renormalizations through a correct definition of the generalized mean field and to calculate real quasi-particle spectra, including the damping and lifetime effects. A careful analysis and detailed presentations of the IGF method will provide an important step to the formulation of the consistent theory of strongly interacting systems and the justification of approximate methods presently used within equation-of-motion approaches. These latter remarks will not be substantiated until next Sections, but it is important to emphasize that the development which follows is not a merely formal exercise but essential for the proper and consistent theory of strongly interacting many-body systems on a lattice.

### 6.2 Spin-Wave Scattering Effects in Heisenberg Ferromagnet

In this Section, we briefly describe , mainly for pedagogical reasons, how the formulation of the quasi-particle picture depends in an essential way on an analysis of the sort introduced in Section 3.1. We consider here the most studied case of a Heisenberg ferromagnet with the Hamiltonian (63) and the equation of motion (64). In an earlier discussion in Sections 4.11 and 5.3, we described the Tyablikov decoupling procedure (65) based on replacing $`S_i^z`$ by $`<S_i^z>`$ in the last term of (64). We also discussed an alternative method of decoupling proposed by Callen (5.3). Both these decoupling procedures retain only the elastic spin-wave scattering effects. But for our purposes, it is essential to retain also the inelastic scattering effects, and therefore, we must carefully identify and separate the elastic and inelastic spin-wave scattering. This is directly related with the correct definition of generalized mean fields. Thus, the purpose of the present consideration is to justify the use of IGF method for the self-consistent theory of spin-wave interactions.
The irreducible part of GF is introduced according to the definition (28) as

(78)
$$^{(ir)}<<(S_i^+S_g^zS_g^+S_i^z)|S_j^{}>>=<<(S_i^+S_g^zS_g^+S_i^z)A_{ig}S_i^+A_{gi}S_g^+|S_j^{}>>$$

Here the unknown quantities $`A_{ig}`$ are defined on the basis of orthogonality constraint (29)

$$<[(S_i^+S_g^zS_g^+S_i^z)^{(ir)},S_j^{}]>=0$$

We have ($`ig`$)

(79)
$$A_{ig}=A_{gi}=\frac{2<S_i^zS_g^z>+<S_i^{}S_g^+>}{2<S^z>}$$

The definition (see eq.(31) ) of a generalized mean field GF $`G^{MF}`$ is given by the equation

(80)
$$\omega G_{ij}^{MF}=2<S^z>\delta _{ij}+\underset{g}{}J_{ig}A_{ig}(G_{ij}^{MF}G_{gj}^{MF})$$

From the Dyson equation in the form (35) we find

(81) $`M_{ij}=(P_{ij})^p=`$
$`<2S^z>^2{\displaystyle \underset{gl}{}}J_{ig}J_{lj}<<(S_i^+S_g^zS_g^+S_i^z)^{(ir)}|((S_i^+S_g^zS_g^+S_i^z)^{(ir)})^{}>>^{(p)}`$
where the proper $`(p)`$ part of the irreducible GF is defined by the equation (34)

$$P_{ij}=M_{ij}+\underset{gl}{}M_{ig}G_{gl}^{MF}P_{lj};M_{ij}=(P_{ij})^p$$

( in the diagrammatic language, this means that it has no parts connected by one $`G^{MF}`$-line). The formal solution of the Dyson equation is of the form (38):

(82) $`G_{ij}(\omega )=`$
$`2<S^z>N^1{\displaystyle \underset{k}{}}\mathrm{exp}[ik(R_iR_j)][\omega \omega (k)2<S^z>M_k(\omega )]^1`$
The spectrum of spin excitations in the generalized mean field approximation is given by

(83)
$$\omega (k)=N^1\underset{ig}{}J_{ig}A_{ig}\{1\mathrm{exp}[ik(R_iR_j)]\}$$

Now it is not difficult to see that the result (83) includes both the simplest spin-wave dispersion law (39) and the result of Tyablikov decoupling (65) as the limiting cases

(84) $`\omega (k)=<S^z>(J_0J_k)+`$
$`(<2S^z>N)^1{\displaystyle \underset{q}{}}(J_qJ_{kq})(\psi _q^++2\psi _q^{zz})`$
where

$$\psi _q^+=\underset{ij}{}<S_i^{}S_j^+>\mathrm{exp}[iq(R_iR_j)]$$

It is seen that due to the correct definition of generalized mean fields (79) we get the spin excitation spectrum in a general way. In the hydrodynamic limit, it leads to $`\omega (k)k^2`$. The procedure is straightforward, and the details are left as an exercise.
Let us remind that till now no approximation has been made. The expressions (81), (82), and (83) are very useful as the starting point for approximate calculation of the self-energy, a determination of which can only be approximate. To do this, it is first necessary to express, using the spectral theorem (24), the mass operator (81) in terms of correlation functions

(85) $`<2S^z>M_k(\omega )=`$
$`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\omega ^{}}{\omega \omega ^{}}}(\mathrm{exp}(\beta \omega ^{})1){\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑t\mathrm{exp}(i\omega ^{}t)`$
$`N^1{\displaystyle \underset{ijgl}{}}J_{ig}J_{lj}\mathrm{exp}[ik(R_iR_j)]`$
$`{\displaystyle \frac{1}{<2S^z>}}<((S_l^+(t)S_j^z(t)S_j^+(t)S_l^z(t))^{(ir)})^{}|(S_i^+S_g^zS_g^+S_i^z)^{(ir)}>^{(p)}`$
This representation is exact, and only the algebraic properties were used to derive it. Thus, the expression for the analytic structure of the single-particle GF ( or the propagator ) can be deduced without any approximation. A characteristic feature of eq.(81) is that it involves the higher-order GFs. A whole hierarchy of equations involving higher-order GFs could thus be rewritten compactly. Moreover, it not only gives a convenient alternative representation, but avoids some of the algebraic complexities of higher-order Green-function theories. Objective of the present consideration is to give a plausible self-consistent scheme of the approximate calculation of the self-energy within the IGF method. To this end, we should express the higher-order GFs in terms of the initial ones, i.e. find the relevant approximate functional form

$$MF[G]$$

It is clear that this can be done in many ways. As a start, let us consider how to express higher-order correlation function in (85) in terms of the low-order ones. We use the following form

(86) $`<((S_l^+(t)S_j^z(t)S_j^+(t)S_l^z(t))^{(ir)})^{}|(S_i^+S_g^zS_g^+S_i^z)^{(ir)}>^{(p)}`$
$`\psi _{jg}^{zz}(t)\psi _{li}^+(t)\psi _{lg}^{zz}(t)\psi _{ji}^+(t)\psi _{ji}^{zz}(t)\psi _{lg}^+(t)+\psi _{li}^{zz}(t)\psi _{jg}^+(t)`$
We find

(87) $`<2S^z>M_k(\omega )=`$
$`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\omega ^{}}{\omega \omega ^{}}}(\mathrm{exp}(\beta \omega ^{})1){\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑t\mathrm{exp}(i\omega ^{}t)`$
$`N^1{\displaystyle \underset{ijgl}{}}J_{ig}J_{lj}\mathrm{exp}[ik(R_iR_j)]`$
$`{\displaystyle \frac{1}{<2S^z>}}\left(\psi _{jg}^{zz}(t)\psi _{li}^+(t)\psi _{lg}^{zz}(t)\psi _{ji}^+(t)\psi _{ji}^{zz}(t)\psi _{lg}^+(t)+\psi _{li}^{zz}(t)\psi _{jg}^+(t)\right)`$
It is reasonable to approximate the longitudinal correlation function by its static value $`\psi _{ji}^{zz}(t)\psi _{ji}^{zz}(0)`$. The transversal spin correlation functions are given by the expression

(88) $`\psi _{ji}^+(t)=`$
$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{d\omega }{2\pi }}[\mathrm{exp}(\beta \omega )1]^1\mathrm{exp}(i\omega t)(2Im<<S_i^+|S_j^{}>>_{\omega +iϵ})`$
After the substitution of eq.( 88) into eq.( 87) for the self-energy, we find an approximate expression in the self-consistent form, which, together with the exact Dyson equation (82), constitute a self-consistent system of equations for the calculation of the GF. As an example, we start the calculation procedure ( which can be made iterative ) with the simplest first ”trial” expression

$$(2Im<<S_i^+|S_j^{}>>_{\omega +iϵ})\delta (\omega \omega (k))$$

After some algebraic transformations we find

(89)
$$<2S^z>M_k(\omega )N^1\underset{q}{}(J_qJ_{kq})^2(\omega \omega (qk))^1\psi _q^{zz}$$

This expression gives a compact representation for the self-energy of the spin-wave propagator in a Heisenberg ferromagnet. The above calculations show that the inelastic spin-wave scattering effects influence the single-particle spin-wave excitation energy

$$\omega (k,T)=\omega (k)+ReM_k(\omega (k))$$

and the energy width

$$\mathrm{\Gamma }_k(T)=ImM_k(\omega (k))$$

Both these quantities are observable, in principle, via the ferromagnetic resonance or inelastic scattering of neutrons. There is no time to go into details of this aspect of spin-wave interaction effects. It is worthy to note only that it is well known that spin-wave interactions in ferromagnetic insulators have a relatively well-established theoretical foundation, in contrast to the situation with antiferromagnets.

## 7 Heisenberg Antiferromagnet at Finite Temperatures

As it is mentioned above, in this article, we describe the foundation of the IGF method, which is based on the equation-of-motion approach. The strength of this approach lies in its flexibility and applicability to systems with complex spectra and strong interaction. The microscopic theory of the Heisenberg antiferromagnet is of great interest from the point of view of application to any novel many-body technique. This is not only because of the interesting nature of the phenomenon itself but also because of the intrinsic difficulty of solving the problem self-consistently in a wide range of temperatures. In this Section, we briefly describe how the generalized mean fields should be constructed for the case of the Heisenberg antiferromagnet, which become very complicated when one uses other many-body methods, like the diagrammatic technique . Within our IGF scheme, however, the calculations of quasi-particle spectra seem feasible and very compact.

### 7.1 Hamiltonian of the Model

The problem to be considered is the many-body quasi-particle dynamics of the system described by the Hamiltonian

(90)
$$H=\frac{1}{2}\underset{ij}{}\underset{\alpha \alpha ^{}}{}J^{\alpha \alpha ^{}}(ij)\stackrel{}{S}_{i\alpha }\stackrel{}{S}_{j\alpha ^{}}=\frac{1}{2}\underset{q}{}\underset{\alpha \alpha ^{}}{}J_q^{\alpha \alpha ^{}}\stackrel{}{S}_{q\alpha }\stackrel{}{S}_{q\alpha ^{}}$$

This is the Heisenberg-Neel model of an isotropic two-sublattice antiferromagnet (the notation is slightly more general than in Section 4.1.2 ). Here $`S_{i\alpha }`$ is a spin operator situated on site $`i`$ of sublattice $`\alpha `$, and $`J^{\alpha \alpha ^{}}(ij)`$ is the exchange energy between atoms on sites $`R_{i\alpha }`$ and $`R_{j\alpha ^{}}`$; $`\alpha ,\alpha ^{}`$ takes two values $`(a,b)`$ . It is assumed that all of the atoms on sublattice $`\alpha `$ are identical, with spin magnitude $`S_\alpha `$. It should be noted that, in principle, no restrictions are placed in the Hamiltonian (90) on the number of sublattices, or the number of sites on a sublattice. What is important is that sublattices are to be distinguished on the basis of differences in local magnetic characteristics rather than merely differences in geometrical or chemical characteristics.
Let us introduce the spin operators $`S_{i\alpha }^\pm =S_{i\alpha }^x\pm iS_{i\alpha }^y`$. Then the commutation rules for spin operators are

$$[S_{i\alpha }^+,S_{j\alpha ^{}}^{}]_{}=2(S_{i\alpha }^z)\delta _{ij}\delta _{\alpha \alpha ^{}};[S_{i\alpha }^{},S_{j\alpha ^{}}^z]_{}=\pm S_{i\alpha }^{}\delta _{ij}\delta _{\alpha \alpha ^{}}$$

For an antiferromagnet, an exact ground state is not known. Neel  introduced the model concept of two mutually interpenetrating sublattices to explain the behaviour of the susceptibility of antiferromagnets. However, the ground state in the form of two sublattices ( the Neel state ) is only a classical approximation. In contrast to ferromagnets, in which the mean molecular field is approximated relatively reasonably by a function homogeneous and proportional to the magnetisation, in ferri- and antiferromagnets, the mean molecular field is strongly inhomogeneous. The local molecular field of Neel  is a more general concept. Here, we present the calculations  of the quasi-particle spectrum and damping of a Heisenberg antiferromagnet in the framework of the IGF method.
In what follows, it is convenient to rewrite (90) in the form

(91)
$$H=\frac{1}{2}\underset{q}{}\underset{\alpha \alpha ^{}}{}I_q^{\alpha \alpha ^{}}(S_{q\alpha }^+S_{q\alpha ^{}}^{}+S_{q\alpha }^zS_{q\alpha ^{}}^z)$$

where

$$I_q^{\alpha \alpha ^{}}=1/2(J_q^{\alpha \alpha ^{}}+J_q^{\alpha ^{}\alpha })$$

It will be shown that the use of ”anomalous averages” which fix the Neel vacuum makes it possible to determine uniquely generalized mean fields and to calculate, in a very compact manner, the spectrum of spin-wave excitations and their damping due to inelastic magnon-magnon scattering processes. A transformation from the spin operators to Bose (or Pauli ) operators is not required.

### 7.2 Quasi-Particle Dynamics of Heisenberg Antiferromagnet

In this section, to make the discussion more concrete, we consider the retarded GF of localized spins defined as $`G^{AB}(tt^{})=<<A(t),B(t^{})>>`$ . Our attention is focused on the spin dynamics of the model. To describe the spin dynamics of the model ( 91) self-consistently, one should take into account the full algebra of relevant operators of the suitable ”spin modes” ( ”relevant degrees of freedom” ) which are appropriate for the case. This relevant algebra should be described by the ’spinor’ $`A=\left(\genfrac{}{}{0pt}{}{S_{ka}^+}{S_{kb}^+}\right)`$, $`B=A^{}`$ , according to the IGF strategy of Section 3.
Once this has been done, we must introduce the generalized matrix GF of the form

(92)
$$\left(\begin{array}{cc}<<S_{ka}^+|S_{ka}^{}>>& <<S_{ka}^+|S_{kb}^{}>>\\ <<S_{kb}^+|S_{ka}^{}>>& <<S_{kb}^+|S_{kb}^{}>>\end{array}\right)=\widehat{G}(k;\omega )$$

To show the advantages of the IGF in the most full form, we carry out the calculations in the matrix form.
To demonstrate the utility of the IGF method, we consider the following steps in a more detailed form. Differentiating the GF $`<<S_{ka}^+|B>>`$ with respect to the first time, $`t`$, we find

(93) $`\omega <<S_{ka}^+|\left({\displaystyle \genfrac{}{}{0pt}{}{S_{ka}^{}}{S_{kb}^{}}}\right)>>_\omega =`$
$`\left\{{\displaystyle \genfrac{}{}{0pt}{}{2<S_a^z>}{0}}\right\}+{\displaystyle \frac{1}{N^{1/2}}}{\displaystyle \underset{q}{}}I_q^{ab}<<S_{kq}^{ab}|B_{ab}>>_\omega `$
$`+{\displaystyle \frac{1}{N^{1/2}}}{\displaystyle \underset{q}{}}I_q^{aa}<<S_{kq}^{aa}|B_{ab}>>_\omega `$
where $`S_{kq}^{ab}=(S_{kq,a}^+S_{qb}^zS_{qb}^+S_{kq,a}^z)`$.
In (93), we introduced the notation

$$B_{ab}=\left\{\genfrac{}{}{0pt}{}{S_{ka}^{}}{S_{kb}^{}}\right\};B_{ba}=\left\{\genfrac{}{}{0pt}{}{S_{kb}^{}}{S_{ka}^{}}\right\}$$

Let us define the irreducible $`(ir)`$ operators as (equivalently, it is possible to define the irreducible GFs)

(94) $`(S_{kq}^{ab}){}_{}{}^{(ir)}=S_{kq}^{ab}A_q^{ab}S_{ka}^++A_{kq}^{ba}S_{kb}^+`$
(95) $`(S_{q\alpha }^z)^{(ir)}=S_{q\alpha }^zN^{1/2}<S_\alpha ^z>\delta _{q,0}`$
The choice of the irreducible parts is uniquely determined by the ”orthogonality” constraint ( 31)

(96)
$$<\left[(S_{kq}^{ab}){}_{}{}^{(ir)},\left(\genfrac{}{}{0pt}{}{S_{ka}^{}}{S_{kb}^{}}\right)\right]>=0$$

From eq.(96) we find that

(97)
$$A_q^{ab}=\frac{2<(S_{qa}^z)^{(ir)}(S_{qb}^z)^{(ir)}>+<S_{qa}^{}S_{qb}^+>}{2N^{1/2}<S_a^z>}$$

By using the definition of the irreducible parts (94), the equation of motion (93) can be exactly transformed to the following form

(98) $`(\omega \omega _{aa})<<S_{ka}^+|B_{ab}>>_\omega +\omega _{ab}<<S_{kb}^+|B_{ab}>>_\omega =`$
$`\left\{{\displaystyle \genfrac{}{}{0pt}{}{2<S_a^z>}{0}}\right\}+<<\mathrm{\Phi }_a^{(ir)}(k)|B_{ab}>>_\omega `$
(99) $`(\omega \omega _{bb})<<S_{kb}^+|B_{ba}>>_\omega +\omega _{ba}<<S_{ka}^+|B_{ba}>>_\omega =`$
$`\left\{{\displaystyle \genfrac{}{}{0pt}{}{2<S_b^z>}{0}}\right\}+<<\mathrm{\Phi }_b^{(ir)}(k)|B_{ba}>>_\omega `$
The following notation was used:

(100) $`\omega _{aa}=((I_0^{aa}I_k^{aa})<S_a^z>+I_0^{ab}<S_b^z>+`$
$`{\displaystyle \underset{q}{}}[(I_q^{aa}I_{kq}^{aa})A_{Nq}^{aa}+I_q^{ab}A_{Nq}^{ab}])`$
(101) $`\omega _{ab}=\left(I_k^{ab}<S_a^z>+{\displaystyle \underset{q}{}}I_{kq}^{ab}A_{Nq}^{ba}\right)`$
(102) $`A_{Nq}^{\alpha \beta }=N^{1/2}A_q^{\alpha \beta }`$
(103) $`\mathrm{\Phi }_a^{(ir)}(k)=`$
$`N^{1/2}{\displaystyle \underset{q}{}}{\displaystyle \underset{\gamma =a,b}{}}I_q^{\alpha \gamma }[S_{kq,a}^+(S_{q\gamma }^z)^{(ir)}S_{q\gamma }^+(S_{kq,a}^z)^{(ir)}]^{(ir)}`$
To calculate the irreducible GFs on the right-hand sides of eqs. (98) and (99), we use the device of differentiating with respect to the second time $`t^{}`$. After introduction of the corresponding irreducible parts into the resulting equations, the system of equations can be represented in the matrix form which can be identically transformed to the standard form (32)

(104)
$$\widehat{G}(k,\omega )=\widehat{G}_0(k,\omega )+\widehat{G}_0(k,\omega )\widehat{P}(k,\omega )\widehat{G}_0(k,\omega )$$

Here we introduced the generalized mean-field (GMF) GF $`G_0`$ and the scattering operator $`P`$ according to the following definitions

(105)
$$\widehat{G}_0=\widehat{\mathrm{\Omega }}^1\widehat{I}$$

(106) $`\widehat{P}=`$
$`{\displaystyle \frac{1}{4<S_a^z>^2}}\left(\begin{array}{cc}<<\mathrm{\Phi }_a^{(ir)}(k)|\mathrm{\Phi }_a^{(ir)}(k)>>& <<\mathrm{\Phi }_a^{(ir)}(k)|\mathrm{\Phi }_b^{(ir)}(k)>>\\ <<\mathrm{\Phi }_b^{(ir)}(k)|\mathrm{\Phi }_a^{(ir)}(k)>>& <<\mathrm{\Phi }_b^{(ir)}(k)|\mathrm{\Phi }_b^{(ir)}(k)>>\end{array}\right)`$
where

(107)
$$\widehat{\mathrm{\Omega }}=\left(\begin{array}{cc}(\omega \omega _{aa})& \omega _{ab}\\ \omega _{ab}& (\omega \omega _{bb})\end{array}\right)$$

The Dyson equation can be written exactly in the form (35) where the mass operator $`M`$ is of the form

(108)
$$\widehat{M}(k,\omega )=(\widehat{P}(k,\omega ))^{(p)}$$

It follows from the Dyson equation that

$$\widehat{P}(k,\omega )=\widehat{M}(k,\omega )+\widehat{M}(k,\omega )\widehat{G}_0(k,\omega )\widehat{P}(k,\omega )$$

Thus, on the basis of these relations, we can speak of the mass operator $`M`$ as the proper part of the operator $`P`$ by analogy with the diagram technique, in which the mass operator is the connected part of the scattering operator. As it is shown in Section 3, the formal solution of the Dyson equation is of the form (38). Hence, the determination of the full GF $`\widehat{G}`$ was reduced to the determination of $`\widehat{G}_0`$ and $`\widehat{M}`$.

### 7.3 Generalized Mean-Field GF

From the definition (105), the GF matrix in the generalized mean-field approximation reads

(109) $`\widehat{G}_0=`$
$`\left(\begin{array}{cc}G_0^{aa}(k,\omega )& G_0^{ab}(k,\omega )\\ G_0^{ba}(k,\omega )& G_0^{bb}(k,\omega )\end{array}\right)={\displaystyle \frac{2<S_a^z>}{det\widehat{\mathrm{\Omega }}}}\left(\begin{array}{cc}(\omega \omega _{aa})& \omega _{ab}\\ \omega _{ab}& (\omega \omega _{bb})\end{array}\right)`$
where

$$det\widehat{\mathrm{\Omega }}=(\omega \omega _{aa})(\omega \omega _{bb})\omega _{aa}\omega _{ab}$$

We find the poles of GF (109) from the equation

$$det\widehat{\mathrm{\Omega }}=0$$

from which it follows that

(110)
$$\omega _\pm (k)=\pm \sqrt{(\omega _{aa}^2(k)\omega _{ab}^2(k))}$$

It is convenient to adopt here the Bogoliubov $`(u,v)`$-transformation notation by analogy with that of Section 4.1.2. The elements of the matrix GF $`G_0(k,\omega )`$ are found to be

(111)
$$G_0^{aa}(k,\omega )=2<S_a^z>\left[\frac{u^2(k)}{\omega \omega _+(k)}\frac{v^2(k)}{\omega \omega _{}(k)}\right]=G_0^{bb}(k,\omega )$$

(112)
$$G_0^{ab}(k,\omega )=2<S_a^z>\left[\frac{u(k)v(k)}{\omega \omega _+(k)}+\frac{u(k)v(k)}{\omega \omega _{}(k)}\right]=G_0^{ba}(k,\omega )$$

where

$`u^2(k)=1/2[(1\gamma _k^2)^{1/2}+1];v^2(k)=1/2[(1\gamma _k^2)^{1/2}1]`$
(113) $`\gamma _k={\displaystyle \frac{1}{z}}{\displaystyle \underset{i}{}}\mathrm{exp}(ikR_i);I_q^{aa}=I_q^{bb}=0`$
The simplest assumption is that each sublattice is s.c. and $`\omega _{\alpha \alpha }(k)=0(\alpha =a,b)`$. Although that we work in the GFs formalism, our expressions (111), (112) are in accordance with the results of the Bogoliubov (u,v)-transformation , but, of course, the present derivation is more general. However, it is possible to say that we diagonalized the generalized mean-field GF by introducing a new set of operators. We used the notation

(114)
$$S_1^+(k)=u_kS_{ka}^++v_kS_{kb}^+;S_2^+(k)=v_kS_{ka}^++u_kS_{kb}^+$$

This notation permits us to write down the results in a compact and convenient form, but all calculations can be done in the initial notation too.
The spectrum of elementary excitations in the GMF approximation for an arbitrary spin $`S`$ is of the form

(115)
$$\omega (k)=Iz<S_a^z>\left[1\frac{1}{N^{1/2}<S_a^z>}\underset{q}{}\gamma _qA_q^{ab}\right]\sqrt{(1\gamma _k^2)}$$

where $`I_q=zI\gamma _q`$, and $`z`$ is the number of nearest neighbors in the lattice. The first term in (115) corresponds to the Tyablikov approximation ( cf.(46)). The second term in (115) describes the elastic scattering of the spin-wave quasi-particles. At low temperatures, the fluctuations of the longitudinal spin components are small, and, therefore, for (115) we obtain

(116)
$$\omega (k)ISz[1C(T)]\sqrt{(1\gamma _k^2)}$$

The function $`C(T)`$ determines the temperature dependence of the spin-wave spectrum

$$C(T)=\frac{1}{2NS^2}\underset{q}{}(<S_{qa}^{}S_{qa}^+>+\gamma _q<S_{qa}^{}S_{qb}^+>)$$

In the case when $`C(T)0`$, we obtain the result of the Tyablikov decoupling for the spectrum of the antiferromagnons

(117)
$$\omega (k)I<S_a^z>z\sqrt{(1\gamma _k^2)}$$

In the hydrodynamic limit, when $`\omega (k)D(T)|\stackrel{}{k}|`$, we can conclude that the stiffness constant $`D(T)=zIS(1C(T))`$ for an antiferromagnet decreases with temperature because of the elastic magnon-magnon scattering as $`T^4`$. To estimate the contribution of the inelastic scattering processes, it is necessary to take into account the corrections due to the mass operator.

### 7.4 Damping of Quasi-Particle Excitations

An antiferromagnet is a system with a complicated quasi-particle spectrum. The calculation of the damping due to inelastic scattering processes in a system of that sort has some important aspects. When calculating the damping, it is necessary to take into account the contributions from all matrix elements of the mass operator $`M`$

$$M=G_0^1G^1$$

It is then convenient to use the representation in which the generalized mean field GF has a diagonal form . In terms of the new operators $`S_1`$ and $`S_2`$, the GF $`G`$ takes the form

$$\stackrel{~}{G}(k;\omega )=\left(\begin{array}{cc}<<\stackrel{~}{S}_1^+(k)|\stackrel{~}{S}_1^{}(k)>>& <<\stackrel{~}{S}_1^+(k)|\stackrel{~}{S}_2^{}(k)>>\\ <<\stackrel{~}{S}_2^+(k)|\stackrel{~}{S}_1^{}(k)>>& <<\stackrel{~}{S}_2^+(k)|\stackrel{~}{S}_2^{}(k)>>\end{array}\right)=\left(\begin{array}{cc}G_{11}& G_{12}\\ G_{21}& G_{22}\end{array}\right)$$

In other words, the damping of the quasi-particle excitations is determined on the basis of a GF of the form

(118)
$$G_{11}(k,\omega )=\frac{2<S_a^z>}{\omega \omega (k)2<S_a^z>\mathrm{\Sigma }(k,\omega )}$$

Here, the self-energy operator $`\mathrm{\Sigma }(k,\omega )`$ is determined by the expression

(119)
$$\mathrm{\Sigma }(k,\omega )=M_{11}(k,\omega )\frac{2<S_a^z>M_{12}(k,\omega )M_{21}(k,\omega )}{\omega +\omega (k)2<S_a^z>M_{22}(k,\omega )}$$

In the case when $`k,\omega 0`$, one can be restricted to the approximation

(120)
$$\mathrm{\Sigma }(k,\omega )M_{11}(k,\omega )=u_k^2M_{aa}+v_ku_k(M_{ab}+M_{ba})+v_k^2M_{bb}$$

It follows from (108) that to calculate the damping, it is necessary to find the GFs $`<<\mathrm{\Phi }_\alpha ^{(ir)}(k)|\mathrm{\Phi }_\beta ^{(ir)}(k)>>`$. As an example, we consider the calculation of one of them. By means of the spectral theorem (25), we can express the GF in terms of the correlation function $`<\mathrm{\Phi }_a^{(ir)}(k)\mathrm{\Phi }_a^{(ir)}(k,t)>`$. We have

(121) $`<<\mathrm{\Phi }_a^{(ir)}(k)|\mathrm{\Phi }_a^{(ir)}(k)>>=`$
$`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\omega ^{}}{\omega \omega ^{}}}(\mathrm{exp}(\beta \omega ^{})1){\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑t\mathrm{exp}(i\omega ^{}t)<\mathrm{\Phi }_a^{(ir)}(k)\mathrm{\Phi }_a^{(ir)}(k,t)>`$
Thus, it is necessary to find a workable ”trial” approximation for the correlation function on the r.h.s. of (121). We consider an approximation of the following form

(122) $`<(S_{qb}^z)^{(ir)}S_{(kq^{})a}^{}S_{(kq^{})a}^+(t)(S_{q^{}b}^z(t))^{(ir)}>`$
$`{\displaystyle \frac{1}{4NS^2}}{\displaystyle \underset{p}{}}\left(\psi _{kp,aa}^+(t)\psi _{q+p,bb}^+(t)\psi _{p,bb}^+(t)+\psi _{kq,ab}^+(t)\psi _{q+p,ab}^+(t)\psi _{p,bb}^+(t)\right)\delta _{q,q^{}}`$
where $`\psi _{q,ab}^+(t)=<S_{qa}^{}S_{qb}^+(t)>`$. By analogy with the diagram technique, we can say that the approximation (122) corresponds to the neglect of the vertex corrections to the magnon-magnon inelastic collisions. Using (122) in (121), we obtain

(123) $`<<\mathrm{\Phi }_a^{(ir)}(k)|\mathrm{\Phi }_a^{(ir)}(k)>>`$
$`{\displaystyle \frac{1}{16NS^4}}{\displaystyle \underset{qp}{}}{\displaystyle \frac{d\omega _1d\omega _2d\omega _3}{\omega \omega _1\omega _2+\omega _3}F(\omega _1,\omega _2,\omega _3)}`$
$`[{\displaystyle \frac{1}{\pi }}ImG_{aa}(kq,\omega _1)][{\displaystyle \frac{1}{\pi }}ImG_{bb}(q+p,\omega _2)][{\displaystyle \frac{1}{\pi }}ImG_{bb}(p,\omega _3)]`$
where

(124)
$$F(\omega _1,\omega _2,\omega _3)=N(\omega _2)[N(\omega _3)N(\omega _1)]+[1+N(\omega _1)]N(\omega _3)$$

Equations (35), (108), and (123) constitute a self-consistent system of equations. To solve this system of equations, we can, in principle, use any convenient initial representation for the GF, substituting it into the right-hand side of eq. (123). The system can then be solved iteratively. To estimate the damping, it is usually sufficient, as the first iteration, to use the simplest single-pole approximation

(125)
$$\frac{1}{\pi }ImG(k,\omega )\delta (\omega \omega (k))$$

As a result, for the damping of the spin-wave excitations we obtain

(126) $`\mathrm{\Gamma }(k,\omega )=2SIm\mathrm{\Sigma }(k,\omega )=`$
$`={\displaystyle \frac{\pi }{N}}(zI)^2(1e^{(\beta \omega )})`$
$`{\displaystyle \underset{qp}{}}N_p(1+N_{q+p})(1+N_{kq})M_{11}(k,p;kq,p+q)\delta (\omega \omega (kq)+\omega (p))`$
The explicit expression for $`M_{11}`$ is given in ref. . In our approach, it is possible to take into account the inelastic scattering of spin waves due to scattering by the longitudinal spin fluctuations too . In general, the correct estimates of the temperature dependence of the damping of antiferromagnons depend strongly on the reduced temperature and energy scales and are rather a nontrivial task. However, under the normal conditions, the damping is weak $`\omega (k)/\mathrm{\Gamma }10^210^3`$, and the antiferromagnons are the well-defined quasi-particle excitations.
In summary, in this Section, we have shown that the IGF method permits us to calculate the spectrum and the damping for a two-sublattice Heisenberg antiferromagnet in a wide range of temperatures in a compact and self-consistent way. At the same time, a certain advantage is that all the calculation can be made in the representation of spin operators for an arbitrary spin $`S`$. The theory we have developed can be directly extended to the case of a large number of magnetic sublattices with inequivalent spins, i.e., it can be used to describe the complex ferrimagnets.
In the framework of our IGF approach, it was shown that the mean fields in an antiferromagnet must include the ”anomalous” averages which represent the local nature of the Neel molecular fields. Thus, the mean field in an antiferromagnet, like the mean field in a superconductor, has a more complicated structure.

## 8 Quasi-Particle Dynamics of Lattice Fermion Models

### 8.1 Hubbard Model. Weak Correlation

The concept of GMFs and the relevant algebra of operators from which GFs are constructed are important for our treatment of electron correlations in solids. It is convenient (and much shorter) to discuss these concepts for weakly and strongly correlated cases separately. First, we should construct a suitable state vector space of a many-body system . The fundamental assumption implies that states of a system of interacting particles can be expanded in terms of states of non-interacting particles . This approach originates from perturbation theory and finds support for weakly interacting many-particle systems. For the strongly correlated case, this approach needs a suitable reformulation, and just at this point, the right definition of the GMFs is vital. Let us consider the weakly correlated Hubbard model (47). In some respect, this case is similar to the ordinary interacting electron gas but with very local singular interaction. The difference is in the lattice ( Wannier ) character of electron states. It is shown below that the usual creation $`a_{i\sigma }^{}`$ and annihilation $`a_{i\sigma }`$ second-quantized operators with the properties

$$a_i^{}\mathrm{\Psi }^{(0)}=\mathrm{\Psi }_i^{(1)};a_i\mathrm{\Psi }^{(1)}=\mathrm{\Psi }^{(0)}$$

$$a_i\mathrm{\Psi }^{(0)}=0;a_j\mathrm{\Psi }_i^{(1)}=0(ij)$$

are suitable variables for description of a system under consideration. Here $`\mathrm{\Psi }^{(0)}`$ and $`\mathrm{\Psi }^{(1)}`$ are vacuum and single-particle states, respectively. The question now is how to describe our system in terms of quasi-particles. For a translationally invariant system, to describe the low-lying excitations of a system in terms of quasi-particles , one has to choose eigenstates such that they all correspond to a definite momentum. For the single-band Hubbard model (47), the exact transformation reads

$$a_{\stackrel{}{k}\sigma }=N^{1/2}\underset{i}{}\mathrm{exp}(i\stackrel{}{k}\stackrel{}{R}_i)a_{i\sigma }$$

Note that for a degenerate band model, a more general transformation is necessary. Then the Hubbard Hamiltonian (47) in the Bloch vector state space is given by

(127)
$$H=\underset{k\sigma }{}ϵ(k)a_{k\sigma }^{}a_{k\sigma }+U/2N\underset{pqrs}{}a_{p+rq\sigma }^{}a_{p\sigma }a_{q\sigma }^{}a_{r\sigma }$$

If the interaction is weak, the algebra of relevant operators is very simple: it is an algebra of a non-interacting fermion system ($`a_{k\sigma },a_{k\sigma }^{},n_{k\sigma }=a_{k\sigma }^{}a_{k\sigma }`$). To calculate of the electron quasi-particle spectrum of the Hubbard model in this limit, let us consider the single-electron GF defined as

(128)
$$G_{k\sigma }(tt^{})=<<a_{k\sigma },a_{k\sigma }^{}>>=i\theta (tt^{})<[a_{k\sigma }(t),a_{k\sigma }^{}(t^{})]_+>$$

The equation of motion for the Fourier transform of GF $`G_{k\sigma }(\omega )`$ is of the form

(129)
$$(\omega ϵ_k)G_{k\sigma }(\omega )=1+U/N\underset{pq}{}<<a_{k+p\sigma }a_{p+q\sigma }^{}a_{q\sigma }|a_{k\sigma }^{}>>_\omega $$

Let us introduce an “irreducible” GF in the following way

$`{}_{}{}^{(ir)}<<a_{k+p\sigma }a_{p+q\sigma }^{}a_{q\sigma }|a_{k\sigma }^{}>>_\omega =`$
(130) $`<<a_{k+p\sigma }a_{p+q\sigma }^{}a_{q\sigma }|a_{k\sigma }^{}>>_\omega \delta _{p,0}<n_{q\sigma }>G_{k\sigma }`$
The irreducible $`(ir)`$ GF in (8.1) is defined so that it cannot be reduced to GF of lower order with respect to the number of fermion operators by an arbitrary pairing of operators or, in other words, by any kind of decoupling. Substituting (8.1) into (129), we obtain

(131) $`G_{k\sigma }(\omega )=G_{k\sigma }^{MF}(\omega )+`$
$`G_{k\sigma }^{MF}(\omega )U/N{\displaystyle \underset{pq}{}}{}_{}{}^{(ir)}<<a_{k+p\sigma }a_{p+q\sigma }^{}a_{q\sigma }|a_{k\sigma }^{}>>_\omega `$
Here we introduced the notation

(132)
$$G_{k\sigma }^{MF}(\omega )=(\omega ϵ(k\sigma ))^1;ϵ(k\sigma )=ϵ(k)+U/N\underset{q}{}<n_{q\sigma }>$$

In this paper, for brevity, we confine ourselves to considering the paramagnetic solutions, i.e. $`<n_\sigma >=<n_\sigma >`$. To calculate the higher-order GF on the r.h.s. of (131), we have to write the equation of motion obtained by means of differentiation with respect to the second variable $`t^{}`$. Constraint (29) allows us to remove the inhomogeneous term from this equation for $`\frac{d}{dt^{}}^{(ir)}<<A(t),a_{k\sigma }^{}(t^{})>>`$.
For the Fourier components, we have

$`(\omega ϵ(k)){}_{}{}^{(ir)}<<A|a_{k\sigma }^{}>>_\omega =<^{(ir)}[A,a_{k\sigma }^{}]_+>+`$
(133) $`U/N{\displaystyle \underset{rs}{}}{}_{}{}^{(ir)}<<A|a_{r\sigma }^{}a_{r+s\sigma }a_{k+s\sigma }^{}>>_\omega .`$
The anticommutator in (8.1) is calculated on the basis of the definition of the irreducible part

$`<[^{(ir)}(a_{k+p\sigma }a_{p+q\sigma }^{}a_{q\sigma }),a_{k\sigma }^{}]_+>=`$
(134) $`<[a_{k+p\sigma }a_{p+q\sigma }^{}a_{q\sigma }<a_{p+q\sigma }^{}a_{q\sigma }>a_{k+p\sigma },a_{k\sigma }^{}]_+>=0`$
If one introduces the irreducible part for the r.h.s. operators by analogy with expression (8.1), the equation of motion (129) takes the following exact form ( cf. eq.(32))

(135)
$$G_{k\sigma }(\omega )=G_{k\sigma }^{MF}(\omega )+G_{k\sigma }^{MF}(\omega )P_{k\sigma }(\omega )G_{k\sigma }^{MF}(\omega )$$

where we introduced the following notation for the operator $`P`$ (35)

(136) $`P_{k\sigma }(\omega )={\displaystyle \frac{U^2}{N^2}}{\displaystyle \underset{pqrs}{}}D_{k\sigma }^{(ir)}(p,q|r,s,;\omega )=`$
$`{\displaystyle \frac{U^2}{N^2}}{\displaystyle \underset{pqrs}{}}(^{(ir)}<<a_{k+p\sigma }a_{p+q\sigma }^{}a_{q\sigma }|a_{r\sigma }^{}a_{r+s\sigma }a_{k+s\sigma }^{}>>_\omega ^{(ir)})`$
To define the self-energy operator according to (34), one should separate the ”proper” part in the following way

$`D_{k\sigma }^{(ir)}(p,q|r,s;\omega )=L_{k\sigma }^{(ir)}(p,q|r,s;\omega )`$
(137) $`+{\displaystyle \frac{U^2}{N^2}}{\displaystyle \underset{r^{}s^{}p^{}q^{}}{}}L_{k\sigma }^{(ir)}(p,q|r^{}s^{};\omega )G_{k\sigma }^{MF}(\omega )D_{k\sigma }^{(ir)}(p^{},q^{}|r,s;\omega )`$
Here $`L_{k\sigma }^{(ir)}(p,q|r,s;\omega )`$ is the “proper” part of GF $`D_{k\sigma }^{(ir)}(p,q|r,s;\omega )`$ which, in accordance with the definition (36), cannot be reduced to the lower-order one by any type of decoupling. We find

(138)
$$G_{k\sigma }=G_{k\sigma }^{MF}(\omega )+G_{k\sigma }^{MF}(\omega )M_{k\sigma }(\omega )G_{k,\sigma }(\omega )$$

Equation (138) is the Dyson equation for the single-particle double-time thermal GF. According to (38), it has the formal solution

(139)
$$G_{k\sigma }(\omega )=[\omega ϵ(k\sigma )M_{k\sigma }(\omega )]^1$$

where the self-energy operator $`M`$ is given by

(140) $`M_{k\sigma }(\omega )={\displaystyle \frac{U^2}{N^2}}{\displaystyle \underset{pqrs}{}}L_{k\sigma }^{(ir)}(p,q|r,s;\omega )=`$
$`{\displaystyle \frac{U^2}{N^2}}{\displaystyle \underset{pqrs}{}}(^{(ir)}<<a_{k+p\sigma }a_{p+q\sigma }^{}a_{q\sigma }|a_{r\sigma }^{}a_{r+s\sigma }a_{k+s\sigma }^{}>>^{(ir)})^{(p)}`$
We wrote explicitly equation (136) for $`P`$ and equation (140) for $`M`$ to illustrate the general arguments of Section 3 and to give concrete equations for determining both the quantities, $`P`$ and $`M`$.
The latter expression (140) is an exact representation ( no decoupling was made till now ) for the self-energy in terms of higher-order GF up to second order in $`U`$ ( for the consideration of higher-order equations of motion, see ref. ). The explicit difference between $`P`$ and $`M`$ follows from the functional form (36). Thus, in contrast to the standard equation-of-motion approach, the calculation of full GF was substituted by the calculation of the mean-field GF $`G^{MF}`$ and the self-energy operator $`M`$. The main reason for this method of calculation is that the decoupling is only introduced into the self-energy operator, as it will be shown in a detailed form below. The formal solution of the Dyson equation (36) determines the right reference frame for the formation of the quasi-particle spectrum due to its own correct functional structure. In the standard equation-of-motion approach, that structure could be lost by using decoupling approximations before arriving at the correct functional structure of the formal solution of the Dyson equation. This is a crucial point of the IGF method.
The energies of electron states in the mean-field approximation are given by the poles of $`G^{MF}`$ . Now let us consider the damping effects and finite lifetimes. To find an explicit expression for the self-energy $`M`$ (140), we have to evaluate approximately the higher-order GF in it. It will be shown below that the IGF method permits one to derive the damping in a self-consistent way simply and much more generally than within other formulations. First, it is convenient to write down the GF in (140) in terms of correlation functions by using the spectral theorem (26)

$`<<a_{k+p\sigma }a_{p+q\sigma }^{}a_{q\sigma }|a_{k+s\sigma }^{}a_{r\sigma }^{}a_{r+s\sigma }>>_\omega =`$
$`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\omega ^{}}{\omega \omega ^{}}}(\mathrm{exp}(\beta \omega ^{})+1){\displaystyle _{\mathrm{}}^+\mathrm{}}\mathrm{exp}(i\omega ^{}t)`$
(141) $`<a_{k+s\sigma }^{}(t)a_{r\sigma }^{}(t)a_{r+s\sigma }(t)a_{k+p\sigma }a_{p+q\sigma }^{}a_{q\sigma }>`$
Further insight is gained if we select the suitable relevant “trial” approximation for the correlation function on the r.h.s. of (8.1). In this paper, we show that the earlier formulations based on the decoupling or/and diagrammatic methods can be obtained from our technique but in a self-consistent way. It is clear that a relevant trial approximation for the correlation function in (8.1) can be chosen in many ways. For example, the reasonable and workable one can be the following “pair approximation” that is especially suitable for a low density of quasi-particles:

$`<a_{k+s\sigma }^{}(t)a_{r\sigma }^{}(t)a_{r+s\sigma }(t)a_{k+p\sigma }a_{p+q\sigma }^{}a_{q\sigma }>^{(ir)}`$
$`<a_{k+p\sigma }^{}(t)a_{k+p\sigma }><a_{q\sigma }^{}(t)a_{q\sigma }><a_{p+q\sigma }(t)a_{p+q\sigma }^{}>`$
(142) $`\delta _{k+s,k+p}\delta _{r,q}\delta _{r+s,p+q}`$
Using (8.1) and (8.1) in (140) we obtain the self-consistent approximate expression for the self-energy operator ( the self-consistency means that we express approximately the self-energy operator in terms of the initial GF, and, in principle, one can obtain the required solution by a suitable iteration procedure ):

(143) $`M_{k\sigma }(\omega )=`$
$`{\displaystyle \frac{U^2}{N^2}}{\displaystyle \underset{pq}{}}{\displaystyle \frac{d\omega _1d\omega _2d\omega _3}{\omega +\omega _1\omega _2\omega _3}}`$
$`\left[n(\omega _2)n(\omega _3)+n(\omega _1)\left(1n(\omega _2)n(\omega _3)\right)\right]g_{p+q\sigma }(\omega _1)g_{k+p\sigma }(\omega _2)g_{q\sigma }(\omega _3)`$
where we used the notation

$$g_{k\sigma }(\omega )=\frac{1}{\pi }ImG_{k\sigma }(\omega +i\epsilon );n(\omega )=[\mathrm{exp}(\beta \omega )+1]^1$$

Equations (143) and (138) constitute a closed self-consistent system of equations for the single-electron GF of the Hubbard model in the weakly correlated limit. In principle, we can use, on the r.h.s. of (143), any workable first iteration-step form of the GF and find a solution by iteration ( see Appendix D ). It is most convenient to choose, as the first iteration step, the following simple one-pole approximation:

(144)
$$g_{k\sigma }(\omega )\delta (\omega ϵ(k\sigma ))$$

Then, using (144) in (143), we get, for the self-energy, the explicit and compact expression

(145)
$$M_{k\sigma }(\omega )=\frac{U^2}{N^2}\underset{pq}{}\frac{n_{p+q\sigma }(1n_{k+p\sigma }n_{q\sigma })+n_{k+p\sigma }n_{q\sigma }}{\omega +ϵ(p+q\sigma )ϵ(k+p\sigma )ϵ(q\sigma )}$$

Formula (145) for the self-energy operator shows the role of correlation effects ( inelastic scattering processes ) in the formation of quasi-particle spectrum of the Hubbard model. This formula can be derived by several different methods, including perturbation theory. Here we derived it from our IGF formalism as a known limiting case. The numerical calculations of the typical behaviour of real and imaginary parts of the self-energy (145) were performed , for the model density of states of the FCC lattice. These calculations and many other (see e.g.,,) show clearly that the conventional one-electron approximation of the band theory is not always a sufficiently good approximation for transition metals like nickel. A more concrete discussion of the numerical calculations and their comparison with experiments deserve a separate consideration and will be considered elsewhere (for a detailed recent discussion, see ).
Although the solution deduced above is a good evidence for the efficiency of the IGF formalism, there is one more stringent test of the method that we can perform. It is instructive to examine other types of possible trial solutions for the six-operator correlation function in the eq.(8.1). The approximation we propose now reflects the interference between the one-particle branch of the spectrum and the collective ones:

$`<a_{k+s\sigma }^{}(t)a_{r\sigma }^{}(t)a_{r+s\sigma }(t)a_{k+p\sigma }a_{p+q\sigma }^{}a_{q\sigma }>^{(ir)}`$
$`<a_{k+s\sigma }^{}(t)a_{k+p\sigma }><a_{r\sigma }^{}(t)a_{r+s\sigma }(t)a_{p+q\sigma }^{}a_{q\sigma }>+`$
$`<a_{r+s\sigma }(t)a_{p+q\sigma }^{}><a_{k+s\sigma }^{}(t)a_{r\sigma }^{}(t)a_{k+p\sigma }a_{q\sigma }>+`$
(146) $`<a_{r\sigma }^{}(t)a_{q\sigma }><a_{k+s\sigma }^{}(t)a_{r+s\sigma }(t)a_{k+p\sigma }a_{p+q\sigma }^{}>`$
It is seen that the three contributions in this trial solution describe the self-energy corrections that take into account the collective motions of electron density, the spin density and the density of “doubles”, respectively. An essential feature of this approximation is that a correct calculation of the single-electron quasi-particle spectra with damping requires a suitable incorporation of the influence of collective degrees of freedom on the single-particle ones. The most interesting contribution comes from spin degrees of freedom, since the correlated systems are often magnetic or have very well developed magnetic fluctuations.
We follow the above steps and calculate the self-energy operator (140) as

$`M_{k\sigma }(\omega )={\displaystyle \frac{U^2}{N}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\omega _1𝑑\omega _2{\displaystyle \frac{1+N(\omega _1)n(\omega _2)}{\omega \omega _1\omega _2}}`$
$`{\displaystyle \underset{i,j}{}}\mathrm{exp}[i\stackrel{}{k}(\stackrel{}{R}_i\stackrel{}{R}_j)]({\displaystyle \frac{1}{\pi }}Im<<S_i^\pm |S_j^{}>>_{\omega _1})`$
(147) $`({\displaystyle \frac{1}{\pi }}Im<<a_{i\sigma }|a_{j\sigma }^{}>>_{\omega _2})`$
where the following notation was used:

$$S_i^+=a_i^{}a_i;S_i^{}=a_i^{}a_i$$

It is possible to rewrite ( 8.1 ) in a more convenient way

(148) $`M_{k\sigma }(\omega )=`$
$`{\displaystyle \frac{U^2}{N}}{\displaystyle \underset{q}{}}{\displaystyle 𝑑\omega ^{}(\mathrm{cot}\frac{\omega \omega ^{}}{2T}+\mathrm{tan}\frac{\omega ^{}}{2T})}`$
$`\left({\displaystyle \frac{1}{\pi }}Im\chi ^\pm (kq,\omega \omega ^{})g_{q\sigma }(\omega ^{})\right)`$
Equations ( 148 ) and ( 138 ) constitute again another self-consistent system of equations for the single-particle GF of the Hubbard model. Note that both the expressions for the self-energy depend on the quasi-momentum; in other words, the approximate procedure does not break the momentum conservation law. The fundamental importance of equations (148) and (143) can be appreciated by examining the problem of the definition of the Fermi surface. It is rather clear, because the poles $`\omega (k,\sigma )=ϵ(k,\sigma )i\mathrm{\Gamma }_k`$ of GF (139) are determined by the equation

$$\omega ϵ(k\sigma )Re[M_{k\sigma }(\omega )]=0$$

It can be shown quite generally that the Luttinger’s definition of the true Fermi surface  is valid in the framework of the present theory. It is worthy to note that for electrons in a crystal where there is a band index, and a quasi-momentum, the definition of the Fermi surface is a little more complicated than the single-band one. Before the single particle energies and Fermi surface are known, one should carry out a diagonalization in the band index.

### 8.2 Hubbard Model. Strong Correlation

Being convinced that the IGF method can be applied successfully to the weakly correlated Hubbard model, we now show that the IFG approach can be extended to the case of an arbitrarily strong but finite interaction. This development incorporates main advantages of the IGF scheme and proves its efficiency and flexibility.
When studying the electron quasi-particle spectrum of strongly correlated systems, one should take care of at least three facts of major importance:

* The ground state is reconstructed radically as compared with the weakly correlated case. This fact makes it necessary to redefine single-particle states. Due to the strong correlation, the initial algebra of operators is transformed into the new algebra of complicated operators. In principle, in terms of the new operators, the initial Hamiltonian can be rewritten as a bilinear form, and the generalized Wick theorem can be formulated. It is very important to stress that the transformation to the new algebra of relevant operators reflects some important internal symmetries of the problem, and nowadays, this way of thinking is formulating in the elegant and very powerful technique of the classification of the integrable models and exactly soluble models (cf.).
* The single-electron GF that describes dynamic properties, should have the two-pole functional structure, which gives in the atomic limit, when the hopping integral tends to zero, the exact two-level atomic solution.
* The GMFs have, in the general case, a very non-trivial structure. The GMFs functional, as a rule, cannot be expressed in terms of the functional of the mean particle densities.

In this section, we consider the case of a large but finite Coulomb repulsion $`U`$ in the Hubbard Hamiltonian (47) . Let us consider the single-particle GF (128) in the Wannier basis

(149)
$$G_{ij\sigma }(tt^{})=<<a_{i\sigma }(t);a_{j\sigma }^{}(t^{})>>$$

It is convenient to introduce the new set of relevant operators

$`d_{i\alpha \sigma }=n_{i\sigma }^\alpha a_{i\sigma },(\alpha =\pm );n_{i\sigma }^+=n_{i\sigma },n_{i\sigma }^{}=(1n_{i\sigma });`$
(150) $`{\displaystyle n_{i\sigma }^\alpha }=1;n_{i\sigma }^\alpha n_{i\sigma }^\beta =\delta _{\alpha \beta }n_{i\sigma }^\alpha ;{\displaystyle \underset{\alpha }{}}d_{i\alpha \sigma }=a_{i\sigma }`$
The new operators $`d_{i\alpha \sigma }`$ and $`d_{j\beta \sigma }^{}`$ have complicated commutation rules, namely,

$$[d_{i\alpha \sigma },d_{j\beta \sigma }^{}]_+=\delta _{ij}\delta _{\alpha \beta }n_{i\sigma }^\alpha $$

The convenience of the new operators follows immediately if one writes down the equation of motion for them

$`[d_{i\alpha \sigma },H]_{}=E_\alpha d_{i\alpha \sigma }+{\displaystyle \underset{ij}{}}t_{ij}(n_{i\sigma }^\alpha a_{j\sigma }+\alpha a_{i\sigma }b_{ij\sigma })`$
(151) $`b_{ij\sigma }=(a_{i\sigma }^{}a_{j\sigma }a_{j\sigma }^{}a_{i\sigma }).`$
It is possible to interpret ,  both contributions to this equation as alloy analogy and resonance broadening corrections. Using the new operator algebra, it is possible identically rewrite GF (149) in the following way

(152)
$$G_{ij\sigma }(\omega )=\underset{\alpha \beta }{}<<d_{i\alpha \sigma }|d_{j\beta \sigma }^{}>>_\omega =\underset{\alpha \beta }{}F_{ij\sigma }^{\alpha \beta }(\omega )$$

The equation of motion for the auxiliary matrix GF

(153)
$$F_{ij\sigma }^{\alpha \beta }(\omega )=\left(\begin{array}{cc}<<d_{i+\sigma }|d_{j+\sigma }^{}>>_\omega & <<d_{i+\sigma }|d_{j\sigma }^{}>>_\omega \\ <<d_{i\sigma }|d_{j+\sigma }^{}>>_\omega & <<d_{i\sigma }|d_{j\sigma }^{}>>_\omega \end{array}\right)$$

is of the following form

(154)
$$(\mathrm{𝐄𝐅}_{ij\sigma }(\omega )𝐈\delta _{ij})_{\alpha \beta }=\underset{li}{}t_{il}<<n_{i\sigma }^\alpha a_{l\sigma }+\alpha a_{i\sigma }b_{il\sigma }|d_{j\beta \sigma }^{}>>_\omega $$

where the following matrix notations was used:

(155)
$$𝐄=\left(\begin{array}{cc}(\omega E_+)& 0\\ 0& (\omega E_{})\end{array}\right);𝐈=\left(\begin{array}{cc}n_\sigma ^+& 0\\ 0& n_\sigma ^{}\end{array}\right).$$

In accordance with the general method of Section 3, we introduce by definition the matrix IGF:

$`𝐃_{il,j}^{(ir)}(\omega )=\left(\begin{array}{cc}<<Z_{11}|d_{j+\sigma }^{}>>_\omega & <<Z_{12}|d_{j\sigma }^{}>>_\omega \\ <<Z_{21}|d_{j+\sigma }^{}>>_\omega & <<Z_{22}|d_{j\sigma }^{}>>_\omega \end{array}\right)`$
(156) $`{\displaystyle \underset{\alpha ^{}}{}}(\left[{\displaystyle \genfrac{}{}{0pt}{}{A_{il}^{+\alpha ^{}}}{A_{il}^\alpha ^{}}}\right][F_{ij\sigma }^{\alpha ^{}+}F_{ij\sigma }^\alpha ^{}]\left[{\displaystyle \genfrac{}{}{0pt}{}{B_{li}^{+\alpha ^{}}}{B_{li}^\alpha ^{}}}\right][F_{lj\sigma }^{\alpha ^{}+}F_{lj\sigma }^\alpha ^{}])`$
Here the notation was used:

$$Z_{11}=Z_{12}=n_{i\sigma }^+a_{l\sigma }+a_{i\sigma }b_{il\sigma };Z_{21}=Z_{22}=n_{i\sigma }^{}a_{l\sigma }a_{i\sigma }b_{il\sigma }$$

It is to be emphasized that the definition (154) is the most important and crucial point of the whole our approach to description of the strong correlation. The coefficients A and B are determined by the orthogonality constraint (31), namely,

(157)
$$<[(𝐃_{il,j}^{(ir)})_{\alpha \beta },d_{j\beta \sigma }^{}]_+>=0$$

After some algebra, we obtain from (157) ($`ij`$)

$`[A_{il}]_{\alpha \beta }=\alpha (<d_{i\beta \sigma }^{}a_{l\sigma }>+<d_{i\beta \sigma }a_{l\sigma }^{}>)(n_\sigma ^\beta )^1`$
$`[B_{li}]_{\alpha \beta }=[<n_{l\sigma }^\beta n_{i\sigma }^\alpha >+\alpha \beta (<a_{i\sigma }a_{i\sigma }^{}a_{l\sigma }a_{l\sigma }^{}>`$
(158) $`<a_{i\sigma }a_{i\sigma }a_{l\sigma }^{}a_{l\sigma }^{}>)](n_\sigma ^\beta )^1`$
As previously, we introduce now GMF GF $`𝐅_{ij\sigma }^0`$; however, as it is clear from (8.2), the actual definition of the GMF GF is very nontrivial. After the Fourier transformation, we get

(159)
$$\left(\begin{array}{cc}F_{k\sigma }^{0++}& F_{k\sigma }^{0+}\\ F_{k\sigma }^{0+}& F_{k\sigma }^0\end{array}\right)=\frac{1}{abcd}\left(\begin{array}{cc}n_\sigma ^+b& n_\sigma ^{}d\\ n_\sigma ^+c& n_\sigma ^{}a\end{array}\right)$$

The coefficients $`a`$, $`b`$, $`c`$, $`d`$ are equal to

$`{\displaystyle \genfrac{}{}{0pt}{}{a}{b}}=\left(\omega E_\pm N^1{\displaystyle \underset{p}{}}ϵ(p)[A^{\pm \pm }(p)B^{\pm \pm }(pq)]\right)`$
(160) $`{\displaystyle \genfrac{}{}{0pt}{}{c}{d}}=N^1{\displaystyle \underset{p}{}}ϵ(p)[A^\pm (p)B^\pm (pq)]`$
Then, using the definition (153), we find the final expression for GMF GF

(161)
$$G_{k\sigma }^{MF}(\omega )=\frac{\omega (n_\sigma ^+E_{}+n_\sigma ^{}E_+)\lambda (k)}{(\omega E_+n_\sigma ^{}\lambda _1(k))(\omega E_{}n_\sigma ^+\lambda _2(k))n_\sigma ^{}n_\sigma ^+\lambda _3(k)\lambda _4(k)}$$

Here we introduced the following notation:

(162) $`{\displaystyle \genfrac{}{}{0pt}{}{\lambda _1(k)}{\lambda _2(k)}}={\displaystyle \frac{1}{n_\sigma ^{}}}{\displaystyle \underset{p}{}}ϵ(p)[A^{\pm \pm }(p)B^{\pm \pm }(pk)]`$
(163) $`{\displaystyle \genfrac{}{}{0pt}{}{\lambda _3(k)}{\lambda _4(k)}}={\displaystyle \frac{1}{n_\sigma ^{}}}{\displaystyle \underset{p}{}}ϵ(p)[A^\pm (p)B^\pm (pk)]`$
$`\lambda (k)=(n_\sigma ^{})^2(\lambda _1+\lambda _3)+(n_\sigma ^+)^2(\lambda _2+\lambda _4)`$
From the equation (161) it is obvious that our two-pole solution is more general than the“Hubbard III”  solution and the Roth solution. Our solution has the correct nonlocal structure and, thus, takes into account the non-diagonal scattering matrix elements more accurately. Those matrix elements describe the virtual “recombination” processes and reflect the extremely complicated structure of single-particle states which virtually include a great number of intermediate scattering processes.
The spectrum of mean-field quasi-particle excitations follows from the poles of the GF (161) and consists of two branches

(164) $`\omega {\displaystyle \genfrac{}{}{0pt}{}{1}{2}}(k)=`$
$`1/2[(E_+E_{}+a_1+b_1)\pm \sqrt{(E_++E_{}a_1b_1)^24cd}]`$
where $`a_1=\omega E_\pm a;b_1=\omega E_\pm b`$. Thus, the spectral weight function $`A_{k\sigma }(\omega )`$ of GF (161) consists of two peaks separated by the distance

(165)
$$\omega _1\omega _2=\sqrt{(Ua_1b_1)^2cd}U(1\frac{a_1b_1}{U})+O(\gamma )$$

For a deeper insight into the functional structure of the solution (161) and to compare with other solutions, we rewrite (161) in the following form

(166)
$$𝐅_{k\sigma }^0(\omega )=\left(\begin{array}{cc}(\frac{a}{n_\sigma ^+}\frac{db^1c}{n_\sigma ^+})^1& \frac{d}{a}(\frac{b}{n_\sigma ^{}}\frac{da^1c}{n_\sigma ^{}})^1\\ \frac{c}{b}(\frac{a}{n_\sigma ^+}\frac{db^1c}{n_\sigma ^+})^1& (\frac{b}{n_\sigma ^{}}\frac{db^1c}{n_\sigma ^{}})^1\end{array}\right)$$

from which we obtain for $`G_\sigma ^{MF}(k,\omega )`$

$`G_{k\sigma }^{MF}(\omega )={\displaystyle \frac{n_\sigma ^+(1+cb^1)}{adb^1c}}+{\displaystyle \frac{n_\sigma ^{}(1+da^1)}{bca^1d}}`$
(167) $`{\displaystyle \frac{n_\sigma ^{}}{\omega E_{}n_\sigma ^+W_{k\sigma }^{}}}+{\displaystyle \frac{n_\sigma ^+}{\omega E_+n_\sigma ^{}W_{k\sigma }^{}}}`$
where

(168) $`n_\sigma ^+n_\sigma ^{}W_{k\sigma }^\pm =N^1{\displaystyle \underset{ij}{}}t_{ij}\mathrm{exp}[ik(R_iR_j)]`$
$`((<a_{i\sigma }^{}n_{i\sigma }^\pm a_{j\sigma }>+<a_{i\sigma }n_{i\sigma }^{}a_{j\sigma }^{}>)+`$
$`(<n_{j\sigma }^\pm n_{i\sigma }^\pm >+<a_{i\sigma }a_{i\sigma }^{}a_{j\sigma }a_{j\sigma }^{}><a_{i\sigma }a_{i\sigma }a_{j\sigma }^{}a_{j\sigma }^{}>))`$
are the shifts for upper and lower splitted subbands due to the elastic scattering of carriers in the Generalized Mean Field. The quantities $`W^\pm `$ are functionals of the GMF. The most important feature of the present solution of the strongly correlated Hubbard model is a very nontrivial structure of the mean-field renormalizations (8.2), which is crucial for understanding the physics of strongly correlated systems. It is important to emphasize that just this complicated form of GMF is only relevant to the essence of the physics under consideration. The attempts to reduce the functional of GMF to a simpler functional of the average density of electrons are incorrect from the point of view of real physics of strongly correlated systems. This physics clearly shows that the mean-field renormalizations cannot be expressed as functionals of the electron mean density. To explain this statement, let us derive the “Hubbard I” solution  (52) from our GMF solution (161). If we approximate (8.2) as

(169)
$$n_\sigma ^+n_\sigma ^{}W^\pm (k)N^1\underset{ij}{}t_{ij}\mathrm{exp}[ik(R_iR_j)]<n_{j\sigma }^\pm n_{i\sigma }^\pm >$$

and make the additional approximation, namely,

$$<n_{j\sigma }n_{i\sigma }>n_\sigma ^2$$

then solution (161) turns into the “Hubbard I” solution (52). This solution, as it is well known, is unrealistic from many points of view.
As to our solution (161) , the second important aspect is that the parameters $`\lambda _i(k)`$ do not depend on frequency, since they depend essentially on elastic scattering processes. The dependence on frequency arises due to inelastic scattering processes which are contained in our self-energy operator. We proceed now with the derivation of the explicit expression for the self-energy.
To calculate a high-order GF on the r.h.s. of (154), we should use the second time variable ($`t^{}`$) differentiation of it again. If one introduces the irreducible parts for the right-hand-side operators by analogy with the expression (8.2), the equation of motion (154) can be rewritten exactly in the following form

(170)
$$𝐅_{k\sigma }(\omega )=𝐅_{k\sigma }^0(\omega )+𝐅_{k\sigma }^0(\omega )𝐏_{k\sigma }(\omega )𝐅_{k\sigma }^0(\omega )$$

Here the scattering operator $`P`$ (36) is of the form

(171)
$$𝐏_{q\sigma }(\omega )=𝐈^1[\underset{lm}{}t_{il}t_{mj}<<𝐃_{il,j}^{(ir)}|𝐃_{i,mj}^{(ir)}>>_\omega ]_q𝐈^1$$

In accordance with the definition (37), we write down the Dyson equation

(172)
$$𝐅=𝐅^0+𝐅^0\mathrm{𝐌𝐅}$$

The self-energy operator $`M`$ is defined by eq. (37). Let us note again that the self-energy corrections, according to (38), contribute to the full GF as additional terms. This is an essential advantage in comparison with the “Hubbard III” solution and other two-pole solutions. It is clear from the form of Roth solution (53) that it includes the elastic scattering corrections only and does not incorporate the damping effects and finite lifetimes.
For the full GF we find, using the formal solution of Dyson equation (36), that it is equal to

$`G_{k\sigma }(\omega )=({\displaystyle \frac{1}{n_\sigma ^+}}(an_\sigma ^+M_{k\sigma }^{++}(\omega ))+{\displaystyle \frac{1}{n_\sigma ^{}}}(bn_\sigma ^{}M_{k\sigma }^{}(\omega ))`$
$`+{\displaystyle \frac{1}{n_\sigma ^+}}(d+n_\sigma ^+M_{k\sigma }^+(\omega ))+{\displaystyle \frac{1}{n_\sigma ^{}}}(c+n_\sigma ^{}M_{k\sigma }^+(\omega )))`$
(173) $`[\mathrm{det}\left((F_{k\sigma }^0(\omega ))^1M_{k\sigma }(\omega )\right)]^1`$
After some algebra, we can rewrite this expression in the following form which is essentially new and, in a certain sense, is the central result of the present theory:

(174)
$$G=\frac{\omega (n^+E_{}+n^{}E_+)L}{(\omega E_+n^{}L_1)(\omega E_{}n^+L_2)n^{}n^+L_3L_4}$$

where

$`L_1(k,\omega )=\lambda _1(k){\displaystyle \frac{n_\sigma ^+}{n_\sigma ^{}}}M_\sigma ^{++}(k,\omega );`$
$`L_2(k,\omega )=\lambda _2(k){\displaystyle \frac{n_\sigma ^{}}{n_\sigma ^+}}M_\sigma ^{}(k,\omega );`$
$`L_3(k,\omega )=\lambda _3(k)+{\displaystyle \frac{n_\sigma ^{}}{n_\sigma ^+}}M_\sigma ^+(k,\omega );`$
$`L_4(k,\omega )=\lambda _4(k)+{\displaystyle \frac{n_\sigma ^+}{n_\sigma ^{}}}M_\sigma ^+(k,\omega );`$
(175) $`L(k,\omega )=\lambda (k)+n_\sigma ^+n_\sigma ^{}(M^{++}+M^{}M^+M^+)`$
Thus, now we have to find explicit expressions for the elements of the self-energy matrix M. To this end, we should use the spectral theorem again to express the GF in terms of correlation functions

(176)
$$M_{k\sigma }^{\alpha ,\beta }(\omega )<D_{mj,\beta }^{(ir)}(t)D_{il,\alpha }^{(ir)}>$$

For the approximate calculation of the self-energy, we propose to use the following trial solution

$`<D^{(ir)}(t)D^{(ir)}><a_{m\sigma }^{}(t)a_{l\sigma }><n_{j\sigma }^\beta (t)n_{i\sigma }^\alpha >`$
$`+<a_{m\sigma }^{}(t)n_{i\sigma }^\alpha ><n_{j\sigma }^\beta (t)a_{l\sigma }>+\beta <b_{mj\sigma }^{}(t)a_{l\sigma }><a_{j\sigma }^{}(t)n_{i\sigma }^\alpha >`$
$`+\beta <b_{mj\sigma }^{}(t)n_{i\sigma }^\alpha ><a_{j\sigma }^{}(t)a_{l\sigma }>+\alpha <a_{m\sigma }^{}(t)a_{i\sigma }><n_{j\sigma }^\beta (t)b_{il\sigma }>`$
$`+\alpha <a_{m\sigma }^{}(t)b_{il\sigma }><n_{j\sigma }^\beta (t)b_{il\sigma }>`$
$`+\alpha \beta <b_{mj\sigma }^{}(t)a_{i\sigma }><a_{j\sigma }^{}(t)b_{il\sigma }>`$
(177) $`+\alpha \beta <b_{mj\sigma }^{}(t)b_{il\sigma }><a_{j\sigma }^{}(t)a_{i\sigma }>`$
It is quite natural to interpret the contributions into this expression in terms of scattering, resonance-broadening, and interference corrections of different types. For example, let us consider the simplest approximation. For this aim, we retain the first contribution in (8.2)

$`[\mathrm{𝐈𝐌𝐈}]_{\alpha \beta }={\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\omega ^{}}{\omega \omega ^{}}}(\mathrm{exp}(\beta \omega ^{})+1)`$
$`{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{dt}{2\pi }}\mathrm{exp}(i\omega ^{}t)N^1{\displaystyle \underset{ijlm}{}}exp[ik(R_iR_j)]t_{il}t_{mj}`$
(178) $`{\displaystyle 𝑑\omega _1n(\omega _1)\mathrm{exp}(i\omega _1t)g_{ml\sigma }(\omega _1)\left(\frac{1}{\pi }ImK_{ij}^{\alpha \beta }(\omega _1\omega ^{})\right)}.`$
Here $`K_{ij}^{\alpha \beta }(\omega )=<<n_{i\sigma }^\alpha |n_{j\sigma }^\beta >>_\omega `$ is the density-density GF. It is worthy to note that the mass operator (8.2) contains the term $`t_{il}t_{mj}`$ contrary to the expression (143) that contains the term $`U^2`$. The pair of equations (8.2) and (172) is a self-consistent system of equations for the single-particle Green function. For a simple estimation, for the calculation of the self-energy (8.2), it is possible to use any initial relevant approximation of the two-pole structure. As an example, we take the expression (52). We then obtain

$`[\mathrm{𝐈𝐌𝐈}]_{\alpha \beta }{\displaystyle \underset{q}{}}|ϵ(kq)|^2K_q^{\alpha \beta }`$
(179) $`[{\displaystyle \frac{n_\sigma }{\omega Uϵ(kq)n_\sigma }}+{\displaystyle \frac{1n_\sigma }{\omega ϵ(kq)(1n_\sigma )}}]`$
In the same way, one can use, instead (52), another initial two-pole solution, e.g. the Roth solution (53), etc.
On the basis of the self-energy operator (8.2) we can explicitly find the energy shift and damping due to inelastic scattering of quasi-particles. This is a great advantage of the present approach.
In summary, in this Section, we obtained the most complete solution to the Hubbard model Hamiltonian in the strongly correlated case. It has correct functional structure, and, moreover, it represents correctly the effects of elastic and inelastic scattering in a systematic and convenient way. The mass operator contains all inelastic scattering terms including various scattering and resonance broadening terms in a systematic way. The obtained solution (174) is valid for all band filling and for arbitrarily strong but finite strength of the Coulomb repulsion. Our solution contains no approximations except those contained in the final calculation of the mass operator. Therefore, we conclude that our solutions to the Hubbard model in the weakly correlated case (139) and in the strongly correlated case (174) describe most fully and self-consistently the correlation effects in the Hubbard model and give a unified interpolation description of the correlation problem. This result is to be contrasted with Hubbard , Roth and many other results in which this interpolation solution cannot be derived within the unified scheme.
It is clear from the present consideration that for the systematic construction of the advanced approximate solutions we need to calculate the collective correlation functions of the electron density and spin density and the density of doubles, but this problem must be considered separately.

### 8.3 Correlations in Random Hubbard Model

In this Section, we apply the IGF method to consider the electron-electron correlations in the presence of disorder to demonstrate the advantage of our approach. The treatment of the electron motion in substitutionally disordered $`A_xB_{1x}`$ transition metal alloys is based upon a certain generalization of the Hubbard model, including random diagonal and off-diagonal elements caused by substitutional disorder in a binary alloy. The electron-electron interaction plays an important role for various aspects of behaviour in alloys, e.g. for the weak localization . The approximation which is used widely for treating disordered alloys is the single-site Coherent Potential Approximation (CPA) . The CPA has been refined and developed in many papers (e.g. , ) and till now is the most popular approximation for the theoretical study of alloys. But the simultaneous effect of disorder and electron-electron inelastic scattering has been considered for some limited cases only and not within the self-consistent scheme.
Let us consider the Hubbard model Hamiltonian (47) on a given configuration of an alloy $`(\nu )`$

(180)
$$H^{(\nu )}=H_1^{(\nu )}+H_2^{(\nu )}$$

where

$`H_1^{(\nu )}={\displaystyle \underset{i\sigma }{}}\epsilon _i^\nu n_{i\sigma }+{\displaystyle \underset{ij\sigma }{}}t_{ij}^{\nu \mu }a_{i\sigma }^{}a_{j\sigma }`$
(181) $`H_2^{(\nu )}={\displaystyle \frac{1}{2}}{\displaystyle \underset{i\sigma }{}}U_i^\nu n_{i\sigma }n_{i\sigma }`$
Contrary to the periodic model (47), the atomic level energy $`\epsilon _i^\nu `$, the hopping integrals $`t_{ij}^{\nu \mu }`$, as well as the intraatomic Coulomb repulsion $`U_i^\nu `$ are here random variables which take the values $`\epsilon ^\nu `$, $`t^{\nu \mu }`$, and $`U^\nu `$, respectively; the superscript $`\nu (\mu )`$ refers to atomic species ($`\nu ,\mu =A,B`$) located on site i(j). The nearest-neighbor hopping integrals were only included .
To unify the IGF method and CPA into a completely self-consistent scheme let us consider the single-electron GF (149) $`G_{ij\sigma }`$ in the Wannier representation for a given configuration $`(\nu )`$. The corresponding equation of motion is of the form (for brevity we omit the superscript $`(\nu )`$ where its presence is clear)

$`(\omega \epsilon _i)<<a_{i\sigma }|a_{j\sigma }^{}>>_\omega =\delta _{ij}+{\displaystyle \underset{n}{}}t_{in}<<a_{n\sigma }|a_{j\sigma }^{}>>_\omega `$
(182) $`+U_i<<n_{i\sigma }a_{i\sigma }|a_{j\sigma }^{}>>_\omega `$
In the present paper, for brevity, we confine ourselves to the weak correlation and the diagonal disorder case. The generalization to the case of strong correlation or off-diagonal disorder is straightforward, but its lengthy consideration preclude us from discussing it this time.
Using the definition (30), we define the IGF for a given (fixed) configuration of atoms in an alloy as follows

(183) $`{}_{}{}^{(ir)}<<n_{i\sigma }a_{i\sigma }|a_{j\sigma }^{}>>=`$
$`<<n_{i\sigma }a_{i\sigma }|a_{j\sigma }^{}>><n_{i\sigma }><<a_{i\sigma }|a_{j\sigma }^{}>>`$
This time, contrary to (159), because of lack of translational invariance we must take into account the site dependence of $`<n_{i\sigma }>`$. Then we rewrite the equation of motion (8.3) in the following form

$`{\displaystyle \underset{n}{}}[(\omega \epsilon _iU_i<n_{i\sigma }>)\delta _{ij}t_{in}]<<a_{n\sigma }|a_{j\sigma }^{}>>_\omega =`$
(184) $`\delta _{ij}+U_i(^{(ir)}<<n_{i\sigma }a_{i\sigma }|a_{j\sigma }^{}>>_\omega )`$
In accordance with the general method of Section 3, we find then the Dyson equation for a given configuration $`(\nu )`$

(185)
$$G_{ij\sigma }(\omega )=G_{ij\sigma }^0(\omega )+\underset{mn}{}G_{im\sigma }^0(\omega )M_{mn\sigma }(\omega )G_{nj\sigma }(\omega )$$

The GMF GF $`G_{ij\sigma }^0`$ and the self-energy operator $`M`$ are defined as

$`{\displaystyle \underset{m}{}}H_{im\sigma }G_{mj\sigma }^0(\omega )=\delta _{ij}`$
$`P_{mn\sigma }=M_{mn\sigma }+{\displaystyle \underset{ij}{}}M_{mi\sigma }G_{ij\sigma }^0P_{jn\sigma }`$
$`H_{im\sigma }=(\omega \epsilon _iU_i<n_{i\sigma }>)\delta _{im}t_{im}`$
(186) $`P_{mn\sigma }(\omega )=U_m(^{(ir)}<<n_{m\sigma }a_{m\sigma }|n_{n\sigma }a_{n\sigma }^{}>>_\omega ^{(ir)})U_n`$
In order to calculate the self-energy operator $`M`$ self-consistently, we have to express it approximately by the lower-order GFs. Employing the same pair approximation as (8.1) (now in the Wannier representation) and the same procedure of calculation, we arrive at the following expression for $`M`$ for a given configuration $`(\nu )`$

$`M_{mn\sigma }^{(\nu )}(\omega )=U_mU_n{\displaystyle \frac{1}{2\pi ^4}}{\displaystyle R(\omega _1,\omega _2,\omega _3)}`$
$`ImG_{nm\sigma }^{(\nu )}(\omega _1)ImG_{mn\sigma }^{(\nu )}(\omega _2)ImG_{mn\sigma }^{(\nu )}(\omega _3);`$
(187) $`R={\displaystyle \frac{d\omega _1d\omega _2d\omega _3}{\omega +\omega _1\omega _2\omega _3}}{\displaystyle \frac{(1n(\omega _1))n(\omega _2)n(\omega _3)}{n(\omega _2+\omega _3\omega _1)}}`$
As we mentioned previously, all the calculations just presented were made for a given configuration of atoms in an alloy. All the quantities in our theory (G, $`G^0`$, P, M) depend on the whole configuration of the alloy. To obtain a theory of a real macroscopic sample, we have to average over various configurations of atoms in the sample. The configurational averaging cannot be exactly made for a macroscopic sample. Hence we must resort to an additional approximation. It is obvious that the self-energy $`M`$ is in turn a functional of $`G`$, namely $`M=M[G]`$. If the process of making configurational averaging is denoted by $`\overline{G}`$, then we have

$$\overline{G}=\overline{G}^0+\overline{G^0MG}$$

A few words are now appropriate for the description of general possibilities. The calculations of $`\overline{G}^0`$ can be performed with the help of any relevant available scheme. In the present work, for the sake of simplicity, we choose the single-site CPA, namely, we take

(188)
$$\overline{G}_{mn\sigma }^0(\omega )=N^1\underset{k}{}\frac{\mathrm{exp}(ik(R_mR_n))}{\omega \mathrm{\Sigma }^\sigma (\omega )ϵ(k)}$$

Here $`ϵ(k)=_{n=1}^zt_{n,0}\mathrm{exp}(ikR_n)`$, $`z`$ is the number of nearest neighbors of the site $`0`$, and the coherent potential $`\mathrm{\Sigma }^\sigma (\omega )`$ is the solution of the CPA self-consistency equations. For the $`A_xB_{1x}`$, we have

$`\mathrm{\Sigma }^\sigma (\omega )=x\epsilon _A^\sigma +(1x)\epsilon _B^\sigma (\epsilon _A^\sigma \mathrm{\Sigma }^\sigma )F^\sigma (\omega ,\mathrm{\Sigma }^\sigma )(\epsilon _B^\sigma \mathrm{\Sigma }^\sigma );`$
(189) $`F^\sigma (\omega ,\mathrm{\Sigma }^\sigma )=\overline{G}_{mm\sigma }^0(\omega )`$
Now, let us return to the calculation of the configurationally averaged total GF $`\overline{G}`$. To perform the remaining averaging in the Dyson equation, we use the approximation

$$\overline{G^0MG}\overline{G}^0\overline{M}\overline{G}$$

The calculation of $`\overline{M}`$ requires further averaging of the product of matrices. We again use the prescription of the factorizability there, namely

$$\overline{M}\overline{(U_mU_n)}\overline{(ImG)}\overline{(ImG)}\overline{(ImG)}$$

However, the quantities $`\overline{U_mU_n}`$ entering into $`\overline{M}`$ are averaged here according to

$`\overline{U_mU_n}=U_2+(U_1U_2)\delta _{mn}`$
$`U_1=x^2U_A^2+2x(1x)U_AU_B+(1x)^2U_B^2`$
(190) $`U_2=xU_A^2+(1x)U_B^2`$
The averaged value for the self-energy is

(191) $`\overline{M}_{mn\sigma }(\omega )`$
$`={\displaystyle \frac{U_2}{2\pi ^4}}{\displaystyle R(\omega _1,\omega _2,\omega _3)Im\overline{G}_{nm\sigma }(\omega _1)Im\overline{G}_{mn\sigma }(\omega _2)Im\overline{G}_{mn\sigma }(\omega _3)}+`$
$`{\displaystyle \frac{U_1U_2}{2\pi ^4}}\delta _{mn}{\displaystyle R(\omega _1,\omega _2,\omega _3)Im\overline{G}_{nm\sigma }(\omega _1)Im\overline{G}_{mn\sigma }(\omega _2)Im\overline{G}_{mn\sigma }(\omega _3)}`$
The averaged quantities are periodic, so we can introduce the Fourier transform of them, i.e.

$$\overline{M}_{mn\sigma }(\omega )=N^1\underset{k}{}\overline{M}_{k\sigma }(\omega )\mathrm{exp}(ik(R_mR_n))$$

and similar formulae for $`\overline{G}`$ and $`\overline{G}^0`$. Performing the configurational averaging of the Dyson equation and Fourier transforming of the resulting expressions according to the above rules, we obtain

(192)
$$\overline{G}_{k\sigma }(\omega )=[\omega ϵ(k)\mathrm{\Sigma }^\sigma (\omega )\overline{M}_{k\sigma }(\omega )]^1$$

where

$`\overline{M}_{k\sigma }(\omega )={\displaystyle \frac{1}{2\pi ^4}}{\displaystyle \underset{pq}{}}{\displaystyle R(\omega _1,\omega _2,\omega _3)N^2Im\overline{G}_{pq\sigma }(\omega _1)Im\overline{G}_{q\sigma }(\omega _2)}`$
(193) $`[U_2Im\overline{G}_{k+p\sigma }(\omega _3)+{\displaystyle \frac{(U_1U_2)}{N}}{\displaystyle \underset{g}{}}Im\overline{G}_{k+pg}(\omega _3)]`$
The simplest way to obtain an explicit solution for the self-energy$`\overline{M}`$ is to start with a suitable initial trial solution as it was done for the periodic case. For a disordered system, it is reasonable to use, as the first iteration approximation the so-called Virtual Crystal Approximation(VCA):

$$\frac{1}{\pi }Im\overline{G}_{k\sigma }^{VCA}(\omega +iϵ)\delta (\omega E_k^\sigma )$$

where for the binary alloy $`A_xB_{1x}`$ this approximation reads

$$\overline{V}=xV^A+(1x)V^B;E_k^\sigma =\overline{\epsilon }_i^\sigma +ϵ(k);$$

$$\overline{\epsilon }_i^\sigma =x\epsilon _A^\sigma +(1x)\epsilon _B^\sigma $$

Note, that the use of VCA here is by no means a solution of the correlation problem in VCA. It is only the use of the VCA for the parametrization of the problem, to start with VCA input parameters. After the integration of (8.3) the final result for the self-energy is

(194) $`\overline{M}_{k\sigma }(\omega )=`$
$`{\displaystyle \frac{U_2}{N^2}}{\displaystyle \underset{pq}{}}{\displaystyle \frac{n(E_{p+q}^\sigma )[1n(E_q^\sigma )n(E_{k+p}^\sigma )]+n(E_{k+p}^\sigma )n(E_q^\sigma )}{\omega +E_{p+q}^\sigma E_q^\sigma E_{k+p}^\sigma }}+`$
$`{\displaystyle \frac{(U_1U_2)}{N^3}}{\displaystyle \underset{pqg}{}}{\displaystyle \frac{n(E_{p+q}^\sigma )[1n(E_q^\sigma )n(E_{k+pg}^\sigma )]+n(E_{k+pg}^\sigma )n(E_q^\sigma )}{\omega +E_{p+q}^\sigma E_q^\sigma E_{k+pg}^\sigma }}`$
It is to be emphasized that the equations (191) - (194) give the general microscopic self-consistent description of inelastic electron-electron scattering in an alloy in the spirit of the CPA. We took into account the randomness not only through the parameters of the Hamiltonian but also in a self-consistent way through the configurational dependence of the self-energy operator.

### 8.4 Electron-Lattice Interaction and MTBA

To understand quantitatively the electrical, thermal, and superconducting properties of metals and their alloys, one needs a proper description of an electron-lattice interaction too ,. A systematic, self-consistent simultaneous treatment of the electron-electron and electron-phonon interaction plays an important role in recent studies of strongly correlated systems. It was argued from different points of view that to understand quantitatively the phenomenon of high-temperature superconductivity one needs a proper inclusion of electron-phonon interaction, too. A lot of theoretical searches for the relevant mechanism of high temperature superconductivity deal with strong electron-phonon interaction models. The natural approach to the description of superconductivity in that type of compounds is the modified tight-binding approximation (MTBA) , . The papers ,, , contain a self-consistent microscopic theory of the normal and superconducting properties of transition metals and strongly disordered binary alloys in the framework of the Hubbard Model (47) and random Hubbard model (180). Here we derive a system of equations for the superconductivity for tight-binding electrons of a transition metal interacting with phonons within the IGF approach. We write the total Hamiltonian of the electron-ion system as the sum

(195)
$$H=H_e+H_i+H_{ei}$$

where $`H_e`$ is the electron part of the Hamiltonian represented by the Hubbard operator (47). The Hamiltonian of an ion subsystem and the operator of electron-ion interaction have the form

(196) $`H_i={\displaystyle \frac{1}{2}}{\displaystyle \underset{n}{}}{\displaystyle \frac{P_n^2}{2M}}+{\displaystyle \frac{1}{2}}{\displaystyle \underset{mn\alpha \beta }{}}\mathrm{\Phi }_{nm}^{\alpha \beta }u_n^\alpha u_m^\beta `$
(197) $`H_{ei}={\displaystyle \underset{\sigma }{}}{\displaystyle \underset{n,ij}{}}V_{ij}^\alpha (\stackrel{}{R}_n^0)a_{i\sigma }^{}a_{j\sigma }u_n^\alpha `$
where

(198)
$$\underset{n}{}V_{ij}^\alpha (\stackrel{}{R}_n^0)u_n^\alpha =\frac{t_{ij}(\stackrel{}{R}_{ij}^0)}{R_{ij}^0}(\stackrel{}{u}_i\stackrel{}{u}_j)$$

Here $`P_n`$ is the momentum operator, $`M`$ is the mass of an ion, and $`u_n`$ is the displacement of the ion from the equilibrium position at the lattice site $`R_n`$.
In a more convenient notation the electron-phonon interaction Hamiltonian in the modified tight-binding approximation reads

(199)
$$H_{ei}=\underset{\nu \sigma }{}\underset{kq}{}V^\nu (\stackrel{}{k},\stackrel{}{k}+\stackrel{}{q})Q_{\stackrel{}{q}\nu }a_{k+q\sigma }^{}a_{k\sigma }$$

where

(200)
$$V^\nu (\stackrel{}{k},\stackrel{}{k}+\stackrel{}{q})=\frac{2iq_0}{(NM)^{1/2}}\underset{\alpha }{}t(\stackrel{}{a}_\alpha )e_\nu ^\alpha (\stackrel{}{q})[\mathrm{sin}\stackrel{}{a}_\alpha \stackrel{}{k}\mathrm{sin}\stackrel{}{a}_\alpha (\stackrel{}{k}\stackrel{}{q})]$$

here $`q_0`$ is the Slater coefficient  having the origin in the exponential decrease of the wave functions of $`d`$-electrons, N is the number of unit cells in the crystal, and M is the ion mass. The quantities $`\stackrel{}{e}_\nu (\stackrel{}{q})`$ are polarization vectors of the phonon modes.
For the ion subsystem, we have

(201)
$$H_i=\frac{1}{2}\underset{q\nu }{}(P_{q\nu }^{}P_{q\nu }+\omega ^2(\stackrel{}{q}\nu )Q_{q\nu }^{}Q_{q\nu })$$

where $`P_{q\nu }`$ and $`Q_{q\nu }`$ are normal coordinates, and $`\omega (q\nu )`$ are acoustical phonon frequencies. It is important to note that in spite of the fact that in Hubbard model (47), the $`d`$\- and $`s(p)`$-bands are replaced by one effective band , the $`s`$-electrons give rise to screening effects and were taken into effects by choosing the proper values of $`U`$ and the acoustical phonon frequencies.

### 8.5 Equations of Superconductivity

To derive the superconductivity equations, we use the IGF method of Section 3 in which the decoupling procedure is carried out only for approximate calculation of the mass operator of the matrix electron GF. According to the arguments of Section 4.3 , eqn.(62), the relevant matrix GF is of the form

(202) $`G_{ij}(\omega )=\left(\begin{array}{cc}G_{11}& G_{12}\\ G_{21}& G_{22}\end{array}\right)=`$
$`\left(\begin{array}{cc}<<a_{i\sigma }|a_{j\sigma }^{}>>& <<a_{i\sigma }|a_{j\sigma }>>\\ <<a_{i\sigma }^{}|a_{j\sigma }^{}>>& <<a_{i\sigma }^{}|a_{j\sigma }>>\end{array}\right)`$
As was discussed in Section 4.4, with this definition, one introduces the so-called anomalous (off-diagonal) GFs which fix the relevant BCS-Bogoliubov vacuum and select proper symmetry broken solutions. Differentiation of $`G_{ij}(tt^{})`$ with respect to the first time gives for the Fourier components of the equations of motion

(203) $`{\displaystyle \underset{j}{}}(\omega \delta _{ij}t_{ij})<<a_{j\sigma }|a_{i^{}\sigma }^{}>>=\delta _{ii^{}}+`$
$`U<<a_{i\sigma }n_{i\sigma }|a_{i^{}\sigma }^{}>>+{\displaystyle \underset{nj}{}}V_{ijn}<<a_{j\sigma }u_n|a_{i^{}\sigma }^{}>>`$
(204) $`{\displaystyle \underset{j}{}}(\omega \delta _{ij}+t_{ij})<<a_{j\sigma }^{}|a_{i^{}\sigma }^{}>>=`$
$`U<<a_{i\sigma }^{}n_{i\sigma }|a_{i^{}\sigma }^{}>>+{\displaystyle \underset{nj}{}}V_{jin}<<a_{j\sigma }^{}u_n|a_{i^{}\sigma }^{}>>`$
Following the general strategy of the IGF method, we separate the renormalization of the electron energy in the Hartree-Fock-Bogoliubov generalized mean field approximation (including anomalous averages) from the renormalization of higher-order due to inelastic scattering. For this, we introduce irreducible parts of the GF in accordance with the definition ( as an example, we take two of the four Green functions)

(205) $`(^{(ir)}<<a_{i\sigma }a_{i\sigma }^{}a_{i\sigma }|a_{i^{}\sigma }^{}>>_\omega )=<<a_{i\sigma }a_{i\sigma }^{}a_{i\sigma }|a_{i^{}\sigma }^{}>>_\omega `$
$`<n_{i\sigma }>G_{11}+<a_{i\sigma }a_{i\sigma }><<a_{i\sigma }^{}|a_{i^{}\sigma }^{}>>_\omega `$
$`(^{(ir)}<<a_{i\sigma }^{}a_{i\sigma }a_{i\sigma }^{}|a_{i^{}\sigma }^{}>>_\omega )=<<a_{i\sigma }^{}a_{i\sigma }a_{i\sigma }^{}|a_{i^{}\sigma }^{}>>_\omega `$
$`<n_{i\sigma }>G_{21}+<a_{i\sigma }^{}a_{i\sigma }^{}><<a_{i\sigma }|a_{i^{}\sigma }^{}>>_\omega `$
From this definition it follows that this way of introducing the IGF broadens the initial algebra of the operators and the initial set of the GFs. This means that “actual” algebra of the operators must include the anomalous terms from the beginning, namely: $`(a_{i\sigma }`$, $`a_{i\sigma }^{}`$, $`n_{i\sigma }`$, $`a_{i\sigma }^{}a_{i\sigma }^{}`$, $`a_{i\sigma }a_{i\sigma }^{})`$. The corresponding initial GF is the form (202). The choice of the irreducible parts of the GF in (205) is specified by the ”orthogonality” constraint (31), which makes it possible to introduce unambiguously the irreducible parts and make the inhomogeneous terms in the equations for them vanish. Using (205) , we rewrite eqs.(203) and (204) in the form

(206) $`{\displaystyle \underset{j}{}}\left((\omega U<n_{j\sigma }>)\delta _{ij}t_{ij}\right)<<a_{j\sigma }|a_{i^{}\sigma }^{}>>=\delta _{ii^{}}`$
$`U<a_{i\sigma }a_{i\sigma }><<a_{i\sigma }|a_{i^{}\sigma }^{}>>+{\displaystyle \underset{j}{}}<<(a_{j\sigma }\rho _{ij\sigma })^{(ir)}|a_{i^{}\sigma }^{}>>`$
(207) $`{\displaystyle \underset{j}{}}\left((\omega +U<n_{j\sigma }>)\delta _{ij}+t_{ji}\right)<<a_{j\sigma }^{}|a_{i^{}\sigma }^{}>>=`$
$`+U<a_{i\sigma }^{}a_{i\sigma }^{}><<a_{i\sigma }|a_{i^{}\sigma }^{}>>{\displaystyle \underset{j}{}}<<(\rho _{ji\sigma }a_{j\sigma }^{})^{(ir)}|a_{i^{}\sigma }^{}>>`$
where

(208)
$$\rho _{ij\sigma }=Un_{j\sigma }\delta _{ij}+\underset{n}{}V_{ijn}u_n(1\delta _{ij})$$

In the representation of the Nambu operators

(209)
$$\psi _{i,\sigma }=\left(\genfrac{}{}{0pt}{}{a_{i\sigma }}{a_{i\sigma }^{}}\right)\psi _{i,\sigma }^{}=(a_{i\sigma }^{},a_{i\sigma })$$

the equation of motion for GF (206) can be represented as

(210) $`{\displaystyle \underset{j}{}}(\omega \tau _0\delta _{ij}t_{ij}\tau _3\mathrm{\Sigma }_{i\sigma }^c)<<\psi _j|\psi _i^{}^{}>>=`$
$`\delta _{ii^{}}\tau _0+{\displaystyle \underset{j}{}}<<(\rho _{ij}\tau _3\psi _j)^{(ir)}|\psi _i^{}^{}>>`$
Here the Hartree-Fock-Bogoliubov elastic Coulomb term (62) is of the form

(211)
$$\mathrm{\Sigma }_{i\sigma }^c=U\tau _3<\psi _{i,\sigma }\psi _{i,\sigma }^{}>\tau _3+\frac{U}{2}(\tau _0+\tau _3)$$

To calculate the irreducible matrix GF in (210), we write down for it the equation of motion with respect to the second time $`t^{}`$ and then separate the irreducible part with respect to the operators on the right-hand-side of the corresponding GF. This gives the Dyson equation in the matrix form

(212)
$$\widehat{G}_{ii^{}}(\omega )=\widehat{G}_{ii^{}}^0(\omega )+\underset{jj^{}}{}\widehat{G}_{ij}^0(\omega )\widehat{M}_{jj^{}}(\omega )\widehat{G}_{j^{}i^{}}(\omega )$$

The generalized mean field GF $`G^0`$ and the mass operator are defined by

(213)
$$\underset{j}{}(\omega \tau _0\delta _{ij}t_{ij}\tau _3\mathrm{\Sigma }_{i\sigma }^c)G_{ji^{}}^0=\delta _{ii^{}}\tau _0$$

(214)
$$M_{kk^{}}=\underset{jj^{}}{}(<<(\rho _{kj}\tau _3\psi _j)^{(ir)}|(\psi _j^{}^{}\tau _3\rho _{j^{}k^{}})^{(ir)}>>)_\omega ^{(p)}$$

The explicit expression for the mass operator (214) is of the form

(215) $`\widehat{M}_{ii^{}}(\omega )=`$
$`{\displaystyle \underset{jj^{}}{}}\left(\begin{array}{cc}(^{(ir)}<<a_j\rho _{ij}|\rho _{j^{}i^{}}a_j^{}^{}>>^{(ir)})^{(p)}& (^{(ir)}<<a_j\rho _{ij}|\rho _{j^{}i^{}}a_j^{}>>^{(ir)})^{(p)}\\ (^{(ir)}<<a_j^{}\rho _{ji}|\rho _{j^{}i^{}}a_j^{}^{}>>^{(ir)})^{(p)}& (^{(ir)}<<a_j^{}\rho _{ji}|\rho _{i^{}j^{}}a_j^{}>>^{(ir)})^{(p)}\end{array}\right)`$
The mass operator (215) describes inelastic scattering of electrons ( the elastic part is contained in $`\mathrm{\Sigma }_{i\sigma }^c`$ ) on fluctuations of the density of a total electron-ion charge in the lattice. To find an approximating expression for the mass operator (215), we adopt the following trial approximation

(216)
$$<\rho _{j^{}i^{}\sigma }(t)a_{j^{}\sigma }^{}(t)a_{j\sigma }\rho _{ij\sigma }>^{(ir)}<\rho _{j^{}i^{}\sigma }(t)\rho _{ij\sigma }><a_{i^{}\sigma }^{}(t)a_{j\sigma }>$$

This approximation was made in the spirit of the approximation of ”two interacting modes” and means ignoring the renormalization of the vertex, i.e., the correlation in the propagation of an electron (hole) and the propagation of charge density fluctuations.
Writing down further spectral representation for the correlation functions in (216), we represent the mass operator by the sum

(217)
$$\widehat{M}_{ii^{}}(\omega )=\widehat{M}_{ii^{}}^1(\omega )+\widehat{M}_{ii^{}}^2(\omega )$$

The first contribution $`M^1`$ has a form characteristic of an interacting electron-phonon system

$`M_{ii^{}}^1(\omega )={\displaystyle \underset{nn^{}}{}}{\displaystyle \underset{jj^{}}{}}V_{ijn}V_{j^{}i^{}n^{}}{\displaystyle \frac{1}{2}}{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\omega _1d\omega _2}{\omega \omega _1\omega _2}}(\mathrm{cot}{\displaystyle \frac{\beta \omega _1}{2}}+\mathrm{tan}{\displaystyle \frac{\beta \omega _2}{2}})`$
(218) $`({\displaystyle \frac{1}{\pi }}Im<<u_n|u_n^{}>>_{\omega _2})({\displaystyle \frac{1}{\pi }}\tau _3Im<<\psi _j|\psi _j^{}^{}>>_{\omega _1}\tau _3)`$
The contribution $`M_{ii^{}}^2`$ has a more complicated structure

(219)
$$M_{ii^{}}^2=\frac{U^2}{2}_{\mathrm{}}^+\mathrm{}\frac{d\omega _1d\omega _2}{\omega \omega _1\omega _2}(\mathrm{cot}\frac{\beta \omega _1}{2}+\mathrm{tan}\frac{\beta \omega _2}{2})\left(\begin{array}{cc}m_{11}& m_{12}\\ m_{21}& m_{22}\end{array}\right)$$

where

$`m_{11}=({\displaystyle \frac{1}{\pi }}Im<<n_i|n_i^{}>>_{\omega _2})({\displaystyle \frac{1}{\pi }}Im<<a_i|a_i^{}^{}>>_{\omega _1})`$
$`m_{12}=({\displaystyle \frac{1}{\pi }}Im<<n_i|n_i^{}>>_{\omega _2})({\displaystyle \frac{1}{\pi }}Im<<a_i|a_i^{}^{}>>_{\omega _1})`$
$`m_{21}=({\displaystyle \frac{1}{\pi }}Im<<n_i|n_i^{}>>_{\omega _2})({\displaystyle \frac{1}{\pi }}Im<<a_i|a_i^{}^{}>>_{\omega _1})`$
$`m_{22}=({\displaystyle \frac{1}{\pi }}Im<<n_i|n_i^{}>>_{\omega _2})({\displaystyle \frac{1}{\pi }}Im<<a_i|a_i^{}^{}>>_{\omega _1})`$
The equations (212) and (217) constitute a self-consistent system of equations for the single-particle GF of the Hubbard model on a vibrating lattice. Note that these equations of superconductivity can be in an obvious way transformed to the standard form of the Eliashberg equations . The numerical calculations of the electron-phonon spectral function $`\alpha ^2(\omega )F(\omega )`$ for a few transition metals were done in ref. . It is worthy to emphasize that in paper  a very detailed microscopic theory of the strong coupling superconductivity in highly disordered transition metal alloys was developed on the basis of the IGF method within the MTBA reformulated approach . The Eliashberg-type strong coupling equations for highly disordered alloys were derived. It was shown that the electron-phonon spectral function in alloys is modified strongly. Thus, the self-consistent system of superconductivity equations obtained in the Wannier representation makes it possible to investigate real transition metals, their alloys, and compounds from a unified point of view.

## 9 Quasi-Particle Dynamics of Anderson Models

### 9.1 Quasi-Particle Dynamics of SIAM

In this Section, we consider the many-body quasi-particle dynamics of the Anderson impurity model at finite temperatures in the framework of the equation-of-motion method. In spite of many theoretical efforts, there is no complete solution of the dynamic problem for the ”simple” Anderson/Hubbard model. One of the main reasons for this is that it has been recognized relatively recently only that the simplicity of the Anderson model manifests itself not in the many-body dynamics ( the right definition of quasi-particles via the poles of GF; see Section 6.1 ) but rather at quite a different level - in the dynamics of two-particle scattering, resulting in the elegant Bethe-ansatz solution (for the relativistic spectrum linear in $`k`$ ), which gives the static characteristics (static susceptibility, specific heat, etc). In this sense, as to the true many-body dynamics, the complete analytic solution of this problem is still quite an open subject. This Section is primarily devoted to the analysis of the relevant many-body dynamic solution of the SIAM and its correct functional structure. We wish to determine which solution actually arises both from the self-consistent many-body approach and intrinsic nature of the model itself. We believe highly that before numerical calculations of the spectral intensity of the Green function at low energy and low temperature it is quite important to have a consistent and closed analytic representation for the one-particle GF of the SIAM and Hubbard model. The paper  clearly shows the importance of the calculation of the GF and spectral densities for SIAM in a self-consistent way. An alternative approach to dynamics of the Anderson model was formulated within a modified version of the Kadanoff-Baym method . Unfortunately, the Neal approach also have certain drawbacks.
A proper many-body description of dynamic correlations is very actual also for the investigation of the dynamics of the many-impurity Anderson model, where standard advanced many-body methods do not work properly in usual formulation. Recently, a lot of efforts were devoted to a better understanding of the static and dynamic properties of the Anderson model in the context of many-impurity case . This field is quite important for the description of magnetic properties of anomalous rare-earth compounds,. The problem of an adequate and consistent description of dynamics of single-impurity and many-impurity Anderson models ( SIAM and MIAM ) and other models of correlated lattice electrons was not yet solved analytically completely . During the last decades, a lot of theoretical papers were published, attacking the Anderson model by many refined many-body analytic methods. Nevertheless, a fully consistent dynamic analytic solution in the closed form for a single-particle propagator of SIAM is still lacking. In this Section, the problem of consistent analytic description of the many-body dynamics of SIAM is discussed in the framework of the equation-of-motion approach for double-time thermodynamic GFs. In addition to the IGF approach, we find a new exact identity relating the one-particle and many-particle GFs. Using this identity, we present a consistent and general scheme for construction of generalized solutions of SIAM. A new approach for the complex expansion for the single-particle propagator in terms of Coulomb repulsion $`U`$ and hybridization $`V`$ is proposed. Using the exact identity, an essentially new many-body dynamic solution of SIAM is derived. This approach offers a new way for the systematic construction of approximate interpolation dynamic solutions of strongly correlated electron systems.

### 9.2 IGF Approach to SIAM

After discussing some of the basic facts about the correct functional structure of the relevant dynamic solution of correlated electron models we are looking for, described in previous Sections, we give a similar consideration for SIAM. It was shown in , using the minimal algebra of relevant operators, that the construction of the GMFs for SIAM is quite nontrivial for the strongly correlated case, and it is rather difficult to get it from an intuitive physical point of view. Let us consider first the following matrix GF

(220)
$$\widehat{G}(\omega )=\left(\begin{array}{cc}<<c_{k\sigma }|c_{k\sigma }^{}>>& <<c_{k\sigma }|f_{0\sigma }^{}>>\\ <<f_{0\sigma }|c_{k\sigma }^{}>>& <<f_{0\sigma }|f_{0\sigma }^{}>>\end{array}\right)$$

Performing the first-time differentiation and defining the irreducible GF

(221) $`(^{(ir)}<<f_{0\sigma }f_{0\sigma }^{}f_{0\sigma }|f_{0\sigma }^{}>>_\omega )=<<f_{0\sigma }f_{0\sigma }^{}f_{0\sigma }|f_{0\sigma }^{}>>_\omega `$
$`<n_{0\sigma }><<f_{0\sigma }|f_{0\sigma }^{}>>_\omega `$
we obtain the following equation of motion in the matrix form

(222)
$$\underset{p}{}\widehat{F}_p(\omega )\widehat{G}_p(\omega )=\widehat{1}+U\widehat{D}^{(ir)}(\omega )$$

where all definitions are rather evident. Proceeding further with the IGF technique, the equation of motion (222) is exactly rewritten in the form of the Dyson equation

(223)
$$\widehat{G}(\omega )=\widehat{G}^0(\omega )+\widehat{G}^0(\omega )\widehat{M}(\omega )\widehat{G}(\omega )$$

The generalized mean field GF $`G^0`$ is defined by

(224)
$$\underset{p}{}F_p(\omega )G_p^0(\omega )=\widehat{I}$$

The explicit solutions for diagonal elements of $`G^0`$ are

(225) $`<<f_{0\sigma }|f_{0\sigma }^{}>>_\omega ^0=(\omega E_{0\sigma }Un_\sigma S(\omega )))^1`$
(226) $`<<c_{k\sigma }|c_{k\sigma }^{}>>_\omega ^0=\left(\omega ϵ_k{\displaystyle \frac{|V_k|^2}{\omega E_{0\sigma }Un_\sigma }}\right)^1`$
where

(227)
$$S(\omega )=\underset{k}{}\frac{|V_k|^2}{\omega ϵ_k}$$

The mass operator, which describes inelastic scattering processes, has the following matrix form

(228)
$$\widehat{M}(\omega )=\left(\begin{array}{cc}0& 0\\ 0& M_{0\sigma }\end{array}\right)$$

where

(229)
$$M_{0\sigma }=U^2(^{(ir)}<<f_{0\sigma }n_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>_\omega ^{(ir)})^{(p)}$$

From the formal solution of the Dyson equation (38) one obtains

(230) $`<<f_{0\sigma }|f_{0\sigma }^{}>>_\omega =\left(\omega E_{0\sigma }Un_\sigma M_{0\sigma }S(\omega )\right)^1`$
(231) $`<<c_{k\sigma }|c_{k\sigma }^{}>>_\omega =\left(\omega ϵ_k{\displaystyle \frac{|V_k|^2}{\omega E_{0\sigma }Un_\sigma M_{0\sigma }}}\right)^1`$
To calculate the self-energy in a self-consistent way, we have to approximate it by lower-order GFs. Let us start by analogy with the Hubbard model with a pair-type approximation (8.1)

(232) $`M_{0\sigma }(\omega )=`$
$`U^2{\displaystyle \frac{d\omega _1d\omega _2d\omega _3}{\omega +\omega _1\omega _2\omega _3}}`$
$`[n(\omega _2)n(\omega _3)+n(\omega _1)(1n(\omega _2)n(\omega _3))]g_{0\sigma }(\omega _1)g_{0\sigma }(\omega _2)g_{0\sigma }(\omega _3)`$
where we used the notation

$$g_{0\sigma }(\omega )=\frac{1}{\pi }Im<<f_{0\sigma }|f_{0\sigma }^{}>>_\omega $$

The equations (223) and (232) constitute a closed self-consistent system of equations for the single-electron GF for SIAM model, but only for weakly correlated limit. In principle, we can use, on the r.h.s. of (232), any workable first iteration-step form of the GF and find a solution by repeated iteration. If we take for the first iteration step the expression

(233)
$$g_{0\sigma }(\omega )\delta (\omega E_{0\sigma }Un_\sigma ),$$

we get, for the self-energy, the explicit expression

(234)
$$M_{0\sigma }(\omega )=U^2\frac{n(E_{0\sigma }+Un_\sigma )(1n(E_{0\sigma }+Un_\sigma ))}{\omega E_{0\sigma }Un_\sigma }=U^2N_\sigma (1N_\sigma )G_\sigma ^0(\omega )$$

where $`N_\sigma =n(E_{0\sigma }+Un_\sigma )`$. This is the well-known ”atomic” limit of the self-energy.
Let us try again another type of the approximation for $`M`$. The approximation which we will use reflects the interference between the one-particle branch and the collective one

$`<f_{0\sigma }(t)f_{0\sigma }^{}(t)f_{0\sigma }(t)f_{0\sigma }^{}f_{0\sigma }f_{0\sigma }^{}>^{(ir)}`$
$`<f_{0\sigma }^{}(t)f_{0\sigma }><n_{0\sigma }(t)n_{0\sigma }>+`$
$`<f_{0\sigma }^{}(t)f_{0\sigma }><f_{0\sigma }^{}(t)f_{0\sigma }(t)f_{0\sigma }^{}f_{0\sigma }>+`$
(235) $`<f_{0\sigma }^{}(t)f_{0\sigma }><f_{0\sigma }(t)f_{0\sigma }(t)f_{0\sigma }^{}f_{0\sigma }^{}>`$
If we retain only the first term in (9.2) and make use of the same iteration as in (233), we obtain

(236)
$$M_{0\sigma }(\omega )U^2\frac{(1n(E_{0\sigma }+Un_\sigma ))}{\omega E_{0\sigma }Un_\sigma }<n_{0\sigma }n_{0\sigma }>$$

If we retain the second term in (9.2), we obtain

$`M_{0\sigma }(\omega )=U^2{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\omega _1𝑑\omega _2{\displaystyle \frac{1+N(\omega _1)n(\omega _2)}{\omega \omega _1\omega _2}}`$
$`({\displaystyle \frac{1}{\pi }}Im<<S_0^\pm |S_0^{}>>_{\omega _1})`$
(237) $`\left({\displaystyle \frac{1}{\pi }}Im<<f_{0\sigma }|f_{0\sigma }^{}>>_{\omega _2}\right)`$
where the following notation was been used:

$$S_0^+=f_0^{}f_0;S_0^{}=f_0^{}f_0$$

It is possible to rewrite (9.2) in a more convenient way now

$`M_{0\sigma }(\omega )=U^2{\displaystyle 𝑑\omega ^{}(\mathrm{cot}\frac{\omega \omega ^{}}{2T}+\mathrm{tan}\frac{\omega ^{}}{2T})}`$
(238) $`\left({\displaystyle \frac{1}{\pi }}Im\chi ^\pm (\omega \omega ^{})g_{0\sigma }(\omega ^{})\right)`$
The equations (223) and (9.2) constitute a self-consistent system of equations for the single-particle GF of SIAM. Note that spin-up and spin-down electrons are correlated when they occupy the impurity level. So, this really improves the H-F theory in which just these correlations were missed. The role of electron-electron correlation becomes much more crucial for the case of strong correlation.

### 9.3 SIAM. Strong Correlation

The simplest relevant algebra of the operators used for the description of the strong correlation has a similar form as for that of the Hubbard model (8.2). Let us represent the matrix GF (220) in the following form

(239)
$$\widehat{G}(\omega )=\underset{\alpha \beta }{}\left(\begin{array}{cc}<<c_{k\sigma }|c_{k\sigma }^{}>>& <<c_{k\sigma }|d_{0\beta \sigma }^{}>>\\ <<d_{0\alpha \sigma }|c_{k\sigma }^{}>>& <<d_{0\alpha \sigma }|d_{0\beta \sigma }^{}>>\end{array}\right)$$

Then we proceed by analogy with the calculations for the Hubbard model. The equation of motion for the auxiliary matrix GF

(240) $`\widehat{F}_\sigma (\omega )=`$
$`\left(\begin{array}{ccc}<<c_{k\sigma }|c_{k\sigma }^{}>>& <<c_{k\sigma }|d_{0+\sigma }^{}>>& <<c_{k\sigma }|d_{0\sigma }^{}>>\\ <<d_{0+\sigma }|c_{k\sigma }^{}>>& <<d_{0+\sigma }|d_{0+\sigma }^{}>>& <<d_{0+\sigma }|d_{0\sigma }^{}>>\\ <<d_{0\sigma }|c_{k\sigma }^{}>>& <<d_{0\sigma }|d_{0+\sigma }^{}>>& <<d_{0\sigma }|d_{0\sigma }^{}>>\end{array}\right)`$
is of the following form

(241)
$$\widehat{E}\widehat{F}_\sigma (\omega )\widehat{I}=\widehat{D}$$

where the following matrix notation was used

(242) $`\widehat{E}=\left(\begin{array}{ccc}(\omega ϵ_k)& V_k& V_k\\ 0& (\omega E_{0\sigma }U_+)& 0\\ 0& 0& (\omega E_{0\sigma }U_{})\end{array}\right)`$
$`\widehat{I}=\left(\begin{array}{ccc}1& 0& 0\\ 0& n_{0\sigma }^+& 0\\ 0& 0& n_{0\sigma }^{}\end{array}\right).`$
$$U_\alpha =\{\begin{array}{cc}U,\hfill & \alpha =+\hfill \\ 0,\hfill & \alpha =\hfill \end{array}$$

Here $`\widehat{D}`$ is a higher-order GF, with the following structure

(243)
$$\widehat{D}(\omega )=\left(\begin{array}{ccc}0& 0& 0\\ D_{21}& D_{22}& D_{23}\\ D_{31}& D_{32}& D_{33}\end{array}\right)$$

In accordance with the general method of Section 3, we by define the matrix IGF:

(244)
$$\widehat{D}^{(ir)}(\omega )=\widehat{D}\underset{\alpha }{}\left(\genfrac{}{}{0pt}{}{A^{+\alpha }}{A^\alpha }\right)(G_\sigma ^{\alpha +}G_\sigma ^\alpha )$$

Here the notation was used:

(245) $`A^{++}={\displaystyle \frac{<(f_{0\sigma }^{}c_{p\sigma }+c_{p\sigma }^{}f_{0\sigma })(n_{0\sigma }n_{0\sigma })>}{<n_{0\sigma }>}}`$
(246) $`A^{}={\displaystyle \frac{<(f_{0\sigma }^{}c_{p\sigma }+c_{p\sigma }^{}f_{0\sigma })(1+n_{0\sigma }n_{0\sigma })>}{<1n_{0\sigma }>}}`$
$`A^+=A^{++},A^+=A^{}`$
The generalized mean-field GF is defined by

(247)
$$\widehat{E}\widehat{F}_\sigma ^0(\omega )\widehat{I}=0;G^0=\underset{\alpha \beta }{}F_{\alpha \beta }^0$$

From the last definition we find that

$`<<f_{0\sigma }|f_{0\sigma }^{}>>_\omega ^0={\displaystyle \frac{<n_{0\sigma }>}{\omega E_{0\sigma }U_+_pV_pA^{++}}}\left(1+{\displaystyle \frac{_pV_pA^+}{\omega E_{0\sigma }U_{}}}\right)`$
(248) $`+{\displaystyle \frac{1<n_{0\sigma }>}{\omega E_{0\sigma }U_{}_pV_pA^{}}}\left(1+{\displaystyle \frac{_pV_pA^+}{\omega E_{0\sigma }U_+}}\right)`$
(249) $`<<c_{k\sigma }|c_{k\sigma }^{}>>_\omega ^0=\left(\omega ϵ_k|V_k|^2F^{at}(\omega )\right)^1`$
where

(250)
$$F^{at}=\frac{<n_{0\sigma }>}{\omega E_{0\sigma }U_+}+\frac{1<n_{0\sigma }>}{\omega E_{0\sigma }U_{}}$$

For $`V_p=0`$, we obtain, from solution (9.3), the atomic solution $`F^{at}`$. The conduction electron GF (249) also gives a correct expression for $`V_k=0`$.

### 9.4 IGF Method and Interpolation Solution

To show explicitly the flexibility of the IGF method, we consider a more extended new algebra of operators from which the relevant matrix GF should be constructed to make the connection with the interpolation solution of the Anderson model. For this aim, let us consider the following equation of motion in the matrix form

(251)
$$\underset{p}{}F(p,k)G_{p\sigma }(\omega )=I+\underset{p}{}V_pD_p(\omega )$$

where $`G`$ is the initial $`4\times 4`$ matrix GF and $`D`$ is the higher-order GF:

(252)
$$G_\sigma =\left(\begin{array}{cccc}G_{11}& G_{12}& G_{13}& G_{14}\\ G_{21}& G_{22}& G_{23}& G_{24}\\ G_{31}& G_{32}& G_{33}& G_{34}\\ G_{41}& G_{42}& G_{43}& G_{44}\end{array}\right)$$

Here the following notation was used

$`G_{11}=<<c_{k\sigma }|c_{k\sigma }^{}>>;G_{12}=<<c_{k\sigma }|f_{0\sigma }^{}>>;`$
$`G_{13}=<<c_{k\sigma }|f_{0\sigma }^{}n_{0\sigma }>>;G_{14}=<<c_{k\sigma }|c_{k\sigma }^{}n_{0\sigma }>>;`$
$`G_{21}=<<f_{0\sigma }|c_{k\sigma }^{}>>;G_{22}=<<f_{0\sigma }|f_{0\sigma }^{}>>;`$
(253) $`G_{23}=<<f_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>;G_{24}=<<f_{0\sigma }|c_{k\sigma }^{}n_{0\sigma }>>;`$
$`G_{31}=<<f_{0\sigma }n_{0\sigma }|c_{k\sigma }^{}>>;G_{32}=<<f_{0\sigma }n_{0\sigma }|f_{0\sigma }^{}>>;`$
$`G_{33}=<<f_{0\sigma }n_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>;G_{34}=<<f_{0\sigma }n_{0\sigma }|c_{k\sigma }^{}n_{0\sigma }>>;`$
$`G_{41}=<<c_{k\sigma }n_{0\sigma }|c_{k\sigma }^{}>>;G_{42}=<<c_{k\sigma }n_{0\sigma }|f_{0\sigma }^{}>>;`$
$`G_{43}=<<c_{k\sigma }n_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>;G_{44}=<<c_{k\sigma }n_{0\sigma }|c_{k\sigma }^{}n_{0\sigma }>>;`$
We avoid to write down explicitly the relevant 16 GFs, of which the matrix GF $`D`$ consist, for the brevity. For our aims, it is enough to proceed forth in the following way.
The equation (251) results from the first-time differentiation of the GF $`G`$ and is a starting point for the IGF approach. Let us introduce the irreducible part for the higher-order GF $`D`$ in the following way

(254)
$$D_\beta ^{(ir)}=D_\beta \underset{\alpha }{}L^{\beta \alpha }G_{\alpha \beta };(\alpha ,\beta )=(1,2,3,4)$$

and define the GMF GF according to

(255)
$$\underset{p}{}\stackrel{~}{F}(p,k)G_{p\sigma }^{MF}(\omega )=I,$$

Then, we are able to write down explicitly the Dyson equation (37) and the exact expression for the self-energy M in the matrix form:

(256)
$$M_{k\sigma }(\omega )=I^1\underset{p,q}{}V_pV_q\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& M_{33}& M_{34}\\ 0& 0& M_{43}& M_{44}\end{array}\right)I^1$$

Here the matrix $`I`$ is given by

$$I=\left(\begin{array}{cccc}1& 0& 0& <n_{0\sigma }>\\ 0& 1& <n_{0\sigma }>& 0\\ 0& <n_{0\sigma }>& <n_{0\sigma }>& 0\\ <n_{0\sigma }>& 0& 0& <n_{0\sigma }>\end{array}\right)$$

and the the matrix elements of M are of the form:

$`M_{33}=(<<A_1^{(ir)}(p)|B_1^{(ir)}(q)>>)^{(p)},M_{34}=(<<A_1^{(ir)}(p)|B_2^{(ir)}(k,q)>>)^{(p)}`$
$`M_{43}=(<<A_2^{(ir)}(k,p)|B_1^{(ir)}(q)>>)^{(p)},M_{44}=(<<A_2^{(ir)}(k,p)|B_2^{(ir)}(k,q)>>)^{(p)}`$
where

$`A_1(p)=(c_{p\sigma }^{}f_{0\sigma }f_{0\sigma }c_{p\sigma }f_{0\sigma }^{}f_{0\sigma });`$
$`A_2(k,p)=(c_{k\sigma }f_{0\sigma }^{}c_{p\sigma }c_{k\sigma }c_{p\sigma }^{}f_{0\sigma });`$
(258) $`B_1(p)=(f_{0\sigma }^{}c_{p\sigma }^{}f_{0\sigma }f_{0\sigma }^{}f_{0\sigma }^{}c_{p\sigma });`$
$`B_2(k,p)=(c_{k\sigma }^{}c_{p\sigma }^{}f_{0\sigma }c_{k\sigma }^{}f_{0\sigma }^{}c_{p\sigma });`$
Since the self-energy $`M`$ describes the processes of inelastic scattering of electrons (c-c , f-f, and c-f types), its approximate representation would be defined by the nature of physical assumptions about this scattering.
To get an idea about the functional structure of our GMF solution (255), let us write down the matrix element $`G_{33}^{MF}`$:

$`G_{33}^{MF}=<<f_{0\sigma }n_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>=`$
$`{\displaystyle \frac{<n_{0\sigma }>}{\omega ϵ_f^{MF}US^{MF}(\omega )Y(\omega )}}+`$
(259) $`{\displaystyle \frac{<n_{0\sigma }>Z(\omega )}{(\omega ϵ_f^{MF}US^{MF}(\omega )Y(\omega ))(\omega E_{0\sigma }S(\omega ))}}`$
(260) $`Y(\omega )={\displaystyle \frac{UZ(\omega )}{\omega E_{0\sigma }S(\omega )}}`$
$`Z(\omega )=S(\omega ){\displaystyle \underset{p}{}}{\displaystyle \frac{V_PL^{41}}{\omega ϵ_p^{MF}}}+`$
(261) $`{\displaystyle \underset{p}{}}{\displaystyle \frac{|V_p|^2L^{42}}{\omega ϵ_p^{MF}}}+S(\omega )L^{31}+{\displaystyle \underset{p}{}}V_pL^{32}`$
Here the coefficients $`L^{41},L^{42},L^{31}`$, and $`L^{32}`$ are certain complicated averages (see definition (254)) from which the functional of the GMF is build. To clarify the functional structure of the obtained solution, let us consider our first equation of motion (251) , before introducing the irreducible GFs (254). Let us put simply, in this equation, the higher-order GF $`D=0`$! To distinguish this simplest equation from the GMF one (255), we write it in the following form

(262)
$$\underset{p}{}F(p,k)G^0(p,\omega )=I$$

The corresponding matrix elements which we are interested in here read

(263) $`G_{22}^0=<<f_{0\sigma }|f_{0\sigma }^{}>>=`$
$`{\displaystyle \frac{1<n_{0\sigma }>}{\omega E_{0\sigma }S(\omega )}}+{\displaystyle \frac{<n_{0\sigma }>}{\omega E_{0\sigma }S(\omega )U}}`$
(264) $`G_{33}^0=<<f_{0\sigma }n_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>={\displaystyle \frac{<n_{0\sigma }>}{\omega E_{0\sigma }S(\omega )U}}`$
(265) $`G_{32}^0=<<f_{0\sigma }n_{0\sigma }|f_{0\sigma }^{}>>=G_{33}^0`$
The conclusion is rather evident. The simplest interpolation solution follows from our matrix GF (252) in the lowest order in $`V`$, even before introduction of GMF corrections, not speaking about the self-energy corrections. The two GFs $`G_{32}^0`$ and $`G_{33}^0`$ are equal only in the lowest order in $`V`$. It is quite clear that our full solution (38) that includes the self-energy corrections is much more richer.
It is worthwhile to stress that our $`4\times 4`$ matrix GMF GF (252) gives only approximate description of suitable mean fields. If we consider more extended algebra, we get the more correct structure of the relevant GMF.

### 9.5 Dynamic Properties of SIAM

To demonstrate clearly the advantages of the IGF method for SIAM, it is worthwhile to emphasize a few important points about the approach based on the equations-of-motion for the GFs. To give a more instructive discussion, let us consider the single-particle GF of localized electrons $`G_\sigma =<<f_{0\sigma }|f_{0\sigma }^{}>>`$. The simplest approximate ”interpolation” solution of SIAM is of the form:

$`G_\sigma (\omega )={\displaystyle \frac{1}{\omega E_{0\sigma }S(\omega )}}+{\displaystyle \frac{U<n_{0\sigma }>}{(\omega E_{0\sigma }S(\omega )U)(\omega E_{0\sigma }S(\omega ))}}=`$
(266) $`{\displaystyle \frac{1<n_{0\sigma }>}{\omega E_{0\sigma }S(\omega )}}+{\displaystyle \frac{<n_{0\sigma }>}{\omega E_{0\sigma }S(\omega )U}}`$
The values of $`n_\sigma `$ are determined through the self-consistency equation

(267)
$$n_\sigma =<n_{0\sigma }>=\frac{1}{\pi }𝑑En(E)ImG_\sigma (E,n_\sigma )$$

The ”atomic-like” interpolation solution (9.5) reproduces correctly the two limits:

$`G_\sigma (\omega )={\displaystyle \frac{1<n_{0\sigma }>}{\omega E_{0\sigma }}}+{\displaystyle \frac{<n_{0\sigma }>}{\omega E_{0\sigma }U}},forV=0`$
(268) $`G_\sigma (\omega )={\displaystyle \frac{1}{\omega E_{0\sigma }S(\omega )}},forU=0`$
The important point about formulas (9.5) is that any approximate solution of SIAM should be consistent with it. Let us remind how to get solution (9.5). It follows from the system of equations for small-$`V`$ limit:

$`(\omega E_{0\sigma }S(\omega ))<<f_{0\sigma }|f_{0\sigma }^{}>>_\omega =1+U<<f_{0\sigma }n_{0\sigma }|f_{0\sigma }^{}>>_\omega ,`$
$`(\omega E_{0\sigma }U)<<f_{0\sigma }n_{0\sigma }|f_{0\sigma }^{}>>_\omega `$
(269) $`<n_{0\sigma }>+{\displaystyle \underset{k}{}}V_k<<c_{k\sigma }n_{0\sigma }|f_{0\sigma }^{}>>_\omega ,`$
$`(\omega ϵ_k)<<c_{k\sigma }n_{0\sigma }|f_{0\sigma }^{}>>_\omega =`$
$`V_k<<f_{0\sigma }n_{0\sigma }|f_{0\sigma }^{}>>_\omega `$
The equations (9.5) are approximate; they include two more terms, treated in the limit of small $`V`$ in paper .
We now proceed further. In paper  the GF $`G`$ was calculated in the limit of infinitely strong Coulomb correlation $`U`$ and for small hybridization $`V`$. The functional structure of the Lacroix solution generalizes the solution (9.5). The starting point is the system of equations:

(270) $`(\omega E_{0\sigma }S(\omega ))<<f_{0\sigma }|f_{0\sigma }^{}>>=1+U<<f_{0\sigma }n_{0\sigma }|f_{0\sigma }^{}>>`$
$`(\omega E_{0\sigma }U)<<f_{0\sigma }n_{0\sigma }|f_{0\sigma }^{}>>=<n_{0\sigma }>+{\displaystyle \underset{k}{}}V_k(<<c_{k\sigma }n_{0\sigma }|f_{0\sigma }^{}>>`$
(271) $`<<c_{k\sigma }f_{0\sigma }^{}f_{0\sigma }|f_{0\sigma }^{}>>+<<c_{k\sigma }^{}f_{0\sigma }f_{0\sigma }|f_{0\sigma }^{}>>)`$
Using a relatively simple decoupling procedure for a higher-order equation of motion, a qualitatively correct low-temperature spectral intensity was calculated. The final expression for $`G`$ for finite $`U`$ is of the form

$`<<f_{0\sigma }|f_{0\sigma }^{}>>={\displaystyle \frac{1}{\omega E_{0\sigma }S(\omega )+US_1(\omega )}}+`$
(272) $`{\displaystyle \frac{U<n_{0\sigma }>+UF_1(\omega )}{K(\omega )(\omega E_{0\sigma }S(\omega )+US_1(\omega ))}}`$
where $`F_1`$, $`S_1`$, and $`K`$ are certain complicated expressions. We write down explicitly the infinite $`U`$ approximate GF :

(273)
$$<<f_{0\sigma }|f_{0\sigma }^{}>>=\frac{1<n_{0\sigma }>F_\sigma (\omega )}{\omega E_{0\sigma }S(\omega )Z_\sigma ^1(\omega )}$$

The following notation was used

(274) $`F_\sigma =V{\displaystyle \underset{k}{}}{\displaystyle \frac{<f_{0\sigma }^{}c_{k\sigma }>}{\omega ϵ_k}}`$
(275) $`Z_\sigma ^1=V^2{\displaystyle \underset{q,k}{}}{\displaystyle \frac{<c_{q\sigma }^{}c_{k\sigma }>}{\omega ϵ_k}}S(\omega )V{\displaystyle \underset{k}{}}{\displaystyle \frac{<f_{0\sigma }^{}c_{k\sigma }>}{\omega ϵ_k}}`$
The functional structure of the single-particle GF (9.5) is quite transparent. The expression in the numerator of (9.5) plays the role of ”dynamic mean field”, proportional to $`<f_{0\sigma }^{}c_{k\sigma }>`$. In the denominator, instead of bare shift $`S(\omega )`$ (227) we have an ”effective shift” $`S^1=S(\omega )+Z_\sigma ^1(\omega )`$. The choice of the specific procedure of decoupling for the higher-order equation of motion specifies the selected ”generalized mean fields” (GMFs) and ”effective shifts”.

### 9.6 Interpolation Solutions of Correlated Models

It is to the point to discuss briefly the general concepts of construction of an interpolation dynamic solution of the strongly correlated electron models. The very problem of the consistent interpolation solutions of the many-body electron models was formulated explicitly by Hubbard in the context of the Hubbard model. Hubbard clearly pointed out one particular feature of consistent theory, insisting that it should give exact results in the two opposite limits of very wide and very narrow bands. The functional structure of a required interpolation solution can be clarified if one considers the atomic (very narrow band) solution of the Hubbard model (47):

(276)
$$G^{at}(\omega )=\frac{1n_\sigma }{\omega t_0}+\frac{n_\sigma }{\omega t_0U}=\frac{1}{\omega t_0\mathrm{\Sigma }^{at}(\omega )}$$

where

(277)
$$\mathrm{\Sigma }^{at}(\omega )=\frac{n_\sigma U}{1\frac{(1n_\sigma )U}{\omega t_0}};t_0=t_{ii}$$

Let us consider the expansion in terms of $`U`$:

(278)
$$\mathrm{\Sigma }^{at}(\omega )n_\sigma U+n_\sigma (1n_\sigma )U^2\frac{1}{\omega t_0}+O(U)$$

The ”Hubbard I” solution (52)can be written as

(279)
$$G_k=\frac{1}{\omega ϵ(k)\mathrm{\Sigma }^{at}(\omega )}=\frac{1}{(G^{at})^1+t_0ϵ(k)}$$

The partial ”Hubbard III” solution, called the ”alloy analogy” approximation is of the form:

(280)
$$\mathrm{\Sigma }(\omega )=\frac{n_\sigma U}{1(U\mathrm{\Sigma }(\omega ))G(\omega )}$$

Equation (280) follows from (277) when one takes into account the following relationship:

(281)
$$\frac{1}{\omega t_0}\frac{1}{1n_\sigma }G(\omega )\mathrm{\Sigma }(\omega )G(\omega )$$

The Coherent Potential Approximation (CPA) provides the basis for physical interpretation of equation (280) which corresponds to elimination of the dynamics of $`\sigma `$ electrons. By analogy with (278), it is possible to expand:

(282)
$$\frac{n_\sigma U}{1(U\mathrm{\Sigma }(\omega ))G(\omega )}n_\sigma U+n_\sigma U(U\mathrm{\Sigma })G^0(\omega \mathrm{\Sigma })+O(U)$$

The solution (273) does not reproduce correctly the U-perturbation expansion for the self-energy:

$`M_\sigma (\omega )U<n_{0\sigma }>+`$
$`U^2{\displaystyle 𝑑E_1𝑑E_2𝑑E_3\frac{n(E_1)n(E_2)(1n(E_3))+(1n(E_1)(1n(E_2))n(E_3)}{\omega E_1E_2+E_3}}`$
(283) $`ImG_\sigma (E_1)ImG_\sigma (E_2)ImG_\sigma (E_3)`$
It can be shown that it is possible , in principle, to find a certain way to incorporate this $`U^2`$ perturbation theory expansion into the functional structure of an interpolation dynamic solution of SIAM in a self-consistent way within the higher-order GFs . The IGF approach with the use of minimal algebra of relevant operators allows one to find an interpolation solution for weak and strong Coulomb interaction $`U`$ and to calculate explicitly the quasi-particle spectra and their damping for both the limits. The U-perturbation expansion (143) is included into the IGF scheme in a self-consistent way. The correct second-order contribution to the local approximation for the Hubbard model is of the form

(284)
$$\stackrel{~}{G}_\sigma \frac{G_\sigma <<n_{0\sigma }|n_{0\sigma }>>}{n_\sigma (1n_\sigma )}$$

The same arguments are also valid for SIAM.

### 9.7 Complex Expansion for a Propagator

We now proceed with analytic many-body consideration. One can attempt to consider a suitable solution for the SIAM starting from the following exact relation derived in paper :

(285) $`<<f_{0\sigma }|f_{0\sigma }^{}>>=g^0+g^0Pg^0`$
(286) $`g^0=(\omega E_{0\sigma }S(\omega ))^1`$
(287) $`P=U<n_{0\sigma }>+U^2<<f_{0\sigma }n_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>`$
The advantage of the equation (285) is that it is a pure identity and does not include any approximation. If we insert our GMF solution (273) into (285), we get an essentially new dynamic solution of SIAM constructed on the basis of the complex (combined) expansion of the propagator in both $`U`$ and $`V`$ parameters and reproducing exact solutions of SIAM for $`V=0`$ and $`U=0`$. It generalizes (even on the mean-field level) the solutions of papers , .
Having emphasized the importance of the role of equation (285) , let us see now what is the best possible fit for higher-order GF in (287). We consider the equation of motion for it:

(288) $`(\omega E_{0\sigma }U)<<f_{0\sigma }n_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>=<n_{0\sigma }>+`$
$`{\displaystyle \underset{k}{}}V_k(<<c_{k\sigma }n_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>+`$
$`<<c_{k\sigma }^{}f_{0\sigma }f_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>><<c_{k\sigma }f_{0\sigma }^{}f_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>)`$
We can think of it as defining new kinds of elastic and inelastic scattering processes that contribute to the formation of generalized mean fields and self-energy (damping) corrections. The construction of suitable mean fields can be quite nontrivial, and to describe these contributions self-consistently, let us consider the equations of motion for higher-order GFs in the r.h.s. of (288)

$`(\omega ϵ_k)<<c_{k\sigma }n_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>=`$
(289) $`V<<f_{0\sigma }n_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>+`$
$`{\displaystyle \underset{p}{}}V(<<c_{k\sigma }f_{0\sigma }^{}c_{p\sigma }|f_{0\sigma }^{}n_{0\sigma }>><<c_{k\sigma }c_{p\sigma }^{}f_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>)`$
$`(\omega ϵ_kE_{0\sigma }+E_{0\sigma })<<c_{k\sigma }f_{0\sigma }^{}f_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>`$
(290) $`=<f_{0\sigma }^{}c_{k\sigma }n_{0\sigma }>`$
$`V<<f_{0\sigma }n_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>+`$
$`{\displaystyle \underset{p}{}}V(<<c_{k\sigma }f_{0\sigma }^{}c_{p\sigma }|f_{0\sigma }^{}n_{0\sigma }>><<c_{k\sigma }c_{p\sigma }^{}f_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>)`$
(291) $`(\omega +ϵ_kE_{0\sigma }E_{0\sigma }U)<<c_{k\sigma }^{}f_{0\sigma }f_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>`$
$`=<c_{k\sigma }^{}f_{0\sigma }f_{0\sigma }^{}f_{0\sigma }>+`$
$`V<<f_{0\sigma }n_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>+`$
$`{\displaystyle \underset{p}{}}V(<<c_{k\sigma }^{}c_{p\sigma }f_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>+<<c_{k\sigma }^{}f_{0\sigma }c_{p\sigma }|f_{0\sigma }^{}n_{0\sigma }>>)`$
Now let us see how to proceed further to get a suitable functional structure of the relevant solution. The intrinsic nature of the system of the equations of motion (9.7) - (291) suggests to consider the following approximation:

(292) $`(\omega ϵ_k)<<c_{k\sigma }n_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>V<<f_{0\sigma }n_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>`$
$`(\omega ϵ_kE_{0\sigma }+E_{0\sigma })<<c_{k\sigma }f_{0\sigma }^{}f_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>><f_{0\sigma }^{}c_{k\sigma }n_{0\sigma }>`$
(293) $`V(<<f_{0\sigma }n_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>><<c_{k\sigma }c_{k\sigma }^{}f_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>)`$
$`(\omega +ϵ_kE_{0\sigma }E_{0\sigma }U)<<c_{k\sigma }^{}f_{0\sigma }f_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>><c_{k\sigma }^{}f_{0\sigma }f_{0\sigma }^{}f_{0\sigma }>+`$
(294) $`V(<<f_{0\sigma }n_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>+<<c_{k\sigma }^{}f_{0\sigma }c_{k\sigma }|f_{0\sigma }^{}n_{0\sigma }>>)`$
It is transparent that the construction of approximations (292) - (9.7) is related with the small-V expansion and is not unique, but very natural. As a result, we find the explicit expression for GF in (287)

(295)
$$<<f_{0\sigma }n_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>\frac{<n_{0\sigma }>F_\sigma ^1(\omega )}{\omega E_{0\sigma }US_1(\omega )}$$

Here the following notation was used

(296) $`S_1(\omega )=S(\omega )`$
$`+{\displaystyle \underset{k}{}}|V|^2({\displaystyle \frac{1}{\omega ϵ_kE_{0\sigma }+E_{0\sigma }}}+{\displaystyle \frac{1}{\omega +ϵ_kE_{0\sigma }E_{0\sigma }U}})`$
(297) $`F_\sigma ^1={\displaystyle \underset{k}{}}(VF_2+V^2F_3)`$
(298) $`F_2={\displaystyle \frac{<c_{k\sigma }^{}f_{0\sigma }f_{0\sigma }^{}f_{0\sigma }>}{\omega +ϵ_kE_{0\sigma }E_{0\sigma }U}}+{\displaystyle \frac{<f_{0\sigma }^{}c_{k\sigma }n_{0\sigma }>}{\omega ϵ_kE_{0\sigma }+E_{0\sigma }}}`$
(299) $`F_3=`$
$`{\displaystyle \frac{<<c_{k\sigma }c_{k\sigma }^{}f_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>}{\omega ϵ_kE_{0\sigma }+E_{0\sigma }}}+{\displaystyle \frac{<<c_{k\sigma }^{}f_{0\sigma }c_{k\sigma }|f_{0\sigma }^{}n_{0\sigma }>>}{\omega +ϵ_kE_{0\sigma }E_{0\sigma }U}}`$
Now one can substitute the GF in (287) by the expression (295). This gives a new approximate dynamic solution of SIAM where the complex expansion both in $`U`$ and $`V`$ was incorporated. The important observation is that this new solution satisfies both the limits (9.5). For example, if we wish to get a lowest order approximation up to $`U^2`$ and $`V^2`$, it is very easy to notice that for $`V=0`$:

$`<<f_{0\sigma }c_{k\sigma }^{}c_{k\sigma }|f_{0\sigma }^{}n_{0\sigma }>>{\displaystyle \frac{<c_{k\sigma }^{}c_{k\sigma }><n_{0\sigma }>}{\omega E_{0\sigma }U}}`$
(300) $`<<c_{k\sigma }c_{k\sigma }^{}f_{0\sigma }|f_{0\sigma }^{}n_{0\sigma }>>{\displaystyle \frac{<c_{k\sigma }c_{k\sigma }^{}><n_{0\sigma }>}{\omega E_{0\sigma }U}}`$
This results in the possibility to find explicitly all necessary quantities and, thus, to solve the problem in a self-consistent way.
In summary, we presented here a consistent many-body approach to analytic dynamic solution of SIAM at finite temperatures and for a broad interval of the values of the model parameters. We used the exact result (285) to connect the single-particle GF with higher-order GF to obtain a complex combined expansion in terms of $`U`$ and $`V`$ for the propagator. To summarize, we reformulated the problem of searches for an appropriate many-body dynamic solution for SIAM in a way that provides us with an effective and workable scheme for constructing of advanced analytic approximate solutions for the single-particle GFs on the level of the higher-order GFs in a rather systematic self-consistent way. This procedure has the advantage that it systematically uses the principle of interpolation solution within the equation-of-motion approach for GFs. The leading principle, which we used here was to look more carefully for the intrinsic functional structure of the required relevant solution and then to formulate approximations for the higher-order GFs in accordance with this structure.
The main results of our IGF study are the exact Dyson equations for the full matrix GFs and a new derivation of the GMF GFs . The approximate explicit calculations of inelastic self-energy corrections are quite straightforward but tedious and too extended for their description. Here we want to emphasize an essentially new point of view on the derivation of the Generalized Mean Fields for SIAM when we are interested in the interpolation finite temperature solution for the single-particle propagator. Our final solutions have the correct functional structure and differ essentially from previous solutions.
Of course, there are important criteria to be met (mainly numerically) , such as the question left open, whether the present approximation satisfies the Friedel sum rule (this question left open in  and  ). A quantitative numerical comparison of self-consistent results (e.g. the width and shape of the Kondo resonance in the near-integer regime of the SIAM) would be crucial too. In the present consideration, we concentrated on the problem of correct functional structure of the single-particle GF itself. In addition to SIAM, it will be instructive to consider sketchy the PAM and TIAM too for completness.

### 9.8 Quasi-Particle Dynamics of PAM

The main drawback of the H-F type solution of PAM (59) is that it ignores the correlations of the ””up” and ”down” electrons. In this Section, we take into account the latter correlations in a self-consistent way using the IGF method. We consider the relevant matrix GF of the form ( cf. (220) )

(301)
$$\widehat{G}(\omega )=\left(\begin{array}{cc}<<c_{k\sigma }|c_{k\sigma }^{}>>& <<c_{k\sigma }|f_{k\sigma }^{}>>\\ <<f_{k\sigma }|c_{k\sigma }^{}>>& <<f_{k\sigma }|f_{k\sigma }^{}>>\end{array}\right)$$

The equation of motion for GF (301) reads

$`\left(\begin{array}{cc}(\omega ϵ_k)& V_k\\ V_k& (\omega E_k)\end{array}\right)\left(\begin{array}{cc}<<c_{k\sigma }|c_{k\sigma }^{}>>& <<c_{k\sigma }|f_{k\sigma }^{}>>\\ <<f_{k\sigma }|c_{k\sigma }^{}>>& <<f_{k\sigma }|f_{k\sigma }^{}>>\end{array}\right)=`$
(302) $`\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)+UN^1{\displaystyle \underset{pq}{}}\left(\begin{array}{cc}0& 0\\ <<A|c_{k\sigma }^{}>>& <<A|f_{k\sigma }^{}>>\end{array}\right)`$
where $`A=f_{k+p\sigma }f_{p+q\sigma }^{}f_{q\sigma }`$. According to eq.(30), the definition of the irreducible parts in the equation of motion (9.8) are defined as follows

$`{}_{}{}^{(ir)}<<f_{k+p\sigma }f_{p+q\sigma }^{}f_{q\sigma }|c_{k\sigma }^{}>>=<<f_{k+p\sigma }f_{p+q\sigma }^{}f_{q\sigma }|c_{k\sigma }^{}>>`$
$`\delta _{p,0}<n_{q\sigma }><<f_{k\sigma }|c_{k\sigma }^{}>>`$
$`{}_{}{}^{(ir)}<<f_{k+p\sigma }f_{p+q\sigma }^{}f_{q\sigma }|f_{k\sigma }^{}>>=<<f_{k+p\sigma }f_{p+q\sigma }^{}f_{q\sigma }|f_{k\sigma }^{}>>`$
$`\delta _{p,0}<n_{q\sigma }><<f_{k\sigma }|f_{k\sigma }^{}>>`$
After substituting these definitions into equation (9.8), we obtain

$`\left(\begin{array}{cc}(\omega ϵ_k)& V_k\\ V_k& (\omega E_\sigma (k))\end{array}\right)\left(\begin{array}{cc}<<c_{k\sigma }|c_{k\sigma }^{}>>& <<c_{k\sigma }|f_{k\sigma }^{}>>\\ <<f_{k\sigma }|c_{k\sigma }^{}>>& <<f_{k\sigma }|f_{k\sigma }^{}>>\end{array}\right)=`$
(303) $`\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)+UN^1{\displaystyle \underset{pq}{}}\left(\begin{array}{cc}0& 0\\ {}_{}{}^{(ir)}<<A|c_{k\sigma }^{}>>& {}_{}{}^{(ir)}<<A|f_{k\sigma }^{}>>\end{array}\right)`$
The following notation was used

$$E_\sigma (k)=E_kUn_\sigma ;n_\sigma =<f_{k\sigma }^{}f_{k\sigma }>$$

The definition of the generalized mean field GF ( which, for the weak Coulomb correlation $`U`$, coincides with the Hartree-Fock mean field ) is evident. All inelastic renormalization terms are now related to the last term in the equation of motion (9.8). All elastic scattering ( or mean field) renormalization terms are included into the following mean-field GF

$$\left(\begin{array}{cc}(\omega ϵ_k)& V_k\\ V_k& (\omega E_\sigma (k))\end{array}\right)\left(\begin{array}{cc}<<c_{k\sigma }|c_{k\sigma }^{}>>^0& <<c_{k\sigma }|f_{k\sigma }^{}>>^0\\ <<f_{k\sigma }|c_{k\sigma }^{}>>^0& <<f_{k\sigma }|f_{k\sigma }^{}>>^0\end{array}\right)=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)$$

It is easy to find that (cf. (225) and (226))

(304) $`<<f_{k\sigma }|f_{k\sigma }^{}>>^0=\left(\omega E_\sigma (k){\displaystyle \frac{|V_k|^2}{\omega ϵ_k}}\right)^1`$
(305) $`<<c_{k\sigma }|c_{k\sigma }^{}>>^0=\left(\omega ϵ_k{\displaystyle \frac{|V_k|^2}{\omega E_\sigma (k)}}\right)^1`$
At this point, it is worthwhile to emphasize a significant difference between both the models, PAM and SIAM. The corresponding SIAM equation for generalized mean field GF (224) reads

$`{\displaystyle \underset{p}{}}\left(\begin{array}{cc}(\omega ϵ_p)\delta _{pk}& V_p\delta _{pk}\\ V_p& \frac{1}{N}(\omega E_{0\sigma }Un_\sigma )\end{array}\right)\left(\begin{array}{cc}<<c_{k\sigma }|c_{k\sigma }^{}>>^0& <<c_{k\sigma }|f_{0\sigma }^{}>>^0\\ <<f_{0\sigma }|c_{k\sigma }^{}>>^0& <<f_{0\sigma }|f_{0\sigma }^{}>>^0\end{array}\right)=`$
(306) $`\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)`$
This matrix notation for SIAM shows a fundamental distinction between SIAM and PAM. For SIAM, we have a different number of states for a strongly localized level and the conduction electron subsystem: the conduction band contains $`2N`$ states, whereas the localized (s-type) level contains only two. The comparison of (9.8) and (9.8) shows clearly that this difficulty does not exist for PAM : the number of states both in the localized and itinerant subsystems are the same, i.e. $`2N`$.
This important difference between SIAM and PAM appears also when we calculate inelastic scattering or self-energy corrections. By analogy with the Hubbard model, the equation of motion (9.8) for PAM can be transformed exactly to the scattering equation of the form ( 36). Then, we are able to write down explicitly the Dyson equation (37) and the exact expression for the self-energy M in the matrix form:

(307)
$$\widehat{M}_{k\sigma }(\omega )=\left(\begin{array}{cc}0& 0\\ 0& M_{22}\end{array}\right)$$

Here the matrix element $`M_{22}`$ is of the form

(308) $`M_{22}=M_{k\sigma }(\omega )=`$
$`{\displaystyle \frac{U^2}{N^2}}{\displaystyle \underset{pqrs}{}}(^{(ir)}<<f_{k+p\sigma }f_{p+q\sigma }^{}f_{q\sigma }|f_{r\sigma }^{}f_{r+s\sigma }f_{k+s\sigma }^{}>>^{(ir)})^{(p)}`$
To calculate the self-energy operator (308) in a self-consistent way, we proceed by analogy with the Hubbard model in Section 8.1. Then we find both expressions for the self-energy operator in form (145) and (148).

### 9.9 Quasi-Particle Dynamics of TIAM

Let us see now how to rewrite the results of the preceeding Sections for the case of TIAM Hamiltonian (60). We again consider the relevant matrix GF of the form (cf.(220) )

(309)
$$\widehat{G}(\omega )=\left(\begin{array}{ccc}G_{11}& G_{12}& G_{13}\\ G_{21}& G_{22}& G_{23}\\ G_{31}& G_{32}& G_{33}\end{array}\right)=\left(\begin{array}{ccc}<<c_{k\sigma }|c_{k\sigma }^{}>>& <<c_{k\sigma }|f_{1\sigma }^{}>>& <<c_{k\sigma }|f_{2\sigma }^{}>>\\ <<f_{1\sigma }|c_{k\sigma }^{}>>& <<f_{1\sigma }|f_{1\sigma }^{}>>& <<f_{1\sigma }|f_{2\sigma }^{}>>\\ <<f_{2\sigma }|c_{k\sigma }^{}>>& <<f_{2\sigma }|f_{1\sigma }^{}>>& <<f_{2\sigma }|f_{2\sigma }^{}>>\end{array}\right)$$

The equation of motion for GF (309) reads

$`{\displaystyle \underset{p}{}}\left(\begin{array}{ccc}(\omega ϵ_p)\delta _{pk}& V_{1p}\delta _{pk}& V_{1p}\delta _{pk}\\ V_{1p}& \frac{1}{N}(\omega E_{0\sigma })& V_{12}\\ V_{2p}& V_{21}& \frac{1}{N}(\omega E_{0\sigma })\end{array}\right)\left(\begin{array}{ccc}G_{11}& G_{12}& G_{13}\\ G_{21}& G_{22}& G_{23}\\ G_{31}& G_{32}& G_{33}\end{array}\right)=`$
(310) $`\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)+`$
$`U\left(\begin{array}{ccc}0& 0& 0\\ <<A_1|c_{k\sigma }^{}>>& <<A_1|f_{1\sigma }^{}>>& <<A_1|f_{2\sigma }^{}>>\\ <<A_2|c_{k\sigma }^{}>>& <<A_2|f_{1\sigma }^{}>>& <<A_2|f_{2\sigma }^{}>>\end{array}\right)`$
The notation is as follows

$$A_1=f_{1\sigma }f_{1\sigma }^{}f_{1\sigma };A_2=f_{2\sigma }f_{2\sigma }^{}f_{2\sigma }$$

In a compact notation, the equation (9.9) has the form (cf. (251))

(311)
$$\underset{p}{}F(p,k)G_{pk}(\omega )=\widehat{I}+UD_p(\omega )$$

We thus have the equatin of motion (311) which is a complete analogue of the corresponding equations for the SIAM and PAM. After introducing the irreducible parts by analogy with the equation (221)

$`{}_{}{}^{(ir)}<<f_{1\sigma }f_{1\sigma }^{}f_{1\sigma }|B>>_\omega =<<f_{1\sigma }f_{1\sigma }^{}f_{1\sigma }|B>>_\omega `$
$`<n_{1\sigma }><<f_{1\sigma }|B>>_\omega `$
$`{}_{}{}^{(ir)}<<f_{2\sigma }f_{2\sigma }^{}f_{2\sigma }|B>>_\omega =<<f_{2\sigma }f_{2\sigma }^{}f_{2\sigma }|B>>_\omega `$
$`<n_{2\sigma }><<f_{2\sigma }|B>>_\omega `$
,
and performing the second-time differentiation of the higher-order GF, and introducing the relevant irreducible parts, the equation of motion (311) is rewritten in the form of Dyson equation (37). The definition of the generalized mean field GF is as follows

$`{\displaystyle \underset{p}{}}\left(\begin{array}{ccc}(\omega ϵ_p)\delta _{pk}& V_{1p}\delta _{pk}& V_{1p}\delta _{pk}\\ V_{1p}& \frac{1}{N}(\omega E_{0\sigma }Un_\sigma )& V_{12}\\ V_{2p}& V_{21}& \frac{1}{N}(\omega E_{0\sigma }Un_\sigma )\end{array}\right)`$
(312) $`\left(\begin{array}{ccc}G_{11}^0& G_{12}^0& G_{13}^0\\ G_{21}& G_{22}^0& G_{23}^0\\ G_{31}^0& G_{32}^0& G_{33}^0\end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)`$
The matrix GF (9.9) describes the mean-field solution of the TIAM Hamiltonian. The explicit solutions for diagonal elements of $`G^0`$ are ( cf. (225))

(313) $`<<c_{k\sigma }|c_{k\sigma }^{}>>_\omega ^0=\left(\omega ϵ_k{\displaystyle \frac{|V_{1k}|^2}{\omega (E_{0\sigma }Un_\sigma )}}\mathrm{\Delta }_{11}(k,\omega )\right)^1`$
(314) $`<<f_{1\sigma }|f_{1\sigma }^{}>>_\omega ^0=(\omega (E_{0\sigma }Un_\sigma )S(\omega ))\mathrm{\Delta }_{22}(k,\omega ))^1`$
(315) $`<<f_{2\sigma }|f_{2\sigma }^{}>>_\omega ^0=(\omega (E_{0\sigma }Un_\sigma )S(\omega ))\mathrm{\Delta }_{33}(k,\omega ))^1`$
Here we introduced the notation

$`\mathrm{\Delta }_{11}(k,\omega )=\left(V_{2k}+{\displaystyle \frac{V_{1k}V_{12}}{\omega (E_{0\sigma }Un_\sigma )}}\right)\left(V_{2k}+{\displaystyle \frac{V_{1k}V_{21}}{\omega (E_{0\sigma }Un_\sigma )}}\right)`$
$`[\omega (E_{0\sigma }Un_\sigma ){\displaystyle \frac{V_{21}V_{12}}{\omega (E_{0\sigma }Un_\sigma )}}]^1`$
$`\mathrm{\Delta }_{22}(k,\omega )=(\lambda _{21}(\omega )+V_{12})(\lambda _{21}(\omega )+V_{21})[\omega (E_{0\sigma }Un_\sigma ){\displaystyle \frac{_p|V_{2p}|^2}{\omega ϵ_p}}]^1`$
$`\mathrm{\Delta }_{33}(k,\omega )=(\lambda _{12}(\omega )+V_{21})(\lambda _{12}(\omega )+V_{12})[\omega (E_{0\sigma }Un_\sigma ){\displaystyle \frac{_p|V_{1p}|^2}{\omega ϵ_p}}]^1`$
(316) $`\lambda _{12}=\lambda _{21}={\displaystyle \underset{p}{}}{\displaystyle \frac{V_{1p}V_{2p}}{\omega ϵ_p}}`$
The formal solution of the Dyson equation for TIAM contains the self-energy matrix

(317)
$$\widehat{M}=\left(\begin{array}{ccc}0& 0& 0\\ 0& M_{22}& M_{23}\\ 0& M_{32}& M33\end{array}\right)$$

where

(318) $`M_{22}=U^2(^{(ir)}<<f_{1\sigma }n_{1\sigma }|f_{1\sigma }^{}n_{1\sigma }>>^{(ir)})^p`$
$`M_{32}=U^2(^{(ir)}<<f_{2\sigma }n_{2\sigma }|f_{1\sigma }^{}n_{1\sigma }>>^{(ir)})^p`$
$`M_{23}=U^2(^{(ir)}<<f_{1\sigma }n_{1\sigma }|f_{2\sigma }^{}n_{2\sigma }>>^{(ir)})^p`$
$`M_{33}=U^2(^{(ir)}<<f_{2\sigma }n_{2\sigma }|f_{2\sigma }^{}n_{2\sigma }>>^{(ir)})^p`$
To calculate the matrix elements (318), the same procedure can be used as it was done previously for the SIAM (9.2). As a result, we find the following explicit expressions for the self-energy matrix elements ( cf.(9.2)

$`M_{22}^{}(\omega )=U^2{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\omega _1𝑑\omega _2{\displaystyle \frac{1+N(\omega _1)n(\omega _2)}{\omega \omega _1\omega _2}}`$
$`({\displaystyle \frac{1}{\pi }}Im<<S_1^{}|S_1^+>>_{\omega _1})`$
(319) $`({\displaystyle \frac{1}{\pi }}Im<<f_1|f_1^{}>>_{\omega _2})`$
$`M_{22}^{}(\omega )=U^2{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\omega _1𝑑\omega _2{\displaystyle \frac{1+N(\omega _1)n(\omega _2)}{\omega \omega _1\omega _2}}`$
$`({\displaystyle \frac{1}{\pi }}Im<<S_1^+|S_1^{}>>_{\omega _1})`$
(320) $`({\displaystyle \frac{1}{\pi }}Im<<f_1|f_1^{}>>_{\omega _2})`$
$`M_{23}^{}(\omega )=U^2{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\omega _1𝑑\omega _2{\displaystyle \frac{1+N(\omega _1)n(\omega _2)}{\omega \omega _1\omega _2}}`$
$`({\displaystyle \frac{1}{\pi }}Im<<S_1^{}|S_2^+>>_{\omega _1})`$
(321) $`({\displaystyle \frac{1}{\pi }}Im<<f_1|f_2^{}>>_{\omega _2})`$
$`M_{23}^{}(\omega )=U^2{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\omega _1𝑑\omega _2{\displaystyle \frac{1+N(\omega _1)n(\omega _2)}{\omega \omega _1\omega _2}}`$
$`({\displaystyle \frac{1}{\pi }}Im<<S_2^+|S_1^{}>>_{\omega _1})`$
(322) $`({\displaystyle \frac{1}{\pi }}Im<<f_1|f_2^{}>>_{\omega _2})`$
where the following notation was used:

$$S_i^+=f_i^{}f_i;S_i^{}=f_i^{}f_i;i=1,2$$

For $`M_{33}`$ we obtain the same expressions as for $`M_{22}`$ with the substitution of index 1 by 2. For $`M_{32}^{}`$ we must do the same. It is possible to say that the diagonal elements $`M_{22}`$ and $`M_{33}`$ describe single-site inelastic scattering processes; off-diagonal elements $`M_{23}`$ and $`M_{32}`$ describe intersite inelastic scattering processes. They are responsible for the specific features of the dynamic behaviour of TIAM ( as well as the off-diagonal matrix elements of the GF $`G^0`$) and, more generally, the cluster impurity Anderson model (CIAM). The nonlocal contributions to the total spin susceptibility of two well formed impurity magnetic moments at a distance $`R`$ can be estimated as

(323)
$$\chi _{pair}<<S_1^{}|S_2^+>>2\chi 12\pi E_F(\frac{\chi }{g\mu _B})^2\frac{\mathrm{cos}(2k_FR)}{(k_FR)^3}$$

In the region of interplay of the RKKY and Kondo behaviour, the key point is then to connect the partial Kondo screening effects with the low temperature behaviour of the total spin susceptibility. As it is known, it is quite difficult to describe such a threshold behaviour analytically. However, progress is expected due to a better understanding of the quasi-particle many-body dynamics both from analytical and numerical investigations.

## 10 Conclusions

In the present paper, we have formulated the theory of the correlation effects for many-particle interacting systems using the ideas of quantum field theory for interacting electron and spin systems on a lattice. The workable and self-consistent IGF approach to the decoupling problem for the equation-of-motion method for double-time temperature Green functions has been presented. The main achievement of this formulation was the derivation of the Dyson equation for double-time retarded Green functions instead of causal ones. That formulation permits to unify convenient analytical properties of retarded and advanced GF and the formal solution of the Dyson equation (38), that, in spite of the required approximations for the self-energy, provides the correct functional structure of single-particle GF. The main advantage of the mathematical formalism is brought out by showing how elastic scattering corrections (generalized mean fields) and inelastic scattering effects (damping and finite lifetimes) could be self-consistently incorporated in a general and compact manner. In this paper, we have thoroughly considered the idealized Anderson and Hubbard models which are the simplest (in the sense of formulation, but not solution) and most popular models of correlated lattice fermions. We have presented here the novel method of calculation of quasi-particle spectra for these and basic spin lattice models, as the most representative examples. Using the IGF method, we were able to obtain a closed self-consistent set of equations determining the electron GF and self-energy. For the Hubbard and Anderson models, these equations give a general microscopic description of correlation effects both for the weak and strong Coulomb correlation, and,thus, determine the interpolation solutions of the models. Moreover, this approach gives the workable scheme for the definition of relevant generalized mean fields written in terms of appropriate correlators.
We hope that these considerations have been done with sufficient details to bring out their scope and power, since we believe that this technique will have application to a variety of many-body systems with complicated spectra and strong interaction. The application of the IGF method to the investigation of nonlocal correlations and quasi-particle interactions in Anderson models  has a particular interest for studying of the intersite correlation effects in the concentrated Kondo system . A comparative study of real many-body dynamics of single-impurity, two-impurity, and periodic Anderson model, especially for strong but finite Coulomb correlation, when perturbation expansion in $`U`$ does not work, is important to characterize the true quasi-particle excitations and the role of magnetic correlations. It was shown that the physics of two-impurity Anderson model can be understood in terms of competition between itinerant motion of carriers and magnetic correlations of the RKKY nature. This issue is still very controversial and the additional efforts must be applied in this field.
The application of the IGF method to the theory of magnetic semiconductors was done in , . As a remarkable result of our approach, let us mention the generalization of the Shastry-Mattis theory for the magnetic polaron to the finite temperatures  . The quasi-particle many-body dynamics of ferromagnetic and aniferromagnetic semiconductors, was studied too. These studies clarified greatly the true nature of carriers in magnetic semiconductors. The application of the IGF method to generalized spin-fermion models that was made in papers , allows one to consider carefully the true nature of carriers in oxides and rare-earth metals. These applications illustrate some of subtle details of the IGF approach and exhibit their physical significance in a representative form.
As it is seen, this treatment has advantages in comparison with the standard methods of decoupling of higher-order GFs within the equation-of-motion approach, namely, the following:

* At the mean-field level, the GF, one obtains, is richer than that following from the standard procedures. The generalized mean fields represent all elastic scattering renormalizations in a compact form.
* The approximations ( the decoupling ) are introduced at a later stage with respect to other methods, i.e. only into the rigorously obtained self-energy.
* Many standard results of the many-particle system theory are reproduced mathematically incomparable more simply.
* The physical picture of elastic and inelastic scattering processes in the interacting many-particle systems is clearly seen at every stage of calculations, which is not the case with the standard methods of decoupling.
* The main advantage of the whole method is the possibility of a self-consistent description of quasi-particle spectra and their damping in a unified and coherent fashion.
* This new picture of interacting many-particle systems on a lattice is far richer and gives more possibilities for the analysis of phenomena which can actually take place. In this sense the approach we suggest produces more advanced physical picture of the quasi-particle many-body dynamics.

Despite the novelty of the IGF techniques introduced above and some (not really big) complexity of the details in its demonstrations, the major conclusions of the present paper can be made intelligible to any reader. The most important conclusion to be drawn from the present consideration is that the GMF for the case of strong Coulomb interaction has quite a nontrivial structure and cannot be reduced to the mean-density functional. This last statement resembles very much the situation with strongly non-equilibrium systems, where only the single-particle distribution function is insufficient to describe the essence of the strongly non-equilibrium state Therefore a more complicated correlation functions are to be taken into account, in accordance with general ideas of Bogoliubov and Mori-Zwanzig. The IGF method is intimately related to the projection method in the sense, that it expresses the idea of “reduced description” of a system in the most general form. This line of consideration is very promising for developing the complete and self-contained theory of strongly interacting many-body systems on a lattice. Our main results reveal the fundamental importance of the adequate definition of Generalized Mean Fields at finite temperatures, that results in a deeper insight into the nature of quasi-particle states of the correlated lattice fermions and spins. We believe that our approach offers a new way for systematic constructions of the approximate dynamic solutions of the Hubbard, SIAM, TIAM, PAM, spin-fermion, and other models of the strongly correlated electron systems on a lattice. The work in this direction is in progress.

## 11 Acknowledgments

I would like to dedicate this article to the memory of the late Professors S.V.Tyablikov, N.N.Bogoliubov, and D.N.Zubarev. Their illuminating and deep remarks, advice, and suggestions were indispensable stimulus for my studies. I express my gratitude to them.

## Appendix A  Appendix . The Gram-Schmidt Orthogonalization Procedure

In this appendix we briefly recall the Gram-Schmidt Orthogonalization Procedure. The Gram-Schmidt orthogonalization procedure is an inductive technique to generate a mutually orthogonal set from any linearly independent set of vectors.
Suppose we have an arbitrary n-dimensional Euclidean space, which means that scalar multiplication has been introduced in some fashion into an n-dimensional linear space. The vectors $`f`$ and $`g`$ are orthogonal if their scalar product is zero

$`(A.1)`$
$$(f,g)=0$$

We now describe the orthogonalization process, which is a means of passing from any linearly independent system of $`k`$ vectors $`f_1,f_2,\mathrm{}f_k`$ to an orthogonal system, also consisting of $`k`$ nonzero vectors. We denote these vectors by $`g_1,g_2,\mathrm{}g_k`$.
Let us put $`g_1=f_1`$, which is to say that the first vector of our system will enter into the orthogonal system we are building. After that, put

$`(A.2)`$
$$g_2=f_2+\alpha g_1$$

Since $`g_1=f_1`$ and the vectors $`f_1`$ and $`f_2`$ are linearly independent, it follows that the vector $`g_2`$ is different from zero for any scalar $`\alpha `$. We choose this scalar from the constraint

$`(A.3)`$
$$0=(g_1,g_2)=\alpha (g_1,g_1)+(g_1,f_2)$$

whence

$`(A.4)`$
$$\alpha =\frac{(g_1,f_2)}{(g_1,g_1)}$$

In other words, we get $`g_2`$ by subtracting from $`f_2`$ the projection of $`f_2`$ onto $`g_1`$. Proceeding inductively, we find

$`(A.5)`$
$$g_n=f_n\underset{j=1}{\overset{n1}{}}\frac{(g_j,f_n)}{(g_j,g_j)}g_j$$

We are left with mutually orthogonal vectors which have the same span as the original set.
Let us consider an important example of a basis $`f_1,f_2,f_3,f_4`$ in a 4-dimensional space and then construct the orthonormal basis of the same space. Next, in the equality $`g_3=f_3+\beta _1g_1+\beta _2g_2`$, chose $`\beta _1`$ and $`\beta _2`$ such that the conditions $`g_3g_1,g_3g_2`$ are fulfilled.
From the equalities

$`(A.6)`$
$$(g_1,g_3)=(g_1,f_3)+\beta _1(g_1,g_1)+\beta _2(g_1,g_2)$$

$`(A.7)`$
$$(g_2,g_3)=(g_2,f_3)+\beta _1(g_1,g_2)+\beta _2(g_2,g_2)$$

we obtain

$`(A.8)`$
$$\beta _1=\frac{(g_1,f_3)}{(g_1,g_1)};\beta _2=\frac{(g_2,f_3)}{(g_2,g_2)}$$

Finally, from the equality $`g_4=f_4+\gamma _1g_1+\gamma _2g_2+\gamma _3g_3`$ we find

$`(A.9)`$
$$\gamma _1=\frac{(g_1,f_4)}{(g_1,g_1)};\gamma _2=\frac{(g_2,f_4)}{(g_2,g_2)};\gamma _3=\frac{(g_3,f_4)}{(g_3,g_3)}$$

Thus, we see that with the choice of $`\alpha ,\beta _1,\beta _2,\gamma _1,\gamma _2,\gamma _3`$ made, the vectors $`g_1,g_2,g_3,g_4`$ are pairwise orthogonal.

## Appendix B Appendix . Moments and Green Functions

It is known that the method of moments  of spectral density is considered sometimes as an alternative approach for describing the many-body quasi-particle dynamics of interacting many-particle systems. The moments technique appears naturally when studying the particle dynamics in many-particle systems in the context of time-dependent correlation functions ( magnetic resonance, liquids, etc.). Qualitatively, a correlation function describes how long a given property of a system persists until it is averaged out by the microscopic motion of particles in the macroscopic system. The time dependence of a particle correlation function sometimes is approximated (at small times) via a power series expansion about the initial time 0.

$`(B.1)`$
$$<A(0)A(t)>=\underset{n=0}{\overset{\mathrm{}}{}}\frac{t^n}{n!}\frac{d^n}{dt^n}<A(0)A(t)>|_{t=0}=$$

$$\underset{n=0}{\overset{\mathrm{}}{}}\frac{(it)^n}{n!}<A(0)[H,[H\mathrm{}[H,A(0)]\mathrm{}]]]>$$

The spectral theorem (26), (27) connects $`A(\omega )`$ and the correlation functions. From the above expression we obtain the moments $`M_n`$ of the spectral density function

$`(B.2)`$
$$M_n=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}𝑑\omega \omega ^nA(\omega )=(1)^n<[[H,[H\mathrm{}[H,A]\mathrm{}]],B]_\eta >$$

So, by definition, the moments are time-independent correlation functions of a combination of the operators. In principle, it is possible to calculate them in a regular way; however, in practice, it is possible to do this only for a first few moments. If the moments $`M_n`$ of a given spectral density form a positive sequence, then GF of appropriate operators is a limit of the sequence

$`(B.3)`$
$$G(E)=\underset{n\mathrm{}}{lim}G_n(E,\gamma )$$

Here the parameter $`\mathrm{}<\gamma <+\mathrm{}`$ and is real. The approximation procedure for GF consists in replacing the $`G(E)`$ by $`G_n(E,\gamma )`$, that depends also on the appropriate choice of the parameter $`\gamma `$. The $`G_n(E,\gamma )`$ have the properties

$`(B.4)`$
$$G_n(E,\mathrm{})=G_{n1}(E,0)$$

and are represented by the fraction

$`(B.5)`$
$$G_n(E,\gamma )=M_0\frac{Q_{n+1}(E)\gamma Q_n(E)}{P_{n+1}(E)\gamma P_n(E)}$$

The polynomials $`P_n`$ are given by the determinant

$`(B.6)`$
$$P_{n1}(E)=\frac{\sqrt{M_0}}{\sqrt{D_{n1}D_n}}\left|\begin{array}{cccc}M_0& M_1& \mathrm{}& M_n\\ M_1& M_2& \mathrm{}& M_{n+1}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ M_{n1}& M_n& \mathrm{}& M_{2n1}\\ 1& E& \mathrm{}& E^n\end{array}\right|$$

$$P_0=1$$

where

$`(B.7)`$
$$D_{n1}=\left|\begin{array}{cccc}M_0& M_1& \mathrm{}& M_n\\ M_1& M_2& \mathrm{}& M_{n+1}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ M_n& M_{n+1}& \mathrm{}& M_{2n}\end{array}\right|$$

$$D_0=D_1=M_0$$

The polynomial $`Q_n(E)`$ (which is of (n-1)-th order in E) is related to the polynomial $`P_n(E)`$ (which is of n-th order in E) via the following relation

$`(B.8)`$
$$Q_n(E)=\frac{1}{2\pi M_0}_{\mathrm{}}^{\mathrm{}}\frac{P_n(E)P_n(\omega )}{E\omega }A(\omega )𝑑\omega $$

It is possible to find a few lowest-order terms

$`(B.9)`$
$$P_0(E)=1;P_1(E)=\frac{E\frac{M_1}{M_0}}{M_2M_0^1}$$

$`(B.10)`$
$$Q_0(E)=0;Q_1=\frac{1}{M_2M_0^1}$$

The expression (B.5) can be represented in the following form

$`(B.11)`$
$$G_n(E,\gamma )=M_0\underset{i=1}{\overset{n+1}{}}\frac{m_i(\gamma )}{EE_i(\gamma )}$$

Here the numbers $`E_i(\gamma )`$ are roots of the equation

$`(B.12)`$
$$P_{n+1}(E)\gamma P_n(E)=0$$

These relations lead to the possibility of practical applications of the moment expansion method. If we know the first $`(2n+2)`$ moments, then the equation (B.12) determines $`(n+1)`$ different roots $`E_i(\gamma )`$. Thus, the spectral density function can be represented by

$`(B.13)`$
$$A(\omega )=2\pi M_0\underset{i=1}{\overset{n+1}{}}m_i\delta (\omega E_i)$$

For example, if we know the moments $`M_0,M_1,M_2`$ then we find, from the equation (B.11), the roots of (B.12)

$`(B.14)`$
$$E_1(\gamma )=M_1M_0^1+\gamma (M_2M_0^1)$$

In this approximation, the GF and corresponding spectral density are represented as

$`(B.15)`$
$$G_0(\gamma )=\frac{M_0}{EE_1(\gamma )};A(\omega )=2\pi M_0\delta (\omega E_1)$$

It is clear that the Tyablikov decoupling approximation (43) corresponds to this approximation within the moment method. An improved decoupling scheme, that conserves the first several frequency moments of the spectral weight function for the Heisenberg and Hubbard models was developed in paper (cf. , ).
It was shown in ref.  that the IGF method permits one to calculate the spectral density for the spin-fermion model in the approximation that preserves the first four moments. This is valid also for the approximation used for the strongly correlated Hubbard model in Section 7.2.
It must be clear from the above consideration that the structure of the obtained solution for single-particle GF depends strongly on the stage at which irreducible parts were introduced . To clarify this, let us consider equation (29) again. Instead of (30), we introduce now the IGFs in the following way

$$\omega G(\omega )=M_0+<<[A,H]_{}A^{}>>_\omega $$

$`(B.16)`$
$$\omega <<[A,H]|A^{}>>=M_1+(^{(ir)}<<[[A,H]H]A^{}>>_\omega )+$$

$$\alpha _1<<A|A^{}>>_\omega +\alpha _2<<[A,H]|A^{}>>_\omega $$

The unknown constants $`\alpha _1`$ and $`\alpha _2`$ are connected by the orthogonality condition

$`(B.17)`$
$$<[[[A,H]H]^{(ir)},A^{}]>=0$$

For illustration, we consider the simplest possibility and write down the following equation

$`(B.18)`$
$$\omega (^{(ir)}<<[[A,H]H]|A^{}>>)=(^{(ir)}<<[[A,H]H]|[H,A^{}]>>)$$

Then by introducing the irreducible parts for the $`\mathrm{𝑟𝑖𝑔ℎ𝑡}`$ operators, we obtain

$`(B.19)`$
$$(^{(ir)}<<[[A,H]H]|A^{}>>)(\omega \alpha _1^{})=(^{(ir)}<<[[A,H]H]|[H,A^{}]>>^{(ir)})$$

It is clear enough that, as a result, we arrive at the following set of equations

$$\omega <<A|A^{}>>_\omega <<[A,H]_{}A^{}>>_\omega =M_0$$

$`(B.20)`$
$$\alpha _1<<A|A^{}>>_\omega +(\omega \alpha _2)<<[A,H]|A^{}>>_\omega =M_1\mathrm{\Phi }$$

where

$`(B.21)`$
$$\mathrm{\Phi }=(^{(ir)}<<[[A,H]H][A,H]^{}>>_\omega ^{(ir)})$$

The solutions of the equations (B.20) are given by

$`(B.22)`$
$$<<A|A^{}>>_\omega =\frac{M_0(\omega \alpha _2)(M_1\mathrm{\Phi })}{\omega (\omega \alpha _2)+\alpha _1}$$

$`(B.23)`$
$$<<[A,H]|A^{}>>_\omega =\frac{\omega (M_1\mathrm{\Phi })+\alpha _1M_0}{\omega (\omega \alpha _2)+\alpha _1}$$

$$\alpha _1M_0+\alpha _2M_1=M_2$$

It is evident that there is similarity between the obtained solutions and the moment expansion method. The structure of equation (B.22) corresponds to the moment expansion (B.11) except for the factor $`\mathrm{\Phi }`$ that should be calculated by considering high-order equations of motion or by some relevant approximation.

## Appendix C Appendix . Projection methods and IGFs

The IGFs method is intimately related to the projection operator method , , that incorporates the idea of ”reduced description” of a system in the most suitable form. The projection operation  makes it possible to reduce the infinite hierarchy of coupled equations to a few relatively simple equations that ”effectively” take into account the essential information about the system that determines the specific nature of the given problem. Projection techniques become standard in the study of certain dynamic processes. Projection operator techniques of Mori-Zwanzig and similar ones  are useful for the derivation of relaxation equations and formulas for transport coefficients in terms of microscopic properties.
This approach was applied to a large variety of phenomena concerning the line-shape problem. It was shown that there is a close relationship between the Mori procedure and the ”classical moment problem” of mathematical analysis.
Let us briefly consider the projection formalism for double-time retarded GFs . Ichiyanagi  constructed the following set of equations for GF (28):

$`(C.1)`$
$$(\frac{d}{dt}i\omega _k)<<A_k(t),A_k^{}(t^{})>>=i\delta (tt^{})<[A_k,A_k^{}]>+F(k,tt^{})$$

$`(C.2)`$
$$(\frac{d}{dt}+i\omega _k)F(k,tt^{})=+i\delta (tt^{})<[K(k),A_k^{}]>+\mathrm{\Pi }(k,tt^{})$$

where $`F(k,tt^{})=<<K(k,t),A_k^{}(t^{})>>`$ and $`\mathrm{\Pi }(k,tt^{})=<<K(k,t),K^{}(k,t^{})>>`$ Here, the definitions were introduced:

$`(C.3)`$
$$i\omega _k=\frac{<[\frac{d}{dt}A_k,A_k^{}]>}{<[A_k,A_k^{}]>},K(k,t)=(1P)A_k(t)$$

$`(C.4)`$
$$PG=<[G,A_k^{}]><[A_k,A_k^{}]>^1A_k$$

The projection operator $`P`$ defined in (C.4) is different from the one introduced by Mori. The main result of paper  is that, using the projection operator, a Dyson equation that determines an irreducible quantity, proper self-energy part, was obtained in the following form:

$`(C.5)`$
$$\left(\omega \omega _k\frac{2\pi }{<[A_k,A_k^{}]>}M(k,\omega )\right)<<A_k|A_k^{}>>_\omega =\frac{<[A_k,A_k^{}]>}{2\pi }$$

Here $`M(k,\omega )`$ is the self-energy , that, in the diagrammatic language, consists of irreducible diagrams.
Our point of view is closely related to that of ref.  and to the development of the this paper by Tserkovnikov in a systematic way . However, our strategy is slightly different in the time evolution aspect. We consider our IGF technique as more convenient from the practical computational point of view.

## Appendix D Appendix . Effective Perturbation Expansion for the Mass Operator

Let us consider a useful example how to iterate the initial ”trial” solution and to get an expansion for the mass operator,. To be concrete, let us consider the calculation of the mass operator for the Hubbard model in Section 8.1. The first iteration for the equation (143) with the trial function (144) have lead us to the expression (145), which we rewrite here in the following form

$`(D.1)`$
$$M_{k\sigma }(\omega )=\frac{U^2}{N^2}\underset{pq}{}\frac{N_{kpq}}{\omega \mathrm{\Omega }_{kpq}}$$

where

$$N_{kpq}=n_{p+q\sigma }(1n_{k+p\sigma }n_{q\sigma })+n_{k+p\sigma }n_{q\sigma }$$

$$\mathrm{\Omega }_{kpq}=ϵ(p+q\sigma )+ϵ(k+p\sigma )+ϵ(q\sigma )$$

Now we are able to calculate the spectral weight function $`g_{k\sigma }(\omega )`$ (74)

$`(D.2)`$
$$g_{k\sigma }(\omega )=\frac{1}{\pi }\frac{\mathrm{\Gamma }_{k\sigma }(\omega )}{[\omega E_{k\sigma }]^2+\mathrm{\Gamma }_{k\sigma }^2(\omega )}$$

We approximate this expression by the following way

$`(D.3)`$
$$g_{k\sigma }(\omega )(1\alpha _{k\sigma })\delta (\omega E_{k\sigma })+\frac{1}{\pi }\frac{\mathrm{\Gamma }_{k\sigma }(\omega )}{[\omega E_{k\sigma }]^2}$$

Here

$$\mathrm{\Gamma }_{k\sigma }(\omega )=\pi \frac{U^2}{N^2}\underset{pq}{}N_{kpq}\delta (\omega \mathrm{\Omega }_{kpq})$$

$$E_{k\sigma }=ϵ(k\sigma )+\mathrm{\Delta }_{k\sigma }$$

$$\mathrm{\Delta }_{k\sigma }=ReM_{k\sigma }(\omega +iϵ)$$

The unknown factor $`(1\alpha _{k\sigma })`$ is determined by the normalization condition

$$_{\mathrm{}}^{\mathrm{}}𝑑\omega g_{k\sigma }(\omega )=1$$

whence

$$\alpha _{k\sigma }=\frac{U^2}{N^2}\underset{pq}{}\frac{N_{kpq}}{\mathrm{\Omega }_{kpq}E_{k\sigma }}$$

Then, using (20), we find for the mean occupation numbers

$`(D.4)`$
$$n_\sigma =\frac{1}{N}\underset{k}{}n(E_{k\sigma })+\frac{U^2}{N^3}\underset{kpq}{}\frac{N_{kpq}}{(\mathrm{\Omega }_{kpq}E_{k\sigma })^2}[n(\mathrm{\Omega }_{kpq})n(E_{k\sigma })]$$

Now we can use the spectral weight function ( D.2) to iterate the equation (143) and to get a perturbation expansion for the self-energy $`M_{k\sigma }`$ in the pair approximation. Instead of the initial trial solution in the form of delta-function (144), we take the expression (D.3). It is easy to check that we get an expansion up to 6th order in U.