text:
# Multi-Instanton Calculus in Supersymmetric Theories

## 1 Introduction

Our understanding of the non–perturbative sector of field and string theories has greatly progressed in recent times. Recently, for the first time, the entire non–perturbative contribution to the holomorphic part of the Wilsonian effective action was computed for $`N=2`$ globally supersymmetric (SUSY) theories with gauge group $`SU(2)`$, using ansätze dictated by physical intuitions . A few years later, a better understanding of non–perturbative configurations in string theory led to the conjecture that certain IIB string theory correlators on an $`AdS_5\times S^5`$ background are related to Green’s functions of composite operators of an $`N=4`$ $`SU(N_c)`$ Super Yang–Mills (SYM) theory in four dimensions in the large $`N_c`$ limit . Although supported by many arguments, these remarkable results remain conjectures and a clear mathematical proof seems to be out of reach at the moment. In our opinion this state of affairs is mainly due to the lack of adequate computational tools in the non–perturbative region. To the extent of our knowledge, the only way to perform computations in this regime in SUSY theories and from first principles is via multi–instanton calculus. Using this tool, many partial checks have been performed on these conjectures, both in $`N=2`$ and $`N=4`$ SUSY gauge theories , , , , . The limits on these computations come from the exploding amount of algebraic manipulations to be performed and from the lack of an explicit parametrization of instantons of winding number greater than two . In order to develop new computational tools that might allow an extension to arbitrary winding number, I revisit instanton computations for $`N=2`$ in the light of the topological theory built out of $`N=2`$ SYM , i.e. the so–called Topological Yang–Mills theory (TYM) .

## 2 Topological Yang–Mills Theory

Here I collect some results which will be relevant to our discussion. I use the same notation of  to which I refer the reader for a detailed exposition of this material. As it is well known , after the twisting procedure, the Lagrangian of $`N=2`$ SYM is invariant under

$`sA=\psi Dc,`$
$`s\psi =[c,\psi ]D\varphi ,`$
$`s\varphi =[c,\varphi ],`$
$`sc={\displaystyle \frac{1}{2}}[c,c]+\varphi .`$ (1)
The BRST operator, $`s`$, defined in (2) is such that $`s^2=0`$ and when the set of equations in (2) is restricted to the solutions of the Euler–Lagrange classical equation (the zero modes) it gives the derivative on the space $`M^+`$. In terms of the parameters of the ADHM construction, these solutions look like

$`A`$ $`=`$ $`U^{}dU,`$ (2)
$`\psi `$ $`=`$ $`U^{}f(d\mathrm{\Delta }^{})U+U^{}(d\mathrm{\Delta })f^{}U,`$ (3)
$`c`$ $`=`$ $`U^{}(s+C)U,`$ (4)
$`\varphi `$ $`=`$ $`U^{}f^{}U+U^{}𝒜U.`$ (5)
Plugging (5) into (2) leads to the action of the operator $`s`$ on the elements of the ADHM construction

$``$ $`=`$ $`s\mathrm{\Delta }+C\mathrm{\Delta }=𝒮\mathrm{\Delta },`$ (6)
$`𝒜`$ $`=`$ $`s\mathrm{\Delta }+C=𝒮,`$ (7)
$`s𝒜`$ $`=`$ $`[C,𝒜],`$ (8)
$`sC`$ $`=`$ $`𝒜CC,`$ (9)
i.e. this is the realization of the BRST algebra on the instanton moduli space. $`C`$ is a connection I must introduce to have a nilpotent $`s`$.

Let us see now how, at the semi–classical level, any correlator which is expressed as a polynomial in the fields, becomes after projection onto the zero–mode subspace, a well–defined differential form on $`M^+`$ . Symbolically

$$fields=_{M^+}\left[(fields)e^{S_{\mathrm{TYM}}}\right]_{zeromode}.$$
(10)
A generic function on the zero–mode subspace will then have the expansion

$`g(\widehat{\mathrm{\Delta }},\widehat{})`$ $`=`$ $`g_0(\widehat{\mathrm{\Delta }})+g_{i_1}(\widehat{\mathrm{\Delta }})\widehat{}_{i_1}+\mathrm{}`$ (11)
$`+`$ $`{\displaystyle \frac{1}{p!}}g_{i_1i_2\mathrm{}i_p}(\widehat{\mathrm{\Delta }})\widehat{}_{i_1}\widehat{}_{i_2}\mathrm{}\widehat{}_{i_p},`$
the coefficients of the expansion being totally antisymmetric in their indices. Now $`\widehat{}_i`$’s and the $`s\widehat{\mathrm{\Delta }}_i`$’s are related by a (moduli–dependent) linear transformation $`K_{ij}`$, which is completely known once the explicit expression for $`C`$ is plugged into the $`\widehat{}_i`$’s:

$$\widehat{}_i=K_{ij}(\widehat{\mathrm{\Delta }})s\widehat{\mathrm{\Delta }}_j.$$
(12)
It then follows that

$`\widehat{}_{i_1}\widehat{}_{i_2}\mathrm{}\widehat{}_{i_p}=`$ (13)
$`=`$ $`ϵ_{j_1\mathrm{}j_p}K_{i_1j_1}K_{i_2j_2}\mathrm{}K_{i_pj_p}s^p\widehat{\mathrm{\Delta }}=`$
$`=`$ $`ϵ_{i_1\mathrm{}i_p}(\mathrm{det}K)s^p\widehat{\mathrm{\Delta }},`$
where $`s^p\widehat{\mathrm{\Delta }}s\widehat{\mathrm{\Delta }}_1\mathrm{}s\widehat{\mathrm{\Delta }}_p`$. I then conclude that

$`{\displaystyle _{M^+}}g(\widehat{\mathrm{\Delta }},\widehat{})={\displaystyle \frac{1}{p!}}{\displaystyle _{M^+}}g_{i_1\mathrm{}i_p}(\widehat{\mathrm{\Delta }})\widehat{}_{i_1}\mathrm{}\widehat{}_{i_p}`$ (14)
$`=`$ $`{\displaystyle _{M^+}}s^p\widehat{\mathrm{\Delta }}|\mathrm{det}K|g_{12\mathrm{}p}(\widehat{\mathrm{\Delta }}).`$
The determinant of $`K`$ naturally stands out as the instanton integration measure for $`N=2`$ SYM theories.

## 3 Conclusions

In this talk i have argued that the results of multi-instanton calculus at the semiclassical level can be easily recovered in the formalism of topological field theories. More benefits come from this reformulation than what I have presented until now. It is for example possible to show that correlators of the type of (10) can be written as a total derivative on the moduli space. The only contributions to these quantities, come from zero-size instantons. Given the peculiar properties of the ADHM construction at the boundary of the moduli space, this might lead to recursion relations among correlators computed at different winding number. To get to this conclusion, we probably have to better understand what is the geometrical role of the action, projected on the subspace of the zero-modes. The connection between moduli spaces of instantons and their construction in terms of D-branes of string theory can be very helpful in this respect.


format: text

text:
# Basis-dependent dynamics of trapped Bose-Einstein condensates and analogies with semi-classical laser theory

## I Introduction

The kinetic theories of Bose-Einstein condensation in a trap can be divided into two general categories. The fully quantum approaches have been based either on a perturbative master equation treatment for the many-body density matrix , or on a single Fokker-Planck equation for the nonequilibrium dynamics of the entire system . Alternative approaches to the nonequilibrium dynamics of trapped Bose-condensed gases are essentially of a perturbative nature, based on a suitably truncated coupled equation of motion hierarchy for normal and anomalous averages . The growth of condensation was first studied numerically by Gardiner and co-workers , based on the theory developed in , and their findings were in qualitative agreement with the description of Kagan, Svistunov and Shlyapnikov . Independent studies by Stoof and co-workers, based on similar approximations, have produced growth curves  which are in very good agreement with those of Gardiner et al. .

The large condensates typically produced in experiments lead to large mean field potentials which significantly dress the single-particle eigenenergies of the trap potential, thus making it convenient to work in a suitably dressed basis. The conventional theoretical picture adopted is thus that of a very large single-mode condensate interacting with a large number of quasiparticles and higher-lying excited atoms, whereas it is not a priori necessary that single-mode condensation will correspond to all experimental conditions . Instead of working in the usual ‘condensed matter’ approach of a condensate and a set of quasiparticle excitations which are, by definition, orthogonal to it (which is suitable for large condensates) we focus our description in terms of a more ‘quantum-optical’ approach, i.e. in terms of individual modes of the trapping potential , for which one must in principle consider the condensate spanning a large number of modes, with the same modes being simultaneously occupied by non-condensate atoms. Within such a picture, this paper explicitly discusses the simplest deviation from single-mode condensation in the context of a bare single-particle basis, in a manner analogous to existing treatments. Formulating the problem in terms of bare particle eigenenergies can only be useful in the case of extremely small, dilute, weakly-interacting condensates, when the mean field effects are small enough that their induced shifts of the single-particle energies can be treated as perturbations. Such an approach is nevertheless beneficial for discussing the inherent multi-mode nature of trapped condensates and their relation to the corresponding (semi-classical) theory of multi-mode photon lasers.

The first part of this paper (Sec. II) reviews in the usual manner the formal development of the coupled equation of motion approach. To simplify the physical picture and bring out the underlying structure in a clear manner, we focus on a very simple system consisting of the three lowest trap eigenstates, for which we derive explicit equations of motion in the Popov approximation  in Sec. III (with corresponding off-equilibrium contributions given in Appendix A, and some further clarifications in Appendix B). By explicitly discussing the interplay of two coupled condensed modes, we show how our treatment reduces to the Hartree-Fock theory for binary condensates  (Sec. IV A1). An important advantage of formulating the theory in terms of bare single-particle eigenstates is that it allows us to discuss in detail the analogy between our equations for the coupled dynamics of two condensed modes, and the corresponding ones arising in the semi-classical treatments of two-mode photon lasers (Sec. IV A2). We show that, the extent to which this analogy can be drawn for the inherent multi-mode nature of a single condensate, depends critically on the choice of the single-particle basis (i.e. bare or dressed by various types of mean fields) in terms of which the analysis is carried out. Our discussion here is distinct from conventional analogies based on two different condensates which are spatially separated , in different spin states , or outcoupled by radiation applied at two different frequencies . Sec. IV B compares our approach to conventional kinetic treatments, where we show explicitly that our theory reduces to the multi-mode kinetic treatment of Walser et al., upon shifting our description to a basis in which the single-particle eigenenergies become dressed by the usual Hartree-Fock-Bogoliubov (HFB) mean field potentials. This suggests that, contrary to an implication in , such dynamic equations (and collisional integrals) are inevitably basis-dependent.

## II The Coupled Equation of Motion Approach

Consider a sufficiently dilute, weakly-interacting partially Bose-condensed trapped gas with a binary-interaction hamiltonian

$$\widehat{H}=\underset{rs}{}\mathrm{\Xi }_{rs}^{Bare}\widehat{a}_r^{}\widehat{a}_s+\frac{1}{2}\underset{rsmn}{}V_{rsmn}\widehat{a}_r^{}\widehat{a}_s^{}\widehat{a}_m\widehat{a}_n$$
(1)
Here $`\widehat{\mathrm{\Xi }}^{bare}=(\mathrm{}^2^2)/(2m)+V_{trap}(𝐫)`$ contains both kinetic energy and trapping potential and $`V_{rsmn}`$ represents the symmetrized 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\}`$, where $`|i=\psi _i(𝐫)`$ denotes a single-particle eigenstate of the trap. 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)`$. We assume the system to be in a symmetry-broken phase and hence express the single-particle operators $`\widehat{a}_i(t)`$ as

$$\widehat{a}_i(t)=\widehat{a}_i+\left(\widehat{a}_i\widehat{a}_i\right)=z_i(t)+\widehat{c}_i(t)$$
(2)
This allows, in general, for a coherent mean-field amplitude $`z_i`$ to form in a number of low-lying trap levels. We can now formulate a non-equilibrium theory for the coupled evolution of condensate mean field amplitudes and fluctuations about these values, based on a suitably truncated hierarchy of coupled equations of motion and appropriate decoupling approximations.

Possibly the most direct approach for studying the dynamics in a closed system is based on solving an appropriate set of such equations self-consistently, in terms of exact interatomic potentials; such an approach has been discussed, for example, in . In a realistic system, the number of trap eigenstates will be very large, making such a procedure computationally very demanding. Since the most interesting dynamics take place in the low-lying levels, in this paper we have chosen to restrict our analysis to such levels (although such treatment will also implicitly yield the behavior of high-lying thermal levels). Accepting a distinction between low- and high-lying levels allows us to adiabatically eliminate all high-lying levels appearing as intermediate states in the equations for the evolution of averages of low-lying states; this procedure is known to lead to the renormalization of the exact (single-vertex) interatomic potential to an effective two-body one (over high-lying states), as discused in . We thus arrive at the situation where the hamiltonian of the system still has the general form of Eq. (1), but with the single-vertex interatomic potential $`V`$ replaced by an effective two-body T-matrix, $`T`$, over high-lying states, with the simultaneous restriction of all bare trap eigenstates being summed over low-lying levels. For sufficiently dilute systems at low temperature, this restricted effective two-body interaction is approximately equal to the full two-body T-matrix, $`T^{2B}`$, giving the scattering of two particles in vacuum. We thus approximate $`T`$ in terms of $`T^{2B}`$, by ensuring that purely two-body effects due to collisions occuring in vacuum are not double-counted (for a detailed discussion of the relation of these effective interactions and the renormalization required to avoid double-counting, the reader is referred to ). We note that it is precisely the quantity $`T^{2B}`$ which corresponds, in three dimensions, to the usual binary s-wave scattering length pseudopotential.

We are interested in working out in a self-consistent manner the evolution of condensate population and incoherent fluctuations about this value, and for consistency we work throughout this paper with averages of two single-particle operators. The Heisenberg equation of motion for such a general product of two operators is given by (setting $`\mathrm{}=1`$)

$`i{\displaystyle \frac{d}{dt}}\widehat{a}_i^{}\widehat{a}_j`$ $`={\displaystyle \underset{r}{\overset{^{}}{}}}\left\{\mathrm{\Xi }_{jr}^{Bare}\widehat{a}_i^{}\widehat{a}_r\mathrm{\Xi }_{ri}^{Bare}\widehat{a}_r^{}\widehat{a}_j\right\}`$ (4)
$`+{\displaystyle \underset{rms}{\overset{^{}}{}}}\left\{T_{jsmr}\widehat{a}_i^{}\widehat{a}_s^{}\widehat{a}_m\widehat{a}_rT_{ismr}^{}\widehat{a}_m^{}\widehat{a}_r^{}\widehat{a}_s\widehat{a}_j\right\}`$
where the primes indicate summation over low-lying levels. Consistent application of second order perturbation theory in the weakly-interacting limit (i.e. when the system can be well described in terms of single-particle wavefunctions) should yield correct expressions for energy level shifts and population damping. To proceed with our treatment, we must determine whether there are any quantities (fluctuations) which evolve faster than others (mean fields), so that the former can be adiabatically eliminated from the first order expressions. For example, in the usual rate equation treatments, one eliminates all off-diagonal normal and anomalous averages in favour of (diagonal) populations; one often speaks of coherences damping out faster than populations due to the coupling of the system to its environment. This gives rise to a set of equations to second order in the effective potentials, coupling populations to populations, in what is often termed the secular approximation .

One could, however, argue that the choice of which low-lying averages can be adiabatically eliminated depends on the basis employed for the description of the system, i.e. essentially on whether the mean field energy shifts are correctly taken into account or not. Starting from bare trap eigenenergies, one would not expect any normal or anomalous averages to be slowly-evolving; hence, the correct second order expression of Eq. (3) in a bare particle basis can be obtained by adiabatic elimination of the entire quantity $`\widehat{a}^{}\widehat{a}^{}\widehat{a}\widehat{a}`$ by means of its respective equation of motion. This is the procedure adopted in the main part of this work, which gives rise to complex equations appearing to contain additional terms when compared directly to similar treatments. However, a subsequent transformation to a dressed single-particle basis (i.e. dressed eigenenergies) shows clearly that such terms drop out from the respective equations of motion for populations in dressed eigenstates, as anticipated. Nonetheless, an approach in terms of bare single-particle eigenenergies allows one to draw important analogies between multi-mode condensation and multi-mode laser theory. Proceeding thus with the treatment in a bare basis, and taking the operator $`\widehat{\mathrm{\Xi }}^{bare}`$ to be diagonal, we obtain

$$i\frac{d}{dt}\widehat{a}_i^{}\widehat{a}_s^{}\widehat{a}_m\widehat{a}_r=(\omega _m+\omega _r\omega _i\omega _s)\widehat{a}_i^{}\widehat{a}_s^{}\widehat{a}_m\widehat{a}_r+F_i(T,\widehat{a},\widehat{a}^{};t)$$
(5)
where $`\omega _i`$ correspond to bare trap energies and $`F_i(T,\widehat{a},\widehat{a}^{};t)`$ defines the collisional evolution of $`\widehat{a}_i^{}\widehat{a}_s^{}\widehat{a}_m\widehat{a}_r`$ in such a basis by

$`F_i(T,\widehat{a},\widehat{a}^{})`$ $`=`$ $`{\displaystyle \underset{lt}{\overset{^{}}{}}}T_{mrlt}\widehat{a}_i^{}\widehat{a}_s^{}\widehat{a}_l\widehat{a}_t{\displaystyle \underset{pq}{\overset{^{}}{}}}T_{pqis}\widehat{a}_p^{}\widehat{a}_q^{}\widehat{a}_m\widehat{a}_r`$ (6)
$`+`$ $`{\displaystyle \underset{plt}{\overset{^{}}{}}}T_{prlt}\widehat{a}_p^{}\widehat{a}_i^{}\widehat{a}_s^{}\widehat{a}_m\widehat{a}_l\widehat{a}_t+{\displaystyle \underset{plt}{\overset{^{}}{}}}T_{pmlt}\widehat{a}_p^{}\widehat{a}_i^{}\widehat{a}_s^{}\widehat{a}_r\widehat{a}_l\widehat{a}_t`$ (7)
$``$ $`{\displaystyle \underset{pql}{\overset{^{}}{}}}T_{pqls}\widehat{a}_p^{}\widehat{a}_q^{}\widehat{a}_i^{}\widehat{a}_l\widehat{a}_m\widehat{a}_r{\displaystyle \underset{pql}{\overset{^{}}{}}}T_{pqli}\widehat{a}_p^{}\widehat{a}_q^{}\widehat{a}_s^{}\widehat{a}_l\widehat{a}_m\widehat{a}_r`$ (8)
Assuming real eigenvalues, we obtain the following exact integral relation

$`{\displaystyle \frac{d}{dt}}\widehat{a}_i^{}\widehat{a}_j=`$ $`i(\omega _j\omega _i)\widehat{a}_i^{}\widehat{a}_ji[{\displaystyle \underset{rms}{\overset{^{}}{}}}T_{jsmr}\widehat{a}_i^{}\widehat{a}_s^{}\widehat{a}_m\widehat{a}_re.c.]`$ (10)
$`\{{\displaystyle \underset{rms}{\overset{^{}}{}}}T_{jsmr}{\displaystyle _{t_0}^t}dt^{^{}}e^{i(\omega _m+\omega _r\omega _s\omega _i)(tt^{^{}})}F_i(T,\widehat{a},\widehat{a}^{};t^{^{}})+e.c.\}`$
where $`e.c.`$ stands for the exchange conjugate (i.e. conjugate expression with labels i and j interchanged). By using the definition $`\widehat{a}_i=z_i+\widehat{c}_i`$, we obtain all second-order collisional terms of our approach. However, for these to be useful, we must express them in terms of a closed system of equations, by imposing suitable approximations. Firstly we decouple averages containing more than two single-particle operators via

$$\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n\widehat{c}_r^{}\widehat{c}_m\widehat{c}_s^{}\widehat{c}_n+\widehat{c}_r^{}\widehat{c}_n\widehat{c}_s^{}\widehat{c}_m+\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_m\widehat{c}_n$$
(11)
and

$`\widehat{c}_p^{}\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_q\widehat{c}_l\widehat{c}_t`$ $`\widehat{c}_p^{}\widehat{c}_q\left(\widehat{c}_r^{}\widehat{c}_l\widehat{c}_s^{}\widehat{c}_t+\widehat{c}_r^{}\widehat{c}_t\widehat{c}_s^{}\widehat{c}_l\right)`$ (12)
$`+`$ $`\widehat{c}_p^{}\widehat{c}_l\left(\widehat{c}_r^{}\widehat{c}_q\widehat{c}_s^{}\widehat{c}_t+\widehat{c}_r^{}\widehat{c}_t\widehat{c}_s^{}\widehat{c}_q\right)+\widehat{c}_p^{}\widehat{c}_t\left(\widehat{c}_r^{}\widehat{c}_q\widehat{c}_s^{}\widehat{c}_l+\widehat{c}_r^{}\widehat{c}_l\widehat{c}_s^{}\widehat{c}_q\right)`$ (13)
$`+`$ $`\widehat{c}_q\widehat{c}_l\left(\widehat{c}_p^{}\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_t+\widehat{c}_p^{}\widehat{c}_s^{}\widehat{c}_r^{}\widehat{c}_t+\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_p^{}\widehat{c}_t\right)`$ (14)
$`+`$ $`\widehat{c}_q\widehat{c}_t\left(\widehat{c}_p^{}\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_l+\widehat{c}_p^{}\widehat{c}_s^{}\widehat{c}_r^{}\widehat{c}_l+\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_p^{}\widehat{c}_l\right)`$ (15)
$`+`$ $`\widehat{c}_l\widehat{c}_t\left(\widehat{c}_p^{}\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_q+\widehat{c}_p^{}\widehat{c}_s^{}\widehat{c}_r^{}\widehat{c}_q+\widehat{c}_r^{}\widehat{c}_s^{}\widehat{c}_p^{}\widehat{c}_q\right)`$ (16)
.

Based on our formulation in terms of effective interactions, we further impose the Markov approximation which takes the quantity $`F_i(T,\widehat{a},\widehat{a}^{})`$ out of the integrand and defines the intermediate propagators. This implies the assumption that the operators $`\widehat{a}_i`$ evolve freely between collisions via $`\widehat{a}_i(t^{^{}})=e^{+i\omega _i(tt^{^{}})}\widehat{a}_i(t)`$. Thus, the second order contributions of Eq. (10) acquire the general form

$$\left[\frac{d}{dt}\widehat{a}_i^{}\widehat{a}_i\right]_{T^2}=\underset{rms}{\overset{^{}}{}}\{\left[T_{ismr}\underset{\mathrm{}}{}\left(_{t_0}^td\tau e^{i(\mathrm{\Delta }\omega _{\mathrm{}})(tt^{^{}})}\right)T_{\mathrm{}}\stackrel{~}{F_i\mathrm{}}(t)\right]+e.c.\}$$
(17)
where the above quantity $`T_{\mathrm{}}\stackrel{~}{F_i\mathrm{}}(t)`$ corresponds to $`F_i`$ defined by Eq. (5), and this notation has been used to indicate that the dotted indices of the second $`T`$ are the same indices as the ones appearing in the exponential of the integrand, as a result of the Markov approximation. Following the notation of Gardiner et al. , we now re-write the above integral as

$$\underset{\eta 0^+}{lim}_0^{\mathrm{}}𝑑\tau e^{\pm i(\mathrm{\Delta }\omega \pm i\eta )\tau }=\delta ^{(p)}(\mathrm{\Delta }\omega )=\pi \delta (\mathrm{\Delta }\omega )\pm i𝒫\left(\frac{1}{\mathrm{\Delta }\omega }\right)$$
(18)
where $`𝒫\left(\frac{1}{\mathrm{\Delta }\omega }\right)`$ corresponds to the principal value of the integral and the upper limit of integration has been approximated by $`(tt_0)\mathrm{}`$, based on the usual assumption that successive collisional events are well separated in time. Having discussed our approximations, we can now write down the coupled rate equations for coherent, incoherent and total populations.

## III Rate Equations for Multi-level Condensation

In this section, we focus on the application of the above methodology to a simple three-level system which we discuss within the Popov approximation , in which anomalous averages of the non-condensate are ignored in the final expressions. We are well aware that this system is rather idealized and do not claim that it will accurately reproduce the entire dynamics of condensed and thermal atoms. In fact, corrections beyond Popov may be significant, as shown in this context by Walser et al.  (see also ), and we defer their explicit discussion to a subsequent paper (but see also Sec. IV. B). In the current paper we are mainly concerned with addressing the simplest departure from single-level condensation, and how this might affect the dynamics of low-lying levels. To this aim, such a system is ideal for a simple comparison with other kinetic theories . Perhaps more significantly, such a small system will enable us to discuss in an explicit manner the desired analogy of the multi-mode nature of a single condensate with multi-mode semiclassical laser theory . Our notation is as follows: The total population $`N_i`$ in level $`i`$ is written as a sum of coherent $`|z_i|^2`$ and incoherent $`n_i=\rho _{ii}`$ populations, via

$$N_i=\widehat{a}_i^{}\widehat{a}_i=\widehat{c}_i^{}\widehat{c}_i+|z_i|^2=n_i+|z_i|^2$$
(19)
with off-diagonal coherent and incoherent matrix elements respectively defined by $`\zeta _{ij}=\left(z_j^{}z_i\right)`$ and $`\rho _{ij}=\widehat{c}_j^{}\widehat{c}_i`$.

### A Evolution of Total Populations

We start with expressions for the coupled evolution of total populations in the three lowest-lying bare levels under the assumption that all of them may exhibit partial condensation. Since we are dealing with a closed system, the evolution of populations in these levels trivially satisfy

$$\frac{dN_0}{dt}=\frac{dN_2}{dt}=\frac{1}{2}\left(\frac{dN_1}{dt}\right)$$
(20)
where

$$\frac{dN_0}{dt}=2i[T_{0211}(\rho _{10}+\zeta _{10})\rho _{12}c.c.]+\{(R_{T_{|0211|^2}}+\stackrel{~}{Q}_{ijji}^{0211})+c.c.\}$$
(21)
Here we have defined the collisional redistribution rate

$`R_{T_{|0211|^2}}=`$ $`2\mathrm{\Gamma }_{0211}\left[(N_0+1)(N_2+1)N_1^2N_0N_2(N_1+1)^2\right]`$ (23)
$`\mathrm{\Gamma }_{0211}(1+N_0+N_2)|z_1|^4+\stackrel{~}{R}_{T_{|0211|^2}}`$
where

$$\mathrm{\Gamma }_{0211}=|T_{0211}|^2\underset{\eta 0^+}{lim}_0^{\mathrm{}}𝑑\tau e^{\pm i(\omega _0+\omega _22\omega _1\pm i\eta )\tau }=|T_{0211}|^2\delta ^{(p)}(\omega _0+\omega _22\omega _1)$$
(24)
and the ‘tilde’ in $`\stackrel{~}{R}_{T_{|0211|^2}}`$ and $`\stackrel{~}{Q}_{ijji}^{0211}`$ denotes contributions which are off-diagonal in $`\rho `$ or $`\zeta `$, whose explicit expressions can be found in Appendix A. In the context of the well-known Boltzmann scattering factors, the appearance of the additional term $`\mathrm{\Gamma }_{0211}|z_1|^4`$ may appear somewhat perplexing. We stress that such contributions only arise because of our choice to formulate our description in terms of a bare basis, and we explicitly show in Sec. IV B that such terms ‘disappear’ in the usual formulation in terms of a self-consistently dressed basis (where they become incorporated in the dressing of the bare trap eigenstates).

Unlike the evolution of total populations, the dynamics of condensed and thermal atoms will also depend on processes of the general form $`T_{ijji}T_{kllk}`$ which will be intrinsically dispersive. For simplicity, we shall henceforth limit our discussion to the case of only two partially condensed levels, namely $`z_0,z_10`$, with level $`2`$ being the lowest purely incoherent level treated here.

### B Evolution of Condensate Populations

The method used to obtain the evolution of condensate populations, is essentially the same as the one used for total populations (Eq. (17)).For ease of our subsequent comparison with the semi-classical laser theory, in this section we directly obtain equations for the modulus $`|z_i|^2`$ of the condensate mean field amplitude $`z_i`$ and not for $`z_i`$ itself. These are given by

$$\frac{d|z_i|^2}{dt}=\underset{rms}{\overset{^{}}{}}\{\left[iT_{ismr}z_i^{}\widehat{a}_s^{}\widehat{a}_m\widehat{a}_r\right]+c.c.\}$$
(25)
where the quantity analogous to $`F_i`$ of Eq. (5) is obtained from the equation of motion of the entire quantity $`z_i^{}\widehat{a}_s^{}\widehat{a}_m\widehat{a}_r`$ and not just $`\widehat{a}_s^{}\widehat{a}_m\widehat{a}_r`$ (see Appendix B). Thus, to second order in the effective interaction, we obtain

$$\frac{d|z_i|^2}{dt}=i\underset{rms}{\overset{^{}}{}}[T_{ismr}z_i^{}\widehat{a}_s^{}\widehat{a}_m\widehat{a}_rc.c.]+[X_{ii}(T^2)+c.c.]$$
(26)
with all $`T^2`$ collisional contributions summarized in $`X_{ii}(T^2)`$, which can be written as

$$X_{ii}(T^2)=R_i^{con}+S_i+E_i+\stackrel{~}{X}_{ii}(T^2)$$
(27)
$`R_i^{con}`$ corresponds to redistribution collisions affecting the level $`i`$ condensate, $`S_i`$ corresponds to self-interactions within level $`i`$, and $`E_i`$ gives the exchange collisions of level $`i`$ with its nearby levels, whereas the term $`\stackrel{~}{X}_{ii}(T^2)`$ corresponds to all off-diagonal contributions. For two condensed levels, in our simplified three-level system, the respective expressions are

$`R_0^{con}=\mathrm{\Gamma }_{0211}\left\{\begin{array}{c}2|z_0|^2\left[(n_2+1)n_1^2n_2(n_1+1)^2\right]\hfill \\ +4|z_0|^2|z_1|^2\left[(n_2+1)n_1n_2(n_1+1)\right]\hfill \\ +|z_0|^2|z_1|^4\left[(n_2+1)n_2\right]\hfill \end{array}\right\}`$ (31)
$`S_0=\lambda _{0000}\left[\begin{array}{c}2|z_0|^2\left[(n_0+1)^2n_0n_0^2(n_0+1)\right]\hfill \\ +|z_0|^4\left[(n_0+1)^2n_0^2\right]\hfill \end{array}\right]`$ (34)
$`E_0`$ $`=`$ $`4\lambda _{0110}\left[\begin{array}{c}|z_0|^2\left[(n_0+1)(n_1+1)n_1n_0n_1(n_1+1)\right]\hfill \\ +|z_0|^2|z_1|^2\left\{\begin{array}{c}\left[(n_0+1)(n_1+1)n_0n_1\right]\hfill \\ +\left[(n_0+1)n_1n_0(n_1+1)\right]\hfill \end{array}\right\}\hfill \end{array}\right]`$ (40)
$`4\lambda _{0220}|z_0|^2\left[(n_0+1)(n_2+1)n_2n_0n_2(n_2+1)\right]`$
and

$`R_1^{con}=2\mathrm{\Gamma }_{0211}\left\{\begin{array}{c}2|z_1|^2\left[(n_1+1)n_0n_2n_1(n_0+1)(n_2+1)\right]\hfill \\ +2|z_0|^2|z_1|^2\left[(n_1+1)n_2n_1(n_2+1)\right]\hfill \\ |z_1|^4\left[(n_0+1)(n_2+1)n_0n_2\right]\hfill \\ |z_0|^2|z_1|^4\left[(n_2+1)n_2\right]\hfill \end{array}\right\}`$ (45)
$`S_1=\lambda _{1111}\left[\begin{array}{c}2|z_1|^2\left[(n_1+1)^2n_1n_1^2(n_1+1)\right]\hfill \\ +|z_1|^4\left[(n_1+1)^2n_1^2\right]\hfill \end{array}\right]`$ (48)
$`E_1`$ $`=`$ $`4\lambda _{0110}\left[\begin{array}{c}|z_1|^2\left[(n_0+1)(n_1+1)n_0n_0n_1(n_0+1)\right]\hfill \\ +|z_0|^2|z_1|^2\left\{\begin{array}{c}\left[(n_0+1)(n_1+1)n_0n_1\right]\hfill \\ +\left[(n_1+1)n_0n_1(n_0+1)\right]\hfill \end{array}\right\}\hfill \end{array}\right]`$ (54)
$`4\lambda _{1221}|z_1|^2\left[(n_1+1)(n_2+1)n_2n_1n_2(n_2+1)\right]`$
where

$$\lambda _{ijji}=|T_{ijji}|^2\underset{\eta 0^+}{lim}_0^{\mathrm{}}𝑑\tau e^{\eta \tau }$$
(55)
while the expressions for the off-diagonal contributions $`\stackrel{~}{X}_{ii}(T^2)`$ can be found in Appendix A.

### C Evolution of Non-condensate Populations

The evolution of the normal averages of levels 0 and 1 is governed by the equation

$$\frac{d}{dt}\widehat{c}_i^{}\widehat{c}_i=\underset{rms}{\overset{^{}}{}}\{iT_{ismr}[\widehat{a}_i^{}\widehat{a}_s^{}\widehat{a}_m\widehat{a}_rz_i^{}\widehat{a}_s^{}\widehat{a}_m\widehat{a}_r]c.c.\}+[Y_{ii}(T^2)+c.c.]$$
(56)
where $`Y_{ii}(T^2)=R_i^{th}S_iE_i+\stackrel{~}{Y}_{ii}(T^2)`$ and $`R_i^{th}`$ are defined by

$`R_0^{th}=\mathrm{\Gamma }_{0211}\left\{\begin{array}{c}2\left[(n_0+1)(n_2+1)n_1^2n_0n_2(n_1+1)^2\right]\hfill \\ +4|z_1|^2\left[(n_0+1)(n_2+1)n_1n_0n_2(n_1+1)\right]\hfill \\ +|z_1|^4\left[(n_0+1)(n_2+1)n_0n_2\right]\hfill \end{array}\right\}`$ (60)
$`R_1^{th}=4\mathrm{\Gamma }_{0211}\left\{\begin{array}{c}\left[(n_1+1)^2n_0n_2n_1^2(n_0+1)(n_2+1)\right]\hfill \\ +|z_0|^2\left[(n_1+1)^2n_2n_1^2(n_2+1)\right]\hfill \\ +|z_1|^2\left[(n_1+1)n_0n_2n_1(n_0+1)(n_2+1)\right]\hfill \\ +|z_0|^2|z_1|^2\left[(n_1+1)n_2n_1(n_2+1)\right]\hfill \end{array}\right\}`$ (65)
The evolution of the non-condensate component of level 2 is given by Eqs. (12)-(14), in the limit $`z_2=0`$. The respective expressions for the off-diagonal contributions $`\stackrel{~}{Y}_{ii}(T^2)`$ can be found in Appendix A.

## IV Interpretation of Rate Equations

We now turn our attention to the interpretation of the above equations, focusing in Sec. A only on the condensate evolution. Ignoring, at first, re-distributional dynamics and off-diagonal contributions, we show how they reduce to the Hartree-Fock coupled nonlinear Schrödinger equations conventionally used for describing two-component condensates (Sec. IV A1). By additionally considering the redistribution dynamics, we further discuss the analogy of a single inherently multi-mode condensate to the semi-classical photon laser theory , with our discussion being limited to only two modes for clarity (Section IV A2). By further considering the evolution of uncondensed, and thus also of total populations of each level, we compare and contrast our approach to the kinetic theory of Walser et al. , which is the formal multi-mode kinetic theory closest to the formalism of this paper (Sec. IV B).

### A Coupled Condensate Evolution

#### 1 Hartree-Fock Theory for Binary Condensates

Let us initially forget about the presence of level 2 in our system and focus on the interactions between the two lowest levels in the trap (i.e. set $`n_2=z_2=0`$ and $`V_{0211},\lambda _{0220},\lambda _{1221}0`$), which amounts to discussing two coupled partially-condensed systems.

Using Eqs. (17)-(18) and (20)-(21) we can deduce the evolution of the ground state condensate mean field $`z_0`$ (as opposed to $`|z_0|^2`$), which to second order reads

$`{\displaystyle \frac{dz_0}{dt}}=i\omega _0z_0`$ $`i\left[T_{0000}i\lambda _{0000}(1+2n_0)\right]\left(2n_0+|z_0|^2\right)z_0+\lambda _{0000}\left(2n_0^2\right)z_0`$ (68)
$`2i\left[T_{0110}2i\lambda _{0110}(1+n_0+n_1)\right]\left(n_1+|z_1|^2\right)z_0+2\lambda _{0110}\left(2n_0n_1\right)z_0`$
$`4\lambda _{0110}\left(n_1n_0\right)|z_1|^2z_0`$
The terms in square brackets can be identified as the corresponding matrix elements of an effective many-body interaction introduced, strictly speaking, only over high-lying levels; this generalized effective interaction can be straightforwardly replaced by the usual many-body T-matrix $`T_{ijji}^{MB}`$ (over all levels), discussed in , by subtracting from the second order expression of Eq. (29) (i.e. terms $`\lambda _{ijji}`$), a term corresponding to the scattering of particles in vacuum, so that one ends up correctly calculating the change in the effective interaction when the pair of atoms collides in a condensed gas (as opposed to the vacuum) . Since our analysis has been carried out in terms of bare particle energies $`\omega _i`$, such a renormalization of the terms appearing in Eq. (29) merely amounts to neglecting the factor of 1 in $`(1+2n_0)`$ and $`(1+n_0+n_1)`$ . We note that the same procedure should be carried out in all second order terms whose scattering amplitude depends on factors of the form $`|T_{ijkl}|^2|z_i|^2|z_j|^2\left[(n_k+1)(n_l+1)n_kn_l\right]`$ (if these are to be written in terms of $`T^{MB}`$), and such terms appear in Eqs. (20)-(24) and (27). After ‘renormalization’, the square bracket in the first line of Eq. (29) takes the form

$$T_{0000}i\lambda _{0000}(1+2n_0)=T_{0000}^{2B}+T_{0000}^{2B}\underset{\eta 0^+}{lim}_0^{\mathrm{}}\frac{d\tau }{i}e^{\eta \tau }\left(2n_0\right)T_{0000}^{2B}=\left[T_{0000}^{MB}\right]_{T^2}$$
(69)
where $`\left[T_{0000}^{MB}\right]_{T^2}`$ corresponds to the second order expression for the many-body T-matrix; a similar identification can be made for the square bracket of the second line of Eq. (29) (where a factor of two arises from the symmetric interchange of atoms $`0`$ and $`1`$ corresponding to direct and exchange Hartree-Fock terms).

If one further ignores the last term of Eq. (29) (which corresponds diagrammatically to ignoring bubble diagrams with respect to the bare many-body ladder diagrams), one obtains the lowest order expression of the general equation

$`i{\displaystyle \frac{dz_0}{dt}}=\omega _0z_0`$ $`+`$ $`T_{0000}^{MB}\left(2n_0+|z_0|^2\right)z_02T_{0000}^{MB}\left[\underset{\eta 0^+}{lim}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{d\tau }{i}}e^{\eta \tau }\right]\left(n_0^2\right)T_{0000}^{MB}z_0`$ (70)
$`+`$ $`2T_{0110}^{MB}\left(n_1+|z_1|^2\right)z_02T_{0110}^{MB}\left[\underset{\eta 0^+}{lim}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{d\tau }{i}}e^{\eta \tau }\right]\left(2n_0n_1\right)T_{0110}^{MB}z_0`$ (71)
with the last term in each line ensuring correct scattering factors for condensate feed collisions from thermal atoms. This expression agrees with the equation for condensate evolution derived by one of us elsewhere . Although this equation is valid when dealing with a single condensed trap level, the somewhat heuristic neglect of $`4\lambda _{0110}\left(n_1n_0\right)|z_1|^2z_0`$ mentioned above suggests that it must be interpreted with some caution when the condensate spans more than one bare trap eigenstates .

We now show explicitly how the above equation reduces to the coupled finite temperature nonlinear Schrödinger equations used for two-component condensation . For this we must first assume that the gas is sufficiently dilute, so that we can further ignore in Eq. (31) the effect of the surrounding medium on binary collisions; this amounts to replacing $`T_{0ii0}^{MB}`$ ($`i=0`$, $`1`$) by the full two-body T-matrix $`T_{0ii0}^{2B}`$, while simultaneously ignoring the ‘kinetic’ (many-body) contributions corresponding to the last term of each line in Eq. (32). By further approximating $`T^{2B}`$ by the usual pseudopotential $`U_{ij}\delta (𝐫𝐫^{^{}})`$ valid in 3D, where $`U_{ij}`$ is parametrized in terms of the scattering length for the collision of an atom in level $`i`$ with an atom in level $`j`$ ($`i`$,$`j=0,1`$) and transforming to coordinate space, we obtain

$$i\frac{\mathrm{\Phi }_0}{t}=H_0^{(0)}\mathrm{\Phi }_0+U_{00}\left(|\mathrm{\Phi }_0|^2+2\stackrel{~}{n}_0\right)\mathrm{\Phi }_0+2U_{01}\left(|\mathrm{\Phi }_1|^2+\stackrel{~}{n}_1\right)\mathrm{\Phi }_0$$
(72)
where $`\mathrm{\Phi }_i`$ and $`\stackrel{~}{n}_i`$ denote the condensate mean field and the non-condensate density of component $`i`$, with $`H_i^{(0)}`$ being the corresponding bare trap hamiltonian (kinetic energy plus trapping potential). By analogy,

$$i\frac{\mathrm{\Phi }_1}{t}=H_1^{(0)}\mathrm{\Phi }_1+U_{11}\left(|\mathrm{\Phi }_1|^2+2\stackrel{~}{n}_1\right)\mathrm{\Phi }_1+2U_{01}\left(|\mathrm{\Phi }_0|^2+\stackrel{~}{n}_0\right)\mathrm{\Phi }_1$$
(73)
Eqs. (32)-(33) correspond to the well-known finite temperature Hartree-Fock equations for two-component condensation .

#### 2 Two-mode Condensation vs. Semi-classical Laser Theory

In the language of quantum optics, a system in which the condensate mean field spans more than one single-particle eigenstate should be analogous to a photon laser in multi-mode operation. Exploring the formal connection between the two, with their similarities and differences, was in fact part of the motivation for this work and the main reason for formulating our approach in terms of bare particle eigenenergies. If such a connection were to be taken literally, one would expect the coupled condensate mean field equations to resemble those of the two-mode laser. Even though we have assumed that level 2 is not itself condensed (and may even be initially unoccupied even by thermal atoms), we should still also consider its presence here; this is because the existence of level 2 is an inherent property of the system that cannot be ignored, since it will affect the evolution of the two lowest condensed modes $`z_0`$, $`z_1`$ via collisional redistribution processes. Combining Eqs. (19)-(24) for the second order contributions to the evolution of condensate amplitudes in a bare single-particle basis, they are found to exhibit the general structure

$$\frac{d|z_0|^2}{dt}=\alpha _0^{RSEH}|z_0|^2\beta _0^S|z_0|^4\theta _{01}^{RE}|z_0|^2|z_1|^2+\xi _{01}^R|z_0|^2|z_1|^4$$
(74)
$$\frac{d|z_1|^2}{dt}=\alpha _1^{RSEH}|z_1|^2\beta _1^{RS}|z_1|^4\theta _{10}^{RE}|z_0|^2|z_1|^2+\xi _{10}^R|z_0|^2|z_1|^4$$
(75)
Focusing initially only on the first three contributions of each of the above coupled equations, we note that they have the same form as those of the two-mode photon laser intensity equations , if we identify $`|z_i|^2`$ with the mode intensity. By analogy, we thus refer to the coefficients appearing in Eqs. (34)-(35) as: the ‘net gain’ coefficient $`\alpha _i`$ of each mode, the ‘self-saturation’ $`\beta _i`$ of each mode and the ‘cross-saturation’ coefficients $`\theta _{ij}`$. In contrast to the optical laser where such equations arise from the polarization of the active medium, the nonlinearity in the case of the atom laser is intrinsic, arising from atom-atom interactions; this is absent in optical lasers, since there are no photon-photon interactions affecting the coherent photon field. Hence, the above coefficients $`\alpha _i`$, $`\beta _i`$ and $`\theta _{ij}`$ depend on incoherent populations $`n_i`$ and collisional rates $`\lambda _{ijji}`$ and $`\mathrm{\Gamma }_{0211}`$. The superscripts $`R`$, $`S`$, $`E`$, $`H`$ used to define the above coefficients stand for Redistribution, Self, Exchange and Higher-level-exchange terms. The above coefficients are respectively defined by

$$\alpha _i^{RSEH}=\alpha _i^R\alpha _i^S\alpha _{ij}^E\alpha _{ik}^H$$
(76)
where

$`\left\{\begin{array}{c}\alpha _0^R=2\mathrm{\Gamma }_{0211}\left[(n_2+1)n_1^2n_2(n_1+1)^2\right]\hfill \\ \alpha _1^R=2\mathrm{\Gamma }_{0211}\left[n_0n_2(n_1+1)(n_0+1)(n_2+1)n_1\right]\hfill \\ \alpha _i^S=2\lambda _{iiii}\left[(n_i+1)^2n_i(n_i+1)n_i^2\right]\hfill \\ \alpha _{ij}^E=4\lambda _{ijji}\left[(n_i+1)(n_j+1)n_j(n_j+1)n_in_j\right]\hfill \\ \alpha _{ij}^H=4\lambda _{ikki}\left[(n_i+1)(n_k+1)n_k(n_k+1)n_in_k\right]\hfill \end{array}\right\}`$ (82)
$$\beta _i^{(R)S}=\left(\beta _i^R\right)+\beta _i^S$$
(83)
where $`\beta _0^R=0`$ and

$`\left\{\begin{array}{c}\beta _1^R=2\mathrm{\Gamma }_{0211}\left[(n_0+1)(n_2+1)n_0n_2\right]\hfill \\ \beta _i^S=\lambda _{iiii}\left[(n_i+1)^2n_i^2\right]\hfill \end{array}\right\}`$ (86)
$$\theta _{ij}^{RE}=\theta _{ij}^R+\theta _{ij}^E$$
(87)
$`\left\{\begin{array}{c}\theta _{01}^R=4\mathrm{\Gamma }_{0211}\left[n_1(n_2+1)(n_1+1)n_2\right]=\theta _{10}^R\hfill \\ \theta _{ij}^E=4\lambda _{ijji}\left\{\left[(n_i+1)(n_j+1)n_1n_j\right]+\left[(n_i+1)n_j(n_j+1)n_i\right]\right\}\hfill \end{array}\right\}`$ (90)
$$\xi _{01}^R=\mathrm{\Gamma }_{0211}\left[(n_2+1)n_2\right]=\frac{1}{2}\xi _{10}^R$$
(91)
Knowledge of their detailed form, enables us to draw important conclusions regarding the signs of the total coefficients $`\alpha _i`$, $`\beta _i`$, $`\theta _{ij}`$ and $`\xi _{ij}`$ of Eqs. (34)-(35). In particular we find, just as in the photon laser, the coefficients $`\beta _i>0`$ always, thus giving rise to self-saturation of the mode intensity $`|z_i|^2`$, whereas coefficients $`\alpha _i`$ and $`\theta _{ij}`$ can be positive or negative (depending on the values of $`n_i`$, $`\lambda _{ijji}`$ and $`\mathrm{\Gamma }_{0211}`$). The net gain coefficient of each mode is given by Eq. (36); here $`\alpha _i^R`$ gives the gain coefficient due to redistribution collisions $`|T_{0211}|^2`$ and can be positive or negative, depending on the relative values of $`n_0`$, $`n_1`$ and $`n_2`$. The remaining contributions to Eq. (36) arise from collisions of an atom in level $`i`$ with an atom in level $`j`$ or $`k`$ (where $`i=0,1`$ and $`j=0,1,2`$, while $`ki,j`$), and leads to saturation of intensity growth since $`\alpha _i^S,\alpha _i^E,\alpha _i^H>0`$ always. The cross-saturation rates $`\theta _{ij}`$ contain contributions from exchange and redistribution collisions and can be positive or negative.

An important difference between Eqs. (34)-(35) and the corresponding ones for two-mode photon lasers is the presence of the higher order cross-saturation terms $`\xi _{ij}|z_0|^2|z_1|^4`$, which arise solely as a result of collisional redistribution. In deriving the two-mode photon laser equations, one conventionally performs a perturbative expansion of the polarization of the active medium in terms of laser intensity. For near-threshold operation and weak laser intensities, the third order perturbative expansion of the medium polarization is usually sufficient, and such a truncation generates at most terms proportional to the square of the laser intensity. Terms of third order in the laser intensity (as well as higher order ones) would indeed arise in two-mode photon laser theory, if one extended the perturbative treatment of the polarization to fifth order (or beyond).

In the case of Bose-Einstein condensation there is no active medium to be polarized, and the nonlinearity of the system is due to intrinsic atom-atom interactions. In this case, the highest order $`|z_0|^{2m}|z_1|^{2n}`$ to be included in the two-mode equations for the ‘intensities’ $`|z_i|^2`$ is based on the number of single-particle operators appearing in the nonlinear interaction term in the hamiltonian of the system, i.e. whether one includes only two-body collisions (as usual, in the dilute limit $`na^31`$), or also three-body collisions (or higher). In the usual case of two-body collisions as defined by Eq. (1), the resulting equations for $`|z_i|^2`$ may contain terms up to order $`(z^{}z)^3`$; so, in general, it would not appear justified to ignore contributions of third order in the mode intensities. We note that the coefficients of the higher order cross-saturation terms appearing in Eqs. (34)-(35) have opposite signs and are given by $`\xi _{10}=2\xi _{01}`$, where $`\xi _{01}=\mathrm{\Gamma }_{0211}>0`$.

The above discussion has been given in terms of bare amplitudes $`|z_i|^2`$, and we should note that the form of these equations changes when the $`|z_i|^2`$ refer to amplitudes of dressed modes. In particular, when the $`|z_i|^2`$ refer to amplitudes in modes dressed only by the condensate mean field (assuming it is physically meaningful to speak of more than a single condensed mode in such a dressed basis), such higher order cross-coupling contributions would not arise, thus yielding a direct analogy with the intensity equations for two-mode photon lasers. The justification for this is given in the next section, where we discuss how the second order collisional integrals become modified, upon shifting our single-particle eigenstates to a basis dressed by mean fields. The ‘exact’ analogy between two-mode condensation and optical lasers arising in this case might suggest that even in a basis dressed by the condensate mean field, one should, in principle, deal with more than one condensed modes in a trapped assembly. This analogy of a single, inherently multi-mode, condensate with the usual semi-classical multi-mode photon laser theory, however, appears to break down when shifting to a basis including higher mean field effects (i.e. those due to uncondensed atoms, anomalous averages, etc.)

In this section, we have presented some similarities and differences between two-mode Bose-Einstein condensation and semi-classical two-mode photon laser theory. A more detailed investigation should discuss such equations in the presence of pumping, evaporative cooling and coherent outcoupling (in a manner analogous to single-mode atom laser models ). The first obvious modification that would occur in this case is that the ‘net gain’ coefficient(s) of the lasing mode(s) should become positive and large, since the contribution of an irreversible evaporative cooling mechanism combined with the redistribution collisions should result in a large flow of particles towards such mode(s). At the same time, of course, the $`\alpha _i`$ coefficients will acquire an additional negative contribution whose value will depend on the rate of outcoupling of atoms from the particular condensed level. This should enable the system to reach a steady state, with the coherent amplitude growth being stabilized both by the outcoupling mechanism, as well as by the collisionally-induced dephasing  due to self ($`\lambda _{iiii}`$) and exchange ($`\lambda _{ijji}`$) interactions. In our formalism, such dephasing can be viewed as destruction of the coherent mean field amplitude $`|z_i|^2`$ in favour of the ‘incoherent populations’ $`n_i`$ (i.e. transfer of population from ‘coherent’ to ‘incoherent’). Inherent two-body and three-body inelastic loss processes, which have recently been shown to be essential for reaching a steady state for an atom laser , will affect the coherent and incoherent populations of a particular mode in a different manner due to the nature of our mean field decorrelation, e.g. $`\widehat{a}_i^{}\widehat{a}_i^{}\widehat{a}_i^{}\widehat{a}_i\widehat{a}_i\widehat{a}_i\left(|z_i|^6+6n_i^3\right)`$, thus creating an essential irreversibility for the condensate mean field $`|z_i|^2`$ to dominate over its corresponding fluctuations $`n_i`$. By only outcoupling the lowest condensed mode of the system, one could, for example, achieve a kind of ‘population inversion’, in the sense that the largest condensate particle number accumulation might occur in some state other than the ground state of the bare trap. In closing this subsection, we note that the quantum nature of the multi-mode atom laser is implicit in the general Fokker-Planck treatment developed by Stoof .

### B Links to Other Kinetic Theories

The expressions given earlier for the evolution of condensed and uncondensed components contain terms in the condensate mean field beyond order $`(z^{}z)`$, and at first sight this appears to be in disagreement with existing kinetic theories . For example, comparing our final expressions for total populations to those of Walser et al.  (in the corresponding limit of diagonal populations and no anomalous averages), we appear to obtain an excess $`\mathrm{\Delta }N`$ for the populations of levels 0 and 2 (whereas level 1 is underestimated by twice that amount), with $`\mathrm{\Delta }N`$ given by

$$\mathrm{\Delta }N=2\mathrm{\Gamma }_{0211}|z_1|^2\left\{\begin{array}{c}4|z_0|^2\left[(n_2+1)n_1n_2(n_1+1)\right]\hfill \\ +|z_0|^2|z_1|^2\left[(n_2+1)n_2\right]\hfill \\ +|z_1|^4\left[(n_0+1)(n_2+1)n_0n_2\right]\hfill \end{array}\right\}$$
(92)
Since all terms of Eq. (43) are proportional to $`|z_1|^2`$, it is clear that our equations will assume the usual form (e.g. as discussed in Walser et al. ) in the limit of extremely weak condensation (i.e. when only the lowest trap eigenstate is occupied by the condensate). The natural question arising then is whether such equivalence remains beyond this simple limit.

The key to understanding this apparent inconsistency is to note that our treatment has so far been given in terms of single-particle eigenfunctions, whose energies are bare (unshifted) ones, as if the trap were completely void. However, since the trap contains a large atomic medium which is condensed, the collisions will actully occur in the presence of the condensate, which forces the eigenenergies to vary in time due to the existing mean field potentials. To account for this, one can thus convert from the simple picture in terms of unshifted energies employed above, to the conventional one in which the effects of the mean fields are included into suitably ‘renormalized’ basis eigenenergies . In the next section we show explicitly that shifting our basis in a suitable manner identically reproduces the kinetic theory of Walser et al. , thus providing an alternative derivation of the latter theory.

To understand the relation of our theory to existing kinetic treatments, we must, in first instance, explain the physical origin of all terms of Eq. (43). We recall from Eq. (3) that all collisional contributions affecting the total population of a particular level $`i`$ arise from terms of general structure $`[_{rms}T_{ismr}\widehat{a}_i^{}\widehat{a}_s^{}\widehat{a}_m\widehat{a}_r+h.c.]`$. Defining in the usual manner  the quantities $`\rho _{ji}=\widehat{c}_i^{}\widehat{c}_j`$ and $`\kappa _{jk}=\widehat{c}_k\widehat{c}_j`$ and focusing initially on the first line of Eq. (43), we note that terms $`|z_0|^2|z_1|^2`$ arise from the elimination of $`\rho `$ appearing in $`a^{}a^{}aa[4\rho (z^{}z)+c.c.]`$, which, among other contributions leads to

$$\frac{dN_i}{dt}=4\underset{rms}{}T_{ismr}z_i^{}\left\{\underset{klt}{}T_{mltk}^{(\delta )}\left[\rho _{ks}(z_l^{}z_t)\right]T_{klts}^{(\delta )}\left[\rho _{mk}(z_l^{}z_t)\right]\right\}z_r+c.c.$$
(93)
Here the notation $`T^{(\delta )}`$ stands for the corresponding two-body t-matrix multiplied by the approximately energy-conserving integral of Eq. (10), and the above term clearly contains a contribution $`\{|T_{0211}|^2z_0^{}\left[(n_1n_2)(z_1^{}z_0)\right]z_1+c.c.\}`$ These latter terms would obviously not arise if the quantity $`\rho `$ was considered essentially constant, as might be more appropriate, for example, for a rapidly-thermalizing ‘hydrodynamic’ system. Similarly, terms of order $`|z_i|^4`$ (i.e. second line of Eq. (43)) arise from the elimination of the pair anomalous average $`\kappa `$ in $`a^{}a^{}aa[\kappa (z^{}z^{})+c.c.]`$, yielding (among other terms)

$$\frac{dN_i}{dt}=\underset{rms}{}T_{ismr}\left(z_i^{}z_s^{}\right)\left\{\underset{pq}{}\left[T_{mrpq}^{(\delta )}+\underset{l}{}\left(T_{mlpq}^{(\delta )}\rho _{rl}+V_{lrpq}^{(\delta )}\rho _{ml}\right)\right]\left(z_pz_q\right)\right\}+c.c.$$
(94)
which leads to $`\{|T_{0211}|^2(z_1^{}z_1^{})(1+n_0+n_2)(z_1z_1)+c.c.\}`$ This proves that our expressions explicitly contain anomalous condensate terms $`(zz)`$ and $`(z^{}z^{})`$, in contrast to the non-condensate anomalous average $`\kappa `$ which has been neglected from all final expressions due to the application of the Popov approximation.

Finally, terms of order $`(z^{}z)^3`$ (last contribution of Eq. (43)) can only arise from elimination of $`z`$ (or $`z^{}`$) wherever it appears in the first order expressions, thus yielding

$`{\displaystyle \frac{dN_i}{dt}}`$ $`={\displaystyle \underset{rms}{}}T_{ismr}z_i^{}\left\{z_s^{}z_m\left[{\displaystyle \underset{pql}{}}T_{rlpq}^{(\delta )}(z_l^{}z_pz_q)\right]+z_s^{}\left[{\displaystyle \underset{pql}{}}T_{mlpq}^{(\delta )}(z_l^{}z_pz_q)\right]z_r\right\}`$ (96)
$`+{\displaystyle \underset{rms}{}}T_{ismr}\left\{z_i^{}\left[{\displaystyle \underset{pql}{}}T_{pqls}^{(\delta )}(z_p^{}z_q^{}z_l)\right]z_mz_r+\left[{\displaystyle \underset{pql}{}}T_{pqli}^{(\delta )}(z_p^{}z_q^{}z_l)\right](z_s^{}z_mz_r)\right\}+c.c.`$
where the products $`(z^{}z^{()}z)`$ appearing immediately after each $`T_{\mathrm{}}^{(\delta )}`$ are now the ones arising from the elimination of a $`z`$ or $`z^{}`$ from the first order expression in $`T(z^{}z^{}zz)`$. This allows us to interpret the term $`|T_{0211}|^2|z_0|^2|z_1|^4`$ of Eq. (43) as arising from the terms $`(z^{}zz)`$ generated in second order expressions, and these could be loosely interpreted as arising either (i) from condensate anomalous averages of the form $`(z_iz_j)`$, or (ii) from off-diagonal coherent terms $`(z_i^{}z_j)`$, with the two statements being equivalent and showing that it would not be justified to include off-diagonal $`(z_i^{}z_j)`$, while ignoring condensate anomalous averages $`(z_iz_j)`$.

In short, the above analysis shows clearly that all terms $`(z^{}z)`$, $`(zz)`$ and $`(z^{}z^{})`$ of second or higher order appearing in Eqs. (13), (17), (26) and hence (43) can only be generated by the adiabatic elimination of quanitites $`z`$, $`\rho `$ and $`\kappa `$ from the corresponding first order expressions. In particular, the terms of Eq. (43) respectively arise from the adiabatic elimination of the quantities $`\rho _{12}`$, $`z_2`$ and $`\kappa _{02}`$. Such elimination is fully justified in a bare basis, which assumes that all quantities evolve in an analogous (rapid) fashion and should thus be treated consistently. However, tranformation to a dressed basis automatically restricts certain quantities to be slowly-evolving (which parameters are slowly-evolving depends explicitly on the choice of unperturbed basis dressing the single-particle eigenstates), so that their elimination can no longer be justified. In the next section we thus extend the above discussion to show explicitly that our treatment yields precisely the second order collisional integrals of Walser et al. , upon renormalizing our single-particle eigenenergies to those dressed by HFB mean fields.

#### 1 Equivalence to Theory of Walser et al.

To show explicitly the link of our treatment to the kinetic theory of Walser et al. , we now focus on the general expression for the normal uncondensed component in an n-level system. To establish exact analogy, in this section we further incorporate the uncondensed anomalous average $`\kappa `$ into our treatment, since this is explicitly present in the expressions of Walser et al. .

After much algebraic manipulation of the second order contributions to Eq. (26), we hence obtain the following result for the second order collisional integrals arising within our bare basis

$`\left[{\displaystyle \frac{d}{dt}}(\rho _{ii})\right]_{T^2}^{Bare}=\left[{\displaystyle \frac{d}{dt}}(\rho _{ii})\right]_{T^2}^{Walser}+({\displaystyle \underset{rms}{}}T_{ismr}\left\{\begin{array}{c}\left(2\rho _{ms}\rho _{ri}+\kappa _{mr}\kappa _{is}^{}\right)\hfill \\ +\left(2\rho _{mi}z_s^{}z_r+z_mz_r\kappa _{is}^{}\right)\hfill \end{array}\right\}+h.c.)`$ (99)
The first term in the above expression has exactly the same form as the second order collisional integrals obtained by Walser et al.  (with the only difference being that in our approach the eigenenergies of Eq. (10) are now explicitly given by their bare values $`\omega _i`$), while the second term yields all the additional terms of our treatment arising when working within a bare-particle basis. The terms within the curly bracket appearing in the latter term denote the corresponding second order contributions arising from the adiabatic elimination of the quantities within the brackets, in the particular basis chosen. Their explicit form is obtained by the formal solution of their respective equations of motion given below, and such second order contributions therefore implicitly include the energy-conserving integral of Eq. (10). In the full hamiltonian of Eq. (1), the equations of motion for the above quantities are given by

$$i\mathrm{}\frac{dz_n}{dt}=\underset{k}{}\mathrm{\Xi }_{nk}^{Bare}z_k+\underset{ijk}{}T_{nijk}\left[z_i^{}z_jz_k+\kappa _{jk}z_i^{}+2\rho _{ji}z_k\right]$$
(100)
$`i\mathrm{}{\displaystyle \frac{d\rho _{ji}}{dt}}`$ $`=`$ $`{\displaystyle \underset{n}{}}\left(\mathrm{\Xi }_{jn}^{Bare}\rho _{ni}\mathrm{\Xi }_{ni}^{Bare}\rho _{jn}\right)+{\displaystyle \underset{r}{}}\left[\eta _{jr}\rho _{ri}\rho _{jr}\eta _{ri}\right]`$ (102)
$`{\displaystyle \underset{r}{}}\left[\kappa _{jr}\mathrm{\Delta }_{ri}^{}\mathrm{\Delta }_{jr}\kappa _{ri}^{}\right]`$
$`i\mathrm{}{\displaystyle \frac{d\kappa _{kj}}{dt}}`$ $`=`$ $`{\displaystyle \underset{n}{}}\left(\mathrm{\Xi }_{jn}^{Bare}\kappa _{nk}+\mathrm{\Xi }_{kn}\kappa _{jn}\right)+{\displaystyle \underset{s}{}}\left[\eta _{ks}\kappa _{sj}+\kappa _{ks}\eta _{sj}^{}\right]`$ (104)
$`+\mathrm{\Delta }_{kj}+{\displaystyle \underset{s}{}}\left[\rho _{ks}\mathrm{\Delta }_{sj}+\mathrm{\Delta }_{ks}\rho _{sj}^{}\right]`$
where we have ignored in the final expressions higher order correlations such as $`\widehat{c}^{}\widehat{c}\widehat{c}`$ and $`\widehat{c}\widehat{c}\widehat{c}`$ which arise in the above expressions (their effect will be discussed elsewhere). In the above expressions, we have further defined $`\eta _{pq}=2_{kl}T_{pklq}\left[z_k^{}z_l+\rho _{lk}\right]`$ and $`\mathrm{\Delta }_{pq}=_{kl}T_{pqkl}\left[z_kz_l+\kappa _{kl}\right]`$. Choosing to describe our system in terms of bare trap eigenenergies essentially amounts to diagonalizing the bare trap hamiltonian $`H_0=_{ij}\mathrm{\Xi }_{ij}^{Bare}\widehat{a}_i^{}\widehat{a}_j`$. In that case, the energies $`\omega _i`$ appearing in the energy-conserving condition (10) correspond to the eigenstates of the harmonic trap.

However, instead of working with a bare basis, one can choose a dressed ‘unperturbed’ basis $`(H_0+H_Q)`$ in which to define the renormalized single-particle eigenenergies, via

$$\widehat{H}=\left(H_0+H_Q\right)+\left(VH_Q\right)$$
(105)
where $`V`$ corresponds to the binary collision hamiltonian defined in Eq. (1) and $`H_Q`$ denotes the quasiparticle hamiltonian

$$H_Q=\frac{1}{2}\underset{pq}{}\left\{h_{pq}\left(\widehat{c}_p^{}\widehat{c}_q+\widehat{c}_q\widehat{c}_p^{}\right)+\left(\mathrm{\Delta }_{pq}\widehat{c}_p^{}\widehat{c}_q^{}+\mathrm{\Delta }_{pq}^{}\widehat{c}_q\widehat{c}_p\right)\right\}$$
(106)
with $`h_{qp}=_{pq}\left(\mathrm{\Xi }_{pq}^{Bare}+\eta _{pq}\right)`$. In this case, one is incorporating part of the collisional evolution of Eqs. (48)-(50) into their ‘free evolution frequencies’ $`\omega _i`$, which now become dressed to $`\stackrel{~}{\omega }_i`$. We thus find that when working in the usual quasiparticle basis, and upon ignoring higher order (triplet) averages, the above equations for $`z`$, $`\rho `$ and $`\kappa `$ simplify to $`(dz_n/dt)=i\stackrel{~}{\omega }_nz_n`$, $`(d\rho _{ji}/dt)=i\left(\stackrel{~}{\omega }_j\stackrel{~}{\omega }_i\right)\rho _{ji}`$ and $`(d\kappa _{jk}/dt)=i\left(\stackrel{~}{\omega }_j+\stackrel{~}{\omega }_k\right)\kappa _{jk}`$ where the $`\stackrel{~}{\omega }_i`$ now correspond to the renormalized eigenenergies in the above chosen basis. Hence (assuming as usual that initial interparticle correlations can be ignored), the quantities $`z`$, $`\rho `$ and $`\kappa `$ lead to no intermediate collisional evolution (to lowest order), which shows clearly that all terms of second or higher order in $`(z^{}z)`$ or $`(zz)`$ indeed vanish when dealing with single-particle eigenenergies dressed by the HFB mean fields; in the latter case our theory conincides with that of Walser et al (who actually formulated their theory in terms of unspecified number-conserving renormalized potentials as can be seen from Eqs. (64)-(67) of  ).

Furthermore, this treatment brings out the implicit basis dependence of the second order collisional integrals of Walser et al. In their treatment, they assume that the evolution of the system can be well parametrized by a restricted set of slowly-varying ‘master’ variables. As such, the possibility of evolution of such variables is absent from their theory, and hence they obtain a basis-independent formulation. Our treatment, however, indicates that their choice of a renormalized single-particle basis is implicit in their particular choice of the slowly-varying quantities. In a single-particle basis where the eigenenergies are dressed by HFB mean fields $`z`$, $`\rho `$ and $`\kappa `$, the collisional integrals cannot include contributions due to the variations of these HFB parameters which are assumed to be static and hence our treatment identically coincides with that of Walser et al., defining at the same time the eigenenergies $`\stackrel{~}{\omega }_i`$ as being dressed by HFB self-consistent potentials. This should be contrasted to the dressing caused simply by number-conserving HF mean fields which is implicit in Eq. (50) of their treatment. However, if one works within a bare basis, there is no a priori reason to assume that any quantity is slowly varying with respect to the others, and this leads to the additional second order contributions we have obtained in this paper (while simultaneously restricting the eigenstates of Eq. (10) to those of the harmonic trap). Hence, we believe our analysis further confirms the statement made in the last sentence of their paper, that their results should be valid when the mean-field induced energy shifts may be neglected during a strong collisional event, and we thus interpret the basis correction terms of our treatment as the evolution of these shifts in a bare basis. Our treatment can be readily generalized to all basis-dependent shifts, by respectively defining other dressed basis as follows: (i) the $`T=0`$ Gross-Pitaevskii basis in which only the condensate parameters $`z`$ evolve slowly and (ii) the Hartree-Fock basis in which both $`z`$ and $`\rho `$ are taken as slowly-varying, whereas the anomalous average $`\kappa `$ is adiabatically eliminated (even over low-lying<sup>*</sup><sup>*</sup>*We remind the reader that the corresponding values of all quantities over high-lying modes have already been implicitly adiabatically eliminated, in favour of the effective two-body interaction $`T`$ over low-lying modes used throughout this work. modes). Finally, in the case of an HFB formulation, such extra terms will only arise from triplet and higher-order averages ignored in Eqs. (48)-(50) and will clearly be of higher (than second) order in the interatomic potential, as also implied by Walser et al. .

## V Conclusions

In this paper we have discussed the dynamics of the lowest-lying single-particle bare trap levels in the weak coupling limit, under the assumption that more than one level exhibit condensation. By comparing our equations to those of the two-mode photon laser, we briefly commented on the similarities and differences of such equations. We note that the usual discussion of two-mode condensation focuses on condensates which are either in different spin states, or under physical separation. Describing each of the two separate condensates by a nonlinear Schrödinger equation (as discussed in Sec. IV), one indeed recovers a direct analogy to the semi-classical theory of two-mode photon lasers (by associating the fields $`|\mathrm{\Phi }_i(r,t)|^2`$ to the laser intensity). However, since the multi-mode nature of the photon laser refers to different modes of the electromagnetic field, we believe that a more direct analogy can be obtained at a more fundamental level, namely within a single trapped condensate (see also ). In particular, in the case of photon lasers, the number of modes and the precise equations governing them depend on the atomic species, the cavity, the intensity and the truncation of the polarization expansion imposed. In the end, the modes into which the field will oscillate result from the strong coupling of the active medium to the cavity, corresponding thus, in some sense, to a dressed basis. In the case of a trapped Bose-Einstein condensate, the number of modes spanned by the condensate is determined by its size, its diluteness (i.e. whether the discussion can be essentially limited to two-body interactions at the particular temperatures and densities), the strength of the interactions and the trap confinement power, whereas the form of the equations additionally depends on the basis in which we choose to describe the system. In this case, the maximum factor $`|z_i|^{2m}|z_j|^{2n}`$ which can be obtained in the equations is set by the complexity of the many-body interactions (i.e. two-body, three-body, etc.), whereas the factor we actually find in our final equations depends on the basis in which the equations are explicitly formulated. Following our renormalization discussion of Sec. V B, we stress that in the case of binary interactions, exact analogy with the semi-classical two-mode photon laser theory is obtained only within a $`T=0`$ Gross-Pitaevskii basis, i.e. a basis dressed only by the condensate mean field potentials via $`H_0^{^{}}=_{rn}\mathrm{\Xi }_{rn}^{Bare}\widehat{a}_r^{}\widehat{a}_n+(1/2)_{rsmn}T_{rsmn}z_r^{}z_s^{}z_mz_n`$, in analogy to the above-mentioned behavior of the photon laser.

Changing our description slightly to explicitly include the mean-field effects on eigenenergies of our single-particle system, we demonstrated that our treatment is consistent with the usual description in which condensed atoms are described in terms of the NLSE, and the evolution of non-condensed atoms is based on the quantum Boltzmann equation . In this case, the single-particle eigenenergies are ultimately effectively dressed by mean fields in the many-body T-matrix approximation, just as was found in the approach of Stoof  (with the many-body corrections arising from suitable inclusion of the anomalous average ). Our approach has been compared in more detail to the treatment of Walser et al , which has been recently shown  to be equivalent to the Kadanoff-Baym Green’s function formalism as applied to trapped Bose gases by Imamovic-Tomasovic and Griffin . In particular, we have discussed how additional contributions which arise in our second order collisional integrals upon choosing a simpler (than HFB) unperturbed basis, depend on the choice of this (bare or partially dressed) basis. This suggests an implicit assumption by Walser et al. that their eigenenergies are dressed by mean fields in the Hartree-Fock-Bogoliubov approximation, which is equivalent to their choice of slowly-evolving master variables. Our method thus yields an alternative microscopic derivation of the theory of Walser et al., based on a coupled equations of motion formalism.

When applying perturbation theory only to second order in the potential, one has to justify why such a truncation is physically realistic. This perturbation theory is essentially a systematic expansion in terms of the diluteness parameter $`\sqrt{na^3}`$. At $`T=0`$ such a treatment can be justified in the limit $`\sqrt{na^3}1`$ . For a large system leading to large mean fields which heavily dress the single-particle eigenenergies from their bare trap values, one conventionally shifts to a description in terms of quasiparticles, thus employing HFB-shifted eigenenergies. In this case, the validity criterion of perturbation theory at finite temperature essentially becomes $`(kT/nU_o)(\sqrt{na^3})1`$ for the homogeneous system , a criterion closely related to the one for the absence of critical fluctuations occuring sufficiently close to the transition point. On the other hand, the main part of this paper has been based on a treatment in terms of bare trap eigenenergies, since this enables a simple and direct analogy with the semi-classical equations for multi-mode photon laser theory. Clearly such a treatment cannot be valid for large, dense condensates, and its validity will be restricted to systems in which the interactions and particle numbers are so small, that the trap eigenenergies become only slighlty perturbed by the mean fields. A minimum (but not necessarily sufficient) criterion here is that the condensates are very weakly-interacting in the sense that $`nU_o\mathrm{}\omega `$. By identifying the additional basis-dependent corrections of our treatment to the second order collisional integrals of Walser et al, we can restate this condition as the requirement that the additional terms are much smaller than the first order contributions. Since the additional terms can be visualized as rates of change of the parameters $`z`$, $`\rho `$ and $`\kappa `$, this criterion essentially reduces to the slow evolution of such quantities. Hence, our bare basis analysis can only be useful in the limit when all mean field potentials evolve very slowly, and do not heavily modify the bare trap eigenenergies. However, our explicit expressions additionally include such slow evolution from a bare basis description, as opposed to existing theories in which such evolution is absent, applicable in the domain where the mean field energy shifts induced (on the already renormalized eigenenergies) during a collision can be neglected. Indeed, by re-formulating our treatment explicitly in terms of HFB eigenenergies, we find contributions of this type arising only in higher orders in the potential, being generated by the careful consideration of triplet and higher order averages, an issue which will be explicitly addressed elsewhere.

## ACKNOWLEDGMENTS

We acknowledge discussions with R. Walser and J. Williams. This work was partially supported by the ULF (Contract No. ERB-FMGECT 950021) and by a grant of the Max-Planck-Institute for Quantum Optics, Garching, where NPP originally developed part of this formalism.

## A Off-Diagonal Contributions

For completeness, we give here all off-diagonal contributions to the condensate / non-condensate populations of the three-level system discussed in the text. In particular, we have for the total population evolution

$`\stackrel{~}{R}_{T_{|0211|^2}}=`$ $`4\mathrm{\Gamma }_{0211}\left\{\begin{array}{c}(N_1N_2)\rho _{01}\rho _{10}+(N_1^TN_2)(\zeta _{01}\rho _{10}+\rho _{01}\zeta _{10})\hfill \\ +(N_1N_0)\rho _{12}\rho _{21}+(N_1^TN_0)(\zeta _{12}\rho _{21}+\rho _{12}\zeta _{21})\hfill \end{array}\right\}`$ (A5)
$`2\mathrm{\Gamma }_{0211}(1+2N_1)\left[\rho _{02}\rho _{20}+\zeta _{02}\rho _{20}+\rho _{02}\zeta _{20}\right]`$
$`8\mathrm{\Gamma }_{0211}\left\{\rho _{01}\rho _{12}\rho _{20}+\zeta _{01}\rho _{12}\rho _{20}+\rho _{01}\zeta _{12}\rho _{20}+\rho _{01}\rho _{12}\zeta _{20}\right\}`$
$`\stackrel{~}{Q}_{ijji}^{0211}=`$ $`\nu _{01}^{(2)}\left\{\begin{array}{c}\rho _{10}\left[4N_0\rho _{12}+4N_0^T\zeta _{12}+2\left(\rho _{10}\rho _{02}+\zeta _{10}\rho _{02}+\rho _{10}\zeta _{02}\right)\right]\hfill \\ +\zeta _{10}\left[4N_0^T\rho _{12}+2N_0^T\zeta _{12}+\left(2\rho _{10}\rho _{02}+\zeta _{10}\rho _{02}+2\rho _{10}\zeta _{02}\right)\right]\hfill \end{array}\right\}`$ (A13)
$`\nu _{21}^{(0)}\left\{\begin{array}{c}\rho _{12}\left[4N_2\rho _{10}+4N_2^T\zeta _{10}+2\left(\rho _{12}\rho _{20}+\zeta _{12}\rho _{20}+\rho _{12}\zeta _{20}\right)\right]\hfill \\ +\zeta _{12}\left[4N_2^T\rho _{10}+2N_2^T\zeta _{10}+\left(2\rho _{12}\rho _{20}+\zeta _{12}\rho _{20}+2\rho _{12}\zeta _{20}\right)\right]\hfill \end{array}\right\}`$
$`+6\vartheta _{01}^{(2)}\left[2N_1\rho _{10}\rho _{12}+2N_1^T\left(\zeta _{10}\rho _{12}+\rho _{10}\zeta _{12}\right)+N_1^{TT}\zeta _{10}\zeta _{12}\right]`$
$`+\epsilon _{01}^{(2)}\left[2\rho _{10}\rho _{12}+2\left(\zeta _{10}\rho _{12}+\rho _{10}\zeta _{12}\right)+\zeta _{10}\zeta _{12}\right]`$
Although $`\stackrel{~}{R}_{T_{|0211|^2}}`$ vanishes in the limit of diagonal $`\rho `$, the same does not apply to $`\stackrel{~}{Q}_{ijji}^{0211}`$ which contains a contribution $`\zeta _{10}\zeta _{12}`$, which will only vanish upon assuming that one of the three levels (e.g. level 2) is fully uncondensed.

In the above expressions, in addition to the total population $`N_i=n_i+|z_i|^2`$, we have defined the ‘population terms’ $`N_i^T=\frac{1}{2}\left(2n_i+|z_i|^2\right)`$ and $`N_i^{TT}=\frac{1}{3}\left(3n_i+|z_i|^2\right)`$. The population $`N_i^T`$ appears familiar from laser physics, as it contains the thermal (chaotic) contribution with a pre-factor of 2 over the corresponding condensed (ordered) term. Such a term will replace the full population term $`N_i`$ whenever the index $`i`$ is multiplied by a condensate amplitude $`z_i`$. By analogy, $`N_iN_i^{TT}`$ whenever the population term appears multiplied by two condensate mean amplitudes $`z_iz_i`$ (i.e. condensate anomalous averages of the form $`zz`$, or $`z^{}z^{}`$, but not normal averages $`z^{}z`$).

In order to keep the notation general, so that it can be easily applicable to n partially condensed levels with next neighbour interactions, we have also defined the following rates

$$\nu _{(i1)i}^{(i\pm 1)}=I_0\left[4V_{(i1)ii(i1)}2V_{(i1)(i\pm 1)(i\pm 1)(i1)}V_{(i1)(i1)(i1)(i1)}\right]V_{(i1)(i\pm 1)ii}$$
(A14)
$$\vartheta _{(i1)i}^{(i\pm 1)}=I_0\left[V_{(i1)ii(i1)}+V_{i(i\pm 1)(i\pm 1)i}V_{iiii}\right]V_{(i1)(i\pm 1)ii}$$
(A15)
$$\epsilon _{(i1)i}^{(i\pm 1)}=I_0\left[2V_{(i1)(i\pm 1)(i\pm 1)(i1)}V_{iiii}\right]V_{(i1)(i\pm 1)ii}$$
(A16)
$$\lambda _{ijji}^{kllk}=V_{ijji}V_{kllk}\underset{\eta 0^+}{lim}_0^{\mathrm{}}𝑑\tau e^{\eta \tau }$$
(A17)
$$\gamma _{ijji}^{0211}=I_0\left(V_{ijji}V_{0211}\right)$$
(A18)
where

$$I_0=\frac{1}{2}\left\{\underset{\eta 0^+}{lim}_0^{\mathrm{}}𝑑\tau e^{\pm i(\omega _0+\omega _22\omega _1\pm i\eta )\tau }+\underset{\eta 0^+}{lim}_0^{\mathrm{}}𝑑\tau e^{\eta \tau }\right\}$$
(A19)
For the condensed components of the two lowest states, we obtain the following off-diagonal contributions

$`\stackrel{~}{X}_{00}(V^2)=`$ $`4|z_0|^2\left\{\begin{array}{c}2\lambda _{0000}^{0110}\rho _{01}\rho _{10}+\left(\mathrm{\Gamma }_{0211}+2\lambda _{0110}^{0220}\right)\rho _{12}\rho _{21}+2\lambda _{0000}^{0220}\rho _{02}\rho _{20}\hfill \\ +2\left(2\gamma _{0110}^{0211}\gamma _{0220}^{0211}\right)\rho _{10}\rho _{12}\hfill \end{array}\right\}`$ (A31)
$`+4|z_1|^2\left[\lambda _{0110}\rho _{01}\rho _{10}+\mathrm{\Gamma }_{0211}\rho _{12}\rho _{21}+2\gamma _{0110}^{0211}\rho _{10}\rho _{12}\right]`$
$`+\zeta _{10}\left\{\begin{array}{c}2|z_0|^2\left[4\lambda _{0000}^{0110}\rho _{01}+\left(4\gamma _{0110}^{0211}2\gamma _{0220}^{0211}\gamma _{0000}^{0211}\right)\rho _{12}\right]\hfill \\ +2|z_1|^2\left[\mathrm{\Gamma }_{0211}\rho _{01}+\left(\gamma _{0110}^{0211}+3\gamma _{1221}^{0211}3\gamma _{1111}^{0211}\right)\rho _{12}\right]\hfill \\ 2\left[2\lambda _{0110}\rho _{01}+\left(2\gamma _{0220}^{0211}+\gamma _{1111}^{0211}\right)\rho _{12}\right]\hfill \\ +4\rho _{01}\left[\mathrm{\Gamma }_{0211}(n_1n_2)2\lambda _{0110}n_1\right]\hfill \\ 4\rho _{12}\left[\gamma _{0000}^{0211}n_0+\left(\gamma _{0110}^{0211}+3\gamma _{1111}^{0211}3\gamma _{1221}^{0211}\right)n_1+\left(4\gamma _{1221}^{0211}\gamma _{2222}^{0211}\right)n_2\right]\hfill \\ 4\left[\left(\lambda _{0211}+2\lambda _{0110}^{0220}\right)\rho _{02}\rho _{21}+\gamma _{0000}^{0211}\rho _{10}\rho _{02}\right]\hfill \\ 2\gamma _{0000}^{0211}\zeta _{10}\rho _{02}\hfill \end{array}\right\}`$
$`\stackrel{~}{X}_{11}(V^2)=`$ $`4|z_1|^2\left\{\begin{array}{c}\left(\mathrm{\Gamma }_{0211}+2\lambda _{1111}^{0110}\right)\rho _{01}\rho _{10}+\left(\mathrm{\Gamma }_{0211}+2\lambda _{1111}^{1221}\right)\rho _{12}\rho _{21}+\left(\mathrm{\Gamma }_{0211}+2\lambda _{0110}^{1221}\right)\rho _{02}\rho _{20}\hfill \\ +2\left(2\gamma _{0110}^{0211}+2\gamma _{1221}^{0211}\gamma _{1111}^{0211}\right)\rho _{10}\rho _{12}\hfill \end{array}\right\}`$ (A45)
$`+4|z_0|^2\left[\lambda _{0110}\rho _{01}\rho _{10}+\mathrm{\Gamma }_{0211}\rho _{12}\rho _{21}+2\gamma _{0110}^{0211}\rho _{10}\rho _{12}\right]`$
$`+\zeta _{10}\left\{\begin{array}{c}2|z_1|^2\left[\left(3\mathrm{\Gamma }_{0211}+4\lambda _{1111}^{0110}\right)\rho _{01}+\left(5\gamma _{0110}^{0211}+5\gamma _{1221}^{0211}3\gamma _{1111}^{0211}\right)\rho _{12}\right]\hfill \\ +2|z_0|^2\left[\left(4\gamma _{0110}^{0211}2\gamma _{0220}^{0211}\gamma _{0000}^{0211}\right)\rho _{12}\right]\hfill \\ 2\left[2\lambda _{0110}\rho _{01}+\left(2\gamma _{0220}^{0211}+\gamma _{1111}^{0211}\right)\rho _{12}\right]\hfill \\ +4\rho _{01}\left[\mathrm{\Gamma }_{0211}(n_1n_2)2\lambda _{0110}n_0\right]\hfill \\ 4\rho _{12}\left[\begin{array}{c}\left(\gamma _{0000}^{0211}+2\gamma _{0220}^{0211}2\gamma _{0110}^{0211}\right)n_0+\left(\gamma _{0110}^{0211}+\gamma _{1221}^{0211}\gamma _{1111}^{0211}\right)n_1\hfill \\ +\left(\gamma _{2222}^{0211}+2\gamma _{0220}^{0211}2\gamma _{1221}^{0211}\right)n_2\hfill \end{array}\right]\hfill \\ 4\left[\left(\mathrm{\Gamma }_{0211}+2\lambda _{0110}^{1221}\right)\rho _{02}\rho _{21}+\left(\gamma _{0000}^{0211}+2\gamma _{0220}^{0211}\right)\rho _{10}\rho _{02}\right]\hfill \\ 2\left(\gamma _{0000}^{0211}+2\gamma _{0220}^{0211}\right)\zeta _{10}\rho _{02}\hfill \end{array}\right\}`$
Although it is exteremely hard to discuss the physical implication of each separate contribution appearing in the above equations, at this point we would like to comment briefly on a carefully selected subset of the above equations, namely

$$\frac{d|z_0|^2}{dt}=8|z_1|^2\{\lambda _{0110}\rho _{01}\rho _{10}+\mathrm{\Gamma }_{0211}\rho _{12}\rho _{21}+2\gamma _{0110}^{0211}\rho _{10}\rho _{12}]$$
(A46)
$$\frac{d|z_1|^2}{dt}=8|z_0|^2\{\lambda _{0110}\rho _{01}\rho _{10}+\mathrm{\Gamma }_{0211}\rho _{12}\rho _{21}+2\gamma _{0110}^{0211}\rho _{10}\rho _{12}]$$
(A47)
Such terms indicate clearly a mechanism of ‘coherent population transfer’, that is growth of condensation in one (bare) level due to the existence of condensation in another (bare) level, which occurs via off-diagonal incoherent couplings $`\rho _{ij}`$. By further ignoring, for simplicity, the redistributional processes ($`V_{0211}0`$) we see a sort of Rabi-like oscillation between the two condensates, via the process

$$\frac{d|z_i|^2}{dt}|V_{ijji}|^2\left[\rho _{ij}\rho _{ji}\right]|z_j|^2$$
(A48)
Such processes do not lead to changes in the total trap level populations. On the contrary, due to the co-existence of coherent and incoherent atoms within the same trap level, this term can be interpreted as follows: If one level initially contains a non-vanishing coherent amplitude, this will be gradually transferred to the other level (even if the other level is initially fully incoherent) without an associated change in the total population of each level (since the incoherent population $`n_i`$ of each level adjusts accordingly to keep $`N_i`$ fixed). Hence, even if we assume that the condensate initially resides only on the bottom trap level, the above contribution generated by nonlinear coupling interactions will tend to give rise to a coherent mean field amplitude in level 1 (and vice versa). Whether this term actually becomes important at all (and in what limits this may be so) depends on how heavily this contribution is overshadowed by all other existing terms in the equations of motion for coherent evolution. Although we would not expect this term to play a significant role in the full population dynamics, this is something which should be confirmed by direct numerical simulation of the full equations given in this paper. We conclude the discussion by noting that such processes may become important in the case of coupled 2-species condensation , or in the creation of non-ground state trapped condensates , whereby the population between such levels is controlled by the application of external fields.

The off-diagonal contributions to the uncondensed dynamics of the two lowest levels are given by

$`Y_{00}(V^2)=`$ $`2\mathrm{\Gamma }_{0211}\left\{\rho _{02}\rho _{20}+2\left[\rho _{12}\rho _{21}n_0\left(\rho _{01}\rho _{10}+\rho _{12}\rho _{21}\rho _{02}\rho _{20}\right)n_1+\rho _{01}\rho _{10}n_2\right]\right\}`$ (A63)
$`2\left[4\mathrm{\Gamma }_{0211}\rho _{02}\rho _{21}\rho _{10}+\nu _{01}^{(2)}\rho _{10}^2\rho _{02}+\nu _{21}^{(0)}\rho _{12}^2\rho _{20}\right]`$
$`+8|z_0|^2\left[\lambda _{0000}^{0110}\rho _{01}\rho _{10}+\lambda _{0110}^{0220}\rho _{12}\rho _{21}+\lambda _{0000}^{0220}\rho _{02}\rho _{20}\right]`$
$`+4|z_1|^2\left[\left(\mathrm{\Gamma }_{0211}\lambda _{0110}\right)\rho _{01}\rho _{10}\mathrm{\Gamma }_{0211}\rho _{02}\rho _{20}\right]`$
$`+\rho _{10}\rho _{12}\left\{\begin{array}{c}2\left(\gamma _{0220}^{0211}\gamma _{1111}^{0211}\right)4\left[\nu _{01}^{(2)}n_03\vartheta _{01}^{(2)}n_1+\nu _{21}^{(0)}n_2\right]\hfill \\ +4|z_0|^2\gamma _{0000}^{0211}+4|z_1|^2\left(\gamma _{0110}^{0211}+3\gamma _{1221}^{0211}3\gamma _{1111}^{0211}\right)\hfill \end{array}\right\}`$
$`+\zeta _{10}\left\{\begin{array}{c}8\lambda _{0000}^{0110}|z_0|^2\rho _{01}\hfill \\ +2|z_1|^2\left[\mathrm{\Gamma }_{0211}\rho _{01}+2\gamma _{0110}^{0211}\rho _{12}\right]\hfill \\ +4\lambda _{0110}\rho _{01}+8\gamma _{0220}^{0211}\rho _{12}\hfill \\ +4\rho _{01}\left[\left(2\lambda _{0110}+\mathrm{\Gamma }_{0211}\right)n_1\mathrm{\Gamma }_{0211}n_2\right]\hfill \\ +4\rho _{12}\left[2\left(\gamma _{0220}^{0211}2\gamma _{0110}^{0211}\right)n_0+4\gamma _{0110}^{0211}n_1+2\gamma _{0220}^{0211}n_2\right]\hfill \\ +4\rho _{02}\rho _{21}\left(2\lambda _{0110}^{0220}\mathrm{\Gamma }_{0211}\right)+8\rho _{10}\rho _{02}\left(\gamma _{0000}^{0211}2\gamma _{0110}^{0211}+\gamma _{0220}^{0211}\right)\hfill \\ +\zeta _{10}\rho _{02}\left(3\gamma _{0000}^{0211}4\gamma _{0110}^{0211}+2\gamma _{0220}^{0211}\right)\hfill \end{array}\right\}`$
and

$`Y_{11}(V^2)=`$ $`4\mathrm{\Gamma }_{0211}\left\{\rho _{02}\rho _{20}+2\left[\rho _{12}\rho _{21}n_0\left(\rho _{01}\rho _{10}+\rho _{12}\rho _{21}\rho _{02}\rho _{20}\right)n_1+\rho _{01}\rho _{10}n_2\right]\right\}`$ (A78)
$`+4\left[4\mathrm{\Gamma }_{0211}\rho _{02}\rho _{21}\rho _{10}+\nu _{01}^{(2)}\rho _{10}^2\rho _{02}+\nu _{21}^{(0)}\rho _{12}^2\rho _{20}\right]`$
$`+4|z_0|^2\left[\lambda _{0110}\rho _{01}\rho _{10}+\mathrm{\Gamma }_{0211}\rho _{12}\rho _{21}\right]`$
$`+4|z_1|^2\left[\left(2\lambda _{0110}^{1111}\mathrm{\Gamma }_{0211}\right)\rho _{01}\rho _{10}+\left(2\lambda _{1111}^{1221}\mathrm{\Gamma }_{0211}\right)\rho _{12}\rho _{21}+\left(2\lambda _{0110}^{1221}+\mathrm{\Gamma }_{0211}\right)\rho _{02}\rho _{20}\right]`$
$`+\rho _{10}\rho _{12}\left\{\begin{array}{c}4\left(2\gamma _{0220}^{0211}\gamma _{1111}^{0211}\right)+8\left[\nu _{01}^{(2)}n_03\vartheta _{01}^{(2)}n_1+\nu _{21}^{(0)}n_2\right]\hfill \\ +8|z_0|^2\left(3\gamma _{0110}^{0211}2\gamma _{0220}^{0211}\gamma _{0000}^{0211}\right)+8|z_1|^2\left(2\gamma _{1111}^{0211}\gamma _{0110}^{0211}\gamma _{1221}^{0211}\right)\hfill \end{array}\right\}`$
$`+\zeta _{10}\left\{\begin{array}{c}2\nu _{01}^{(2)}|z_0|^2\rho _{12}\hfill \\ +2|z_1|^2\left[\left(4\lambda _{0110}^{1111}\mathrm{\Gamma }_{0211}\right)\rho _{01}+\left(3\gamma _{1111}^{0211}\gamma _{0000}^{0211}\gamma _{1221}^{0211}\right)\rho _{12}\right]\hfill \\ +4\lambda _{0110}\rho _{01}+2\left(3\gamma _{1111}^{0211}2\gamma _{0220}^{0211}\right)\rho _{12}\hfill \\ +4\rho _{01}\left[2\lambda _{0110}n_03\mathrm{\Gamma }_{0211}n_1+3\mathrm{\Gamma }_{0211}n_2\right]\hfill \\ +4\rho _{12}\left[\left(6\gamma _{0110}^{0211}2\gamma _{0220}^{0211}\gamma _{0000}^{0211}\right)n_05\vartheta _{01}^{(2)}n_1+\left(6\gamma _{1221}^{0211}2\gamma _{0220}^{0211}\gamma _{2222}^{0211}\right)n_2\right]\hfill \\ +4\rho _{02}\rho _{21}\left(2\lambda _{0110}^{1221}+3\mathrm{\Gamma }_{0211}\right)+4\rho _{10}\rho _{02}\left(8\gamma _{0110}^{0211}2\gamma _{0220}^{0211}\gamma _{0000}^{0211}\right)\hfill \\ +8\zeta _{10}\rho _{02}\gamma _{0110}^{0211}\hfill \end{array}\right\}`$
## B Derivation of Condensate Evolution

The second order contribution to $`|z_i|^2`$ of Eq. (17) can be found by adiabatically eliminating the quantity $`z_i^{}\widehat{a}_s^{}\widehat{a}_m\widehat{a}_r`$ via

$$i\frac{dz_i}{dt}=\omega _iz_i+\underset{pql}{\overset{^{}}{}}T_{ipql}\widehat{a}_p^{}\widehat{a}_q\widehat{a}_l$$
(B1)
and

$$i\frac{d}{dt}\widehat{a}_s^{}\widehat{a}_m\widehat{a}_r=(\omega _m+\omega _r\omega _s)\widehat{a}_s^{}\widehat{a}_m\widehat{a}_r+f_{\widehat{a}^{}\widehat{a}\widehat{a}}$$
(B2)
with

$`f_{\widehat{a}^{}\widehat{a}\widehat{a}}`$ $`={\displaystyle \underset{lt}{}}T_{mrlt}\widehat{a}_s^{}\widehat{a}_l\widehat{a}_t`$ (B4)
$`+{\displaystyle \underset{plt}{}}T_{prlt}\widehat{a}_p^{}\widehat{a}_s^{}\widehat{a}_m\widehat{a}_l\widehat{a}_t+{\displaystyle \underset{plt}{}}T_{pmlt}\widehat{a}_p^{}\widehat{a}_s^{}\widehat{a}_r\widehat{a}_l\widehat{a}_t{\displaystyle \underset{pql}{}}T_{pqls}\widehat{a}_p^{}\widehat{a}_q^{}\widehat{a}_l\widehat{a}_m\widehat{a}_r`$
Hence, the quantity $`F_{\widehat{a}^{}\widehat{a}\widehat{a}}`$ corresponding to $`F`$ of Eq. (8), now becomes

$$F_{\widehat{a}^{}\widehat{a}\widehat{a}}=z_i^{}f_{\widehat{a}^{}\widehat{a}\widehat{a}}\underset{pql}{}T_{pqli}\widehat{a}_p^{}\widehat{a}_q^{}\widehat{a}_l\widehat{a}_s^{}\widehat{a}_m\widehat{a}_r$$
(B5)
## C Second Order Collisional Integrals

For completeness, we give here the second order collisional integrals of Walser et al.  as generated by our approach.

$`\left({\displaystyle \frac{d\rho }{dt}}\right)^{Walser}`$ $`=`$ $`2{\displaystyle \underset{rsmn}{}}{\displaystyle \underset{pqlt}{}}T_{rsmn}T_{pqlt}^{(\delta )}`$ (C1)
$`\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 \\ +\left[\left(\rho _{mp}+\delta _{mp}\right)\left(\rho _{nq}+\delta _{nq}\right)\left(z_s^{}z_t\right)\rho _{lr}\rho _{mp}\rho _{rq}\left(z_s^{}z_t\right)\left(\rho _{lr}+\delta _{lr}\right)\right]\hfill \\ +2\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 \\ +2\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 \\ +2\left[\left(z_p^{}z_m\right)\kappa _{nt}\kappa _{qs}^{}\rho _{lr}\left(z_p^{}z_m\right)\kappa _{rt}\kappa _{qs}^{}\left(\rho _{lr}+\delta _{lr}\right)\right]\hfill \\ +2\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 \\ +2\kappa _{qi}^{}\left[\begin{array}{c}\left[\left(\rho _{mp}+\delta _{mp}\right)\kappa _{rl}\rho _{ts}\rho _{mp}\kappa _{rl}\left(\rho _{ts}+\delta _{ts}\right)\right]\hfill \\ +\left[\left(z_p^{}z_m\right)\kappa _{rl}\rho _{ts}\left(z_p^{}z_m\right)\kappa _{rl}\left(\rho _{ts}+\delta _{ts}\right)\right]\hfill \\ +\left[\left(\rho _{mp}+\delta _{mp}\right)\kappa _{rl}\left(z_s^{}z_t\right)\rho _{mp}\kappa _{rl}\left(z_s^{}z_t\right)\right]\hfill \\ +\left[\left(\rho _{mp}+\delta _{mp}\right)\left(z_rz_l\right)\rho _{ts}\rho _{mp}\left(z_rz_l\right)\left(\rho _{ts}+\delta _{ts}\right)\right]\hfill \end{array}\right]\hfill \end{array}\right\}`$ (C15)
As explained in the text

$$T_{pqlt}^{(\delta )}=𝑑t^{^{}}e^{i\left(\omega _l+\omega _t\omega _p\omega _q\right)(tt^{^{}})}=\pi \delta (\mathrm{\Delta }\omega )iP\left(\frac{1}{\mathrm{\Delta }\omega }\right)$$
(C16)
Here $`\mathrm{\Delta }\omega =\left(\omega _l+\omega _t\omega _p\omega _q\right)`$ and the $`\omega _i`$ denotes the eigenenergy of level $`i`$ in the particular basis chosen for the analysis of the system. In our original formulation, these correspond to bare trap eigenenergies, whereas in order to establish exact analogy with the collisional integrals of Walser et al. , these should be replaced by their corresponding values dressed by normal and anomalous HFB mean fields.

In the original formulation of Walser et al., their basis is left unspecified, and their resulting expressions are assumed to be valid for any single-particle basis dressed by, at most, number-conserving mean fields (i.e. no anomalous averages), a conclusion not supported by our analysis.