学术源码, 模型和能带

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

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_1x方向的跃迁项 (hopping),H_2y方向的跃迁项 (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 次浏览

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

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

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

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

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

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

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

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

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

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

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

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

发表评论

您的电子邮箱地址不会被公开。