Exploring the growth of correlations in a quasi one-dimensional trapped Bose gas
EExploring the growth of correlations in a quasi one-dimensional trappedBose gas
M. Eckart, R. Walser and W.P. Schleich
Institute of Quantum Physics, Ulm University, D-89069 Ulm, Germany ∗ (Dated: November 18, 2018)Phase correlations, density fluctuations and three-body loss rates are relevant for manyexperiments in quasi one-dimensional geometries. Extended mean-field theory is used toevaluate correlation functions up to third order for a quasi one-dimensional trapped Bosegas at zero and finite temperature. At zero temperature and in the homogeneous limit, wealso study the transition from the weakly correlated Gross-Pitaevskii regime to the stronglycorrelated Tonks-Girardeau regime analytically. We compare our results with the exact Lieb-Liniger solution for the homogeneous case and find good agreement up to the cross-overregime. ∗ Electronic address: [email protected] a r X i v : . [ c ond - m a t . m e s - h a ll ] D ec Contents
I. Introduction II. Lieb-Liniger theory for bosons in one dimension III. Extended mean-field theory for bosons in one dimension
IV. Analytic solution for the stationary HFB-equations in the homogeneous system at zerotemperature r (cid:28) ˜ ξ r (cid:29) ˜ ξ V. Numerical results for trapped atoms at zero and finite temperature
VI. Conclusions and outlook Acknowledgments A. Higher transcendental functions
B. Deformation of the integration contour in the complex plane References I. INTRODUCTION
Low dimensional physics has always been in the focus of interest as reduced dimensionalitygoes hand in hand with an increase of quantum fluctuations. Thus, by reducing the volume ofaccessible phase space, the effect of quantum fluctuations of the remaining degrees of freedommust be enhanced in order to comply with the fundamental Heisenberg uncertainty principle.Currently, this fact has stimulated many fascinating experiments in the context of ultracoldgases [1, 2, 3, 4, 5, 6, 7, 8], which explore various aspects of the geometric transitions. Serendip-itously, also most exactly solvable models of field theory are one-dimensional [9] and rest onthe celebrated Bethe-ansatz invented in the 1930ies. In the context of atomic Bose gases, todaymost prominent are the spatially homogeneous models of hard- and soft-core bosons on a stringof M. Girardeau [10] as well as E. Lieb and W. Liniger [11, 12]. One of the many interestingquestions which can be explored, is the cross-over from the weakly correlated Gross-Pitaevskiiregime ( γ (cid:28) ) to the strongly correlated Tonks-Girardeau regime [13, 14] ( γ (cid:29) ). Thereby,one commonly uses the Lieb-Liniger (LL) parameter γ , cf. (5), to measure the relative strength ofkinetic to repulsive self-energy in the dilute Bose gas.Today, we also have alternative tools available: First, there are the exact few-body calcula-tions, i. e. multi-channel time-dependent Hartree-Fock (MCTHF) or configuration interaction (CI)methods, which originate from atomic, molecular and nuclear physics. While originally designedfor fermionic energy structure calculation, they are nowadays also applied to few-boson systems( ≈ ˆ a x , which is the single particledensity n x at position x and the conjugate phase quadrature correlation function g (1) x,y . The fluc-tuations about the mean density are measured with the second order density-density correlationfunction g (2) x,y n x = (cid:104) ˆ a † x ˆ a x (cid:105) , g (1) x,y = (cid:104) ˆ a † y ˆ a x (cid:105)√ n x n y , g (2) x,y = (cid:104) ˆ a † x ˆ a † y ˆ a y ˆ a x (cid:105) n x n y . (1)In here, (cid:104) . . . (cid:105) = Tr { . . . ρ } denotes an average over the state of the system described by the many-body density operator ρ . Such second order correlation functions have been measured experimen-tally [4, 5, 26, 27, 28, 29], while the third order density-density-density correlation g (3) x,y,z = (cid:104) ˆ a † x ˆ a † y ˆ a † z ˆ a z ˆ a y ˆ a x (cid:105) n x n y n z , (2)became observable only recently [3, 30] via the three-body recombination rate [31].Theoretically, much attention has been directed towards second order correlation functions[32, 33, 34, 35, 36, 37, 38, 39], while less is known about the third order correlation function. Thissituation has been rectified recently in [40, 41] where the diagonal behaviour of this correlationfunction was calculated in the framework of Lieb-Liniger theory. This is where extended mean-field theory is useful, because we can calculate arbitrary orders of the correlation function andwill present calculations of the diagonal and off-diagonal behaviour of the third order correlationfunction at zero as well as finite temperature. However, the extended mean-field (EMF) approachis restricted to values of the correlation parameter γ ≤ , because any mean-field theory is knownto fail in the strongly correlated regime.This paper is organized as follows: In section II, we briefly review the central ideas of Lieb-Liniger theory [11]. This celebrated solution of the one-dimensional homogeneous Bose gas isan ideal benchmark for the extended mean-field theory [36, 42, 43, 44], whose basic conceptsare summarized in section III. In section IV, we specialize the kinetic equations to a quasi one-dimensional homogeneous situation at zero temperature for which analytical solutions can befound and compare correlation functions with the LL-predictions. The extension to inhomoge-neous, harmonically trapped systems at finite temperatures is studied numerically in section V. II. LIEB-LINIGER THEORY FOR BOSONS IN ONE DIMENSION
Lieb-Liniger theory based on the Bethe ansatz [9] describes a one-dimensional homogeneousgas of N bosons on a ring of length L . It is one of very few exactly solvable problems in many-body physics and provides a solution for every value of the correlation parameter γ . Even ininhomogeneous trapped systems this is very useful, if we can make the local density approxi-mation. In the language of second quantization, the starting point for Lieb-Liniger theory is thefollowing Hamiltonian ˆ H = (cid:90) L d x ˆ a † x (cid:18) − (cid:126) m ∂ ∂x + g a † x ˆ a x (cid:19) ˆ a x , (3)where m denotes the mass of a boson and the creation and annihilation operators satisfy the usualbosonic commutation relation. With the help of the Hellmann-Feynman theorem [45], one canobtain the diagonal part ( x = y = 0 ) of the translation invariant second order correlation function g (2) x,y = g (2) LL , introduced in (1), as L g (2) LL n = dE dg = (cid:104) Ψ | d ˆ Hdg | Ψ (cid:105) (4)by differentiating the ground state energy E with respect to the coupling constant g . Here, | Ψ (cid:105) represents the ground state and n = N/L denotes the linear particle density. It was shown byLieb and Liniger [11] that the ground state energy only depends on the dimensionless correlationparameter γ . It is basically the ratio of the repulsive mean-field energy gn to the kinetic energy (cid:126) / md at an average distance d = 1 /n . Another length scale of the problem is the healinglength ξ , which equates the kinetic energy of a wave function at scale ξ to the mean-field energy γ = mg (cid:126) n , ξ = (cid:126) √ mng . (5)We call bosons weakly correlated for γ (cid:28) (Gross-Pitaevskii regime) and strongly correlated for γ (cid:29) (Tonks-Girardeau regime).In terms of this parameter, the ground state energy and second order correlation function E = N (cid:126) n m e ( γ ) , g (2) LL = e (cid:48) ( γ ) , (6)are given in terms of the solutions of the Lieb-Liniger equations e ( γ ) = γ λ ( γ ) (cid:90) − d x h ( x, γ ) x (7) h ( x, γ ) = 12 π + 1 π (cid:90) − d y λ ( γ ) h ( y, γ ) λ ( γ ) + ( y − x ) , λ ( γ ) = γ (cid:90) − d x h ( x, γ ) . (8) FIG. 1: The Lieb-Liniger correlation function versus the correlation parameter γ . In subplot a), we depictthe second order correlator g (2) LL (solid line), the Gross-Pitaevskii approximation g (2) LL,GP (dashed dottedline) and the Tonks-Girardeau approximation g (2) LL,T G (dotted line), while subplot b) shows the third ordercorrelator g (3) LL (solid line) and the approximations g (3) LL,GP (dashed dotted line) and g (3) LL,T G (dotted line).
In the weakly correlated Gross-Pitaevskii limit γ → , as well as in the strongly correlatedTonks-Girardeau regime γ → ∞ , one obtains for the correlation function [35] g (2) LL,GP = 1 − √ γπ , for γ (cid:28) , g (2) LL,T G = 4 π γ , for γ (cid:29) . (9)A comparison of these approximations with the exact solution is presented in figure 1. The va-lidity of these results has recently been tested experimentally over a wide range of the correlationparameter [5] and will be used to probe the extended mean-field approach presented in the nextsection. −1.0 −0.5 0.0 0.5 1.0 0.0 1.0 2.0 3.0 4.0 5.0 γ = 0.01 γ = 0.1 γ = 1 γ = 10 γ = 100 x h ( x , γ ) FIG. 2: The auxiliary Lieb-Liniger function h ( x, γ ) , as a function of the variable x for five different valuesof the correlation parameter γ . It is also possible to derive an exact result for the diagonal part of the third order correlationfunction within the framework of Lieb-Liniger theory, however the task is considerably moredifficult. The exact result is derived in [41] by introducing a new function ˜ e ( γ ) which has the form ˜ e ( γ ) = γ λ ( γ ) (cid:90) − d x h ( x, γ ) x (10)and with the help of this function one obtains g (3) LL = 3˜ e (cid:48) ( γ ) − e ( γ ) − e ( γ ) e (cid:48) ( γ )2 γ + (cid:16) γ (cid:17) e (cid:48) ( γ ) + 9 e ( γ ) − e ( γ ) γ . (11)A comparison of the exact result for the third order correlation function with the approximationsin the Gross-Pitaevskii and the Tonks-Girardeau regime [35] is presented in figure 1 g (3) LL,GP = 1 − √ γπ , for γ (cid:28) , g (3) LL,T G = 16 π γ , for γ (cid:29) . (12)The auxiliary function h ( x, γ ) of (8) is depicted for various values of the correlation parameter infigure 2. III. EXTENDED MEAN-FIELD THEORY FOR BOSONS IN ONE DIMENSIONA. Time-dependent Hartree-Fock-Bogoliubov equations
The evolution of a weakly interacting dilute gas of bosons in three dimensions can be describedby a Hamiltonian ˆ H = (cid:90) d xy ˆ a † x (cid:20) H xy + 12 V bin ( x − y )ˆ a † y ˆ a y (cid:21) ˆ a x , (13) H xy = (cid:104) x | p m + V ext ( x ) | y (cid:105) , V ext ( x ) = 12 mω x + 12 mω ⊥ ( y + z ) , (14)where H is the single-particle energy in an external potential V ext and V bin is the two-particlepotential. As we are interested in the quasi one-dimensional limit, we will consider a cigar shapedtrapping configuration (angular frequencies ω and ω ⊥ ) with a large aspect ratio β . The energy andlength scales will be set by the transverse oscillator β = ω ⊥ /ω (cid:29) , a ⊥ = (cid:112) (cid:126) /mω ⊥ , ε ⊥ = (cid:126) ω ⊥ . (15)Conceptually, it is straight forward in the extended mean-field theory to use real, finite rangebinary interaction potentials and obtain proper two-body T matrices including many-body correc-tions. However, for convenience, we will use the pseudo-potential approximation in here V bin ( x − y ) = 4 π (cid:126) a s m δ x − y , (16)where a s denotes the s-wave scattering length. In order to compactify the notation we will inter-changeably use a subscript notation also for continuous functions.Extended mean-field theory uses a reduced state description based on a set of master variables { i ∈ I | γ i } . Basically, this implies the existence of a well separated hierarchy of time, energy andlength scales [46, 47] and leads to a rapid attenuation of correlation functions. Mathematicallyspeaking, it allows for a selfconsistent expansion of the full many-body density matrix ρ in termsof a perturbation series of simple many-body density matrices σ ( i ) , which depend parametricallyon the master variables ρ = σ (0) { γ } + σ (1) { γ } + O ( V bin ) . (17)This non-perturbative series in terms of the interaction potential V bin has been introduced first byChapman and Enskog in the context of kinetic theory of gases [48]. In addition to the simple seriesexpansion, we impose a selfconsistency constraint such that the operators ˆ γ i , corresponding to thec-number master variables γ i , fulfill γ i = (cid:104) ˆ γ i (cid:105) = Tr[ˆ γ i ρ ] = Tr[ˆ γ i σ (0) { γ } ] . (18)As far as the master variables are concerned, we choose the mean-field α x , the normal fluctua-tions of the single-particle density ˜ f and the fluctuations of the anomalous two-particle correlationfunction ˜ m such that (cid:104) ˆ a x (cid:105) = α x f ( c ) x , y = α x α ∗ y , m ( c ) x , y = α x α y , (19) (cid:104) ˆ a † y ˆ a x (cid:105) = f ( c ) x , y + ˜ f x , y , (cid:104) ˆ a y ˆ a x (cid:105) = m ( c ) x , y + ˜ m x , y . (20)A Gaussian operator σ (0) { γ } is compatible with the requirements of (18,19,20). In turn, this im-plies the factorizability of multi operator products (Wick’s theorem) and also yields non-Gaussiancorrections by calculating the contribution of σ (1) { γ } .By studying the coordinate transformation properties of the fluctuations [49], one finds that theaverages ˜ f and ˜ m are components of a positive semi-definite, generalized density matrix G ≥ .Thus, the system is described by a row vector χ , containing the mean-field α as well as its complexconjugate, and by the density matrix Gχ x = α x α ∗ x , G x , y = ˜ f x , y ˜ m x , y ˜ m ∗ x , y δ x , y + ˜ f ∗ x , y . (21)It can be shown from a Cauchy-Schwartz inequality that at T = 0 , the generalized density matrix G obeys an idem-potency relation Gσ G + G = 0 , σ = − . (22)Starting with the Heisenberg equation of motion, it is straightforward to derive the equations ofmotion for χ and G . In order to obtain higher-order correlation functions within the present ap-proximation scheme [42, 50], like g (2) or g (3) of (1,2), one has to evaluate the Gaussian Tr[ . . . σ (0) { γ } ] as well as the non-Gaussian contributions Tr[ . . . σ (1) { γ } ] . However it is clear that the Gaussian con-tribution will dominate for weak correlations. Thus, we have evaluated in here only the Gaussiancontributions. But already for γ ≈ , deviations from that can be noticed.0 B. Reduction to a quasi one-dimensional, stationary configuration
In a very prolate trap, the transverse motion in the directions y and z is effectively frozen outand only amplitudes proportional to the ground state ϕ ( y, z, t ) = e − ( y + z ) / a ⊥ − iω ⊥ t √ πa ⊥ (23)need to be considered. By projecting all three-dimensional functions onto the longitudinal axis x → x , one obtains the time-dependent Hartree-Fock-Bogoliubov equations (THFB) for χ x and G x,x (cid:48) i (cid:126) ∂ t χ x = Π x χ x + O ( V bin ) , i (cid:126) ∂ t G x,x (cid:48) = Σ x G x,x (cid:48) − h. c. + O ( V bin ) , (24)with the following abbreviations for the single particle Hamiltonian and the self energies Π x = Π N Π A − Π ∗ A − Π ∗ N , Σ x = Σ N Σ A − Σ ∗ A − Σ ∗ N , (25) Π N = H x + gf ( c ) x,x + 2 g ˜ f x,x , Σ N = H x + 2 gf ( c ) x,x + 2 g ˜ f x,x , (26) Σ A = gm ( c ) x,x + g ˜ m x,x , Π A = g ˜ m x,x , H x = − (cid:126) m ∂ x + 12 β mω ⊥ x . (27)In the course of the dimensional reduction, we had to introduce an effective one-dimensionalcoupling constant g = 2 (cid:126) ω ⊥ a s , which is in agreement with previous derivations [36, 51, 52]. Inorder to obtain the stationary solution for the time-independent fields χ x and G x,x (cid:48) , we make theansatz χ ( t ) = e − i (cid:126) µtσ χ, G ( t ) = e − i (cid:126) µtσ G e i (cid:126) µtσ (28)which introduces the chemical potential µ and employs the Pauli matrix σ of (22). This ansatz im-plies that the normal fluctuations ˜ f x,x (cid:48) ( t ) become time-independent, while the anomalous fluctua-tions ˜ m x,x (cid:48) ( t ) oscillate with twice the chemical potential. The properties of the resulting stationaryHFB equation will be investigated in section IV in the case of a homogeneous gas of bosons andin section V for a harmonic trapping potential, both at zero and for finite temperatures.The fact that the eigenvalue in the resulting stationary equations is indeed the chemical potential[49] can be seen from a variation of the total energy functional E ( α, ˜ f , ˜ m ) = (cid:104) ˆ H (cid:105) given by E = (cid:90) d x d y δ ( x − y )[ α ∗ y ( H x + g | α x | ) α x + ( H y + g ˜ f x,y ) ˜ f x,y ]+ g (cid:90) d x [(2 | α x | ˜ f x,x + α x ˜ m ∗ x,x + 12 | ˜ m x,x | ) + h. c. ] + O ( V bin ) , (29)with the constraint that the number of particles N = (cid:82) d x ( f ( c ) x,x + ˜ f x,x ) .1 IV. ANALYTIC SOLUTION FOR THE STATIONARY HFB-EQUATIONS IN THEHOMOGENEOUS SYSTEM AT ZERO TEMPERATURE
For the calculations that are presented in the following sections, we use standard parametersfor Rb in natural units of length a ⊥ and energy ε ⊥ m = 1 . · − kg , a s = 5 . · − m ,ω ⊥ = 2 π · Hz , ω = 2 π · Hz ,a ⊥ = 3 . · − m , ˜ g = 2 a s /a ⊥ = 3 . · − . (30)To reach the homogeneous case, we have to decrease the influence of the external trappingpotential in (27) by weakening it to the limit β (cid:29) , where we can neglect it practically. There-fore, the equilibrium state should possess the same translation symmetry as the generator of thedynamics. Consequently, we can assume that the mean field is space independent α x = α and thedensity matrix only depends on relative differences r = x − x (cid:48) χ x = χ, G x,x (cid:48) = G x − x (cid:48) = (cid:90) ∞−∞ d k π e − ikr G k . (31)Translation invariant systems are best described in Fourier-space, which was introduced above. Wecan also choose the mean-field to be real-valued by a suitable phase rotation. This is a consequenceof the global number conservation that is built into the dynamical HFB equations [53]. As themean-field is real-valued f ( c ) = m ( c ) = α , so are the fluctuations ˜ f = ˜ f x,x and ˜ m = ˜ m x,x andwith these assumptions, the normalization constraint reads n = NL = f ( c ) + ˜ f , (32)where n is the linear particle density on a length L . Furthermore the THFB equations (24) simplifysignificantly to µ = ˜ g ( f ( c ) + 2 ˜ f + ˜ m ) , k − µσ ) G k − h. c. (33)From the equation for the chemical potential, it is clear that energy and length scales emerge.It will be beneficial to introduce such scales for coherence ( k c , ξ c , ω c ), the pairing correlations( ˜ k, ˜ ξ, ˜ ω ), and their weighted sums and differences as k c = ξ − c = √ ω c , ˜ k = ˜ ξ − = √ ˜ ω, k ± = ξ ±− = √ ω ± (34) ω c = 4˜ gf ( c ) , ˜ ω = − g ˜ m , ω ± = ω c ± ˜ ω . (35)2In particular, in the Gross-Pitaevskii regime one can assume that ω c (cid:29) ˜ ω . With these definitions,one finds that the self energy is simply a × matrix in k -space with eigenvalues ω k Σ k − µσ = 12 k + ω + ω − − ω − − k − ω + , ω k = 12 (cid:112) ( k + ω c )( k + ˜ ω ) . (36)Now, we can finally evaluate the density matrix part of the HFB equations (33). Moreover, we alsohave to consider the idem-potency relation of (22). It holds for the vacuum state at zero tempera-ture and one obtains another, now quadratic relation between normal and anomalous fluctuationsin k -space ˜ m k = − ω − k + ω + (cid:18)
12 + ˜ f k (cid:19) , ˜ f k = ˜ m k − ˜ f k . (37)The system of equations can be solved point wise in k -space and leads to two solutions. One ofwhich has to be rejected on physical grounds. Thus, we find ˜ m k = − ω − ω k , ˜ f k = k + ω + ω k − . (38)The high momentum tail of the correlation functions is responsible for the short scale behaviourin real space. In the limit k → ∞ the leading terms are ˜ m k ∼ − ω − k , ˜ f k ∼ ˜ m k . (39) A. Diagonal contributions of normal and anomalous fluctuations
In order to obtain the diagonal part of the translation invariant correlation functions ˜ m r and ˜ f r at r = 0 , we have to evaluate the inverse Fourier transform of (31) ˜ m = − ω − π (cid:90) ∞−∞ dkω k , ˜ f = (cid:90) ∞−∞ dk π ˜ f k . (40)Serendipitously, this can be done exactly in terms of elliptic integrals [54] ˜ m = − k c π (cid:18) m f ( c ) (cid:19) K (cid:18) m f ( c ) (cid:19) , (41) ˜ f = ˜ m − k c π (cid:20) E (cid:18) m f ( c ) (cid:19) − K (cid:18) m f ( c ) (cid:19)(cid:21) , (42)where K and E are the complete elliptic integral of the first kind and second kind, respectively.Basic definitions are given in A 1.3 f ( c ) ˜ f − ˜ m FIG. 3: Diagonal part of the anomalous fluctuations − ˜ m (left scale, solid line) and normal fluctuations ˜ f (right scale, dashed line) versus mean field density f ( c ) . The asymptotic approximations for ˜ m (dasheddotted line) and ˜ f (dotted line) according to (43) agree well for the considered parameter range. The scaling properties of the correlation functions are most relevant for a physical insight.Thus, we can study the Gross-Pitaevskii regime of weakly correlated bosons, where | ˜ m /f ( c ) | (cid:28) and use a series expansion for the elliptic integral (1 − x ) K (1 − x ) ≈ ln(4 / √ x ) . With thisapproximation we get ˜ m = − (cid:112) ˜ gf ( c ) π W (cid:18) π (cid:113) f ( c ) / ˜ g (cid:19) , ˜ f = − (cid:112) ˜ gf ( c ) π − ˜ m . (43)In this explicit formula, we had to introduce the Lambert- W -function, which is defined in A 2 andan excellent asymptotic expansion is given in terms of logarithms [55].The approximations for the fluctuations are compared to exact numerical calculations in figure3 and give a good agreement. We will use these approximations in the following sections toevaluate the ground state energy and correlation functions.4 B. Off-diagonal contribution of normal and anomalous fluctuations
1. Short length scale behaviour: r (cid:28) ˜ ξ A rather simple, yet surprisingly efficient insight into the short range behaviour of the off-diagonal of the fluctuations can be obtained by using an iteration scheme for their Fourier trans-forms, which has its origin in (37). Starting with ˜ f (0) k = 0 and using the recursion relation ˜ m ( i +1) k = − ω − k + ω + (cid:18)
12 + ˜ f ( i ) k (cid:19) , ˜ f ( i +1) k = ( ˜ m ( i +1) k ) − ( ˜ f ( i ) k ) (44)we get a rapid convergence towards the exact results. It is remarkable that even with the inverseFourier transforms of low orders of this iteration scheme, we get a functional behaviour for ˜ f ( r ) and ˜ m ( r ) , which is equivalent to their exact behaviour for short ranges. However in contrast tothe exact form of the Fourier transforms of the fluctuations in (38), it is possible to perform theinverse Fourier transform analytically in every order of the iteration scheme. A closer look revealsthat the dependence of the fluctuations on r has to be of the form ˜ m ( i ) ( r ) = e − k + | r | P ( i ) ( r ) , ˜ f ( i ) ( r ) = e − k + | r | Q ( i ) ( r ) (45)where P ( i ) ( r ) and Q ( i ) ( r ) are polynomials in r of order i − and i +1 − respectively. Conse-quently the length scale on which the correlations decay is given by ξ + = 1 k + = 1 (cid:112) g ( f ( c ) − ˜ m ) . (46)In the Gross-Pitaevskii regime, where f ( c ) (cid:29) ˜ m , we recover the healing length ξ ≈ / (cid:112) gf ( c ) ,already introduced at the beginning in (5).The short range behaviour of the anomalous fluctuation and the normal fluctuation is depictedin figure 4. There, we compare the th order result of the iteration scheme to the exact numericalevaluation of the inverse Fourier transform. We assumed N = 100 particles, distributed over alength of L = 90 a ⊥ . This length was chosen such that the density in the homogeneous case issimilar to the density in the center of the trapped system, which will be discussed in section V. Oneobtains a good agreement between the approximation and the exact results in the regime where r (cid:28) ˜ ξ ≈ . with ξ ≈ . . At the origin we note that the anomalous fluctuation shows thetypical cusp whereas the normal fluctuation has a smooth behaviour and consequently a vanishingfirst derivative at r = 0 .5 −3 −2 −1 −3 −2 −1 − ˜ m ( r ) ˜ f ( r ) rr a)b) FIG. 4: Off-diagonal part of the normal and anomalous fluctuations versus the distance r . In subplot a)we depict the normal fluctuation ˜ f ( r ) (solid line), its short range approximation ˜ f (4) ( r ) (dashed dottedline) and the long range approximation (dotted line) according to (49). Subplot b) shows the anomalousfluctuation − ˜ m ( r ) (solid line), the short range approximation − ˜ m (4) ( r ) thereof (dashed dotted line) andthe long range approximation (dotted line) according to (48).
2. Long length scale behaviour: r (cid:29) ˜ ξ In order to get an approximation for the fluctuations in this regime, we start with the Fouriertransform of the anomalous fluctuation in (38) and note for further consideration that the Fouriertransform of the modified Bessel function of the second kind K ( c | r | ) is given by π/ √ k + c .Thus the convolution property of the Fourier transform yields ˜ m ( r ) = − ω − π (cid:90) ∞−∞ K ( k c | r (cid:48) | ) K (˜ k | r (cid:48) + r | ) dr (cid:48) . The modified Bessel function of the second kind K ( r ) is presented in figure 5. K ( r ) divergeslogarithmically at the origin and decreases exponentially for large arguments K ( r ) ∼ − γ e + ln r + O ( r ) , r → e − r (cid:2)(cid:112) π r + O ( r − / ) (cid:3) , r → ∞ (47)where γ e ≈ . denotes Euler’s constant.6 −4 −2 K ( r ) r FIG. 5: The modified Bessel function of the second kind K ( r ) versus r . In the Gross-Pitaevskii regime where f ( c ) (cid:29) | ˜ m | , we get from r (cid:29) ˜ ξ that also r (cid:29) ξ c andtherefore the first Bessel function in the integral closely resembles a δ -function. Thus we obtainthe result ˜ m ( r ) = − k c π (cid:18) m f ( c ) (cid:19) K (cid:16) ˜ k | r | (cid:17) ≈ − k c π K (cid:16) ˜ k | r | (cid:17) , (48) ˜ f ( r ) = k c πk − [ k c K (˜ k | r | ) − ˜ k K (˜ k | r | )] . (49)An alternative derivation of this result with the help of complex integration is given in B. In figure4 the asymptotic form of the fluctuations in terms of the Bessel functions is compared to an exactnumerical simulation with the parameters that were mentioned in the previous subsection. Ineither case the semi-logarithmic plot reveals an exponential decay for large distances r (cid:29) ˜ ξ , inagreement with the asymptotic behaviour of the Bessel functions. C. Comparison to Lieb-Liniger theory
Having all the ingredients at hand to calculate correlation functions, we are now ready for aquantitative comparison with the results of Lieb-Liniger theory which provides an exact solutionfor the behaviour of the second and third order correlation function.As outlined in section III A, we have to obtain the values of the correlation functions frommultiple operator averages. If ˆ o is such a general operator then (cid:104) ˆ o (cid:105) = Tr[ˆ o ( σ (0) { γ } + σ (1) { γ } + O (˜ g ))] (50)While we have already evaluated the Gaussian and non-Gaussian averages for the multinomial op-erator averages [42], it is clear that the Gaussian contribution will dominate for weak correlations.7Therefore we will focus in here on the Gaussian contribution and disregard the non-Gaussiancontributions in the following explicit expressions of order two and one, respectively g (1) x,y = f ( c ) x,y + ˜ f x,y √ n x n y + O (˜ g ) , (51) g (2) x,y = 1 + 2 (cid:60) ( f ( c ) x,y ˜ f y,x + m ( c ) x,y ∗ ˜ m y,x ) + ˜ f x,y ˜ f y,x + ˜ m x,y ˜ m ∗ y,x n x n y + O (˜ g ) , (52) g (3) x,y = 1 + 2 n x n y (cid:104) (cid:60) ( f ( c ) x,y ˜ f y,x + m ( c ) x,y ∗ ˜ m y,x ) + ˜ f x,y ˜ f y,x + ˜ m x,y ˜ m ∗ y,x (cid:105) (53) + 1 n x n y (cid:104) f ( c ) x,x ˜ f y,y + ˜ f x,x ˜ f y,y + 2 ˜ f x,y ˜ f y,x + 2 ˜ m ∗ x,y ˜ m y,x (cid:105) + 1 n y (cid:104) (cid:60) (cid:0) m ( c ) y,y ˜ m ∗ y,y (cid:1) + ˜ m y,y ˜ m ∗ y,y + f ( c ) y,y ˜ f y,y (cid:105) + 4 n x n y (cid:60) (cid:104) ˜ f x,y ( ˜ m y,y ˜ m ∗ y,x + m ( c ) y,y ˜ m ∗ y,x ) + f ( c ) x,y ( ˜ m y,y ˜ m ∗ y,x + ˜ f y,y ˜ f y,x )+ m ( c ) x,y ( ˜ f y,y ˜ m ∗ y,x + ˜ m ∗ y,y ˜ f y,x ) (cid:105) + O (˜ g ) , where n x = f ( c ) x,x + ˜ f x,x denotes the total density. This way we can calculate the full diagonaland off-diagonal behaviour of the correlation functions. It works equally well for the trapped andhomogeneous case.In figure 6 we see a comparison of approximations and exact numerical results within extendedmean-field theory as well as Lieb-Liniger theory. In either case we observe a good agreementbetween our results and Lieb-Liniger theory. However as γ increases the deviation from the exactresult grows. We attribute this deviation to the non-Gaussian contributions that have been dropped.Another relevant quantity is the ground state energy of the system. By comparing the value ofthe energy functional (29) with the Lieb-Liniger ground state energy for a range of the correlationparameter γ , we obtain figure 7. In particular, we plot the relative deviation of the ground stateenergies. This is to be compared with deviations from a simple mean-field approach neglectingfluctuations and for the Bogoliubov method in the GP-regime which includes excitations of themean-field. The latter approach results in [56] e ( γ ) LL,GP = γ − π γ . (54)In either case, we present all the results in the form of a normalized deviation from the dimen-sionless ground state energy e ( γ ) from Lieb-Liniger theory given by (6).8 FIG. 6: Second and third order correlation functions versus the correlation parameter γ . In subplot a) wecompare the exact results from Lieb-Liniger theory, g (2) LL (solid line), and the approximation in the GP-regime, g (2) LL,GP (dashed dotted line) to g (2) x,x ≡ g (2)0 , calculated with an extended mean-field theory. Wedepict exact results (dashed line) using (41) and (42) as well as approximated results (dotted line) using(43). In subplot b) we depict the same comparison for the third order correlation function. −4 −3 −2 −1 −8−4048 γ e ( γ ) / e ( γ ) LL − i n % FIG. 7: Relative deviations from the dimensionless ground state energy of Lieb-Liniger theory as a functionof γ . Results from an extended mean-field approach (solid line), the simple mean-field theory withoutfluctuations (dotted line) and the Bogoliubov approach (dashed dotted line) are compared. γ ≈ the maximum deviation of our results isless than 4 % and we obtain reliable results throughout the region of interest, i.e. γ ≤ . However,this appears to be the limit for a quasi one-dimensional extended mean-field theory and differentapproaches have to be used in the strongly correlated regime.0 FIG. 8: Coherent single particle density matrix f ( c ) x,y versus x and y . As the ground state is real valued, thecoherent part of the pairing field m ( c ) x,y is also represented in this figure. V. NUMERICAL RESULTS FOR TRAPPED ATOMS AT ZERO AND FINITE TEMPERATUREA. The zero temperature limit for a trapped gas
In the previous section we have studied the homogeneous case. In here, this will be extendedto harmonically trapped systems and we present correlation functions up to third order. First of allwe depict the spatial shape of the master variables ˜ f , ˜ m and of the quantities f ( c ) , m ( c ) , which areessential for the calculation of the correlation functions. The plots show numerical simulations fora particle number of N = 10 in a trap with standard parameters for Rb according to (30).The coherent contribution to the single particle density matrix f ( c ) x,y in figure 8 has off-diagonallong range order and extends over the complete system. As the Hamiltonian for a one-dimensionaltrap is real-valued, so is the ground state solution α x . Hence, the coherent contribution of thepairing field m ( c ) x,y is identical to f ( c ) x,y and shown in figure 8.In contrast to the coherent contributions, the normal fluctuation ˜ f x,y in figure 9 and the anoma-lous fluctuation ˜ m x,y in figure 10 are primarily localized along the diagonal.1 FIG. 9: Normal fluctuations ˜ f x,y versus x and y . The coherence in the off-diagonal direction is only of short range and the negativity of thepairing field is an indication of a reduced likelihood of finding two particles at the same location.The behaviour of the first, second and third order correlation function is presented in figures11, 12 and 13. A generic feature of all three correlation functions is that they become morepronounced for smaller particle numbers. For the first order correlation function the diagonalhas to be identical to one and the deviation in the off-diagonal is fairly small as expected fora coherent system. However, the second order density-density correlation is a more sensitiveprobe as this correlation function is less than one, thus exhibits non-classical behaviour. Thisanti-bunching is particularly strong for smaller particle numbers when we approach the Tonks-Girardeau regime of a fermionized Bose gas and the correlation function vanishes eventually.Recently this effect has been investigated in a number of experiments, e.g. [5, 6, 14] and confirmsthe theoretical predictions. The same statements apply to the third order correlation function andit can be observed that the deviation from one is even more pronounced. This also implies thatthe third order correlation function [3, 30] is the most sensitive probe for quantum aspects of thefield. In addition we notice values which are clearly below one for | y | (cid:29) and x (cid:54) = y , because inthis case g (3) x,y ≈ g (2) y,y . This can easily be seen by looking at (52,53) and taking into account that allterms with off-diagonal contributions of the fluctuations in g (3) x,y are negligible for x (cid:54) = y .2 FIG. 10: Anomalous fluctuations − ˜ m x,y versus x and y .FIG. 11: First order correlation function g (1) x,y versus x and y . FIG. 12: Second order correlation function g (2) x,y versus x and y .FIG. 13: Third order correlation function g (3) x,y versus x and y . −4 −3 −2 −1 γ −4 −3 −2 −1 γ g ( ) , g ( ) , a)b) FIG. 14: Correlation functions versus the correlation parameter γ . In subplot a) we depict g (2)0 , and comparesimulation results for the trapped case (circles) with analytic calculations obtained with Lieb-Liniger theory(solid line). In subplot b) we depict g (3)0 , and again compare simulation results for the trapped case (circles)with analytic calculations obtained with Lieb-Liniger theory (solid line). In both comparisons the circlesoriginate for particle numbers N ranging from N = 10 on the left hand side to N = 10 on the right handside. B. Behaviour in the center of the trap
In figure 14 we compare the results of our simulations for the second and third order correlationfunction with Lieb-Liniger theory. In contrast to the comparison in subsection IV C an externalpotential is now included in the calculations with the extended mean-field theory whereas thetheoretical curve is for a homogeneous gas of bosons. Our simulations are for particle numbersranging from N = 10 − and we only used the values of the correlation functions in the centerof the trap for the comparison. Compared to subsection IV C the results in the trapped case deviateslightly more from the exact results originating from the homogeneous Lieb-Liniger theory but thequalitative behaviour is very similar.5 FIG. 15: The diagonal second order correlation function g (2) x,x versus γ for various particle numbers. Withan increasing value of the correlation parameter γ the circles correspond to points further outwards fromthe origin. We compare results for the trapped case (circles) with analytic calculations for the second ordercorrelation function obtained with LDA-Lieb-Liniger theory (solid line). Plots a) to c) are for N = 10 , N = 10 and N = 10 (from left to right) and plots d) to f) are for N = 10 , N = 10 and N = 10 (fromleft to right). C. Diagonal behaviour in the local density approximation
The local density approximation (LDA) is a frequently employed approximation scheme totransfer results of homogeneous systems to spatially trapped gases. It is assumed that a smoothvariation of the density profile can be incorporated by an adiabatic adjustment of a locally uniformgas. The LDA uses a local effective chemical potential [38] µ ( x ) = µ − V ( x ) = µ − mω x , (55)where µ denotes the global equilibrium chemical potential. In order for the LDA to be applicable,it is thus necessary that the short-range correlation length is much smaller than the characteristicinhomogeneity length.6 FIG. 16: The diagonal third order correlation function g (3) x,x versus γ for various particle numbers. Withan increasing value of the correlation parameter γ the circles correspond to points further outwards fromthe origin. We compare results for the trapped case (circles) with analytic calculations for the third ordercorrelation function obtained with LDA-Lieb-Liniger theory (solid line). Plots a) to c) are for N = 10 , N = 10 and N = 10 (from left to right) and plots d) to f) are for N = 10 , N = 10 and N = 10 (fromleft to right). In this context, we want to compare the diagonal behaviour of our numerically calculated cor-relation functions to theoretical predictions. By definition, the first order correlation function isidentical to one along the diagonal and our data behaves accordingly. For the second and thirdorder correlation function we will compare our results with the predictions from Lieb-Liniger the-ory in the LDA. Naturally the LDA works best in the center of the trap. It can not be expectedto work in regions where the density drops rapidly and the inhomogeneity length is very small inthese regions.In the Gross-Pitaevskii regime the chemical potential µ connects the density n to the correlationparameter γ , via µ ( x ) = gn ( x ) , γ ( x ) = mg (cid:126) n ( x ) . (56)7In our simulations we tune the particle number in the trap, which decreases γ for an increasingnumber of particles. Qualitatively one can expect that the inhomogeneous correlation functions arehigher than the homogeneous results because in the LDA the external potential leads to a smallerchemical potential and according to (56) also to a smaller density compared to the homogeneouscase. Due to the monotonous decrease of the correlation functions, there is a tendency of the in-homogeneous values to be shifted to larger γ values. All the features that have just been describedcan be seen in figures 15 and 16, where we plotted the correlation functions for particle numbersranging from N = 10 − and restricted the plotted regions to the Thomas-Fermi radius.8 FIG. 17: Finite temperature first order correlation function g (1) x,y versus x and y , for N = 100 and T =10 (cid:126) ω/k B . D. The finite temperature result for a trapped gas
The zero-temperature results of the previous section can be extended easily to account for finitetemperature effects [36, 49]. One obtains an equilibrium solution for the density matrix G of thethermal system (21) from the eigenstates of the selfenergy matrix (25), according to the Bose-Einstein distribution. We present results in the present section for a particle number of N = 100 and a temperature T = 10 (cid:126) ω/k B .The main thermal effect is a strong increase of the fluctuations at the edge of the trap at thecost of a reduction of the condensate density [57]. This effect is clearly seen by comparing thefirst order correlation function in figure 17 to the zero temperature result in figure 11. At finitetemperatures, we also obtain a reduction of first order coherence. Consequently, this leads toa situation where the gas is almost thermalized at the edge of the trap, whereas it is coherentin the center. The suppression of density fluctuations, also known as anti-bunching, is also lesspronounced at finite temperature. This can be seen by comparing figures 18 and 19 to figures 12and 13, which give the zero temperature results.9 FIG. 18: Finite temperature second order correlation function g (2) x,y versus x and y , for N = 100 and T = 10 (cid:126) ω/k B .FIG. 19: Finite temperature third order correlation function g (3) x,y versus x and y , for N = 100 and T =10 (cid:126) ω/k B . xx g ( ) x , − x g ( ) x , − x a)b) FIG. 20: Off-diagonal first g (1) x, − x (a) and second order correlation function g (2) x, − x (b) versus x . The individualcurves correspond to temperatures k B T = 0 (cid:126) ω (smallest value for x = 0 ) to (cid:126) ω (largest value for x = 0 )with increments of (cid:126) ω . For a thermal gas of noninteracting bosons, one finds g (2)0 , = 2! and g (3)0 , = 3! . It can be seenthat these values are attained at the edge of the trap where fluctuations dominate. In figure 19 wealso notice a value of two for | y | (cid:29) and x (cid:54) = y , because we again have g (3) x,y ≈ g (2) y,y = 2 in thiscase.In figure 20, we present the off-diagonal of the first and second-order correlation function g (1 , x, − x versus x for temperatures from k B T = 0 − (cid:126) ω with increments of (cid:126) ω . It can be noticedthat correlations are strongly attenuated with increasing temperature. Looking at the off-diagonalof g (2) x, − x in figure 20, we see a reduction of the anti-bunching dip in the center with increasingtemperatures, however it is still present for high values of the temperature. It can be understoodqualitatively from the stronger increase of fluctuations with the temperature at the edge of the trap.Thus, the anti-bunching dip in the center of the trap remains visible even at finite temperatures.1 VI. CONCLUSIONS AND OUTLOOK
We have presented a detailed study of an extended mean-field theory which shows that it is afeasible approach for the description of a weakly correlated gas of bosons in a quasi-1D setup.This approach agrees well with exact predictions of zero-temperature Lieb-Liniger theory and caneasily be applied to spatially inhomogeneous systems at finite temperature. We have not analyzedthe finite temperature theory of Yang-Yang due its complexity.There are many relevant applications for using this extended mean-field theory in different geo-metrical configurations like a double-well potential [58, 59]. Yet another extension of our approachis the dimensional crossover out-of-equilibrium where in general the increase of available phasespace volume leads to a decrease of correlations. An evaluation of the such correlation functionsis work in progress.
Acknowledgments
The authors acknowledge the support by the German Research Foundation via SFB/TRR 21which is a collaboration of the Universities of Stuttgart, T¨ubingen, Ulm and the Max Planck Insti-tute for Solid State Research in Stuttgart.2
APPENDIX A: HIGHER TRANSCENDENTAL FUNCTIONS1. Complete elliptic integrals
Following the definitions and the notation in [54], the complete elliptic integral of the first kindreads K ( m ) = (cid:90) π/ d θ (cid:112) − m sin θ , (A1)where the parameter ≤ m ≤ . For the calculation of ˜ m in section IV A, we encounter anintegral of the form I = (cid:90) ∞ d k √ k + ω c √ k + ˜ ω = K ( m ) k c (A2)where ω c > ˜ ω . We can show that the evaluation of this integral leads to the complete ellipticintegral of the first kind by making the substitution k = k c cot θ and using m = ( ω c − ˜ ω ) /ω c .Similarly the complete elliptic integral of the second kind is defined as E ( m ) = (cid:90) π/ (cid:112) − m sin θ d θ . (A3)In order to calculate ˜ f in section IV A we end up with the integral I = (cid:90) ∞ d k √ k + ω c − √ k + ˜ ω √ k + ˜ ω (A4)after separating the constant contribution which leads to the previously discussed integral. Thesame substitution as above k = k c cot θ simplifies the integral to the form I = k c (cid:90) π/ − (cid:112) − m sin θ sin θ (cid:112) − m sin θ d θ = − k c ( E ( m ) − K ( m )) . (A5)
2. The Lambert-W function
The Lambert- W -function is implicitly defined by the solution of the transcendental equation[55] z = W e W . (A6)In the case of large arguments z (cid:29) , one can use an asymptotic expansion W ( z ) = (cid:96) − (cid:96) + (cid:96) /(cid:96) + . . . , (A7)with (cid:96) = ln z and (cid:96) = ln ln z .3 FIG. 21: Integration contour for the evaluation of (B1).
APPENDIX B: DEFORMATION OF THE INTEGRATION CONTOUR IN THE COMPLEXPLANE
The results of subsection IV B 2 can be derived alternatively with the help of complex inte-gration [60]. If we take the inverse Fourier transform of the anomalous fluctuation of (38), weget ˜ m ( r ) = − ω − π (cid:90) ∞−∞ e − ikr √ k + ω c √ k + ˜ ω dk . (B1)Using the substitutions k = ˜ kz , r (cid:48) = ˜ kr and b = k c / ˜ k this equation reduces to ˜ m ( r ) = − ω − π (cid:90) ∞−∞ e − ir (cid:48) z √ z + b √ z + 1 dz ˜ k . (B2)For the evaluation of this integral we make a branch cut between − i and − ib and choose the pathof integration as can be seen in Fig. 21.The contributions from C and C vanish if the contour is moved to infinity and the contribu-tions from C and C cancel each other. The integrals along the semicircles around − ib and thecircle around − i tend to zero if the radius tends to zero. Due to the branch cut the contributionsfrom C and C are equal. Thus, the integral to be solved reads I = (cid:90) ∞−∞ e − ir (cid:48) z dz √ z + 1 √ z + b = (cid:90) − i − ib e − ir (cid:48) z dz √ z + 1 √ z + b (B3)4and by changing the variable of integration ( z = − i − iy ), taking into account that b (cid:29) in theGross-Pitaevskii regime, we get I ≈ e − r (cid:48) (cid:90) ∞ e − yr (cid:48) dy (cid:112) y + y (cid:112) b − − y − y . (B4)As we are looking for an approximation for r (cid:48) = ˜ kr (cid:29) , we notice that only small values of y play an important role for the evaluation of the integral.Hence we neglect the expression y + y in the second term in the denominator whichyields I ≈ e − r (cid:48) √ b (cid:90) ∞ e − yr (cid:48) (cid:112) y + y dy = 2 √ b K ( r (cid:48) ) (B5)and as ˜ m is an even function in r we get the final result ˜ m ( r ) ≈ − ω − πk c K (˜ k | r | ) ≈ − k c π K (˜ k | r | ) . (B6) [1] A. G¨orlitz et al. , Phys. Rev. Lett. , 130402 (2001).[2] H. Moritz, T. Stoferle, M. Kohl, and T. Esslinger, Phys. Rev. Lett. , 250402 (2003).[3] B. Laburthe Tolra et al. , Phys. Rev. Lett. , 190401 (2004).[4] D. Hellweg et al. , Phys. Rev. Lett. , 010406 (2003).[5] T. Kinoshita, T. Wenger, and D. S. Weiss, Phys. Rev. Lett. , 190406 (2005).[6] C.-S. Chuu et al. , Phys. Rev. Lett. , 260403 (2005).[7] A.H. van Amerongen et al. , Preprint arXiv:0709.1899 , (2007).[8] S. Hofferberth et al. , Nature , 324 (2007).[9]
The Many-Body Problem: An Encyclopedia of Exactly Solved Models in One Dimension , edited byD. C. Mattis, World Scientific, Singapore (1995).[10] M. Girardeau,
J. Math. Phys. , 516 (1960).[11] E. Lieb and W. Lininger, Phys. Rev. , 1605 (1963).[12] C. N. Yang and C. P. Yang,
J. Math. Phys. , 1115 (1969).[13] K. K. Das, M. D. Girardeau, and E. M. Wright, Phys. Rev. Lett. , 110402 (2002).[14] B. Paredes et al. , Nature , 277 (2004).[15] B. D. Esry,
Phys. Rev. A , 1147 (1997). [16] O. E. Alon, A. I. Streltsov, K. Sakmann, and L. S. Cederbaum, Europhysics Letters , 8 (2004).[17] S. Z¨ollner, H.-D. Meyer, and P. Schmelcher, Phys. Rev. A , 063611 (2006).[18] A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, Phys. Rev. A , 063626 (2006).[19] S. R. White, Phys. Rev. Lett. , 2863 (1992).[20] D. Jaksch et al. , Phys. Rev. Lett. , 3108 (1998).[21] D. van Oosten, P. van der Straten, and H. T. C. Stoof, Phys. Rev. A , 053601 (2001).[22] M. Rigol and A. Muramatsu, Phys. Rev. Lett. , 240403 (2005).[23] C. Kollath, U. Schollw¨ock, J. von Delft, and W. Zwerger, Phys. Rev. A , 031601 (2004).[24] M. Aizenman et al. , Phys. Rev. A , 023612 (2004).[25] P. D. Drummond and P. Deuar, J. Opt. B: Quantum Semiclass. Opt. , S281 (2003).[26] M. Yasuda and F. Shimizu, Phys. Rev. Lett. , 3090 (1996).[27] P. Bouyer and M. Kasevich, prepint (1996).[28] B. Saubam´ea et al. , Phys. Rev. Lett. , 3146 (1997).[29] A. Perrin et al. , Phys. Rev. Lett. , 150405 (2007).[30] E. A. Burt et al. , Phys. Rev. Lett. , 337 (1997).[31] Y. Kagan, B. V. Svistunov, and G. V. Shlyapnikov, JETP Lett. , 209 (1985).[32] D. S. Petrov, G. V. Shlyapnikov, and J. T. M. Walraven, Phys. Rev. Lett. , 3745 (2000).[33] M. Olshanii and V. Dunjko, Phys. Rev. Lett. , 090401 (2003).[34] D. M. Gangardt and G. V. Shlyapnikov, New J. Phys. , 79.1 (2003).[35] D. M. Gangardt and G. V. Shlyapnikov, Phys. Rev. Lett. , 010401 (2003).[36] R. Walser, Opt. Comm. , 107 (2004).[37] N. M. Bogoliubov, C. Malyshev, R. K. Bullough, and J. Timonen,
Phys. Rev. A , 023619 (2004).[38] K. V. Kheruntsyan, D. M. Gangardt, P. D. Drummond, and G. V. Shlyapnikov, Phys. Rev. A , 053615(2005).[39] G. E. Astrakharchik and S. Giorgini, J. Phys. B , S1 (2006).[40] V. V. Cheianov, H. Smith, and M. B. Zvonarev, Phys. Rev. A , 051604 (2006).[41] V. V. Cheianov, H. Smith, and M. B. Zvonarev, cond-mat 0602468 (2006).[42] R. Walser, J. Williams, J. Cooper, and M. Holland, Phys. Rev. A , 3878 (1999).[43] M. Holland, J. Park, and R. Walser, Phys. Rev. Lett. , 1915 (2001).[44] M. Holland, S. Kokkelmans, M. Chiofalo, and R. Walser, Phys. Rev. Lett. , 120406 (2001).[45] R. P. Feynman, Phys. Rev. , 340 (1939). [46] A. I. Akhiezer and S. V. Peletminskii, Methods of Statistical Physics , Pergamon Press Ltd., Oxford,England (1981).[47] D. Zubarev, V. Morozov, and G. R¨opke,
Statistical Mechanics of Nonequilibrium Processes , AkademieVerlag, Berlin (1997).[48] S. Chapman and T. G. Cowling,
The Mathematical Theory of Non-Uniform Gases , Cambridge Uni-versity Press, Cambridge (1970).[49] J.-P. Blaizot and G. Ripka,
Quantum Theory of Finite Systems , The MIT Press, Cambridge, Mas-sachusetts (1986).[50] J. Wachter, R. Walser, J. Cooper, and M. Holland,