拓扑不变量, 学术

BHZ模型的自旋陈数和Z2不变量(附Python、Matlab代码)

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

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)]

一、 自旋陈数

因为自旋向上和自旋向下是两个相互独立的子体系(对应的陈数刚好相反),自旋是守恒的,因此可以定义自旋陈数。

陈数的计算方法为:陈数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

[4] Quantum Spin Hall Effect in Three Dimensional Materials: Lattice Computation of Z2 Topological Invariants and Its Application to Bi and Sb

[5] First-principles calculation of Z2topological invariants within the FP-LAPW formalism

4,884 次浏览

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

27 thoughts on “BHZ模型的自旋陈数和Z2不变量(附Python、Matlab代码)”

  1. 你好,请问一下对于BHZ模型,如何构造时间反演算符T,使得满足时间反演对称性TH_kT^{-1}=H_{-k}^*且T^2=-1

  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
    那么这个偏函数里求得dd,应该是8X8的矩阵的行列式吗

  3. [eigenvalue, index]=sort(diag(eigenvalue), 'ascend');
    vector_new_2 = vector(:, index(2));
    vector_new_1 = vector(:, index(1));
    请问楼主上面的代码是取了导带来计算吗,取价带3,4来算的话会不会Z2是相反数啊

    1. 取的是价带。ascend是上升的意思,最前面两个带是价带。另外两条带我没考虑过,你可以算算。理论上如果考虑所有的带,Z2应该等于0。

      1. 请问这里取的两条带必须要是不同自旋贡献的兼并的两条带吗,要算Z2是不是模型必须要有这个自由度啊

        1. 如果算自旋陈数,那么需要两个自旋是独立的。算Z2没有这个要求。
          至于Z2是否必须考虑自旋,目前我看到的都是有考虑,且需要满足时间反演对称性。其他情况我不是很了解,没深入研究过,你可以查下文献看看。

  4. 作者可以出一篇运用wannier function center的计算Z2不变量的文章吗?然后比较说明一下和用贝里曲率算法的比较说明,这两个方法的对比,看了文献,并不是非常的理解。谢谢

    1. 这方面的文章我看的也不多,目前不是很了解。之后有空的时候可能会考虑。你可以多找些文献看看,看不同文献中的描述。

      1. 嗯嗯,是的。对于简单的BHZ正方晶格的积分选取是比较显而易见的。那如果对于非正方晶格比如三角晶格或六角晶格,就会出现缺角啊,还是说根据话高对称点之间的能带图的取法来选取呢?

      1. 十分感谢博主,对于一个初次接触拓扑的研究生帮助甚大,希望可以多交流。

  5. 博主您好,我想请问一下,对于这样的二能带系统,里面的A具体是怎么定义的,能展开说说吗?因为对于二能带,A(联络)应该是矩阵形式,那么对应于U,U的具体形式是怎么样的呢?看程序看的不太明显…麻烦了

      1. 能带有两个指标很明显A的分量为A11,A12,A21,A22,数字为能带指标啊…为啥不是矩阵呢

        1. 是有形成一个矩阵,不同能带上的波函数是正交的,非对角线上的矩阵元(内积值)应该为零。U最后是取了一个行列式,具体公式你可以参照参考资料[2]的PPT,或者其他参考资料。我对这方面也没有仔细研究过,只是简单算了下。

          1. 参考文献3中的6式,内积的两个波函数分别是两个相似变换矩阵(2n列向量排列形成),两个方阵相乘再求行列式,
            应该是这样不,感觉和代码里不一样

发表评论

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

Captcha Code