Time-evolution algorithm

The macroscopic dielectric function, $\epsilon_{ij}(\omega)$, measures how a given dielectric medium reacts when subject to an external electric field. From $\epsilon_{ij}(\omega)$ one can extract several optical properties such as absorption, optical conductivity, reflectance. However, it is important that the interacting electrons and holes are taken into account. This makes the evaluation of the macroscopic dielectric function more involved, since it goes beyond the single-particle level, working either at the Bethe-Salpeter or time-dependent density-functional theory level.

Within VASP, users can select two different methods for how $\epsilon_{ij}(\omega)$ is computed. The first is based on the eigendecomposition of the electron-hole Hamiltonian, $H^\mathrm{exc}$. It allows for the evaluation of $\epsilon_{ij}(\omega)$ by first obtaining the eigenvalues and eigenvectors of $H^\mathrm{exc}$ and it is based on the Bethe-Salpeter equation or the Casida equation. The second method transforms the mathematical expression of $\epsilon_{ij}(\omega)$ into a time-dependent integral. By propagating in time the dipolar moments and then applying a Fourier transform, it can compute $\epsilon_{ij}(\omega)$.

The advantage of the later method in comparison to the former is related to the cost, with the time-dependent integral being $O(N^2)$, while the eigendecomposition has a cost of $O(N^3)$, where $N$ is the rank of $H^\mathrm{exc}$. This means that for very large numbers of bands or k-points, the time-dependent formalism is cheaper than the eigendecomposition method.

Below we present a brief description of the method, from its theoretical support to how calculations should be performed, with the relevant approximations needed in the two-particle Hamiltonian.

The macroscopic-dielectric function as a time-dependent integral

The starting point is that one can re-write $\epsilon_{ij}(\omega)$ as a time-dependent integral. It starts from its expression, given by

$\epsilon^M(\omega)=1+\frac{4 \pi}{\Omega_0} \sum_\lambda\left|\sum_{c v \mathbf{k}} \mu_{c v \mathbf{k}} A_{c v \mathbf{k}}^\lambda\right|^2\left[\frac{1}{\omega+E_\lambda+\mathrm{i} \eta}-\frac{1}{\omega-E_\lambda+\mathrm{i} \eta}\right]$,

where $\mu_{v c \mathbf{k}}^j=\frac{\left\langle c \mathbf{k}\left|v_j\right| v \mathbf{k}\right\rangle}{\varepsilon_c(\mathbf{k})-\varepsilon_v(\mathbf{k})}$ is the dipolar moment associated to the the conduction $c$ and valence band $v$, and k-point $k$, $\lambda$ is the index of the eigenstate of $H^\mathrm{exc}$, with $A^\lambda$ and $E_\lambda$ being the associated eigenvector and eigenvalue. This equation can be brought into operational form,

$\epsilon^M(\omega)=1+\frac{4 \pi}{\Omega_0}\left\langle\mu\left|\left[\frac{1}{\omega+\mathrm{i} \eta+\hat{H}^{\mathrm{exc}}}-\frac{1}{\omega+\mathrm{i} \eta-\hat{H}^{\mathrm{exc}}}\right]\right| \mu\right\rangle$

by using the spectral decomposition $\left[\hat{H}^{\mathrm{exc}}-\omega\right]^{-1}=\sum_\lambda \frac{\left|A_\lambda\right\rangle\left\langle A_\lambda\right|}{E_\lambda-\omega}$. Then, one can bring the new expression of $\epsilon(\omega)$ into a time-dependent integral, by using

$\frac{1}{\omega+\mathrm{i} \eta-\hat{H}^{\mathrm{exc}}}|\mu\rangle=-\mathrm{i} \int_0^{\infty} e^{-\mathrm{i}\left(\omega-\hat{H}^{\mathrm{exc}}+\mathrm{i} \eta\right) t}|\mu\rangle=-\mathrm{i} \int_0^{\infty} e^{-\mathrm{i}(\omega+\mathrm{i} \eta) t} e^{\mathrm{i} \hat{H}^{\mathrm{exc}}t}|\mu\rangle$,

and recognising that $e^{\mathrm{i} \hat{H}^{\mathrm{exc}}t}|\mu\rangle = |\xi(t)\rangle$ is the exponential form of a time-dependent equation. These considerations allow the expression of $\epsilon_{ij}(\omega)$ to be written as

$\epsilon_{ij}(\omega)=\delta_{ij}-\frac{4\pi e^2}{\Omega}\int_0^{\infty} \mathrm{d} t \sum_{c,v,\mathbf{k}}\left(\langle\mu^j_{cv\mathbf{k}}| \xi^i_{cv\mathbf{k}}(t)\rangle+ \mathrm{c.c.}\right) e^{-\mathrm i(\omega-\mathrm i \delta) t}$,

The fundamental aspect behind this transformation is that the new, time-dependent vector $\left.\mid \xi^j(t)\right\rangle$ follows the equation

$\mathrm i \frac{\mathrm d}{\mathrm d t}\left|\xi^j(t)\right\rangle=\hat{H}(t)\left|\xi^j(t)\right\rangle,$

with the initial vector elements given by $\left|\xi^j(0)\right\rangle=\left|\mu^j\right\rangle$. It is worthy to note that $H(t)$ is split into two parts, $H(t) = H_0 + V(t)$, where $H_0$ is the ground-state Hamiltonian, and $V(t)$ is the time-dependent perturbation.

The delta-like perturbation

explain how the system is perturbed with the delta-potential

The many-body terms in the hamiltonian

Screened two-particle interaction

Ladder diagrams from many-body perturbation theory

explain what different components are included in H (LFXC,LHARTREE,LADDER, or none)