Impact of Charge Transitions at Atomic Defect Sites on Electronic Device Performance

Chapter 3 Computational Methods

This chapter outlines the computational methods used in this thesis to investigate the properties of atomic defects in semiconductors and insulators. In particular, density functional theory (DFT), molecular dynamics (MD) and machine learning interatomic potential (MLIP) are discussed.

3.1 Density functional theory (DFT)

The quantum mechanical description of electrons in a material is governed by the Schrödinger equation. Solving this equation for a solid state material system is analytically impossible and computationally infeasible already for a system size of a few atoms. Every electron and ion in the system interacts with each other and different quantum mechanical effects have to be considered, leading to coupled high-dimensional partial equations. Therefore, such a problem has to be solved approximately. An established approach that is computationally efficient while still giving an accurate description of the electronic structure is DFT. This modeling method has been proven for decades to reliably reproduce and predict material properties based on first-principles calculations. The theorems and approximations that led to this theory are outlined in the following.

3.1.1 Hohenberg-Kohn theorems

A significant step towards the development of DFT was made by Hohenberg and Kohn in 1964 with the formulation of two groundbreaking theorems [169], that demonstrated that the ground state properties of an atomic system are uniquely defined by the electron density $\rho$. Following the Born-Oppenheimer approximation [106] as described in Section 2.1.2, the electronic part of the Hamiltonian $H_{\mathrm{el}}$ can be separated to define a Schrödinger equation for the electrons as given in Eq. (2.4). $H_{\mathrm{el}}$ includes three terms: the kinetic energy of the electrons $T$, a potential for the interaction between the electrons $V_{\mathrm{ee}}$ and a term for an external potential $V_{\mathrm{ex}}$, representing all interactions of the electrons with positively charged nuclei as well as possible other external effects such as an applied electric field. $V_{\mathrm{ex}}$ thus depends on the positions of the nuclei ${\mathit{R}}_j$, but also on their total number $M$ and charge $Z$. For a system with $N$ electrons, $H_{\mathrm{el}}$ is given in atomic units by

$$\begin{array}{}\text{(3.1)}& H_{\mathrm{el}}=T+V_{\mathrm{ee}}+V_{\mathrm{ex}}=\sum _i^N\frac{{\mathrm{\nabla }}_{{\mathit{r}}_i}^2}{2}+\sum _{i\ne j}^N\frac{1}{|{\mathit{r}}_j-{\mathit{r}}_i|}+\sum _{i,j}^{N,M}V({\mathit{R}}_i,{\mathit{r}}_i),\end{array}$$

where ${\mathrm{\nabla }}_{{\mathit{r}}_i}$ is the Laplace operator.

The first Hohenberg and Kohn theorem states that the external potential $V_{\mathrm{ex}}$ is a unique functional of the electron density $\rho (\mathit{r})$, which consequently means that the full many particle ground state is a unique functional of $\rho (\mathit{r})$ [169]. This was proven by a reductio ad absurdum argument, showing that two different external potentials cannot yield the same ground state electron density ${\rho }_0$.

A functional of a certain parameter will be denoted with square brackets in the following. Because the total energy of the system in the ground state is a functional of ${\rho }_0$, the same must be true for all of its individual elements.

$$\begin{array}{}\text{(3.2)}& E_0[{\rho }_0]=E_{\mathrm{ex}}[{\rho }_0]+T[{\rho }_0]+E_{\mathrm{ee}}[{\rho }_0]=\int {\rho }_0(\mathit{r})V_{\mathrm{ex}}\mathrm{d}\mathit{r}+E_{\mathrm{HK}}[{\rho }_0].\end{array}$$

Determining ${\rho }_0$ is therefore sufficient to obtain the full many particle ground state energy and, by extension, all ground state properties of the system. In the above equation, the universal functional $E_{\mathrm{HK}}[{\rho }_0]=T[{\rho }_0]+E_{\mathrm{ee}}[{\rho }_0]$ was separated and contains both classical and non-classical contributions [170]. Universal here means that the expressions for $T$ and $E_{\mathrm{ee}}$ are independent of the ionic system for the same number of electrons.

The exact form of the functional $E_{\mathrm{HK}}[{\rho }_0]$ is still to be discovered. It includes not only the Coulomb repulsion between the electrons which can easily be expressed, but also all non-classical contributions such as the exchange energy due to the Pauli exclusion principle and effects of Coulomb correlation. If the expression of the functional $E_{\mathrm{HK}}[{\rho }_0]$ would be known exactly, the Schrödinger equation could be solved for any system. As a first step, at least the classical Coulomb part can be separated from $V_{\mathrm{ee}}[\rho ]$ to make the functional more accessible [170].

$$\begin{array}{}\text{(3.3)}& E_{\mathrm{ee}}[\rho ]=\frac{1}{2}\int \int \frac{\rho ({\mathit{r}}_1)\rho ({\mathit{r}}_2)}{{\mathit{r}}_{12}}\mathrm{d}{\mathit{r}}_1\mathrm{d}{\mathit{r}}_2+E_{\mathrm{ncl}}[\rho ],\end{array}$$

where the first term corresponds to the Hartree energy and the non-classical many-body correlation effects such as self-interaction, exchange and Coulomb correlation are all packed in the non-classical $E_{\mathrm{ncl}}[\rho ]$.

The second Hohenberg and Kohn theorem states that only for the true ground state density ${\rho }_0$, $E_{\mathrm{HK}}({\rho }_0)$ gives the minimum energy of the system $E_0$. This was proven by applying the variational principle, showing that for slight deviations of ${\rho }_0$ the total energy as a function of the electron density increases. Mathematically, this can be written as $⟨{\psi }_0|H_{\mathrm{el}}|{\psi }_0⟩=E_0[{\rho }_0]\le E[\tilde{\rho }]$ which is true for any trial density $\tilde{\rho }\ge 0$, $\int \tilde{\rho }(\mathit{r})d\mathit{r}=N$ [170].

The Hohenberg-Kohn theorems give the basic principle of DFT and motivate the search for ${\rho }_0$. As soon as ${\rho }_0$ is obtained, the Hamiltonian is uniquely determined and different material properties can be derived from ${\rho }_0$. Because ${\rho }_0$ only depends on three spatial coordinates $\mathit{r}$ (${\rho }_0(\mathit{r})$) in contrast to the 3$N$ coordinates of the complete electron wave function, DFT significantly simplifies the computational challenge and allows for the investigation of more complex systems such as defective solids. Still, to make the solution of a system as described by Eq. (3.2) computationally feasible, certain additional methods and approximations are commonly applied which will be described in the following.

3.1.2 Kohn-Sham equations

Due to the electron-electron interactions in $V_{\mathrm{ee}}$, the electronic part of the Hamiltonian is not further separable. However, the complex problem can be simplified by an approach that was proposed by Kohn and Sham in 1965 [171]. The main idea is to introduce non-interacting particles ${\varphi }_i$ that generate the same electron density as the given system with interacting particles

$$\begin{array}{}\text{(3.4)}& {\rho }_0(\mathit{r})=\sum _i^N|{\varphi }_i(\mathit{r}){|}^2.\end{array}$$

By construction, ${\rho }_0$ of the fictitious system is equal to the electron density of the real interacting system ${\rho }_s(r)$. This is achieved by the introduction of an effective external potential $v_{\mathrm{eff}}(\mathit{r})$ that simulates all effects of the former interacting electrons. In addition to the separation of the interaction functional as given in Eq. (3.3), the functional of the kinetic energy $T[\rho ]$ in Eq. (3.2) can now be split up in a term that only considers the kinetic energy of the non-interacting particles $T_{\mathrm{s}}[\rho ]$ and an additional term representing the remaining correlation interactions $T_{\mathrm{C}}[\rho ]$.

$$\begin{array}{}\text{(3.5)}& T[\rho ]=T_{\mathrm{s}}[\rho ]+T_{\mathrm{C}}[\rho ]=-\frac{1}{2}\sum _i^N⟨{\varphi }_i|{\mathrm{\nabla }}^2|{\varphi }_i⟩+T_{\mathrm{C}}[\rho ].\end{array}$$

As a result, also the full functional $E[\rho ]$ is separable into classical terms (kinetic energy, Hartree energy and external potential energy) and the exchange-correlation energy $E_{\mathrm{XC}}$, which must be approximated. In this exchange-correlation energy, all quantum mechanical effects of self-interaction corrections, exchange and correlation are incorporated as well as a portion of the kinetic energy [172].

$$\begin{array}{}\text{(3.6)}& E_{\mathrm{XC}}[\rho ]=T_{\mathrm{C}}[\rho ]+E_{\mathrm{ncl}}[\rho ].\end{array}$$

The resulting Schrödinger-like equations for the $N$ single particles can be derived by applying the variational principle to minimize $E[\rho (\mathit{r})]$ in the space of non-interacting orbitals ${\varphi }_i$ under the condition $⟨{\varphi }_i|{\varphi }_j⟩={\delta }_{ij}$ [172, 173, 171].

$$\begin{array}{}\text{(3.7)}& (-\frac{{\mathrm{\nabla }}^2}{2}+\int \frac{\rho ({\mathit{r}}_2)}{r_{12}}d{\mathit{r}}_2+V_{\mathrm{XC}}({\mathit{r}}_1)-\sum _A^M\frac{Z_{\mathrm{A}}}{r_{1\mathrm{A}}}){\varphi }_i=(-\frac{{\mathrm{\nabla }}^2}{2}+v_{\mathrm{eff}}(\mathit{r})){\varphi }_i={ϵ}_i{\varphi }_i(\mathit{r}).\end{array}$$

These equations include all interaction functionals in an effective potential and are given for the non-interacting particles with wave functions ${\varphi }_i$ and energies ${ϵ}_i$. The effective potential includes the Hartree potential for the electrostatic interaction between the electrons, the exchange-correlation potential and the external potential due to electrostatic interactions of the electrons with the nuclei.

During the derivation, every component that can not be solved directly was packed into the unknown $V_{\mathrm{XC}}$. Because the exact form of $V_{\mathrm{XC}}$ is not available, this term is simply defined as the functional derivative of the corresponding to the exchange-correlation energy $E_{\mathrm{XC}}$

$$\begin{array}{}\text{(3.8)}& V_{\mathrm{XC}}=\frac{\delta E_{\mathrm{XC}}}{\delta \rho }.\end{array}$$

If this potential where known, the Kohn-Sham equation would lead to the correct ground state energy and electron density of the Schrödinger equation [172]. Since this is not the case, significant efforts have been made to reasonably approximate the exchange-correlation effects, which will be outlined in the next section.

3.1.3 Exchange-correlation functionals

The foundation of many approximations for the exchange-correlation energy is based on the treatment of the electron density as a locally homogeneous electron gas. This concept was already introduced in the original work of Kohn and Sham [171] and is nowadays known as the local density approximation (LDA). The corresponding functional is given by

$$\begin{array}{}\text{(3.9)}& E_{\mathrm{XC}}^{\mathrm{LDA}}[\rho ]=\int \rho (\mathit{r}){ϵ}_{\mathrm{XC}}(\rho (\mathit{r}))\mathrm{d}\mathit{r},\end{array}$$

where ${ϵ}_{\mathrm{XC}}(\rho )$ is the exchange and correlation energy per electron of a uniform electron gas with density $\rho$, which can be split in an exchange and correlation part ${ϵ}_{\mathrm{XC}}(\rho )={ϵ}_{\mathrm{X}}(\rho )+{ϵ}_{\mathrm{C}}(\rho )$. Furthermore, an analytical expression for the exchange part is given by ${ϵ}_{\mathrm{X}}(\rho (\mathit{r}))\propto \rho (\mathit{r}{)}^{1/3}$ and highly accurate numerical simulations based on quantum Monte-Carlo simulations exist for the calculation of ${ϵ}_{\mathrm{X}}$ [173]. Even though this approach of a homogeneous electron gas is quite simplistic, it has been proven to effectively capture certain molecule properties, such as equilibrium structures and harmonic frequencies [174].

However, LDA cannot account for rapidly varying electron densities, which description is crucial for analyzing localized charges at defect sites – one of the key focuses of this thesis. Additionally, LDA is known to overestimate binding energies and underestimate electronic band gaps. A significant improvement in the treatment of $E_{\mathrm{XC}}[\rho ]$ can be achieved by including the gradient of the electron density $\mathrm{\nabla }\rho$ in the definition of the functional.

Such a treatment is known as the generalized gradient approximation (GGA) and is generally given by

$$\begin{array}{}\text{(3.10)}& E_{\mathrm{XC}}^{\mathrm{GGA}}[\rho ]=\int f(\rho (\mathit{r}),\mathrm{\nabla }\rho (\mathit{r}))\mathrm{d}\mathit{r},\end{array}$$

where $f$ is a function depending on the electron density and its gradient. Various ways to parametrize $f(\rho (\mathit{r}),\mathrm{\nabla }\rho (\mathit{r})$ have been proposed in literature to improve the physical description compared the LDA. The GGA functional that was employed for some initial cell and geometry optimizations in this thesis is PBE, which are the initials of the names of its developers [175].

$$\begin{array}{}\text{(3.11)}& E_{\mathrm{XC}}^{\mathrm{PBE}}[\rho ]=\int F(\rho ,\mathrm{\nabla }\rho ){ϵ}_{\mathrm{XC}}(\rho )\mathrm{d}\mathit{r}.\end{array}$$

This functional uses only fundamental constants as parameters and the dependency on the gradient of $\rho$ is incorporated through the enhancement factor $F(\rho ,\mathrm{\nabla }\rho )$.

While GGA tends to provide a more accurate description of the total energies, energy barriers and structural energy differences when compared to LDA [175], it still lacks an accurate description of the electronic properties of a solid state system. In particular, GGA is known to considerably underestimate the electronic band gap of insulators and semiconductors [176]. Because predicting the band gap as accurately as possible is crucial for analyzing the charge trapping properties of defects in electronic devices as described in Section 2.1.4, more advanced models are needed.

A significant improvement of the description of electronic properties can be achieved by the employment of hybrid functionals, which incorporate calculations of the exact Hartree-Fock exchange $E_{\mathrm{X}}^{\mathrm{HF}}$ into $E_{\mathrm{XC}}$ [177]. Hereby, the exchange-correlation energy is split in a term describing the exchange energy of the non-interacting reference system $E_{\mathrm{X}}^{\mathrm{HF}}$ and an exchange-correlation term of the fully-interacting system $E_{\mathrm{XC}}^{\mathrm{GGA}}$. It is often considered sufficient to only calculate the exchange part exactly and approximate the correlation part $E_{\mathrm{C}}^{\mathrm{PBE}}[\rho ]$, as exchange is by far the dominant component in $E_{\mathrm{XC}}$ [177]. $E_{\mathrm{X}}^{\mathrm{HF}}$ can be calculated as the exchange energy of the Slater determinant formed from the Kohn-Sham orbitals, in a similar fashion as it is done in Hartree-Fock theory [177, 178].

$$\begin{array}{}\text{(3.12)}& E_{\mathrm{X}}^{\mathrm{HF}}=-\frac{1}{2}\sum _{ij}^N\int \int {\varphi }_i({\mathit{r}}_1){\varphi }_j({\mathit{r}}_2)g_{12}{\varphi }_j({\mathit{r}}_1){\varphi }_i({\mathit{r}}_2)\mathrm{d}{\mathit{r}}_1\mathrm{d}{\mathit{r}}_2,\end{array}$$

where $g_{12}$ denotes the Coulomb operator.

The idea behind hybrid functionals is to obtain a partially interacting system, where both $E_{\mathrm{X}}^{\mathrm{HF}}$ and $E_{\mathrm{XC}}^{\mathrm{GGA}}$ depend on the same electron density $\rho$. Determining only a small portion of $E_{\mathrm{X}}$ from exact Hartree-Fock exchange improves the accuracy of $E_{\mathrm{XC}}$ compared to a full Hartree-Fock treatment of exchange [179]. The exact composition of $E_{\mathrm{XC}}$ depends on the respective implementation. For example, the hybrid functional PBE0 is given in its basic form as

$$\begin{array}{}\text{(3.13)}& E_{\mathrm{XC}}^{\mathrm{PBE}0}[\rho ]=\alpha E_{\mathrm{X}}^{\mathrm{HF}}[\rho ]+(1-\alpha )E_{\mathrm{X}}^{\mathrm{PBE}}[\rho ]+E_{\mathrm{C}}^{\mathrm{PBE}}[\rho ],\end{array}$$

with a standard mixing parameter $\alpha =0.25$. In this thesis, a variant of the PBE0 potential was used [180], which employs a truncated Coulomb operator in Eq. (3.12) with a cutoff radius $R_{\mathrm{c}}$ for more efficient computation.

$$\begin{array}{}\text{(3.14)}& g_{12}=\{\begin{array}{ll}\frac{1}{r_{12}},& r_{12}\le R_{\mathrm{c}},\\ 0,& r_{12}>R_{\mathrm{c}}\end{array}\end{array}$$

The mixing parameter $\alpha$ of the hybrid functional can be adjusted to match the experimentally determined band gap of the investigated material [181]. A quantitative comparison between electronic energies calculated with PBE and
PBE0_TC_LRC is presented in Appendix D.

3.1.4 Solving the equation

With all approximations applied and all functionals defined, the task is now to solve Eq. (3.7) as accurately and as efficiently as possible. Since the potentials depend on the electron density, which itself is directly related to the Kohn-Sham orbitals by Eq. (3.4), the equation has to be solved in a self-consistent manner. Starting with an educated initial guess for $\rho$, such as the superposition of atomic charge densities, the initial effective potential $v_{\mathrm{eff}}(\rho )$ is defined. Subsequently, the Kohn-Sham equations including $v_{\mathrm{eff}}(\rho )$ are solved to obtain the single-particle orbitals ${\varphi }_i$ and the according energies ${ϵ}_i$. The new electron density is calculated from these orbitals and the effective potential updated for the new $\rho$. These steps are repeated until the input and output electron densities become equal within a convergence criteria specified for the desired precision.

Over the last decades, several software packages have been developed to tackle this problem in a numerically stable way. Although each of these packages offers its own approach for this challenge, the obtained values and properties of most recent codes only vary slightly [182]. A main difference between the packages lies in their approach to expand the Kohn-Sham orbitals in an appropriate basis set

$$\begin{array}{}\text{(3.15)}& {\varphi }_i=\sum _{\beta }{c}_{\beta }{u}_{\beta }\end{array}$$

which is necessary to allow for efficient numerical calculations.

Two established examples for possible basis sets are plane waves and Gaussian type orbitals (GTOs) [183]. Plane waves satisfy the periodic boundary conditions of bulk materials, offer an orthonormal and complete basis set and are well-behaved in reciprocal space by construction. The size of the basis set only depends on the system size and the cutoff energy. However, since plane waves are delocalized, they are not well suited to describe localized orbitals such as trapped charges at defect sites. Consequently, a large large basis set is required for such investigations which makes the calculations expensive. Additionally, the description of core electrons with plane waves is computationally expensive, which is why the fast varying wave functions near the core region are often replaced with pseudopotentials [184, 185].

GTOs, on the other hand, are centered at the nucleus and provide a compact description of localized orbitals. This is particularly helpful for calculations in amorphous systems, which lack any internal periodicity. Furthermore, GTOs considerably speed up the evaluation of most integrals as their simple mathematical form often allows for an analytical solution. On the downside, GTOs are not orthogonal and can thus suffer from numerical issues such as over-completeness or linear dependency.

The code utilized for this thesis is CP2K [186], which attempts to utilize the benefits of both approaches by combining them into a hybrid basis set called the Gaussian plane wave (GPW) method. The hybrid basis set consists of plane waves and Gaussian functions [187] and is employed in the Quickstep code [188] as implemented in CP2K. In this approach, the GTOs are used for expansion of the wave functions, while the electron density is modeled in a plane wave basis set.

For all calculations in this thesis, core electrons are not explicitly treated but replaced with pseudopotentials to considerably speed up the calculations [185]. Furthermore, to accelerate the computation of $E_{\mathrm{X}}^{\mathrm{HF}}$, the auxiliary density matrix method was employed [189]. CP2K initially solves the Kohn-Sham equations only at the $\mathrm{\Gamma }$ point, corresponding to wave vector $k=0$. This was considered sufficient for all investigations in this thesis, as the Brillouin zone becomes very small in reciprocal space for large supercells and is therefore adequately represented by the $\mathrm{\Gamma }$ point.

In summary, DFT is an established tool to model the electronic density of the ground state of many-body atomic systems, capturing both classical and quantum-mechanical effects. When combined with sophisticated models for approximating the exchange and correlation effects, DFT mostly gives accurate predictions of structural and electronic properties. However, even though various approximations are applied, solving the Schrödinger equation remains complex and computationally demanding and other methods may be required for specific problems.