模型和能带, 科学

BBH高阶拓扑绝缘体模型(附Python代码)

BBH模型(Benalcazar-Bernevig-Hughes Model)的参考文献为:Wladimir A. Benalcazar, B. Andrei Bernevig, and Taylor L. Hughes, "Quantized electric multipole insulators", Science Volume 357:61-66 (2017).

模型:

\begin{aligned}  h=&[\gamma+\lambda cos(k_x)]\Gamma_4+\lambda sin(k_x)\Gamma_3\\ &+[\gamma+\lambda cos(k_y)]\Gamma_2+\lambda sin(k_y)\Gamma_1+\delta \Gamma_0 \end{aligned}

其中,

\Gamma_0=\tau_3 \sigma_0=\begin{pmatrix}1 & 0 & & \\0 & 1 & & \\& & -1 & 0 \\& & 0 & -1\end{pmatrix}

\Gamma_1=-\tau_2\sigma_1=-\begin{pmatrix}& & 0 & -i \\& & -i & 0 \\0 & i & & \\i & 0 & &\end{pmatrix}

\Gamma_2=-\tau_2\sigma_2=-\begin{pmatrix}& & 0 & -1 \\& & 1 & 0 \\0 & 1 & & \\-1 & 0 & &\end{pmatrix}

\Gamma_3=-\tau_2\sigma_3=-\begin{pmatrix}& & -i & 0 \\& & 0 & i \\i & 0 & & \\0 & -i & &\end{pmatrix}

\Gamma_4=\tau_1\sigma_0=\begin{pmatrix}& & 1 & 0 \\& & 0 & 1 \\1 & 0 & & \\0 & 1 & &\end{pmatrix}

以上用到泡利矩阵的张量积:泡利矩阵以及泡利矩阵的张量积

所以哈密顿量展开写为:

\begin{aligned}h&=\begin{pmatrix}\delta  &  0  & \gamma+\lambda cos(k_x)+i\lambda sin(k_x) & \gamma+\lambda cos(k_y)+i \lambda sin(k_y)  \\0   &   \delta  & -[\gamma+\lambda cos(k_y)]+i \lambda sin(k_y) & \gamma+\lambda cos(k_x)-i \lambda sin(k_x) \\\gamma+\lambda cos(k_x)-i \lambda sin(k_x) & -[\gamma+\lambda cos(k_y)]-i \lambda sin(k_y) & -\delta & 0 \\\gamma+\lambda cos(k_y)-i \lambda sin(k_y) & \gamma+\lambda cos(k_x)+i \lambda sin(k_x) & 0 & -\delta\end{pmatrix}\\&=\begin{pmatrix}\delta & 0 & \gamma+\lambda exp(ik_x) & \gamma+\lambda exp(ik_y) \\0 & \delta & -\gamma-\lambda exp(-i k_y) & \gamma+\lambda exp(-i k_x) \\\gamma+\lambda exp(-i k_x) & -\gamma-\lambda exp(i k_y) & -\delta & 0 \\\gamma+\lambda exp(-ik_y) & \gamma+\lambda exp(i k_x) & 0 & -\delta\end{pmatrix}\end{aligned}

文献中给出的对应示意图:

元胞为(1,2,3,4)。把哈密顿量反傅里叶变换后,可以知道:

  • 元胞内部:3到1的跃迁项为\gamma;2到4的跃迁项为\gamma;4到1的跃迁项为\gamma;2到3的跃迁项为-\gamma
  • 元胞之间:1到3的跃迁项为\lambda;4到2的跃迁项为\lambda;1到4的跃迁项为\lambda;3到2的跃迁项为-\lambda

即:元胞内部跃迁项绝对值为\gamma,元胞间的跃迁项绝对值为\lambda。此外,2和3之间的跃迁要加个负号。

计算局域态密度(LDOS)的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/8557
"""

import numpy as np
from math import *


def hamiltonian(Nx, Ny): 
    delta = 1e-3
    gamma = 1e-3
    lambda0 = 1
    h = np.zeros((4*Nx*Ny, 4*Nx*Ny))
    # 元胞内部跃迁
    for x in range(Nx):
        for y in range(Ny):
            h[x*Ny*4+y*4+0, x*Ny*4+y*4+0] = delta
            h[x*Ny*4+y*4+1, x*Ny*4+y*4+1] = delta
            h[x*Ny*4+y*4+2, x*Ny*4+y*4+2] = -delta
            h[x*Ny*4+y*4+3, x*Ny*4+y*4+3] = -delta

            h[x*Ny*4+y*4+0, x*Ny*4+y*4+2] = gamma
            h[x*Ny*4+y*4+0, x*Ny*4+y*4+3] = gamma
            h[x*Ny*4+y*4+1, x*Ny*4+y*4+2] = -gamma
            h[x*Ny*4+y*4+1, x*Ny*4+y*4+3] = gamma
            h[x*Ny*4+y*4+2, x*Ny*4+y*4+0] = gamma
            h[x*Ny*4+y*4+2, x*Ny*4+y*4+1] = -gamma
            h[x*Ny*4+y*4+3, x*Ny*4+y*4+0] = gamma
            h[x*Ny*4+y*4+3, x*Ny*4+y*4+1] = gamma

    # y方向上的元胞间跃迁
    for x in range(Nx):
        for y in range(Ny-1):
            h[x*Ny*4+y*4+0, x*Ny*4+(y+1)*4+3] = lambda0
            h[x*Ny*4+(y+1)*4+1, x*Ny*4+y*4+2] = -lambda0
            h[x*Ny*4+y*4+2, x*Ny*4+(y+1)*4+1] = -lambda0
            h[x*Ny*4+(y+1)*4+3, x*Ny*4+y*4+0] = lambda0

    # x方向上的元胞间跃迁
    for x in range(Nx-1):
        for y in range(Ny):
            h[x*Ny*4+y*4+0, (x+1)*Ny*4+y*4+2] = lambda0
            h[(x+1)*Ny*4+y*4+1, x*Ny*4+y*4+3] = lambda0
            h[(x+1)*Ny*4+y*4+2, x*Ny*4+y*4+0] = lambda0
            h[x*Ny*4+y*4+3, (x+1)*Ny*4+y*4+1] = lambda0
    return h
    

def main():
    Nx = 10
    Ny = 10
    fermi_energy = 0
    h = hamiltonian(Nx, Ny)
    green = np.linalg.inv((fermi_energy+1e-6j)*np.eye(h.shape[0])-h)
    
    x_array = []
    y_array = []
    DOS = []
    for x in range(Nx):
        for y in range(Ny):
            x_array.append(x*2+2)
            y_array.append(y*2+2)
            DOS.append(-np.imag(green[x*Ny*4+y*4+0, x*Ny*4+y*4+0])/pi)

            x_array.append(x*2+1)
            y_array.append(y*2+1)
            DOS.append(-np.imag(green[x*Ny*4+y*4+1, x*Ny*4+y*4+1])/pi)

            x_array.append(x*2+1)
            y_array.append(y*2+2)
            DOS.append(-np.imag(green[x*Ny*4+y*4+2, x*Ny*4+y*4+2])/pi)

            x_array.append(x*2+2)
            y_array.append(y*2+1)
            DOS.append(-np.imag(green[x*Ny*4+y*4+3, x*Ny*4+y*4+3])/pi)
    DOS = DOS/np.sum(DOS)
    Plot_2D_Scatter(x_array, y_array, DOS, xlabel='x', ylabel='y', title='BBH Model', filename='BBH Model')


def Plot_2D_Scatter(x, y, value, xlabel='x', ylabel='y', title='title', filename='a'): 
    import matplotlib.pyplot as plt
    fig = plt.figure()
    ax = fig.add_subplot(111)
    plt.subplots_adjust(bottom=0.2, right=0.8, left=0.2) 
    for i in range(np.array(x).shape[0]):
        ax.scatter(x[i], y[i], marker='o', s=1000*value[i], c=[(1,0,0)])
    ax.set_title(title, fontsize=20, fontfamily='Times New Roman')
    ax.set_xlabel(xlabel, fontsize=20, fontfamily='Times New Roman') 
    ax.set_ylabel(ylabel, fontsize=20, fontfamily='Times New Roman') 
    ax.tick_params(labelsize=15)  # 设置刻度值字体大小
    labels = ax.get_xticklabels() + ax.get_yticklabels() 
    [label.set_fontname('Times New Roman') for label in labels] # 设置刻度值字体
    # plt.savefig(filename+'.jpg', dpi=300) 
    plt.show()
    plt.close('all')


if __name__ == '__main__':
    main()

计算结果为:

换一种画图方式(Plot 3D Surface):

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

import numpy as np
from math import *


def hamiltonian(Nx, Ny): 
    delta = 1e-3
    gamma = 1e-3
    lambda0 = 1
    h = np.zeros((4*Nx*Ny, 4*Nx*Ny))
    # 元胞内部跃迁
    for x in range(Nx):
        for y in range(Ny):
            h[x*Ny*4+y*4+0, x*Ny*4+y*4+0] = delta
            h[x*Ny*4+y*4+1, x*Ny*4+y*4+1] = delta
            h[x*Ny*4+y*4+2, x*Ny*4+y*4+2] = -delta
            h[x*Ny*4+y*4+3, x*Ny*4+y*4+3] = -delta

            h[x*Ny*4+y*4+0, x*Ny*4+y*4+2] = gamma
            h[x*Ny*4+y*4+0, x*Ny*4+y*4+3] = gamma
            h[x*Ny*4+y*4+1, x*Ny*4+y*4+2] = -gamma
            h[x*Ny*4+y*4+1, x*Ny*4+y*4+3] = gamma
            h[x*Ny*4+y*4+2, x*Ny*4+y*4+0] = gamma
            h[x*Ny*4+y*4+2, x*Ny*4+y*4+1] = -gamma
            h[x*Ny*4+y*4+3, x*Ny*4+y*4+0] = gamma
            h[x*Ny*4+y*4+3, x*Ny*4+y*4+1] = gamma

    # y方向上的元胞间跃迁
    for x in range(Nx):
        for y in range(Ny-1):
            h[x*Ny*4+y*4+0, x*Ny*4+(y+1)*4+3] = lambda0
            h[x*Ny*4+(y+1)*4+1, x*Ny*4+y*4+2] = -lambda0
            h[x*Ny*4+y*4+2, x*Ny*4+(y+1)*4+1] = -lambda0
            h[x*Ny*4+(y+1)*4+3, x*Ny*4+y*4+0] = lambda0

    # x方向上的元胞间跃迁
    for x in range(Nx-1):
        for y in range(Ny):
            h[x*Ny*4+y*4+0, (x+1)*Ny*4+y*4+2] = lambda0
            h[(x+1)*Ny*4+y*4+1, x*Ny*4+y*4+3] = lambda0
            h[(x+1)*Ny*4+y*4+2, x*Ny*4+y*4+0] = lambda0
            h[x*Ny*4+y*4+3, (x+1)*Ny*4+y*4+1] = lambda0
    return h
    

def main():
    Nx = 10
    Ny = 10
    fermi_energy = 0
    h = hamiltonian(Nx, Ny)
    green = np.linalg.inv((fermi_energy+1e-6j)*np.eye(h.shape[0])-h)
    DOS = np.zeros((Ny*2, Nx*2))
    for x in range(Nx):
        for y in range(Ny):
            DOS[y*2+1, x*2+1] = -np.imag(green[x*Ny*4+y*4+0, x*Ny*4+y*4+0])/pi
            DOS[y*2+0, x*2+0] = -np.imag(green[x*Ny*4+y*4+1, x*Ny*4+y*4+1])/pi
            DOS[y*2+1, x*2+0] = -np.imag(green[x*Ny*4+y*4+2, x*Ny*4+y*4+2])/pi
            DOS[y*2+0, x*2+1] = -np.imag(green[x*Ny*4+y*4+3, x*Ny*4+y*4+3])/pi
    DOS = DOS/np.sum(DOS)
    Plot_3D_Surface(np.arange(1, 2*Nx+0.001), np.arange(1, 2*Ny+0.001), DOS, xlabel='x', ylabel='y', zlabel='DOS', title='BBH Model', filename='BBH Model')


def Plot_3D_Surface(x, y, matrix, xlabel='x', ylabel='y', zlabel='z', title='title', filename='a'): 
    import matplotlib.pyplot as plt
    from matplotlib import cm
    from matplotlib.ticker import LinearLocator
    fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
    plt.subplots_adjust(bottom=0.1, right=0.65) 
    x, y = np.meshgrid(x, y)
    if len(matrix.shape) == 2:
        surf = ax.plot_surface(x, y, matrix, cmap=cm.coolwarm, linewidth=0, antialiased=False) 
    elif len(matrix.shape) == 3:
        for i0 in range(matrix.shape[2]):
            surf = ax.plot_surface(x, y, matrix[:,:,i0], cmap=cm.coolwarm, linewidth=0, antialiased=False) 
    ax.set_title(title, fontsize=20, fontfamily='Times New Roman')
    ax.set_xlabel(xlabel, fontsize=20, fontfamily='Times New Roman') 
    ax.set_ylabel(ylabel, fontsize=20, fontfamily='Times New Roman') 
    ax.set_zlabel(zlabel, fontsize=20, fontfamily='Times New Roman') 
    # ax.set_zlim(-1, 1)  # 设置z轴的范围
    ax.zaxis.set_major_locator(LinearLocator(2)) # 设置z轴主刻度的个数
    ax.zaxis.set_major_formatter('{x:.2f}')   # 设置z轴主刻度的格式
    ax.tick_params(labelsize=15)  # 设置刻度值字体大小
    labels = ax.get_xticklabels() + ax.get_yticklabels() + ax.get_zticklabels()
    [label.set_fontname('Times New Roman') for label in labels] # 设置刻度值字体
    cax = plt.axes([0.80, 0.15, 0.05, 0.75]) # color bar的位置 [左,下,宽度, 高度]
    cbar = fig.colorbar(surf, cax=cax)  # color bar
    cbar.ax.tick_params(labelsize=15) # 设置color bar刻度的字体大小
    for l in cbar.ax.yaxis.get_ticklabels():
        l.set_family('Times New Roman')
    # plt.savefig(filename+'.jpg', dpi=300) 
    plt.show()
    plt.close('all')


if __name__ == '__main__':
    main()

计算结果为:

2,408 次浏览

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

15 thoughts on “BBH高阶拓扑绝缘体模型(附Python代码)”

    1. 是不是通过热力学量平均值定义从而和Green function联系上的那种?如Abrikosov的P96,(10.25)

      1. 我也不是很了解,我猜电荷密度可能是需要对态密度进行能量上的积分。如果考虑有限温度,也是要考虑热力学平均。可以查查教材或文献,看有没有这个概念的定义。

  1. 这个局部态密度LDOS的物理含义是什么呀?
    看程序里面是对(λ*I-H)求逆矩阵,然后将逆矩阵的每个元素按照原子排列分配给2*Nx*2*Ny个原子,最后对这些原子绘图
    之前看到有个对二维弹性板(固体声学)做逆向设计的文章,也提到了LDOS,在那篇文章里面是把归一化后每个点的面外位移的平方定义为了LDOS
    那么是不是弹性力学里面的结构矩阵也能定义成这个形式?

    1. LDOS=local density of states,表示的是态密度在实空间的分布。波函数模的平方也能体现出电子在空间上的分布情况,归一化后是等价的,参考数值验证“波函数模平方分布”和“格林函数计算出的态密度分布”的关系(附Python代码)

      你说的弹性力学中的LDOS我不是很了解,估计也类似于波函数模的平方吧。是否能写出这种形式需要具体考虑或推导,或者查阅相关的文献。

  2. 冒昧打扰一下楼主,想请问一下高阶拓扑相是怎么定义的?您可否一下定义高阶的文章?

  3. 楼主,您好!我有一个一直困扰的问题没有解决,想请教一下您:
    在算矩阵本征值时,有时候python的结果误差比较大(和Fortran比较),得到的结果是不准确的。初步判断应该是精度不够的问题。想请问您,这该如何解决?

    1. 本征值应该都是一样的才对。看在python中矩阵是否定义为复数,例如np.zeros((x,y), dtype=’complex’)。

      在内容一样的情况下,检查是很容易的。在某个语句前后分别输出结果,一步步对比验证。

  4. 冒昧打扰一下楼主,请问如果把四极矩的晶格形状正方形变成正六边形,那么角态会还集中在六个角吗,哈密顿量该怎么改写呢?

    1. 对称性不一样,具体的形式和结果我也不确定。可能有这方面的文献吧,可以找一找看。

      1. 哦好的,谢谢您,我有查到了相关的论文,“Higher-Order Topological Insulator in Twisted Bilayer Graphene”,讲述了石墨烯晶格中的情况,角态和这个类似。
        还想再请问您,这个BBH模型中它后面讲述的八极矩的程序该怎么在您的四极矩的基础上改进,我看您第二段程序是换了一种画图方式,内容还是四极矩的。

        1. 实空间的哈密顿量需要重新写。八级矩一般是要画三维的分布图,可以参考:使用matplotlib画图。计算输出的数据格式不唯一,根据画图程序的要求调整。

发表评论

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