Electron-energy-loss spectrum

One of the many ways with which is possible to probe neutral excitations in a material is by injecting electrons into the sample. These are called electron-energy-loss spectroscopy experiments, where the incoming electron can create bound electron-hole pairs (i.e. excitons), plasmons, or even higher-order multi-pair excitations.

The incoming electron acts as an external potential, $V_\mathrm{ext}(\mathbf r', t')$, which induces a charge density in the material, $\rho_\mathrm{ind}(\mathbf r, t)$. Within linear-response theory these two quantities can be related by the reducible polarisability function, $\chi$, via a Green-Kubo relation

$\rho_\mathrm{ind}(\mathbf r, t) = \int \mathrm d^3r'\mathrm d t \chi(\mathbf r, t,\mathbf r', t')V_\mathrm{ext}(\mathbf r', t').$

If the external potential is taken as proportional to a plane-wave of momentum $\mathbf q$, then the electron energy-loss spectrum (EELS) can be taken from the imaginary part of the inverse dielectric function, $\epsilon^{-1}(\mathbf q,\omega)$, since $\epsilon^{-1} = 1 + v\chi$

$\mathrm{EELS}(\mathbf q,\omega) = -\mathrm{Im}\left[\epsilon^{-1}(\mathbf q,\omega)\right].$

Inclusion of local fields

In general, the microscopic quantity $\epsilon^{-1}$ is a function of two coordinates, i.e. $\epsilon^{-1} := \epsilon^{-1}(\mathbf r , \mathbf r', \omega)$. This has important consequences on inhomogeneous systems where a homogeneous, constant, external electric field can induce fluctuations at the interatomic scale, and thus create microscopic fields. A direct consequence from the inhomogeneous character of the system is that in reciprocal space $\epsilon$ has to be written as $\epsilon_{\mathbf G, \mathbf G'}(\mathbf q, \omega)$, where $\mathbf G$ is a reciprocal lattice vector. The microscopic fields are then the $\mathbf G \neq 0$ components of the tensor.

From $\epsilon^{-1} = 1 + v\chi$ it is possible to see that a problem arises when $\mathbf q \to 0$, i.e. the optical limit. In reciprocal space this equation becomes

$\epsilon^{-1}_{\mathbf G, \mathbf G'}(\mathbf q, \omega) = \delta_{\mathbf G, \mathbf G'} + \frac{4\pi}{|\mathbf q + \mathbf G|^2}\chi_{\mathbf G, \mathbf G'}(\mathbf q, \omega)$

where $v(\mathbf q + \mathbf G) = 4\pi/|\mathbf q + \mathbf G|^2$ is the Coulomb potential. At $\mathbf q= 0$, all components without microscopic fields are divergent. To circumvent this issue, the evaluation of $\epsilon^{-1}$ is replaced the Coulomb potential with

$\bar v(\mathbf q + \mathbf G) = \left\{ \begin{array}{ll} 0, & \mathbf G=0 \\ 4\pi/|\mathbf q + \mathbf G|^2, & \mathbf G\neq 0 \end{array} \right.,$

leaving the $v_0\chi_{00}$ component to be dealt with separately and then added at the end.

Micro-macro connection and relation to measured quantities

It is important to note that the actual measured quantity, $\epsilon_\mathrm{M}(\mathbf q, \omega)$, does not depend on the microscopic fields. To connect both the microscopic and macroscopic quantities, an averaging procedure is taken out, so that

$\epsilon_M(\mathbf q,\omega) = \frac{1}{\epsilon_{\mathbf G = 0, \mathbf G'=0}^{-1}(\mathbf q,\omega)}.$

Since VASP computes the macroscopic function, EELS can be extracted from the final result via

$\mathrm{EELS}(\omega) = -\mathrm{Im}\left[\frac{1}{\epsilon_M(\mathbf q,\omega)}\right]$

Note that the inclusion of local fields and the connection to the macroscopic observable must be considered regardless of the level of approximation to the polarisability function, $\chi$. Within VASP, this can be done at several levels of approximation, which are discussed in the next section.

Computing EELS with VASP

EELS from density functional theory (DFT)

The simplest calculation that yields a the macroscopic dielectric function is a ground state calculation using DFT, with the tags NBANDS and LOPTICS. Using bulk silicon as an example with the following INCAR

SYSTEM = Si
NBANDS = 48
ISMEAR = 0 ; SIGMA = 0.1
ALGO = N
LOPTICS = .TRUE.
CSHIFT = 0.4

the macroscopic dielectric function can be extracted from the vasprun.xml file at the end of the calculation by first running the following script

awk 'BEGIN{i=0} /<dielectricfunction comment="density-density">/,\
                /<\/dielectricfunction>/ \
                 {if ($1=="<r>") {a[i]=$2 ; b[i]=($3+$4+$5)/3 ; c[i]=$4 ; d[i]=$5 ; i=i+1}} \
     END{for (j=0;j<i/2;j++) print a[j],b[j],b[j+i/2]}' vasprun.xml > dielectric_function.DFT.dat

which writes the trace of the tensor to a file called dielectric_function.DFT.dat. Plotting the EEL spectrum can be done with visualisation software like gnuplot, using the instruction

p 'dielectric_function.DFT.dat' u 1:($2/($2**2+$3**2)) w l lc rgb 'black' t 'DFT'

Imaginary part of the inverse macroscopic dielectric function for bulk Si from DFT

The main feature in the EEL spectrum is peak at close to 17.5 eV, which would correspond to the plasmon energy in bulk Si.

The results form this DFT calculation should be analysed with the following considerations in mind. Firstly, this calculation is performed only for $\mathbf q = 0$, with no local fields taken into account. Because of that, the expression for $\epsilon_{\mathrm M}$ is often called the independent particle random-phase approximation dielectric function. Also, there are no interactions between electrons and holes taken into account. To account for both local-field effects and the electron-hole interaction, approximations beyond DFT must be taken into account.

Before continuing, it is important to remark that plotting the EELS can also be done using py4vasp as python module. Taking the next script as an example, where the vaspout.h5 had its name changed to vaspout.DFT.h5 in order to preserve the information for latter plots,

import py4vasp
import numpy as np
import matplotlib.pyplot as plt

calc = py4vasp.Calculation.from_file('vaspout.DFT.h5')

energies = calc.dielectric_function.read('')['energies']
epsilon  = calc.dielectric_function.read('')['dielectric_function'][1][1][:]

eels = np.imag(epsilon)/(np.real(epsilon)**2 + np.imag(epsilon)**2)

plt.plot(energies, eels, label='DFT')
plt.ylabel(r'Im[$\epsilon^{-1}_M$]')
plt.xlabel('Energy loss [eV]')
plt.legend()
plt.show()

the EEL spectrum can be easily visualised

EEL spectrum from DFT plotted with py4vasp.

It is important to note that for the case of bulk Silicon, since the material is centro-symmetric, only one component of the 3x3 tensor. For non-symmetric systems, it is better to take the trace of the dielectric function.

EELS from time-dependent density functional theory (TDDFT)

VASP can compute the macroscopic dielectric function form TDDFT calculations using different functionals. Such calculations not only account for the electron-hole interaction beyond the independent particle level, but also include local-fields and can be performed at points in the Brillouin zone other than $\Gamma$.

The evaluation of $\epsilon_{\mathrm M}$ is performed via Casida's equation, set with ALGO=TDHF in the INCAR. The following INCAR will be used as an example, where the electron-hole interaction is included using a model dielectric function.

SYSTEM    = Si
ISMEAR    = 0 
SIGMA     = 0.1
NBANDS    = 48     
ALGO      = TDHF
CSHIFT    = 0.4
IBSE      = 0
NBANDSO   = 4       ! number of occupied bands
NBANDSV   = 8       ! number of unoccupied bands
LHARTREE  = .TRUE.
LADDER    = .TRUE.
LFXC      = .FALSE.
LMODELHF  = .TRUE. 
AEXX      = 0.083
HFSCREEN  = 1.22

at the end of the calculation the dielectric function can be extracted from the vasprum.xml file. Results are shown at the end of the next subsection, in order to compare the difference between the TDDFT and the many-body perturbation theory approaches.

EELS from many-body perturbation theory (MBPT)

The addition of electron-hole interactions to the dielectric function can be done at MBPT level, with ALGO=BSE. The caveat is that VASP requires a previous GW calculation in order to generate the WFULLxxxx.tmp files, where the dielectric screening is stored.

Like in the TDDFT case, users can select how to build $\epsilon_\mathrm{M}$ via the tag IBSE. In the INCAR file used as an example below, the exact diagonalisation solver is employed

SYSTEM    = Si
ISMEAR    = 0 
SIGMA     = 0.1
NBANDS    = 48     
ALGO      = BSE
CSHIFT    = 0.4
IBSE      = 2
NBANDSO   = 4       ! number of occupied bands
NBANDSV   = 8       ! number of unoccupied bands
LHARTREE  = .TRUE.
LADDER    = .TRUE.

and the same number of virtual and occupied states are used, but now no information is given about any parameters with which to model the electron-hole interaction.

Extracting the data at the end of the calculation can be done for both MBPT and TDDFT using the following script

awk 'BEGIN{i=0} /<dielectricfunction>/,\
                /<\/dielectricfunction>/ \
                 {if ($1=="<r>") {a[i]=$2 ; b[i]=($3+$4+$5)/3 ; c[i]=$4 ; d[i]=$5 ; i=i+1}} \
     END{for (j=0;j<i/2;j++) print a[j],b[j],b[j+i/2]}' vasprun.xml > dielectric_function.BSE.dat

Mind: If both TDDFT and MBPT calculations are run in the same directory, VASP will overwrite the vasprun.xml and vaspout.h5 files, making it impossible to compare both calculations. Please be sure to either change the name of this file after the calculation is done, and before a new one begins, or perform both calculations in different directories.

Imaginary part of the inverse macroscopic dielectric function for bulk Si from both BSE and TDHF at $\Gamma$-point.

Plotting both results can be done in gnuplot using

p 'dielectric_function.BSE.dat' u 1:($2/($2**2+$3**2)) w l lc rgb 'black' t 'BSE', \
'dielectric_function.TDHF.dat' u 1:($2/($2**2+$3**2)) w l lc rgb 'red' t 'TDHF'

which results in the following figure. Both spectra look very similar, the differences coming from the treatment of the electron-hole interaction in both TDDFT and BSE case. By varying the range separation parameters, TDDFT results can be brought up to match almost exactly those from BSE, but this will require more computing time, so it is advisable to test the best values of AEXX and HFSCREEN in lower settings, away from convergence.

However, while the use of a model (diagonal) dielectric function is well suited for bulk Silicon, for systems with lower symmetry the off-diagonal elements that are ignored here will be important, and can lead to larger discrepancies in the results from TDDFT and MBPT.

To plot the EEL spectrum with py4vasp, it should be noted that the type of calculation has to be specified when reading the dielectric function. It is simply a matter of editing the following lines

calc_bse  = py4vasp.Calculation.from_file('vaspout.BSE.h5')
calc_tdhf = py4vasp.Calculation.from_file('vaspout.TDHF.h5')

epsilon_bse  = calc_bse.dielectric_function.read('BSE')['dielectric_function']
epsilon_tdhf = calc_tdhf.dielectric_function.read('TDHF')['dielectric_function']

which when used with the rest of the script will produce the following plot for EELS.

EELS plot using py4vasp for bulk Si from both BSE and TDHF at $\Gamma$-point.

Calculations at finite momentum

Warning: With VASP, finite momentum calculations at TDDFT or BSE level, i.e. ALGO=TDHF or BSE, must always use ANTIRES=2 regardless of the solver, functional, or approximation used for the electron-hole interaction. Otherwise results will be unphysical!

The macroscopic dielectric function can be evaluated at other points in the Brillouin zone. First, it is important to check the point list in the OUTCAR or IBZKPT from the ground-state calculation. In the OUTCAR the points are listed in this section

 Subroutine IBZKPT returns following result:
 ===========================================

 Found     16 irreducible k-points:

 Following reciprocal coordinates:
            Coordinates               Weight
  0.000000  0.000000  0.000000      1.000000
  0.166667  0.000000  0.000000      8.000000
  0.333333  0.000000  0.000000      8.000000
...

with the first three entries being the reduced coordinates, the fourth entry being the weight of the respective point. The same information is also repeated in the IBZKPT file, where the points are listed as

Automatically generated mesh
      16
Reciprocal lattice
    0.00000000000000    0.00000000000000    0.00000000000000             1
    0.16666666666667    0.00000000000000    0.00000000000000             8
    0.33333333333334    0.00000000000000    0.00000000000000             8
...

The point at which $\epsilon_M(\mathbf q,\omega)$ is going go be computed is then selected in the INCAR file with the tag KPOINT_BSE. The syntax is

 KPOINT_BSE = index_of_k-point  n1 n2 n3

where the first entry is the index of the point on the list found in the OUTCAR or IBZKPT files, and $n_i,\, i=1,2,3$ optional integer arguments. If present, these three indices can be used the evaluate the dielectric function at a k point outside of the first Brillouin zone corresponding to

$\bold{k} + n_{1} \bold{b}_{1}+ n_{2} \bold{b}_{2} + n_{3} \bold{b}_{3}.$

Still using bulk silicon as an example, the following INCAR evaluates $\epsilon_M(\mathbf q,\omega)$ at the second k-point on the list

SYSTEM = Si
NBANDS = 48
NBANDSO = 4 ; NBANDSV = 8
ISMEAR = 0 ; SIGMA = 0.1
EDIFF = 1E-8
ALGO = BSE
LADDER = .TRUE.
LHARTREE = .TRUE.
ANTIRES = 2
CSHIFT = 0.4
IBSE = 2
KPOINT_BSE = 2
OMEGAMAX = 30

Extracting the dielectric function from the vasprun.xml file is done in a similar way as $\Gamma$-point calculations, but with a small caveat that the dielectric function now only has one single component. Using the script

awk 'BEGIN{i=0} /<dielectricfunction>/,\
                /<\/dielectricfunction>/ \
                 {if ($1=="<r>") {a[i]=$2 ; b[i]=$3 ; c[i]=$4 ; d[i]=$5 ; i=i+1}} \
     END{for (j=0;j<i/2;j++) print a[j],b[j],b[j+i/2]}' vasprun.xml > dielectric_function.TDHF.dat

and plotting the EELS result can also be done with gnuplot.

Imaginary part of the inverse macroscopic dielectric function for bulk Si from both BSE and TDHF for KPOINT_BSE=2. In this example it corresponds to the coordinates $(0.5,0.0,0.0)$ in reciprocal space.

Much like calculations at $\Gamma$, results from model TDDFT are very close to those from MBPT. The reduction in intensity for spectra away from $\Gamma$ is normal, due to the fact that matrix elements are now away from the dipole approximation.

To plot the same results with py4vasp, the new shape of the dielectric function has to be taken into account, and so no indexing is added to the epslion array

import py4vasp
import numpy as np
import matplotlib.pyplot as plt

calc_bse  = py4vasp.Calculation.from_file('vaspout.BSE.h5')
calc_tdhf = py4vasp.Calculation.from_file('vaspout.TDHF.h5')

energies = calc_bse.dielectric_function.read('BSE')['energies']
epsilon_bse  = calc_bse.dielectric_function.read('BSE')['dielectric_function']
epsilon_tdhf = calc_tdhf.dielectric_function.read('TDHF')['dielectric_function']

eels_bse  = np.imag(epsilon_bse )/(np.real(epsilon_bse )**2 + np.imag(epsilon_bse )**2)
eels_tdhf = np.imag(epsilon_tdhf)/(np.real(epsilon_tdhf)**2 + np.imag(epsilon_tdhf)**2)

plt.plot(energies, eels_bse,  '', label='BSE')
plt.plot(energies, eels_tdhf, '', label='TDHF')
plt.ylabel(r'Im[$\epsilon^{-1}_M$]')
plt.xlabel('Energy loss [eV]')
plt.legend()
plt.show()

which will generate the following plot

EELS plot for KPOINT_BSE=2 using py4vasp.