關於量子蒙地卡羅退火法 - 政大學術集成
75
0
0
全文
(2) 立. 政 治 大. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. 2. i n U. v. DOI:10.6814/NCCU201901066.
(3) 致謝 能夠完成論文，首先最為感謝的是我的指導教授林瑜琤老師。無論 在學業還是在生活上，老師總是給予我相當大的幫助以及鼓勵。在程 式以及論文上，更是投入許多時間用心盡力地指導我，使我有能力完 成這篇論文。再來感謝張太乙學長用心維護所上的機器，使我能夠安 心且能全力地計算論文所需的數據。感謝行政人員筱嘉姐、楊于廷學 姐以及恆毅在這段期間幫我許多忙。也感謝所上的老師、同學以及學 弟妹的支持與鼓勵。最後感謝我的父母，在我求學這一路上以來，總 是支持並鼓勵我，讓我能無後顧之憂地完成學業。. 立. 政 治 大. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. i. i n U. v. DOI:10.6814/NCCU201901066.
(4) 立. 政 治 大. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. ii. i n U. v. DOI:10.6814/NCCU201901066.
(5) 摘要 本論文檢驗以蒙地卡羅模擬退火來探討平衡態相變點定標分析之可 能性。以量子易辛模型為例，我們分別探討動力學指數為 z = 1 的量 子臨界點，具 z = ∞ 的無序量子臨界點，及 KosterlitzThouless (KT) 相變。應用有限溫度隨機級數展開法及基態投射演算法，我們考慮的 退火路徑涵蓋降溫、降橫場（量子擾動項）及同時降溫及降場三種情 形。我們的計算結果顯示對於 z = 1 量子臨界點，上述後兩類量子退 火過程在緩慢改變參數下均能正確反應臨界點位置及臨界指數。通過 KT 相變的退火過程亦可找出吻合理論的定標行為。唯 z = ∞ 的量子 臨界點為退火過程的瓶頸，似乎任意緩慢的退火速率均很難突破這個 瓶頸來達到無序系統近似靜態的極限。. 立. 政 治 大. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. v. 關鍵字：模擬退火、隨機級數展開量子蒙地卡羅演算法、零溫投射蒙 地卡羅演算法、非均質量子易辛鏈、三角反鐵磁、KosterlitzThouless 相變. Ch. engchi. iii. i n U. DOI:10.6814/NCCU201901066.
(6) 立. 政 治 大. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. iv. i n U. v. DOI:10.6814/NCCU201901066.
(7) Abstract This thesis examines the use of quantum Monte Carlo simulated anneal ing in the study of finitesize scaling for equilibrium phase transitions. For. 治 政 大 with z = ∞, and the exponent z = 1, a disordered quantum critical point 立 KosterlitzThouless (KT) transition approached through various annealing. quantum Ising models, we study quantum critical points with the dynamic. ‧ 國. 學. protocols in quantum Monte Carlo simulations using the stochastic series ex pansion method and a zerotemperature projector method. We demonstrate. ‧. that annealing by decreasing a transverse field at zero temperature, or by de. Nat. sit. y. creasing the temperature and the transverse field simultaneously can correctly. er. io. capture the critical scaling behaviors at z = 1 quantum critical points and the KT transition, if the rate of change is sufficiently slow. However, the z = ∞. al. n. v i n quantum critical point is anCannealing and our approaches fail to U h e n gbottleneck i h c reach the quasistatic limit of the random quantum Ising chain.. Keywords: simulated annealing, stochastic series expansion method, zero temperature projector method, random quantum Ising chain, triangular Ising antiferromagnet, KosterlitzThouless transition. v. DOI:10.6814/NCCU201901066.
(8) 立. 政 治 大. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. vi. i n U. v. DOI:10.6814/NCCU201901066.
(9) Contents 致謝. i. 摘要. 立. Abstract. ‧ 國. 模型概述. 1. 量子相變 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 1.2. 量子易辛模型 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 2. 1.3. 有限尺度定標 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 4. 1.4. 無序效應 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. y. sit. n. er. io. al. Ch. engchi. i n U. v. 隨機級數展開量子蒙地卡羅方法. 1. 7 11. 隨機級數展開法之推導 . . . . . . . . . . . . . . . . . . . . . . . . . .. 11. 2.1.1. 局域更新 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 15. 2.1.2. 叢集更新 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 19. 2.2. 無序系統的 SSE 方法 . . . . . . . . . . . . . . . . . . . . . . . . . . .. 20. 2.3. 零溫投射量子蒙地羅方法 . . . . . . . . . . . . . . . . . . . . . . . . .. 22. 基態投射法 . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 23. 2.1. 2.3.1 3. vii. 1.1. Nat. 2. v. ‧. 1. iii. 學. Contents. 政 治 大. 模擬結果 3.1. 27. 自旋鐵磁鏈及方晶格鐵磁 . . . . . . . . . . . . . . . . . . . . . . . . .. 27. 平衡態模擬 . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 27. 3.1.1. vii. DOI:10.6814/NCCU201901066.
(10) 模擬退火與 KibbleZurek 機制 . . . . . . . . . . . . . . . . . .. 30. 3.1.3. 模擬退火 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 34. 3.1.4. 零溫模擬退火 . . . . . . . . . . . . . . . . . . . . . . . . . . .. 46. 3.2. 三角反鐵磁 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 47. 3.3. 無序自旋鏈 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 52. 3.3.1. 有限溫度模擬退火 . . . . . . . . . . . . . . . . . . . . . . . .. 52. 3.3.2. 零溫模擬退火 . . . . . . . . . . . . . . . . . . . . . . . . . . .. 57. 結論. 61. 參考文獻. 立. 政 治 大. 63. ‧. ‧ 國. 學. io. sit. y. Nat. n. al. er. 4. 3.1.2. Ch. engchi. viii. i n U. v. DOI:10.6814/NCCU201901066.
(11) Chapter 1 模型概述. 立. 量子相變. 學. ‧ 國. 1.1. 政 治 大. 量子相變（quantum phase transition，QPT）描述了由量子起伏（quantum fluc tuation）導致量子多體系統在基態上的驟變 [1]。與古典相變不同之處在於量子. ‧. 相變發生在絕對零溫，由與溫度無關的量子起伏之參數 λ 所引發，例如：以量. sit. y. Nat. 子易辛模型為例，垂直於易辛軸的橫磁場。藉由橫磁場的變化，令量子易辛模. io. er. 型從有磁性的鐵磁體（ferromagnet）受到量子擾動的破壞而變成磁性消失的順磁 體（paramagnet）。磁性消失的相變為連續（二階）相變，為沒有相共存（phase. al. n. v i n Ch coexistense）的相變。其特徵在於關聯長度（correlation length）ξ 會在臨界點發散： engchi U ξ ∝ |λ − λc |−ν. (1.1). |λ − λc | 定義至量子臨界點 λc 之距離，指數 ν 則是關聯長度臨界指數且 ν > 0。 所謂關聯長度，是指熱力學系統某處的狀態發生變化時，最遠可影響某一處的 狀態變化的距離。而發散的關聯長度代表著系統每一處的狀態均緊密關聯 ( 無論 空間上的的實際距離 ) ，也意味著系統的結構並不重要，也因此產生了普適性 ( universality ) 的概念。不同的系統的臨界點只要能對應一系列相同的臨界指數，即 歸屬於同一的普適類（universality class） 。 系統基態至最低激發態的能隙（energy gap）∆ 隨著靠近量子臨界點而呈冪次. 1. DOI:10.6814/NCCU201901066.
(12) 方消失： ∆ ∝ |λ − λc |zν. (1.2). 這裡定義另一臨界指數 z，稱為動力學指數 ( dynamic exponent ) ，根據式 (1.1) 以 及式 (1.2) 可得： ∆ ∝ ξ −z. (1.3). 又弛豫時間 ( relaxtion time ) 與能隙成反比關係 ξτ ∝ ∆−1 ，可得： ξτ ∝ ξ z ,. (1.4). 政 治 大. 這意味在量子相變處時間與空間的關係是密不可分的。. 學. 1.2. ‧ 國. 立. 量子易辛模型. ‧. 在此介紹本論文將所探討的一典型用以描述磁性物質及量子相變的模型：量子. σ ˆiz σ ˆjz − h. X. ⟨ij⟩. n. al. X. Ch. i. engchi. σ ˆix ,. sit. io. ˆ = −J H. er. Nat. y. 易辛模型。其哈密頓（Hamiltonian）算符定義為：. i n P U. σ ˆ z,x 為庖立矩陣（Pauli matrix）之 z 或 x 的分量，而. (1.5). v. 的下標 ⟨ij⟩ 表示為只考慮. 相鄰的晶格點， J 為相鄰兩自旋的鍵結交互作用強度， J > 0 為鐵磁性； J < 0 則為反鐵磁性，而 h 為垂直易辛軸（z 軸）之橫磁場強度。本論文將探討的易辛 模型為一維鐵磁自旋鏈、二維方晶格鐵磁、二維的三角反鐵磁模型。上述鐵磁系 統於絕對零溫均存在一介於鐵磁態及順磁態的臨界點，其中二維模型的鐵磁態延 伸至有限溫度（T > 0）直至居里溫度（Curie temperature）。三角反鐵磁模型具有 所謂幾何挫折性（geometric frustration），因為三角晶格任意兩邊的相鄰自旋雖然 可以滿足反鐵磁性的排列，但第三個自旋無論朝上或朝下均會與其中一相鄰自旋 一起違背反鐵磁性的排列。不同於建立在簡單方晶格的模型，三角反鐵磁模型具 有豐富的相態。. 2. DOI:10.6814/NCCU201901066.
(13) 0.5 Tc,2(h) Tc,1(h) QCP. PM phase 0.4 0.3. T/J. KT phase. 0.2 0.1. Clock phase. 0 0. 0.2. 0.4. 0.6. 1. 0.8. h/hc. 立. 政 治 大. ‧. ‧ 國. 學. Figure 1.1: 三角反鐵磁隨溫度 T 及橫場 h 變化之相態圖。KosterlitzThouless (TK) 相變由兩弧形的 相邊界（虛線與實線）標出，這是根據 [2] 的式 (46) 及 Tc,1 = 4/9Tc,2 所描繪出來的。這兩相變線 界線界定低溫具磁性的時鐘態，兩線間的 TK 相及高溫的順磁態。兩條相邊界交會於量子臨界點 T = 0, h/hc = 1 ，及無量子起伏的古典 h = 0 處。. 三角反鐵磁模型在沒有橫場的作用下（h = 0）為無量子擾動的古典模型，因為. sit. y. Nat. 具有幾何挫折性所以在任何溫度（T > 0）皆呈現無磁性的無序態 [3]，在絕對零 D. E. er. io. z 溫時呈臨界態，這裡基態的自旋關聯函數（correlation function）C (r) ≡ σiz σi+r. 隨距離 r 逐漸以 r−2 冪次方遞減 [2]。絕對零溫時，在弱橫場的作用下基態反而呈. al. n. v i n C h h ，零溫基態的磁序將被破壞而變為順磁體。 現磁序 [4]，直至強度大於一臨界值 engchi U 零溫基態的磁有序態稱為時鐘相態（clock phase），因其序參數（order parameter） c. 的分佈呈六點旋轉對稱；偵測此相態的序參數可定義為一複數的 XY 序參數 [4]： √ meiθ ≡ m1 + m2 ei(4π/3) + m3 ei(−4π/3) / 3. (1.6). 其中 mα , α = 1, 2, 3 為三角晶格上三子晶格的個別 z 軸磁化量期望值。三角反鐵磁 模型在有限溫度下已知具有三種相態 [2, 4, 5]，分別為處於高溫的順磁態、中溫的 KosterlitzThouless（KT）相態，以及低溫的有序時鐘相態。其中 KT 相態與高溫 無序態及低溫的有序態分別由兩條在 h − T 相圖中的相邊界區分（圖1.1）。KT 相 態為一臨界相態（critical phase），在此相態中關聯函數 C (r) 會隨距離 r 以冪次方 的形式遞減： C (r) ∝ r−η ，而臨界指數 η 又會隨不同溫度而變化 [6]。在低溫的 3. DOI:10.6814/NCCU201901066.
(14) KT 有序相變 Tc,1 (h) 處， η = 1/9 ；在高溫無序相邊界 Tc,2 (h) 上， η = 1/4 ，與 二維的 XY 模型相同 [7]。而在量子臨界點處則為三維 XY 普適類的 [2, 4]。. 1.3. 有限尺度定標. 多體系統的相變發生在熱力學極限下，也就是當系統的粒子數 N = ∞ 時 。 而在計算模擬上，我們一般只能考慮有限尺度的系統。並藉有限尺度定標分析 （finitesize scaling analysis）將觀察量推導至熱力學極限 [8–10]。 在熱力學極限下，關聯長度 ξ 發散於臨界點（式 (1.1)），許多物理量隨與臨界. 政 治 大. 點距離 δ 呈現冪次方消失（或發散）的行為:. 立. (1.7). 定義了相關的臨界指數 κ 。結合式 (1.1) ，上式可寫為：. 學. ‧ 國. O ∼ |δ|+κ. ‧. (1.8). sit. y. Nat. O ∼ ξ −κ/ν. n. al. er. io. 在有限的尺度下，當 ξ → L 時，觀察量 O 會快速地收斂至某一值 O∗ ，且 O∗ 與系 統尺度 L 呈下式關係：. C hO. e(L)n ∼g Lc h i. ∗. −κ/ν. i n U. v. (1.9). 值 O∗ 對應的臨界點位置偏離熱力學極限下的真實臨界點位置，此偏差距離也與 系統尺度 L 有關： |δL | ∼ L−1/ν .. (1.10). 另外，量子臨界點發生於絕對零度 β ≡ 1/T = ∞，所以我們也需知道在有限溫度 所量測的觀察值如何推至零溫極限。事實上，溫度倒數 β 可視為量子系統的另一 個維度，又稱虛數時間；原因在於固定溫度下求觀察量期望值所需的 Boltzmann 因子（Boltzmann’s factor）e−β H ，若將其中 β 以虛數時間 it 替換，其形式與時間 ˆ. 演化算符 e−itH 相同。如此，虛數時間維度的長度取決於溫度倒數 β 值。有限系 ˆ. 統及有限溫度下所測量的物理量如何逼近熱力學極限及零溫極限的關係由式 (1.2). 4. DOI:10.6814/NCCU201901066.
(15) 定義的動力學指數透過 β ∼ Lz. (1.11). 描述。根據定標化理論 [9, 10]，我們可將式 (1.9)、式 (1.10) 及式 (1.11) 更廣義地 寫成以下定標形式： . . ˜ δL1/ν , β/Lz , O (δ, L, β) = L−κ/ν O. (1.12). ˜ 稱為定標函數（scaling function）。式 (1.12) 也可改寫為 O. 政 治 大. ˜ (ξ/L, β/Lz ) . O (δ, L, β) = L−κ/ν O. 立. (1.13). 在計算上，若 z 值已知，我們可固定 β/Lz 比值：β/Lz = constant；或在未知 z 值. ‧ 國. 學. 時，對固定系統尺度 L 針對不同 β 檢驗測量值的收斂情形。如此我們可在作定標 分析時不考慮定標函數的第二個變量。. ‧. io. n. al. sit. . m (δ, L) ∼ L−βm /ν m ˜ δL1/ν ,. (1.14). er. Nat. . y. 以本論文的觀察量磁化量為例，其定標形式可寫為：. i n U. v. βm 為磁化量的臨界指數1 。為了減少有限尺度定標分析涉及太多待定的臨界指. Ch. engchi. 數，我們在計算上常先分析一無因次的量，使得式 (1.12) 中的指數 κ = 0，常用 的無因次量為序參數的所謂 Binder 比值（Binder ratio）或 Binder 累積量（Binder cumulant），其為磁化量（序參數） m 的二階矩（secondorder moment）與四階矩 （fourthorder moment）的組合 [11, 12]： ". #. 3 ⟨m4 ⟩ g= 1− , 2 3⟨m2 ⟩2. (1.15). 式中的括弧符號 ⟨x⟩ 表示為 x 的期望值。根據上式的定義，在遠離臨界點的有序 態得 g → 1，因為這裡的磁化量分佈為位於 m = ±1 的雙峰型；而遠離臨界點的 無序態，磁化量分佈傾向以零為中心的高斯分佈（常態分佈），對 ⟨m2 ⟩ 和 ⟨m4 ⟩ 做 1. 我們用 βm 來取代磁化量慣用的臨界指數符號 β，以與溫度倒數 β 作區別。. 5. DOI:10.6814/NCCU201901066.
(16) 高斯積分，可得 Binder 比值 g → 0 。由於四階矩因次與二階矩因次的平方抵銷， 故無因次量的定標形式具簡單的形式： . . g = g˜ δL1/ν .. (1.16). 如此一來，對應不同的 L 的 Binder 比值對 δ 作圖，便會交於 δ = 0 處，即在臨界 點時，Binder 比值與尺度 L 大小無關。 對於三角反鐵磁，我們好奇的是它在 KT 相態的性質，根據 [4, 13] 我們觀測其 磁化量平方以及磁化率（susceptibility）。將磁化率定義為：. 政 治 大 D. 立. E. χ = L2 m2 /T .. ‧ 國. χ ∝ ξ 2−η ,. m ∝ ξ −η/2. (1.18). ‧. ξ ∝ e−at ,. 學. 藉由 KT 理論及定標分析理論可知 [13, 14, 14]：. (1.17). Nat. sit. y. 這裡 a 為非普適性常數（nonuniversal constant）， t = (T − Tc ) /Tc 為約化溫度. er. al. n. 形式表示為：. io. （reduced temperature），且 Tc 為臨界溫度。故可將磁化量與磁化率的有限尺度定標. Ch. engchi. m (ξ, L) = L−η/2 m ˜ KT (ξ/L) ,. i n U. v. χ (ξ, L) = L2−η χ˜KT (ξ/L) .. (1.19). 帶入關聯長度 ξ 與約化溫度 t 的關係，我們可將式 (1.19) 重寫為 . m (L) = L−b m ˜ KT L−1 e−at. . χ (L) = Lc χ˜KT L−1 e−at. . . (1.20). (1.21). 這裡我們參考 [4,13] 定義兩個指數 b = η/2 和 c = 2 − η 。式 (1.20) 適用於 T < Tc,1 而式 (1.21) 適用於 T > Tc,2 。至於磁化率（定義於式 (1.17)）的定標形式在零溫. 6. DOI:10.6814/NCCU201901066.
(17) 下量子相變點為 . χ (δ, L) = L2−2βm /ν χ˜ δL1/ν. . (1.22). 這可從磁化率的定義（式 (1.17） ）及磁化量的定標形式 (1.14) 看出。這個量子相變 點屬 3DXY 普適類，對應的臨界指數為 βm ≈ 0.3485 、 ν ≈ 0.67155 [15]。. 1.4. 無序效應. 通常實際的物質難免都會有雜質或是缺陷。這些不均質的性質使得無序（非均 質）系統和有序（均質）系統具有不同的物理性質，尤其於臨界點時兩者之間差. 政 治 大. 異性會十分明顯。在這裡我們只關注與時間無關的雜質或缺陷，也聚焦於無序交. 立. 互耦合的效應。. ‧ 國. 學. 在均質系統中，熱力學觀察量（例如：磁化率）只在相變點發生奇點，偏離相 變點奇點則消失；然而，R. B. Griffiths [16] 以無序的古典易辛模型作為例子，提. ‧. 出在非均質的情形下奇點可能會發生在非臨界的區域，原因在於即使是非均質的 系統，其順磁相內亦隨機存在局域有序鐵磁性的區塊，這些區塊產生的機率很小，. y. Nat. sit. 卻能使磁化率產生奇點，我們稱這個狀態為 Griffiths phase ，而在非臨界區域發生. n. al. er. io. 奇點的現象稱為 Griffiths singularities 。雖然此現象在理論上可被預期，但在古典. i n U. v. 的系統並不明顯，且在實驗上也難以被實現。同樣的無序效應在量子相變點附近. Ch. engchi. 卻能十分明顯，且在實驗上已被實現出 [17–19]。這是因為在局域鐵磁性區塊產生 極小的能隙，一體積為 V = ℓd 的區塊之能隙為 ε ∼ exp(−aℓd )，在量子系統這意 味極長的弛豫時間，而造成相對大的動力學指數 z。又這些局域區塊發生的機率 極小，且機率隨區塊體積呈指數下降：Pℓ ∼ exp(−bℓd )，這個極小機率加上亦為指 數小的能隙造成一冪次方的態密度（density of states）ρ(ε) ∼ ε−1+d/z 。冪次方態密 度進而可造成一些發散的物理量，如磁化率 χ(T ) ∼ T d/z−1 發散在 z > d 及 T = 0 情況下 [20]。 在這裡我們將探討無序的一維量子易辛鐵磁模型。在系統中相鄰自旋的鍵結 交互作用強度都不盡相同，且每個晶格點受到的橫場強度也不同。我們將哈密頓 （式 (1.5)）改寫為： ˆ =− H. X. z − ˆi+1 Ji σ ˆiz σ. X. hi σ ˆix. (1.23). i. i. 7. DOI:10.6814/NCCU201901066.
(18) 其中我們將自旋交互作用強度 Ji 設為均勻介於 (0, 1] 間的亂數，而橫場 hi 為均勻 介於 (0, h] 的亂數。. h1. h2. J1. h3. h4. J2. h5. J3. J4. Figure 1.2: 無序一維量子易辛鏈示意圖。紅色圓點表示自旋，水藍色長條大小表示自旋所受橫場的 強度，而相鄰自旋的鍵結交互強度大小則以 Ji 表示。. 立. 政 治 大. 量子易辛鐵磁模型具顯著的 Griffiths 效應，且動力學指數 z 隨愈靠近量子臨界. ‧ 國. 學. 點而增大，對定義於式 (1.23) 的一維系統而言，量子臨界點發生於 [21, 22]. ‧. ln h = ln J ,. (1.24). Nat. sit. y. 這裡 x 表示求無序樣本的鍵結/橫場強度 x 平均，即所謂的無序平均（disorder. er. io. average） 。根據強無序重整群（strongdisorder renormalization group）的分析 [21,22]，. al. n. 動力學指數在量子臨界點可達到無窮大 z = ∞，這也是說弛豫時間與式 (1.4) 不. i n C 同，不再是冪次方的成長，而是呈指數成長 h e n g c[23]： hi U ξτ ∝ eα. v. √. ξ. (1.25). 其中 α 為常數。對有限尺度系統而言，上式的關係式可換成虛數時間的長度 β 與 系統尺度 L 的關係： ′. β ∝ eα. √. L. (1.26). 如此，對於無序系統的有限尺度定標形式可廣義地寫為： . ˜ δL1/ν , ln (β) /L1/2 O (δ, L, β) = L−κ/ν O. 8. . (1.27). DOI:10.6814/NCCU201901066.
(19) 根據式 (1.24) 這裡定義控制參數到臨界點的距離為 [21, 22] δ ≡ ln h − ln J .. (1.28). 由強無序重整群分析 [21, 22]，已知一維無序量子易辛模型的關聯長度及磁化量相 √ 關臨界指數 ν = 2、βm = 2 − (1 + 5)/2 ≈ 0.38，這裡的關聯長度及磁化量均是無 序平均量。. 立. 政 治 大. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. 9. i n U. v. DOI:10.6814/NCCU201901066.
(20) 立. 政 治 大. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. 10. i n U. v. DOI:10.6814/NCCU201901066.
(21) Chapter 2 隨機級數展開量子蒙地卡羅方法. 治 政 本章節將描述隨機級數展開量子蒙地卡羅方法（Stochastic Series Expansion 大 立 quantum Monte Carlo method）[24]，簡稱 SSE ，這是一個高效率且被廣泛應用的 ‧ 國. 學. 於有限溫度系統的量子蒙地卡羅演算法。此方法的基礎建立於 1960 年代針對海森 堡鐵磁模型的 Handscomb 方法 [25]。主要是由 Anders Sandvik 將之改良並發展至. ‧. 目前的通用形式 [24, 26]，可廣泛應用於量子多體問題中各種晶格模型。此章節將. sit. y. Nat. 簡介 SSE 的方法以及應用在本論文所探討的 1/2 自旋的量子易辛模型（Quantum. io. er. Ising model）[27, 28]。這裡我們也同時簡介與 SSE 方法息息相關但用以探討絕對 零溫基態問題的零溫投射法。. n. al. 2.1. Ch. engchi. i n U. v. 隨機級數展開法之推導. D E. . . ˆ 中物理量 O 的期望值為 O ˆ = Tr ρˆO ˆ ，其中密度矩陣（density 在量子系統 H matrix）ˆ ρ 在固定溫度 T = β −1 下的表示式為. ρˆ =. 1 ˆ Tr e−β H , Z. (2.1). 其中歸一化常數即所謂配分函數（partition function） ： . . ˆ −β Hˆ . Z = Tr Oe. (2.2). 11. DOI:10.6814/NCCU201901066.
(22) 期望值因而表示為： ˆ −β Hˆ . ˆ = 1 Tr Oe O Z. D E. (2.3). 我們可利用泰勒展開式（Taylor expansion）將 e−β H 展開： ˆ. e−β H = ˆ. ∞ X (−β)n. n!. n=0. ˆn H. (2.4). 若選擇 |x⟩ 為基底，例如 N 個 1/2 自旋模型的標準基底 |x⟩ = |σ1 σ2 · · · σN ⟩（其中 σi 為 σ ˆiz 的本徵向量 |↑⟩ 或 |↓⟩），配分函數可重新表示為：. Z=. ∞ X βn X D. . n. E. 治 n! 政 大 ˆ x| −H. n=0. 立. |x .. (2.5). {x}. 學. （α = 1, 2, · · · , n − 1），得： ∞ X βn X X n=0. n!. ···. {x0 } {x1 }. X D {xn−1 } |. ˆ 1 x0 |−H|x. ED. E. D. ˆ 0 ˆ 2 · · · xn−1 |−H|x x1 |−H|x {z. }. (2.6). sit. Nat. n個. E. y. Z=. ‧. ‧ 國. P 再利用基底 {|x⟩} 的完備性，將上式插入 n − 1 個單位算符 Iˆα = {xα } |xα ⟩ ⟨xα |. . n. er. io. 其中 |xα ⟩ 的下標僅標註所插入單位算符 Iˆα 的順序。為了滿足矩陣跡數（Trace） . . a l 兩 端 的 狀 態 ⟨x | 和 |x ⟩ i 需 v 相 同。 接 著， 我 們 將 n Ch 式 (1.5) 的哈密頓改寫成各項局域算符的組合： engchi U n. ˆ 的 形 式， 必 須 注 意 −H. 0. ˆ =− H. XX t. 0. ˆ t,a , H. (2.7). a. 其中 t 表示算符的種類，例如鍵結算符或橫場算符等；而 a 標示算符位於晶格模 型的位置，例如晶格位置 i 或鍵結位置 (ij)。將上式代入式 (2.6) ，將 Z 表示為. Z=. n D ∞ X βn X X Y n=0. n!. E. ˆ tα ,aα |xα , xα−1 |H. (2.8). {xα } {Sn } α=1. ˆ tα ,aα ）乘積的指標序列：Sn = [t1 , a1 ] , [t2 , a2 ] , · · · [tn , an ] 這裡 Sn 為一串 n 個局域算符（H 。 SSE 方法即採用隨機的方式（蒙地卡羅方法）更新算符序列 Sn 。一初始狀態 12. DOI:10.6814/NCCU201901066.
(23) |x0 ⟩ 將被算符序列 Sn 中的算符依序作用，故 Sn 可視為 D 維自旋模型的一額外的 維度，構成一 D + 1 維的模擬空間。在實際計算中，算符數目 n 為一浮動的值， 我們選擇一大於 n 的高冪次 M 來截斷展開式，並同時引入 M − n 個單位算符 Iˆ 於算符序列中，如此我們可有一固定長度的算符序列 SM 。M − n 個單位算符 Iˆ 可 有多種排入算符序列 SM 的方式，這些排列組合方式之數目為： M n. !. =. M! . n! (M − n)!. (2.9). 故我們需將配分函數調整為. Z=. 政 治 大. M D E β n (M − n)! X X Y ˆ tα ,aα |xα . xα−1 |H M! {xα } {SM } α=1. 立. (2.10). P. ⟨O⟩ =. Ox,SM Wx,SM {x,SM } Wx,SM. {x,SM }. P. ‧. ‧ 國. 學. 比照式 (2.3) 觀察量期望值，在以 |x⟩ 為基底的 SSE 方法中我們可將其改寫為：. Nat. al. n. Wx,SM =. M D E β n (M − n)! Y ˆ tα ,aα |xα . xα−1 |H M! α=1. er. io. sit. y. 其中 {xα , SM } 為狀態向量—算符的組態，對應的權重 Wx,SM 為：. Ch. engchi U. v ni. (2.11). 在 蒙 地 卡 羅 演 算 法 中， 我 們 將 依 權 重 Wx,SM 來 採 樣， 此 為 所 謂 重 要 性 採 樣 （importance sampling）。 ˆ tα ,aα |xα ⟩ 為負數將產生不合 由上式可明顯看出，若其中一個矩陣元素 ⟨xα−1 |H 理的負機率問題（sign problem ），使得演算法不適用。對於我們這裡所探討的模 型而言，負機率問題將不會出現。以下我們討論如何在量子易辛模型的哈密頓算 符加上一些常數項來避免負號機率。我們將式 (1.23) 中的量子易辛哈密頓算符改 寫為. ˆ =− H. X i. ˆ −1,i − H. X i. 13. ˆ 0,i − H. X. ˆ 1,(ij) H. (2.12). ⟨ij⟩. DOI:10.6814/NCCU201901066.
(24) ˆ t,a 分別為： 上式各項局域算符 H ˆ −1,i = hi σ H ˆix. (2.13). ˆ 0,i = hi H. (2.14). ˆ 1,(ij) = Jij σ H ˆiz σ ˆjz ± Jij .. (2.15). 如此在 (2.15) 加入了常數項使得對應之矩陣元素不為負，如: 第二項為加號為鐵磁（Jij > 0）鍵結項： E. D. D. E. 政 治 大. ˆ 1,(ij) |↑↑ = ↓↓|H ˆ 1,(ij) |↓↓ = 2Jij ↑ ↑ |H. 立. (2.16). D. E. D. 學. ‧ 國. 第二項為減號為反鐵磁（Jij < 0）鍵結項： E. ˆ 1,(ij) |↑↓ = ↓↑|H ˆ 1,(ij) |↓↑ = −2Jij ↑ ↓ |H. (2.17). ‧. ˆ 0,i ，將使各晶格點算符之各非零矩陣元素的值相同，這方便我 又引入常數算符 H. y. Nat. D. E. al. ˆ −1,i |↓ ↑|H. D. E. D. E. ˆ −1,i |↑ = ↓|H. er. io. sit. 們作所謂叢集更新時方便使用（將後續討論） [28]。晶格點項非零矩陣元素如下：. n. v i n ˆC |↑ h e =n g↓| cHˆ h i|↓ U= h ↑| H. D. E. 0,i. = hi. 0,i. (2.18). i. ˆ 0,0 = Iˆ 表示單位算符。由上述各非零矩陣元素可見，僅橫場算符 另外我們以 H ˆ −1,i 為所謂的非對角算符 ( offdiagonal operator )；鍵結算符 H ˆ 1,(ij) 、常數算符 H ˆ 0,i 及單位算符 H ˆ 0,0 則為對角算符 ( diagonal operator )。 H 在 SSE 量子蒙地卡羅計算中我們主要對算符序列 SM 內各項算符（包含單位 算符）作抽樣及更新，藉此進而更新自旋組態。更新方式主要分成兩大類型，第 一類為依 Metropolis 法則 [29] 的局域組態的更新（local update），另一類為類似 SwendsenWang 方法 [30] 的叢集更新（cluster update）；前者主要針對對角算符 ˆ 0,i 和 H ˆ 1,(ij) 以及單位算符 H ˆ 0,0 的更新，後者則是針對非對角算符 H ˆ −1,i 的更新， H 運算時，每一步的蒙地卡羅步驟都包含著這兩種類型的更新。以下我們分別以兩 小節描述這類型的更新法則。. 14. DOI:10.6814/NCCU201901066.
(25) 局域更新. 2.1.1. 這裡我們以一具有 N 個自旋、 Nb 個鍵結之量子易辛鐵磁模型為例，且考慮所 有鍵結強度為定值 Jij = J > 0, ∀i, j 及橫場強度亦為定值 hi = h > 0 , ∀i 的情形。 我們從頭檢視算符序列 SM 中的每一個算符，並以下列幾點規則作算符序列的 更新： ˆ −1,i ，則更新該算符所對應之晶格位置 1. 若遇到一橫場項算符 (非對角算符) H 其作用後的自旋狀態，即自旋的翻轉。 ˆ 0,i 或 H ˆ 1,(ij) ，則以下列機率（其推導可看第 16 頁）移 2. 若遇到一對角算符 H. 政 治 大. 除它（即改變成一單位算符）：. 立 P. #. ". remove. M −n+1 ,1 = min β [hN + (2J) Nb ]. (2.19). ‧ 國. 學. 無論此移除動作是否被接受，我們都繼續檢驗算符序列 SM 中的下一個算. ‧. 符。. Nat. sit. n. er. io. ˆ 0,i ： (a) 以以下機率選擇一常數算符 H. al. C hP. . y. ˆ 0,0 ，則嘗試插入一對角算符。 3. 若出現一單位算符 H. i n hN U. . v. i + 2JN e nHˆg c =h hN. select. 0,i. (2.20) b. ˆ 1,(ij) ： 或以以下機率選擇插入一鍵結算符 H . . . . ˆ 1,(ij) = 1 − P H ˆ 0,i = Pselect H. 2JNb hN + 2JNb. (2.21). ˆ 0,i （或鍵結算符 H ˆ 1,(ij) ）在系統哈密頓 上述機率呈現的正是常數算符 H 算符（式 (2.12)）中所有的對角算符出現的比例及矩陣元素大小的比重 （式 (2.18) 及式 (2.17)）。 (b) 一旦選擇好預加入的算符類型，以下列機率（其推導可看第 18 頁） ： ". Pinsert. #. β [hN + (2J) Nb ] = min ,1 M −n 15. (2.22). DOI:10.6814/NCCU201901066.
(26) 考慮接受加入此算符，並隨機選擇置入的晶格位置；若是接受的算符 ˆ 1,(ij) ，則置入的鍵結位置 (i, j) 相鄰的兩自旋必須為平行 為鍵結算符 H （相同方向）狀態，否則拒絕此鍵結算符的加入。 上述的步驟必須重複進行至有一對角算符被置入為止。 公式推導 以下我們推導式 (2.19) 中的接受機率。首先需滿足所謂的細緻平衡（detailed balance）的條件： P (A → B) W (A) = P (B → A) W (B) ,. 立. 政 治 大. (2.23). 這裡 P (A → B) 為從 A 組態轉移至另一 B 組態的機率， W (A) 為組態 A 的權重. Pselect (B) 及接受這個轉移的機率 Paccept (A → B)，即：. (2.24). sit. y. Nat. P (A → B) = Paccept (A → B)Pselect (B) .. ‧. ‧ 國. 學. （式 2.11）。其中轉移機率 P (A → B) 又可分成兩部份，分別為選擇組態 B 的機率. er. io. 雷同的符號定義適用於反向的 P (B → A)。如此，式 (2.23) 可改寫為. n. al. Paccept (A → B)Pselect. i n C (B)W (A) A)P h e=nPg c(Bh → i U. v. select (A)W (B) .. accept. (2.25). 根據 Metropolis 的更新法則 [29]，以下的接受機率滿足上述細緻平衡條件： ". P accept. #. W (B) Pselect (A) (A → B) = min ,1 W (A) Pselect (B). (2.26). 針對我們的計算，我們考慮以下機種情況： ˆ 0,i ）並以一單位算符取代 1. 組態 A 至組態 B 的轉移為移除一個常數算符（H ˆ ˆ 表達 W (B), 而 W (H ˆ 0,i ) 表達轉移前的 （I）的對應組態轉移。故我們以 W (I) 組態 A 權重。則接受這個移除的機率為 ". ˆ 0,i → I) ˆ = min Premove (H 16. #. ˆ Pselect (H ˆ 0,i ) W (I) ,1 ˆ 0,i ) Pselect (I) ˆ W (H. (2.27). DOI:10.6814/NCCU201901066.
(27) 權重 W 相關比值可參考式 (2.11) 得出: D. β n−1 (M −(n−1))! QM ˆ α=1 xα−1 |Htα ,aα |xα M! D E β n (M −n)! QM ˆ ′ ′ x | H |x ′ ′ ′ t ,a α α α α =1 α −1 M!. ˆ W (I) = ˆ 0,i ) W (H. D. E. E. ˆ 0,0 |xα M − n + 1 xα−1 |H D E = ˆ 0,i |xα′ β xα′ −1 |H =. (2.28). M −n+11 , β h. ˆ 牽涉到 n − 1 非單位算符；第二個等式是將相同 其中第一個等式中的 W (I) 的矩陣元素抵消後的結果。選擇組態的機率則為：. 立. 政 治 h大 , hN + (2J) N. ˆ 0,i ) = Pselect (H. b. (2.29). ˆ = 1. Pselect (I). ‧ 國. 學. 故式 (2.27) 之取代常數算符的接受機率可得：. ‧. ". al. n 此正為式 (2.19)。. (2.30). er. io. sit. y. Nat. ˆ 0,i Premove (H. #. h ˆ = min M − n + 1 1 · ,1 → I) β h hN + (2J) Nb " # M −n+1 = min ,1 β [hN + (2J) Nb ]. Ch. engchi. i n U. v. ˆ 1,(ij) 的情形： 2. 我們再看以單位算符取代鍵結算符 H . ˆ 1,(i,j) Premove (H. . ˆ ˆ ˆ = min W (I) Pselect (H1,(ij) ) , 1 → I) ˆ 1,(ij) ) Pselect (I) ˆ W (H. (2.31). 其權重比值： ˆ W (I) = ˆ 1,(ij) ) W (H. D. β n−1 (M −(n−1))! QM ˆ α=1 xα−1 |Htα ,aα |xα M! D E β n (M −n)! QM ˆ ′ =1 xα′ −1 |Ht ′ ,a ′ |xα′ α α α M!. D. E. ˆ 0,0 |xα M − n + 1 xα−1 |H D E = ˆ 1,(ij) |xα′ β xα′ −1 |H =. E. (2.32). M −n+1 1 β 2J. 17. DOI:10.6814/NCCU201901066.
(28) 第二個等式也是將相同的矩陣元素抵消的結果。又選擇機率鍵結算符的機率 ˆ 1,(ij) ) 為： Pselect (H ˆ 1,(ij) ) = Pselect (H. 2J . hN + (2J) Nb. (2.33). 故式 (2.31) 取代鍵結算符的接受機率與上述 1. 中取代常數算符的接受機率 一樣，同為式 (2.19)。也就是說，被取代的對角算符種類並不會影響其被接 受移除的機率。 ˆ t,a , t = 0/1）被接受的機率，也 3. 接著考慮以對角算符取代單位算符（Iˆ → H 就是非單位算符數目提升 (n → n + 1) 的情況。. 治 ˆ ) P 政 W (H ˆ ˆ (I → H ) = min ˆ 大 W (I) P 立 ". Pinsert. t,a. t,a. (2.34). 學. 其權重比值為：. β n+1 (M −(n+1))! M! β n (M −n)! M!. ˆ t,a |xα ⟩ ⟨xα−1 |H ˆ 0,0 |x′α ⟩ ⟨xα′ −1 |H ˆ t,a |xα ⟩ β ⟨xα−1 |H = ˆ 0,0 |x′ ⟩ M − n ⟨xα′ −1 |H. ˆ t,a ) W (H = ˆ W (I). Nat. α. n. er. io. al. ˆ t,a |xα ⟩ β ⟨xα−1 |H . = M −n 1. sit. y. ‧. ‧ 國. #. ˆ ,1 . ˆ select (Ht,a ) select (I). Ch. engchi. i n U. v. D. E. ˆ 0,i |xα = 若是以常數算符取代單位算符，上式最後等式中的矩陣元素為 xα−1 |H D. E. ˆ 1,(ij) |xα = 2J h ；而若以鍵結算符取代單位算符，則此矩陣元素為 xα−1 |H ˆ 0,i ) （式 (2.29)）或選擇鍵結算符的機率 。帶入選擇常數算符的機率 Pselect (H ˆ 1,(ij) )（式 (2.33)），式 (2.34) 的接受插入對角算符的機率成為： Pselect (H ". #. ˆ t,a ) = min β [hN + (2J) Nb ] , 1 . Pinsert (Iˆ → H M −n. (2.35). 這個結果不分插入的算符是常數算符或鍵結算符。. 18. DOI:10.6814/NCCU201901066.
(29) Figure 2.1: 此示意圖為自旋數目 N = 4 及算符序列長度 M = 10 的 SSE 方法模擬空間。箭頭方 向為算符序列和自旋狀態之傳遞方向，圓點為自旋狀態，其中實心圓點為自旋朝下，而空心圓點 為自旋朝上，左右兩端自旋狀態相同滿足了矩陣跡數的形式。方塊則是為算符總類，空心正方形 ˆ 0,i ，實心正方形為橫場算符 H ˆ −1,i ，而橫跨兩自旋的實心長方形則是為鍵結算符 為常數算符 H ˆ H1,(ij) 。. 立. ‧ 國. 學. 2.1.2. 政 治 大. 叢集更新. ‧. 所謂的叢集是指在 D + 1 維度的模擬空間中一群自旋與算符的集合。我們以下. y. Nat. 列方式建構一叢集：從一算符出發，沿晶格空間的方向及算符傳遞方向納入其它. al. er. io. sit. 自旋與算符於叢集內，並以下述兩個法則決定其成長與停止：. n. ˆ −1,i 或常數算符 H ˆ 0,i 。 1. 叢集終止於晶格點算符，即橫場算符 H. Ch. engchi. i n U. v. ˆ 1,(ij) ，四端連結之自旋會被納入同一叢集，且會沿著四端 2. 若遇到鍵結算符 H 方向繼續成長（如圖）。 當叢集建構完畢後，我們以 1/2 的機率翻轉叢集內所有的自旋狀態與算符種類， 其翻轉步驟包含： • 自旋翻轉。 ˆ −1,i 變常數算符 H ˆ 0,i ；常數算符 H ˆ 0,i 變橫場算符 H ˆ −1,i 。 • 橫場算符 H ˆ 0,i 其矩陣元素值與 兩算符之所以可以直接互換，這是因為我們加入的常數算符 H ˆ −1,i 的矩陣元素值相同；又因為矩陣元素值相同，使得變換前後狀態 橫場算符 H 的權重和機率相同，且翻轉機率都是 1/2 ，滿足了細緻平衡。故我們可以在叢集 內將兩算符直接互換。 19. DOI:10.6814/NCCU201901066.
(30) 政 治 大. 立. ‧ 國. 學. (a) 建構叢集. (b) 叢集更新. ‧. Figure 2.2: SSE 方法的叢集更新示意圖。虛線框住的範圍包還有在方框邊界的算符都屬於同一叢 ˆ 1,(ij) 所連結的四端自旋並沿著四端成長方向的自旋也納入同一叢集。圖 ( a )：建 集，鍵結算符 H 構了三個叢集。圖 ( b )：以 1/2 的機率更新叢集後（最右邊的叢集則未成功翻轉），叢集內自旋翻 ˆ −1,i 跟常數 轉（實心圓點變空心圓點，空心圓點變實心圓點），算符種類也跟著轉換（橫場算符 H ˆ ˆ 算符 H0,i ） ，以此更新了非對角算符 H−1,i 。. er. io. sit. y. Nat. a l 方法 無序系統的 SSE n. 2.2. Ch. engchi. i n U. v. 基本上無序系統的 SSE 方法與有序系統的 SSE 方法步驟都相同，但與有序系 統不同的是，無序系統的每個相鄰自旋鍵結強度以及自旋所受的橫場強度不盡相 同（常數算符的值是對應同位置的橫場強度）。故我們必須考慮在進行算符序列的 更新上，無論是接受算符的移除或置入的機率，或者是置入哪種算符的機率，都 是否會因為矩陣元素的值不同而導致機率必須修正? 我們可以參考第 16 至 18 頁 的式 (2.27 2.35) 的公式推導，可推得： 1. 移除對角算符之機率： ". Premove. #. ˆ Pselect (H ˆ t,a ) W (I) = min ,1 ˆ t,a ) Pselect (I) ˆ W (H " # M −n+1 = min ,1 β [hsum + 2Jsum ] 20. (2.36). DOI:10.6814/NCCU201901066.
(31) 其中 hsum =. P i. hi 、 Jsum =. P. ⟨ij⟩. Jij ，因為這裡場及鍵結強度已不是定值，. 故分別得改寫為隨晶格位置變化的值： hi 、 Jij 。 2. 接受對角算符置入之機率： ". #. ˆ ˆ ˆ t,a ) = min W (Ht,a ) Pselect (I) , 1 Pinsert (H ˆ Pselect (H ˆ t,a ) W (I) # " β [hsum + 2J sum ] ,1 = min M −n. (2.37). 另外，當我們以對角算符取代常數算符時，由於各算符的矩陣元素值不盡相同， ˆ 0,i 或是鍵結算符 H ˆ 1,(ij) 的機率也不能再只有分成兩個比 選擇置入常數算符 H. 政 治 大. 重 ( 式 (2.20)、式 (2.21) )；而且給各算符不同的值的同時，也代表著給定其位置. 立. 了，故不能再跟均質系統相同，先決定好算符種類再隨機安排位置。因此我們依. ‧ 國. 學. 照每個對角算符的權重大小來決定置入算符的種類 [27]：. ‧. 1. 首先，先將對角算符賦予編號 (1, 2, · · · M ) ，並將其權重 P (i) 歸一化且累加. n. al. i. P (i) P (i). Ch. sit. io. Pc (k) = PMi. (2.38). er. Nat. Pk. y. 起來，形成累積機率 Pc 並暫存在一個表中：. i n U. v. 其中 k = [1, 2 · · · , M ] ， Pc (k) 為對應編號 k 的對角算符之累積機率。. engchi. 2. 接著我們產生一隨機數字 0 ≤ R ≤ 1 ，然後從表中尋找最小的編號 k 使得 Pc (k) ≥ R 。而找尋最小的編號 k 的方式我們使用一個叫做二元搜尋 ( binary search）的方法。 所謂二元搜尋（binary search）也稱作區間分半搜尋（halfinterval search） ，是一種 在有序的陣列中尋找某一特定元素的搜尋方法。搜尋的過程是由中間元素開始， 若中間元素是我們要找的特定元素則停止搜尋；若是特定元素大於或小於中間的 話，則在大於或小於的那一半搜尋，而且也是從那一半的中間元素比較，重覆步 驟直到找到特定元素為止。就好像遊戲終極密碼一樣，猜錯數字就提示在比正確 數字大或小，不停的縮小範圍直到猜對數字為止。. 21. DOI:10.6814/NCCU201901066.
(32) 立. 政 治 大. Nat. io. sit. y. 零溫投射量子蒙地羅方法. er. 2.3. ‧. ‧ 國. 學. Figure 2.3: 二元搜尋算符之示意圖。橫數軸代表式 (2.38) 定義的累積機率 Pc (k), k = 1, 2, · · · 8，其 ˆ t ,a ⟩。隨機產生一介於 (0, 1] 的亂 間隔長度取決於各對角算符的比重，即其矩陣元素值的大小 ⟨H i i 數，例如 R = 0.57，在此示意圖中，可知 R 介於線段 [0.46, 0.66] 間，即對應編號 k = 5 的算符 ˆ 1,0 。二元搜尋過程中，先確定 R 落在全數軸的右半部，即 [Pc (4), Pc (8)] 那一區域，接著判斷它 H 落於那區域的左半側（Pc (4) < R < Pc (6)），依此進行至找到對應的區間及其算符。. 在第 1 章開始提到，量子相變是發生在絕對零度 T = 0 的情形下，其反應系統. n. al. Ch. i n U. 在基態的驟變。在基態物理量 O 的期望值表示為：. engchi. O=. v. 1 ⟨ϕ0 |O|ϕ0 ⟩ Z. (2.39). 其中 ϕ0 為基態的向量（波函數）。而歸一化常數 Z 則為： Z = ⟨ϕ0 |ϕ0 ⟩. (2.40). 繼討論有限溫度的 SSE 方法之後，我們這裡討論一解量子多體系統零溫狀態的蒙 地卡羅方法主要利用高次方的哈密頓算符來投射出基態向量 [28, 31, 32]。. 22. DOI:10.6814/NCCU201901066.
(33) 基態投射法. 2.3.1. 政 治 大. Figure 2.4: 此示意圖為自旋數目 N = 4 及算符數 2n = 10 的零溫投射方法模擬空間。箭頭方向為 算符乘積和自旋狀態之傳遞方向，圓點為自旋狀態，其中實心圓點為自旋朝下，而空心圓點為自 ˆ 0,i ，實心正方形為橫場算符 H ˆ −1,i ，而 旋朝上。方塊則是為算符總類，空心正方形為常數算符 H ˆ 橫跨兩自旋的實心長方形則是為鍵結算符 H1,(ij) 。與 SSE 不同的地方在於算符傳遞方向的兩端為 開放邊界條件，也就是左右兩端的狀態不必相同，而期望值矩陣元 ⟨xl |O|xr ⟩ 的計算沿著算符乘積 方向做測量並做平均。. 立. ‧. ‧ 國. 學. 首先取一任意狀態的向量 |Ψtr ⟩ 當試探態向量（trial state vector），令能量本徵. Nat. n. Ch. X. ck |ϕk ⟩. k. engchi U. er. io. al. |Ψtr ⟩ =. sit. y. 向量 {|ϕn ⟩ , n = 0, 1, 2, · · · } 為基底向量，試探態向量則可展開表示為：. (2.41). v ni. ˆ n 作用於試探態 其中 ck 為展開係數。而投射法是將一高次方（ n ）哈密頓算符 H 向量 |Ψtr ⟩ 上，並且給予一常數 C ，使得最低能量本徵值的絕對值 |E0 − C| 大於 其於本徵值的絕對值 [31]，因此可表示為： . ˆ C −H. n. ". |Ψ tr ⟩ = c0 (C − E0 ). n. . c1 C − E 1 |ϕ0 ⟩ + c0 C − E 0. n. #. |ϕ1 ⟩ + · · ·. (2.42). 由於在上式中的能量比值 |(C − Ek )/(C − E0 )| < 1 ，所以當 n → ∞ 下，其值趨 近於零，可得 . ˆ C −H. n. n→∞. |Ψ tr ⟩ −−−→ c0 (C − E0 )n |ϕ0 ⟩. (2.43). ˆ n 的作用下，任意不與基態向 若 c0 = ⟨ϕ0 |Ψtr ⟩ ̸= 0。故在高次方哈密頓算符 H 量垂直的試探態向量可趨近於基態向量。因此我們可將在基態的歸一化常數 23. DOI:10.6814/NCCU201901066.
(34) Z = ⟨ϕ0 |ϕ0 ⟩ 改寫為： D. . ˆ Z = Ψtr | C − H. n . ˆ C −H. n. |Ψtr. E. (2.44). 其中哈密頓算符以量子易辛模型為例，其形式如同式 (2.12) ，為晶格點算符與鍵 結算符的相加形式。為了簡化符號我們引入一算符乘積 P 來替代上式的展開項， . ˆ C −H. n. X. =. PˆSn , PˆSn =. n Y. ˆ t ,a H j j. (2.45). j=1. Sn. ˆ t.a 為哈密頓中的各項局域算符，定義於式 (2.13)(2.15)，而 Sn 為一串 n 算 其中 H. 政 治 大. 符乘積。式 (2.39) 的物理量 O 期望值在投射方法上可被改寫為：. 立. D. P. ‧ 國. ⟨O⟩ = P. Ψtr |PˆS†′ OPˆSn |Ψtr. D. ′ Sn Sn. n. Ψtr |PˆS†′ PˆSn |Ψtr. E. E. n. 學. ′ Sn Sn. (2.46). ‧. 接 著， 我 們 選 擇 |x⟩ 作 為 基 底， 並 在 式 (2.46) 引 入 基 底 完 備 性 構 成 的 單 位 算. y. E. io. 寫為. sit. ˆ tα ,aα |xα 、 ⟨xα−1 |O|xα ⟩ 。我們用上述定義的符號將歸一化常數 Z 進一步 xα−1 |H. er. D. Nat. P 符 Iˆ = {xα } |xα ⟩ ⟨xα | ， 如 此， 我 們 期 望 值 的 計 算 包 含 著 下 列 矩 陣 元 計 算：. n. a l X X Y Dx |Hˆ |x iEv Z= n Ch U engchi P 2n. α−1. tα ,aα. α. (2.47). {xα } S2n α=1. 我們也已知歸一化常數 Z 可表示為權重的相加 Z =. {x,S2n }. Wx,S2n ， 比 對. 式 (2.47) 可得權重 W ：. Wx,S2n =. 2n D Y. ˆ tα ,aα |xα xα−1 |H. E. (2.48). α=1. 如 此， 算 符 乘 積. Q2n α=1. ˆ tα ,aα 我 們 沿 著 算 符 乘 積 傳 遞 方 向 進 行 期 望 值 矩 陣 元 H. ⟨xα−1 |O|xα ⟩ 的測量並做平均。而哈密頓算符的矩陣元素乘積也影響算符更新的機 率，哈密頓算符的形式（式 (2.12)）與其矩陣元素的值（式 (2.17)、式 (2.18)）便 如同我們在 SSE 方法上的設計。 零溫投射法在算符更新的程序上與 SSE 方法大致相同，包含著針對對角算符 更新的局域更新以及針對非對角算符更新的叢集更新。不同的地方在於零溫投射 24. DOI:10.6814/NCCU201901066.
(35) 法並沒有像 SSE 方法引入 M − n 個單位算符，因此我們會在這兩大更新上作小修 改： 局域更新 針對一具有 N 個自旋、 Nb 個鍵結之量子易辛鐵磁模型，我們檢驗具有 2n 算符的 乘積 S2n 中每一個的算符，並以下列幾點規則作自旋狀態或算符序列的更新： ˆ −1,i ，則更新該算符所對應之晶格位 1. 若遇到一橫場項算符（非對角算符）H 置其作用後的自旋狀態，即自旋的翻轉。 ˆ 0,i 或 H ˆ 1,(ij) ， 我 們 先 拿 到 原 來 的 對 角 算 符， 並 2. 若 遇 到 一 對 角 算 符 H. 政 治 大. ˆ 0,i ； 或 以 式 (2.21) 的 機 率 置 入 一 鍵 以 式 (2.20) 的 機 率 置 入 常 數 算 符 H. 立. ˆ 1,(ij) ，這時置入的鍵結位置 (i, j) 相鄰的兩自旋必須為平行（相同 結算符 H. ‧ 國. 學. 方向）狀態，否則拒絕此鍵結算符的加入。. sit. y. Nat. 叢集更新. ‧. 上述的步驟必須重複進行至有一對角算符被置入為止。. al. er. io. 如同 SSE 方法的叢集更新法則，唯一的差別在於零溫投射法在模擬空間算符傳遞. n. 方向的兩端為開放邊界條件，因此這裡決定叢集的成長與終止有三法則：. Ch. engchi. i n U. v. ˆ −1,i 或常數算符 H ˆ 0,i 。 1. 叢集終止於單位算符，即橫場算符 H 2. 叢集終止於算符傳遞方向兩端邊界的自旋狀態。 ˆ 1,(ij) ，四端連結之自旋會被納入同一叢集，且會沿著四端 3. 若遇到鍵結算符 H 方向繼續成長。 當叢集建構完畢後，剩下就跟 SSE 方法的一樣，以 1/2 的機率翻轉叢集內所有的 自旋狀態與算符種類，其翻轉步驟也包含： • 自旋翻轉。 ˆ −1,i 變常數算符 H ˆ 0,i ；常數算符 H ˆ 0,i 變橫場算符 H ˆ −1,i 。 • 橫場算符 H 而無序系統的零溫投射法，其算符的更新法則，就如同第2.2節所述的方法一樣。 25. DOI:10.6814/NCCU201901066.
(36) 立. 政 治 大. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. 26. i n U. v. DOI:10.6814/NCCU201901066.
(37) Chapter 3 模擬結果. 治 政 我們測試蒙地卡羅退火法對探討以下幾類相變問題的可行性或效率： （一) 具 大 立 z = 1 的量子臨界點；（二）KosterlitzThouless (KT) 相態；（三）具無窮大動力學 ‧ 國. 學. 指數 z = ∞ 的無序臨界點。所應用的方法為第 2 章所描述的有限溫度隨機級數展 開量子蒙地卡羅方法（SSE）及零溫投射法。. ‧ y. sit. al. n. 平衡態模擬. er. io. 3.1.1. 自旋鐵磁鏈及方晶格鐵磁. Nat. 3.1. Ch. engchi. i n U. v. 為利用有限溫度 SSE 方法來探測絕對零溫量子相變，我們先針對固定系統尺 度 L 逐步降低溫度（提高 β）來檢驗期望值的收斂情形。若期望值於 β ≥ β ∗ 收斂 至一定值，我們便可將 β ∗ 及其以上的 β 值視為模擬上的等效零溫。這裡我們考慮 均質的量子易辛模型（式 (1.5)）定義於一維度及二維方晶格。 圖 3.1 與圖 3.2 分別展示對於 J = 1 的一維量子易辛鐵磁模型的 Binder 比值 g 、 磁化量平方 m2 在已知的量子臨界點 hc = 1 時對應不同的溫度 β −1 的關係，這裡 模擬的總蒙地卡羅步數為 200000 。從圖中可觀察到，對不同系統尺度 L 的期望值 數據的曲線在足夠低溫時 ( β ∗ ≈ 512 ) 會逐漸收斂至一定值。除此之外，我們將 β 視為虛數時間，利用 β ∼ Lz 的關係 ( z = 1 ) 與臨界定標關係 [10]，使對應不同的 L 的數據坐落於一曲線。這裡的使用臨界指數 ν = 1 以及 βm = 0.125 [33, 34]。. 27. DOI:10.6814/NCCU201901066.
(38) 0.6. 0.6. g. 0.8. g. 0.8. 0.4. 0.2. 0. 0.4. L = 16 L = 32 L = 64 L = 128 L = 256 1. 4. 16. 64. 256. 1024. z=1 0.2. 0 -3 10. 4096. β. -2. -1. 10. 0. 10. 10. β/L. (a). 1. z. 10. L = 16 L = 32 L = 64 L = 128 L = 256 2. 10. 3. 10. (b). Figure 3.1: (a)：一維量子易辛鐵磁模型在臨界點 hc = 1 時，不同溫度 β −1 下的 Binder 比值。當 β & 512 ， g 值逐漸趨於一定值。(b)：利用 β ∼ Lz 的關係及取動力學指數 z = 1，可以使對應不 同的 L 的數據坐落於一曲線。. 政 治 大. 立. 2 2βm /ν. io. 4. 16. 64. β (a). al. 256. n. 1. 1024. Ch. sit. 0.1. 0.5. L = 16 L = 32 L = 64 L = 128 L = 256. L = 16 L = 32 L = 64 L = 128 L = 256. y. Nat. 0.2. 0 -3 10. 4096. v ni. -2. 10. engchi U. er. 0.3. 1. mL. m. 2. 0.4. ν =1 z =1 βm = 0.125. ‧. ‧ 國. 0.5. 0. 1.5. 學. 0.6. -1. 0. 10. 10. β/L. 1. z. 10. 2. 10. 3. 10. (b). Figure 3.2: (a)：一維量子易辛鐵磁模型在 hc = 1 時，不同溫度 β −1 下的磁化量平方 m2 。當 β & 512 ， m2 值逐漸趨於一定值。(b)：利用 β ∼ Lz 的關係以及磁化量平方的有限尺度定標形 式，可以使對應不同的 L 的數據坐落於一曲線。這裡的使用臨界指數理論值 [33, 34]: ν = 1 以及 βm = 0.125。. 我們固定在足夠低 ( β ≥ β ∗ ) 的溫度 β = 1024 下，藉由改變橫場來偵測量子相 變點 hc 。圖 3.3(a) 表示著在 β = 1024 時 Binder 比值對橫場強度的關係，圖中可 觀察到對應不同的尺度大小 L 的曲線 ( L ≥ 32 ) 交會於 h = 1 。 Binder 比值與磁 化量平方的定標分析（圖 3.3(b) 與圖 3.4(b)）亦顯示量子相變點 hc = 1 。. 28. DOI:10.6814/NCCU201901066.
(39) 1. 1. L=8 L = 16 L = 32 L = 64 L = 128 L = 256. 0.9 0.8 0.8 0.6. g. g. 0.7. 0.4. 0.6. L=8 L = 16 L = 32 L = 64 L = 128 L = 256. 0.2 0 0.9. hc = 1 ν = 1. 0.5 0.4 1. 0.95. 0.3 -5. 1.1. 1.05. -4. -3. -2. h. -1. 0. (h-hc )L. (a). 1 1/ν. 2. 3. 4. 5. (b). Figure 3.3: (a)：一維量子易辛鐵磁模型在 β = 1024 時， Binder 比值對橫場強度作圖。當 hc = 1 時，對應不同的 L 的曲線會交會在 hc = 1 ，而 hc = 1 正是此模型零溫量子相變點。(b)：有限尺 度定標作圖。這裡的使用關聯長度臨界指數的理論值 ν = 1 。. 政 治 大. h (a). hc = 1 ν = 1. 0.6 0.4. al. 1.05. Ch. 1.1. -5. -4. engchi U. y. βm = 0.125. sit. 1. er. 2 2βm /ν. 1. 1.2. 0.8. n. 0.95. io. 0 0.9. L=8 L = 16 L = 32 L = 64 L = 128 L = 256. 1.4. Nat. 0.2. L=8 L = 16 L = 32 L = 64 L = 128 L = 256. 1.6. mL. m. 0.4. 1.8. ‧. 2. 0.6. 學. ‧ 國. 立. v ni. -3. -2. -1. 0. (h-hc )L. 1 1/ν. 2. 3. 4. 5. (b). Figure 3.4: (a)：一維量子易辛鐵磁模型在 β = 1024 時，磁化量平方對橫場強度作圖。(b)：有限尺 度定標作圖。這裡的使用的臨界指數為 ν = 1 以及 βm = 0.125 。. 為了後續關於二維量子易辛模型的討論，這裡我們也探測同樣在 J = 1 二維量 子易辛鐵磁模型對不同 β 的期望值收斂情形。固定橫場在量子臨界點 hc = 3.044 下 [35]，作 Binder 比值（圖 3.5(a)） 、磁化量平方（圖 3.6(a)）對不同 β 的關係。當 β ≥ 256 時，數據曲線逐漸收斂至一定值。這裡模擬的總蒙地卡羅步數為 200000 ，參考 [36] 所使用的臨界指數為 ν = 0.63 以及 βm = 0.327 ，卻無法使對應不同的 L 的數據完美坐落於一曲線。故這裡使用我自己嘗試出來的臨界指數為 ν = 0.637 以及 βm = 0.32 。擬合情況差的原因可能在於蒙地卡羅步數不夠多或者系統尺度 不夠大所造成的。 29. DOI:10.6814/NCCU201901066.
(40) 0.6 0.6 0.5 0.4. g. g. 0.4 0.3. L=8 L = 16 L = 32 L = 64. 0.2 0.1 0. 1. 2. 4. 8. 16. 32. 0.2. 0 -2 10. 64 128 256 512 1024. β. z=1. -1. 0. 10. 1. 10. β/L. (a). z. 10. L=8 L = 16 L = 32 L = 64 2. 10. (b). Figure 3.5: (a)：二維量子易辛鐵磁模型在 hc = 3.044 時，不同溫度 β −1 下的 Binder 比值。當 β & 256 ， g 值逐漸趨於一定值。(b)：利用 β ∼ Lz 的關係及 z = 1，可以使對應不同的 L 的數據 坐落於一曲線。. 政 治 大. 立. ‧ 國. 1.2 m. 64. β. 0.6. al. (a). 256. Ch. 0.4 -2 10. 1024. z =1 βm = 0.32 ν = 0.637. sit. 0.8. n. 16. 1. v ni -1. 10. engchi U. er. io 4. L=8 L = 16 L = 32 L = 64. y. Nat. 0.05. 1. m 2L 2β /ν. 0.1. 1.4. ‧. m2. 0.15. 0. 1.6. L=8 L = 16 L = 32 L = 64. 學. 0.2. 0. 1. 10. β/L. z. 10. 2. 10. (b). Figure 3.6: (a)：二維量子易辛鐵磁模型在 hc = 3.044 時，不同溫度 β −1 下的磁化量平方 m2 。當 β & 256 ， m2 值逐漸趨於一定值。(b)：利用 β ∼ Lz 的關係以及磁化量平方的有限尺度定標形 式，並使用 ν = 0.637 及 βm = 0.32，此與其他文獻 [36] 得出的 ν = 0.63 以及 βm = 0.327 有些微 差距。. 3.1.2. 模擬退火與 KibbleZurek 機制. 接著，我們要討論一種非平衡態的模擬過程，稱之為模擬退火（simulated annealing）[37]，是藉由將控制參數（例如：溫度、磁場。）從一起始值緩慢地調 降至零，使退火過程保持近似平衡態的最佳化方法。在計算過程中，調降控制參 數（統稱「退火」）的速率對於在通過臨界點時系統缺陷的產生有關鍵性的影響。 此類情形的發生與 Thomas W. B. Kibble 與 Wojciech H. Zurek 提出的 KibbleZurek 30. DOI:10.6814/NCCU201901066.
(41) 1 0.6 0.8 0.5. 0.4 0.2. L=8 L = 16 L = 32 L = 64 L = 128 L = 256. 0. 0.4. m. g. 2. 0.6. 0.3. L=8 L = 16 L = 32 L = 64 L = 128 L = 256. 0.2 0.1. -2. 0. 10. 0. 2. 10. -2. 10. u. 0. 10. 2. 10. 10. u. (a). (b). Figure 3.7: 對一維量子易辛鐵磁模型固定橫場 h = 1 並降溫度，作 (a)：Binder 比值對下降速率圖。 (b)：磁化量平方對下降速率圖。. 政 治 大 機制 ( KibbleZurek meshanism 立 ) 有關，其討論早期宇宙與凝態系統淬煉通過一臨. ‧ 國. 統 [40, 41]。. 學. 界點時缺陷與退火速率之關係 [38, 39]。KibbleZurek ( KZ ) 機制也擴展至量子系. ‧. 以二維易辛模型為例，臨界溫度 Tc 區分高溫時的無序態及低溫時的對稱破缺 有序態。在臨界點附近弛豫時間發散，即系統至平衡態所需的時間發散，此稱為. y. Nat. ′. sit. 臨界慢化現象 (critical slowing down)。弛豫時間 ξτ 發散情形可由動力學指數 z ′ 描. n. al. er. io. 述：ξτ ∼ ξ z ，此處先以 z ′ 標註可能與量子動力學指數 z （定義於式 (1.2)）不同的. i n U. v. 指數。KZ 理論的定量描述如下：弛豫時間 ξτ 隨著靠近臨界溫度的關係為. Ch. engchi ′. ξτ ∝ |T − Tc |−z ν. (3.1). 在退火過程中，當參數 T 以速率 u 且以線性方式隨時間 T 下降時： T (T) = T (0) − uT. (3.2). 根據 KZ 理論，這裡至臨界點之上剩餘時間 τ = |T − Tc | /u 須至少與弛豫時間 ξτ 相當，系統才得以保持近似平衡的退火過程。由此我們定義一特徵速率 uKZ ，稱 為 KibbleZurek ( KZ ) 速率，來界定退火速率慢 ( u . uKZ ) 或快 ( u > uKZ ) 。由 τ≡. ′ |T − Tc | ∼ ξτ ∼ |T − Tc |−z ν uKZ. 31. (3.3). DOI:10.6814/NCCU201901066.
(42) 我們得出 KZ 速率與臨界點距離 |T − Tc | 之關係 ′. uKZ ∼ |T − Tc |1+z ν. (3.4). 或與關聯長度 ξ 之關係 ′. uKZ ∼ ξ −(z +1/ν). (3.5). 上述的討論是基於熱力學極限下 L → ∞ 的情形下，因此對於有限大小系統，由 式 (3.5) 可得 KZ 速率與系統長度 L 的關係：. 政 治 大 ′. uKZ ∼ L−(z +1/ν). 立. (3.6). 學. ‧ 國. 所以非平衡態的物理量有限尺度定標形式 (式 (1.12) ) 可擴充表示為： . ˜ δL1/ν , uL−(z′ +1/ν) O (δ, L, u) = L−κ/ν O. . (3.7). ‧. y. Nat. 我們這裡對量子易辛模型進行 SSE 蒙地卡羅退火。在橫場固定在量子臨界值. sit. 下 hc /J = 1 ( 令 J = 1 )。我們藉緩慢線性降溫至等效零溫 β ∗ ，也就是逼近量子臨. n. al. er. io. 界點，來探討上述的 KZ 機制。若起始溫度為 Ti ，終止溫度為 T ∗ ， Smc = T 定義 為蒙地卡羅時間（步驟），這裡速率定義為. Ch. engchi. i n U. T ∗ −1 − Ti −1 u= Smc. v. (3.8). 我們設 Ti −1 = 0.2 ，而 T ∗ −1 = 1024 ，我們從 Smc = 21 至 Smc = 220 來調慢速率。 在系統大小有限的情況下，我們期待觀察量與速率的關係可由雷同於式 (3.7) 的定 標形式描述 . ˜ uL−(z′ +1/ν) O (L, u) = L−κ/ν O. . (3.9). ˜ 的第一個變 因為我們降溫至等效量子臨界點（T ∗ , hc ），因此式 (3.7) 定標函數 O 數消失。 我們將模擬得出的 Binder 比值以及磁化量平方與速率之關係呈現在圖 3.7。 依據式 (3.6)，我們作這兩個觀察量對速率之定標圖於圖 3.8 及圖 3.9，這裡我們 32. DOI:10.6814/NCCU201901066.
(43) 1.2. 1.2. 1. 1 0.8. g. g. 0.8 L=8 L = 16 L = 32 L = 64 L = 128 L = 256. 0.6 0.4 0.2. L=8 L = 16 L = 32 L = 64 L = 128 L = 256. 0.6. z′ = 1 ν =1. 0.4. 0. 2. 10. 4. 10. uL. 6. 10. 10. 0.2 -2 10. z′ = 0.3 ν =1. 0. 2. 10. z′+1/ν. uL. (a). 10 z′+1/ν. 4. 10. (b). Figure 3.8: Binder 比值模擬退火之定標圖，分別以不同的動力學指數 z 作定標，擬合的臨界指數為 ν = 1 。(a)：量子力學動力學指數 z ′ = 1 。(b)： SwendsenWang 演算法的臨界動力學指數 z ′ = 0.3 。. 政 治 大. 立. mL. 2 2βm /ν. ‧ 國. 1. 10. 0. 10. 0.5. 2. 10. (a). 4. al. 10. 6. 10. -2. 10. 0. Ch. 2. 10. v ni. z′+1/ν. n. uL. sit. io. -2. y. Nat. 0.5. er. 2 2βm /ν. 1. 1.5. βm = 0.125 z′ = 0.3 ν =1. ‧. mL. 1.5. L=8 L = 16 L = 32 L = 64 L = 128 L = 256. 2. βm = 0.125 z′ = 1 ν =1. 學. L=8 L = 16 L = 32 L = 64 L = 128 L = 256. 2. uL. 10 z′+1/ν. 4. 10. (b) U engchi Figure 3.9: 磁化量平方模擬退火之定標圖，分別以不同的動力學指數 z 作定標，擬合的臨界指數為. ν = 1 以及 βm = 0.125 。(a)：量子力學動力學指數 z ′ = 1 。(b)： SwendsenWang 演算法的臨界動 力學指數 z ′ = 0.3 。. 僅考慮緩慢速率部份。採用的動力學指數分別為量子動力學指數 z ′ = z = 1 和 SwendsenWang 演算法於二維古典易辛模型的動力學指數 z ′ = 0.3 [42]。從定標圖 來看， SwendsenWang 演算法的臨界動力學指數 z ′ = 0.3 的定標圖（圖3.9(b)）相 對於量子力學動力學指數 z ′ = 1 的圖（圖3.9(a)）擬合品質比較好，雖無法確定 z ′ = 0.3 即為 SSE 方法最適切的動力學指數。只能推斷 SSE 蒙地卡羅退火法呈現 的是演算法相關的動態行為，而非反應量子哈密頓動力學。也就是說，上節所討 論的以 z = 1 作隨 β 變化的平衡態觀察值之定標行為無法在此以退火法重現。. 33. DOI:10.6814/NCCU201901066.
(44) 3.1.3. 模擬退火. 雖然量子蒙地卡羅退火法不等效於可反應量子動力學的量子退火法（quantum annealing）。 但 考 慮 有 限 尺 度 的 系 統， 我 們 應 可 藉 蒙 地 卡 羅 退 火 法 極 度 緩 慢 改變參數的方式作近似平衡的計算。若下降速率慢於臨界速率 u < uKZ ，或 ˜ 的第 u ∼ L−α , α > z + 1/ν ，對應至非平衡的定標型式 (3.7) ，式 (3.7) 定標函數 O 二個變數將消失或趨近一常數。不同於前一節的討論，我們這裡在退火過程中時 時作觀察量的測量，若系統確能達近似平衡，得出的觀察量與參數的關係應該可 以在臨界點附近呈現正確的平衡態有限尺度定標形式（式 (1.12)）。 為討論量子臨界點的有限尺度定標，我們調降的參數為橫場（量子參數），但. 政 治 大. 亦同時調降溫度，我們將如此的退火路徑於圖 3.10 中的路徑 (i) 表示。另外，我. 立. 們也考慮固定低溫只調降橫場的退火路徑，以路徑 (ii) 表示。. ‧. ‧ 國. 學. paramagnet. sit. y. Nat. n. al. er. io. ferromagnet. Ch. engchi. i n U. v. Figure 3.10: 量子模擬退火的路徑示意圖。其中 hc 與 Tc 分別為量子與古典的臨界點，而紅色區域 為磁化量有序的鐵磁體，其餘區域為無序的順磁體。以一維模型而言，Tc = 0，即鐵磁體僅侷限 於 T = 0 軸上。我們循兩種路徑去逼近量子臨界點。. 路徑 (i) ：同時降溫降橫場 一維 首先我們對一維量子易辛鐵磁模型進行同時降溫及降橫場的退火模擬。退火過 程如下：從起始溫度 βi = 0.2 及起始橫場 hi = 5 隨蒙地卡羅步驟降至終止溫度 βf = 1024 及零場 hf = 0。若總蒙地卡羅步數為 Ssm ，每步調降的溫度及橫場距離 為： dβ =. βf − βi Smc 34. (3.10). DOI:10.6814/NCCU201901066.
(45) dh =. hi − hf Smc. (3.11). 圖 3.11(a) 與 圖 3.13(a) 分 別 用 兩 個 不 同 長 度 的 蒙 地 卡 羅 步 數 Smc = 214 及 Smc = 217 作 Binder 比值對橫場強度關係圖。我們觀察到雖然數值曲線呈現一些小 震盪，但足以看出它們與平衡態的計算結果大致吻合。尤其是更緩慢調降參數的 Ssm = 217 結果之有限尺度定標更佳擬合平衡態臨界相變的理論值；雖然在有限尺 度定標圖 3.13(b) 中我們拋棄較小系統尺度 ( L < 32 ) 的數值，但平衡態計算的結 果圖 3.3(b) 同樣尺度的數值亦呈現不甚理想的擬合品質（注意：那裡僅顯示很小 的參數範圍：0.95 ≤ h ≤ 1.05）。 1. ‧ 國. 14 L=8 Smc = 2 L = 16 L = 32 N = 1000 s L = 64 L = 128 L = 256 h=1. 0.6. hc = 1. 1.01 1.02 1.03 1.04 1.05. 0.4 -5. -4. -3. h. -2. -1. 0. (h-hc )L. n. al. er. io. (a). sit. 1. y. ν =1. Nat. 0.4 0.95 0.96 0.97 0.98 0.99. 14. Smc = 2. Ns = 1000. g. g. 0.8. L = 32 L = 64 L =128 L = 256. ‧. 0.6. 1. 學. 0.8. 立. 政 治 大. Ch. engchi U. v ni. 1 1/ν. 2. 3. 4. 5. (b). Figure 3.11: 一維量子易辛鐵磁模型線性降溫和降橫場的模擬退火，模擬參數為 Smc = 214 以及 Ns = 1000 。(a)：Binder 比值對橫場強度圖。(b)：Binder 比值的有限尺度定標圖。這裡採用的零 溫量子相變點為 hc = 1 ，臨界指數為 ν = 1 。. 35. DOI:10.6814/NCCU201901066.
(46) 0.8. 2. L = 32 L = 64 L = 128 L = 256. 1.8 1.6. 2. 2 2βm /ν. 0.6. 0.2. L=8 L = 16 L = 32 L = 64 L =128 L = 256. 1 0.8. 14. Smc= 2. 0.6. Ns = 1000. 0.4. 0 0.95 0.96 0.97 0.98 0.99. 14. Smc = 2 Ns = 1000 hc = 1 ν =1 βm = 0.125. 1.2. mL. m. 0.4. 1.4. 1. 0.2 -5. 1.01 1.02 1.03 1.04 1.05. -4. -3. -2. h. -1. 0. (h-hc )L. (a). 1 1/ν. 2. 3. 4. 5. (b). Figure 3.12: 一維量子易辛鐵磁模型線性降溫和降橫場的模擬退火，模擬參數為 Smc = 214 以及 Ns = 1000 。(a)：磁化量平方對橫場強度圖。(b)：磁化量平方的有限尺度定標圖。這裡採用的零 溫量子相變點為 hc = 1 ，臨界指數為 ν = 1 以及 βm = 0.125 。. io 1. h (a). al. n. 0.4 0.95 0.96 0.97 0.98 0.99. 1.01 1.02 1.03 1.04 1.05. Ch. 0.4 -5. engchi. y. 0.6. sit. L=8 17 L = 16 S = 2 mc L = 32 L = 64 L = 128 Ns = 1000 L = 256 h=1. 17. Smc = 2 Ns = 1000 hc = 1 ν =1. L = 32 L = 64 L = 128 L = 256. er. g. g. 0.8. Nat. 0.6. 1. ‧. 0.8. 學. 1. ‧ 國. 立. 政 治 大. -4. i n U -3. v. -2. -1. 0. (h-hc )L. 1 1/ν. 2. 3. 4. 5. (b). Figure 3.13: 一維量子易辛鐵磁模型線性降溫和降橫場的模擬退火，模擬參數為 Smc = 217 以及 Ns = 1000 。(a)：Binder 比值對橫場強度圖。(b)：Binder 比值的有限尺度定標圖。這裡使用的零 溫量子相變點為 hc = 1 ，臨界指數為 ν = 1 。. 36. DOI:10.6814/NCCU201901066.
(47) 0.8. 2. 1.6. 2. 2 2βm /ν. 0.6. L=8 L = 16 L = 32 L = 64 L = 128 L = 256. 17. 1.2. mL. m. 0.4. 0.2. L = 32 L = 64 L = 128 L = 256. 0.8 17. Smc = 2 Ns = 1000. 0 0.95 0.96 0.97 0.98 0.99. 1. Smc = 2 Ns = 1000 hc = 1 ν =1 βm = 0.125. 0.4 1.01 1.02 1.03 1.04 1.05. -4. -5. -3. -2. h. -1. 0. (h-hc )L. (a). 1 1/ν. 2. 3. 4. 5. (b). Figure 3.14: 一維量子易辛鐵磁模型線性降溫和降橫場的模擬退火，模擬參數為 Smc = 217 以及 Ns = 1000 。(a)：磁化量平方對橫場強度圖。(b)：磁化量平方的有限尺度定標圖。這裡使用的零 溫量子相變點為 hc = 1 ，臨界指數為 ν = 1 以及 βm = 0.125 。. 政 治 大. 立. ‧ 國. 學. 另外，我們也嘗試調降參數後再逆向調升參數，也就是說，到達 (βf , hf ) 時再 以相同的距離 dβ 、 dh 上升溫度及橫場值。因為平衡態觀察量之期望值應唯一，. ‧. 所以若調降參數過程的數據和調升參數過程的數據可重合，則表示退火過程（或 升火過程）系統已達近似平衡。針對前面討論的的例子，我們確實觀察到退火速. y. Nat. sit. 率過快的 Smc = 214 下降過程及上升過程大 L 的數據呈現叉開的行為（圖 3.15(a)） ，. n. al. er. io. 而達近似平衡態的 Smc = 217 數據則在於雙向過程重合一起（圖 3.15(b) ）；後者的. i n U. v. 雙向數據在臨界點附近的定標分析圖中亦呈現理想的擬合（圖 3.16 及圖 3.17） 。 1. Ch. engchi 1. 0.8. 0.8 14. g. g. Smc = 2. Ns = 400 0.6. 0.4 0.98. 0.6. L = 256 Up L = 256 Down L = 8 Up L = 8 Down 1. 0.99. 1.01. 1.02. 0.4 0.98. 17. Smc = 2 Ns = 400 L = 256 Up L = 256 Down L = 8 Up L = 8 Down 1. 0.99. h. 1.01. 1.02. h. (a). (b). Figure 3.15: 一維量子易辛鐵磁模型線性升降溫度和橫場的模擬退火及回火，取 Ns = 400 樣本，分 別顯示上升和下降過程的 Binder 比值的圖。(a)：採用蒙地卡羅步數 Smc = 214 。(b)：Smc = 217 。. 37. DOI:10.6814/NCCU201901066.
(48) 1. 1 17. Smc = 2 Ns = 800. L = 32 L = 64 L = 128 L = 256. 0.6. 0.8. 17. Smc = 2 Ns = 800 hc = 1 ν =1. g. g. 0.8. L=8 L = 16 L = 32 L = 64 L = 128 L = 256 h=1. 0.6. 0.4 0.95 0.96 0.97 0.98 0.99. 1. 1.01 1.02 1.03 1.04 1.05. 0.4 -5. -4. -3. -2. h. -1. 0. (h-hc )L. (a). 1 1/ν. 2. 3. 4. 5. (b). Figure 3.16: 一維量子易辛鐵磁模型線性升降溫度和橫場的模擬退火及回火，參數為 Smc = 217 以 及 Ns = 800。(a)：Binder 比值對橫場強度圖，虛線標示 hc = 1 。(b)：Binder 比值的有限尺度定 標圖。採用零溫量子相變點為 hc = 1 ，及臨界指數 ν = 1 。. 0.6. 1. h (a). 1.01 1.02 1.03 1.04 1.05. Ch. 0.2 -5. y. er. al. 17. Smc = 2 Ns = 800 hc = 1 ν =1 βm = 0.125. 0.4. n. 0 0.95 0.96 0.97 0.98 0.99. 1 0.8. io. L=8 L = 16 17 L = 32 L = 64 Smc = 2 L = 128 N = 800 s L = 256. 1.2. sit. 2 2βm /ν. Nat. 0.2. L = 32 L = 64 L = 128 L = 256. 1.6 1.4. mL. m. 0.4. 2 1.8. ‧. 2. 0.6. 學. 0.8. ‧ 國. 立. 政 治 大. i n U. -4. engchi. -3. v. -2. -1. 0. (h-hc )L. 1 1/ν. 2. 3. 4. 5. (b). Figure 3.17: 一維量子易辛鐵磁模型線性升降溫度和橫場的模擬退火及回火，參數為 Smc = 217 以 及 Ns = 800。(a)：磁化量平方對橫場強度圖。(b)：磁化量平方的有限尺度定標圖。採用零溫量子 相變點為 hc = 1 ，臨界指數為 ν = 1 以及 βm = 0.125 。. 以上我們所討論的均為等間隔 dβ、dh （線性）調整參數的結果。以下我們討 論以非等間隔的方式調整參數的情形。我們採用： dβ =. ln (βf ) − ln (βi ) Smc. (3.12). β = exp ((βi ) + dβ × i) , i = 0, 1, 2, · · · Smc 如此一來，在高溫高橫場時， β 下降（上升）速率較快；而在相對低溫低橫場時， 38. DOI:10.6814/NCCU201901066.
(49) β 下降（上升）速率較慢。比起線性方式下降（上升），非線性方式下降（上升） 能使我們在量子臨界點附近及有序態（h < hc ）有較長的模擬時間。 1 14. 0.9. 0.9. 0.8. 0.8. 0.7 0.6 0.5. L=8 L = 16 L = 32 L = 64 L = 128 L = 256 h=1. Smc = 2 Ns = 400 hc = 1 ν =1. 0.7. g. g. 1. 14. Smc = 2 Ns = 400. 0.6. L = 32 L = 64 L = 128 L = 256. 0.5. 0.4 0.95 0.96 0.97 0.98 0.99. 1. 0.4 -5. 1.01 1.02 1.03 1.04 1.05. h (a). 立. -4. -3. -2. -1. 0. (h-hc )L. 政 治 大. 1 1/ν. 2. 3. 4. 5. (b). 0.1. ‧ 國. 0 0.95 0.96 0.97 0.98 0.99. 1. Ch. 1 0.8. engchi 0.6. 1.01 1.02 1.03 1.04 1.05. h. 0.4 -5. y. sit. 2 2βm /ν. al. Smc = 2 Ns = 400 hc = 1 ν =1 βm = 0.125. er. L=8 L = 16 L = 32 S = 214 mc L = 64 L = 128 N = 400 s L = 256. 14. 1.2. mL. 2. 0.2. 1.4. n. 0.3. L = 32 L = 64 L = 128 L = 256. 1.6. io. m. 0.4. 1.8. Nat. 0.5. ‧. 0.6. 學. Figure 3.18: 一維量子易辛鐵磁模型非線性降溫度和線性降橫場的模擬退火，參數為 Smc = 214 以 及 Ns = 400。(a)：Binder 比值對橫場強度圖。(b)：Binder 比值的有限尺度定標圖。這裡使用的零 溫量子相變點為 hc = 1 ，臨界指數為 ν = 1 。. i n U. -4. -3. v. -2. -1. 0. (h-hc )L. (a). 1 1/ν. 2. 3. 4. 5. (b). Figure 3.19: 一維量子易辛鐵磁模型非等間隔降溫度和線性降橫場的模擬退火，參數為 Smc = 214 以及 Ns = 400。(a)：磁化量平方對橫場強度圖。(b)：磁化量平方的有限尺度定標圖。這裡使用的 零溫量子相變點為 hc = 1 ，臨界指數為 ν = 1 以及 βm = 0.125 。. 圖 3.20 顯示測試不同 Smc 對 L = 256 及 L = 8 在調降參數和調升參數過程的 Binder 比值的圖。其中 Smc = 218（圖 3.21）為足以達到近似平衡的蒙地卡羅步數， 其定標分析的擬合結果（圖 3.22 與 3.23）也顯示與平衡態結果相似的理想品質。 整體來看，非線性的降溫方法所得到的數據曲線較線性降溫的數據曲線光滑，我 們可推斷數據曲線的小震盪主要來自高溫。 39. DOI:10.6814/NCCU201901066.
(50) 1. 1 14. 15. Smc = 2 Ns = 400. Smc = 2 Ns = 400. g. 0.8. g. 0.8. 0.6. 0.4 0.98. 0.6. L = 256 Up L = 256 Down L = 8 Up L = 8 Down. 1. 0.99. 1.01. 0.4 0.98. 1.02. L = 256 Up L = 256 Down L = 8 Up L = 8 Down 1. 0.99. h (a) 1 16. 17. Smc = 2 Ns = 400. Smc = 2 Ns = 400. 政 治 大. 0.8. 0.8. g. g. 立. 0.6. 0.99. ‧ 國. L = 256 Up L = 256 Down L = 8 Up L = 8 Down. L = 256 Up L = 256 Down L = 8 Up L = 8 Down. 學. 0.4 0.98. 1.02. (b). 1. 0.6. 1.01. h. 1. 1.01. 0.4 0.98. 1.02. 1. 0.99. 1.01. 1.02. h. ‧. h. (c). (d). Nat. y. sit. n. al. er. io. Figure 3.20: 一維量子易辛鐵磁模型非等間隔升降溫度和線性橫場的模擬退火及回火，分別採用 Ns = 400 樣本，上升和下降過程的 Binder 比值。採用的蒙地卡羅步數分別為 (a)：Smc = 214 。 (b)： Smc = 215 。(c)： Smc = 216 。(d)： Smc = 217 。. 1. Ch. engchi. i n U. v. 18. Smc = 2 Ns = 400. g. 0.8. 0.6. 0.4 0.98. L = 256 Up L = 256 Down L = 8 Up L = 8 Down 1. 0.99. 1.01. 1.02. h Figure 3.21: 一維量子易辛鐵磁模型非等間隔升降溫度和等間隔升降橫場的模擬退火及回火，採用 Ns = 400 ， Smc = 218 ，在上升和下降過程的 Binder 比值。. 40. DOI:10.6814/NCCU201901066.
相關文件
In an Ising spin glass with a large number of spins the number of lowest-energy configurations (ground states) grows exponentially with increasing number of spins.. It is in
Abstract In this paper, we consider the smoothing Newton method for solving a type of absolute value equations associated with second order cone (SOCAVE for short), which.. 1
Conjugate gradient method
Step 4 If the current bfs is not optimal, then determine which nonbasic variable should become a basic variable and which basic variable should become a nonbasic variable to find a
Keywords: Fuzzy Delphi method (FDM); Analytic network process (ANP); Collective residence environment sustainable development; Housing community
Therefore, this study based on GIS analysis of road network, such as: Nearest Neighbor Method, Farthest Insertion Method, Sweep Algorithm, Simulated Annealing
【Keywords】Life-City; Multiple attribute decision analysis (MADA); Fuzzy Delphi method (FDM); Fuzzy extented analytic hierarchy process
Keywords: Ant Colony System (ACS), Farthest Insertion Method, Vehicle Routing Problem (VRP)... 第一章
