Vortices in rotating Bose-Einstein condensates confined in homogeneous traps
aa r X i v : . [ m a t h - ph ] A p r Vortices in rotating Bose-Einstein condensatesconfined in homogeneous traps
T.Rindler-Daller
Fakult¨at f¨ur Physik, Universit¨at Wien, Boltzmanngasse 5, 1090 Vienna, Austria
Abstract
We investigate analytically the thermodynamical stability of vortices in the groundstate of rotating 2-dimensional Bose-Einstein condensates confined in asymptoti-cally homogeneous trapping potentials in the Thomas-Fermi regime. Our startingpoint is the Gross-Pitaevskii energy functional in the rotating frame. By estimatinglower and upper bounds for this energy, we show that the leading order in energyand density can be described by the corresponding Thomas-Fermi quantities andwe derive the next order contributions due to vortices. As an application, we con-sider a general potential of the form V ( x, y ) = ( x + λ y ) s/ with slope s ∈ [2 , ∞ )and anisotropy λ ∈ (0 ,
1] which includes the harmonic ( s = 2) and ’flat’ ( s → ∞ )trap, respectively. For this potential, we derive the critical angular velocities for theexistence of vortices and show that all vortices are single-quantized. Moreover, wederive relations which determine the distribution of the vortices in the condensatei.e. the vortex pattern. Key words:
Static properties of condensates; thermodynamical, statistical and structuralproperties, Tunneling, Josephson effect, Bose-Einstein condensates in periodicpotentials, solitons, vortices and topological excitations
PACS: : 03.75.Hh, 03.75.Lm
Many efforts have been made in understanding ultra-cold quantum gases,especially since the experimental achievement of Bose-Einstein condensates(BECs) in 1995. A particular interesting subject is the study of rotating BECs.When the trap is subjected to an external rotation the condensate does not ⋆ Journal Ref.: Physica A 387, 1851-1874 (2008)
Preprint submitted to Elsevier otate like a solid body. Instead, beyond a critical angular velocity quantizedvortices appear manifesting the genuine quantum character of the system. In-deed, vortices in BECs were observed in 1999 for the first time (see Refs. [29]and [27,28]). Theoretical studies were already presented before (see e.g. [32]for one of the earliest papers on the subject) and have since then grown toa substantial branch of its own (see e.g. [3,8,11,12,15,16,17,25,36]). A generaltreatment of BECs can be found in the textbooks of [30] and [31].Most of the theoretical studies have been undertaken in the framework ofthe Gross-Pitaevskii (GP) theory whose validity as an approximation of thequantum mechanical many-body ground state was established in [22] for thenon-rotating case and in [24] for rotating systems. Particular attention hasbeen put on the so-called Thomas-Fermi (TF) regime of strong coupling. Thisis especially true for the study of vortex structures (see the monograph [1]). In[18,19] a rigorous analysis of vortices for BECs in harmonic anisotropic trappotentials was achieved for a GP-type functional in the TF limit. A previousanalysis was developed in [35] in the context of superfluids. The methodol-ogy of those papers originates from [7] where a rigorous analysis of vorticesin Ginzburg-Landau models of vanishing magnetic field in the regime whichcorresponds to the TF limit was developed. In [33] and [34], general resultson symmetry breaking which are not limited to the TF regime were proven intraps of arbitrary shape.Within the GP theory, the properties of vortices are determined by two phys-ical parameters apart from the external trap, namely angular velocity andinteraction strength between the particles. In this paper, we consider theground state of rotating 2D Bose-Einstein condensates which are trapped inasymptotically homogeneous anisotropic potentials rotating with angular ve-locity Ω. The aim of this investigation consists of deducing analytically theGross-Pitaevskii energy and density in presence of vortices and deriving theirproperties in the TF regime. We consider thermodynamical conditions forvortex existence, i.e. we are looking for angular velocities which reduce thetotal energy in such a way that vortices are energetically favoured to appear.This work was originally inspired by the papers of [4] and [8] which consideranisotropic harmonic potentials. There and for instance in Refs. [15,25,32], itwas established by numerical methods that vortices are single-quantized. Weshow here by analytical estimates, in particular, that this is true for a verylarge class of trapping potentials. In fact, the majority of studies uses nu-merical and variational methods for a limited number of trap potentials (e.g.harmonic or harmonic-plus-quartic) whereas we derive analytical formulae fora very large class of potentials. Thereby, we try to present the analysis in sucha way that both the physical ideas and mathematical estimates are broughtout in a clear way.This paper is organized as follows: In
Section 2 , we state the setting andpresent the main result. We decompose the condensate wave function in avortex-free part and a vortex-carrying part. This allows a splitting of the un-2erlying energy functional in separate contributions which can be estimatedsubsequently. In
Section 3 , we study the leading asymptotics of the energyand density. In
Section 4 , we justify a model for the structure and number ofvortex cores which is compatible with the considered order of magnitude ofthe angular velocities.
Sections 5-7 contain lower and upper bound estimatesof the vortex-carrying energy contributions in terms of the winding numberof the vortices and the coupling parameter. In
Section 8 , we specify an ex-ternal potential which is of a general anisotropic homogeneous form. For thispotential, we deduce the critical angular velocity for the appearance of one ora finite number of vortices. The leading orders of the energy in presence ofvortices are calculated and it is shown that all vortices have winding numberone, i.e. they are all single-quantized. Furthermore, we deduce relations whichdetermine the distribution of the vortices in the condensate, i.e. the vortexpattern. Finally, in
Section 9 we present the conclusions.
Our starting point is the 2D Gross-Pitaevskii energy functional in the referenceframe rotating (uniformly) with ˜ Ω (see e.g. [8,30,31]): E GP [ u ] = Z R " ~ m |∇ u | + V ( r ) | u | + N g | u | − iu ∗ ~ ˜ Ω · ( ∇ u × r ) . (1)Indeed, it is only meaningful to consider this reference frame as far as the(temporal) stability of structures is concerned which appear due to the ro-tation: The external trap is time-independent and the states are stationarywith respect to that frame (see e.g. Ref. [8]). The function u ( r ) is a complexfield (the complex conjugate is denoted as u ∗ ). We write the associated po-lar decomposition as u = | u | e iS u where | u | is proportional to the density ofcondensed particles with normalization Z R | u | = 1and S u is the phase function. The minimizer of (1) is called the order parameteror ’wave function of the condensate’ in the rotating frame. The external trappotential is denoted by V and N is the number of particles with mass m .The third term in (1) describes the effective interaction between the particleswhere the coupling constant in 2D is given by g = √ π ~ a/ ( mh ) with the3D scattering length a and the thickness h of the system in the stronglyconfined direction which we choose to be the z -axis, so that the system iseffectively 2D (in the x - y -plane). We denote × as the vector product in R , r = ( x, y,
0) and ˜ Ω = (0 , , ˜Ω) is the angular velocity vector assuming thatthe gas rotates around the z -axis. An important parameter, consisting of the3cattering length and a density, is given by the ’healing length’ ξ . It is definedoriginally by setting ~ / (2 mξ ) = 2 πρa ~ /m where the r.h.s is the energy perparticle for gases in a box in the limit of dilute systems ρa → ρ , so ξ = 1 / √ πaρ . For inhomogeneous and rotating systems, the healinglength may be defined accordingly by using an appropriate (mean) value forthe density. In particular, the healing length determines the effective radius ofa vortex core in rotating systems.In 2 dimensions, the ratio between the healing length and the characteristiclength of the system L , which is set by the external trap or box respectively,is ε ∼ ξ L = ~ √ πN gm where we introduce the dimensionless parameter ε . In this paper, we will beconcerned with the TF limit where this ratio tends to zero (meaning physicallythat 0 < ξ ≪ L or 0 < ε ≪ ε -dependent factor (see also [9]):Suppose V is homogeneous of order s , i.e. V ( γ r ) = γ s V ( r ) for γ >
0. Werescale the energy functional (1) by setting r = k r ′ and u ( r ) = u ′ ( r ′ ) /k with N g/ ~ / (4 ε m ) and k = ( ~ / (4 ε m )) / ( s +2) . Then we have E GP [ u ] = 1 k Z R " ~ m |∇ ′ u ′ | + k s +2 V ( r ′ ) | u ′ | + N g | u ′ | − i ~ u ′∗ k ˜ Ω · ( ∇ ′ u ′ × r ′ ) d r ′ (2)with R | u ′ | = 1. Choosing ~ = 1 = m and inserting k , (2) becomes E GP [ u ] =(16 ε ) / ( s +2) E GP ′ [ u ′ ] with the energy on the r.h.s. (omitting the primes) E GP [ u ] = Z R " |∇ u | + | u | ε ( V + | u | ) − iu ∗ Ω ( ε ) · ( ∇ u × r ) (3)and the scaled angular velocity Ω ( ε ) is related to the original unscaled one byΩ( ε ) = ˜Ω / (16 ε ) / ( s +2) . (4)For brevity, we will also write Ω but it should be kept in mind that Ω dependson ε after scaling. In the forthcoming, we study the functional in (3) whichcan be also written in the following form E GP [ u ] = Z R " | ( ∇ − i ( Ω × r )) u | + | u | ε ( V + | u | ) −
12 Ω r | u | (5)and r := | r | . Critical points of E GP [ u ] are solutions of the following associatedEuler-Lagrange equation, called Gross-Pitaevskii equation∆ u = u ε ( V + 2 | u | − ε µ GP ) + 2 i ( Ω × r ) · ∇ u (6)4here the GP chemical potential µ GP is fixed by the normalization. Denotinga minimizer of (3) as u ε , it is given by µ GP = E GP [ u ε ] + 14 ε Z R | u ε | . (7)The corresponding amplitude squared | u ε | will be referred to as Gross-Pitaevskiidensity. Inserting u = | u | e iS u into (6) results in hydrodynamic-like relationsfor the density and the velocity:∆ | u | − | u | ( ∇ S u ) + 2 | u | ( Ω × r ) · ∇ S u − | u | ε ( V + 2 | u | − ε µ GP ) = 0 , ∇ · [ | u | ( ∇ S u − Ω × r )] = 0 . The GP functional (3) for Ω = 0 decribes the gas without rotation E GP [ f ] = Z R "
12 ( ∇ f ) + f ε ( V + f ) (8)with f a real, positive function. The minimizer of (8) will be denoted as f ε .The normalization condition R R f = 1 fixes the associated chemical potential ν GP which is given by ν GP = E GP [ f ε ] + 14 ε Z f ε (9)and which is of the order 1 /ε . The functional (8) tends for ε → E TF [ ρ ] = 14 ε Z R ρ ( V + ρ ) , (10)which is a functional for the density ρ = f alone. It can be shown (see e.g.Ref. [22]) that it has a unique positive minimizer, the Thomas-Fermi density, ρ TF = 12 [4 ε µ TF − V ] + =: 12 [ µ − V ] + (11)where [ . ] + denotes the positive part and µ := 4 ε µ TF . The TF chemical po-tential µ TF (or µ respectively) is determined by Z D ρ TF = 1 (12)where D = { ( x, y ) ∈ R : ρ TF > } is the Thomas-Fermi domain whose shape depends on the external potential V . Moreover, µ TF = E TF [ ρ TF ] + 14 ε Z ( ρ TF ) , /ε whereas µ := 4 ε µ TF is of the order of a constantindependent of ε . In the TF regime where ε is small, vortex cores are small compared to thecharacteristic length scale of the system, producing narrow ’holes’ which ef-fectively shrink as ε →
0. It is argued in
Section 4 that vortices appear at acritical angular velocity of the order Ω ≃ C | ln ε | with C a positive constant(independent of ε ) depending on the external trap. Explicit expressions for C will be determined in the forthcoming analysis (see also Refs. [4] and [8] forthe harmonic trap case).In the minimization of (8), i.e. (3) with Ω = 0, one considers all functions inthe subspace of angular momentum zero and the density profile is given by f ε .Considering (3) with Ω > ≤ C | ln ε | asymptot-ically, the overall density can still be described by the vortex-free density f ε ingood approximation. However, in a non-isotropic potential V there appears aphase S (depending on V ), i.e. the vortex-free function is then more generally f e iS . Since this function has no vortex, the phase S is non-singular and (6)gives ∆ f = − f ∇ S · [2( Ω × r ) − ∇ S ] + f ε ( V + 2 f − ε ˜ ν GP ) (13)and ∇ · [ f ( ∇ S − Ω × r )] = 0 (14)where ˜ ν GP is the associated chemical potential. A solution without vortex isa minimizer of the problem min {E GP [ f e iS ] : f e iS ∈ H with f > , R f = 1 } (see also [18]) with E GP [ f e iS ] = Z R "
12 ( ∇ f ) + f ε ( V + f ) + 12 f [( ∇ S ) − ∇ S · ( Ω × r )] . (15)Later on in this paper we are going to consider external traps of the form V ( x, y ) = ( x + λ y ) s/ (16)with slope s ∈ [2 , ∞ ) and λ ∈ (0 ,
1] describing the anisotropy. It is a fairly gen-eral potential which includes also the important special cases of the harmonic( s = 2) and flat ( s → ∞ ) trap which are extensively used in experiments. Thecorresponding phase to this potential is S = λ − λ + 1 Ω xy (17)which vanishes for the isotropic case λ = 1. This expression for S was alsodeduced for the harmonic trap in Refs. [4] and [8]. Note, however, that it is6ot dependent on the slope parameter s . We also see that the terms in (15)involving ∇ S are at most of the order | ln ε | for Ω ≤ | ln ε | and hence of muchlower order than the remaining part described by (8) which is ∼ /ε .We now decompose the order parameter u of (3) into the vortex-free part f e iS and a part which carries the vorticity. A similar splitting can be found in Refs.[4,5,8,20] and more recently in Refs. [6,14,37]. Writing u = | u | e iS u = f e iS v = f | v | e i ( S + S v ) with | u | = f | v | and S u = S + S v , the contribution v = | v | e iS v accounts for the presence of vortices. In a vortex point, the amplitude vanishes,i.e. | u | = | v | = 0 since f = 0 and the phase fulfills the usual circulationcondition which is a quantization condition because u (resp. v ) is a complexfield: I C ∇ S u · τ = I C ( ∇ S v + ∇ S ) · τ = 2 πd + 0since S has no singularity and τ is a unit tangent vector to the curve C encircling the vortex with winding number d . Without the presence of vortices,there would be u = f e iS with density | u | = f and the phase S u would benon-singular. Inserting the decomposition u = f e iS v in the energy functional(3) results in the following splitting (see also [4,5]) where the integrals are over R : The first term becomes Z |∇ ( f e iS v ) | = Z (cid:20) f |∇ v | + 12 | v | [( ∇ f ) + f ( ∇ S ) ] + 14 ∇ ( f ) · ∇| v | + 12 f ∇ S · ( iv ∇ v ∗ − iv ∗ ∇ v ) (cid:21) , the second one is simply Z | f e iS v | ε ( V + | f e iS v | ) = Z f | v | ε ( V + f | v | )and for the rotation term we get − Z if e − iS v ∗ Ω · ( ∇ ( f e iS v ) × r ) = Z if v ∗ ∇ v · ( Ω × r ) − Z f | v | ∇ S · ( Ω × r ) . Putting the terms together and separating the vortex-free part of the energy(15), we have E GP [ u ] = E GP [ f e iS ] + Z ( | v | − × "
12 ( ∇ f ) + 12 f ( ∇ S ) − f ∇ S · ( Ω × r ) + V f ε ++ Z ∇ ( f ) · ∇| v | + Z f |∇ v | + Z f ε ( | v | −
1) ++
Z (cid:20) f ∇ S ( iv ∇ v ∗ − iv ∗ ∇ v ) + if v ∗ ∇ v · ( Ω × r ) (cid:21) . (18)7he third term of this expression becomes Z ∇ ( f ) · ∇| v | = Z ∇ ( f ) · ∇ ( | v | −
1) = − Z
14 ( | v | − f )= − Z
12 ( | v | − f ∆ f − Z
12 ( | v | − ∇ f ) = Z ( | v | − " f ∇ S · ( Ω × r ) − f ( ∇ S ) − f ε ( V + 2 f − ε ˜ ν GP ) −
12 ( ∇ f ) , where we used (13) for ∆ f .Moreover, for the fifth term in (18) we use the identity Z f ε ( | v | −
1) = Z f ε ( | v | −
1) + Z f ε (1 − | v | ) . Inserting the last two equations into (18) we get the following splitting of thefunctional in (3) E GP [ u ] = E GP [ f e iS ]+ Z R " f |∇ v | + f ε (1 − | v | ) − Z R if v ∗ ∇ v · ( ∇ S − Ω × r )=: E GP [ f e iS ] + G f [ v ] − R f [ v ] (19)where we used ˜ ν GP R f ( | v | −
1) = 0 because of the normalization conditionsand the last term in (18) was written in a more convenient form using
Z (cid:20) f ∇ S · ( iv ∇ v ∗ − iv ∗ ∇ v ) + f iv ∗ ∇ v · ( Ω × r ) (cid:21) = Z f ∇ S · ( iv, ∇ v ) − Z " f ( Ω × r ) · ( iv, ∇ v ) + f i ( Ω × r ) · ∇ ( | v | ) = Z f ( iv, ∇ v ) · ( ∇ S − Ω × r ) = − Z f iv ∗ ∇ v · ( ∇ S − Ω × r )where ( u, v ) := ( uv ∗ + u ∗ v ) /
2. The terms apart from the vortex-free energy in(19) describe the contribution of the vorticity to the energy: The second term G f [ v ] looks formally like a ’weighted’ Ginzburg-Landau (GL) energy functionalwithout magnetic field and accordingly will be called GL-type energy in theforthcoming and R f [ v ] is the rotation energy.Using the splitting (19), vortices of u (if present) are vortices of v and theyare described via the functionals G f [ v ] − R f [ v ]. We have the following main result:
Main result:
Let u ε be a minimizer of (3) and f ε a minimizer of (15) for8 in (16) and S in (17) and under the normalization constraints. Let C and δ be positive constants independent of ε with 0 < δ ≪ o (1) denotea quantity which goes to zero as ε → n ≥ n = C [ | ln ε | + ( n −
1) ln | ln ε | ] =: Ω + C ( n −
1) ln | ln ε | (20)with C := s + 2 sµ /s λ , we have the following results:i) If Ω ≤ Ω − C δ ln | ln ε | and ε sufficiently small, then u ε has no vortices in D \ ∂ D and the Gross-Pitaevskii energy is E GP [ u ε ] = E GP [ f ε e iS ] + C. (21)ii) If Ω n + C δ ln | ln ε | ≤ Ω ≤ Ω n +1 − C δ ln | ln ε | for n ≥
1, then, for ε sufficiently small, u ε has n vortices with winding number one located in r , ..., r n ∈ D \ ∂ D , r i = ( x i , y i ) , i = 1 , .., n . Setting ˜ r i = (˜ x i , ˜ y i ) with ˜ x i = x i √ Ω , ˜ y i = y i λ √ Ω, the configuration (˜ r , ..., ˜ r n ) minimizes the function w ( a , .., a n ) = − πµ X i = j ln[( X i − X j ) + λ − ( Y i − Y j ) ] + πµ λ n X i =1 ( X i + Y i ) −− π ln Ω4Ω s/ n X i =1 ( X i + Y i ) s/ with a i = ( X i , Y i ) , i = 1 , .., n and the Gross-Pitaevskii energy is E GP [ u ε ] = E GP [ f ε e iS ] + π µn | ln ε | − s (1 + λ )( s + 2) µ /s Ω ! ++ π µn ( n −
1) ln Ω + w (˜ r , .., ˜ r n ) + C + o (1) . (22)The proof is split into several estimates which are shown in the following sec-tions. There, positive constants are denoted by C (sometimes carrying primes)and they may change from line to line. In this section, we show the leading asymptotics for the GP energy and density.We will see, in particular, that it is not affected by vortices whose influencecan only be seen in the next lower order. The leading term in the energy comesfrom the TF contribution in (10) which is ∼ /ε whereas vortices contributeto the order Ω ∼ | ln ε | (see also Section 4 ). However, the determination of the9recise expressions in (22) requires a more detailed analysis which is carriedout in
Sections 5-8 .For the following estimates, we introduce the function b ( r ) := 12 ( µ − V ( r )) (23)whose positive part is the TF density, i.e. [ b ( r )] + := ρ TF . Estimate 1:
For Ω( ε ) satisfying C V | ln ε | ≤ Ω( ε ) < C | ln ε | where C V de-pends on the parameters of the external potential V , C > C V , and for ε sufficiently small, E GP [ u ε ] = E TF [ ρ TF ] + C | ln ε | (24)and Z R (cid:16) | u ε | − ρ TF (cid:17) = o (1) . Proof:
This can be shown similar as Prop. 2.3 in [9]. The lower bound can be triviallyobtained by neglecting the first positive term in (5) E GP [ u ε ] ≥ E TF [ ρ TF ] − C Ω( ε ) . The upper bound can be obtained by using E GP [ u ε ] ≤ E GP [ u ε ] | Ω=0 and q ρ TF as a trial function, E GP [ u ε ] | Ω=0 ≤ E TF [ ρ TF ] + C | ln ε | . Concerning the density asymptotics, we estimate the following: using the neg-ativity of b ( r ) outside the TF domain, we have Z R ( | u ε | − ρ TF ) ≤ Z R (cid:16) | u ε | − b ( r ) | u ε | + ( ρ TF ) (cid:17) . On the other hand, we deduce4 ε E TF [ | u ε | ] = Z R h | u ε | + | u ε | V i = Z R | u ε | − Z R b ( r ) | u ε | + µ, that is Z R ( | u ε | − ρ TF ) ≤ ε E TF [ | u ε | ] + Z R ( ρ TF ) − µ = 4 ε ( E TF [ | u ε | ] − E TF [ ρ TF ]) ≤ Cε | ln ε | and the last inequality follows from (24). Thus, the GP density approachesthe TF density for ε → Estimate 1 . (A similar result is also truefor higher angular velocities as is shown in [9]). In the same way, we see thefollowing result which is used to show
Estimate 3 below Z D (cid:16) | u ε | − ρ TF (cid:17) + Z R \D | u ε | = 4 ε E TF [ | u ε | ] + Z D ( ρ TF ) − µ ε ( E TF [ | u ε | ] − E TF [ ρ TF ]) ≤ Cε | ln ε | , so that Z R \D | u ε | ≤ Cε | ln ε | . (25)Now we return to the non-rotating ground state described by (8). We havethe following point-wise estimate for f ε within the TF domain: Estimate 2:
Let f ε be a minimizer of (8) under the normalization constraint.It is the unique positive solution of∆ f = f ε ( V + 2 f − ε ν GP ) in R (26)with the chemical potential ν GP in (9). If ε is sufficiently small, then | q ρ TF ( r ) − f ε ( r ) | ≤ Cε / q ρ TF ( r ) (27)for r ∈ D in := { r ∈ R : V ( r ) ≤ µ − ε / } . That is, we may replace thevortex-free density f ε by the Thomas-Fermi density ρ TF within a region al-most as large as the Thomas-Fermi domain making only an error of order o (1). Proof:
As is shown in Ref. [22], there exists a unique minimizer for the functional (8).Since each minimizer fulfills (26) (which is the corresponding Euler-Lagrangeequation) and E GP [ f ε ] = E GP [ | f ε | ] the positivity of the minimizer f ε fol-lows. Now we look at (27). It can be shown similar as in Refs. [2,5] byusing suitable sub- and supersolutions: We consider a disc B δ ( r ) around r ∈ D ′ := { r ∈ R : V ( r ) ≤ µ − t, t > } with radius δ < t and con-struct a subsolution w ( r ) = √ ρ tanh q with r ∈ B δ ( r ), q := δ −| r − r | δε and ρ := min B δ ( r ρ TF . Using w ( r ), we see that ∆ w ≥ w ε (cid:16) V + 2 w − ε ν GP (cid:17) is fulfilled since 4 ε ν GP > µ + 2[ ρ tanh q − ρ TF ] for ε sufficiently small. On ∂B δ ( r ) there is | r − r | = δ and w | ∂B δ = 0 < f ε . So w is a subsolution for(26) in B δ ( r ) and √ ρ − f ε ( r ) ≤ √ ρ − w ( r ) = 2 √ ρe − δ/ε e − δ/ε ≤ √ ρe − δ/ε ≤ q ρ TF e − δ/ε . (28)Since ρ TF is smooth in D \ ∂ D , we can approximate ρ TF ( r ) by ρ , making asmall error of the order o (1). So, q ρ TF ( r ) − f ε ( r ) q ρ TF ( r ) ≤ q ρ TF ( r ) − √ ρ q ρ TF ( r ) + 2 √ ρ q ρ TF ( r ) e − δ/ε ≤ C δ √ t + e − δ/ε ! . must be chosen such that e − δ/ε is exponentially small. We choose δ = ε / and t = ε / as in [2]. Likewise we construct a supersolution p ( r ) = √ m coth[arcoth( q M / m) + √ mq] with r ∈ B δ ( r ), m := max B δ ( r ρ TF , M =max D ρ TF and q as above. Using p ( r ), we see that ∆ p ≤ p ε (cid:16) V + 2 p − ε ν GP (cid:17) is fulfilled since 4 ε ν GP ≤ µ + 2( m coth (arcoth( q M / m) + √ mq) − ρ TF ) for ε sufficiently small. On ∂B δ ( r ) there is | r − r | = δ and p | ∂B δ = √ M ≥ f ε . So p ( r ) is a supersolution for (26) in B δ ( r ). Proceeding as before, we get f ε ( r ) − q ρ TF ( r ) q ρ TF ( r ) ≤ p ( r ) − q ρ TF ( r ) q ρ TF ( r ) ≤ C δ √ t + e − δ/ε ! . Choosing δ and t appropriately again, we get (27) for any r ∈ D in .It is intuitively clear that only the vortex-free density f ε and not the ’full’GP density | u ε | can satisfy a pointwise estimate as above: u ε may have vor-tices whereas ρ TF carries no vorticity at all. However, what can be shown isthe fact that, in the TF regime, where ε → | u ε | is exponentially smalloutside the TF domain (see also Refs. [18] and [9]): Estimate 3:
For r ∈ Θ ε := { r ∈ R : V ( r ) > µ + ε / } and ε sufficientlysmall, there is | u ε ( r ) | ≤ Cε / | ln ε | / exp b ( r ) Cε / ! , where b ( r ) is defined in (23). Proof:
By using (6) we have −
12 ∆ | u ε | = −|∇ u ε | − V ε | u ε | − | u ε | ε +2 µ GP | u ε | − i ( u ε ( Ω × r ) ·∇ u ∗ ε + u ∗ ε ( Ω × r ) ·∇ u ε ) . The estimate 2Ω( ε ) | iu ∗ ε ∇ u ε × r | ≤ |∇ u ε | + Ω( ε ) | r | | u ε | leads to −
12 ∆ | u ε | ≤ (cid:20) Ω( ε ) | r | ε − V − | u ε | + 2 µ GP ε (cid:21) | u ε | ε in R . From (7) and
Estimate 1 follows ε µ GP = ε E GP [ u ε ] + 14 Z R | u ε | ≤ ε E TF [ ρ TF ] + o (1) + 14 Z R | u ε | = ε µ TF + o (1) + 14 Z R ( | u ε | − ( ρ TF ) ) ≤ ε µ TF + o (1) . −
12 ∆ | u ε | ≤ (cid:20) Ω( ε ) | r | ε − V ε µ TF + o (1) − | u ε | (cid:21) | u ε | ε ≤ C b ( r ) ε | u ε | < ε where b ( r ) < − ε / . That is | u ε | fulfills − ε ∆ | u ε | − C ′ b ( r ) | u ε | ≤ ε . (29)So | u ε | is subharmonic in Θ ε for ε sufficiently small. That means, there is forall r = | r | with B ̺ ( r ) ⊂ Θ ε that | u ε ( r ) | ≤ π̺ Z B ̺ ( r ) | u ε | ≤ √ π̺ (cid:18)Z r ∈ Θ ε | u ε | (cid:19) / ≤ C̺ ε / | ln ε | / using (25). If we now take r ∈ Σ ε := { r ∈ R : V ( r ) ≥ µ + ε / } and choose ̺ = ε / we get | u ε ( r ) | ≤ Cε / | ln ε | / so that | u ε ( r ) | → ε for ε →
0. Moreover, from (29) it follows that | u ε | is a subsolution of − ∆ w + C ′′ ε − / w = 0 in Σ ε w = Cε / | ln ε | / on ∂ Σ ε . (30)On the other hand, one can verify that˜ u = Cε / | ln ε | / exp b ( r ) Cε / ! is a supersolution of (30). Therefore 0 ≤ | u ε ( r ) | ≤ ˜ u for r ∈ Θ ε .So, since | u ε | is exponentially small in ε outside of the TF domain D , theabove energy splitting (19) can be put now into the form E GP [ u ε ] = E GP [ f ε e iS ]+ Z D " f ε |∇ v ε | + f ε ε (1 − | v ε | ) − Z D if ε v ∗ ε ∇ v ε · ( ∇ S − Ω × r )+ o (1)(31)=: E GP [ f ε e iS ] + G f [ v ε ] − R f [ v ε ] + o (1)where v ε = u ε /f ε e iS and o (1) → ε →
0. Thus in the following, it sufficesto restrict our considerations to the Thomas-Fermi domain D .13 Vorticity in the Thomas-Fermi regime
We know from experiments that angular momentum is quantized in the formof vortices when the gas is subjected to an external rotation. Hence we mayapproximate the vorticity field by N v isolated point vortices. However, it is adifficult task in general to prove the validity of this approximation rigorouslyfrom more basic properties. It has been shown in the work of [19] for BECs inharmonic anisotropic traps that the vorticity is indeed concentrated in a finite(independent of ε ) number of vortex cores if one assumes that the angularvelocity is bounded by Ω ≤ C | ln ε | asymptotically. This has been achieved byusing a number of technical vortex core constructions. We will not generalizethese methods to the more general traps considered here but instead we liketo argue by physical reasoning how the number of vortices scales with Ω( ε ).In experiments Ω and ε are independent parameters. Usually, the interactionbetween the particles is tuned and afterwards Ω is increased (independentlyof ε ) beyond the critical value. So in principle one could study the wholeparameter domain spanned by Ω and (here) positive ε . However, we restrictedto the TF regime where the scaled angular velocity Ω = Ω( ε ) depends on ε insuch a way that for ε →
0, Ω → ∞ and hence we cover only a fraction of thepossible parameter domain. Which dependence of Ω( ε ) may occur in the TFregime ? There are essentially three regimes in Ω for non-harmonic traps whereinteresting effects appear (see [9]), namely Ω ∼ | ln ε | , Ω ∼ /ε, Ω ≫ /ε (thefirst regime also applies to harmonic traps). One may ask for a connectionbetween different vortex core sizes, the magnitude of Ω( ε ) and the kind ofdefects appearing in the condensate. For Ω ∼ | ln ε | , one may deduce similarestimates for the vortex energy using core sizes of the order σ = κε or σ = ε α with constants κ, α > ∼ /ε there appears a ’hole’ around theorigin and the core size of the vortices itself is of the order √ ε . For even largervelocities Ω ≫ /ε , the condensate is expelled to a small layer at the boundaryand there remains a ’giant vortex’ state filling out almost all of the condensate. One may argue that vortices with larger core radii, say e.g. σ ∼ / | ln ε | couldin principle exist at lower angular velocities of the order Ω ∼ ln | ln ε | . However,the characteristic length, where perturbations of the condensate wave function aresmoothed out, is given by the healing length ξ or ε respectively. Hence we expectthe cores to be of the order ε and larger cores are not stable in the setting describedhere. Moreover, for angular velocities of the order | ln ε | we may also not expectthe appearance of pathological cases like non-isolated vortices forming dense 1-dimensional structures because they would have a much higher energy than wouldbe favourable at this order of Ω. It seems that the underlying equations are tooregular to support such kinds of defects even at much higher angular velocities.
14n the regime of large vorticity one may also consider a kind of correspondenceprinciple for a large number of vortices which is argued by Feynman [13] inthe context of rotating superfluid He: a dense lattice of uniform distributedvortices should ’mimic’ solid-body rotation on average, although the flow isstrictly irrotational away from the vortex cores. The circulation around aclosed contour C which encloses a large number of vortices N v is Γ = H C ∇ S u · τ = 2 πdN v for vortices with winding number d . On the other hand, if thevortex lattice mimics solid-body rotation there is Γ = 2Ω A where A is thearea enclosed by the contour C . In this approximation, the vortex densityper area is n v = N v /A = Ω / ( πd ) and the area per vortex is 1 /n v = πd/ Ωand so decreases with increasing Ω. The crucial ingredient in this argumentis the assumption of a uniform distribution of vortices. But this is justifiedonly if the number of vortices is very large, i.e. if Ω is very large which meansfor the TF regime that, indeed, N v , Γ, and Ω have to increase as ε → A is finitely large, i.e. bounded from belowby a positive constant (independent of ε ). Actually, A is the whole condensatedomain and the contour C is the boundary of that domain. So, from Γ = 2Ω A we see that Γ ≃ Ω, i.e. the circulation is of the same order than the angularvelocity if the vortices are distributed uniformly.On the other hand, considering the case that N v and Γ respectively can bebounded from above by a finite constant (independent of ε ), then vortices cannot be distributed uniformly but instead they form a polygonal lattice (see e.g.the pictures in [28]). So the above argument concerning solid-body rotationgives only an upper bound for Γ. However, this bound is still quite goodin experimental realizations as is demonstrated in [10]. In order to estimateroughly the order of magnitude of Ω for a finite number of vortices to appear,one may calculate the GP energy of a single vortex. This has been done mostoften in the approximation of a homogeneous system (see e.g. Refs. [30,31])or for a condensate in harmonic traps (see e.g. Refs. [26,4]). In any case,the leading contribution comes from the angular kinetic energy. This can bealready seen heuristically by considering a vortex of circulation Γ = 2 πd andcore radius σ ∼ ε which is located at the origin of a flat trap with radius R . Writing the vortex in the form v ( r, θ ) = ̺ ( r ) e iθd with ̺ ( r ) ∼ r d if 0 ≤ r ≤ εR and ̺ ( r ) ∼ R − if εR ≤ r ≤ R , the kinetic energy is then R |∇ v | ∼ R − ( d | ln ε | + C ), whereas the rotation term gives − iv ∗ Ω ( ε ) · ( ∇ v × r ) = − d Ω( ε ). Thus one may expect a vortex of winding number d to appear when Z ( |∇ v | − iv ∗ Ω ( ε ) · ( ∇ v × r )) ∼ d R | ln ε | − d Ω( ε ) < ε ) > C dR | ln ε | and the constant C is fixed by the external potential accordingly. Hence onevortex or a finite number of them are favourable to exist if the angular velocityis of the order Ω( ε ) ≃ C | ln ε | . Furthermore, from Γ < A ≤ C we get For | ln ε | ≪ Ω( ε ) ≪ /ε the number of vortices is no longer bounded as ε → ≤ C/ Ω, that is the vortices are enclosed within a disc centered at the originhaving a radius of the order r v ≤ C √ Ω ≃ C ′ q | ln ε | . (32)So with regard to the above discussion, we model the vorticity in terms of afinite number of vortices within the TF domain denoting their positions as r i = ( x i , y i ) ∈ D \ ∂ D , i = 1 , .., n , n ∈ N . The vortex cores are modelledas non-overlapping discs B i = B ( r i , σ ) with core radius σ ∼ ε , all containedwithin D : ¯ B ( r i , σ ) ⊂ D for all i, ¯ B ( r i , σ ) ∩ ¯ B ( r j , σ ) = ∅ for all i = j (33)assuming that | r i − r j | > σ . Otherwise, their energy would surpass the order of | ln ε | and would hence not be favourable for the angular velocities consideredhere. In a vortex point r i , the condensate wave function vanishes | u | ( r i ) = | v | ( r i ) = 0 ∀ i and the circulation condition can be written as Z ∂B i ∇ S u · τ = Z ∂B i ∇ S v · τ = Z ∂B i ∂S v ∂τ = 2 πd i ∀ i, (34)where τ is the unit tangent vector to B i and d i is the degree of the vortex in r i . The domain outside the cores is denoted as˜ D := D \ [ i B i . In that region, there holds | u | → f , i.e. | v | → ≤ | v | ≤ − o (1) in B ( r i , σ ) (35)and | v | = 1 − o (1) in ˜ D (36)where o (1) goes to zero for σ → ε → o (1) depends on the steepness of the radial falloff of the vortex core profile.For the core radii we are going to use, namely σ = ε α , α >
0, the error due tothe core profile is negligible within the orders considered. G f [ v ]In this section, we consider the functional G f [ v ε ] = Z D " f ε |∇ v ε | + f ε ε (1 − | v ε | ) (37)16hich is part of the energy splitting (31). Because of (27), we can replace f ε by q ρ TF in (37) and the error is of the order o (1). Using the polar decomposition v ε = | v ε | e iS vε , (37) is equivalent to G f [ v ε ] = Z D " ρ TF ∇| v ε | ) + | v ε | ( ∇ S v ε ) ] + ( ρ TF ) ε (1 − | v ε | ) − o (1) . (38)Minimizing this at fixed ρ TF results in − ∆ | v ε | + | v ε | ( ∇ S v ε ) − ( ρ TF ) ε | v ε | (1 − | v ε | ) = 0and ∇ · h | v ε | ∇ S v ε i = 0 . (39)We have the following estimate: Estimate 4:
Let f ε be a minimizer of (15), u ε a minimizer of (3) and v ε = u ε /f ε e iS . Let σ = Cε α with constants C, α > v ε satisfy (33) - (36)in presence of vortices in r i having winding numbers d i , i = 1 , .., n . Then, for ε sufficiently small and Ω ≤ C | ln ε | asymptotically, the GL-type energy canbe bounded from below by G f [ v ε ] ≥ π | ln σ | n X i =1 d i ρ TF ( r i ) + π ln σε n X i =1 | d i | ρ TF ( r i ) −− π X i = j d i d j ln | r i − r j | ρ TF ( r i ) + o (1) . (40) Proof:
First we are going to estimate G f [ v ε ] in the vortex-free domain where | v ε | =1 − o (1). Then (38) reduces to G f [ v ε ] | ˜ D = 12 Z ˜ D ρ TF ( ∇ S v ε ) = 12 Z ˜ D ρ TF V (41)up to an error of order o (1) and we use the relation ∇ S v ε = V where V is the(linear) superfluid velocity of the condensate (see (39)). Indeed, it is shown inRef. [23] that BECs are 100 % superfluid in their ground state. Minimizingthe functional with respect to V gives ρ TF ∇ · V + ∇ ρ TF · V = 0 in ˜ D . (42)Because of the circulation condition (34)2 πd i = Z ∂B i ∇ S v ε · τ = Z ∂B i V · τ = Z B i ∇ × V · d o , there is ∇ × V ( r ) = 2 πd i δ ( r − r i ) in B i (43)17enoting the oriented surface element as d o . Since the cores B i are small andthe TF density ρ TF is smooth in D \ ∂ D , ρ TF is nearly constant within them,and we may approximate ∇ · V = 0 in B i . (44)This allows to define a stream function ψ ε , which is the dual to the phase S v ε ,so ∇ S v ε = ∇ × ψ ε . Using this form for V together with (43) and (44), thestream function becomes ψ ε ( r ) = − d i ln | r − r i | in B i . (45)Now we calculate the integral12 Z ˜ D ρ TF V = 12 Z ˜ D ρ TF ( V × ∇ ψ ε ) · d o = 12 Z ˜ D ∇ × [ ρ TF ψ ε V ] · d o −− Z ˜ D ψ ε ∇ ρ TF · V − Z ˜ D ψ ε ρ TF ( ∇ · V )= 12 Z ∂ ˜ D ρ TF ψ ε V · τ − Z ˜ D ψ ε [ ∇ ρ TF · V + ρ TF ( ∇ · V )]= 12 Z ∂ D ρ TF ψ ε V · τ + 12 n X i =1 Z ∂B i ρ TF ψ ε V · τ = 12 n X i =1 ρ TF ( r i ) Z ∂B i ψ ε V · τ + o (1) , (46)where we used Stokes theorem, (42) and ρ TF = 0 on ∂ D . Using again the factthat the cores B i are small and ρ TF is smooth in D \ ∂ D , we replace it by itsvalue in the core center ρ TF ( r i ) and the error is of the order o (1).If we insert ψ ε from (45) we get12 X i ρ TF ( r i ) Z ∂B i ψ ε V · τ = − X i ρ TF ( r i ) d i Z ∂B i ln | r − r i | V · τ −− X i ρ TF ( r i ) X j = i d j Z ∂B i ln | r − r j | V · τ + o (1) . (47)The first term on the r.h.s. describes the ’diagonal part’ ( i = j ) and the secondone the ’non-diagonal part’ ( i = j ). Using ln | r − r i | = ln σ and | r − r j | = | r i − r j | + o (1) for r ∈ ∂B i , there simply remains in each case the circulationcondition which gives 2 πd i . So we have12 X i ρ TF ( r i ) Z ∂B i ψ ε V · τ = − π ln σ X i d i ρ TF ( r i ) − π X i = j d i d j ln | r i − r j | ρ TF ( r i ) . (48)Since σ ≪ − ln σ can be replaced by | ln σ | and we recover the first andthird term in (40). In this result, we see the familiar energy dependence onthe winding number squared d and the logarithmic divergence due to thevortex cores | ln σ | (see also the approach in Ref. [26] for the harmonic trap18ase). If there are more vortices present than one, their interaction energy ismodelled by W ( r , .., r n ) = − π X i = j d i d j ln | r i − r j | ρ TF ( r i ) (49)in (48). It has the form of a Coulombian interaction in 2-dimensional systemswhere vortices with the same sign of the winding number repel each other andvortices with opposite sign attract each other. In Ref. [7], the analogue to thisfunction is called renormalized energy because it remains after the core energy,which is the leading order, is separated: W < | ln σ | as long as | r i − r j | > σ .We also see that W is bounded from below by a constant if all winding num-bers have the same sign.It remains to find a lower bound of G f [ v ε ] in the vortex cores B i : G f [ v ε ] | S i B i = X i Z B i " ρ TF |∇ v ε | + ( ρ TF ) ε (1 − | v ε | ) − o (1) ≥ X i Z B i ρ TF |∇ v ε | − o (1) = X i ρ TF ( r i )2 Z B i |∇ v ε | − o (1) (50)where we used again the fact that the TF density ρ TF varies only of the order o (1) in the small discs B i . The integral over B i can be estimated as follows: Z B i |∇ v ε | ≥ Z B i \ B ε |∇ v ε | ≥ Z σε Z π r (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∂v∂φ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) rdrdφ with polar coordinates ( r, φ ) on the annulus, B ε a disc with radius ε centeredat r i , and we use the polar decomposition v ( r, φ ) = | v | ( r ) e id i φ for a vortexwith winding number d i in the disc B i with radius σ . Using R π | ∂v∂φ | ≥ π | d i | and Cauchy-Schwartz inequality, we get Z σε Z π r (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∂v∂φ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) drdφ ≥ π | d i | ln σε . (51)Combining (47), (48), (50) and (51), we complete the proof of (40). R f [ v ]The estimate for the rotation term R f [ v ε ] = Z D if v ∗ ε ∇ v ε · ( ∇ S − Ω × r ) (52)19n (31) proceeds similar as in [35]. Because of (27), we replace f ε by q ρ TF in(52) and the error is of the order o (1). Then we get from (14) ∇ · [ ρ TF ( ∇ S − Ω × r )] = 0 in ˜ D (53)from where we see that there is a real function χ ( x, y ) satisfying ρ TF ( ∇ S − Ω × r ) = Ω ∇ ⊥ χ (54)with ∇ ⊥ χ = ( − ∂ y χ, ∂ x χ ), ∂ x = ∂∂x , ect. We impose the accompanying bound-ary condition χ = 0 on ∂ D . To determine the auxiliary function χ we use (54)and rewrite it as ( ∇ S − Ω × r ) ⊥ = Ω ρ TF ∇ χ where r ⊥ = ( − y, x ) if r = ( x, y ). Applying the operator ∇ , we get ∂ x ( ∂ y S − Ω x ) + ∂ y ( − ∂ x S − Ω y ) = Ω ∇ · ∇ χρ TF ! or ∇ · ∇ χρ TF ! = − . (55)We have the following estimate for the rotation term: Estimate 5:
Let f ε be a minimizer of (15), u ε a minimizer of (3), v ε = u ε /f ε e iS and χ the solution of (55). Let σ = Cε α with constants C, α > v ε satisfy (33) - (36) in presence of vortices in r i having winding numbers d i , i = 1 , .., n . Then for ε sufficiently small and Ω ≤ C | ln ε | asymptotically, therotation energy is R f [ v ε ] = 2 π Ω n X i =1 d i χ ( r i ) + o (1) . (56) Proof:
We can see that the contribution in the vortex cores becomes small: (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) n X i =1 Z B i iρ TF v ∗ ε ∇ v ε · ( ∇ S − Ω × r ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ n X i =1 Z B i | ρ TF || v ε ||∇ v ε ||∇ S − Ω × r |≤ nµ Ω (cid:18)Z B i | v ε | (cid:19) / (cid:18)Z B i |∇ v ε | (cid:19) / ≤ Cσ | ln ε | / ≤ o (1)where we used that | R B i |∇ v ε | | ≤ C | ln ε | which will be shown in (59).For estimating the rotation term outside of the vortex discs, we use again | v ε | = 1 − o (1) to get 20 ˜ D iρ TF v ∗ ε ∇ v ε · ( ∇ S − Ω × r ) = Ω Z ˜ D iv ∗ ε ∇ v ε · ∇ ⊥ χ = − Ω Z ˜ D ∇ S v · ∇ ⊥ χ = − Ω Z ˜ D ∇ ⊥ · ( χ ∇ S v ) + Ω Z ˜ D χ ∇ ⊥ · ∇ S v = − Ω Z ∂ ˜ D χ ∇ S v · τ = − Ω Z ∂ D χ ∂S v ∂τ + Ω n X i =1 Z ∂B i χ ∂S v ∂τ = Ω n X i =1 χ ( r i ) Z ∂B i ∂S v ∂τ + o (1) = 2 π Ω n X i =1 d i χ ( r i ) + o (1) (57)where we used (54), ∇ ⊥ · ∇ S v = 0, χ = 0 on ∂ D and (34). Furthermore, sincethe cores B i are small and χ is smooth in D \ ∂ D , we replace it by its valuein the core center χ ( r i ) and the error is of the order o (1). We thus arrive at (56).The expressions (40) and (56) for the vortex contributions suggest that thelowest vortex energy is attained for vortices with winding number d i = 1 forall i , which will be explicitly shown in Section 8.2 . In estimating (40), we used ∇ × V = 2 πd i δ ( r − r i ) in the vortex core B i , which is valid for any vortexcore radius σ . On the other hand, due to the characteristic scale ε the core isnot much larger than σ ∼ ε α , α >
0. Then, the gradient term of G f [ v ] outside the core dominates over the contribution of the core itself. Concerning theinteraction energy, one could first of all ask which terms of G f [ v ] − R f [ v ] inthe splitting of the functional E GP [ u ] in (19) will contribute to the interactionbetween vortices. We have just seen that the rotation energy outside of vortexcores is of the ’diagonal’ form given in (57), whereas it is of order o (1) inthe cores. So the interaction must be modelled by the GL-type energy G f [ v ].In addition, the interaction is only relevant in the domain outside the cores.There we have | v | ≃ G f [ v ] plays the signif-icant role. In particular, the form of the core and interaction energy in (40)was deduced by minimizing (41) with respect to ∇ S v ε = V . The core energydominates as long as σ ≪ | r i − r j | for all i = j , i.e. as long as the vortex coresize is much smaller than the distance between vortices. This is vastly fulfilledin the regime ε → ≃ C | ln ε | since then there is σ ≃ Cε α , α > | r i − r j | ≥ C/ q | ln ε | which will be shown in Section 8.3 . Remark:
Our analysis may be compared with the works of [4] and [18]. Westart from the original energy functional (1) and rescale it to arrive at (2) and(3) respectively. In [18,19], a functional is used, motivated in [4] and justifiedby the normalization condition, which already in the beginning looks like aGL-type functional. The so encountered additional term is ’thrown away’ be-cause it does not contribute to the vortex energy. But one has to keep in mindthat the true leading order is 1 /ε which can be explicitly seen in (3). In theestimate of the lower bound of (37), we consider equations (41) and (42). Thebehaviour of V in the discs is derived from the quantization condition (34)and is given in (45) in terms of the stream function. These ingredients are used21n the estimate of (46). Instead, methods of Ref. [7] are adapted in [19] to thefunctional studied by considering a ’linear problem’ as in [7]. By introducing asuitable function, it gives eventually the interaction energy between vortices.Concerning the forthcoming estimates, we will benefit from inequality (62).The corresponding equality for s = 2 is used in [4,18]. In [18,19], it is assumedthat Ω ≤ C | ln ε | asymptotically. Otherwise, there is no a priori assumptionon the fine structure of vorticity. We have argued in Section 4 that the as-sumption of a non-zero circulation which is bounded from above by a naturalnumber independent of ε leads to an angular velocity of the order Ω ∼ | ln ε | .So, our assumptions on the vortex fine structure, concerning number and sizeof vortex cores, are actually compatible with this order of Ω. G f [ v ] − R f [ v ]The energy without vortex is always larger or equal to E GP [ u ε ] and it can beused as a trial function for the whole energy (see also the proof of Estimate1 ): E GP [ u ε ] = E GP [ f ε e iS ] + G f [ v ε ] − R f [ v ε ] + o (1) ≤ E GP [ f ε e iS ] + C + o (1)so G f [ v ε ] − R f [ v ε ] ≤ C + o (1). A more precise upper bound is obtained asfollows: Estimate 6:
Let f ε be a minimizer of (15), u ε a minimizer of (3) and v ε = u ε /f ε e iS . Then, for ε sufficiently small and Ω ≤ C | ln ε | asymptotically, thevortex energy G f [ v ε ] − R f [ v ε ] can be bounded from above by G f [ v ε ] − R f [ v ε ] ≤ π | ln ε | k X i =1 d i ρ TF ( r i ) − π Ω k X i =1 d i χ ( r i ) + W ( r , .., r k ) + C + o (1)(58)where d i ≥ i , W ( r , .., r k ) from (49) and i, j = 1 , .., k . Proof:
We fix k ≥ , k ∈ N vortex positions r , ..., r k in D , each is center of a discwith fixed radius R > D and do not overlap i.e. ¯ B ( r i , R ) ⊂ D and ¯ B ( r i , R ) ∩ ¯ B ( r j , R ) = ∅ for all i = j . We use the trial function ˆ v = | ˆ v | e i ˆ S v with | ˆ v | = 1 in ˜ D andˆ v ( r ) = r − r i | r − r i | for | r − r i | ≥ ε r − r i ε otherwise22n the discs B i ( r i , R ). Since |∇ ˆ v | = ( ∇ ˆ S v ) in ˜ D and the phase is smooth andbounded outside the discs with finite size, there is G f [ˆ v ] | ˜ D = Z ˜ D ρ TF ∇ ˆ S v ) ≤ C. However, in the discs there is G f [ˆ v ] | S i B i = k X i =1 ρ TF ( r i )2 Z B i |∇ ˆ v | + k X i =1 ( ρ TF ) ( r i )4 ε Z B i (1 − | ˆ v | ) + o (1)where Z B i |∇ ˆ v | − π = Z B i \ B ε |∇ ˆ v | = Z B i \ B ε | r − r i | = Z π Z Rε rdrdφr − rε cos φ + ε = π ln( r − ε ) | Rε ≤ π ln r | Rε = 2 π | ln ε | + 2 π ln R, (59)and the other contributions are at most of the order of a constant. With theabove trial function, the rotation energy in ˜ D is the same as in (56) apart fromthe sum running from i = 1 to k . The contribution inside the vortex discs issimply (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) k X i =1 Z B ( r i ,R ) iρ TF ˆ v ∗ ∇ ˆ v · ( ∇ S − Ω × r ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ k X i =1 Z B ( r i ,R ) | ρ TF || ˆ v ||∇ ˆ v ||∇ S − Ω × r |≤ kµ Ω (cid:18)Z B i | ˆ v | (cid:19) / (cid:18)Z B i |∇ ˆ v | (cid:19) / ≤ kµ Ω( π ( R − ε ) + 3 πε / / (2 π | ln ε | + 2 π ln R + 4 π ) / ≤ Cε | ln ε | / . Then, to recover (58) we finally use the fact that W ( r , .., r k ) can be boundedfrom below by a constant. In the above estimates (24), (40), (56) and (58), the external trap potential V enters via the Thomas-Fermi density ρ TF which was not specified until now.These estimates are valid as long as V satisfies (asymptotical) homogeneity(see also the remark at the end of this section). As an application, we considernow the potential in (16). The associated TF density is ρ TF ( x, y ) = 12 (cid:16) µ − ( x + λ y ) s/ (cid:17) . (60)23rom (12) we have µ = ( s +2 s λπ ) s/ ( s +2) . For s → ∞ , µ → λ/π , hence µ isalways smaller than one. The auxiliary function χ is determined from (55) to χ ( x, y ) = 11 + λ " s + 2 ( x + λ y ) ( s +2) / − µ x + λ y ) + s s + 2) µ ( s +2) /s . (61)It can be estimated from above in terms of the TF density ρ TF by χ ( x, y ) ≤
11 + λ s /s s + 2 ( ρ TF ( x, y )) (2+ s ) /s , (62)where strict equality only holds for the harmonic trap s = 2 ! This upper boundwill be very useful in Section 8.2 where the winding numbers of vortices arederived. The phase S can be determined by inserting (61) and (60) into (54)and is already given in (17). The upper bound of the energy using (31) and (58) is E GP [ u ε ] ≤ E GP [ f ε e iS ] + π | ln ε | k X i =1 d i ρ TF ( r i ) − π Ω k X i =1 d i χ ( r i ) + C + o (1) . For a trial function having one vortex with winding number d = 1 at the originwe get E GP [ u ε ] ≤ E GP [ f ε e iS ] + π µ | ln ε | − πsµ ( s +2) /s (1 + λ )( s + 2) Ω + C + o (1) . (63)The energy E GP [ u ε ] will be smaller than the vortex-free energy E GP [ f ε e iS ] ifthe r.h.s. of (63) is smaller or equal to E GP [ f ε e iS ] − o (1). Equivalently, theangular velocity must fulfill Ω ≥ Ω + C + o (1)where Ω = s + 2 sµ /s λ | ln ε | =: C | ln ε | . (64)(Equation (64) has to be multiplied by (16 ε ) / ( s +2) in order to obtain theunscaled angular velocity ˜Ω , see (4)). So for Ω ≥ Ω + C + o (1), minimizersof E GP [ u ] will have vortices, or in other terms: Ω is the leading order in theangular velocity where the solution with one vortex having d = 1 starts to beglobally thermodynamically stable. We may also see the following: Consider( x, y ) ∈ D , let s > δE = π | ln ε | ρ TF ( x, y ) − π Ω χ ( x, y ) + C .Then one can see that at the origin ∇ ( δE )(0 ,
0) = and ∆( δE )(0 ,
0) =24 πµ Ω >
0, i.e. (0 ,
0) is a local minimum for all
Ω. However, (0 ,
0) is a globalminimum for δE (0 , < + C + o (1) < Ω with Ω in (64). However,for the harmonic trap s = 2 one has ∇ ( δE )(0 ,
0) = and ∆( δE )(0 ,
0) = − π | ln ε | (1 + λ ) + 4 πµ Ω. For ∆( δE )(0 , <
0, the origin is a local maximum;for ∆( δE )(0 , > local thermodynamical stability which is Ω > λ µ | ln ε | = Ω / in (64) with s = 2 which was also shown in Ref. [4]. Estimate 7: If σ = ε α , < α < ε sufficiently small and Ω ≤ Ω + CF ( ε )with F ( ε ) of lower order than | ln ε | , then d i = 1 for all i . Proof:
We use σ = ε α with 0 < α < W ≥ C andΩ ≤ C | ln ε | + C F ( ε )with C from (64) and C is another positive constant. This upper bound forΩ is suggested by (64) and F ( ε ) is assumed to be of lower order than | ln ε | .In the next section, we will see that F ( ε ) = ln | ln ε | . So π | ln σ | n X i =1 d i ρ TF ( r i ) + π ln σε n X i =1 d i ρ TF ( r i ) − π Ω n X i =1 d i χ ( r i ) ≤ C or π | ln σ | X i ( d i − d i ) ρ TF ( r i ) + π | ln ε | X i d i ρ TF ( r i ) ≤ C + 2 π Ω X i d i χ ( r i ) ≤ C + 2 π Ω1 + λ ss + 2 µ /s X i d i ρ TF ( r i ) ≤ C + ( C | ln ε | + C F ( ε )) 2 π λ ss + 2 µ /s X i d i ρ TF ( r i ) ≤ C + π | ln ε | X i d i ρ TF ( r i ) + C ′ F ( ε ) X i d i ρ TF ( r i ) . (65)Here we used (62) and ( ρ TF ) ( s +2) /s ≤ (cid:16) µ (cid:17) /s ρ TF . We also note that the lastinequality in (65) is valid only for C from (64).So we see that (65) reduces to n X i =1 ( d i − d i ) ρ TF ( r i ) ≤ o (1)25or ε sufficiently small. Therefore, if the vortices are not located at the bound-ary of the Thomas-Fermi domain where ρ TF vanishes, there must be d i = 1for all i for sufficiently small ε . n vortices Since all vortices have winding number one, we see from the lower and upperbounds of E GP [ u ε ] in (40), (56) and (58) that they coincide in their orders(up to a constant). By applying the transformation ˜ r i = (˜ x i , ˜ y i ) with ˜ x i = x i √ Ω , ˜ y i = y i λ √ Ω, the vortex interaction energy W ( r , .., r n ) in (49) can bedecomposed as follows: W ( r , .., r n ) = − π X i = j ln | r i − r j | ρ TF ( r i )= − π X i = j ln( | x i − x j | + | y i − y j | )( µ − ( x i + λ y i ) s/ )= π µn ( n −
1) ln Ω − π s/ X i (˜ x i + ˜ y i ) s/ − πµ X i = j ln (cid:18) (˜ x i − ˜ x j ) + 1 λ (˜ y i − ˜ y j ) (cid:19) + π s/ X i = j ln (cid:18) (˜ x i − ˜ x j ) + 1 λ (˜ y i − ˜ y j ) (cid:19) (˜ x i + ˜ y i ) s/ = π µn ( n −
1) ln Ω − π µ X i = j ln (cid:18) | ˜ x i − ˜ x j | + 1 λ | ˜ y i − ˜ y j | (cid:19) + o (1) . (66)The first term contains the relevant order ln Ω whereas the remainder is of theorder of a constant. The rotation term (56) becomes − π Ω n X i =1 χ ( r i ) = − π Ω1 + λ X i " s + 2 ( x i + λ y i ) ( s +2) / − µ x i + λ y i ) + s s + 2) µ ( s +2) /s = − πsn (1 + λ )( s + 2) µ ( s +2) /s Ω + π Ω µ λ X i ( x i + λ y i ) − π Ω(1 + λ )( s + 2) X i ( x i + λ y i ) ( s +2) / = − πsn (1 + λ )( s + 2) µ ( s +2) /s Ω + πµ λ X i (˜ x i + ˜ y i ) − π (1 + λ )( s + 2) 1Ω s/ X i (˜ x i + ˜ y i ) ( s +2) / (67)26here in the last step the variable transformation was applied.The remaining part of the lower bound of G f [ v ε ] in (40) is (with d i = 1 ∀ i ) π | ln ε α | X i ρ TF ( r i ) + π ln( ε α ε ) X i ρ TF ( r i ) = π | ln ε | X i ρ TF ( r i )for small ε and using σ = ε α . With u ε and f ε as above, we thus recover (22)for the Gross-Pitaevskii energy in presence of n vortices E GP [ u ε ] = E GP [ f ε e iS ] + π µn | ln ε | − s (1 + λ )( s + 2) µ /s Ω ! ++ π µn ( n −
1) ln Ω + w (˜ r , ..., ˜ r n ) + C + o (1) (68)with w (˜ r , ..., ˜ r n ) = − πµ X i = j ln (cid:18) (˜ x i − ˜ x j ) + 1 λ (˜ y i − ˜ y j ) (cid:19) ++ πµ λ X i (˜ x i + ˜ y i ) − π ln Ω4Ω s/ X i (˜ x i + ˜ y i ) s/ where we put all terms proportional to Ω − m , m > o (1) since Ω ≤ C | ln ε | asymptotically.A necessary condition for the minimizing configuration to have more thanone vortex is min U n E GP [ u ] ≤ min U E GP [ u ] for n ≥ U n is the set of functions with n vortices having winding number oneeach and U is the set of functions with one vortex at the origin with windingnumber one. Using (69), we first want to deduce a rough estimate for thecritical angular velocity Ω n for n vortices to appear. To this aim, we neglectin (68) the term coming from the interaction and we take the energy with allvortices close to the origin, i.e. we approximate ρ TF ( r i ) ≈ µ/ i , whichis a more stringent condition on the l.h.s. of (69). Indeed, from (32) we expectthe vortices to be near to the origin. Hence π | ln ε | µ ( n − π µn ( n −
1) ln Ω n + πsµ ( s +2) /s (1 + λ )( s + 2) Ω ≤ πsnµ ( s +2) /s (1 + λ )( s + 2) Ω n + C, and 1 + λ s + 2 sµ /s | ln ε | n − n + 1 + λ s + 2 sµ /s n −
12 ln Ω n + 1 n Ω ≤ Ω n + C. Using (64) and the fact that Ω ≤ Ω n for n ≥
2, we have the estimateΩ + 1 + λ s + 2 sµ /s n −
12 ln Ω ≤ Ω n + C + C n −
12 ln | ln ε | + C n −
12 ln C − C ≤ Ω n with C and Ω from (64). We thus see that the critical angular velocity hasto be at least of the order C | ln ε | + C ′ ln | ln ε | (we neglect the constant term).This order of magnitude for Ω is assumed in Ref. [19] from the outset beforethe number of vortices is rigorously derived. Using now the ansatzΩ = Ω + C ν ( ε ) ln | ln ε | (70)with ( k −
1) + δ ≤ ν ( ε ) ≤ k − δ for an integer k ≥ < δ ≪ ε , we see the following: Inserting this form of Ω inour energy estimate (68) we get for the upper bound E GP [ u ε ] ≤ E GP [ f ε e iS ] + π µk | ln ε | − πskµ /s (1 + λ )( s + 2) Ω + π µk ( k −
1) ln Ω + C = E GP [ f ε e iS ] − π kµν ( ε ) ln | ln ε | + π µk ( k −
1) ln | ln ε | + π µk ( k −
1) ln C + C or G f [ v ε ] − R f [ v ε ] ≤ − π µkν ( ε ) ln | ln ε | + π µk ( k −
1) ln | ln ε | ++ π µk ( k −
1) ln C + C (71)respectively. Considering the case k = 0 (no vortices), i.e. ν ( ε ) in (70) satisfies − δ ≤ ν ( ε ) ≤ − δ , we see from (70) and (64) that Ω ≤ Ω − C δ ln | ln ε | and G f [ v ε ] − R f [ v ε ] = C , showing (21).Considering k = 1 (one vortex), i.e. δ ≤ ν ( ε ) ≤ − δ and comparing its energywith the lower bound of G f [ v ε ] and R f [ v ε ] we have − π µν ( ε ) ln | ln ε | + C ≥ π | ln ε | n X i =1 ρ TF ( r i ) − π Ω n X i =1 χ ( r i ) + C ≥ π n X i =1 ρ TF ( r i )( − ν ( ε ) ln | ln ε | + C ) ≥ − π µnν ( ε ) ln | ln ε | so 1 − o (1) ≤ n, that is, there is at least one vortex forΩ + C δ ln | ln ε | ≤ Ω ≤ Ω + C (1 − δ ) ln | ln ε | = Ω − C δ ln | ln ε | (using (20)) and ε sufficiently small. Now we compare the lower and upperbound of the energy if k > G f [ v ε ] − R f [ v ε ] ≥ − π µnν ( ε ) ln | ln ε | + π µn ( n −
1) ln | ln ε | ++ π µn ( n −
1) ln C + C. (72)Comparison of (71) and (72) gives − nν ( ε ) + 12 n ( n −
1) + o (1) ≤ − kν ( ε ) + 12 k ( k −
1) + o (1) . Assuming now that n ≤ k −
1, we have ν ( ε )( k − n ) ≤
12 ( k − n )( k + n −
1) + o (1) , so ( k −
1) + δ ≤ ν ( ε ) ≤
12 ( k + n −
1) + o (1) ≤ k − o (1)which is a contradiction for ε sufficiently small, since δ is a fixed constant.On the other hand, assuming n ≥ k + 1 we have ν ( ε )( n − k ) ≥
12 ( n − k )( k + n −
1) + o (1) , so k − δ ≥ ν ( ε ) ≥
12 ( k + n −
1) + o (1) ≥ k + o (1)which is again a contradiction for ε sufficiently small. So we see that there areexactly n ≡ k vortices for ε sufficiently small and from this follows G f [ v ε ] −R f [ v ε ] = − π µnν ( ε ) ln | ln ε | + π µn ( n −
1) ln | ln ε | + π µn ( n −
1) ln C + C and Ω n + C δ ln | ln ε | ≤ Ω ≤ Ω n +1 − C δ ln | ln ε | by using (20). This completes the proof of the main result stated in the endof Section 2 . Special cases:
For the harmonic trap s = 2 there isΩ n = 1 + λ µ [ | ln ε | + ( n −
1) ln | ln ε | ] . From (4), the unscaled angular velocity is then˜Ω n = 2 µ (1 + λ ) ε [ | ln ε | + ( n −
1) ln | ln ε | ]29hich may be compared to the results in Refs. [4] and [8]. For the flat trap s → ∞ , there is Ω n = 1 + λ | ln ε | + ( n −
1) ln | ln ε | ]and the same is true for ˜Ω n since ˜Ω n → Ω n for s → ∞ . (We put in mind thatfor the flat trap, the scaled energy functional converges to the original one,i.e. E GP ′ [ u ′ ] → E GP [ u ], see Section 2 ). By comparing the first critical angularvelocities, we see the following: For the flat trap ˜Ω ∼ | ln ε | , resembling thecorresponding velocity for the rotating bucket and this is no surprise since theflat trap approximates the bucket. On the other hand, there is ˜Ω ∼ ε | ln ε | for the harmonic trap which is much smaller. The ratio of the unscaled firstcritical angular velocities is thus˜Ω ( s → ∞ )˜Ω ( s = 2) = µ ε . The minimization of the energy in (68) with respect to the coordinates ˜ r i =(˜ x i , ˜ y i ) determines the distribution of the vortices in the condensate and there-fore the resulting pattern which appears for a given number n of vortices. Theenergy in (68) is minimal with respect to the coordinates if w (˜ r , ..., ˜ r n ) isminimal. Setting ∇ w = , we obtain πµ X i = j ˜ x i − ˜ x j (˜ x i − ˜ x j ) + λ − (˜ y i − ˜ y j ) + πs ln Ω2Ω s/ X i ˜ x i (˜ x i + ˜ y i ) s/ − = 2 πµ λ X i ˜ x i (73)and πµ λ X i = j ˜ y i − ˜ y j (˜ x i − ˜ x j ) + λ − (˜ y i − ˜ y j ) + πs ln Ω2Ω s/ X i ˜ y i (˜ x i + ˜ y i ) s/ − = 2 πµ λ X i ˜ y i . (74)Multiplying (73) and (74) with ˜ x i and ˜ y i respectively and adding them togethergives X i (˜ x i + ˜ y i ) = 1 + λ n ( n − λ µ s ln ΩΩ s/ X i (˜ x i + ˜ y i ) s/ . (75)On the other hand, multiplying them with ˜ y i and − λ ˜ x i respectively andadding them gives(1 − λ ) X i ˜ x i ˜ y i = 1 + λ µ s ln ΩΩ s/ (1 − λ ) X i ˜ x i ˜ y i (˜ x i + ˜ y i ) s/ − . (76)30he relations (75) and (76) are constraints for the (non-dimensionalized) vor-tex positions. They simplify considerably for the harmonic trap s = 2 (andonly for this trap!). They were already deduced in Ref. [4] and we only statethem for completeness: P i ˜ x i = P i ˜ y i = 0 and X i (˜ x i + ˜ y i ) = n ( n − (cid:16) λ − ln ΩΩ µ (cid:17) ,
21 + λ − ln ΩΩ µ ! (1 − λ ) X i ˜ x i ˜ y i = 0 . For the anisotropic case λ = 1, the last relation leads to P i ˜ x i ˜ y i = 0. For n = 2vortices, one already sees that ˜ x = − ˜ x and the same for the ˜ y -coordinates.Similarly one can proceed for n > s > ε →
0, (75) and (76) reduce to n X i =1 (˜ x i + ˜ y i ) = 1 + λ n ( n − o (1)and (1 − λ ) n X i =1 ˜ x i ˜ y i = o (1)where o (1) ∼ ln ΩΩ s/ . Remarkably, the distribution of vortices in anharmonictraps with s > Remark:
The analysis in the foregoing sections holds generally for asymptot-ically homogeneous traps according to Def.1.1 in Ref. [22] which is as follows: V is asymptotically homogeneous of order s > U with U ( r ) = 0 for r = 0 such that γ − s V ( γ r ) − U ( r )1 + | U ( r ) | → γ → ∞ and the convergence is uniform in r . U is clearly uniquely determined andhomogeneous of order s , i.e. U ( γ r ) = γ s U ( r ) for all γ ≥ V itself is homogeneous, there is V ≡ U . But if V for instanceis a harmonic-plus-quartic potential, U contains only the quartic contribution.Consider for example the following trap V ( x, y ) = ( x + λ y )[1 + ζ ( x + λ y )] (77)with ζ ∈ (0 ,
1) independent of ε describing the degree of anharmonicity. Thistrap is asymptotically homogeneous of order s = 4. But since V is not homo-geneous, equation (3) is not exactly right. However, using the above definitionfor asymptotically homogeneous potentials, (3) can be used if V ( r ) is replaced31y U ( r ) + o (1). So only the leading contribution, i.e. the asymptotically ho-mogeneous one in the potential is ’visible’. In order to see the anharmoniccontribution of (77) in our regime, one would have to introduce an additionalscaling parameter, i.e. ζ would have to depend on ε (see for instance Ref. [2]). In this paper, we studied the Gross-Pitaevskii (GP) energy and density forBose-Einstein condensates confined in asymptotically homogeneous traps whichare subjected to an external rotation in the Thomas-Fermi (TF) limit whenthe coupling parameter goes to infinity. We derived by analytical estimatesthe leading order of the GP energy and density, which are given by the corre-sponding TF quantities, and the next orders due to vortices. In deriving thecontributions of the vortices, we estimated the relation between the vortexcore sizes and the considered magnitude of angular velocity. As an example,we considered a very general anisotropic homogeneous potential for which wecalculated the critical angular velocities for a finite number of vortices to-gether with the associated GP energy. We have shown that all vortices insidethe Thomas-Fermi domain are single-quantized and arranged in a polygonallattice whose shape can be deduced explicitly by a few simple constraints sat-isfied by the vortex positions. In fact, the results may be used to comparewith experiments when the latter involve asymptotically homogeneous trapsin the TF regime. In this paper, we considered the above trap for the rea-son of explicitness and because it incorporates the harmonic and flat trap forwhich most experimental results are available, but any trap potential satisfy-ing asymptotical homogeneity could be used.
Acknowledgments
The author thanks Jakob Yngvason and Michele Correggi for helpful dis-cussions. Financial support by the Austrian Science Fund FWF under grantP17176-N02 is gratefully acknowledged.
References [1] Aftalion A.,
Vortices in Bose-Einstein condensates , Birkh¨auser (2006)[2] Aftalion A.,Alama S.,Bronsard L., Arch.Rational Mech.Anal. , no.2, 247(2005)[3] Aftalion A., Danaila I., Phys.Rev. A , 023603 (2003)[4] Aftalion A.,Du Q., Phys.Rev.A , 063603 (2001)
5] Andr´e N., Shafrir I., Calc.Var.Part.Diff.Equ. , 191 (1998)[6] Baym G.,Pethick C.J., Phys.Rev.A , 043619 (2004)[7] Bethuel F., Brezis H., Helein F., Ginzburg-Landau vortices , Progress inNonl.Diff. Equ. and their Appl., Vol. 13, Birkh¨auser (1994)[8] Castin Y.,Dum R., Eur.Phys.J. D , 399 (1999)[9] Correggi M., Rindler-Daller T., Yngvason J., J.Math.Phys. , 042104 (2007)[10] Dalibard J., in ’Dynamics and thermodynamics of systems with long-range interactions’, Eds.: Dauxois T., Ruffo S., Arimondo E., Wilkens M.,Lect.Not.Phys., Springer (2002)[11] Fetter A.L., Phys.Rev. A , 063608 (2001)[12] Fetter A.L., Svidzinsky A.A., J.Phys.:Condens.Matter , R135 (2001)[13] Feynman R.P., in Progress in Low Temp.Phys. , 17, ed.C.J.Gorter,Amsterdam, North-Holland (1955)[14] Fischer U., Baym G., Phys.Rev.Lett. , 140402 (2003)[15] Garcia-Ripoll J.J., Perez-Garcia V.M., Phys.Rev. A , 6, 4864 (1999)[16] Garcia-Ripoll J.J., Perez-Garcia V.M., Phys.Rev. A , 041603 (2001)[17] Garcia-Ripoll J.J., Molina-Terriza G., Perez-Garcia V.M., Torner L.,Phys.Rev.Lett. , 14, 140403 (2001)[18] Ignat R., Millot V., J.Funct.Anal. , no.1, 260 (2006)[19] Ignat R., Millot V., Rev.Math.Phys. , no.2, 119 (2006)[20] Kavoulakis G.M.,Baym G., New J.Phys. , 51 (2003)[21] Lieb E.H., Seiringer R., Yngvason J., Phys.Rev.A , 043602 (2000)[22] Lieb E.H., Seiringer R., Yngvason J., Comm.Math.Phys. , 17 (2001)[23] Lieb E.H., Seiringer R., Yngvason J., Phys.Rev. B , 134529 (2002)[24] Lieb E.H., Seiringer R., Comm.Math.Phys. , 505 (2006)[25] Lundh E., Phys.Rev.A, , 043604 (2002)[26] Lundh E., Pethick C.J., Smith H., Phys.Rev. A , 2126 (1997)[27] Madison K.W., Chevy F., Wohlleben W., Dalibard J., Phys.Rev.Lett. , 806(2000)[28] Madison K.W., Chevy F., Wohlleben W., Dalibard J., Mod.Opt. , 2715 (2000)[29] Matthews M.R., Anderson B.P. Haljan P.C., Hall D.S., Wiemann C.E., CornellE.A., Phys.Rev.Lett. , 2498 (1999)
30] Pethick C.J., Smith H.,
Bose-Einstein condensation of dilute gases , CambridgeUniv.Press (2001)[31] Pitaevskii L.P., Stringari S.,
Bose-Einstein condensation , Oxford Science Publ.,Oxford (2003)[32] Rokhsar D.S., Phys.Rev.Lett. ,12, 2164 (1997)[33] Seiringer R., Commun.Math.Phys. , 491 (2002)[34] Seiringer R., J.Phys. A , 9755 (2003)[35] Serfaty S., ESAIM, Cont.Opt.Calc.Var. , 201 (2001)[36] Simula T.P., Virtanen M.M., Salomaa M.M., Phys.Rev. A , 033614 (2002)[37] Watanabe G., Gifford S.A., Baym G., Pethick C.J., Phys.Rev. A , 063621(2006), 063621(2006)