# Vortex pairing in two-dimensional Bose gases

###### Abstract

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.

###### pacs:

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.

## Ii Formalism

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

(1) |

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

(2) |

with the spatial extent of the mode ^{1}^{1}1For 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:

(3) |

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

(4) |

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 ^{2}^{2}2The 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

(5) |

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

(6) |

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

(7) |

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,

(8) | ||||

(9) |

with . The fields are thus defined by

(10) | ||||

(11) |

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

(14) |

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

(15) |

where the projection operator

(16) |

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

(17) | ||||

(18) |

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

(19) |

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

(20) |

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

(21) |

where

(22) |

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

(23) |

with corresponding atom number ; total number is simply

(24) |

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,

(25) | ||||

(26) |

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

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

(27) |

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

(28) |

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,

(29) |

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.

#### 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

(30) |

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

(31) |

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.

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 .

#### 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

(32) |

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

(33) |

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.

#### 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 ^{3}^{3}3For
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.

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 ^{4}^{4}4We 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.

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.

###### Acknowledgements.

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

(34) |

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:

(35) |

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

(36) |

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

(37) |

where . Rearranging, this is

(38) |

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

(39) |

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

(40) |

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:

(41) |

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:

(42) |

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

(43) | ||||

(44) |

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)

(45) |

This integral is of the general form

(46) |

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

(47) | ||||

(48) |

(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:

(49) |

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

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

(50) | ||||

(51) |

Expanding this expression to first order in yields

(52) |

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

(53) | ||||

(54) |

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 :

(55) | ||||

(56) |

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

(57) |

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:

(58) |

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

(59) | ||||

(60) |

which in combination with Eq. (54) yields

(61) | ||||

(62) |

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

(63) | ||||

(64) | ||||

(65) |

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:

(66) |

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:

(67) |

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

(68) |

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

(69) |

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

(70) |

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

(71) |

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