text:
# Limitations of squeezing due to collisional decoherence in Bose-Einstein condensates

## I Introduction

Squeezing and entanglement in coherent many particle systems have become important issues in experiments with Bose-Einstein condensates (BECs) as well as in the theoretical work related to them. Achieving squeezed and/or entangled states of a large number of ultracold atoms would have an enormous impact on the research field which includes quantum information technology (QIT) , the development of new time standards  and the test of quantum state reduction . In particular, whereas techniques for entangling photons  are already available for the transmission of quantum information , schemes for the local storage and processing are still missing. A promising step in the desired direction represent the experiments which allow for photon pulses being slowed and even stored for a considerably long time in cold and ultracold, condensed atomic gases  (for a recent review cf. ). These experiments are based on electromagnetically induced transparancy  and use the possibility of dark states, which for the case of many-particle systems involving a condensate may be described as entangled states of atoms and photons.

In  we have studied the limitations of relative number squeezing in entangled states of atoms and photons produced by Bragg-scattering of light off a condensate. Entanglement between photons and atoms has been proposed and studied by different authors before . We took into account interactions of the condensate atoms with non-condensed atoms to the degree that the quasiparticles resulting from the diagonalization of the quadratic part of the Hamiltonian remained free. We have shown that the amount of squeezing in such systems underlies important limits imposed by the collisional interactions of the atoms. In this paper we would like to extend this work and study the effect of quasiparticle interactions which lead to collisional damping of excitations. We will find additional constraints on the squeezing considered in .

Moreover, we will consider the case where the non-equilibrium matter wave excitations consist of atoms which are in a different internal state than the condensate atoms. This is e.g. relevant in the case of a dark state involving a condensate, a classical coupling laser and quantized probe and atomic modes, which allows to slow and stop light pulses in atomic matter . The quantized atomic mode is in general in a different internal state than the classical (condensate) mode and the decoherence time is governed by the spin-conserving collisions between the atoms of the different modes. In this paper we derive the collisional damping rate of an excited atomic mode for such a system with two internal degrees of freedom, in an analogous way as in the single level case. In  we have studied the limitations due to this collisional damping of the light delay and storage times in experiments based on electromagnetically induced transparency.

Mean field theory and its extensions have been widely used to describe stationary properties of homogeneous and trapped dilute, weakly interacting Bose-Einstein condensates at finite temperatures . In these approaches, the non-condensate component is treated as a static thermal bath and its dynamical coupling to the condensate oscillations is neglected. This implies that damping of out-of-equilibrium-excitations is not accounted for. Time dependent mean field approaches have been introduced  to describe the damping in trapped condensates at finite temperatures . The calculated damping rates agree with the results from perturbation theory . In  damping rates in trapped condensates at finite temperature have been calculated explicitly. In order to describe the dynamics of a decaying excitation consistently at large times it is necessary to take into account the connection between dissipation and fluctuations. To this end it is convenient to use the Heisenberg-Langevin approach (cf. e.g. ) to treat the time development of the system . In this paper we would like to recall this procedure as applied to the dissipation dynamics of condensates with atoms in a single internal state and extend it to the treatment of the two internal levels case. Moreover, the Heisenberg-Langevin equations will enable us to derive the time evolution of expectation values of products of two and more quasiparticle operators which we need to study the squeezing described above.

In the basic, quadratic approximation to the Hamiltonian of the system, including a quartic collisional contact interaction, mean field theory yields the ground state wave function and the spectrum of elementary excitations, which may be understood as non-interacting quasiparticles. In this theory, the Heisenberg equations of motion of the Fock operators assume the form

$$\dot{\widehat{\beta }}_i(t)=i\omega _i\widehat{\beta }_i(t),$$
(1)
where $`\widehat{\beta }_i(t)`$ is the annihilation operator of the quasiparticle mode labelled by $`i`$, which satisfies Bose equal-time commutation relations with its hermitian conjugate, $`[\widehat{\beta }_i(t),\widehat{\beta }_j^{}(t)]=\delta _{ij}`$. $`\omega _i`$ is the corresponding eigenfrequency.

Beyond this free field theory the atomic collisions induce interactions of the quasiparticles. The purpose of this paper is to consider the impact of these interactions on the dynamic behaviour of a BEC system which is initially in a state out of thermal equilibrium. Specifically we study the decay behaviour of a single quasiparticle mode whose initial population is different from the equilibrium one. As this quasiparticle mode is interacting with the bath of other modes the system will evolve back to a stationary equilibrium state. The dynamics of such processes may be described using a time-dependent mean field approach based on the finite temperature extension of the Gross-Pitaevski equation . This allows to describe oscillations and damping of the condensate mean field which initially is assumed to be slightly displaced from its stationary form.

As an alternative to the mean-field approach one may adhere for the first to operators and find the equations of motion for the quasiparticle operators. In the basis of the non-interacting quasiparticles of the quadratic theory these will be in general coupled non-linear equations. For the physical situations we are considering here, we are interested in approximations which allow to write the equation for the excited mode $`q`$ in the form

$$\dot{\widehat{\beta }}_q(t)=i(\omega _q+\delta \omega _q\frac{i}{2}\gamma _q)\widehat{\beta }_q(t),$$
(2)
where $`\delta \omega _q`$ is an energy shift and $`\gamma _q`$ a decay width. As was shown using the extended time dependent mean-field theory , the decay width arises from the coupling of the excited mode to the bath of the other modes which are assumed to be in thermal equilibrium and to have a broad range of suffiently closely spaced frequencies. Eqn. (2) may be used e.g. to calculate the time dependence of the mean population of mode $`q`$, which results as an exponential decay with rate $`\gamma _q`$. It is also clear, however, that the operator which evolves in time as described by (2) does not obey the equal time commutation relation, but commutes with $`\widehat{\beta }_q^{}(t)`$ for $`t\mathrm{}`$. Moreover, according to (2) the population of the mode $`q`$ decays to zero, i.e. not to the correct equilibrium value at finite temperatures. These problems are absent in the mean-field approach, where only the deviation from the stationary condensate field decays and where no operator ordering problems occur.

It is well known in quantum optics that these problems are resolved by the Heisenberg-Langevin theory of dissipative dynamics. Let us briefly summarize its results. The equation of motion (2) for $`\widehat{\beta }_q`$ contains a further term. In the Markov approximation, the equation may be written in the form

$$\dot{\widehat{\beta }_q}(t)=i(\omega _q+\delta \omega _q\frac{i}{2}\gamma _q)\widehat{\beta }_q(t)\widehat{F}_q(t),$$
(3)
where $`\widehat{F}_q`$ is a Langevin quantum noise operator which has a zero mean value. The Langevin operator in the Markov approximation obeys the commutation relation

$$[\widehat{F}_q(t),\widehat{F}_q^{}(t^{})]=\gamma _q\delta (tt^{}).$$
(4)
This ensures that the equal time commutation relation for $`\widehat{\beta }_q(t)`$ is conserved in time,

$$[\widehat{\beta }_q(t),\widehat{\beta }_q^{}(t)]=1.$$
(5)
In the Markov approximation the noise operator is delta-correlated in time, i.e. the two-time correlation function is given by

$$\widehat{F}_q^{}(t)\widehat{F}_q(t^{})=\gamma _qn_q^0\delta (tt^{}),$$
(6)
where $`n_q^0`$ is the equilibrium value of the mode population. Eqn. (6) means, that the random forces described by $`\widehat{F}_q`$ have a short memory. The equation of motion for the quasiparticle number operator $`\widehat{n}_q=\widehat{\beta }_q^{}\widehat{\beta }_q`$ becomes

$$\dot{\widehat{n}_q}(t)=\gamma _q\left[\widehat{n}_q(t)n_q^0\right]\widehat{F}_q^{}(t),$$
(7)
where $`\widehat{F}_q^{}`$ is the quantum noise operator defined by

$$\widehat{F}_q^{}\widehat{\beta }_q^{}\widehat{F}_q+\widehat{F}_q^{}\widehat{\beta }_q\widehat{\beta }_q^{}\widehat{F}_q+\widehat{F}_q^{}\widehat{\beta }_q,$$
(8)
which by definition has a zero mean value. To prove eqn. (7), it has to be shown expicitly that $`\widehat{\beta }_q^{}\widehat{F}_q+\widehat{F}_q^{}\widehat{\beta }_q=n_q^0`$. Eqn. (7) shows that the quasiparticle population of mode $`q`$ decays, as expected, exponentially to the equilibrium value $`n_q^0`$. Before the population can decay to zero the fluctuations introduced by the coupling to the bath modes force it to settle at the thermal equilibrium value.

Our paper is organized as follows: In section II we review the Heisenberg-Langevin approach to quasiparticle damping for the single internal level case. We use a linear response treatment of the operator time evolution analogous to the one used in  for the mean field fluctuation and derive both (4) and (6) in the Markov approximation. In section III we extend the approach to the case of two internal levels and derive the Beliaev- and Landau-type decay widths as well as the corresponding noise operators. In section IV we study time evolution equations of the form (7) for normal ordered products of two and more quasiparticle operators. We introduce a quantum regression theorem which allows to use the basic equation (2) for the operator to calculate the time evolution of the operator products. We use this to calculate the implications of quasiparticle damping for the time evolution of the atom-photon squeezing produced by Bragg scattering of photons off condensates as described in . Section V contains our conclusions.

## II Heisenberg-Langevin equation for systems of atoms in a single internal state

The Hamiltonian of the interacting system of non-relativistic bosonic atoms in the presence of an external trapping potential $`V(𝐫)`$, in the grand canonical description, has the form

$`\widehat{K}\widehat{H}\mu \widehat{N}`$ $`=`$ $`{\displaystyle d^3r\left\{\widehat{\psi }^{}(𝐫,t)\left(\frac{\mathrm{}^2\mathbf{}^2}{2m}+V(𝐫)\mu \right)\widehat{\psi }(𝐫,t)+\frac{g}{2}\widehat{\psi }^{}(𝐫,t)\widehat{\psi }^{}(𝐫,t)\widehat{\psi }(𝐫,t)\widehat{\psi }(𝐫,t)\right\}}.`$ (9)
We have chosen the contact potential approximation for the interatomic potential, $`U(𝐫)=g\delta (𝐫)`$, where $`g`$ is given in terms of the s-wave scattering length by $`g=4\pi \mathrm{}^2a/m`$. The field operators $`\widehat{\psi }(𝐫,t)`$ and $`\widehat{\psi }^{}(𝐫,t)`$ obey the usual equal-time commutation relations. Using these, the equation of motion for the field operator follows to be

$`i\mathrm{}{\displaystyle \frac{}{t}}\widehat{\psi }(𝐫,t)`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{}^2\mathbf{}^2}{2m}}+V(𝐫)\mu \right)\widehat{\psi }(𝐫,t)+g\widehat{\psi }^{}(𝐫,t)\widehat{\psi }(𝐫,t)\widehat{\psi }(𝐫,t).`$ (10)
In the usual treatment of the time evolution of a Bose condensate away from, but close to equilibrium a time dependent condensate wave function is defined by

$$\mathrm{\Phi }(𝐫,t)=\widehat{\psi }(𝐫,t),$$
(11)
where the avarage $`\mathrm{}`$ is with respect to the non-equilibrium initial state of the system. The field operator is then decomposed into a condensate and a non-condensate component,

$$\widehat{\psi }(𝐫,t)=\widehat{\mathrm{\Phi }}(𝐫,t)+\stackrel{~}{\psi }(𝐫,t),$$
(12)
and the non-condensate part is defined by $`\stackrel{~}{\psi }(𝐫,t)=0`$ . The condensate part is usually taken to be a c-number, but we will keep its operator character as far as its non-equilibrium contribution is concerned:

$$\widehat{\mathrm{\Phi }}(𝐫,t)=\mathrm{\Phi }_0(𝐫)+\widehat{\varphi }(𝐫,t).$$
(13)
Here $`\mathrm{\Phi }_0(𝐫)=\widehat{\psi }(𝐫)_0`$ is the stationary expectation value w.r.t. the equilibrium state of the BEC and assumed to be real. In the usual mean field treatment  the condensate part $`\widehat{\mathrm{\Phi }}`$ is replaced by its non-equilibrium expectation value $`\mathrm{\Phi }=\widehat{\mathrm{\Phi }}`$, and therefore $`\widehat{\varphi }`$ by the fluctuation field $`\delta \mathrm{\Phi }=\widehat{\varphi }`$, which is assumed to be small compared to $`\mathrm{\Phi }_0`$. In our further procedure we will also assume that the non-equilibrium expectation value of $`\widehat{\varphi }`$ remains small compared to $`\mathrm{\Phi }_0`$, in order to justify an analogous linear response approach to the time-evolution of the field operator $`\widehat{\varphi }`$.

We now define the fluctuation operators of the normal and anomalous densities as

$`\delta \widehat{\stackrel{~}{n}}`$ $`=`$ $`\stackrel{~}{\psi }^{}\stackrel{~}{\psi }\stackrel{~}{n}^0,`$ (14)
$`\delta \widehat{\stackrel{~}{m}}`$ $`=`$ $`\stackrel{~}{\psi }\stackrel{~}{\psi }\stackrel{~}{m}^0,`$ (15)
where we have, as we will do in the following, suppressed the space and time arguments of the operators and densities. The non-equilibrium expectation values of the operators (14,15) are the fluctuation functions $`\delta \stackrel{~}{n}`$, $`\delta \stackrel{~}{m}`$ defined in , which describe the deviation from the equilibrium expectation values $`\stackrel{~}{n}^0=\stackrel{~}{\psi }^{}\stackrel{~}{\psi }_0`$, $`\stackrel{~}{m}^0=\stackrel{~}{\psi }\stackrel{~}{\psi }_0`$.

Inserting eqns. (12)–(15) into (10) we obtain, to linear order in $`\widehat{\varphi }`$, $`\stackrel{~}{\psi }`$:

$`i\mathrm{}{\displaystyle \frac{}{t}}\left(\stackrel{~}{\psi }+\widehat{\varphi }\right)`$ $`=`$ $`\left(\widehat{K}_0+2g\left[n_0+\stackrel{~}{n}^0\right]\right)\stackrel{~}{\psi }+g\left[n_0+\stackrel{~}{m}^0\right]\stackrel{~}{\psi }^{}`$ (18)
$`+\left(\widehat{K}_0+2g\left[n_0+\stackrel{~}{n}^0+\mathrm{\Phi }_0\stackrel{~}{\psi }^{}+\mathrm{\Phi }_0\stackrel{~}{\psi }\right]\right)\widehat{\varphi }+g\widehat{\varphi }^{}\left[n_0+\stackrel{~}{m}^0+2\mathrm{\Phi }_0\stackrel{~}{\psi }\right]`$
$`+g\mathrm{\Phi }_0\left[2\delta \widehat{\stackrel{~}{n}}+\delta \widehat{\stackrel{~}{m}}\right],`$
where we have defined $`\widehat{K}_0=\mathrm{}^2\mathbf{}^2/(2m)+V\mu `$, and where $`n_0=|\mathrm{\Phi }_0|^2`$ is the condensate density. In deriving (18) we have also used that the equilibrium condensate wave function $`\mathrm{\Phi }_0`$ satisfies the generalized stationary Gross-Pitaevskii equation

$$\left(\widehat{K}_0+g\left[n_0+2\stackrel{~}{n}^0+\stackrel{~}{m}^0\right]\right)\mathrm{\Phi }_0=0,$$
(19)
and we have neglected the fluctuation operators

$$\delta (\stackrel{~}{\psi }^{}\stackrel{~}{\psi }\stackrel{~}{\psi })\stackrel{~}{\psi }^{}\stackrel{~}{\psi }\stackrel{~}{\psi }2\stackrel{~}{n}^0\stackrel{~}{\psi }\stackrel{~}{m}^0\stackrel{~}{\psi }^{}.$$
(20)
In the case, where the system is in equilibrium, the r.h.s. (18) reduces to its first line. Then, eqn. (18) may be solved by introducing the Bogoliubov expansion

$$\stackrel{~}{\psi }(𝐫,t)=\underset{i}{}\left[u_i(𝐫)\widehat{\beta }_i(t)+v_i^{}(𝐫)\widehat{\beta }_i^{}(t)\right]$$
(21)
of the field operator in terms of quasiparticle mode operators. In equilibrium, (18) is linear in the field operators, one has $`\widehat{\beta }_i(t)=\widehat{\beta }_i(t_0)\mathrm{exp}(i\omega _i(tt_0))`$ and the mode functions $`u_i`$, $`v_i`$ are the solutions of the Bogoliubov-de Gennes equations

$$\left(\begin{array}{cc}(𝐫)& (𝐫)\\ ^{}(𝐫)& ^{}(𝐫)\end{array}\right)\left(\begin{array}{c}u_i(𝐫)\\ v_i(𝐫)\end{array}\right)=\mathrm{}\omega _i\left(\begin{array}{c}u_i(𝐫)\\ v_i(𝐫)\end{array}\right),$$
(22)
where $`(𝐫)=\widehat{K}_0(𝐫)+2g[n_0(𝐫)+\stackrel{~}{n}^0(𝐫)]`$, $`(𝐫)=g[n_0(𝐫)+\stackrel{~}{m}^0(𝐫)]`$. A normalization of the mode functions $`u_i`$, $`v_i`$ by

$$d^3r\left[u_i^{}(𝐫)u_j(𝐫)v_i^{}(𝐫)v_j(𝐫)\right]=\delta _{ij}$$
(23)
ensures that the quasiparticle operators $`\widehat{\beta }_i`$, $`\widehat{\beta }_i^{}`$ satisfy bosonic commutation relations.

Using (21) we can expand the fluctuation operators (14), (15) in terms of the operators

$`\widehat{f}_{ij}`$ $`=`$ $`\widehat{\beta }_i^{}\widehat{\beta }_j\delta _{ij}n_i^0,`$ (24)
$`\widehat{g}_{ij}`$ $`=`$ $`\widehat{\beta }_i\widehat{\beta }_j`$ (25)
(the notation using $`f`$ and $`g`$ follows ), and find

$`\delta \widehat{\stackrel{~}{n}}`$ $`=`$ $`{\displaystyle \underset{ij}{}}\left\{\left[u_i^{}u_j+v_i^{}v_j\right]\widehat{f}_{ij}+u_iv_j\widehat{g}_{ij}+u_i^{}v_j^{}\widehat{g}_{ij}^{}\right\},`$ (26)
$`\delta \widehat{\stackrel{~}{m}}`$ $`=`$ $`{\displaystyle \underset{ij}{}}\left\{2v_i^{}u_j\widehat{f}_{ij}+u_iu_j\widehat{g}_{ij}+v_i^{}v_j^{}\widehat{g}_{ij}^{}\right\}.`$ (27)
We now proceed with solving eqn. (18) in the non-equilibrium case. For simplicity, we assume that the initial state is, in the basis of the quasiparticle modes defined by (21), a quasiparticle number state. We assume that a single quasiparticle mode is populated much stronger than the thermal equilibrium population, $`n_q(t=0)=\widehat{\beta }_q^{}(0)\widehat{\beta }_q(0)n_i^0=\widehat{\beta }_i^{}\widehat{\beta }_i_0=[\mathrm{exp}(\mathrm{}\omega _i/k_BT)1]^1`$ ($`iq`$).

We write the fluctuation operator $`\widehat{\varphi }`$ in the form

$$\widehat{\varphi }(𝐫,t)=u_q(𝐫)\widehat{\beta }_q(t)+v_q^{}(𝐫)\widehat{\beta }_q^{}(t).$$
(28)
In order to justify the linear response treatment, the particle population is still assumed to be much lower than that of the condensate mode, $`\widehat{\varphi }^{}(𝐫,0)\widehat{\varphi }(𝐫,0)=d^3r([|u_q(𝐫)|^2+|v_q(𝐫)|^2]n_q(t)+|v_q(𝐫)|^2)n_0`$.

Note that in (28) we assume $`\widehat{\psi }`$ to be defined by (21), but with the sum over $`i`$ excluding $`q`$. Following (14) and (15) the sums over $`i`$ and $`j`$ in the expansions (26) and (27) then also exclude $`q`$. However, we will include the terms proportional to $`\stackrel{~}{\psi }\widehat{\varphi }`$, $`\stackrel{~}{\psi }^{}\widehat{\varphi }`$ and $`\widehat{\varphi }^{}\stackrel{~}{\psi }`$ in (18) into the term $`g\mathrm{\Phi }_0[2\delta \widehat{\stackrel{~}{n}}+\delta \widehat{\stackrel{~}{m}}]`$ which means that in (26), (27) the sums over $`i`$ and $`j`$ also include $`q`$, except the combination $`i=j=q`$, which is neglected in linear response.

Now, inserting (21) and (28) in (18) and using (22) and (23), we find the equation of motion for $`\widehat{\beta }_q(t)`$:

$`\dot{\widehat{\beta }_q}(t)`$ $`=`$ $`i\omega _q\widehat{\beta }_q(t){\displaystyle \frac{i}{\mathrm{}}}g{\displaystyle \underset{ij}{}}\left\{2A_{q,ij}^{}\widehat{f}_{ij}(t)+B_{q,ij}^{}\widehat{g}_{ij}(t)+\stackrel{~}{B}_{q,ij}^{}\widehat{g}_{ij}^{}(t)\right\},`$ (29)
where

$`A_{q,ij}`$ $`=`$ $`{\displaystyle d^3r\mathrm{\Phi }_0\left[u_qA_{ij}^{}+v_qA_{ji}\right]},`$ (30)
$`B_{q,ij}`$ $`=`$ $`{\displaystyle d^3r\mathrm{\Phi }_0\left[u_qB_{ij}^{}+v_qC_{ij}^{}\right]},`$ (31)
$`\stackrel{~}{B}_{q,ij}`$ $`=`$ $`{\displaystyle d^3r\mathrm{\Phi }_0\left[u_qC_{ij}+v_qB_{ij}\right]}.`$ (32)
and

$`A_{ij}`$ $`=`$ $`u_i^{}u_j+v_i^{}v_j+v_i^{}u_j,`$ (33)
$`B_{ij}`$ $`=`$ $`u_iu_j+u_iv_j+v_iu_j,`$ (34)
$`C_{ij}`$ $`=`$ $`v_iv_j+u_iv_j+v_iu_j.`$ (35)
The time evolution equations for the density and pair operators $`\widehat{f}_{ij}`$, $`\widehat{g}_{ij}`$ may be derived in an analogous way as in . Inserting (12) and (13) into (9), using (14), (15), (20), (21), and (22), and keeping only terms up to linear order in the fluctuation operators, one finds:

$`\widehat{K}`$ $`=`$ $`{\displaystyle \underset{i}{}}\mathrm{}\omega _i\widehat{\beta }_i^{}\widehat{\beta }_i`$ (40)
$`+g{\displaystyle }d^3r\mathrm{\Phi }_0[2\widehat{\varphi }^{}\stackrel{~}{\psi }^{}\stackrel{~}{\psi }+\widehat{\varphi }^{}\stackrel{~}{\psi }\stackrel{~}{\psi }+\mathrm{h}.\mathrm{c}.]`$
$`+g{\displaystyle d^3r\left[4\delta \stackrel{~}{n}\stackrel{~}{\psi }^{}\stackrel{~}{\psi }+\delta \stackrel{~}{m}\stackrel{~}{\psi }^{}\stackrel{~}{\psi }^{}+\delta \stackrel{~}{m}^{}\stackrel{~}{\psi }\stackrel{~}{\psi }\right]}`$
$`+g{\displaystyle d^3r\left[\delta (\stackrel{~}{\psi }^{}\stackrel{~}{\psi }^{}\stackrel{~}{\psi }\stackrel{~}{\psi })+\mathrm{\Phi }_0\left(\delta (\stackrel{~}{\psi }^{}\stackrel{~}{\psi }^{}\stackrel{~}{\psi })+\delta (\stackrel{~}{\psi }^{}\stackrel{~}{\psi }\stackrel{~}{\psi })\right)\right]}`$
$`g{\displaystyle }d^3r\mathrm{\Phi }_0[(2\stackrel{~}{n}^0+\stackrel{~}{m}^0)\widehat{\varphi }+\mathrm{h}.\mathrm{c}.],`$
where

$`\delta (\stackrel{~}{\psi }^{}\stackrel{~}{\psi }^{}\stackrel{~}{\psi }\stackrel{~}{\psi })`$ $`=`$ $`\stackrel{~}{\psi }^{}\stackrel{~}{\psi }^{}\stackrel{~}{\psi }\stackrel{~}{\psi }4\stackrel{~}{\psi }^{}\stackrel{~}{\psi }\stackrel{~}{\psi }^{}\stackrel{~}{\psi }\stackrel{~}{\psi }^{}\stackrel{~}{\psi }^{}\stackrel{~}{\psi }\stackrel{~}{\psi }\stackrel{~}{\psi }\stackrel{~}{\psi }\stackrel{~}{\psi }^{}\stackrel{~}{\psi }^{},`$ (41)
$`\delta (\stackrel{~}{\psi }^{}\stackrel{~}{\psi }^{}\stackrel{~}{\psi })`$ $`=`$ $`\stackrel{~}{\psi }^{}\stackrel{~}{\psi }^{}\stackrel{~}{\psi }2\stackrel{~}{n}^0\stackrel{~}{\psi }^{}\stackrel{~}{m}^0\stackrel{~}{\psi },`$ (42)
$`\delta \stackrel{~}{n}`$ $`=`$ $`\stackrel{~}{\psi }^{}\stackrel{~}{\psi }\stackrel{~}{n}^0=\delta \widehat{\stackrel{~}{n}},`$ (43)
$`\delta \stackrel{~}{m}`$ $`=`$ $`\stackrel{~}{\psi }\stackrel{~}{\psi }\stackrel{~}{m}^0=\delta \widehat{\stackrel{~}{m}}.`$ (44)
Then, using (21), (24), (25) the equations of motion for $`\widehat{f}_{ij}`$ and $`\widehat{g}_{ij}`$ read ($`i,jq`$)

$`{\displaystyle \frac{}{t}}\widehat{f}_{ij}`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}[\widehat{K},\widehat{\beta }_i^{}\widehat{\beta }_j]`$ (45)
$`=`$ $`i(\omega _j\omega _i)\widehat{f}_{ij}2{\displaystyle \frac{i}{\mathrm{}}}g(n_i^0n_j^0){\displaystyle d^3r\mathrm{\Phi }_0\left[A_{ij}^{}\widehat{\varphi }+A_{ji}\widehat{\varphi }^{}\right]},`$ (46)
$`{\displaystyle \frac{}{t}}\widehat{g}_{ij}`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}[\widehat{K},\widehat{\beta }_i\widehat{\beta }_j]`$ (47)
$`=`$ $`i(\omega _j+\omega _i)\widehat{g}_{ij}2{\displaystyle \frac{i}{\mathrm{}}}g(1+n_i^0+n_j^0){\displaystyle d^3r\mathrm{\Phi }_0\left[B_{ij}^{}\widehat{\varphi }+C_{ij}^{}\widehat{\varphi }^{}\right]},`$ (48)
In deriving the time evolution (45,47) of the fluctuation operators from (40), the terms proportional $`\delta \stackrel{~}{n}`$ and $`\delta \stackrel{~}{m}`$ which describe the coupling to the fluctuations of the bath modes have been neglected. Also the terms arising from the fluctuations $`\delta (\stackrel{~}{\psi }^{}\stackrel{~}{\psi }^{}\stackrel{~}{\psi }^{}\stackrel{~}{\psi })`$, $`\delta (\stackrel{~}{\psi }^{}\stackrel{~}{\psi }^{}\stackrel{~}{\psi })`$ have not been included. Note also, that in the cases $`i=q`$ or $`j=q`$, for which the differential equations (45) and (47) take the same form, further terms in $`\widehat{K}`$ have to be taken into account. These terms are of second order in $`\widehat{\varphi }`$, $`\widehat{\varphi }^{}`$ and have been neglected in (40). We will in the following assume these special cases to be included. Equations (45) and (47) are then the analogues of eqns. (26) and (27) of .

We can integrate (45) and (47) formally and insert the solutions back into (29).

$`\dot{\widehat{\beta }_q}(t)`$ $`=`$ $`i\omega _q\widehat{\beta }_q(t)`$ (52)
$`{\displaystyle \frac{i}{\mathrm{}}}g{\displaystyle \underset{ij}{}}\left\{2A_{q,ij}^{}e^{i(\omega _j\omega _i)t}\widehat{f}_{ij}(0)+B_{q,ij}^{}e^{i(\omega _j+\omega _i)t}\widehat{g}_{ij}(0)+\stackrel{~}{B}_{q,ij}^{}e^{i(\omega _j+\omega _i)t}\widehat{g}_{ij}^{}(0)\right\}`$
$`2(g/\mathrm{})^2{\displaystyle \underset{ij}{}}{\displaystyle }dt\{2(n_i^0n_j^0)|A_{q,ij}|^2e^{i(\omega _i\omega _j)(tt^{})}`$
$`+(1+n_i^0+n_j^0)[|B_{q,ij}|^2e^{i(\omega _i+\omega _j)(tt^{})}|\stackrel{~}{B}_{q,ij}|^2e^{i(\omega _i+\omega _j)(tt^{})}]\}.`$
The terms involving the time integrals then lead, in the Markov approximation, to an energy shift $`\delta \omega _q`$ and a decay width $`\gamma _q`$. The energy shift may be calculated in the Popov approximation ($`\stackrel{~}{m}^0=0`$) , where it vanishes for $`q0`$, i.e. leaves the spectrum gapless as required by the Hugenholtz-Pines theorem . For temperatures $`k_BT\mu `$, however, the anomalous density $`\stackrel{~}{m}^0`$ can not be neglected in the calculations. A self consistent solution of the Hartree-Fock-Bogoliubov equations including $`\stackrel{~}{m}^0`$, together with a renormalization of the coupling $`g`$ results in a gapless excitation spectrum .

We are here primarily interested in the damping of the excited quasiparticle mode. The decay width $`\gamma _q`$ follows, analogously to  to be

$$\gamma _q=\gamma _L(q)+\gamma _B(q),$$
(53)
$`\gamma _L`$ and $`\gamma _B`$ are known as the Landau and Beliaev decay widths  and are found from (52) to be:

$`\gamma _L(q)`$ $`=`$ $`8\pi (g/\mathrm{})^2{\displaystyle \underset{ij}{}}|A_{q,ij}|^2(n_i^0n_j^0)\delta (\omega _q+\omega _i\omega _j),`$ (54)
$`\gamma _B(q)`$ $`=`$ $`4\pi (g/\mathrm{})^2{\displaystyle \underset{ij}{}}|B_{q,ij}|^2(1+n_i^0+n_j^0)\delta (\omega _q\omega _i\omega _j).`$ (55)
To derive (54,55) we have chosen the Markov approximation. This means that we assume that the sum over $`i`$ and $`j`$ over the appearing products of coupling functions (30), multiplied by the energy exponentials, give functions of $`\tau =tt^{}`$ which are sharply peaked around $`\tau =0`$. The conditions for this assumption to hold are generally that the energies of the bath modes are closely spaced compared to the energy of the decaying mode and that the couplings are sufficiently well behaved at large energies. The time dependence of the sharply peaked function may then be approximated as Markovian, i.e. by a delta-distribution, and integration over time yields the energy conservation delta-distributions in (54,55).

The first term on the r.h.s. of (52) is the Langevin noise operator corresponding to the decay of the quasiparticle mode. Hence, the equation of motion for $`\widehat{\beta }_q`$ may be written in the form (3), where the decay width $`\gamma _q`$ is given by (53)–(55) and the Langevin noise operator is

$`\widehat{F}_q(t)`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}g{\displaystyle \underset{ij}{}}\{2A_{q,ij}^{}e^{i(\omega _j\omega _i)t}\widehat{f}_{ij}(0)`$ (58)
$`+B_{q,ij}^{}e^{i(\omega _j+\omega _i)t}\widehat{g}_{ij}(0)`$
$`+\stackrel{~}{B}_{q,ij}^{}e^{i(\omega _j+\omega _i)t}\widehat{g}_{ij}^{}(0)\}`$
Using this expression we will now show that in the Markov approximation $`\widehat{F}_q`$ obeys the commutation relation (4) and its two-time correlation function is given by (6).

To this end we need the commutation relations

$`[\widehat{f}_{ij}(0),\widehat{f}_{kl}^{}(0)]`$ $`=`$ $`\delta _{ik}\delta _{jl}(n_i^0n_j^0),`$ (59)
$`[\widehat{g}_{ij}(0),\widehat{g}_{kl}^{}(0)]`$ $`=`$ $`\left(\delta _{il}\delta _{jk}+\delta _{ik}\delta _{jl}\right)(1+n_i^0+n_j^0).`$ (60)
In deriving these we have neglected on the right hand sides terms proportional to the fluctuation operators $`\widehat{f}_{ij}`$. Using (59,60) and the fact, that mixed commutators of $`\widehat{f}_{ij}`$ and $`\widehat{g}_{kl}`$ only yield fluctation operators $`\widehat{g}_{ij}`$ we arrive at

$`[\widehat{F}_q(t)e^{i\omega _qt},\widehat{F}_q^{}(t^{})e^{i\omega _qt^{}}]`$ $`=`$ $`2(g/\mathrm{})^2{\displaystyle \underset{ij}{}}\{2(n_i^0n_j^0)|A_{q,ij}|^2e^{i(\omega _q+\omega _i\omega _j)(tt^{})}`$ (63)
$`+(1+n_i^0+n_j^0)[|B_{q,ij}|^2e^{i(\omega _q\omega _i\omega _j)(tt^{})}`$
$`|\stackrel{~}{B}_{q,ij}|^2e^{i(\omega _q+\omega _i+\omega _j)(tt^{})}]\}.`$
In the Markov approximation the sum over $`i,j`$ is replaced by a function proportional to $`\delta (tt^{})`$. Integration of (63) and comparison with (53)–(55) then shows that the coefficient of this delta-distribution is given by $`\gamma _q`$, as stated in (4).

Eqn. (6) may be proven in a similar way, using

$`\widehat{f}_{ij}^{}(0)\widehat{f}_{kl}(0)`$ $`=`$ $`\delta _{ik}\delta _{jl}n_j^0(n_i^0+1),`$ (64)
$`\widehat{g}_{ij}^{}(0)\widehat{g}_{kl}(0)`$ $`=`$ $`\left(\delta _{il}\delta _{jk}+\delta _{ik}\delta _{jl}\right)n_i^0n_j^0.`$ (65)
In deriving these relations we have made use of Wick’s theorem. Taking into account still the respective energy matchings in the different terms, the population factors become

$`\delta (\omega _q+\omega _i\omega _j)n_j^0(n_i^0+1)`$ $`=`$ $`\delta (\omega _q+\omega _i\omega _j)n_q^0(n_i^0n_j^0),`$ (66)
$`\delta (\omega _q\omega _i\omega _j)n_i^0n_j^0`$ $`=`$ $`\delta (\omega _q\omega _i\omega _j)n_q^0(1+n_i^0+n_j^0).`$ (67)
where we used the Planck distribution $`n_i^0=\widehat{\beta }_i^{}\widehat{\beta }_i_0=[\mathrm{exp}(\mathrm{}\omega _i/k_BT)1]^1`$. Together with (60) we arrive at

$`\widehat{F}_q^{}(t)\widehat{F}_q(t^{})e^{i\omega _q(tt^{})}`$ $`=`$ $`n_q^0[\widehat{F}_q(t),\widehat{F}_q^{}(t^{})]e^{i\omega _q(tt^{})}.`$ (68)
where the r.h.s. is given by (63). This yields, in the Markov approximation, the identity (6) for the correlation function.

## III Heisenberg-Langevin equation for systems of atoms in two internal states

In the following we would like to extend the procedure reviewed in sec. II to the case with two internal atomic degrees of freedom, where the atoms in the excited mode $`q`$ are in a different internal state than the condensate atoms. Such situations are important e.g. in the context of the dynamics of dark states involving Bose-Einstein condensates which are the basis for experiments using electromagnetically induced transparency  for the slowing and stopping of light pulses in atomic matter . The decay of the dark state induced by atomic collisions limits the maximum time for which a light pulse may be stored in collective atomic excitations of condensates .

The interactions which cause the quasiparticle decay will be mainly collisions between atoms in the excited mode and condensate atoms, such that they do not change the respective internal states. Let us call the two internal levels of the atoms $`b`$ and $`c`$, the Bose-condensed atoms being in internal state $`|b`$. The Hamiltonian in the grand canonical ensemble then has the form

$`\widehat{K}`$ $``$ $`\widehat{H}\mu \widehat{N}=\widehat{K}_b+\widehat{K}_c+\widehat{K}_{bc},`$ (69)
$`\widehat{K}_b`$ $`=`$ $`{\displaystyle d^3r\left\{\widehat{\psi }_b^{}\widehat{K}_{0,b}\widehat{\psi }_b+\frac{g_{bb}}{2}\widehat{\psi }_b^{}\widehat{\psi }_b^{}\widehat{\psi }_b\widehat{\psi }_b\right\}},`$ (70)
$`\widehat{K}_c`$ $`=`$ $`{\displaystyle d^3r\left\{\widehat{\psi }_c^{}\left(\widehat{K}_{0,c}+\omega _{cb}\right)\widehat{\psi }_c+\frac{g_{cc}}{2}\widehat{\psi }_c^{}\widehat{\psi }_c^{}\widehat{\psi }_c\widehat{\psi }_c\right\}},`$ (71)
$`\widehat{K}_{bc}`$ $`=`$ $`{\displaystyle d^3rg_{bc}\widehat{\psi }_c^{}\widehat{\psi }_b^{}\widehat{\psi }_c\widehat{\psi }_b},`$ (72)
where $`\widehat{K}_{0,\nu }=\mathrm{}^2\mathbf{}^2/(2m)+V_\nu (𝐫)\mu `$, $`\nu =b,c`$, and $`\omega _{cb}=\omega _c\omega _b`$ is the frequency difference of the internal levels. $`g_{\mu \nu }`$, $`\mu ,\nu =b,c`$ are the collisonal coupling constants for the respective processes being related to the respective scattering lengths by $`g_{\mu \nu }=4\pi \mathrm{}^2a_{\mu \nu }/m`$. As for the single level case we split the field operators for the atoms in level $`b`$ into mean-field and fluctuation parts,

$$\widehat{\psi }_b(𝐫,t)=\mathrm{\Phi }_0(𝐫)+\stackrel{~}{\psi }_b(𝐫,t),$$
(73)
and the operator for the atoms in $`c`$ into fluctuation and excited modes:

$$\widehat{\psi }_c(𝐫,t)=\stackrel{~}{\psi }_c(𝐫,t)+\widehat{\varphi }_c(𝐫,t),$$
(74)
With these, we find the equation of motion for the excitation operator $`\widehat{\varphi }_c`$ in an analogous way as in the single level case:

$`i\mathrm{}{\displaystyle \frac{}{t}}\widehat{\varphi }_c`$ $`=`$ $`\left(\widehat{K}_{0,c}+2g_{cc}\stackrel{~}{n}_c^0+g_{bc}\left[n_0+\stackrel{~}{n}_b^0\right]\right)\widehat{\varphi }_c+g_{cc}\stackrel{~}{m}_c^0\widehat{\varphi }_c^{}`$ (76)
$`+g_{bc}\mathrm{\Phi }_0\left[\delta \widehat{\stackrel{~}{n}}_{bc}+\delta \widehat{\stackrel{~}{m}}_{bc}\right],`$
where $`\stackrel{~}{n}_\nu ^0=\stackrel{~}{\psi }_\nu ^{}\stackrel{~}{\psi }_\nu _0`$, $`\stackrel{~}{m}_\nu ^0=\stackrel{~}{\psi }_\nu \stackrel{~}{\psi }_\nu _0`$ are the mean equilibrium density and anomalous density of the non-condensed fraction in $`\nu =b,c`$. The fluctuation operators $`\delta \widehat{\stackrel{~}{n}}_{bc}=\stackrel{~}{\psi }_b^{}\stackrel{~}{\psi }_c`$, $`\delta \widehat{\stackrel{~}{m}}_{bc}=\stackrel{~}{\psi }_b\stackrel{~}{\psi }_c`$ are defined as in (14,15), with their equilibrium mean values though being zero.

In order to derive from (76) the equation of motion for the Fock operator of the excited mode we use the decomposition (21) for $`\stackrel{~}{\psi }_b`$ and a particle mode decomposition for atoms in the internal state $`c`$:

$`\stackrel{~}{\psi }_c(𝐫,t)`$ $`=`$ $`{\displaystyle \underset{i}{}}\stackrel{~}{\psi }_{c,i}(𝐫)\widehat{c}_i(t).`$ (77)
Here the mode functions $`\stackrel{~}{\psi }_{c,i}`$ are solutions of the equation

$$^c(𝐫)\stackrel{~}{\psi }_{c,i}(𝐫)=\mathrm{}\omega _i^c\stackrel{~}{\psi }_{c,i}(𝐫),$$
(78)
where $`^c(𝐫)=\widehat{K}_{0,c}(𝐫)+\omega _{cb}+2g_{cc}\stackrel{~}{n}_c^0+g_{bc}\left[n_0+\stackrel{~}{n}_b^0\right]`$, and where we have chosen the approximation $`^c(𝐫)=g_{cc}\stackrel{~}{m}_c^0=0`$. Note that the eigenfunctions then are not quasiparticle modes since the zero mode in $`c`$ is not macroscopically populated. The equation of motion for $`\widehat{c}_q`$ then reads

$`\dot{\widehat{c}}_q(t)`$ $`=`$ $`i\omega _q^c\widehat{c}_q(t){\displaystyle \frac{i}{\mathrm{}}}g_{bc}{\displaystyle \underset{ij}{}}\left\{A_{q,ij}^{bc}\widehat{f}_{ij}^{bc}(t)+B_{q,ij}^{bc}\widehat{g}_{ij}^{bc}(t)\right\},`$ (79)
where

$`A_{q,ij}^{bc}`$ $`=`$ $`{\displaystyle d^3r\mathrm{\Phi }_0\left[\stackrel{~}{\psi }_{c,q}(u_i+v_i)\stackrel{~}{\psi }_{c,j}^{}\right]},`$ (80)
$`B_{q,ij}^{bc}`$ $`=`$ $`{\displaystyle d^3r\mathrm{\Phi }_0\left[\stackrel{~}{\psi }_{c,q}(u_i^{}+v_i^{})\stackrel{~}{\psi }_{c,j}^{}\right]},`$ (81)
and

$`\widehat{f}_{ij}^{bc}`$ $`=`$ $`\widehat{\beta }_i^{}\widehat{c}_j,`$ (82)
$`\widehat{g}_{ij}^{bc}`$ $`=`$ $`\widehat{\beta }_i\widehat{c}_j.`$ (83)
In order to solve (79) we need, as in the single level case, the time evolution of the operators $`\widehat{f}_{ij}^{bc},\widehat{g}_{ij}^{bc}`$. To this end we write down the relevant terms of the Hamiltonian $`\widehat{K}`$ which contribute to the commutator with these operators. Using (73) and (74) as well as (77) and (78) they read

$`\widehat{K}_c+\widehat{K}_{bc}`$ $`=`$ $`{\displaystyle \underset{i}{}}\mathrm{}\omega _i\widehat{c}_i^{}\widehat{c}_i+g_{cc}{\displaystyle d^3r\delta (\stackrel{~}{\psi }_c^{}\stackrel{~}{\psi }_c^{}\stackrel{~}{\psi }_c\stackrel{~}{\psi }_c)},`$ (87)
$`+g_{bc}{\displaystyle }d^3r\mathrm{\Phi }_0[\stackrel{~}{\psi }_b^{}\widehat{\varphi }_c^{}\stackrel{~}{\psi }_c+\widehat{\varphi }_c^{}\stackrel{~}{\psi }_c\stackrel{~}{\psi }_b+\mathrm{h}.\mathrm{c}.]`$
$`+g_{bc}{\displaystyle }d^3r[\delta \stackrel{~}{n}_b\stackrel{~}{\psi }_c^{}\stackrel{~}{\psi }_c+\delta \stackrel{~}{n}_c\stackrel{~}{\psi }_b^{}\stackrel{~}{\psi }_b+(\delta \stackrel{~}{n}_{bc}\stackrel{~}{\psi }_c^{}\stackrel{~}{\psi }_c+\delta \stackrel{~}{m}_{bc}^{}\stackrel{~}{\psi }_c\stackrel{~}{\psi }_b+\mathrm{h}.\mathrm{c}.)]`$
$`+g_{bc}{\displaystyle }d^3r[\delta (\stackrel{~}{\psi }_c^{}\stackrel{~}{\psi }_b^{}\stackrel{~}{\psi }_c\stackrel{~}{\psi }_b)+\mathrm{\Phi }_0(\delta (\stackrel{~}{\psi }_c^{}\stackrel{~}{\psi }_c\stackrel{~}{\psi }_b)+\mathrm{h}.\mathrm{c}.)],`$
where

$`\delta (\stackrel{~}{\psi }_c^{}\stackrel{~}{\psi }_b^{}\stackrel{~}{\psi }_c\stackrel{~}{\psi }_b)`$ $`=`$ $`\stackrel{~}{\psi }_c^{}\stackrel{~}{\psi }_b^{}\stackrel{~}{\psi }_c\stackrel{~}{\psi }_b\stackrel{~}{\psi }_c^{}\stackrel{~}{\psi }_c\stackrel{~}{\psi }_b^{}\stackrel{~}{\psi }_b\stackrel{~}{\psi }_b^{}\stackrel{~}{\psi }_b\stackrel{~}{\psi }_c^{}\stackrel{~}{\psi }_c(\stackrel{~}{\psi }_c^{}\stackrel{~}{\psi }_b\stackrel{~}{\psi }_b^{}\stackrel{~}{\psi }_c\stackrel{~}{\psi }_c^{}\stackrel{~}{\psi }_b^{}\stackrel{~}{\psi }_c\stackrel{~}{\psi }_b+\mathrm{h}.\mathrm{c}.),`$ (88)
$`\delta (\stackrel{~}{\psi }_c^{}\stackrel{~}{\psi }_c^{}\stackrel{~}{\psi }_b)`$ $`=`$ $`\stackrel{~}{\psi }_c^{}\stackrel{~}{\psi }_c^{}\stackrel{~}{\psi }_b\stackrel{~}{n}_c^0\stackrel{~}{\psi }_b,`$ (89)
$`\delta \stackrel{~}{n}_c`$ $`=`$ $`\stackrel{~}{\psi }_c^{}\stackrel{~}{\psi }_c,`$ (90)
and $`\delta \stackrel{~}{n}_{bc}=\delta \widehat{\stackrel{~}{n}}_{bc}`$, $`\delta \stackrel{~}{m}_{bc}=\delta \widehat{\stackrel{~}{m}}_{bc}`$. The operator $`\delta (\stackrel{~}{\psi }_c^{}\stackrel{~}{\psi }_c^{}\stackrel{~}{\psi }_c\stackrel{~}{\psi }_c)`$ and the functions $`\delta \stackrel{~}{n}_b`$ are defined as in (41) and (43) respectively. From the terms on the right hand side of (87) only the quadratic one and the ones in the second line will contribute to the leading order time evolution of $`\widehat{f}_{ij}^{bc}`$ and $`\widehat{g}_{ij}^{bc}`$. The equations of motion for these operators are then found to be

$`{\displaystyle \frac{}{t}}\widehat{f}_{ij}^{bc}`$ $`=`$ $`i(\omega _j^c\omega _i^b)\widehat{f}_{ij}^{bc}{\displaystyle \frac{i}{\mathrm{}}}g_{bc}(n_i^bn_j^c)A_{q,ij}^{bc}\widehat{c}_q(t),`$ (91)
$`{\displaystyle \frac{}{t}}\widehat{g}_{ij}^{bc}`$ $`=`$ $`i(\omega _j^c+\omega _i^b)\widehat{g}_{ij}^{bc}{\displaystyle \frac{i}{\mathrm{}}}g_{bc}(n_i^bn_j^c)B_{q,ij}^{bc}\widehat{c}_q(t),`$ (92)
where we have neglected terms of higher order in the fluctuation opertors and used $`\widehat{\varphi }_q=\stackrel{~}{\psi }_{c,q}\widehat{c}_q`$. $`n_i^b`$ and $`n_j^c`$ are the equilibrium populations of modes $`i`$ and $`j`$ in internal states $`b`$ and $`c`$, respectively. Analogously to the single level case we then obtain, in the Markov approximation, the equation of motion for $`\widehat{c}_q`$ in the form (3). We find that the decay rate $`\gamma _q=\gamma _L(q)+\gamma _B(q)`$ has Landau- and Beliaev-like contributions (53) which read

$`\gamma _L(q)`$ $`=`$ $`2\pi (g_{bc}/\mathrm{})^2{\displaystyle \underset{ij}{}}|A_{q,ij}^{bc}|^2(n_i^0n_j^0)\delta (\omega _q+\omega _i\omega _j),`$ (93)
$`\gamma _B(q)`$ $`=`$ $`2\pi (g_{bc}/\mathrm{})^2{\displaystyle \underset{ij}{}}|B_{q,ij}^{bc}|^2(1+n_i^0+n_j^0)\delta (\omega _q\omega _i\omega _j),`$ (94)
The Langevin noise operator is given by

$`\widehat{F}_q(t)`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}g_{bc}{\displaystyle \underset{ij}{}}\{A_{q,ij}^{bc}e^{i(\omega _j^c\omega _i^b)t}\widehat{f}_{ij}^{bc}(0)`$ (96)
$`+B_{q,ij}^{bc}e^{i(\omega _j^c+\omega _i^b)t}\widehat{g}_{ij}^{bc}(0)\}.`$
The proofs that (96) obeys the commutator (4) and eqn. (6) for the two-time correlation function are conducted in analogy to the single mode case, eqns. (59)–(68).

Our result (93,94) for the damping rate $`\gamma _q=\gamma _L(q)+\gamma _B(q)`$ may e.g. be used to calculate the loss of coherence in a dark state used to delay and store light pulses. In  we have calculated the damping rate for a homogeneous condensate and found that the Beliaev rate for small momenta $`q`$ of the decaying mode is proportional to $`q^5`$ as in the case where the internal atomic state is the same as that of the condensate atoms. The rate differs only by a constant numerical factor: $`\gamma _C(q)_{T=0}\mathrm{}q^5/96m\pi n_0`$ compared to $`\gamma _C(q)_{T=0}3\mathrm{}q^5/320m\pi n_0`$ in the single level case (assuming $`a_{bc}=a_{bb})`$.

## IV Dissipative dynamics of higher correlation functions: Atom-photon squeezing

We would now like to apply our formalism in order to calculate the dissipative dynamics of expectation values of two and four quasiparticle operators. This will enable us to compute the time evolution of the relative number squeezing between atoms and photons produced in photon scattering off condensates. Such schemes have been proposed in . In extending this work in  we have considered the production of atom-photon squeezing using a finite homogeneous Bose-Einstein condensate of interacting atoms at zero temperature, illuminated by laser light with given frequency $`\mathrm{\Omega }_l`$ detuned slightly from an atomic transition, cf. fig.1.

The laser photons with momentum $`k_l`$ predominantly interact with the condensate atoms ($`k=0`$) and may transfer them in a centre of mass mode with momentum $`\mathrm{\Delta }k`$ (see fig. 1) . Such processes yield a state with relative number squeezing in a particular mode pair, i.e. the variance of the relative occupation number is suppressed as compared to the classical case. We computed the squeezing parameter

$$\xi _3:=\left[\mathrm{\Delta }(n_an_b)\right]^2/\widehat{n}_a+\widehat{n}_b,$$
(97)
which measures the relative number variance $`[\mathrm{\Delta }(n_an_b)]^2=(\widehat{n}_a\widehat{n}_b)^2\widehat{n}_a\widehat{n}_b^2`$ in relation to its classical limit given by the sum of mean occupations $`\widehat{n}_a+\widehat{n}_b=\widehat{a}^{}\widehat{a}+\widehat{b}^{}\widehat{b}`$, of the scattered photon mode and the recoiling atom mode. (The index ’$`3`$’ of $`\xi _3`$ refers to the analogy of the relative number squeezing with spin squeezing. $`\xi _3`$ measures the variance of the 3-component $`\widehat{J}_3=\frac{1}{2}(\widehat{n}_a\widehat{n}_b)`$ compared to half the total value $`J/2=\frac{1}{4}\widehat{n}_a+\widehat{n}_b`$ of an abstract angular momentum. This is defined by the Fock operators $`\widehat{a}`$, $`\widehat{b}`$, using the Schwinger representation of angular momentum. For the issue of spin squeezing and the definition (97) cf. . It will be instructive to study also the squeezing of $`\widehat{J}_1=\frac{1}{2}(\widehat{a}^{}\widehat{b}+\widehat{b}^{}\widehat{a})`$ and $`\widehat{J}_2=\frac{1}{2i}(\widehat{a}^{}\widehat{b}\widehat{b}^{}\widehat{a})`$ for which we define $`\xi _i:=[\mathrm{\Delta }J_i]^2/(J/2)`$, $`i=1,2`$.

In the classical (uncorrelated) limit $`\xi _3=1`$, whereas for maximum squeezing we have $`\xi _3=0`$. There are two reasons why such perfect squeezing is spoilt in an experiment: Collisional atomic interactions and photon (rescattering) processes different from those which produce the squeezed state. In  we focused on the limitations to the maximum squeezing due to those atomic interactions which are taken into account in the basic Hartree-Fock-Bogoliubov approximation, in the Hamiltonian quadratic in the particle operators. There we have shown that “unwanted” collisions amongst the particles present profound limitations to the achievable squeezing, with correlations between atomic modes in the initial state being the main limiting factor at small times.

Here we would like to extend our work to include interactions between quasiparticles described by the third order terms in the Hamiltonian considered in section II. The main effect at this order of the approximation is collisional loss from the excited modes which will induce decay of the population of the recoiling atomic mode. This will also affect the time evolution of the correlation functions $`(\mathrm{\Delta }n_a)^2`$, $`(\mathrm{\Delta }n_b)^2`$, and $`\widehat{n}_a\widehat{n}_b=\widehat{n}_a\widehat{n}_b\widehat{n}_a\widehat{n}_b`$, which enter the squeezing parameter (97). We will neglect energy shifts at this order.

In order to be able to compute the dynamics of the occupation numbers and the four-operator expectation values including dissipation, we will in general have to find equations of motion of the type (7). Taking into account the interaction with the laser, the equations for the different expectation values will be coupled. The task would then be to find the equilibrium expectation values, to which the observables will evolve at large times, as e.g. $`n_q^0`$ in the case of eqn. (7). As we are here interested in the zero-temperature case only, we expect these equilibrium limits to be zero in general. It should then be justified to take the simpler route and use a regression theorem, i.e. the basic equations for the Fock operators, like (2), in conjunction with Wick’s theorem, to calculate the time dependence of the expectation values. Such a procedure is readily justified in the case without decay. We will have to show, however, that such a regression theorem may also be used for the dissipative case, and give precise rules for its application.

We start with the dynamic equation for the operators involved in the process to be studied. In the homogeneous limit, the quasiparticle modes are denoted by their momenta. We recall that the fluctuation field operator (21) may be expanded in terms of particle Fock operators, $`\stackrel{~}{\psi }(r,t)=V^{1/2}_k\widehat{b}_k\mathrm{exp}(ikr)`$, and the particle operators may be expressed by quasiparticle ones through $`\widehat{b}_k=u_k\widehat{\beta }_k+v_k\widehat{\beta }_k^{}`$ (cf. (21)). The quasiparticle transformation coefficients are given by $`u_k=(1\alpha _k^2)^{1/2}`$, $`v_k^2=u_k^21`$, with $`\alpha _k=1+\overline{k}^2\overline{k}(2+\overline{k}^2)^{1/2}`$, $`\overline{k}=|k|/k_0`$, $`k_0=\sqrt{8\pi an_0}`$, $`n_0`$ being the condensate density. The Bogoliubov quasiparticle frequencies are $`\omega _k=(\mathrm{}k_0^2/2m)\overline{k}\sqrt{2+\overline{k}^2}`$.

We assume two-photon resonance, where the difference $`\mathrm{\Delta }\mathrm{\Omega }=\mathrm{\Omega }_l\mathrm{\Omega }_s`$ of the frequencies of the laser, $`\mathrm{\Omega }_l`$, and of the scattered photon, $`\mathrm{\Omega }_s`$, equals the frequency of the recoiling atom mode, $`\mathrm{\Delta }\mathrm{\Omega }=\omega _{\mathrm{\Delta }k}`$, with $`\mathrm{\Delta }k=k_lk_s`$ being the recoil momentum. If we neglect non-resonant terms, the coupled equations of motion read

$`{\displaystyle \frac{d}{dt}}\left(\begin{array}{c}\stackrel{~}{\beta }_{\mathrm{\Delta }k}(t)\\ \stackrel{~}{a}_s^{}(t)\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}\gamma _{\mathrm{\Delta }k}/2& i\stackrel{~}{\mathrm{\Omega }}\\ i\stackrel{~}{\mathrm{\Omega }}& 0\end{array}\right)\left(\begin{array}{c}\stackrel{~}{\beta }_{\mathrm{\Delta }k}(t)\\ \stackrel{~}{a}_s^{}(t)\end{array}\right)\left(\begin{array}{c}\stackrel{~}{F}_{\mathrm{\Delta }k}(t)\\ 0\end{array}\right),`$ (106)
$`\dot{\stackrel{~}{\beta }}_{\mathrm{\Delta }k}(t)`$ $`=`$ $`(\gamma _{\mathrm{\Delta }k}/2)\stackrel{~}{\beta }_{\mathrm{\Delta }k}(t)\stackrel{~}{F}_{\mathrm{\Delta }k}(t).`$ (107)
Here we have chosen an interaction picture by defining $`\stackrel{~}{\beta }_q(t)=\widehat{\beta }_q(t)\mathrm{exp}(i\omega _qt)`$, $`\stackrel{~}{a}_s(t)=\widehat{a}_s(t)\mathrm{exp}(i\mathrm{\Omega }_st)`$, $`\stackrel{~}{F}_q(t)=\widehat{F}_q(t)\mathrm{exp}(i\omega _qt)`$, and $`\stackrel{~}{\mathrm{\Omega }}=\mathrm{\Omega }(t)\mathrm{exp}(i\mathrm{\Omega }_lt)`$. We denote by $`\stackrel{~}{\mathrm{\Omega }}=2|𝐝_{ge}𝓔^{(+)}(k_1)|g_{ge}(k_2)/\mathrm{\Delta }`$ the two photon Rabi frequency, given in terms of the dipole momentum $`𝐝_{ge}`$ between the internal ground and excited states, of the positive frequency part of the electric field describing the laser mode, the coupling $`g_{ge}(k_2)`$ between the scattered photon mode and the atom, and the detuning $`\mathrm{\Delta }`$ from the excited state. Neglecting the decay width $`\gamma _{\mathrm{\Delta }k}`$ of the quasiparticle mode as well as the corresponding Langevin operator $`\widehat{F}_{\mathrm{\Delta }k}`$, we may straightforwardly calculate the time evolution of the particle occupation numbers $`\widehat{n}_{k_s}^a`$, $`\widehat{n}_{\pm \mathrm{\Delta }k}^b`$, and of the squeezing parameter $`\xi _3=[(\mathrm{\Delta }n_s^a)^2+(\mathrm{\Delta }n_{\mathrm{\Delta }k}^b)^22\widehat{n}_s^a\widehat{n}_{\mathrm{\Delta }k}^b]/\widehat{n}_s^a+\widehat{n}_{\mathrm{\Delta }k}^b`$, using their expressions in terms of quasiparticle operators given above:

$`\widehat{n}_s^a`$ $`=`$ $`\widehat{a}_s\widehat{a}_s^{}1,`$ (108)
$`\widehat{n}_{\pm \mathrm{\Delta }k}^b`$ $`=`$ $`u_{\mathrm{\Delta }k}^2\widehat{\beta }_{\pm \mathrm{\Delta }k}^{}\widehat{\beta }_{\pm \mathrm{\Delta }k}+v_{\mathrm{\Delta }k}^2\left(\widehat{\beta }_{\mathrm{\Delta }k}^{}\widehat{\beta }_{\mathrm{\Delta }k}+1\right),`$ (109)
$`(\mathrm{\Delta }n_s^a)^2`$ $`=`$ $`\widehat{n}_s^a\left(\widehat{n}_s^a+1\right),`$ (110)
$`(\mathrm{\Delta }n_{\mathrm{\Delta }k}^b)^2`$ $`=`$ $`\widehat{n}_{\mathrm{\Delta }k}^b\left(\widehat{n}_{\mathrm{\Delta }k}^b+1\right),`$ (111)
$`\widehat{n}_s^a\widehat{n}_{\mathrm{\Delta }k}^b`$ $`=`$ $`u_{\mathrm{\Delta }k}^2|\widehat{a}_s\widehat{\beta }_{\mathrm{\Delta }k}|^2.`$ (112)
To be able to include collisional loss, however, a more careful approach is required. In appendix A we show that the dynamics of the expectation values appearing on the right hand sides of (108,109,112) is given by the system of equations

$`\left(\begin{array}{c}\stackrel{~}{\beta }_{\mathrm{\Delta }k}^{}\stackrel{~}{\beta }_{\mathrm{\Delta }k}(t)n_{\mathrm{\Delta }k}^0\\ \stackrel{~}{a}_s\stackrel{~}{a}_s^{}(t)+n_{\mathrm{\Delta }k}^0\\ \stackrel{~}{a}_s\stackrel{~}{\beta }_{\mathrm{\Delta }k}(t)\mathrm{c}.\mathrm{c}.\end{array}\right)`$ $`=`$ $`\mathrm{exp}\left[\left(\begin{array}{ccc}\gamma _{\mathrm{\Delta }k}& 0& i\stackrel{~}{\mathrm{\Omega }}\\ 0& 0& i\stackrel{~}{\mathrm{\Omega }}\\ 2i\stackrel{~}{\mathrm{\Omega }}& 2i\stackrel{~}{\mathrm{\Omega }}& \gamma _{\mathrm{\Delta }k}/2\end{array}\right)(tt_0)\right]\left(\begin{array}{c}\stackrel{~}{\beta }_{\mathrm{\Delta }k}^{}\stackrel{~}{\beta }_{\mathrm{\Delta }k}(t_0)n_{\mathrm{\Delta }k}^0\\ \stackrel{~}{a}_s\stackrel{~}{a}_s^{}(t_0)+n_{\mathrm{\Delta }k}^0\\ \stackrel{~}{a}_s\stackrel{~}{\beta }_{\mathrm{\Delta }k}(t_0)\mathrm{c}.\mathrm{c}.\end{array}\right),`$ (122)
$`\stackrel{~}{\beta }_{\mathrm{\Delta }k}^{}\stackrel{~}{\beta }_{\mathrm{\Delta }k}(t)n_{\mathrm{\Delta }k}^0`$ $`=`$ $`\mathrm{exp}[\gamma _{\mathrm{\Delta }k}(tt_0)]\left(\stackrel{~}{\beta }_{\mathrm{\Delta }k}^{}\stackrel{~}{\beta }_{\mathrm{\Delta }k}(t_0)n_{\mathrm{\Delta }k}^0\right).`$ (123)
Note that the operator ordering in the expectation values of the density operators is crucial. This ordering may of course be changed in the final result using the commutation relations. Thus (122) shows that applying the basic equations of motion (106) for the photon operators, at temperature $`T=0`$, where the equilibrium quasiparticle population $`n_q^0`$ vanishes, is only possible for their anti-normal-ordered products. More generally, products of the operators $`\stackrel{~}{a}_s`$, $`\stackrel{~}{\beta }_{\mathrm{\Delta }k}`$, and their respective hermitian conjugates have to be normal ordered in the sense, that $`\stackrel{~}{a}_s`$ has to be interpreted as a creation operator and $`\stackrel{~}{a}_s^{}`$ as an annihilator, whereas the quasiparticle operators are taken in the usual sense. This is due to the fact that our basic Hamiltonian in the approximation which leads to the equations (106) is not excitation number conserving due to its pair production term $`\stackrel{~}{a}_s^{}\stackrel{~}{\beta }_{\mathrm{\Delta }k}^{}`$. For this reason we have already written in eqns. (108), (109), and (112) the crucial operator products in the required ordering.

The next step would be to derive corresponding equations for products of more operators, which can be done straightforwardly along the same lines. However, as we are for the first interested in the zero-temperature case, we can restrict the effort by following the operator ordering rule stated above and using the basic time evolution equations (106), (107). In appendix A we show that this rule, applied to products of an even number of operators, in the Markov approximation leads to the correct time evolution of the expectation values at $`T=0`$.

In the rest of this paper we present the results for the time evolution of the occupation numbers and squeezing in the same system as considered in , but including quasiparticle damping. We consider specifically a homogeneous condensate of $`N=10^6`$ <sup>23</sup>Na atoms within a volume $`V=10^8`$cm<sup>3</sup>, which we choose such that the condensate density $`n_0=10^{14}`$cm<sup>-3</sup> is of the order which is typically observed in experiments. With $`a=2.8`$nm we thus have $`k_0=2.6510^6`$m<sup>-1</sup>, $`\mathrm{}k_0^2/2m=(2\pi )\mathrm{\hspace{0.17em}1.5}`$kHz. The time dependence of the occupation numbers and the squeezing parameter is shown in Figs. 2 and 3. These results have been discussed in detail in . We plot the logarithm of the occupation numbers $`\widehat{n}_{k_2}^a`$, $`\widehat{n}_{\mathrm{\Delta }k}^b`$, $`\widehat{n}_{\mathrm{\Delta }k}^b`$, for $`\overline{\mathrm{\Delta }k}=\mathrm{\Delta }k/k_0=5`$, i.e. for coupling into the particle regime of the Bogoliubov spectrum, where $`\omega _k\mathrm{}k^2/2m`$. We have chosen a two photon Rabi-frequency of $`\stackrel{~}{\mathrm{\Omega }}=1`$s<sup>-1</sup>. Initially the atomic occupation is $`\widehat{n}_{\mathrm{\Delta }k}^b(t=0)=v_{\mathrm{\Delta }k}^2`$ due to collisional interactions. Only for $`t\stackrel{>}{}0.1`$ms the photon and atom number occupations are of the same magnitude, $`\widehat{n}_{k_2}^a\widehat{n}_{\mathrm{\Delta }k}^b`$. For $`\widehat{n}_{\pm \mathrm{\Delta }k}^b`$ approaching $`N_0=10^6`$, the calculation yields non-physical results since we have neglected the depletion of $`N_0`$. The range of times where the calculation is valid may however be increased by choosing a smaller $`\stackrel{~}{\mathrm{\Omega }}`$.

The effects of the quasiparticle damping due to the finite width $`\gamma _{\mathrm{\Delta }k}`$ become visible for $`t\stackrel{>}{}\mathrm{\hspace{0.17em}1}`$ms. The curves shown correspond to a width of $`\gamma _{5k_0}=2.810^3\omega _{5k_0}=7.010^2`$s<sup>-1</sup>. The damping affects all occupation numbers shown, also that of the photon mode, since the photon production rate depends on the occupation of the recoiling atomic mode. A hypothetical stronger decay rate would result in a time evolution, where $`\widehat{n}_{\pm \mathrm{\Delta }k}^b`$ approach the initial value $`v_{\mathrm{\Delta }k}^2`$ for large times and $`\widehat{n}_{k_2}^a`$ grows at a considerably lower rate.

Fig. 3 shows the squeezing parameter $`\xi _3`$. For times $`t\stackrel{<}{}10\mu `$s it is close to $`1`$, i.e. there is essentially no squeezing compared to the classical case of independent modes. This is because the relative atom-photon number variance is still dominated by the initial “thermal” variance of the atom mode population, $`[\mathrm{\Delta }n_{\mathrm{\Delta }k}^b(t=0)]^2=\widehat{n}_{\mathrm{\Delta }k}^b(\widehat{n}_{\mathrm{\Delta }k}^b+1)=v_{\mathrm{\Delta }k}^2(v_{\mathrm{\Delta }k}^2+1)`$, which, in the particle regime where $`v_{\mathrm{\Delta }k}^21`$, resembles the variance of a coherent state. (For recoil momenta in the quasiparticle regime, $`\mathrm{\Delta }k<k_0`$ the initial value of $`\xi _3`$ would be considerably larger than 1.)

For $`t\stackrel{>}{}10\mu `$s the squeezing starts to improve due to the atom-photon pair production, but only until $`t0.1`$ms. It becomes clear that the quasiparticle damping causes the squeezing to deteriorate already at a considerably earlier time and larger value of $`\xi _3`$ compared to the loss-less case. Decoherence due to quasiparticle damping is therefore a further important limiting factor.

In the context of spin squeezing, we found analytically that within the restricted Fock space of the modes $`\pm \mathrm{\Delta }k`$, $`k_2`$, we have $`\widehat{J}_1=\widehat{J}_20`$ for all times. Thus the mean abstract spin vector $`\widehat{𝐉}`$ points in 3-direction. Fig. 3 shows, that although the squeezing of all three spin components becomes worse for large times compared to the case of coherent modes, there is a relative squeezing of up to 2 orders of magnitude of the 3-component compared to the 1- and 2-components which are of equal size

## V Conclusions

In this paper we have studied the dynamics of interacting condensed Bose-gases at finite temperature. We have set up a system of equations of motion for the operators and operator products which describe the fluctuations around the mean field and used the linear response approximation to describe the time evolution of the system. We have focused on the dissipation in the time evolution of systems which initially are out of equilibrium, specifically for initial states with a single excited (quasiparticle) eigenmode. As is well known, the decay of this mode is caused predominantly by interaction processes of the quasiparticles which involve a condensate atom and which can be classified as Landau or Beliaev processes. We have reviewed the derivation of the general expressions of the quasiparticle decay rates for a trapped condensate at finite temperature and extended this to the case where the atoms in the initially excited mode are in a different internal state than the condensate atoms. The interactions taken into account in this two-level case are collisions between atoms in different internal states which conserve the internal quantum number. We calculated the coupling functions which enter the expressions of the Beliaev and Landau decay rates and found the same infrared $`q^5`$-proportionality of the Beliaev rate as for the single internal level case.

In order to gain a consistent description of the dissipative dynamics of densities and higher correlation functions in the long time limit we used the Heisenberg-Langevin approach. In this way the time evolution equations of correlation functions including their correct equilibrium limit may be determined straightforwardly. In this paper we are interested in the zero temperature dynamics of expectation values, which yield the number squeezing between different Fock modes of a system. We were able to introduce and justify a regression theorem which we used to compute higher correlation functions using the basic equations of motion for the Fock operators only.

Specifically, we considered Bragg-scattering of photons off a finite homogeneous Bose-Einstein condensate of sodium atoms and studied the relative occupation number squeezing between the scattered photon and the recoiling atomic mode. The parameter which quantifies the squeezing measures the variance of the number difference of the modes compared to its classical limit which is given by the sum of the mean occupations. Hence our calculations required to determine the evolution of densities as well as mean values of products of four operators. Our numerical results show that there is a substantial limitation of squeezing due to the collisional interactions between the atoms. Compared to the approximation where quasiparticle interactions have been neglected the achievable squeezing is, in addition, considerably reduced by decoherence due to dissipation.

## Acknowledgments

I am especially grateful to Prof. Keith Burnett who has helped much to make my stay in Oxford possible and joyful. I would like to thank him, as well as Dr. Thorsten Köhler, Dr. Samuel Morgan, and David C. Roberts for many lively discussions and valuable suggestions. This work has been supported by the Alexander von Humboldt-Foundation, and the European Community under contract no. HPMF-CT-1999-0023.

## A Operator ordering rule

In this appendix we derive the time evolution equations (122), (123), for the expectation values appearing on the right hand sides of (108,109). Moreover we justify our procedure, that, for $`T=0`$ and in the Markov approximation we may compute expectation values of normal ordered products of more than two operators by using the basic equations (106,107), without taking into account the Langevin operators. Normal ordering here means, that $`\stackrel{~}{a}_s`$ is treated as a creation operator, and $`\stackrel{~}{a}_s^{}`$ as an annihilator, whereas the quasiparticle operators are treated in the usual way.

In order to obtain the dynamical equations for the expectation values of the required products of two operators, we need to solve the following equations which are derived from (106,107).

$`{\displaystyle \frac{d}{dt}}\left(\begin{array}{c}\stackrel{~}{\beta }_{\mathrm{\Delta }k}^{}\stackrel{~}{\beta }_{\mathrm{\Delta }k}\\ \stackrel{~}{a}_s\stackrel{~}{a}_s^{}\\ \stackrel{~}{a}_s\stackrel{~}{\beta }_{\mathrm{\Delta }k}\mathrm{c}.\mathrm{c}.\end{array}\right)`$ $`=`$ $`\left(\begin{array}{ccc}\gamma _{\mathrm{\Delta }k}& 0& i\stackrel{~}{\mathrm{\Omega }}\\ 0& 0& i\stackrel{~}{\mathrm{\Omega }}\\ 2i\stackrel{~}{\mathrm{\Omega }}& 2i\stackrel{~}{\mathrm{\Omega }}& \gamma _{\mathrm{\Delta }k}/2\end{array}\right)\left(\begin{array}{c}\stackrel{~}{\beta }_{\mathrm{\Delta }k}^{}\stackrel{~}{\beta }_{\mathrm{\Delta }k}\\ \stackrel{~}{a}_s\stackrel{~}{a}_s^{}\\ \stackrel{~}{a}_s\stackrel{~}{\beta }_{\mathrm{\Delta }k}\mathrm{c}.\mathrm{c}.\end{array}\right)\left(\begin{array}{c}\stackrel{~}{F}_{\mathrm{\Delta }k}^{}\stackrel{~}{\beta }_{\mathrm{\Delta }k}+\mathrm{c}.\mathrm{c}.\\ 0\\ \stackrel{~}{F}_{\mathrm{\Delta }k}^{}\stackrel{~}{a}_s^{}\mathrm{c}.\mathrm{c}.\end{array}\right),`$ (A13)
$`{\displaystyle \frac{d}{dt}}\stackrel{~}{\beta }_{\mathrm{\Delta }k}^{}\stackrel{~}{\beta }_{\mathrm{\Delta }k}`$ $`=`$ $`\gamma _{\mathrm{\Delta }k}\stackrel{~}{\beta }_{\mathrm{\Delta }k}^{}\stackrel{~}{\beta }_{\mathrm{\Delta }k}(\stackrel{~}{F}_{\mathrm{\Delta }k}^{}\stackrel{~}{\beta }_{\mathrm{\Delta }k}+\mathrm{c}.\mathrm{c}.).`$ (A14)
The expectation values involving the Langevin operator determine the equilibrium values of the observables. Using the formal integration of (106), together with eqn. (6) for the two time correlation function, as well as eqns. (24), (25), (58) (which show that the expectation values of the Langevin operator with the Fock operators at time $`t_0`$ vanish), we find: $`\stackrel{~}{F}_{\pm \mathrm{\Delta }k}^{}\stackrel{~}{\beta }_{\pm \mathrm{\Delta }k}=\gamma _{\mathrm{\Delta }k}n_{\mathrm{\Delta }k}^0`$, $`\stackrel{~}{F}_{\mathrm{\Delta }k}^{}\stackrel{~}{a}_s^{}=0`$. This leads, upon integration of (A13), (A14), to the final result for the time evolution, eqns. (122), (123).

The crucial condition for the application of the regression theorem is the operator ordering in the expectation values of the density operators. The general rule is that products of the operators $`\stackrel{~}{a}_s`$, $`\stackrel{~}{\beta }_{\mathrm{\Delta }k}`$, and their respective hermitian conjugates have to be normal ordered. This is meant in the sense, that $`\stackrel{~}{a}_s`$ has to be interpreted as a creation operator and $`\stackrel{~}{a}_s^{}`$ as an annihilator, whereas the quasiparticle operators are taken in the usual sense. In the following we sketch the proof that using the basic time evolution equations (106,107) and Wick’s theorem for products of an even number of operators, in the Markov approximation leads to the correct time evolution of the expectation values at $`T=0`$.

The time evolution equation of a (in the above stated sense) normal ordered product of $`n`$ Fock operators, using the chain rule and eqns. (106,107), involves terms containing the Langevin operator. It is the expectation values of these terms which we have to consider and which need to vanish for the desired procedure to be justified. Upon insertion of the formal solutions of eqns. (106,107) for the Fock operators in these terms, their value is determined by expectation values of mixed products of Langevin operators at different times and Fock operators at time $`t_0`$. The normal ordering of the initially considered operators carries through to a normal ordering of these mixed products. Moreover, the terms involve time integrals over those Langevin operators which entered through the integrated eqns. (106), (107). For the case of four operators which we need to compute squeezing, these integrals are e.g. of the form

$`{\displaystyle _{t_0}^t}𝑑t^{}𝑑t^{\prime \prime }𝑑t^{\prime \prime \prime }\stackrel{~}{F}^{}(t)\stackrel{~}{F}^{}(t^{})\stackrel{~}{F}(t^{\prime \prime })\stackrel{~}{F}(t^{\prime \prime \prime }),`$ (A15)
$`{\displaystyle _{t_0}^t}𝑑t^{}\stackrel{~}{F}^{}(t)\widehat{\alpha }^{}(t_0)\widehat{\alpha }(t_0)\stackrel{~}{F}(t^{}),`$ (A16)
where $`\widehat{\alpha }`$ is either of the operators $`\stackrel{~}{a}_s^{}`$, $`\stackrel{~}{\beta }_{\pm \mathrm{\Delta }k}`$. When inserting (58) we see that, at $`T=0`$, only those terms can be different from zero, which contain an operator $`\widehat{g}_{ij}(t_0)`$ at the very left and one $`\widehat{g}_{kl}^{}(t_0)`$ at the very right. (If there is an $`\stackrel{~}{a}_s(t_0)`$ at the very left and/or an $`\stackrel{~}{a}_s^{}(t_0)`$ at the very right, the expectation value may be broken down using Wick’s theorem, and the problem reduced to the case of $`n2`$ or less operators.)

The final step is to show that the remaining non-vanishing terms are zero in the Markov approximation. This is done as before by assuming that the coupling functions (3032) which multiply the integrals (A15) vary only slowly over the range of modes, such that the sums over the mode indices result in functions which are strongly peaked at a single time and may be written as proportional to delta-distributions in time. Integration over the time argument of one of the Langevin operators at the very right or left then immeditately yields, that the remaining terms are all proportional to $`\delta (\omega _{\mathrm{\Delta }k}+\omega _i+\omega _j)`$. Since we consider $`|\mathrm{\Delta }k|0`$ only, the argument of the distribution is positive. The expectation values in the equations of motion, which involve Langevin operators, therefore vanish in the Markov approximation at $`T=0`$, for the same reason as this is the case, with $`n_q^0=0`$, for (68).


format: text

text:
# The phase diagram of QCD with two and four flavors: results with HYP fermions Talk presented by A. Hasenfratz

## 1 INTRODUCTION

In this talk we present preliminary results of the two and four flavor QCD finite temperature phase diagram on $`N_T=4`$ temporal lattices using a hypercubic blocked (HYP), staggered fermion action. Because of the large lattice spacing these simulations are more relevant to illustrate the effect of the HYP action than to extract physical values. Our goal here is to show how the improved flavor symmetry of the HYP fermions influences the phase diagram.

In any fermionic action the thin link gauge connections can be replaced by some sort of smeared links. Smearing partially removes short range vacuum fluctuations (dislocations) improving the topological properties of the gauge field configuration and improving the chiral properties of the fermionic action. At the same time smearing can distort the short range properties of the configuration introducing new type of lattice artifacts. An optimal smearing is a transformation where short range distortions are minimal while most of the dislocations are removed. The hypercubic blocking was designed to be very compact, and its parameters are optimized non-perturbatively to remove most of the vacuum fluctuations. Interestingly the non-perturbatively optimized parameters remove the perturbative tree level flavor symmetry violating terms of the action.

The HYP blocking can be combined with any fermionic action. Here we consider staggered fermions. In  we showed in the quenched approximation that flavor symmetry with HYP staggered fermions is improved by about an order of magnitude relative to thin link staggered fermions. The improved flavor symmetry is especially important for physical quantities that are sensitive to chiral symmetry. The finite temperature phase transition restores chiral symmetry at vanishing quark mass. Theoretical predictions  for this transition are based on the full chiral and flavor symmetry of the fermionic action. Lattice simulations studying the nature of the finite temperature phase transition at small quark masses could therefore be very sensitive to chiral symmetry violations of the lattice action. Since the finite temperature phase transition occurs at large lattice spacing, $`a0.30.15`$fm on typical temporal lattices $`N_T=48`$, the use of chirally improved actions could be essential there.

Before presenting our numerical results we mention a puzzle about the flavor dependence of the finite temperature phase transition. Phenomenological instanton model calculations predict a strong flavor dependence for the chiral transition temperatures . For $`N_f=2`$ flavors the chiral transition is expected to be second order at $`T_c150`$MeV with crossover at finite quark mass values. For $`N_f=3`$ flavors the transition is predicted to be first order, $`T_c100`$MeV in the chiral limit while for $`N_f=4`$ flavors the chiral transition, if it exists at all, is at a very small temperature. At finite quark mass both the $`N_f=3`$ and 4 flavor cases show first order phase transition that persists up to some critical quark mass value and the phase transition temperature depends strongly on the quark mass. Lattice simulations observe only a weak flavor dependence, especially for $`N_f=4`$ flavors. Instanton models rely on plausible but unproven assumptions yet they predict the zero temperature chiral behavior of QCD very successfully. It is puzzling why their prediction is different from lattice calculations in the finite temperature case.

## 2 $`N_F=4`$ FLAVOR SIMULATIONS

We performed simulations on $`N_T=4`$ temporal lattices with $`N_S=8,10`$ and 16 spatial size. The details of the simulations can be found in . To see the effect of smearing first we compare thin link, one level of APE and HYP smeared actions at approximately matched quark masses. In quenched simulations we found that the three actions are matched at quark masses $`am_{\mathrm{HYP}}=0.1,`$ $`am_{\mathrm{APE1}}=0.08`$, and $`am_{\mathrm{thin}}=0.06`$. Dynamical simulations show that once the quark masses are matched, the lattice spacing at fixed gauge coupling $`\beta `$ remains approximately matched also.

As it can be seen in figure 1 the strong first order transition observed with thin link action weakens and disappears as smearing is introduced at constant physical quark mass. One might wonder if the HYP action might show a first order transition at some smaller gauge coupling. We cannot exclude this possibility without actual simulations, but even if there is a first order phase transition, it occurs at a much smaller temperature than the phase transition of the thin link action suggesting that it is not a physically relevant transition.

We did not find first order signals even at smaller quark masses close to the physical light quark mass indicating that our simulations were all above the end point of the first order transition line. To identify the crossover it is traditional to look at the susceptibility of the quark condensate. We found only very broad peaks in the susceptibility. However the peak of the susceptibility does not necessarily identify the crossover transition. It is more sensitive to nearby critical points like the end point of the first order line. If the transition line has a strong quark mass-temperature dependence, this could be very different from the actual crossover region. Because of that concern we decided to use the chiral condensate at fixed gauge coupling to identify the crossover temperature. At small quark masses we expect the condensate to depend linearly on the quark mass,

$$\frac{<\overline{\psi }\psi >}{m_q}=\frac{\mathrm{\Sigma }}{m_q}+c$$

where $`\mathrm{\Sigma }`$ and $`c`$ are constants independent of the quark mass. We distinguish three different cases. If all quark masses are in the chirally symmetric phase, $`\mathrm{\Sigma }=0`$. If all quark masses are in the chirally broken phase, $`\mathrm{\Sigma }>0`$ is the chiral condensate. If at fixed gauge coupling the smaller quark mass points are in the chirally symmetric while the larger quark mass points are in the chirally broken phase, the formers predict $`\mathrm{\Sigma }=0`$ while the latter ones predict a non-zero condensate at the transition point. $`\mathrm{\Sigma }`$ however is not this value but its linear extrapolation to zero quark mass and could take any value, even negative ones. Figure 2 illustrates such an extrapolation at $`\beta =4.8`$, 4.9 and 5.0. The $`\beta =5.0`$ points are all in the chirally symmetric phase, the other two are mixed. From the plot we can read off the transition mass at fixed gauge coupling. The negative slope at larger quark masses indicates that the condensate at the crossover is small.

The predicted crossover region is plotted in figure 3. At $`\beta =5.0`$ the crossover occurs above $`am_q=0.1`$ corresponding to a temperature around $`220`$MeV. At $`\beta =4.95,`$ $`am_t=0.04`$ the crossover temperature is approximately $`T=170`$MeV. Whether this temperature changes significantly as we move to smaller quark masses has to be investigated in the future.

In any case our data indicate that the first order end point occurs at a fairly small quark mass value (and possibly small temperatures).

## 3 $`N_f=2`$ FLAVOR SIMULATIONS

We have preliminary results of the two flavor system on $`8^3\times 4`$ lattices at several quark masses. In the simulations we approximate the square root of the fermionic determinant with a finite polynomial. Some of the details of the algorithm are discussed in . The phase diagram does not show significant deviation from the thin link action results. Figure 4 shows $`<\overline{\psi }\psi >/am_q`$ as the function of the gauge coupling for $`am_q=0.01`$ and $`am_q=0.04`$. The figure suggests a phase transition around $`\beta =5.2`$. Further data at quark mass $`am_q=0.02,0.03`$ and 0.06 support this conclusion. The phase transition for the two flavor system on $`N_T=4`$ is a crossover, it occurs around $`\beta =5.2`$ and depends only weekly on the quark mass. The susceptibility of the chiral condensate predicts the same value.