学术源码, 拓扑不变量

陈数Chern number的计算(Kubo公式,附Python代码)

另外几篇关于计算陈数的方法:

本篇通过Kubo公式计算陈数。这个方法跟高效法、Wilson loop方法一样,不需要对波函数求导,因此数值计算过程不依赖于波函数相位。

从定义出发,贝里联络(Berry connection)公式为:

A_n(\bm{R})=i\langle n(\bm{R})|\frac{\partial}{\partial \bm{R}}|n(\bm{R})\rangle

贝里曲率(Berry curvature)公式:

\Omega_n(\bm{R})= \nabla_ {\bm{R}} \times  A_n(\bm{R})

展开写为[1]:

\begin{aligned}\Omega_{\mu\nu}^{n}(\bm{R}) &=\frac{\partial}{\partial R_\mu} A_\nu^n(\bm{R})- \frac{\partial}{\partial R_\nu} A_\mu^n(\bm{R}) \\&= i[\frac{\partial}{\partial R_{\mu}}\langle n(\bm{R})|\frac{\partial}{\partial R_{\nu}}|n(\bm{R})\rangle-(\nu  \leftrightarrow \mu)]  \end{aligned}

进一步推导(由“向同学”提供):

继续推导需要用到以下等式:

\langle n'|\frac{\partial H}{\partial \bm{R}}|n\rangle = \langle n'| \frac{\partial}{\partial \bm{R}} |n\rangle(E_n-E_{n'}) for n'\neq n

该等式的证明如下[2]。

薛定谔方程:H|n\rangle =E_n|n\rangle

两边求导:\frac{\partial H}{\partial \bm{R}}|n\rangle+H\frac{ \partial |n\rangle}{ \partial \bm{R} }=\frac{\partial E}{\partial  \bm{R} }|n\rangle+E_n\frac{  \partial |n\rangle }{ \partial  \bm{R} }

两边左乘\langle n'|\langle n'| \frac{\partial H}{\partial \bm{R}}|n\rangle+ \langle n'| H\frac{ \partial |n\rangle}{ \partial \bm{R} }= \langle n'| \frac{\partial E_n}{\partial  \bm{R} }|n\rangle+ \langle n'| E_n\frac{  \partial |n\rangle }{ \partial  \bm{R} }

整理后: \langle n'| \frac{\partial H}{\partial \bm{R}}|n\rangle+ E_{n'}\langle n'|\frac{ \partial }{ \partial \bm{R} }|n\rangle=E_n \langle n'|\frac{  \partial  }{ \partial  \bm{R} }|n\rangle 。其中,考虑了 n'\neq n ,因此\langle n'| n \rangle = 0

继续整理,得到: \langle n'|\frac{\partial H}{\partial \bm{R}}|n\rangle = \langle n'| \frac{\partial}{\partial \bm{R}} |n\rangle(E_n-E_{n'}) 。证毕!

接下来继续推导贝里曲率公式,考虑贝里曲率中的第一项。

中间插入单位算符 \sum\limits_{n'} |n'\rangle \langle n'|,得到: \sum\limits_{n'} \langle \frac{\partial}{\partial R_{\mu}} n| n'\rangle \langle n'| \frac{\partial}{\partial R_{\nu}}|n\rangle。由于贝里曲率后面还有一项(\nu  \leftrightarrow \mu) ,对于n' =n的情况两项相减为零,因此这一项可以直接写为: \sum\limits_{n'\neq n} \langle \frac{\partial}{\partial R_{\mu}} n| n'\rangle \langle n'| \frac{\partial}{\partial R_{\nu}}|n\rangle

代入上面被证明过的表达式,得到: \sum\limits_{n'\neq n} \frac{ \langle n|\frac{\partial H}{\partial R_{\mu}}|n'\rangle   \langle n'|\frac{\partial H}{\partial R_{\nu}}|n\rangle  }{(E_n-E_{n'})^2}

贝里曲率公式:

\Omega_{\mu\nu}^{n}(\bm{R}) = i \sum\limits_{n'\neq n} \frac{ \langle n|\frac{\partial H}{\partial R_{\mu}}|n'\rangle   \langle n'|\frac{\partial H}{\partial R_{\nu}}|n\rangle-(\nu  \leftrightarrow \mu)  }{(E_n-E_{n'})^2}

以这篇“陈数Chern number的计算(定义法,附Python/Matlab代码)”中的哈密顿量为例子,计算陈数,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/16148
"""

import numpy as np
from math import *
import time


def hamiltonian(kx, ky):  # one QAH model with Chern number = 2
    t1 = 1.0
    t2 = 1.0
    t3 = 0.5
    m = -1.0
    matrix = np.zeros((2, 2), dtype=complex)
    matrix[0, 1] = 2*t1*cos(kx)-1j*2*t1*cos(ky)
    matrix[1, 0] = 2*t1*cos(kx)+1j*2*t1*cos(ky)
    matrix[0, 0] = m+2*t3*sin(kx)+2*t3*sin(ky)+2*t2*cos(kx+ky)
    matrix[1, 1] = -(m+2*t3*sin(kx)+2*t3*sin(ky)+2*t2*cos(kx+ky))
    return matrix


def main():
    start_time = time.time()
    n = 200       # integration
    delta = 1e-6  # derivation
    chern_number = 0 
    for kx in np.arange(-pi, pi, 2*pi/n):
        for ky in np.arange(-pi, pi,2*pi/n):
            H = hamiltonian(kx, ky)
            eigenvalue, eigenvector = np.linalg.eig(H)
            vector_0 = eigenvector[:, np.argsort(np.real(eigenvalue))[0]]
            vector_1 = eigenvector[:, np.argsort(np.real(eigenvalue))[1]]
            eigenvalue = np.sort(np.real(eigenvalue))

            H_delta_kx = hamiltonian(kx+delta, ky)-hamiltonian(kx, ky) 
            H_delta_ky = hamiltonian(kx, ky+delta)-hamiltonian(kx, ky)

            berry_curvature = 1j*(np.dot(np.dot(np.dot(np.dot(np.dot(vector_0.transpose().conj(), H_delta_kx/delta), vector_1), vector_1.transpose().conj()), H_delta_ky/delta), vector_0)- np.dot(np.dot(np.dot(np.dot(np.dot(vector_0.transpose().conj(), H_delta_ky/delta), vector_1), vector_1.transpose().conj()), H_delta_kx/delta), vector_0))/(eigenvalue[0]-eigenvalue[1])**2

            chern_number = chern_number + berry_curvature*(2*pi/n)**2
    chern_number = chern_number/(2*pi)
    print('Chern number = ', chern_number)
    end_time = time.time()
    print('运行时间(min)=', (end_time-start_time)/60)


if __name__ == '__main__':
    main()

计算结果:Chern number = (-1.9999999999798186+1.9718230813235673e-19j)

参考资料:

[1] Quantized Hall Conductance in a Two-Dimensional Periodic Potential

[2] Berry phase effects on electronic properties

[3] https://phyx.readthedocs.io/en/latest/TI/Lecture%20notes/2.html

[4] https://mp.weixin.qq.com/s/VJxeLOIyBoJQ1-x-P4Pwtg

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

2,277 次浏览

16 thoughts on “陈数Chern number的计算(Kubo公式,附Python代码)”

  1. 写了个matlab版的,结果对的上:

    %kubo公式,特点是不用加波函数规范
    function []=main_kubo_chen()
    tic
    n=200;%布里渊区积分离散参数
    delta=1e-6;%求导步长,应该是要比积分步长小
    c_num=0;%陈数初始化
    band_index=1;%算哪条陈数的能带,本征值从低到高
    N=2;%哈密顿量矩阵的阶数,也是波函数个数,根据实际情况写
    kx0=-pi:2*pi/n:pi;
    ky0=kx0;%在第一布里渊区积分
    
    for ii=1:length(kx0)
        for jj=1:length(ky0)
            kx=kx0(ii);
            ky=ky0(jj);
            H=hami(kx,ky);
            [eigen_vecor,eigen_value]=eig(H);
            
             [e_value,index]=sort(real(diag(eigen_value)));
             wave_fun_matrix=eigen_vecor(:,index);%按照本征值从小到大排列波函数
             
             D_H_kx = (hami(kx+delta, ky)-hami(kx, ky))/delta;%矩阵,求导
             D_H_ky = (hami(kx, ky+delta)-hami(kx, ky))/delta;%矩阵
             
             berry_cur=0;%贝利曲率初始化
             
             wave_fun0=wave_fun_matrix(:,band_index);%|n>
             value0=e_value(band_index);%En
             for uu=1:N
                 if uu==band_index
                     continue
                 end
                 
                 temp_wave_fun=wave_fun_matrix(:,uu);%|n'>
                 temp_value=e_value(uu);%En'
                 
                 a1=wave_fun0'*D_H_kx*temp_wave_fun*temp_wave_fun'*D_H_ky*wave_fun0;
                 a2=wave_fun0'*D_H_ky*temp_wave_fun*temp_wave_fun'*D_H_kx*wave_fun0;
                 
                 berry_cur=berry_cur+1j*(a1-a2)/(value0-temp_value)^(2);
             end
             %贝利曲率计算完毕,berry_cur
             c_num=c_num+berry_cur*((2*pi/n)^2);
        end
    end
    
    c_num=c_num/(2*pi);
    
    disp(['第',num2str(band_index),'条能带陈数为:',num2str(c_num)])
    toc
    end
    
    
    
    
    function  [out]=hami(kx,ky)
    %一个测试哈密顿量
    t1 = 1.0;
    t2 = 1.0;
    t3 = 0.5;
    m = -1.0;
    matrix =zeros(2);
    matrix(1, 2) = 2*t1*cos(kx)-1j*2*t1*cos(ky);
    matrix(2,1) = 2*t1*cos(kx)+1j*2*t1*cos(ky);
    matrix(1,1) = m+2*t3*sin(kx)+2*t3*sin(ky)+2*t2*cos(kx+ky);
    matrix(2,2) = -(m+2*t3*sin(kx)+2*t3*sin(ky)+2*t2*cos(kx+ky));
    out=matrix;
    end

        1. 这个久保公式法在有能带兼并的时候咋办,En-En’=0,发散了,数值是缺失值NAN

          1. 我可以先试着加微扰打开带隙吧。或者把价带波函数当成整体来算?这个我还没验证过。

    1. 你好,这个matlab代码可以画3维的贝里曲率图吗?z轴是贝里曲率,x、y轴分别是kx、ky那种,谢谢~

        1. 可否具体出一个画3维贝里曲率图形的matlab代码,我修改了一番这个评论区的求陈数的代码,没有画出来贝里曲率图,感谢指导~

  2. 好像计算陈数还有个kubo公式的方法,威尔逊环等方法,我见过的就有四五种。

发表评论

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