劉 燕, 張素英（山西大學理論物理研究所,山西太原 030006）

環(huán)形勢阱中旋轉(zhuǎn)玻色愛因斯坦凝聚體的基態(tài)

劉 燕, 張素英
（山西大學理論物理研究所,山西太原 030006）

應(yīng)用托馬斯－費米近似和虛時演化數(shù)值方法研究環(huán)形勢阱中旋轉(zhuǎn)玻色愛因斯坦凝聚體的基態(tài)密度分布.當增加其旋轉(zhuǎn)角頻率,或者增加環(huán)形勢阱的寬度及相應(yīng)的中心高度,凝聚體基態(tài)密度分布均從渦旋晶格相轉(zhuǎn)變?yōu)榫逌u旋相.當旋轉(zhuǎn)角頻率為零時,增加環(huán)形勢阱的寬度及相應(yīng)的中心高度,凝聚體基態(tài)密度分布從一個圓盤變?yōu)閳A環(huán).解析結(jié)果與數(shù)值結(jié)果相互吻合.

托馬斯－費米近似；中心洞；巨渦旋

0 引言

近年來,有關(guān)冷原子研究的快速發(fā)展以及實驗技術(shù)的不斷提高實現(xiàn)了對各種不同形式的束縛外勢的調(diào)控,如諧振子勢、光晶格勢以及環(huán)形勢等[1－2].對不同束縛外勢阱中的旋轉(zhuǎn)玻色愛因斯坦凝聚體（BEC）的理論研究也取得了很大的進展.對于諧振子勢阱中的旋轉(zhuǎn)玻色愛因斯坦凝聚體,當旋轉(zhuǎn)角頻率Ω趨近于或者大于勢阱束縛頻率ω⊥時,凝聚體的托馬斯－費米近似半徑發(fā)散[3－4].若在諧振子勢的基礎(chǔ)上加入一個四次勢,則可以克服這一缺陷使托馬斯－費米近似半徑在Ω＞ω⊥時繼續(xù)收斂,且隨著旋轉(zhuǎn)角頻率的增大,凝聚體的基態(tài)密度分布呈現(xiàn)出各種新奇的渦旋拓撲結(jié)構(gòu)[5－10].改變凝聚體粒子間的相互作用強度,拓撲結(jié)構(gòu)也會發(fā)生相應(yīng)的變化[11－17].除此之外,凝聚體的基態(tài)密度分布還和束縛勢阱的結(jié)構(gòu)有關(guān).對于旋轉(zhuǎn)環(huán)形勢阱中的玻色愛因斯坦凝聚體,當粒子間互作用強度和勢阱旋轉(zhuǎn)角頻率一定時,調(diào)節(jié)環(huán)形勢阱的寬度及相應(yīng)的中心高度,凝聚體的基態(tài)密度分布必然也會出現(xiàn)新的渦旋相.我們考慮單分量玻色愛因斯坦凝聚體束縛在如下形式的環(huán)形勢阱中且勢阱以角速度Ω＝Ωz繞z軸旋轉(zhuǎn)：

勢阱的最低點坐標為（±a0,0）,中心最高點坐標為（0,/2）,其中a0＝,r表示徑向坐標,V0,r0是無量綱化的常數(shù),可以調(diào)整勢阱的寬度及相應(yīng)的中心高度.我們首先應(yīng)用托馬斯－費米近似理論求解上述環(huán)形勢阱中的玻色愛因斯坦凝聚體在如下三種情況下的基態(tài)密度分布情況：①勢阱的寬度及相應(yīng)的中心高度一定,增大其旋轉(zhuǎn)角頻率；②勢阱的旋轉(zhuǎn)角頻率一定,增大其寬度及相應(yīng)的中心高度；③令勢阱旋轉(zhuǎn)角頻率為零,增大其寬度及相應(yīng)的中心高度.進一步運用虛時演化數(shù)值方法模擬這三種情況下凝聚體的基態(tài)密度分布情況,并將其同解析結(jié)果進行比較.

此外,我們還可以保持旋轉(zhuǎn)角頻率及勢阱寬度不變,數(shù)值模擬凝聚體的基態(tài)密度分布隨勢阱中心高度增加的變化.但是,如何解析求解量子化渦旋數(shù)目與旋轉(zhuǎn)角頻率及勢阱結(jié)構(gòu)的關(guān)系還需進一步研究.

1 理論模型

絕對零度下受外勢束縛的玻色愛因斯坦凝聚體可以用一個波函數(shù)ψ來描述,它滿足如下非線性薛定諤方程

其中as是一個正的s波散射長度,描述凝聚體粒子間的排斥相互作用.如果外勢阱繞z軸以Ω ＝Ωz的角速度旋轉(zhuǎn),則波函數(shù)ψ滿足如下非線性薛定諤方程

其中Lz＝r×p＝－i h－r×Δ＝－i h－（x?y－y?x）是角動量算符的z分量.波函數(shù)ψ滿足如下歸一化條件

N為總粒子數(shù).將波函數(shù)分解為ψ（x,y,t,z）＝?（x,y,t）φ（z）,令＝ψa0/,＝ω⊥t,＝r/a0,＝Ω/ω⊥,＝Lz/h－,在方程（3）中作變量代換并將波浪線去掉,則得到如下準二維無量綱化的非線性薛定諤方程：

其中g(shù)＝4πηNas,η＝∫d z｜φ（z）｜4∫dz｜φ（z）｜2,波函數(shù)?滿足歸一化條件

無量綱化之后的外勢阱的示意圖及沿x軸的截面如圖1所示.

圖1 環(huán)形勢示意圖及沿x軸的截面圖Fig.1 Schematic of annular potential and cross section of potential at y＝0

相應(yīng)地,與上述非線性薛定諤方程對應(yīng)的拉格朗日函數(shù)為

其中,

化學式μ為拉格朗日乘子.

2 托馬斯－費米近似

其中相位梯度 Δθ＝v＝Ω×r表示超流速.在托馬斯－費米近似中,通常忽略曲率密度上述自由能函數(shù)對求變分可得如下托馬斯－費米密度

2.1 出現(xiàn)中心洞的臨界條件

由（10）式中托馬斯－費米密度為零可得

或者

2.2 有中心洞的渦旋晶格的性質(zhì)

根據(jù)歸一化條件（6）可得

下面我們分三種情況分別討論旋轉(zhuǎn)環(huán)形勢阱中玻色愛因斯坦凝聚體的基態(tài)密度分布隨旋轉(zhuǎn)角頻率及勢阱結(jié)構(gòu)的變化情況.

情況1 勢阱寬度及相應(yīng)的中心高度一定,增大其旋轉(zhuǎn)角頻率

結(jié)合（12）,（15）,（16）式,可得

當Ω→Ωh時,d＝R2－R1≈R2,R≈d/2,隨著Ω的增大,圓環(huán)平均半徑R增大,而寬度d在減小.

情況2 勢阱旋轉(zhuǎn)角頻率一定,增大其寬度及相應(yīng)的中心高度

結(jié)合（13）,（15）,（16）式,可得

只有r0＞時,才能滿足μ＜/2.由式（14）和式（20）得

當r0→時,d＝R2－R1≈R2,R≈d/2,隨著r0的增大,R增大,d減小.

情況3 令勢阱旋轉(zhuǎn)角頻率為零,增大其寬度及相應(yīng)的中心高度

若旋轉(zhuǎn)角頻率為零,則由公式（13）可知出現(xiàn)中心洞的臨界條件為

此時

當r0→時,d＝R2－R1≈,R＝d/2,隨著r0的增大,R增大,d減小.

上述三種情況下無量綱化的環(huán)形勢阱V（r）＝V0（r2－r0）2/2束縛的玻色愛因斯坦凝聚體的托馬斯費米近似半徑在任何參數(shù)條件下都是收斂的.而無量綱化的環(huán)形勢阱V（r）＝V0（r－r0）2/2束縛的玻色愛因斯坦凝聚體只有在Ω2＜V0的條件下才能滿足托馬斯費米近似半徑收斂[18].

3 數(shù)值模擬

運用虛時演化數(shù)值方法模擬無量綱化的兩維非線性薛定諤方程（5）描述的凝聚體的基態(tài)密度分布并同解析結(jié)果進行比較[19－20].首先,我們選取g＝100,V0＝0.5,r0＝3,對于不同的旋轉(zhuǎn)角頻率Ω,凝聚體的基態(tài)密度分布如圖2所示.當Ω＝1.4時,凝聚體基態(tài)密度分布由圓盤變?yōu)閳A環(huán),隨著旋轉(zhuǎn)角頻率的增大,圓環(huán)平均半徑增大,寬度變窄,量子化渦旋數(shù)目增多且最終全部包含在中心洞的內(nèi)部形成巨渦旋；其次,選取g＝100,V0＝0.5,Ω＝2,對于不同的r0,凝聚體的基態(tài)密度分布如圖3所示,當r0＝1.5時,凝聚體基態(tài)密度分布由圓盤變?yōu)閳A環(huán),增大r0,凝聚體基態(tài)密度分布變化情況同增大旋轉(zhuǎn)角頻率Ω相似；最后,考慮強相互作用下凝聚體的基態(tài)密度分布隨r0的變化情況,取g＝1 000,V0＝0.5,Ω＝2,對于不同r0凝聚體的基態(tài)密度分布如圖4所示,當r0＝7時,密度分布由圓盤變?yōu)閳A環(huán),增大r0,圓環(huán)平均半徑增大,寬度變窄,量子化渦旋數(shù)目增多且一部分包含在中心洞的內(nèi)部形成巨渦旋,一部分呈鏈狀均勻的排列在環(huán)的內(nèi)部.圖5顯示,在形成巨渦旋的三種不同條件下,凝聚體基態(tài)密度分布的解析結(jié)果同數(shù)值結(jié)果相互吻合.另外,我們考慮角頻率Ω＝0時凝聚體基態(tài)密度分布隨r0增大的變化情況,圖6顯示g＝100,V0＝0.5,r0＝7和g＝1 000, V0＝0.5,r0＝12兩種情況下,凝聚體的基態(tài)密度分布都沒有出現(xiàn)量子化渦旋,但都出現(xiàn)了中心洞,且均與解析結(jié)果相一致.

圖2 g＝100時不同Ω對應(yīng)凝聚體的基態(tài)密度分布和相位（第一行為密度分布,第二行為相位）Fig.2 Ground state density profiles and phases for different angular velocities with g＝100（The first row is ground state density profiles.The second row is corresponding phase plots.）

圖3 g＝100時,不同r0對應(yīng)凝聚體的基態(tài)密度分布和相位（第一行為密度分布,第二行為相位）Fig.3 Ground state density profiles and phases for different r0with g＝100（The first row is the ground state density profiles.The second row is corresponding phase plots.）

圖4 g＝1 000時,不同r0對應(yīng)凝聚體的基態(tài)密度分布和相位（第一行為密度分布,第二行為相位）Fig.4 Ground state density profiles and phases for different r0with g＝1 000（The first row is ground state density profiles.The second row is corresponding phase plots.）

4 結(jié)論

通過托馬斯－費米近似和虛時演化方法求解可知,當粒子間相互作用強度和環(huán)形勢阱結(jié)構(gòu)一定時,增大勢阱旋轉(zhuǎn)角頻率到某一臨界值Ωh或當粒子間相互作用強度和勢阱旋轉(zhuǎn)角頻率一定時,增大控制環(huán)形勢阱寬度及相應(yīng)中心高度變化的參量r0到某一臨界值,凝聚體基態(tài)密度分布均從圓盤變?yōu)閳A環(huán).繼續(xù)增大Ω或者r0,圓環(huán)平均半徑增大,而寬度變窄,量子化渦旋數(shù)目增多.當粒子間相互作用較弱時,量子化渦旋會隨著

圖5 三種不同條件下,凝聚體基態(tài)密度分布沿x軸的截面（a）g＝100,V0＝0.5,Ω＝2.5,r0＝3,（b）g＝100,V0＝0.5,Ω＝2,r0＝4.5,（c）g＝1 000,V0＝0.5,Ω＝2,r0＝10Fig.5 Cross sections of ground state density profiles of condensates along x axis, （a）g＝100,V0＝0.5,Ω＝2.5,r0＝3,（b）g＝100,V0＝0.5,Ω＝2,r0＝4.5,（c）g＝1 000,V0＝0.5,Ω＝2,r0＝10

圖6 Ω＝0時凝聚體的基態(tài)密度分布（a）,（d）、相位（b）,（e）及密度沿x軸截面的分布（c）,（f）：（a）,（b）,（c）g＝100,V0＝0.5,r0＝7；（d）,（e）,（f）g＝1 000,V0＝0.5,r0＝12Fig.6 Ground state density profiles（a）,（d）,phase plots（b）,（e）and cross sections of density profiles of condensates along x axis（c）,（f）withΩ＝0：（a）,（b）,（c）g＝100,V0＝0.5,r0＝7；（d）,（e）,（f）g＝1 000,V0＝0.5,r0＝12

Ω或者r0的增大全部包含在中心洞中形成巨渦旋,當粒子間相互作用較強時,量子化渦旋隨著Ω或者r0的增大一部分包含在中心洞中形成巨渦旋,一部分呈鏈狀均勻排列在環(huán)的內(nèi)部,這里我們只給出了其隨r0增大的變化情況.當旋轉(zhuǎn)角頻率為零時,增大r0到某一臨界值,凝聚體基態(tài)密度分布也會從圓盤變?yōu)閳A環(huán),但不會出現(xiàn)量子化渦旋.凝聚體基態(tài)密度分布的解析結(jié)果同數(shù)值結(jié)果相互吻合.可見,旋轉(zhuǎn)是產(chǎn)生量子化渦旋的直接因素.當旋轉(zhuǎn)角頻率為零時,沒有量子化渦旋；當旋轉(zhuǎn)角頻率不為零時,增大旋轉(zhuǎn)角頻率或者增大勢阱寬度及相應(yīng)的中心高度,量子化渦旋的數(shù)目都會增多,且最終形成巨渦旋.

[1] Bretin V,Stock S,Seurin Y,Dalibard J.Fast rotation of a Bose-Einstein condensate[J].Phys Rev Lett,2004,92：050403.

[2] Ryu C,Andersen M F,Clade P,Natarajan V,Helmerson K,PhillipsW D.Observation of persistent flow of a Bose-Einsteincondensate in a toroidal trap[J].Phys Rev Lett,2004,99：260401.

[3] Haljan PC,Coddington I,Engels P,Cornell E A.Driving Bose-Einstein-condensate vorticity with a rotating normal cloud [J].Phys Rev Lett,2001,87：210403.

[4] Stringari S.Phase diagram of quantized vortices in a trapped Bose-Einstein condensed gas[J].Phys Rev Lett,1999,82：4371.

[5] Fetter A L.Rotating vortex lattice in a Bose-Einstein condensate trapped in combined quadratic and quartic radial potentials [J].Phys Rev A,2001,64：063608.

[6] Lundh E.Multiply quantized vortices in trapped Bose-Einstein condensates[J].Phys Rev A,2002,65：043604.

[7] Kasamatsu K,Tsubota M,Ueda M.Giant hole and circular superflow in a fast rotating Bose-Einstein condensate[J].Phys Rev A,2002,66：053606.

[8] Jackson A D,Kavoulakis G M,Lundh E.Phase diagram of a rotating Bose-Einstein condensate with anharmonic confinement [J].Phys Rev A,2004,69：053619.

[9] Fetter A L,Jackson B,StringariS.Rapid rotation ofa Bose-Einstein condensate in a harmonic plus quartic trap[J].Phys Rev A,2005,71：013605.

[10] Kim JK,Fetter A L.Dynamics of rapidly rotating Bose-Einstein condensates in a harmonic plus quartic trap[J].Phys Rev A, 2005,72：023619.

[11] Fu H,Zaremba E.Transition to the giant vortex state in a harmonic-plus-quartic trap[J].Phys Rev A,2006,73：013614.

[12] Aftalion A,Mason P.Rotation ofa Bose-Einstein condensate held under a toroidal trap[J].Phys Rev A,2010,81：023607.

[13] Jeng B W,Wang Y S,Chien C S.A two-parameter continuation algorithm for vortex pinning in rotating Bose-Einstein condensates[J].Comput Phys Comm,2013,184：493－508.

[14] Song T L,Ma Y L.Finite vortex numbers and symmetric vortex structures in a rotating trapped Fermi gas in the BCS-BEC crossover[J].Eur Phys JB,2011,82：303－311.

[15] Song T L,Ma C R,Ma Y L.Vortex lattices in a rotating Fermisuperfluid in the BCS-BEC crossover withmany Landau levels [J].Ann Phys,2012,327：1933－1947.

[16] HuaW,Liu X.Improved shootingmethod for Gross-Pitaevskii equation：Interference of three Bose-Einstein condensates[J].Chinese JComput Phys,2011,28（6）：922－926.

[17] Che L,Luo X,Su M.Interference of three Bose-Einstain condensates in one-dimensional time-dependent Gross-Pitaevskii equation[J].Chinese JComput Phys,2011,28（3）：463－468.

[18] Jin J J,Zhang S Y,Han W,Wei Z F.The ground states and spin textures of rotating two-component Bose-Einstein condensates in an annular trap[J].JPhys B：AtMol Opt Phys,2013,46：075302

[19] Antoine X,Bao W,Besse C.Computationalmethods for the dynamics of the nonlinear Schr?dinger equations[J].Comput Phys Comm,2013,184：2621－2633.

[20] Dalfovo F,Stringari S.Bosons in anisotropic traps：Ground state and vortices[J].Phys Rev A,1996,53：2477.

Ground States of Rotating Bose-Einstein Condensates in an Annular Trap

LIU Yan, ZHANG Suying
（Institute of Theoretical Physics,Shanxi University,Taiyuan 030006,China）

Thomas-Fermi approximation（TFA）and imaginary-time propagation method are used to study ground states of rotating Bose-Einstein condensates in an annular trap.Ground state density profiles of condensates experience a transition from vortex lattice phase to giant vortex phasewith increase ofangular frequency orwith increase ofwidth and center heightof trap potential.Particularly, ground state density profiles change from a disc shape into an annulus shape with the increase of width and center height of trap potential,when angular frequency is zero.Finally,comparison between ground state density profiles,obtained by analyticalmethod and numericalmethod ismade.They coincide with each other.

Thomas-Fermi approximation；central hole；giant vortex

1001-246X（2015）06-0744-07

O469

A

2014－12－05；

2015－02－13

國家自然科學基金（91430109）,高等學校博士學科點專項科研基金（20111401110004）及山西省自然科學基金（2014011005-3）資助項目

劉燕（1988－）,女,碩士生；研究方向：計算物理,E-mail：zhangsy＠sxu.edu.cn