本問題集來自身邊計(jì)算勢力腦的自發(fā)涌現(xiàn)。此問題搞凝聚態(tài)物理計(jì)算的容易遇到或想到，本科初學(xué)凝聚態(tài)導(dǎo)論的過程中也會(huì)去思考，不是第一次發(fā)現(xiàn)此問題了，也不是第一次思考這一問題了。總之，我想這是一個(gè)新手普遍會(huì)遇到的問題，而非一次鉆牛角尖式的發(fā)問。作者沒有系統(tǒng)學(xué)習(xí)過量子多體理論，本文的一些說法一定是存在某些問題的，請辨證閱讀！

類似問題：https://www.zhihu.com/question/450788784

文題中“初步回答”的含義是，僅僅在本文把問題表述出來，作者沒有做到清晰表述，所以程度上只能歸類到初步。

一、問題集描述（把問題清晰表述出來，往往就是答案本身）

請看以下有關(guān)聯(lián)的模糊表述：

1. 規(guī)定玻戈留玻夫變換后的新算符和原算符的玻色費(fèi)米型一致（線性組合并不能保證此條件，比如費(fèi)米算符的線性組合不一定依然滿足費(fèi)米反對易關(guān)系，需要額外規(guī)定）；

2. 可以直接對角化費(fèi)米子體系（在粒子數(shù)表象下）獲取能譜；

3. “不能”直接對角化玻色子體系（在粒子數(shù)表象下）獲取能譜；

請移步https://jsci.cnu.edu.cn/HTMLfile/2021/1/202101003.html

4. 產(chǎn)生湮滅算符的玻色對易關(guān)系及費(fèi)米反對易關(guān)系與波函數(shù)交換對稱及交換反對稱的關(guān)聯(lián)；

4i. *無關(guān)內(nèi)容：波函數(shù)交換反對稱與泡利不相容原理的等同；

5. 任意N個(gè)費(fèi)米子體系的波函數(shù)（交換反對稱），可用Slatter行列式構(gòu)筑（實(shí)際上是Dirac先提出），行列式維度N；

請移步https://blog.sciencenet.cn/blog-100379-869543.html

6. 對于N個(gè)格點(diǎn)，M個(gè)費(fèi)米子(M<N)占據(jù)這N個(gè)格點(diǎn)，在粒子數(shù)表象下有C_N^M個(gè)基底，基底可寫為類似Slatter行列式的形式，也可以用費(fèi)米子產(chǎn)生算符作用在|0,0,...,0>上表達(dá)，再利用費(fèi)米子的產(chǎn)生湮滅算符的反對易關(guān)系，可以構(gòu)造出交換反對稱的性質(zhì)?？梢杂眠@樣表述的基底夾住哈密頓量，寫出費(fèi)米子體系哈密頓量的所有矩陣元，之后可直接對角化；

7.?能夠構(gòu)造一個(gè)粒子數(shù)不守恒的基底，由二項(xiàng)式定理，知2^N=C_N^0+C_N^1+...+C_N^N，也可以由排列組合問題的基本方法，知每一個(gè)格點(diǎn)只有兩種可能（0態(tài)或1態(tài)），故N個(gè)格點(diǎn)共2^N個(gè)態(tài)（無相互作用費(fèi)米子體系的直積態(tài)表示）；

8.?無相互作用玻色子天然適配“單粒子基底”的直積態(tài)（線性空間，對易），而費(fèi)米子體系則需要加入額外的結(jié)構(gòu)（暫且自行命名為行列式直積態(tài)），來構(gòu)造出交換反對稱的特性（線性空間+某種規(guī)則或稱結(jié)構(gòu)，反對易），比如從粒子數(shù)表象01到自旋表象↑↓的約當(dāng)-魏格納變換就可以實(shí)現(xiàn)這一額外引入的結(jié)構(gòu)；

請移步https://www.zhihu.com/question/450788784

9. ……

二、用Jordan-Wigner變換對角化費(fèi)米子體系哈密頓量

以下圖片內(nèi)容及附錄的程序來自緊扣的dagger的通訊作者，C_{jl}^\ddagger，他是真正的double dagger，也是問題的提供人，沒有他，本篇專欄將不會(huì)被寫成。

他計(jì)算了一個(gè)4格點(diǎn)的僅考慮最近鄰hoping的體系，體系粒子數(shù)不守恒體系（每個(gè)格點(diǎn)至多激發(fā)1個(gè)費(fèi)米子，故總粒子數(shù)至多4個(gè)）。

結(jié)果如下：

當(dāng)然也可繼續(xù)計(jì)算6格點(diǎn)的能譜，會(huì)看到tight-binding的余弦曲線特征，隨格點(diǎn)增加逐漸顯露：

附錄：C_{jl}^\ddagger的程序源碼

經(jīng)double dagger同意，源碼可以貼出。源碼沒有經(jīng)過任何優(yōu)化，且某些關(guān)鍵行被注釋掉了，需要自行恢復(fù)。

'''

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()

? ??

? ??