The Charge Transport (CHARGE) solver is a physicsbased electrical simulation tool for semiconductor devices, which selfconsistently solves the system of equations describing the electrostatic potential (Poisson’s equation) and density of free carriers (the driftdiffusion equations). The driftdiffusion 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 selfconsistent algorithm used in CHARGE, starting from the nonlinear Poisson and driftdiffusion equations. CHARGE solves the driftdiffusion 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 driftdiffusion 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 driftdiffusion 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 = pn+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 materialspecific 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 driftdiffusion and Poisson’s equations on an unstructured finiteelement 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 userspecified model or collection of models that describe its behavior. In particular, specialized multicoefficient 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 steadystate and timevarying 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 steadystate. Steadystate simulations can be used to examine the system’s behavior at a fixed operating point, and are also useful when extracting smallsignal 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 timedependent behavior of the component can then be used to directly evaluate its largesignal timedomain response or extract largesignal 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 driftdiffusion equations). Poisson’s equation and the driftdiffusion equations are secondorder 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
$$pn+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 fieldinduced 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 driftdiffusion). The CHARGE transport solver will find a selfconsistent solution to the coupled, nonlinear system of equations.
\begin{array}{c}{\nabla \cdot\left(\varepsilon_{r} \nabla V\right)\frac{q}{\varepsilon_{0}}(pn+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 steadystate, 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 steadystate 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 firstorder 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 nonlinear, a full time domain simulation can be performed.
Coupled Heat Charge Transport
To extend the charge transport solver to account for selfheating effects, the driftdiffusion 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 FermiDirac statistics are used. In the absence of a temperature gradient, these equations reduce to the classic driftdiffusion equations.
In a coupled heat and charge transport (CHCT) simulation, the driftdiffusion equations defined above are solved selfconsistently 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, finiteelement 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
Quantity  Description  Units  Unit description 

Length 
Length units in semiconductor models. This is reflected in the semiconductor device literature, and in the parameter coefficients for the material models. 
cm 
Centimeter 
Energy 
The electron energy E is related to the local electrostatic potential (voltage) as E = qV. All energies (and voltages) are referenced from the (equilibrium) Fermi level of an electrical contact in the system. 
eV 
Electron volts 
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. 103108.
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)