Real Space Hamiltonian to k Space Hamiltonian
一维实空间哈密顿量
实空间哈密顿量
对于一维最理想的模型，其哈密顿量为:
\[H = \sum_{i}\varepsilon_i c_i^{\dagger}c_i + \sum_{<i,j>}t_{ij}c_i^{\dagger}c_j + h.c.
\]
若只考虑最近邻近似，考虑最简单的情况: \(\varepsilon_i = a, t_{i,i+1}=b\), 写成矩阵形式为:
\[H = \begin{pmatrix}
a & b & & & & &\\
b & a & b & & & &\\
& b & a & b & & &\\
& & \cdots & \cdots & \cdots & &\\
& & & b & a & b & \\
& & & & b & a & b\\
& & & & & b & a\\
\end{pmatrix} _{n\times n}
\]
很明显，该哈密顿量的解析解为:
\[E= a+ 2b*\cos(k), k=\frac{i*\pi}{n+1}, i =1,2,3,...
\]
能谱图为:
考虑 \(\varepsilon_i = a, t_{i,i+1}=b, t_{i+1,i}=c\), 写成矩阵形式为:
\[H = \begin{pmatrix}
a & b & & & & &\\
c & a & b & & & &\\
& c & a & b & & &\\
& & \cdots & \cdots & \cdots & &\\
& & & c & a & b & \\
& & & & c & a & b\\
& & & & & c & a\\
\end{pmatrix} _{n\times n}
\]
该哈密顿量的解析解为:
\[E= a+ 2\sqrt{bc}*\cos(k), k=\frac{i*\pi}{n+1}, i =1,2,3,...
\]
其能谱图为:
当 N 增加时，能谱图为:
\(k\) 空间哈密顿量
将上述实空间矩阵变换到 \(k\) 空间：
\[H(k)=U^{\dagger}HU
\]
矩阵 \(U\) 的形式为:
\[U = \sqrt{\frac{2}{n+1}}\sin(k_p*r_m)\\
k_p = \frac{p\pi}{n+1}. (p = 1,2,3,..,n)\\
r_m = m. (m = 1,2,3,..,n)
\]
这里的 \(p\) 是实空间指标，\(m\) 是 \(k\) 空间指标.
当 \(b=c\) 时, 实空间矩阵变换后的 \(k\) 空间矩阵是对角矩阵, 对角元即为本征值.
当 \(b\neq c\)时, \(k\) 空间矩阵不是对角矩阵.
a = 0, b = 2, c =3. 不难发现 k 空间矩阵并不是完全的对角矩阵.
二维实空间哈密顿量
考虑简单的二维哈密顿量
\[H = \frac{\hbar^2}{2m}(k_x^2+k_y^2)
\]
应用有限差分方法, \(k_x = -i\partial_x, k_y = -i\partial_y\), \(x,y\) 方向网格密度为 \(N_x, N_y\).
\(t_x = t_y = \frac{\hbar^2}{2ma^2}\), 实空间哈密顿量可写为
\[H = \begin{pmatrix}
H00 & H01 & & & \\
H10 & H00 & H01 & & \\
& \cdots & \cdots & \cdots &\\
& & H10 & H00 & H01\\
& & & H10 & H00\\
\end{pmatrix}_{N_x*N_y\times N_x*N_y}
\]
其中
\[H00 = \begin{pmatrix}
2t_x + 2t_y & -t_x & & & \\
-t_x &2t_x + 2t_y & -t_x & & \\
& \cdots & \cdots & \cdots &\\
& & -t_x &2t_x + 2t_y & -t_x\\
& & & -t_x &2t_x + 2t_y \\
\end{pmatrix}_{N_x\times N_x},\\
H01 = \begin{pmatrix}
-t_y & & & & \\
&-t_y && & \\
& & \cdots &&\\
& & &-t_y& \\
& & &&-t_y\\
\end{pmatrix}_{N_x\times N_x}
\]
将实空间哈密顿量变换到 \(k\) 空间
\[H_k = U^{\dagger}HU\\
U = \frac{2}{\sqrt{(N_x+1)*(N_y+1)}}\sin(k_L*r_M)
\]
这里的[1]
\[\sin(k_L*r_M) = \sin(k_p x_m)\sin(k_q y_n)\\
k_L = (\frac{p\pi}{N_x+1},\frac{q\pi}{N_y+1}). (p = 1,2,3,..,N_x, q = 1,2,3,..,N_y)\\
r_M = (x_m,y_n). (m = 1,2,3,..,N_x, n = 1,2,3,..,N_y)\\
M = (m-1)N_y + n
\]
详细参数见附件python文件, 下图为放大后的 HK 矩阵元.
参考文献
[1] Mincheol Shin. Full-quantum simulation of hole transport and band-to-band tunneling in nanowires using the kp method. J. Appl. Phys. 106, 054505 (2009).
[2] https://physics.stackexchange.com/questions/674136/numerically-transforming-hamiltonian-into-k-space
附件
一维哈密顿量, 点击查看代码
import numpy as np
import cmath
import time
from matplotlib import pyplot as plt
'''
For any complex a, b, c
[a, b, 0, 0, 0, 0, 0 ]
[c, a, b, 0, 0, 0, 0 ]
[0, c, a, b, 0, 0, 0 ]
[0, 0, c, a, b, 0, 0 ]
..
[0, 0, 0, 0, c, a, b ]
[0, 0, 0, 0, 0, c, a ]
'''
def Ham_R_space(a,b,c,N):
H = np.zeros((N,N),complex)
for j in range(N):
H[j][j] = a
if j + 1 < N:
H[j][j+1] = b
H[j+1][j] = c
# H[0][N-1] = c
# H[N-1][0] = b
return H
def Ham_k_space(a,b,c,N):
HRspace = Ham_R_space(a=a,b=b,c=c,N=N)
p = np.linspace(1,N,N)
y = np.linspace(1,N,N)
U = np.empty((N,N), complex)
for i in range(N):
U[i][:] = (2/(N+1))**0.5*np.sin(p[i]*np.pi/(N+1)*y)
Hkspace = np.dot(np.dot(U.conjugate().T,HRspace),U)
return Hkspace
def analytic_solution(a,b,c,N):
res = np.zeros(N,complex)
for i in range(N):
k = (i+1)*np.pi/(N+1)
res[i] = a + 2*cmath.sqrt(b*c)*np.cos(k)
# print('解析解\n',res)
return res
def main(a,b,c,N):
p = np.linspace(1,N,N)
HRspace = Ham_R_space(a=a,b=b,c=c,N=N)
plt.matshow(HRspace.real)
plt.colorbar()
eig_R = np.linalg.eigvals(HRspace)
# print('严格解:\n',sorted(eig_R))
Hkspace = Ham_k_space(a=a,b=b,c=c,N=N)
plt.matshow(Hkspace.real)
plt.colorbar()
eig_K = np.linalg.eigvals(Hkspace)
plt.figure()
plt.plot(np.diagonal(Hkspace).real,linewidth=7, label = 'diag(Hkspace)')
plt.plot(sorted(eig_R,reverse=True),linewidth=4, label = 'eig(HRspace)')
plt.plot(sorted(eig_K,reverse=True),linewidth=2, label = 'eig(Hkspace)')
res = analytic_solution(a=a,b=b,c=c,N=N)
# print('\n解析解:\n',sorted(res))
plt.plot(res,'.',label = 'analytical solution of HRspace')
plt.xlabel('$k$')
plt.ylabel('Energy')
plt.legend(loc='best')
plt.show()
if __name__ == '__main__':
main(a=0,b=2,c=3,N=200)
二维哈密顿量, 点击查看代码
import numpy as np
import time,cmath
from matplotlib import pyplot as plt
'''
[[A, B, 0, 0, 0, 0, 0]
[B*, A, B, 0, 0, 0, 0]
[0, B*, A, B, 0, 0, 0]
[0, 0, B*, A, B, 0, 0]
...
[0, 0, 0, 0, B*, A, B]
[0, 0, 0, 0, 0, B*, A]]
A = B =
[[x1, x2, 0, 0, 0, 0, 0] [[x4, x5, 0, 0, 0, 0, 0]
[x3, x1, x2, 0, 0, 0, 0] [x6, x4, x5, 0, 0, 0, 0]
[ 0, x3, x1, x2, 0, 0, 0] [ 0, x6, x4, x5, 0, 0, 0]
[ 0, 0, x3, x1, x2, 0, 0] [ 0, 0, x6, x4, x5, 0, 0]
... ...
[ 0, 0, 0, 0, x3, x1, x2] [ 0, 0, 0, 0, x6, x4, x5]
[ 0, 0, 0, 0, 0, x3, x1]] [ 0, 0, 0, 0, 0, x6, x4]]
'''
## 解析解为
def ham_R_space(x1,x2,x3,x4,x5,x6,Nx=4,Ny=4):
H = np.zeros((Nx*Ny,Nx*Ny),dtype="complex")
# H = [[ 0 for i in range(Nx*Ny)] for i in range(Nx*Ny)]
for i in range(Nx):
for j in range(Ny):
index = Nx*j+i
## x1
H[index][index] = x1
## x4
## (i, j) --> (i, j+1)
if j + 1 < Ny:
H[index][index+Nx] = x4
H[index+Nx][index] = x4.conjugate()
## x2 matrix
## (i, j) --> (i+1, j)
if i + 1 < Nx:
H[index][index+1] = x2
H[index+1][index] = x3
## x5
## (i, j) --> (i+1, j+1)
if i + 1 < Nx and j + 1 < Ny:
H[index][index+Nx+1] = x5
H[index+Nx+1][index] = x5.conjugate()
## x6
## (i, j) --> (i+1, j-1)
if i > 0 and j + 1 < Ny:
H[index][index+Nx-1] = x6
H[index+Nx-1][index] = x6.conjugate()
eig = np.linalg.eigvals(H)
# print('严格解:')
# print(sorted(eig))
return H
def ham_K_space(x1,x2,x3,x4,x5,x6,Nx=4,Ny=4):
HRspace = ham_R_space(x1=x1,x2=x2,x3=x3,x4=x4,x5=x5,x6=x6,Nx=Nx,Ny=Ny)
U = np.zeros((Nx*Ny,Nx*Ny),dtype="complex")
# H = [[ 0 for i in range(Nx*Ny)] for i in range(Nx*Ny)]
for p in range(Nx):
for q in range(Ny):
index1 = Nx*q+p
for i in range(Nx):
for j in range(Ny):
index2 = Nx*j+i
kx = (p+1)*np.pi/(Nx+1)
ky = (q+1)*np.pi/(Ny+1)
U[index1][index2] = 2/((Nx+1)*(Ny+1))**0.5*np.sin(kx*(i+1))*np.sin(ky*(j+1))
# test = np.dot(U.conjugate().T,U)
# print('U test: ',test)
Hkspace = np.dot(np.dot(U.conjugate().T,HRspace),U)
return Hkspace
def hamiltonian_analyze1(x1,x2,x3,x4,x5,x6,Nx=4,Ny=4):
### 该解析方法只适用于 x4, x5, x6 均为实数，且 x5 = x6
res = []
for i in range(1,Nx+1):
for j in range(1, Ny+1):
kx = i*np.pi/(Nx+1)
ky = j*np.pi/(Ny+1)
H00 = x1 + x4*np.exp(1j*ky) + x4.conjugate()*np.exp(-1j*ky)
H01 = x2 + x5*np.exp(1j*ky) + x6.conjugate()*np.exp(-1j*ky)
H10 = x3 + x6*np.exp(1j*ky) + x5.conjugate()*np.exp(-1j*ky)
tmp = H00 + 2*cmath.sqrt(H01*H10)*np.cos(kx)
res.append(tmp)
print('\n解析解:\n',sorted(np.array(res)))
def hamiltonian_analyze2(x1,x2,x3,x4,x5,x6,Nx=4,Ny=4):
### 该解析方法只适用于 x2=x3, 且 x5=x6, 且不能为 负数
### 或者在 x5 = x6 = 0 的情况下，该解析解有效!
res = []
for i in range(1,Nx+1):
for j in range(1, Ny+1):
kx = i*np.pi/(Nx+1)
ky = j*np.pi/(Ny+1)
H00 = x1 + 2*cmath.sqrt(x2*x3)*np.cos(kx)
H01 = x4 + 2*cmath.sqrt(x5*x6)*np.cos(kx)
H10 = x4.conjugate() + 2*cmath.sqrt(x5.conjugate()*x6.conjugate())*np.cos(kx)
tmp = H00 + 2*cmath.sqrt(H01*H10)*np.cos(ky)
res.append(tmp)
print('\n解析解:\n',sorted(np.array(res)))
def main():
x1 = 1
x2 = 2
x3 = 2
x4 = 1
x5 = 2
x6 = 2
H = ham_R_space(x1=x1,x2=x2,x3=x3,x4=x4,x5=x5,x6=x6,Nx=20,Ny=20)
HK = ham_K_space(x1=x1,x2=x2,x3=x3,x4=x4,x5=x5,x6=x6,Nx=20,Ny=20)
plt.matshow(HK.real)
plt.colorbar()
plt.matshow(H.real)
plt.colorbar()
# hamiltonian_analyze1(x1=x1,x2=x2,x3=x3,x4=x4,x5=x5,x6=x6,Nx=4,Ny=4)
# hamiltonian_analyze2(x1=x1,x2=x2,x3=x3,x4=x4,x5=x5,x6=x6,Nx=4,Ny=4)
plt.show()
main()
[2023/10/11] 补充
实空间哈密顿量矩阵 \(H_R\) 通过 \(U\) 矩阵变换到 \(K\) 空间哈密顿量矩阵\(H_K\), 下面这两种结果并不完全一致。
1), \(H_K = U^{\dagger}H_RU\)
通过U矩阵展开相乘得到 \(H_K\) 矩阵
2), \(H_K = k_p^2\delta_{p,p'}\)
假设网格点无限多，实空间是连续函数，然后对x积分得到的结果。
原因是，第一种情况计算出来的本征值形式为 \(2t+2t\cos(k_p)\), 第二种情况计算出来的是 \(tk_p^2\), 这两种结果只在 \(k_p\) 很小时近似成立。
浙公网安备 33010602011771号