半导体器件建模仿真-Scharfetter-Gummel方法

在半导体器件仿真中，需要求解耦合的非线性偏微分方程，即漂移-扩散模型。其核心方程之一是载流子的电流连续性方程。对于电子电流密度公式为：

\[J_n = q \mu_n n E + q D_n \nabla n \]

其中 \(E\) 是电场强度: \(E = -\nabla \psi\)， $ D_n$ 是电子扩散系数，通过爱因斯坦关系 $D_n = \mu_n V_T $ 与迁移率关联，$ V_T = kT/q $ 是热电压。

例如，对于一维稳态连续性方程 \(dJ_n/dx=qR\)（\(R\) 为净复合率），得到：

\[\frac{d}{dx} \left( \mu_n V_T \frac{dn}{dx} - \mu_n n \frac{d\psi}{dx} \right) = R \]

这是一个关于电子浓度 $ n $ 和电势 $ \psi $ 的方程。

在半导体内部，尤其是在空间电荷区（如PN结附近），电势 \(\psi\) 的变化非常剧烈。电子浓度 \(n\) 与电势呈玻尔兹曼关系（在热平衡下）：\(n \propto \exp(\psi / V_T)\)。这意味着即使电势变化几个 \(V_T\)（室温下约 26 mV），浓度就会发生数量级的变化。

如果用中心差分来离散化电流密度项，例如：

\[J_{i+1/2} \approx q\mu_n V_T \frac{n_{i+1} - n_i}{\Delta x} - q\mu_n \frac{\psi_{i+1} - \psi_i}{\Delta x} \cdot \frac{n_i + n_{i+1}}{2} \]

当电势差 \(\Delta\psi = \psi_{i+1} - \psi_i\) 很大（正或负）时，这种离散化会严重违背物理上的指数关系，导致数值解出现非物理振荡（不稳定性）、不收敛，或要求网格极其精细才能得到正确结果。

Scharfetter-Gummel 核心思想

1969年，Scharfetter 和 Gummel 提出了一种精妙的离散化方法，其核心思想是：

在两个相邻的网格点之间的单元（即网格边）上，假设电流密度 ( J ) 和电场 ( E ) 为常数，然后求解关于载流子浓度的一阶常微分方程。

这个假设将问题局部化，并允许我们推导出两个网格点浓度之间的精确关系，该关系自然地包含了指数变化。

考虑相邻网格点 \(i\) 和 \(i+1\) 之间的边中点\(i+1/2\)。在边单元上，假设 $ J_n$ 和 \(E = -d\psi/dx\) 为常数，记为 $ J_{i+1/2} $ 和 $E_{i+1/2} $。那么电流方程变为：

\[J_{i+1/2} = -q \mu_n n E_{i+1/2} + q D_n \frac{dn}{dx} \]

这是一个关于 \(n(x)\) 的一阶线性常微分方程。

将方程改写为：

\[\frac{dn}{dx} - \frac{E_{i+1/2}}{V_T} n = \frac{J_{i+1/2}}{q D_n} \]

这是一个标准的一阶线性方程。其积分因子为 \(\exp\left( \int \frac{E}{V_T} dx \right) = \exp\left( \frac{E_{i+1/2} x}{V_T} \right)\)。注意，$ E_{i+1/2} = -(\psi_{i+1} - \psi_i) / \Delta x = -\Delta\psi / \Delta x $。

在区间 \([x_i, x_{i+1}]\) 上积分，利用边界条件 \(n(x_i) = n_i\)， \(n(x_{i+1}) = n_{i+1}\)，经过运算可以得到 \(J_{i+1/2}\) 的表达式：

\[J_{i+1/2} = q \mu_n V_T \cdot \frac{1}{\Delta x} \left[ B\left( -\frac{\Delta\psi}{V_T} \right) n_i - B\left( \frac{\Delta\psi}{V_T} \right) n_{i+1} \right] \]

其中 $ B(\cdot)$ 是伯努利函数：

\[B(x) = \frac{x}{e^x - 1} \]

这就是著名的 Scharfetter-Gummel 公式。它表达了边上的电流 \(J_{i+1/2}\) 如何通过相邻节点的浓度 \(n_i, n_{i+1}\) 和电势差 \(\Delta\psi\) 来计算。

数值特性：

• 稳定性：无论 \(\Delta\psi\) 多大（正或负），伯努利函数 \(B(x)\) 总是正值且平滑的。当 \(x \to +\infty\) 时， \(B(x) \to 0\)；当 \(x \to -\infty\) 时，\(B(x) \to -x\)；当 \(x \to 0\) 时， \(B(x) \to 1\)。这保证了数值通量不会出现非物理的爆炸或振荡。
• 一致性：当网格很细（\(\Delta\psi \to 0\)）时，SG 格式退化为标准的中心差分格式。

Scharfetter-Gummel 方法 是一种针对对流-扩散型方程（其中对流项系数远大于扩散项，即高场情况）的特殊空间离散格式。其精髓在于在网格边上局部求解常数通量/常数场假设下的精确解，并以此构造数值通量。SG 方法是所有现代半导体器件仿真软件（如 Silvaco TCAD, Synopsys Sentaurus）的基石。它使得在相对较粗的网格上稳定、准确地模拟具有强内建电场和偏压的器件（如 PN 结、MOSFET）成为可能。

示例：PN结

考虑一个简单的PN结，网格点 \(i\) 和 \(i+1\)。假设：

• 热电压 \(V_T = 0.0259 \, \text{V}\) (300K)
• 电子迁移率 \(\mu_n = 1000 \, \text{cm}^2/\text{V·s}\)
• 网格间距 \(\Delta x = 0.1 \, \mu\text{m}\)
• 电势差：\(\psi_{i+1} - \psi_i = 0.26 \, \text{V}\)(这很大，等于 \(10V_T\))
• 边界浓度：\(n_i = 1 \times 10^{15} \, \text{cm}^{-3}\)，\(n_{i+1} = 1 \times 10^{10} \, \text{cm}^{-3}\)

电势差0.26V出现在PN结耗尽区。根据玻尔兹曼关系，电子浓度应满足：

\[\frac{n_{i+1}}{n_i} = \exp\left(-\frac{\Delta\psi}{V_T}\right) = \exp(-10) \approx 4.54 \times 10^{-5} \]

中心差分电流密度离散为：

\[J_{i+1/2}^{\text{center}} = q\mu_n V_T \frac{n_{i+1} - n_i}{\Delta x} - q\mu_n \frac{\psi_{i+1} - \psi_i}{\Delta x} \cdot \frac{n_i + n_{i+1}}{2} \]

代入数值：

\[\frac{V_T}{\Delta x} = \frac{0.0259}{10^{-5}} = 2590 \text{V/cm}\\ \frac{\Delta\psi}{\Delta x} = \frac{0.26}{10^{-5}} = 2.6 \times 10^4 \text{V/cm}\\ q\mu_n = 1.6 \times 10^{-19} \times 1000 = 1.6 \times 10^{-16} \]

扩散项：

\[D_{\text{diff}} = q\mu_n V_T \frac{n_{i+1} - n_i}{\Delta x} = 1.6 \times 10^{-16} \times 2590 \times (10^{10} - 10^{15}) = (4.144 \times 10^{-13}) \times (-9.999 \times 10^{14}) = -4.144 \times 10^{2} \, \text{A/cm}^2 \]

漂移项：

\[D_{\text{drift}} =- q\mu_n \frac{\Delta\psi}{\Delta x} \cdot \frac{n_i + n_{i+1}}{2} = -1.6 \times 10^{-16} \times 2.6 \times 10^4 \times \frac{10^{15} + 10^{10}}{2} = -(4.16 \times 10^{-12}) \times (5.00005 \times 10^{14}) = -2.08 \times 10^{3} \, \text{A/cm}^2 \]

总电流：

\[J_{\text{center}} = -414.4 - 2080 = -2494.4 \, \text{A/cm}^2 \]

SG电流密度离散为：

\[J_{i+1/2}^{\text{SG}} = q\mu_n V_T \cdot \frac{1}{\Delta x} \left[ B\left(-\frac{\Delta\psi}{V_T} \right) n_i - B\left( \frac{\Delta\psi}{V_T} \right) n_{i+1} \right] \]

其中 \(\Delta\psi/V_T = 10\)。

伯努利函数：

\[B(10) = \frac{10}{e^{10} - 1} = \frac{10}{22026 - 1} \approx 4.54 \times 10^{-4} \\ B(-10) = \frac{-10}{e^{-10} - 1} = \frac{-10}{4.54 \times 10^{-5} - 1} \approx \frac{-10}{-0.9999546} \approx 10.000455 \]

代入：

\[J_{\text{SG}} = 1.6 \times 10^{-16} \times 2590 \times [10.000455 \times 10^{-4} \times 10^{15} - 4.54 \times 10^{10}]\approx 1.0\text{A/cm}^2 \]

SG方法结果：这个值非常小 (接近零), 基本符合物理预期, 小误差可能来自数值舍入和有限的网格精度.

总结

中心差分方法在强电场下完全失效：

• 无法满足热平衡条件
• 产生巨大的非物理电流
• 需要极其精细的网格才能得到勉强可用的结果

Scharfetter-Gummel方法具有显著优势：

• 严格保持热平衡条件
• 在粗网格上也能稳定准确
• 物理基础坚实，数学特性优良
• 已成为半导体器件仿真的标准方法

SG方法的成功在于放弃了简单线性近似，转而采用基于物理的指数拟合，从而能够正确处理半导体中载流子的指数分布特性。