A Computation of the correlation function

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


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.

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


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


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


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


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


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


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


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


with and


Here and in the following we denote


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


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


We define the remainder function


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


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., .

Figure 1: Selection of correlation functions for and in terms of . Black data represents the sum from Eq. (6) (points) and the analytic result from Eq. (15) (solid line). The LCA from Eq. (19) is indicated by the dashed orange lines. Upper panel. Central correlation function (upper curve) decaying algebraically with exponent , and trap-averaged correlation function (lower curve) decaying with strictly larger exponent . Long-dashed red lines correspond to with and , respectively. The inset shows the same data on a linear scale. Lower panel. Radial correlation function for fixed (from left to right). Note that the peaks heights are . The inset shows the angular correlation function for the same values of (from top to bottom). The long-dashed blue lines show , which deviates considerably.

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


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


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

Figure 2: Scaling exponent extracted from for (blue, dark red, orange). The solid lines give the result of applying Eq. (15) to the whole cloud. As a first estimate of the influence of the superfluid transition we replace the phase correlations in the outer regions of the cloud by an exponential decay with correlation length . This yields an increase of for small (dashed lines). The red triangles correspond to the QMC data from Ref. Murthy et al. (2015) for the quasi-2D gas with and , with determined from the scaling of . The error bars estimate the fitting error of , which is substantial due to its small value.

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


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


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)


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


appearing in the phase-only action


We introduce dimensionless variables and to obtain the eigenvalue problem


We employ a polar coordinate Ansatz to arrive at


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


we eventually arrive at with


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


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


such that


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


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


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


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




Since is quadratic and we have


according to the rules of Gaussian integration.

We proceed by computing


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


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


for a real variable . This eventually yields


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


with . For the largest term in the -sum is


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


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


for the largest value of . Thus


In the following we choose


with Euler’s constant . For large and this implies


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