能带论
金属电子论不考虑电子电子相互作用和电子离子相互作用；考虑这些相互作用，可以得到能带论。
对于晶体中的多体问题，其哈密顿量为：
能带论做出了以下假设：
- 绝热近似：处理电子问题时，认为离子实在格点上固定不动；在上述假设中，电子和离子实被分开，即电声无耦合。
- 单电子近似：用平均场来代替电子电子间和电子离子间的相互作用；
- 周期场近似：平均场是周期性势场。这使得电子的波函数成为布洛赫函数。
这样，能带论的核心问题就是求解以下周期势场的单电子问题：
这和我们在量子力学中求解的周期势场问题是相同的。一种严格解法是平面波法，两种近似解法为近自由电子近似和紧束缚近似。如何选择两种近似方法，得看材料中的电子表现更像波还是粒子，进一步，是在周期势场中穿梭还是被束缚在原子附近。
Bloch定理
周期势场中的哈密顿量是正格矢的周期函数，那么波函数也是吗？并非。考虑周期场中单电子的波函数为$\psi_k^n(\vec{r})$，显然平移后的波函数为：
考虑玻恩-冯卡门边界条件：
或者写为分量的形式：
将波函数拆开，根据Bloch定理，其相位函数可以确定，振幅函数同样也是一个周期函数：
这种波函数$u_k^n(\vec{r})$称为布洛赫波函数。
在描述布洛赫波函数的时候，我们用到了量子数$n$和$k$。量子数$n$表示该电子波函数在晶体原子间距变大的时候，退化为孤立原子的第n个本征态；量子数$k$表示该电子波函数的波矢，不过并不与电子的动量成比例。
能谱结构
- 写出能量本征值方程，并对其作规范变化： 记$\hat H_k=-\frac{\hbar^2}{2m}(\nabla-i\vec{k})^2+V(r)$，其仍然是一个厄米算符。
- 对于确定的$\vec{k}$，由Sturm-Liouville定理，有无数个分立的能量本征值$E_n(\vec{k})$。
- 对于确定的$n$，能量本征值$E_n(\vec{k})$是一个周期函数，即$E_n(\vec{k})=E_n(\vec{k}+\vec{K_h})$，其中$\vec{K_h}$是倒格矢。这说明能带的计算可以在一个布里渊区内完成。
- 显然，由于周期性，确定的$n$对应的能量本征值$E_n(\vec{k})$是一个能带（有上下界）；
- 由于玻恩-冯卡门边界条件，$\vec{k}$的取值是离散的，所以能带在$\vec{k}$轴上是分立的。但是，由于宏观晶体的尺寸较大，$\vec{k}$的取值是准连续的。
平面波法
利用平面波分解可以严格求解周期势场中的单电子薛定谔方程问题。将单电子波函数用倒格矢展开：
其中，平面波展开的系数为：
最终波函数可以写为：
代入波动方程中，并在左侧作用$\langle \vec{k}+\vec{K_{h’}}|$得到：
其中，对于势能的计算来自于：
这正是周期势的傅里叶展开：
一般取$\bar{V}=0$。
至此，我们问题转化为线性齐次方程组的问题。利用行列式为0，可以求解能量本征值的方程：
近自由电子近似
尽管平面波是严格求解周期势场中单电子薛定谔方程的方法，但在高阶行列式中，求解往往是难以收敛的。考虑势场的空间变化十分微弱，我们使用量子力学中的微扰论来处理。这时候电子的行为十分接近自由电子，这就是近自由电子近似：
零级近似
如果做零级近似，即$\Delta V=0$，代入平均场，电子的本征解即为平面波解：
一级近似的非简并微扰
将波函数写为：
依照微扰理论，可以计算出一阶微扰：
代入，即可得到能量本征值的二阶近似解：
一般来说，由于周期势很弱，电子本征波函数和电子本征能量和自由电子相差不远，但可见在某些简并情况下发散。
简并微扰
观察上述能量表达式，可知当$\vec{k}^2=(\vec{k}+\vec{K_h})^2$时，能量本征值是简并的。此时，其他散射波都可以忽略，最后剩下两个方程：
其中耦合的波函数为：
可以通过行列式：
解出：
布里渊区与能隙
上述非简并微扰下，当$|\vec{k}|=|\vec{k}+\vec{K_h}|$时，计算发散。这个条件等同于$\vec{k}$落在倒格矢的垂直平分面上。倒格矢的垂直平分面依次围成布里渊区。
为什么会在布里渊区的边界上发生能隙现象呢？入射波遭布拉格衍射后，一般来说相互抵消，不影响入射波；但在布里渊区边界上，入射波和衍射波叠加形成驻波。以下是一个简单的图像：
当$k=-\frac{\pi n}{a}$时，$k’=k+\frac{2\pi n}{a}=-k$，波函数$|\psi_k^0\rangle=\frac{1}{\sqrt{L}}e^{ikx}$是简并的，此时需要考虑简并微扰。
左边同乘$\langle\psi^0_k|,\langle\psi^0_{k’}|$，有：
其中$\langle\psi^0_{k’}|\Delta V|\psi^0_k\rangle=\langle k’|\Delta V|k\rangle=V_n$，由解存在性，有：
解得：
相应的本征波函数为：
对于非简并微扰，入射波和衍射波不耦合，但在简并微扰中，我们通过对角化提取出了驻波。这意味着能量更高的波函数喜欢出现在两个原子的中间，能量较低的波函数喜欢出现在原子上，这正是能隙的成因。
近似简并微扰
实际问题中，在布里渊区边界的附近，简并微扰也是近似有效的：
我们讨论其中的两种退化情形：
当$|E_k^0-E_{k’}^0|\gg|V_n|$时：
如果再做一次近似，得到：
这说明，当距离布里渊边界较远的时候，能量曲线为抛物线，可以忽略微扰。
当$|E_k^0-E_{k’}^0|\ll|V_n|$时：
代入$k=-\frac{n\pi}{a}(1-\Delta),k’=\frac{n\pi}{a}(1+\Delta)$，有：
其中$T_n=\frac{\hbar^2n^2\pi^2}{2ma^2}$，这意味着在距离布里渊边界较近的时候，$\Delta\rightarrow0$，能量曲线为带有能隙的抛物线。可以看到，能量差值是以布里渊边界为零点的抛物线函数，所以能看到能隙是弯弯的。
紧束缚近似
万尼尔函数
在平面波方法中，我们将布洛赫波函数按平面波（复数）展开。实际上，由于布洛赫波函数是倒空间的周期函数，所以也可以用正空间的周期函数展开：
其中$a_n(\vec{R}_l,\vec{r})$称为万尼尔函数：
可见万尼尔函数只是宗量$\vec{r}-\vec{R}_l$的函数（这说明万尼尔函数是以$\vec{R}_l$为中心的局域函数），记：
不难证明万尼尔函数是正交完备的：
- 正交性：
- 完备性：
紧束缚近似假设
紧束缚近似的基本假设为：晶体中每个原子的势场对电子都有较强的束缚，使得电子的行为接近孤立原子中的电子。这时候，其他原子对该电子的吸引势被看作微扰。
在紧束缚近似假设下，可以选择孤立原子的定域波函数$\varphi_n(\vec{r}-\vec{R}_l)$作为万尼尔函数：
于是周期势场中的单电子波函数可以用上述函数展开：
指标n表示不同轨道类型的波函数，如s和p轨道。该方法之所以是一种近似，是因为不同孤立原子的波函数之间是有重叠的，但在紧束缚假设中，重叠的部分很少，可以近似正交：
代入周期势场的单电子薛定谔方程：
其中$V(\vec{r})=\sum_l U(\vec{r}-\vec{R}_l)$。左乘$\varphi_n^*(\vec{r}-\vec{R}_{l})$并对$\vec{r}$积分，利用正交性得到：
只取最邻近项：
讨论以上公式的物理意义。当原子的间距$\vec{R}_l$较大的时候，电子表现出孤立原子的特征：
即能带是由简并态的直线进化而来的。这点从下述讨论的三维简单立方晶体的能带结构中可以看出。
考虑晶体结构，简并态被打破，这是由于其他原子对电子的吸引势。交叠积分$J(\vec{R}_l-\vec{R}_{l’})$其实代表了第l’个原子跃迁到第l个原子时的跃迁矩阵元：
可以用微扰的角度看待$J(0)$项。
如图所示，$V(\vec{r})-U(\vec{r}-\vec{R}_{l})$是一个在$\vec{r}=\vec{R}_{l}$处的微扰势场，所以当$l=l’$时，$J(0)$可以通过微扰论得到：
正交平面波方法
从上述的讨论中，我们发现，离子实区域外的电子受到周期势场的作用，更像一个波；而离子实区域内的电子受到孤立原子势场的作用，更像一个粒子。对于两种电子，分类讨论可以解决平面波方法的收敛问题。
若内层有M个电子态，他们仍可以使用相应的独立原子轨道为拟合：
对于外层的非定域电子（离子实内振荡，离子实外平滑），如果使用平面波，就无法避免高能平面波的引入。为了解决这一问题，考虑将此前的M个原子轨道作为基函数的一部分，线性组合得到：
这样可以有效减少平面波的引入。考虑外层波函数与内层波函数的正交性：
得到新的基组叫做正交平面波（Orthogonal Plane Wave，OPW）基组：
这样，外层波函数就可以表示为：
新的正交平面波基组自身反而不正交了，这导致在计算能量本征值时需要引入重叠矩阵的概念。
考虑基组$\{|\phi_i\rangle\}$，薛定谔方程指出存在本征能量和波函数：
其中$|\psi\rangle=\sum_i c_i|\phi_i\rangle$，代入后有：
左乘$\langle \phi_j|$，有：
- 对于正交基组，$\langle \phi_j|\phi_i\rangle=\delta_{ij}$，有：
- 对于非正交基组，$\langle \phi_j|\phi_i\rangle=S_{ij}$，有：这说明在非正交基组下，需要求解的是带有重叠矩阵的广义本征值问题。
赝势方法
将正交平面波展开代入薛定谔方程：
其中：
变为矩阵方程：
行列式为0：
其具有和平面波方法矩阵方程同样的形式。以下是对比：
- 平面波方法：
- 正交平面波方法：
将正交平面波方法导出的方程进行变形：
相比于真实的薛定谔方程，我们将正交平面波基组替换为平面波基组，将真实势能$V(\vec{r})$替换为赝势$U(\vec{r})$。这代表我们通过改变势能为一个更加平滑的赝势来简化问题，得到了平滑解。对于只关心离子实外能带的情况，不失为一种简化。
能带电子密度
能带论除了考虑能带的色散关系，还得考虑能带上的电子密度。考虑能带的准连续性，可以用delta函数的筛选性表示为：
其中$N_n(E)$为能带$n$上能量为$E$的电子数密度。
一般可以表示为更为实用的形式：
其中$dS_E$为能带$n$上能量为$E$的$\vec{k}$空间的面积元，该积分在等能面上积分，通常需要数值计算的帮助。
考虑不同能带的交叠，最终的电子密度为：
自由电子的能态密度
自由电子的能谱为：
能带电子密度为：
同样的道理，可以计算其他维度的能态密度，这里仅给出正比关系：
具体计算过程见 近独立粒子的最概然分布 。这里复述结果（未计入自旋简并度）：
维度 态密度 $D(\epsilon)$ 1 $\dfrac{L}{2\pi\hbar}\sqrt{\dfrac{2m}{\epsilon}}$ 2 $\dfrac{S}{2\pi} \dfrac{m}{\hbar^2}$ 3 $\dfrac{V}{4\pi^2}(\dfrac{2m}{\hbar^2})^\frac32\sqrt{\epsilon}$
紧束缚近似三维简单立方晶体的能带电子密度
考虑紧束缚近似下的简单立方晶体，能带色散关系为：
沿高对称路径，通过数值计算可以得到能带色散关系和能带电子密度的图像如下：
它是由简并态的直线进化而来的：
绘图代码为：import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
import matplotlib.font_manager as fm
# 1. 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei'] # 设置中文字体为黑体
plt.rcParams['axes.unicode_minus'] = False # 解决负号显示异常问题
# 参数设置
J1 = 1.0 # hopping积分 (单位: eV)
a = 1.0 # 晶格常数
E0 = 0.0 # 能带中心
Nk = 300 # 每个方向的k点数量
Nbins = 300 # 能量分箱数
# 生成k空间网格 (第一布里渊区: -π/a ≤ k ≤ π/a)
k = np.linspace(-np.pi/a, np.pi/a, Nk)
kx, ky, kz = np.meshgrid(k, k, k, indexing='ij')
# 计算能带能量 E(k)
def energy(kx, ky, kz, E0, J1, a):
return E0 - 2*J1*(np.cos(kx*a) + np.cos(ky*a) + np.cos(kz*a))
E = energy(kx, ky, kz, E0, J1, a)
# 计算态密度 N(E)
E_min, E_max = E0 - 6*J1, E0 + 6*J1 # 能量范围
E_bins = np.linspace(E_min, E_max, Nbins)
N = np.zeros_like(E_bins)
# 数值积分: 统计每个能量区间的态密度贡献
for i in range(Nbins-1):
mask = (E >= E_bins[i]) & (E < E_bins[i+1])
N[i]= np.sum(mask) # 统计每个能量区间的k点数量
# 归一化 (考虑k空间体积和网格密度)
N = (N-N.min()) / (N.max()-N.min())
# ========================================
# 创建包含两个子图的图形
# ========================================
plt.figure(figsize=(12, 5), dpi=100)
# ========================================
# 子图1: 能带结构图 (左)
# ========================================
plt.subplot(1, 2, 1)
# 定义高对称路径 Γ→X→M→Γ→R
k_path = np.array([
[0, 0, 0], # Γ
[np.pi/a, 0, 0], # X
[np.pi/a, np.pi/a, 0], # M
[0, 0, 0], # Γ
[np.pi/a, np.pi/a, np.pi/a] # R
])
# 在路径上采样
n_points = 100
k_interp = []
for i in range(len(k_path)-1):
for t in np.linspace(0, 1, n_points):
k_interp.append(k_path[i] + t*(k_path[i+1] - k_path[i]))
k_interp = np.array(k_interp)
# 计算能带
E_band = energy(k_interp[:,0], k_interp[:,1], k_interp[:,2], E0, J1, a)
# 绘制能带
x_axis = np.arange(len(k_interp))
plt.plot(x_axis, E_band, 'r-', lw=2)
# 标记高对称点
high_sym_pos = [0, n_points, 2*n_points, 3*n_points, 4*n_points-1]
high_sym_labels = ['Γ', 'X', 'M', 'Γ', 'R']
for pos, label in zip(high_sym_pos, high_sym_labels):
plt.axvline(x_axis[pos], color='gray', ls=':', alpha=0.5)
# 设置x轴刻度
plt.xticks([x_axis[pos] for pos in high_sym_pos], high_sym_labels)
plt.xlim(0, len(x_axis)-1)
# 坐标轴标签和标题
plt.xlabel('高对称点', fontsize=12)
plt.ylabel(r'能量 $E-J_0-nJ_1$ (eV)', fontsize=12)
plt.title('(a) 能带结构 (Γ→X→M→Γ→R路径)', y=1.02, fontsize=12)
plt.grid(True, alpha=0.3)
# ========================================
# 子图2: 能带电子态密度曲线 (右)
# ========================================
plt.subplot(1, 2, 2)
# 绘制态密度曲线
plt.plot(E_bins[:-1], N[:-1], 'k-', lw=2, label='态密度')
# 在Γ点附近进行√E拟合 (能带底部附近)
# 找到能带最小值附近的区域 (最低10%的能量范围)
E_range = E_max - E_min
fit_range = 0.1 * E_range # 取最低10%能量范围进行拟合
mask_fit = (E_bins < (E_min + fit_range)) & (E_bins > E_min) # 拟合区域
# 准备拟合数据
x_fit = E_bins[mask_fit][:-1] - E_min # 相对于Γ点的能量
y_fit = N[mask_fit][:-1] # 对应的态密度
# 定义√E拟合函数
def sqrt_func(x, a):
return a * np.sqrt(x)
# 进行拟合 (使用非负最小二乘)
from scipy.optimize import curve_fit
popt, pcov = curve_fit(sqrt_func, x_fit, y_fit, p0=[1.0])
# 生成拟合曲线
E_range = E0 - E_min
x_plot = np.linspace(0, E_range, 100)
y_plot = sqrt_func(x_plot, *popt)
# 绘制拟合曲线
plt.plot(x_plot + E_min, y_plot, 'b--', lw=2,
label=f'√E拟合')
# 标记Γ点位置
plt.axvline(E_min, color='gray', ls=':', alpha=0.5)
high_sym_pos = [0, Nbins//3, 2*Nbins//3-1, 3*Nbins//3-2]
high_sym_labels = ['Γ(n=-6)', 'X(n=-2)', 'M(n=2)','R(n=6)']
x_axis_dos = E_bins[:-1]
for pos, label in zip(high_sym_pos, high_sym_labels):
plt.axvline(x_axis_dos[pos], color='gray', ls=':', alpha=0.5)
# 设置x轴刻度
plt.xticks([x_axis_dos[pos] for pos in high_sym_pos], high_sym_labels)
# 坐标轴标签和标题
plt.xlabel(r'能量 $E-J_0-nJ_1$ (eV)', fontsize=12)
plt.ylabel(r'电子态密度 $N(E)$', fontsize=12)
plt.title('(b) 能带电子态密度曲线', y=1.02, fontsize=12)
plt.ylim(0, None)
plt.grid(True, alpha=0.3)
plt.legend(loc='upper left')
# 调整子图间距
plt.tight_layout()
plt.show()
这里用到了能带色散关系的求导：
对于带底和带顶进行特殊讨论：
- 带底：当$k_x=k_y=k_z=0$时：其中$m^*_-=\frac{\hbar^2}{2a^2J(1)}$，这说明带底的有效质量为正。
- 带顶：当$k_x=k_y=k_z=\frac{\pi}{a}$时：其中$m^*_+=-\frac{\hbar^2}{2a^2J(1)}$，这说明带顶的有效质量为负。
带顶和带底的色散关系都类似自由电子，所以在上述图像中，的确可以看到$\sqrt{E}$形式的态密度曲线。可以理论计算出态密度公式：
紧束缚近似下石墨烯的能带电子密度
先利用二次量子化后的紧束缚模型计算色散关系：
其中$t$为跃迁积分，$\langle i,j\rangle$表示最近邻原子对。利用二次量子化的傅里叶展开和正交性公式：
代入后得到：
则哈密顿简化为：
将其写成矩阵形式：
代入三个原子的坐标：
求解矩阵的特征值：
这里的$a$是晶格常数。当
import numpy as np |
