打卡11月完成的第二个小任务,仔细研读了汪忠老师这篇Edge States and Topological Invariants of Non-Hermitian Systems文章,将基本内容都重复了一下,最主要的是学习计算了一下非厄SSH模型的拓扑不变量.
{:.info}

前言

非厄密是最近关注度比较高的一个领域,自己最近正好有空,想自己寻找一个合适的方向进一步学习.老实讲虽然现在的导师是做超导的,但是我这两年主要的关注点都在拓扑上面,其实对超导的很多内容也不是很清楚,所以就想反正拓扑掌握的比较多,干脆就在学习学习非厄密中的拓扑,我对这个方向还是非常有兴趣的,觉得它的理论看起来会更漂亮一些.

Non-Hermitian Winding number

我就不重复文章的内容了,这里主要整理一下自己在计算时候遇到的一些问题,因为之前一直在做的是厄密问题,所以在算非厄密的时候多多少少有一些坑,不过我都踩过了.

非厄密与厄密最大的区别就在于这个时候的变量(只从函数的角度看)是$\beta$不再是$k$,所以本征态不再是Bloch形式的拓展态,会变成随着深入体系存在局域在边界上的态,也就是非厄密趋肤效应.我这里就是班门弄斧而已,想搞懂这个问题,还请移步文章.这个时候的$\beta$和原来的$k$之间也会存在着变换关系$\beta=re^{ik}$,由这个$\beta$所确定的就是广义布里渊区(GBZ).所以在计算的时候,也一定是参照这个GBZ进行的.

在非厄密的时候是存在两种基矢的，分别为左矢$\vert u_L\rangle$和右矢$\vert u_R\rangle$.

在对$H(\beta)$进行对角化的时候可以表示为$H(\beta)=TJT^{-1}$,$J$是对角的,每个元素就是本征值,T则是每个本征值所对应的本征矢量

这里特别要注意T中的结果是$H(\beta)$的右矢,而$T^{-1\dagger}$中的则是左矢.所以通常如果使用程序计算矩阵$H(\beta)$的本征矢量和本征值之后,其实你的到的是$T$,所以可以通过求逆后再进行厄密共轭然后得到对应的左矢.这是和厄密情况完全不同的,因为厄密情况下$T$和$T^{-1\dagger}$都是幺正矩阵,它们两个是完全相同的,所以在厄密问题中,我从来没有注意过,页没出现这个情况,在非厄米的时候一定要注意.
{:.warning}

这里$\tilde{u_R}=\sigma_z\rvert u_R\rangle$,最后得到的$Q$是off-diagonal的形式

这里就可以计算非厄密的winding number了

这里在进行数值计算的时候又有坑,(\ref{wind})中的$q^{-1}$就是矩阵$Q$中的$q^{-1}$,公式中的$dq$在数值求和的时候,就是相邻两个点上计算得到的$Q$中的$q$的差值,这样最后就可以得到正确的结果.

png

还有一些其它的结果,比如在$k$时候计算得到的开边界能谱以及非厄密趋肤效应的波函数分布,还有文章中的其它几个图,都进行了重复

png{:width=”330px”,:height=”495px”}png{:width=”330px”,:height=”495px”}

png{:width=”330px”,:height=”495px”}png{:width=”330px”,:height=”495px”}

png{:width=”330px”,:height=”495px”}

代码

能带计算

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
# Import package which is necessary
#import Pkg
#Pkg.add("LinearAlgebra")
#Pkg.add("BenchmarkTools")
#Pkg.add("PyPlot")
using LinearAlgebra,PyPlot,DelimitedFiles
# =========================================================
#--------------------------------
function matrixSet(xn::Int64,t1::Float64,t2::Float64,t3::Float64,gam::Float64)
    ham = zeros(ComplexF64,xn*2,xn*2)
    #----------------------------
    sigx = zeros(Float64,2,2)
    sigx[1,2] = 1.0
    sigx[2,1] = 1.0
    sigy = zeros(ComplexF64,2,2)
    sigy[1,2] = -1im
    sigy[2,1] = 1im
    #-----------------------------
    for k in 0:xn-1
        if k == 0    # First line     
           for m1 in 1:2
                for m2 in 1:2
                    ham[m1,m2] = t1*sigx[m1,m2] + 1im*gam/2.0*sigy[m1,m2]
                    ham[m1,m2 + 2] = (t2 + t3)/2.0*sigx[m1,m2] - 1im*(t2 - t3)/2.0*sigy[m1,m2]  
                end
            end
        elseif k == xn-1
            for m1 in 1:2
                for m2 in 1:2
                    ham[k*2 + m1,k*2 + m2] = t1*sigx[m1,m2] + 1im*gam/2.0*sigy[m1,m2]
                    ham[k*2 + m1,k*2 + m2 - 2] = (t2 + t3)/2.0*sigx[m1,m2] + 1im*(t2 - t3)/2.0*sigy[m1,m2]  
                end
            end
        else
            for m1 in 1:2
                for m2 in 1:2
                    ham[k*2 + m1,k*2 + m2] = t1*sigx[m1,m2] + 1im*gam/2.0*sigy[m1,m2]
                    #   right hopping
                    ham[k*2 + m1,k*2 + m2 + 2] = (t2 + t3)/2.0*sigx[m1,m2] - 1im*(t2 - t3)/2.0*sigy[m1,m2]
                    #   left hopping
                    ham[k*2 + m1,k*2 + m2 - 2] = (t2 + t3)/2.0*sigx[m1,m2] + 1im*(t2 - t3)/2.0*sigy[m1,m2]  
                end
            end
        end
    end
    #----------------------------
    return ham
end
# ======================================================================
function eigHam(k::Float64)
    t1 = 1.0 + 0im
    t2 = 1.0 + 0im
    t3 = 0
    gam = 1.0
   #-------------------
    dx = t1 + (t2 + t3)*cos(k)
    dy = (t2 - t3)*sin(k)
    return sqrt(dx^2 + (dy + 1im*gam/2.0)^2)
end
# ========================================================================
function edgeBand(xn::Int64)
    en::Int64 = 100;
    valSet = zeros(ComplexF64,2*en + 1,xn*2)
    valSetre = zeros(ComplexF64,2*en + 1,xn*2)
    valSetim = zeros(ComplexF64,2*en + 1,xn*2)
    tlist = []
    for m1 in -en:en
        t1 = -3.0*m1/en
        append!(tlist,t1)
        ham = matrixSet(xn,t1,1.0,0.0,4.0/3.0)
        val = eigvals(ham)
        repart = map(real,val)
        impart = map(imag,val)
        val = map(abs,val)
        sort!(val)
        sort!(repart)
        sort!(impart)
        valSetre[m1 + en + 1,:] = repart
        valSetim[m1 + en + 1,:] = impart
        valSet[m1 + en + 1,:] = val
    end
    PyPlot.figure()
    PyPlot.plot(tlist,valSet)
    xlabel("t1")
    ylabel("|E|")
    savefig("absE.png",bbox_inches="tight")
    PyPlot.figure()
    PyPlot.plot(tlist,valSetre)
    xlabel("t1")
    ylabel("Re(E)")
    savefig("ReE.png",bbox_inches="tight")
    PyPlot.figure()
    PyPlot.plot(tlist,valSetim)
    xlabel("t1")
    ylabel("Im(E)")
    savefig("ImE.png",bbox_inches="tight")
    PyPlot.show()
end
# ======================================================================
function wave(xn::Int64)
    ham = matrixSet(xn,1.0,1.0,0.0,4.0/4.0)
    PyPlot.figure()
    val,vec = eigen(ham)
    PyPlot.plot(1:xn,map(abs,vec[1:xn,:]))
    PyPlot.show()
end
# =======================================================================
function main()
    edgeBand(40)
    wave(40)
end
# ========================================================================
main()

Winding number计算

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
function pauli()
    # 构建Pauli矩阵
    sigx = zeros(Float64,2,2)
    sigy = zeros(ComplexF64,2,2)
    sigz = zeros(Float64,2,2)
    sigm = zeros(ComplexF64,2,2)
    sigp = zeros(ComplexF64,2,2)
    #----
    sigx[1,2] = 1.0
    sigx[2,1] = 1.0
    #----
    sigy[1,2] = -1im
    sigy[2,1] = 1im
    #----
    sigz[1,1] = 1
    sigz[2,2] = -1
    #----
    sigm = (sigx - 1im*sigy)/2.0
    sigp = (sigx + 1im*sigy)/2.0
    return sigp,sigm
end
# =========================================================
function hamset(k::Number,tv::Float64)::Matrix{ComplexF64}
    # 哈密顿量构建
    t1::Float64 = tv
    t2::Float64 = 1.0
    gam::Float64 = 4.0/3.0
    r::Float64 = sqrt(abs((t1 - gam/2)/(t1 + gam/2)))
    sigp,sigm = pauli()
    ham = zeros(ComplexF64,2,2)
    ham = (t1 - gam/2.0 + r*exp(1im*k)*t2)*sigm + (t1 + gam/2.0 + 1/r*exp(-1im*k)*t2)*sigp
    return ham
end
# ===========================================================
function Qmat(tt::Float64,tv::Float64)
    occ::Int64 = 1
    sigz = zeros(Float64,2,2)
    sigz[1,1] = 1
    sigz[2,2] = -1
    h1 = hamset(tt,tv)
    vecR = eigvecs(h1)
    vecL = inv(vecR)'
    q1 = vecR[:,occ]
    q11 = sigz*vecR[:,occ]
    q2 = vecL[:,occ]
    q22 = sigz*vecL[:,occ]
    Q = q11*q22' - q1*q2'
    return Q
end
# ================================================================
function winding(tv::Float64)
    re1::ComplexF64 = 0 + 0im
    kn::Int64 = 150
    qmlist = []
    qlist = []
    dq = []
    for k in 0:kn
        k = 2*pi/kn*k
        Q = Qmat(k,tv)
        append!(qlist,Q[1,2])
        append!(qmlist,Q[2,1])
    end
    #--------------------------
    dq = qlist[2:end] - qlist[1:end-1]
    for k in 1:length(qmlist)-1
        re1 += qmlist[k]*dq[k]
    end
    return re1*1im/(2*pi)
end
# ==================================================================
function main()
    tlist = []
    wlist = []
    for tv in -3:0.01:3
        append!(tlist,tv)
        append!(wlist,winding(tv))
    end
    wlist = map(real,wlist)
    #wlist = map(floor,wlist)
    plot(tlist,wlist)
    xlabel("t1")
    ylabel("Winding")
    savefig("TopoInv.png",bbox_inches="tight")
end
# ===============================
main()

公众号

相关内容均会在公众号进行同步，若对该Blog感兴趣，欢迎关注微信公众号。
{:.info}

QR Code

Email

yxliphy@gmail.com