BHZ模型哈密顿量与准一维体系的能带图（附Python代码）

发布时间：2019年11月27日2022年3月29日

最近修改时间：2022年3月29日

BHZ是一个量子自旋霍尔效应 (QSH) 的模型。BHZ模型是以文章“Quantum Spin Hall Effect and Topological Phase Transition in HgTe Quantum Wells”的作者姓名B. Andrei Bernevig, Taylor L. Hughes, Shou-Cheng Zhang的缩写命名的。

BHZ模型哈密顿量：

$H_{\mathrm{eff}}(k_x, k_y)= \begin{pmatrix} H(k) & 0 \\ 0 & H^{*}(-k) \\ \end{pmatrix}$

$H(k)= \varepsilon(k)+d_i(k) \sigma_i$

$d_1+id_2=A[\sin(k_x)+i\sin(k_ y)]$

$d_3=-2B[2-(M/2B)-\cos(k_x)-\cos(k_y)]$

$\varepsilon(k)=C-2D[2-\cos(k_x)-\cos(k_y)]$

通过反傅里叶变换，得到实空间的参数（参考：离散格子的傅里叶变换和反傅里叶变换）：

$H_0= \begin{pmatrix} E_s & 0 & 0 & 0 \\ 0 & E_p & 0 & 0 \\ 0 & 0 & E_s & 0 \\ 0 & 0 & 0 & E_p \\ \end{pmatrix}$

$H_1= \begin{pmatrix} V_{ss} & V_{sp} & 0 & 0 \\ -V_{sp}^* & V_{pp} & 0 & 0 \\ 0 & 0 & V_{ss} & V_{sp}^* \\ 0 & 0 & -V_{sp} & V_{pp} \\ \end{pmatrix}$

$H_2= \begin{pmatrix} V_{ss} & iV_{sp} & 0 & 0 \\ iV_{sp}^* & V_{pp} & 0 & 0 \\ 0 & 0 & V_{ss} & -iV_{sp}^* \\ 0 & 0 & -iV_{sp} & V_{pp} \\ \end{pmatrix}$

$E_s= C+M-4(D+B)/a^2$

$E_p= C-M-4(D-B)/a^2$

$V_{ss}=(D+B)/a^2$

$V_{pp}=(D-B)/a^2$

$V_{sp}=-iA/(2a)$

其中， $H_0$ 是在位能 (on-site energy)， $H_1$ 是 $x$ 方向的跃迁项 (hopping)， $H_2$ 是 $y$ 方向的跃迁项 (hopping)。

计算BHZ模型的准一维情况的能带，Python程序如下：

"""
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/2327
"""

import numpy as np
import matplotlib.pyplot as plt
from math import *   # 引入sqrt(), pi, exp等
import cmath  # 要处理复数情况，用到cmath.exp()
import functools  # 使用偏函数functools.partial()


def get_terms(A, B, C, D, M, a):
    E_s = C+M-4*(D+B)/(a**2)
    E_p = C-M-4*(D-B)/(a**2)
    V_ss = (D+B)/(a**2)
    V_pp = (D-B)/(a**2)
    V_sp = -1j*A/(2*a)
    H0 = np.zeros((4, 4))*(1+0j)  # 在位能 (on-site energy)
    H1 = np.zeros((4, 4))*(1+0j)  # x方向的跃迁 (hopping)
    H2 = np.zeros((4, 4))*(1+0j)  # y方向的跃迁 (hopping)
    H0[0, 0] = E_s
    H0[1, 1] = E_p
    H0[2, 2] = E_s
    H0[3, 3] = E_p

    H1[0, 0] = V_ss
    H1[1, 1] = V_pp
    H1[2, 2] = V_ss
    H1[3, 3] = V_pp
    H1[0, 1] = V_sp
    H1[1, 0] = -np.conj(V_sp)
    H1[2, 3] = np.conj(V_sp)
    H1[3, 2] = -V_sp

    H2[0, 0] = V_ss
    H2[1, 1] = V_pp
    H2[2, 2] = V_ss
    H2[3, 3] = V_pp
    H2[0, 1] = 1j*V_sp
    H2[1, 0] = 1j*np.conj(V_sp)
    H2[2, 3] = -1j*np.conj(V_sp)
    H2[3, 2] = -1j*V_sp
    return H0, H1, H2


def BHZ_model(k, A=0.3645/5, B=-0.686/25, C=0, D=-0.512/25, M=-0.01, a=1, N=100):  # 这边数值是不赋值时的默认参数
    H0, H1, H2 = get_terms(A, B, C, D, M, a)
    H00 = np.zeros((4*N, 4*N))*(1+0j)  # 元胞内，条带宽度为N
    H01 = np.zeros((4*N, 4*N))*(1+0j)  # 条带方向元胞间的跃迁
    for i in range(N):
        H00[i*4+0:i*4+4, i*4+0:i*4+4] = H0  # a:b代表 a <= x < b
        H01[i*4+0:i*4+4, i*4+0:i*4+4] = H1
    for i in range(N-1):
        H00[i*4+0:i*4+4, (i+1)*4+0:(i+1)*4+4] = H2
        H00[(i+1)*4+0:(i+1)*4+4, i*4+0:i*4+4] = np.conj(np.transpose(H2))
    H = H00 + H01 * cmath.exp(-1j * k) + H01.transpose().conj() * cmath.exp(1j * k)
    return H


def main():
    hamiltonian0 = functools.partial(BHZ_model, N=50)  # 使用偏函数，固定一些参数
    k = np.linspace(-pi, pi, 300)  # 300
    plot_bands_one_dimension(k, hamiltonian0)


def plot_bands_one_dimension(k, hamiltonian, filename='bands_1D'):
    dim = hamiltonian(0).shape[0]
    dim_k = k.shape[0]
    eigenvalue_k = np.zeros((dim_k, dim))  # np.zeros()里要用tuple
    i0 = 0
    for k0 in k:
        matrix0 = hamiltonian(k0)
        eigenvalue, eigenvector = np.linalg.eig(matrix0)
        eigenvalue_k[i0, :] = np.sort(np.real(eigenvalue[:]))
        i0 += 1
    for dim0 in range(dim):
        plt.plot(k, eigenvalue_k[:, dim0], '-k')  # -.
    # plt.savefig(filename + '.jpg')  # plt.savefig(filename+'.eps')
    plt.show()


if __name__ == '__main__':
    main()

计算结果为：

参数来源于文献[2]。可以看到，在能隙的地方出现了一个交叉的边缘态（交叉点略靠近导带），这个边缘态是受到拓扑保护的。

这里选的M参数等于-0.01，体系是拓扑的。如果M改为0.002，拓扑边缘态消失，体系是平庸非拓扑的。

参考资料：

[1] BHZ模型原始文献“Quantum Spin Hall Effect and Topological Phase Transition in HgTe Quantum Wells”

[2] 文献“Numerical study of the topological Anderson insulator in HgTe/CdTe quantum wells”

【未经同意请勿转载。引用请注明出处：https://www.guanjihuan.com】

5,522 次浏览

Published by guanjihuan

View all posts by guanjihuan

19 thoughts on “BHZ模型哈密顿量与准一维体系的能带图（附Python代码）”

匿名说道：

2022年6月19日 11:07

关博士您好，请问做反傅里叶变换时是不是要先知道实空间的晶格排列呢？

回复
1. guanjihuan说道：
  
  2022年6月19日 14:03
  
  哈密顿量在倒空间的形式已经蕴含着晶格的信息，比较常见的是方格子。需要注意的是：这个晶格排列只是一个有效的模型，不是实际材料的晶格排列，这还是有区别的。实际晶格排列的材料的能带和波函数是通过第一性原理计算给出。傅里叶变换和反傅里叶变换可参考：离散格子的傅里叶变换和反傅里叶变换。
  
  回复

老师您好，我还有个疑问，是后面在重写程序时发现的。
1，重写时我发现在def BHZ_model函数中，交换哈密顿量量中的子矩阵H2似乎不影响计算，比如改成以下形式
for i in range(N-1):
H00[i*4+0:i*4+4, (i+1)*4+0:(i+1)*4+4] = np.conj(np.transpose(H2))
H00[(i+1)*4+0:(i+1)*4+4, i*4+0:i*4+4] = H2
所以这个H2的含义是指向y方向的下一个原子跃迁还是“被”y的反方向下一个原子跃迁呢？
2，为了验证我的猜测，我对程序稍作了修改，将条带的方向旋转90度，计算了大小为N*2的带子，然后发现需要把np.conj(np.transpose(H2))与哈密顿矩阵的H2交换位置才能出现边界态，否则图像看起来是个平庸绝缘体。请问在将这个H2拼接进H矩阵时有什么注意事项吗？

def BHZ_model(k, A=0.3645/5, B=-0.686/25, C=0, D=-0.512/25, M=-0.01, a=1, N=50):  # 这边数值是不赋值时的默认参数
    H0, H1, H2 = get_terms(A, B, C, D, M, a) 
    # 修改成沿y方向无限长，x方向宽度为N的带子，仿真其中的一段(包含N*2个元素)
    # 实空间标号顺序是从左到右，从下到上，共N*2
    H00 = np.zeros((4*N*2, 4*N*2))*(1+0j)  # 元胞内，条带宽度为N
    H01 = np.zeros((4*N*2, 4*N*2))*(1+0j)  # 条带方向元胞间的跃迁
    for i in range(2*N): # 元胞内部的哈密顿量，包含占位能、上下两段各自的x方向，和之间的y方向
        H00[i*4+0:i*4+4, i*4+0:i*4+4] = H0  # a:b代表 a <= x < b
    for i in range(N-1):
        H00[i*4+0:i*4+4, (i+1)*4+0:(i+1)*4+4] = np.conj(np.transpose(H1))
        H00[(i+1)*4+0:(i+1)*4+4, i*4+0:i*4+4] = H1
        H00[(i+N)*4+0:(i+N)*4+4, (i+N+1)*4+0:(i+N+1)*4+4] = np.conj(np.transpose(H1))
        H00[(i+N+1)*4+0:(i+N+1)*4+4, (i+N)*4+0:(i+N)*4+4] = H1
        H00[i*4+0:i*4+4, (i+N)*4+0:(i+N)*4+4] = np.conj(np.transpose(H2))
        H00[(i+N)*4+0:(i+N)*4+4, i*4+0:i*4+4] = H2
    
    # 元胞与元胞之间的哈密顿量，可能是周期性原因要额外补上Bloch的相位exp(ik*Rn)? Rn在这里是晶格矢量定值，已翻倍为a*2
    for i in range(N-1):
        H01[i*4+0:i*4+4, (i+N)*4+0:(i+N)*4+4] = H2 * cmath.exp(-1j * k * 2*a)
        H01[(i+N)*4+0:(i+N)*4+4, i*4+0:i*4+4] = np.conj(np.transpose(H2)) * cmath.exp(1j * k * 2*a)
    
    H = H00 + H01
    return H
def main():
    hamiltonian0 = functools.partial(BHZ_model, N=25)  # 使用偏函数，固定一些参数
    k = np.linspace(-pi/2, pi/2, 300)  # 300  # 沿着y方向的原胞宽度翻倍了，k似乎应该减半
    plot_bands_one_dimension(k, hamiltonian0)

guanjihuan说道：

2021年7月27日 14:50

只要写的时候全程固定一个方向，应该是对结果没有影响的。如果是关于x方向的跃迁相反，那么kx方向的能带也要相反。

回复
1. 匿名说道：
  
  2021年7月27日 16:33
  
  谢谢老师，我按照您说的，只交换元胞间的哈密顿子矩阵方向，结果也正确了o(￣▽￣)o
  
  回复

匿名说道：

2021年5月7日 11:21

请问老师。如果哈密顿量里出现sinkx乘sinky咋办，没办法把他俩先分开写到不同矩阵里，再做逆变换，，

回复
1. guanjihuan说道：
  
  2021年5月7日 21:08
  
  参考这篇：离散格子的傅里叶变换和反傅里叶变换
  
  回复
  1. 匿名说道：
    
    2021年5月7日 22:22
    
    谢谢老师！
    
    回复
匿名说道：

2020年12月28日 19:54

请问大神这里的逆傅里叶变换是不是很复杂，一个矩阵元都8项了，解析好难推

回复
1. guanjihuan说道：
  
  2020年12月29日 09:17
  
  找到规律就可以直接写了，不用重新推导。可以参考离散格子的傅里叶变换和反傅里叶变换中“傅里叶变换常见的形式”部分。
  
  回复
  1. 匿名说道：
    
    2020年12月29日 09:32
    
    多谢多谢
    
    回复
匿名说道：

2020年12月27日 16:25

非常感谢

回复
匿名说道：

2020年12月22日 18:44

请问博主画能带为什么要把实空间的形式写出来，不可以直接把K空间的哈密顿量对角化得到能带吗

回复
1. 匿名说道：
  
  2020年12月22日 21:55
  
  请问是需要把实空间的哈密顿量再做一次傅里叶变换，然后再取一个方向（周期性）的哈密顿量对角化画出能带吗
  
  回复
2. guanjihuan说道：
  
  2020年12月22日 23:42
  
  已知的哈密顿量是H(kx,ky)，是一个二维体系的情况，直接对角化后是二维体系的能带E(kx,ky)。如果需要算有限宽度的条带，即准一维的情况，那么需要先知道实空间的跃迁项，把条带的最小元胞找出来，在一个方向傅里叶变换，例如x方向，对角化后得到E(kx)。
  
  回复
  1. 匿名说道：
    
    2020年12月25日 22:07
    
    多谢多谢
    
    回复
  2. 匿名说道：
    
    2020年12月26日 21:20
    
    请问这个二维的哈密顿量直接画出来能带是不是应该有一个锥啊，我对角化画出来不管取d>64埃还是d<64埃的参数都是一个锥。无法区分拓扑平庸和不平庸
    
    回复
    1. 匿名说道：
      
      2020年12月26日 21:24
      
      说错了，不是一个锥，是两个山谷（valley），或者是两把对顶的伞面那样
      
      回复
      1. guanjihuan说道：
        
        2020年12月26日 21:51
        
        从二维的能带图中是无法看出是否是拓扑的，因为拓扑绝缘体也是绝缘体，能带是有带隙的。拓扑的信息包含在波函数中，而不在本征值中。只有当存在边缘时，由于体边对应关系（bulk-boundary corresponding），才会有交叉状的边缘态。
        
        回复

Published by guanjihuan

19 thoughts on “BHZ模型哈密顿量与准一维体系的能带图（附Python代码）”

发表评论 取消回复

发表评论取消回复