对角化费米子体系和玻色子体系的异同（初步回答）
本问题集来自身边计算势力脑的自发涌现。此问题搞凝聚态物理计算的容易遇到或想到，本科初学凝聚态导论的过程中也会去思考，不是第一次发现此问题了，也不是第一次思考这一问题了。总之，我想这是一个新手普遍会遇到的问题，而非一次钻牛角尖式的发问。作者没有系统学习过量子多体理论，本文的一些说法一定是存在某些问题的，请辨证阅读！
类似问题：https://www.zhihu.com/question/450788784
文题中“初步回答”的含义是，仅仅在本文把问题表述出来，作者没有做到清晰表述，所以程度上只能归类到初步。
一、问题集描述（把问题清晰表述出来，往往就是答案本身）
请看以下有关联的模糊表述：
1. 规定玻戈留玻夫变换后的新算符和原算符的玻色费米型一致（线性组合并不能保证此条件，比如费米算符的线性组合不一定依然满足费米反对易关系，需要额外规定）；
2. 可以直接对角化费米子体系（在粒子数表象下）获取能谱；
3. “不能”直接对角化玻色子体系（在粒子数表象下）获取能谱；
请移步https://jsci.cnu.edu.cn/HTMLfile/2021/1/202101003.html
4. 产生湮灭算符的玻色对易关系及费米反对易关系与波函数交换对称及交换反对称的关联；
4i. *无关内容：波函数交换反对称与泡利不相容原理的等同；
5. 任意N个费米子体系的波函数（交换反对称），可用Slatter行列式构筑（实际上是Dirac先提出），行列式维度N；
请移步https://blog.sciencenet.cn/blog-100379-869543.html
6. 对于N个格点，M个费米子(M<N)占据这N个格点，在粒子数表象下有C_N^M个基底，基底可写为类似Slatter行列式的形式，也可以用费米子产生算符作用在|0,0,...,0>上表达，再利用费米子的产生湮灭算符的反对易关系，可以构造出交换反对称的性质。可以用这样表述的基底夹住哈密顿量，写出费米子体系哈密顿量的所有矩阵元，之后可直接对角化；
7. 能够构造一个粒子数不守恒的基底，由二项式定理，知2^N=C_N^0+C_N^1+...+C_N^N，也可以由排列组合问题的基本方法，知每一个格点只有两种可能（0态或1态），故N个格点共2^N个态（无相互作用费米子体系的直积态表示）；
8. 无相互作用玻色子天然适配“单粒子基底”的直积态（线性空间，对易），而费米子体系则需要加入额外的结构（暂且自行命名为行列式直积态），来构造出交换反对称的特性（线性空间+某种规则或称结构，反对易），比如从粒子数表象01到自旋表象↑↓的约当-魏格纳变换就可以实现这一额外引入的结构；
请移步https://www.zhihu.com/question/450788784
9. ……
二、用Jordan-Wigner变换对角化费米子体系哈密顿量
以下图片内容及附录的程序来自紧扣的dagger的通讯作者，C_{jl}^\ddagger，他是真正的double dagger，也是问题的提供人，没有他，本篇专栏将不会被写成。
他计算了一个4格点的仅考虑最近邻hoping的体系，体系粒子数不守恒体系（每个格点至多激发1个费米子，故总粒子数至多4个）。
结果如下：
当然也可继续计算6格点的能谱，会看到tight-binding的余弦曲线特征，随格点增加逐渐显露：
附录：C_{jl}^\ddagger的程序源码
经double dagger同意，源码可以贴出。源码没有经过任何优化，且某些关键行被注释掉了，需要自行恢复。
'''
ED for spinless free fermion model in Fock space
'''
import numpy as np
from copy import deepcopy
a = np.array([[0,0],[1,0]])
a_dag = np.array([[0,1],[0,0]])
op_n = a_dag@a
I = np.array([[1,0],[0,1]])
F = np.array([[1,0],[0,-1]])
def FreeFermion_ED(N):
"""
[ Note ]
For fermions, the wavefunction should be antisymmetric under particle exchange,
therefore the Hilbert space is not simply the kron product of each site's Hilbert space.
Extra structure should be considered.
Here we use the Jordan-Wigner transformation to map the fermion operators to spin operators.
Namely,
c_j = \prod_{i<j} \sigma^z_i \sigma^-_j,
c_j^\dagger = \prod_{i<j} \sigma^z_i \sigma^+_j,
where
\sigma^z = |0><0| - |1><1|,
\sigma^- = |0><1|, \sigma^+ = |1><0|.
In this program, we use `F` to represent \sigma^z, `c` to represent \sigma^-, and `c_dag` to represent \sigma^+.
"""
dim = 2**N
H_tol = np.zeros((dim, dim),dtype=np.float64)
# add term `c_dag c`
op_list = [a_dag, a] + [I for _ in range(N-2)]
op_list_boundary = [a] + [F for _ in range(N-2)] + [a_dag]
for i in range(N):
if i == N-1:
# the cross-boundary term `c_{N-1}^\dagger c_0` should be carefully considered.
op_list = op_list_boundary
H_sub = op_list[0]
for j in range(1, N):
H_sub = np.kron(H_sub, op_list[j])
H_tol += H_sub
op_list = [op_list[-1]] + op_list[:-1] # shift to right by one site
H_tol += H_tol.T # add term `c_cdag`
w, v = np.linalg.eigh(H_tol)
print(w)
return w, v
def measure_occupy(N, i, state_vec):
'''
measure the occupation number of site i
'''
op_list = [I for _ in range(N)]
op_list[i] = op_n # 名字取得不好，下次一定
measure_op_mat = op_list[0]
for j in range(1, N):
measure_op_mat = np.kron(measure_op_mat, op_list[j])
occupy_num = np.einsum('i,ij,j->', state_vec, measure_op_mat, state_vec)
return occupy_num
def count_degency(w, thershold):
'''
count the degeneracy of each energy level
'''
w_origin = deepcopy(w)
w_unique = []
w_degency = []
while w_origin.size > 0:
next_w_origin = []
w_unique.append(w_origin[0])
w_degency.append(0)
for item in w_origin:
if(np.isclose(w_origin[0], item, atol=thershold)):
w_degency[-1] += 1
else:
next_w_origin.append(item)
w_origin = np.array(next_w_origin)
print("unique : ", w_unique)
print("degency: ", w_degency)
return w_unique, w_degency
if __name__ == '__main__':
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(color_codes=True)
N = 6
w, v = FreeFermion_ED(N)
# particle measurement
# occupy_distribution = np.zeros(N)
# plt.figure()
# for i in range(2**N):
# for j in range(N):
# occupy_distribution[j] = measure_occupy(N, j, v[:,i])
# plt.subplot(4,4,i+1)
# plt.plot(occupy_distribution, 'o-')
# plt.text(0,2.2,'$%d,\ \langle \sum_i \hat n_i \\rangle = %f$' % (i, np.sum(occupy_distribution)))
# plt.ylim(-0.1,3.5)
# if i%4 == 0:
# plt.ylabel('$\langle \hat n_i \\rangle$')
# if int(i/4) == 3:
# plt.xlabel('site')
# plt.show()
grid = plt.GridSpec(1, 4, wspace=0.5, hspace=0.5)
fig = plt.figure(figsize=(10,5))
plt.subplot(grid[0,0])
sns.rugplot(y=w, height=1)
plt.xticks([])
plt.ylabel("Energy Level (of the whole system)")
w_unique, w_degency = count_degency(w, 1e-10)
for index in range(len(w_unique)):
plt.text(0.09, w_unique[index], str(w_degency[index])+" fold")
k_list = np.array([2*np.pi/N*i for i in range(N)])
plt.subplot(grid[0,1:4]).set_xticks([i/N for i in range(N)])
plt.xlabel('k / 2$\pi$')
plt.ylabel('Single Particle Energy Spectrum')
plt.plot(k_list/(2*np.pi), 2*np.cos(k_list), 'o-')
plt.show()
# plt.subplot(122).set_xticks([i for i in range(N)])
# plt.xlabel('site')
# plt.plot(occupy_distribution, 'o-', label='occupy distribution (for %d),\n$\sum_i n_i = %f$' % (level,np.sum(occupy_distribution)))
# plt.legend()
