Quantum Fluctuations of Vortex-Lattice State in Ultrafast Rotating Bose Gas Qiong Li, Bo Feng, and Dingping Li Department of Physics, Peking University, Beijing, 100871, China Quantum fluctuations in an ultrafast rotating Bose gas at zero temperature are investigated. We calculate the condensate density perturbatively to show that no condensate is present in the thermodynamic limit. The excitation from Gaussian fluctuations around the mean field solution causes infrared divergences in loop diagrams, nevertheless,in calculating theatom numberdensity, 1 the correlation functions and the free energy, we find the sum of the divergences in the same loop 1 ordervanishesand obtain finitephysical quantities. Thelong-range correlation is explored and the 0 algebraic decay exponentfor thesingle-particle correlation function is obtained. The atom number 2 density distribution is obtained at the one-loop level, which illustrates the quantum fluctuation n effectstomelt themeanfieldvortex-lattice. Bythenon-perturbativeGaussian variationalmethod, a we locate thespinodal point of the vortex-latticestate. J 4 PACSnumbers: 03.75.Hh,03.75.Lm,05.30.Jp,05.30.Rt 1 ] I. INTRODUCTION s a g Theappearanceofvortexexcitationsinresponsetorotationisacharacteristicfeatureofsuperfluid[1–3]. Sincethe - t discoveryofBose-Einsteincondensation(BEC)inatomicgases[4–6],muchworkhasbeendevotedtothepropertiesof n rotatinggaseouscondensatesintraps,andthesedevelopmentshavebeenreviewedin[7,8]. Astherotationfrequency a Ω increases, more and more vortices occur, from a single one to several ones, then they form an Abrikosov lattice, u q i.e., a triangular array with a surface density nυ = mπ~Ω(m is the mass of the condensed atoms) [9–18]. In the frame . co-rotating with harmonic traps, rotation effects are described by a centrifugal force, which effectively reduces the t a transverse harmonic trapping, and a Coriolis force which has the same mathematical structure as the Lorentz force m that an electron experiences in a uniform magnetic field. In the fast rotation regime, when the centrifugal force - almost cancels the transverse confinement, the energy levels of the single particle part of the Hamiltonian organize d into Landau levels with spacing 2~Ω, and the condensates expand radically, leading to a very dilute atom number n density, which ensures the mean interaction energy smaller than Landau level spacing, so that the cold atoms are o confined in the lowest Landau level (LLL) single particle orbit. c [ For rotating bosonic atoms in the LLL, the filling fraction, i.e., the ratio of the number of atoms to the number of vortices,isthe parametercontrollingthenatureofthe system. Athighfillingfractions,the condensateisinthemean 1 field quantum-Hall regime and forms an ordered vortex lattice ground state [19–24]. As the filling fraction decreases, v there is a zero temperature phase transition from a triangular vortex-lattice to strongly correlated vortex-liquid[25] 7 4 and the melting point is located by various approaches [25–27] to be approximately at the filling fraction 6 10. ∼ 7 The experiment done by V. Schweikhard etc. [28] created an ordered vortex lattice in the mean-field quantum-Hall 2 regime, and provided evidence that the elastic shear strength of the vortex lattice drops substantially as the BEC . enters the mean field quantum-Hall regime. 1 0 In the ultrafast rotation limit when the transverse confinement is exactly canceled by the centrifugal force, the 1 condensates expand to be a two dimensional configuration and atoms are frozen in the lowest energy level in the z 1 direction (assuming a strong confinement in the z direction). For such a two dimensional system, J. Sinova etc. [26] v: findthatthe solutiontothe Gross-Pitaevskii(GP)equationisanAbrikosovtriangularvortex-latticeandtheintegral i for the fraction of atoms outside the condensate diverges logarithmically with the system size, which implies that no X BEC occurs in the thermodynamic limit even at zero temperature. r Fluctuation effects in the transverse plane in an ultrafast rotation limit are important and one has to go beyond a meanfieldtreatment[19]. Inthepresentpaper,wearedevotedtoquantumfluctuationsintheultrafastrotationlimit and focus on the thermodynamic limit at zero temperature. Basedon analyticalcalculations in the perturbative framework,we show that the condensate density is zero in the thermodynamic limit, namely, there is no BEC. The single-particle correlation function and the density fluctuation correlation function are shown to fall off as an inverse power of the separation distance in the large distance limit, whichindicatesanalgebraiclong-rangeorderandthealgebraicdecayexponentisobtained. Theatomnumberdensity distribution is obtainedat the one-looplevel, which illustrates the quantumfluctuation effects to melt the mean field vortex-lattice. By loop expansion around the mean field solution, we calculate the free energy density up to two-loop. The mean field solution is an Abrikosov triangular vortex-lattice and the excitation from Gaussian fluctuations has a quadratic dispersionatsmallwavevectors[26]. Wefindthatthequadraticdispersioncausesinfrareddivergencesinthetwo-loop 2 diagrams, nevertheless, the sum of the divergences vanishes and the two-loop contribution to the free energy density is finite. We also study the model by non-perturbative Gaussian variational method. The free energy density calculated by the perturbationtheoryandthatbythe Gaussianvariationalmethodcoincide verywellatlargefilling fractions. The vortex lattice solution exists only when the filling fraction ν is greater than a certain value ν , which is found to be s about 1.1. The point ν =ν is the so called spinodal point. Between the spinodal point and the melting point is the s meta-stable vortex-lattice state. The paper is organized as follows: In section II, we formulate the model. In section III, we explore the long-range correlations. In section IV, we calculate the free energy density perturbatively up to two-loop. In section V, we use the non-perturbative Gaussian variational method to study the model. In section VI, we give a summary and the conclusions. II. THE MODEL ForasystemofN bosonicatomsinanaxisymmetricharmonictrap(withtrapfrequenciesω andω )rotatingwith z angular velocity Ωe , the Hamiltonian in the rotating frame is ⊥ z N (p mΩzˆ r )2 1 1 N H = [ i− × i + m(ω2 Ω2)(x2+y2)+ mω2z2]+g δ(r r ), (1) 2m 2 ⊥− i i 2 z i i− j i=1 i<j=1 X X where g = 4π~m2as is the strength of the hard core repulsive interactions, with as the s-wave scattering length. The centrifugal force effectively reduces the radial confinement and the Coriolis force is equivalent to the Lorentz force exerted on a particle with charge Q by a magnetic field B = 2mΩe . In the present paper we are confined to the Q z ultrafast rotation limit by setting Ω = ω , and assume the axial confinement is so strong that atoms are frozen in the lowest harmonic state in the z directi⊥on. Consequently, what we concern is essentially a two dimensional system of charged bosonic atoms experiencing an effective magnetic field in the z direction, described by the Hamiltonian N (p mΩzˆ r )2 N H = i− × i +g δ(r r ). (2) i j 2m − i=1 i<j=1 X X The kinetic part of the Hamiltonian has equally spaced Landau levels (with spacing 2~Ω) and the interaction part is a small perturbation, ng 2~Ω (n is the mean number density of the atoms). Without thermal fluctuations at zero ≪ temperature,allatomsareconfinedintheLLL.Inthe LLLsubspace,thekinetic partoftheHamiltonianisquenched and represented by ~Ω, hence we have the grand Hamiltonian in second quantized form (with the unit ~=1) g Hˆ µNˆ = d2r (Ω µ)Ψ (r)Ψ(r)+ Ψ (r)Ψ (r)Ψ(r)Ψ(r) , (3) † † † − − 2 Z h i and the grand-canonicalpartition function in the functional formalism (β,µ)= [Ψ Ψ]e S[Ψ∗,Ψ], (4) ∗ − Z D Z where the action S[Ψ ,Ψ] takes the form ∗ β 1 dτ d2r Ψ∗(r,τ)(∂τ µ+Ω)Ψ(r,τ)+ g Ψ(r,τ)4 , (5) − 2 | | Z0 Z (cid:20) (cid:21) in which β = 1 and µ is the chemical potential. By variable rescaling : 1 (µ Ω) = a ,τ = 1 τ ,β = kBT mgΩ − µ mgΩ ′ 1 β ,r= 1 r,Ψ=√2mΩΨ, the partition function simplifies mgΩ ′ √2mΩ ′ ′ β′ ′(β′,aµ)= [Ψ′∗Ψ′]exp dτ′ d2r′ Ψ′∗(r′,τ′)(∂τ′ aµ)Ψ′(r′,τ′)+ Ψ′(r′,τ′)4 . (6) Z D − − | | Z Z0 Z (cid:2) (cid:3) With all primes omitted, we have 3 β (β,a )= [Ψ Ψ]exp dτ d2r Ψ (r,τ)(∂ a )Ψ(r,τ)+ Ψ(r,τ)4 . (7) µ ∗ ∗ τ µ Z D − − | | Z Z0 Z (cid:2) (cid:3) Note that the kinetic part is absorbed in a shift of the effective chemical potential a , and only the interaction part µ is relevant. Hereafter we will work with the rescaled variables. Using Landau gauge A = ( By,0) for the effective magnetic field B = 2mΩe , we have the LLL magnetic Bloch − Q z representation [29, 30] 2 ϕk(r)=318 ∞ exp 2π y √3n √3kxd +i π(n2 n)+2πnx + √3nkyd+kxx . (8) −√3 d − 2 − 4π ! 2 − d 2 ! n= X−∞ ϕ(r), i.e. ϕk(r) with k=0, is a superposition of the lowest Landau levels of the kinetic part of the Hamiltonian and corresponds to an Abrikosov triangular vortex-lattice with lattice spacing d. In the rescaling length unit, d equals √4π3 and the primitive vectors of the vortex-lattice are d1 = √4π3(1,0) , d2 = √4π3(12,√23). ϕk(r) describes the tqransverseoscillationsofthevortex-lattice,andtheprimitivevecqtorsofthereciprocqallatticeared˜ = 4π(√3, 1), 1 √3 2 −2 d˜2 = √4π3(0,1). The vortex cores, i.e. the zero point of the function ϕk(r) with k = k1d˜1 +k2d˜2,qare uniformly distribquted at sites (n k )d +(n + 1 +k )d , where n and n are integers. Note that both the lattice cell and 1− 2 1 2 2 1 2 1 2 the Brillouin zone have an area of 2π, and the vortex number density is n = 1 in the rescaling units. υ 2π III. THE LONG-RANGE CORRELATIONS As is known [26], for the ultrafast rotating two dimensional Bose gas, there is no BEC in the thermodynamic limit even at zero temperature. In spite of that, we will calculate U(1) invariant quantities, like the atom number density, the free energy and the correlation functions, in the perturbative framework, and show that the infrared divergences are canceled, similar to the method used in two dimensional non-linear σ model [31, 32]. In two dimensional O(N) non-linear σ model, David in [32] proved that, using the “wrong” spontaneously broken symmetry phase, any O(N) invariant observable has an infrared finite weak coupling perturbative expansion. In this paper, though we will only show some U(1) invariant quantities are also free of infrared divergences at most to two loops in the perturbative framework, we believe that it is true to all orders similar to the two dimensional O(N) non-linear σ model. This method was extensively used to study the vortex lattice in type II superconductors, for example, in Ref. [33, 34]. Inthis section,We calculatethe condensatedensity perturbatively andshowitis zerointhe thermodynamiclimit. By calculating the single-particle correlationfunction and the density fluctuation correlationfunction, we obtain the algebraicdecayexponent. Theatomnumberlocaldensityisalsocalculated,whichshowsthatatlargefillingfractions, the numberdensityretainsthe vortex-latticeconfiguration,while atsmallfilling fractions,quantumfluctuations tend to smooth away the vortex-lattice. A. The atom number density In a usual fashion [35–38], we separate the field as the condensate part and the fluctuation part Ψ(r,τ)=√n ϕ(r)+ψ(r,τ), (9) c wheren ,thecondensatenumberdensity,isarealnumberminimizingthefreeenergyandthefluctuationpartψ(r,τ) c can be expanded as 1 ψ(r,τ)= ψkmϕk(r)e−2iθke−iωmτ . (10) √Aβ k BZ m X∈ X In powers of ψk∗m and ψkm, we divide the action S into four parts 4 S = βA a n +β n2 , 0 − µ c A c S2 = (cid:0)(−iωm−aµ+(cid:1)4ncβk)ψk∗mψkm+nc|γk| ψk∗mψ−∗k−m+ψkmψ−k−m , p X(cid:2) (cid:0) (cid:1)(cid:3) 1 S3 = 2√nc√Aβ ψp∗1ψp∗2ψp3Pp1p2p30+c.c. , (11) p1X,p2,p3(cid:0) (cid:1) 1 S = ψ ψ ψ ψ P , 4 Aβ p∗1 p∗2 p3 p4 p1p2p3p4 p1,pX2,p3,p4 where p≡(k,ωm) , ψp ≡ψkm and p ≡ k BZ m . The quadratic part can be diagonalized as ∈ P P P S2 = (−iωm+ǫ(k))a∗kmakm, (12) p X by the following Bogoliubov transformation: akm = ukψkm+υkψ∗k m, (13) − − a∗km = ukψk∗m+υkψ k m, (14) − − and inversely, ψkm = ukakm−υka∗−k−m, (15) ψk∗m = uka∗km−υka−k−m, (16) where 1 ǫ (k) 0 uk = +1 , s2 ǫ(k) (cid:18) (cid:19) 1 ǫ (k) 0 υk = 1 , (17) s2 ǫ(k) − (cid:18) (cid:19) and ǫ0(k) = ( aµ+4ncβk), − ǫ(k) = (−aµ+4ncβk)2−4n2c|γk|2. (18) q Note that the free energy density depends on two parameters, a and n , nevertheless, n is related to a by the µ c c µ F constraint ∂F(∂nncc,aµ) =0 and hence nc is renormalized order by order. In the zero-looporder, = a n +β n2, (19) Fc − µ c A c and the condensate density equals a n(0) = µ , (20) c 2β A which is also the total number density in the zero-loop order, denoted as n . Now we want to calculate the total 0 number density and the condensate density to one-loop. The total number density n is given by 1 1 A d2r Ψ†(r)Ψ(r) =nc+ Aβ hψk∗mψkmi (21) Z p (cid:10) (cid:11) X 5 where nc is the condensate density and A1β phψk∗mψkmi is the density of the atoms outside the condensate. First we need to know the one-loop correction to the condensate density, denoted as n(1). Up to one-loop order, the free P c energy density takes the form 1 d2k F0+1(aµ,nc)=−aµnc+βAn2c + 4π 2π (−aµ+4ncβk)2−4n2c|γk|2. (22) ZBZ q Minimizing with respect to n leads to 0+1 c F aµ 1 d2k βk( aµ+4ncβk) nc γk 2 n = − − | | , (23) c 2βA − 2πβA ZBZ 2π (−aµ+4ncβk)2−4n2c|γk|2 q from which we see that the one-loop correction to n equals c n(1) = 1 d2k2βk(2βk−βA)−|γk|2 . (24) c −4πβA ZBZ 2π (2βk βA)2 γk 2 − −| | q The density of the atoms outside the condensate is given by 1 1 d2kE (k) 0 Aβ hψk∗mψkmi1 loop = 4π 2π E(k) , (25) p − ZBZ X where E (k) and E(k) are defined as 0 E0(k) = 2βk βA, (26) − E(k) = (2βk βA)2 γk 2. (27) − −| | q The total number density, n0+1, is equalto n(c0)+n(c1)+ A1β phψk∗mψkmi1 loop . The filling fraction, n/nυ, is given by 2πn at the one-loop level. − 0+1 P Obviously, n(c1) and A1β phψk∗mψkmi1 loop both contain infrared divergences, but in the sum the divergences are canceled and n(c1)+ A1β Pphψk∗mψkmi1 lo−op , which is the one-loop correction to the total number density, equals − P 1 d2k n = E(k) 0.023. (28) 1 −4πβ 2π ≃− A ZBZ (0) (1) (0) (1) The condensate density, n =n +n , is infrareddivergent,n +n =c(1 αlnL), where c and α are positive c c c c c constants and L is the system size. If we can calculate n to all loops, n −c 1 αlnL+ 1(αlnL)2+ cexp( αlnL) = cL α. When we take L , n 0. Tcherefore, in thecth≃ermod−ynamic lim2it, there will··b·e n≃o − c conden−sate, as is also shown by J. Sinova et→c.∞[26]. F→or a finite system, there is a fi(cid:0)nite infrared cutoff 1, and(cid:1)the condensate density, n =cL α , will be finite. ∼ L c − Up to one-loop, the local density n(r) is given by 1 d2kE (k) Ψ†(r)Ψ(r) =(n(c0)+n(c1))ϕ∗(r)ϕ(r)+ 4π 2π E0(k) ϕ∗k(r)ϕk(r), (29) ZBZ (cid:10) (cid:11) where n(0)ϕ (r)ϕ(r) is the mean field local density, denoted as n (r). Obviously, in Eq. (29) n(1)and the last term c ∗ 0 c both contain infrared divergences, nevertheless, by arranging them properly we see the divergences are canceled and the one-loop correction to the total local density is obtained as 1 d2kE (k) n1(r) = n(c1)ϕ∗(r)ϕ(r)+ 4π 2π E0(k) ϕ∗k(r)ϕk(r) ZBZ 1 d2kE (k) = n1ϕ∗(r)ϕ(r)+ 4π 2π E0(k) (ϕ∗k(r)ϕk(r)−ϕ∗(r)ϕ(r)), (30) ZBZ in which the second term is free of divergences. 6 0 0.4 0.8 1.2 1.5 |ϕ(r)|2 ν = 10 ν = 6 ν = 2 ν = 1.5 FIG. 1: The mean field atom numberdensity distribution forms a perfect triangular lattice, described by thefunction |ϕ(r)|2. Quantumfluctuationstendtosmooththemeanfieldvortex-lattice. Asthefillingfractionlowers,quantumfluctuationsincrease and the vortex-latticebecomes more smooth. Note that the mean field local density, n (r)=n(0)ϕ (r)ϕ(r) , forms a vortex-lattice, and the one-loop correction, 0 c ∗ n (r) , includes quantum fluctuations. In order to explore quantum fluctuation effects on the atom number density 1 configuration, we plot the mean field density distribution, n0(r) = ϕ (r)ϕ(r), and the total density distribution, n0 ∗ n0(r)+n1(r),atdifferent filling fractions. FromFig. 1, onefinds thatquantumfluctuations tend to smooththe vortex- n0+n1 lattice. Atlargefillingfractions,the numberdensitystillretainsthe vortex-latticeconfiguration,whileatsmallfilling fractions, the vortex-lattice is smoothed away. B. The single-particle correlation function Up to one-loop, the single-particle correlation function Ψ (r )Ψ(r ) is equal to † 1 2 (cid:10) (cid:11) Ψ†(r1)Ψ(r2) = n(c0)+n(c1) ϕ∗(r1)ϕ(r2)+ ψ†(r1)ψ(r2) (cid:10) (cid:11) = (cid:16)n(0)+λn(1(cid:17)) ϕ (r )ϕ(r )+(cid:10)λ ψ (r )ψ((cid:11)r ) , (31) c c ∗ 1 2 † 1 2 (cid:16) (cid:17) (cid:10) (cid:11) where λ is used to keep trace of the loop order and should be set to 1 in the end. Similar to Eq. (29), the above equation can be arranged as 1 d2kE (k) n(0)+λn(1)+λ 0 ϕ (r )ϕ(r )+λf(r ,r ) c c 4π 2π E(k) ∗ 1 2 1 2 (cid:18) ZBZ (cid:19) = n(0)+λn ϕ (r )ϕ(r )+λf(r ,r ), (32) c 1 ∗ 1 2 1 2 (cid:16) (cid:17) where f(r1,r2) is given by41π BZ d22πkEE0((kk))[ϕ∗k(r1)ϕk(r2)−ϕ∗(r1)ϕ(r2)] and free of divergences. Now we define the relative coordinate r r r , and the center coordinate R r1+r2. We take R=0 for R ≡ 1 − 2 ≡ 2 simplicity and analyze the large r limit. Based on numerical calculations, we find | | lim f(r ,r ) 0.18ln rϕ (r )ϕ(r ) (33) 1 2 ∗ 1 2 r ≃− | | | |→∞ Therefore lim Ψ (r )Ψ(r ) n(0)+λn ϕ (r )ϕ(r ) 0.18λln rϕ (r )ϕ(r ) r † 1 2 ≃ c 1 ∗ 1 2 − | | ∗ 1 2 | |→∞(cid:10) (cid:11) (cid:16) (cid:17) = n(0)+λn ϕ (r )ϕ(r )(1 αln r) (34) c 1 ∗ 1 2 − | | (cid:16) (cid:17) where α = 0.18λ , and to one loop order α 0.18λ = 0.18. Of course this result is not physical if we only include n(c0)+λn1 ≈ n(c0) n(c0) one-loop correction, as limr Ψ†(r1))Ψ(r2) . We shall include all order in perturbation theory to get the physical result. We argue, i→n∞all order, similar t→o th−e∞calculation in Ref. [34], (cid:10) (cid:11) 7 1 rlim Ψ†(r1)Ψ(r2) ∝ ϕ∗(r1)ϕ(r2) 1−α′ln|r|+ 2(α′ln|r|)2+... (35) →∞ (cid:18) (cid:19) (cid:10) (cid:11) = ϕ∗(r1)ϕ(r2) r−α′ | | To one loop, α = α = 0.18. similar calculations can be done in usual BKT (Berezinsky, Kosterlitz and Thouless) ′ n(c0) phase transition systems, and the algebraic decay exponent can be obtained correctly [39, 40]. One will wonder in usual BKT phase transition systems, the phase transition is continuous. In the vortex lattice phase, from Eq. (35), the correlation decays algebraically, and the rotational symmetry is broken as the factor ϕ (r )ϕ(r ) in Eq. (35) is notrotationallyinvariant. In the vortexliquid phase, the correlationdecays exponentially, ∗ 1 2 andtherotationalsymmetryisunbroken. Therefore,thephasetransitionisthespontaneousbreakingoftherotational symmetry. As in usual solid to liquid phase transition, the phase transition is a first order melting transition. C. The density fluctuation correlation function In the following, we will calculate the density fluctuation correlation to one-loop, δnˆ(r )δnˆ(r ) = [nˆ(r ) nˆ(r ) ][nˆ(r ) nˆ(r ) ] 1 2 1 1 2 2 h i h −h i −h ii = nˆ(r )nˆ(r ) nˆ(r ) nˆ(r ) . (36) 1 2 1 2 h i−h ih i As shown in subsection IIIA, there is no condensate and thus the system retains the U(1) symmetry, therefore, only the U(1) gauge invariant terms in the contractions remain nˆ(r )nˆ(r ) nˆ(r ) nˆ(r ) 1 2 1 2 h i−h ih i = Ψ (r )Ψ(r ) Ψ(r )Ψ (r ) + Ψ (r )Ψ(r )Ψ (r )Ψ(r ) . (37) † 1 2 1 † 2 † 1 1 † 2 2 c (cid:10) (cid:11)(cid:10) (cid:11) (cid:10) (cid:11) The term Ψ (r )Ψ(r ) is investigated in the last subsection. In the larger r limit, the first term in Eq. (37) † 1 2 1 2 | − | falls off as an inverse power of the separation distance r r , 1 2 (cid:10) (cid:11) | − | lim Ψ†(r1)Ψ(r2) Ψ(r1)Ψ†(r2) ϕ(r1)2 ϕ(r2)2 r1 r2 −2α. (38) r r ∼| | | | | − | | 1− 2|→∞ (cid:10) (cid:11)(cid:10) (cid:11) ThesecondterminEq. (37)isthetwo-bodyconnectedGreen’sfunctionanditishardtocalculatenon-perturbatively. But we speculate that this term will not alter the asymptotic behavior of δnˆ(r )δnˆ(r ) in the larger r limit, 1 2 1 2 h i | − | andhencethe densityfluctuationcorrelationfunctionalsodecaysalgebraicallywiththe distanceinthe largedistance limit. Theresultwillindicate,neartheBragpeak,thestructurefunctionS(Q+k)whereQbelongstothereciprocal lattice and k is small will have a scaling limk 0S(Q+k) k−(2−2α). Similarcalculationcanbe done forvortexl→attice intype∝II|su|perconductorsandthe densityfluctuationcorrelation will be shown to have similar behavior. For the vortex lattice in type II superconductors near the Brag peak, the structure function S(Q+k) will have a scaling limk 0S(Q+k) k−η (details will be published elsewhere). → ∝| | IV. THE FREE ENERGY DENSITY CALCULATION BY LOOP EXPANSION Inthissection,weshallcalculatethefreeenergydensitybyloopexpansionuptotwo-loopandshowthecancelation of infrared divergences. A. Mean-field contribution In the saddle point approximation, δS[Ψ ,Ψ] δS[Ψ ,Ψ] ∗ ∗ =0, =0, (39) δΨ δΨ ∗ 8 with the LLL constraint, one obtains [3, 41] a Ψ (r)= µ ϕ(r), (40) 0 2β A r whereϕ(r)=ϕk=0(r) , βA = 21π celld2r|ϕ(r)|4. Obviously,the saddlepointapproximationisequivalentto the mean field GP equation. We have the mean field contribution to the free energy density R a2 µ (a )= , (41) 0 µ F −4β A and the mean field contribution to the number density ∂ (a ) a c µ µ n (a )= F = . (42) 0 µ − ∂a 2β µ A B. One-loop correction Following the loop expansion procedure presented in [42], we set a Ψ(r,τ)= µ ϕ(r)+ψ(r,τ), (43) 2β A r in which aµ ϕ(r) is the mean field part and ψ(r,τ) is the higher order corrections, and then expand S[Ψ ,Ψ] in 2βA ∗ powers ofqψ∗(r,τ) and ψ(r,τ). In the magnetic Bloch representation, ψ(r,τ) can be expanded as 1 ψ(r,τ)= ψkmϕk(r)e−2iθke−iωmτ (44) √Aβ k BZ m X∈ X where A is the area of the sample, ωm = 2πβm is the bosonic Matsubara frequency, and θk is defined in Eq. (46). In powers of ψk∗m and ψkm, we divide the action into four parts a2 µ S = Aβ , 0 −4βA! a 1 a µ µ S2 = p [(cid:18)−iωm+ βA (2βk−βA)(cid:19)ψk∗mψkm+ 2βA|γk| ψk∗mψ−∗k−m+ψkmψ−k−m ], X (cid:0) (cid:1) 2a µ S3 = β Aβ ψp∗1ψp∗2ψp3Pp1p2p30+c.c. , (45) r A p1X,p2,p3(cid:0) (cid:1) 1 S4 = Aβ ψp∗1ψp∗2ψp3ψp4Pp1p2p3p4, p1,pX2,p3,p4 where p≡(k,ωm) , ψp ≡ψkm , p ≡ k BZ m and ∈ 1P P P βk = d2rϕ∗k(r)ϕ∗(r)ϕk(r)ϕ(r), 2π Zcell 1 γk = d2rϕ∗(r)ϕ∗(r)ϕk(r)ϕ k(r), 2π Zcell − eiθk = γk , (46) γk | | d2r Pp1p2p3p4 = δm1+m2,m3+m4 2π ϕ∗k1(r)ϕ∗k2(r)ϕk3(r)ϕk4(r)e2i(θk1+θk2−θk3−θk4). Zcell In the one-loop approximation, we keep only the quadratic part of the action and diagonalize it as 9 S2 = (−iωm+ǫ(k))a∗kmakm, (47) p X where a ǫ(k)= µ (2βk βA)2 γk 2, (48) β − −| | A q by the following Bogoliubov transformation: akm = ukψkm+υkψ∗k m, (49) − − a∗km = ukψk∗m+υkψ k m, (50) − − and inversely, ψkm = ukakm−υka∗−k−m, (51) ψk∗m = uka∗km−υka−k−m, (52) where uk = 12 EE0((kk)) +1 ,υk = 12 EE0((kk)) −1 , E0(k) and E(k) are defined in Eqs. (26) and (27). By Taylor r r expandingβk and(cid:16)γk ,wefi(cid:17)ndtheexcit(cid:16)ationǫ(k)(cid:17)hasaquadraticdispersionatsmallwavevectors,i.e. limk 0ǫ(k) k2, whichis consis|ten|twithpreviousresults [26,43,44]. The one-loopcontributionto the freeenergydensit→y, (a ∼), 1 µ F takes the form 1 1 −Aβ ln D[a∗a]exp− (−iωm+ǫ(k))a∗kmakm Z p X 1 1 1 = ǫ(k)+ ln 1 e βǫk . (53) − A 2 β − k BZ(cid:20) (cid:21) X∈ (cid:0) (cid:1) By setting the area A to infinity and the temperature T to zero, we obtain a (a )= µ E(k) , (54) 1 µ k F 4πβ h i A where d2k, means average over the Brillouin zone. The one-loop correction to the number density, h···ik ≡ BZ 2π n (a ), is equal to 1 µ R ∂ (a ) 1 F1 µ = E(k) , (55) k − ∂a −4πβ h i µ A which is consistent with the result obtained in the last section as in Eq. (28). C. Two-loop correction By the Bogoliubov transformation shown in Eqs. (51) and (52), we switch to the field a∗km, akm, and write the cubic and quartic part as 2a 1 µ S3 =rβA √Aβ p1X,p2,p3[ap1ap2a∗p3(Λp1p2p3 −Λ′p1p2p3)+ap1ap2ap3(Q′p1p2p3 −Qp1p2p3)]+c.c. (56) and 1 S4 = {Aβ [ap1ap2ap3ap4Pp1p2−p3−p4υk1υk2uk3uk4 p1,pX2,p3,p4 −2a∗p1ap2ap3ap4(Pp1−p4p2p3uk1uk2uk3υk4 +Pp2p3p1−p4υk1υk2υk3uk4)]+c.c.} 1 +Aβ a∗p1a∗p2ap3ap4[4Pp1−p4p3−p2uk1uk3υk2υk4 (57) p1,pX2,p3,p4 +Pp1p2p3p4uk1uk2uk3uk4 +Pp3p4p1p2υk1υk2υk3υk4], 10 FIG. 2: The two-loop Feynman diagrams. where ap akm and ≡ 1 p1p2p3 =3(P0−p1p2p3uk2uk3υk1 +P0−p2p3p1uk1uk3υk2 +P0−p3p2p1uk1uk2υk3), Q 1 ′p1p2p3 =3(Pp3p2−p10υk2υk3uk1 +Pp1p3−p20υk1υk3uk2 +Pp1p2−p30υk1υk2uk3), (58) and Q Λp1p2p3 = P0p1p3−p2υk1υk3uk2 +P0p2p3−p1υk2υk3uk1 +P0p3p2p1uk1uk2uk3, Λ′p1p2p3 = P−p2p3p10uk1uk3υk2 +P−p1p3p20uk2uk3υk1 +Pp1p2p30υk1υk2υk3. (59) The two-loop contribution to the free energy density, (a ), takes the form 2 µ F 1 [a ,a]exp (S +S +S ) 1 1 ∗ 2 3 4 ln D − = S S S , (60) 4 3 3 − βA [a ,a]exp S βA h i− 2h i (cid:20) R D ∗ − 2 (cid:21)2 loop (cid:18) (cid:19) − where denotes the sum of all thRe connected Feynman diagrams with G = 1 as a propagator. The two-loohp··F·eiynman diagrams are depicted in Fig. 2. The contribution from thepdiag−raiωmm+“ǫ(k”) equals ∞ 1 β2A2 4Pp1−p2p1−p2u2k1υk22 +4Pp1−p1p2−p2uk1υk1uk2υk2 +2Pp1p2p1p2 u2k1u2k2 +υk21υk22 Gp1Gp2 (61) pX1,p2(cid:2) (cid:0) (cid:1)(cid:3) 1 1 = 4π2βA (h|γk|ukυkik)2+ 8π2 βk1−k2 u2k1 +υk21 u2k2 +υk22 k1,k2; (62) (cid:10) (cid:0) (cid:1)(cid:0) (cid:1)(cid:11) the contribution from the diagram “ ” equals ⊖ 12a 2 µ −βAβ2A2 ′p1p2p3 − p1p2p3 Gp1Gp2Gp3 3 p1X,p2,p3(cid:12)(cid:12)(cid:12)Q Q (cid:12)(cid:12)(cid:12) 2 1 = −π2 < ′k1k2h−k1−k2i− k1k2h−k1−k2i E(k1)+E(k2)+E( k1+k2 ) >k1,k2; (63) (cid:12) (cid:12) h i and the contribution from t(cid:12)hQe diagram “ - Q”equals (cid:12) (cid:12) (cid:12) (cid:13)(cid:13) 8a µ Λ Λ Λ Λ G G G −β β2A2 p1p2p1 − ′p1p2p1 ∗p3p2p3 − ′p∗3p2p3 p1 p2 p3 A p1X,p2,p3(cid:0) (cid:1)(cid:0) (cid:1) = −4π21β βk(u2k+υk2)−|γk|ukυk k 2. (64) A (cid:0)(cid:10)(cid:2) (cid:3)(cid:11) (cid:1) We have set the temperature T to zero and the area A to infinity in the end. The notation h···ik1,k2 ≡ d2k1 d2k2, means average over the Brillouin zone, and k +k represents the reduced wave vector in the BZ 2π BZ 2π h 1 2i Brillouin zone. Obviously, ∂ (a ) = 0, and hence the two-loop correction to the atom number density is zero, R R ∂aµF2 µ n (a ) = 0. By Taylor expansion, one can confirm that each of the three contributions in Eqs. (62), (63) and (64) 2 µ has infrared divergences. However, all the divergences are exactly canceled if they are summed up. By numerical integration, we find that (a )=0 (65) 2 µ F It is amazing that it precisely amounts to zero, and the reason is still under investigations.