Equation of motion truncation scheme based on partial orthogonalization

We introduce a general scheme to consistently truncate equations of motion for Green's functions. Our scheme is guaranteed to generate physical Green's functions with real excitation energies and positive spectral weights. There are free parameters in our scheme akin to mean field parameters that may be determined to get as good an approximation to the physics as possible. As a test case we apply our scheme to a two-pole approximation for the 2D Hubbard model. At half-filling we find an insulating solution with several interesting properties: it has low expectation value of the energy and it gives upper and lower Hubbard bands with the full non-interacting bandwidth in the large U limit. Away from half-filling, in particular in the intermediate interaction regime, our scheme allows for several different phases with different number of Fermi surfaces and topologies.


I. INTRODUCTION
Green's function methods are widely used to study many-body systems and they represent a natural framework that connects microscopical details of a theory with its macroscopical properties. 1 The attempt to self-consistently determine these quantities has a long history and it is still remains one of the central paradigms in the study of strongly correlated systems. From the most recent DMFT, 2 where they are used to fix the mapping of a lattice model onto an impurity one, to the older equation of motion approach. 3,4 In the latter method, given an interacting Hamiltonian, an extensively growing chain of coupled equations are derived. 5,6 For few-body systems it is possible to use various implementations of this method to obtain the single particle Green's function exactly. 7 However, in order to study thermodynamical properties of an interacting system a truncation procedure able to approximately decouple this extensively growing system of coupled equations plays a crucial role. Early attempts in the construction of truncation schemes explored arbitrary truncation schemes and decoupling schemes of Tyablikov-type. 3,4 Despite some successful applications these decoupling schemes often led to violation of the analytical structure of the Green's functions, predicting imaginary poles and negative spectral weight for the single particle Green's function. Despite these difficulties Hubbard in his pioneering work, 8 managed to find a useful decoupling for a two-pole approximation for the Hubbard model. This decoupling (Hubbard-I) is still often used in treating strongly correlation in presence of local interactions, especially in studies of quantum systems out of equilibrium, 9 and multi-orbital systems. 10 Almost a decade after these early works Roth developed a universal decoupling scheme able to enforce correct analytical properties for approximated Green's functions. 11 This decoupling scheme is now called the Roth procedure, and often relies on parameters that can not be determined within the scheme itself, making unavoidable ulterior approximations. For this reason this method is often regarded as an uncontrolled approximation, which severely limits its applicability. The works of Mancini and Avella et al. 12 show that the Roth procedure leads to violations of other physical principles such as the Pauli principle and that it is possible to constrain some, if not all of the unknown decoupling parameters, by enforcing such physical requirements. Despite much progress in finding easy extendable decoupling schemes, 13 the possibility to systematically check what are the approximations involved in the decoupling still remains a neglected aspect.
In this paper we present a decoupling scheme based on a partial orthogonalization of the operators involved, where the relation between the true Green's function and the approximate one can readily be obtained. The paper is organized as follows: in Sec. II we provide a general discussion of the formalism, we clarify the role of the Hermiticity of the E-matrix and we present our decoupling scheme based on the partial orthogonalization of the operators. In Sec. III we apply our scheme to a two-pole approximation of the Hubbard model making evident the relationship between the approximate and the true Green's function. In Sec. IV we analyze the global sum rules that should be respected in the two-pole approximation of the Hubbard model and we present a variational scheme as a guiding principle for the determination of the unknown orthogonalization parameters. In Sec. V we provide numerical results at half-filling and in Sec. VI we give analytical formulas that are useful to understand the Green's function. In Sec. VII and Sec. VIII we discuss numerical results for hole doping in the strongand intermediate-coupling regimes respectively. Finally, in Sec. IX we provide some conclusions and an outlook.

II. A SCHEME FOR THE TRUNCATION OF THE EOM
A. Formalism review We will mainly use the notation of Tserkovnikov in the following. 5 For completeness we briefly review what we will need for this paper. Let us first assume we have a set of fermionic operators {Â i } M i=1 closed under the commutation with the hamiltonian for some evolution matrix Then the equation of motion (EOM) for the Green's function matrix gives here the normalization matrix N Consequently the Green's function, viewed as a matrix becomes where the second form is obtained making use of the fact thatĤ is Hermitean. For these two forms to be consistent we have the condition that which will be of crucial importance in the developments below. Finally one may calculate averages of bilinear of all of the operators involved using the formula where the contour encircles the real axis.
B. A partial orthogonalization scheme for the truncation of the EOM As shown in the previous work this framework gives exact results if the set of operator {Â i } M i=1 are closed under the commutation with the Hamiltonian. 7 In an extend many-body system the number of operators necessary to close the equation of motion exactly will typically grow exponentially with the size of the system, making a direct application of this scheme unfeasible.
To produce a truncation scheme capable of producing physical Green's functions, it is important to notice that in a Hermitean theory the average of the operators involved in the dynamics and their evolution are not independent. In particular as noticed by Roth 11 the matrix needs to be Hermitean. Here . . . indicate some average over exact eigenstates of the theoryĤ. K is the full evolution matrix of the operators and N is the normalization matrix introduced above. Using the fact that the matrix N is Hermitean by construction (i.e., it holds for averages in any state) this gives the same consistency condition as Eq. (5) above. This condition together with the fact that N has to be positive definite guarantees that the Green's function posseses real poles and positive spectral weight. 11 When the hierarchy of the evolution of an operatorÂ 1 is considered at most one new operator is generated in each step, i.e., etc. until the EOM closes and no new operators are generated. Note thatÂ 2 is not unique since one can add a part ofÂ 1 to it, and similarly for the other higherÂ's. In any event K is only non-zero on the first upper diagonal and below. Let us now truncate the EOM at the q-th operator. A brute force truncation of the matrices involved gives and the corresponding N trunc Now we note that differs from the corresponding sub-block of the full E only in the last row, through the coupling of K q,q+1 to the (q + 1)-th column of the full N matrix. Therefore an arbitrary truncation of the equation of motion is going to generate an evolution that in general does not satisfy the condition in Eq. (5), leading to a potentially unphysical approximation for the Green's function.
In this paper we propose to restore the Hermiticity of E trunc adding to the first operator not considered explicitly in the dynamicsÂ q+1 a linear combination of the operatorsÂ 1 , . . . ,Â q Most of the λ parameters will be fixed by demanding that This partial orthogonalization procedure ensures that E trunc is Hermitean, because it makes it identical to the corresponding block of E except for the last element on the diagonal E qq which is not fixed by our procedure. The Roth procedure corresponds to also orthogonalizing with respect toÂ † q . This gives q equations for q unknowns, and therefore also fixes the value of E qq , whereas in our scheme we have q − 1 equations for q unknowns, leaving E qq arbitrary. We will use this additional freedom to make sure that our approximation fulfills other physically relevant criteria such as Pauli principle constraints or sum rules.
In the next section we are going to elucidate this procedure by applying it to a two-pole approximation of the Hubbard model. In particular it will be evident that the effect of this procedure is a non-unique modification of the last row of K trunc . This arbitrariness can be exploited to enforce global sum rules for the Green's functions and open up the possibility of using different criteria to fix the free parameters λ i not fixed by Eq. (13).
A last remark on this scheme is that despite the freedom in the choice of the parameters λ i one can always write the residual Green's functions not considered explicitly in the dynamics, making transparent the approximation involved in this truncation of the equation of motion.

III. APPLICATION TO THE HUBBARD MODEL IN A TWO-POLE APPROXIMATION
In this section we will apply our scheme to a two-pole approximation to the Green's function in the Hubbard model. Let us consider the Hubbard hamiltonian We will denote the total number of sites with N s , c kσ indicates fermion operator with spin σ and n iσ = c † iσ c iσ . Let us consider the first three operators that appear in the equation of motion hierarchŷ where we have introduced Let us first do a brute force truncation of the evolution after two operators. The truncated evolution becomes and the respective N matrix becomes withn ↓ = n i↓ which is independent of the site index i. In this case E trunc = K trunc N trunc is not Hermitean (except in special cases such as U = 0, k = 0 orn ↓ = 0) and this leads to an unphysical approximation for the Green's function for some range of the parameters.
Let us now apply our scheme to this particular problem, in this case we need to determine λ 1k , λ 2k such that Evaluating the anticommutator averages we obtain As already anticipated in Sec. II, the values of λ 1k and λ 2k are not uniquely determined by this procedure. Without any loss of generality let us eliminate λ 1k writing Using this we can writê where A 3k |c † k = 0. At this point the equation of motion for the operatorÂ 1k can be rewritten as (B is here arbitrary) and forÂ 2k Consequently the new evolution given by the partial orthogonalization procedure is The physical condition in Eq. (5) is now satisfied for this evolution for any choice of the model parameters λ 2k ,n ↓ , U, k . The approximate Green's function for the truncated theory becomes Assuming no spin symmetry breaking, the parameter n ↓ can be determined self-consistently, by applying the fermionic characterization of the spectral theorem stated in Eq. (6) to G 11 , obtaining: To see that we can always write the residual Green's function highlighting the approximation involved in the truncation of the equation of motion let us analyze the special case where λ 2k = k . The equation of motion of the Green's function with B † = c † k becomes: Recalling that A 1k = c k↑ we find that the conventional fermion Green's function may be written exactly as From this it is clear that truncating the equation of motion implies that the term on the last line is neglected, making the approximation evident. Moreover we note that A 3k |c † k↑ does not contain poles at k and k + U (since double poles in the original Green's function are not allowed) and its total spectral weight is vanishing (since A 3k |c † k↑ = 0). We may also note that excitations at k and U + k appears in the exact thermal Green's function (although their weight may be exponentially small) since they are exact energy differences between states with charge 1 and 0 and 2N s − 1 and 2N s respectively.

IV. THE GLOBAL CONSTRAINTS ON THE TWO-POLE APPROXIMATION OF THE HUBBARD MODEL
As noticed and stressed by Mancini and Avella 12 the Roth procedure does not ensure that global sum rules such as those related to the Pauli principle and Ward identities are satisfied. In the context of a two-pole approximation of the Hubbard model, these violations can be related to global constraint between averages. In particular the average double occupancy of the system can be evaluated in two inequivalent ways At the operatorial level these two ways of writing the averages are equal. Consequently when we evaluate these averages using the spectral theorem and the effective evolution, we have to make sure that This constraint is very important, because it removes a fundamental ambiguity related to the determination of the energy in the Roth scheme. In particular we can notice that in the previously studied solution, where we used λ 1k = 0 and λ 2k = k the constraint in Eq. (31) is automatically satisfied, because the argument of the integral is identically 0 for every k, making the solution suitable for unambiguous physical interpretation. On a physical level this choice of the parameters makes the evolution diagonal in the two Hubbard operators which are orthogonal by construction.
A. Variational determination of the orthogonalization parameter As previously stated the determination of the orthogonalization parameters λ 2k plays a crucial role. Different values for this parameters gives different approximations to the true Green's function, all of them are physical in the sense that the spectral weights are positive and the excitation energies real, which is a fundamental requirement. On the other hand different values of this parameter may correspond to quite different physics. In some sense λ 2k may be viewed as a kind of mean field parameter, in the sense of variational mean field theory. 14 Any choice for λ 2k is allowed and gives physical results, but we want to determine the parameter to approximate the physics in the "best" possible way. The definition of "best" is however not unique, since approximations do not get everything correctly. Depending on what one choose to optimize different approximations will result.
It may be reasonable to demand that the solution posses the full lattice symmetry (i.e., assuming unbroken lattice symmetry). Then the evolution matrix K(k) may be expanded in terms of proper basis functions with full lattice symmetry. The simplest non-trivial possibility is to take the ansatz for λ 2k to be where a 0 and a 1 are some real k-independent constants. This may be viewed as the first two terms in a locality expansion. Let us also note that this is exactly the form for λ 2k that is obtained in the Roth procedure in the two-pole approximation in the Hubbard model. 15 If we further assume unbroken spin symmetry the average Free energy of the system (i.e. including the chemical potential term in the energy) may be evaluated using To fix the parameters a 0 and a 1 we propose a zero temperature scheme based on minimizing the free energy. In particular we are going to use ∆(a 0 , a 1 ) = 0, min to fix a 0 and a 1 . This is a constrained minimization problem and may be studied with standard methods in several ways. We can for example first fix a 1 and then try to solve the equation ∆(a 0 , a 1 ) = 0 for a 0 . This may in general have more than one solution so it is crucial in this scheme to always check the number of roots of ∆(a 0 ). In addition the parametern ↓ will be determined self-consistently.

V. NUMERICAL RESULTS FOR THE HALF FILLED CASE
In this section we are going to report some numerical results for the half filled case for a square lattice 100×100 at T = 0. Throughout we will measure energies in units of t, which amounts to setting t = 1. Half filling is obtained by taking µ = U/2. In Sec. VI below an analytical treatment of the half-filled case will be presented as well.
In particular we are going to report the results obtained for two possible set of parameter a 0 , a 1 which satisfy the constraint Eq. (31): the a 0 = 0, a 1 = 1 case and the a 0 , a 1 obtained by the variational scheme presented in Sec. IV A. From Eq. (31) it is possible to notice that for for a 0 = 0 we have ∆(0, a 1 ) = 0 independently on the value of a 1 and this is the only possible root, as can be seen in Fig. 1 (here we report ∆(a 0 ) only for a particular value of a 1 but the situation is the same for other values of a 1 ). To carry out the Free energy minimization carefully, it is important to have a sketch of the Free energy landscape as a function a 1 , since we will put a 0 = 0. A representative curve can be seen in Fig. 2, and we notice that F (a 1 ) posses a global minima for negative values of a 1 . In particular after carry out the constrained minimization numerically we found that the minimum of the free energy is reached for a 1 = −3, a 0 = 0, independently on the coupling strength U . At this point we are going to compare the avarage energy and double occupancy obtained for the two choices of the decoupling parameter a 1 = 1, a 0 = 0 and a 1 = −3, a 0 = 0 against the benchmark results gathered from Le Blanc et al., 16 reported respectively in Tab. I and Tab II. From Tab. I it is possible to notice both decouplings a 1 = 1 and a 1 = −3 predicts energies that may be lower than the exact methods. Consequently this scheme is not variational in the usual sense. This is expected since we are approximating the Green's function. In fact, despite that we know analytically the term neglected in the Green's function Eq. (29), the approximate Green's function properties are determined self-consistently within the truncated theory which can be different from the original one. From Tab. II we can notice that in the case a 1 = 1 the double occupancy drop to zero for U ≥ 8. In the other case a 1 = −3 double occupancy is predicted to be of the order (t/U ) 2 , which agree at least in order of magnitude with the benchmark results.
It is important to highlight that despite the crudeness of the two-pole approximation, this scheme independently on the value of parameter a 1 is capable to capture the effect of the correlation predicting a double occupancy that is significantly reduced from the mean field value n ↑ n ↓ = 1/4. The two possible choice of parameters a 1 =, a 0 = 0 and a 1 = −3, a 0 = 0 predict big differences at the level of predicted observables however. In particular for the case a 0 = 0, a 1 = 1 the truncated theory posses two energy bands shifted rigidly by U (i.e., independently of the momentum), as may be seen in Eq. (29) and in Fig. 3. With this choice of parameters the occupations of all the k-points in the first Brillouin zone are half occupied for U > 8t and for U < 8t we have a formation of fully occupied region around the Γ point surrounded by a region of half occupied k-points, and an empty region close to M point. The half-filled region in between shrinks as the interaction strength is decreased as can be seen in Fig. 4. We can also notice that in the limit U → 0, we recover the diamond-shaped Fermi surface for free fermions on the square lattice at half-filling.
To capture the metallic or insulting behavior of the solution one should in principle evaluate the conductivity or the charge-charge correlation function, which is in principle unaccessible with the operators used here. However we can have an indication on the metallic or in- sulating behavior of the system by analyzing density of states. Let us first consider the case a 1 = 1. In this case we can see in Fig. 5 that for U > 8t the density of states posses an hard gap and there is a formation of separated lower and upper Hubbard bands, which is a signature of an insulating phase. For U < 8t the lower and upper Hubbard bands overlap giving rise to a gapless density of state which is an indication of a metallic phase. Consequently for the choice of parameter a 1 = 1, a 0 = 0, we can notice that U = 8 represent a critical value of the interaction above which the system is in an insulating state and below which the system is in a metallic state. A radically different behavior is predicted by the choice of parameters a 1 = −3, a 0 = 0. In this case the system posseses two bands that repel with increasing interaction strength and there is always a small gap between the two of these parameters the occupation in first Brillouin zone is characterized by the presence of an almost fully occupied region around the Γ point which changes continuously to a low but non-zero occupation at the corner of the first Brillouin zone (the M point). As the interaction is decreased the almost fully occupied region around the Γ becomes increasingly occupied and the corner of the Brillouin zone get increasingly depleted. From the figures it looks like a Fermi surface is formed at a U = 0.1 along the high symmetry vector M Γ as it is possible to notice in Fig. 7. There is however a small gap that is not seen on this scale, this becomes clear in the analytic treatment in Sec. VI. In the limit of U → 0 also in this case we recover the diamond-shaped Fermi surface for a free electron gas on square lattice. As we did for the case a 1 = 1, a 0 = 0 we can also study the density of state in order to get an indication on the phase of the system. In this case there is always a gap between the upper and lower Hubbard band, but the gap is very small for small U as can be seen in Fig. 8 characterize the possible phases of the system for this choice of parameters it is going to be beneficial a study of the F (U ). In fact, the solution with a 1 = −3 always has a lower expectation value of the Free energy than the solution at a 1 = 1. Moreover, since a 1 = −3 is insulating, our scheme indicates that the insulator is stable at half-filling.

A. Relation to the two-pole approximation of Avella and collaborators
We can relate our approach to that of Avella et al by comparing the associated E-matrices. 15 The relation between their parameters (∆ and p) and ours (a 0 and a 1 ) are given by The issue of the determination of these parameters is discussed at length in Ref. 15. We note that they choose to determine ∆ self-consistently from the Green's function, and fix p so that Pauli principle is satisfied. This is different from our procedure where a 0 (and therefore ∆) is determined so that the Pauli principle is satisfied, in the next step we fix a 1 to minimize expectation value of the Free energy. Comparing the results we also have two classes of solutions, but the parameters obtained are not identical.

VI. ANALYTICAL RESULTS
Since the two-pole approximation involves 2 × 2 matrices everything may be evaluated exactly. A straightforward calculation gives (dropping k indexes on k and λ 2k and other parameters for brevity) where the poles are located at (37) The other parameters that are related to the weight of the poles are . (38c) Using this we may calculate many quantities of interest, such as the density of spin-up electrons the Pauli principle constraint (∆ = 0) and average double occupancy as well as the average kinetic energy (of two spin species)

A. Simplifying assumptions -insulator
It is possible to find the solution with a 1 = −3 obtained in the numerical study above analytically. In this subsection we present this solution is some detail since it provides an interesting zeroth order approximate Green's function at half-filling.
Let us assume that U is sufficiently large and chemical potential sufficiently small so that n −k = 1 and n +k = 0 for all k. We must then have Thenn ↑ =n ↓ = 1/2 solves Eq. (39). With this choice and therefore any λ 2 = a 1 will satisfy the Pauli principle constraint. The other parameters then become Using this me may write down expressions for average double occupancy and average kinetic energy Minimizing Ĥ 0 + U D we find a minimum at a 1 = −3, with the energy being The gain in energy due to hopping is increased with respect to more conventional approaches, such as antiferromagnetic mean field. Let us also note that the band structure for this solution is in the large-U limit we therefore get giving us two Hubbard bands with the full bare noninteracting bandwidth. Note however that the sign of the kinetic term is opposite to what if would be in the non-interacting case. The solution a 0 = 0, a 1 = 1 in the same region has D = 0 and Ĥ 0 = 0 so is always higher in energy than a 0 = 0, a 1 = −3. This agrees with our numerical findings.

VII. HOLE DOPED CASE IN THE STRONG COUPLING REGIME
In this section we are going to apply our scheme in the hole doped case for an interaction strength larger than the bandwidth namely U = 12. From the Free energy plots in Fig. 9 it is clear that upon hole doping (decreasing the chemical potential) the minima around a 1 = 1 is pushed down in energy with respect to the one near a 1 = −3, until it becomes the global one below a critical value near µ = 1.4. From an analysis of the DOS in Fig. 10, it is possible to notice that the solution around a 1 = −3 is characterized by the presence of an hard gap and the system is predicted to be insulating up to µ = 1.4. On the other hand the solution around a 1 = 1 is characterized by a smaller gap and when µ < 1.4, becomes the global minima of the Free energy. The DOS found in Fig. 11 suggests the formation of a metallic phase and it is possible to notice a spectral weight transfer from high energy states to the low energy ones.
Another interesting feature of Fig. 9 for µ ≤ 1.4 is that the Free energy as a function of a 1 features a discontin-  In this case the lower band get attracted to the upper one, pushing it above the chemical potential for certain momenta. This results in the formation of unoccupied kpoints in the first Brillouin and two Fermi surfaces that may be seen in Fig. 12. This may be viewed as a Lifshitz transition. 17,18 This solution is however not energetically favorable, and is not likely stabilized without additional interactions.

VIII. HOLE DOPED CASE IN THE INTERMEDIATE COUPLING REGIME
In this section we are going to apply our scheme to the hole doped case for an interaction comparable to the bandwidth, namely U = 4. The Free energy plots in Fig. 13 indicate that the situation is more involved in this case compared to the one obtained in the strong coupling limit of Sec. VII. There appears three local minima: one in the region a 1 ∈ [−4, −3], one in the region a 1 ∈ [0, 1], and one in the region a 1 ∈ [3,4]. In our discussion below we will call these minima m 1 , m 2 , and m 3 . For µ > 0.3 m 1 is the global minimum. Upon hole doping we can see that the local minimum m 3 is pushed down in energy and the minimum m 2 gets formed. For µ < 0.3 the minimum in m 3 becomes the global one until for µ < −0.3 the minimum in m 2 becomes the global minimum.
The character of the solution m 1 may elucidated from Fig. 14. It is characterized by an insulating gap and predicts the system to be half filled for µ ∈ [0. 3,2]. This solution is also characterized by the absence of a Fermi surface, and should be viewed as being in the same phase  as the corresponding half-filled solution studied above with a 1 = −3, and also the corresponding strong coupling solution. The solution m 3 may be characterized by studying Fig. 15, it is gapless which suggests a metallic state. There is moreover two sharp Fermi surfaces with discontinuities in the occupation numbers and a partially depleted "ring" around the Γ-point is formed. This interesting state is not present in the strong interaction solution.
The solution m 2 is also characterized by the absence of a gap as can be seen in Fig. 16. There is one sharp Fermi surface and consequently, in contrast to the solution m 3 , there is no formation of a depleted ring around the Γpoint. As in the strong coupling case above there exists discontinuities in some curves in Fig. 13. In particular for µ = 0.3 there is a discontinuity in F (a 1 ) around a 1 = 1.8 that is barely visible in the figure. The origin of this is a Lifshitz type transition where the Fermi surface change topology passing from a connected to a non-connected one as can be seen in Fig. 17. The discontinuity near a 1 ≈ −3.2 at µ = 0.3 is of the same type as the considered above, see Fig. 12.

IX. CONCLUSIONS AND OUTLOOK
In the context of the Green's function equation of motion method, we disclose the dependency between the algebra of the operators and their evolution, stressing that the Hermiticity of the E-matrix is a fundamental relation that all physical theories must satisfy. We also realized that for an arbitrary truncation the Hermiticity of the E-matrix is generally violated which leads to unphysical approximation for the Green's function.
To overcome this type of problem a novel truncation scheme for the equations of motion based on a partial orthogonalization was developed, in the context of the hierarchy of the operators. The main outcome of this procedure is an approximation for the fermionic Green's function, which can in principle be extended to an arbi-trary number of poles.
We applied this truncation scheme to a two-pole approximation for the Hubbard Model showing that the Hubbard-I and Mancini results can be obtained as a particular choices of a much wider range of decoupling possibilities. We introduced a variational procedure to determine the partial orthogonalization parameter(s). By employing it we analyzed a set of possible solutions for the two-pole approximation for the Hubbard model and we show that independently of the choice of the orthogonalization parameter both the atomic limit and the noninteracting limit are obtained as special cases for the halffilled case. Furthermore the solutions obtained, suggests the presence of a Mott metal-insulator transition both in the large coupling limit and in the intermediate one.
In the latter case we also find the presence of three competing solutions: one with an insulating character and two with metallic ones, characterized by different occupations in the first Brillouin zone and different number of Fermi surfaces. We want to stress that the variational procedure proposed to fix the parameters in this paper is not the only option available and in principle whatever decoupling parameters which satisfy the algebra constraint should be considered valid. Despite that, this method allows a transparent way to determine the part of the Green's function that is neglected from the original theory and constrain its total spectral weight. This enables further refinements of the approximate Green's function, where the effect of the neglected part can be incorporated in the theory using an adequate form of the self-energy.
In the end it is important to recall that this scheme can be applied both in the study of fermionic and bosonic systems. Various application of this novel decoupling scheme also in case of broken symmetries are planned for future works.

ACKNOWLEDGMENTS
Funding from the Knut and Alice Wallenberg Foundation and the Swedish research council Vetenskapsrådet is gratefully acknowledged.