Rapidly rotating Bose-Einstein condensates in anharmonic potentials

Rapidly rotating Bose-Einstein condensates confined in anharmonic traps can exhibit a rich variety of vortex phases, including a vortex lattice, a vortex lattice with a hole, and a giant vortex. Using an augmented Thomas-Fermi variational approach to determine the ground state of the condensate in the rotating frame - valid for sufficiently strongly interacting condensates - we determine the transitions between these three phases for a quadratic-plus-quartic confining potential. Combining the present results with previous numerical simulations of small rotating condensates in such anharmonic potentials, we delineate the general structure of the zero-temperature phase diagram.


Introduction
Gaseous Bose-Einstein condensates of alkali-metal atoms are fertile systems in which to study the behaviour of superfluids under rotation. For rotation frequencies above c1 , the critical frequency for creating a single vortex, the properties of these systems resemble those of liquid helium under typical experimental conditions. States with a small number of vortices were first observed in the experiments of [1,2]. As increases, more and more vortices form, and indeed it has become possible to create lattices with up to ∼ 200 vortices [3]- [5]. The physics of rapidly rotating gases is different in harmonic and anharmonic traps. In a harmonic trap, the rotation rate is limited by the radial trap frequency ω, which sets the scale of c1 . As the angular frequency of rotation of the gas approaches the transverse trapping frequency, the centrifugal force approaches the restoring force exerted by the trap and the atoms become more and more weakly contained. On the other hand, the regime > ω becomes accessible in anharmonic trapping potentials that rise more steeply than quadratically.
Condensates rotating in anharmonic trapping potentials are expected to have a rich structure of vortex states, as described in [6]- [9], reviewed below. These states include an array of singly quantized vortices (which becomes a triangular lattice when the number of vortices is large); a vortex array with a finite radius 'hole' in the centre-opened by the centrifugal force-in which the condensate density is zero but the vorticity non-zero; and giant (multiply quantized) vortices 4 . Our aims in this paper are twofold. First we study, via a variational approach, the states in an anharmonic quadratic-plus-quartic radial potential [11], where ρ is the (cylindrical) radial coordinate, m is the atomic mass, d = (1/mω) 1/2 is the oscillator length (h = 1), and λ is a dimensionless measure of the anharmonicity. We consider a container infinitely long and cylindrically symmetric and rotating about the z axis, with a density per unit length σ = N/Z in the axial (z) direction; the particles are assumed to interact via s-wave scattering with scattering length a s > 0. The variational approach is adopted from [9], with the density profile of the gas given by an augmented Thomas-Fermi approximation which explicitly takes into account the important kinetic energy of flow about the individual vortices and the modification of the interaction energy by the density variations due to the vortices. We then combine the present calculation with the results of the earlier studies to delineate the phase diagram of a rotating Bose condensate in an anharmonic trap, when the total vorticity is 1 (i.e. well above the critical rotation frequency c1 ). As we shall see, the phase diagram that consistently accounts for these results has the schematic structure in the -σ a s plane shown in figure 1. At very small interaction strength or density, multiply quantized vortices are always preferable, as found in [6], while at larger interaction strength at small the system forms an array of singly quantized vortices. At larger the array, as shown in [8], develops a hole in the centre. The phase diagram contains a triple point between the three distinct regions. We infer the existence of a triple point from the fact that in the multiply quantized vortex regime the order parameter is topologically different from that 51.3  in the vortex lattice regime, while the vortex lattice with hole is distinguished from the vortex lattice by a vanishing average density near the rotation axis.

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
Lundh [6], by numerical solution of the Gross-Pitaevskii equation in two dimensions, determined the ground state vortex structures of small total vorticity 10, in a rapidly rotating condensate in radial trapping potentials ∼ρ n (n > 2), as well as in the potential (1). His calculations, which are in the parameter range corresponding to the lower left of figure 1, show that for small dimensionless interaction strength, σ a s , the system prefers multiply quantized vortices, as illustrated in figure 1, while above a critical (λ-dependent) interaction strength, (σ a s ) c ∼ 1, the ground state consists of an array of singly quantized vortices. The fact that a giant vortex state is favoured at small σ a s can be understood from the fact that in the limit of a non-interacting gas trapped in a potential that rises faster than quadratically, the dependence of the single-particle energies, (ν), on angular momentum,hν, rises faster than linearly [6,12]. Thus, in this limit, in the ground state of a fixed number of bosons in the frame rotating at given , all the particles occupy the single angular momentum state that minimizes (ν) − ν. Condensates in the potential (1) were also studied via numerical simulation by Kasamatsu et al [8], at fixed parameters λ = 1/2 (corresponding to their k = 1) and σ a s = 250/8π. At low the system contains an array of singly quantized vortices. The value of σ a s assumed in this calculation well exceeds the critical coupling (σ a s ) c found by Lundh, and is above the triple point in figure 1. As they increase the rotation frequency to 2.14ω, the vortex array begins to develop, as a consequence of the centrifugal force, a hole in the centre. At still higher rotational velocity, ∼ 3.2-3.5, the system undergoes a transition to a single multiply quantized giant vortex. Simula et al [13] have also done numerical simulations of stable multiply quantized vortices in the presence of a pinning potential.
Similar behaviour was found in [9] for a condensate contained in a cylindrical hard-walled bucket, by constructing an augmented Thomas-Fermi solution of the Gross-Pitaevskii equation. At low rotation speed the vortices form a uniform lattice. However, at a critical rotation frequency, h ∼ 1/2mξ 0 R-where ξ 0 = (8πna s ) −1/2 is the coherence or healing length, andn is the average density-the lattice begins to develop a hole in the centre. Finally, at a second critical frequency,

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
G ∼ 1/2mξ 2 0 , the system makes a transition to a giant vortex state. Fetter [7], in considering condensates in the potential (1), finds a critical rotation speed at which the density in the centre of the trap goes to zero; however, the onset of the hole is obscured by the particular approximation employed following factorization of the condensate wavefunction.
While detailed comparison of our calculation with the numerical simulations of [6] and [8] is hindered by the fact that the number of vortices in the simulations is not large enough to form a uniform lattice, as we assume (see, e.g., figure 1(e) of [8] and figure 6 of [6]), our results are qualitatively consistent with and provide a simple theoretical understanding of the numerical results. At large σ a s , where our approach is valid, we see that with increasing (> ω) the system indeed undergoes a transition from a vortex lattice to a vortex lattice with a hole around the axis of rotation. Furthermore, our calculations show that as increases the system eventually undergoes a transition to a single giant vortex, for σ a s up to at least 350, the largest value for which we studied this transition.
We first describe the three states of the rapidly rotating gas-a vortex lattice, both with and without a hole along the axis, in section 2.1, and a giant vortex, in section 2.2-and calculate the augmented Thomas-Fermi solutions for the lowest energy in the rotating frame at given . In section 3 we determine the transitions between different phases for the particular anharmonicity λ = 1/2, and the less steep λ = 1/8. In section 3 we summarize our results, and examine the limits of validity of our approach in section 4.

Vortex lattice
We first consider a uniform lattice of singly quantized vortices in which the average density varies slowly on the scale of the inter-vortex spacing. In the Thomas-Fermi regime, in which U 0 n , where U 0 = 4πh 2 a s /m, one can approach the vortex state by integrating out the condensate structure on small length scales, in the sense of the renormalization group, and including its effects in an effective long-wavelength Hamiltonian [9]. We write the order parameter as ψ(r) = s(r) √ n(r), which is the product of a rapidly varying factor, s(r) = e iϕ(r) |s(r)|, which vanishes at each vortex core, and has a phase ϕ which wraps by 2π around each vortex, times a real slowly varying envelope function, √ n(r). If we choose |s| 2 to average to unity over each unit cell of the lattice, then n(r) is the smoothed density profile of the system. The kinetic energy in the laboratory frame can be written in terms of s and n as Integrating the last term by parts, and using the fact that |s| 2 averages to unity over each unit cell of the lattice, we have The We write the final term as a sum over integrations about individual vortices (labelled by i), and replace the unit cells by cylindrical Wigner-Seitz cells of radius , with π 2 = n v . In the neighbourhood of a given vortex, δv 1/mρ c , where ρ c is the radial coordinate measured from the vortex line and |s| is approximately radially symmetric about the vortex line. Thus where the integration is over the unit cell containing the given vortex i. Since n varies slowly over a given cell, we may write the contribution of a given cell to the energy as d 3 r n(r)n v A i , where The A i are effectively independent of the cell, and in the following we shall assume the A i to have a uniform value which we write as πa/m, where (recalling with the integration over the central cell. Neglecting the slowly varying gradient of the smoothed density, we may write the total kinetic energy as The interaction energy (U 0 /2) d 3 r |ψ(r)| 4 is greater than (U 0 /2) d 3 r n(r) 2 , where n is the smoothed density profile. The difference is given by a multiplicative renormalization of U 0 by the factor b = n 2 / n 2 , which is of order unity [9]. In addition, the total angular momentum along the z axis is given by [15] L z = d 3 r n(r)( v mρ 2 + 1); (9) the two terms are the angular momentum of solid body rotation, and the contributions of the angular momenta of the individual cells. Thus the energy of the lattice in the rotating frame, E L = E − L z , in the Thomas-Fermi regime is We denote the smoothed density of the lattice array, with cylindrical symmetry about the trap axis, by n L (ρ).
To evaluate a and b we assume the simple variational ansatz of [9], that the order parameter ramps up linearly over a distance ξ from the cores of the vortices, and then flattens out, where ζ = (ξ/ ) 2 is the fractional area of the vortex core in the unit cell. The effective core radius ξ is in general different from the Gross-Pitaevskii healing length ξ 0 . With equation (11),

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
Minimization of E L with respect to n L (ρ) yields the Thomas-Fermi profile, with µ the chemical potential. Since the density of the gas at the centre of the cloud, ρ = 0, cannot be negative, the system begins to develop a hole when n L (ρ = 0) falls to zero, or from equation (13), µ−(a +1) v + = 0. For smaller , this quantity is positive and the inner radius of the cloud, R 1 , vanishes, while for larger , R 1 = 0. We note that when a hole develops in a uniform vortex lattice, the hole must, for the lattice to remain uniform, contain the same vorticity, πR 2 1 n v = mR 2 1 v (where n v is the number of vortices per unit area) as would be present were the hole filled with condensate. Equation (13) implies that the inner and outer radii of the cloud, R 1 and R 2 , respectively, where the density vanishes, are given by the solutions of where here lengths with a tilde are measured in units of d, e.g.R i = R i /d, i = 1, 2, and energies or frequencies with a tilde are measured in units of ω.
Integrating the density over the volume, n L d 3 r = N, we find Minimization of E L with respect to v yields v = − N(a + 1)/I, where I = d 3 r mρ 2 n L (ρ) is the moment of inertia. Equation (16), with (13), implies that For a given value of the dimensionless parameters σ a s , λ and ζ , we solve equations (14), (15) and (17) to find R 1 , R 2 , µ and v ; then from equations (10) and (13) we construct the energy E L ( , ζ ) in the rotating frame, as a function of and ζ , and finally vary ζ to minimize E L . A useful check on the numerical results is the relation ∂E L /∂ = −L z , implied by equation (10). We give the results in section 3 after discussing the description of the giant vortex state.

Giant vortex
The order parameter of the giant vortex state, whose vorticity is localized along the rotation axis, is where n G (ρ) is the density profile of the cloud, ρ and φ are cylindrical polar coordinates, and ν is the winding number. The angular momentum per particle along the z axis is L z /N = ν, while the energy of the system in the rotating frame in the Thomas-Fermi limit is

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
Setting the variational derivative δE G /δn G (ρ) equal to the chemical potential, µ, we find the density profile of the cloud: where µ = µ + ν . The inner and outer radii of the cloud, R 1 and R 2 , solve In addition, the constraint n G d 3 r = N implies that Minimization of E G with respect to ν yields and thus 4σ a s˜ = ν 2μ lnR For given values of the dimensionless parameters σ a s and λ, we solve equations (21), (22) and (24), and express R 1 , R 2 , µ and in terms of ν. We then numerically solve for E G from equations (19) and (20) as a function of ν, and thus implicitly as a function of . The relation ∂E G /∂ = −L z , implied by equation (19), again provides a check on the numerical results.

Results
We have limited the calculations of the vortex phases here to the cases of intermediate anharmonicity, λ = 1/2, and weak anharmonicity, λ = 1/8, which well illustrate the richness of the phase diagram. Figure 2 shows the phase diagram for λ = 1/2. Relatively small and relatively large σ a s favour the formation of a vortex lattice. In the opposite limit, large and small σ a s favour the formation of a giant vortex. Finally, at intermediate (for σ a s in the range considered) the system prefers a vortex lattice with a hole. We do not see evidence in our calculations for the triple point shown in figure 1, since comparison of the numerical results of [6] and [8] indicates that it must occur at a small interaction strengths where the present approach fails. However, the existence of the triple point is required by the fact that in the multiply quantized vortex regime the order parameter is topologically different from that in the vortex lattice regime, while the vortex lattice with hole is distinguished from the vortex lattice by a vanishing average density near the axis. Quite generally, the transition from a vortex lattice to a vortex lattice with a hole is continuous, second order, while the transition from a vortex lattice with or without a hole to a giant vortex is discontinuous, first order. For σ a s = 200, the vortex lattice develops a hole for /ω ≈ 3.90, and then makes a transition to a giant vortex for /ω ≈ 5.48. At this point the giant vortex has a vorticity ν 168. For σ a s = 550, corresponding to the parameters of the multivortex MIT experiment [3], we find that the vortex lattice develops a hole for /ω ≈ 4.50, and that even up to the highest value of /ω ≈ 8.4 that we examined, this solution has lower energy than the giant vortex. For such high values of σ a s and we observe that the two solutions (vortex lattice with a hole and giant vortex) look macroscopically very similar. The question remains open whether for all large interaction strengths a giant vortex state will eventually have lower energy than a vortex lattice with a hole.
For comparison, we have also calculated the phases for interaction strength σ a s = 250/8π ≈ 10 considered in [8], although it is well below the limits of validity of the present approach. Then the transition to a giant vortex occurs for /ω ≈ 1.5, compared with the value ∼3.2 found in [8]; however ν (taken to be continuous in the calculation) is between 2 and 3. For such low values of ν our approach is clearly not valid, and the discrepancy should not be taken seriously.
The particular case of weak anharmonicity λ = 1/8, with σ a s = 90 also shows the three phases: for /ω ≈ 2.30 the vortex lattice develops a hole in the middle, while for /ω ≈ 2.40 it makes a transition to a giant vortex. Figure 3 shows the gentle crossing of E G ( ) (solid curve), and E L ( ) (dashed curve), while figure 4 shows the angular momentum per particle L z /N (of the lowest state), the negative of the slope of the lower of the two curves of figure 3, as a function of . Since the curves cross at a finite but small angle, L z /N is discontinuous at the point of the transition. Figure 5 shows R 1 and R 2 as a function of for the three solutions.

Limits of validity of the present approach
The present approach makes two crucial assumptions, first that the profile of the gas can be described by the Thomas-Fermi approximation, with the vortex corrections to the kinetic and interaction energies included, and second that the array of singly quantized vortices is a uniform lattice. The validity of the Thomas-Fermi approach requires that the kinetic energy per particle, ∼1/2m(R 2 − R 1 ) 2 , associated with localizing the atoms within the radii of the profile of the cloud R 1 and R 2 , be much less thannU 0 = 1/2mξ 2 0 , wheren ∼ N/π(R 2 2 − R 2 1 )Z is the average   density of the cloud, and ξ 0 is the coherence length. Thus, validity of Thomas-Fermi requires or equivalently, ξ 0 R 2 − R 1 . For R 1 = 0, the condition (25) is simply 8σ a s 1. In addition to condition (25), the basic features of the lowest state have to be captured by the trial functions being considered. For small values of σ a s and /ω, the numerical solutions [6,8] have too small a vorticity for the vortices to arrange in a triangular lattice. Rather, in [8] at /ω = 2.5 for example, the vortex hole, containing four units of vorticity, is surrounded by a ring of 10 vortices, instead of a developed lattice. A reasonable estimate for the lowest value of N v at which the calculation is valid is ∼50, which typically (e.g. at λ = 1/2) requires /ω 3.5.

Summary
In this study we have examined, by means of a variational method, the lowest state of a Bose-Einstein condensate rotating in a quadratic-plus-quartic confining potential, in the limit of fast rotation. We find that as the frequency of rotation, the density of the gas and the strength of the quartic term in the trapping potential are varied, the system undergoes transitions among three phases: a vortex lattice with and without a hole, and a giant vortex. The phase diagram, whose structure we have mapped for a specific choice of the coefficient of the anharmonic term, is quite rich. This structure is not limited to the specific form of the anharmonic potential considered here. The augmented Thomas-Fermi approach is reliable in the limit of large densities and high angular frequencies. Many problems remain, including determining the details of the critical point in figure 1, which requires going beyond Thomas-Fermi; the variation of the phase diagram as the coefficient of the quartic term in the confining potential is varied; and finally, the structure of the phase diagram at finite temperatures.