# 利用 Python 展示 Tight-Binding 模型(1D 與 2D band structure)
1D band structure
---
### 自由粒子能量公式的物理意義
公式如下：
$$
E(k) = \frac{1}{2}k^2
$$
這是自由粒子（free particle)在量子力學中的能量與波向量（$k$）的關係式。在物理上，這個公式可以從以下幾個角度來理解：
### 1. 它從哪裡來?(公式推導)
原本的動能公式是：
$$
E = \frac{p^2}{2m}
$$
而在量子力學裡，動量 $p$ 跟波向量 $k$ 的關係是：
$$
p = \hbar k
$$
把這個代進去，就變成：
$$
E(k) = \frac{\hbar^2 k^2}{2m}
$$
如果用**自然單位**或**簡化單位制**（讓 $\hbar = 1$, $m = 1$），就變成你看到的這個：
$$
E(k) = \frac{1}{2}k^2
$$
### 2. 它代表什麼？
這個公式說明的是：
> 粒子的**能量大小與它的波向量平方成正比**。
- $k$ 越大，代表粒子「震動得越快」、「波長越短」→ 動量越大
- 所以 $E(k)$ 就越高，對應到比較「快的粒子」
### 3. 在能帶圖裡的意義
在固態物理的能帶結構圖中，這個公式畫出來是個**拋物線**。當我們在晶格中考慮電子的行為時，這就是最簡單的近似模型，用來模擬電子在不受外力情況下的能量分布。
### 4. 另外補充 $k$ 是什麼？
它描述一個波（像電子波、光波、聲波）在空間中的震動頻率和方向。
數學上定義是：
$$
\vec{k} = \frac{2\pi}{\lambda}
$$
其中：
- $\lambda$：波長
- $\vec{k}$：波向量，是一個有方向的量，單位是 $\mathrm{m^{-1}}$（1/長度）
## 接下來簡單介紹程式與物理意義
```
import numpy as np
import matplotlib.pyplot as plt
```
#### 引入numpy 以及 matplotlib
NumPy（Numerical Python）是Python數學與科學計算的基礎工具庫在tight-binding 模型中會遇到的計算(例如:建立向量、對波向量取內積、矩陣哈密頓量運算等)
```
k_array = np.linspace(-4*np.pi, 4*np.pi, 201)
```
建立一個波向量（wave vector）陣列，從 −4𝜋到 4𝜋均勻取 201 點。這代表我們考慮多個週期的 𝑘 空間。
```
E_array = 1/2 * k_array**2
```
這就是自由粒子的能量公式，此圖將會是一個開口向上的對稱拋物線，代表在沒有晶格位能的情況下粒子的能量分佈。
```
plt.xlabel("wave vector ")
plt.ylabel("Energy")
```
設定$x$軸為波向量${k}$，$y$軸為能量$𝐸$
```
plt.xlim([-2.5*np.pi, 2.5*np.pi])
plt.ylim([0,65])
```
侷限範圍，使圖集中在第一布里淵區附近（-2.5$\pi$ ~ 2.5$\pi$），能量上限為 65，讓圖形不會太亂。
```
plt.plot(
k_array, E_array,'black',
k_array+2*np.pi,E_array, 'black',
k_array-2*np.pi,E_array, 'black',
k_array+4*np.pi,E_array, 'black',
k_array-4*np.pi,E_array, 'black',
k_array+6*np.pi,E_array, 'black',
k_array-6*np.pi,E_array, 'black',
k_array+8*np.pi,E_array, 'black',
k_array-8*np.pi,E_array, 'black',
```
上方程式碼將自由粒子能量公式間格2$\pi$做位移，分別做到+-8$\pi$。
```
k_array, np.cos(k_array)+1, 'blue',
k_array ,4.5*np.cos(k_array)+11, 'blue',
k_array, -8*np.cos(k_array)+31, 'blue',
k_array, 11*np.cos(k_array)+60, 'blue',
```
加入了四條藍色曲線，形狀為$cos({k})$模擬能帶的 dispersion（色散關係）。這些曲線的形狀和間距模擬了不同能帶之間的能隙與對稱性。
```
[np.pi, np.pi],[0,65],'gray',
[-np.pi, -np.pi],[0,65],'gray',
```
(x=$\pi$,0)到(x=$\pi$,65)畫一條灰色線，以及(x=-$\pi$,0)到(x=-$\pi$,65)也畫一條灰線，就完整標示出第一布里淵區為k=[-$\pi$ ~ $\pi$]。
> Dispersion 意指能量 $E$ 隨波向量 $k$ 的變化，描述電子在晶格中可擁有的能量範圍與運動特性。白話問就是:$$ 在不同波向量下，電子有哪些能量可選 $$
注意:
本圖中能量 $E(k)$ 是以簡化單位制表示（即 $\hbar = 1$, $m = 1$）由上方推導區得知，因此數值為無單位的規模化結果，用來強調能量與波向量的**關係**，而非實際物理量。
2D band structure
---
### 以石墨烯(Graphene)為例
藍色&紅色皆為碳(Carbon)不同色只是為了方便觀察。
### 接下來模擬典型的 能帶色散關係圖（band dispersion relation）
```
a0 = 1.418 # [a0] = Å，代表碳原子之間的距離
```
定義 graphne 的晶格常數（原子間距，如上圖紅球和藍球的距離），單位是 Å
```
k_points = 301
```
決定${k}$空間中要取幾個點來掃描波向量$k_x$, $k_y$(分別為x及y方向的波向量)值越大，解析度越高，圖就越平滑。但過高可能導致電腦當機。
```
kx_list = np.linspace(-np.pi/a0, np.pi/a0, k_points)
ky_list = kx_list
```
建立一個 $k_x$,$k_y$ 皆分別為[-$\pi$ ~ $\pi$]的範圍取301個點!這樣就涵蓋第一布里淵區（對稱的波向量區域）。
```
energy_list = []
```
建立一個空的energy_list 來存所有能量本徵值。每個 k 點會有兩個能量（上下兩條能帶）
```
for kx in kx_list:
for ky in ky_list:
```
這兩層迴圈讓每個$k_x$,$k_y$的組合都能計算到一次能量!也就是完全掃過剛建立的空間網格。
```
ham_matrix = np.array([
[0,np.exp(1j*np.dot([kx,ky],[0,a0]))+np.exp(1j*np.dot([kx,ky],[a0*np.sqrt(3)/2,-a0/2]))+np.exp(1j*np.dot([kx,ky],[-a0*np.sqrt(3)/2,-a0/2]))],
[np.exp(1j*np.dot([kx,ky],[0,-a0]))+np.exp(1j*np.dot([kx,ky],[-a0*np.sqrt(3)/2,a0/2]))+np.exp(1j*np.dot([kx,ky],[a0*np.sqrt(3)/2,a0/2])),0]
])
```
這裡有點攏長,最主要的是建立一個2×2 的 Hamiltonian（哈密頓量矩陣)，在 tight-binding 模型裡，只有最近鄰（nearest neighbors）的原子之間有非零的跳躍機率。因此
$$
H(\vec{k}) =
\begin{bmatrix}
0 & f(\vec{k}) \\
f^*(\vec{k}) & 0
\end{bmatrix}
$$
- 左上、右下是對角項 → 設為 0，代表 藍-藍、紅-紅 不互相跳遷（模型假設）。
- 右上$f(\vec{k})$意旨從 A → B 的跳遷
- 左下$f^*(\vec{k})$意旨是從 B → A 的跳遷，是右上的複共軛（因為 Hamiltonian 必須是 Hermitian）
### 什麼是 $f(\vec{k})$ ?
$$f(k) = e^{i \vec{k} \cdot \vec{\delta}_1} + e^{i \vec{k} \cdot \vec{\delta}_2} + e^{i \vec{k} \cdot \vec{\delta}_3}$$
這三個 𝛿 是 graphene 中 A → B 的三個最近鄰位移向量。
其中：
```
[0,
np.exp(1j*np.dot([kx, ky], [0, a0])) +
np.exp(1j*np.dot([kx, ky], [a0*np.sqrt(3)/2, -a0/2])) +
np.exp(1j*np.dot([kx, ky], [-a0*np.sqrt(3)/2, -a0/2]))
]
```
這三項就是 石墨烯(Graphene)中 A → B 的三個跳遷方向！
- 第一項:向上跳,[0,$a_0$]
- 第二項:往右下跳一格,$\left[ \frac{a_0 \sqrt{3}}{2}, -\frac{a_0}{2} \right]$
- 第三項:往左下跳一格,$\left[ -\frac{a_0 \sqrt{3}}{2},\ -\frac{a_0}{2} \right]$
#### 也就是:
$$
\vec{\delta}_1 = \left( 0,\ a_0 \right),\quad
\vec{\delta}_2 = \left( \frac{a_0 \sqrt{3}}{2},\ -\frac{a_0}{2} \right),\quad
\vec{\delta}_3 = \left( -\frac{a_0 \sqrt{3}}{2},\ -\frac{a_0}{2} \right)
$$
#### 接著分別套入delta位置
```
np.exp(1j * np.dot([kx, ky], delta))
```
意思是:
- 把波向量$k$ =($k_x$,$k_y$)跟空間位移向量$\vec{\delta}$做內積
- 再套入相位因子$e^{i \vec{k} \cdot \vec{\delta}}$
- 模擬晶格中波動因跳遷所帶來的干涉效應
```
[np.exp(1j*np.dot([kx, ky], [0, -a0])) +
np.exp(1j*np.dot([kx, ky], [-a0*np.sqrt(3)/2, a0/2])) +
np.exp(1j*np.dot([kx, ky], [a0*np.sqrt(3)/2, a0/2])), 0]
```
和上方雷同,但皆差一個負號為(-$\vec{\delta}$)
```
energy_list.append(np.linalg.eigvalsh(ham_matrix))
```
用 numpy 的 eigvalsh() 計算這個 2x2 Hermitian 矩陣的eigenvalue（也就是上下兩條能帶）。然後把它加入能量列表中。
```
kx_list, ky_list = np.meshgrid(kx_list, ky_list)
```
將一維的 $k_x$, $k_y$ 陣列變成二維網格，對應到倒空間中的每一個 k 點座標。
```
energy_list = np.array(energy_list)
```
把能量 list 轉成 numpy 陣列以便後續操作
```
energy_list = np.reshape(energy_list, (k_points, k_points, 2))
```
把能量陣列reshape 成 (301, 301, 2)，對應 (kx, ky) 格點網格，每點有兩條能帶。
```
plot3D = plt.axes(projection = '3d')
```
開啟一個 3D 繪圖座標軸！
```
plot3D.set_xlabel('kx [1/a]')
plot3D.set_ylabel('ky [1/a]')
plot3D.set_zlabel('Energy')
```
定 x, y, z 軸的標籤名稱（$k_x$, $k_y$, $能量$）。
```
plot3D.view_init(10,15)
```
設定 3D 視角（仰角 10 度，方位角 15 度），依照個人需求改成你喜歡的角度!
```
plot3D.plot_surface(kx_list, ky_list, energy_list[:,:,0],
cmap='winter', linewidth = 1, antialiased = True,
rstride = 1, cstride = 1)
```
畫出第一條能帶（最低能態），用藍色 colormap 當配色。
```
plot3D.plot_surface(kx_list, ky_list, energy_list[:,:,1],
cmap='autumn', linewidth = 1, antialiased = True,
rstride = 1, cstride = 1)
```
畫出第二條能帶（較高能態），用紅橘色。這樣就能清楚區分兩條能帶！
然後把它加入能量列表中。
利用 tight-binding 模型模擬 Graphene 的能帶結構，透過掃描倒空間的波向量 ($k_x$, $k_y$)，建立 2x2 Hamiltonian 並計算能量eigenvalue本徵值。
程式模擬的是 Graphene 的能帶結構，利用最簡單的 tight-binding 模型，只考慮最近鄰跳遷，建立出一個 2x2 的 Hamiltonian，並計算每個 $k$ 點的能量本徵值。畫出來的上下兩條能帶展示了 石墨烯(Graphene) 的色散關係，特別是在 K 點會出現特徵性的 Dirac cone 結構。
Sign in with Wallet
Connect another wallet