Vortex pairing in two-dimensional Bose gases

Vortex pairing in two-dimensional Bose gases

Christopher J. Foster The University of Queensland, School of Mathematics and Physics, ARC Centre of Excellence for Quantum-Atom Optics, Qld 4072, Australia    P. Blair Blakie The Jack Dodd Centre for Quantum Technology, Department of Physics, University of Otago, Dunedin, New Zealand.    Matthew J. Davis The University of Queensland, School of Mathematics and Physics, ARC Centre of Excellence for Quantum-Atom Optics, Qld 4072, Australia

Recent experiments on ultracold Bose gases in two dimensions have provided evidence for the existence of the Berezinskii-Kosterlitz-Thouless (BKT) phase via analysis of the interference between two independent systems. In this work we study the two-dimensional quantum degenerate Bose gas at finite temperature using the projected Gross-Pitaevskii equation classical field method. While this describes the highly occupied modes of the gas below a momentum cutoff, we have developed a method to incorporate the higher momentum states in our model. We concentrate on finite-sized homogeneous systems in order to simplify the analysis of the vortex pairing. We determine the dependence of the condensate fraction on temperature and compare this to the calculated superfluid fraction. By measuring the first order correlation function we determine the boundary of the Bose-Einstein condensate and BKT phases, and find it is consistent with the superfluid fraction decreasing to zero. We reveal the characteristic unbinding of vortex pairs above the BKT transition via a coarse-graining procedure. Finally, we model the procedure used in experiments to infer system correlations [Hadzibabic et al., Nature 441, 1118 (2006)], and quantify its level of agreement with directly calculated in situ correlation functions.

03.75.Hh, 03.75.Lm

I Introduction

At low temperatures a three-dimensional (3D) Bose gas can undergo a phase transition to a Bose-Einstein condensate. In contrast, thermal fluctuations prevent a two-dimensional (2D) Bose gas from making a phase transition to an ordered state, in accordance with the Mermin-Wagner-Hohenberg theorem Mermin and Wagner (1966); Hohenberg (1967). However, the 2D Bose gas supports topological defects in the form of vortices, and in the presence of interactions can instead undergo a Berezinskii-Kosterlitz-Thouless (BKT) Berezinskii (1971); Kosterlitz and Thouless (1973); Posazhennikova (2006) transition to a quasi-coherent superfluid state. The BKT transition was first observed in liquid helium thin films Bishop and Reppy (1978), however, more recently, evidence for this transition has been found in dilute Bose gases Stock et al. (2005); Hadzibabic et al. (2006); Krüger et al. (2007); Schweikhard et al. (2007); Cladé et al. (2009) (also see Görlitz et al. (2001); Smith et al. (2005)).

Ultracold gases have proven to be beautiful systems for making direct comparisons between experiment and ab initio theory. Experiments in the 2D regime present a new challenge for theory as strong fluctuations invalidate mean-field theories (e.g. see Prokof’ev et al. (2001); Prokof’ev and Svistunov (2002); Posazhennikova (2006); Bisset and Blakie (2009); Mathey and Polkovnikov (2009); Sato et al. (2009); Gies and Hutchinson (2004); Schumayer and Hutchinson (2007)), and only recently have quantum Monte Carlo Holzmann and Krauth (2008); Holzmann et al. (2009) and classical field (c-field) Simula and Blakie (2006); Simula et al. (2008); Bisset et al. (2009) methods been developed that are directly applicable to the experimental regime.

In this paper we study a uniform Bose gas of finite spatial extent and parameters corresponding to current experiments. To analyze this system we use the projected Gross-Pitaevskii equation (PGPE), a c-field technique suited to studying finite temperature Bose fields with many highly occupied modes. We develop a technique for extracting the superfluid density based on linear response properties, and use this to understand the relationship between superfluidity and condensation in the finite system.

With this formalism we then examine two important applications: First, we provide a quantitative validation of the interference technique used in the ENS experiment to determine the nature of two-point correlation in the system. To do this we simulate the interference pattern generated by allowing two independent 2D systems to expand and interfere. Then applying the experimental fitting procedure to analyze the interference pattern we can extract the inferred two point correlations, which we can then compare against the in situ correlations that we calculate directly. Second, we examine the correlations between vortices and antivortices in the system to directly quantify the emergence of vortex-antivortex pairing in the low temperature phase. A similar study was made by Giorgetti et al. using a semiclassical field technique Giorgetti et al. (2007). We find results for vortex number and vortex pair distributions consistent with their results, and we show how a coarse graining procedure can be used to reveal the unpaired vortices in the system.

We now briefly outline the structure of this paper: In Sec. II we review the 2D Bose gas and relevant BKT physics. In Sec. III we outline the c-field technique and how it is specialized to describing a uniform, but finite, 2D Bose gas. In Sec. IV we present our main results, before concluding.

Ii Formalism

Here we consider a dilute 2D Bose gas described by the Hamiltonian


where is the atomic mass, , and is the bosonic field operator.

We take the two-dimensional geometry to be realized by tight confinement in the direction that restricts atomic occupation to the lowest mode. The dimensionless 2D coupling constant is


with the spatial extent of the mode 111For example, for tight harmonic confinement of frequency we have . and the s-wave scattering length. We will assume that so that the scattering is approximately three-dimensional Petrov et al. (2000), a condition well-satisfied in the ENS and NIST experiments Stock et al. (2005); Hadzibabic et al. (2006); Krüger et al. (2007); Cladé et al. (2009). For reference, the ENS experiment reported in Hadzibabic et al. (2006) had , whereas in the NIST experiments Cladé et al. (2009).

In contrast to experiments we focus here on the uniform case; no trapping potential in the plane is considered. We perform finite sized calculations corresponding to a square system of size with periodic boundary conditions. Working in the finite size regime simplifies the simulations and is more representative of current experiments. We note that the thermodynamic limit corresponds to taking while keeping the density, , constant.

ii.1 Review of BKT physics

The BKT superfluid phase has several distinctive characteristics, which we briefly review.

ii.1.1 First order correlations

Below the BKT transition the first-order correlations decay according to an inverse power law:


Systems displaying such algebraic decay are said to exhibit quasi-long-range order Chaikin and Lubensky (1995). This is in contrast to both the high temperature (disordered phase) in which the correlations decay exponentially, and long-range ordered case of the 3D Bose gas in which for .

ii.1.2 Superfluid density

Nelson and Kosterlitz Nelson and Kosterlitz (1977) found that the exponent of the algebraic decay is related to the ratio of the superfluid density and temperature. To within logarithmic corrections


where is the superfluid density and is the thermal de Broglie wavelength. Furthermore, Nelson and Kosterlitz showed that this ratio converges to a universal constant as the transition temperature, , is approached from below: (i.e., ). Thus, the superfluid fraction undergoes a universal jump from to as the temperature decreases through .

ii.1.3 Vortex binding transition

Another important indicator of the BKT transition is the behavior of topological excitations, which are quantized vortices and antivortices in the case of a Bose gas. A single vortex has energy which scales with the logarithm of the system size. At low temperatures this means that the free energy for a single vortex is infinite (in the thermodynamic limit), and vortices cannot exist in isolation. As originally argued in Kosterlitz and Thouless (1973), the entropic contribution to the free energy also scales logarithmically with the system size, and will dominate the free energy at high temperatures allowing unbound vortices to proliferate. This argument provides a simple estimate for the BKT transition temperature.

Although unbound vortices are thermodynamically unfavored at , bound pairs of counter-rotating vortices may exist since the total energy of such a pair is finite 222The vortex-antivortex pair energy depend on the pair size rather than the system size.. This leads to a distinctive qualitative characterization of the BKT transition: as the temperature increases through pairs of vortices unbind.

ii.1.4 Location of the BKT transition in the dilute Bose gas

While the relation between the superfluid density and temperature at the transition is universal, the total density, , at the transition is not. General arguments Popov (1983); Kagan et al. (1987); Fisher and Hohenberg (1988) suggest that the transition point for the dilute uniform 2D Bose gas is given by


where is a constant. Prokofév, Ruebenacker and Svistunov Prokof’ev et al. (2001); Prokof’ev and Svistunov (2002) studied the homogeneous Bose gas using Monte Carlo simulations of an equivalent classical model on a lattice. Using an extrapolation to the infinite-sized system, they computed a value for the dimensionless constant, . By inverting Eq. (5), we obtain the BKT critical temperature for the infinite system


We use the superscript to indicate that this result holds in the thermodynamic limit.

Iii Method

iii.1 c-field and incoherent regions

We briefly outline the PGPE formalism, which is developed in detail in Ref. Blakie et al. (2008). The Bose field operator is split into two parts according to


where is the coherent region c-field and is the incoherent field operator (see Blakie et al. (2008)). These fields are defined as the low and high energy projections of the full quantum field operator, separated by the cutoff wave vector . In our theory this cutoff is implemented in terms of the plane wave eigenstates of the time-independent single particle Hamiltonian, that is,


with . The fields are thus defined by


where the are Bose annihilation operators, the are complex amplitudes, and the sets of quantum numbers defining the regions are

C (12)
I (13)

iii.1.1 Choice of C region

In general, the applicability of the PGPE approach to describing the finite temperature gas relies on an appropriate choice for , so that the modes at the cutoff have an average occupation of order unity. In this work we choose an average of five or more atoms per mode using a procedure discussed in appendix A. This choice means that all the modes in C are appreciably occupied, justifying the classical field replacement . In contrast the I region contains many sparsely occupied modes that are particle-like and would be poorly described using a classical field approximation. Because our 2D system is critical over a wide temperature range, additional care is needed in choosing C. Typically strong fluctuations occur in the infrared modes up to the energy scale . Above this energy scale the modes are well described by mean-field theory (e.g. see the discussion in Kashurnikov et al. (2001); Prokof’ev et al. (2001)). For the results we present here, we have


for simulations around the transition region and at high temperature. At temperatures well below , the requirement of large modal occupation near the cutoff competes with this condition and we favor the former at the expense of violating Eq. (14).

iii.1.2 PGPE treatment of C region

The equation of motion for is the PGPE


where the projection operator


formalizes our basis set restriction of to the C region. The main approximation used to arrive at the PGPE is to neglect dynamical couplings to the incoherent region Davis et al. (2001a).

We assume that Eq. (15) is ergodic Davis et al. (2001b), so that the microstates {} generated through time evolution form an unbiased sample of the equilibrium microstates. Time averaging can then be used to obtain macroscopic equilibrium properties. We generate the time evolution by solving the PGPE with three adjustable parameters: (i) the cutoff wave vector, , that defines the division between C and I, and hence the number of modes in the C region; (ii) the number of C region atoms, ; (iii) the total energy of the C region, . The last two quantities, defined as


are important because they represent constants of motion of the PGPE (15), and thus control the thermodynamic equilibrium state of the system.

iii.1.3 Obtaining equilibrium properties for the C region

To characterize the equilibrium state in the C region it is necessary to determine the average density, temperature and chemical potential, which in turn allow us to characterize the I region (see Sec. III.2). These and other C region quantities can be computed by time-averaging, e.g, the average C region density is given by


where is a set of times (after the system has been allowed to relax to equilibrium) at which the field is sampled. We typically use 2000 samples from our simulation to perform such averages over a time of s. Another quantity of interest here is the first order correlation function, which we calculate directly via the expression


Derivatives of entropy, such as the temperature () and chemical potential () can be calculated by time averaging appropriate quantities constructed from the Hamiltonian (17) using the Rugh approach Rugh (1997). The detailed implementation of the Rugh formalism for the PGPE is rather technical and we refer the reader to Refs. Davis and Morgan (2003); Davis and Blakie (2005) for additional details of this procedure.

A major extension to the formalism of the PGPE made in this work is the development of a method for extracting the superfluid fraction, , from these calculations. For this we use linear response theory to relate the superfluid fraction to the long wavelength limit of the second order momentum density correlations. An extensive discussion of this approach, and the numerical methods used to implement it, are presented in appendix D.

iii.2 Mean-field treatment of I region

Occupation of the I region modes, , accounts for about 25% of the total number of atoms at temperatures near the phase transition. We assume a time independent state for the I region atoms defined by a Wigner function Naraschewski and Glauber (1999), allowing us to calculate quantities of interest by integrating over the above-cutoff momenta, Bezett et al. (2008); Davis and Blakie (2006).

Our assumed Wigner function corresponds to the self-consistent Hartree-Fock theory as applied in Davis and Blakie (2006). In two dimensions this is




is the Hartree-Fock energy, is the I region density, and is the chemical potential (shifted by the mean-field interaction with the I region atoms). Note that the average densities are constant in the uniform system, so has no explicit dependence, however, we include this variable for generality when defining the associated correlation function.

The I region density appearing in Eq. (22) is given by


with corresponding atom number ; total number is simply


An analytic expression for and simplified procedure for numerically calculating the first order correlation function of the I region atoms, , can be obtained by taking integrals over the phase space. These results are discussed in appendix B.

iii.3 Equilibrium configurations with fixed and

Generating equilibrium classical fields with given values of and is straightforward since the PGPE simulates a microcanonical system (see appendix A.3). However, we wish to simulate systems with a given temperature and total number. As described in the preceding two sections these can only be determined after a simulation has been performed. In appendix A we outline a procedure for estimating values of and for desired values of and based on a root finding scheme using a Hartree-Fock-Bogoliubov analysis for the initial guess.

Iv Results

We choose simulation parameters in analogy with the Paris experiment of Hadzibabic et al. Hadzibabic et al. (2006). This experiment used an elongated atomic cloud of approximately Rb atoms, with a spatial extent (Thomas-Fermi lengths) of 120 m and 10 m along the two loosely trapped and directions. The tight confinement in the direction was provided by an optical lattice.

Although our simulation is for a uniform system, we have chosen similar parameters where possible. Our primary simulations are for a system in a square box with m, with Rb atoms. We also present results for systems with m and m at the same density in order to better understand finite-size effects. All simulations are for the case of corresponding to the experimental parameters reported in Hadzibabic et al. (2006).

The cutoff wave vector varied with temperature to ensure appropriate occupation of the highest modes (see Sec. III.1.1). For the 100 m system, the number of C region modes ranged between 559 at low temperatures to 11338 at the highest temperature studied.

iv.1 Simulation of expanded interference patterns between two systems

In order to make a direct comparison with the experimental results of Hadzibabic et al. (2006), we have generated synthetic interference patterns and implemented the experimental analysis technique. Our simulated imaging geometry is identical to that found in Hadzibabic et al. (2006), with expansion occurring in the -direction. The interference pattern is formed in the - plane via integration of the density along the -direction (“absorption imaging”).

Our algorithm for obtaining the interference pattern due to our classical field is very similar to that presented in Hadzibabic et al. (2004). Our above cutoff thermal cloud is taken into account separately. We consider a pair of fields from different times during the simulation, chosen such that the fields can be considered independent. The 3D wavefunction corresponding to each field is reconstructed by assuming a harmonic oscillator ground state in the tight-trapping direction. These two reconstructed fields are spatially separated by m, corresponding to the period of the optical lattice in Hadzibabic et al. (2004).

Given this initial state, we neglect atomic interactions and only account for expansion in the tightly-trapped direction. This yields a simple analytical result for the full classical field at later times. The contribution of the above-cutoff atoms is included by an incoherent addition of intensities. The result is integrated along the -direction to simulate the effect of absorption imaging with a laser beam, that is,


Rather than integrate the full field along the -direction, we use only a slice of length m in keeping with the experimental geometry of Ref. Hadzibabic et al. (2006).

The interference patterns, , generated this way contained fine spatial detail not seen in the experimental images. To make a more useful comparison to experiment it is necessary to account for the finite optical imaging resolution by applying a Gaussian convolution in the - plane with standard deviation 3 m. Had ().

In accordance with the Paris experiment, we use a 22 ms expansion time to generate interference patterns for quantitative analysis (see Sec. IV.3.2). To obtain characteristic interference images for display in Hadzibabic et al. (2006), the experiments used a shorter 11 ms expansion Had (). We exhibit examples of interference patterns at various temperatures in Fig. 1, for this shorter expansion time. These images show a striking resemblance to the results presented in Ref. Hadzibabic et al. (2006).

Figure 1: Synthetic interference patterns generated from the 50 m grid by simulation of the experimental procedure of Ref. Hadzibabic et al. (2006). (a) At low temperatures, , the interference fringes are straight. (b) Just below the transition temperature, , the fringes become wavy due to decreased spatial phase coherence. Phase dislocations become common at temperatures above the transition, (c) , and (d) . These “zipper patterns” indicate the presence of free vortices. (e) When simulation of the finite imaging resolution is disabled, the zipper patterns from the field in sub figure (d) are no longer clearly visible; the high-frequency details obscure the phase information without providing obvious additional information about the existence of vortex pairs.

iv.2 Condensate and superfluid fractions

For a 2D Bose gas in a box we expect a nonzero condensate fraction due to the finite spacing of low-energy modes. A central question is whether we can observe a distinction between the crossover due to Bose condensation and that due to BKT physics. To address this question we have computed both the condensate and superfluid fractions from our dynamical simulations.

The condensate fraction in a homogeneous system is easily identified as the average fractional occupation of the lowest momentum mode. This is directly available from our simulations as a time average of the mode of the classical field,


Extracting the superfluid fraction from dynamical classical field simulations provides a more difficult challenge. For this we use linear response theory to relate the superfluid fraction to the long wavelength limit of the second order momentum density correlations. Details concerning the technique are given in appendix D.

Figure 2: (color online). Condensate fraction (solid dots) and superfluid fraction (crosses) as functions of temperature for the 100 m grid. The transition temperature in the thermodynamic limit, Prokof’ev et al. (2001), is shown as a vertical dot-dashed line. The vertical dashed line shows our estimate for the transition temperature in the finite system. The thick solid line is the condensate fraction for an ideal Bose gas in the grand canonical ensemble with the same number of atoms and periodic spatial domain. The superfluid fraction becomes negative in places because the extrapolation of the momentum correlations to is sensitive to statistical noise at high temperature (see appendix D.2 for details).
Figure 3: (color online). Detail of the superfluid fraction near the transition temperature. Solid dots represent the calculation based on momentum correlations as described in appendix D. Results for the largest and smallest grids are shown (left and right, respectively). The data for the 100 m grid is omitted for clarity, but lies between the curves shown as expected. Open circles represent the calculation of the superfluid fraction from the associated fitted values for the decay coefficient , via Eq. (4). The open circles terminate where the power law fitting procedure fails.

Figure 2 compares the results for the superfluid and condensate fractions computed on the 100 m grid. These results are qualitatively similar to the results for the larger and smaller grids. In particular, we note that there is no apparent separation between temperatures at which the superfluid and condensate fractions fall to zero. Also shown in Fig. 2 is the condensate fraction for the ideal Bose gas confined to an identical finite-size box in the grand canonical ensemble. The large shift between ideal and computed transition temperatures indicates the effect of interactions in the 2D system. Because the average system density is uniform, this large shift is to due to critical fluctuations (also see Kashurnikov et al. (2001)).

In our calculations we identify the transition temperature, , as where the superfluid fraction falls off most rapidly (i.e., the location of steepest slope on the versus graph; see Fig. 2). As the system size increases, this transition temperature moves toward the value for an infinite-sized system, Prokof’ev et al. (2001). This effect is illustrated by the behavior of the superfluid fraction in Fig. 3.

iv.3 First order correlations — algebraic decay

Algebraic decay of the first order correlations, as described by Eq. (3), is a characteristic feature of the BKT phase. Above the BKT transition, the first order correlations should revert to the exponential decay expected in a disordered phase.

The normalized first order correlation function, is defined by


where is the unnormalized first order correlation function Naraschewski and Glauber (1999).

iv.3.1 Direct calculation of

In the PGPE formalism the C and I contributions to the correlation function are additive Bezett et al. (2008), that is,


where and are defined in Eqs. (20) and (49), respectively. It is interesting to note that and individually display an oscillatory decay behavior — originating from the cutoff — an effect which correctly cancels when the two are added together.

Having calculated , we obtain the coefficient by fitting the algebraic decay law, Eq. (3), using nonlinear least squares; sample fits are shown in Fig. 4. The fit is conducted over the region between 10 and 40 de Broglie wavelengths. The short length scale cutoff is to avoid the contribution of the non-universal normal atoms, for which the thermal de Broglie wavelength sets the appropriate decay length. The long distance cutoff is chosen to be small compared to the length scale , to avoid the effect of periodic boundary conditions on the long range correlations.

The quality of the fitting procedure, and the breakdown of expression (3) at the BKT transition can be observed by adding an additional degree of freedom to the fitting function. In particular, at each temperature we fit the quadratic and extract the parameter ( is discarded). The abrupt failure of the fits can be observed in the inset of Fig. 5 as a sudden increase in the value of — an effect which is in excellent agreement with the value of as estimated from the superfluid fraction.

Figure 4: (color online). Sample fits to the algebraic decay of at various temperatures, ranging from below to above the transition. High temperatures correspond to curves at the bottom of the figure which have rapid falloff of with distance. Fits are shown on a log-log scale in the inset to emphasize the failure of a power law in describing the behavior of at high temperature.
Figure 5: (color online). Comparison of two methods for determining the algebraic decay coefficient for the first order correlation function . The line with circle markers represents direct fits to . These fits fail at the transition temperature as shown by the sharply diverging value of in the inset. The filled points represent the values obtained from a simulation of the experimental analysis procedure of Hadzibabic et al. (2006), described in Sec. IV.3.2. Horizontal dotted lines at 0.25 and 0.5 correspond to the expected values of just below and above the transition, respectively Hadzibabic et al. (2006). The vertical line is the BKT transition temperature, as estimated from the superfluid fraction calculated in Sec. IV.2.

iv.3.2 Calculation of via interference patterns

So far a direct probe of the in situ spatial correlations has not been possible, although important progress has been made by the NIST group Cladé et al. (2009). In the experiments of Hadzibabic et al. Hadzibabic et al. (2006) a scheme proposed by Polkovnikov et al. Polkovnikov et al. (2006) was used to infer these correlations from the “waviness” of interference patterns produced by pair of quasi-2D systems (see Sec. IV.1). In this section we simulate the experimental data analysis method, and compare inferred predictions for the correlation function against those we can directly calculate. This allows us to characterize the errors associated with this technique arising from finite size effects and finite expansion time.

To make this analysis we follow the procedure outlined in Hadzibabic et al. (2006). We fit our numerically generated interference patterns (see Sec. IV.1) to the function


where is a Gaussian envelope in the -direction, is the interference fringe contrast, is the fringe spacing and is the phase of the interference pattern in the -direction.

Defining the function


the nature of spatial correlations is then revealed by the manner in which decays with . In particular, we identify the parameter , defined by Polkovnikov et al. (2006). For an infinite 2D system in the superfluid regime () (i.e. corresponds to the algebraic decay of correlations). For , where correlations decay exponentially, is equal to .

Fitting to the algebraic decay law we can determine . A comparison between inferred from the interference pattern and obtained directly from is shown in Fig. 5. Both methods give broadly consistent predictions for when , however our results show that there is a clear quantitative difference between the two schemes, and that underestimates the coefficient of algebraic decay in the system (i.e. using in Eq. (4) would overestimate the superfluid density). Near and above transition temperature, where the fits to fail, we observe that converges toward . The agreement between and in the low temperature region improves as the size of the grid is increased.

iv.4 Vortices and pairing

The simplest description of the BKT transition is that it occurs as a result of vortex pair unbinding: At vortices only exist in pairs of opposite circulation, which unbind at the transition point to produce free vortices that destroy the superfluidity of the system. However, to date there are no direct experimental observations of this scenario, and theoretical studies of 2D Bose gases have been limited to qualitative inspection of the vortex distributions. In the c-field approach vortices and their dynamics are clearly revealed, unlike other ensemble-based simulation techniques where the vortices are obscured by averaging. This gives us a unique opportunity to investigate the role of vortices and pairing in a dilute Bose gas.

We detect vortices in the c-field microstates by analyzing the phase profile of the instantaneous field (see appendix C). An example of a phase profile of a field for is shown in Fig. 6(a). The vortex locations reveal a pairing character, that is, the close proximity of pairs of positive (clockwise) and negative (counterclockwise) vortices relative to the average vortex separation. An important qualitative feature of our observed vortex distributions is that at high temperatures, pairing does not disappear from the system entirely. Indeed, most vortices at high temperature could be considered paired or grouped in some manner, as shown in Fig. 6(b). Perhaps this is not surprising, since positive and negative vortices have a logarithmic attraction, and we observe them to create and annihilate readily in the c-field dynamics. However, this does indicate that the use of pairing to locate the transition may be ambiguous, and we examine this aspect further below.

Figure 6: (color online). Phase profile of a c-field with vortices indicated. Vortices with clockwise (white ) and anticlockwise (black ) circulation. The phase of the classical field is indicated by shading the background between dark blue (phase 0) and light yellow (phase ). (a) Distinctive pairing below the transition at (b) A “vortex plasma” above the transition at .

It is also of interest to measure the number of vortices, , present in the system as a function of temperature (see Fig. 7). At the lowest temperatures the system is in an ordered state, and the energetic cost of having a vortex is prohibitive. As the temperature increases there is a rapid growth of vortex population leading up to the transition point followed by linear growth above .

Figure 7: Total number of vortices (dots) and number of unpaired vortices (circles) as a function of temperature near the transition. While at the transition temperature is already very high, becomes nonzero only close to the transition, providing clear evidence of vortex unbinding at work. The inset shows the variation in the total number over the full temperature range of the simulations. Above the transition temperature the growth in the number of vortices becomes linear with temperature.

iv.4.1 Radial vortex density

The most obvious way to characterize vortex pairing is by defining a pair distribution function for vortices of opposite sign. Adopting the notation of Giorgetti et al. (2007), this is


where is the vortex density function which consists of a sum of delta-spikes,

for positive vortices at positions . We use the analogous definition for . The associated dimensionless two-vortex correlation function is


The angular average of can be calculated directly from the detected vortex positions using a binning procedure on the pairwise distances , and is shown in Fig. 8.

These results quantify the effect discussed earlier: Positive and negative vortices show a pairing correlation which does not disappear above . The characteristic size of this correlation, given by twice the width of the peak feature in Fig. 8, is m (taking full width half maximum).

The shape of our pairing peak is qualitatively similar to that described in Giorgetti et al. (2007). However, in contrast to their results the width does not appear to change appreciably with temperature. Additional simulations show that increasing the interaction strength causes the peak to become squarer and wider. It is clear that while the pair size and strength revealed in does not change appreciably as the transition is crossed, the amount of pairing relative to the background uncorrelated vortices changes considerably. This background of uncorrelated vortices is given by the horizontal plateau at large as shown in the inset.

Figure 8: (color online). Angular average of the two-vortex pair distribution functions for vortices of opposite sign. Three temperatures centered about the transition are shown: dot markers nK , ; circle markers nK , ; cross markers nK , . The vertical dotted line shows the value of the healing length at . The main plot shows normalized by the positive vortex density; comparable magnitudes for the peaks near show that vortex pairing remains important over the range of temperatures studied, not only below the transition. The inset shows in the natural dimensionless units for which as .

iv.4.2 Revealing unpaired vortices with coarse-graining

The function clearly indicates the existence of vortex pairing in the system. However, it does not provide a convenient way to locate the transition temperature, since a large amount of pairing exists both below and above the transition: The expected number of neighbors for any given vortex — roughly, the area of the pairing peak of shown in Fig. 8 — does not change dramatically across the transition. is the expected density of positive vortices.

We desire a quantitative observation of vortex unbinding at the transition and have therefore investigated several measures of vortex pairing 333For example, the Hausdorff distance (see, e.g., Papadopoulos (2005)) between the set of positive vortices and the set of negative vortices.. However, measures based directly on the full set of vortex positions seem to suffer from the proliferation of vortices at high temperature — an effect which tends to wash out clear signs of vortex unbinding. With this in mind, we have developed a procedure for measuring the number of unpaired vortices in our simulations, starting from the classical field rather than the full set of vortex positions.

The basis of our approach for detecting unpairing is to coarse-grain the classical field by convolution with a Gaussian filter of spatial width (standard deviation) . This removes all vortex pairs on length scales smaller than than . Figure 9 shows the count of remaining vortices as a function of filter width, along with some examples of coarse-grained fields. For , the number of remaining vortices levels off and only decreases slowly with increasing . Ultimately the number of remaining vortices goes to zero as .

Setting the filter width to be larger than the characteristic pairing distance, , yields a coarse-grained field from which the pairs have been removed, but unpaired vortices remain. In our simulations we have m; we take the vortices which remain after coarse-graining with a Gaussian of standard deviation m to give an estimate of the number of unpaired vortices, . Figure 7 shows that becomes nonzero only near the transition, in contrast to which is nonzero well below . The sharp increase in at is a quantitative demonstration of vortex unbinding at work.

Figure 9: (color online). The coarse-graining procedure: number of vortices as a function of filter width for a temperature near the transition. The smooth curve is an average over many realizations of the field, whereas the stepped curve shows typical behavior of the number for a single field. Insets show the coarse-grained fields for various filter widths; the transformation removes vortex-antivortex pairs which are separated by approximately less than the standard deviation of the filter. In this example unpaired vortices remain at m.

In the experiment of Ref. Hadzibabic et al. (2006), the fraction of interference patterns with dislocations (e.g., see Figs. 1(c) and (d)) was measured. While isolated vortices are clearly identified by interference pattern dislocations, a lack of spatial resolution in experiments means that this type of detection method obscures the observation of tightly bound vortex pairs. The experimental resolution of 3 m is broadly consistent with the scale of the coarse-graining filter (i.e., m). With this in mind, we introduce the quantity , defined as the probability of observing an unpaired vortex in a m control volume at a given temperature 444We choose a fixed control volume with m in order to compare results between simulations with different grid sizes.. For the 50 m grid we have simply .

Computing from our simulations yields the results shown in Fig. 10. Our results show a dramatic jump in at a temperature that is consistent with the transition temperature determined from the superfluid fraction calculation presented in Sec. IV.2.

Figure 10: (color online). Comparison of vortex unpairing measures. The dots are our pairing measure based on coarse-graining the field. Circles represent the pairing as determined by the number of dislocations in the simulated interference patterns. This was the same method used in the experimental analysis of Hadzibabic et al. (2006) and coincides remarkably well with our coarse-graining based measure. Both curves are consistent with the vertical line showing the transition temperature as determined from the superfluid fraction calculation in Sec. IV.2. The inset shows the calculated coarse-grained pairing measure for all three grid sizes, along with vertical lines showing the estimates for derived from the superfluid fraction calculations.

From the definition, we expect that should be close to the experimentally measured frequency of dislocations. To demonstrate this relationship, we have simulated interference patterns (as described in Sec. IV.1) and detected dislocations using the experimental procedure of Ref. Hadzibabic et al. (2006): A phase gradient was considered to mark a dislocation whenever radm. From this we can compute the probability of detecting at least one dislocation as a function of temperature. As shown in Fig. 10, the results of this procedure compare very favorably with our measure of pairing based on . We note that inhomogeneous effects in experiments probably broaden the jump in appreciably compared to our homogeneous results.

V Conclusion

In this paper we have used c-field simulations of a finite-sized homogeneous system in order to investigate the physics of the 2D Bose gas in a regime corresponding to current experiments. We have directly computed the condensate and superfluid fractions as a function of temperature, and made comparisons to the superfluid fraction inferred by the first order correlation function, and using the interference scheme used in experiments. Our results for these quantities provide a quantitative test of the interference scheme for a finite system.

A beautiful possibility is the direct experimental observation of vortex-antivortex pairs, their distribution in the system, and hence a quantitative measurement of their unbinding at the BKT transition. We have calculated the vortex correlation function across the transition and provided a coarse-graining scheme for distinguishing unpaired vortices. These results suggest that the dislocations observed in experiments, due to limited optical resolution, provide an accurate measure of the unpaired vortex population and accordingly are a strong indicator of the BKT transition.

We briefly discuss the effect that harmonic confinement (present in experiments) would have on our predictions. The spatial inhomogeneity will cause the superfluid transition to be gradual, occurring first at the trap center where the density is highest, in contrast to our results where the transition occurs in the bulk. So far the superfluid fraction for the trapped system has been determined by using the universality result for the critical density in the homogeneous gas Prokof’ev and Svistunov (2002) in combination with the local density approximation Holzmann and Krauth (2008); Bisset et al. (2009). It would be interesting to be able to compute the superfluid fraction independently as we have done here; however it is not clear how to do so.

Bisset et al. Bisset et al. (2009) used an extension of the c-field method for the trapped 2D gas to examine and found similar results for the onset of algebraic decay of correlations at the transition. Their analysis was restricted to the small region near the trap center where the density is approximately constant; we expect the results of our vortex correlation function and the coarse-graining scheme should similarly be applicable to the trapped system in the central region. Except in very weak traps, the size of this region is relatively small and will likely prove challenging to measure experimentally.

Our results for the homogeneous gas emphasize the clarity with which ab initio theoretical methods can calculate quantities directly observable in experiments, such as interference patterns. This should allow direct comparisons with experiments, providing stringent tests of many-body theory.

The authors are grateful for several useful discussions with Keith Burnett and Zoran Hadzibabic. CJF and MJD acknowledge financial support from the Australian Research Council Centre of Excellence for Quantum-Atom Optics. PBB is supported by FRST contract NERF-UOOX0703.

Appendix A Simulation using the PGPE

Here we outline our procedure for determining the properties of the C region and the steps used to create initial states for the PGPE solver. The C region itself is characterized by the cutoff momentum , while the initial states are characterized by the energy and number . We want to obtain values of these three properties which are consistent with a specified temperature and total number of atoms .

a.1 Hartree-Fock-Bogoliubov analysis

To generate an initial estimate of the C region parameters we solve the self-consistent Hartree-Fock-Bogoliubov (HFB) equations in the so-called Popov approximation Griffin (1996) to find an approximate thermal state for the system at a temperature . The resulting state is a Bose Einstein distribution of quasiparticles interacting only via the mean-field, expressed in terms of the quasiparticle amplitudes and .

Occupations for the C region field may be computed directly from the quasiparticle occupations via


where is the Bose Einstein distribution and is the quasiparticle energy which is obtained by solving the Bogoliubov-de Gennes equations self consistently  Griffin (1996). This allows us to compute the cutoff as the maximum value of consistent with sufficient modal occupation:


We choose for the sufficient occupation condition on the C region modes.

The number of atoms below the cutoff may be computed directly from the sum of the condensate number and the number of C region excited state atoms, :


For the total energy below the cutoff, we use the expression


where . Rearranging, this is


The expression (38) differs from Eq. (22) of Griffin (1996) as we have retained the zeroth order (constant) terms that are required to match the energy scale of the HFB analysis to the zero point of energy in the classical field simulations.

a.2 Initial conditions for fixed total number

A simple comparison between simulations at varying temperatures can only be carried out if the total number of atoms is fixed. This presents a problem in our simulations: although the number of atoms and energy of the C region can be directly specified (see Sec. A.3), we may only determine the total number after performing a simulation. This is because the number of atoms in the I region depends on the temperature and chemical potential which are calculated by ergodic averaging of the C region simulations.

Formally, this may be stated as a root-finding problem: solve


with initial guess provided by the solution to the HFB analysis in Sec. A.1. Although both and affect the total number , we choose to fix to the initial guess and to vary until the desired total number is found.

We note that evaluating the function is very computationally expensive and difficult to fully automate since it involves a simulation and several steps of analysis. For this reason we use a nonstandard root finding procedure: For the first iteration we simulate three energies about the initial guess such that the results crudely span ; these three simulations can be performed in parallel which significantly reduces the time to a solution. A second guess was obtained by quadratic fitting of as a function of which gives accurate to within about 5% of . An addition iteration using the same interpolation method takes to within 0.3%, which we consider sufficient.

We note that changing during the root finding procedure means we have no direct control over the final temperature of each specific simulation. In our case this is not a problem because we only require a range of temperatures spanning the transition. In principle one could solve for a given temperature by allowing to vary in addition to .

a.3 Initial conditions for given and

We compute initial conditions for the C region field in a similar way to Davis et al. (2002). Using the representation for the C region given by Eq. (10), the task is to choose appropriate values for the . As a first approximation, choose the smallest value for a momentum cutoff such that the field with coefficients


has energy greater than . Here is chosen so that the field has normalization corresponding to atoms, and is a randomly chosen phase which is fixed for each mode at the start of the procedure. The random phases allow us to generate many unique random initial states at the same energy.

By definition, the field defined by (40) has energy slightly above the desired energy. This problem is solved by mixing it with the lowest energy state:


using a root finding procedure to converge on the desired energy . The scheme generates random realizations of a non-equilibrium field with given and which are then simulated to equilibrium before using ergodic averaging for computing statistics.

Appendix B I region integrals

Our assumed self-consistent Wigner function (Sec. III.2) for the I region atoms takes a particularly simple form in the homogeneous case:


The above-cutoff density may then be found by direct integration:


with the thermal de Broglie wavelength.

In a similar way, the assumed Wigner function allows any desired physical quantity to be estimated via a suitable integral. A particular quantity of interest in the current work is the first order correlation function, which can be obtained from the Wigner function as Naraschewski and Glauber (1999)


This integral is of the general form


for constants and . Noting that depends only on the length, of , and transforming to polar coordinates , we have


(see (Gradshteyn and Ryzhik, 2000, p902) for the Bessel function identity).

Thus we obtain in terms of a one dimensional integral which may be performed numerically:


Appendix C Vortex detection

The defining feature of a “charge-” vortex is that the phase of the complex field changes continuously from to around any closed curve which circles the vortex core. We express our field on a discrete grid in position space; the aim of vortex detection is then to determine which grid plaquettes (that is, sets of four adjacent grid points) contain vortex cores.

To obtain the phase winding about a plaquette, first consider the phase at two neighboring grid points A and B. We are interested in the unwrapped phase difference between the grid points; unwrapping ensures that the phase is continuous between A and B. (In the discrete setting such continuity is poorly defined; the best we can do is to correct for the possibility of phase jumps by adding or subtracting factors of so that .) The unwrapped phase differences around a grid plaquette tell us a total phase change where is the winding number or “topological charge”.

Due to the necessity of unwrapping the phase, a four-point grid plaquette cannot unambiguously support vortices with charge larger than one. Luckily, such vortices are energetically unfavorable in 2D Bose gases (Pitaevskii and Stringari, 2003, p.83) so we need only concern ourselves with detecting vortices with winding number in this work. The positions obtained from a given run of our vortex detection algorithm are the labeled and for winding numbers and , respectively.

Appendix D Superfluid fraction

One of the important characteristics of the BKT transition is the presence of superfluidity, even in the absence of conventional long range order. In the following we describe a method to calculate the superfluid fraction from our classical field description. The method is attractive because it makes use of momentum correlations which may be extracted directly from our equilibrium simulations without any need to introduce additional boundary conditions or moving defects.

d.1 Superfluid fraction via momentum density correlations

Our derivation is based on the procedure presented in (Forster, 1975, p.214), (see also Baym (1969) and (Pitaevskii and Stringari, 2003, p.96)). The central idea is to establish a relationship between i) the autocorrelations of the momentum density in the simulated ensemble and ii) the linear response of the fluid to slowly moving boundaries; (i) is a quantity we can calculate, while (ii) is related to the basic properties of a superfluid via a simple thought experiment.

To connect the macroscopic, phenomenological description of superfluidity with our microscopic theory, we make use of the thought experiment shown schematically in Fig. 11(b): Consider an infinitely long box, containing superfluid, and accelerate the box along its long axis until it reaches a small velocity . Due to viscous interactions with the walls, such a box filled with a normal fluid should have a momentum density at equilibrium of . The notation denotes an expectation value in the ensemble with walls moving with velocity .

Because superfluids are nonviscous, the observed value for the momentum density in a superfluid is less than the value expected for a classical fluid. In Landau’s two-fluid model we attribute the observed momentum density, , to the “normal fraction” where is the normal fluid density. The superfluid fraction remains stationary in the lab frame, even at equilibrium and makes up the remaining mass with density .

In order to apply the usual procedures of statistical mechanics to the thought experiment, we consider two frames: the “lab frame” in which the walls move with velocity in the -direction and the “wall frame” in which the walls are at rest.

Figure 11: Thought experiment used in deriving the superfluid density. The walls move with velocity in the direction. To begin with, we imagine that the superfluid sits in a box of dimensions as shown in a). We later take the limit as the box walls recede to infinity to get the thermodynamic limit d). The order of the limits is critically important: the path b) leads to superflow while the path c) results in the entire fluid moving along with the walls.

Assuming that the fluid is in thermal equilibrium with the walls, the density matrix in the grand canonical ensemble is given by the usual expression where is the Hamiltonian of the system in the wall frame and . A Galilean transformation relates to the Hamiltonian in the lab frame, , where is the total momentum, is the total mass and is the momentum density operator at point . The expectation value for the momentum density in the presence of moving walls is then given by the expression


Expanding this expression to first order in yields


where all the expectation values on the right hand side are now taken in the equilibrium ensemble with the walls at rest. Since in our equilibrium ensemble, this simplifies to


where is a dyad [i.e., a rank-two tensor; the outer product of and ].

To make further progress, we wish to take the limit as the system gets very large (write “”). To this end, we first consider some properties of the correlation functions in the infinite system. The infinite system is homogeneous, which implies that for any , where indicates an average in the infinite system. As a consequence, we may express the correlations — in the infinite system — in terms of the Fourier transform in the relative coordinate :


where all the important features of the correlations are now captured by the tensor


Because of the isotropy of the fluid in the infinite system, obeys the transformation law for any orthogonal matrix . This implies that may be decomposed into the sum of longitudinal and transverse parts:


where , and is the identity. The transverse and longitudinal functions and are scalars which depend only on the length .

We now return our attention to the finite system. If the finite box is large then the momentum correlations in the bulk will be very similar to the values for the infinite system. Therefore, when and are far from the boundaries, we may approximate


which in combination with Eq. (54) yields


Here we have defined the nascent delta function which has the property as .

We are now in a position to carry out the limiting procedure to increase the box size to infinity. However, care must be taken because the simple expression is not well defined without further qualification of the limiting process .

To resolve this subtlety we must insert a final vital piece of physical reasoning. Let us assume for simplicity that is directed along the -direction, and the box is aligned with the and axes with dimensions . As shown in Fig. 11, there are two possibilities for taking the limits, representing different physical situations.

On the one hand [Fig. 11(b)], we may take the limit first, which gives us an infinitely long channel in which superfluid can remain stationary while only the normal fraction moves with the walls in the -direction. We have


where we use the fact that can be decomposed into the product with as . Employing the decomposition of given in Eq. (58) allows the density of the normal fraction to be related to the transverse component of evaluated at zero:


On the other hand [Fig. 11(c)] we may take the limit first, resulting in an infinitely long channel — with velocity perpendicular to the walls — in which the entire body of the fluid must move regardless of the superfluidity. In a similar way to the previous paragraph, , and making use of the decomposition in Eq. (58), the total density is related to the longitudinal component of the correlations:


With these expressions, the normal fraction may finally be expressed directly as


while the superfluid fraction is . Thus, we have expressed the superfluid and normal fractions in terms of a correlation function which can be directly computed from our simulation results.

d.2 Numerical procedure

To determine the superfluid fraction for our system, we need to estimate the tensor of momentum density correlations from the simulation results. For our finite system constrained to a periodic simulation box, we may compute the momentum correlations only at discrete grid points. The discrete analogue of Eq. (57) leads to the expression


where the constant of proportionality is not important to the final result, and are the discrete Fourier coefficients of over our simulation box.

The momentum density operator is given by


which may be derived by considering the continuity equation for the number density, . For a given classical field Eq. (10), the Fourier coefficients of may be written as


where is the area of the system. Computing a value for all at each time step, we then evaluate via the usual ergodic averaging procedure using Eq. (69).

Having evaluated , we are left with performing the decomposition into longitudinal and transverse parts. For this, simply note that Eq. (58) implies , and , where is a unit vector perpendicular to .

Values for and may be collected for all angles as a function of , and a fitting procedure used to perform the extrapolation ; this procedure is illustrated in Fig. 12. At low temperatures, the extrapolation is quite reliable, but becomes more difficult near the transition where sampling noise increases and changes rapidly near . Without a known functional form, we settled for a quadratic weighted least squares fit of and versus . A weighting of was used to counteract the fact that the density of samples of vs scales proportionally with