BHZ模型的能带图参考:BHZ模型哈密顿量与准一维体系的能带图(附Python代码)。
BHZ模型哈密顿量:
一、 自旋陈数
因为自旋向上和自旋向下是两个相互独立的子体系(对应的陈数刚好相反),自旋是守恒的,因此可以定义自旋陈数。
陈数的计算方法为:陈数Chern number的计算(高效法,附Python/Matlab代码)。
代码如下:
"""
This code is supported by the website: https://www.guanjihuan.com
The newest version of this code is on the web page: https://www.guanjihuan.com/archives/5778
"""
import numpy as np
import matplotlib.pyplot as plt
from math import * # 引入pi, cos等
import cmath
import time
def hamiltonian1(kx, ky): # Half BHZ for spin up
A=0.3645/5
B=-0.686/25
C=0
D=-0.512/25
M=-0.01
matrix = np.zeros((2, 2))*(1+0j)
varepsilon = C-2*D*(2-cos(kx)-cos(ky))
d3 = -2*B*(2-(M/2/B)-cos(kx)-cos(ky))
d1_d2 = A*(sin(kx)+1j*sin(ky))
matrix[0, 0] = varepsilon+d3
matrix[1, 1] = varepsilon-d3
matrix[0, 1] = np.conj(d1_d2)
matrix[1, 0] = d1_d2
return matrix
def hamiltonian2(kx, ky): # Half BHZ for spin down
A=0.3645/5
B=-0.686/25
C=0
D=-0.512/25
M=-0.01
matrix = np.zeros((2, 2))*(1+0j)
varepsilon = C-2*D*(2-cos(-kx)-cos(-ky))
d3 = -2*B*(2-(M/2/B)-cos(-kx)-cos(-ky))
d1_d2 = A*(sin(-kx)+1j*sin(-ky))
matrix[0, 0] = varepsilon+d3
matrix[1, 1] = varepsilon-d3
matrix[0, 1] = d1_d2
matrix[1, 0] = np.conj(d1_d2)
return matrix
def main():
start_clock = time.perf_counter()
delta = 0.1
for i0 in range(2):
if i0 == 0:
hamiltonian = hamiltonian1
else:
hamiltonian = hamiltonian2
chern_number = 0 # 陈数初始化
for kx in np.arange(-pi, pi, delta):
print(kx)
for ky in np.arange(-pi, pi, delta):
H = hamiltonian(kx, ky)
eigenvalue, eigenvector = np.linalg.eig(H)
vector = eigenvector[:, np.argsort(np.real(eigenvalue))[0]] # 价带波函数
H_delta_kx = hamiltonian(kx+delta, ky)
eigenvalue, eigenvector = np.linalg.eig(H_delta_kx)
vector_delta_kx = eigenvector[:, np.argsort(np.real(eigenvalue))[0]] # 略偏离kx的波函数
H_delta_ky = hamiltonian(kx, ky+delta)
eigenvalue, eigenvector = np.linalg.eig(H_delta_ky)
vector_delta_ky = eigenvector[:, np.argsort(np.real(eigenvalue))[0]] # 略偏离ky的波函数
H_delta_kx_ky = hamiltonian(kx+delta, ky+delta)
eigenvalue, eigenvector = np.linalg.eig(H_delta_kx_ky)
vector_delta_kx_ky = eigenvector[:, np.argsort(np.real(eigenvalue))[0]] # 略偏离kx和ky的波函数
Ux = np.dot(np.conj(vector), vector_delta_kx)/abs(np.dot(np.conj(vector), vector_delta_kx))
Uy = np.dot(np.conj(vector), vector_delta_ky)/abs(np.dot(np.conj(vector), vector_delta_ky))
Ux_y = np.dot(np.conj(vector_delta_ky), vector_delta_kx_ky)/abs(np.dot(np.conj(vector_delta_ky), vector_delta_kx_ky))
Uy_x = np.dot(np.conj(vector_delta_kx), vector_delta_kx_ky)/abs(np.dot(np.conj(vector_delta_kx), vector_delta_kx_ky))
F = cmath.log(Ux*Uy_x*(1/Ux_y)*(1/Uy))
# 陈数(chern number)
chern_number = chern_number + F
chern_number = chern_number/(2*pi*1j)
if i0 == 0:
chern_number_up = chern_number
else:
chern_number_down = chern_number
spin_chern_number = (chern_number_up-chern_number_down)/2
print('Chern number for spin up = ', chern_number_up)
print('Chern number for spin down = ', chern_number_down)
print('Spin chern number = ', spin_chern_number)
end_clock = time.perf_counter()
print('CPU执行时间(min)=', (end_clock-start_clock)/60)
if __name__ == '__main__':
main()
计算结果为:

二、Z2不变量
如果自旋不守恒,一般情况下自旋陈数失去了物理意义(但可以通过特殊的处理方式来定义自旋陈数,这里不进行讨论)。在时间反演对称性下可以定义Z2不变量,计算方法见参考资料[1-5]。
计算BHZ模型Z2不变量的代码(仍以自旋守恒的体系为例):
"""
This code is supported by the website: https://www.guanjihuan.com
The newest version of this code is on the web page: https://www.guanjihuan.com/archives/5778
"""
import numpy as np
import matplotlib.pyplot as plt
from math import * # 引入pi, cos等
import cmath
import time
def hamiltonian(kx, ky): # BHZ模型
A=0.3645/5
B=-0.686/25
C=0
D=-0.512/25
M=-0.01
matrix = np.zeros((4, 4))*(1+0j)
varepsilon = C-2*D*(2-cos(kx)-cos(ky))
d3 = -2*B*(2-(M/2/B)-cos(kx)-cos(ky))
d1_d2 = A*(sin(kx)+1j*sin(ky))
matrix[0, 0] = varepsilon+d3
matrix[1, 1] = varepsilon-d3
matrix[0, 1] = np.conj(d1_d2)
matrix[1, 0] = d1_d2
varepsilon = C-2*D*(2-cos(-kx)-cos(-ky))
d3 = -2*B*(2-(M/2/B)-cos(-kx)-cos(-ky))
d1_d2 = A*(sin(-kx)+1j*sin(-ky))
matrix[2, 2] = varepsilon+d3
matrix[3, 3] = varepsilon-d3
matrix[2, 3] = d1_d2
matrix[3, 2] = np.conj(d1_d2)
return matrix
def main():
start_clock = time.perf_counter()
delta = 0.1
Z2 = 0 # Z2数
for kx in np.arange(-pi, 0, delta):
print(kx)
for ky in np.arange(-pi, pi, delta):
H = hamiltonian(kx, ky)
eigenvalue, eigenvector = np.linalg.eig(H)
vector = eigenvector[:, np.argsort(np.real(eigenvalue))[0]] # 价带波函数1
vector2 = eigenvector[:, np.argsort(np.real(eigenvalue))[1]] # 价带波函数2
H_delta_kx = hamiltonian(kx+delta, ky)
eigenvalue, eigenvector = np.linalg.eig(H_delta_kx)
vector_delta_kx = eigenvector[:, np.argsort(np.real(eigenvalue))[0]] # 略偏离kx的波函数1
vector_delta_kx2 = eigenvector[:, np.argsort(np.real(eigenvalue))[1]] # 略偏离kx的波函数2
H_delta_ky = hamiltonian(kx, ky+delta)
eigenvalue, eigenvector = np.linalg.eig(H_delta_ky)
vector_delta_ky = eigenvector[:, np.argsort(np.real(eigenvalue))[0]] # 略偏离ky的波函数1
vector_delta_ky2 = eigenvector[:, np.argsort(np.real(eigenvalue))[1]] # 略偏离ky的波函数2
H_delta_kx_ky = hamiltonian(kx+delta, ky+delta)
eigenvalue, eigenvector = np.linalg.eig(H_delta_kx_ky)
vector_delta_kx_ky = eigenvector[:, np.argsort(np.real(eigenvalue))[0]] # 略偏离kx和ky的波函数1
vector_delta_kx_ky2 = eigenvector[:, np.argsort(np.real(eigenvalue))[1]] # 略偏离kx和ky的波函数2
Ux = dot_and_det(vector, vector_delta_kx, vector2, vector_delta_kx2)
Uy = dot_and_det(vector, vector_delta_ky, vector2, vector_delta_ky2)
Ux_y = dot_and_det(vector_delta_ky, vector_delta_kx_ky, vector_delta_ky2, vector_delta_kx_ky2)
Uy_x = dot_and_det(vector_delta_kx, vector_delta_kx_ky, vector_delta_kx2, vector_delta_kx_ky2)
F = np.imag(cmath.log(Ux*Uy_x*np.conj(Ux_y)*np.conj(Uy)))
A = np.imag(cmath.log(Ux))+np.imag(cmath.log(Uy_x))+np.imag(cmath.log(np.conj(Ux_y)))+np.imag(cmath.log(np.conj(Uy)))
Z2 = Z2 + (A-F)/(2*pi)
print('Z2 = ', Z2) # Z2数
end_clock = time.perf_counter()
print('CPU执行时间(min)=', (end_clock-start_clock)/60)
def dot_and_det(a1, b1, a2, b2): # 内积组成的矩阵对应的行列式
x1 = np.dot(np.conj(a1), b1)
x2 = np.dot(np.conj(a2), b2)
x3 = np.dot(np.conj(a1), b2)
x4 = np.dot(np.conj(a2), b1)
return x1*x2-x3*x4
if __name__ == '__main__':
main()
计算结果为:

在自旋守恒的情况下,自旋陈数和Z2不变量是一致的。
Matlab代码为(内容同上):
% This code is supported by the website: https://www.guanjihuan.com
% The newest version of this code is on the web page: https://www.guanjihuan.com/archives/5778
clear;clc;
delta=0.1;
Z2=0;
for kx=-pi:0.1:0
for ky=-pi:0.1:pi
[V1,V2]=get_vector(Hamiltonian(kx,ky));
[Vkx1,Vkx2]=get_vector(Hamiltonian(kx+delta,ky)); % 略偏离kx的波函数
[Vky1,Vky2]=get_vector(Hamiltonian(kx,ky+delta)); % 略偏离ky的波函数
[Vkxky1,Vkxky2]=get_vector(Hamiltonian(kx+delta,ky+delta)); % 略偏离kx,ky的波函数
Ux = dot_and_det(V1, Vkx1, V2, Vkx2);
Uy = dot_and_det(V1, Vky1, V2, Vky2);
Ux_y = dot_and_det(Vky1, Vkxky1, Vky2, Vkxky2);
Uy_x = dot_and_det(Vkx1, Vkxky1, Vkx2, Vkxky2);
F=imag(log(Ux*Uy_x*(conj(Ux_y))*(conj(Uy))));
A=imag(log(Ux))+imag(log(Uy_x))+imag(log(conj(Ux_y)))+imag(log(conj(Uy)));
Z2 = Z2+(A-F)/(2*pi);
end
end
Z2= mod(Z2, 2)
function dd = dot_and_det(a1,b1,a2,b2) % 内积组成的矩阵对应的行列式
x1=a1'*b1;
x2=a2'*b2;
x3=a1'*b2;
x4=a2'*b1;
dd =x1*x2-x3*x4;
end
function [vector_new_1, vector_new_2] = get_vector(H)
[vector,eigenvalue] = eig(H);
[eigenvalue, index]=sort(diag(eigenvalue), 'ascend');
vector_new_2 = vector(:, index(2));
vector_new_1 = vector(:, index(1));
end
function H=Hamiltonian(kx,ky) % BHZ模型
A=0.3645/5;
B=-0.686/25;
C=0;
D=-0.512/25;
M=-0.01;
H=zeros(4,4);
varepsilon = C-2*D*(2-cos(kx)-cos(ky));
d3 = -2*B*(2-(M/2/B)-cos(kx)-cos(ky));
d1_d2 = A*(sin(kx)+1j*sin(ky));
H(1, 1) = varepsilon+d3;
H(2, 2) = varepsilon-d3;
H(1, 2) = conj(d1_d2);
H(2, 1) = d1_d2 ;
varepsilon = C-2*D*(2-cos(-kx)-cos(-ky));
d3 = -2*B*(2-(M/2/B)-cos(-kx)-cos(-ky));
d1_d2 = A*(sin(-kx)+1j*sin(-ky));
H(3, 3) = varepsilon+d3;
H(4, 4) = varepsilon-d3;
H(3, 4) = d1_d2 ;
H(4, 3) = conj(d1_d2);
end
额外说明:代码严格来说可能还需要处理下高对称点和高对称线的情况,具体可以看参考资料中的说明。目前这个例子计算结果没有出现问题,估计是这些位置的数值刚好为零。
参考资料:
[1] Z2拓扑不变量与拓扑绝缘体
[2] Time reversal polarization and a Z2 adiabatic spin pump
[3] Chern number and Z2 invariant
[5] First-principles calculation of Z2topological invariants within the FP-LAPW formalism
【说明:本站主要是个人的一些笔记和代码分享,内容可能会不定期修改。为了使全网显示的始终是最新版本,这里的文章未经同意请勿转载。引用请注明出处:https://www.guanjihuan.com】
你好,请问一下对于BHZ模型,如何构造时间反演算符T,使得满足时间反演对称性TH_kT^{-1}=H_{-k}^*且T^2=-1
大概跟这个差不多,这个算符矩阵是二维的:时间反演算符 Time Reversal Operator。
如果这个算符矩阵和二阶的单位矩阵进行张量积就可以变成四维的,但要注意张量积的前后顺序,有一个顺序是对的,刚好和BHZ模型的分块形式对应上。张量积参考这篇:泡利矩阵以及泡利矩阵的张量积。
请问如果我要算的哈密顿有四条价带,function dd = dot_and_det(a1,b1,a2,b2) % 内积组成的矩阵对应的行列式
x1=a1'*b1;
x2=a2'*b2;
x3=a1'*b2;
x4=a2'*b1;
dd =x1*x2-x3*x4;
end
那么这个偏函数里求得dd,应该是8X8的矩阵的行列式吗
是2x2矩阵的行列式。内积后是一个数,不是矩阵。
[eigenvalue, index]=sort(diag(eigenvalue), 'ascend');
vector_new_2 = vector(:, index(2));
vector_new_1 = vector(:, index(1));
请问楼主上面的代码是取了导带来计算吗,取价带3,4来算的话会不会Z2是相反数啊
取的是价带。ascend是上升的意思,最前面两个带是价带。另外两条带我没考虑过,你可以算算。理论上如果考虑所有的带,Z2应该等于0。
请问这里取的两条带必须要是不同自旋贡献的兼并的两条带吗,要算Z2是不是模型必须要有这个自由度啊
如果算自旋陈数,那么需要两个自旋是独立的。算Z2没有这个要求。
至于Z2是否必须考虑自旋,目前我看到的都是有考虑,且需要满足时间反演对称性。其他情况我不是很了解,没深入研究过,你可以查下文献看看。
多谢指导
请问bhz模型的陈数是两个自旋陈数之和(-1)+(1)=0吗
不是之和,是相减后除以2。
作者可以出一篇运用wannier function center的计算Z2不变量的文章吗?然后比较说明一下和用贝里曲率算法的比较说明,这两个方法的对比,看了文献,并不是非常的理解。谢谢
这方面的文章我看的也不多,目前不是很了解。之后有空的时候可能会考虑。你可以多找些文献看看,看不同文献中的描述。
请问这里面的积分空间是根据什么来选取的呢?
如果是计算Z2不变量,在这个方法中是选取一半的布里渊区。
嗯嗯,是的。对于简单的BHZ正方晶格的积分选取是比较显而易见的。那如果对于非正方晶格比如三角晶格或六角晶格,就会出现缺角啊,还是说根据话高对称点之间的能带图的取法来选取呢?
一样的,可以通过程序把积分限制半个布里渊内,不管是什么图形。
如何还是按照离散矩阵的形式来积分吗?
可以参考这一篇中的积分方法:Haldane模型中陈数的计算(附Python代码)
作者可以出一个matlab的求Z2不变量的参考程序吗?谢谢
在博文的最后更新了matlab代码,供参考。内容是一样的。
十分感谢博主,对于一个初次接触拓扑的研究生帮助甚大,希望可以多交流。
博主您好,我想请问一下,对于这样的二能带系统,里面的A具体是怎么定义的,能展开说说吗?因为对于二能带,A(联络)应该是矩阵形式,那么对应于U,U的具体形式是怎么样的呢?看程序看的不太明显…麻烦了
公式参考这篇中给的文献:陈数Chern number的计算(高效法,附Python/Matlab代码)。A或U都是标量吧,不是矩阵,内积的结果是一个值。
能带有两个指标很明显A的分量为A11,A12,A21,A22,数字为能带指标啊…为啥不是矩阵呢
是有形成一个矩阵,不同能带上的波函数是正交的,非对角线上的矩阵元(内积值)应该为零。U最后是取了一个行列式,具体公式你可以参照参考资料[2]的PPT,或者其他参考资料。我对这方面也没有仔细研究过,只是简单算了下。
参考文献3中的6式,内积的两个波函数分别是两个相似变换矩阵(2n列向量排列形成),两个方阵相乘再求行列式,
应该是这样不,感觉和代码里不一样