Local correlations in a strongly interacting 1D Bose gas

We develop an analytical method for calculating local correlations in strongly interacting 1D Bose gases, based on the exactly solvable Lieb-Liniger model. The results are obtained at zero and finite temperatures. They describe the interaction-induced reduction of local many-body correlation functions and can be used for achieving and identifying the strong-coupling Tonks-Girardeau regime in experiments with cold Bose gases in the 1D regime.


I. INTRODUCTION
Recently, one-dimensional (1D) Bose gases have been created in long cylindrical traps by tightly confining the transverse motion of particles to zero-point oscillations [1,2,3]. In the present stage, one of the main goals is to achieve the strong-coupling Tonks-Girardeau (TG) regime [4,5], where due to repulsion between particles the 1D Bose gas acquires fermionic properties. These studies have revived an interest in dynamical and correlation properties of 1D gases. In most cases, correlations and dynamics of trapped atomic gases in the 1D regime can be investigated [6,7,8,9,10,11,12] by using the Lieb-Liniger model [5], which assumes that particles interact via a delta-function repulsive potential. Of particular importance are local 3-body and 2-body correlations [11,12] as they govern the rates of inelastic processes, such as 3-body recombination and photoassociation in pair collisions. The measurement of these rates provides a way to identify the strong-coupling Tonks-Girardeau regime, and the intrinsic decay process of recombination in 3-body collisions is crucial for the stability of the gas.
The Lieb-Liniger model for a uniform system of N bosons belongs to the class of integrable models of statistical physics [13]. The ground state and spectrum of elementary excitations for this model were found by Lieb and Liniger [5] and a theory for finding thermodynamic functions at finite temperatures was constructed by C.N. Yang and C.P. Yang [14]. These quantities are determined by relatively simple integral equations for the distribution of quantum numbers (quasi-momenta or rapidities). On the other hand, the problem of correlation properties is far from being completely resolved, except for some correlation functions in the limiting cases of weak and strong interactions. For example, the case of infinitely strong interactions is to a certain extent equivalent to that of free fermions and the interactions play the role of Pauli principle [4]. In this limit, any correlation function of the density is given by the corresponding expression for fermions, and the one-particle density matrix is the determinant of a certain integral operator [15,16]. The expressions for the one-body and two-body correlations for an arbitrary interaction strength were obtained by using the Inverse Scattering Method [17]. However, closed analytical results can be found only as perturbative expansions in the limiting cases of weak and strong interactions [18,19,20].
In this paper, we calculate local many-body correlations, that is the expectation value of a product of 2m field operators (m = 2, 3, . . .) at zero space and time separation. The locality allows us to obtain closed analytical results at zero and finite temperatures in the limiting cases of strong and weak interactions. The paper is organized as follows. In section II we present the Hamiltonian and introduce relevant parameters of the problem. Section III contains the main result of the paper, the derivation of local correlation functions at zero and finite temperature in the strongcoupling regime. In Section IV we use this result for the particular case of two-body and three-body correlations. For the sake of completeness, local correlations for the weak-coupling regime, following from the Bogoliubov approach, are presented in Section V. We conclude in Section VI.

II. INTERACTING BOSE GAS IN 1D
We consider the system of N bosons on a 1D ring of length L in the thermodynamic limit N, L → ∞, with a fixed density n = N/L. The particles interact via a repulsive delta-function interaction potential and are described by the Hamiltonian Here M is the particle mass, x j is the coordinate of the j-th particle, and c = mg/h 2 , with g > 0 being the coupling constant. The quantity c is the inverse interaction length for the 2-body problem with the delta-function potential gδ(x). In other words, the wave function of two particles decreases on a distance scale 1/c as they approach each other. The Hamiltonian (1) is diagonalized by means of the Bethe Ansatz [5]. The many-body wave function is symmetric with respect to permutation of particle coordinates, and in the domain 0 < x 1 < . . . < x N < L the wave function reads In Eq. (2) the sum runs over N ! permutations P acting on N indices of quantum numbers (quasi-momenta) k j . These quantum numbers are chosen such that k 1 < k 2 < . . . < k N . They determine the amplitudes of each term in Eq. (2) and their quantization in a finite system follows from the requirement of periodicity of the wave function (2). The corresponding eigenenergy is E = k 2 j , and N (c) in Eq. (2) is the normalization constant. The latter has been calculated in Refs. [22,23].
The key parameter of the system is the ratio of the mean interparticle separation 1/n to the interaction length 1/c: At sufficiently low temperatures, for γ ≫ 1 the gas is in the strong-coupling regime: the many-body wave function strongly decreases at interparticle distances much smaller than 1/n. In the extreme case γ = ∞ the amplitudes (3) are determined by the sign of the permutation: So, the wave function (2) vanishes when at least two particle coordinates coincide. In this limit any correlation function of the density operators can be calculated by using the theory of free fermions, since in the domain 0 < x 1 < x 2 < . . . < x N < L the wave function (2) is the Slater determinant constructed out of the plane waves with momenta k j . In the case of large but finite γ, the calculations can be performed perturbatively as expounded in the next Section for the local density correlations.
In the opposite limit, γ ≪ 1, the gas is in the weak-coupling regime. In this case the ground state of the system is well described by the Gross-Pitaevskii mean-field theory, and correlation functions at T = 0 can be found on the basis of the Bogoliubov-Popov approach (see [21]). This approach also covers a significant part of the weak-coupling regime at finite temperatures.
For finite temperatures, aside from 1/c and 1/n, one has another length scale, the thermal de Broglie wavelength Λ = (2πh 2 /mT ) 1/2 . The latter is especially important at temperatures T > T d , where T d =h 2 n 2 /2m is the temperature of quantum degeneracy. The relation between the three length scales determines various finite-temperature regimes of the 1D Bose gas, which we discuss below in the limits of strong and weak interactions.

III. LOCAL CORRELATIONS IN THE STRONG-COUPLING LIMIT
We now consider the strong coupling limit γ ≫ 1 and calculate the local m-particle correlation function where Ψ † (x),Ψ(x) are the creation and annihilation operators of bosons. At T = 0 the expectation value . . . is taken with respect to the ground state of the system, and at finite temperatures we assume the average in the grandcanonical ensemble. We first present the derivation at zero temperature and then generalize the method to the case of finite temperatures. For T = 0, in first quantization the correlation function g m (6) reads where the ground state wave function Φ (γ) 0 is given by by the Bethe Ansatz expression (2). In the strong-coupling limit the amplitudes a(P ) (3) can be expanded in inverse powers of the interaction strength The quasi-momenta are close to their ground state values at γ = ∞, given by the uniform Fermi-Dirac distribution in the interval [−k F , k F ], with the Fermi momentum k F = πn. One then sees that in fact Eq. (8) gives the expansion in inverse powers of γ. For γ ≫ 1, we extract the leading behavior of the wave function Φ γ 0 at m coinciding points by appropriately symmetrizing the amplitudes a(P ) for each permutation: where the sum runs over m! permutations p of numbers 1, 2, . . . , m, and is the Vandermonde determinant. Up to the leading term in 1/c, the ground state wave function at m coinciding points is therefore given by This allows us to express Φ Returning to Eq. (7), one sees that the local correlation function is given by derivatives of the m-particle correlation function of free fermions at x 1 = . . . = x m = y 1 = . . . = y m = 0: The operators ψ † , ψ are now fermionic field operators with canonical anticommutation relations. Using Wick's theorem we express the expectation value of these operator as a product of one-particle Green functions: The one-particle Green function is given by where N 0 (k) is the ground state Fermi-Dirac distribution. Then, from Eqs. (15), (14) and (13) we obtain where ∆ m (∂ x ) stands for ∆ m (∂ x1 , . . . , ∂ xm ). This expression can be further simplified as the Vandermonde determinant is a totally antisymmetric function of its variables. Therefore, reordering y pj arguments cancels the factor (−1) p . Using this fact, acting with derivatives on the corresponding arguments in the exponents, and noting that all permutations give the same result, we arrive at the equation: The developed method was briefly outlined in Ref. [11] for 3-particle local correlations. It is readily generalized to the case of finite temperatures. In the grand canonical ensemble the expression for the m-particle local correlation function g m (6) reads: where Z is the partition function, z = exp(µ/T ) is the fugacity, µ is the chemical potential, and the index α labels eigenstates of the system. Expanding the eigenfunctions Φ (γ) α at coinciding points exactly as above in the case of the ground state wave function, we express the leading contribution to g m through derivatives of the finite-temperature correlation function of free fermions at x 1 = x 2 = . . . = x m = y 1 = . . . = y m = 0: The finite-temperature m-fermion correlation function is again calculated with the use of Wick's theorem. This gives Eq. (14) in which G(x i , y Pi ) are now finite-temperature Green functions for free fermions. We then proceed in the same way as at T = 0 and arrive at Eq. (17), with N 0 (k) replaced by the Fermi-Dirac distribution at finite temperatures: At temperatures T ≫ T d the characteristic momentum of particles is the thermal momentum k T ∼ 1/Λ. Therefore, the small parameter for the expansion of the amplitudes a(P ) in Eq. (8) is 1/Λc. Thus, one must satisfy the inequality Λc ≫ 1, which requires temperatures From this point on we consider together the cases of zero and finite temperature. In the thermodynamic limit, expressing the momenta in terms of the size of the Fermi zone, k = 2k F x = 2πnx, we have where the m-fold integral I m (Λn, z) depends on the ratio of the de Broglie wavelength to the mean interparticle separation and on the fugacity z (or chemical potential), and is given by where N (x) is given by Eq. (20) with k = 2πnx. In order to keep the density constant one must fix z as a function of Λn by the normalization condition Under this condition the integral in Eq. (23) becomes a function of Λn only. For calculating this function we employ the method of orthogonal polynomials used in Random Matrix theory [24]. Consider a set of polynomials P j (x) = x j + . . . for j = 0, 1, 2, . . ., orthogonal with the weight N (x) on an infinite interval: Then the value of the integral (23) is expressed through the normalization coefficients h j as follows: At T = 0 we have Λn = ∞. In this limit the distribution function N (x) is uniform in the interval −1/2 < x < 1/2, and is zero otherwise. Polynomials P j (x) can be expressed in terms of Jacobi polynomials [25] and we obtain From Eq. (22) we then obtain the local correlation function where the coefficients A m satisfy the recurrence relation Although general expressions for the normalization coefficients h j in Eq. (25) are not known, for sufficiently small m these coefficients can be calculated straightforwardly as they are related to the moments of the distribution function N (x). For example, one finds We finally consider high temperatures T ≫ T d that still satisfy Eq. (21). In this case the fugacity is z = Λn, and the distribution function is Gaussian: N (x) = z exp (−πΛ 2 n 2 x 2 ). The polynomials P j (x) are Hermite polynomials, and we find , where the recurrence relation for the coefficients B m reads B m+1 = (m + 1)Γ(m + 2) B m , and B 1 = 1.

IV. THREE-BODY AND TWO-BODY CORRELATIONS
In this Section we use the general results of Section III for the particular case of three-body local correlations in the strong-coupling limit. Two-body correlations have been discussed in Ref. [12] and are expounded here for completeness.
At temperatures T ≪ T d we have Λn ≫ 1, and the local correlation functions g 2 and g 3 are close to their zero temperature values. Thus, the system remains in the Tonks-Girardeau regime. For calculating g 2 and g 3 we use Eqs. (22) and (26), with normalization coefficients h 1 and h 2 following from Eq. (30). The quantities < x 2 > and < x 4 > at finite T are obtained on the basis of the Sommerfeld expansion: (33) This gives the following expressions: The results of Eqs. (34) and (35) have a clear physical explanation. For T ≪ T d the characteristic momentum of particles is of the order of k F = πn. Fermionic correlations are present at interparticle distances x > ∼ 1/c, since for smaller distances the correlations do not change. Therefore, g 2 is nothing else than the pair correlation function for free fermions at a distance x ∼ 1/c [12]. This function is of the order of (k F /c) 2 ∼ (π/γ) 2 , which agrees with Eq. (34). With the same arguments, one finds that the 3-body correlation function is g 3 ∼ (k F /c) 2 ∼ (π/γ) 6 , which coincides with the result of Eq. (35).
For temperatures in the interval T d ≪ T ≪ γ 2 T d , where 1/c ≪ Λ ≪ 1/n, we use Eq. (32) directly and obtain: This regime can be called the regime of "high-temperature fermionization" [12]. The gas is no longer degenerate, but due to strong interaction between particles the correlations still have a fermionic character. With regard to local correlations, the main difference from the low temperature case is related to the value of the characteristic momentum of particles. It is now of the order of 1/Λ, instead of k F at T ≪ T d . Accordingly, the two-body and three-body local correlation functions are g 2 ∼ (1/Λc) 2 and g 3 ∼ (1/Λc) 6 , as described by Eqs. (36) and (37).

V. WEAK-COUPLING LIMIT
In the weak-coupling limit, where γ ≪ 1, one can rely on the mean-field approach. The mean-field interaction energy per particle is proportional to ng and it is reasonable to introduce the correlation length l c =h/ √ mng. At temperatures T ≪ T d , over a wide range of parameters the correlation length l c ≪ l φ , where l φ is the phase coherence length. Then the equilibrium state is a quasicondensate, that is the state in which density fluctuations are suppressed but the phase still fluctuates [26]. A review of earlier studies of 1D Bose gases in the weak-coupling limit may be found in Ref. [21]. The zero temperature phase coherence length is always larger than l c , and at finite T one has l φ =h 2 n/mT . Therefore, the condition l c ≪ l φ can be rewritten as One then sees that the quasicondensate regime requires sufficiently low temperatures or relatively large interaction between particles. Due to the presence of phase coherence on a long distance scale, local correlation functions in the quasicondensate regime are close to n 2 . Deviations from this value can be calculated straightforwardly by using the Bogoliubov approximation [27]. For g 2 the results are known (see [12,28] and references therein), and g 3 at T = 0 is given in Ref. [11]. Here we complete the picture and present general results for g m . They are obtained in the same way as g 2 in Ref. [12].
We represent the bosonic field operator as a sum of a macroscopic component Ψ 0 and a small component Ψ ′ describing finite-momentum excitations. To be more precise, the component Ψ 0 contains the contribution of excitations with momenta k < ∼ k 0 ≪ l −1 c , whereas Ψ ′ includes the contribution of larger k. At the same time, the momentum k 0 is chosen such that most of the particles are contained in the part Ψ 0 . This picture is along the lines of Ref. [21], and the momentum k 0 drops out of the answer as the main contribution of the excitation part Ψ ′ to local correlation functions is provided by excitations with k ∼ l −1 c . Confining ourselves to the terms quadratic in Ψ ′ and taking then into account that |Ψ 0 | 2m = n m − mn m−1 Ψ ′ † Ψ ′ , the m-particle local correlation function is given by The normal and anomalous averages, Ψ ′ † Ψ ′ and Ψ ′ Ψ ′ , can be calculated by using the same Bogoliubov transformation for Ψ ′ as in the 3D case: where a k and a † k are annihilation and creation operators of elementary excitations, u k , v k = (ε k ± E k )/2 √ ε k E k , ε k = E 2 k + 2ngE k is the Bogoliubov excitation energy, and E k =h 2 k 2 /2m. This immediately gives g m n m = 1 + m(m − 1) withÑ k = [exp (ε k /T ) − 1] −1 being the equilibrium occupation numbers for the Bogoliubov excitations. The integral term in Eq. (41) contains the contribution of both vacuum and thermal fluctuations. The former is determined by excitations with k ∼ l −1 c , and at T = 0 we find At temperatures T ≪ ng ∼ γT d , thermal fluctuations provide only a small correction (π √ γ/24)(T /γT d ) 2 to the result of Eq. (42). For temperatures T ≫ ng ∼ γT d (but still T ≪ √ γT d ), the contribution of thermal fluctuations to the integral term in Eq. (41) is the most important. It comes from excitation energies ε k ∼ ng and, hence, one may putÑ (k) = T /ε k in the integrand. This yields One clearly sees from Eq. (43) that g m increases with decreasing interaction strength or increasing temperature. The ratio T /T d √ γ is simply l c /l φ . For l c ∼ l φ or T ∼ T d √ γ, the gas leaves the quasicondensate regime and Eq. (43) is no longer valid. At a temperature and interaction strength satisfying the condition T ≫ T d √ γ the gas is in the decoherent regime where the interaction between particles is much less important. Note that for a very small interaction strength this regime is present even at temperatures far below T d . The cross-over from the quasicondensate to the decoherent regime was discussed in Refs. [12,28] on the basis of calculations for g 2 . In the decoherent regime the local correlation function is close to the result for an ideal Bose gas following from Wick's theorem: g m /n m = m!.
The corrections to this value can be calculated perturbatively [12]. The lower is T , the smaller is γ at which the gas enters the decoherent regime. For T = 0 this transition is discontinuous and occurs at zero interaction strength.

VI. CONCLUDING REMARKS
From a general point of view, local correlations of the 1D system are less universal than the long-wave behavior. Local correlation functions depend on a particular model used for calculations. Our approach leads to physically transparent analytical results for g m of the 1D Bose gas in the Lieb-Liniger model.
The use of this model for trapped atomic gases in the 1D regime is justified by the short-range character of interatomic interaction: the characteristic radius of interactions is much smaller than any length-scale of the Lieb-Liniger model. For an atomic gas in an infinitely long cylindrical trap, harmonically confined with frequency ω 0 in the transverse direction to zero-point oscillations, the coupling constant g is expressed through the 3D scattering length a and reads [29] g = 2h 2 a/M l 2 0 , where l 0 = h/M ω 0 is the amplitude of zero point oscillations. Then, to describe the system by the 1D Hamiltonian (1), it is sufficient to satisfy the inequalities l 0 << 1/n, Λ T [11,12].
The reduction of g 3 in the strong-coupling limit, which follows from our results, is important in two aspects. First, it indicates a possibility of achieving this limit at a high gas density (large number of particles), since the rate of decay due to 3-body recombination is proportional to g 3 and will be strongly suppressed. Second, the measurement of 3-body losses of particles can be used for identifying the strong-coupling Tonks-Girardeau regime and the regime of high-temperature fermionization. The identification of these regimes can also be provided by the measurement of photoassociation in pair atomic collisions as the rate of this process is proportional to the two-particle local correlation function g 2 . In this respect, it is worth mentioning that recent Monte Carlo calculations of g 2 at T = 0 for the number of particles as low as 100 [30] agree with our results obtained in the thermodynamic limit.