二維隨機格點伊辛模型中相變的蒙特卡羅模擬

劉 洋，趙 薇，周恒為，黃以能

（伊犁師范大學物理科學與技術學院，新疆凝聚態相變與微結構實驗室，新疆伊寧 835000）

0 引言

為了描述順磁-鐵磁相變，Lenz和Ising提出了現在一般稱之為伊辛模型（Ising Model，IM）［1-3］的著名理論模型［4］.IM的兩個主要假設為：1）伊辛自旋（spin）假設，即自旋只能處于向上或向下兩個狀態，這里的自旋表示是材料晶體點陣上的原子、分子、離子中電子的自旋磁矩和軌道磁矩的總磁矩，準確地說是贗自旋（Pseudo-spin），一般簡稱為自旋［5］；2）最近鄰相互作用假設，即晶體點陣上的自旋，只有在最近鄰的情況下，才存在相互作用.后來的研究表明，IM模型也是描述順電-鐵電相變最為成功的模型［6］.

Klein-Brout提出的長程相互作用隨機格點（random-site）伊辛自旋模型［7］，是為了描述磁性金屬固溶體材料（CuφMn1-φ、AuφFe1-φ等）而提出的.Klein-Brout 模型的3 個主要假設為：1）自旋和自旋空位（spin vacant）假設，即晶體點陣的任一格點上由自旋（如AuφFe1-φ中的Fe原子）或自旋空位占據，自旋空位定義為無自旋或自旋為0（如AuφFe1-φ中的Au原子）；2）點陣上自旋隨機分布假設，即在保證每個點陣格點上的平均濃度為1-φ條件下，自旋在點陣上隨機分布；3）自旋之間長程相互作用假設，即模型體系中，任意一對自旋之間的相互作用為RKKY勢［8-10］.

Binder等［5］用IM的最近鄰自旋相互作用假設，對Klein-Brout模型的長程相互作用進行簡化，即得到現在一般所謂的隨機格點伊辛模型（random-site-IM，RSIM）.Zhang-Huang［6］發現，在RSIM中再引入隨機場（random-field，RF），即RFRSIM，可以對弛豫鐵電性（relaxor-ferroelectricity）進行較好的描述.

目前，雖然已經獲得一維RSIM（1D-RSIM）體系熱平衡態的精確解，但是2D-RSIM、3D-RSIM仍未獲得相應的精確解［5，6］.已有的蒙特卡羅（Monte Carlo，MC）模擬結果［11-14］，則存在所模擬的自旋空位濃度φ點既不夠多、系統性也不夠足等問題.針對該問題，本文用MC方法［15］，對2D-RSIM體系的比熱隨溫度變化，特別是所表現的相變行為，隨φ的演化，進行了細致的計算機模擬.

1 二維隨機格點伊辛模型與蒙特卡羅方法介紹

2D-RSIM模型的哈密頓量H為

其中，σi,j表示二維點陣中第i行第j列格點上的自旋，σi,j=1或-1，表示自旋的狀態；J為最近鄰自旋之間的相互作用能常數；φ為自旋空位濃度(0 ＜φ≤1)；r是0到1之間的隨機數；ri,j為隨機函數，當r＜φ時=0；r≥φ時=1；n為沿著晶軸方向的自旋數目；n→∞表示熱力學極限（thermodynamic limit）；μ0為真空磁導率；μ為自旋磁矩.對伊辛自旋體系，H就是體系中任意一個自旋取1或-1的自旋構型（spin configuration）的能量.

另外，上述自旋體系處于可以交換能量的溫度為T的熱?。╤eat bath）中，即自旋體系和熱浴共同構成一個正則系綜（canonical ensemble）.

當模型體系處于熱平衡時，基于Boltzmann-Gibbs 統計得，體系中單個自旋的平均內能U、平均比熱C為

將方程（3A）（3B）（3C）代入方程（2B），原則上可以嚴格計算2D-RSIM的熱平衡態的C，但是如上文所述，迄今為止仍然未獲得其精確解.

MC模擬是一種以概率和統計理論為基礎，用計算機實現統計抽樣，以獲得問題近似解的方法.下文具體介紹利用MC方法，對2D-RSIM體系中，晶體點陣、點陣上自旋隨機分布、體系中自旋翻轉、體系的熱平衡自旋構型、體系的C等模擬過程.

因為熱力學實際系統包含巨大的粒子數，所以往往可以看作是一個無限大的系統.但是當進行模擬時，因為計算機內存和運行時間的限制，只能選擇有限體系.本文選擇模擬2D-RSIM體系中所包含的自旋數目n2(1-φ)=2002.由于有限體系必然存在邊界，本文選用的是自由性邊界條件，即邊界上自旋與模擬體系之外不存在相互作用.

二維正方點陣的模擬：生成晶格常數為a、原胞數目為n×n的二維正方點陣，如圖1所示.

圖1 二維隨機格點伊辛模型（2D-RSIM）示意圖（空位濃度φ=0.3），其中藍色實線表示正方點陣，紅色圓點表示自旋

點陣上隨機自旋的模擬：對應每個格點位置，生成0至1之間的隨機數r，若r≥φ，則在格點上模擬放置自旋；若r＜φ，表示沒有自旋（如圖1）.

自旋翻轉的模擬：本文選擇自旋翻轉的Glauber 動力學［5，6，16］，即自旋構型的變化由單個自旋的翻轉導致，并且單個自旋翻轉概率（自旋每次嘗試翻轉過程中，實現由σi,j向-σi,j翻轉的概率）為

其中，Fi,j為σi,j的局域場：

由于不同文獻中對Markov過程［15］的一個MC步（MC Step）的表述不盡相同，本文選擇的一個MC步過程為：對圖1所示的2D-RSIM體系中每個自旋，依次使用自旋翻轉的Glauber動力學（方程4），即利用MATLAB的隨機數函數產生0到1之間的隨機數r;當p＞r時，自旋發生偏轉；當p≤r時，自旋不發生偏轉.

熱平衡自旋構型的模擬：本文選擇自旋全部向上的自旋構型為模擬的初始構型，即σi,j=1.體系每經過一個MC步后，一般會得到一個新的自旋構型，對第k步體系的一個確定自旋構型的能量為

因為本文模擬的系綜是正則系綜，所以選擇隨MC步數變化的Hk作為熱平衡的判據參量，即一定的MC步數內Hk的平均值，基本不隨MC 步數變化，就認為模擬體系趨于溫度為T的熱平衡態.如圖2 所示，當φ=0.0 時，對T=2.50J=1.11Tc（Tc將在下文討論），k＞102步以后，Hk的平均值基本不隨k變化，即可以認為體系趨于熱平衡態；對T=2.26J=1.00Tc，當k＞103步以后，體系趨于熱平衡態；對T=2.03J=0.90Tc，T=1.83J=0.81Tc，T=1.42J=0.63Tc，當k＞20步以后，體系趨于熱平衡態.

圖2 系列溫度（T）下，蒙特卡羅模擬的二維隨機格點伊辛模型（2D-RSIM）（自旋空位濃度φ=0.0）的能量（Hk）隨蒙特卡羅步數的變化結果

圖3的相應模擬結果表明，對不同的φ，當T=0.81Tc時，k＞103以后，體系基本都趨于熱平衡態.

圖3 特定約化溫度（T Tc=0.81，Tc為相變溫度）下，蒙特卡羅模擬的系列自旋空位濃度φ的二維隨機格點伊辛模型（2D-RSIM）的能量（Hk）隨蒙特卡羅步數的變化結果

由于本文重點模擬2D-RSIM的相變特征，即研究溫區主要在Tc附近，k＞103步后認為體系基本趨于熱平衡是可行的.但是，考慮到不同φ以及不同T下，體系到達熱平衡所需的MC步數不同，步數越大模擬結果會越好.因此本文選擇的到達熱平衡的MC步數（kE）為50 000步，總的MC步數（kF）為150 000.

2 模擬結果與分析討論

圖4為系列自旋空位濃度φ、自旋數目為n×n的2D-RSIM中，單個自旋的平均比熱C隨溫度T變化的計算機模擬結果.當φ=0.0時，2D-RSIM簡化為2D-伊辛模型（2D-IM）.由于2D-IM的C的精確解已經由昂薩格（Onsager）給出［3，17］，這里作為驗證本文模擬結果的依據.首先，C的精確解表明，2D-IM的相變溫度Tc=2.27J，本文的模擬結果Tc=2.26J（圖4a），比精確解低0.44%，微小的偏差源于本文模擬的是n×n=200×200的有限體系的邊界效應；其次，C的精確解還表明，其在Tc附近對數發散，即C∝-ln| |T-Tc，除去極為靠近Tc點C值較小外，本文的模擬結果也給出了基本相同的結果（圖4a）.因此，本文對2D-IM的計算機MC模擬結果（方程4-6）與精確解結果（方程2-3）偏差較小，可以預期本文對2D-RSIM的模擬結果也具有較高的可信度.

圖4 系列自旋空位濃度φ、自旋數目為n× n的2D-RSIM中，單個自旋的平均比熱C隨溫度T變化的計算機模擬所得結果.其中，黑線：φ=0.0，n →∞的精確解；紅色圓點：φ=0.0，n×n=200×200；藍色圓點：φ=0.1，n×n=211×211；洋紅色圓點：φ=0.2，n×n=224×224；海軍藍色圓點：φ=0.3，n× n=239×239；紫色圓點：φ=0.4，n× n=258×258；橄欖綠色圓點：φ=0.5，n× n=283×283；深青色圓點：φ=0.6，n× n=316×316；皇家藍圓點：φ=0.7，n× n=365×365；橙色圓點：φ=0.8，n× n=447×447；紫羅蘭圓點：φ=0.9，n× n=633×633

圖4的結果還表明，對所有的φ值，C隨著T變化明顯存在一個峰，并且隨著φ變大，峰形由尖銳變得越來越圓滑，并且在不斷變寬.雖然按照早期的Ehrenfest相變分類法［18］，圓滑的比熱峰表示體系中不存在相變；但是隨著研究的深入，盡管仍然存在爭議［19，20］，現在一般認為體系中發生了彌散相變（diffuse phase transi‐tion）［6］或者初始相變（incipient phase transition）［19］.

本文作者認為2D-RSIM中發生了彌散相變，并將C的峰值對應的溫度稱為相變溫度Tc（圖4a）.由圖5可見，Tc隨φ的增大總體上單調降低，細致分析則可以發現存在3 個Tc隨φ線性變化的區域（I 區、II 區、III區），其中I區下降斜率最大，III區下降斜率最小.上述Tc隨φ變化的3個線性區域是本文首次發現，對應的微觀機制仍需深入研究.

圖5 計算機模擬所得的二維隨機格點伊辛模型（2D-RSIM）的相變溫度Tc隨自旋空位濃度φ的變化結果