The RPA and RPA+BSE Dielectric functions

The macroscopic dielectric function $\epsilon_M(\omega)$ is usually calculated at the level of the random phase approximation (RPA). Here we will go through the background of how we can go beyond the RPA using the Bethe-Salpeter equation (BSE) for the polarization. A tutorial using Questaal to calculate $\epsilon_M$ at BSE level is included on this page. We will introduce quantities and methods here that are used in the program, and relate to these in the tutorial given.

Introduction

The RPA macroscopic dielectric function can be calculated in lmf using the OPTICS category in the ctrl file; see this page for a detailed description. Whilst the dielectric function in the form produced by lmf is useful for analysis there are two crucial effects missing. The first is the effect of local fields, as we will discuss below; and the second is the inclusion of excitonic (electron-hole interaction) effects. The latter can have a huge impact on both the quantitative and qualitative description of $\epsilon(\omega)$. The difference between including and neglectting the electron-hole interaction can be demonstrated when comparing the absorption calculated, Im{$\epsilon(\omega)$}, with that determined experimentally. Excitonic and local-field effects can be included through use of the Bethe-Salpeter equation for the polarisation.

Dielectric Function

In the GW formalism, the self energy is calculated from the Green’s function and screened Coulomb interaction $\Sigma=iGW$. The screened Coulomb interaction is calculated using the bare Coulomb interaction and the dielectric function $\epsilon(\mathbf r,\mathbf r';\omega)=\delta(\mathbf r,\mathbf r')-\int d\mathbf r" v(\mathbf r-\mathbf r")P(\mathbf r",\mathbf r';\omega)$. The RPA (or independent particle) polarisation assumes a sum over independent transitions: $$P_{RPA}(\mathbf r,\mathbf r';\omega)=\sum_{n_1n_2}(f_{n_2}-f_{n_1})\frac{\psi_{n_2}^{*}(\mathbf r)\psi_{n_1}(\mathbf r)\psi_{n_1}^{*}(\mathbf r')\psi_{n_2}(\mathbf r')}{(\varepsilon_{n_2}-\varepsilon_{n_1})-\omega-i\eta}$$

The role of local fields

Now, the macroscopic dielectric function is calculated from $$\epsilon_M(\omega)=\lim_{\mathbf q\rightarrow 0}\frac{1}{\epsilon^{-1}(\mathbf q;\omega)|_{\mathbf G=\mathbf G'=0}}$$

i.e., to determine $\epsilon_M(\omega)$, we calculate the dielectric matrix in a particular basis, this matrix is inverted to get $\epsilon^{-1}$, which must then be projected into the plane wave basis. The macroscopic contribution $\mathbf G=\mathbf G'=0$ is taken and finally we assume the optical limit: $\mathbf q\rightarrow 0$. If we assume $\epsilon_M(\omega)$ is just $\epsilon_{\mathbf G=\mathbf G'=0}(\mathbf q;\omega)$ then we effectively neglect local field effects. To include local field effects we can, however, calculate $\epsilon_M(\omega)$ from a modified response function $\bar{P}$ (W. Hanke, 1978, Adv. Phys. 27, 287) $$\epsilon_M(\omega)=\lim_{\mathbf q\rightarrow 0}(1-\frac{4\pi}{|\mathbf q|^2}\bar{P}_{\mathbf G=\mathbf G'=0}(\mathbf q;\omega))$$

which means that we do not have to calculate the whole dielectric matrix and then invert it to extract just one single matrix element. Instead, we just have to calculate one single matrix element of the modified response function $\bar{P}$, where $\bar{P}=P+P\bar{v}\bar{P}$, and $\bar{v}$ is the Coulomb interaction with $v_{\mathbf G=0}(\mathbf q)$ set to zero. Local field effects are included in the generation of the polarisation in the GW(QSGW) calculation and will also be included when we calculate $\epsilon_M(\omega)$ using the Bethe-Salpeter equation. Comparing the effect of local field corrections on the imaginary part of the dielectric function is presented in, for example, Fig. 11 Phys. Rev. B76, 165106 (2007).

Ladder Diagrams

Excitonic effects (or ladder diagrams in the language of Feynman diagrams) can be included though the Bethe-Salpeter equation for the 4-point polarisation: $$\bar{P}^4=P_{RPA}^4+P_{RPA}^4K^4\bar{P}^4$$

with $K^4(\mathbf r,\bar{\mathbf r},\mathbf r',\bar{\mathbf r}';\omega)=\bar{v}(\mathbf r,\mathbf r')-W(\mathbf r,\bar{\mathbf r};\omega)$, where $\bar{v}$ was introduced above, and $W$ is the screened Coulomb interaction calculated using the RPA polarization. The 4-point RPA polarisation is $P_{RPA}^4(\mathbf r,\bar{\mathbf r},\mathbf r',\bar{\mathbf r}')=P_{RPA}(\mathbf r,{\mathbf r}')\delta(\mathbf r-\bar{\mathbf r})\delta(\mathbf r'-\bar{\mathbf r}')$. We express $$\bar{P}(\mathbf r,\bar{\mathbf r},\mathbf r',\bar{\mathbf r}')=\sum_{n_1n_2n_3n_4}\bar{P}_{n_1n_2n_3n_4}\psi_{n_1}(\mathbf r)\psi_{n_2}^{*}(\bar{\mathbf r})\psi_{n_3}^{*}(\mathbf r')\psi_{n_4}(\bar{\mathbf r}')$$

where $n_i$ contains the band index and $\mathbf k$-point, and in the basis of eigenfunctions $\bar{P}_{n_1n_2n_3n_4}=\int d\mathbf r d\bar{\mathbf r} d\mathbf r' d\bar{\mathbf r}' \psi_{n_1}^{*}(\mathbf r)\psi_{n_2}(\bar{\mathbf r})\psi_{n_3}(\mathbf r')\psi_{n_4}^{*}(\bar{\mathbf r}')\bar{P}(\mathbf r,\bar{\mathbf r},\mathbf r',\bar{\mathbf r}')$ (easily shown using the completeness relation for the eigenfunctions). The polarization in the eigenfunction basis can then be expressed as $$\bar{P}_{n_1n_2n_3n_4}(\omega)=(f_{n_4}-f_{n_3})(H-\omega)^{-1}_{n_1n_2n_3n_4}$$

with $$H_{n_1n_2n_3n_4}(\omega)=(\varepsilon_{n_2}-\varepsilon_{n_1})\delta_{n_1n_3}\delta_{n_2n_4}-(f_{n_4}-f_{n_3})K_{n_1n_2n_3n_4}(\omega)$$

$K$ is constructed from $\bar{v}$ and $W$ in the mixed product basis that is used in the QSGW routine to calculate the self energy. If we take just the diagonal part of $H$ then $\bar{P}$ will reduce to the RPA polarisation. In the code, $H$ is constructed assuming a static interaction kernel and adopting the Tamm-Dancoff approximation (Onida. Rev. Mod. Phys. 74, 601, 2002) and this is then diagonalized to represent $(H-\omega)^{-1}$ in the spectral form. Finally, the macroscopic dielectric function is $$\epsilon_M(\omega)=1-\lim_{\mathbf q \rightarrow 0}\frac{4\pi}{\Omega_{cell}N_k} \sum_{n_1n_2\mathbf k n_3n_4\mathbf k'}\sum_{\lambda}(f_{n_4\mathbf k'+\mathbf q}-f_{n_3\mathbf k'})\hat{\mathbf q}\cdot\langle \psi_{n_1\mathbf k}|\mathbf r|\psi_{n_2\mathbf k +\mathbf q}\rangle^{*}\frac{A_{(n_1\mathbf k)(n_2\mathbf k+\mathbf q),\lambda}A_{\lambda,(n_3\mathbf k')(n_4\mathbf k'+\mathbf q)}^{\dagger}}{E_{\lambda}-\omega+ i\eta}\hat{\mathbf q}\cdot\langle \psi_{n_3\mathbf k'}|\mathbf r|\psi_{n_4\mathbf k'+\mathbf q}\rangle$$

with $E_{\lambda}$, $\mathbf A_{\lambda}$ the eigenvalues and eigenfunctions of $H$, $\eta$ a Lorentzian broadening parameter that sets the height and width of the Lorentzian type peaks, and $\hat{\mathbf q}$ is a unit vector in the direction of the perturbing field. The matrix elements $\langle \psi_{n_1}|\mathbf r|\psi_{n_2}\rangle$ are determined using lmf (see the tutorial as to how these are determined for the BSE $\epsilon_M(\omega)$).