4. Operators

Often it is interesting to evaluate the matrix elements of certain position-dependent operators in the eigenbasis of a Hamiltonian. The operators module facilitates this calculation.

We consider the following time-independent Schrödinger equation

(4.1)\[H_0 \left|n\right> = E_n\left|n\right>,\]

where \(H_0\) is an unperturbed Hamiltonian, \(E_n\) are the eigenenergies, and \(\left|n\right>\) are the eigenstates. The goal of the Schrödinger solver is to solve for the eigenenergies, \(E_n\), and the eigenfunctions, \(\left<\mathbf{r} |n\right>\), of Eq. (4.1) for all \(0 \le n < N_0\), where \(N_0\) is a solver parameter that is set by the num_states attribute of the Schrödinger solver. The matrix elements of an operator, \(U(\mathbf{r})\), are given by

(4.2)\[\left(\hat{U}\right)_{nm} = \left<n|U|m\right>\]

where we have used the notation \(\hat{O}\) to refer to a matrix representation of the operator \(O\) and \(n, m \le N_0\) run over the eigenstates of \(H_0\). The matrix elements can be evaluated using the Operator class. For example, the matrix elements of the position operator \(x\) can be obtained by running

from device.operators import Operator

def pos_x_func(x, y, z):
    return x

pos_x = Operator(d, pos_x_func)
x_matrix = pos_x.get_operator_matrix()

The Operator constructor takes as input a device, d, and a function of position, pos_x_func. The get_operator_matrix method outputs a 2D numpy array whose matrix elements are given by Eq. (4.2). In this specific example, pos_x_func represents the position operator \(x\) and x_matrix will contain the matrix elements of \(x\).


The get_operator_matrix_element can be used to compute individual matrix element.


In the Schrödinger theory section , we saw that QTCAD devices are equipped with various setters that allow the Schrödinger solver to account for effects such spin–orbit coupling, Zeeman effect, orbital–magnetic coupling, and strain, when solving for the eigenstates and eigenenergies of a confined system. While these effects can be included directly at the level of the Schrödinger solver, it is sometimes interesting to consider perturbations which are accounted for after an unperturbed Hamiltonian has been diagonalized.

We again consider the following time-independent Schrödinger equation

(4.3)\[H_0 \left|n\right> = E_n\left|n\right>,\]

where \(H_0\) is an unperturbed Hamiltonian, \(E_n\) are the eigenenergies, and \(\left|n\right>\) are the eigenstates. In its eigenbasis, the matrix form of \(H_0\) has matrix elements

(4.4)\[\left(\hat{H}_0\right)_{nm} = E_n\delta_{nm},\]

where \(n, m \le N_0\) run over the eigenstates of \(H_0\). We now consider a perturbing term, \(U\). In the eigenbasis of \(H_0\), the matrix elements of \(\hat{U}\) are given by

(4.5)\[\left(\hat{U}\right)_{nm} = \left<n|U|m\right>\]

Diagonalizing the \(N_0 \times N_0\) matrix, \(\hat{H}_0 + \hat{U}\) will lead eigenvectors, \(\vec{\alpha}_i\), and eigenvalues, \(E_{i}\). The corresponding states

(4.6)\[\left|\psi_i\right> = \sum_{n=1}^{N_0}\alpha_i^n \left|n\right>,\]

where \(\alpha_i^n\) is the \(n^{\mathrm{th}}\) component of \(\vec{\alpha}_i\), approximate the eigenstates of the full Hamiltonian, \(H_0 + U\), while the eigenvalues, \(E_{i}\), approximate the eigenenergies.

The procedure outlined above for getting approximate eigenenergies and eigenfunctions for Hamiltonians of the type \(H_0 + U\) has been implemented in QTCAD via the Hamiltonian Operator class. The constructor for this class takes as input any object that possesses energies, eigenfunctions, and mesh attributes (this includes device objects and other Hamiltonian Operator objects) as well as some parameters which define the perturbation. These objects are equipped with a solve method which constructs and diagonalizes the Hamltonian, \(\hat{H}_0 + \hat{U}\). The eigenenergies and eigenfunctions are stored in the attributes, energies and eigenfunctions, respectively. In addition, the matrix \(\hat{U}\) is stored in the U_mat attribute. Currently, there are 3 distinct perturbations implemented within QTCAD:

Zeeman Operator

The Zeeman effect describes the coupling of magnetic fields to the spin degree of freedom of a particle. Like all operators, the Zeeman operator constructor takes in an object that possesses energies, eigenfunctions, and mesh attributes. In addition, the magnetic field, B, and the g-factor, g must also be specified (g is an optional argument and if it is not specified, the first input must be a device object which will supply g to the perturbation automatically).

from device.operators import Zeeman

Z = Zeeman(d, B, g=None)


For electrons the Zeeman Hamiltonian is assumed to be of the form [Win03]

(4.7)\[U = H_Z = \frac{\mu_B}{\hbar} \mathbf S \cdot \overleftrightarrow{\mathbf g} \cdot \mathbf B,\]

where \(\mu_B = \frac{e\hbar}{2m_e}\) is the Bohr magneton (\(e > 0\) is the electron charge and \(m_e\) is the electron rest mass), \(\overleftrightarrow{\mathbf g}\) is the effective electron g-tensor (\(3\times3\)) in the system under consideration, \(\mathbf B\) is the magnetic field, and \(\mathbf S\) is the electron spin operator. Both the magnetic field, \(\mathbf B\), and the g-tensor, \(\overleftrightarrow{\mathbf g}\), can be position dependent. Consequently, B can be given as a function of the spatial coordinates (\(x\), \(y\), and \(z\)) which outputs a 1D \(1\times 3\) numpy array, a 1D numpy array of length 3 giving a constant magnetic field over all space, or a 2D numpy array where the first index runs over the global nodes of the mesh over which the device is defined and the second over the components of the magnetic field. Similarly, g can be given as a function that outputs a 2D \(3\times 3\) numpy array, a 2D numpy array of dimensions \(3\times 3\) giving a constant g-tensor over all space, or a 3D numpy array where the first index runs over the global nodes and the second and third indices over the components of the g-tensor.

In the case where the confined carriers are electrons, and no spin-orbit coupling or magnetic field is set throughout the device, the eigenfunctions will have no spin index. The Zeeman perturbation will therefore add a spin index to the states before computing the matrix elements of \(H_Z\).


The Zeeman perturbation class generalizes the Zeeman for holes within the Luttinger-Kohn formalism [see Eq. (3.30)],

(4.8)\[U = H_Z = - 2 \frac{\mu_B}{\hbar} \mathbf J \cdot \overleftrightarrow{\mathbf g}_{\kappa} \cdot \mathbf B - 2 \frac{\mu_B}{\hbar} \mathbf {\cal J} \cdot \overleftrightarrow{\mathbf g}_{q} \cdot \mathbf B.\]

Here, \(\mathbf J\) is the vector of spin-3/2 matrices [its components are given in Eq. (3.15)] and \(\mathbf {\cal J} = (J_x^3, J_y^3, J_z^3)\). In the standard Luttinger-Kohn formalism,

(4.9)\[\begin{split}\overleftrightarrow{\mathbf g}_{\kappa} = \left(\begin{array}{lll} \kappa & 0 & 0\\ 0 & \kappa & 0\\ 0 & 0 & \kappa\\ \end{array} \right),\end{split}\]


(4.10)\[\begin{split}\overleftrightarrow{\mathbf g}_{q} = \left(\begin{array}{lll} q & 0 & 0\\ 0 & q & 0\\ 0 & 0 & q\\ \end{array} \right),\end{split}\]

where, \(\kappa\) and \(q\) are material parameters. In general however, the Zeeman perturbation permits arbitrary \(\overleftrightarrow{\mathbf g}_{\kappa}\) and \(\overleftrightarrow{\mathbf g}_{q}\). The \(g\)-tensors can be specified through the g keyword argument by providing it with a 2-tuple, (g_kappa, g_q), or a single tensor, g_kappa, in which case we take g_q=0. In the case where the default value is taken for g, values for \(\kappa\) and \(q\) will automatically be generated from the device input and the tensors, \(\overleftrightarrow{\mathbf g}_{\kappa}\) and \(\overleftrightarrow{\mathbf g}_{q}\) will be taken to have the standard Luttinger–Kohn form [Eqs. (4.9) and (4.10)]. The magnetic field is specified in the same way as for electrons.

Gate Operator

Sometimes, we wish to understand how a change on a boundary condition can affect our system. The Gate operator is appropriate in this context. In addition to a device-like object, the constructor for this object has two arguments. The first is a string or list of strings which identifies boundaries over which we wish to change the boundary conditions. The second is a scalar or list of scalars that specifies the new biases to apply to these gates.

from device.operators import Gate

G = Gate(d, gate, V, phys_d=None, base_bias=None)

There is also an important keyword argument, phys_d, which specifies the physical device over which the boundary conditions are defined. By default the physical device is chosen to be d itself (the first argument). However, in the case where d is another operator or a subdevice, phys_d may be different than d. We can also specify a base_bias which tells the Gate operator the biases applied to the gates in the default configuration, before the change in boundary condition. Internally, the Gate operator will compute a potential \(V(\mathbf r)\) associated with these biases. If base_bias is set to None, \(V(\mathbf r)\) is taken to be the potential stored in phys_d.

This operator will compute the potential energy, \(\tilde V(\mathbf r)\), respecting the new boundary conditions. We can write a potential-energy perturbation as

(4.11)\[U = \Delta V(\mathbf r) = \tilde{V}(\mathbf r) - V(\mathbf r).\]

The solve method will diagonalize \(\hat{H_0}+\hat{U}\) where \(U\) is given by Eq. (4.11).

Hamiltonian Operator

Finally, we can also compute matrix elements of arbitrary operators defined over all space. The Hamiltonian Operator can be used for this purpose

from device.operators import HamiltonianOperator

H = HamiltonianOperator(d, U)

For an \(n\)-band system (e.g. \(n=1\) for electrons) U should be an \(n\times n\) matrix. Here, U is either a function of space or a numpy array defined over all global nodes. For this specific type of operator, we have

\[U = U(\mathbf r),\]

and the solve method will diagonalize \(\hat{H_0}+\hat{U}\).