Exploring the growth of correlations in a quasi one-dimensional trapped Bose gas

Phase correlations, density fluctuations and three-body loss rates are relevant for many experiments in quasi one-dimensional geometries. Extended mean-field theory is used to evaluate correlation functions up to third order for a quasi one-dimensional trapped Bose gas at zero and finite temperature. At zero temperature and in the homogeneous limit, we also study the transition from the weakly correlated Gross-Pitaevskii regime to the strongly correlated 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-over regime.

Currently, this fact has stimulated many fascinating experiments in the context of ultracold gases [1,2,3,4,5,6,7,8], which explore various aspects of the geometric transitions. Serendipitously, also most exactly solvable models of field theory are one-dimensional [9] and rest on the celebrated Bethe-ansatz invented in the 1930ies. In the context of atomic Bose gases, today most prominent are the spatially homogeneous models of hard-and soft-core bosons on a string of M. Girardeau [10] as well as E. Lieb and W. Liniger [11,12]. One of the many interesting questions which can be explored, is the cross-over from the weakly correlated Gross-Pitaevskii regime (γ 1) to the strongly correlated Tonks-Girardeau regime [13,14] (γ 1). Thereby, one commonly uses the Lieb-Liniger (LL) parameter γ, cf. (5), to measure the relative strength of kinetic to repulsive self-energy in the dilute Bose gas.
Today, we also have alternative tools available: First, there are the exact few-body calculations, i. e. multi-channel time-dependent Hartree-Fock (MCTHF) or configuration interaction (CI) methods, which originate from atomic, molecular and nuclear physics. While originally designed for fermionic energy structure calculation, they are nowadays also applied to few-boson systems (≈ 10-100 particles) in arbitrary trap geometries [15,16,17,18]. Second, there is now the possibility to prepare atomic gases in optical lattices and this opens up the rich methodology of the Bose-Hubbard model and density matrix renormalization group methods [19,20,21,22,23,24].
Third, there are stochastic multi-mode trajectory simulations [25] that also successfully address the same questions.
Irrespective of the choice of method, all need to predict experimentally accessible observables in terms of correlation functions, spatial averages or Fourier transforms, thereof. Most relevant are obviously the lowest order moments of the bosonic field operatorâ x , which is the single particle density n x at position x and the conjugate phase quadrature correlation function g (1) x,y . The fluctuations about the mean density are measured with the second order density-density correlation function g (2) x,y n x = â † xâx , g (1) x,y = â † yâx √ n x n y , g (2) x,y = â † xâ † yâyâx n x n y .
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. This situation has been rectified recently in [40,41] where the diagonal behaviour of this correlation function was calculated in the framework of Lieb-Liniger theory. This is where extended meanfield theory is useful, because we can calculate arbitrary orders of the correlation function and will present calculations of the diagonal and off-diagonal behaviour of the third order correlation function at zero as well as finite temperature. However, the extended mean-field (EMF) approach is restricted to values of the correlation parameter γ ≤ 1, because any mean-field theory is known to 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 is an ideal benchmark for the extended mean-field theory [36,42,43,44], whose basic concepts are summarized in section III. In section IV, we specialize the kinetic equations to a quasi onedimensional homogeneous situation at zero temperature for which analytical solutions can be found and compare correlation functions with the LL-predictions. The extension to inhomogeneous, 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 homogeneous gas of N bosons on a ring of length L. It is one of very few exactly solvable problems in manybody physics and provides a solution for every value of the correlation parameter γ. Even in inhomogeneous trapped systems this is very useful, if we can make the local density approximation. In the language of second quantization, the starting point for Lieb-Liniger theory is the following HamiltonianĤ where m denotes the mass of a boson and the creation and annihilation operators satisfy the usual bosonic commutation relation. With the help of the Hellmann-Feynman theorem [45], one can obtain the diagonal part (x = y = 0) of the translation invariant second order correlation function g (2) x,y = g by differentiating the ground state energy E 0 with respect to the coupling constant g. Here, |Ψ 0 represents the ground state and n = N/L denotes the linear particle density. It was shown by Lieb and Liniger [11] that the ground state energy only depends on the dimensionless correlation parameter γ. It is basically the ratio of the repulsive mean-field energy gn to the kinetic energy 2 /2md 2 at an average distance d = 1/n. Another length scale of the problem is the healing length ξ, which equates the kinetic energy of a wave function at scale ξ to the mean-field energy We call bosons weakly correlated for γ 1 (Gross-Pitaevskii regime) and strongly correlated for In terms of this parameter, the ground state energy and second order correlation function LL = e (γ), are given in terms of the solutions of the Lieb-Liniger equations h(x, γ) = 1 2π   In the weakly correlated Gross-Pitaevskii limit γ → 0, as well as in the strongly correlated Tonks-Girardeau regime γ → ∞, one obtains for the correlation function [35] g (2) LL,T G = 4π 2 3γ 2 , for γ 1.
A comparison of these approximations with the exact solution is presented in figure 1. The validity of these results has recently been tested experimentally over a wide range of the correlation parameter [5] and will be used to probe the extended mean-field approach presented in the next section. It is also possible to derive an exact result for the diagonal part of the third order correlation function within the framework of Lieb-Liniger theory, however the task is considerably more difficult. The exact result is derived in [41] by introducing a new functionẽ(γ) which has the form and with the help of this function one obtains A comparison of the exact result for the third order correlation function with the approximations in the Gross-Pitaevskii and the Tonks-Girardeau regime [35] is presented in figure 1 g (3) LL,T G = 16π 6 15γ 6 , for γ 1.
The auxiliary function h(x, γ) of (8) is depicted for various values of the correlation parameter in figure 2.

A. Time-dependent Hartree-Fock-Bogoliubov equations
The evolution of a weakly interacting dilute gas of bosons in three dimensions can be described by a Hamiltonian where H is the single-particle energy in an external potential V ext and V bin is the two-particle potential. As we are interested in the quasi one-dimensional limit, we will consider a cigar shaped trapping configuration (angular frequencies ω and ω ⊥ ) with a large aspect ratio β. The energy and length scales will be set by the transverse oscillator Conceptually, it is straight forward in the extended mean-field theory to use real, finite range binary interaction potentials and obtain proper two-body T matrices including many-body corrections. However, for convenience, we will use the pseudo-potential approximation in here where a s denotes the s-wave scattering length. In order to compactify the notation we will interchangeably 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 and length scales [46,47] and leads to a rapid attenuation of correlation functions. Mathematically speaking, it allows for a selfconsistent expansion of the full many-body density matrix ρ in terms of a perturbation series of simple many-body density matrices σ (i) , which depend parametrically on the master variables This non-perturbative series in terms of the interaction potential V bin has been introduced first by Chapman and Enskog in the context of kinetic theory of gases [48]. In addition to the simple series expansion, we impose a selfconsistency constraint such that the operatorsγ i , corresponding to the c-number master variables γ i , fulfill As far as the master variables are concerned, we choose the mean-field α x , the normal fluctuations of the single-particle densityf and the fluctuations of the anomalous two-particle correlation functionm such that A Gaussian operator σ (0) {γ} is compatible with the requirements of (18,19,20). In turn, this implies the factorizability of multi operator products (Wick's theorem) and also yields non-Gaussian corrections by calculating the contribution of σ (1) {γ} . By studying the coordinate transformation properties of the fluctuations [49], one finds that the averagesf andm are components of a positive semi-definite, generalized density matrix G ≥ 0 . Thus, the system is described by a row vector χ, containing the mean-field α as well as its complex conjugate, and by the density matrix G It can be shown from a Cauchy-Schwartz inequality that at T = 0, the generalized density matrix G obeys an idem-potency relation Starting with the Heisenberg equation of motion, it is straightforward to derive the equations of motion for χ and G. In order to obtain higher-order correlation functions within the present approximation scheme [42,50], like g (2) or g (3) of (1,2), one has to evaluate the Gaussian Tr[. . . σ {γ} ]. However it is clear that the Gaussian contribution will dominate for weak correlations. Thus, we have evaluated in here only the Gaussian contributions. But already for γ ≈ 1, deviations from that can be noticed.

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 out and only amplitudes proportional to the ground state 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 with the following abbreviations for the single particle Hamiltonian and the self energies In the course of the dimensional reduction, we had to introduce an effective one-dimensional coupling constant g = 2 ω ⊥ a s , which is in agreement with previous derivations [36,51,52]. In order to obtain the stationary solution for the time-independent fields χ x and G x,x , we make the ansatz which introduces the chemical potential µ and employs the Pauli matrix σ 3 of (22). This ansatz implies that the normal fluctuationsf x,x (t) become time-independent, while the anomalous fluctuationsm x,x (t) oscillate with twice the chemical potential. The properties of the resulting stationary HFB equation will be investigated in section IV in the case of a homogeneous gas of bosons and in 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) = Ĥ given by with the constraint that the number of particles N = dx(f (c) x,x +f x,x ).

IV. ANALYTIC SOLUTION FOR THE STATIONARY HFB-EQUATIONS IN THE HOMOGENEOUS SYSTEM AT ZERO TEMPERATURE
For the calculations that are presented in the following sections, we use standard parameters for 87 Rb in natural units of length a ⊥ and energy ε ⊥ m = 1.4432 · 10 −25 kg, a s = 5.8209 · 10 −9 m, To reach the homogeneous case, we have to decrease the influence of the external trapping potential in (27) by weakening it to the limit β 1, where we can neglect it practically. Therefore, the equilibrium state should possess the same translation symmetry as the generator of the dynamics. Consequently, we can assume that the mean field is space independent α x = α and the density matrix only depends on relative differences Translation invariant systems are best described in Fourier-space, which was introduced above. We can also choose the mean-field to be real-valued by a suitable phase rotation. This is a consequence of the global number conservation that is built into the dynamical HFB equations [53]. As the mean-field is real-valued f (c) = m (c) = α 2 , so are the fluctuationsf 0 =f x,x andm 0 =m x,x and with these assumptions, the normalization constraint reads where n is the linear particle density on a length L. Furthermore the THFB equations (24) simplify significantly to 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 In particular, in the Gross-Pitaevskii regime one can assume that ω c ω. With these definitions, one finds that the self energy is simply a 2 × 2 matrix in k-space with eigenvalues ω k Now, we can finally evaluate the density matrix part of the HFB equations (33). Moreover, we also have to consider the idem-potency relation of (22). It holds for the vacuum state at zero temperature and one obtains another, now quadratic relation between normal and anomalous fluctuations in k-spacem The system of equations can be solved point wise in k-space and leads to two solutions. One of which has to be rejected on physical grounds. Thus, we find The high momentum tail of the correlation functions is responsible for the short scale behaviour in real space. In the limit k → ∞ the leading terms arẽ

A. Diagonal contributions of normal and anomalous fluctuations
In order to obtain the diagonal part of the translation invariant correlation functionsm r andf r at r = 0, we have to evaluate the inverse Fourier transform of (31) Serendipitously, this can be done exactly in terms of elliptic integrals [54] where K and E are the complete elliptic integral of the first kind and second kind, respectively.
Basic definitions are given in A 1.
In this explicit formula, we had to introduce the Lambert-W -function, which is defined in A 2 and an excellent asymptotic expansion is given in terms of logarithms [55].
The approximations for the fluctuations are compared to exact numerical calculations in figure   3 and give a good agreement. We will use these approximations in the following sections to evaluate the ground state energy and correlation functions. forms, which has its origin in (37). Starting withf (0) k = 0 and using the recursion relatioñ we get a rapid convergence towards the exact results. It is remarkable that even with the inverse Fourier transforms of low orders of this iteration scheme, we get a functional behaviour forf (r) andm(r), which is equivalent to their exact behaviour for short ranges. However in contrast to the exact form of the Fourier transforms of the fluctuations in (38), it is possible to perform the inverse Fourier transform analytically in every order of the iteration scheme. A closer look reveals that the dependence of the fluctuations on r has to be of the form where P (i) (r) and Q (i) (r) are polynomials in r of order 2 i − 2 and 2 i+1 − 3 respectively. Consequently the length scale on which the correlations decay is given by In the Gross-Pitaevskii regime, where f (c) m 0 , we recover the healing length ξ ≈ 1/ 2gf (c) , already introduced at the beginning in (5).
The short range behaviour of the anomalous fluctuation and the normal fluctuation is depicted in figure 4. There, we compare the 4th order result of the iteration scheme to the exact numerical evaluation of the inverse Fourier transform. We assumed N = 100 particles, distributed over a length of L = 90a ⊥ . This length was chosen such that the density in the homogeneous case is similar to the density in the center of the trapped system, which will be discussed in section V. One obtains a good agreement between the approximation and the exact results in the regime where r ξ ≈ 10.5 with ξ ≈ 3.88. At the origin we note that the anomalous fluctuation shows the typical cusp whereas the normal fluctuation has a smooth behaviour and consequently a vanishing first derivative at r = 0.

Long length scale behaviour: r ξ
In order to get an approximation for the fluctuations in this regime, we start with the Fourier transform of the anomalous fluctuation in (38) and note for further consideration that the Fourier transform of the modified Bessel function of the second kind K 0 (c|r|) is given by π/ √ k 2 + c 2 .
Thus the convolution property of the Fourier transform yields The modified Bessel function of the second kind K 0 (r) is presented in figure 5. K 0 (r) diverges logarithmically at the origin and decreases exponentially for large arguments where γ e ≈ 0.5772 denotes Euler's constant. In the Gross-Pitaevskii regime where f (c) |m 0 |, we get from r ξ that also r ξ c and therefore the first Bessel function in the integral closely resembles a δ-function. Thus we obtain the resultm An alternative derivation of this result with the help of complex integration is given in B. In figure   4 the asymptotic form of the fluctuations in terms of the Bessel functions is compared to an exact numerical simulation with the parameters that were mentioned in the previous subsection. In either case the semi-logarithmic plot reveals an exponential decay for large distances r ξ , in agreement 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 a quantitative comparison with the results of Lieb-Liniger theory which provides an exact solution for 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 from multiple operator averages. Ifô is such a general operator then While we have already evaluated the Gaussian and non-Gaussian averages for the multinomial operator averages [42], it is clear that the Gaussian contribution will dominate for weak correlations.
Therefore we will focus in here on the Gaussian contribution and disregard the non-Gaussian contributions in the following explicit expressions of order two and one, respectively x,y +f x,y √ n x n y + O(g 2 ), (51) g (2) x,y = 1 + x,y * m y,x ) +f x,yfy,x +m x,ym * y,x x,xfy,y +f x,xfy,y + 2f x,yfy,x + 2m * x,ymy,x y,ym * y,y +m y,ym * y,y + f (c) y,yfy,y where n x = f (c) x,x +f x,x denotes the total density. This way we can calculate the full diagonal and off-diagonal behaviour of the correlation functions. It works equally well for the trapped and homogeneous case.
In figure 6 we see a comparison of approximations and exact numerical results within extended mean-field theory as well as Lieb-Liniger theory. In either case we observe a good agreement between our results and Lieb-Liniger theory. However as γ increases the deviation from the exact result 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 of the energy functional (29) with the Lieb-Liniger ground state energy for a range of the correlation parameter γ, we obtain figure 7. In particular, we plot the relative deviation of the ground state energies. This is to be compared with deviations from a simple mean-field approach neglecting fluctuations and for the Bogoliubov method in the GP-regime which includes excitations of the mean-field. The latter approach results in [56] In either case, we present all the results in the form of a normalized deviation from the dimensionless ground state energy e(γ) from Lieb-Liniger theory given by (6).
LL (solid line), and the approximation in the GPregime, g LL,GP (dashed dotted line) to g (2) x,x ≡ g (2) 0,0 calculated with an extended mean-field theory. We depict 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. The results show a clear improvement over simple mean-field theory and it also improves on the Bogoliubov method. Up to the cross-over at γ ≈ 1 the maximum deviation of our results is less than 4% and we obtain reliable results throughout the region of interest, i.e. γ ≤ 1. However, this appears to be the limit for a quasi one-dimensional extended mean-field theory and different approaches have to be used in the strongly correlated regime.
FIG. 8: Coherent single particle density matrix f (c) x,y versus x and y. As the ground state is real valued, the coherent 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 TEMPERATURE
A. The zero temperature limit for a trapped gas In the previous section we have studied the homogeneous case. In here, this will be extended to harmonically trapped systems and we present correlation functions up to third order. First of all we depict the spatial shape of the master variablesf ,m and of the quantities f (c) , m (c) , which are essential for the calculation of the correlation functions. The plots show numerical simulations for a particle number of N = 10 2 in a trap with standard parameters for 87 Rb according to (30).
The coherent contribution to the single particle density matrix f (c) x,y in figure 8 has off-diagonal long range order and extends over the complete system. As the Hamiltonian for a one-dimensional trap is real-valued, so is the ground state solution α x . Hence, the coherent contribution of the pairing 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 fluctuationf x,y in figure 9 and the anomalous fluctuationm x,y in figure 10 are primarily localized along the diagonal. The coherence in the off-diagonal direction is only of short range and the negativity of the pairing 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 figures 11, 12 and 13. A generic feature of all three correlation functions is that they become more pronounced for smaller particle numbers. For the first order correlation function the diagonal has to be identical to one and the deviation in the off-diagonal is fairly small as expected for a coherent system. However, the second order density-density correlation is a more sensitive probe as this correlation function is less than one, thus exhibits non-classical behaviour. This anti-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 confirms the theoretical predictions. The same statements apply to the third order correlation function and it can be observed that the deviation from one is even more pronounced. This also implies that the third order correlation function [3,30] is the most sensitive probe for quantum aspects of the field. In addition we notice values which are clearly below one for |y| 1 and x = y, because in this case g x,y ≈ g y,y . This can easily be seen by looking at (52,53) and taking into account that all terms with off-diagonal contributions of the fluctuations in g (3) x,y are negligible for x = y. x,y versus x and y. x,y versus x and y. x,y versus x and y.
x,x versus γ for various particle numbers. With

C. Diagonal behaviour in the local density approximation
The local density approximation (LDA) is a frequently employed approximation scheme to transfer results of homogeneous systems to spatially trapped gases. It is assumed that a smooth variation of the density profile can be incorporated by an adiabatic adjustment of a locally uniform gas. The LDA uses a local effective chemical potential [38] µ where µ 0 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 characteristic inhomogeneity length. x,x versus γ for various particle numbers. With In this context, we want to compare the diagonal behaviour of our numerically calculated correlation functions to theoretical predictions. By definition, the first order correlation function is identical to one along the diagonal and our data behaves accordingly. For the second and third order correlation function we will compare our results with the predictions from Lieb-Liniger theory in the LDA. Naturally the LDA works best in the center of the trap. It can not be expected to work in regions where the density drops rapidly and the inhomogeneity length is very small in these regions.
In the Gross-Pitaevskii regime the chemical potential µ connects the density n to the correlation parameter γ, via In our simulations we tune the particle number in the trap, which decreases γ for an increasing number of particles. Qualitatively one can expect that the inhomogeneous correlation functions are higher than the homogeneous results because in the LDA the external potential leads to a smaller chemical potential and according to (56) also to a smaller density compared to the homogeneous case. Due to the monotonous decrease of the correlation functions, there is a tendency of the inhomogeneous values to be shifted to larger γ values. All the features that have just been described can be seen in figures 15 and 16, where we plotted the correlation functions for particle numbers ranging from N = 10 0 − 10 5 and restricted the plotted regions to the Thomas-Fermi radius.
x,y versus x and y, for N = 100 and T = 10 ω/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 finite temperature effects [36,49]. One obtains an equilibrium solution for the density matrix G of the thermal 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 ω/k B .
The main thermal effect is a strong increase of the fluctuations at the edge of the trap at the cost of a reduction of the condensate density [57]. This effect is clearly seen by comparing the first order correlation function in figure 17 to the zero temperature result in figure 11. At finite temperatures, we also obtain a reduction of first order coherence. Consequently, this leads to a situation where the gas is almost thermalized at the edge of the trap, whereas it is coherent in the center. The suppression of density fluctuations, also known as anti-bunching, is also less pronounced at finite temperature. This can be seen by comparing figures 18 and 19 to figures 12 and 13, which give the zero temperature results.
x,−x (a) and second order correlation function g (2) x,−x (b) versus x. The individual curves correspond to temperatures k B T = 0 ω (smallest value for x = 0) to 10 ω (largest value for x = 0) with increments of 2 ω.
For a thermal gas of noninteracting bosons, one finds g 0,0 = 2! and g 0,0 = 3!. It can be seen that these values are attained at the edge of the trap where fluctuations dominate. In figure 19 we also notice a value of two for |y| 1 and x = y, because we again have g x,y ≈ g (2) y,y = 2 in this case.
In figure 20, we present the off-diagonal of the first and second-order correlation function g (1,2) x,−x versus x for temperatures from k B T = 0 − 10 ω with increments of 2 ω. It can be noticed that correlations are strongly attenuated with increasing temperature. Looking at the off-diagonal of g (2) x,−x in figure 20, we see a reduction of the anti-bunching dip in the center with increasing temperatures, however it is still present for high values of the temperature. It can be understood qualitatively 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.

VI. CONCLUSIONS AND OUTLOOK
We have presented a detailed study of an extended mean-field theory which shows that it is a feasible 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 can easily be applied to spatially inhomogeneous systems at finite temperature. We have not analyzed the finite temperature theory of Yang-Yang due its complexity.
There are many relevant applications for using this extended mean-field theory in different geometrical configurations like a double-well potential [58,59]. Yet another extension of our approach is the dimensional crossover out-of-equilibrium where in general the increase of available phase space volume leads to a decrease of correlations. An evaluation of the such correlation functions is work in progress.
APPENDIX A: HIGHER TRANSCENDENTAL FUNCTIONS

Complete elliptic integrals
Following the definitions and the notation in [54], the complete elliptic integral of the first kind reads K(m) = π/2 where the parameter 0 ≤ m ≤ 1. For the calculation ofm 0 in section IV A, we encounter an integral of the form where ω c >ω. We can show that the evaluation of this integral leads to the complete elliptic integral 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 In order to calculatef 0 in section IV A we end up with the integral after separating the constant contribution which leads to the previously discussed integral. The same substitution as above k = k c cot θ simplifies the integral to the form

The Lambert-W function
The Lambert-W -function is implicitly defined by the solution of the transcendental equation [55] z = W e W .
In the case of large arguments z 1, one can use an asymptotic expansion W (z) = 1 − 2 + 2 / 1 + . . . , with 1 = ln z and 2 = ln ln z. The results of subsection IV B 2 can be derived alternatively with the help of complex integration [60]. If we take the inverse Fourier transform of the anomalous fluctuation of (38), we Using the substitutions k =kz, r =kr and b 2 = k 2 c /k 2 this equation reduces tõ For the evaluation of this integral we make a branch cut between −i and −ib and choose the path of integration as can be seen in Fig. 21.
The contributions from C 1 and C 6 vanish if the contour is moved to infinity and the contributions from C 2 and C 5 cancel each other. The integrals along the semicircles around −ib and the circle around −i tend to zero if the radius tends to zero. Due to the branch cut the contributions from C 3 and C 4 are equal. Thus, the integral to be solved reads and by changing the variable of integration (z = −i − iy), taking into account that b 1 in the Gross-Pitaevskii regime, we get I ≈ 2e −r ∞ 0 e −yr dy 2y + y 2 b 2 − 1 − 2y − y 2 . (B4) As we are looking for an approximation for r =kr 1, we notice that only small values of y play an important role for the evaluation of the integral.
Hence we neglect the expression 1 + 2y + y 2 in the second term in the denominator which yields and asm is an even function in r we get the final result