Matlab 蒙特卡洛模拟⼆维伊⾟模型相变过程
⼀、什么是伊⾟模型
伊⾟（Ising）模型是描述磁系统相变最简单的模型，但模型⾥⾃旋之间的相互作⽤赋予了它奇妙的特性，最有趣的就是对称性破缺。这⼀模型可以被推⼴⽤于研究连续的量⼦相变、基本粒⼦超弦理论、动⼒学临界⾏为等，甚⾄被认为可以描述深林⽕灾、交通拥堵、舆论传播等
社会经济现象。
如图，每个格点的⽅向只有向上或向下两者状态，但临近的⾃旋之间有相互作⽤，⽽且点阵可以是⼀维、⼆维、三维、甚⾄更⾼维，这两个特点让伊⾟模型的严格求解成为了世纪难题。为了定量描述这个系统的能量，我们假设第个格点的⾃旋为，只能取+1或-1，如果相邻两个格点同⽅向，则它们相互作⽤的能量更⼩，设为，如果为反⽅向，则为，称为耦合系数，通常为正值，代表铁磁系统，如果为负值，则代表反铁磁系统。如果外磁场的强度为B，格点的⾃旋磁矩为，那么可以写出这个体系的哈密顿量：
⾃发对称性破缺
我们先定性地了解⼀下这个系统的性质，令外磁场零，当温度时，体系为了保持能量最低，所⽤的格点趋向于同⽅向，系统整体要么向下，要么向上，呈现强磁性。当温度时，系统热运动占主导地位，格点⽅向呈现随机性，系统整体不带磁性，从上或从下观察体系，呈现出对称性，或者说⽆法通过系统磁性区分上或下。现在再考虑，当温度T从逐渐降温，那么系统必定存在某个温度，⾼于此温度时系统⽆磁性，低于此温度时，系统磁性逐渐加强。这个温度就是临界温度，也是相变点，系统从对称磁体转变为⾮对称磁体，⽽这就是
对称性破缺，因为这种破缺不是外界扰动（如外加磁场）引起的，⽽是由内部的关联作⽤⼒造成的，所以称之为⾃发的对称性破缺。上图就是对⾃发对称性破缺的定性描述，当逐渐降温到时，系统磁性开始出现分化，要么向下要么向上，最终平均的磁化强度趋向于+1或-1。
我们感兴趣的问题主要有两个：第⼀，不同维度、不同分布的格点，其临界温度是多少；第⼆，附近，随温度以什么样的幂指数趋近于：
i s i s i −J J J J μH =−J s s −<ij >∑i j μB s i =1∑N
i T →0T →∞∞T c T c s T c T c s T αT c ∼s (T −c T )α
统计物理的思路
我们先考虑简单的9个格点的例⼦，实际格点数的量级为
。
假设格点耦合强度，那么这个9格点体系的能量为：
第⼀个格点有两个相邻格点，右边的与其同向，耦合能为-1，下⾯的与其反向，耦合能为1；第⼆个格点有三个相邻格点，左和下与其同向，耦合能为-2，右边与其反向，耦合能为1。类似的可以计算其余格点的耦合能。最后除以2是因为每个相互作⽤重复计算了⼀次。⽽这⼀特定分布以某概率P出现：
对于9个格点的体系，总共有种分布，每种分布的出现概率于其总能量有关，所以分布概率满⾜归⼀化条件。给定不同温度，我们可以计算出不同温度下平均磁化率的数学期望，得到的曲线。
但是对于粒⼦为量级的格点，可能的分布有种，根本不可能统计出结果。
⼀维的情况可以通过数学上的处理，最终可以提出N，并得到，也就是说⼀维的伊⾟模型不会有⾃发
的对称性破缺，这是因为⼀维的格点只有两个相邻格点，相互作⽤太弱，不⾜以对抗热运动，始终表现为整体0磁化率的对称状态。
⼆维的情况下，如果⽤平均场近似的⽅法（具体可以参考林宗涵的热⼒学与统计物理），基本思想是将相互作⽤的耦合能转化为外磁场强度，这就可以⽤近独⽴的模型来计算配分函数，进⽽得到所有的统计量，获得的临界温度为 ，平均磁化率：1944年，昂萨格推导出了⼆维伊⾟模型的严格解，临界温度，平均磁化率：
⼆、⼆维伊⾟模型的精确解
matlab求导平均磁化率 这⾥只列出⼆维伊⾟模型精确解的结论，推导过于复杂，定义，则平均磁化率：令可以解得：，令，⼩量泰勒展开化简可以得到：
当时，容易得到结果为1，表⾯所有的格点同向，平均磁化率为1。
1029J =1E =[(−1+1)+(−2+1)+(2)+(2−1)+(−3+1)+(−1+2)+(−1+1)+(−1+2)+2]/2=2
P ∝e −E /kT
2=9512s −s T 102921029
T =c 0T =c k 2J ∼s (T −c T )        (T →1/2T )
c −T =c k 2.269J ∼s (T −c T )        (T →1/8T )
c −s
β≡kT 1
=s {0,[1−sinh(2βJ )],−41/8T >T c
T ≤T c
sinh(2βJ )=0T =c ≈k ln(1+)22J
k 2.269J T =T −c δT =s [4⋅cosh()⋅]=kT c 2J kT c 2J T c δT 1/8  1.224[1−]∼T c T
1/8(T −c T )1/8
T →c 0
假设耦合强度等于1开尔⽂温度下的热运动能量，即，做出
曲线如下：
平均能量和⽐热
J J =k −s T
另⼀个能够判断相变的参数是⽐热，这⾥的表⽰单个格点的平均能量，定性来看，当温度趋于0时，所有格点同向，，当温度趋于⽆穷时，格点⽅向随机，某格点四周平均有两个同向和两个⽅向，。其曲线如下(令
)：
这⾥的⽐热在临界温度时会突然增⼤，表⾯临界温度附近变化有个⼩温度变化，需要吸收极⼤能量，这也很符合相变的特点。具体的计算公式和matlab代码详见附录A。
虽然⼤体了解了相变的过程，以及理论上的精确解，我们能否通过实验的⽅式，验证这⼀结论呢？借助计算机模拟这⼀过程来验证结果呢？
三、⼆维伊⾟模型模拟
因为不可能遍历所有的格点组合，我们只能利⽤采样的⽅式去计算平均能量，采样的条件应该是体系在某个温度下已经平衡。 计算机模拟的基本思想是，⾸先随机给定⼀种分布，在特定温度下，让体系趋向平衡，再在这个平衡体系中采样求平均。
假设体系有的格点，初始时同⼀分布，相当于温度很低；
我们设定⼀个希望“加热”到的平衡温度，接下来是模拟最关键的地⽅，如何改变格点的分布以趋向于设定温度？
1、我们先任意选择⼀个格点，计算改变这个格点的能量变化，因为体系出现的概率正⽐于，那么格点变化前后两个体系出现的概率⽐为，或者说格点改变的概率与不改变的概率⽐为：，那么格点需要改变的概率为，因此我们产⽣⼀个(0,1)概率随机数，如果它⼩于则选择改变格点；
2、重复1的步骤，直到体系相对平稳，这个过程称为马尔可夫链，之后在平稳的体系下采样若⼲次并做统计平均，获得平均能量，平均磁化率，以及描述能量⽅差的⽐热。
3、设定另⼀个希望考察的温度，重复1、2步骤，获得该温度下的统计参数，并和理论值⽐较。
我们同样假设，选取格点数为。临界温度点附近，马尔可夫链长取5万次，采样数为25万次；其他温度点马尔可夫链长1万次，采样数为5万次。这是因为临界温度附近的涨落很⼤，需要更长的时间趋向平衡，需要更多的统计样本获得较准确平均值。详细的代码及解释可以参看附录B。
模拟平均磁化率 C =v dT dE
E =E −4/2=−2=E 0J =k 20×20T 0T 0ΔE e −E /kT e −ΔE /kT 1:e −ΔE /kT e −ΔE /kT e −ΔE /kT E s C v J =k 20×20−s T
可以看出平均磁化率在临界温度附近很不稳定，这是因为临界相变时涨落很⼤的缘故，⾼温时的磁化
率不是严格的零，可能与格点数少和马尔可夫链较短有关系，如何确定呢？通过平均磁化率求⽐较困难，⼀般是通过⽐热发散的位置确定，参考下⼀部分。选取T=1.7~2.2的16个数据点，拟合曲线：
得到，这和理论值相当接近。
模拟平均能量 T c T c C v T ≈c    2.3ln ∼s αln(T −c T )
α≈0.12,R =20.891/8=0.125−E T