7. Schrödinger–Poisson solver
In previous sections, we explained the theories underdinning QTCAD’s nonlinear Poisson and Schrödinger solvers separately. In this approach, the confinement potential is first found within the Thomas–Fermi approximation by solving the nonlinear Poisson equation; Schrödinger’s equation may then be solved to find the electronic structure of singleelectron or singlehole bound states in this potential (if any).
Notably, this approach neglects the backaction of the charges that are confined in the nanostructure on the electric potential. For fewelectron and fewhole systems in quantum dots, this backaction is negligible in many cases. However, in some circumstances, quantummechanically confined charges may significantly alter the band profile of the structure. This is particularly important for systems in which quantum confinement only exists along certain directions; e.g. in quantum wells (confinement along one direction, translational invariance along two directions) or in nanowires (confinement along two directions, translational invariance along one direction). In such cases, quantum confinement is only partial, and many carriers can occupy the quantized levels (called subbands), forming 1D or 2D electron or hole gases that can significantly influence band bending. These effects may be accounted for with selfconsistent Schrödinger–Poisson simulations.
Differential equations
In Schrödinger–Poisson simulations, we solve the nonlinear Poisson and Schrödinger equations selfconsistently.
In the first step of this iterative procedure, the electric potential \(\varphi\) (and the resulting conduction and valence band edges, see Relationship between electric potential and band edges and Band alignment in heterostructures) is called the “initial guess” and is typically taken to be the solution of the nonlinear Poisson equation for a classical, 3D electron (or hole) gas.
A selfconsistent loop is then entered in which inner (nonlinear Poisson) and outer (Schrödinger) iteration schemes are implemented:
In each outer iteration step, Schrödinger’s equation is solved using the electric potential \(\varphi_{\mathrm {prev.}}\) from the previous outer iteration step.
After each outer iteration, an inner loop is entered in which the nonlinear Poisson’s equation is solved to find the next electric potential \(\varphi_{\mathrm {next}}\), and in which the statistical model giving the charge density takes quantum confinement into account.
Inner and outer iteration schemes are alternated until variations of the electric potential between outer iterations are below a certain tolerance, \(\varphi < \varphi_{\mathrm{tol}}\). At that point, convergence (or selfconsistency) is reached.
Schrödinger equations
For quantummechanically confined electrons, we consider the following Schrödinger equation:
For holes, we instead use:
Confinement potential
The confinement potential in (7.1) and (7.2) is given by the equations introduced in Band alignment in heterostructures, to which we may optionally add an exchange–correlation term \(V_\mathrm{xc}(\mathbf r)\) derived from densityfunctional theory. More precisely, the total confinement potential is given by
where \(V(\mathbf r)\) and \(V_\mathrm{ext}(\mathbf r)\) have the same meaning as in Band alignment in heterostructures. When activated, in QTCAD, the exchange–correlation term is obtained from the Hedin–Lundqvist parameterization of the local density approximation (LDA) functional, i.e. (for electrons) [GNM+13, HL71, TR04]
where \(m^\ast_\mathrm{xc}\) is a scalar effective mass which is specific to this parameterization. As suggested in [TR04], when the inverse effective mass tensor is anisotropic, an average effective mass may be used.
Note
QTCAD currently does not support exchange–correlation for holes. When holes are used as the confined carriers, the exchange–correlation potential is set to zero automatically.
Poisson’s equation
We now remark that the term \(V(\mathbf r)\) in the total confinement potential given by (7.3) depends on the electrostatic potential \(\varphi\) (see Band alignment in heterostructures) which is obtained from the solution of the nonlinear Poisson equation
Carrier densities
When the confined carriers are chosen to be electrons (by setting the
conf_carriers
of the
Device
object to "e"
), the hole
density \(p\) is given by its classical expression,
(1.3), and vice versa.
For electrons in three dimensions, the confined carrier density is calculated quantummechanically from
where \(n_F(x)\equiv 1/(1+\mathrm e^x)\) is the Fermi–Dirac distribution, \(E_i\) and \(F_i(\mathbf r)\) are the \(i\)th quantized energy level and envelope function of the nanostructure, respectively, and \(g_{q,e}\) is the total (spin and valley) degeneracy of the quantized energy levels for electrons.
For holes in three dimensions, the confined carrier density is given by
where \(F_{n,i}(\mathbf r)\) is the projection on band \(n\) of the envelope function of the \(i\)th eigensolution of the coupled envelope function equations with energy \(E_i\). Note that unlike Eq. (7.6), Eq. (7.7) does not have a degeneracy factor; possible degeneracies are instead considered through the sum over \(i\).
Note that QTCAD enforces the following normalization conditions for electrons and holes:
Carrier densities in translationallyinvariant systems
Certain structures display translational invariance along one or two direction(s) and quantum confinement along the other(s). We will note positions within coordinate systems with translational invariance using \(\mathbf r_\parallel\), and those without invariance with \(\mathbf r_\perp\). For example, a nanowire has translational invariance along one direction (e.g. \(\mathbf{\hat z}\)), and confinement along two directions (e.g. \(\mathbf{\hat x}\) and \(\mathbf{\hat y}\)). By contrast, a quantum well has translational invariance along two directions (e.g. \(\mathbf{\hat y}\) and \(\mathbf{\hat z}\)), and confinement along one direction (e.g. \(\mathbf{\hat x}\)).
For structures with translational invariance, the envelope functions take the form
where \(\mathbf k_\parallel\) is a wavevector defined in the translationinvariant direction(s). To calculate the confined carrier densities, it is then necessary to integrate over all possible values of this wavevector. For electrons with quadratic bands, assuming an \((x,y,z)`\) coordinate system that coincides with the principal axes of the effective mass tensor (i.e. in which the inverse effective mass tensor is diagonal), this yields new analytic expressions for the quantummechanically confined electron density, which is now given by [GNM+13]
where \(\eta_i\equiv(E_FE_i)/k_B T\), and the expression for \(\mathcal N\) depends on dimensionality.
For 1D confinement (i.e. for two translationinvariant directions), we have [GNM+13]
where \(m_\perp^\ast \equiv \sqrt{m_y m_z}\) is an average effective mass in the translationinvariant directions.
For 2D confinement (i.e. for one translationinvariant direction), we rather have [GNM+13]
where \(m_\perp^\ast \equiv m_z\) is the average effective mass along the translationinvariant direction.
Finally, in (7.11) and (7.12), we have introduced the complete Fermi–Dirac integral of order \(q\), defined by
where \(\Gamma(x)\) is the Gamma function. Note that \(\mathcal F_0(\eta)=\ln[1+\exp(\eta)]\).
For holes, the calculation of the confined carrier density under translation invariance is complicated by band mixing, under which the analytic expressions given above for electrons may not be valid. In this case, the integral over \(\mathbf k_\parallel\) should be done numerically for each subband, which is a computationally intensive overhead. Consequently, when solving Schrödinger–Poisson for holes in 1D, bandmixing terms are neglected (band offdiagonal matrix elements of the Hamiltonian are set to \(0\)). The leading corrections to the carrier densities and eigenenergies are \({\cal{O}}(\epsilon^2)\),
where the proportionality constant, \(\alpha\), depends on material parameters [Luttinger parameters see Eqs. (3.4), (3.5), and (3.6)] and physical constants, the maximum is taken over all occupied states, and \(\Delta\) represents the heavyhole lighthole splitting. Thus, this approach is welljustified the limit of tight 1D quantum confinement (large enough band splittings) and small carrier density (the maximal value of \(\mathrm{max}\left(k_{\parallel}\right)\) over all occupied states must be small enough).
In the absence of band mixing, each state \(i\) will have a unique “band character”, e.g. \(n=1/2\) for the (lighthole) state that transforms like a state with \(J_z=1/2\) units of angular momentum along the axis of quantization [as a reminder, for holes, \(n\) runs over the two heavyhole and the two lighthole states, see Eq. (3.15)]. We define \(n_i\) to be the character of eigenstate \(i\). The carrier density is given by
where
where \(m_{\perp, n}^\ast \equiv \sqrt{m_{y, n} m_{z, n}}\) is an average effective mass of a band of character \(n\) in the translationinvariant directions.
Note
2D Schrödinger–Poisson calculations for holes are not yet supported.
Note
For 1D and 2D simulations, it is important to respect the following
conventions regarding the geometry. For systems with translational invariance
along a single direction, this direction is the \(z\) axis, by convention,
and a 2D geometry must be drawn in Gmsh or devicegen in the \((x,y)\)
plane. For systems with translational invariance along two directions, these
directions are the \((y,z)\) plane, by default, and a 1D geometry must
be drawn in Gmsh along the \(x\) axis.
While the conventions for the geometries defined using Gmsh must always be
respected, the translationallyinvariant directions can be modified via the
ti_directions
attribute of the Schrödinger
solver
.
The solvers will reinterpret the directions defined in Gmsh to respect the
choice for translationallyinvariant directions. For example, if a 1D
geometry is drawn along the \(x\) axis in Gmsh, imported into QTCAD,
and
ti_directions
is set to ["x", "y"]
, the QTCAD mesh will be rotated so that the
\(x\) axis in Gmsh corresponds to the \(z\) axis in QTCAD.
Carrier densities in nearlytranslationallyinvariant systems
There is a separate class of systems which is characterized by quantum confinement along one or two dimensions without pure translational invariance along the other(s). Although these systems display no pure translational invariance, we will assume that they are nearly translationally invariant. This concept will be made clearer with the example below.
One type of system which can display near translational invariance is the quantum well. Quantum wells are typically formed near heterostructure interfaces where the difference in band gaps and band offsets creates a confinement potential for electrons and/or holes along the growth axis. In this basic form, the quantum well possesses translational invariance along the directions orthogonal to the growth direction. However, the twodimensional particle (electron or hole) gases confined to quantum wells are often manipulated via electric potentials applied to gates to generate additional confinement and create systems like gated quantum dots or quantum point contacts. These additional electric potentials break the pure translational symmetry of the system, leaving behind systems with 3D confinement (quantum dots) and 2D confinement (quantum point contacts) and a surrounding particle gas which we will assume is only locally translationally invariant. To analyze these systems we employ the quantumwell solver.
The quantumwell solver can only be applied to a 3D
device
(device defined over a 3D
mesh
).
The user must also specify a set of coordinates, \(\{(x_j, y_j)\}\), which
are used to generate related 1D systems, indexed by \(j\).
The potential for system \(j\) is obtained by taking a linecut of the full
3D potential along the confinement direction, \(z\), for fixed values
\(x_j\) and \(y_j\).
The 1D Schrödinger equation is then solved for each of these 1D systems to
obtain eigenstates \(F^j_{\perp,i}(z)\) and energies \(E_i^j\).
For electrons, the carrier densities along the different linecuts, \(j\),
are obtained via Eqs. (7.10) and
(7.11)
where \(\eta_i\equiv(E_FE^j_i)/k_B T\). For holes, Eqs. (7.15) and (7.16) are used instead,
These carrier densities are then interpolated throughout the 3D device to obtain the full 3D carrier densities. We have therefore implicitly assumed that the system has neartranslational invariance in the vicinity of points \(\{(x_j, y_j)\}\) and that the potential along the \(x\) and \(y\) directions varies slow enough so that the interpolation of the carrier density between points \(\{(x_j, y_j)\}\) is accurate.
Note
Similar to other QTCAD solvers, the quantumwell
solver
will set the charge to 0 in
dot regions
.
Once the quantumwell solver
has
converged, a subdevice
composed of
the dot region can be passed to the Schrödinger
solver
to get the quantumdot
singleparticle eigenstates and eigenenergies.
These singleparticle eigenstates and eigenenergies can then be passed to the
manybody solver
to get the
manybody spectrum of the dot.
Therefore, the quantumwell solver
can be used to understand e.g. how the electrons/holes confined to a gated
quantum dot (3D confinement) are affected by the surrounding 2D electron/hole
gas confined to the quantum well (1D confinement).
Note
As is the case for holes in translationallyinvariant systems, when using the
quantumwell solver
to study holes
confined in only 1 direction, bandmixing terms are neglected (band
offdiagonal matrix elements of the Hamiltonian are set to \(0\)).
Boundary conditions
The boundary conditions in Schrödinger–Poisson calculations are those described in the Schrödinger and Poisson theory sections for the envelope functions and electric potential, respectively, without modification.
Relevant device attributes
Here, we list Device
attributes that
are uniquely used by the Schrödinger–Poisson solver. For attributes that
are related to individual Schrödinger
and Poisson solvers, see the tables
in the appropriate theory sections.
Parameter 
Symbol 
QTCAD name 
Unit 
Default 
Setter 
Total degeneracy of confined electrons 
\(g_{q,e}\) 

– 
4.0 
Through 
Exchange–correlation effective mass 
\(m_\mathrm{xc}\) 

kg 
0.258\(m_e\) 
Through 