Including ladder diagrams in W

The screened Coulomb interaction W 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. We correct the full polarization matrix at all k points in the IBZ and this is used to calculate the screened Coulomb interaction and hence the correlation part of the self-energy. A tutorial using Questaal to include ladder diagrams in W 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

Questaal performs QSGW calculations using the script lmgwsc, which in turn calls further scripts; see the output from lmgwsc to see the set of scripts run and their output. The screened Coulomb interaction is calculated in the program hx0fp0 at the level of the random phase approximation (RPA) and this is used to determine the screened part of the self-energy. Local field effects are included in hx0fp0, however, excitonic effects are not. Excitonic effects, or ladder diagrams/vertex contributions, in W can have a huge impact on the electronic structure, and the neglection of these perhaps contribute the most to the bandgap overestimation of QSGW.

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}$$

where the spin, k and band index are subsumed into the indices $n_1$ and $n_2$.

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: $$P(1,2)=P_{RPA}(1,2)-\int P_{RPA}(1,1,3,4)W(3,4;\omega=0)P(3,4,2,2){\mathrm d}(34)$$

with $W$ the screened Coulomb interaction calculated using the RPA polarization. Note that the index (1,2) means $(\mathbf r_1 t_1\sigma_1,\mathbf r_2 t_2 \sigma_2)$. The 4-point polarisation is such that $P(\mathbf r,\mathbf r')=P(\mathbf r,\mathbf r,{\mathbf r}',\mathbf r')$. The Bethe-Salpeter equation can be derived by combining the Hedin’s equations for the polarization and vertex and assuming the GW approximation for the interaction kernel, $\delta \Sigma/\delta G=iW$. We express the functions in transition space, such that $$P(\mathbf r,\bar{\mathbf r},\mathbf r',\bar{\mathbf r}')=\sum_{n_1n_2n_3n_4}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, again, $n_i$ contains the band and spin indices and $\mathbf k$-point, and in the basis of eigenfunctions $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}')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 $$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})W_{n_1n_2n_3n_4}(\omega)$$

where $W_{n_1n_2n_3n_4}$ is calculated from the RPA $W$ in the mixed product basis used in questaal.
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: $$P_{I,J}(\mathbf q,\omega)=\sum_{n_1n_2n_3n_4}\sum_{\sigma}\sum_{\mathbf k,\mathbf k'}\langle \widetilde{M}_I^{\mathbf q}\psi_{n_1\mathbf k}^{\sigma}|\psi_{n_2\mathbf k+\mathbf q}^{\sigma}\rangle\sum_{\lambda}\frac{A_{n_1n_2\mathbf k,\lambda}^{\sigma}(\mathbf q)A_{\lambda ,n_3n_4\mathbf k'}^{\sigma\dagger}(\mathbf q)}{E_{\lambda}^{\sigma}-\omega\pm i\eta} \langle \psi_{n_1\mathbf k'+\mathbf q}^{\sigma}|\psi_{n_3\mathbf k'}^{\sigma}\widetilde{M}_J^{\mathbf q}\rangle,$$

where $A^{\sigma}(\mathbf q)$ and $E^{\sigma}(\mathbf q)$ are the eigenvectors and eigenvalues of the matrix $H$.

The broadening $\eta$ determines the width and height of the peaks; in Questaal, the default is 0.01 Hartree. Finally, the number of states considered in the BSE must be considered. The matrix $H$ has the shape $N\times N$, where $N$ is the product of the number of spin channels, the number of $k$-points in the FBZ, the number of valence states, and the number of conduction states. Due to memory limitations, we cannot consider transitions between all the states. We pick an energy window about the Fermi level and consider states corresponding to transitions within this window. Transitions not considered in the BSE will be included in the polarization at the level of the RPA.

Finally, once we have the polarization, we can calulate the screened Coulomb interaction, and hence the screened part of the self-energy $$W_{\mu,\nu}(\mathbf q,\omega)=\sqrt{V_{\mu}(\mathbf q)}\left[1-V(\mathbf q)P(\mathbf q,\omega)\right]_{\mu,\nu}^{-1}\sqrt{V_{\nu}(\mathbf q)}$$