磁阱中超冷玻色氣體的 BEC 相變與 Kibble - Zurek 動力學（下）

壹. 古典場模型與數值模擬

GP 方程式是平均場近似下的理論，它僅能描述凝聚體波函數在絕對零度下的同調演化。我們採用古典場方法來處理 BEC 在非零溫環境下的演化，包含了玻色氣體中同調與非同調的交互作用過程。古典場方法的原理簡述如下^[3]。

根據給定的能量閥值E_cut，如圖三（a）所示，系統可分成兩部分；（一）古典場區（C-field region）： 由所有能量低於E_cut 的模態組成，包含了 BEC 及其他高佔有數的低能量模態。（二）非同調區（incoherent region）：由能量高於 E_cut 之低佔有數激發態所組成（可視為一處於熱平衡溫度為 T 的熱庫）。系統的演化主要發生在古典場區，與 GP 方程有類似的形式。這也是此方法最大的優點，不僅在數值計算上容易處理並且可以適用在大部分的實驗參數範圍。

如圖三（b）、（c）所示，非同調區與古典場區粒子間的碰撞模式可分為：（i）非同調區中兩粒子碰撞後，其中之一進入古典場區；（ii）一古典場與一非同調區粒子碰撞，古典場區粒子數不變但佔有數重新分配。由於效應（ii）在古典場佔有數極大時可忽略，故我們僅考慮（i），此時古典場可以由一隨機投影GP 方程（stochastic projected GP equation，SPGPE）描述^[3,4]，

\[d\psi=\hat{\mathcal{P}}\left\{\frac{i+\gamma}{\hbar}\left[\mu-\hat{\mathcal{L}}\right]\psi\right\}dt+dW，
\hspace{1cm}(1)\]

其中投影算符 \( \hat{\mathcal{P}}\) 排除非古典場模態進行演化，而

\[\hat{\mathcal{L}}=-\frac{\hbar\nabla^2}{2M}+V_{\mathrm{trap}}\left(\mathbf{r}\right)+g\left|\psi\right|^2。 \hspace{1cm}(2)\]

\(dW\left(\mathbf{r},t\right)\)為經散射進入古典區的粒子所帶進之複數高斯白噪音，它滿足以下相干關係

\[\left\langle d W^\ast\left(\mathbf{r},t\right)dW\left(\mathbf{r}\prime,t\right)\right\rangle=\frac{2\gamma k_BT}{\hbar}\delta\left(\mathbf{r}-\mathbf{r}^\prime\right)dt，  \hspace{1cm}(3)\]

其中 \(\gamma \) 為古典區粒子的增加率，\(\mu\) 為化學勢，\(k_B\) 為波茲曼常數。依據“ 各態歷經假說”（ergodic assumption），可利用短時間平均來計算密度矩陣\(\rho\left(\mathbf{r},\mathbf{r}^\prime,t\right)=\left\langle\psi^\ast\left(\mathbf{r}^\prime,t\right)\psi\left(\mathbf{r},t\right)\right\rangle \)。再依據 Penrose － Onsager 定理， 亦即 BEC 模態（或是 PO 模態）對應到密度矩陣之最大本徵值的解，我們便可以計算 BEC 在隨機噪訊環境下的演化。

根據 Trento 團隊的實驗參數^[1,2]，我們考慮於 \( -\tau_Q < t < \tau_Q\) 期間化學勢與溫度呈線性變化：

取\(\mu_f=\mathrm{22}\hbar \omega_{⊥}，T_0=500 \hspace{1mm}\rm nK\) 與\(\Delta T = 290 \hspace{1mm}\rm nK\) ；位能勢\(V_{\mathrm{trap}}\left(\mathbf{r}\right)=\left({\mathrm{M}}/{\mathrm{2}}\right)\left[\omega_x^2x^2+\omega_\bot^2\left(y^2+z^2\right)\right] \)的振盪頻率，\(\omega_x=\mathrm{2}\pi×13 \) Hz，\(\omega_\bot=\mathrm{2}\pi×131.4 \) Hz。被捕捉在此位能阱的原子團呈現狹長（\(\omega_\bot\approx10\omega_x \)）的雪茄型，屬三維系統，故可觀測到的缺陷為渦漩線。我們透過 SPGPE 模擬來探索該系統相變點之動力學行為^[4,5]。系統的初態設定為一熱平衡的玻色氣體，末態為一內含 BEC 的超冷玻色氣（凝聚體的粒子數約佔原有氣體粒子數的75%）。圖四中，我們考慮\(\tau_Q= 9 \sim 720 \hspace{1mm}\rm ms\)，\(\gamma=0.005\) 。根據不同的 \( \tau Q\)，我們歸納出以下三種降溫模式：（a）緩慢降溫—系統在通過相變點後，逐一生成分散的 BEC，但無大量的缺陷生成。因為熱擾動影響，缺陷將會離開系統，故系統有很大的機會觀測不到缺陷；（b）典型降溫—系統生成數個拓樸缺陷，隨後因為熱漲落與缺陷間的交互作用，僅有少數能夠存活長時間。（c）快速降溫—系統出現大量的缺陷，即便缺陷的平均生命期有限，但依然能夠有數個能夠存活的非常久。

貳. 線性近似與臨界暴漲

系統粒子數由化學勢 \(\mu\) 決定。對於均勻系統，當 \(\mu > 0\) 時，系統將發生凝聚態現象，故可以藉由μ 來做為相變調控參數。當系統尚未發生凝聚態現象時，即 \(\mu > 0\)  ( \(t < 0\) ) ，古典場波函數僅由無序的相位構成且 \( \psi \sim 0\)（及其期望值），故我們能夠忽略式（2）中的非線性項。對SPGPE 式（1）做傅立葉變換得到\[\frac{d\widetilde{\psi}}{dt}=\left[\frac{\left(\gamma+i\ \right)}{\hbar}\frac{\hbar^2k^2}{2M}+\frac{\gamma}{\hbar}\mu\left(t\right)\right]\widetilde{\psi}+\widetilde{\zeta}，\hspace{1cm} (5)\] 其中 \(\widetilde{\psi} \) 為古典場波函數的傅立葉轉換，而噪訊項 \(\widetilde{\xi}\) 則遵守相干條件，

\[\left\langle{\widetilde{\zeta}}^\ast\left(\mathbf{k},t\right)\widetilde{\zeta}\left(\mathbf{k}^\prime,t\right)\right\rangle\ =\frac{2\gamma k_BT}{\hbar}\delta\left(\mathbf{k}-\mathbf{k}^\prime\right)\delta\left(t-t\prime\right)。
\hspace{1cm}(6)\]

隨機方程式（5）於動量空間中有精確解且僅躁訊項有主要貢獻，其精確解形式為
\[\widetilde{\psi}\left(\mathbf{k},t\right)=\int_{-\infty}^{t}{dt^\prime\widetilde{\zeta}\left(\mathbf{k},t^\prime\right)e^{-\frac{\gamma\hbar k^2}{2M}\left(t-t^\prime\right)+\frac{\gamma\mu_f}{2\hbar\tau_Q}\left(t^2-t^{\prime2}\right)}} \hspace{1cm}(7)\] 

古典場波函數之動量分布（譜函數）為可量測之物理量，其定義為
 \[f\left(\mathbf{k},t\right)=\left\langle{\widetilde{\psi}}^\ast\left(\mathbf{k},t\right)\widetilde{\psi}\left(\mathbf{k},t\right)\right\rangle。 \hspace{1cm}(8)\] 



我們可以假定當 \(t > 0\)，系統發生BEC 相變，故上述的解僅適用於接近相變點。利用高斯積分將其近似為
\[f\left(\mathbf{k},t\right)\approx\frac{2\sqrt\pi\gamma k_BT}{\hbar}\hat{t}e^{\left({t}/{\hat{t}}-\hat{\xi}k^2\right)^2}， \hspace{1cm}(9)\] 
其中
 \[\hat{t}=\left(\frac{\hbar\tau_Q}{\gamma\mu_f}\right)^{-\frac{1}{2}}， \hat{ξ}=τ_Q^{1/4}  \left(\frac{4M^2μ_f}{γℏ^3} \right)^{-\frac{1}{4}}。 \hspace{1cm}(10)\] 

由式（9），我們發現譜函數滿足標度不變性（scaling invariant）
 \[f\left(\mathbf{k},t\right)=\hat{t}F\left(\frac{t}{\hat{t}},\hat{\xi}\mathbf{k}\right)， \hspace{1cm}(11)\] 
另外，\({\hat{t}}/{{\hat{\xi}}^3}\propto{\hat{\xi}}^{-1} \)符合式（5）的標度行為假設。

上述標度行為適用於\({\ \tau}_Q\gg\hat{t}\  \)之模擬，因其不受式（4）之降溫參數停止變化之影響。此外要明確地知道相變點的位置，在考慮位能阱及粒子交互作用的情況下，我們知道相變點不會發生在 \(\mu = 0\)。我們透過研究在不同時刻下給定 \(\mu\) 與 \(T\) 之古典場平衡態性質發現相變發生在 \(t_c ≈ (0.188 ± 0.013)τ_Q \) ^[4]，於此我們將對於數值計算中之 \(t\) 平移至 \(t_c\)。在引入相變點位移後，我們發現模擬數值計算滿足式（11）。圖五中以 \(k = 0\) 為例，將譜函數經根據式（11）進行標度轉換後，譜函數至\(\left(t-t_c\right)/\hat{t}\approx\mathrm{1.25}\  \)左右皆遵守標度律。

由式（9），譜函數在通過相變點後會隨時間增加呈指數發散，其特徵時間為 \( \hat{t}\)，但不同波數的發散時刻會隨著 \(k\) 增加而延後。由譜函數我們回推對應的系統密度為
\[\left\langle\left|\psi\left(\mathbf{r},t\right)\right|^2\right\rangle\approx\frac{\gamma k_BT}{16\pi^{3/2}\hbar}\frac{\hat{t}}{{\hat{\xi}}^3}\left(\frac{\hat{t}}{t}\right)^\frac{3}{2}{e^{\left(t/\hat{t}\right)}}^2 \hspace{2cm}，(12)\]

本近似預期可精確至 \(t\sim\hat{t}\)，因隨著時間演化密度成指數發散，非線性項將不可忽略，系統也將不再臨界發散，故可以被激發的最高動量波數\(~\ {\hat{\xi}}^{-1}\)將決定系統的最小尺度，即相干長度之尺度。我們發現與KZ 圖像不同是，缺陷之相干長度並非由相變點前之凍結時域決定，而是由過相變點後，系統所能激發的最大波數決定。故此處的 \( \hat{t}\) 與KZ圖像之凍結時域大小相符但僅發生在相變後，且恰巧對應到 \(zν = 1\) 其中 \(z = 2\) 之平均場預測。

參. 降溫相圖

以線性近似分析SPGPE，我們可以更進一步利用局域密度近似來分析Trento 團隊的實驗。為利於討論與計算，我們依然假定相變點發生在 \(\mu = 0 \)。在陷阱中，空間中各區間發生相變的條件為
\[\mu\left(t\right)-V_{\mathrm{trap}}\left(\mathbf{r}\right)>0，\hspace{2cm} (13)\]
可推知系統中相變波前尺寸為
\[a_{x/\bot}\left(t\right)=\sqrt{\frac{2\mu_f}{M\omega_{x/\bot}^2}\frac{t}{\tau_Q}}，\hspace{2cm} (14)\]
與其對應的波前速率
\[v_{x/\bot}=\frac{da_{x/\bot}}{dt}=\sqrt{\frac{\mu_f}{2M\omega_{x/\bot}^2\tau_Qt}}，\hspace{2cm} (15)\]
比較系統發生相變後，聲速地平線波前（相干長度增長）增長的速率
\[\hat{v}\approx\frac{\hat{\xi}}{\hat{t}}=\left(\frac{\gamma^3\hbar\mu_f}{4M^2}\right)^\frac{1}{4}\tau_Q^{-\frac{1}{4}}，\hspace{2cm} (16)\]
當聲速地平線擴展速度小於相變波前時，我們可得勻相軸長為
\[{\hat{a}}_{x/\bot}=\left(\frac{4\mu_f^3}{M^2\omega_{x / \bot}^8\gamma\hbar\tau_Q^3}\right)^{\frac{1}{4}\ }\hspace{2cm} (17)\]
反之當聲速地平線擴展速度大於相變波前時，即當 \( a_x (t)/ \hat{a}_x\)或 \(a_{\perp}(t)/\hat{a}_{\perp}\) 大於1 時，勻相近似將會失準。我們整理出以下幾個降溫模式

\( \tau_Q\lesssim\hat{t} \)，          準瞬間降溫
 \( \hat{t}<\tau_Q<\frac{\mu_f}{\gamma\hbar\omega_\bot^2} \)，    勻向降溫
\(\frac{\mu_f}{\gamma\hbar\omega_\bot^2}<\tau_Q<\frac{\mu_f}{\gamma\hbar\omega_x^2}\)，長軸勻向降溫
\(\frac{\mu_f}{\gamma\hbar\omega_x^2}<\tau_Q\)。  非勻向降溫

圖六為上述不同降溫模式所對應的相圖。我們的降溫範圍剛好處於準瞬間降溫與勻向降溫區間，其中前節所述之” 快速降溫”對應到準瞬間降溫，而典型與緩慢降溫屬於勻向降溫，可以解釋即便我們在線性近似之預測所中忽略的外加陷阱，卻可以依然適用於我們的數值計算中。想要探討長軸勻向與非勻向降溫動力學，根據本參數需要極緩慢的降溫，可預期會有更加複雜有趣的結果。


肆. 後相變時期

最後，我們比較Trento 實驗進行KZ 缺陷數目量測的結果[1]。該實驗使用飛行時間（Time of flight）量測法， 觀測在特定時刻凝聚體由位能阱釋放後膨脹的影像。他們僅能在凝聚體生成後約250 ms 左右進行有效的缺陷數測量。由上述理論與模擬，在臨界暴漲期間，即便已過了相變點，PO 模態尚未能夠顯示完整的凝聚體[4]，故可以預期實驗上可觀測到凝聚體的時間是已過了臨界暴漲期。為了便於與實驗觀測比較，我們模擬實驗的量測方式，以PO 模態粒子數達最終平衡態粒子數的5% 時為實驗可能量測到凝聚體生成的時間點，標記為 \(t_bec\)。在圖七我們比較 \(t − t_bec\) 在50 ms 與200 ms 所能觀測到缺陷數目與實驗量測之結果。因實驗僅能控制溫度隨時間之變化，故我們將可觀測之渦旋數 \(N_{\mu}\) 對 \(dT/dt \) 作圖比較，發現數值計算結果整體行為與實驗吻合，實驗中所量測到的標度律為 \(N_{\mu}\propto \tau Q^{−7/3}\)，符合受侷限系統之預測[6]。即便因為熱漲落，數值模擬在後期所能觀測到的粒子數較實驗結果為少，但整體趨勢一致。

快速降溫出現極值之平台，可以理解為因降溫過程生成大量之缺陷，在SPGPE 模擬中我們發現因受交互作用影響、熱漲落而將缺陷排出於系統外，導致其數目快速地減少，故僅存數個可觀測的缺陷。這些缺陷因距離遙遠，猶如相互隔離的缺陷，即便在熱漲落下，依舊可以存活相當長的時間。而較緩的降溫則因本身能夠生成的缺陷數變不多，且如上述相互隔離的缺陷可以相當地長壽，故能夠吻合KZ 圖像預測。總而言之，想要簡單地驗證KZ 機制，需要具有穩定不受熱漲落影響的缺陷系統[7]。

然而即便這些隔離的長壽缺陷，最終還是會受到熱漲落而離開系統。根據模擬與實驗顯示，單一的渦漩最長可存活超過 1 秒，而接著系統將演化至一有限溫度、處於平衡態的凝聚體，如圖三（a）所示。

伍. 總結

本文利用古典場近似探索一受侷限之超冷玻色氣體經過 BEC 相變點之動力學行為。我們發現缺陷數並非如 KZ 圖像所描述，由被凍結的相位決定，而是在通過相變點後，系統將發生臨界暴漲，此時所能激發之最大動量模態，或是最小系統尺度決定。其伴隨著密度暴增，將離開暴漲時域，而此時域恰與 KZ 圖像之凍結時域相符。透過在 SPGPE 線性分析中引入局域密度近似，得出勻向－非勻向降溫相圖，而本文探索的參數剛好處於勻向區。

透過 SPGPE 數值模擬（這是首宗引入化學勢與溫度參量的數值模擬）與分析古典場中的 PO 凝聚體模態，我們得以更進一步探究凝聚態在相變後之動力學行為，得知了生成的拓樸缺陷會因為彼此的交互作用而排出系統，導致有大量的自發缺陷無法順利在實驗中被觀測到，數值模擬結果與實驗有相當好地吻合，並了解在系統如何逐漸演化至平衡態。

本研究由英國 Newcastle 大學、義大利 Trento 大學、波蘭 Jagiellonian 大學與台灣國立彰化師範大學四方合作而成，論文詳見[4,5]。

參考文獻：

[1] S. Donadello, S. Serafini, T. Bienaimé, F. Dalfovo, G. Lamporesi and G. Ferrari, “Creation and counting of defects in a temperature-quenched Bose-Einstein conden-sate”, Phys. Rev. A 94, 023528（2016）.

[2] S. Donadello, S. Serafini, M. Tylutki, L. P. Pitaevskii, F. Dalfovo, G. Lamporesi and G. Ferrari, “Observation of Solitonic Vortices in Bose-Einstein Condensates” Phys. Rev. Lett. 113, 065302（2014）.

[3] 蘇士煒、劉翼綱、郭西川， “ 冷原子不冷－非零溫玻色－愛因斯坦凝聚體的動力學，” 物理雙月刊，第34 卷

3 期，184（2012）。

[4] I.-K. Liu, S. Donadello, G. Lamporesi, G. Ferrari, S.-C. Gou, F. Dalfovo and N. P. Proukakis, “Dynamical equilibration across a quenched phase transition in a trapped quantum gas”, Comm. Phys. 1, 24（2018）.

[5] I.-K. Liu, J. Dziamaga, S.-C. Gou, F. Dalfovo and N. P. Proukakis, “Kibble-Zurek dynamics in a trapped ultracold Bose gas” Phys. Rev. Research 2, 033183（2020）.

[6] A. del Campo, A. Retzker, and M. B. Plenio, “The inhomogeneous Kibble-Zurek mechanism: Vortex nucleation during Bose-Einstein condensation”, New Journal of Physics 13, 083022（2010）.

[7] S.-W. Su, S.-C. Gou, A. Bradley, O. Fialko, J. Brand, “Kibble-Zurek scaling and its breakdown for spontaneous generation of Josephson vortices in Bose-Einstein condensates”, Phys. Rev. Lett. 110, 215302（2013）.