{'fraction_non_alphanumeric': 0.03394770443568534, 'fraction_numerical': 0.011797035474379615, 'mean_word_length': 4.299646954986761, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Self-Consistent Random Phase Approximation - Application to the Hubbard Model for finite number of sites

Within the 1D Hubbard model linear closed chains with various numbers of sites are considered in Self Consistent Random Phase Approximation (SCRPA). Excellent results with a minimal numerical effort are obtained for 2+4n sites cases, confirming earlier results with this theory for other models. However, the 4n sites cases need further considerations. SCRPA solves the two sites problem exactly. It therefore contains the two electrons and high density Fermi gas limits correctly.

laws (Ward identies), Goldstone theorem, and restoration of spontaneously broken symmetries. Though there exist respectable general theories (see e.g. [7,8]), any practical attempt to go beyond this basic HF-RPA scheme conserving these properties turned out to be technically extremely demanding and no well accepted general and practical extension has emerged so far. Nevertheless, standard RPA has also quite serious shortcomings and it is desirable to overcome them. One of the most prominent is its violation of the Pauli principle, often paraphrased as the "quasiboson approximation". It is most critical for only moderately collective modes or when the self interaction of the gas of quantum fluctuations becomes important as in ultra small finite quantum systems. Since a couple of years two of the present authors and collaborators have been working on a non linear extension of RPA [9] which has shown surprisingly accurate results in a number of non-trivial models [10]. It is called Self -Consistent -RPA (SCRPA) and can be obtained from minimising an energy weighted sum rule. Therefore s-RPA which is perturbative in the sense that it sums a certain class of diagrams (the bubbles) is upgraded in SCRPA to a nonperturbative variational theory though it is in general not of the Raleigh-Ritz type. A strong bonus of this extension of s-RPA is that it generally preserves its positive features as conservation laws and restoration of symmetries as well as numerical tractability, since it leads to equations of the Schrödinger type [11]. In this paper we want to apply this theory to the Hubbard model for strongly correlated electrons. Because of its necessarily increased numerical complexity over s-RPA, we first want to consider finite clusters in reduced dimensions. Before going into the details, let us very briefly repeat the main ideas of SCRPA.
One way of presentation is to outline its strong analogy with the Hartree-Fock-Bgoliubov (HFB) approach to interacting boson fields b † , b. The HFB canonical transformation reads The amplitudes u, v can be determined [12] from minimising the following mean energy (energy weighted sum rule) where H is the usual many body Hamiltonian with two body interactions and the groundstate |0 is supposed to be the vacuum to the quasiboson operators q ν , i.e.
With this scheme and the usual orthonormalisation conditions for the amplitudes u, v, which allows the inversion of (1), one derives standard HFB theory [6] with no need to construct |0 explicitly. Of course, in this way the fact that HFB is a Raleigh-Ritz variational theory is not manifest but the scheme has the advantage to be physically transparent and to lead to the final equations with a minimum of mathematical effort. For SCRPA we follow exactly the same route. We replace in (1) the ideal boson operators by fermion pair operators of the particle-hole (ph) type and form an ansatz for a general transformation of ph -Fermion pairs with |ν = Q † ν |0 an excited state of the spectrum. In analogy with (2) we minimise a mean excitation energy with |0 , in analogy with (3), the vacuum to the operators Q ν , i.e.
Q ν |0 = 0 (6) and arrive at equations of the usual RPA type [6] with Here we supposed to work in a single particle basis which diagonalises the density matrix (natural orbits) : 0|a † k a k ′ |0 ≡ n k δ kk ′ (9) and therefore the n k 's are the occupation numbers. For H with a two body interaction, (8) only contains correlation functions of the a † a and a † a a † a types and, since (6) admits the usual RPA orthonormalisation relations for the amplitudes X , Y [6], the relation (4) can be inverted and with (6) the correlation functions in (8) be expressed by X , Y. However, to be complete, occupation numbers n k = 0|a † k a k |0 and two body correlation functions with other index combinations than two times particle and two times hole need extra considerations. That will be done in the main text. This is, in short, the SCRPA scheme which, as HFB theory, is obviously non linear, since the elements A and B in (7) become functionals of the X and Y amplitudes. We want to point out that no bosonisation of Fermion pairs is operated at any stage of the theory.
We want to apply this scheme to the Hubbard model of strongly correlated electrons which is one of the most wide spread models to investigate strong electron correlations and high T c superconductivity. Its Hamiltonian is given by where c † iσ , c jσ are the electron creation and destruction operators at site 'i' and then iσ = c + iσ c iσ are the number operators for electrons at site 'i' with spin projection σ. As usual t is the nearest neighbour hopping integral and U the on site coulomb matrix element. In this exploratory work, we will limit ourselves to the simplest cases possible, i.e. we will consider closed chains in one dimension with increasing number of sites at half filling, starting with the 2 -sites problem. It will turn out that the next case of 4 -sites is a configuration with degeneracies which cause problems in SCRPA, as do all 4n (n = 1, 2, 3, . . .) configurations in 1D. We therfore will postpone the treatment of these cases to future work and directly jump to the case of 6 -sites and only shortly outline at the end why the 4 -sites case is unfavorable and how the problem can eventually be cured. In this work we will stop with the 6 -sites case considering it as sufficiently general to be able to extrapolate to the more electron case. In this way one may hope to approach the thermodynamic limit in increasing the number of sites as much as possible. Let us mention that an earlier attempt to solve SCRPA in 1D in the thermodynamic limit in a strongly simplified version of SCRPA, the so-called renormalised RPA (r-RPA), produced interesting results [13].
In detail our paper is organised as follows : in section (II) we present the two sites case with its exact solution. In section (III) we outline the 6 -sites case with a detailed discussion of the results and in section (IV) we present the difficulties encountered in the 4 -sites case and how, eventually, one can overcome them. Finally in section (V) we give our conclusions together with some perspectives of this work.

II. THE TWO SITES PROBLEM
In this section we will apply the general formalism of SCRPA outlined in the introduction to the two sites problem at half filling, i.e. two electrons with periodic boundary conditions. This case may seem trivial, the fact, however, is that such popular many body approximations as standard RPA (s-RPA), GW [14], Gutzwiller wave function [15], the TPSC approch by Vilk and Tremblay [16], etc.. do not yield very convincing results in this study case, whereas it has recently been shown that SCRPA solves two body problems exactly [10,11,17]. We again will briefly demonstrate this here for the two sites problem.
First we will transform (10) into momentum space. With the usual transformation to plane waves c j,σ = 1 √ N k ak ,σ e −ikx j this leads to the standard expression for a zero range two body interaction : wherenk ,σ = a † k,σ ak ,σ is the occupation number operator of the mode (k, σ) and the single particle energies are cos (k d ) with the lattice spacing set to unity.
For our further considerations it is convenient to transform (11) to HF quasiparticles operators via (we switch to 1D) where h and p are momenta below and above the Fermi momentum, respectively, so that b k,σ |HF = 0 for all k where |HF is the Hartree-Fock groundstate in the plane wave basis. For the two sites problem with periodic boundary conditions we then write, after normal ordering, the Hamiltonian (11) in the following way with and where we introduced the abreviation "1" and "2" for the two momenta k 1 = 0 and k 2 = −π of the system, respetively. The HF groundstate is |HF = b 1,↑ b 1,↓ |vac and the corresponding energy is given by The RPA excitation operator corresponding to (4) can, because of rotational invariance in spin space, be separeted according to spin singlet (S = 0, charge) and spin triplet (S = 1) excitations. The latter still can be divided into spin longitudinal (S = 1, m s = 0) and spin transverse (S = 1, m s = ±1) excitations. Let us first consider the charge and spin longitudinal sectors. For later convenience we will not separate them and write for the corresponding RPA operator where K ± σ = J ± σ / 1 − M σ , M σ =ñ 1σ +ñ 2σ , and the mean values . . . are always taken with respect to the RPA vacuum Because of the orthonormality relations one can invert (18) to obtain With (21) the expectation value of (26) gives For the calculation of the correlation functions which appear on the right hand side of (27) one commutes the destructors Q ν to the right and uses (6), yielding again correlation functions M σ M σ ′ . One then obtains a closed linear system of equations for the latters. Details are given in Appendix (A). The SCRPA matrix elements can be expressed in the following way With our previous relations (23), (24) and (27) we can entirely express the elements of (28) and (29) by the RPA -amplitudes and therefore we have a completely closed system of equation for the amplitudes X , Y. With the orthonormality relations (20) we further more have and, therefore, the SCRPA equation can be written in the following form (31) The system (31) has the two positive The SCRPA equation (31) can be solved numerically by iteration leading, as expected, to the exact result. This latter fact can also be seen analytically in noticing that by symmetry Therefore the 4 × 4 equation (31) decouples into two 2 × 2 equations corresponding to charge (ch) and spin (sp).
Then we see that the exact groundstate wave function which contains only up to 2p − 2h excitations is the exact vacuum to the RPA operators, i.e Q ch(sp) |RP A = 0 under the condition that We therfore can express the SCRPA equations by the single parameter φ and obtain the solution analytically (up to the solution a nonlinear equation for φ). The solution agrees for all quantities with the exact result. For example the groundstate energy is given by This expression can either be derived directly from H using (33) and (34) or one uses a generalisation of the standard RPA expression for the groundstate energy [6]: (36) It is straightforward to verify that expressions (35) and (36) are identical.
The standard RPA expression are recovered from (31) in replacing in all expectation values the RPA groundstate by the uncorrelated HF determinant. In Fig.1 we compare the standard RPA with SCRPA and exact results for the excitation energies and in Fig.2 the corresponding groundstate energies together with the HFvalues are shown. From these figures one should especially appreciate the long way SCRPA has gone from s-RPA to recover the exact result. For instance it is clearly seen that the instability of s-RPA at U = 2 is, as expected for such a small system, an artefact and is completely washed out by the self consistent treatment of quantum fluctuation contained in the SCRPA approach. Without explicit demonstration let us also mention that SCRPA in the spin transverse channel with as well as in the particle-particle channel with also gives the exact solution for the two sites problem. How the pp-SCRPA works can be seen in ref [10] where for the pairing problem the two particle problem is also solved exactly. The fact that SCRPA solves the two sites problem exactly is non trivial, since other well known many body approaches [14,15,16], as already mentioned, so far failed to obtain this limit correctly.

III. THE SIX SITES PROBLEM
After this positive experience with the two sites problem we next will consider the 1-dimensional 6 -sites case, as for the 4 -sites case problems appear needing particular considerations to be outlined in section IV. We again consider the plane wave transformation explained in section II with the corresponding Hamiltonian in momentum space (11). In the first Brillouin zone −π ≤ k < π we have for N = 6 the following wave numbers With the HF transformation such that b k,σ |HF = 0 for all k, we can write the hamiltonian in the following way (normal order with respect to b † , b) with the following definition of operators and The level scheme is shown in Fig.3. The hole states are ε π π/3 π/3 π/3 π/3 labeled h = {1, 2, 3} and the particle states p = {4, 5, 6}. The HF groundstate is We see that the Hamiltonian for 6 -sites has largely the same structure as the one for 2 -sites. It is only augmented by H |q|= π 3 + H |q|= 2π 3 which contains the Soperators on which we will comment below.
There are three different absolute values of momentum transfers as shown in Table I. Since the momentum transfer |q| is a good quantum number, the RPA equations are block diagonal and can be written down for each |q| -value separately. For exemple for |q| = π 3 we have the following RPA operator for charge and longitudinal spin excitations where and M pσ,hσ =ñ p,σ +ñ h,σ .
We write this RPA operator in short hand notation as again with the properties The matrix elements in the SCRPA equation Since the SCRPA equations have the same mathematical structure as standard RPA, one also has equivalent orthonormality relations etc.. , in analogy to eqs. (20) of the 2 -sites case. This allows to invert (47) and to calculate the expectation values which will appear in (49a) and (49b) in complete analogy to (23). The missing expectation values M i can be expressed by the X , Y amplitudes in observing that J ± i , and J 0 i = 1 2 (M i −1) form, as in the 2 -sites case, an SU 2 Lie algebra for spin-1 2 particles. Using the Casimir relation one again obtains M i = 2 J + i J − i and thus We also will need expectation values of (50)). Those can again be calculated following the same procedure as outlined in (27) and Appendix (A).
In order to solve the SCRPA equations we now practically have prepared all we need. Nontheless, at this point we have to discuss a limitation of our RPA ansatz (44) which is not absolutely necessary but which turned out to be convenient for numerical reasons. The fact is that our RPA ansatz is restricted to ph and hp configurations, as this is also the case in standard RPA. In the latter case this is a strict consequence of the use of HF occupation numbers n 0 p , n 0 h with values zero or one, respectively. In the SCRPA case with a correlated groundstate the occupation numbers are different from zero and one and a priori there is no formal reason not to include into the RPA operator also pp and hh configuration of the Such terms are usually called scattering or anomalous terms [19]. With rounded occupation numbers the SCRPA equations (at T = 0) are formally and mathematically equivalent to standard RPA equations at finite temperature where also pp and hh components are to be included, in principle [18]. The inclusion of those scattering terms [18,19] (the S-terms in (39)) usually is of little quantitative consequence [11], entails, however, the important formal property that, as for standard RPA, the energy weighted sum rule is fullfilled exactly [11,19]. Inspite of this desirable feature, we had to refrain from the inclusion of the scattering configurations in this work because the factors 1 − M i by which the SCRPA matrix is divided (see eqs. (49a) and (49b)), can become very small in these cases and this perturbed the convergence process of the iterative solution of the SCRPA equations. Though we do not exclude that a more adequate numerical procedure could be found to stabilise the iteration cycle, we decided to postpone such an investigation, because, as already mentioned and as will be shown later, the influence of the scattering terms is, as found already in other studies [11], very small. We will shortly come back to this discussion when presenting the results for the energy weighted sum rule below. As a consequence and for consistency we then also will have to disregard the S-terms of the Hamiltonian (remember that also in standard RPA these terms do not contribute). Under these conditions we then obtain a completely closed system of SCRPA equations. For completeness we give some examples of SCRPA matrix elements which correspond to the ansatz (44) for |q| = π . . .
The other matrix elements can be elaborated along the same lines. Of course in the approximation where the expectation values in (51a) and (51b) are evaluated with the HF groundstate the usual matrix elements of standard RPA are recovered. We should also mention that in expressions (51a) and (51b) expectation values as for example J − 1↑,4↑ J + 4↓,1↓ which involve momentum transfers other than the one under consideration (|q 3 | = π 3 in the specific example) must be discarded. That this implicit channel coupling cannot be taken into account without deteriorating the quality of the SCRPA solutions is an empirical law which has been established quite sometime ago [20]. It is part of the decoupling scheme and it is intuitively understandable that, since each channel is summing specific correlations, one can not mix the channels implicitly without perturbing the balance of the minimisation procedure which is done channel by channel. It can also be noticed that, neglecting the S-terms in H, the channel coupling disappears.
We here give for the transfer |q| = π 3 the totality of the elements of matrix SCRPA A and B just as it is was used in the numerical calculation. For others transfers there will be analogous expressions. Indeed with the following abbreviations the elements of matrices A and B are given by Let us add that the matrices A and B are symmetric and that the expectation values . . . in (52a) and (52b) can be expressed in an analogous way as the expectation values (23) and (27) by the amplitudes X , Y. The structure of the self consistent matrix elements (52a) and (52b) is also quite transparent : the bare interaction which survives in the limit of standard RPA is renormalised, i.e. screened, by two body correlation functions which are calculated self consistently. The general structure of the scheme is in a way similar to the one proposed by Tremblay and coworkers [16], however, the details of the expressions and the spirit of derivation are different. One can also interpret our theory as a mean field theory of quantum fluctuations as this was done in refs. [9]  Let us now come to the presentation of the results. In Fig. (4, 5 and 6) we display the excitation energies in the three channels |q| = π, 2π 3 and π 3 as a fonction of U/t. The exact values are given by the continuous lines, the SCRPA ones by crosses and the ones corresponding to standard RPA by the broken lines. We see that in all three cases SCRPA results are excellent and strongly improve over standard RPA. As expected, this is particularly important at the phase transition points where the lowest root of standard RPA goes to zero, indicating the onset of a staggered magnetisation on the mean field level. It is particularly interesting that SCRPA allows to go beyond the mean field instability point. However, contrary to the two sites case where SCRPA, in the plane wave basis, solved the model for all values of U , here at some values U > U cr the system "feels" the phase transition and SCRPA stops to converge and also deteriorates in quality. Up to these values of U SCRPA shows very  In Fig. 7 we show the groundstate energy (see eq. (36)) as a function of U . In addition to exact, SCRPA, and s-RPA values we also show the HF energy. Again we see that SCRPA is in excellent agreement with the exact solution. Standard RPA is also good for low values of U but strongly deteriorates close to the lowest phase transition point which occurs in the |q| = π channel at U = 12 t 5 . The HF energies, on the contrary, deviate quite strongly from the exact values. The reader certainly has remarked that our RPA ansatz (44) has so far not separated charge and spin excitations. In the 2 -sites problem this was automatically and exactly the case. However, here, since we did not consider the S -operators neither in the Hamiltonian nor in the RPA operator, spin symmetry is violated. On the other hand this permits to evaluated the importance of the S -operators. Normally the eigenvectors of the RPA matrix should be such that for charge (ch) excitations the operators J + ph↑ + J + ph↓ and J − ph↑ + J − ph↓ can be factored whereas for spin (sp) excitations the combinations J + ph↑ − J + ph↓ and J − ph↑ − J − ph↓ hold. Because of our violation of spin symmetry this factorisation is not exact. To have a measure of this violation we plot in Fig. 8 the ratio For exact spin symmetry r should be zero. From Fig. 8 we see that the violation is on the level of a fraction of one percent. This, therefore justifies, a posteriori having neglected the scattering terms (S -terms) in the Hamiltonian and RPA operator. A further indication that Sterms are not important comes from the energy weighted sum rule. We know that the sum rule including the S -terms is fullfilled in SCRPA [13,19]. However, neglecting them gives a slight violation. Considering the exact relation with we trace in Fig. 9 the ratio ξ = R−L R . Again we see that the violation is on the level of a fraction of one percent, confirming the very small influence of the scattering terms.
A further quantity which crucially tests the ground state correlations are the occupation numbers. We have no direct access to them, however, we will use the socalled Catara approximation for their evaluation [21] : We show these quantities in Figs. 10 and 11 in comparison with the exact values and the ones of standard RPA. We again see the excellent performance of SCRPA. Concluding this section we can say that the expectation we had from the 2 -sites case, with its exact solution, have very satisfactorily also been fullfilled in the 6 -sites case. However, in spite of the very good performance of SCRPA, there is the limitation that SCRPA, in the symmetry conserving basis of plane waves used here, can not be employed in the strong U limit. One also may wonder how the extension to cases with sites number 2 + 4n with n > 1 works. For such cases it does not make sense any more to elaborate the Hamiltonian in its detailed form as given in eq.(40). This explicit expression was only given to make clear the detailed internal structure of the approach for a definite example. In the general case with many sites one would just take the form (11) of the Hamiltonian, calculate the double commutators as needed in (8) and then express the resulting correlation functions by the X -and Y-amplitudes. That such a program is feasable in terms of analytic work and numerical execution was demonstrated in our earlier work on the multilevel pairing model [10] where cases up to hundred levels were treated. However, this number was not considered of an upper limit. Though the present For each approximation, s-RPA and SCRPA and exact solution, the occupation numbers are represented like k = 0, π 3 , − π 3 . Let us notice that the modes k = π 3 and k = − π 3 are degenerate.
model is slightly more complicated, we think that a generalisation to the case of many sites is perfectly possible. It needs, however, some investment which is planned for the future. This also concerns the D = 2 case. Another question to ask is whether the degradation of the SCRPA results going from the N = 2 to the N = 6 case does not go on considering N = 10, 14, etc? One again may cite the experience with the multilevel pairing model [10] where also the N = 2 case turned out to be exact in SCRPA but not the other cases. However, all N > 2 cases showed more or less the same degrees of accuracy: excellent results of SCRPA up to the phase transtion point and deterioration beyond. Since this behavior has also been found in simpler models [12] we think that this is a generic feature of SCRPA and that this behavior will also translate to the case of the present model.
Another problem for further work is how to continue the present theory into the strong coupling regime. Of course, there exists the possiblilty to perform SCRPA in the symmetry broken basis, but details and how to match with the summetry unbroken phase must still be worked out. Also the inclusion of higher order operators, as will shortly be discussed in the next section, may be an intersting direction in this respect.

IV. FOUR SITES PROBLEM
A. The symmetry unbroken case The problem of the 4 -sites case is easily located in regarding the level scheme of Fig. 12 (see also ref [22] deeling with the attractive Hubbard model in 1D). We see that the Fermi energy coincides with the second level which is half filled. The uncorrelated groundstate is therefore degenerate and excitations with momentum transfer |q| = π cost no energy. On the other hand for excitations with |q| = π 2 there is no problem. The corresponding RPA operator is given by In Fig. 13 we show the results of s-RPA and SCRPA, together with the exact solution. We see that the lower excitation is still very well reproduced by SCRPA, whereas for the second excited state SCRPA only reduces the difference of s-RPA to exact by half. The real problem shows up for the transfer |q| = π. The corresponding operator is The standard RPA produces a doubly degenerate zero  mode independent of U as seen on Fig. 14. As compared with the exact solution, we see that these two zero modes approximate two very low lying exact solutions. Unfortunately, because of these modes at low energy the SCRPA could not be stabilised. The only possibility consisted in excluding the components K ± 32,↑ and K ± 32,↓ in the RPA operator. Then the self consistency was achieved without problem and the result is shown in Fig. 14. The result of SCRPA is half way in between s-RPA and the exact solution. On the other hand, because of the omission of the two lower states, the groundstate energy can not correctly be calculated in SCRPA. Therefore, for the 4 -sites problem in the symmetry unbroken basis (plane waves), the SCRPA cannot fully account for the situation.

B. Symmetry broken basis
An analysis of the HF solution shows that, as soon as U = 0, the plane wave state becomes unstable and the system prefers a staggered magnetisation. The general HF transformation can be written as with u = cos(ϑ) and v = sin(ϑ) e i ϕ . The minimisation of the groundstate energy with shows that ϕ = 0 for any value of U and the angle ϑ is obtained from The occupation numbers are given by n 1,↑ = n 3,↑ = n 2,↓ = n 4,↓ = 1 2 1 + sin 2 (ϑ) and shown in Fig. 15 which illustrates the spontaneous symmetry breaking for any value of U . For U → ∞ we have a perfect antiferromagnet.
We can now perform a SCRPA calculation in the symmetry broken basis. The RPA operators are given by with σ = ± 1 2 . We also have two other excitation operators and In Figs. 16 and 17 we give the results. The most striking feature is that s-RPA and SCRPA are very close and that the error with respect to the exact solution does not become greater than 25% for any value of U . Though the improvement of SCRPA over s-RPA is very small in each channel, at the end in the groundstate energy this sums to a more substantial correction in the right direction for the ground state energy. This is shown in Fig. 17 as a function of atan( U t ). We see that HF, s-RPA and SCRPA become exact for U = 0 and U → ∞. In between SCRPA deviates e.g. by 8% from the exact result at U ≃ 6 (atan( U t ) ≃ 1.4) whereas this deviation is 20% for s-RPA.
Concluding this section on the 4 -sites case at half filling we can say that in the symmetry unbroken basis SCRPA is unable to account for some low lying excitations and therefore fails to reproduce the groundstate energy as well. In the symmetry broken basis SCRPA gives very little correction over s-RPA. However, the maximum error is not greater than 25% for all values of U for the excited states and the groundstate energy in SCRPA whereas this is 30% for standard RPA. This may be an interesting result in view of the importance of the so-called 'plaquettes' (see e.g. ref [23]) in high T c superconductivity. Nevertheless, even though one plaquette (4 -sites) Groundstate energies in HF, standard RPA, SCRPA, and exact solution as a function of atan( U t ) for four sites with spin projection ms = 0 in the symmetry broken basis. may reasonably be described, the present approach can not account for the situation of many plaquettes in interaction which is the real situation in 2D. For the future it is therefore very interesting to develope an extension of the present SCRPA which not only gives the exact solution for the 2 -sites case but equally for the 4 -sites case. Such a generalisation is possible in including into the RPA operator in addition to the Fermion pair operators also quadruples of Fermion operators. This is a general principle and it has already been demonstrated to hold true in the case of the simpler Lipkin model [24]. One could call such an extension a second SCRPA in analogy to the well known standard second RPA which involves in addition to the ph configurations also 2p − 2h ones. In the case of many plaquettes this second SCRPA would then constitute a self consistent mean field theory for plaquettes.

V. DISCUSSION, CONCLUSIONS AND OUTLOOK
In this work a many body approach which has essentially been developed in the nuclear physics context in recent years [9] has been applied to the Hubbard model for finite number of sites. The theory is an extension of standard RPA, called Self Consistent RPA (SCRPA), which aims to correct its well known deficiencies as the quasiboson approximation with its ensuing violation of the Pauli principle and its perturbation theoretical aspect. Of course the appealing features of RPA, as for instance fullfillment of sum rules, restoration of broken symmetries, Goldstone theorem, numerical practicability and physical transparency should be kept as much as possible. That this is indeed the case with SCRPA has in the past been demonstrated with applications to several non trivial models [10] as for instance the many level pairing (Richardson) model [10] and the 3-level Lipkin model [11]. SCRPA can be derived by minimising an energy weighted sum rule and it is therefore a non perturbative variational approach though it is in general not of the Raleigh Ritz type. The resulting equations are a non linear version of the RPA type which can be interpreted as the mean field equations of interacting quantum fluctuations. Though the SCRPA equations are of the Schrödinger type, their non linearity non the less makes their numerical solution quite demanding. We therefore thought it indicated to begin with applications to the Hubbard model restricting them to low dimensional cases given by a finite number of sites where exact diagonalisation can easily be obtained. We then logically started out considering the two sites case (with periodic boundary conditions), increasing the number of sites by steps of two, i.e N = 2, 4, 6, . . .. To our satisfaction SCRPA solves the 2 -sites problem exactly for any value of U . This, as a matter of fact, did not come entirely as a surprise, since the same happened already with the pairing problem for two fermions [10] and indeed it can be shown that SCRPA solves a general two body problem exactly [17]. It is nontheless worth pointing out that other respectable many body theories fail in the two particle case, apart from the low U -limit.
In the four sites problem at half filling SCRPA failed. This, as in all 4n (n = 1, 2, 3, . . .) cases, presents the particular problem that the system is unstable with respect to the formation of staggered magnetisation for any finite value of U and this prevented the SCRPA solution to exist in the plane wave basis for particular values of the momentum transfer |q|. At the end of the paper we indicated that extending the present RPA ansatz of ph pairs to include quadruples of fermion operators can solve not only the two electrons but also the four electrons case exactly. This is particularly interesting in view of the fact that the 4 -sites case (plaquette) may be very important for the explanation of high T c superconductivity, in considering the many plaquette configurations in 2D [23]. In this work we jumped directly to the six sites problem which, as all 2 + 4n cases, causes no particular difficulties in SCRPA, even in the symmetry unbroken basis of plane waves. Of course, in the case of 6 -sites, SCRPA is not exact any more. However, it is shown that the results are still excellent for all quantities considered : excited states, groundstate, and occupation numbers. Contrary to the two sites case, the SCRPA solutions in the plane wave basis cannot be obtained for all values of U . Somewhere after the point where, as a function of U , the first mean field instability shows up, the SCRPA also starts to deteriorate and in fact does not converge any longer. Often the mean field critical value of U is by passed by 20% up to 50% in SCRPA, still staying excellent. However, to go into the strong U -limit we have to introduce the above mentioned quadruple fermion operators or perform a SCRPA calculation in the symmetry broken basis [12]. Such investigations shall be left for the future. We also gave arguments why we think that, going to the N > 6 cases, the precision we found for N = 6 will not deteriorate. We therfore think that our formalism will allow to find precise results for system sizes where an exact diagonalisation becomes prohibitive. Problems in 2D with closed shell configurations probably also can and shall be considered with the present formalism. Also, as shown in [10], the extension to finite temperatures is possible.
We also should mention that in this work we neglected the so-called scattering terms of the form a † p a p ′ or a † h a h ′ , that is fermion ph operators where either both indices are above or both below the Fermi level. In standard RPA those configurations automatically decouple from the ph and hp spaces. However, in SCRPA with its rounded occupation numbers, there is formally no reason not to include them. As a matter of fact, as shown in earlier work [11,19], to assure the fullfillment of the f-sum rule and the restoration of broken symmetries, these scattering terms must be taken into account. In the present case, as well as in earlier studies, the scattering terms seem to be almost linearly dependent with the ordinary ph and hp configurations. This fact induced difficulties with the iteration procedure, since they correspond to very small eigenvalues of the norm matrix. Though we do not exclude the possibility that this difficulty could be mastered with a more refined numerical algorithm, we finally refrained from persuing this effort, since we could show that the influence of the scattering terms on the results is only on the level of a fraction of percent and also the f-sum rule is only violated on this order.
In short we showed that SCRPA, as in previous models, performs excellently in the symmetry unbroken regime of the Hubbard model. However, the high U -limit and the 4n sites cases need further developements.
Finally, one can express the correlation function according to the amplitudes RPA, M i and of M i M j as .
with the average value N σNσ ′ = N σ + N σ ′ , which gives us Thus for our case, there is the relation