# Efficient computation of optical excitations in two-dimensional materials
with the Xatu code

Alejandro José Uría-Álvarez  Juan José Esteve-Paredes  M. A. García-Blázquez
Juan José Palacios  Departamento de Física de la Materia Condensada,
Universidad Autónoma de Madrid, 28049 Madrid, Spain  Instituto Nicolás
Cabrera, Condensed Matter Physics Centre (IFIMAC), 28049 Madrid, Spain

###### Abstract

Here we describe an efficient numerical implementation of the Bethe-Salpeter
equation to obtain the excitonic spectrum of semiconductors. This is done on
the electronic structure calculated either at the simplest tight-binding level
or through density funcional theory calculations based on local orbitals. We
use a simplified model for the electron-electron interactions which considers
atomic orbitals as point-like orbitals and a phenomenological screening. The
optical conductivity can then be optionally computed within the Kubo
formalism. Our results for paradigmatic two-dimensional materials such as hBN
and MoS2, when compared with those of more sophisticated first-principles
methods, are excellent and envision a practical use of our implementation
beyond the computational limitations of such methods.

###### keywords:

Exciton, Bethe-Salpeter Equation, Optics, Many-Body Physics, Localized
Orbitals, Tight-Binding

††journal: Computer Physics Communications

PROGRAM SUMMARY

Program Title: Xatu
CPC Library link to program files: (to be added by Technical Editor)
Developer’s repository link:
https://github.com/alejandrojuria/xatu
Code Ocean capsule: (to be added by Technical Editor)
Licensing provisions: GPLv3
Programming language: C++, Fortran, Python
Nature of problem: The exciton spectrum is obtained as the solution of the
Bethe-Salpeter equation for insulators and semi-conductors. Constructing the
equation involves determining the screening of the electrostatic interaction
and then determining the matrix elements of the interaction kernel, which are
computationally-intensive tasks, specially if one takes a purely ab-initio
approach.
Solution method: The Bethe-Salpeter equation can be efficiently set up and
solved assuming that the basis of the reference electronic structure
calculation, obtained either from tight-binding models or density functional
theory with actual localized orbitals, corresponds to point-like localized
orbitals. This, in addition to using an effective screening instead of
computing the dielectric constant, allows to obtain the interaction kernel at
very low computational cost and, thereof, the exciton spectrum as well as the
light absorption of materials.
Additional comments including restrictions and unusual features: The code
requires using at least C++11, given that it uses version-specific features.
All linear algebra routines have been delegated to the Armadillo library.

##  1 Introduction

Bound electron-hole pairs, namely excitons, are known to be largely
responsible for the most prominent features of the optical response of
semiconductors near the band edge[1], particularly for low dimensional
materials[2, 3]. This includes, of course, absorption and photoluminescence,
but also the photovoltaic response, where substantial efforts, both
experimentally and theoretically, are being made for energy-harvesting real-
life applications [4]. On the theory side, multiple ways to describe excitons
have been developed varying in accuracy and sophistication [5], from an
effective two-body description [6, 7] and configuration interaction [8, 9, 10]
to many-body perturbation theory (MBPT) [11, 12, 13] or time-dependent
techniques [14, 15, 16, 17]. MBPT itself can be purely electronic or include
electron-phonon interactions [18, 19, 20]. The current standard for exciton
calculations is GW-BSE: the GW approximation [21] is used to correct the
density functional theory (DFT) electronic band structure, specifically the
gap in insulators and semiconductors [22]. The resulting band structure is
then used to compute the exciton spectrum with the Bethe-Salpeter equation
(BSE). The accuracy of the MBPT approach has prompted the development of
several software applications for the calculation of first-principles many-
body excitations [23, 24, 25, 26].

Being first-principles calculations, one can seek quantitative agreement with
experiments, at the cost of computational time. Alternatively, one can seek
less costly, qualitative comparison through an effective description of the
interactions. Our code Xatu is intended for such purpose. While the base
electronic structure can be computed at any degree of fidelity, from the
simplest tight-binding (TB) model to the more sophisticated GW approximation,
electron-hole interactions are taken into consideration through a simplified
model where orbitals are considered to be point-like along with
phenomenological models for screening. In principle, any band structure can be
used as the starting point as long as it comes from a localized orbitals basis
code; plane waves-based calculations are out of scope since they go against
the nature of the approximation used for the interaction. Ultimately, these
two approximations result in a considerable reduction of the computational
cost.

Beyond the intrinsic speed-up coming from the calculation scheme itself, Xatu
has been written mainly in C++, and is designed to be as efficient and general
as possible while keeping its usability relatively simple. It targets a wide
range of systems in the landscape of computational tools for optical
excitations, from those that simply are out of the range of first-principles
ones because of the complexity of the unit cell, to those that require a quick
iteration, while obtaining qualitative and sometimes even quantitative
agreement with experiments.

##  2 Exciton theory

###  2.1 The Bethe-Salpeter equation

Here we review the basic aspects of our theoretical approach, highlighting the
simplifications and analogies with respect to the standarized GW-BSE method
(see e.g. [27]). From a quantum chemistry perspective, for the description of
excitons we consider the exact, non-relativistic electronic Hamiltonian of the
solid of interest:

$H=H_{0}+V=\sum_{i,j}t_{ij}c^{\dagger}_{i}c_{j}+\frac{1}{2}\sum_{i,j,k,l}V_{ijkl}c^{\dagger}_{i}c^{\dagger}_{j}c_{l}c_{k},$
(1)

where the indices include orbital and position degrees of freedom; we restrict
to basis of localized orbitals. $H_{0}$ describes the kinetic and ion-electron
interaction terms and $V$ is the electrostatic interaction between electrons.
Diagonalization of $H_{0}$ yields a Bloch eigenbasis $\ket{n\mathbb{k}}$ with
energies $\varepsilon_{n\mathbb{k}}$, which here will correspond to insulating
or semi-conducting materials. The interaction term in 1 contains

$\displaystyle V_{ijkl}$ $\displaystyle=\braket{i,j}{V}{k,l}$
$\displaystyle=\int
d\mathbb{r}d\mathbb{r}^{\prime}\varphi^{*}_{i}(\mathbb{r})\varphi^{*}_{j}(\mathbb{r}^{\prime})V(\mathbb{r},\mathbb{r}^{\prime})\varphi_{k}(\mathbb{r})\varphi_{l}(\mathbb{r}^{\prime})$
(2)

where $V(\mathbb{r},\mathbb{r}^{\prime})$ is the two-body interaction. This
can be the bare Coulomb interaction or some alternative interaction to take
into account dimensionality or screening. Since the non-interacting
Hamiltonian $H_{0}$ describes insulating materials, it is usually a good
approximation to take the ground state for the interacting Hamiltonian $H$ as
the Fermi sea:

$\ket{GS}=\prod_{n,\mathbb{k}}^{\varepsilon_{n\mathbb{k}}\leq\varepsilon_{F}}c^{\dagger}_{n\mathbb{k}}\ket{0}$
(3)

where $\ket{0}$ denotes the state with zero electrons, and $\varepsilon_{F}$
is the Fermi energy. Then, an electron-hole pair of center-of-mass momentum
$\mathbb{Q}$ between the conduction band $c$ and the valence band $v$, and
located at momentum $\mathbb{k}$ is defined as:

$\ket{v,c,\mathbb{k},\mathbb{Q}}=c^{\dagger}_{c\mathbb{k}+\mathbb{Q}}c_{v\mathbb{k}}\ket{GS}$
(4)

meaning that one electron of momentum $\mathbb{k}$ from the valence bands is
promoted to the conduction bands with momentum $\mathbb{k}+\mathbb{Q}$. Note
that even though we denote these states as electron-hole pairs, we are not
actually using hole quasiparticle operators, but simply refer to the hole as
the absence of an electron in the Fermi sea. We will stick to the electron
picture throughout this work, unless specified otherwise. These electron-hole
pairs will serve as the basis for the exciton states, $\ket{X}_{\mathbb{Q}}$:

$\displaystyle\ket{X}_{\mathbb{Q}}$
$\displaystyle=\sum_{v,c,\mathbb{k}}A_{vc}^{\mathbb{Q}}(\mathbb{k})\ket{v,c,\mathbb{k},\mathbb{Q}}$
$\displaystyle=\sum_{v,c,\mathbb{k}}A_{vc}^{\mathbb{Q}}(\mathbb{k})c^{\dagger}_{c\mathbb{k}+\mathbb{Q}}c_{v\mathbb{k}}\ket{GS}$
(5)

Therefore, the exciton is expressed as a linear combination of electron-hole
pairs over different bands and momenta. Note that $\mathbb{Q}$ serves as a
good quantum number for the exciton states, since the interaction is momentum-
conserving. The interaction only mixes electron-hole pairs with the same net
momentum, which is $\mathbb{Q}$. This can be seen by computing explicitly a
general interaction matrix element, $V_{ijkl}$.

Next, we determine the $A_{vc}^{\mathbb{Q}}(\mathbb{k})$ coefficients that
minimize the expectation value $\braket{X}{H}{X}_{\mathbb{Q}}$:

$\frac{\delta E[X]}{\delta X}=\frac{\delta}{\delta
X}\left[\frac{\braket{X}{H}{X}_{\mathbb{Q}}}{\braket{X}{X}_{\mathbb{Q}}}\right]=0$
(6)

Performing this derivative explicitly is equivalent to the problem of
diagonalizing the Hamiltonian represented in the basis of electron-hole pairs:

$\sum_{v^{\prime},c^{\prime},\mathbb{k}^{\prime}}H_{vc,v^{\prime}c^{\prime}}(\mathbb{k},\mathbb{k}^{\prime},\mathbb{Q})A_{v^{\prime}c^{\prime}}^{\mathbb{Q}}(\mathbb{k}^{\prime})=E_{X}A_{vc}^{\mathbb{Q}}(\mathbb{k})$
(7)

where
$H_{vc,v^{\prime}c^{\prime}}(\mathbb{k},\mathbb{k}^{\prime},\mathbb{Q})=\braket{v,c,\mathbb{k},\mathbb{Q}}{H}{v^{\prime},c^{\prime},\mathbb{k}^{\prime},\mathbb{Q}}$.
The expansion in electron-hole pairs of the exciton is actually an ansatz: we
obtain exact eigenstates of the Hamiltonian restricted to a partition of the
Hilbert space, $PHP$, where $P$ is a projector over the single electron-hole
pairs.

$PHP=\sum_{\begin{subarray}{c}v,c,\mathbb{k}\\\
v^{\prime},c^{\prime},\mathbb{k}^{\prime}\end{subarray}}H_{vc,v^{\prime}c^{\prime}}(\mathbb{k},\mathbb{k}^{\prime},\mathbb{Q})\ket{v,c,\mathbb{k},\mathbb{Q}}\bra{v^{\prime},c^{\prime},\mathbb{k}^{\prime},\mathbb{Q}}$
(8)

In fact, if we only consider charge-conserving excitations, we could represent
the Hamiltonian in the following way:

$H=\bigoplus^{N_{e}}_{n=0}P_{n}HP_{n}+C,\ \text{with }$ (9)

where

$\displaystyle P_{n}=\sum_{\begin{subarray}{c}\\{c_{i}\\},\\{v_{i}\\}\\\
\\{c^{\prime}_{i}\\},\\{v^{\prime}_{i}\\}\end{subarray}}\ket{\\{c_{i}\\},\\{v_{i}\\}}\bra{\\{c^{\prime}_{i}\\},\\{v^{\prime}_{i}\\}}\
$

and

$\displaystyle\ket{\\{c_{i}\\},\\{v_{i}\\}}=\prod^{n}_{i=1}c^{\dagger}_{c_{i}}\prod^{n}_{i=1}c_{v_{i}}\ket{GS}$
(10)

$N_{e}$ is the total number of electrons, $C$ the coupling between the
different excitation sectors, and $P_{n}$ is the projector over the n-th
electron-hole pairs sector. If instead of using the Bloch states from $H_{0}$
we formulate the problem in terms of the Hartree-Fock (HF) solution to 1, then
the coupling between the Fermi sea and the single-pair sector, $P_{0}HP_{1}$,
is exactly zero according to Brillouin’s theorem[28, 29]. As we will mention
later, we will assume that this always holds even when the ground state has
not been calculated in the HF approximation. The same, however, is not true
for $P_{0}HP_{2}$ or $P_{1}HP_{2}$, i.e., the interaction couples the ground
state and the one electron-hole pair sector with the two electron-hole pairs
sector. Thus, the proposed ground state and the exciton states are never exact
but approximate eigenstates. Given that the material is insulating, we expect
the coupling to be weak due to the energy differences, which justifies the
ansatz. Keeping with the exact diagonalization approach, one could try to
diagonalize the Hamiltonian including more excitation sectors. Although
possible in principle, it becomes quickly unfeasible since the Hilbert space
in many-body systems grows exponentially and, in this case, the eigenstates
would involve a mixture of excitations, losing the interpretation as a bound
electron-hole pair.

Going back to 7, we compute next the Hamiltonian matrix elements in the
$H_{0}$ basis, which are given in terms of the single particle energies and
the interaction matrix elements:

$\displaystyle
H_{vc,v^{\prime}c^{\prime}}(\mathbb{k},\mathbb{k}^{\prime},\mathbb{Q})=$
$\displaystyle\delta_{\mathbb{k}\mathbb{k}^{\prime}}\delta_{vv^{\prime}}\big{[}\varepsilon_{c\mathbb{k}+\mathbb{Q}}\delta_{cc^{\prime}}+\Sigma_{cc^{\prime}}(\mathbb{k}+\mathbb{Q},\mathbb{k}^{\prime}+\mathbb{Q})\big{]}$
$\displaystyle-\delta_{\mathbb{k}\mathbb{k}^{\prime}}\delta_{cc^{\prime}}\big{[}\varepsilon_{v\mathbb{k}}\delta_{vv^{\prime}}+\Sigma_{v^{\prime}v}(\mathbb{k}^{\prime},\mathbb{k})\big{]}-(D-X)_{vcv^{\prime}c^{\prime}}(\mathbb{k},\mathbb{k}^{\prime},\mathbb{Q})$
(11)

where

$\begin{split}D_{vc,v^{\prime}c^{\prime}}(\mathbb{k},\mathbb{k}^{\prime},\mathbb{Q})&=V_{c\mathbb{k}+\mathbb{Q},v^{\prime}\mathbb{k}^{\prime},c^{\prime}\mathbb{k}^{\prime}+\mathbb{Q},v\mathbb{k}}\\\
X_{vc,v^{\prime}c^{\prime}}(\mathbb{k},\mathbb{k}^{\prime},\mathbb{Q})&=V_{c\mathbb{k}+\mathbb{Q},v^{\prime}\mathbb{k}^{\prime},v\mathbb{k},c^{\prime}\mathbb{k}^{\prime}+\mathbb{Q}}\end{split}$
(12)

and

$\begin{split}\Sigma_{nm}(\mathbb{k},\mathbb{k}^{\prime})=\sum_{j,\mathbb{k}^{\prime\prime}}^{\text{occ}}\big{(}V_{n\mathbb{k},j\mathbb{k}^{\prime\prime},m\mathbb{k}^{\prime},j\mathbb{k}^{\prime\prime}}-V_{n\mathbb{k},j\mathbb{k}^{\prime\prime},j\mathbb{k}^{\prime\prime},m\mathbb{k}^{\prime}}\big{)}\end{split}$
(13)

$D$, $X$ correspond to the direct and exchange interactions between the
electron-hole pair, whereas $\Sigma$ is the self-energy coming from the
interaction of the electron/hole with the Fermi sea. At this point we could
obtain the exciton spectrum diagonalizing (2.1). Instead, it is more
convenient to solve first for the ground-state of (1) at the mean-field level,
i.e. in the HF approximation [30]. If we now write (2.1) in the HF band basis,
we obtain:

$\displaystyle(\varepsilon_{c\mathbb{k}+\mathbb{Q}}-\varepsilon_{v\mathbb{k}})A_{vc}^{\mathbb{Q}}(\mathbb{k})+\sum_{v^{\prime},c^{\prime},\mathbb{k}^{\prime}}K_{vc,v^{\prime}c^{\prime}}(\mathbb{k},\mathbb{k}^{\prime},$
$\displaystyle\mathbb{Q})A_{v^{\prime}c^{\prime}}^{\mathbb{Q}}(\mathbb{k}^{\prime})$
$\displaystyle=E_{X}A_{vc}^{\mathbb{Q}}(\mathbb{k})$ (14)

where $\varepsilon_{n\mathbb{k}}$ are now the HF quasiparticle energies, and
$K=-(D-X)$ is the interaction kernel. Thus, the self-energies are now
incorporated into the quasiparticle energies instead. Note that the Fermi sea
energy has been set to zero, so that exciton energies can be compared directly
with the gap of the system. This is the standard form of the Bethe-Salpeter
equation for excitons using the Tamm-Dancoff approximation (TDA) [31, 32], and
it defines the starting point for any exciton calculation. The main difference
with MBPT comes from the interaction kernel, which there involves a
dynamically screened interaction, usually in the random-phase approximation
[33, 34, 35]. The determination of the dielectric constant is a
computationally intensive task [23], which we avoid setting instead an
effective static screening.

So far we have seen that it is more convenient to pose the exciton problem in
terms of the HF basis, as it simplifies the problem and allows to decouple
excitation sectors. In practice, we do not address the problem of determining
the mean-field solution to (1). Instead, we start directly from equation (2.1)
assuming that the initial band structure, which is already known, verifies it.
Namely, for tight-binding band structures we drop the self-energy terms
assuming that we are using a HF solution. Alternatively, if the band structure
comes from DFT or MBPT (e.g. GW approximation), then we also remove the self-
energy terms since the quasiparticle energies already include self-energy
corrections (although they do not cancel exactly with those from (2.1)). Thus,
from now on we regard the starting band structure as the non-interacting
Hamiltonian $H_{0}$.

###  2.2 Interaction matrix elements

With Eq. (2.1) established, a practical expression for the interaction matrix
elements (2.1) remains to be obtained. The single-particle states, using a
basis of localized orbitals, can be written as:

$\varphi_{n\mathbb{k}}(\mathbb{r})=\frac{1}{\sqrt{N}}\displaystyle\sum_{\mathbb{R}}e^{i\mathbb{k}\cdot\mathbb{R}}\sum_{i,\alpha}C^{n\mathbb{k}}_{i\alpha}\phi_{\alpha}(\mathbb{r}-\mathbb{R}-\mathbb{t}_{i})$
(15)

where $\\{\phi_{\alpha}\\}$ denote the orbitals located at the atom $i$ of the
motif and $N$ is the number of unit cells of the system. As mentioned before,
this wavefunction may correspond to that of a tight-binding model (meaning
that the spatial nature of the orbitals is ignored and are typically
considered orthonormal), or a DFT calculation with a local orbital basis set,
which are in general non-orthogonal. While the origin of the single-particle
states can be different, for the actual calculation of the interactions we
will treat them on the same footing, approximating them as point-like
orthonormal orbitals. Depending on how we treat the interaction, different
working expressions for the matrix elements can be obtained. For instance, we
address first the direct term, which is given by:

$\displaystyle
D_{vc,v^{\prime}c^{\prime}}(\mathbb{k},\mathbb{k}^{\prime},\mathbb{Q})$
$\displaystyle=\int
d\mathbb{r}d\mathbb{r}^{\prime}\varphi^{*}_{c\mathbb{k}+\mathbb{Q}}(\mathbb{r})\varphi^{*}_{v^{\prime}\mathbb{k^{\prime}}}(\mathbb{r}^{\prime})V(\mathbb{r},\mathbb{r}^{\prime})\varphi_{c^{\prime}\mathbb{k^{\prime}}+\mathbb{Q}}(\mathbb{r})\varphi_{v\mathbb{k}}(\mathbb{r}^{\prime})$
(16)

We substitute the single-particle Bloch states (15) in Eq. (16). Expanding
each term, we end up having to evaluate the same four-body integral, but now
between the orbitals that compose each state:

$\int
d\mathbb{r}d\mathbb{r}^{\prime}\phi^{*}_{\alpha}(\mathbb{r})\phi^{*}_{\beta}(\mathbb{r}^{\prime})V(\mathbb{r},\mathbb{r}^{\prime})\phi_{\gamma}(\mathbb{r})\phi_{\delta}(\mathbb{r}^{\prime})$
(17)

At this point, there are two ways to compute the present four-body integral:
we can evaluate directly the interaction in real space, or, instead, use its
Fourier series to work in reciprocal space. In both cases we consider point-
like orbitals centered at $\mathbb{R}+\mathbb{t}_{i}$:

$\phi_{\alpha}(\mathbb{r}-\mathbb{R}-\mathbb{t}_{i})\phi_{\beta}(\mathbb{r}-\mathbb{R^{\prime}}-\mathbb{t}_{j})\approx\delta_{\alpha\beta}\delta(\mathbb{r}-\mathbb{R}-\mathbb{t}_{i})\delta_{ij}\delta_{\textbf{R},\textbf{R}^{\prime}}.$
(18)

Integrating in real space, after simplifying the resulting deltas, we obtain
the following expression for the direct term $D$:

$\displaystyle
D_{vc,v^{\prime}c^{\prime}}(\mathbb{k},\mathbb{k}^{\prime},\mathbb{Q})$
$\displaystyle=\frac{1}{N}\sum_{ij}\sum_{\alpha\beta}(C_{i\alpha}^{c\mathbb{k}+\mathbb{Q}})^{*}(C_{j\beta}^{v^{\prime}\mathbb{k}^{\prime}})^{*}C_{i\alpha}^{c^{\prime}\mathbb{k}^{\prime}+\mathbb{Q}}C_{j\beta}^{v\mathbb{k}}V_{ij}(\mathbb{k}^{\prime}-\mathbb{k})$
(19)

where

$V_{ij}(\mathbb{k}^{\prime}-\mathbb{k})=\sum_{\mathbb{R}}e^{i(\mathbb{k}^{\prime}-\mathbb{k})\mathbb{R}}V(\mathbb{R}-(\mathbb{t}_{j}-\mathbb{t}_{i})).$
(20)

Here $V_{ij}(\mathbb{k}^{\prime}-\mathbb{k})$ can be regarded as a lattice
Fourier transform centered at $\mathbb{t}_{j}-\mathbb{t}_{i}$. Since it is
defined as a sum over lattice vectors and not an integral, one cannot use the
shift property from the Fourier transform. Attempting to do so would result in
breaking the spatial symmetries of the Hamiltonian. Then, the direct term can
be interpreted as the weighted average of the Fourier transform of the
interaction between the electron and the hole, over all positions and
orbitals. The exchange term is computed analogously:

$\displaystyle
X_{vc,v^{\prime}c^{\prime}}(\mathbb{k},\mathbb{k}^{\prime},\mathbb{Q})$
$\displaystyle=\frac{1}{N}\sum_{ij}\sum_{\alpha\beta}(C_{i\alpha}^{c\mathbb{k}+\mathbb{Q}})^{*}(C_{j\beta}^{v^{\prime}\mathbb{k}^{\prime}})^{*}C_{i\alpha}^{v\mathbb{k}}C_{j\beta}^{c^{\prime}\mathbb{k}^{\prime}+\mathbb{Q}}V_{ij}(\mathbb{Q})$
(21)

In case that there is only one atom in the motif, then expressions (2.2),
(2.2) simplify even further since the interaction decouples from the tight-
binding coefficients, yielding:

$\begin{split}&D_{vc,v^{\prime}c^{\prime}}(\mathbb{k},\mathbb{k}^{\prime},\mathbb{Q})=\frac{1}{N}V(\mathbb{k}^{\prime}-\mathbb{k})(U^{\dagger}_{\mathbb{k}+\mathbb{Q}}U_{\mathbb{k}^{\prime}+\mathbb{Q}})_{cc^{\prime}}(U^{\dagger}_{\mathbb{k}}U_{\mathbb{k}^{\prime}})_{v^{\prime}v}\\\
&X_{vc,v^{\prime}c^{\prime}}(\mathbb{k},\mathbb{k}^{\prime},\mathbb{Q})=\frac{1}{N}V(\mathbb{Q})(U^{\dagger}_{\mathbb{k}+\mathbb{Q}}U_{\mathbb{k}})_{cv}(U^{\dagger}_{\mathbb{k}^{\prime}}U_{\mathbb{k}^{\prime}+\mathbb{Q}})_{v^{\prime}c^{\prime}}\end{split}$
(22)

where $U_{\mathbb{k}}$ is the unitary matrix that diagonalizes the Bloch
Hamiltonian $H(\mathbb{k})$ [36]. The evaluation of these expressions is much
faster than the corresponding ones (2.2) and (2.2) for a general case.
Additionally, for $\mathbb{Q}=0$, the exchange term (2.2) becomes exactly
zero, which is not true in general, although it is usually neglected. As
mentioned before, for DFT band structures we evaluate the interaction using
the same point-like approximation, performing first a Löwdin orthogonalization
of the basis. This allows to improve the TB descriptions, incorporating fine
details to the quasiparticle dispersion along the BZ. In such treatments, our
interaction matrix elements are an approximation to the true ones involving
ab-initio orbitals. Given that in DFT the orbitals are known (e.g. gaussian-
type basis in the CRYSTAL [37] code), one could, in principle, evaluate the
integrals (17) exactly for a closer ab-initio calculation of excitons.

The previous calculation corresponds to the evaluation of the interaction
matrix elements in real space. An alternative approach consists of writing the
interaction as its Fourier series before evaluating (17) [17, 16]:

$V(\mathbb{r}-\mathbb{r^{\prime}})=\frac{1}{N}\sum_{\mathbb{q}}V(\mathbb{q})e^{i\mathbb{q}\cdot(\mathbb{r}-\mathbb{r^{\prime}})}$
(23)

where

$V(\mathbb{q})=\frac{1}{V_{cell}}\int_{\Omega}V(\mathbb{r})e^{-i\mathbb{q}\cdot\mathbb{r}}d\mathbb{r}$
(24)

$\Omega=NV_{cell}$ denotes the volume of the crystal. Usually, one takes
$\Omega\to\infty$ meaning that we can evaluate the integral analytically, this
is, $V(\mathbb{q})$ becomes the Fourier transform of the potential. Note,
however, that $\mathbb{q}$ is not restricted to the first Brillouin Zone (BZ),
and $V(\mathbb{q})$ is not periodic in the BZ. Therefore, in principle one has
to sum over $\mathbb{q}\in\text{BZ}$, but also over reciprocal vectors
$\mathbb{G}$, i.e.:

$V(\mathbb{r}-\mathbb{r^{\prime}})=\frac{1}{N}\sum_{\mathbb{q}\in\text{BZ}}\sum_{\mathbb{G}}V(\mathbb{q}+\mathbb{G})e^{i(\mathbb{q}+\mathbb{G})\cdot(\mathbb{r}-\mathbb{r^{\prime}})}$
(25)

The evaluation of the integral is done in the same way, although in this case
there is a plane wave instead of the electrostatic interaction. This approach
is particularly useful when using a plane wave basis, since it allows to
evaluate the four-body integrals exactly without need for approximation (18).
The interaction matrix elements $D$, $X$ are now given by:

$\begin{split}&D_{vc,v^{\prime}c^{\prime}}(\mathbb{k},\mathbb{k}^{\prime},\mathbb{Q})=\frac{1}{N}\sum_{\mathbb{G}}V(\mathbb{k}-\mathbb{k}^{\prime}+\mathbb{G})I^{\mathbb{G}}_{c\mathbb{k}+\mathbb{Q},c^{\prime}\mathbb{k}^{\prime}+\mathbb{Q}}(I^{\mathbb{G}}_{v\mathbb{k},v^{\prime}\mathbb{k}^{\prime}})^{*}\\\
&X_{vc,v^{\prime}c^{\prime}}(\mathbb{k},\mathbb{k}^{\prime},\mathbb{Q})=\frac{1}{N}\sum_{\mathbb{G}}V(\mathbb{Q}+\mathbb{G})I^{\mathbb{G}}_{c\mathbb{k}+\mathbb{Q},v\mathbb{k}}(I^{\mathbb{G}}_{c^{\prime}\mathbb{k^{\prime}}+\mathbb{Q},v^{\prime}\mathbb{k}^{\prime}})^{*}\end{split}$
(26)

where

$I^{\mathbb{G}}_{n\mathbb{k},m\mathbb{k}^{\prime}}=\sum_{i\alpha}(C_{i\alpha}^{n\mathbb{k}})^{*}C_{i\alpha}^{m\mathbb{k^{\prime}}}e^{i(\mathbb{k}-\mathbb{k}^{\prime}+\mathbb{G})\cdot\mathbb{t}_{i}}$
(27)

Usually $V(\mathbb{q})$ decays fast enough so it suffices to sum only over
$\mathbb{G}=\mathbb{0}$ for the excitons to converge in energy. Xatu allows to
use the interactions evaluated in real-space (expressions (2.2), (2.2)), or in
reciprocal space (expressions (LABEL:reciprocal_interaction)). They are
benchmarked in section 4; by default we use the interactions in real-space
since the calculation converges faster with the number of points in the BZ
mesh, it can be used rigorously for finite systems such as ribbons, and the
current implementation performs on par with the reciprocal one.

Once the interaction kernel is determined, Eq. (2.1) can be solved to obtain
the exciton energies and wavefunctions, i.e., the coefficients
$A_{vc}^{\mathbb{Q}}(\mathbb{k})$. These can be used to compute different
quantities. For instance, given that the exciton is written as a linear
combination of electron-hole pairs with well-defined $\mathbb{k}$ quantum
number, we can define the probability density of finding the exciton in a
specific pair in $\mathbb{k}$-space as:

$|\psi_{X}(\mathbb{k})|^{2}=\sum_{v,c}|A_{vc}^{\mathbb{Q}}(\mathbb{k})|^{2}$
(28)

which is the straightforward definition since all electron-hole pairs are
orthonormal to each other.

###  2.3 Spinful excitons

If the single-particle basis includes spin, one can also compute the expected
value of the total spin of the exciton, $S_{z}^{T}=S_{z}^{e}+S_{z}^{h}$. Given
that we are using a fully electronic description of the exciton, we need to
specify the electrons whose spin we want to measure. To this purpose, we write
the total spin operator in second quantization as:

$S^{T}_{z}=\sum_{c^{\prime},c,\mathbb{k}}\sigma_{c^{\prime}\mathbb{k}+\mathbb{Q},c\mathbb{k}+\mathbb{Q}}c^{\dagger}_{c^{\prime}\mathbb{k}+\mathbb{Q}}c_{c\mathbb{k}+\mathbb{Q}}-\sum_{v^{\prime},v,\mathbb{k}}\sigma_{v^{\prime}\mathbb{k},v\mathbb{k}}c_{v^{\prime}\mathbb{k}}c^{\dagger}_{v\mathbb{k}}$
(29)

where $\sigma_{nm}=\braket{n}{S_{z}}{m}$. The labels
$c,c^{\prime},v,v^{\prime}$ refer exclusively to the conduction and valence
bands used in the definition of the excitons. Note that the second term, which
corresponds to the spin of the hole, has a minus sign. This is because holes,
when described as quasiparticles, have opposite momentum and spin of the
corresponding electronic state, i.e.
$h^{\dagger}_{n,-\mathbb{k},-\sigma}=(-1)^{\sigma}c_{n\mathbb{k}\sigma}$, for
states below the Fermi energy, $\varepsilon_{n\mathbb{k}}<\varepsilon_{F}$
[38]. These $h$ operators describe creation/annihilation of holes in terms of
their electronic counterpart. Although we keep $\mathbb{k}$ the same (since we
are still in the electronic picture), we already incorporate this minus sign
to give a correct description of the total spin of the exciton. As we will see
later, this sign change is also necessary to retrieve the known singlet and
triplet states when summing angular momentum. The two pictures are equivalent,
and all the previous calculations can be reproduced in the electron-hole
picture. The expected value of the total spin is then given by:

$\displaystyle\braket{X}{S^{T}_{z}}{X}=\left[\sum_{v,c,\mathbb{k}}\sum_{c^{\prime}}A^{\mathbb{Q}}_{vc}\right.$
$\displaystyle(\mathbb{k})(A^{\mathbb{Q}}_{vc^{\prime}}(\mathbb{k}))^{*}\sigma_{c\mathbb{k}+\mathbb{Q},c^{\prime}\mathbb{k}+\mathbb{Q}}$
$\displaystyle\left.-\sum_{v^{\prime}}A^{\mathbb{Q}}_{vc}(\mathbb{k})(A^{\mathbb{Q}}_{v^{\prime}c}(\mathbb{k}))^{*}\sigma_{v\mathbb{k},v^{\prime}\mathbb{k}}\right]$
(30)

If $[H_{0},S_{z}]=0$, then the spin projection $S_{z}$ is also a good quantum
number for the Bloch states. Therefore, they can be written now as
$\ket{n\mathbb{k}\sigma}$, or in real space as
$\varphi_{n\mathbb{k}}(\mathbb{r})\chi_{\sigma}$, where $\chi_{\sigma}$
denotes the spin part of the state. This means that the spin operator $S_{z}$
is diagonal, $\sigma_{nm}=\sigma_{n}\delta_{nm}$, which allows us to simplify
expression (2.3):

$\braket{S_{z}^{T}}=\sum_{v,c,\mathbb{k}}|A^{\mathbb{Q}}_{vc}(\mathbb{k})|^{2}(\sigma_{c}-\sigma_{v})$
(31)

Another consequence of having the spin well-defined is that it also propagates
to the electron-hole pairs that serve as a basis for the exciton states, i.e.
$\ket{\tilde{v},\tilde{c},\mathbb{k},\mathbb{Q}}=c^{\dagger}_{\tilde{c}\mathbb{k}+\mathbb{Q}}c_{\tilde{v}\mathbb{k}}\ket{GS}$,
where $\tilde{v}=(v,\sigma_{v})$, $\tilde{c}=(c,\sigma_{c})$. In principle, we
allow the spin of the electron and the hole to be different,
$\sigma_{c}\neq\sigma_{v}$. Taking into account the spin in the computation of
the interaction matrix elements, we arrive at constraints on which electron-
hole pairs interact. Then the direct and exchange terms read:

$\begin{split}D_{\tilde{v}\tilde{c},\tilde{v}^{\prime}\tilde{c}^{\prime}}(\mathbb{k},\mathbb{k}^{\prime},\mathbb{Q})&=\delta_{\sigma_{c}\sigma_{c^{\prime}}}\delta_{\sigma_{v}\sigma_{v^{\prime}}}D_{vc,v^{\prime}c^{\prime}}(\mathbb{k},\mathbb{k}^{\prime},\mathbb{Q})\\\
X_{\tilde{v}\tilde{c},\tilde{v}^{\prime}\tilde{c}^{\prime}}(\mathbb{k},\mathbb{k}^{\prime},\mathbb{Q})&=\delta_{\sigma_{c}\sigma_{v}}\delta_{\sigma_{v^{\prime}}\sigma_{c^{\prime}}}X_{vc,v^{\prime}c^{\prime}}(\mathbb{k},\mathbb{k}^{\prime},\mathbb{Q})\end{split}$
(32)

which can be directly obtained by substituting the single-particle states,
since the spin part is not mixed with the orbital part of the states (i.e.
$\ket{n\mathbb{k}\sigma}=\ket{n\mathbb{k}}\otimes\ket{\sigma}$). At this
point, we can arrange the electron-hole pairs into four groups depending on
their spin:

$\\{\ket{++},\ket{--},\ket{+-},\ket{-+}\\}_{e}=\\{\ket{\sigma_{c}\sigma_{v}}\\}_{e}$

The $e$ subindex is used to denote that this corresponds to the electronic
picture. Then the Hamiltonian represented in terms of the spin groups, and
taking into account (32) becomes:

$H=\begin{pmatrix}H_{0}-D+X&X&0&0\\\ X&H_{0}-D+X&0&0\\\ 0&0&H_{0}-D&0\\\
0&0&0&H_{0}-D\end{pmatrix}$ (33)

where $H_{0}$, $D$, $X$ are blocks which include matrix elements corresponding
to different electron-hole pairs but same spin group. If we now take into
account that the hole in its quasiparticle representation must have spin
opposite of that of the electron vacancy, then our states are
$\\{\ket{+-},\ket{-+},\ket{++},\ket{--}\\}_{eh}$, where $eh$ denotes electron-
hole picture. Therefore, the exciton spectrum would be composed of groups of
three triplet states and one singlet state, as when adding angular momenta. If
instead we turn off the exchange interaction, then every state should have at
least four-fold degeneracy. Any additional degeneracy would come from spatial
symmetries of the Hamiltonian, in particular from the irreducible
representations of the little group at $\mathbb{Q}$ (see A).

###  2.4 Real-space wavefunction

Plotting the probability density (28) is useful to extract some information
about the exciton such as the wavefunction type ($s$, $p$, etc, following the
hydrogenic model). The same can be argued for its real-space wavefunction,
$\psi_{X}(\mathbb{r}_{e},\mathbb{r}_{h})$. However, obtaining it is not as
straightforward as the $\mathbb{k}$ wavefunction. To do so, first we define
the field operators as:

$\psi^{\dagger}(\mathbb{r})=\sum_{n\mathbb{k}}\varphi^{*}_{n\mathbb{k}}(\mathbb{r})c^{\dagger}_{n\mathbb{k}},\
\psi(\mathbb{r})=\sum_{n\mathbb{k}}\varphi_{n\mathbb{k}}(\mathbb{r})c_{n\mathbb{k}}$
(34)

where $\varphi_{n\mathbb{k}}(\mathbb{r})$ are the single-particle states in
coordinate representation. Then, we can define the amplitude or real space
wavefunction of the exciton in the following way:

$\psi_{X}(\mathbb{r}_{e},\mathbb{r}_{h})=\braket{GS}{\psi(\mathbb{r}_{e})\psi^{\dagger}(\mathbb{r}_{h})}{X}$
(35)

This definition is motivated by the fact that
$\varphi_{n\mathbb{k}}(\mathbb{r})=\braket{GS}{\psi(\mathbb{r})}{n\mathbb{k}}$.
Before computing the amplitude, it is convenient to switch to the electron-
hole picture. The field operator written in terms of electron and hole
operators is:

$\psi(\mathbb{r})=\sum_{c\mathbb{k}}\varphi_{c\mathbb{k}}(\mathbb{r})c_{c\mathbb{k}}+\sum_{v\mathbb{k}}\varphi_{v\mathbb{k}}(\mathbb{r})h^{\dagger}_{v\mathbb{-k}}\equiv\psi_{e}(\mathbb{r})+\psi^{\dagger}_{h}(\mathbb{r})$
(36)

where $\psi_{e}(\mathbb{r})$, $\psi_{h}(\mathbb{r})$ are the annihilation
field operator for the electrons, and holes respectively. Since we are
switching from the electronic to the electron-hole picture, the same has to be
done for the exciton state,
$\ket{X}=\sum_{v,c,\mathbb{k}}A_{vc}^{\mathbb{Q}}(\mathbb{k})c^{\dagger}_{c\mathbb{k}+\mathbb{Q}}h^{\dagger}_{v,-\mathbb{k}}\ket{0}$.
Evaluating the exciton amplitude in terms of the electron and hole field
operators, we obtain:

$\displaystyle\psi_{X}(\mathbb{r}_{e},\mathbb{r}_{h})$
$\displaystyle=\braket{GS}{\psi_{e}(\mathbb{r}_{e})\psi_{h}(\mathbb{r}_{h})}{X}$
$\displaystyle=\sum_{v,c,\mathbb{k}}A_{vc}^{\mathbb{Q}}(\mathbb{k})\varphi_{c\mathbb{k}+\mathbb{Q}}(\mathbb{r}_{e})\varphi^{*}_{v\mathbb{k}}(\mathbb{r}_{h})$
(37)

To obtain the first equality note that there are four cross terms containing
electron and hole field operators. Two of them are zero, since they they move
around the electron or the hole [e.g.
$\psi_{e}(\mathbb{r}_{e})\psi^{\dagger}_{e}(\mathbb{r}_{h})$], meaning that
the final state is still orthonormal to the ground state. There is a third
term consisting on creation of an electron and a hole,
$\psi^{\dagger}_{e}(\mathbb{r}_{e})\psi^{\dagger}_{h}(\mathbb{r}_{h})$. This
term is also zero because we assume that our ground state is the Fermi sea,
meaning that it does not contain excited electrons. If this were the case,
then the exciton could also consist on deexcitations or antiresonant
transitions. This is known as the Tamm-Dancoff approximation, and it is also
usually present in GW-BSE. To obtain the final expression for the
wavefunction, it remains to substitute the expression of the field operators.
One recovers the electron-hole pairs states of the exciton basis (up to a sign
from operator permutation), and from orthonormality it results in expression
(2.4).

At this point, to be able to plot the exciton real space wavefunction, we
still need to evaluate (2.4) in terms of the single-particle states
$\varphi_{n\mathbb{k}}(r)$. Since the exciton wavefunction depends on both the
position of the electron and the hole, first we need to fix the position of
either of them to be able to plot the wavefunction. Since we assume the
orbitals are point-like, both the electron and hole can only be localized at
the atomic positions, so we will evaluate the wavefunction and the probability
density at these points only.

We set the electron to be located at cell $\mathbb{R}_{e}$ and atom
$\mathbb{t}_{m}$ of the motif, $\mathbb{r}_{e}=\mathbb{R}_{e}+\mathbb{t}_{m}$,
while the hole is at position $\mathbb{r}_{h}=\mathbb{R}_{h}+\mathbb{t}_{n}$.
Using the point-like approximation (18), the probability density of finding
the electron at a given position with the hole fixed reads:

$|\psi_{X}(\mathbb{R}_{e}+\mathbb{t}_{m},\mathbb{R}_{h}+\mathbb{t}_{n})|^{2}=\sum_{\alpha\beta}|\psi^{\alpha\beta}_{X}(\mathbb{R}_{e}+\mathbb{t}_{m},\mathbb{R}_{h}+\mathbb{t}_{n})|^{2}$
(38)

where

$\displaystyle|\psi_{X}^{\alpha\beta}(\mathbb{R}_{e}+\mathbb{t}_{m},$
$\displaystyle\mathbb{R}_{h}+\mathbb{t}_{n})|^{2}=$
$\displaystyle\frac{1}{N^{2}}\sum_{v,c,\mathbb{k}}\sum_{v^{\prime},c^{\prime},\mathbb{k}^{\prime}}A_{vc}^{\mathbb{Q}}(\mathbb{k})(A_{v^{\prime}c^{\prime}}^{\mathbb{Q}}(\mathbb{k}^{\prime}))^{*}e^{i(\mathbb{k}-\mathbb{k}^{\prime})\cdot(\mathbb{R_{e}}-\mathbb{R_{h}})}$
$\displaystyle\cdot
C^{c,\mathbb{k}+\mathbb{Q}}_{m\alpha}(C^{c^{\prime},\mathbb{k}^{\prime}+\mathbb{Q}}_{m\alpha})^{*}(C^{v,\mathbb{k}}_{n\beta})^{*}C^{v^{\prime},\mathbb{k}^{\prime}}_{n\beta}$
(39)

For both the reciprocal and the real-space probability densities, one could
expect them to have the symmetries of the crystal, since $[H,C]=0$, where $C$
is any symmetry operator from the space group. However, if the states are
degenerate, then they are not necessarily eigenstates of the symmetry
operators and consequently the associated densities will not be invariant
under symmetry transformations. Still, in this case it is possible to define a
probability density that preserves the symmetry of the crystal for each
degenerate subspace:

$|\psi_{X}(\mathbb{r},\mathbb{r}_{h})|^{2}=\sum_{n}|\psi_{X}^{(n)}(\mathbb{r},\mathbb{r}_{h})|^{2}$
(40)

where the index $n$ runs over exciton states degenerate in energy. An
analogous expression holds for the $\mathbb{k}$ wavefunction. It is always
good practice to check that the resulting probability densities preserve the
symmetry of the crystal to ensure that the calculation was done correctly. The
proof for the invariance under symmetry operations of (40) is given in A.

###  2.5 Optical conductivity and light absorption

As an example of a post-processing calculation, we investigate here the
interaction of the material with an incident linearly-polarized electric
field. We elaborate below how to compute the optical response by means of the
exciton eigenfunctions.

For a sufficiently low-intensity and linearly-polarized homogeneus electric
pulse, $\mathbb{\varepsilon}(t)$, the induced current per unit frequency in
the bulk of the material can be written
$J_{a}=\sum_{\beta}\sigma_{ab}(\omega)\tilde{\varepsilon}_{b}(\omega)$, where
the linear optical conductivity reads [39]

$\begin{split}\sigma_{ab}(\omega)=&\frac{\pi
e^{2}\hbar}{V}\sum_{k}^{N_{X}}\frac{1}{E_{k}}\bigg{[}V_{k}^{a}(V_{k}^{b})^{\ast}\bigg{]}\delta(\hbar\omega-
E_{k})\end{split}$ (41)

Here, $N_{X}$ is the number of exciton states, $E_{k}$ is the energy of the
$k$-th excited state, $V_{k}^{a}=\braket{GS}{\hat{v}^{a}}{X_{k}}$ is the
velocity matrix element (VME) transition to the ground state and $V$ is the
volume of the solid under periodic boundary conditions. In the equation above,
only excitons with $\mathbb{Q}=0$ are considered, as finite momentum excitons
cannot be achieved by light incidence. Thus, we drop $\mathbb{Q}$ from the
notation, and instead specify the excitation index $k$ in the exciton
coefficients, $A^{k}_{vc}(\mathbb{k})$. The exciton VME is found to be

$V_{k}^{a}=\sum_{cv\mathbb{k}}A_{vc}^{k}(\mathbb{k})v^{a}_{vc}(\mathbb{k}),$
(42)

where
$v^{a}_{vc}(\mathbb{k})\equiv\braket{v\mathbb{k}}{\hat{v}^{a}}{c\mathbb{k}}=i\hbar^{-1}\braket{v\mathbb{k}}{[H_{0},\hat{r}^{a}]}{c\mathbb{k}}$
($H_{0}$ is the non-interacting or mean-field Hamiltonian). With light
polarized along the $a$ direction, an exciton is dark (or bright) if
$V_{k}^{a}=0$ ($V_{k}^{a}\neq 0$). In general, the brightness of an exciton
and its contribution to Eq. (41) is dictated by the selection rules of
$A_{vc}^{k}(\mathbb{k})$ and $v^{a}_{vc}(\mathbb{k})$ over the Brillouin zone.
The calculation of VMEs has to be worked out taking into account the
underlying local basis of our approach. It is found [40, 41] (we simplify the
notation here by doing $i\alpha\rightarrow\alpha$ for the rest of the section)

$\begin{split}&\braket{n\mathbb{k}}{\hat{\mathbb{v}}}{n^{\prime}\mathbb{k}}=\\\
&\sum_{\alpha\alpha^{\prime}}(C^{n\mathbb{k}}_{\alpha})^{*}C^{n^{\prime}\mathbb{k}}_{\alpha^{\prime}}\nabla_{\mathbb{k}}H_{\alpha\alpha^{\prime}}(\mathbb{k}),\\\
&+i\sum_{\alpha\alpha^{\prime}}(C^{n\mathbb{k}}_{\alpha})^{*}C^{n^{\prime}\mathbb{k}}_{\alpha^{\prime}}\Big{[}\epsilon_{n}(\mathbb{k})\mathbb{\xi}_{\alpha\alpha^{\prime}}(\mathbb{k})-\epsilon_{n^{\prime}}(\mathbb{k})\mathbb{\xi}_{\alpha^{\prime}\alpha}^{\ast}(\mathbb{k})\Big{]}.\end{split}$
(43)

with
$\mathbb{\xi}_{\alpha\alpha^{\prime}}(\mathbb{k})=i\braket{u_{\alpha\mathbb{k}}}{\nabla_{\mathbb{k}}u_{\alpha^{\prime}\mathbb{k}}}$
the Berry connection between Bloch basis states. After some algebra, it reads

$\mathbb{\xi}_{\alpha\alpha^{\prime}}(\mathbb{k})=\sum_{\mathbb{R}}e^{i\mathbb{k}\cdot\mathbb{R}}\braket{\alpha\mathbb{0}}{\hat{\mathbb{r}}}{\alpha^{\prime}\mathbb{R}}+i\nabla_{\mathbb{k}}S_{\alpha\alpha^{\prime}}(\mathbb{k}).$
(44)

In the case of an underlying non-orthonormal local orbital basis, the overlap
matrix
$S_{\alpha\alpha^{\prime}}(\mathbb{k})\equiv\braket{\alpha\mathbb{k}}{\alpha^{\prime}\mathbb{k}}$
is accounted and makes the Berry connection above non-hermitian. Instead, one
has
$\mathbb{\xi}_{\alpha\alpha^{\prime}}(\mathbb{k})=\mathbb{\xi}_{\alpha^{\prime}\alpha}^{\ast}(\mathbb{k})+i\nabla_{\mathbb{k}}S_{\alpha\alpha^{\prime}}(\mathbb{k})$.
Eq. (43) allows to evaluate the optical matrix elements by means of the non-
interacting hamiltonian plus position matrix elements between the local
orbitals.

In the case of an orthogonal basis set, as in tight-binding models, the
overlap matrix is an unitary matrix at all points of the Brillouin zone. In
this case VMEs read

$v^{a}_{vc}(\mathbb{k})=\sum_{\alpha\alpha^{\prime}}(C^{n\mathbb{k}}_{\alpha})^{*}C^{n^{\prime}\mathbb{k}}_{\alpha^{\prime}}\Bigg{[}\frac{\partial
H_{\alpha\alpha^{\prime}}(\mathbb{k})}{\partial
k_{a}}+iH_{\alpha\alpha^{\prime}}(\mathbb{k})(t^{a}_{\alpha^{\prime}}-t^{a}_{\alpha})\Bigg{]}$
(45)

This expression is sometimes known as the “diagonal tight-binding
approximation (TBA)” in ab-initio calculations involving the maximally-
localized Wannier functions [42, 43], where the inter-orbital position matrix
elements are discarded. Eq. (41) can thus be evaluated and is implemented in
our code.

Additionaly, Eq. (41) can be compared with the frequency dependent expression
for its non-interacting counterpart. In the limit of no correlations, it
reduces to

$\displaystyle\sigma_{ab}(\omega)=\frac{\pi
e^{2}\hbar}{V}\sum_{cv\mathbb{k}}\frac{1}{\varepsilon_{c\mathbb{k}}-\varepsilon_{v\mathbb{k}}}$
$\displaystyle\bigg{[}v_{cv}^{a}(\mathbb{k})v_{vc}^{b}(\mathbb{k})\bigg{]}$
$\displaystyle\cdot\delta(\hbar\omega-[\varepsilon_{c\mathbb{k}}-\varepsilon_{v\mathbb{k}}])$
(46)

From the frequency-dependent optical conductivity one can obtain related
quantities of interest. For instance, the ratio of absorbed incident flux per
unit frequency and unit length (considering vaccum surroundings) is [44]
$S(\omega)=\sigma(\omega)/c\epsilon_{0}$, also called absorbance. Note that
this quantity is ill-defined for 2D lattice systems. In such case, all the
absorbance is assumed to occur at $z=0$ reference plane of the material.

##  3 Implementation

The programming languages of choice for the implementation of the exciton
theory were C++ and Fortran, which are the usual options for heavy numerical
computations. In the case of C++, to facilitate manipulation of matrices we
use the library Armadillo [45], on top of the usual libraries for linear
algebra (BLAS and LAPACK). The core of the code was written in C++, except the
post-diagonalization calculation of the optical conductivity being written in
Fortran. This routine is wrapped inside the C++ library. The software was
designed with a hybrid approach in mind: previous packages such as DFT codes
require the preparation of input files, which are then fed to the program and
result in some output files which may be post-processed to extract
information. We propose to use the same scheme, i.e. to prepare an input file
which describes the system where we want to compute the excitons (namely the
Hamiltonian $H_{0}$), and another one with the description of the excitons
(participating bands, $\mathbb{k}$ mesh, etc). However, there is an
alternative usage, which is employing directly the exciton API defined to
built the program. This a common approach, where one builds a library to
expose some functionality to the user (e.g. Python libraries). Therefore, one
can define some system and script the computation of excitons using the API.
This is advised whenever we are interested in performing some manipulation of
the excitons, and not only obtaining the spectrum or the absorption. There is
a third approach, consisting on using the system files to leverage the
definition of the system to other programs (e.g. DFT), and then use the API
instead of the exciton configuration file. The different forms to use the code
will be reviewed in B. The CLI option parsing has been done using the header-
only library TCLAP, which is distributed with this package.

###  3.1 Complexity analysis

Next we will discuss the numerical implementation of the exciton computation
and related quantities. Solving the Bethe-Salpeter equation (2.1) amounts to
diagonalizing the corresponding matrix $PHP$. Diagonalization is done using
the standard linear algebra libraries, meaning that the main problem is
constructing $PHP$ as fast as possible. Consider a system formed by $N$ unit
cells in total (meaning $\sqrt{N}$ along each Bravais vector for a two
dimensional system). To treat the interaction rigorously, one has to compute
the excitons on a BZ mesh with the same number of $\mathbb{k}$ points as unit
cells, due to the periodic boundary conditions. Therefore, one has to compute
$N^{2}$ matrix elements, and each of them requires computing the lattice
Fourier transform, which involves summations over the $N$ unit cells. This has
to be done for all possible band pairs $B$, so a naive implementation of (2.1)
would have $\mathcal{O}(N^{3}B^{2})$ time complexity, on par with matrix
diagonalization algorithms. Note that each interaction matrix element also
requires knowing the tight-binding coefficients
$\\{C^{n\mathbb{k}}_{i\alpha}\\}$. If the dimension of the Bloch Hamiltonian
is $M$, then diagonalizing on the fly the system for each element of the BSE
would result in time $\mathcal{O}(N^{2}B^{2}(N+M^{3}))$.

The easiest way to reduce the time complexity of the BSE construction is to
increase the space complexity, i.e. to precompute and store quantities that
appear multiple times, instead of computing them on the fly. This can be done
for the Bloch Hamiltonian eigenvectors. Before constructing $PHP$, we
diagonalize $H(\mathbb{k})$ $\forall\mathbb{k}\in\text{BZ}$, and store the
eigenvectors. At this point, if were to store all eigenvectors, the spatial
complexity would go from $\mathcal{O}(1)$ to $\mathcal{O}(NM^{2})$. Since we
only need the eigenvectors corresponding to the bands that participate in the
exciton formation, it suffices to store only those, meaning that the spatial
complexity would be $\mathcal{O}(NMB)$, i.e. we have to store $N$ matrices of
size $M\times B$. Accessing directly the eigenvectors results in a time
complexity of $\mathcal{O}(N^{3}B^{2}+NM^{3})$.

The same could be done for the lattice Fourier transform
$V_{ij}(\mathbb{k}-\mathbb{k}^{\prime})$. Since it depends on the difference
between between two $\mathbb{k}$ points, we could simply store $V_{ij}$ for
each pair of $\mathbb{k}$ points. This implies high spatial complexity
$\mathcal{O}(N^{2})$, but overall it does not report any speed advantage,
since precomputing this would be of order $\mathcal{O}(N^{3})$. However, it is
possible to reduce the time cost of the algorithm: as long as the $\mathbb{k}$
point mesh covers the whole BZ uniformly (as given by Monkhorst-Pack), then we
can map the $\mathbb{k}$ point difference back to a single $\mathbb{k}$ point
using the periodicity of $V_{ij}(\mathbb{k}-\mathbb{k}^{\prime})$:

$\displaystyle\forall\mathbb{k},\mathbb{k}^{\prime}\in\text{BZ,
}\exists\mathbb{G}\in\text{Reciprocal lattice, }$
$\displaystyle\mathbb{k}^{\prime\prime}\in\text{BZ}$ $\displaystyle\text{ s.t.
}\mathbb{k}-\mathbb{k}^{\prime}=\mathbb{G}+\mathbb{k}^{\prime\prime}$ (47)

Therefore, it suffices to compute and store $V_{ij}(\mathbb{k})$
$\forall\mathbb{k}\in\text{BZ}$. Then, when initializing the matrix elements
of $PHP$, one has to find the vector $\mathbb{k}^{\prime\prime}$ such that it
verifies (3.1). The time complexity now is $\mathcal{O}(N^{2})$, which is a
reduction of an order of magnitude. The space complexity is also reduced,
being now $\mathcal{O}(N)$.

With this, the algorithm for determining $PHP$ has time order
$\mathcal{O}(N^{2}B^{2}+NM^{3})$, and the memory requirements are
$\mathcal{O}(N+NMB)=\mathcal{O}(NMB)$. As we will see, this allows for very
fast computation of the BSE matrix, meaning that the main bottleneck lies in
the diagonalization, as it often happens. In some cases we might be interested
in the whole spectrum, but usually it suffices to determine the lowest energy
eigenstates. To address this, the code includes a custom implementation of the
Davidson algorithm, which is suited to obtain the ground state of quantum
chemistry Hamiltonians [46].

So far the discussion was focused on how to reduce the complexity of the
algorithm, but it is equally important to comment on how to perform the actual
computation of the matrix elements. The big O notation neglects all constant
factors, which is fine for theoretical considerations, but might have a
considerable impact on the real behaviour of the code. The general strategy
followed was to vectorize all calculations to make use of the highly optimized
and parallel existing linear algebra routines. The remaining parts that do not
allow vectorization, such as matrix element initialization in $PHP$, were all
parallelized with OpenMP. Currently, all the parallelism is shared-memory and
distributed parallelism might be implemented in the future.

For instance, consider the direct interaction term which requires computing
expression (2.2). Supposed that the lattice Fourier transform of the
interaction is already computed for all motif combinations $i,j$ and for all
$\mathbb{k}$ points, we basically have to sum over tight-binding coefficients
multiplied by the interaction. Given that the Bloch eigenstates are already
stored as columns in matrices, we want to write this as matrix-vector
products. Specifically, we can use $V_{ij}$ as a bilinear form, so with a
well-defined matrix $\tilde{V}$ the direct term can be written as:

$D_{vc,v^{\prime}c^{\prime}}(\mathbb{k},\mathbb{k}^{\prime},\mathbb{Q})=C^{T}_{cc^{\prime}}\tilde{V}(\mathbb{k}^{\prime}-\mathbb{k})C_{v^{\prime}v}$
(48)

where

$\tilde{V}=V(\mathbb{k}-\mathbb{k}^{\prime})\otimes\mathbb{I}_{n},\ \
C_{nm}=C^{*}_{n}\odot C_{m}$ (49)

$\odot$ denotes element-wise array product, $C_{n}$ is the vector of
coefficients corresponding to state $\ket{n}$ and $\mathbb{I}_{n}$ denotes a
square matrix of ones of dimension $n$, $n$ being the number of orbitals per
atom. Note that this expression is only valid if all atoms have the same
number of orbitals. Otherwise, one must take into account the different number
of orbitals per chemical species when performing the Kronecker´s products. The
exchange term $X$ can be computed in an analogous way. Note that this assumes
that the order of the single-particle basis is
$\\{\ket{i}\otimes\ket{\alpha}\otimes\ket{\sigma}\\}$, i.e. for each atomic
position, we run over orbitals, and for each orbital we run over spin. This is
also relevant for the computation of the spin of the excitons, since it
follows this convention.

As we mentioned at the beginning, to compensate for the lack of screening of
the theory, one typically uses the Rytova-Keldysh potential [47, 48] instead
of the bare Coulomb potential in the context of two-dimensional materials.
However, both interactions diverge at $r=0$. We regularize this divergence by
setting $V(0)=V(a)$ [36], where $a$ denotes the lattice parameter. Currently,
the code only implements the Keldysh potential, given by:

$V(r)=\frac{e^{2}}{8\varepsilon_{0}\bar{\varepsilon}r_{0}}\left[H_{0}\left(\frac{r}{r_{0}}\right)-Y_{0}\left(\frac{r}{r_{0}}\right)\right]$
(50)

where $\bar{\varepsilon}=(\varepsilon_{m}+\varepsilon_{s})/2$, with
$\varepsilon_{s}$, $\varepsilon_{m}$ being the dielectric constants of the
substrate and the embedding medium (usually vacuum) respectively, and $r_{0}$
the effective screening length. Those three parameters have to be specified
for all calculations. $H_{0}$, $Y_{0}$ are Struve and Bessel functions of
second kind respectively.

Also, since the interaction decays quickly, we employ a radial cutoff, such
that for distances $r>R_{c}$ we take the interaction to be zero. Then, the
effective interaction is:

$\tilde{V}(r)=\left\\{\begin{array}[]{cc}V(a)&\text{if }r=0\\\ V(r)&\text{ if
}r<R_{c}\\\ 0&\text{else}\end{array}\right.$ (51)

where $R_{c}$ is the cutoff radius. The cutoff has two purposes: first, it
enforces the crystal symmetries in the transformed potential (as a function of
$\mathbb{k}$). Secondly, it allows to compute the summation over lattice
positions faster. Instead of evaluating the potential over all lattice
positions, we restrict the sum to the lattice positions where we know the
potential is different from zero. As for the interactions computed using the
Fourier series of the potential, we set $V(\mathbb{q}=0)=0$ to remove the long
wavelength divergence.

Lastly, it is also worth mentioning how to compute the probability of finding
the electron on a given spatial position (2.4). Since this requires two
summations over $\mathbb{k},\mathbb{k}^{\prime}$, its cost would be
$\mathcal{O}(N^{2})$. To obtain the whole wavefunction, a priori we have to
evaluate this over each position in the crystal, meaning that the cost would
be $\mathcal{O}(N^{3})$. However, this would be the worse case scenario in
which the exciton is strongly delocalized in real-space. Usually, it will
suffice to compute the real-space wavefunction on a contour of the hole
position, for a few unit cells only. To actually compute the probability, we
want to use the fact that we are storing the exciton coefficients as vectors.
First, note that (2.4) can be written as:

$\displaystyle|\psi_{X}^{\alpha\beta}(\mathbb{t}_{n}+\mathbb{R}_{e},$
$\displaystyle\mathbb{t}_{m}+\mathbb{R}_{h})|^{2}$
$\displaystyle=\left|\frac{1}{N}\sum_{v,c,\mathbb{k}}A_{vc}^{\mathbb{Q}}(\mathbb{k})e^{i\mathbb{k}\cdot(\mathbb{R_{e}}-\mathbb{R_{h}})}C^{c,\mathbb{k}+\mathbb{Q}}_{m\alpha}(C^{v,\mathbb{k}}_{n\beta})^{*}\right|^{2}$
(52)

which already reduces the complexity down to $\mathcal{O}(N)$. Then, the
probability is computed as $||A\odot C||^{2}$, where $A$ is the vector of
exciton coefficients that incorporates the exponential terms and $C$ are the
tight-binding coefficients arranged such that they match the electron-hole
pair ordering of the exciton.

##  4 Examples

So far we have discussed the theory underlying the code and its numerical
implementation. Therefore, it remains to show actual examples of the
capabilities of the code. One context where excitons are relevant is
valleytronics: materials with honeycomb structure which exhibit the band gap
at the $\mathbb{K},\mathbb{K}^{\prime}$ points of the Brillouin zone (the
”valleys”), and whose optical excitations can be tuned according to the valley
[49, 50]. The materials most commonly used for this purpose are transition
metal dichalcogenides (TMDs), with formula $\text{WX}_{2}$, where W is the
transition metal and S some chalcogenide. Another similar material that has
become highly relevant is hexagonal boron nitride (hBN), although in this case
due to its good properties as an insulating substrate [51]. These materials
have become the prototypical examples to test the capabilities of an exciton
code, and have been studied extensively. We will characterize the excitons in
both hBN and $\text{MoS}_{2}$, i.e. obtain the exciton spectrum for
$\mathbb{Q}=0$, show the associated wavefunctions and compute the optical
conductivity. We will also show how a simple strain model of hBN can be used
to break some crystal symmetries and modify the excitonic ground state. All
the calculations shown are done with the real-space approach to the
interaction matrix elements and neglecting the exchange term, unless specified
otherwise.

###  4.1 hBN

Monolayer hexagonal boron nitride has a large quasi-particle band gap, with
ab-initio calculations predicting a value of $6-8$ eV depending on the method
[52]. As we will see, the band structure of hBN is relatively flat along the
$\mathbb{M}-\mathbb{K}$ path in the Brillouin zone. This, in conjunction with
small screening results in excitons that are strongly delocalized in
reciprocal space, but are tightly bound in real space.

This material can be described easily with a minimal 2-band tight-binding
model [53], equivalent to graphene but with opposite onsite energies for each
atom of the motif. The tight-binding model for hBN reads:

$H=\sum_{i}\frac{\Delta}{2}(c^{\dagger}_{i}c_{i}-d^{\dagger}_{i}d_{i})+\sum_{\begin{subarray}{c}<i,j>\\\
i\neq j\end{subarray}}\left[tc^{\dagger}_{i}d_{j}+\text{h.c.}\right]$ (53)

where $c^{\dagger}(d^{\dagger})$ denote creation operators for B (N) atoms.
The indices $i,j$ run over unit cells, and the summation over $<i,j>$ spans
only the first-neighbours. The parameters are $t=-2.3$ eV, $\Delta/2=3.625$
eV, and the corresponding system file can be found in the code repository
under the folder /models.

Figure 1: (a) hBN lattice and (b) band structure of the tight-binding model.

From a Slater-Koster perspective, hBN is described by $\text{p}_{z}$ orbitals.
Taking the model to be spin-polarized for simplicity, there are only two bands
and there must be only one electron per unit cell (half-filling) for it to be
an insulator. Once the model is defined and the system file appropriately
constructed, we can begin setting the parameters of the calculation. First, we
need to specify the constants that appear in the Keldysh potential in Eq.
(50). These parameters determine the strength of the electrostatic interaction
and consequently affect the exciton binding energies. Here we follow previous
works to set these quantities [53], but sometimes we will be interested in
exploring the effect of tuning the dielectric constants, or instead we will
want to set them to reproduce known experimental results. Nevertheless, values
for typical substrates can be found in literature and $r_{0}$ can also be
estimated from ab-initio calculations [6, 7].

The other parameters of the exciton file are related to the convergence of the
excitons themselves. Varying the number of $\mathbb{k}$ points in the mesh,
$N_{k}$, one obtains the convergence curves shown in Fig. (2). The convergence
has been done with both the default interactions (in real-space) and with
reciprocal interactions (Fig. (2a)). For reciprocal interactions energies
converge much slower than the real-space counterpart, on top of requiring
summing over several $\mathbb{G}$ reciprocal cells. In materials with highly
localized excitons in $\mathbb{k}$ space, it usually suffices to take only
$\mathbb{G}=0$ (e.g. $\text{MoS}_{2}$). However, we will see later that hBN
excitons are highly delocalized in reciprocal space, which is why the
interaction can see neighbouring reciprocal unit cells.

Figure 2: (a) Convergence of the ground state and first and second excited
states as a function of the number of $\mathbb{k}$ points for hBN, computed
both with interactions in real and reciprocal space. The reciprocal space
calculations have to be converged also with respect to the number of
reciprocal lattice vectors included, $N_{G}$, with $N_{G}=25$ in this case.
(b) Measured calculation time as a function of $N$. For both the real and
reciprocal space calculations, the asymptotic behaviour is
$\mathcal{O}(N^{3})$. However, at small values of $N$ the required time is
partially dominated by the BSE matrix initialization, which in both cases
scales as $\mathcal{O}(N^{2})$.

After checking convergence, we can start studying the exciton themselves. The
energies of the first 8 states and their degeneracies are given in table (1).
To make sense of the degeneracies, one has to check the character table of the
point group of the material: hBN has the crystallographic point group
$D_{3h}$, with both one- and two-dimensional irreducible representations.
Since the symmetry operations and their action on single-particle states are
specific to each problem, the code does not address the problem of identifying
the irreducible representation of each exciton, nor labeling them in terms of
symmetry eigenvalues. Instead, we only check that the $\mathbb{Q}-$excitonic
wavefunctions have the allowed degeneracies and (40) is invariant under the
little group at $\mathbb{Q}$.

$n$ | Energy (eV) | Binding energy (eV) | Degeneracy
---|---|---|---
1 | 5.3357 | -1.9143 | 2
2 | 6.0738 | -1.1762 | 1
3 | 6.1641 | -1.0859 | 2
4 | 6.1723 | -1.0777 | 1
5 | 6.3511 | -0.8989 | 2
Table 1: Exciton spectrum from the tight-binding model for hBN computed with
$N_{k}=60^{2}$. The binding energy $E_{b}$ is defined as $E_{b}=E_{X}-\Delta$,
where $\Delta$ is the gap of the system.

The $\mathbb{k}$ probability densities of the first eight excitonic states,
grouped by degenerate levels, are shown in Fig. (3). Each energy level has the
symmetry of the lattice, as expected since we are plotting (40). The
additional symmetry in this case is due to time-reversal symmetry and the fact
that $\mathbb{Q}=0$ is a time-reversal invariant momenta. We see that the
wavefunctions peak at the valleys, although they also spread over the
$\mathbb{K}-\mathbb{M}-\mathbb{K}^{\prime}$ paths. This means that the
excitons are formed by strongly interacting electron-hole pairs in
$\mathbb{k}$ space, which explains why we need to sum over several reciprocal
cells when using the reciprocal interactions. As for the shape of excitons, we
find the common pattern: the first state peaks is $s$-like in the sense that
it does not have nodes. The next state would be $p$-like and so on. Note that
the hydrogen analogy only concerns the shape of the wavefunctions, and not the
energy spectrum, which in general differs from the hydrogen series.

Figure 3: Plot of the $\mathbb{k}$ exciton probability densities in TB hBN,
for the first 5 energy levels (from left to right, top to bottom). For each
level, we actually show
$|\Psi(\mathbb{k})|^{2}=\sum_{n}|\psi_{n}(\mathbb{k})|^{2}$ for
$\mathbb{Q}=0$, where the index $n$ runs over degenerate states.

Since the excitons are delocalized in reciprocal space, we expect them to be
strongly localized in real space. The real space densities of each degenerate
level are shown in Fig. (4). The hydrogenic picture makes more sense when
looking at the real-space wavefunction, since it can be understood then as the
problem of two interacting opposite sign charges. The spectrum and the
degeneracies do not match that of hydrogen, but the wavefunctions behave
radially as we would expect.

Figure 4: Plot of the real-space exciton probability densities in TB hBN with
$\mathbb{Q}=0$, for the first 5 energy levels (from left to right, top to
bottom). The red dot shows the position of the hole. For each level we plot
the sum of the probability densities of each state of the degenerate subspace,
$|\Psi(r_{e},r_{h})|^{2}=\sum_{n}|\psi_{n}(r_{e},r_{h})|^{2}$.

In hBN the spin-orbit coupling is small and it suffices to compute the
excitons as a spinless system, in particular given that we are also neglecting
the exchange interaction. If we consider a spinful system, without exchange
again, we obtain exactly the same energy levels but now four-fold degenerate
(on top of the previous spatial degeneracy). The same stands for both types of
wavefunctions.

Our study of the exciton spectrum in hBN concludes with the calculation of the
optical conductivity [16, 39, 52], which reflects the light absorbance from a
source up to a constant factor. So far we have not discussed which excitons of
the spectrum are bright or dark. This can be seen through the calculation of
the optical oscillator strengths within Eq. (41), which determine the
transition rate for photon emission. The frequency-dependent conductivity of
monolayer hBN is given in Fig. (5). Electron-hole interactions move the
spectral power from the continuum to pronounced sub-band gap peaks. Attending
to Table (1) and Fig. (3), we see that non-degenerate excitons with mainly $s$
character are bright. The relative height of the peaks can be understood by
looking at the magnitude of the wavefunctions near the $K$ and $K$’ points.
All bright excitons can be excited with linearly polarized light along two
orthonormal polarization directions, giving rise to an isotropic conductivity
consistent with the $D_{3h}$ point group of the material.

Figure 5: Optical conductivity of monolayer hBN as a function of the incident
energy. We compare the conductivity obtained with the Kubo formula in the
independent particle approximation (IPA), and with the BSE, which shows a
dramatic change from the inclusion of excitons.

It is of interest to check the validity of the results against a more refined
description of the band structure of the material. This can be done with the
code by using a local orbital-based DFT calculation as the starting
Hamiltonian, instead of using a parametrized tight-binding model. The exciton
energies will depend on the gap as estimated from the functional used, but we
expect to get similar wavefunctions and conductivity. Since we consider
several orbitals for each chemical species now, we have multiple valence and
conduction bands so we should converge the exciton with respect to the number
of bands as well. It is a proper check to do, but in this case the different
bands are well separated, so their effect should be negligible.

The DFT band structure and the wavefunctions of the ground state exciton are
shown in Fig. (6). One could use standard LDA functionals, but here we opt for
a hybrid functional (HSE06[54] in this case), which is efficiently implemented
in CRYSTAL [37]. This type of functional yields a better estimation of the
single-particle gap due to a different treatment of the exchange-correlation
term. For both LDA (not shown) and hybrid functionals such as HSE06, the
wavefunctions closely resemble those obtained with TB models. For instance, we
observe the same sublattice polarization present in the TB real-space
densities with the HSE06 calculation (Fig. 6c). The energy spectrum shows the
same degeneracies, although the positions of some of the levels are exchanged.

Figure 6: (a) DFT band structure of monolayer hBN as obtained with the HSE06
functional. (b) Reciprocal and (c) real-space probability densities of the
$\mathbb{Q}=0$ ground state exciton with $N_{k}=60^{2}$, $N_{c}=N_{v}=1$. The
DFT calculation involved a basis size of 36. We have run successfully exciton
calculations on different systems with varying basis sizes, from 8 to 92.

To illustrate the applicability of the code beyond standard cases, we now
study the effect of strain on the exciton spectrum. If we apply some uniaxial
in-plane strain along the $x$ axis, the point group of the material will
change to $C_{2v}$ (with rotation axis along $x$). The degeneracy of the
ground state came from the spatial symmetries, meaning that it should be
broken for any strain value, given that all irreducible representations of
$C_{2v}$ are of dimension 1. Therefore, we can study the energy splitting of
the ground state as a function of the applied strain.

The strain model used is fairly straightforward. Based on the original tight-
binding model, we now consider the hopping parameters to have an exponential
dependence on the distance:

$t(r)=t_{0}e^{-a(r-r_{0})},$ (54)

where $a$ is some decay length, $t_{0}$ the original value of the hopping and
$r_{0}$ the reference length. Additionally, the distortion of the lattice due
to strain is taken to affect only bonds parallel to the strain. A rigorous
approach would have to implement appropriate distortion of all atomic
positions according to the stress tensor [55], but for our purposes this
simple model suffices. This is illustrated in Fig. (7)

Figure 7: Schematic of the distortion of the hBN lattice due to the
application of uniaxial strain along the $x$ axis. Note that in reality all
bonds should be distorted, due to the phenomena of Poisson contraction.

The procedure to study the exciton spectrum as a function of strain is as
follows: we generate different system files (i.e. different Hamiltonians) for
different values of the strain, which translates into different atomic
positions. Then, we run the exciton simulation for each system file, storing
the energies. As we expected, now all states are non-degenerate because of the
symmetry group $C_{2v}$. We can plot the ground state splitting as a function
of strain, which is shown in Fig. 8(a). In Fig.8(b) we show the conductivity
for some finite value of the strain. The response is no longer isotropic due
to the lattice symmetry breaking caused by strain, where the exciton peaks
shift for both light polarizations.

Figure 8: (a) Energy splitting of the ground state exciton as a function of
strain in hBN. We observe that the splitting is linear on strain, for small
values. (b) Frequency-dependend conductivity on strained hBN,
$\varepsilon=0.1$. The dashed line shows the position of the ground state for
$\varepsilon=0$.

###  4.2 $\text{MoS}_{2}$

To conclude the examples section, we also analyze the exciton spectrum of
$\text{MoS}_{2}$. Same as hBN, in monolayer form this material presents itself
in a honeycomb lattice, although it is not planar. Instead, it is formed by
three layers of composition S-Mo-S respectively. The description of the band
structure of $\text{MoS}_{2}$ requires a more complex model, which is why we
use it to showcase the code. We use a Slater-Koster tight-binding model [56],
where each chemical species has a different set of orbitals (Mo has $d$
orbitals, and S only $p$ orbitals). This, together with the non-negligible
spin-orbit coupling results in a more complex band structure than that of hBN.
Both the lattice and the band structure can be found in Fig. (9).

Figure 9: (a) Crystal and (b) tight-binding band structure of $\text{MoS}_{2}$. $n$ | Energy (eV) | Binding energy (eV) | Degeneracy
---|---|---|---
1 | 1.7673 | -0.3527 | 2
2 | 1.7797 | -0.3403 | 2
3 | 1.9105 | -0.2095 | 2
4 | 1.9232 | -0.1968 | 2
Table 2: Exciton spectrum from the tight-binding model for MoS2 computed with
$N_{k}=40^{2}$, $N_{v}=N_{c}=2$. This model has a direct gap at $\mathbb{K}$
of 2.12 eV, used to compute the shown exciton binding energies.

After checking convergence with the number of $\mathbb{k}$ points and the
number of bands, we obtain the spectrum shown in table (2). In this case, the
point group of the material is again $D_{3h}$ and the irreducible
representations realized by the wavefunctions at $\mathbb{Q}=0$ are compatible
with the character table of the group [57]. As before, to ensure that the
excitons were computed correctly we can plot the total densities to ensure
that they have the expected symmetries.

In Fig. (10a) we show the reciprocal probability density of the first energy
level. As opposed to hBN, we observe that the states are strongly localized at
the valleys. Resolving the degeneracy by labeling each exciton with the
$C_{3}$ eigenvalues would result in each exciton localized in a different
valley [36]. This shows that at least the low energy spectrum of
$\text{MoS}_{2}$ can be studied at one valley instead of the whole BZ [58].
This allows to get a more precise description of the exciton since one can use
a more refined mesh. The wavefunction for the exciton obtained at one valley
can be seen in Fig. (10b), using the code feature to reduce the BZ mesh by
some integer factor. For higher excited states this does not hold since the
states become more extended across the BZ, reaching both valleys.

Figure 10: (a) Probability density of the ground state exciton in MoS2
obtained over the full BZ with $N_{k}=60^{2}$ for $\mathbb{Q}=0$. (b) Ground
state exciton computed in a contour of the $\mathbb{K}$ valley with
$N_{k}=30^{2}$ with a reduction factor of 2. Both calculations were done with
$N_{c}=N_{v}=2$.

Since the excitons are very localized in reciprocal space, they should be
delocalized in real-space, meaning that the radius of the exciton should be
large (e.g. compared to that of hBN). To complete the characterization of the
excitons, we calculate the optical conductivity as shown in Fig. (11) While
the exciton energies converge quickly with $N_{k}$, it is usually necessary to
include more $\mathbb{k}$ points in the calculation of the optical
conductivity in order to smooth unphysical oscillations derived from the
discrete mesh. As it can be seen, the shape of the spectrum matches previous
tight-binding studies [59, 36] and agrees well with ab-initio results [60]. At
low energies, the optical conductivity of MoS2 presents the characteristic A
and B exciton peaks, that are understood considering the main spin-allowed
electron-hole excitations at the $\mathbb{K}$ and $\mathbb{K}^{\prime}$
points. The split of $\sim$$100$ meV between such peaks reflects the effect of
SOC in TMD materials [61]. At higher energies, the main feature of the spectra
is a pronounced peak similar to the non-interacting case but red-shifted in
energy. The excitons giving rise to such peak are often called “C” excitons
and were fully characterized in Ref. [59], already showing the potential of
tight-binding methods for studying new exciton physics.

Figure 11: Optical conductivity of MoS2 with and without excitons. The BSE
calculation was done with $N_{k}=34^{2}$, $N_{v}=2$ and $N_{c}=6$. The first
two peaks correspond to the A B excitons at the valleys, while the rest of the
conductivity can be regarded as a shift of the non-interacting one.

##  5 Conclusions

We have developed a software package that allows to solve the Bethe-Salpeter
equation constructed from either tight-binding models or DFT calculations
based on localized orbitals. By considering orbitals as point-like, the
computation of the interactions becomes drastically simplified. Together with
an effective screening, this results in a fast determination of the BSE
matrix. More specifically, our real-space implementation of the interaction
matrix elements is shown to be faster and more precise than its reciprocal-
space counterpart, which is the formulation more commonly used.

As in GW-BSE approximations, the starting band structure plays a crucial role
in determining the resulting exciton spectrum. Therefore it is key to select
the best possible functional (typically hybrids) or the most accurate tight-
binding models that capture the most prominent features of the band structure.
Then, by choosing appropriately the screening parameters, it is possible to
reproduce the results of GW-BSE or similar first-principles codes at a
fraction of the computational cost.

The Xatu code currently provides all the tools needed to extract and
characterize the exciton spectrum, either using the binary or via its API.
Nevertheless, the package is still under development, as new functionalities
and optimizations are added. Our future plans include giving support for
distributed parallelism to enable bigger system sizes and calculation of
different excitation types such as trions or biexcitons. The code is currently
aimed at the description of 2D materials, but it can support 0D and 3D
systems. Since the Keldysh potential is only adequate for 2D systems, we will
implement additional potentials suitable for different dimensionalities. We
also plan to add the possibility of performing exact calculations of the
interaction matrix elements when using Gaussian-based DFT codes to compute the
band structure. Currently we provide an interface with the CRYSTAL code [37],
and ideally more interfaces to community codes will be added over time, such
as SIESTA [62] or Wannier90 [63]. The project has been released under an open-
source license and as such community contributions are welcome and encouraged.

Note added: Upon completion of this work we became aware of a very recent
submission which also addresses the problem of determining the exciton
spectrum from Wannier-based tight-binding models [64]. Nevertheless, we
believe that the thorough characterization of excitons we provide here can be
advantageous and complementary to other different tools.

##  6 Acknowledgments

The authors acknowledge financial support from Spanish MICINN (Grant Nos.
PID2019-109539GB-C43 & TED2021-131323B-I00), María de Maeztu Program for Units
of Excellence in R&D (GrantNo.CEX2018-000805-M), Comunidad Autónoma de Madrid
through the Nanomag COST-CM Program (GrantNo.S2018/NMT-4321), Generalitat
Valenciana through Programa Prometeo (2021/017), Centro de Computación
Científica of the Universidad Autónoma de Madrid, and Red Española de
Supercomputación.

##  Appendix A Symmetry

Throughout the document, we have mentioned and used several times the fact
that each individual exciton state does not necessarily have the symmetries of
the lattice, but it is the sum of the square amplitude of them who does. In
this appendix we will give a proof of this statement: Given
$|\Psi|^{2}=\sum_{n}|\psi_{n}|^{2}$, where $\psi_{n}$ denotes the wavefunction
of the exciton states on some degenerate subspace of $PHP$, and given some
symmetry operation $C$ such that $[H,C]=0$, then

$C|\Psi|^{2}=|\Psi|^{2}.$ (55)

First we have the consider the action of the symmetry operator $C$ on an
exciton state. Given that the eigenstates of a degenerate subspace of $H$ are
not in general eigenstates of $C$, the most general action is to mix the
degenerate states, i.e.:

$C\psi_{n}=\sum_{i}\alpha_{in}\psi_{i}$ (56)

The coefficients $\alpha_{in}$ are the matrix elements of $C$. To prove (55),
we need to know the action of $C$ on the squared amplitude, $C|\psi_{n}|^{2}$.
So first we want to prove the following property:

$C|\psi_{n}|^{2}=|C\psi_{n}|^{2}$ (57)

This can be proven using the action of the symmetry operation on the
coordinate of the wavefunction, i.e. $C\psi_{n}(x)=\psi_{n}(C^{-1}x)$:

$\displaystyle C|\psi_{n}|^{2}(x)$
$\displaystyle=|\psi_{n}|^{2}(C^{-1}x)=\psi_{n}(C^{-1}x)\psi_{n}^{*}(C^{-1}x)$
$\displaystyle=C\psi_{n}(x)C\psi_{n}^{*}(x)=|C\psi_{n}|^{2}(x)$ (58)

where we have also used that $C\psi_{n}^{*}(x)=(C\psi_{n})^{*}(x)$. This last
identity can be proved conjugating the action of the symmetry on the
coordinates:

$(C\psi_{n})^{*}(x)=\psi_{n}^{*}(C^{-1}x)=C\psi_{n}^{*}(x)$ (59)

This enables us to compute $C|\psi_{n}|^{2}$ in terms of a expansion on the
states of the degenerate subspace:

$C|\psi_{n}|^{2}=|C\psi_{n}|^{2}=\left|\sum_{i}\alpha_{in}\psi_{i}\right|^{2}=\sum_{i,j}\alpha_{in}\alpha^{*}_{jn}\psi_{i}\psi^{*}_{j}$
(60)

Finally, with this expression we can prove the symmetry invariance of
$|\Psi|^{2}=\sum_{n}|\psi_{n}|^{2}$. To do so, we act with the symmetry
operation $C$ on $|\Psi|^{2}$:

$\displaystyle C|\Psi^{2}|$
$\displaystyle=\sum_{n}C|\psi_{n}|^{2}=\sum_{n}\left[\sum_{i,j}\alpha_{in}\alpha^{*}_{jn}\psi_{i}\psi^{*}_{j}\right]$
$\displaystyle=\sum_{ij}\left[\sum_{n}\alpha_{in}\alpha^{*}_{jn}\right]\psi_{i}\psi^{*}_{j}=\sum_{i}|\psi_{i}|^{2}=|\Psi|^{2}$
(61)

where we have used that $C$ is unitary, i.e.
$\sum_{n}\alpha_{in}\alpha^{*}_{jn}=\delta_{ij}$. This proves that the sum of
the squared amplitude of each degenerate state is invariant under the symmetry
operations.

On a different note, back in the examples we used the character table of the
point group of the solid to justify the observed state degeneracies, and in
the previous proof we also considered some general symmetry $C$ such that
$[H,C]=0$. For the abstract, unrepresented Hamiltonian, given any operation
$C$ of the point group of the solid, it is true that $[H,C]=0$. However, we
are not working with the total Hamiltonian, but with a sector of it. So one
must actually look for symmetry operations that commute with $PHP$:

$[PHP,C]=0$ (62)

Since the sectors of electron-hole pairs of different momentum are
disconnected, we can define
$\tilde{H}(\mathbb{Q})=P_{\mathbb{Q}}HP_{\mathbb{\mathbb{Q}}}$, where
$P_{\mathbb{Q}}$ is the projector over electron-hole pairs of $\mathbb{Q}$
total momentum. This Hamiltonian is analogous to the Bloch Hamiltonian
$H(\mathbb{k})$, and it can be shown that it transforms in the same way:

$C^{-1}\tilde{H}(\mathbb{Q})C=\tilde{H}(C^{-1}\mathbb{Q})$ (63)

meaning that for $\mathbb{Q}=0$ the symmetry group is the crystallographic
point group, but for $\mathbb{Q}\neq 0$ the Hamiltonian is invariant only
under symmetry operations of the little group of $\mathbb{Q}$, whose
irreducible representations thus dictate the (unitary) transformation
properties of the $\mathbb{Q}-$excitonic wavefunctions.

Proof that $H(\mathbb{Q})$ transforms as the Bloch Hamiltonian $H(\mathbb{k})$
under symmetry operations:

$\displaystyle C^{-1}H(\mathbb{Q})C$
$\displaystyle=C^{-1}P_{\mathbb{Q}}CC^{-1}HCC^{-1}P_{\mathbb{Q}}C$
$\displaystyle=C^{-1}P_{\mathbb{Q}}CHC^{-1}P_{\mathbb{Q}}C$ (64)

where we have used that $[H,C]=0$. So we only need to see how the projectors
transform under the symmetry operation to determine how $H(\mathbb{Q})$
transforms.

$\displaystyle
C^{-1}P_{\mathbb{Q}}C=\sum_{\mathbb{k},v,c}C^{-1}c^{\dagger}_{c,\mathbb{k+Q}}c_{v\mathbb{k}}\ket{GS}\bra{GS}c^{\dagger}_{v\mathbb{k}}c_{c,\mathbb{k+Q}}C$
(65)

Inserting identities, we can transform each creation/annihilation operator
according to
$C^{-1}c^{\dagger}_{n\mathbb{k}}C=c^{\dagger}_{n,C^{-1}\mathbb{k}}$, up to an
arbitrary phase that is cancelled in (65). From (3) it follows that the Fermi
sea is invariant under point group operations, i.e. $C\ket{GS}=\ket{GS}$
(again, up to an arbitrary phase that is cancelled), we arrive to the
following expression:

$\displaystyle C^{-1}P_{\mathbb{Q}}C$
$\displaystyle=\sum_{\mathbb{k},v,c}c^{\dagger}_{c,C^{-1}\mathbb{k}+C^{-1}\mathbb{Q}}c_{v,C^{-1}\mathbb{k}}\ket{GS}\bra{GS}c^{\dagger}_{v,C^{-1}\mathbb{k}}c_{c,C^{-1}\mathbb{k}+C^{-1}\mathbb{Q}}$
(66)

From the C-invariance of the BZ in $\mathbb{k}-$space, we arrive at the final
expression for the transformed projector:

$\displaystyle
C^{-1}P_{\mathbb{Q}}C=\sum_{\mathbb{k},v,c}\ket{v,c,\mathbb{k},C^{-1}\mathbb{Q}}\bra{v,c,\mathbb{k},C^{-1}\mathbb{Q}}=P_{C^{-1}\mathbb{Q}}$
(67)

Therefore, the projected exciton Hamiltonian $H(\mathbb{Q})$ also transforms
in the same way:

$\displaystyle C^{-1}H(\mathbb{Q})C=H(C^{-1}\mathbb{Q})$ (68)

Likewise, the application of time-reversal yields
$T^{-1}H(\mathbb{Q})T=H(-\mathbb{Q})$ whenever it is a symmetry of the system.

##  Appendix B Usage

The installation instructions can be found at the repository
https://github.com/alejandrojuria/xatu, so it will not be discussed here. The
code has been developed with a hybrid approach in mind: one can resort to
configuration files to run the program, in analogy with DFT codes, or instead
program both the non-interacting system and run the exciton simulation using
the provided API. First we are going to discuss its usage with configuration
files. The basic usage as a CLI program is described by:

⬇

xatu [OPTIONS] systemfile [excitonfile]

The executable always expects one file describing the system where we want to
compute the excitons, and then another file specifying the parameters of the
simulation. Their content is addressed in the next sections. The executable
can also take optional flags, generally to tune the output of the simulation.
By default, running the program without additional flags prints the exciton
energies, without writing the results to any file.

⬇

-h (--help)

Used to print a help message with the usage of the executable and a list of
all possible flags that may be passed. The simulation is not performed (even
in presence of configuration files).

⬇

-s (--states) nstates

The number of states specified with this flag is also used for any of the
output flags. By default, the number of states is 8 (i.e. if the flag is not
present).

⬇

-p (--precision) decimals

One can specify the number of decimals used when printing the exciton
energies. This is relevant to detect state degeneracy without inspecting
manually the states. Defaults to 6 decimals if not present.

⬇

-d (--dft) [ncells]

This flag is used to indicate that the systemfile provided corresponds to a
CRYSTAL output file, instead of following the standarized format. DFT
calculations usually involve several unit cells to determine the Bloch
Hamiltonian, so the optional value ncells can be passed to specify how many we
want to take into account. Otherwise all of them are read and used.

⬇

-eck (--energy, \--eigenstates, \--kwf)

The optional flags -e, -c, -k, -r are used to specify which exciton output is
written to file. -e writes the energies, -c writes the eigenvectors, -k writes
the reciprocal density. Note that they can be combined instead of being
written separately (e.g. -ek instead of -e -k).

⬇

-r (--rswf) [holeIndex] [-r ncells]

Used to write the real-space probability densities to a file. One can give the
index of the atom where the hole is located (defaults to first atom of the
motif). It can be used a second time to specify the number of unit cells where
we want to compute the amplitude (e.g. -r 1 -r 10 fixes the hole at the second
atom of the motif, and uses 10 unit cells along each axis).

⬇

-s (--spin)

Computes the total spin of the excitons, and writes it to a file. This assumes
that the single-particle basis includes spin without performing any check, so
incorrect usage could result in wrong results or runtime errors.

⬇

-a (--absorption)

Computes the optical conductivity (which reflects the absorption of light up
to a constant factor) as a function of frequency using the exciton spectrum,
and saves the result to a file. A file named ”kubo_w.in” with the adequate
format (shown below) must be present in the working directory.

⬇

-m (--method) diag | davidson | sparse

Choose method to obtain the eigenstates of the BSE. By default, full
diagonalization is used. If the Davidson or sparse (Lanczos) method is
selected, then it is used to compute the number of states specified before.

⬇

-b (--bands) kpointsfile

To check that the system file was written correctly, one can use this option
to diagonalize the Bloch Hamiltonian on the $\mathbb{k}$ points specified on a
file, and write the energy bands to a file. No exciton calculation is
performed.

###  B.1 Structure of a system file

The system configuration files contain all the information needed to
characterize completely the material of study: it provides the lattice vectors
and motif positions, which is required for the real-space evaluation of the
excitons. Then, we have the number of orbitals of each unique chemical
species, which is needed to compute the matrix elements correctly, and the
filling, which determines which bands participate in the formation of the
exciton. Finally, the file contains the matrices needed to build the Bloch
Hamiltonian, this is, the Fock matrices $H(\mathbb{R})$ and their
corresponding Bravais vectors $\mathbb{R}$. The Bloch Hamiltonian is then
reconstructed as:

$H(\mathbb{k})=\sum_{\mathbb{R}}H(\mathbb{R})e^{i\mathbb{k}\cdot\mathbb{R}}$
(69)

Note that even though one has to provide the orbitals of each species, the
specific type of orbital is not needed since the interaction is computed using
the point-like approximation. A system file is specified using labels for each
block. Blocks always begin with the block delimiter #, followed by a label. A
block is then defined as all the content between two consecutive block
delimiters. The expected content for each label will be discussed next. Any
line containing ! is regarded as a comment, and empty lines are skipped.

  * #

BravaisLattice: Basis vectors of the Bravais lattice. The number of vectors
present is also used to determine the dimensionality of the system. The
expected format is one vector per line, x y z.

  * #

Motif: List with the positions and chemical species of all atoms of the motif
(unit cell). The chemical species are specified with an integer index, used
later to retrieve the number of orbitals of that species. The expected format
is one atom per line, x y z index.

  * #

Orbitals: Number of orbitals of each chemical species present. The position of
the number of orbitals for each species follows the indexing used in the motif
block. This block expects one or more numbers of orbitals, the same as the
number of different species present, n1 [n2 ...].

  * #

Filling: Total number of electrons in the unit cell. Required to identify the
Fermi level, which is the reference point in the construction of the excitons.
Must be an integer number.

  * #

BravaisVectors: List of Bravais vectors $\mathbb{R}$ that participate in the
construction of the Bloch Hamiltonian (69). Expected one per line, in format x
y z.

  * #

FockMatrices: Matrices $H(\mathbb{R})$ that construct the Bloch Hamiltonian
$H(\mathbb{k})$. The matrices must be fully defined, i.e., they cannot be
triangular, since the code does not use hermiticity to generate the Bloch
Hamiltonian. The Fock matrices given must follow the ordering given in the
block BravaisVectors. The matrices can be real or complex, and each one must
be separated from the next using the delimiter &. In case the matrices are
complex, the real and imaginary parts must be separated by a space, and the
complex part must carry the imaginary number symbol (e.g. $1.5$ $-2.1j$). Both
$i$ and $j$ can be used.

  * #

[OverlapMatrices]: In case that the orbitals used are not orthonormal, one can
optionally provide the overlap matrices $S(\mathbb{R})$. The overlap in
$\mathbb{k}$ space is given by:

$S(\mathbb{k})=\sum_{\mathbb{R}}S(\mathbb{R})e^{i\mathbb{k}\cdot\mathbb{R}}$

This is necessary to be able to reproduce the bands, which come from solving
the generalized eigenvalue problem
$H(\mathbb{k})S(\mathbb{k})\Psi=ES(\mathbb{k})\Psi$. This will be specially
necessary if the system was determined using DFT, since in tight-binding we
usually assume orthonormality. This block follows the same rules as
FockMatrices: each matrix $S(\mathbb{R})$ must be separated with the delimiter
&, and they must follow the order given in BravaisVectors.

Several examples of valid system files are provided in the code repository,
under the folder /models.

###  B.2 Structure of an exciton file

The purpose of the system file was to specify completely the system where we
want to compute the excitons. Then, the exciton file is used to describe the
excitons themselves: number of points in the mesh or submesh, bands that
participate and center-of-mass momentum for example, as well as some
additional flags. The idea is to keep the functionality as orthogonal as
possible between files. With one system file, we can test for the convergence
of the excitons with the number of kpoints, or with the number of bands
modifying the exciton file only. Finally, we have the runtime options of the
program, which in general do not affect the energy and modify the output
exclusively. The philosophy is to maximize the reproducibility and facilitate
tracking of the experiments.

The exciton files are built following the rules of the system files. They are
composed of blocks, starting with #. Each block has a label, which determines
the expected content of the block. Next we provide a list of the possible
parameters used in the construction of an exciton file:

  * #

Label: Used to specify the name of the files containing the output of the
program. The files will be named [Label].eigval, [Label].kwf, etc.

  * #

Bands: Number of bands above and below the Fermi level. The minimum value is
1, to describe one conduction band and one valence band (i.e. only one
combination of bands).

  * #

[BandList]: As an alternative to Bands, one can specify a list with the
indices of the bands that compose the exciton. 0 is taken as the last valence
band, meaning that 1 would be the first conduction band, -1 is the second
valence band and so on. This option can be used to generate asymmetric
combinations of bands. It overrides the Bands block.

  * #

Ncells: Number of points in one direction of the Brillouin zone, or
equivalently number of unit cells along one axis. The same number of points is
taken along all directions.

  * #

[Submesh]: Used to specify a submesh of the Brillouin zone. Takes a positive
integer $m$, which divides the BZ along each axis by that factor. The
resulting area is meshed with the number of points specified in the Ncells
block. This option can become memory intensive (it scales as
$\mathcal{O}(m^{d})$, $d$ the dimension).

  * #

[ShiftMesh]: In case that we are using a submesh, then probably we also want
to shift the meshed area to center it at the gap, where the exciton peaks.
Takes a vector with its components, kx ky kz.

  * #

Dielectric: The Keldysh interaction requires setting the dielectric constants
of substrate $\epsilon_{s}$, the medium $\epsilon_{m}$ and the screening
length $r_{0}$, which involves the dielectric constant of the material. This
block expects three values, es em r0.

  * #

[TotalMomentum]: One can optionally specify the total or center-of-mass
momentum $\mathbb{Q}$ of the exciton. By default, it is taken to be zero,
unless this block is specified. It expects a vector in form qx qy qz.

  * #

[Reciprocal]: If present, the interaction matrix elements are computed in
reciprocal space instead of direct space, which is the default. It takes an
integer argument to specify the number of reciprocal cells to sum over, nG.

  * #

[Exchange]: Flag to turn on the exchange interaction. By default computations
neglect the exchange, and use only the direct term. It has to be set to true
or false.

  * #

[Scissor]: Used to specify a scissor shift of the bands to correct the gap.
This optional field takes a single value, shift

As it can be seen, a minimum exciton simulation only requires specifying the
number of bands, the number of $\mathbb{k}$ points and the dielectric
constants. The modification of any of the present parameters is expected to
result in a variation of the exciton results (energies, wavefunctions,
conductivity), which is why all this parameters have been delegated to the
same file.

###  B.3 Absorption file

For the calculation of the Kubo conductivity, one needs to provide a separate
input file named kubo_w.in in the working folder. This file is used to specify
all parameters relative to the conductivity calculation, namely the desired
energy interval, the point sampling and the broadening to be used, as well as
the output files. Its format is as follows:

    #initial frequency (eV)
    5
    #frequency range (eV)
    8
    #number of frequency points (integer)
    300
    #broadening parameter (eV)
    0.05
    #type of broadening
    lorentzian
    #output kubo name files
    kubo_hBN_sp.dat
    kubo_hBN_ex.dat

Do note that as opposed to the previous configuration files, the name of each
section starting with $\\#$ is not actually relevant for the parsing; the
program always expects the same fields to be present in the file, and in the
same order as presented here. For the broadening, three different options are
allowed (’lorentzian’, ’exponential’, ’gaussian’)

###  B.4 As a library

So far we have discussed a more streamlined usage of the package. In some
cases, however, the user could benefit from accessing directly the results of
an exciton calculation, instead of having to dump it to a file to postprocess
it later. To enable this possibility, the package has been also designed as a
library, meaning that one can import the classes and functions defined in the
API, and use them to build some extra functionality. Some use cases could be
scenarios with exciton interactions, such as exciton-exciton interactions or
exciton-polaritons.

To do so, the package provides a header file which defines a namespace. Within
the namespace we have access to all the exciton functionality, which is
completely documented. For instructions on how to build the documentation, we
refer to the project repository where the most up-to-date information will be
present. Additionally, some usage examples can be found under the root
directory in the folder /main.

The outline for a general exciton simulation is the following: one first has
to create a System object, which can be done with a system file.
Alternatively, one can define a subclass that inherits from System, and use it
to implement the desired behaviour (namely the Bloch Hamiltonian). Then this
System is passed on to the Exciton class, which we configure with the desired
parameters. The interacting Hamiltonian is initialized and solved, returning a
Result object which contains the eigenvalues and eigenvectors. With this, now
we can compute some observables, or instead use these states to perform some
other calculations out of the scope of the code.

###  B.5 Output

To conclude the usage section, we will describe briefly the structure of the
output files. Here we describe how each file is written so the user can write
their own custom routines; we also provide some example Python scripts under
the folder /plot.

  * 1.

Energy: The energy file has in the first line the total number of energies
written in the file. The second line contains all the energies, separated by a
tabulation, e1 e2 ... en. All energies are written, including degenerate
levels. All the exciton energies are given with respect to the Fermi sea
energy. To obtain the binding energy, one must substract the gap from the
exciton energy. Units are $[E]=$ eV.

  * 2.

States: The first line contains the dimension of the BSE matrix $n$, i.e. the
number of different electron-hole pairs. The next $n$ lines specify the
valence, conduction bands of each electron-hole pair and their $\mathbb{k}$
point, kx ky kz v c. Afterwards, each line specifies completely the
coefficients of each exciton state. The format per line is: Re(A1) Im(A1)
Re(A2) Im(A2) ....

  * 3.

Reciprocal probability density: For the reciprocal density, on each line we
specify the coordinates of the k point and the associated probability, kx ky
kz P. Each state is separated from the next by a delimiter #. Units are
$[k]=$Å-1.

  * 4.

Real-space probability density: The first line has the coordinates of the
hole, hx hy hz. The following ones each have the coordinates of one atomic
position, and the probability of finding the electron: x y z P. Densities for
different states are separated by #. Units are $[x]=$Å.

  * 5.

Spin: On each line we write the index of the current exciton, and next the
total spin projection, the hole and the electron spin, n St Sh Se. Spin units
are $[S_{z}]=\hbar$.

  * 6.

Absorption: Both conductivities with and without exciton effects are computed
and written to two different files. Each row contains the following columns:
$\omega,\ \sigma_{xx},\ \sigma_{xy},\ \sigma_{xz}$, $\sigma_{yx},\
\sigma_{yy},\ \sigma_{yz},\ \sigma_{zx},\ \sigma_{zy},\ \sigma_{zz}$. Units
are $[\omega]=$ eV, $[\sigma_{ij}]=e^{2}/\hbar$.

## References

  * [1] W. Hanke, L. J. Sham, Many-particle effects in the optical spectrum of a semiconductor, Phys. Rev. B 21 (1980) 4656–4673.  doi:10.1103/PhysRevB.21.4656.
URL https://link.aps.org/doi/10.1103/PhysRevB.21.4656

  * [2] A. Chernikov, T. C. Berkelbach, H. M. Hill, A. Rigosi, Y. Li, B. Aslan, D. R. Reichman, M. S. Hybertsen, T. F. Heinz, Exciton binding energy and nonhydrogenic rydberg series in monolayer ${\mathrm{ws}}_{2}$, Phys. Rev. Lett. 113 (2014) 076802.  doi:10.1103/PhysRevLett.113.076802.
URL https://link.aps.org/doi/10.1103/PhysRevLett.113.076802

  * [3] P. Cudazzo, I. V. Tokatly, A. Rubio, Dielectric screening in two-dimensional insulators: Implications for excitonic and impurity states in graphane, Phys. Rev. B 84 (2011) 085406.  doi:10.1103/PhysRevB.84.085406.
URL https://link.aps.org/doi/10.1103/PhysRevB.84.085406

  * [4] T. Mueller, E. Malic, Exciton physics and device application of two-dimensional transition metal dichalcogenide semiconductors, npj 2D Materials and Applications 2 (1) (2018) 29.  doi:10.1038/s41699-018-0074-2.
URL https://doi.org/10.1038/s41699-018-0074-2

  * [5] M. F. C. M. Quintela, J. C. G. Henriques, L. G. M. Tenório, N. M. R. Peres, Theoretical methods for excitonic physics in 2d materials, physica status solidi (b) 259 (7) (2022) 2200097.  arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/pssb.202200097, doi:https://doi.org/10.1002/pssb.202200097.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/pssb.202200097

  * [6] E. Prada, J. V. Alvarez, K. L. Narasimha-Acharya, F. J. Bailen, J. J. Palacios, Effective-mass theory for the anisotropic exciton in two-dimensional crystals: Application to phosphorene, Phys. Rev. B 91 (2015) 245421.  doi:10.1103/PhysRevB.91.245421.
URL https://link.aps.org/doi/10.1103/PhysRevB.91.245421

  * [7] T. C. Berkelbach, M. S. Hybertsen, D. R. Reichman, Theory of neutral and charged excitons in monolayer transition metal dichalcogenides, Phys. Rev. B 88 (2013) 045318.  doi:10.1103/PhysRevB.88.045318.
URL https://link.aps.org/doi/10.1103/PhysRevB.88.045318

  * [8] A. Franceschetti, H. Fu, L. W. Wang, A. Zunger, Many-body pseudopotential theory of excitons in inp and cdse quantum dots, Phys. Rev. B 60 (1999) 1819–1829.  doi:10.1103/PhysRevB.60.1819.
URL https://link.aps.org/doi/10.1103/PhysRevB.60.1819

  * [9] M. Dvorak, P. Rinke, Dynamical configuration interaction: Quantum embedding that combines wave functions and green’s functions, Phys. Rev. B 99 (2019) 115134.  doi:10.1103/PhysRevB.99.115134.
URL https://link.aps.org/doi/10.1103/PhysRevB.99.115134

  * [10] M. Bieniek, K. Sadecka, L. Szulakowska, P. Hawrylak, Theory of excitons in atomically thin semiconductors: Tight-binding approach, Nanomaterials 12 (9) (2022).  doi:10.3390/nano12091582.
URL https://www.mdpi.com/2079-4991/12/9/1582

  * [11] G. Csanak, H. Taylor, R. Yaris, Green’s function technique in atomic and molecular physics, Vol. 7 of Advances in Atomic and Molecular Physics, Academic Press, 1971, pp. 287–361.  doi:https://doi.org/10.1016/S0065-2199(08)60363-2.
URL https://www.sciencedirect.com/science/article/pii/S0065219908603632

  * [12] G. Strinati, Application of the green’s functions method to the study of the optical properties of semiconductors, La Rivista del Nuovo Cimento (1978-1999) 11 (12) (1988) 1–86.  doi:10.1007/BF02725962.
URL https://doi.org/10.1007/BF02725962

  * [13] M. Rohlfing, S. G. Louie, Electron-hole excitations and optical spectra from first principles, Phys. Rev. B 62 (2000) 4927–4944.  doi:10.1103/PhysRevB.62.4927.
URL https://link.aps.org/doi/10.1103/PhysRevB.62.4927

  * [14] G. Onida, L. Reining, A. Rubio, Electronic excitations: density-functional versus many-body green’s-function approaches, Rev. Mod. Phys. 74 (2002) 601–659.  doi:10.1103/RevModPhys.74.601.
URL https://link.aps.org/doi/10.1103/RevModPhys.74.601

  * [15] C. Attaccalite, M. Grüning, A. Marini, Real-time approach to the optical properties of solids and nanostructures: Time-dependent bethe-salpeter equation, Phys. Rev. B 84 (2011) 245110.  doi:10.1103/PhysRevB.84.245110.
URL https://link.aps.org/doi/10.1103/PhysRevB.84.245110

  * [16] E. Ridolfi, P. E. Trevisanutto, V. M. Pereira, Expeditious computation of nonlinear optical properties of arbitrary order with native electronic interactions in the time domain, Phys. Rev. B 102 (2020) 245110.  doi:10.1103/PhysRevB.102.245110.
URL https://link.aps.org/doi/10.1103/PhysRevB.102.245110

  * [17] G. Cistaro, M. Malakhov, J. J. Esteve-Paredes, A. J. Uría-Álvarez, R. E. F. Silva, F. Martín, J. J. Palacios, A. Picón, Theoretical approach for electron dynamics and ultrafast spectroscopy (edus), Journal of Chemical Theory and Computation 19 (1) (2023) 333–348.  doi:10.1021/acs.jctc.2c00674.
URL https://doi.org/10.1021/acs.jctc.2c00674

  * [18] G. Antonius, S. G. Louie, Theory of exciton-phonon coupling, Phys. Rev. B 105 (2022) 085111.  doi:10.1103/PhysRevB.105.085111.
URL https://link.aps.org/doi/10.1103/PhysRevB.105.085111

  * [19] H.-Y. Chen, D. Sangalli, M. Bernardi, Exciton-phonon interaction and relaxation times from first principles, Phys. Rev. Lett. 125 (2020) 107401.  doi:10.1103/PhysRevLett.125.107401.
URL https://link.aps.org/doi/10.1103/PhysRevLett.125.107401

  * [20] Y.-h. Chan, J. B. Haber, M. H. Naik, J. B. Neaton, D. Y. Qiu, F. H. da Jornada, S. G. Louie, Exciton lifetime and optical line width profile via exciton–phonon interactions: Theory and first-principles calculations for monolayer mos2, Nano Letters 23 (9) (2023) 3971–3977.  doi:10.1021/acs.nanolett.3c00732.
URL https://doi.org/10.1021/acs.nanolett.3c00732

  * [21] M. S. Hybertsen, S. G. Louie, Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies, Phys. Rev. B 34 (1986) 5390–5413.  doi:10.1103/PhysRevB.34.5390.
URL https://link.aps.org/doi/10.1103/PhysRevB.34.5390

  * [22] L. Hedin, On correlation effects in electron spectroscopies and the gw approximation, Journal of Physics: Condensed Matter 11 (42) (1999) R489.  doi:10.1088/0953-8984/11/42/201.
URL https://dx.doi.org/10.1088/0953-8984/11/42/201

  * [23] J. Deslippe, G. Samsonidze, D. A. Strubbe, M. Jain, M. L. Cohen, S. G. Louie, Berkeleygw: A massively parallel computer package for the calculation of the quasiparticle and optical properties of materials and nanostructures, Computer Physics Communications 183 (6) (2012) 1269–1289.  doi:https://doi.org/10.1016/j.cpc.2011.12.006.
URL https://www.sciencedirect.com/science/article/pii/S0010465511003912

  * [24] D. Sangalli, A. Ferretti, H. Miranda, C. Attaccalite, I. Marri, E. Cannuccia, P. Melo, M. Marsili, F. Paleari, A. Marrazzo, G. Prandini, P. Bonfà, M. O. Atambo, F. Affinito, M. Palummo, A. Molina-Sánchez, C. Hogan, M. Grüning, D. Varsano, A. Marini, Many-body perturbation theory calculations using the yambo code, Journal of Physics: Condensed Matter 31 (32) (2019) 325902.  doi:10.1088/1361-648X/ab15d0.
URL https://dx.doi.org/10.1088/1361-648X/ab15d0

  * [25] C. Vorwerk, B. Aurich, C. Cocchi, C. Draxl, Bethe–salpeter equation for absorption and scattering spectroscopy: implementation in the exciting code, Electronic Structure 1 (3) (2019) 037001.
  * [26] E. Perfetto, G. Stefanucci, Cheers: a tool for correlated hole-electron evolution from real-time simulations, Journal of Physics: Condensed Matter 30 (46) (2018) 465901.  doi:10.1088/1361-648X/aae675.
URL https://dx.doi.org/10.1088/1361-648X/aae675

  * [27] R. M. Martin, L. Reining, D. M. Ceperley, Interacting Electrons: Theory and Computational Approaches, Cambridge University Press, 2016.
  * [28] A. Szabo, N. S. Ostlund, Modern quantum chemistry: introduction to advanced electronic structure theory, Courier Corporation, 2012.
  * [29] E. K. Gross, E. Runge, Many-particle theory (1986).
  * [30] I. Ozfidan, M. Korkusinski, A. D. Güçlü, J. A. McGuire, P. Hawrylak, Microscopic theory of the optical properties of colloidal graphene quantum dots, Phys. Rev. B 89 (2014) 085310.  doi:10.1103/PhysRevB.89.085310.
URL https://link.aps.org/doi/10.1103/PhysRevB.89.085310

  * [31] S. Hirata, M. Head-Gordon, Time-dependent density functional theory within the tamm–dancoff approximation, Chemical Physics Letters 314 (3) (1999) 291–299.  doi:https://doi.org/10.1016/S0009-2614(99)01149-5.
URL https://www.sciencedirect.com/science/article/pii/S0009261499011495

  * [32] A. Dreuw, M. Head-Gordon, Single-reference ab initio methods for the calculation of excited states of large molecules, Chemical Reviews 105 (11) (2005) 4009–4037.  doi:10.1021/cr0505627.
URL https://doi.org/10.1021/cr0505627

  * [33] M. S. Hybertsen, S. G. Louie, First-principles theory of quasiparticles: Calculation of band gaps in semiconductors and insulators, Phys. Rev. Lett. 55 (1985) 1418–1421.  doi:10.1103/PhysRevLett.55.1418.
URL https://link.aps.org/doi/10.1103/PhysRevLett.55.1418

  * [34] S. L. Adler, Quantum theory of the dielectric constant in real solids, Phys. Rev. 126 (1962) 413–420.  doi:10.1103/PhysRev.126.413.
URL https://link.aps.org/doi/10.1103/PhysRev.126.413

  * [35] N. Wiser, Dielectric constant with local field effects included, Phys. Rev. 129 (1963) 62–69.  doi:10.1103/PhysRev.129.62.
URL https://link.aps.org/doi/10.1103/PhysRev.129.62

  * [36] F. Wu, F. Qu, A. H. MacDonald, Exciton band structure of monolayer ${\mathrm{mos}}_{2}$, Phys. Rev. B 91 (2015) 075310.  doi:10.1103/PhysRevB.91.075310.
URL https://link.aps.org/doi/10.1103/PhysRevB.91.075310

  * [37] R. Dovesi, A. Erba, R. Orlando, C. M. Zicovich-Wilson, B. Civalleri, L. Maschio, M. Rérat, S. Casassa, J. Baima, S. Salustro, B. Kirtman, Quantum-mechanical condensed matter simulations with crystal, WIREs Computational Molecular Science 8 (4) (2018) e1360.  arXiv:https://wires.onlinelibrary.wiley.com/doi/pdf/10.1002/wcms.1360, doi:https://doi.org/10.1002/wcms.1360.
URL https://wires.onlinelibrary.wiley.com/doi/abs/10.1002/wcms.1360

  * [38] A. L. Fetter, J. D. Walecka, Quantum theory of many-particle systems, Courier Corporation, 2012.
  * [39] T. G. Pedersen, Intraband effects in excitonic second-harmonic generation, Phys. Rev. B 92 (2015) 235432.  doi:10.1103/PhysRevB.92.235432.
URL https://link.aps.org/doi/10.1103/PhysRevB.92.235432

  * [40] J. J. Esteve-Paredes, J. J. Palacios, A comprehensive study of the velocity, momentum and position matrix elements for Bloch states: Application to a local orbital basis, SciPost Phys. Core 6 (2023) 002\.  doi:10.21468/SciPostPhysCore.6.1.002.
URL https://scipost.org/10.21468/SciPostPhysCore.6.1.002

  * [41] C.-C. Lee, Y.-T. Lee, M. Fukuda, T. Ozaki, Tight-binding calculations of optical matrix elements for conductivity using nonorthogonal atomic orbitals: Anomalous hall conductivity in bcc fe, Phys. Rev. B 98 (2018) 115115.  doi:10.1103/PhysRevB.98.115115.
URL https://link.aps.org/doi/10.1103/PhysRevB.98.115115

  * [42] G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, J. Ibañez-Azpiroz, H. Lee, J.-M. Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. Nohara, Y. Nomura, L. Paulatto, S. Poncé, T. Ponweiser, J. Qiao, F. Thöle, S. S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. A. Mostofi, J. R. Yates, Wannier90 as a community code: new features and applications, Journal of Physics: Condensed Matter 32 (16) (2020) 165902.  doi:10.1088/1361-648X/ab51ff.
URL https://dx.doi.org/10.1088/1361-648X/ab51ff

  * [43] J. Ibañez Azpiroz, S. S. Tsirkin, I. Souza, Ab initio calculation of the shift photocurrent by wannier interpolation, Phys. Rev. B 97 (2018) 245143.  doi:10.1103/PhysRevB.97.245143.
URL https://link.aps.org/doi/10.1103/PhysRevB.97.245143

  * [44] M. Wu, S. Chen, S. Camp, K. J. Schafer, M. B. Gaarde, Theory of strong-field attosecond transient absorption, Journal of Physics B: Atomic, Molecular and Optical Physics 49 (6) (2016) 062003.  doi:10.1088/0953-4075/49/6/062003.
URL https://dx.doi.org/10.1088/0953-4075/49/6/062003

  * [45] C. Sanderson, R. Curtin, Armadillo: a template-based c++ library for linear algebra, Journal of Open Source Software 1 (2) (2016) 26.  doi:10.21105/joss.00026.
URL https://doi.org/10.21105/joss.00026

  * [46] E. R. Davidson, The iterative calculation of a few of the lowest eigenvalues and corresponding eigenvectors of large real-symmetric matrices, Journal of Computational Physics 17 (1) (1975) 87–94.  doi:https://doi.org/10.1016/0021-9991(75)90065-0.
URL https://www.sciencedirect.com/science/article/pii/0021999175900650

  * [47] N. S. Rytova, Screened potential of a point charge in a thin film, arXiv e-prints (2018) arXiv:1806.00976arXiv:1806.00976.
  * [48] L. V. Keldysh, Coulomb interaction in thin semiconductor and semimetal films, Soviet Journal of Experimental and Theoretical Physics Letters 29 (1979) 658.
  * [49] K. F. Mak, D. Xiao, J. Shan, Light–valley interactions in 2d semiconductors, Nature Photonics 12 (8) (2018) 451–460.  doi:10.1038/s41566-018-0204-6.
URL https://doi.org/10.1038/s41566-018-0204-6

  * [50] D. Xiao, G.-B. Liu, W. Feng, X. Xu, W. Yao, Coupled spin and valley physics in monolayers of ${\mathrm{mos}}_{2}$ and other group-vi dichalcogenides, Phys. Rev. Lett. 108 (2012) 196802.  doi:10.1103/PhysRevLett.108.196802.
URL https://link.aps.org/doi/10.1103/PhysRevLett.108.196802

  * [51] C. R. Dean, A. F. Young, I. Meric, C. Lee, L. Wang, S. Sorgenfrei, K. Watanabe, T. Taniguchi, P. Kim, K. L. Shepard, J. Hone, Boron nitride substrates for high-quality graphene electronics, Nature Nanotechnology 5 (10) (2010) 722–726.  doi:10.1038/nnano.2010.172.
URL https://doi.org/10.1038/nnano.2010.172

  * [52] F. Zhang, C. S. Ong, J. W. Ruan, M. Wu, X. Q. Shi, Z. K. Tang, S. G. Louie, Intervalley excitonic hybridization, optical selection rules, and imperfect circular dichroism in monolayer $h\text{$-$}\mathrm{BN}$, Phys. Rev. Lett. 128 (2022) 047402.  doi:10.1103/PhysRevLett.128.047402.
URL https://link.aps.org/doi/10.1103/PhysRevLett.128.047402

  * [53] T. Galvani, F. Paleari, H. P. C. Miranda, A. Molina-Sánchez, L. Wirtz, S. Latil, H. Amara, F. m. c. Ducastelle, Excitons in boron nitride single layer, Phys. Rev. B 94 (2016) 125303.  doi:10.1103/PhysRevB.94.125303.
URL https://link.aps.org/doi/10.1103/PhysRevB.94.125303

  * [54] A. V. Krukau, O. A. Vydrov, A. F. Izmaylov, G. E. Scuseria, Influence of the exchange screening parameter on the performance of screened hybrid functionals, The Journal of Chemical Physics 125 (22) (2006) 224106.  arXiv:https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.2404663/13263224/224106\\_1\\_online.pdf, doi:10.1063/1.2404663.
URL https://doi.org/10.1063/1.2404663

  * [55] V. M. Pereira, A. H. Castro Neto, N. M. R. Peres, Tight-binding approach to uniaxial strain in graphene, Phys. Rev. B 80 (2009) 045401.  doi:10.1103/PhysRevB.80.045401.
URL https://link.aps.org/doi/10.1103/PhysRevB.80.045401

  * [56] E. Ridolfi, D. Le, T. S. Rahman, E. R. Mucciolo, C. H. Lewenkopf, A tight-binding model for mos2 monolayers, Journal of Physics: Condensed Matter 27 (36) (2015) 365501.  doi:10.1088/0953-8984/27/36/365501.
URL https://dx.doi.org/10.1088/0953-8984/27/36/365501

  * [57] G. L. Bir, G. E. Pikus, Symmetry and strain-induced effects in semiconductors, Wiley New York, 1974.
  * [58] M. Bieniek, L. Szulakowska, P. Hawrylak, Band nesting and exciton spectrum in monolayer ${\mathrm{mos}}_{2}$, Phys. Rev. B 101 (2020) 125423\.  doi:10.1103/PhysRevB.101.125423.
URL https://link.aps.org/doi/10.1103/PhysRevB.101.125423

  * [59] E. Ridolfi, C. H. Lewenkopf, V. M. Pereira, Excitonic structure of the optical conductivity in ${\mathrm{mos}}_{2}$ monolayers, Phys. Rev. B 97 (2018) 205409.  doi:10.1103/PhysRevB.97.205409.
URL https://link.aps.org/doi/10.1103/PhysRevB.97.205409

  * [60] D. Y. Qiu, F. H. da Jornada, S. G. Louie, Optical spectrum of ${\mathrm{mos}}_{2}$: Many-body effects and diversity of exciton states, Phys. Rev. Lett. 111 (2013) 216805.  doi:10.1103/PhysRevLett.111.216805.
URL https://link.aps.org/doi/10.1103/PhysRevLett.111.216805

  * [61] Y. Li, A. Chernikov, X. Zhang, A. Rigosi, H. M. Hill, A. M. van der Zande, D. A. Chenet, E.-M. Shih, J. Hone, T. F. Heinz, Measurement of the optical dielectric function of monolayer transition-metal dichalcogenides: ${\mathrm{mos}}_{2}$, $\mathrm{Mo}\mathrm{S}{\mathrm{e}}_{2}$, ${\mathrm{ws}}_{2}$, and $\mathrm{WS}{\mathrm{e}}_{2}$, Phys. Rev. B 90 (2014) 205422.  doi:10.1103/PhysRevB.90.205422.
URL https://link.aps.org/doi/10.1103/PhysRevB.90.205422

  * [62] A. García, N. Papior, A. Akhtar, E. Artacho, V. Blum, E. Bosoni, P. Brandimarte, M. Brandbyge, J. I. Cerdá, F. Corsetti, R. Cuadrado, V. Dikan, J. Ferrer, J. Gale, P. García-Fernández, V. M. García-Suárez, S. García, G. Huhs, S. Illera, R. Korytár, P. Koval, I. Lebedeva, L. Lin, P. López-Tarifa, S. G. Mayo, S. Mohr, P. Ordejón, A. Postnikov, Y. Pouillon, M. Pruneda, R. Robles, D. Sánchez-Portal, J. M. Soler, R. Ullah, V. W.-z. Yu, J. Junquera, Siesta: Recent developments and applications, The Journal of Chemical Physics 152 (20) (2020) 204108.  arXiv:https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0005077/16745439/204108\\_1\\_online.pdf, doi:10.1063/5.0005077.
URL https://doi.org/10.1063/5.0005077

  * [63] G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, J. Ibañez-Azpiroz, H. Lee, J.-M. Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. Nohara, Y. Nomura, L. Paulatto, S. Poncé, T. Ponweiser, J. Qiao, F. Thöle, S. S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. A. Mostofi, J. R. Yates, Wannier90 as a community code: new features and applications, Journal of Physics: Condensed Matter 32 (16) (2020) 165902.  doi:10.1088/1361-648X/ab51ff.
URL https://dx.doi.org/10.1088/1361-648X/ab51ff

  * [64] A. C. Dias, J. F. Silveira, F. Qu, Wantibexos: A wannier based tight binding code for electronic band structure, excitonic and optoelectronic properties of solids, Computer Physics Communications 285 (2023) 108636.  doi:https://doi.org/10.1016/j.cpc.2022.108636.
URL https://www.sciencedirect.com/science/article/pii/S0010465522003551