1. Poisson solvers
QTCAD contains two Poisson solvers, a linear Poisson solver (which solves the Poisson equation usually introduced in basic electricity classes), and a nonlinear Poisson solver which selfconsistently accounts for the carrier densities in classical reservoirs.
Linear Poisson solver
The linear Poisson solver solves the partial differential equation
where \(\varepsilon\) is the dielectric permittivity of the medium, \(\varphi\) is the electric potential, and \(\rho\) is the (free) charge density.
Relevant device attributes
Parameter 
Symbol 
QTCAD name 
Unit 
Default 
Setter 
Permittivity 
\(\varepsilon\) 

F/m 
11.8 \(\varepsilon_0\) 

Charge density 
\(\rho\) 

C/\(\mathrm m^3\) 
0 
Variable 
Symbol 
QTCAD name 
Unit 
Potential 
\(\varphi\) 

V 
See Parameter setters for instructions on how to employ device parameter setters.
Another way to set the device permittivity attribute is through materials.
See qtcad.device.materials
for a description of materials attributes,
which include permittivity.
Note
The default dielectric permittivity value is that of silicon,
Relevant solver attributes
The solver attributes are set by instantiating a
SolverParams
object with the desired attribute values and inputting it in the
solver_params
keyword argument of the solver constructor.
The default values of the relevant linear poisson solver attributes are given
in the SolverParams
API reference.
We remark that the same workflow is used to configure the parameters of all other QTCAD solvers.
Nonlinear Poisson solver
Differential equation
The nonlinear Poisson solver solves the partial differential equation
where \(\varepsilon\) is the dielectric permittivity of the medium and \(\varphi\) is the electric potential.
In the above equations, the total charge density \(\rho\) is decomposed into two contributions. The first contribution, \(\rho_\mathrm{semi}\), arises from charges in a semiconductor at equilibrium:
where \(e\) is the elementary charge (\(e>0\)), \(p\) and \(n\) are the hole and electron densities, respectively, and \(N_+\) and \(N_\) are the densities of ionized donors and acceptors, respectively.
In general, \(p\) and \(n\) depend, themselves, on the variable \(\varphi\) for which the equation is solved, making this Poisson equation nonlinear. In incomplete ionization models, \(N_+\) and \(N_\) may also depend on \(\varphi\).
The second contribution to the total charge in
Eq. (1.2), \(\rho_0\), is a fixed background volume
charge density that may be set by the user with the method
set_vol_charge_density
to model additional charges whose
density is independent of \(\varphi\). This may be used, for example, to
model charged defects.
Carrier statistics
In QTCAD, carrier statistics may be calculated using two main models: Fermi–Dirac statistics and Maxwell–Boltzmann statistics.
While electron and hole populations in semiconductors are fundamentally determined by Fermi–Dirac statistics, at room temperature and weak doping, Maxwell–Boltzmann statistics may provide a sufficiently good approximation.
The carrier statistics used when solving the nonlinear Poisson equation is
determined by the statistics
attribute of objects of the Device
class.
This attribute may have three possible values: "MB"
(Maxwell–Boltzmann),
"FD"
(Fermi–Dirac), or "FD_approx"
(approximate Fermi–Dirac). For
example, to set the value of this attribute to "MB"
for a hypothetical
Device
object called dvc
:
dvc.statistics = "MB"
The default value of the statistics
attribute is "FD_approx"
.
In a bulk semiconductor at equilibrium, assuming isotropic and parabolic bands, the electron and hole densities are given by [AM76]
where \(E_C\) and \(E_V\) are the conduction and valence band edges, respectively, \(E_F\) is the Fermi energy, \(T\) is the device temperature (assumed to be uniform), \(k_B\) is the Boltzmann constant, and \(F_{\frac{1}{2}}(x)\) is the complete Fermi–Dirac integral of order \(\frac{1}{2}\), which is explicitly given by
Finally, in Eq. (1.3), we have introduced the electron and hole effective densities of states
where \(g_{c(v)}\) is the total electron (hole) band degeneracy (including both spin and valley degeneracies), \(\hbar\) is the reduced Planck constant, and \(m_{c(v)}\) is the electron (hole) densityofstates effective mass.
Note
Eq. (1.3) may be applied to conduction electrons in silicon even if its bands are anisotropic. This can be done simply by taking \(g_c = 2\times 6 =12\) (2 for spin degeneracy, 6 for valley degeneracy) and \(m_c = (m_\ell m_t^2)^{1/3}\), where \(m_\ell\) and \(m_t\) are the longitudinal and transverse effective masses, respectively.
Note
In QTCAD finiteelement simulations, equilibrium is assumed throughout the
simulation domain by default. As a consequence, the Fermi level (chemical
potential) is typically uniform. We further choose the arbitrary zero of
our energy scale so that in the default behavior, the Fermi level is zero
throughout the device, \(E_F\equiv 0\). However, the Fermi energy
may be modified manually to depart from this default equilibrium behavior
using the device attribute setter method
set_fermi_lvl
.
Though arbitrary functions of space may be used in principle, in practice,
it is often useful to modify the Fermi energy only in specific regions
forming the source and the drain. Though phenomenological, this approach
is appropriate when the source and drain are individually at
quasiequilibrium, and only weakly interacting with a quantum dot system
through tunnel couplings.
We now distinguish the three possible choices for the statistics
device
attribute.
When statistics
is chosen to be "FD"
, the Fermi–Dirac integral,
Eq. (1.4), is calculated directly, numerically. This
is the most accurate approach; however, it imposes a computational overhead
that can significantly slow down calculations.
When statistics
is chosen to be "FD_approx"
(which it is, by default),
the Fermi–Dirac integral, Eq. (1.4), is approximated
analytically by [BB78]
The maximum relative error introduced by Eq. (1.6) over all possible values of \(x\) is 0.4 %, which is negligible in most practical cases.
Finally, when statistics
is chosen to be "MB"
, the Fermi–Dirac integral
is approximated by
Relationship between electric potential and band edges
The electric potential is related to the band edges through [GNM+13]
where \(\chi\) and \(E_g\) are the electron affinity and band gap of the material. In addition, we have introduced \(\varphi_F\), a constant arbitrary reference potential. This potential does not affect the physics of the device, since a constant shift in electric potential corresponds to a gauge transformation that has no impact on electric fields. By default, \(e\varphi_F= E_w^\mathrm{Si}\), where \(E_w^\mathrm{Si}\) is the work function of intrinsic silicon at the device temperature. Changing the reference potential has no impact on the physics, but may affect the convergence of the solvers.
Note
If the structure under consideration contains several different materials (i.e. if it is a heterostructure), the parameters of Eq. (1.8) have a spatial dependence. It is this spatial dependence that determines band alignment across heterojunctions. The case of constant \(\varphi_F\) corresponds to Anderson’s rule for band alignment in heterostructures. For more details on band alignment, see Band alignment in heterostructures.
Dopant ionization models
QTCAD supports two ionization models: complete ionization and incomplete ionization.
The default ionization model in a device that is input to the nonlinear
Poisson solver is specified by setting its ionization
attribute through
the method
set_ionization_model
.
This attribute may have either of the following two values: "complete"
or
"incomplete"
. By default, ionization
is complete
.
For example, to set the ionization
attribute of a hypothetical
Device
object called dvc
to "incomplete"
:
dvc.set_ionization_model("incomplete")
Under complete ionization, all donors and acceptors are assumed to be ionized, meaning that \(N_+=N_D\) and \(N_=N_A\), where \(N_D\) and \(N_A\) are the donor and acceptor densities, respectively.
Under incomplete ionization, assuming thermal equilibrium of the dopants with the carriers, the ionized dopant concentrations are given by [Gru10]
where \(g_D\) (\(g_A\)) is the donor (acceptor) level degeneracy and \(E_d\) (\(E_a\)) is the donor (acceptor) binding energy.
The default ionization model may be overridden when defining regions
(see Regions and boundaries) using the optional argument
ionization_model
of the
new_region
method of the
Device
class. For example,
assuming that the default ionization model of a
Device
object called dvc
has been previously set to incomplete
, complete ionization
may be imposed in a region called, e.g., "source"
, using the command
dvc.new_region("source", mt.Si, ndoping=1e20*1e6,
ionization_model="complete")
where we assume that the materials
module
has already been imported as mt
. This is particularly useful at very
low temperatures at which weakly doped regions will follow the incomplete
ionization model, but at which regions doped at a concentration much higher
than the critical Mott dopant concentration may still be completely ionized.
Relevant device attributes
Parameter 
Symbol 
QTCAD name 
Unit 
Default 
Setter 
Permittivity 
\(\varepsilon\) 

F/m 
11.8 \(\varepsilon_0\) 

Electron effective mass 
\(m_c\) 

kg 
0.321 \(m_e\) 
Through 
Electron degeneracy 
\(g_c\) 

– 
12 
Through 
Hole effective mass 
\(m_v\) 

kg 
0.549 \(m_e\) 
Through 
Hole degeneracy 
\(g_v\) 

– 
2 
Through 
Electron affinity 
\(\chi\) 

J 
4.05 eV 
Through 
Band gap 
\(E_g\) 

J 
1.12 eV 
Through 
Donor density 
\(N_D\) 

\(\mathrm m^{3}\) 
0 

Donor binding energy 
\(E_d\) 

J 
46 meV 

Donor level degeneracy 
\(g_D\) 

– 
2 

Acceptor density 
\(N_A\) 

\(\mathrm m^{3}\) 
0 

Acceptor binding energy 
\(E_a\) 

J 
44 meV 

Acceptor level degeneracy 
\(g_A\) 

– 
4 

Total charge density 
\(\rho\) 

\(\mathrm{C/m}^{3}\) 
0 
Through dopant densities, background density, and solver iterations 
Background charge density 
\(\rho_0\) 

\(\mathrm{C/m}^{3}\) 
0 

Fermi energy 
\(E_F\) 

eV 
0 

Reference potential 
\(\varphi_F\) 

V 
See note 

Temperature 
\(T\) 

K 
300.0 
Note
The default device attributes are those of silicon.
The default dopant parameters are those of phosphorus (for donors) or boron (for acceptors) in silicon.
The default reference potential is the silicon work function at the device temperature.
More information about materials can be found in
materials
.
Special care must be taken when defining device parameters using setters. Device attribute setters must be called before setting boundary conditions, because those attributes may enter the formulas used to calculate boundary values.
Boundary Conditions
Here, we describe the boundary conditions that may be used in combination with both linear and nonlinear Poisson solvers.
Dirichlet boundaries (API reference
)
On Dirichlet boundaries, the value of the electric potential is simply set to be equal to the applied potential \(\varphi_\mathrm{bias}\), i.e.
Gate boundaries (API reference
)
Gate boundaries correspond to metallic surfaces in contact with an insulator. In such cases, alignment of the metal, oxide, and semiconductor vacuum levels leads to an electric potential given by [GNM+13]
where \(E_F\) is the Fermi energy at the boundary nodes, \(\varphi_F\) is the reference potential (see Relationship between electric potential and band edges), \(e\) the elementary charge, and \(E_w\) is the work function of the metallic gate.
Schottky boundaries (API reference
)
Schottky boundaries correspond to interfaces between a metallic gate and a semiconductor. At such boundaries, the electric potential is given by
where \(\Phi_B\) is the Schottky barrier height (in the absence of any applied potential) and \(\chi\) is the electron affinity of the semiconductor. The Schottky barrier height is defined as the difference between the interfacial conduction band edge \(E_C\) and the Fermi level \(E_F\).
Ohmic boundaries (API reference
)
Ohmic conditions may be applied to boundaries in contact with a semiconductor. At Ohmic boundaries, the semiconductor is assumed to be at a chargeneutral equilibrium with its environment, leading to the condition [GNM+13]:
where the carrier densities \(p\) and \(n\) are given by Eq. (1.3), while the formulas for the ionized dopant densities \(N_+\) and \(N_\) are given in Dopant ionization models.
Since these quantities depend on electric potential \(\varphi\), Eq. (1.13) results in a nonlinear equation that is solved for \(\varphi\), giving the value of the electric potential at Ohmic boundaries.
Note that in QTCAD, equilibrium is assumed throughout the device. As a consequence, no external potential may be applied at Ohmic boundaries, since this would violate the equilibrium condition.
Frozen boundaries (API reference
)
For semiconductors at very low temperatures, such that \(k_B T/e \ll E_{d(a)}\), where \(E_{d(a)}\) is the donor (acceptor) binding energy, and under the incomplete ionization model, the solution of (1.13) asymptotically tends to
In other words, the Fermi level lies between the dopant level and the conduction (valence) band edge. Using Relationship between electric potential and band edges, this leads to the following boundary condition for the electric potential:
where the band edges are given by Eq. (1.14).
In this low temperature regime and under the incomplete ionization model, we expect dopants to be completely frozen out, hence the name of the boundary condition.
Surface charges (API reference
)
An arbitary background surface charge density \(\sigma\) may be added to any surface of a device. The surface may be at a boundary, or it may separate two regions, in which case it may be called an interface.
Charged defects at boundaries or interfaces often play an important role in device physics.
At an interface between two regions 1 and 2, surface charges impose the following continuity relation between the displacement field \(\mathbf D\equiv\varepsilon\nabla\varphi\) in each region [Jac99]:
where \(\mathbf{\hat n}\) is the unit vector normal to the interface, pointing from region 1 to region 2, and \(\mathbf D_i\) is the displacement field in region \(i\), infinitesimally close to the interface.
At a boundary, the surface charge density imposes a Neumann boundary condition on the potential, i.e., it sets the potential gradient through:
where \(\mathbf{\hat n}\) is the unit vector normal to the boundary, pointing inside the simulation domain, and \(\mathbf D\) is the displacement field infinitesimally close to the boundary yet still inside the simulation domain. In the TCAD literature, such a boundary condition is sometimes called a floating metal gate boundary condition.
Natural boundaries
All boundaries at which no specific boundary condition is specified by the user are called natural boundaries. For Poisson’s equation (linear or nonlinear), such natural boundary conditions correspond to
where \(\mathbf{\hat n}\) is the unit vector normal to the surface, pointing outside the simulation domain. In other words, natural boundaries set the component of the displacement field which is normal to the surface to zero.
Point Charges
Another feature available for QTCAD’s linear and nonlinear Poisson solvers is
the possibility of adding point charges.
Point charges are added to any region of the device using the method
add_point_charges
.
With the inclusion of point charges, the Poisson equation becomes
subject to the essential boundary conditions
where the \(\Sigma_i\) are the boundary regions over which the potential is forced to be \(\varphi_i\) (according to the Dirichlet, Gate, Schottky, Ohmic, or Frozen boundaries defined above) and
is the charge density due to point charges. To avoid numerical instabilities associated with the \(\delta\)function, point charges are spread over a small region (typically on the order of \(0.1\) nm) around their position. In other words, the \(\delta\)functions of Eq. (1.21) are approximated by a Gaussian
or a uniform distribution
where \(\sigma\) is the radius of the ball over which the point charge is spread.
Nonlinear Poisson solver with point charges
Because point charges produce a diverging potential, \(\varphi\), they can lead to situations where the conductionband edge drops very far below the Fermi level or where the valenceband edge climbs very far above the Fermi level. Both situations would lead to an unphysical density of electrons (in the former case) or holes (in the later case) in gases accumulating near the point charges if the nonlinear Poisson solver is directly applied. In reality, we expect the point charges observed experimentally in semiconductor devices to be quantized and not associated to some continuum of charge (gas). Therefore we treat point charges separately when solving the nonlinear Poisson equation. Specifically, we perform two calculations. The first calculation is the standard nonlinear Poisson calculation without point charges, where we solve the following equation selfconsistently:
subject to the essential boundary conditions of Eq. (1.20)
Once that has converged, we perform the second calculation, a linear Poisson calculation with the point charges
subject to the modified essential boundary conditions
Once both these equations have been solved, we add the solutions together to obtain the final potential, \(\varphi = \varphi_0 + \varphi_{\delta}\). Inserting \(\varphi = \varphi_0 + \varphi_{\delta}\) into the Poisson equation, we can verify that \(\varphi\) is a solution to Eq. (1.19):
Moreover, the boundary conditions of Eq. (1.20) are also satisfied [see Eqs. (1.27)]:
By separating the calculation over two steps, we avoid the unphysical accumulation of carriers near the point charges since \(\rho_{\mathrm{semi}}\) is computed independent of the point charges.
Relevant device attributes
Parameter 
Symbol 
QTCAD name 
Unit 
Default 
Setter 
Point charge density 
\(\rho_{\delta}\) 

\(\mathrm{C m}^{3}\) 
0 

Point charge radius 
\(\sigma\) 

\(\mathrm m\) 
– 

Point charge positions 
\(\mathbf{r}_j\) 

\(\mathrm m\) 
– 

Point charge values 
\(Q_j\) 

\(\mathrm C\) 
– 

Potential due to point charges 
\(\varphi_{\delta}\) 

\(\mathrm V\) 
0 
– 