7. Schrödinger–Poisson solver

In previous sections, we explained the theories underdinning QTCAD’s non-linear Poisson and Schrödinger solvers separately. In this approach, the confinement potential is first found within the Thomas–Fermi approximation by solving the non-linear Poisson equation; Schrödinger’s equation may then be solved to find the electronic structure of single-electron or single-hole bound states in this potential (if any).

Notably, this approach neglects the back-action of the charges that are confined in the nanostructure on the electric potential. For few-electron and few-hole systems in quantum dots, this back-action is negligible in many cases. However, in some circumstances, quantum-mechanically 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 self-consistent Schrödinger–Poisson simulations.

Differential equations

In Schrödinger–Poisson simulations, we solve the non-linear Poisson and Schrödinger equations self-consistently.

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 non-linear Poisson equation for a classical, 3D electron (or hole) gas.

A self-consistent loop is then entered in which inner (non-linear 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 non-linear 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 self-consistency) is reached.

Schrödinger equations

For quantum-mechanically confined electrons, we consider the following Schrödinger equation:

(7.1)\[\begin{split}H(\mathbf r) F(\mathbf r) &= E F(\mathbf r),\\ H(\mathbf r) F(\mathbf r) &= V_\mathrm{conf}(\mathbf r)F(\mathbf r) - \frac{\hbar^2}2\nabla\cdot \mathbf M^{-1}_e \cdot \nabla F(\mathbf r).\end{split}\]

For holes, we instead use:

(7.2)\[\sum_{m} \left[V_\mathrm{conf}(\mathbf r) \delta_{nm} + \sum_{\alpha\beta} k_\alpha D^{\alpha\beta}_{nm} k_\beta\right] F_m(\mathbf r) = E F_n(\mathbf r).\]

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 density-functional theory. More precisely, the total confinement potential is given by

(7.3)\[V_\mathrm{conf}^\mathrm{xc}(\mathbf r) = V(\mathbf r) + V_\mathrm{ext}(\mathbf r) + V_\mathrm{xc}[n],\]

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]

(7.4)\[\begin{split}V_\mathrm{xc}[n]&=-\frac{e^2}{4\pi^2\varepsilon}[3\pi^2n(\mathbf r)]^{\frac13} \left[1+0.7734 x\ln\left(1+\frac1x\right)\right],\\ x &= \frac1{21}\left(\frac{4\pi n(\mathbf r)b^3}{3}\right)^{-\frac13},\\ b &= \frac{4\pi\varepsilon\hbar^2}{m^\ast_\mathrm{xc} e^2},\end{split}\]

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

(7.5)\[\nabla\cdot\left(\varepsilon\nabla\varphi\right)=-e(p-n+N_+-N_-).\]

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 quantum-mechanically from

(7.6)\[n(\mathbf r) = g_{q,e} \sum_i n_F\left(\frac{E_i-E_F}{k_B T}\right) |F_i(\mathbf r)|^2,\]

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

(7.7)\[p(\mathbf r) = \sum_i n_F\left(\frac{E_i-E_F}{k_B T}\right) \sum_n|F_{n,i}(\mathbf r)|^2,\]

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:

(7.8)\[\begin{split}\int d\mathbf r |F_i(\mathbf r)|^2 &= 1,\qquad\mbox{(electrons)}\\ \sum_n\int d\mathbf r |F_{n,i}(\mathbf r)|^2 &= 1. \qquad\mbox{(holes)}\end{split}\]

Carrier densities in translationally-invariant 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

(7.9)\[F(\mathbf r) = \mathrm e^{i\mathbf k_\parallel\cdot \mathbf r_\parallel} F_\perp(\mathbf r_\perp),\]

where \(\mathbf k_\parallel\) is a wavevector defined in the translation-invariant 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 quantum-mechanically confined electron density, which is now given by [GNM+13]

(7.10)\[n(\mathbf r_\perp) = \sum_i \mathcal N\left(\eta_i\right) |F_{\perp,i}(\mathbf r_\perp)|^2,\]

where \(\eta_i\equiv(E_F-E_i)/k_B T\), and the expression for \(\mathcal N\) depends on dimensionality.

For 1D confinement (i.e. for two translation-invariant directions), we have [GNM+13]

(7.11)\[\mathcal N\left(\eta\right) = \frac{g_{q,e}}{2}\frac{m_\perp^\ast k_B T}{\pi\hbar^2} \mathcal F_0(\eta_i),\]

where \(m_\perp^\ast \equiv \sqrt{m_y m_z}\) is an average effective mass in the translation-invariant directions.

For 2D confinement (i.e. for one translation-invariant direction), we rather have [GNM+13]

(7.12)\[\mathcal N\left(\eta\right) = \frac{g_{q,e}}{2} \left(\frac{2 m_\perp^\ast k_B T}{\pi\hbar^2}\right)^{\frac12} \mathcal F_{-\frac12}(\eta),\]

where \(m_\perp^\ast \equiv m_z\) is the average effective mass along the translation-invariant direction.

Finally, in (7.11) and (7.12), we have introduced the complete Fermi–Dirac integral of order \(q\), defined by

(7.13)\[\mathcal F_q(x)= \frac1{\Gamma(q+1)} \int_0^\infty dt \frac{t^q}{\mathrm e^{t-x}+1},\]

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, band-mixing terms are neglected (band off-diagonal matrix elements of the Hamiltonian are set to \(0\)). The leading corrections to the carrier densities and eigenenergies are \({\cal{O}}(\epsilon^2)\),

(7.14)\[\epsilon \lesssim \alpha \frac{\mathrm{max}\left(k_{\parallel}\right)}{\Delta},\]

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 heavy-hole light-hole splitting. Thus, this approach is well-justified 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 (light-hole) 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 heavy-hole and the two light-hole states, see Eq. (3.15)]. We define \(n_i\) to be the character of eigenstate \(i\). The carrier density is given by

(7.15)\[p(\mathbf r_\perp) = \sum_{i} \mathcal N_i\left(\eta_i\right) |F_{\perp,i,n_i}(\mathbf r_\perp)|^2,\]

where

(7.16)\[\mathcal N_i\left(\eta\right) = \frac{m_{\perp, n_i}^\ast k_B T}{2\pi\hbar^2} \mathcal F_0(\eta),\]

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 translation-invariant 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 translationally-invariant 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 translationally-invariant 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 nearly-translationally-invariant 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 two-dimensional 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 quantum-well solver.

The quantum-well 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)

(7.17)\[n(x_j, y_j, z) = \sum_i \mathcal N\left(\eta^j_i\right) |F^j_{\perp,i}(z)|^2,\]

where \(\eta_i\equiv(E_F-E^j_i)/k_B T\). For holes, Eqs. (7.15) and (7.16) are used instead,

(7.18)\[p(x_j, y_j, z) = \sum_{i} \mathcal N_i\left(\eta^j_i\right) |F^j_{\perp,i,n_i}(z)|^2.\]

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 near-translational 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 quantum-well solver will set the charge to 0 in dot regions. Once the quantum-well solver has converged, a subdevice composed of the dot region can be passed to the Schrödinger solver to get the quantum-dot single-particle eigenstates and eigenenergies. These single-particle eigenstates and eigenenergies can then be passed to the many-body solver to get the many-body spectrum of the dot. Therefore, the quantum-well 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 translationally-invariant systems, when using the quantum-well solver to study holes confined in only 1 direction, band-mixing terms are neglected (band off-diagonal 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}\)	gqc	–	4.0	Through materials
Exchange–correlation effective mass	\(m_\mathrm{xc}\)	mc_xc	kg	0.258\(m_e\)	Through materials