Time-dependent density-functional theory

Time-dependent density-functional theory (TDDFT) is an extension of density-functional theory 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, for a fixed initial state, the time-dependent external potential [math]\displaystyle{ v_\mathrm{ext}(\mathbf r, t) }[/math] is uniquely determined by the time-dependent density [math]\displaystyle{ n(\mathbf r, t) }[/math]. As a consequence, all physical observables are functionals of the density and the initial state, 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).

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) + \frac{\delta(t_1-t_2)}{|\mathbf r_1 - \mathbf r_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], yields the Dyson equation

[math]\displaystyle{ \chi(1,2) = \chi_0(1,2) + \chi_0(1,3)\left[\frac{\delta(t_3-t_4)}{|\mathbf r_3 - \mathbf r_4|} + f_\mathrm{xc}(3,4)\right]\chi(4,2). }[/math]

Casida equation formalism for TDDFT

It is often useful to rewrite this equation 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) = \delta(t_3-t_4)/|\mathbf r_3 - \mathbf r_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[\frac{\delta(t_5-t_6)}{|\mathbf r_5 - \mathbf r_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 equatio for [math]\displaystyle{ \chi(1,2) }[/math] to 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'|\frac{1}{|\mathbf r_1-\mathbf r_2|}|vc'\rangle - \langle cv'|f_\mathrm{xc}(1,2)|c'v\rangle, }[/math]

[math]\displaystyle{ B_{vc}^{v'c'} = \langle vv'|\frac{1}{|\mathbf r_1-\mathbf r_2|}|cc'\rangle - \langle vv'|f_\mathrm{xc}(1,2)|c'c\rangle. }[/math]

Here, [math]\displaystyle{ v }[/math] and [math]\displaystyle{ c }[/math] denote valence and conduction states, respectively, and [math]\displaystyle{ \varepsilon_v }[/math] and [math]\displaystyle{ \varepsilon_c }[/math] are the corresponding eigenvalues. The [math]\displaystyle{ A }[/math] matrix describes the resonant (excitations) and anti-resonant (de-excitation) transitions, while the [math]\displaystyle{ B }[/math] matrix deals with the coupling between both. Due to the presence of 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 calculatios at finite momentum q the full BSE equaiton must be solved.

Connection to the dielectric function

The physical observable measured in optical experiments is the macroscopic dielectric function [math]\displaystyle{ \varepsilon_M(\omega) }[/math], which is related to the density-density response function [math]\displaystyle{ \chi }[/math] by

[math]\displaystyle{ \varepsilon^{-1}_{\mathbf G, \mathbf G'}(\mathbf q, \omega) = \delta_{\mathbf G, \mathbf G'} + \frac{4\pi}{|\mathbf q + \mathbf G|\,|\mathbf q + \mathbf G'|} \chi_{\mathbf G, \mathbf G'}(\mathbf q, \omega). }[/math]

The macroscopic dielectric function follows from [math]\displaystyle{ \varepsilon_M(\omega) = 1/\varepsilon^{-1}_{\mathbf G = \mathbf G' = \mathbf 0}(\mathbf q \to \mathbf 0, \omega) }[/math], where the inversion of the full [math]\displaystyle{ \varepsilon^{-1}_{\mathbf G, \mathbf G'} }[/math] matrix automatically includes local-field effects (the off-diagonal [math]\displaystyle{ \mathbf G \ne \mathbf G' }[/math] components). Equivalently, the eigenvalues [math]\displaystyle{ \omega_\lambda }[/math] and eigenvectors [math]\displaystyle{ X_\lambda }[/math] of the Casida equation give the excitation energies and oscillator strengths that determine the position and intensity of peaks in the optical absorption spectrum [math]\displaystyle{ \mathrm{Im}\,\varepsilon_M(\omega) }[/math].

Approximation hierarchy: IPA, RPA, and TDDFT

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): If both the Coulomb kernel [math]\displaystyle{ 1/|\mathbf r - \mathbf r'| }[/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 \approx \chi_0 }[/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): 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)\frac{\delta(t_3-t_4)}{|\mathbf r_3 - \mathbf r_4|}\chi(4,2). }[/math]

This approximation includes the long-range Coulomb interaction between electrons and holes (electron-hole attraction and repulsion), 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, local-field effects, and continuum resonances, but typically fails to describe bound excitons in semiconductors and insulators.

• Full TDDFT: Retaining both the Coulomb kernel and a non-trivial exchange-correlation kernel [math]\displaystyle{ f_\mathrm{xc} }[/math] yields the complete TDDFT response. The choice of [math]\displaystyle{ f_\mathrm{xc} }[/math] determines whether bound excitons, excitonic continuum enhancement, and other many-body effects are accurately captured (see § Approximations for the exchange-correlation kernel).

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. Instead of constructing and diagonalizing the full excitonic Hamiltonian in transition space, this 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. This narrow pulse excites all possible valence-to-conduction transitions simultaneously, while the constant displacement field replicates the long-wavelength limit [math]\displaystyle{ \mathbf q \to 0 }[/math].

The time-dependent wavefunctions are expanded as

[math]\displaystyle{ \left|\phi_i(t)\right\rangle=\left\{\left|\phi_i^0\right\rangle+\lambda\sum_{a \in \mathrm{unocc} .} c_{a i}(t)\left|\phi_a^0\right\rangle\right\} e^{-i \varepsilon_i t}, }[/math]

so that changes in an initial occupied state [math]\displaystyle{ |\phi_i^0\rangle }[/math] are captured by time-dependent contributions from unoccupied states [math]\displaystyle{ |\phi_a^0\rangle }[/math]. The time-dependent coefficients [math]\displaystyle{ c_{ai}(t) }[/math] are propagated forward in time starting from [math]\displaystyle{ c_{ai}(0)=0 }[/math] by repeatedly updating the time-dependent Hamiltonian [math]\displaystyle{ H[\phi_i(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_i(t_n)\right\rangle=\left\{\left|\phi_i^0\right\rangle+\lambda\sum_{a \in \mathrm{unocc.}} c_{ai}(t_n)\left|\phi_a^0\right\rangle\right\} e^{-i \varepsilon_i t_n}. }[/math]
2. Update the time-dependent Hamiltonian [math]\displaystyle{ H[\phi_i(t_n)] }[/math].
3. Calculate the change in the time-dependent coefficients: ::[math]\displaystyle{ \delta c_{ai}(t_n) = \left\langle\phi_a^0\right| H[\phi_i(t_n)] \left|\phi_i(t_n)\right\rangle. }[/math]
4. Compute the coefficients at the next time step: ::[math]\displaystyle{ c_{ai}(t_{n+1}) = c_{ai}(t_{n-1}) + 2\mathrm i \, \Delta t \, \delta c_{ai}(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.

Dielectric function as a time-dependent integral

The connection to the macroscopic dielectric function follows from rewriting [math]\displaystyle{ \varepsilon_M(\omega) }[/math] as a time-dependent integral. Starting from the operator form

[math]\displaystyle{ \varepsilon^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, }[/math]

where [math]\displaystyle{ \hat H^\mathrm{exc} }[/math] is the excitonic Hamiltonian and [math]\displaystyle{ \mu_{c v \mathbf k}^j = \langle c \mathbf k | v_j | v \mathbf k\rangle/(\varepsilon_c(\mathbf k)-\varepsilon_v(\mathbf k)) }[/math] is the dipole moment associated with the transition [math]\displaystyle{ v \to c }[/math] at [math]\displaystyle{ \mathbf k }[/math], the resolvent can be expressed as a time integral using the identity

[math]\displaystyle{ \frac{1}{\omega+\mathrm{i} \eta-\hat{H}^{\mathrm{exc}}}|\mu\rangle=-\mathrm{i} \int_0^{\infty} \mathrm d t \, e^{-\mathrm{i}(\omega+\mathrm{i} \eta) t} \, e^{\mathrm{i} \hat{H}^{\mathrm{exc}}t}|\mu\rangle. }[/math]

Identifying [math]\displaystyle{ |\xi(t)\rangle = e^{\mathrm i \hat H^\mathrm{exc} t}|\mu\rangle }[/math] as the time-evolved dipole vector, the dielectric function takes the form

[math]\displaystyle{ \varepsilon_{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 \eta) t}. }[/math]

The vector [math]\displaystyle{ |\xi^j(t)\rangle }[/math] satisfies a Schrödinger-like equation of motion

[math]\displaystyle{ \mathrm i \frac{\mathrm d}{\mathrm d t}\left|\xi^j(t)\right\rangle=\hat{H}^\mathrm{exc}(t)\left|\xi^j(t)\right\rangle, }[/math]

with the initial condition [math]\displaystyle{ |\xi^j(0)\rangle = |\mu^j\rangle }[/math]. In practice, [math]\displaystyle{ |\xi^j(t)\rangle }[/math] is propagated in time and its projections onto the dipole vectors are accumulated at each step. Because all operations are of matrix-vector type, this approach scales as [math]\displaystyle{ O(N^2) }[/math] compared to [math]\displaystyle{ O(N^3) }[/math] for the Casida eigendecomposition, making it significantly more efficient for systems requiring large numbers of bands or k points.

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}^\mathrm{adia}(\mathbf r, \mathbf r', t-t') = \delta(t-t') \left.\frac{\delta v_\mathrm{xc}[n](\mathbf r)}{\delta n(\mathbf r')}\right|_{n=n_0(\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 PBE gives APBE, and so on. All TDDFT kernels discussed below are adiabatic. The adiabatic approximation works well for single excitations dominated by transitions between Kohn-Sham orbitals, but fails for processes that require true memory effects, such as double excitations or strong nonadiabatic dynamics.

Local exchange-correlation kernel

VASP can include the local exchange-correlation kernel, [math]\displaystyle{ f_\mathrm{xc} }[/math], in TDDFT calculations:

[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 often lack the long-range component (which goes as [math]\displaystyle{ -1/q^2 }[/math], where [math]\displaystyle{ q }[/math] is the momentum difference between the electron and the hole). When using them in periodic or extended systems, they will likely 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 the ground-state calculation, 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}}\left(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3, \mathbf{r}_4\right) = -\frac{\delta^2 E_{\mathrm{x}}^{\text{exact}}}{\delta \rho(\mathbf{r}_1,\mathbf{r}_3) \, \delta \rho(\mathbf{r}_2,\mathbf{r}_4)}, }[/math]

where [math]\displaystyle{ E_{\mathrm{x}}^{\text{exact}} }[/math] is the exact-exchange energy. Because the exact-exchange energy is naturally a functional of the one-particle density matrix [math]\displaystyle{ \rho(\mathbf{r},\mathbf{r}') }[/math] rather than of the density alone, taking the second functional derivative yields a kernel that depends on four spatial coordinates rather than two. This four-point structure is analogous to that of the BSE kernel and is the origin of the long-range [math]\displaystyle{ -1/q^2 }[/math] behavior, allowing the kernel to properly describe bound electron-hole pairs in solids and large molecules.

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}}^{\text{NQ}} \sim \chi_0^{-1} \left[ G\, G\, W\, G\, G \right] \chi_0^{-1}, }[/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, evaluating it requires computing the screened interaction [math]\displaystyle{ W }[/math], so it is essentially as expensive as a BSE calculation and does not provide a computational advantage over BSE itself.

TDDFT versus 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. The Coulomb part is identical in both methods; 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: a two-electron integral over the Coulomb interaction weighted by the orbital products. 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 function, [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 and frequency-dependent.
• 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 but does not depend on the system's dielectric response.

In this sense, hybrid-functional TDDFT can be viewed as a BSE-like calculation in which the screening is approximated by a system-independent constant or model function. 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, material-independent screening.

Related tags and articles

How-to

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

Theory

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