CHARGE solver introduction

The Charge Transport (CHARGE) solver is a physics-based electrical simulation tool for semiconductor devices, which self-consistently solves the system of equations describing the electrostatic potential (Poisson’s equation) and density of free carriers (the drift-diffusion equations). The drift-diffusion model is an established and robust method that will produce accurate results for a wide range of semiconductor devices under common operating conditions.

Charge Transport Solver Physics

This section will introduce the basic mathematical and physics formalism behind the self-consistent algorithm used in CHARGE, starting from the non-linear Poisson and drift-diffusion equations.  CHARGE solves the drift-diffusion equations for electrons and holes (carriers)

\begin{aligned} \mathbf{J}_{n} &=q \mu_{n} n \mathbf{E}+q D_{n} \nabla n \\ \mathbf{J}_{p} &=q \mu_{p} p \mathbf{E}-q D_{p} \nabla p \end{aligned}

where Jn,p is the current density (A/cm2), q is the positive electron charge, μn,p is the mobility, E is the electric field, Dn,p is the diffusivity, and n and p are the densities of the electrons and holes, respectively (the subscripts n and p indicate quantities that are specific to the carrier type). Each carrier (electron or hole) moves under the influence of two competing processes: drift due to the applied electric field, and random thermal diffusion due to the gradient in the density. These processes are represented in the drift-diffusion equations as the sum of two terms.

The mobility μn,p  describes the ease with which carriers can move through the semiconductor material, and is related to the diffusivity Dn,p  through the Einstein relation,

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

where kB is the Boltzmann constant. The mobility is a key property of the material, and may be modeled as a function of temperature, impurity (doping) concentration, carrier concentration, and electric field. For further details, please see the section on material models.

To solve the drift-diffusion equations, the electric field must be known. To determine the electric field, Poisson's equation is solved:$$-\nabla\cdot\left(\varepsilon\nabla V\right)=q\rho$$

where ε is the dielectric permittivity, V the electrostatic potential (E = \(-\nabla V \) ) and ρ the net charge density,

$$\rho = p-n+C$$

which includes the contribution \(C\) from the ionized impurity density.

Finally, the auxiliary continuity equations are required to account for charge conservation

\begin{array}{l}{\frac{\partial n}{\partial t}=\frac{1}{q} \nabla \cdot \mathbf{J}_{n}-R_{n}} \\ {\frac{\partial p}{\partial t}=-\frac{1}{q} \nabla \cdot \mathbf{J}_{p}-R_{p}}\end{array}

where Rn,p is the net recombination rate (the difference between the recombination rate and generation rate). The physical processes associated with the material are assumed to act equivalently when applied to electrons or holes, and as a result, R = Rn = Rp. The recombination and generation processes are important factors in the material-specific calculation of carrier behavior. Multiple recombination and generation processes are modeled, which may depend on temperature, impurity (doping) concentration, carrier concentration, electric field, and current density. For further details, please see the section on material models.

CHARGE discretizes and solves the drift-diffusion and Poisson’s equations on an unstructured finite-element mesh in one and two dimensions. The simulation region is partitioned into multiple domains along boundaries between materials with unique physical descriptions. The materials used in the simulation may be categorized as insulators, semiconductors, or conductors; each type of material has an associated user-specified model or collection of models that describe its behavior. In particular, specialized multi-coefficient models are provided for semiconductors that describe the fundamental properties, mobility, and recombination processes inherent to that material.

The system of equations solved by CHARGE admits both a steady-state and time-varying result. By enforcing the condition

$$\frac{\partial n}{\partial t}=\frac{\partial p}{\partial t} = 0$$

in the continuity equations, the carrier density and electrostatic potential can be solved at steady-state. Steady-state simulations can be used to examine the system’s behavior at a fixed operating point, and are also useful when extracting small-signal parameters for a component (e.g. for frequency response analysis). Alternately, by specifying an initial condition for the carrier density and electrostatic potential, the equations can be solved in a sequence of discrete times. The time-dependent behavior of the component can then be used to directly evaluate its large-signal time-domain response or extract large-signal AC parameters.

Boundary conditions are very important in an accurate semiconductor device simulation. Two categories of boundary condition are present in CHARGE: those that relate to the electrostatic potential (Poisson’s equation) and those that relate to the carrier densities (the drift-diffusion equations). Poisson’s equation and the drift-diffusion equations are second-order partial differential equations (PDE), and each requires that the solution be explicitly specified for at least one location. This is known as a Dirichlet boundary condition. For the electrostatic potential, the Dirichlet condition takes the form of a boundary (internal or external) with a fixed voltage specified,

$$V(x) = V_1$$

as is typical of an electrical contact. For the carrier densities, the majority carrier concentration is set to its equilibrium value at the interface between a contact and the semiconductor, such that

$$p-n+C=0$$

At internal boundaries between two domains that are not contacts, the boundary conditions are determined from the physical properties of the interface. In the case of the electrostatic potential, the electric flux density must be continuous across the boundary in the absence of a surface charge. For the electron and hole densities, boundary conditions between the semiconductor and adjacent materials may be specified in terms of a surface recombination current density, which will default to zero (no carrier flux across the boundary) for insulators or an infinite recombination velocity (forcing the carrier density to its equilibrium value) for contacts. At the external (open) boundaries of the simulation domain, homogenous Neumann boundary conditions are applied: the electric field normal to the boundary is set to zero as is the surface recombination current density. Physically, this corresponds to an insulating boundary across which no charge can flow.

At boundaries between different semiconductor materials and between contacts and semiconductors, the solver also accounts for thermionic emission and tunneling over the barrier as well as field-induced barrier lowering based on the model in reference [1]. However, certain limitations apply:

• Tunneling constant is calculated internally from material properties and can't be adjusted.

• Tunneling model is local and derived assuming triangular barriers, representative of those at Schottky contacts. The model may not accurately capture thin barriers and barriers of different shapes, such as thin oxides, quantum wells, and other semiconductor heterostructures.

Small Signal AC Analysis

Small signal AC analysis gives a frequency domain solution to the linearized charge transport equations (Poisson and drift-diffusion). The CHARGE transport solver will find a self-consistent solution to the coupled, non-linear system of equations.

\begin{array}{c}{-\nabla \cdot\left(\varepsilon_{r} \nabla V\right)-\frac{q}{\varepsilon_{0}}(p-n+C)=F_{V}(V, n, p)=0} \\ {-\frac{\partial n}{\partial t}+\nabla \cdot \mathbf{J}_{n}-R=-\frac{\partial n}{\partial t}+F_{n}(V, n, p)=0} \\ {\frac{\partial p}{\partial t}+\nabla \cdot \mathbf{J}_{p}+R=\frac{\partial p}{\partial t}+F_{p}(V, n, p)=0}\end{array}

At steady-state, the time derivatives \(\frac{\partial n}{\partial t}\) and \(\frac{\partial p}{\partial t}\) are zero, and solution \(\bar{x} = (\bar{V}, \bar{n}, \bar{p})\) is determined. The function \(F_{V,n,p}\) can be expanded in a Taylor series around the operating point \(\bar{x}\) with a perturbation \(\widetilde{x} e^{j\omega t}\)(in response to an applied perturbation \(q=\bar{q}+\widetilde{q} e^{j\omega t}\) and truncated at first order, e.g.

\begin{aligned} F_{V}(x) &=F_{V}(\overline{x})+\left.\frac{\partial F_{V}}{\partial x}\right|_{x}(x-\overline{x})+\ldots \\ & \approx F_{V}(\overline{x})+\left(\left.\frac{\partial F_{V}}{\partial V}\right|_{\overline{x}} \widetilde{V}+\left.\frac{\partial F_{V}}{\partial n}\right|_{\overline{x}} \widetilde{n}+\left.\frac{\partial F_{V}}{\partial p}\right|_{\overline{x}} \widetilde{p}\right) e^{j \omega t} \end{aligned}

The steady-state residual \(F_{V}(\bar{x})\) is zero. Incorporating the time derivatives, we obtain the system of equations

$$\left.\frac{\partial F_{V}}{\partial V}\right|_{\overline{x}} \widetilde{V}+\left.\frac{\partial F_{V}}{\partial n}\right|_{\overline{x}} \widetilde{n}+\left.\frac{\partial F_{V}}{\partial p}\right|_{\overline{x}} \widetilde{p}=\widetilde{q}_{V}$$

$$-j \omega \tilde{n}+\left.\frac{\partial F_{n}}{\partial V}\right|_{\overline{x}} \widetilde{V}+\left.\frac{\partial F_{n}}{\partial n}\right|_{\overline{x}} \widetilde{n}+\left.\frac{\partial F_{n}}{\partial p}\right|_{\overline{x}} \widetilde{p}=\widetilde{q}_{n}$$

$$j \omega \widetilde{p}+\left.\frac{\partial F_{p}}{\partial V}\right|_{\overline{x}} \widetilde{V}+\left.\frac{\partial F_{p}}{\partial n}\right|_{\overline{x}} \widetilde{n}+\left.\frac{\partial F_{p}}{\partial p}\right|_{\overline{x}} \widetilde{p}=\widetilde{q}_{p}$$

where we assume that the perturbations to the carrier densities at the terminals are zero (\(\widetilde{q}_{n,p}=0\)) and \(\widetilde{q}_{n,p}\) represents the applied small amplitude signal. This system of equations is solved for the (complex) small amplitude responses \(\widetilde{V}\), \(\widetilde{n}\) and \(\widetilde{p}\).

The terminal currents are calculated as the sum of their components, including the displacement and carrier currents. The small amplitude displacement current density is

$$\mathbf{J}_{d}=-j \omega \varepsilon \nabla \widetilde{V}$$

The electron and hole currents are found by assuming a first-order expansion around the operating point, e.g.

$$\tilde{\mathbf{J}}_{n}(\widetilde{x}) e^{j \omega t}=\mathbf{J}_{n}(x)-\mathbf{J}_{n}(\overline{x})=\left[\left.\frac{\partial \mathbf{J}_{n}}{\partial \varphi_{n}}\right|_{\overline{x}} \widetilde{\varphi}_{n}+\left.\frac{\partial \mathbf{J}_{n}}{\partial \varphi_{p}}\right|_{\overline{x}} \widetilde{\varphi}_{p}+\left.\frac{\partial \mathbf{J}_{n}}{\partial u}\right|_{\overline{x}} \widetilde{u}\right] e^{j \omega t}$$

The net currents are resolved as the integrated flux into the terminals.

The small signal analysis is a computationally effective way of resolving the frequency domain response of a semiconductor device under conditions where the small amplitude/linear approximation is valid. For large signal response when the behavior of the semiconductor device is non-linear, a full time domain simulation can be performed.

Coupled Heat Charge Transport

To extend the charge transport solver to account for self-heating effects, the drift-diffusion equations are adapted to account for the influence of a gradient in the temperature,

\begin{array}{l}{\mathbf{J}_{n}=\mu_{n}\left(q \mathbf{E}_{n}+k \nabla T\right) n+\mu_{n} k T \nabla n} \\ {\mathbf{J}_{p}=\mu_{p}\left(q \mathbf{E}_{p}+k \nabla T\right) p+\mu_{p} k T \nabla p}\end{array}

where \(\mu\)is the carrier mobility, \(k\) is the Boltzmann constant, \(q\) is the elementary charge, and where the driving fields are described as

\begin{array}{l}{\mathbf{E}_{n}=\frac{1}{q} \nabla E_{c}-\frac{k T}{q} \nabla \log N_{c}+\frac{3}{2} \frac{k T}{q} \nabla \log T=\frac{1}{q} \nabla E_{c}-\frac{3}{2} \frac{k T}{q} \nabla \log \left(\frac{m_{n} *}{m_{0}}\right)} \\ {\mathbf{E}_{p}=\frac{1}{q} \nabla E_{v}-\frac{k T}{q} \nabla \log N_{v}-\frac{3}{2} \frac{k T}{q} \nabla \log T=\frac{1}{q} \nabla E_{v}+\frac{3}{2} \frac{k T}{q} \nabla \log \left(\frac{m_{p} *}{m_{0}}\right)}\end{array}

In these equations,\(E_c\) and \(E_v\) are the conduction band and valence band energies, respectively;  \(N_c\) and \(N_v\)  are the electron and hole effective densities of states, respectively; and\({m_n}^*\) and \({m_p}^*\) are the electron and hole effective masses, respectively, scaled by the electron rest mass\(m_0\). An additional term is added when Fermi-Dirac statistics are used. In the absence of a temperature gradient, these equations reduce to the classic drift-diffusion equations.

In a coupled heat and charge transport (CHCT) simulation, the drift-diffusion equations defined above are solved self-consistently with Poisson’s equation and the heat transport equation (see heat transport solver physics description),

$$c_p \rho \frac{\partial T}{\partial t} -\nabla\cdot\left(k\nabla T\right) = Q$$

where the applied heat \(Q\)(W/m3) is now generated from the combination of Joule heating from the carrier currents and from carrier recombination:

$$Q=Q_n+Q_p+Q_R$$

In this description of heating,

\begin{array}{r}{Q_{n, p}=\mathbf{J}_{n, p} \cdot \mathbf{E}_{n, p}} \\ {Q_{R}=q\left(E_{g}+3 k T\right) R}\end{array}

where \(E_g\) is the bandgap energy and \(R\) is the net recombination (recombination less generation).

The boundary conditions for the heat transport equation are set by the temperatures at the contacts, which are assumed to be constant and uniform. All other thermal boundaries are assumed to be insulating.

Meshing

CHARGE uses an unstructured, finite-element mesh, like the one shown in the following screenshot. The solution to the system of equations used to determine the physical quantities of interest is estimated from the discrete formulation of those equations. Consequently, it's important to understand that of the fundamental simulation quantities (material properties, geometry information, electrostatic potential, and carrier concentrations) are calculated at each mesh vertex.

A finer mesh (with shorter edge lengths and smaller elements) will better approximate the exact solution to the system of equations, but at a substantial cost. As the mesh features become smaller, the simulation time and memory requirements will increase.  CHARGE provides a number of tools, including the automatic and guided mesh refinement, that allow you to obtain accurate results, while minimizing computational effort.

Units and normalization

References

1. K. Matsuzawa, K. Uchida, and A. Nishiyama, "A unified simulation of Schottky and ohmic contacts", IEEE Transactions on Electron Devices, Vol. 47, No. 1, January 2000, pp. 103-108.

Selected Application Examples

• Getting Started Examples
• Vertical Photodetector
• Traveling Wave Modulator
• Thin film GaAs solar cell
• GaAs PIN diode
• MOSFET

Video resources

• Charge overview (YouTube.com)
• More videos (lumerical.com)