%\addbibresource{ref.bib}

\chapter{Band Structure Methods} \label{chap:theory}

A fundamental understanding of a material’s band structure is essential for predicting and tailoring its electronic and optical properties. In this chapter, we introduce the theoretical and computational methods used to model and analyse the electronic structure of semiconductors, with particular emphasis on first-principles approaches, particularly density functional theory and the GW approximation. Key concepts such as plane-wave basis sets, pseudopotentials, and the virtual crystal approximation are addressed, while we also discuss semi-empirical models, including the Kane model and the $\mathbf{k.p}$ method, which we use for calculating more accurate band structures.



\section{Density Functional Theory} 

Density functional theory (DFT) is a highly useful computational approach that enables the quantitative analysis of various material properties in a purely theoretical manner, without the need for experimental data. DFT methods offer the capability to model materials that are challenging to study through experimentation, or as a means to independently validate experimental findings in less explored material systems. The overall aim of DFT is to solve the time-independent Schroedinger equation:
\begin{equation} \label{eq:Ht}
H \ket{\Psi} = E\ket{\Psi}
\end{equation}

In condensed matter, many properties of materials can be described theoretically by a Hamiltonian ($H$) for a system of nuclei and their electrons:
 \begin{equation} \label{eq:H}
H = -\frac{\hbar^2}{2m_e} \sum_{i} \nabla_{i}^{2}  - \sum_{I} \frac{\hbar^2}{2M_I}  \nabla_{I}^{2}  - \sum_{i,I} \frac{Z_I e^2}{| \textbf{r}_i-\textbf{R}_I |} + \frac{1}{2}  \sum_{i\neq j} \frac{e^2}{| \textbf{r}_i-\textbf{r}_j |}+ \frac{1}{2}  \sum_{I\neq J} \frac{Z_I Z_J e^2}{| \textbf{R}_I-\textbf{R}_J |}
\end{equation}

For such a Hamiltonian in Eq. \ref{eq:Ht}, $E$ represents the eigen-energy of the many-body wave function (\textbf{$\Psi$}) for the system of electrons and nuclei. In Eq. \ref{eq:H}, upper-case subscripts denote nuclei with mass $M_I$ and coordinates $R_I$ with charge $eZ_I$, and lower-case subscripts represents electrons with mass $m_e$ and position $r_i$. The first and second terms represent the kinetic energies of the electrons and nuclei, respectively. The third, fourth and fifth terms represent the electrons-nuclei interactions, electron-electron interactions and the interaction between nuclei, respectively. 

We employ the Born-Oppenheimer approximation, whereby the nuclei masses are assumed to be infinitely heavier than the electrons, thus the electrons respond more rapidly to changes in their surroundings than nuclei. The electrons can therefore be considered to be moving in the field of fixed nuclei. In essence, we are saying that the velocity of the nuclei compared to the electrons is negligible, and therefore the kinetic energy of the nuclei is now zero and the nuclei-nuclei potential energy (which we will now denote as $V_{nn}$) is a constant. This reduces Eq. \ref{eq:H} to Eq. \ref{eq:H1} below:
 \begin{equation} \label{eq:H1}
H = -\frac{\hbar^2}{2m_e} \sum_{i} \nabla_{i}^{2} - \sum_{i,I} \frac{Z_I e^2}{| \textbf{r}_i-\textbf{R}_I |} + \frac{1}{2}  \sum_{i\neq j} \frac{e^2}{| \textbf{r}_i-\textbf{r}_j |}+V_{nn}
\end{equation}

The system can be thought of as a set of $N$ electrons evolving under the influence of an external potential $V_{ext}$, such that:
\begin{equation}
    V_{ext}(\mathbf{r}_i)=\frac{Z_I e^2}{| \mathbf{r}_i-\mathbf{R}_I |}
\end{equation} 

Finally, by defining $U(\mathbf{r}_i,\mathbf{r}_j)=\frac{1}{2} \frac{e^2}{| \mathbf{r}_i-\mathbf{r}_j |}$, we can express the $N$-electron Hamiltonian as:
\begin{equation} \label{eq:H2}
H = \sum_{i}^N \left( -\frac{\hbar^2}{2m_e}  \nabla_{i}^{2} + V_{ext}(\mathbf{r}_i) \right)- \sum_{i\neq j} U(\mathbf{r}_i,\mathbf{r}_j) 
\end{equation}

Finding the eigenvalues and wave functions in such a many-body system is near-impossible for periodic crystals. It is more practical to instead calculate the electronic density, which contains all physically meaningful properties of the crystal. This is computationally possible thanks to the development of DFT, which allows one to consider a many-body interacting electron gas as a single particle moving in an effective potential $V_{eff}$.  It relies on two theorems, which were proved by Hohenberg and Kohn \cite{hohenberg1964inhomogeneous}:


\begin{enumerate}[label=\textbf{Theorem \arabic*:}, ref=\textbf{Theorem \arabic*}, leftmargin=2.7cm]
    \item The external potential $V_{ext}$ acting on the electrons is uniquely determined by the ground-state electron density $n(\mathbf{r})$. Thus, $V_{ext}$ and all other ground state properties can be determined by $n(\mathbf{r})$.
    \item For any external potential $V_{ext}$, there exists a universal energy functional $E[n]$ in terms of the density $n(\mathbf{r})$ whose global minimum is the exact ground state energy of the many-body system.
\end{enumerate}



Given these theorems, Kohn and Sham \cite{kohn1965self} proposed solving the many-body problem by replacing the system with an easier-to-solve auxiliary problem, whereby the ground-state density of the original interacting system is the same as that of non-interacting particles. This leads to the Kohn-Sham equation, which is the non-interacting Schroedinger equation:
 \begin{equation} \label{eq:KS}
H_{aux}\psi_{i}(\textbf{r}) = \epsilon_{i}(\textbf{r})\psi_{i}(\textbf{r}) 
\end{equation}

where $\psi_{i}(\textbf{r}) $ are the single particle orbitals. The auxiliary Hamiltonian $H_{aux}$ is given by 
\begin{equation} \label{Haux}
H_{aux} = \frac{-\hbar^{2}}{2m_{e}}\nabla^{2}+V_{eff}(\textbf{r})
\end{equation}

where $V_{eff}$, known as the Kohn-Sham potential, is an effective one-particle potential felt only by electrons and yields the same density:
\begin{equation} \label{eq:den}
    n(\textbf{r})=\sum\limits_{i=1}^N | \psi_{i}(\textbf{r})|^2
\end{equation}

as the interacting system for electrons at position $\textbf{r}$. If we define the Coulomb interaction energy of the electron density interacting with itself, known as the Hartree potential, $V_H$, as:
\begin{equation}
    V_H \left[ n(\textbf{r}) \right] = e^{2} \int d^{3}r' \frac{n(\textbf{r})}{| \textbf{r}-\textbf{r'} |}
\end{equation}


% define correlation term
and also account for many-body effects with the exchange-correlation functional energy $E_{XC}$, which essentially corresponds to the energy difference between the fully-interacting many-electron system and a non-interacting system of electrons with the same charge density $n(\textbf{r})$, we can define the effective potential as:
\begin{equation} \label{Veff}
V_{eff}(\textbf{r}) = V_{ext}(\textbf{r}) +V_H \left[ n(\textbf{r}) \right] +V_{XC}[n](\textbf{r})
\end{equation}

where $V_{XC}(\textbf{r})=\frac{\delta E_{XC}}{\delta n(\textbf{r})}$ is the potential that defines the exchange and correlation contributions to the single electron equations. A solution to the Kohn-Sham (KS) equations (Eq. \ref{eq:KS}, \ref{Haux}, \ref{eq:den} and \ref{Veff}) is found by performing a self-consistent iterative method:


\begin{enumerate}
    \item The initial trial electron density $\textit{n(\textbf{r})}$ is defined. 
    \item Using this electron density, solve the KS equations to find the single-particle wavefunctions $\psi_{i}(\textbf{r})$.
    \item Calculate the electron density defined by the single-particle wavefunctions $\psi_i$ from step (2), such that $n_{KS}(\textbf{r})=2\sum_i\psi_{i}^*(\textbf{r})\psi_{i}(\textbf{r}) $
    \item If $n_{KS}(\textbf{r})=n(\textbf{r})$, this is the ground state electron density which can be used to calculate the total energy. If $n_{KS}(\textbf{r})\neq n(\textbf{r})$, $n(\textbf{r})$ is updated and the iterative process starts again from step (2).
\end{enumerate}


Ultimately, we are solving Equations \ref{eq:KS}, \ref{Haux}, \ref{eq:den} and \ref{Veff} iteratively until one achieves self-consistency between the density and the potential $V_{eff}$, with the condition that the ground-state density of the original interacting system is equal to that of some exactly soluble system of independent-particles. Typically in any DFT calculation, a self-consistent calculation is initially carried out to determine the electron density of the system. This is followed by a non-self-consistent calculation where the band structure and other material properties are determined.






\subsection{Local Density Approximation}

The exchange-correlation energy $E_{XC}$ is unknown, and therefore must be defined. Only the minimum energy value of the exchange-correlation energy functional has physical meaning. At the minimum, the Kohn-Sham energy functional is equal to the ground-state energy of the system of electrons. Several approaches exist for the approximation of $E_{XC}$, such as the local density approximation (LDA), the generalised gradient approximation (GGA), and the hybrid XC functional. We implement the LDA, whereby the exchange-correlation energy per particle at any particular point in space is approximated as the exchange-correlation energy per particle of a homogenous electron gas with a constant electron density $n$ at all points in space equivalent to the density at this same point. It is a local approximation, where $E_{XC}$ of the LDA form can be written as:
\begin{equation} \label{eq:Elda}
    E_{XC}^{LDA}[n(\textbf{r})] = \int d\textbf{r} \ n(\textbf{r}) \ \epsilon_{XC}[n(\textbf{r})]
\end{equation}

where $\epsilon_{XC}[n(\textbf{r})]$ is the exchange-correlation energy per particle of a homogeneous electron gas of density $n(\textbf{r})$, and $\epsilon_{XC}[n(\textbf{r})]$ is multiplied by a weighted probability $n(\textbf{r})$ that there exists an electron at position $\mathbf{r}$. The exchange-correlation energy is linearly separated into exchange and correlation terms:
\begin{equation}
    \epsilon_{XC}[n(\textbf{r})] = \epsilon_{X}[n(\textbf{r})] + \epsilon_{C}[n(\textbf{r})]
\end{equation}

The exchange energy arises from the Pauli exclusion principle; accounting for the fact that two electrons with the same spin cannot occupy the same quantum state at the same position simultaneously. It is typically derived exactly for a homogeneous electron gas. If we take a homogenous electron gas of $N_e$ electrons interacting in a volume $V$ which has a positive background
charge (keeping the overall charge of the system neutral), the Wigner-Seitz radius $r_S$ is the radius of a sphere whose volume is that of a single electron, such that:
\begin{equation}
    \frac{V}{N_e}=\frac{1}{n}=\frac{4\pi r_S^3}{3}
\end{equation}

Isolating $r_S$ in terms of $n$, the exchange energy per electron can therefore be expressed as:
\begin{equation} \label{eq:ex}
    \epsilon_{X}=\frac{E_{X}}{N_e}=-\frac{1}{4\pi\epsilon_0}\frac{3e^2}{4\pi}\left(\frac{9\pi}{4}\right)^{\frac{1}{3}}\frac{1}{r_S}
\end{equation}

Applying Eq. \ref{eq:ex} to Eq. \ref{eq:Elda} give us Eq. \ref{ExLDA} below, which shows the LDA exchange energy, a more thorough derivation of which can be found in Perdew and Zunger \cite{perdew1981self}. 
\begin{equation} \label{ExLDA}
     E_{X}^{LDA}[n(\textbf{r})] = - \frac{3}{4} \left( \frac{3}{\pi}\right) ^{1/3} \int n(\textbf{r})^{4/3} d\textbf{r}
\end{equation}


Although the integrand of Eq. \ref{ExLDA} allows for an analytical solution for the exchange energy contributions $\epsilon_{X}$, there are no such expressions for $\epsilon_{C}$, the correlation contributions. However, it can be parameterised using highly accurate numerical quantum Monte-Carlo simulations of the homogeneous electron gas. The prevailing parameterisation for the correlation functional is the one introduced by Ref. \cite{perdew1981self} and applied to homogeneous electron gases at different electron densities by Ceperley and Alder \cite{ceperley1980ground}.


The LDA provides precise results for many materials, particularly structural properties such as lattice constants and bond lengths. Although the LDA is computationally efficient and forms the basis for many DFT calculations, it underestimates band gaps in semiconductors, often predicting materials to be more metallic than they are. The correct way to correct for this underestimation is discussed later in the chapter.


% e_c: It includes the effects of electron-electron correlations, electron-electron interactions, and the dynamical behaviour of electrons in the system. 




\subsection{Structural Relaxation}


The results in this thesis rely heavily on alloy scattering effects, which arise when substitutional impurities are introduced into the crystal lattice. These impurities cause local distortions by displacing nearby host atoms, which in turn alters the material's band structure and other properties. In such calculations, accounting for structural relaxation is crucial for accurately capturing these effects.

Structural relaxation is an iterative process used to find the equilibrium atomic positions by minimizing the forces acting on atoms and the stress on the unit cell. This involves calculating the total energy, the forces and the stress tensor for a given configuration, then adjusting the atomic positions until the forces and energy converge below specified thresholds. This is possible thanks to the Hellmann-Feyman theorem \cite{feynman1939forces}, which states that the force on an atom $A$ can be defined as:
\begin{equation}
    F_A=\frac{\partial E}{\partial \mathbf{R}_A}=-\bra{\psi} \frac{\partial H}{\partial \mathbf{R}_A}\ket{\psi}
\end{equation}









\subsection{Plane Waves}

Plane waves are a natural basis to represent single-particle wavefunctions $\psi_i$. In this thesis, we are concerned with the calculation of properties in periodic solids. Bloch's theorem \cite{bloch1929quantenmechanik} provides a framework for understanding how electrons behave in periodic crystal lattices. It states that the wave function $\psi_{nk}(\Vec{r})$ for a bulk material can be expressed as the product of a plane wave term $e^{i\Vec{k}.\Vec{r}}$ and a function $u_{nk}(\Vec{r})$ with the same periodicity as the lattice:
\begin{equation}
    \psi_{nk}(\Vec{r}) = \frac{1}{\sqrt{N_{cell}}}e^{i\Vec{k}.\Vec{r}} u_{nk}(\Vec{r})
\end{equation}
where $N_{cell}$ is the number of primitive cells, $n$ is the $n^{th}$ state corresponding to wavevector $\Vec{k}$, and $\Vec{r}$ is the position vector. If $\Vec{R}$ is a Bravais lattice vector, then:
\begin{equation}
    u_{nk}(\Vec{r}+\Vec{R})=u_{nk}(\Vec{r})
\end{equation}

This implies that states $\Vec{k}$ which lie beyond the first Brillouin zone (BZ) can be translated back into the BZ by adding a suitable reciprocal lattice vector $\Vec{G}$. Electron wave functions of the form $\psi_{nk}(\Vec{r})$ are called Bloch states. If we plot the energy states $E_{nk}$ in a periodic solid as a function of their corresponding wavevectors $\Vec{k}$, this gives the band structure of the solid (which we show in Chapter \ref{chap:gesn}). A proof of Bloch's theorem can be found in Ashcroft and Mermin \cite{ashcroft2022solid}. 
% DO I SHOW BZ IMAGE HERE?

The periodic component $u_{nk}$ can be Fourier transformed and expressed in terms of plane waves of the reciprocal lattice vectors $\mathbf{G}$:
\begin{equation}
    u_{nk}(\textbf{r})=\frac{1}{\sqrt{\Omega_{cell}}} \sum_\mathbf{G} c_{nk} \ e^{i \mathbf{G}.r}
\end{equation}

where $\Omega_{cell}$ is the primitive cell volume, and \textbf{G} satisfies $\mathbf{G}\cdot\mathbf{R}=2\pi \times integer$. Using Bloch's theorem, we can now define a wavefunction of a periodic solid as:
\begin{equation}
    \psi_{nk}(\textbf{r})=\frac{1}{\sqrt{N_{cell}\Omega_{cell}}} \sum_\mathbf{G} c_{nk} \ e^{i (k+\mathbf{G}).r}
\end{equation}


As a result, Eq. \ref{eq:KS}, the Kohn-Sham Equation, can be rewritten as:
\begin{equation} \label{eq:KS1}
    \sum_{\mathbf{G'}}\left[ \frac{\hbar^2}{2m_e}|\mathbf{k}+\mathbf{G} |^2\delta_{\mathbf{G}\mathbf{G'}} +V_{eff}(\mathbf{G}-\mathbf{G'})\right]c_{nk}(\mathbf{G})=\epsilon_{nk} \ c_{nk}(\mathbf{G}) 
\end{equation}

where the Fourier transform of $V_{eff}$ is:
\begin{equation}
    V_{eff}(\mathbf{G})=\frac{1}{\sqrt{\Omega_{cell}}}\int_{\Omega_{cell}}V_{eff}(\mathbf{r})\ e^{-i\mathbf{G}.r}d\mathbf{r}
\end{equation}

The eigenvalues (eigen-energies) $\epsilon_{nk}$ in Eq. \ref{eq:KS1} are discrete, but become a continuum of bands when the $\mathbf{k}$-points lie close to each other. For large systems, $\mathbf{k}$ is a quasi-continuous variable, but the periodic part of the electronic wavefunctions ($u_{nk}$) will be almost identical at $\mathbf{k}$ points very close to each other. This allows the wavefunction to be represented over a small region of reciprocal space by a single $\mathbf{k}$-point. The integration over $\mathbf{k}$-space required to calculate the density can therefore be reduced to a sum over a $\mathbf{k}$-point grid. We can therefore calculate the total energy of the system from electronic states computed on a finite $\mathbf{k}$-point grid, provided it is fine enough to capture variations in
the electron density. The most common method for generating this finite grid is
the Monkhorst and Pack scheme \cite{monkhorst1976special}.





% https://nanohub.org/resources/1524/download/empiricalpseudopotentialmethod_word.pdf
% https://arxiv.org/pdf/1212.6011.pdf
\subsection{Pseudopotentials}

However, the use of plane waves becomes problematic when dealing with regions near the nuclei. Core electrons are relatively unaffected by an atom's chemical environment, with the majority of physical properties of solids being solely dependent on its valence electrons. Using plane-waves to expand the electron wavefunctions which include both the core and valence electrons is excessively time consuming. Moreover, since valence states are orthogonal to core states, valence electron wavefunctions will have a large number of nodes and would require a large number of plane-waves to describe them. It is therefore practical to replace the strong Coulomb potential of the nucleus with a pseudopotential which accounts for the screening and repulsion of the core electrons, therefore replacing the valence electron wave-functions with pseudo-wave-functions. Unlike valence electron wave-functions, the pseudo-wave-functions vary smoothly in the core region. Pseudopotentials essentially reduce the computational cost of solving the Schrödinger equation by approximating electron behaviour in a way that reduces the need to consider the behaviour of the nucleus and inner electrons.

%The exclusion principle demands that valence states must be orthogonal to core states, seeing as the atomic wave-functions are eigenstates of the Hamiltonian. Because of this, the wave-functions of the valence electrons will have a large number of nodes and would require a large number of plane-waves to describe them. Also, due to the core states being localised around the nucleus, the valence states will only preserve this orthogonality with the core electrons by oscillating rapidly in the core region. This leads to a high kinetic energy for the valence electrons in the core region, roughly cancelling out the high potential energy coming from the strong Coulomb potential, so the valence electrons are much more weakly bound than the core electrons. 

\bigskip

Since the crystal potential is periodic, so too is the pseudopotential, and it can be expanded into a Fourier series over the reciprocal lattice. Electronic-structure software employing pseudopotentials (such as ABINIT, VASP or Quantum ESPRESSO) will usually incorporate relativistic effects on core electrons in an approximate manner by using scalar-relativistic pseudopotentials. This is a computationally efficient and quite reliable approach. The most popular pseudopotentials are norm-conserving. We use Hartwigsen-Goedecker-Hutter (HGH) type pseudopotentials \cite{hartwigsen1998relativistic}, which are LDA pseudopotentials and are norm-conserving.

% https://nanohub.org/resources/1524/download/empiricalpseudopotentialmethod_word.pdf
% https://arxiv.org/pdf/1212.6011.pdf







\subsection{ABINIT}

ABINIT \cite{ab} is an open-source software package used to perform first-principles electronic structure calculations on materials. These first-principles calculations are based on DFT, using pseudopotentials and plane-wave basis sets to represent the electronic structures of materials. For band structure calculations on ABNIIT, a self-consistent field (SCF) calculation is performed whereby the density is outputted, followed by a second calculation which uses this density in a non-SCF band structure calculation. ABINIT is capable of predicting the electronic structures, crystal structures, optical properties, and vibrational properties of materials, and can also take advantage of parallel computing. Its user-friendly interface and extensive documentation makes it one of the most used packages in DFT. In this thesis, we will use ABINIT within the DFT framework to solve the KS equations iteratively for a system of atoms from first principles. 


\bigskip

ABINIT also includes a post-processor code called "Cut 3-Dimensional" (Cut3d), which analyses any file produced by ABINIT that contains 3D real space data, such as potential files, density files and wavefunction files. In the case of wavefunction files, analysis of only one wave function at a time may be performed, where the user must specify the k-point, band number and (where necessary) the spinor-component. This was necessary for several calculations in this work, including determining the momentum matrix elements, which are used in the calculation of optical properties, and wavefunction overlaps, which are essential in computing scattering matrix elements.



\subsubsection{Convergence Testing in ABINIT}

In DFT calculations using plane-wave basis sets, as implemented in ABINIT, the accuracy of the computed results is closely tied to the choice of the plane-wave energy cut-off parameter, denoted as \( E_{\text{cut}} \). The plane-wave energy cut-off defines the maximum kinetic energy of the plane waves included in the basis set used to expand the electronic wavefunctions. Since DFT calculations involve solving the Kohn-Sham equations self-consistently, it is crucial to ensure that the truncation of the basis set does not compromise the reliability of the results. The plane-wave basis set employed in ABINIT is truncated by restricting the kinetic energy of the included plane waves to a value below  $E_{\text{cut}}$. In other words, if we recall Eq. \ref{eq:KS1}, the kinetic energy term has upper bound:
\begin{equation}
    \frac{\hbar^2}{2m_e}|\mathbf{k}+\mathbf{G} |^2<E_{cut}
\end{equation}

for all $\mathbf{k}$ in the Brillouin Zone. A higher cut-off energy allows the inclusion of more plane waves, resulting in a more accurate description of the electronic wavefunctions and, consequently, more precise total energies, forces, and stress tensors. However, an excessively high cut-off can lead to unnecessary computational cost without significant improvement in accuracy, since plane-waves with higher energies contribute less to the wavefunction. Therefore, it is essential to perform convergence tests to determine the minimum \( E_{\text{cut}} \) that yields converged results within an acceptable tolerance for the system under investigation. Typically, the total energy or lattice parameter is calculated at incrementally increased values of \( E_{\text{cut}} \) until the differences between successive calculations are within a chosen tolerance. 



Performing plane-wave pseudopotential calculations also requires defining the k-point grid to sample the Brillouin zone effectively. The \textit{ngkpt} variable in ABINIT is used to define a Monkhorst-Pack k-point grid, which defines the number of k-points along each axis in the reciprocal space grid. For instance, a value of {4x4x4} creates a grid with four k-points along each of the three lattice directions (giving a total of $4^3=64$ points). It is once again essential to establish clear convergence criteria in a similar manner to \( E_{\text{cut}} \). The accurate sampling of k-points plays a crucial role in ensuring the accuracy and efficiency of these computations.


\section{GW Method} 

Whilst DFT gives an accurate total energy for the many-body system, and for many systems treated in this work the bands produced are a very good approximation for calculating different material properties, the energy eigenvalues of the KS Hamiltonian are not the correct excitation energies of the system. The eigenvalues produced are those that correspond to a one-electron potential that has the same total density as the many-body system. More precisely, the exchange-correlation functionals describing the interactions between electrons often fail to capture the electron-electron interactions accurately, which leads to a bandgap underestimation. Therefore, properties such as the band gap are not well produced and require another method of calculation to obtain true accuracy. 


%\bigskip
One such method is GW, which operates by finding the response of a system that is perturbed by transitioning an electron to higher energy states. These energy transitions correspond to the observable energies. The name GW stems from the Green's function $G$, which describes the propagation of an electron in a solid, and the screened Coulomb interaction $W$ which represents the interaction between electrons, taking into account the screening effect due to other electrons. 

Using ABINIT, we implement the single-shot GW method ($G_0W_0$), which begins with a zeroth-order approximation (a DFT calculation) of the wavefunctions and energies. This is followed by a calculation of the screening (KS susceptibility $\chi^0$ and then inverse dielectric matrix $\epsilon^{-1}$), such that:
\begin{equation}
    W(\mathbf{r}, \mathbf{r}'; \omega) = \int d\mathbf{r}'' \, \epsilon^{-1}(\mathbf{r}, \mathbf{r}''; \omega) \, v(\mathbf{r}'' - \mathbf{r}')
\end{equation}
where $v$ is the bare Coulomb interaction. The KS band structure and $W$ are then used to evaluate the matrix elements of the self-energy $\Sigma=iGW$, finally obtaining the quasi-particle corrections. In this process, we check for convergence of the screening with respect to both the number of bands and plane-waves, as well as convergence of the self-energy with respect to the number of bands.

GW is very computationally expensive, and is best used to find bandgaps for a small number of k-points rather than to produce an entire band structure. While the GW method captures the different many-body effects at the $\Gamma$ and $L$ points in GeSn, a simple scissors shift applies a uniform correction to all states and cannot correct the relative difference between these gaps, which DFT miscalculates. An excellent review of the mathematics behind GW can be found in Ref. \cite{GW}.







\section{Virtual Crystal Approximation} 

The virtual crystal approximation (VCA) is a simple and yet effective theoretical framework used for studying alloys. Its purpose is due to the fact that electronic structures become more complex when impurities are introduced into the crystal lattice, as in any alloy. The VCA is an interpolation of the behaviours of the constituent atoms, and can be used in determining many properties of an alloy, including the bulk modulus or lattice constant.

In electronic structure calculations, we implement the alchemical potential approximation (APA), which is essentially the VCA applied to potentials. The APA works by treating each atom in the system as an "alloy" by averaging the potentials $V$ of the individual constituents. In the case of a GeSn alloy of Sn concentration $x$:
\begin{equation*}
    V_{GeSn}=(1-x)V_{Ge}+xV_{Sn}
\end{equation*}
where $V_{Ge}$ and $V_{Sn}$ are the potentials of Ge and Sn respectively. This approximation essentially smooths out the abrupt changes in the potential energy landscape created by the atomic disorder, making it simpler to solve the quantum mechanical equations that describe the electronic structure. By solving the Schrödinger equation with the effective potential and composition, one can calculate the electronic band structure of the virtual crystal. Pseudopotentials are often used in the APA.

However, as we will see in later chapters, the APA does not work when calculating elastic scattering parameters, since it doesn't capture certain effects such as local atomic fluctuations or short-range order. Methods beyond the APA are required to find the elastic scattering contributions.




\section{Kane Model}  \label{sec:Kane}

Given the underestimation of the band energies using the LDA and the high computational cost of the GW method, an alternative approach is needed to accurately model the band structure dispersion near the bandgap for multiple $\mathbf{k}$-points. 

A suitable theoretical framework for such an issue is the Kane model \cite{kane1957band} which, in bulk (3D) semiconductors, provides a compact yet effective way to capture the energy-momentum relationship near the Brillouin zone centre (the $\Gamma$ point). The Kane model centres around constructing a Kane Hamiltonian which couples the material's conduction and valence bands through momentum matrix elements. A four-band Kane model is commonly used (eight-band when SOC is considered), which considers one conduction band state and three valence band states (light-hole, heavy-hole, and split-off bands). In this case, the model is constructed as a $4\times 4$ matrix Hamiltonian (below) describing the coupling between these states, taking into account the $3D$ crystal momentum: 

\[
\mathbf{H}(k_x, k_y, k_z) = \begin{bmatrix}
E_c & b k_x & b k_y & b k_z \\
b k_x & E_v & 0 & 0 \\
b k_y & 0 & E_v & 0 \\
b k_z & 0 & 0 & E_v
\end{bmatrix}
\]

%In order to accurately determine $E_k$, the conduction band energies around the $\Gamma$-point, we first use DFT to find the band dispersion around the $\Gamma$-point. We then fit the DFT eigen-energies to a $3D$-Kane model \cite{kane1957band}, the Hamiltonian $\textbf{H}$ for which is described generally using the matrix below for direct gap $E_g$ and parameters $a$ and $b$:


In this Hamiltonian, the conduction band energy $E_c=a (k_x^2 + k_y^2 + k_z^2) + E_g$, where $E_g$ is the direct gap, and the valence band energy $E_v=-a (k_x^2 + k_y^2 + k_z^2)$, where the parameter $a\approx\frac{\hbar^2}{2m^*}$. The parameter $b=\frac{\hbar^2}{2m^*}|p|$ represents the coupling between the conduction and valence bands, typically involving the momentum matrix elements \textbf{p}. The calculated DFT eigen-energies are fitted to the eigenvalues of this matrix, whereby we solve for $a$ and $b$.  

Unlike the scissors shift, which applies a rigid bandgap correction, the Kane model accounts for band coupling between the conduction and valence bands and resulting non-parabolicity, providing a more physically accurate description of the electronic structure in narrow-gap semiconductors. This is essential for describing the effective mass and mobility of electrons in 3D systems.



\section{k.p Theory} % NEED TO FINISH THIS< INCLUDE MATHEMATICS

\textbf{k.p} theory is based on time-independent perturbation theory, employed to study how the electronic structure and energy levels of a material change in response to external perturbations (such as an external electric field or strain) with various levels of refinement. For example, the number of basis states (bands) can be varied
(e.g. single-band, 8-band) \cite{chaudhuri2021multiscale}. 

\textbf{k.p} theory operates by expanding the electronic wave functions and Hamiltonian (i.e. calculating the eigenvalues and eigenvectors) around a reference $\mathbf{k}$-point ($\mathbf{k_0}$), usually a high-symmetry point in the BZ such as the $\Gamma$-point. The resulting Hamiltonian is a sum of the free-electron term $H_0$ and the $\mathbf{k}$-dependent term $H_k$, such that:
\begin{align}
  H_0&=\frac{\mathbf{p}^2}{2m}+V  \nonumber \\ 
  H_k&=\frac{\hbar}{m}(\mathbf{k}.\mathbf{p})+\frac{\hbar^2\mathbf{k}^2}{2m}
\end{align}
where $\mathbf{k}$ is a simple vector that characterises the periodicity of the wavefunction, and the momentum operator $\mathbf{p}\neq \hbar \mathbf{k}$. Thus, the band structure over the entire BZ can be extrapolated from the BZ centre energy gaps and optical matrix elements \cite{peter2010fundamentals}. Following Schulz and O'Reilly \cite{ScOR2017_CRC}, applying this Hamiltonian to a Bloch wavefunction, second-order perturbation theory can be used to evaluate the variation of the energy levels $E_n(\mathbf{k})$ at band $n$ with wavevector $\mathbf{k}$ close to the reference point $\mathbf{k_0}$, shown in Eq. \ref{eq:En} below:
\begin{equation} \label{eq:En}
    E_n(\mathbf{k})=E_n(\mathbf{k_0})+\frac{\hbar}{m}\mathbf{k}.\mathbf{p}_{nn}+\frac{\hbar^2\mathbf{k}^2}{2m}+\frac{\hbar^2}{m^2}\sum_{n\neq n'}\frac{|\mathbf{k}.\mathbf{p}_{nn'}|^2}{E_n(\mathbf{k_0})-E_{n'}(\mathbf{k_0})}
\end{equation}
where $\mathbf{p}_{nn'}$ denotes the momentum matrix element between the $nth$ and $n'th$ states at $\mathbf{k_0}$.


The \textbf{k.p} method is particularly useful in estimating the effective masses of bands by perturbatively calculating the energy dispersion about the known BZ point. We exploit the fact that all of the important information in GeSn can be extracted from the $\Gamma$ point. Using this, we can use a \textbf{k.p} model to correct the DFT band structure for GW using the momentum matrix elements at $\Gamma$. A detailed description of how \textbf{k.p} theory can be found in Ref. \cite{kodama2016kp}.










