# 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

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

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\).

Note

The
`get_operator_matrix_element`

can be used to compute individual matrix element.

## Perturbations

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

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

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

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

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)
Z.solve()
```

#### Electrons

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

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\).

#### Holes

The `Zeeman`

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

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,

and

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)
G.solve()
```

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

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)
H.solve()
```

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

and the `solve`

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