BBH模型的Wilson loop及Nested Wilson loop计算
这里整理一下高阶拓扑的BBH模型的一些拓扑性质的计算,包括边界态,Wilson loop以及Nested Wilson loop计算.
{:.info}
Wilson loop
这个计算没有什么好说的,在其他模型里面也已经做过很多次了,比如BHZ模型Wilson loop计算这篇博客,废话不多说直接上代码1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104import numpy as np
import matplotlib.pyplot as plt
def Pauli():
s0 = np.array([[1,0],[0,1]])
sx = np.array([[0,1],[1,0]])
sy = np.array([[0,-1j],[1j,0]])
sz = np.array([[1,0],[0,-1]])
return s0,sx,sy,sz
#--------------------------------------------------------------
def hamset(kx,ky):
s0 = np.zeros([2,2],np.complex)
sx = np.zeros([2,2],np.complex)
sy = np.zeros([2,2],np.complex)
sz = np.zeros([2,2],np.complex)
ham = np.zeros([4,4],np.complex)
gamx = 0.5
gamy = 0.5
lamx = 1.0
lamy = 1.0
s0,sx,sy,sz = Pauli()
ham = (gamx + lamx*np.cos(kx))*np.kron(sx,s0) + lamx*np.sin(kx)*np.kron(-sy,sz) + (gamy + lamy*np.cos(ky))*np.kron(-sy,sy) + lamy*np.sin(ky)*np.kron(-sy,sx)
return ham
#--------------------------------------------------------------
def WilsonLoop():
nkx = 100
nky = 100
kxlist = np.linspace(-np.pi, np.pi, nkx)
kylist = np.linspace(-np.pi, np.pi, nky)
anglist = []
for ky in kylist:
wave1 = []
wave2 = []
for kx in kxlist:
eigval, eigvec = np.linalg.eigh(hamset(kx, ky))
if kx != np.pi:
wave1.append(eigvec[:, 0]) # 第一个占据态波函数
wave2.append(eigvec[:, 1]) # 第二个占据态波函数
else:
# 首位波函数相同,消除规范
wave1.append(wave1[0])
wave2.append(wave2[0])
Wan = np.eye(2,dtype=complex)
for i0 in range(nkx - 1):
F = np.zeros((2, 2), dtype = complex)
F[0, 0] = np.dot(wave1[i0 + 1].transpose().conj(), wave1[i0]) # 两个占据态波函数之间的交叠
F[1, 1] = np.dot(wave2[i0 + 1].transpose().conj(), wave2[i0])
F[0, 1] = np.dot(wave1[i0 + 1].transpose().conj(), wave2[i0])
F[1, 0] = np.dot(wave2[i0 + 1].transpose().conj(), wave1[i0])
Wan = np.dot(F, Wan)
eigval, eigvec = np.linalg.eig(Wan)
ang = np.log(eigval)/2/np.pi/1j
for i0 in range(2):
if np.real(ang[i0]) < 0:
ang[i0] += 1
ang = np.sort(ang)
anglist.append(ang.real)
return anglist
#------------------------------------------------
def WilsonLoop2():
nkx = 100
nky = 100
kxlist = np.linspace(-np.pi, np.pi, nkx)
kylist = np.linspace(-np.pi, np.pi, nky)
hn = hamset(0,0).shape[0]
Nocc = int(hn/2)
wave = np.zeros([hn,hn,nkx],np.complex)
anglist = []
for ky in kylist:
ix = 0
for kx in kxlist: # 计算沿着kx方向的Wilson loop
eigval, eigvec = np.linalg.eigh(hamset(kx, ky)) # 求解哈密顿量的本征矢量
if kx != np.pi:
wave[:,:,ix] = eigvec[:,:] # 存储所有的本征波函数,用来后面计算Wilson loop
ix += 1
else:
# 首尾波函数相同,消除规范
wave[:,:,nkx - 1] = wave[:,:,0]
Wan = np.eye(Nocc,dtype = complex)
F = np.zeros((Nocc, Nocc), dtype = complex) # Wannier Hamiltonian
for i0 in range(nkx - 1):
for i1 in range(Nocc): # 直接通过循环,只遍历占据态的波函数
for i2 in range(Nocc):
F[i1,i2] = np.dot(wave[:,i1,i0 + 1].transpose().conj(),wave[:,i2,i0])
Wan = np.dot(F, Wan)
eigval, eigvec = np.linalg.eig(Wan)
ang = np.log(eigval)/(2*np.pi*1j)
for i0 in range(2):
if np.real(ang[i0]) < 0:
ang[i0] += 1
ang = np.sort(ang)
anglist.append(ang.real)
return anglist
#------------------------------------------------------
def main():
a1 = WilsonLoop()
a2 = WilsonLoop2()
fig = plt.figure()
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
ax1.plot(a1)
ax2.plot(a2)
#--------------------------------------------------------
if __name__ == '__main__':
main()
这里采用了两种写代码的方式,第一种就是老老实实按照占据态的计算,而第二种则是将所有都整理到一起,需要用到占据态的时候在直接提取占据态对应的波函数来计算.
边界态
这部分同样没有太多需要讲解的,直接参考快速格林函数方法计算Chern绝缘体边界态这篇博客,直接上代码1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141import numpy as np
import matplotlib.pyplot as plt
import os
import time
import seaborn as sns
def Pauli():
s0 = np.array([[1,0],[0,1]])
sx = np.array([[0,1],[1,0]])
sy = np.array([[0,-1j],[1j,0]])
sz = np.array([[1,0],[0,-1]])
return s0,sx,sy,sz
#------------------------------------------------
def hamset(ki):
s0 = np.zeros([2,2],np.complex)
sx = np.zeros([2,2],np.complex)
sy = np.zeros([2,2],np.complex)
sz = np.zeros([2,2],np.complex)
H00 = np.zeros([4,4],np.complex)
H01 = np.zeros([4,4],np.complex)
gamx = 0.5
gamy = 0.5
lamx = 1.0
lamy = 1.0
s0,sx,sy,sz = Pauli()
H00 = (gamx + lamx*np.cos(ki))*np.kron(sx,s0) + lamx*np.sin(ki)*np.kron(-sy,sz)
H01 = (gamy + lamy/2.0)*np.kron(-sy,sy) + lamy/(2*1j)*np.kron(-sy,sx)
return H00,H01
# -----------------------------------------------------------------------
def Iteration(omega,ki):
err = 1e-16
eta = 0.01
iternum = 200
# H00 = np.zeros((2,2),np.complex128)
# H01 = np.zeros((2,2),np.complex128)
H00,H01 = hamset(ki)
epsiloni = H00
epsilons = H00
epsilons_t = H00
alphai = H01
betai = H01.T.conjugate()
omegac = omega + eta*1j
# s0,sx,sy,sz = Pauli()
s0 = np.eye(4)
for iter in range(iternum):
g0dem = omegac*s0 - epsiloni
g0 = np.linalg.inv(g0dem)
mat1 = np.dot(alphai,g0)
mat2 = np.dot(betai,g0)
g0 = np.dot(mat1,betai)
epsiloni = epsiloni + g0
epsilons = epsilons + g0
g0 = np.dot(mat2,alphai)
epsiloni = epsiloni + g0
epsilons_t = epsilons_t + g0
g0 = np.dot(mat1, alphai)
alphai = g0
g0 = np.dot(mat2,betai)
betai = g0
real_temp = np.sum(np.concatenate(np.abs(alphai)))
if (real_temp < err):
break
GLLdem = np.dot(omegac,s0) - epsilons
GLL = np.linalg.inv(GLLdem)
# GLL = epsilons
GLL = np.sum(np.concatenate(np.abs(GLL)))
GRRdem = np.dot(omegac,s0) - epsilons_t
GRR = np.linalg.inv(GRRdem)
GRR = np.sum(np.concatenate(np.abs(GRR)))
# GRR = epsilons_t
GBdem = np.dot(omegac,s0) - epsiloni
GB = np.linalg.inv(GBdem)
GB = np.sum(np.concatenate(np.abs(GB)))
# GB = epsiloni
return GLL,GRR,GB
#------------------------------------------------------------
def surface():
nx = 100
max_omg = 3
re = np.zeros((len(range(-nx,nx))**2,5))
con = 0
ix = -1
iy = -1
GLL = np.zeros((len(range(-nx,nx)),len(range(-nx,nx))))
GRR = np.zeros((len(range(-nx,nx)),len(range(-nx,nx))))
GB = np.zeros((len(range(-nx,nx)),len(range(-nx,nx))))
for i0 in range(-nx,nx):
kx = np.pi*i0/nx
for i1 in range(-nx,nx):
omg = max_omg*i1/nx
re1,re2,re3 = Iteration(omg,kx)
re[con,0] = kx
re[con,1] = omg
re[con,2] = re1
re[con,3] = re2
re[con,4] = re3
GLL[iy,ix] = np.log(re1)
GRR[iy,ix] = np.log(re2)
GB[iy,ix] = np.log(re3)
con += 1
iy += 1
ix += 1
iy = 0
# np.savetxt("GLL.dat", [kilist,re1list], fmt="%15.10f")
# np.savetxt("GRR.dat", [kilist,re1list], fmt="%15.10f")
# np.savetxt("GB.dat", [kilist,re1list], fmt="%15.10f")
np.savetxt("density.dat",re , fmt="%15.10f")
return GLL,GRR,GB
#-----------------------------------------------------------------------
def main():
os.chdir(os.getcwd())
tstart = time.time()
GLL,GRR,GB = surface()
tend = time.time()
#print(tend - tstart)
# 绘图
sns.set()
ax = sns.heatmap(GB)
plt.show()
#-----------------------------------------------------------------------
if __name__ == '__main__':
main()
Nested Wilson Loop
关于Nested Wilson loop的计算可以参考Electric multipole moments, topological multipole moment pumping, and chiral hinge states in crystalline insulators 这篇文章，这里主要说一下自己在学习计算的时候踩过的坑。
与计算Wilson loop相同，这里最主要的仍然是找到一个Wannier band basis，也就是文章的中的公式
其实在做计算的时候，最让人困扰的不过是公式中的一大堆符号对应的到底是什么，这里就来讲这个公式拆解开，一步一步的讲清楚里面的含义。这里假设你已经知道为什么要计算Nested Wilson loop，
我在这里就简单的阐述一下。首先要是体系的Wilson loop计算对应的Wannier哈密顿量的能带是有能隙的，也就是说你的体系是4带模型，那么当占据态是2条能带的时候，每个占据态能带会对应着一个Wannier center，
比如BHZ模型的两条占据态能带对应的Wannier band就是相互交叉的，而且因为Wilson loop与边界态之间的拓扑等价性，TI是有边界态的，所以其对应的Wilson loop在形状上就与边界态类似。而对于高阶拓扑相，
首先就是要使得边界态打开能隙，那么相对应的就是其Wilson loop计算得到的Wannier center随着某个动量参数的演化是不会相互交叉的，这一点在上面BBH模型中已经计算过了，所以此时就可以对某一个单独的Wannier band
计算它的Nested Wilson loop，所以首先第一步就是必须要明白什么样的情况下，是需要计算体系的Nested Wilson loop。
这里的$\sum_{n=1}^\text{Nocc}$不用讲太多，是需要对占据态进行求和，但是这个$n$其实表示的只是说哈密顿量的占据态，也就是说对于$\rvert u^n_\mathbf{k}\rangle$而言，这是哈密顿量的占据态波函数，$n$表示占据态其实是对$\rvert u^n_\mathbf{k}\rangle$
而言的，虽然$[v_{x,\mathbf{k}}^j]^n$中同样存在这个$n$，但是在这个地方$n$代表的不是占据态，在这里$j$代表的才是你选择的是哪一个Wannier band来计算其对应的Nested Wilson loop，也就是这里$j$代表的是你选择的是那个占据的Wannier band，而$n$在这里
表示的是一个Wannier哈密顿量本征矢量的第$n$个分量。假如$H_\text{Wann}$是Wannier哈密顿量，其满足
那么这里的$[v_{x,\mathbf{k}}^j]^n$表示的就是这个本征矢量的第$n$个分量，$j$则表示你选择是哪个本征值对应的本征矢量，也就是选择了哪一个Wannier band。这里的$x$则表示你在做Wilson loop的时候，是沿着哪个方向进行的，即就是讲上面公式中的$H_\text{Wann}$替换成你
构建的那个Wilson loop的哈密顿量就可以。
至于$\rvert u^n_\mathbf{k}\rangle$就很简单了，它表示的就是你的哈密顿量的本征态，当然了在计算的时候，还是要选择正确的占据态才可以。下面直接上代码，在其中同样做了注释
1 | import numpy as np |
完整的代码可以点击这里下载.
参考
公众号
相关内容均会在公众号进行同步，若对该Blog感兴趣，欢迎关注微信公众号。
{:.info}
|
yxliphy@gmail.com |
wechat
alipay
