A Computation of the correlation function

# Quasi-long-range order in trapped two-dimensional Bose gases

## Abstract

We study the fate of algebraic decay of correlations in a harmonically trapped two-dimensional degenerate Bose gas. The analysis is inspired by recent experiments on ultracold atoms where power-law correlations have been observed despite the presence of the external potential. We generalize the spin wave description of phase fluctuations to the trapped case and obtain an analytical expression for the one-body density matrix within this approximation. We show that algebraic decay of the central correlation function persists to lengths of about 20% of the Thomas–Fermi radius. We establish that the trap-averaged correlation function decays algebraically with a strictly larger exponent weakly changing with trap size and find indications that the recently observed enhanced scaling exponents receive significant contributions from the normal component of the gas. We discuss radial and angular correlations and propose a local correlation approximation which captures the correlations very well. Our analysis goes beyond the usual local density approximation and the developed summation techniques constitute a powerful tool to investigate correlations in inhomogeneous systems.

###### pacs:
05.70.Jk, 64.60.an, 67.85.-d

With the achievement of quantum degeneracy in ultracold atom gases a new experimental platform for studying fundamental questions related to phase transitions and critical phenomena has appeared. In particular, correlation functions can be measured through interference and time-of-flight techniques Polkovnikov et al. (2006); Hadzibabic et al. (2006); Cladé et al. (2009); Tung et al. (2010); Plisson et al. (2011); Singh and Mathey (2014); Murthy et al. (2014). A crucial aspect of the experiments, however, consists in the absence of translation invariance due to the presence of the external trapping potential. Consequently, correlations between points r and are not uniquely determined by the value of . This effect should be unimportant for short-range correlations, and, in fact, thermodynamic properties of ultracold gases are often very well-captured by a local density approximation. The situation changes drastically for a system at a critical point, where the correlation length diverges and thus competes with the inhomogeneity of the trapping potential.

Whereas the experimental preparation of critical systems typically requires a highly fine-tuned setup, they are rather effortlessly realized in two-dimensional (2D) systems whose low-energy excitations can be mapped onto an XY model. Examples apart from 2D ultracold quantum gases Hadzibabic and Dalibard (2011); Desbuquois et al. (2012); Ries et al. (2015); Fletcher et al. (2015); Murthy et al. (2015) are given by thin Helium films Bishop and Reppy (1978), layered magnets Dürr et al. (1989); Ballentine et al. (1990), and 2D exciton-polariton condensates Roumpos et al. (2012). To quantify correlations in these systems we introduce the one-body density matrix , where is the creation operator for a particle at point r. For the spatially homogeneous case we then have with some function due to translation and rotation invariance. The Mermin–Wagner–Hohenberg theorem Mermin and Wagner (1966); Hohenberg (1967) forbids long-range order at any finite temperature such that . However, in the XY model an ordered low-temperature phase with infinite correlation length exists, where correlations decay with a temperature-dependent scaling exponent as for large . Since is well below unity, this behavior of correlations is named quasi-long-range order (QLRO). Above a transition temperature, vortices disorder the system and QLRO is lost Holzmann and Krauth (2008); Choi et al. (2013). The described picture has been developed by Berezinskii and Kosterlitz, Thouless (BKT) Berezinskii (1971, 1972); Kosterlitz and Thouless (1973); Kosterlitz (1974).

The presence of the trap raises several conceptual questions towards the validity of the BKT picture in ultracold quantum gas experiments: On which length scales is QLRO visible? If samples bigger than the state of Texas are needed to reach the thermodynamic limit in the homogeneous system Bramwell and Holdsworth (1994), do we need ultracold atomic clouds covering the whole continent for the local density approximation to hold? Does the formula for a Bose gas with mass , temperature , and superfluid density imply a locally varying scaling exponent due to the inhomogeneous density? Can the 2D XY-model explain the large scaling exponent recently observed in experiment and Quantum Monte Carlo (QMC) computations Murthy et al. (2015)? Can we obtain the correlation function of the trapped system from the homogeneous one by a generalization of the local density approximation?

In this Rapid Communication, we address and answer all of these questions within the spin wave approximation for the phase fluctuations of the trapped gas. The main results are (1) the generalization of the textbook spin wave theory to a discrete collective mode spectrum, (2) the derivation of explicit analytical expressions for the one-body density matrix and first-order correlation functions that are suited for practical applications, and (3) the demonstration that the trap-averaged correlation function decays algebraically with an increased exponent. Our calculations provide a simple picture to the large exponents observed in experiment and simulation. We limit the derivation to the key steps and invite the mathematically interested reader to consult the detailed Supplemental Material (SM) SOM ().

Spin wave theory. At sufficiently high temperatures, quantum corrections are small and the macroscopic properties of the trapped gas are described by the classical action

 S[Φ]=1T∫d2r[ |∇Φ|22M+(−μ+V(r))|Φ|2+g2|Φ|4] (1)

for the complex bosonic field . Herein, and are the chemical potential and coupling constant, and we use units such that . We assume an isotropic and harmonic trapping potential . The action is stationary towards small changes in if is satisfied. We neglect gradient terms of the fields and obtain a Thomas–Fermi (TF) density profile given by with and TF radius . By neglecting deviations of from the stationary configuration we then arrive at the action for the phase given by

 Sph[θ]=n02MT∫Rλd2r (1−r2R2)(∇θ)2. (2)

We explicate the length scales that delimit the validity of our description, given by the TF radius at large distances, and a microscopic scale , which is given by the thermal wavelength in cold atom experiments, or by the lattice spacing in spin models. For typical 2D ultracold quantum gases we have and , such that is large.

After partial integration the phase-only action can be written as quadratic form with self-adjoint operator

 DR=−(1−r2R2)∇2+2r⋅∇R2 (3)

and scalar product . The normalized eigenfunctions of for are given by

 θnm(r)=√2n+|m|+1πs|m|/2P(|m|,0)n(1−2s)eimϕ (4)

with , , , , and Jacobi and Legendre polynomials and , respectively Abramowitz and Stegun (1964); SOM (). The energy levels read

 εnm=2|m|+4n(n+|m|+1). (5)

These energies match collective superfluid modes of harmonically trapped 2D Bose gases Boettcher et al. (2011)1. In fact, spin waves are related to the superfluid velocity by means of , i.e., they describe the same physical excitations Edwards et al. (1996); Stringari (1996); De Rossi et al. (2016).

Since the phase-only action is quadratic, though with a nontrivial discrete spectrum, correlation functions are obtained via Gaussian integration. Within the spin wave approximation we have with and

 ⟨Δθ(r,r′)2⟩=MTn0∑n,m ′ |θnm(r)−θnm(r′)|2εnm. (6)

In the homogeneous limit () we have and the eigenfunctions are plane waves with energies , and

 ⟨Δθ(r,r′)2⟩hom =MTn0∫λ−1d2q(2π)2|θq(r)−θq(r′)|2εq ≃η0log(|r−r′|2λ2) for |r−r′|≫λ, (7)

with Herbut (2007); Altland and Simons (2010). Due to and for large and small , we recover this formula for . The homogeneous case is recalled in the SM SOM (). In the homogeneous system, spin wave theory breaks down for Kosterlitz and Thouless (1973); Kosterlitz (1974). In a trap it is applicable to the superfluid core in the center of the atomic cloud. We estimate its radius from as . For () we have (), which is sufficiently large to observe the effects studied here. Much of the elegance of BKT theory stems from replacing the integral in Eq. (7) by the logarithm, which is valid for all relevant separations, but not exact.

Evaluation of the sum. Just as the momentum integration in the homogeneous case is limited by , also Eq. (6) has an ultraviolet cutoff on possible values of , indicated by the prime. In fact, besides , the values of and are limited from above according to . In the following we choose summation cutoffs with Euler’s constant SOM (). We can then write

 ⟨Δθ(r,r′)2⟩ =η02[F(N)0(s,s)−2F(N)0(s,s′)+F(N)0(s′,s′)] +η0M∑m=1[F(Nm)m(s,s)−2cos(mΔϕ)F(Nm)m(s,s′) +F(Nm)m(s′,s′)] (8)

with and

 F(N)0(s,s′) =N∑n=12n+1n(n+1)Pn(1−2s)Pn(1−2s′), F(Nm)m(s,s′) =Nm∑n=02n+m+1m2+n(n+m+1)(ss′)m/2 ×P(m,0)n(1−2s)P(m,0)n(1−2s′). (9)

Here and in the following we denote

 r=(rcosϕrsinϕ), s=r2R2, Δϕ=ϕ−ϕ′, s>=max(s,s′), s<=min(s,s′), ^λ=λR. (10)

Equation (8) can be implemented numerically and thus allows the computation of the phase fluctuations in a trap for arbitrary values of . Note in particular that . If is large, however, analytical approximation formulas can be derived from the limits.

Evaluating Eq. (9) for relies on the observation that for a self-adjoint operator with eigenvalues and eigenfunctions , the function is the Green function of the operator due to the completeness of the . Accordingly, if we know the Green function by some other means, we also know the result of the summation Englis and Peetre (1998); Gustavsson (2001). Now, for fixed the operator is of Sturm–Liouville-type and thus its Green function is given by , where are zero modes of the operator. With this approach we find

 F(∞)0(s,s′) =−1−log[s>(1−s<)], (11) F(∞)m(s,s′) =1m(s)m/2um(s<)vm(s>) (12)

for . We have and with hypergeometric function and , such that and Abramowitz and Stegun (1964); SOM (). The key properties of the functions and are and

 limm→∞um(s)=limm→∞vm(s)=1√1−s. (13)

We define the remainder function

 R(s,s′,Δϕ) =M0∑m=11m(s)m/2cos(mΔϕ) (14) ×(um(s<)vm(s>)−1√(1−s)(1−s′)),

which converges quickly in and is well-approximated by summing up the first terms. The careful limit in Eq. (8) is presented in the SM SOM ().

Correlation functions. The phase fluctuations due to Eq. (8) for are given by

 ⟨Δθ(r,r′)2⟩ ≃η02(1+11−s>)log(¯s>^λ2(1−s>))−η02(1−11−s<)log(¯s<^λ2(1−s<)) +η0√(1−s)(1−s′)log(Δ¯r2¯s>)+η0(R(s,s,0)+R(s′,s′,0)−2R(s,s′,Δϕ)) (15)

with . The latter two definitions give a proper way to treat the singularities that occur when . Their origin is already visible in the homogeneous limit (7): Whereas the integral is well-defined for all , the logarithm expression requires a meaningful procedure of how to treat the ultraviolet singularities. The remainder term involving the -functions in Eq. (15) can be neglected for most cases of interest. The homogeneous limit is recovered for , i.e., .

The central correlation function derived from Eq. (15) is given by

 g1(r)=n0(1−s)12[1+η(s)](rλ)−η(s)e−η02R(s,s,0) (16)

with . The factor is close to unity and may be neglected. The algebraic decay with for small agrees with the central correlation functions obtained in Refs. Petrov et al. (2000); Crecchi and Vicari (2011); Ceccarelli et al. (2013). The angular correlations for read

 ⟨Δθ(r,r,Δϕ)2⟩ =η01−slog(2s(1−cosΔϕ)^λ2(1−s)) +2η0[R(s,s,0)−R(s,s,Δϕ)] (17)

for . (This corresponds to . For the correlations vanish.) In Fig. 1 we plot a selection of correlation functions.

Applications. To obtain a useful estimate of the trapped correlation function from the homogeneous result, we write Eq. (7) as

 ⟨Δθ(r,r′)2⟩hom=2G(r,r′)−G(r,r)−G(%r′,r′) (18)

with , the Green function of . In each separately we apply a local density approximation and obtain

 ⟨Δθ(r,r′)2⟩LCA =η02[11−s+11−s′]log(1^λ2) +η0√(1−s)(1−s′)log(|r−r′|2R2). (19)

We name this particular procedure local correlation approximation (LCA), and observe from Fig. 1 that it works sufficiently well to approximate the actual result. In contrast, the approximation fails considerably.

The Fourier transform of the momentum distribution defines the trap-averaged correlation function Cohen-Tannoudji and Guéry-Odelin (2011); Murthy et al. (2015)

 ¯g1(r)=∫d2x ρ(x,r+x). (20)

In a translation invariant setting we have , whereas both functions differ in the trapped case. As shown in Figs. 1 and 2, decays algebraically for with an increased exponent -, which weakly depends on . In the SM SOM () we analytically estimate , which supports the exceptionally weak -dependence. Although for infinitely large systems, it would require absurd values of .

The spin wave description only applies to the superfluid core of the cloud. Hence, applying Eq. (15) to the whole cloud becomes less accurate for increasing values of and underestimates . To estimate the effect of the normal gas on the trap-averaging we set for with approximate core radius SOM (). Assuming the correlation length in the normal gas to coincide with should give the right order of magnitude of the effect. Computing with this modification of , we find an increased scaling exponent which is slightly above the spin wave prediction for . Increasing further, we obtain exponents up to for small , whereas the algebraic decay disappears for larger in the typical fitting range . Since most experiments are in the first regime, an enhanced scaling exponent is likely to be measured. Neglecting the spin wave part and only keeping exponential correlations, no algebraic decay of is observed.

In order to relate to the enhanced exponents found in Ref. Murthy et al. (2015), we compare to the QMC results where is determined by the decay of the central correlation function Holzmann et al. (2010)2. As shown in Fig 2, we find good qualitative agreement with the QMC results when including the exponential decay of correlations in the normal component. Hence we identify the superfluid core, a significant normal component, and a ratio of as a possible explanation for the large universal transition exponent . Although it is possible to include more quantitative knowledge of the quasi-2D equation of state Holzmann et al. (2008, 2010); Prokof’ev and Svistunov (2002); Boettcher et al. (2016); Wu et al. (2015) into the LCA approximation, it is clear that our qualitative findings are robust against these refinements. A full quantitative analysis of the BKT transition in the trapped gas is left for future work.

In conclusion, we have generalized the spin wave theory of phase fluctuations towards the trapped 2D Bose gas, analytically computed correlation functions, and applied this to obtain the trap-averaged correlation function. Our analysis provides a natural explanation for the unexpected large exponents observed in Ref. Murthy et al. (2015). It should be possible to apply the method to other dimensions and polytropic equations of state, just as it is the case for superfluid hydrodynamic modes Fuchs et al. (2003); Heiselberg (2004); Boettcher et al. (2011). We proposed an LCA prescription which can be further developed by considering these other mode spectra. An improved understanding of correlations in inhomogeneous systems is also mandatory for conducting experiments on universality in disordered systems Allard et al. (2012); Beeler et al. (2012); Carleo et al. (2013) and situations far from equilibrium Mathey et al. (2011); Altman et al. (2015); Dagvadorj et al. (2015).

Acknowledgments. We gratefully acknowledge inspiring discussions with A. Böttcher, R. Ganesh, I. F. Herbut, P. A. Murthy, and M. M. Scherer. IB acknowledges funding by the European Research Council (ERC) under Advanced Grant No. 290623 and the German Research Foundation (DFG) under Grant No. BO 4640/1-1. MH was supported by ANR-12-BS04-0022-01, ANR-13-JS01-0005-01 and DIM Nano’K from Région Île-de-France.

Supplemental Material

## Appendix A Computation of the correlation function

In this section we derive the formula for the phase fluctuations in detail. The presentation provides all the intermediate steps that have been left out in the main text for reasons of brevity.

We start by solving the eigenvalue problem for the operator

 DR=−(1−r2R2)∇2+2r⋅∇R2 (S1)

appearing in the phase-only action

 Sph[θ]=n02MT∫Rλd2r (1−r2R2)(∇θ)2=n02MT(θ,DRθ). (S2)

We introduce dimensionless variables and to obtain the eigenvalue problem

 ^Dθε=εθε. (S3)

We employ a polar coordinate Ansatz to arrive at

 −(1−^r2)(d2d^r2+1^rdd^r−m2^r2)θεm+2^rdd^rθεm=εθεm. (S4)

Note that since is invariant under , the eigenvalues only depends on . With the new coordinate

 s=^r2=r2/R2 (S5)

we eventually arrive at with

 ^D(m)=−dds[4s(1−s)dds]+m21−ss. (S6)

Note that this operator is of Sturm–Liouville form, see Eq. (S114). The eigenfunctions of are easily seen to be of the form with being the Jacobi polynomials, see Eq. (S98). The corresponding eigenvalues are given by

 εnm=2|m|+4n(n+|m|+1) (S7)

with such that . For the solution reduces to the Legendre polynomial . In App. D.1 we collect relevant properties of the Legendre and Jacobi polynomials.

We define the normalized eigenfunctions

 θnm(r)=√2n+|m|+1πs|m|/2P(|m|,0)|n(1−2s)eimϕ (S8)

such that

 ^Dθnm=R2DRθnm=εnmθnm, (S9)

and, due to Eq. (S95), . Due to the completeness of the Jacobi polynomials, every regular phase configuration can be expanded according to

 θ(r)=∑n,manmθnm(r). (S10)

Since is real we have . We arrive at the regular phase-only action

 Missing dimension or its units for \hskip (S11)

The smallest and largest length scales and introduced in Eq. (S2) manifest in lower and upper boundaries on the summation (indicated by a prime). The implications of this are addressed below in detail.

We now consider the one-body density matrix

 ρ(r,r′)=⟨Φ∗(r)Φ(r′)⟩. (S12)

Note that since this correlation function coincides with the connected correlation function . We write and assume to be regular. Approximating by the stationary solution we arrive at

 ρ(r,r′) =√n(r)n(r′)⟨e−i(θ(r)−θ(r′))⟩ (S13)

with

 ⟨O⟩=∫Dθ O e−Sph[θ]∫Dθ e−Sph[θ]. (S14)

Since is quadratic and we have

 ⟨e−i(θ(r)−θ(r′))⟩=e−12⟨[θ(r)−θ(r′)]2⟩ (S15)

according to the rules of Gaussian integration.

We proceed by computing

 ⟨Δθ(r,r′)2⟩=⟨[θ(r)−θ(r′)]2⟩ (S16)

with the Gaussian action . We write the functional measure as and find

 ⟨[θ(r)−θ(r′)]2⟩ =∫Dθ (θ(r)−θ(%r′))(θ(r)−θ(r′))e−Sph[θ]∫Dθ e−Sph[θ] =∑n,m∑n′,m′(θnm(r% )−θnm(r′))(θn′m′(r)−θn′m′(r′)) ×∏ν,μ∫daνμ anman′m′e−12γνμaνμaν,−μ∫daνμ e−12γνμaνμaν,−μ (S17)

with . In the product, all terms except one yield unity. The nontrivial one is due to

 ∫da a2e−12aγa∫da e−12aγa=1γ (S18)

for a real variable . This eventually yields

 ⟨Δθ(r,r′)2⟩=MTn0∑n,m ′ |θnm(r)−θnm(r′)|2εnm. (S19)

The macroscopic description of the trapped Bose gas in terms of the phase-only action is restricted to length scales . This implies an ultraviolet energy cutoff according to . The energies are thus restricted according to

 0<εnm<εmax∼(Rλ)2 (S20)

with . For the largest term in the -sum is

 N=N0=−1+√1+εmax2≃√εmax2∼(Rλ). (S21)

The largest value of , denoted by , is attained for and is found to be

 M=εmax2∼(Rλ)2. (S22)

If , the possible values of are limited from above by . Given we have

 2m+4Nm(Nm+m+1)<εmax (S23)

for the largest value of . Thus

 Nm =−m+12+√εmax+m2+12 =−m+12+√2M+m2+12. (S24)

In the following we choose

 N=e−γ(Rλ), M=e−γ(Rλ)2 (S25)

with Euler’s constant . For large and this implies

 2HN=log(1^λ2), HM=log(1^λ2) (S26)

with harmonic number . We will see below that fixing and according to Eqs. (S26) leads to convenient expressions for the phase fluctuations by trading for . Note, however, that there is no such that Eqs. (S25) can be true at the same time. On the other hand, due to , Eqs. (S25) are close approximations to , i.e., with a proportionality factor of unity in .

We now write the summation formula like in the main text as

 ⟨Δθ(r,r′)2⟩ =η02[F(N)0(s,s)−2F(N)0(s,s′)+F(N)0(s′,s′)] +η0M∑m=1[F(Nm)m(s,s)−2cos(mΔϕ)F(Nm)