5. Tight-binding Schrödinger solver

The tight-binding (TB) model is an atomistic approach to describe the electronic structure of atomic systems. In this model, quantum states and the Hamiltonian operator are expressed in a basis of orbitals centered on the atoms of the system. The TB Hamiltonian captures the energy associated with an electron lying on one of these orbitals, as well as the energy associated with the hopping of an electron between orbitals centered around neighboring atoms.

While the TB model was originally developed to describe the bandstructure of crystals [SK54], it has since been applied to describe the electronic properties of large, non-periodic systems such as quantum dots [LJonssonW+01] and nanoscale transistors [PHLG20]. Thanks to its explicit modeling of all individual atoms of a system, the TB model can capture intrinsically atomistic effects such as the confinement of electrons over a small number of atoms, atomistic disorder arising from random alloying and rough heterointerfaces, and strain fluctuations over nanoscale distances; such effects are difficult to properly capture in continuum models like the effective mass and \(\mathbf{k}\cdot\mathbf{p}\) models. In addition, thanks to its local modeling of electrons, Hamiltonian matrices that arise in the TB model tend to be very sparse compared to those arising in first-principles atomistic models like density functional theory. This allows for efficient numerical modeling of systems with up to millions of atoms, provided that TB model parameters are available.

The Solver class of the qtcad.atoms.schrodinger module implements a TB Schrödinger equation solver to resolve the electronic structure of quantum dots. It is a flexible solver which may take into account a variety of relevant effects, such as spin–orbit coupling, external electric potentials (including those computed via a finite-element Poisson solver, such as the nonlinear Poisson solver and the quantum-well solver, thereby enabling multiscale simulations), and magnetic field effects (including the Zeeman effect and magnetic orbital effects). The TB basis set used in this solver can be chosen from three different options with varying degrees of accuracy and computational cost.

The constructor method of the Solver class takes as input atoms, an Atoms object, which contains the atomic structure of the system. The solve method solves the TB Schrödinger equation by:

constructing the TB Hamiltonian matrix \(H\) of the system described by atoms;
performing a partial diagonalization of \(H\) to obtain a few eigenpairs of interest through an efficient, sparse eigensolver.

Here, we describe how the TB Hamiltonian matrix is constructed in QTCAD. In Basis sets, we describe which TB basis sets are available. In Block structure of the Hamiltonian, we describe the structure of the TB Hamiltonian; it may be partitioned into so-called onsite and offsite blocks, which we respectively describe in Onsite blocks and Offsite blocks. In Spin–orbit coupling, we describe how to include spin–orbit coupling in the TB Hamiltonian. In Treatment of surface atoms, we describe an effective boundary condition on the surface atoms of an atomic structure to avoid midgap surface states. In External electric potential and multiscale simulations, we describe how to include external electric potentials in the TB Hamiltonian, which are needed to perform multiscale simulations combining a finite-element Poisson solver with the TB model. Finally, in Magnetic effects, we describe how to include magnetic effects in the TB Hamiltonian, including the Zeeman effect and magnetic orbital effects.

Basis sets

The TB model basis set consists of atomic orbitals centered on the atoms of the system. A basis state may thus be denoted by the ket \(|n,i\rangle\), where \(n\) indexes the atom and \(i\) indexes the type of atomic orbital.

While each atom has, in principle, an infinite number of atomic orbitals, only a finite number of them are considered in the TB model. Generally, only valence orbitals (namely those in the outermost shell of the atom, energetically close to the Fermi energy) are included in the model, as they are sufficient to describe the relevant chemistry and physics of interatomic orbital hybridization. Three sets of atomic orbital types are available in QTCAD:

the \(sp^3\) atomic orbitals set, which consists of an \(s\) orbital, a \(p_x\) orbital, a \(p_y\) orbital, and a \(p_z\) orbital;
the \(sp^3s^{\star}\) atomic orbitals set, which additionally includes an \(s^{\star}\) orbital;
the \(sp^3d^5s^{\star}\) atomic orbitals set, which additionally includes a \(d_{yz}\) orbital, a \(d_{xz}\) orbital, a \(d_{xy}\) orbital, a \(d_{x^2-y^2}\) orbital, and a \(d_{3z^2-r^2}\) orbital.

The orbital labels \(s\), \(p_x\), and so so on, refer to the symmetry of the atomic orbitals under rotations along the relevant coordinate axes, and follow the convention of Ref. [SK54]. We note that the \(s^{\star}\) label refers to an orbital with the same continuous spherical symmetry as the \(s\) orbital, but at a higher energy level.

The \(sp^3\) TB model, which has \(N_O=4\) orbitals per atom, is the lightest computationally, but unfortunately predicts that silicon has a direct bandgap [SK54]. This incorrect prediction can be corrected in the \(sp^3s^{\star}\) TB model, which has \(N_O=5\) orbitals per atom, at minimal additional computational cost [KBB+00, VHD83]. The current state of the art is the \(sp^3d^5s^{\star}\) TB model, which has \(N_O=10\) orbitals per atom, as it most accurately describes the bandstructure (in terms of bandgap, valleys, effective masses, and so on) of semiconductors and how it varies under strain [NRT+09].

When spin–orbit coupling and/or the Zeeman effect are considered, spin states must be added to the basis set, which doubles the number of basis states. In this case, the basis states may thus be denoted by the ket \(|n,i,\sigma\rangle\), where \(\sigma\) indexes the spin state. These spin states are the spin-up and spin-down states along the \(z\) axis, and are respectively denoted by \(\uparrow\) and \(\downarrow\).

Overall, the total number of basis states in the TB model is

(5.10)\[\mathrm{Number~of~basis~states} = N \times N_O \times N_S\,,\]

where \(N\) is the number of atoms in the system, \(N_O\) is the number of atomic orbitals per atom, and \(N_S\) is the number of spin states per orbital.

Finally, we note that although atomic orbitals on neighboring atoms generally have a non-zero overlap (physically speaking), the TB model assumes an orthonormal basis set, so that

(5.11)\[\langle n,i,\sigma | n^{\prime},i^{\prime},\sigma^{\prime} \rangle = \delta_{n,n^{\prime}} \delta_{i,i^{\prime}} \delta_{\sigma,\sigma^{\prime}}\,.\]

This orthonormality condition is achieved using Löwdin’s method [Lowdin50], which preserves the symmetry of the atomic orbitals.

The TB basis set used to construct the TB Hamiltonian matrix is indirectly set through the set_tb_params method of the Atoms class; the TB basis set is automatically set according to the TB model parameter set by the user and defined in the qtcad.atoms.materials module.

Block structure of the Hamiltonian

To construct the TB Hamiltonian matrix of an atomic structure, we first need to define the ordering of the TB basis states:

first, the basis states are ordered according to atom index;
second, the basis states for each atom are ordered according to atomic orbital type, with the \(s\) orbital first, followed by the \(p_x\), \(p_y\), and \(p_z\) orbitals, then followed (if applicable) by the \(d_{yz}\), \(d_{xz}\), \(d_{xy}\), \(d_{x^2-y^2}\), and \(d_{3z^2-r^2}\) orbitals, and finally (if applicable) followed by the \(s^{\star}\) orbital;
third, the basis states for each orbital are ordered according to spin (if applicable), with the \(\uparrow\) spin first and \(\downarrow\) spin second.

In this ordering, the TB Hamiltonian matrix \(H\) is an \(\left(N N_O N_S\right) \times \left(N N_O N_S\right)\) matrix which may be partitioned into \(\left(N_O N_S\right) \times \left(N_O N_S\right)\) blocks:

(5.12)\[\begin{split}H = \begin{pmatrix} H_{11} & H_{12} & \cdots & H_{1N} \\ H_{21} & H_{22} & \cdots & H_{2N} \\ \vdots & \vdots & \ddots & \vdots \\ H_{N1} & H_{N2} & \cdots & H_{NN} \end{pmatrix}\,,\end{split}\]

where \(H_{nn^{\prime}}\), for \(n,n^{\prime}=1,\ldots,N\), are \(\left(N_O N_S\right) \times \left(N_O N_S\right)\) blocks. For \(n=n^{\prime}\), the block \(H_{nn}\) is referred to as an “onsite block”. For \(n\neq n^{\prime}\), the block \(H_{nn^{\prime}}\) is referred to as an “offsite block”.

QTCAD uses a nearest-neighbor TB model, so that the offsite block \(H_{nn^{\prime}}\) is non-zero if and only if the atoms indexed by \(n\) and \(n^{\prime}\) are nearest neighbors. In materials where the atoms are tetrahedrally-bonded (such as silicon, germanium, and silicon–germanium alloys), each atom (excluding surface atoms) has exactly four nearest neighbors, meaning that, for a given atom \(n\), \(H_{nn^{\prime}}\) is non-zero for only four values of \(n^{\prime}\), thereby making the TB Hamiltonian matrix relatively sparse.

Onsite blocks

The onsite block \(H_{nn}\) of the TB Hamiltonian is parametrized in terms of TB model parameters which only depend on the chemical species of the atom indexed by \(n\). These parameters are available in the qtcad.atoms.materials module, specifically in OneAtom objects.

In this section, we assume that \(N_S=1\); the effects of spin–orbit coupling and of the Zeeman effect will be described later.

We also note that only the \(sp^3d^5s^{\star}\) TB model can capture strain effects.

\(sp^3\) model

In the \(sp^3\) TB model, the onsite block is given by:

(5.13)\[\begin{split}H_{nn} = \begin{pmatrix} E_s & 0 & 0 & 0 \\ 0 & E_p & 0 & 0 \\ 0 & 0 & E_p & 0 \\ 0 & 0 & 0 & E_p \end{pmatrix}\,,\end{split}\]

where \(E_s\) and \(E_p\) are the onsite energies of the \(s\) and \(p\) orbitals of the atom index by \(n\), respectively.

\(sp^3s^{\star}\) model

Similarly, in the \(sp^3s^{\star}\) TB model, the onsite block is given by:

(5.14)\[\begin{split}H_{nn} = \begin{pmatrix} E_s & 0 & 0 & 0 & 0 \\ 0 & E_p & 0 & 0 & 0 \\ 0 & 0 & E_p & 0 & 0 \\ 0 & 0 & 0 & E_p & 0 \\ 0 & 0 & 0 & 0 & E_{s^{\star}} \end{pmatrix}\,,\end{split}\]

where \(E_{s^{\star}}\) is the onsite energy of the \(s^{\star}\) orbital of the atom indexed by \(n\).

\(sp^3d^5s^{\star}\) model

The onsite block has a more complex structure in the \(sp^3d^5s^{\star}\) TB model, as this model describes strain effects which couple orbitals [NRT+09]. Specifically:

(5.15)\[\begin{split}H_{nn} = \begin{pmatrix} H_{ss} & H_{sp} & H_{sd} & H_{ss^{\star}} \\ H_{ps} & H_{pp} & H_{pd} & H_{ps^{\star}} \\ H_{ds} & H_{dp} & H_{dd} & H_{ds^{\star}} \\ H_{s^{\star}s} & H_{s^{\star}p} & H_{s^{\star}d} & H_{s^{\star}s^{\star}} \end{pmatrix}\,.\end{split}\]

The sub-block \(H_{ss}\) is given by:

(5.16)\[H_{ss} = E_s + \alpha_s \frac{\Delta \Omega}{\Omega_0} \,,\]

where \(\alpha_s\) is a parameter that describes the dependence of the \(s\) orbital energy on hydrostatic strain and \(\frac{\Delta \Omega}{\Omega_0}\) is the relative change in unit cell volume around the atom indexed by \(n\), which, in this model, is approximated as:

(5.17)\[\frac{\Delta \Omega}{\Omega_0} = \frac{3}{4} \sum_{n^{\prime}} \frac{d_{nn^{\prime}}-d_0}{d_0}\,.\]

Here, the sum runs over the nearest neighbors of the atom, \(d_{nn^{\prime}}\) is the distance between the atoms indexed by \(n\) and \(n^{\prime}\), and \(d_0\) is the equilibrium distance between the atoms indexed by \(n\) and \(n^{\prime}\).

The sub-block \(H_{pp}\) is given by:

(5.18)\[\begin{split}H_{pp} = \left(E_p + \alpha_p \frac{\Delta \Omega}{\Omega_0}\right) I_3 + \sum_{n'} \beta_p\left(d_{nn'}\right) \begin{pmatrix} l^2 - \frac{1}{3} & lm & ln \\ ml & m^2 - \frac{1}{3} & mn \\ nl & nm & n^2 - \frac{1}{3} \end{pmatrix} \,,\end{split}\]

where \(\alpha_p\) is a parameter that describes the dependence of the \(p\) orbitals energies on hydrostatic strain and where \(I_3\) is the identity matrix. The variables \(l\), \(m\), and \(n\) are defined through:

(5.19)\[\mathbf{d}_{nn^{\prime}} = d_{nn^{\prime}} \left(l, m, n\right)\,.\]

Note that their dependence on \(n\) and \(n^{\prime}\) is implicit, for notational simplicity. In addition,

(5.20)\[\beta_p(d_{nn^{\prime}}) = \beta_p^0 + \beta_p^1 \frac{d_{nn^{\prime}} - d_0}{d_0}\]

captures uniaxial and biaxial strain effects.

Similar expressions hold for the other sub-blocks of \(H_{nn}\); their exact forms are detailed in Ref. [NRT+09].

Offsite blocks

The offsite block \(H_{nn^{\prime}}\) of the TB Hamiltonian is parametrized in terms of TB model parameters which only depend on the chemical species of the atoms indexed by \(n\) and \(n^{\prime}\). These parameters are available in the qtcad.atoms.materials module, specifically in TwoAtoms objects.

As in the previous section, we assume that \(N_S=1\).

Bond integrals and Slater–Koster rules

At the most basic level, the offsite block \(H_{nn^{\prime}}\) is parametrized in terms of bond integrals, which are labeled by \(V_{ii^{\prime}\zeta}\), where \(i\) labels an orbital on atom \(n\), \(i^{\prime}\) labels an orbital on atom \(n^{\prime}\), and \(\zeta\) labels the symmetry of the bond; \(\zeta\) can be either \(\sigma\), \(\pi\), or \(\delta\).

From there, the elements of the offsite block are given by the Slater–Koster rules, which harness the symmetry of the atomic orbitals to express the matrix elements in terms of the bond integrals and the bond direction parameters \(l\), \(m\), and \(n\). The Slater–Koster rules are too numerous to be exhaustively listed here; interested reader may consult Ref. [SK54] for a complete list. For illustration purposes, they include the following:

(5.21)\[\begin{split}\left\langle n,s | H | n^{\prime},s \right\rangle = V_{ss\sigma}\,, \\ \left\langle n,s | H | n^{\prime},p_x \right\rangle = l V_{sp\sigma}\,.\end{split}\]

Harrison exponents

In the \(sp^3d^5s^{\star}\) model, these equations are modified to include a bond-length dependence, so as to capture the decay of the overlap between atomic orbitals on neighboring atoms as a function of the distance that separates them. Specifically, the bond integrals acquire a bond-length dependence:

(5.22)\[V_{ii^{\prime}\zeta} \left(d_{nn^{\prime}} \right) = V_{ii^{\prime}\zeta} \left(d_0 \right) \left(\frac{d_0}{d_{nn^{\prime}}} \right)^{\eta_{ii^{\prime}}}\,,\]

where \(\eta_{ii^{\prime}}\) is a so-called Harrison exponent [FH79].

We note that the Harrison exponents are only available for the \(sp^3d^5s^{\star}\) model.

Spin–orbit coupling

To account for spin–orbit coupling in the TB Hamiltonian, the first step is to extend the TB basis states to include a spin component. According to the ordering of basis states described in Block structure of the Hamiltonian, this can be achieved by taking the Kronecker product of the spin-less TB Hamiltonian with the \(2\times 2\) identity matrix \(I_2\). For illustration purposes, this operation results in the following generalization of Eq. (5.13) for the onsite block in the \(sp^3s^{\star}\) model:

(5.23)\[\begin{split}\begin{align} H_{nn} &= \begin{pmatrix} E_s & 0 & 0 & 0 & 0 \\ 0 & E_p & 0 & 0 & 0 \\ 0 & 0 & E_p & 0 & 0 \\ 0 & 0 & 0 & E_p & 0 \\ 0 & 0 & 0 & & E_{s^{\star}} \end{pmatrix} \otimes I_2 \\ &= \begin{pmatrix} E_s & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & E_s & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & E_p & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & E_p & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & E_p & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & E_p & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & E_p & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & E_p & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & E_{s^{\star}} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & E_{s^{\star}} \end{pmatrix}\,. \end{align}\end{split}\]

The second step is to include the spin–orbit coupling term in the TB Hamiltonian. In QTCAD, spin–orbit coupling only affects the onsite blocks of the TB Hamiltonian, specifically the diagonal sub-blocks associated with the \(p\) orbitals. The matrix describing spin–orbit coupling is modeled via a single parameter \(\lambda\) and is given by:

(5.24)\[\begin{split}H_{\mathrm{SOC}} = \lambda \mathbf{L} \cdot \mathbf{S} = \lambda \begin{pmatrix} 0 & 0 & -i & 0 & 0 & 1 \\ 0 & 0 & 0 & i & -1 & 0 \\ i & 0 & 0 & 0 & 0 & -i \\ 0 & -i & 0 & 0 & -i & 0 \\ 0 & -1 & 0 & i & 0 & 0 \\ 1 & 0 & i & 0 & 0 & 0 \end{pmatrix}\,,\end{split}\]

where \(\mathbf{L}\) is the vector of orbital angular momentum operators for \(p\) orbitals and \(\mathbf{S}\) is the vector of spin operators.

To construct the TB Hamiltonian, \(H_{\mathrm{SOC}}\) is then simply added to the onsite diagonal sub-blocks associated with the \(p\) orbitals. For example, in the \(sp^3s^{\star}\) model, the onsite block is then given by:

(5.25)\[\begin{split}H_{nn} = \begin{pmatrix} E_s & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & E_s & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & E_p & 0 & -i\lambda & 0 & 0 & \lambda & 0 & 0 \\ 0 & 0 & 0 & E_p & 0 & i\lambda & -\lambda & 0 & 0 & 0 \\ 0 & 0 & i\lambda & 0 & E_p & 0 & 0 & -i\lambda & 0 & 0 \\ 0 & 0 & 0 & -i\lambda & 0 & E_p & -i\lambda & 0 & 0 & 0 \\ 0 & 0 & 0 & -\lambda & 0 & i\lambda & E_p & 0 & 0 & 0 \\ 0 & 0 & \lambda & 0 & i\lambda & 0 & 0 & E_p & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & E_{s^{\star}} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & E_{s^{\star}} \end{pmatrix}\,,\end{split}\]

We note that spin–orbit coupling is only available for the \(sp^3s^{\star}\) and \(sp^3d^5s^{\star}\) models.

Enabling spin–orbit coupling may be done via the set_soc method of the Atoms class.

Treatment of surface atoms

The method to construct the TB Hamiltonian described so far does not differentiate between bulk atoms, which are tetrahedrally-bonded to their four nearest neighbors, and surface atoms, which are bonded to fewer than four atoms. This is problematic, as surface atoms have dangling bonds, which, without proper treatment, result in midgap, surface states [HTJ+16].

There exists two common approaches to treat surface atoms in the TB model:

by hydrogen passivation, which consists in saturating the dangling bonds of surface atoms with hydrogen atoms;
by raising the dangling bond energies so that no surface states is present in the bandgap.

Hydrogen passivation suffers from three major drawbacks:

it requires additional TB parameters, such as onsite energies for hydrogen atoms and bond integrals between hydrogen atoms and the other atoms;
it increases the size of the TB Hamiltonian matrix;
in certain cases, it can lead to incorrect predictions; for example, it has been shown to lead to an overestimation of the performance of nanoscale transistors [HTJ+16].

As such, QTCAD only supports the second approach, which is described in Ref. [LOVAK04]. Succinctly, this approach only affects the onsite blocks of the TB Hamiltonian matrix and involves the following steps:

detect surface atoms;
if the atom has two (three) nearest neighbors, construct a set of four \(sp^3\) hybridized orbitals, two (three) of which are aligned with the physical bonds, and two (one) of which is determined by the tetrahedral bonding symmetry;
this results in a basis set of \(sp^3\) hybridized orbitals, \(\left\{|sp^3_a\rangle, |sp^3_b\rangle, |sp^3_c\rangle, |sp^3_d\rangle\right\}\);
raise the energy of these hyrbidized orbitals by constructing the matrix \(\delta_{sp^3} I_4\), where \(\delta_{sp^3}\) is a parameter that describes the dangling bond energy shift and \(I_4\) is the \(4\times 4\) identity matrix;
transform this matrix into the original orbital basis \(\left\{|s\rangle, |p_x\rangle, |p_y\rangle, |p_z\rangle\right\}\);
add the resulting matrix to the diagonal sub-block associated with the \(sp^3\) orbitals of the appropriate onsite block of the TB Hamiltonian matrix.

In general, setting \(\delta_{sp^3}\) to a value larger than \(5~\mathrm{eV}\) is sufficient to eliminate all midgap, surface states; it may be set through the energy_shift_db attribute of the SolverParams class for the TB Schrödinger Solver.

External electric potential and multiscale simulations

To model realistic quantum dot devices, it is necessary to include the effect of an external electric potential in the TB Hamilonian. Such an external potential may arise from the presence of a gate voltage, and may lead to confinement of the wavefunction in the quantum dot.

In QTCAD, the external electric potential \(\varphi\left(\mathbf{r}\right)\) may be set in two ways:

from a user-specified function, which is passed to the set_potential method of the Atoms object;
from the results of a finite-element Poisson solver, which is done by passing an appropriate Device object to the set_potential_from_device method of the Atoms object.

The following diagonal matrix is then added to the TB Hamiltonian matrix:

(5.26)\[\begin{split}H_{\varphi} = -e \begin{pmatrix} \varphi\left(\mathbf{r}_1\right) I_{N_O N_S} & 0 & \cdots & 0 \\ 0 & \varphi\left(\mathbf{r}_2\right) I_{N_O N_S} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \varphi\left(\mathbf{r}_N\right) I_{N_O N_S} \end{pmatrix}\,,\end{split}\]

where \(e\) is the elementary charge, \(\mathbf{r}_i\) is the position of the atom indexed by \(i\), and \(I_{N_O N_S}\) is the \(\left(N_O N_S\right) \times \left(N_O N_S\right)\) identity matrix. The computation of \(\varphi\left(\mathbf{r}_i\right)\) is done either via direct evaluation of the user-specified function or via interpolation of the finite-element Poisson solver solution. In the latter case, the simulation may be considered to be a multiscale simulations, since it would combine the finite-element method to resolve the electrostatics of the device on micrometer lengthscales with the TB method to resolve the quantum electronic structure of the device on atomic lengthscales.

We note that unlike the finite-element Schrödinger solver, the TB Schrödinger solver does not explicitly model holes. Only electrons are considered in the TB model. Holes may be indirectly modeled by studying electrons near the top of the valence band.

In addition, we stress that the external potential \(\varphi\left(\mathbf{r}\right)\) in Eq. (5.26) is the electric potential, not the confinement potential. In particular, it does not include any band offset; such phenomena are indirectly captured through the orbital onsite energies (\(E_s\), \(E_p\), and so on) of the TB model.

Finally, we note that care should be taken when setting the external electric potential to ensure that it properly confines charge. When modeling electrons, \(\varphi\left(\mathbf{r}\right)\) should be concave down; when modeling holes, \(\varphi\left(\mathbf{r}\right)\) should be concave up.

Note

The gauge for the electric potential and energy in TB simulations performed in the qtcad.atoms module is the default gauge of the finite-element simulations performed in the qtcad.device module, and is described in Band alignment in heterostructures. Specifically, the zero of energy is the Fermi energy of intrinsic, unstrained, bulk silicon at the device temperature (if the potential is set through the set_potential_from_device method) or at absolute zero (if the potential is set through the set_potential method). In the latter case, the zero of energy is then at the middle of the bandgap of intrinsic, unstrained, bulk silicon.

Magnetic effects

Magnetic fields are likewise important to consider for the modeling of realistic spin qubits, as these devices are often operated in the presence of an external magnetic field to achieve spin splitting and/or for quantum control.

Magnetic fields affect the TB Hamiltonian in two ways:

they couple to the spin and orbital angular momentum around each atomic site via the Zeeman effect by modifying the onsite blocks of the TB Hamiltonian;
they couple to the orbital degrees of freedom via the minimal coupling Hamiltonian through a Peierls substitution by modifying the offsite blocks of the TB Hamiltonian.

These effects are described in more detail below.

Zeeman effect

The Zeeman Hamiltonian of an electron, which in the context of the TB model leads to an additional contribution to the onsite blocks of the TB Hamiltonian matrix, is given by the following expression:

(5.27)\[H_Z = \frac{\mu_B}{\hbar} \left(g_s \mathbf S + g_l\mathbf L\right) \cdot \mathbf B\,.\]

Here, \(\mu_B = \frac{e\hbar}{2m_e}\) is the Bohr magneton (\(m_e\) is the electron rest mass), \(g_s=2\) is the electron \(g\)-factor for the spin, \(g_l=1\) is the electron \(g\)-factor for the orbital angular momentum, \(\mathbf S\) and \(\mathbf L\) are the spin and orbital angular momentum of the electron, and \(\mathbf B\) is the magnetic field. The electron spin operator can be written using the Pauli matrices:

(5.28)\[\mathbf S = \frac{\hbar}{2} \boldsymbol{\sigma}\,,\]

where the components of the vector of matrices \(\boldsymbol{\sigma}\) are given by:

(5.29)\[\begin{split}\sigma_x = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, \quad \sigma_y = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix}, \quad \sigma_z = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}\,.\end{split}\]

The orbital angular momentum operator is less trivial to express and depends on the TB basis being used. For example, in the \(sp^3d^5s^{\star}\) TB model, the basis is given by \(\{s, p_x, p_y, p_z, d_{yz}, d_{xz}, d_{xy}, d_{x^2-y^2}, d_{3z^2-r^2}, s^{\star}\}\), where \(s\), \(p\), and \(d\) refer to the total angular momentum around each atomic site and correspond to total angular momentum \(\ell=0, 1, 2\), respectively, and the subscripts refer to the transformation properties of the basis under symmetry operations of the crystal. Since the states with \(\ell=0\) are spherically symmetric and invariant under all symmetry operations, the \(s\) and \(s^{\star}\) orbitals possess no subscript. The TB basis orbitals can be expressed in terms of the eigenstates of the orbital angular momentum matrix \(L_z\) along \(z\) as follows:

(5.30)\[\begin{split}\begin{align} s &= U_0 Y_{0,0} = Y_{0,0}\,, \\ \begin{pmatrix} p_x \\ p_y \\ p_z \end{pmatrix} &= U_1 \begin{pmatrix} Y_{1,1} \\ Y_{1,-1} \\ Y_{1, 0} \end{pmatrix} = \begin{pmatrix} -\frac{1}{\sqrt{2}} & 0 & \frac{1}{\sqrt{2}} \\ \frac{i}{\sqrt{2}} & 0 & \frac{i}{\sqrt{2}} \\ 0 & 1 & 0 \end{pmatrix} \begin{pmatrix} Y_{1,1} \\ Y_{1,-1} \\ Y_{1, 0} \end{pmatrix}\,, \\ \begin{pmatrix} d_{yz} \\ d_{xz} \\ d_{xy} \\ d_{x^2-y^2} \\ d_{3z^2-r^2} \end{pmatrix} &= U_2 \begin{pmatrix} Y_{2,2} \\ Y_{2,1} \\ Y_{2,0} \\ Y_{2, -1} \\ Y_{2, -2} \end{pmatrix} = \begin{pmatrix} 0 & \frac{i}{\sqrt{2}} & 0 & \frac{i}{\sqrt{2}} & 0 \\ 0 & -\frac{1}{\sqrt{2}} & 0 & \frac{1}{\sqrt{2}} & 0 \\ -\frac{i}{\sqrt{2}} & 0 & 0 & 0 & \frac{i}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} & 0 & 0 & 0 & \frac{1}{\sqrt{2}} \\ 0 & 0 & 1 & 0 & 0 \end{pmatrix} \begin{pmatrix} Y_{2,2} \\ Y_{2,1} \\ Y_{2,0} \\ Y_{2, -1} \\ Y_{2, -2} \end{pmatrix}\,, \\ s^{\star} &= U_0 Y_{0, 0}^{\star} = Y_{0, 0}^{\star}\,, \end{align}\end{split}\]

where \(Y_{\ell,m}\) is the simultaneous eigenstates of total angular momentum with quantum number \(\ell\) and \(z\)-component of angular momentum with quantum number \(m\), and \(U_{\ell}\) is a unitary basis transformation matrix.

To complete the description of the coupling between the magnetic field and the onsite orbital angular momentum, we also need to specify the orbital angular momentum operators in the \(Y_{\ell,m}\) basis. For \(\ell=0\):

(5.31)\[\mathbf L^0 = 0\,.\]

For \(\ell=1\), the orbital angular momentum operator is given by:

(5.32)\[\begin{split}\begin{align} L^1_x &= \frac{\hbar}{\sqrt{2}}\begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 0 \end{pmatrix}\,, \\ L^1_y &= \frac{\hbar}{\sqrt{2}}\begin{pmatrix} 0 & -i & 0 \\ i & 0 & -i \\ 0 & i & 0 \end{pmatrix}\,, \\ L^1_z &= \hbar\begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & -1 \end{pmatrix}\,. \end{align}\end{split}\]

For \(\ell=2\), the orbital angular momentum operator is given by:

(5.33)\[\begin{split}\begin{align} L^2_x &= \frac{\hbar}{2}\begin{pmatrix} 0 & 2 & 0 & 0 & 0 \\ 2 & 0 & \sqrt{6} & 0 & 0 \\ 0 & \sqrt{6} & 0 & \sqrt{6} & 0 \\ 0 & 0 & \sqrt{6} & 0 & 2 \\ 0 & 0 & 0 & 2 & 0 \end{pmatrix}\,, \\ L^2_y &= \frac{\hbar}{2i}\begin{pmatrix} 0 & 2 & 0 & 0 & 0 \\ -2 & 0 & \sqrt{6} & 0 & 0 \\ 0 & -\sqrt{6} & 0 & \sqrt{6} & 0 \\ 0 & 0 & -\sqrt{6} & 0 & 2 \\ 0 & 0 & 0 & -2 & 0 \end{pmatrix} \,, \\ L^2_z &= \hbar\begin{pmatrix} 2 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 & -2 \end{pmatrix} \,. \end{align}\end{split}\]

We now have all the elements to write the orbital angular momentum operator in the \(sp^3d^5s^{\star}\) basis:

(5.34)\[\begin{split}\mathbf L = \begin{pmatrix} U^\dagger L_x U \\ U^\dagger L_y U \\ U^\dagger L_z U \end{pmatrix} \,,\end{split}\]

where

(5.35)\[\begin{split}U = \begin{pmatrix} U_0 & 0 & 0 & 0 \\ 0 & U_1 & 0 & 0 \\ 0 & 0 & U_2 & 0 \\ 0 & 0 & 0 & U_0 \end{pmatrix} \,,\end{split}\]

(5.36)\[\begin{split}L_j = \begin{pmatrix} L_j^0 & 0 & 0 & 0 \\ 0 & L_j^1 & 0 & 0 \\ 0 & 0 & L_j^2 & 0 \\ 0 & 0 & 0 & L_j^0, \end{pmatrix} \,,\end{split}\]

and \(j\;\in\;\{x,y,z\}\).

A similar procedure can be followed to compute the Zeeman Hamiltonian for the \(sp^3\) and \(sp^3s^{\star}\) TB models.

The Zeeman Hamiltonian obtained via this procedure is then added to all the onsite blocks of the TB Hamiltonian matrix. The Zeeman effect can be enabled in a simulation by calling the set_Zeeman method of the Atoms object. The magnetic field can be set through the set_Bfield method.

Magnetic orbital effects—Peierls substitution

In addition to the Zeeman effect, magnetic fields also couple to the orbital degrees of freedom via the minimal coupling Hamiltonian. This is done, effectively, by replacing the momentum operator \(\mathbf p\) in the TB Hamiltonian as follows:

(5.37)\[\mathbf p \rightarrow \mathbf p + e \mathbf A\left(\mathbf{r}\right)\,.\]

Here, \(\mathbf A\left(\mathbf{r}\right)\) is the vector potential associated with the magnetic field \(\mathbf B\):

(5.38)\[\mathbf B = \nabla \times \mathbf A\left(\mathbf{r}\right)\,.\]

Currently, QTCAD only supports uniform magnetic fields, which can be set through the set_Bfield method of the Atoms object. For a uniform magnetic field, the vector potential, in the symmetric gauge, is given by:

(5.39)\[\mathbf A\left(\mathbf{r}\right) = \frac{1}{2} \mathbf B \times \mathbf{r}\,.\]

Note that physical observables are invariant under the choice of gauge.

In the TB model, Eq. (5.37) is incorporated by modifying the offsite blocks \(H_{nn^{\prime}}\) of the TB Hamiltonian as outlined in Ref. [GV95]:

(5.40)\[H_{nn^{\prime}} \rightarrow H_{nn^{\prime}} e^{ \frac{i}{\hbar} \left({\mathbf r_{n^{\prime}}} - {\mathbf r_{n}}\right) \cdot \frac{1}{2}\left[ \mathbf A(\mathbf r_{n^{\prime}}) + \mathbf A(\mathbf r_{n}) \right] } \,.\]

This modification of the offsite blocks is known as the Peierls substitution with the multiplicative phase factor known as the Peierls phase. The Peierls substitution is automatically applied when the magnetic field set through the set_Bfield method is non-zero.