伊辛模型铁磁相变的研究

姚贵斌 物基一班 2014301020066

摘要:利用蒙特卡罗算法，对二维伊辛模型下的铁磁相变作了模拟和研究，对其磁化强度、平均能量、平均热容、对外场的敏感性作了模拟。
关键词： 铁磁相变 伊辛模型 蒙特卡罗算法
代码

背景介绍

伊辛模型

伊辛模型是模拟铁磁性物质的结构，解释从铁磁性到顺磁性这类相变现象的一种粗略的模型。它的优点在于，用统计物理方法，对二维情形求得了数学上严格的解。这就使得铁磁性物质相变的大致特征，获得了理论上的描述。
这个模型所研究的系统是由N个阵点排列成n维周期性点阵，这里n=1,2,3。点阵的几何结构可以是立方的或六角形的,每个阵点上都赋予一个取值+1或－1的自旋变数i，如果i=+1,即第N个阵点的自旋向上；如i=-1，即第N个阵点的自旋向下并且认为只是最近邻的自旋之间有相互作用。点阵的位形用一组自旋变数（i=1，2，…N,）来确定。

图为一维伊辛模型示意图。箭头向上表示自旋为+1，箭头向下表示自旋为-1.模型的总能量为

$$E=-J\Sigma_{<ij>}s_is_j-\mu H\Sigma_is_i$$
式中：
（1）$s_i$表示第i个自旋的取值（只能取+1或-1），$J>0$为常数，与相互作用强度有关；
（2）$\Sigma_{<ij>}s_is_j$表示只对相邻最近的两自旋求和；
（3）$H$为磁场强度，$\mu$为每个自旋的磁矩。
伊辛模型采用的是正则系综，即体系粒子数不变并与一恒温热源达到平衡。对正则系综，其处于某一状态$\alpha$的概率为
$$P_\alpha=\frac{1}{Z}e^{-E_\alpha/k_BT}$$
其中$Z=\Sigma_\alpha e^{-E_\alpha/k_BT}$为配分函数，$E_\alpha$为状态$\alpha$的能量，$T$为热源温度，$k_B$为玻尔兹曼常数。

二维情况

用伊辛模型描述二维体系的翻转过程。单个自旋的能量：

$$E=-(J\Sigma_{<ij>}s_j)s_i-\mu Hs_i$$
当其翻转（即自旋与原来相反）后，能量变为：
$$E'=-(J\Sigma_{<ij>}s_j)(-s_i)-\mu H(-s_i)=-E$$
故翻转所需能量为$E_{flip}=-2E$
按照体系能量总是最低的原则，当$E_{flip}<0$时，必定翻转；当$E_{flip}>0$时，其反转的概率应为
$$P_{flip}\sim e^{-E_{flip}/k_BT}$$
此过程可以用蒙特卡罗算法实现。
二位方格情况下，设无外磁场，即$H=0$，记方格上第i列第j行的自旋为$s(i,j)$,则对内部自旋，与之相距最近的有四个。
$$E=-[s(i-1,j)+s(i+1,j)+s(i,j-1)+s(i,j+1)]s_i$$
对于边界上的自旋，采用周期性边界条件:对最左边一列，令其左边为最右边一列；对最下边一行，令其下面一行为最上边一行；其余同理。

对相变的模拟

采用不同的单位制。以$J$作为能量的单位，以$J/k_B$作为温度的单位。
四种温度下，$10\times 10$二维方格自旋体系的总磁化强度

$$M=\Sigma_i s_i$$
随时间变化的情况，令初始时所有的自旋均向上。具体变化情况如下图所示。

性质随温度的变化

温度从0-5单位内，平均磁化强度随温度的变化

温度在2.2单位附近时，总磁化强度变化很大以致看起来不“连续”。即临界温度大致为2.2$J/k_BT.$

同样范围内每个自旋的平均能量随温度的变化

随着温度升高，每个自旋的平均能量一直增加。增加速度先增加后降低，此时的临界温度大致为2.4$J/k_BT.$温度达到5$J/k_BT$后，平均能量仍小于0.

单个自旋平均热容量随温度的变化

右图中取样点数与左图相同，只是间隔变小。
可见温度趋于零时的热容也趋于零；在临界温度前，随着温度升高，热容升高且升高加快，达到一很大的极值后再降低且降低减缓。
临界温度在2.3$J/k_BT$与2.4$J/k_BT$之间。

$\chi=\frac{dM}{dH}$随温度的变化(习题8.6)

上面三图中取样点数相同。
由左上图可见$\chi$可以达到很大（在临界温度时），且几乎是突然地，“陡峭”地增加的。到最大值以后再逐渐减小。
临界温度大致为2.2$J/k_BT$，略小于2.2$J/k_BT$。

不同边界条件

上述模拟均采用周期性边界条件。将有限边界条件（即边界之外视为0）与周期边界条件作比较（习题8.10）：

其中图例PBC指周期边界条件（Periodic Boundary Condition），LBC指有限边界条件（Limited Boundary Condition）左图与右图点列来自两次分开的模拟
可见有限边界条件对应的临界温度比周期边界条件对应的临界温度要小，也更小于理论值2.27$J/k_BT$.从总体趋势来看，LBC的图像相对PBC为向左平移了，有限边界条件会降低体系的“耐热能力”。

总结

实验中某些由物理量的模拟结果得到的临界温度与理论值2.27仍有差距。可能原因有：
1. 蒙特卡罗算法的性质就决定了它得到的是近似值不是精确结果（取样点越多越接近精确值），因此误差无法避免。
2. 程序中在先进行了一定多的时间以使其在计算开始时是平衡状态。可能由于温度在接近临界温度时达到平衡的过程变慢，导致实际计算开始时不是平衡状态而带来误差。
3. 误差不影响总体趋势。

参考资料

1. 计算物理第二版
2. 百度百科词条伊辛模型