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 electron-hole pairs, 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 $\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, the final result can be linked to EELS 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.

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

Time-propagation algorithm

VASP can employ a time-propagation algorithm to evaluate $\epsilon_{\mathrm M}$ using different exchange-correlation functionals to model the electron-hole interaction. This calculation requires a previous step where the ground state is computed including extra empty states. Afterwards, VASP will integrate a time-dependent equation, using the functional specified in the INCAR.

The following INCAR file can be used as an example, where a range-separated hybrid functional (HSE06) is used. The system in case is the same (i.e. bulk silicon), and the time-propagation is done using 2000 time-steps. Notice that the same number of bands (NBANDS) is used here as in the previous step, but the total number of virtual states (NBANDSV is much lower. As it is explained in the page related to the time-propagation algorithm, only a few virtual states are needed in this calculation, and, in fact, this method does not propagate them in time.

SYSTEM = Si
ALGO = TIMEEV
NBANDS = 48
NBANDSO = 4
NBANDSV = 8
ISMEAR = 0 ; SIGMA = 0.05
IEPSILON = 1 
NELM = 2000
CSHIFT = 0.1
OMEGAMAX = 20
LHARTREE = .TRUE.
LADDER = .TRUE.
LFXC = .FALSE.
LHFCALC = .TRUE.
LMODELHF = .TRUE.
AEXX = 0.088
HFSCREEN = 1.26

Once the calculation is over, the macroscopic dielectric function can be extracted from the vasprun.xml file, using a very similar script to the one previously employed

awk 'BEGIN{i=0} /<dielectricfunction comment="time-propagation">/,\
                /<\/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.TIMEEV.dat

Casida equation

Another method that can employ hybrid functionals to evaluate $\epsilon_{\mathrm M}$ is the use of Casida's equation, set with ALGO=TIMEEV in the INCAR. Here an eigenvalue equation is employed instead of a time-dependent differential equation. If the following INCAR is used as an example

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 using

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.TDHF.dat

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. Otherwise results will be unphysical!