Linear Response

From OctopusWiki

Jump to: navigation, search

WARNING, This page is work in progress

1 Static Response
2 Dynamic Response

[edit]

Static Response

We will consider a Kohn-Sham equation

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

and a perturbative 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 dependency in lambda, for the zeroth order we get the fundamental state equation.

[edit]

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)|\psi^{(1),i}_k> = -P_cH^{(1),i}|\psi^{(0)}_k>$

where $P c$ is a proyector onto the unoccupied space. As the equation is linear, this implies that the perturbations only have components in the unoccupied space.

[edit]

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\ .$

[edit]

Relation to the derivatives

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

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.

[edit]

2n+1 theorem

[edit]

Electric field

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

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

[edit]

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)$

[edit]

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)$

[edit]

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 $ω$ (and period $T = 2π / ω$ ).

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 tipically 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}\,$

[edit]

Density

From the expansion for the wavefunctions we can calculate the variations 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)\ .$

[edit]

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}\,$

[edit]

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 or 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 fundamental 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.

[edit]

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\qquad \mbox{for k in the occupied space}$

$<\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$

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

[edit]

Component of the perturbation on the occupied space

The perturbation can be separated in two perperdincular 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 swaping 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 proyector in the RHS of the Steinheimer equation. This is not necessarily true is we want to calculate other things that requiere 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)}>\ .$

[edit]

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 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\ .$