Dielectric screening

The dielectric response function is a central concept in the understanding and calculation of many phenomena in condensed matter systems. It is essential in understanding the response of a material to electromagnetic probes. In fact, in performing GW calculations, we need dielectric function as an input to calculate the screening Coulomb interaction. Moreover, constructing exciton kernels in Bethe-Salpeter equation (BSE) calculations also involves the dielectric function. In BerkeleyGW, we explicitly calculate the dielectric function in random-phase approximation (RPA) with the binary epsilon.flavor.x (The binaries sigma_real.x  is for systems with inversion symmetry and sigma_cplx.x for systems without inversion symmetry). Dynamical screening effects can be included either from the generalized plasmon-pole (GPP) model or explicit finite-frequency calculations.

Random-phase Approximation (RPA)

Random-phase approximation (RPA) to the polarizability means the summation of ring-like Feynman diagrams in calculating the polarizability. [1,2] As a first-order approximation to the dielectric screening effects, RPA is famous for describing the plasmon excitation and correlation energy in homogeneous interacting electron gas, also known as Lindhard function there. [1,2] We calculate the RPA dielectric function using the famous Adler-Wiser equation of polarizability, [3]

\chi_{GG'}(q,\omega)=\sum_{n,n',k} (f(\varepsilon_{n',k+q})-f(\varepsilon_{n,k})) \frac{\langle n, k | e^{-i(q+G)\cdot r}|n',k+q\rangle \langle n',k+q|e^{i(q+G')\cdot r'}|n,k\rangle}{\varepsilon_{n',k+q}-\varepsilon_{n,k}} ,

where f(\varepsilon) is the occupation of each quasi-particle states | n,k\rangle.

The dielectric function is then given by,

\epsilon={\bf 1}-{\bf V}\cdot {\bf \chi},

where {\bf V} is the matrix elements of bare Coulomb interaction and {\bf 1} is the identity matrix.

epsilon is a standalone executable that computes either the static or dynamic RPA polarizability and corresponding inverse dielectric function from input electronic eigenvalues and eigenvectors computed in a suitable meanfield code. In the code, we define the matrix elements M as,[4]

M(G,q,(n,n',k))=M_{nn'}(k,q,G) \cdot \frac{1}{\sqrt{\varepsilon_{n,k+q}-\varepsilon_{n',k}}}.

There is a clear problem in directly computing \epsilon_{00}(q = 0) due to the fact that the Fourier-transformed Coulomb potential V(|q+G|)  diverges as q+G \rightarrow 0 except in the case of box-type truncation schemes. For semiconducting systems, due to orthogonality, the matrix elements themselves go to 0 with the form |M_{nn}'(k, q, G = 0)| \propto |q|. Thus \epsilon (q \rightarrow 0) contains a non-trivial q^2/q^2 limit. The epsilon code has implemented a simpler scheme, however, in which we numerically take the limit as q \rightarrow 0 by evaluating \epsilon_{00}(q_0) at a small but finite q_0 usually taken as approximately \frac{1}{1000}th of the Brillouin zone, in one of the periodic directions. For semiconducting systems, where \epsilon_{00}(q = 0) \rightarrow C, it is sufficient to construct a separate k-grid for the conduction and valence bands shifted by the small vector q_0 in order to compute M'_{nn}(k, q_0, G = 0), where n is a valence and n' a conduction band, and to evaluate the correct limiting q^2/q^2 ratio. For metals, however, intra-band transitions have |M_{nn}'(k, q, G = 0)| \propto C, yielding \epsilon_{00}(q \rightarrow 0) \propto C '/q^2 . In this case, the two-k-grid treatment is insufficient, because the proportionality coefficient C' depends sensitively on the density of states (DOS) at the Fermi energy. Therefore a k-grid sampling of the same spacing as q_0 is required, although fewer conduction bands are necessary in the sum since \epsilon(q \rightarrow 0) is dominated by intra-band transitions. Thus we typically calculate \epsilon(q \rightarrow 0) using a single fine wavefunction grid by using the smallest q consistent with the grid.

The resultant dielectric functions are stored in eps0mat.h5 (\epsilon(q_0\rightarrow 0)) and epsmat.h5 \epsilon(q_1) files.

Plasmon-pole Model

To extend the dielectric matrix to finite frequencies, we propose a generalized plasmon-pole (GPP) model to include dynamical dielectric screening effects. Since the evaluation of the self-energy operator generally involves a sum over frequencies in the screened Coulomb interaction, the fine details of the frequency dependence of the dielectric function should not be important. Thus for the present purposes, the GPP model is valid if it represents the average features of \mathrm{Re} \, \epsilon^{-1}(\omega) for all the important elements of the dielectric matrix. The form of the dynamical dielectric function is given by, [3]

\mathrm{Im} \epsilon^{-1}_{GG'}(q,\omega)=A_{GG'}(q)\left [ \delta (\omega - \tilde{\omega}_{GG'}(q)) - \delta (\omega + \tilde{\omega}_{GG'}(q)) \right ],


\mathrm{Re} \epsilon^{-1}_{GG'}(q,\omega) = 1 + \frac{\Omega^2_{GG'}(q)}{\omega^2 - \tilde{\omega}^2_{GG'}(q)}.

And the effective bare plasma frequency \Omega_{GG'}(q) is defined as,

\Omega^2_{GG'}(q)=\omega^2_p \frac{(q+G)\cdot(q+G')}{|q+G|^2} \frac{\rho(G-G')}{\rho(0)},

where \omega_p is the conventional plasma frequency proportional to the average electron density.

The mode frequency \tilde{\omega}_{GG'}(q) and the amplitude A are determined by the Kramer-Kronig relation and a generalized f-sum rule relating the imaginary part of the many-body dielectric matrix to the plasma frequency and the crystalline charge density. Please refer to Eqn. 28&29 in Ref. 3.


[1] A. L. Fetter, and J. D. Walecka. Quantum theory of many-particle systems. Courier Corporation, 2012.

[2] G. D. Mahan,  Many-particle physics. Springer Science & Business Media, 2013.

[3] M. S. Hybertsen and S. G. Louie, Phys. Rev. B 34, 5390 (1986).

[4]  J. Deslippe, G. Samsonidze, D. A. Strubbe, M. Jain, M. L. Cohen, and S. G. Louie, Comput. Phys. Commun. 183, 1269 (2012), 1111.4429.