一、什么是伊辛模型
伊辛(Ising)模型是描述磁系统相变最简单的模型,但模型里自旋之间的相互作用赋予了它奇妙的特性,最有趣的就是对称性破缺。这一模型可以被推广用于研究连续的量子相变、基本粒子超弦理论、动力学临界行为等,甚至被认为可以描述深林火灾、交通拥堵、舆论传播等社会经济现象。
如图,每个格点的方向只有向上或向下两者状态,但临近的自旋之间有相互作用,而且点阵可以是一维、二维、三维、甚至更高维,这两个特点让伊辛模型的严格求解成为了世纪难题。为了定量描述这个系统的能量,我们假设第
i
i
i个格点的自旋为
s
i
s_i
si,
s
i
s_i
si只能取+1或-1,如果相邻两个格点同方向,则它们相互作用的能量更小,设为
−
J
-J
−J,如果为反方向,则为
J
J
J,
J
J
J称为耦合系数,通常为正值,代表铁磁系统,如果
J
J
J为负值,则代表反铁磁系统。如果外磁场的强度为B,格点的自旋磁矩为
μ
\mu
μ,那么可以写出这个体系的哈密顿量:
H
=
−
J
∑
<
i
j
>
s
i
s
j
−
μ
B
∑
i
=
1
N
s
i
H=-J\sum_{<ij>}s_i s_j-\mu B \sum_{i=1}^N s_i
H=−J<ij>∑sisj−μBi=1∑Nsi
自发对称性破缺
我们先定性地了解一下这个系统的性质,令外磁场零,当温度
T
→
0
T\rightarrow 0
T→0时,体系为了保持能量最低,所用的格点趋向于同方向,系统整体要么向下,要么向上,呈现强磁性。当温度
T
→
∞
T\rightarrow \infty
T→∞时,系统热运动占主导地位,格点方向呈现随机性,系统整体不带磁性,从上或从下观察体系,呈现出对称性,或者说无法通过系统磁性区分上或下。现在再考虑,当温度T从
∞
\infty
∞逐渐降温,那么系统必定存在某个温度
T
c
T_c
Tc,高于此温度时系统无磁性,低于此温度时,系统磁性逐渐加强。这个温度就是临界温度,也是相变点,系统从对称磁体转变为非对称磁体,而这就是对称性破缺,因为这种破缺不是外界扰动(如外加磁场)引起的,而是由内部的关联作用力造成的,所以称之为自发的对称性破缺。
上图就是对自发对称性破缺的定性描述,当逐渐降温到
T
c
T_c
Tc时,系统磁性开始出现分化,要么向下要么向上,最终平均的磁化强度
s
‾
\overline{s}
s趋向于+1或-1。
我们感兴趣的问题主要有两个:第一,不同维度、不同分布的格点,其临界温度
T
c
T_c
Tc是多少;第二,
T
c
T_c
Tc附近,
s
‾
\overline{s}
s随温度
T
T
T以什么样的幂指数
α
\alpha
α趋近于
T
c
T_c
Tc:
s
‾
∼
(
T
c
−
T
)
α
\overline{s}\sim(T_c-T)^{\alpha}
s∼(Tc−T)α
统计物理的思路
我们先考虑简单的9个格点的例子,实际格点数的量级为
1
0
29
10^{29}
1029。
假设格点耦合强度
J
=
1
J=1
J=1,那么这个9格点体系的能量为:
E
=
[
(
−
1
+
1
)
+
(
−
2
+
1
)
+
(
2
)
+
(
2
−
1
)
+
(
−
3
+
1
)
+
(
−
1
+
2
)
+
(
−
1
+
1
)
+
(
−
1
+
2
)
+
2
]
/
2
=
2
E=[(-1+1)+(-2+1)+(2)+(2-1)+(-3+1)+(-1+2)+(-1+1)+(-1+2)+2]/2=2
E=[(−1+1)+(−2+1)+(2)+(2−1)+(−3+1)+(−1+2)+(−1+1)+(−1+2)+2]/2=2
第一个格点有两个相邻格点,右边的与其同向,耦合能为-1,下面的与其反向,耦合能为1;第二个格点有三个相邻格点,左和下与其同向,耦合能为-2,右边与其反向,耦合能为1。类似的可以计算其余格点的耦合能。最后除以2是因为每个相互作用重复计算了一次。而这一特定分布以某概率P出现:
P
∝
e
−
E
/
k
T
P\varpropto e^{-E/kT}
P∝e−E/kT
对于9个格点的体系,总共有
2
9
=
512
2^9=512
29=512种分布,每种分布的出现概率于其总能量有关,所以分布概率满足归一化条件。给定不同温度,我们可以计算出不同温度下平均磁化率的数学期望
s
‾
\overline{s}
s,得到
s
‾
−
T
\overline{s}-T
s−T的曲线。
但是对于粒子为 1 0 29 10^{29} 1029量级的格点,可能的分布有 2 1 0 29 2^{10^{29}} 21029种,根本不可能统计出结果。
一维的情况可以通过数学上的处理,最终可以提出N,并得到 T c = 0 T_c=0 Tc=0,也就是说一维的伊辛模型不会有自发的对称性破缺,这是因为一维的格点只有两个相邻格点,相互作用太弱,不足以对抗热运动,始终表现为整体0磁化率的对称状态。
二维的情况下,如果用平均场近似的方法(具体可以参考林宗涵的热力学与统计物理),基本思想是将相互作用的耦合能转化为外磁场强度,这就可以用近独立的模型来计算配分函数,进而得到所有的统计量,获得的临界温度为
T
c
=
2
J
k
T_c= \frac{2J}{k}
Tc=k2J,平均磁化率:
s
‾
∼
(
T
c
−
T
)
1
/
2
(
T
→
T
c
−
)
\overline{s}\sim(T_c-T)^{1/2} ~~~~~~~~(T\rightarrow T_c^-)
s∼(Tc−T)1/2 (T→Tc−)
1944年,昂萨格推导出了二维伊辛模型的严格解,临界温度
T
c
=
2.269
J
k
T_c=\frac {2.269J}{k}
Tc=k2.269J,平均磁化率:
s
‾
∼
(
T
c
−
T
)
1
/
8
(
T
→
T
c
−
)
\overline{s}\sim(T_c-T)^{1/8} ~~~~~~~~~(T\rightarrow T_c^-)
s∼(Tc−T)1/8 (T→Tc−)
二、二维伊辛模型的精确解
平均磁化率 s ‾ \overline s s
这里只列出二维伊辛模型精确解的结论,推导过于复杂,定义
β
≡
1
k
T
\beta\equiv \frac{1}{kT}
β≡kT1,则平均磁化率:
s
‾
=
{
0
,
T
>
T
c
[
1
−
sinh
(
2
β
J
)
−
4
]
1
/
8
,
T
≤
T
c
\overline{s} =\begin{cases} 0,&{T>T_c} \\ [1-\sinh(2\beta J)^{-4}]^{1/8}, &{T\leq T_c}\end{cases}
s={0,[1−sinh(2βJ)−4]1/8,T>TcT≤Tc
令
sinh
(
2
β
J
)
=
0
\sinh(2\beta J)=0
sinh(2βJ)=0可以解得:
T
c
=
2
J
k
ln
(
1
+
2
)
≈
2.269
J
k
T_c=\frac{2J}{k\ln(1+\sqrt{2})}\approx \frac{2.269J}{k}
Tc=kln(1+2)2J≈k2.269J,令
T
=
T
c
−
δ
T
T=T_c-\delta T
T=Tc−δT,小量泰勒展开化简可以得到:
s
‾
=
[
4
⋅
2
J
k
T
c
cosh
(
2
J
k
T
c
)
⋅
δ
T
T
c
]
1
/
8
=
1.224
[
1
−
T
T
c
]
1
/
8
∼
(
T
c
−
T
)
1
/
8
\overline{s}=[4·\frac{2J}{kT_c}\cosh(\frac{2J}{kT_c})·\frac{\delta T}{T_c}]^{1/8}=1.224[1-\frac{T}{T_c}]^{1/8}\sim(T_c-T)^{1/8}
s=[4⋅kTc2Jcosh(kTc2J)⋅TcδT]1/8=1.224[1−TcT]1/8∼(Tc−T)1/8
当
T
c
→
0
T_c \rightarrow 0
Tc→0时,容易得到结果为1,表面所有的格点同向,平均磁化率为1。
假设耦合强度
J
J
J等于1开尔文温度下的热运动能量,即
J
=
k
J=k
J=k,做出
s
‾
−
T
\overline{s}-T
s−T曲线如下:
平均能量和比热
另一个能够判断相变的参数是比热
C
v
=
d
E
‾
d
T
C_v=\frac{d \overline{E}}{dT}
Cv=dTdE,这里的
E
‾
\overline{E}
E表示单个格点的平均能量,定性来看,当温度趋于0时,所有格点同向,
E
‾
=
−
4
/
2
=
−
2
\overline{E}=-4/2=-2
E=−4/2=−2,当温度趋于无穷时,格点方向随机,某格点四周平均有两个同向和两个方向,
E
‾
=
0
\overline{E}=0
E=0。其曲线如下(令
J
=
k
J=k
J=k):
这里的比热在临界温度时会突然增大,表面临界温度附近变化有个小温度变化,需要吸收极大能量,这也很符合相变的特点。具体的计算公式和matlab代码详见附录A。
虽然大体了解了相变的过程,以及理论上的精确解,我们能否通过实验的方式,验证这一结论呢?借助计算机模拟这一过程来验证结果呢?
三、二维伊辛模型模拟
因为不可能遍历所有的格点组合,我们只能利用采样的方式去计算平均能量,采样的条件应该是体系在某个温度下已经平衡。 计算机模拟的基本思想是,首先随机给定一种分布,在特定温度下,让体系趋向平衡,再在这个平衡体系中采样求平均。
- 假设体系有 20 × 20 20\times 20 20×20的格点,初始时同一分布,相当于温度很低;
- 我们设定一个希望“加热”到的平衡温度
T
0
T_0
T0,接下来是模拟最关键的地方,如何改变格点的分布以趋向于设定温度
T
0
T_0
T0?
1、我们先任意选择一个格点,计算改变这个格点的能量变化 Δ E \Delta E ΔE,因为体系出现的概率正比于 e − E / k T e^{-E/kT} e−E/kT,那么格点变化前后两个体系出现的概率比为 e − Δ E / k T e^{-\Delta E/kT} e−ΔE/kT,或者说格点改变的概率与不改变的概率比为: 1 : e − Δ E / k T 1:e^{-\Delta E/kT} 1:e−ΔE/kT,那么格点需要改变的概率为 e − Δ E / k T e^{-\Delta E/kT} e−ΔE/kT,因此我们产生一个(0,1)概率随机数,如果它小于 e − Δ E / k T e^{-\Delta E/kT} e−ΔE/kT则选择改变格点;
2、重复1的步骤,直到体系相对平稳,这个过程称为马尔可夫链,之后在平稳的体系下采样若干次并做统计平均,获得平均能量 E ‾ \overline E E,平均磁化率 s ‾ \overline s s,以及描述能量方差的比热 C v C_v Cv。
3、设定另一个希望考察的温度,重复1、2步骤,获得该温度下的统计参数,并和理论值比较。
我们同样假设 J = k J=k J=k,选取格点数为 20 × 20 20\times 20 20×20。临界温度点附近,马尔可夫链长取5万次,采样数为25万次;其他温度点马尔可夫链长1万次,采样数为5万次。这是因为临界温度附近的涨落很大,需要更长的时间趋向平衡,需要更多的统计样本获得较准确平均值。详细的代码及解释可以参看附录B。
模拟平均磁化率 s ‾ − T \overline{s}-T s−T
可以看出平均磁化率在临界温度附近很不稳定,这是因为临界相变时涨落很大的缘故,高温时的磁化率不是严格的零,可能与格点数少和马尔可夫链较短有关系,如何确定 T c T_c Tc呢?通过平均磁化率求 T c T_c Tc比较困难,一般是通过比热 C v C_v Cv发散的位置确定 T c ≈ 2.3 T_c\approx 2.3 Tc≈2.3,参考下一部分。
选取T=1.7~2.2的16个数据点,拟合曲线:
ln
s
‾
∼
α
ln
(
T
c
−
T
)
\ln\overline{s}\sim \alpha \ln(T_c-T)
lns∼αln(Tc−T)
得到
α
≈
0.12
,
R
2
=
0.89
\alpha \approx 0.12,R^2=0.89
α≈0.12,R2=0.89,这和理论值
1
/
8
=
0.125
1/8=0.125
1/8=0.125相当接近。
模拟平均能量 E ‾ − T \overline{E}-T E−T
在临界温度附近进行了较密集的温度取点,而且加长了马尔可夫链,但是仍能看到较大的涨落,可以通过这个现象来确定临界温度 T c ≈ 2.3 T_c\approx 2.3 Tc≈2.3。
最后,让我们欣赏一下格点从同一分布到临界温度的变化过程吧,为了便于观察,选择40X40的格点,马尔可夫链长为1万:
初始同一分布的格点,逐渐趋向于高温
T
=
5
T=5
T=5下的平衡态,格点最后呈现随机分布:
初始同一分布的格点,温度降低至
T
=
3
T=3
T=3的平衡态,格点最后呈现小块状:
初始无序分布的格点,温度降低到临界温度
T
=
2.3
T=2.3
T=2.3时,格点最后呈现的块状增大。
初始无序分布的格点,温度降低到临界温度以下
T
=
2
T=2
T=2,格点最后呈现的大的块状,说明已经发生了明显的相变。
初始无序分布的格点,温度降低更多至
T
=
1
T=1
T=1,格点越来越趋向于同一分布。
附录
A、平均能量和比热的精确解:(前方高能)
定义
β
≡
1
/
k
T
\beta \equiv 1/kT
β≡1/kT
E
‾
=
−
J
cot
(
2
β
J
)
×
[
1
+
2
π
A
⋅
B
(
λ
)
]
\overline E=-J\cot(2\beta J)\times[1+\frac{2}{\pi}A·B(\lambda)]
E=−Jcot(2βJ)×[1+π2A⋅B(λ)]
{
A
≡
2
tanh
(
2
β
J
)
2
−
1
B
(
λ
)
≡
∫
0
π
/
2
d
ϕ
1
−
λ
2
sin
ϕ
2
λ
≡
2
sinh
(
2
β
J
)
cosh
(
2
β
J
)
2
\begin{cases} A\equiv 2\tanh(2\beta J)^2-1\\ B(\lambda)\equiv \int_0^{\pi/2} \frac{d\phi}{\sqrt{1-\lambda^2\sin\phi ^2}}\\ \lambda\equiv \frac{2\sinh(2\beta J)}{\cosh(2\beta J)^2} \end{cases}
⎩⎪⎪⎨⎪⎪⎧A≡2tanh(2βJ)2−1B(λ)≡∫0π/21−λ2sinϕ2dϕλ≡cosh(2βJ)22sinh(2βJ)
帝国主义都是纸老虎,我们仔细发现,只要确定了温度
T
T
T或
β
\beta
β,那么可以依次确定
λ
,
B
(
λ
)
,
A
,
E
‾
\lambda,B(\lambda),A,\overline E
λ,B(λ),A,E,也就是平均能量和温度是一一对应的,最后通过求导得到比热
C
v
=
d
E
‾
/
d
T
C_v=d \overline E/d T
Cv=dE/dT。
%matlab code
clear;clc;
beta=0.1:0.01:1; %温度从10到1
for i=1:1:size(beta,2);%遍历beta
lambda=2.*sinh(2.*beta(i))./cosh(2.*beta(i)).^2;%计算lambda
phi=linspace(0,pi/2,1000);%求B的积分参数
b=1./sqrt(1-lambda.^2.*sin(phi).^2);
B=trapz(phi,b);%积分B(lambda)
A=2*tanh(2.*beta(i)).^2-1;
e_bar(i)=-coth(2.*beta(i)).*(1+2/pi.*A.*B); %每个格点的平均能量
end
%%
plot(1./beta,e_bar,'k','LineWidth',2);hold on;%E-T曲线
beta1=beta(2:end)/2+beta(1:end-1)/2;
cv=-beta1.^2.*diff(e_bar)./diff(beta); %对能量求导得到比热
plot(1./beta1,cv,'r','LineWidth',2);
B、蒙特卡洛马尔可夫链模拟二维伊辛模型相变过程
具体细节详见代码:
%matlab code
clear;clc;
n=10000; %马尔可夫链长度1万
ns=20; %20*20的格点
beta_mc=(0.1:0.01:0.4); %温度从10到2.5,链长1万,样品长5万
%T_mc=(2.1:0.01:2.4); %第三批模拟温度设定,临界温度附近取点更密集,还要调整n=50000
%beta_mc=1./T_mc;
tic; %计时用,n=10000时,通常需要跑一分多钟
for jj=1:1:size(beta_mc,2)
X=sign(rand(ns,ns)); %所有格点方向一致,相当于从0度开始升温
%马尔可夫链长度为5万次
for j=1:1:n
%随机选取一个格点,行列存储在index[1,2]
index=unidrnd(ns,1,2);
% 利用周期性边界条件,分别计算格点上下左右四个点行列坐标
tmp1=rem(index(1),ns)+1;tmp2=rem(index(1)+1,ns)+1;tmp3=rem(index(1)-1,ns)+1;
tmp4=rem(index(2),ns)+1;tmp5=rem(index(2)+1,ns)+1;tmp6=rem(index(2)-1,ns)+1;
% 计算改变格点方向后的能量变化
cen=X(tmp1,tmp4);right=X(tmp1,tmp5);left=X(tmp1,tmp6);
up= X(tmp2,tmp4);down= X(tmp3,tmp4);
deE=2*cen*(right+left+up+down);
% 判断是否改变格点
if rand<exp(-deE*beta_mc(jj))
X(tmp1,tmp4)=-X(tmp1,tmp4);
end
end
% 取样5万次,平衡时同样需要判断是否改变格点
for j=1:1:5*n
index=unidrnd(ns,1,2);
% 利用周期性边界条件,分别计算格点上下左右四个点行列坐标
tmp1=rem(index(1),ns)+1;tmp2=rem(index(1)+1,ns)+1;tmp3=rem(index(1)-1,ns)+1;
tmp4=rem(index(2),ns)+1;tmp5=rem(index(2)+1,ns)+1;tmp6=rem(index(2)-1,ns)+1;
% 计算改变格点方向后的能量变化
cen=X(tmp1,tmp4);right=X(tmp1,tmp5);left=X(tmp1,tmp6);
up= X(tmp2,tmp4);down= X(tmp3,tmp4);
deE=2*cen*(right+left+up+down);
% 判断是否改变格点
if rand<exp(-deE*beta_mc(jj))
X(tmp1,tmp4)=-X(tmp1,tmp4);
end
%计算一种特定分布时的平均磁化率
m(j)=abs(mean(mean(X)));
%计算一种特定分布时的平均能量
Xt1=X;Xt1(1,:)=[];Xt1=[Xt1; X(1,:)];
Xt2=X;Xt2(:,1)=[];Xt2=[Xt2, X(:,1)];
e(j)=-mean(mean(X.*Xt1+X.*Xt2));
end
% 特定温度下的统计量
m_bar(jj)=mean(m);
e_bar(jj)=mean(e);
cv_bar(jj)=beta_mc(jj)^2*ns^2*std(e)^2;
end
toc;
% 作图观察
figure(1);
plot(1./beta_mc,m_bar,'ko');
figure(2);
plot(1./beta_mc,e_bar,'ko');
figure(3);
plot(1./beta_mc,cv_bar,'ro');
20X20格点,从同一分布开始升温,分布进行了三批温度选择:
- β = ( 0.1 : 0.01 : 0.4 ) \beta=(0.1:0.01:0.4) β=(0.1:0.01:0.4) 31个温度点, T ∈ [ 2.5 , 10 ] T\in[2.5,10] T∈[2.5,10]马尔可夫链长度为1万次,采样5万次。
- T = ( 0.1 : 0.1 : 2.1 ) T=(0.1:0.1:2.1) T=(0.1:0.1:2.1) 21个温度点,$马尔可夫链长度为1万次,采样5万次。
- T = ( 2.1 : 0.01 : 2.4 ) T=(2.1:0.01:2.4) T=(2.1:0.01:2.4) 31个温度点,马尔可夫链长度为5万次,采样25万次。