R Balian
email: roger.balian@cea.fr

M Vénéroni

Correlation functions from a unified variational principle: trial Lie groups

Keywords: numbers: 05.30.-d, 04.25.-g, 24.60.-k, 67.10.-j Many-body systems in equilibrium and non-equilibrium, Variational principle, Generating functional, Correlations, Backward Heisenberg equation, Lie groups Contents I. Introduction

Time-dependent expectation values and correlation functions for many-body quantum systems are evaluated by means of a unified variational principle. It optimizes a generating functional depending on sources associated with the observables of interest. It is built by imposing through Lagrange multipliers constraints that account for the initial state (at equilibrium or off equilibrium) and for the backward Heisenberg evolution of the observables. The trial objects are respectively akin to a density operator and to an operator involving the observables of interest and the sources. We work out here the case where trial spaces constitute Lie groups. This choice reduces the original degrees of freedom to those of the underlying Lie algebra, consisting of simple observables; the resulting objects are labeled by the indices of a basis of this algebra. Explicit results are obtained by expanding in powers of the sources. Zeroth and first orders provide thermodynamic quantities and expectation values in the form of mean-field approximations, with dynamical equations having a classical Lie-Poisson structure. At second order, the variational expression for two-time correlation functions separates-as does its exact counterpart-the approximate dynamics of the observables from the approximate correlations in the initial state. Two building blocks are involved: (i) a commutation matrix which stems from the structure constants of the Lie algebra; and (ii) the second-derivative matrix of a free-energy function. The diagonalization of both matrices, required for practical calculations, is worked out, in a way analogous to the standard RPA. The ensuing structure of the variational formulae is the same as for a system of non-interacting bosons (or of harmonic oscillators) plus, at non-zero temperature, classical gaussian variables. This property is explained by mapping the original Lie algebra onto a simpler Lie algebra. The results, valid for any trial Lie group, fulfill consistency properties and encompass several special cases: linear responses, static and time-dependent fluctuations, zero-and high-temperature limits, static and dynamic stability of small deviations.

I. INTRODUCTION

Variational methods have proved their flexibility and efficiency in many domains of physics, chemistry and applied mathematics, in particular when no small parameter allows perturbative approaches. In physics one may wish, for systems of fermions, bosons or spins, to evaluate variationally various quantities, such as thermodynamic properties, expectation values, fluctuations or correlation functions of some observables of interest. One may deal with the ground state, with equilibrium at finite temperature, or with time-dependences in off-equilibrium situations.

We want, moreover, to perform these evaluations consistently. Suppose for instance that we have optimized the free energy of a system by determining variationally its approximate state within some trial class; nothing tells us that this state is also suited to a consistent evaluation of other properties. Is it possible, remaining in the same trial class, to optimize some other quantity than the free energy, for instance a statistical fluctuation?

The wanted properties can be of different types. For instance, one may be interested in both the expectation values and the correlations of some set of basic observables. Consistency then requires the simultaneous optimization of these quantities. In order to evaluate them in a unified framework it appears natural, as in probability theory and statistical mechanics, to rely on characteristic functions, or more generally for time-dependent problems, on functionals that generate correlation functions. As in field theory, time-dependent sources ξ j (t) are associated with the basic observables Q j . Expansion of the generating functional in powers of these sources will supply expectation values and correlation functions of the set Q j in a consistent fashion.

Our strategy, therefore, will be the variational optimization of a generating functional ψ{ξ} for (connected) correlation functions. To face this problem we will rely on a general method [1,2] allowing the systematic construction of a variational principle that optimizes some wanted quantity. In this procedure, the equations that characterize this quantity are implemented as constraints by means of Lagrange multipliers.

The desired generating functional is expressed as ψ{ξ} ≡ ln TrA(t i ) D in terms of two entities, the state D in the Heisenberg picture and the "generating operator"

A(t) ≡ T exp[i ∞ t dt ′ j ξ j (t ′ ) Q H j (t ′ , t
)] taken at the initial time t = t i . The operators Q H j (t ′ , t) entering A(t) are the observables of interest in the Heisenberg picture, t ′ being the running time and t the reference time at which Q H j (t, t) reduces to the observable Q j in the Schrödinger picture. The variational determination of ψ{ξ} together with its expansion at successive orders in the sources ξ j (t ′ ) provides the various desired outcomes: namely, at zeroth order (for A(t) = I) thermodynamic potentials if D is a (non-normalized) Gibbs state, or ground state energy in the zero-temperature limit; at first order expectation values; at second order correlation functions for an initial off-equilibrium state D, etc... Static correlations within the state D will be obtained from sources located at the origin t i of times.

Implementing the variational principle requires the use of formally simple equations which characterize the two ingredients D and A(t i ) of the generating functional ψ{ξ}, and which will be taken as constraints accounted for by Lagrange multipliers. The state D = exp(-βK) will be characterized by Bloch's equation for D(τ ) = exp(-τ K); the associated Lagrange multiplier is an operator depending on an imaginary time τ . In order to characterize the second ingredient A(t i ) of ψ{ξ}, we have defined A(t) for an arbitrary initial time t. The observables Q H j (t ′ , t) entering the "generating operator" A(t) then satisfy a backward Heisenberg equation [Eq. (10)] in terms of the reference time t which will eventually be fixed at t i . This backward Heisenberg equation plays a crucial role as it produces for A(t) the formally simple differential equation (11); the associated Lagrange multiplier is a time-dependent operator. The equations for the density operator D(τ ) and for the generating operator A(t) are complemented by the boundary conditions D(0) = I and A(+∞) = I, where I is the unit operator. The "time" τ of D(τ ) varies forward, the time t of A(t) backward.

Unrestricted variations of the trial operators D(τ ) and A(t), and of their associated multipliers, yield the exact generating functional, the stationarity conditions being the exact dynamical equations for D(τ ) and A(t). These operators should be restricted within a trial subspace to make the evaluations feasible, and then the resulting equations become coupled. Their solution will be simplified by expansion in powers of the sources ξ j (t), yielding a tractable, unified and consistent treatment.

In this article, we choose as trial subspace a Lie group of operators. This will be the sole approximation. The formalism will be developed for an arbitrary Lie group (Secs. III-V). Explicit calculations are then allowed for a sufficiently simple underlying Lie algebra.

For arbitrary systems and for any trial Lie group, mean-field like approximations are recovered at zeroth and first orders in the sources (Sec. IV) for thermodynamic quantities and static or dynamic expectation values (for instance, selecting for fermions at finite temperature the Lie algebra of single-particle operators yields back in this general framework the static and time-dependent Hartree-Fock theories).

However, at second order in the sources, the formalism generates non trivial results for fluctuations, for static correlations and for two-time correlation functions (Sec. V). Remarkably, new variational approximations come out for these quantities, although the trial operators belong to the same simple class as the one that provides standard results for the expectation values (exponentials of single-particle operators in the example of fermions). Within the Lie group the trial operators adapt themselves to each question being asked so as to optimize the answer -while the use of the generating functional ensures the consistency of the results thus obtained.

The second half of this article (Secs. VI-X) is devoted to the properties of the outcomes of the present variational theory with trial Lie groups. In particular, the static and dynamic stabilities are related to each other (Sec. IX); a quasi-bosonic structure (Sec. VIII) and a classical structure (Sec. X) are exhibited for any system; various consistency properties are reviewed (Sec. VI).

Some past works are related to the present variational approach, which generalizes and unifies them within a natural framework. For fermions at finite temperature, the optimization of expectation values has been shown to lead to the static HF and dynamic TDHF equations [3], while variational expressions for fluctuations [4] and correlation functions [5] have been derived. In particular, the large fluctuations observed in heavy ion nuclear reactions have thus been correctly reproduced [6][7][8][9][10][11], as recalled in Sec. VI B 5. Applications of the variational principle of Sec. II C have been made to bosonic systems [12], to Bose condensation [13,14], to field theory in φ 4 , including two-time correlation functions [15,16], and to restoration of broken particle-number invariance for paired fermions at finite temperature [17]. Let us also mention an application to cosmology [18] and an extension to control theory [19,20]. In the case of small fluctuations around the mean-field trajectory for fermions, the variational principle leads to time-dependent RPA corrections similarly to other approaches assuming either initial sampling of quantum fluctuations [21,22] or directly solving time-dependent RPA [23] (for a recent review, see [START_REF] Lacroix | Progress of time-dependent nuclear reaction theory[END_REF]).

The main results are recapitulated in Sec. XI, before the conclusion in Sec. XII.

II. A UNIFIED VARIATIONAL APPROACH A. Generating functional

We consider a quantum physical system prepared at the initial time t i in a state represented by a density operator of the form D ∝ exp(-βK) in the Hilbert space H . The dynamics is generated by the Hamiltonian H, possibly time-dependent. If the system is in canonical (or grand canonical) equilibrium, one has K = H (or K = H -µN ); for dynamical problems K needs not commute with H. Ground state problems are traited by letting β → ∞.

Partition functions will be evaluated as Tr D for the unnormalized state

D = e -β K . (1)
The normalized state will be denoted as D = D/Tr D. We are mainly interested in expectation values, fluctuations and correlation functions of some set of observables denoted as Q S j in the Schrödinger picture. We will work in the Heisenberg picture in which the observables

Q H j (t f , t i ) = U † (t f , t i ) Q S j U (t f , t i ) (2)
depend on two times, the initial reference time t i and the final running time t f . In the unitary evolution operator

U (t f , t i ) = T e -i t f ti dtH(t) (3)
T denotes time ordering from right to left. In order to generate consistently the desired quantities, we associate with each observable Q j a time-dependent source ξ j (t) and we introduce the generating operator

A(t) ≡ T e i ∞ t dt ′ j ξ j (t ′ ) Q H j (t ′ , t) . (4)
Then, the generating functional

ψ{ξ} ≡ ln TrA(t i ) D , (5)
which depends on the functions ξ j (t), encompasses the quantities of interest. In particular the partition function Tr D, the expectation values

Q j t = Tr Q H j (t, t i ) D (6)
at the time t, and the two-time correlation functions

C jk (t ′ , t ′′ ) = Tr T Q H j (t ′ , t i ) Q H k (t ′′ , t i ) D -Q j t ′ Q k t ′′ , (7)
are obtained as functional derivatives with respect to the sources according to the successive terms of the expansion of ψ{ξ},

ψ{ξ} = ln Tr D + i ∞ ti dt j ξ j (t) Q j t -1 2 ∞ ti dt ′ ∞ ti dt ′′ jk ξ j (t ′ ) ξ k (t ′′ ) C jk (t ′ , t ′′ ) + . . . . (8)
Variances ∆Q 2 j (t), hence fluctuations ∆Q j (t), are found as C jj (t, t); static expectation values and correlations in the state D are found by letting both times equal to t i . Linear responses are also covered by the formalism.

B. The constraints

In order to optimize simultaneously all the quantities embedded in the generating functional ψ{ξ}, use is made of a general procedure [1,2] inspired by the Lagrange multiplier method. The purpose is to construct an expression whose stationary value provides the quantity we are looking for, namely here ψ{ξ} = ln Tr A(t i )D. To implement in this expression the quantities D = e -βK and A(t), we will characterize them by equations regarded as constraints.

To characterize the state D, we introduce a trial "time"-dependent operator D(τ ) compelled to satisfy the initial condition D(0) = I and the Bloch equation

dD(τ ) dτ + K D(τ ) = 0 , (9)
where the imaginary "time" τ runs from 0 to β. One recovers D = exp(-βK) from the solution of ( 9) for τ = β.

Let us turn to the generating operator A(t) defined by ( 4) in terms of the observables Q H j (t ′ , t). In this form, A(t) looks difficult to handle and we wish to characterize it by a formally simple equation that can be taken as a constraint. To this aim, the operators Q H j (t ′ , t) are regarded as functions of the initial time t rather than of the final running time t ′ which is kept fixed at t ′ = t f . They thus satisfy the backward Heisenberg equation [3]

dQ H j (t f , t) dt = - 1 i [Q H j (t f , t), H] . (10)
This differential equation, together with its final boundary condition

Q H j (t f , t f ) = Q S j (t f ) at t = t f , is equivalent to the definition (2) of Q H j (t f , t).
Contrary to the standard forward Heisenberg equation (a differential equation in terms of t ′ ), the backward equation (10) holds even when the observable Q S j and/or the Hamiltonian H depend on time in the Schrödinger picture. Note that in the backward Eq.( 10), H is written in the Schrödinger picture if it is time-dependent.

In the present context, the forward equation for Q H j (t ′ , t) would be of no help in dealing with the definition (4) of the generating operator A(t) whereas the backward Heisenberg equation (10) readily provides the differential equation [5]

dA(t) dt + 1 i [A(t), H] + iA(t) j ξ j (t) Q S j = 0 , (11)
which, together with the boundary condition A(+∞) = I, is equivalent to (4). The generating functional ψ{ξ} defined by (5) involves A(t i ), and this operator will be found by letting t run backward in Eq. (11) from +∞ to t i , with the final condition A(+∞) = I. Here again we shall simulate the operator A(t), solution of the exact equation (11), by a trial operator A(t) which will satisfy (11) approximately.

In order to optimize the generating functional ψ{ξ} = ln TrA(t i ) D, a variational expression [Eq. ( 12) below] will be constructed, which relies on the equations (9) for D and (11) for A(t). These equations are regarded as constraints with which Lagrange multipliers will be associated. We denote the Lagrange multiplier accounting for the Bloch equation ( 9) by A(τ ), an operator depending on the "time" τ ; we denote the Lagrange multiplier accounting for Eq.( 11) by D(t), a time-dependent operator. These notations are inspired by the duality between observables A and states D at the root of the algebraic formulation of quantum mechanics [START_REF] Emch | Algebraic Methods in Statistical Mechanics and Quantum Field Theory[END_REF][START_REF] David | The formalisms of quantum mechanics: an introduction[END_REF] where expectation values are expressed as scalar products Tr AD. Here, D(τ ) (for 0 ≤ τ ≤ β) and the multiplier D(t) (for t ≥ t i ) are state-like quantities, whereas the multiplier A(τ ) and the operator A(t) are observable-like quantities.

C. The variational principle

The implementation of the constraint (9) for D(τ ) by introducing the Lagrange multiplier A(τ ), and of the constraint (11) for A(t) by introducing of the Lagrange multiplier D(t), results in the variational expression [5] Ψ{A, D} = ln Tr A(t

= t i )D(τ = β) (12)
- β 0 dτ Tr A(τ ) dD(τ ) dτ + KD(τ ) [Tr A(τ )D(τ )] -1 + ∞ ti dt Tr   dA(t) dt + 1 i [A(t), H] + i A(t) j ξ j (t) Q S j   D(t)[Tr A(t)D(t)] -1 ,
where normalizing denominators are included for convenience. Together with the mixed boundary conditions

A(t = +∞) = I , D(τ = 0) = I , (13)
Ψ{A, D} should be made stationary with respect to the four time-dependent operators D(τ ), A(τ ), A(t), D(t) (with 0 ≤ τ ≤ β and t i ≤ t ≤ +∞). The stationarity conditions include the additional continuity relations

D(τ = β) = D(t = t i ) and A(τ = β) = A(t = t i ) , (14)
another argument for the notation. (In view of this continuity one might replace τ by a complex time t = t i + i(β -τ ) so as to rewrite the two integrals of ( 12) as a single integral on a Keldysh-like contour [5,27].) For unrestricted variations of A and D the stationary value of Ψ is the required generating functional ψ{ξ}. It is attained for values of A and D that let the two square brackets of ( 12) vanish, so that we recover the evolution equations ( 9) and (11), the solutions of which are D(τ ) = exp (-τ K) and A(t) = A(t).

The data of the problem, K, H and the observables Q S j , are operators in the Hilbert space H . They all enter explicitly the variational expression Ψ{A, D}. Simpler variational principles (VPs) derive from (12) in two special circumstances. If the initial state D is workable, the first integral over τ should be omitted. In this case, the variational principle ( 12) can be viewed for ξ j (t) = 0 as a transposition of the Lippmann-Schwinger VP [28] from the Hilbert space, with duality between bras and kets, to the Liouville space, with duality between observables and states [3]. For static problems, the last integral over t is irrelevant [17]. Classical problems enter the same framework, with the replacement of the Hilbert space by the phase space, traces by integrals and commutators by Poisson brackets.

As usual, practical exploitation of the above variational approach relies on restricting the trial spaces so that the expression (12) of Ψ{A, D} can be explicitly worked out. The denominators Tr AD have been introduced so as to let the functional Ψ{A, D} be invariant under time-dependent changes of normalization of D and A. This allows us to select a "gauge", that is, to fix at each time the traces of D and A in such a way that the stationarity conditions take the form

Tr δA(τ ) dD(τ ) dτ + KD(τ ) = 0 , (0 ≤ τ ≤ β) (15)
Tr

dA(τ ) dτ -A(τ )K δD(τ ) = 0 , (0 ≤ τ ≤ β) (16)
Tr

  dA(t) dt + 1 i [A(t), H] + iA(t) j ξ j (t)Q S j   δD(t) = 0, (t ≥ t i ) (17)
Tr δA(t)

  dD(t) dt - 1 i [H, D(t)] -i j ξ j (t)Q S j D(t)   = 0, (t ≥ t i ) ( 18
)
in the restricted space for D(τ ), A(τ ), A(t), D(t) and the corresponding space for their infinitesimal variations. In agreement with the boundary conditions ( 13) and ( 14) relating the sectors τ and t, equations ( 15) and ( 18) for D should be solved forward in time, with τ running from 0 to β and t running from t i to ∞, whereas Eqs.( 17) and ( 16) for A should be solved backward. We obtain the stationary value of Ψ as

ψ{ξ} ≃ ln Tr A(t) D(t) = ln Tr A(τ ) D(τ ) (19)
for arbitrary t or τ , as a consequence of the stationarity conditions written for δA ∝ A, δD ∝ D. In the following the restricted choice of the trial space will imply that the allowed variations δA around A depend on A, and likewise for δD, so that the forward or backward equations (15)(16)(17)(18) are coupled. Practical solutions will take advantage of their expansion in powers of the sources ξ j (t).

The variational procedure has duplicated the dynamical equations, introducing Eqs.( 16) and ( 18), besides the approximate Bloch equation (15) and the approximate equation (17) for the generating operator A(t). While the formalism was set up in the Heisenberg picture, the stationarity condition (18) reduces, in the absence of sources and for unrestricted variations of A(t) and D(t), to the Liouville-von Neumann equation

dD(t) dt = 1 i [H, D(t)] . (20)
The Lagrange multiplier matrix D(t) thus behaves as a time-dependent density operator in the Schrödinger picture. However, such an interpretation does not hold in the presence of sources, in which case D(t) is not even hermitian.

III. LIE GROUP AS TRIAL SPACE

A. Parametrizations and entropy

From now on, we specialize the trial space for the operators A and D involved in the variational principle (12), assuming it to be endowed with a Lie group structure. This will be our sole approximation. The Lie group is generated by a Lie algebra {M} of operators acting on the Hilbert space H , a basis of which is denoted as M α . This algebra is characterized by the structure constants Γ γ αβ entering the commutation relations

[M α , M β ] = i Γ γ αβ M γ . (21)
(We use throughout the convention of summation over repeated indices.) These constants are antisymmetric and satisfy the Jacobi identity

Γ ǫ αβ Γ δ ǫγ + Γ ǫ βγ Γ δ ǫα + Γ ǫ γα Γ δ ǫβ = 0 .
It is convenient to include in the algebra the unit operator I, denoted as M 0 . A seminal example is provided, for a many-fermion problem, by the algebra of fermionic single-particle operators

M α = a † µ a ν [with α ≡ (ν, µ) = 0]
. Other examples are given in Sec. XI A, such as, for condensed bosons, the set of creation and annihilation operators and of their products in pairs.

The operators A and D are then parametrized, at each time τ or t, according to

A = e L α Mα , D = e J α Mα . ( 22
)
The parameters L α and J α are functions of the times τ or t; their sets will be denoted as {L} and {J}. They will be complex due to the presence of the sources ξ j (t), and to a possible non-hermiticity of the operators M α . [In the fermionic example, the operators D have the nature of non-normalized independent-particle states; for bosonic problems they would encompass coherent states.] The results will be conveniently expressed by writing D as D = Z D, where Z denotes the normalization factor

Z{J} ≡ Tr D = Tr e J α Mα , (23)
and where D is a normalized operator. Instead of the set {J}, the operator D can alternatively be parametrized by Z and by the set {R}, defined by

R α ≡ Tr M α D = ∂ ln Z ∂J α (24)
and including R 0 = 1. [In the example of fermions, for α = (ν, µ) = 0, the set R α ≡ R νµ = Tr a † µ a ν D are the Wick contractions associated with the independent-particle trial operator D, so that a covariant vector {R} can also be regarded as a single-particle density matrix. Contravariant vectors such as {L} or {J} are then regarded as matrices with switched indices, so as to produce the usual expressions in Hilbert space for operators and scalars, such as

L α M α = a † µ L µν a ν or J α R α = J µν R νµ .
] The converse equations of [START_REF] Lacroix | Progress of time-dependent nuclear reaction theory[END_REF], which relate the set {J} to the set {R}, Z, involve the von Neumann entropy

(k B = 1) S{R} ≡ -Tr D ln D ≡ ln Tr D - Tr D ln D Tr D = ln Z -J α R α , (25)
and they read

J α = - ∂S{R} ∂R α (α = 0) , J 0 = ln Z -S{R} - α =0 J α R α . (26)
The Legendre transform ( 25)-( 26) from ln Z{J} to the von Neumann entropy S{R} stems from the exponential form of D as function of the parameters {J}. The situation is the same as in thermodynamics where a Legendre transform relates thermodynamic potentials and entropy when going from intensive to extensive variables, due to the the Boltzmann-Gibbs form of the equilibrium states. The Jacobian matrix of the transformation relating the two parametrizations of D is the entropic matrix

S αβ ≡ ∂ 2 S{R} ∂R α ∂R β = - ∂J α ∂R β , S -1 αβ = - ∂ 2 ln Z{J} ∂J α ∂J β = - ∂R α ∂J β (α, β = 0) , (27)
which is the (negative) matrix of second derivatives of S{R}.

As a remark, we note that a metric ds 2 can be defined in the full space of density operators D by ds 2 = -d 2 S{ D} [29]. The quantity -δR α S αβ δR β can indeed be interpreted as the square of the distance, within the trial subset of density operators, between the state D parametrized by {R} and the state D + δ D parametrized by {R + δR}. The matrix -S can thus be regarded as a metric tensor, and the relation δJ α = -S αβ δR β (α, β = 0) as a correspondence between covariant and contravariant coordinates.

B. Symbols and images

We shall have to deal with quantities of the form Tr QD, where D will be some element of the Lie group and Q some operator in the full Hilbert space H , not necessarily belonging to the Lie algebra. In order to take care of such operators it appears convenient to represent them by means of two useful tools.

Let us first introduce the symbol q{R} of Q, a scalar which depends both on the operator Q (in the Schrödinger picture) and on a normalized running element D of the Lie group. This symbol is defined by

q{R} ≡ Tr Q D , (28)
as function of the parameters R α (with α = 0) that characterize D. If Q belongs to the Lie algebra, this function is linear since R α = Tr M α D itself is the symbol of M α . Otherwise q{R} is non-linear. If D were a density operator, a symbol would be an expectation value but here D is an arbitrary normalized element of the Lie group, not necessarily hermitian. The symbol q{R} of the operator Q may be viewed as a generalization, for any Lie group and for mixed states, of the expectation value of Q in a coherent state [START_REF] Perelomov | Generalized Coherent States and Their Applications[END_REF], this value being regarded as a function of the parameters that characterize this state.

Let us then introduce a second object associated with an operator Q, its image Q{R}, an element of the Lie algebra depending again both on Q and on the running operator D. The image Q{R} of Q is constructed by requiring that both operators Q{R} and Q should be equivalent, in the sense that

Tr Q{R} D = Tr Q D = q{R} , (29)
and that

Tr Q{R} δD = Tr Q δD [START_REF] Perelomov | Generalized Coherent States and Their Applications[END_REF] for any infinitesimal variation δD around D within the Lie group, in particular for δD ∝ M α D, or δD ∝ D M α , or δD ∝ ∂ D/∂R α . Let us show that these conditions are sufficient to determine uniquely the image Q{R} associated with a given Q. Since it must belong to the Lie algebra, Q{R} is parametrized by a set of coordinates Q α {R} according to

Q{R} = Q α {R} M α . (31)
For α = 0, these coordinates Q α {R} are determined by inserting (31) into [START_REF] Perelomov | Generalized Coherent States and Their Applications[END_REF], by taking δD ∝ (∂ D/∂R β )δR β and by using the definition (28) of the symbol of Q, which yields

Q α {R} = ∂q{R} ∂R α , (α = 0) . (32)
The coordinate Q 0 {R} is obtained from (29). Altogether, the image Q{R} of Q defined by their equivalence (29), [START_REF] Perelomov | Generalized Coherent States and Their Applications[END_REF] is related to the symbol q{R} of Q through

Q{R} = q{R} M 0 + (M α -R α M 0 ) ∂q{R} ∂R α (α = 0) . (33)
An operator Q belonging to the Lie algebra coincides with its image. Otherwise, its coordinates Q β depend on the R α 's, so that Q{R} is an effective operator simulating Q in the Lie algebra, for a state close to D. The dependence of Q α {R} on {R} arises from the occurrence of D in the equivalence relation (29), [START_REF] Perelomov | Generalized Coherent States and Their Applications[END_REF].

C. The commutation matrix C and the entropic matrix S

We shall have to handle products of two operators of the Lie algebra. When solving the dynamical equations we shall in particular encounter commutators [M α , M β ]. Their symbol defines the commutation matrix

C αβ {R} = 1 i Tr [M α , M β ] D = Γ γ αβ R γ , (34)
expressed in terms of the structure constants Γ γ αβ (for γ = 0, we have R 0 = 1). The matrix C will play a crucial role. Using ln D = J γ M γ , the vanishing of Tr[M β , ln D] D = 0 is expressed by

C βγ J γ = Γ δ βγ R δ J γ = 0 . ( 35
)
Taking the derivatives of this identity with respect to R α and using the relation (27) between the two parametrizations {J} and {R} of D, one obtains for the product C S the relations

Γ α βγ J γ = C βγ S γα , (α = 0) (36)
Γ 0 βγ J γ = -C βγ S γδ R δ .
These equations help to show that the automorphism of the Lie algebra engendered by the element D λ of the group can be expressed in the two equivalent forms

M α -R α → D -λ (M α -R α )D λ = e i C S λ γ α (M γ -R γ ) . ( 37
)
(This is proved by evaluating the derivative with respect to λ of the left-hand-side, using [START_REF] Cloizeaux | Many-Body Physics[END_REF] then integrating from 0 to λ.)

The property (37) will be exploited later on (Sec. V E and Appendix A). We use it here to find the expression of the symbol of the product M α M β , or equivalently of the correlation Tr

M α M β D -R α R β = Tr(M α -R α )M β D (with R α = Tr M α D).
To this aim, we start from the expression (27) of S -1 , and evaluate explicitly therein the derivatives with respect to {J}:

-S -1 αβ = ∂ 2 ln Z{J} ∂J α ∂J β = ∂ ∂J α Tr e J γ Mγ M β Tr e J γ Mγ = Tr 1 0 dλ D1-λ (M α -R α ) Dλ M β , (38)
where we made use of the first-order expansion in the shift {δJ} of

exp (J γ M γ + δJ γ M γ ) ≈ D + δJ γ 1 0 dλ D 1-λ M γ D λ , (39)
with D = exp J γ M γ . We recognize in the r.h.s. of [START_REF] Mackay | Hamiltonian Dynamical Systems[END_REF] the Kubo correlation of M α and M β in the normalized state D. By means of (37) the integration over λ can be performed explicitly in [START_REF] Mackay | Hamiltonian Dynamical Systems[END_REF], which yields

-(S -1 ) αβ = Tr 1 0 dλ e i C S λ γ α (M γ -R γ ) M β D (40) = e i C S -I i C S γ α Tr (M γ -R γ ) M β D .
Hence the ordinary correlations between operators of the Lie algebra in any element D of the Lie group are found to be given by

Tr M α M β D -R α R β = - i C S e i C S -I S -1 αβ . (41)
D. The variational formalism in the restricted space

With the above tools in hand, it is possible to implement the Lie-group form (22) of the trial operators into the variational expression (12) of Ψ{A, D} by taking as variables the parameters {L}, {R}, Z that characterize A and D at the times τ or t. These operators A and D appear within Ψ{A, D} through products AD and DA. Such products belong to the Lie group and are characterized by the parameters R AD The initial state operator K = -β -1 ln D, the Hamiltonian H and the observables Q S j enter Ψ{A, D} through traces of the form Tr QA D/Tr A D and Tr QD A/Tr D A, where Q stands for K, H or Q S j . We are thus led to introduce, for any operator D parametrized by {R}, the symbols

k{R} ≡ Tr K D , h{R} ≡ Tr H D , q j {R} ≡ Tr Q S j D (42)
of K, H and Q S j . These symbols occur within Ψ{A, D} for values of {R} equal to {R AD } or {R DA }. The variational expression Ψ{A, D}, when specialized to a Lie group, then takes the form

Ψ{A, D} = ln Tr A(t = t i ) D(τ = β) (43)
-

β 0 dτ ∂ ln Z AD ∂R α dR α dτ + d ln Z dτ + i k{R DA } + ∞ ti dt   ∂ ln Z AD ∂L α dL α dt + 1 i h{R DA } -h{R AD } + i j ξ j (t) q j {R DA }   ,
which should be regarded as a functional of the original trial parameters {L}, {R} taken at times τ and t, and Z(τ ) = Tr D(τ ).

The stationarity conditions (15)(16)(17)(18), obtained by functional derivation with respect to these parameters, now read

dD(τ ) dτ + K{R DA } D(τ ) = 0 , (0 ≤ τ ≤ β) (44)
dA(τ ) dτ -A(τ ) K{R DA } = 0 , (0 ≤ τ ≤ β) (45)
dA(t) dt + 1 i A(t) H{R DA } -H{R AD } A(t) + i j ξ j (t) A(t) Q j {R DA } = 0 , (t ≥ t i ) ( 46
)
dD(t) dt + 1 i D(t) H{R AD } -H{R DA } D(t) -i j ξ j (t) Q j {R DA } D(t) = 0 , (t ≥ t i ) . ( 47
)
These equations involve the images K{R}, H{R} and Q j {R} issued, according to the general relation [START_REF] Ring | The Nuclear Many-Body Problem[END_REF], from the derivatives of the corresponding symbols k{R}, h{R} and q j {R} defined by (42). The stationarity conditions ( 44)-( 47) should be solved with the boundary conditions ( 13) and ( 14). When they are satisfied, Z AD is constant in τ and t, and the optimal estimate for the generating functional ψ{ξ} reduces to the first term ln Z AD of ( 43), in agreement with (19). For most Lie groups, solving the coupled equations (44)(45)(46)(47) is hindered by the need of expressing explicitly the parameters {R DA }, {R AD }, Z AD of AD and DA in terms of those ({L}, {R} and Z) of A and D. However, we are interested in the first terms of the expansion of ψ{ξ} in powers of the sources ξ j (t). Accordingly we shall only need to express, as functions of the parameters {R} and Z of an arbitrary element D of the Lie group, the following ingredients: (i) the symbols (42) of the operators K, H and Q S j , and (ii) the entropy function [START_REF] Emch | Algebraic Methods in Statistical Mechanics and Quantum Field Theory[END_REF] which allows us to relate the sets {R} and {J}. This is feasible for many Lie groups. [In the example of the algebra of single-fermion operators a † µ a ν , this is achieved by Wick's theorem.] Thus, explicit solutions of the equations of motion will be found at the first few orders in the sources.

IV. ZEROTH AND FIRST ORDERS

A. Thermodynamic quantities

At zeroth-order in the sources {ξ}, the quantity of interest is the partition function Tr e -βK , or the "generalized free energy"

F ≡ -β -1 ln Tr e -β K = -β -1 ψ{ξ = 0} (48) (k B = 1, T = β -1
) which reduces to the standard free energy for K = H, or to the grand potential for

K = H -µN . It is variationally approximated by -β -1 Ψ{A (0) , D (0) } = -β -1 ln Tr A (0) (t = t i )D (0) (τ = β)
, where the upper index denotes the order in the sources {ξ}.

For {ξ} = 0, the stationarity condition (46) yields A (0) (t) = I for t ≥ t i , hence, from Eq.( 14),

A (0) (τ = β) = I and F ≃ -β -1 ln Tr D (0) (τ = β). Thus D (0) (τ = β) ≡ D (0)
appears as an approximation, variationally suited to the evaluation of thermodynamic quantities, of the exact state D = e -β K .

To obtain D (0) (τ = β) we have to solve the first two stationarity conditions (44)(45) with

D (0) (0) = A (0) (β) = I. We make the Ansatz D(τ )A(τ ) = D (0) , (49)
where

D (0) ≡ D (0) (τ = β
) is a constant operator, still to be determined and characterized by its parameters R

α ≡ Tr M α D(0) for α = 0 and Z (0) ≡ Tr D (0) . The image K{R DA } is then the constant operator K{R (0) }, so that we can solve Eqs.( 44) and ( 45) in the form

D (0) (τ ) = [D (0) ] τ /β , A (0) (τ ) = [D (0) ] (β-τ )/β , where D (0) is determined by the equation ln D (0) = -β K{R (0) } . ( 50
)
The operator equation (50) provides D (0) (τ = β) = D (0) . More explicitly, in the basis {M} of the Lie algebra, the components α = 0 of (50) read J (0)α = -β K α {R (0) } in terms of the coordinates J (0)α of ln D (0) ≡ J (0)α M α and the coordinates K α {R (0) } of K{R (0) } [defined as in [START_REF] Ornstein | Proc. Acad. Sci. Amsterdam[END_REF][START_REF] Broomfield | [END_REF][START_REF] Ring | The Nuclear Many-Body Problem[END_REF]]. This yields the self-consistent equations

∂S{R (0) } ∂R (0) α = β ∂k{R (0) } ∂R (0) α (51)
which determine the parameters {R (0) } of D(0) . The component α = 0 of (50) yields

J (0)0 = k{R (0) } - R (0) α ∂k{R (0) }/∂R (0)
α , where Eq.( 33) has been used for K. Together with ( 26) and (51), this equation provides for the sought free energy F the alternative expressions

F ≃ -β -1 ln Tr e -β K{R (0) } = -β -1 ln Z (0) = f {R (0) } , (52)
where we defined the free-energy function f {R} through

f {R} ≡ k{R} -T S{R} . ( 53
)
The above self-consistent equations determine a hermitian operator D(0) which can be interpreted as an approximation for the exact density operator D (whereas the trial operator D occurring in the presence of sources is not hermitian). The general relation [START_REF] Thouless | [END_REF] entails

Γ γ αβ K β {R (0) } R (0) γ = 0 , (54)
an equation which is equivalent to Tr M γ [ln D (0) , D(0) ] = 0. [In the Hartree-Fock example, (54) expresses the commutation of the single-particle density matrix with the effective hamiltonian matrix.]

Our variational principle provides solutions for ψ{ξ = 0} = -β F that are not maxima but only stationary values of Ψ{A, D}; it relies on the Bloch equation rather than on the maximization of the entropy under constraints. However, it turns out that the above result (51) coincides with the outcome of the standard maximum entropy procedure. Indeed the latter amounts to minimizing the left-hand side of the Bogoliubov inequality

Tr K D + β -1 Tr D ln D ≥ -β -1 ln Tr e -β K ≡ F (55)
with respect to the normalized trial density operator D. When D in (55) is restricted to the Lie group, the left-hand side reduces to f {R}, where {R} is the set characterizing D. The best estimate for F is thereby the minimum of the free-energy function f {R}, which requires that the equations (51) are satisfied. The equivalence between the two variational approaches provides a criterion for selecting the best solution {R (0) } of the self-consistent stationarity conditions (51) when they have several solutions, namely, the one for which f {R (0) } is the absolute minimum of f {R}.

The standard relations S = -∂F/∂T ≃ -Tr D(0) ln D(0) = S{R (0) } and K = F + T S ≃ k{R (0) } are satisfied as usual, so that the approximation is thermodynamically consistent. Thermodynamic coefficients are obtained by derivation, which introduces the matrix (α, β = 0)

F αβ = ∂ 2 f {R (0) } ∂R (0) α ∂R (0) β = K αβ -T S αβ = ∂ 2 k{R (0) } ∂R (0) α ∂R (0) β -T ∂ 2 S{R (0) } ∂R (0) α ∂R (0) β . ( 56
)
In particular, the heat capacity is found for K = H as

C ≃ β K α {R (0) }(F -1 ) αβ K β {R (0) } . ( 57
)
The positivity of the stability matrix F at the minimum of f {R}, entailed by the inequality (55), will play an important role below. [For the fermionic single-particle Lie algebra a † µ a ν , we recover the thermal HF approximation, either under the form (50) or through the minimization of f {R}.]

B. Expectation values

Expectation values are obtained by expanding the generating functional ψ{ξ} at first order in the sources ξ j (t) [Eq. (8)]. These sources occur both directly in Ψ{A, D}, as exhibited by the last term of (43), and indirectly through the values of A and D at the stationarity point where ψ{ξ} = Ψ{A, D}. At this point, however, the partial derivatives of Ψ{A, D} with respect to A and D vanish; one is left with the explicit derivative ∂Ψ{A (0) , D (0) }/∂ξ j (t) taken at the zeroth-order point {A (0) ,D (0) }, so that

Q j t = 1 i ∂ψ{ξ} ∂ξ j (t) ξ=0 ≃ Tr Q S j D (0) (t) A (0) (t) Tr A (0) (t) D (0) (t) (58)
involves only the zeroth-order approximations A (0) (t) and D (0) (t).

We have seen that A (0) (t) = I for t ≥ t i , and that D(0) (t i ), equal to D(0) (τ = β) ≡ D(0) , is given self-consistently by (50). Hence, the static expectation value Q j ti in the state D is variationally expressed by

Tr Q S j D = Tr Q S j e -βK Tr e -βK ≃ Tr Q S j D(0) ≡ q j {R (0) } = Q α j {R (0) } R (0) α , (59)
that is, by the symbol (28) of the observable Q S j of interest. For dynamical problems, the expectation value (58) involves D(0) (t), provided by the last stationarity condition (47) taken for ξ j (t) = 0 and for A (0) (t) = I. Thus, D(0) (t) is determined by the zeroth order self-consistent equation

d D(0) (t) dt = 1 i H{R (0) (t)}, D(0) (t) , (60)
with the initial condition D(0) (t i ) = D(0) . [The norm Z (0) (t) is constant and equal to Z (0) .] In (60) the operator

H{R (0) (t)} is the image taken for R (0) α (t) = Tr M α D(0) (t)
of the Hamiltonian H. The time-dependent expectation values (58) are therefore found as

Q j t ≡ Tr Q H j (t, t i ) D ≃ Tr Q S j D(0) (t) ≡ q j {R (0) (t)} . ( 61
)
As noted at the end of Sec. III, the variational equation ( 60) for the Lagrange multiplier D (0) (t) comes out as an approximation for the Liouville-von Neumann equation (20).

In coordinate form, the equations of motion for the variables R (0)

α (t) = M α t that parametrize D(0) (t) are found from (60) as dR (0) α (t) dt = C αβ {R (0) (t)} ∂h{R (0) (t)} ∂R (0) β (t) , (62)
with the initial conditions {R (0) (t i )} = {R (0) }. [For the single-fermion Lie algebra, we recover the time-dependent Hartree-Fock (TDHF) approximation for the single-particle density matrix R (0) νµ (t), with the static HF solution as initial condition. Indeed, the multiplication in (62) by C αβ {R} produces for this algebra a commutation with R (0) νµ (t). The usual single-particle effective HF Hamiltonian comes out as the image H{R (0) (t)}. The current use of TDHF to evaluate expectation values is thus given a variational status.] Equations ( 62) can be rewritten in the alternative form

dR (0) α dt = L γ α R (0) γ (t) , (63)
where the matrix

L γ α {R (0) (t)} = Γ γ αβ H β {R (0) (t)} (64)
plays the role of an effective Liouvillian acting in the Lie algebra; the quantities H β {R} are the coordinates H β {R} ≡ ∂h{R}/∂R β [Eq.( 32)] of the image H{R} of H.

The dynamics of the density operator D(0) (t) takes therefore the classical form (63) in terms of the scalar variables R (0) α parametrizing D(0) (t). This classical structure will be analyzed in Sec. X. Non-linearity occurs through the restriction of the trial space: the image H{R} depends on the variables {R (0) (t)} when H does not belong to the Lie algebra.

V. CORRELATION FUNCTIONS

Variational approximations for the two-time correlation functions C jk (t ′ , t ′′ ) are provided by the second-order term of the expansion (8) of the generating functional ψ{ξ} in powers of the sources ξ j (t). However, we have seen that, owing to the stationarity of Ψ, zeroth order was sufficient to determine the variational approximation for expectation values. Likewise, it is sufficient here to expand up to first order the first derivative of the generating functional ψ{ξ} according to

1 i ∂ψ{ξ} ∂ξ k (t ′′ ) ≈ Q k t ′′ + i ∞ ti dt ′ j ξ j (t ′ ) C jk (t ′ , t ′′ ) (65) ≃ 1 i ∂Ψ{A, D} ∂ξ k (t ′′ ) = Tr Q S k D(t ′′ ) A(t ′′ ) Tr D(t ′′ ) A(t ′′ ) ,
where the stationarity of Ψ{A, D} with respect of A and D has been used. We can thus obtain C jk (t ′ , t ′′ ) from the first-order contribution to the r.h.s. of (65), and only A (0) , D (0) , A (1) and D (1) have to be determined from the variational equations (44)(45)(46)(47).

A. The approximate backward Heisenberg equation

One building block for correlation functions will be the quantity Q H k (t ′′ , t) defined for t ′′ > t by expanding the trial operator A(t), which enters (65), up to first order as

A(t) ≈ A (0) (t) + A (1) (t) ≡ I + i ∞ t dt ′′ k ξ k (t ′′ ) Q H k (t ′′ , t) . ( 66
)
Comparison with the expansion of the exact generating operator A(t) defined by (4) shows that Q H k (t ′′ , t) simulates the Heisenberg observable Q H k (t ′′ , t). Since A(t) belongs to the Lie group, Q H k (t ′′ , t) belongs to the Lie algebra and can be expressed as Q Hα k (t ′′ , t)M α . The stationarity condition (46) with respect to A(t), expanded up to first order, determines the coordinates Q Hα k (t ′′ , t). For α = 0 these obey for t ≤ t ′′ the equations

dQ Hα k (t ′′ , t) dt = -Q Hβ k (t ′′ , t) (L + C H) α β (α, β = 0) , (67)
which appear as the reduction in the Lie algebra of the backward Heisenberg equation (10) for Q H k (t ′′ , t). The matrix L is the effective Liouvillian (64), C is the commutation matrix defined by [START_REF] Blaizot | Quantum Theory of Finite Systems[END_REF] and H is the second-derivative matrix of the symbol h{R} of the Hamiltonian H,

H αβ {R} = ∂ 2 h{R} ∂R α ∂R β . ( 68
)
These three matrices are functions of {R}; in (67) they are taken at the point

R α = R (0)
α (t) determined by the zeroth-order equations (63). Equations (67) should be solved backward in time from the final boundary condition

Q Hα k (t ′′ , t ′′ ) = Q α k {R (0) (t ′′ )} = ∂q k {R (0) (t ′′ )} ∂R (0) α (t ′′ ) , (69)
where

Q α k {R}M α = Q k {R} is the image of the observable Q k while q k {R}
is its symbol. The exact Ehrenfest equation for Tr M α D(t) [issued from the Liouville-von Neumann Eq.( 20)] and the exact Heisenberg equation for Q H k (t ′′ , t ′′ ) involve the same Liouvillian. However, in the variational treatment, the corresponding approximate Ehrenfest equation ( 63) and the approximate Heisenberg equation (67) differ; the latter contains the corrective term CH in addition to the effective Liouvillian L. This difference arises because the variational equations for the set R (0) α (t) occur at zeroth order in the sources whereas those for the set Q Hα k (t ′′ , t) occur at first order. From the stationarity condition (46) for A(t) one also finds the component Q H0 k (t ′′ , t) of A (1) . This gives an alternative expression for the time-dependent expectation value of Q k at a time t,

Q k t ≃ q k {R (0) (t)} = Q Hα k (t, t ′′ ) R (0) α (t ′′ ) = Tr Q H k (t, t ′′ ) D(0) (t ′′ ) (70)
(including α = 0 with R (0) 0 = 1), which holds for any intermediate time t ′′ , and thus interpolates the Schrödinger picture for t ′′ = t [Eq.( 61)] and the Heisenberg picture for t ′′ = t i , as does the exact expression

Q k t = Tr Q H k (t, t ′′ )[U (t ′′ , t i ) D U † (t ′′ , t i )] . (71)
B. Bypassing D (1) (t) for real times t

The correlation functions that we want to determine through (65) depend on the combination D (1)

(t ′′ ) + D (0) (t ′′ ) A (1) (t ′′ ), since A (0) = I.
The time-dependence of A (1) (t ′′ ) has been expressed through Eqs.(66) and (67). We still need the other ingredient, the first-order operator D (1) (t ′′ ) whose evolution is not simple. However, using the coupled equations of motion for D(t) and A(t), we will show below that the time t ′′ occurring in (65) can be shifted down to t i .

From Eqs.( 46) and ( 47), one can derive for the products D(t)A(t) and A(t)D(t) the uncoupled equations

d(DA) dt = 1 i H{R DA }, DA + i j ξ j (t) Q j {R DA }, DA , (72)
d(AD) dt = 1 i H{R AD }, AD . (73)
These equations cannot be fully solved because the boundary conditions on A and D occur at different times [Eq.( 13)]. Nevertheless, Eq.(72) will happen to be sufficient for our purpose.

In parallel with Eq.(66), we parametrize at first order the product D(t)A(t) by R

(1)

jβ (t, t ′ ) according to R DA β (t) = Tr M β D(t) A(t) Tr D(t) A(t) ≈ Tr M β [ D(0) (t) + δ D(t)] (74) = R (0) β (t) + i ∞ ti dt ′ j ξ j (t ′ ) R (1) jβ (t, t ′ ).
We have in (74) denoted as δ D(t) the first-order variation of DA/Tr DA around D(0) (t). The quantity C jk (t ′ , t ′′ ) that we wish to evaluate now satisfies from (65) and (74) the relation

+ i ∞ ti dt ′ j ξ j (t ′ ) C jk (t ′ , t ′′ ) = Tr Q S k δ D(t ′′ ) . (75)
Since δ D(t ′′ ) is a variation around D(0) (t ′′ ) within the Lie group, the image property [START_REF] Perelomov | Generalized Coherent States and Their Applications[END_REF] allows us to replace the operator 74) and (75), we can get rid of the sources ξ j (t ′ ), so that

Q S k by its image Q k {R (0) (t ′′ )} = Q β k {R (0) (t ′′ )}M β . Comparing Eqs.(
C jk (t ′ , t ′′ ) = R (1)
jβ (t ′′ , t ′ ) Q β k {R (0) (t ′′ )} . ( 76
)
Remember that the coordinates

Q β k {R (0) (t ′′ )} of the image of Q S k are related to the symbol q k {R (0) (t ′′ )} ≡ Tr Q S k D(0) (t ′′ ) by Q β k {R (0) (t ′′ )} = ∂q k {R (0) (t ′′ )} ∂R (0) β (t ′′ ) . (77)
The time-dependence of R

jβ (t, t ′ ) is found by expanding Eq.( 72) at first order. Noting that TrDA does not depend on time, one obtains

dR (1) jβ (t, t ′ ) dt = (L + C H) γ β R (1) jγ (t, t ′ ) + i C βγ Q γ j {R (0) (t)} δ(t -t ′ ), (78)
where again L, C and H depend on time through {R (0) (t)}. The same kernel L + C H as in the approximate backward Heisenberg equation ( 67) is recovered; it will also be encountered in the context of small deviations (Sec. IX B). [For fermionic systems, L + C H is the time-dependent RPA matrix issued from the kernel L of the TDHF equation ( 63).] As a consequence of the duality between the kernels of Eqs.( 67) and (78), we obtain the identity

d dt [R (1)
jβ (t, t ′ ) Q Hβ k (t ′′ , t)] = -i Q γ j {R (0) (t)} C γβ Q Hβ k (t ′′ , t) δ(t -t ′ ) . (79)
Since Q Hβ k (t ′′ , t) vanishes for t ′′ < t, the r.h.s. of (79) disappears if t ′′ < t ′ . One can therefore evaluate C jk (t ′ , t ′′ ) for t ′ ≥ t ′′ from (76) by noting that the product R

jβ (t, t ′ ) Q Hβ k (t ′′ , t) does not depend on t for t ′ ≥ t ′′ > t > t i . Using the boundary condition Q Hβ k (t ′′ , t ′′ ) = Q β k {R (0) (t ′′
)}, we then shift t ′′ down to t i in (76), which yields

C jk (t ′ , t ′′ ) = R (1)
jβ (t i , t ′ ) Q Hβ k (t ′′ , t i ) , (t ′ ≥ t ′′ ) . ( 80
)
The continuity condition D(t i )A(t i ) = D(β)A(β) entails that this correlation function (80) is equal to

C jk (t ′ , t ′′ ) = R (1) jβ (β, t ′ ) Q Hβ k (t ′′ , t i ) (t ′ ≥ t ′′ ) (81)
in terms of the boundary value at τ = β of {R

(1) j (τ, t ′ )}. One needs therefore to determine, only in the interval 0 ≤ τ ≤ β, the set {R (1) j (τ, t ′ )} defined through the first-order parametrization of D(τ )A(τ ) by

R DA β (τ ) = Tr M β D(τ ) A(τ ) Tr D(τ ) A(τ ) ≈ R (0) β + i ∞ ti dt ′ j ξ j (t ′ )R (1)
jβ (τ, t ′ ) , (82)
where

R (0) β = Tr M β D(0) .
The evaluation of D(t) at times t > t i has been bypassed, and it remains to solve the coupled equations ( 44) and (45) for A (1) (τ ) and D (1) (τ ) in the range 0 ≤ τ ≤ β, with the boundary conditions D (1) (τ = 0) = 0 and

A (1) (τ = β) = i ∞ ti dt ′ j ξ j (t ′ ) Q Hα j (t ′ , t i ) M α (83)
issued from ( 13) and (66), respectively. For t ′′ ≥ t ′ , the symmetry C jk (t ′ , t ′′ ) = C kj (t ′′ , t ′ ) can be checked by means of the r.h.s. of Eq.(79).

C. Two-time correlation functions

The linear structure of the boundary condition (83) for A (1) (τ ), together with the boundary condition D (1) (τ = 0) = 0 for D (1) (τ ), entail the occurrence of an overall factor Q Hα j (t ′ , t i ) in A (1) (τ ) and D (1) (τ ), hence in R

jβ (τ, t ′ ), and finally in the correlation function C jk (t ′ , t ′′ ). We had already acknowledged the explicit dependence of C jk (t ′ , t ′′ ) on Q Hβ k (t ′′ , t i ). Therefore, for any trial Lie group, the two-time correlation functions have the general form (α, β = 0)

C jk (t ′ , t ′′ ) ≃ Q Hα j (t ′ , t i ) B αβ Q Hβ k (t ′′ , t i ) , (t ′ > t ′′ ) . ( 84
)
Beside the approximate Heisenberg observables given by the dynamical equations (67) and the boundary condition (69), a matrix B αβ appears, which depends only on the zeroth and first order contributions to D(β) and A(β).

Noticeably, the expression (84) has the same factorized structure as the exact formula (7). In spite of the coupling between the stationarity conditions in the sectors t ≥ t i and 0 ≤ τ ≤ β, Eq.( 84) [as Eq.( 7)] displays separately two types of ingredients: The result B of the optimization of the initial state is disentangled from that of the dynamics, embedded in the approximate Lie-algebra Heisenberg observables

Q H j (t ′ , t i ) and Q H k (t ′′ , t i ).

D. The correlation matrix B

The last task consists in determining explicitly the correlation matrix B αβ . We have noted above that, through the boundary condition (83), the factor Q Hα j (t ′ , t i ) occurs in R

jβ (τ, t ′ ). By introducing new (α-indexed) operators A

jβ (τ, t ′ ) as R (1) jβ (τ, t ′ ) = Q Hα j (t ′ , t i ) R (1) αβ (τ ) . (85)
The quantities R

αβ (τ ) defined by ( 82) and (85) will then be given by R

αβ (τ ) = Tr [D (0) (τ ) A (1) α (τ ) + D (1) α (τ ) A (0) (τ )](M β -R (0) β ) / Tr D (0) , (1)
and the correlation matrix B αβ will be found as

B αβ = R (1)
αβ (τ = β), that is,

B αβ = Tr M α (M β -R (0) β ) D(0) + Tr D (1)
α (β)(M β -R (0) β ) Tr D (0) . ( 87
)
It remains to determine the coupled operators D

α (τ ) and A

α (τ ). The solution is explicitly worked out in Appendix A.1. It had already been given in special cases, for fermionic systems [5], extended BCS theory [17], φ 4 quantum field [15].

The result thus found for the correlation matrix B is

B = i C F 1 1 -exp(-i β C F) F -1 (88) = 1 2 i C F coth 1 2 i β C F F -1 + 1 2 i C .
It involves the commutation matrix C defined by [START_REF] Blaizot | Quantum Theory of Finite Systems[END_REF] and the positive matrix F = K -T S of second derivatives (56) of the trial free energy. Both are taken for the parameters {R (0) } that were determined at zeroth order. For vanishing eigenvalues of iCF, the coefficient of F -1 is meant as β -1 . [For the single-particle fermionic Lie group, the matrix iCF is the static RPA kernel.] By letting t ′ -0 = t ′′ = t i and Q S j = M α , Q S k = M β in (84), one identifies B αβ as the variational approximation for the correlations, in the initial state D, of the operators M α that span the Lie algebra:

Tr M α M β D -Tr M α D Tr M β D ≃ B αβ , (α, β = 0) . ( 89
)
The variational expressions (84), (88) [together with Eqs.(67)] for the correlation functions issued from the use of a Lie-group trial space are the most important outcomes of the present variational approach. We comment below the features of this expression, and work out its consequences in the forthcoming sections.

E. Status of the result for static correlations; Kubo correlations

The optimization of thermodynamic quantities and expectation values (Sec. IV) has resulted in a mean-field type of approximation, with the mere replacement of the exact state D by the zeroth-order contribution D (0) to the trial object D, as in Tr

M α D ≃ Tr M α D(0) = R (0)
α . However, the optimized approximation (88),(89) that we found for the correlations B αβ of M α and M β does not follow from such a simple replacement. Evaluated in the state D(0) by the formula [START_REF] Weinstein | Fluids and Plasmas : Geometry and Dynamics[END_REF], such correlations, instead of (88), would be given by

Tr M α M β D(0) -R (0) α R (0) β = - i C S exp[i C S] -1 1 S -1 αβ , (90)
where S and C are evaluated from ( 27) and ( 34) for {R} = {R (0) }. Contrary to B, the naive expression (90) is not variationally optimized. While the first term of the variational expression (87) of B provides a contribution equal to (90), its second term, arising from D

α (β), introduces a correction which substitutes S -βK = -βF to S . The matrix K takes into account effects coming from the part of the operator K that lies outside the Lie algebra. [When {M} is the fermionic single-particle algebra, the left-hand-side of (90) is the Fock term from Tr (a

† µ a ν )(a † σ a τ ) D(0) since R (0) α R (0)
β is the Hartree term. The full matrix B involves an RPA kernel, in which the matrix K is the effective two-body interaction.] Although the approximations for expectation values and correlations functions stem from the same variational principle, they appear intrinsically different. Not only B cannot be expressed from D (0) , but moreover there is no density operator approximating D in the original Hilbert space H that would produce the optimized correlations B in the same form as (89). While D(0) can be interpreted as a state, the trial operator D has no perturbative status and is just a calculational tool involving the sources. In the expression (87) of B, the operator D (1) α , which depends on M α , is not a correction to D(0) .

Nevertheless, we will show in Sec. VII that, through a mapping of the original Hilbert space H into a new space H and of the Lie algebra {M} into a reduced algebra {M}, the matrix elements of B can be interpreted as exact correlations between the operators {M} in an effective state D. The quantities {R (0) } will also appear as exact expectation values of {M} in D. The variational results, for both expectation values and correlations, will thus be unified through a modification of the Lie algebra.

One may also wonder about the origin of the Bose-like factor exhibited, for arbitrary Lie groups, by the expression (88) of B. A clue will be given in Secs. VII and VIII where an algebra of Bose operators arises from the mapped Lie algebra {M}.

An alternative understanding of the structure of the matrix B can be reached by relating it to the matrix of Kubo correlations, even though the latter have less direct physical relevance than ordinary correlations. Kubo correlations between the operators M α and M β are defined for the exact state D by

Tr 1 β β 0 dτ e τ K M α e -τ K M β D -Tr M α D Tr M β D . (91)
We have seen in Sec. III C that, if the state D were bluntly replaced in (91) by the element D(0) of the Lie group, Kubo correlations would be directly related by [START_REF] Mackay | Hamiltonian Dynamical Systems[END_REF] to the matrix S. We show in Appendix A.2 that a variational approximation B K for the Kubo correlations (91) reads

B K = 1 β (F -1 ) = (β K -S) -1 . ( 92
)
The naive approximation [START_REF] Mackay | Hamiltonian Dynamical Systems[END_REF], namely B K ≃ -S -1 , obtained by replacing D by D(0) in (91), is thus variationally corrected by inclusion in (92) of the term β K. Likewise, the naive approximation (90) for the ordinary correlations in the state D results from the variational expression (88) for B by omission of K within F = K -T S.

Once B K has been obtained in the form (92), it is possible to recover from it the expression (88) for B. In the special case K = 0, for which β F = -S, this was achieved in Sec. III C. We saw there that the factor [D (0) ] -τ /β (M α -R α )[D (0) ] τ /β entering the Kubo correlation (91) is expressed as e i C S τ /β γ α (M γ -R γ ) through the automorphism (37) of the Lie algebra. Integration of e i C S τ /β between 0 and β yields

B K = e i C S -I i C S B , (93)
in agreement with [START_REF] Weinstein | Fluids and Plasmas : Geometry and Dynamics[END_REF]. In the general case K = 0, it is shown likewise in Appendix B that

B(B K ) -1 = i β C F 1 1 -exp(-i β C F) . ( 94
)
This ratio stems therefore from a property inherent to the Lie group underlying our variational approach, namely the exponentional form of the automorphism [START_REF] Mermin | [END_REF]. The Bose-like factor that enters B arises from this property, and from the simplicity of B K . [The quasi-boson structure of the static thermal RPA for fermions [17] appears as a special case, see Sec. VIII.]

Thermodynamic quantities and expectation values were determined in Sec. IV; they depend on the mean-field images of K and H within the Lie algebra, namely, K (0) which self-consistently (Sec. IV A) determines {R (0) } and H which governs (Sec. IV B) the time-dependence of {R (0) (t)}. More elaborate ingredients are required for the evaluation of static and dynamic correlation functions: The matrix K enters B through F while the matrix H (together with the mean-field Liouvillian L) governs the time-dependence of the approximate Heisenberg observables Q H j (t ′ , t). In the rest of this article we will analyse the properties of the approximate correlation functions [expressed by Eqs.( 84),( 88) and (67)] and of the dynamics [expressed by Eqs.(62)].

VI. PROPERTIES OF APPROXIMATE CORRELATION FUNCTIONS AND FLUCTUATIONS

In this section we review some consequences of the variational expressions found above for the correlation functions, encompassing special cases and conservation laws.

A. Time-dependent correlation functions in an equilibrium initial state

We first consider the special case of an equilibrium initial state D = exp(-βK) for which the operator K is equal to the Hamiltonian H, plus possibly some constants of motion such as the particle number for a grand-canonical equilibrium. One can then generate simply the dynamics by using K instead of H in the backward Heisenberg equation (10).

The exact expectation values Q j t then do not evolve. Their approximation (61) depends on time through {R (0) (t)} governed by (63). In the effective Liouvillian L defined by (64), the coordinates H β of the image of the Hamiltonian are replaced by those of the image K of K, that is, K γ {R} = ∂k{R}/∂R γ . Using the self-consistent equilibrium condition (50) for D(0) and the identity [START_REF] Cloizeaux | Many-Body Physics[END_REF] for the product CS, we find

L α β {R (0) } = Γ α βγ K γ {R (0) } = -T Γ α βγ J (0)γ = -T C βγ S γα . ( 95
)
The identity (54) then entails, as expected, that R (0)

α (t) = R (0)
α , and hence that the approximation Q j t does not depend on time.

For two observables Q j and Q k that do not commute with K, the exact two-time correlation function C jk (t ′ , t ′′ ) defined by (7) depends only on the time difference t ′ -t ′′ . This property is not obvious for the approximation (84), which involves two factors depending separately on t ′ and t ′′ . To elucidate this point, let us solve for H = K the approximate backward Heisenberg equation (67) for the image

Q H j (t ′ , t) = Q H j (t ′ -t) with the boundary condition Q H j (0) = Q j {R (0)
}. This equation involves the kernel L + CH which, according to (95), takes the simple form

C(K -T S) = CF.
Hence, the Heisenberg equation (67) simplifies into dQ Hα

j (t ′ , t)/dt = -Q Hβ j (t ′ , t) (C F) α β
, where C and F are evaluated for {R (0) }; it is readily solved as

Q Hα j (t ′ , t i ) = Q β j {R (0) } e C F (t ′ -ti) α β (96)
(with α, β = 0) in terms of the boundary condition Q β j {R (0) } = ∂q j {R (0) }/∂R (0)

β . In order to evaluate the correlation function (84), use is also made of

Q Hα k (t ′′ , t i ) = e -F C (t ′′ -ti) α β Q β k {R (0) } , (97)
and of the relation

F -1 e -F C (t ′′ -ti) = e -C F (t ′′ -ti) F -1 . ( 98
)
The explicit dependence on time of the correlation functions is then given for t ′ > t ′′ by

C jk (t ′ , t ′′ ) ≃ Q α j {R (0) } e C F (t ′ -t ′′ ) i C F 1 1 -exp(-i β C F) F -1 αβ Q β k {R (0) } . ( 99
)
The identity of H and K has led to the occurrence of the same matrix C F in the dynamical equation ( 96) and in the matrix B [Eq.( 88)] that accounts for the correlations in the initial state. As a consequence, only t ′ -t ′′ appears in (99), as it should. This property would not have been satisfied if the correlation matrix B were naively evaluated, as in Eq.( 90), by replacing in (89) the exact state D by D(0) .

For j = k, C jj (t ′ -t ′′ ) provides the variational approximation for the autocorrelation function of the observable Q j in the equilibrium state D ∝ exp(-βK).

B. Other special cases

Commutators and linear responses

The antisymmetric part of B, namely i C/2, is simple and depends only on the zeroth order in the sources, not on F. As a consequence the approximation B αβ -B βα for the expectation value Tr [M α , M β ] D = i Tr Γ γ αβ M γ D of the commutator [M α , M β ] is obtained as i C αβ {R (0) }, a property in agreement with the expectation value R (0) γ of M γ found at first order in the sources.

Linear responses are expectation values of commutators, and therefore involve only this antisymmetric part of C jk (t ′ , t ′′ ). Alternatively, they can be evaluated directly from the variational expression Ψ{A, D} by including a time-dependent perturbation in the Hamiltonian. Then, the response of Q j to a perturbation Q k appears as an expectation value, and is hence directly obtained at first order in the sources, without the occurrence of K which enters the symmetric part of C jk (t ′ , t ′′ ). Both approaches yield

χ jk (t ′ , t ′′ )= (1/i )θ(t ′ -t ′′ )[C jk (t ′ , t ′′ ) -C kj (t ′′ , t ′ )] ≃ θ(t ′ -t ′′ ) Q Hα j (t ′ , t) C αβ {R (0) (t)} Q Hβ k (t ′′ , t) , ( 100
)
where t is an arbitrary time in the interval t i ≤ t ≤ t ′′ and θ the usual step function. In particular, by letting t = t i in (100), the responses are variationally expressed in the Heisenberg picture in terms of the matrix C{R (0) } and of the Heisenberg observables Q H j (t ′ , t i ) and Q H k (t ′′ , t i ) given by the approximate backward Heisenberg equations (67). For H = K, D ∝ exp(-βK), the response (100) in an equilibrium state depends only on the time difference and takes the form

χ jk (t ′ , t ′′ )= (1/i ) θ(t ′ -t ′′ ) Tr D e i H (t ′ -t ′′ )/ Q j e -i H (t ′ -t ′′ ) , Q k ≃ θ(t ′ -t ′′ ) Q α j {R (0) } e C F (t ′ -t ′′ ) C αβ Q β k {R (0) } . ( 101
)

Static correlations; classical limit

Static correlations between observables Q j and Q k in the state D ∝ e -βK are variationally obtained by letting t ′ -0 = t ′′ = t i in (84), which yields

C jk (t i + 0, t i ) = Q α j {R (0) } B αβ Q β k {R (0) } . (102)
The matrix B of correlations between the operators {M} is here saturated by the coordinates Q α j {R (0) } = ∂q j {R (0) }/∂R (0) α of the images in the Lie algebra of the considered observables. For an Abelian algebra, the matrix C vanishes and the ratio (94) reduces to unity. The ordinary and Kubo correlations are identical, and the matrix B simplifies into B K = (β F) -1 . In the high-temperature limit β → 0 and in the classical limit → 0, the ratio (94) also tends to I, so that

B → (βF) -1 .
If Q j and Q k are commuting conserved observables, the occurrence of the commutation matrix C in this ratio also reduces it to I.

For a Curie-Weiss model of interacting spins σ j = ±1 at equilibrium, the Weiss mean-field expressions for the thermodynamic properties and the expectation values σ j are recovered, while the Ornstein-Zernike approximation [START_REF] Ornstein | Proc. Acad. Sci. Amsterdam[END_REF] for correlations is recovered from B = (β F) -1 .

For the fermionic single-particle Lie algebra, the matrix B is the variational approximation for the correlations of the operators a † µ a ν and a † σ a τ :

Tr (a † µ a ν )(a † σ a τ ) D -Tr a † µ a ν D Tr a † σ a τ D ≃ B νµ,τ σ . ( 103
)
If D were replaced by an independent-particle state D, the expectation value Tr (a † µ a ν )(a † σ a τ ) D would be given by Wick's theorem and Tr a † µ a ν D Tr a † σ a τ D would be the Hartree term so that the left-hand side of (103) would reduce to the Fock term Tr(a † µ a τ D) Tr(a ν a † σ D). For a more general state D, the expression (88) of B involves the static RPA kernel iCF, which takes into account not only the Fock term but also, through K, effects of the interactions present in K.

Fluctuations

The static fluctuation ∆Q j of the observable Q j in the state D is variationally given by (102) where k = j, that is,

∆Q 2 j = Q α j {R (0) } B αβ Q β j {R (0) } . (104)
This fluctuation ∆Q j (t) evolves in time according to (84) with k = j, t ′ = t ′′ = t, that is,

∆Q 2 j (t) ≃ Q Hα j (t, t i ) B αβ Q Hβ j (t, t i ) , (105)
which involves only the symmetric part of (88). For H = K, the fluctuation is time-independent. For H = K, examples of time-dependences that are not properly accounted for by non-variational mean-field approximations are given at the end of Sec. VI A and in Sec. VI C.

Initial state in the Lie group

In case the exact density operator D belongs to the trial Lie group, the operator K belongs to the Lie algebra and coincides with its image, D(0) equals D; the matrix K vanishes, F reduces to -T S and B to the trivial form (90). [For fermions, B reduces to the Fock terms.] The quantities Q Hα j (t ′ , t i ) and Q Hβ k (t ′′ , t i ) are in this case the only ingredients, apart from D(0) , that enter C jk (t ′ , t ′′ ).

Zero-temperature limit

The present formalism encompasses ground-state properties, found by letting β → ∞. This limit entails simplifications. [For instance, in the fermionic case, the parameters {R (0) } of D(0) constitute a matrix satisfying [R (0) ] 2 = R (0) .] The number of vanishing eigenvalues of the commutation matrix C increases. While S{R (0) } tends to 0 as T → 0, the quantities β -1 ∂S/∂R (0) α , β -1 S and CF remain finite, due to the singularity of the von Neumann entropy [START_REF] Emch | Algebraic Methods in Statistical Mechanics and Quantum Field Theory[END_REF] for vanishing eigenvalues of D. The resulting simplifications of the correlation matrix B will be exhibited below in the diagonalized form (135) of B.

A further simplification occurs if the initial state lies in the Lie group. For fermions (possibly with pairing) the initial state is then a Slater determinant (or a BCS state). In this case, it has been shown [6] that one can by-pass the solution of the equations (67) for Q Hα j (t ′ , t i ), of the RPA type. The proof relies on the fact that these equations (67) for Q Hα j (t ′ , t i ) involve the same kernel as the dynamical equations (165) for small deviations δR (0) α (t), and that the latter equations can in practice be worked out by expansion of the simpler time-dependent mean-field equations (62) for R (0) α (t). Two-time correlation functions C jk (t ′ , t ′′ ) and time-dependent fluctuations ∆Q j (t) can thus be evaluated by running the existing TDHF (or TDHFB) codes alternatively forward and backward, with appropriate shifts in the boundary conditions. (For another derivation, see [4,11,[START_REF] Broomfield | [END_REF].) This technique, variationally consistent, has been successfully applied [7][8][9][10] to describe, in nuclear systems, the fluctuations of single-particle observables which were severely underestimated by the conventional use of TDHF; for a review, see [11].

C. Images of Heisenberg operators; conservation laws

It has already been noted that the equations of motion (67) for the coordinates Q Hα j (t ′ , t) of A (1) (t) appear as a variational counterpart of the backward Heisenberg equation (10). To be more precise, let us write the time dependence of the Lie algebra operator

Q H j (t ′ , t) = Q Hα j (t ′ , t) M α . Using the relation (L + C H) α β = Γ α βγ H γ {R (0) (t)} + C βγ {R (0) (t)} H γα {R (0) (t)} = ∂(C βγ H γ )/∂R (0)
α (t) and the definition (33) of images, one can rewrite the equations of motion for Q Hα j (t ′ , t) (including α = 0) as

dQ H j (t ′ , t) dt = Image of - 1 i [Q H j (t ′ , t), H] = - 1 i [Q H j (t ′ , t), H] , (106)
the image H of the Hamiltonian H being evaluated with respect to {R (0) (t)}. While the equations (67) were written in terms of the coordinates Q Hα j of Q H j , the introduction of images gives these equations a simple operator form: The right-hand side is simply the image of the r.h.s. of the exact backward Heisenberg equation (10) for time-dependent observables. The boundary condition Q H j (t ′ , t ′ ) = Q S j is also the image of the boundary condition for the Heisenberg observable Q H j (t ′ , t). Up to now, no specific assumptions have been made about the data. Let us now consider an observable Q S j that belongs to the Lie algebra and commutes with H. This conservation law, together with the equation of motion (106), shows that the operator Q H j (t ′ , t) is constant and equal to Q S j . Hence, the approximate expectation value Q j t , evaluated through (70) for t ′ = t i , is constant. Less trivially, the fluctuation ∆Q j (t), evaluated through (105), is also constant as it should. This property, which arises naturally through the use of the approximate backward Heisenberg equation, was not granted. [For instance, in the time-dependent Hartree-Fock approximation for fermions, fluctuations evaluated from D (0) (t) through Wick's theorem are not constant for a conserved single-particle observable.] Another conservation property holds for two observables Q S j and Q S k of the Lie algebra such that

Q S j is conserved, [Q S j , H] = 0, and that [Q S k , H] = i Q S j .
This occurs for instance if Q S k = X is a coordinate of the center of mass of the system, and Q S j = P/m the associated velocity operator (m is the total mass). Then, Eq.( 106) implies that, as the exact Heisenberg operators, the approximate ones satisfy P H (t ′ , t) = P S and X H (t ′ , t) = X S + (t ′ -t)P S /m (see Sec. 5.5 of [4]). Hence, not only X t and P t , but also the fluctuations and correlations of X and P produced by the general formula (84) have the proper time dependence. In particular, ∆P is constant and ∆X 2 (t) satisfies d ∆X 2 (t)/dt = ∆P 2 /m 2 ; we thus acknowledge that the present approximation for the fluctuations accounts for the exact spreading of the wave packet [whereas the TDHF approximation produces a constant width].

VII. MAPPED LIE ALGEBRA AND MAPPED HILBERT SPACE

The purpose of this section is to rewrite in a unified form the variational approximations (59),(61) for expectation values and (84),( 88),(99) for correlation functions . To this aim, we will rely on a correspondence that associates with the Lie algebra {M} a simpler Lie algebra {M}. It will turn out that the various expressions found above can be re-expressed as traces over a single effective density operator D acting in a new mapped space H rather than in the original Hilbert space H .

A. Unifying the approximate expectation values and correlations in the initial state

The optimization of the expectation value M α = Tr D M α has yielded the approximation

M α app = R (0) α = Tr D(0) M α ; in contrast the optimization of M α M β = Tr D M α M β has yielded M α M β app = B αβ + R (0) α R (0)
β that cannot be expressed as a trace over a density operator in the Hilbert space H . We wish to map the set {M} acting in H onto a new set {M} acting in a new space H , and to introduce in H an effective density operator D so as to re-express our approximations in terms of D. Namely, we wish the exact expectation values over D in the mapped space H , denoted as M α map and M α M β map , to coincide with the corresponding variational approximations in the original space H , denoted as M α app and M α M β app , according to

M α map ≡ Tr M α D = M α app = R (0) α , (107)
M α M β map ≡ Tr M α M β D = M α M β app = B αβ + R (0) α R (0) β . ( 108
)
Going from the space H to H will be a price to pay for expressing both expectation values and correlations in terms

[M α , M β ] = i Γ γ αβ R (0) γ M 0 = i C αβ . (109)
(From now on we shall most often drop, as in the end of (109), the unit operator M 0 .) Multiplication of any number of operators M α generates an enveloping algebra, the space of representation of which defines H . The structure of this space will be cleared up in Sec. VIII by setting C αβ into a canonical form.

In order to satisfy the conditions (107) on expectation values M α app the sought effective density operator D should depend only on the differences M α -R

α . As regards the conditions (108), we remember that correlations are often generated in statistical mechanics by a probability distribution having the form of an exponential of the free energy regarded as a function of the running variables. For instance an energy distribution is the product e -β F (E) = e -β E+S(E) of the Boltzmann-Gibbs exponential e -β E by the level density, the exponential e S(E) of the entropy, which accounts for the discarded variables. Another example was provided in our approximation by the thermodynamic coefficients given by the second derivative (56) of the free energy f {R} at {R} = {R (0) }. Here, likewise, the reduction from {M} to {M} suggests to rely on an effective free-energy operator rather than on an effective Hamiltonian. We therefore replace, in the free-energy function f {R} ≡ k{R} -T S{R} defined by (53), the variables {R} by the corresponding new operators {M}. The operators M α fluctuate around R (0) α , and this leads us to expand the operator f {M} as

f {M} ≈ f {R (0) } + 1 2 (M α -R (0) α ) F αβ (M β -R (0) β ) + ... , (110)
where the first-order term is absent thanks to the stationarity of f {R} at {R (0) }.

We thus guess that the distribution D governing the operators {M}, which is expected to yield the identities (107)-( 108), should have in the mapped space H the exponential form

D ≡ e -β F Tr e -β F , F ≡ 1 2 (M α -R (0) α ) F αβ (M β -R (0) β ) , (111)
where F behaves as a kind of free-energy operator. The inclusion of an entropic contribution in the exponent of D accounts for the elimination of degrees of freedom associated with the replacement of the original Hilbert space H by the space H that involves a smaller set of observables, those generated by the operators {M} and their products. The non-negativity of the operator F is ensured by that of the matrix F. More concrete interpretations of F will be given in Sec. VIII by Eqs.( 139) or (144).

While

M α map = R (0)
α is evident, the surmise (108) is proved in Appendix B. It is first shown there that the Kubo correlations of the operators {M} in the state D are given by 1

β β 0 dτ e τ F (M α -R (0) α ) e -τ F (M β -R (0) β ) map = 1 β F -1 αβ . ( 112
)
The ordinary correlation matrix of the operators {M} is then derived from (112) and shown to coincide with the matrix B defined by ( 88):

(M α -R (0) α )(M β -R (0) β ) map = i C F 1 1 -exp[-i β C F] F -1 αβ = B αβ . ( 113
)
Thus, the operator D can be regarded as a substitute in H to the exact state D ∝ exp(-βK) that satisfies the anticipated identities ( 107) and ( 108): We can interpret the matrix elements of B as exact correlations of the operators {M} in the state D. At first order in the sources, the variational treatment amounted to replace the operators {M} by their expectation values {R (0) }. Here, at second order, we reproduce the variational approximation B for correlations by keeping contributions of lowest order in the deviation {M -R (0) }, both as regards the reduction (109) of the algebra of {M} into that of {M} and as regards the expansion (110) of the free energy operator.

B. Heisenberg dynamics of the mapped Lie algebra

Let us extend the above results to the time-dependent correlation functions. We restrict here to the case where H = K, as in Sec. VI A. As already seen, the solution of the approximate backward Heisenberg equation is then generated by the kernel iCF according to Eqs.(96) for the coordinates α = 0 [and to Eq.( 70) for α = 0]. This time-dependence keeps the Heisenberg operators M H α (t ′ , t) in the Lie algebra {M} since they are given by

M H α (t ′ , t) -R (0) α = e C F (t ′ -t) β α (M β -R (0) β ) , (114)
or equivalently by

d M H α (t ′ , t) d t = -(C F) β α (M H β (t ′ , t) -R (0) β ) . ( 115
)
The equations (115) constitute a linear set. We are thus led to define, in the mapped space H , the time-dependent operators M H α (t ′ , t) by the corresponding equation

d M H α (t ′ , t) d t = -(C F) β α (M H β (t ′ , t) -R (0) β ) , (116)
with the boundary condition M H α (t, t) = M α . Moreover, from the definition (111) of the operator F and from the mapped algebra (109), it follows that

[M α , F ] = i (C F) β α (M β -R (0) β ) . ( 117
)
Hence the equations of motion of the set {M H (t ′ , t)} defined by (116) read alternatively

d M H α (t ′ , t) d t = - 1 i M H α (t ′ , t), F . (118)
One recognizes the structure of an exact backward Heisenberg equation (10) in the space H , where the free-energy operator F plays the role of a Hamiltonian. As in the static case, the mapping of the operators {M} onto {M} replaces approximate properties in H by exact ones in H , here the approximate dynamical equation (115) in the space H by (118) where the commutator is restored in the space H . Accordingly, the mapped Heisenberg operators M H α (t ′ , t) are given by

M H α (t ′ , t) = e i F (t ′ -t)/ M α e -i F (t ′ -t)/ , (119)
a mere unitary transformation in the space H , whereas the transformation (114) of the set {M} in the space H is not unitary. The linearity of Eqs.(116), or equivalently the quadratic nature of F , exhibits harmonic-oscillator dynamics in the mapped space. This will be made more precise in Secs. VIII D and VIII E. In Eq.(118) the effective "Hamiltonian" F in the space H should not be confused with an approximation for the Hamiltonian H = K of the original problem. Its expressions (111) and ( 56) are related to the free-energy function f {R} and its deviations, rather than to the original Hamiltonian H. The contribution to F of the entropic term -β -1 S, through F = K -β -1 S, is essential. This contribution is the only one left if K belongs to the Lie algebra, and it remains finite at zero temperature.

C. Unified formulation of the variational expressions

Our mapping has provided a formalism (Sec. VII A) in which both optimized expectation values and correlations in the state D are generated, in the mapped space H , as traces over the effective density operator D ∝ exp(-βF ). This operator depends on the temperature both explicitly and through the kernel F of the operator F .

An arbitrary observable Q j acting in the original Hilbert space is now represented by its mapped image given, according to [START_REF] Ring | The Nuclear Many-Body Problem[END_REF], by

Q j = q j {R (0) } M 0 + (M α -R (0) α M 0 ) ∂q j {R (0) } ∂R (0) α (α = 0) . ( 120
)
Then, the optimized expectation value of a single operator Q j is

Q j app = q j {R (0) } = Q j map ≡ Tr Q j D , (121)
and that of the product

Q j Q k is Q j Q k app = Q j Q k map ≡ Tr Q j Q k D , (122)
both being directly obtained from the effective state D in the space H .

In the case H = K considered in Secs. VI A and VII B, the time-dependence of the Heisenberg operator associated with Q j is governed by the backward Heisenberg equation (118), which yields

Q H j (t) = e i F t/ Q j e -i F t/ . (123)
Hence, the variational two-time correlation functions take in the mapped space H the simple form

C jk (t ′ , t ′′ ) = T Q H j (t ′ -t ′′ ) Q k map -Q j map Q j map . (124)
The same operator F [Eq.( 111)] occurs both in the density operator D and as an effective Hamiltonian in (123). It depends on the following ingredients: R (0) α defined self-consistently by Eq.( 51), C αβ {R (0) } given by Γ γ αβ R

γ , F given by (56) and D by (111).

VIII. THE EIGENMODES AND THEIR INTERPRETATION

In this section, we consider only initial states at equilibrium, in which case the dynamics is governed by the Hamiltonian H = K. The variational expressions of correlation functions obtained in Sec. V then involve functions of the sole matrices C and F. Practical evaluations rely on their diagonalization. The eigenvalues and eigenvectors thus obtained will enable us to interpret the above results.

A. Diagonalization of the evolution kernel and the correlation matrix

The product iCF of the commutation matrix C and the stability matrix F is the kernel which governs the evolution, as exhibited by the dynamical equations (96) of Sec. VI A. It also occurs in the expression (88) of the correlation matrix B. The diagonalization of the matrix iCF, and the study of its properties, appear therefore appropriate.

The Lie algebra is globally hermitian, namely, if the operators M α of the basis are not individually hermitian, they come in conjugate pairs (e.g., a † µ a ν and a † ν a µ in the fermionic example). We denote as α the index of the operator M α ≡ M † α , which may or may not differ from M α . The change α → α of all the indices in the quantities R α , Q α , C αβ or F αβ transforms them into their complex conjugates. Hence, iC, which is antisymmetric, and F, which is symmetric, are equivalent to hermitian matrices. Moreover, F is non negative. (Only the case of a strictly positive matrix F is examined here; vanishing eigenvalues are considered in the end of Sec. IX C.) Altogether, the matrix iCF is equivalent to the antisymmetric hermitian matrix iF 1/2 CF 1/2 , so that its right and left eigenvectors, respectively denoted as ψ and φ, constitute a complete biorthonormal basis in the space α, while its eigenvalues are real and either vanish or come in opposite pairs.

Using the above properties, one can classify into three subsets the eigenvectors of iCF, denoted with indices n, -n and p, respectively:

(i) To the positive eigenvalues Ω n > 0 are associated the right and left eigenvectors ψ n and φ n defined by (α,

β = 0) i C αγ F γβ ψ n β = Ω n ψ n α , (125)
φ α n i C αγ F γβ = φ β n Ω n , φ β n = (ψ n α ) * F αβ . (126)
(ii) Taking the complex conjugate of (125) provides the right and left eigenvectors ψ -n and φ -n associated with the negative eigenvalue -Ω n :

ψ -n α = (ψ n α ) * , φ α -n = (φ α n ) * . (127)
(iii) The eigenvectors ψ p and φ p associated with vanishing eigenvalues are given by

i C αγ F γβ ψ p β = 0 , φ α p C αβ = 0 , (128)
ψ p α = (ψ p α ) * , φ β p = (φ β p ) * = (ψ p α ) * F αβ . ( 129
)
The biorthonormality of the sets {ψ} and {φ} is equivalent to the orthonormalization relations for the right eigenvectors {ψ} expressed by

φ β n ψ n ′ β = (ψ n α ) * F αβ ψ n ′ β = δ nn ′ , φ β p ψ n ′ β = (ψ p α ) * F αβ ψ p ′ β = δ pp ′ , ψ n α F αβ ψ n ′ β = ψ n α F αβ ψ p β = (ψ n α ) * F αβ ψ p β = 0 . (130)
Likewise, the closure property is equivalent to

(F -1 ) αβ = n [ψ n α (ψ n β ) * + (ψ n α ) * ψ n β ] + p ψ p α ψ p β . (131)
The matrix iCF is thus diagonalized as

i(C F) β α = n Ω n [ψ n α φ β n -(ψ n α ) * (φ β n ) * ] , (132)
a sum involving only the modes n associated with Ω n > 0. The matrices F and C can then be expressed in the form

F αβ = n [(φ α n ) * φ β n + φ α n (φ β n ) * ] + p φ α p φ β p , (133)
i C αβ = n Ω n [ψ n α (ψ n β ) * -(ψ n α ) * ψ n β ] , (134)
which for F includes also terms associated with the vanishing eigenvalues of iCF.

The above diagonalization appears as a generalization, for an arbitrary Lie group and at non-zero temperature, of the standard diagonalization of the RPA matrix for fermion systems [START_REF] Ring | The Nuclear Many-Body Problem[END_REF][START_REF] Blaizot | Quantum Theory of Finite Systems[END_REF]. The reality of the eigenvalues Ω n , which follows, as shown above, from the positivity of the matrix F is well known in that case [START_REF] Thouless | [END_REF][START_REF] Cloizeaux | Many-Body Physics[END_REF][START_REF] Mermin | [END_REF]. The present diagonalization is also similar to the diagonalization of a quadratic Hamiltonian of boson operators [START_REF] Blaizot | Quantum Theory of Finite Systems[END_REF].

Inserting the expansions (133) and (134) of F and C into the expression (88) of the correlation matrix B yields

B αβ = n ψ n α Ω n 1 -exp(-β Ω n ) (ψ n β ) * (135) + n (ψ n α ) * Ω n exp(β Ω n ) -1 ψ n β + 1 β p ψ p α ψ p β .
This expression involves only the eigenvalues and right eigenvectors of iCF, normalized according to (130). In the zero-temperature limit, only the first term of (135) survives, and B reduces to

B αβ = n ψ n α Ω n (ψ n β ) * , (136)
where the parameters Ω n and ψ n α are found from the limit β → ∞ of iCF. The quasi-scalars do not contribute. In the high-temperature limit or in the classical limit, B = (βF)

-1 is given in diagonalized form by (131). The same holds for the Kubo correlations.

B. Diagonalization of the effective free-energy operator

The above diagonalization of the matrix F will help us to give an interpretation of the effective Hamiltonian F defined in the mapped space H [Eq.(111) of Sec. VII A]. Let us introduce, besides the unit operator M 0 = I, the following new basis for the mapped algebra {M} (α = 0):

A n ≡ φ α n (M α -R (0) α ) √ Ω n , A -n = A † n , (137)
Y p ≡ φ α p (M α -R (0) α ) β = (Y p ) † .
The pair of operators A n , A -n is associated with each mode Ω n , and the single operator Y p with each vanishing eigenvalue of iCF. Conversely, the original operators of the mapped Lie algebra are decomposed over the modes n and p according to

M α -R (0) α = n Ω n [ψ n α A n + (ψ n α ) * A † n ] + p ψ p α Y p / β . (138)
The eigenvectors ψ appear as the amplitudes of M α on the different modes.

The free-energy operator F takes in the new basis the diagonal form

F = n Ω n 1 2 A † n A n + A n A † n + p 1 2β Y 2 p . (139)
The commutation relations of the operators A n , A † n and Y p follow from those ([M α , M β ] = i C αβ ) of the operators {M}, from the diagonal form (134) of the commutation matrix C αβ and from the biorthogonality (130) of the amplitudes ψ and φ. The transformation (137) thus implies

[A n , A † n ′ ] = δ nn ′ , ( 140
) [A n , A n ′ ] = [A n , Y p ] = [Y p , Y p ′ ] = 0 .

C. A bosonic and scalar algebra

The mapped algebra {M} is spanned, according to (138), by two sets of operators. On the one hand, the commutation relations (140) express that this algebra has a symplectic sector n: the set {A, A † } is simply an algebra of bosonic operators, with single-boson states labelled by the index n. On the other hand, in the sector p, the set {Y} commute with all other operators of the mapped algebra, and hence each Y p can be regarded as a scalar random variable.

The form (139) of the free-energy operator F in the space H is suggestive. It can be identified with a Hamiltonian of non-interacting bosons, with single-particle states n having the energy Ω n , plus a classical quadratic part with a coefficient 1/2β for each variable Y p .

The evaluation of expectation values ... map involves traces over the space H with the weight D ∝ exp(-βF ). As regards the bosonic part of (139), H contains a Fock space and the trace is as usual a summation over the occupation numbers of each single-boson state. As regards the scalar part, the trace is meant as the integration p dY p . Thus, the expectation values over D ∝ exp(-βF ) of single operators A n , A † n or Y p vanish, while we find for pairs

A † n A n map = 1 e β Ωn -1 = A n A † n map -1 , (141)
Y p 2 map = M 2 0 map = 1 ;
all other expectation values of pairs vanish. Each mode n yields a Bose factor associated in canonical equilibrium with the energy Ω n ; no chemical potential occurs for the distribution D ∝ exp(-βF ). For the scalar variables, the variance 1 arising from the Gaussian weight exp(-Y p 2 /2) within D provides an expectation value 1/2β for the term Y p 2 /2β of F , in agreement with the equipartition theorem of classical statistical mechanics. We are considering in this section the dynamics (119) for H = K, and we have seen that F then plays the role of a Hamiltonian in the space H . The diagonal form (139) of F indicates that the mapped operators A n oscillate as

A H n (t ′ , t) = e -iΩn(t ′ -t) A n , (142)
while the scalars Y p remain constant. Thus, the operators A n , A † n behave as bosonic annihilation and creation operators in all respects: commutation relations (140), average occupation in canonical equilibrium (141) and dynamics (142).

D. Interpretation in terms of oscillators

The mapped bosonic Fock space can equivalently be regarded as a space of oscillators, each single bosonic state n corresponding to an oscillator mode with frequency Ω n . The operators A n , A † n are thus replaced by the position and momentum operators

X n = /2 (A n + A † n ) , P n = /2 (A n -A † n )/i , (143)
satisfying the canonical commutation relations [X n , P n ′ ] = i δ nn ′ . The effective Hamiltonian F then takes the form

F = n Ω n (P n 2 + X n 2 ) + 1 2β p Y p 2 . ( 144
)
It now describes uncoupled harmonic oscillators, plus a quadratic energy associated with the classical variables Y p . In this alternative interpretation, we have

X 2 n map = P 2 n map = 2 coth β Ω n 2 , (145)
and the dynamics describes harmonic oscillators with frequency Ω n .

E. Quasi-bosons and quasi-scalars

From now on, we return to the original space H . The correspondence between the algebras {M} and {M} leads us to perform on the original Lie algebra the same transformation as (137)-(138), which introduces for a given D (0) a new basis A n , A † n , Y p (besides M 0 ) in the algebra {M}. The commutators between these operators, issued from the Lie structure (21), are not simple. However, the expectation values of these commutators, evaluated as traces over D(0) , have the same structure as the commutators (140) of the corresponding mapped operators

A n , A † n , Y p , that is, [A n , A † n ′ ] app = δ nn ′ , ( 146
) [A n , A n ′ ] app = [A n , Y p ] app = [Y p , Y p ′ ] app = 0 .
The operators A n , A † n can thus be termed "quasi-boson" annihilation and creation operators, while the operators Y p can be termed "quasi-scalars".

The matrix i C αβ was defined as the expectation value [M α , M β ] app in the original basis of the Lie algebra. Going to the new basis amounts to diagonalize C according to (134). In this new basis, i C takes a form involving only 2 × 2 diagonal blocks 0 1 -1 0 and zeros, in agreement with the above relations.

The variational expectation values in the state D(0) of the operators A n , A † n and Y p vanish,

A n app = A † n app = Y p app = 0 , (147)
while M 0 app = 1. The approximate expectation values of pairs of operators,

A † n A n app = 1 e β Ωn -1 = A n A † n app -1 , (148)
Y p 2 app = M 2 0 app = 1 ,
are the same as the equations (141) which were exact in the space H ; again all other expectation values of pairs vanish. We recognize, as in the mapped space, the Bose factor of a mode n for the quasi-boson operators A n , A † n , and the unit variance for the Gaussian quasi-scalar variables Y p .

The change of basis {M} → {A, A † , Y} is issued from the diagonalization of Sec. VIII A. This change leads to re-express the correlation matrix B as

B αβ = M α M β app -R (0) α R (0) β = n ψ n α Ω n A n A † n app (ψ n β ) * (149) + n (ψ n α ) * Ω n A † n A n app ψ n β + 1 β p ψ p α Y p 2 app ψ p β .
The Bose-like factors that appeared in the expression (135) of B now come out as expectation values (148) of pairs of quasi-boson operators A n and A † n ′ , while the last term of ( 135) is consistent with the variance equal to 1 of the quasi-scalar operators Y p . The coefficients arise from the change of basis (138) that expresses the operators M α -R (0) α on the basis {A, A † , Y}.

More generally, to write the optimized correlation of arbitrary observables Q j and Q k , one should express their images on the basis A, A † , Y according to

Q j = q j {R (0) } M 0 + n [Q n jb A n + (Q n jb ) * A † n ] + p Q p js Y p , (150)
where the new, bosonic and scalar, coordinates are (α = 0)

Q n jb = Ω n ψ n α Q α j {R (0) } (151)
Q p js = 1 √ β ψ p α Q α j {R (0) } , (we recall that Q α j {R (0) } ≡ ∂q j {R (0) }/∂R (0)
α ). One thus finds for the correlations

Q j Q k app -Q j app Q k app = C jk (t i , t i ) (152) = n coth β Ω n 2 Re Q n jb (Q n kb ) * + i Im Q n jb (Q n kb ) * + p Q p js Q p ks .
The antisymmetric part of the correlations,

1 2 Q j Q k -Q k Q j app = n i Im Q n jb (Q n kb ) * , ( 153
) agrees with [Q j , Q k ] = i Q α j {R (0) } C αβ Q β k {R (0)
}. An alternative interpretation of the symmetric part is obtained by replacing, as in Sec. VIII D, quasi-bosons by quasi-oscillators with variables X n , P n . Writing the image of Q j in the basis {X, P, Y} then produces the factors ( /2) coth(β Ω n /2) which are the expectation values X 2 n app = P 2 n app . Since H = K here, the approximate Heisenberg operators in the space H follow the same evolution as (142) in the mapped space H , that is,

A H n (t ′ , t) = e -iΩn(t ′ -t) A n , Y H p (t ′ , t) = Y p . (154)
Hence, the two-time correlation function C jk (t ′ , t ′′ ) is given, when the initial state is at equilibrium, by

C jk (t ′ , t ′′ ) = n Q n jb e -i Ωn(t ′ -t ′′ ) 1 -e -β Ωn (Q n kb ) * (155) + (Q n jb ) * e i Ωn(t ′ -t ′′ ) e β Ωn -1 Q n kb + p Q p js Q p ks , (t ′ > t ′′ ) ,
which is the explicit form for the general equation ( 124). It exhibits the coefficients of the images Q j and Q k on the quasi-boson and quasi-scalar basis, the Bose factors and the boson dynamics.

F. Linear response and excitation energies

The linear response (100) is directly found from ( 101) and (150); it does not involve the quasi-scalar contributions nor the Bose factor. Its dissipative part, defined through a Fourier transform with respect to t ′ -t ′′ , comes out as

χ ′′ jk (ω) = π n Q n jb δ(ω -Ω n ) (Q n kb ) * -(Q n jb ) * δ(ω + Ω n ) Q n kb . (156)
Thus, the present approximation yields the Ω n 's as resonance frequencies. Consistency properties such as the Kramers-Kronig dispersion relations and the Kubo fluctuation-dissipation relations are satisfied. At the zero-temperature limit, the exact expression of χ ′′ jk (ω) is given, in terms of the ground state |0 of K and of the excited states |exc with excitation energies E exc , as

χ ′′ jk (ω) = π exc [ 0|Q j |exc exc|Q k |0 δ( ω -E exc ) (157) -0|Q k |exc exc|Q j |0 δ( ω + E exc )] .
Comparison with (156) shows that the positive eigenvalues Ω n of i CF can be identified as approximations for the excitation energies E n of some set of states, labelled as |n . The amplitude Q n jb of the image Q j over the quasi-boson annihilation operator A n appears as an approximation for the matrix element 0|Q j |n . In particular, the quasi-boson operators A n satisfy approximately

0|A n |n ′ = δ nn ′ , (158)
as if the state |n were obtained by creating a quasi-boson.

IX. SMALL DEVIATIONS

We consider in this section static and dynamic changes brought in by a shift of the initial conditions.

A. Static deviations

The change in the equilibrium properties arising from a variation δK of the operator K that defines the state D∝ exp(-βK) is accounted for in the variational treatment by the shift δ D(0) of D(0) , parametrized by the set

δR (0) α = Tr δ D(0) M α = Tr[ D(0) + δ D(0) ](M α -R (0) α ) . ( 159
)
Denoting by δK the image of δK, the first-order change in the self-consistent equations (51) provides

S αβ δR (0) β = β δK α + β K αβ δR (0) β , (160)
that is, using the definition (56) of the matrix F,

δR (0) α ≡ -(F -1 ) αβ δK β (α, β = 0) . (161)
The resulting change ∆F of the free energy F ≃ minf {R} = min[k{R} -T S{R}], expanded up to second order in δK, is found as

∆F ≈ δk + 1 2 δR (0) α F αβ δR (0) β + δK α δR (0) α (162) = δk - 1 2 δR (0) α F αβ δR (0) β = δk - 1 2 δK α (F -1 ) αβ δK β .
In ( 160)-( 162), the matrices K, S, F, the image δK and the symbol δk of δK are taken at {R (0) }. In particular, the expression (57) of the heat capacity is recovered by taking δK = Kδβ/β. This shift takes a suggestive form if F αβ is replaced by its diagonalized form (133). Using ( 159) and (147), and denoting as

δR (0) nb = Tr [ D(0) + δ D(0) ] A n = Tr δ D(0) A n , (163) δR
∆F ≈ δk - 1 2 δR (0) α F αβ δR (0) β (164) = δk - n Ω n |δR (0) nb | 2 - p 1 2β (δR (0) ps ) 2 .
Each oscillator mode provides a contribution Ω n weighted by the amplitude |δR

α (t) dt = (L + C H) β α δR (0) β (t) , (165)
obtained by varying R α (t i ) are given by (159). [For the fermionic singleparticle Lie-algebra, (165) is the dynamical RPA equation issued from the TDHF equation.] The kernel L + CH governing the forward equation (165) for the small deviations turns out to be the dual of the kernel governing the backward equation (67) for the approximate Heisenberg observables Q H j (t ′ , t). The shift δ Q j t , variationally given by δ Q

j t = Q α j {R (0) (t)}δR (0) α (t), is therefore equal to δ Q j t ≃ Q Hα j (t,
{ {R α , R β } } = C αβ {R} ≡ Γ γ αβ R γ (170)
and

{ {g 1 , g 2 } } = ∂g 1 ∂R α C αβ {R} ∂g 2 ∂R β . (171)
One can readily check that the above rules are satisfied, owing in particular to the properties of the structure constants Γ γ αβ , namely, antisymmetry and Jacobi identity. The Poisson structure (170)-( 171) is thus generated by the Lie algebra of the set {M}.

The equations of motion (62) for {R (0) (t)} can now be rewritten as

dR (0) α (t) dt = { {R (0) α (t), h{R (0) (t)}} } . (172)
These quantum variational equations are therefore identified with classical equations involving the brackets (171) and governed by a classical Hamiltonian h{R (0) (t)}, the symbol of the quantum Hamiltonian H. The relation (61) expresses the time-dependent expectation value of the observable Q j t as a function of {R (0) (t)}; together with (172), it implies that Q j t evolves according again to the classical dynamical equation

d Q j t dt = { {q j {R (0) (t)}, h{R (0) (t)}} } , (173)
which involves the symbols of Q j and H. The variational approach, together with the introduction of symbols, thus generate approximate dynamics of expectation values that have a classical structure, generated by the Lie-Poisson bracket (171) and by a Hamiltonian. The symbol h{R (0) (t)} is obviously a constant of the motion. Moreover, the von Neumann entropy S{R (0) (t)} defined by ( 25) is also a constant of the motion. Indeed, using (26) then [START_REF] Thouless | [END_REF], one finds

{ {R (0) β (t), S{R (0) (t)}} } = C βγ {R (0) } ∂S{R (0) (t)} ∂R (0) γ (t) = -C βγ {R (0) (t)} J γ(0) (t) = 0 , (174)
which implies that S{R (0) (t)} remains constant during the evolution (172) of {R (0) (t)}. Lie-Poisson structures for dynamical equations issued from a variational principle have been recognized in cases such as the Vlasov equation [40], time-dependent Hartree-Fock equations [START_REF] Balian | [END_REF], time-dependent Hartree equations for bosons [12] and for φ 4 field theory [15,16].

Proposals of non-linear extensions of quantum mechanics [42] have suggested a formulation in terms of a Poisson structure [43]. Here it is the restriction of the algebra of observables to the trial Lie algebra which produces a Poisson structure within standard quantum mechanics.

B. Canonical variables

We have seen in Secs. VIII B and VIII D that the diagonalization of CF generates in the mapped space H a linear transformation (137),(143) of the operators {M}, which produces pairs of canonically conjugate operators X n , P n and scalars Y p . This corresponds in the original Hilbert space to a construction of quasi-oscillator operators X n , P n and quasi-scalars Y p (Sec. VIII E). Accordingly, the expectation values of their small deviations, defined by

δX n (t) ≡ /2 [δR (0) nb (t) + δR (0) * nb (t)] ≡ /2 Tr δ D(0) (t) (A n + A † n ) , δP n (t) ≡ /2 [δR (0) nb (t) -δR (0) * nb (t)]/i ≡ /2 Tr δ D(0) (t) (A n -A † n )/i ,
evolve according to (169) as conjugate variables of classical harmonic operators with frequency Ω n , while the quantities δY p ≡ δR for the energy (when K is the Hamiltonian):

K ≃ k{R (0) } , (181)
and for the expectation values of the observables Q j in the state D ∝ exp (-βK) (Sec. IV B):

Q j ≡ Tr Q S j e -βK Tr e -βK ≃ q j {R (0) } . (182)
Thermodynamic coefficients, such as the specific heat [Eq.( 57)], are found from the second derivatives

F αβ = ∂ 2 f {R (0) } ∂R (0) α ∂R (0) β = K αβ -T S αβ ≡ ∂ 2 k{R (0) } ∂R (0) α ∂R (0) β -T ∂ 2 S{R (0) } ∂R (0) α ∂R (0) β , (183)
which appear as a matrix in the space of indices α.

While the above expressions appear as mere extensions to arbitrary Lie groups of standard mean-field results, the approach has also provided variational expressions for the correlations in the state D ∝ exp (-βK) (Secs. V D and V E). For correlations between the elements of the Lie algebra, we have obtained (α, β = 0)

M α M β -M α M β ≃ B αβ , B = i C F 1 1 -exp(-i β C F) F -1 , (184)
where the matrix C is defined by

C αβ = Γ γ αβ R (0)
γ . More generally, the correlations between two observables Q j and Q k have been found as

Q j Q k -Q j Q k ≃ ∂q k {R (0) } ∂R (0) α B αβ ∂q k {R (0) } ∂R (0) α , (185)
in terms of the matrix B and of the symbols of the operators Q j and Q k . Fluctuations ∆Q j result for j = k. In the classical or high-temperature limit, B reduces to (βF) -1 .

C. Time-dependent quantities

Here D ∝ exp (-βK) is the exact initial state at the time t i , and the subsequent evolution is governed by the Hamiltonian H. The above results hold as variational approximations at the time t i . For later times, the optimization of the expectation value Q j t yields (Sec. IV B)

Q j t ≡ Tr e -βK e i H t/ Q S j e -i H t/ Tr e -βK ≃ q j {R (0) (t)} , (186)
where the time-dependent expectation values R

α (t) of the operators M α are given by the equations

dR (0) α (t) dt = L γ α {R (0) (t)} R (0) γ (t) (187)
with the initial conditions R (0)

α (t i ) = R (0)
α . The effective Liouvillian L is expressed self-consistently in terms of the symbol h{R} of the Hamiltonian by

L γ α {R} = Γ γ αβ ∂h{R} ∂R β . (188)
(In the presence of non-vanishing structure constants Γ 0 αβ , R

0 (t) should be replaced by 1 in the term γ = 0 of (187).) We have shown (Sec. X) that the dynamical equations (187) and the resulting ones for Q j t have a classical structure, where C αβ {R} = Γ γ αβ R γ appears as the Lie-Poisson tensor associated with the classical coordinates {R} and where h{R} behaves as a classical Hamiltonian.

Here again, non standard variational results have been obtained (Sec. V C) for two-time correlation functions. Their expression involves an approximation (Sec. V A) for the observables in the Heisenberg picture defined by

Q H j (t ′ , t) ≡ e i H(t ′ -t)/ Q S j e -i H(t ′ -t)/ , (189)
where t ′ is the usual final running time and t is a reference time at which Q H j (t + 0, t) reduces to the Schrödinger observable Q S j . This approximation Q

H j (t ′ , t) ≃ Q H j (t ′ , t) ≡ Q Hα j (t ′ , t)M α is characterized by the differential equations dQ Hα j (t ′ , t) dt = -Q Hβ j (t ′ , t) (L + C H) α β (α, β = 0) , (190)
in terms of the reference time t which runs backward from t ′ to the initial time t i . The boundary condition is given by

Q Hα j (t ′ , t ′ -0) = ∂q j {R (0) (t ′ )} ∂R (0) α (t ′ ) . (191)
In Eq.(190), the three matrices L defined by (188), C αβ {R} = Γ γ αβ R γ and H, which denotes the matrix of second derivatives of the symbol h{R} with respect to the variables R α , are evaluated at the point R = R (0) (t).

The optimized expression for the two-time correlation functions C jk (t ′ , t ′′ ) then reads (Secs. V C and V D)

C jk (t ′ , t ′′ ) ≡ Tr e -β K Q H j (t ′ , t i ) Q H k (t ′′ , t i ) Tr e -β K -Q j t ′ Q k t ′′ ≃ Q Hα j (t ′ , t i ) B αβ Q Hβ k (t ′′ , t i ) (t ′ > t ′′ ) , (192)
which involves both the approximate correlations (184) at the initial time and the approximate Heisenberg operators given by Eqs. ( 190) and (191).

D. Properties of the results and consequences

Special cases (Sec. VI) include time-dependent fluctuations, obtained from C jj (t + 0, t), and linear responses

χ jk (t ′ , t ′′ )= (1/i )θ(t ′ -t ′′ )[C jk (t ′ , t ′′ ) -C kj (t ′′ , t ′ )] ≃ θ(t ′ -t ′′ ) Q Hα j (t ′ , t) C αβ {R (0) (t)} Q Hβ k (t ′′ , t) , (193)
where t is arbitrary in the interval t i ≤ t ≤ t ′′ . If the initial state D ∝ e -βK is in equilibrium, with a dynamics generated by H equal to K (or to K plus a constant of motion), the approximate expectation values Q j t remain constant. The equations (191) are solved as

Q Hα j (t ′ , t i ) = Q β j {R (0) } e C F (t ′ -ti) α β . (194)
Hence, two-time correlation functions are variationally given in this case by

C jk (t ′ , t ′′ ) ≃ Q α j {R (0) } e C F (t ′ -t ′′ ) i C F 1 1 -exp(-i β C F) F -1 αβ Q β k {R (0) } , (195)
an expression which, as it should, depends only on the time difference t-t ′ . Derived within a unified framework, these results satisfy consistency properties and conservation laws. For instance, the approximation preserves the relation between the static stability of a thermodynamic equilibrium state (f {R} minimum) and the dynamical stability of the motions (187) around it. Other special cases are also considered in Sec. VI, including classical, high-temperature and zero-temperature limits, as well as Kubo correlations. Results for small deviations are presented in Sec. IX.

In Secs. VII and VIII C-VIII F, we have exhibited, for any Lie group, an attractive interpretation for the static correlations (184) and for the two-time functions (194)-(195). To this aim, the Lie algebra {M} has been mapped into a simpler Lie algebra {M}, the operators M α of which are linear combinations of creation and annihilation bosonic operators A n and A † n (or equivalently of position and momentum operators X n and P n of harmonic oscillators), and of scalar random variables Y p . Then, the dynamics (194) is mapped onto a dynamics governed by an effective Hamiltonian

F = n Ω n 1 2 A † n A n + A n A † n + p 1 2β Y 2 p (196) = n Ω n (P n 2 + X n 2 ) + 1 2β p Y p 2 ,
which describes uncoupled modes. The above correlations (184) and (195) take for an arbitrary Lie group the form of expectation values of pairs of quasi-boson operators and of scalars in a canonical equilibrium state with Hamiltonian F . The explicit expression

C jk (t ′ , t ′′ ) = n Q n jb e -i Ωn(t ′ -t ′′ ) 1 -e -β Ωn (Q n kb ) * (197) + (Q n jb ) * e i Ωn(t ′ -t ′′ ) e β Ωn -1 Q n kb + p Q p js Q p ks , (t ′ > t ′′ )
derived in Sec. VIII relies on the diagonalization (132) of the matrix i C F, which enters the formalism at several places. The coefficients Q n jb and Q p js appear as coordinates, in the new basis of the Lie algebra, of the observable Q j , their explicit form being given by Eqs.(151) in terms of the eigenvalues ±Ω n and eigenvectors of the matrix i C F (which satisfy the equations ( 130), ( 133) and ( 134)). Special cases of (197) are the static correlations written in the diagonalized form (135) and the linear responses, which satisfy the Kramers-Kronig dispersion relations and the Kubo fluctuation-dissipation relation.

XII. CONCLUSION

We have presented above a variational approach to the determination of various physical quantities pertaining to many-body systems. Although our scope has been formal, the generality and flexibility of the treatment appear well suited to many specific problems. The results, listed in Sec. XI, have been obtained by merging several ingredients:

(i)The evaluation of a generating functional (Sec. II A) has allowed the simultaneous optimization of different quantities such as the (static or time-dependent) expectation values and correlations of the observables of interest: they are obtained by expanding in powers of the sources the functional ψ{ξ} ≡ ln TrA(t i ) D. Deriving the correlation functions as second-order contributions in the sources of this generating functional provides for them non-trivial approximations, even for a simple restricted trial space. As an example, for a system of interacting fermions, the use of independent-particle trial objects leads to standard mean-field theories for expectation values but to expressions of the form (192) for correlation functions. It is the dependence of the trial objects on the sources which leads to such elaborate results. Moreover, approximating all quantities in a unique framework preserves consistency properties.

(ii)The variational principle for the optimization of the generating functional is built by means of a general method (Sec. II C). The object to be optimized is characterized by some simple equations, regarded as constraints on the ingredients, namely, on the initial state D ∝ exp(-βK) and on the time-dependent observables of interest Q H j (t ′ , t). Lagrange multipliers are then associated with these constraints.

(iii) For dynamical problems, the Heisenberg observables Q H j (t ′ , t) enter the generating functional, together with the sources ξ j (t), through the "generating operator"

A(t) ≡ T e i ∞ t dt ′ j ξ j (t ′ ) Q H j (t ′ , t) . (198)
The operator A(t) is characterized by the differential equation ( 11) expressing dA(t) dt , which plays the role of a constraint in the variational principle. This Eq.( 11) has been derived as a consequence of the backward Heisenberg equation

dQ H j (t ′ , t) dt = - 1 i [Q H j (t ′ , t), H] . (199)
We recall that, while the standard (forward) Heisenberg equation describes the variation of the Heisenberg observable Q H j (t ′ , t) as function of the running time t ′ , (199) is a differential equation with respect to a reference time t which runs backward from t ′ to the initial time t i . We have explained in Sec. II B why this backward dynamics is the suitable one.

(iv) As regards the initial state D = e -βK , we have characterized it by the Bloch equation

dD(τ ) dτ + K D(τ ) = 0 , (200)
where τ runs from 0 to β. We thus deal with finite temperature. Ground state problems are treated in the limit β → 0; in fact, this procedure turns out to be more convenient than the direct implementation of ground states.

(v) The variational equations associated with the above constraints on A(t) and D(τ ) have been worked out by taking a Lie group as trial space for A(t), D(τ ) and for their associated Lagrange multipliers (Sec. III), thus replacing operators in Hilbert space by functions, their symbols. The approximate states that occur in mean-field theories belong to such Lie groups, for instance static and dynamic states in Hartree-Fock approximations for fermions, Hartree-Fock-Bogoliubov states for fermions with pairing, or coherent states for bosons. However, the trial operators D used here depend on the sources; they are not overall approximations for the exact state D, but serve only to optimize the generating functional (and hence all quantities of interest).

We have recalled in Sec. XI the outcomes of the above approach. Some well-known approximate results have been recovered in special cases, for instance, static and dynamic mean-field or RPA treatments which now appear within a general unified framework based on the use of a Lie group as trial space. Moreover, the new variational results (184), (185), (192) and (195) have been derived for correlations and fluctuations. In particular, we have shown for H = K that quasi-bosons (or quasi-oscillators) come out for any Lie group, and not only in the usual context of zero-temperature RPA for fermions. Let us add a few comments.

In order to determine the generating functional ψ{ξ} ≡ ln TrA(t i ) D, we have characterized A(t i ) and D = e -βK by introducing the functions A(t) and D(τ ) = e -τ K determined from the simple differential equations ( 11) and (200). To build a tractable variational principle, we have been led to introduce the trial time-dependent quantities D(τ ) and A(t) whereas we need only their boundary values D(β) and A(t i ). Moreover, the number of variables is doubled by the introduction of the Lagrange multipliers A(τ ), associated with the equation for D(τ ), and D(t), associated with the equation for A(t). These quantities are not coupled by the exact equations of motion. However, in a restricted trial space, the stationarity conditions entail a coupling between them which allows a better optimization of the generating functional ψ{ξ}.

In spite of their resemblance with density operators and exponentials of observables, the trial objects D(τ ), A(t), A(τ ), D(t) are only computational tools for the evaluation of expectation values and correlation functions. In particular the trial quantity D(t), which looks formally like a density operator in the Schrödinger picture, adapts itself to the question asked, namely to the value of ψ{ξ}. It thus depends on the sources (and is not even hermitian). It cannot be interpreted as an approximate state of the system, and in fact there exists no approximate density operator in the original Hilbert space that would produce both the optimized expectation values (186) and correlations (184), ( 185), (192).

Since all the results have been derived from the same variational principle, they naturally satisfy some consistency properties fulfilled by exact quantities. For instance, if an observable Q j belongs to the Lie algebra and commutes with H, exact conservation laws express the constancy in time of the expectation value Q j t and of the fluctuation ∆Q j (t). These two properties are ensured by the variational approximations (Sec. VI C). In contrast, a naive evaluation of fluctuations based on the approximation Tr Q j 2 D(t) ≃ Tr Q j 2 D(0) (t) where D(0) (t) evolves according to (60), would provide an unphysical time dependence: The Schrödinger picture for D(0) (t) is variationally suited to the evaluation of expectation values, but not of correlation functions or fluctuations, which involve here the approximate Heisenberg operators Q H j (t, t ′ ). The variational approximation (195) fulfils another consistency property, as discussed in Sec. VI A, namely that any correlation function C jk (t ′ , t ′′ ) depends only on the time difference t ′ -t ′′ when H = K. We have also stressed that the approach unifies static and dynamical properties, generalizing theorems well known for fermion systems [START_REF] Thouless | [END_REF][START_REF] Cloizeaux | Many-Body Physics[END_REF][START_REF] Mermin | [END_REF] : When the thermodynamic equilibrium is stable, i.e., when the trial free energy f {R} is minimum at {R} = {R (0) }, the matrix F of second derivatives is positive, and this implies that the eigenvalues of iCF are real. The latter property ensures dynamical stability: a small deviation of {R} around {R (0) } is never amplified (Sec. IX C).

The approximate equations ( 182) and (186) determine the expectation values Q j t at the initial time t i and for arbitrary t, express these quantities in terms of the expectation values R (0) α (t) = Tr M α D(0) (t), which are the symbols of the elements M α of the Lie algebra. Thus, at first order in the sources, these elements M α are replaced by the c-numbers R (0) α (t) as if the algebra were replaced by a commutative one. Accordingly, the approximate dynamical equations then have a classical Lie-Poisson structure (Sec. X). At second order in the sources, we have shown (Sec. VII) that in the initial state the fluctuations and correlations B of the operators M α have the same form as if the Lie algebra were replaced by a mapped algebra in which the commutators [M α , M β ] are replaced by c-numbers (109), which are their expectation values [START_REF] Blaizot | Quantum Theory of Finite Systems[END_REF]. We have also seen that this amounts to replace the deviations (138) by linear combinations of operators of non-interacting quasi-bosons (or of harmonic oscillators) together with quasi-scalars. This property extends for H = K to two-time correlations. We thus acknowledge two successive modifications of the Lie algebra: For first-order quantities, the operators M α are approximated by scalars, whereas for second-order quantities it is their commutators which are approximated by scalars.

Although formal and technical, the present approach is systematic, flexible, and consistent owing to the optimization of the generating functional. It allows a variational evaluation of correlation functions, which lie beyond the realm of standard mean-field theories. It generalizes some known approaches to arbitrary Lie groups, setting them within a unified framework. It may be applied to various questions of statistical physics and field theory, static or dynamic,

α=

  Tr M α A D/Tr A D, R DA α = Tr M α DA/Tr DA, Z AD = Z DA = Tr AD, which should be expressed in terms of our basic variables {L}, {R}, Z by relying on the group properties.

  α (τ ), again determined by the coupled equations (44) and (45) but with the boundary conditions A (1) α (β) = M α and D (1) α (β) = 0, one can explicitly factorize R (1)

  of a unique effective state D. The first step consists in replacing, in the original Lie structure [M α , M β ] = i Γ γ αβ M γ , the operator M γ on the right side by the c-number R (0) γ = Tr M γ D(0) , the expectation value of M γ in the state D(0) of H . This procedure associates with the original Lie algebra {M} in H a reduced Lie algebra {M} characterized by the simpler commutation relations

  nb = Tr δ D(0) A † n , δR (0) ps = Tr δ D(0) Y p the variations of the symbols of the quasi-boson and quasi-scalar operators A n , A † n and Y p , one finds from the variation of f {R} the form

2 ,

 2 while each quasi-scalar mode brings in the equipartition contribution (2β) -1 .B. Dynamic deviationsLet us turn to the change in the time-dependent expectation values induced by a shift D(0) → D(0) + δ D(0) around the equilibrium state D(0) . The resulting small deviations δR

  α (t) are governed by the set of equations d δR

α

  (t) in Eq.(62); the initial conditions δR

( 0 )

 0 ps ≡ Tr δ D(0) (t)Y p remain constant since { {Y p , h} } = 0 for any h. A classical symplectic structure thus

  ′′ . This independence on t ′′ agrees with the expression (70) ofQ j t . { {f, g} } = -{ {g, f } }; (iii) Jacobi identity: { {{ {f, g} }, k} } + { {{ {k, f } }, g} } + { {{ {g, k} }, kf } } = 0; (iv) Leibniz derivation rule: { {f g, k} } = f { {g, k} } + { {f, k} }g .We consider here functions g{R} of the variables {R}, and define a Poisson structure through

	for any intermediate time t (i) bilinearity; (ii) antisymmetry:	t ′′ ) δR (0) α (t ′′ )	(166)

Acknowledgments A. Derivation of the correlation matrices B and B K 1. Ordinary correlations 2. Kubo correlations B. The variational results and the effective state References Acknowledgments

During the long elaboration of this work, Marcel Vénéroni constantly manifested his deep gratitude to Marie-Thérèse Commault for her efficient and friendly assistance in setting the manuscript into shape. Denis Lacroix from the IPN at Orsay has kindly revised the references.

C. Dynamic and static stability

We now specialize to the case H = K, for which an initial shift {δR (0) } generates a deviation {δR (0) (t)} of {R (0) (t)} around the fixed value {R (0) }. The linearized equations (165) then involve the constant kernel L + C H = C F, and can be solved as δR (0) α (t) = e (t-ti) C F β α δR (0)

Here, as in the backward equations ( 114),(115), the dynamics is governed by the product C F, so that the result (167) can alternatively be written as

Within the variational treatment we recover for small deviations the equivalence between Schrödinger and Heisenberg pictures, already exhibited in Eq.( 70).

If the matrix F is positive, one can change the basis so as to diagonalize iC F according to (132) and use the new variables δR

δR (0) ps (t) = δR (0) ps .

These evolutions merely reflect the motion (154) of the Heisenberg operators {A H , Y H }. Again the modes are decoupled, quasi-bosons (or oscillators) undergo pure oscillations and quasi-scalars are static.

The above sinusoidal form of the dynamics entails Lyapunov stability. This property is defined, according to [START_REF] Mackay | Hamiltonian Dynamical Systems[END_REF], as follows: "The equilibrium point x 0 is said to be Lyapunov stable if given any neighborhood U of x 0 , there is a subneighborhood V of x 0 such that if x lies in V then its orbit remains in U forever." In other words, the trajectories tend uniformly to {R (0) } as their initial point {R (0) (t i )} tends to {R (0) }. Hence, the Lyapunov stability of the linearized motion is ensured if all eigenvalues of F are positive, that is, if the approximate free energy associated with D(0) is a local minimum of f {R}. Such a property is well known for fermions, both at zero [START_REF] Thouless | [END_REF] and at finite temperature [START_REF] Cloizeaux | Many-Body Physics[END_REF][START_REF] Mermin | [END_REF], in which case the minimization of the Hartree-Fock (free) energy entails the reality of the RPA modes. We find here, in a general variational context for any trial Lie group, that the static stability of the approximate "state" D (0) implies the dynamic stability of motions δR (0) (t) around it.

The matrix F may have vanishing eigenvalues. This occurs, for instance, if a continuous invariance is broken by the approximation D(0) for the equilibrium state; in this case f {R} is minimum for a continuous set of solutions {R (0) }. Since F is then not invertible, iF 1/2 CF 1/2 is no longer defined and it is not ascertained that the matrix iCF is diagonalizable (some right eigenvectors may be missing). If iCF is diagonalizable, all its vanishing eigenvalues yield a constant contribution to the set δR (0) (t) so that Lyapunov stability is still ensured. However, if it is not, contributions behaving as powers of t come out, so that the dynamics is unstable (for a more detailed discussion, see Appendix B of [START_REF] Balian | [END_REF]). Such a behavior may be associated with Goldstone modes. For instance, if the original problem is translationally invariant and if this invariance is broken by a localized solution D(0) , a small shift in the initial conditions produces an instability, characterized by a boost at some constant velocity.

X. CLASSICAL STRUCTURE OF LARGE AMPLITUDE MOTION

We have written in Sec. IV B several equivalent dynamical equations for the time-dependent expectation values {R (0) (t)} of the operators {M}. We now analyse the structure of these equations, which will turn out to have a classical form for any Lie group.

A. Poisson structure

Let us consider {R} as a set of classical dynamical variables, and let us show that the tensor C αβ {R} ≡ Γ γ αβ R γ can be regarded as the generator for this set of a Poisson structure issued from the Lie algebra (21). We first recall the definition of Poisson structures [40,[START_REF] Weinstein | Fluids and Plasmas : Geometry and Dynamics[END_REF]. Consider some set of dynamical variables {x} parametrizing points on a manifold, and functions f, g, h, ... of these variables. A Poisson structure is a mapping { {f, g} } = k from a pair of functions f, g to a third function k, which obeys the following rules: appears for linearized motions, with ordinary Poisson brackets {δX n , δP n ′ } = δ nn ′ , {δX n , δX n ′ } = {δP n , δP n ′ } = 0 and Hamiltonian (1/2) n Ω n (P n 2 + X n 2 ) . The above canonical classical structure pertains to the dynamics of small amplitude motions of {R (0) (t)} = {R (0) }+ {δR (0) (t)} around equilibrium {R (0) }. We have seen however (Sec. X A) that a more elaborate Poisson structure occurs for large amplitude motions. In this case, the bracket { {R α , R β } } ≡ C αβ {R} depends on the dynamical variables, whereas for linearized motions we had simply to diagonalize the constant matrix C αβ {R (0) } according to (134). Moreover, the motion is no longer harmonic.

Nevertheless, one can re-express the classical dynamics of Sec. X A in terms of canonical variables by means of a non-linear change of variables in the space {R}. Indeed, a "splitting theorem" [40] states than an arbitrary Poisson structure can locally be split into symplectic components and invariant components: There exist independent (nonlinear) functions ξ n , π n , η p of the coordinates {R} such that their brackets (171) reduce to

In general such a transformation does not exist globally but only locally in some neighborhood around each point of the trajectory of the point {R (0) (t)} in the space {R}. The dynamical variables ξ n and π n are canonically conjugate in the elementary sense (their brackets { { , } } reduce to ordinary Poisson brackets). These variables are dynamically coupled and their motion (173) is governed by Hamilton's equations, while the variables η p (the Casimir invariants) are structurally conserved in the flow.

The construction of the set {ξ, π, η} is not simple and does not result from the mere diagonalization of C αβ {R} at each point, contrary to the linearized dynamics. In fact, the reduction of the Lie-Poisson structure to the form (175) has been achieved only in special cases; relevant to the present work are the Vlasov equation for which a rigorous proof has been given [40], the time-dependent Hartree-Fock theory at zero [16,[44][45][46][47] and finite [START_REF] Balian | [END_REF] temperatures.

C. Stability of equilibrium and of non-linearized motions

In Sec. IX C we have shown that small amplitude motions around thermodynamic equilibrium {R (0) } governed by the linearized equations for {R (0) (t)} (with H = K) are Lyapunov stable if the trial free-energy function f {R} is minimum at {R (0) }. Let us show that this stability property also holds for motions around equilibrium governed by the non-linearized equations of motion (62).

To this aim, we rely on the classical form (172) of this equation written in the Poisson formalism. We note moreover that, owing to the property (174) of the entropy function S{R (0) (t)}, the addition of -β -1 S{R (0) (t)} to the classical Hamiltonian h{R (0) (t)} = k{R (0) (t)} does not affect the dynamics, so that the free-energy function also generates the large amplitude trajectories of {R (0) (t)}:

The time-dependence of f {R (0) (t)} is given, according to (173), by df {R (0) (t)}/dt = { {f {R (0) (t)}, f {R (0) (t)}} } = 0. Hence f {R (0) (t)} remains constant along any trajectory. On the other hand, if {R (0) } is a stable static equilibrium point, f {R} has an isolated minimum equal to f {R (0) } = F . These two conditions ensure Lyapunov stability according to the Lagrange-Dirichlet theorem (see for instance [START_REF] Weinstein | Fluids and Plasmas : Geometry and Dynamics[END_REF]): the trajectory {R (0) (t)} uniformly tends to the point

Thus, for non-linearized as well as for linearized motions of {R (0) (t)}, the present variational approximations preserve the following property of the exact dynamics near an equilibrium state: the static stability, f {R} minimum at {R} = {R (0) }, entails the Lyapunov stability of the dynamics in some neighborhood.

XI. R ÉSUM É OF OUTCOMES AND DIRECTIONS FOR USE

We have dwelt at length on the derivation of the results of a variational approach based on the principles of Sec. II and on the restriction of trial spaces to Lie groups. We summarize below some of the formal outcomes thus obtained; they are sufficient for practical applications to specific many-body problems.

A. General features

The initial step consists in selecting, among the whole set of operators in Hilbert space, a Lie algebra spanned by a set of operators {M} labeled by the index α. (We include in this Lie algebra the unit operator I, denoted as M 0 .) We gave above as seminal example a system of fermions for which the set {M} encompasses the single-particle operators a † µ a ν , the index α denoting the pair (ν, µ). The approach can be applied to many other systems. For fermions with pairing, the Lie algebra {M} includes the operators a µ a ν and a † ν a † µ (with µ > ν); here the index α denotes not only the pair (ν, µ) but also distinguishes between the operators a † µ a ν , a µ a ν and a † ν a † µ . For bosons with possible condensation, one can take as basis {M} for the Lie algebra the operators I, a µ , a † µ , a † µ a ν , a µ a ν (µ ≥ ν) and a † ν a † µ (µ ≥ ν); coherent states arise from the inclusion of the operators a µ and a † µ . If such a bosonic system is translationally invariant, the algebra can be reduced to

where k denotes the momentum of single-particle states). The formalism also applies to other sub-algebras, or to other systems such as mixtures of fermions and bosons, spin systems or quantum fields as in [15], the only restriction being the Lie algebraic structure of the set {M}. In all such cases the algebra {M} is characterized by the structure constants Γ γ αβ entering the commutation relations

As we wish to derive static or dynamic properties of some observables Q j of interest, we have relied on the variational evaluation of a generating functional (Secs. II A and II C), so as to deal simultaneously and consistently with all such properties. This has entailed the introduction of two sets of trial objects, namely, trial "generating operators" A that depend on the observables Q j and their associated (time-dependent) sources, and trial density operators D. Both A and D are elements of the Lie group generated by the chosen Lie algebra {M}, that is, exponentials of elements of this algebra. We focused on expectation values and pair correlations of the observables Q j ; they are found by expansion of the generating functional in powers of the sources. We summarize below the results thus obtained.

We parametrize the trial operators D, which behave as non-normalized density operators, by their normalization Z = Tr D and by the numbers R α = Tr M α D associated with the operators M α (α = 0), where D denotes the normalized operator Z -1 D. The exponential form of the operator D is suited to investigate systems at non-zero temperature; ground state problems are dealt with by taking a zero-temperature limit. We may thus approximately answer questions about an equilibrium (unnormalized) state D = exp(-βK); for a system in canonical equilibrium, K is the Hamiltonian. Other choices of K allow us to deal with non-equilibrium problems, D = exp(-βK) being then the initial state.

An essential tool consists in the representation of the observables Q j , of the operator K entering D = exp(-βK) and of the Hamiltonian H by their symbols q j {R}, k{R} and h{R} within the Lie group (Sec. III B). Defined for any operator Q and an arbitrary element D of the Lie group as

a symbol is a function of the scalar variables R α (α = 0) that parametrize D. The practical implementation of the approach is based on the possibility of evaluating explicitly such symbols q{R}. (In the fermionic and bosonic examples, this is feasible owing to Wick's theorem.) One also needs to express in terms of the variables R α the entropy function S{R} ≡ -Tr D ln D.

B. Static quantities

The generalized free energy F ≡ -β -1 ln Tr exp( -β K) and the thermodynamic quantities issued from it are obtained (Sec. IV A) by looking for the minimum of the trial free-energy function f {R} ≡ k{R} -T S{R} (where k{R} is the symbol of K); the values of R α = R (0) α at this minimum are given by ∂f {R} ∂Rα = 0. The resulting, consistent, approximations are then, for the free energy:

for the entropy:

at zero or non-zero temperature, at equilibrium or off equilibrium.

Appendix A: Derivation of the correlation matrices B and B K

Ordinary correlations

As indicated in Sec. V D, the determination of the correlation matrix B involves the solution at first order of the coupled equations ( 44), (45) for D(τ ), A(τ ) in the range 0 ≤ τ ≤ β, with the boundary conditions D

where we denote for shorthand K{R (0) } ≡ K (0) , have been derived in Sec. IV A. At first order, we have parametrized the combination D(τ ) A(τ ) by {R (1)

The matrix element B αβ was then identified with

β ) Tr D (0) .

A preliminary task is to find the τ -dependence of R

. This is done through the equation

a consequence of ( 44), (45) equivalent to

This equation, since the boundary conditions D

α (τ = 0) = 0 and A

α (τ = β) = M α are imposed at different times τ , cannot fully determine the quantity R (1) αβ (τ ) but it will provide its dependence on τ . Expanding (A4) up to first order involves the kernel

which depends self-consistently on D

α (τ ) and A

α (τ ) through {R (1) (τ )}. Using (A5), then ( 35),( 36) and (56), one obtains

which is readily solved, in terms of the still unknown quantities R

(1)

We now turn to the determination of D

α (τ ) which enters (A2) for τ = β. (We shall not need A

α (τ ) for 0 ≤ τ < β.) The first order contribution to (44) takes the form

Using the boundary condition D

α (0) = 0, we can solve (A8) as

The integrand involves the transform of the element M γ of the Lie algebra by the element e K (0) (τ ′ -τ ) of the group, with K (0) = M β K (0)β . This transform is also an element of the algebra given by the automorphism [START_REF] Mermin | [END_REF] for D = D (0) and λ = (τ -τ ′ )/β as

Inserting (A7) and (A10) into (A9) leads to

We note that the bracket in Eq.(A11) is the derivative with respect to τ ′ of

where (e x -1)/x is defined by continuity as 1 for x = 0. Thus, by integration of (A11), we obtain for τ = β:

We have thus found the expression of D

α (β) to be inserted into (A2) so as to obtain B αβ . The resulting two terms in (A2) involve a correlation in the state D(0) ∝ exp[-βK (0) ] between two operators of the basis {M} of the Lie algebra, M α and M β , or M β and M γ , respectively. Such correlations have been evaluated in Sec. III C (within replacement of D by D (0) ) and are provided by Eq.( 41). Inserting both terms in (A2), we find

and hence

(A14)

Kubo correlations

In order to find a variational approximation for the Kubo correlation

we first replace A(t) by I in the generating functional ln TrA(t i ) D and introduce, instead of the sources ξ j (t) entering the exponent of A, small sources δJ α entering terms -β -1 δJ α M α added to K in the exponent of D = exp(-βK). We now deal with a single exponential operator instead of a product of two, and Kubo correlations are the second-order terms in {δJ} in the expansion of ln Tr D. The variational approximation B K for Kubo correlations is then obtained through the formalism of Sec. IV A within replacement of the image K by K -β -1 δJ α M α . This yields ln Tr exp(-βK + δJ α M α ) ≃ -βf {R} + δJ α R α where {R} = {R (0) + δR} is given by the stationarity condition β ∂f {R (0) + δR}/∂R α = β F αβ δR β = δJ α . As in the evaluation of the thermodynamic coefficients, the second-order terms in {δJ} are finally found as

Here as in the case of B, this approximation for (A15) is obtained by replacing -S by β F = β K -S in the naive expression [START_REF] Mackay | Hamiltonian Dynamical Systems[END_REF] written for D = D(0) .

where the operators M ′ α in the mapped space H are defined through the shift (we drop M 0 ):

These operators obey the same commutation relations as (109), so that the replacement of {M} by {M ′ } does not modify the trace. Hence, we find Inserting in (B1) and identifying with (B4), we recover at first order

At second order, we obtain in the space H the Kubo correlations of the operators {M}:

In order to derive therefrom the ordinary correlations we proceed as in Sec. III C. We note that

where we used the expression (111) of F and the commutation relations (109). Integration over τ then yields

and hence, through a new integration as in (40),

By combining (112) and (B10) we find that the ordinary correlations of the operators {M} are given by

so that we recover here the matrix B of Eq.( 88), now derived as an exact correlation in the mapped space H .