Time-dependent density-functional theory

Time-dependent density-functional theory (TDDFT) is an extension of density-functional theory (DFT) to systems with time-varying external potentials, enabling the computation of excited-state properties and response functions. TDDFT calculations can be based on ground-state electronic structures obtained from DFT, hybrid functionals, or even GW approximations.

The theoretical foundation of TDDFT is the Runge-Gross theorem, which is the time-dependent analog of the Hohenberg-Kohn theorem of density-functional theory. It states that the time-dependent density [math]\displaystyle{ n(\mathbf r, t) }[/math] is uniquely determined by the time-dependent external potential [math]\displaystyle{ v_\mathrm{ext}(\mathbf r, t) }[/math]. As a consequence, all physical observables are functionals of the density, and the interacting many-electron problem can be mapped onto an auxiliary system of non-interacting electrons reproducing the same density (the time-dependent Kohn-Sham system).

Linear-response formalism

In the linear-response regime, the external potential is split into a static part and a small time-dependent perturbation, [math]\displaystyle{ v_\mathrm{ext}(\mathbf r, t) = v(\mathbf r) + \delta v(\mathbf r, t) }[/math]. The induced density variation [math]\displaystyle{ \delta n }[/math] is then related to the perturbation through the density-density response function [math]\displaystyle{ \chi }[/math],

[math]\displaystyle{ \delta n(1) = \int \mathrm d 2 \, \chi(1,2) \, \delta v(2), }[/math]

where [math]\displaystyle{ 1 \equiv (\mathbf{r}_1, t_1) }[/math] and [math]\displaystyle{ 2 \equiv (\mathbf{r}_2, t_2) }[/math] denote space-time coordinates. Equivalently, [math]\displaystyle{ \chi }[/math] is the functional derivative [math]\displaystyle{ \chi(1,2) = \delta n(1)/\delta v_\mathrm{ext}(2) }[/math]. Within the Kohn-Sham scheme, the density response of the interacting system equals that of an auxiliary non-interacting system responding to an effective perturbation [math]\displaystyle{ \delta v_\mathrm{KS} = \delta v + \delta v_\mathrm{H} + \delta v_\mathrm{xc} }[/math], which defines the independent-particle response

[math]\displaystyle{ \chi_0(1,2) = \frac{\delta n(1)}{\delta v_\mathrm{KS}(2)}. }[/math]

Applying the chain rule to [math]\displaystyle{ \chi(1,2) }[/math] and using

[math]\displaystyle{ \frac{\delta v_\mathrm{KS}(1)}{\delta v_\mathrm{ext}(2)} = \delta(1,2) + v(1,2) + f_\mathrm{xc}(1,2), }[/math]

with the exchange-correlation kernel [math]\displaystyle{ f_\mathrm{xc}(1,2) = \delta v_\mathrm{xc}(1)/\delta n(2) }[/math] and the bare Coulomb interaction [math]\displaystyle{ v(1,2) = \delta(t_1-t_2)/|\mathbf r_1 - \mathbf r_2| }[/math], yields the Dyson equation

[math]\displaystyle{ \chi(1,2) = \chi_0(1,2) + \chi_0(1,3)\left[v(3,4) + f_\mathrm{xc}(3,4)\right]\chi(4,2). }[/math]

Approximation hierarchy

The response function formalism above can be systematically simplified by neglecting different interaction terms in the Dyson equation for [math]\displaystyle{ \chi }[/math]:

• Independent-particle approximation (IPA) (LHARTREE=.FALSE. and LFXC=.FALSE.): If both the Coulomb interaction [math]\displaystyle{ v }[/math] and the exchange-correlation kernel [math]\displaystyle{ f_\mathrm{xc} }[/math] are neglected, the full response function reduces to the independent-particle response,

[math]\displaystyle{ \chi(1,2) = \chi_0(1,2). }[/math]

In this limit, optical spectra are computed from non-interacting Kohn-Sham transitions without any electron-hole interaction or local-field effects.

• Random-phase approximation (RPA) (LHARTREE=.TRUE. and LFXC=.FALSE.): If [math]\displaystyle{ f_\mathrm{xc} }[/math] is neglected but the Coulomb kernel is retained, the Dyson equation becomes

[math]\displaystyle{ \chi(1,2) = \chi_0(1,2) + \chi_0(1,3)\,v(3,4)\,\chi(4,2). }[/math]

This approximation includes the long-range Coulomb interaction, but omits exchange-correlation contributions beyond those already present in the Kohn-Sham eigenvalues used to construct [math]\displaystyle{ \chi_0 }[/math]. RPA captures plasmons and local-field effects, but typically fails to describe bound excitons in semiconductors and insulators.

• Full TDDFT (LHARTREE=.TRUE. and LFXC=.TRUE.; additionally LADDER=.TRUE. if a hybrid functional is used): Retaining both the Coulomb interaction and an exchange-correlation kernel [math]\displaystyle{ f_\mathrm{xc} }[/math] yields the complete TDDFT response,

[math]\displaystyle{ \chi(1,2) = \chi_0(1,2) + \chi_0(1,3)\left[v(3,4) + f_\mathrm{xc}(3,4)\right]\chi(4,2). }[/math]

The choice of [math]\displaystyle{ f_\mathrm{xc} }[/math] determines whether bound excitons and other many-body effects are accurately captured.

Casida equation formalism for TDDFT

The Casida approach (ALGO=TDHF) solves for the excitation spectrum by diagonalizing the excitonic Hamiltonian in the transition space. It is often useful to rewrite the Dyson equation for [math]\displaystyle{ \chi }[/math] in terms of a four-point response function [math]\displaystyle{ L(1,2,3,4) }[/math], which describes the response of the system to a two-particle perturbation. The four-point function is related to the two-point density response [math]\displaystyle{ \chi }[/math] by

[math]\displaystyle{ \chi(1,2) = \int \mathrm d 3 \, \mathrm d 4 \, L(1,3,2,4) \, v(3,4), }[/math]

where [math]\displaystyle{ v(3,4) }[/math] is the bare Coulomb interaction. The four-point response function satisfies a Dyson-like equation

[math]\displaystyle{ L(1,2,3,4) = L_0(1,2,3,4) + L_0(1,2,5,6) \left[v(5,6) + f_\mathrm{xc}(5,6)\right] L(5,6,3,4), }[/math]

where [math]\displaystyle{ L_0 }[/math] is the independent-particle four-point response. This four-point formulation is particularly useful when comparing Casida TDDFT with the Bethe-Salpeter equation, as both can be expressed in terms of analogous four-point functions.

Working within the frequency domain and using a basis set that considers transitions from valence to conduction states at the same k point (transition basis) makes it possible to recast the equation for [math]\displaystyle{ \chi(1,2) }[/math] into an eigenvalue problem

[math]\displaystyle{ \left( \begin{matrix} A & B \\ -B^* & -A^* \end{matrix} \right) \left( \begin{matrix} X \\ Y \end{matrix} \right) = \omega \left( \begin{matrix} X \\ Y \end{matrix} \right), }[/math]

which is also known as the Casida equation. The [math]\displaystyle{ A }[/math] and [math]\displaystyle{ B }[/math] matrices are given by

[math]\displaystyle{ A_{vc}^{v'c'} = (\varepsilon_v-\varepsilon_c)\delta_{vv'}\delta_{cc'} + \langle cv'|v|vc'\rangle - \langle cv'|f_\mathrm{xc}|c'v\rangle, }[/math]

[math]\displaystyle{ B_{vc}^{v'c'} = \langle vv'|v|cc'\rangle - \langle vv'|f_\mathrm{xc}|c'c\rangle. }[/math]

Here, [math]\displaystyle{ v }[/math] and [math]\displaystyle{ c }[/math] denote valence and conduction states, respectively, [math]\displaystyle{ \varepsilon_v }[/math] and [math]\displaystyle{ \varepsilon_c }[/math] are the corresponding eigenvalues, and the two-electron integrals [math]\displaystyle{ \langle \cdot|\cdot|\cdot\rangle }[/math] involve the bare Coulomb interaction or the exchange-correlation kernel as indicated. The [math]\displaystyle{ A }[/math] matrix describes the resonant (excitation) and anti-resonant (de-excitation) transitions, while the [math]\displaystyle{ B }[/math] matrix describes the coupling between them. Due to this coupling, the Casida matrix is non-Hermitian.

Tamm-Dancoff approximation

The non-Hermitian structure of the Casida equation arises from the coupling between resonant and anti-resonant transitions through the off-diagonal block [math]\displaystyle{ B }[/math]. A common simplification is the Tamm-Dancoff approximation (TDA), in which the coupling block is set to zero, [math]\displaystyle{ B = 0 }[/math]. The eigenvalue problem then reduces to the Hermitian form

[math]\displaystyle{ A \, X_\lambda = \omega_\lambda \, X_\lambda. }[/math]

The TDA significantly reduces the computational cost and guarantees real, positive excitation energies, but it neglects the mixing between excitations and de-excitations. For optical absorption spectra of solids it is usually a good approximation, but for calculations at finite momentum q the full equation must be solved.

Connection to the dielectric function

The macroscopic dielectric function is obtained by using the eigenvalues [math]\displaystyle{ \omega_\lambda }[/math] and eigenvectors [math]\displaystyle{ A^\lambda }[/math] of the Casida equation in the spectral representation

[math]\displaystyle{ \varepsilon_M(\omega) = 1 + \frac{4\pi}{\Omega} \sum_\lambda \left|\sum_{cv\mathbf k} \mu_{cv\mathbf k} A_{cv\mathbf k}^\lambda\right|^2 \left[\frac{1}{\omega + \omega_\lambda + \mathrm i\eta} - \frac{1}{\omega - \omega_\lambda + \mathrm i\eta}\right], }[/math]

where [math]\displaystyle{ \mu_{cv\mathbf{k}}^j=\frac{\langle c\mathbf{k}|v_j|v\mathbf{k}\rangle}{\varepsilon_c(\mathbf{k})-\varepsilon_v(\mathbf{k})} }[/math] is the dipole matrix element associated with the transition from valence state [math]\displaystyle{ v }[/math] to conduction state [math]\displaystyle{ c }[/math] at k-point [math]\displaystyle{ \mathbf k }[/math], [math]\displaystyle{ v_j }[/math] is the velocity operator along direction [math]\displaystyle{ j }[/math], and [math]\displaystyle{ \omega_\lambda }[/math] are the excitation energies. The eigenvalues give the peak positions and the eigenvectors determine the oscillator strengths (peak intensities) in the optical absorption spectrum [math]\displaystyle{ \mathrm{Im}\,\varepsilon_M(\omega) }[/math].

Time-evolution or real-time TDDFT

An alternative to solving the Casida equation is to compute the frequency-dependent response via real-time propagation of the Kohn-Sham orbitals (ALGO=TIMEEV). The starting point is the time-dependent Kohn-Sham equation,

[math]\displaystyle{ \mathrm i \frac{\partial}{\partial t}\left|\phi_{v\mathbf k}[n(t)]\right\rangle = \left[-\frac{\nabla^2}{2} + V_{\mathrm H}[n(t)] + V_{\mathrm{xc}}[n(t)] + V_\mathrm{ext}(t)\right]\left|\phi_{v\mathbf k}[n(t)]\right\rangle, }[/math]

where all potentials are functionals of the time-evolving density [math]\displaystyle{ n(\mathbf r, t) }[/math].

Instead of constructing and diagonalizing the full Hamiltonian in transition space as done in Casida formalism, the time-evolution method applies a delta-like perturbation

[math]\displaystyle{ V_\mathrm{ext}(\mathbf r, t) = \lambda \, \mathbf r\cdot \mathbf D \, \delta(t) }[/math]

to the ground-state system, where [math]\displaystyle{ \lambda }[/math] is a small perturbation parameter and [math]\displaystyle{ \mathbf D }[/math] is the electric displacement field.

Propagation of the time-dependent coefficients

The time-dependent wavefunctions are expanded as

[math]\displaystyle{ \left|\phi_{v\mathbf k}(t)\right\rangle=\left\{\left|\phi_{v\mathbf k}^0\right\rangle+\lambda\sum_{c \in \mathrm{unocc} .} c_{cv\mathbf k}(t)\left|\phi_{c\mathbf k}^0\right\rangle\right\} e^{-\mathrm i \varepsilon_{v\mathbf{k}} t}, }[/math]

so that changes in an initial occupied state [math]\displaystyle{ |\phi_{v\mathbf k}^0\rangle }[/math] are captured by time-dependent contributions from unoccupied states [math]\displaystyle{ |\phi_{c\mathbf k}^0\rangle }[/math]. The time-dependent coefficients [math]\displaystyle{ c_{cv\mathbf k}(t) }[/math] are propagated forward in time starting from [math]\displaystyle{ c_{cv\mathbf k}(0)=0 }[/math] by repeatedly updating the time-dependent Hamiltonian [math]\displaystyle{ H[\phi_{v\mathbf k}(t)] }[/math].

In essence, the time-evolution algorithm works as follows:

1. Set up states at time step [math]\displaystyle{ t_n }[/math]: [math]\displaystyle{ \left|\phi_{v\mathbf k}(t_n)\right\rangle=\left\{\left|\phi_{v\mathbf k}^0\right\rangle+\lambda\sum_{c \in \mathrm{unocc.}} c_{cv\mathbf k}(t_n)\left|\phi_{c\mathbf k}^0\right\rangle\right\} e^{-\mathrm i \varepsilon_{v\mathbf{k}} t_n}. }[/math]
2. Update the time-dependent Hamiltonian [math]\displaystyle{ H[\phi_{v\mathbf k}(t_n)] }[/math].
3. Calculate the change in the time-dependent coefficients: [math]\displaystyle{ \delta c_{cv\mathbf k}(t_n) = \left\langle\phi_{c\mathbf k}^0\right| H[\phi_{v\mathbf k}(t_n)] \left|\phi_{v\mathbf k}(t_n)\right\rangle. }[/math]
4. Compute the coefficients at the next time step: [math]\displaystyle{ c_{cv\mathbf k}(t_{n+1}) = c_{cv\mathbf k}(t_{n-1}) + 2\mathrm i \, \Delta t \, \delta c_{cv\mathbf k}(t_n). }[/math]

Here, [math]\displaystyle{ \Delta t }[/math] is the time step chosen for the propagation. Updating the Hamiltonian with the new states is essential, as it is a functional of the time-dependent density.

The same approximation hierarchy described in the Approximation hierarchy section applies to the time-evolution approach. The tags LHARTREE, LFXC, and LADDER control which interaction terms are included in the time-dependent Hamiltonian: LHARTREE controls the Coulomb (Hartree) energy, LFXC controls the local exchange-correlation kernel, and LADDER controls the nonlocal exchange contribution from hybrid functionals. At each time step, the corresponding energy terms in the Hamiltonian are updated based on the time-evolving density.

Connection to the dielectric function

The dielectric function is obtained by propagating the time-dependent coefficients [math]\displaystyle{ c_{cv\mathbf k}(t) }[/math] and accumulating their overlap with the dipole vector [math]\displaystyle{ |\mu\rangle }[/math] at each time step, followed by a Fourier transform:

[math]\displaystyle{ \varepsilon_{ij}(\omega) = \delta_{ij} - \frac{4\pi e^2}{\Omega} \int_0^\infty \mathrm d t \sum_{cv\mathbf k} \left(\langle\mu^j_{cv\mathbf k}|c^i_{cv\mathbf k}(t)\rangle + \mathrm{c.c.}\right) e^{-\mathrm i(\omega - \mathrm i\eta)t}. }[/math]

Here, [math]\displaystyle{ \mu_{cv\mathbf{k}}^j }[/math] is the dipole matrix element defined in the connection to the dielectric function above.

Two-point Dyson equation for TDDFT

Another approach to TDDFT (ALGO=CHI) is to solve the two-point Dyson equation directly and find the density-density response function [math]\displaystyle{ \chi }[/math]. This approach allows one to calculate the full [math]\displaystyle{ \chi_{\mathbf q}(\mathbf G,\mathbf G',\omega) }[/math] and is used for finding the Coulomb interaction screened via the RPA dielectric function, W. However, it requires inverting the Dyson equation at every frequency, which becomes too costly for calculating spectra with good resolution. Since this approach does not invoke the Tamm-Dancoff approximation, the RPA dielectric function obtained this way is equivalent to the Casida RPA dielectric function.

Approximations for the exchange-correlation kernel

The exchange-correlation kernel [math]\displaystyle{ f_\mathrm{xc} }[/math] is a key quantity in TDDFT that approximates the interaction between electrons and holes. The choice of kernel depends on the exchange-correlation functional used in the ground-state calculation.

Adiabatic approximation

In its exact form, the exchange-correlation kernel [math]\displaystyle{ f_\mathrm{xc}(\mathbf r, \mathbf r', t-t') }[/math] is nonlocal in time, meaning that [math]\displaystyle{ v_\mathrm{xc} }[/math] at time [math]\displaystyle{ t }[/math] depends on the density at all earlier times [math]\displaystyle{ t' < t }[/math]. This memory dependence is intractable in practice. The adiabatic approximation replaces the time-nonlocal kernel by the instantaneous functional derivative of a ground-state exchange-correlation functional,

[math]\displaystyle{ f_\mathrm{xc}(\mathbf r, \mathbf r', t-t') = \delta(t-t') \frac{\delta v_\mathrm{xc}[n](\mathbf r)}{\delta n(\mathbf r')}, }[/math]

evaluated at the ground-state density [math]\displaystyle{ n_0 }[/math]. The kernel is therefore frequency-independent in this approximation. Combining the adiabatic approximation with the local-density approximation gives the adiabatic LDA (ALDA), with GGA gives AGGA, and so on. All TDDFT kernels discussed below are adiabatic.

Local exchange-correlation kernel

The local exchange-correlation kernel, [math]\displaystyle{ f_\mathrm{xc} }[/math], in TDDFT calculations is given by

[math]\displaystyle{ f_{\mathrm{xc}}^{\text {loc }}\left(\mathbf{r}, \mathbf{r}^{\prime}\right)=\frac{\delta^2 E_{\mathrm{xc}}^{\mathrm{DFT}}}{\delta n(\mathbf{r}) \delta n\left(\mathbf{r}^{\prime}\right)}, }[/math]

where [math]\displaystyle{ E_{\mathrm{xc}}^{\mathrm{DFT}} }[/math] is the local or semilocal exchange-correlation functional (e.g., LDA or PBE).

These local kernels lack the long-range component (which goes as [math]\displaystyle{ -1/q^2 }[/math]). In periodic or extended systems, they fail to properly reproduce the binding energies of electron-hole pairs.

The ALDA and APBE kernels work well for metallic systems and for optical properties where excitonic effects are not important (such as plasmon frequencies). However, they fail to describe bound excitons in semiconductors and insulators, where the long-range electron-hole interaction is crucial for determining excitation energies and binding energies.

Exchange-correlation kernel from exact exchange

When hybrid functionals or Hartree-Fock exchange are used in TDDFT, the exchange-correlation kernel includes a contribution from exact exchange. This nonlocal exchange contribution naturally provides the [math]\displaystyle{ -1/q^2 }[/math] long-range behavior in the kernel, which is essential for capturing excitonic effects. The exact-exchange contribution to [math]\displaystyle{ f_\mathrm{xc} }[/math] takes the form

[math]\displaystyle{ f_{\mathrm{x}}^{\text{exact}}(\mathbf{r}, \mathbf{r'}) = \frac{\delta^2 E_{\mathrm{x}}^{\text{exact}}}{\delta^2 n(\mathbf{r},\mathbf{r'})}, }[/math]

where [math]\displaystyle{ E_{\mathrm{x}}^{\text{exact}} }[/math] is the exact-exchange energy.

Alternatively, the screened exchange interaction potential [math]\displaystyle{ W(\mathbf r,\mathbf r';\omega) }[/math] from many-body perturbation theory can be used. This treats the electron-hole interaction by including the ladder diagrams from many-body perturbation theory, providing an alternative route to the correct long-range physics.

Nanoquanta kernel

The nanoquanta kernel is an exchange-correlation kernel derived from many-body perturbation theory that maps the Bethe-Salpeter equation (BSE) onto a TDDFT-like equation. Rather than being a simple analytical form, it is constructed by requiring that the TDDFT response function reproduces the BSE response function. This leads to a kernel of the form

[math]\displaystyle{ f_{\mathrm{xc}}^{NQ}(1,2)=\int \mathrm d 3456 \, \chi_0^{-1}(1,3) \, G(3,4) \, G(5,3) \, W(4,5) \, G(4,6) \, G(6,5) \, \chi_0^{-1}(6,2), }[/math]

where [math]\displaystyle{ \chi_0 }[/math] is the independent-particle response function, [math]\displaystyle{ G }[/math] is the single-particle Green's function, and [math]\displaystyle{ W }[/math] is the screened Coulomb interaction from GW theory. The four Green's functions [math]\displaystyle{ G }[/math] contract the two-point density indices on the outside with the four-point structure of [math]\displaystyle{ W }[/math] inside, effectively mapping the four-point BSE kernel onto a two-point TDDFT kernel. Because [math]\displaystyle{ W }[/math] contains the proper [math]\displaystyle{ -1/q^2 }[/math] long-range behavior, the resulting nanoquanta kernel inherits this feature and is therefore capable of describing bound excitons in semiconductors and insulators.

By construction, the nanoquanta kernel yields optical spectra of comparable accuracy to BSE calculations, including bound excitons and continuum excitonic enhancement. However, computing it is as expensive as a BSE calculation and does not provide a computational advantage over BSE itself.

See the LFXC tag page for details on how the [math]\displaystyle{ f_\mathrm{xc} }[/math] kernels are implemented in VASP and what approximations are made.

TDDFT compared to BSE

Casida TDDFT and BSE

The Casida formulation of TDDFT and the Bethe-Salpeter equation (BSE) share essentially the same mathematical structure: both are non-Hermitian eigenvalue problems of the form

[math]\displaystyle{ \left(\begin{matrix} A & B \\ -B^* & -A^* \end{matrix}\right) \left(\begin{matrix} X \\ Y \end{matrix}\right) = \omega \left(\begin{matrix} X \\ Y \end{matrix}\right), }[/math]

where the matrix elements of [math]\displaystyle{ A }[/math] and [math]\displaystyle{ B }[/math] are built from valence-to-conduction transitions and contain a Hartree (Coulomb) contribution together with a term that goes beyond the random-phase approximation (RPA). The Coulomb part is identical in both methods and the difference lies entirely in how the beyond-RPA contribution is described:

• In TDDFT, the beyond-RPA term is the exchange-correlation kernel [math]\displaystyle{ f_\mathrm{xc} }[/math].
• In BSE, the beyond-RPA term is the screened Coulomb interaction [math]\displaystyle{ W }[/math].

When [math]\displaystyle{ f_\mathrm{xc} }[/math] includes a fraction of exact exchange (as in hybrid or range-separated hybrid functionals), both formalisms involve integrals of the same form and the remaining difference is how that interaction is screened:

• In BSE, the bare Coulomb interaction [math]\displaystyle{ v }[/math] is screened by the inverse dielectric tensor, [math]\displaystyle{ W = \varepsilon^{-1} v }[/math], so that the screening is determined by the actual electronic response of the system and is in general nonlocal.
• In Casida TDDFT with hybrid functionals, the screening is replaced by a constant prefactor [math]\displaystyle{ c_\mathrm{x} }[/math], the fraction of exact exchange (e.g., 0.25 for PBE0). For range-separated hybrids, this prefactor becomes a simple function of [math]\displaystyle{ |\mathbf q + \mathbf G| }[/math] that mimics a screened interaction.

In this sense, hybrid-functional TDDFT can be viewed as a BSE-like calculation in which the screening is approximated by a model function or a constant. This makes TDDFT considerably cheaper than BSE, because it avoids the need for a preceding GW or screening calculation, but at the cost of using an approximate, diagonal screening model.

Related tags and articles

Tags

ALGO, LHARTREE, LFXC, LADDER

How-to

Time-dependent density-functional theory calculations, Plotting exciton wavefunction

Theory

Time-evolution algorithm, Bethe-Salpeter equation, Dielectric properties