Coupled Self-Consistent RPA Equations for Even and Odd Particle Numbers. Tests with Solvable Models

TL;DR

This paper develops coupled self-consistent RPA equations for even and odd particle numbers, enabling self-consistent calculation of occupation probabilities, and demonstrates their effectiveness with solvable models like Lipkin and Hubbard.

Contribution

Introduces a novel coupled equation framework for even and odd particle numbers based on the same correlated vacuum, extending SCRPA methods.

Findings

Very good results in Lipkin model tests

Accurate occupation probabilities for odd particle numbers

Unified approach improves consistency between even and odd cases

Abstract

Coupled equations for even and odd particle number correlation functions are set up via the equation of motion method. For the even particle number case this leads to self-consistent RPA (SCRPA) equations already known from the literature. From the equations of the odd particle number case the single particle occupation probabilities are obtained in a self-consistent way. This is the essential new procedure of this work. Both, even and odd particle number cases are based on the same correlated vacuum and, thus, are coupled equations. Applications to the Lipkin model and the 1D Hubbard model give very good results.

Tables

Table 1. Table 1: The various momentum transfers in the 6 sites case.

$q=\pm \frac{2\pi }{3}$	$q=\pm \pi$	$q=\pm \frac{\pi }{3}$
$51\to q_{51}=+\frac{2\pi }{3}$	$61\to q_{61}=-\pi$	$42\to q_{42}=-\frac{\pi }{3}$
$63\to q_{63}=+\frac{2\pi }{3}$	$52\to q_{52}=+\pi$	$53\to q_{53}=+\frac{\pi }{3}$
$41\to q_{41}=-\frac{2\pi }{3}$	$43\to q_{43}=-\pi$
$62\to q_{62}=-\frac{2\pi }{3}$

Equations

$$q_{\mu }^{†}=\sum _h{x}_h^{\mu }{a}_h+\sum _{p{p}^{\prime }{h}^{\prime }}{U}_{p{p}^{\prime }{h}^{\prime }}^{\mu }{a}_{h^{\prime }}^{†}{a}_{p^{\prime }}{a}_p,$$

$$q_{\rho }^{†}=\sum _p{x}_p^{\rho }{a}_p^{†}+\sum _{p^{\prime }{h}^{\prime }h}{U}_{p^{\prime }{h}^{\prime }h}^{\rho }{a}_h^{†}{a}_{h^{\prime }}^{†}{a}_{p^{\prime }}.$$

$$q_{\mu }|Z⟩=0,$$

$$|Z⟩$$

$$\sum _h{x}_h^{\mu \ast }{Z}_{p{p}^{\prime }h{h}^{\prime }}$$

$${\lambda }_{\mu }$$

$$Q_{\nu }^{†}=\sum _{ph}{X}_{ph}^{\nu }{a}_p^{†}{a}_h-Y_{ph}^{\nu }{a}_h^{†}{a}_p.$$

$${\Omega }_{\nu }$$

$$Q|Z⟩=0.$$

$$\displaystyle\left(\begin{array}[]{cc}A&B\\ -B^{*}&-A^{*}\end{array}\right)\left(\begin{array}[]{c}X^{\nu}\\ Y^{\nu}\end{array}\right)=\Omega_{\nu}\left(\begin{array}[]{c}X^{\nu}\\ Y^{\nu}\end{array}\right)$$

$$\sum _{ph}\left(|X_{ph}^{\nu }{|}^2-|Y_{ph}^{\nu }{|}^2\right)=1$$

$$A_{ph,p^{\prime }{h}^{\prime }}$$

$$B_{ph,p^{\prime }{h}^{\prime }}$$

$$H=\sum _{k{k}^{\prime }}{t}_{k{k}^{\prime }}{a}_k^{†}{a}_{k^{\prime }}+\frac{1}{4}\sum _{klmn}{\bar{v}}_{klmn}{a}_k^{†}{a}_l^{†}{a}_n{a}_m.$$

$${\bar{v}}_{klmn}=⟨kl|v|mn⟩-⟨kl|v|nm⟩.$$

$$\sum _{h^{\prime }}{ϵ}_{h{h}^{\prime }}{x}_{h^{\prime }}^{\mu }+\sum _{p{p}^{\prime }{h}^{\prime }}{\mathcal{C}}_{h,p{p}^{\prime }{h}^{\prime }}{U}_{p{p}^{\prime }{h}^{\prime }}^{\mu }={\lambda }_{\mu }{x}_h^{\mu }$$

$$\sum _{h^{\prime }}{\mathcal{C}}_{p{p}^{\prime }h,h^{\prime }}^{\ast }{x}_{h^{\prime }}^{\mu }+\sum _{p_2{p}_1{h}_1}{\mathcal{D}}_{p{p}^{\prime }h,p_2{p}_1{h}_1}{U}_{p_2{p}_1{h}_1}^{\mu }={\lambda }_{\mu }{U}_{p{p}^{\prime }h}^{\mu }$$

$$\left(\begin{array}{cc}ϵ& \mathcal{C}\\ {\mathcal{C}}^{†}& \mathcal{D}\end{array}\right)\left(\begin{array}{c}{x}^{\mu }\\ U^{\mu }\end{array}\right)={\lambda }_{\mu }\left(\begin{array}{c}{x}^{\mu }\\ U^{\mu }\end{array}\right)$$

$${ϵ}_{h{h}^{\prime }}$$

$${\mathcal{C}}_{p{p}^{\prime }{h}^{\prime },h}^{\ast }$$

$${\mathcal{D}}_{p{p}^{\prime }{h}^{\prime },p_2{p}_1{h}_1}$$

$${\mathcal{N}}_{p{p}^{\prime }{h}^{\prime }}$$

$$a_p^{†}{a}_h=\sqrt{n_h-n_p}\sum _{\nu }{X}_{ph}^{\nu }{Q}_{\nu }^{+}+Y_{ph}^{\nu }{Q}_{\nu }$$

$$n_h=⟨a_h^{†}{a}_h⟩=\sum _{\mu }|⟨\{a_h^{†},q_{h,\mu }^{†}\}⟩{|}^2=\sum _{\mu }|x_h^{\mu }{|}^2$$

$$n_p=⟨a_p^{†}{a}_p⟩=\sum _{\rho }|x_p^{\rho }{|}^2$$

$${\mathcal{G}}_h(\omega )=\frac{1}{\omega -{ϵ}_h-{\mathcal{M}}_h(\omega )}$$

$${\lambda }_{\alpha }-{ϵ}_h-{\mathcal{M}}_h({\lambda }_{\alpha })=0$$

$${\mathcal{M}}_h=\sum _{p{p}^{\prime }{h}^{\prime },p_2{p}_1{h}_1}{\mathcal{C}}_{h,p{p}^{\prime }{h}^{\prime }}^{\ast }(\omega -\mathcal{D}{)}_{p{p}^{\prime }{h}^{\prime },p_2{p}_1{h}_1}^{-1}{\mathcal{C}}_{p_2{p}_1{h}_1,h}$$

$${\mathcal{G}}_h(\omega )=\sum _{\alpha }\frac{r_{\alpha }}{\omega -{\lambda }_{\alpha }}$$

$$r_{\alpha }=\frac{1}{1-{\mathcal{M}}^{\prime }{|}_{\omega ={\lambda }_{\alpha }}}$$

$$i{\mathcal{G}}_h(t-t^{\prime })$$

$$i{\mathcal{G}}_p(t-t^{\prime })$$

$$n_h$$

Full text

I Introduction

Developments of Many-Body approaches for strongly correlated systems is an active field of research. In the past we developed an RPA theory which goes beyond the standard one and which is based on a correlated ground state. To lowest order this leads to RPA equations where the single particle (s.p.) occupation numbers are not the uncorrelated (Hartre-Fock (HF)) ones but correlated ones which are obtained in a self-consistent way from the RPA solution. It is known as ’re-normalized’ RPA (r-RPA) DSD05 ; SchaferP . However, in general the RPA equations contain additionally vertex (e.g., screening) corrections which also can be obtained self-consistently from the RPA solution. The whole procedure has been dubbed Self-Consistent RPA (SCRPA) jemai13 and references in there. SCRPA can also be seen as a sub-product of an even more general approach which is the Time-Dependent Density Matrix (TDDM) theory based on a decoupling of the BBGKY hierarchy of one, two, etc. correlation functions SchuToh ; tddm-scrpa . In the past, there was always a certain debate how to include the single particle (s.p.) occupations into the SCRPA scheme. Since the latter are s.p. quantities and the RPA gives rise to two body correlations, it was not completely evident in which way to close the system of equations. However, already Rowe promoted the so-called particle number operator method to obtain self-consistent s.p. occupations Row68 . Later this was further elaborated by F. Catara Catara and it has become known as the ’Catara method’ since then. This method mostly works quite well jemai05 but also fails more or less in particular cases jemai13 . It is for this reason that in this work we elaborate a different scheme which seems to be more natural since it is also based on the equation of motion (eom) method and employs the so-called odd particle number RPA (o-RPA) recently proposed by one of the authors plus collaborators TohSchu13 . The latter can also be reformulated as a Dyson equation for the s.p. Green’s function (GF) with a self-energy obtained from the eom for the $2$particle-$1$hole (2p-1h) and 2h-1p GF. Since both, the SCRPA and o-RPA will be based on the same correlated vacuum, naturally even and odd particle number channels become coupled. We may coin this scheme eo-SCRPA. We will apply this approach to two exactly solvable model cases: the Lipkin model and the 1D Hubbard model. In both cases the results turn out to be promising.

The paper is organized as follows: in section II we present the general theory. In section III, applications to the Lipkin and Hubbard models are given. Section IV contains the conclusions and details of our procedures are presented in the Appendices.

II General Theory

As mentioned in the introduction, we will base our approach on the coupling of the even and odd particle number eom. The latter will be obtained from the following ansatz

[TABLE]

where the indices ’p, h’ refer to s.p. states ’above’ and ’below’ the Fermi surface, respectively, and $a^{\dagger}_{k},a_{k}$ are the fermion creation and annihilation operators. The equation for the s.p. basis in which the equations will be worked out will be given below. Our ’quasi-particle’ operator in (1) has the good quality that its destructor exactly kills the so-called Coupled-Cluster-Doubles (CCD) wave function:

[TABLE]

with

[TABLE]

where $|\rm HF\rangle$ is the Hartree-Fock (HF) Slater determinant and the amplitudes must full-fill the following relations

[TABLE]

The coefficients will be determined from the minimisation of a sum rule for the average s.p. energy

[TABLE]

and equivalently for $q^{\dagger}_{\rho}$. The even particle number equation relies on the usual RPA excitation operator RingSchuck

[TABLE]

The $X,~{}Y$ amplitudes are the solutions of another sum-rule defining an average excitation energy of the even systems

[TABLE]

The destruction operator $Q_{\nu}$ does not exactly kill the CCD ground state without introducing a generalization jemai13 but it kills it to very good approximation as studies of model cases have shown jemai13 . Therefore the even particle number equation (SCRPA) is the only point where the approach is not entirely consistent though the theory remains very performent as we will see below in the Application section. We, thus, will henceforth always suppose that

[TABLE]

The SCRPA equation corresponding to the minimisation of the mean excitation energy $\Omega_{\nu}$ in (7) can be written as

[TABLE]

with the normalisation of the amplitude given as usual by

[TABLE]

and

[TABLE]

where $\langle...\rangle=\langle Z|...|Z\rangle/\langle Z|Z\rangle$. The SCRPA equations are well documented in the literature DSD05 ; jemai05 ; tddm-scrpa ; jemai13 and we will not repeat their explicit form here. Let us simply say that $A$ and $B$ are functional of one and two particle density matrices when the Hamiltonian of the system is given by

[TABLE]

The first part of the Hamiltonian represents the kinetic energy and the second part the two body interaction with the anti-symmetrized matrix element

[TABLE]

From the minimisation of the sum-rule in eq.(5), we obtain two coupled equations

[TABLE]

or written as a matrix eigenvalues equation

[TABLE]

with

[TABLE]

where we supposed that hitherto we work in the Mean-Field (MF) basis with diagonal s.p. MF energies $\epsilon_{h},\epsilon_{p}$. The matrices ${{\cal C}}$ and ${{\cal D}}$ in (II) are obtained from the minimisation of the mean s.p. energy given in (5)

[TABLE]

Equation (II) is essentially already given in TohSchu13 . However, the way we solve this equation and, thus, couple it to the SCRPA of (15) is novel. Let us briefly describe the procedure. The coefficients ${{\cal C}}$ and ${{\cal D}}$ contain two and three body correlation functions. In particular they contain p-h operators which are given by the inversion of (6) valid because the $X,Y$ amplitudes in (15) form a complete orthonormal set of states

[TABLE]

and its hermitian conjugate. All other correlation functions which do not contain those ph operators and which are not of the $pp^{\prime}h$ type shall be discarded since they are supposedly less important. Commuting the destructor $Q$ to the right until they hit and kill the vacuum state $|Z\rangle$ leads to expressions of diverse correlation functions which only contain s.p. occupations $n_{h}$ and $n_{p}$ and RPA amplitudes $X,Y$. We want to call the coupled equations (15-17) and (20-22) ’even-odd SCRPA’ (eo-SCRPA). One may find more details in the Application section below. This procedure to obtain the s.p. occupation numbers is the essential new point of this work. It is clear that in this way eqs. (15) and (II) become coupled. In our earlier publications the s.p. occupation numbers appearing in the SCRPA equations have always been obtained in a different, in our opinion less natural way. We should also say the the formal expressions of SCRPA are not altered, only the way how the s.p. occupation probabilities in there are calculated is new.

It may be helpful at this point to discuss for instance the matrix ${{\cal D}}$ in (20,22) a little more and give a graphical representation. From the double commutator in ${{\cal D}}$, we retain only those terms where a particle state of the triplet operator on the right connects to the interaction and the same for the triplet on the left. In doing so, what is left from the interaction is a density operator $a^{{\dagger}}a$ for which we will make the diagonal approximation. Of course anti-symmetrization of the two particle indices will be fully respected. We then can make a graphical representation of the interaction process contained in ${{\cal D}}$ as shown in Fig. 1.

After the diagonalization of the matrix which implies a self-consistency on the occupancies, we can find the occupation numbers $n_{h}$ as

[TABLE]

and

[TABLE]

where the summation extends over all the amplitudes where $\lambda_{\mu}<E_{F}$ for hole state (or $\lambda_{\rho}>E_{F}$ for particle state). These simple expressions stem from the fact that, e.g., $a^{\dagger}_{h}$ commutes with $a^{\dagger}_{h}a_{p}a_{p^{\prime}}$. Please notice that these occupation numbers enter also the $A$ and $B$ matrices in Eq.(16). Again details of the procedure will become more clear in the applications we will give below.

Another way to find the same results for the occupancies is to define the Green Function (GF)

[TABLE]

from where we find the resonances as

[TABLE]

The mass operator ${{\cal M}}_{h}$ is obtained from the eq.(20), eliminating the amplitude $U$,

[TABLE]

The solution of (27) has obviously the same eigenvalues as (20) and then ${{\cal G}}_{h}$ can be written as

[TABLE]

where

[TABLE]

and ${{\cal M}}^{\prime}_{h}=d{{\cal M}}_{h}(\omega)/d\omega=-{{\cal C}}^{\dagger}(\omega-{{\cal D}})^{-2}{{\cal C}}$. We can easily check that $\sum_{\alpha}r_{\alpha}=1$ (the sum over all residua) and we can write the Green function dependent on time as

[TABLE]

Thus, we can find the s.p. occupation probabilities as

[TABLE]

Once we have the GF’s, we can calculate the ground state energy in the usual way via FeterW71

[TABLE]

In order to test our idea, we chose two models where we know the exact solution. The first application concerns the Lipkin model as an orientation to nuclear physics. The second one focuses on solid state physics where the Hubbard model is chosen.

III Applications

III.1 The Lipkin Model

The single-particle space of the Lipkin model consists of two fermion levels, each of which has a N-fold degeneracy see, e.g.,RingSchuck . The upper (lower) level has the energy of $\frac{e}{2}$ ($-\frac{e}{2}$). The Hamiltonian of the Lipkin model is given by

[TABLE]

with $e$ is the inter-shell spacing, $V$ is the coupling constant and

[TABLE]

with $2J_{0}=\hat{n}_{1}-\hat{n}_{0}$, $\hat{n}_{i}=\sum_{m}c^{\dagger}_{im}c_{im}$ and $N$ is the number of particles equivalent to the degeneracies of the shells. We consider the odd excitation operator as (1), that is

[TABLE]

with the minimisation of the sum rule in (5). Based on the solution of the SCRPA equations jemai11 ; RingSchuck ; jemai13 with the definition of the pair excitation operator as $Q^{+}=(XJ_{+}-YJ_{-})/\sqrt{\langle-2J_{0}\rangle}$ (with $J_{+}=\sqrt{\langle-2J_{0}\rangle}(XQ^{+}-YQ)$), we obtain the $X,~{}Y$ amplitudes as being the solutions of SCRPA equations. From the minimisation of expression (5), we obtain a $2\times 2$ matrix eigenvalue equation (see appendix A), with the Hamiltonian and norm matrices,

[TABLE]

where we define the elements of the two matrices as

[TABLE]

and the corresponding secular equation

[TABLE]

where the eigenvalues $\lambda$ are given in App. (57). In the above equations (39) the correlation functions are expressed by the RPA amplitudes $X,Y$ in the way it is described in section (II) and App. A. The correlation functions which contain quadratic forms of occupation number operators as $\langle J_{0}J_{0}\rangle$ in above equation can in principle be expressed by the RPA amplitudes as well but leading to heavier expressions. Usually, we, therefore will employ the factorization approximation leading in the present case to $\langle J_{0}J_{0}\rangle\simeq\langle J_{0}\rangle^{2}$ what mostly turns out to be quite satisfactory. However, in the case of the Lipkin model one also can use the Casimir relation to close the system of equations, see App. A where also more details of the procedure are given. The results are shown in Figs. 2 - 7. They concern in the order: i) the expectation value $\langle J_{0}\rangle$ of the difference of populations in upper and lower level, ii) the square of this quantity, iii) the first excitation energy, iv) the correlation energy, and v) the excitation energy between the system with $N\pm 1$ and $N$ particles, as $\lambda_{\pm}=E^{N\pm 1}_{\alpha}-E^{N}_{0}$. All quantities are very well reproduced throughout couplings up to the critical value $\chi=\chi_{\rm crit.}$ where the standard RPA breaks down and the system wants to change to the ’deformed’ basis. However, even values slightly beyond $\chi_{\rm crit.}=1$ are still quite acceptable. All quantities for $N=2$ are reproduced exactly. By some lucky accident the occupancies even for $N=4$ come out to be exact (as shown in Figs. 2, 5 and 7). In Fig.2, Fig.3, and Fig.4, in the panels with $N=10$ and $N=20$, we also show the results of the Catara method Catara for the calculation of the occupation numbers and first excited state (as a reminder, let us mention that using the Catara method for the occupation numbers has been named the SCRPA method in the past; we keep the same name while getting the occupations from the selfconsistent odd RPA). One can thus appreciate the important improvement obtained with the method of the present work where even and odd RPA’s are coupled.

III.2 The Hubbard Model

The Hubbard model is widely used to deal with the physics of strongly correlated electrons. Since the model can be solved exactly in one dimension (1D) and for small cluster sizes, it is very useful for theoretical investigations jemai05 . To be precise, our ”Hubbard model” is a 6-site system at half filling with periodic boundary condition, described by the usual Hamiltonian jemai05 ; Hubbard :

[TABLE]

Here, $\hat{n}_{i\sigma}=c_{i\sigma}^{\dagger}c_{i\sigma}$, $c_{i\sigma}^{\dagger}$ and $c_{i\sigma}$ are the creation and annihilation operators for an electron at site $i$ with spin $\sigma$, $U$ is the on-site (spin-independent) interaction, $-t$ is the hopping term of the kinetic energy. The eigenstates of the system can be expressed as linear combinations of Slater determinants. The Hamiltonian is rewritten in plane wave basis,

[TABLE]

with the transformation

[TABLE]

where $\hat{n}_{\bf{k},\sigma}=a^{\dagger}_{\bf{k},\sigma}a_{\bf{k},\sigma}$, $\varepsilon_{\bf{k}}=-2t\cos\left(ka\right)$, which are, respectively, the number operator of particles of the mode $({\bf{k}},\,\sigma)$ and the energies of one particle on a lattice with a the parameter of the lattice which is taken as $a=1$. For a problem with $N$sites, the condition of periodicity is given by $c_{N+1,\sigma}=c_{1,\sigma}$. This implies that $e^{-ik\,N}=1$, hence the values taken by $k$ will be $k=\frac{2\,\pi}{N}\,n$. In addition, the first Brillouin zone is defined on the field where $-\pi\leqslant k<\pi$, which gives us the values of $n$ as $\frac{-N}{2}\leqslant n<\frac{N}{2}$.

For the six sites, we have the possible states with the following wave vectors:

[TABLE]

and with the kinetic energies (see Fig.8), respectively,

[TABLE]

The transfer wave vector($q_{ph}=k_{p}-k_{h}$) takes the possible values as shown in the Table 1.

At this point we proceed exactly as in the case of the Lipkin model: The excitation operator for the even system is given by

[TABLE]

with $~{}K^{\pm}_{ph\sigma}=J^{\pm}_{ph\sigma}/\sqrt{N_{ph\sigma}}$, $~{}J^{+}_{ph\sigma}=a^{{\dagger}}_{p\sigma}a_{h\sigma}$, $~{}N_{ph\sigma}=n_{h\sigma}-n_{p\sigma}$. With the inversion

[TABLE]

we can calculate the mean values needed for the matrix elements of the SCRPA equations for the even particle number case

[TABLE]

where we replaced the ”ph” operators by the RPA creation and destruction operators from the inversion (23) and then commute the $Q$ operators to the right until they kill the ground state. All matrices become functional of the occupancies $n_{h}$ and $n_{p}$ and $X,Y$ amplitudes in analogy to what was the case in the Lipkin model and, thus, the diagonalisation process implies at the same time an iteration on the occupancies and the amplitudes.

For the odd particle number case, we make again the following ansatz

[TABLE]

From there, we can, as outlined in the general section II and as just now for the case of the Lipkin model calculate the occupation numbers. For more details, see App. B. The results for the occupation numbers are again very satisfying, see Fig. 9. Also the excitation energies of the even particle number system, see Fig. 10 are very well reproduced. In Fig. 11 we show the ground state energies for the exact case compared to the eo-SCRPA solution.

One should notice that there is barely an improvement using the eo-SCRPA versus the standard SCRPA because the latter produced already excellent results. So, we do not show the old SCRPA results again. It is not quite clear why there is this difference between the Lipkin and Hubbard models. Probably the fact that in Lipkin model, contrary to the Hubbard model, one uses collective ph operators makes it more difficult to fullfill the Pauli principle. So the performance of one or the other approach seems to depend on the situation.

IV Conclusion

In this work, we coupled even and odd particle number RPA self-consistently. Both systems are based on the same correlated RPA ground state. From the odd system, we get the occupation numbers, odd particle excitation energies, and the ground state energies whereas from the even SCRPA equations we get the excitation energies of the even system and transition probabilities. To make things clear, we should mention again that the SCRPA employed here has the same mathematical structure as the one used before jemai05 , only the single particle occupation probabilities are now calculated via the odd selfconsistent RPA whereas they were obtained before via the so-called Catara method Catara . Both even and odd systems are coupled through non-linear equations which both contain the RPA amplitudes $X,Y$ and the s.p. occupation numbers $n_{k}$ in a non-linear way. We called this system of equations ’even-odd’ SCRPA (eo-SCRPA). Applications to the Lipkin model and a six sites Hubbard ring at half filling gave very satisfying results for all quantities. The equations are relatively complex due to their non-linearity but they should be solvable with modern computers for realistic problems such as, e.g., the calculation of collective states in nuclei. The equations to be solved seem not to be of higher numerical complexity than, e.g., the Brueckner Hartree-Fock equations which have been solved a number of times for nuclei. The coupling of even and odd RPA’s has a couple of advantages: it gives richer results, i.e., excitation energies of even and odd particle number systems; there is a natural way how to obtain the ground state energy via the s.p. Green’s function and, last but not least, the results seem to be promising.

V Acknowledgements

We are grateful for long-standing collaboration on SCRPA with D. Delion, J. Dukelsky, and M. Tohyama.

Appendix A Equation of Motion for odd particle number operator for Lipkin Model

We consider the odd excitation operator as

[TABLE]

and the coefficients will be determined from minimisation of expression (5). Based on the solution of the SCRPA equations with the definition of the pair excitation operator as $Q^{\dagger}=(XJ_{+}-YJ_{-})/\sqrt{\langle-2J_{0}\rangle}$ (with $J_{+}=\sqrt{\langle-2J_{0}\rangle}(XQ^{\dagger}-YQ)$), the $X,~{}Y$ amplitudes are the solutions of the SCRPA equations with $H$ of the Lipkin Hamiltonian (34). From the minimisation of (5), we obtain a $2\times 2$ matrix eigenvalue equation. The norm matrix is given by

[TABLE]

where we have used the inversion (23) and the killing condition $Q|0\rangle=0$. For the first Hamiltonian element we have

[TABLE]

and for the off diagonal elements

[TABLE]

And the anti-commutator for ${{\cal H}}_{11}$ is given by

[TABLE]

with, for example

[TABLE]

where we have again used the inversion (23) and the killing condition $Q|0\rangle=0$.

The correlation functions which contain quadratic forms of occupation number operators as $\langle J_{0}J_{0}\rangle$ in above equation can in principle be expressed by the RPA amplitudes as well but leading to heavier expressions. Usually, we, therefore will employ the factorization approximation leading in the present case to $\langle J_{0}J_{0}\rangle\simeq\langle J_{0}\rangle^{2}$. However, in the Lipkin model one also can use the Casimir relation

[TABLE]

Then all matrix elements ${{\cal H}}_{ij}$ become functions of $\langle J_{0}\rangle$ and the RPA amplitudes $X,Y$. The eigenvalue problem can therefore be solved leading to a self-consistency problem for $\langle J_{0}\rangle$ and the RPA amplitudes which are obtained from the SCRPA equations (17) jemai13 . The occupation numbers are then given by

[TABLE]

where $\lambda_{\pm}$ are the eigenvalues of the $2\times 2$ matrix problem,

[TABLE]

with $\beta=e-VXY(N-4)-VXY(N-4)(1+2Y^{2})\frac{\langle J_{0}\rangle}{n_{11}}$. Thus,

[TABLE]

Appendix B Equation of Motion for Hubbard Model

For the Hubbard model (41) we define the odd excitation operator as in Eq.(1),

[TABLE]

with $J^{+}_{ph-}=a^{\dagger}_{p-}a_{h-}$ and $\sigma=\uparrow,\downarrow=+,-$. Remembering the notations for the occupation probabilities

[TABLE]

we have $n_{k_{2}\sigma}=n_{k_{3}\sigma}$, $n_{k_{4}\sigma}=n_{k_{5}\sigma}$, $n_{k_{2}\sigma}=1-n_{k_{3}\sigma}$ and $n_{k_{1}\sigma}=1-n_{k_{6}\sigma}$. This gives

[TABLE]

The term without interaction $H_{0}=\sum_{k\sigma}\varepsilon_{k}\hat{n}_{k\sigma}$ is given by

[TABLE]

with ${{\cal N}}_{p^{\prime}ph}=\langle(1-\hat{n}_{p^{\prime}+})(-2J^{0}_{ph,-})\rangle+\langle J^{+}_{ph,-}J^{-}_{hp,-}\rangle$. The term in the Hamiltonian for the transfer $q=0$, $H_{q=0}=\frac{U}{6}\sum_{kk^{\prime}}\hat{n}_{k+}\hat{n}_{k^{\prime}-}$ leads to

[TABLE]

with $\sum_{k}\hat{n}_{k\sigma}=\sum_{p}\hat{n}_{p\sigma}+\sum_{h}\hat{n}_{h\sigma}=3$ in the half-filled case. Now let us calculate the elements ${{\cal C}}_{p^{\prime}ph}$ for the first row (or column) as

[TABLE]

The elements of the matrix except the first row (or column) are given as follows

[TABLE]

In the following, as already discussed several times, we retain from (65) only those terms where the particle states of the left and right triple operators in ${{\cal D}}$ connect to the interaction. The remaining density operator from the interaction is approximated by its diagonal form. This leads to expressions evaluated in (66) below. First let us discuss what kind of terms we are neglecting in this way. It should be noted that the terms of type $\langle J^{\pm}_{ph}J^{\pm}_{p^{\prime}h^{\prime}}J^{\pm}_{p^{\prime\prime}h^{\prime\prime}}\rangle=0$, $\langle J^{\pm}_{ph}S_{p_{1}p_{2}}J^{\pm}_{p^{\prime\prime}h^{\prime\prime}}\rangle$ are probably small (with $S_{p_{1}p_{2}}=a^{\dagger}_{p_{1}}a_{p_{2}}$ for $p_{1}\neq p_{2}$) and $\langle J^{\pm}_{ph}S_{h_{1}h_{2}}J^{\pm}_{p^{\prime\prime}h^{\prime\prime}}\rangle$ also small (with $S_{h_{1}h_{2}}=a^{\dagger}_{h_{1}}a_{h_{2}}$ for $h_{1}\neq h_{2}$) in eq.(65). As shown in jemai13 , the term $\langle SJ\rangle=0$ and $\langle SS\rangle$ are small. Only the terms non-zero in eq.(65) like $\langle J^{\pm}_{ph+}n_{k\pm}J^{\pm}_{p^{\prime}h^{\prime}-}\rangle$ which can be calculated as shown in (66) are kept. With the short hand notation $ph\sigma\equiv i$, $k\sigma\equiv k$, $\hat{N}_{i}=\hat{n}_{h\sigma}-\hat{n}_{p\sigma}$ and $N_{i}=n_{h\sigma}-n_{p\sigma}$, we can evaluate the following terms

[TABLE]

where we used $Q^{\dagger}_{\nu}=\sum_{l}(X^{\nu}_{l}J^{+}_{l}+Y^{\nu}_{l}J^{-}_{l})N^{-1/2}_{l}$ and the commutators

[TABLE]

This entails, $\langle[Q_{\nu},Q^{\dagger}_{\nu^{\prime}}]\rangle=1$ and

[TABLE]