最近准备认真学一下半导体器件仿真，学习前，准备先自己仿真一下最简单的一维 PN 结器件，来熟悉一下半导体器件背后的物理原理，以及最重要的 TCAD 软件背后的仿真原理，打好器件仿真的基础。

问题概述

首先，这是一个需要通过数学求解的物理问题。因此，我们首先需要明确，这个问题的输入（也就是问题的已知量）有哪些，以及问题的输出（即需要求解的物理量）有哪些。

输入

一些器件参数和材料参数，主要包括掺杂浓度、器件尺寸、材料本征参数（比如介电常数等）。

输出

主要包括载流子浓度（电子浓度和空穴浓度）和电势。为什么是这三个？这是因为通过这三个物理量，我们就可以求出其他的物理量，比如电场。

注：为什么不求电场，而是电势呢？这是因为电场是矢量，求解它不仅需要求大小，更需要考虑方向，这会是问题十分麻烦，因此通过引入电势，来帮助我们计算问题。因此，电势的引入，是充当一个中间量的作用，来帮助我们求出我们所关注的电场。

物理方程

明确了相关物理量之后，我们就需要列出这个问题相关的微分方程。这是一个电学问题，因此，我们从 Maxwell 方程组入手。Maxwell 方程如下：

$$
\begin{aligned}
&\nabla \times \mathbf{E}=-\frac{\partial \mathbf{B}}{\partial t} \\
&\nabla \times \mathbf{H}=\mathbf{J}+\frac{\partial \mathbf{D}}{\partial t} \\
&\nabla \cdot \mathbf{D}=\rho \\
&\nabla \cdot \mathbf{B}=0
\end{aligned}
$$

泊松方程

接下来我们对其进行处理，首先我们假设介质是线性、各向同性的，可以得到如下关系：

$$
\mathbf{D}=\varepsilon\mathbf{E},\quad \mathbf{B}=\mu\mathbf{H}
$$

同时，利用电场与电势的关系，将电场方程转化为电势方程：

$$
\mathbf{E}=-\nabla \varphi
$$

首先，我们对 $\nabla \cdot \mathbf{D}=\rho$ 进行处理：

$$
\nabla \cdot \mathbf{D}=\nabla \cdot(\varepsilon \mathbf{E})
$$

假设介质是均匀的，即介电常数不随位置变化，可以得到：

$$
\nabla \cdot \mathbf{D}=\nabla \cdot(\varepsilon \mathbf{E})=\varepsilon \nabla\cdot (-\nabla\varphi)=-\varepsilon \nabla^{2}\varphi=\rho
$$

可以得到：

$$
\nabla^{2} \varphi=-\frac{\rho}{\varepsilon}
$$

这个方程称为泊松方程。可以发现，这个方程的得出是有一定的条件的，即介质是线性、各向同性、均匀的。

电流连续性方程

接下来我们对 $\nabla \times \mathbf{H}=\mathbf{J}+\frac{\partial \mathbf{D}}{\partial t}$ 进行处理。等式两边同时去散度，可以得到左边为 0，即：

$$
0=\nabla\cdot (\nabla \times \mathbf{H})=\nabla\cdot \mathbf{J}+ \frac{\partial}{\partial t}(\nabla \cdot \mathbf{D})=\nabla\cdot \mathbf{J}+ \frac{\partial \rho}{\partial t}
$$

所以可以得到：

$$
\nabla\cdot \mathbf{J}+ \frac{\partial \rho}{\partial t}=0
$$

这个方程被称为电流连续性方程。可以发现这个方程的物理意义还是比较明显的。第一项表示单位体积中流出的电流，表示单位时间内流出的电荷，第二项表示体电荷密度的变化量。即单位体积内流出的电流等于单位时间内该体积中净电荷密度的减少。

半导体器件中的电流连续性方程

在半导体中，共有两种载流子，分别为电子和空穴。所以电流密度也是由两部分组成，分别是电子电流 $\mathbf{J_n}$ 和空穴电流 $\mathbf{J_p}$。所以电流连续性方程可以写成：

$$
\nabla\cdot \mathbf{J_n}+\nabla\cdot \mathbf{J_p}-q\frac{\partial n}{\partial t}+q\frac{\partial p}{\partial t}=0
$$

所以

$$
\nabla\cdot \mathbf{J_p}+q\frac{\partial p}{\partial t}=-\nabla\cdot \mathbf{J_n}+q\frac{\partial n}{\partial t}
$$

这个式子的物理意义为：单位时间、单位体积内流出以及增加的正电荷数等于单位时间、单位体积内流入以及减少的负电荷数。而电子和空穴的变化，便是由于半导体内的产生和复合中心。我们用 $G$ 表示产生率，用 $R$ 表示复合率。那么上式可以写成：

$$
\nabla\cdot \mathbf{J_p}+q\frac{\partial p}{\partial t}=-\nabla\cdot \mathbf{J_n}+q\frac{\partial n}{\partial t}=q(G-R)
$$

化简可得：

$$
\frac{\partial p}{\partial t}=-\frac{1}{q}\nabla \cdot \mathbf{J_p}+G-R
$$

$$
\frac{\partial n}{\partial t}=\frac{1}{q}\nabla \cdot \mathbf{J_n}+G-R
$$

这两个式子便分别是空穴连续性方程和电子连续性方程。

半导体器件中的物理方程

至此，我们就从 Maxwell 方程出发，得到了半导体器件背后的物理方程，如下：

$$
\begin{aligned}
\nabla^{2} \varphi&=-\frac{\rho}{\varepsilon} \\
\frac{\partial p}{\partial t}&=-\frac{1}{q}\nabla \cdot \mathbf{J_{p}}+G-R\\
\frac{\partial n}{\partial t}&=\frac{1}{q}\nabla \cdot \mathbf{J_n}+G-R
\end{aligned}
$$

可以发现，上式中的未知量有很多，显然是求不出来解的。因此，我们希望将上面的 $\rho$、$\mathbf{J_n}$、$\mathbf{J_p}$、$G$、$R$ 转换成 $n$ 、$p$、$\varphi$。这就涉及到许多的物理模型了。我们需要知道的是，物理模型并不是唯一的，且准确度不相同，有些模型考虑的因素很多，有些考虑的比较少。因此，在仿真时，我们需要根据实际需求，选择不同的物理模型。但是，无论物理模型有多么复杂，我们始终需要明确它的对象，并且物理模型再复杂，上面的三个方程是不会改变的。不同的物理模型，也就是对“将上面的 $\rho$、$\mathbf{J_n}$、$\mathbf{J_p}$、$G$、$R$ 转换成 $n$ 、$p$、$\varphi$”这个过程，进行不同的转换。因此，搞懂上面的物理原理，可以极大的帮助我们理清 TCAD 仿真中的流程和各种设置。

引入物理模型

$\rho$

首先我们来看体电荷密度，这涉及到体电荷模型，在本次仿真中，我们假设杂质完全电离。此时，可以得到：

$$
\rho=q(N_{d}-N_{a}+p-n)
$$

其中，$N_d$ 为施主杂质浓度，$N_a$ 为受主杂质浓度，$p$ 为空穴浓度，$n$ 为电子浓度。

$G，R$

在这次的问题中，我们不考虑产生率，只考虑复合。即 $G=0$。在 TCAD 有许多与产生有关的模型，比如碰撞电离，还可以指定产生率等参数。重点说一下复合模型。在 TCAD 中，有许多复合模型，这是因为半导体器件中存在许多复合机制，比如辐射复合、SRH 复合、俄歇复合、表面复合等等。在这里，我们只考虑 SRH 复合。

在半导体物理中，学过 SRH 复合的复合率的表达式，即：

$$
R=\frac{np-n_{i}^{2}}{\tau_{p}(n+n_i)+\tau_{n}(p+p_{i})}
$$

其中，这个式子涉及到本征载流子浓度 $n_i$。在半导体中存在禁带变窄效应等负效应，因此，在仿真中，可以根据实际问题，利用相关的物理模型。

$\mathbf{J_{n}},\mathbf{J_{p}}$

有关电流密度，依旧有许多物理模型，主要包括漂移扩散模型、热力学模型、流体力学模型、密度梯度模型，具体就不介绍了。在这里，我们使用最简单的漂移扩散模型。根据半导体物理的知识，可以写出：

$$
\begin{aligned}
\mathbf{J_{n}}&=q\mu_{n}n\mathbf{E}+qD_{n}\nabla n \\
\mathbf{J_{p}}&=q\mu_{p}p\mathbf{E}-qD_{p}\nabla p
\end{aligned}
$$

其中，$\mu$ 为迁移率，$D$ 为扩散系数，它们满足爱因斯坦关系式：

$$
D_{n}=\frac{kT}{q}\mu_{n},\quad D_{p}=\frac{kT}{q}\mu_{p}
$$

有关迁移率，由于半导体中存在许多的散射机构，这也涉及到许多的迁移率模型。在这里，为简单起见，我们就把迁移率当成一个常数。

总结

综上，我们可以得到微分方程组：

$$
\begin{aligned}
&\nabla^{2} \varphi=-\frac{q(N_{d}-N_{a}+p-n)}{\varepsilon} \\
&\frac{\partial p}{\partial t}=-\frac{1}{q}\nabla \cdot \mathbf{J_{p}}-R\\
&\frac{\partial n}{\partial t}=\frac{1}{q}\nabla \cdot \mathbf{J_n}-R\\
&R=\frac{np-n_{i}^{2}}{\tau_{p}(n+n_i)+\tau_{n}(p+p_{i})}\\
&\mathbf{J_{n}}=-q\mu_{n}n\nabla \varphi+qD_{n}\nabla n \\
&\mathbf{J_{p}}=-q\mu_{p}p\nabla \varphi-qD_{p}\nabla p
\end{aligned}
$$

可以发现，上面的微分方程组的未知量为：$\varphi,p,n,R,\mathbf{J_{n}},\mathbf{J_{p}}$，共有六个未知量，而同时也有六个方程，因此是可以求解的。

一维 PN 结对应的微分方程组

建立如上图所示的坐标系。我们假设准静态情况，即 $\frac{\partial p}{\partial t}=\frac{\partial n}{\partial t}=0$，令 $N(x)=N_d(x)-N_a(x)$ 上面的微分方程组可以化为如下形式：

$$
\begin{aligned}
&\frac{\mathrm{d}^{2} \varphi(x)}{\mathrm{d} x^{2}} =-\frac{q(N(x)+p(x)-n(x))}{\varepsilon} \\
&\frac{1}{q}\frac{\mathrm{d} J_{p}(x)}{\mathrm{d} x} +R=0\\
&\frac{1}{q}\frac{\mathrm{d} J_{n}(x)}{\mathrm{d} x} -R=0\\
&R=\frac{n(x)p(x)-n_{i}^{2}}{\tau_{p}(n(x)+n_i)+\tau_{n}(p(x)+p_{i})}\\
&J_{n}(x)=-q\mu_{n}n(x)\frac{\mathrm{d} \varphi(x)}{\mathrm{d} x} +qD_{n}\frac{\mathrm{d} n(x)}{\mathrm{d} x} \\
&J_{p}(x)=-q\mu_{p}p(x)\frac{\mathrm{d} \varphi(x)}{\mathrm{d} x} -qD_{p}\frac{\mathrm{d} p(x)}{\mathrm{d} x}
\end{aligned}
$$

边界条件

求解微分方程，需要知道边界条件。在 TCAD 中，也有许多的物理模型，对应不同的边界条件，比如欧姆接触、肖特基接触等。在这里，我们考虑欧姆接触。假设 PN 结两端加了一个电压 $V_A$。在边界处，需要满足电中性条件和热平衡条件，即：

$$
\begin{aligned}
&N(0)+p(0)-n(0)=0 \\
&N(L)+p(L)-n(L)=0
\end{aligned}
$$

$$
\begin{aligned}
&n(0)p(0)=n_{i}^{2} \\
&n(L)p(L)=n_{i}^{2}
\end{aligned}
$$

电势边界条件为：

$$
\begin{aligned}
&\varphi(0)=\frac{kT}{q}ln(\frac{n(0)}{n_i})\\
&\varphi(L)=V_{A}-\frac{kT}{q}ln(\frac{p(L)}{n_i})
\end{aligned}
$$

归一化

接下来，理论上，就可以通过计算机求解微分方程了。但是我们可以注意到，对于载流子浓度等这些物理量，往往都是十几个数量级，而对于器件尺寸等物理量，往往都很小。它们之间的数据跨度是很大的。对于计算机而言，是以比特为单位存储数据的，如果数据跨度很大的话，很容易出现溢出现象，因此需要进行归一化操作。归一化操作是将物理量除以一个适当的比例系数，缩小数值范围，这个比例系数称为归一化因子。归一化因子的选取具有一定的任意性，只要保证模型能够成立就行。

引入本征得拜长度：

$$
L_{D}=\sqrt{\frac{\varepsilon kT}{q^{2}n_{i}}}
$$

可以得到常见的归一化系数表：

通过归一化操作，我们可以得到归一化后的微分方程组：

$$
\begin{aligned}
&\frac{\mathrm{d}^{2} \varphi}{\mathrm{d} x^{2}} =-(N+p-n) \\
&\frac{\mathrm{d} J_{p}}{\mathrm{d} x} =-R\\
&\frac{\mathrm{d} J_{n}}{\mathrm{d} x} =R\\
&R=\frac{np-1}{\tau_{p}(n+1)+\tau_{n}(p+1)}\\
&J_{n}=-\mu_{n}n\frac{\mathrm{d} \varphi}{\mathrm{d} x} +\mu_{n}\frac{\mathrm{d} n}{\mathrm{d} x} \\
&J_{p}=-\mu_{p}p\frac{\mathrm{d} \varphi}{\mathrm{d} x} -\mu_{p}\frac{\mathrm{d} p}{\mathrm{d} x}
\end{aligned}
$$

边界条件为：

$$
\begin{aligned}
& N(0)+p(0)-n(0)=0 \\
& N(L)+p(L)-n(L)=0 \\
& n(0)p(0)=1 \\
& n(L)p(L)=1 \\
& \varphi(0)=ln(n(0))\\
& \varphi(L)=V_{A}-ln(p(L))
\end{aligned}
$$

这种微分方程一般是没有解析解的，因此求解需要用到数值求解方法，求解方法主要有差分法和有限元法。在之后的文章中，会简单介绍一下利用差分法，如何求解这个方程。