半导体器件建模仿真-DD-Poisson耦合_Boltzmann统计

本篇文章主要介绍基于Boltzmann统计的漂移扩散与Poisson方程耦合的半导体器件建模仿真物理过程。

核心物理公式

对于半经典半导体器件建模仿真，核心物理公式包括: 载流子浓度分布, Poisson方程, 连续性方程, 漂移扩散方程. 下面将分别介绍，

Boltzmann 统计方程

在 Boltzmann 统计规律下, 载流子浓度公式为:

\[n(r) = n_ie^{\frac{E_{Fn}(r)-E_{i}(r)}{kT}}=n_ie^{\frac{\phi_{i}(r)-\phi_{n}(r)}{kT/q}}\\ p(r) = n_ie^{\frac{E_{i}(r)-E_{Fp}(r)}{kT}}=n_ie^{\frac{\phi_{p}(r)-\phi_{i}(r)}{kT/q}} \]

其中 \(n_i\) 是本征载流子浓度, \(E_{Fn}(r)\) 和 \(E_{Fp}(r)\) 分别对应电子和空穴的准费米能级，\(\phi_{n}(r)\) 和 \(\phi_{p}(r)\) 分别对应电子和空穴的准费米势。

Poisson 方程

描述电势与空间电荷分布的关系:

\[\nabla\cdot(\epsilon\nabla\phi) = -\rho=-q(p-n+N_D^{+}-N_A^{-}) \]

其中 \(\epsilon\) 是介电常数, \(N_D^+, N_A^-\) 分别是电离施主浓度和受主浓度.

连续性方程

描述载流子的产生、复合、流入和流出.

• 电子连续性方程：

\[\frac{\partial n}{\partial t}=\frac{1}{q}\nabla\cdot J_n + (G-R) \]

• 空穴连续性方程:

\[\frac{\partial p}{\partial t}=-\frac{1}{q}\nabla\cdot J_p + (G-R) \]

其中 \(G, R\) 分别是产生率和复合率, \(J_n, J_p\) 是电子和空穴电流密度.

漂移扩散方程

\[J_n = q\mu_nnE+qD_n\nabla n\\ J_p = q\mu_pnE-qD_p\nabla p \]

其中 \(\mu_n, \mu_p\) 分别是电子和空穴的迁移率, \(D_n, D_p\) 分别是电子和空穴的扩散系数. 由爱因斯坦关系与迁移率关联: \(D_n = kT/q \mu_n, D_p = kT/q \mu_p\).

将连续性方程转换为对准费米势的方程

由于电场 \(E=-\nabla\phi\), 所以电流密度方程也可写成:

\[J_n = -q\mu_nn\nabla\phi+qD_n\nabla n\\ J_p = -q\mu_pn\nabla\phi-qD_p\nabla p \]

对于 Boltzmann 统计分布，

\[n = n_ie^{\frac{q(\phi-\phi_{n})}{kT}}\\ p = n_ie^{\frac{q(\phi_{p}-\phi)}{kT}} \]

对 \(n, p\) 取梯度, 得

\[\nabla n =\frac{qn}{kT}(\nabla\phi-\nabla\phi_n)\\ \nabla p =\frac{qn}{kT}(\nabla\phi_p-\nabla\phi) \]

化简可得非常简洁的形式:

\[J_n = -q\mu_nn\nabla\phi_n\\ J_p = -q\mu_pp\nabla\phi_p \]

对于稳态仿真, 电子和空穴的连续性方程为:

\[\nabla\cdot(q\mu_nn\nabla\phi_n) = q(G-R)\\ \nabla\cdot(q\mu_pp\nabla\phi_p) = -q(G-R) \]

---

耦合求解流程

未知变量:
\(\phi\): 静电势, \(\phi_n\): 电子准费米势, \(\phi_p\): 空穴准费米势.

依赖变量:
\(n\): 电子浓度, \(p\): 空穴浓度.

控制方程:
1, Poisson 方程: \(\nabla\cdot(\epsilon\nabla\phi) = -q(p-n+N_D^{+}-N_A^{-})\).
2, 电子连续性方程: \(\nabla\cdot(\mu_nn\nabla\phi_n) = G-R\).
3, 空穴连续性方程: \(\nabla\cdot(\mu_pp\nabla\phi_p) = -(G-R)\).

下面介绍两种迭代流程：

一. Gummel 迭代

这是一种解耦迭代方法, 按顺序求解单个方程, 每次使用其他变量的最新值. 它稳健性好, 对初始猜测要求较低, 但收敛速度可能较慢 (尤其在高偏压时).

流程如下:
1, 初始化 (\(k\)=0)
设定初始猜测值, 一般从热平衡状态开始:

• \(\phi^0 = V_{app}\) 在接触处施加偏压，内部通过求解平衡Poisson方程获得.
• \(\phi_n^0 = \phi_p^0 =0\), 平衡时准费米势能处处为零.
• \(n^0,p^0\) 由上面三个变量通过 Boltzmann方程求解.

2, 外循环迭代 (\(k\) --> \(k+1\))

• 步骤 \(A\): 求解Poisson方程, 更新 \(\phi\)
  固定 \(n^k,p^k\)， 求解关于 \(\phi\) 的方程:

\[\nabla \cdot (\epsilon \nabla \phi^{k+1}) = -q(p^k - n^k + N_D^+ - N_A^-) \]

• 步骤 \(B\): 更新载流子浓度
  使用新的 \(\phi^{k+1}\) 和旧准费米势 \(\phi_n^{k},\phi_p^{k}\) 更新载流子浓度:

\[\tilde{n} = n_i \exp\left(\frac{q(\phi^{k+1} - \phi_n^k)}{k_B T}\right)\\ \tilde{p} = n_i \exp\left(\frac{q(\phi_p^k - \phi^{k+1})}{k_B T}\right) \]

• 步骤 \(C\): 求解连续性方程, 更新 \(\phi_n, \phi_p\).
  电子连续性方程: 固定 \(\phi^{k+1}\) 和 \(\tilde{p}\) (用于计算净复合率\(R\) )，将方程重写为以 \(\phi_n\) 为变量的形式：

\[\nabla \cdot \left[ \mu_n n_i e^{q(\phi^{k+1} - \phi_n)/(k_B T)} \nabla \phi_n \right] = -R(\phi^{k+1}, \phi_n, \phi_p^k) \]

这是一个关于 \(\phi_n\) 的非线性方程, 通常用一次牛顿迭代线性化后求解, 求解得到 \(\phi_n^{k+1}\).
空穴连续性方程，也类似地:

\[\nabla \cdot \left[ \mu_p n_i e^{q(\phi_p - \phi^{k+1})/(k_B T)} \nabla \phi_p \right] = R(\phi^{k+1}, \phi_n^{k+1}, \phi_p) \]

求解得到 \(\phi_p^{k+1}\).

• 步骤 \(D\): 最终更新载流子浓度
  使用完全更新后的势 \(\phi^{k+1}, \phi_n^{k+1}, \phi_p^{k+1}\) 计算最终载流子浓度：

\[n^{k+1} = n_i \exp\left(\frac{q(\phi^{k+1} - \phi_n^{k+1})}{k_B T}\right) \\ p^{k+1} = n_i \exp\left(\frac{q(\phi_p^{k+1} - \phi^{k+1})}{k_B T}\right) \]

备注: 这里步骤 \(A\) 和步骤 \(C\) 都需要内循环迭代, 该内循环迭代可采用 Newton 法.

3, 收敛性检查:
计算两次迭代间所有节点上变量的最大变化量或残差范数:

\[\max( |\phi^{k+1}-\phi^k|, |\phi_n^{k+1}-\phi_n^k|, |\phi_p^{k+1}-\phi_p^k| ) < \text{Tolerance} \]

若满足收敛条件，则退出迭代，进入后处理（计算电流等）。若不满足，则 \(k = k+1\)，返回步骤 A 继续迭代。有时会采用阻尼因子 (如 \(\phi_{new} = \omega*\phi_{new} + (1-\omega)*\phi_{old}\)) 来稳定迭代过程。

二. Newton 迭代

Newton 迭代属于全耦合迭代方法，同时求解所有变量的修正量。它收敛速度快（二阶收敛），但对初始猜测要求高，且每次迭代需要构建和求解更大的矩阵，对内存要求大。

流程如下:
1, 初始化 (\(m\)=0): 同Gummel法.
2, 牛顿迭代循环 (\(m\) --> \(m+1\)):

• 步骤 \(A\): 构建耦合系统的残差向量 \(F(U)\).
  • 将三个未知数在每个离散网格点上的值堆叠成一个大向量 \(U = [\phi, \phi_n, \phi_p]^T\)。
  • 将三个控制方程离散化后，其残差也构成一个向量 \(F(U) = [F_{poisson}, F_{continuity-n}, F_{continuity-p}]^T\)。
  • 理想解满足 \(F(U) = 0\)。

\[\nabla \cdot (\epsilon \nabla \phi) + q(p - n + N_D^{+} - N_A^{-}) = 0\\ \nabla \cdot (\mu_n n \nabla \phi_n) - (G - R) =0\\ \nabla \cdot (\mu_p p \nabla \phi_p) +(G - R)=0 \]

• 步骤 \(B\)： 构建雅可比（Jacobian）矩阵 \(J\)
  \(J\) 是 \(F\) 对 \(U\) 的导数矩阵，\(J_{ij} = \partial F_i /\partial U_j\)。它是一个大型稀疏矩阵，结构如下：

\[J = \begin{bmatrix} \frac{\partial F_P}{\partial \phi} & \frac{\partial F_P}{\partial \phi_n} & \frac{\partial F_P}{\partial \phi_p} \\ \frac{\partial F_n}{\partial \phi} & \frac{\partial F_n}{\partial \phi_n} & \frac{\partial F_n}{\partial \phi_p} \\ \frac{\partial F_p}{\partial \phi} & \frac{\partial F_p}{\partial \phi_n} & \frac{\partial F_p}{\partial \phi_p} \end{bmatrix} \]

其中非零块体现了方程间的耦合。例如，\(\partial F_P/\partial\phi_n\) 项源于 Poisson 方程中的 \(-n(\phi, \phi_n)\)。

• 步骤 \(C\)： 求解线性系统，得到更新量
  求解方程：

\[J(U^m) \cdot \delta U = -F(U^m) \]

这是一个大型稀疏线性系统，常用直接法（如LU分解，适用于中小规模问题）或迭代法（如GMRES，适用于大规模问题）求解。得到更新量 \(\delta U = [\delta\phi, \delta\phi_n, \delta\phi_p]^T\)。

• 步骤 \(D\)：更新变量

\[U^{m+1} = U^m + \lambda \delta U \]

• 步骤 \(E\)：更新依赖变量 \(n, p\)
  用新的 \(U^{m+1}\) 通过 Boltzmann 关系更新 \(n^{m+1}, p^{m+1}\)。

3, 收敛性检查：
检查残差范数 \(||F(U^{m+1})||\) 或更新量范数 \(||\delta U||\) 是否小于容差。若收敛则停止，否则继续下一轮牛顿迭代。

Newton法具体实现

其核心在于构建描述所有变量在所有离散点之间耦合关系的全雅可比矩阵 \(J\). 设离散后有 \(N\) 个网格点, 则未知向量 \(U\) 的维度为 \(3N\):

\[U = [\phi_1, ..., \phi_N, \phi_{n,1}, ..., \phi_{n,N}, \phi_{p,1}, ..., \phi_{p,N}]^T \]

残差向量 \(F\) 与之对应：

\[F = [F_{P,1}, ..., F_{P,N}, F_{n,1}, ..., F_{n,N}, F_{p,1}, ..., F_{p,N}]^T \]

雅可比矩阵 \(J\) 是一个 \(3N \times 3N\) 的大型稀疏块矩阵。

下面将展开每个物理块 \(\partial F_X/\partial U_Y\) 的具体含义, 考虑一个格点 \(i\) 及其邻点 \(j\).

1. Poisson 方程残差 \(F_{P,i}\) 的导数

离散后的 \(F_{P,i}\) 来源于 \(\nabla\cdot(\epsilon\nabla\phi) + q(p - n + N_D - N_A) = 0\).

• (1) 对 \(\phi_j\) 的导数 (\(\partial F_{P,i} / \partial\phi_j\)):
  当 \(j = i\) (对角元) 包含两部分:
  • 扩散项贡献： \(\partial /\partial\phi_i [\nabla\cdot(\epsilon\nabla\phi)]_i\). 这部分是纯数学的离散拉普拉斯算子对自身节点的导数，通常是一个与节点周围网格几何相关的常数。
  • 电荷密度贡献： \(q (\partial p_i/\partial\phi_i - \partial n_i/\partial\phi_i)\)。利用 Boltzmann 关系：

\[\frac{\partial n_i}{\partial \phi_i} = \frac{q}{kT} n_i, \quad \frac{\partial p_i}{\partial \phi_i} = -\frac{q}{kT} p_i \]

所以此项为 \(-q^2/(kT) (n_i + p_i)\).

当 \(j \neq i\): 主要来自 \([\nabla\cdot(\epsilon\nabla\phi)]_i\) 项，即离散格式中连接节点 \(i\) 和 \(j\) 的系数.

• (2) 对 \(\phi_{n,j}\) 的导数 (\(\partial F_{P,i} / \partial\phi_{n,j}\)):
  仅当 \(j = i\) 时有值，因为 \(n_i\) 依赖于 \(\phi_{n,i}\).

\[\frac{\partial F_{P,i}}{\partial \phi_{n,i}} = -q \frac{\partial n_i}{\partial \phi_{n,i}} = -q \left( -\frac{q}{k T} n_i \right) = \frac{q^2}{k T} n_i \]

• (3) 对 \(\phi_{p,j}\) 的导数 (\(\partial F_{P,i} / \partial\phi_{p,j}\)):
  同理，仅当 \(j = i\) 时有值，因为 \(p_i\) 依赖于 \(\phi_{p,i}\).

\[\frac{\partial F_{P,i}}{\partial \phi_{p,i}} = q \frac{\partial p_i}{\partial \phi_{p,i}} = q \left( \frac{q}{k T} p_i \right) = \frac{q^2}{k T} p_i \]

2. 电子连续性方程残差 \(F_{n,i}\) 的导数

离散后的 \(F_{n,i}\) 来源于 \(\nabla\cdot(\mu_n n \nabla\phi_n) + R = 0\)。

• (1) 对 \(\phi_j\) 的导数 (\(\partial F_{n,i} / \partial\phi_j\)):
  当 \(j = i\) (对角元) 包含两部分:
  • 漂移-扩散电流项: \(n = n_i exp(q(\phi - \phi_n)/kT)\)，所以 \(\mu_n n\nabla\phi_n\) 通过 \(n\) 依赖于 \(\phi\)。该项的导数比较复杂，包含 \(\partial n/\partial\phi\) 与 \(\nabla\phi_n\) 的乘积项。
  • 复合率项 \(R\)： \(R\) 通常是 \(n\) 和 \(p\) 的函数，而 \(n, p\) 都依赖于 \(\phi\)。例如，对于SRH复合：

\[R \propto \frac{np - n_i^2}{n+p+...} \Rightarrow \frac{\partial R}{\partial \phi} = \frac{\partial R}{\partial n}\frac{\partial n}{\partial \phi} + \frac{\partial R}{\partial p}\frac{\partial p}{\partial \phi} \]

它体现了电势变化如何直接影响电子电流的平衡.

• (2) 对 \(\phi_{n,j}\) 的导数 (\(\partial F_{n,i} / \partial\phi_{n,j}\)):
  当 \(j = i\) (对角元)
  • 电流项的主贡献: 来自 \(\nabla\cdot(\mu_n n \nabla\phi_n)\) 对 \(\phi_{n,i}\) 的导数。离散后，这部分会给出一个类似“离散扩散算子”的结构，但其系数与 \(n\) 相关。
  • 载流子浓度贡献： 同样来自 \(\partial n/\partial\phi_n\) 对电流项和复合项的影响。
  • 复合率贡献： \(\partial R/\partial\phi_n = (\partial R/\partial n)(\partial n/\partial\phi_n)\).

当 \(j \neq i\), 主要来自离散化后 \(\nabla\cdot(\mu_n n \nabla\phi_n)\) 项中连接 \(i\) 和 \(j\) 的系数，这些系数本身也依赖于 \(n_i\) 和 \(n_j\)（即依赖于 \(\phi_i, \phi_{n,i}, \phi_j, \phi_{n,j}\)），因此计算非常繁复。这是Newton法矩阵构建中最复杂的部分之一。

• (3) 对 \(\phi_{p,j}\) 的导数 (\(\partial F_{n,i} / \partial\phi_{p,j}\)):
  主要通过复合率项 \(R\) 耦合。通常仅当 \(j = i\) 时有显著值：

\[\frac{\partial F_{n,i}}{\partial \phi_{p,i}} = \frac{\partial R}{\partial \phi_{p,i}} = \frac{\partial R}{\partial p_i} \frac{\partial p_i}{\partial \phi_{p,i}} = \frac{\partial R}{\partial p_i} \cdot \left( \frac{q}{k_B T} p_i \right) \]

3. 空穴连续性方程残差 \(F_{p,i}\) 的导数

其结构与电子方程完全对称，只需交换 \(n ↔ p\) 和 \(\phi_n ↔ \phi_p\) 的角色，并注意复合率 \(R\) 的符号。

两种迭代对比

Gummel迭代采取“解耦-顺序求解”策略，逻辑清晰，内存消耗相对小，更稳健。而Newton迭代则是“全耦合一次性求解”策略，收敛速度快，是商用TCAD软件在大多数情况下的默认选择。但每次迭代计算量大（需构建和求解大的Jacobian矩阵），且对初始值敏感。

在实际仿真中，常采用混合策略：开始时用几轮Gummel迭代让解进入牛顿法的收敛区域，然后切换到牛顿迭代进行快速精确求解。整个流程的自动化（包括网格生成、离散化、矩阵求解、收敛判断）由半导体器件仿真软件（如Sentaurus TCAD, Silvaco Atlas, COMSOL等）完成。