Green's Function Formalism for Highly Correlated Systems

We present the Composite Operator Method (COM) as a modern approach to the study of strongly correlated electronic systems, based on the equation of motion and Green's function method. COM uses propagators of composite operators as building blocks at the basis of approximate calculations and algebra constrains to fix the representation of Green's functions in order to maintain the algebraic and symmetry properties.


Introduction
The Green's function method is a very convenient formalism in condensed matter physics, and many progresses have been achieved in the last fifty years. When applied to interacting systems, such an approach is usually based on the hypothesis that the interaction among the particles is weak and can be treated in the framework of some perturbation schemes. In this line of thinking a consolidated scheme has been constructed, mostly based on diagrammatic expansions, Wick's theorem, Dyson equation, and so on. However, in the last few decades new materials with unconventional properties have been discovered. It is believed that the origin of such anomalous behaviors is generally due to strong electronic correlations in narrow conduction bands [1]. In this line of thinking many analytical methods have been developed for the study of strongly correlated electron systems [2]. The main difficulties are connected with the absence of any obvious small parameter in the strong coupling regime and with the simultaneous presence of itinerant and atomic aspects. The concept that breaks down is the existence of the electrons as particles with some well-defined and intrinsic properties. The presence of interaction modifies the properties of the particles: what are observed are new particles with new peculiar properties entirely determined by the dynamics and by the boundary conditions. These new objects appear as the final result of the modifications imposed by the interactions on the original particles and contain, by the very beginning, the effects of correlations. The choice of new fundamental particles, whose properties have to be self-consistently determined by dynamics, symmetries and boundary conditions, becomes relevant.
As a simple example, let us consider an atomic system described by the Hamiltonian ϕ σ denotes an Heisenberg electronic field with spin σ =↑, ↓, satisfying canonical anticommutation relations; µ is the chemical potential and V is the strength of the interaction. This model is exactly solvable in terms of the operators which are eigenoperators of the Hamiltonian Due to the presence of the interaction, the original electrons ϕ σ are no more observables and new stable elementary excitations, described by the field operators ξ and η, appear. Due to the V -interaction, two sharp features develop in the band structure: the energy level E = −µ of the bare electron splits in the two levels E 1 = −µ and E 2 = V − µ. The bare electron reveals itself to be precisely the wrong place to start. A perturbative solution will never give the band splitting.
On the basis of this evidence one can be induced to move the attention from the original fields to the new fields generated by the interaction. The operators describing these excitations, once they have been found, can be written in terms of the original ones and are known as composite operators.
The convenience of developing a formulation to treat composite excitations as fundamental objects has been noticed for the many-body problem of condensed matter physics since long time. Recent years have seen remarkable developments in many-body theory in the form of an assortment of techniques that may be termed composite particle methods. The beginnings of these types of techniques may be traced back to the work of Bogolubov [3] and later to that of Dancoff [4]. The work of Zwanzig [5], Mori [6] and Umezawa [7] has to be mentioned. Closely related to this work is that of Hubbard [8][9][10], Rowe [11], Roth [12] and Tserkovnikov [13,14]. The slave boson method [15][16][17], the spectral density approach [18,19] and the composite operator method (COM) [20][21][22][23][24][25][26][27] are also along similar lines. This large class of theories is founded on the conviction that an analysis in terms of elementary fields might be inadequate for a system dominated by strong interactions.
All these approaches are very promising because all the different approximation schemes are constructed on the basis of interacting particles: some amount of the interaction is already present in the chosen basis and permits to overcome the problem of finding an appropriate expansion parameter. However, one price must be paid. In general, the composite fields are neither Fermi nor Bose operators, since they do not satisfy canonical (anti)commutation relations, and their properties, because of the inherent definition, must be self-consistently determined. They can only be recognized as fermionic or bosonic operators according to the number, odd or even, of the constituting original electronic fields. New techniques of calculus have to be developed in order to treat with composite fields. In developing perturbation calculations, where the building blocks are now the propagators of composite fields, the consolidated scheme (diagrammatic expansions, Wick's theorem, . . . ) gives rise to very complicated approaches [28] whose application is far from being easy. The formulation of the Green's function method must be revisited and new frameworks of calculations have to be formulated.

Green's Function and Equation of Motion Formalism
Let us consider a system of N e interacting Wannier-electrons residing on a Bravais lattice of N sites, spanned by the vectors R i = i. We ignore the presence of magnetic impurities and restrict the analysis to single-band electron models. The generalization of the formalism to more complex systems is straightforward [see for example [29][30][31][32][33]]. The system is enclosed in a finite but macroscopically large volume V and is supposed to be in a thermodynamical equilibrium state at a finite temperature T. In a second quantization scheme this system is described by a certain Hamiltonian describing, in complete generality, the free propagation of the electrons and all the interactions among them and with external fields. ϕ(i) denotes an Heisenberg electronic field [ i = (i, t)] satisfying canonical anticommutation relations derived from the Pauli principle Any physical property of the system can be connected to the expectation value of a specific operator. The expectation value of an arbitrary operator A = A[ϕ(i)] can be computed, for the grand canonical ensemble, by means of where the trace implies a sum over a complete set of states in the Hilbert space.
is the total number operator, β = (k B T ) −1 is the inverse temperature, µ is the chemical potential which is fixed in order to give the chosen average number of particles N e = N .
To evaluate the expectation value A , it is possible to use the equation of motion in order to derive one or more equations for this quantity or, better, for the corresponding Green's functions, as explained below. However, the equation of motion (2.4) generates higher-order operators and more and more complex equations are needed. The traditional approximation schemes, often based on perturbative calculations, use as building blocks the noninteracting propagators. The mean-field formulation, which corresponds to a linearization of the equation of motion (2.4), also belongs to this category. On the hypothesis that the original fields are not a good basis, we choose a set of composite fields {ψ(i)} in terms of which a perturbation scheme will be constructed. Firstly, we choose the set ψ(i) according to the physical properties we want to study. Roughly, the properties of electronic systems can be classified in two large classes: single particle properties, described in terms of fermionic propagators, and response functions, described in terms of bosonic propagators. These two sectors, fermionic and bosonic, are not independent but interplay each other, and a fully self-consistent solution usually requires that both sectors are simultaneously solved. Once the sector, fermionic or bosonic, has been fixed, we have several criteria for the choice of the new basis. In constructing the composite fields no recipe can be given without thinking to its drawbacks, but many recipes can assure a correct and controlled description of relevant aspects of the dynamics. One can choose: the higher order fields emerging from the equations of motion (i.e., the conservation of some spectral moments is assured), the eigenoperators of some relevant interacting terms (i.e., the relevant interactions are treated exactly), the eigenoperators of the problem reduced to a small cluster, . . .
Let ψ(i) be a n-component field We do not specify the nature, fermionic or bosonic, of the set {ψ(i)}. In the case of fermionic operators it is intended that we use the spinorial representation The dynamics of these operators is governed by the given Hamiltonian H = H[ϕ(i)] and can be written as It is always possible to decompose the source J(i) under the form where the linear term represents the projection of the source on the basis {ψ} and is calculated by means of the equation Here η = ±1; usually, it is convenient to take η = 1 (η = −1 ) for a fermionic (bosonic) set ψ(i) (i.e., for a composite field constituted of an odd (even) number of original fields) in order to exploit the canonical anticommutation relations of {ψ(i)}; but, in principle, both choices are possible. Accordingly, we define · · · denotes the quantum statistical average over the grand canonical ensemble, according to Eq. (2.3). Hereafter, unless otherwise specified, time and space translational invariance will be considered. The action of the derivative operator ε(−i∇) on ψ(i) is defined in momentum space where k runs over the first Brillouin zone. The constraint (2.9) gives after defining the normalization matrix and the m-matrix Since the components of ψ(i) contain composite operators, the normalization matrix I(k) is not the identity matrix and defines the spectral content of the excitations In fact, the use of composite operators has the advantage of describing crossover phenomena as the phenomena in which the weight of some operator is shifted to another one.
It is worth noting that the normalization matrix I(k) and the m(k) matrix are the lowest order generalized spectral moments M (p) (k) which are defined as where Coming back to our original problem, the evaluation of the expectation value A , it is possible to use the equation of motion i ∂ ∂t ψ(i) = [ψ(i), H] in order to derive equations for A . However, the correlation functions satisfy homogeneous equations. A convenient generalization of the concept of correlation functions is furnished by the Green's functions (GF) which have some advantages in the construction and solution of the equations determining them. In particular, the two-time Green's functions contain most of the relevant information on the properties of the system: expectation values of the observables, excitation spectrum, response to external perturbation, and so on. Different types of GF can be defined; for statistical systems it is better to consider the real time thermodynamic Green's functions where the averaging process of the Heisenberg operators is performed over the grand canonical ensemble.
By considering the two-time thermodynamic Green's functions [34][35][36], let us define the causal function .17) the retarded and advanced functions By means of the Heisenberg equation (2.7) and using the decomposition (2.8), the Green's function where the derivative operator Λ(∂ i ) is defined as and the propagator G Q 0 (i, j) is defined by the equation By introducing the Fourier transform equation (2.19) in momentum space is written as where the self-energy Σ Q * (k, ω) has the expression Next, we introduce the irreducible self-energy Σ Q (k, ω) by means of the definition and can be formally solved as The formal definition (2.24) of self-energy must be manipulated to avoid any flowing on tautology. By noting that We have constructed a generalized perturbative approach designed for formulations using composite fields. Equation (2.27) is a Dyson-like equation and may represents the starting point for a perturbative calculation in terms of the propagator G Q 0 (k, ω). Contrarily to the usual perturbation schemes, the calculation of the "free propagator" G Q 0 (k, ω) is not an easy task and large part of this article will be dedicated to this problem. Then, the attention will be given to the calculation of the self-energy Σ Q (k, ω), and some approximate methods will be presented. It should be noted that the computation of the two quantities G Q 0 (k, ω) and Σ Q (k, ω) are intimately related. The total weight of the self-energy corrections is bounded by the weight of the residual source operator δJ(i). According to this, it can be made smaller and smaller by increasing the components of the basis ψ(i) [e.g. by including higher-order composite operators appearing in δJ(i)]. The result of such a procedure will be the inclusion in the energy matrix of part of the self-energy as an expansion in terms of coupling constants multiplied by the weights of the newly includes basis operators. In general, the enlargement of the basis leads to a new self-energy with a smaller total weight. However, it is necessary pointing out that this process can be quite cumbersome and the inclusion of fully momentum and frequency dependent self-energy corrections can be necessary to effectively take into account low-energy and virtual processes. According to this, one can chose a reasonable number of components for the basic set and then use another approximation method to evaluate the residual dynamical corrections.

Equations of motion
In the previous Section we have constructed a generalized perturbative approach based on a Dyson equation designed for formulations using composite fields. Two quantities appear in the Dyson equation (2.27): the "free propagator" G Q 0 (k, ω) and the self-energy Σ Q (k, ω). By postponing to next Sections the problem of computing the self-energy, in this Section we concentrate on the calculation of the Green's functions G Q 0 (k, ω) which constitute the building blocks of the perturbation scheme we are trying to formulate. For the sake of simplicity we will drop the sub index 0 in the definition of G Q 0 (k, ω). One fundamental aspect in a Green's function formulation is the choice of the representation. The knowledge of the Hamiltonian and of the operatorial algebra is not sufficient to completely specify the GF. The GF refer to a specific representation (i.e., to a specific choice of the Hilbert space) and this information must be supplied as a boundary condition to the equations of motion that alone are not sufficient to completely determine the GF. As well known, the same system can exist in different phases according to the external conditions; the existence of infinite inequivalent representations [37][38][39] where the equations of motions can be realized, allows us to pick up, among the many possible choices, the right Hilbert space appropriate to the physical situation under study. The use of composite operators leads to an enlargement of the Hilbert space by the inclusion of some unphysical states. As a consequence of this, it is difficult to satisfy a priori all the sum rules and, in general, the symmetry properties enjoined by the system under study. In addition, since the representation where the operators are realized has to be dynamically determined, the method clearly requires a process of self-consistency.
From this discussion it is clear that fixing the representation is not an easy task and requires special attention. In the literature the properties of the GF are usually determined by starting from the knowledge of the representation. Owing to the difficulties above discussed we cannot proceed in this way. Therefore, we will derive the general properties of the GF on the basis of the two elements we have: the dynamics, fixed by the choice of the Hamiltonian (2.1), and the algebra, fixed by the choice of the basic set (2.5). The problem of fixing the representation will be considered in the next Sections.
Let ψ(i) be a n-component field satisfying linear equations of motion with the energy matrix ε(i, j) defined by (2.12)-(2.14). If the fields ψ(i) are eigenoperators of the total Hamiltonian, the equations of motion (3.1) are exact. There are many non-trivial realistic systems for which it is possible to obtain a complete set of eigenoperators of the Hamiltonian (for instance see [40][41][42][43] where the dependence on the parameter η has been explicitly introduced. As mentioned in Section 2, the set ψ(i) can be fermionic or bosonic and the parameter η generally takes the value η = 1 (η = −1) for a fermionic (bosonic) set ψ(i). The three Green's functions G C , G R and G A satisfy the same equation of motion which alone is not sufficient and must be supplemented by other equations. Indeed, the GF are determined by solving a first order differential equation of motion, thereby the GF are given only within an arbitrary constant of integration. The retarded and advanced GF can be completely determined because the factor θ[±(t i − t j )] provides the boundary condition: The determination of the causal GF is not so immediate. The most general solution of equation (3.2) is where ω l (k) are the eigenvalues of the ε(k), σ (l) (k) are defined by Ω(k) is the n × n matrix, whose columns are the eigenvectors of ε(k). We note that while the spectral density matrix σ (l) (k) is completely determined by the energy ε(k) and normalization I (η) (k) matrices, the matrix g (l,η)Q (k) is not fixed by the equations of motion and must be determined by means of the boundary conditions. P represents the principal value. By recalling the retarded and advanced nature of Then, the retarded and advanced functions are completely determined in terms of the matrices ε(k) and I (η) (k). As well known, as functions of ω the G R,A(η) (k, ω) are analytic in the upper and lower half-planes, respectively. The determination of g (l,η)C (k) requires more work. From the definitions (2.16), (2.17) and (2.18) we can derive the following exact relations where there appear the two correlation functions The definition of C ψ † ψ (i, j) must be interpreted at level of matrix elements: C ψ † ψ;ba (i, j) = ψ † b (j)ψ a (i) ; no scalar product is intended, neither in the spin space. These functions are defined for all real times and the equations of motion they obey contain no inhomogeneous terms involving a delta function. Indeed, by means of the equations of motion (3.1) the Fourier transform of these correlation functions satisfy the homogeneous equations These equations tell us that the Fourier transforms of the correlation functions are zero unless the frequency ω is equal to one of the energy levels ω l (k) of the system. The solutions of (3.9) have the general form with the momentum-dependent Fourier components c (l) ψψ † (k) and c (l) ψ † ψ (k) to be determined.
We now recall the Kubo-Martin-Schwinger (KMS) relation where A(t) and B(t) are Heisenberg operators at time t. This relation implies that the Fourier transforms C ψψ † (k, ω) and C ψ † ψ (k, ω) are related and the η−commutator [ψ(i), ψ † (j)] η can be expressed in terms of the correlation function as By putting (3.14) into (3.6) and (3.7) and by taking into account Eqs.
The solution of Eqs. (3.15) and (3.16) is remarkably different according to the value of the parameter η and we shall treat separately the two cases.

Fermionic fields
For the case of fermionic fields it is convenient to choose η = 1. Then, the solution of Eqs. (3.15) and (3.16) is where f F (ω) is the Fermi distribution function: f F (ω) = 1 e βω +1 . By recalling that all the fermionic energies are shifted by the chemical potential, the locus in the k-space, defined by ω l (k) = 0, will define the Fermi surface. By looking at (3.15b) we can see that the imaginary part of the causal GF vanishes on the Fermi surface. In the fermionic case the right procedure of calculation is to start from the retarded (advanced) GF, and then to compute the other Green's functions and correlation functions by means of the relations We note the dispersion relations By introducing the spectral function we can establish the spectral representation

Bosonic fields
For the case of bosonic fields it is convenient to choose η = −1. For any given momentum k we can always write Obviously, A(k) can also be the empty set (i.e., A(k) = ∅ and B(k) = N). For l ∈ B(k) the solution of (3.15) and (3.16) is For l ∈ A(k) the Fourier coefficients c (l) ψψ † (k) and g (l,−1)C (k) cannot be determined from Eqs. (3.16). It is convenient to introduce the function Γ(k) By considering that from (3.16) and (3.10) we must distinguish two situations: 2. If l∈A(k) σ (l,−1) (k) = 0 but finite, then in order to satisfy (3.35) C ψψ † (k, ω) must have a singularity of the type 1 ω in the limit ω → 0. In fact It is clear from (3.40) and (3.41) that the situation where l∈A(k) σ (l,−1) (k) = 0 leads to a situation in which for l ∈ A(k) the Fourier coefficients c (l) ψψ † (k) and c (l) ψ † ψ (k) diverge as [βω l (k)] −1 . Since the correlation function in direct space must be finite, at finite temperature this is admissible only in the thermodynamic limit and if the dispersion relation ω l (k) is such that the divergence in momentum space is integrable and the corresponding correlation function in real space remains finite. For finite systems and for infinite systems where the divergence is not integrable we must have l∈A(k) σ (l,−1) (k) = 0. The calculation of the spectral density matrices σ (l,−1) (k) it not a simple dynamical problem, but requires the self-consistent calculation of some expectation values, where the boundary condition and the choice of the representation play a crucial role. A finite value of l∈A(k) σ (l,−1) (k) is generally related to the presence of long-range order and the previous statement is nothing but the Mermin-Wagner theorem [44].
Summarizing, by using (3.5) and by putting (3.36) and (3.37) into (3.15), for finite systems and T = 0, we get the following general expressions for the GF and correlations function It is possible to have σ (l,−1) (k) = 0 f or l ∈ A(k) only for infinite systems and if the divergence is integrable.
The previous formulas show that when zero-energy modes are present a zerofrequency singularity appears in the correlation function C ψψ † (k, ω) and in the imaginary part of the causal function G C(−1) (k, ω). Such singularity does not contribute to the retarded and advanced GF. Then, in the bosonic case the right procedure of calculation is to start from the causal GF, and compute the other GF by means of the relations We note the dispersion relations Also in the bosonic case we can introduce a spectral function However, the zero-frequency function Γ(k) does not contribute to ρ (−1) (k, ω) and a spectral representation can be established only for the retarded (advanced) GF For the bosonic causal GF a spectral representation exists only when Γ(k) = 0:

Sum rules and some useful relations
Coming back to a generic value of Γ(k), we note that from the definition (3.4) the following sum rule can be derived This is a particular case of a general sum rule. From the results (3.21), (3.22), (3.45) and (3.46) we obtain Recalling the expression (2.15) of the generalized spectral moment M (p,η) (k) we immediately have Some interesting results can be obtained by noting that the correlation function C ψψ † (i, j) = ψ(i)ψ † (j) and the energy matrix ε(i, j) do not depend on η. As mentioned above, once we have chosen a basic set {ψ(i)}, fermionic or bosonic, it is a only a question of convenience to chooseη = 1 or η = −1. Let us consider the case of a bosonic set {ψ(i)} and let us suppose to perform two series of calculations: one with η = −1 and another with η = 1. Then, it is immediate to obtain the following relations: 1. by equating (3.21) and (3.45) we obtain We see that the general structure of the GF is remarkably different according to the statistics. For fermionic composite fields (i.e., when it is natural to choose η = 1) all the Green's functions and correlation functions are completely determined. The zero-frequency function Γ(k), defined on the Fermi surface ω l (k) = µ, contributes to the spectral function ρ (+1) (k, ω) (see 3.28), it is directly related to the spectral density functions σ (l,+1) (k) by means of equation (3.58), and its calculation does not require more information. Also, it does not contribute to the imaginary part of the causal GF. For bosonic composite fields (i.e., when it is natural to choose η = −1) the retarded and advanced GF are completely determined, but the causal GF and the correlation function depend on the zero-frequency function Γ(k), defined on the surface ω l (k) = 0. It is now clear that the causal and retarded (advanced) GF contain different information and that the right procedure of calculation is controlled by the statistics. In particular, in the case of bosonic fields (i.e., for η = −1) one must start from the causal function and then use (3.34) to compute the other GF. On the contrary, for fermionic fields (i.e., for η = 1) the right procedure for computing the correlation function requires first the calculation of the retarded (advanced) function and then the use (3.23), (3.24) and (3.25) to compute the other GF. Moreover, it is worth noting that Γ(k) is undetermined within the bosonic sector (i.e., η = −1). It is true that Γ(k) could be computed by considering an anticommutating algebra: remaining in the bosonic sector we make the choice η = 1 and Γ(k) can be calculated by (3.58) or equivalently by means of the following relation Γ(k) = 1 2 lim ω→0 ωG C(+1) (k, ω) which can be easily obtained from (3.44). However, the calculation of the σ (l,+1) (k) requires the calculation of the normalization matrix I (+1) (k) that, for bosonic fields, generates unknown momentum dependent correlation functions whose determination can be very cumbersome as requires, at least in principle, the self-consistent solution of the integral equations connecting them to the corresponding Green's functions. In practice, also for simple, but anyway composite, bosonic fields the Γ(k) remains undetermined and other methods should be used. Similar methods, like the use of the relaxation function [45], would lead to the same problem.
Actually, all issues related to Γ(k) have a natural playground in dealing with the ergodicity of the dynamics under investigation. More detail on this topic can be found in a manuscript in this same issue [46].
The formulation given in this Section needs some modifications in the case of zero temperature. In particular, Eqs. (3.15) and (3.16) are not applicable and we must proceed in a different way. After a straightforward derivation [26], it is immediate to see that the limit T → 0 of the expressions (3.21), (3.22), (3.45) and (3.46) gives the right result.

A self-consistent scheme
As stressed in Section 1, in the study of highly interacting systems, where traditional perturbative calculations in terms of the noninteracting fields fail, a way to reconcile the powerful perturbation theory with the presence of complex and/or strong interactions is to describe the system in terms of a new set of fields, composite operators, generated by the interactions themselves. These field operators undoubtedly constitute a better starting point: they appear as the final result of the modifications imposed by the interactions on the original particles and contain, by the very beginning, the effects of the correlations. Once a choice of composite fields has been made, the relative Green's function formalism can be set up, as illustrated in Section 2, where a generalized Dyson equation has been derived. On this basis one can construct a non-standard perturbation formalism where the basic ingredients are the propagators of a subset of the fundamental basis, satisfying linear equations of motion [cfr. (3.1)]. By means of the equations of motion and by using the boundary conditions related to the definitions of the various Green's functions we have been able to derive explicit expressions for these latter [cfr. (3.21), (3.22), (3.45) and (3.46)]. However, these expressions can only determine the functional dependence; the knowledge of the GF is not fully achieved yet. The reason is that the algebra of the field is not canonical. As a consequence, the inhomogeneous terms I (η) (k) in the equations of motion (3.2) and the energy matrix ε(k) contain some unknown static correlation functions, correlators, that have to be self-consistently calculated. Three serious problems arise with the study of the Green's functions: (a) the calculation of some parameters expressed as correlation functions of field operators not belonging the chosen basis; (b) the appearance of some zero-frequency constants (ZFC) and their determination; (c) the problem of fixing the representation where the Green's functions are formulated.
In the Composite Operator Method [26,27] (COM) the three problems (a), (b) and (c) are not considered separately but they are all connected in one self-consistent scheme. The main idea is that fixing the values of the unknown parameters and of the ZFC implies to put some constraints on the representation where the GF are realized. As the determination of this representation is not arbitrary, it is clear that there is no freedom in fixing these quantities. They must assume values compatible with the dynamics and with the right representation. Which is the right representation? From the algebra it is possible to derive several relations among the operators. We will call algebra constraints (AC) all possible relations among the operators dictated by the algebra. This set of relations valid at microscopic level must be satisfied also at macroscopic level, when expectations values are considered. Also, we note that, in general, the Hamiltonian has some symmetry properties (i.e. rotational invariance in coordinate and spin space, phase invariance, gauge invariance,. . . ). These symmetries generate a set of relations among the matrix elements: the Ward-Takahashi identities [47,48] (WT). Now, certainly the right representation must be the one where the relations among the operators and the conservation laws are maintained when expectation values are taken; in other words, all the AC and WT are preserved. By imposing these conditions we obtain a set of self-consistent equations that will fix the unknown correlators, the ZFC and the right representation at the same time. Several equations can be written down, according to the different symmetries we want to preserve. A large class of self-consistent equations is given by the following equation where the l.h.s. is fixed by the AC, the WT and the boundary conditions compatible with the phase under investigation and in the r.h.s. the correlation function C ψψ † (k, ω) is computed by means of the equation of motion, as illustrated in Section 3. Equations (4.1) generate a set of self-consistent equations which determine the unknown parameters (i.e., ZFC and unknown correlators) and, consequently, the proper representation [26,27], avoiding the problem of uncontrolled and uncontrollable decoupling.

Approximation schemes
The generalized Dyson equation (2.27) is an exact equation and permits, in principle, once the normalization matrix I(i, j) [cfr. Eq. (2.13)], the m-matrix m(i, j) [cfr. Eq. (2.14)] and the propagator B(i, j) [cfr. Eq. (2.24)] are known, in the framework of the self-consistent scheme outlined in Sections 3 and 4, the calculation of the various Green's functions. However, for most of the physical systems of interest the calculation of the propagator B(i, j) is a very difficult task and some approximations are needed. Various approximate schemes have been proposed. We will summarize some of them.

The n-pole approximation
The simplest approximation is based on completely neglecting the dynamical part Σ(k, ω) [cfr. (2.28)] of the self-energy. This approximation is largely used in the literature [2, 6, 8-14, 19-25, 49-62] and is called pole-approximation. In this approximation we only need the knowledge of the normalization matrix and the m-matrix. The constraint (2.9) produces a physics of the solution totally extraneous to the complementary physical space, orthogonal to that spanned by the multiplet ψ(i). This approximation, or assumption to a larger extent, consists in retaining that one can neglect finite life-time effects (i.e., the dynamical part of the self-energy) paying attention to the choice of a proper extended operatorial basis, with respect to which the self-energy corrections have a small total weight. Indeed, the total weight of the corrections is bounded by the thermal average (2.25) involving the residual source δJ(i). It is worth noting [63] that the n-pole structure of the various GF corresponds to a Dyson-like equation where the self-energy components Σ Q ab (k, ω) have a (n − 1)-pole structure.

Self-consistent Born approximation
In order to improve the approximation one needs to take into account selfenergy corrections by developing some methods to calculate the effects of Σ(k, ω). In the self-consistent Born approximation (SCBA), or non-crossing approximation, the many-particle Green's functions, appearing in the expression of Σ(k, ω) [see (2.29)], are calculated by assuming that the fermionic and bosonic modes propagate independently.
By recalling the results given in Section 2, the knowledge of the self-energy requires the calculation of the higher-order propagator B Q (i, j) = Q[δJ(i)δJ † (j)] . In order to illustrate the approximation, let us consider the case where the basic set {ψ(i)} is of a fermionic type. Then, typically we have to calculate GF of the form where F (i) and B(i) are fermionic and bosonic field operators, respectively. By means of the spectral representation (3.29) we can write where H C (i, j) = T [B(i)F (i)F † (j)B † (j)] is the causal function. In the SCBA we approximate This approximation has been used in many works (for instance see [64][65][66][67]. By assuming that the system is ergodic we can use the spectral representations (3.30) and (3.54) to obtain Use of (5.2)-(5.4) leads to where d is the dimensionality of the system, a is the lattice constant and Ω B is the volume of the Brillouin zone.
It should be noted that the SCBA can also be applied to the correlation function. We start from the expression where H(i, j) = B(i)F (i)F † (j)B † (j) is the correlation function. In the SCBA we approximate H(i, j) = B(i)F (i)F † (j)B † (j) ≈ F (i)F † (j) B(i)B † (j) (5.8) Then, by proceeding in the similar way we arrive to the same expression (5.6).

Two-site resolvent approach
In this subsection we will consider an approximation scheme [68,69], where the dynamical part Σ(k, ω) [cfr. (2.29)] of the self-energy is estimated by a two-site approximation in combined use with the resolvent method [70]. In this approximation the higher order propagator (2.25) is written as where B Q 0 (ω) is related to level transitions on equal site The Green's function (2.28) takes the form G Q (k, ω) = 1 ω − ε(k) + t 2 V (ω)α(k) I(k) (5.12) where V (ω) has to be calculated from the definition (2.29) by making use of approximation (5.9). We give now a brief sketch of the calculation. B Q 0 (ω) and B Q 1 (ω) are computed by expressing δJ(i) in terms of transitions among the two-site levels δJ = nm a nm Φ † n Φ m (5.13) where {Φ m } is the complete set of operators for two-site levels. By means of the non crossing approximation [70,71], the propagator Q[Φ † n (t i )Φ m (t i )Φ † n (t j )Φ m (t j )] is expressed in terms of the resolvent R nm (t i − t j ) = Q[Φ n (t i )Φ † m (t j )] R where the subscript R indicates the reservoir system, i.e., the part H R of the Hamiltonian other than that concerned with two sites. The calculation of the resolvent brings to a modification of the original two-site levels by the surroundings. In this scheme effects of time delay in the local correlations are treated trough the time-dependent modifications of the two-site level transitions, which are included as time-dependent local effects in the electron self-energy.

Conclusions
In this article we have illustrated an approach to the study of highly correlated electronic systems, based on the equation of motion and Green's function method. Such an approach is based on two main ideas: (i) propagators of composite operators as building blocks at the basis of approximate calculations; (ii) use of algebra constrains to fix the representation of the GF in order to maintain the algebraic and symmetry properties. This formalism has been applied to the study of several models of highly interacting systems, and we refer the interested readers to Ref. [27] for an exhaustive bibliography.