Time-dependent Hartree-Fock method (response of density matrices)
There are many different ways to derive the linear response time-dependent HF (LR-TDHF) or LR-TDDFT equation. This derivation is based on the linear response of HF density matrices to an external perturbation following the derivations in the famous paper. The perturbation and the response of density matrices have the same form (same frequency) based on the linear response theory.
For any single-electron wave function, the time-dependent Schrödinger equation can be written as
\[i\hbar\frac{\partial}{\partial t}|\phi_{p}\rangle=\hat{f}(t)|\phi_{p}\rangle\]If we expand \(\phi_{p}\) with a set of time-independent basis functions \(\{|\mu\rangle\}\), we have
\[\sum_{\mu}i\hbar|\mu\rangle\frac{\partial}{\partial t}c_{\mu p}(t)=\sum_{\mu}\hat{f}(t)c_{\mu p}(t)|\mu\rangle\]Left project basis function and we have
\[i\hbar SC^{\prime}(t)=F(t)C(t)\]and its complex conjugate
\[-i\hbar C^{^{\dagger}\prime}(t)S=C^{^{\dagger}}(t)F(t)\]We have
\[\begin{aligned}i\hbar\frac{\partial}{\partial t}D(t)=i\hbar\left(\frac{\partial}{\partial t}C(t)\tilde{I}C^{\dagger}(t)\right)= & i\hbar\left(\frac{\partial}{\partial t}C(t)\right)\tilde{I}C^{\dagger}(t)+i\hbar C(t)\tilde{I}\left(\frac{\partial}{\partial t}C^{\dagger}(t)\right)\\ = & S^{-1}F(t)D(t)-D(t)F(t)S^{-1} \end{aligned}\]For a given orthonormal basis set, e.g., the unperturbed single-electron wave function, we have (in a.u.)
\[i\frac{\partial}{\partial t}D(t)=\left[F(t),D(t)\right]\]Expand the above equation with perturbation series.
\[D(t) =D^{(0)}+D^{(1)}(t)+\cdots\] \[F(t) =F^{(0)}+F^{(1)}(t)+\cdots\]and collect the first-order terms
\[i\frac{\partial}{\partial t}D^{(1)}(t)=F^{(0)}D^{(1)}(t)+F^{(1)}(t)D^{(0)}-D^{(0)}F^{(1)}(t)-D^{(1)}(t)F^{(0)}\]With the perturbation Hamiltonian and the response of density matrices written in frequency-domain
\[g_{pq}=\frac{1}{2}\left[f_{pq}e^{-i\omega t}+f_{qp}^{*}e^{i\omega t}\right]\] \[D_{pq}^{(1)}=\frac{1}{2}\left[d_{pq}e^{-i\omega t}+d_{qp}^{*}e^{i\omega t}\right]\]We note that \(F^{(1)}\neq g\) since it has another first-order contribution from the two-electron interaction.
\[F_{pq}^{(1)}=g_{pq}+\sum_{rs}\langle pr||qs\rangle D_{rs}^{(1)}\]Collect the terms with \(e^{-i\omega t}\)
\[\begin{aligned}\omega d_{pq}= & \sum_{r}F_{pr}^{(0)}d_{rq}+\sum_{r}\left[f_{pr}+\sum_{mn}\langle pm||rn\rangle d_{mn}\right]D_{rq}^{(0)}\\ - & \sum_{r}d_{pr}F_{rq}^{(0)}-\sum_{r}D_{pr}^{(0)}\left[f_{rq}+\sum_{mn}\langle rm||qn\rangle d_{mn}\right] \end{aligned}\]Noticing that \(DD=C\tilde{I}C^{\dagger}C\tilde{I}C^{\dagger}=D\), we have
\[D^{(0)}D^{(1)}+D^{(1)}D^{(0)}=D^{(1)}\]with canonical orbitals, \(F_{pq}^{(0)}=\epsilon_{p}\delta_{pq},D_{pq}^{(0)}=\delta_{pq}\delta_{pi}\). We thus have the following restrictions on \(D^{(1)}\)
\[\left[\begin{array}{cc} D_{ii}^{(1)} & D_{ia}^{(1)}\\ 0 & 0 \end{array}\right]+\left[\begin{array}{cc} D_{ii}^{(1)} & 0\\ D_{ai}^{(1)} & 0 \end{array}\right]=\left[\begin{array}{cc} D_{ii}^{(1)} & D_{ia}^{(1)}\\ D_{ai}^{(1)} & D_{aa}^{(1)} \end{array}\right]\]As a result, we can only have non-zero \(D_{ai}^{(1)}\) and \(D_{ia}^{(1)}\).
\[\omega d_{ai}=\left(\epsilon_{a}-\epsilon_{i}\right)d_{ai}+f_{ai}+\sum_{bj}\langle ab||ij\rangle d_{jb}+\sum_{bj}\langle aj||ib\rangle d_{bj}\] \[\omega d_{ia}=\left(\epsilon_{i}-\epsilon_{a}\right)d_{ia}-f_{ia}-\sum_{bj}\langle ib||aj\rangle d_{jb}-\sum_{bj}\langle aj||ib\rangle d_{bj}\]Take \(f\rightarrow0\) (zero-frequency limit) and rename \(x_{ai}=d_{ai},y_{ai}=d_{ia}\)
\[\left[\begin{array}{cc} A & B\\ B^{*} & A^{*} \end{array}\right]\left[\begin{array}{c} X\\ Y \end{array}\right]=\omega\left[\begin{array}{cc} I & 0\\ 0 & -I \end{array}\right]\left[\begin{array}{c} X\\ Y \end{array}\right]\]where
\[A_{ai,bj}=\left(\epsilon_{a}-\epsilon_{i}\right)\delta_{ab}\delta_{ij}+\langle aj||ib\rangle\] \[B_{ai,bj}=\langle ab||ij\rangle\]