Simulation objects
CHARGE solver - Simulation object
General tab
The charge transport solver contains several settings:
SIMULATION REGION: Select the simulation region the solver will use.
NORM LENGTH (UM): The length of the device in the direction perpendicular to the plane of the simulation; any normalizations length will be with respect to this value.
SOLVER MODE: The color of the simulation region will change depending on which option is picked. Also, the available options in the contacts table will change accordingly. Please visit the Solvers section for more information.
STEADY STATE - dc simulations
TRANSIENT - time dependent simulations
SSAC - Small signal AC simulations
TEMPERATURE DEPENDENCE: The temperature dependence can be considered in the CHARGE solver, below is the choice of temperature dependence. Please visit the Solvers section for more information.
ISOTHEMAL - An electrical simulation mode with a uniform temperature applied to all contacts and the simulation domain. Temperature-dependent material properties are updated with parameters calculated for the supplied simulation temperature.
NON-ISOTHEMAL - An electrical simulation mode that allows a spatially varying temperature map added to a region specified by the user. In regions where no temperature information is present, the simulation region temperature is used. Contact temperatures are determined from the average temperature on the contact surface. Heating is not calculated, and no updates to the temperature profile are made.
COUPLED - A self-consistent CHARGE-HEAT coupled simulation mode that solves Poisson's equations, the drift-diffusion equations and the heat transport equations. The contact temperatures are specified as boundary conditions (fixed temperature on a surface), and heat sources can be applied (e.g. from an external import). Heating is calculated from the ohmic loss (\( \overrightarrow{J}\cdot\overrightarrow{E} \)) and from recombination.
SIMULATION TEMPERATURE (K): The temperature in Kelvin at which the simulation will be done. This is only effective when Isothermal Temperature dependence is selected.
Continue from previous solutions
ENABLE CONTINUATION: If this option is checked, the user has the choice to specify a file that has already been run with a solution. This solution will be used as the starting point for the Newton solver. In simulations where the starting point isn't at zero volts, this can be very useful to save time.
Mesh tab
Global Mesh Constraints
The global mesh settings section contains several settings:
MINIMUM EDGE LENGTH: The minimum length of an edge of a triangle (in 2D) or tetrahedron (3D) used in the mesh.
MAXIMUM EDGE LENGTH: The maximum length of an edge of a triangle (in 2D) or tetrahedron (3D) used in the mesh.
Auto Refinement Settings
The auto refinement settings section contains several settings:
MAX REFINE STEPS: The automatic refinement proceeds in multiple stages, creating a quality triangulation and refining the mesh according to the change in doping density and, if present, optical generation rate. This setting limits the number of refinements at each stage, and corresponds to the number of vertexes that can be added to the mesh at each stage.
SENSITIVITY: This setting controls the threshold at which the mesh will be refined due to the gradient in the doping density or optical generation rate. The default value will roughly correspond to a limit of a factor of 2 change in the doping density or generation rate over the span of an element in the mesh.
Geometry Options
DEFLECTION TOLERANCE (UM): Controls how curved surfaces are broken up into multiple linear segments. A smaller deflection tolerance will force the geometry builder to break up a curved surface into smaller segments. Only available in 2D simulations and 3D simualtions with vertex insertion volume meshing.
MAX PLC EDGE LENGTH (UM): Forces the solver to break long objects into smaller segments in order to help with the geometry building with long and narrow structures. This is only available in 3D simulations with vertex insertion volume meshing.
Mesh Options
TRIANGLE QUALITY: The quality of mesh triangle, defined the minimum angle used in the triangular mesh cells. The higher the angle, the higher the quality of the triangle.
VOLUME MESHING: The meshing algorithm used by the solver for 3D simulations
Vertex insertion: This method forms vertices within boundaries of the volume to generate a 3D mesh and refines the mesh in a volume based on variations in parameters such as doping and optical generation across the volume. This is the default option and suitable for most cases.
Hybrid: In addition to vertex insertion, this method uses a surface meshing technique which can result in a better mesh than vertex insertion only method for thin geometry features. It should only be used if vertex insertion method fails generating the mesh and the accuracy of the results should be verified.
Transient tab
Transient Simulation Controls
The transient simulation controls section contains several settings:
MIN TIME STEP: The smallest time step that will be used in the transient simulation (used as the initial time step)
MAX TIME STEP: The largest time step that will be used in the transient simulation.
ABS LTE LIMIT: If the absolute error is less than this value then the solver increases the time step.
REL LTE LIMIT: If the relative error is less than this value then the solver increases the time step.
TIME INTEGRATION: Set the time stepping scheme
BDF2 - Backward Differentiation Formula (2nd order) method.
TR-BDF2 - Trapezoidal Rule/Backward Differentiation Formula (2nd order) method.
Downsampling
DOWN SAMPLE MODE: The type of downsampling to use, None, for no downsampling, count, for specifying the number of point to downsample at, and interval to specify the downsample step size.
Global Source Shutter
The global source shutter will apply a shutter to all the optical generation (source) objects in the simulation.
SHUTTER MODE: disabled, for no shutter, step on and step off for step functions, pulse on and pulse off for a pulse with on and off times. The time shutter function will be plotted as the option is chosen for ease of use. The on and off times can then be specified.
SHUTTER TON: Sets the on time of the shutter.
SHUTTER TOFF: Sets the off time of the shutter.
SHUTTER TSLEW: Sets a slew in the stepping of the shutter. This feature is useful in helping the solver to converge where a sharp step in the shutter tend to makes the solver diverge.
SHUTTER SLEW FUNCTION: linear for a linear transition between off and on states, log for an exponential change between the off and on states.
SHUTTER SLEW CUTOFF: The shutter uses the logarithmic function to step from this value to 1 (for turning on) and from 1 to this value (for turning off). The stepping from zero to the slew cutoff value (and back) is abrupt (within one time step).
SOURCE POWER SCALING: This option can be used to scale the amplitude of the source. By default the value is set to '1' which means no scaling.
Small Signal AC tab
PERTURBATION AMPLITUDE (V): Amplitude of the SSAC perturbation
SMALL SIGNAL AC OPTICAL GENERATION: Adds a small signal AC component to the import optical generation rate object (OGR) for frequency response analysis of illuminated devices. It is a global option that applies to all such objects in the simulation. If there is an import OGR object and this option is checked the OGR object will be used as a small signal source in the small signal AC simulation (CHARGE solver mode ssac).
FREQUENCY SPACING: It has three modes for frequency spacing of the SSAC
SINGLE - A single frequency
LINEAR - Linear sampling for a range of frequency
LOG - Logarithm sampling for a range of frequency
Results tab
This tab contains a list of all the spatial results that can be recorded throughout the simulation. One can pick to enable or disable one or more of the results to save memory as needed. The tab also contains basic information about the available results such as their units and descriptions. Please see the results section below for more information on the available results.
Advanced tab
Note: Changes to advanced settings should be made only in case the simulation does not converge or the convergence is too slow.
Non-Linear Solver Settings
SOLVER TYPE: the charge transport solver simultaneously solves the equations for the electrostatic potential, charge and heat transport (if licensed and configured). The solutions to these equations must be self-consistent, i.e. the charge calculated from the drift-diffusion equations satisfies the Poisson equation, and vice versa. Two common approaches are used to solve this system of equations: Gummel's method and Newton's method.
GUMMEL: decouples the charge problem from the electrostatic potential problem at each step. As both equations are non-linear, they must in turn be solved using an iterative or direct method. Using Gummel's method, the electrostatic potential is first solved holding the charge fixed. Next, this solution to the electrostatic potential is used as a fixed input to the charge equations, which are updated. This process continues until the solution is self-consistent. Gummel's method is stable and efficient for devices where the currents are small and the variations in the charge distribution are not too extreme (the charge and electrostatic potential are weakly coupled problems).
NEWTON: the system of equations for the electrostatic potential, charge and heat transport (if licensed and configured) are solved simultaneously using Newton's method. Newton's method requires a good initial guess to aid convergence (see initialization), which may diverge if the initial guess is too far from the solution. However, Newton's method can handle devices where the variations in charge density are large.
DC UPDATE MODE: self-consistent (default) solves the system of equations for the electrostatic potential (Poisson equation), charge density (drift-diffusion), and heat transport. The electrostatic option will only solve Poisson's equation, using the carrier density as a fixed input; the charge option will only solve the drift-diffusion equations, using the electrostatic potential as a fixed input. Both the electrostatic option and the charge option require an existing solution that provides the fixed input (see continuation).
FERMI STATISTICS: Electrons and holes are fermions, and therefore obey Fermi-Dirac statistics. At a finite temperature, the energy distribution of the electrons is described by the Fermi-Dirac function,
$$ f(E)=\frac{1}{1+\exp\left(\frac{E-E_{f}}{kT}\right)} $$
where \( k \) is the Boltzmann constant, \( T \) is the temperature, and \( E_{f} \) is the Fermi energy. When \( \left|E_{c,v}-E_{f}\right|>>kT \), the Fermi-Dirac distribution can be approximated by the Boltzmann distribution,
$$ f(E)=\exp\left(-\frac{E-E_{f}}{kT}\right) $$
When the net doping density is non-degenerate in a semiconductor device, the Fermi energy is located within the band gap. When the condition \( \left|E_{c,v}-E_{f}\right|>>kT \) is satisfied (typically \( \left|E_{c,v}-E_{f}\right|>>3kT \)), the Fermi-Dirac distribution can be replaced with the Boltzmann distribution.
The carrier density is calculated from the integrated product of the Fermi-Dirac distribution (probability of occupancy) and the density of states (available states to occupy). For electrons, the equation is
$$ n=\int_{-\infty}^{+\infty} f(E)N(E)dE $$
When the Fermi-Dirac distribution is used, this integral does not have an analytic solution, and must be approximated numerically. However, for non-degenerate conditions, the Fermi-Dirac distribution is approximated by the Boltzmann distribution, and the preceding equation reduces to
$$ n=N_{c}\exp\left(-\frac{E-E_{F}}{kT}\right) $$
which describes the electron density at equilibrium (\( N_{c} \) is the effective density of states and is a constant related to the specific properties of the semiconductor).
MULTITHREADING: If enabled, the user can limit the number of threads used by the non-linear solver, otherwise, the solver will use all available processors.
Gummel's Method
These settings apply when the SOLVER TYPE is set to "gummel".
USE GLOBAL ITERATION LIMIT: use the global iteration limit (see Convergence Controls) to limit the number of iteration that each solver (electrostatics, charge, and heat) will attempt before determining that it has failed to converge.
DDS ITERATION LIMIT: iteration limit for the drift-diffusion solver when it is executed independently as part of a Gummel's method simulation. Takes a value between 1 and 10000.
POISSON ITERATION LIMIT: iteration limit for the electrostatic solver when it is executed independently as part of a Gummel's method simulation. Takes a value between 1 and 10000.
HT ITERATION LIMIT: iteration limit for the heat transport solver when it is executed independently as part of a Gummel's method simulation. Takes a value between 1 and 10000.
Convergence Controls
GLOBAL ITERATION LIMIT: limit the total number of iterations that can be performed as part of a non-linear solve. This may apply to the number of Newton iterations when the SOLVER TYPE is specified as "newton" or the number of Gummel iterations (each consisting of sub-iterations by specific solvers) when the SOLVER TYPE is "gummel."
CONVERGENCE METHOD: choice of "limit update," "line search cubic," or "line search bank rose". If the first option is selected, the update calculated in a Newton iteration will be limited to a maximum value (see Update Limiting) to prevent divergence when the estimate is far from the solution. If a "line search" option is selected, the corresponding line search method is used to control the size of the update:
For cubic: Numerical Recipes in C: The Art of Scientific Computing, Cambridge University Press (ISBN 0-521-43108-5)
For bank rose: R.E. Bank and D.J. Rose, Numer. Math. 37, 279-295 (1981).
NR RATIO: Newton-Richardson ratio, value between 0 and 1. This parameter defines the limiting value for the ratio of the residuals calculated at sequential Newton iterations. If the residual norm is decreased by more than this ratio \( (||r^{(k+1)}||< \) NR RATIO \( ||r^{(k)}||) \), the Jacobian is not refactored. If 0, the option is disabled.
GRADIENT MIXING: choice of "disabled", "fast", and "conservative." To aid convergence when gradient quantities change discontinuously between iterations, the quantity from the next iteration can be gradually mixed with the converged value from the previous iteration. This can improve convergence for simulations with field-dependent material properties.
Update Limiting
These parameters apply when the CONVERGENCE METHOD is chosen to be "limit update".
USE DEFAULT UPDATE LIMITS: when selected, used default values for maximum update size. The default values for each limit are listed below.
DDS MAX UPDATE: The maximum change that will be applied to the quasi-Fermi level estimates at the next non-linear iteration. This value is in multiples of the thermal voltage \(\left(V_{th}=\frac{kT}{q}\right)\), and the default is 5.
POISSON MAX UPDATE: The maximum change that will be applied to the electrostatic potential estimate at the next non-linear iteration. This value is in multiples of the thermal voltage \(\left(V_{th}=\frac{kT}{q}\right)\), and the default is 20.
HT MAX UPDATE: The maximum change that will be applied to the temperature estimate at the next non-linear iteration. This value units of Kelvins (K), and the default is 100.
Convergence Criteria
CONVERGENCE CRITERIA: choice of "update," "residual" or "any." If "update" is selected the simulation will be deemed to be converged when the maximum update is less than the absolute or relative tolerance. When "residual" is selected, the simulation is deemed to be converged when the residual norm is less than the absolute tolerance. When "any" is selected, the simulation is deemed to be converged when any of the convergence criteria are met.
UPDATE ABS TOL specifies the absolute convergence tolerance compared with the maximum update between non-linear iteration steps:
$$ x^{(k+1)}_{i}-x^{k}_{i}<\epsilon_{a}, $$
where \( x \) is the scaled variable, \( i \) is a vertex index in the mesh, and \( \epsilon_{a} \) is the absolute tolerance.
UPDATE REL TOL specifies the relative convergence tolerance compared with the maximum update between non-linear iteration steps:
$$ \frac{x^{(k+1)}_{i}-x^{(k)}_{i}}{x^{(k)}_{i}}<\epsilon_{r}, $$
where \( \epsilon_{r} \) is the absolute tolerance. If \( x^{(k)}_{i}=0 \) the relative error is set to 1 (i.e. the criterion becomes \( 1<\epsilon_{r} \) ).
The total convergence criterion for update is then:
$$ \left|\left|\min\left(\frac{x^{(k+1)}-x^{(k)}}{\epsilon_{a}},\frac{x^{(k+1)}-x^{(k)}}{\epsilon_{r}x^{(k)}}\right)\right|\right|_{\infty}<1. $$
RESIDUAL ABS TOL specifies the absolute convergence tolerance compared with the Euklidean norm (i.e. second norm) of the residual at each non-linear iteration step. In this case, the residual is scaled such that
the units of entries in the electrostatic part are fC/um in 2D and fC in 3D
the units of entries in the drift-diffusion part are mA/um in 2D and mA in 3D
the units of entries in the heat transport part are mW/um in 2D and mW in 3D
Length in the 2D cases refers to the norm length.
Initialization can be used to generate an improved initial guess when the operating point of the device is far from equilibrium. The initialization algorithm proceeds by incrementally increasing the applied bias voltage from equilibrium and computing an incomplete solution at each step.
ENABLE INITIALIZATION: enable the initialization algorithm
INIT STEP SIZE: specifies the size of the bias voltage increment used during initialization. Larger values result in faster initialization calculations but may generate a lower quality initial guess.
Results returned
Once the simulation region has been meshed, the CHARGE simulation object returns the mesh grid as an unstructured dataset which shows the returns the charge concentration, area of the mesh cells, and the unique ID of the partitioned regions.
Doping intrinsic and extrinsic \( cm^{-3} \)
Area of the finite element mesh cell
Distinct ID number given to partitioned regions
After the simulation is run the grid object will be overwritten, but the grid can be viewed from any of the new results. If the simulation is transient or ssac the results will be a function of time. In steady state each result will have the electrical boundary conditions as parameters. Depending on the the simulation type the CHARGE object will return a slightly different set of results, and choosing whether to save these or not can be done from the results tab. Each result is returned as an unstructured dataset.
bandstructure
$$eV$$
Conduction band
Efn
Electron quasi Fermi level
Efp
Hole quasi Fermi level
Intrinsic energy
Valence band
Evac
Vacuum potential
xfrac
a.u.
Alloy composition fraction
$$ A/cm^2$$
Electron current density
Hole current density
free_charge
$$ cm^{-3} $$
Net free charge density (p-n)
space_charge
Total charge density (ND-NA+p-n)
Electron density
Hole density
Net doping density
Acceptor doping density
Donor doping density
$$V/m$$
$$ V$$
Electrostatic potential
$$m$$
Electron diffusion length (SRH)
Hole diffusion length (SRH)
$$cm^2/V \cdot s$$
$$cm^2/V\cdot s$$
Hole mobility
Goptext
$$ 1/(cm^3 \cdot s) $$
External optical generation rate
Auger recombination rate
Rbbt
Band-to-band tunneling rate
Rii
Impact ionization recombination rate
Rnet
Net recombination rate
Ropt
Radiative recombination rate
Rsrh
SRH recombination rate
Rstim
Stimulated recombination rate
Rsurf
Surface recombination rate
taun
$$ s $$
Electron lifetime (SRH)
taup
Hole lifetime (SRH)
$$ W/m^3 $$
$$ K $$
Boundary Conditions
Each electrical boundary condition is returned as matrix dataset with the same name as the boundary condition set in the folder. The results will be a function of time for ssac/transient simulations or with the set points as parameters for steady state simulations.
"electrical_BC_name"
$$ A $$
Electron current
Hole current
Displacement current
$$ V $$
Source voltage
Contact voltage
Like the spatial charge results the thermal results will have time as a parameter for transient simulations, and with the boundary conditions as parameters for steady state simulations. Rather than having a result for each thermal boundary condition, as in the charge results, the thermal solver returns one matrix dataset called boundaries with the results from each boundary as an attribute.
electrostatic
$$ J/(K\cdot kg) $$
$$ W/(K \cdot m) $$
rho
$$ kg/m^3 $$
Mass density
P_"BC_name"
$$ W $$
A_"BC_name"
$$ m^3 $$
I_"elec_BC_name"
Current (electrical BC only)
CHARGE Product Reference Manual
CHARGE solver introduction
Planar silicon solar cell
CMOS image sensor - Angular response 3D