Coupled Self-Consistent RPA Equations for Even and Odd Particle Numbers. Tests with Solvable Models
aa r X i v : . [ nu c l - t h ] A ug Coupled Self-Consistent RPA Equations for Even and Odd Particle Numbers.Tests with Solvable Models.
M. Jema¨ı
Laboratoire des mat´eriaux avanc´es et ph´enom`enes quantiques, FST,Universit´e Tunis El-Manar 2092 El-Manar, Tunis, Tunisia. andISSATM, Universit´e de Charthage, Avenue de la R´epublique P.O. Box 77 - 1054 Amilcar, Tunis, Tunisia. ∗ P. Schuck
Institut de Physique Nucl´eaire d’Orsay, Universit´e Paris-Sud, CNRS–IN2P315, Rue Georges Clemenceau, 91406 Orsay Cedex, France. andUniv. Grenoble Alpes, CNRS, LPMMC, 38000 Grenoble, France † (Dated: August 16, 2019)Coupled equations for even and odd particle number correlation functions are set up via theequation 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 particlenumber case the single particle occupation probabilities are obtained in a self-consistent way. Thisis the essential new procedure of this work. Both, even and odd particle number cases are basedon the same correlated vacuum and, thus, are coupled equations. Applications to the Lipkin modeland the 1D Hubbard model give very good results. PACS numbers: 21.60.-n, 21.60.Fw, 71.10.-w, 75.10.Jm
I. INTRODUCTION
Developments of Many-Body approaches for stronglycorrelated systems is an active field of research. In thepast we developed an RPA theory which goes beyond thestandard one and which is based on a correlated groundstate. To lowest order this leads to RPA equations wherethe single particle (s.p.) occupation numbers are notthe uncorrelated (Hartre-Fock (HF)) ones but correlatedones which are obtained in a self-consistent way fromthe RPA solution. It is known as ’re-normalized’ RPA(r-RPA) [1, 2]. However, in general the RPA equationscontain additionally vertex (e.g., screening) correctionswhich also can be obtained self-consistently from theRPA solution. The whole procedure has been dubbedSelf-Consistent RPA (SCRPA) [3] and references in there.SCRPA can also be seen as a sub-product of an even moregeneral approach which is the Time-Dependent DensityMatrix (TDDM) theory based on a decoupling of theBBGKY hierarchy of one, two, etc. correlation func-tions [4, 5]. In the past, there was always a certaindebate how to include the single particle (s.p.) occu-pations into the SCRPA scheme. Since the latter ares.p. quantities and the RPA gives rise to two body cor-relations, it was not completely evident in which way toclose the system of equations. However, already Rowepromoted the so-called particle number operator methodto obtain self-consistent s.p. occupations [6]. Later thiswas further elaborated by F. Catara [7] and it has be-come known as the ’Catara method’ since then. This ∗ Electronic address: [email protected] † Electronic address: [email protected] method mostly works quite well [8] but also fails moreor less in particular cases [3]. It is for this reason thatin this work we elaborate a different scheme which seemsto be more natural since it is also based on the equa-tion of motion (eom) method and employs the so-calledodd particle number RPA (o-RPA) recently proposed byone of the authors plus collaborators [9]. The latter canalso be reformulated as a Dyson equation for the s.p.Green’s function (GF) with a self-energy obtained fromthe eom for the 2particle-1hole (2p-1h) and 2h-1p GF.Since both, the SCRPA and o-RPA will be based on thesame correlated vacuum, naturally even and odd parti-cle number channels become coupled. We may coin thisscheme eo-SCRPA. We will apply this approach to twoexactly solvable model cases: the Lipkin model and the1D Hubbard model. In both cases the results turn outto be promising.The paper is organized as follows: in section II wepresent the general theory. In section III, applicationsto the Lipkin and Hubbard models are given. Section IVcontains the conclusions and details of our procedures arepresented in the Appendices.
II. GENERAL THEORY
As mentioned in the introduction, we will base our ap-proach on the coupling of the even and odd particle num-ber eom. The latter will be obtained from the followingansatz q † µ = X h x µh a h + X pp ′ h ′ U µpp ′ h ′ a † h ′ a p ′ a p ,q † ρ = X p x ρp a † p + X p ′ h ′ h U ρp ′ h ′ h a † h a † h ′ a p ′ . (1)where the indices ’p, h’ refer to s.p. states ’above’ and’below’ the Fermi surface, respectively, and a † k , a k are thefermion creation and annihilation operators. The equa-tion for the s.p. basis in which the equations will beworked out will be given below. Our ’quasi-particle’ op-erator in (1) has the good quality that its destructor ex-actly kills the so-called Coupled-Cluster-Doubles (CCD)wave function: q µ | Z i = 0 , (2)with | Z i = exp X pp ′ hh ′ Z pp ′ hh ′ a † p a † p ′ a h ′ a h | HF i , (3)where | HF i is the Hartree-Fock (HF) Slater determinantand the amplitudes must full-fill the following relations X h x µ ∗ h Z pp ′ hh ′ = U µ ∗ pp ′ h ′ , X p x ρ ∗ p Z pp ′ hh ′ = U ρ ∗ p ′ h ′ h . (4)The coefficients will be determined from the minimisationof a sum rule for the average s.p. energy λ µ = 12 h{ q µ , [ H, q † µ ] }ih{ q µ , q † µ }i = 1 h |{ q µ , q † µ }| i X α ( E N +1 α − E N ) |h | q µ | α i| (5)and equivalently for q † ρ . The even particle number equa-tion relies on the usual RPA excitation operator [10] Q † ν = X ph X νph a † p a h − Y νph a † h a p . (6)The X, Y amplitudes are the solutions of another sum-rule defining an average excitation energy of the evensystemsΩ ν = 12 h | [ Q ν , [ H, Q † ν ]] | ih | [ Q ν , Q † ν ] | i = 1 h | [ Q ν , Q + ν ] | i X µ ( E µ − E ) |h | Q ν | µ i| . (7)The destruction operator Q ν does not exactly kill theCCD ground state without introducing a generalization[3] but it kills it to very good approximation as studiesof model cases have shown [3]. Therefore the even par-ticle number equation (SCRPA) is the only point wherethe approach is not entirely consistent though the the-ory remains very performent as we will see below in theApplication section. We, thus, will henceforth alwayssuppose that Q | Z i = 0 . (8) The SCRPA equation corresponding to the minimisationof the mean excitation energy Ω ν in (7) can be writtenas (cid:18) A B − B ∗ − A ∗ (cid:19) (cid:18) X ν Y ν (cid:19) = Ω ν (cid:18) X ν Y ν (cid:19) (9)with the normalisation of the amplitude given as usualby X ph (cid:0) | X νph | − | Y νph | (cid:1) = 1 (10)and A ph,p ′ h ′ = h [ a † h a p , [ H, a † p ′ a h ′ ]] i√ n h − n p √ n h ′ − n p ′ ,B ph,p ′ h ′ = − h [ a † h a p , [ H, a † h ′ a p ′ ]] i√ n h − n p √ n h ′ − n p ′ (11)where h ... i = h Z | ... | Z i / h Z | Z i . The SCRPA equationsare well documented in the literature [1, 3, 5, 8] and wewill not repeat their explicit form here. Let us simplysay that A and B are functional of one and two particledensity matrices when the Hamiltonian of the system isgiven by H = X kk ′ t kk ′ a † k a k ′ + 14 X klmn ¯ v klmn a † k a † l a n a m . (12)The first part of the Hamiltonian represents the kineticenergy and the second part the two body interaction withthe anti-symmetrized matrix element¯ v klmn = h kl | v | mn i − h kl | v | nm i . From the minimisation of the sum-rule in eq.(5), we ob-tain two coupled equations X h ′ ǫ hh ′ x µh ′ + X pp ′ h ′ C h,pp ′ h ′ U µpp ′ h ′ = λ µ x µh (13) X h ′ C ∗ pp ′ h,h ′ x µh ′ + X p p h D pp ′ h,p p h U µp p h = λ µ U µpp ′ h or written as a matrix eigenvalues equation (cid:18) ǫ CC † D (cid:19) (cid:18) x µ U µ (cid:19) = λ µ (cid:18) x µ U µ (cid:19) (14)with ǫ hh ′ = h{ a h , h H, a † h ′ i }i = ǫ h δ hh ′ (15)where we supposed that hitherto we work in the Mean-Field (MF) basis with diagonal s.p. MF energies ǫ h , ǫ p .The matrices C and D in (14) are obtained from the min-imisation of the mean s.p. energy given in (5) C ∗ pp ′ h ′ ,h = h n a † h ′ a p ′ a p , h H, a † h io i p N pp ′ h ′ D pp ′ h ′ ,p p h = h n a † h ′ a p ′ a p , (cid:2) H, a † p a † p a h (cid:3)o i p N pp ′ h ′ p N p p h N pp ′ h ′ = h n a † h ′ a p ′ a p , a † p a † p ′ a h ′ o i (16)Equation (14) is essentially already given in [9]. How-ever, the way we solve this equation and, thus, coupleit to the SCRPA of (9) is novel. Let us briefly describethe procedure. The coefficients C and D contain two andthree body correlation functions. In particular they con-tain p-h operators which are given by the inversion of (6)valid because the X, Y amplitudes in (9) form a completeorthonormal set of states a † p a h = p n h − n p X ν X νph Q + ν + Y νph Q ν (17)and its hermitian conjugate. All other correlation func-tions which do not contain those ph operators and whichare not of the pp ′ h type shall be discarded since they aresupposedly less important. Commuting the destructor Q to the right until they hit and kill the vacuum state | Z i leads to expressions of diverse correlation functionswhich only contain s.p. occupations n h and n p and RPAamplitudes X, Y . We want to call the coupled equations(9-11) and (14-16) ’even-odd SCRPA’ (eo-SCRPA). Onemay find more details in the Application section below.This procedure to obtain the s.p. occupation numbers isthe essential new point of this work. It is clear that inthis way eqs. (9) and (14) become coupled. In our ear-lier publications the s.p. occupation numbers appearingin the SCRPA equations have always been obtained in adifferent, in our opinion less natural way. We should alsosay the the formal expressions of SCRPA are not altered,only the way how the s.p. occupation probabilities inthere are calculated is new.It may be helpful at this point to discuss for instancethe matrix D in (14,16) a little more and give a graph-ical representation. From the double commutator in D ,we retain only those terms where a particle state of thetriplet operator on the right connects to the interactionand the same for the triplet on the left. In doing so,what is left from the interaction is a density operator a † a for which we will make the diagonal approximation.Of course anti-symmetrization of the two particle indiceswill be fully respected. We then can make a graphicalrepresentation of the interaction process contained in D as shown in Fig. 1. After the diagonalization of the ma-trix which implies a self-consistency on the occupancies,we can find the occupation numbers n h as n h = h a † h a h i = X µ |h{ a † h , q † h,µ }i| = X µ | x µh | (18) D p p h,p ′ p ′ h ′ → ✲ ✲ p p h (cid:0)(cid:0)✒❅❅❘ (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)✒ r ❅❅❅❅❅❘✲ (cid:0)(cid:0)✒❅❅❘ p ′ p ′ h ′ − ( p ↔ p ) − ( p ′ ↔ p ′ ) FIG. 1: Schematic representation of the three body interac-tion D . The wiggly line stands for the (collective) ph-modes.The full dot represents the two body interaction. and n p = h a † p a p i = X ρ | x ρp | (19)where the summation extends over all the amplitudeswhere λ µ < E F for hole state (or λ ρ > E F for particlestate). These simple expressions stem from the fact that,e.g., a † h commutes with a † h a p a p ′ . Please notice that theseoccupation numbers enter also the A and B matrices inEq.(10). Again details of the procedure will become moreclear in the applications we will give below.Another way to find the same results for the occupanciesis to define the Green Function (GF) G h ( ω ) = 1 ω − ǫ h − M h ( ω ) (20)from where we find the resonances as λ α − ǫ h − M h ( λ α ) = 0 (21)The mass operator M h is obtained from the eq.(14),eliminating the amplitude U , M h = X pp ′ h ′ ,p p h C ∗ h,pp ′ h ′ ( ω − D ) − pp ′ h ′ ,p p h C p p h ,h (22)The solution of (21) has obviously the same eigenvaluesas (14) and then G h can be written as G h ( ω ) = X α r α ω − λ α (23)where r α = 11 − M ′ | ω = λ α (24)and M ′ h = d M h ( ω ) /dω = −C † ( ω −D ) − C . We can easilycheck that P α r α = 1 (the sum over all residua) and wecan write the Green function dependent on time as i G h ( t − t ′ ) = − θ ( t − t ′ ) X α ( λ α The single-particle space of the Lipkin model consistsof two fermion levels, each of which has a N-fold degener-acy see, e.g.,[10]. The upper (lower) level has the energyof e ( − e ). The Hamiltonian of the Lipkin model is givenby H = eJ − V (cid:0) J + J − (cid:1) (28)with e is the inter-shell spacing, V is the coupling con-stant and J = 12 N X m =1 (cid:16) c † m c m − c † m c m (cid:17) ,J + = N X m =1 c † m c m , J − = ( ˆ J + ) † , (29)with 2 J = ˆ n − ˆ n , ˆ n i = P m c † im c im and N is thenumber of particles equivalent to the degeneracies of theshells. We consider the odd excitation operator as (1),that is q † µ = 1 N X m x µ m c m + U µ m J + c † m (30)with the minimisation of the sum rule in (5). Based onthe solution of the SCRPA equations [3, 10, 11] with thedefinition of the pair excitation operator as Q + = ( XJ + − Y J − ) / p h− J i (with J + = p h− J i ( XQ + − Y Q )), weobtain the X, Y amplitudes as being the solutions ofSCRPA equations. From the minimisation of expression(5), we obtain a 2 × H ij = (cid:18) H H H H (cid:19) and N = (cid:18) n n n n (cid:19) (31) χ < - J > / N sRPASCRPAeo-SCRPAExactN=4 χ < - J > / N sRPASCRPAeo-SCRPAExactN=10 χ < - J > / N sRPASCRPAeo-SCRPAExactN=20 χ < - J > / N sRPASCRPAeo-SCRPAExact N=100 FIG. 2: The difference between occupation number of thetwo level in Lipkin model, normalized by N as a function of χ = V ( N − 1) for N = 4 , , , N = 4. Also, we present the results ofCatara method for N = 10. χ < J > / N sRPAeo-SCRPASCRPAExactN=4 χ < J > / N sRPASCRPAeo-SCRPAExactN=10 χ < J > / N sRPASCRPAeo-SCRPAExactN=20 χ < J > / N sRPASCRPAeo-SCRPAExactN=100 FIG. 3: Same as Fig.2 but for the square of the differencebetween occupation number of the two level in Lipkin model,normalized by N . χ Ω sRPASCRPAeo-SCRPAExact N=4 χ Ω sRPASCRPAeo-SCRPAExactN=20 χ Ω sRPASCRPAeo-SCRPAExactN=100 χ Ω sRPASCRPAeo-SCRPAExactN=200 FIG. 4: Same as Fig.2 but for the first excited state for N = 4 , , , N → ∞ limit. Also, we present the results of Cataramethod for N = 20. χ -0.6-0.4-0.20 E c o rr eo-SCRPAExactN=4 χ -0.8-0.6-0.4-0.20 E c o rr eo-SCRPAExactN=10 χ -1.2-1-0.8-0.6-0.4-0.20 E c o rr eo-SCRPAExact N=20 χ -2.5-2-1.5-1-0.50 E c o rr eo-SCRPAExactN=100 FIG. 5: The correlation energy as a function of χ = V ( N − N = 4 , , , 100 with eo-SCRPA (blue dashed line)compared to the exact solution (full black line). Note againthat for N = 4 the exact result is obtained with our approach. χ -2-1.8-1.6-1.4-1.2-1-0.8-0.6 λ + eo-SCRPAExact N=4 χ -0.8-0.7-0.6-0.5 λ + eo-SCRPAExactN=10 FIG. 6: Excitation energy between the system N + 1 and N particles as a function of χ = V ( N − 1) for N = 4 , 10 witheo-SCRPA (blue dashed line) (A9) compared to the exactsolution (full black line) λ + = E N +1 α − E N . where we define the elements of the two matrices as n = 1 N X m h{ c m , c † m }i n = n = 1 N X m h{ c m , J + c † m }i n = 1 N X m h{ c m J − , J + c † m }iH = 1 N X m h{ c m , [ H, c † m ] }iH = H = 1 N X m h{ c m J − , [ H, c † m ] }iH = 1 N X m h{ c m J − , [ H, J + c † m ] }i (32)and the corresponding secular equationdet (cid:26)X i ′ j ′ N − / ii ′ H i ′ j ′ N − / j ′ j − λI (cid:27) = 0 (33) χ λ - eo-SCRPAExactN=4 χ λ - eo-SCRPAExactN=10 FIG. 7: Same as Fig. 6 but for the excitation energy betweenthe system N − N particles as a function of χ = V ( N − N = 4 , 10. Note that for N = 4 the exact result λ − = E N − α − E N is obtained with our approach eq. (A9). where the eigenvalues λ are given in App. (A9). In theabove equations (33) the correlation functions are ex-pressed by the RPA amplitudes X, Y in the way it isdescribed in section (II) and App. A. The correlationfunctions which contain quadratic forms of occupationnumber operators as h J J i in above equation can inprinciple be expressed by the RPA amplitudes as wellbut leading to heavier expressions. Usually, we, there-fore will employ the factorization approximation leadingin the present case to h J J i ≃ h J i what mostly turnsout to be quite satisfactory. However, in the case of theLipkin model one also can use the Casimir relation toclose the system of equations, see App. A where alsomore details of the procedure are given. The results areshown in Figs. 2 - 7. They concern in the order: i) theexpectation value h J i of the difference of populations inupper and lower level, ii) the square of this quantity, iii)the first excitation energy, iv) the correlation energy, andv) the excitation energy between the system with N ± N particles, as λ ± = E N ± α − E N . All quantitiesare very well reproduced throughout couplings up to the ε F ππ/3π/3π/3π/3 k=k =2k=k =−2k=k =−k=k =−k=k = k=k =0 FIG. 8: Hatree Fock States at U = 0 for the chain with6 sites at half filling and projection of spin m s = 0. Theoccupied states are represented by the full arrows and thosenot occupied are represented by the dashed arrows. critical value χ = χ crit . where the standard RPA breaksdown and the system wants to change to the ’deformed’basis. However, even values slightly beyond χ crit . = 1are still quite acceptable. All quantities for N = 2 arereproduced exactly. By some lucky accident the occu-pancies even for N = 4 come out to be exact (as shownin Figs. 2, 5 and 7). In Fig.2, Fig.3, and Fig.4, in thepanels with N = 10 and N = 20, we also show the resultsof the Catara method [7] for the calculation of the occu-pation numbers and first excited state (as a reminder, letus mention that using the Catara method for the occu-pation numbers has been named the SCRPA method inthe past; we keep the same name while getting the occu-pations from the selfconsistent odd RPA). One can thusappreciate the important improvement obtained with themethod of the present work where even and odd RPA’sare coupled. B. The Hubbard Model The Hubbard model is widely used to deal with thephysics of strongly correlated electrons. Since the modelcan be solved exactly in one dimension (1D) and for smallcluster sizes, it is very useful for theoretical investiga-tions [8]. To be precise, our ”Hubbard model” is a 6-sitesystem at half filling with periodic boundary condition,described by the usual Hamiltonian [8, 13]: H = − t X h i,j i ,σ c † iσ c jσ + U X i,σ ˆ n i,σ ˆ n i, − σ . (34)Here, ˆ n iσ = c † iσ c iσ , c † iσ and c iσ are the creation and an-nihilation operators for an electron at site i with spin σ , U is the on-site (spin-independent) interaction, − t isthe hopping term of the kinetic energy. The eigenstatesof the system can be expressed as linear combinationsof Slater determinants. The Hamiltonian is rewritten inplane wave basis, H = X k σ ε k ˆ n k σ + U N X kk ′ q σ a † k σ a k + q σ a † k ′ − σ a k ′ − q − σ (35)with the transformation c j,σ = 1 √ N X k a k ,σ e − i k x j , (36)where ˆ n k ,σ = a † k ,σ a k ,σ , ε k = − t cos ( ka ), which are,respectively, the number operator of particles of the mode( k , σ ) and the energies of one particle on a lattice with athe parameter of the lattice which is taken as a = 1. For aproblem with N sites, the condition of periodicity is givenby c N +1 ,σ = c ,σ . This implies that e − ik N = 1, hence thevalues taken by k will be k = πN n . In addition, the firstBrillouin zone is defined on the field where − π k < π ,which gives us the values of n as − N n < N .For the six sites, we have the possible states with thefollowing wave vectors: k = 0 , k = − k = π , k = − k = 2 π , k = − π (37)and with the kinetic energies (see Fig.8), respectively, ε k = − ε k = 2 t, ε k = ε k = − ε k = − ε k = t. (38)The transfer wave vector( q ph = k p − k h ) takes the possiblevalues as shown in the Table I. q = ± π q = ± π q = ± π → q = + π → q = − π → q = − π → q = + π → q = + π → q = + π → q = − π → q = − π → q = − π TABLE I: The various momentum transfers in the 6 sites case. At this point we proceed exactly as in the case of theLipkin model: The excitation operator for the even sys-tem is given by Q † ν = X phσ X νphσ K + phσ − Y νphσ K − hpσ (39)with K ± phσ = J ± phσ / p N phσ , J + phσ = a † pσ a hσ , N phσ = n hσ − n pσ . With the inversion J − hpσ = p N phσ X ν (cid:0) X νphσ Q ν + Y νphσ Q † ν (cid:1) J + phσ = p N phσ X ν (cid:0) Y νphσ Q ν + X νphσ Q † ν (cid:1) . (40)we can calculate the mean values needed for the matrixelements of the SCRPA equations for the even particlenumber case h J + p ′ h ′ σ ′ J − hpσ i = p N p ′ h ′ σ ′ N phσ X ν Y νp ′ h ′ σ ′ Y νphσ , h J − h ′ p ′ σ ′ J + phσ i = p N p ′ h ′ σ ′ N phσ X ν X νp ′ h ′ σ ′ X νphσ , h J + p ′ h ′ σ ′ J + phσ i = p N p ′ h ′ σ ′ N phσ X ν Y νp ′ h ′ σ ′ X νphσ , h J − h ′ p ′ σ ′ J − hpσ i = p N p ′ h ′ σ ′ N phσ X ν X νp ′ h ′ σ ′ Y νphσ ,(41) k sRPAeo-SCRPAExact n k n k =n k FIG. 9: Occupation numbers as function of the interaction U/t for various values of the momenta k = − π , k = − π/ k = 2 π/ k = 2 π/ k = − π/ n k = 1 − n k and n k = n k = 1 − n k = 1 − n k where we replaced the ”ph” operators by the RPA cre-ation and destruction operators from the inversion (17)and then commute the Q operators to the right until theykill the ground state. All matrices become functional ofthe occupancies n h and n p and X, Y amplitudes in anal-ogy to what was the case in the Lipkin model and, thus,the diagonalisation process implies at the same time aniteration on the occupancies and the amplitudes.For the odd particle number case, we make again thefollowing ansatz q † h,µ = x µh a h + + X p ′ ph U µp ′ ph a † p ′ + J + ph − q † p,ρ = x ρp a † p + + X p ′ h ′ h U ρp ′ h ′ h a † h + J − h ′ p ′ − . (42)From there, we can, as outlined in the general section IIand as just now for the case of the Lipkin model calculatethe occupation numbers. For more details, see App. B.The results for the occupation numbers are again verysatisfying, see Fig. 9. Also the excitation energies of theeven particle number system, see Fig. 10 are very well re-produced. In Fig. 11 we show the ground state energiesfor the exact case compared to the eo-SCRPA solution.One should notice that there is barely an improvementusing the eo-SCRPA versus the standard SCRPA becausethe latter produced already excellent results. So, we donot show the old SCRPA results again. It is not quiteclear why there is this difference between the Lipkin andHubbard models. Probably the fact that in Lipkin model,contrary to the Hubbard model, one uses collective ph op-erators makes it more difficult to fullfill the Pauli prin-ciple. So the performance of one or the other approach Ω /t sRPA eo-SCRPA Exact |q|= π/3 Ω /t sRPAeo-SCRPAExact |q|=2 π/3 Ω /t sRPAeo-SCRPA Exact |q|= π FIG. 10: Same as Fig.9 but for the energies of excited statesfor different channels | q | = π/ 3, 2 π/ π . seems to depend on the situation. E G S /t sRPAeo-SCRPAExact FIG. 11: Same as Fig.9 but for the ground state energy. IV. CONCLUSION In this work, we coupled even and odd particle num-ber RPA self-consistently. Both systems are based onthe same correlated RPA ground state. From the oddsystem, we get the occupation numbers, odd particle ex-citation energies, and the ground state energies whereasfrom the even SCRPA equations we get the excitationenergies of the even system and transition probabilities.To make things clear, we should mention again that theSCRPA employed here has the same mathematical struc-ture as the one used before [8], only the single particleoccupation probabilities are now calculated via the oddselfconsistent RPA whereas they were obtained beforevia the so-called Catara method [7]. Both even and oddsystems are coupled through non-linear equations whichboth contain the RPA amplitudes X, Y and the s.p. oc-cupation numbers n k in a non-linear way. We called thissystem of equations ’even-odd’ SCRPA (eo-SCRPA). Ap-plications to the Lipkin model and a six sites Hubbardring at half filling gave very satisfying results for all quan-tities. The equations are relatively complex due to theirnon-linearity but they should be solvable with moderncomputers for realistic problems such as, e.g., the cal-culation of collective states in nuclei. The equations tobe solved seem not to be of higher numerical complexitythan, e.g., the Brueckner Hartree-Fock equations whichhave been solved a number of times for nuclei. The cou-pling of even and odd RPA’s has a couple of advantages:it gives richer results, i.e., excitation energies of even andodd particle number systems; there is a natural way howto obtain the ground state energy via the s.p. Green’sfunction and, last but not least, the results seem to bepromising.0 V. ACKNOWLEDGEMENTS We are grateful for long-standing collaboration onSCRPA with D. Delion, J. Dukelsky, and M. Tohyama. Appendix A: Equation of Motion for odd particlenumber operator for Lipkin Model We consider the odd excitation operator as q † µ = 1 N X m x µ m c m + U µ m J + c † m (A1)and the coefficients will be determined from minimisa-tion of expression (5). Based on the solution of theSCRPA equations with the definition of the pair exci-tation operator as Q † = ( XJ + − Y J − ) / p h− J i (with J + = p h− J i ( XQ † − Y Q )), the X, Y amplitudes arethe solutions of the SCRPA equations with H of the Lip-kin Hamiltonian (28). From the minimisation of (5), weobtain a 2 × n = 1 N X m h{ c m , c † m }i = 1 n = n = 1 N X m h{ c m , J + c † m }i = 0 n = 1 N X m h{ c m J − , J + c † m }i = − N ( N − (cid:0) Y (cid:1) h J i + 2 N h J J i (A2)where we have used the inversion (17) and the killingcondition Q | i = 0. For the first Hamiltonian elementwe have H = 1 N X m h{ c m , h H, c † m i }i = − e H = H = 1 N X m h{ c m J − , [ H, c † m ] }i = − e N X m h{ c m J − , c † m }i− VN X m h{ c m J − , J + c † m }i = − V n (A4)And the anti-commutator for H is given by H = 1 N X m h{ c m J − , [ H, J + c † m ] }i (A5)= 3 e n + V (2 − N )[ h J − J − i + h J J − J − i ]= 3 e n − V XY (2 − N )[(1 + 2 Y ) h J i + h J i ] with, for example h J J − J − i = − XY h J i − XY h J J i (A6)where we have again used the inversion (17) and thekilling condition Q | i = 0.The correlation functions which contain quadraticforms of occupation number operators as h J J i in aboveequation can in principle be expressed by the RPA ampli-tudes as well but leading to heavier expressions. Usually,we, therefore will employ the factorization approximationleading in the present case to h J J i ≃ h J i . However,in the Lipkin model one also can use the Casimir relation h J J i = N ( N + 2) + 4 h J i − h J + J − i (A7)Then all matrix elements H ij become functions of h J i and the RPA amplitudes X, Y . The eigenvalue prob-lem can therefore be solved leading to a self-consistencyproblem for h J i and the RPA amplitudes which are ob-tained from the SCRPA equations (11) [3]. The occupa-tion numbers are then given by n = N λ − − H /n λ − − λ + and n = N − n (A8)where λ ± are the eigenvalues of the 2 × λ ± = − e β ± p β + V n (A9)with β = e − V XY ( N − − V XY ( N − Y ) h J i n .Thus, h− J i = n − n = 2 n − N (A10) Appendix B: Equation of Motion for HubbardModel