Localised electrostatic charges in non-uniform dielectric media: Theory

For examples, see Calculating the formation energies of charged defects.


The purpose of this section is to enable us to calculate the formation energy of charged defects. In the Zhang-Northrup formula, this is given by

\[E^f[X^q] = E[X^q] - E_0 - \sum_i\mu_in_i + q (\epsilon_v + \epsilon_F)\]

In this formula, \(X\) labels the type of defect and \(q\) its charge state, i.e. the net charge contained in some volume surrounding the defect. \(q\) is defined such that \(q=-1\) for an electron. \(E[X^q]\) is the total energy of the sample with the defect, and \(E_0\) the energy of the pristine (bulklike) sample.

These quantities are usually calculated in a supercell approach with periodic boundary conditions. That is, we create the defect we are interested in, and place it in an environment containing many repetitions of the pristine unit cell. In the infinite supercell limit, we can describe the properties of the isolated defect.

When we employ periodic boundary conditions, our calculation includes spurious electrostatic interactions between localised charge distribution of the defect state, and all its periodically repeated images. These long-ranged interactions mean that the convergence with respect to supercell size is slow.

To accelerate this convergence, we employ an electrostatic correction scheme, and write

\[E^f[X^q]_{\mathrm{corrected}} = E^f[X^q]_{\mathrm{uncorrected}} - E_{\mathrm{periodic}} + E_{\mathrm{isolated}} + q\Delta V\]

This assumes that DFT describes the bonding and energies of the system well, but contains errors in the electrostatics. The correction thus consists in subtracting the spurious interactions between the periodic images (the term included in the DFT calculation) and adding in the energy of an isolated charge distribution in the given dielectric environment. Finally, the potentials between the charged and neutral states must be aligned to the same reference.

Implementing this scheme consists of the following steps, which will be described in further detail in the sections below. First, we determine the \(z\)-dependent dielectric function of the 2D layers. Knowing this, we use a model charge distribution to emulate the behaviour of the defect charge state, and find the potential associated with this model distribution by solving the Poisson equation. We do this twice: first with periodic boundary conditions, and then with zero boundary conditions. The electrostatic energy associated with the charge distribution in a given potential can then be found from the usual formula,

\[U = \frac{1}{2}\int_{\Omega} \rho V.\]

Defining the dielectric response of 2D layers

In two dimensions, the bulk dielectric response is poorly defined, and some model must be used.

The approach used here, following Ref [1] is to assume that the dielectric function of the isolated layer is isotropic in plane, and varies only in the \(z\) direction. Additionally, the screening is assumed to follow the density distribution of the system so that we can write

\[\varepsilon^{i}(z) = k^i\cdot n(z) + 1,\]

Where \(i\) varies over “in-plane” and “out-of-plane”, and \(n\) is the in-plane averaged density of the system. The normalization constants, \(k^i\), are chosen such that

\[\begin{split}\frac{1}{L} \int \mathrm{d} z\, \varepsilon^{\parallel}(z) &= \varepsilon^{\parallel}_{\mathrm{DFT}} \\ \frac{1}{L} \int \mathrm{d} z\, \left(\varepsilon^{\perp}(z)\right)^{-1} &= \left(\varepsilon^{\perp}_{\mathrm{DFT}}\right)^{-1}\end{split}\]

Calculating the energy of the periodic images

Once we have the dielectric response, we proceed by solving the poisson equation

\[\nabla \cdot \boldsymbol{\varepsilon}(z) \odot \nabla \phi(\mathbf r) = -\rho(\mathbf r).\]

Here \(\odot\) is the elementwise (Hadamard) product, and \(\boldsymbol{\varepsilon} = (\varepsilon^{\parallel}, \varepsilon^{\parallel}, \varepsilon^{\perp})\). \(\rho(\mathbf r)\) is a model charge distribution that we use to describe the charge distribution of the defect state. For convenience, a Gaussian is often chosen, so that

\[\rho(\mathbf r) = \frac{1}{\left(\sqrt{2\pi}\sigma\right)^3} e^{-(\mathbf r - \mathbf r_0)^2/(2\sigma^2)}.\]

In terms of the coordinate axes, the poisson equation reduces to

\[\varepsilon^{\parallel}(z)\frac{\partial^2 V}{\partial x^2} + \varepsilon^{\parallel}(z)\frac{\partial^2 V}{\partial y^2} + \varepsilon^{\perp}(z)\frac{\partial^2 V}{\partial z^2} + \frac{\partial \varepsilon^{\perp}}{\partial z} \frac{\partial V}{\partial z} = -\rho(\mathbf r)\]

Since we wish to solve this equation with periodic boundary conditions, we Fourier transform the above equation, giving

\[\varepsilon^{\parallel}(G_z) * \left[(G_x^2 + G_y^2)V(\mathbf G)\right] + \varepsilon^{\perp}(G_z) * \left[G_z^2 V(\mathbf G)\right] + \left[G_z \varepsilon^{\perp}(G_z)\right] * \left[ G_z V\mathbf(G)\right] = \rho(\mathbf G),\]

where \(*\) denotes a convolution along the \(z\) axis. Writing out these convolutions, we finally arrive at the expression

\[\rho_{G_x,G_y,G_z} = \sum_{G_z'} \left[\varepsilon^\parallel_{G_z - G_z'}\left(G_x^2 + G_y^2\right) + \varepsilon^{\perp}_{G_z - G_z'}G_zG_z'\right] V_{G_x, G_y, G_z'}.\]

For each value of \((G_x, G_y)\), we can thus calculate the corresponding potential \(V_{G_z}\) through a matrix inversion, and use that to calculate the energy of the model charge distribution.

We can also use the potential \(V_{G_z}\) to calculate the alignment term, \(\Delta V\). We can Fourier transform this to get real-space potential of the model charge distribution. If we have described the electrostatics of the system well, this potential should be similar to the true potential of the defect charge distribution, up to a constant shift. Defining

\[\Delta V(\vec{r}) = V(\vec{r}) - [V^{X^q}_\mathrm{el}(\vec{r}) - V^{0}_\mathrm{el}(\vec{r}) ],\]

We set

Calculating the energy of the isolated system

We start as before, with the Poisson equation, but since we would like to describe the energy of the isolated defect, we do not impose periodic boundary conditions and Fourier transform. Instead, following Ref. [2] we can exploit the in-plane symmetry of the problem and expand \(\phi\) using cylindrical Bessel functions.

\[\phi(\mathbf r) = \int_0^\infty \mathrm{d}k'\, 2qe^{-k'^2\sigma^2/2} \varphi_{k'}(z) J_0(\rho k')\]

Inserting this into the above equation and using the orthogonality relation \(\int \rho\mathrm{d}\rho J_0(\rho k)J_0(\rho k') = \delta(k - k') / k\) we find that \(\varphi_k\) must obey the Poisson equation

\[-\frac{\partial}{\partial z}\left(\varepsilon^{\perp}(z) \frac{\partial \varphi_k(z)}{\partial z}\right) + k^2\varepsilon^{\parallel}(z) \varphi_k(z) = \frac{1}{\sqrt{2\pi}\sigma}e^{-\left(z - z_0\right)^2/\left(2\sigma^2\right)},\]

where \(z_0\) is the center of the gaussian density along the \(z\) direction. The normalization of \(\varphi_k\) defined above was chosen precisely so that the right hand side of this equation is a normalized gaussian along the \(z\) direction.

We solve this equation by separating the response into two components: The bulk response, describing the screening far away from the material, and the remaining \(z\) -dependent response close to the system. We thus define \(\Delta \varepsilon^i(z) = \varepsilon^i(z) - \varepsilon^{i}_{\mathrm{bulk}}\) and the Green’s function of the bulk response \(\hat K = (-\varepsilon^{\perp}_{\mathrm{bulk}} \frac{\partial^2}{\partial z^2} + k^2\varepsilon^{\parallel}_{\mathrm{bulk}})^{-1}\). As an implementation detail, we note that for 2D materials, the bulk response is generally 1.

The equation for \(\varphi_k\) can then be written as

\[\hat{K}^{-1} \varphi_k - \frac{\partial}{\partial z}\left(\Delta\varepsilon^{\perp} \frac{\partial \varphi_k}{\partial z}\right) + k^2\Delta\varepsilon^{\parallel}\varphi_k = \frac{1}{\sqrt{2\pi}\sigma}e^{-\left(z - z_0\right)^2/\left(2\sigma^2\right)}.\]

Only the first term on the left hand side is affected by the boundary conditions on \(\varphi_k\). We can solve this by Fourier transforming along the \(z\) axis and wigner-seitz truncating the Green’s function, which yields the following equation

\[\sum_{G_z} D_{G_zG_z'} \left(\varphi_k\right)_{G_z} = e^{-i G_z'z_0 - G_z'^2\sigma^2/2},\]

with the matrix \(D\) given by

\[\frac{1}{L}D_{G_z'G_z} = \frac{\varepsilon^{\parallel}_{\mathrm{b}}k^2 + \varepsilon^{\perp}_{\mathrm{b}}G_z^2}{1 - e^{-kL/2}\cos(G_zL/2)}\delta_{G_zG_z'} + \Delta\varepsilon^{\parallel}_{G_z - G_z'}k^2 + \varepsilon^{\perp}_{G_z - G_z'}G_zG_z'.\]

Finding \(\varphi_k\) is thus just a simple matrix inversion. Once we have solved the poisson equation, we calculate the total energy.

\[\begin{split}U &= \frac{1}{2} \int_{\Omega} \rho(\mathbf r) \phi(\mathbf r) \\ &= q^2 \int k \mathrm{d}k e^{-k^2 \sigma^2} U_k,\end{split}\]


\[U_k \int \mathrm{d}z\, \varphi_k(z) \frac{1}{\sqrt{2\pi}\sigma}e^{-\left(z - z_0\right)^2/\left(2\sigma^2\right)}\]

Using the solution to the Poisson equation, this reduces to

\[U_k = \sum_{G_z,G_z'} e^{i(G_z - G_z')z_0 - (G_z^2 + G_z'^2)\sigma^2 / 2} \left(D^{-1}\right)_{G_zG_z'},\]

With \(D\) defined as above.