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 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 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) + \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]

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): 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(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): 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, 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: Retaining both the Coulomb kernel 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[\frac{\delta(t_3-t_4)}{|\mathbf r_3 - \mathbf r_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

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) = \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 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'|\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 (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 = E_\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 + E_\lambda + \mathrm i\eta} - \frac{1}{\omega - E_\lambda + \mathrm i\eta}\right], }[/math]

where [math]\displaystyle{ \mu_{cv\mathbf k} }[/math] are the dipole matrix elements. The eigenvalues give the excitation energies (peak positions) and the eigenvectors determine the oscillator strengths (peak intensities) in the optical absorption spectrum [math]\displaystyle{ \mathrm{Im}\,\varepsilon_M(\omega) }[/math].

Alternatively, the macroscopic dielectric function [math]\displaystyle{ \varepsilon_M(\omega) }[/math] 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 limit 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). Using the Casida eigenvectors, this can be written as

[math]\displaystyle{ \varepsilon_M(\mathbf{q},\omega)= 1+2\lim_{\mathbf{q}\rightarrow 0}v(q)\sum_{\lambda} \left|\sum_{c,v,\mathbf k}\langle c\mathbf{k}|e^{\mathrm i\mathbf{qr}}|v\mathbf{k}\rangle X_\lambda^{cv\mathbf{k}}\right|^2 \left(\frac{1}{\omega_\lambda - \omega - \mathrm i\delta}\right). }[/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. The starting point is the time-dependent Kohn-Sham equation,

[math]\displaystyle{ \mathrm i \frac{\partial}{\partial t}\left|\phi_i[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_i[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. 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].

Propagation of the time-dependent coefficients

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.

The dielectric function as a time-dependent integral

The key theoretical insight behind the time-evolution approach is that the macroscopic dielectric function can be recast as a time-dependent integral. Starting from 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 + E_\lambda + \mathrm i\eta} - \frac{1}{\omega - E_\lambda + \mathrm i\eta}\right], }[/math]

where [math]\displaystyle{ \mu_{cv\mathbf k} }[/math] are the dipole matrix elements, [math]\displaystyle{ A^\lambda }[/math] are the eigenvectors and [math]\displaystyle{ E_\lambda }[/math] are the eigenvalues of the excitonic Hamiltonian [math]\displaystyle{ H^\mathrm{exc} }[/math], one can write this in operator form as

[math]\displaystyle{ \varepsilon_M(\omega) = 1 + \frac{4\pi}{\Omega} \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]

Using the identity

[math]\displaystyle{ \frac{1}{\omega + \mathrm i\eta - \hat{H}^\mathrm{exc}}|\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 \, \mathrm d t, }[/math]

and recognizing that [math]\displaystyle{ |\xi(t)\rangle = e^{\mathrm i \hat{H}^\mathrm{exc} t}|\mu\rangle }[/math] satisfies the time-dependent equation

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

with the initial condition [math]\displaystyle{ |\xi(0)\rangle = |\mu\rangle }[/math], the dielectric function becomes

[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}|\xi^i_{cv\mathbf k}(t)\rangle + \mathrm{c.c.}\right) e^{-\mathrm i(\omega - \mathrm i\eta)t}. }[/math]

The dielectric function is therefore obtained by propagating the time-dependent vector [math]\displaystyle{ |\xi(t)\rangle }[/math] and accumulating its overlap with the dipole vector [math]\displaystyle{ |\mu\rangle }[/math] at each time step, followed by a Fourier transform. Since all operations in this scheme are of the matrix-vector type, the computational cost scales as [math]\displaystyle{ O(N_\mathrm{rank}^2) }[/math] rather than [math]\displaystyle{ O(N_\mathrm{rank}^3) }[/math] for the eigendecomposition approach.

Two-point Dyson equation for TDDFT

Another approach to TDDFT 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

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 lack the long-range component (which goes as [math]\displaystyle{ -1/q^2 }[/math], where [math]\displaystyle{ q }[/math] is the momentum transfer between the electron and the hole). 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 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 n(\mathbf{r}_1,\mathbf{r}_3) \, \delta n(\mathbf{r}_2,\mathbf{r}_4)}, }[/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}}^{(2)}(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.

Time-evolution TDDFT and BSE

The time-evolution approach can also be combined with BSE-level interactions. In this case, the Hamiltonian used during propagation includes the screened Coulomb interaction [math]\displaystyle{ W }[/math] instead of (or in addition to) an exchange-correlation kernel. This yields optical spectra equivalent to solving the full BSE eigenvalue problem, but avoids explicit construction and diagonalization of the excitonic Hamiltonian. The computational advantage is that the time-evolution approach scales more favorably with system size than direct diagonalization of the BSE matrix.

Related tags and articles

How-to

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

Theory

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