R Balian

Spht

VARIATIONAL APPROACH TO MULTI-TIME CORRELATION FUNCTIONS

Variational approach to multi-time correlation functions

R. Balian, M.

INTRODUCTION AND OUTLINE

Multi-time correlation functions are common tools in several branches of physics. In many-body problems, in equilibrium and non-equilibrium statistical mechanics and field theories, one deals with causal Green's functions or with linear response functions which are closely connected with two-time correlation functions. For weakly interacting systems the calculation of these functions relies on perturbation theory.

When the interactions are strong, alternative approaches are needed. Among them, variational methods have often proved fruitful. The most familiar ones lead to mean-field approximations which provide a simple understanding of many phenomena; nevertheless, they are ill-suited to the evaluation of correlation functions. In the present paper we propose a more elaborate, but still workable, variational approach. Like any other variational treatment, it is less systematic than perturbative expansions. Through its flexibility, it is expected, however, to lend itself to consistent evaluations of multi-time correlation functions for strongly interacting systems. Variational aspects also exist in diagrammatic methods, where a preliminary step often consists in choosing the vacuum by minimisation of the classical action. Our procedure will be different; the quantity that we shall make stationary will be the generating functional itself.

The physical problems that we have in mind are embodied in the following question. At some initial time t 0 , possibly equal to -∞, the state of the system is characterized by a given density operator D (t 0 ) . In statistical mechanics, or in non-zero temperature field theories, it may represent thermal equilibrium as well as a non-equilibrium situation. It may also be the ground state, in which case D (t 0 ) reduces to a projection operator. For atomic or nuclear collisions, it may describe the two incoming fragments boosted towards each other. We wish to determine the time-correlation functions of some set of operators denoted in the Schrödinger picture by Q j . These may be single field operators, or they may be composite. For instance, in thermo-or hydrodynamics, they may represent local density or local current operators. In sects. 4, 5, 7 and 8 we shall be mainly interested in single-particle operators Q j that are bilinear in fermionic creation and annihilation operators. In the Heisenberg picture, the observables Q j of interest become

Q H j (t , t 0 ) = U † (t , t 0 ) Q j U (t , t 0 ) . (1.1)
We have denoted by t the usual running time, and the evolution operator U is defined through dU (t , t 0 ) dt = -i HU (t , t 0 ) , U (t 0 , t 0 ) = 1, (

with h = 1. The multi-time causal functions, which involve the T -product of an arbitrary number of operators Q j , are obtained as usual from a generating functional. We shall find it convenient to define the latter as

ϕ{ξ} ≡ ln TrA (t 0 ) D (t 0 ) , (1.3)
where the operator A(t), a functional of the sources ξ j (t ) , is

A(t) ≡ T exp    i ∞ t dt j ξ j (t ) Q H j (t , t)    . (1.4)
This operator is unitary and its inverse A † would produce the generating functional for the anticausal functions. To lowest orders, the expansion of ϕ in powers of the ξ's yields the single-time expectation values

Q j (t ) ≡ Tr Q H j (t , t 0 ) D (t 0 ) , (1.5)
and the two-time causal functions

C jk (t , t ) ≡ Tr T Q H j (t , t 0 ) Q H k (t , t 0 ) D (t 0 ) -Q j (t ) Q k (t ) , (1.6)
according to

ϕ{ξ} = i ∞ t0 dt j ξ j (t ) Q j (t ) -1 2 ∞ t0 dt dt jk ξ j (t ) ξ k (t ) C jk (t , t ) + ... . (1.7)
Higher order terms generate the cumulants of the three-time, four-time,... causal functions associated with the operators Q j . Of course, equal-time correlation functions can be obtained as special cases. It may be convenient not to normalize D (t 0 ) . A zeroth order term then arises in the expansion (1.7); for instance, if we take D (t 0 ) = exp(-βH) in statistical mechanics, -ϕ{0}/β is the free energy.

Our strategy will be based on the use of a variational expression admitting the generating functional (1.3) as its stationary value. We shall rely on a general method [1,2] that allows for a systematic construction of variational principles specifically suited to the optimization of a desired quantity. A prerequisite is the knowledge of some simple set of equations characterizing this desired quantity, which is here ϕ. It is in order to obtain such an equation for A (t 0 ) that we have introduced a running time t in our definition (1.4).

Indeed, we shall see in sect. 2 that even though it is a complicated operator, A(t) satisfies the simple linear differential equation

d dt A(t) = i[A(t), H] -iA(t) j ξ j (t)Q j . (1.8)
As this equation is supplemented by the boundary condition A(+∞) = 1, it must be solved backward in the time t. For a given set of sources, this yields A (t 0 ) .

The simplicity of the backward differential equation (1.8) for A(t) can be explained by the following remarks. Usually, the Heisenberg operator (1.1) is regarded as a function of the running time t , while t 0 is just a fixed initial reference time at which the Heisenberg and Schrödinger operators coincide. Let us instead, in U (t , t 0 ) and hence in Q H j (t , t 0 ) , regard t as fixed and t 0 as the new running time. In terms of this variable, now denoted as t, the evolution operator U (t , t) satisfies dU (t , t) dt = iU (t , t) H, U (t , t ) = 1, (1.9) to be compared with the usual equation of motion (1.2). Likewise, we can characterize the Heisenberg operator (1.1) as being the solution of the "backward Heisenberg equation" [3]

d dt Q H j (t , t) = i Q H j (t , t) , H . (1.10)
As the boundary condition at t = t is Q H j (t , t ) = Q j , and as t is posterior to t, the time t runs backward in (1.10). This equation (1.10) should be compared to the usual, forward, Heisenberg equation. In case neither Q j nor H depend explicitly on time in the Schrödinger picture, the latter equation, which results from (1.1) and (1.2), reads

d dt Q H j (t , t) = -i Q H j (t , t) , H , (1.11)
with the boundary condition Q H j (t, t) = Q j at the initial time. In the general case, however, the forward Heisenberg equation suffers from two defects. (i) For a time-dependent Hamiltonian H (t ) , the operator which enters this equation is not H (t ) itself, but its expression in the Heisenberg picture, H H (t , t) ≡ U † (t , t) H (t ) U (t , t) .

(1.12) (ii) If, moreover, the observable Q j (t ) depends explicitly on time in the Schrödinger picture, we must include its derivative, ∂Q H j (t , t) ∂t ≡ U † (t , t) dQ j (t ) dt U (t , t) , (1.13) in the right-hand side of the equation. Altogether, the conventional Heisenberg equation in the general case has the form

d dt Q H j (t , t) = -i Q H j (t , t) , H H (t , t) + ∂Q H j (t , t) ∂t , (1.14)
with the boundary condition The existence of two alternative equations of motion, (1.14) and (1.15), is analogous to the situation in the theory of continuous Markov processes. There also, the transition probability can be obtained either as the solution of the forward Kolmogorov equation (better known in physics as the Fokker-Planck equation) or as the solution of the backward Kolmogorov equation [START_REF] Gardiner | Handbook of stochastic methods[END_REF]. Moreover, although these equations are in principle equivalent, the backward Kolmogorov equation possesses simpler mathematical properties [START_REF] Karlin | A first, and a second course in stochastic processes[END_REF] and turns out to be more convenient for some physical applications such as fragmentation [START_REF] Cvitanovic | [END_REF].

Q H j (t, t) = Q(t
Returning to our problem, we note that our operator of interest (1.4) depends on the "observation" times t through the sources ξ j (t ) , and that it is moreover integrated over t . The conventional Heisenberg equation (1.14) is therefore of no use to us. However, apart from being a functional of the sources ξ j (t ) , the operator (1.4) appears as an explicit function of the initial time t, both through the Heisenberg operators Q H j (t , t) and through the lower integration bound. As we shall see in sect. 2, the two terms of (1.8) readily follow from this twofold dependence on t, the first term being a direct consequence of the backward equation (1.15).

Following the general scheme for building variational principles [1,2], we shall treat the equation (1.8) for A(t) as a set of constraints, with which we associate Lagrangian multipliers (sect. 3). We thus find an expression, eq. (3.2), containing two variational objects: A(t), which simulates A(t) itself, and the Lagrangian multiplier D(t), which accounts for the equation of motion for A(t). The stationary value of this expression yields the generating functional (1.3). By restricting the trial classes of A(t) and D(t) so as to make the calculations feasible, one can find variational approximations for this generating functional.

The original expression (1.3) that we want to evaluate depends only on one unknown quantity, the operator A (t 0 ) . Nevertheless, in order to achieve a variational formulation, we have been led to introduce two time-dependent trial operators, A(t) and D(t). While A(t) runs backward in time, D(t) runs forward.

Such a doubling of variables, with two different arrows of time, is of course reminiscent of old approaches to time-dependent problems [7 -9]. However, our present problem is stated more naturally in the Liouville space than in the underlying Hilbert or Fock space. For this reason, our two variational objects, A(t) and D(t), which are dual to each other in the Liouville space, are of different natures: while the former is akin to an observable, we shall see that the latter is akin to a density operator.

An example, drawn from statistical mechanics of interacting fermions, is given in sect. 4. We take exponentials of single-particle (s.p.) operators as trial classes for the operators A(t) and D(t) in Fock Fock space, where H generates both the equations for the states and for the observables, the kernels of the two corresponding approximate equations differ. Moreover, each of these two equations appears variationally suited to the evaluation of a specific quantity. In sect. 5 we extend to response functions the results of sect.

4 on causal functions (see eq. ( 5.3)). The main objective there is to verify that our variational treatment satisfies a number of consistency requirements. We previously checked in sect. 4 that it satisfies some conservation laws which are violated by standard mean-field approaches.

In practical cases, the exact initial state D (t 0 ) is intractable and the above procedure must be supplemented with some approximation for D (t 0 ) . For instance, in equilibrium statistical mechanics, the canonical distribution cannot be handled directly for interacting particles. More generally, whether the system is in equilibrium or not, the initial state usually involves intricate correlations. However, if we write it as

D (t 0 ) ≡ exp(-β H), (1.16)
the operator β H is in general a single-particle or a two particle operator, much simpler than D (t 0 ) itself.

Indeed, D (t 0 ) most often describes a situation of, at least partial or local, equilibrium. It is thus obtained by maximizing the von Neumann entropy -Tr D ln D subject to constraints imposed by macroscopically controlled data. This procedure leads to the form (1.16), where β H is a linear combination of the simple observables associated with these data.

There exist several standard approaches [10], perturbative or variational, which allow one to deal approximately with states of the type (1.16). One class of methods assumes that D (t 0 ) can be generated dynamically in the form

D (t 0 ) = U (t 0 , t 1 ) D (t 1 ) U † (t 0 , t 1 ) , (1.17)
by starting from a simple, uncorrelated, state D (t 1 ) at a time t 1 which will eventually tend to -∞. The correlations existing in D (t 0 ) are thus built up by the evolution from t 1 to t 0 and they can be evaluated approximately from (1.17). However, such an approach is suitable only for infinite systems and it is not fully under control. It requires that dissipation takes place between t 1 and t 0 , and that the data characterizing D (t 0 ) are (exact or approximate) constants of the motion, which are introduced through D (t 1 ) . It might be used in the present context by taking the uncorrelated state D (t 1 ) as our initial condition in sect. 3, and letting t 0 ≡ t 1 tend to -∞. Another class of methods by-passes the difficulties associated with the approach to equilibrium or to quasi-equilibrium. One starts directly from the explicit expression (1.16), and looks for some manageable approximation for it.

We shall adopt an approach of the latter type. However, we shall not replace D (t 0 ) by an a priori approximation in our variational treatment of sect. 3. Instead, we shall again determine the approximation for D (t 0 ) variationally by optimizing the generating functional (1.3) with respect to both A (t 0 ) and D (t 0 ) .

The equations that D (t 0 ) should satisfy are then treated as constraints. We are thus led to merge into a single variational principle, eq. (6.5), the two sets of equations (1.8) and (6.3) which respectively characterize the two operators A (t 0 ) and D (t 0 ) entering the generating functional (1.3). We achieve this task in sect.

6. Here again, following the general procedure [1,2], we introduce Lagrangian multipliers to enforce the equations characterizing A(t) and D(t). We thus obtain a variational expression that has by construction the generating functional (1.3) as its stationary value, and where the objects A and D play symmetric roles.

This unifying variational principle is again exemplified with systems of interacting fermions in sect. 7,

where we choose once more exponentials of s.p. operators for the various trial spaces. Here not only the TDHF and the backward RPA equations, but also the static HF and RPA equations come out naturally from the variational evaluation of the two-time functions, given by eq. (7.18). Finally, in sect. 8 we apply these results to systems in equilibrium. It is shown that the variational treatment of sects. 6 and 7, where the state D (t 0 ) and the operator A (t 0 ) appear on the same footing, provides approximations for the two-time functions which depend, as they should, only on the time difference and not on the initial time t 0 . Moreover these approximations satisfy the fluctuation-dissipation relation, another consistency property of this fully variational approach.

The general method that we propose is contained in sects. 2, 3 and 6, which can be read independently of the rest. All the other sections constitute an application to a problem of many-body physics, worked out both for its own sake and for illustrative purposes. We have omitted some intermediate steps in derivations that were close to those already published in refs. [2,3,11], where details and a larger bibliography can be found. Nevertheless, the present paper is self-contained, and it subsumes several partial results established in these references in a less direct way.

A DIFFERENTIAL EQUATION WHICH GENERATES T-PRODUCTS

As already indicated in the introduction, our first step in the evaluation of the generating functional

(
A(t) ≡ T exp   i ∞ t dt j ξ j (t ) Q H j (t , t)   = 1 + i ∞ t dt j ξ j (t ) Q H j (t , t) - t >t >t dt dt jk ξ k (t ) Q H k (t , t) ξ j (t ) Q H j (t ,
A(t) = i[A(t), H] -iA(t) j ξ j (t)Q j , (2.2)
with A(+∞) = 1.

Note that the entire dependence of A(t) upon the sources is generated solely by the last term of (2.2).

This simple feature can be traced back to our definition (1.4) or (2.1) of A(t), where we have introduced the reference time of the Heisenberg picture as a lower bound of integration at which the sources are turned on.

Note also that the time-ordering in the definition of A(t) is reflected in the ordering of the operators A(t) and Q j in the last term of (2.2).

The same idea can easily be extended to the generation of anticausal functions by means of the operator

A † (t) ≡ T exp    -i ∞ t dt j ξ j (t ) Q H j (t , t)    , (2.3)
where T denotes the anti T -product. Provided that ξ j Q j is Hermitean, A † (t) is both the inverse and the Hermitean conjugate of A(t). It satisfies the differential equation

d dt A † (t) = i A † (t), H + i j ξ j (t)Q j A † (t), (2.4)
which differs from (2.2) through the ordering of the operators in the last term. By combining (2.2) and (2.4), we obtain the Heisenberg-like equations

d dt A † A = i[A † A, H - j ξ j Q j ], (2.5)

VARIATIONAL PRINCIPLE FOR THE DYNAMICS OF THE GENERATING FUNC-TIONAL

We are seeking an expression which reduces at its stationary point to e ϕ = Tr A (t 0 ) D (t 0 ) , where A (t 0 ) is given by the solution of the differential equation (2.2). Since A(t) is complicated, we shall eventually approximate it by a more tractable expression, A(t), to be determined variationally. In this subsection, we assume that D (t 0 ) is simple enough so that Tr A (t 0 ) D (t 0 ) can be explicitly calculated in the trial class for A.

There exists a systematic procedure for constructing variational principles, and hence variational approximations, especially designed to evaluate some quantity of interest [1,2]. Let us briefly recall this method.

We want to determine the value f 0 ≡ f {x α 0 } taken by a given function f {x α } of a set of variables x α at some point {x α 0 } ; this point is determined as the solution of a given set of equations g β {x α } = 0. In spite of the fact that the equations g β {x α } = 0 fully determine {x α 0 } , we use a Lagrangian multiplier method.

Namely, we associate with them a set of multipliers λ β and introduce the functional

Φ {x α , λ β } ≡ f {x α } - β λ β g β {x α } (3.1)
of the two sets of variables x α and λ β . For unrestricted variations of the point {x α , λ β } , the stationary value of Φ is seen to be f 0 . Variational approximations are obtained by restricting the trial spaces for the unknown quantities {x α } and for the multipliers {λ β } to subspaces such that (3.1) can be explicitly handled. The doubling of variables thus allows us to find the desired quantity f 0 through a definite variational scheme.

Detailed explanations about the procedure can be found in [1,2].

For the problem at hand the desired quantity f 0 is Tr D (t 0 ) A (t 0 ) , and the variables x α are the matrix elements of a trial time-dependent operator A(t) at each time t ≥ t 0 . The index α thus stands for both the time and the pair of indices labelling the matrix A(t) in a fixed basis of the Fock space. The point {x α 0 } represents here the actual operator A(t) defined by (2.1). The equations g β {x α } = 0 which determine it are the matrix elements at each time of the differential equation (2.2), and the label β, like α, denotes both time and matrix indices. The set of Lagrangian parameters λ β constitute a time-dependent operator (for t ≥ t 0 ) that we denote as D(t). Notice that the set {x α 0 } describes A(t) at all times t, whereas the desired quantity e ϕ = Tr A (t 0 ) D (t 0 ) involves only part of this set, namely A (t 0 ) . The expression (3.1) now takes the form

Φ{A(t); D(t)} = Tr A (t 0 ) D (t 0 ) + Φ dyn , (3.2a) Φ dyn ≡ Tr ∞ t0 dt D(t)   dA(t) dt -i[A(t), H] + iA(t) j ξ j (t)Q j   , (3.2b)
where Φ dyn accounts for the dynamics of A(t). The functional Φ depends on the two time-dependent trial operators A(t) and D(t) (for t ≥ t 0 ). It is also a functional of the sources ξ j (t). It should be made stationary with respect to variations of A(t) and D(t) subject to the boundary condition

A(+∞) = A(+∞) = 1. (3.3)
We need not assume that D (t 0 ) has unit trace, since multiplying D (t 0 ) by z 0 simply adds the irrelevant constant lnz 0 to ϕ.

The stationarity conditions with respect to D(t), for any t ≥ t 0 , are

Tr δD(t)   dA(t) dt -i[A(t), H] + iA(t) j ξ j (t)Q j   = 0. (3.4)
In order to vary A(t) we perform an integration by parts in (3.2). We obtain, for t > t 0 ,

Tr δA(t)   dD(t) dt + i[H, D(t)] -i j ξ j (t)Q j D(t)   = 0, (3.5)
and, for t = t 0 ,

Tr δA (t 0 ) [D (t 0 ) -D (t 0 )] = 0. (3.6)
For unrestricted variations of D(t), the conditions (3.4) reproduce, as expected from the general expression (3.1), the constraints g β {x} = 0 which are here the equations of motion (2.2); together with the boundary condition (3.3) these equations imply that A (t 0 ) , at the stationary point, equals A (t 0 ) as defined by (1.4).

For unrestricted variations of A(t), the equations (3.5) and (3.6) reduce to

d dt D(t) = -i[H, D(t)] + i j ξ j (t)Q j D(t), (3.7)
D (t 0 ) = D (t 0 ) . (3.8)
If there are no sources, (3.7) is just the Liouville-von Neumann equation; hence the Lagrangian multiplier D(t) coincides, at the stationary point, with the density operator D(t) of the system in the Schrödinger picture. Note that this picture has come out from our variational principle, even though we have adopted the Heisenberg point of view. If however ξ j (t) = 0, (3.7) is the same equation (2.4) as the one for A † , except for the boundary condition (3.8) which replaces A † (∞) = 1. Hence, its solution is

D(t) = A † (t)e -iH(t-t0) A (t 0 ) D (t 0 ) e iH(t-t0) . (3.9)
When the stationarity conditions are satisfied, the desired characteristic functional ϕ{ξ} is related to the stationary value of Φ by

ϕ{ξ} = ln Φ st . (3.10)
If D(t) were to vary freely, the stationarity conditions (3.5) and (3.6) with respect to A(t) would be pointless. For practical applications, however, we have to restrict the trial spaces of D(t) and A(t) so as to be able to calculate (3.2) explicitly. In this case, the best estimate for ϕ is still given by (3.10) when D(t)

and A(t) vary in the restricted spaces. The equations arising from (3.4) then give an approximation for the time-dependence of A(t), and hence for A (t 0 ) . As illustrated by the example of sect. 4, the variations δD(t) usually depend upon D(t), and the approximate equations of motion for A(t) therefore involve D(t). We thus need to determine D(t) within the same approximation scheme, and the equations ( It also results from (3.4) for δD ∝ D that the integrand of (3.2) vanishes, and hence that the approximate stationary value of Φ is

e ϕ = Φ st = Tr A (t 0 ) D (t 0 ) = Tr A(t)D(t), (3.12)
where t is arbitrary. On the other hand, the time-derivatives of A(t) and D(t) obviously belong to the class of variations δA(t) and δD(t), respectively. Writing the equations (3.4) and (3.5) for such variations and adding them, we obtain

d dt Tr D(t)[A(t), H] = Tr    j ξ j (t)Q j d dt (D(t)A(t))    . (3.13)
In practice, an approximate variational scheme for evaluating causal correlation functions for a given initial state D (t 0 ) can thus be developed as follows. After having chosen the relevant operators Q i , we introduce the functional (3.2). We select for A(t) and D(t) variational spaces which make the calculation of (3.2) feasible in terms of the trial parameters. (Of course these spaces should be large enough so that the resulting approximation does not leave out the most important features of the exact solution.) The optimization of the (time-dependent) trial parameters then yields an explicit form for eqs. (3.3) to (3.6).

The one-time, two-time,... correlation functions come out by expanding the result (3.10) in powers of the sources ξ j (t), according to (1.7). For actual calculations, we shall also expand the stationarity conditions (3.3) to (3.6), and solve them iteratively. The complications due to their coupling will thus be overcome.

Moreover, we shall find it convenient, when working out the expansion of (3.10), to take advantage of the stationarity of Φ. In particular, although at its stationarity point the functional (3.2) depends on the sources both explicitly and through A(t) and D(t), only the explicit dependence contributes to the first derivative of the generating functional with respect to the ξ j (t). We have therefore at the stationarity point of (3.2)

δϕ δξ j (t) = ie -ϕ TrD(t)A(t)Q j = i Tr D(t)A(t)Q j Tr D(t)A(t) , (3.14)
where the last expression follows from (3.12). It has been shown in a similar context (see sects. 3.1 and 5.3 of ref. [11]) that considerable simplifications in the expansion of ϕ can arise from equations such as (3.14).

We shall meet with similar simplifications in the next section.

APPLICATION: TWO-TIME CAUSAL FUNCTIONS AND THE BACKWARD RPA

We now illustrate the formalism of sect. 3 by working out an extension of the mean-field approximation for a many-body system of fermions. The Hamiltonian, in second quantized form, is

H = αβ B αβ c † α c β + 1 4 αβγδ αβ | V | γδ c † α c † β c γ c δ , (4.1)
where B includes the kinetic energy and a possible s.p. potential, and where the matrix elements of the two-body potential (which we shall also denote as V 12 ) are antisymmetrized. We shall be interested in oneand two-time correlation functions of some given operators Q j . The latter can be s.p. operators of the form

Q ≡ αβ Q αβ c † α c β , (4.2)
but they can also be more complicated, including, for instance, the Hamiltonian Q 0 ≡ H.

We assume in the present section that D (t 0 ) is a given independent-particle state, characterized by its trace z 0 and the contraction matrix ρ 0 in the s.p. space:

ρ 0αβ ≡ 1 z 0 Tr D (t 0 ) c † β c α . (4.3)
Our purpose is to determine the generating functional (1.3) by means of the variational principle (3.2).

The exact expression (1.4) for A (t 0 ) suggests that we take as trial space for A(t) the set of exponentials of s.p. operators. This set would contain the exact solution A(t) if all the Q k had the form (4.2) and if the particles did not interact. Likewise, the exact solution of (3.7), (3.8) for non-interacting fermions suggests to take for D(t) the same type of trial space, as in any mean-field approximation. However, since for small values of the sources A(t) is expected to remain close to the unit operator and D(t) close to the time-dependent Hartree-Fock (TDHF) state, we parametrize differently A(t) and D(t). Namely, we set

A(t) ≡ exp   -(t) - αβ L αβ (t) c † α c β   , (4.4)
while we characterize D(t), like D (t 0 ) , by its trace z(t) and its contraction matrix

ρ αβ (t) = 1 z(t) Tr D(t) c † β c α . (4.5)
The variational parameters entering (3.2) are therefore (t), L αβ (t), z(t) and ρ αβ (t), with the boundary conditions

(∞) = L αβ (∞) = 0, z (t 0 ) = z 0 , ρ αβ (t 0 ) = ρ 0αβ , (4.6)
arising from (3.3) and (3.8).

With this choice for the trial spaces, we are able to write explicitly the functional (3.2), even for finite values of the sources, by taking advantage of the algebraic properties of exponentials of s.p. operators [3,11].

Indeed, the products A(t)D(t) and D(t)A(t) are also exponentials of s.p. operators, characterized by their trace in Fock space,

w(t) ≡ Tr AD = Tr DA = z exp -+ tr ln 1 -ρ 1 -e -L , (4.7)
and by the contraction matrices

ρ αβ (t) ≡ 1 w(t) Tr DA c † β c α = 1 e L (ρ -1 -1) + 1 , (4.8)
σ αβ (t) ≡ 1 w(t) Tr AD c † β c α = 1 (ρ -1 -1) e L + 1 . (4.9)
We denoted by tr the trace in the s.p. space, and omitted the t-dependence of , L, z, ρ. Wick's theorem then allows us to write out the s.p. reduction of (3.2) in the form

Φ{ (t), L(t); z(t), ρ(t)} = w (t 0 ) + +∞ t0 dt w(t)   - d dt + tr 1 (ρ -1 -1) + e -L de -L dt + iE (σ ) -iE (ρ ) + i j ξ j (t)F j (ρ )   . (4.10)
The quantity E(ρ) is the Hartree-Fock energy associated with H,

E(ρ) ≡ tr Bρ + 1 2 tr 1 tr 2 V 12 ρ 1 ρ 2 , (4.11)
and F j (ρ) is likewise obtained from the operator Q j by means of Wick's theorem in terms of the contractions ρ αβ ; in the case of the operator (4.2), we have F j (ρ) = tr Q j ρ. The last term of (4.10), which arises from the sources, is unsymmetrical with respect to the two sets of contractions (4.8) and (4.9). This lack of symmetry reflects the fact that we are dealing with causal correlation functions. Would we wish to evaluate anticausal correlations, replacing the T -product by the anti-T -product everywhere would replace F j (ρ ) by F j (σ ) in (4.10). In the following we shall need to write explicitly the variations of the functional Φ. This will introduce the matrices W (ρ ) , W (σ ) and Q j (ρ ) defined by

W αβ (ρ ) ≡ ∂E (ρ ) ∂ρ βα , Q jαβ (ρ ) ≡ ∂F j (ρ ) ∂ρ βα ; (4.11 )
these matrices, which are the s.p. reductions of H and Q j , respectively, depend on ρ and L through the definitions (4.8) and (4.9) of ρ and σ , and hence on time. In the special case (4.2), Q j (ρ) is just the given matrix Q j ; in the general case, Q j (ρ) depends on ρ and F j (ρ) differs from tr Q j (ρ)ρ.

The stationarity of Φ with respect to z(t) is expressed by the vanishing of the integrand of (4.10), which will yield (t) through a quadrature. Its stationarity with respect to (t) is equivalent to (3.11), that is, to dw/dt = 0, and it provides z(t). Hence, (3.12) reads here ϕ = ln w(t) (4.12)

for any value of t. By varying ρ and L we find with the above definitions [3]

i de L dt = {W (ρ ) - j ξ j Q j (ρ )}e L -e L W (σ ) , (4.13)
i dρ dt = (1 -ρ){W (ρ ) - j ξ j Q j (ρ )}ρ -ρW (σ ) (1 -ρ), (4.14)
to be solved with the boundary conditions L(∞) = 0, ρ (t 0 ) = ρ 0 . We are interested in the dependence of ϕ upon the sources ξ j (t); it will be given by (4.12) and (4.7), where the variational parameters ( , L, z, ρ)

depend on the ξ s through the stationarity conditions. As consequences of (4.13) and (4.14), ρ and σ satisfy the TDHF-like equations

i dρ dt = [W (ρ ) - j ξ j Q j (ρ ) , ρ ], (4.15)
i dσ dt = [W (σ ) , σ ] , dE (σ ) dt = 0, (4.16)
which, on account of the relation (3.9) between D(t) and A † (t), correspond to (2.5) and (2.6). Finally, eq.

(3.14), together with (4.8), (4.12) and Wick's theorem, reduces to

δϕ δξ j (t) = iF j [ρ (t)] . (4.17)
According to our general strategy, we now expand this whole set of equations in powers of the sources, in order to determine the successive terms of (1.7). To lowest order, when there are no sources, we have at any time L (0) = 0, and hence ρ (0) = ρ (0) = σ (0) . The equation of motion for ρ (0) is

i dρ (0) dt = W ρ (0) , ρ (0) , ρ (0) (t 0 ) = ρ 0 , (4.18)
and (4.17) provides

Q j (t) = F j ρ (0) (t) . (4.19)
Thus, the best variational answer for the expectation value of Q j at the time t is given by the TDHF equation (supplemented by Wick's theorem if Q j is not of the s.p. type (4.2)).

To next order, we shall obtain from (4.17) the two-time causal correlations C jk . This will require the first-order correction to ρ , defined as

ρ (t) -ρ (0) (t) ∼ ρ (1) (t) ≡ -i ∞ t0 dt k ρ k (t, t ) ξ k (t ) . (4.20)
Indeed, the identification of (1.7) with (4.17) yields

C jk (t , t ) = -tr Q j (t ) ρ k (t , t ) , (4.21)
where Q j (t ) stands for Q j ρ (0) (t ) . The equation of motion for ρ k results from the expansion of (4.15).

This introduces the RPA kernel R associated at the time t with ρ (0) (t), and defined by [11] W ρ (0) + δρ , ρ (0) + δρ

αβ -W ρ (0) , ρ (0)
i d dt ρ k (t, t ) = R(t)ρ k (t, t ) -i Q k (t), ρ (0) (t) δ (t -t ) . ( 4
L(t) ∼ L (1) (t) ≡ -i ∞ t0 dt j ξ j (t ) L j (t , t) . (4.24)
Inserting the expansion of (4.8) and (4.9) into the equation of motion (4.13) for L, we find to lowest order

i d dt L j (t , t) = W ρ (0) , L j + tr 2 V 12 L j , ρ (0) 2 -iQ j δ (t -t ) = -L j (t , t) R(t) -iQ j (t)δ (t -t ) , (4.25)
an equation decoupled from (4.23). In the right-hand side, we have recognized the dual of the RPA kernel (4.22). Since L j (t , ∞) = 0, we have L j (t , t) = 0 for t > t and, hence, L j (t , t -0) = Q j . Thus, L j is given for t < t by the solution of

i d dt L j (t , t) = -L j R, L j (t , t ) = Q j (t ) . (4.26)
We note that the same kernel R occurs in (4.23) and (4.25). This gives rise to the equation

d dt tr L j (t , t) ρ k (t, t ) = -tr Q j (t ) ρ k (t , t ) δ (t -t ) -tr [L j (t , t ) , Q k (t )] ρ (0) (t ) δ (t -t ) , (4.27)
which expresses that tr L j ρ k does not depend on t except for jumps at t = t and t = t . Moreover, according to (4.21), the first term in the right side of (4.27) is simply C jk (t , t ) δ (t -t ) . We can therefore obtain this two-time causal function by integrating (4.27) from t = +∞, where L j vanishes, down to t 0 , where ρ k is easily deduced from the boundary condition ρ (t 0 ) = ρ 0 . Indeed, from (4.8), (4.20) and (4.24), we obtain

ρ k (t 0 , t ) ∼ -ρ (t 0 ) L k (t , t 0 ) [1 -ρ (t 0 )] = -ρ 0 L k (t , t 0 ) (1 -ρ 0 ) . (4.28)
For t > t , the last term of (4.27) is irrelevant since L j (t , t ) = 0, and the integration over t is trivial. We thus find, using (4.28),

C jk (t , t ) = -tr L j (t , t 0 ) ρ k (t 0 , t ) = tr L k (t , t 0 ) (1 -ρ 0 ) L j (t , t 0 ) ρ 0 . ( 4
L kαβ (t , t 0 ) c † α c β L jγδ (t , t 0 ) c † γ c δ . (4.30)
The equation (4.29) encompasses some partial results that were obtained previously [11] for t = t . A convenient formula was also derived long ago [12] in the special case when ρ 2 0 = ρ 0 , j = k and t = t .

Numerical calculations [13 -15] relying on it have shown the relevance of the present ideas for the evaluation of fluctuations in heavy-ion collisions.

(t) is tr Q j ρ (0) (t),
and not tr L j (t, t 0 ) ρ 0 , as might have been erroneously inferred from (4.30) and from the exact expression

Q j (t) = Tr D 0 Q H j (t, t 0 ) .
In the TDHF equation, the exact Hamiltonian H is replaced by the effective s.p.

Hamiltonian W depending on time through ρ (0) . In the backward RPA-like eq. ( 4.25), the exact Hamiltonian H entering the Heisenberg equation (1.10) is replaced by the kernel R, which again depends on time through ρ (0) but differs from W when the particles interact. These different guises of H are the price we pay for the s.p. reduction and for the decoupling of the equations.

We can also introduce H, not as the generator of time-translations, but as an operator Q 0 ≡ H on the same footing as the other Q j 's. We then identify the expectation value F 0 (ρ) with the HF energy E(ρ), and the s.p. reduction Q 0 (ρ) of Q 0 in the Schrödinger picture with W (ρ). It is easy to check from the equation of motion (4.18) of ρ (0) that W (t) ≡ W ρ (0) (t) evolves according to

i d dt W (t) = -W R, (4.31)
for a time-independent Hamiltonian H. On the other hand, the s.p. reduction L 0 (t , t) of the Heisenberg operator Q H 0 (t , t) is defined by the backward RPA equation (4.26), and comparison with (4.31) shows that L 0 (t , t) = θ (t -t) W (t). Hence, the s.p. reduction W of the Hamiltonian H is the same in the Schrödinger as in the Heisenberg picture. This property reflects the identity of the exact H in the two pictures when it is time-independent (eq. (1.12)). It holds in spite of the time-dependence of W and provides a further test of consistency for our approximate equations of motion. Furthermore, the replacement of Q j by Q 0 = H in our variational expression for the two-time functions C jk (t , t ) simply replaces L j (t , t 0 ) by W (ρ 0 ) in (4.29). We thus find an approximation for the time correlation function between the Hamiltonian and another operator Q k , which, satisfactorily, does not depend on the time t , even though W ρ (0) (t ) does.

The above property is related to energy conservation. Let us also consider the consistency of our approximation under other symmetry laws, assuming that Q j is conserved, i.e., that it commutes with H.

In such a case, the exact causal function C jk (t , t ) does not depend on t , although it still depends on t .

We thus wish to examine the time-dependence of the matrix L j (t , t 0 ) , defined by (4.26). We shall specialize to a s.p. operator Q j of the form (4.2), and shall prove that in this case the solution of (4.26) is simply

L j (t , t) = Q j . (4.32)
To this aim, we first express how the kernel R, defined in the s.p. Liouville space by (4.22), acts in Fock space. Let us denote by δρ and by Q a pair of arbitrary test matrices with which we saturate R. According to (4.22), we have

(Q | R | δρ) = tr Q δ[W (ρ), ρ], (4.33)
where the r.h.s. involves a change of ρ from ρ (0) to ρ (0) + δρ. On the other hand, if we denote by D the normalized independent-particle density operator corresponding to ρ, we have

tr Q[W (ρ), ρ] = Tr D[Q, H], (4.34)
where Q is the s.p. operator (4.2) associated with the test matrix Q. Altogether, we obtain

(Q | R | δρ) = Tr δD[Q, H]. (4.35)
Applying now (4.35) to Q j , with an arbitrary δρ, we find that Q j R = 0, and thence we justify (4.32). The approximation (4.29) for two-point functions (for t > t ) thus reduces to

C jk (t , t ) = tr L k (t , t 0 ) (1 -ρ 0 ) Q j ρ 0 . (4.36)
It is obviously independent of t , as it should, at least for a s.p. conserved quantity Q j such as momentum or angular momentum.

EXTENSION TO LINEAR RESPONSES

The various two-time functions can be obtained, for the causal functions from the variational expression (4.29), and for the anticausal functions from a similar expression. In particular, if a time-dependent perturbation of the form k ζ k (t)Q k , acting after the preparation time t 0 , is added to the Hamiltonian, the expectation value of Q j at the time t is for small ζ's obtained by expanding (1.5) in the form

Q j (t ) = Q j (t ) 0 + ∞ t0 dt k χ jk (t , t ) ζ k (t ) + ..., (5.1)
where Q j (t ) 0 is the value when there is no perturbation, and where χ jk (t , t ) is the linear response of

Q j to Q k .
As well-known the latter quantity is given by the retarded commutator

χ jk (t , t ) = -i θ (t -t ) Tr Q H j (t , t 0 ) , Q H k (t , t 0 ) D (t 0 ) , (5.2)
a formula which can be derived directly, and also variationally [3]. From (4.29) and the similar formula for anti-T -products, we find for (5.2) the approximation

χ jk (t , t ) = -iθ (t -t ) tr [L j (t , t 0 ) , L k (t , t 0 )] ρ 0 , (5.3)
while Q j (t ) 0 is variationally given by (4.19).

Like the approximation (4.29) for the correlation function (1.6), the approximate response (5.3) can be interpreted as a s.p. reduction of (5.2), where the backward Heisenberg equation (1.10) for the operators Q H j in Fock space is changed into the backward RPA for the matrices L j . We have noted in sect. 4 that different s.p. reductions occur, depending on the quantity of interest, but that they satisfy several consistency properties of the one-and two-time functions. We have just found on variational grounds still another s.p.

expression, (5.3), suited to linear responses. It is then natural to wonder whether (5.3) also satisfies some general consistency requirements which should be obeyed by any reasonable approximation for response functions.

(i) The exact expression (5.2) remains unchanged if the reference time t 0 entering both the state and the operators is shifted forward to any time t satisfying t ≥ t ≥ t 0 , the state D (t 0 ) being replaced by its Schrödinger value at this time t.

We thus wish to perform such a shift in the approximation (5.3). The key point is the following: whereas for Rδρ. We thus obtain

i d dt L k (t , t) , ρ (0) (t) = R(t) L k (t , t) , ρ (0) (t) -i Q k (t), ρ (0) (t) δ (t -t ) , (5.4)
where, as in sect. 4, we wrote Q k (t) for Q k ρ (0) (t) . We now combine (5.4) with the equation (4.25) for L j (t , t) to find for t > t :

i d dt tr L j (t , t) L k (t , t) , ρ (0) (t) = -i tr L j (t , t) Q k (t), ρ (0) (t) δ (t -t ) .
(5.5)

Hence, the approximation (5.3) for the linear response can equivalently be written as

χ jk (t , t ) = -iθ (t -t ) tr [L j (t , t) , L k (t , t)] ρ (0) (t), (5.6)
for any t such that t > t ≥ t 0 . This expression is the s.p. equivalent of property (i): the time t 0 of (5.3) has been replaced in (5.6) by any later time t (smaller than t ), while the initial contraction matrix ρ 0 has been replaced by ρ (0) (t) which follows the TDHF evolution. This confirms that, in our variational s.p.

χ jk (t , t -0) = -i F [Qj ,Q k ] ρ (0) (t ) . (5.7)
On the other hand, we can let t + 0 = t = t -0 in (5.6) and obtain the equal time commutator in the approximate form

χ jk (t , t -0) = -i tr [Q j (t ) , Q k (t )] ρ (0) (t ) . ( 5
Q j (t ) = F j [ρ(t )] = F j ρ (0) (t ) + tr Q j ρ (0) (t ) ρ(1) (t ), (5.9)
where ρ ≡ ρ (0) + ρ( 1) is itself expanded up to first-order in the ζ's (rather than in the sources as in sect. 4).

The matrix ρ(t) evolves according to the TDHF equation associated with the Hamiltonian H, namely

i dρ dt = W , ρ , W ≡ W [ρ(t)] + k ζ k (t)Q k [ρ(t)] , (5.10)
with the initial condition ρ (t 0 ) = ρ 0 .

To check this consistency property, we need to push the time t in (5.6) up to t . The existence for Due to the jump in Lk , the commutator satisfies the differential equation

L k (t , t) of the jump -Q k (t) at t = t ,
i d dt Lk (t , t) , ρ (0) (t) = R Lk , ρ (0) + i Q k , ρ (0) δ (t -t ) , (5.12)
analogous to (5.4); it vanishes for t < t , in particular for t = t 0 . In order to check the property (iii) we insert (5.11) into (5.1) and compare the result with (5.9). It remains to show that the first-order contribution (in the ζ's) of the solution ρ of (5.10) can be identified with

ρ(t ) -ρ (0) (t ) ∼ ρ(1) (t ) = -i ∞ t0 dt k ζ k (t ) Lk (t , t ) , ρ (0) (t ) . (5.13)
Let us therefore expand (5.10). From the definition of R, we find

i dρ (1) (t) dt = Rρ (1) + k ζ k (t)Q k (t), ρ (0) (t) , (5.14)
with ρ(1) (t 0 ) = 0. Comparison of the differential equations (5.12) and (5.14), and of the associated boundary conditions, completes the proof of (5.13) and hence of the property (iii).

(iv) The consistency between the causal functions and the response functions is related to the symmetry of the former in the interchange of the times t , t and of the observable indices j, k. This symmetry was not obvious on the expression (4.21), and we have only derived the expression (4.29) of C jk (t , t ) for t > t .

Let us now turn to the case t > t > t 0 . When bringing t back to t 0 in (4.21), we meet with the jump of L j ρ k at t = t , which is described by the last term of (4.27). We thus have, for t = t -0,

C jk (t , t ) = -tr L j (t , t) ρ k (t, t ) + tr [L j (t , t) , Q k (t)] ρ (0) (t). (5.15)
As a function of t, in the range t > t ≥ t 0 , the first term of (5.15) is constant according to (4.27) while its second term is constant according to (5.5). (The latter is directly related to χ jk since it accounts for the change in the time-ordering of Q H j and Q H k when t becomes larger than t .) Letting t = t 0 in (5.15) and using (4.28), we finally find for t > t

C jk (t , t ) = tr L j (t , t 0 ) ρ 0 L k (t , t 0 ) (1 -ρ 0 ) + tr [L j (t , t 0 ) , L k (t , t 0 )] ρ 0 = tr L j (t , t 0 ) (1 -ρ 0 ) L k (t , t 0 ) ρ 0 . (5.16)
This result complements (4.29), and the expected symmetry is verified.

THE INITIAL STATE PROBLEM: A VARIATIONAL APPROACH

The variational principle of sect. 3 is manageable only if the initial state D (t 0 ) has a sufficiently simple form. In the illustrative application of sects. 4 and 5, we assumed that D (t 0 ) had an independent-particle form; this made the explicit calculation of Tr A (t 0 ) D (t 0 ) feasible. Unfortunately, in realistic problems (which involve the ground state of many-body systems, the Boltzmann-Gibbs distribution in equilibrium statistical mechanics, or the initial state for non-equilibrium situations) not only the evolution but also the exact density operator D (t 0 ) are complicated. As already mentioned in the introduction, a possible approach consists in replacing this state D (t 0 ) by some simple approximation chosen a priori. For instance, if we wish to apply the formalism of sect. 4 to grand canonical equilibrium, we might replace the actual distribution D (t 0 ) = exp(-βH + αN ) by its standard mean-field approximation; the s.p. density matrix ρ 0 entering our formulae would then refer to the static HF density.

This procedure, however, is unsatisfactory. Though the HF density ρ 0 is variational, it optimizes only the thermodynamic potential. It has no reason to provide for D (t 0 ) the variational approximation that is best suited for the evaluation of the generating functional (1.3), even if we remain within the s.p. trial class. More generally we wish to construct, by means of the general method of refs. [1,2], a single variational principle for e ϕ = Tr A (t 0 ) D (t 0 ) in which both A (t 0 ) and D (t 0 ) will be replaced by trial objects A (t 0 ) and D (t 0 ) , to be determined approximately so as to optimize ϕ. This variational principle should encompass the determination of A (t 0 ) by means of the equation of motion (2.2) as above and, simultaneously, the determination of D (t 0 ) by means of another appropriate equation.

As was already noted, the exact state D (t 0 ) need not be normalized and we shall write it as exp(-β H).

Usually, in contrast to D (t 0 ) , H will be a rather simple operator. In ground state problems, H is identical with the Hamiltonian H, and β is infinite. For a grand canonical state, H is H -µN and β is the inverse temperature. For non-equilibrium problems, the logarithm -β H of the initial state D (t 0 ) is some operator which does not commute with H. For atomic, molecular or nuclear collisions, it includes, in addition to the Hamiltonians of the two incoming fragments, a term -v • P where v is the relative velocity and P the relative momentum operator [13 -15].

In order to apply the general scheme recalled in sect. 3, we need a simple constraint forcing the variational state D (t 0 ) to lie as close as possible to exp(-β H). One or other variant of this technique can be used. A first possibility, proposed in ref. [11] for the evaluation of equal-time correlations, consists in imposing on M ≡ -lnD (t 0 ) the constraint

M -β H = 0. ( 6.1)
The variables {x α } of sect. 3 then include, in Fock space, both the matrix elements of A(t) at each time and those of M, while {x α 0 } stands for both A(t) and β H. The constraints g β include the equations (2.2) and (6.1) for A(t) and D (t 0 ) ; the Lagrangian multipliers {λ β } include D(t), as above, plus a matrix B associated with (6.1). The generating functional e ϕ is thus given as the stationary value of

Tr A (t 0 ) e -M + Tr B(M -β H) + Φ dyn , (6.2)
where Φ dyn is the same as (3.2b).

Rather than extending this procedure to the present problem of the evaluation of multi-time functions, we propose here an alternative method which leads to a more symmetric and more convenient variational formulation. In analogy with the representation of A(t) as the solution of the differential equation (2.2), we build up the exponential form exp(-β H) of D (t 0 ) by means of the Bloch equation d du exp(-u H) + Hexp(-u H) = 0, (

where u runs from 0 to β. We introduce therefore a trial operator intended to approximate exp(-u H) and denoted as D(t), where the complex time t ≡ t 0 + i(β -u) runs from t 0 + iβ to t 0 while D(t) runs from 1 to D (t 0 ) (fig. 1). The equation ( 6.3) will be our simple constraint determining D(t) variationally. We shall regard the segment t 0 + iβ, t 0 as the continuation of t 0 , ∞ and denote this segment by t 0 + iβ < t < t 0 .

Such a replacement of temperature variables by imaginary time variables is commonly used in perturbation theory [10]; it will also prove convenient in our variational context. The use of the same notation D(t), that represents for t > t 0 the Lagrangian multiplier associated with the equation for A(t) as in the previous variational expression (3.2), and that also represents for t = t 0 + i(β -u) the solution of (6.3), will be made clearer by eqs. (6.8) and (6.9); it should not be a source of confusion.

To account for the equation (6.3), which, together with the boundary condition

D (t 0 + iβ) = 1, (6.4)
characterizes D(t) for t 0 + iβ < t < t 0 , we introduce a Lagrangian multiplier operator depending on t. For the sake of symmetry, we denote here again this operator by the same notation A(t) as the operator defined for t > t 0 by (2.2). Hence, we now introduce the functional

Ψ = Tr A (t 0 + 0) D (t 0 + i0) -Tr t0 t0+iβ dt A(t) dD(t) dt + i HD(t) + Φ dyn , (6.5)
where Φ dyn is the same as in (3.2b), and where A(t) and D(t) are constrained to satisfy the boundary conditions A(∞) = 1 and D (t 0 + iβ) = 1. The stationary value of Ψ under arbitrary variations of the operators A(t) and D(t), with t lying either on the segment t 0 + iβ < t < t 0 or on the half-line t > t 0 of the complex plane, provides e ϕ . Our notation allows us to match the two integrals of (3.2b) and (6.5) into a single one along the complex L-shaped contour of fig. 1, after an integration by parts. Note, however, that the rôles of the operators A(t) and D(t) are interchanged along the two sections of this contour: the unknown quantities {x α } are D(t) for t 0 + iβ < t ≤ t 0 and A(t) for t ≥ t 0 , while the Lagrangian multipliers {λ β } are A(t) for t 0 + iβ < t < t 0 , and D(t) for t > t 0 . (We might have written iD(t) H, instead of i HD(t), in the integral of (6.5); this alternative choice would not affect the final results of sects. 7 and 8.)

For restricted trial spaces, the stationarity conditions for t > t 0 are the same as (3.4), (3.5) The stationarity of (6.5) with respect to A (t 0 + 0) and to D (t 0 + i0) yields the matching conditions

Tr δA (t 0 + 0) [D (t 0 + 0) -D (t 0 + i0)] = 0, (6.8)
Tr δD (t 0 + i0) [A (t 0 + 0) -A (t 0 + i0)] = 0. (6.9)

The solution of (6.8), (6.9) simply expresses that both D(t) and A(t) are continuous at t 0 . Our choice of notations anticipated this property.

A(t) for t > t 0 , and hence on ξ(t). It is owing to this dependence that we obtain consistently non-trivial approximate correlation functions from a restricted class for A(t) and D(t), such as the independent-particle class of sect. 7. The situation is reminiscent of broken invariances. For instance, the BCS approximation permits a description of superconductivity by means of a simple variational state in which the conservation of the particle number is broken, in contrast to the exact state.

APPLICATION: VARIATIONAL ASPECT OF THE STATIC RPA

We now apply the formalism of sect. 6 to the fermion system of sect.

where the last term Φ dyn is the same integral as in (4.10). We have denoted the Wick expression for Tr D(t)A(t) H by Ē (ρ ) and shall write W (ρ ) for ∂ Ē/∂ρ . In a grand canonical equilibrium, Ē(ρ) reduces to E(ρ) -µ trρ and W (ρ) to W (ρ) -µ.

The stationarity conditions for ρ and L are the same as (4.14) and (4.13) on the half-line t > t 0 ; for t 0 + iβ < t < t 0 , they take the form

i dρ dt = (1 -ρ) W (ρ ) ρ, (7.2)
i de L dt = W (ρ ) e L . (7.3)
The other stationarity conditions are the vanishing of the integrand of (7.1), the constancy of w(t), and the continuity of ρ, z, L and l across t 0 . We shall make use below of the equations

i dρ dt = W (ρ ) , ρ , d Ē (ρ ) dt = 0, (7.4)
dσ dt = 0, (7.5) which hold for t 0 + iβ < t < t 0 and are consequences of (4.8), (4.9), (7.2) and (7.3). We shall also use equation (4.17) for ∂ϕ/∂ξ, which remains valid.

To zeroth order in the source terms, the solution of the equations (7.2) and (7.3), with the boundary conditions ρ (t 0 + iβ) = 1 2 (a consequence of D (t 0 + iβ) = 1) and L (t 0 ) = 0, is

ρ (0) (t 0 + iβ -iu) = 1 e u W0 + 1 , (7.6)
L (0) (t 0 + iβ -iu) = (β -u) W0 , (7.7)
where W0 is defined by the static HF equation

W0 = W (ρ 0 ) , ρ 0 = 1 e β W0 + 1 . (7.8)
The quantity ρ (0) is a constant, equal to ρ 0 for t 0 + iβ ≤ t ≤ t 0 , while ρ (0) = ρ (0) obeys the TDHF equation (4.18) for t > t 0 . Hence, Q j (t ) is approximated variationally by (4.18), (4.19). This approximation consists in using Wick's theorem with contractions given as follows: the s.p. density matrix ρ (0) (t) evolves between t 0 and t according to the TDHF equation associated with H, whereas the initial condition ρ (0) (t 0 ) = ρ 0 is the static HF density (7.8) associated with H.

To next order, we have worked out in sect. 4 the equations for t > t 0 and the results obtained there still hold. We again introduce the RPA kernel R which involves the TDHF s.p. density ρ (0) (t). From R we obtain L j (t , t 0 ) and L k (t , t 0 ) as above by running the backward RPA equation (4.26). However, a new feature appears now at this order since our approximate initial state ρ (t 0 ) is no longer the HF density ρ 0 but involves the shift:

ρ (t 0 ) -ρ 0 ∼ ρ (1) (t 0 ) ≡ -i ∞ t0 dt k ρ k (t 0 , t ) ξ k (t ) . (7.9)
The correction ρ (1) (t 0 ) is introduced by the variational method in the process of optimizing the two-point functions, while the HF density ρ 0 optimizes only the expectation values Q j (t ) . Its existence implies that the quantity ρ k (t 0 , t ) entering the expressions (4.29) or (5.15) for C jk (t , t ) is no longer given by (4.28) but by

ρ k (t 0 , t ) = ρ k (t 0 , t ) -ρ 0 L k (t , t 0 ) (1 -ρ 0 ) . (7.10)
The additional term ρ k (t 0 , t ) should in principle be derived from the expansion of the equations (7.2) and (7.3) that couple ρ(t) and L(t) for t 0 + iβ ≤ t ≤ t 0 . Their solution is complicated by the mixed boundary conditions,

ρ k (t 0 + iβ, t ) = 0, L (t 0 ) = -i k dt ξ k (t ) L k (t , t 0 ) , (7.11)
which refer to both ends of the time interval, and by the fact that L(t) is not small for t ≡ t 0 + i(β -u), 0 ≤ u ≤ β when the sources ξ are small, as is seen from (6.10) or (7.7). However, we shall show that the explicit determination of (7.10) does not require the full solution of these equations (7.2) and (7.3). First, the expansion of (7.4) yields the equation for ρ k (t, t ) along

t 0 + iβ ≤ t ≤ t 0 , i dρ k dt = Rρ k , (7.12)
where R is the RPA kernel (4.22) associated with H, W0 and ρ 0 . (It reduces to R, associated with H, W ρ (0) and ρ (0) , for grand canonical equilibrium.) The solution of (7.12) then provides

ρ k (t 0 , t ) = e -β Rρ k (t 0 + iβ, t ) . (7.13)
On the other hand, since ρ (t 0 + iβ) = 1 2 commutes with L (t 0 + iβ) , we have

ρ k (t 0 + iβ, t ) = σ k (t 0 + iβ, t ) . (7.14)
We now take advantage of (7.5) and write

σ k (t 0 + iβ, t ) = σ k (t 0 , t ) . (7.15)
Finally, comparing the expansions of (4.8) and (4.9) at t = t 0 yields

σ k (t 0 , t ) = ρ k (t 0 , t ) -[L k (t , t 0 ) , ρ 0 ] , (7.16)
and hence we find from eqs. (7.13) to (7.16)

ρ k (t 0 , t ) = - 1 e β R -1 [L k (t , t 0 ) , ρ 0 ] . (7.17)
We shall not consider here the question of the vanishing eigenvalues of R, in spite of the importance of this problem in situations with broken invariance.

We are now in position to write the expression of the two-time function which accounts for the correction (7.9) optimizing the approximate initial state. We denote this new result as Cjk (t , t ) so as to distinguish it from the approximation C jk (t , t ) given by (4.29), (5.16), obtained by assuming the initial state to be uncorrelated and characterized by the s.p. density ρ 0 . We find Cjk (t , t ) by inserting (7.10) into the first line of (4.29) for t > t , or into (5.15) for t > t . We thus obtain our final expression

Cjk (t , t ) = tr L j 1 e β R -1 [L k , ρ 0 ] , t > t , (7.18a)
= tr L j 1 1 -e -β R [L k , ρ 0 ] , t > t , (7.18b)
where L j and L k stand for L j (t , t 0 ) and L k (t , t 0 ) , respectively. This result, also expressed by

Cjk (t , t ) = 1 2 tr L j coth 1 2 β R ∓ 1 [L k , ρ 0 ] , t < > t , (7.18)
encompasses all the previous mean-field approximations that we wrote for correlation functions of interacting fermions. Here, both the evolution and the initial state are determined variationally. In addition to the HF equation (7.8) for ρ 0 and to the backward dynamical RPA used to evaluate L j (t , t 0 ) and L k (t , t 0 ) , the formula (7.18) involves the static RPA kernel R associated with ρ 0 . A new, doubly variational, aspect of the RPA thus emerges naturally from our approach: the RPA occurs both statically and dynamically in our variational approximation for the two-time functions. Remember however that, within the present variational framework, the single-time functions, as well as the thermodynamic functions for grand canonical equilibrium (or the ground state energy), are given consistently by the static and dynamic HF approximations. Moreover, we look here for a stationary value of Ψ, not for a minimum or a maximum. These features contrast with efforts to relate the static RPA to the Rayleigh-Ritz principle [16].

Let us conclude this section with a few additional remarks. Note first that pure state problems correspond to a limit β -→ ∞ in (7.18). On the other hand, for static problems, to evaluate the initial correlations between the observables Q j in the state D (t 0 ) , we let in the above result t = t ± 0 = t 0 ; in this case, L j (t , t 0 ) and L k (t , t 0 ) should simply be replaced in (7.18) by Q j and Q k . For arbitrary t and t the expression (7.18) accounts both for the dynamics at times later than t 0 , through the equations of motion The first part,

R0αβ,γδ ≡ W0αδ δ γβ -δ αδ W0γβ , (7.20)
describes the commutator with W0 , as R0 acts according to

R0 • = W0 , • = -• R0 . (7.21)
Using this definition together with (7.8), we can readily check that the expression (4.29), (5.16) can also be written as

C jk (t , t ) = 1 2 tr L j coth 1 2 β R0 ∓ 1 [L k , ρ 0 ] , t < > t . (7.22)
Comparison with (7.18) shows that the two-time function C jk associated with an independent-particle state D (t 0 ) is recovered by replacing R by R0 in Cjk . Thus, the contribution to Cjk arising from the correlations between particles in the initial state is described by the second part R -R0 of (7.19). To write this second part, it is convenient to introduce [17] two simpler superoperators: (i) an antisymmetric one,

Cαβ,γδ ≡ i (ρ 0αδ δ γβ -δ αδ ρ 0γβ ) , (7.23)
which describes in the Liouville space the commutator with ρ 0 ,

C• = -i [•, ρ 0 ] = -• C , (7.24)
and (ii) a symmetric superoperator In the equation (4.26) for L j (t , t 0 ) , the shift t → t + δt therefore amounts to a change t 0 → t 0 -δt of the initial time, and it produces a change in L j equal to

Vαβ,γδ ≡ ∂ 2 Ē (ρ 0 ) ∂ρ 0βα ∂ρ 0δγ , ( 7
δL j = -iL j Rδt = 1 β L j CGδt, (8.1)
where we used (7.28). The symmetry of G and the antisymmetry of C imply that (8.1) is also equal to δL j = -β -1 GCL j δt. We can thus write the shift in L k as

δL k = - 1 β GCL k δt. ( 8
δC jk = δt 2β L j CG cotg 1 2 CG 0 ∓ i C -cotg 1 2 CG 0 ∓ i CGC L k = 1 2 δt L j CV, cotg 1 2 CG 0 C L k . (8.3)
The explicit form of (8.3),

δC jk = iδt tr 12 V 12 [(1 -ρ 0 ) L j ρ 0 ] 1 [(1 -ρ 0 ) L k ρ 0 ] 2 -[ρ 0 L j (1 -ρ 0 )] 1 [ρ 0 L k (1 -ρ 0 )] 2 , (8.4)
shows that δC jk does not vanish in general. In contrast to the exact two-time function, the approximation (4.29) is affected by the shift δt of the two times t and t with respect to t 0 . Thus, the approach of sects.

3 and 4 entails a spurious dependence on the initial time when the exact equilibrium state is replaced by its standard HF approximation. (Related pathologies occur in TDHF calculations for heavy-ion collisions; in particular, they do not reproduce the spreading of the wave packet and are not compatible with angular momentum analysis.)

Let us show that this defect of the standard mean-field theories is cured by the approach of sects. 6 and 7, which incorporates a consistent variational determination of the approximate state at the time t 0 . Indeed, the same calculation as in (8.3), starting from (7.29) rather than from (7.22), yields

δ Cjk = δt 2β L j |CG(cotg 1 2 CG ∓ i)C -(cotg 1 2 CG ∓ i)CGC|L k = 0, (8.5)
owing to the identity of the two kernels R and R and to the symmetry properties of C and G.

We can apprehend the above consistency of the approximation (7.18), or (7.29), for the two-time function by interpreting this equation in perturbation theory, with W as an unperturbed Hamiltonian. Then G 0 = Ḡ0 describes a particle-hole propagator. According to (7.28), the factor cothβ R/2 in (7.18) or (7.29) describes an iteration of the hole-particle interaction V with the propagator G 0 . Such an iteration is represented by chain (or ring) diagrams along the segment t 0 + iβ, t 0 . Likewise, the solution of the equation of motion (4.26) for L j or L k , which involves the kernel R, is represented by chain diagrams for real times t > t 0 . The identity of R and R allows us to match the two-time sheets, thus resulting in a single type of diagrams over the whole interval t 0 + iβ < t < ∞.

This consistency can also be exhibited algebraically, rather than diagrammatically, as we now demonstrate. Let us rewrite in a simpler way the approximation (7.18) for the causal functions in the case when D (t 0 ) describes a grand canonical equilibrium. It is then convenient to introduce the function of τ ≡ t -t , defined for both τ > 0 and τ < 0 by

Tr Q H j (t , t 0 ) Q H k (t , t 0 ) D (t 0 ) -Q j Q k ≡ Q j (τ )Q k = Q j Q k (-τ ) . (8.6)
(The expectation value Tr Q H j (t , t 0 ) D (t 0 ) ≡ Q j is here time-independent.) This function is a building block for all the two-time functions. Its Fourier transform satisfies the relation The fluctuation-dissipation theorem readily follows by rewriting (8.7) in terms of the expectation values of the commutator and of the anticommutator of Q j (τ ) and Q k . A variational approximation to (8.6) is provided by (7.18b) for τ > 0. We write it out more explicitly by solving formally the equation (4.26) for L j (t , t 0 ) . Since R = R is now independent of time, we have for t > t 0 L j (t , t 0 ) = Q j e -iR(t -t0) . (8.8)

Taking then advantage of the fact that our approximate expression Cjk (t , t ) depends only on τ = t -t , we let t = t 0 , and the insertion of (8.8) into (7.18b) yields for τ > 0

Q j (τ )Q k = tr Q j e -iRτ 1 1 -e -βR [Q k , ρ 0 ] . (8.9)
Likewise, the insertion of (8.8) into (7.19a) yields for τ < 0

Q k Q j (τ ) = tr Q j e -iRτ 1 e βR -1 [Q k , ρ 0 ] . (8.10)
In fact, both (8.9) and (8.10) hold for arbitrary signs of τ, as can be seen by considering anticausal rather than causal functions.

The form (8.9), (8.10) of our approximate results makes it easy to verify that they satisfy some general properties of the exact two-time functions. The fact that (8.9), (8.10) depend only on the time difference t -t is obvious, although it was not in the original formula (7.18). Moreover, by changing τ into -τ and interchanging j and k in (8.10), we should find (8.9); the check is easy, though not completely trivial.

Our approximation also fulfils the fluctuation-dissipation relation, or equivalently the identity (8.7). To verify this, we note that the eigenvalues of R are real provided the HF state is stable [17]. Thus, the Fourier transforms of (8.9) and (8.10) simply involve 2πδ(ω -R) instead of e -iRτ . Hence, the Bose-like factors of (8.9) and (8.10) can be replaced by 1 -e -βω -1 = e βω e βω -1 -1 and by e βω -1 -1 , respectively, which yields (8.7).

Finally, the invariance under translation of the initial time can be shown to hold for our whole set of approximate multi-time functions. For both variational principles (3.2) and (6.5), the shift δt changes all the sources ξ j (t ) into ξ j (t -δt) and hence produces a change in the stationary value of Φ or Ψ equal to δe ϕ = -iδt Tr where we made use of (3.13). In the approximation of sects. we find that Ē [ρ (t 0 )] = Ē [σ (t 0 )] . The equality of tr ρ and tr σ completes the proof that δϕ, still given by (8.12), vanishes. Thus, the approximate multi-time correlation functions do not show any spurious dependence on the initial time t 0 , provided they are consistently evaluated through the approach of sect. 6

where both the operator A (t 0 ) and the state D (t 0 ) are determined variationally.

CONCLUSION

In sects. 4 and 7 we have derived variationally an approximate expression, eq. (7.18), for two-time correlation and response functions in systems of interacting fermions. Our formalism is suited to both equilibrium and non-equilibrium situations; it accounts simultaneously for the initial state and for the dynamics.

We made use of techniques from many-body physics and field theory. One of our approximations consisted in choosing as trial states the same independent-particle class as in the conventional mean-field theories.

The latter theories, however, yield trivial results for correlation functions; in fact, they are not intended for this purpose, since the variational criterion determining their approximate state is the optimization of only the thermodynamic functions.

The reason why we find more elaborate results, in spite of the independent-particle form of the trial states, lies in our use of a new variational principle (sects. 3 and 6), constructed in accordance with a quite general scheme. This principle directly provides the quantity of interest, namely the generating functional for the multi-time correlation functions, as the stationary value of the expression (6.5). This expression contains a trial time-dependent density-like operator, and also its dual, a time-dependent observable-like operator. Within a restricted trial space, the stationarity conditions lead to coupled equations for these dual quantities. An explicit iterative solution is obtained by expanding in powers of the sources, so as to evaluate the one-time, two-time, ... functions. It is through these features that our variational principle partly overcomes the restriction imposed on our trial space.

FIGURE CAPTION

Figure 1: The complex t-plane. The sources ξ j (t ) and the Hamiltonian H refer to real times t > t 0 , while H ≡ -ln D 0 /β refers to times between t 0 + iβ and t 0 . The equations (3.4) and (6.7) for A(t) should be solved backward in time from t = +∞ down to t 0 then to t 0 + iβ, while the equations (6.6) and (3.5) for D(t) run from t 0 + iβ to t 0 then to +∞.

  space. The variational outcome for the expectation values(1.5) is that of the usual time-dependent Hartree-Fock (TDHF) theory. However, for the two-time functions (1.6), the result, eq. (4.29), is more elaborate and involves a "backward RPA equation", (4.26), related to the standard time-dependent random-phase approximation. This backward RPA equation is associated with small deviations around the forward TDHF equation, and its kernel depends on the solution of the latter. Our variational scheme thus exhibits the TDHF equation as the s.p. reduction of the Liouville-von Neumann equation, and the backward RPA as the s.p. reduction of the backward Heisenberg equation(1.10). In contrast to what happens in the complete

2 )

 2 in terms of a single time t, in spite of the complicated time-dependence of the sources: (i) the properties of the backward Heisenberg equation (1.10) that were discussed in the introduction; and (ii) the replacement, in the definition (1.4), of the initial reference time (possibly equal to -∞) by a running time t.

  3.5),(3.6) become relevant since they are now coupled to(3.3) and(3.4). Altogether, for restricted variations of A(t) and D(t), we should solve the set of coupled differential equations (3.4) and (3.5), with the boundary conditions(3.3) and(3.6). Since (3.3) refers to t = +∞ and (3.6) to t = t 0 , the equation for A(t) should be solved backward in time, while the equation for D(t) should be solved forward in time. We have assumed that D (t 0 ) has a tractable form. It is then natural to suppose that it belongs to the trial class for D(t), so that the solution of (3.6) is given by(3.8). We shall return to this point in sect. 6.The equations (3.4) and (3.5) entail some simple consequences when the trial spaces are such that variations δA proportional to A (δA ∝ A) and δD ∝ D are allowed[3]. In this case, by combining (3.4) and (3.5), we find that d dt Tr D(t)A(t) = 0. (3.11)

R

  αβ,δγ δρ γδ . (4.22) This kernel acts as a superoperator in the Liouville space, the vectors δρ αβ of which are identified with matrices in the s.p. space and are labelled by a pair of indices αβ; we have found it convenient to write the products such as (4.22) with a twist of indices, since the scalar product of x αβ and y αβ in the Liouville space is (x | y) ≡ tr xy = x βα y αβ . It then results from (4.15) that

  .23) The equation (4.23) is not sufficient for determining (4.21); indeed, ρ does not satisfy a simple boundary condition, since in its definition (4.8) we have L(∞) = 0 and ρ (t 0 ) = ρ 0 . Nevertheless, let us write the equation of motion for L j , defined by

ρ ( 0 )

 0 (t) evolves according to the unperturbed TDHF equation (4.18) and L k (t , t) according to the backward RPA (4.26), the commutator of these two matrices follows the forward RPA-like equation (4.23). We check this by making use of the explicit form (4.25) for L k R, of the Jacobi identity and of the explicit form(4.22)

( 4 .

 4 26), and for the correlations already existing in the state D (t 0 ) , through R. In order to analyze the latter point, we split R according to (4.22) into two parts, R = R0 + R -R0 .(7.19)

  .25) which acts on an arbitrary vector δρ according to Vδρ = δ W and reduces to αγ | V | δβ when H reduces to the Hamiltonian (4.1) plus, possibly, s.p. operators. We then have R -R0 = i C V.(7.26)

  dτ e iωτ Q j (τ )Q k = e βω +∞ -∞ dτ e iωτ Q k Q j (τ ) .(8.7)

  t)A(t)) = -iδt Tr [D (t 0 ) , A (t 0 )] H,(8.11)

3 and 4

 4 this change had no reason to vanish: in the example of interacting fermions, it is equal toδϕ = iδt {E [σ (t 0 )] -E [ρ (t 0 )]} .(8.12)However, according to the variational determination of the state D (t 0 ) of sects. 6 and 7, the quantitiesE [ρ (t 0 )] and E [σ (t 0 )] are now related to Ē byE [ρ (t 0 )] = Ē [ρ (t 0 )] + µ tr ρ (t 0 ) , (8.13) E [σ (t 0 )] = Ē [σ (t 0 )] + µ tr σ (t 0 ) ,(8.14)since H = H -µN . The equations (7.4),(7.5) imply that Ē (ρ ) and Ē (σ ) do not change between t 0 + iβ and t 0 . From ρ (t 0 + iβ) = σ (t 0 + iβ)

  The expression (4.29) or (4.30) holds for t > t , and the ordering of the matrices L k and L j reflects that of the T -product in the definition (1.6) of C jk . For t > t , we can obtain C jk by interchanging the operators Q j (t ) and Q k (t ) . We shall check in the end of sect. 5 that this result also follows from (4.21) .19) require bringing t 0 to t by means of the forward TDHF equation(4.18). In particular, if Q j is a s.p. operator of the form (4.2), the consistent variational answer for Q j

	Thus, although the Schrödinger and Heisenberg equations of motion are equivalent for the exact density
	operators and observables, this is no longer true for their s.p. reductions, the TDHF and the backward
	RPA equations. The evaluation of the two-time functions (4.29) involves bringing back from t (or t ) to t 0

and (4.27). For anti-T -products we should reverse the ordering of the matrices in (4.29), consistently with the replacement of F j (ρ ) by F j (σ ) in (4.10). Moreover, (4.29) is consistent with an alternative variational treatment suited to the equal-time correlation functions (ref.

[11]

, sect. 5). Another consistency requirement is satisfied: if all the observables Q j are s.p. ones and commute with H, our approximation gives the exact result since the exact operator A(t) then belongs to our trial class. Altogether, our variational treatment involves two different approximate s.p. objects, with which are associated two different approximate equations of motion. On the one hand, ρ (0) (t) appears as the s.p. reduction of the full density operator D(t). It obeys the TDHF equation (4.18), which is the variational s.p. counterpart of the exact Liouville-von Neumann equation. The quantity ρ (0) (t) was shown to be variationally suited to the evaluation of the single-time expectation values (4.19). On the other hand, the Heisenberg operators Q H j (t , t) are simulated here by the s.p. operators αβ L jαβ (t , t) c † α c β . The backward RPA equation (4.26) obeyed by the latter operators arises as the variational approximation for the backward Heisenberg equation (1.10) when one wishes to evaluate the two-time functions (4.30). the time t in L j (t , t) (or L k (t , t)) by means of the backward RPA equation (4.26), while the single-time functions (4

  approximation, TDHF plays the rôle of the Liouville-von Neumann equation while the backward RPA for L j and L k plays the rôle of the backward Heisenberg equation for Q H j and Q H

k . Note that, whereas the expectation value

(4.19)

ought to be written in the Schrödinger picture and the causal correlation function (4.29) in the Heisenberg picture, the linear response can be expressed in either picture by letting t = t or t = t 0 in (5.6). Actually, eq. (5.6) with t = t had already been derived variationally (ref.

[3]

, eq. (7.25)), but without introducing a generating functional.

(ii) If the two times t and t tend towards each other, the exact expression (5.2) is just the expectation value, at the single time t = t + 0, of the commutator -i [Q j , Q k ] . We have already constructed for this quantity two alternative approximations, which we need to compare. On the one hand, if we regard [Q j , Q k ] as a single operator, we are led through (4.19) to express its expectation value in terms of the contraction matrix ρ (0) by means of Wick's theorem, namely

  αβ (t ) and the second to δ αβ , all other contractions lying within Q j or within Q k . Let us then classify the various terms of (5.7) according to the number of contractions which join Q j and Q k . This number is even. The terms containing no such contraction cancel out since we have a commutator. The terms with two such contractions are obviously in one-to-one correspondence with the terms of (5.8). The two expressions (5.7) and (5.8) thus differ only if (5.7) involves four, or six, ... contractions between Q j and Q k . Altogether, the two approaches provide, satisfactorily, the same result for the equal-time commutator whenever one at least of the operators Q j or Q k is a s.p. operator.(iii) The two-time response functions arise from expanding (5.1) in powers of the ζ k 's up to first-order.Alternatively, we can try to evaluate directly Q j (t ) , without carrying out this expansion, by regarding this quantity as a single-time expectation value reached after an evolution governed between the times t 0 and t by the Hamiltonian H ≡ H + ζ k (t)Q k . This very question was already answered in sect. 4 by eq.

	According to its definition (4.11 ) , Q (0)	.8)

j (t ) is obtained from the full operator Q j by contracting all its creation and annihilation operators except for one pair. Hence (5.8) involves two contractions astride Q j and Q k , one equal to ρ (4.19), which is valid even if the Hamiltonian H is time-dependent. To be consistent we thus need to show that (5.1) and (5.3) can equivalently be obtained as the first two terms of the expansion in powers of the ζ's of

  [START_REF] Gardiner | Handbook of stochastic methods[END_REF]. We first take the operator H characterizing the initial state D (t 0 ) = exp(-β H) as arbitrary, and shall specialize to grand canonical equilibrium in sect. 8 by letting H = H -µN . In order to build for D (t 0 ) and A (t 0 ) an approximation that is suited to the evaluation of Tr D (t 0 ) A (t 0 ) , we rely on the variational principle (6.5), and take the same trial class for A(t) and D(t) as in (4.4), (4.5) and(4.6). This will provide us with a new extension of the mean-field approximation. We use the same notations as in sect. 4, in particular the definitions (4.7) and (4.8) of w(t) and ρ (t), but now the times can lie on either part, t 0 + iβ < t < t 0 or t > t 0 , of the integration

	contour of fig. 1. We obtain
	Ψ{l(t), L(t); z(t), ρ(t)} = w (t 0 )
	-	t0 t0+iβ	dt w	1 z	dz dt	+ tr	1 ρ	ρ -1	1 1 -ρ	dρ dt	+ i Ē (ρ ) + Φ dyn ,

  .2)Inserting now (8.1) and (8.2) into the expression (7.22), equivalent to (4.29), for the two-time function and making use of the definitions (7.24) and (7.28), we find the variation of C jk under the time-shift δt :

Although the symmetry introduced between A and D and between the two parts of the contour is convenient, it is solely formal. Note that, when expanding the equations in powers of the sources, A(t) will be close to 1 only for t > t 0 . For t = t 0 + i(β -u), the exact solution of (6.7) is

and it has a finite non-trivial limit when the sources are turned off, in which case A (t 0 ) = 1.

As in sect. 3, it can be convenient to choose trial spaces such that the variations δD(t) ∝ D(t) and δA(t) ∝ A(t) are allowed. The approximation for ϕ resulting from (3.4), (3.5), (6.6) and (6.7) is then given by ϕ ln Tr A(t)D(t), (6.11) for arbitrary t on the complex integration path, since (3.11) holds everywhere. Moreover, we can use again (3.14) to simplify the calculation of one-and two-time functions.

The functional (6.5) is our basic tool for evaluating multi-time correlation functions. It accounts not only for the time evolution, as in sect. 3, but also for our knowledge of the initial state D (t 0 ) . The approximation D (t 0 ) which follows from our scheme is expected to be an improvement on a priori approximations of the same form. Indeed, we construct it variationally by selecting within the trial class that state D (t 0 ) which optimizes the generating functional ϕ, our quantity of interest. For a practical calculation, after having made an ansatz about the quantities A(t) and D(t) as in sect. 3, we need to write the differential equations (3.4), (3.5), (6.6), (6.7) for A(t) and D(t) along the contour t 0 + iβ < t < t 0 , t > t 0 , with the mixed boundary conditions D (t 0 + iβ) = 1, A(∞) = 1. In order to solve these equations, we can again proceed by recursion while expanding the generating functional (6.11) or (3.14) in powers of the sources, so as to find eventually the one-time, two-time, ... functions.

We shall illustrate this procedure by the example of sect. 7. The explicit form of the stationarity conditions will show us that the approximate initial state D (t 0 ) depends on the sources, in contrast to the exact D (t 0 ) . This fact may look surprising since the sources ξ(t) will be turned on at later times. Actually, the approximate equations for A(t) and D(t), derived by restricting the trial spaces and then by optimizing the variational expression (6.5) for ϕ, are coupled. Such a coupling lets D (t 0 ) depend on the values of

The above formal expressions are useful to check the symmetry of (7.18) in the interchange of (j, t ) and (k, t ) . If we introduce the s.p. entropy S = -tr ρ 0 lnρ 0 -tr (1 -ρ 0 ) ln (1 -ρ 0 ) and its associated metric tensor [18] Ḡ0αβ,γδ ≡ -

we obtain [11,17]

where Ḡ is symmetric and C is antisymmetric. By rewriting the causal function (7.18) in the Liouville space as

we exhibit its symmetry.

The second term of (7.22) or (7.29) accounts for the time-ordering in the definitions of C jk or Cjk .

Only this term survives in the evaluation of a commutator. Hence, the linear response has in the present approximation the same form as (5.3). Even though the static HF density (7.8) is not variationally suited to the evaluation of two-time causal functions, it turns out to optimize not only the single-time expectation values but also the retarded linear responses. All the consistency properties of sect. 5 are therefore satisfied by the expression (7.18). We shall moreover see in sect. 8 that (7.18) fulfils the fluctuation-dissipation theorem.

A SPECIAL CASE: MULTI-TIME FUNCTIONS AT EQUILIBRIUM

Until now the initial state D (t 0 ) = exp(-β H) was arbitrary. Let us now specialize to thermodynamic equilibrium, taking H = H -µN . In this case Q j (t) should be time-independent, and the two-time function Cjk (t , t ) should depend only on the time difference t -t . The first property is easily seen from (4.19), since the TDHF density ρ (0) (t) arising from the HF equilibrium state ρ 0 remains constant when W0 , ρ 0 = 0 and W (ρ) = W (ρ) + µ.

However, as we shall see, the approximation (4.29) for the two-time function, C jk (t , t ) , with the HF choice for ρ 0 , is not satisfactory. Indeed let us evaluate the change of (4.29) under a shift δt of both times t and t , keeping t 0 fixed. Since ρ (0) (t) = ρ 0 does not depend on time, the kernel R = R also remains constant.

As expected, the HF and TDHF approximations are recovered for single-time functions, but our expression (7.18) for two-time correlation functions involves in addition two RPA-type kernels, one associated with the dynamics, the other with the initial state. There exist many paths which lead to the RPA; here, both forms of the RPA are obtained mechanically within a consistent independent-particle scheme, and they acquire together a variational status -although they correspond to a stationary point and not to an extremum.

We have demonstrated in sects. 4, 5 and 8 that our approach satisfies several requirements imposed by consistency, invariance or conservation laws, which are violated by standard mean-field approximations.

This supports the idea that a significant improvement is provided by the unified treatment of sects. 3 and 6.

Indeed, we gained in consistency because our variational principle was specifically designed for the evaluation of generating functionals: this linked the optimization of the dynamics with that of the initial state.

Additional support for this variational approach is furnished by numerical calculations based on earlier partial results [12]. In particular, eq. ( 4.29) has been applied to heavy ion collisions in the special case when t = t and ρ 2 0 = ρ 0 . These computations [13][14][15], which led to a better agreement with experiment, have revealed the importance of the corrections introduced by the present method.

Although we have dwelt on an illustrative example from many-body theory, our method is by no means restricted to systems of interacting fermions or to collisions of complex objects. The problem stated at the beginning of the introduction is quite general, and the variational formalism that we set up in sects. 2, 3 and 6 should prove useful in domains other than statistical mechanics. For instance, variational treatments akin to ours have been devised in quantum mechanics. In particular, variational approximate solutions of one-dimensional time-dependent Schrödinger equations have been numerically tested [19]. On the other hand, simple models of quantum mechanics have been investigated as a first step towards a non-perturbative approach to quantum field theory [20]. Other applications have also been worked out for off-equilibrium quantum field theories in a cosmological context [21]. Being systematic and general, the formalism of sects.

3 and 6 opens new prospects in these directions, as it relies on the generating functional, a basic tool in field theory, and as it is tailored to non-zero temperature situations [10].

This formalism can also be translated to deal with classical statistical mechanics. The operators Q j are then replaced by random variables, the traces by integrals over the many-particle phase space, and the backward Heisenberg equation (1.10) by "backward Hamilton equations" written with Poisson brackets, while D (t 0 ) is now the initial density in phase. On may thus imagine applications to topics such as plasma physics, turbulence, or the large scale structure of the Universe (in which case the quantities of interest are the correlations between the positions of the galaxies). An independent-particle approach corresponding to that of sects. Finally, the flexibility of the approach allows for a combination of variational and perturbative approximations (see ref. [20] and sect. 2.4 of ref. [2]), a question that we have not touched upon in this paper.

We wish to thank Hubert Flocard and John Gillespie for their critical readings of our manuscript.


J A De Jesús Osuna-Coutiño

Abiel-Aguilar-González

Miguel Arias-Estrada

PhD Forum: GPU-based Visual Odometry for Autonomous Vehicle Applications

Keywords: Camera pose, embedded systems, GPU

One popular task under several computer vision applications is camera pose estimation under video sequences. In previous work, several camera pose estimations approaches have been developed and several algorithms have been proposed. Unfortunately, most previous formulations iterative behavior and depend on nonlinear optimizations applied over some geometric constraint and this limits the performance under embedded applications. For real-time embedded applications, another approach, more efficient in terms of computational size and processing speed could be reached via hardware-based solutions, for example GPU-based solutions. In this work, we present a GPU-based solution for the camera pose problem. As early formulation we focused our algorithm for an autonomous vehicle application. Preliminary results are encouraging and show the viability of the proposed approach.

INTRODUCTION

In several computer vision applications, camera pose estimation across video sequences appears repeatedly under several contexts: visual servoing [START_REF] Chaumette | Visual servoing[END_REF], Simultaneous Localization and Mapping (SLAM) [START_REF] Mur-Artal | Orb-slam: a versatile and accurate monocular slam system[END_REF], augmented reality [START_REF] Hong | A 27 mw reconfigurable marker-less logarithmic camera pose estimation engine for mobile augmented reality processor[END_REF], etc. In previous work, primary interest was for an accurate/generalized solution for the camera pose problem [START_REF] Hartley | In defense of the eight-point algorithm[END_REF][START_REF] Longuet-Higgins | A computer algorithm for reconstructing a scene from two projections[END_REF]. Unfortunately, these solutions have iterative behavior and requires relatively high computational resources. Then, several tasks have limitations for embedded applications. In this work, we believe that a smart camera that delivers camera pose at video cadence could be an important contribution for current embedded computer vision applications. In general, several computer vision applications such as SLAM, augmented reality, etc., could reach onboard camera pose processing.

CAMERA POSE ESTIMATION BASED ON 2D PIXEL DISPLACEMENTS

As early formulation, we consider an autonomous vehicle scenario. In this case we suppose that all possible camera movements are linear dependent of the x, y pixel displacements whiting the scene. i.e., given the input image divided in four quadrants (as illustrated in Fig. 1) and, consider four different cases: turn left, turn right, forward and backward; for the turn left case: all pixels within all the quadrants must have right displacement under the horizontal axe and similar displacement for the vertical axe, see Fig. 1a. For the turn right case: all pixels for all the quadrants must have left displacement under the horizontal axe and similar displacement for the vertical axe, as shown in Fig. 1b. For the forward case: pixel displacements are to outside the image, see Fig. 1c. Finally, for the backward case: pixel displacements are to the image center, as illustrated in Fig. 1d.

PRELIMINARY RESULTS

As early validation, we implemented our algorithm in a MatLab R2017a code that uses the KITTI dataset [START_REF] Geiger | Vision meets robotics: The kitti dataset[END_REF] as input sequences. First, feature points (Harris corners) are tracked along the video sequence. Then, several 2D pixel displacements (considering four different quadrants, as shown in Fig. 1) are correlated with its corresponding camera movement (taken from the ground truth in the KITTI dataset). Finally, unknown camera movements were estimated by comparing 2D pixel displacements in the input image with the previously stored pixel displacements and its corresponding camera movement. In Fig. 2, preliminary results are shown. In this case, pixel displacements and camera movements for 800 random frames for the first video sequence from the KITTI dataset were stored, then, camera pose for the third and fourth sequences were computed by comparing the pixel displacements for the input images with the pixel displacements and camera movements previously stored.

CONCLUSIONS AND FUTURE WORK

In this work, we have demonstrated that given an autonomous vehicle scenario camera movements are linear dependent with the 2D pixel displacements within a video sequence. As work in progress, we are generalizing our algorithmic approach and implementing it inside the NVIDIA Jetson TX1 board.

Figure 1 :Figure 2 :

 12 Figure 1: Given an autonomous vehicle scenario. Four different cases: turn left, turn right, forward and backward have to fulfill with all possible movements from a single moving camera installed over the vehicle. In such scenario, pixel displacements inside the images (red arrows) have different directions for each case.