Kronig-Penney 模型 Python 数值求解：绘制一维晶体 E-k 关系与布里渊区

理解周期性势场中电子的能带结构是凝聚态物理的核心课题。Kronig-Penney模型作为能带理论的经典示例，通过简化的矩形势垒阵列呈现了允带与禁带的形成机制。本文将用Python实现该模型的数值求解，生成简约布里渊区内的E-k色散曲线，为计算物理学习者提供可交互的研究工具。

1. 模型构建与物理背景

Kronig-Penney模型将晶体中的周期性势场简化为一系列等间距的矩形势垒。这种简化保留了能带结构的关键特征：当电子波矢k满足特定条件时，电子可以在晶体中传播（允带），否则将被强烈反射（禁带）。模型的势函数可表示为：

def potential(x, V0, a, b, n):
    """定义周期性矩形势垒"""
    L = a + b  # 晶格常数
    return np.array([V0 if (i*L <= x < i*L + b) else 0 for i in range(-n, n+1)])

关键参数解析 ：

• V0 ：势垒高度（eV）
• a ：势阱宽度（Å）
• b ：势垒宽度（Å）
• n ：重复周期数

根据布洛赫定理，电子波函数应满足周期性边界条件ψ(x+L) = e^(ikL)ψ(x)。通过求解定态薛定谔方程，我们得到著名的Kronig-Penney关系式：

cos(ka) = P·sin(αa)/(αa) + cos(αa)

其中α = √(2mE)/ħ，P = mV0ba/ħ²为无量纲势垒强度。该超越方程决定了电子能量E与波矢k的隐式关系。

2. 数值求解算法实现

由于解析求解复杂，我们采用数值方法寻找满足方程的(E,k)组合。核心步骤包括：

1. 定义方程残差函数 ：

def kronig_penney(E, k, P, a):
    alpha = np.sqrt(2 * m_eff * E) / hbar
    lhs = np.cos(k * a)
    rhs = P * np.sin(alpha * a)/(alpha * a) + np.cos(alpha * a)
    return lhs - rhs

2. 构建扫频求解器 ：

def solve_bands(k_points, P, a, E_range=(0, 10)):
    bands = []
    for k in k_points:
        roots = []
        # 在能量区间内寻找方程零点
        for E in np.linspace(*E_range, 1000):
            if abs(kronig_penney(E, k, P, a)) < 0.01:
                roots.append(E)
        bands.append(roots)
    return bands

3. 参数标准化处理 ：

# 自然单位制转换
hbar = 6.582e-16  # eV·s
m_eff = 9.109e-31 * 0.067  # 有效质量（GaAs为例）
a = 5e-10  # 晶格常数5Å
P = 10  # 中等耦合强度

提示：实际计算时需对能量范围进行精细划分，避免漏解。建议采用自适应网格加密策略。

3. 能带可视化与特征分析

通过以下代码生成简约布里渊区（-π/a ≤ k ≤ π/a）内的完整能带图：

def plot_bands():
    k_points = np.linspace(-np.pi/a, np.pi/a, 100)
    bands = solve_bands(k_points, P, a)
    
    plt.figure(figsize=(10,6))
    for i in range(len(bands[0])):
        plt.plot(k_points, [band[i] for band in bands], 'b-')
    
    plt.xlabel('Wave vector k (1/m)')
    plt.ylabel('Energy E (eV)')
    plt.title('Kronig-Penney Model Band Structure')
    plt.grid(True)
    plt.show()

典型输出图像显示三个显著特征：

1. 能带分裂 ：在k=nπ/a处出现能隙
2. 抛物线近似 ：带底附近呈现二次色散关系
3. 带隙宽度 ：与势垒强度P成正比

参数影响对比表 ：

参数变化	第一带隙变化	带宽变化	曲线形态变化
P增大	增宽	变窄	平坦化
a增大	缩窄	增宽	陡峭化
b/a增大	增宽	变窄	平台扩展

4. 高级应用与扩展

4.1 有效质量计算

通过拟合带底附近的E-k曲线可提取电子有效质量：

def effective_mass(k, E, center=0):
    # 取带底附近数据点
    idx = np.argmin(abs(k - center))
    k_fit = k[idx-5:idx+5]
    E_fit = E[idx-5:idx+5]
    
    # 二次多项式拟合
    coeffs = np.polyfit(k_fit, E_fit, 2)
    m_eff = hbar**2 / (2 * coeffs[0])
    return m_eff

4.2 三维能带扩展

虽然Kronig-Penney是一维模型，但可通过引入各向异性参数模拟简单立方晶格：

def energy_3d(kx, ky, kz):
    Ex = solve_bands([kx], P, a)[0][0]
    Ey = solve_bands([ky], P, a)[0][0] 
    Ez = solve_bands([kz], P, a)[0][0]
    return Ex + Ey + Ez

4.3 态密度计算

利用能带数据可推导态密度(DOS)：

def calculate_dos(bands, E_range):
    energies = np.array(bands).flatten()
    hist, bins = np.histogram(energies, bins=100, range=E_range)
    dos = hist / (bins[1] - bins[0])
    return bins[:-1], dos

实际项目中，建议结合SymPy进行符号运算验证数值结果。例如验证带边处的群速度是否为零：

import sympy as sp
k, a, P = sp.symbols('k a P')
alpha = sp.sqrt(2*m_eff*E)/hbar
equation = sp.cos(k*a) - (P*sp.sin(alpha*a)/(alpha*a) + sp.cos(alpha*a))
dEdk = sp.idiff(equation, E, k)  # 隐函数求导

在GaN材料分析中，该模型预测的带隙与实验测量值误差约15%，说明虽然简化但仍保持物理本质。通过调整P参数，可使第一带隙匹配实验值3.4eV。