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.
• MQW DOMAINS: To enable CHARGE-MQW coupled simulation, select the name of the geometry structure group that contains all the quantum wells and barriers in the multi-quantum well stack. This defines the layers where MQW equations should be solved and coupled to the CHARGE equations.

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.
• DISABLE INTERFACE MESH AUTOREFINEMENT: Enabling this option will disable automatic refinement of the mesh at interfaces between different partitioned domains. This will usually lead to a smaller number of mesh vertices and faster simulation, which is useful in large geometry devices, such as planar layered devices with very large width or area, where the variation of electrical quantities, such as electric field and carrier density, is across larger distances.

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.

Advanced Options

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

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

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.

MQW tab

This tab provides the user interface for the MQW solver options. These options are similar to the user interface options in the standalone MQW solver object. For more details, please check the MQW solver documentation page. There are the following differences:

• The MQW layer table is auto-populated from the geometry information defined in MQW domains under General tab. Only the strain and valence band offset columns are editable.
• The layer table is auto-partitioned by default, such that there is one quantum well per partition.
• The Parameters tab is removed since the parameters are provided by coupling to CHARGE.

Advanced tab

• OVERRIDE PARTITION TABLE: If this option is checked the partition table becomes editable.
• SPATIAL REFINEMENT: MQW is a 1D solver (one-dimensional quantum confinement). When coupled to 2D CHARGE, MQW has to be solved at different lateral locations (normal to the growth direction). This slider enables adjusting the trade-off between simulation speed and simulation accuracy. Moving the slider to the right enables higher accuracy at the expense of higher computational cost by solving MQW at more closely spaced lateral locations. Moving the slider to the left enables faster simulation at the expense of lower accuracy by solving MQW on a coarser lateral grid.

Results returned

Mesh

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.

Spatial

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.