Bose-Einstein Condensation from a Rotating Thermal Cloud: Vortex Nucleation and Lattice Formation

# Bose-Einstein Condensation from a Rotating Thermal Cloud: Vortex Nucleation and Lattice Formation

A. S. Bradley ARC Centre of Excellence for Quantum-Atom Optics, School of Physical Sciences, University of Queensland, Brisbane, QLD 4072, Australia.    C. W. Gardiner Jack Dodd Centre for Quantum Technology, Department of Physics, University of Otago, PO Box 56, Dunedin, New Zealand.    M. J. Davis ARC Centre of Excellence for Quantum-Atom Optics, School of Physical Sciences, University of Queensland, Brisbane, QLD 4072, Australia.
July 21, 2019
###### Abstract

We develop a stochastic Gross-Pitaveskii theory suitable for the study of Bose-Einstein condensation in a rotating dilute Bose gas. The theory is used to model the dynamical and equilibrium properties of a rapidly rotating Bose gas quenched through the critical point for condensation, as in the experiment of Haljan et al. [Phys. Rev. Lett., 87, 21043 (2001)]. In contrast to stirring a vortex-free condensate, where topological constraints require that vortices enter from the edge of the condensate, we find that phase defects in the initial non-condensed cloud are trapped en masse in the emerging condensate. Bose-stimulated condensate growth proceeds into a disordered vortex configuration. At sufficiently low temperature the vortices then order into a regular Abrikosov lattice in thermal equilibrium with the rotating cloud. We calculate the effect of thermal fluctuations on vortex ordering in the final gas at different temperatures, and find that the BEC transition is accompanied by lattice melting associated with diminishing long range correlations between vortices across the system.

###### pacs:
03.75.Kk,03.75.Lm,05.70.Fh,05.10.Gg,

## I Introduction

Rotating Bose-Einstein condensation was first realized by Haljan et al.at JILA Haljan2001 (), who grew large vortex lattices by evaporatively cooling a rotating thermal cloud. The evaporation is performed in a prolate trap, along the axis of rotation so that escaping atoms make large axial excursions but do not stray far from the rotation axis. Escaping atoms carry away significant energy but little angular momentum, increasing the angular momentum per trapped atom. Using this cooling mechanism, which requires typically s of evaporative cooling and an extremely cylindrical trap, very large vortex lattices containing more than 300 vortices have been created Engels2004 (). This nucleation mechanism is distinct from more studied vortex formation mechanisms Matthews1999 (); Leonhardt2002 (); Madison2000 (); Hodby2002 (); Abo-Shaeer2001 (); Inouye2001 (); Anderson2001a (); Dutton2001a (); Scherer2007 () because, rather than manipulating a condensate to generate vortices, vortices can be formed at the onset of condensation as the rotating gas reaches the critical temperature () for Bose-Einstein condensation. This mechanism thus has connections with the Kibble-Zurek Kibble1976 (); Zurek1985 () and Berezinskii-Kosterlitz-Thouless Berezinskii1971 (); Kosterlitz1973 (); Hadzibabic2006a () mechanisms, as it is fundamentally caused by spontaneous symmetry breaking and thermal fluctuations near a phase transition.

In this work we study the dynamics of vortex lattice formation during rotating evaporative cooling using a stochastic Gross-Pitaevskii equation (SGPE) Stoof1999 (); SGPEI (); SGPEII (). Our central aim is to include rotation, thermal fluctuations, and two-body collisions, by separating the system into a condensate band coupled to a (possibly rotating) thermal reservoir. The condensate band is described by a generalized mean field theory known as the truncated Wigner method Graham1973 (); Steel1998 (); Sinatra2001 () that includes the projected classical field method Davis2001b (); Davis2006a (); Simula2006a () as a special case in the regime of high temperatures, while the non-condensate band formalism is based on quantum kinetic theory QKV ().

We also provide expressions for the dissipation rates of the theory under the condition that the high-energy component is approximately Bose-Einstein distributed, and details of our numerical algorithm that is central to the utility of the SGPE approach. These are the first SGPE simulations of vortex formation in a trapped Bose gas.

## Ii System and equations of motion

The trapped system is divided into a condensate band of states with energy beneath the cut-off , and the remaining non-condensate band of high-energy states. The formal separation of the Hamiltonian into bands according to energy allows different methods to applied in the two bands. The condensate band is treated using the truncated Wigner phase space method Graham1973 (); Steel1998 () which is valid when the modes are significantly occupied Blakie2007a (). Within this approach the condensate band includes the condensate and the most highly occupied excitations. The non-condensate band is comprised of many weakly occupied modes and is treated as approximately thermalised; the effect of this band is to generate dissipative terms in the equation of motion for the condensate band. The high-energy states play the role of a grand canonical thermal reservoir for the low energy condensate band. In this section we provide a minimal overview of the formalism and simply state the final equations of motion to be used in this work. Refs. SGPEI (); SGPEII () contain the essential formalism, and Ref. SGPEIII () gives necessary modifications for the rotation system, along with a number of corrections to the derivation of the equations of motion.

### ii.1 Hamiltonian and separation of the rotating system

The cold-collision regime for bosonic atoms of mass is characterized by the interaction parameter , where is the s-wave scattering length Dalfovo1999 (). The second-quantized many-body Hamiltonian for a dilute Bose gas confined by a trapping potential in the rotating frame defined by angular frequency is , where

 ^Hsp=∫d3x^ψ†(x)(−ℏ2∇22m−ΩLz+V(x,t))^ψ(x), (1)

and

 ^HI=u2∫d3x^ψ†(x)^ψ†(x)^ψ(x)^ψ(x). (2)

The orthogonal projectors that separate the system are defined with respect to the single particle basis of the trap. This provides a convenient basis because (a) The many body spectrum becomes well approximated at high-energy by a single particle spectrum, allowing the cutoff to be imposed in this basis Davis2001b (), and (b) the numerical method used to solve the final equations of motion uses the single particle basis Blakie2005a ().

The condensate band projector takes the form

 P≡−∑n|n⟩⟨n|, (3)

where the bar denotes the high-energy cut-off and represents all quantum numbers required to index the eigenstates of the single particle Hamiltonian. To minimize notation we will use both for the operator and for its action on wave functions; the meaning will always be clear from the context. The field operator is decomposed as

 ^ψ(x)=^ϕ(x)+^ψNC(x)≡P^ψ(x)+Q^ψ(x), (4)

where the non-condensate field describes the high-energy thermal modes, and the condensate band is described by the field . In the spatial representation the eigensolutions of , denoted , generate the action of the projector on an arbitrary wave function

 PΨ(x)=−∑nYn(x)∫d3zY∗n(z)Ψ(z). (5)

The condensate band field theory generated by such a decomposition has a nonlocal commutator

 [^ϕ(x),^ϕ†(y)]=δC(x,y) (6)

where . In the harmonic oscillator representation of a one dimensional system, the approximate delta function has a position dependent width, being narrowest in the center of the trap. This reflects the fact that trap wave-functions become smoother near their semi-classical turning points. Note that since , ; put another way, this means that is a true Dirac-delta function for any condensate band field.

There are two physical processes arising from the interactions between condensate and non-condensate bands, referred to as growth and scattering SGPEI (); SGPEII (); SGPEIII (), and shown schematically in Figure 2. The growth terms correspond to two-body collisions between condensate and non-condensate atoms whereby particle transfer can take place. The scattering describes inter-band collisions that conserve the populations in each band but involve a transfer of energy and momentum between bands.

The non-condensate band is described semi-classically in terms of the Wigner function

 F(x,K)≡∫d3v⟨^ψ†NC(x+v2)^ψNC(x−v2)⟩eiK⋅v. (7)

In this work the trapping potential assumes the cylindrical form

 V(x,t)≡m2(ω2rr2+ω2zz2), (8)

where and are the radial and axial trapping frequencies respectively. The non-condensate band Wigner function for a Bose gas confined in a cylindrically symmetric trap in equilibrium in the rotating frame at frequency , with chemical potential , is

 F(x,K)=1exp[(ℏω(x,K)−μ)/kBT]−1, (9)

where

 ℏω(x,K)=ℏ2K22m−ℏΩ⋅(x×K)+V(x),=ℏ22m(K−mΩ×x/ℏ)2+Veff(x), (10)

defines the effective potential as , and is the distance from the symmetry axis of the trap. Since is the rotating frame momentum, we can work in terms of this variable, defining

 ℏω(x,K)=ℏωr(x,Kr)≡ℏ2K2r2m+Veff(x). (11)

Hereafter we drop the subscript for brevity since we always describe the system in this frame. In the semiclassical approximation is only significant in the non-condensate region , which is the region of phase space where .

### ii.2 Stochastic Gross-Pitaevskii equation

Following the theory in Refs. SGPEI (); SGPEII (); SGPEIII (), the equation of motion for the density operator describing the total system

 iℏ∂^ρ∂t=[^H,^ρ], (12)

is reduced to a master equation for the condensate band density operator, which is the trace of over the non-condensate band degrees of freedom . The master equation is then mapped to a Fokker-Planck equation for the Wigner distribution describing the condensate band. In the truncated Wigner approximation, an equivalent stochastic equation of motion can be obtained that allows the dynamics to be efficiently simulated numerically QN (). Formally, there is a correspondence between symmetrically ordered moments of the field operator and ensemble averages over many trajectories of the projected field that evolves according to the SGPE to be defined below SGPEIII (). For example, for the density

 ⟨α∗(x)α(x)⟩W=12⟨^ϕ†(x)^ϕ(x)+^ϕ(x)^ϕ†(x)⟩, (13)

where denotes stochastic averaging over trajectories of the SGPE.

To state the SGPE we require some expressions from mean field theory, in particular the GPE obeyed by in the absence of any thermal atoms in the non-condensate band. The mean field Hamiltonian for corresponding to is

 HGP=∫d3xα∗(x)(−ℏ2∇22m−ΩLz+V(x))α(x)+u2∫d3x|α(x)|4, (14)

with atom number given by .

The growth terms involve the Gross-Pitaevskii equation found in the usual way as

 iℏ∂α(x)∂t=¯δHGP¯δα∗(x)≡PLGPα(x)=P{(−ℏ2∇22m−ΩLz+V(x)+u|α(x)|2)α(x)} (15)

that defines the evolution operator . To obtain this equation of motion we have used the projected functional calculus SGPEII () that plays the same role here as ordinary functional calculus in non-projected field theory.

The scattering process couples to the divergence of the condensate band current

 jGP(x)≡iℏ2m([∇α∗(x)]α(x)−α∗(x)∇α(x))−(Ω×x)|α(x)|2, (16)

where the second line is a rigid body rotation term arising from the transformation to the rotating frame. The limit gives the laboratory frame velocity field , corresponding to the irrotational system mimicking rigid body rotation. Persistent currents around the vortex cores are responsible for maintaining irrotationality.

#### ii.2.1 Full non-local SGPE

We finally arrive at the Stratonovich stochastic differential equation (SDE) SGPEIII ()

 (S)dα(x,t)=−iℏPLGPα(x)dt+P{G(x)kBT(μ−LGP)α(x)dt+dWG(x,t)+iℏα(x)kBT∫d3x′M(x−x′)∇⋅jGP(x′)dt+iα(x)dWM(x,t)}, (17)

where the growth noise is complex, and the scattering noise is real. The noises are Gaussian, independent of each other, and have the non-zero correlations

 ⟨dW∗G(x,t)dWG(x′,t)⟩ = 2G(x)δC(x,x′)dt, (18) ⟨dWM(x,t)dWM(x′,t)⟩ = 2M(x−x′)dt. (19)

The growth and scattering amplitudes, and , are calculated in Appendix A for the case where the non-condensate band is described by a thermal Bose-Einstein distribution.

The two validity conditions for this equation are (i) high temperature: where is the characteristic system frequency, and (ii) significant occupation of modes in the condensate band, necessitated by the truncated Wigner method. The first line of Eq. 17 is the PGPE developed by Davis et al Davis2001b () that, in essence, is the GPE formally confined to a low energy subspace of the condensate band. The condensate growth terms of the second line are closely connected to phenomenological dissipative theories of BECs based on the GPE Choi1998 (); Penckwitt2002 (); Tsubota2002 () but here have a specific form given in Appendix A. The remaining terms describe scattering processes that transfer energy and momentum between the two bands; these terms where first derived in Ref. SGPEII (), and have not appeared in previous SGPE theories.

Irrespective of the form of and , the equilibrium state for the Wigner distribution of the system has the grand canonical form SGPEII ()

 Ws∝exp(μNGP−HGPkBT), (20)

that is generated by the stochastic evolution of Eq. (17). Solutions of the SGPE also known to satisfy a set of generalized Ehrenfest relations which in the steady state reduce to a statement of the fluctuation-dissipation theorem for the grand canonical Bose gas Bradley2005b ().

#### ii.2.2 Simple growth SGPE

A simpler equation that reaches the same equilibrium state can be obtained by noting that the scattering terms in Eq. (17) do not directly change the condensate band population. Such terms have been included in a kinetic theory of Bose-Einstein condensation QKIV (), where they were found to play a role in the early stages of condensate growth by causing atoms in many modestly occupied states to merge more rapidly into a single highly occupied state which then dominates the growth process. Within the SGPE theory, the scattering terms generate an effective potential from density fluctuations. This can be seen by noting that upon neglecting the growth and scattering terms, that are relatively small compared to the mean field evolution, the continuity equation for the condensate band gives . Consequently, since is purely real, the deterministic part of the scattering terms generate an effective potential that will tend to smooth out density fluctuations.

Upon neglecting these terms the resulting stochastic equation is much easier to integrate numerically as it does not require the non-local noise correlations of the full equations. The simulations in this work are performed by numerically integrating the simple growth SDE

 dα(x)=−iℏPLGPα(x)dtP{γkBT(μ−LGP)α(x)dt+dWγ(x,t)}, (21)

using the mixed quadrature pseudo-spectral method detailed in Appendix B. The rate is given by Eq. (44), and the noise is complex Gaussian with non-vanishing correlator

 ⟨dW∗γ(x,t)dWγ(x′,t)⟩=2γδC(x,x′)dt. (22)

## Iii Simulations of Rotating Bose-Einstein condensation

The central question of interest for this work regards the formation of vortices during condensation. In particular, does a vortex-free condensate form and then vortices enter it from the edge, or does condensation proceed into a state with vortices already present? That this is non-trivial to answer is easily seen from the single particle energy spectrum of the cylindrically symmetric harmonic trap in three dimensions

 ϵnlm=ℏωr(2n+|l|+1)−ℏΩl+ℏωz(m+1/2), (23)

where are the radial, angular and axial quantum numbers. In the absence of interactions BEC will form in the ground state mode with energy . Vortices arise from occupation of states with nonzero angular momentum and the positive angular momentum part of the spectrum behaves as leading to near-degeneracy of positive angular momentum modes as . States with angular momentum are increasingly easy to excite as the rate of rotation increases and in the interacting gas there is a trade-off between the interaction and rotational contributions to the energy due to the presence of a vortex that is known to favor the creation of vortices in the condensate once the system is rotating faster than the critical angular frequency, , for energetic stability of a vortex Fetter2001 (). For a rapidly rotating BEC transition (), we can expect that vortices will play a dominant role at all stages in the BEC growth process. It is then possible that atoms condense into a vortex-filled state. In general the vortex configuration can be disordered. It is not required, for example, that the vortices minimize the system energy by forming a regular Abrikosov lattice at all stages of condensate growth.

### iii.1 Properties of the rotating ideal Bose gas and the non-condensate band

An important effect of rotation in Bose-Einstein condensation is the reduction of  Stringari1999 (); Coddington2004a (), that, along with the slow one-dimensional evaporative cooling, accounts for the long evaporation time in the JILA experiments. The transition temperature for the ideal Bose gas is modified by the effective radial frequency in the rotating frame, , so that

 TC=ToC(1−Ω2ω2r)1/3, (24)

where is the nonrotating transition temperature for trapped atoms.

In this work we model the non-condensate region assuming the high-energy atoms are quickly thermalised into a Bose-Einstein distribution, that we approximate as non-interacting since the density is usually small at high energies. We can find analytical expressions for the properties of the noninteracting gas, including the effect of the cutoff, by defining an incomplete Bose-Einstein function

 gν(z,y)≡1Γ(ν)∫∞ydxxν−1∞∑l=1(ze−x)l,=∞∑l=1zllνΓ(ν,yl)Γ(ν), (25)

where is the incomplete Gamma function. In analogy with the reduction to an ordinary Gamma function , we have , reducing to the ordinary Bose-Einstein function. We then find for the non-condensate band density, for example

 nNC(x)=∫RNCd3k(2π)3F(x,k)=1(2π)3∞∑n=1eβ[μ−Veff(x)]∫∞KR(x)dkk2e−nβℏ2k2/2m, (26)

where . Changing integration variable to gives

 nNC(x)=g3/2(eβ[μ−Veff(x)],βℏ2KR(x)2/2m])/λ3dB, (27)

where is the thermal de Broglie wave length. Setting , we recover the standard form for the semi-classical particle density of the ideal gas Dalfovo1999 (). The total number in the non-condensate band is given by

 NNC=g3(eβμ,βER)/(βℏ¯ω)3, (28)

where is the geometric mean frequency in the rotating fame. In this way the usual semi-classical expressions can be generalized to include a cutoff in terms of the incomplete Bose-Einstein function.

The angular momentum and its fluctuations can also be calculated. We now write , and the variance per particle as . In rotating frame coordinates the angular momentum is , which in equilibrium gives the rigid body relation between angular momentum and rotational inertia . From this expression we find the angular momentum is

 ⟨Lz⟩=2ℏΩω⊥g4(eβμ,βER)(βℏω′)4, (29)

where . Similarly, the mean squared angular momentum is

 ⟨L2z⟩=2ℏ2g5(eβμ,βER)(βℏ~ω)5(1+4Ω2ω2⊥), (30)

where . In Figure 3 we give an example of the the variation of , and the angular momentum per particle together with its fluctuations as a function of the angular frequency of rotation. The suppression of by the rotation is seen to be associated with diverging angular momentum and angular momentum fluctuations near the critical frequency for stable rotation in the harmonic trap.

### iii.2 Numerical procedure

The long cooling time and the large changes in length scales in the axial and radial dimensions during cooling makes modeling the full condensation dynamics challenging. However, as the gas cools and spins up, the system becomes progressively more two-dimensional. We model the system by considering the two dimensional projected subspace of the quasi-two dimensional system. We thus use the SGPE of Eq. 21 with rescaled interaction parameter , where in this work we take the thickness through the oblate system as , obtained by projecting onto the ground state of the -dimension. Note that a consequence of the approximation is that the growth rate is unchanged by reduction to the two dimensional effective equation.

We take as our starting point the gas after it has been evaporatively cooled to a rapidly rotating state above the transition temperature. We then simulate the dynamics of a rapid quench below threshold. The procedure is summarized as follows:

• Initial states are generated by sampling the grand canonical ensemble for the ideal Bose gas for initial parameters , and .

• To approximate a rapid quench, the temperature and chemical potential of the non-condensate band are altered to the new values and , while the rotation of the cloud is held fixed at , consistent with axial evaporation. The condensate band field is then evolved according to the SGPE.

Specifically, we aim to model a possible run of the JILA experiment consisting of a system of atoms held in a magnetic trap with frequencies , initially at temperature with chemical potential , and rotating with angular frequency . We examine a range of final quench temperatures with constant chemical potential . We note the ratio , where a value approaching unity indicates the onset of the mean field quantum Hall regime. The ratio , indicates the system is oblate but not highly two-dimensional. We have found that simulations for parameters that are highly oblate are very challenging numerically, so we consider as our first application of the method these more modest parameters. For the results presented here we use a cut-off value of , giving a condensate band containing single particle states. To check the cut-off independence of our results we have also carried out simulations for (not shown here). We find qualitative agreement for single trajectory dynamics and quantitative agreement (at the level of a few percent) for the ensemble averaged condensate fraction shown in Figure 10. We are thus confident that our results are cut-off independent. For computational reasons we choose constant dimensionless damping for all simulations. Although we have derived a form for the damping rate given in Appendix A, the values of arising from Eq. (44) for our parameters are orders of magnitude smaller then our chosen value, so we have chosen a fixed damping such that the growth rate is comparable to the JILA experimental time of 50s. This affects the specific time-scale for growth, but not the qualitative features of the simulations which merely require that the damping is much slower than any other time-scale of the system. We emphasize that the equilibrium properties are independent of .

Since our numerical simulation method employs a projected eigenbasis of the associated linear Schrödinger problem, the initial Bose-Einstein distribution is easily sampled. The grand canonical Wigner distribution for a non-interacting trapped Bose gas distributed over a set of single particle states with frequencies is

 W({αnl,α∗nl})=−∏nl1π(Nnl+1/2)exp{−|αnl|2Nnl+1/2}, (31)

where , the modes are defined in Eq. 55, and

 Nnl=1e(ℏωnl−μ0)/kBT0−1. (32)

The energy cut-off is denoted by the upper limit of the product notation. The single particle spectrum of the condensate band is restricted by the cutoff to

 ℏωnl≡ℏωr(2n+|l|−Ωl+1)≤ER. (33)

The distribution is sampled as , where are complex Gaussian random variates with moments , and . We simulate the time dynamics of the quenched system using the numerical algorithm described in Appendix B.

We note here that application of the Penrose-Onsager criterion for obtaining the BEC Penrose1956 (); Blakie2005a () is hindered in this system because each run of the SGPE theory leads to vortices in different locations, and ensemble averaging the trajectories generates spurious results for . Specifically, the condensate number is significantly reduced below the condensate number that we find using a different technique of short-time averaging single trajectories. Generalizing the Penrose-Onsager criterion for BEC to multi-mode mixed states is not pursued in this work. Here we have taken the view that a given trajectory samples a sub-ensemble of the full range of quantum evolution, consistent with spontaneously broken symmetry in the initial conditions. We present single trajectory results of our simulations of the SGPE that give a qualitative indication of what could be seen in individual experimental measurements of the atomic density time series.

### iii.3 Initial states and dynamics

In Figure 4 we show the density of the condensate band for the thermal gas above the BEC transition. Individual samples of the Wigner function are seen to contain phase singularities. We also construct and diagonalize the one body density matrix, according the the Penrose-Onsager criterion for BEC. By averaging samples from the Wigner function we find that the largest eigenvalue is and the associated eigenvector contains several phase singularities. Using samples washes out the singularities and gives the noninteracting ground state; the point here is that the high degree of degeneracy means that we must average over very many samples to remove the singularities. Put another way, this means that individual trajectories always contain many vortices, and the motion of these vortices destroys long range coherence associated with BEC.

In Figure 5 we plot the condensate number as a function of time after the quench for a range of different quench temperatures. The condensate is calculated by constructing the one body density matrix for the sub-ensemble via short time averaging individual trajectories Blakie2005a (). We have found that by averaging over samples taken uniformly over an interval of radial trap periods, we obtain sufficiently good statistics to determine the condensate number , given by the largest eigenvalue of . The main point to note from this result is that at higher temperatures the absolute condensate number is slightly suppressed. However, as we will see in the next section, the condensate fraction is strongly suppressed, due to the increasing atom number in the non-condensate band for higher temperatures.

A typical time series is shown in Fig. 6, for the lowest temperature considered in Fig. 5, . The initial sample contains many phase singularities, including a few that are circulating counter to the cloud rotation. Counter-circulating vortices only tend to arise in the initial state, are thermodynamically very energetic, and are quickly annihilated by positive vortices during post-quench dynamics. The condensate, shown in Fig. 7, forms with many vortices already present, but in a highly disordered configuration. The vortices then assemble into a regular triangular lattice.

At low temperatures the equilibrium state of the system consists of a an Abrikosov vortex lattice with small thermal fluctuations of the vortices about their average positions. The time dynamics for the highest temperature of Fig. 5, , are shown in Fig. 8. We see that growth proceeds into a disordered state for which periodic order remains frustrated by significant thermal fluctuations.

### iii.4 Equilibrium states

We now consider equilibrium properties of the rotating system, in the grand canonical ensemble generated by evolution according to the SGPE. In Fig. 10 we plot the condensate fraction of the final state for the temperatures used in Fig. 5. The total atom number in the system is calculated as , where is given by Eq. (28) evaluated for our chosen energy cutoff, and is the mean total number of atoms in the condensate band.

We note that the ideal gas relationship between condensate fraction and temperature, , is unaltered by rotation, provided the rotating frame transition temperature, given by Eq. 24, is used in place of the non-rotating transition temperature. The rotating ideal gas curve is plotted for comparison, and we see that the condensate fraction predicted by the SGPE simulations is suppressed relative to the ideal gas result. Although we are unable to simulate the nonrotating system with our two-dimensional algorithm, since the ideal gas behavior is universal with respect to rotation we also compare with the experimental data for the non-rotating Bose gas Ensher1996a (). Our simulations show a condensate fraction significantly suppressed relative to the experimental data for the non-rotating system.

This difference presumably comes about from two effects. The most obvious is the loss of long range order caused by thermal motion of vortices that are absent from the non-rotating system. A less obvious and possibly less significant effect is that the reduction of condensate number by interactions is more pronounced for a rotating gas, due to the exclusion of atoms from vortex cores. Our highest temperature simulation, at , is right on the edge of Bose-Einstein condensation, and this is reflected in the lack of long range order seen in Fig. 8.

The loss of long range order can be seen more clearly by looking at the relative positions of vortices in the condensate wave-function at different temperatures Engels2003a (). In Figure 11 we plot the condensate wave function, along with the histogram of the separation of each pair of vortices in the system. We only consider vortices with radius . For the lowest temperature () the histogram is indistinguishable from the ground state histogram (not shown). The presence of many peaks for large separations reflects the long range spatial order in the system. As temperature increases, order is reduced starting at large vortex separation. At the highest temperature () we see that only two histogram peaks clearly persist, corresponding to nearest and next-to nearest neighbor order being preserved.

## Iv Conclusions

We have developed and applied the SGPE theory to the rotating Bose gas. The formalism presented here can be equally applied to the non-rotating case and in that limit extends the prior SGPE theory SGPEI (); SGPEII (); SGPEIII () by providing expressions for the dissipation rates. We have developed an efficient algorithm for numerically integrating the SGPE in a rotating reference frame and applied it to the case of quenching a rotating Bose gas below threshold, modeling the rapid-quench limit of the experiment at JILA Haljan2001 ().

Our simulations indicate that many vortices co-rotating with the thermal cloud are trapped in the emerging condensate, reflecting the large angular momentum of the system. In our calculations these primordial vortices are nucleated from vortices in the individual samples of the initial Wigner distribution and from fluctuations arising from the condensate growth terms in the SGPE. The process can be regarded as stimulated vortex production, in contrast with the non-rotating case where vortices can appear spontaneously but the total vorticity in the BEC must be zero on the average Drummond1999 ().

For low temperatures the vortices order into a regular Abrikosov lattice as the system approaches thermal equilibrium with the quenched thermal cloud. At sufficiently high temperatures lattice crystallization is frustrated by thermal fluctuations of the positions of vortex cores. This disorder is reflected in a loss of long range correlations in the relative positions of vortices that are most evident as the system approaches .

We note that the concept of rotation has not been investigated in the context of the Kibble-Zurek mechanism, which explains vortex nucleation during the phase transition in terms of spontaneous symmetry breaking and critical slowing down Anglin1999 (). Future quantitative studies of dynamics at the transition will allow for a deeper understanding of the condensation process, in particular the role of spontaneous symmetry breaking as a vortex nucleation mechanism in rotating Bose-Einstein condensation.

We gratefully acknowledge B. P. Anderson and M. K. Olsen for their critical reading of this manuscript, and P. Jain, P. B. Blakie and M. D. Lee for helpful discussions about this work.

This research was supported by the Australian Research Council Center of Excellence for Quantum Atom Optics, by the New Zealand Foundation for Research Science and Technology under the contract NERF-UOOX0703: Quantum Technologies, and by the Marsden Fund under contracts PVT202 and UOO509.

## Appendix A Calculation of dissipation rates

Reservoir interaction rates closely connected to those required here have been previously computed within the context of the quantum kinetic theory (QKT) of BEC QKV (). We now derive the form of the growth and scattering rates including the effects of the cutoff, the external trapping potential and the rotation of the thermal cloud.

### a.1 Growth

#### a.1.1 General form

From Eq. (84) of Ref. SGPEII (), the growth amplitude is

 G(x)≡∫d3vG(+)(x,v,0)=G1(x)+G2(x) (34)

where

 G1(x)≡u2(2π)5ℏ2∫∫∫RNCd3K1d3K2d3K3×F(x,K1)F(x,K2)Δ123(0,0), (35)
 G2(x)≡u2(2π)5ℏ2∫∫∫RNCd3K1d3K2d3K3×F(x,K1)F(x,K2)F(x,K3)Δ123(0,0), (36)

and conserves energy and momentum. Note that effect of working in the rotating frame is only to introduce the effective potential term into the cut-off and Wigner functions in , which is otherwise unaltered. We will take the reservoir to be in equilibrium with Bose-Einstein distribution at chemical potential and temperature . Integration over the momentum -function reduces to

 G1(x)=u2(2π)5ℏ28mπ2ℏ∞∑p,q=1eβ(μ−Veff(x))(p+q)∫∞KR(x)K1dK1×∫∞KR(x)K2dK2e−(pK21+qK22)βℏ2/2m×∫1−1dzδ(z−mVeff(x)/ℏ2K1K2), (37)

where Eq. (11) gives the expression for the momentum cutoff

 KR(x)=√2m[ER−Veff(x)]/ℏ2 (38)

Changing to variables , , and taking account of the delta-function, we obtain

 G1(x)=u2(2π)5ℏ28mπ2ℏ(mβℏ2)2∞∑p,q=1eβ(μ−Veff(x))(p+q)pq×∫∞smin(x)dse−s∫∞tmin(x,s)dte−t, (39)

where , and . We now show that the rate is composed of two regions: an inner region where the rate is spatially invariant, and an outer limit where there is a weak spatial dependence.

#### a.1.2 Inner region

From the from of the lower limits we find that wherever , . Proceeding similarly for we find the same condition to hold, so that in the inner region the growth rates are spatially invariant, taking the final form

 Gin1 = 4mπℏ3(akBT)2[ln(1−eβ(μ−ER))]2, (40) Gin2 = 4mπℏ3(akBT)2e2β(μ−ER) (41) ×∞∑r=1erβ(μ−2ER)(Φ[eβ(μ−ER),1,r+1])2,

where the Lerch transcendent is defined as . These terms are superficially similar to the condensate growth rates previously obtained using the quantum kinetic theory QKPRLII (), although the form of the cutoff dependence is different in the SGPE theory. Note that this form is independent of any collisional contribution to the effective potential.

#### a.1.3 Outer region

In the region where the rates depend on position, albeit quite weakly. At the edge of the condensate band where , we find

 Gout1 = 4mπℏ3(akBT)2βER (42) ×∞∑p,q=1eβ(μ−ER)(p+q)√pqK1(β√pq) Gout2 = 4mπℏ3(akBT)2βER, (43) ×∞∑p,q,r=1eβ(μ−ER)(p+q)+β(μ−2ER)r√(p+r)(q+r) ×K1(β√(p+r)(q+r)),

where is the modified Bessel function of the second kind.

In Fig. 12 we compute the variation of with radius for a spherical parabolic potential. For illustration we have also neglected the collisional contribution to the effective potential. To show the maximum variation of the total rate we also plot (inset) the ratio versus . For all data we use which is typically a good approximation for the energy at which the spectrum of the many body system becomes single particle-like SGPEII (). The rates are spatially invariant except for a small region near the edge of the condensate band where the scattering phase space in the harmonic trap is modified by the energy cutoff. Since the variation with radius is usually not significant, for our simulations we will use the approximate growth rate

 G(x)≈γ≡Gin1+Gin2, (44)

which is exact in the inner region and a good approximation elsewhere.

### a.2 Scattering

Our starting point of the scattering amplitude calculation, that is related to the actual scattering via Fourier transform, is Eq. (86) of Ref. SGPEII ()

 ~M(x,k,ϵ)=4πu2ℏ2∫RNCd3K1(2π)3∫RNCd3K2(2π)3Δ12(k,ϵ)×F(x,K1)[1+F(x,K2)], (45)

where . Upon evaluating the momentum delta-function we find

 ~M(x,k,0)=4mu2(2π)5ℏ3∫RNCd3KF(x,K)[1+F(x,K)] ×δ(2{K−mℏΩ×x}⋅k−k2). (46)

Using the notation , where

 eK(x)≡ℏ2K22m+Veff(x), (47)

in the variable we have

 ~M(x,k,0) = 4mu2(2π)5ℏ3∫RNCd3qδ(2q⋅k−k2), (48) ×fT,μ(eq(x))[1+fT,μ(eq(x))],

where . Writing to distinguish the term with one or two BE-distributions as before, and choosing the -component of along , we find, for the term involving one BE-distribution

 ~M1(x,k)=4mu2(2π)5ℏ3π|k|∞∑p=1epβ(μ−Veff(x))×∫∞KR(x)K1dK1e−pβℏ2K21/2mΘ(|k|/2K1≤1). (49)

The theta function generates the effective cutoff . Evaluating the integral gives

 ~M1(x,k)=4mu2(2π)5ℏ2πmβℏ2|k|×eβ(ℏ2K2min(x)/2m+Veff(x)−μ)(eβ(ℏ2K2min(x)/2m+Veff(x)−μ)−1)2, (50)

which can be further simplified as follows. We want the condition for , which can be expressed as . The restriction of is determined by the requirement that , where