The RPA and RPA+BSE Dielectric functions

The macroscopic dielectric function ϵM(ω)\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 ϵM\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.


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 Σ=iGW\Sigma=iGW. The screened Coulomb interaction is calculated using the bare Coulomb interaction and the dielectric function ϵ(r,r;ω)=δ(r,r)dr"v(rr")P(r",r;ω)\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:

PRPA(r,r;ω)=n1n2(fn2fn1)ψn2(r)ψn1(r)ψn1(r)ψn2(r)(εn2εn1)ωiη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

ϵM(ω)=limq01ϵ1(q;ω)G=G=0\epsilon_M(\omega)=\lim_{\mathbf q\rightarrow 0}\frac{1}{\epsilon^{-1}(\mathbf q;\omega)|_{\mathbf G=\mathbf G'=0}}

i.e., to determine ϵM(ω)\epsilon_M(\omega), we calculate the dielectric matrix in a particular basis, this matrix is inverted to get ϵ1\epsilon^{-1}, which must then be projected into the plane wave basis. The macroscopic contribution G=G=0\mathbf G=\mathbf G'=0 is taken and finally we assume the optical limit: q0\mathbf q\rightarrow 0. If we assume ϵM(ω)\epsilon_M(\omega) is just ϵG=G=0(q;ω)\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 ϵM(ω)\epsilon_M(\omega) from a modified response function P¯\bar{P} (W. Hanke, 1978, Adv. Phys. 27, 287)

ϵM(ω)=limq0(14πq2P¯G=G=0(q;ω))\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 P¯\bar{P}, where P¯=P+Pv¯P¯\bar{P}=P+P\bar{v}\bar{P}, and v¯\bar{v} is the Coulomb interaction with vG=0(q)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 ϵM(ω)\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:


with K4(r,r¯,r,r¯;ω)=v¯(r,r)W(r,r¯;ω)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 v¯\bar{v} was introduced above, and WW is the screened Coulomb interaction calculated using the RPA polarization. The 4-point RPA polarisation is PRPA4(r,r¯,r,r¯)=PRPA(r,r)δ(rr¯)δ(rr¯)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

P¯(r,r¯,r,r¯)=n1n2n3n4P¯n1n2n3n4ψn1(r)ψn2(r¯)ψn3(r)ψn4(r¯)\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 nin_i contains the band index and k\mathbf k-point, and in the basis of eigenfunctions P¯n1n2n3n4=drdr¯drdr¯ψn1(r)ψn2(r¯)ψn3(r)ψn4(r¯)P¯(r,r¯,r,r¯)\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




KK is constructed from v¯\bar{v} and WW 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 HH then P¯\bar{P} will reduce to the RPA polarisation. In the code, HH 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ω)1(H-\omega)^{-1} in the spectral form. Finally, the macroscopic dielectric function is

ϵM(ω)=1limq04πΩcellNkn1n2kn3n4kλ(fn4k+qfn3k)q^ψn1krψn2k+qA(n1k)(n2k+q),λAλ,(n3k)(n4k+q)Eλω+iηq^ψn3krψn4k+q\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λE_{\lambda}, Aλ\mathbf A_{\lambda} the eigenvalues and eigenfunctions of HH, η\eta a Lorentzian broadening parameter that sets the height and width of the Lorentzian type peaks, and q^\hat{\mathbf q} is a unit vector in the direction of the perturbing field. The matrix elements ψn1rψn2\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 ϵM(ω)\epsilon_M(\omega)).

Edit This Page