拓扑不变量, 学术

SSH模型中的Wilson loop(附Python代码)

这是之前的一篇关于SSH模型的博文:SSH模型的哈密顿量、能带图和卷绕数(附Python代码)。除了卷绕数(winding number)这个不变量,我们还可以计算Wilson loop,进而得到电极化(polarization),或称为电偶极矩(dipole moment)。本篇内容主要来源于文献[1]。

这里写出SSH模型哈密顿量:

h(k)=\begin{pmatrix}0 & \gamma+\lambda e^{-ik} \\\gamma+\lambda e^{ik} & 0\end{pmatrix}

该哈密顿量满足空间反演对称性(inversion symmetry):

\hat{\mathcal{I}}h(k)\hat{\mathcal{I}}^{-1}=h(-k)

其中,空间反演对称性算符为\hat{\mathcal{I}}=\sigma_1

同时也满足手征对称性(chiral symmetry):

\hat{\Pi}h(k)\hat{\Pi}^{-1}=-h(k)

其中,手征对称性算符为\hat{\Pi}=\sigma_3

可证明:在空间反演对称性或手征对称性下,偶极矩均可量子化到0或者\frac{1}{2}(mod 1)。证明见参考资料[1]。

SSH模型同时满足以上这两种对称性,因此偶极矩是量子化的。当|\gamma|>|\lambda|时,极化p=0;当|\gamma|<|\lambda|时,极化p=\frac{1}{2}

如果在SSH模型的哈密顿量中加入\sigma_3项,会同时打破模型的空间反演对称性和手征对称性,这时候偶极矩不再量子化。

说明1:这里添加了固定波函数规范的代码,在本篇博文后面给出了更好的代码,不需要对波函数进行规范。

说明2:固定规范函数可通过开源项目Guan进行调用:https://py.guanjihuan.com。此外,固定规范的方法介绍也可参考这篇:非简并波函数和简并波函数的固定规范

一、当|\gamma|<|\lambda|时,代码如下:

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

import numpy as np
import cmath
from math import *


def hamiltonian(k):  # SSH模型哈密顿量
    gamma = 0.5
    lambda0 = 1
    h = np.zeros((2, 2), dtype=complex)
    h[0,0] = 0
    h[1,1] = 0
    h[0,1] = gamma+lambda0*cmath.exp(-1j*k)
    h[1,0] = gamma+lambda0*cmath.exp(1j*k)
    return h


def main():
    Num_k = 100
    k_array = np.linspace(-pi, pi, Num_k)
    vector_array = []
    for k in k_array:
        vector  = get_occupied_bands_vectors(k, hamiltonian)   
        vector_array.append(vector)
        # vector_array.append(vector*cmath.exp(1j*np.random.uniform(0, pi))) # 随机相位测试

    # 波函数固定一个规范
    vector_sum = 0
    for i0 in range(Num_k):
        vector_sum += np.abs(vector_array[i0])
    index = np.argmax(np.abs(vector_sum))
    for i0 in range(Num_k):
        vector_array[i0] = find_vector_with_fixed_gauge_by_making_one_component_real(vector_array[i0], index=index)

    # # 波函数固定一个规范
    # import guan
    # vector_array = guan.find_vector_array_with_fixed_gauge_by_making_one_component_real(vector_array, precision=0.005)

    # 计算Wilson loop
    W_k = 1
    for i0 in range(Num_k-1):
        F = np.dot(vector_array[i0+1].transpose().conj(), vector_array[i0])
        W_k = np.dot(F, W_k)
    nu = np.log(W_k)/2/pi/1j
    # if np.real(nu) < 0:
    #     nu += 1
    print('p=', nu, '\n')
    

def get_occupied_bands_vectors(x, matrix):  
    matrix0 = matrix(x)
    eigenvalue, eigenvector = np.linalg.eig(matrix0) 
    vector = eigenvector[:, np.argsort(np.real(eigenvalue))[0]] 
    return vector


def find_vector_with_fixed_gauge_by_making_one_component_real(vector, precision=0.005, index=None):
    if index == None:
        index = np.argmax(np.abs(vector))
    sign_pre = np.sign(np.imag(vector[index]))
    for phase in np.arange(0, 2*np.pi, precision):
        sign =  np.sign(np.imag(vector[index]*cmath.exp(1j*phase)))
        if np.abs(np.imag(vector[index]*cmath.exp(1j*phase))) < 1e-9 or sign == -sign_pre:
            break
        sign_pre = sign
    vector = vector*cmath.exp(1j*phase)
    if np.real(vector[index]) < 0:
        vector = -vector
    return vector


if __name__ == '__main__':
    main()

计算结果为:

p= (-0.5+0.009257772639491353j)

二、当|\gamma|<|\lambda|时:

def hamiltonian(k):  # SSH模型哈密顿量
    gamma = 1.5
    lambda0 = 1
    h = np.zeros((2, 2), dtype=complex)
    h[0,0] = 0
    h[1,1] = 0
    h[0,1] = gamma+lambda0*cmath.exp(-1j*k)
    h[1,0] = gamma+lambda0*cmath.exp(1j*k)
    return h

计算结果为:

p= (-6.308806949364591e-17+0.003169518547424918j)

三、在SSH模型的哈密顿量中加入\sigma_3项:

def hamiltonian(k):
    gamma = 0.5
    lambda0 = 1
    delta = 0.3
    h = np.zeros((2, 2), dtype=complex)
    h[0,0] = delta
    h[1,1] = -delta
    h[0,1] = gamma+lambda0*cmath.exp(-1j*k)
    h[1,0] = gamma+lambda0*cmath.exp(1j*k)
    return h

计算结果为:

p= (0.328981519117191+0.007898978069732571j)

由评论区Kabu Xiong提供的资料[5],指出了Wilson loop对相位是不依赖的,是gauge invariant,这是因为在Wilson loop 中,波函数都出现了两次,而且是相互复数共轭,相位刚好相消。

然而在程序中,我们把vector_array.append(vector)修改成vector_array.append(vector*cmath.exp(1j*np.random.uniform(0, pi))),也就是每个波函数加一个随机的相位。这时候只有加了波函数规范的代码,结果才是正确的。

原因在于-pipi两个点的波函数应该是要相位相同。如果两个点是独立算的,那么相位是不能完全相消的。因此,如果在程序中不重新算pi点的波函数,直接引用-pi点的波函数,这时候是不需要对波函数进行规范。

代码如下(建议采用这种更好方式):

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

import numpy as np
import cmath
from math import *


def hamiltonian(k):
    gamma = 0.5
    lambda0 = 1
    delta = 0
    h = np.zeros((2, 2), dtype=complex)
    h[0,0] = delta
    h[1,1] = -delta
    h[0,1] = gamma+lambda0*cmath.exp(-1j*k)
    h[1,0] = gamma+lambda0*cmath.exp(1j*k)
    return h


def main():
    Num_k = 100
    k_array = np.linspace(-pi, pi, Num_k)
    vector_array = []
    for k in k_array:
        vector  = get_occupied_bands_vectors(k, hamiltonian)   
        if k != pi:
            vector_array.append(vector)
        else:
            vector_array.append(vector_array[0])

    # 计算Wilson loop
    W_k = 1
    for i0 in range(Num_k-1):
        F = np.dot(vector_array[i0+1].transpose().conj(), vector_array[i0])
        W_k = np.dot(F, W_k)
    nu = np.log(W_k)/2/pi/1j
    # if np.real(nu) < 0:
    #     nu += 1
    print('p=', nu, '\n')
    

def get_occupied_bands_vectors(x, matrix):  
    matrix0 = matrix(x)
    eigenvalue, eigenvector = np.linalg.eig(matrix0) 
    vector = eigenvector[:, np.argsort(np.real(eigenvalue))[0]] 
    return vector


if __name__ == '__main__':
    main()

这时,即使把vector_array.append(vector)修改成vector_array.append(vector*cmath.exp(1j*np.random.uniform(0, pi))),也就是每个波函数加一个随机的相位,即使没有波函数规范的代码结果都是正确的。

参考资料:

[1] Wladimir A. Benalcazar, B. Andrei Bernevig, and Taylor L. Hughes, “Electric multipole moments, topological multipole moment pumping, and chiral hinge states in crystalline insulators", Phys. Rev. B 96, 245115 (2017): P8-P13

[2] S. Franca, J. van den Brink, and I. C. Fulga, "An anomalous higher-order topological insulator", Phys. Rev. B 98, 201114(R) (2018).

[3] 凝聚态中的拓扑(四):从TKNN到Z2拓扑绝缘体

[4] Band topology in classical waves: Wilson-loop approach to topological numbers and fragile topology

[5] 和YZZ同学的讨论

4,130 次浏览

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

评论说明:
(1)在保留浏览器缓存的前提下,目前支持72小时自主修改或删除个人评论。如果自己无法修改或删除评论,可再次评论或联系我。如有发现广告留言,请勿点击链接,博主会不定期删除。
(2)评论支持Latex公式。把latexpage作为标签放在任何位置,评论中的公式可正常编译,示例:
$Latex formula$  [latexpage]

32 thoughts on “SSH模型中的Wilson loop(附Python代码)”

  1. 感谢博主的分享,想把代码改成matlab学习编程思路。
    vector = eigenvector[:, np.argsort(np.real(eigenvalue))[0]] 这句代码的意思是要提取每次计算的特征向量吗?不明白。每次计算都会有两个特征向量。

      1. clear
        clc
        close all

        k0=linspace(-pi,pi,100);
        for i=1:100
        str=zeros(2,1);
        k=k0(i);
        h1=HH(k);
        [V,D]=eig(h1);
        vector=V(:,1); %取能量低的一条带
        if k~=pi
        VV(:,i)=vector;
        else
        VV(:,i)=V(:,1);
        end
        end
        W_k=1;
        for j=1:99
        F=dot(VV(:,j+1)',VV(:,j));
        ff(:,j)=F;
        W_k=dot(F,W_k);
        end
        nu=log(W_k)/(2*pi*1i);

        function h=HH(k)
        gamma=0.5;
        lambda0=1;
        delta=0;
        h=zeros(2);
        h(1,1)=delta;
        h(2,2)=-delta;
        h(1,2)=gamma+lambda0*exp(-1i*k);
        h(2,1)=gamma+lambda0*exp(1i*k);
        end

        结果没对,不知道是哪里出了问题
        W_k=1.32804668806497e-33 + 5.66626118539546e-49i;
        nu=6.79052539801457e-17 + 12.0482836587981i;

        1. 在matlab中,dot(A, B)已经包含了对A的厄密共轭, 即:dot(A, B)=A'*B。所以不需要额外再加一撇。

  2. 对于计算wilson loop那一块,初始化为1,然后做i0循环。 我想请问 计算下来的F是一个数还是一个矩阵,就涉及到F后面的两个波函数内积的情况,是行矢量点乘列矢量还是列矢量点乘行矢量。 后面是调用np.dot我理解上是矩阵乘法,但是初始化的是一个数。所以我想请问,F是矩阵还是数。 这个问题我可以用博主的代码自己打开看看,先在这里问一下。

  3. 老师您好,请问对于一些简并能带计算Berry Curvature的时候,需要考虑non-Abelian的Berry Curvature,从定义上虽然能做但会不会波函数的随机相位在数值上会有问题呢?这方面最近有数值方面的工作吗

    1. 简并的时候比较麻烦,如果需要波函数连续,可能是要对简并的波函数基做幺正变换。目前我还没做过测试。

  4. 博主您好,请问为什么把本征矢乘以一个,能使模最大分量变实数的一个相因子,就可以达到统一规范的目的呢?(从你的代码中看出来似乎是这个逻辑)

    1. 因为波函数乘exp(i*phase)都是方程的解,所以每一个k点的波函数计算可能存在不同的规范。把波函数某个分量变成实数,可以起到固定规范的作用,从而达到波函数连续,这是因为乘一个exp(i*phase)都会导致该分量为非实数。而至于为什么取最大分量,是为了避开波函数分量为0的情况,从而使得代码运行不容易出错。
      该方法可以参照文献“An anomalous higher-order topological insulator - An anomalous higher-order topological insulator”的补充材料中的说明。
      另外再强调下:Wilson loop的计算可以不考虑波函数的规范任意性(前提是-pi和pi的位置用同一个波函数)。

      1. 这个分量在各个波函数中也不必是同一个位置的是吧?谢谢您的回复和提供的文献,我再学习下这篇文献。

    2. 局部使得波函数最大分量变实数可能还不一定使得波函数连续,需要考虑整个布里渊区,目前代码已经完成补充修改。我猜是存在两个k点的最大分量处于不同位置的情况,尤其是临界点附近(即存在两个最大分量),所以需要全程固定一个指标(index),虽然旧代码在这个案例中没出现错误。
      但这种方案仍存在一点问题,需要小心存在局部k点该分量为0左右。即使该分量整个布里渊区是最大的,但不排除某些k点为0+或0-。对于这种情况,需要提高相位查找的精度来解决。

  5. 请问这个波函数的规范是什么含义呢?
    之前计算一个周期晶格的不变量,跟wilson loop非常像,计算过程相同,但是波函数用的是bloch函数的周期部分(即补上相位exp(ikr)以后的部分),积分路径是k空间上的若干个小平行四边形,结果吻合得很好。中间好像没有涉及到对波函数的规范,所以不太明白规范的具体意义

    1. 不变量和可观测量都是和规范的无关的(gauge independent)。波函数相位的选取主要是在数值计算过程中会产生,例如求导时相邻的两个波函数需要连续。这里计算wilson loop也不需要选取规范,前提是把-pi和pi两个点选取为同一个波函数(这两个点本身就是同一个状态)。

    1. 如果带不相交或者简并,照样算就可以了。如果存在多个带相交或者简并,会麻烦一些,目前我还没考虑。

    1. 感谢提供的资料。Wilson loop是规范不变的。代码稍微调整下,是不需要对波函数进行规范,博文已经更新。

  6. 作者您好!我想请问一下,如何用python循环出SSH模型的实空间中的哈密顿?谢谢!

  7. 作者你好,我想问一下,“ # 波函数固定一个规范”那一步是为什么呀,我用你的代码把这一步去掉后,也能得到相同的结果。如果不规范会有什么问题吗

    1. 如果把vector_array.append(vector)修改成vector_array.append(vector*cmath.exp(1j*np.random.uniform(0, 1))),也就是每个波函数加一个随机的相位。这时候只有加了波函数规范的代码,结果才是正确的。

      1. 哦哦,明白些了。还有一个问题,就是我想问一下那个波函数规范的公式原型是怎样的呀,我从代码上分析不出来,也没有看到相关的资料说明。谢谢!

        1. 目前没看到什么公式,只要使得波函数连续即可。可以自己写代码进行判断。

          固定一个规范的目的也是使得波函数连续。这里代码选取的规范是使波函数的最大分量(不为零的分量)的虚部为零。也可以选取其他规范。

    2. 博文更新了。因为wilson loop是gauge invariant,所以适当调整下代码,是可以不需要对波函数进行规范。

发表回复

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