text:
# 1 Introduction

## 1 Introduction

The anomalous $`J/\psi `$ suppression observed by the NA50 Collaboration  in Pb–Pb collisions at the CERN SPS can be naturally interpreted as $`J/\psi `$ dissolution by the Quark–Gluon Plasma (QGP). This phenomenon was predicted by Matsui and Satz  more than 15 years ago. The interpretation of the anomalous $`J/\psi `$ suppression is somewhat complicated by the fact that there are also other contributions to $`J/\psi `$ suppression. Nuclear absorption (NA) (sometimes referred to as ”pre-resonance absorption” or ”Gerschel–Hüfner mechanism”)  is probably responsible for most of $`J/\psi `$ suppresion in nuclear collisions induced by lighter ions. Disintegration of $`J/\psi `$ by collisions with secondary hadrons  (also referred to as ”co-mover” interaction) is another possible contribution. In order to make the $`J/\psi `$ suppression by the gas of secondary hadrons sufficiently large, in some approaches  cross-section ($`\sigma _{co}`$) for the process hadron+$`J/\psi D\overline{D}+X`$ has to be of about $`5`$ mb, depending also on the value of the cross-section ($`\sigma _a`$) for nuclear absorption. In the approach based on the Dual Parton Model (DPM)  one needs $`\sigma _{co}=0.6`$ mb, for $`\sigma _a=6.7`$ mb and  $`\sigma _{co}=1`$ mb for $`\sigma _a=4.5`$ mb. It has been argued  that, at energies corresponding to the thermal motion of hadrons, the cross-sections for the disintegration of $`J/\psi `$ by hadrons are in fact an order of magnitude smaller. The issue is not yet definitely clarified.

There is one key difference between the suppression of $`J/\psi `$ by the QGP and by other mechanisms. For $`J/\psi `$ dissolution by QGP a rapid onset of the anomalous suppression is expected or at least possible, whereas for other mechanisms one expects a smooth dependence of $`J/\psi `$ suppression on the impact parameter $`b`$ and on the total transverse energy $`E_T`$ of a nuclear collision. Models by Blaizot and Ollitrault (BO)  and by Kharzeev, Lourenço, Nardi and Satz (KLNS)  do show a rather rapid onset of the dependence of the anomalous $`J/\psi `$ suppression on $`E_T`$ or on $`b`$. Plots of experimental data on $`J/\psi `$ suppression, in particular those with the nuclear absorption ”subtracted”  do show such an abrupt onset.

Although the anomalous $`J/\psi `$ suppression is a rather spectacular Phenomenon, indicating a rapid onset of a new mechanism of $`J/\psi `$ Suppression, there are a few mechanisms that are able to obtain some agreement with the basic features of the data. Apart from the suppression by QGP  these include the disintegration of $`J/\psi `$ by co-movers  and Monte Carlo microscopic models . In order to see why it is possible that rather different models lead to a similar overall $`J/\psi `$ suppression, it would be interesting to see not only the resulting total $`J/\psi `$ suppression but also contributions of individual mechanisms involved in this total suppression. The purpose of the present paper is to propose an alternative way of looking at the data on $`J/\psi `$ suppression, and in particular on the onset of anomalous $`J/\psi `$ suppression. Our approach is based on the additive decomposition of contributions of different mechanisms responsible for the disintegration of $`J/\psi `$’s produced by a hard process in nuclear collisions. The method can be used provided that the time sequence of contributing mechanisms is known or assumed. In some cases there are good reasons to believe that two mechanisms are at work in the same time period. This would concern, for instance, the nuclear absorption and possible depletion of gluon structure functions  during the first part of the nuclear collision. In such a case both mechanisms should be considered as one stage.

In the next section we shall describe the method; in Sect. 3 we shall present a few illustrative examples. In Sect. 4 we shall discuss the relationship of the proposed way of looking at the data with the standard procedure used by the NA50 Collaboration, and comments and conclusions will be presented in Sect. 5.

## 2 Separation of different contributions to $`J/\psi `$ disintegration

We shall consider here for simplicity the case when only two mechanisms of $`J/\psi `$ suppression are present: the nuclear absorption and the dissolution of $`J/\psi `$ by QGP, the latter taken from a point of view very close to that of the BO  and KLNS  models.

We shall use the following notation

$`N_{prod}^{J/\psi }(E_T)=`$ the total number of $`J/\psi `$ produced in collisions at given $`E_T`$,

$`N_{NA}^{J/\psi }(E_T)=`$ the number of $`J/\psi `$’s disintegrated by NA,

$`N_{QGP}^{J/\psi }(E_T)=`$ the number of $`J/\psi `$’s dissolved by QGP,

$`N_{exp}^{J/\psi }(E_T)=`$ the number of surviving $`J/\psi `$’s, measured in experiment

$`N_{other}^{J/\psi }(E_T)=`$ the numbers of $`J/\psi `$’s disintegrated by other mechanisms.

The balance between input and output requires

$$N_{prod}^{J/\psi }(E_T)=N_{exp}^{J/\psi }(E_T)+N_{NA}^{J/\psi }(E_T)+N_{QGP}^{J/\psi }(E_T)+N_{other}^{J/\psi }(E_T).$$
(1)
We suggest that all terms in Eq. (1) be calculated and plotted separately, or in combinations, for a given A-B collision, for all values of $`E_T`$. In this sense our approach tries to give an answer to the question ”Where are all the $`J/\psi `$’s gone?” or corresponds to ”an accountant’s look at the data on $`J/\psi `$ suppression”.

Admittedly, the suggestion we are making is rather trivial, but we do hope that looking at the data in this way can make the results of a particular model more transparent and can help to avoid some inconsistencies.

In specific applications Eq.(1) can be rewritten in such a way that more reliably known terms are put on one side, thus providing constraints for the more model-dependent and less reliably known ones. We can for instance put on one side the best known terms $`N_{prod}^{J/\psi }(E_T)N_{exp}^{J/\psi }(E_T)`$ or even $`N_{prod}^{J/\psi }(E_T)N_{exp}^{J/\psi }(E_T)N_{NA}^{J/\psi }(E_T)`$, and try to describe or fit these expressions by a particular model. For instance in the latter case by $`N_{QGP}^{J/\psi }(E_T)+N_{other}^{J/\psi }(E_T)`$.

The purpose of the present paper is not to attempt a detailed analysis of the data. We wish just to describe the method and to give a few examples as illustrations. For that purpose we shall take nuclei as homogeneous hard spheres with radii $`R_A=1.2A^{1/3}`$fm, we shall neglect the energy fluctuations, including the ”knee” in multiplicity distribution at $`E_T100`$GeV and we shall assume here a strong correlation between $`b`$ and $`E_T`$:

$$E_T(b)=0.325GeVN_w(b),$$
(2)
where $`N_w(b)`$ is the number of interacting (”wounded”) nucleons. With these simplifications the quantities entering Eq. (1) can be expressed in the following way:

$$N_{prod}^{J/\psi }(E_T)=\sigma _{nn}^{J/\psi }N_{coll}(E_T),$$
(3)
where $`\sigma _{nn}^{J/\psi }`$ is the $`J/\psi `$ production cross-section in the average nucleon–nucleon collision and $`N_{coll}(E_T)`$ is the number of nucleon–nucleon collisions at a given value of $`E_T`$. In a rather simplified case we have

$$N_{coll}(E_T)=_0^{R_A}\frac{sds}{\sigma _{nn}}_0^{2\pi }𝑑\theta \rho _A\sigma _{nn}2L_A(s)\rho _A\sigma _{nn}2L_B(b,s,\theta )=$$

$$_0^{R_A}\frac{sds}{\sigma _{nn}}_0^{2\pi }𝑑\theta _{L_A(s)}^{L_A(s)}\rho _A\sigma _{nn}𝑑z_A_{L_B(s,\theta )}^{L_B(b,s,\theta )}\rho _B\sigma _{nn}𝑑z_B,$$
(4)
where

$$L_A(s)=\sqrt{R_A^2s^2},L_B(b,s,\theta )=\sqrt{R_B^2b^2s^2+2bscos(\theta )},$$

when the expression under the square–root in $`L_B`$ is negative, there is no tube-on-tube collision and the contribution vanishes. Nuclear densities are denoted as $`\rho _A,\rho _B`$ (we are working in the approximation of nuclei as hard spheres, therefore $`\rho _A,\rho _B`$ are constants), $`b`$ is the impact parameter, $`s`$ is the distance from the centre of the A-nucleus and $`\theta `$ is the angle between $`\stackrel{}{b}`$ and $`\stackrel{}{s}`$. A nuclear collision is taken as the sum of tube-on-tube collisions, with lengths of tubes $`2L_A`$ and $`2L_B`$, and $`z_A`$,$`z_B`$ specifying the coordinate of nucleons within both colliding tubes: $`z_A`$ varies from $`L_A`$ to $`L_A`$, $`z_B`$ from $`L_B`$ to $`L_B`$, and both $`z_A`$ and $`z_B`$ increase in the direction of motion of the tubes in the c.m. frame of the nucleon–nucleon collision. The nucleon-nucleon non-diffractive cross-section $`\sigma _{nn}`$ is taken as $`30`$mb. Experimental results on the survival probability $`S(E_T)`$ of $`J/\psi `$ enter the expression

$$N_{exp}^{J/\psi }(E_T)=S(E_T)N_{prod}^{J/\psi }(E_T).$$
(5)
Nuclear absorption of $`J/\psi `$ is given by the expression

$$N_{NA}^{J/\psi }(E_T)=\sigma _{nn}^{J/\psi }_0^{R_A}\frac{sds}{\sigma _{nn}}_0^{2\pi }𝑑\theta _{L_A(s)}^{L_A(s)}\rho _A\sigma _{nn}𝑑z_A_{L_B(s,\theta )}^{L_B(b,s,\theta )}\rho _B\sigma _{nn}𝑑z_B$$

$$(1e^{\rho _A\sigma _a[z_A+L_A(s)]}e^{\rho _B\sigma _a[z_B+L_B(b,s,\theta )]}),$$
(6)
where $`\sigma _a`$ is the cross-section describing the absorption of $`J/\psi `$ by nucleons.

The integral in Eq. (6) can be decomposed into two separate integrals, the former giving the total number of $`J/\psi `$’s produced and the latter being equal to the number of $`J/\psi `$’s that survive the nuclear absorption. The term $`N_{QGP}^{J/\psi }(E_T)`$ gives the number of $`J/\psi `$’s destroyed by the QGP; it is given as

$$N_{QGP}^{J/\psi }(E_T)=\sigma _{nn}^{J/\psi }_0^{R_A}\frac{sds}{\sigma _{nn}}_0^{2\pi }𝑑\theta _{L_A(s)}^{L_A(s)}\rho _A\sigma _{nn}𝑑z_A_{L_B(s,\theta )}^{L_B(b,s,\theta )}\rho _B\sigma _{nn}𝑑z_B$$

$$e^{\rho _A\sigma _a[z_A+L_A(s)]}e^{\rho _B\sigma _a[z_B+L_B(b,s,\theta )]}\mathrm{\Theta }(\kappa \kappa _{crit}).$$
(7)
Following KLNS , see also Ref., we have introduced

$$\kappa =\frac{\rho _A\sigma _{nn}2L_A(s).\rho _B\sigma _{nn}2L_B(b,s,\theta )}{\rho _A\sigma _{nn}2L_A(s)+\rho _B\sigma _{nn}2L_B(b,s,\theta )}.$$
(8)
The parameter $`\kappa _{crit}`$ specifies the onset of the QGP formation. This parameter has to be determined by the data in such a way that

$$N_{prod}^{J/\psi }(E_T)=N_{exp}^{J/\psi }(E_T)+N_{NA}^{J/\psi }(E_T)+N_{QGP}^{J/\psi }(E_T),$$
(9)
where we have assumed that the contribution $`N_{other}^{J/\psi }`$ can be neglected. Equation (9) simply says that out of all $`J/\psi `$’s produced some survived, others were disintegrated by the nuclear absorption and others by QGP.

The multiplicative factor $`\sigma _{nn}^{J/\psi }`$ is present in all of the terms in Eq. (9) and in what follows we shall leave it out.

The time ordering of the different processes enters Eq. (7). It is most likely that the QGP is formed after the colliding nuclei have passed through one another and the QGP can thus destroy only those $`J/\psi `$’s that survived the nuclear absorption.

## 3 Illustrative examples

Figures showing the experimental data in the way they were used by NA50 Collaboration  are well known. In these figures they plot the survival probability of $`J/\psi `$ as

$$S(E_T)=\frac{J/\psi _{measured}}{J/\psi _{produced}},$$
(10)
inferred from direct data on $`J/\psi `$ over Drell–Yan pair production. The coefficient $`S(E_T)`$ contains in a multiplicative way probabilities that $`J/\psi `$ has survived nuclear absorption, QGP dissolution and other possible mechanisms of $`J/\psi `$ disintegration. In the way of looking at data proposed here we plot separately the surviving $`J/\psi `$’s (as given by experiment), those disintegrated by nuclear absorption, and those dissolved by the QGP and possibly also by other mechanisms.

In Fig. 1 we plot as an illustration

the terms $`N_{prod}^{J/\psi }`$, $`N_{exp}^{J/\psi }`$ (survivors), and $`N_{NA}^{J/\psi }`$, leaving out the common factor $`\sigma _{nn}^{J/\psi }`$ in Eqs. (3)–(7). The anomaly is not easily visible, since it appears as a rapidly growing difference $`N_{prod}^{J/\psi }N_{exp}^{J/\psi }`$.

In order to make the anomaly visible, we plot in Fig. 2 the expression $`\mathrm{\Delta }N`$ = $`N_{prod}^{J/\psi }N_{exp}^{J/\psi }N_{NA}^{J/\psi }`$ with a nuclear absorption cross-section $`\sigma _a=0.7fm^2`$. The shape of $`\mathrm{\Delta }N`$ indicates a presence of a $`J/\psi `$ dissolving mechanism with threshold between $`E_T=30`$GeV and $`E_T=40`$GeV. Note that values of $`\mathrm{\Delta }N`$ are slightly negative for the lowest three $`E_T`$ points, which is connected with the fact that in the standard presentation of data calculations based on the NA mechanism are below the lowest $`E_T`$ data.

In Fig. 3 we show the opening of the space for

other contributions, most probably for $`J/\psi `$ suppression by a hadron gas, when $`\sigma _a`$ in nuclear absorption becomes smaller.

The true consistency check of the phenomenological description of data is provided by Eq. (1) with $`N_{exp}^{J/\psi }`$ given by Eq. (5) and $`S(E_T)`$ taken from experimental data on $`J/\psi `$ survival. In the illustrative case discussed above we have $`N_{other}^{J/\psi }=0`$ and Eq. (1) reduces to

$$N_{prod}^{J/\psi }(E_T)=\mathrm{\Sigma }(E_T)N_{exp}^{J/\psi }(E_T)+N_{NA}^{J/\psi }(E_T)+N_{QGP}^{J/\psi }(E_T).$$
(11)
In Fig. 4 we compare $`N_{prod}^{J/\psi }(E_T)`$ and $`\mathrm{\Sigma }(E_T)`$. The

agreement is quite reasonable, in view of simplifying assumptions that we made. Agreement with data is worse for the points with the highest values of $`E_T`$ where a large part of $`E_T`$ is given by fluctuations  and our simplified model is not applicable. Note that the agreement visible in Fig. 4 is non-trivial since $`N_{exp}^{J/\psi }(E_T)=S(E_T)N_{prod}^{J/\psi }`$ is given by experimental data on $`S(E_T)`$.

The information contained in Fig. 1 and Fig. 4 can be compressed into

a single Fig. 5, where we plot $`N_{NA}^{J/\psi }(E_T)`$; $`N_{NA}^{J/\psi }(E_T)`$\+ $`N_{QGP}^{J/\psi }(E_T)`$ and $`N_{NA}^{J/\psi }(E_T)`$\+ $`N_{QGP}^{J/\psi }(E_T)`$\+ $`N_{exp}^{J/\psi }(E_T)`$; and for comparison with the last expression also $`N_{prod}^{J/\psi }(E_T)`$.

## 4 Common points of the new way of looking at data and the standard procedure

We shall first show that the way of plotting experimental data and results of phenomenological calculations is in some aspects equivalent to the standard way introduced and used by the NA50 Collaboration.

To make the argument transparent, suppose that a certain number $`N_0`$ of $`J/\psi `$’s is passing successively through three obstacles (media) and can be disintegrated by each of them. The probability that it passes (survives) through the first obstacle is denoted as $`P_1`$. The number of $`J/\psi `$’s passing through is $`N_0P_1`$, the number of those that fail to pass is $`N_{F1}=N_0(1P_1)`$. In the same way the number of $`J/\psi `$’s passing through the second obstacle is $`N_0P_1P_2`$ and the number of those that passed the first obstacle, but failed to pass the second one is $`N_{F2}=N_0P_1(1P_2)`$. Finally, the number of $`J/\psi `$’s passing throught the third obstacle is $`N_0P_1P_2P_3`$ and the number of those that failed to pass the third obstacle is $`N_{F3}=N_0P_1P_2(1P_3)`$. The obvious identity

$$N_0=N_0(1P_1)+N_0P_1(1P_2)+N_0P_1P_2(1P_3)+P_1P_2P_3N_0$$
(12)
can be rewritten as

$$N_0=N_{F1}+N_{F2}+N_{F3}+SN_0.$$
(13)
where

$$S=P_1P_2P_3=\frac{N_0N_{F1}N_{F_2}N_{F3}}{N_0}$$
(14)
gives the overall survival probability. When comparing the results of the calculations with the data in the usual way the comparison is between experimentally measured values of S and the probability $`P_1P_2P_3`$ calculated from phenomenological models. The alternative way used here consists in proceeding according to Eq. (12), plotting individual terms $`N_{F1},N_{F2},N_{F3},SN_0`$, with $`S`$ taken from experiment and verifying the validity of Eq.(12).

In a real situation, each $`J/\psi `$ is born at a certain value of the parameters $`s,\theta `$ (when considering nuclear collisions at a fixed value of $`b`$, and at a certain value of $`b,s,\theta `$ when integrating also over $`b`$. The probability of passing though different media depends on values of these parameters. Equations (11–13) change only a little. Lumping $`s,\theta `$, or $`b,s,\theta `$ into a single parameter $`x`$ and introducing probabilities $`P_1(x)`$, $`P_2(x)`$ and $`P_3(x)`$, we have in an obvious notation:

$$N_0=n_0(x)𝑑x$$

$$N_{P1}=n_0(x)P_1(x)𝑑x,N_{F1}=n_0(x)[1P_1(x)]𝑑x$$

$$N_{P2}=n_0(x)P_1(x)P_2(x)𝑑x,N_{F2}=n_0(x)P_1(x)[1P_2(x)]𝑑x$$

$$N_{P3}=n_0(x)P_1(x)P_2(x)P_3(x)𝑑x,N_{F3}=n_0(x)P_1(x)P_2(x)[1P_3(x)]𝑑x$$

and one can write again Eqs. (11)–(13).

When calculating the total $`J/\psi `$ suppression via an expression corresponding to Eq. (14), we get:

$$S(E_T)=\frac{N_{prod}^{J/\psi }N_{GH}^{J/\psi }N_{QGP}^{J/\psi }}{N_{prod}^{J/\psi }}.$$
(15)
Using Eqs. (4)–(6) and the identity $`1\mathrm{\Theta }(\kappa \kappa _{crit})=\mathrm{\Theta }(\kappa _{crit}\kappa )`$, we obtain from Eq. (15) the standard expression valid for $`J/\psi `$ suppression by nuclear absorption and by the QGP

$$S^{J/\psi }(E_T)=_0^{R_A}\frac{sds}{\sigma _{nn}}_0^{2\pi }𝑑\theta _{L_A(s)}^{L_A(s)}\rho _A\sigma _{nn}𝑑z_A_{L_B(s,\theta )}^{L_B(b,s,\theta )}\rho _B\sigma _{nn}𝑑z_B$$

$$e^{\rho _A\sigma _a[z_A+L_A(s)]}e^{\rho _B\sigma _a[z_B+L_B(b,s,\theta )]}\mathrm{\Theta }(\kappa _{crit}\kappa ).$$
(16)
## 5 Comments and conclusions

We have desribed here an alternative way of plotting the data and model calculations on $`J/\psi `$ suppression and presented a few illustrative examples. A more detailed analysis of $`J/\psi `$ suppression along this path would certainly require using more realistic Woods–Saxon nuclear densities, a more realistic relationship between $`b`$ and $`E_T`$, detailed analysis of $`J/\psi `$ suppression in collisions induced by lighter ions at $`200`$GeV per nucleon and a detailed evaluation of experimental errors. We are of the opinion that such an undertaking, although rather tedious, might help us to understand why models based on rather different assumptions lead to roughly similar results.

For models considering only nuclear absorption and $`J/\psi `$ interaction with co-movers, as in Refs.  the present scheme requires only minimal modifications. Equation (9) should be rewritten as

$$N_{prod}^{J/\psi }(E_T)=N_{exp}^{J/\psi }(E_T)+N_{NA}^{J/\psi }(E_T)+N_{co}^{J/\psi }(E_T).$$
(17)
where the last term corresponds to the number of $`J/\psi `$’s suppressed by the interaction with co-movers (with the gas of secondary hadrons). The term $`N_{prod}^{J/\psi }(E_T)`$ is given by Eqs. (3) and (4), $`N_{exp}^{J/\psi }(E_T)`$ by Eq.(5), $`N_{NA}^{J/\psi }(E_T)`$ by Eq. (6) and $`N_{co}^{J/\psi }(E_T)`$ is calculated, in analogy to Eq. (7), as

$$N_{co}^{J/\psi }(E_T)=\sigma _{nn}^{J/\psi }_0^{R_A}\frac{sds}{\sigma _{nn}}_0^{2\pi }𝑑\theta _{L_A(s)}^{L_A(s)}\rho _A\sigma _{nn}𝑑z_A_{L_B(s,\theta )}^{L_B(b,s,\theta )}\rho _B\sigma _{nn}𝑑z_B$$

$$e^{\rho _A\sigma _a[z_A+L_A(s)]}e^{\rho _B\sigma _a[z_B+L_B(b,s,\theta )]}\left(1\mathrm{exp}(𝑑\tau v\sigma _{co}\rho _{co})\right).$$
(18)
In Eq.(18), $`\sigma _{co}`$ is the cross-section for the disintegration of $`J/\psi `$ by the interaction with co-movers, $`\rho _{co}(\tau )`$ is the density of co-movers as a function of the proper time $`\tau `$ of $`J/\psi `$. For more details see Refs. .

For Monte Carlo microscopic models such as the one in Ref., an analytic recipe for separating different contributions cannot be written down, but the separation is most likely possible as well.

The alternative way of plotting and looking at data might be useful, but it will not resolve the physics problems. The question of whether the QGP is responsible for anomalous $`J/\psi `$ suppression observed by the NA50 Collaboration will probably and hopefully be resolved in one of the following ways (or by combination of some of them):

* Accurate determination of $`\sigma _a`$ by data on $`J/\psi `$ suppression at $`E_{Lab}=160`$GeV per nucleon. Using this value of $`\sigma _a`$ together with data on $`S(E_T)`$ in Pb–Pb interactions might reveal a threshold of ”missing contribution”, when plotted as in Fig. 2 above.
* Combining data on anomalous $`J/\psi `$ suppression with information on other signatures, including the data  on transverse momenta of surviving $`J/\psi `$’s and their analyses  and possible future data on nucleon number dependence on anomalous $`J/\psi `$ suppression. In this respect see e.g. Ref.. The explanation of data on anomalous $`J/\psi `$ suppression by QGP will only be generally accepted provided that the explanation will lead to experimentally verified predictions concerning the onset as function of A,B and $`E_T`$.
* Additional information on pre-resonance interaction with nucleons, see e.g. Refs. (”inverse experiment” and antiproton interactions in nuclei).
* Making use of recent results of lattice calculations , which are becoming more and more accurate and lead to the critical temperature for the phase transition of about $`173`$MeV. This represents a constraint for models with $`J/\psi `$ suppression by a gas of hadrons, since the energy density of the hadron gas higher than $`ϵ_{QGP}(T_c)`$ indicates the transition of the hadron gas into the QGP.

Acknowledgements
One of the authors (J.P.) is indebted to the CERN Theory Division for the hospitality extended to him. We would like to thank F. Karsch, C.Salgado and B. Tomášik for comments and helpful discussions.


format: text

text:
# Self-consistent quantum kinetics of condensate and non-condensate via a coupled equations of motion formalism

## 1 Introduction

The formulation of a consistent kinetic theory for the description of systems exhibiting Bose-Einstein condensation has been a very active area of research, since the very early days of the theoretical study of the behaviour of superfluid liquid helium . Such pioneering work for homogeneous systems has led to the establishment of the essential theoretical framework, which has since been extended to the case of trapped condensates. Due to the extreme weakness of their interactions, in comparison to those of liquid helium, the recently formed Bose-Einstein condensates in alkali gases are ideal systems for testing the validity of such theories. There are at the moment various non-equilibrium approaches to the dynamics of dilute, weakly-interacting, trapped Bose-Einstein condensates. Nonetheless, despite recent progress in the field, there exists to date no uniquely accepted non-equilibrium theory for the coupled dynamics of condensate and thermal cloud. On the one hand, we find the theories of Stoof  and Gardiner-Zoller et al. , based respectively on a non-perturbative Fokker-Planck equation for the non-equilibrium dynamics of the gas, and a master equation for the many-body density matrix. These have been compared to experiments , yielding very good agreement between them and with existing data for the issue of onset of condensation , and they further appear to be in reasonable qualitative agreement with earlier work of Kagan, Svistunov and Shlyapnikov . Other types of kinetic theories (except that of Castin and Dum ) assume the existence of a mean field for condensed atoms and deal perturbatively with more complex correlations around the mean field, based on suitable decoupling approximations, as first discussed by Kirkpatrick and Dorfmann . Such kinetic treatments have been recently discussed in the context of trapped gases by Proukakis and Burnett , Zaremba-Nikuni-Griffin  and Walser-Williams-Cooper-Holland . In a recent paper, Wachter et al.  have shown the equivalence of the non-equilibrium Green’s function approach originally proposed by Kadanoff and Baym , as applied to inhomogeneous systems by Imamovic-Tomasovic and Griffin , to the second-order gapless kinetic theory of Walser et al. . Central to this lies the fact that these theories yield the same second order damping rates, originally obtained in Beliaev’s pioneering work , and later extended by Popov , and more recently by Fedichev and Shlyapnikov .

Based on a hierarchy of coupled equations of motion formalism for the condensate mean field and fluctuations around it, the author and Burnett  have established an equation for condensate dynamics which extends beyond Hartree-Fock-Bogoliubov theory by the inclusion of triplet correlations; these correlations are related to the source field generating the condensate in the Green’s functions approach, as discussed, for example, in . Discussion of these equations has led to the ab initio introduction of an effective interaction in them , in a manner analogous to early diagrammatic work in the field . By focusing on static forms of these equations, Proukakis et al.  have further discussed gapless many-body theories, extending also beyond the bare t-matrix, establishing at the same time relations with other approaches . Extending such work, Morgan  has formulated a gapless, number-conserving, self-consistent second order perturbation theory for the excitations of inhomogeneous condensates, which extends beyond the quasiparticle basis. Morgan has discussed at length the issues of gaplessness and infrared and ultraviolet divergences which have plagued the study of both homogeneous  and inhomogeneous condensates from their very beginning.

In this paper, we extend the coupled equation of motion kinetic approach presented by the author elsewhere. In particular, in our preceeding work  we discussed how one formulates equations beyond the simplest HFB theory, by discussing the equations of motion of correlations of up to three single-particle operators. Despite containing all the essential dynamics for the condensate mean field, that treatment only dealt with correlations of four single-particle operators in their mean-field approximation, hence not being able to generate the dynamics of the non-condensate in a consistent manner. To achieve this, one must additionally look into the dynamic off-equilibrium contribution of four-operator correlations, which drive such correlations away from their respective equilibrium mean values due to interactions with the surrounding particles. Doing this above the transition point, the author  successfully obtained the classical Boltzmann equation, discussing at the same time how the actual interatomic potential becomes upgraded to an effective interaction. In this paper, we extend our earlier treatment by consistently treating, to second order in the effective (bare t-matrix) interaction, both the dynamics of the condensate and the non-condensate. This is based on the application of second order perturbation theory, starting from the usual Hartree-Fock-Bogoliubov (HFB), or quasiparticle basis, in much the same manner as the theory of Morgan . Working in such a basis, we explicitly re-derive the second order theory of Walser et al. , which is of a different nature, formulated essentially as a perturbative expansion of a many-body density operator around a self-consistently adjusted Gaussian and its correlations. This equivalence was anticipated, since both theories rely on the assumption of a few slowly-varying quantities, given in our case, by the HFB order parameters. By further noting that the theory of Walser et al.  was recently shown  to be equivalent to existing theories based on Green’s functions , we hence additionally establish a connection with such theories. Furthermore, since our perturbative treatment deals with the same hamiltonian as the equilibrium theory of Morgan , we further believe that this paper establishes an indirect link between the static self-consistent approach of Morgan and the dynamic theory of Walser et al. We therefore consider this work as another step towards the formulation of a consistent theory for the study of both static and dynamic processes in finite temperature Bose-Einstein condensates.

The layout of this paper is as follows. Firstly we present our formalism, giving explicit expressions for the treatment of the hamiltonian of the system and the corresponding energy functional. We then introduce, in Sec. III, the generalised matrix notation for the condensate and quasiparticle propagators and self energies, in the usual manner. In Sec. IV, we give explicitly the first and second order evolution of these generalised propagators (in terms of the effective interatomic potential). Sec. V shows clearly the equivalence of our treatment to other kinetic theories (focusing on the recently proposed theory of Walser et al.), and we conclude this paper with some general remarks (Sec. VI). In the appendices, we explain how to deal with more complex normal and anomalous averages for both non-equilibrium contributions (via their respective equations of motion given in Appendix A) and their corresponding equilibrium values (for which we give the required decoupling approximations in Appendix B).

## 2 Formalism

We begin our treatment with the usual binary interaction Hamiltonian,

$$\widehat{H}=\underset{rn}{}\mathrm{\Xi }_{rn}\widehat{a}_r^{}\widehat{a}_n+\frac{1}{2}\underset{rsmn}{}V_{rsmn}\widehat{a}_r^{}\widehat{a}_s^{}\widehat{a}_m\widehat{a}_n$$
(1)
which should incorporate most of the interesting physical processes occuring in dilute, weakly-interacting Bose-Einstein condensed gases. Here $`\mathrm{\Xi }=(\mathrm{}^2^2)/(2m)+V_{trap}(𝐫)`$ corresponds to the unperturbed hamiltonian in a harmonic trap, and $`V_{rsmn}`$ represents the symmetrised form of the interaction potential between a pair of particles, defined by $`V_{rsmn}=\frac{1}{2}\left\{rs|\widehat{V}|mn+rs|\widehat{V}|nm\right\}`$. Here $`|i=\psi _i(𝐫)`$ denotes a single-particle eigenstate of the trap, and the single-particle operators $`\widehat{a}_i`$ are related to the Bose field operator $`\widehat{\mathrm{\Psi }}(𝐫,t)`$ via $`\widehat{\mathrm{\Psi }}(𝐫,t)=_i\psi _i(𝐫)\widehat{a}_i(t)`$. Extending our earlier treatment regarding the ab initio introduction of an effective interaction , Morgan  has shown that, as long as high-lying states are adiabatically eliminated, and one is only interested in states up to a certain cut-off , the above hamiltonian can be equivalently written in terms of an effective re-summed two-body interaction. In fact, by numerically ensuring the independence of this effective interaction on the cut-off, this can be essentially replaced by the usual two-body t-matrix, often approximated by a local pseudopotential in three dimensions . All our subsequent expressions will be given in terms of $`V`$, bearing in mind that this essentially corresponds to such a re-summed two-body effective interaction. Our formalism relies on the existence of symmetry breaking, and hence we express the single-particle operators $`\widehat{a}_i`$ as

$$\widehat{a}_i=\widehat{a}_i+(\widehat{a}_i\widehat{a}_i)=z_i+\widehat{c}_i$$
(2)
One can formulate an infinite hierarchy of coupled equations of motion for the condensate mean field and fluctuations around this field. To make this problem tractable, certain approximations are required. Firstly, we note that for the very dilute gases we are dealing with, most of the interesting physics should be already apparent by their second order expressions (in the effective re-summed interatomic interaction), since such expressions contain both energy shifts and irreversible damping processes. Our theory further requires consistent decoupling approximations (discussed in Appendix B). Due to the large number of atoms in existing condensates and the mean field potentials they generate, the single-particle eigenstates of the harmonic trap are not necessarily the best states for a perturbative expansion. Hence, starting from an HFB basis, we incorporate all remaining effects as a perturbation. To this aim, we write our hamiltonian as

$$\widehat{H}=H_0+\widehat{H}_1^{^{}}+\widehat{H}_2^{^{}}+\widehat{H}_3^{^{}}+\widehat{H}_4^{^{}}$$
(3)
where

$$H_0=\underset{rn}{}\mathrm{\Xi }_{rn}z_r^{}z_n+\frac{1}{2}\underset{rsmn}{}V_{rsmn}z_r^{}z_s^{}z_mz_n$$
(4)
is the unpeturbed hamiltonian in a trap, additionally dressed by the condensate mean field, and the remaining operators are given by

$`\widehat{H}_1^{^{}}`$ $`={\displaystyle \underset{rn}{}}\mathrm{\Xi }_{rn}\left[\widehat{c}_r^{}z_n+z_r^{}\widehat{c}_n\right]+{\displaystyle \underset{rsmn}{}}V_{rsmn}\left[\widehat{c}_r^{}\left(z_s^{}z_m+\widehat{c}_s^{}\widehat{c}_m\right)z_n+z_r^{}\left(z_s^{}z_m+\widehat{c}_s^{}\widehat{c}_m\right)\widehat{c}_n\right]`$ (5)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{rsmn}{}}V_{rsmn}\left[\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_mz_n+\widehat{c}_r^{}z_s^{}\widehat{c}_m\widehat{c}_n\right]`$
$`\widehat{H}_2^{^{}}`$ $`={\displaystyle \underset{rn}{}}\mathrm{\Xi }_{rn}\widehat{c}_r^{}\widehat{c}_n+{\displaystyle \underset{rsmn}{}}V_{rsmn}\left[2\widehat{c}_r^{}\left(z_s^{}z_m+\widehat{c}_s^{}\widehat{c}_m\right)\widehat{c}_n\widehat{c}_r^{}\widehat{c}_n\widehat{c}_s^{}\widehat{c}_m\right]`$ (6)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{rsmn}{}}V_{rsmn}\left[\widehat{c}_r^{}\widehat{c}_s^{}\left(z_mz_n+\widehat{c}_m\widehat{c}_n\right)+\left(z_r^{}z_s^{}+\widehat{c}_r^{}\widehat{c}_s^{}\right)\widehat{c}_m\widehat{c}_n\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n\right]`$
$`\widehat{H}_3^{^{}}`$ $`={\displaystyle \underset{rsmn}{}}V_{rsmn}\left[\left(\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m2\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\right)z_n\right]`$ (7)
$`+{\displaystyle \underset{rsmn}{}}V_{rsmn}\left[z_r^{}\left(\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n2\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n\right)\right]`$
$`\widehat{H}_4^{^{}}`$ $`={\displaystyle \frac{1}{2}}{\displaystyle \underset{rsmn}{}}V_{rsmn}\left[\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n4\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n\right]`$ (8)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{rsmn}{}}V_{rsmn}\left[2\widehat{c}_r^{}\widehat{c}_n\widehat{c}_s^{}\widehat{c}_m+\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n\right]`$
The preceeding operator-dependent hamiltonians ($`\widehat{H}_i^{^{}}`$) have been denoted by a prime, indicating that these are not the usual hamiltonians one would obtain from equation (1) upon writing $`\widehat{a}_i=z_i+\widehat{c}_i`$, but they have been modified by the application of the following mean field approximations (see e.g. )

$$\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\widehat{c}_r^{}\widehat{c}_m\widehat{c}_s^{}+\widehat{c}_s^{}\widehat{c}_m\widehat{c}_r^{}+\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m$$
(9)
$$\widehat{c}_r^{}\widehat{c}_m\widehat{c}_n\widehat{c}_r^{}\widehat{c}_m\widehat{c}_n+\widehat{c}_r^{}\widehat{c}_n\widehat{c}_m+\widehat{c}_r^{}\widehat{c}_m\widehat{c}_n$$
(10)
$`\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n`$ $`\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n+\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_n\widehat{c}_m+\widehat{c}_s^{}\widehat{c}_r^{}\widehat{c}_m\widehat{c}_n+\widehat{c}_s^{}\widehat{c}_r^{}\widehat{c}_n\widehat{c}_m+\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n+\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n`$ (11)
$`\left(2\widehat{c}_r^{}\widehat{c}_n\widehat{c}_s^{}\widehat{c}_m+\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n\right)`$
Such mean field approximations are conventionally performed in HFB treatments aimed at reducing the hamiltonian of the system to the usual quadratic form $`\left(\widehat{H}_0+\widehat{H}_1^{^{}}+\widehat{H}_2^{^{}}\right)`$, ignoring completely contributions beyond the mean field labelled above as $`\left(\widehat{H}_3^{^{}}+\widehat{H}_4^{^{}}\right)`$, since, by definition, $`\left(\widehat{H}_3^{^{}}+\widehat{H}_4^{^{}}\right)=0`$ when taking averages over HFB eigenstates. Such a hamiltonian can be diagonalised by a transformation to quasiparticle operators , which mix the shifted single-particle creation and annihilation operators $`\widehat{c}_i^{()}`$. In such a theory, the effects of thermal atoms are treated only in an approximate manner, and this is known to lead to an unphysical gap in the spectrum of elementary excitations in the homogeneous limit . This gap can be attributed to the inconsistent treatment of interactions between condensed and uncondensed components of the system and, more specifically, to the fact that the above approximation for three-particle operators is not physically justified . For consistency we note here that the approximation on four operators arises from an extension of Wick’s theorem  (which is rigorous for equilibrium averages) to the non-equilibrium case.

To avoid such complications, we will hence work with an “unperturbed” hamiltonian consisting of the usual quadratic HFB hamiltonian $`\left(H_0+\widehat{H}_1^{^{}}+\widehat{H}_2^{^{}}\right)`$, but also keep the terms $`\left(\widehat{H}_3^{^{}}+\widehat{H}_4^{^{}}\right)`$ such that our hamiltonian is still exact, as also done by Morgan . As we will show explicitly in this paper, this implies that all first-order contributions (assuming no steady state anomalous correlations beyond the pair correlation) arise from the HFB hamiltonian $`\left(H_0+\widehat{H}_1^{^{}}+\widehat{H}_2^{^{}}\right)`$, whereas the remaining parts of the hamiltonian $`\left(\widehat{H}_3^{^{}}+\widehat{H}_4^{^{}}\right)`$ are solely responsible for all higher order collisional processes, damping rates, etc. In this paper, we will focus on the application of the above hamiltonian to the dynamical domain, somewhat similarly to the approach of Zaremba-Nikuni-Griffin . We start our treatment from a suitably generalised version of our earlier kinetic equations, and gradually build upon them, which leads to the explicit re-derivation of the kinetic theory of Walser et al. .

### 2.1 Some Definitions

Applying perturbation theory on the HFB hamiltonian $`\left(H_0+\widehat{H}_1^{^{}}+\widehat{H}_2^{^{}}\right)`$ amounts to assuming that the correlations of up to two single-particle operators can be thought of as slowly-varying, whereas all other (higher order) correlations will be treated perturbatively (beyond their corresponding equilibrium contributions). We thus define the lowest order normal and anomalous averages of single-particle operators as

$$\rho _{ji}=\widehat{c}_i^{}\widehat{c}_j$$
(12)
$$\kappa _{ji}=\widehat{c}_i\widehat{c}_j$$
(13)
These quantities, along with the condensate mean field $`z_i`$ form the set of order parameters of the HFB basis . This is analogous to the assumption made by Walser et al.  that such quantities form the set of relevant (master) variables for describing the evolution of the system on a coarse-grained timescale. For a binary hamiltonian, as already noted in , a consistent treatment of the dynamic evolution of such quantities, to second order in the (effective) interatomic potential, requires an in-depth analysis of the evolution of all correlations, normal and anomalous, containing up to four single-particle shifted operators. As is well-known, higher order correlations containing an even number of shifted single-particle operators do not vanish at equilibrium, but actually acquire a steady state value consistent with Wick’s theorem. This states that, at thermal equilibrium, such higher order averages can be decomposed into products of lower order averages, which in our case correspond to the slowly-varying HFB order parameters $`z`$, $`\rho `$ and $`\kappa `$. To proceed further with our coupled equation of motion methodology, we must hence treat such quantities in two steps : firstly we split them into their respective ‘equilibrium’ and ‘non-equilibrium’ components; subsequently, we decompose all equilibrium higher-order averages into products of HFB parameters, while at the same time perturbatively eliminating their rapidly evolving parts which drive them out of equilibrium. The non-equilibrium quantities we will be concerned with are defined by

$$\sigma _{rsmn}=\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n^{(0)}=\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n\left(\rho _{mr}\rho _{ns}+\rho _{nr}\rho _{ms}+\kappa _{mn}\kappa _{rs}^{}\right)$$
(14)
and

$$\xi _{rsmn}=\widehat{c}_r^{}\widehat{c}_s\widehat{c}_m\widehat{c}_n\widehat{c}_r^{}\widehat{c}_s\widehat{c}_m\widehat{c}_n^{(0)}=\widehat{c}_r^{}\widehat{c}_s\widehat{c}_m\widehat{c}_n\left(\rho _{sr}\kappa _{mn}+\rho _{mr}\kappa _{ns}+\rho _{nr}\kappa _{ms}\right)$$
(15)
This appearance of ‘equilibrium’ values is not encountered by triplet correlations (or higher order odd single-particle operator averages). These are assumed to vanish identically at equilibrium (due to our choice of the HFB basis as a good basis for describing the equilibrium properties of the system) and can therefore be treated perturbatively. In this case, we define

$$\lambda _{rmn}=\widehat{c}_r^{}\widehat{c}_m\widehat{c}_n$$
(16)
$$\gamma _{smn}=\widehat{c}_s\widehat{c}_m\widehat{c}_n$$
(17)
### 2.2 Energy functional

Having defined our notation, we can now explicitly write down the energy functional $`\widehat{E}`$ of our system based on the hamiltonian (1). Its exact form is given by

$$\widehat{E}=\widehat{H}=\widehat{E}_{HFB}^{(1)}+\widehat{E}^{(2)}$$
(18)
Here we have defined the HFB energy functional in the usual manner

$`\widehat{E}_{HFB}^{(1)}`$ $`\left(H_0+\widehat{H}_1^{^{}}+\widehat{H}_2^{^{}}\right)={\displaystyle \underset{rn}{}}\mathrm{\Xi }_{rn}\left(z_r^{}z_n+\rho _{nr}\right)`$ (19)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{rsmn}{}}V_{rsmn}\left[z_r^{}z_s^{}z_mz_n+4\rho _{nr}z_s^{}z_m+2\rho _{nr}\rho _{ms}+\kappa _{rs}^{}z_mz_n+z_r^{}z_s^{}\kappa _{mn}+\kappa _{rs}^{}\kappa _{mn}\right]`$
and the additional ‘perturbative’ beyond-HFB functional

$`\widehat{E}^{(2)}`$ $``$ $`\left(\widehat{H}_3^{^{}}+\widehat{H}_4^{^{}}\right)={\displaystyle \frac{1}{2}}{\displaystyle \underset{rsmn}{}}V_{rsmn}\left[2\lambda _{mrs}^{}z_n+2\lambda _{smn}z_r^{}+\sigma _{rsmn}\right]`$ (20)
Consideration of $`\widehat{E}_{HFB}^{(1)}`$ leads to the well-known (reversible) first order expression for the dynamics of the HFB order parameters. The second term, $`\widehat{E}^{(2)}`$ has no lowest order contribution in the HFB basis (since the HFB basis, by definition, does not allow for a non-vanishing steady state value for the triplets) and is hence conventionally ignored.

In this paper we demonstrate that it is precisely this second contribution to the energy functional which generates all second (and higher order) collisional contributions and should, thus, not be ignored . In fact, we will show explicitly that the second order energy functional discussed by Walser et al.  is nothing but our functional $`\left(\widehat{E}_{HFB}^{(1)}+\widehat{E}^{(2)}\right)`$ when treated to second order, thus justifying its particular form. To obtain the desired higher order contributions in a perturbative expansion, we thus require the equations of motion of $`\lambda `$ and $`\sigma `$, as given in Appendix A.

Making use of these expressions and re-arranging, we obtain for the second order expression of the energy functional

$`\widehat{E}^{(2)}`$ $`=(i){\displaystyle \underset{rsmn}{}}{\displaystyle \underset{pqlt}{}}V_{rsmn}\stackrel{~}{V}_{pqlt}`$ (29)
$`\times \left\{\begin{array}{c}\left[\left(\rho _{mp}+\delta _{mp}\right)\left(\rho _{nq}+\delta _{nq}\right)\rho _{ts}\rho _{lr}\rho _{mp}\rho _{nq}\left(\rho _{ts}+\delta _{ts}\right)\left(\rho _{lr}+\delta _{lr}\right)\right]\hfill \\ +2\left[\left(\rho _{mp}+\delta _{mp}\right)\left(z_q^{}z_n\right)\rho _{ts}\rho _{lr}\rho _{mp}\left(z_q^{}z_n\right)\left(\rho _{ts}+\delta _{ts}\right)\left(\rho _{lr}+\delta _{lr}\right)\right]\hfill \\ +2\left[\left(\rho _{mp}+\delta _{mp}\right)\left(\rho _{nq}+\delta _{nq}\right)\left(z_s^{}z_t\right)\rho _{lr}\rho _{mp}\rho _{nq}\left(z_s^{}z_t\right)\left(\rho _{lr}+\delta _{lr}\right)\right]\hfill \\ +4\left[\left(\rho _{mp}+\delta _{mp}\right)\kappa _{nt}\kappa _{qs}^{}\rho _{lr}\rho _{mp}\kappa _{nt}\kappa _{qs}^{}\left(\rho _{lr}+\delta _{lr}\right)\right]\hfill \\ +4\left[\left(\rho _{mp}+\delta _{mp}\right)\kappa _{nt}\kappa _{qs}^{}\left(z_r^{}z_l\right)\rho _{mp}\kappa _{nt}\kappa _{qs}^{}\left(z_r^{}z_l\right)\right]\hfill \\ +4\left[\left(z_p^{}z_m\right)\kappa _{nt}\kappa _{qs}^{}\rho _{lr}\left(z_p^{}z_m\right)\kappa _{nt}\kappa _{qs}^{}\left(\rho _{lr}+\delta _{lr}\right)\right]\hfill \\ +4\left[\left(\rho _{mp}+\delta _{mp}\right)\left(z_nz_t\right)\kappa _{qs}^{}\rho _{lr}\rho _{mp}\left(z_nz_t\right)\kappa _{qs}^{}\left(\rho _{lr}+\delta _{lr}\right)\right]\hfill \\ +4\left[\left(\rho _{mp}+\delta _{mp}\right)\kappa _{nt}\left(z_q^{}z_s^{}\right)\rho _{lr}\rho _{mp}\kappa _{nt}\left(z_q^{}z_s^{}\right)\left(\rho _{lr}+\delta _{lr}\right)\right]\hfill \end{array}\right\}`$
where $`\delta _{ij}`$ is the usual Cronecker delta.

We stress that the above expression only corresponds to the second order contribution of the energy functional $`\widehat{E}`$, and not to the exact energy functional given by equations (18)-(20). In particular, we note that this second order expression additionally contains higher order averages such as $`\gamma `$, $`\lambda `$, $`\sigma `$ and $`\xi `$ which play the crucial role of generating all higher order contributions to the energy functional (but do not modify the system evolution to second order, just as $`\widehat{E}^{(2)}`$ has no effect on first order expressions for condensate / non-condensate evolution).

In the above expression, we have defined the approximately energy-conserving matrix element

$$\stackrel{~}{V}_{pqlt}=𝑑t^{^{}}e^{i\left(\stackrel{~}{\omega }_l+\stackrel{~}{\omega }_t\stackrel{~}{\omega }_p\stackrel{~}{\omega }_q\right)(tt^{^{}})}=\pi \delta (\mathrm{\Delta }\stackrel{~}{\omega })iP\left(\frac{1}{\mathrm{\Delta }\stackrel{~}{\omega }}\right)$$
(30)
where $`\mathrm{\Delta }\stackrel{~}{\omega }=\left(\stackrel{~}{\omega }_l+\stackrel{~}{\omega }_t\stackrel{~}{\omega }_p\stackrel{~}{\omega }_q\right)`$ and $`\stackrel{~}{\omega }_i`$ denotes the energy of level $`i`$ dressed by HFB mean field potentials (i.e. expectation value in the quasiparticle hamiltonian $`\left(H_0+\widehat{H}_1^{^{}}+\widehat{H}_2^{^{}}\right)`$). We note that the particular structure of the energy-conserving indices in the above expression for the energy functional (i.e. that they correspond precisely to the second matrix element $`V`$ of $`\widehat{E}_2`$) is a direct consequence of the Markov approximation which has been additionally imposed here (see section IV B). However, the Markov approximation only affects the indices appearing in the quasiparticle energies of the exponential in the integrand, and not the general structure of the energy functional.

Apart from the first line of equation (29) containing the usual Boltzmann scattering rates for normal (uncondensed) averages $`\rho `$, we find terms involving condensate populations $`(z^{}z)`$, as well as condensed ($`zz`$) and uncondensed ($`\kappa `$) anomalous averages. The appealing feature of this energy functional (which, admittedly, is not apparent from its first order expression of equation (20)) is its explicit symmetry with respect to condensate and non-condensate contributions. For example, additionally to the usual ‘classical’ Boltzmann rates (first line of equation (29)), we obtain all similar contributions in which one of the normal uncondensed averages $`\rho `$ is replaced by $`(z^{}z)`$. Even more appealing is the fact that this further holds when considering anomalous averages. Hence, in addition to terms $`\left[(\rho +1)\kappa \kappa ^{}\rho \right]`$, $`\widehat{E}^{(2)}`$ contains terms $`\left[(\rho +1)\kappa (zz)^{}\rho \right]`$, $`\left[(\rho +1)(zz)\kappa \rho \right]`$, $`\left[(\rho +1)\kappa \kappa ^{}(z^{}z)\right]`$ and $`\left[(z^{}z)\kappa \kappa ^{}\rho \right]`$, etc. We note that it is precisely the choice of the HFB basis for our unperturbed hamiltonian which prohibits, in the final expressions, the appearance of multiple products of normal or anomalous condensate averages beyond their simplest forms $`(z^{}z)`$, $`(zz)`$ or $`(z^{}z^{})`$. This simplification would not arise if we were working with a simpler unperturbed hamiltonian, such as the one describing bare trap eigenstates, the Gross-Pitaevskii basis (dressed only by the condensate mean field) or the Hartree-Fock basis (additionally dressed by the normal non-condensate average). Finally, as already pointed out by Walser et al. , although the in and out rates for the normal uncondensed component differ due to the process of bosonic enhancement (modifying $`\rho `$ to $`(\rho +1)`$), the mean field $`(z^{}z)`$ is never bosonically enhanced, and can thus be thought of as a classical field.

## 3 Generalised matrix notation

In order to present our formalism in the most general manner and establish a straighforward link with existing theories, we will henceforth work with generalised density matrices and hamiltonians, in a manner which clearly distinguishes between first order HFB and higher order perturbative results. We hence define the generalised condensate matrix $`R_c`$ ($`2n\times 1`$) and the generalised quasiparticle density matrix $`R_e`$ ($`2n\times 2n`$) by

$`R_c=\left(\begin{array}{c}z\\ z^{}\end{array}\right)R_e=\left(\begin{array}{c}\rho \\ \kappa ^{}\end{array}\begin{array}{c}\kappa \\ (\rho ^{}+\text{1}\text{1})\end{array}\right)`$ (37)
where 11 is the unity matrix of the n-dimensional Fock space. These quantities can be thought of as the propagators for the condensed and uncondensed parts of the system, and are analogous to the corresponding generalised propagators conventionally employed in Green’s functions formulations . We note that, although $`R_e`$ includes all of the effects of a quasiparticle basis, these are explicitly given in terms of correlations of single-particle operators, with an explicit transformation to quasiparticle basis being beyond the scope of this paper (see e.g. ).

Correspondingly, we define the generalised condensate and quasiparticle hamiltonians as

$`H_c=\left(\begin{array}{c}h^{(c)}\\ (\mathrm{\Delta }^{(c)})^{}\end{array}\begin{array}{c}\mathrm{\Delta }^{(c)}\\ (h^{(c)})^{}\end{array}\right)H_e=\left(\begin{array}{c}h\\ \mathrm{\Delta }^{}\end{array}\begin{array}{c}\mathrm{\Delta }\\ h^{}\end{array}\right)`$ (46)
which can be thought of as the generalised self-energies for the above Green’s functions.

These generalised HFB matrix hamiltonians are already known from variational approaches, in which they can be readily obtained by variation of the energy functional $`\widehat{E}_{HFB}^{(1)}`$ with respect to the corresponding HFB propagators. In particular, one finds

$`\begin{array}{c}h_{ij}\frac{\widehat{E}_{HFB}^{(1)}}{\rho _{ji}(t)}=h_{ij}^c+_{kl}V_{iklj}\left(z_k^{}z_l\right)=i|\widehat{\mathrm{\Xi }}|j+2_{kl}V_{iklj}\left[z_k^{}z_l+\rho _{lk}\right]\\ \mathrm{\Delta }_{ij}\frac{\widehat{E}_{HFB}^{(1)}}{\kappa _{ji}^{}(t)}=\mathrm{\Delta }_{ij}^c+_{kl}V_{ijkl}\left(z_kz_l\right)=_{kl}V_{ijkl}\left[z_kz_l+\kappa _{kl}\right]\end{array}`$ (49)
For consistency, and ease of subsequent comparison to the expressions of Walser et al. , in this paper we have defined the generalised hamiltonians $`H_c`$ and $`H_e`$ in a slightly different manner from those of our earlier treatment . This difference is only a matter of notation and arises from how the Pauli matrix $`\sigma ^{(3)}=\left(\begin{array}{c}\text{1}\text{1}\hfill \\ 0\hfill \end{array}\begin{array}{c}0\\ \text{1}\text{1}\end{array}\right)`$ is incorporated in our equations<sup>1</sup><sup>1</sup>1We also note that, despite the different notation, our previous presentation contained a minor error, in that the left hand sides of equations (32) and (2.32) of  and  respectively, should contain an additional factor of $`\eta =\sigma ^{(3)}`$).. We also point out the useful identities $`h_{ij}^{}=h_{ji}`$, $`\mathrm{\Delta }_{ij}=\mathrm{\Delta }_{ji}`$, $`\rho _{ji}=\rho _{ij}^{}`$ and $`\kappa _{kj}=\kappa _{jk}`$, whereas all higher order correlations are symmetric with respect to the interchange of indices labelling operators of the same ‘type’ (i.e. creation $`\widehat{c}^{}`$ or annihilation $`\widehat{c}`$).

## 4 Coupled evolution of condensate and quasiparticle propagators

In this section we discuss the dynamic evolution of the condensate and quasiparticle propagators defined above, by employing the Heisenberg equation of motion (setting $`\mathrm{}=1`$ for simplicity)

$$i\frac{d}{dt}\widehat{O}=[\widehat{O},\widehat{H}]$$
(50)
for the mean value of a general operator $`\widehat{O}`$.

### 4.1 Exact first order evolution

To first order in the effective potential, the equations of motion for the condensed and uncondensed propagators acquire the following exact form

$`i{\displaystyle \frac{dR_c}{dt}}=H_cR_c+J`$ (51)
$$i\frac{dR_e}{dt}=\left(HRRH^{}\right)+K$$
(52)
In each of the above equations, the first terms correspond to the HFB contributions, arising entirely from the hamiltonian $`\left(H_0+\widehat{H}_1^{^{}}+\widehat{H}_2^{^{}}\right)`$. The static form of such equations can be readily derived variationally from $`\widehat{E}_{HFB}^{(1)}`$ . Our approach extends beyond such treatments, in that it further generates the ‘non-equilibrium’ matrices $`J`$ and $`K`$ for condensate and non-condensate evolution. These terms, arising solely from $`\widehat{E}^{(2)}`$, depend on averages of higher order than the HFB parameters, and are thus conventionally set to zero at equilibrium. In this work, these matrices are treated perturbatively, thus generating terms of higher order in the potential, which can still depend on the non-vanishing equilibrium quantities $`z`$, $`\rho `$ and $`\kappa `$. This shows that such contributions, and hence the beyond-HFB energy functional $`\widehat{E}^{(2)}`$ cannot be ignored even in static treatments, a point discussed clearly in the related work of Morgan . To further establish a connection to kinetic treatments involving a damped nonlinear Schrödinger equation for the condensate mean field at finite temperatures, we note that it is precisely the presence of such additional contributions (labelled here by the kinetics matrix $`J`$) which lead to the damping term. In fact, by carefully identifying the physical origin of such a damping term (in a manner somewhat analogous to the present treatment), Zaremba et al.  and Bijlsma et al.  have respectively performed detailed studies of the hydrodynamic regime and the onset of condensation. It further appears that the kinetic matrix $`J`$ is also taken into account in treatments based on linearised equations of motion, either explicitly , or implicitly .

Coming back to our formalism, the exact form of these kinetic matrices is given by

$`J=\left(\begin{array}{c}L\\ L^{}\end{array}\right)`$ (55)
$`K=\left(\begin{array}{c}(M\stackrel{~}{M}^{})\\ (N+\stackrel{~}{N})^{}\end{array}\begin{array}{c}(N+\stackrel{~}{N})\\ (M\stackrel{~}{M}^{})^{}\end{array}\right)`$ (60)
where

$$L_{ji}=\underset{smn}{}V_{jsmn}\lambda _{smn}$$
(61)
$$M_{ji}=\underset{smn}{}V_{jsmn}[\sigma _{ismn}+\lambda _{imn}z_s^{}+\lambda _{mis}^{}z_n+\lambda _{nis}^{}z_m]$$
(62)
$$N_{ji}=\underset{smn}{}V_{jsmn}[\xi _{ismn}+\lambda _{sim}z_n+\lambda _{sin}z_m+\gamma _{imn}z_s^{}]$$
(63)
and $`\stackrel{~}{N}`$ represents the transpose of matrix $`N`$. Note that the above definitions follow the corresponding definitions of , which they generalise by additional inclusion of quartic terms $`\sigma `$ and $`\xi `$ required for the consistent treatment of the non-condensate kinetics.

### 4.2 Second order collisional integrals

To proceed further and derive all relevant second order collisional integrals, we merely need to derive the corresponding equations of motion for the above kinetic matrices $`J`$ and $`K`$. This requires the evolution of higher order normal and anomalous single-particle operator correlations, as given in Appendix A, and the application of consistent decoupling approximations, given in Appendix B, for truncating the infinitely coupled equation of motion hierarchy. Our choice to apply perturbation theory starting from an HFB unperturbed basis, implies that the only effect of the hamiltonian $`\left(H_0+\widehat{H}_1^{^{}}+\widehat{H}_2^{^{}}\right)`$, beyond yielding the first order evolution, is to define the renormalised unperturbed eigenenergies $`\stackrel{~}{\omega }_i`$ which appear in all higher-order expressions. In particular, we note that these eigenenergies dressed by HFB mean fields are implicit in the second-order collisional evolution of condensate / non-condensate propagators via the approximately energy-conserving matrix element $`\stackrel{~}{V}_{pqlt}`$ of equation (30). To keep our subsequent notation compact, we henceforth suppress the free evolution of matrices $`J`$ and $`K`$ (see also Appendix A). Moreover, in giving their respective lowest order collisional evolution, we additionally impose the Markov approximation on the slowly-evolving HFB order parameters $`z`$, $`\rho `$ and $`\kappa `$ appearing on the right hand sides of equations (70)-(75). This relies on the assumption that such quantities evolve significantly only over many complete collisional events, and thus decay much slower than all other interparticle correlations. We expect this to be a valid assumption for the dilute, weakly-interacting systems under consideration. (A detailed discussion of non-Markovian kinetic treatments which also deal explicitly with initial correlations in a system can be found in .) Within the above simplifications, we thus obtain

$`i{\displaystyle \frac{dJ}{dt}}=\left(\begin{array}{c}\mathrm{\Gamma }_z\\ \left[\mathrm{\Gamma }_z^{}\right]^{}\end{array}\begin{array}{c}\mathrm{\Gamma }_z^{}\\ \left[\mathrm{\Gamma }_z\right]^{}\end{array}\right)\left(\begin{array}{c}z\\ z^{}\end{array}\right)`$ (70)
$`i{\displaystyle \frac{dK}{dt}}=\left(\begin{array}{c}\left[\mathrm{\Gamma }_\rho \rho +\mathrm{\Gamma }_\kappa \kappa ^{}+I_\rho \right]\\ \left[\mathrm{\Gamma }_\rho \kappa +\mathrm{\Gamma }_\kappa \left(\rho +\text{1}\text{1}\right)^{}+I_\kappa \right]^{}\end{array}\begin{array}{c}\left[\mathrm{\Gamma }_\rho \kappa +\mathrm{\Gamma }_\kappa \left(\rho +\text{1}\text{1}\right)^{}+I_\kappa \right]\\ \left[\mathrm{\Gamma }_\rho \rho +\mathrm{\Gamma }_\kappa \kappa ^{}+I_\rho \right]^{}\end{array}\right)+h.c.`$ (75)
where $`h.c.`$ denotes the hermitian conjugate and we have used a compact matrix notation, writing, for example, $`\left[\left(\mathrm{\Gamma }_\rho \right)\rho \right]_{ji}=_l\left(\mathrm{\Gamma }_\rho \right)_{jl}\rho _{li}`$, and similarly for the other contributions.

We find that the rates $`\mathrm{\Gamma }_\rho `$ and $`\mathrm{\Gamma }_\kappa `$ appearing in the above equation can be calculated variationally from the explicit second order expression for the energy functional of equation (29), as

$$\left(\mathrm{\Gamma }_\rho \right)_{rl}=i\left(\frac{\widehat{E}^{(2)}}{\rho _{lr}}\right)$$
(76)
$$\left(\mathrm{\Gamma }_\kappa \right)_{sq}=i\left(\frac{\widehat{E}^{(2)}}{\kappa _{qs}^{}}\right)$$
(77)
It is also interesting to note that

$$\left(I_\rho \right)_{rl}=i\left(\frac{\widehat{E}^{(2)}}{\left(\rho _{lr}+\delta _{lr}\right)}\right)$$
(78)
whereas we have not been able to obtain a correspondingly simple expression for $`I_\kappa `$.

For the condensate rates $`\mathrm{\Gamma }_{z^{()}}`$ we find

$$i\left(\frac{\widehat{E}^{(2)}}{z_i^{}}\right)=\underset{l}{}\left(\mathrm{\Gamma }_z\right)_{il}z_l+\underset{q}{}\left(\mathrm{\Gamma }_z^{}\right)_{iq}z_q^{}=\left[\frac{dz_i}{dt}\right]_{V^2}$$
(79)
where the subscript $`V^2`$ has been used to indicate second-order contributions to the evolution of the condensate mean field. Unlike the case of the condensate propagator, the presence of the non-vanishing matrices $`I_\rho `$, $`I_\kappa `$ in equation (75) indicate that, in this case, one may not be able to define variationally a simple generalised hamiltonian including all rates (as was done earlier for the first order expressions).

Below we define explicitly the form of all ($`n\times n`$) matrices appearing in the above equations, in terms of their respective elements

$`\left(\mathrm{\Gamma }_z\right)_{il}=2{\displaystyle \underset{smn}{}}{\displaystyle \underset{pqlt}{}}V_{jsmn}\stackrel{~}{V}_{pqlt}\left\{\begin{array}{c}\left[\left(\rho _{mp}+\delta _{mp}\right)\left(\rho _{nq}+\delta _{nq}\right)\rho _{ts}\rho _{mp}\rho _{nq}\left(\rho _{ts}+\delta _{ts}\right)\right]\hfill \\ +2\left[\left(\rho _{mp}+\delta _{mp}\right)\kappa _{nt}\kappa _{qs}^{}\rho _{mp}\kappa _{nt}\kappa _{qs}^{}\right]\hfill \end{array}\right\}`$ (82)
$`\left(\mathrm{\Gamma }_z^{}\right)_{iq}=2{\displaystyle \underset{smn}{}}{\displaystyle \underset{pqlt}{}}V_{jsmn}\stackrel{~}{V}_{pqlt}\left\{2\left[\left(\rho _{mp}+\delta _{mp}\right)\kappa _{nt}\rho _{ls}\rho _{mp}\kappa _{nt}\left(\rho _{ls}+\delta _{ls}\right)\right]\right\}`$ (83)
$`\left(\mathrm{\Gamma }_\rho \right)_{jl}=2{\displaystyle \underset{smn}{}}{\displaystyle \underset{pqlt}{}}V_{jsmn}\stackrel{~}{V}_{pqlt}\left\{\begin{array}{c}\left[\left(\rho _{mp}+\delta _{mp}\right)\left(\rho _{nq}+\delta _{nq}\right)\rho _{ts}\rho _{mp}\rho _{nq}\left(\rho _{ts}+\delta _{ts}\right)\right]\hfill \\ +2\left[\left(\rho _{mp}+\delta _{mp}\right)\left(z_q^{}z_n\right)\rho _{ts}\rho _{mp}\left(z_q^{}z_n\right)\left(\rho _{ts}+\delta _{ts}\right)\right]\hfill \\ +\left[\left(\rho _{mp}+\delta _{mp}\right)\left(\rho _{nq}+\delta _{nq}\right)\left(z_s^{}z_t\right)\rho _{mp}\rho _{nq}\left(z_s^{}z_t\right)\right]\hfill \\ +2\left[\left(\rho _{mp}+\delta _{mp}\right)\kappa _{nt}\kappa _{qs}^{}\rho _{mp}\kappa _{nt}\kappa _{qs}^{}\right]\hfill \\ +2\left[\left(z_p^{}z_m\right)\kappa _{nt}\kappa _{qs}^{}\left(z_p^{}z_m\right)\kappa _{nt}\kappa _{qs}^{}\right]\hfill \\ +2\left[\left(\rho _{mp}+\delta _{mp}\right)\left(z_nz_t\right)\kappa _{qs}^{}\rho _{mp}\left(z_nz_t\right)\kappa _{qs}^{}\right]\hfill \\ +2\left[\left(\rho _{mp}+\delta _{mp}\right)\kappa _{nt}\left(z_q^{}z_s^{}\right)\rho _{mp}\kappa _{nt}\left(z_q^{}z_s^{}\right)\right]\hfill \end{array}\right\}`$ (91)
$`\left(\mathrm{\Gamma }_\kappa \right)_{jq}=2{\displaystyle \underset{smn}{}}{\displaystyle \underset{pqlt}{}}V_{jsmn}\stackrel{~}{V}_{pqlt}\left\{\begin{array}{c}2\left[\left(\rho _{mp}+\delta _{mp}\right)\kappa _{nt}\rho _{ls}\rho _{mp}\kappa _{nt}\left(\rho _{ls}+\delta _{ls}\right)\right]\hfill \\ +2\left[\left(\rho _{mp}+\delta _{mp}\right)\kappa _{nt}\left(z_s^{}z_l\right)\rho _{mp}\kappa _{nt}\left(z_s^{}z_l\right)\right]\hfill \\ +2\left[\left(z_p^{}z_m\right)\kappa _{nt}\rho _{ls}\left(z_p^{}z_m\right)\kappa _{nt}\left(\rho _{ls}+\delta _{ls}\right)\right]\hfill \\ +2\left[\left(\rho _{mp}+\delta _{mp}\right)\left(z_nz_t\right)\rho _{ls}\rho _{mp}\left(z_nz_t\right)\left(\rho _{ls}+\delta _{ls}\right)\right]\hfill \end{array}\right\}`$ (96)
$`\left(I_\rho \right)_{ji}`$ $`=2{\displaystyle \underset{smn}{}}{\displaystyle \underset{pqlt}{}}V_{jsmn}\stackrel{~}{V}_{pqlt}\left(\delta _{li}\right)`$ (99)
$`\times \left\{\begin{array}{c}\rho _{mp}\rho _{nq}\left(\rho _{ts}+\delta _{ts}\right)+2\left(z_p^{}z_m\right)\rho _{nq}\left(\rho _{ts}+\delta _{ts}\right)+\rho _{mp}\rho _{nq}\left(z_s^{}z_t\right)\hfill \\ +2\rho _{mp}\kappa _{nt}\kappa _{qs}^{}+2\left(z_p^{}z_m\right)\kappa _{nt}\kappa _{qs}^{}+2\rho _{mp}\left(z_nz_t\right)\kappa _{qs}^{}+2\rho _{mp}\kappa _{nt}\left(z_q^{}z_s^{}\right)\hfill \end{array}\right\}`$
$`\left(I_\kappa \right)_{ji}`$ $`=2{\displaystyle \underset{smn}{}}{\displaystyle \underset{pqlt}{}}V_{jsmn}\stackrel{~}{V}_{pqlt}\left(\delta _{li}\right)`$ (102)
$`\times \left\{\begin{array}{c}2\rho _{mp}\kappa _{nt}\left(\rho _{ls}+\delta _{ls}\right)+2\rho _{mp}\kappa _{nt}\left(z_s^{}z_l\right)+2\rho _{mp}\left(z_nz_t\right)\left(\rho _{ls}+\delta _{ls}\right)\hfill \\ +2\left(z_p^{}z_m\right)\kappa _{nt}\left(\rho _{ls}+\delta _{ls}\right)+\kappa _{ml}\kappa _{nt}\kappa _{ps}^{}+2\kappa _{ml}\left(z_nz_t\right)\kappa _{ps}^{}+\kappa _{ml}\kappa _{nt}\left(z_p^{}z_s^{}\right)\hfill \end{array}\right\}`$
Formal integration of the equations of motion (70)-(75) (taking into account their suppressed ‘free’ evolution in the HFB basis) and substitution into equations (51)-(52) generates all second order collisional processes, as discussed in more detail in the next section. Before proceeding, we give at this point the most compact form of our equations to second order as

$`i\left({\displaystyle \frac{dR_c}{dt}}\right)=\left({\displaystyle \frac{\widehat{E}}{R_c}}\right)R_c`$ (103)
$`i\left({\displaystyle \frac{dR_e}{dt}}\right)_{ji}=\left[{\displaystyle \underset{l}{}}\left({\displaystyle \frac{\widehat{E}}{(R_e)_{lj}}}\right)\left(R_e\right)_{li}i\left(I_{ji}\right)\right]h.c.`$ (104)
where the matrix $`I`$ is defined as

$`I=\left(\begin{array}{c}I_\rho \\ I_\kappa ^{}\end{array}\begin{array}{c}I_\kappa \\ I_\rho ^{}\end{array}\right)`$ (109)
## 5 Link to other kinetic theories

Having discussed the extent to which a variational methodology can be useful, let us now show explicitly how our equations reproduce those of Walser et al. , as written in generalised matrix form in their follow-up paper . By thinking in terms of forward and backward scattering rates incorporating terms $`I_\rho `$ and $`I_\kappa `$, the equations of motion for the non-equilibrium matrix $`K`$ can be straightforwardly re-expressed as

$`i{\displaystyle \frac{dK}{dt}}=`$ $`\left(\begin{array}{c}\left(\mathrm{\Gamma }_\rho I_\rho \right)\\ I_\kappa ^{}\end{array}\begin{array}{c}\left(\mathrm{\Gamma }_\kappa +I_\kappa \right)\\ I_\rho ^{}\end{array}\right)\left(\begin{array}{c}\rho \\ \kappa ^{}\end{array}\begin{array}{c}\kappa \\ \left(\rho ^{}+\text{1}\text{1}\right)\end{array}\right)`$ (118)
$``$ $`\left(\begin{array}{c}I_\rho \\ \left[\mathrm{\Gamma }_\kappa +I_\kappa \right]^{}\end{array}\begin{array}{c}I_\kappa \\ \left[\mathrm{\Gamma }_\rho I_\rho \right]^{}\end{array}\right)\left(\begin{array}{c}\left(\rho +\text{1}\text{1}\right)\\ \kappa ^{}\end{array}\begin{array}{c}\kappa \\ \rho ^{}\end{array}\right)+h.c.`$ (127)
Following the notation of Walser et al. , we now define forward and backward collision operators as

$`\mathrm{\Gamma }^<=\left(\begin{array}{c}I_\rho \\ \left[\mathrm{\Gamma }_\kappa +I_\kappa \right]^{}\end{array}\begin{array}{c}I_\kappa \\ \left[\mathrm{\Gamma }_\rho I_\rho \right]^{}\end{array}\right)`$ (132)
and $`\mathrm{\Gamma }^>=\sigma _1\left(\mathrm{\Gamma }^<\right)^{}\sigma _1`$. The corresponding generalised quasiparticle self-energies are $`R_e^>=R_e`$ as defined in equation (23) and $`R_e^<=\sigma _1\left(R_e^>\right)^{}\sigma _1`$. We thus obtain

$`i{\displaystyle \frac{dK}{dt}}=[(\mathrm{\Gamma }^<R_e^<\mathrm{\Gamma }^>R_e^>)+h.c.]`$ (133)
yielding for the quasiparticle propagator the following second order expression

$`{\displaystyle \frac{dR_e^>}{dt}}=iH_eR_e^>+\left(\mathrm{\Gamma }^<R_e^<\mathrm{\Gamma }^>R_e^>\right)+h.c.`$ (134)
Careful examination of the above equation shows clearly that our approach generates exactly the same forward and backward collision operators for the non-condensate propagator, as discussed in .

Although our methodology generates the forward and backward collision operators for the non-condensate propagator in a rather straightforward manner, the situation becomes somewhat more complicated when considering the evolution of the condensed component. The reason is that our methodology does not yield such expressions directly, but instead it provides us with the sum of forward and backward rates, via

$`\left({\displaystyle \frac{dR_c}{dt}}\right)=i\left({\displaystyle \frac{\widehat{E}}{R_c}}\right)R_c=i\left({\displaystyle \frac{\widehat{E}_{HFB}^{(1)}}{R_c}}+{\displaystyle \frac{\widehat{E}^{(2)}}{R_c}}\right)R_c=iH_cR_c+YR_c`$ (135)
where

$`Y=\left(\begin{array}{c}\mathrm{\Gamma }_z\\ \left[\mathrm{\Gamma }_z^{}\right]^{}\end{array}\begin{array}{c}\mathrm{\Gamma }_z^{}\\ \left[\mathrm{\Gamma }_z\right]^{}\end{array}\right)`$ (140)
thus leaving their particular definition somewhat arbitrary. A natural definition corresponds to identifying these rates in terms of positive and negative contributions of $`\mathrm{\Gamma }_{z^{()}}`$ (denoted respectively by $`\pm ve[\mathrm{\Gamma }_{z^{()}}]`$), thus defining the forward and backward scattering rates for the condensate propagator as

$`Y^<=\left(\begin{array}{c}ve\left[\mathrm{\Gamma }_z\right]\\ \left(+ve\left[\mathrm{\Gamma }_z^{}\right]\right)^{}\end{array}\begin{array}{c}ve\left[\mathrm{\Gamma }_z^{}\right]\\ \left(+ve\left[\mathrm{\Gamma }_z\right]\right)^{}\end{array}\right)`$ (145)
and

$`Y^>=\sigma _1\left(Y^<\right)^{}\sigma _1=\left(\begin{array}{c}+ve\left[\mathrm{\Gamma }_z\right]\\ \left(ve\left[\mathrm{\Gamma }_z^{}\right]\right)^{}\end{array}\begin{array}{c}+ve\left[\mathrm{\Gamma }_z^{}\right]\\ \left(ve\left[\mathrm{\Gamma }_z\right]\right)^{}\end{array}\right)`$ (150)
which clearly satisfy $`Y=\left(Y^<Y^>\right)`$.

In comparing our final expressions to those of Walser et al. , we note that our combined evolution equations agree entirely with the first formulation of the theory of Walser et al. . In their subsequent work , where they have explicitly identified forward and backward rates as described in this section, their expressions for condensate rates differ from those we have given above in that their corresponding expression for $`\mathrm{\Gamma }_z^{}`$ of equation (83) is

$`\left(\mathrm{\Gamma }_z^{}\right)_{iq}=2{\displaystyle \underset{smn}{}}{\displaystyle \underset{pqlt}{}}V_{jsmn}\stackrel{~}{V}_{pqlt}\left\{\begin{array}{c}2\left[\left(\rho _{mp}+\delta _{mp}\right)\kappa _{nt}\rho _{ls}\rho _{mp}\kappa _{nt}\left(\rho _{ls}+\delta _{ls}\right)\right]\hfill \\ +\left[\kappa _{ml}\kappa _{nt}\kappa _{ps}^{}\kappa _{ml}\kappa _{nt}\kappa _{ps}^{}\right]\hfill \end{array}\right\}`$ (153)
with the second line of the above expression, which identically cancels itself, being generated additionally to what appears in our expression of equation (83). The appearance of such a contribution is physically appealing, since it implies that the collision processes occuring in condensate and non-condensate are of the same basic structure , indicating that processes $`\mathrm{\Gamma }^>`$ can be generated from $`Y^>`$ by functional differentiation, as done in Green functions’ techniques . Although the net contribution of such a term is zero, its presence modifies the physical expressions of forward and backward collision integrals, i.e. it adds a particular collisional diagram to each. In our treatment, it is clear that such an extra ‘self-cancelling’ contribution cannot be generated in the first place. This is because the adiabatic elimination procedure yields for the second order expression for $`z`$ terms proportional, at most, to correlations of five operators, whereas these additional contributions would require correlations of seven single-particle operators, and hence a treatment of three-particle interactions , which extends beyond the binary hamiltonian employed here. It is therefore somewhat perplexing that such terms are found in the second formulation of the theory of Walser et al. . However, if one is looking for a theory which satisfies such symmetries as mentioned above, due to specific physical considerations, we can indeed add and subtract these (or other similar) contributions to the forward and backward condensate collision integrals, obtaining exactly the same contributions as those discussed by Walser et al. in  and by Wachter et al. . Such treatment will, however, not be a rigorous derivation, since its last stage will rely on additional physical arguments (e.g. gapless or conserving ) arising from alternative approaches.

We can now summarize the entire findings of our perturbative HFB treatment in the following systems of equations

$`\left\{\begin{array}{c}i\frac{dR_c}{dt}=H_cR_c+J\\ i\frac{dR_e^>}{dt}=H_eR_e^>R_e^>H_e^{}+K\end{array}\right\}`$ (156)
$`\left\{\begin{array}{c}i\frac{dJ}{dt}=\left(Y^<Y^>\right)R_c\\ i\frac{dK}{dt}=[(\mathrm{\Gamma }^<R_e^<\mathrm{\Gamma }^>R_e^>)+h.c.]\end{array}\right\}`$ (159)
Combining these results yields for the propagators of the system the following second order evolution

$`\left\{\begin{array}{c}\frac{dR_c}{dt}=iH_cR_c+\left(Y^<Y^>\right)R_c\\ \frac{dR_e^>}{dt}=iH_eR_e^>+\left(\mathrm{\Gamma }^<R_e^<\mathrm{\Gamma }^>R_e^>\right)+h.c.\end{array}\right\}`$ (162)
## 6 Conclusions

By applying second order perturbation theory from an HFB basis and making no further approximations on the binary interaction hamiltonian, we have generalised an earlier quantum kinetic approach to Bose-Einstein condensation presented by the author and Burnett elsewhere , by treating in a self-consistent manner the dynamics of both condensed and non-condensed components of the system. This involves a detailed consideration of the evolution of correlations of up to four shifted single-particle operators. Our methodology can be summarised as follows: Firstly we obtain the exact equations of motion for the condensate and non-condensate propagators in the binary interaction hamiltonian. These depend on matrices expressed in terms of higher order correlations, which can be approximately treated by separating off their ‘equilibrium’ and ‘non-equilibrium’ contributions (somewhat analogously to ). The ‘equilibrium’ contributions are decomposed into lower order averages in terms of the HFB order parameters (thus yielding the usual first order evolution), whereas the non-equilibrium matrices contain all contributions of order higher than those obtained from the HFB energy functional. In this paper, we have focused on the lowest order expression for the evolution of such non-equilibrium matrices. Substitution of these results into the exact equations of motion for the condensate and non-condensate propagators yields their respective evolution to second order in the effective interatomic potentials and these are shown to be identical to those of Walser et al. . Such an equivalence was of course anticipated, since they are both based on essentially the same underlying assumption of slowly-varying quantities. In the work of Walser et al. , these quantities correspond to the ‘master’ variables determining the evolution of the gas on a coarse-grained timescale, whereas in the treatment presented here, these quantities are implicit in our choice of the HFB hamiltonian as the unperturbed hamiltonian for applying perturbation theory. Our approach has provided the exact form of the energy functional (equations (18)-(20)), which readily enables us to extend this treatement to higher orders. This may be useful, for example, to test a key underlying assumption of all such second order theories for dilute, weakly-interacting Bose-condensed gases, that the higher order terms (in the perturbative expansion) are significantly smaller. This should be valid sufficiently far from the critical region, when $`(kT/nU_o)(\sqrt{na^3})1`$  where $`a`$ is the s-wave scattering length, $`n`$ the condensate density and $`U_o`$ the effective two-body interatomic potential (and can thus be justified by the fact that the problem has been formulated in terms of an unperturbed quasiparticle basis and that the interactions have been treated by an effective resummed bare t-matrix).

Furthermore, our work is based on the same idea as the second order theory of excitations developed by Morgan  and recently applied numerically by the Oxford group , thus showing the link between these two approaches. We should also point out that our approach is formally connected to the detailed kinetic treatment of Zaremba-Nikuni-Griffin  (whose analysis focuses on the hydrodynamic regime). Moreover, in a recent paper , the theory of Walser et al.  has been shown to be equivalent to the Green’s function methodology originally developed by Kadanoff and Baym  and recently applied to the study of trapped condensates by Imamovic-Tomasovic and Griffin . This latter approach is in turn formally related to the Schwinger-Keldysh method for treating non-equilibrium many-body systems, which has been employed by Stoof . We stress that the work we have presented in this paper includes all relevant quasiparticle effects, despite being expressed in terms of single-particle operators; to generalise this, one could explicitly transform to quasiparticle operators. Carrying out such a transformation within the framework of the kinetic theory of Walser et al. , Wachter has shown  that the collisional integrals simplify enormously, acquiring the usual form of Boltzmann-like factors, thus enabling a direct physical interpretation (which is rather concealed in the lengthy expressions in terms of single-particle operators). Similar work has been performed by Imamovic-Tomasovic and Griffin who recently described such a quasiparticle kinetic equation . Here we should further point out that the static theory of Morgan  was formulated explicitly in terms of quasiparticles, and so where the related linear response discussions of Rusch et al.  and Giorgini . The links of the current approach independently to those of Walser et al. (and hence Imamovic-Tomasovic and Griffin) and Morgan suggest that these time-dependent theories may correspond to some form of time-dependent generalisation of the gapless excitation theory of Morgan, although a more direct link remains to be established. Based on these considerations we hence believe that this paper represents yet another contribution towards the long-sought goal of a universally accepted kinetic theory for dilute, weakly-interacting systems exhibiting Bose-Einstein condensation.

I am grateful to Peter Lambropoulos for detailed discussions. I would also like to thank Reinhold Walser for in-depth discussions during the First International Workshop on the Theory of Quantum Gases and Quantum Coherence held in Salerno, Italy, and Keith Burnett for various general discussions on quantum kinetic theory during the past few years. This work was partially supported by the ULF (Contract No. ERB-FMGECT 950021).

## Appendix A Equations of motion for multiple shifted single-particle operator correlations

In this appendix we give expressions for the equations of motion of correlations of up to four shifted single-particle operators. Our expressions are only explicit for the effect of the reduced hamiltonian $`\left(\widehat{H}_3^{^{}}+\widehat{H}_4^{^{}}\right)`$, with the effect of the HFB hamiltonian $`\left(H_0+\widehat{H}_1^{^{}}+\widehat{H}_2^{^{}}\right)`$ being implicit in the renormalised eigenenergies $`\stackrel{~}{\omega }_i`$ appearing below (and also in the approximately energy-conserving matrix elements $`\stackrel{~}{V}_{pqlt}`$ of equation (30)). Explicit expressions in terms of the entire hamiltonian of the system $`\widehat{H}`$ of equation (1) (including also the ones needed to derive the first order HFB-basis contributions to the coupled dynamics, i.e. first terms in each of equations (51)-(52)) can be found in . We start by noting that the HFB quantities $`z,\rho ,\kappa `$ will, to lowest order in the potential, have no contribution in the hamiltonian $`\left(\widehat{H}_3^{^{}}+\widehat{H}_4^{^{}}\right)`$, since the quantities found on their right-hand-sides acquire no equilibrium mean value within the HFB approximation. In particular, we find (making use of symmetries wite respect to summed indices)

$$i\frac{d}{dt}\left(z_i\right)=\stackrel{~}{\omega }_iz_i+\underset{smn}{}V_{ismn}\lambda _{smn}$$
(163)
$`i{\displaystyle \frac{d}{dt}}\left(\rho _{ji}\right)`$ $`=\left(\stackrel{~}{\omega }_j\stackrel{~}{\omega }_i\right)\rho _{ji}+{\displaystyle \underset{smn}{}}V_{jsmn}\left[\sigma _{ismn}+\left(2\lambda _{mis}^{}z_n+\lambda _{imn}z_s^{}\right)\right]`$ (164)
$`{\displaystyle \underset{smn}{}}V_{mnsi}\left[\sigma _{mnsj}+\left(2\lambda _{mjs}z_n^{}+\lambda _{jmn}^{}z_s\right)\right]`$
$`i{\displaystyle \frac{d}{dt}}\left(\kappa _{ji}\right)`$ $`=\left(\stackrel{~}{\omega }_i+\stackrel{~}{\omega }_j\right)\kappa _{ji}+{\displaystyle \underset{smn}{}}V_{jsmn}\left[\xi _{simn}+\left(2\lambda _{sim}z_n+\gamma _{imn}z_s^{}\right)\right]`$ (165)
$`+{\displaystyle \underset{smn}{}}V_{ismn}\left[\xi _{sjmn}+\left(2\lambda _{sjm}z_n+\gamma _{jmn}z_s^{}\right)\right]`$
where the entire first order evolution of condensate / non-condensate propagators has now been summed up in the renormalised HFB eigenenergies $`\stackrel{~}{\omega }_i`$.

For the non-HFB quantities, we find the following evolution

$`i{\displaystyle \frac{d}{dt}}\left(\lambda _{smn}\right)`$ $`=\left(\stackrel{~}{\omega }_m+\stackrel{~}{\omega }_n\stackrel{~}{\omega }_s\right)\lambda _{smn}+{\displaystyle \underset{lt}{}}V_{mnlt}\left(2\rho _{ls}z_t\right)`$ (166)
$`+{\displaystyle \underset{plt}{}}V_{pnlt}\left[2\rho _{mp}\rho _{ls}z_t+2\kappa _{ml}\kappa _{ps}^{}z_t+2\rho _{ls}\kappa _{mt}z_p^{}\right]`$
$`+{\displaystyle \underset{plt}{}}V_{pmlt}\left[2\rho _{np}\rho _{ls}z_t+2\kappa _{nl}\kappa _{ps}^{}z_t+2\rho _{ls}\kappa _{nt}z_p^{}\right]`$
$`{\displaystyle \underset{pql}{}}V_{pqls}\left[2\rho _{mp}\rho _{nq}z_l+2\left(\rho _{mp}\kappa _{nl}+\rho _{np}\kappa _{ml}\right)z_q^{}\right]`$
$`i{\displaystyle \frac{d}{dt}}\left(\gamma _{imn}\right)`$ $`=\left(\stackrel{~}{\omega }_i+\stackrel{~}{\omega }_m+\stackrel{~}{\omega }_n\right)\gamma _{imn}`$ (167)
$`+2{\displaystyle \underset{lt}{}}\left[V_{mnlt}\left(\kappa _{il}z_t\right)+V_{inlt}\left(\kappa _{ml}z_t\right)+V_{imlt}\left(\kappa _{nl}z_t\right)\right]`$
$`+2{\displaystyle \underset{plt}{}}V_{pnlt}\left[\rho _{ip}\kappa _{ml}z_t+\rho _{mp}\kappa _{il}z_t+\kappa _{il}\kappa _{mt}z_p^{}\right]`$
$`+2{\displaystyle \underset{plt}{}}V_{pmlt}\left[\rho _{ip}\kappa _{nl}z_t+\rho _{np}\kappa _{il}z_t+\kappa _{il}\kappa _{nt}z_p^{}\right]`$
$`+2{\displaystyle \underset{plt}{}}V_{pilt}\left[\rho _{mp}\kappa _{nl}z_t+\rho _{np}\kappa _{ml}z_t+2\kappa _{ml}\kappa _{nt}z_p^{}\right]`$
$`i{\displaystyle \frac{d}{dt}}\left(\sigma _{ismn}\right)`$ $`=\left(\stackrel{~}{\omega }_m+\stackrel{~}{\omega }_n\stackrel{~}{\omega }_i\stackrel{~}{\omega }_s\right)\sigma _{ismn}`$ (168)
$`+{\displaystyle \underset{lt}{}}V_{mnlt}\left(2\rho _{li}\rho _{ts}\right){\displaystyle \underset{pq}{}}V_{pqis}\left(2\rho _{mp}\rho _{nq}\right)`$
$`+2{\displaystyle \underset{plt}{}}V_{pnlt}\left(\rho _{mp}\rho _{li}\rho _{ts}+\rho _{ls}\kappa _{mt}\kappa _{pi}^{}+\rho _{li}\kappa _{mt}\kappa _{ps}^{}\right)`$
$`+2{\displaystyle \underset{plt}{}}V_{pmlt}\left(\rho _{np}\rho _{li}\rho _{ts}+\rho _{ls}\kappa _{nt}\kappa _{pi}^{}+\rho _{li}\kappa _{nt}\kappa _{ps}^{}\right)`$
$`2{\displaystyle \underset{pql}{}}V_{pqli}\left(\rho _{mq}\rho _{np}\rho _{ls}+2\rho _{mp}\kappa _{nl}\kappa _{qs}^{}\right)`$
$`2{\displaystyle \underset{pql}{}}V_{pqls}\left(\rho _{mq}\rho _{np}\rho _{li}+2\rho _{mp}\kappa _{nl}\kappa _{qi}^{}\right)`$
$`i{\displaystyle \frac{d}{dt}}\left(\xi _{simn}\right)`$ $`=\left(\stackrel{~}{\omega }_i+\stackrel{~}{\omega }_m+\stackrel{~}{\omega }_n\stackrel{~}{\omega }_s\right)\xi _{simn}`$ (169)
$`+{\displaystyle \underset{lt}{}}\left[V_{mnlt}\left(2\rho _{ls}\kappa _{it}\right)+V_{inlt}\left(2\rho _{ls}\kappa _{mt}\right)+V_{imlt}\left(2\rho _{ls}\kappa _{nt}\right)\right]`$
$`+2{\displaystyle \underset{plt}{}}V_{pnlt}\left[\kappa _{ml}\kappa _{it}\kappa _{ps}^{}+\rho _{mp}\rho _{ls}\kappa _{it}+\rho _{ip}\rho _{ls}\kappa _{mt}\right]`$
$`+2{\displaystyle \underset{plt}{}}V_{pmlt}\left[\kappa _{nl}\kappa _{it}\kappa _{ps}^{}+\rho _{np}\rho _{ls}\kappa _{it}+\rho _{ip}\rho _{ls}\kappa _{nt}\right]`$
$`+2{\displaystyle \underset{plt}{}}V_{pilt}\left[\kappa _{nl}\kappa _{mt}\kappa _{ps}^{}+\rho _{np}\rho _{ls}\kappa _{mt}+\rho _{mp}\rho _{ls}\kappa _{nt}\right]`$
$`2{\displaystyle \underset{pql}{}}V_{pqls}\left[\rho _{mp}\rho _{nq}\kappa _{li}+\rho _{mp}\rho _{iq}\kappa _{nl}+\rho _{ip}\rho _{nq}\kappa _{ml}\right]`$
We stress that equations (166)-(169) given above are not exact, even within our consistent decoupling approximations, because we have systematically ignored from their right hand sides contributions which depend on ‘non-equilibrium’ correlations of more than two shifted single-particle operators (i.e. $`\lambda `$, $`\gamma `$, $`\sigma `$, or $`\xi `$). Such contributions will ultimately yield the evolution of the generalised condensate / quasiparticle propagators beyond second order in the potential, in the same manner that equations (163)-(165) will also yield contributions to third or higher orders. We believe this is fully equivalent to the statement made in Appendix A of Walser et al. , that terms leading to higher order contributions are consistently neglected. We could, of course, have very straightforwardly included such terms in the equations of motion (166)-(169). However, since we do not deal with such terms any further, this would only introduce unnecessary complexity.

## Appendix B Decoupling approximations

As mentioned in the text, correlations of even numbers of shifted single-particle operators, can be thought of as consisting of two contributions , one part being due to the mean value at equilibrium when decomposed by Wick’s theorem into products of lower averages (here referring to binary averages) and the remaining term yielding the rate of change of this quantity from its equilibrium value. For example $`\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n=\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n^0+\sigma _{rsmn}`$ with the ‘non-equilibrium’ quantities (here $`\sigma `$) treated perturbatively, by means of their respective equations of motion given in Appendix A.

Since the HFB basis, by definition, does not allow for nonzero equilibrium values of triplet correlations, we note that, by analogy, any correlations of odd products of shifted single-particle operators will vanish to lowest order (since their final decomposed forms will always depend on such triplet products as discussed in ) and will yield non-vanishing contributions only at the next order in the potential, via their respective equations of motions discussed in Appendix A. The decoupling approximations we have employed for the ‘equilibrium’ quantities are as follows

$$\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n^0=\rho _{mr}\rho _{ns}+\rho _{nr}\rho _{ms}+\kappa _{mn}\kappa _{rs}^{}$$
(170)
At equilibrium, this would arise directly from Wick’s theorem, which is the main motivation for such a decoupling approximation here. The remaining decorrelations are performed analogously, namely

$$\widehat{c}_p^{}\widehat{c}_q\widehat{c}_l\widehat{c}_t=\rho _{qp}\kappa _{lt}+\rho _{lp}\kappa _{qt}+\rho _{tp}\kappa _{ql}$$
(171)
$$\widehat{c}_p\widehat{c}_q\widehat{c}_l\widehat{c}_t=\kappa _{pq}\kappa _{lt}+\kappa _{pl}\kappa _{qt}+\kappa _{pt}\kappa _{ql}$$
(172)
$`\widehat{c}_p^{}\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_q\widehat{c}_l\widehat{c}_t`$
$`=\left\{\begin{array}{c}\rho _{qp}\left(\rho _{lr}\rho _{ts}+\rho _{tr}\rho _{ls}\right)+\rho _{lp}\left(\rho _{qr}\rho _{ts}+\rho _{tr}\rho _{qs}\right)+\rho _{tp}\left(\rho _{qr}\rho _{ls}+\rho _{lr}\rho _{qs}\right)\hfill \\ +\kappa _{ql}\left(\kappa _{pr}^{}\rho _{ts}+\kappa _{ps}^{}\rho _{tr}+\kappa _{rs}^{}\rho _{tp}\right)+\kappa _{qt}\left(\kappa _{pr}^{}\rho _{ls}+\kappa _{ps}^{}\rho _{lr}+\kappa _{rs}^{}\rho _{lp}\right)\hfill \\ +\kappa _{lt}\left(\kappa _{pr}^{}\rho _{qs}+\kappa _{ps}^{}\rho _{qr}+\kappa _{rs}^{}\rho _{qp}\right)\hfill \end{array}\right\}`$ (176)
$`\widehat{c}_p^{}\widehat{c}_q^{}\widehat{c}_m\widehat{c}_r\widehat{c}_l\widehat{c}_t`$
$`=\left\{\begin{array}{c}\rho _{mp}\left(\rho _{rq}\kappa _{lt}+\rho _{lq}\kappa _{rt}+\rho _{tq}\kappa _{rl}\right)+\rho _{rp}\left(\rho _{mq}\kappa _{lt}+\rho _{lq}\kappa _{mt}+\rho _{tq}\kappa _{ml}\right)\hfill \\ +\rho _{lp}\left(\rho _{rq}\kappa _{mt}+\rho _{mq}\kappa _{rt}+\rho _{tq}\kappa _{mr}\right)+\rho _{tp}\left(\rho _{rq}\kappa _{lm}+\rho _{mq}\kappa _{rl}+\rho _{lq}\kappa _{mr}\right)\hfill \\ +\kappa _{pq}^{}\left(\kappa _{rm}\kappa _{lt}+\kappa _{lm}\kappa _{rt}+\kappa _{tm}\kappa _{rl}\right)\hfill \end{array}\right\}`$ (180)
## References