二維激子極化激元凝聚中渦旋疊加態穩態及動力學特性研究

陳海軍, 任 元, 王 華, 湯國志, 劉 通

(中國人民解放軍戰略支援部隊航天工程大學, 北京 101416)

1 引 言

實驗上實現了冷原子玻色愛因斯坦凝聚(BEC)以后, 與之相關的理論研究和實驗研究層出不窮[1-4]. 實現冷原子BEC需要相當低的溫度, 通常通過磁光阱和蒸發冷卻技術來完成[5]. 低溫下實現冷原子BEC, 給BEC的工程應用帶來了限制. 近年來人們發現半導體微腔中的激子極化激元凝聚系統(exciton-polariton condensates)有望在室溫下實現BEC[6], 給BEC的應用從實驗室走向工程實踐提供了可能. 激子極化激元BEC和冷原子BEC具有許多相同的性質, 例如在其中可以形成穩定的孤立子, 約瑟夫森振蕩, 宏觀自捕獲和空間干涉等[7].

實現BEC以后, BEC超流特性的研究引起了足夠重視, BEC中量子化渦旋的出現是BEC具有超流特性的直接證據[8]. 量子化渦旋的存在形式有單個渦旋, 渦旋陣列和正反渦旋疊加態.渦旋疊加態是由具有相同軌道角動量量子數的正反渦旋疊加形成渦旋疊加態. 2006年Liu M等人研究了單組分冷原子BEC中正反渦旋疊加態的基本性質, 并提出了一種產生正反渦旋疊加態的有效方案[9]. 2012年Thanvanthri S等人研究了基于原子BEC的正反渦旋疊加態在超穩物質波陀螺中的應用, 并首次提出了渦旋干涉陀螺的概念[10]. 2015年Padhi A N B等人研究了激子極化激元凝聚體系中單渦旋在旋轉參考系中形成渦旋陣列中渦旋個數隨系統參數的變化規律[11]. 2016年Dai W D等人進一步研究了扁平勢和無序勢中渦旋疊加態的穩態結構及其Sagnac效應, 理論闡明了渦旋干涉陀螺中Sagnac干涉相位, 軌道角動量量子數和半導體微腔旋轉角速率三者之間的關系[12].

本文利用分步Crank-Nicolson方案的虛時和實時有限差分方法[13]求解耗散系統的Gross-Pitaevskii(GP)方程, 找出了諧振子勢阱和高斯型勢壘相結合的復合阱中的激子極化激元凝聚中正反渦旋疊加態的穩態結構, 并直觀地驗證這種穩態結構在半導體微腔旋轉下的穩定性, 給出了穩定的渦旋疊加態和半導體微腔的旋轉角速率之間的關系, 確定了渦旋疊加態在旋轉參考系中具有Sagnac效應. 另外, 研究了泵浦光寬度和增益對系統中形成渦旋陣列動力學過程的影響.

本文的結構安排如下: 第二部分給出了描述耗散系統的無量綱化形式的GP方程, 系統的能量, 化學勢, 角動量和自由能和數值計算方法. 第三部分討論了正反渦旋疊加態的穩態結構. 第四部分給出了半導體微腔旋轉下正反渦旋疊加態的穩定性及疊加態旋轉角速率和半導體微腔旋轉角速率之間的定量關系. 第五部分討論了泵浦光寬度和增益對渦旋陣列形成動力學的影響.

2 理論模型和數值計算方法

2.1 GP方程和系統參數

(1)

(2)

聯合方程(1)和(2)得到

(3)

在零溫極限下, 采取Bogoliubov近似, 把場算符用波函數替代, 得到描述玻色體系動力學行為的Gross-Pitaevskii方程.

(4)

其中g=4π?2a/m表示粒子之間的兩體相互作用強度.

對于弱相互作用的激子極化激元體系仍然可以用方程(4)描述其動力學過程. 由于激子極化激元系統中存在泵浦和損耗, 是非平衡系統, 因此要在方程(4)的基礎上加入耗散項和增益項.另外, 為了研究半導體微腔旋轉時體系的動力學特性, 還要引入旋轉項. 方程(4)是三維方程, 我們考慮的是準二維體系, 因此要在凍結z方向自由度的基礎上對方程進行約化處理得到二維方程. 另外, 為了數值計算方便, 我們選取體系的特征長度, 時間和能量

t0=1/ω～ps

E0=?ω～0.66 meV

(5)

對方程進行無量綱化處理. 最終得到的無量綱化GP方程為

(6)

為了研究系統的穩態結構和渦旋演化動力學過程, 數值計算過程中可以根據系統參數的變化趨勢來判斷系統是否達到穩態, 包括系統能量, 耗散能和化學勢.在單渦旋向渦旋陣列的演化過程中, 自由能的變化尤為重要, 它們的表達式分別是

▽Ψ|2+

(7)

2.2 數值計算方法

由于激子極化激元凝聚系統存在粒子數的增加和損耗, 而虛時演化方法每一步數值計算對波函數進行了歸一化, 相當于忽略了粒子數的損耗和增加, 因此不能直接利用虛時演化方法求解激子極化激元系統的穩態解. 我們采用的方法是先忽略方程(6)中的虛數項, 利用虛時演化方法找出不存在泵浦, 損耗和增益的情況下, 系統的穩態解, 然后利用實時演化方法, 驗證加上泵浦, 損耗和增益時, 穩態解能夠穩定存在所對應的泵浦, 損耗和增益值, 這樣就找出了特定參數下激子極化激元系統中穩定的渦旋疊加態.

諧振子系統的基態解形式是高斯型, 因此數值計算過程中, 我們選用如下形式的高斯型渦旋疊加態函數作為初始波函數.

[(x+iy)l1+(x-iy)l2]

(8)

w是初始波包的寬度,l1和l2是渦旋量子數, 根據其不同取值可以構造出單渦旋態和正反渦旋疊加態, 尤其當l1=l2=l時, 所形成的正反渦旋疊加態呈現對稱的“花瓣”狀, 花瓣的個數為2l.

3 渦旋疊加態的穩態結構

圖1是虛時演化結果. 只存在諧振子勢阱時, 系統的穩態解應為高斯型函數, 當給定正反渦旋疊加態時, 最終也會演化為高斯型. 為了阻止這種演化, 形成穩定的正反渦旋疊加態, 在諧振勢阱中央加入一個高斯型勢壘阻止這種演化, 表達式為

圖1 虛時演化方法計算穩態結構的過程, l=2, g=1, wc=1.2, wb=1.2, b0=10.5Fig. 1 The process of calculating steady-state structure imaginary-time propagation method, l=2, g=1, wc=1.2, wb=1.2, b0=10.5

b0exp[-(x2+y2)/(2wb)]

b0和wb分別表示勢壘的高度和寬度, 這樣所形成的簡諧勢阱和勢壘疊加的囚禁結構類似于文獻[22]中提出的墨西哥帽子勢阱.t=0時,m=2, 花瓣的個數為4, (a)給出初始時刻體系的密度分布|Ψ(x,y)|2(紅的部分密度大, 藍的部分密度小, 下同), 渦旋量子數是2, 呈現花瓣數為4的對稱狀態. 為了數值計算的穩定性, 非線性相互作用項g隨時間逐漸增加. (b)給出非線性相互作用項達到預定值時系統的密度分布. (c)表示經過一定時間演化所形成的穩態結構, 判斷系統達到穩態的標準是系統參數(能量, 化學勢, 自由能)不隨時間進行變化, 保證在一定的誤差范圍內波動. (d)(e)(f)是和上述密度分布所對應的相位分布. 可以看出, 給定初態后, 經過很短時間的

演化系統就行形成了穩態結構, 穩態結構和初態結構不同, 這是由于諧振子勢阱和中心勢壘所形成的環狀囚禁區域, 使得粒子沿著圓環方向分布, 形成穩定的花瓣結構.

圖1是無虛數項時用虛時演化方法的到的穩態結果, 保持所有參數取值不變, 當加上虛數項時(泵浦、損耗和增益), 這時虛時演化結果一般不再是系統的穩態結構.為了使無虛數項時的穩態結構在加上虛數項時仍然是體系的穩態結構, 必須使虛數部分ρ-γ-η|Ψ|2總體效應為0, 由此結論可以得出|Ψ|2=(ρ-γ)/η, 同時|Ψ|2還受制于方程(6)中其它部分的制約.我們的計算方案是給定損耗和飽和增益項, 然后根據系統參數(7)隨時間演化趨勢調整泵浦功率, 當泵浦功率取某一合適值時, 泵浦, 損耗和增益項互相平衡, 系統參數維持穩定狀態不隨時間變化, 此時的泵浦功率為臨界泵浦功率.

在計算過程中, 粒子之間的兩體相互作用強度g=1, 損耗和增益飽和項取η=γ=0.001, 泵浦光寬度wp=6. 圖2第一行和第二行分別給出穩定正反渦旋疊加態結構的密度和相位分布.前三列表示無虛數項的情況; 后三列表示存在泵浦, 損耗和增益的情況, 對應的臨界泵浦功率分別是,

p1=0.00115442,p2=0.00222305,

p3=0.00273581,

(9)

第一, 四列渦旋量子數l=2, 勢壘寬度和高度分別是wb=1.2,b0=10.5; 第二, 四列m=3,wb=2.6,b0=400；第三, 六列m=4,wb=2.6,b0=1200.

圖2 l=2, 3, 4三種情況下有無GP方程中虛數項對應的穩態結構密度和相位分布Fig. 2 Steady-State structure density and phase distribution in GP equation without imaginary and imaginary items for l=2, 3, 4

可以看出隨著渦旋量子數增加, 離心勢能逐漸增加, 正反渦旋疊加態半徑和能量隨之增加;相同渦旋量子數條件下, 有無虛數項所得到穩態結果完全一致, 說明泵浦, 損耗和增益飽和三項只要達到平衡, 不會對穩定狀態的結構有所影響;渦旋量子數越大, 需要的泵浦光強度越大, 勢壘強度和寬度也隨之增加.

4 旋轉情況

4.1 渦旋疊加態的Sagnac效應驗證

為了研究激子極化激元凝聚體系中的超流特性, 渦旋動力學是其中重要的研究內容. 我們考慮了正反渦旋疊加態(l=2)的旋轉動力學特征. 由于正反渦旋疊加態的總角動量為0, 當體系旋轉時, 正反渦旋會產生附加的相位差, 導致密度分布所形成的花瓣整體旋轉.

圖3給出了l=2的正反渦旋疊加態在半導體微腔角速率Ω=0.1ω時隨時間演化的動力學過程, 第一, 二行分別表示不同時刻的密度分布和相位分布(a)-(f). 可以看出, 一方面, 當凝聚體隨著體系一起旋轉時, 初始給定的正反渦旋疊加態的密度和相位分布整體形狀不隨時間變化, 說明渦旋疊加態是長時間穩定的(計算時長為500); 另一方面, (g)圖表示密度分布沿x軸的截面圖隨時間變化情況, 可以看出截面圖隨時間分布周期性地交替變化, 說明疊加態的轉動是勻速的, 為了討論旋轉速度的具體大小, 我們考慮了1/4周期內的旋轉情況, 從(h)圖可以看出, 疊加態1/4周期大約是16 ps, 則時間周期是64, 根據T=2π/Ω=2π/0.1≈63, 和我們數值計算結果一致.

圖3 l=2的穩定渦旋疊加態隨時間的旋轉過程Fig. 3 The rotational process of a stable vortex superposition state over time for l=2

4.2 渦旋疊加態的旋轉角速率

為了進一步闡明渦旋疊加態的旋轉角速率, 半導體微腔角速率和渦旋軌道角動量量子數三者之間的關系, 我們進一步在不同軌道角動量量子數與半導體微腔旋轉速率的條件下進行對比計算, 結果如圖4所示. 從圖4(a)可知, 疊加態的轉動周期Ti隨著半導體微腔速度Ω的增加而減小, 且與渦旋軌道角動量l沒有關系. 圖4(b)進一步給出了疊加態旋轉速率Ωi和Ω之間的關系, 可以看出, 在不同l取值下, Ωi和Ω始終相等. 以上分析表明, 疊加態BEC和慣性空間之間是相對穩定的.

事實上, 根據圖3(g), 我們可以定義干涉圖樣變化頻率為fd=1/△Td, 其中△Td表示在微腔角速率為Ω的條件下實現某一時刻渦旋密度分布與靜態條件下的渦旋密度分布完全一致所用的最短時間. 圖5給出了不同軌道角動量量子數下花瓣重合速度(渦旋密度變化速率)Ωd與微腔角速率Ω之間的關系. 從圖5可以看出, 在給定l的條件下, 渦旋密度變化速率Ωd與微腔角速率Ω成正比例關系, 且其斜率滿足關系Ωd/Ω=2l, 即在軌道角動量量子數為2,3,4的情況下, 斜率分別為4,6,8. 因此密度分布的相位變化是

Δφd=ΩdTd=2/ΩTd

(10)

圖4 疊加態旋轉周期和速率與微腔旋轉角速率之間關系Fig. 4 Relationship between rotation period and rate of superposition state and the rotational angular rate of microcavity

文獻[11]穩定情況下的密度分布為

[1+cos(2|l|φ)]|Ψ(|l|,ρ,z)|2

(11)

旋轉相位Δφd后, 密度分布應為

[1+cos(2|l|φ+ΩdΔT))]|Ψ(|l|,ρ,z)|2

(12)

與文獻[11]中的結論一致.

圖5 密度變化頻率隨微腔角速率變化Fig. 5 Variation of density change frequency with the angular velocity of microcavity

5 單渦旋態的動力學演化過程

圖6 旋轉參考系中渦旋晶格隨泵浦光寬度的變化, l=1, g=10, p=2, η=γ=1, wc=1.2Fig. 6 The change of vortex array in rotating reference system with the width of pumped light, l=1, g=10, p=2, η=γ=1, wc=1.2

單渦旋態的動力學演化過程如圖6和7所示, 給出了復數項中高斯型泵浦光寬度wp和增益飽和項η對動力學演化過程的影響. 我們取Ω=0.95ω研究旋轉帶來的動力學過程, 根據超流的特性, 當體系的旋轉速率超過某一臨界值時, 體系中會出現三角形排列的渦旋晶格, 我們的計算結果也是如此.

圖7 旋轉參考系中渦旋晶格隨增益飽和項的變化, l=1, g=1, p=2, γ=1, wp=5.5, wc=1.2Fig. 7 The change of vortex array with gain saturation in rotating reference system, l=1, g=1, p=2, γ=1, wp=5.5, wc=1.2

圖6給出的是激光寬度對動力學過程的影響, 從上到下每一行對應的泵浦光寬度分別是wp=4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 第一列對應t=500時的粒子密度分布, 第二列是和第一列對應的相位分布. 主要結果有: (1) 由于激子極化激元凝聚是激子和光子的耦合疊加態, 渦旋區域隨著泵浦光寬度的增加而增加. (2) 渦旋是從凝聚體邊緣逐漸進入中心, 進入時間隨著激光寬度增加而提前,wp=4.5時在我們計算的時間范圍內沒有渦旋進入凝聚體,wp=7.0時在很短時間內就有渦旋進入凝聚體. (3) 渦旋陣列呈三角形排列, 渦旋個數隨泵浦光寬度增加而增加. (4) 初始狀態的中心渦旋隨著時間演化會突然消失, 進而以小渦旋形式代替.

圖7是動力學過程隨飽和增益項η的變化, 從上到下每行對應的飽和增益項分別是η=0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 每一列的的含義和圖6一致, 可以看出: (1) 渦旋從邊緣進入凝聚體的時間隨η的增加而提前, 并且渦旋個數隨η增加而減少. (2) 渦旋區域大小沒有明顯變化. (3) 中心渦旋隨著時間增加會分裂成小渦旋, 并偏離中央區域, 分裂時間沒有明顯規律, 圖6當wp=4.5時, 在計算時間內沒有渦旋進入凝聚體系, 中央渦旋也沒有產生分裂, 圖7中當η=1.3時, 雖然有渦旋較早進入凝聚體系, 但是在計算的時間內, 中央渦旋沒有產生分裂. 從動力學過程中可以看出, 雖然能量, 化學勢等物理參數不隨時間變化, 但是最終的渦旋晶格狀態是動態變化的.

6 結 論

通過虛時和實時演化方法求解耗散GP方程研究了諧振子勢阱和高斯型勢壘疊加形成的類墨西哥帽子勢阱中正反渦旋疊加態的穩態結構, 旋轉動力學過程以及泵浦光寬度和增益對單渦旋演化成渦旋陣列過程的影響. 由于耗散GP方程中包含泵浦, 損耗和增益飽和項, 不能直接利用虛時演化方法求解系統的穩態解, 可以利用虛實演化相結合的方法研究穩態結構和實時方法經過長時間演化研究系統的穩態結構和動力學過程.

在此復合勢阱中可以形成不同軌道角動量量子數的穩定疊加態, 且其可以在旋轉條件下保持穩定. 渦旋疊加態形成的干涉圖樣相對于慣性空間是穩定的, 且渦旋密度分布的變化速率正比于半導體微腔的旋轉速率和角動量量子數, 這與基于Sagnac效應的理論分析結果一致. 單渦旋態在高速旋轉情況下會形成渦旋晶格, 渦旋晶格是超流體的特征之一, 證明了激子極化激元凝聚具有超流性, 并且證明泵浦光寬度和增益飽和項對渦旋陣列結構有明顯影響.