模型和能带, 学术

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

发布时间：

最近修改时间：2022年10月21日

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), dtype=complex)  # 在位能 (on-site energy)
    H1 = np.zeros((4, 4), dtype=complex)  # x方向的跃迁 (hopping)
    H2 = np.zeros((4, 4), dtype=complex)  # 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), dtype=complex)  # 元胞内，条带宽度为N
    H01 = np.zeros((4*N, 4*N), dtype=complex)  # 条带方向元胞间的跃迁
    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”

7,777 次浏览

【说明：本站主要是个人的一些笔记和代码分享，内容可能会不定期修改。为了使全网显示的始终是最新版本，这里的文章未经同意请勿转载。引用请注明出处：https://www.guanjihuan.com】

头像

Published by guanjihuan

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

  1. 关老师，您好，我想请教一个三维二阶模型（宋志达教授的文章（(d − 2)-Dimensional Edge States of Rotation Symmetry Protected Topological States）的能带图。模型的哈密顿量为H=(M-t(cos(kx)+cos(ky)+cos(kz)))τ0σzs0+Asin(kx)τ0σxsx+Asin(ky)τ0σxsy+Asin(kz)τ0σxsz+Δ(cos(kx)-cos(ky))τyσys0 #（其中τ，σ，s都是泡利矩阵）

    def SONG_model(k, t=1, A=1, M=2, Δ=0.2，a=1, N=100): # 这边数值是不赋值时的默认参数
    H0, H1, H2,H3 = get_terms(t, A, Δ, M, a) # H0表示在位能，H1表示x方向跃迁，H2表示y方向跃迁，H3表示z方向跃迁。只有z方向具有周期性，在x和y方向限定宽度

    H00 = np.zeros((8*N, 8*N), dtype=complex) # 元胞内，宽度为N
    H01 = np.zeros((8*N, 8*N), dtype=complex) # z方向元胞间的跃迁
    for i in range(N):
    H00[i*8+0:i*8+8, i*8+0:i*8+8] = H0 # a:b代表 a <= x < b
    H01[i*8+0:i*8+8, i*8+0:i*8+8] = H3
    for i in range(N-1):
    H00[i*8+0:i*8+4, (i+1)*8+0:(i+1)*8+8] = H1 #x方向
    H00[(i+1)*8+0:(i+1)*8+8, i*8+0:i*8+8] = np.conj(np.transpose(H1))

    H00[i*8+0:i*8+4, (i+1)*8+0:(i+1)*8+8] = H2 #y方向
    H00[(i+1)*8+0:(i+1)*8+8, i*8+0:i*8+8] = np.conj(np.transpose(H2))
    H = H00 + H01 * cmath.exp(-1j * k) + H01.transpose().conj() * cmath.exp(1j * k)
    return H
    我的问题是：同时限定x和y的宽度时，这个哈密顿量可以这样写吗？我感觉有些错误，请关老师指点一下，期待您的回复

    回复

      1. 关老师，您讲的对，这个宽度是限定x和y的宽度。关于自旋我是做直积使哈密顿量变为8*8的矩阵。我仿照您发的BBH模型写出了H00，但是z方向的跃迁H01感觉有错误，请您指点
        clear
        close all
        clc
        Nx=5;
        Ny=5;
        a=1;
        d=0.2;
        t=1;
        M=2;
        H0=zeros(8,8);%在位能
        H1=zeros(8,8);%x方向的跃迁
        H2=zeros(8,8);%y方向的跃迁
        H3=zeros(8,8);%z方向的跃迁

        H0(1,1)=M;
        H0(2,2)=M;
        H0(3,3)=-M;
        H0(4,4)=-M;
        H0(5,5)=M;
        H0(6,6)=M;
        H0(7,7)=-M;
        H0(8,8)=-M;

        H1(1,1)=-t/2;
        H1(2,2)=-t/2;
        H1(3,3)=t/2;
        H1(4,4)=t/2;
        H1(5,5)=-t/2;
        H1(6,6)=-t/2;
        H1(7,7)=t/2;
        H1(8,8)=t/2;
        H1(1,4)=-1i*a/2;
        H1(2,3)=-1i*a/2;
        H1(3,2)=-1i*a/2;
        H1(4,1)=-1i*a/2;
        H1(5,8)=-1i*a/2;
        H1(6,7)=-1i*a/2;
        H1(7,6)=-1i*a/2;
        H1(8,5)=-1i*a/2;
        H1(1,7)=-d/2;
        H1(2,8)=-d/2;
        H1(3,5)=d/2;
        H1(4,6)=d/2;
        H1(5,3)=d/2;
        H1(6,4)=d/2;
        H1(7,1)=-d/2;
        H1(8,2)=-d/2;

        H2(1,1)=-t/2;
        H2(2,2)=-t/2;
        H2(3,3)=t/2;
        H2(4,4)=t/2;
        H2(5,5)=-t/2;
        H2(6,6)=-t/2;
        H2(7,7)=t/2;
        H2(8,8)=t/2;
        H2(1,4)=-a/2;
        H2(2,3)=a/2;
        H2(3,2)=-a/2;
        H2(4,1)=a/2;
        H2(5,8)=-a/2;
        H2(6,7)=a/2;
        H2(7,6)=-a/2;
        H2(8,5)=a/2;
        H2(1,7)=d/2;
        H2(2,8)=d/2;
        H2(3,5)=-d/2;
        H2(4,6)=-d/2;
        H2(5,3)=-d/2;
        H2(6,4)=-d/2;
        H2(7,1)=d/2;
        H2(8,2)=d/2;

        H3(1,1)=-t/2;
        H3(2,2)=-t/2;
        H3(3,3)=t/2;
        H3(4,4)=t/2;
        H3(5,5)=-t/2;
        H3(6,6)=-t/2;
        H3(7,7)=t/2;
        H3(8,8)=t/2;
        H3(1,3)=-1i*a/2;
        H3(2,4)=1i*a/2;
        H3(3,1)=-1i*a/2;
        H3(4,2)=1i*a/2;
        H3(5,7)=-1i*a/2;
        H3(6,8)=1i*a/2;
        H3(7,5)=-1i*a/2;
        H3(8,6)=1i*a/2;

        h00=zeros(8*Nx*Ny,8*Nx*Ny);%元胞内，宽度为Nx*Ny
        h01=zeros(8*Nx*Ny,8*Nx*Ny);%z方向，元胞间的跃迁
        for x=0:Nx-1
        for y=0:Ny-1
        h00(x*Ny*8+y*8+1:x*Ny*8+y*8+8,x*Ny*8+y*8+1:x*Ny*8+y*8+8)=H0;%在位能
        h01(x*Ny*8+y*8+1:x*Ny*8+y*8+8,x*Ny*8+y*8+1:x*Ny*8+y*8+8)=H3;%元胞间z方向跃迁

        for x=0:Nx-2
        for y=0:Ny-2
        h00(x*Ny*8+y*8+1:x*Ny*8+y*8+8,(x+1)*Ny*8+y*8+1:(x+1)*Ny*8+y*8+8)=H1;%元胞内x跃迁
        h00((x+1)*Ny*8+y*8+1:(x+1)*Ny*8+y*8+8,x*Ny*8+y*8+1:x*Ny*8+y*8+8)=H1';%元胞内x跃迁
        h00(x*Ny*8+y*8+1:x*Ny*8+y*8+8,x*Ny*8+(y+1)*8+1:x*Ny*8+(y+1)*8+8)=H2;%元胞内y跃迁
        h00(x*Ny*8+(y+1)*8+1:x*Ny*8+(y+1)*8+8,x*Ny*8+y*8+1:x*Ny*8+y*8+8)=H2';%元胞内y跃迁
        end
        end
        end
        end
        kz=linspace(-pi,pi,100);
        for i=1:100
        hh=h00+h01*exp(1i*kz(i))+h01'*exp(-1i*kz(i));
        e(:,i)=sort(eig(hh));
        end
        plot(kz,e,'b')

        回复

        2. 可能没太大问题吧。具体的你可以自己计算验证下，看结果是否合理以及是否符合物理或预期

          回复

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

    回复

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

      回复

  3. 老师您好，我还有个疑问，是后面在重写程序时发现的。
    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)

    回复

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

      回复

      1. 谢谢老师，我按照您说的，只交换元胞间的哈密顿子矩阵方向，结果也正确了o(￣▽￣)o

        回复

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

    回复

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

    回复

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

    回复

    1. 请问是需要把实空间的哈密顿量再做一次傅里叶变换，然后再取一个方向（周期性）的哈密顿量对角化画出能带吗

      回复

    2. 已知的哈密顿量是H(kx,ky)，是一个二维体系的情况，直接对角化后是二维体系的能带E(kx,ky)。如果需要算有限宽度的条带，即准一维的情况，那么需要先知道实空间的跃迁项，把条带的最小元胞找出来，在一个方向傅里叶变换，例如x方向，对角化后得到E(kx)。

      回复

      2. 请问这个二维的哈密顿量直接画出来能带是不是应该有一个锥啊，我对角化画出来不管取d>64埃还是d<64埃的参数都是一个锥。无法区分拓扑平庸和不平庸

        回复

        1. 说错了，不是一个锥，是两个山谷（valley），或者是两把对顶的伞面那样

          回复

          1. 从二维的能带图中是无法看出是否是拓扑的，因为拓扑绝缘体也是绝缘体，能带是有带隙的。拓扑的信息包含在波函数中，而不在本征值中。只有当存在边缘时，由于体边对应关系（bulk-boundary corresponding），才会有交叉状的边缘态。

            回复

      3. 博主您好！不知道您有没有读过沈顺清老师的这篇Finite Size Effects on Helical Edge States in a Quantum Spin-Hall System，文中直接解析求解了BHZ上自旋部分哈密顿 H(kx, -i∂y)↑，也能根据参数取值的不同得到体态或边缘态，请问他这样做跟您这样做的区别是什么呢？他这样做是不是只能得到费米能附近的两条带而不能像您这样得到多条带？

        回复

        1. BHZ模型是量子自旋霍尔效应。一半的BHZ模型相当于量子反常霍尔效应。
          如果不存在非对角块，那么从能带上来看两者是没什么区别，只是能带条数减半了。

          回复

发表评论

您的电子邮箱地址不会被公开。 必填项已用*标注