標籤：

Hofstadter模型中的拓撲（Julia）

07-14

這是我為了學院里的量子力學討論班準備的第二期的數值習題內容。

上一次的單體數值zhuanlan.zhihu.com/p/60 里我們講到了轉角石墨烯怎麼算，我在回答zhihu.com/question/3213 中提到了轉角體系里有趣的Hofstadter Butterfly，我們這一節來好好講講Hofstadter Model的模擬。

我會在討論班中推導：

  • 從Landau規範下的磁場中運動的電子，離散化得到Hofstadter模型。
  • 介紹量子相空間，其自身希爾伯特空間的內積可以誘導出一個自然的辛結構（也叫做Berry Curvature）和一個黎曼結構（也叫做Fubini-Study Metric）。
  • 從辛結構出發介紹不變數理論，主要是複線叢的陳數，介紹從U(1)上的和樂計算陳數的Wilson Loop公式以及TKNN公式；叢黎曼結構中介紹不變數理論，主要是歐拉不變數。

我們在這個討論班中僅僅只需要單體量子物理的知識就可以了，但我們依然會使用產生湮滅算符的記號，然而我們可以用單體量子態來等效這些產生湮滅算符。

下面的版本也可以參考我的博客（因為沒有備案，可能國內用戶訪問緩慢physcai.com/2019/05/15/）

Hofstadter Model and its Topology(Julia Version)

About the author:

Jia-Qi Cai

Huazhong University of science and Technology, Wuhan

Email: caidish at uw.edu

#import essential package
using LinearAlgebra, Statistics, Compat, Plots, LaTeXStrings, SparseArrays

In this notebook, we show how to use Julia to simulate quantum Hall effect. We mainly focus on Hofstadter model, the Hamiltonian of which is simply: H = sum_{m,n} - t(c^dagger_{m+1,n}c_{m,n} + c_{m,n+1}^dagger c_{m,n}e^{i 2pi Phi m} + h.c.)

Energy Band

When the flux is commensurable, i.e, Phi = frac{p}{q} , the model have well difined Bloch excitation, which is band electron with q bands. The eigen equations are: -2tcos(k_x + 2pi Phi n)psi_n - t(e^{-ik_y}psi_{n-1} + e^{ik_y}psi_{n+1}) = E_{k_x,k_y}psi_n where n = 1,...,q and the boudary condition is psi_{j+q}= psi_{j} . Solving the linear system above will give the band structure of Hofstadter model

function TorusHamiltonian(kx,ky,t,p,q)::Array{Complex{Float64},2}
Phi = p/q
diagL = -t*exp(-1*im*ky)*ones(q-1)
diagD = complex(-2*t*cos.(kx*ones(q) + 2*pi*Phi*collect(0:1:q-1)))
diagR = -t*exp(+1*im*ky)*ones(q-1)
H = convert(Array,Tridiagonal(diagL,diagD,diagR))
# We use the method Tridiagonal to setup the matrix, and convert it to normal
# matrix to define the boundary condition
H[q,1] = -t*exp(+1*im*ky)
H[1,q] = -t*exp(-1*im*ky)
# The bouandary condition is set.
return H
end;

Lets try to use the Hamiltonian defined in the BZ (geometrically a torus) to generate the band.

p = 1
q = 6
t = 1
kx_step = 100
ky_step = 100
kx_list = LinRange(-pi/q,pi/q,kx_step)
ky_list = LinRange(-pi ,pi ,ky_step)
能帶 = zeros(ky_step,kx_step,q)
for (m,kx) = enumerate(kx_list)
for (n,ky) = enumerate(ky_list)
能帶[n,m,:] = eigvals(TorusHamiltonian(kx,ky,t,p,q))
end
end
#result = permutedims(result,(2,1,3))
能帶 = reshape(能帶,(ky_step,kx_step*q))
plot(ky_list,能帶,legend=false,xlabel = L"k_y",ylabel = L"frac{E}{t}")

I deliberately use Unicode "能帶"(energy band in Chinese) for the result to show the capability of Julia. Please see the winding number protected Dirac fermion near the charge neutral point.

Edge state and Chern number

In this section we show have to calculate the edge state spectrum and the Chern number of the band.

Edge state

Lets see what will happen if we cut the torus into a cylinder, namely only k_y is defined, the eigen equation for this case is: -t(psi_{m+1}(k_y) + psi_{m-1}(k_y)) -2t cos(k_y - 2piPhi m) psi_m(k_y) = Epsi_m(k_y)

where m = 1,...,N and N is the length of the cylinder, with a hard boundary condition psi_0 = psi_{N+1} = 0

function CylindricHamiltonian(ky,N,t,Phi)::Array{Complex{Float64},2}
diagL = -t*ones(N-1)
diagD = -2*t*cos.(ky*ones(N) - 2*pi*Phi*collect(1:1:N))
diagR = -t*ones(N-1)
H = convert(Array,Tridiagonal(diagL,diagD,diagR))
return H
end;
CylindricHamiltonian (generic function with 1 method)
p = 1
q = 6
Phi = p/q
t = 1
ky_step = 100
N = 100
ky_list = range(-pi,stop = pi,length=ky_step)
result = zeros(ky_step,N)
for (m,ky) = enumerate(ky_list)
result[m,:] = eigvals(CylindricHamiltonian(ky,N,t,Phi))
end
plot(ky_list,result,legend=false,xlabel = L"k_y",ylabel = L"frac{E}{t}")

Calculate the Chern number: TKNN method

Here the TKNN formula of Berry curvature (and then Chern number) is used to compute the Chern number.

function ChernNumberByTKNN(H,Hkx,Hky,N,xv,yv,Nx,Ny)
#H is the function of (kx,ky) to give the Hamiltonian.
#Hkx is the kx derivative of H and Hky is the ky derivative of H
#N denotes the number of band. xv,yv define the boundary of the Brioullion zone
#Nx,Ny is the meshes for Chern number
lower_x = minimum(xv)
higher_x = maximum(xv)
lower_y = minimum(yv)
higher_y = maximum(yv)
x_list = collect(range(lower_x,stop = higher_x,length = Nx))
dx = x_list[2] - x_list[1]
y_list = collect(range(lower_y,stop = higher_y,length = Ny))
dy = y_list[2] - y_list[1]
B = zeros(Nx,Ny,N)
C = zeros(N)
for (numX,kx) = enumerate(x_list)
for (numY,ky) = enumerate(y_list)
F = eigen(H(kx,ky))
Wave = F.vectors
Eigen = F.values
p = sortperm(Eigen)
Eigen = Eigen[p]
Wave = Wave[:,p]
Okx = Hkx(kx,ky)
Oky = Hky(kx,ky)
for n = 1:N
for m = 1:N
if m!=n
ketn = Wave[:,n]
ketm = Wave[:,m]
En = Eigen[n]
Em = Eigen[m]
#println(ketn)
#println(ketm)
#println(Okx)
#println(Oky)
#println(imag(ketn*Okx*ketm * ketm*Oky*ketn - ketn*Oky*ketm * ketm*Okx*ketn)/(Em-En)^2)
B[numX,numY,n] += -imag(ketn*Okx*ketm * ketm*Oky*ketn - ketn*Oky*ketm * ketm*Okx*ketn)/(Em-En)^2
end
end
end
end
end
for i = 1:N
C[i] = 1/(2*pi)*sum(B[:,:,i])*dx*dy
end
return C
end;
ChernNumberByTKNN (generic function with 1 method)
function TorusHkx(kx,ky,t,p,q)::Array{Complex{Float64},2}
Phi = p/q
diagL = zeros(Complex,q-1)
diagD = convert(Array{Complex,1},2*t*sin.(kx*ones(q) + 2*pi*Phi*collect(0:1:q-1)))
diagR = zeros(Complex,q-1)
H = convert(Array,Tridiagonal(diagL,diagD,diagR))
# We use the method Tridiagonal to setup the matrix, and convert it to normal
# matrix to define the boundary condition
H[q,1] = 0
H[1,q] = 0
# The bouandary condition is set.
return H
end

function TorusHky(kx,ky,t,p,q)::Array{Complex{Float64},2}
Phi = p/q
diagL = im*t*exp(-1*im*ky)*ones(q-1)
diagD = complex(0*ones(q))
diagR = -im*t*exp(+1*im*ky)*ones(q-1)
H = convert(Array,Tridiagonal(diagL,diagD,diagR))
# We use the method Tridiagonal to setup the matrix, and convert it to normal
# matrix to define the boundary condition
H[q,1] = -im*t*exp(+1*im*ky)
H[1,q] = im*t*exp(-1*im*ky)
# The bouandary condition is set.
return H
end;
p = 1
q = 6
t = 1
kx_step = 100
ky_step = 100
kx = [-pi,pi]
ky = [-pi,pi]
H = (kx,ky)->TorusHamiltonian(kx,ky,t,p,q)
Hkx=(kx,ky)->TorusHkx(kx,ky,t,p,q)
Hky=(kx,ky)->TorusHky(kx,ky,t,p,q)
C = ChernNumberByTKNN(H,Hkx,Hky,q,kx,ky,kx_step,ky_step)/q
for n = 1:q
println("The $n-th Chern Number is: " * string(C[n]))
end
The 1-th Chern Number is: 1.0187478474466722
The 2-th Chern Number is: 1.0445831127776701
The 3-th Chern Number is: -2.0633309602243317
The 4-th Chern Number is: -2.063330960224346
The 5-th Chern Number is: 1.044583112777662
The 6-th Chern Number is: 1.0187478474466733

Now we can see the band chern number is [-1,-1,4,-1,-1] for five band.

Calculate the Chern number: wilson loop method

We then want to calculate the Chern number of each band. The method is to use wilson loop. If you are not familiar with this method, please read the paper: arxiv.org/abs/cond-mat/.

function ChernNumberByWilsonLoop(H,n,xv,yv,Nx,Ny)
#H is the function of (kx,ky) to give the Hamiltonian.
#n denotes n-th band. xv,yv define the boundary of the Brioullion zone
#Nx,Ny is the meshes for Chern number
lower_x = minimum(xv)
higher_x = maximum(xv)
lower_y = minimum(yv)
higher_y = maximum(yv)
x_list = collect(range(lower_x,stop = higher_x,length = Nx))
dx = x_list[2] - x_list[1]
y_list = collect(range(lower_y,stop = higher_y,length = Ny))
dy = y_list[2] - y_list[1]
append!(x_list,[higher_x + dx,higher_x + 2*dx])
append!(y_list,[higher_y + dy,higher_y + 2*dy])

U_x_list = zeros(Complex,Nx+1,Ny+1)
U_y_list = zeros(Complex,Nx+1,Ny+1)

W_list = zeros(Complex,Nx,Ny)
for (num_y,ky) = enumerate(y_list[1:Ny+1])
Wave1 = sortWave(H(x_list[1],ky),n)
for num_x = 1:Nx+1
Wave2 = sortWave(H(x_list[num_x + 1],ky),n)
U_x_list[num_x,num_y] = Wave2*Wave1;
Wave1 = Wave2
end
end
for (num_x,kx) = enumerate(x_list[1:Nx+1])
Wave1 = sortWave(H(kx,y_list[1]),n)
for num_y = 1:Ny+1
Wave2 = sortWave(H(kx,y_list[num_y + 1]),n)
U_y_list[num_x,num_y] = Wave2*Wave1;
Wave1 = Wave2
end
end
for num_x = 1:Nx
for num_y = 1:Ny
W_list[num_x,num_y] = U_x_list[num_x,num_y]*U_y_list[num_x+1,num_y]*conj(U_x_list[num_x,num_y+1])*conj(U_y_list[num_x,num_y])
end
end
B = angle.(W_list)
return 1/(2*pi)*sum(B[:])
end

function sortWave(M,n)
F = eigen(M)
Wave = F.vectors
Eigen = F.values
p = sortperm(Eigen)
Eigen = Eigen[p]
Wave = Wave[:,p]
return Wave[:,n]
end
sortWave (generic function with 1 method)
p = 1
q = 6
t = 1
kx_step = 100
ky_step = 100
kx = [-pi,pi]
ky = [-pi,pi]
H = (kx,ky)->TorusHamiltonian(kx,ky,t,p,q)
for n = 1:q
println("The $n-th Chern Number is: " * string(ChernNumberByWilsonLoop(H,n,kx,ky,kx_step,ky_step)/q))
end
The 1-th Chern Number is: -1.018781432213558
The 2-th Chern Number is: -1.0437354586033607
The 3-th Chern Number is: 5.062544010724829
The 4-th Chern Number is: 5.06254401072483
The 5-th Chern Number is: -1.0437354586033607
The 6-th Chern Number is: -1.018781432213558

Finite size wave function and edge state

So now weve already known that there is nontrivial topology and its corresponding edge state, we may compute a real finite size lattice and see the wave funciton in the topological gap. The single particle Hamiltonian is: H = sum_{m,n}-t (left|m+1,n ight angle leftlangle m,n ight|+e^{i2piPhi m}left|m,n+1 ight angle leftlangle m,n ight|+ h.c.)

where left|m,n ight angle = left|m ight angle otimes left|n ight angle and otimes denotes the kronecker tensor.

# we first test the indexing of the kronecker

a = [1 0 0 0]#ket{1}
b = [0 0 1 0]#ket{3}
c = kron(a,b)#ket{3,1}
function index(m,n,M,N)
index = m + M*(n-1)
end
println(c[index(3,1,4,4)])

a1 = [0 1 0 0]#ket{2}
b1 = [0 0 0 1]#ket{4}
c1 = kron(a1,b1)#ket{4,2}
println(c1[index(4,2,4,4)])

Mat = c1*c #bra{4,2}*ket{3,1}
println(Mat[index(4,2,4,4),index(3,1,4,4)])
1
1
1

The code below generates Hilbert basis via dense matrix manipulation, which is SUPER SLOW and should never be used.

#function that give the ket{i,j}
function basis_dense(i,j,M,N)
res = zeros(Complex,M*N,1)
res[index(i,j,M,N),1] = 1
return res
end

# function that give the bra{m,n}*ket{i,j}
function mat_dense(m,n,i,j,M,N)
res = zeros(Complex,M*N,M*N)
res[index(m,n,M,N),index(i,j,M,N)] = 1
return res
end
mat_dense (generic function with 1 method)

we may use the sparse matrix to give the Hamiltonian and the matrix

println(sparse([3,2],[3,2],[1.0im,2im],5,5))
println(spzeros(Complex{Float64}, 5,5))
[2, 2] = 0.0+2.0im
[3, 3] = 0.0+1.0im
5×5 SparseMatrixCSC{Complex{Float64},Int64} with 0 stored entries
function basis_sp(i,j,M,N)
res = sparse([index(i,j,M,N)],[1],[1 + 0.0im],M*N,1)
end

function mat_sp(m,n,i,j,M,N)
res = sparse([index(m,n,M,N)],[index(i,j,M,N)],[1],M*N,M*N)
end
mat_sp (generic function with 1 method)

Now we now what is sparse matrix, lets construct our Hamiltonian by sort all the indexes and the value.

function FiniteHamiltonian(t,p,q,M,N)::SparseMatrixCSC{Complex{Float64},Int64}
dim = M*N
row_list = []
col_list = []
val_list = []
for m = 1:M-1
for n = 1:N
append!(row_list,index(m+1,n,M,N))
append!(col_list,index(m,n,M,N))
append!(val_list,-t*(1+0.0im))

append!(row_list,index(m,n,M,N))
append!(col_list,index(m+1,n,M,N))
append!(val_list,-t*(1+0.0im))
end
end
for m = 1:M
for n = 1:N-1
append!(row_list,index(m,n+1,M,N))
append!(col_list,index(m,n,M,N))
append!(val_list,-t*exp(1im*2*pi*p/q*m))

append!(row_list,index(m,n,M,N))
append!(col_list,index(m,n+1,M,N))
append!(val_list,-t*exp(-1im*2*pi*p/q*m))
end
end
H = sparse(row_list,col_list,val_list,M*N,M*N)
return H
end
FiniteHamiltonian (generic function with 1 method)
t = 1
p = 1
q = 3
M = 30
N = 31
H = FiniteHamiltonian(t,p,q,M,N)
F = eigen(Array(H))
Wave = F.vectors
Eigen = F.values
perm = sortperm(Eigen)
Eigen = Eigen[perm]
Wave = Wave[:,perm]
p1 = plot(Eigen,marker =(:c,1),xlabel = "Number of Eigenvalues",ylabel = L"frac{E}{t}")

Phi = p/q
ky_step = 100
Nky = 100
ky_list = range(-pi,stop = pi,length=ky_step)
result = zeros(ky_step,Nky)
for (m,ky) = enumerate(ky_list)
result[m,:] = eigvals(CylindricHamiltonian(ky,Nky,t,Phi))
end
p2 = plot(ky_list,result,legend=false,xlabel = L"k_y")

plot(p1,p2,legend = false,size=(800,300))

We found that there exists "in-gap" state in finite size eigen value(left) comparing to edge state in cylindric(right). Are they really edge states? We may plot those in gap state.

edgewave = Wave[:,(Eigen.>1.2).&(Eigen.<1.8)]
stateToPlot = edgewave[:,1]
NumToPlot = norm.(stateToPlot)
NumToPlot = reshape(NumToPlot,(M,N))
heatmap(NumToPlot)

Its exactly an edge state!!!!

練習1：分別使用上面提供的TorusHamiltonian，CylindricHamiltonian，和FiniteHamiltonian來計算Hofstadter Butterfly譜，也就是在不同的 Phi 下體系的能譜圖。

  • 對於Torus，令 q = 137,p=1,...,137 ，在布里淵區內選取 400 個點，對角化 H(k_x,k_y,frac{p}{q }) 壓縮成一個 frac{p}{q} 對應的譜 E(frac{p}{q}) 。
  • 對於CylindricHamiltonian，修改它的定義，使之直接接受 Phi 作為參數，取圓柱長度為 N=100 , k_y 採樣100個點，對角化 H(k_y,frac{p}{q })壓縮成 E(frac{p}{q})
  • 對於FiniteHamiltonian，修改它的定義，使之直接接受 Phi 作為參數，對角化 H(frac{p}{q }) 得到 E(frac{p}{q})

練習 2***，計算 H(k_x,k_y) 對應的量子度規 g_{xx},{g_{xy}},g_{yy} ，求出其歐拉不變數，它的物理意義是什麼？

參考文獻:

1、Topological Insulators and Topological Superconductors by Bernevig&Hughes

2、Artificial Gauge Fields with Ultracold Atoms in Optical Lattices, PhD thesis of Monika Aidelsburge

3、journals.jps.jp/doi/10.

4、sciencedirect.com/scien, XG Wen&A.Zee

推薦閱讀：

※AM227 物理中的數值方法
※玻璃隕石成因的幾種繆談，不要讓價格迷失雙眼
※通過聲波精確列印液滴
※詳解離心機的基本原理、類型及操作應用
※常溫狀態水的表面張力是多少？

TAG:量子物理 | 量子多體理論 | 物理學 |