蒙特卡洛方法模擬二維伊辛模型的相變分析

詹明浩 張威威 梁炳清 余寶明 蔣藍(lán)毅

（1.北部灣大學(xué)，廣西 欽州 535011；2.常州大學(xué)，江蘇 常州 213000）

0 前言

在研究系統(tǒng)處于平衡態(tài)時(shí)的熱力學(xué)性質(zhì)時(shí)，通常通過對(duì)配分函數(shù)求導(dǎo)得到系統(tǒng)的熱力學(xué)函數(shù)[1]。在物質(zhì)發(fā)生相變時(shí)，相變點(diǎn)的某些熱力學(xué)量會(huì)發(fā)生突變或者出現(xiàn)無窮尖峰，必須通過建立體現(xiàn)系統(tǒng)核心特性的模型，得到其在相變點(diǎn)的宏觀特性，才能從單一的配分函數(shù)表達(dá)式中同時(shí)描述相和相的轉(zhuǎn)變。從20世紀(jì)30年代中葉開始，統(tǒng)計(jì)物理學(xué)對(duì)統(tǒng)計(jì)模型的大量研究已經(jīng)形成了一個(gè)專門的領(lǐng)域，誕生了許多模型，在解釋鐵磁性的問題上，伊辛模型體現(xiàn)了它的優(yōu)勢(shì)[2]。伊辛模型盡管很成熟、非常經(jīng)典，但是現(xiàn)在也需要不斷地對(duì)其進(jìn)行驗(yàn)證。而對(duì)于蒙特卡洛方法，不斷將它與新的科學(xué)技術(shù)結(jié)合，使之獲得新的生命力并得到發(fā)展，嘗試性地取得一些微小的進(jìn)步是非常具有意義的。綜上所述，用蒙特卡洛方法（MC）模擬二維伊辛（lsing）模型的相變過程就顯得尤為重要。

1 蒙特卡洛方法與伊辛模型

1.1 伊辛模型

考慮1個(gè)一維的鐵磁體，定義這個(gè)自旋粒子，用（↑↓）分別表示自旋的方向。將該粒子排列在規(guī)則的晶格上。20世紀(jì)20年代楞次（Wilhelm Lenz）教授就是由此創(chuàng)立了伊辛模型，一維伊辛模型的示意圖如圖1所示。

圖1 一維伊辛模型示意圖

圖1中每一個(gè)黑點(diǎn)代表一個(gè)自旋，它們之間存在最近鄰的相互作用。一維伊辛模型給出的結(jié)果如下：一維的伊辛模型只有一個(gè)相，不存在相變。該結(jié)果正是1925年楞次教授的博士生伊辛的成果，因此該模型也被稱楞次 -伊辛模型。

1944年該理論出現(xiàn)了一個(gè)重要進(jìn)展，美國(guó)科學(xué)家昂薩格教授解出了2D伊辛模型的解析。

圖2（a）是8×8的1個(gè)二維伊辛模型晶格，黑點(diǎn)代表自旋，圖2（b）是以中心自旋n為例的相互作用的示意圖，n1，2，3，4是不同位置的最近鄰自旋，其相互作用關(guān)系用{i，j}來表示。在這里，筆者給出最初昂薩格的解[3]，即Tc=2.269。

圖2 二維伊辛模型示意圖

1.2 蒙特卡洛方法的基本思想

伊辛模型系統(tǒng)各格點(diǎn)上的自旋取值抽樣采用了Metroplis 方法[4]。這對(duì)于伊辛模型來說雖然不是最高效的采樣方法，但卻是最易于實(shí)現(xiàn)的方法。該文采用MATLAB軟件在20×20的晶格中進(jìn)行模擬。蒙卡方法的基本思想如下：對(duì)于某一問題的隨機(jī)模擬，通過問題的解來設(shè)定參數(shù)，根據(jù)現(xiàn)實(shí)情況來建模。以多次隨機(jī)模擬來求得問題的概率。設(shè)隨機(jī)變量的數(shù)學(xué)期望是E（ξ），對(duì)ξ進(jìn)行N次模擬得到各個(gè)ξ值，可得到其算術(shù)平均值ξN如公式（1）所示。

式中：ξ為ξN的計(jì)算結(jié)果。

當(dāng)N→∞時(shí)，。

首先，產(chǎn)生隨機(jī)數(shù)r1，r2，且它們范圍是0～1的數(shù)。

其次，將r1和r2對(duì)應(yīng)于直角坐標(biāo)系上的x、y值，如果r2≤e-r1，那么就可多設(shè)一個(gè)未知數(shù)v來表示觀察的次數(shù)，符合統(tǒng)計(jì)規(guī)則，v加1，N也加1；如果r2＞e-r1，那么不符合統(tǒng)計(jì)規(guī)則，則v加0，N加1。

再次，多次重復(fù)上述步驟。

最后，得出結(jié)果。

2 二維伊辛模型的蒙特卡洛模擬

比熱容、能量與溫度T的關(guān)系如公式（2）所示。

圖3為2D伊辛模型的比熱容Cv與平均能量E隨溫度的變化，上部分圖的熱容量Cv與下部分圖的平均能量E相對(duì)應(yīng)。比熱可以作為判斷相變的標(biāo)準(zhǔn)，且其判斷結(jié)果比較精確，它在相變時(shí)作為熱力學(xué)量會(huì)出現(xiàn)尖峰，E表示單個(gè)格點(diǎn)的平均能量，當(dāng)溫度趨于0時(shí)，所有格點(diǎn)同向。

圖3 比熱容Cv與平均能量E隨溫度T的變化

當(dāng)溫度趨于無窮時(shí)，格點(diǎn)方向隨機(jī)，某格點(diǎn)四周平均有2個(gè)同向和2個(gè)方向，E=0。通過模擬得到圖4。圖4中模擬的比熱在臨界溫度時(shí)會(huì)突然增大，需要吸收極大的能量，這也很符合相變的特點(diǎn)，在這之后比熱逐漸下降。可以看出在溫度低于相變臨界溫度的情況下，熱容量隨溫度的升高而增大；在高于相變臨界溫度時(shí)，隨溫度的升高而減小。當(dāng)溫度為2.273 K時(shí)，熱容圖像由升到降出現(xiàn)相變，這與相變臨界溫度理論值2.269 K吻合，該結(jié)果對(duì)二維伊辛模型的進(jìn)一步研究具有深遠(yuǎn)的意義，同時(shí)也是下一步研究的基礎(chǔ)。

3 磁場(chǎng)及耦合強(qiáng)度對(duì)二維伊辛模型相變點(diǎn)臨界溫度的影響

3.1 結(jié)果與討論

3.1.1 磁場(chǎng)對(duì)伊辛模型的影響

如圖5所示，在不同磁場(chǎng)作用下相變點(diǎn)存在以下變化，見表1。

表1 2D伊辛模型在不同磁場(chǎng)B作用下系統(tǒng)的相變臨界點(diǎn)

在之前的基礎(chǔ)上，通過計(jì)算機(jī)模擬了耦合強(qiáng)度J=1，對(duì)磁場(chǎng)B=0、B=0.3、B=0.5、B=0.7以及B=0.9進(jìn)行計(jì)算，用磁感應(yīng)強(qiáng)度B來表示磁場(chǎng)強(qiáng)度，磁場(chǎng)強(qiáng)度H的單位為A/m，換算關(guān)系為B=4π×10-7N/A2·H，得到圖像如圖5所示的關(guān)系圖，上部分圖的熱容量Cv與下部分圖的平均能量E相對(duì)應(yīng)。

從圖5可以看出，磁場(chǎng)對(duì)伊辛模型有作用，使系統(tǒng)相變溫度出現(xiàn)了延遲，磁場(chǎng)絕對(duì)值越大，相變臨界溫度就越高。當(dāng)存在外磁場(chǎng)時(shí)，相變臨界點(diǎn)隨著外磁場(chǎng)的增大向高溫方向移動(dòng)，這是由于外磁場(chǎng)使自旋趨于有序排列，而熱運(yùn)動(dòng)使自旋趨于無序排列；隨著外場(chǎng)的增大，磁場(chǎng)的作用逐漸增強(qiáng)，無序區(qū)間向高溫方向壓縮，而有序區(qū)間則向高溫方向擴(kuò)大，導(dǎo)致臨界溫度增高；因此，當(dāng)外磁場(chǎng)足夠強(qiáng)時(shí)，系統(tǒng)幾乎在整個(gè)溫度空間內(nèi)處于有序狀態(tài)且會(huì)出現(xiàn)相變臨界的現(xiàn)象。

3.1.2 耦合強(qiáng)度對(duì)伊辛模型的影響

不同耦合強(qiáng)度下系統(tǒng)熱容量和平均能量的關(guān)系如圖4所示。

圖4 在不同耦合強(qiáng)度J作用下系統(tǒng)比熱容Cv與平均能量E隨溫度的變化

由圖5可知，在不同耦合強(qiáng)度作用下相變點(diǎn)存在以下變化，見表2。

圖5 不同磁場(chǎng)B作用下系統(tǒng)比熱容Cv與平均能量E隨溫度T的變化

表2 2D伊辛模型在不同耦合強(qiáng)度作用下系統(tǒng)的相變臨界點(diǎn)

這里討論了耦合強(qiáng)度對(duì)伊辛模型的影響，取外磁場(chǎng)B=0，分別取耦合強(qiáng)度J=0.8、J=0.9、J=1、J=1.25以及J=1.5進(jìn)行計(jì)算，計(jì)算結(jié)果如圖5所示，上部分圖的熱容量Cv與下部分圖的平均能量E相對(duì)應(yīng)。

從圖5可以看出，在不同耦合強(qiáng)度下，系統(tǒng)的能量總是呈上升的趨勢(shì)，隨著耦合強(qiáng)度的增大，能量曲線向右平移，相變臨界值增加。得出結(jié)果為2.381，與預(yù)期相變點(diǎn)2.269存在誤差，之所以存在誤差，可能是由以下2個(gè)原因?qū)е碌模?）實(shí)驗(yàn)本身誤差。2）相鄰局域的溫度不同，相互之間有耦合作用。

耦合強(qiáng)度在一定波動(dòng)范圍內(nèi)對(duì)二維伊辛模型的相變存在影響，耦合強(qiáng)度越大，相變點(diǎn)臨界溫度也越高。如果耦合強(qiáng)度變化過大，將會(huì)使模擬失敗，這違反了基本自旋相互作用的假設(shè)。

3.2 模擬結(jié)論

在模擬過程中還發(fā)現(xiàn)了如下問題：當(dāng)程序接近相變臨界點(diǎn)時(shí)，Metroplis方法會(huì)收斂得很慢，單個(gè)自旋的翻轉(zhuǎn)統(tǒng)計(jì)會(huì)導(dǎo)致如果系統(tǒng)準(zhǔn)備進(jìn)入完全自旋向上或自旋向下時(shí)，突然接受了某個(gè)粒子的翻轉(zhuǎn)，那么又會(huì)從頭進(jìn)行一輪新的統(tǒng)計(jì)，這就會(huì)導(dǎo)致在亞穩(wěn)態(tài)停留的時(shí)間過長(zhǎng)，甚至出現(xiàn)無法收斂到穩(wěn)態(tài)而一直停留在亞穩(wěn)態(tài)的情況。希望在以后條件允許的情況下可以對(duì)算法的這個(gè)問題進(jìn)行改進(jìn)。

4 結(jié)語

在用MC模擬二維伊辛模型的相變的問題上，前人有非常豐富的經(jīng)驗(yàn)，該文通過查閱文獻(xiàn)資料，重現(xiàn)了MC模擬二維伊辛模型的相變過程，并對(duì)該過程進(jìn)行模擬、探索、對(duì)比與分析，豐富了該理論的模擬過程，從而使該理論得到發(fā)展。

同時(shí)，該文還重新模擬這個(gè)古老又經(jīng)典的模型，是對(duì)真理的不斷驗(yàn)證，還可以為其他科技工作者給予參考。

進(jìn)一步研究耦合強(qiáng)度與溫度關(guān)系中磁化率和磁化強(qiáng)度的規(guī)律對(duì)該模型具有重要意義；運(yùn)用平均場(chǎng)理論計(jì)算合理的磁場(chǎng)梯度，在適當(dāng)?shù)鸟R爾可夫鏈長(zhǎng)和步數(shù)下進(jìn)行模擬；對(duì)采樣算法進(jìn)行改進(jìn)，優(yōu)化相變臨界點(diǎn)的收斂速度。