COMSOL Blog

Intro to Density-Gradient Theory for Semiconductor Device Simulation

by Chien Liu

November 27, 2019

As semiconductor technology progresses to smaller and smaller device dimensions, the effect of quantum confinement becomes more and more important. The density-gradient theory provides a computationally efficient method to include quantum confinement in the conventional drift-diffusion formulation commonly used for semiconductor device physics simulation. In the first post of a two-part blog series, we will briefly review the theory, with a focus on the equations implemented in the Semiconductor Module. Interested readers may consult Ref. 1 for a comprehensive review.

Electrostatics and Charge Carrier Conservation

We start by reviewing the fundamental equations common to the conventional drift-diffusion theory and the density-gradient (DG) theory. In the next section, we will discuss the difference between them.

Semiconductor device physics simulation concerns the transport of charge carriers driven by the electric field (drift) and the carrier concentration gradient (diffusion). The electric field \mathbf{E} (V/m) is given by the electrostatic equations under the quasistatic assumption:

(1)

\nabla \cdot (\epsilon \mathbf{E})=\rho \qquad \mathbf{E}= – \nabla V

where\epsilon is the permittivity (F/m), \rho is the charge density (C/m³), and V is the electric potential (V).

The charge density, \rho, is given by the contributions from the hole, electron, and ionized donor and ionized acceptor concentrations (1/m³) (~p, n, and ~N_d^{+} and N_a^{-}, respectively).

(2)

\rho=q(p-n+N_d^{+}-N_a^{-})

where q is the elementary charge (C).

The electron and hole concentrations (n and ~p) are governed by conservation laws given in the form of current continuity equations:

(3)

\frac{\partial n}{\partial t}&=&\frac{-\nabla\cdot\mathbf{J}_n}{-q}-Re_n+Ge_n

\frac{\partial p}{\partial t}&=&\frac{-\nabla\cdot\mathbf{J}_p}{+q}-Re_p+Ge_p

where t is the time (s); \mathbf{J}_n and \mathbf{J}_p are the electron and hole current densities (A/m²); and Re_n, Ge_n, Re_p, Ge_p are the recombination and generation rates for the electrons and holes (1/m³/s).

The electron and hole current densities, \mathbf{J}_n and \mathbf{J}_p, can be expressed in terms of the quasi-Fermi levels, E_{fn} and E_{fp} (V):

(4)

\mathbf{J}_n &=& q~n~\mu_n \nabla E_{fn} + q~n ((E_c-E_{fn})\mu_n + Q_n) \nabla T /T \\

\mathbf{J}_p &=& -q~p~\mu_p \nabla E_{fp} – q~p ((E_v-E_{fp})\mu_p + Q_p) \nabla T /T

where \mu_n and \mu_p are the electron and hole mobilities (m²/V/s), E_c and E_v are the conduction band and valance band edges (V), Q_n and Q_p are the electron and hole nonequilibrium contributions to the thermal diffusion coefficient (m²/s), and T is the temperature (K).

The conduction band edge, E_c, and the valance band edge, E_v, are related to the electric potential, V; electron affinity, \chi; and band gap, E_g:

(5)

E_c = -V-\chi \qquad E_v = E_c – E_g

Note that in the implementation for the Semiconductor Module, all of the energy level variables (E_{fn},~E_{fp},~E_c,~E_v,~\chi,~E_g) are scaled by the elementary charge to take on the unit of electric potential (V).

The Equations of States for the Drift-Diffusion and DG Theories

Both the conventional drift-diffusion theory and the DG theory obey the electrostatics and current conservation laws listed above. The difference between them lies in the equation of states of the electron and hole gases.

In the drift-diffusion theory, the charge carrier concentrations n and ~p are related to the quasi-Fermi levels E_{fn} and E_{fp} via the following equations:

(6)

n&=&N_c F_{1/2}(\frac{E_{fn}-E_c}{k_B T/q})\\

p&=&N_v F_{1/2}(\frac{-E_{fp}+E_v}{k_B T/q})

where N_c and N_v are the conduction band and valance band effective density of states (1/m³), F_{1/2} is the Fermi–Dirac integral, and k_B is the Boltzmann constant (J/K).

In contrast, the DG theory adds a contribution from the gradients of the concentrations to the equation of states via the quantum potentials V^{DG}_n and V^{DG}_p (V):

(7)

n&=&N_c F_{1/2}(\frac{E_{fn}-E_c+V^{DG}_n}{k_B T/q})\\

p&=&N_v F_{1/2}(\frac{-E_{fp}+E_v+V^{DG}_p}{k_B T/q})

where the quantum potentials, V^{DG}_n and V^{DG}_p, are defined in terms of the density gradients:

(8)

\nabla \cdot \left( \mathbf{b}_n \nabla \sqrt{n} \right) &\equiv& \frac{\sqrt{n}}{2} V^{DG}_n\\

\nabla \cdot \left( \mathbf{b}_p \nabla \sqrt{p} \right) &\equiv& \frac{\sqrt{p}}{2} V^{DG}_p

The DG coefficients, \mathbf{b}_n and \mathbf{b}_p (V m²), are given by the inverse of the DG effective mass tensors, \mathbf{m}_n and \mathbf{m}_p (kg):

(9)

\mathbf{b}_n &=&\frac{\hbar^2}{12 q} \left[\mathbf{m}_n\right]^{-1}\\

\mathbf{b}_p &=&\frac{\hbar^2}{12 q} \left[\mathbf{m}_p\right]^{-1}

where \hbar is the Planck constant.

Solution Strategy

For the conventional drift-diffusion formulation, it is possible to explicitly insert the equation of states (6) into the fundamental equations, (1)~(5), and solve for the three dependent variables: V, E_{fn}, and E_{fp}.

For the DG formulation, the equation of states (7)~(8) forms implicit relations for the charge carrier concentrations. In this situation, it is necessary to introduce new dependent variables to solve the implicit equations. Following Ref. 1, we use the Slotboom variables \phi_n and \phi_p (V) as the additional dependent variables:

(10)

n&\equiv& exp(\frac{\phi_n}{k_B T/q})\\

p&\equiv& exp(\frac{\phi_p}{k_B T/q})

The quantum potentials can now be expressed in terms of the Slotboom variables:

(11)

V^{DG}_n &=& + E_c – E_{fn} + k_B T/q \left[ log(F_{1/2})\right]^{-1} \left( \frac{\phi_n}{k_B T/q} – log(N_c) \right)\\

V^{DG}_p &=& -E_v + E_{fp} + k_B T/q \left[ log(F_{1/2})\right]^{-1} \left( \frac{\phi_p}{k_B T/q} – log(N_v) \right)

where \left[ log(F_{1/2})\right]^{-1} is the inverse function of log(F_{1/2}(\cdot)).

The DG formulation then solves for the five dependent variables, V, E_{fn}, E_{fp}, \phi_n, and \phi_p, using the fundamental equations (1)~(5) and the equation of states (7)~(11). It is evident that this approach does not add an excessive computational burden as compared to other, more sophisticated quantum mechanical methods. Therefore, the DG theory provides an efficient alternative for engineering purposes.

Recombination Rates

Because the gradient of the carrier concentrations enters the equation of states in the DG theory, the equilibrium concentrations are no longer only a function of the equilibrium Fermi level, thus the recombination rates involve more complex evaluations (Ref. 2).

For explicit traps, the electron and hole recombination rates r_e and r_h (1/m³/s) are now computed with the formula based on the quasi-Fermi level differences:

(12)

r_e &=& n~C_n~N_t (1-f_t) (1-e^{\frac{E_{ft}-E_{fn}}{k_B T/q}})\\

r_h &=& p~C_p~N_t~f_t (1-e^{\frac{E_{fp}-E_{ft}}{k_B T/q}})

where C_n and C_p are the average capture probabilities for the charge carriers (m³/s), N_t is the trap density (1/m³), f_t is the trap occupancy (1), and E_{ft} is the quasi-Fermi level of the trap (V).

The trap occupancy is given by Fermi–Dirac statistics:

(13)

f_t = \frac{1}{1+\frac{1}{g_D}~e^{\frac{E_t-E_{ft}}{k_B T/q}}} \qquad \mbox{or equivalently} \qquad \frac{1-f_t}{f_t} = \frac{1}
{g_D}~e^{\frac{E_t-E_{ft}}{k_B T/q}}

where g_D is the degeneracy factor (1) and E_t is the trap energy level (V).

For the direct, Auger, and Shockley–Read–Hall (SRH) recombinations, the rate expressions may look familiar; however, the underlying definitions of the various parameters are more complex. For example, the SRH recombination rates R_n and R_p (1/m³/s) are:

(14)

R_n=R_p=\frac{n~p-n_{eq}^{DG}~p_{eq}^{DG}}{\tau_p(n+n_1)+\tau_n(p+p1)}

where \tau_n and \tau_p are the electron and hole lifetimes (s).

The electron and hole equilibrium concentrations, n_{eq}^{DG} and p_{eq}^{DG} (1/m³), now become:

(15)

n_{eq}^{DG} &\equiv& N_c F_{1/2}(\frac{V_{eq,adj}-E_c+V_n^{DG}}{k_B T/q})\\

p_{eq}^{DG} &\equiv& N_v F_{1/2}(\frac{-V_{eq,adj}+E_v+V_p^{DG}}{k_B T/q})

where V_{eq,adj} is the equilibrium Fermi level (V).

Note the appearance of the quantum potentials V^{DG}_n and V^{DG}_p in the expressions.

The parameters n_1 and p_1 (1/m³) are no longer constants, even in simple cases, and are computed using the original definitions:

(16)

n_1 &\equiv& n~e^{\frac{E_{t}-E_{fn}}{k_B T/q}}\\

p_1 &\equiv& p~e^{\frac{E_{fp}-E_{t}}{k_B T/q}}

Note the dependence on the carrier concentrations n and p, which implicitly depend on the concentration gradients.

Boundary Conditions for the Slotboom Variables

In most cases, the boundary conditions for the Slotboom variables, \phi_n and \phi_p, are simple natural boundary conditions (Ref. 3):

(17)

\mathbf{n} \cdot (\mathbf{b}_n \nabla \sqrt{n}) = 0 \qquad \mathbf{n} \cdot (\mathbf{b}_p \nabla \sqrt{p}) = 0

In the case of a boundary representing an abrupt potential barrier (for example, a silicon-oxide interface), Jin et. al. in Ref. 4 suggested using the Wentzel–Kramers–Brillouin (WKB) approximation to obtain the boundary condition:

(18)

\mathbf{n} \cdot (\mathbf{b}_n \nabla \sqrt{n}) = -\frac{b_{n,ox}}{d_n}\sqrt{n}

where b_{n,ox} is the b_{n} coefficient in the oxide (V m²) and d_n is the penetration depth into the oxide (m), given by:

(19)

b_{n,ox}=\frac{\hbar^2}{12~q~m_{n,ox}^\star} \qquad d_n = \frac{\hbar}{\sqrt{2~q~m_{n,ox} \Phi_{n,ox}}}

where m_{n,ox}^\star and m_{n,ox} are the effective masses in the oxide (kg) and \Phi_{n,ox} is the potential barrier height (V).

Options for Heterojunctions

Currently, in the conventional drift-diffusion formulation, the Semiconductor Module offers two options for heterojunctions: Continuous quasi-Fermi level and Thermionic emission.

The former can be easily extended to the DG formulation: Just let the quasi-Fermi level and the Slotboom variable be continuous across the heterojunction, which is automatic for Lagrange shape functions. This mimics the continuous nature of the quantum mechanical wave function, although it should be treated as phenomenology at best (Ref. 1).

The latter option assumes the thermionic emission process dominates and allows the quasi-Fermi level and the Slotboom variable to be discontinuous across the heterojunction. The same formula as the drift-diffusion theory is used for the thermionic current density, yielding similar results.

Final Remarks

The discussion above gives an overview of the DG formulation available in the Semiconductor Module. It should be emphasized that only the DG confinement theory is implemented, not the DG tunneling theory (Ref. 1). This implementation provides a computationally efficient option to include the effect of quantum confinement in device physics simulations for engineering purposes. In a follow-up blog post, we will use three example models to demonstrate the power of this simulation approach.

References

M.G. Ancona, “Density-gradient theory: a macroscopic approach to quantum confinement and tunneling in semiconductor devices,” Journal of Computational Electronics, vol. 10, p. 65, 2011.
M.G. Ancona, Z. Yu, R.W. Dutton, P.J. Vande Voorde, M. Cao, and D. Vook, “Density-Gradient Analysis of MOS Tunneling,” IEEE Transactions On Electron Devices, p. 2310, Vol. 47, No. 12, December 2000.
M.G. Ancona, D. Yergeau, Z. Yu, and B.A. Biegel, “On Ohmic Boundary Conditions for Density-Gradient Theory”, Journal of Computational Electronics 1: 103–107, 2002.
S. Jin, Y.J. Park, and H.S. Min, “Simulation of Quantum Effects in the Nano-scale Semiconductor Device,” Journal of Semiconductor Technology and Science, vol. 4, no. 1, p. 32, 2004.

Comments (0)

EXPLORE COMSOL BLOG

Intro to Density-Gradient Theory for Semiconductor Device Simulation

Electrostatics and Charge Carrier Conservation