%\part{Theory}

\chapter{Electronic structure theory}
\label{chp:ElecStruc}

\section{Introduction}
\label{sec:ElecStrucIntro}

The radiative and nonradiative properties of, for instance, the nitride QWs targeted in this work will depend on their underlying electronic structure. In general, the realistic description of the electronic structure of a material or nanostructure is a formidable task. Determining the electronic structure essentially amounts to solving the \sch~equation for the system of electrons including their interaction with each other and the lattice. Generally, there are around $10^{23}$ valence electrons involved in bonding in a cubic centimetre of a solid, presenting a significant, if not insurmountable, computational task~\cite{QuantumTheoryOfSolidsEoin}.

However, methods to reduce the computational complexity of this problem while still capturing the essential physics have been developed. Approaches to determining the electronic structure are broadly divided into two categories, \emph{ab initio} and semi-empirical methods. \emph{Ab initio} (or first-principles) methods target solutions to the \sch~equation that, at least in principle, do not require the fitting of model parameters to experimental results. Semi-empirical models, in contrast, will require the fitting of input parameters to either experimental or first-principles results. In this chapter we will discuss state-of-the-art \emph{ab initio} methods used in electronic structure calculations of III-N materials and to develop semi-empirical models. Subsequently, we will outline our tight-binding model in more detail. Before proceeding, however, we first revise some foundational concepts underlying the description of electronic states which will be useful in our subsequent discussions.

\section{Basic considerations}
\label{sec:ElecStrucBasicConsiderations}

To begin, we will look at the so-called reciprocal lattice and some preliminary approximations often applied when determining the electronic structures of solids. Next, we will introduce the idea of an energy band and Bloch's theorem, including some consequences of this theorem.

\subsection{Reciprocal space}

In Chapter~\ref{chp:LEDs} above we described the wurtzite crystal structure and some of the experimentally relevant directions and crystallographic planes, e.g., the $c$-plane, the $a$-plane and the $m$-plane. These planes are categorised by Miller-Bravais indices which are related to the reciprocal of the intersections of the plane with the Cartesian axes. Experimentally, information about the crystallographic planes can be extracted by diffraction measurements, and can be used to determine information about the underlying material properties, such as lattice constants~\cite{UeYo1994}. Diffraction is related to the so-called ``reciprocal'' lattice, which in turn is related to the spatial frequency of points in the real space lattice. 

Given an ideal real space lattice with primitive lattice vectors $\mathbf{a}_{j}$, we recall that any point in the lattice, $\mathbf{R}$, is described by
\begin{equation*}
    \mathbf{R} = n_{1}\mathbf{a}_{1} + n_{2}\mathbf{a}_{2} + n_{3}\mathbf{a}_{3} \,\, ,
\end{equation*}
with $n_{i} \in \mathbb{Z}$. The associated reciprocal lattice is described by vectors $\mathbf{G}$ such that~\cite{Kittel,HaugAndKochThird}
\begin{equation}
    e^{i\mathbf{G}\cdot\mathbf{R}} = 1 \,\, ,
    \label{eq:reciprocalCondition}
\end{equation}
where $\mathbf{G} = m_{1}\mathbf{b}_{1} + m_{2}\mathbf{b}_{2} + m_{3}\mathbf{b}_{3}$ and $m_{i} \in \mathbb{Z}$. The condition given in Eq.~\eqref{eq:reciprocalCondition} requires that $\mathbf{b}_{i}\cdot\mathbf{a}_{j} = 2\pi\delta_{ij}$, with $i,j \in \{1,2,3\}$~\cite{Kittel}. Thus, (with $i\neq j, k$) $\mathbf{b}_{i}$ is perpendicular to both $\mathbf{a}_{j}$ and $\mathbf{a}_{k}$, and we can write 
\begin{equation*}
    \mathbf{b}_{i} = C(\mathbf{a}_{j}\times\mathbf{a}_{k}) \,\, ,
\end{equation*}
where $C$ is some constant. Making use of the requirement that $\mathbf{a}_{i}\cdot\mathbf{b}_{i} = 2\pi$,
\begin{equation*}
    \mathbf{a}_{i}\cdot\mathbf{b}_{i} = \mathbf{a}_{i}\cdot(C(\mathbf{a}_{j}\times\mathbf{a}_{k})) = 2\pi \,\, ,
\end{equation*}
and so
\begin{equation*}
    C = \frac{2\pi}{\mathbf{a}_{i}\cdot(\mathbf{a}_{j}\times\mathbf{a}_{k})} \,\, .
\end{equation*}
Hence, in 3D, the vectors $\mathbf{b}_{i}$ are given by~\cite{Kittel,HaugAndKochThird}
\begin{equation*}
    \mathbf{b}_{i} = 2\pi\frac{\mathbf{a}_{j}\times\mathbf{a}_{k}}{\mathbf{a}_{i}\cdot\left(\mathbf{a}_{j}\times\mathbf{a}_{k}\right)} \,\, .
\end{equation*}
The vectors $\mathbf{b}_{i}$ define the reciprocal lattice and span the reciprocal space. They can be thought of as a discrete 3D Fourier transform of the real space (sometimes called direct) lattice, defined by the vectors $\mathbf{a}_{i}$. An arbitrary vector in the reciprocal space is often represented by the vector $\mathbf{k}$.

\begin{figure}[t!]
\centering
    \includegraphics[width=0.75\columnwidth]{Figures/HEXFBZ.eps}
    \caption{Schematic of the first Brillouin zone for a hexagonal lattice spanned by the reciprocal lattice vectors $\mathbf{b}_{i}$. High symmetry points, including the $\Gamma$-point at $\mathbf{k} = \mathbf{0}$ are indicated.}
    \label{fig:HexFBZ}
\end{figure}

In the real space lattice it is possible to define a primitive unit cell as the parallelepiped defined by the primitive lattice vectors, $\mathbf{a}_{i}$, which is the minimum repeatable volume that generates the lattice. This primitive unit cell is called the Wigner-Seitz cell. Similarly, in the reciprocal lattice one may define such a unit cell using the primitive reciprocal lattice vectors, $\mathbf{b}_{i}$. In the reciprocal lattice, this unit cell is called the first Brillouin zone (FBZ) and is depicted for a hexagonal real space lattice in Fig.~\ref{fig:HexFBZ}~\cite{Kittel}. These ideas will be useful in the next section, as we return to the overall problem at hand, solving the \sch~equation.

\subsection{Electronic band structures}

As previously stated, in order to determine the electronic structure of a crystal, we must solve the \sch~equation to find the energy levels and wave functions of the electrons. To do so, the form of the Hamiltonian operator, $\hat{H}$ must be known. In principle, the full Hamiltonian for a system of electrons in a crystal lattice must take into account interactions between the electrons as well as interactions between the ions of the lattice, and the interactions between electrons and ions of the lattice~\cite{AdvancedSolidStatePhysics}:
\begin{align}
    \hat{H} =& \overbrace{-\frac{\hbar^{2}}{2m}\sum_{i}^{n}\nabla_{i}^{2}}^{\hat{T}_{\text{e}}}
    \overbrace{-\frac{\hbar^{2}}{2M}\sum_{\alpha}^{N}\nabla_{\alpha}^{2}}^{\hat{T}_{\text{I}}}
    \overbrace{-\frac{1}{4\pi\epsilon_{0}}\sum_{i}^{n}\sum_{\alpha}^{N}\frac{Ze^{2}}{|\mathbf{r}_{i}-\mathbf{R}_{\alpha}|}}^{\hat{V}_{\text{eI}}}\nonumber \\
    &\underbrace{+\frac{1}{8\pi\epsilon_{0}}\sum_{i}^{n}\sum_{j\neq i}^{n}\frac{e^{2}}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}}_{\hat{V}_{\text{ee}}}
    \underbrace{+\frac{1}{8\pi\epsilon_{0}}\sum_{\alpha}^{N}\sum_{\beta\neq\alpha}^{N}\frac{(Ze)^{2}}{|\mathbf{R}_{\alpha}-\mathbf{R}_{\beta}|}}_{\hat{V}_{\text{II}}}\,\, .
    \label{eq:fullHamiltonian}
\end{align}
For simplicity, we have assumed here that there is only one species of ion in the lattice, with proton number $Z$. Furthermore, we consider $n$ electrons of mass $m$, at positions $\mathbf{r}_{i}$, and $N$ ions of mass $M$ at positions $\mathbf{R}_{\alpha}$. Electrons are indexed by Latin labels, and ions by Greek labels. Here we have neglected relativistic effects, e.g., spin-orbit coupling, which is of secondary importance in III-N semiconductors~\cite{QuantumTheoryOfSolidsEoin}. The Hamiltonian can thus be represented as 
\begin{equation}
    \hat{H} = \hat{T}_{\text{e}} + \hat{T}_{\text{I}} + \hat{V}_{\text{eI}} + \hat{V}_{\text{ee}} + \hat{V}_{\text{II}} \,\, ,
    \label{eq:HamiltonianOperatorRepresentation}
\end{equation}
where $\hat{T}_{\text{e}}$ is the kinetic energy operator for the electrons, $\hat{T}_{\text{I}}$ is the kinetic energy operator for the ions, $\hat{V}_{\text{eI}}$ is the potential operator for the attractive electron-ion interaction, $\hat{V}_{\text{ee}}$ is the potential operator for the repulsive electron-electron interaction and $\hat{V}_{\text{II}}$ is the potential operator for the repulsive ion-ion interaction. 

Given that the mass of the ions is much larger than the mass of the electrons (i.e. $M \gg m$), we employ here the Born-Oppenheimer approximation to decouple the electronic system from the ionic part~\cite{BoOp1927}. Thus, we neglect the $\hat{T}_{\text{I}}$ term and recast the problem as that of a number of electrons moving in an effective potential due to the ions. To do so, we combine the $\hat{V}_{\text{eI}}$ and $\hat{V}_{\text{II}}$ terms into $\hat{V}_{\text{eff}}$, leaving a separate interaction potential term that describes the interaction of the electrons with each other~\cite{AdvancedSolidStatePhysics}. Thus, Eq.~\eqref{eq:HamiltonianOperatorRepresentation} is reduced to~\cite{ModernManyParticlePhysics}
\begin{equation}
    \hat{H} = \hat{T}_{\text{e}} + \hat{V}_{\text{eff}} + \hat{V}_{\text{ee}} = \sum_{i}^{n}\left(-\frac{\hbar^{2}}{2m}\nabla_{i}^{2} + V_{\text{eff}}\left(\mathbf{r}_{i}\right)\right) + \sum_{i}^{n}\sum_{j\neq i}^{n}V\left(\mathbf{r}_{ij}\right) \,\, .
    \label{eq:manyBodyHamiltonian}
\end{equation}

For the purpose of discussion, for now we will neglect the Coulomb interaction between electrons and assume that the impact of these interactions is accounted for in the effective potential. In this simplified picture, the effective potential is static, and so we must solve the time independent \sch~equation,
\begin{equation}
    \label{eq:Schrod}
    \hat{H}\Psi_{n} = \text{E}_{n}\Psi_{n} \,\, ,
\end{equation}
to find the energy levels, $\text{E}_{n}$, and wave function, $\Psi_{n}$, of the system. When the interactions between electrons are neglected, the wave function for the system of electrons, $\Psi_{n}$, can be written as the product of the single particle wave functions, $\psi_{n,i}$, that individually solve Eq.~\eqref{eq:Schrod}. The energy of the system, $\text{E}_{n}$ is then the sum of the single particle energies, $\text{E}_{n,i}$. Even in this simplified picture, solving the single particle \sch~equation over the length scale of the entire solid is a significant challenge. However, the effective potential above will have the periodicity of the crystal lattice, thus we may make use of Bloch's theorem~\cite{Kittel,Bloch1929}. 

Bloch's theorem applies to perfectly crystalline materials in which the translational invariance of the potential is not violated. In such a material, for some arbitrary real space vector $\mathbf{r}$, and a lattice vector $\mathbf{R}$, defined above, we would have
\begin{equation*}
    V_{\text{eff}}\left(\mathbf{r} + \mathbf{R}\right) = V_{\text{eff}}\left(\mathbf{r}\right) \,\, .
\end{equation*}
This being the case, Bloch's theorem states that the single particle electronic wave functions, $\psi_{n}$, will be the product of a plane wave, $e^{i\mathbf{k}\cdot\mathbf{r}}$, with wave vector $\mathbf{k}$ and a function with the same periodicity as the lattice (called a Bloch factor), $u_{\mathbf{k}}(\mathbf{r})$, i.e.,
\begin{equation*}
    \psi_{n,\mathbf{k}}\left(\mathbf{r}\right) = e^{i\mathbf{k}\cdot\mathbf{r}} u_{n,\mathbf{k}}\left(\mathbf{r}\right)\,\, ,
\end{equation*}
where, for a given $\mathbf{k}$,
\begin{equation*}
    u_{n,\mathbf{k}}\left(\mathbf{r} + \mathbf{R}\right) = u_{n,\mathbf{k}}\left(\mathbf{r}\right) \,\, .
\end{equation*}
Thus
\begin{equation*}
    \psi_{n,\mathbf{k}}\left(\mathbf{r} + \mathbf{R}\right) = e^{i\mathbf{k}\cdot\mathbf{R}}\psi_{n,\mathbf{k}}\left(\mathbf{r}\right)\,\, ,
\end{equation*}
i.e., the electronic wave function is periodic with the periodicity of the lattice, up to a phase factor of $e^{i\mathbf{k}\cdot\mathbf{R}}$, and, as a result, the charge density, $|\psi_{n,\mathbf{k}}|^{2}$, will be the same at $\mathbf{r}$ and $\mathbf{r}+\mathbf{R}$. A wave function that fulfills Bloch's theorem, i.e., has the form of a plane wave times a Bloch factor, is called a Bloch wave function. Thus, if we know the nature of the (Bloch) wave function within a Wigner-Seitz cell, we could extend that knowledge to encompass the whole lattice, while taking the plane wave phase factor into account. From Bloch's theorem, a wave vector $\mathbf{k}$ can be associated with each energy state, $\text{E}_{n}$. So, for a given value of $\mathbf{k}$, we will find a (discrete) set of allowed energies, $\text{E}_{n}$ by solving the \sch~equation (Eq.~\eqref{eq:Schrod}) assuming a Bloch wave function. 

\begin{figure}[t!]
\centering
    \includegraphics[width=\columnwidth]{Figures/TBGaNBandStructure.eps}
    \caption{Band structure of wurtzite GaN determined in a tight-binding model. Tight-binding parameters from Ref.~\cite{CaSc2013local} have been used.}
    \label{fig:TBBandsGaN}
\end{figure}

Furthermore, due to the periodicity of the lattice, it can be shown that the wave vectors, $\mathbf{k}$, are not uniquely determined. To illustrate this, we take the approach used in Ref.~\cite{HodmasElectronicBands} and consider a 1D periodic lattice with a unit cell of length $L$, thus
\begin{equation*}
    \psi_{n,k}(x) = e^{ikx}u_{n,k}(x) \,\, ,
\end{equation*}
and $u_{n,k}(x) = u_{n,k}(x + L)$. Multiplying by a plane wave with the periodicity of the lattice, $e^{i2\pi mx/L}$, where $m$ is any integer, and its complex conjugate, we have
\begin{align*}
    \psi_{n,k}(x) &= e^{ikx}e^{i\frac{2\pi m}{L}x}e^{-i\frac{2\pi m}{L}x}u_{n,k}(x) \\
    &= e^{i(k+\frac{2\pi m}{L})x}\left(e^{-i\frac{2\pi m}{L}x}u_{n,k}(x)\right)\,\, .
\end{align*}
The term $2\pi m/L$ can be identified as a 1D reciprocal lattice vector, denoted $G_{m}$, as it satisfies Eq.~\eqref{eq:reciprocalCondition}. Thus, the term $\exp(-iG_{m}x)u_{k}(x)$ also has the periodicity of the lattice and so the properties of $\psi_{n,k}$ are equivalent at $k$ and $k+G_{m}$ (up to a phase factor): e.g., for the energy, $\text{E}_{n}(k)$, of the state $\psi_{n}$ at a given point $k$ in reciprocal space, we have that $\text{E}_{n}(k+G_{m}) = \text{E}(k)$. This allows us to restrict our analysis to the interval $-\frac{\pi}{L} < k < \frac{\pi}{L}$, which corresponds to the FBZ in this 1D example. This is the so-called ``reduced zone scheme'' where any value of $k$ outside the FBZ is translated back to the FBZ by suitable choice of $G_{m}$. The variation of $\text{E}_{n}$ with $k$ gives rise to so-called energy ``bands'' for a given value of $n$. A plot of all the energy bands (i.e. for all values of $n$) gives rise to the ``band structure'' of a material. However, visualising 3D reciprocal space is generally more difficult than in the 1D example above. Thus, we select certain ``high symmetry'' points in the FBZ to gain insight into the 3D band structure. For example, the point $\mathbf{k}=\mathbf{0}$ at the centre of the FBZ is a high symmetry point often referred to as the $\Gamma$-point. Other so-called high symmetry points in the first Brillouin zone include A, L, M, H and K (see, for instance, Ref.~\cite{KoSa1983}, and Fig.~\ref{fig:HexFBZ} above). An example of a band structure is shown in Fig.~\ref{fig:TBBandsGaN} for GaN. Moving along the $x$-axis of Fig.~\ref{fig:TBBandsGaN} corresponds to moving along a straight line in reciprocal space between high symmetry points. 

Above, we briefly discussed the independent particle model, in which the interaction between electrons is neglected and the solution wave function is a combination of the single particle (i.e. independent particle) wave functions. We now turn to discuss more advanced models that aim to include the electron-electron interaction.

%However, due to the Fermionic nature of electrons, the constructed wave function, $\Psi$, must be anti-symmetric, which can be achieved by constructing a Slater determinant from the single particle wave functions, $\psi_{i}$, via~\cite{ModernManyParticlePhysics}
%\begin{equation*}
%    \Psi(\mathbf{r}_{1},\mathbf{r}_{2},\cdots,\mathbf{r}_{N}) =  
%    \frac{1}{\sqrt{N!}}\det
%    \begin{pmatrix}
%    \psi_{1}(\mathbf{r}_{1}) & \psi_{2}(\mathbf{r}_{1}) & \cdots & \psi_{N}(\mathbf{r}_{1}) \\
%    \psi_{1}(\mathbf{r}_{2}) & \psi_{2}(\mathbf{r}_{2}) & \cdots & \psi_{N}(\mathbf{r}_{2}) \\
%    \vdots & \cdots & \ddots & \vdots \\
%    \psi_{1}(\mathbf{r}_{N}) & \psi_{2}(\mathbf{r}_{N}) & \cdots & \psi_{N}(\mathbf{r}_{N})
%    \end{pmatrix}
%    \,\, .
%\end{equation*}

\section{\textit{Ab initio} methods}
\label{sec:ElecStrucAbInitio}

In the single particle picture described above we noted that the overall solution wave function of Eq.~\eqref{eq:Schrod}, $\Psi_{n}$, would be the product of the single particle wave functions, $\psi_{n,i}$. However, due to the Fermionic nature of electrons, $\Psi_{n}$ must be anti-symmetric, which can be achieved by constructing a Slater determinant from the single particle wave functions. For a system of $N$ electrons at positions $\mathbf{r}_{\alpha}$, with $\alpha \in \{1, \ldots , N\}$~\cite{ModernManyParticlePhysics}:
\begin{equation}
    \label{eq:SlaterDet}
    \Psi_{n}(\mathbf{r}_{1},\mathbf{r}_{2},\ldots,\mathbf{r}_{N}) =  
    \frac{1}{\sqrt{N!}}\det
    \begin{pmatrix}
        \psi_{n,1}(\mathbf{r}_{1}) & \psi_{n,2}(\mathbf{r}_{1}) & \cdots & \psi_{n,N}(\mathbf{r}_{1}) \\
        \psi_{n,1}(\mathbf{r}_{2}) & \psi_{n,2}(\mathbf{r}_{2}) & \cdots & \psi_{n,N}(\mathbf{r}_{2}) \\
        \vdots & \cdots & \ddots & \vdots \\
        \psi_{n,1}(\mathbf{r}_{N}) & \psi_{n,2}(\mathbf{r}_{N}) & \cdots & \psi_{n,N}(\mathbf{r}_{N})
    \end{pmatrix}
    \,\, .
\end{equation}
However, the true many-body wave function is a linear combination of Slater determinants. As such, we must account for the quantum mechanical nature of the wave function including exchange and correlation energies and effects. Generally, the accurate calculation of solutions to the \sch~equation in the presence of interacting electrons is a significant challenge. Even so, the groundbreaking work of Hohenberg and Kohn~\cite{HoKo1964}, and the subsequent work of Kohn and Sham~\cite{KoSh1965}, paved the way for density functional theory, and frameworks based on this approach that aim to treat both exchange and correlation effects accurately. This is the subject of our next section.

\subsection{Density functional theory}

Density functional theory (DFT) has found widespread use in electronic structure calculations and has been applied to a wide variety of materials, including III-N materials~\cite{KiSt2015,BrGe2000,MuWa2001,KiRi2011}. Indeed, it is DFT that underlies the current state-of-the-art approach to determining the electronic structure of the optoelectronically relevant III-N semiconductors, the so-called Heyd-Scuseria-Ernzerhof (HSE) hybrid functional DFT~\cite{CaSc2013local,HeSc2003}. However, in order to discuss the strengths and weaknesses of this approach we must briefly discuss the underlying ideas of DFT and the HSE hybrid functional. To do so, we follow the approach taken in Ref.~\cite{HodmasElectronicBands}.

The foundation of DFT is the proof by Hohenberg and Kohn~\cite{HoKo1964} that the ground state energy of a many particle system is uniquely determined by a functional of the ground state electron density. The primary benefit of this approach is that it avoids the need to determine the total wave function described above in Eq.~\eqref{eq:SlaterDet}, $\Psi_{n}$. From Eq.~\eqref{eq:SlaterDet} we can see that $\Psi_{n}$ depends on the three spatial coordinates of every electron, and as mentioned, there are around $10^{23}$ valence electrons involved in bonding in every cubic centimetre of a solid. Thus, the total wave function of such a system would depend on around $3\times10^{23}$ variable spatial coordinates. Instead, the ground state electron density, $\rho(\mathbf{r})$, depends only on three spatial variables. Proof of the theorem of Hohenberg and Kohn can be found in the original paper in Ref.~\cite{HoKo1964} and is also discussed in Refs.~\cite{HodmasElectronicBands} and~\cite{ModernManyParticlePhysics}. 
Although Hohenberg and Kohn's theorem provided a general framework for the determination of the ground state electron density, $\rho_{\text{GS}}$, and energy, $\text{E}_{\text{GS}}$, in principle, one still must determine the functional form of the Hamiltonian. For instance, how the kinetic energy, $T$, depends on the electron density, $\rho$, in general. 

An important approach developed to overcome this was the method of Kohn and Sham~\cite{KoSh1965}. In Kohn and Sham's method, the interacting many particle system is mapped to an effective single particle system, where the kinetic energy contribution is given by the kinetic energy of the non-interacting single particle system, $T_{0}$, with an effective potential, $V_{\text{eff}}$, where part of $V_{\text{eff}}$ takes into account exchange and correlation effects. The energy functional, $\text{E}[\rho]$, of the system can then be decomposed as follows~\cite{ModernManyParticlePhysics};
\begin{align}
    \text{E}[\rho] &= T_{0}[\rho] + V_{\text{eff}}[\rho] \,\, , \nonumber \\
    &= T_{0}[\rho] + \int d\mathbf{r}\rho(\mathbf{r})\left[V_{\text{ext}}(\mathbf{r}) + U(\mathbf{r})\right] + \text{E}_{\text{xc}}[\rho] \,\, ,
    \label{eq:KohnShamEq}
\end{align}
where $T_{0}[\rho]$ is the kinetic energy of the non-interacting system with density $\rho$, $V_{\text{ext}}$ is the external potential in which the particles move, $U(\mathbf{r})$ is the Hartree contribution to the interaction, which would be the Coulomb potential for electrons;
\begin{equation*}
    U(\mathbf{r}) = \frac{e^{2}}{2}\int d\mathbf{r}\int d\mathbf{r}' \frac{\rho(\mathbf{r})\rho(\mathbf{r}')}{|\mathbf{r}-\mathbf{r}'|}\,\, ,
\end{equation*}
and finally, the term $\text{E}_{\text{xc}}[\rho]$ is the so-called exchange-correlation functional. All of the terms in Eq.~\eqref{eq:KohnShamEq} except $\text{E}_{\text{xc}}$ can be treated exactly, i.e., the form of their functionals can be determined. Furthermore, even though the kinetic energy of the non-interacting system will not generally be equal to the kinetic energy of the interacting system, it is assumed that corrections to the kinetic energy due to the interaction will be accounted for in the, as yet undetermined, exchange-correlation functional. 

One of the first approximations explored for $\text{E}_{\text{xc}}$ was the local density approximation (LDA), in which an inhomogeneous electron density is approximated locally by a homogeneous electron gas with a given density~\cite{KoSh1965,Martin2020ElecStruc}. LDA DFT has been shown to work remarkably well when applied to some materials, while systematically failing for others~\cite{HodmasElectronicBands,Perdew1985,CoMo2012}. One of the most well known issues with LDA is its underestimation of the band gap of many semiconductors and insulators, the so-called band gap problem~\cite{Perdew1985}. 

Beyond LDA, the generalised gradient approximation (GGA) was developed, in which the exchange-correlation functional includes terms that depend on the gradient of the density~\cite{PeBu1996GGA,PeBu1996Simple}. However, even in GGA approaches, the band gap problem persisted. To overcome this, the idea to combine Hartree-Fock theory and LDA/GGA functional-based approaches was introduced by Becke \emph{et al.}~\cite{HodmasElectronicBands,Becke1993HF}. In Hartree-Fock theory, exchange effects can be taken into account exactly, thus Becke \emph{et al.} introduced the idea of ``hybrid functionals'' which contain contributions from the Hartree-Fock exchange functional and an LDA/GGA exchange-correlation functional.

Many different approaches to approximate the exchange-correlation functional are available in the literature~\cite{HodmasElectronicBands}. As mentioned at the beiginning of this section, an approach that has attracted a lot of attention recently is due to Heyd \emph{et al.}, who introduced the Heyd-Scuseria-Ernzerhof (HSE) screened hybrid functional, which is a combination of the Hartree-Fock exchange functional and a so-called PBE0 functional, with a linear mixing parameter $a$~\cite{HeSc2003,PeEr1996}. One of the primary benefits of the HSE hybrid functional is that it makes use of the screening of the long range part of the Coulomb potential to significantly reduce computational cost~\cite{PaMa2006}. 

The HSE hybrid functional has been applied to a wide range of material systems with great success~\cite{WrKu2009,MaHo2008,HoGa2011,WuSt2010}, including III-N semiconductors~\cite{MoVa2010,MoMi2011,CaSc2013local,DeRi2009}. However, the introduction of the mixing parameter, $a$, represents a step away from the \emph{first principles} origins of DFT and can introduce ambiguity in certain situations. For instance, the $a$ parameter that describes certain defect geometries of a system might differ from the $a$ that describes the band gap accurately~\cite{CaAl2009PRB,CaAl2009Physica}. Furthermore, and of particular relevance to this work, is the ambiguity in choosing $a$ parameters for alloys of, for instance, III-N semiconductors. While the $a$ parameters of the constituent binary III-N materials may be well determined by fitting, there is no rigorous prescription for assigning an $a$ parameter to an alloy of these materials, other than semi-empirical fitting~\cite{MoMi2011}.

A further limitation of DFT methods within the Kohn-Sham framework is that the computational complexity of the problem scales with the cube of the number of atoms considered~\cite{DFTAPracticalIntroduction}. Thus, even though the HSE functional represents a significant improvement in terms of computational cost relative to, for instance PBE0 hybrid functional methods~\cite{PaMa2006}, typical HSE simulations are generally on the order of 10s to, at most, 100s of atoms~\cite{MoMi2011,CaAl2009PRB,CaSc2012}. Although implementations of linear-scaling DFT are currently being developed, so far they have only been applied to a small subset of material science problems~\cite{ScPa2019,PrAa2020,NaBa2020}. A very recently developed state-of-the-art linear scaling DFT software, the CONQUEST code~\cite{NaBa2020}, may provide a promising route to realising DFT-based simulations of, for instance, mesoscopic nanostructures. However, any linear-scaling DFT approach would require significant research efforts to be applied to and optimised for, e.g., III-N-based systems. Further work would also be required to accurately describe a semiconductor alloy in such a framework.

As we have discussed above in Sec.~\ref{sec:LEDsIIINMaterials}, the need to include many 1000s of atoms in models of III-N heterostructures arises from the length scales over which electrons and holes localise. Typical localisation lengths for holes in, for instance, InGaN QW systems, are on the order of 1 nm~\cite{WaGo2011}, and so simulation supercells must have larger dimensions to capture these effects. Generally, it is important to balance accuracy with computational cost, and computationally less intense semi-empirical methods have been developed that accurately reproduce important features of bulk band structures. Semi-empirical approaches can be designed to reproduce experimental or DFT results by fitting model input parameters to data, allowing the semi-empirical model to be extended to larger systems. These methods are the topic of the next section.

\section{Semi-empirical methods}
\label{sec:ElecStrucSemiEmp}

Semi-empirical methods have been widely used in electronic structure calculations of solids generally~\cite{MultibandEffMassIntro} and due to their lower computational cost (compared to \emph{ab initio} methods) have found widespread use in the analysis of many optoelectronically relevant materials, nanostructures and devices~\cite{HodmasElectronicBands,MultibandEffMassIntro,AufderMaurReview}. Broadly, these methods can be divided into two categories, continuum-based and atomistic. Continuum-based approaches essentially assume a uniform repeated unit cell describes a material, while atomistic approaches include internal degrees of freedom, i.e., parameters that may vary within unit cells. In this section we will discuss the ideas underlying a class of continuum-based models, so-called \kp theories, widely used in the analysis of not only conventional III-V semiconductors, but also III-Ns, despite most of these continuum-based models not taking carrier localisation effects into account, a key component in modelling III-Ns. We will then discuss recent efforts to include carrier localisation in these models, as well as the limitations of these approaches. Finally, we will turn to the atomistic tight-binding approach upon which the work here is based. 

\subsection{Effective mass approximation and \kp theories}

In our discussion of \emph{ab initio} methods above, we saw that, for instance, DFT approaches can be used to determine the band structure. We also saw that for a given energy band, $\text{E}(\mathbf{k})$, we can restrict our analysis to the first Brillouin zone (FBZ). DFT approaches allow the determination of $\text{E}(\mathbf{k})$ throughout the entire FBZ. However, certain properties of a system may be strongly dependent on a specific region of the FBZ, for example the optical properties of GaAs systems are determined by the behaviour of $\text{E}(\mathbf{k})$ around $\mathbf{k}=\mathbf{0}$ as it is a direct band gap semiconductor~\cite{HodmasElectronicBands}. In principle, the behaviour of any given energy band at a given $\mathbf{k}$-point will depend on all the other bands, however, the influence of bands that are further away (in energy) will be negligible. In a first step, we explicitly analyse only one band, and encompass modifications to this band due to the coupling between it and the other bands via an effective mass approximation (EMA). In such an approach, the energy dispersion relation, $\text{E}(\mathbf{k})$, around $\mathbf{k}=\mathbf{0}$ is approximated to be parabolic~\cite{HodmasElectronicBands};
\begin{equation*}
    \text{E}(\mathbf{k}) = \frac{\hbar^{2}k^{2}}{2m^{*}} \,\, ,
\end{equation*}
where $k=|\mathbf{k}|$, $\hbar$ is the reduced Planck constant and $m^{*}$ is the so-called ``effective'' mass of the carrier. The effective mass is implicitly modified by the other bands and can be determined from, for example, experimental measurements or DFT band structures~\cite{HaugAndKochThird}.

\begin{figure}[t!]
\centering
    \includegraphics[width=\columnwidth]{Figures/TBGaNBandStructureNearGamma.eps}
    \caption{Band structure of wurtzite GaN determined in a tight-binding model around the $\Gamma$-point. Tight-binding parameters from Ref.~\cite{CaSc2013local} have been used. The figure highlights the more complicated valence band structure when compared to the conduction band.}
    \label{fig:TBBandsGaNGamma}
\end{figure}

Single band EMAs are often applied to bands around $\mathbf{k}=\mathbf{0}$, the $\Gamma$-point, both for valence and conduction bands when dealing with direct gap III-V materials and heterostructures. Figure~\ref{fig:TBBandsGaNGamma} shows an illustrative band structure of GaN near the $\Gamma$-point. Compared to the top of the valence band, the bottom of the conduction band is energetically isolated, and so a single band EMA may be a suitable description. However, the more complicated structure of the top of the valence band, consisting of a number of valence ``sub-bands'' energetically close together, requires ideally an approach that takes into account more bands. Thus, to improve the accuracy of the EMA we can take more bands into account and include the coupling between these bands explicitly. Doing so will also provide a better description of energetically higher and lower lying conduction and valence bands, respectively, which are important for processes such as Auger recombination, in contrast to radiative recombination, where the structure of the band edges (i.e. near the band gap) is more important~\cite{HodmasElectronicBands}. To do so, we begin by considering the single particle Hamiltonian of an electron in a periodic solid;
\begin{equation}
    \label{eq:kp_base_hamiltonian}
    \hat{H} = -\frac{\hbar^{2}}{2m}\nabla^{2} + V(\mathbf{r}) \,\, ,
\end{equation}
where $V(\mathbf{r}) = V(\mathbf{r} + \mathbf{R})$ is the periodic potential arising from the ions of the lattice. We will seek to construct the Hamiltonian in a basis that allows us to account for couplings between bands. To do so, we follow the approach laid out in Ref.~\cite{HodmasElectronicBands}.

The basic assumption underlying this multi-band EMA is that we already have knowledge of the energy levels, $\text{E}_{n}(\mathbf{k})$, and periodic Bloch factors, $u_{n,\mathbf{k}}(\mathbf{r})$, at some $\mathbf{k}$-point, which we denote by $\mathbf{k}_{0}$. By Bloch's theorem (Sec.~\ref{sec:ElecStrucBasicConsiderations}), the electronic wave functions, $\Psi_{n,\mathbf{k}}(\mathbf{r})$, can be expressed as a plane wave times this periodic function, $u_{n,\mathbf{k}}(\mathbf{r})$, such that $\Psi_{n,\mathbf{k}}\left(\mathbf{r}\right) = e^{i\mathbf{k}\cdot\mathbf{r}}u_{n,\mathbf{k}}\left(\mathbf{r}\right)$. Thus, the \sch~equation at $\mathbf{k}_{0}$ is
\begin{equation}
    \label{eq:kp_bloch_hamiltonian_kZero}
    \hat{H}\Psi_{n,\mathbf{k}_{0}} = \hat{H}\left(e^{i\mathbf{k}_{0}\cdot\mathbf{r}}u_{n,\mathbf{k}_{0}}\left(\mathbf{r}\right)\right) = \text{E}_{n}(\mathbf{k}_{0})\left(e^{i\mathbf{k}_{0}\cdot\mathbf{r}}u_{n,\mathbf{k}_{0}}\left(\mathbf{r}\right)\right) \,\, .
\end{equation}
We now seek to find $\text{E}_{n}(\mathbf{k})$ for some $\mathbf{k}$-point close to $\mathbf{k}_{0}$, at which point the \sch~equation is
\begin{equation}
    \label{eq:kp_bloch_hamiltonian_kNonzero}
    \hat{H}\Psi_{n,\mathbf{k}} = \hat{H}\left(e^{i\mathbf{k}\cdot\mathbf{r}}u_{n,\mathbf{k}}\left(\mathbf{r}\right)\right) = \text{E}_{n}(\mathbf{k})\left(e^{i\mathbf{k}\cdot\mathbf{r}}u_{n,\mathbf{k}}\left(\mathbf{r}\right)\right) \,\, .
\end{equation}
We introduce a dependence on $\mathbf{k} - \mathbf{k}_{0}$ by rewriting $e^{i\mathbf{k}\cdot\mathbf{r}} = e^{i(\mathbf{k}-\mathbf{k}_{0}+\mathbf{k}_{0})\cdot\mathbf{r}}$ and rearranging;
\begin{equation}
    \label{eq:kpDisplacedHamiltonian}
    \hat{H}e^{i(\mathbf{k}-\mathbf{k}_{0})\cdot\mathbf{r}}\left(e^{i\mathbf{k}_{0}\cdot\mathbf{r}}u_{n,\mathbf{k}}(\mathbf{r})\right) = \text{E}_{n}(\mathbf{k})e^{i(\mathbf{k}-\mathbf{k}_{0})\cdot\mathbf{r}}\left(e^{i\mathbf{k}_{0}\cdot\mathbf{r}}u_{n,\mathbf{k}}(\mathbf{r})\right) \,\, .
\end{equation}
We now multiply both sides of Eq.~\eqref{eq:kpDisplacedHamiltonian} by $e^{-i(\mathbf{k}-\mathbf{k}_{0})\cdot\mathbf{r}}$;
\begin{align*}
    \left[e^{-i(\mathbf{k}-\mathbf{k}_{0})\cdot\mathbf{r}}\hat{H}e^{i(\mathbf{k}-\mathbf{k}_{0})\cdot\mathbf{r}}\right]&\left(e^{i\mathbf{k}_{0}\cdot\mathbf{r}}u_{n,\mathbf{k}}(\mathbf{r})\right) = \\&\left[e^{-i(\mathbf{k}-\mathbf{k}_{0})\cdot\mathbf{r}}\text{E}_{n}(\mathbf{k})e^{i(\mathbf{k}-\mathbf{k}_{0})\cdot\mathbf{r}}\right]\left(e^{i\mathbf{k}_{0}\cdot\mathbf{r}}u_{n,\mathbf{k}}(\mathbf{r})\right) \,\, .
\end{align*}
The resulting modified \sch~equation includes a Hamiltonian that depends on $\mathbf{q} = \mathbf{k}-\mathbf{k}_{0}$. We denote this new Hamiltonian by $\hat{H}_{\mathbf{q}}$, such that
\begin{equation*}
    \hat{H}_{\mathbf{q}} = e^{-i\mathbf{q}\cdot\mathbf{r}}\hat{H}e^{i\mathbf{q}\cdot\mathbf{r}} \,\, .
\end{equation*}
This $\mathbf{q}$-dependent Hamiltonian acts on Bloch wave functions, $\Phi_{n,\mathbf{k}}(\mathbf{r})$, where
\begin{equation*}
    \Phi_{n,\mathbf{k}}(\mathbf{r}) = e^{i\mathbf{k}_{0}\cdot\mathbf{r}}u_{n,\mathbf{k}}(\mathbf{r}) \,\, .
\end{equation*}
With some manipulation to expand the $\nabla^{2}(\exp{(i\mathbf{q}\cdot\mathbf{r})}\Phi_{n,\mathbf{k}}(\mathbf{r}))$ term we can look at how $\hat{H}_{\mathbf{q}}$ acts on $\Phi_{n,\mathbf{k}}(\mathbf{r})$. It can be shown that
\begin{align*}
    \hat{H}_{\mathbf{q}}\Phi_{n,\mathbf{k}}(\mathbf{r}) &= \left[e^{-i\mathbf{q}\cdot\mathbf{r}}\left(-\frac{\hbar^{2}}{2m}\nabla^{2} + V(\mathbf{r})\right)e^{i\mathbf{q}\cdot\mathbf{r}}\right]\Phi_{n,\mathbf{k}}(\mathbf{r}) \,\, , \\
    &= \left[-\frac{\hbar^{2}}{2m}\nabla^{2} + V\left(\mathbf{r}\right) + \frac{\hbar^{2}}{m}\mathbf{q}\cdot\frac{1}{i}\nabla + \frac{\hbar^{2}q^{2}}{2m}\right]\Phi_{n,\mathbf{k}}(\mathbf{r}) \,\, , \\
    &= \left[\hat{H} + \frac{\hbar^{2}}{m}\mathbf{q}\cdot\frac{1}{i}\nabla + \frac{\hbar^{2}q^{2}}{2m}\right]\Phi_{n,\mathbf{k}}(\mathbf{r}) \,\, , \\
    &= \left[\hat{H} + \hat{H}'\right]\Phi_{n,\mathbf{k}}(\mathbf{r}) \,\, ,
\end{align*}
where $q = |\mathbf{q}|$. Using Eq.~\eqref{eq:kp_base_hamiltonian} to identify $\hat{H}$, the above can be written as $\hat{H}$ plus a correction, $\hat{H}'$, that is given by
\begin{equation}
    \label{eq:kpName}
    \hat{H}' = \frac{\hbar^{2}}{m}\mathbf{q}\cdot\frac{1}{i}\nabla + \frac{\hbar^{2}q^{2}}{2m} = \frac{\hbar}{m}\mathbf{q}\cdot\mathbf{p} + \frac{\hbar^{2}q^{2}}{2m} \,\, ,
\end{equation}
where we have used the definition of the momentum operator $\mathbf{p} = -i\hbar\nabla$ to recast $\hat{H}'$ in terms of $\mathbf{q}$ and $\mathbf{p}$. As mentioned above, the optical properties of the direct band gap semiconductors studied here are determined by the behaviour of the band structure around the $\Gamma$-point, thus, this point is often of most interest~\cite{HodmasElectronicBands}. When $\mathbf{k}_{0}=\mathbf{0}$, $\mathbf{q}=\mathbf{k}$, hence the name of this class of theories derives from the resulting \kp term in Eq.~\eqref{eq:kpName}. However, we note that this method can be also be applied to other choices of $\mathbf{k}_{0}$.

Having identified the correction to $\hat{H}$ when studying the band structure at $\mathbf{k}$ with $\mathbf{k} \neq \mathbf{k}_{0}$, we can treat these extra terms as perturbations. To do so, we apply second-order perturbation theory to derive an expression for the variation of $\text{E}_{n}(\mathbf{k})$ near $\mathbf{k}_{0} = \mathbf{0}$, i.e., near the $\Gamma$-point~\cite{QuantumTheoryOfSolidsEoin}. As an illustration, we consider a non-degenerate band, where, for the $n^{\text{th}}$ such band, we have
\begin{equation}
    \label{eq:kp_SOPT_Ek}
    \text{E}_{n}(\mathbf{k}) = \text{E}_{n}(\mathbf{0}) + \frac{\hbar}{m}\mathbf{k}\cdot\mathbf{p}_{nn} + \frac{\hbar^{2}k^{2}}{2m} + \frac{\hbar^{2}}{m^{2}}\sum_{n' \neq n} \frac{|\mathbf{k}\cdot\mathbf{p}_{nn'}|^{2}}{\text{E}_{n}(\mathbf{0})-\text{E}_{n'}(\mathbf{0})} \,\, .
\end{equation}
where
\begin{equation*}
    \mathbf{p}_{nn'} = \bra{u_{n,\mathbf{0}}}\mathbf{p}\ket{u_{n',\mathbf{0}}} \,\, ,
\end{equation*}
and $\mathbf{p}$ is the momentum operator as before. The momentum matrix elements, $\mathbf{p}_{nn'}$, contribute corrections to the energy of a given state near $\Gamma$ due to the presence of other energy states. It can be shown, due to symmetry, that in diamond crystal structures $\mathbf{p}_{nn} = 0$~\cite{HodmasElectronicBands}, and in III-N systems this term is generally small relative to the quadratic term in Eq.~\eqref{eq:kp_SOPT_Ek}~\cite{HodmasElectronicBands,DuBr2000}. It is also possible to expand the wave function $\Phi_{n,\mathbf{k}}(\mathbf{r})$ near $\Gamma$ in terms of the known wave functions, $\Psi_{n,\mathbf{k}_{0}}(\mathbf{r})$, at $\Gamma$ via~\cite{HaugAndKochThird};
\begin{equation*}
    \Phi_{n,\mathbf{k}}(\mathbf{r}) = \Psi_{n,\mathbf{k}_{0}}(\mathbf{r}) + \frac{\hbar}{m}\sum_{n'\neq n}\frac{\mathbf{k}\cdot\mathbf{p}_{nn'}}{\text{E}_{n}(\mathbf{0})-\text{E}_{n'}(\mathbf{0})}\Psi_{n',\mathbf{k}_{0}}(\mathbf{r}) \,\, .
\end{equation*}

In practice, the sum over $n'$ in Eq.~\eqref{eq:kp_SOPT_Ek} will be dominated by contributions from bands close in energy to the target band, limiting the number of bands we need to consider explicitly. L\"{o}wdin provided a method to categorise bands into ``A'' bands, which will be included explicitly, and ``B'' bands, the impact of which will be included implicitly via the effective mass~\cite{Lowdin1951}. Thus, Eq.~\eqref{eq:kp_SOPT_Ek} becomes
\begin{align}
    \label{eq:kpAB}
    \text{E}_{n}(\mathbf{k}) = &\text{E}_{n}(\mathbf{0}) + \frac{\hbar^{2}k^{2}}{2m} \nonumber \\
    &+ \frac{\hbar^{2}}{m^{2}}\sum_{n'\in A \neq n} \frac{|\mathbf{k}\cdot\mathbf{p}_{nn'}|^{2}}{\text{E}_{n}(\mathbf{0})-\text{E}_{n'}(\mathbf{0})} + \frac{\hbar^{2}}{m^{2}}\sum_{n'\in B} \frac{|\mathbf{k}\cdot\mathbf{p}_{nn'}|^{2}}{\text{E}_{n}(\mathbf{0})-\text{E}_{n'}(\mathbf{0})}\,\, .
\end{align}
In general, Eq.~\eqref{eq:kpAB} can be written as~\cite{HaugAndKochThird}
\begin{equation*}
    \text{E}_{n}(\mathbf{k}) = \text{E}_{n}(\mathbf{0}) + \underbrace{\frac{\hbar^{2}}{m^{2}}\sum_{n'\in A \neq n} \frac{|\mathbf{k}\cdot\mathbf{p}_{nn'}|^{2}}{\text{E}_{n}(\mathbf{0})-\text{E}_{n'}(\mathbf{0})}}_{\text{Explicit A Bands}} + \underbrace{\sum_{i\in\{1,2,3\}}\frac{\hbar^{2}k_{i}^{2}}{2}\left(\frac{1}{m_{i}^{*}}\right)}_{\text{Implicit B Bands}} \,\, ,
\end{equation*}
where the sum over $i$ is over the three directions of reciprocal space, $k_{i}$ is the magnitude of $\mathbf{k}$ in the $i$-direction and $m_{i}^{*}$ is the effective mass in that direction. The effective mass in direction $i$, which is modified by the presence of the B bands, is given by
\begin{equation*}
    \frac{1}{m_{i}^{*}} = \frac{1}{m} + \frac{2}{m^{2}}\sum_{n'\in B}\frac{|\mathbf{\hat{k}}_{i}\cdot\mathbf{p}_{nn'}|^{2}}{\text{E}_{n}(\mathbf{0}) - \text{E}_{n'}(\mathbf{0})}\,\, ,
\end{equation*}
where $\mathbf{\hat{k}}_{i}$ is a unit vector in the $i$-direction of reciprocal space.

The theories described above can also be extended to degenerate bands, as would be required, for instance, for the topmost valence bands in the band structure in Fig.~\ref{fig:TBBandsGaNGamma}. This can be achieved within the framework of degenerate perturbation theory~\cite{MultibandEffMassIntro}. Multi-band approaches have also been developed that take into account spin-orbit coupling, such as the models of Luttinger and Kohn~\cite{ChCh1996}, Kane~\cite{HaugAndKochThird} and Pidgeon and Brown~\cite{PiBr1966}; for III-Ns, it turns out that this effect is generally of secondary importance~\cite{QuantumTheoryOfSolidsEoin}. Multi-band \kp theories of increasing levels of complexity have been developed, generally named after how many bands are treated explicitly~\cite{HodmasElectronicBands}.

However, when considering semiconductor alloys, the above analysis should, in principle, break down due to the presence of alloy disorder. Despite this, it has been found that for many conventional III-V semiconductor alloys (e.g. AlGaAs) the alloy disorder represents a weak perturbation to the system and interpolating material parameters can give a reasonable description of the electronic structure~\cite{HarrisonQWWD}. Such an approximation is another example of a virtual crystal approximation (VCA), as mentioned in Sec.~\ref{sec:LEDsIIINMaterials} to determine, for instance, the lattice parameters of an InGaN alloy via a linear interpolation scheme. Conceptually, the problem of solving the \sch~equation in a disordered, aperiodic potential, is converted to a picture in which the electrons move in a fictitious ideal ``virtual'' crystal. The atoms of this crystal then have the interpolated properties of the alloy constituents, and the model above is applied to an ``effective'' periodic unit cell. VCAs form the basis of many \kp descriptions of semiconductor alloys underlying bulk materials and/or heterostructures.

Multi-band \kp approaches have been applied to the analysis of many nanostructures of realistic sizes, beyond the reach of standard DFT~\cite{HodmasElectronicBands,HarrisonQWWD,MultibandEffMassIntro,QuantumTheoryOfSolidsEoin}. Using VCAs, 1D descriptions of QWs can be targeted as the confining potential is still periodic in the growth plane. In this case, if the growth plane is coplanar with the $x$-$y$ plane, then $k_{x}$ and $k_{y}$ can be used as ``good'' quantum numbers and we need only to solve a particle in a box problem along the growth direction. However, as we discussed in Sec.~\ref{sec:LEDsIIINMaterials}, in III-N alloys, alloy disorder leads to a breakdown of the VCA description due to variations in the electronic structure even within a unit cell. As a result, $k_{x}$ and $k_{y}$ are no longer good quantum numbers and standard continuum-based approaches like those discussed above may not provide a good description of III-N alloys~\cite{AufderMaurReview}. Thus, a number of modified continuum-based approaches have been developed.

\subsection{Modified continuum-based approaches}

To develop modified continuum-based models the confining potential is described in a fully 3D picture. To obtain the confining potential one often starts from a 3D ``compositional'' map of, e.g., cation species in the alloy~\cite{WuWa2003,JoTe2017,WaGo2011,YaSh2014,PiLi2017LLT2,WuSh2012}. In this section, as an example, we will consider the case of InGaN alloys and QWs. Having determined information about the In and Ga atom distributions from, for example, experimental data (i.e. atom clustering or random distribution), a digital alloy is generated with either In or Ga atoms at a given position. However, that resolution may be too fine for a continuum-based model. Thus, a map is generated on a larger length scale that gives the average In/Ga concentration in a given region. Often a Gaussian smoothing procedure is used to determine the local In or Ga content.

In order to determine the confining potential, the position dependence of the local energy band edges are inferred from a VCA-like description. At every point in the simulated material, the conduction and valence band edge energies are determined from the variation of the band gap with In content given by V\'{e}gard's law (see Sec.~\ref{sec:LEDsIIINMaterials}):
\begin{equation}
    \text{E}_{\text{g}}(\text{In$_{x}$Ga$_{1-x}$N}) = x\text{E}_{g}(\text{InN}) + (1-x)\text{E}_{\text{g}}(\text{GaN}) + bx(1-x) \,\, ,
    \label{eq:VegMCBA}
\end{equation}
where $x$ is the In composition at a given point and $b$ is the quadratic bowing parameter. Equation~\eqref{eq:VegMCBA} is widely used~\cite{WeNa2020} and can be refined based on the work in Ref.~\cite{CaSc2013local}. Fluctuations in strain and built-in field can also be accounted for based on the local alloy content (see, for example, Ref.~\cite{PiLi2017LLT2}). Similar VCA-like interpolation schemes, if known, can be used for all material parameters (such as relative permittivity or effective mass~\cite{PiLi2017LLT2}). The 3D potential map is then used as input for solving the \sch~equation. Most often, modified continuum-based models use EMAs. However, even when working in a single band EMA, when targeting finer and finer resolutions of the effective confining potential map, as well as when modelling larger regions of the structure (for instance in device simulations~\cite{OdLu2020}), the computational expense can be significant. Much of the computational expense stems from solving the eigenvalue problem associated with the \sch~equation. Recently, a numerical approach to avoid explicitly solving the \sch~equation for random energy landscapes has been developed, the so-called localisation landscape theory.

In general, localisation landscape theory (LLT) is a method that allows the estimation of a local ground state wave function in a disordered potential energy landscape without solving the \sch~equation. Instead, in LLT, an upper bound for the amplitude of the ground state wave function is determined from the solution to a system of linear equations. For more details on LLT, and a derivation of the so-called landscape equation,
\begin{equation}
    \label{eq:landscapeEq}
    \hat{H}\ket{u} = \identity\,\, ,
\end{equation}
see Ref.~\cite{FiMa2012}. The mathematical object $\ket{u}$ can be used to approximate the ground state wave function and by taking the inverse of $\braket{r|u} = u(r)$, an effective confining potential can be found that describes the energy landscape~\cite{FiPi2017LLT1}. It is possible to construct the higher lying, excited states using $\ket{u}$ as well~\cite{FiPi2017LLT1}.

So far, LLT has been applied to a number of topics using a single band EMA, i.e., by using a single band effective mass Hamiltonian, $\hat{H}$, in Eq.~\eqref{eq:landscapeEq}. Very recently, the method has been established for single band tight-binding Hamiltonians, with some theoretical results being reported~\cite{WaZh2021,ArFi2022}. The numerical efficiency of the method is a clear advantage when determining properties of disordered systems that require larger simulation cells, such as the determination of the density of states or absorption edges~\cite{PiLi2017LLT2,ChKe2020}. On top of this, the effective confining potential, $1/u(r)$, can be used as input to larger scale transport calculations~\cite{OdCh2021,LiPi2017LLT3}, thus, device simulations.

As already indicated above, so far LLT has only been developed for single band models and it is not clear how to extend this scheme to multi-band models. For InGaN systems, since the valence band structure is similar between InN and GaN systems (i.e. the highest valence bands in both InN and GaN are mostly $p_{x}$/$p_{y}$-like in orbital character~\cite{RiWi2008}), a single band approach may be a reasonable approximation to describe localised states in a disordered energy landscape. However, for AlGaN structures, the valence band ordering is different (as the highest valence band in AlN is mostly $p_{z}$-like in character~\cite{JaBa2002}) and changes with Al content. Thus, the single band EMA may present a significant shortcoming and fail to describe important physics that is also relevant for device applications, such as the polarisation characteristics of emitted light~\cite{FiSc2022}.

Furthermore, some processes, such as Auger recombination, may require a description of the electronic structure of the full first Brillouin zone, beyond the reach of a single band framework. Thus, modified continuum-based approaches have been developed in which a multi-band \kp theory is used in conjunction with the extracted 3D confining potential map described above~\cite{WuWa2003,JoTe2017,WaGo2011,YaSh2014,WuSh2012}. However, due to the underlying eigenvalue problem, solving the \sch~equation with a multi-band EMA (with the 3D potential map as input) can also be numerically expensive. Furthermore, for both LLT and multi-band approaches, a Gaussian smoothing of the compositional map is used~\cite{JoTe2017,WaGo2011,PiLi2017LLT2}. This smoothing procedure introduces another degree of freedom, which is related to the softening/interpolation scheme, and may ``wash out'' carrier localisation effects by averaging out large variations in, for instance, confining potential.

Overall, in this section we have seen that there are a number of modified continuum-based approaches that can be applied to determine the electronic structure of semiconductors. However, there are trade-offs that must be made between, for instance, single band versus multi-band descriptions or simulation resolution versus numerical expense. All these modified continuum-based approaches have in common that they cannot resolve atomistic features and can even be computationally expensive when trying to use a fine numerical resolution. A model that bypasses this shortcoming is the tight-binding approach, which allows for an atomistic resolution as well as single band or multi-band descriptions. In the next section we discuss tight-binding models in general and the model underpinning this work.

\subsection{Tight-binding}
\label{sec:TheoryTightbinding}

In the tight-binding (TB) approach, we assume that the electrons in a system are tightly bound, or strongly localised, to the ions in the lattice. The wave functions of the extended bulk, or heterostructure, carrier states are then constructed as linear combinations of these strongly localised atomic-like orbitals~\cite{SlKo1954}. Tight-binding models of varying levels of sophistication have been developed and can be categorised by a number of parameters. For instance, in this work we consider a TB model based on four atomic-like orbitals: $s$, $p_{x}$, $p_{y}$ and $p_{z}$~\cite{ScSc2006}. Such a TB model is said to have an $sp^{3}$ basis but examples of models using, for instance, an $sp^{3}d^{5}s^{*}$ basis can be found~\cite{DiPe2020,LuSc2006,LoPe2014}. 

Given the strongly localised nature of the atomic-like orbitals, it is expected that the interaction between orbitals on nearby atoms decays rapidly, and so TB models with differing ``cutoff distances'' have been used, referred to as nearest neighbour, second nearest neighbour, etc.~\cite{JeDo1989,YaNa1995}. These two assumptions, that a small number of atomic-like orbitals can describe much of the band structure, and that the overlap of atomic-like orbitals decreases rapidly with increasing distance between atomic sites, are the main assumptions of the TB method~\cite{StefanSchulzThesis}. In this work we use a nearest neighbour TB model. Although spin-orbit coupling effects can be taken into account in the TB approach (see, e.g., Ref.~\cite{StefanSchulzThesis}), as we have mentioned before, in III-N materials these effects are generally of secondary importance~\cite{QuantumTheoryOfSolidsEoin}. We begin by discussing the fundamental ideas underlying the TB approach before turning to methods used to include strain and polarisation effects with an atomistic resolution.

\subsubsection{Basic considerations}

To discuss the TB method we will first construct the Hamiltonian for the system. To do so we follow the approach of Ref.~\cite{StefanSchulzThesis}. For illustration, we will consider the case of a ideal bulk crystal so that we may apply Bloch's theorem, and subsequently will discuss the modification of the model for alloys. We begin by considering the strongly localised basis states, $\ket{\mathbf{R}_{n},\alpha,\nu}$, describing the atomic-like orbital of orbital character $\nu$, with $\nu\in\{s,p_{x},p_{y},p_{z}\}$, in the unit cell at position $\mathbf{R}_{n}$, in which there may be anions and cations of differing species labeled by $\alpha$, and $n$ runs over the unit cells of the lattice. Assuming, in a first step, an isolated atom located at spatial position $\mathbf{R}_{l}$, the \sch~equation reads;
\begin{equation*}
    \hat{H}_{\text{at}}\ket{\mathbf{R}_{l},\alpha,\nu} = \text{E}_{\text{at}}^{\alpha,\nu}\ket{\mathbf{R}_{l},\alpha,\nu}\,\, ,
\end{equation*}
where the atomic Hamiltonian is
\begin{equation*}
    \hat{H}_{\text{at}} = \frac{\mathbf{\hat{p}}^{2}}{2m} + V_{\text{at}}(\mathbf{R}_{l},\alpha)\,\, .
\end{equation*}
Here, $V_{\text{at}}(\mathbf{R}_{l},\alpha)$ is the atomic potential of the atom $\alpha$ at $\mathbf{R}_{l}$. Due to the presence of the other atoms in the lattice, the electronic structure of such an atom will be modified, with the associated single particle Hamiltonian;
\begin{equation}
    \label{eq:TBSPHamDef}
    \hat{H}_{\text{SP}} = \hat{H}_{\text{at}}(\mathbf{R}_{l},\alpha) + \hat{V}_{\text{lattice}}(\mathbf{R}_{l},\alpha)\,\, ,
\end{equation}
where
\begin{equation*}
    \hat{V}_{\text{lattice}}(\mathbf{R}_{l},\alpha) = \sum_{m\neq l}\sum_{\alpha'} V_{\text{at}}(\mathbf{R}_{m},\alpha')
\end{equation*}
is the lattice periodic potential due to the other atoms in the lattice. For a single particle wave function characterised by wave vector $\mathbf{k}$, $\ket{\mathbf{k}}$, we are required to solve
\begin{equation}
    \label{eq:TBSPHam}
    \hat{H}_{\text{SP}}\ket{\mathbf{k}} = \text{E}(\mathbf{k})\ket{\mathbf{k}} \,\, . 
\end{equation}
The bulk single particle wave functions can then be represented as linear combinations of Bloch-type wave functions constructed from the atomic orbitals:
\begin{align}
    \label{eq:TBSPWF}
    \ket{\mathbf{k}} &= \sum_{\alpha,\nu}c_{\alpha,\nu}(\mathbf{k})\left(\sqrt{\frac{V_{0}}{V}}\sum_{n}e^{i\mathbf{k}\cdot(\mathbf{R}_{n}+\mathbf{\Delta}_{\alpha})}\ket{\mathbf{R}_{n},\alpha,\nu}\right) \,\, , \nonumber \\
    &= \sum_{\alpha,\nu} c_{\alpha,\nu}(\mathbf{k})\ket{\mathbf{k},\alpha,\nu} \,\, ,
\end{align}
where $c_{\alpha,\nu}(\mathbf{k})$ are the expansion coefficients, $V_{0}$ is the volume of a unit cell, $V$ is the volume of the full system, $n$ runs over all unit cells as before and the position of the atom $\alpha$ within unit cell $n$ is given by $\Delta_{\alpha}$. It can be shown that this ansatz for $\ket{\mathbf{k}}$ satisfies Bloch's theorem~\cite{StefanSchulzThesis}, thus we restrict our analysis to $\mathbf{k}$ being in the first Brillouin zone. Our next step will be to determine the expansion coefficients. To do so, we substitute our ansatz for $\ket{\mathbf{k}}$ given in Eq.~\eqref{eq:TBSPWF} into Eq.~\eqref{eq:TBSPHam} and multiply on the left by the bra vector $\bra{\mathbf{k},\alpha',\nu'}$. After some manipulation we are left with the following matrix equation:
\begin{equation}
    \label{eq:TBSchMatrixEquation}
    \sum_{\alpha,\nu}H^{\text{SP}}_{\alpha',\nu';\alpha,\nu}(\mathbf{k})c_{\alpha,\nu}(\mathbf{k}) = \text{E}(\mathbf{k})\sum_{\alpha,\nu}S_{\alpha',\nu';\alpha,\nu}(\mathbf{k})c_{\alpha,\nu}(\mathbf{k}) \,\, ,
\end{equation}
where the Hamiltonian matrix elements, $H_{\alpha',\nu';\alpha,\nu}^{\text{SP}}(\mathbf{k})$, and overlap matrix elements, $S_{\alpha',\nu';\alpha,\nu}(\mathbf{k})$, are given by
\begin{equation*}
    H_{\alpha',\nu';\alpha,\nu}^{\text{SP}}(\mathbf{k}) = \frac{V_{0}}{V}\sum_{m,n}e^{i\mathbf{k}\cdot(\mathbf{R}_{m} + \Delta_{\alpha'} - \mathbf{R}_{n} - \Delta_{\alpha})}\bra{\mathbf{R}_{m},\alpha',\nu'}\hat{H}_{\text{SP}}\ket{\mathbf{R}_{n},\alpha,\nu} \,\, ,
\end{equation*}
and
\begin{equation*}
    S_{\alpha',\nu';\alpha,\nu}(\mathbf{k}) = \frac{V_{0}}{V}\sum_{m,n}e^{i\mathbf{k}\cdot(\mathbf{R}_{m} + \Delta_{\alpha'} - \mathbf{R}_{n} - \Delta_{\alpha})}\braket{\mathbf{R}_{m},\alpha',\nu'|\mathbf{R}_{n},\alpha,\nu} \,\, . 
\end{equation*}
The atomic-like orbitals are not generally mutually orthogonal to each other when calculating the overlap of wave functions localised on different atoms in the lattice. Thus, the condition that
\begin{equation*}
    \braket{\mathbf{R}_{m},\alpha',\nu'|\mathbf{R}_{n},\alpha,\nu} = \delta_{n,m}\delta_{\alpha',\alpha}\delta_{\nu',\nu} \,\, ,
\end{equation*}
for nearest neighbour sites will not necessarily be fulfilled. However, the orthogonalisation procedure of L\"{o}wdin permits the transformation of the TB basis into an orthogonal basis while preserving the localised nature of the basis orbitals and their symmetry~\cite{Lowdin1951}.

Proceeding now under the assumption that we have used a L\"{o}wdin transformation to work with an orthogonal basis, the elements of the overlap matrix, $S$, become
\begin{equation*}
    S_{\alpha',\nu';\alpha,\nu} = \delta_{\mathbf{R}_{m},\mathbf{R}_{n}}\delta_{\alpha',\alpha}\delta_{\nu',\nu} \,\, ,
\end{equation*}
and overall Eq.~\eqref{eq:TBSchMatrixEquation} reduces to
\begin{equation}
    \label{eq:TBSchMatReduced}
    \sum_{\alpha,\nu}H^{\text{SP}}_{\alpha',\nu';\alpha,\nu}(\mathbf{k})c_{\alpha,\nu}(\mathbf{k}) = \text{E}(\mathbf{k})c_{\alpha',\nu'}(\mathbf{k}) \,\, .
\end{equation}
In order to discuss some of the features of the solutions to Eq.~\eqref{eq:TBSchMatReduced}, we consider a simplified situation in which there is only one atomic species and one atomic-like orbital in the basis. Thus, the sum over $\alpha$ and $\nu$ is eliminated in our discussion here:
\begin{equation*}
    H_{\alpha,\nu;\alpha,\nu}^{\text{SP}}(\mathbf{k}) = \bra{\mathbf{k},\alpha,\nu}\hat{H}_{\text{SP}}\ket{\mathbf{k},\alpha,\nu} = \text{E}(\mathbf{k})\,\, .
\end{equation*}
And,
\begin{align}
    \bra{\mathbf{k},\alpha,\nu}\hat{H}_{\text{SP}}\ket{\mathbf{k},\alpha,\nu} &= \frac{V_{0}}{V}\sum_{m,n}e^{i\mathbf{k}\cdot(\mathbf{R}_{m}-\mathbf{R}_{n})}\bra{\mathbf{R}_{m},\alpha,\nu}\hat{H}_{\text{SP}}\ket{\mathbf{R}_{n},\alpha,\nu}\,\, , \nonumber \\
    &= \frac{V_{0}}{V}\sum_{m,n}e^{i\mathbf{k}\cdot(\mathbf{R}_{m}-\mathbf{R}_{n})}I_{mn} \,\, .
    \label{eq:TBHamMatInts}
\end{align}
Due to the localised nature of the basis states, the $I_{mn}$ terms in Eq.~\eqref{eq:TBHamMatInts} become exponentially small for large values of $|\mathbf{R}_{m}-\mathbf{R}_{n}|$. Referring back to our definition for $\hat{H}_{\text{SP}}$ in Eq.~\eqref{eq:TBSPHamDef}, and considering now a specific atomic site, $m$, we can approximate the integrals, $I_{mn}$ by
\begin{align}
    \label{eq:TBMatrixElements}
    I_{mn} &= \delta_{m,n}\left[\bra{\mathbf{R}_{m},\alpha,\nu}\hat{H}_{\text{at}}\ket{\mathbf{R}_{n},\alpha,\nu} + \sum_{l\neq n}\bra{\mathbf{R}_{m},\alpha,\nu}V_{\text{at}}(\mathbf{R}_{l},\alpha)\ket{\mathbf{R}_{n},\alpha,\nu}\right] \nonumber \\
    &+ \delta_{m,n\pm 1}\left[\sum_{l}\bra{\mathbf{R}_{m},\alpha,\nu}V_{\text{at}}(\mathbf{R}_{l},\alpha)\ket{\mathbf{R}_{n},\alpha,\nu}\right] \nonumber \\
    &+ \delta_{m,n\pm 2}\left[\sum_{l}\bra{\mathbf{R}_{m},\alpha,\nu}V_{\text{at}}(\mathbf{R}_{l},\alpha)\ket{\mathbf{R}_{n},\alpha,\nu}\right] + \ldots \\
    &= \delta_{m,n}\widetilde{\text{E}} + \delta_{m,n\pm 1}\lambda + \ldots \nonumber \,\, .
\end{align}
In our nearest neighbour model we truncate the expression after the $\delta_{m,n\pm 1}$ term. The term $\widetilde{\text{E}}$ is called the on-site energy, and can be interpreted as the atomic energy modified by the presence of the other atoms in the lattice. The $\lambda$ term involves integrals of orbitals located at different atomic sites and is called a hopping integral/matrix element, as it gives the probability amplitude of an electron moving from one site to the other~\cite{StefanSchulzThesis}. 

The contributions to the sums in Eq.~\eqref{eq:TBMatrixElements} involve so-called ``two centre'' and ``three centre'' integrals depending on $\mathbf{R}_{m}$, $\mathbf{R}_{n}$ and $\mathbf{R}_{l}$. A two centre integral involves $\mathbf{R}_{m} = \mathbf{R}_{l} \neq \mathbf{R}_{n}$. A three centre integral involves $\mathbf{R}_{m} \neq \mathbf{R}_{l} \neq \mathbf{R}_{n}$. The two centre approximation of Slater and Koster proposes to discard three centre integrals as they are relatively smaller than the two centre integrals~\cite{StefanSchulzThesis}. Slater and Koster also evaluated expressions for the hopping matrix elements depending on the orbital and bond type, while taking into account the directional cosines between the atomic sites under consideration~\cite{SlKo1954}. 

The remaining question is now how to determine the quantities $\widetilde{\text{E}}$ and $\lambda$, which we call the TB matrix elements. Generally, there are two approaches to determining the TB matrix elements in Eq.~\eqref{eq:TBMatrixElements}. One approach is to assume an explicit form for the basis states and calculate the matrix elements, for instance using \emph{ab initio} methods. This leads to density functional tight-binding (DFTB) approaches~\cite{StefanSchulzThesis}. Alternatively, the TB matrix elements can be treated as free parameters and fitted to band structure data, either experimental or determined via \emph{ab initio} methods. The latter approach is often referred to as empirical or semi-empirical TB. In this work we use such a semi-empirical TB method that has been fitted to HSE hybrid functional DFT band structures of bulk InN, GaN and AlN and has also been benchmarked against experimental data for alloyed systems and heterostructures~\cite{CaSc2013local,TaCa2016}. \blue{The HSE functional HSE06 was utilised for the fitting with an exact exchange parameter $\alpha = 0.25$ and screening parameter $\mu = 0.2$. The cutoff energy for plane waves was 600 eV and a $6\times6\times4$ $\Gamma$-centred Monkhorst $\Vec{k}$-point mesh was used. Further details can be found in Ref.~\cite{CaSc2013local}.}

Throughout this chapter we have stressed the need for an atomistic resolution when modelling an alloy, and above we have discussed the application of a TB model to determine the band structures of III-N bulk materials. However, extending the model to wurtzite III-N alloys is reasonably intuitive. For cation sites (e.g. In, Ga and Al) the nearest neighbour environment always consists of four N anions, in the absence of defects. Thus, in principle (i.e. in the absence of strain effects), assigning a value to the on-site matrix element for cation sites is straightforward, and the bulk value would be appropriate. The InN/GaN band offset can be taken into account by shifting the InN on-site matrix element by the valence band offset, $\Delta \text{E}_{\text{VB}} = 0.62$ eV, as determined via HSE DFT calculations~\cite{MoMi2011}. For the anion sites, however, the nearest neighbour environment can have varying numbers of, e.g., In and Ga atoms. Thus, we must take into account the difference in the on-site matrix element for N in InN compared to GaN~\cite{CaSc2013local}. A common approach is to take a weighted average of the bulk values, weighted in the ratio of nearest neighbour atomic species~\cite{OrLi2002,LiPo1992,BoKh2007}.  

To determine the hopping matrix elements, in the absence of strain effects, the relevant bulk value would be appropriate. However, strain effects will modify the hopping matrix elements due to variations in bond lengths and bond angles. Furthermore, strain effects will also modify the on-site energy due to variations in the local potential due to the surrounding lattice. In the next section we first discuss how we will include these strain effects in the TB model described above, before turning to look at how we determine the strain state of the system. Lastly, we will see how we factor the impact of the built-in polarisation field into our TB model.

\subsubsection{Including strain effects}

In principle, incorporating strain effects in a TB model would require some knowledge of how both the on-site and hopping matrix elements change with variations of the nearest neighbour bond lengths and bond angles. TB models have been developed with this approach but require the fitting of a large number of parameters~\cite{BoKl2002,NiRi2009}. However, it has been shown, at least for an $sp^{3}$ basis, that strain can be included as an on-site correction to the TB Hamiltonian such that~\cite{BoKl2002,JaSc1998}
\begin{equation*}
    \hat{H} = \hat{H}_{0} + \hat{H}_{\text{strain}} \,\, .
\end{equation*}
$\hat{H}_{0}$ is the Hamiltonian of the unstrained system and $\hat{H}_{\text{strain}}$ is a site diagonal Pikus-Bir Hamiltonian for the atomic site $i$ such that~\cite{CaSc2013local,WiSc2006,ScBa2010}
\begin{equation}
    \label{eq:OnSiteStrainCorr}
    \hat{H}_{\text{strain}}^{i} = 
    \begin{pmatrix}
        S_{s} & 0 & 0 & 0 \\
        0 & S_{xx} & S_{xy} & S_{xz} \\
        0 & S_{xy} & S_{yy} & S_{yz} \\
        0 & S_{xz} & S_{yz} & S_{zz}
    \end{pmatrix} \,\, .
\end{equation}
In Eq.~\eqref{eq:OnSiteStrainCorr},
{\allowdisplaybreaks
\begin{align*}
    &S_{s} = a_{ct}(\varepsilon_{11} + \varepsilon_{22}) + a_{cp}\varepsilon_{33} \,\, , \\
    &S_{xx} = (D_{2}+D_{4})(\varepsilon_{11}+\varepsilon_{22}) + D_{5}(\varepsilon_{11}-\varepsilon_{22}) + (D_{1} + D_{3})\varepsilon_{33} \,\, , \\
    &S_{yy} = (D_{2}+D_{4})(\varepsilon_{11}+\varepsilon_{22}) - D_{5}(\varepsilon_{11}-\varepsilon_{22}) + (D_{2} + D_{3})\varepsilon_{33} \,\, , \\
    &S_{zz} = D_{2}(\varepsilon_{11}+\varepsilon_{22}) + D_{1}\varepsilon_{33} \,\, , \\
    &S_{xy} = 2D_{5}\varepsilon_{12} \,\, , \\
    &S_{xz} = \sqrt{2}D_{6}\varepsilon_{13} \,\, , \\
    &S_{yz} = \sqrt{2}D_{6}\varepsilon_{23} \,\, ,
\end{align*}}where the $D_{i}$ terms are the valence band deformation potentials, extracted from HSE DFT~\cite{YaRi2009} and $a_{cp}$ and $a_{ct}$ are the conduction band deformation potentials, related to the band gap deformation potentials, $a_{1}$ and $a_{2}$, given in Ref.~\cite{VuMe2003}, by $a_{1} = a_{cp} - D_{1}$ and $a_{2} = a_{ct} - D_{2}$. The term $S_{s}$ would give, for example, the change in energy of the conduction band edge state with strain. As above, for anion sites, a weighted average in the nearest neighbour atomic species ratio is used. Using this approach, the deformation potentials for the highest valence band and lowest conduction band states are included directly from the HSE DFT calculations, avoiding further fitting of parameters~\cite{CaSc2013local}. 

We note that the deformation potentials are fitted in the vicinity of the $\Gamma$-point~\cite{YaRi2009}, and so may not accurately reproduce the variation of the band structure with strain away from $\Gamma$. However, alloy disorder breaks the translational symmetry of the lattice, meaning we can no longer label quantum states by $\mathbf{k}$, and it is assumed in our analysis below that $\mathbf{k}=\mathbf{0}$ anyway. Moreover, even at $\Gamma$, the energy shift of higher/lower lying states may not be correctly described by these band edge deformation potentials. Regardless, this approach is used in, for instance, 8-band \kp models~\cite{WiSc2006} in the literature and, furthermore, although knowledge of the higher lying states is required to calculate Auger rates, overall, they are of secondary important.

To implement this approach we require some information about the local strain field, which will determine the on-site correction. As discussed above in Sec.~\ref{sec:LEDsIIINMaterials}, strain arises in a material when there is a deformation of the crystal. From an atomistic perspective, local strain fluctuations are caused by the displacement of atoms within a given unit cell as well as overall deformations of the unit cells in the crystal. Thus, we must determine how the atoms in the lattice will be displaced in order to minimise the elastic energy of the system. 

To evaluate the elastic energy of the system we require some knowledge of the interactions between atoms. A number of methods have been developed to model the interatomic interactions~\cite{StWe1985,JoCh1924,MuPo1962,Keating1966}. Here we use the approach developed by Martin~\cite{Martin1970}, which is an example of a valence force field (VFF) model and is adapted from the work in Ref.~\cite{Keating1966}. Further details are given in Ref.~\cite{TaCa2019}. In this VFF model, a specific form is chosen for the interatomic interaction potential that takes into account changes in energy due to changes in bond geometries (discussed in detail below). Overall, when assuming this form for the interatomic potential, the elastic/potential energy at a given atomic site, $V_{i}$, is given by
\begin{flalign}
    \label{eq:VFFMartin}
    %\begin{split}
    V_{i} &= \frac{1}{2}\sum_{j\neq i}\overbrace{\frac{1}{2}k_{r}\Delta r_{ij}^{2}}^{\text{bond stretching}}&&\nonumber\\
    &+ \sum_{j\neq i}\sum_{k\neq i, k>j}\overbrace{\left\{\frac{1}{2}k_{\theta}^{i}r_{ij}^{0}r_{ik}^{0}\Delta\theta_{ijk}^{2}\right.}^{\text{bond bending}}+\overbrace{k_{r\theta}^{i}\Delta\theta_{ijk}\left(r_{ij}^{0}\Delta r_{ij}+r_{ik}^{0}\Delta r_{ik}\right)}^{\text{bond-angle}}+ \overbrace{\left.k_{rr}^{i}\Delta r_{ij}\Delta r_{ik}\right\}}^{\text{bond-bond}}&&\nonumber\\
    &+ \sum_{j'\neq i}\underbrace{\frac{e^{2}}{4\pi\epsilon_{0}\epsilon_{r}}\frac{Z_{i}^{*}Z_{j'}^{*}}{r_{ij'}}}_{\text{Attractive Coulomb}}- \frac{1}{2}\sum_{j\neq i}\underbrace{\frac{1}{4}\alpha_{M}\frac{e^{2}}{4\pi\epsilon_{0}\epsilon_{r}}\frac{Z_{i}^{*}Z_{j}^{*}}{r_{ij}^{2}}\Delta r_{ij}}_{\text{Repulsive Coulomb screening}}\,\, .&&
    %\end{split}
\end{flalign}
In this expression $k_{r}$, $k_{\theta}^{i}$, $k_{r\theta}^{i}$ and $k_{rr}^{i}$ are the force constants related to changes in bond length (bond stretching term), changes in bond angle (bond bending term), the coupling between changes in bond length and changes in bond angle (bond-angle term) and the coupling between changes in bond lengths (of different bonds, bond-bond term), respectively. In the unstrained bulk system, the bond length between atomic site $i$ and $j$ is given by $r_{ij}^{0}$, while the angle between the bonds connecting atomic site $i$ to atomic sites $j$ and $k$ is given by $\theta_{ijk}^{0}$. The bond lengths and bond angles in the strained system are given by $r_{ij}$ and $\theta_{ijk}$, respectively, with the associated variations $\Delta r_{ij} = r_{ij}-r_{ij}^{0}$ and $\Delta\theta_{ijk} = \theta_{ijk}-\theta_{ijk}^{0}$ (see Fig.~\ref{fig:VFFTerms} below). The elementary charge is denoted by $e$, the permittivity of vacuum is $\epsilon_{0}$ and the relative permittivity is $\epsilon_{r}$, $Z_{i}^{*}$ is the effective charge for atom $i$, and finally, $\alpha_{M}$ is the Madelung constant. We will discuss each individually below.

\begin{figure}[t!]
\centering
    \begin{tabular}{|c|c|}
    \hline
    \includegraphics[width=0.45\columnwidth]{Figures/VFFBondStretching.eps} & \includegraphics[width=0.45\columnwidth]{Figures/VFFBondBending.eps} \\
    \hline
    \includegraphics[width=0.45\columnwidth]{Figures/VFFBondAngle.eps} & \includegraphics[width=0.45\columnwidth]{Figures/VFFBondBond.eps} \\
    \hline
    \end{tabular}
    \caption{Schematics of the origins of the a) bond stretching, b) bond bending, c) bond-angle and d) bond-bond terms in Eq.~\eqref{eq:VFFMartin}. In all four diagrams, the top configuration represents the unstrained, ideal bulk situation, while the bottom is an example of a modification to bond length or bond angle that may occur in an alloy situation, wherein the system will be strained. Atomic sites are labeled by $i$, $j$ and $k$. Equilibrium (i.e. unstrained bulk) bond lengths and bond angles are denoted by $r_{ij}^{0}$ and $\theta_{ijk}^{0}$, respectively. Alloy (strained) bond lengths and bond angles are denoted by $r_{ij}$ and $\theta_{ijk}$, respectively.}
    \label{fig:VFFTerms}
\end{figure}

Figure~\ref{fig:VFFTerms} shows schematics of the origin of the bond stretching, bond bending, bond-angle and bond-bond terms. The bond stretching term in Eq.~\eqref{eq:VFFMartin}, due to the situation depicted in Fig.~\ref{fig:VFFTerms} a), is reasonably straightforward, and tracks the change in energy due to changes in bond lengths. The bond bending term, seen in Fig.~\ref{fig:VFFTerms} b), accounts for changes in energy due to changes in the angle between bonds. In Fig.~\ref{fig:VFFTerms} c) we see that changes in bond length may induce changes in bond angle, thus the bond-angle term takes these contributions into account. Lastly, Fig.~\ref{fig:VFFTerms} d) depicts a situation in which the change in length of one bond induces a change in length of a different bond, which is included via the bond-bond term in Eq.~\eqref{eq:VFFMartin}. In each of the situations depicted in Fig.~\ref{fig:VFFTerms}, only one possible deformation is shown in the strained case as an example. For instance in Fig.~\ref{fig:VFFTerms} a), $r_{ij} > r_{ij}^{0}$, but depending on the local environment, it could also be the case that $r_{ij} < r_{ij}^{0}$ in the strained system.

Each sum in Eq.~\eqref{eq:VFFMartin}, except the sums involving Coulomb terms, runs over the four nearest neighbour bonds and six nearest neighbour angles in a local tetrahedron, in the case of the wurtzite structures considered here, illustrated in Fig.~\ref{fig:LocalTet}. In the attractive Coulomb term, the sum over $j$ runs over all atoms in the crystal, thus taking into account the potential due to the rest of the lattice. The repulsive Coulomb screening term in Eq.~\eqref{eq:VFFMartin} is included to screen the Coulomb interaction in the nearest neighbour tetrahedron, a step which is required to preserve the stability of the crystal~\cite{Martin1970}. The force constants and effective charges are fitted to HSE DFT results to accurately reproduce elastic constants and internal strain in bulk InN, GaN and AlN. The model also reproduces the non-ideal $c/a$ ratio of wurtzite nitrides, as well as the internal parameter $u$. For more details see Refs.~\cite{DanielTannerThesis} and~\cite{ScCa2015}. The model has been implemented in the molecular dynamics software package LAMMPS, in which the potential energy of the structure is minimised in order to determine the relaxed atomic positions~\cite{ThAk2022}. In the implementation of an alloy in this framework, when a mix of, for instance, InN and GaN bonds are involved, weighted averages of the bulk material force constants are taken, otherwise the bulk parameters are used.

\begin{figure}[t!]
\centering
    \includegraphics[width=0.5\columnwidth]{Figures/figLocalTet.eps}
    \caption{Local tetrahedron of the cation (blue filled circle) atomic site located at position $\mathbf{r}_{0}$ formed by nearest neighbours located at vectors $\mathbf{d}^{\alpha}$ from atom 0. Only the vector to atom 1, $\mathbf{d}^{1}$, is shown. Bonds are depicted by solid lines. $\mathbf{d}^{1}$ is colinear with the bond between atom 0 and atom 1. The edges of the local tetrahedron, $\mathbf{R}_{ij}$, are depicted by dashed arrows. }
    \label{fig:LocalTet}
\end{figure}

Finally, with the relaxed atomic positions determined, a local strain tensor can be constructed using the method in Ref.~\cite{PrKi1998}, see also Ref.~\cite{CaSc2013local}. The local strain tensor at $\mathbf{r}$, $\varepsilon_{ij}(\mathbf{r})$, is given by
\begin{equation}
    \begin{pmatrix}
        \varepsilon_{xx} & \varepsilon_{yx} & \varepsilon_{zx} \\
        \varepsilon_{xy} & \varepsilon_{yy} & \varepsilon_{zy} \\
        \varepsilon_{xz} & \varepsilon_{yz} & \varepsilon_{zz} 
    \end{pmatrix} = 
    \begin{pmatrix}
        R_{12,x}^{0} & R_{23,x}^{0} & R_{34,x}^{0} \\
        R_{12,y}^{0} & R_{23,y}^{0} & R_{34,y}^{0} \\
        R_{12,z}^{0} & R_{23,z}^{0} & R_{34,z}^{0} 
    \end{pmatrix}^{-1}
    \begin{pmatrix}
        R_{12,x} & R_{23,x} & R_{34,x} \\
        R_{12,y} & R_{23,y} & R_{34,y} \\
        R_{12,z} & R_{23,z} & R_{34,z}
    \end{pmatrix}
    -\identity \,\, ,
    \label{eq:strainTensor}
\end{equation}
where, for example, $\mathbf{R}_{12} = (R_{12,x},R_{12,y},R_{12,z})^{T}$ is the edge of the local tetrahedron between atoms 1 and 2 in Fig.~\ref{fig:LocalTet}. With the local strain tensor determined, we can include the effect of strain on the TB Hamiltonian determined above via Eq.~\eqref{eq:OnSiteStrainCorr}. However, as discussed above in Sec.~\ref{sec:LEDsIIINMaterials}, strain effects will also lead to piezoelectric polarisation effects, and we seek a method to treat local fluctuations in piezoelectric polarisation that arise from local fluctuations in strain. We now turn to discuss the method developed to achieve this.

\subsubsection{Including polarisation field effects}

As discussed above in Sec.~\ref{sec:LEDsIIINMaterials}, strain in wurtzite III-N QWs lead to piezoelectric polarisation effects. Indeed, in an InGaN alloy the difference in lattice constants between InN and GaN is $\sim$ 10\%, thus we expect relatively large local fluctuations may occur depending on the local atomic environment. With the model that can describe local fluctuations in strain established above, we are in a position to determine the associated local fluctuations in piezoelectric polarisation via the approach detailed in Ref.~\cite{CaSc2013local}. The central result of that work is that the $i^{\text{th}}$ component of the total polarisation at an atomic site can be separated into a macroscopic and local contribution, given by, for the atomic site labeled $0$;
\begin{equation}
    \label{eq:LocalPol}
    P_{i} = \underbrace{\sum_{j=1}^{6}e_{ij}^{(0)}\varepsilon_{j}}_{\text{macroscopic}} + \underbrace{P_{i}^{\text{SP}}-\frac{e}{V_{0}}\frac{Z_{i}^{0}}{N_{\text{coor}}^{0}}\left[\mu_{i}-\sum_{j=1}^{3}\left(\delta_{ij}+\varepsilon_{ij}\right)\mu_{j,0}\right]}_{\text{local}} \,\, .
\end{equation}
In this expression $e_{ij}^{(0)}$ are the piezoelectric coefficients determined by a ``clamped-ion'' calculation in which the internal degrees of freedom are not allowed to relax~\cite{BeFi1997}. The macroscopic strain is denoted by $\varepsilon_{j}$ in Voigt notation. The $i^{\text{th}}$ component of the spontaneous polarisation is denoted by $P_{i}^{\text{SP}}$. The volume of the local tetrahedron is $V_{0}$, $Z_{i}^{0}$ is the Born effective charge of atom 0, $N_{\text{coor}}^{0}$ is the number of nearest neighbours of atom 0. The $i^{\text{th}}$ component of the bond asymmetry parameter, $\mu_{i}$, is given by 
\begin{equation*}
    \mu_{i} = \sum_{\alpha=1}^{N_{\text{coor}}^{0}} d_{i}^{\alpha} \,\, ,
\end{equation*}
where
\begin{equation*}
    d_{i}^{\alpha} = \sum_{j=1}^{3}\left(\delta_{ij} + \varepsilon_{ij}\right)d_{j,0}^{\alpha} + t_{i}^{\alpha} - t_{i}^{0} \,\, ,
\end{equation*}
and $d_{j,0}^{\alpha}$ is the $j^{\text{th}}$ component of the vector from the atom under consideration, atom 0, to each of the nearest neighbours for the unstrained system (see Fig.~\ref{fig:LocalTet}). The $i^{\text{th}}$ component of the internal strain vector for atom $\alpha$ is $t_{i}^{\alpha}$, hence $t_{i}^{0}$ corresponds to atom 0. Summing these vectors over all the nearest neighbours gives a measure of the bond asymmetry, $\boldsymbol{\mu}$, from which the piezoelectricity in strained systems arises. Determining this vector for an unstrained binary wurtzite system, $\boldsymbol{\mu}_{0}$, we would see that only the $z$ component is nonzero, while for an unstrained zinc blende binary structure, every component is zero~\cite{CaSc2012AlloyFlucsNanostrucs}. The components of $\boldsymbol{\mu}_{0}$ contribute to the local part of Eq.~\eqref{eq:LocalPol} via $\mu_{j,0}$.

A number of approximations have been employed in the derivation of this equation in Ref.~\cite{CaSc2013local}. These include 1) truncating the description of the piezoelectric response to strain to first order in both the macroscopic and internal strain, 2) assuming the spontaneous polarisation is constant throughout the crystal for binaries, 3) that the off-diagonal components of the Born effective charge tensor are zero, 4) using a spherical approximation for the Born effective charge of the nearest neighbours of the atom under consideration and 5) neglecting the strain dependence of the Born effective charge. The impact of these approximations on the results is discussed at length in Ref.~\cite{CaSc2013local}.

Having determined the polarisation vector, we must now determine the associated potential, as described above in Sec.~\ref{sec:LEDsIIINMaterials}. Previously, we used Poisson's equation (Eq.~\eqref{eq:Poisson}) to determine the potential in a continuum-based approach. In the atomistic approach here, however, the relaxed atomic positions make an irregular grid, presenting a further technical challenge. Instead, a point-dipole method is used, in which the potential at each atomic site is determined as the sum of all the dipoles in the lattice~\cite{CaSc2013local}.

Generally, polarisation is defined to be dipole moment per unit of volume, and so here, the $P_{i}$ in Eq.~\eqref{eq:LocalPol} can be considered a density of dipole moments. A natural discretisation of the strained, irregular grid determined by the relaxed atomic positions is the irregular tetrahedron formed by the nearest neighbours of a given atomic site. This volume, determined by the atomic positions $\mathbf{r}_{j}$, $j\in\{1,2,3,4\}$ as in Fig.~\ref{fig:LocalTet}, can be determined to be
\begin{equation*}
    V_{1234} = \frac{1}{6}|(\mathbf{r}_{1}-\mathbf{r}_{4})\cdot\left[(\mathbf{r}_{2}-\mathbf{r}_{4})\times(\mathbf{r}_{3}-\mathbf{r}_{4})\right]| \,\, .
\end{equation*}
Labeling the grid points by their positions $\mathbf{R}_{i}$ (and ensuring that the volume of the full simulation cell is considered~\cite{CaSc2013local}), the dipole moment, $\mathbf{p}_{i}$, in the corresponding volume, $V_{i}$, is given by
\begin{equation*}
    \mathbf{p}_{i} = \mathbf{P}_{i}V_{i} \,\, .
\end{equation*}
For instance, the point dipole at $\mathbf{R}_{i}$ will generate an $\mathbf{r}$-dependent potential, $\phi_{i}(\mathbf{r})$, given by the multipole expansion~\cite{Zangwill};
\begin{equation*}
    \phi_{i}(\mathbf{r}) = \frac{1}{4\pi\epsilon_{0}\epsilon_{r}}\frac{\mathbf{p}_{i}\cdot\left(\mathbf{r}-\mathbf{R}_{i}\right)}{|\mathbf{r}-\mathbf{R}_{i}|^{3}} \,\, .
\end{equation*}
Thus, the total potential at a given atomic site will be the sum of all the contributions at that point. For example, the polarisation potential, $\varphi_{\text{pol}}$, at the atomic site at position $\mathbf{r}_{i}$, will be given by
\begin{equation}
    \label{eq:PolAtAtomicSite}
    \varphi_{\text{pol}}(\mathbf{r}_{i}) = \sum_{j\neq i}\phi_{j}(\mathbf{r}_{i}) \,\, .
\end{equation}

To include the polarisation in the calculation of the electronic structure in our TB framework, we use a site-diagonal correction, similar to the approach for strain corrections described above. For the atomic site at $\mathbf{r}_{i}$, the correction to the Hamiltonian is
\begin{equation*}
    \hat{H}_{\text{pol}}(\mathbf{r}_{i}) = 
    \begin{pmatrix}
        \varphi_{\text{pol}}(\mathbf{r}_{i}) & 0 & 0 & 0 \\
        0 & \varphi_{\text{pol}}(\mathbf{r}_{i}) & 0 & 0 \\
        0 & 0 & \varphi_{\text{pol}}(\mathbf{r}_{i}) & 0 \\
        0 & 0 & 0 & \varphi_{\text{pol}}(\mathbf{r}_{i})
    \end{pmatrix}
    \,\, .
\end{equation*}
Finally, the full Hamiltonian is 
\begin{equation*}
    \hat{H} = \hat{H}_{\text{SP}} + \hat{H}_{\text{strain}} - q\hat{H}_{\text{pol}} \,\, .
\end{equation*}

\blue{In summary, the following is an overview of the method used to set up the TB Hamiltonian for a strained, alloy disordered III-N alloy or heterostructure: 
\begin{enumerate}
    \item TB parameters are determined by fitting to HSE-DFT band structures (DFT implemented in VASP, see Ref.~\cite{CaSc2013local} for further details).
    \item Large scale disorderd supercells set up with random seeding of In and Ga atoms in accordance with experimental data and usual practice in semi-empirical calculations (implemented in custom code).
    \item Equilibrium atomic positions are determined (via atomistic VFF) by minimising the elastic energy for each supercell generated in Step 2 above. The VFF model is implemented in the numerically efficient and massively parallelised software package LAMMPS~\cite{ThAk2022}.
    \item The local strain tensor of the relaxed supercell is determined on the basis of Eq.~\eqref{eq:strainTensor} (implemented in custom code).
    \item The built-in polarisation field is determined via the procedure outlined in Sec.~\ref{sec:ElecStrucSemiEmp}.3 (implemented in custom code).
    \item TB Hamiltonian is constructed including strain and polarisation fields as outlined above (implemented in custom code).
    \item Diagonalising this TB Hamiltonian returns the energy eigenvalues and eigenstates of the disordered system (a Jacobi-Davidson method, JDQR, is used to calculate a selected number of eigenstates and eigenvalues of large sparse matrices in a numerically efficient manner~\cite{SlVa2000}).
\end{enumerate}}

The electronic structure model established allows us to determine the energies and wave functions that are required to calculate the radiative and Auger recombination rates. Gaining insight into these rates will ultimately enable us to target questions about the IQE and the influence of carrier localisation. We proceed in the next chapter to describe a framework for achieving this, as well as discussing the screening of the polarisation field as carrier density increases.

\chapter{Calculating radiative and Auger recombination rates in a tight-binding framework}
\label{chp:RecombinationRates}

In the preceeding chapter we established that determining the electronic structure to an atomistic resolution while including strain, polarisation and carrier localisation effects is crucial for the analysis of $c$-plane III-N-based QWs. We now proceed within the TB framework detailed above, seeking a way to determine the radiative and Auger recombination rates entering the ABC model discussed above in Sec.~\ref{sec:LEDsModellingEfficiency}. Doing so will allow us to gain insight into the competition between these processes, and subsequently their impact on the IQE. The framework described below can be applied generally, but we will highlight where the approach is specific to III-N systems or the TB model employed here. We will first discuss the determination of the rate of radiative recombination, before turning to the rate of Auger recombination, and finally, modifications to the framework to take into account carrier density dependent built-in field screening effects in III-N QWs. 

\section{Calculation of the radiative recombination rate}
\label{sec:RecombinationRadiative}

To gain insight into the IQE, as determined in the ABC model, we must determine the coefficient of spontaneous radiative recombination, $B$. The $B$ coefficient is related to the rate of spontaneous radiative recombination per unit volume per time, $r_{\text{sp}}$ via $B=r_{\text{sp}}/(np)$, where $n$ is the electron and $p$ is the hole carrier density. We begin by deriving an expression for the total spontaneous radiative recombination rate in the active layer, $R_{\text{sp}}$. We mainly follow Ref.~\cite{RidleyQuantumProcesses}, but have modified the approach to connect it to our TB framework. To begin, we consider the transition rate between two discrete energy levels, $\text{E}_{i}$ for an electron state in the conduction band and $\text{E}_{f}$ for a hole state in the valence band. 

We apply Fermi's Golden Rule to derive an expression for the transition probability per unit time (rate), $W_{i\rightarrow f}$, between such states~\cite{HarrisonQWWD,RidleyQuantumProcesses,GrAsIntroToQuantumOptics2010}:
\begin{equation}
    \label{eq:FermisGoldernRule}
    W_{i\rightarrow f} = \frac{2\pi}{\hbar}|\bra{i}\hat{H}\ket{f}|^{2}\rho(\text{E})\delta(\hbar\omega - \text{E}_{i} - \text{E}_{f}) \,\, ,
\end{equation}
where $\rho(\text{E})$ is the density of photon states of energy $\text{E}=\hbar\omega$ per unit volume per unit energy, the $\delta$-function ensures that the emitted photon obeys energy conservation and the modulus square of the interaction matrix element, $|\bra{i}\hat{H}\ket{f}|^{2} = |H_{if}|^{2}$, is given by
\begin{equation}
    \label{eq:RadIntMatEl}
    |H_{if}|^{2} = \frac{e^{2}}{m^{2}}\frac{2\hbar}{\epsilon_{0}n_{r}^{2}\omega}\frac{1}{4}|\bra{i}\mathbf{a}\cdot\mathbf{p}\ket{f}|^{2} \,\, .
\end{equation}
In Eq.~\eqref{eq:RadIntMatEl}, $e$ is the electron charge, $m$ is the electron mass, $\epsilon_{0}$ is the vacuum permittivity, $n_{r}$ is the refractive index of the material, $\omega$ is the frequency of the photon, $\mathbf{a}$ is a unit vector defining the light polarisation, $\mathbf{p}$ is a vector of the momentum operator for each direction, $\mathbf{k}$ is the photon wave vector, $\mathbf{r}$ is a position vector and we have made use of the dipole approximation~\cite{RidleyQuantumProcesses}. The dipole approximation assumes that the extent of the carrier wave function is much smaller than the wavelength of emitted light, thus $\mathbf{k}\cdot\mathbf{r} \ll 1$ and so $e^{i\mathbf{k}\cdot\mathbf{r}} \approx 1$. It is a commonly used approximation given the wavelength of UV (100 -- 280 nm) and visible (280 -- 700 nm) light in comparison to the extent of carriers ($\sim$ 1 nm). The quantity $\rho(\text{E})$ is derived in detail in Ref.~\cite{RidleyQuantumProcesses} and in general obtained by counting the number of emitted photon modes per volume per unit energy:
\begin{equation}
    \label{eq:RidleyDOS}
    \rho(\text{E}) = \frac{n_{r}^{3}\omega^{2}}{\pi^{2}c^{3}\hbar} \,\, ,
\end{equation}
where $c$ is the speed of light in vacuum. The expression in Ref.~\cite{RidleyQuantumProcesses} is given in terms of the speed of light in a medium, $v_{l}$, which is related to $c$ and $n_{r}$ by $v_{l}=c/n_{r}$. Lastly, we have included a factor of two to account for left and right circular polarised light. Thus, overall the expression in Eq.~\eqref{eq:FermisGoldernRule} becomes
\begin{equation}
    \label{eq:RadLinearPol}
    W_{i\rightarrow f} = \frac{4e^{2}n_{r}\omega}{\epsilon_{0}\hbar m^{2}c^{3}} \frac{1}{4\pi}|\bra{i}\mathbf{a}\cdot\mathbf{p}\ket{f}|^{2}\delta(\hbar\omega - \text{E}_{i} - \text{E}_{f}) \,\, .
\end{equation}
So far we have only considered linearly polarised light propagating in one direction. To take into account all possible polarisations over propagation directions, we average Eq.~\eqref{eq:RadLinearPol} over a sphere. Calculating the average over solid angles of $|\mathbf{a}\cdot\mathbf{p}_{if}|^{2}$ we find
\begin{equation}
    \label{eq:AveragOverSolidAngles}
    \langle|\mathbf{a}\cdot\mathbf{p}_{if}|^{2}\rangle_{\text{av}} = \frac{4\pi}{3}|p_{if}|^{2} \,\, ,
\end{equation}
where $p_{if}$ are the so-called momentum matrix elements;
\begin{equation*}
    |p_{if}|^{2} = (\bra{i}\hat{p}_{x}\ket{f})^{2} + (\bra{i}\hat{p}_{y}\ket{f})^{2} + (\bra{i}\hat{p}_{z}\ket{f})^{2} \,\, .
\end{equation*}
Finally, using Eq.~\eqref{eq:AveragOverSolidAngles} in Eq.~\eqref{eq:RadLinearPol} we find that
\begin{equation}
    \label{eq:RadTransitionProb}
    W_{i\rightarrow f} = \frac{4e^{2}n_{r}\omega}{\epsilon_{0}\hbar m^{2}c^{3}} \frac{|p_{if}|^{2}}{3}\delta(\hbar\omega - \text{E}_{i} - \text{E}_{f}) \,\, .
\end{equation}

Following Ref.~\cite{RidleyQuantumProcesses}, to determine the total spontaneous radiative recombination rate we must sum over all possible energies of emitted photons. Thus we take the integral of Eq.~\eqref{eq:RadTransitionProb} over all possible values of $\hbar\omega$:
\begin{align}
    \label{eq:WijTotal}
    W_{i\rightarrow j} &= \int_{0}^{\infty} d(\hbar\omega) \frac{4e^{2}n_{r}\omega}{\epsilon_{0}\hbar m^{2}c^{3}} \frac{|p_{if}|^{2}}{3}\delta(\hbar\omega - \text{E}_{i} - \text{E}_{f}) \,\, , \nonumber \\
    &= \frac{1}{2}\int_{-\infty}^{\infty} d(\hbar\omega) \frac{4e^{2}n_{r}\omega}{\epsilon_{0}\hbar m^{2}c^{3}} \frac{|p_{if}|^{2}}{3}\delta(\hbar\omega - \text{E}_{i} - \text{E}_{f}) \,\, , \nonumber \\
    &= \frac{2}{3}\frac{n_{r} e^{2}}{m^{2}c^{3}\epsilon_{0}}\frac{(\text{E}_{i} - \text{E}_{j})}{\hbar^{2}}|p_{ij}|^{2} \,\, . 
\end{align}

Determining the momentum matrix elements requires knowledge of the wave functions of states $\ket{i}$ and $\ket{f}$. In our case we will use the TB wave functions and employ the method developed in Ref.~\cite{ScSc2006}. To do so we recast the momentum matrix elements as dipole matrix elements, $d_{ij} = \bra{i}\mathbf{a}\cdot\mathbf{\hat{r}}\ket{j}$, where $\mathbf{a}$ is the light polarisation vector as before and $\mathbf{\hat{r}}$ is a vector of the position operators. To do so, as an illustration we first look component-wise at the commutator of the Hamiltonian, $\hat{H}_{0}$, of an electron in a periodic potential, $V$, and the position operator for the $x$-direction. It can be shown that in general
\begin{equation*}
    \left[\hat{H}_{0},\hat{x}\right] = \left[\frac{\mathbf{\hat{p}}\cdot\mathbf{\hat{p}}}{2m} + V,\hat{x}\right] = -\frac{i\hbar}{m}\hat{p}_{x} \,\, ,
\end{equation*}
given that $[\hat{A}\hat{B},\hat{C}] = \hat{A}[\hat{B},\hat{C}] + [\hat{A},\hat{C}]\hat{B}$ and $\left[\hat{p}_{x},\hat{x}\right] = -i\hbar$. Applying this result to our system above, described by the Hamiltonian $\hat{H}$, we can find that 
\begin{equation*}
    -i\frac{\hbar}{m}\bra{i}\hat{p}_{x}\ket{j} = \bra{i}\left[\hat{H},\hat{x}\right]\ket{j} = (\text{E}_{i} - \text{E}_{j})\bra{i}\hat{x}\ket{j} \,\, ,
\end{equation*}
where $\ket{i}$ and $\ket{j}$ are eigenstates of $\hat{H}$. Thus,
\begin{equation*}
    \bra{i}\hat{p}_{x}\ket{j} = \frac{\left(\text{E}_{i}-\text{E}_{j}\right)}{-i\hbar}m\bra{i}\hat{x}\ket{j} \,\, .
\end{equation*}
Similar expressions hold for $\hat{y}$ and $\hat{z}$ and so overall we have that 
\begin{equation}
    \label{eq:DipMomRel}
    |p_{ij}|^{2} = \frac{\left(\text{E}_{i}-\text{E}_{j}\right)}{\hbar^{2}}^{2}m^{2}|d_{ij}|^{2} \,\, .
\end{equation}
Substituting Eq.~\eqref{eq:DipMomRel} back into Eq.~\eqref{eq:WijTotal} we find that the transition rate between states $i$ and $j$ is
\begin{equation*}
    W_{i\rightarrow j} = \frac{2}{3}\frac{n_{r} e^{2}}{c^{3}\epsilon_{0}}\frac{(\text{E}_{i} - \text{E}_{j})^{3}}{\hbar^{4}}|d_{ij}|^{2} \,\, .
\end{equation*}

To account for the temperature and carrier density dependence of the involved transition rates, the above expression is weighted by the factor $P_{ij}(n,p,T)=f^{e}_{i}(n,T)f^{h}_{j}(p,T)$, where $f^{e}_{i}$ and $f^{h}_{j}$ are Fermi factors that account for the occupation of states as temperature and carrier density vary. These factors are determined via Fermi-Dirac statistics~\cite{QuantumTheoryOfSolidsEoin}:
\begin{eqnarray}
    \label{eq:FermiDiracStats}
    f^{e}_{i}(n,T) = f^{e}(\text{E}_{i},n,T) = \frac{1}{1+\text{exp}[(\text{E}_{i}-\text{E}_{fn}(n))/k_{B}T]}\,\, , \nonumber \\
    f^{h}_{j}(p,T) = f^{h}(\text{E}_{j},p,T) = \frac{1}{1+\text{exp}[(\text{E}_{fp}(p)-\text{E}_{j})/k_{B}T]}\,\, .
\end{eqnarray}
The quasi-Fermi levels, $\text{E}_{fn}(n)$ and $\text{E}_{fp}(p)$ for electrons and holes, respectively, depend on the electron, $n$, and hole, $p$, carrier densities, while $k_{B}$ is the Boltzmann constant and $T$ is the system temperature. To determine the quasi-Fermi levels the electron and hole carrier densities are used in conjunction with the volume of the active region to find the number of carriers available to undergo recombination, $N$ for electrons (or $P$ for holes). Then, the expression
\begin{equation*}
    g(\text{E}_{fn}) = \int f^{e}(\text{E},\text{E}_{fn}) d\text{E} - N
\end{equation*}
is minimised to determine $\text{E}_{fn}$ for electrons and similarly for holes by replacing $\text{E}_{fn}$ by $\text{E}_{fp}$ and $N$ by $P$. Thus, each $f^{\alpha}$ in Eq.~\eqref{eq:FermiDiracStats} is a function of carrier density, $n$ (or $p$), and temperature, $T$. In what follows we assume that the material is electrically neutral (i.e. that $n = p$) as is commonly done in other works~\cite{KiSt2015,tHvO1983,DaYo2019scaling}.

Lastly, we sum over the available conduction and valence band states, taking into account the two possible spin states of the electrons, to find the total spontaneous radiative recombination rate:
\begin{equation}
    \label{eq:TotalSpontRadRate}
    R_{\text{sp}} = 2\frac{2}{3}\frac{ e^{2}}{\epsilon_{0}c^{3}\hbar^{4}}\sum_{i,j} P_{ij}(n,T)n_{r}(\lambda_{ij})(\Delta \text{E}_{ij})^{3} |d_{ij}|^{2} \,\, ,
\end{equation}
where $\Delta \text{E}_{ij} = \text{E}_{i} - \text{E}_{j}$ and we have included the wavelength dependence of the refractive index, $n_{r}(\lambda_{ij})$, with $\lambda_{ij} = hc/\Delta \text{E}_{ij}$. The wavelength/emission energy dependence of the refractive index is evaluated using a Sellmeier type description~\cite{PaWa2017};
\begin{equation*}
    n_{r}(\lambda_{ij}) = \sqrt{A + \frac{B\lambda_{ij}^{2}}{(\lambda_{ij}^{2}-C^{2})}} \,\, ,
\end{equation*}
where the Sellmeier coefficients, $A$, $B$, and $C$, are material dependent. In the study of InGaN/GaN QWs below, we have used the Sellmeier coefficients of the GaN barriers, as has been done previously for InGaN/GaN systems in Ref.~\cite{PaWa2017}. For GaN the values of these coefficients are $A = 5.15$, $B = 0.35$ and $C = 339.8$ nm~\cite{AnNa2003}. For the study of AlGaN/AlN QWs we have used values for AlN reported in Ref.~\cite{AnNa2003}: $A = 1.00$, $B = 3.12$ and $C = 138.0$ nm. 

Once $R_{\text{sp}}$ has been determined, the rate of spontaneous emission per unit volume ($r_{\text{sp}}^{\text{3D}}$) or unit (in-plane) area ($r_{\text{sp}}^{\text{2D}}$) of the well can be calculated via $r_{\text{sp}}^{\text{3D}} = R_{\text{sp}}/V_{QW}$ or $r_{\text{sp}}^{\text{2D}} = R_{\text{sp}}/A_{QW}$. In order to compare our results with reported experimental data we calculate both the sheet (2D) radiative recombination coefficient, $B_{\text{2D}}= r_{\text{sp}}^{\text{2D}}/np$ where $n$ and $p$ are the sheet carrier densities for electrons ($n$) and holes ($p$) in units of cm$^{-2}$, and the volume (3D) radiative recombination coefficient, $B_{\text{3D}}= r_{\text{sp}}^{\text{3D}}/np$ where $n$ and $p$ are the volume carrier densities for electrons and holes in units of cm$^{-3}$. \blue{A widely made assumption in the literature when modelling the recombination coefficients is that of charge neutrality in the well, i.e., that $n = p$~\cite{KiSt2015,tHvO1983,DaYo2019scaling}. We note that full LED structures in general employ multi-quantum well (MQW) systems which may lead to the situation of inhomogeneous current spreading between the wells under \emph{electrical injection}, i.e., $n \neq p$. However, resonant optical excitation of (In,Ga)N QWs is often applied to gain insight into the fundamental properties of recombination processes. In such experiments the wavelength of the excitation light pulse is tuned to the band gap of a well so that carriers are directly excited within the well and not the barrier material. Under these conditions the assumption that $n = p$ is expected to be valid and provides an ideal set up to compare theory and experiment~\cite{ShMu2007,BiNi2013}.}

\blue{As previously mentioned, the TB energies and wave functions determined in Sec.~\ref{sec:TheoryTightbinding} are used to calculate the dipole matrix elements, or, in general, interaction matrix elements. In empirical atomistic models such as the TB framework employed here, the alloy ``properties'' (e.g. gradient in composition, random vs. clustered alloys, ``bleeding'' of atoms into the barrier material) are input to the calculation of the electronic structure. Thus, in general one requires some information about the alloy distribution which is usually based on experimental data. In line with experimental atom probe tomography (APT) data, the cation species in InGaN and AlGaN QWs are randomly distributed~\cite{OlBe2010,SmKa2003,RiMa2016}. This approach is used elsewhere in the literature in empirical atomistic studies of random alloys~\cite{ChLi2010,UsOR2014,SmWa1989,Farkas2021}, and also in modified continuum-based approaches~\cite{JoTe2017}.}

\blue{Building on such an assumption, our previous theoretical studies provided good agreement with experimental and other literature theoretical data on, for instance, low temperature FWHM values of InGaN QWs over a wide composition range~\cite{TaDa2020}. Given that a large FWHM value at low temperatures in directly connected to fluctuations in the band gap of the well (thus, is affected by local alloy configurations), it is important to capture configurations that result in, for instance, strong carrier localisation effects and thus low band gap values. As such, one may find ``extreme'' alloy configurations that significantly deviate from, for example, special quasi-random structures often considered in DFT calculations~\cite{WeFe1990}. The significance of these ``extreme'', or potentially, ``outlier'', configurations can be judged in several ways; (i) by analysing their total energy (e.g. are these configurations more unlikely than others), (ii) by their statistical significance, e.g. how often these configurations are found when studying a large number of random alloy configurations and (iii) by comparison with experimental studies. To give an example for (iii), if an ``extreme'' configuration leads to, for instance, faster radiative recombination in comparison to other configurations, then such a configuration can be expected to dominate the recombination processes, as it out performs the other processes and can provide a characteristic feature in, e.g., time resolved PL studies. Thus, even when one can study only a limited number of configurations, comparison with experimental data can shed light onto the importance of the different configurations randomly generated in the theoretical model. Therefore, we follow here the widely used assumption from the literature and investigate the impact of the alloy microstructure on the results by carrying out the calculations described above for a number of alloy configurations. For each configuration, the targeted carrier density and temperature are used to determine the quasi-Fermi level, as described above. Subsequently, the $B$ coefficient is averaged over the number of investigated configurations.}

\section{Calculation of Auger recombination rate}
\label{sec:RecombinationAuger}

\begin{figure}[t!]
    \centering
    \begin{tabular}{cc}
        \includegraphics[width=0.45\columnwidth]{augerSchemaEEH.eps} & \includegraphics[width=0.45\columnwidth]{augerSchemaHHE.eps}
    \end{tabular}
    \caption{a) Schematic representation of the direct electron-electron-hole ($eeh$) Auger recombination process between electron states $\mathbf{1}$ and $\mathbf{2}$, and hole state $\mathbf{3}$. State $\mathbf{4}$ is an unoccupied electron state that the electron in state $\mathbf{2}$ is scattered into. b) Schematic representation of the direct hole-hole-electron ($hhe$) Auger recombination process between electron state $\mathbf{1}$ and hole states $\mathbf{2}$ and $\mathbf{3}$. State $\mathbf{4}$ is an unoccupied hole state that the hole in state $\mathbf{2}$ is scattered into.}
    \label{fig:augerRecombSchema}
\end{figure}

As described above in Sec.~\ref{sec:LEDsModellingEfficiency}, and illustrated schematically in Fig.~\ref{fig:augerRecombSchema}, Auger recombination is a three particle process. Several Auger processes exist, including direct and indirect (phonon-assisted)~\cite{KiSt2015}, as well as trap/defect-assisted Auger recombination~\cite{ZhTu2022,DaYo2019TAAR}. We focus here on direct Auger recombination in the presence of alloy disorder (so-called alloy-enhanced direct Auger). It has been found theoretically that alloy disorder can increase the rate of direct Auger recombination in bulk InGaN, which is essentially negligible in the absence of alloy disorder~\cite{KiSt2015}. The alloy-enhanced Auger recombination rate was also found to be larger in bulk InGaN compared to the phonon-assisted Auger rate in the band gap range of 1.9 -- 3.0 eV, thus we expect that alloy-enhanced Auger recombination plays an important role in InGaN QWs too. For the study of AlGaN QWs with band gaps outside this range below, a similar analysis of the competition between alloy-enhanced direct Auger and indirect Auger recombination in bulk AlGaN or AlGaN QWs has not been reported. However, given the observation of carrier localisation effects in AlGaN QWs (see Sec.~\ref{sec:LEDsIIINMaterials}), we expect that the alloy-enhanced Auger process may also play an important role. Trap-assisted Auger recombination has recently been explored as an explanation for the large overall (i.e. sum of direct, indirect, alloy-enhanced, trap-assisted contributions) Auger rates observed experimentally but the nature and origin of the ``trap'' states has yet to be determined~\cite{ZhTu2022,DaYo2019TAAR}. Taking all of this together, we focus in the following on alloy-enhanced Auger, and other processes may be targeted in future studies.

As an illustration, in the follow analysis we focus on the $eeh$ Auger process, noting that similar arguments hold for the $hhe$ Auger process. In $eeh$ Auger recombination an electron and hole pair recombine and interact with another electron, but rather than emitting a photon another electron is scattered into an energetically higher lying excited state. In order to calculate the rate of Auger recombination, we consider the interaction of the electrons via the Coulomb force to be a perturbation to the system. Building on Fermi's Golden Rule, introduced above, the transition rate between initial state $i$ and final state $j$ is determined. We mainly follow the approach of Ref.~\cite{KiSt2015} but highlight changes due to the presence of alloy disorder. 

Overall, the Auger recombination rate is given by
\begin{equation}
    R = 2 \frac{2\pi}{\hbar} \sum_{\textbf{1234}} P(n,T)\left|M_{\textbf{1234}} \right|^{2} \delta\left( \text{E}_{\textbf{1}} + \text{E}_{\textbf{2}} - \text{E}_{\textbf{3}} - \text{E}_{\textbf{4}} \right)\,\, , \label{eq:RateAuger}
\end{equation}
where the factor of 2 accounts for spin. In a semiconductor material where translational invariance is preserved, the bold indices denote both band index and $\mathbf{k}$-point. Given that random alloy fluctuations lead to strong carrier localization effects in InGaN/GaN QW systems~\cite{WaGo2011,ScCa2015,DaSc2016}, and similar effects have been observed in AlGaN/AlN QWs~\cite{FrNi2020PLSShape,FrNi2020PLBroadening}, the alloy microstructure breaks the translational symmetry of the system so that the wave vector $\mathbf{k}$ is no longer a ``good'' quantum number. The $\delta$-function assures energy conservation. The pre-factor
\begin{equation*}
    P(n,T)=f^{e}_{\textbf{1}}(n,T)f^{e}_{\textbf{2}}(n,T)f^{h}_{\textbf{3}}(p,T)f^{e}_{\textbf{4}}(n,T) \,\, ,
\end{equation*}
accounts for state occupations, and ensures that transitions only occur between occupied and empty states; $f$ denotes the occupation number according to Fermi-Dirac statistics, as described above. The Coulomb interaction between carriers in different states, $\psi_\textbf{m}$, is described by the matrix elements $M_{\textbf{1234}}$ such that
\begin{eqnarray*}
    \left| M_{\textbf{1234}} \right|^{2} = \left|M_{\textbf{1234}}^{d}- M_{\textbf{1234}}^{x} \right|^{2}+ \left| M_{\textbf{1234}}^{d} \right|^{2} + \left|M_{\textbf{1234}}^{x} \right|^{2}\,\, ,
\end{eqnarray*}
where $M^{d}_{\textbf{1234}}$ ($M^{x}_{\textbf{1234}}$) are the direct (exchange) screened Coulomb matrix elements:
\begin{eqnarray*}
    M_{\textbf{1234}}^{d} &=& \langle \psi_{\textbf{1}}\psi_{\textbf{2}}\left| W \right| \psi_{\textbf{3}} \psi_{\textbf{4}} \rangle\,\, ,\\
    M_{\textbf{1234}}^{x} &=& \langle \psi_{\textbf{1}}\psi_{\textbf{2}} \left| W \right| \psi_{\textbf{4}}\psi_{\textbf{3}} \rangle\,\, .
\end{eqnarray*}
The direct, screened Coulomb interaction (noting that similar arguments hold for the exchange term) may thus be written as:
\begin{equation}
    M_{\textbf{1234}}^{d} =\int\int d^3rd^3r'\psi^{*}_{\textbf{1}}(\mathbf{r})\psi^{*}_{\textbf{2}}(\mathbf{r}')\psi_{\textbf{3}}(\mathbf{r}')\psi_{\textbf{4}}(\mathbf{r})\frac{e^{2}}{4\pi\epsilon_{0}\epsilon_{\infty}}\frac{e^{-\alpha|\mathbf{r}-\mathbf{r}'|}}{|\mathbf{r}-\mathbf{r}'|}\,. 
    \label{eq:CME}
\end{equation}
The wave functions, $\psi_\textbf{m}$, required for evaluating the above matrix elements are obtained from our TB model.  Furthermore, when evaluating $M_{\textbf{1234}}^{\lambda}$, with $\lambda=d$ or $\lambda=x$ for direct or exchange matrix elements respectively, we assume a position independent, and constant, background dielectric function~\cite{HaugAndKochThird}. For the study of InGaN/GaN QWs below, we take the GaN barrier value, $\epsilon_{\infty}(\text{GaN}) = 5.5$, as has been done in Ref.~\cite{KiSt2015}. This approximation is often applied in Auger matrix element calculations~\cite{ZeAn1995,AnZe1997}. In the study of 50\% Al AlGaN/AlN QWs we use a 50/50 ratio of AlN and GaN of $\epsilon_{\infty}(\text{Al}_{0.5}\text{Ga}_{0.5}\text{N}) = 5.0$, based on a dielectric constant of $\epsilon_{\infty}(\text{AlN}) = 4.64$ for AlN. This value for AlN is determined from the average of the values reported in Ref.~\cite{BeFi1998DielectricConstants} for two theoretical investigations ($\epsilon_{\infty}^{\text{theo, a}}(\text{AlN}) = 4.61$; $\epsilon_{\infty}^{\text{theo, b}}(\text{AlN}) = 4.62$) and an experimental investigation ($\epsilon_{\infty}^{\text{exp}}(\text{AlN}) = 4.68$).  

In Ref.~\cite{KiSt2015}, the inverse screening length $\alpha$ in Eq.~\eqref{eq:CME} is determined separately for electrons ($\alpha_e$), and holes ($\alpha_h$), by using the Debye-H{\"u}ckel equation: 
\begin{eqnarray*}
    \alpha_{e} = \sqrt{\frac{4 \pi n e^{2} }{ \epsilon_{0} \epsilon_{\infty} k_{B} T } } \,\, , \\
    \alpha_{h} = \sqrt{\frac{4 \pi p e^{2} }{ \epsilon_{0} \epsilon_{\infty} k_{B} T } } \,\, ,
\end{eqnarray*}
where $k_B$ is the Boltzmann constant, $T$ is the temperature, $\epsilon_0$ is the vacuum permittivity, $n$ ($p$) is the electron (hole) carrier density and $e$ is the elementary charge. The total inverse free-carrier screening length is then given by $\alpha^2=\alpha^2_e+\alpha^2_h$, which is used in Eq.~\eqref{eq:CME}.

To calculate $M_{\textbf{1234}}^{\alpha}$, we largely follow the procedure detailed in Ref.~\cite{ScSc2006} to obtain Coulomb matrix elements from TB wave functions. The TB wave functions may be written as:
\begin{eqnarray}
    \psi_{\textbf{m}}(\mathbf{r})=\sum_{\mathbf{R},\sigma}c^{\textbf{m}}_{\mathbf{R},\sigma}\phi_{\mathbf{R},\sigma}(\mathbf{r})\,\,
    \label{eq:TBWF}
\end{eqnarray} 
where $\phi_{\mathbf{R},\sigma}(\mathbf{r})$ are the localised, atomic-like orbitals, $\sigma=\left\{s,p_x,p_y,p_z\right\}$, located at a lattice site, $\mathbf{R}$, and $c^{\textbf{m}}_{\mathbf{R},\sigma}$ are the expansion coefficients of the orbital $\sigma$ at site $\mathbf{R}$. Equipped with this information, the matrix elements are evaluated as
\begin{eqnarray*}
    \nonumber 
    M_{\textbf{1234}}^{d}
    &=&\sum_{\mathbf{R},\mathbf{R}'}\sum_{\alpha,\beta}c^{\textbf{1}*}_{\mathbf{R},\alpha}c^{\textbf{2}*}_{\mathbf{R}',\beta}c^{\textbf{3}}_{\mathbf{R}',\beta}c^{\textbf{4}}_{\mathbf{R},\alpha}W(\mathbf{R}-\mathbf{R'}) \,\, ,
\end{eqnarray*}
with
\begin{eqnarray*}
    W(\mathbf{R}-\mathbf{R}')=\frac{e^{2}}{4\pi\epsilon_{0}\epsilon_{\infty}}\frac{e^{-\alpha|\mathbf{R}-\mathbf{R}'|}}{|\mathbf{R}-\mathbf{R}'|}\,
    ,\,\,\text{for}\quad\mathbf{R}\neq\mathbf{R}'
\end{eqnarray*}
and
\begin{eqnarray*}
    W(\mathbf{0})=\frac{1}{v^2_\text{tet}}\int_{\text{tet}}\frac{e^{2}}{4\pi\epsilon_{0}}\frac{1}{|\mathbf{r}-\mathbf{r}'|}\approx W_0\, ,\,\text{when}\,\mathbf{R}=\mathbf{R}'\,.
\end{eqnarray*}

The on-site energy, $W_0$, can be calculated quasi-analytically, as described in Ref.~\cite{ScSc2006}, and in more detail in Appendix~\ref{chp:AppendixDeterminationOfOnsiteEnergies}, by integrating over a sphere of volume equal to that of the local tetrahedron formed  by the nearest neighbors of an atom in a wurtzite unit cell. This approximation is widely used in determining the Coulomb matrix elements from TB wave functions~\cite{ScSc2006,ShCh2005,ScBa2012,BaSc2013}, building on the finding that in general the dominant contributions to the Coulomb matrix elements stem from the long-ranged part; therefore, the localised orbitals, $\phi_{\mathbf{R},\sigma}(\mathbf{r})$, underlying the TB model can be treated, in a first approximation, as point charges, and thus the explicit structure of the localised orbitals is of secondary importance~\cite{ScSc2006,ShCh2005,BaSc2005,Zielinski2012}. The impact of the on-site Coulomb energy, $W_0$, on the Auger rate is also discussed in more detail in Appendix~\ref{chp:AppendixDeterminationOfOnsiteEnergies}.

Finally, and following previous work on calculating the Auger rate in III-N materials~\cite{JoTe2017,KiYa2012}, the fourth state ($\textbf{4}$) in Eq.~\eqref{eq:RateAuger}, into which the carriers can scatter, is always considered to be unoccupied. In addition to the pre-evaluated large number of electron and hole states ($\geq 100$ for each carrier type) calculated for the results presented below, we have also performed further calculations to investigate the impact of this energetically higher/lower lying (fourth) state on the results. Firstly, our calculations revealed a continuum of states (energy separation $\leq 0.5$ meV) in the range of the magnitude of the band gap above (below) the conduction (valence) ``band'' edge. Testing Auger rates with different ``final'' fourth states from this continuum revealed that the specific choice of the fourth state was of secondary importance to the results when compared to the spread in values introduced by alloy fluctuations. Thus, a pre-calculated continuum state is used for each alloy composition investigated (10\%, 15\% and 25\% for the study of InGaN QWs and 50\% for the study of AlGaN QWs). Additionally, despite the existence of such a continuum of states, the energy conservation requirement in the Auger rate calculation, given by the $\delta$-function in Eq.~\eqref{eq:RateAuger}, is never exactly numerically satisfied. Therefore, the $\delta$-function was replaced by a Gaussian function and the impact of the broadening parameter on the results was tested. These investigations showed that a broadening parameter of 1 meV is sufficient to achieve converged results, i.e., choosing a smaller value did not change the results.

Equipped with this information, we obtain the Auger coefficients from the Auger rate, $R$, using \mbox{$C=R/(Vn^{3})$} (see Sec.~\ref{sec:LEDsModellingEfficiency}). This is repeated for the two possible Auger recombination processes, $eeh$ and $hhe$, also discussed above. Here, $V$ is the volume of the well, and $n$ is the carrier density. As was assumed for the determination of the radiative recombination rate, for all Auger recombination calculations we let $n = p$, i.e., we assume that the structure is electrically neutral, as is done in other works~\cite{KiSt2015,tHvO1983,DaYo2019scaling}. Unlike in the case of the radiative recombination rate above, we do not calculate sheet (2D) Auger recombination rates but only the volume (3D) Auger recombination rate, as this is the quantity most often reported in the experimental literature~\cite{GaDr2012,PiprekReview2020,PiLi2017LLT2}. The total Auger rate is then calculated as $C_{\text{tot}} = C_{eeh} + C_{hhe}$.

As we have seen in our discussion of the calculation of the radiative and Auger recombination rates above, we expect that both rates will depend on temperature and carrier density. However, as we have discussed above in Sec.~\ref{sec:LEDsIIINMaterials}, as carrier density is increased in III-N heterostructures, the polarisation field will be screened, and as a result the electronic structure may also change. In the next section we will discuss this effect and a method we have developed to take this change into account.

\section{Polarisation field screening}
\label{sec:RecombinationScreening}

As we have previously discussed, $c$-plane III-N heterostructures exhibit strong built-in polarisation fields. These fields lead to the spatial separation of charge carriers on opposite QW interfaces, which is part of the so-called quantum confined Stark effect (QCSE) (see for example Fig.~\ref{fig:BandsFlatTilted} in Sec.~\ref{sec:LEDsIIINMaterials}). However, the separation of electron and hole charge densities leads to an electric field pointing in the direction opposite to the built-in polarisation field. The potential arising from this ``new'' field counteracts the polarisation field, diminishing it, thus the introduction of charge carriers to the system ``screens'' the polarisation field. The screening of the polarisation field will depend on the charge (i.e. carrier) density, and can change the electronic structure of the QW significantly (e.g. reducing the spatial separation of carriers/QCSE).

For \emph{low} carrier densities in III-N QWs, the TB model outlined above is sufficient to determine the electronic structure, as the impact of screening is essentially negligible. As a result, it is possible to use the obtained single particle TB states and energies to calculate, e.g., radiative recombination rates~\cite{McTa2020}, excitonic effects~\cite{ScTa2015}, or Auger recombination rates at low carrier densities~\cite{McKi2022}. However, at \emph{higher} carrier densities this is no longer the case as the electrostatic built-in field will be noticeably screened. Theoretical and experimental works investigating InGaN/GaN QWs indicate that, for those systems, the screening effect becomes important at carrier densities on the order of \mbox{$n=1\times10^{19}$} cm$^{-3}$ and above~\cite{DaYo2020review,MuHa2022,HaCh2022}. Similarly, it has been suggested in theoretical studies of AlGaN QW systems that screening effects become relevant above $n=1\times10^{19}$ cm$^{-3}$~\cite{PaCh1998}. Thus, a self-consistent approach is required to account for the screening of the built-in field and the change in the electronic structure.

However, as established above, to capture carrier localisation effects in III-N QWs, very large simulation supercells and an atomistic resolution are required, as provided by the TB model employed here. In general it is numerically extremely demanding to carry out such a self-consistent \sch-Poisson calculation, but especially so in a TB framework, since it has to be performed for each carrier density and for a large number of electronic states when targeting, for instance, an investigation of Auger recombination rates in large supercells. To keep the information about alloy-induced carrier localization effects, but at the same time take built-in field screening effects into account, we proceed as follows. 

We start from a one dimensional Schr\"{o}dinger equation describing a wurtzite III-N QW grown along the $c$-axis:
\begin{equation}
    \left[-\frac{\hbar^{2}}{2m_{b}}\frac{\partial^{2}}{\partial z^{2}}+U^{b}(z)\right]\psi(z)=\text{E}_b\psi(z)\,\, ,
    \label{eq:1dS}
\end{equation}
where the confining potential $U^{b}(z)$ for electrons, $b=e$, or holes, $b=h$, is given by~\cite{NiGa2005} 
\begin{equation}
    U^{b}(z)=U^{b}_{0}(z)+U^{b}_{p}(z)+U^{b}_{\text{scr}}(z)\,\, .
\end{equation}
The different contributions to $U^{b}(z)$ are the bare (band offset) confining potential, $U^{b}_{0}(z)$, the intrinsic (spontaneous and piezoelectric related) electrostatic potential, $U^{b}_{p}(z)$, and the screening potential $U^{b}_{\text{scr}}(z)$, respectively. Here, $U^{b}_{\text{scr}}(z)$ arises from the spatial separation of electron and hole wave functions in the well due to the presence of the intrinsic built-in field. Following Haug and Koch~\cite{HaugAndKochThird}, the screening potential is calculated via:
\begin{equation}
    U^{e,h}_{\text{scr}}(z)=\pm\frac{e^{2}n_{\text{sys}}}{2\epsilon_{0}\epsilon_{r}}\int dz'\left[|\psi^{e}(z)|^{2}-|\psi^{h}(z)|^{2}\right]|z-z'|\, .
    \label{eq:ScP}
\end{equation}
Eqs.~\eqref{eq:1dS} -~\eqref{eq:ScP} are evaluated self-consistently for a given carrier density, $n_{\text{sys}}$, within the well. We use band offsets, average strain fields, spontaneous and piezoelectric coefficients as well as effective masses from our TB model as input for Eqs. \eqref{eq:1dS} - \eqref{eq:ScP} for the studied $c$-plane III-N QW systems. Solving these equations self-consistently allows us to determine $U^{e,h}_{\text{scr}}(z)$ for a given carrier density, $n_{\text{sys}}$. As a result, one can determine by how much the intrinsic electrostatic built-in field is reduced at a given $n_{\text{sys}}$. This information is then fed back into the TB model by adjusting/scaling the built-in potential across the well and diagonalising the Hamiltonian for each carrier density $n_{\text{sys}}$. In doing so, we obtain, for a given $n_{\text{sys}}$, the (screened) TB wave functions and energies for the studied supercell, while still including local alloy contributions. The obtained wave functions and energies then form the input for the radiative and Auger recombination rate calculations detailed above. 

We note that a similar screening approach has been used by Auf der Maur \emph{et al.}~\cite{AuPe2016}, to study the impact of alloy fluctuations on the radiative recombination rate in $c$-plane InGaN/GaN QWs.  Although we will return to the discussion of the impact of polarisation field screening on the electronic structure below, for InGaN and AlGaN QW systems we find that (i) we obtain an onset of the built-in potential screening at carrier densities similar to experimental and theoretical studies~\cite{DaYo2020review,PaCh1998} and (ii) the expected behavior of the electron and hole wave function separation is observed (i.e. that the separation decreases) with increasing carrier density. Thus, we expect that the above detailed approach is sufficient to capture the main effects of the built-in field screening on the results. Future studies can target a self-consistent 3D approach. Having established the theoretical framework used to determine the electronic structure in $c$-plane III-N QWs, along with the radiative and Auger recombination rates, we now present results for, in a first step, InGaN/GaN QWs, before turning to AlGaN/AlN QWs. 