is the volume of the grid in $k$-space over which the calculation is
performed. For the Landau-level mode, the volumetric offset density is instead
given by

$n_{\mathrm{offset}}=n_{z}[(n_{\mathrm{max}}+1)+(n_{\mathrm{max}}+2)]\frac{eB}{h}\frac{1}{n_{z}\Delta
z},$ (135)

where $n_{\mathrm{max}}$ is the maximum Landau level index (set by the command
line argument llmax) and $B$ is the magnetic field.

The full-band envelope approach not only considers the active carriers in the
quantum well but also results in interface dipoles and side bumps at the model
edges (see Fig. 18). Note that the algorithm described in 3.11.1 keeps the
total carrier density in the system constant, which includes these features.
Hence the carrier density inside the quantum well will be different from the
one specified by the argument cardens.

The implementation of the full-band approach is done by the SelfConSolver
derived classes SelfConSolverFullDiag and SelfConSolverLLFullDiag for
dispersion and Landau-level mode, respectively. By default, the full-band
approached is used for the self-consistent Hartree method. To use the naive
method, for gaining speed at the expense of becoming less physically sound,
one may set the configuration value `selfcon_full_diag=false`.

Figure 18: (a) Self-consistently calculated Hartree potential for a structure
with $70\,\mathrm{nm}$ thick HgTe layer and (Hg,Cd)Te barriers. The total
carrier density has been set as $0.004\,\mathrm{nm}^{-2}$. The boundary
conditions are determined by imposing an electric field strength of
$0.48\,\mathrm{mV}\,\mathrm{nm}^{-1}$ inside the bottom barrier. This
calculation has been performed with the default setting
selfcon_full_diag=true.

####  3.11.3 Interpretation of the results

Due to the combination of diagonalization and integration (in the steps
solving the Schrödinger and the Poisson equation, respectively), the iterative
dynamics of $V_{\mathrm{H}}^{(i)}$ as function of iteration step $i$ is
difficult to characterize and to predict. The convergence behaviour can depend
sensitively on input parameters such as boundary conditions and the initial
time step. If the solver converges, the result only represents a single
solution, out of possibly many solutions. Whether the result is the most
natural solution can only be determined by critical scrutiny.

When one relies on the self-consistent method, it is advisable to run the
calculation multiple times with slightly different settings, in order to see
if the algorithm is stable against such variations. The terminal output from
kdotpy can also be a useful information source, especially for tracking
convergence behaviour. To allow for systematic analysis of the iterative
dynamics, the self-consistent solvers can provide debug output, with the
internal state (carrier density $\rho(z)$ and Hartree potential
$V_{\mathrm{H}}(z)$) being written to a CSV file at each iteration. A thorough
systematic analysis of the self-consistent Hartree method is beyond the scope
of this work and will be addressed elsewhere.

##  4 Installation

###  4.1 Prerequisites

kdotpy is packaged as a Python package that is intended to be installed using
PIP. For installing kdotpy, you must have a working Python installation
(version 3.9.0 or above) and an up-to-date version of PIP. When PIP installs
kdotpy, it automatically checks if the packages numpy, scipy, and matplotlib
are installed and have compatible versions. Packages that kdotpy depends on
optionally may be installed manually with PIP.

If your Python installation already contains a variety of packages, the
dependencies of the installed packages could potentially conflict with those
of kdotpy. This problem may be avoided by installing kdotpy in a virtual
environment, which is like a ‘sandbox’ where kdotpy and its dependencies do
not interfere with other packages. The virtual environment can be created with
python3 -m venv _directoryname_. After activating the virtual environment, you
may now install kdotpy following any of the installation methods below.
Consult the documentation of the Python venv module for more information555See
https://docs.python.org/3/library/venv.html..

Installation methods 2 and 3 below rely on git, which you must have installed
in order for these methods to work. If you do not wish to use git, you may
wish to use one of the other methods.

###  4.2 Download and installation

We list four methods for installing kdotpy. You can choose any one of them,
but combining them is not recommended.

  1. 1.

_Using PIP, from the Python Package Index (PyPI):_

        python3 -m pip install kdotpy

  2. 2.

_Using PIP, from our Gitlab repository:_

        python3 -m pip install<EMAIL_ADDRESS>

  3. 3.

_Using git clone and installing it from your local copy:_

        git clone https://git.physik.uni-wuerzburg.de/kdotpy/kdotpy.git
    python3 -m pip install ./kdotpy

Note that this clones the files into the subdirectory kdotpy of your current
working directory.

  4. 4.

_Using a manual download and installing it from your local copy:_
If you prefer to avoid the PyPI and git, then you can download a .zip or
.tar.gz file from our Gitlab repository at https://git.physik.uni-
wuerzburg.de/kdotpy/kdotpy [42]. Download the preferred file format from the
dropdown menu under the ‘Code’ button. Unpack the files into an empty
directory. Then use

        python3 -m pip install ./directoryname

to install kdotpy, replacing ./directoryname by the appropriate directory name
where you unpacked the files.

You may obtain version 1.0.0, the exact version discussed in this article,
also from the Gitlab repository of this journal.

###  4.3 Installation for active developers

If you wish to make modifications to the kdotpy code, then it is recommended
to use PIP’s -e option (short for \--editable). This option dynamically links
the package installation directory to the source files. If you change the
source, you thus do not have to re-install the package. To use an editable
install, use

    git clone https://git.physik.uni-wuerzburg.de/kdotpy/kdotpy.git
    python3 -m pip install -e ./kdotpy

If you intend to actively contribute to the project, you must change to a non-
protected branch in order to be able to push your changes to the repository.
(Where you are allowed to push depends on your membership role in our Gitlab
project. The master branch is protected for all except the project owner. The
journal repository shall not be used for development purposes.)

###  4.4 Testing the installation

If the install was successful, then entering

    kdotpy version

on the command line should return the version number. If not, then something
went wrong during the installation and you may wish to try one of the other
methods. (This assumes that if you use a virtual environment, it has been set
up and activated correctly.)

For testing the functionality of kdotpy, you can use

    kdotpy test

for running the standardized tests. Note that the tests may take a few minutes
to complete. They will generate output in the subdirectory test relative to
where you start kdotpy test.

##  5 Usage

###  5.1 General remarks

kdotpy is designed as a standalone application. If you have followed the
installation instructions above, you can simply run kdotpy from the command
line, followed by the ‘subprogram’ label and further arguments. The first
argument is always the subprogram, but the order of the following arguments is
usually unimportant. (For an overview of the subprograms, see Sec. 3.1.2.) You
can run kdotpy from any folder.

Alternatively, you can also use python3 -m kdotpy followed by the subprogram
and further arguments. You can also import specific functions or submodules of
kdotpy from other Python scripts or the interactive interpreter using from
kdotpy import _function_. As we have designed kdotpy primarily as a command-
line tool, such imports are not recommended for normal use.

Below, we provide detailed tutorials that describe how to set up a calculation
in kdotpy step by step, and how to interpret the output. The curious user may
also find the tests defined by kdotpy test useful; the command lines for these
standardized tests can be extracted using kdotpy test showcmd.

###  5.2 Example calculation: Basic dispersion

####  5.2.1 Introduction

This basic tutorial describes a systematic manner to compose the command line
from scratch, adding the necessary arguments step by step. As an example, we
will calculate the dispersion of a $7\,\mathrm{nm}$ HgTe quantum well: We
assume the substrate is $\mathrm{Cd}_{0.96}\mathrm{Zn}_{0.04}\mathrm{Te}$ and
the barriers are $\mathrm{Hg}_{0.32}\mathrm{Cd}_{0.68}\mathrm{Te}$. We will
calculate the dispersion along the diagonal axis (110) and try to answer the
following questions about the maxima of the dispersion at finite momentum in
the valence band (also known as _camel back_):

  * •

Where is the finite-momentum maximum in momentum and energy?

  * •

Do we have a direct or an indirect gap?

  * •

What is the orbital character at the finite-momentum maximum?

####  5.2.2 Setting it up step by step

  1. 1.

We first have to determine the geometry and/or mode, defined by the number of
translationally invariant dimensions. In this example, we have a ‘2D’
geometry, with momentum coordinates $(k_{x},k_{y})$. Thus, the appropriate
subprogram to use here is kdotpy 2d.

        kdotpy 2d

  2. 2.

Next, we determine the number of orbitals in the model. This is usually eight
orbitals, for which we use the option 8o. Alternatively 6o for six orbitals
may be used. Let us also specify that we do not want to use the axial
approximation by entering noax.

        kdotpy 2d 8o noax

  3. 3.

Let us now enter the substrate material and ‘layer stack’. We use msubst for
the substrate material, mlayer for the layer materials and llayer for the
layer thicknesses. We also enter the resolution of the discretization in the
$z$ direction with zres; $0.25\,\mathrm{nm}$ is a good value to start with.

        kdotpy 2d 8o noax msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 7 10 zres 0.25

Note that we have taken the barriers to be 10 nm thick. Usually, even if the
barriers are much larger in reality, setting the thickness to
$10\,\mathrm{nm}$ gives an accurate representation of the dispersion while
reducing calculation time. Note also how the materials are entered.

  4. 4.

We define the parametrization of the momentum. Here, we take 60 points along
the chosen momentum axis, symmetric around zero, with `k -0.6 0.6 / 120`. The
(110) axis has a $45$ degree angle with the $k_{x}$ axis, hence we enter `kphi
45`.

        kdotpy 2d 8o noax msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 7 10 zres 0.25 k -0.6 0.6 / 120 kphi 45

If we would run this command, we already get a result, `dispersion.pdf`. It
looks promising, but kdotpy indicates that there are a few problems:

        Calculating bands (k=0)...
    Warning (band_type): Unable to determine band character and/or number of
    nodes for 50 eigenstates.
    Possible causes: spin degeneracy not broken, nonzero potential,
    one-dimensional geometry, etc.
    1 / 1
    ERROR (estimate_charge_neutrality_point): Failed, because E+ and/or E-
    bands are missing

  5. 5.

The assignment of band characters has failed because the spin degeneracy has
not being broken. The band characters are used to determine the position of
the charge neutrality point and the band indices, and if this fails this can
lead to subsequent problems, especially with post-processing functions that
rely on the band indices, like density of states. We avoid these problems by
splitting the degeneracy with split 0.01. The value $0.01$ (in meV) hardly
ever needs to be changed.

        kdotpy 2d 8o noax msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 7 10 zres 0.25 k -0.6 0.6 / 120 kphi 45 split 0.01

  6. 6.

With erange -80 0 we zoom in on the energy range of our interest. Moreover,
let us choose a colour scale representing the orbital degree of freedom of the
eigenstates. We set this by obs orbitalrgb. We also include the figure legend
with legend. We add subband character labels with char.

        kdotpy 2d 8o noax msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 7 10 zres 0.25 k -0.6 0.6 / 120 kphi 45 split 0.01 erange -80 0
    obs orbitalrgb legend char

  7. 7.

We choose a label for the filenames with out -7nm and the target folder with
outdir data-qw.

        kdotpy 2d 8o noax msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 7 10 zres 0.25 k -0.6 0.6 / 120 kphi 45 split 0.01 erange -80 0
    obs orbitalrgb legend char out -7nm outdir data-qw

Note that if outdir is omitted, the files will end up in the subfolder data if
it exists, and in the current folder otherwise.

  8. 8.

Finally, we add the post-processing option extrema, which will gives us useful
information about the extrema.

        kdotpy 2d 8o noax msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 7 10 zres 0.25 k -0.6 0.6 / 120 kphi 45 split 0.01 erange -80 0
    obs orbitalrgb legend char out -7nm outdir data-qw extrema

Other than the subprogram kdotpy 2d, which must be at the beginning of the
command line, the other command-line arguments may be added in any order.

####  5.2.3 Results and interpretation

Figure 19:  (a) Dispersion of a $7\,\mathrm{nm}$ HgTe quantum well as shown by
dispersion-7nm.pdf. The colours indicate orbital character, see the legend. At
$k=0$, the band character labels are shown, where ‘$\pm$’ indicates that the
states are approximately two-fold degenerate. (b) The properties of the
extrema in extrema-7nm.csv. The maxima of the band with bindex = -1 are
highlighted (highlighting added manually, for illustration). (c) In
dispersion-7nm.csv, we find the observable values for all eigenstates that
kdotpy has calculated. The state closest to the maximum at finite momentum has
been highlighted. We read off the expectation values for the observables
gamma6, gamma8h, and gamma8l to gain information on the orbital character at
this point.

After running this command, we will find the following files in the target
output folder data-qw:

  * •

`dispersion-7nm.pdf`: The dispersions $E(k)$ for all calculated eigenstates,
see Fig. 19(a). The colours visualize the values of the observable(s) chosen
with the option obs.

  * •

`dispersion-7nm.csv`: A csv file where all data points (momentum, energy) are
given with the values of the observables. The points are ordered by momentum.

  * •

`dispersion-7nm.byband.csv`: A csv file with the band energies for each curve
of `dispersion-7nm.pdf`, i.e., ordered by band.

  * •

`extrema-7nm.csv`: A csv file listing all the extremal values.

  * •

`output-7nm.xml`: The XML data file (see Sec. 3.10.2) that can be read again
by the subprograms `kdotpy merge` and `kdotpy compare` for the purpose of
replotting and comparing data sets and that contains a large amount of
metadata that can be used for diagnostics and reproducing the results later.
Discarding this file is strongly discouraged.

With these results, we can answer the questions above:

  * •

_Where is the finite-momentum maximum in momentum and energy?_ — In the plot,
we see that these maxima appear in the band with character $\mathrm{E}1\pm$ at
$k=0$. We find the following data in extrema.csv, see Fig. 19(b). We find
three maxima of the $\mathrm{E}1+$ band, one at $k=0$ (highlighted in blue)
and two at finite $k$ (highlighted in red), namely, at $k=\pm 0.463\
\mathrm{nm}^{-1}$ and $E=-40.1\,\mathrm{meV}$. Since we have calculated along
the (110) axis, with kphi 45, these points are located at $\mathbf{k}=\pm
0.463\,\mathrm{nm}^{-1}\times(\cos 45^{\circ},\sin
45^{\circ})=\pm(0.327,0.327)\,\mathrm{nm}^{-1}$.

  * •

_Do we have a direct or an indirect gap?_ — The maximum of the $\mathrm{E}1+$
band at $k=0$ lies at $E=-37.2\,\mathrm{meV}$. The maxima at finite $k$ are
lower, so the global maximum is at $k=0$. The value bindex = -1 confirms that
this is the highest valence band state. The $\mathrm{H}1-$ state with bindex =
1 is the lowest conduction band state. It only has one minimum, at $k=0$ and
$E=-19.7\,\mathrm{meV}$. We find that the gap is direct (at $k=0$) and its
size is $\lvert\Delta\rvert=17.5\,\mathrm{meV}$.

  * •

_What is the orbital character at the finite-momentum maximum?_ — We extract
this information from dispersion-7nm.csv that contains the values of the
observables. The maximum is located at $k=\pm 0.463\,\mathrm{nm}^{-1}$, so we
check the states at the nearest value $k=0.46\,\mathrm{nm}^{-1}$ (as a first
approximation). The band with bindex = -1 is the correct one. We could also
identify it by energy. The orbital character is determined by the probability
densities in the $\lvert\Gamma_{6},\pm\frac{1}{2}\rangle$,
$\lvert\Gamma_{8},\pm\frac{3}{2}\rangle$, and
$\lvert\Gamma_{8},\pm\frac{1}{2}\rangle$ states, which is encoded by the
observables gamma6, gamma8h and gamma8l, respectively. We thus find $0.8\%$
$\Gamma_{6}$ (‘electron or s orbital’), $50.4\%$ $\Gamma_{8\mathrm{H}}$
(‘heavy hole’), and $48.7\%$ $\Gamma_{8\mathrm{L}}$ (‘light hole’). We confirm
in the dispersion plot that the colour at this position is approximately cyan,
a mixture of $50\%$ green and $50\%$ blue.

###  5.3 Example calculation: Landau levels

####  5.3.1 Introduction

We set up a Landau-level calculation for the same structure as in Sec. 5.2. We
aim to calculate the Landau levels up to $B=10\,\mathrm{T}$ and to answer the
following questions, centred around the properties of the ‘lowest Landau
levels’, i.e., those that border the band gap:

  * •

Which Landau levels are the ‘lowest’ ones?

  * •

What is the critical field $B_{c}$ where the inversion of the lowest LLs is
undone?

  * •

What orbital character do the states have at this crossing?

We set up the calculation step by step, reusing some of the steps of Sec. 5.2.

####  5.3.2 Setting it up step by step

  1. 1.

First, we determine the subprogram. The geometry is 2D, and we use the Landau
level formalism. We thus use

        kdotpy ll

  2. 2.

As before, we use the eight-orbital model. In order to speed up the
calculation, let us use the axial approximation. For kdotpy ll, this is
achieved by omitting noax.

        kdotpy ll 8o

  3. 3.

The layer stack parameters are the same as in Sec. 5.2.

        kdotpy ll 8o msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%\
    llayer 10 7 10 zres 0.25

  4. 4.

Next, we specify the grid for the magnetic field values. Let is choose 100
points with the upper limit $10\,\mathrm{T}$. We could choose evenly spaced
points (every $0.1\,\mathrm{T}$ with b 0 10 / 100, but for LL calculations,
quadratic stepping is recommended. This is achieved by the ‘double slash’
notation b 0 10 // 100

        kdotpy ll 8o msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 7 10 zres 0.25 b 0 10 // 100

  5. 5.

Like in Sec. 5.2, we zoom in at the energy range erange -80 0 and we use split
0.01. This time, we also enter the diagonalization parameters explicitly. We
choose the maximum LL index to be 20 using nll 20. For each LL index, we aim
to get approximately 12 (subband) states. For this we need about 240
eigenvalues: neig 240. We also put the target energy for the diagonalization
somewhere near the gap, which we know is at approximately $-30\,\mathrm{meV}$.
Note, however, that if we put the target energy at that value, the
$\mathrm{E}2\pm$ subbands might be out of range. These subbands are required
for the determination of the charge neutrality point to work properly, so we
put the target energy higher up, i.e., we put targetenergy 0. (For the same
reason, we aim for 12 subband states instead of the minimally required number
of 8.)

        kdotpy ll 8o msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 7 10 zres 0.25 b 0 10 // 100 split 0.01 erange -80 0 nll 20
    neig 240 targetenergy 0

  6. 6.

For colouring the states, let us choose a combination of Landau level index
and total angular momentum. We set this by obs llindex.jz. Again, we include
the figure legend with legend and the band characters with char.

        kdotpy ll 8o msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 7 10 zres 0.25 b 0 10 // 100 split 0.01 erange -80 0 nll 20
    neig 240 targetenergy 0 obs llindex.jz legend char

  7. 7.

We choose a label for the filenames `out -7nm-landau` and the folder where the
files should go `outdir data-landau`.

        kdotpy ll 8o msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 7 10 zres 0.25 b 0 10 // 100 split 0.01 erange -80 0 nll 20
    neig 240 targetenergy 0 obs llindex.jz legend char out -7nm-landau
    outdir data-landau

####  5.3.3 Results and interpretation

Figure 20:  (a) The Landau level spectrum, bdependence-7nm-landau.pdf. The
colours indicate the Landau level index $n$ and the shade (brighter/darker)
indicates the sign of $\langle J_{z}\rangle$, see the colour legend. (b)
Screenshot of bdependence-7nm-landau.byband.csv which provides the energies of
the ‘lowest’ Landau levels $\lvert n,b\rangle$ with $(n,b)=(-2,1)$ and
$(0,-1)$ as a function of the magnetic field $B$ (column A, labelled bz). (c)
Screenshot of bdependence-7nm-landau.csv which contains the expectation values
for all relevant observables for all states. In (b) and (c), some columns were
hidden and the coloured highlights were added manually to lay emphasis on the
lowest Landau level states near the crossing.

After running this command, we get the following files in the folder `data-
landau`:

  * •

`bdependence-7nm-landau.pdf`: A plot of the Landau level spectrum, i.e., the
eigenenergies of the states $\lvert n,b\rangle$ as function of the magnetic
field, see Fig. 20(a).

  * •

`bdependence-7nm-landau.csv`: A csv file where all data points (magnetic
field, energy) are given with the values of the observables. The points are
ordered by magnetic field value.

  * •

`bdependence-7nm-landau.byband.csv`: A csv file with the band energies for
each curve, i.e., ordered by LL band.

  * •

`output-7nm-landau.xml`: The XML file with all parameters and data.

We can extract sufficient information to answer the questions above:

  * •

_Which Landau levels are the ‘lowest’ ones?_ — Let us examine the plot of Fig.
20(a). The lowest LL index formally is $-2$, but the actual lowest index
differs by subband character. Let us focus on the $\mathrm{H}1\pm$ and
$\mathrm{E}1\pm$ subbands, as these are the most relevant for the physics. The
‘$\pm$’ indicates the expectation value $\langle J_{z}\rangle$ at zero
momentum and zero magnetic field. This degree of freedom is encoded in the
colour shading, see the legend: dark colours represent the $+$ states
($\langle J_{z}\rangle>0$) and bright colours the $-$ states ($\langle
J_{z}\rangle<0$).

    * –

For $\mathrm{H}1+$, dark colours emanating from ‘$\mathrm{H}1\pm$’, we find
that the lowest index is $1$ (red colour); $0$, $-1$, and $-2$ are missing.

    * –

For $\mathrm{H}1-$, bright colours emanating from ‘$\mathrm{H}1\pm$’, the
baby-blue LL can be identified as index $-2$, so all indices are present for
$\mathrm{H}1-$.

    * –

For $\mathrm{E}1+$, dark colours emanating from ‘$\mathrm{E}1\pm$’, the lowest
index is $0$ (dark green).

    * –

For $\mathrm{E}1-$, bright colours emanating from ‘$\mathrm{E}1\pm$’, the
lowest index is $-1$ (bright orange green).

These observations can be confirmed from the data in the csv files. The lowest
LLs near the gap are thus ($\mathrm{E}1+,n=0$) and ($\mathrm{H}1-,n=-2$).

  * •

_What is the critical field $B_{c}$ where the inversion of the lowest LLs is
undone?_— In principle we could estimate this from Fig. 20(a), but we use one
of the csv files in order to answer this question more precisely: In
`bdependence-7nm-landau.byband.csv`, see Fig. 20(b), the bands are ordered in
columns labelled by $(n,b)$ in the first row, where $n$ is the LL index and
$b$ is the band index. The band indices are counted for each LL separately,
which makes intuition a bit tricky for the lowest LLs ($n=-2,-1,0$), but it is
guaranteed that the LLs at the crossing have $(n,b)=(0,-1)$ and $(-2,1)$
respectively. The respective columns are highlighted. We find the crossing
between $4.624\,\mathrm{T}$ and $4.761\,\mathrm{T}$. Linear interpolation may
be used as to obtain a more accurate value.

  * •

_What orbital character do the states have at this crossing?_ — From
`bdependence-7nm-landau.csv`, depicted in Fig. 20(c), we extract the orbital
character from the respective states at $B=4.624\,\mathrm{T}$ (which lies
closer to the actual crossing than $B=4.761\,\mathrm{T}$. The desired
information is encoded in the observables gamma6, gamma8h, and gamma8l for the
two respective states there. The lowest LL coming from the $\mathrm{E}1+$
band, $(n,b)=(0,-1)$, is a mixture of $49.4\%$ $\Gamma_{6}$, $46.7\%$
$\Gamma_{8\mathrm{L}}$, and $3.3\%$ $\Gamma_{8\mathrm{H}}$. We thus find that
it is predominantly a $J_{z}=\pm\frac{1}{2}$ state. The one from the
$\mathrm{H}1-$ band is purely $\Gamma_{8\mathrm{H}}$, and has $\langle
J_{z}\rangle=-\frac{3}{2}$ exactly.

###  5.4 Example calculation: Magneto-transport

####  5.4.1 Introduction

In the previous example calculation, we have set up a basic LL calculation.
With some postprocessing in kdotpy, we are able to extract much more useful
information and simulate magneto-transport experiments. The key ingredients
here are density of states (DOS) and Chern numbers (or Berry curvature). We
will set up an example simulation, where we will address the following
questions:

  * •

Which Landau levels are filled at the constant carrier density $2\times
10^{11}\,\mathrm{cm}^{-2}$?

  * •

Do we see plateaus in the Hall resistance $R_{xy}$?

  * •

Are there Shubnikov-de Haas (SdH) oscillations at low magnetic fields?

  * •

Can we simulate a Landau fan such that it resembles experimental data?

####  5.4.2 Setting up the simulation step by step

We take the the command line from Sec. 5.3 as a starting point:

    kdotpy ll 8o msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 7 10 zres 0.25 b 0 10 // 100 split 0.01 erange -80 0 nll 20
    neig 240 targetenergy 0 obs llindex.jz legend char out -7nm-landau
    outdir data-landau

We modify and add command line arguments as follows:

  1. 1.

We change some parameters slightly for this calculation. We choose erange -60
40 instead of erange -80 0, in order to give a better view of the data. We
increase the number of eigenvalues from neig 240 to neig 300, so that we
obtain a sufficient number of states to view this energy range. We also change
the output filenames and output folder.

        kdotpy ll 8o msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 7 10 zres 0.25 b 0 10 // 100 split 0.01 erange -60 40 nll 20
    neig 300 targetenergy 0 obs llindex.jz legend char out -7nm-hall
    outdir data-hall

  2. 2.

In order to calculate the density of states, we add the options dos and
localdos for calculating the total and local density of states, respectively.
For the DOS, we use Gaussian broadening as in Ref. [3], with a broadening
width (standard deviation) $\sigma(B)=\sigma_{1}\sqrt{B\,[\mathrm{T}]}$ with
$\sigma_{1}=0.5\,\mathrm{meV}$ (see also Appendix B.6). The appropriate
command is broadening gauss 0.5 sqrt, but this may be shortened to broadening
0.5; for kdotpy ll, gauss and sqrt are default settings for the broadening.

        kdotpy ll 8o msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 7 10 zres 0.25 b 0 10 // 100 split 0.01 erange -60 40 nll 20
    neig 300 targetenergy 0 obs llindex.jz legend char out -7nm-hall
    outdir data-hall dos localdos broadening 0.5

  3. 3.

For calculating the Hall conductance, we need to calculate the Chern numbers
and set its broadening correctly. We add chern to calculate the Chern numbers.

        kdotpy ll 8o msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 7 10 zres 0.25 b 0 10 // 100 split 0.01 erange -60 40 nll 20
    neig 300 targetenergy 0 obs llindex.jz legend char out -7nm-hall
    outdir data-hall chern dos localdos broadening 0.5

  4. 4.

We replace `broadening 0.5` by `broadening 0.5 10%` to set the broadening of
the Hall conductance as function of energy to $10\%$ of that of the carrier
density (integrated density of states).

        kdotpy ll 8o msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 7 10 zres 0.25 b 0 10 // 100 split 0.01 erange -60 40 nll 20
    neig 300 targetenergy 0 obs llindex.jz legend char out -7nm-hall
    outdir data-hall chern dos localdos broadening 0.5 10%

  5. 5.

For convenience, one may also use the command shortcut hall, which is
equivalent to chern dos localdos broadening 0.5 10%.

        kdotpy ll 8o msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 7 10 zres 0.25 b 0 10 // 100 split 0.01 erange -60 40 nll 20
    neig 300 targetenergy 0 obs llindex.jz legend char out -7nm-hall
    outdir data-hall hall

This command line is equivalent to the one in step 4.

For this calculation, the following configuration options are relevant:

  * •

`dos_unit`: The plots in Fig. 21 have been produced with `dos_unit=cm`, so
that densities are expressed in units of $\mathrm{cm}^{-2}$ with an
appropriate power of ten.

  * •

`dos_energy_points`: The default energy resolution used for DOS calculations
(etc.) is `dos_energy_points=1000`, which is usually is a good compromise
between accuracy and computation time. If the value is insufficiently large,
one may see visible artifacts, like jagged equal-density lines where one would
expect smooth curves. If that is the case, one should increase this value.

  * •

`berry_ll_simulate`: By default (`berry_ll_simulate=false`), the plots related
to Hall conductance are extracted from the calculated Chern numbers (Berry
curvature), stored as the observable chern. When the number of states, set by
command-line argument neig, is not sufficient, some values are inaccurate or
incorrect; see Sec. 3.5.3 for a detailed discussion. In bad cases, this
renders the data to be partially or completely unusable. One can mitigate this
problem by setting `berry_ll_simulate=true`. In this case, the output
functions use a simulated Chern number, which is exactly $1$ for all Landau-
level eigenstates and is stored as the observable chernsim. This value
coincides with the calculated Chern number (up to numerical error) for almost
all states, except for states the edges if the spectrum, which should usually
be mistrusted anyway. If `berry_ll_simulate=true`, the output file names will
be changed by insertion of sim or simul in order to be able to easily spot
whether the calculated or simulated Chern numbers have been used. It is
recommended to run kdotpy with both settings and to compare the results at
least once in order to get an intuition for the implications of this setting
and for the physics behind it.

The configuration values may be set permanently with

    kdotpy config ’dos_unit=cm;dos_energy_points=1000;berry_ll_simulate=true’

or used once by appending

    config ’dos_unit=cm;dos_energy_points=1000;berry_ll_simulate=true’

to the command line that starts kdotpy ll.

####  5.4.3 Results and interpretation

The above commands produce many pdf and csv files, and here we shall only
discuss the ones relevant for answering the questions above. In most cases,
the pdf is accompanied by a csv file which contains exactly the data shown in
the figure.

Figure 21:  (a) In the file bdependence-density.pdf, equal-density contours
are overlaid on top of the magnetic field depdendence plot, cf. Fig. 20(a).
The thick line near $-20\,\mathrm{meV}$ is the charge neutral point $n=0$.
Contours are drawn for each multiple of $1\times 10^{11}\,\mathrm{cm}^{-2}$
for $\lvert n\rvert\leq 15\times 10^{11}\,\mathrm{cm}^{-2}$. The medium-thick
lines are at odd multiples of $5\times 10^{11}\,\mathrm{cm}^{-2}$ and the very
thick lines are at multiples of $10\times 10^{11}\,\mathrm{cm}^{-2}$. (b) The
plot in rxy-constdens-7nm-hall.pdf shows the simulated Hall resistance
$R_{xy}\equiv 1/\sigma_{\mathrm{H}}$ as function of magnetic field $B$ for the
constant density $2\times 10^{11}\,\mathrm{cm}^{-2}$. The dashed line
indicates the classical Hall slope given by $R_{xy}=B/ne$. (c) The file dos-
constdens-7nm-hall.pdf visualizes the density of states as function of
magnetic field for a constant carrier density. (d) In dsigmah-dn-7nm-hall.pdf,
we plot the derivative $d\sigma_{\mathrm{H}}/dn$ of the Hall conductance with
respect to the carrier density $n$. This type of plot can easily be compared
to plots obtained from experiments which measure Hall conductance as function
of magnetic field and gate voltage. The figures (b) and (c) are single pages
in a multipage PDF unless cardens has been used in order to specify a single
density. For this figure, we have used the configuration settings dos_unit=cm,
dos_energy_points=1000, and berry_ll_simulate=false.

  * •

_Which Landau levels are filled at the constant carrier density $2\times
10^{11}\,\mathrm{cm}^{-2}$?_— The answer to this question of course depends on
the magnetic field. It is encoded in the file bdependence-density.pdf, see
Fig. 21(a), which is generated when kdotpy ll is used with the dos option.
Equal-density contours are overlaid onto the magnetic-field dependence plot,
(LL energies as function of magnetic field $B$). In Fig. 21(a), we find the
equal-density curve for $n=2\times 10^{11}\,\mathrm{cm}^{-2}$ to be at
approximately $20\,\mathrm{meV}$ for low magnetic fields. A different
visualization may be obtained by adding cardens 0.002 to the command line, in
which case only the contour line for the given carrier density
($0.002\,\mathrm{nm}^{-2}=2\times 10^{11}\,\mathrm{cm}^{-2}$) is shown.

  * •

_Do we see plateaus in the Hall resistance $R_{xy}$?_— The Hall resistance is
shown in the plot rxy-constdens-7nm-hall.pdf, see Fig. 21(b). We find plateaus
corresponding to the filling factors $\nu=1,2,\ldots$. At higher filling
fractions (lower magnetic fields), the plateau transitions are more washed out
due to the broadening. Here, the curve tends more towards the classical value
$R_{xy}=B/ne$, indicated by the dashed line. This type of visualization can be
compared directly to experimental results. For completeness, kdotpy also
generates the plot sigmah-constdens-7nm-hall.pdf, which shows the Hall
conductance $\sigma_{\mathrm{H}}$ as function of magnetic field for constant
carrier density. In both cases, these output files correspond to the density
specified by cardens, if this argument is provided on the command-line.
Otherwise, the output will be multi-page PDFs with multiple densities, namely
integer multiples of $1\times 10^{11}\,\mathrm{cm}^{-2}$.

  * •

_Are there Shubnikov-de Haas (SdH) oscillations at low magnetic fields?_ — The
constant-carrier-density curves in Fig. 21(c) clearly show oscillations with
magnetic field. The plot shows density of states, which is thought to
correlate with the longitudinal resistance $R_{xx}$. The red markers at the
bottom indicate the values of $1/B$ equal to integer multiples of $e/hn$,
where $n$ is the carrier density. The minima of the density of states thus
align well with the multiples of $e/hn$ for $B\gtrsim 1\,\mathrm{T}$. At lower
fields, the resolution of $B$ values is insufficient to make any conclusive
statement. This may easily be mitigated by an additional calculation at lower
fields, for example with b 0 1 // 100.

  * •

_Can we simulate a Landau fan such that it resembles experimental data?_ — In
magneto-transport experiments, Landau fans are obtained by measuring $R_{xy}$
while sweeping the gate voltage for many values of the magnetic field. This
yields $R_{xy}$ as function of magnetic field $B$ and gate voltage
$V_{\mathrm{g}}$. In simulations, however, the natural quantity on the
‘vertical axis’ is energy. The DOS calculation in kdotpy bridges this gap: The
relation of carrier density $n$ as function of energy $E$ can be reversed as
to obtain a spectrum as function of density. The latter often correlates
approximately linearly with gate voltage $V$, so that comparison between the
two is physically meaningful. In Fig. 21(d), we visualize the plateau
transitions by taking the derivative $d\sigma_{\mathrm{H}}/dn$ in order to
highlight the plateau transitions. This method is also customarily applied to
experimental data.

###  5.5 Example calculation: Optical transitions

####  5.5.1 Introduction

Instead of extracting information on magneto-transport experiments in
postprocessing, we can also simulate magneto-optical transition spectra in
kdotpy. We will use an example simulation to answer the following questions:

  * •

How large is the band gap of the system?

  * •

What are the involved LL states for a specific transition?

  * •

What information regarding band ordering can we extract from the transition
spectrum?

####  5.5.2 Setting up the simulation step by step

Again, we take the the command line from Sec. 5.3 as a starting point:

    kdotpy ll 8o msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 7 10 zres 0.25 b 0 10 // 100 split 0.01 erange -80 0 nll 20
    neig 240 targetenergy 0 obs llindex.jz legend char out -7nm-landau
    outdir data-landau

We modify and add command line arguments as follows:

  1. 1.

As for the magneto-transport analysis, we slightly adjust some parameters. We
choose erange -80 50 for a better view of the data. We increase the number of
eigenvalues from neig 240 to neig 300 to get a sufficient number of states to
view in this energy range. We also change the output filenames and output
folder.

        kdotpy ll 8o msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 7 10 zres 0.25 b 0 10 // 100 split 0.01 erange -80 50 nll 20
    neig 300 targetenergy 0 obs llindex.jz legend char
    out -7nm-optical-transitions outdir data-optical-transitions

  2. 2.

In order to calculate all optical transitions, we add the option transitions.

        kdotpy ll 8o msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 7 10 zres 0.25 b 0 10 // 100 split 0.01 erange -80 50 nll 20
    neig 300 targetenergy 0 obs llindex.jz legend char
    out -7nm-optical-transitions outdir data-optical-transitions transitions

  3. 3.

It is also advisable to filter transitions by state occupancy, since a
transition can only be observed if the initial state is (partially) occupied
while the final state is (partially) unoccupied. Thus, we add the option
cardens 0, assuming a charge neutral sample. (There is also the option to
calculate filtered transitions for multiple carrier densities at once by using
the input `cardens # # / #`, analogous to the range input for b.)
Additionally, we also use the options dos and broadening 0.5. dos enables
calculation of the (electro-)chemical potential, which is used to determine
the occupation of states, while broadening applies a Gaussian square-root
broadening to the LL states (shortened, see Sec. 5.4). Omitting broadening
often leads to abrupt jumps in the (electro-)chemical potential, leading to
unphysical gaps along specific transition features. Note, that broadening only
is applied to DOS-related quantities, transition spectra will always be delta
peaks.

Per default, these options suppress the output of all possible transitions.

        kdotpy ll 8o msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 7 10 zres 0.25b 0 10 // 100 split 0.01 erange -80 50 nll 20
    neig 300 targetenergy 0 obs llindex.jz legend char
    out -7nm-optical-transitions outdir data-optical-transitions transitions
    cardens 0 dos broadening 0.5

  4. 4.

Last, we change some configuration values using config. With
`transitions_max_deltae=80` we set the upper limit of the energy axis in the
transition plot to $80\,\mathrm{meV}$, while `transitions_min_amplitude=3e11`
filters out all transitions with transition matrix elements smaller than the
given threshold, keeping only the most prominent transitions.

        kdotpy ll 8o msubst CdZnTe 4% mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 7 10 zres 0.25 b 0 10 // 100 split 0.01 erange -80 50 nll 20
    neig 300 targetenergy 0 obs llindex.jz legend char
    out -7nm-optical-transitions outdir data-optical-transitions transitions
    cardens 0 dos broadening 0.5
    config ’transitions_max_deltae=80;transitions_min_amplitude=3e11’

The quantity that will be plotted in the transitions spectrum can be changed
by using the configuration option `plot_transitions_quantity`. We will use the
default value rate, plotting the rate density as defined in Eq. (121). More
information on this topic can be found in Sec. 3.8.4.

####  5.5.3 Results and interpretation

Using the commands above produces additional output files. This includes
bdependence-transitions-7nm-optical-transitions.pdf, where all filtered
transitions are drawn as vertical lines into the LL fan chart, and
transitions-filtered-7nm-optical-transitions.pdf, together with its csv file,
showing the filtered optical transitions spectrum. We will address the
question given above by discussing the transitions spectrum in Fig. 22(a).

Figure 22:  (a) The filtered transitions spectrum, transitions-filtered-7nm-
optical-transitions.pdf. The strength of each individual transition is
indicated by the size and colour of its respective symbol. The blue arrow
points towards the transition feature that extrapolates to the band gap energy
at small magnetic fields, while the black arrow indicates the feature for
which we want to find the involved LL states. Both were added in post-
production and are normally not include in the plot. The screenshots on the
right are taken from the transitions-filtered-7nm-optical-transitions.csv
file. In (b) the data is filtered by the band indices B1 and B2, taking the
values $-1$ and $1$, respectively. For (c) data was filtered for magnetic
fields bx $\geq 6\,\mathrm{T}$. The Landau level spectrum corresponding to
this figure is the one shown in Fig. 21.

  * •

_How large is the band gap of the system?_ — The interband transition between
the highest valence and the lowest conduction subband at low magnetic fields
extrapolates to the energetic value of the direct band gap, see e.g. Fig.
21(a). We suspect, that the transition feature marked with the blue arrow in
Fig. 22(a) is this specific transition. To confirm this, we use the
transitions-filtered-7nm-optical-transitions.csv file and filter the initial
state B1 by highest valence subband and the final state B2 by lowest
conduction subband. Per definition, these subbands, respectively, have the
band index $-1$ and $1$. The blue and green columns in Fig. 22(b) show these
indices, while the red column shows the transition energies. The lowest
calculated magnetic field bx is $0.001\,\mathrm{T}$ and the corresponding
transition energies deltaE are $>17.5\,\mathrm{meV}$, confirming that this
feature extrapolates to the same band gap as discussed in Sec. 5.2. For zero
magnetic field all LLs in a subband are degenerate, thus, all transition
energies would be identical (up to the artificial offset added by split).

  * •

_What are the involved LL states for a specific transition?_ — We analyse the
transition marked with the black arrow in Fig. 22(a) regarding which subbands
and LL indices are involved. We use the same output file transitions-
filtered-7nm-optical-transitions.csv and filter the magnetic field values bx
to be $>6\,\mathrm{T}$, see Fig. 22(c), so that we can clearly distinguish
this transition from other transitions at higher energies than shown in the
plot. The marked columns are the initial (blue) and the final (green) LL
(light) and band indices (dark). As can be seen, this transition is an
intraband transition inside the subband with band index $1$, which is the
inverted subband H1. The respective LL indices of initial and final state are
$-2$ and $-1$, hence the polarization of this transition corresponds to
$O_{+}$ [see Sec. 3.8].

  * •

_What information regarding band ordering can we extract from the transition
spectrum?_ — Continuing the analysis for the same transition, we discuss why
it can only be observed at magnetic fields $\gtrapprox 5\,\mathrm{T}$. For
that, we use the LL fan chart shown in the previous tutorial section Fig. 21,
where the same sample was simulated. The relevant (electro-)chemical potential
for the current case is the thick line near $-20\,\mathrm{meV}$, the charge
neutral point. At magnetic fields slightly below $5\,\mathrm{T}$ the lowest LL
band of H1 and E1 cross, reverting to a trivial band order, while
simultaneously changing from completely unoccupied to fully occupied in the
case of the H1 state (vice versa for the E1 state). This enables the
corresponding intraband transition in the H1 subband, which we see in the
transition spectrum. Consequently, we can indirectly identify the critical
magnetic field, where the lowest LL states revert to trivial band order in the
magneto-optical transition spectrum. The abruptness of the crossing is a
consequence of the axial approximation, which we used here for illustrative
purposes. In a more realistic picture, the non-axial and bulk-inversion
asymmetric terms lead to an anticrossing instead.

###  5.6 Example calculation: Dispersions and wave functions of a 3D TI

####  5.6.1 Introduction

In thick layers of a band inverted material like HgTe, so called _three-
dimensional topological insulators_ , the transport properties are dominated
by topological surface states. At an interface between materials with inverted
and normally ordered bands, the band inversion induces eigenstates with a
linear (two-dimensional Dirac) dispersion, confined to the interface [43]. The
physics is markedly different from those of narrow quantum wells (see Sec.
5.2), where instead the strong confinement in the $z$ direction leads to the
inversion of subbands.

In this usage example, we investigate the key features of a 3D topological
insulator by analyzing the dispersions and wave functions. We try to answer
the following questions:

  * •

Can we identify the Dirac point and the surface states in the dispersion?

  * •

How well are the surface states confined?

  * •

How much is the overlap between the surface states?

  * •

What is the orbital content of the wave functions?

####  5.6.2 Setting it up step by step

  1. 1.

Despite the system being called ‘3D topological insulator’, the appropriate
geometry is still 2D, because there is confinement in the $z$ direction. We
take the previous example, Sec. 5.2, and adjust the layer stack parameters
msubst and llayer. We also use a different momentum and energy range. We raise
the number of eigenvectors with neig 100, because there are more subbands due
to weaker confinement in the $z$ direction. We choose the observable z in
order to see whether we can locate the surface states at the top and bottom
interface.

        kdotpy 2d 8o noax msubst CdTe mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 70 10 zres 0.25 k -0.3 0.3 / 120 kphi 0 split 0.01 erange -100 40
    neig 100 obs z legend char out -70nm outdir data-3dti

  2. 2.

We obtain the dispersion plot of Fig. 23(a). At $E\approx-90\,\mathrm{meV}$,
we observe the linear dispersions characteristic of a Dirac point. Between
valence band and conduction bands (negative and positive energies,
respectively), there are states bridging the gap. The gray colour indicates
that the expectation values $\langle z\rangle\approx 0$ for all states, which
seemingly suggests absence of surface character. In order to diagnose this
counterintuitive result, let us plot plot the wave functions at
$k=0.14\,\mathrm{nm}^{-1}$ using plotwf separate 0.14.

        kdotpy 2d 8o noax msubst CdTe mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 70 10 zres 0.25 k -0.3 0.3 / 120 kphi 0 split 0.01 erange -100 40
    neig 100 obs z legend char out -70nm outdir data-3dti plotwf separate 0.14

  3. 3.

We obtain two files with wave functions, i.e., at $k=-0.14\,\mathrm{nm}^{-1}$
and at $k=0.14\,\mathrm{nm}^{-1}$. In the latter file, let us look for one of
the states at $E\approx 13.5\,\mathrm{meV}$, see Fig. 23(b). The wave function
contains symmetric and antisymmetric components for different orbitals. In
order to break the degeneracy due to the mirror symmetry in $z$ direction, let
us add a weak potential difference between top and bottom surface with vinner
0.1. This argument leads to an electrostatic potential consistent with a
constant displacement field $\mathbf{D}=\epsilon(z)\mathbf{E}(z)$, such that
the potential difference between top and bottom surface is
$0.1\,\mathrm{meV}$.

        kdotpy 2d 8o noax msubst CdTe mlayer HgCdTe 68% HgTe HgCdTe 68%
    llayer 10 70 10 zres 0.25 k -0.3 0.3 / 120 kphi 0 split 0.01 erange -100 40
    neig 100 obs z legend char out -70nm outdir data-3dti plotwf separate 0.14
    vinner 0.1

####  5.6.3 Results and interpretation

Figure 23:  (a) Dispersion of a $70\,\mathrm{nm}$ thick HgTe layer on a CdTe
substrate, dispersion-70nm.pdf. The colours indicate the expectation value
$\langle z\rangle$. Due to mirror symmetry $z\to-z$, $\langle z\rangle=0$ for
all states. (b) A wave function at $k=0.14\,\mathrm{nm}^{-1}$ and $E\approx
13.5\,\mathrm{meV}$. The solid curves indicate the real parts, the dashed
curves the imaginary parts for each orbital component. (c) The dispersion with
a degeneracy breaking potential (argument vinner 0.1). The surface character
of the states in the gap and near the Dirac point is indicated by red and blue
colour. States with bulk character are gray. (d) One of the surface states at
$k=0.14\,\mathrm{nm}^{-1}$ and $E\approx 13.5\,\mathrm{meV}$. Clearly, the
degeneracy is resolved. The angles on the right-hand side indicate the complex
arguments in each orbital component.

After applying vinner 0.1 in order to break the mirror symmetry, we obtain
`dispersion-70nm.pdf` and `wfs-70nm_0.140_0.pdf` as shown in Fig. 23(c) and
(d). With these plots and the corresponding csv files, we can answer the
questions above:

  * •

_Can we identify the Dirac point and the surface states in the dispersion?_ —
In Fig. 23(c), we find two pairs of linearly dispersing states, crossing at
$E\approx-89\,\mathrm{meV}$ at zero momentum. The subband labels are
$\mathrm{E}1\pm$ and $\mathrm{L}1\pm$. The states emanating from the Dirac
point have a distinct surface character is indicated by the red and blue
colour (where the latter is poorly visible because the dispersions almost
coincide).

These four Dirac states can be described by a massless two-dimensional Dirac
Hamiltonian. We note that the Dirac model for the surface states is distinct
from the BHZ model [39], which has a different set of basis states, namely the
$\mathrm{E}1\pm$ and $\mathrm{H}1\pm$ subbands; see also the example in Sec.
5.2.

The surface states seem to disappear into the bulk valence band for energies
$-70\,\mathrm{meV}\leq E\leq-10\,\mathrm{meV}$ approximately, due to
hybridization with the bulk valence band. In the bulk gap
($0\,\mathrm{meV}\leq E\leq 20\,\mathrm{meV}$ approximately), the surface
states cross the gap from valence to conduction band.

  * •

_How well are the surface states confined?_ — The level of confinement can be
extracted from the wave functions. In this example, we investigate one of the
surface states in the bulk gap, at $k=0.14\,\mathrm{nm}^{-1}$ and $E\approx
13.5\,\mathrm{meV}$. From the wave functions shown in Fig. 23(d), we estimate
a characteristic ‘size’ of about $10$–$20\,\mathrm{nm}$.

For a more rigorous value, we can take the file `{wfs-70nm_0.140_0.1.csv`
(depending on the configuration settings, this might be inside a tar or zip
archive) and calculate
$\lvert\psi\rvert^{2}(z)=\sum_{i}\lvert\psi_{i}\rvert^{2}(z)$, where $i$ are
the orbitals. We determine the characteristic length $\xi$ by fitting
$\psi(z)\sim\mathrm{e}^{-z/\xi}$. From approximating
$\xi^{-1}\approx\partial_{z}\ln(\lvert\psi\rvert^{2})$ by a discrete
derivative around $z=0$, we find $\xi\approx 7\,\mathrm{nm}$. Note that the
apparent discrepancy with the estimate below comes from the difference between
$\psi_{i}(z)$ and $\lvert\psi\rvert^{2}(z)$.

  * •

_How much is the overlap between the surface states?_ — In order to estimate
the overlap between the surface states at top and bottom surface, we calculate
the integral of $\lvert\psi\rvert(z)$ for $z>0$, for one of the surface
states. We find $\int_{0}^{z_{\mathrm{max}}}\lvert\psi\rvert^{2}(z)\approx
0.01$ by summing up the appropriate values in `wfs-70nm_0.140_0.1.csv`. We
finally point out that the overlap and the confinement length are momentum
dependent, so that for a complete picture of the physics, this analysis has to
be repeated for other momenta.

  * •

_What is the orbital content of the wave functions?_ — The data in Fig. 23(d)
visualizes the wave function as $\psi(z)\sum_{i}\psi_{i}(z)\lvert i\rangle$,
where $\lvert i\rangle$ are the orbital basis functions (see Eq. 10). The
curves being solid or dashed indicate the $\psi_{i}(z)$ being real or
imaginary. In many cases, the functions $\psi_{i}(z)$ can be factorized as a
complex constant $\mathrm{e}^{\mathrm{i}\phi_{i}}$ times a real function
$f_{i}(z)$. The angles in degrees in the right hand side of Fig. 23(d)
indicate the complex arguments $\psi_{i}$. Thus, the wave function for this
state can be written as

$\psi(z)=f_{1,2}(z)(-\lvert\mathbf{1}\rangle-\mathrm{i}\lvert\mathbf{2}\rangle)+f_{3,6}(z)(\lvert\mathbf{3}\rangle-\mathrm{i}\lvert\mathbf{6}\rangle)+f_{4,5}(z)(\mathrm{i}\lvert\mathbf{4}\rangle-\lvert\mathbf{5}\rangle)+f_{7,8}(z)(\mathrm{i}\lvert\mathbf{7}\rangle+\lvert\mathbf{8}\rangle)$
(136)

in terms of four real functions $f_{i,i^{\prime}}(z)$, where components $i$,
$i^{\prime}$ with identical shapes of the envelope function have been taken
together pairwise.

###  5.7 Auxiliary scripts

####  5.7.1 kdotpy merge

The subprogram kdotpy merge imports multiple XML files and merges the result
into a single data set and plots the result in a single dispersion or
magnetic-field dependence plot. (Further postprocessing is not supported in
version 1.0.0 of kdotpy.) The merging is done ‘horizontally’, i.e., data sets
with different values of $k$ or $B$) as well as ‘vertically’, i.e., with data
at the same values of $k$ and $B$, but for example with respect to different
target energies. The command-line syntax may be as simple as

    kdotpy merge data/output.1.xml data/output.2.xml

Usually, additional options are provided, for example setting the output
filenames and directory.

    kdotpy merge bandalign -25 2 out .merged
    outdir data-merged -- data/output.1.xml data/output.2.xml

The `--` (double hyphen) symbol separates the input file names from the other
options. It may be omitted if there is no potential confusion between the
other options and the input file names. The command line option bandalign can
be used to do or redo band alignment on the combined data set. Manually
assisted realignment as described in Sec. 3.6.4 is also supported.

For data integrity, kdotpy merge compares the physical parameters associated
with the source data files (stored in the <parameter> tag). If values differ,
it will show a warning as to inform the user that the data sets might not fit
together. There is also a check on whether the grids of the data files are
compatible, i.e., whether the combination of the grids can also be represented
as a VectorGrid.

One of the important use cases of kdotpy merge is the ability to split up
computationally demanding computations into more manageable pieces. For
example, strip calculation with kdotpy 1d (which has a high RAM footprint) can
be split up across multiple nodes of a high-performance computing cluster.
Another boon of kdotpy merge is that incomplete data sets can be completed
with just the missing data and using the already existing data to save
calculation time. This situation may occur, for example, when one out of many
cluster job fails, or when in hindsight the value of $n_{\mathrm{eig}}$ was
set too low to cover the desired certain energy interval. In Refs. [4] and
[6], we have used this strategy for the strip calculations with kdotpy 1d,
(with kdotpy versions v0.42 and v0.64, respectively 666The actual versions
used were development versions preceding v0.42 (for Ref. [4]) and v0.64 (for
Ref. [6]).. For the dispersion figures, $3$–$10$ data sets were merged.

####  5.7.2 kdotpy compare

The subprogram kdotpy compare is similar to kdotpy merge, but it can take an
arbitrary number of data files and merge them into multiple data sets, which
can then be compared in a single plot. The command-line syntax is similar to
kdotpy merge, however data sets must be separated using vs, as in the
following example,

    kdotpy compare data/output.1a.xml data/output.1b.xml vs data/output.2.xml

In this example, the first data set is constructed by merging the data from
data/output.1a.xml and data/output.1b.xml, and the second data set is the
single file data/output.2.xml. The two data sets are visualized with different
markers and/or colours. It is also possible to enter more than two data sets.
The syntax for extra options is similar to that of kdotpy merge.

kdotpy compare only handles one-dimensional dispersions and magnetic field
dependence. It is primarily used for quick comparisons between two related
systems, for example to study the effect of enabling or disabling a term in
the Hamiltonian. It is also useful as a debugging tool.

####  5.7.3 kdotpy batch

The auxiliary subprogram kdotpy batch is a tool for running any of the
calculation subprograms in batch mode, where many calculations are run in
sequence, but with a variation in one or more parameters. The philosophy is
similar to a shell (bash) script, but kdotpy batch may have a slightly more
intuitive way of iterating over the parameters and it can handle runs in
parallel.

The behaviour is best illustrated from an example command line,

    kdotpy batch [opts] @x 5.0 9.0 / 4 @y [CdTe, HgTe, HgCdTe 68%] do kdotpy 2d
    8o noax llayer 10 @x 10 mater HgCdTe 68% @y HgCdTe 68% out -@0of@@ ...

This kdotpy batch command runs iteratively the kdotpy subprogram after the
argument do. Here, kdotpy 2d is run fifteen times, with `@x` replaced by the
values 5.0, 6.0, 7.0, 8.0, and 9.0, and `@y` replaced by CdTe, HgTe, and
HgCdTe 68%. The command syntax for ranges is the same as for the grid
components, see Sec. 3.2.3. The special notations `@0` and `@@` denote a
counter and the total value respectively, i.e., `@0of@@` yields 1of15, 2of15,
etc. All ‘@ variables’ defined before do may appear after do an arbitrary
number of times. Curly brackets may be used (e.g., `@{x}`) if variables appear
inside a longer string.

Multiple processes may run in parallel with the cpu and proc arguments, to be
specified before do. These options also take into account the number of CPUs
of each process, if it is given by the argument cpu after do. Another useful
argument is dryrun, which prints the command lines that result from the ‘`@`
substitutions’, but does not execute the tasks.

In a batch calculation, the console output of the calculation scripts is
written to the filesstdout.1.txt, stdout.2.txt, etc. and stderr.1.txt,
stderr.2.txt, etc., where the extension can be adjusted with the configuration
values `batch_stdout_extension` and `batch_stderr_extension`. The kdotpy batch
subprogram itself writes to the standard stdout and stderr streams. At the end
of the batch run, kdotpy batch provides a status message that summarizes how
many tasks were completed successfully.

####  5.7.4 kdotpy config

Whereas the user may edit the configuration file ~/.kdotpy/kdotpyrc directly,
we also provide a command line tool kdotpy config, which can do the following
simple operations on the configuration:

  * •

kdotpy config list or kdotpy config show: Show all non-default configuration
values.

  * •

kdotpy config all: Show all configuration values.

  * •

kdotpy config help _key_ : Show information on _key_ (taken from the help
file; multiple _key_ arguments possible).

  * •

kdotpy config set _key_ =_value_ : Set _key_ to a new value (multiple _key_
=_value_ pairs possible).

  * •

kdotpy config reset _key_ : Set _key_ to its default value (multiple _key_
arguments possible).

  * •

kdotpy config file: Print the full path of the configuration file.

  * •

kdotpy config edit: Open the configuration file in an editor. (Use environment
variable VISUAL or EDITOR to set your command line editor.)

  * •

`kdotpy config fig_lmargin=10 fig_rmargin=10 fig_hsize` (no operation
specified): Set the values for `fig_lmargin` and `fig_rmargin`, show all three
values.

Multiple _key_ and/or _key_ =_value_ can also be combined with a semicolon ;,
for example

    kdotpy config ’fig_lmargin=10;fig_rmargin=10;fig_hsize’

(In bash one must use single quotes because of the special meaning of ;.) The
syntax is thus analogous to the config argument for the calculation
subprograms, for example kdotpy 2d.

####  5.7.5 kdotpy test

The package includes a set of standardized tests that can be run with the
command kdotpy test. The purpose of these tests is to catch errors during code
development and to catch performance issues. Each standardized test is defined
in kdotpy-test.py by a command line. Success or failure is determined by the
exit code. The output (by default written into the subdirectory test relative
to the current working directory) is not checked for errors and has to be
inspected manually, if desired. After all tests have been attempted, a summary
is printed to the terminal indicating success or failure and the time it took
to run.

Currently, in kdotpy version 1.0, there are 19 standardized tests, labelled
with the following test ids: `1d`, `2d_cartesian`, `2d_offset`, `2d_orient`,
`2d_polar`, `2d_qw`, `2d_qw_2`, `2d_qw_bia`, `2d_selfcon`, `batch`, `bulk`,
`bulk_3d`, `bulk_ll`, `ll_axial`, `ll_bia`, `ll_legacy`, `compare_2d`,
`compare_ll`, and `merge`. The full test suite may be run by typing

    kdotpy test

To run only a selection (of one or more tests), one may add the test ids as
arguments, for example

    kdotpy test 2d_qw ll_axial

to run the specified ones only. The test suite may be interrupted with Ctrl-C.

The test suite supports the following extra commands:

  * •

kdotpy test list: Shows all test ids on the terminal.

  * •

kdotpy test showcmd _testids_ : Shows the command lines without running the
tests. The command lines may be copied, modified, and run manually for the
purpose of testing and debugging.

  * •

kdotpy test verbose _testids_ : Runs the tests in verbose mode, by adding the
argument verbose to their command lines.

  * •

kdotpy test python3.9 _testids_ : Uses a specific Python command (in this
example python3.9) for running the test commands. It is recommended to use a
virtual environment for this purpose.

In these commands, _testids_ may be any sequence of test ids, or it may be
omitted to apply the command to all 19 tests.

The tests are defined in kdotpy-test.py in a way that is compatible with the
pytest package for automated testing. We have included automated testing with
pytest into the CI/CD workflow in our Gitlab project.

##  6 Conclusion and outlook

This publication marks the public release of kdotpy as an open source software
project. We release it with the expectation that kdotpy is beneficial to the
research community. Conversely, the open source nature encourages input from
the community as to improve and extend the project even further. It is our
hope that interaction with the community will lead to new insights and new
directions for the future development of kdotpy. We encourage our users to
actively participate in the discussions in our Gitlab repository [42] or via
other channels.

At the time of writing, we are already in the process of designing new
features for future releases. Improved strain handling is scheduled to be
included in an upcoming release in the near future. The treatment of
electrostatics is currently undergoing continuous improvement, where in
particular we are reconsidering the boundary conditions appropriate for
several experimental setups. We will also add the necessary infrastructure for
treatment of dopants in the material. We also expect to gain more experience
with materials other than the built-in ones. In general, the future directions
of kdotpy will also be influenced by input from new experimental results from
transport and spectroscopy and by feedback and suggestions from our users.

## Acknowledgements

We thank Domenico Di Sante, Giorgio Sangiovanni, Björn Trauzettel, Florian
Goth, and Fakher Assaad for feedback and support at various stages of the
project.

##### Author contributions

W.B. designed the software package. W.B., F.B., C.B., and M.H. developed the
source code, with W.B. acting as lead developer. W.B., F.B., C.B., L.B., M.H.,
and M.S. performed code review and benchmarks. W.B., C.F., S.S., L.-X.W.,
H.B., and L.W.M. provided input on the physics of quantum transport. F.B.,
C.B., L.B., M.H., M.S., and T.K. provided input on the physics of optical
spectroscopy. W.B., J.B., and E.M.H. developed the theoretical model at the
early stages of the project. W.B., T.K., and L.W.M. managed the project. The
writing of the paper was led by W.B. with input from all authors.

##### Funding information

We acknowledge financial support from the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) in the project SFB 1170 _ToCoTronics_
(Project ID 258499086) and in the Würzburg-Dresden Cluster of Excellence on
Complexity and Topology in Quantum Matter _ct.qmat_ (EXC 2147, Project ID
39085490).

##  Appendix A Implementation details

###  A.1 Parallelization

####  A.1.1 Basic method

The Python programming language was not designed with parallelization in mind.
The global interpreter lock (GIL) is often considered (see, e.g., Ref. [44])
as a major obstacle for running parallelized Python code. However, the
multiprocessing module avoids this obstacle by running Python code in separate
parallel subprocesses, such that each process is affected only by its own GIL.

In kdotpy, the multiprocessing module is used to parallelize iterative tasks
of the form

    data = [f(x, *f_args, **f_kwds) for x in vals]

where each call to the function f is an independent task. The main loop in
kdotpy is the iterated diagonalization of the Hamiltonian over the momentum
$k$ and magnetic field $B$ values of the grid. Due to the uniform nature of
these tasks, this iteration lends itself well for parallelization.

In kdotpy, the function `parallel_apply()` in the parallel submodule
facilitates parallelization over the diagonalization function by implementing
the necessary boilerplate code around the multiprocessing.Pool object. The
following code illustrates the recipe implemented in `parallel_apply()`
(simplified compared to the actual code for the sake of illustration):

    pool = multiprocessing.Pool(processes = num_processes)
    output = [
        pool.apply_async(f, args=(x,) + f_args, kwds = f_kwds) for x in vals
    ]
    while True:
       jobsdone = sum(1 for x in output if x.ready())
       if jobsdone >= n:
           break
       print(f"{jobsdone} / {n}")
       time.sleep(poll_interval)
    data = [r.get() for r in output]
    pool.close()
    pool.join()

Here, the while loop, that runs in the parent process, is used for process
monitoring by printing the number of completed tasks to standard output 777In
the actual code, an instance of the Progress class handles the progress
counter. In the code example here, we have replaced this by a print statement
for clarity.. This loop runs one iteration every `poll_interval` seconds in
the main process until all tasks are completed. The result data is collected
before the pool is joined and closed.

Unfortunately, multiprocessing does not handle keyboard interrupts (pressing
Ctrl-C to abort the program) and other signals properly out of the box. Thus,
`parallel_apply()` contains additional code in order to handle these
asynchronous events gracefully: Once a worker process is aborted or
terminated, the parent process of kdotpy detects it and aborts the full
calculation, including the remaining worker processes. If this is not done
correctly, worker processes can become _orphaned_ (i.e., detached from the
parent process) and continue working even when the parent process no longer
exists. These problems are mitigated by using a custom signal handler with
`parallel_apply()`.

####  A.1.2 Advanced method

The Task-Model framework is also built around the process pools of the
multiprocessing module, but with a higher level of abstraction, that allows a
more fine-grained control of how tasks are run. The kdotpy submodules tasks
and models define the following classes that constitute this framework:

  * •

Model _X_ class: This class stores the ‘recipe’, analogous to a
diagonalization function. For the different ‘models’ (1D, 2D, LL, etc.), a
separate class is derived from the ModelBase class. It stores the steps of the
recipe as a list of functions

        self.steps = [
        self.load_ddp,
        self.construct_ham,
        self.solve_ham,
        self.post_solve
    ]

where `self.load_ddp` enables loading of a DiagDataPoint object, and
`self.construct_ham`, `self.solve_ham`, and `self.post_solve` apply the
construction of the Hamiltonian, application of the eigensolver, and the
processing of eigenstates (i.e., the steps listed in Sec. 3.4.2).

  * •

Task class: This class contains the functions `self.worker_func`,
`self.callback`, and `self.error_callback` as well as references to the
process or thread pool. The `Task.run()` method calls `pool.apply_async()`,
where pool is either a process or a thread pool.

  * •

TaskManager class: The single instance of this class is responsible for
initializing, managing, and joining the process or thread pool, somewhat
similar to `parallel_apply()`. The primary job of the TaskManager is to send
the tasks in the queue to the available CPU or GPU workers and to join if all
tasks are completed. The handling of asynchronous events is also done by the
TaskManager. The TaskManager class is derived from PriorityQueue from the
queue module of Python.

The Task-Model framework has several benefits compared to `parallel_apply()`.
Firstly, the Model _X_ classes can be derived from. The derived class does not
need to redefine all steps, but only those one where it differs from its
parent class. Secondly, the data points are split into subtasks, so that each
subtasks can be executed where this is done most optimally (CPU or GPU). The
added flexibility of the Task-Model framework comes with a small price: As a
result of more overhead, the Task-Model framework is marginally slower than
`parallel_apply()`, but the difference is often barely noticeable.

###  A.2 Density of states

####  A.2.1 Triangular IDOS element

In Eq. (101), we have determined the fraction of the interval
$[k_{i},k_{i+1}]$ where a linear function is below a certain value $E$. For
two dimensions, the analogous problem would be considering a linearly
interpolated function on a triangle in the two dimensional plane. Let us label
the vertices $1,2,3$, and assume the function values are $e_{1,2,3}$ with
$e_{1}\leq e_{2}\leq e_{3}$ (without loss of generality). The momentum
coordinates of the vertices are irrelevant. The interpolated function is given
by

$e(u,v)=e_{1}+(e_{2}-e_{1})u+(e_{3}-e_{1})v=(1-u-v)e_{1}+u\,e_{2}+v\,e_{3}$
(A.1)

where $(u,v)$ are scaled coordinates with $0\leq u\leq 1$, $0\leq v\leq 1$,
and $u+v\leq 1$, such that the vertices of the triangle correspond to
$(u,v)=(0,0)$, $(1,0)$, and $(0,1)$. We are interested in the fraction $f(E)$
of the triangle where $e(u,v)<E$. This can be found by finding the
intersection points of the line $e(u,v)=E$ with the sides of the triangle and
calculating the area of the resulting polygon. The fraction $f(E)$ is this
area divided by $\frac{1}{2}$ (the area of the full triangle), and is given by

$f(E)=\left\\{\begin{array}[]{ll}1&\text{if }E\geq e_{3}\\\
1-\frac{(E-e_{3})^{2}}{(e_{3}-e_{2})(e_{3}-e_{1})}&\text{if }e_{2}\leq
E<e_{3}\\\ \frac{(E-e_{1})^{2}}{(e_{2}-e_{1})(e_{3}-e_{1})}&\text{if
}e_{1}\leq E<e_{2}\\\ 0&\text{if }E<e_{1}\end{array}\right.$ (A.2)

The conditions are such that the denominators are nonzero; for example, if
$e_{2}=e_{1}$, then the condition $e_{1}\leq E<e_{2}$ is never fulfilled.

The implementation in `triangle_idos_element()` takes as input the three-
dimensional array $e_{i,j,v}$, where $i,j$ label the elementary squares in
momentum space, and $v=1,2,3$ labels the three vertices. It applies the
following recipe.

  * •

The momentum space indices $(i,j)$ are flattened to $I$. The result is a two-
dimensional array $e_{I,v}$ of shape $(n_{k},3)$, where
$n_{k}=n_{k_{x}}n_{k_{y}}$ is the number of elementary squares in momentum
space.

  * •

Sort along the last axis so that $e_{I,1}\leq e_{I,2}\leq e_{I,3}$ is
satisfied for all $I$.

  * •

Evaluate $f^{(1)}=(E-e_{I,1})^{2}/(e_{I,2}-e_{I,1})(e_{I,3}-e_{I,1})$ and
$f^{(2)}=1-(E-e_{I,3})^{2}/(e_{I,3}-e_{I,2})(e_{I,3}-e_{I,1})$ (see Eq. (A.2))
for all $I$ and for all $E$ in the energy range. At the points where the
denominator vanishes, substitute $0$. In both cases, the resulting arrays f1
and f2 are two-dimensional arrays of shape $(n_{k},n_{E})$, where $n_{E}$ is
the number of energy values in the energy range.

  * •

Evaluate the conditions $E<e_{1}$, $E<e_{2}$, and $E<e_{3}$ as the two-
dimensional boolean arrays cond0, cond1, and cond2. Then evaluate $f_{I}(E)$
as

        f = np.where(cond0, zeros,
        np.where(cond1, f1, np.where(cond2, f2, ones))
    )

where `zeros = np.zeros_like(f1)` and `ones = np.ones_like(f1)`. The output is
an array of shape $(n_{k},n_{E})$, like all input arrays.

  * •

Subtract $1$ if the band is hole-like.

  * •

For all momenta indices $I$, where any $e_{I,v}$ ($v=1,2,3$) is undefined
(NaN), substitute $f_{I}(E)=0$ for all $E$.

####  A.2.2 Tetrahedral IDOS element

For three dimensions, we follow similar ideas to arrive at an expression for
$f(E)$. Let the energy values at the vertices of an elementary tetrahedron be
$e_{1,2,3,4}$, where $e_{1}\leq e_{2}\leq e_{3}\leq e_{4}$. The interpolated
function is given by
$e(u,v,w)=e_{1}+(e_{2}-e_{1})u+(e_{3}-e_{1})v+(e_{4}-e_{1})w$, with $0\leq
u,v,w\leq 1$ and $u+v+w\leq 1$. The fraction of the tetrahedra (with volume
$\frac{1}{6}$ that satisfies $e(u,v,w)<E$ is

$f(E)=\left\\{\begin{array}[]{ll}1&\text{if }E\geq e_{4}\\\
1-\frac{(E-e_{4})^{3}}{(e_{4}-e_{3})(e_{4}-e_{2})(e_{4}-e_{1})}&\text{if
}e_{3}\leq E<e_{4}\\\
\frac{(E-e_{1})^{3}}{(e_{2}-e_{1})(e_{3}-e_{1})(e_{4}-e_{1})}-\frac{(E-e_{2})^{3}}{(e_{2}-e_{1})(e_{3}-e_{2})(e_{4}-e_{2})}&\text{if
}e_{2}\leq E<e_{3}\text{ and }e_{1}<e_{2}\\\
\frac{(E-e_{1})^{2}}{(e_{3}-e_{1})(e_{4}-e_{1})}\left(1+\frac{e_{3}-E}{e_{3}-e_{1}}+\frac{e_{4}-E}{e_{4}-e_{1}}\right)&\text{if
}e_{2}\leq E<e_{3}\text{ and }e_{1}=e_{2}\\\
\frac{(E-e_{1})^{3}}{(e_{2}-e_{1})(e_{3}-e_{1})(e_{4}-e_{1})}&\text{if
}e_{1}\leq E<e_{2}\\\ 0&\text{if }E<e_{1}\end{array}\right.$ (A.3)

In `tetrahedral_idos_element()`, the arrays defining $f(E)$ are calculated
analogous to those in `triangle_idos_element()` but with some additional
intermediate steps to avoid recalculation of the products of $(e_{w}-e_{v})$
in the denominators. The conditions $E<e_{1}$, $E<e_{2}$, $E<e_{3}$, and
$E<e_{4}$ are applied as

    f = np.ones_like(f1)
    f[cond3] = f3[cond3]
    f[cond2] = f2[cond2]
    f[cond1] = f1[cond1]
    f[cond0] = 0.0

which is equivalent to a sequence of calls to np.where(), but with a better
performance for the larger arrays in three momentum dimensions.

Note: In the expression for $f(E)$ with $e_{2}\leq E<e_{3}$ and $e_{1}<e_{2}$,
it is possible to factor out $e_{2}-e_{1}$. Doing this yields

$f^{(2)}(E)=\frac{\left(e_{31}e_{41}+e_{41}(e_{3}-E)+e_{31}(e_{4}-E)\right)(E-e_{1})^{2}+e_{21}(E-e_{1})^{3}-3e_{21}e_{31}e_{41}(E-e_{1})+e_{21}^{2}e_{31}e_{41}}{e_{32}e_{42}e_{31}e_{41}},$
(A.4)

using the short-hand notation $e_{vw}=e_{v}-e_{w}$. Unlike Eq. (A.3), there is
no factor $e_{2}-e_{1}$ in the denominator, hence it is valid for
$e_{1}=e_{2}$ as well. If we substitute $e_{1}=e_{2}$, we retrieve the
expression for $e_{2}\leq E<e_{3}$ and $e_{1}=e_{2}$ of Eq. (A.3).

###  A.3 Charge neutrality point in symbolic Landau level mode

In dispersion mode and Landau level mode full, the charge neutrality point
always lies between the bands with band indices $-1$ and $1$, by definition.
In the Landau level mode sym, each state is characterized by a pair $(n,b)$ of
Landau level index $n=-2,-1,0,\ldots$ and band index $b$, so that finding the
charge neutrality point requires an extra step: For each state, we determine a
_universal band index_ $u$, equivalent to ordinary band indices in full Landau
level mode, with the property that the charge neutrality point lies between
the states with $u=-1$ and $u=+=1$.

The algorithm for finding the universal band indices is implemented in
`DiagDataPoint.get_ubindex()` and proceeds as follows. The states in
`DiagDataPoint` are sorted by eigenvalue. We define an array with ones for all
states with band index $b>0$ and an array with ones where $b<0$,

    pos = np.where(
        bindex_sort > 0,
        np.ones_like(bindex_sort),
        np.zeros_like(bindex_sort)
    )
    neg = 1 - pos

The arrays are summed cumulatively from below and above, respectively,

    npos = np.cumsum(pos)
    nneg = neg.sum() - np.cumsum(neg)

The array of universal band indices is then simply the difference between the
two cumulative sums

    ubindex = npos - nneg
    ubindex[ubindex <= 0] -= 1

where the nonpositive values have to be decreased by one in order to make sure
that the integers $u$ are nonzero. By definition, the resulting array is an
increasing sequence as function of energy, and identical to that obtained from
regular band indices in the Landau level mode full.

###  A.4 Solving Poisson’s equation

As already mentioned in Sec. 3.2.7, a way to incorporate electrostatic
potentials into the Hamiltonian is by parsing a set of boundary conditions.
This set is used to solve Poisson’s equation within the function
`solve_potential()` in potential.py, for ‘static’ as well as self-consistent
potentials (see Secs. 3.2.7 and 3.11, respectively). For the remainder of this
section we will refer to this potential as Hartree potential $V_{\mathrm{H}}$.
It should be emphasized that $V_{\mathrm{H}}$ parametrizes the potential
energy of the carriers in the electrostatic potential, not the electric
potential itself 888These two quantities differ by a factor of $-e$, which is
$-1$ in the chosen system of units. The Hartree potential $V_{\mathrm{H}}$
satisfies Poisson’s equation

$\partial_{z}\left[\varepsilon(z)\,\partial_{z}V_{\mathrm{H}}(z)\right]=\frac{e}{\varepsilon_{0}}\rho(z)$
(A.5)

[identical to Eq. (132)] with dielectric function/constant $\varepsilon(z)$,
elementary charge $e$, vacuum permittivity $\varepsilon_{0}$ and spatial
charge density $\rho(z)$. In the self-consistent case, as described in Secs.
3.11 and 3.7.8, $\rho(z)$ is extracted from all eigenvectors, whereas in the
‘static’ case $\rho(z)=0$.

The generalized solution can be calculated by integrating twice over $z$, with
generally different lower integration limits $z_{a}$ for the first integration
and $z_{b}$ for the second. The first integral over Eq. (A.5) yields

$\varepsilon(z)\,\partial_{z}V_{\mathrm{H}}(z)-\varepsilon(z_{a})\,\partial_{z}V_{\mathrm{H}}(z_{a})=\frac{e}{\varepsilon_{0}}\int_{z_{a}}^{z}dz^{\prime}\rho(z^{\prime})=\frac{e}{\varepsilon_{0}}\left[\mathcal{I}_{\rho}(z)-\mathcal{I}_{\rho}(z_{a})\right],$
(A.6)

where we use the shorthand notation
$\mathcal{I}_{\rho}(z)\equiv\int_{0}^{z}dz^{\prime}\rho(z^{\prime})$. As this
equation is valid for all $z$, it follows that

$\mathcal{C}(z)\equiv\varepsilon(z)\,\partial_{z}V_{\mathrm{H}}(z)-\frac{e}{\varepsilon_{0}}\mathcal{I}_{\rho}(z)$
(A.7)

is constant. We evaluate this constant at $z=z_{a}$, divide by
$\varepsilon(z)$ and integrate a second time to find the Hartree potential,

$\displaystyle V_{\mathrm{H}}(z)=\;$ $\displaystyle
V_{\mathrm{H}}(z_{b})+\frac{e}{\varepsilon_{0}}\left(\int_{0}^{z}dz^{\prime}\frac{\mathcal{I}_{\rho}(z^{\prime})}{\varepsilon(z^{\prime})}-\int_{0}^{z_{b}}dz^{\prime}\frac{\mathcal{I}_{\rho}(z^{\prime})}{\varepsilon(z^{\prime})}\right)+\mathcal{C}(z_{a})\left(\int_{0}^{z}dz^{\prime}\frac{1}{\varepsilon(z^{\prime})}-\int_{0}^{z_{b}}dz^{\prime}\frac{1}{\varepsilon(z^{\prime})}\right)$
(A.8)

We have written the integrals $\int_{z_{b}}^{z}$ explicitly as
$\int_{0}^{z}-\int_{0}^{z_{b}}$, which reflects the implementation in the
code.

In order to find a unique solution, two boundary conditions need to be given,
which fix $V_{\mathrm{H}}$ or $\partial_{z}V_{\mathrm{H}}$ (which can be
interpreted as electric field) at specific $z$ coordinates. The function
`solve_potential()` takes the boundary conditions as a set of keyword
arguments, namely a subset of values v1, v2, v3, dv1, dv2, and v12 and
coordinates z1, z2, and z3. Out of the former set, only certain combinations
of values may be set to a numerical value (the other being None); any other
combination will cause a ValueError exception to be raised. The following
combinations are valid:

  1. 1.

dv1 and v1, with z1 as coordinate: $\partial_{z}V_{\mathrm{H}}(z_{1})=\partial
V_{1}$ and $V_{\mathrm{H}}(z_{1})=V_{1}$. The $z$ coordinate for both boundary
conditions is the same, thus $z_{a}=z_{b}=z_{1}$.

$\displaystyle V_{\mathrm{H}}(z)=\;$ $\displaystyle
V_{1}+\frac{e}{\varepsilon_{0}}\left(\int_{0}^{z}dz^{\prime}\frac{\mathcal{I}_{\rho}(z^{\prime})}{\varepsilon(z^{\prime})}-\int_{0}^{z_{1}}dz^{\prime}\frac{\mathcal{I}_{\rho}(z^{\prime})}{\varepsilon(z^{\prime})}\right)$
$\displaystyle+\left(\varepsilon(z_{1})\,\partial
V_{1}-\frac{e}{\varepsilon_{0}}\mathcal{I}_{\rho}(z_{1})\right)\left(\int_{0}^{z}dz^{\prime}\frac{1}{\varepsilon(z^{\prime})}-\int_{0}^{z_{1}}dz^{\prime}\frac{1}{\varepsilon(z^{\prime})}\right).$
(A.9)

  2. 2.

dv2 and v2, with z2 as coordinate: Analogous to case 1, using
$z_{a}=z_{b}=z_{2}$ instead.

  3. 3.

dv1 and v2, with z1 and z2 as coordinates:
$\partial_{z}V_{\mathrm{H}}(z_{1})=\partial V_{1}$ and
$V_{\mathrm{H}}(z_{2})=V_{2}$. The $z$ coordinates for both boundary
conditions are different, thus $z_{a}=z_{1}$ and $z_{b}=z_{2}$.

$\displaystyle V_{\mathrm{H}}(z)=\;$ $\displaystyle
V_{2}+\frac{e}{\varepsilon_{0}}\left(\int_{0}^{z}dz^{\prime}\frac{\mathcal{I}_{\rho}(z^{\prime})}{\varepsilon(z^{\prime})}-\int_{0}^{z_{2}}dz^{\prime}\frac{\mathcal{I}_{\rho}(z^{\prime})}{\varepsilon(z^{\prime})}\right)$
$\displaystyle+\left(\varepsilon(z_{1})\,\partial
V_{1}-\frac{e}{\varepsilon_{0}}\mathcal{I}_{\rho}(z_{1})\right)\left(\int_{0}^{z}dz^{\prime}\frac{1}{\varepsilon(z^{\prime})}-\int_{0}^{z_{2}}dz^{\prime}\frac{1}{\varepsilon(z^{\prime})}\right).$
(A.10)

  4. 4.

v1 and dv2, with z1 and z2 as coordinates: Analogous to case 3, using
$z_{a}=z_{2}$ and $z_{b}=z_{1}$ instead.

  5. 5.

v1 and v2, with z1 and z2 as coordinates: $V_{\mathrm{H}}(z_{1})=V_{1}$ and
$V_{\mathrm{H}}(z_{2})=V_{2}$. Only potential values are given. The solution
is given by Eq. (A.8), but we cannot evaluate the constant
$\mathcal{C}(z_{a})$ directly, because the derivative
$\partial_{z}V_{\mathrm{H}}(z_{a})$ is unknown. By setting $z_{b}=z_{1}$ and
$z=z_{2}$ (or vice versa), we find that

$\displaystyle\mathcal{C}(z_{a})=\;$
$\displaystyle\frac{V_{\mathrm{H}}(z_{2})-V_{\mathrm{H}}(z_{1})-\frac{e}{\varepsilon_{0}}\left(\int_{0}^{z_{1}}dz^{\prime}\frac{\mathcal{I}_{\rho}(z^{\prime})}{\varepsilon(z^{\prime})}-\int_{0}^{z_{1}}dz^{\prime}\frac{\mathcal{I}_{\rho}(z^{\prime})}{\varepsilon(z^{\prime})}\right)}{\int_{0}^{z_{2}}dz^{\prime}\frac{1}{\varepsilon(z^{\prime})}-\int_{0}^{z_{1}}dz^{\prime}\frac{1}{\varepsilon(z^{\prime})}}.$
(A.11)

By substitution of $\mathcal{C}(z_{a})$ into Eq. (A.8), we thus find

$\displaystyle V_{\mathrm{H}}(z)=\;$ $\displaystyle
V_{1}+\frac{e}{\varepsilon_{0}}\left(\int_{0}^{z}dz^{\prime}\frac{\mathcal{I}_{\rho}(z^{\prime})}{\varepsilon(z^{\prime})}-\int_{0}^{z_{1}}dz^{\prime}\frac{\mathcal{I}_{\rho}(z^{\prime})}{\varepsilon(z^{\prime})}\right)$
$\displaystyle+\frac{V_{2}-V_{1}-\frac{e}{\varepsilon_{0}}\int_{z_{1}}^{z_{2}}dz^{\prime}\frac{\mathcal{I}_{\rho}(z^{\prime})}{\varepsilon(z^{\prime})}}{\int_{z_{1}}^{z_{2}}dz^{\prime}\frac{1}{\varepsilon(z^{\prime})}}\left(\int_{0}^{z}dz^{\prime}\frac{1}{\varepsilon(z^{\prime})}-\int_{0}^{z_{1}}dz^{\prime}\frac{1}{\varepsilon(z^{\prime})}\right).$
(A.12)

  6. 6.

v12 and v3, with z1, z2, and z3 as coordinates:
$V_{\mathrm{H}}(z_{2})-V_{\mathrm{H}}(z_{1})=V_{12}$ and
$V_{\mathrm{H}}(z_{3})=V_{3}$. We substitute $\mathcal{C}(z_{a})$ from Eq.
(A.11) into Eq. (A.8) with $z_{b}=z_{3}$, and obtain

$\displaystyle V_{\mathrm{H}}(z)=\;$ $\displaystyle
V_{3}+\frac{e}{\varepsilon_{0}}\left(\int_{0}^{z}dz^{\prime}\frac{\mathcal{I}_{\rho}(z^{\prime})}{\varepsilon(z^{\prime})}-\int_{0}^{z_{3}}dz^{\prime}\frac{\mathcal{I}_{\rho}(z^{\prime})}{\varepsilon(z^{\prime})}\right)$
$\displaystyle+\frac{V_{12}-\frac{e}{\varepsilon_{0}}\int_{z_{1}}^{z_{2}}dz^{\prime}\frac{\mathcal{I}_{\rho}(z^{\prime})}{\varepsilon(z^{\prime})}}{\int_{z_{1}}^{z_{2}}dz^{\prime}\frac{1}{\varepsilon(z^{\prime})}}\left(\int_{0}^{z}dz^{\prime}\frac{1}{\varepsilon(z^{\prime})}-\int_{0}^{z_{3}}dz^{\prime}\frac{1}{\varepsilon(z^{\prime})}\right).$
(A.13)

In `solve_potential()` the NumPy arrays densz [$\rho(z)$] and epsilonz
[$\varepsilon(z)$] first are integrated incrementally, yielding the arrays

    int_densz = integrate_arr(densz) * dz
    int_invepsilonz = integrate_arr(1. / epsilonz) * dz

which represent the integrals $\mathcal{I}_{\rho}(z)$ and
$\int_{0}^{z}dz^{\prime}1/\varepsilon(z^{\prime})$ for each individual value
of $z$. Subsequently,
$\int_{0}^{z}dz^{\prime}\mathcal{I}_{\rho}(z^{\prime})/\varepsilon(z^{\prime})$
is calculated as

    int_dens_over_epsz = integrate_arr(int_densz / epsilonz) * dz

The value of the integrals at the chosen $z$ coordinates ($z_{1}$, $z_{2}$,
$z_{3}$) are evaluated by using np.interp(), e.g.,

    int_dens_over_eps_z1, int_dens_over_eps_z2, int_dens_over_eps_z3 = \
        np.interp([z1, z2, z3], zval, int_dens_over_epsz)

where zval is the array of $z$ coordinates where the other arrays are defined.
A similar evaluation is done for `epsilonz`, `int_densz`, `int_invepsilon`,
followed by

    int_dens_over_epsz_1_2 = int_dens_over_eps_z2 - int_dens_over_eps_z1
    int_invepsilonz_1_2 = int_invepsilon_z2 - int_invepsilon_z1

Depending on the combination of boundary conditions, the corresponding
solution is selected from the six cases given above. Taking case 1 (dv1 and
v1) as example, the Hartree potential is then evaluated as

    int_const = epsilon_z1 * dv1 - eovereps0 * int_dens_z1
    vz = v1 + eovereps0 * (int_dens_over_epsz - int_dens_over_eps_z1)
            + int_const * (int_invepsilonz - int_invepsilon_z1)

where `int_const` is the integration constant $\mathcal{C}(z_{a})$ from Eq.
(A.7). The implementations for the other cases are analogous.

The boundary conditions are chosen automatically by parsing command line
arguments like vinner, vouter, vsurf or efield. In self-consistent
calculations the automatically determined boundary conditions can be manually
overwritten by using the command line argument potentialbc.

On a side note, one could also imagine solving the Poisson equation as a
matrix equation involving the vector $V_{\mathrm{H}}(z_{i})$ on the grid of
coordinates $z_{i}$, obtained by replacing the derivative operator by the
appropriate matrix (cf. Sec. 2.2.1). This results in a matrix-vector equation,
which can be solved by matrix inversion. In the continuum limit (grid
resolution $\Delta z\to 0$), this method converges to the same solution as the
numerical integration method detailed above. Our benchmarks show that the
numerical integration method outperforms the matrix method within the context
of kdotpy, with smaller numerical errors at finite resolutions $\Delta z$. We
thus decided to include the numerical integration approach only.

###  A.5 Extrema

####  A.5.1 Nine-point extremum solver

The function `nine_point_extremum_solver()` takes a $3\times 3$ grid of
momentum values $(k_{x,i^{\prime}},k_{y,j^{\prime}})$ values and a $3\times 3$
grid of energy values $e_{i^{\prime},j^{\prime}}$ as inputs. The indices
$i^{\prime}$ and $j^{\prime}$ are taken as $i^{\prime}=i-1,i,i+1$ and
$j^{\prime}=j-1,j,j+1$, where $i,j$ is a location where an extremum is
detected: We say there is a minimum (maximum) at $i,j$ if the values

$\displaystyle
e_{i+1,j},e_{i-1,j},e_{i,j+1},e_{i,j-1},e_{i+1,j+1},e_{i-1,j+1},e_{i+1,j-1},e_{i-1,j-1}$
(A.14)

are all $>e_{i,j}$ ($<e_{i,j}$). If this is the case, the
`nine_point_extremum_solver()` locates the extremum more precisely by fitting
Eq. (125) to the input data. Assume that the momenta are aligned on an equally
spaced cartesian grid, and define $\Delta
x=k_{x,i+1}-k_{x,i}=k_{x,i}-k_{x,i-1}$ and $\Delta
y=k_{y,j+1}-k_{y,j}=k_{y,j}-k_{y,j-1}$. First calculate the coefficients of
the Hessian matrix,

$\displaystyle a$ $\displaystyle=\frac{e_{i+1,j}-2e_{i,j}+e_{i-1,j}}{2(\Delta
x)^{2}},$ $\displaystyle c$
$\displaystyle=\frac{e_{i+1,j+1}-e_{i-1,j+1}-e_{i+1,j-1}+e_{i-1,j-1}}{4\,\Delta
x\,\Delta y},$ (A.15) $\displaystyle b$
$\displaystyle=\frac{e_{i,j+1}-2e_{i,j}+e_{i,j-1}}{2(\Delta y)^{2}}.$

We note that $c$ is the only coefficient where the $e_{i\pm 1,j\pm 1}$ appear;
the other linear combinations of $e_{i\pm 1,j\pm 1}$ are not considered. We
calculate the auxiliary variables

$X=\frac{e_{i+1,j}-e_{i-1,j}}{2\Delta x},\qquad
Y=\frac{e_{i,j+1}-e_{i,j-1}}{2\Delta y}$ (A.16)

and use them to find the location $(k_{x,0},k_{y,0})$ of the extremum in
momentum space,

$k_{x,0}=k_{x,i}+\frac{cY-2bX}{4ab-c^{2}},\qquad
k_{y,0}=k_{y,i}+\frac{cX-2aY}{4ab-c^{2}}.$ (A.17)

The denominator $4ab-c^{2}$ is the determinant of the Hessian matrix. If the
location defined by Eq. (A.17) lies outside of the rectangle
$[k_{x,i-1},k_{x,i+1}]\times[k_{y,j-1},k_{y,j+1}]$, raise the warning _“Poorly
defined extremum found”_ and use the location $(k_{x,i}-X/2a,k_{y,j}-Y/2b)$
instead, i.e., effectively setting $c=0$. The energy value $f_{0}$ of the
extremum is

$f_{0}=e_{i,j}-a(k_{x,i}-k_{x,0})^{2}-b(k_{y,j}-k_{y,0})^{2}-c(k_{x,i}-k_{x,0})(k_{y,j}-k_{y,0}).$
(A.18)

The function `nine_point_extremum_solver()` returns the result as $f_{0}$,
$(k_{x,0},k_{y,0})$, $(a,b,c)$.

The band masses are calculated from the eigenvalues of the Hessian matrix,
constructed from $(a,b,c)$. For cartesian coordinates, the eigenvalues are
equal to $\lambda_{\pm}=\frac{1}{2}(a+b\pm\sqrt{(a-b)^{2}+c^{2}})$. For polar
coordinates, $k_{x}$ and $k_{y}$ should be interpreted as $k_{r}=\lvert
k\rvert$ and $k_{\phi}$ respectively. In order to obtain band masses for
cartesian coordinates, the elements of the Hessian matrix must be rescaled as

$h_{\mathrm{polar}}=\begin{pmatrix}2a&c/k_{r}\\\
c/k_{r}&2b/k_{r}^{2}\end{pmatrix},$ (A.19)

where the factor $1/k_{r}$ comes from the Jacobian of the coordinate
transformation between polar and cartesian coordinates (assuming the angular
coordinate is in radians). The band masses are then obtained from the
eigenvalues of $h_{\mathrm{polar}}$.

####  A.5.2 Nineteen-point extremum solver

The function `nineteen_point_extremum_solver()` for three dimensions
effectively applies the method of the nine-point extremum solver in the
$k_{x}k_{y}$, $k_{x}k_{z}$, and $k_{y}k_{z}$ planes separately. The elements
of the Hessian matrix $h$ [the $3\times 3$ matrix in Eq. (127)] are calculated
as

$\displaystyle a$
$\displaystyle=\frac{e_{i+1,j,l}-2e_{i,j,l}+e_{i-1,j,l}}{2(\Delta x)^{2}},$
$\displaystyle d$
$\displaystyle=\frac{e_{i+1,j+1,l}-e_{i-1,j+1,l}-e_{i+1,j-1,l}+e_{i-1,j-1,l}}{4\,\Delta
x\,\Delta y},$ $\displaystyle b$
$\displaystyle=\frac{e_{i,j+1,l}-2e_{i,j,l}+e_{i,j-1,l}}{2(\Delta y)^{2}},$
$\displaystyle e$
$\displaystyle=\frac{e_{i+1,j,l+1}-e_{i-1,j,l+1}-e_{i+1,j,l-1}+e_{i-1,j,l-1}}{4\,\Delta
x\,\Delta z},$ (A.20) $\displaystyle c$
$\displaystyle=\frac{e_{i,j,l+1}-2e_{i,j,l}+e_{i,j,l-1}}{2(\Delta z)^{2}},$
$\displaystyle f$
$\displaystyle=\frac{e_{i,j+1,l+1}-e_{i,j-1,l+1}-e_{i,j+1,l-1}+e_{i,j-1,l-1}}{4\,\Delta
y\,\Delta z}.$

The coefficients $d$, $e$, and $f$ involve the twelve pink-coloured points in
Fig. 9(c), but for each group of four, only one linear combination is
considered. We note that the corner points $e_{i\pm 1,j\pm 1,l\pm 1}$
[coloured black in Fig. 9(c)] are included in the input arguments, but not
considered.

If the Hessian matrix is non-singular (in `nineteen_point_extremum_solver()`,
we use the condition $\lvert\det h\rvert>10^{-6}$), we find the momentum
location as

$\begin{pmatrix}k_{x,0}\\\ k_{y,0}\\\
k_{z,0}\end{pmatrix}=\begin{pmatrix}k_{x,i}\\\ k_{y,j}\\\
k_{z,l}\end{pmatrix}+h^{-1}\begin{pmatrix}-X\\\ -Y\\\ -Z\end{pmatrix}$ (A.21)

where $h^{-1}$ is the inverse of the Hessian matrix and
$X=(e_{i+1,j,l}-e_{i-1,j,l})/2\Delta x$, $Y=(e_{i,j+1,l}-e_{i,j-1,l})/2\Delta
y$, and $Z=(e_{i,j,l+1}-e_{i,j,l-1})/2\Delta z$ [analogous to Eq. (A.16)]. If
$\lvert\det h\rvert\leq 10^{-6}$ or the point lies outside of the box
$[k_{x,i-1},k_{x,i+1}]\times[k_{y,j-1},k_{y,j+1}]\times[k_{z,l-1},k_{z,l+1}]$,
the _“Poorly defined extremum found”_ warning is raised and the location
replaced by $(k_{x,i}-X/2a,k_{y,j}-Y/2b,k_{z,l}-Z/2c)$. The energy value
$f_{0}$ of the extremum is given by

$f_{0}=e_{i,j,l}-\begin{pmatrix}\tilde{k}_{x}&\tilde{k}_{y}&\tilde{k}_{z}\end{pmatrix}\begin{pmatrix}2a&d&e\\\
d&2b&f\\\ e&f&2c\end{pmatrix}\begin{pmatrix}\tilde{k}_{x}\\\ \tilde{k}_{y}\\\
\tilde{k}_{z}\end{pmatrix},$ (A.22)

with
$(\tilde{k}_{x},\tilde{k}_{y},\tilde{k}_{z})=(k_{x,i},k_{y,j},k_{z,l})-(k_{x,0},k_{y,0},k_{z,0})$.
The `nineteen_point_extremum_solver()` function returns the result as $f_{0}$,
$(k_{x,0},k_{y,0},k_{z,0})$, $(a,b,c,d,e,f)$. The band masses are calculated
from the eigenvalues of the Hessian matrix. If the coordinate system is
cylindrical or spherical, the appropriate coordinate transformation is applied
prior to finding the eigenvalues.

##  Appendix B Reference

###  B.1 Units and constants

####  B.1.1 Units

In kdotpy, physical quantities are represented by numerical data types without
the explicit specification of units. For this reason, it is necessary to agree
on a system of units such that calculations that physical quantities can be
done intuitively, possibly without the need of conversion factors. Whereas a
system where units are treated explicitly would be more fail-safe, it would
create a lot of overhead which may reduce performance.

We choose a system of units such that physical quantities in context of solid-
state physics have reasonable values. The basic units are:

  * •

Length: $\mathrm{nm}$

  * •

Time: $\mathrm{ns}$

  * •

Energy: $\mathrm{meV}$ (_not_ $\mathrm{eV}$)

  * •

Voltage: $\mathrm{mV}$ (_not_ $\mathrm{V}$)

  * •

Temperature: $\mathrm{K}$

  * •

Electric charge: $\mathrm{e}$

  * •

Magnetic field: $\mathrm{T}$

With _magnetic field_ we mean _magnetic flux density_ , more accurately
speaking. The unit for magnetic flux density in this unit system is defined as
Tesla, $\mathrm{T}=\mathrm{V}\,\mathrm{s}/\mathrm{m}^{2}$. This unit is
‘incompatible’ with the combination of units for voltage, length, and time,
which would be $\mathrm{mV}\,\mathrm{ns}/\mathrm{nm}^{2}$. The conversion,
given by $1\,\mathrm{T}=10^{-6}\,\mathrm{mV}\,\mathrm{ns}/\mathrm{nm}^{2}$ is
done internally by kdotpy when it interprets magnetic field values.

From the basic units, we derive units for other physical quantities, for
example,

  * •

Momentum (more appropriately wave vector) $k$: $\mathrm{nm}^{-1}$

  * •

Density (such as particle/carrier density): $\mathrm{nm}^{-1}$,
$\mathrm{nm}^{-2}$, or $\mathrm{nm}^{-3}$ depending on dimensionality

  * •

Density of states: $\mathrm{nm}^{-1}\,\mathrm{meV}^{-1}$,
$\mathrm{nm}^{-2}\,\mathrm{meV}^{-1}$, or
$\mathrm{nm}^{-3}\,\mathrm{meV}^{-1}$

  * •

Charge density: $e\,\mathrm{nm}^{-1}$, $e\,\mathrm{nm}^{-2}$, or
$e\,\mathrm{nm}^{-3}$

  * •

Electric field: $\mathrm{mV}/\mathrm{nm}$

  * •

Velocity: $\mathrm{nm}/\mathrm{ns}=\mathrm{m}/\mathrm{s}$

Inputs and outputs use the basic and derived units if not explicitly stated
otherwise.

####  B.1.2 Physical constants

The following physical constants are defined in physconst.py.

  * •

`m_e = 0.510998910e9` in $\mathrm{meV}$: Electron mass in energy equivalents
($E=m_{\mathrm{e}}c^{2}$).

  * •

`e_el = 1.6021766208e-19` in $\mathrm{C}$: Elementary charge $e$. We define
this value to be positive.

  * •

`cLight = 299792458.` in $\mathrm{nm}/\mathrm{ns}=\mathrm{m}/\mathrm{s}$:
Speed of light $c$.

  * •

`hbar = 6.582119514e-4` in $\mathrm{meV}\,\mathrm{ns}$: Reduced Planck
constant $\hbar$.

  * •

`hbarm0 = hbar**2 * cLight**2 / m_e / 2` in $\mathrm{meV}\,\mathrm{nm^{2}}$:
$\hbar^{2}/2m_{e}=38.09982350(23)\,\mathrm{meV}\,\mathrm{nm^{2}}$. The
approximate value $38\,\mathrm{meV}\,\mathrm{nm}^{2}$ is useful for making
estimates.

  * •

`eoverhbar = 1e-6 / hbar` in $1/(\mathrm{T}\,\mathrm{nm}^{2})$: $e/\hbar$
defined such, that eoverhbar times magnetic field in T yields a density in
$\mathrm{nm}^{-2}$. This takes into account the conversion factor $10^{-6}$ in
$\mathrm{T}=10^{-6}\,\mathrm{mV}\,\mathrm{ns}/\mathrm{nm}^{2}$. The product
eoverhbar * A, with the vector potential $\mathbf{A}$ in
$\mathrm{T}\,\mathrm{nm}$, yields a value in $\mathrm{nm}^{-1}$, as
appropriate for a momentum quantity. This appears in the Peierls substitution,
for example. The product eoverhbar * B, with the magnetic field $B$ in
$\mathrm{T}$, yields a quantity in units of $\mathrm{nm}^{-2}$, cf.
$l_{B}^{2}=\hbar/eB$

  * •

`muB = 5.7883818012e-2` in $\mathrm{meV}/\mathrm{T}$: Bohr magneton
$\mu_{B}=e\hbar/2m_{e}$.

  * •

`kB = 8.6173303e-2` in $\mathrm{meV}/\mathrm{K}$: Boltzmann constant $k_{B}$.

  * •

`eovereps0 = 1.80951280207e4` in $\mathrm{mV}\,\mathrm{nm}$:
$e/\varepsilon_{0}$, electron charge divided by permittivity constant (also
called vacuum permittivity).

  * •

`gg = 2` (dimensionless): Gyromagnetic ratio.

  * •

`r_vonklitzing = 25812.8074555` in $\Omega$ (ohm): Von Klitzing constant
$R_{\mathrm{K}}$, resistance value corresponding to one quantum of
conductance.

These values have been extracted from the NIST Reference on Constants, Units,
and Uncertainty, revision 2014, see Ref. [45] and online at
https://physics.nist.gov/cuu/Constants/index.html. The revised NIST values
from 2018 may deviate slightly (but negligibly) from the values listed here
and in physconst.py.

###  B.2 Material parameters

####  B.2.1 Chemistry

  * •

compound: The chemical formula of the compound, for example HgTe.

  * •

elements: Comma-separated list of the elements, for example Hg, Te. If
elements is not given, it is determined automatically from compound (and vice
versa).

  * •

composition: Numbers for each of the elements in the chemical formula, that
indicate their (stoichiometric) proportions in the compound. The numbers must
be separated by commas and may be a function of the variables x, y, and/or z.
For example, for a crystal with molecular formula Hg1-xCdxTe, use compound =
HgCdTe and composition = 1-x, x, 1.

####  B.2.2 Special commands

  * •

copy: Copy all valid parameters of another material into the present one. The
value is the id of the source material. For example

        [mat2]
    copy = mat1

copies all parameters from material mat1 into material mat2. Each parameter
may be subsequently overwritten manually.

  * •

linearmix: Define a material as a linear combination of two others, where the
parameters are linearly interpolated. The value is a 3-tuple of the form mat1,
mat2, var, where mat1 and mat2 are the source materials and var is the
interpolation variable (x, y, or z). For example

        [HgCdTe]
    linearmix = HgTe, CdTe, x

makes material HgCdTe a linear combination of HgTe and CdTe. Each parameter
$p$ is interpolated as
$p_{\mathrm{HgCdTe}}=(1-x)\,p_{\mathrm{HgTe}}+x\,p_{\mathrm{CdTe}}$. Each
parameter may be subsequently overwritten manually.

####  B.2.3 Band energies

  * •

Ev: ‘Valence’ band energy in meV. This is the energy of the $\Gamma_{8}$
orbitals at $k=0$ for the _unstrained bulk_ material. For inverted materials,
Ev indicates the energy of the $\Gamma_{8}$ orbitals, not that of the actual
valence band.

  * •

Ec: ‘Conduction’ band energy in meV. This is the energy of the $\Gamma_{6}$
orbitals at $k=0$ for the _unstrained bulk_ material. For inverted materials,
Ec indicates the energy of the $\Gamma_{6}$ orbitals, not that of the actual
conductance band.

  * •

`delta_so`: ‘Split-off’ energy difference $\Delta_{\mathrm{SO}}$ in meV
between the $\Gamma_{8}$ and $\Gamma_{7}$ orbitals in the unstrained bulk
material. To put it more precisely, the energy of the $\Gamma_{8}$ orbital at
$k=0$ is `Ev - delta_so`.

####  B.2.4 Quadratic terms

  * •

P: Kane matrix element $P=(\hbar/m_{0})\langle S|p_{x}|X\rangle$ in meV nm.
Note: The definition of $P$ differs between references by a sign and/or a
factor of $\mathrm{i}$. Here, we take a positive real value for P.

  * •

gamma1: Luttinger parameter $\gamma_{1}$, dimensionless value. This Luttinger
parameter describes the spherically isotropic component of the (bare) band
mass of the $\Gamma_{8}$ orbitals.

  * •

gamma2: Luttinger parameter $\gamma_{2}$, dimensionless value. This Luttinger
parameter describes the component proportional to
$k_{x}^{2}+k_{y}^{2}-2k_{z}^{2}$ for the $\Gamma_{8}$ orbitals. This component
is axially, but not spherically symmetric.

  * •

gamma3: Luttinger parameter $\gamma_{3}$, dimensionless value. This Luttinger
parameter describes the component proportional to $k_{x}^{2}-k_{y}^{2}$ for
the $\Gamma_{8}$ orbitals. This component breaks axial symmetry.

  * •

F: Band mass parameter for the $\Gamma_{6}$ orbitals. The bare band mass term
is $\frac{\hbar^{2}}{2m_{0}}(2F+1)(k_{x}^{2}+k_{y}^{2}+k_{z}^{2})$. In $2F+1$,
the $1$ is the contribution of free electrons and $2F$ is the contribution
from perturbative corrections from remote bands in $\mathbf{k}\cdot\mathbf{p}$
theory.

  * •

kappa: See Section _Magnetic couplings_ below.

####  B.2.5  Bulk-inversion asymmetry

  * •

`bia_c`: Linear bulk-inversion asymmetry coefficient $C$ in meV nm-1.

  * •

`bia_b8p`: Quadratic bulk-inversion asymmetry coefficient $B_{8v}^{+}$ in meV
nm-2.

  * •

`bia_b8m`: Quadratic bulk-inversion asymmetry coefficient $B_{8v}^{-}$ in meV
nm-2.

  * •

`bia_b7`: Quadratic bulk-inversion asymmetry coefficient $B_{7v}$ in meV nm-2.

####  B.2.6 Magnetic couplings

  * •

ge: Gyromagnetic factor $g_{e}$ for the $\Gamma_{6}$ orbitals. The magnetic
coupling is also known as the Zeeman term and is equal to
$g_{\mathrm{e}}\mu_{\mathrm{B}}\mathbf{B}\cdot\mathbf{S}$. The value of
$g_{e}$ contains contributions from the free electron (equal to $2$) and
perturbative corrections from remote bands.

  * •

kappa: Gyromagnetic factor $\kappa$ for the $\Gamma_{8}$ and $\Gamma_{7}$
orbitals. The magnetic coupling is of the form
$-2\mu_{\mathrm{B}}\kappa\mathbf{J}\cdot\mathbf{B}$ in the ($\Gamma_{8}$,
$\Gamma_{8}$) block of the Hamiltonian and analogous couplings in the
($\Gamma_{8}$, $\Gamma_{7}$) and ($\Gamma_{7}$, $\Gamma_{7}$) blocks. The
coefficient kappa also appears as a non-magnetic term containing a commutator
$[\kappa,k_{z}]$, which contributes only at material interfaces.

  * •

q: Coefficient $q$ of the $\Gamma_{8}$ magnetic coupling
$-2\mu_{\mathrm{B}}q\mathbf{\mathcal{J}}\cdot\mathbf{B}$, where
$\mathbf{\mathcal{J}}=(J_{x}^{3},J_{y}^{3},J_{z}^{3})$.

  * •

`exch_yNalpha`: Exchange energy parameter in meV for the $\Gamma_{6}$
orbitals. The paramagnetic exchange interaction is relevant for paramagnetic
materials like (Hg,Mn)Te. The interaction strength typically depends linearly
on the Mn concentration ($y$ in Ref. [3], hence the y in the parameter name).
In the materials file, this dependence must be given explicitly, e.g., by
setting `exch_yNalpha = y * 400`. The given dependence may be chosen non-
linear if desired. See Ref. [3] for more details on the physics of this
paramagnetic coupling in (Hg,Mn)Te.

  * •

`exch_yNbeta`: Exchange energy parameter in meV for the $\Gamma_{8}$ and
$\Gamma_{7}$ orbitals.

  * •

`exch_g`: The gyromagnetic factor $g_{\mathrm{ex}}$ in the argument of the
Brillouin function. For (Hg,Mn)Te, this quantity is usually written as
$g_{\mathrm{Mn}}$ see Eq. (58) and Ref. [3]. The default value is 2.

  * •

`exch_TK0`: The offset temperature $T_{0}$ in K that appears in the
denominator of the argument of the Brillouin function. The default value is
$10^{-6}$. See Sec. 2.5.2.

####  B.2.7 Electrostatic properties

  * •

`diel_epsilon`: Dielectric constant $\varepsilon_{\mathrm{r}}$, also known as
relative permittivity. The value is dimensionless. The present version of
kdotpy implements a constant value only, i.e., no frequency dependence is
considered.

  * •

`piezo_e14`: Piezo-electric constant $e_{14}$ in $e$ nm-2. Note that this is
the only independent component of the piezo-electric tensor for the zincblende
crystal structure. To convert from units of C m-2, multiply the value by
$10^{18}e$, for example `piezo_e14 = 0.035 * 1e18 * e_el`.

####  B.2.8 Lattice constant and strain

  * •

a: Lattice constant in nm of the unstrained material.

  * •

`strain_C1`: Strain parameter (deformation potential) $C_{1}$. The
corresponding strain term is $C_{1}\mathop{\mathrm{tr}}\epsilon$ acting on the
$\Gamma_{6}$ orbitals. In Refs. [2, 3], the notation $C$ is used.

  * •

`strain_Dd`: Strain parameter (deformation potential) $D_{d}$. The
corresponding strain term is $D_{d}\mathop{\mathrm{tr}}\epsilon$ acting on the
$\Gamma_{8}$ and $\Gamma_{7}$ orbitals. In Refs. [2, 3], the notation $a$ is
used.

  * •

`strain_Du`: Strain parameter (deformation potential) $D_{u}$. The
corresponding strain term is a linear combination of the diagonal entries
$\epsilon_{ii}$ of the strain tensor, and acts on the $\Gamma_{8}$ and
$\Gamma_{7}$ orbitals. In Refs. [2, 3], the notation $b=-\frac{2}{3}D_{u}$ is
used.

  * •

`strain_Duprime`: Strain parameter (deformation potential) $D^{\prime}_{u}$.
The corresponding strain term is a linear combination of the off-diagonal
entries $\epsilon_{ij}$ (shear components) of the strain tensor, and acts on
the $\Gamma_{8}$ and $\Gamma_{7}$ orbitals. In Refs. [2, 3],
$d=-\frac{2}{\sqrt{3}}D^{\prime}_{u}$ is used.

###  B.3 Lattice orientation

The $\mathrm{SO}(3)$ rotation matrix $R$ can be parametrized in different
ways, for example directly in terms of the device coordinate axis in terms of
the lattice direction, or in terms of Euler angles. The command line option
orientation is the most generic way to input a rotation. By giving up to three
angles or directions, the crystal lattice can be rotated to any possible
orientation with respect to the device coordinate system:

  * •

orientation #ang: Rotation around $z$, equivalent to stripangle.

  * •

orientation #ang #ang: Tilt the $z$ axis, then rotate around the $z$ axis.

  * •

orientation #ang #ang #ang: Euler rotation $z,x,z$, i.e., rotate around the
$c$ crystal axis, tilt the $z$ axis, and rotate around the $z$ axis.

  * •

orientation #dir: Longitudinal direction $x$, equivalent to stripdir.

  * •

orientation - #dir: Growth direction $z$.

  * •

orientation #dir #dir: Longitudinal and growth direction $x$, $z$.

  * •

orientation - #dir #dir: Transversal and growth direction $y$, $z$.

  * •

orientation #dir #dir #dir: Longitudinal, transversal, and growth direction
$x$, $y$, $z$.

  * •

orientation #dir #ang or orientation #ang #dir: Growth direction $z$ and
rotation around $z$.

Angles #ang are entered an as explicit floating point number containing a
decimal sign . or with the degree symbol (∘ or d). Directions #dir are
triplets either of digits without separators and possibly with minus signs
(e.g., 100, 111, 11-2, -110) or of numbers separated by commas without spaces
(e.g., 1,1,0 or 10,-10,3). If multiple directions are given as arguments, they
must be pairwise orthogonal, which is tested by calculating the inner products
between them.

The command line arguments stripangle and stripdir are simplified options for
rotation of a strip geometry around the $z$ axis.

###  B.4 Custom eigensolvers

As discussed in Sec. 3.4.1, the diagonalization involves two components: The
actual eigensolver and an efficient LU solver for the shift-and-invert step.
kdotpycan be configured to use custom solvers provided by external packages
for both steps. In this section, we list them and discuss their strengths and
weaknesses.

####  B.4.1 Eigensolvers

  * •

eigsh from SciPy: The default implementation that ships with the SciPy module.
It interfaces to the external ARPACK package [38] that implements the
Lanczos/Arnoldi algorithm.
Pro: It is stable and there are no convergence problems. It comes with SciPy,
so it does not require any steps to set up.
Contra: A single calculation can only make use of one CPU core (multi-core
usage possible by parallelizing over $k$ or $B$ points). The memory bandwidth
may be a bottleneck for large problems, due to gemm operations (matrix-matrix
products with general matrices in BLAS [46]) during eigenvector
orthonormalization, which is the slowest operation per Lanzcos/Arnoldi
iteration.

  * •

eigsh from CuPy [47]: Also a Lanzcos/Arnoldi algorithm, written in pure
Python. The CuPy library mimics, extends and replaces the NumPy/SciPy
packages. Where possible and efficient, the computational workload is done on
a CUDA capable GPU. Performance is best on larger problem sizes. Depending on
the availability of single and double precision GPU cores, using single
precision floating-point arithmetic may be faster, at the expense of lower
precision. Using CuPy requires a CUDA capable GPU and a compatible version of
the CUDA library 999The CUDA API, provided by Nvidia, is an extension of C/C++
that facilitates computation on GPUs that support it..
Pro: It is much faster for large problems compared to the CPU eigsh solver,
for example about 6 times in non-axial LL mode, depending on setup and problem
parameters. The GPUs typically have higher memory bandwidth. TensorCores can
drastically speed up the gemm operation.
Contra: The RAM limits on GPU are generally tighter than those on CPU. The
solver can fail to converge, especially if single precision floating-point
arithmetic is used. This can be mitigated by switching to higher precision or
falling back to the CPU solver (done automatically by default), but this leads
to a reduction in speed. Setting up and configuring this solver requires more
effort.

  * •

FEAST solver: Either provided via Intel MKL (version 2.1) or by manual
compilation of a shared library (version 4.0).
Pro: Can be more efficient than default eigsh; no additional LU solver is
required.
Contra: This solver can be unreliable with densely clustered eigenstates,
which is frequently the case for kdotpy.

  * •

eigh from jax: Converts the sparse matrix to a dense matrix, then uses eigh
from the Google JAX library to solve on accelerators like GPUs using highly
optimized instructions. Scales across multiple accelerators. Currently
implemented as fixed double precision.
Pro: Extremely fast when a large fraction of eigenvalues are requested (in
fact, the algorithm does a full diagonalization internally).
Contra: Very memory intensive on GPUs. The entire Hamiltonian (including
solutions) needs to fit into GPU memory.

####  B.4.2 LU solvers

  * •

SuperLU: The default implementation shipped with current SciPy module.
Pro: It performs well and does not require additional effort to set up.
Contra: Can fail with a memory allocation error for very large matrix sizes.

  * •

UMFPACK: Uses external library shipped with scikit/UMFPACK. Previously, this
was used as default by SciPy, but this was changed due to license issues. If
UMFPACK is installed (scikit libraries and Python packages), SciPy chooses it
automatically for use with eigsh.
Pro: It is stable and efficient, and usually fast.
Contra: In some cases, it may rather be slower than SuperLU.

  * •

PARDISO: Provided by the Intel MKL library. It requires the Intel MKL
libraries to be installed along with a suitable Python package that links to
it.
Pro: It is fast, for example about 3 to 4 times with kdotpy 1d, efficient and
stable.
Contra: The official Python package is no longer maintained. It may require
more work to set up.

####  B.4.3 Performance considerations

Optimizing solver speed is connected to the problem to be solved, as well as
to the computation environment used for the calculation. For the very
different workloads created by various kdotpy problem sets and the large range
of hardware that kdotpy could be run on (from small dual core, low RAM PCs
over powerful GPU accelerated workstations up to multiple HPC cluster nodes),
it is not possible to identify a single eigensolver and a set of optimal
parameters that satisfies every scenario.

Automatic settings in kdotpy always try to use settings for highest stability,
maximum speed and minimal RAM requirements as a rough guideline (in order of
priority). However in most scenarios, performance data is not available and
has to be found by user experience. The following hardware parameters are
important aspects to that affect eigensolver performance:

  * •

CPU/GPU core clock speed: This parameter is mostly fixed. For some systems one
might be able to increase the speed by overclocking without affecting other
limiting factors, but this typically gives minor speed boosts on the order of
few percent.

  * •

CPU instruction sets: Vector and matrix BLAS operations profit much from
modern CPU instructions such as AVX. Intel MKL downthrottles AMD CPUs by not
making use of some of those features, but this behaviour can be disabled. An
older version of MKL might perform best in such cases, but should be
benchmarked against current versions.

On CPU cores: While Python itself has many pitfalls in terms of efficient
threading and multi core usage, prominently due to its global interpreter lock
(GIL), we can still make good use of additional CPU cores by solving multiple
matrices in parallel processes by parallelization over $k$ or $B$ points. As
long as we do not run into other limits listed below, this yields almost
linear performance improvements. Some of the external libraries (all ones
mentioned above, except ARPACK) can also use more cores through
multithreading, however this higher core usage is not always rewarded with
speed gains. A high degree of parallelization can be energy-inefficient due to
the extra overhead. On hardware with inadequate cooling, pushing the CPU to
its limits for an extensive time may also lead to performance degradation from
thermal throttling.

  * •

RAM size: The maximum amount of RAM in the system sets a hard limit on the
number of matrices that can be solved in parallel. Each independent problem
requires a similar amount of RAM space during calculation, that can not be
shared. Once the threshold is exceeded, the whole process is likely to crash.

  * •

CPU cache size and RAM bandwidth: This is the amount of data the CPU can
request from RAM (per second). For some calculations, it is possible, that
data for CPU operations can not be stored efficiently enough in the CPU cache
(cache miss) and a lot of data has to fetched from RAM. This is a limitation
for large matrix-matrix-multiplications (BLAS function gemm), as the
operation’s data cannot fit in cache. In non-axial LL mode, this can be
observed when requesting many eigenvalues (large $n_{\mathrm{eig}}$ for thick
layers, as the orthogonalization of eigenvectors (as part of the Lanczos
algorithm) is basically done using two gemm calls. Running more simultaneous
processes does not increase total solution speed, as it slows down each
process due to shared RAM access.

###  B.5 Observables

####  B.5.1 Spatial observables

The following observables are defined as integrals over spatial coordinates.
In the discrete coordinates used in kdotpy, they are implemented as
summations, e.g., $\sum_{i}\psi^{*}_{i}z_{i}\psi_{i}$ for the expectation
value $\langle z\rangle$.

  * •

y: Expectation value $\langle y\rangle$, where $y$ is the coordinate
transverse to the strip.

  * •

y2: Expectation value $\langle y^{2}\rangle$.

  * •

sigmay: Standard deviation of y, $\sigma_{y}=\sqrt{\langle
y^{2}\rangle-\langle y\rangle^{2}}$.

  * •

z: Expectation value $\langle z\rangle$, where $z$ is the coordinate parallel
to the growth direction.

  * •

z2: Expectation value $\langle z^{2}\rangle$.

  * •

sigmaz: Standard deviation of $z$, $\sigma_{y}=\sqrt{\langle
z^{2}\rangle-\langle z\rangle^{2}}$.

  * •

zif: Expectation value $\langle z_{\mathrm{if}}\rangle$, where
$z_{\mathrm{if}}$ is the (signed) distance to the nearest interface.

  * •

zif2: Expectation value $\langle z_{\mathrm{if}}^{2}\rangle$.

  * •

sigmazif: Standard deviation of $z_{\mathrm{if}}$,
$\sigma_{z_{\mathrm{if}}}=\sqrt{\langle z_{\mathrm{if}}^{2}\rangle-\langle
z_{\mathrm{if}}\rangle^{2}}$.

  * •

well: Probability density inside the well.

  * •

wellext: Probability density inside the well and the adjacent 2 nm of
barriers.

  * •

interface: Probability density less than 1 nm from each interface. This
quantity cannot be used to compare samples of different size.

  * •

interface10nm: Probability density less than 10 nm from each interface. This
quantity cannot be used to compare samples of different size.

  * •

custominterface[]: Probability density up to a length in nm, set by the
command-line argument custominterfacelength.

  * •

interfacechar: ‘Interface character’. Probability density less than 1 nm from
each interface divided by what this quantity would be for a uniform
probability density, i.e., for a normalized wave function with constant
magnitude given by $|\psi|^{2}\equiv 1/V$, where $V$ is the volume (size) of
the complete sample. If $\Omega$ denotes the domain near the interfaces (e.g.,
less than 1 nm away) and $V_{\mathrm{if}}=\mathrm{vol}(\Omega)$ its size, then
the interface character is $(V/V_{\mathrm{if}})\int_{\Omega}|\psi(r)|^{2}dr$.
Regardless of sample size, values larger than 1 indicate strong interface
character.

  * •

interfacechar10nm: Interface character, like interfacechar, but with an
interface region of 10 nm instead of 1 nm.

  * •

custominterfacechar[]: Interface character, like interfacechar, with the size
of the interface region set by the command-line argument
custominterfacelength.

####  B.5.2 Inverse participation ratio

Generically, an inverse participation ratio (IPR) for a probability
distribution is defined in terms of the fourth moment (kurtosis) scaled by the
square of the second moment (variance). For a wave function $\psi(r)$ in terms
of spatial coordinate(s) $r$, we define

$\mathrm{IPR}=\frac{\left(\int|\psi(r)|^{2}dr\right)^{2}}{\int|\psi(r)|^{4}dr}=\left(\int|\psi(r)|^{4}dr\right)^{-1},$
(B.1)

where the second equality assumes normalization of the wave function. The
resulting IPR has dimensions of length (for a one-dimensional spatial
coordinate). kdotpy provides IPR observables for three choices of coordinates:

  * •

ipry : IPR over the coordinate $y$. If necessary, we first integrate over the
coordinate z in order to obtain $|\psi(y)|^{2}$. The dimensionful variety has
units of length. For the dimensionless variety, divide by the width $w$ in $y$
direction.

  * •

iprz : IPR over the coordinate $z$. If necessary, we first integrate over the
coordinate y in order to obtain $|\psi(z)|^{2}$. The dimensionful variety has
units of length. For the dimensionless variety, divide by the thickness $d$ in
$z$ direction.

  * •

ipryz : IPR over the coordinates $(y,z)$. The dimensionful variety has units
of area (length squared). For the dimensionless variety, divide by $wd$.

####  B.5.3 Internal degrees of freedom

The following observables are defined in terms of internal degrees of freedom,
i.e., (combinations of) orbitals. This includes spin operators.

  * •

sx, sy, sz: Spin expectation value $\langle S_{i}\rangle$ ($i=x,y,z$). This is
proper spin in units of $\hbar$, so that the range of possible values is
$[-\frac{1}{2},\frac{1}{2}]$.

  * •

jx, jy, jz: Expectation values $\langle J_{i}\rangle$ ($i=x,y,z$) of the total
angular momentum. The range of possible values is
$[-\frac{3}{2},\frac{3}{2}]$.

  * •

yjz: Expectation values $\langle yJ_{z}\rangle$, i.e., the product of the
coordinate $y$ and the angular momentum $J_{z}$. This is roughly chirality for
edge states.

  * •

gamma6, gamma7, gamma8. Expectation values $\langle P_{\Gamma_{i}}\rangle$
($i=6,7,8$) of the projection to the $\Gamma_{i}$ orbitals. This is the
probability density in these orbitals.

  * •

gamma8h: Expectation values $\langle P_{\Gamma_{8H}}\rangle$ of the projection
to the orbital $|\Gamma_{8},j_{z}=\pm 3/2\rangle$. ‘H’ stands for heavy hole.

  * •

gamma8l: Expectation values $\langle P_{\Gamma_{8L}}\rangle$ of the projection
to the orbital $|\Gamma_{8},j_{z}=\pm 1/2\rangle$. ‘L’ stands for light hole.

  * •

orbital: Difference of gamma6 and gamma8. The value ranges between $-1$ for
purely $\Gamma_{8}$ states and $1$ for purely $\Gamma_{6}$ states.

  * •

jz6, jz7, jz8. Expectation values $\langle J_{z}P_{\Gamma_{i}}\rangle$
($i=6,7,8$) of the angular momentum in $z$ direction in the $\Gamma_{i}$
orbitals.

  * •

split: Expectation value of the ‘artificial split’ Hamiltonian, which can be
encoded as $\mathrm{sgn}(J_{z})$.

  * •

orbital[j]: The squared overlaps of the eigenstates within orbital number $j$,
where $j$ runs from 1 to $n_{\mathrm{orb}}$ (the number of orbitals, i.e., 6
or 8). Available for kdotpy 2d only, if the command line option orbitalobs is
given.

####  B.5.4 Plain parity and ‘isoparity’ operators

The following observables implement several plain parity operators (a purely
spatial operation) as well as ‘isoparity’ operators, which combine a spatial
and a spin component [48]. The latter are given as actions of the appropriate
representations of the relevant point group, namely the double group of
$T_{d}$ or one of its subgroups.

  * •

pz: Expectation value of the parity operator in the $z$ coordinate, i.e.,
$\langle\mathcal{P}_{z}\rangle=\int\psi^{*}(z)\psi(-z)dz$.

  * •

px, py: Expectation value of the parity operator in the coordinate $x$ or $y$,
i.e., $\langle\mathcal{P}_{x}\rangle$, $\langle\mathcal{P}_{y}\rangle$

  * •

pzy: Expectation value of the combined parity operator
$\langle\mathcal{P}_{z}\mathcal{P}_{y}\rangle$

  * •

isopz: Expectation value of isoparity $\langle\tilde{\mathcal{P}}_{z}\rangle$:
this is the parity operator combined with the diagonal matrix $Q$ in orbital
space, given by $Q=\mathrm{diag}(1,-1,1,-1,1,-1,-1,1)$ in the basis given by
Eq. 10. This is a conserved quantum number under a few assumptions [48].

  * •

isopx, isopy: Expectation value of in-plane isoparity
$\langle\tilde{\mathcal{P}}_{x}\rangle$,
$\langle\tilde{\mathcal{P}}_{y}\rangle$. This is the parity operator in x or
y, combined with the appropriate matrix $Q_{x}$, $Q_{y}$ in orbital space,
which can be found from the appropriate representations of the point group
[48].

  * •

isopzy: Expectation value of the combined isoparity operator
$\langle\tilde{\mathcal{P}}_{z}\tilde{\mathcal{P}}_{y}\rangle$

  * •

isopzw: Expectation value of a modified isoparity operator like isopz, that
acts only in the quantum well layer. Formally speaking, if the well ranges
from $z=z_{\mathrm{min}}$ to $z_{\mathrm{max}}$, then it maps $z^{\prime}=z-c$
($c=(z_{\mathrm{min}}+z_{\mathrm{max}})/2$ being the center of the well) to
$-z^{\prime}$. While isopz is not a conserved quantity for asymmetric
geometries (e.g., a well layer and two barriers with unequal thickness),
isopzw can remain almost conserved in that case. Due to incomplete confinement
in the well region, the eigenvalues may deviate significantly from $\pm 1$.

  * •

isopzs: Expectation value of a modified isoparity operator like isopz, that
acts only in a symmetric region around the centre of the quantum well layer.
It is like isopzw, but $z_{\mathrm{min}}$ and $z_{\mathrm{max}}$ are chosen
such that the range fits inside the stack, while keeping
$c=(z_{\mathrm{min}}+z_{\mathrm{max}})/2$ at the centre of the well. It tends
to have eigenvalues closer to $\pm 1$, because it generally covers the
probability density in (a part of) the barrier layers as well.

####  B.5.5 Hamiltonian terms

The following observables are expectation values of the individual terms in
the Hamiltonian.

  * •

hex: Exchange energy (expectation value of the exchange Hamiltonian).

  * •

hex1t: Exchange energy at $1\,\mathrm{T}$. Regardless of the actual value of
the $B$ field, give the expectation value of the exchange Hamiltonian
evaluated at $B_{z}=1\,\mathrm{T}$.

  * •

hexinf: Exchange energy in the large-field limit. Regardless of the actual
value of the $B$ field, give the expectation value of the exchange Hamiltonian
evaluated for $B_{z}\to\infty$. In this limit, the Brillouin function is
saturated at its maximum absolute value.

  * •

hz: Zeeman energy (expectation value of the exchange Hamiltonian).

  * •

hz1t: Zeeman energy at $1\,\mathrm{T}$. Regardless of the actual value of the
$B$ field, give the expectation value of the Zeeman Hamiltonian evaluated at
$B_{z}=1\,\mathrm{T}$.

  * •

hstrain: Expectation value of the strain Hamiltonian.

####  B.5.6 Landau levels

The following are observables based on the Landau level index $n$.

  * •

llindex: In the LL mode sym or legacy, this is the (conserved) LL index $n$.
The minimal value is $-2$. The lowest indices $n=-2,-1,0$ are incomplete,
i.e., for these indices, not all orbitals contribute.

  * •

llavg: In the LL mode full, the Hamiltonian is not diagonal in the basis of LL
indexed by $n$. This observable returns the expectation value $\langle
n\rangle$ in this basis.

  * •

llbymax: In the LL mode full, write the wave function in the basis of LL
indexed by $n$, and provide the value $n$ with the highest probability
density.

  * •

llmod4: In the LL mode full, analogous to llavg, returns the expectation value
$\langle n\mod 4\rangle$. While llavg is not a conserved quantity if axial
symmetry is broken, this observable is usually conserved if bulk inversion
symmetry is not broken.

  * •

llmod2: In the LL mode full, analogous to llavg, returns the expectation value