Vortex proliferation in the Berezinskii-Kosterlitz-Thouless regime on a
two-dimensional lattice of Bose-Einstein condensates
V. Schweikhard, S. Tung, and E. A. Cornell[*]
JILA, National Institute of Standards and Technology and University of Colorado,
and Department of Physics, University of Colorado, Boulder, Colorado 80309-0440
(Dated: October 30, 2018)
We observe the proliferation of vortices in the Berezinskii-Kosterlitz-Thouless regime on a two-
dimensional array of Josephson-coupled Bose-Einstein condensates. As long as the Josephson (tun-
neling) energy J exceeds the thermal energy T , the array is vortex-free. With decreasing J/T ,
vortices appear in the system in ever greater numbers. We confirm thermal activation as the vortex
formation mechanism and obtain information on the size of bound vortex pairs as J/T is varied.
PACS numbers: 03.75.Lm, 03.75.Gg, 74.50.+r, 74.81.Fa, 67.90.+z
One of the defining characteristics of superfluids is
long-range phase coherence [1], which may be destroyed
by quantum fluctuations, as in the Mott-insulator tran-
sition [2, 3], or thermal fluctuations, e.g. in one-
dimensional Bose gases [4, 5] and in a double-well system
[6]. In two dimensions (2D), Berezinskii [7], Kosterlitz
and Thouless [8] (BKT) developed an elegant description
of thermal phase fluctuations based on the unbinding of
vortex-antivortex pairs, i.e. pairs of vortices of opposite
circulation. The BKT picture applies to a wide variety
of 2D systems, among them Josephson junction arrays
(JJA), i.e. arrays of superfluids in which phase coherence
is mediated via a tunnel coupling J between adjacent
sites. Placing an isolated (free) vortex into a JJA is ther-
modynamically favored if its free energy F = E−TS ≤ 0.
In an array of period d the vortex energy diverges with ar-
ray size R as E ≈ J log(R/d) [9], but may be offset by an
entropy gain S ≈ log(R/d) due to the available ≈ R2/d2
sites. This leads to a critical condition (J/T )crit ≈ 1 in-
dependent of system size, below which free vortices will
proliferate. In contrast, tightly bound vortex-antivortex
pairs are less energetically costly and show up even above
(J/T )crit. The overall vortex density is thus expected to
grow smoothly with decreasing J/T in the BKT crossover
regime.
Transport measurements, both in continuous superflu-
ids [10, 11] and superconducting JJA [12] have confirmed
the predictions of BKT, without however directly observ-
ing its microscopic mechanism, vortex-antivortex unbind-
ing. A recent beautiful experiment [13] in a continuous
2D Bose gas measured the phase-phase decay function
through the BKT cross-over, and saw evidence for ther-
mal vortex formation. For related theoretical studies see
e.g. [14]. In this work we present more detailed vortex-
formation data, collected in a 2D array of BECs with ex-
perimentally controllable Joephson couplings. The sys-
tem was studied theoretically in [15].
Our experiment starts with production of a partially
Bose-condensed sample of 87Rb atoms in a harmonic, ax-
ially symmetric magnetic trap with oscillation frequen-
cies {ωx,y, ωz} = 2π{6.95, 15.0}Hz. The number of con-
densed atoms is kept fixed around 6 × 105 as the tem-
perature is varied. We then transform this system into
4.7µm
. . ..
. . ..
... ..
Nwell
FIG. 1: Experimental system. (a) 2D optical lattice intensity
profile. A lattice of Josephson-coupled BECs is created in the
white-shaded area. The central box marks the basic building
block of our system, the double-well potential shown in (b).
The barrier height VOL and the number of condensed atoms
per well, Nwell, control the Josephson coupling J , which acts
to lock the relative phase ∆φ. A cloud of uncondensed atoms
at temperature T induces thermal fluctuations and phase de-
fects in the array when J < T . (c) Experimental sequence: A
BEC (i) is loaded into the optical lattice over 10 s, suppress-
ing J to values around T. We allow 2 s for thermalization. To
probe the system, we ramp off the lattice on a faster timescale
tr [23] and take images of the recombined condensate. When
J is reduced below T (ii)-(vii), vortices (dark spots) appear
as remnants of the thermal fluctuations in the array.
a Josephson junction array, as illustrated in Fig. 1. In
a 10 s linear ramp, we raise the intensity of a 2D hexag-
onal optical lattice [16] of period d = 4.7µm in the x-y
plane. The resulting potential barriers of height VOL be-
tween adjacent sites [Fig.1(b)] rise above the condensate’s
chemical potential around VOL ≈ 250− 300Hz, splitting
it into an array of condensates which now communicate
http://arxiv.org/abs/0704.0289v3
only through tunneling. This procedure is adiabatic even
with respect to the longest-wavelength phonon modes of
the array [18, 19] over the full range of VOL in our exper-
iments. Each of the ≈ 190 occupied sites (15 sites across
the BEC diameter 2×RTF ≈ 68µm [20]) now contains a
macroscopic BEC, withNwell ≈ 7000 condensed atoms in
each of the central wells at a temperature T that can be
adjusted between 30−70nK. By varying VOL in a range
between 500Hz and 2 kHz we tune J between 1.5µK
and 5nK, whereas the “charging” energy Ec, defined in
[1], is on the order of a few pK, much smaller than both
J and T . In this regime, thermal fluctuations of the rela-
tive phases ∆φTh ≈
T/J are expected, while quantum
fluctuations ∆φQ ≈ (Ec/4J)1/4 are negligible [1].
The suppression of the Josephson coupling greatly sup-
presses the energy cost of phase fluctuations in the x-y
plane, between condensates, J [1− cos(∆φ)], compared to
the cost of axial (z) phase fluctuations inside the con-
densates [21]. As a result, axial phase fluctuations re-
main relatively small, and each condensate can be ap-
proximated as a single-phase object [22].
After allowing 2 s for thermalization, we initiate our
probe sequence. We first take a nondestructive thermom-
etry image in the x-z plane, from which the temperature
T and, from the axial condensate size Rz, the number
of condensed particles per well, Nwell, is obtained (see
below). To observe the phase fluctuations we then turn
down the optical lattice on a time-scale tr [23], which
is fast enough to trap phase winding defects, but slow
enough to allow neighboring condensates to merge, pro-
vided their phase difference is small. Phase fluctuations
are thus converted to vortices in the reconnected conden-
sate, as has been observed in the experiments of Scherer
et al. [24]. We then expand the condensate by a factor
of 6 and take a destructive image in the x-y plane.
Figure 1(c) illustrates our observations: (ii)-(vii) is a
sequence of images at successively smaller J/T (measured
in the center of the array [25]). Vortices, with their cores
visible as dark “spots” in (iii)-(vii), occur in the BEC
center around J/T = 1. Vortices at the BEC edge ap-
pear earlier, as here the magnetic trap potential adds
to the tunnel barrier, suppressing the local J/T below
the quoted value. That the observed “spots” are indeed
circulation-carrying vortices and antivortices is inferred
from their slow ≈ 100ms decay after the optical lattice
ramp-down, presumably dominated by vortex-antivortex
annihilation. From extensive experiments on vortices in
our system we know that circulation-free “holes” fill so
quickly due to positive mean field pressure, that they
do not survive the pre-imaging expansion. Vortices with
identical circulation would decay by dissipative motion
to the BEC edge, in our trap over & 10 s.
To investigate the thermal nature of phase fluctua-
tions, we study vortex activation while varying J at
different temperatures. For a quantitative study, accu-
rate parameter estimates are required. The Josephson-
coupling energy J is obtained from 3D numerical sim-
ulations of the Gross-Pitaevskii equation (GPE) for
0.1 1 10
0.1 1 10 100
 cold: 30<T<40nK
 hot: 55<T<70nK
10 100 1000
J [nK]
 cold: 30<T<40nK
 hot: 55<T<70nK
10 100 1000
J [nK]
(b) (c)
FIG. 2: Quantitative study of vortices. The areal density of
vortices is quantified by the plotted D defined in the text. D is
extracted only from the central 11% of the condensate region
[circle in inset (a)] to minimize effects of spatial inhomogene-
ity. (b) D vs J for two datasets with distinct “cold” and “hot”
temperatures. Each point represents one experimental cycle.
The increase in D with decreasing J . 100nK signals the
spontaneous appearance of vortices, while the “background”
D . 0.01 for J & 200nK is not associated with vortices.
Vortices clearly proliferate at larger J for the “hot” data,
indicating thermal activation as the underlying mechanism.
The large scatter in D at low J is due to shot noise on the
small average number of vortices in the central condensate
region. (c) same data as in (b), but averaged within bins of
size ∆[log(J)] = 0.15. Error bars of D are standard errors.
(d) same data as (b), but plotted vs J/T . “Cold” and “hot”
datasets almost overlap on what appears to be a universal
vortex activation curve, as confirmed by averaging [inset (e)],
clearly revealing the underlying competition of J and T .
the central double-well system [6, 26] [Fig.1(b)], self-
consistently including mean-field interactions of both
condensed and uncondensed atoms [27]. A use-
ful approximation for J in our experiments is [25]:
J(VOL, Nwell, T ) ≈ Nwell × 0.315nK exp[Nwell/3950 −
VOL/244Hz](1+0.59T/100nK). The finite-T correction
to J arises from both the lifting-up of the BEC’s chem-
ical potential and the axial compression by the thermal
cloud’s repulsive mean field, but does not take into ac-
count the effects of phase fluctuations on J (in condensed-
matter language, we calculate the bare J). Nwell is de-
termined by comparison of the experimentally measured
Rz, to Rz(VOL, Nwell, T ) obtained from GPE simula-
tions. Both experimental and simulated Rz are obtained
from a fit to the distribution of condensed and uncon-
densed atoms, to a Thomas-Fermi profile plus mean-field-
modified Bose function [27]. In determination of all J
values, there is an overall systematic multiplicative un-
certainty ∆J/J = ×
1.6, dominated by uncertainties in
the optical lattice modulation contrast, the absolute in-
tensity calibration, and magnification in the image used
to determine Nwell. In comparing J for “hot” and “cold”
clouds (see Fig. 2) there is a relative systematic error of
15% associated with image fitting and theory uncertain-
ties in the thermal-cloud mean-field correction to J .
The qualitative results of our work are consistent
whether we use an automated vortex-counting routine
or count vortices by hand, but the former shows signs of
saturation error at high vortex density, and the latter is
vulnerable to subjective bias. As a robust vortex-density
surrogate we therefore use the “roughness” D of the con-
densate image caused by the vortex cores. Precisely, we
define D as the normalized variance of the measured col-
umn density profile from a fit to a smooth finite-T Bose
profile [27], with a small constant offset subtracted to
account e.g. for imaging noise. To limit spatial inhomo-
geneity in J , caused by spatially varying condensate den-
sity and optical lattice intensity, to < 10%, D is extracted
only from the central 11% of the condensate area which
contains 20 lattice sites [Fig. 2(a)]. Comparison to au-
tomated vortex-counts shows that D is roughly linear in
the observed number of vortices, irrespective of the sign
of their circulation, with a sensitivity of ≈ 0.01/vortex.
Figure 2 shows results of our quantitative study. In
Fig. 2(b), we plot D vs J for two datasets with dis-
tinct temperatures. At large J & 200nK a background
D . 0.01 is observed, that is not associated with vortices,
but due to residual density ripples remaining after the
optical lattice ramp-down. Vortex proliferation, signaled
by a rise of D above ≈ 0.01, occurs around J ≈ 100nK
for “hot” BECs and at a distinctly lower J ≈ 50nK for
“cold” BECs [confirmed by the averaged data shown in
Fig. 2(c)], indicating thermal activation as the vortex
formation mechanism. Plotting the same data vs J/T
in Fig. 2(d) shows collapse onto a universal vortex acti-
vation curve, providing strong evidence for thermal ac-
tivation. A slight residual difference becomes visible in
the averaged “cold” vs “hot” data [Fig.2(e)], perhaps be-
cause of systematic differences in our determination of J
at different temperatures.
The vortex density D by itself provides no distinction
between bound vortex-antivortex pairs and free vortices.
In the following we exploit the flexibility of optical po-
tentials to distinguish free or loosely bound vortices from
tightly bound vortex-antivortex pairs. A “slow” optical
lattice ramp-down allows time for tightly bound pairs
to annihilate before they can be imaged. By slowing
down the ramp-down duration τ [inset of Fig. 3 (a)],
0.1 1 10
ramp-down timescale τ
 τ = 6 ms
 τ = 36 ms
0 20 40 60
smallest surviving pairs
ramp-down timescale τ [ms]
      (J/T<<1)
      (J/T>>1)
0 20 40 60
ramp-down timescale τ [ms]
1.3kHz
(c) I II & III larger pairs & free vortices
FIG. 3: (a) Vortex density D probed at different optical lat-
tice ramp-down timescales τ . A slow ramp provides time
for tightly bound vortex-antivortex pairs to annihilate, allow-
ing selective counting of loosely bound or free vortices only,
whereas a fast ramp probes both free and tightly bound vor-
tices. A fit to the vortex activation curve determines its mid-
point (J/T )50%, its 27% − 73% width ∆(J/T )27−73, and the
limiting values D< (D>) well below (above) (J/T )50%. (b)
A downshift in (J/T )50% is seen for slow ramps, consistent
with the occurrence of loosely bound or free vortices at lower
J/T only. (c) Mapping between ramp-down timescale τ and
estimated size of the smallest pairs surviving the ramp (up-
per axis). The difference D< − D> measures the number of
vortices surviving the ramp (right axis). Comparison to sim-
ulated vortex distributions yields a size estimate of the small-
est surviving pairs (upper axis). Inset: smallest possible pair
sizes in a hexagonal array, I: d/
3, II: d, III: 2d/
we therefore selectively probe vortices at increasing spa-
tial scales. This represents an attempt to approach the
“true” BKT vortex unbinding crossover that is comple-
mentary to transport measurements employed so success-
fully in superconductive and liquid Helium systems.
Figure 3(a) shows vortex activation curves, probed
with two different ramp-down times. Two points are
worth noticing: First, a slow ramp compared to a fast
one shows a reduction of the vortex density D< in arrays
with fully randomized phases at low J/T . The difference
directly shows the fraction of tightly bound pairs that an-
nihilate on the long ramp. Second, a slower ramp shows
vortex activation at lower (J/T )50%, confirming that free
or very loosely bound vortices occur only at higher T
(lower J). Specifically, the data clearly show a range
around J/T ≈ 1.4 where only tightly bound pairs exist.
Figure 3(b) quantitatively shows the shift of (J/T )50%
from 1.4 to 1.0 with slower ramp time. We can make
a crude mapping of the experimental ramp-down time-
scale to theoretically more accessible vortex-antivortex
pair sizes as follows: In Fig. 3(c), we see the decrease of
the saturated (low-J/T ) vortex density D< with increas-
ing ramp timescale τ . The right axis shows the inferred
number of vortices that survived the ramp. We compare
this number of surviving vortices to simulations [28] of
a 20-site hexagonal array with random phases. In these
simulations we find, on average, a total of 10 vortices,
6 of which occur in nearest-neighbor vortex-antivortex
pairs [configuration I in Fig. 3(c)], 1.7 (0.4) occur in
configuration II (III) respectively, and 1.9 occur in larger
pairs or as free vortices. Experimentally≈ 11 vortices are
observed for the fastest ramps, in good agreement with
the expected total number of vortices. For just some-
what slower ramps of τ ≈ 5ms, only 3 vortices survive,
consistent with only vortices in configuration II & III or
larger remaining (indicated in Fig. 3, top axis) [29]. For
τ & 30ms ramps less than 2 vortices remain, according
to our simulations spaced by more than 2d/
3. Thus we
infer that ramps of τ ≈ 30ms or longer allow time for
bound pairs of spacing . 2d/
3 to decay before we ob-
serve them. The downward shift of (J/T )50% in Fig. 3(b)
thus tells us that loosely bound pairs of size larger than
3, or indeed free vortices, do not appear in quan-
tity until J/T ≤ 1.0, whereas more tightly bound vortex
pairs appear in large number already for J/T ≤ 1.4.
A further interesting observation concerns the width
of the vortex activation curve. The relative width, de-
termined from fits to data such as the ones shown in
Fig. 3(a), is ∆(J/T )27−73/(J/T )50% ≈ 0.3, independent
of ramp-down duration. This width is neither as broad
as in a double-well system [6, 30], where the coherence
factor rises over a range ∆(J/T )27−73/(J/T )50% ≈ 1.4,
nor as broad as expected from our simulations [28] of
an array of uncoupled phases, each fluctuating inde-
pendently with ∆φRMS =
T/J , for which we find
∆(J/T )27−73/(J/T )50% ≈ 0.85. Presumably collective
effects in the highly multiply connected lattice narrow
the curve. On the other hand, the width is 3 times larger
than the limit due to spatial inhomogeneity in J , suggest-
ing contributions to the width due to finite-size effects
or perhaps revealing the intrinsically smooth behavior of
vortex activation in the BKT regime.
In conclusion, we have probed vortex proliferation in
the BKT regime on a 2D lattice of Josephson-coupled
BECs. Allowing variable time for vortex-antivortex pair
annihilation before probing the system provides a time-
to-length mapping, which reveals information on the size
of pairs with varying J/T . We acknowledge illuminating
conversations with Leo Radzihovsky and Victor Gurarie.
This work was funded by NSF and NIST.
[*] Quantum Physics Division, National Institute of Stan-
dards and Technology.
[1] A. Leggett, Rev. Mod. Phys. 73, 307 (2001).
[2] M. Greiner et al., Nature 415, 39 (2002).
[3] D. Jaksch et al., Phys. Rev. Lett. 81, 003108 (1998).
[4] S. Richard et al., Phys. Rev. Lett. 91, 010405 (2003);
[5] D. Hellweg et al., Phys. Rev. Lett. 91, 010406 (2003).
[6] R. Gati et al., Phys. Rev. Lett. 96, 130404 (2006).
[7] V. Berezinskii, Sov. Phys.-JETP 32, 493 (1971); 34, 610
(1972).
[8] J. Kosterlitz, D. Thouless, J. Phys. C 6, 1181 (1973).
[9] M. Tinkham, Introduction to Superconductivity,
McGraw-Hill, Inc., New York (1996).
[10] G. Agnolet et al., Phys. Rev. B 39, 8934 (1989).
[11] A. Safonov et al., Phys. Rev. Lett. 81, 4545 (1998)
[12] D. Resnick et al., Phys. Rev. Lett. 47, 1542 (1981).
[13] Z. Hadzibabic et al., Nature 441, 1118 (2006).
[14] T. Simula and P. Blakie, Phys. Rev. Lett. 96, 020404
(2006).
[15] A. Trombettoni et al., New J. Phys. 7, 57 (2005).
[16] Three circularly polarized laser beams (λ = 810.1nm) in-
tersect in a tripodlike configuration, with θ = 6.6◦ angles
to the z-axis. Calculation of the optical dipole potential
[17] includes counterrotating terms and interaction with
both the D1 and D2 lines, as well as the ‘fictitious mag-
netic field’ due to the circular polarization. The tilted
bias field of the TOP trap makes P ≈ 0.5[17].
[17] R. Grimm et al., Adv. At. Mol. Opt. Phys. 42, 95 (2000).
[18] J. Javanainen, Phys. Rev. A 60, 4902 (1999).
[19] K. Burnett et al., J. Phys. B 35, 1671 (2002).
[20] To avoid radial flows during VOL ramp-up, RTF is kept
constant by balancing the lattice-enhanced mean field
pressure with radial confinement due to the optical lattice
envelope, by the choice of a 67µm 1/e2 intensity waist.
Axial (z) confinement is due to the magnetic trap alone.
[21] D. Petrov et al., Phys. Rev. Lett. 87, 050404 (2001).
[22] In the axial condensate region between z =
−Rz/3,+Rz/3, where according to our 3D GPE
simulations 85% of the tunnel current is localized
and hence the relative phase is measured, axial phase
fluctuations [21] vary between ≈ 600mrad (“cold” data
in Fig. 2) and 800mrad (“hot”) in the regime J/T ≈ 1.
[23] Within a dataset, the ramp-down rate is kept fixed, tr =
τ × VOL/1.3 kHz, τ = 18ms if not otherwise indicated.
[24] D. Scherer et al., Phys. Rev. Lett. 98, 110402 (2007).
[25] J is averaged over junctions within the central 11% of
the array area, from which all quantitative experimental
results are extracted.
[26] D. Ananikian and T. Bergeman, Phys. Rev. A 73, 013604
(2006).
[27] M. Naraschewski and D. Stamper-Kurn, Phys. Rev. A
58, 2423 (1998).
[28] Following [24], in simulations we count a vortex if all
three phase differences in an elemental triangle of junc-
tions are ∈ (0, π), or if all are ∈ (−π, 0).
[29] The very short annihilation time of configuration I pairs
[ τ < 5ms in Fig. 3(c) ] is not unexpected: their spacing
3 = 2.8µm is comparable to the diameter of a vortex
core in the bulk condensate (after VOL rampdown).
[30] L. Pitaevskii and S. Stringari, Phys. Rev. Lett. 87,
180402 (2001).

  We observe the proliferation of vortices in the
Berezinskii-Kosterlitz-Thouless regime on a two-dimensional array of
Josephson-coupled Bose-Einstein condensates. As long as the Josephson
(tunneling) energy J exceeds the thermal energy T, the array is vortex-free.
With decreasing J/T, vortices appear in the system in ever greater numbers. We
confirm thermal activation as the vortex formation mechanism and obtain
information on the size of bound vortex pairs as J/T is varied.

Introduction to Superconductivity,
McGraw-Hill, Inc., New York (1996).
[10] G. Agnolet et al., Phys. Rev. B 39, 8934 (1989).
[11] A. Safonov et al., Phys. Rev. Lett. 81, 4545 (1998)
[12] D. Resnick et al., Phys. Rev. Lett. 47, 1542 (1981).
[13] Z. Hadzibabic et al., Nature 441, 1118 (2006).
[14] T. Simula and P. Blakie, Phys. Rev. Lett. 96, 020404
(2006).
[15] A. Trombettoni et al., New J. Phys. 7, 57 (2005).
[16] Three circularly polarized laser beams (λ = 810.1nm) in-
tersect in a tripodlike configuration, with θ = 6.6◦ angles
to the z-axis. Calculation of the optical dipole potential
[17] includes counterrotating terms and interaction with
both the D1 and D2 lines, as well as the ‘fictitious mag-
netic field’ due to the circular polarization. The tilted
bias field of the TOP trap makes P ≈ 0.5[17].
[17] R. Grimm et al., Adv. At. Mol. Opt. Phys. 42, 95 (2000).
[18] J. Javanainen, Phys. Rev. A 60, 4902 (1999).
[19] K. Burnett et al., J. Phys. B 35, 1671 (2002).
[20] To avoid radial flows during VOL ramp-up, RTF is kept
constant by balancing the lattice-enhanced mean field
pressure with radial confinement due to the optical lattice
envelope, by the choice of a 67µm 1/e2 intensity waist.
Axial (z) confinement is due to the magnetic trap alone.
[21] D. Petrov et al., Phys. Rev. Lett. 87, 050404 (2001).
[22] In the axial condensate region between z =
−Rz/3,+Rz/3, where according to our 3D GPE
simulations 85% of the tunnel current is localized
and hence the relative phase is measured, axial phase
fluctuations [21] vary between ≈ 600mrad (“cold” data
in Fig. 2) and 800mrad (“hot”) in the regime J/T ≈ 1.
[23] Within a dataset, the ramp-down rate is kept fixed, tr =
τ × VOL/1.3 kHz, τ = 18ms if not otherwise indicated.
[24] D. Scherer et al., Phys. Rev. Lett. 98, 110402 (2007).
[25] J is averaged over junctions within the central 11% of
the array area, from which all quantitative experimental
results are extracted.
[26] D. Ananikian and T. Bergeman, Phys. Rev. A 73, 013604
(2006).
[27] M. Naraschewski and D. Stamper-Kurn, Phys. Rev. A
58, 2423 (1998).
[28] Following [24], in simulations we count a vortex if all
three phase differences in an elemental triangle of junc-
tions are ∈ (0, π), or if all are ∈ (−π, 0).
[29] The very short annihilation time of configuration I pairs
[ τ < 5ms in Fig. 3(c) ] is not unexpected: their spacing
3 = 2.8µm is comparable to the diameter of a vortex
core in the bulk condensate (after VOL rampdown).
[30] L. Pitaevskii and S. Stringari, Phys. Rev. Lett. 87,
180402 (2001).
