# Electronic polarizability in Hohenber and Kohn DFT paper

I try to reproduce Eq.(31) in Hohenberg and Kohn's paper Inhomogeneous Electron Gas (Phys. Rev. 136, B864 (1964)). My understanding to this model is that the external positive charge background(Eq.(27) in the paper) plays the rule as a perturbation to the electron density and in the Hamiltonian we have a perturbation term:

$$\delta H = -\int \frac{n_{ext}(\boldsymbol{r^{'}})}{|\boldsymbol{r}-\boldsymbol{r^{'}}|} \text d^3\boldsymbol{r^{'}}.$$

What I did is:

Firstly, I write $$\delta H$$ into second quantization form and do Fourier transformation(the definition of $$a(\boldsymbol{q})$$ and $$\rho_{\boldsymbol{q}}$$ are given by Eq.(27) and Eq.(30) in the paper): $$\delta H = -\int \frac{n_{ext}(\boldsymbol{r^{'}})}{|\boldsymbol{r}-\boldsymbol{r^{'}}|} \hat{\psi}^{\dagger}(\boldsymbol{r})\hat{\psi}(\boldsymbol{r})\text d^3\boldsymbol{r^{'}}\text d^3\boldsymbol{r}\\ =\sum_{\boldsymbol{q},\boldsymbol{k}}\frac{4\pi a(\boldsymbol{q})}{\boldsymbol{q}}\hat{c}^{\dagger}_{\boldsymbol{k}-\boldsymbol{q}}\hat{c_{\boldsymbol{k}}}\\ =\sum_{\boldsymbol{q}}\frac{4\pi a(\boldsymbol{q})}{\boldsymbol{q}}\rho_{\boldsymbol{q}}$$

Then, I also express the electron charge density operator in second quantization momentum operator: $$\hat{n}(\boldsymbol{r}) = \hat{\psi}^{\dagger}(\boldsymbol{r})\hat{\psi}(\boldsymbol{r}) = \sum_{\boldsymbol{q^{'}}}e^{-i\boldsymbol{q^{'}}\boldsymbol{r}}\rho_{\boldsymbol{q^{'}}}.$$

Then I write down the 1st order perturbation of the ground state $$|0\rangle$$: $$|0^{(1)}\rangle = \sum_{i\neq 0} \frac{\langle i|\delta H|0\rangle}{E_0-E_i}$$

Finally, I calculate the expectation value of the density operator for the perturbed ground state $$|0\rangle + |0^{(1)}\rangle$$ to first order. The perturbation to the electron density is: $$\sum_{n\neq 0}\frac{\langle i|\delta H |0^{(1)}\rangle\langle 0^{(1)}|\hat{n}|i\rangle + \langle i|\hat{n} |0^{(1)}\rangle\langle 0^{(1)}|\delta H|i\rangle}{E_0-E_i}......(a)$$

Substitute the expression of $$\delta H$$ and $$\hat{n}$$ to the equation (a), I found $$b_1(\boldsymbol{q})$$ in the paper can be written as: $$b_1(\boldsymbol{q^{'}})=-\sum_{\boldsymbol{q}} \frac{4\pi a(\boldsymbol{q})}{\boldsymbol{q}} \sum_{i\neq 0}\frac{\langle i |\rho_{\boldsymbol{q}}|0\rangle\langle 0|\rho_{\boldsymbol{q^{'}}}|i\rangle + \langle i |\rho_{\boldsymbol{q^{'}}}|0\rangle\langle 0|\rho_{\boldsymbol{-q}}|i\rangle}{E_0-E_i}$$

I stopped here as I cannot find any way to go from equation above to Eq.(31) in the paper. If I can show that $$\sum_{i\neq 0}\frac{\langle i |\rho_{\boldsymbol{q}}|0\rangle\langle 0|\rho_{\boldsymbol{q^{'}}}|i\rangle}{E_0-E_i} = 0$$ except $$q = -q^{'}$$, then the problem resolved. However, I don't think this is obvious.

In first quantization: \begin{align} \delta\hat{H}&=-\int_{V}\frac{n_{ext}(\mathbf{r}')}{\vert \mathbf{r}-\mathbf{r}'\vert}\mathrm{d}\mathbf{r'}\equiv\hat{f}^{(1)}(\mathbf{r}) \end{align} In second quantization it reads: \begin{align} \delta\hat{H}&=\int_V\mathrm{d}\mathbf{r}\hat{\Psi}^\dagger(\mathbf{r})\hat{f}^{(1)}(\mathbf{r})\hat{\Psi}(\mathbf{r})=\int_V\mathrm{d}\mathbf{r}\sum_{\mathbf{k}}\frac{1}{\sqrt{V}}e^{-i\mathbf{k\cdot\mathbf{r}}}c^\dagger_{\mathbf{k}}~\hat{f}^{(1)}(\mathbf{r})~\sum_{\mathbf{k}'}\frac{1}{\sqrt{V}}e^{i\mathbf{k'\cdot\mathbf{r}}}c_{\mathbf{k}'}\\&=\sum_{\mathbf{k},\mathbf{k}'}c^\dagger_{\mathbf{k}}c_{\mathbf{k}'}\frac{1}{V}\int_V\mathrm{d}\mathbf{r}e^{-i (\mathbf{k}-\mathbf{k}')\cdot\mathbf{r}}\hat{f}^{(1)}(\mathbf{r})\\&=\sum_{\mathbf{k},\mathbf{k}'}c^\dagger_{\mathbf{k}}c_{\mathbf{k}'}\frac{1}{V}\int_V\mathrm{d}\mathbf{r}e^{-i (\mathbf{k}-\mathbf{k}')\cdot\mathbf{r}}\int_{V}\mathrm{d}\mathbf{r}'\sum_{\mathbf{q}}\bigg(-\frac{\lambda}{V}\bigg)\frac{a(\mathbf{q})e^{-i\mathbf{q}\cdot\mathbf{r}'}}{\vert\mathbf{r}-\mathbf{r}'\vert}\\&= -\frac{\lambda}{V^2}\sum_{\mathbf{k},\mathbf{k}',\mathbf{q}}c^\dagger_{\mathbf{k}}c_{\mathbf{k}'}a(\mathbf{q})\int_{V\times V}\mathrm{d}\mathbf{r}\mathrm{d}\mathbf{r}'\frac{e^{-i(\mathbf{k}-\mathbf{k}')\cdot\mathbf{r}}e^{-i\mathbf{q\cdot\mathbf{r}'}}}{\vert\mathbf{r}-\mathbf{r}'\vert} \end{align} Next, we define $$$$\begin{bmatrix} \mathbf{t}\\\mathbf{T} \end{bmatrix}=\underbrace{\frac{1}{\sqrt{2}}\begin{pmatrix} 1 & -1 \\1 &1 \end{pmatrix}}_{\det = 1}\begin{bmatrix} \mathbf{r}\\\mathbf{r}' \end{bmatrix}$$$$ So that \begin{align} \delta\hat{H}&= -\frac{\lambda}{V^2}\sum_{\mathbf{k},\mathbf{k}',\mathbf{q}}c^\dagger_{\mathbf{k}}c_{\mathbf{k}'}a(\mathbf{q})\frac{1}{2}\int_{V\times V}\mathrm{d}\mathbf{t}\mathrm{d}\mathbf{T}\frac{e^{-\frac{i}{\sqrt{2}}(\mathbf{k}-\mathbf{k}')\cdot(\mathbf{t}+\mathbf{T})}e^{-\frac{i}{\sqrt{2}}\mathbf{q\cdot(\mathbf{T}-\mathbf{t})}}}{\sqrt{2}\vert\mathbf{t}\vert} \\&= -\frac{\lambda}{V^2}\sum_{\mathbf{k},\mathbf{k}',\mathbf{q}}c^\dagger_{\mathbf{k}}c_{\mathbf{k}'}a(\mathbf{q})\frac{1}{2}\int_{V\times V}\mathrm{d}\mathbf{t}\mathrm{d}\mathbf{T}\frac{e^{-\frac{i}{\sqrt{2}}(\mathbf{k}-\mathbf{k}'-\mathbf{q})\cdot\mathbf{t}}e^{-\frac{i}{\sqrt{2}}((\mathbf{k}-\mathbf{k}')+\mathbf{q})\cdot\mathbf{T}}}{\sqrt{2}\vert\mathbf{t}\vert} \end{align} Short moment to remember that \begin{align} \frac{1}{2\pi}\int_{-\infty}^{\infty}e^{-i\xi x}\mathrm{d}x=\delta(\xi) \end{align} to perform the integral in $$\mathbf{T}$$ yielding \begin{align} \delta\hat{H}&= -\frac{\lambda}{V^2}\sum_{\mathbf{k},\mathbf{k}',\mathbf{q}}c^\dagger_{\mathbf{k}}c_{\mathbf{k}'}a(\mathbf{q})\frac{1}{2}\int_{V}\mathrm{d}\mathbf{t}\frac{e^{-\frac{i}{\sqrt{2}}(\mathbf{k}-\mathbf{k}'-\mathbf{q})\cdot\mathbf{t}}\delta(\sqrt{2}^{-1}(\mathbf{k}-\mathbf{k}'+\mathbf{q}))}{\sqrt{2}\vert \mathbf{t}\vert} ~8\pi^3 \\&=-\frac{\lambda}{V^2}\sum_{\mathbf{k},\mathbf{k}',\mathbf{q}}c^\dagger_{\mathbf{k}}c_{\mathbf{k}'}a(\mathbf{q})\frac{1}{2}\int_{V}\mathrm{d}\mathbf{t}\frac{e^{-\frac{i}{\sqrt{2}}(\mathbf{k}-\mathbf{k}'-\mathbf{q})\cdot\mathbf{t}}\delta(\mathbf{k}-\mathbf{k}'+\mathbf{q})}{\sqrt{2}^{-3}\sqrt{2}\vert \mathbf{t}\vert}8\pi^3\\&\Rightarrow \frac{1}{V}\sum_{\mathbf{k}}\mapsto\int\frac{\mathrm{d}\mathbf{k}}{8\pi^3} \\&=-\frac{\lambda}{V}\sum_{\mathbf{k}',\mathbf{q}}a(\mathbf{q})\frac{1}{2}\int_{V}\mathrm{d}\mathbf{t}\int\frac{\mathrm{d}\mathbf{k}}{8\pi^3}c^\dagger_{\mathbf{k}}c_{\mathbf{k}'}\frac{e^{-\frac{i}{\sqrt{2}}(\mathbf{k}-\mathbf{k}'-\mathbf{q})\cdot\mathbf{t}}\delta(\mathbf{k}-(\mathbf{k}'-\mathbf{q}))}{\vert \mathbf{t}\vert}32\pi^3\\&=-\frac{\lambda}{V}\sum_{\mathbf{k}',\mathbf{q}}c^\dagger_{\mathbf{k}'-\mathbf{q}}c_{\mathbf{k}'}a(\mathbf{q})\int_{V}\mathrm{d}\mathbf{t}\frac{e^{\frac{i}{\sqrt{2}}(2\mathbf{q})\cdot\mathbf{t}}}{\vert \mathbf{t}\vert}4 \\&=-\frac{\lambda}{V}\sum_{\mathbf{k}',\mathbf{q}}c^\dagger_{\mathbf{k}'-\mathbf{q}}c_{\mathbf{k}'}a(\mathbf{q})\frac{1}{2}\int_{V}\mathrm{d}\mathbf{t}\frac{e^{i\sqrt{2}\mathbf{q}\cdot\mathbf{t}}}{\vert\mathbf{t}\vert}4\\&=-\frac{\lambda}{V}\sum_{\mathbf{k}',\mathbf{q}}c^\dagger_{\mathbf{k}'-\mathbf{q}}c_{\mathbf{k}'}a(\mathbf{q})\frac{16\pi}{4\vert \mathbf{q}\vert^2}=-\frac{\lambda}{V}\sum_{\mathbf{q}}\rho_\mathbf{q}a(\mathbf{q})\frac{4\pi}{\vert \mathbf{q}\vert^2}=\delta\hat{H} \end{align} For the last line see e.g. here

Next consider (see, e.g., Gasiorowicz) $$$$\vert\Psi_0^{appr}\rangle=\vert\psi_0\rangle+\sum_{l\neq 0}\frac{\langle\psi_l\vert\delta\hat{H}\vert\psi_0\rangle}{\epsilon_l-\epsilon_0}\vert\psi_l\rangle$$$$

The density operator in second quantization is $$$$\hat{n}(\mathbf{r})=\hat{\Psi}^\dagger(\mathbf{r})\hat{\Psi}(\mathbf{r})=\frac{1}{V}\sum_{\mathbf{k}}\sum_{\mathbf{k}'}e^{-i\underbrace{(\mathbf{k}-\mathbf{k}')}_{=: \mathbf{q}}\mathbf{r}}c^\dagger_\mathbf{k}c_\mathbf{k'}$$$$ which implies $$\mathbf{k}=\mathbf{k}'+\mathbf{q}$$ $$$$\hat{n}(\mathbf{r})=\hat{\Psi}^\dagger(\mathbf{r})\hat{\Psi}(\mathbf{r})=\frac{1}{V}\sum_{\mathbf{q}}\sum_{\mathbf{k}'}e^{-i\mathbf{q}\mathbf{r}}c^\dagger_{\mathbf{k}'+\mathbf{q}}c_\mathbf{k'}=\frac{1}{V}\sum_{\mathbf{q}}e^{-i\mathbf{q}\mathbf{r}}\rho_{-\mathbf{q}}$$$$

So next we compute $$\langle\Psi_0^{approx}\vert\hat{n}(\mathbf{r})\vert\Psi_0^{approx}\rangle$$ to leading order in $$\lambda$$: \begin{align} \langle\Psi_0^{approx}\vert\hat{n}(\mathbf{r})\vert\Psi_0^{approx}\rangle&=\langle\psi_0\vert\hat{n}(\mathbf{r})\vert\psi_0\rangle+2\sum_{l\neq 0}\langle\psi_l\vert\frac{\langle\psi_l\vert\delta\hat{H}\vert\psi_0\rangle}{\epsilon_l-\epsilon_0}\cdot\hat{n}(\mathbf{r})\cdot\vert\psi_0\rangle \\&=n_0(\mathbf{r})+2\sum_{l\neq 0}\frac{\langle\psi_0\vert\delta\hat{H}\vert\psi_l\rangle\langle\psi_l\vert\hat{n}(\mathbf{r})\vert\psi_0\rangle}{\epsilon_l-\epsilon_0} \\&=n_0(\mathbf{r})+2\sum_{l\neq 0}\frac{1}{V}\sum_{\mathbf{q}}e^{-i\mathbf{q}\mathbf{r}}\frac{\langle\psi_0\vert\delta\hat{H}\vert\psi_l\rangle\langle\psi_l\vert\rho_{-\mathbf{q}}\vert\psi_0\rangle}{\epsilon_l-\epsilon_0} \\&=n_0(\mathbf{r})-2\sum_{l\neq 0}\frac{\lambda}{V^2}\sum_{\mathbf{q}}\sum_{\mathbf{q}'}\frac{4\pi}{q'^2}e^{-i\mathbf{q}\mathbf{r}}a(\mathbf{q}')\frac{\langle\psi_0\vert\rho_{\mathbf{q}'}\vert\psi_l\rangle\langle\psi_l\vert\rho_{-\mathbf{q}}\vert\psi_0\rangle}{\epsilon_l-\epsilon_0} \\&= n_0(\mathbf{r})-2\sum_{l\neq 0}\frac{\lambda}{V^2}\sum_{\mathbf{q}}\frac{4\pi}{q^2}e^{-i\mathbf{q}\mathbf{r}}a(\mathbf{q})\frac{\langle\psi_0\vert\rho_{\mathbf{q}}\vert\psi_l\rangle\langle\psi_l\vert\rho_{-\mathbf{q}}\vert\psi_0\rangle}{\epsilon_l-\epsilon_0} \end{align}

The key in the last step is that the excited state $$\vert\psi_l\rangle$$ appears in both factors. Therefore the excitement operator $$\hat{\rho}_{\mathbf{q}'}$$ yields nonzero values only if $$\mathbf{q}=\mathbf{q'}$$.

In other words, the term to the right is exciting one particle to an orbital distant $$+\mathbf{q}$$, whereas the term to the left is de-exciting a particle to an orbital distant $$-\mathbf{q}'$$; on the other hand, the index $$l$$ for the excited state is “locked” (equal in both terms) and therefore the only case in which neither one of the two inner products is zero is if $$\mathbf{q}=\mathbf{q}'$$.

Finally, we notice that it seems that we are off HK's result by volume factor. However, Mermin, who quotes them, seems to have the same prefactor in agreement with what derived here.