区域 III（x > a）：V = 0

通过在 x=0 和 x=a 处的边界条件（波函数及其导数的连续性），可以解出各系数的关系。透射系数定义为透射波与入射波的概率流之比：

几个重要的物理结论：

1. T ≠ 0
2. ：即使 E < V0，透射系数也不为零！
3. 指数衰减
4. ：当势垒较厚（κa >>0）时，

1. 势垒越厚、越高，透射越难
2. ：T随 a 和 V0 指数下降

定态分析给出了透射概率，但看不到"过程"。要真正理解隧穿，我们需要模拟一个波包撞击势垒的动态演化。这里采用波包而不用平面波是因为平面波在空间上无限延展无法表示一个从左边飞来的“粒子”

含时薛定谔方程

初始条件：一个从左边飞向势垒的高斯波包：

• x0：波包初始位置
• σ：波包宽度
• k0：中心波矢，决定波包的平均动量

直接用显式差分求解含时薛定谔方程会遇到稳定性问题。Crank-Nicolson格式是更好的选择：

含时薛定谔方程的时间演化算符是

直接计算指数很难，我们用有限差分近似。

显式格式（前向Euler）：

问题：不稳定，概率不守恒。

Crank-Nicolson格式取时间中点的平均：

整理成：

注意到 A = B*，所以 |Aψ|=|Bψ|这保证了概率守恒（幺正性）。

哈密顿量

离散化后变成三对角矩阵：

所以 A 和 B 也是三对角的矩阵。

Python代码实现

下面是完整的动态隧穿模拟代码：

import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.animation import FuncAnimationfrom scipy.linalg import solve_banded# ============ 物理参数（采用原子单位 ℏ = m = 1）============# 这样薛定谔方程变成：i ∂Ψ/∂t = -0.5 ∂²Ψ/∂x² + V·Ψ# 空间网格N = 1000 # 网格点数x_min, x_max = -50, 50 # 空间范围x = np.linspace(x_min, x_max, N)dx = x[1] - x[0]# 时间参数dt = 0.5 # 时间步长N_steps = 600 # 总步数# 势垒参数V0 = 0.5 # 势垒高度barrier_width = 3.0 # 势垒宽度barrier_start = -barrier_width / 2barrier_end = barrier_width / 2# 构建势能V = np.zeros(N)V[(x >= barrier_start) & (x <= barrier_end)] = V0# 初始波包参数x0 = -20.0 # 初始位置sigma = 3.0 # 波包宽度k0 = 0.7 # 中心波矢（对应动能 E = k0²/2 = 0.245 < V0）# 平均能量E_avg = 0.5 * k0**2print(f"势垒高度 V0 = {V0:.3f}")print(f"波包平均能量 E = {E_avg:.3f}")print(f"E/V0 = {E_avg/V0:.2%}")print(f"经典预测：{'可以越过' if E_avg > V0 else '无法越过'}")print(f"量子预测：有透射概率！")# ============ 初始化波函数 ============def gaussian_wavepacket(x, x0, sigma, k0): """归一化高斯波包""" norm = (1 / (2 * np.pi * sigma**2))**0.25 envelope = np.exp(-(x - x0)**2 / (4 * sigma**2)) phase = np.exp(1j * k0 * x) return norm * envelope * phasepsi = gaussian_wavepacket(x, x0, sigma, k0)# 验证归一化norm_check = np.trapz(np.abs(psi)**2, x)print(f"初始归一化检验：∫|Ψ|² dx = {norm_check:.6f}")# ============ Crank-Nicolson 矩阵构建 ============# (1 + i·dt·H/2)·Ψ^{n+1} = (1 - i·dt·H/2)·Ψ^n# 其中 H = -0.5·d²/dx² + V# 系数alpha = 1j * dt / (4 * dx**2)beta = 1j * dt / 2# 构建三对角矩阵（banded storage format for scipy.linalg.solve_banded）# 矩阵 A = 1 + i·dt·H/2 的对角元和非对角元main_diag_A = 1 + 2*alpha + beta * Vupper_diag_A = -alpha * np.ones(N-1)lower_diag_A = -alpha * np.ones(N-1)# 矩阵 B = 1 - i·dt·H/2main_diag_B = 1 - 2*alpha - beta * Vupper_diag_B = alpha * np.ones(N-1)lower_diag_B = alpha * np.ones(N-1)# 将 A 存储为 banded 格式AB = np.zeros((3, N), dtype=complex)AB[0, 1:] = upper_diag_A # 上对角线AB[1, :] = main_diag_A # 主对角线AB[2, :-1] = lower_diag_A # 下对角线def apply_B(psi): """计算 B·Ψ""" result = main_diag_B * psi result[:-1] += upper_diag_B * psi[1:] result[1:] += lower_diag_B * psi[:-1] return resultdef crank_nicolson_step(psi): """单步 Crank-Nicolson 演化""" # 右端向量 b = apply_B(psi) # 求解 A·Ψ^{n+1} = b psi_new = solve_banded((1, 1), AB, b) return psi_new# ============ 时间演化并记录 ============print("\n开始时间演化...")psi_history = [psi.copy()]prob_density_history = [np.abs(psi)**2]# 记录透射和反射概率x_barrier_center = 0T_history = [] # 透射概率（势垒右侧）R_history = [] # 反射概率（势垒左侧）save_interval = 3 # 每隔几步保存一帧for step in range(N_steps): psi = crank_nicolson_step(psi) # 计算透射和反射概率 prob_density = np.abs(psi)**2 T = np.trapz(prob_density[x > barrier_end], x[x > barrier_end]) R = np.trapz(prob_density[x < barrier_start], x[x < barrier_start]) T_history.append(T) R_history.append(R) if step % save_interval == 0: psi_history.append(psi.copy()) prob_density_history.append(prob_density)print(f"演化完成，共 {len(psi_history)} 帧")# 最终结果T_final = T_history[-1]R_final = R_history[-1]total = T_final + R_finalprint(f"\n最终结果：")print(f" 透射概率 T = {T_final:.4f} ({T_final*100:.2f}%)")print(f" 反射概率 R = {R_final:.4f} ({R_final*100:.2f}%)")print(f" 总概率 T+R = {total:.4f}")# ============ 解析透射系数对比 ============def analytical_transmission(E, V0, a): """解析透射系数（单能量）""" if E >= V0: return 1.0 # 简化处理 kappa = np.sqrt(2 * (V0 - E)) sinh_term = np.sinh(kappa * a)**2 T = 1 / (1 + V0**2 * sinh_term / (4 * E * (V0 - E))) return TT_analytical = analytical_transmission(E_avg, V0, barrier_width)print(f" 解析透射系数（中心能量）T_analytical = {T_analytical:.4f}")# ============ 动画可视化 ============fig, axes = plt.subplots(2, 1, figsize=(12, 9), gridspec_kw={'height_ratios': [3, 1]})# 主图：波函数演化ax1 = axes[0]ax1.set_xlim(x_min, x_max)ax1.set_ylim(0, 0.25)ax1.set_xlabel('Position x', fontsize=12)ax1.set_ylabel('|Ψ(x)|²', fontsize=12)ax1.set_title('Quantum Tunneling Through a Barrier', fontsize=14)# 画势垒barrier_height_display = 0.2ax1.fill_between([barrier_start, barrier_end], 0, barrier_height_display, alpha=0.3, color='gray', label='Barrier')ax1.axvline(barrier_start, color='gray', linestyle='--', alpha=0.5)ax1.axvline(barrier_end, color='gray', linestyle='--', alpha=0.5)# 初始化线条line_psi, = ax1.plot([], [], 'b-', linewidth=2, label='|Ψ|²')line_real, = ax1.plot([], [], 'g-', linewidth=1, alpha=0.5, label='Re(Ψ)×0.3')ax1.legend(loc='upper right')ax1.grid(True, alpha=0.3)# 文本标注text_time = ax1.text(0.02, 0.95, '', transform=ax1.transAxes, fontsize=11)text_prob = ax1.text(0.02, 0.88, '', transform=ax1.transAxes, fontsize=11)# 副图：透射/反射概率随时间变化ax2 = axes[1]ax2.set_xlim(0, N_steps * dt)ax2.set_ylim(0, 1.1)ax2.set_xlabel('Time', fontsize=12)ax2.set_ylabel('Probability', fontsize=12)time_array = np.arange(N_steps) * dtax2.plot(time_array, T_history, 'b-', linewidth=1.5, label='Transmitted')ax2.plot(time_array, R_history, 'r-', linewidth=1.5, label='Reflected')ax2.axhline(T_analytical, color='b', linestyle='--', alpha=0.5, label=f'Analytical T={T_analytical:.3f}')ax2.legend(loc='right')ax2.grid(True, alpha=0.3)# 时间指示线time_indicator, = ax2.plot([], [], 'k-', linewidth=2)def init(): line_psi.set_data([], []) line_real.set_data([], []) time_indicator.set_data([], []) return line_psi, line_real, time_indicator, text_time, text_probdef animate(frame): psi_frame = psi_history[frame] prob = prob_density_history[frame] line_psi.set_data(x, prob) line_real.set_data(x, 0.3 * np.real(psi_frame) + 0.1) # 偏移显示 t = frame * save_interval * dt text_time.set_text(f't = {t:.1f}') # 当前透射/反射概率 if frame * save_interval < len(T_history): T_curr = T_history[frame * save_interval] R_curr = R_history[frame * save_interval] text_prob.set_text(f'T = {T_curr:.3f}, R = {R_curr:.3f}') # 时间指示线 time_indicator.set_data([t, t], [0, 1.1]) return line_psi, line_real, time_indicator, text_time, text_probani = FuncAnimation(fig, animate, init_func=init, frames=len(psi_history), interval=30, blit=True)plt.tight_layout()# 保存动画print("\n保存动画中...")ani.save('quantum_tunneling.gif', writer='pillow', fps=30)print("动画已保存为 quantum_tunneling.gif")plt.show()# ============ 静态对比图：不同能量的隧穿 ============fig2, axes2 = plt.subplots(2, 2, figsize=(14, 10))# 不同波矢（能量）的对比k_values = [0.5, 0.7, 0.9, 1.1]for idx, k0_test in enumerate(k_values): ax = axes2[idx // 2, idx % 2] E_test = 0.5 * k0_test**2 # 重新模拟 psi_test = gaussian_wavepacket(x, x0, sigma, k0_test) # 存储几个关键时刻 snapshots = [] snapshot_times = [0, 150, 300, 450] psi_temp = psi_test.copy() for step in range(max(snapshot_times) + 1): if step in snapshot_times: snapshots.append((step * dt, np.abs(psi_temp)**2)) psi_temp = crank_nicolson_step(psi_temp) # 画势垒 ax.fill_between([barrier_start, barrier_end], 0, 0.2, alpha=0.2, color='gray') # 画不同时刻的波函数 colors = plt.cm.viridis(np.linspace(0.2, 0.9, len(snapshots))) for (t, prob), color in zip(snapshots, colors): ax.plot(x, prob, color=color, linewidth=1.5, label=f't={t:.0f}') ax.set_xlim(-40, 40) ax.set_ylim(0, 0.2) ax.set_xlabel('Position x') ax.set_ylabel('|Ψ|²') status = "E > V₀ (经典允许)" if E_test > V0 else "E < V₀ (经典禁止)" ax.set_title(f'k₀={k0_test}, E={E_test:.3f}, V₀={V0}\n{status}') ax.legend(loc='upper right', fontsize=8) ax.grid(True, alpha=0.3)plt.tight_layout()plt.savefig('tunneling_comparison.png', dpi=150)plt.show()# ============ 透射系数 vs 能量曲线 ============fig3, ax3 = plt.subplots(figsize=(10, 6))E_range = np.linspace(0.01, 1.0, 500)T_range = [analytical_transmission(E, V0, barrier_width) for E in E_range]ax3.plot(E_range, T_range, 'b-', linewidth=2)ax3.axvline(V0, color='r', linestyle='--', label=f'V₀ = {V0}')ax3.axhline(1, color='gray', linestyle=':', alpha=0.5)# 标记波包平均能量ax3.axvline(E_avg, color='g', linestyle='--', alpha=0.7)ax3.scatter([E_avg], [analytical_transmission(E_avg, V0, barrier_width)], color='g', s=100, zorder=5, label=f'Wave packet E = {E_avg:.3f}')ax3.set_xlabel('Energy E', fontsize=12)ax3.set_ylabel('Transmission Coefficient T', fontsize=12)ax3.set_title('Transmission Through Square Barrier', fontsize=14)ax3.set_xlim(0, 1)ax3.set_ylim(0, 1.1)ax3.legend()ax3.grid(True, alpha=0.3)# 添加注释ax3.annotate('Classical\nforbidden', xy=(0.25, 0.1), fontsize=11, ha='center')ax3.annotate('Classical\nallowed', xy=(0.75, 0.8), fontsize=11, ha='center')plt.tight_layout()plt.savefig('transmission_vs_energy.png', dpi=150)plt.show()

运行代码生成的动画展示了波包隧穿的完整过程：

观察图像易得：

• 透射波包的形状与入射波包相似，但振幅小得多
• 透射概率（蓝线）在波包穿过后趋于稳定值
• 反射 + 透射 = 1，概率守恒

隧穿现象的根源是量子力学的波动性。在势垒内部，波函数不是变成零，而是变成指数衰减的形式。只要势垒不是无限厚，波函数在另一侧就不会完全消失。

这和光学中的全内反射类似：光从玻璃射向空气，入射角大于临界角时会全反射。但如果在空气侧很近的地方放另一块玻璃，光可以"隧穿"过去。

从解析公式可以看出影响透射的因素有

1. 势垒宽度 a
2. ：T∝exp(-2κa)，宽度翻倍，透射指数衰减
3. 势垒高度 V0
4. ：κ∝V0-E，高度增加，衰减加快
5. 粒子能量 E
6. ：能量越接近 V0，κ 越小，透射越容易
7. 粒子质量 m
8. ：κ正∝m，质量越大，隧穿越难

这解释了为什么电子容易隧穿而质子很难——质量差1836倍！下图展示了不同能量下的隧穿情况。

量子隧穿效应的一个重要应用就是扫描隧道显微镜。

扫描隧道显微镜（STM）用一根尖锐的金属针尖靠近样品表面，针尖和表面之间有真空势垒。当距离只有几个埃时，电子可以隧穿，产生可测量的隧穿电流。隧穿电流对距离极其敏感：电流 I ∝ exp(-2κd)。距离变化1埃，电流可以变化一个数量级。这使得STM能够分辨单个原子！

如果这个动画让你对量子隧穿有了新的理解，欢迎分享给同样在量子世界里"穿墙"的朋友