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 non-linear Poisson solver which self-consistently accounts for the carrier densities in classical reservoirs.

Linear Poisson solver

The linear Poisson solver solves the partial differential equation

(1.1)\[-\nabla\cdot\left(\varepsilon\nabla\varphi\right)=\rho,\]

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\)	eps	F/m	11.8 \(\varepsilon_0\)	set_permittivity
Charge density	\(\rho\)	rho	C/\(\mathrm m^3\)	0	set_charge_density

Variable	Symbol	QTCAD name	Unit
Potential	\(\varphi\)	phi	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.

Non-linear Poisson solver

Differential equation

The non-linear Poisson solver solves the partial differential equation

(1.2)\[ \begin{align}\begin{aligned}&-\nabla\cdot\left(\varepsilon\nabla\varphi\right) = \rho,\\&\rho = \rho_\mathrm{semi}+\rho_0.\end{aligned}\end{align} \]

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:

\[\rho_\mathrm{semi} = e(p-n+N_+-N_-),\]

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 non-linear. 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 non-linear 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]

(1.3)\[n = N_c F_{\frac{1}{2}}\left(\frac{E_F-E_C}{k_B T}\right),\qquad p = N_v F_{\frac{1}{2}}\left(\frac{E_V-E_F}{k_B T}\right),\]

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

(1.4)\[F_{\frac{1}{2}}(x) = \frac{2}{\sqrt \pi}\int_0^\infty dt \frac{\sqrt t}{1+\mathrm{e}^{t-x}}.\]

Finally, in Eq. (1.3), we have introduced the electron and hole effective densities of states

(1.5)\[N_c = \frac{g_c}8\left(\frac{2k_BTm_c}{\pi\hbar^2}\right)^{3/2},\qquad N_v = \frac{g_v}8\left(\frac{2k_BTm_v}{\pi\hbar^2}\right)^{3/2},\]

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) density-of-states 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 finite-element 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 quasi-equilibrium, 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]

(1.6)\[\begin{split}F_{\frac{1}{2}}(x)&\approx \left(\mathrm e^{-x}+\frac{3\sqrt\pi}4 v^{-3/8}\right)^{-1},\\ \textrm{where}\quad v &= x^4+50+33.6x\left\{1-0.68\exp\left[-0.17(x+1)^2\right]\right\}.\end{split}\]

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

(1.7)\[F_{\frac{1}{2}}(x) \approx \mathrm e^x.\]

Relationship between electric potential and band edges

The electric potential is related to the band edges through [GNM+13]

(1.8)\[\begin{split}E_C &= -e(\varphi-\varphi_F)-\chi,\\ E_V &= -e(\varphi-\varphi_F)-\chi-E_g,\end{split}\]

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 non-linear 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]

(1.9)\[\begin{split}N_+ &= \frac{N_D}{1+g_D\exp[(E_F-E_C+E_d)/k_BT]},\\ N_- &= \frac{N_A}{1+g_A\exp[(E_V-E_F+E_a)/k_BT]},\end{split}\]

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\)	eps	F/m	11.8 \(\varepsilon_0\)	set_permittivity
Electron effective mass	\(m_c\)	mc	kg	0.321 \(m_e\)	Through materials
Electron degeneracy	\(g_c\)	gc	–	12	Through materials
Hole effective mass	\(m_v\)	mv	kg	0.549 \(m_e\)	Through materials
Hole degeneracy	\(g_v\)	gv	–	2	Through materials
Electron affinity	\(\chi\)	chi	J	4.05 eV	Through materials
Band gap	\(E_g\)	Eg	J	1.12 eV	Through materials
Donor density	\(N_D\)	ND	\(\mathrm m^{-3}\)	0	new_region or set_donor_density
Donor binding energy	\(E_d\)	Ed	J	46 meV	set_donor_binding_nrg
Donor level degeneracy	\(g_D\)	gD	–	2	set_donor_degen
Acceptor density	\(N_A\)	NA	\(\mathrm m^{-3}\)	0	new_region or set_acceptor_density
Acceptor binding energy	\(E_a\)	Ea	J	44 meV	set_acceptor_binding_nrg
Acceptor level degeneracy	\(g_A\)	gA	–	4	set_acceptor_degen
Total charge density	\(\rho\)	rho	\(\mathrm{C/m}^{-3}\)	0	Through dopant densities, background density, and solver iterations
Background charge density	\(\rho_0\)	rho_0	\(\mathrm{C/m}^{-3}\)	0	set_vol_charge_density
Fermi energy	\(E_F\)	EF	eV	0	set_fermi_lvl
Reference potential	\(\varphi_F\)	phi_F	V	See note	align_bands
Temperature	\(T\)	T	K	300.0	set_temperature

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 non-linear 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.

(1.10)\[\varphi = \varphi_\mathrm{bias}.\]

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]

(1.11)\[\varphi = \varphi_\mathrm{bias} - E_F/e + \varphi_F-E_w/e,\]

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

(1.12)\[\varphi = \varphi_\mathrm{bias} - E_F/e + \varphi_F -(\Phi_B+\chi)/e,\]

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 charge-neutral equilibrium with its environment, leading to the condition [GNM+13]:

(1.13)\[p - n + N_+ - N_- = 0,\]

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 non-linear 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

(1.14)\[\begin{split}E_C - E_F &= \frac{E_d}2 + \frac{k_B T}2\log\left( \frac{g_D N_c}{N_D}\right) \qquad &\mbox{(n doping)},\\ E_V - E_F &= -\frac{E_a}2 - \frac{k_B T}2\log\left( \frac{g_A N_v}{N_A}\right) \qquad &\mbox{(p doping)}.\end{split}\]

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:

(1.15)\[\begin{split}\varphi &= \varphi_\mathrm{bias}+\varphi_F-(\chi+E_C)/e \qquad &\mbox{(n doping)},\\ \varphi &= \varphi_\mathrm{bias}+\varphi_F-(\chi+E_V+E_g)/e \qquad &\mbox{(p doping)},\end{split}\]

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]:

(1.16)\[\left(\mathbf D_2-\mathbf D_1\right)\cdot \mathbf{\hat n}=\sigma,\]

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:

(1.17)\[\mathbf{\hat n}\cdot \mathbf D =\sigma,\]

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 non-linear), such natural boundary conditions correspond to

(1.18)\[\mathbf{\hat n}\cdot \left(\varepsilon \nabla \varphi\right) =0,\]

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 non-linear 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

(1.19)\[-\nabla\cdot\left(\varepsilon\nabla\varphi\right) = \rho + \rho_{\delta},\]

subject to the essential boundary conditions

(1.20)\[\varphi(\mathbf{r}) = \varphi_i, \quad \forall \mathbf{r} \in \Sigma_i.\]

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

(1.21)\[\rho_{\delta} (\mathbf{r}) = \sum_j Q_j \delta(\mathbf{r} - \mathbf{r}_j)\]

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

(1.22)\[\delta (\mathbf{r}) \approx \frac{1}{(2\pi)^{3/2}\sigma^3} e^{-\frac{r^2}{2\sigma^2}}\]

or a uniform distribution

(1.23)\[\begin{split}\delta (\mathbf{r}) &\approx\left\{ \begin{array}{ll} \frac{3}{4\pi\sigma^3} &\textrm{if}\;r<\sigma,\\ 0 &\textrm{otherwise}. \end{array} \right.\end{split}\]

where \(\sigma\) is the radius of the ball over which the point charge is spread.

Non-linear Poisson solver with point charges

Because point charges produce a diverging potential, \(\varphi\), they can lead to situations where the conduction-band edge drops very far below the Fermi level or where the valence-band 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 non-linear 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 non-linear Poisson equation. Specifically, we perform two calculations. The first calculation is the standard non-linear Poisson calculation without point charges, where we solve the following equation self-consistently:

(1.24)\[-\nabla\cdot\left(\varepsilon\nabla\varphi_{0}\right) = \rho_{\mathrm{semi}} + \rho_0,\]

subject to the essential boundary conditions of Eq. (1.20)

(1.25)\[\varphi_0(\mathbf{r}) = \varphi_i, \quad \forall \mathbf{r} \in \Sigma_i.\]

Once that has converged, we perform the second calculation, a linear Poisson calculation with the point charges

(1.26)\[-\nabla\cdot\left(\varepsilon\nabla\varphi_{\delta}\right) = \rho_{\delta}\]

subject to the modified essential boundary conditions

(1.27)\[\varphi_{\delta}(\mathbf{r}) = 0, \quad \forall \mathbf{r} \in \Sigma_i.\]

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):

(1.28)\[ \begin{align}\begin{aligned}-\nabla\cdot\left(\varepsilon\nabla\varphi\right) &= -\nabla\cdot\left(\varepsilon\nabla\left[\varphi_0 + \varphi_{\delta}\right]\right)\\&= -\nabla\cdot\left(\varepsilon\nabla\varphi_0\right) -\nabla\cdot\left(\varepsilon\nabla\varphi_{\delta}\right)\\&= \rho_{\mathrm{semi}} + \rho_0 + \rho_{\delta}.\end{aligned}\end{align} \]

Moreover, the boundary conditions of Eq. (1.20) are also satisfied [see Eqs. (1.27)]:

(1.29)\[\varphi(\mathbf{r}) = \varphi_0(\mathbf{r}) + \varphi_{\delta}(\mathbf{r}) = \varphi_i + 0, \quad \forall \mathbf{r} \in \Sigma_i.\]

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}\)	rho_pc	\(\mathrm{C m}^{-3}\)	0	add_point_charges
Point charge radius	\(\sigma\)	point_charge_radius	\(\mathrm m\)	–	add_point_charges
Point charge positions	\(\mathbf{r}_j\)	point_charge_pos	\(\mathrm m\)	–	add_point_charges
Point charge values	\(Q_j\)	point_charge_val	\(\mathrm C\)	–	add_point_charges
Potential due to point charges	\(\varphi_{\delta}\)	phi_pc	\(\mathrm V\)	0	–