空間各向異性與無序性之 (3+1)維量子海森堡模型的蒙地卡羅研究

全文

(1)國立臺灣師範大學理學院物理學系博士論文 Department of Physics College of Science. National Taiwan Normal University Ph.D. Dissertation. 「空間各向異性與無序性之 (3+1) 維量子海森堡模型」的蒙地卡羅研究 Monte Carlo studies of (3+1)- dimensional quantum Heisenberg models with spatial anisotropy and disorder. 高銘佐 Ming-Tso Kao. 指導教授：江府峻博士 Advisor: Fu-Jiun Jiang, Ph.D.. 中華民國 103 年 7 月 July, 2014.

(2) ii.

(3) 誌謝首先，致上最深的感謝給我的指導教授江府峻博士，因為他的指導，鼓勵，耐心與支持，本論文才得以完成。感謝各位口試委員老師的建議與指教，本校的高賢忠教授，張明哲教授，政大林瑜琤教授，清大的陳柏中教授，中山大學陳宗緯教授。…也感謝我認識的朋友和我親愛的家人，漫長的研究生生涯，感謝你們的鼓勵和支持，才能完成博士學位。希望完成這個里程碑，對我個人，或是家庭，甚至是整個社會，都能有正面積極的貢獻和意義。. iii.

(4) iv.

(5) Acknowledgements I would like to express my deep gratitude to my supervisor , Prof. FuJiun Jiang, for his advice, encouragement, enduring patience and constant support. Without his consistent instruction, the thesis could not have reached its present form. I’m glad to thank my thesis committees, Prof. Hsien-Chung Kao, Prof. Ming-Che Chang, Prof. Yu-Cheng Lin, Prof. Po-Chung Chen, Prof. Tsung-Wei Chen and my supervisor, Prof. Fu–Jiun Jiang. They teach me a lots …. Additionally, I want to thank my friends and my family. Without your encouragement and support, I cannot finish my PhD degree.. v.

(6) vi.

(7) 摘要本論文主要是使用蒙地卡羅方法 (Monte Carlo Method) 來模擬研究 (3+1) 維量子海森堡模型 (quantum Heisenberg model)。特別是我們探討了空間各向異性 (spatial anisotropy) 與無序性 (disorder) 對此模型特性之影響。研究空間各向異性量子海森堡模型的動機是想要針對 dimerization 類別的海森堡模型，定量上去探討在量子臨界點附近 (quantum critical √ point) 新建立的普適關係 (universal relation)，即 TN / c3 ∝ Ms 。其中 TN 是 Néel temperature ，c 是自旋波速 (spin wave velocity) 及 Ms 是交錯磁化密度 (staggered magnetization density)。我們所作的模擬結果與 Sushkov [31] 藉由級數展開法 (series expansion) 所得到的結果是一致的。另外對無序性的研究，我們計算三維鍵結無序 (bond disorder) 量子海森堡模型的 TN 和 Ms ，方法是引進兩個參數，即隨機耦合強度 D 和隨機機率 P ，來描述反鐵磁交換耦合 (exchange couplings) Jij 的隨機性。D 和 P 的值皆在 0 和 1 之間，每個交換耦合強度為 Jij (1 + D) 或 Jij (1 − D) 的機率分別為 P 及 (1 − P ) 。我們發現對這種無序性模型在靠近乾淨系統附近，用平均交換耦合強度 J 歸一化的 TN (即 TN /J) 和交錯磁化密度 Ms 之間也呈現一種線性關係。. 關鍵字：量子海森堡模型，空間各向異性，無序性，蒙地卡羅模擬. vii.

(8) viii.

(9) Abstract In this thesis we simulate the three-dimensional quantum Heisenberg model using first principles Monte Carlo method. We focus on the effects of spatially anisotropy and random-exchange disorder. Our motivation is to investigate √ quantitatively the newly established universal relation TN / c3 ∝ Ms near the quantum critical point (QCP) associated with dimerization. Here, TN , c and Ms are the Néel temperature, the spin wave velocity and the staggered magnetization density, respectively. Our Monte Carlo results agree nicely with the corresponding results determined by the series expansion method. As for the random-exchange disorder, the randomness for the antiferromagnetic exchange couplings Jij (bond disorder) for any two nearest neighbour spin ⟨ij⟩ is introduced by two parameters D and P . Specifically, given a set of 0 < P < 1 and 0 < D < 1, the probability that each antiferromagnetic coupling takes the value Jij (1 + D) (Jij (1 − D)) is P (1 − P ). Remarkably, in contrast to the scenario of the dimerized systems that the linear relation between TN and Ms appears close to a quantum critical point at which the antiferromagnetism is destroyed, for the model considered here the Néel temperatures, when being normalized properly, scale linearly with the staggered magnetization density near the data associated with the clean system. Our study also confirms that in three dimensions the antiferromagnetism is robust against the employed bond disorder. Keywords:Quantum Heisenberg model, spatial anisotropy, disorder, Monte Carlo simulations. ix.

(10) x.

(11) Contents. 誌謝. iii. Acknowledgements. v. 摘要. vii. Abstract. ix. 1 Introduction 導論. 1. 2 Motivation 研究動機. 7. 3 Methods 研究方法 3.1. 3.2. 11. 蒙地卡羅方法 (Monte Carlo Method) . . . . . . . . . . . . . . . . . . .. 11. 3.1.1. 細致平衡與遍歷性 (Detailed balance and Ergodicity) . . . . . .. 12. 3.1.2. Metropolis 演算法 . . . . . . . . . . . . . . . . . . . . . . . . .. 14. 3.1.3. 隨機級數展開法 SSE . . . . . . . . . . . . . . . . . . . . . . .. 14. 3.1.4. Measurement in SSE . . . . . . . . . . . . . . . . . . . . . . . .. 23. 誤差分析 Error analysis . . . . . . . . . . . . . . . . . . . . . . . . . .. 23. 3.2.1. 數據分級 binning . . . . . . . . . . . . . . . . . . . . . . . . .. 24. 3.2.2. Jackknife 分析 . . . . . . . . . . . . . . . . . . . . . . . . . . .. 24. 3.2.3. Bootstrap 分析 . . . . . . . . . . . . . . . . . . . . . . . . . . .. 25. 4 Numerical results 數值結果 4.1. 27. 3d ladder 空間各向異性之海森堡模型 . . . . . . . . . . . . . . . . . . xi. 27.

(12) 4.2. 4.1.1. Néel 溫度 TN 的計算 . . . . . . . . . . . . . . . . . . . . . . .. 28. 4.1.2. 交錯磁化強度密度 Ms 的計算 . . . . . . . . . . . . . . . . . .. 30. 4.1.3. 自旋波速 c 的計算 . . . . . . . . . . . . . . . . . . . . . . . . .. 30. 4.1.4. 理論預測與蒙地卡羅的比較 . . . . . . . . . . . . . . . . . . .. 32. disorder 無序性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 35. 4.2.1. 不同 (D，P ) 所對應之 Ms 和 ρs 的計算 . . . . . . . . . . . . .. 38. 4.2.2. 不同 (D，P ) 所對應之 TN 的計算 . . . . . . . . . . . . . . . .. 39. 4.2.3. TN 和 MS 之線性關係的驗證 . . . . . . . . . . . . . . . . . . .. 41. 5 Discussion and Conclusions 討論與結論. 45. Bibliography. 47. xii.

(13) List of Figures 1.1. 釔鋇銅氧，本圖摘自 http://upload.wikimedia.org/wikipedia/commons/0/06/Ybco.jpg . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 1.2. 高溫超導摻滲電子電洞相圖，本圖摘自 http://for538.wmi.badwmuenchen.de//projects/P4_crystal_growth/. 1.3. 1. . . . . . . . . . . . . . . . .. 2. 鉈銅氯 Thallium copper chloride (T lCuCl3 ) 隨著 g 的變化而可能所在的態，本圖摘自參考文獻 [29] . . . . . . . . . . . . . . . . . . . . . .. 5. 2.1. 3D dimerized model ，粗的 bond 表示較強的 coupling . . . . . . . . .. 7. 2.2. T lCuCl3 實驗與模擬的比較，摘自參考文獻 [31] . . . . . . . . . . .. 8. 3.1. 四種在 bonds 的作用下可能的 spin 狀態 . . . . . . . . . . . . . . . . .. 18. 3.2. propagated state . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 19. 3.3. loop update 改變 bond 的種類及翻轉其所對應的 spin . . . . . . . . . .. 21. 3.4. loop update 跨過邊界的情形 . . . . . . . . . . . . . . . . . . . . . . . .. 21. 3.5. single spin flip update . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 22. 4.1. ( 3 + 1 ) 維空間各向異性的 ladder-dimer 量子海森堡模型 . . . . . . .. 27. 4.2. 蒙地卡羅在 J ′ /J = 3.5 時所得到的 ρs L (上圖) 及 ρs3 L (下圖) 之數據，虛線僅是引導作用，盒子尺寸為 L = 8, 12, 16, · · · , 32, 36 和 40，其中，ρs L = (ρs1 + ρs2 )L/2。 . . . . . . . . . . . . . . . . . . . . . .. 4.3. 29. 由蒙地卡羅計算所得出的 ⟨(m2s )⟩ 對 L 的關係圖，從上而下對應的 J ′ /J 依序是 2.5, 3.0, 3.125, 3.25,· · · ,3.75 和 3.875，實線僅是引導作用。 31. 4.4. 從蒙地卡羅計算所得到的 Ms 的值。實線則是從文獻 [34] 重製而得到的。 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii. 31.

(14) 4.5. 圖上顯示 TN 對 J ′ /J 藉由蒙地卡羅 (實心圓點)，級數展開 (實心方塊)，及場論計算 (實線) 所得到的關係圖。場論及數值展開的結果是. 4.6. 從文獻 [31] 估計和重製，dashed 及 dotted 線僅是引導作用。 . . . . . √ 圖上顯示 TN / c3 對 Ms 的關係圖。實線是用 a + bMs 形式擬合. 33. 的結果，虛線是用 b1 Ms + d1 Ms log(Ms ) 的擬合結果，正方形點是 J ′ /J = 1.0 的結果，J ′ /J = 1.0 的結果是參考文獻 [51] 和 [52] 的計算得來的。 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 34. 4.7. β- doubling scheme 的示意圖，摘自文獻 [42] . . . . . . . . . . . . . .. 36. 4.8. 對於不同 L ，D 和 P ，S(π, π) 收斂到其基態值的示意圖 . . . . . .. 36. 4.9. 不同的 D 和 P ，其平均交錯結構因子 S(π, π) 對 L 的關係圖，虛線僅是引導作用 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 37. 4.10 不同的 D 值，其 Ms 對 P 的關係圖。在固定 D 之下，可看出 Ms 在 P ∼ 2 附近有最小值，虛線僅是引導作用 . . . . . . . . . . . . . . . .. 37. 4.11 不同的 D 值，其 ρs 對 P 的關係圖，虛線僅是引導作用 . . . . . . .. 38. 4.12 不同的 D 值，其 ρs /J 對 P 的關係圖，虛線僅是引導作用 . . . . . .. 39. 4.13 ρs L for L = 16, 20, 24, · · · , 36, 40 ，上圖為 (D, P )=(0.70, 0.20)，下圖為 (D, P )=(0.50, 0.30)，虛線僅是引導作用 . . . . . . . . . . . . . . .. 40. 4.14 TN /J (上圖) 及 TN (下圖) 對 Ms 之關係圖，虛線僅是引導作用 . . .. 42. 4.15 (3+1) 維 ladder-dimer 量子海森堡模型之 TN /J 對 Ms 之關係圖，虛線僅是引導作用 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 43. 4.16 (3+1) 維 ladder-dimer 量子海森堡模型及 P = 0.3 之 TN /J 對 Ms 之關係圖，虛線及實線僅是引導作用 . . . . . . . . . . . . . . . . . . .. xiv. 43.

(15) List of Tables 4.1. 利用 finite-size scaling 所得到的 TN 的值。所有擬合的 χ2 /DOF 皆小於 1.6 。 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 4.2. 30. 對各種耦合 J ′ /J ，藉由 winding number squared 所求出的自旋波速 c 的數值。 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. xv. 32.

(16) xvi.

(17) Chapter 1 Introduction 導論本章我們將介紹在此論文中所研究的模型 –海森堡模型 (Heisenberg model)，特別是我們會著重於說明此模型的重要性及它與實際材料間的關聯性 [1]，同時我們也會介紹與此論文有關的觀念，如在二階相變臨界點附近的物理現象 [2]。. 本圖摘自 http://upload.wikimedia.org/wikipedia/commons/0/06/Ybco.jpg Figure 1.1: 釔鋇銅氧，一般熟知的陶瓷材料是反鐵磁材料，它在摻滲電洞或電子至某一濃度 (稱為臨界密度) 且降溫至某一低溫以下 (稱為臨界溫度 Tc )，其性質會從反鐵磁轉變為超導 1.

(18) [1]。在圖1.1中，可以清楚看到 Y Ba2 Cu3 O7 此一高溫超導體的幾何結構。雖然從圖上看出它是三維的晶格結構，但實際上，從實驗上的數據顯示層與層之間的交互作用很弱，而主要強的交互作用是落在銅氧 (CuO) 平面上 [3] [4]。因為銅原子的位能很大，電子在層與層之間發生穿隧效應的機率非常小。所以這材料的物理基本上是由銅氧平面所決定，也就是說雖然它們有三維晶格結構，而實際上的物理卻是二維的。因此在文獻上幾乎所有與高溫超導有關的研究，都集中於探討二維的模型 [1] [5]−[12]。. Figure 1.2: 高溫超導摻滲電子電洞相圖，本圖摘自 http://for538.wmi.badwmuenchen.de//projects/P4_crystal_growth/ 經過多年的研究，雖然對於高溫超導的機制並沒有完全弄清楚，但一般來說，圖1.2為人們所相信的高溫超導對於摻滲 (doping) 的相圖。對於原來在室溫下是反鐵磁的材料，在適當的條件，摻入電子或電洞之後，隨著摻入電子或電洞的濃度增加，其反鐵磁長程有序會被破壞掉，而濃度到某個臨界點時，這些材料會發生相變而成為高溫超導體。特別是在這些高溫超導體內並無 long-range 的 order ，而只有 strong short range (及 intermediate range) 的 correlations [1]。當然 doping 的濃度不同，會造成許多不同的態 (phase)，這正是此類型高溫超導材料豐富多樣的特性，因此吸引許多科學家去研究，這也是我們對反鐵磁材料有濃厚研究興趣的理由之一。另外當然也因為它與高溫超導有密切的關聯，希望有機會也可以進一步明白高溫超導的機制。而在與高溫超導相關的模型之中，最基本也是最重要的是所謂的 Hubbard model 2.

(19) (如果只考慮 hole doped 的情形下則為 t-J model)。它的 Hamiltonian 的形式是:. H = −t. ∑. (ˆ a†iσ a ˆjσ. +. a ˆ†jσ a ˆiσ ). +U. ⟨ij⟩σ. N ∑. (ni↑ ni↓ ). i=1. 其中，ˆ a† 和 a ˆ 是 fermion 的 creation operator 和 annihilation operator ，n 是 fermion number operator ，t 是 hopping term ，U 是 on site repulsion potential ，⟨ij⟩ 表示最近鄰，σ 表示 spin 的正或負 (上或下)，N 是 lattice sites 的總數。由於在 Hubbard model 中，ˆ a† 及 a ˆ 都是 fermion operators ，他們會滿足 {ˆ a†iσ , a ˆjσ′ } = δij δσσ′ 及 {ˆ a†iσ , a ˆ†jσ′ } = {ˆ aiσ , a ˆjσ′ } = 0 ，所以基本上會有 sign problem ，因此蒙地卡羅方法不適合用來研究這個模型。這裡 Hubbard model 是比較 general 的理論，對週期性的 potential 有很好的近似，可以用來解釋導電的機制，如導體、絕緣體等。另外它的 Hamiltonian 有 hole-electron symmetry ，當 Hubbard model project 到 hole ，我們就得到 t-J model (t-J model 的 Hamiltonian 的形式和上式的差別，僅在第二項) [1]。前面提到高溫超導可以藉由反鐵磁材料摻入電子或電洞而得到，而對於反鐵磁材料，其相關模型是 (量子) Heisenberg model 。它的 Hamiltonian 如下：. H=J. ∑. S⃗i · S⃗j. ij. 其中 J 是 anitferromagnetic coupling ，S⃗i 是 spin-half operator。從 Hubbard model 出發，在 half-filling 及 U ≫ t 的條件下，用微擾展開可以證明 Hubbard model 會和 Heisenberg model 一致 [13] [14] [15]。比較 Hubbard model 和 Heisenberg model 之間的差別，Hubbard model 有 hopping term ，但有 Pauli exclusion principle 的限制，所以晶格上的 site 可以是空的或者是有一對 spin 上下的電子，但是 Heisenberg model 要求晶格上的 site 只能有一個電子，spin 上或 spin 下。因為 Heisenberg model 的形式很簡單，所以為了解釋新的實驗結果，物理學家會從和 Heisenberg model 相關的模型出發。比如說，針對不同的初始條件，可以加上磁場，或是改變 coupling 的強弱，或是排列的形式等等 [16] [17]。這也是為什麼我們對 Heisenberg model 很有興趣。另外對於 spin-1/2 的 Heisenberg model 在 non-frustrated lattices 上沒有 sign problem ，雖然它的 Hamiltonian 仍都是 fermion operators ，但可證明其所產生的負號會相消，所以我們可以用蒙地卡羅方法來研究它 [18] [19]。 3.

(20) 總結以上論述，為什麼在過去 20 多年間 Heisenberg-type models 無論在數值上或是解析上都被廣泛的研究，有兩個簡單的理由，其一是在正方晶格 (square lattice) 上的量子 (spin-1/2) 海森堡模型可定性的，甚至定量的描述了未摻雜雜質的銅氧化物 (undoped cuprate)，而這些 undoped cuprate 材料在摻雜電荷載子之後有可能變成高溫超導體，這正是物理學家們至今對海森堡模型仍有著濃厚的研究興趣的一個最主要原因。此外，可應用於此類模型的數種高效率蒙地卡羅演算法的提出 [18] [20] [21] 及高速的中央處理器的進展，使得我們對「非自旋阻挫態」(nonfrustrated) 的海森堡模型可以有非常精確的數值結果 [1]，因此這類模型非常適合用來測試新的想法或者驗證理論預測的結果，而這也是 Heisenberg-type models 仍廣受歡迎的另一個因素。至於蒙地卡羅方法是常見的數值模擬方法，簡單說就是使用隨機亂數抽樣去作統計平均後得到所要計算的物理量，其內容會在第三章介紹。另外在本論文所提到的蒙地卡羅方法主要是參考 Anders Sandvik 教授所提出的隨機級數展開法 (Stochastic series expansion (SSE) method) 的演算法 [22] [23] [24] [25] ，也會在3.1.3 節中提及其主要的內容。在模擬計算海森堡模型時，需要給定晶格的大小，通常我們也會使用某種邊界條件，如週期邊界條件，來增加統計的多變性。但在模擬計算中，所選取晶格尺寸大小會影響所算出來的值，而在有限制大小的模擬數據中，我們可找出其相應的 bulk 性質，這與所謂的「有限尺度理論」(Finite-size scaling theory) [2] [26] [27] [28] 是有相關的。有限尺度理論是研究臨界點附近的有限尺度效應。在熱力學系統中，可使用序參量 (order parameter) 及幾個臨界指數 (critical exponents) 來描述系統的臨界現象。其中，關聯長度 (correlation length) 是指，熱力學系統中某處的原子或分子的狀態如果發生變化，距離多遠的原子或分子會受到影響。用二維的 Ising model 為例子，在臨界溫度 Tc 以下，自旋都會是同一個值 (+1 或是 −1)，在 Tc 以上，自旋會以數個不同的小區域為單位，有的是正號或是負號，當 T → ∞，可以說系統的關聯長度趨近於零，自旋的正負號完全是隨機的，反之在 Tc 時，關聯長度則會趨近於無限大。因此模擬計算中可達到的尺寸大小有時可能會影響所得到的結論。而我們也會在稍後的論文探討這一點。量子相變是指在絕對零度時，不同量子態之間的一個相變，特別是這個相變 4.

(21) 是由量子效應所控制的。Subir Sachdev 教授特別從絕對零度出發，引進可調參數 (tuning parameter) g-factor 去描述量子的臨界現象 [29] [30]。在絕對零度，沒有熱擾動 (thermal fluctuation)，但由於海森堡測不准原理 (Heisenberg’s uncertainty principle) 造成量子擾動 (quantum fluctuation) 使得電子內的自旋組態 (spin configuration) 有不同的變化而有不同的相態 (phase)。因此，在量子臨界點，相變似乎僅發生在絕對零度，然而實際上即使在非零溫下，仍有量子擾動的影響。如超導、鐵磁 (ferromagnetic) 和反鐵磁 (antiferromagnetic) 等不同的量子相態 (quantum phase)，也存在於極低溫之下。. Figure 1.3: 鉈銅氯 Thallium copper chloride (T lCuCl3 ) 隨著 g 的變化而可能所在的態，本圖摘自參考文獻 [29] Sachdev 特別提到，強關聯電子系統可以有許多不同的量子相態，並且可以藉由改變外在的參數 (external parameter) 來改變這些相態。這些相態及外在參數可以是，固體材料所受到的壓力，外加的磁場強度，或藉由摻雜改變電子的密度等等。Sachdev 提出 Thallium copper chloride (T lCuCl3 ) (鉈銅氯) 這個例子來說明可調參數 g-factor 這個概念。在 Thallium copper chloride 中，Cu2+ 離子會是在 spin-1/2 態下，在一般的大氣壓下，自旋喜歡配對成 dimers 且形成單重態 (singlet) 鍵結 √ (| ↑↓⟩ − | ↓↑⟩) / 2。在圖1.3中，紅色實線表示較強的耦合 J ，且 J > 0 ，綠色虛線是較弱的耦合 J/g，且 g > 1 ，T lCuCl3 的基態會隨著不同大小的 g 值的變化，而有不同的行為。對 T lCuCl3 而言，g 與壓力成反比，當壓力增加，g 變小，T lCuCl3 會經歷量子相變成為 Néel 反鐵磁態 [29]。而在 dimerized Heisenberg 模型中模擬的結果，在量子臨界點 g = gc 附近會有 Néel 反鐵磁有序態的相變點 [31]，這也是為什麼人們想研究 dimerized Heisenberg model 在臨界點附近的行為，藉此而尋找是否有普適關係存在。 5.

(22) 6.

(23) Chapter 2 Motivation 研究動機以往由於反鐵磁二維的 (量子) Heisenberg model 與高溫超導模型有密切的關係 [1]，而高溫超導主要是由 CuO 平面上的交互作用決定它的物理特性，所以主要的研究都集中在二維的海森堡模型，相對而言三維的海森堡模型比較少人研究。即使如此，自然界也存在和三維 Heisenberg model 相關的材料，如在上一章提過的鉈銅氯 T lCuCl3 即是一個例子。這個材料的實驗數據，顯示在圖2.2，特別是它會因為壓力的變化而有不同的 phase ，當壓力增加到大約在 1.07 kbar，其反鐵磁的 order 會出現 [29] [32] [33]。. Figure 2.1: 3D dimerized model ，粗的 bond 表示較強的 coupling Susukov [31] 這一組研究學者，建議可以用三維 dimerized 反鐵磁的 Heisenberg model 來解釋實驗的結果。這三維的 dimerized model 在 z 方向是較強與較弱的 couplings 交錯出現，在其他兩個方向是較弱的 couplings ，如圖 2.1 所示。在他們 7.

(24) 的計算中，當強弱 couplings 的 ratio 大約在 4.013 附近，會對應到 T lCuCl3 的臨界壓力 1.07 kbar。. Figure 2.2: T lCuCl3 實驗與模擬的比較，摘自參考文獻 [31]. 在圖2.2 顯示 Susukov 作的數值計算和實驗室作的實驗結果有相當的吻合度。圖2.2 中的點是實驗的結果，而實線及虛線為 [31] 所得到的結果。由圖中我們可看出 [31] 所得到的 spin gap ∆x ，∆y 和 ∆z 在靠近量子臨界點 (QCP) 附近與實驗結果是一致的。他們也計算了 Néel 溫度 TN ，雖然數值與實驗有些差距，但兩者之間的趨勢還是符合。在這裡為什麼要特別考慮 TN 呢？這是因為 [31] 中預測了一個與 TN 有關的 relation ，其關係式如下： T √N = c3. √. 12 Ms 5. (2.1). 這個關係式是由場論計算得到的，上式中，c 是 spin wave velocity ，Ms 是 staggered magnetization density 。這裡 universal 指的是與 microscopic 的細節無關，任何適用的 (dimerized) models ，其 TN ，c ，Ms 都會滿足這個關係。這個關係式從實驗的角度來說也是重要的，比如說從兩個已知量可利用上式求得第三個物理量的值。所以我們的動機就是想要檢驗這個關係式，特別是用蒙地卡羅的方法來檢驗。除了前面提的關係式之外，還有以下兩個關係式，是由 Jin 和 Sandvik 所提出 8.

(25) 的 [34]， TN = AMs , J¯ TN = A′ Ms T⋆. (2.2) (2.3). 上式中，J¯ 是 antiferromagnetic couplings 的平均強度，而 T ⋆ 是測量 uniform susceptibility 最大值所對應的溫度。嚴格來說，Jin 和 Sandvik 所發現的是更廣義的 universal relations 。特別是在靠近 critical points 附近，這些更廣義的 universal relations 會簡化為 (2.2) 及 (2.3) 式。這些關係式的重要性在於，若是實驗上，J¯ 不太容易測量，而我們知道 TN 和 Ms 的實驗數據，就可以透過這個線性關係式得到 J¯ 的數值，也就是說對實驗而言，這些關係式有它們的重要性。目前為止，(2.2) 及 (2.3) 式對 dimerized 引起的相變在可適用的範圍內，都被驗證是對的 [34]。因此，我們所要作的研究就是使用蒙地卡羅模擬的方法驗證，在靠近 QCP 附近，(2.1) 式這線性關係對 3D 的 dimerized Heisenberg model 是否正確 [35]。除此之外，並非材料都會像 T lCuCl3 一樣，可以用規則排列的 dimerized Heisenberg model 來描述。為了模擬比較貼近真實的系統，一種方法是採用 antiferromagnetic coupling disorder (或稱 bond disorder ) 作適當的安排。在實際上的材料內不同晶格對的 couplings (bonds) 的強度可能會在某一平均值附近變化。所以我們可以引進 bond disorder ，即根據某種機率分佈，對不同的 Jij 選擇它的強度。理論上，有人主張對二維以上 bond disorder 不會破壞反鐵磁長程有序。這個預測在二維的情況是被確定的，詳細來說在二維的情況，box-like 及 singular distribution 的 bond disorder 都不會破壞長程有序 [30] [36] [37] [38]。所以我們想要了解這些線性關係，如 TN /J¯ = AMs 在有 bond disorder 的情況，是不是也成立。另外對於三維反鐵磁 Heisenberg model 的情況下，我們也想明白 bond disorder 是否會破壞反鐵磁長程有序 [39]。我們會在第 4 章詳細的呈現數值模擬的結果。. 9.

(26) 10.

(27) Chapter 3 Methods 研究方法這章主要回顧在使用蒙地卡羅模擬計算的方法來研究海森堡模型時所需要的背景知識，包括蒙地卡羅模擬，我們所用到的演算法，以及統計誤差分析。特別是本章的題材主要是參照參考文獻 [18] [19] [20] [28]。通常蒙地卡羅方法可粗略分成兩類，一種是所要求解的問題具有內在隨機性，所以可藉由電腦可重複快速的運算能力來直接模擬這種過程，如在核物理中，研究中子在反應爐中的傳輸過程，即是此類問題的一個代表例子。另一種是所考慮的問題可轉化成某種隨機分布的函數，透過隨機取樣及隨機事件出現的機率等概念來求解。這裡所探討的問題是偏向後者。. 3.1 蒙地卡羅方法 (Monte Carlo Method) 統計物理可以用來處理多粒子系統，特別是對系統內的微觀狀態，我們可利用配分函數 (partition function) 的觀念來計算物理量的平均值，例如對物理量 A 其統計的結果是. ∑ −Ei /kT e A ⟨A⟩ = ∑i −Ei /kT . ie. 其中 i 是指系統所有可能的態，e−Ei /kT 是其所對應的波茲曼權重 (Boltzmann weight) ，k 是 Boltzmann constant ，而 Ei 是系統在 i 此一微觀態的能量。蒙地卡羅方法是最重要的一種可用來解決統計物理問題的數值方法 [19] [20]。用電腦來描述真實系統時，不可能去計算粒子間所有的交互作用，蒙地卡羅模擬的精神就是經過隨機亂數適當的取樣，讓計算出來的答案，符合統計力學的要求， 11.

(28) 可以較快速且較準確的逼近真實的結果。在數學上，產生亂數，是指從一給定的數集合中，不按順序隨機選取其中的數字。例如擲骰子，假設是公平的骰子，那麼六面中，擲出任一面的機率都是均等，也就是被選到的機率都是相同的 1/6 。這樣的情形，我們稱為骰子的每一面的權重 (weight) 都相等。重複擲骰子數次，一一記錄，其結果應該是均勻亂數的數列。在模擬計算中一般常利用電腦產生均勻分布在 [0,1] 之間的 (偽) 亂數。在過去常見的亂數產生器 (random number generator) 是用賭場輪盤以機械的方式產生的，這正是此數值方法以摩洛哥大賭城 -蒙地卡羅為名的含意。蒙地卡羅方法是基於大數法則，意思是當丟骰子次數越多，其平均值就越趨近於理論值，因此需要大量重複的抽樣，隨現代電腦軟硬體的進步，及演算法來處理取樣技巧的發展，提供了蒙地卡羅模擬發揮的舞台。而所謂的演算法是指，在有限的步驟解決數學問題的程序。蒙地卡羅的演算法有許多種，而其中最重要及最基礎的是 Metropolis 演算法。蒙地卡羅可以用來解決統計力學的問題，在一個統計力學適用的系統中，其物理觀測量的期望值 ⟨f ⟩β ，系統所可能的態 σ 的機率 Pβ (σ) ，及配分函數 Zβ 有以下關係：. Pβ (σ) =. e−βH(σ) Zβ. H(σ) = ⟨σ|H|σ⟩ ∑ Zβ = e−βH(σ) ⟨f ⟩β =. ∑. σ. f (σ)Pβ (σ). σ. 其中 β 是溫度倒數，H(σ) 為其 Hamiltonian 在態 σ 的期望值。我們可以透過正確的 Boltzmann weight e−βH(σ) 來產生 (generate) 正確的態 (configurations) ，其細節將於下一小節介紹。. 3.1.1. 細致平衡與遍歷性 (Detailed balance and Ergodicity). 蒙地卡羅在實際上操作，是利用偽隨機亂數對微觀系統作取樣，以趨近達到平衡態的真實系統，所以需要設計適當的演算法。基本上這演算法需要滿足的條件是，細緻平衡 (detailed balance) 與各態歷經 (ergodicity) [20]，在介紹這兩個概念之前， 12.

(29) 我們必須先提一下馬可夫鍊 (Markov chain)，這個概念因俄羅斯數學家安德烈 ·馬可夫 (Andrey Markov) 而得名。馬可夫鏈是指一個時間離散的數學系統，描述從某一狀態跳到另一狀態的過程。它是沒有記憶的隨機過程。在馬可夫鏈的每一步，系統根據機率分布，可以從一個狀態變到另一個狀態，也可以保持當前狀態。狀態的改變叫做 transition ，與不同的狀態改變相關的機率叫做 transition probability 。若用一個序列 {c1 , c2 , c3 , · · · } 來描述 Markov chain 的過程 [19]，則 Markov chain 可表示為 c1 → c2 → c3 · · · 其中 c1 有一定的機率可以到不同的態 c2 ，或是維持在 c1 ，也就是說 c2 可以是 c1 。此過程中，在給定當前資訊的情況下，只有 ci 狀態可用來預測 ci+1 ，過去的 ci−1 狀態對於預測未來 ci+1 狀態是無關的。例如隨機漫步就是馬可夫鏈的例子。隨機漫步中每一步的狀態是在二維平面上的點，每一步可以移動到任何一個相鄰的點，其移動機率可假設都是相同的 (無論之前漫步路徑是如何的)。這樣可以預期當系統平衡時，隨機漫步應該要滿足各態歷經和細緻平衡。. 各態歷經說明系統中對任兩個可能的狀態 A 和 B ，即其在 Boltzmann 分布中機率都不為零的情況，則當 A 到 B 作 transition 時 (A → B)，其 transition probability P (A → B) 必不為零，即 P (A → B) ̸= 0 或是說時間夠久，transition 一定可以到樣本空間中的任一狀態。但實際上演算法需要設計能夠在有限的步驟達到各態歷經。. 另外，細致平衡的意思就是從 A 狀態改變到 B 狀態的機率，與 B 狀態改變到 A 狀態的機率，兩者的比值等於總權重 W (B) 和總權重 W (A) 的比值， W (B) P (A → B) = P (B → A) W (A) 其中. W (A) W (B). = S(A/B) e−(E(A)−E(B)/kT ，E(A) 及 E(B) 分別為狀態 A 及 B 所對應的能. 量，而 S(A/B) 則是與選擇的機率相關的權重，k 是 Boltzmann 常數。 13.

(30) 3.1.2 Metropolis 演算法在討論所謂的取樣法之前，先說明一下權重的重要性 [40] [41]。當我們用蒙地卡羅法做簡單取樣來計算積分，若該函數為常數函數，則無論取樣數多少，準確度皆為百分之百；但若函數為 δ 函數，則無論取樣數多少，準確度幾乎皆為零。即，若在積分區間內，函數為一平滑函數，則簡單取樣的誤差很小；但若函數變動很大，則誤差就會很大，而且收斂也會很慢。因此取樣時可依據一個函數 f (x) 的量值大小來決定。這就是需要考慮權重 (weight) 的原因之一。權重函數 (weight function) W (x) ，是用來決定每個 x 被選中的機率，另外 W (x) 在積分範圍內要歸一化，即 ∫ W (x)dx = 1 。 Metropolis 演算法是 Nicholas Metropolis 先生在 1953 年命名的。下列我們將簡單說明如何在計算積分時應用 Metropolis 演算法，即如何產生一連串限制在積分範圍內的亂數 x (x′ ) ，使的 x (x′ ) 的分佈滿足權重 W (x) (W (x′ )) 的要求。 1. 由均勻亂數產生一任意 x， 2. 再產生另外一個亂數 x′ ， 3. 兩者的機率權重比值為 P (x→x′ ) =. W (x′ ) W (x). 4. 若 P ≥ 1 則決定接受 x′ 為新的 x 值，若 P < 1 則再產生另一亂數 y ，若 y < P 則決定接受 x′ 為新的 x 值，不然保持舊的 x 值 5. 如此即可得到一連串的 x ，重複此步驟 N 次，將每次所得到的 x 所對應的 f (x)/W (x) 求平均，即為所求但新舊的 weight 的比值有可能超過 1 ，所以可以將步驟 3 改寫成 ( ′. P (x→x ) = min. ) W (x′ ) ,1 , W (x). 我們可以驗證 Metropolos 演算法會滿足 ergodicity 和 detailed balance [20] [21]。. 3.1.3. 隨機級數展開法 SSE. 這節所提的隨機級數展開法 (Stochastic series expansion) ，縮寫成 SSE ，主要是參照參考文獻 [18] [22]−[25] [28] [42] [43]。 14.

(31) 在量子統計下，可以用 partition function ，Z ，計算在溫度 T 之下 operator A 的期望值 ⟨A⟩ ， Z = T r{e−βH },. ⟨A⟩ =. (3.1). 1 T r{Ae−βH }, Z. (3.2). 其中 β = 1/T ，H 為所考慮的 Hamiltonian 且設所有自然常數單位為 unity 1 。所以最重要的問題是我們如何計算 partition function 中的指數函數，一個簡單的想法就是對它用泰勒展開：. e. −βH. ∑ (−β)n (−1) (−1)2 2 = (1 + (βH) + (βH) + · · · ) = Hn 1! 2! n! n=0 ∞. 如果我們進一步選用適當的基底. Z=. ∑ α. |α⟩ ，則 Z 可以寫成，. ∞ ∑∑ βn α. n=0. (3.3). n!. ⟨α|(−H)n |α⟩,. (3.4). 或者寫得精確一點，因為 operator H 彼此之間不能交換，所以我們在每 2 個 H 之 ∑ 間再插入一組所選取的基底 i |αi ⟩⟨αi |，則我們會有. Z=. ∞ ∑ n=0. (−1). nβ. n. n!. ∑. ⟨α0 |H|αn−1 ⟩ · · · ⟨α2 |H|α1 ⟩⟨α1 |H|α0 ⟩,. (3.5). {α}n. 提醒一下 (3.1)-(3.5) 式是一般 general 的情況，以下我們會針對在正方晶格上的量子 Heisenberg model 應用以上所引入的觀念。這樣的應用是由 Prof. Anders Sandvik 所發展出來的，在文獻上稱為 Stochastic Series Expansion (SSE) 方法。量子 Heisenberg model 的 Hamiltonian 為. H=J. ∑. Si · Sj. (3.6). ⟨i,j⟩. 其中 ⟨i, j⟩ 表示是最近鄰的 sites，J 是反鐵磁 coupling ，且 J > 0 ，Si 為在 site i 的 spin-1/2 operator 。SSE 本身並沒有限制，但為了在泰勒展開後都是正值，避免 sign 15.

(32) problem 。這裡要求必須是 bipartite，非阻挫 (non-frustrated) 的晶格點，bipartite 是像西洋棋盤一樣黑白相間排列的規則來將正方晶格分為兩類 sub-lattices 。這裡我們將以二維正方晶格為例來說明 SSE 的精神及其 implementation ，我們假設晶格大小為 N = Lx Ly ，且其晶格有週期邊界條件。此外，我們選擇 z 軸為 quantized 的方向，即 Siz = ± 21 。則自旋 S = 1/2 的 N 個 spin 晶格其 basis 可表示為， z |α⟩ = |S1z , S2z , · · · , SN ⟩,. Siz = ±. 1 2. (3.7). 注意到，利用升降 (ladder) 算符，Hamiltonian 可寫成，. H=J. B [ ∑. z z Si(b) Sj(b). b=1. ) 1( + − − + + S S + Si(b) Sj(b) 2 i(b) j(b). ] (3.8). 這裡 b = 1, · · · , B，是表示連接兩個最相鄰的 spin sites i, j ，B 表示 couplings 的總數目，對週期邊界條件的二維正方晶格，couplings 的總數 B 為 2N 。由 (3.8) 式，我們可看到 Hamiltonian 可分成兩類 operators ，即 diagonal operators (H1,b ) 和 off-diagonal operators (H2,b )，. H1,b =. 1 z z − Si(b) Sj(b) , 4. H2,b =. ) 1( + − − + Si(b) Sj(b) + Si(b) Sj(b) 2. (3.9). 藉由 diagonal 及 off-diagonal operators ，Hamiltonian ( 3.8) 式可改寫為. H=−. B ∑ 2 ∑. Ha,b. (3.10). b=1 a=1. 注意到這裡 diagonal operators 前面我們多加了常數 1/4 ，使得其所對應的矩陣元素恆為正值，這樣的 shift 不會影響到我們要計算觀測量的期望值。，另外 off-operator 前面多了負號，這可看成將 x 和 y 的 spin operators 轉了 180 度，可以證明這樣的情況並不會改變系統的物理，特別是其為正定的特性 (以下我們會簡單解釋正定這一特性，正因為這一特性，我們可以應用 Monte Carlo method 來研究在 bipartite lattices 上的 Heisenberg model)。注意到，對新的 Hamiltonian 的每一項，無論是 diagonal 或是 off-diagonal 需要作用在一對其 S z 值為不同的 spin 上才有不為零的. 16.

(33) 矩陣元素。的確由 ( 3.9) 式，我們有 ⟨↑i(b) ↑j(b) |H1,b | ↑i(b) ↑j(b) ⟩ = 0,. ⟨↑i(b) ↑j(b) |H2,b | ↑i(b) ↑j(b) ⟩ = 0. ⟨↓i(b) ↓j(b) |H1,b | ↓i(b) ↓j(b) ⟩ = 0,. ⟨↓i(b) ↓j(b) |H2,b | ↓i(b) ↓j(b) ⟩ = 0. 及 1 1 ⟨↑i(b) ↓j(b) |H1,b | ↑i(b) ↓j(b) ⟩ = , ⟨↓i(b) ↑j(b) |H2,b | ↑i(b) ↓j(b) ⟩ = , 2 2 1 1 ⟨↓i(b) ↑j(b) |H1,b | ↓i(b) ↑j(b) ⟩ = , ⟨↑i(b) ↓j(b) |H2,b | ↓i(b) ↑j(b) ⟩ = , 2 2 由上式我們可明白在 diagonal operators 加上常數 1/4 ，會使得 diagonal 及 offdiagonal 的矩陣元素都是是 1/2 ，這會大大簡化我們的計算。將 Hamiltonian( 3.10) 式代回一般情況下的 partition function (3.4) 式，. Z=. ∑. ⟨α|e−βH |α⟩ =. α. ∑ α. ⟨α|1 +. (−1) (−1)2 (βH) + (βH)2 + · · · |α⟩ 1! 2!. 則我們會得到. Z=. ∑. )] )]2 ⟩ [ ( B 2 [ ( B 2

(34)

(35) ∑∑ ∑∑ (−1)2 (−1)

(36)

(37) α

(38) 1 + β − Ha,b + β − Ha,b + · · ·

(39) α 1! 2! b=1 a=1 b=1 a=1. ⟨. α. 在實際應用上，在展開的 partition function 中的高次方項的貢獻非常小，所以通常我們會選擇一個很大的正整數 M 來 truncate 泰勒展開，這個正整數 M 稱為 cut-off 。另外，為了方便計算，小於 M 次方的部份，會補上 unit operators 即 H0,0 = 1 ，使得展開中的每一項其 order 均為相同的 M 次方。考慮這些因素後，上式可改寫為一連串的算符的乘積， ∑ ∑ β n (M − n)! Z= M! α {Hab }. 為了方便討論，我們稱. ∏M p=1.

(40) ⟩ ⟨

(41) M

(42)

(43) ∏

(44)

(45) Ha(i),b(i)

(46) α α

(47)

(48)

(49). (3.11). i=1. Ha(p),b(p) 為 (M ) 算符串。這裡要注意的是在執行演算 17.

(50) 法時，M 的初始值一開始設定不會太大，但重複蒙地卡羅運算時會適當的調整 M 的值，使得要計算的精確度和誤差在控制之內。而 M 算符串中被 p 個算符作用之後的 propagated state |α(p)⟩ 記為， |α(p)⟩ =. p ∏. Ha(j),b(j) |α⟩. (3.12). j=1. 由於 propagated state 必須滿足邊界週期性條件，即 |α(M )⟩ = |α(0)⟩，因此在展開時，對每一個 site 我們會有偶數個與其相關的 off-diagonal operators ，而 diagonal operators 的數目並無任何限制。如果我們將算符串表示為 {Ha,b } ，並將隨機採樣後的樣本稱為 configuration 並記為 (α, {Ha,b }) ，則 (α, {Ha,b }) 指的是不同狀態搭配不同算符串。而 configuration 的集合則是 configuration space 。其隨機採樣的機率是 P = W /Z ，即各自的 weight 除以 partition function ，對於一個包含 n 個非 unit operators 且是所允許的算符串其所對應的 configuration 的 (Boltzmann) weight W (α, {Hab }) 為 ( )n β (M − n)! W (α, {Hab }) = 2 M!. (3.13). 在 propagating 的過程，configuration 所允許 operators 作用的情況，可用四種圖示 (圖 3.1) 來表示：. Figure 3.1: 四種在 bonds 的作用下可能的 spin 狀態其中綠色正立三角形點表示 spin up (+) ，藍色倒立三角形點表示 spin down (−)，可將作用前後總共四個 spin states 看作 bond 的四支腳 (leg)，依序標記腳位編號 0 ， 1 ，2 ，3 為左上，右上，左下，右下。. 18.

(51) 圖 3.2 表示的是某一個 configuration 的 propagating 過程的圖像。其中空白的 states 表示 unit operators ，而 diagonal 及 off-diagonal 則是用 bond 來標示，紅色白色相間棒代表 diagonal ；空心虛線棒代表 off-diagonal 。圖 3.2中將左圖簡化為右圖，沒有變化的部份僅以直線標示。另外除了顯示初末狀態之外也標示了在 bonds 作用前後的 spin 狀態。. Figure 3.2: propagated state 而在 propagated states 中，每個 operator ，包含 unit operator ，都有四支腳。這裡可以設計儲存連接腳位的表，從第一個 propagated state 的 bond 依照先後順序開始編號，之後在所謂的 loop update 中就可以透過這個連接腳位的表來做 update。SSE 是一個非常有效率的演算法的原因，在於它保持 propagated state 初末狀態的邊界週期性，藉由 diagonal update，插入或移除 bond，然後連接 bond 並做 loop update 。這些過程可以一次同時改變許多的 spins ，而非針對單一的 spin site 去考慮。以下我們要簡單說明這些步驟。. 19.

(52) Diagonal update Diagonal update 的目的是在 propagated states 中改變 diagonal bond 的數目。最簡單的方法就是將 diagonal operator 換成 unit operator ，或是將 unit operator 換成 diagonal operator ，這兩個過程我們稱為對 diagonal operator 的 removal 及 diagonal operator 的 insertion 。接著我們將說明如何執行這兩個過程，如在蒙地卡羅 important sampling 所進行的步驟一樣，我們從一些任意允許的 configurations (也就是非零的 weight W ) 來產生 Markov chains 。用的演算法可為 Metropolis 方法。在 Metropolis 方法中，接受一舊的暫存 configuration (其所對應的 weight 記為 Wold ) 改變到一個新的暫存 configuration (其所對應的 weight 為 Wnew ) 機率是， ( P = min. ) Wnew ,1 Wold. 所以我們要計算新舊 weight 的比值，weight 僅與 propagated states 中非 unit operators 的 operators 的數目 n 有關，我們可從式 (3.13) 中算出插入一個 bond 時新舊 weight 的比值是： W (n + 1) β/2 = W (n) M −n 注意到，當插入 bond 時，我們選擇任一 bond b 的機率都是均等，即為 1/B 。另一方面，移除 bond 時我們總是選取已存在此 propagated state 的 bond 。所以插入一個 bond 接受的機率是 ( P (n → n + 1) = min. Bβ/2 ,1 M −n. ). 這裡要注意的是，要求插入 bond 僅允許發生在其所連接的 sites 的 spin 是相反的。同樣的，移除 bond 的 weight 的比值及移除 bond 的機率也可以從式 (3.13) 中算出，分別是 M −n+1 W (n − 1) = W (n) β/2 ( P (n → n − 1) = min 20. ) M −n+1 ,1 Bβ/2.

(53) Loop update. Figure 3.3: loop update 改變 bond 的種類及翻轉其所對應的 spin. Figure 3.4: loop update 跨過邊界的情形在 loop update 主要的作用是在 configurations 中，將 diagonal operators 和 offdiagonal operators 互換及將 loop 上的 spins 做翻轉。由於 propagated states 滿足週期性邊界條件，在一個 configuration 中，作用在某一個 site 上的 off-diagonal operators 的數目必須是偶數。我們可在圖 3.3，特別是橘色線所經過的部份，看到 loop update ，確實可以一次同時改變許多的 spins (及 bond 的種類)，因此可以提高改變 21.

(54) configuration 的效率。在圖 3.4 中，可以看到 loop 的路徑跨越了 boundary states 可以改變 boundary state |α⟩。也就是說在 loop update 的這個步驟中，實際的物理狀態確實被改變了。. Figure 3.5: single spin flip update 在圖 3.5 中顯示了在 loop update 中，即使 configuration 中沒有 operator 作用到的單一 spin ，也有機會被翻轉使得 configuration 改變。最後由於 diagonal operators 和 off-diagonal operators 的 matrix element 都是相同的 1/2 ，所以決定是否 flip 每個 loop 的機率就是 1/2 。這裡僅介紹 SSE 重要的精神，若要明白更詳細的理論推導，可以參考文獻 [18] [22]−[25] [28] [42] [43]。. 22.

(55) 3.1.4 Measurement in SSE 對 3D ladder-dimer 及 bond disordered quantum Heisenberg model 的研究，我們需要計算相關的物理觀測量，如 spin stiffness ρs ，staggered magnetization Ms ，spin wave velocity c 等。如何計算這些物理量，會在第四章中詳細介紹。由於我們所考慮的物理量 ρs 及 c 與 winding numbers squared ⟨W 2 ⟩ 有關，所以在這裡我們先說明在 SSE algorithm 下，這些 winding numbers squared ⟨Wi2 ⟩ ，⟨Wt2 ⟩ 的定義及測量 [18] [26] [46] [47]。當中對 spatial winding numbers squared ，我們有， ⟨Wi2 ⟩. ( )2 1 + − =⟨ (N − Ni ) ⟩ Li i. (3.14). 其中，Li 是 lattice 在 i 方向的長度，i ∈ 1, 2, 3 代表空間三個方向。Ni+ 與 Ni− 分別代表 propagated states 中 i 方向的 Sa+ Sb− 以及 Sa− Sb+ operators 的數目，這裡 {a, b} 代表在 i 方向上 bond 連接相鄰近的 sites 。藉由求出 spatial winding number squared，我們可以得到 spin stiffness 。另外在模擬時，在 propagated state 中 configuration 的演化，可以將它看成時間的維度，或者說是 β（溫度）的維度，它與空間一樣也有一個 winding number 要考慮，稱為 temporal winding number squared ⟨Wt2 ⟩，其定義為. ⟨Wt2 ⟩ = ⟨. ( N ∑. )2 Siz. ⟩. (3.15). i=1. 這裡 Siz 是在 i 位置的 spin 沿 z 方向的分量。在 spin operator 中，將 z 軸選擇為測量方向較為方便，因為在空間中無論座標軸如何選擇，僅有一個座標軸所對應的 spin 分量有較簡潔的形式。. 3.2 誤差分析 Error analysis 由於我們所用的方法是蒙地卡羅模擬，所以需要考慮誤差分析，在這節我們將簡單介紹 binnning ，Jackknife 及 bootstrap 等重要的概念 [20]。 23.

(56) 3.2.1. 數據分級 binning. 數據分級 (binning) 是預先處理數據，用來減少誤差估算偏差的技巧 [18] [48]。這是一個被廣泛使用的技巧。比如說在處理影像，binning 就是將一群畫素變成一個畫素，比如使用 2 乘 2 的 binning ，以四個畫素代表一個單一較大的畫素，用來描述所有的畫素。這種聚集的減化，以減少數據的數量，方便於分析，雖然丟掉一些資訊，但也可以減少資料的雜訊。由於我們不知道所產生數據之間的關聯性，我們在處理蒙地卡羅中的原始資料數據時，為了減少統計誤差，也使用了 binning 的方法。 Binning 的概念如下，比如說有 100 個數據，a1 , a2 , · · · , a100 ，可以作 2 個 1 bin ，即 a′1 = (a1 + a2 )/2, · · · , a′50 = (a99 + a100 )/2 ，這樣 a′1 與 a′2 所代表原始數據的間距也拉大，我們可以以此精神繼續作下去，如 4 個為 1 bin，a′′1 = (a1 + a2 + a3 + a4 )/4, · · · ，越多數據作 1 bin 可將數據間隔拉大。對每個 binning 的過程，我們計算其所對應的誤差，待誤差不再隨 bin 更多的數據而增加，即表示所對應的物理量之間的關聯性就可以忽略，如此可避免隨機亂數造成突然數值的起伏及減少連續測量物理量所造成的鄰近數據的相關性。Binning 是一種可有效減低雜訊 (noise) 的方法。. 3.2.2 Jackknife 分析將時間序列分割最簡單的想法是：將時間級數等分，然後對每個序列計算所需要之函數的值。統計上，Jackknife 是一種重取抽樣 resampling 的技巧。Jackknife 技巧是有系統逐次的扣掉原始數據集合內其中一個數據，加總平均後 resample 產生新的數據，並由新的數據計算平均及誤差。Jackknife 是由 Maurice Quenouille 在 1949 年發展的，John. W. Tukey 在 1958 擴展了該技術，涵義是說，像童子軍的折刀，它是一種“粗略的”工具，但可以解決各種問題，甚至有時比專門設計的工具更有效率。特別是對於計算要用到不同物理量的延伸物理量的誤差時，非常可靠。Jackknife 可視為 bootstrap 的線性近似 ( bootstrap 將於下一小節介紹)。例如有 M 個序列的函數，. f (x1 , y1 ), f (x2 , y2 ), f (x3 , y3 ), · · · , f (xM , yM ), 24.

(57) 原來的平均是. M 1 ∑ f (xi , yi ) U0 = M i=1. Jackknife 之後第 1 個數據，就是原數據第 1 個函數去掉不要加總之後的平均，為 1 ∑ U1 = f (xi , yi ) M − 1 i=2 M. 同理第 j 個新數據是原數據第 j 個函數不要加總之後的平均： 1 ∑ f (xi , yi ) Uj = M − 1 i=1 M. i̸=j. 在作完 Jackknife 之後的新數據集合，其平均會趨近於原來函數的平均。Jackknife 的平均 U ，期望值 ⟨U ⟩ 及誤差 ∆U 分別為 M 1 ∑ U= Ui M i=1. ⟨U ⟩ ≈ U0 − (M − 1)(U − U0 ) v u M uM − 1 ∑ ∆U = t (Ui − U )2 M i=1. 3.2.3 Bootstrap 分析 Bootstrap 原意是提靴帶，自己抬自己的腳走路，用在電腦軟硬體系統，指的是系統如何辨認硬體週邊，傳遞訊息，進行開機的過程；在演算法的關聯是一種 resampling 的技巧，bootstrap 由 Bradley Efron 於 1979 年發表。當樣本的 sampling distribution 不是常態分佈 (the normal distribution)，為了改善誤差，則以 boostrap resampling 等方法來分析。它採用隨機將樣本重新抽樣的方式 (random sampling with replacement)。對於原始數據的資料量不夠多的情況，使用 bootstrap 的效果很好。Bootstrap 即是在原來樣本數中反覆抽取自己的樣本，來計算標準差，信賴區間和偏差。舉例說明，如原始樣本數為 (x1 , x2 , x3 , x4 , x5 ) ，則 resample 後新樣本可能是 x⋆1 = (x2 , x3 , x5 , x3 , x1 ) ，x⋆2 = (x4 , x1 , x2 , x3 , x3 ) ，· · · 。重覆這個過程，比方說到 100 次，然後計算 resampling 的平均值。如將 x⋆i 的平均值記為 µ∗i ，則我們的 25.

(58) 樣本數就有 µ⋆1 , µ⋆2 , · · · , µ⋆100 。然後就可以對新的樣本空間 {µ⋆1 , µ⋆1 , · · · , µ⋆100 } 進行相關的計算及分析。注意到其選取的方式是隨機且機率相等，其中有可能是重覆選取，有的卻沒有抽樣到，但在取樣次數足夠多時，新樣本的平均值與真正期望值之間的差距就可以忽略不計。舉一極端的例子說明其好處，例如某次考試有 4 個 100 分，1 個 0 分，我們可以想像考 0 分是因為失常，所以可手動會把它去掉，在大量數據的處理時，使用 bootstrap resampling 的技巧就會自動稀釋類似 0 分所造成的這種偏差。我們在計算 TN 時，就使用了 bootstrap resampling 的這個技巧。. 26.

(59) Chapter 4 Numerical results 數值結果 4.1 3d ladder 空間各向異性之海森堡模型我們所研究的是三維 ladder-dimer 量子海森堡模型 (quantum Heisenberg model) [35] ，此模型的漢米頓算符 (Hamilton operator) 可表示如下：. H=. ∑. J S⃗x · S⃗y +. ∑. J ′ S⃗x′ · S⃗y′. (4.1). ⟨x′ y ′ ⟩. ⟨xy⟩. 其中 J(J ′ ) 是連接最近鄰自旋 (nearest neighboring spin) ⟨xy⟩(⟨x′ y ′ ⟩) 的反鐵磁交換耦合常數 (antiferromagnetic exchange coupling)。此模型的示意圖請參照圖4.1。為了研究由於 dimerization 引起的在臨界點 (critical point) 附近的 TN 和 Ms 之間的普適關係 (universal behaviour)，在我們模擬的計算中，我們測量了自旋剛性係數 (spin. J J’ Figure 4.1: ( 3 + 1 ) 維空間各向異性的 ladder-dimer 量子海森堡模型 27.

(60) stiffnesses ) ρsi ，其定義為. ρsi =. 1 ⟨Wi2 ⟩, i ∈ 1, 2, 3, βL1 L2 L3. (4.2). 其中 β 是溫度倒數，Li 指的是在 i 方向上盒子尺寸的長度，⟨Wi2 ⟩ 是在 i 方向上的 winding number 的平方。此外，Ms 的值是由記錄 ⟨(mzs )2 ⟩ 所得到。其中 ∑ m ⃗ s = L1 L12 L3 x (−1)x1 +x2 +x3 S⃗x ，這裡 mzs 指的交錯磁化 (staggered magnetization) m ⃗ s 的 z 分量。如同前面所提到的，我們使用了 SSE (stochastic series expansion with operator-loop update) 演算法 [25]，在我們的計算中，我們使用了許多不同 J ′ /J ，溫度倒數和 box-size L ，另外在整個計算中 J 設為 1 。需注意到的是，對此模型而言其 dimerization 所引發的 QCP 是在 (J ′ /J)c ∼ 4.0 [49]，所以我們進行計算的範圍為 2.5 ≤ J ′ /J ≤ 4.0 。接著在這節中我們會陸續說明所得到的數據及結果。. 4.1.1 Néel 溫度 TN 的計算當 T > TN 時，反鐵磁的長程有序會被破壞掉，我們在計算某個固定 J ′ /J 值時其所對應 Néel 溫度 TN 的方法如下，對固定 J ′ /J 值，我們取 L = 8, 12, 16, · · · , 36, 40 及不同的 T 值。而 TN 就可以從 finite-size scaling 的方法得到。具體的說，在靠近 TN 及對 i ∈ 1, 2, 3 的觀測量 ρsi L ，不同 L 對 T 的曲線應該會交錯在 TN 。有趣的是，對我們所考慮的 J ′ /J 值 (J ′ /J = 2.5, 3.0, 3.25, · · · , 3.875) ，適用於這些觀測量的 scaling ansatz 的修正項在 L ≥ 20 是可忽略的。換句話說，我們可以用 leading 的 scaling ansatz 來求 TN 。這裡我們所用的 scaling ansatz 型式為 g(x) ，其中 g 是一個參數為 x 的平滑函數，而 x 包含一個線性因子 (T − TN )/TN ，在實際應用上通常會對 g(x) 作泰勒展開。以 J ′ /J = 3.5 為例 (參看圖 4.2)，從 ρs L 這個物理量中，我們得到 TN = 0.7751(2) ，其中，ρs L = (ρs1 + ρs2 )L/2 (圖 4.2中的上圖)。應用相同的步驟，從 ρs3 L 中，我們得到 TN = 0.7750(2) (圖 4.2中的下圖)。值得注意的是從這兩個不同的觀測量所得到的 TN 相當一致。在其他 J ′ /J 值所對應的 TN ，也是利用上述方法得到，在表 4.1我們列出所有結果。另外在做 fitting 時我們使用了類似 bootstrap resampling 的技巧來估算 TN 及其所對應的誤差。 28.

(61) 0.6 0.5. ρsL. 0.4 0.3 0.2 0.1 0 0.75. 0.76. 0.77. 0.78. 0.79. 0.8. 0.78. 0.79. 0.8. T 1.2 1. ρs3L. 0.8 0.6 0.4 0.2 0.75. 0.76. 0.77. T Figure 4.2: 蒙地卡羅在 J ′ /J = 3.5 時所得到的 ρs L (上圖) 及 ρs3 L (下圖) 之數據，虛線僅是引導作用，盒子尺寸為 L = 8, 12, 16, · · · , 32, 36 和 40，其中， ρs L = (ρs1 +ρs2 )L/2。. 29.

(62) Observable ρs L ρs3 L ρs L ρs3 L ρs L ρs3 L ρs L ρs3 L. J ′ /J. 2.5 2.5 3.0 3.0 3.25 3.25 3.375 3.375. TN. J ′ /J. TN. 1.0014(2) 3.5 0.7751(2) 1.0014(2) 3.5 0.7750(2) 0.9317(2) 3.625 0.7087(3) 0.9316(2) 3.625 0.7086(3) 0.8690(2) 3.75 0.6197(2) 0.8689(2) 3.75 0.6193(3) 0.8270(2) 3.875 0.4853(3) 0.8269(2) 3.875 0.4849(4). Table 4.1: 利用 finite-size scaling 所得到的 TN 的值。所有擬合的 χ2 /DOF 皆小於 1.6 。. 4.1.2. 交錯磁化強度密度 Ms 的計算. 我們可從觀測量 ⟨(mzs )2 ⟩ 來計算 Ms 。具體的說，對有限晶格大小的 ⟨(mzs )2 ⟩ 做外 √ 插可得到 ⟨(mzs )2 ⟩(∞) ，而 Ms 可以由 Ms = 3(mzs )2 (∞) 進一步得到。利用此方法來計算 Ms ，我們需要在絕對零度的 ⟨(mzs )2 ⟩ 的值。我們在 J ′ /J= 2.5, 3.0, 3.125, 3.25, 3.375, 3.5, 3.625, 3.75, 3.875 ，計算了兩組數據來當作測試，一是 L = 20 及 βJ = 20；另一是 L = 20 及 βJ = 40 。在我們所考慮的耦合 J ′ /J ，這兩組不同溫度倒數 β 所得到的 ⟨(m2s )⟩ 的值是相當一致的。所以在 βJ = L 的條件下進行模擬的計算所得到的 Ms 值應該是可信的，參考文獻 [34] 有提到類似的計算。在圖 4.3中，所顯示的是在所考慮的 J ′ /J 下，L = 6, 8, 10,..., 32, 36, 40 所對應的 ⟨(m2s )⟩ 的值。我們做外插法時所使用的公式是 a + b/L + c/L2 + d/L3 。圖 4.4 的圓點為我們所得到的 Ms 值。另外在圖 4.4中的實線是從參考文獻 [31] 中重製而來的。從圖 4.4中，我們可以清楚看到蒙地卡羅的結果和 [31] 的結果吻合度相當好。. 4.1.3. 自旋波速 c 的計算. 有許多方法可以用來計算自旋波速 c 。這裡我們使用 winding number squared 的這個想法。具體的說，對每個 J ′ /J 我們調整 L1 /L3 的比例，使得三個空間方向的 winding numbers squared 大約是相同的值 [44] [50]。然後去調 β 讓在 i ∈ 1, 2, 3 這三個方向都可以達到 ⟨Wt2 ⟩ ∼ ⟨Wi2 ⟩ 的這個條件，這裡 ⟨Wt2 ⟩ 是溫度方向 (temporal) 的 winding number squared 。一旦達到此條件，數值 c 就可以從以下不等式估算， L/β2 ≤ c ≤ L/β1 ，上式中 L = (L1 L2 L3 )1/3 ，β1 (β2 ) 是可使得 ⟨Wi2 ⟩ ≤ ⟨Wt2 ⟩ ， 30.

(63) 0.06. 0.04. z 2. < (ms ) >. 0.05. 0.03 0.02 0.01 0 0. 10. 20. 30. 40. L Figure 4.3: 由蒙地卡羅計算所得出的 ⟨(m2s )⟩ 對 L 的關係圖，從上而下對應的 J ′ /J 依序是 2.5, 3.0, 3.125, 3.25,· · · ,3.75 和 3.875，實線僅是引導作用。. 0.4 0.35 0.3. Ms. 0.25 0.2 0.15 0.1 0.05 0. 2.9. 3. 3.1. 3.2. 3.3. 3.4. 3.5. 3.6. 3.7. 3.8. 3.9. 4. J’/J Figure 4.4: 從蒙地卡羅計算所得到的 Ms 的值。實線則是從文獻 [34] 重製而得到的。. 31.

(64) J ′ /J. L1 2.5 22 2.5 36 3.0 32 3.0 44 3.25 12 3.25 18 3.25 24 3.375 12 3.375 24 3.5 34. L3 28 46 46 58 16 24 32 16 32 46. J ′ /J. c/J. 2.215(8) 3.5 2.215(9) 3.625 2.282(13) 3.625 2.283(11) 3.75 2.317(7) 3.75 2.317(8) 3.875 2.317(11) 3.875 2.335(12) 4.0 2.334(13) 4.0 2.347(12) 4.0. L1 46 22 34 22 44 16 32 16 32 42. L3 62 30 46 30 60 22 44 22 44 58. c/J. 2.348(10) 2.360(12) 2.360(13) 2.376(12) 2.378(11) 2.391(7) 2.389(8) 2.408(13) 3.408(15) 2.401(10). Table 4.2: 對各種耦合 J ′ /J ，藉由 winding number squared 所求出的自旋波速 c 的數值。. i ∈ 1, 2, 3 ( (⟨Wi2 ⟩ ≥ ⟨Wt2 ⟩) ，i ∈ 1, 2, 3 ) 這個條件成立的最大 (最小) 之溫度導數。在各向同性的情況，也就是 J ′ /J = 1.0 ，自旋波速理論預測 c ∼ 1.9091J [51]。值得一提的是，在我們測試的模擬中，當 J ′ /J = 1.0，L1 = L2 = L3 = 20 ，βJ = 10.476 (由 L/β ∼ 1.9091J 及 L = 20 等條件可得出 βJ = 10.476) 時，這三個空間與溫度的 winding numbers squared 的比值大約是 0.994。這確認了用這個方法去求 c 的想法是正確且十分有效的。另外，對所考慮的每個 J ′ /J 值，在估算 c 時，我們至少考慮了兩組盒子尺寸，表 4.2列出了對 J ′ /J=2.5, 3.0, 3.25, 3.375, 3.5, 3.625, 3.75, 3.875, 4.0 我們所得到的 c 的數值。表 4.2的結果顯示對所考慮的耦合 J ′ /J ，其對應的 c 值已經收斂。即使我們部分的 c 沒有真正收斂，我們可以預期這偏差必定極小。因此這種系統不確定的因素並不會對我們想研究的 TN 和 Ms 之間的普適關係有太大的影響。. 4.1.4. 理論預測與蒙地卡羅的比較. 參考文獻 [31] 用相應的場論預測在 QCP 附近，TN 和 Ms 之間普適關係是 √ TN =. 12c3 Ms 5. (4.3). 注意到原本參考文獻 [31] 中，應用在空間各向異性系統的預測關係中的 c3 是 c1 c2 c3 ，這裡 ci 指的是在 i 方向自旋波速。另一方面，我們在前一節中有詳細描述如何測量 c 。乍看之下，我們所想驗證的 TN 和 Ms 的關係式似乎與原版的有 32.

(65) 1.2 1. TN. 0.8 0.6 0.4 0.2 0. 2.8. 3. 3.2. 3.4. 3.6. 3.8. 4. J’/J Figure 4.5: 圖上顯示 TN 對 J ′ /J 藉由蒙地卡羅 (實心圓點)，級數展開 (實心方塊)，及場論計算 (實線) 所得到的關係圖。場論及數值展開的結果是從文獻 [31] 估計和重製，dashed 及 dotted 線僅是引導作用。. 些不同。而實際上，我們所驗證的關係式和原版是等效的。這可由下列敘述來理解。要計算相對應在各方向 i ∈ 1, 2, 3 中的自旋波速 ci ，在模擬中我們需要有一個方法，而此方法對任何空間方向不會有所偏好。在此條件下，最自然的選擇條件是要求使 ⟨W12 ⟩ ，⟨W22 ⟩ 及 ⟨W32 ⟩ 有相同的值。要使這個條件成立的方法是調整所模擬中盒子尺寸的比值。一旦達成這條件，相對應盒子尺寸邊長 Li 的 ci 可以藉 ci (Li ) = Li /β 的關係給出，其中 β 是溫度倒數。所以對每個方向 i ∈ 1, 2, 3 的 ⟨Wi2 ⟩ 皆會與 ⟨Wt2 ⟩ 相等。而此方法所求的 ci 的收斂性是以由計算盒子尺寸 Li 及 2Li 所得到的自旋波速來決定。如果在 Li 及 2Li 所計算的自旋波速相吻合，則此 ci 即是收斂的。注意到文獻 [31] 原預測的式子中有一項是 c1 c2 c3 ，對一體積是 L1 L2 L3 的有限晶格，c1 (L1 )c2 (L2 )c3 (L3 ) 可由 L1 L2 L3 /β 3 給出，而這正是這裡使用的觀測量 c3 的定義。所以這裡定義的 c3 在熱力學極限下，可以逼近正確的值 c1 c2 c3 。也就是說我們研究的進行方式原則上與 [31] 原始預測的關係是等效的。注意在實際的模擬中，這三個空間 winding numbers squared 取相同值的條件在離散的有限晶格僅可大略近似，這也是為什麼我們使用這方法計算 c 時，需要考慮 2 個不同的溫度倒數來估算 ci 的值。對 Ms 來說，在圖4.4中可以看到，我們計算的結果與文獻 [31] 數值的結果十分吻合，另外在圖4.5可以看出，我們蒙地卡羅計算的 TN 與從參考文獻 [31] 級數展開方法所得到的結果，在接近 QCP J ′ /J ∼ 4.0 時相當一致。最後方程式 (4.3) 意味 33.

(66) 0.4 0.35. TN / c. 1.5. 0.3 0.25 0.2 0.15 0.1 0.05 0 0. 0.05. 0.1. 0.15. 0.2. 0.25. 0.3. 0.35. 0.4. 0.45. Ms √ Figure 4.6: 圖上顯示 TN / c3 對 Ms 的關係圖。實線是用 a + bMs 形式擬合的結果，虛線是用 b1 Ms + d1 Ms log(Ms ) 的擬合結果，正方形點是 J ′ /J = 1.0 的結果， J ′ /J = 1.0 的結果是參考文獻 [51] 和 [52] 的計算得來的。 √ √ 著 TN / c3 對 Ms 的曲線圖應該是線性的。在圖 4.6，我們考慮 TN / c3 對 Ms 的關 √ 係圖。的確從圖 4.6中可看出，TN / c3 定性上是與 Ms 成線性關係。在圖 4.6中，對 √ J ′ /J = 2.5, 3.0,· · · , 3.875 用 a + bMs 去擬合 TN / c3 可以得到 a = 0.0117(20) 。這與期望值 a = 0 僅有非常小的差別。我們將這點偏差歸咎為修正項的效應。由於所得到 a = 0.0117(20) 跟零相差很小，我們可以預期對所考慮的參數 J ′ /J ，此修正 √ 項的效應應可忽略。有趣的是，從擬合得到的 b 的值只有預期 12/5 的一半。這需要更進一步的研究探討為何有此差別。雖然在 Ms 值小時的區域有些偏差，我們發現方程式 (4.3) 在使用我們所定義的 c 時確實是正確的。如果我們從量子色動力學 (quantum chromodynamics) 的手性微擾 (chiral perturbation) 理論的精神下來推演，修正項的型式應為 Ms log(Ms ) 。值得一提的，如果我們使用 b1 Ms + d1 Ms log(Ms ) 來擬合，(在圖 4.6中的虛線)，則擬合的結果非常好。特別是在我們實際所計算的範圍，這兩種擬合的公式可以說都和我們的數據非常吻合。所以在這裡，對我們 √ 所考慮的 3D 的 ladder-dimer 模型，我們驗證了在 QCP 附近，TN / c3 = AMs 這線性關係是成立的。. 34.

(67) 4.2 disorder 無序性在這一節中我們會描述及介紹所引入的反鐵磁耦合無序性 (antiferromagnetic coupling disorder) 及所測量的物理量 (observables)。這裡所考慮的三維量子海森堡模型 (quantum Heisenberg model) 可以用以下漢米頓算符 (Hamilton operator) 來表示： H=. ∑. Jij S⃗i · S⃗j. (4.4). ⟨ij⟩. 其中 Jij 是連接最近鄰自旋 (nearest neighbour spin) ⟨ij⟩ 的反鐵磁交換耦合常數， S⃗i 表示作用在標記位址 (site) i 上的 spin-1/2 算符。為簡單起見，在這節我們會稱 Jij 為 (反鐵磁) 鍵結 (bond) 。為了研究鍵結無序性 (bond disorder) 的效應對基態 (ground state) 性質的影響以及想瞭解 Néel 溫度 TN 及交錯磁化密度 (staggered magnetization density) Ms 之間的普適關係在無序系統是否會成立 [45]，我們引進兩個參數，D 和 P ，來描述鍵結 Jij 的隨機性。參數的 D 和 P 的值皆在 0 和 1 之間，即 0 ≤ D < 1 和 0 ≤ P < 1，而 bond disorder 的引入如下：每個交換耦合強度會有 P 的機率為 Jij (1 + D) ；(1 − P ) 的機率為 Jij (1 − D) (我們在模擬中，設 Jij = 1)。這裡需強調的是在我們作蒙地卡羅模擬的過程中，每個隨機組態都是由自己的隨機種子 (random seed) 所產生的。所以對不同組態所測量到的物理量彼此之間的關聯性應是極小的，另外我們執行模擬計算的過程如下：對任何固定的一組 (D,P )，晶格中每一個 bond 根據 (D,P ) 機率的要求，來決定每個 bond 的耦合強度。然後做 thermalization ，平衡後再作測量 (measurement)，例如做 200 次測量，接著作平均得到一組數據，重複這樣的過程來產生幾百個或幾千個數據，最後再作數據分析 [39]。舉例來說，固定一組 (D,P )，對第 1 個 bond，這次的耦合強度是 J(1 + D)，下一次的強度可能是 J(1 − D) 。藉由改變 D 和 P ，我們可研究反鐵磁如何隨 D 和 P 變化而減弱，進而也可探討 TN 和 Ms 值的普適 (線性) 關係是否對無序系統也成立。 Ms 的計算已在前一節所介紹，即 Ms =. √ S(π, π)(∞) ，其中 S(π, π) 是零溫的. staggered structure factor [53]，其在立方晶格的定義為. S(π, π) = 3. ⟨( 1 ∑ )⟩ i1 +i2 +i3 z 2 (−1) S i L3 i 35.

(68) Figure 4.7: β- doubling scheme 的示意圖，摘自文獻 [42]. D = 0.65, P = 0.20, L = 16 D = 0.75, P = 0.15, L = 14 D = 0.99, P = 0.10, L = 18 D = 0.90, P = 0.20, L = 12. 0.2. S(π,π). 0.15. 0.1. 0.05. 0. 1. 10. 100. βJ. 1000. 10000. Figure 4.8: 對於不同 L ，D 和 P ，S(π, π) 收斂到其基態值的示意圖. S(π, π) 與前一節所用的 ⟨(mzs )2 ⟩ 定義上的差別，只相差一個常數的倍數。除了 Ms 以外，另一個與長程 (long-range) 反鐵磁有序相關的物理量是 spin stiffness ρs 。也就是說，當 L → ∞ 的時候，ρs 的值不為零的話，則表示仍有反鐵磁的存在。而 ρs 的定義與測量方式也於上一節介紹過。為了要更有效率的得到所測量的物理量其基態的值，我們使用 β 倍增 (β- doubling) 的技巧，可參照圖4.7及參考文獻 [42]。接者我們將介紹我們在三維反鐵磁中鍵結無序性對其基態的效應的計算結果。 36.

(69) D = 0.4, P = 0.2 D = 0.5, P = 0.4 D = 0.75, P = 0.3 D = 0.95, P = 0.4 D = 0.9, P = 0.2 D = 0.95, P = 0.2 D = 0.99, P = 0.2. 0.3. S(π,π). 0.25 0.2 0.15 0.1 0.05 0 0. 5. 10. 15. 20. 25. 30. L Figure 4.9: 不同的 D 和 P ，其平均交錯結構因子 S(π, π) 對 L 的關係圖，虛線僅是引導作用. 0.5. Ms. 0.4. 0.3. D = 0.5 D = 0.75 D = 0.95 D = 0.99. 0.2. 0.1 0. 0.1. 0.2. 0.3. 0.4. 0.5. 0.6. 0.7. 0.8. 0.9. 1. P Figure 4.10: 不同的 D 值，其 Ms 對 P 的關係圖。在固定 D 之下，可看出 Ms 在 P ∼ 2 附近有最小值，虛線僅是引導作用. 37.

(70) D = 0.5 D = 0.75 D = 0.95 D = 0.99. 0.3 0.25. ρs. 0.2 0.15 0.1 0.05 0 0. 0.1. 0.2. 0.3. 0.4. 0.5. 0.6. 0.7. 0.8. 0.9. 1. P Figure 4.11: 不同的 D 值，其 ρs 對 P 的關係圖，虛線僅是引導作用. 4.2.1. 不同 (D，P ) 所對應之 Ms 和 ρs 的計算. 為了瞭解鍵結隨機性如何使反鐵磁減弱，我們對數百個無序組態作相關物理量的測量。特別是每個無序組態的產生都有自己開始的隨機種子 (seed)。因此這些從不同組態得到的物理量之間的關聯 (correlations) 應是在控制之下。此外我們也用了保守的誤差估計以減少了統計誤差 (statistical fluctuations) 可能造成的影響。先前提到的 β-倍增 (β- doubling) 方法可以有效率的用來計算基態性質。在圖 4.8可以看出，對不同組 D 和 P 的無序平均交錯結構因子 S(π, π) 如何收斂到其基態值。在得到 S(π, π) 的值之後 (在有限大小晶格中 L 與 S(π, π) 的關係，如圖 4.9所示) ， Ms 的 bulk 值則可用 a0 + a1 /L + a2 /L2 + a3 /L3 多項式的形式去擬合 (fitting) 而得到。在我們的模擬計算中，這裡盒子的尺寸 L 大部分是從 L = 4 到 L = 20 ，某些點最多到 L = 28 。圖 4.10中展示了對某些考慮的參數 D 和 P ，其所對應的 Ms 的值。對自旋剛性係數我們也使用了一個類似的分析，其結果在圖 4.11。值得注意的是，對每個 D 而言，ρs 的最小值發生時所對應的 P 值和 Ms 最小值發生時的 P 值並不相同，只有在考慮 ρs /J 時，其最小值所對應的 P 值才會和 Ms 最小值發生時的 P 值相吻合。在 ρs /J 中，J 為反鐵磁耦合強度的平均值。在我們的計算中，我們用了 J = (1 + D)P + (1 − D)(1 − P ) 這個定義。圖 4.10及 4.12顯示對我們所考慮的這種 disorder 的系統來說，ρs /J 才是衡量反鐵磁比較正確的物理量，而非 ρs 。 38.

(71) 0.2. ρs / J. 0.15. 0.1. D = 0.5 D = 0.75 D = 0.95 D = 0.99. 0.05. 0. 0. 0.1. 0.2. 0.3. 0.4. 0.5. 0.6. 0.7. 0.8. 0.9. 1. P Figure 4.12: 不同的 D 值，其 ρs /J 對 P 的關係圖，虛線僅是引導作用. 4.2.2. 不同 (D，P ) 所對應之 TN 的計算. 我們在計算不同的 D 和 P 值所對應的 Néel 溫度 TN ，與前一小節 4.1.1 節的方法類似：對固定的 (D,P )，在靠近 TN 附近，對 ρs L 這個物理量，其不同的 L 對 T 的曲線會交會在 TN 。只是這裡所用的 fitting 公式, 多了一個修正項，其形式為 (1 + b0 L−ω )(b + b1 tL1/ν + b2 (tL1/ν )2 + · · · ) 。其中 bi 中，i = 0, 1, 2, · · · 等都是常數而 t = (T − TN )/TN 。這裡用來計算 TN 的觀測量是 ρs L 。但圖4.12告訴我們說 ρs /J 才是適合描述此一系統的物理量，注意到對固定的 (D,P ) 而言，其 J 都是相同的固定值，因此由 ρs L 計算所得到的 TN ，應與從 ρs L/J 所得到的 TN 是一致的。的確，我們記錄了幾個由模擬所得到的 J 值，而這些 J 值與其理論值的誤差遠小於 ρs L 的誤差。所以我們所求出的 TN 是可靠的。在圖 4.13中所得到的交點，即 TN 的值在上圖是 0.34381(13) ，而下圖是 0.6628(2)。有趣的是從 fitting 所得到的 ν 的平均值大約是 0.705 。這與 O(3) universality class 的結果 ν ∼ 0.71 是一致的。要注意是，我們用的是可用來研究乾淨系統的 finite size scaling 的公式。雖然這裡我們只對 TN 的值有興趣，但從理論的觀點來看 ν 的值是否真的與 O(3) 的預測一致，或是我們得到的結果僅是巧合，在未來是一個有趣而值得研究的主題。 39.

(72) 0.24. D = 0.70, P = 0.20. ρs L. 0.2. 0.16. 0.12. 0.08 0.336. 0.34. 0.344. 0.348. 0.352. T/J 0.48. ρsL. 0.4. D = 0.50, P = 0.30. 0.32. 0.24. 0.16. 0.652. 0.656. 0.66. 0.664. 0.668. 0.672. T/J Figure 4.13: ρs L for L = 16, 20, 24, · · · , 36, 40 ，上圖為 (D, P )=(0.70, 0.20)，下圖為 (D, P )=(0.50, 0.30)，虛線僅是引導作用. 40.

(73) 4.2.3. TN 和 MS 之線性關係的驗證. 對三維 dimerized 量子海森堡模型，在靠近 QCP 有幾個 TN 和 Ms 之間相關的普適 (線性) 關係。舉例來說，在 4.1 節中，我們驗證了在臨界點附近 TN /c1.5 ∝ Ms 此一線性關係。此外，在 QCP 附近，如 TN /T ∗ 和 TN /J 對 Ms 也有類似的關係式存在 [34]。在這節中，我們將探討對無序性系統是否也有其類似的線性關係，在前面 4.2.1 小節中我們計算 ρs 的研究中，似乎暗示了考慮 TN /J 這一物理量是一個很好的選擇。從圖4.12，我們注意到當 D 增加時，Ms 的最小值發生在 P ∼ 0.2 。因此我們針對固定在 P = 0.2 ，及不同的 D ，考慮 TN 對 Ms 的關係，來驗證兩者之間是否有線性關係。在圖 4.14 下圖中，當 P = 0.2 時，我們看不出來 TN 與 Ms 有任何明顯的關係。令人意外的是，如果我們考慮的是 TN /J 對 Ms 的關係，則在乾淨系所對應的點附近，確實可觀察到 TN /J 對 Ms 有一個明顯的線性關係，這可從圖 4.14 上圖中得到驗證。為了避免可能發生的偏差，圖 4.15中，我們也考慮了 ladder-dimer model 的 TN /J 對 Ms 的關係當作對照。的確 ladder-dimer model 的 TN /J 對 Ms 在 QCP 附近有線性的行為，但在靠近各向同性點附近，並沒有線性關係。在圖 4.16中，我們將 ladder-dimer model 和固定 P = 0.3 的 bond disorder 的數據作一比較，可以更清楚看到對 bond disorder 的線性關係是比較接近各向同性點。而 ladder-dimer model 的 TN /J 和 Ms 的線性關係是發生於 QCP 附近。也就是說，靠近乾淨系統 TN /J 和 Ms 之間的線性關係是我們所考慮的這個無序系統的特性。. 41.