Variational extensions of BCS theory

R. Balian et al./Physics Reports 317 (1999) 251}358

VARIATIONAL EXTENSIONS OF BCS THEORY

R. BALIAN , H. FLOCARD
, M. VED NED RONI
CEA/Saclay, Service de Physique TheH orique, 91191 Gif-sur-Yvette Cedex, France
Groupe de Physique TheH orique, Institut de Physique NucleH aire, 91406 Orsay Cedex, France

AMSTERDAM } LAUSANNE } NEW YORK } OXFORD } SHANNON } TOKYO

251

Physics Reports 317 (1999) 251}358

Variational extensions of BCS theory R. Balian , H. Flocard
*, M. VeH neH roni
CEA/Saclay, Service de Physique The& orique, 91191 Gif-sur-Yvette Cedex, France
Groupe de Physique The& orique, Institut de Physique Nucle& aire, 91406 Orsay Cedex, France Received November 1998, editor: E. BreH zin Contents 1. Introduction: statistical mechanics of "nite systems 1.1. The two types of equilibrium data 1.2. Inference in statistical mechanics 1.3. Two di$culties of mean-"eld approximations 1.4. Broken symmetries in in"nite and "nite systems 1.5. Outline 2. The variational principle 2.1. A general method for the construction of variational principles 2.2. The action-like functional and the stationarity conditions 2.3. General properties of variational approximations 2.4. Cumulants of conserved observables 2.5. Connection with the standard variational principle 2.6. Variational principle for DK (u) alone 2.7. Comments 3. Extended "nite-temperature HFB approximation for characteristic functions 3.1. Generalities and notations 3.2. The variational ansaK tze and the reduced action functional 3.3. The coupled equations 3.4. Formal results 4. Expansion of the extended HFB approximation: #uctuations and correlations 4.1. The HFB approximation recovered

254 254 255 256 257 259 260 261 262 265 269 270 271 272 273 274 275 278 279 281 282

4.2. Fluctuations and correlations of conserved observables 4.3. Fluctuations and correlations of non-conserved observables 4.4. Diagrammatic interpretation 4.5. Comments 5. Projected extension of the thermal HFB approximation 5.1. Generalities and notations 5.2. The variational ansaK tze and the associated functional 5.3. The coupled equations 5.4. Unbroken PK -invariance 6. Projection on even or odd particle number 6.1. Characteristic function in the N-parity projected HFB approximation 6.2. The partition function in the N-parity projected HFB approximation 6.3. The N-parity projected BCS model: generalities 6.4. The N-parity projected BCS model: limiting cases 6.5. The N-parity projected BCS model: a numerical illustration 6.6. Discussion 7. Summary and perspectives Acknowledgements Appendix A. Geometric features of the HFB theory A.1. The reduced Liouville space A.2. Expansion of the HFB energy, entropy and grand potential

* Corresponding author. e-mail: #[email protected]. 0370-1573/99/$ - see front matter  1999 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 0 - 1 5 7 3 ( 9 8 ) 0 0 1 3 4 - 3

283 288 293 297 298 298 301 303 305 310 310 313 316 321 328 339 341 347 348 348 348

R. Balian et al. / Physics Reports 317 (1999) 251}358 A.3. The HFB grand potential and the RPA equation A.4. Riemannian structure of the HFB theory A.5. Lie}Poisson structure of the time-dependent HFB theory

350 351

Appendix B. Liouville formulation of the projected "nite-temperature HFB equations References

253

354 357

352

Abstract A variational principle is devised which optimizes the characteristic function at thermodynamical equilibrium. The Bloch equation is used as a constraint to de"ne the equilibrium state, and the trial quantities are an unnormalized density operator and a Lagrangian multiplier matrix which is akin to an observable. The conditions of stationarity yield for the latter a Bloch-like equation with an imaginary time running backwards. General conditions for the trial spaces are given that warrant the preservation of thermodynamic relations. The connection with the standard minimum principle for thermodynamic potentials is discussed. We apply our variational principle to the derivation of equations which are tailored for (i) the consistent evaluation of #uctuations and correlations and (ii) the restoration through projection of broken symmetries. When the trial spaces are chosen to be of the independent-quasi-particle type, we obtain an extension of the Hartree}Fock}Bogoliubov theory which optimizes the characteristic function. The expansion of the latter in powers of its sources yields for the #uctuations and correlations compact formulae in which the RPA kernel emerges variationally. Variational expressions for thermodynamic quantities or characteristic functions are also obtained with projected trial states, whether an invariance symmetry is broken or not. In particular, the projection on even or odd particle number is worked out for a pairing Hamiltonian, which leads to new equations replacing the BCS ones. Qualitative di!erences between even and odd systems, depending on the temperature ¹, the level density and the strength of the pairing force, are investigated analytically and numerically. When the single-particle level spacing is small compared to the BCS gap D at zero temperature, pairing correlations are e!ective, for both even and odd projected cases, at all temperatures below the BCS critical temperature ¹ . There exists a crossover temperature ¹ such that odd}even e!ects disappear for  " ¹ such that ¹ (¹(¹ . Below ¹ , the free-energy di!erence between odd and even systems decreases "  " quasi-linearly with ¹. The low temperature entropy for odd systems has the Sackur}Tetrode form. When the level spacing is comparable with D, pairing in odd systems is predicted to take place only between two critical temperatures, thus exhibiting a reentrance e!ect.  1999 Elsevier Science B.V. All rights reserved. PACS: 05.30.Fk; 21.60.-n; 74.20.Fg; 74.80.Bj; 74.25.Bt Keywords: Variational principles; Fluctuations; Broken symmetry restoration; Projection; Even}odd e!ects; Heavy nuclei; Superconducting grains

254

R. Balian et al. / Physics Reports 317 (1999) 251}358

1. Introduction: statistical mechanics of 5nite systems 1.1. The two types of equilibrium data Let us consider a "nite physical system in thermodynamic equilibrium. We want to evaluate variationally the thermodynamic functions of this system along with the expectation values, #uctuations and correlations of some set of observables QK . The equilibrium state of A the system is de"ned from our knowledge of conserved quantities such as the number of particles, the energy, the total angular momentum, the parity. However, because our system of interest is "nite, these equilibrium data should be handled with care. Indeed, depending on the circumstances, their values can be given in two di!erent ways, either with certainty or statistically, and this allows us to classify the conserved quantities into two types. Take for instance a "nite, closed system thermalized with a heat bath. Among the variables that determine its state, the particle number N is given with certainty; we shall call N a conserved quantity of type I. On the other hand, the energy is known only statistically as only its expectation value is "xed; let us call it a conserved quantity of type II. The Gibbs argument, which identi"es a heat bath with a large set of copies of the system, entails that this system is in a canonical Boltzmann}Gibbs equilibrium characterized by a Lagrangian multiplier b de"ning the inverse temperature. Apart from N, our information is speci"ed by b or equivalently by the expectation value of the Hamiltonian. Take now an open system which communicates with heat and particle baths. In this case, both the energy and the particle number are conserved quantities of type II, and the Gibbs argument leads to the grand canonical distribution. For other equilibrium data one has likewise to recognize, according to the circumstances, whether they are of type I or type II. Note that, aside from the truly conserved quantities, the type II equilibrium data may also include nearly conserved quantities (or order parameters) which are here assumed to have all been identi"ed beforehand. For macroscopic systems, this distinction between types I and II is immaterial because the relative #uctuations that occur for observables of type II are negligible in the thermodynamic limit. Canonical and grand canonical distributions then become equivalent. When dealing with xnite systems, however, it is important to specify carefully which statistical ensemble is relevant to the physical situation under consideration. This problem is encountered in various systems of current interest. In nuclear physics, depending on the experimental conditions and/or the accuracy of the description, we may have to describe nuclei having either a well-de"ned or a #uctuating particle number, angular momentum or energy. Similar situations may occur for atoms or molecules. Another example is that of metallic or atomic clusters which can be prepared with a given particle number [1]. Let us also mention mesoscopic rings for which a proper description of the conductivity properties requires that there be no #uctuations (i.e., a canonical distribution) in the number of electrons [2]. Likewise, the description of mesoscopic superconducting islands [3] or metallic grains [4}6] demands a distribution with a given parity of the number of electrons. The recent discovery of Bose}Einstein condensates in dilute alkali atomic gases is currently arousing great interest [7]; in this case also, the theoretical formulation requires us to use states with well-de"ned numbers of particles.

R. Balian et al. / Physics Reports 317 (1999) 251}358

255

1.2. Inference in statistical mechanics A standard formal method to assign a density operator DK to a system characterized by data of both types proceeds as follows [8]. The type I data, referring to some conserved observables, determine a Hilbert space, or a subspace H (for instance the subspace with a given number of  particles) in which the density matrix DK operates. On the other hand, type II data are the expectation values 1XK 2 of the remaining conserved observables which act in H . As a conseI  quence, they impose on DK the constraints 1XK 2"Tr XK DK /Tr DK , (1.1) I  I  where the trace Tr is meant to be taken on the Hilbert subspace H . As we shall see, it is   convenient to leave DK unnormalized or, equivalently, not to include the unit operator among the set XK . Within the subspace H the state DK is determined, up to a normalization factor, by I  maximizing the von Neumann entropy Tr DK ln DK #ln Tr DK , (1.2) S,!   Tr DK  under the constraints (1.1). Associating a Lagrangian multiplier j with each datum 1XK 2, one I I recovers the generalized Boltzmann}Gibbs distribution






(1.3) DK Jexp ! j XK , I I I where j and 1XK 2 are related to each other through I I R 1XK 2"! ln Tr exp ! j XK . (1.4) I  I I Rj I I In what follows, we take the form (1.3) of the density operator for granted, irrespective of the mechanism that led to it. Equilibrium distributions of the form (1.3) can be the outcome of some transport equation, the detailed nature of which is beyond the scope of this work. The form (1.3) can also be justi"ed by various theoretical arguments. One of them is the above-mentioned maximum entropy criterion which has its roots in information theory [9]. Alternatively, the Gibbs argument provides a direct statistical derivation of the Boltzmann}Gibbs distribution. An extension to non-commuting observables XK , relying as in the classical case on Laplace's principle of I indi!erence, has been worked out in Ref. [10]. Formally, the assignment of the density operator (1.3) to the system solves our problem since the thermodynamic functions, expectations values, #uctuations, correlations of the observables and all their cumulants are conveniently obtained from the characteristic function






u(m),ln Tr AK (m)DK , (1.5)  where the operator AK (m), depending on the sources m associated with the observables QK , is de"ned as A A






AK (m),exp ! m QK . A A A

(1.6)

256

R. Balian et al. / Physics Reports 317 (1999) 251}358

Indeed, by expanding Eq. (1.5) in powers of the sources, 1 u(m)"u(0)! m 1QK 2# m m C !2 (1.7) A A A B AB 2 A AB one gets at "rst order the expectation values 1QK 2 of the observables QK . At second order appear A A their correlations C ,1QK QK #QK QK 2!1QK 21QK 2 (1.8) AB  A B B A A B and #uctuations *Q,C , and in higher orders appear the cumulants of higher rank. Moreover, A AA since we chose not to enforce the normalization of DK , the generalized thermodynamic potential is given by the zeroth-order term:






u(0)"ln Tr DK "ln Tr exp ! j XK "S! j 1XK 2 .   I I I I I I

(1.9)

1.3. Two dizculties of mean-xeld approximations Unfortunately, in most physical cases, the calculation above cannot be worked out explicitly. Indeed, the statistical operator (1.3) is tractable only if all the operators XK are of the single-particle I type. For systems of interacting particles, this is no longer the case since, in general, the data 1XK 2 I include the average energy 1HK 2. Recourse to some approximation scheme then becomes unavoidable. Many approaches start with the replacement of the density operator (1.3) by one with the form DK Je\+K ,

(1.10)

where M K is a single-particle operator M aRa , or more generally a single quasi-particle operator GH GH G H which includes also pairs aRaR and a a . (In this article, we shall mainly deal with systems of G H G H fermions; for the boson case, terms linear in the operators aR and a should be added to M K .) G G However, the form (1.10) entails two dizculties, both of which will be adressed in this article. (i) The "rst di$culty concerns the optimum selection of the independent-quasi-particle approximate state (1.10). One often uses a variational criterion, and this is an approach that we shall follow also. [Conventional many-body perturbation theory is not suited to the problems we have in mind since no small parameter is available.] The choice of this criterion, however, raises a crucial question. A standard option, which leads to the Hartree}Fock (or Hartree}Fock}Bogoliubov) solution, consists in choosing the minimization of the free energy or, at zero temperature, of the energy. This procedure is clearly adapted to the evaluation of the thermodynamic quantities that are given by u(0), the value of Eq. (1.5) when the sources m vanish (see Eq. (1.9)). However, the resulting A approximate state is not necessarily suited to an optimization of the characteristic function u(m) for non-zero values of the sources m . For restricted trial spaces such as (1.10) (or (1.11)), a variational A method designed to evaluate u(m), rather than u(0), is expected to yield an optimal independentparticle density operator which, owing to its dependence on the intensity of the sources, is better adapted to our purpose than the (m-independent) Hartree}Fock solution. (ii) The second di$culty is more familiar. Let us for instance consider a system for which the particle number is of type I and the energy of type II. It should be described by a canonical equilibrium distribution DK Jexp(!bHK ) built up in the N-particle Hilbert space H . Since our 

R. Balian et al. / Physics Reports 317 (1999) 251}358

257

system is "nite, the results are expected to di!er from those of a grand canonical equilibrium. However, perturbative and variational approaches are often worked out more conveniently in the entire Fock space H. In particular, even an independent-particle state such as Eq. (1.10) where M K " M aRa resides in this Fock space H. Thus the simple ansatz (1.10) requires leaving the GH GH G H space H to which the exact canonical state belongs; as a consequence, the standard variational  treatment introduces spurious components into the trial state. Unphysical thermodynamic yuctuations of the particle-number thus arise from the form of the trial state (1.10), which would be better suited to a grand canonical than to the canonical equilibrium being approximated. These #uctuations should be eliminated, and it is natural to improve the approximation by introducing a projection onto states with the right number of particles. This is achieved by replacing Eq. (1.10) by the more elaborate ansatz (1.11) DK JPK e\+K PK , , , where PK is the projection onto the subspace H with the exact particle number N. The ,  characteristic function (1.5) can now be evaluated by taking the trace Tr in the full Fock space H. Likewise, when the exact state has a non-#uctuating angular momentum or spatial parity, any trial state of the form (1.10) introduces thermodynamic #uctuations for these quantities, and adequate projections PK are required so as to suppress these #uctuations. However, if the operator M K commutes with PK , it is not necessary in Eq. (1.11) to insert the projection PK on both sides of the exponential. 1.4. Broken symmetries in inxnite and xnite systems Projections turn out to be very useful in another circumstance, namely when a symmetry is broken by the operator M K . Since in this case M K does not commute with PK , the projection should appear on both sides as in Eq. (1.11). Although projection matters mostly for "nite systems, let us "rst recall, for the sake of comparison, the situation in the in"nite systems. Consider a set of interacting particles which, in the thermodynamic limit, presents a spontaneously broken symmetry. One can think of a Bose liquid or a superconducting electron gas where the NK -invariance is broken. In such a system an essential property of the exact statistical operator DK is the long-range order displayed by the associated n-point functions. For an interacting Bose #uid at low temperature, if we denote by tR(r) the creation operator at the point r, expectation values over DK such as 1tR(r )t(r )2, 1tR(r )tR(r )t(r )t(r )2 factorize in terms of a single function u(r) as       u(r )Hu(r ), u(r )Hu(r )Hu(r )u(r ) when the various points all move apart from one another. On the       other hand, since for any size of the system its exact density operator DK commutes with NK , as does the Hamiltonian HK , the quantity 1t(r)2 vanishes. In order to give a meaning to the factors u(r), a standard procedure consists of adding to HK small sources which do not commute with NK : for instance terms in a and aR for bosons, or in aRaR and a a for fermions. The resulting equilibrium G G G H G H density operator DK displays therefore an explicit broken invariance. Despite the fact that such sources do not exist in nature, the remarkable properties of DK make it an almost inevitable substitute for DK . Indeed, in the limit when one "rst lets the size of the system grow to in"nity and then lets the sources tend to zero, DK is equivalent to DK as regards all the observables that commute with NK , such as tR(r )t(r ) or tR(r )tR(r )t(r )t(r ); however, the expectation value 1t(r)2 over       DK does not any longer vanish and is interpreted as an order parameter. Furthermore, DK satis"es

258

R. Balian et al. / Physics Reports 317 (1999) 251}358

the clustering property for all operators, even those which do not commute with NK ; in particular, we now have 1tR(r )t(r )2 P1tR(r )2 1t(r )2 when "r !r "PR so that 1tR(r)2 can be identi      "ed with u(r). The clustering property of DK as well as the equivalence of DK and DK for the physical observables were recognized for the pairing of electrons as early as the inception of the BCS theory [11]. Thus, for su$ciently large systems, the same physical results are obtained not only for di!erent canonical ensembles, but also with or without the inclusion of sources that break invariances. The situation is di!erent for a "nite system governed by the same interactions as in the in"nite case above, under similar conditions of density and temperature. Here again the expectation value 1t(r)2 is zero in the exact equilibrium state, but now no analogue to DK can be devised. Indeed, we can no longer de"ne an order parameter by means of factorization of n-point functions since the relative coordinates cannot grow to in"nity. Nor can we introduce sources which would produce non-vanishing expectation values 1t(r)2 without changing the physical n-point functions, as this requires the sources to tend to zero only after the thermodynamic limit has been taken. Even the translational invariance of a nucleus, an atom or a cluster cannot be broken in the exact ground state or in thermal equilibrium states if the Hamiltonian includes no external forces and commutes with the total momentum. Returning to the NK -invariance, for a nucleus with pairing or for a superconducting mesoscopic island or metallic grain, it is no longer without physical consequences to break the exp(i NK ) gauge invariance. This not only introduces spurious o!-diagonal elements in NK for the density matrix DK , but also a!ects observable matrix elements in the diagonal blocks in NK . Thus, the breaking of invariances does not have the same status for "nite and in"nite systems. In the latter case, it stands as a conceptual ingredient of the theoretical description which correctly implements the clustering property of the n-point functions. For "nite systems it is also essential, but now as a basic tool for building approximation schemes which economically account for some main features of the correlations. Both in perturbative and variational treatments, the starting point is most often an approximate density operator of the form (1.10) which does not commute with conserved observables such as NK . This yields tractable and sometimes remarkably accurate approximations, for which the clustering of the correlation functions arises naturally from Wick's theorem. Remember the success of the BCS theory in describing pairing correlations between nucleons in "nite nuclei. (One may wonder why nuclear physicists, who had acknowledged the existence of pairing e!ects long before 1957, did not anticipate the BCS theory; it may be that, too obsessed by the conservation of the particle number, they were inhibited in developing models which would break the NK -invariance.) Nonetheless, more realistic descriptions of "nite systems require the restoration of the broken invariances. A state of the form (1.10) involves #uctuations in the particle number that must be eliminated, whether they arise from a grand canonical description, or from unphysical o!-diagonal elements of NK connecting Hilbert spaces with di!erent particle numbers. A natural procedure to account for pairing correlations and at the same time enforce the constraint on the particle number consists in using a projection. The use of Eq. (1.11) arises naturally in this context. Likewise, breaking translational or rotational invariance introduces spurious #uctuations of the linear or angular momentum, as well as spurious o!-diagonal elements in the trial density operator (1.10). To remedy these defects we need again a trial density operator of the form (1.11) involving the appropriate projection PK on both sides.

R. Balian et al. / Physics Reports 317 (1999) 251}358

259

The discussion above pertains to the breaking of a conserved quantity of type I, that is, one which is known with certainty , such as NK in a canonical ensemble. For a quantity of type II, known only on average, such as NK in a grand canonical ensemble, the invariance breaking again introduces o!-diagonal elements of NK , which are spurious for a "nite system. One can suppress them, while letting NK free to #uctuate, by means of a di!erent type of projection, exhibited below in Eq. (5.11). The problem of the restoration of broken symmetries has already a long history. Within the nuclear context, in the zero-temperature limit, a general review of the litterature (up to 1980) may be found in [12]. An early reference on parity projection at "nite temperature is [13]. More recent works include Refs. [14}17], where additional references can be found. For a Monte-Carlo approach to the temperature dependence of pair correlations in nuclei, see Refs. [18,19]. In condensed matter physics, recent experimental studies about isolated mesoscopic metallic islands [3,20] or small metallic grains [4}6] have motivated several theoretical works about modi"ed BCS theories with well de"ned number parity [21}28].

1.5. Outline The two problems (i) and (ii) that we have stated in Section 1.3 will be worked out within a general variational setting described in Section 2. We shall write an action-like functional, Eq. (2.5), which yields as its stationary value the characteristic function (1.5). This corresponds to the static limit of a variational principle proposed elsewhere to evaluate multi-time correlation functions [29]. The approach uses a method [30,31], recalled in Section 2.1, which allows the construction of variational principles (and variational approximations) tailored to the optimal evaluation of the quantity of interest (here the characteristic function). This systematic and very general method encompasses several known variational principles, including the Lippmann}Schwinger formalism [32] (the quantity of interest being the S-matrix). An important ingredient is the incorporation in the trial functional of the Lagrangian multipliers associated with the primary equations, in our case the Bloch equation determining the canonical density operator (1.3). As a result, the number of variational degrees of freedom is doubled. Here, the variational quantities entering our functional (2.5) consists of two matrices: one akin to a density operator and the other to the observable AK (m) de"ned in Eq. (1.6). Within this framework, one can state some general conditions on the trial spaces that ensure the preservation of the thermodynamic identities by the variational approximations. Starting in Section 3, we apply the variational principle of Section 2 to systems of fermions for which pairing e!ects are signi"cant. Section 3 (as well as Section 4) deals with grand canonical equilibrium. The two variational operators introduced in Section 2 are taken as exponentials (3.4) of single-quasi-particle operators and the resulting expression of the approximate action-like functional is given in Section 3.2. The equations which express the stationarity and determine the optimal characteristic function, are derived in Section 3.3. They couple, with mixed boundary conditions, the variational parameters characterizing the two trial matrices and they generalize the Hartree}Fock}Bogoliubov (HFB) approximation [33]. In Section 4, the coupled equations obtained in Section 3 are expanded in powers of the sources m . The HFB approximation is recovered at the zeroth order, while the "rst order gives back the A

260

R. Balian et al. / Physics Reports 317 (1999) 251}358

usual HFB formula for the average values of observables (Section 4.1). Fluctuations and correlations of observables are provided by the second order. In the ensuing formulas for conserved and non-conserved observables, the RPA (random-phase approximation) kernel plays a central ro( le. In Section 4.4, through a diagrammatic analysis, we show the correspondence between our results and the summation of the RPA (bubble) diagrams. Section 5.1 addresses the projection problem stated in (ii) in Section 1.3. The functional and the coupled equations that result from variational ansaK tze of the type (1.11) are derived in Sections 5.2 and 5.3, respectively. Special attention is devoted to the case where the variational approximation does not break the PK -invariance (Section 5.4). Section 6 specializes the formalism of Section 5 to the projection on even or odd particle number. The more detailed study of this relatively simple case is motivated by its relevance to the description of heavy nuclei and to the interpretation of recent experiments on isolated superconducting islands [3] or small metallic grains [5]. Equations are obtained in Section 6.3 which should replace the BCS ones. Ensuing results are compared to those of Ref. [22]. In Section 6.4 we analyze some limiting situations with a special emphasis on the low temperature limit. Numerical applications are presented in Section 6.5; they illustrate physical situations encountered in nuclei, metallic islands or aluminium grains. In the last section we summarize our main results and present some suggestions for other applications and for extensions. Appendix A introduces our notation for the Liouville representation of the reduced fermion quasi-particle space. In particular, it gives a proof of the factorization of the RPA kernel (used in Section 4) into the stability matrix and the matrix whose elements are the constants associated with the Lie}Poisson structure of the HFB equations. In Appendix B, using the Liouville-space notation, we reexpress the equations derived in Section 5 in a fashion which exhibits their formal similarity with those obtained in Section 3. The length of this paper is partly a consequence of our wish to allow several reading options. Section 2 presents and discusses the variational principle for characteristic functions at thermodynamic equilibrium which underlies all the other sections. It can be read for itself since its scope is more general than the applications to fermion systems explored in the following Sections. The contents of Sections 5 and 6, which deal with projections, are largely independent of that of Section 4 which focusses on correlations. The reader primarily interested by the variational projection method can therefore bypass Section 4, except for the brief Section 4.1. If he is mostly concerned by the projection on even or odd particle-number, and ready to admit a few formulas demonstrated in Section 5, he can directly jump to Section 6. If his interest is focused upon the improvements over BCS introduced by parity-number projection, he may even begin with Section 6.3. Alternatively, the reader mainly interested by the variational evaluation of correlations and #uctuations in the grand canonical formalism can skip Sections 5 and 6. 2. The variational principle Our purpose in this section is to write a variational expression adapted to the evaluation of the characteristic function u(m) de"ned in Eq. (1.5), exp u(m),Tr AK (m)DK . 

(2.1)

R. Balian et al. / Physics Reports 317 (1999) 251}358

261

The operator AK (m), de"ned in Eq. (1.6), involves the observables of interest QK and depends on the A associated sources m . The density operator DK describes thermodynamic equilibrium. It has matrix A elements in the Hilbert space H de"ned by the conserved quantities of type I. Therefore, the trace  Tr in Eq. (2.1) is taken on this Hilbert space H . According to the Gibbs argument, the operator   DK is an exponential of the conserved quantities XK of type II that are given on average (Eq. (1.3)). I From now on we assume that the Hamiltonian HK is one of these conserved quantities, with b as its conjugate variable. We thus write DK in the form DK "e\@)K ,

bKK , j XK . I I I

(2.2)

We have not introduced the unit operator among the set +XK ,. Then u(0) does not vanish but its I exponential is identi"ed (Eq. (1.9)) with the partition function Tr exp(!bKK ), for instance the  canonical partition function for KK ,HK , or the grand canonical partition function for KK ,HK !kNK . This will allow us to evaluate variationally not only the expectation values and the cumulants of the set +QK ,, but also the thermodynamic quantities. A According to the circumstances, the space H can be the full Fock space, as in Sections 3 and 4,  or a more restricted one such as the N-particle space for canonical equilibrium. Sections 5 and 6 resort in principle to the latter case; nonetheless, the introduction, as in Eq. (1.11), of a projection onto H makes it possible to use again the trace Tr in the full Fock space.  2.1. A general method for the construction of variational principles We shall rely on a methodical procedure for devising variational principles, and hence variational approximations, aimed at the evaluation of some quantity of interest. We brie#y recall the method; for a more detailed exposition see, for instance, Refs. [30,31]. Given a function f +xG, of a set of variables +xG,, we wish to estimate the value f +xG ,, the quantity  of interest, at a point +xG , which is in principle completely determined by the set of equations  g +xG,"0. In practice, we are concerned by the (common) situation where either the equations H g +xG,"0 are untractable and cannot a!ord the exact solution +xG , and/or where f +xG, cannot be H  calculated at the point +xG ,.  It is natural to associate with each equation g +xG,"0 a variable - H akin to a Lagrangian H multiplier, although we do not look here for a stationary value of f under constraints since our constraints are in su$cient number to determine f +xG ,. Then, a variational procedure suited to an  approximate evaluation of f +xG , rests upon the introduction of the functional  I+xG, - H,,f +xG,! - H g +xG, , H H

(2.3)

which involves both the variables xG and the additional quantities - H. The stationary value I of  the functional I+xG, - H, is obviously equal to the exact value f +xG , for unrestricted variations of the  point +xG, - H,. In fact, in this unrestricted case, letting I+xG, - H, be stationary with respect to the set +- H, ensures that g "0 and hence that xG"xG , so that I+xG , - H,"f +xG ,; it is therefore unnecessH     ary to write the stationarity conditions with respect to the set +xG,.

262

R. Balian et al. / Physics Reports 317 (1999) 251}358

With respect to all the variables +xG,- H,, the stationary point +xG , - H, is a saddle point (Fig. 1).   When we restrict the trial spaces of these variables, the stationarity of I provides a point +xG, - H, which can be close to the exact saddle point +xG , - H, if these spaces have been judiciously chosen.   Then, since the surface I+xG, - H, is yat near its saddle point, the di!erence between f +xG , and the  stationary value of I in the restricted space will be much smaller than the error on +xG , - H,.   Several known variational principles enter the framework above [30,31]. Let us, for instance, recall how the Lippmann}Schwinger formulation can be retrieved. The quantity f +xG , of interest is  here the transition amplitude 1W"U(t )2 from the state "U2 at the initial time t to the state "W2 at   the "nal time t . The variables +xG, are the components of the ket "U(t)2 for t (t4t (the function    f +xG , depends only on the part t"t of them). The constraints g +xG,"0 that determine "U(t)2 are   H the components of the SchroK dinger equation (R/Rt#iHK )"U(t)2"0. Denoting as 1W(t)" the associated Lagrangian multipliers - H, we "nd the scattering amplitude as the stationary value of



I+U(t) , W(t),"1W"U(t )2! 

R

R



dt1W(t)



R #iHK U(t)2 , Rt

where "U(t)2 is subject to the initial boundary condition "U(t )2""U2. The additional conditions  RI/RxG"0 provide the SchroK dinger equation for 1W(t)" which must be solved backward in time with the "nal condition 1W(t )""1W". In the original work of Lipmann}Schwinger [32], the  evolution of "U(t)2 is expressed in the interaction representation, and the times t and t are taken   to !R and #R, respectively. When the elimination of one set of variables by means of the exact conditions RI/RxG"0 is feasible, a new variational principle may emerge, under suitable conditions, such that the quantity of interest f +xG , is obtained as a maximum or a minimum, and not merely as a stationary value. This  is visualized in the schematic one-dimensional case of Fig. 1 by the dashed line. We shall encounter such situations in Sections 2.5 and 2.6. 2.2. The action-like functional and the stationarity conditions In our problem, the quantity of interest f +xG , is the characteristic function (2.1), which depends  on the matrix elements of the operator DK "exp(!bKK ). We wish to characterize DK by means of some set of simple constraints that will be identi"ed with the equations g +xG ,"0 of the general H  method. To this aim we remark that a set of linear equations is provided by the Bloch equation. More precisely, we choose to construct the canonical form DK "exp(!bKK ) as the solution of the symmetrized equation 1 dDK # [KK DK #DK KK ]"0 , du 2

(2.4)

where DK (u) is a u-dependent operator (04u4b). This equation is subject to the initial condition DK (0)"1) where 1) denotes the unit operator in the subspace H of the Fock space. The variables  +xG, to be determined are thus the u-dependent matrix elements of DK (u) for 0(u4b, while the constraints g +xG,"0 are identi"ed with the matrix elements of Eq. (2.4). H

R. Balian et al. / Physics Reports 317 (1999) 251}358

263

Fig. 1. A geometric picture of the variational method. The x-axis stands for the variables +xG,, the y-axis for the variables +- H,. The level lines of the functional (2.3) show the existence of the saddle point +xG , - H,, the altitude of which is the   quantity f +xG , of interest. The level lines corresponding to this altitude include the straight line +xG"xG ,, which visualizes   the stationarity of I+xG,- H, with respect to the unrestricted set +- H,. On the other hand, the dashed line visualizes the stationarity of the functional I+xG, - H, with respect to the unrestricted set +xG,, and it meets the line +xG"xG , at the  saddle point +xG , - H,. When +xG, - H, is restricted to some trial set, only part of the surface I+xG, - H, is explored.   Stationarity of I with respect to both +xG, and +- H, in this subset provides a variational approximation in which the error on the altitude of the saddle point is of second order with respect to the deviation of +xG, - H, from the exact point +xG , - H,. Elimination of half the variables by means of RI/RxG"0 leads to a new functional represented by the dashed   line, which may provide f +xG , as a maximum or a minimum under some convexity conditions. 

The exact solution of Eq. (2.4) is exp(!uKK ), the matrix elements of which, for any u, constitute the set +xG ,. Here the quantity of interest f +xG ,"ln Tr AK (m)DK depends only on a subset of the   variables +xG ,, namely on the matrix elements of DK (b). The Lagrangian multipliers - H associated  with the constraints (2.4) are now the matrix elements of a u-dependent operator AK (u), and the functional (2.3) takes the form I+DK (u), AK (u),,Tr AK (m)DK (b)  @ dDK (u) 1 ! du Tr AK (u) # [KK DK (u)#DK (u)KK ] . (2.5)  du 2  Thus, the characteristic function (2.1) is the stationary value of the action-like functional (2.5) under arbitrary variations of the trial operators AK (u) and DK (u), the latter being subject to








264

R. Balian et al. / Physics Reports 317 (1999) 251}358

the boundary condition DK (0)"1) .

(2.6)

The speci"cs of our problem enter the functional (2.5) through the temperature 1/b , the operator KK and through the operator AK (m) (a functional of the observables and their sources). A variational approximation for exp u(m) will be generated by restricting the trial class for AK (u) and DK (u) and by looking for the associated stationary value of Eq. (2.5). The error will be of second order in the di!erence between the optimum trial operators and the exact ones. The doubling of the number of variational degrees of freedom, which include not only the trial state DK (u) but also the Lagrangian multiplier operator AK (u), is only bene"cial when approximations are made; indeed, as already indicated in Section 2.1 for the general procedure, for unrestricted variations of AK (u) the stationarity conditions on AK (u) alone yield the Bloch equation (2.4) and are su$cient to ensure that I+DK , AK , provides the exact value for exp u(m). When DK (u) and AK (u) are restricted to vary within some subset, the stationarity conditions of the functional I+DK , AK , with respect to AK (u) are read directly from Eq. (2.5): Tr dAK 






dDK 1 # [KK DK #DK KK ] "0 . du 2

(2.7)

As regards the variation with respect to DK (u), an integration by parts allows us to rewrite the functional (2.5) in the form I+DK (u), AK (u),"Tr AK (0)#Tr DK (b)(AK (m)!AK (b))   @ dAK (u) 1 # duTr DK (u) ! [AK (u)KK #KK AK (u)] .  2 du  From Eq. (2.8), for 04u(b, one obtains the stationarity conditions:











dAK 1 ! [AK KK #KK AK ] "0 , Tr dDK  du 2



(2.8)

(2.9)

and for u"b: Tr dDK (b)(AK (m)!AK (b))"0 . (2.10)  In Eqs. (2.7) and (2.9) we have omitted the u-dependence of DK , of AK , and of the allowed variations dDK and dAK . By construction of the functional (2.5), Eq. (2.7) gives back the symmetrized Bloch equation (2.4) for unrestricted variations of AK . For unrestricted variations of DK , Eqs. (2.9) and (2.10) yield (in the Hilbert space H ) the Bloch-like equation for AK (u)  dAK 1 ! [AK KK #KK AK ]"0 (2.11) du 2 in the interval 04u(b, and the boundary condition AK (b)"AK (m)

(2.12)

R. Balian et al. / Physics Reports 317 (1999) 251}358

265

at the "nal value u"b. The variable u is thus evolving backward from b towards 0 in Eq. (2.11), the solution of which is (2.13) AK (u)"exp[(u!b)KK ]AK (m)exp[(u!b)KK ] .   The exact equation (2.11) merely duplicates the Bloch equation (2.4). However, in any physical application, the trial spaces for DK (u) and AK (u) have to be restricted to make the calculations tractable. Then, in general, Eqs. (2.7) and (2.9) are coupled and must be solved simultaneously. Moreover, we saw that the boundary conditions (2.6) on DK (0) and (2.10) on AK (b) imply that Eq. (2.7) should be solved starting from u"0 up to u"b while Eq. (2.9) should be solved in the opposite direction. 2.3. General properties of variational approximations The stationary value I of the functional (2.5) depends on an ensemble of parameters g. Among  them are obviously the source parameters m associated with the observables QK [the dependence is A A through the boundary condition (2.10)] and the temperature 1/b. The parameters g also include those, denoted as u, which can occur in the operator KK . In particular an external xeld u, coupled to an observable JK of the system, gives rise to a term !uJK in the Hamiltonian HK . For instance, u may represent a magnetic "eld coupled to the component JK of the magnetic moment operator along the "eld. Moreover, bu may be one of the Lagrangian multipliers j of Eq. (2.2). For instance, for a I grand canonical ensemble where KK "HK !kNK , the chemical potential k is regarded as one of the parameters u; likewise, for a system in uniform rotation, u may be an angular velocity associated with a component JK of the total angular momentum. The symbol g will therefore stand in this subsection for the set +m, b, u,. By construction, our quantity of interest expu(m) is equal to the stationary value I of the  functional (2.5) (`sa stands for stationary); we shall use any one of the three notations ePK"I (g),I +AK (m), b, KK ,,I+DK , AK , g, . (2.14)   E E The second notation speci"es some of the g parameters on which this stationary solution depends (parameters of the u-type are implicit within KK ). The last notation emphasizes that the optimized result depends on the parameters g both directly and through the optimal solution +DK , AK , E E determined within our variational class by the conditions dI+DK , AK , g,/dDK "0 and dI+DK , AK , g,/dAK "0. Let us now state with these notations a general property of variational principles that we shall often use: Owing to the stationarity conditions, the full xrst-order derivative of I with respect to  any g reduces to its partial derivative evaluated at the stationary point:



d R d I (g), I+DK , AK , g," I+DK , AK , g, . E E D dg Rg dg  K DK E AK AK E

(2.15)

The contributions to the derivative occuring through DK and AK cancel due to the stationarity. The E E relation (2.15) holds not only for the exact solution but also within any trial class for AK and DK . Another general property holds whatever the variational space for DK (u) and AK (u). The derivatives dDK /du and dAK /du belong to the set of allowed variations dDK and dAK , respectively. Writing Eqs. (2.7) and (2.9) for these variations and subtracting one equation from the other, we obtain

266

R. Balian et al. / Physics Reports 317 (1999) 251}358

along the stationary path the conservation law: d Tr KK [AK (u)DK (u)#DK (u)AK (u)]"0 . E E E E du 

(2.16)

2.3.1. Characteristic function We shall always use for the operator AK a trial space such that it contains the unit operator and that variations dAK proportional to AK (dAK JAK ) are allowed. We then "nd from Eq. (2.7) that the integrand of the functional (2.5) vanishes, so that the stationary value ePK"I (g)"Tr AK (m)DK (b) (2.17)   E does not involve the integral term. The exact relation (2.1) is thus preserved by the variational approximation. Moreover, since we do not normalize DK , its trial space will allow variations dDK proportional to DK (dDK JDK ). Another helpful property then emerges (which of course holds for the exact solution). Indeed, substituting dAK JAK in Eq. (2.7) and dDK JDK in Eq. (2.9) and combining these equations, we obtain d Tr AK (u)DK (u)"0 . E du  E

(2.18)

Using dDK JDK in Eq. (2.10), we "nd ePK"Tr AK (m)DK (b)"Tr AK (b)DK (b) , (2.19)  E  E E even when AK does not belong to the variational space for AK . The relations (2.18) and (2.19) imply that the approximation (2.17) for exp u is equal to Tr AK (u)DK (u) for any value u, and in particular  E E that ePK"Tr AK (0) .  E

(2.20)

2.3.2. Partition function and expectation value of observables From Eq. (2.17) it follows (when dAK JAK ) that the partition function Tr DK (b), where the  K subscript m"0 is shorthand for g,+m"0, b, u,, is given by the stationary value of the functional (2.5) in the absence of sources (AK "1) ): eP,I +1) , b, KK ,"Tr DK (b) . (2.21)   K Let us now consider arbitrary values of the source parameters m . Using Eq. (2.15) where g stands A for one of these parameters, we "nd d dAK (m) du(m) , I +AK (m), b, KK ,"Tr DK (b) .  E dm  dm dm A A A According to Eqs. (1.7) and (2.21), in the limit where the m's vanish, Eq. (2.22) yields ePK



du Tr QK DK (b) Tr QK DK (b) ! ,1QK 2"  A K "  A K . A dm I +1) , b, KK , Tr DK (b) A K   K

(2.22)

(2.23)

R. Balian et al. / Physics Reports 317 (1999) 251}358

267

This result tells us that the calculation of the expectation value of any observable requires only the solution of Eqs. (2.7) and (2.9) for AK "1) . Eq. (2.23) shows therefore that the density operator DK (b) suited to the calculation of the partition function is also suited to the variational evaluation of K the expectation values of the observables QK . A The knowledge of the m-dependence of DK (b) is only needed for the calculation of the higherE order cumulants. The use of Eq. (2.22) will allow us to gain one order when expanding u in powers of the sources. 2.3.3. The energy To "nd the derivative of I with respect to the inverse temperature b, we rewrite the functional  (2.5) as










du Tr AK DK (u)d(u!b)!AK (u) 





dDI (u) 1 # [KK DK (u)#DK (u)KK ] h(b!u) . du 2

(2.24)  We then apply Eq. (2.15) with g replaced by b. Noting that dDK (u)/du" belongs to the variational S@ class allowed for dDK (b) and using Eq. (2.10), we obtain I"

ePK

d 1 du(m) , I "! Tr KK [AK (b)DK (b)#DK (b)AK (b)] ,  E E E E db 2  db

(2.25)

in which AK (b) can be replaced by AK when the variational space for AK contains AK . When AK "1) , Eq. E (2.25) becomes d Tr KK DK (b) du(0) K "1KK 2,E , "! ln Tr DK (b)"  (2.26) !  K db Tr DK (b) db  K where we have used Eq. (2.23) and introduced the notation E for the equilibrium energy in the canonical case (KK ,HK ), or for the expectation value 1HK !kNK 2 in the grand-canonical case. Thus, in the framework of our variational method (for any trial space allowing dAK JAK ), the derivative with respect to b of the approximate thermodynamic potential (AK (m)"1) ) provides the same result as the expectation value of KK obtained from Eq. (2.23) with QK "KK (i.e. AK (m)"exp !mKK ). 2.3.4. Thermodynamic potential and entropy We can also de"ne an approximate thermodynamic potential F ,!b\u(0) , % and an approximate entropy

(2.27)

du(0) S,b(E!F )"u(0)!b % db Tr KK DK (b) K , "ln Tr DK (b)#b   K Tr DK (b)  K which, owing to Eq. (2.26), satisfy the relation RF S"! % . R¹

(2.28)

(2.29)

268

R. Balian et al. / Physics Reports 317 (1999) 251}358

Our approximation will thus be consistent with the fundamental thermodynamic relations between the partition function, the temperature, the energy and the entropy. However, the approximate entropy (2.28) is not necessarily linked by the von Neumann relation (1.2) to the variational quantity DK (b) which optimizes F . K % Note that the entropy does not directly enter the functional (2.5), in contrast with the usual maximum principle for the thermodynamic potential, which requires the preliminary calculation of the entropy associated with the trial density operator. In fact, within our formalism, the status of the entropy S is that of a subordinate quantity: it can be calculated, using Eq. (2.28), from the b-dependence of the stationary value of the functional (2.5) with AK "1) . This is advantageous since we can thus bypass the explicit calculation of Tr DK ln DK , which is unfeasible when the trial state  DK is, for instance, a projected independent particle state. 2.3.5. External xelds If we use the relation (2.15) where g stands for the external "eld u in a term !uJK entering the operator KK , we have



d 1 @ du du Tr JK [DK (u)AK (u)#AK (u)DK (u)] . (2.30) , I"  E E E E  2 du du  When AK "1) , i.e. when m "0, we shall see in Section 2.5 that under some rather general conditions A (satis"ed for instance in Sections 4.1, 5.4.2, 6.2 and 6.3), the products DK (u)AK (u)"AK (u)DK (u) do E E E E not depend on the variable u. Then Eq. (2.30) reduces to eP

1 d Tr JK DK (b) K "1JK 2 . ln Tr DK (b)"  (2.31)  K b du Tr DK (b)  K In this case, as for Eq. (2.26), our variational approximation is consistent with the thermodynamic identities expressing the expectation values of observables JK as derivatives of the thermodynamic potential. 2.3.6. Evolution versus u of the operators AK DK and DK AK In Sections 3 and 4, AK (u) and DK (u) will be chosen within the trial class of exponentials of quadratic forms of creation and annihilation operators. Apart from the properties used above, namely that dAK and dDK belong to the classes of AK and DK respectively, we note that in this case the quantities dAK AK \, AK \ dAK , dDK DK \ or DK \ dDK span the set of quadratic forms.  Let us more generally consider the case when the trial classes for DK and AK are such that DK \dDK and dAK AK \ belong to a common set of operators. Letting then XK "DK \ dDK "dAK AK \ in Eqs. (2.7) and (2.9) and subtracting one equation from the other, we "nd for any XK within the considered set






(2.32)






(2.33)

Tr XK 

1 dAK DK # [AK DK , KK ] "0 . 2 du

Likewise, if we let >K "AK \ dAK "dDK DK \, we "nd Tr >K 

dDK AK 1 ! [DK AK , KK ] "0 . du 2

R. Balian et al. / Physics Reports 317 (1999) 251}358

269

We shall encounter below these equations (Eqs. (3.52) and (5.41)), with arbitrary quadratic forms of creation and annihilation operators for XK and >K . For unrestricted variations of DK and AK , the operators XK and >K are arbitrary and Eqs. (2.32) and (2.33) imply that AK DK and DK AK both satisfy exact equations of the Liouville-von Neumann type (with a Hamiltonian $iKK /2): 1 dAK DK " [KK , AK DK ], 2 du

dDK AK 1 "! [KK , DK AK ] . du 2

(2.34)

Note however that the boundary conditions (2.6) and(2.12) do not involve DK and AK evaluated at the same argument but DK (0) and AK (b) at the initial and "nal limits of the evolution interval. A large arbitrariness is involved in the solution of the equations (2.34) with these mixed boundary conditions, and they are not su$cient to entail the form exp (!uKK ) for DK (u) and the form (2.13) for AK (u). However the simple structure of Eqs. (2.34) makes them useful. 2.4. Cumulants of conserved observables The variational principle (2.5) provides in general rather complicated expressions for the correlations, and more generally for the higher order cumulants of the observables QK . In this A section we consider a single observable QK which is conserved, that is, which commutes with the operator KK . (The extension to several observables commuting with KK , but not necessarily among themselves, is readily performed by replacing mQK by m QK .) In the exact case, one has the obvious A A A relation (2.35) ePK"Tr e\K/K e\@)K "Tr e\@)K >K/K @ ,   which tells that the exponential of the characteristic function (1.5) associated with the observable QK is equal to the partition function corresponding to the shifted operator KK (m),KK #mQK /b .

(2.36)

We shall now show that the relation (2.35) remains valid in variational approximations such that the trial sets for the operators AK and DK are invariant under left or right multiplication by exp (!jQK ), where j is any c-number and QK the considered conserved observable. This invariance property will be satis"ed in Sections 3 and 4, where AK and DK are taken as exponentials of quadratic forms, provided QK is a conserved single quasi-particle observable. We shall encounter the relation (2.35) in a di!erent context in Sections 5 and 6. Let us indeed consider a trial set +AK (u), DK (u), satisfying the boundary conditions (2.6) and (2.12) with AK (m)"exp (!mQK ). We note that the two operators AK (u), DK (u) de"ned by AK (u)"eSK/K @AK (u)eSK/K @,

DK (u)"e\SK/K @DK (u) e\SK/K @ ,

(2.37)

belong to our trial set and obey the boundary conditions (2.6) and (2.12) with AK "1) . Moreover the value of the functional (2.5) calculated with AK (u), DK (u), AK "1) and KK  is equal to the value calculated with AK (u), DK (u), AK (m)"exp (!mQK ) and KK . By sweeping over the trial sets +AK (u), DK (u), for AK (m)"exp (!mQK ) and +AK (u), DK (u), for AK "1) , we obtain exp u+e\K/K , b, KK ,"I +e\K/K , b, KK ,"I +1) , b, KK #mQK /b, ,  

(2.38)

270

R. Balian et al. / Physics Reports 317 (1999) 251}358

which, within our variational approach, is the counterpart of the relation (2.35). As a consequence, for trial spaces satisfying the invariance above, all the cumulants at temperature 1/b of a conserved observable can be obtained from a calculation of the partition function for the shifted operator KK  de"ned in Eq. (2.36). In particular, the mean value 1QK 2 and the variance *QK are given by





1 R u+1) , b, KK #jQK , "! , b Rj K H R 1 R u+e\K/K , b, KK , u+1) , b, KK #jQK , *QK , " , (2.39) Rm b Rj K H which are consistent with exact thermodynamical relations. An example of special interest concerns the particle-number operator NK in the grand-canonical equilibrium (KK ,HK !kNK ). To obtain the cumulants of NK (i.e. AK "exp (!mNK )) it su$ces to consider the partition function in which k is replaced by k!j with j,m/b. Then Eqs. (2.39) become R 1QK 2,! u+e\K/K , b, KK , Rm



1 R u+1) , b, HK !kNK ,, 1NK 2" b Rk



1 R *NK " 1NK 2 . b Rk

(2.40)

These equations, valid in any variational space that satis"es the aforementioned invariance, coincide with the usual thermodynamic relations for the expectation value and the #uctuations of the particle-number operator. Similar relations arise for any conserved observable such as the momentum or the angular momentum. 2.5. Connection with the standard variational principle While for DK (u) the exact solution exp (!uKK ) has a trivial, exponential dependence on u, the exact solution (2.13) for AK (u) depends on u in a more complicated fashion since ln AK (u) is a linear function of u only if the operators AK and KK commute, as was the case in Section 2.4. For restricted variational spaces, due to the coupling between the di!erential equations (2.7) and (2.9), we expect in general both the approximate solutions DK (u) and AK (u) to display a non-trivial u-dependence. In fact, this feature endows the method with a versatility which partly compensates for the unavoidable restrictions of the trial spaces. Nevertheless, for m "0 (AK "1) ), that is, when we only want a variational evaluation of the A thermodynamic potential u(0), a simpli"cation occurs whenever the trial spaces satisfy the two properties: (i) the trial sets for the operators DK (u) and AK (u) are identical; (ii) if DK belongs to this set, so does cDK ? where c and a are any positive constants. These properties are obviously satis"ed for independent quasi-particle density operators (Sections 3 and 4). They shall also hold in Sections 5.4, 6.2}6.5, where we shall consider projections of these density operators onto subspaces of the Fock space, because the power a of the operator DK is well-de"ned in a basis where it is diagonal, even when DK has vanishing eigenvalues. Thus, within our variational space, we can make the ansatz DK (u)"e\S)K , 

AK (u)"e\@\S)K  , 

(2.41)

R. Balian et al. / Physics Reports 317 (1999) 251}358

271

where the trial operator KK does not depend on u. The trial choice (2.41) for AK (u) is suggested by   the simple form taken by Eq. (2.13) when AK "1) . As a consequence, we have AK (u)DK (u)"DK (u)AK (u)"DK (b)"e\@)K ,ZDK , (2.42)      where we de"ne Z as the normalization Tr e\@)K . The variational expression (2.5) then reduces to  I+KK ,"Z(1!ln Z)!Z Tr DK (ln DK #bKK ) . (2.43)   By maximizing Eq. (2.43) with respect to Z, we obtain the restricted variational quantity ln I+DK ,"ln Z"!Tr DK ln DK !b Tr DK KK ,   where we recognize the von Neumann entropy

(2.44)

S"!Tr DK ln DK (2.45)  associated with the normalized DK . The standard variational principle for thermodynamic potentials is thereby recovered since we now "nd that u(0) is the maximum with respect to DK (under the constraint Tr DK "1) of the right-hand side S!b1KK 2 of Eq. (2.44). This, however, takes place only  for thermodynamic potentials or expectation values (2.23), and for particular trial spaces which satisfy the conditions (i) and (ii). We shall encounter in Sections 5.2 and 5.3 a special case of interest in which our variational principle cannot be reduced to the standard one, even for AK "1) , because the condition (ii) is violated: indeed, since in that case a symmetry is broken, the trial density operator is of the form (1.11) and involves on both sides a projection which does not commute with the operator M K . Here the entropy (2.28) resulting from our variational approximation can be identi"ed with Eq. (2.45). It also coincides with Eq. (1.2) where DK is replaced by DK (b), and it is the Legendre transform  of the approximation for ln Z. One notes that the variational principle (2.5), which involves a saddle-point, has been transformed into the maximum principle (2.43). We already noted in Section 2.1 and Fig. 1 that the elimination of part of the variables xG and - H in the functional (2.3) by means of the conditions RI/RxG"0 may, under suitable conditions, lead to a principle of maximum. Here, the constraints on DK (u) and AK (u) imposed by the ansatz (2.41) also reduce the number of independent variables, though in a di!erent fashion. This restriction on the variational space turns out to be su$cient to ensure that u(0) is obtained as a maximum of I+KK ,.  The analysis of Section 2.4 showed that for an operator QK commuting with KK , the characteristic function u(m) is equal to the thermodynamic potential associated with the shifted operator KK "KK #mQK /b when the trial spaces are invariant under multiplication by exp (!jQK ). In this case the cumulants of any order, and in particular the correlations, can be deduced from the standard variational principle (2.44) with KK replaced by KK . 2.6. Variational principle for DK (u) alone The ansatz (2.41), which restricts the parametrization of the two u-dependent variational quantities DK (u) and AK (u) to a single operator KK , is not suited to the evaluation of the characteristic 

272

R. Balian et al. / Physics Reports 317 (1999) 251}358

function u(m) for m O0 (and therefore of cumulants of non-conserved quantities). In the followA ing Sections we shall therefore exploit the freedom allowed by the full set of the variables DK (u) and AK (u). Nevertheless the form of the exact solution (2.13) for AK (u) suggests another possibility: instead of letting AK (u) and DK (u) vary independently, one can constrain them according to AK (u)"DK (b!u)AK DK (b!u) ,

(2.46)

where AK stands for AK (m). Inserting this trial form into the functional (2.5) leads to the variational expression








@ dDK (u) 1 # [KK DK (u)#DK (u)KK ] , (2.47) du Tr DK (b!u)AK DK (b!u)  du 2  in terms of a single variational quantity DK (u) obeying the boundary condition (2.6). Had we not symmetrized the Bloch equation, we would have found the alternative functional Tr AK DK (b)! 








@ dDK (u) #KK DK (u) , (2.48) du Tr AK DK (b!u)  du  whose stationary value with respect to DK (u) yields the wanted quantity Tr AK DK . A special case of this variational principle has already been proposed and applied in connection with the Monte Carlo method [34]; it is obtained from Eq. (2.48) in the case of canonical equilibrium KK "HK by replacing AK by the dyadic "R21R", where R is a point in the 3N-dimensional space. The stationary value of Eq. (2.48) is then the density matrix itself in the R-representation. If we further specialize DK (u) to depend exponentially on u, as does the exact solution, in the form Tr AK DK (b)! 

DK (u)"e\MK S@ ,

(2.49)

both functionals (2.47) and (2.48) become








@ (2.50) du eMK S@ KK e\MK S@ .  We thus recover the variational expression (3.20) of Ref. [31]. In Eq. (2.50) the trial quantity is the operator MK . For unrestricted variations of MK , the stationarity condition is MK "bKK and the stationary value is the characteristic function Tr AK e\@)K . It has been shown in [31] that, under  rather general conditions, the functional (2.50) is maximum for MK "bKK . Tr AK e\MK 1) #MK ! 

2.7. Comments All subsequent sections are based on the variational functional (2.5) and on the resulting equations (2.7) and (2.9) which express its stationarity in restricted trial spaces. This approach presents several advantages. Within the same variational scheme, the optimization of the characteristic function u(m) provides approximate thermodynamic quantities, average values, correlations or #uctuations and (at least in

R. Balian et al. / Physics Reports 317 (1999) 251}358

273

principle) higher cumulants. We have seen in Sections 2.3 and 2.4 how the results are consistent with general requirements, in particular with thermodynamic relations. We have also seen how the stationarity of Eq. (2.5) entails through Eq. (2.22) cancellations in the expansion of the characteristic function in powers of the sources. This will greatly simplify the forthcoming calculations. We have relied on the Bloch equation to characterize the density operator of the system and to devise the action-like functional (2.5). The simplicity of the Bloch equation is re#ected by this functional, which will allow, for instance, to use for the operator DK (u) trial forms including projections (Sections 5 and 6). Such trial forms are di$cult to implement in the usual variational principle for thermodynamic potentials, since the associated von Neumann entropy cannot be written explicitly. In the present formalism, the entropy can be evaluated a posteriori through Eqs. (2.27) and (2.29). It is only under special conditions that the elimination of AK (u) from the action (2.5) results in the standard variational principle of minimum thermodynamic potential (Section 2.5). The variational principle (2.5) involves a stationarity rather than a maximum or minimum property which, admittedly, would yield a better control. For instance, the existence and uniqueness of the solution for the coupled equations (2.7) and (2.9) with boundary conditions (2.6) and (2.10) are not warranted for an arbitrary choice of the trial spaces. On the other hand, we feel that the internal consistency allowed by the inclusion of sources is an important feature, even if the variational approximation only involves a saddle point. This is exampli"ed by Section 4.5. From now on, we apply the variational expression (2.5) to systems of fermions. However, the formalism of Section 2 is general and, owing to its #exibility, can be adapted to a wider class of problems.

3. Extended 5nite-temperature HFB approximation for characteristic functions We now consider systems of interacting fermions for which pairing e!ects are important. In Section 3.1, after some notations and de"nitions, we recall some properties of operators of the independent-quasi-particle type (Eq. (3.4)) which play a central ro( le throughout this article. In Section 3.2 we take the form (3.4) as the trial class for both the operators DK (u) and AK (u) that enter the functional (2.5). Inserting this choice in (2.5) we derive the explicit form of the reduced action (Eq. (3.32)). In Section 3.3 we write the self-consistent coupled equations (3.36)}(3.37), (3.42)}(3.43) which express the stationarity of the reduced functional (3.32). These equations are suited to the optimization of the characteristic function (2.1); they account for pairing correlations and constitute a generalization of the Hartree}Fock}Bogolyubov equations. In this section (as in Section 4), we deal with grand canonical equilibrium, letting KK ,HK !kNK and taking traces in the full Fock space. The formalism will also apply if KK includes additional constraints of the single-quasi-particle type. For instance, for a nucleus rotating around the z-axis with angular velocity u, the constraint is !uJK where JK is the third component of the angular momentum. For a "nite system, one may also want to control its position in space by including in the set XK the center-of-mass operator, or, if it is deformed, its orientation by including the I quadrupole-moment tensor.

274

R. Balian et al. / Physics Reports 317 (1999) 251}358

3.1. Generalities and notations For the sake of conciseness, it is convenient to introduce 2n-dimensional column and line vectors de"ned as


 a  $

a L , aR"(aR 2aR a 2a ) ,  L  L aR  $

a"

(3.1)

aR L where aR and a are the creation and annihilation operators associated with a "xed single-particle G G basis (i"1, n). If we denote by a the jth component +j"1, 2n, of the column vector a, the H fermionic anticommutation rules are given by either one of the relations [a , aR ] "d , [a , a ] "p , H HY > HHY H H > HHY where p is the 2n;2n matrix




p"



0

1

1

0

(3.2)

.

(3.3)

As in the HFB theory, let us introduce operators with the form ¹K "exp (!l!aRLa) , (3.4)  where l is a c-number and L is a 2n;2n matrix. As a consequence of the anticommutation relation (3.2), we can restrict ourselves to matrices L which satisfy the relation pLp"!L2 ,

(3.5)

where L2 denotes the transposed matrix of L. Instead of l and L one can alternatively use the c-number q and the 2n;2n matrix T de"ned by q"e\J,

T"e\L ,

(3.6)

to specify the operator (3.4). Another equivalent parametrization involves the normalization f, f,Tr ¹K "e\J[det(1#e\L)]"q exp+tr ln(1#T), ,  L and the generalized reduced density matrix R (or matrix of the contractions),







1 Tr a ¹K aR H HY" L R , HHY e #1 Tr ¹K

(3.7)



T " . (3.8) T#1 HHY HHY In Eq. (3.7), the symbol tr denotes a trace in the space of 2n;2n matrices while Tr is the trace in L the Fock space. The property (3.5) implies the relations: pT\p"T2 ,

pRp"1!R2 .

(3.9)

R. Balian et al. / Physics Reports 317 (1999) 251}358

275

As a consequence, the 2n;2n matrix R can be written in the form




R"

o i



i \ , 1!o2

(3.10)

> where the n;n matrices o, i and i correspond to the normal and abnormal contractions: \ > o "Tr a ¹K aR/Tr ¹K , GH G H i "Tr a ¹K a /Tr ¹K , (3.11) \GH G H i "Tr aR¹K aR/Tr ¹K . >GH G H The matrices i and i are antisymmetric. When the operator ¹K is hermitian, one has \ > o"oR, i "iR and the matrix R"RR reduces to the usual HFB form. > \ Operators of the type (3.4) allow the use of a generalized Wick theorem, even when they are not hermitian [35]. They also form a non-abelian multiplicative group, the algebra of which is characterized by the relations (3.20)}(3.24). These simple algebraic properties (more details can be found in [35] and the Appendix of [36]) will be used extensively below.



3.2. The variational ansaK tze and the reduced action functional Let us come back to our problem. We assume that the observables QK are operators of the A single-quasi-particle type (this restriction is not essential and can be removed at the cost of some formal complications (see [29])): QK "q #aRQ a , A A A  so that the operator AK (m) (see Eq. (1.6)) belongs to the class (3.4). We thus write it as AK (m),¹K  ,exp (!l  !aRL  a) , 

(3.12)

(3.13)

with l  , m q , L  , m Q . (3.14) A A A A A A Now we choose the trial spaces for the operators DK (u) and AK (u) entering the functional (2.5). This choice determines our variational approximation. We take as trial sets: DK (u)"¹K B (u) ,

AK (u)"¹K ? (u) ,

(3.15)

with (3.16) ¹K B (u),exp(!l B (u)!aRL B (u)a) ,  ¹K ? (u),exp(!l ? (u)!aRL ? (u)a) . (3.17)  The 2n;2n matrices L B (u) and L ? (u) (constrained to satisfy the relation (3.5)) and the c-numbers l B (u) and l ? (u) constitute a set of independent variational parameters. Using Eqs. (3.6)}(3.8), we de"ne the coe$cients q ? , q B , the normalization factors f ? , f B , the matrices T ? , T B and the reduced density matrices R ? , R B associated with the operators ¹K ? and ¹K B .

276

R. Balian et al. / Physics Reports 317 (1999) 251}358

The boundary condition (2.6) on DK (0) implies the relations: l B (0)"0 ,

L B (0)"0 ,

(3.18)

or equivalently q B (0)"1 ,

T B (0)"1 ;

f B (0)"2L ,

R B (0)" . (3.19)  The products DK (u)AK (u) and A K (u)DK (u) occur in the functionals (2.5) and(2.8). Taking advantage of the group properties [35] of the operators (3.4), we write AK DK "¹K ?B ,

DK AK "¹K B? ,

(3.20)

with ¹K ?B "exp (!l ?B !aRL ?B a) , (3.21)  ¹K B? "exp (!l B? !aRL B? a) , (3.22)  where the u-dependence of the c-numbers l ?B and l B? , and of the matrices L ?B and L B? has not been indicated; using the notation (3.6), one has the relations: e\J ?B "q ?B "e\J B? "q B? "e\J ? \J B "q B q ? ,

(3.23)

e\L ?B "T ?B "T ? T B ,

(3.24)

e\L B? "T B? "T B T ? .

According to Eqs. (3.7) and (3.23)}(3.24), the normalization of DK AK and AK DK , which we denote in shorthand by Z, is given at each `timea u by Z,Tr AK DK "f ?B "Tr DK AK "f B?

"e\J ? >J B  exp +tr ln(1#e\L ? e\L B ),  L "q ? q B exp +tr ln(1#T ? T B ), . (3.25)  L As we already know from Section 2.3, the normalization Z does not depend on u at the stationary point and it is a variational approximation for the quantity exp u that we are looking for (see below Eqs. (3.48)}(3.50)). It thus plays the ro( le of a generalized partition function which reduces to the grand canonical partition function when the sources m are set to zero. A We can now write the functional (2.5) (which provides the characteristic function (2.1)) in terms of the variational parameters q B , T B , q ? , T ? . The term Tr AK dDK /du is obtained by di!erentiating q B and T B in Eq. (3.25) with respect to u, while holding q ? and T ? "xed, which yields Tr AK






d lnq B 1 dT B

dDK "Z # tr T B \R B?

. L du 2 du du

(3.26)

The last two terms in the integrand of the functional (2.5) are supplied by a straightforward application of the generalized Wick theorem. They involve the operator KK ,HK !kNK where the Hamiltonian HK of the system is assumed to have the form L 1 L < aRaRa a , HK " B aRa # GHIJ G H J I GH G H 4 GHIJ GH

(3.27)

R. Balian et al. / Physics Reports 317 (1999) 251}358

277

and where L (3.28) NK " aRa . G G G The matrix elements B of the one-body part of HK and the antisymmetrized matrix elements < of GH GHIJ the two-body interaction satisfy the relations B "BH , < "!< "!< "
Tr KK DK AK "Z E+R B? , ,

(3.30)

with 1 (3.31) E+R," (B !kd ) o # < (2o o !i i ) . GHIJ IG JH >GH \IJ GH GH HG 4 GHIJ GH When o"oR and i "iR , the expression (3.31) is formally identical to the HFB energy. However, > \ the quantities E+R ?B , and E+R B? , are not necessarily real since R ?B and R B? do not have to be hermitian. Altogether, for the trial choices (3.16) and (3.17), the functional (2.5) reads I+AK , DK ,"Z  !



@

du Z

 dlnq B 1 dT B 1 ; # tr T B \R B?

# [E+R B? ,#E+R ?B ,] , du 2 L du 2






(3.32)

where use has been made of the relation T B \R B? "R ?B T B \ .

(3.33)

According to Eq. (3.25), the boundary term Z  is given by Z  ,Tr AK (m)DK (b),Tr ¹K  DK (b)"q  q B exp +tr ln(1#T  T B (b)), , (3.34)  L where the c-number q  and the matrix T  , which specify the characteristic function we are looking for, are expressed in terms of the sources m by Eqs. (3.12)}(3.14) and (3.6). Similarly the A form (2.8) of I reads



I+AK , DK ,"Z  !Z(b)#Z(0)#




;

@



du Z



dlnq ? 1 dT ? 1 # tr T ? \R ?B

! [E+R ?B ,#E+R B? ,] . du 2 L du 2

(3.35)

278

R. Balian et al. / Physics Reports 317 (1999) 251}358

In the functionals (3.32) and (3.35), the u-dependence has been omitted in the variational parameters q B , T B , and q ? , T ? , as well as in the normalization Z given by Eq. (3.25) and in the contractions R B? , R ?B associated with T B? , T ?B , respectively; the latter quantities are related to T ? and T B through Eq. (3.24). 3.3. The coupled equations The next step is to write the equations expressing the stationarity of the functional (3.32) or (3.35). We derive them directly from these expressions, but could also have used the general equations (2.7) and (2.9)}(2.10). Requiring the stationarity of (3.32) with respect to q ? (u) (in the functional (3.32), q ? (u) only appears as a proportionality factor in Z(u)) gives us the equation of evolution for l B (u)"!ln q B (u): dT B 1 dl B 1 " tr T B \R B?

# [E+R B? ,#E+R ?B ,] , 2 L du 2 du

(3.36)

which we shall solve explicitly below (Eq. (3.46)). We now require the stationarity of the functional (3.32) with respect to T ? (u), which appears through R ?B (u), R B? (u) and Z(u). Eq. (3.36) implies that we can disregard in this calculation the variation of Z(u). If we note that the variations of R ?B (u) and R B? (u) resulting from those of T ? (u) satisfy dR B? "T B dR ?B T B \ as a consequence of Eq. (3.33), we readily "nd the di!erential equation dT B

1 "! [H+R B? ,T B #T B H+R ?B ,] , du 2

(3.37)

where H+R, has the same formal de"nition, (3.38) dE+R,"tr H+R, dR ,  L as the HFB Hamiltonian, although H+R, is not necessarily hermitian. For the Hamiltonian (3.27), the 2n;2n matrix H takes the form




H+R,"

h

D



\ , D !h2 >

(3.39)

with h "(B !kd )# < o , GH GH GH GIHJ JI IJ 1 D " < i , \GH 2 GHIJ \IJ IJ 1 D " i < . >GH 2 >IJ IJGH IJ

(3.40)

R. Balian et al. / Physics Reports 317 (1999) 251}358

279

As a consequence of Eq. (3.29), the matrices D et D are antisymmetric. This entails for H the \ > same relation as Eq. (3.5): pHp"!H2 ,

(3.41)

which is consistent with pdRp"!dR2 and Eq. (3.38). Likewise, the stationarity conditions of the functional (3.35) with respect to q B (u) and T B (u) provide, for 04u(b, the di!erential equations dT ? 1 dl ? 1 " tr T ? \R ?B

! [E+R ?B ,#E+R B? ,] , 2 L du 2 du

(3.42)

for l ? (u)"!ln q ? (u) and dT ? 1 " [H+R ?B ,T ? #T ? H+R B? ,] , du 2

(3.43)

for T ? (u). Finally, by rendering the functional (3.35) stationary with respect to q B (b) and T B (b) (see Eqs. (3.13) and (3.6)), one obtains the boundary conditions for AK (u): q ? (b)"q  ,

T ? (b)"T  .

(3.44)

The main equations to be solved are the di!erential equations (3.37) and (3.43) for T B (u) and T ? (u) which are coupled through the matrices H+R ?B , and H+R B? ,. Moreover, one has to deal with the mixed boundary conditions T B (0)"1 , T ? (b)"T  .

(3.45)

All our variational quantities depend upon the sources m (see Eqs. (3.12)}(3.14) and(3.6)) through A the boundary conditions (3.44). One notes that the evolution equation (3.43) of T ? can be deduced from Eq. (3.37) for T B by an overall change in sign of the right hand side and by the exchange a  d, as could have been inferred from the comparison between the forms (3.32) and(3.35) of the functional I. (The same is true for the equations (3.36) and (3.42) governing l B and l ? if one takes account of Eqs. (3.37) and (3.43).) Once Eqs. (3.37) and (3.43) are solved, the c-numbers l B (u) and l ? (u) are obtained by integration of Eqs. (3.36) and (3.42) whose r.h.s. depend only on T B and T ? . Indeed, using Eqs. (3.37) and (3.43) with the boundary conditions l B (0)"0 and l ? (b)"l  , we "nd



l B (u)"

S

du +!tr [R B? H+R B? ,#R ?B H+R ?B ,]#[E+R B? ,#E+R ?B ,], ,   L

 l ? (u)"l  #l B (b)!l B (u) .

(3.46) (3.47)

3.4. Formal results From the equations for T B , T ? , l B , l ? , together with the de"nition (3.25) and the relation (3.38), one can readily verify the conservation laws (2.18) and (2.16) which now

280

R. Balian et al. / Physics Reports 317 (1999) 251}358

read, respectively: dZ "0, du

d (E+R ?B ,#E+R B? ,)"0 . du

(3.48)

One can also verify that the integrands of the functionals (3.32) or (3.35) vanish. As was shown in Section 2, these properties directly result from the feature that variations dAK JAK and dDK JDK are allowed in the trial classes (3.16)}(3.17). We are looking eventually for the characteristic function u+AK (m), b, KK , which is given by the logarithm of the stationary value of the functional (3.32). This value is provided by the boundary term Z  de"ned in Eq. (3.34). In this formula, the matrix T B (b) and the c-number q B (b)"exp [!l B (b)] are furnished by the solution of the coupled Eqs. (3.37), (3.43) and by (3.46), respectively. Since, according to Eq. (3.48), the normalization Z (de"ned in (3.25)) does not depend on u we have the relation u(m),u+AK (m), b, KK ,"ln Z  "ln Z(u) ,

04u4b .

(3.49)

In particular, Z can be evaluated at u"0 (see Eq. (2.20)); from Eqs. (3.25) and (3.19) we "nd ln Z"!l ? (0)#tr ln(1#T ? (0)) , (3.50)  L where l ? (0)"l  #l B (b) is explicitly given by Eqs. (3.46)}(3.47) once the coupled equations (3.37) and (3.43) which yield T ? (0) have been solved. Let us add a few remarks about the properties of these last equations. In the derivations of this section (as well as in Section 2), the hermiticity of the operators DK (u) and AK (u), as de"ned in Eqs. (3.15)}(3.17), was nowhere assumed. If, for some value of u, the operators DK and AK are hermitian, i.e. if T B "T B R, T ? "T ? R, q B "q B H and q ? "q ? H, it follows that R ?B "R B? R and hence that H+R ?B ,"H+R B? ,R. When the operator AK (m) is hermitian, the inspection of Eqs. (3.37), (3.43) and (3.45)}(3.47) shows that these equations preserve the hermiticity of T B and T ? and the reality of q B and q ? . In this case, the stationary value Z given by (3.50) is also real. If AK (m) is not hermitian, as would happen for the characteristic functions with imaginary m's, the boundary condition (3.44) prohibits AK (u) and DK (u) from being hermitian. Combining Eqs. (3.37) and (3.43), one obtains the evolution equations for the products T ?B "T ? T B and T B? "T B T ? , and hence for the HFB-like contraction matrices R ?B and R B? associated (through Eq. (3.8)) with the operators ¹K ?B and ¹K B? . The calculation relies on the property dR"(1!R)dT(1!R) ,

(3.51)

and it leads to dR ?B 1 dR B?

1 " [H+R ?B ,, R ?B ] , "! [H+R B? ,, R B? ] . du 2 du 2

(3.52)

These equations, the single-quasi-particle reductions of Eqs. (2.34), involve commutators. They are reminiscent of the time-dependent Hartree}Fock}Bogoliubov equations (see Section A.5 of Appendix A), $i u replacing the time t. However, here H+R ?B , and H+R B? , are in general not hermitian and u is real. A consequence of the simple form (3.52) is that the eigenvalues of R ?B and

R. Balian et al. / Physics Reports 317 (1999) 251}358

281

R B? remain constant with u, in contrast to those of R ? and R B . It also follows from Eqs. (3.52) and(3.38) that each of the quantities E+R ?B , and E+R B? , is conserved, a "ner result than Eq. (3.48). In spite of the form of Eqs. (3.52), one should not infer from them that the evolutions of the matrices R ?B and R B? are decoupled. Actually the boundary conditions (3.19) and(3.44) on R B (0) and R ? (b) do not entail any explicit boundary conditions for R ?B and R B? . Moreover, Eqs. (3.52) are not equivalent to Eqs. (3.37) and (3.43) since R ? and R B are not determined uniquely by the knowledge of R ?B and R B? .

4. Expansion of the extended HFB approximation: 6uctuations and correlations In this section, starting from the characteristic function (3.49) as expressed by either one of the formulae (3.34) or (3.50), we derive the thermodynamic quantities, the expectation values of the observables, and their #uctuations and correlations. In principle this demands the determination of u(m) and its expansion up to second order in the sources. However, we shall see that the full solution of the coupled equations (3.37), (3.43), (3.46) and (3.47) for "nite values of the sources m is not A required. Indeed all the quantities of interest can be obtained through a prior expansion of the coupled equations in powers of the sources m , followed by a solution of the resulting equations A order by order. Moreover, owing to the variational nature of the approach, we shall need to expand the coupled equations to "rst order only. This procedure turns out to overcome the di$culties associated with the mixed boundary conditions. As expected, the usual HFB equations are recovered in the case AK "1) where the sources m vanish (Section 4.1.1). In agreement with Eq. (2.23) these equations are not only relevant for the A thermodynamic potential and for the ensuing quantities, but also for the calculation of the average values of single-quasi-particle operators (Section 4.1.2). In Sections 4.2 and 4.3 we work out the next order in the sources; this yields an explicit variational expression for the correlations C between two single-quasi-particle observables AB QK and QK (and for the #uctuations when c"d). Two di!erent situations are discussed. The "rst one A B corresponds to the case where the observables QK or QK commute with the operator KK . One then A B recovers the formula (4.27) which involves the HFB stability matrix F . When QK and QK do not  A B commute with KK , one is led to the more general formula (4.63) which involves both the stability matrix F and the RPA kernel K . In Sections 4.2.3 and 4.3.3 we discuss the e!ect of the vanishing   of some eigenvalues of F and of K that may in particular occur when an invariance is broken. We   shall also evaluate the Kubo correlations, Eq. (4.37), and give a diagrammatic interpretation (Section 4.4) of the results obtained here variationally. Although the trial operator DK describes independent quasi-particles, its dependence on the sources yields expressions for the cumulants of order two and beyond which di!er from what might have been naively inferred from the standard mean-"eld HFB approximation. Indeed, the latter is only suited to the evaluation of expectation values; in particular the approximations that it provides for correlations and #uctuations violate the usual thermodynamic relations. This defect is cured in the present fully variational approach, as illustrated in Section 4.2.1 by the example of the #uctuations of the particle number in the grand canonical ensemble. In the present section, the increasing orders of the expansion with respect to the sources m will be A indicated by lower indices (0, 1,2) assigned to the relevant quantities.

282

R. Balian et al. / Physics Reports 317 (1999) 251}358

4.1. The HFB approximation recovered 4.1.1. Thermodynamic quantities To zeroth order in the sources (i.e., AK (m)"1) ) our problem amounts to evaluating the partition function Z"Tr DK (b). We expect that the coupled equations of Section 3.3 will yield the standard  HFB approximation for Z, since the trial spaces for DK (u) and AK (u) satisfy the conditions required in the discussion of Section 2.5. Let us verify this point. If we associate with the contraction matrix R the single-particle Hamiltonian H de"ned by   1 1!R , R , H , (4.1) bH ,ln  e@ #1  R  the usual optimum HFB matrix R is given, at temperature 1/b, by the self-consistent equation  H "H+R , , (4.2)   where H+R, was de"ned in Eqs. (3.39)}(3.40). Because the relation (4.2) entails the commutation of H+R , with R , Eqs. (3.52) are satis"ed by the u-independent solution   R ?B (u)"R B? (u)"R . (4.3)    In the equations (3.37) and (3.43) for T B (u) and T ? (u), the e!ective Hamiltonians do not depend   on u since H+R ?B ,"H+R B? ,"H . Moreover, the boundary conditions (3.45) reduce to    T B (0)"1, T ? (b)"T  "1. We thus "nd    H  (4.4) T B (u)"e\S , T ? (u)"e\@\SH .   Eqs. (4.3)}(4.4) are the single-quasi-particle reductions of the general relations (2.42) and (2.41). One also notes that R is equal to R B (b), the contraction matrix associated with the state DK (b).    Using the solution (4.4) of the coupled equations, we can now calculate l ? (0) from Eqs. (3.46)}(3.47), and hence ln Z from Eq. (3.50). The integrations over u, which can be performed explicitly, yield ln Z from which, according to Eqs. (2.27) and (2.28), we de"ne the thermodynamic grand potential F and the entropy S+R , %  ln Z"S+R ,!bE+R ,,!bF +R , . (4.5)   %  The entropy is given by (4.6) S+R,"!tr [R ln R#(1!R)ln(1!R)] .  L It is thus the entropy of a quasi-particle gas characterized by the contraction matrix R. Owing to the symmetry relation (3.9), the two terms under the trace in Eq. (4.6) are equal. The HFB energy E+R ,"1HK !kNK 2 is given by the expression (3.31).  As expected from the discussion of Section 2.5, Eq. (4.5) coincides with the outcome of the variational principle (2.44), which reads here (4.7) ln Z"ln Tr DK (b)"MaxR[S+R,!bE+R,]"!b Min F +R, . %  The stationarity condition of Eq. (4.7) is obtained from Eq. (3.38) and from the "rst-order variation 1!R 1 dR dS+R," tr ln R 2 L

(4.8)

R. Balian et al. / Physics Reports 317 (1999) 251}358

283

of the reduced entropy (4.6); it thus gives back the usual self-consistency condition (4.2). Therefore, at zeroth order in the sources, we have recovered the standard self-consistent HFB approximation for the thermodynamic potential. 4.1.2. Expectation values of observables The variational approximation for the expectation values 1QK 2 in grand canonical equilibrium is A in principle given by the "rst-order terms in the sources m . However, the knowledge of the A zeroth-order solution (4.1),(4.2) is here su$cient since, according to Eq. (2.23) and to the de"nitions (3.12) and (3.8), we have 1 Tr QK DK (b) A  "q # tr Q R . (4.9) 1QK 2" A 2 L A  A Tr DK (b)  For the expectation value of the operator KK , we can apply Eq. (2.26) which was obtained by derivation with respect to the parameter b. Using Eqs. (3.30) at zeroth order and (4.3), we "nd R 1KK 2"E+R ,"! ln Z ,  Rb

(4.10)

where, as in Eq. (4.9), the expectation value is taken over DK (b).  For the calculation of 1NK 2, we can use either (4.9) or the discussion in Sections 2.3 and 2.4. From Eq. (2.31) (or (2.40)) we get 1 R n 1 ln Z , 1NK 2"tr o" # tr NR " L  L b Rk 2 2

(4.11)

where, as in Eq. (3.12), the operator NK can be represented by the c-number n/2 and the 2n;2n matrix




N"

1

0

0 !1

.

(4.12)

The formulas (4.9)}(4.11) are those given by the standard approximation, which consists in evaluating expectation values in the HFB state DK (b). Our variational procedure thus does not  yield anything new when the characteristic function is expanded to "rst order in the sources. This is not the case at the next order. 4.2. Fluctuations and correlations of conserved observables 4.2.1. The particle number operator As a "rst example of a conserved observable (i.e., which commutes with KK ), we consider the particle number NK and we evaluate its variance *NK . A naive method would consist in taking the di!erence between the expectation value 1NK 2 and the square 1NK 2, both evaluated by means of the Wick's theorem applied to the state DK (b). This would yield  1 (4.13) 1NK 2!1NK 2" tr N(1!R )NR ,   2 L where R is the contraction matrix (4.1)}(4.2) corresponding to DK (b).  

284

R. Balian et al. / Physics Reports 317 (1999) 251}358

However, the state DK (b) is variationally "tted to the evaluation of the thermodynamic potential  F , not to the evaluation of *NK . In agreement with our general philosophy, we should rather rely % on a variational principle suited to the evaluation of the characteristic function u. The variational approximation for the #uctuation *NK  is then supplied by the second derivative of u with respect to the source associated with the observable NK . The result is expected to di!er from Eq. (4.13), since the optimal values of the trial operators DK (u) and AK (u) are thus dependent on the sources m . A Actually, under conditions on the trial classes that are obviously satis"ed here, we have already determined in Section 2.4 the result of this evaluation: 1 RF 1 R1NK 2 %, "! *NK " b Rk b Rk

(4.14)

where F and 1NK 2 are the outcomes of the variational calculation at zeroth and "rst orders in the % sources, given by Eqs. (4.5) and (4.11), respectively. The identities (4.14), satis"ed by our approximation for *NK , agree with exact thermodynamic relations which are violated by the naive result (4.13). Inserting Eq. (4.11) into (4.14), we have 1 RR *NK " tr N  . L 2b Rk

(4.15)

The derivative RR /Rk is given by Eq. (4.2) which de"nes R self-consistently; the parameter   k enters this equation through the contribution !kN to H+R ,. The expression (4.15) reduces to  (4.13) only for a system of non-interacting particles, with [N, R ]"0, and is otherwise non trivial  due to self-consistency. We shall write it explicitly below (Eq. (4.27)) in a more general context. The result (4.15) arising from of our variational treatment is more satisfactory than the naive expression (4.13). We noted that it meets the thermodynamic identities (4.14). Moreover, for a system whose pairing correlations do not disappear at zero temperature, Eq. (4.13) has a spurious "nite limit as b goes to R, whereas the variational estimate (4.15) for the particle-number dispersion vanishes in this limit as 1/b, with a coe$cient equal to the density of single-particle states at the Fermi level. It is therefore consistent with the expectation that, for "nite systems with or without pairing, results from the grand-canonical and canonical ensembles must coincide at zero temperature. 4.2.2. Characteristic function for conserved single-quasi-particle observables The ideas above are readily extended to correlations between any pair of single-quasi-particle observables QK and QK commuting with the operator KK . Indeed, it follows from the discussion of A B Section 2.4 that the evaluation of the characteristic function, and therefore of the cumulants, can be reduced in this case to the calculation of a thermodynamic potential when the trial spaces for AK and DK are invariant with respect to the transformation (2.37). Within the spaces (3.16) and (3.17) considered in this section, this invariance is satis"ed for observables QK of the type (3.12). A The identity (2.38) then shows us that the variational evaluation of the characteristic function u(m) amounts to the variational evaluation of the partition function associated with the shifted Hamiltonian KK "KK # m QK /b. The calculation is thus the same as in Section 4.1.1, except for A A A the change of KK into KK , and Eq. (4.7) is replaced by





u(m)"MaxR S+R,!bE+R,! m (q #tr Q R) . A A  L A A

(4.16)

R. Balian et al. / Physics Reports 317 (1999) 251}358

285

The optimum matrix R is characterized by the extremum condition (4.16), which yields the self-consistent HFB equations 1 H+R,,H+R,# m Q . (4.17) A A b A The characteristic function (4.16) depends on the m's both directly and through the matrix R. From its "rst derivatives we recover the expectation values (4.9). The second derivatives, taken at m"0, yield the correlations between the QK 's. In order to write their explicit expression, setting A R"R #dR, we shall expand S+R,, E+R, and F +R, around R up to second order. To this  %  aim, we introduce a condensed notation. Up to now, we have regarded R , H as well as other quantities denoted by a script HHY HHY capital (Q, L, T, dR, H , etc.) as 2n;2n matrices. We "nd it convenient to consider  them alternatively as vectors in the Liouville space, with the pair (jj) playing the ro( le of a single index (see Appendix A, Section A.1). Thus, the "rst-order variations (4.8) of S+R, and Eq. (3.38) of E+R, appear now as scalar products that we shall denote by the colon sign : which stands for tr .  L Likewise the second-order variations will generate matrices in the Liouville space, with two pairs of indices (j j). With this notation we can write, for R"R #dR, the expansions of the HFB entropy S+R,,  energy E+R, and grand potential F +R,, de"ned by Eqs. (4.6), (3.31) and (4.5), respectively, as (see % Appendix, Section A.2) 1 , R" H +R , e@ Y Y #1

S+R,KS+R ,#bH : dR#dR : S : dR ,     where bH stands for ln(1!R )/R as in Eq. (4.1);    E+R,KE+R ,#H+R , : dR#dR : E : dR ,     where H+R , is de"ned by Eq. (3.38);  F +R,KF +R ,#dR : F : dR , % %    where we made use of the stationary condition (4.2) and where 1 F ,E ! S   b 

(4.18)

(4.19)

(4.20)

(4.21)

is the stability matrix around the HFB solution. In Eqs. (4.18)}(4.20), the second derivatives S ,  E and F , as well as the "rst derivatives bH and H+R ,, are functions of the HFB density R . In      the j-basis of Section 3.1, S , E and F are 4-index tensors; we regard them here as (symmetric)    matrices in the Liouville space (jj). The explicit expression of S is given in Eq. (A.10).  HHYIIY Contraction of such a matrix with a vector proceeds as 1 dR . (4.22) (S : dR) , S  HHYIIY IYI  HHY 2 IIY In Appendix A (Section A.4), we endow !S with a geometric interpretation. It appears naturally  as a Riemannian metric tensor which de"nes a distance between two neighbouring states characterized by the contraction matrices R and R #dR.  

286

R. Balian et al. / Physics Reports 317 (1999) 251}358

The correlation C between the two observables QK and QK (de"ned as in Eqs. (1.7)}(1.8)) is AB A B obtained from the expansion of the expression (4.16) around R :  u(m)Kln Z! m 1QK 2!Min R[b dR : F : dR# m Q : dR] . A A B   A A A A The extremum condition,

(4.23)

b F : dR"! m Q ,  A A A determines dR as function of the m's. From Eq. (4.23), one has

(4.24)





Ru R C " "!Q : . dR AB Rm Rm A Rm A B K B K Inserting into Eq. (4.25) the formal solution of Eq. (4.24),



R 1 dR "! F\ : Q , B Rm b  K B one obtains "nally 1 C " Q : F\ : Q . AB b A  B

(4.25)

(4.26)

(4.27)

We shall discuss the problem of the vanishing eigenvalues of F in Section 4.2.3. The formula (4.27)  encompasses the #uctuations *QK , obtained for c"d. A Remark. The above results have been obtained indirectly by using the identity (2.38). We could also have found them directly by solving the coupled equations of Section 3.3. This is feasible, but not straightforward. Indeed, having de"ned R and H+R, by Eq. (4.17), we can "rst check that the solution of Eqs. (3.52) is given by R ?B (u)"e\SL  @ R eSL  @, R B? (u)"eSL  @ R e\SL  @ .

(4.28)

The proof uses the identity (4.32) below where U is replaced by exp uL  . We can then solve the coupled equations (3.37) and (3.43), together with the required boundary conditions T B (0)"1 and T ? (b)"exp(!L  ), which results in T B (u)"eSL  @ e\SHY+RY, eSL  @ , T ? (u)"e\SL  @ e\@\SHY+RY, e\SL  @ .

(4.29)

We "nally recover the characteristic function (4.16) by means of (3.46)}(3.47) and (3.49)}(3.50). If L  commutes with R, and therefore with H+R,, the expressions (4.28) and (4.29) simplify. Otherwise, when some QK -invariances are broken, the u-dependence in the solution (4.29) of the A evolution equations (3.37) and (3.43) is no longer as simple as it is for m"0 or for the exact solution

R. Balian et al. / Physics Reports 317 (1999) 251}358

287

[given by DK (u)"exp(!uKK ) and by Eq. (2.13) for AK (u)]. In particular ln T B (u) is not proportional here to u, but contains contributions in u, u,2 Although, in the case of commuting observables, we were able to indirectly solve the coupled equations (3.37) and (3.43), it would not have been easy to guess a priori the complicated u-dependence of the solutions (4.28)}(4.29). 4.2.3. Broken invariances Since equilibrium is described by the minimum of the grand potential F +R,, given by Eq. (4.20), % the eigenvalues of the matrix F de"ned in Eq. (4.21) are non-negative. We have, however, to  consider the case of vanishing eigenvalues for which the inversion of F in Eqs. (4.24), (4.26) raises  a problem. The relation (4.49) below shows that such vanishing eigenvalues of F give rise to the  usual spurious modes of the RPA. A characteristic situation where this happens is that of broken invariance. For instance, in the mean-"eld description of pairing correlations, the particle-number NK symmetry is broken. For convenience, in this Subsection, we shall use this operator as the representative of all broken symmetries, although the whole discussion applies to any conserved single-quasi-particle observable. By de"nition of a conserved observable, NK satis"es the commutation relation [KK , NK ]"0. When the NK -invariance is not broken, NK commutes with the approximate density operator DK (b) and we  have both [R , N]"0 and [H , N]"0, where the matrix N has been de"ned in Eq. (4.12).   However, when the NK -invariance is broken, N does not commute with R and H despite the   commutation of NK with KK in the Fock space. In the HFB theory, this non-commutation manifests itself by the occurence of non-zero o!-diagonal elements in Eqs. (3.10) and (3.39). Let us now start our analysis of the meaning of Eq. (4.27) by recalling how the NK -invariance, whether or not it is broken, is re#ected on our variational objects. We can associate with NK a continuous set of unitary transformations ;K ,exp ieNK which leave KK invariant. They are represented in the 2n-dimensional j-space of Section 3.1 by unitary matrices U,e CN .

(4.30)

From the de"nition E+R,"Tr DK KK /Tr DK and the commutation [KK , NK ]"0, we "nd Tr ;K DK ;K \KK E+R," "E+UR U\, , Tr ;K DK ;K \

(4.31)

for any R satisfying the symmetry condition (3.9). We also have S+R,"S+UR U\, since S+R,, de"ned by Eq. (4.6), depends only on the eigenvalues of R. Thus, if F +R, reaches its minimum for % some R satisfying (4.2) and if the NK -invariance is broken, the transformation UR U\ generates   a set of distinct solutions which give the same value F +UR U\,"F +R ,. This implies that the %  %  second derivative F of F +R, has a vanishing eigenvalue. To "nd the associated eigenvector, we  % note that the "rst-order variation of Eq. (4.31) with respect to R, together with the de"nition (3.38) of H+R,, yields UH+R,U\"H+URU\, .

(4.32)

Hence, letting R"R , taking an in"nitesimal transformation U and using dH"E : dR (a   consequence of the de"nition (4.19) of E ), we "nd  [N, H ]"E : [N, R ] . (4.33)   

288

R. Balian et al. / Physics Reports 317 (1999) 251}358

The same argument holds when E+R, is replaced by F +R,, but then the "rst derivatives of % F vanish at R"R and Eq. (4.33) is replaced by %  F : [N, R ]"0,  

[N, R ] : F "0 .  

(4.34)

These equations express that the commutator [N, R ], which does not vanish when the NK  invariance is broken by R , is a right or left eigenvector of the matrix F for the eigenvalue 0.   We are now in position to discuss the expression (4.27) for the correlation C of two conserved AB observables QK and QK when some invariance (here the NK -invariance) is broken. Depending on the A B nature of the observables, we can distinguish three cases. E (1) If both QK and QK are such that A B Q : [N, R ]"1[QK , NK ]2"0 , 

(4.35)

in particular if QK and QK commute with NK (as well as with KK ), their expectation values (4.9) are A B well de"ned, although R is not. Indeed, all the density matrices UR U\ for which F +R, is   % minimum, will provide the same value for 1QK 2. The condition (4.35) ensures that the equation A (4.24) for RdR/Rm " has a "nite solution, that we can denote !b\F\ : Q . Due to Eq. (4.34) B K  B this solution is not unique; it is only de"ned within the addition of ie [N, R ], where e is  arbitrary. The condition (4.35), however, implies that this addition is irrelevant. Hence the correlation (4.27) is xnite and well-dexned, in spite of the vanishing eigenvalue (4.34) of F . This  conclusion holds in particular for the #uctuation *N K  obtained when QK "QK "NK . A B E (2) If QK , but not QK , satisxes Eq. (4.35), the expectation value 1QK 2 depends on the speci"c B A A solution R of Eq. (4.2). This feature is characteristic of a situation with broken NK -invariance.  Moreover, the additive arbitrary contribution ie[N, R ] to RdR/Rm " also contributes to Eq.  B K (4.27), so that the correlation between QK and QK is xnite, but ill-dexned. A B E (3) Finally, if both QK and QK violate Eq. (4.35), Eq. (4.24) has no "nite solution. The #uctuation of A B QK , in particular, is inxnite, consistently with the fact that 1QK 2 is itself ill-de"ned. A A These conclusions hold without any changes for the breaking of invariances other than NK . 4.3. Fluctuations and correlations of non-conserved observables In this section we consider the general case of single-quasi-particle observables QK which do not A commute with the operator KK . 4.3.1. Kubo correlations When the observables QK do not commute with KK , the expression (4.16) is still the variational A approximation for ln Tr exp[!bKK ! m QK ]. By expansion in powers of the m 's, this expression A A A A generates cumulants of the Kubo type. In particular, while its zeroth and "rst-order contributions coincide with Eqs. (4.7) and (4.9), respectively, the second-order contribution, Tr e\@)K @ du eS)K QK e\S)K QK  A B!1QK 21QK 2 , C) "C) , A B AB BA b Tr e\@)K

(4.36)

R. Balian et al. / Physics Reports 317 (1999) 251}358

289

yields the Kubo correlations between the QK 's, which are thus given variationally by Eq. (4.27), A that is 1 C) " Q : F\ : Q . B AB b A 

(4.37)

4.3.2. Standard correlations: derivation In order to obtain our variational approximation for the true correlations C " AB 1QK QK #QK QK 2!1QK 21QK 2 between QK and QK , we have to return to the coupled equations of B A A B A B  A B Section 3.3, facing their u-dependence and expanding them in powers of the sources m . Section 4.1 A gave the zeroth-order solution, which happened to furnish both the thermodynamical potential and the expectation values. Fluctuations and correlations are provided by the next order that we now work out. The "nal result is expressed by Eq. (4.63). We start from the fundamental relation (2.22) that involves the stationary operator DK (b) and the operator AK (m) de"ned in Eq. (1.6); more explicitly it reads: Tr DK (b) dv e\T BKB/K BQK e\\T BKB/K B Ru A  "! . (4.38) Tr DK (b) e\ BKB/K B Rm A As discussed in Section 4.1.2, Eq. (4.38) taken at m"0 yields 1QK 2 as the expectation value (4.9) of A QK in the state DK (b) already determined in Section 4.1.1. The correlations C and the #uctuations A  AB *QK "C are obtained by expanding Eq. (4.38) up to "rst order in the m's which appear both A AA explicitly and through the expansion DK (b),DK (b)#DK (b)#2,DK (b)# m DK (b)#2    B B B of the stationary state DK (b). For shorthand we denote with an index d the coe$cient

(4.39)



RDK (b) (4.40) DK (b), B Rm B K of m in the "rst-order contribution to DK (b). Note that, thanks to the variational nature of the A method, only this "rst-order contribution is needed to evaluate the wanted second derivatives of the characteristic function. By using the formalism of Sections 3.1 and 3.2, we can rewrite Eq. (4.38) as

 

 Ru "!q ! tr R ?B (b) dv e\TL  Q eTL 

A  L A Rm  A  "!q ! tr R B? (b) dv eTL  Q e\TL  . (4.41) A  L A  The correlation C is the derivative of this expression with respect to m , evaluated at m"0. AB B Contributions arise here from the dependences in m of both R ?B (b) (or R B? (b)) and L  , B m Q . As in Eq. (4.39), we expand the matrices R G (m) and T G (m) around m"0, introducing the A A A same notations R G ,RR G /Rm " and T G ,RT G /Rm " as in Eq. (4.40) (the symbol [i] stands A A K A A K

290

R. Balian et al. / Physics Reports 317 (1999) 251}358

for [a], [d], [ad] or [da]). We obtain C "!tr Q R ?B (b)#tr [Q , Q ]R AB  L A B  L B A  "!tr Q R B? (b)!tr [Q , Q ]R  L A B  L B A  "!tr Q +R ?B (b)#R B? (b), . (4.42)  L A B B We now work out the coupled equations of Section 3.3 so as to "nd an explicit expression for Eq. (4.42). First we note, from the boundary conditions (3.44), that we have T ? (b)"1, T ? (b)"!Q ,  B B and hence, using Eqs. (3.8), (3.24) and (4.3):

(4.43)

R ?B (b)"R B (b)!(1!R ) Q R , B B  B  (4.44) R B? (b)"R B (b)!R Q (1!R ) . B B  B  The insertion of Eq. (4.44) into the expression (4.42) of C yields AB (4.45) C "tr [Q R Q (1!R )#Q (1!R )Q R ]!tr Q R B (b) .  A  B   L A B AB  L A  B The "rst term in the right-hand side is what would emerge from a naive calculation using Wick's theorem with respect to the state DK (b). The departure from this result, given by the last term, arises  from the variational nature of the calculation which modi"es DK (b) and hence the associated contraction matrix R B (b) when sources are present. Rather than evaluating R B (b) we shall transform Eq. (4.42) so as to bring the `timea u down B to 0 and thus "nd a more explicit expression for the correlation C . To this aim we expand AB the equations of motion (3.52) for R ?B and R B? around the lowest-order contribution (4.3). This leads to dR B?

1 dR ?B 1 B "! K : R B? , B " K : R ?B , 2  B du 2  B du

(4.46)

where the kernel K is de"ned, for R"R #dR, by   K : dR,d[H+R,, R]  RH+R, "[H+R ,, dR]# dR , R . (4.47)  IIY  RR IIY RR IIY The (2n;2n);(2n;2n) tensor K is the usual RPA tensor associated with the HFB Hamiltonian  H+R, and customarily introduced to describe small deviations around the equilibrium state R .  As in Eq. (4.22), the colon symbol in K : dR stands for half a trace (with a twist in the indices)  involving the last pair of indices of K and the pair of indices of dR; equivalently, K : dR can be   considered as a scalar product in the 2n;2n vector space of dR. It is shown in the Appendix (Section A.3) that K is related simply to the matrix F which, according to Eq. (4.20), characterizes   the second-order variations of the thermodynamic potential F +R, around R . More precisely if %  we de"ne, for any matrix M, the tensor C as the commutator with the HFB density matrix R ,   C : M,[R , M] , (4.48)  







R. Balian et al. / Physics Reports 317 (1999) 251}358

291

the matrix K is (in the Liouville space) the product  K "!C F ,!C : F , (4.49)      again de"ned as in Eq. (4.22). The occurrence of the simple relation (4.49) and its meaning are explained in Sections A.3 and A.5 of Appendix A. The boundary condition T B (0)"1 implies R ?B (0)"R B? (0)"R ? (0),Y . B B B A formal solution of Eqs. (4.46), together with the boundary condition (4.50), gives

(4.50)

R B? (u)"e\SK : Y , R ?B (u)"eSK : Y, B B where Y is still unknown. The correlation (4.42) can now be written as

(4.51)

C "!Q : cosh b K : Y . AB A   With the aim to "nd Y, using Eq. (4.44) and again Eq. (4.51), we can obtain

(4.52)

(4.53) R ?B (b)!R B? (b)"[R , Q ]"2 sinh bK : Y . B B  B   If K had no vanishing eigenvalue, we could eliminate Y between Eqs. (4.52) and (4.53) and write  (4.54) C "!Q : coth b K : [R , Q ] .    B AB  A Unfortunately, due both to the vanishing eigenvalues of C associated with all one-quasi-particle  operators that commute with aRH a and to the possible vanishing eigenvalues of F associated    with broken invariances (see Section 4.2.3), the operator K "!C F has no inverse.    In order to determine the unknown quantity Y,R ? (0) entering Eq. (4.52), we return to the B equation (3.43) which governs the evolution of T ? (u). Eq. (4.43) gives the boundary condition on B T ? (b) while, in agreement with Eq. (3.51), T ? (0) is related to Y through B B T ? (0)"(1!R )\Y(1!R )\ . (4.55) B   We now relate T ? (b) to T ? (0) by expanding Eq. (3.43) to "rst order in the sources around the B B zeroth-order solution (4.1)}(4.4). This yields the di!erential equation 2

dT ?

B !H T ? #T ? H "(E : R ?B )e\@\SH#e\@\SH(E : R B? ) ,  B B   B  B du

(4.56)

where the u-dependence has been omitted in T ? , R ?B and R B? . We used the expansion B B B H+R #dR,KH #E : dR, entailed by the de"nition (4.19) of the matrix E . The quantities     R ?B and R B? are given by Eq. (4.51) in terms of Y and as functions of u. In order to transform the B B left-hand side into a derivative, we multiply Eq. (4.56) by exp (b!u)H on the left and on the   right. As a result, terms of the form eTHMe\TH appear on the r.h.s. of this equation. Using the fact (see Appendix, Section A.3) that the matrix S C in the Liouville space describes the commutation   with bH , we rewrite these terms as  (4.57) eT HMe\T H"eT@ SC : M . Relying now on Eq. (4.51) we note that, after this transformation, the two terms on the r.h.s. of Eq. (4.56) depend on u only through exponentials of S C and K . In addition, these terms can be   

292

R. Balian et al. / Physics Reports 317 (1999) 251}358

explicitly integrated over u by means of the identity d S C  K S C  K  K  K [(e! @\S e! @\S @   : E : e! S "G2 @  !1)C\e!S #C\(e!S !e!@ )] ,    du

(4.58)

which is checked by using Eqs. (4.21) and (4.49). Although C has no inverse, the r.h.s. is  nevertheless well de"ned because C\ is always multiplied either on the left side by S C or on the    right side by C F when the exponentials are expanded. This amounts to extending to operators   the natural convention (e?V!e@V)x\"a!b for x"0. Integrating the transformed di!erential equation for T ? (u), with the boundary condition T ? (b)"!Q given in Eq. (4.43), we "nd B B B  H  H   e@\S T ? (u) e@\S #Q B B  @\S S C  K S C  K    )C\eS !(1!e\ @\S "[(1!e  @ @  )C\e\S    #2C\(sinh bK !sinh uK )] : Y . (4.59)      For u"0, Eqs. (4.55) and (4.59) provide the required equation for Y. We can simplify the result by rewriting Eq. (4.55) as e@H T ? (0) e@H"4 cosh[bH ] Y cosh[bH ] B     (4.60) "!2 (sinh S C )C\ : Y ,     an equality that is derived by using the expressions (A.11) and (A.19) of S and C in the basis where   H is diagonal. Again the r.h.s. of (4.60) is well de"ned, despite the vanishing eigenvalues of C . We   eventually "nd Q "2C\ sinh(bK ) : Y . (4.61) B    Assuming that F\ exists, Eq. (4.61) can be inverted in the Liouville space; using the relation (4.49),  this gives 1 K  Y"! F\ : Q . B 2 sinh b K   

(4.62)

4.3.3. Standard correlations: xnal expression It remains only to insert Eq. (4.62) into Eq. (4.52), which yields (4.63) C "Q : (K coth b K )F\ : Q .    B AB  A  We recall that F is the stability matrix of the HFB theory, de"ned by Eq. (4.20), and that K is the   RPA matrix, de"ned by Eq. (4.47). These two matrices are related to each other, and to the HFB reduced density operator R , through Eqs. (4.48) and (4.49).  The "nal result (4.63) is our variational approximation for the correlations C or for the AB yuctuations *QK . The symmetry C "C is a simple consequence of the relation K "!C F , A AB BA    whereas this symmetry was not obvious from our starting equations (4.42), (4.45) or (4.54). The formula (4.63) gives a meaning to the formal expression (4.54) which, due to the vanishing eigenvalues of the commutator C in K "!C F (Eq. (4.49)), was not de"ned. As for     the vanishing eigenvalues of F , our discussion of Section 4.2.3 still holds for Eq. (4.63) as well as 

R. Balian et al. / Physics Reports 317 (1999) 251}358

293

for Eq. (4.27), and all its conclusions remain valid when QK and QK are arbitrary single-quasiA B particle operators. Finally, when one of the observables, say QK , is a conserved quantity (i.e., [QK , KK ]"0), we "nd A A from Eq. (4.34) that Q : K ,!Q : C F "0 , (4.64) A  A   even if the QK -invariance is broken (if it is not, we have Q : C "0). Then the left-action of the A A  operator K coth b K on Q amounts to a multiplication by the c-number 2/b, so that from Eq.    A (4.63) we recover the result (4.27) of Section 4.2.2 for the correlations when one of the two observables QK or QK is conserved. A B 4.4. Diagrammatic interpretation The HFB e!ective Hamiltonian H and the RPA kernel K which occur in the results above also   appear, as well known, in perturbative expansions. The former is associated with instantaneous self-energy insertions, the latter with the iteration of a pair of quasi-particle or quasi-hole lines. It is thus instructive to interpret in terms of diagrams the expressions obtained here variationally. A partial result was already obtained in [37] in the special case when (i) there is no broken invariance and (ii) the observables QK and QK are conserved. Moreover that derivation disregarded A B the di$culty caused by the vanishing eigenvalues of K . We sketch below a more complete proof.  Let us "rst consider the perturbative expansion of the generalized grand potential ln Z,ln Tr e\@)K Y ,

(4.65)

where 1 KK ,HK !kNK # m QK A A b A includes source terms. By splitting KK  into

(4.66)

KK "KK  #
294

R. Balian et al. / Physics Reports 317 (1999) 251}358

Fig. 2. A diagram contributing to ln Z in the HFB approximation. Each line represents the propagator associated with KK  which includes the sources. Such diagrams sum up into the grand potential associated with an independent  quasi-particle system governed by the self-consistent HFB Hamiltonian.

the matrix element Q . The same partial summation which, for ln Z, eliminated all the diagrams AH HY of the type of Fig. 2, now leaves for 1QK 2 the single diagram where Q is closed onto itself by the A A self-consistent propagator which accounts for all insertions. For equal times, this propagator reduces to R . Thus, the expression (4.9) for 1QK 2 is again represented by the sum of the  A contributions of the diagrams with instantaneous self-energy insertions only. A new feature occurs if we expand ln Z up to second order in the sources m . We then obtain the A Kubo correlation (4.36) as



R ln Z . C) " AB Rm Rm K A B

(4.68)

Two dots Q and Q are now inserted in each diagram, but they may or may not lie on the same A B closed loop. In the "rst case, if we restrict ourselves to the diagrams retained above for ln Z, the summation of the self-energy insertions gives rise only to a non-interacting contribution, evaluated with the self-consistent Hamiltonian H . In the second case, this summation leads to the diagrams  of the type of Fig. 3, which exhibit the irreducible skeleton attached to the dots Q and Q . A B From the topology of the primitive diagrams of Fig. 2, we thus obtain for C) the set of bubble (or AB open ring, or chain) diagrams (Fig. 3). These diagrams start from Q and end up at Q . They involve A B any number n"0, 1, 2,2 of vertices
R. Balian et al. / Physics Reports 317 (1999) 251}358

295

Fig. 3. Bubble diagram contributing to the correlation between the observables QK and QK . A pair of lines represents the A B two-quasi-particle propagator C expressed by Eqs. (4.69), (4.70). A vertex represents the two-body interaction, written as a tensor E . For Kubo correlations C) the `timea u associated with QK is integrated over the interval [0, b]; it is "xed at  AB A 0 for standard correlations C . AB

bubble diagrams. In the latter case, the interaction
296

R. Balian et al. / Physics Reports 317 (1999) 251}358

provides 1 C , (uO0) C(u)"! b\C S #u   

(4.71)

C(0)"!bS\ . 

(4.72)

The expressions (4.70) and (4.72) for the limiting value C(0) are consistent with the value (A.26) of the metric tensor !S\ for R "R .  I H Bubble diagrams iterate a vertex associated with the two-body interaction < by means of the propagator C(u). Regarded as a matrix in the 2n;2n Liouville space, the interaction < can be identi"ed with the matrix E (see Eq. (4.19)) of the second derivatives of E+R,, as seen from the  expression (3.31). Finally, the initial and "nal vertices (Fig. 3) associated with the observables QK and QK (the correlation of which we want to calculate) give rise to factors Q and Q , regarded A B A B here as vectors of the Liouville space. Altogether, the contribution of the n-th order bubble diagram to the Kubo correlation C) , AB evaluated with DK (b) as the unperturbed density matrix, is obtained by iterating n times !E with   the propagator C. However, we have to associate a `timea u, running from 0 to b, with the "nal observable QK of the chain of bubbles. This amounts to keeping only the value u"0 in all the A propagators C. We thus obtain, after summation over n, the contribution of the bubble diagrams to C) : AB  C) 
" Q : b\[!C(0)E ]LC(0) : Q AB A  B L 1 :Q . (4.73) "b\Q : B A C(0)\#E  It then su$ces to take into account Eq. (4.72), together with the de"nition E !b\S of F , to    identify C) 
 with our variational expression (4.37). AB For the standard correlations C , we have now to set to zero the `timesa u associated with AB QK and QK . This is accounted for by a summation over u"2nib\ m where m3Z. We thus have A B  C
" Q : b\ [!C(u)E ]LC(u) : Q AB A  B L S 1 C(u) : Q . "b\Q : B A 1#C(u)E  S

(4.74)

The ordering of QK and QK is relevant only in the evaluation of the zeroth-order term, where the A B summand behaves as !C /u for uP$R as seen from Eq. (4.71); the average over both orderings  as in Eq. (1.8) is achieved by taking a principal part in the summation over u for uP$R. By using the expression (4.71) of C(u) we "nd for uO0: 1 1 1 C(u)"! C "! C , b\C S #u!C E  u#K  1#C(u)E      

(4.75)

R. Balian et al. / Physics Reports 317 (1999) 251}358

297

where we utilized the relations (4.21) and (4.49) which introduces K . The summation over  u"2nib\ m is then performed by means of 1 1 1 1 1  " coth bz! , (4.76) b u#z 2 2 bz S where  is meant as a principal part for uP$R, the value u"0 being excluded. The S expression (4.76) also holds for z"0 where it is continuous. We still have to add the contribution u"0, already written in Eq. (4.73). We thus obtain (4.77) C
"Q : ([! coth bK #(bK )\]C #[bF ]\) : Q .      B AB A  The combination of the "rst two terms is de"ned even if K has vanishing eigenvalues, since these  give no contribution to the bracket. The last term may diverge for vanishing eigenvalues of F , in  agreement with the discussion of Section 4.2.3, but it otherwise cancels the second term. Finally, using again Eq. (4.49) so as to exhibit the "nite factor (coth bK )K , we "nd for the contribution of    the bubble diagrams, C
"Q : ( coth bK )K F\ : Q , AB A      B the same expression as our variational result (4.63).

(4.78)

4.5. Comments The formula (4.63), obtained by inserting the trial choice (3.15)}(3.17) in the variational principle (2.5), covers the general case where the single-quasi-particle observables QK and QK do not commute with the A B operator KK . Within the same variational framework, we have also found the simpler expression (4.37) or (4.27), either for the Kubo correlations or in the case where one of the two observables is conserved. We may regard all these formulae as quantal extensions of the Ornstein-Zernike approximation. Indeed, the latter can be derived [38] from a classical version of the variational principle (2.44) where KK is replaced by the operator KK  which includes the sources (see Eq. (2.36)). The stability matrix F of the HFB theory (Eq. (4.21)) as well as the associated RPA kernel K (Eq.   (4.47)) have emerged naturally in a variational frame (which, however, is not of the Rayleigh-Ritz type) without any other assumption than the independent quasi-particle trial choice (3.15)}(3.17) for the operators DK (u) and AK (u). In the present formulation, the matrices F and K have   automatically come out in a speci"c way, entailed by the variational approach; this should be contrasted with the various fashions in which the so-called RPA approximation can be implemented either in static or in dynamic many-body problems. In particular, both matrices F and  K appear in the general expression (4.63) while only F enters Eqs. (4.27) and (4.37).   The results of Section 4.1 show that, in our variational scheme, the evaluation of the thermodynamic quantities and of the expectation values 1QK 2 requires only the use of the standard HFB A theory, without any RPA-type corrections, in contrast to the evaluation of the correlations and #uctuations. Far from being an inconsistency, this disparity is needed to ful"l the thermodynamic relations, as illustrated by Section 4.2.1. These would be violated, either by a naive use of HFB to evaluate the #uctuations as in Eq. (4.13), or in apparently more sophisticated approaches where thermodynamic quantities (such as the energy) or expectation values 1QK 2 would include RPAA type corrections.

298

R. Balian et al. / Physics Reports 317 (1999) 251}358

We have given in Section 4.4 the diagrammatic interpretation of our variational results. Thermodynamic quantities and expectation values 1QK 2 are represented by the sum of all the A diagrams (of the type of Fig. 2) in which any self-energy insertion is instantaneous; this sum reduces to the term without interactions when the unperturbed Hamiltonian is chosen to be the HFB self-consistent one aRH a. On the other hand, with the same choice of the unperturbed Hamil  tonian, our variational approximation for the correlations is represented by the sum of the bubble diagrams (of the type of Fig. 3). These diagrams, which also appear in the RPA, involve the interaction, not only through H but also directly. Let us stress again that they do not  contribute to the variational expressions for ln Z and 1QK 2. This di!erence in the diagrammatic A expansions is in fact logical, as shown in Section 4.4. Indeed, it is a consistent scheme to retain, for the perturbative expansion of the characteristic function ln Z, the set of diagrams of the type of Fig. 2 and then to expand the outcome in powers of the sources. At the zeroth and "rst order in these sources, this yields the HFB result for ln Z and 1QK 2. At the second order, in this perturbative A treatment as well as in our variational one, the bubble diagrams of Fig. 3 thus emerge from the HFB diagrams of Fig. 2 contributing to ln Z. Note again that the systematic insertion of bubble diagrams (a procedure sometimes called `RPA approximationa or `standard RPAa) in the calculation of all the physical quantities (such as the energy) would not be consistent in the framework of generating functions, and would not preserve thermodynamic consistency, as discussed for instance in Ref. [39]. The diagrammatic interpretation of the variational approximation raises also the question of its status with respect to the Pauli principle. Indeed, the subset of diagrams like the one in Fig. 2 for ln Z (or ln Z) or the one in Fig. 3 for the correlations violate the Pauli principle in the intermediate states, since exchanging the end points of two lines may lead to diagrams which have a di!erent topology and which are not included in the initially retained subseries. Nevertheless, the Pauli principle is properly taken into account at each step of our variational approach. Thus, deciding whether an approximation violates the Pauli principle or does not is an ambiguous question, as the answer depends on the path which leads to this approximation. Another variational principle, based on the maximization of the entropy rather than on the Bloch equation (2.4) for the characterization of the state, has been previously used [37] to evaluate the correlations C . While the formula obtained for the Kubo correlations was the same as Eq. AB (4.37), the result for ordinary correlations di!ered from Eq. (4.63). Its perturbative interpretation also involved the set of bubble diagrams. However, starting from second order, only the u"0 contribution to each diagram was included (see Section 4.1 of [37]). In contrast, our result (4.63) includes the summation over all values u"2nib\m. This appears as another advantage of the present variational principle based on Bloch's equation.

5. Projected extension of the thermal HFB approximation 5.1. Generalities and notations In Sections 3 and 4 we have dealt with equilibriums in a grand-canonical ensemble. However, as discussed in the Introduction, many problems involve data of type I that prohibit #uctuations for the associated conserved operators. For instance, in the case of the canonical ensemble, the state is

R. Balian et al. / Physics Reports 317 (1999) 251}358

299

required to have a well-de"ned particle number N . Trial states of the type (3.4) that were used in  Sections 3 and 4, namely exponentials ¹K of single-quasi-particle operators, violate this requirement. We want, nonetheless, to keep the convenience associated with the operators ¹K and the Fock space. To resolve this dilemma, in conjunction with the variational functional (2.5) we shall use density operators DK involving projections, such as in Eq. (1.11). For instance, a natural choice for a trial state with a well-dexned particle number is (5.1) DK "PK ¹K B PK  , , , where PK  is the projection , 1 L PK " (5.2) dh e\ F, e F,K , , 2n  onto the subspace with the particle number N . As in the preceding sections, the variational  parameters associated with the choice (5.1) are contained in the operator ¹K B which belongs to the class (3.4). We shall consider successively the two types of situations that were analyzed in the Introduction. (i) The "rst type (Sections 5.2 and 5.3) corresponds to the case where the choice (5.1) involves an operator ¹K B which breaks a symmetry of the Hamiltonian. An example is that of a system in canonical equilibrium when pairing sources are added. Then, the operator ¹K B does not commute with PK , although the exact state DK "e\@)K does. In this case, projections are necessary on both , sides of ¹K B . (ii) In the second type (Section 5.4), no symmetry is broken by the operator ¹K B . This situation is illustrated by the thermal Hartree-Fock approximation. Then, the operator ¹K B does not violate the N-invariance and DK can be de"ned with only one projection on either side of ¹K B ; this projection is needed because ¹K B acts in the unrestricted Fock space. Whatever the case, projected trial states are not easy to handle in the form (5.1). In particular, they do not satisfy directly Wick's theorem which has been an important tool in Sections 3 and 4. In order to retain the simple properties of the exponential operators ¹K , we wish to express each projection as a sum of such exponentials. This is feasible when the projection is associated with some underlying group of operators ¹K E having single-quasi-particle operators as generators. In this case the projection can be expressed as a weighted sum of the operators ¹K E over the elements g of this group. For instance, the projection PK  in Eq. (5.2) is built from the group elements , ¹K E ,e F,K (g,h) which are indexed by the parameter 04h(2n. As a consequence, our projected trial states will appear as sums of exponential operators of the ¹K type, each of which is easy to handle. Likewise, if we only wish to eliminate the odd (or even) particle-number components, we should introduce a projection on even (g"1) or odd (g"!1) numbers, which can be expressed as the sum



PK "(1#ge L,K ) , g"$1 E  over the two elements of the group e F,K with h"0 and h"n. The projection on a given value M of the z-component of the angular momentum,



1 L dh e\ F+ e F(K X , PK " + 2n 

(5.3)

(5.4)

300

R. Balian et al. / Physics Reports 317 (1999) 251}358

is similar to Eq. (5.2). To project in addition on a given value of the total angular momentum JK , we should introduce the rotation operators ¹K E "RK (X) in the Fock space; here the group index g denotes the three Euler angles X, and RK (X) has again the form of an exponential of single-particle operators. Integrating over the Euler angles with the rotation matrices D( (X) yields the operator +) 2 J#1 dX D(H (X)RK (X)" "cJM21cJK" , (5.5) PK ( " +) +) 8n A where "cJM2 denotes a basis in the Fock space labelled by J, M and another set of quantum numbers c. When K"M, Eq. (5.5) de"nes the projection on given values of J and M. The projections PK over the spaces with a given spatial parity can be expressed as ! (5.6) PK "(1$e L,K \) , !  where the single-particle operator





1 NK , dr [tR(r)!tR(!r)][t (r)!t (!r)] , (5.7) \ 4 N N N N N counts the number of particles lying in the single-particle states which are odd in the exchange r, pP!r, p (the factor  in Eq. (5.7) accounts both for the normalization of the basic antisymmet ric states [tR(r)!tR(!r)]"02/(2 and for the double-counting in the integration over r). N N In order to encompass formally all these situations, we shall denote by



(5.8) PK " ¹K E , E the decomposition of our relevant projection as a linear combination of the group operators ¹K E . The symbol  2 denotes a weighted sum or integral over the elements of the group. In particular, E gH. The operators for Eq. (5.2), it stands for (2n)\Ldhe\ F,2 while for (5.3), it stands for   H  ¹K E have the form (3.4) of an exponential of a single-quasi-particle operator, which we shall denote as !l E !aRL E a. The coe$cients l E and L E depend on the group index g; for instance, in Eq.  (5.2), we have l F "!ihn/2 and L F "!ihN (04h(2n), and in Eq. (5.3), l H "!ijnn/2 and L H "!ijnN ( j"0, 1), where N is the 2n;2n matrix de"ned by (4.12). As in Section 3.1, we denote the 2n;2n matrix exp(!L E ) as T E and the scalar exp(!l E ) as q E . We shall have to deal with products of two to four ¹K -operators, such as ¹K E ¹K B ¹K EY or ¹K E ¹K B ¹K EY ¹K ? ; from the algebra properties (3.23)}(3.24), these products have still the same ¹K -form and they are characterized by matrices such as T EBEY ,T E T B T EY ,

(5.9)

and scalars q EBEY "q E q B q EY . There is a contraction matrix R EBEY associated with the matrix T EBEY through Eq. (3.8); R EBEY , as T EBEY , depends on the ordering of the original operators. The equations that we shall derive in the next Subsections can also handle cases which are more general than pure projections. For instance, in the description of rotational bands of axial nuclei, it is often a good approximation to assume that the quantum number K associated with the

R. Balian et al. / Physics Reports 317 (1999) 251}358

301

component of the angular momentum on an internal axis of inertia has a well-de"ned value K .  Then rotational states with given J and M values can be sought with the trial form (5.10) DK "PK ( ¹K B PK (  , )+ +) where the operators PK ( de"ned by Eq. (5.5) are not projections when MOK . Thermal +)  properties of triaxial nuclei, which involve sums over the K quantum number, can also be handled by our method. Another example is the case where a symmetry (such as NK ) is broken. Suppose that we still want to work in the grand canonical ensemble but also want to suppress the o!-diagonal contributions between subspaces with di!erent values of N. We are led to consider operators of the form

  

L 1 L PK DK PK " d

d  eG(\(Y,e\G(,K DK eG(Y,K , , , 4n   , , 1 L " d e\G(,K DK eG(,K , 2n  which are also amenable to our treatment.

(5.11)

5.2. The variational ansaK tze and the associated functional Now we must choose the trial classes for the operators DK (u) and AK (u) to be inserted into the functional (2.5) which, at its stationary point, supplies the characteristic function of interest. For the operator DK (u) we already made the choice





¹K E ¹K B (u)¹K EY , ¹K EBEY (u) , EEY EEY where ¹K B (u) is de"ned by Eq. (3.16). For the operator AK (u) we keep the same choice, DK (u)"PK ¹K B (u)PK "

AK (u)"¹K ? (u) ,

(5.12)

(5.13)

as in Eq. (3.17). Again, as in Section 3.2, the quantities l B (u), L B (u), l ? (u), L ? (u) form a set of independent variational parameters. The exact density operator DK satis"es PK DK "DK PK "DK , and when de"ned in the full Fock space the operator KK commutes with PK . Therefore, it makes no di!erence if we project DK (u) or AK (u) in the functional (2.5); the choice DK (u)"¹K B (u), AK (u)"PK ¹K ? (u)PK is equivalent to (5.12)}(5.13). The results of Section 3.2 are readily adapted through a mere replacement of ¹K B by ¹K EBEY and a summation over the group indices g and g. Thus the normalization Z[ , given by



Z[ "TrAK (u)DK (u)"

ZEEY(u) EEY

(5.14)

with (5.15) ln ZEEY(u)"![l ? (u)#l B (u)#l E #l EY ]#tr ln[1#T ? (u)T E T B (u)T EY ] ,  L replaces Z as de"ned by Eq. (3.25). We recall that, according to Eq. (2.18), Z[ does not depend on u when AK (u) and DK (u) are stationary solutions.

302

R. Balian et al. / Physics Reports 317 (1999) 251}358

The action-like functional (3.32), is now replaced by

 
 @

I"Z[  !



du

dl B (u) 1 dT B 1 # tr T B \R BEY?E

# [E+R BEY?E , ZEEY ! L du 2 du 2 EEY

#E+R ?EBEY ,] .

(5.16)

The boundary term Z[  "TrAK DK (b) is obtained from Eqs. (5.14)}(5.15) by letting u"b in l B (u) and T B (u) and replacing l ? (u) and L ? (u) by the given quantities l  and L  (see Eq. (3.14)). We have used the analogue of the property (3.33) to transform T E \R EBEY? T E into R BEY?E . Moreover, from the commutation [KK , ¹K E ]"0 and in agreement with Eq. (4.31), we "nd that for any R the energy (3.31) satis"es E+T E RT E \,"E+R, ,

(5.17)

a property that we have also used to write E+R EBEY? , instead of E+R BEY?E , (similarly, one has E+R ?EBEY ,"E+R EY?EB ,). The expression (5.16) exhibits the fact that we might as well have projected AK (u) instead of DK (u). Indeed, taking into account the property (5.17), this expression does not change under the replacement of gdg by d and of a by gag. Actually, we can readily check that if one projects AK instead of DK , Eq. (3.32) is directly transposed into Eq. (5.16) through the replacement of a by gag. Likewise, by an obvious extension of Eq. (3.35) where d is replaced by gdg, we "nd for the functional I the alternative form:

 


I"Z[  !Z[ (b)#Z[ (0)#

@

du

dl ? (u) 1 dT ?

# tr T ? \R ?EBEY

ZEEY ! du 2 L du EEY

 1 ! [E+R BEY?E ,#E+R ?EBEY ,] . 2

(5.18)

The expressions (5.16) and (5.18) generalize (3.32) and (3.35) through the averaging over g and g with the weight ZEEY/Z[ . They lead us to introduce the matrices



RM B ? ,Z[ \

ZEEYR BEY?E , EEY



RM ? B ,Z[ \

ZEEYR ?EBEY , EEY

(5.19)

and the pseudo-energies



EM B?,Z[ \

ZEEYE+R BEY?E , , EEY



EM ?B,Z[ \

ZEEYE+R ?EBEY , , EEY

(5.20)

(remember that E+R EBEY? ,"E+R BEY?E , and that E+R ?EBEY ,"E+R EY?EB ,). With these notations we can rewrite Eqs. (5.16) and (5.18) in a more condensed fashion, formally closer to the functionals

R. Balian et al. / Physics Reports 317 (1999) 251}358

303

(3.32) and (3.35):




I"Z[  !

@





dl B 1 dT B 1 du Z[ ! # tr T B \RM B ?

# [EM B?#EM ?B] , L du 2 du 2




I"Z[  !Z[ (b)#Z[ (0)#

@



dl ? 1 dT ? 1 du Z[ ! # tr T ? \RM ? B

! [EM B?#EM ?B] , du 2 L du 2  (5.21)

where RM B ? , RM ? B , EM B? and EM ?B replace R B? , R ?B , E+R B? , and E+R ?B , while Z[ replaces Z. An important change between the present functional and the one of Section 3.2 lies in that EM B? di!ers in general from E+RM B ? , (as well as from E+RM B ? ,) when the operator KK contains a two-body interaction. It implies that the variational outcome Z[  yielded by the functional (5.21) will di!er from that which would be obtained by projecting the optimum state found in Section 3; this holds even for the case (AK "1) ) of thermodynamic potentials. 5.3. The coupled equations As in Section 3.3 we "nd the equations governing the evolution of l B (u) and T B (u) by requiring that the functional (5.16) be stationary with respect to l ? (u) and T ? (u). The resulting di!erential equation for l B (u), dT B 1 dl B 1 " tr T B \RM B ?

# [EM B?#EM ?B] , 2 L du 2 du

(5.22)

is similar to Eq. (3.36), with RM B ? , EM B? and EM ?B replacing R B? , E+R B? , and E+R ?B ,. The di!erential equation for T B (u) di!ers from Eq. (3.37) by an explicit summation over the group indexes. It reads:







dT B 1 # [H BEY?E T B #T B H EY?EB ] ZEEYT E (1!R BEY?E ) du 2 EEY ;(1!R EY?EB )T EY "0 ,

(5.23)

with kEYE , H G ,H+R G ,# B 1!R G

(5.24)

where i stands for dgag or gagd and where the c-number kEYE is de"ned by B dT B

1 #E+R BEY?E ,!EM B?#E+R EY?EB ,!EM ?B . kEYE, tr T B \(R BEY?E !RM B ? ) B du 2 L (5.25)

 





To derive the equation (5.23) we made use of Eqs. (3.51), (5.22) and of the fact that H+R, satis"es for any R the relation H+T E RT E \,"T E H+R,T E \ ,

(5.26)

304

R. Balian et al. / Physics Reports 317 (1999) 251}358

a consequence of the de"nition (3.38) and of Eq. (5.17). Note that Eq. (5.22) does not imply the vanishing of the integrand of Eq. (5.16) before summation over g and g. Likewise the di!erential equation for l ? (u) is obtained by varying the functional (5.18) with respect to l B (u): dl ? 1 dT ? 1 " tr T ? \RM ? B

! [EM ?B#EM B?] , du 2 L du 2

(5.27)

while the equation for T ? (u) is provided by varying Eq. (5.18) with respect to T B (u) for 04u(b:





ZEEYT EY (1!R ?EBEY ) EEY



dT ? 1 ! [H ?EBEY T ? #T ? H EBEY? ] (1!R EBEY? )T E "0 , du 2 (5.28)

where H ?EBEY and H EBEY? are de"ned according to Eq. (5.24) with i standing for agdg or gdga and with the c-number kEYE replaced by kEEY given by B ? dT ?

1 #E+R ?EBEY ,!EM ?B#E+R EBEY? ,!EM B? . kEEY, !tr T ? \(R ?EBEY !RM ? B ) L ? du 2 (5.29)









Note that l ? #l B disappears from Eqs. (5.23) and (5.28) if we divide them by Z[ . These equations thus couple only the matrices T B and T ? . Once they are solved, l B and l ? are obtained by integration of Eqs. (5.22) and (5.27). In Appendix B, using the Liouville formulation and the notation (5.19)}(5.20), Eqs. (5.23) and (5.28) are rewritten in the forms (B.10)}(B.11) which look similar to Eqs. (3.37) and (3.43). Finally, in agreement with Eq. (2.19), the stationarity of Eq. (5.18) with respect to l B (b) and T B (b) provides the equations Z[ (b)"Z[  ,

R BEY?E (b)"R BEYE (b) ,

(5.30)

which are satis"ed by the same boundary conditions as in Section 3: l ? (b)"l  ,

T ? (b)"T  .

(5.31)

The boundary conditions at u"0 are l B (0)"0, T B (0)"1 as in Section 3. Again, all our variational quantities (l ? (u), l B (u), T ? (u), T B (u)) depend upon the sources m through the A boundary conditions (5.31). Despite their somewhat complex form, the coupled equations (5.23) and (5.28) entail some simple consequences. First, using Eqs. (5.22) and (5.27), one notes that ZEEY varies with u according to d ln ZEEY"kEYE!kEEY B ? du 1 d tr ln(1!R BEY?E )!c (u) , "! 2 du L

(5.32)

R. Balian et al. / Physics Reports 317 (1999) 251}358

305

where the function c (u) is given by the average



1 d c (u)"! ZEEY tr ln(1!R BEY?E ) . L du 2Z[ EEY Summation over g and g of Eq. (5.32) leads to dZ[ "0 , du

(5.33)

(5.34)

which was expected from Eq. (2.18). Thus the ultimate quantity of interest, the characteristic function u provided by the stationary value of ln I, is equal to ln Z[ (u) for any u. One can also check from Eqs. (5.23), (5.28) and (5.32) that d (EM B?#EM ?B)"0 , du

(5.35)

which is a special case of the general relation (2.16). The relations (5.34) and (5.35) generalize Eq. (3.48). When the operator AK is hermitian, the coupled equations (5.23) and (5.28) admit solutions with hermitian matrices T B (u) and T ? (u), and the resulting c-numbers l B (u) and l ? (u) are real. This is seen by grouping in all the equations the terms corresponding to ¹K E , ¹K EY respectively with the associated terms corresponding to ¹K EY \, ¹K E \, as is shown at the end of Appendix B. 5.4. Unbroken PK -invariance The approach above is general and covers situations in which the trial operators ¹K B (u) and ¹K ? (u) do not commute with the projection PK . It applies for instance to "nite systems of fermions in canonical equilibrium in which pairing correlations are operative. In this case the operator KK denotes the Hamiltonian HK alone without the term !kNK , the projection PK is given by (5.2) and the independent quasi-particle trial state ¹K B (u) as well as the trial operator ¹K ? (u) do not commute with NK , nor with the operations ¹K F "e F,K of the associated group. We specialize below to the simpler case when PK commutes with ¹K B (u). 5.4.1. Notation and general formalism As discussed in Section 5.1 and in the Introduction, there exist situations where projection is required while the PK -invariance is not broken by the variational approximation. In these cases ¹K B (u) and ¹K ? (u) commute with all the group elements ¹K E and therefore with PK . Only one projection is then needed on either side of ¹K B (u) and, instead of Eq. (5.12), the trial state can be written as



DK (u)"PK ¹K B (u)PK "¹K B (u)PK "PK ¹K B (u)" ¹K E ¹K B (u) . E A single summation over the group index g is su$cient in Eq. (5.14):



Z[ "TrAK (u)DK (u)" ZE(u) , E

(5.36)

(5.37)

306

R. Balian et al. / Physics Reports 317 (1999) 251}358

while Eq. (5.15) becomes (5.38) ln ZE(u)"![l ? (u)#l B (u)#l E ]#tr ln [1#T ? (u)T E T B (u)] .  L Likewise, the double averaging in Eq. (5.19) or (5.20) is replaced by a single one over g with the weight ZE/Z[ ; for instance, we have





RM ?B,Z[ \ ZER ?BE . (5.39) RM B?,Z[ \ ZER B?E , E E The expressions of the functional I and of the coupled equations are readily derived from those obtained in Sections 5.2 and 5.3 by noting that the matrix T E commutes with T B (u) and T ? (u) and that the factor T EY can everywhere be replaced by T  "1. Since g disappears and the ordering of g with respect to a and d is irrelevant, there remain only two (instead of four) di!erent contraction matrices: R B?E "R BE? "R EB? and R ?BE "R ?EB "R E?B . The functional (5.16) then simpli"es to

 


I"Z[  !

@



dl B (u) 1 dT B 1 # tr T B \R B?E

# [E+R B?E ,#E+R ?BE ,] . du ZE ! du 2 L du 2 E  (5.40)

As before, the stationary value of ln I (our variational approximation for the characteristic function) is equal to ln Z[ which does not depend on u. Accordingly, the coupled equations of Section 5.3 take a simpler form. Note, however, that because in general kE (Eq. (5.25)) di!ers from kE (Eq. (5.29)), the matrix H EB? is di!erent from B ? H B?E and H ?EB di!erent from H ?EB (Eq. (5.24)). By combining the forms taken by Eqs. (5.23), (5.28) and (5.32) in this commutative case and by using Eq. (3.51), we "nd that the averaged contraction matrices (5.39) satisfy 1 d RM ?B" du 2Z[



ZE[H+R ?BE , , R ?BE ] ,

E 1 d RM B?"! du 2Z[



ZE[H+R B?E , , R B?E ] .

(5.41)

E This was expected from the general equations (2.32) and(2.33), since their validity conditions are satis"ed when the operators ¹K B (u) and ¹K ? (u) commute with PK . While Eqs. (3.52), which are appropriate to the unprojected variational treatment, had the same form as the time-dependent HFB equation (with an imaginary time), Eqs. (5.41) deal with the larger set of matrices R B?E and R ?BE and involve a summation over the group index g. The apparent decoupling of Eqs. (5.41) is again "ctitious because the boundary conditions cannot be expressed explicitly in terms of RM ?B and RM B?. 5.4.2. Partition function, entropy In the rest of this section, still assuming that the ¹K E -invariance is not broken, we focus on the case AK "1) for which our variational principle provides an approximation for the exact

R. Balian et al. / Physics Reports 317 (1999) 251}358

307

partition function (5.42) Z "Tr PK e\@)K . . We have shown in Section 2.4 that by shifting appropriately the operator KK , the variational principle (2.5) used with AK "1) also yields the characteristic function and therefore the cumulants of the one-body operators which commute both with KK and with the projection PK . As seen in Section 2.5, the boundary condition AK (b)"AK "1) allows us to restrict the trial classes (5.12) and (5.13) to AK (u)"e\@\S)K  ; (5.43) DK (u)"PK e\S)K ,   the subscript 0 refers to the zero values of the sources m . The trial operator KK is a (u-independent) A  single-quasi-particle operator commuting with PK , (5.44) KK "h #aRH a ,     which depends on the parameters h , a real number, and H , a 2n;2n hermitian matrix obeying   (3.41) and commuting with T E . The choice (5.43) entails DK (u)AK (u)"PK e\@)K "AK (u)DK (u)"DK (b) ,      and the u-integration becomes trivial in the functional (2.5), which reduces to



I"Tr PK e\@)K (1#bKK !bKK )" ZE [1#bh #btr H RE !bE+RE ,] ,    L     E

(5.45)

(5.46)

with ln ZE "!bh !l E #tr ln [1#e\@HT E ] , (5.47)    L where the trial quantities h and H must be determined variationally. In the generalized   contraction matrices R T E

T E

 " , (5.48) RE , H  e@ #T E 1#R (T E !1)  the matrices H , T E and R all commute; for the unity element of the group   (¹K  "1) ; l  "0, T  "1), the matrix R given by Eq. (5.48) reduces to the contraction matrix  R of the usual "nite-temperature HFB formalism (Eq. (4.1)). As in Eq. (2.43), we recognize in the  simple expression (5.46) the standard variational principle associated with the maximization of the entropy compatible with the constraint on 1KK 2. If one introduces the quantities







1 1 RM , ZE RE , EM , ZE E+RE , , Z[ , ZE ,      Z [ Z [ E E E the functional (5.46)}(5.47) takes the more compact form I"Z[ [1#bh #btr H RM !bEM ] .   L 

(5.49)

(5.50)

308

R. Balian et al. / Physics Reports 317 (1999) 251}358

Within the selected variational space, the best approximation for the exact partition function Z (Eq. (5.42)) is supplied by the maximum of I since the r.h.s. of Eq. (5.46) is always smaller than . Z . The maximization with respect to the number h yields .  (5.51) h "!tr H RM #EM .   L  This equation ensures that, according to the discussion of Section 2.5, the maximum of I reduces to Z[ . Requiring now the functional I to be maximum with respect to the matrix H provides 



ZE RE (H !HE)(1!RE )"0 ,     E

(5.52)

with kE , HE,H+RE ,#  1!RE  where

(5.53)

kE,!tr H (RE !RM )#E+RE ,!EM . (5.54)  L    Eq. (5.52), together with the de"nitions (5.53) and (5.54), determines the trial matrix H to which  RE is related through Eq. (5.48). The j-basis where H and RE are diagonal (with eigenvalues    H and RE , respectively) is thus de"ned by the equations H H



ZE RE H+RE , (1!RE )"0, jOk  H  HI I E while the diagonal eigenvalues of the matrix H are given by 



(5.55)



ZE +RE H+RE , (1!RE )#RE (E+RE ,!EM ), , H " (M\)  H  HH H H  I IH E H where M is the symmetric matrix

(5.56)



M , ZE +d RE (1!RE )#RE (RE !RM ), . (5.57) HI  HI H H  H I I E Therefore the projection generates a more complicated relation (Eqs. (5.52) and (5.48)) between H and the RE 's than the relation between H and R (Eqs. (4.1) and (4.2)) given by the grand     canonical HFB formalism. Eqs. (5.52)}(5.54) can also be obtained directly from the evolution Eqs. (5.23) and (5.28) when one takes into account the simpli"cations induced by the choice (5.43)}(5.44) and by the commutation properties resulting from Eqs. (5.41). Indeed, kE is then the u-independent common value of the quantities kE"kE"kE introduced in Eqs. (5.25) and (5.29), while the matrices H G de"ned in B ? (5.24) also become u-independent and all reduce to H ?EB "H ?EB "H B?E "H EB? "HE. Several consistency properties can be derived from the formulation above. When the functional (5.50) is stationary with respect to the number h , the approximate entropy S[ associated with the  approximate normalization Z[ and projected HFB energy EM by the thermodynamic relation S[ "ln Z[ #bEM ,

(5.58)

R. Balian et al. / Physics Reports 317 (1999) 251}358

309

(a specialisation of Eq. (2.28) to the restricted variational spaces de"ned by Eqs. (5.43) and (5.44)) is also related to the approximate density operator DK (b)"PK exp(!bKK ) by the von Neumann   relations (1.2) and (2.45). Using Eqs. (5.47), (5.49) and (5.51), one obtains the more explicit form S[ "!tr [RM ln R #(1!RM ) ln(1!R )]#ln    L



e\J E [det+1#R (T E !1),] .  E

(5.59)

As expected, S[ reduces to the HFB entropy S+R , (Eq. (4.6)) when the group of operators T E is  restricted to its unity element. Furthermore, in accordance with Eq. (2.26), the stationarity of Eq. (5.46) with respect to the matrix H ensures that the equilibrium energy, !R ln Z[ /Rb, is equal to EM .  According to the discussion of Section 2.3, this guarantees that, in agreement with thermodynamics, the projected entropy S[ also satis"es the relation (2.29) where F ,!¹ ln Z[ . % Remark. Through the ansatz (5.43) we have restricted the trial class by choosing an exponential u-dependence for the operators DK (u) and AK (u) and by relating their exponents. Therefore, it is not obvious that the maximum of I as given by Eq. (5.46) is equal to the stationary value of the functional (5.40) which is de"ned within the more general class (5.12)}(5.13) of operators whose u-dependences are not limited a priori. In the light of Section 2.5, the commutation of PK with ¹K ? and ¹K B is su$cient to ensure this equality. Indeed, if we associate a projection PK to AK (u) as well as to DK (u), the conditions (i) and (ii) of Section 2.5 are both satis"ed. Nevertheless let us give a more explicit proof along the line already followed in Section 4.1 for the HFB case. We have to check that the ansatz (5.43) with h and H given respectively by Eqs. (5.51)   and (5.52) satis"es the set of coupled equations (5.22) and (5.23), or (5.27) and (5.28) which express the unrestricted stationarity. We "rst recall that, for the choice (5.43), the matrix R B?E (u) is u-independent and given by Eq. (5.48). As l B (u)"uh and l ? (u)"(b!u)h , the functional   (5.40) is identi"ed with Eq. (5.46) and both Eqs. (5.22) and (5.27) are identical with Eq. (5.51). Eq. (5.23) then becomes



ZE RE (H !HE!eSHHEe\SH) (1!RE )"0 .       E

(5.60)

Moreover, when H is solution of Eq. (5.52) (which implies Eq. (5.55) in the basis where H is   diagonal) one has the commutation relation



H , 



E



ZE RE HE(1!RE ) "0 ,   

(5.61)

so that Eq. (5.60) is satis"ed for all values of u. This last property can also be demonstrated by means of Eqs. (5.41). Likewise, within the replacement of u by b!u, the stationarity condition (5.28) also has the form (5.60) and this is a consequence of Eq. (5.52). Therefore, in the case AK "1) and unbroken PK -invariance, the result of the optimization of the functional (5.16) is identical to the approximation Z[ of the exact parition function Z (Eq. (5.42)) obtained by the maximization of . the functional (5.46).

310

R. Balian et al. / Physics Reports 317 (1999) 251}358

6. Projection on even or odd particle number In this section we apply the formalism of Section 5 to a speci"c problem with the purpose to illustrate the usefulness of the method. In su$ciently small systems of fermions, such as atomic nuclei or metallic grains, for which the Hamiltonian favors pairing correlations, there exists a qualitative diwerence between the properties of systems with even and with odd numbers of fermions, whereas systems with N, N$2, N$4 2 particles display strong similarities. At zero temperature the HFB approximation, which is known to provide a good description of pairing correlations in nuclei, involves a sum of components with the same particle-number parity. On the other hand, at non-zero temperature, any number, odd and even, of quasi-particle con"gurations are allowed by the standard HFB solution. As a consequence the HFB thermal equilibrium state has both even and odd particle-number components. A realistic description of "nite systems with "xed particle number N would require  a full projection as de"ned by Eqs. (5.1)}(5.2). (See Chapter 11 of Ref. [12] for the litterature, starting with Refs. [40}42], about the particle-number projection at zero temperature.) Nevertheless, the similar properties of systems with the same particle-number parity suggest that a "rst signi"cant step towards this goal consists in projecting the independent quasi-particle trial state onto the space with even or odd particle (or equivalently quasi-particle) number by means of the projection PK de"ned in Eq. (5.3). E Using the formulation of Section 5 we want therefore to approximate variationally the characteristic function u (m) of a system with a given parity of the particle number (or N-parity). In E particular, for AK "1) , we shall thus obtain an approximation for the number-parity projected grand-partition function (6.1) Z , e@I, Tr e\@&K , E , , where the prime indicates that the sum runs only over even (g"1) or odd (g"!1) values of N. We still enforce the particle-number average by means of the chemical potential k entering the operator KK "HK !kNK . 6.1. Characteristic function in the N-parity projected HFB approximation In the specialization of the equations of Section 5 to the present case, the sum  in the general E formula (5.8) reduces to a sum of two terms: PK "(¹K  #g¹K  ),(1) #ge L,K ) .  E  In the 2n;2n representation, we have simply

(6.2)

l  "0, T  "1; l  "!inn, T  "e LN"!1 . (6.3)  The commutation of T  and T  with any matrix T of the algebra introduced in Section 3.1 is obvious. This re#ects the fact that all operators of the single-quasi-particle type (3.12), or of the trial forms ¹K B or ¹K ? given by Eq. (3.4), commute with ¹K  and ¹K  and with the projection PK because E they involve only an even number of creation or annihilation operators. Thus, although a BCS type

R. Balian et al. / Physics Reports 317 (1999) 251}358

311

of state (3.16) breaks the particle-number invariance eGF,K , it does not break the N-parity invariance eGL,K , and we can work with the variational spaces DK (u)"PK ¹K B (u), AK (u)"¹K ? (u) , (6.4) E where only one projection PK is needed. E Therefore, in the present section we shall not need the general formalism of projection developed in Sections 5.2 and 5.3; it shall be su$cient to rely on Section 5.4 which is suited to unbroken invariances. The reader who has chosen to skip Section 5 is requested to accept Eqs. (5.37)}(5.40). With the replacement of  by a sum and the simpli"cations associated with Eq. (6.3), they furnish E the expression of the action functional corresponding to the N-parity projection. Let us "rst transpose the results obtained for characteristic functions in Section 5.4.1. The matrices R ?BE (u) and R B?E (u) which occur in the functional (5.40) can here be written explicitly: R G "R G ,

R G

, R G " 2R G !1

(6.5)

with i"ad or da and with R ?B or R B? de"ned as in Section 3.2. The weight factors ZE, given by Eq. (5.38), take now the compact form (6.6) Z"e\J ? \J B +det (1!R ?B ),\ , L Z"rZ, r"+det [N(1!2R ?B )], , (6.7) L and hence the normalization Z[ (Eq. (5.37)) reads E (6.8) Z[ "Z(1#gr) . E  According to the general formalism of Section 2, the evaluation of the approximate characteristic function u(m)"ln Z[ "ln Tr AK (m)DK (b) requires the solution of the coupled equations which deterE mine DK (u) and AK (u). To this aim we will use Eqs. (5.23)}(5.25), (5.28)}(5.29) and(5.32)}(5.33) with the following simpli"cations: (i) suppression of the integration on the group index g which is everywhere set to 0 (¹K EY "¹K  "1) ), (ii) replacement of  by the sum with j"0 or 1, (iii) E H commutation of T H with T ? and T B whenever it is required. Thanks to these simpli"cations we can express the evolution equations in a more explicit form than in Section 5.3 since the group elements T E are now trivial. For instance, Eq. (5.23) for T B (u) writes: dT B

1 1 1 dT B

!gr "! [H M B?T B #T B H M ?B] , B B 2R B? !1 du 2R ?B !1 2 du

(6.9)

where we have introduced the generalized HFB Hamiltonians H M ?B and H M B? de"ned by B B 1 1 H M ?B,H ?B !gr H ?B

, B 2R ?B !1 2R ?B !1 1 1 H B?

, H M B?,H B? !gr B 2R B? !1 2R B? !1

(6.10)

312

R. Balian et al. / Physics Reports 317 (1999) 251}358

in terms of the matrices given in Eq. (5.24). Taking into account that the quantities kEYE de"ned in B Eq. (5.25) now satisfy the relation k"!grk B B

(6.11)

with

   

 

1 R B? (1!R B? ) dT B

k" 2tr T B \ B L 2(1#gr) 2R B? !1 du #E





R B?

R ?B

!E+R B? ,#E !E+R ?B , , 2R B? !1 2R ?B !1

(6.12)

one can also express the matrices H M ?B and H M B? as B B





1 R G

1 2grk B , H M G"H+R G ,!gr H # B 2R G !1 2R G !1 2R G !1 2R G !1

(6.13)

where i stands as above for ad or da. Likewise, the evolution equation (5.28) for T ? (u) can be put in a form similar to Eq. (6.9). In Eqs. (6.9)}(6.10) and (6.13), the subscript d in the matrices H M G accounts for the fact that the c-numbers k and k and therefore the matrices H M G and B ? B B H M G are in general di!erent. ? For g"g"0, the specialization of Eqs. (5.32) and (5.33) to the N-parity projection yields a compact form for the evolution equation of ln Z,





1 dR B?

d ln Z gr "! tr , du 1#gr L 2R B? !1 du

(6.14)

in terms of R B? (or of R ?B ). This equation could also have been obtained from the general property dZ[ /du"0, using Eq. (6.7) and the derivative of r drawn from Eq. (6.8). To complete the E set of equations which may be useful for the calculation of the characteristic function, we can rewrite Eqs. (5.41) in terms of the matrices R B? and R ?B . For the evolution of R B? we get 1 R B? (R B? !1) d ln Z dR B?

1 dR B?

!gr #2 "![H M B?, R B? ] ,  B 2R B? !1 du 2R B? !1 2R B? !1 du du (6.15) where d ln Z/du is given by Eq. (6.14). A similar equation, with da replaced by ad and a plus sign before the commutator in the right-hand side, holds for R ?B . These equations for R B? and R ?B are coupled through the boundary conditions, as were Eqs. (3.52) which they generalize. Remark. Due to the square roots occuring in Eqs. (6.6)}(6.7), we must carefully "x the determinations which control the signs of Z and Z. In Eq. (6.6) we can rely on the fact that, for small values of the sources, the contraction matrix R ?B is nearly real, lying between 0 and 1, which makes the determination unambiguous. In Eq. (6.7) we have introduced the factor N to account for the term !l  of Eq. (5.38), but both the matrices N and 1!2R ?B have positive and negative eigenvalues;

R. Balian et al. / Physics Reports 317 (1999) 251}358

313

moreover, they do not commute when pairing takes place. In order to overcome this di$culty and to de"ne unambiguously the sign of r, we proceed by continuity, starting from the unperturbed state without pairing. In this case, when N and R ?B are simultaneously diagonalized, the two diagonal blocks of N(1!2R ?B ) are identical, due to Eq. (3.9), and provide twice the same factors in Eq. (6.7). We then choose the determination of r as r"“ (1!2R ?B ) , (6.16) H H where R ?B are the diagonal elements of R ?B lying in the block such that the eigenvalues of N in H Eq. (4.12) equal #1. Hence, r is positive if the number of eigenvalues of the single-particle unperturbed Hamiltonian that are smaller than the chemical potential k is even, while r is negative if this number is odd. In the low temperature limit, for AK "1) and arbitrary k, one obtains for the grand potential F "!b\ ln Z[ the expected result, namely the minimum of (E!k1NK 2) for % E either even N (g"#1) or odd N (g"!1). When the interactions are switched on and pairing takes place, we note that within the determinant of Eq. (6.7) the particle number matrix N can be replaced by the matrix N (associated with the quasi-particle number) which commutes with R ?B . We can thus de"ne r by continuity, starting from Eq. (6.16) and keeping the same form (6.16) where the R ?B are now the eigenvalues of R ?B associated with the subspace where N equals #1. H 6.2. The partition function in the N-parity projected HFB approximation If one is only interested in the evaluation of the N-parity projected grand partition function Z (obtained for AK "1) ), or more generally of the characteristic function for conserved singleE particle observables QK (obtained for AK "exp(! m QK ) with [QK , KK ]"0), it is not necessary to A A A A A solve the coupled equations of Section 6.1, which are suited to the evaluation of the most general characteristic function. Simpli"cations occur as in Section 5.4.2 by relating to each other the trial forms of DK and AK . We therefore introduce a trial operator KK of the form (5.44) which depends on the c-number  h and the 2n;2n matrix H ; for the restricted variational spaces, we take   DK (u)"PK e\S)K ,  E

AK (u)"e\@\S)K  . 

(6.17)

The matrices T B (u) and T ? (u) are therefore equal to exp(!uH ) and exp(!(b!u)H ),   respectively, while the u-independent contraction matrices R ?B and R B? are expressed in terms of   H by  1 "R . R ?B (u)"R B? (u)"    1#e@H

(6.18)

We shall also need the specialization of the matrices RE , de"ned in Eq. (5.48), to the simple  projection group associated with PK : E R 1  " R"R , R" .    2R !1 1!e@H 

(6.19)

314

R. Balian et al. / Physics Reports 317 (1999) 251}358

As in Section 5.4.2, the c-number h is given by Eq. (5.51) where RM and EM , de"ned by Eqs. (5.49),  E E take now the explicit forms 1 RM " (R #gr R) , (6.20) E 1#gr     1 (E+R ,#gr E+R,) . (6.21) EM "    E 1#gr  In these equations the quantity r "Z/Z is given by    b  b r " det N tanh H (6.22) "“ tanh H ,   2 2 H H where N stands for the quasi-particle-number operator and H denote the eigenvalues of the H matrix H . The sign of r is obtained by continuity according to Eq. (6.16).   The matrix H is the solution of the self-consistent Eq. (5.52) which, after a suitable reorganiza tion, can be written










b M . H 1!gr coth H "H   2 

(6.23)

Using Eqs. (6.19) and (5.53) to de"ne H from R ,R and H from R, the matrix H M is given by    b b H M "H!gr coth H H coth H . (6.24)  2  2 









In terms of the matrices R and R, the operator H M can be expressed as   b b b H M "H+R ,!gr coth H H+R, coth H !2 gr kcoth H ,     2  2  2 











(6.25)

in which the quantity k (see Eq. (5.54) and note that k"!gr k) is  1 1 H  #E+R,!E+R , . k" tr (6.26)   1#gr 2 L sinh bH   In this projected variational scheme, Eq. (6.23) replaces the standard HFB self-consistent equation H "H+R ,. The logarithm of the exact N-parity projected grand-partition function   Z (Eq. (6.1)) is approximated by ln Z[ , the stationary value of the functional (5.46), as E E H (6.27) ln Z[ "tr ln[1#e\@ ]#tr [bH RM ]!bEM #ln[(1#gr )] .  L  E E   E  L Equivalently, this quantity can be written as






(6.28) ln Z[ "ln Z!gr b k#ln[(1#gr )] ,  E   where ln Z is the logarithm of the HFB partition function, S+R ,!bE+R ,, which is associated   with the independent-quasi-particle contraction matrix R (see Eq. (4.5)). 

R. Balian et al. / Physics Reports 317 (1999) 251}358

315

The entropy S of the exact projected state is approximated by the specialization to the present E case of the quantity S[ given in Eq. (5.59): E S[ "!tr [RM ln R #(1!RM )ln(1!R )]#ln[(1#gr )] . (6.29) E  L E  E    After some rearrangement, S[ takes the form E bH gr  #ln[(1#gr )] ,  tr (6.30) S[ "S+R ,!   E  2(1#gr ) L sinh bH   where S+R , is the usual HFB entropy given by Eq. (4.6). As stressed in Section 5.4.2, the Eqs.  (6.23)}(6.26) satis"ed by H ensure that the approximate self-consistent projected entropy S[ obeys  E the standard thermodynamic relations (2.28) and (2.29) with F "!b\ ln Z[ . % E Since our trial operators DK (u) and AK (u) satisfy the validity conditions of Eq. (2.23), the density matrix which optimizes the partition function can be used to obtain the variational approximation for the expectation value of any single-quasi-particle operator QK of the type (3.12). For instance, the transposition of Eq. (2.23) gives (6.31) 1QK 2"q#tr QRM ,  L E which generalizes the HFB formula (4.9). Altogether, to obtain the variational approximation Z[ for the N-parity projected grandE partition function Z de"ned by Eq. (6.1), one "rst determines the matrix H (jointly with the E  contraction matrix RM ) from the self-consistent Eq. (6.23) which contains (through H M ) the cE numbers r and k, given by Eqs. (6.22) and (6.26), respectively. One then evaluates Z[ by means of  E Eq. (6.28). Approximations for the thermodynamic quantities such as the entropy (6.29)}(6.30), or for the expectation values (6.31) of single-quasi-particle observables, in particular of NK , involve the averaged contraction matrix RM de"ned in Eq. (6.20); this matrix is positive since, from Eqs. (6.18), E (6.19) and (6.22), its diagonal eigenvalues are given by 1!gr coth @H   H . (6.32) RM "R EH H 1#gr  According to the general discussion of Section 2.4, our variational approximation for the #uctuation of NK is given by the thermodynamic relation 1 RRM 1 R1NK 2 E. " tr N *NK " 2b L Rk b Rk

(6.33)

The derivative RRM /Rk is deduced from RH /Rk, itself evaluated from Eq. (6.23). E  For large systems, r is exponentially small while k is proportional to the volume, as shown by  Eqs. (6.22) and (6.26); these behaviours imply that the projection does not modify the HFB results. However, for su$ciently small systems and at su$ciently low temperature, the factors entering r  (Eq. (6.22)) may be close to $1, except near the Fermi surface; the terms involving r may thus be  signi"cant compared to the dominant terms, although "r " is always less than 1.  The results of this section could have been obtained directly from the formalism developed in Section 6.1 by using the speci"c u-dependence of T B and the commutation of H with H M . The  latter property, already established in the Remark of Section 5.4.2 for a more general case, can also

316

R. Balian et al. / Physics Reports 317 (1999) 251}358

be proven from Eqs. (6.14) and (6.15) once the u-independence (6.18) of the contraction matrix R B? is acknowledged. 6.3. The N-parity projected BCS model: generalities We now use the formalism of the previous section by specializing to the BCS-type pairing. In this case, the n creation operators aR can be arranged (possibly by means of a unitary transformation) in G a set of couples +(aR, aR ), 14p4n/2, such that the Hamiltonian (3.27) takes the form N N (6.34) HK " e (aRa #aR a  )! G aRaR a  a . NO N N O O N N N N N NO N This may encompass models for electrons in a superconductor (p then denotes the momentum and a spin pointing upwards, while p denotes the opposite momentum and spin), or for nucleons in a deformed nuclear shell (the couple momentum-spin is then replaced by the component of the total angular momentum along the principal axis of inertia). For the sake of simplicity we have not introduced a magnetic "eld, which would lift the degeneracy between the states p and p . We have also assumed that the diagonal matrix elements G are zero (a non-zero value of G entails a small NN NN state-dependent shift of the single-particle energy e , which leads to more complicated formulas but N does not modify the conclusions). Moreover, we consider the usual case where the quantities G are real and symmetric in the exchange pq. NO Remark. The reader who wishes to begin with this Section 6.3 without having followed the derivations of Sections 5, 6.1 and 6.2 can do so if he is ready to adopt the de"nitions and accept the expressions and equations that we now recapitulate. Our variational space was introduced in Section 6.2, Eq. (6.17), where the projection PK was given in Eq. (6.2) and where the trial operator E KK (see Eq. (5.44)) involved a c-number h (which will play no ro( le in the following Sections) and   a 2n;2n hermitian matrix H . The contraction matrices R and R were de"ned in Eqs.    (6.18)}(6.19) from the matrix H which is itself determined self-consistently by Eqs. (6.23) and  (6.25). The matrix RM and the scalar quantities r , EM and k were given in Eqs. (6.20), (6.21), (6.22) E  E and (6.26), respectively. The approximate projected partition function Z[ , entropy S[ and expectaE E tion values 1QK 2 were given in Eqs. (6.28), (6.30) and (6.31). 6.3.1. The variational space Within the BCS scheme the operator KK can be put into the familiar diagonal quasi-particle form 






KK " h ! e # e (bRb #bR b  ) ,  N N N N N N  N N by means of the canonical transformation





(6.35)

b "u a #v aR , a "u b !v bR , N N N N N N N N N N (6.36) b  "u a  !v aR , a  "u b  #v bR . N N N N N N N N N N (Note that two di!erent sign conventions are currently used in the literature for this transformation.) The parameters u and v can be chosen real and non-negative, with u#v"1. The N N N N variational quantities are thus e , v and h . N N 

R. Balian et al. / Physics Reports 317 (1999) 251}358

317

All the 2n;2n matrices that enter our formalism can then be decomposed into a number n/2 of 2;2 symmetric matrices labelled by p and corresponding to the operators a and aR , plus another N N set of matrices labelled by p which di!ers from the previous one through a change in sign of v . In N particular the p-components of the matrix H (Eq. (5.44)) are  H "e ; , N N N

(6.37)

where ; is the real symmetric orthogonal matrix N




u!v N ;, N N 2u v N N



2u v N N . v!u N N

(6.38)

Accordingly, the p-components of the contraction matrix R (Eq. (6.18)) read 




aRa 1 R " Tr ¹K B (b) N N N Tr ¹K B (b) aRaR N N



aa N N "(1!; t ) ,  N N a  aR N N

(6.39)

where we have introduced for shorthand the notation t ,tanh b e . N  N

(6.40)

Likewise, the matrix R (Eq. (6.19)) has the p-components  R "(1!; t\) . N  N N

(6.41)

The energy E+R ,, associated with the operator KK "HK !kNK and evaluated for Eq. (6.39), is  now given by E+R ," (e !k)[1!(u!v)t ]! G u v u v t t .  N N N N NO N N O O N O N NO

(6.42)

The corresponding matrix H +R , has the p-components N 




e !k H +R ," N N  D N

D N k!e



,

(6.43)

N

with D, G u v t . N NO O O O O

(6.44)

By substituting t\ to t in Eqs. (6.42) and (6.43) we obtain the energy E+R, and the matrix O O  H+R,. In particular, D is replaced by  N D, G u v t\ . N NO O O O O

(6.45)

318

R. Balian et al. / Physics Reports 317 (1999) 251}358

The c-numbers r (Eq. (6.22)) and k (Eq. (6.26)) are respectively equal to  r "“ t ,  N N 1 k" (t\!t )+e !(e !k)(u!v)!u v (D#D), . N N N N N N N N N N 1#gr  N

(6.46) (6.47)

6.3.2. The variational equations The trial parameters e and v must be determined from the self-consistent matrix equation N N (6.23), (6.25). By means of the canonical transformation (6.36), let us write this equation in the basis where H is diagonal. From the diagonal element of Eq. (6.23) we thus obtain  e (1!gr t\)#2gr kt\ N N  N " (u!v)(e !k)(1!gr t\)#2u v (D!gr Dt\) , (6.48) N N N N N N N  NN and from the o!-diagonal element: 0"2u v (e !k)(1#gr t\)!(u!v)(D#gr Dt\) . (6.49) N N N N N N N  NN Together with the preceding de"nitions of D, D, t , r and k, Eqs. (6.48) and (6.49) are the N N N  self-consistent equations which determine the variational parameters e , u and v . N N N Let us further work out these equations. By introducing the weighted averaged pairing gap D#gr Dt\  NN , D , N (6.50) EN 1#gr t\ N and making use of the relation (u!v)#(2u v )"1, we can simplify Eq. (6.49) into N N N N 2u v u!v 1 N N" N N" , (6.51) D e !k e EN N N where e has the familiar BCS form N e ,((e !k)#D . (6.52) N N EN Except for the de"nition of D , Eq. (6.51) is the same as the usual BCS equation for the coe$cients EN u , v of the canonical transformation (6.36). As regards Eq. (6.48), it appears as an extension of the N N BCS relation e "(u!v)(e !k)#2u v D ; actually, we can rewrite Eq. (6.48) as N N N N N N N 2gr 2gr k   e "(u!v) (e !k)#2u v D ! (D!D) ! , (6.53) N N N N N N EN t!rt\ N N t !gr t\ N N N N and hence, by using Eq. (6.51) to eliminate u and v , as N N 2gr k D (D!D) 2gr  EN N N .  e "e ! ! (6.54) N N t !gr t\ t!rt\ e N N N N N The equations above reduce to the usual BCS ones if g is replaced by zero. Then both D and N D become identical with the BCS gap D , while D becomes irrelevant. However here we have EN !1 N N





R. Balian et al. / Physics Reports 317 (1999) 251}358

319

g"1 for even systems, g"!1 for odd systems, so that several quantities D, D and D now occur. N N EN The quantities e are the energies which appear in the diagonal form of the density operator. Eqs. N (6.52) and (6.54) show that these e have no longer the BCS form. We can eliminate u and v from N N N the de"nitions (6.44), (6.45) and (6.50) of D, D and D . These quantities thus satisfy the equations N N EN D (t\!t ) 1 EO O O , (6.55) D!D" G NO N N 2 e O O 1 D t 1#gr t\ t\ N O . EO O D " G (6.56) EN 2 NO e 1#gr t\ N O O Eq. (6.56) exhibits the projection correction to the so-called BCS gap equation, while Eq. (6.55) obviously has no BCS counterpart. Finally, by using (6.53) and the identity u v (Dt !D t\)"0, we can write the expression (6.47) for k in a more explicit form: N N N NN N N t\!1 t\!t 1#gr t\ D (D!D) 2gr N N N EN N N .  k 1# "! N (6.57) 1!gr t\ 2(1#gr ) 1!gr t\ e 1#gr N  N N  N N Altogether, we have to solve the self-consistent Eqs. (6.54) and(6.56) for the two sets e and D . N EN The quantity e is indirectly involved through t ,tanh(be /2). The two c-numbers r and k are N N N  expressed by Eqs. (6.46) and(6.57), respectively. The di!erence D!D, which enters (6.54) and N N (6.57), can be eliminated by means of Eq. (6.55). The coe$cients u and v are given explicitly by N N Eq. (6.51).





6.3.3. Comments In Eq. (6.56) the last fraction, which takes care of the projection, can also be written as 1#(t\!1) f (gr t\) with f (x)"x/(1#x). The positive term (t\!1) di!ers signi"cantly from O N O zero only when be (1. For an even-number projection (g"1), the function f is positive and O varies between 0 and 1/2 while for an odd-number projection (!1(gr (0) f is negative and  varies from !R to zero. The following conclusions can therefore be drawn: E in the sums of Eqs. (6.56)}(6.57), the modi"cations due to projection a!ect predominantly the terms corresponding to the single-particle states near the Fermi surface. E the even-number projected gap is larger than the BCS one while the odd-number gap is smaller. E larger di!erences with BCS are expected for the projection on an odd number of particles. These features will be con"rmed below, analytically (Section 6.4) and numerically (Section 6.5). In particular the fact that the odd-N projection produces more important e!ects than the even-N projection, especially at low temperatures, can be understood as follows. Consider the ordinary grand canonical equilibrium. The associated grand partition function Z"Z #Z is the sum of > \ the projected functions Z and Z de"ned by Eq. (6.1). In grand canonical equilibrium, the > \ probability of "nding an odd number of particles is Z /(Z #Z ), and we obtained for it \ > \ the approximation Z[ /(Z[ #Z[ ) where Z[ is given by Eq. (6.28). In the BCS grand canonical \ > \ E equilibrium, the ground state is built by putting either zero or two particles in each level (p, p ). Hence, in the zero-temperature limit, Z[ /Z[ tends to 0. More generally, at low temperature, Z[ is \ > \ smaller than Z[ since ln Z then behaves as !b(E!kN). The probability for odd N is then smaller > than that for even N, and it even vanishes as ¹P0.

320

R. Balian et al. / Physics Reports 317 (1999) 251}358

A technical consequence is the following. When building up an approximation for an odd-N system, one must carefully enforce the elimination of the even-N components. This was achieved in Section 6.3.2 since all our calculations for g"!1 are completely decoupled from those for g"#1. In contrast, Ref. [43] is based on the approximate evaluation of Z (by the BCS approximation) and of the staggered partition function denoted as Z , which is Z[ !Z[ in our  > \ language. Two gap functions, D and D , are thus introduced, satisfying equations similar to Eqs.  (6.44) and (6.45), respectively. The self-consistent equations which yield the values of D and D are  then derived by optimizing independently Z"Z[ #Z[ and Z "Z[ !Z[ . Since at low > \  > \ temperature Z and Z are both dominated by Z[ , it is not surprising that D and D are close to  >  each other and become equal at ¹"0. The evaluation of Z[ as a di!erence between the \ independently optimized values of Z and Z is however not reliable. In our approach, we  optimize directly the partition function Z[ for odd-N systems. This leads to the coupled Eqs. (6.50), \ (6.55) and (6.56) for D, D and D , whose results signi"cantly di!er from those of Ref. [43]. They N N \N underline the importance, within a variational approach, to optimize the quantity of interest itself. When the system becomes su$ciently large, both r and r k become small. Hence the above set   of coupled equations reduces to the standard BCS equations. The projection on the particlenumber parity therefore produces e!ects only for su$ciently small systems. 6.3.4. Thermodynamic quantities Once Eqs. (6.54) and (6.56) are solved, the thermodynamics is governed by the Eqs. (6.30), (6.21) and (6.28) which yield S[ , EM and Z[ , respectively. To compute the approximation S[ given by Eq. E E E E (6.30) to the exact N-parity projected entropy S , we take into account the twofold degeneracy of E the eigenvalues (e!@CN#1)\ of the contraction matrix R . We "nd that S[ di!ers from the  E independent quasi-particle entropy S , !1 b b (6.58) S ,S+R ," 2 ln 2 cosh e !be tanh e , N N !1  2 2 N N according to








gr  be (t\!t )#ln[(1#gr )] . S[ "S ! (6.59) E !1 1#gr N N N    N It results immediately from Eqs. (6.31), (6.39) and (6.51) that the expectation value of the particle number is given in terms of the chemical potential by





e !k t #gr t\ N N . (6.60) 1NK 2" 1! N 1#gr e  N N This expression is also the derivative with respect to k of F "!ln Z[ /b"EM !S[ /b, which is our % E E E approximation to the logarithm of the exact N-parity projected grand-potential !ln Z /b. The E expectation value 1HK !kNK 2 is obtained from Eqs. (6.21) and (6.42) as





t #gr t\ N 1HK !kNK 2,EM " (e !k) 1!(u!v) N N N E N 1#gr  N t t #gr t\t\ N O , ! G u v u v N O NO N N O O 1#gr  NO

(6.61)

R. Balian et al. / Physics Reports 317 (1999) 251}358

321

while, using Eqs. (6.50), (6.51), (6.60) and (6.61), the energy 1HK 2 is given by





e (e !k)#D t #gr t\ N  EN N . (6.62) 1HK 2" e ! N N N e 1#gr N  N Here again, the deviation from the BCS energy is exhibited by the two terms that contain r .  In accordance with the general properties of our variational scheme, the thermodynamic relations (2.28) and (2.29) are satis"ed by the approximate entropy S[ . The expectation values of the E single-particle observables are given by Eq. (6.31). The speci"c heat, as expected, is equal both to !bdS[ /db and to the derivative of 1HK 2 with respect to the temperature. In Eq. (6.62), the E temperature appears directly in t "tanh(be /2), and indirectly through e , D , r and also through N N N EN  k since 1NK 2 should be kept "xed in the derivative. Finally, the general arguments at the end of Section 6.2 show that the characteristic function for any single-quasi-particle operator commuting with KK , such as NK , is variationally obtained by adding a source term to the operator KK . We "nd therefore that the #uctuation satis"es as it should 1 R1NK 2 , *N" b Rk

(6.63)

where in Eq. (6.60) the derivative should be taken both explicitly and through e , D and r . N EN  A somewhat similar extension of the BCS theory, also involving a projection on the spaces with a given parity of the particle number, has been worked out in Ref. [22]. The results agree with ours for the partition function Z[ . The comparison between our variational Eqs. (6.48)}(6.49) and those E of Ref. [22] appears more di$cult. In particular it is not clear to us why the entropy found in Ref. [22] does not satisfy, in contrast to ours, the thermodynamic relations (2.27) and (2.28). On the other hand, the general discussion of Section 2 as well as Eq. (6.63) show that #uctuations and correlations calculated with our method must de"nitely be di!erent; for instance we "nd here that the #uctuations of the particle number tend to zero with the temperature, as expected from general grounds. 6.4. The N-parity projected BCS model: limiting cases At zero temperature, the BCS (and HFB) solutions for the ground and excited states are exact eigenstates of the number-parity operator, and their eigenvalues $1 depend on the number of quasi-particles added to the even BCS vacuum. We thus expect the ¹"0 limit of the even-N projected solution (g"#1) of Eq. (6.56) to coincide with the BCS ground state. We shall also see below that this same equation produces the lowest one-quasi-particle BCS state when g"!1. This odd-N BCS solution with lowest energy assigns a speci"c role to the twofold-degenerate individual level whose energy e is closest to the Fermi energy k. In this section we use the index p"0 to label any quantity pertaining to this level. For odd-N systems, we shall later on denote by the index p"1 the level with second smallest quasi-particle energy. The ¹"0 odd-N ground state is obtained by creating one p"0 quasi-particle on the BCS vacuum, that is by e!ecting in all sums the substitution (u , v )P(v ,!u ) for only one of the two quasi-particle states 0 or 0 associated     with e . Under this transformation the product u v changes sign and the contributions of these   

322

R. Balian et al. / Physics Reports 317 (1999) 251}358

two quasi-particle states to the gap cancel exactly. Thus for pO0, the odd-N gap equation is D G \O , (6.64) NO e O O O$ both in the odd-N BCS and projected cases. The only formal di!erence with the usual BCS equation, which at ¹"0 describes the even-N systems, is the suppression of the q"0 term in the sum. Nevertheless, through self-consistency, depending on the level density at the Fermi surface, this di!erence may lead to a signi"cant decrease of the pairing gaps that we shall investigate below both analytically and numerically. Since creating a quasi-particle with index p"0 in the BCS vacuum amounts to creating exactly one particle on the level p"0 or p "0 , the approximation above is usually referred in nuclear physics as the blocking approximation. Note that Eq. (6.64) should be solved with the constraint v"N!1, on the average occupancy of the levels N$ N pO0. To describe odd-N systems at "nite temperature, one could think of the following naive generalization of Eq. (6.64); we shall refer to it as BCS1. One of the two states at the Fermi surface (p"0 or p "0 ) is again assumed to be occupied exactly by one particle while the remaining N!1 particles are a!ected by the pairing correlations. The occupations of the pO0 levels are then calculated according to the standard "nite-temperature BCS procedure. The temperature dependence enters the equation for the pO0 gaps only through the usual thermal factor. These prescriptions lead to the equation: 1 D " \N 2

D !1 Ot . G (6.65) NO e N O O O$ In this BCS1 approximation, the unpaired particle which occupies the level at the Fermi surface forbids the formation of the pair (0, 0 ). In contrast, at "nite ¹, our odd-N projected equation (6.56) does not compel any speci"c level to be blocked by the unpaired particle. As the temperature grows, the statistical #uctuations make it possible that other levels than p"0 be concerned. We shall see in Section 6.5.1 (Fig. 4) that this can lead to sizeable di!erences with Eq. (6.65). However, we show below how, at ¹"0, Eq. (6.64) is recovered from Eq. (6.56) when g"!1. D

1 " !1N 2

6.4.1. The extreme low temperature regime In this section, we consider the situation where the single-particle level spacing at the Fermi surface is at the same time larger than the temperature and smaller than the ¹"0 pairing gap. Nuclei and small metallic grains with pairing provide physical illustrations of such a situation. We defer to the next Section the analysis of the case when the single-particle level density is much larger than b and such that it is legitimate to replace discrete sums on level indices by integrals weighted by the level density. In Section 6.5.4, we shall consider still another situation in which the single-particle level spacing and the BCS pairing gap are almost equal. Let us "rst recall the low-temperature limits of some quantities and equations of interest. For large values of b, the usual BCS gap equation becomes D

D 1 !1 O (1!2e\@CO) . " G NO e !1 N 2 O O

(6.66)

R. Balian et al. / Physics Reports 317 (1999) 251}358

323

In Eq. (6.66) the terms proportional to exp(!be ) lead to a decrease of the gaps versus the O temperature ¹"1/b. In the same limit, the BCS entropy (6.58) is approximated by "2 be e\@CN , !1 N N while the particle-number average (6.60) is given by S



(6.67)



e !k (1!2e\@CN) . (6.68) " 1! N !1 e N N In the N-parity projected extension of BCS, the quantity r (Eq. (6.46)) plays a central ro( le. In the  low-temperature limit considered here, its value tends towards 1 and one has 1NK 2

ln r "!4 e\@CN , (6.69)  N up to terms equal to, or smaller than exp(!3bD ). (However, for some odd-N systems, it can E happen that r vanishes with the temperature because one term of the product exactly vanishes; see  Section 6.5.4.) Even particle number. Using Eq. (6.69), the expansion of Eq. (6.56) for g"#1 yields






1 D (6.70) D " G >O 1!2e\@CO!4 e\@CO>CP . NO e >N 2 O O P P$NO At ¹"0 Eqs. (6.66) and (6.70) are identical. In Eq. (6.70), the terms which control the lowtemperature corrections to the gaps are also negative. However, they are now of the order of, or smaller than, exp(!2bD ) while they were of the order of exp(!bD ) for BCS (near ¹"0, > !1  D KD ). This re#ects the fact that any excitation with "xed N-parity requires the creation of > !1  at least two quasi-particles and therefore an energy larger than 2D . Therefore the decrease of the > pairing gaps with increasing temperature will be slower than for the BCS case. The low temperature behaviour of the entropy (6.59) is given by (6.71) S[ "2 be e\@CN#2 b(e #e ) e\@CN>CO . N O > N NON$O N Like for all other thermodynamic quantities, the ¹-dependence is governed by terms smaller than exp(!2bD ). For instance, for the average particle number, the limit of large b is > e !k 1NK 2" 1! N 1!2e\@CN!4 e\@CN>CO . (6.72) e N N O O$N The analysis of Eqs. (6.54), (6.55) and (6.57) shows that, at ¹"0, the energy e of a state in the N vicinity of the Fermi surface is smaller than its BCS counterpart. For instance, assuming that all the matrix elements G (pOq) have the same value GI , one "nds NO D D D GI D > (1!d ) 3 >N! > #d > . (6.73) e "e ! N N e N N 6 e e e  N  

















324

R. Balian et al. / Physics Reports 317 (1999) 251}358

In fact, in the numerical analysis of Section 6.5, we shall "nd that e and e di!er by only few N N percent. Odd particle number. For N odd (g"!1) and ¹ small, we must treat separately the gaps D associated with the quasi-particles pO0 and the gap D associated with the lowest \N \ quasi-particle energy. Let us "rst consider Eq. (6.56) for pO0. At low temperature both the numerator and the denominator in the right-hand-side fraction vanish as exp(!be ). After dividing out by this  common factor, we obtain the well-conditioned equation:










D D D \O# G \!G \O e\@CO\C . G (6.74) NO e N e NO e O  O O O$ As expected, at ¹"0, Eq. (6.74) is identical with Eq. (6.64). As bPR, Eq. (6.56) automatically singles out the level with the lowest quasi-particle energy and generates the odd particlenumber BCS ground state. The blocking approximation thus emerges naturally in the zerotemperature limit. The low-temperature correction to the gap D is positive and depends on terms \N which are proportional to the exponential of the di!erence e !e . Therefore, starting from the O  value provided by Eq. (6.64), one expects the gap to "rst increase with temperature as exp(!b(e !e )).   The case p"0 is slightly more complicated. One "nds that the temperature corrections involve no longer e but the second smallest quasi-particle energy e , the corresponding gap D and the   \ degeneracy g :  1 D " \N 2










D 1 D 1 D D 1 \# \ e\@CO\C . G \O! G \O!G (6.75) D "! G  O O  \ e 2 e g e e 2  O  O  O O$ The ¹"0 limit of this equation provides an expression of D in terms of the D (qO0). \ \O However, D does not enter the equation (6.64) which, at ¹"0, determines all the other gaps, and \ it is thus relevant only at non-zero temperature. Despite their formally di!erent expressions, we shall "nd in the numerical application performed in Section 6.5 that the values of D and D are \ \N equal within a few percent. The low-temperature entropy (6.59) behaves as (6.76) S[ "ln 2# b (e !e ) e\@CN\C . \ N  N N$ When bPR, the limit of S[ is equal to ln 2. This value re#ects the twofold degeneracy of the \ ground state. For the average number of particles, the low temperature limit of Eq. (6.60) yields





e !k 1! N (1!2e\@CN\C) . (6.77) 1NK 2"1# e N N N$ In contrast to the corresponding BCS formula (6.68), the temperature correction depends on the di!erences e !e . N 

R. Balian et al. / Physics Reports 317 (1999) 251}358

325

Using Eqs. (6.54), (6.55) and (6.57) one "nds that, at ¹"0, e is larger than e . If one assumes N N again that all the matrix elements G (pOq) are equal to the same value GI , one obtains NO D GI D \ . (6.78) e "e # (1!d ) \N#d N e N e N N 2 N  To derive this formula, we have taken g "2 as in the numerical examples presented in Section 6.5.  In these applications, we shall also "nd that the magnitude of e !e is small (a few percent). N N





6.4.2. The condensed matter low temperature regime In the physics of mesoscopic superconductors, one often considers a situation in which the twofold degenerate unperturbed energies e are su$ciently dense near the Fermi surface to allow N the replacement of a sum over p by an integral over e,e !k with a weight w equal to the level N $ density. One also assumes that the pairing matrix elements are constant and equal to GI over an interval K around the Fermi energy k, and that they vanish outside. The width of this interval, taken of the order of the Debye energy, is assumed to be much larger than the pairing gap. It must be "nite, however, to ensure the convergence of the integrals. The ordinary BCS gap equation at zero temperature



K



w de $ , (6.79) \K(e#D then de"nes a p-independent quantity D that we shall take in this section as a reference energy. In this "eld of physics, the `low-temperaturea regime actually refers to an intermediate range of temperatures such that the conditions 2 " GI

1 1 ; ;b;w $ K D

(6.80)

are valid. Under these conditions Eqs. (6.66) and (6.67) give, after elimination of the cut-o! K by means of Eq. (6.79), analytical expressions for the low temperature dependence of the BCS pairing gap: D

!1

 

(b)"D 1!



2n e\@D , bD

(6.81)

and of the BCS entropy: (6.82) "2 w D (2nbD e\@D . !1 $ We shall check below (Section 6.5.2) that, in the limit w D<1, the pairing gap D obtained with $ EN the N-parity projection does not depend much on p and that it is almost identical with D (b). This !1 will be illustrated by the numerical results plotted in Fig. 5. Moreover, we shall see that the di!erence between e and e (e)"((e !k)#D K(e#D is small. Under the conditions N N EN (6.80), the approximation S




  






2nD be 3 D e\@CN\DKe\@CD 1# K d(e) 1# # d(e) , b 8D 8bD 2b

(6.83)

326

R. Balian et al. / Physics Reports 317 (1999) 251}358

holds. Hence the expression (6.69) becomes ln r "!4 w  $








2n D 3 e\@D 1# . b 8bD

(6.84)

Even particle number. The gap equation (6.70) provides an analytical expression for the lowtemperature dependence of the even-N gap D (b): > w D (b)"D 1!4n $ e\@D . (6.85) > b





The extrapolation to zero temperature of D is equal to the BCS gap D. Under the same conditions, > one can rewrite the expression (6.71) as (6.86) S[ "8n (w D) e\@D . > $ The BCS entropy (6.82) is larger than S[ . Indeed, the disorder decreases when odd-particle > components are removed. Odd particle number. Here we cannot use Eqs. (6.74), (6.76) as a starting point since they are only valid when w 4b (while the validity of Eqs. (6.70), (6.71) rests on the condition 1/D;b). We must $ directly implement the conditions (6.80) and Eqs. (6.83)}(6.84) in the general equations (6.56) and (6.59). For the (p-independent) odd-N gap D (b), Eq. (6.56) can be reformulated as \ K w de 1 1  2 $ " ! # . (6.87) D (b) 2bD (b) GI K (e#D (b) \ \ \  \ The limit of this equation, in the regime (6.80), yields







1 1 D (b)"D 1! # . (6.88) \ 2w D 4bw D $ $ At zero temperature the odd-N gap D is therefore smaller than the BCS gap by the quantity \ 1/2w . This property has been shown to be related with the blocking of the p"0 level located at $ the Fermi surface, as described by the odd-N gap Eqs. (6.64) or (6.74). In the temperature range de"ned by the conditions (6.80), the gap grows linearly with temperature. In contrast, in the very low temperature regime investigated in Section 6.4.1 (b5w ), the gap grows exponentially as $ exp(!b(e !e ))Kexp(!b/2wD).   $ In the regime de"ned by the inequalities (6.80) the odd-N projected thermodynamic quantities have a larger temperature variation than the corresponding BCS ones. This results from the fact that the excitation energies are then of the order of the inverse of wD and therefore much smaller $ than w\ and therefore D. $ For the entropy (6.59) we obtain






1 1 2n wD $ . S[ "ln 2# # ln \ b 2 2

(6.89)

It is the expression of the Sackur-Tetrode entropy (see Chap. 7, p. 323 of Ref. [8]) for one classical particle of mass D, with an internal twofold degeneracy, enclosed in a one-dimensional box

R. Balian et al. / Physics Reports 317 (1999) 251}358

327

of length 2n w . This result is natural since the elementary excitations concern only a $ single unpaired particle, with energy e KD#e/2D; if we identify e"e !k with the momentum N N of a particle in a box, the identi"cation of the level densities leads to the size 2n w for this $ box. While elementary excitations of even-N systems require energies of at least 2D, the e!ective spectrum for odd-N systems is nearly continuous, with the form m/2wD where m is an $ integer. This occurence of gapless elementary excitations is consistent with recent experiments on the conductance of odd-N ultrasmall aluminium grains [6]. The analogy with the Sakur-Tetrode entropy shows also that in the variational projected solution, contrary to the BCS1 one, the unpaired particle does not remain on the single-particle level of energy e "k. The projected  solution is a coherent superposition of con"gurations in which this particle explores all singleparticle levels such that b"e !k"41. The fact that we recover a formula of classical thermoN dynamics is another indication that, strictly speaking, the regime associated with the conditions (6.80) is not a low temperature regime. Actually, Eq. (6.80) implies that the temperature is su$ciently large so that wD/b<1, which ensures the positivity of the entropy S[ as given by $ \ Eq. (6.89). Ordering of entropies. The inequalities S[ 'S 'S[ implied by Eqs. (6.82), (6.86) and (6.89) \ !1 > are easily understood. When the system is constrained to have an even particle number, we have seen that its lowest excited states are obtained from the ground state by the creation of two quasi-particles, which requires at least an energy 2D. For an odd system, the excitation of the (gapless) unpaired last particle is much easier, consistently with S[ 'S . \ !1 At non-zero temperature the BCS state is a weighted sum over even and odd components. The concavity of the entropy implies that S is larger than the smallest of the quantities S[ and S[ , !1 > \ that is S[ . The inequality S (S[ is made possible by the very small weight of the odd > !1 \ components entering the BCS state at low temperatures. The smallness of this weight also explains why, as already noted in Section 6.3.3, the odd-N projection a!ects the BCS results more than the even one. Indeed, we have met several qualitative di!erences between the properties of odd-systems and their BCS description, while BCS explains fairly well the even systems. This point will also be exhibited by the numerical applications of Section 6.5, especially on Figs. 4 and 9 which show the gaps at the Fermi surface as a function of temperature. The same considerations on the entropy also hold in the extreme low temperature regime of Section 6.4.1. Indeed, Eqs. (6.67), (6.71) and (6.76) again imply S[ 'S 'S[ . \ !1 > 6.4.3. Higher temperatures High temperatures relevant for superconductivity are of the order of the critical temperature, that is ¹ "0.57D. General conclusions on the behaviour of the equations can only be obtained  when the density of single-particle states at the Fermi surface is signi"cantly larger than one, which entails w D<1. In this limit, the quantity “ t tends towards 0 for any value of p and q. Then, $ P$N O P Eq. (6.56) is equivalent to:





D 1 (6.90) D " G EOt 1#g(1!t) “ t . NO e O O P EN 2 O O P$N O This di!ers from the BCS equation by the factor within the braces, which deviates from 1 by a vanishingly small correction. The projected gap values converge (from above when g"1 or from

328

R. Balian et al. / Physics Reports 317 (1999) 251}358

below when g"!1) towards the BCS values as one reaches the critical temperature. In particular, the odd solution converges towards the standard BCS solution and not towards BCS1. In Section 6.5.4, we present a numerical case which does not pertain to the limit w D<1 and $ yields qualitatively di!erent results. 6.4.4. Large systems In large systems at non-zero temperatures, the product r decreases as the inverse of an  exponential of the volume of the system. Then, an iterative solution of Eqs. (6.54),(6.56), (6.62) and (6.57) is attainable by starting from the BCS solution. Since in several places r appears multiplied  by coth(be /2), the e!ect of the projection is surmised to be especially sizeable near the Fermi N surface. In this large system limit, the approximate projected entropy (6.59), for both values of g, di!ers from the BCS entropy by a constant: S[ "S !ln 2 . (6.91) E !1 This result shows that, for a given non-zero temperature, the even-N and odd-N grand-canonical states become equivalent when the size of the system increases. With respect to the BCS grand canonical state, the N-parity projection reduces the number of relevant microscopic states by a factor 2, thence the term !ln 2 in Eq. (6.91). The same feature occurs for a non-interacting system. Eq. (6.91) also shows that for large systems, pairing interactions a!ect the entropy similarly in BCS and N-parity projected grand-canonical states. The quantity k increases linearly with the volume, so that e !e is expected to be dominated N N by the second term in the r.h.s. of Eq. (6.54). Then, at the lowest non-trivial order, the quasi-particle energy can be expressed as D (t\!t ) D (t\!t ) gr t\ OY . EO O O EOY OY (6.92) e Ke #  N G OOY N N e e 2 OY O OOY The self-consistency and the constraint on the particle number moreover shift k and D away from EN their BCS values. 6.5. The N-parity projected BCS model: a numerical illustration In this section we estimate the importance of the corrections to BCS introduced by the number-parity projection. To this aim, we have performed three schematic calculations that exemplify situations encountered in the physics of (i) superdeformed nuclear states around mass A"190, (ii) superconducting mesoscopic metallic islands and (iii) ultrasmall metallic grains. 6.5.1. Heavy nuclei We "rst consider an application to a "eld of nuclear structure physics which, recently, has been the subject of an intense experimental exploration. Superdeformed nuclear states are generally populated via fusion reactions between two heavy ions. Once the colliding nuclei have fused and a few neutrons have evaporated, the compound nucleus decays with a small probability (a few percent) into superdeformed metastable con"gurations via the emission of statistical c-rays. In

R. Balian et al. / Physics Reports 317 (1999) 251}358

329

heavy nuclei, the angular momentum dependence of the moment of inertia observed in superdeformed bands gives strong evidence for the presence of pairing correlations. Most of the presently available data concern the properties of superdeformed ground states or of the lowest excitations. However, it is expected that the next generation of c-arrays, in Europe and in the US, will soon allow investigations of the statistical behaviour of the spectrum and of the properties of heated nuclei. In view of the general scope of the present work, rather than performing a realistic nuclear calculation (see for instance Ref. [44]) that would require the solution of the full extended HFB equations written in Section 6.2, we prefer to investigate a simple model in which the nuclear context outlined above is used as a guidance for selecting the relevant parameter values. Thus we consider twofold degenerate single-particle levels of energies e , regularly spaced with a density N w "2 MeV\. Our active pairing space includes 21 levels and extends 5 MeV above and 5 MeV $ below the Fermi surface, so that the cut-o! K equals 10 MeV. The matrix elements of the interaction are taken to have the form G "GI (1!d ) with GI "0.23 MeV. For this value of GI , NO NO the BCS pairing gap D at zero temperature is 1 MeV, consistent with the data for the even heavy nuclei. In fact, it is partly because D is close to 1 MeV in medium and heavy nuclei that the MeV turns out to be a natural energy unit in nuclear dynamics despite the fact that binding energies are several orders of magnitude larger. For these values of the parameters w , K and GI (or equivalently D), we have solved the N-parity $ projected coupled equations (6.54) } (6.57) with g"$1 as function of the temperature. This provides us with the projected pairing gap D de"ned by Eq. (6.50). Although our interaction is EN very simple, it leads to state-dependent gaps for the BCS as well as for the projected cases. The numerics, however, shows that this state dependence is weak: we "nd that at all temperatures the variation of the gap D versus p never exceeds 9%. We thus content ourselves with plots of the EN ¹-dependent value of the gap associated with the quasi-particle having the lowest energy in either the even or odd system. Hereafter, we denote this quantity D where D "D for N-even projection E E > and D "D for N-odd projection; we denote moreover by D (¹) the ¹-dependent BCS gap E \ !1 associated with the lowest-energy quasi-particle (D (¹),D ). !1 E As a reference energy in all the "gures of the present section, both for the temperature in abscissa and for the ¹-dependent projected gaps D in ordinates, we take the zero-temperature BCS gap E associated with the lowest energy quasi-particle which we denote D (D,D (¹"0)). !1 We "rst consider (upper part of Fig. 4) the even-number case with 1NK 2"20. (Here, 1NK 2 should not to be confused with the numbers of neutrons and protons which, for a heavy nucleus, are of the order of 120 and 80, respectively. 1NK 2 corresponds to the particles within the active pairing space; it does not take into account the nucleons of the fully occupied levels whose energies are lower than k!K/2.) We have plotted the temperature dependences of the BCS gap D and of the projected !1 gap D . As expected from the discussion of Eq. (6.56) in Section 6.3.3, the even-N projected gap is > always larger than the BCS one. The two curves meet at ¹"0 and also at the critical temperature ¹ "0.57D, consistently with the analysis of Section 6.4.3.  The lower part of Fig. 4 presents results for the odd-number case with 1NK 2"21. It compares the temperature-dependent gaps obtained with our odd-N projected method, with BCS and with BCS1. As explained in Section 6.4, in the BCS1 option the level p"0 (with e "k"0) is always  occupied by exactly one nucleon and therefore not available for pairing. We have therefore chosen to show in all three cases the gap D (g"0, $1) associated with the lowest-energy quasi-particle E

330

R. Balian et al. / Physics Reports 317 (1999) 251}358

Fig. 4. Comparison of BCS and number-parity projected gaps at the Fermi surface versus temperature for w D"2, $ w K"20 and 1NK 2"20 or 21. These parameter values are illustrative of the physics of heavy nuclei. The upper $ (respectively lower) part of the "gure corresponds to an even (respectively odd) average number of particles. With respect to BCS, the pairing gap is enhanced by even-N projection and (more strongly) reduced by odd-N projection, except near ¹ . For the BCS1 curve, see Sections 6.4 and 6.5.1. The odd-N projected gap increases with ¹ up to 0.4D. 

whose index is distinct from 0. The odd-N projected curve interpolates the BCS and BCS1 ones. In agreement with the discussion in the introduction of Sections 6.4 and 6.4.1, this curve meets the BCS1 curve at ¹"0 where they both yield the correct result. As discussed in Section 6.4.3, the projected curve terminates at the same critical temperature as BCS. The failure of BCS at low temperature is associated with the small weight of the odd-N components in the BCS state (see Sections 6.3.3 and 6.4.2). The BCS1 approximation, by enforcing the odd parity at zero temperature, approaches at low temperature the projected state better than does BCS. On the other hand, when the temperature increases, the even and odd components tend to have equal weights while the di!erences between even and odd systems fade out. In that limit the BCS1 approximation becomes inadequate since it prevents the formation of a pair in the level e . This inhibits a correct  description of pairing correlations which, at high temperature, take place coherently on all levels in the energy range of the gap for the odd as well as for the even systems.

R. Balian et al. / Physics Reports 317 (1999) 251}358

331

When the temperature grows, the odd-N projected gap begins to increase as expected from (6.74)}(6.88), then reaches a maximum, decreases and converges asymptotically towards the BCS curve near the critical temperature. At low temperatures, there is a signi"cant di!erence between the usual BCS gap D and D . In nuclear physics, this di!erence should be taken into account when !1 \ mass di!erences between odd and even neighbour nuclei are used to extract pairing gap values. In the numerical example of Fig. 4, the condition w D<1 is not satis"ed. It is therefore not $ surprising that the results do not comply with Eqs. (6.85) and (6.88) derived under the conditions (6.80). When gaps and temperatures are measured in D units, the BCS curve is universal up to small di!erences resulting from the details of the spectrum. For instance, the critical temperature always comes out equal (with an accuracy better than 1%) to the value 0.57D which is derived for a continuous spectrum with constant level density. In contrast, as anticipated in Section 6.4, the N-parity projected curves show more diversity. In our model, their characteristics are mostly determined by the product w D and to a lesser degree by the product w K and the ratio K/D. The $ $ quantity w D, `the number of Cooper pairsa, is the number of twofold degenerate single-particle $ levels in an energy interval equal to the gap, while w K is the number of levels in the active pairing $ space. In the nuclear physics case of Fig. 4 we have taken w D"2 and K/D"10, the number w K $ $ of levels being 21. The data corresponding to condensed matter physics are expected to be quite di!erent. 6.5.2. Mesoscopic metallic islands For superconducting materials, K is of the order of the Debye energy while D is of the order of the critical temperature, so that K/D is larger than above. On the other hand, depending on the size of the metallic islands or grains, such as those considered in recent experiments [3}5,45], the product w D can vary between 0.5 and several thousands. $ Let us "rst consider the case w D<1, which corresponds to most of the experimental studies $ performed on superconducting metallic islands. In this limit the e!ects of the N-parity projection on the ¹-dependent pairing gaps D are weak. In order to be able to visualize them, we have ! considered in Fig. 5 not too large a value of w D, taking w D"10 with K/D"20 and hence $ $ w GI "0.34. Thus, the active space comprises 201 twofold-degenerate levels. The average particle $ numbers of the even and odd systems are 200 and 201, respectively. The BCS curve and the even-N projected curve D are almost indistinguishable. In agreement with the discussion of the projected > gap equation in Section 6.3.3, the odd-N projected gap di!ers from the BCS value more than the even-N. At low temperature, it increases and is represented with good accuracy by the linear ¹-dependence of Eq. (6.88). In Fig. 6 we have plotted the entropy of the odd system calculated in four di!erent ways: with or without (i.e. GI "0) pairing interaction and with or without projection on number parity. The odd-N projected entropy S[ tends to ln 2 when ¹P0; above the critical temperature it joins the \ non interacting odd-N projected entropy. In contrast, the BCS entropy vanishes at small temperature; above the critical temperature, it joins the grand-canonical entropy of the non-interacting system. For ¹ larger than 0.35D the odd-N projection lowers the BCS entropy, according to S !S[ " ln 2. In contrast, at small ¹, one has S[
332

R. Balian et al. / Physics Reports 317 (1999) 251}358

Fig. 5. Comparison of BCS and N-parity projected gaps at the Fermi surface versus temperature for w D"10, $ w K"200 and 1NK 2"200 (upper part) or 201 (lower part). Such parameter values are illustrative of superconducting $ mesoscopic islands. The e!ects are qualitatively the same as in Fig. 4 but weaker. Even}odd e!ects disappear above ¹ "0.3D. "

almost undistinguishable from the expression 2n S" w ¹ , 3 $

(6.93)

obtained by assuming a continuous spectrum with constant level density. The only e!ect of the odd-N projection is again a lowering of the entropy by ln 2. In Fig. 7 we show the specixc heat C "d1HK 2/d¹"¹dS/d¹ for the odd system. With or without 4 projection, above the critical temperature ¹ "0.57D, the speci"c heat of the interacting system  follows the linear form (6.93). Both the critical temperature and the magnitude of the discontinuity of C are the same for the BCS and odd-N projected curves. The low temperature behaviour, 4 however, is a!ected by the projection. The enlarged scale inset in Fig. 7 shows that the BCS speci"c heat drops to zero at low temperature, while the odd-N projected curve approaches the value   expected for a single classical particle moving in a one-dimensional box (see Section 6.4.2). As mentioned in Section 6.3.3, the odd-N speci"c heat C takes this value  only over the intermediate 4 

R. Balian et al. / Physics Reports 317 (1999) 251}358

333

Fig. 6. Entropy of an odd-number system for w D"10, w K"200 and 1NK 2"201. The solid curve corresponds to the $ $ variational solution with projection on an odd number of particles. Around ¹"0.1D it has a Sakur}Tetrode form and it tends to ln 2 for ¹P0. The dashed curve (BCS) is obtained with the standard "nite-temperature BCS formalism. The dashed-dotted and dotted curves are obtained by dropping the interactions. The latter unprojected non-interacting case coincides with the entropy of the Fermi gas.

range of temperatures where the quasi-particle spectrum e is well approximated by a quadratic N function of e !k. In Section 6.5.4 we shall investigate in more details its zero-temperature limit. In N the absence of pairing interaction C does not depend on whether odd-N projection is or is not 4 e!ected, and it is again given by the linear function (6.93). 6.5.3. The free-energy diwerence We now consider the di!erence dF between the free energies of the odd and even systems, a quantity which is experimentally accessible [45]. In our formalism this di!erence can be reached by means of the approximate free energies obtained by the Legendre transform, F (1NK 2)"F (k)#k1NK 2 , (6.94) E %E of the N-parity projected BCS grand potentials F "!b\ln Z[ given by Eqs. (2.27) and (6.28). %E E The free energies F do not pertain to canonical ensembles with non-#uctuating particle number, E but to N-parity projected grand canonical ensembles. Nevertheless, we expect that dF does not di!er signi"cantly from F (1NK 2"N#1)!F (1NK 2"N). \ > At low temperatures, the experimental projected free-energy di!erence dF decreases linearly (see Fig. 5 of Ref. [45]) and tends to zero for a temperature lower than the critical temperature. Let us see how this behaviour can be understood within our formalism. First, using the Legendre transform (6.94), we numerically calculate dF from the odd-even grand potential di!erence for given values of ¹ and k, 1 dF ,F (k)!F (k)"! (ln Z[ !ln Z[ ) , \ > % %\ %> b

(6.95)

334

R. Balian et al. / Physics Reports 317 (1999) 251}358

Fig. 7. Speci"c heat of an odd-number system for w D"10, w K"200 and 1NK 2"201. Symbols and abbreviations are $ $ the same as in Fig. 6. The solid curve reaches the classical value for one degree of freedom near ¹"0.1D while the BCS value tends to zero, as exhibited by the inset. The projected and unprojected curves with no interaction are identical and have the linear behaviour of the Fermi gas.

with Z[ given by Eq. (6.28). Fig. 8 shows the variation with temperature of dF, for the same E parameters w D"10 and K/D"20 as in Section 6.5.2 (we have taken as our unit scale D the value $ for the odd system 1NK 2"201, which di!ers from the even system value by less than 0.1%). In agreement with experiment, one observes that the projected and BCS curves are identical for temperatures above 0.4D while, at lower temperatures, dF decreases almost linearly from a zerotemperature value which, by extrapolation, can be estimated to lie close to D. In the BCS case, for a "nite system in grand-canonical equilibrium, the free energy F (1NK 2) is !1 a function of the expectation value 1NK 2 which may vary continuously from 1NK 2"N"200 to 1NK 2"N#1"201. Moreover, RF /R1NK 2"k(1NK 2) is the BCS chemical potential as func!1 tion of 1NK 2 at the considered temperature. Accordingly, we can obtain the di!erence dF "F (N#1)!F (N) by integration of k(1NK 2) between N and N#1. In the !1 !1 !1 present model with N#1 equidistant twofold degenerate single-particle levels, for symmetry reasons, k(N#1) is equal to the energy of the level at mid-spectrum (101th level) which we take equal to zero. The Fermi level k(1NK 2) lies halfway between two consecutive levels for even 1NK 2, and coincides with a level for odd 1NK 2 provided the occupation number is close to 1 (resp. 0) at the bottom (resp. top) of the active space. Finally, when w ¹<1, or w D<1, k(1NK 2) is expected to $ $ vary smoothly, without oscillations. Hence we "nd 1 (1NK 2!N!1) , (6.96) k(1NK 2)K 2w $ from which we derive the BCS di!erence of free energies: ,> 1 k(1NK 2) d1NK 2K! . (6.97) dF " !1 4w , $



R. Balian et al. / Physics Reports 317 (1999) 251}358

335

Fig. 8. Di!erence between the free energies of neighbouring odd (1NK 2"201) and even (1NK 2"200) systems as function of temperature for the same parameter values w D"10 and w K"200 as for Figs. 5}7. The solid and dashed curves $ $ correspond to self-consistent calculations with and without projection on N-parity, respectively. Above the crossover temperature ¹ "0.3D, the even}odd e!ect disappears, as for D in Fig. 5. Below ¹ the projected di!erence dF decreases " E " nearly linearly. It lies close to the dashed-dotted curve which displays the low-temperature perturbative approximation (6.99).

At this level of approximation one "nds that dF is temperature independent, a result supported !1 by the numerical evaluation, as illustrated by Fig. 8. It is also possible to obtain an analytic expression of dF for temperatures in the range (6.80). We "rst calculate dF perturbatively around the BCS solution by means of Eqs. (6.28) and (6.84). If we % neglect the contribution r (k #k ), we obtain  \ > 1 1#r ¹ ¹ dF K1! ln 8n(w D) %K ln . (6.98) $ bD 1!r 2D D D  In order to relate dF to dF by means of Eq. (6.94), we have in principle to control the variation of % k associated with the addition of one particle and the change of N-parity. However, we note that the derivative with respect to k of the approximate expression (6.98) is practically zero. (Indeed, RD/Rk" "0 for symmetry reasons.) This means that the relation between 1NK 2 and k is I the same for the even-N and odd-N projected solutions. Hence this relation is also the same as for BCS (Eq. (6.96)). Since k vanishes for 1NK 2"N#1, the di!erence dF coincides with % F (N#1)!F (N#1), as shown by Eq. (6.94). The extra contribution F (N#1)!F (N) to dF \ > > > is then evaluated by integration of RF /R1NK 2"k(1NK 2) as in Eq. (6.97). We "nd: > 1 ¹ 1 KD! ! ln(8nD w ¹) . (6.99) dF"dF ! $ % 4w 4w 2 $ $ This expression is represented on Fig. 8 by the curve labelled `proj. pert.a. It is very close to the full solution of the projected equations up to ¹"0.3D.






336

R. Balian et al. / Physics Reports 317 (1999) 251}358

Thus, our results agree with the observation [45] that, when the spacing between single-particle levels is small compared to the gap (w D<1), N-parity ewects disappear above some temperature $ ¹ which lies below the critical temperature. Actually, Fig. 8 shows that for w D"10 a sharp " $ crossover takes place around ¹ "0.3D (that is, around 0.5 ¹ ) between the N-parity dependent "  and independent regimes. Below this value, projection strongly modi"es dF, while parity e!ects disappear above. Figs. 5}7 exhibit the identity of the various BCS, even-N and odd-N projected results above ¹"0.4D. More generally, for arbitrary values of w D much larger than 1, the di!erence dF is given at $ low temperature by the approximation (6.99), at higher temperature by the value !1/(4w ) $ associated with the same trivial parity e!ect than for a non-interacting system or for the BCS approximation. The crossover between the two regimes takes place around the solution of the equation 1 ¹ "K . D 8n ¹ " ln w D $ D



(6.100)

This crossover temperature is also visible on Fig. 5, which shows that D (¹) is the same for BCS and E for the two projected cases in the range ¹ (¹(¹ , where ¹ K0.3D for w D"10. "  " $ When the size of the system decreases, w D also decreases. The crossover temperature ¹ thus $ " rises; from Eq. (6.100) one "nds that it reaches the critical temperature ¹ "0.57D for w DK1.5.  $ Thus, if the distances between levels become comparable to the gap, one can expect that odd-even e!ects will take place over the entire range of temperatures where pairing occurs. Physical systems which satisfy this condition seem presently available since one can now prepare and study ultrasmall metallic grains [4,5]. 6.5.4. Ultrasmall metallic grains As an example of modelisation of ultrasmall metallic grains, we have chosen the parameter values w DK1.1 and w K"100. These are obtained for w GI "0.23, and for an active space $ $ $ having 101 twofold degenerate levels. For the odd (resp. even) system we take 1NK 2"101 (resp. 100). In Fig. 9 we present the gap D as function of the temperature for odd and even systems. The E critical temperature associated with the even-N projection comes out somewhat larger than the BCS value. Even more striking is the behaviour of the odd-N projected curve. Then the gap is zero at ¹"0 but takes a non-zero value above some xrst transition temperature up to a second one, lower than the BCS critical temperature. This peculiar behaviour can be explained as follows. At ¹"0, the unpaired particle occupies one of the two single-particle states e "k with a probability exactly  equal to one. This blocks the formation of pairs utilizing the two degenerate states with energy e and creates an e!ective single-particle spectrum with a distance between the levels at the Fermi  surface which is twice as large as D. Consequently, pairing is inhibited. When the temperature grows, the bachelor particle statistically populates other single-particle states, as indicated in the discussion above on the entropy (6.89). This allows pairing correlations on the strategic level at the Fermi surface (e "k), provided the temperature is not too low. One thus observes a reentrance  ewect, with two transition temperatures: one at which the pairing switches on and, at a higher temperature, the usual one at which it switches o!.

R. Balian et al. / Physics Reports 317 (1999) 251}358

337

Fig. 9. Comparison of BCS and N-parity projected gaps at the Fermi surface versus temperature for w D+1.1, $ w K"100, 1NK 2"100 (upper part) or 101 (lower part of the "gure). Such parameters are illustrative of extremely small $ superconducting grains. The critical temperature is raised by even-N projection, lowered by odd-N projection. In the latter case, pairing disappears below a new, lower critical temperature 0.15D, thus displaying a reentrance phenomenon.

The entropy of the even system 1NK 2"100 is shown in Fig. 10. Projection strongly reduces the entropy below the critical point. This projection e!ect is seen to be more important for the interacting system. At small ¹, the four entropies vanish exponentially. As expected the decrease rate of the projected entropy is the fastest. For ¹'0.3D, the grand-canonical entropy of the non-interacting system agrees with the linear form of Eq. (6.93). It bends down below ¹"0.3D under the e!ect of the spectrum discreteness; in the domain ¹'0.3D, projection again induces only a lowering by ln 2. Fig. 11 displays the entropy of the odd system 1NK 2"101. When the pairing interactions are present, the BCS formalism, which only constrains the average number of particles, yields the result 0 instead of the right zero-temperature limit ln 2. This correct limit is obtained by the odd-N projection. When there is no interaction, for ¹'0.3D, the grand-canonical value (dotted curve) is again described by Eq. (6.93). However the entropy tends to ln 4 at small ¹. This comes from the fact that, at ¹"0, the Fermi energy exactly coincides with the energy of the last occupied level

338

R. Balian et al. / Physics Reports 317 (1999) 251}358

Fig. 10. Entropy of an even-number system for w D+1.1, w K"100, 1NK 2"100. The symbols and notations are the $ $ same as in Fig. 6. Projection strongly reduces the entropy.

Fig. 11. Entropy of an odd-number system for 1NK 2"101 and the same parameters w D+1.1, w K"100 as in Figs. $ $ 9 and 10. The symbols and notations are the same as in Fig. 6. The projection lowers the entropy above ¹ , but raises it at  lower temperatures with the limit ln 2 at ¹P0.

whose average occupation is equal to 1/2. When this value is inserted in the standard formula for the grand-canonical entropy of a non-interacting fermion system (while taking into account the twofold level degeneracy), one obtains ln 4 instead of the expected value ln 2. Projection on odd-N parity decreases the entropy by ln 2 and restores the correct limit. Note that this is not a trivial

R. Balian et al. / Physics Reports 317 (1999) 251}358

339

result since the exact particle number is not restored by the projection; only its parity is. In Fig. 11 the reentrance phenomenon manifests itself by the fact that the entropy of the interacting system is larger than that of the non-interacting system over the temperature interval 0.14(¹/D(0.34; otherwise the e!ect of the interactions is rather small. The original features associated with the N-parity projected mean-"eld approximation when w D is close to one are even more visible on the diwerence of the specixc heats of neighbouring odd $ and even systems. In Fig. 12 one sees that the projected curve shows discontinuities for the even-system critical temperature, and also for the two odd-system critical temperatures. The comparison with the same curve with w D"10 shown in the inset of Fig. 12, which corresponds to $ the second derivative of the curves on Fig. 8, displays the striking qualitative change arising from the variation of w D over one order of magnitude. $ 6.6. Discussion We shall resume in the general conclusion the various results obtained above. We focus here on the change in the pairing properties when the size of the systems decreases. The dimensionless parameter characterizing this size appears to be w D, which takes values 10 in Fig. 5 (metallic $ islands), 2 in Fig. 4 (heavy nuclei) and 1.1 in Fig. 9 (aluminum grains). From this standpoint the pairing properties of heavy nuclei (Section 6.5.1) are intermediate between those of metallic

Fig. 12. Di!erence dC between the speci"c heats of neighbour odd and even systems 1NK 2"101 and 1NK 2"100 4 (multiplied by the ratio of the zero-temperature BCS gap D over ¹) as function of temperature for the same parameters w D+1.1, w K"100 as in Figs. 9}11. The solid and dashed curves correspond to self-consistent calculations with and $ $ without projection on particle-number parity, respectively. The critical temperature for the even system and the two critical temperatures for the odd system show up as discontinuities. The curve di!ers qualitatively from that of the inset, which for 1NK 2"200 and 1NK 2"201 corresponds to the parameters w D"10, w K"200 of Figs. 3}8. In this inset, the $ $ dashed, dashed-dotted and dotted curves are indistinghishable from the temperature axis.

340

R. Balian et al. / Physics Reports 317 (1999) 251}358

mesoscopic islands (Section 6.5.2) and those of ultrasmall aluminium grains (Section 6.5.4). We recall that D denotes the zero-temperature BCS gap associated with the lowest-energy quasiparticle, that is, D,D (¹"0). !1 Our analysis shows that even-N systems are sensitive to the parameter w D at temperatures not $ too far from the BCS critical temperature, and that their projected gap D (¹) is slightly stronger > than the BCS gap D (¹). The discrepancy with BCS is stronger for odd-N systems; as measured !1 by D (¹)/D, it signi"cantly increases when the size (i.e., w D) decreases. The di!erence is especially \ $ important when the temperature is much lower than D, since the p"0 level is completely blocked at ¹"0. Accordingly, for su$ciently low temperatures, the ratio D /D increases with ¹ as \ a consequence of the thermal excitation which depopulates this p"0 level (Figs. 4 and 5). The zero-temperature value of D /D decreases with w D: thus, for "xed D, odd nuclei or odd supercon\ $ ducting metallic grains display a weaker and weaker pairing gap D as their size decreases. When \ w D becomes su$ciently small, the pairing fades out, "rst at ¹"0. Then, for still smaller odd $ systems, the temperature domain over which pairing correlations are e!ective shrinks (Fig. 9) and it eventually disappears at a "nite temperature for odd-N systems. In contrast, for even systems whatever their size, D /D remains equal to 1 at zero temperature. > Of course, as suggested long ago by Anderson [46], superconductivity disappears at any temperature for a "xed coupling constant when the size of the system } directly re#ected on the level spacing } is su$ciently small (the BCS gap itself then vanishes at ¹"0 and hence at all temperatures due to the discreteness of the spectrum). As already pointed out by several authors [24}28] this disappearance (as characterized by the projected gap D ) takes place for larger sizes E when N is odd than when N is even. In both the BCS and N-parity projected results, the sharp transition from the superconducting to the normal phase arises because of the invariance-breaking mean-"eld nature of the approximation. In reality, the "nite size of the system prevents any sharp transition, and more elaborate treatments are required to determine to which extent the discontinuities that were found in the present section are smoothed out. Various methods have been devised to study this broadening. Several e!ects were already taken into account in Ref. [47], including #uctuations of the order parameter in the Ginzburg-Landau theory and quasi-particle contributions associated with discrete spectra. A revival of such ideas can be found in [26], where path integral techniques are used to account for the #uctuations, and in [16,17], where the static path approximation is combined with the projection on the particle-number, or on the N-parity, with applications to nuclear physics. It does not seem easy to merge these techniques with our variational approach so as to render the treatment of Section 6 more realistic around the transition lines. Moreover, some of them do not apply at low temperatures. Already proposed in [47], the projection on a well-de"ned particle number seems better suited to our purpose. One can surmise that it should account, at least partly, for the #uctuations and "nite-size e!ects which smooth out the transition. This conjecture is supported by the work of Ref. [14], where the rounding-o! of the BCS transition produced by the projection on the particle number was quantitatively exhibited on an exactly soluble model. The approximation is based on the minimization of the free energy after projection. However, the use at non-zero temperature of this standard variational principle requires the evaluation of the entropy of the projected independent quasi-particle state. The di$culty of such a calculation led the authors of Ref. [14] to treat only a degenerate BCS Hamiltonian of the form (6.34) with e "0 N

R. Balian et al. / Physics Reports 317 (1999) 251}358

341

for all p (with a constant G ). Another indication may be given in the zero-temperature nuclei by NO the smoothness of the superconducting to normal phase transition which is experimentally observed when the angular momentum increases and theoretically explained through the projection method (see for instance Refs. [12,44]). Returning to our problem and assuming that variational approximations in the canonical ensemble can describe the most important e!ects of the transition broadening, we can imagine to apply the method of Section 5 with the projection on the particle number (Eq. (5.2) rather than on the N-parity as we did in Section 6. This approach seems promising since our variational principle (5.16) avoids the explicit calculation of the trial entropy. It introduces instead the two u-dependent trial quantities DK (u) and AK (u) whose evolution must be determined over the interval [0, b]. The fact that the variation of I is performed after projection should, as suggested by the results of Ref. [14], entail the products u v to be non-zero at any temperature, thus producing N N a smooth transition. The projection on a well-de"ned particle-number should have a negligible incidence on the results obtained for mesoscopic metallic islands where 1NK 2 is large and hence has small relative #uctuations. In nuclei, on the other hand, the smoothing of the superconducting to normal phase transition due to the "nite number of neutrons and protons should be signi"cant; nevertheless the increase of D with ¹ occuring in odd-N systems (Fig. 4) may be real. For very small metallic \ grains, for which the smoothing e!ect is expected to be the largest, it is an open question whether or not the reentrance phenomenon, as it is for instance revealed in Fig. 12 by a maximum in the curve for the di!erence of speci"c heats, will survive.

7. Summary and perspectives A general variational method for evaluating the properties of "nite systems at thermodynamic equilibrium has been proposed and applied to "nite fermion systems with pairing correlations. Extensions of the standard mean-"eld theories were thus obtained for two types of problems: (1) the variational evaluation not only of thermodynamic quantities, but also of expectation values, #uctuations or correlations between observables and more generally of characteristic functions; (2) the restoration of some constraints on the physical space of states or/and of some broken symmetries violated by the variational AnsaK tze. To reach this double goal, we relied on an action-like functional (2.5) which is designed to yield as its stationary value the relevant characteristic function. Two dual trial quantities enter this functional. The operator DK (u) is akin to a statistical operator evolving from 1) to the actual density operator DK "e\@)K

(7.1)

as u goes from 0 to b, while the operator AK (u) is akin to the observable AK de"ned in Eq. (1.6). Such a doubling of the number of variational degrees of freedom is distinctive of a general class of variational principles built to optimize some quantity of interest de"ned by given constraints. Here the constraint is the (symmetrized) Bloch equation (2.4) satis"ed by DK (u), while AK (u) appears as a Lagrangian multiplier associated with this constraint. For unrestricted trial spaces, the stationar-

342

R. Balian et al. / Physics Reports 317 (1999) 251}358

ity of the functional (2.5) expresses the Bloch equation (2.4) for DK (u); in addition, it yields the equation (2.11) for AK (u) which is the Bloch analogue of a (backward) Heisenberg equation. These two equations are decoupled and duplicate each other. However, when the trial spaces are restricted (as is unavoidable for any application), the two equations generally become coupled, with mixed boundary conditions involving both ends of the integration range [0, b]. The fundamental thermodynamic relations are preserved by the variational approximations under rather general conditions on the choice of the trial spaces (Section 2.3); these conditions are satis"ed by the AnsaK tze that have been used throughout the present article. For conserved observables, the characteristic functions, and therefore all the cumulants, are obtained from the approximate partition function via a simple shift of the operator KK (Section 2.4). Moreover, in the limit when the sources m entering the characteristic function vanish, and provided the trial spaces satisfy two A additional conditions stated in Section 2.5, our variational principle reduces to the standard maximization principle for thermodynamic potentials. Beginning with Section 3, we used the functional (2.5) to study systems of fermions in which pairing correlations are e!ective. For the variational operators DK (u) and AK (u) we choose the trial forms (7.2) ¹K B (u),exp(!l B (u)!aRL B (u)a) ,  ¹K ? (u),exp(!l ? (u)!aRL ? (u)a) , (7.3)  where a de"ned in Eq. (3.1) denotes the 2n operators of creation aR or annihilation a ; the trial H H quantities are the scalars l B and l ? and the 2n;2n matrices L B and L ? . The generalized Wick theorem then allows to write explicitly the resulting functional as Eq. (3.32). The stationarity conditions are expressed by the set of coupled equations (3.37), (3.43), (3.46) and (3.47) which constitute an extension of the Hartree}Fock}Bogoliubov approximation. The equation (3.37) proceeds forwards from 0 to b with respect to the variable u while Eq. (3.43) evolves backwards from b to 0, since the boundary conditions are "xed at u"0 for DK and at u"b for AK . In Section 4 thermodynamic quantities, #uctuations and correlations have been evaluated via the expansion of the coupled equations in powers of the sources m . The thermodynamic quantities A are given by the zeroth order term in the expansion while expectation values 1QK 2 of single-quasiA particle observables (7.4) QK "q #aRQ a A A A  are given by the "rst order. In both cases, one recovers (Section 4.1) the standard HFB results as expected from the general discussion of Section 2. The #uctuations and correlations C ,1QK QK #QK Q 2!1QK 21QK 2 are given by the second B A A B AB  A B order. In the case where at least one of the two observables QK and QK commute with the operator A B KK , we obtain (Eq. (4.27)) 1 C " Q : F\ : Q . AB b A  B

(7.5)

The 2n;2n matrices Q and Q de"ned in Eq. (7.4) are regarded in Eq. (7.5) as vectors in the A B Liouville space, and the colon symbol stands as a scalar product in this space, or equivalently as

R. Balian et al. / Physics Reports 317 (1999) 251}358

343

half the trace (tr ) in the original space of the 2n;2n matrices (Section 4.2.2 and Eq. (A.1)). In the L Liouville space, the stability matrix F is de"ned in Eq. (4.20) as the second derivative with respect  to the reduced density R of the HFB grand potential F around its minimum, the HFB solution % R . If both observables QK and QK do not commute with KK , Eq. (7.5) expresses their Kubo  A B correlation (Section 4.3.1). It then appears as a generalization to the non-commutative case of the Ornstein}Zernike formula, here derived from a general variational principle. The ordinary correlation of two observables QK and QK that do not commute with KK is covered A B by the more general formula (Eq. (4.63)) C "Q : (K coth b K )F\ : Q .    B AB  A 

(7.6)

The matrix K is the usual RPA kernel written in the Liouville space (Eq. (4.47) or (A.15)). The HFB  stability matrix F is related to the dynamical RPA matrix K through Eq. (4.49), or (A.22)}(A.23).   Notwithstanding the independent-quasi-particle nature of the operators DK (u) and AK (u), we have obtained non-trivial approximate expressions for the correlations. In particular, formulas (7.5) or (7.6) incorporate long-range e!ects. This was made possible because, through the optimization of the functional (2.5), the variational quantities DK (u) and AK (u) acquire a dependence on the sources m . A Despite the coupling inherent to the stationarity equations (3.37) and (3.43), we have found the explicit expressions (7.5) and (7.6) for the correlations. In this context, we have found a simple u-dependence for the solutions of these equations to "rst order in the sources. Although the "nal formulas (7.5) or (7.6) are compact, their derivation was complicated by the existence of zero eigenvalues in the kernel F , due to some broken invariance, and in the kernel K (Section 4.3.2).   Indeed, for the formulas (7.5) and (7.6) to be meaningful, one must specify how to handle these zero eigenvalues. We showed in Sections 4.2.3 and 4.3.3 that the precise meaning to be given to these formulas depends on the commutation, or non-commutation, of the observables QK and QK with the A B conserved single-particle operator associated with the broken invariance. Analyzed in perturbation theory, the formulas (7.5) and (7.6) were shown to correspond to a summation of bubble diagrams (Section 4.4). Let us stress that the RPA kernel comes out from our variational approach without any additional assumption beyond the choice (7.2)}(7.3) of the trial AnsaK tze. In this way, the RPA acquires a genuine variational status (albeit not in the Rayleigh}Ritz sense). These various aspects of the RPA are discussed in Section 4.5. Sections 5 and 6 are devoted to the second type of problems posed in the Introduction, namely to those situations which require projections. Though the operator KK , and hence the exact density operator (7.1), commute with any conserved quantity, the corresponding invariance can be broken by the variational approximation if the commutation is not preserved by the trial operator DK (u). In this case, a projection PK over an eigenspace of the conserved quantity is required on both sides of DK to restore the symmetry. On the other hand, even when the invariance is not broken, a projection is needed if the trial operators act in a space which is wider than the physical Hilbert space H . (An  example is the extension of the thermal Hartree}Fock approximation to a system with a well de"ned particle number.) Sections 5 and 6 cope with both situations. In the "rst case where the invariance is broken, we replace the trial expression (7.2) for the density operator by DK (u)"PK ¹K B (u)PK ,

(7.7)

344

R. Balian et al. / Physics Reports 317 (1999) 251}358

where ¹K B (u) still has the independent quasi-particle form (7.2); for AK (u) we keep the ansatz (7.3). In order to handle the product (7.7) we take advantage of the fact that the projection PK can be written as a sum of group elements ¹K E which have an independent quasi-particle form, so that Eq. (7.7) is also a sum of terms of the ¹K -type. The resulting functional (Eq. (5.16)) and coupled stationarity equations (Eqs. (5.22), (5.23), (5.27) and (5.28)) are written in Sections 5.2 and 5.3. Although these coupled equations are noticeably more complicated than the unprojected ones, the number of variational parameters, and therefore the dimension of the numerical problem, do not increase. The case of unbroken PK -invariance, where the trial operators ¹K B (u) and ¹K ? (u) commute with the group elements ¹K E , is discussed in Section 5.4. Since PK now commutes with ¹K B (u), the ansatz (7.7) reduces to DK (u)"PK ¹K B (u)"¹K B (u)PK ,

(7.8)

which entails several simpli"cations in the action-like functional (Eq. (5.40)) and in the stationarity equations. We have examined in some detail (Section 5.4.2) the limiting case AK (m)"1) which admits the explicit solution (5.43) for DK (u) and AK (u). This leads to further simpli"cations in the functional (Eq. (5.50)) and in the self-consistent equations (5.51)}(5.54). Several thermodynamic consistency properties were veri"ed. As an example the projected entropy, energy and grand potential satisfy the usual thermodynamic identities (Eqs. (2.26)}(2.29)). In Section 6 the theoretical framework elaborated in Section 5.4 was applied to study the e!ects of the particle-number parity in a "nite superconducting system at thermodynamic equilibrium. In the trial density operator (7.8), PK given by Eq. (5.3) is the projection on either even or odd particle numbers. This problem is of interest for recent or future experiments on mesoscopic superconducting islands, small metallic grains or heavy nuclei. The explicit form of the coupled equations (see Eq. (6.9)) was written in Section 6.1. The parity-projected grand partition function (or more generally the characteristic function for conserved single-particle operators QK ) was evaluated in A Section 6.2, where Eq. (6.23), together with the de"nitions (6.25)}(6.26), was shown to replace the standard HFB self-consistent equation (4.2). Section 6.3 deals with a system (electrons in a superconductor or nucleons in a deformed nucleus) governed by a BCS model Hamiltonian, where the non-interacting part involves twofold degenerate levels p with energies e , and the interaction scatters pairs from one level to another. The BCS N gap is now replaced by the quantities D de"ned in Eq. (6.50), and the usual BCS equation is EN replaced by the number-parity projected self-consistent equation (6.56), 1#gr t\t\ D t 1 N O , EO O (7.9) D " G NO((e !k)#D 1#gr t\ EN 2 N O O EO in which G is a pairing matrix element of the BCS Hamiltonian (6.34). The notation t stands for NO N t ,tanhbe , (7.10) N  N where the quasi-particle energy e , given by Eq. (6.54), di!ers from the standard BCS expression N ((e !k)#D . The number r is de"ned (Eq. (6.46)) as N EN  r ,“ t .  N N

(7.11)

R. Balian et al. / Physics Reports 317 (1999) 251}358

345

The formulae (7.9)}(7.11) depend on the parity of the particle number N through the value of g (g"#1 for even N or g"!1 for odd N). The e!ects of the projection appear indirectly through the occurence of e in Eq. (7.10), and directly through the last fraction of Eq. (7.9). N As discussed in Section 6.3.3, this fraction is larger than one for even-N systems, smaller for odd-N systems. As a consequence, the even-N projected gap is larger than the BCS gap, while the odd-N projected gap is smaller, in agreement with the idea that pairing is inhibited in small odd-N systems. The analysis of this fraction reveals, moreover, that the odd-N projection di!ers from BCS more than the even-N. In Section 6.4 we considered some limiting situations. At zero temperature the BCS and even-N projected gaps are equal. However, while the BCS quantities deviate from their zero-temperature limit by terms of the form exp(!bD), the deviations are of the form exp(!2bD) for the even-N projection. The di!erences are more striking between BCS and the odd-N projection. At zero temperature, the BCS gap is larger than the odd-N projected gap. When the temperature grows, the latter starts to increase before reaching a maximum. The odd-N entropy tends to ln 2 when ¹ tends to zero, re#ecting the twofold degeneracy of the ground state; in Section 6.5.4 we commented upon the fact that this result is not completely trivial. In Section 6.4.2, we had shown that, at temperatures larger than the single-particle level spacing and smaller than the pairing gap, the odd-N entropy coincides with the Sackur-Tetrode entropy of a single-particle with a mass D moving in a one-dimensional box. This is consistent with gapless elementary excitations in odd-N systems. When the single-particle level density is su$ciently large, the critical temperatures are equal for the BCS and the odd or even projected solutions. Moreover, at non-zero temperature, for a large system or when the level spacing becomes small with respect to the pairing gap at the Fermi surface, Eq. (7.9) reduces to the usual BCS equation. However, the value of the projected (odd or even) entropy is lower than the BCS value by ln 2. Note that the entropy is not involved in the action functional (2.5), in contrast with the standard variational principle for thermodynamic potentials. Here, the entropy comes out as a by-product of our variational approximation; it is obtained in the end through the di!erentiation (2.28). In order to push further the comparison between the BCS and the projected results we presented in Section 6.5 schematic calculations which exemplify three di!erent physical situations. A crucial parameter turns out to be the product w D where D is the zero-temperature BCS gap and w the $ $ single-particle level density at the Fermi surface. This parameter, which indicates how many single-particle levels lie within the energy range of the gap, is sometimes called the `number of Cooper pairsa. When w D is large, the projection e!ects are weak. Fig. 5 illustrates the case w D"10. The BCS $ $ and even-N projected gaps are almost undistinguishable. The relative di!erence between the BCS and odd-N gaps is 1/2w D at ¹"0. The di!erence between the free energies of the odd and even $ systems can be estimated perturbatively from the BCS solution. The outcome (Eq. (6.99)), equal to D!1/4w at ¹"0, decreases almost linearly with ¹ and displays a cross-over towards the $ asymptotic value !1/4w at a temperature ¹ smaller than the critical temperature (see Fig. 8). $ " Between ¹ and ¹ , the gap function is the same whatever the N-parity. These results are "  consistent with experiments on mesoscopic superconducting islands which show that, when w D1, the odd}even e!ects disappear at some temperature below ¹ [45]. $  The value w DK2 corresponds roughly to the case of heavy nuclei (Fig. 4). In the even system, $ the di!erence between the projected and BCS gaps becomes more appreciable than in Fig. 5. In the

346

R. Balian et al. / Physics Reports 317 (1999) 251}358

odd system, the projected gap at zero temperature is reduced by about 30% with respect to BCS; this result can also be obtained from a BCS calculation in which the pair formation is blocked for the level that coincides with the Fermi surface. As the temperature grows, so does the odd-N projected gap. The BCS, even and odd projected gaps merge again at higher temperature where they all display the same critical behaviour. Fig. 9, where w DK1.1, corresponds to a situation which could be encountered in ultrasmall $ Aluminium grains [4}6]. In the even-N projected state, the critical temperature is higher than for the BCS state. In the odd system a new phenomenon occurs: there is no gap at zero temperature. However, the phase diagram (¹, D) obtained by projecting on odd number exhibits a reentrance ewect: when the temperature increases, a gap appears at a "rst transition temperature, reaches a maximum and disappears at a second transition temperature smaller than the BCS one. This behaviour can be understood as follows: at ¹"0, at least one particle remains unpaired and it fully occupies one of the two single-particle states at the Fermi surface, forbidding the formation of a pair with components on this strategic level. When the temperature increases, the unpaired particle partially frees up this level, which then becomes available for the formation of Cooper pairs. According to Fig. 12, the existence of two transitions in odd ultrasmall metallic grains could be experimentally detected by speci"c heat measurements of odd and even systems. However, in more elaborate treatments going beyond the projection on the N-parity, the transitions will be smoothed out by additional "nite-size e!ects [47]; further studies are required to ensure that the reentrance phenomenon is not completely washed away by these e!ects. Besides its relevance for "nite superconducting systems, one can imagine using the variational method introduced in this paper for studying other physical problems. In particular the formulation of Section 5.4, which treats situations with unbroken PK -invariance, is applicable in various problems. Consider for instance a "nite fermion system at equilibrium without pairing correlations. In this case the exponents in the trial quantities (7.2) and (7.3) need only contain single-particle terms in aRa which commute with NK . The projection PK has the form (5.2). G H A variational method to evaluate the canonical partition function is then furnished by Section 5.4.2, which constitutes a projected extension (with a well-de"ned number of particles) of the thermal Hartree}Fock approximation. It su$ces in Eqs. (5.43)}(5.44) to take a trial operator KK of the  single-particle type. This procedure also provides, along the lines of Section 2.4, a consistent variational estimate for the characteristic function of any conserved single-particle observable (commuting with both NK and HK ). For other single-particle observables (commuting with NK only) one should resort to the method of Section 5.4.1 where the trial operators ¹K B (u) and ¹K ? (u) have a more complicate u-dependence. Likewise, for "nite systems without odd-multipole deformations, the projection over a given space-parity can be performed by implementing the operator (5.6) in Sections 5.4.1 and 5.4.2. In case the symmetry is broken (Sections 5.2 and 5.3), the functional (5.16) takes a more intricate form, due to the occurence of two projections. One could then make use of a further approximation by using the scheme described in Ref. [48] where an integral such as (5.2) on the group elements is replaced by a truncated sum. This would generate approximate expressions for the canonical thermodynamic functions and correlations which should be easier to handle numerically than the general expressions (5.23)}(5.25) and (5.28)}(5.29). In the variational methods described above, one could imagine implementing trial quantities DK (u) and AK (u) with a more general form than Eqs. (7.2) and(7.3). For instance Ref. [49] has

R. Balian et al. / Physics Reports 317 (1999) 251}358

347

demonstrated, in the related context of time-dependent Hartree}Fock equations, the feasibility of extensions where T B (u) is multiplied by some polynomial in creation and annihilation operators. Using the methods developed here, one could obtain extensions of the HFB formalism, or of the projected HFB formalism, which would take into account a wider class of correlations. The variational setting of Section 2 could easily be adapted to problems of classical statistical mechanics. A "rst possibility would be to start from the classical version of the functional (2.5) in which traces, density operators and observables are replaced by integrals, distribution functions and random variables in phase space. The commutative nature of these quantities implies that a characteristic function would be obtained as a maximum. Alternatively, one could start from the quantal functional (2.5) and then apply a semi-classical limiting process, for instance through the use of the Wigner transform. The latter procedure would yield in addition the "rst quantum corrections to the classical limit. This could be of interest for heavy nuclei or for su$ciently large clusters of atoms and molecules. Minor modi"cations are needed to deal with systems of interacting bosons. Besides the change of the anticommutation rules (3.2) into commutation ones, linear terms in a and aR should be added to the observables (3.12) and to the exponents of the trial operators (7.2) and (7.3). This should result in a systematic variational approach of Bose}Einstein condensation in "nite interacting systems, a phenomenon which has become of direct experimental interest with the recent achievement of condensates in dilute atomic gases at ultra-cold temperatures. It is worthwhile to note that here again the projection on the right number of particles is an important requirement. Another physical problem where a projection is essential arises in particle physics, when mean-"eld-like approximations violate the colour invariance. One may imagine restoring it by means of a projection on colour singlets along the lines of Section 5. In conclusion we would like to stress the consistent character of the present method, besides its generality and its #exibility. We have veri"ed the agreement of our approximation scheme with thermodynamics. We have recovered naturally the RPA (the complete series of bubble diagrams). Finally, although we have not touched upon the subject in this already long article, we point out the complementarity of the present static approach with the variational treatment of dynamical problems of Ref. [29]. Indeed, the use of the Bloch equation as a constraint for the initial state combines coherently with the use of the (backward) Heisenberg equation as a dynamical constraint; a variational principle comes out which (within restricted trial spaces) optimizes both the initial state and the evolution of the system. This comprehensive static and dynamic variational method should provide a convenient tool for evaluating time-dependent correlations, and hence transport properties, either in the grand-canonical formalism as in Section 4 or with projections as in Sections 5 and 6.

Acknowledgements We are indebted to M.H. Devoret, D. Este`ve and N. Pavlo! for very instructive discussions, to S. Creagh and N. Whelan for reading parts of our manuscript, and to M.T. Commault and M. Trehin for their help in the preparation of the "gures. One of us (M.V.) gratefully acknowledges the support of the Alexander von Humbolt-Stiftung and the hospitality of the Physics Department of the Technische UniversitaK t MuK nchen.

348

R. Balian et al. / Physics Reports 317 (1999) 251}358

Appendix A. Geometric features of the HFB theory A.1. The reduced Liouville space The Liouville representation of quantum mechanics (see for instance [51], Section 3) relies on the idea that the algebra of observables, when regarded as a vector space, can be spanned by a basis of operators. Any observable is then written as a linear combination of the operators of the chosen basis, with coe$cients interpreted as (covariant) coordinates. Correlatively, for any state, the expectation values of the operators of the basis are regarded as the (contravariant) coordinates of this state. In the HFB context, it is su$cient to consider the subalgebra of quadratic forms of creation and annihilation operators a (equal to either a and aR). With the notation (3.12), keeping H G G aside the trivial c-number q, we are led to regard the coe$cients Q as the covariant coordinates of HHY the operator  aRQ a . The factor  and the constraint (3.5) imposed on Q are consistent with   HHY H HHY HY the duplication aRa "!a aR#d "! p aRa p #d occuring in the operator basis. H HY HY H HHY IIY HYI I IY IYH HHY We shall thus consider Q both as a matrix in the 2n-dimensional j-space of creation and HHY annihilation operators (Section 3.1) and as a vector in the 2n;2n-dimensional reduced Liouville space of observables. Likewise, since in the HFB approximation we focus on density operators of the form (3.4) and since they are characterized by the contraction matrices R de"ned by Eqs. (3.8) and(3.10)}(3.11), HHY the contravariant coordinates of these density operators reduce to the set R , with the constraint HHY (3.9). Again, R is both a matrix in the j-space and a vector in the reduced Liouville space of states. HHY As in the full space of observables and states, the expectation value of an operator aRQa, in  a state characterized by its contractions R , is given here by the scalar product of the vectors HHY R and Q in the Liouville space. We denote this scalar product as (A.1) Q : R"tr QR" Q R .  HHY HYH  L HHY In Eq. (A.1) the colon symbol : indicates therefore both a summation on a twisted pair of indices and a multiplication by the factor  accounting for the duplication of coordinates both in R and Q.  Among others, we shall use in this Appendix the following two matrices in Liouville space: "2d d , R "2p p . (A.2) HHYIIY HIY IHY HHYIIY HI IYHY They are invariant under the exchange (jj)  (kk) (symmetry in Liouville space). One has also R"1. Moreover, the relation (3.9) implies that the variation dR of a contraction matrix satis"es 1

R : dR"pdR2p"!dR .

(A.3)

A.2. Expansion of the HFB energy, entropy and grand potential The HFB energy, de"ned by E+R,"Tr DK KK /Tr DK (including therefore the term !k1NK 2), is given by Wick's theorem as the function (3.31) of the Liouville vector R . Its partial derivatives HHY have been de"ned in Eqs. (3.38)}(3.40) as the 2n;2n matrix H+R, constrained by (3.41). In the Liouville space, in agreement with Eq. (3.38), H appears as the gradient of E, that is, a covariant HHY

R. Balian et al. / Physics Reports 317 (1999) 251}358

349

vector which, with the notation (A.1), reads dE+R,"H+R, : dR .

(A.4)

Using Eqs. (3.7) and (3.8), the HFB entropy is evaluated as S+R,,!Tr DK ln DK /Tr DK !ln Tr DK "!tr [R l n R#(1!R) ln (1!R)] , (A.5)  L where the two terms under the trace are equal. Similarly to Eq. (A.4), its gradient, written by means of Eq. (3.9) so as to satisfy the same identity (3.41) as H, is given by dS+R,"L : dR ,

(A.6)

with L,ln [(1!R)/R]. The second derivatives of E+R, and S+R, de"ne the twice covariant tensors E and S entering the expansions E+R#dR,"E+R,#H+R, : dR#dR : E : dR , (A.7)  S+R#dR,"S+R,#L : dR#dR : S : dR#2 . (A.8)  Both matrices E and S are symmetric (for instance, E "E ). Moreover, they satisfy HHYIIY IIYHHY R : E : R"E , R : S : R"S , (A.9) so that E : dR obeys the same symmetry relation (A.3) as dR. The explicit form of the second derivative E of E+R, comes out from Eq. (3.31) and dH"E : dR. When R and dR are HHY HHY hermitian matrices in the j-space, the matrix S of the second variations of S+R, is always negative in Eq. (A.8). As shown in [51], it can be expressed as

 

S "!2 HHYIIY



dv



1 R#v(1!R)



1 R#v(1!R)

,

(A.10)

HIY IHY  which, in a basis where R is diagonal (R "R d ), reduces to HHY H HHY L !L H I. (A.11) S "2d d HHYIIY HIY IHY R !R H I The ratio is equal to !1/R (1!R ) when R "R . H H I H The standard HFB variational principle (Section 4.1) provides the grand potential at grand canonical equilibrium as the minimum of 1 F +R,"E+R,! S+R, , % b

(A.12)

which is attained for R"R . From Eqs. (A.4) and (A.6), the stationarity condition yields the  self-consistent equation 1 H+R ," L ,  b 

(A.13)

350

R. Balian et al. / Physics Reports 317 (1999) 251}358

which de"nes the HFB solution R , in agreement with Eq. (4.1). Around this equilibrium, the  expansion of the grand potential F +R, reads % (A.14) F +R #dR,"F +R ,#dR : F : dR#2 ,  %  %   where F ,E !S is a positive matrix, like !S .    @  A.3. The HFB grand potential and the RPA equation The expansion of Eq. (3.52) around R ?B "R (or around R B? "R ) leads to Eq. (4.46) which   involve the RPA kernel K de"ned in Eq. (4.47) as  RH+R, dR , R . (A.15) K : dR,[H+R ,, dR]# IIY    RR IIY RR IIY In Liouville space the matrix K is nonsymmetric; it depends on R and on the parameters of KK .   We show here that K is directly related to the matrix F which enters Eq. (A.14), or Eq. (4.20).   To this aim, we "rst introduce in the Liouville space a new matrix which generates the commutators with R. Its action on any vector M is de"ned by





C : M"!M : C"[R, M] .



(A.16)

Accordingly, its expression as a matrix is C "d R !d R #p (Rp) !p (pR) , (A.17) HHYIIY IHY HIY HIY IHY HYIY IH IH HYIY where the last two terms yield the same contribution as the "rst two when applied to a vector M satisfying Eq. (A.3) (i.e., R : M"!M). The matrix C is antisymmetric, and it moreover satis"es C : R"R : C"!C .

(A.18)

In a basis where R is diagonal, it reduces to C "2d d (R !R ) , (A.19) HHYIIY HIY IHY H I as it is obvious from Eq. (A.16). The commutators with the matrix L can then be represented in Liouville space by the product C : S which acts according to C : S : dR"S : C : dR"!dR : C : S"[L, dR] .

(A.20)

Indeed, in a basis where R is diagonal, this identity, analogous to Eq. (A.16), is a direct consequence of Eqs. (A.11) and (A.19). Note that the matrices C and S commute. The "rst term of Eq. (A.15) is thus proportional to Eq. (A.20), in which C and S have to be taken at the stationary point where the relation L ,ln (1!R )/R "bH+R , is satis"ed. To write the     second term, we note that RH dR "E : dR , IIY RR IIY IIY

(A.21)

R. Balian et al. / Physics Reports 317 (1999) 251}358

351

a consequence of Eq. (A.7). Hence, using Eq. (A.16) and then adding Eq. (A.20), we "nd for any R and dR:





RH+R, 1 1 dR , R # [L, dR]"!C : E : dR# C : S : dR IIY RR b b IIY IIY "!C : F : dR .

(A.22)

The sought for expression of the kernel (A.15) is then simply K "!C : F , (A.23)    where the matrices C and F are evaluated at the equilibrium point R"R . The RPA kernel K is     thus directly related to the matrix F of the second derivatives of the grand potential.  At the minimum (attained at R ) of the trial grand potential F +R,, the matrix F is non %  negative. Therefore the eigenvalues of K , which are the same as those of !F : C : F in the     space of Liouville vectors satisfying Eq. (A.3), are real. We thus recover the well-known consistency between the stability of the variational HFB state and the nature of the associated RPA dynamics in real time. Moreover, due to the antisymmetry of C , the non-vanishing eigenvalues of K come in   opposite pairs. A.4. Riemannian structure of the HFB theory By relying on the very existence of the von Neumann entropy S(DK ), it has been shown that the set of density operators DK can be endowed with a natural Riemannian structure [51]. The basic idea is to introduce a metric such that the distance ds between two neigbouring states DK and DK #dDK is de"ned by ds"!dS .

(A.24)

With this metric the general projection method of Nakajima}Zwanzig appears as an orthogonal projection in the space of states. A reduction of the exact description through projection on a subset of simpler states, in our case the HFB states of the form (3.4), induces from Eq. (A.24) a natural metric on this subset. It is thus natural to de"ne !S, the positive symmetric matrix describing the second variations of the entropy (A.8), as a metric tensor in the reduced 2n;2n Liouville space of states. This de"nition allows us to introduce distances between neighbouring HFB states, and also to establish a correspondance between the covariant and contravariant components of the vectors of Section A.1. More precisely, this canonical mapping relates the variations of the observables aRLa to the variations of those states which are their exponentials; it reads:  S : dR"dL . (A.25) We can also introduce the twice contravariant metric tensor !S\ in the subspace of single quasi-particle observables by inverting S in Liouville space (S\ : S"1). In a basis where R and L are diagonal, the matrix elements of S\ are R !R H I , S\ "2d d HHYIIY HIY IHY L !L H I

(A.26)

352

R. Balian et al. / Physics Reports 317 (1999) 251}358

where the ratio is !R (1!R ) when R "R . In an arbitrary basis, it can be written as H H I H L L  eS e\S S\ "!2 du L . (A.27) HHYIIY e #1 eL#1 HIY IHY  Finally, we can introduce as in Eq. (3.7) the reduced thermodynamic potential

 









1 1 ln f+L,"ln Tr exp ! aRLa " tr ln (1#e\L) . 2 2 L

(A.28)

The expression of the "rst-order term in the expansion in powers of dL, 1 ln f+L#dL,"ln f+L,!R : dL! dL : S\ : dL#2 , 2

(A.29)

exhibits the fact that the entropy (A.8) is the Legendre transform of the reduced grand potential (A.28), while the second-order term is seen to describe a distance between two neighbouring observables. This last feature is the counterpart in the reduced HFB description of a structure already recognized [51] in the full Liouville space. All these properties, on which we shall not dwell here, confer a geometric nature to the reduced HFB description and set it into the general projection method [51]. A.5. Lie}Poisson structure of the time-dependent HFB theory We shall now introduce another geometric type of structure, of a symplectic type. This will allow us to regard the reduced energy E+R, or the grand potential F +R, as a classical Hamiltonian for % the dynamics associated with the HFB theory. The evolution of the variables R is governed in the time-dependent mean-xeld HFB theory by HHY the equation i

dR "[H+R,, R] , dt

(A.30)

written in the 2n;2n matrix representation. In order to analyze this equation, we "rst note that the basic operators aaR of the Liouville representation are the generators of a Lie group with the algebra [a aR , a aR ]"d a aR !d a aR #p a a !p aR aR . (A.31) H HY I IY IHY H IY HIY I HY HYIY I H IH HY IY Following the arguments of Ref. [50], where a similar formulation was set up for the Hartree}Fock case, we note that the contractions R are through Eq. (3.8) in one-to-one correspondence with HHY the operators a aR . Regarding the R 's as a set of classical dynamical variables, we introduce H HY HHY among them a Poisson structure; indeed, the Lie-algebra relation (A.31) suggests to characterize the Lie}Poisson bracket between the dynamical variables by the Poisson tensor i +R , R ,"C , (A.32) HHY IIY HHYIIY where each C is the linear function of the R's inferred from the right-hand side of Eq. (A.31). As a matrix, C is readily identi"ed with the matrix introduced in Eqs. (A.16), (A.17) which described, in

R. Balian et al. / Physics Reports 317 (1999) 251}358

353

the Liouville space, the commutation with R. While the metric of Section A.4 was characterized by the symmetric Riemann tensor !S, the present Poisson structure is characterized by the antisymmetric Poisson tensor C. For more details on Poisson structures, see for instance Ref. [52]. In contrast to the standard Poisson brackets between canonically conjugate position and momentum variables, which are c-numbers, the Lie}Poisson brackets (A.32) depend on the dynamical variables R . In the Liouville representation, the bracket of two functions f and g of the HHY R's is evaluated by saturating their gradients with the tensor C according to Rg Rf . i + f, g," : C : RR RR

(A.33)

The compatibility of Eq. (A.33) with the de"nition (A.32) results from the property (A.18) of the tensor C and from the relation 1 RR HHY" (1!R) , (A.34) HHYIIY 2 RR IIY which itself is a consequence of Eq. (3.8). Using the rules (A.33) and (A.34), we identify the r.h.s. of Eq. (A.30) as a Lie}Poisson bracket since the time-dependent mean-"eld HFB equation reads dR 1 " [H+R,, R] i

dt 1 RE 1 RE 1 RE RR " : C" : C : (1!R)" :C: , i RR 2i RR i RR RR

(A.35)

that is dR "+E+R,, R, . dt

(A.36)

This equation has thus the form of a classical dynamical equation of motion; the HFB energy E+R, plays the ro( le of a classical Hamiltonian, while a Lie}Poisson structure is associated with the dynamical variables R by Eq. (A.32). Since the gradient (A.6) of S+R,, as a matrix in the j-space, commutes with R, the relation (A.16) implies that the Lie}Poisson bracket +S+R,, R, vanishes. Hence, we can alternatively regard the grand potential (A.12), instead of E+R,, as a classical Hamiltonian and write dR "+F +R,, R, . % dt

(A.37)

This remark is useful when the time-dependent mean-"eld equation (A.30) is used in the small amplitude limit to describe approximately collective excitations around the grand-canonical equilibrium state. By expanding R as R #dR, Eq. (A.30) then reduces to the RPA equation  associated with the HFB approximation, ddR 1 " K : dR , dt i 

(A.38)

354

R. Balian et al. / Physics Reports 317 (1999) 251}358

where the RPA kernel K is expressed by Eq. (A.15). We can now identify the r.h.s. of Eq. (A.38) with  the "rst-order term +dF +R,, R , in the expansion of Eq. (A.37) in powers of dR around R . Using %   the expansion (A.14) and the expression (A.32) of the Lie}Poisson bracket, we "nd that Eq. (A.38) is equivalent to i ddR/dt"!C : F : dR , (A.39)   and thus recover the expression (A.23) for K . This new derivation of Eq. (A.23) provides us with  a simple interpretation of the RPA equation (A.39) as an equation for small motions in ordinary classical dynamics. Indeed, small changes of canonically conjugate variables q , p around a minG G imum of a Hamiltonian H+q, p, are governed by !d.& +q, q, +q, p, .& .& dq d dq .N K! .O.O .O.N " , (A.40) dt dp d.& +p, q, +p, p, .& .& dp .N .N.O .N.N where one recognizes the same multiplicative structure as in the r.h.s. of Eq. (A.39), within the replacement of the ordinary Poisson bracket by Eq. (A.32) and of H+q, p, by F +R,. % For small deviations around an arbitrary TDHFB trajectory, the time-dependent RPA kernel K is de"ned by
















 

RH+R, ddR "K : dR,[H+R,, dR]# dR , R . (A.41) IIY RR dt IIY IIY Since we no longer lie at the minimum of F +R,, this variation (A.41) involves not only the % expansion of F but also of R, which appears both directly and through the Poisson tensor C. We % thus get, in the product i

K : dR"H : dC!C : E : dR "!C : F : dR#(H!b\L) : dC ,

(A.42)

an additional term proportional to the deviation between bH+R, and L"ln[(1!R)/R]. Apart from this term, which arises from the dependence of the structure constants C on the variables R, the present classical dynamics di!ers from an ordinary Hamiltonian dynamics through the presence of vanishing eigenvalues in the tensor C. The Poisson structure generated by (A.32) therefore di!ers from the usual symplectic structure associated with canonically conjugate variables because some combinations of the dynamical variables R have a vanishing bracket with any variable. This property re#ects the existence of structural constants of the motion, which never vary whatever the ewective Hamiltonian in Eq. (A.36). In particular, the eigenvalues of R and S+R, always remain constant.

Appendix B. Liouville formulation of the projected 5nite-temperature HFB equations Using the Liouville-space formalism introduced in Appendix A, it is possible to recast the coupled equations (5.23) and(5.28) governing the evolution of T B (u) and T ? (u), de"ned in Eqs. (3.6) and(3.16)}(3.17), so as to make these equations formally similar to Eqs. (3.37) and (3.43).

R. Balian et al. / Physics Reports 317 (1999) 251}358

355

To this aim, we "rst introduce two operators T ? and T B whose action on a Liouville vector Q is given by T ? : Q,T ? QT ? \ ,

T B : Q,T B QT B \ .

(B.1)

They satisfy the relations R : T ? : R"T ? , T ? \"T ? 2,

R : T B : R"T B , T B \"T B 2 ,

(B.2) (B.3)

where the superscript T denotes the matrix transposition in the Liouville space. Next, we de"ne the matrix



MM B ? , HHYIIY

ZE EY+2([1!R EBEY? ]T E ) (T E \R EBEY? ) HIY IHY E EY #(R EBEY? !RM B ? ) (R BEY?E !RM B ? ) , , HHY IIY and the two Liouville vectors



GM B ? ,

E EY

(B.4)

ZE EY+(1!R EBEY? )H+R EBEY? ,R EBEY? #(E+R EBEY? ,!EM B?)(R EBEY? !RM B ? ), , (B.5)



GM ? B ,

ZE EY+(R ?EBEY H+R ?EBEY ,(1!R ?EBEY )#(E+R ?EBEY ,!EM ?B)(R ?EBEY !RM ? B ), .

E EY After multiplication to the right by T ? , Eq. (5.23) can be reformulated as MM B ? :






dT B

1 T B \ # (GM B ? #T ? \ : GM ? B )"0 . du 2

(B.6)

By introducing in the Liouville space the matrix MM ? B ,T ? : MM B ? : T B ,

(B.7)

namely



MM ? B , HHYIIY

ZE EY+2(R ?EBEY T EY \) (T EY [1!R ?EBEY ]) HIY IHY E EY #(R ?EBEY !RM ? B ) (R EY ?EB !RM ? B ) , , HHY IIY

(B.8)

and the vectors H M B ? "MM B ? \ : GM B ? ,

H M ? B "MM ? B \ : GM ? B ,

(B.9)

we can rewrite Eq. (B.6) as 1 dT B

"! [H M B ? T B #T B H M ? B ] , 2 du a form which resembles closely Eq. (3.37).

(B.10)

356

R. Balian et al. / Physics Reports 317 (1999) 251}358

In the same way Eq. (5.28) can be transformed into dT ? 1 " [H M ? B T ? #T ? H M B ? ] , 2 du

(B.11)

where the matrices H M ? B and H M B ? in the single-particle space are obtained from Eqs. (B.4)}(B.5) and(B.8)}(B.9) by the exchanges a  d and g  g. One major di!erence, however, between the coupled equations (B.10)}(B.11) and Eqs. (3.37) and (3.43) is that the former involve four di!erent quantities (H M B ? , H M ? B , H M ? B , H M B ? ) instead of two (H+R B? ,, H+R ?@ ,) for the latter. The matrices MM ? B and MM B ? required for the calculation of the quantities H M ? B and H M B ? appearing in Eq. (B.11) are related to MM ? B and MM B ? through MM ? B "MM ? B 2 ,

MM B ? "MM B ? 2 .

(B.12)

From the de"nition (B.4)}(B.5) and (B.7) of the matrices MM B ? and MM ? B and of the Liouville vectors GM B ? and GM ? B , one establishes the relations R : MM B ? : R"MM B ? , R : GM B ? "!GM B ? ,

R : MM ? B : R"MM ? B , R : GM ? B "!GM ? B .

(B.13) (B.14)

One deduces that, as the reduced HFB Hamiltonian H+R,, the vectors in Liouville space H M B ? and H M ? B verify the relation (3.41) which in the Liouville formulation becomes: R:H M B ? "!H M B ? ,

R:H M ? B "!H M ? B .

(B.15)

Through the exchanges a  d and g  g one extends the results (B.13)}(B.15) to the matrices MM ? B , MM B ? and vectors GM ? B , GM B ? , H M ? B , H M B ? . When the operator AK (de"ned in Eq. (1.6)) is hermitian, Eqs. (B.10), (B.11), (5.23) and(5.28), have hermitian solutions. Indeed, in such a case these equations preserve the hermiticity of T ? and T B and the reality of l ? and l B . To prove it, one "rst notices that the relations T E R"T E\ , T ? R"T ? , T B R"T B result in R EBEY? R"R ?EY\BE\ ,

(B.16)

ZE EYH"ZEY\E\ .

(B.17)

Then, since ( )H" \ \, one sees from Eq. (B.17) that Z is real. The use of this property along E EY EY E with the equality E+RR,"E+R,H in Eqs. (5.19) and (5.20) implies RM B ? R"RM ? B ,

RM B ? R"RM ? B ,

EM B?H"EM ?B .

(B.18) (B.19)

Finally, on formulae (B.4)}(B.5) and (B.8), one readily checks the relations MM B ? H "MM ? B

, MM B ? H "MM ? B

, HYHIYI HHYIIY HYHIYI HHYIIY GM B ? R"GM ? B , GM B ? R"GM ? B .

(B.20) (B.21)

They yield the equations H M B ? R"H M ? B ,

H M B ? R"H M ? B ,

which entail the hermiticity of dT B /du and dT ? /du.

(B.22)

R. Balian et al. / Physics Reports 317 (1999) 251}358

357

References [1] T.P. Martin (Ed.), Large Clusters of Atoms and Molecules, NATO ASI Series E, vol. 313, Kluwer Academic Publishers, Dordrecht, 1996. [2] H. Bouchiat, in: E. Akkermans, G. Montambaux, J.-L. Pichard, J. Zinn-Justin (Eds.), Proc. 1994 Les Houches Summer School on Mesoscopic Quantum Physics, North-Holland, Amsterdam, 1995, pp. 99}146. [3] M.H. Devoret, D. Esteve, C. Urbina, in: E. Akkermans, G. Montambaux, J.-L. Pichard, J. Zinn-Justin (Eds.), Proc. 1994 Les Houches Summer School on Mesoscopic Quantum Physics, North-Holland, Amsterdam, 1995, pp. 605}658. [4] D.C. Ralph, C.T. Black, M. Tinkham, Phys. Rev. Lett. 74 (1995) 3241. [5] C.T. Black, D.C. Ralph, M. Tinkham, Phys. Rev. Lett. 76 (1996) 688. [6] D.C. Ralph, C.T. Black, M. Tinkham, Phys. Rev. Lett. 78 (1997) 4087. [7] A. Gri$n, D.W. Snoke, S. Stringari (Eds.), Bose}Einstein Condensation, Cambridge University Press, Cambridge, UK, 1995. [8] R. Balian, From Microphysics to Macrophysics, vol. I, Chapt. 4, Springer, Berlin, 1991. [9] E.T. Jaynes, Phys. Rev. 106 (1957) 620; 108 (1957) 171. [10] R. Balian, N. Balazs, Ann. Phys. (NY) 179 (1987) 97. [11] J. Bardeen, L.N. Cooper, J.R. Schrie!er, Phys. Rev. 106 (1957) 162. [12] P. Ring, P. Schuck, The Nuclear Many-Body Problem, Springer, New York, 1980. [13] K. Tanabe, K. Sugarawa-Tanabe, H.J. Mang, Nucl. Phys. A 357 (1981) 20. [14] C. Essebag, J.L. Egido, Nucl. Phys. A 552 (1993) 205. [15] R. Rossignoli, P. Ring, Ann. Phys. (NY) 235 (1994) 350. [16] R. Rossignoli, N. Canosa, J.L. Egido, Nucl. Phys. A 605 (1996) 1. [17] R. Rossignoli, N. Canosa, P. Ring, Phys. Rev. Lett. 80 (1998) 1853. [18] K. Langanke, D.J. Dean, P.B. Radha, S.E. Koonin, Nucl. Phys. A 602 (1996) 244. [19] G. Puddu, Nuovo Cimento 111 A (1998) 311. [20] M.T. Tuominen, J.M. Hergenrother, T.S. Tighe, M. Tinkham, Phys. Rev. Lett. 69 (1992) 1997. [21] D.V. Averin, Yu.V. Nazarov, Phys. Rev. Lett. 69 (1992) 1993. [22] B. JankoH , A. Smith, V. Ambegaokar, Phys. Rev. B 50 (1994) 1152. [23] D.S. Golubev, A.D. Zaikin, Phys. Lett. A 195 (1994) 380. [24] J. von Delft, A.D. Zaikin, D.S. Golubev, W. Tichy, Phys. Rev. Lett. 77 (1996) 3189. [25] R.A. Smith, V. Ambegaokar, Phys. Rev. Lett. 77 (1996) 4962. [26] K.A. Matveev, A.I. Larkin, Phys. Rev. Lett. 78 (1997) 3749. [27] F. Braun, J. von Delft, D.C. Ralph, M. Tinkham, Phys. Rev. Lett. 79 (1997) 921. [28] A. Mastellone, G. Falci, R. Fazio, Phys. Rev. Lett. 80 (1998) 4542. [29] R. Balian, M. VeH neH roni, Nucl. Phys. B 408 (1993) 445. [30] E. Gerjuoy, A.R.P. Rau, L. Spruch, Rev. Mod. Phys. 55 (1983) 725. [31] R. Balian, M. VeH neH roni, Ann. Phys. (NY) 187 (1988) 29. [32] B.A. Lippmann, J. Schwinger, Phys. Rev. 79 (1950) 469. [33] N.N. Bogolyubov, Uspekhi Fiz. Nauk. 67 (1959) 549 [translation: Sov. Phys. Usp. 67 (2) (1959) 236] [34] P.A. Whitlock, M.H. Kalos, G.V. Chester, J. Stat. Phys. 32 (1983) 389. [35] R. Balian, E. BreH zin, Nuovo Cimento B 64 (1969) 37. [36] R. Balian, M. VeH neH roni, Ann. Phys. (NY) 164 (1985) 334. [37] R. Balian, M. VeH neH roni, Ann. Phys. (NY) 216 (1992) 351. [38] E. BreH zin, J.C. Le Guillou, J. Zinn-Justin, in: C. Domb, M.S. Green (Eds.), Phase Transitions and Critical Phenomena, vol. 6, Academic Press, London, 1976. [39] T. Troudet, N. Vinh Mau, D. Vautherin, Nucl. Phys. A 458 (1986) 460. [40] B.F. Bayman, Nucl. Phys. 15 (1960) 33. [41] A.K. Kerman, R.D. Lawson, M.H. Macfarlane, Phys. Rev. 124 (1961) 162. [42] K. Dietrich, H.J. Mang, J.H. Pradal, Phys. Rev. 135 (1964) B22. [43] D.V. Averin, Yu.V. Nazarov, Physica B 203 (1994) 310.

358 [44] [45] [46] [47] [48] [49] [50] [51] [52]

R. Balian et al. / Physics Reports 317 (1999) 251}358 B. Gall, H. Flocard, P. Bonche, J. Dobaczewski, P.H. Heenen, Zeit. Phys. A 348 (1994) 89. P. Lafarge, P. Joyez, D. Esteve, C. Urbina, M.H. Devoret, Letters to Nature 365 (1993) 422. P.W. Anderson, J. Phys. Chem. Solids 11 (1959) 26. B. MuK hlschlegel, D.J. Scalapino, R. Denton, Phys. Rev. B. 6 (1972) 1767. H. Flocard, N. Onishi, Ann. Phys. (NY) 254 (1996) 225. H. Flocard, Ann. Phys. (NY) 191 (1989) 382. R. Balian, M. VeH neH roni, Ann. Phys. (NY) 195 (1989) 324. R. Balian, Y. Alhassid, H. Reinhardt, Phys. Rep. 131 (1986) 1. A. Weinstein, in: Fluids and Plasmas: Geometry and Dynamics, Contemporary Mathematics, vol. 28, Amer. Math. Soc., Providence, RI, 1984.