Linear Response

Static Response

We will consider a Kohn-Sham equation

$$ (1)\ H^{(0)} |\psi^{(0)}_k> = \epsilon^{(0)}_k | \psi^{(0)}_k> $$

and a perturbing potential of the form

$$ V_{ext}=V_{ext}(\lambda_i)\ . $$

$$ (2)\ H|\psi_k> = \epsilon_k |\psi_k> $$

Then we can expand the Hamiltonian, the wavefunctions and the eigenvalues in terms of $\lambda_i,$

$$ H = H+\lambda_i H^{(1),i}+\lambda_i\lambda_j H^{(2),ij}+\cdots $$

$$ |\psi_k> = |\psi^{(0)}_k>+\lambda_i|\psi^{(1),i}_k>+\lambda_i\lambda_j|\psi^{(2),ij}>+\cdots $$

$$ \epsilon_k = \epsilon^{(0)}_k+\lambda_i\epsilon^{(1),i}_k+\lambda_i\lambda_j\epsilon^{(2),ij}+\cdots $$

Putting all the expansions until second order in equation (2) we get

$$ \left(H^{(0)}-e^{(0)}\right)|\psi_k^{(0)}>+ \lambda_i\left(H^{(1),i}-\epsilon^{(1),i}\right)|\psi_k^{(0)}>+ $$

$$ \lambda_i\lambda_j\left(H^{(0)}-\epsilon^{(0)}\right)|\psi_k^{(2),ij}>+ \lambda_i\left(H^{(1),i}-\epsilon^{(1),i}\right)|\psi_k^{(0)}>+ $$

$$ \lambda_i\lambda_j\left(H^{(1),i}-\epsilon^{(1),i}\right)|\psi_k^{(1),i}>+ \lambda_i\lambda_j\left(H^{(2),ij}-\epsilon^{(2),ij}\right)|\psi_k{(0)}> + O(\lambda^3)=0 \ , $$

we can separate terms according to the dependence on lambda; for the zeroth order we get the ground-state equation.

First Order

For the first order we get the Sternheimer equation

$$ \left(H^{(0)}-\epsilon^{(0)}\right)|\psi^{(1),i}_k> = -\left(H^{(1),i}-\epsilon^{(1),i}\right)|\psi^{(0)}_k> $$

The RHS can be written as

$$ \left(H^{(0)}-\epsilon^{(0)}\right)P_c|\psi^{(1),i}_k> = -P_cH^{(1),i}|\psi^{(0)}_k> $$

where $P_c$ is a projector onto the unoccupied space. We need only the components in the unoccupied space for many applications.

Second Order

Doing the same for the second-order perturbation we get

$$ \left(H^{(0)}-\epsilon^{(0)}\right)|\psi^{(2),ij}>+ \left(H^{(1),i}-\epsilon^{(1),i}\right)|\psi^{(1),j}>+ \left(H^{(2),ij}-\epsilon^{(2),ij}\right)|\psi^{(0)}>=0\ . $$

Relation to the derivatives

Now we will find the relations between perturbations as defined here and the derivatives of the magnitudes with respect to $\lambda$.

We can see from the second-order equation it’s not symmetrical with respect to the index exchange, so $|\psi^{(2),ij}>\neq|\psi^{(2),ji}>$. To relate the perturbation terms to the derivatives (which are symmetrical) we must take this in account, taking the average over all index permutations, as well as the $\frac1{n!}$ factor (which cancels the average factor), so:

$$ \left.\frac{\partial A}{\partial\lambda_i}\right|_{\lambda_i=0} = A^{(1),i} $$

$$ \left.\frac{\partial^2 A}{\partial\lambda_i\lambda_j}\right|_{\lambda_i=0,\lambda_j=0} = A^{(2),ij} + A^{(2),ji} $$

$$ \left.\frac{\partial^3 A}{\partial\lambda_i\lambda_j\lambda_k}\right|_{\lambda_i=0,\lambda_j=0,\lambda_k=0} = A^{(3),ijk} + A^{(3),jki} + A^{(3),kij} + A^{(3),ikj} + A^{(3),kji} + A^{(3),jik} $$

An alternative way to do this is to symmetrize the equations, and then solve only for one component.

2n+1 theorem

Electric field

For the case of a uniform static electric field $\vec{u}$

$$ V_{ext}(u_i)=u_ir^i\, $$

Current density

$$ j=\frac1{2i}\sum_{k}\left(\psi_k^{*}\nabla\psi_k-\psi_k\nabla\psi_k^{*}\right) $$

$$ j^{(1)}= \frac1{2i}\sum_{m}\left( \psi_k^{(0)*}\nabla\psi_k^{(1)}-\psi_k^{(0)}\nabla\psi_k^{(1)*}+ \psi_k^{(1)*}\nabla\psi_k^{(0)}-\psi_k^{(1)}\nabla\psi_k^{(0)*} \right) $$

ELF

The Electron Localization Function

$$ D(r)=\sum_i\left|\nabla\psi(r)\right|^2 -\frac14\frac{\left|\nabla\rho(r)\right|^2}{\rho(r)}+\frac{j^2(r)}{\rho(t)} $$

The normalization is

$$ f_{ELF}=\frac1{1+\left[D(r)/D_0(r)\right]^2} $$

with

$$ D_0(r)=\frac35\left(6\pi^2\right)\rho^{5/3} $$

Putting our expansion to first order:

$$ D(r)=\sum_i\left|\nabla\psi^{(0)}(r)+\nabla\psi^{(1)}\right|^2 -\frac14\frac{\left|\nabla\rho^{(0)}(r)+\nabla\rho^{(1)}(r)\right|^2}{\rho^{(0)}(r)+\rho^{(1)}(r)} +\frac{\left(j^{(0)}+j^{(1)}\right)^2}{\rho^{(0)}+\rho^{(1)}} $$

$$ D(r)=\sum_i\left[ \nabla\psi^{(0)*}(r)\cdot\nabla\psi^{0}(r)+ \nabla\psi^{(1)}(r)\cdot\nabla\psi^{0}(r)+ \nabla\psi^{(0)*}(r)\cdot\nabla\psi^{1}(r) \right] $$

$$ -\frac14 \left(\frac1{\rho^{(0)}(r)}-\frac{\rho^{(1)}(r)}{\left[\rho^{(0)}(r)\right]^2}\right) \left[ \nabla\rho^{(0)}(r)\cdot\nabla\rho^{(0)}(r)+ \nabla\rho^{(1)}(r)\cdot\nabla\rho^{(0)}(r)+ \nabla\rho^{(0)}(r)\cdot\nabla\rho^{(1)}(r) \right] $$

$$ +\frac{\left(j^{(0)}\right)^2+2j^{(0)}j^{(1)}}{\rho^{(0)}} -\frac{\rho^{(1)}\left(j^{(0)}\right)^2}{\left(\rho^{(0)}\right)^2} $$

Identifying the form of the ELF in the zeroth-order terms, we get the first-order pertubation ELF

$$ D^{(1)}(r)= \sum_i\left[ \nabla\psi^{(1)*}\cdot\nabla\psi^{0}+ \nabla\psi^{(0)*}\cdot\nabla\psi^{1} \right] -\frac12 \frac{ \nabla\rho^{(0)}} {\rho^{(0)}(r)}\cdot\nabla\rho^{(1)} +\frac14\frac{\left|\nabla\rho^{(0)}\right|^2}{\left(\rho^{(0)}\right)^2}\rho^{(1)} +2\frac{j^{(0)}j^{(1)}}{\rho^{(0)}} -\frac{\rho^{(1)}\left(j^{(0)}\right)^2}{\left(\rho^{(0)}\right)^2} $$

And the normalization:

There is a contribution from the change in $D_0$:

$$ D_0^{(1)}=6\pi^2\left(\rho^{(0)}\right)^{2/3}\rho^{(1)} $$

and with this

$$ f_{ELF}^{(1)}= -2\left(f_{ELF}^{(0)}\right)^2 \frac{D^{(0)}}{D_0^{(0)}} \left(\frac{D^{(1)}}{D^{(0)}}-\frac{D^{(0)}}{D_0^{(0)}}\frac{D^{(1)}_0}{D_0^{(0)}}\right) $$

Dynamic Response

We will consider a time-independent ground state system described by a Kohn-Sham equation

$$ (1) \ H^{(0)} |\psi^{(0)}_k> = \epsilon^{(0)}_k | \psi^{(0)}_k>\ , $$

and we will apply a perturbative potential $V_{\omega}(\vec{r},t)$ which is periodic in time with frequency $\omega$ (and period $T=2\pi/\omega$).

In particular we will consider

$$ V_{\omega}=u_lr^l\left(ae^{i\omega{t}}+be^{-i\omega{t}}\right)\ , $$

where $a$ and $b$ are some constant coefficients.

As the perturbed system is time-dependent, we must use the TDDFT Kohn-Sham equations

$$ (1) \ H |\psi_k> = i\frac{\partial}{\partial{t}}| \psi_k>\ , $$

Due to the form of the potential, we can take the following first-order expansion of the wavefunctions

$$ |\psi_k(t)> = e^{-i\epsilon_k{t}}\left[|\psi^{(0)}_k>+u_je^{i\omega{t}}|\psi_k^{(+1),j}>+u_je^{-i\omega{t}}|\psi_k^{(-1),j}>\right]\ . $$

for the density

$$ \rho(\vec{r},t)=\rho^{(0)}(\vec{r})+u_je^{i\omega{t}}\rho^{(+1),j}(\vec{r})+u_je^{-i\omega{t}}\rho^{(-1),j}(\vec{r}) $$

and for the hamiltonian (including the perturbative potential)

$$ H=H^{(0)}+H^{(1)} (u_je^{i\omega{t}}\rho^{(+1),j}(\vec{r})+u_je^{-i\omega{t}}\rho^{(-1),j}(\vec{r}))+u_lr^l\left(ae^{i\omega{t}}+be^{-i\omega{t}}\right)\, $$

where the $H^{(1)}$ operator is the variation of the hamiltonian with respect to the density, and typically includes the Exchange-Correlation and Hartree terms.

$$ H=H^{(0)} +u_je^{ i\omega{t}}\left(H^{(1)}\rho^{(+1),j}(\vec{r})+ar^j\right) +u_je^{-i\omega{t}}\left(H^{(1)}\rho^{(-1),j}(\vec{r})+br^j\right)\ , $$

$$ H=H^{(0)} +u_je^{ i\omega{t}}V^{(+1),j} +u_je^{-i\omega{t}}V^{(-1),j}\, $$

Density

From the expansion for the wavefunctions we can calculate the variation of the density

$$ \rho(\vec{r},t)=\sum_k^{N/2}\psi_k^*\psi_k=\sum_k\left[\psi_k^{*(0)}\psi^{(0)}_k+ u_je^{i\omega{t}}\psi_k^{*(0)}\psi_k^{(+1),j}+ u_je^{-i\omega{t}}\psi_k^{*(0)}\psi_k^{(-1),j}+ u_je^{-i\omega{t}}\psi_k^{*(+1),j}\psi_k^{(0)}+ u_je^{i\omega{t}}\psi_k^{*(-1),j}\psi_k^{(0)}\right] +O(u^2) $$

This gives us

$$ \rho^{(+1),j}=\sum_k^{N/2} \left( \psi_{k}^{*(0)} \psi_{k}^{(+1),j} + \psi_{k}^{*(-1),j} \psi_{k}^{(0)} \right) $$

and

$$ \rho^{(-1),j}=\sum_k^{N/2}\left[ \psi_k^{*(0)}\psi_k^{(-1),j}+\psi_k^{*(+1),j}\psi_k^{(0)}\right] . $$

But we can see that

$$ \rho^{*(-1),j}=\rho^{(+1),j} \ . $$

Using this, we can simplify equation (2):

$$ \rho=\rho^{(0)}+u_je^{i\omega{t}}\rho^{(+1),j}+\left(u_je^{i\omega{t}}\rho^{(+1),j}\right)^* $$

and finally

$$ \rho=\rho^{(0)}+2u_j\cos(\omega{t})\Re\left(\rho^{(+1),j}\right)-2u_j\sin(\omega{t})\Im\left(\rho^{(+1),j}\right)\ . $$

Dynamic current density

$$ J(\vec{r},t)=J^{(0)}(\vec{r})+u_je^{i\omega{t}}J^{(+1),j}(\vec{r})+u_je^{-i\omega{t}}J^{(-1),j}(\vec{r}) $$

$$ J^{(+1),j}=\frac{1}{2i}\sum_{k}\left[ \psi_k^{*(0)}\nabla\psi^{(+1)}+\psi_k^{*(-1)}\nabla\psi_k^{(0)} -\psi_k^{(0)}\nabla\psi^{*(-1)}-\psi_k^{(+1)}\nabla\psi_k^{*(0)} \right] $$

$$ J^{(-1),j}=\frac{1}{2i}\sum_{k}\left[ \psi_k^{*(0)}\nabla\psi^{(-1)}+\psi_k^{*(+1)}\nabla\psi_k^{(0)} -\psi_k^{(0)}\nabla\psi^{*(+1)}-\psi_k^{(-1)}\nabla\psi_k^{*(0)} \right] $$

$$ J^{*(+1),j}=J^{(-1),j}\, $$

Equations for the perturbations

If we take the expansions for the wavefunctions and the hamiltonian and we put them into the TDDFT equation we obtain the equation that the variations obey:

$$ (H^{(0)}-\epsilon_k^{(0)}\pm\omega)|\psi_k^{(\pm1),l}>=-V^{(\pm1),l}|\psi_k^{(0)}>\ . $$

This time there is no projector, so the perturbation has a component onto the occupied space.

As the operator in this equation is the same ground state hamiltonian, we can write the solution in terms of the eigenfunctions (occupied and unoccupied) of the ground state

$$ |\psi_k^{(\pm1),l}>=-\sum_{m=N/2+1}^{\infty}\frac{<\psi_m^{(0)}|V^{(\pm1),l}|\psi_k^{(0)}>}{\epsilon_m-\epsilon_k\pm\omega}|\psi_m^{(0)}>\ , $$

altough to do this numerically we would have to calculate the unoccupied orbitals.

Normalization

One of the conditions we have to ask is that the perturbed wavefunctions should be normalized; this is satisfied if

$$ <\psi_k^{(+1),l}|\psi_k^{(0)}>+<\psi_k^{(0)}|\psi_k^{(-1),l}>=0 $$

using the exact solution for $|\psi_k^{(\pm1),l}>$ we know that

$$ <\psi_k^{(0)}|\psi_k^{(\pm1),l}>=-\frac{<\psi_k^{(0)}|V^{(\pm1),l}|\psi_k^{(0)}>}{\pm\omega} $$

Using this, the normalization condition becomes

$$ \frac{<\psi_k^{(0)}|V^{*(+1),l}|\psi_k^{(0)}>}{-\omega}+\frac{<\psi_k^{(0)}|V^{(-1),l}|\psi_k^{(0)}>}{\omega}=0 $$

$$ <\psi_k^{(0)}|V^{(-1),l}-V^{*(+1),l}|\psi_k^{(0)}>=0 $$

$$ <\psi_k^{(0)}|H^{(1),l}\rho^{(-1),l}+br^l-H^{(1),l}\rho^{*(+1)}-a^*r^l|\psi_k^{(0)}>=0 $$

$$ (b-a^*)<\psi_k^{(0)}|r^l|\psi_k^{(0)}>=0 $$ for k in the occupied space.

So $b=a^*$ which implies that the external potential must be real.

Component of the perturbation on the occupied space

The perturbation can be separated in two perpendicular components

$$ |\psi_k^{(\pm1),l}>=|\psi_k^{u\,(\pm1),l}>+|\psi_k^{o\,(\pm1),l}> $$

The part in the occupied space is

$$ |\psi_k^{o\,(\pm1),l}>=-\sum_m^{N/2}\frac{<\psi_m^{(0)}|V^{(\pm1),l}|\psi_k^{(0)}>}{\epsilon_m-\epsilon_k\pm\omega}|\psi_m^{(0)}>\ , $$

We will now calculate the contribution to the density perturbation of this part, we have that

$$ \rho^{o\,(+1),l}=\sum_k^{N/2}\left[\psi_k^{*(0)}\psi_k^{o\,(+1),l}+\psi_k^{*\,o\,(-1),l}\psi_k^{(0)}\right]\ , $$

using the explicit expresion

$$ \rho^{o\,(+1),l}=\sum_k^{N/2}\sum_m^{N/2}\left[ \psi_k^{*(0)}\frac{<\psi_m^{(0)}|V^{(+1),l}|\psi_k^{(0)}>}{\epsilon_m-\epsilon_k+\omega}\psi_m^{(0)}+ \frac{<\psi_k^{(0)}|V^{*(-1),l}|\psi_m^{(0)}>}{\epsilon_m-\epsilon_k-\omega}\psi_m^{*(0)} \psi_k^{(0)}\right]\ , $$

using that

$$ V^{*(-1),l}=H^{1}\rho^{*(-1)}+b^*r=H^{1}\rho^{(+1)}+ar=V^{(+1),l}\, $$

and swapping the indexes in the second term

$$ \rho^{o\,(+1),l}=\sum_k^{N/2}\sum_m^{N/2}\left[ \frac{<\psi_m^{(0)}|V^{(+1),l}|\psi_k^{(0)}>}{\epsilon_m-\epsilon_k+\omega}\psi_k^{*(0)}\psi_m^{(0)}+ \frac{<\psi_m^{(0)}|V^{(+1),l}|\psi_k^{(0)}>}{-\epsilon_m+\epsilon_k-\omega}\psi_k^{*(0)}\psi_m^{(0)}\right]\ , $$

so finally

$$ \rho^{o\,(+1),l}=0\ . $$

This means that at least for the first order, the part in the occupied space doesn’t contribute to the density and we can put a projector in the RHS of the Steinheimer equation. This is not necessarily true if we want to calculate other things that require the perturbations of the orbitals and not the perturbations of the density.

So finally the equations we have to solve are

$$ (H^{(0)}-\epsilon_k^{(0)}\pm\omega)|\psi_k^{(\pm1),l}>=-P_cV^{(\pm1),l}|\psi_k^{(0)}>\ . $$

Polarizability

We will now calculate the polarizability in terms of the perturbation of the density, first we take the dipole moment and we replace the expresion for the density

$$ p_l(t)=\int\rho(\vec{r},t)r_l\,d^3r= \int\rho^{(0)}(\vec{r})r_l\,d^3r+ u_m\int\left[2\cos(\omega{t})\Re\left(\rho^{(+1),m}(\vec{r})\right) -2\sin(\omega{t})\Im\left(\rho^{(+1),m}(\vec{r})\right)\right]r_l\,d^3r\ , $$

the second term can be identified as the time dependent polarizability

$$ \alpha_{lm}(t)=\int\left[2\cos(\omega{t})\Re\left(\rho^{(+1),m}(\vec{r})\right) -2\sin(\omega{t})\Im\left(\rho^{(+1),m}(\vec{r})\right)\right]r_l\,d^3r\ . $$

Now, we take the fourier transform, because this is the frequency-dependent polarizability we are interested in

$$ \alpha_{lm}(\omega)=\int_0^T\frac{dt}{T}e^{i\omega{t}}\int\left[2\cos(\omega{t})\Re\left(\rho^{(+1),m}(\vec{r})\right) -2\sin(\omega{t})\Im\left(\rho^{(+1),m}(\vec{r})\right)\right]r_l\,d^3r\ . $$

The time integral can be done explicity giving

$$ \alpha_{lm}(\omega)=\int\left[\Re\left(\rho^{(+1),m}(\vec{r})\right) -i\Im\left(\rho^{(+1),m}(\vec{r})\right)\right]r_l\,d^3r $$

and finally this can be written as

$$ \alpha_{lm}(\omega)=\int\rho^{*(+1),m}(\vec{r})r_l\,d^3r\ . $$