Electron Gas Susceptibility LaTeX Source

\documentclass[aps,12pt,prb]{revtex4}

\begin{document}

\title{Current-Current Response of Electron Gas}

\author{R. B. Laughlin}

\affiliation{Stanford}

\date{17 May 08}

\maketitle

\section{Density Response Definition}

\noindent
Taking $\omega$ to be shorthand for $\omega + i \eta$, we have

\begin{displaymath}
<\rho_q \rho_{-q}>_{\omega} = 2 \sum_{|k| < k_f \; |k+q| > k_f}
\biggl\{ \frac{1}{\hbar \omega - (E_{k+q} - E_k)}
+ \frac{1}{- \hbar \omega - (E_{k+q} - E_k)} \biggr\}
\end{displaymath}

\begin{displaymath}
= 2 \sum_{|k| < k_f}
\biggl\{ \frac{1}{\hbar \omega - (E_{k+q} - E_k)}
+ \frac{1}{- \hbar \omega - (E_{k+q} - E_k)} \biggr\}
\end{displaymath}

\noindent
This last step is a trick.  You can formally sum over transitions for 
which $|k+q| < k_f$ because the anti-resonant contribution for such 
transitions exactly cancels the resonant contributions.

\section{Density Response Computation}

\noindent
Taking

\begin{displaymath}
E_k = \frac{\hbar^2 k^2}{2m}
\; \; \; \; \; \; \; \; \; 
E_f = \frac{\hbar^2 k_f^2}{2m}
\; \; \; \; \; \; \; \; \; 
x = \frac{q}{k_f}
\; \; \; \; \; \; \; \; \; 
y = \frac{\hbar \omega}{E_f}
\end{displaymath}

\noindent
we then have

\begin{displaymath}
<\rho_q \rho_{-q}>_{\omega} = 2 (\frac{L}{2 \pi})^3
\int_{|k| < k_f}
\biggl\{ \frac{1}{\hbar \omega - (E_{k+q} - E_k)}
+ \frac{1}{- \hbar \omega - (E_{k+q} - E_k)} \biggr\} \; d^3k
\end{displaymath}

\begin{displaymath}
= 2 (\frac{L}{2 \pi})^3 \frac{2 \pi k_f^3}{E_f} \int_{-1}^1 d\mu
\int_0^1 \xi^2 d\xi
\biggl\{ \frac{1}{y - (x^2 + 2x\xi \mu)}
+ \frac{1}{- y - (x^2 + 2 x \xi \mu)} \biggr\}
\end{displaymath}

\begin{displaymath}
= 2 (\frac{L}{2 \pi})^3 \frac{\pi k_f^3}{x E_f}
\int_0^1 \xi d\xi \biggl\{
\ln \biggl[\frac{y - (x^2 - 2x\xi)}{y - (x^2 + 2x\xi)}\biggr] +
\ln \biggl[\frac{-y - (x^2 - 2x\xi)}{-y - (x^2 + 2x\xi)}\biggr] \biggr\}
\end{displaymath}

\begin{displaymath}
= - 2 (\frac{L}{2 \pi})^3 \frac{\pi k_f^3}{E_f} \biggl\{1
+ \biggl[\frac{(y-x^2)^2-(2x)^2}{8x^3} \biggr] \ln
\biggl[ \frac{y - x^2 + 2x}{y - x^2 - 2x} \biggr]
\end{displaymath}

\begin{displaymath}
+ \biggl[\frac{(-y-x^2)^2-(2x)^2}{8x^3} \biggr] \ln
\biggl[ \frac{-y - x^2 + 2x}{-y - x^2 - 2x} \biggr] \biggr\} \; \; \; .
\end{displaymath}

\section{Integration Details}

\noindent
The $\xi$ intergral is elementary.  First we evaluate

\begin{displaymath}
\int_0^1 \xi \ln(y - x^2 +2x\xi) \; d\xi =
\frac{1}{2x} \int_0^1 \biggl\{ (2x\xi + y -x^2)
\ln (2x\xi + y -x^2) 
- (y - x^2) \ln (2x\xi + y -x^2) \biggr\} \; d\xi
\end{displaymath}

\begin{displaymath}
= \frac{1}{4x^2} \biggl\{ \frac{(2x\xi + y -x^2)^2}{2}
\biggl[ \ln (2x\xi + y -x^2) - \frac{1}{2} \biggr] 
\end{displaymath}

\begin{displaymath}
- (y - x^2) (2x\xi + y - x^2)
\biggl[ \ln (2x\xi + y -x^2) - 1 \biggr] \biggr\}
\bigg|_0^1
\end{displaymath}

\begin{displaymath}
= \frac{1}{4x^2} \biggl\{ \frac{(y -x^2 + 2x)^2}{2}
\biggl[ \ln (y -x^2 + 2x) - \frac{1}{2} \biggr] 
- (y - x^2) (y - x^2 +2x) \biggl[ \ln (y -x^2 + 2x) - 1 \biggr]
\end{displaymath}

\begin{displaymath}
+ \frac{(y - x^2)^2}{2}\biggl[ \ln (y -x^2) - \frac{1}{2}
\biggr] \biggr\} \; \; \; .
\end{displaymath}

\noindent
Then we subtract this expression from itself with $x \longrightarrow 
-x$:

\begin{displaymath}
\int_0^1 \xi \ln \biggl[
\frac{y - x^2 +2x\xi}{y - x^2 - 2 x\xi} \biggr] \; d\xi =
\frac{1}{4x^2} \biggl\{
2x(y-x^2) - \biggl[\frac{(y-x^2)^2-(2x)^2}{2} \biggr] \ln
\biggl[ \frac{y - x^2 + 2x}{y - x^2 - 2x} \biggr] \biggr\}
\end{displaymath}

\noindent
Finally we add this expression from itself with $y \longrightarrow 
-y$:

\begin{displaymath}
\int_0^1 \xi \ln \biggl\{
\biggl[ \frac{y - x^2 +2x\xi}{y - x^2 - 2 x\xi} \biggr] +
\biggl[ \frac{- y - x^2 +2x\xi}{- y - x^2 - 2 x\xi} \biggr] \biggr\}
 \; d\xi
\end{displaymath}

\begin{displaymath}
= x
+ \biggl[\frac{(y-x^2)^2-(2x)^2}{8x^2} \biggr] \ln
\biggl[ \frac{y - x^2 + 2x}{y - x^2 - 2x} \biggr]
+ \biggl[\frac{(-y-x^2)^2-(2x)^2}{8x^2} \biggr] \ln
\biggl[ \frac{-y - x^2 + 2x}{-y - x^2 - 2x} \biggr] \; \; \; .
\end{displaymath}

\section{Tests}

\noindent
Let's now do a couple of tests.  First of all, we should obviously have

\begin{displaymath}
\lim_{\omega \rightarrow 0} \lim_{q \rightarrow \infty}
<\rho_q \rho_{-q}>_{\omega}
= - 2 \sum_{|k| < k_f} \frac{2}{E_q} = - \frac{2N}{E_q} \; \; \; .
\end{displaymath}

\noindent
That works out, since

\begin{displaymath}
\lim_{x \rightarrow \infty}
- 2 (\frac{L}{2 \pi})^3 \frac{\pi k_f^3}{E_f} \biggl\{1
+ \biggl[\frac{x^4-4x^2}{4x^3} \biggr] \ln
\biggl[ \frac{x^2 - 2x}{x^2 + 2x} \biggr] \biggr\}
= - 2 (\frac{L}{2 \pi})^3 \frac{\pi k_f^3}{E_f} \frac{8}{3 x^2}
\; \; \; .
\end{displaymath}

\noindent
We also have

\begin{displaymath}
\lim_{x \rightarrow 0}
- 2 (\frac{L}{2 \pi})^3 \frac{\pi k_f^3}{E_f} {\rm Re} \biggl\{1
+ \biggl[\frac{x^4-4x^2}{4x^3} \biggr] \ln
\biggl[ \frac{x^2 - 2x}{x^2 + 2x} \biggr] \biggr\}
= - 4 (\frac{L}{2 \pi})^3 \frac{\pi k_f^3}{E_f} 
\; \; \; .
\end{displaymath}

\noindent
This is indeed consistent with the Thomas-Fermi screening wavevector
$q_{TF}$, per

\begin{displaymath}
\lim_{\omega \rightarrow 0} \lim_{q \rightarrow 0}
- \frac{4 \pi e^2}{L^3 q^2} {\rm Re} (<\rho_q \rho_{-q}>_{\omega}) = 
\frac{2 k_f^3 e^2}{\pi E_f} = q_{TF}^2 \; \; \; .
\end{displaymath}

\section{Branch Cuts}

\noindent
The imaginary part of the Lindhard function is easiest to work out
by common sense.  It is

\begin{displaymath}
{\rm Im} \biggl\{ \ln \biggl[ \frac{y - x^2 + 2x}{y - x^2 - 2x} \biggr]
\biggr\} = 
\left[ \begin{array}{rcl}
- \pi &  \; \; \; \; \; \; \; \; & x^2 - 2x < {\rm Re}(y) < x^2 + 2x \\
0 & & {\rm otherwise}
\end{array} \right]
\end{displaymath}

\begin{displaymath}
{\rm Im} \biggl\{ \ln \biggl[ \frac{- y - x^2 + 2x}{- y - x^2 - 
2x} 
\biggr] \biggr\} = 
\left[ \begin{array}{rcl}
\pi &  \; \; \; \; \; \; \; \; & x^2 - 2x < {\rm Re}(- y) < x^2 + 2x \\
0 & & {\rm otherwise}
\end{array} \right]
\end{displaymath}

\noindent
Note that this rule works for both $x > 2$ and $x < 2$.  In these
formulas, x is real and positive by definition but y can be either 
positive or negative.

\section{Current Response Definition}

\begin{displaymath}
_\omega = 2 \sum_{|k| < k_f \; |k+q| > k_f}
( \frac{\hbar}{m})^2 
(k^\mu + \frac{q^\mu}{2}) (k^\nu + \frac{q^\nu}{2})
\end{displaymath}

\begin{displaymath}
\times \biggl\{ \frac{1}{\hbar \omega - (E_{k+q} - E_k)}
+ \frac{1}{- \hbar \omega - (E_{k+q} - E_k)} \biggr\}
+ \frac{N}{m} \delta^{\mu \nu}
\end{displaymath}

\begin{displaymath}
= 2 \sum_{|k| < k_f}
( \frac{\hbar}{m})^2 
(k^\mu + \frac{q^\mu}{2}) (k^\nu + \frac{q^\nu}{2})
\biggl\{ \frac{1}{\hbar \omega - (E_{k+q} - E_k)}
+ \frac{1}{- \hbar \omega - (E_{k+q} - E_k)} \biggr\}
+ \frac{N}{m} \delta^{\mu \nu}
\end{displaymath}

\noindent
By symmetry the current-current response tensor splits into longitudinal
and transverse parts, per

\begin{displaymath}
_\omega = 
\frac{q^\mu q^\nu}{q^2} _\omega 
+ (\delta^{\mu \nu} - \frac{q^\mu q^\nu}{q^2})
_\omega 
\end{displaymath}

\noindent
We don't have to compute the longitudinal component, since it's
related to the density-density response function through

\begin{displaymath}
q^2 _\omega =
\sum_{\mu \nu} q^\mu _\omega q^\nu = 
\omega^2 <\rho_q \rho_{-q}>_{\omega}
\end{displaymath}

\section{Current Response Computation}

\begin{displaymath}
_{\omega} - \frac{N}{m}
\end{displaymath}

\begin{displaymath}
= 2 (\frac{L}{2 \pi})^3
(\frac{\hbar}{m})^2
\int_{|k| < k_f} \biggl[k^2 - (\frac{k \cdot q}{q})^2 \biggr]
\biggl\{ \frac{1}{\hbar \omega - (E_{k+q} - E_k)}
+ \frac{1}{- \hbar \omega - (E_{k+q} - E_k)} \biggr\} \; d^3k
\end{displaymath}

\begin{displaymath}
= 2 (\frac{L}{2 \pi})^3 \frac{2 \pi k_f^3}{E_f} 
(\frac{\hbar k_f}{m})^2 \int_{-1}^1
(1 - \mu^2) d\mu \int_0^1 \xi^2 d\xi
\biggl\{ \frac{1}{y - (x^2 + 2x\xi \mu)}
+ \frac{1}{- y - (x^2 + 2 x \xi \mu)} \biggr\}
\end{displaymath}

\begin{displaymath}
= 2 (\frac{L}{2 \pi})^3 \frac{\pi k_f^3}{x E_f}
(\frac{\hbar k_f}{m})^2 \int_0^1 \xi d\xi \biggl\{
\biggl[ 1 - (\frac{y - x^2}{2x\xi})^2 \biggr]
\ln \biggl[\frac{y - (x^2 - 2x\xi)}{y - (x^2 + 2x\xi)}\biggr]
\end{displaymath}

\begin{displaymath}
+ \biggl[ 1 - (\frac{-y - x^2}{2x\xi})^2 \biggr]
\ln \biggl[\frac{-y - (x^2 - 2x\xi)}{-y - (x^2 + 2x\xi)}\biggr]
- \frac{2x}{\xi} \biggr\}
\end{displaymath}

\begin{displaymath}
= - 2 (\frac{L}{2 \pi})^3 \frac{\pi k_f^3}{E_f} 
(\frac{\hbar k_f}{m})^2 \biggl\{1
+ \biggl[\frac{(y-x^2)^2-(2x)^2}{8x^3} \biggr] \ln
\biggl[ \frac{y - x^2 + 2x}{y - x^2 - 2x} \biggr]
\end{displaymath}

\begin{displaymath}
+ \biggl[\frac{(-y-x^2)^2-(2x)^2}{8x^3} \biggr] \ln
\biggl[ \frac{-y - x^2 + 2x}{-y - x^2 - 2x} \biggr] + f_x(y) \biggr\} \; 
\; \; .
\end{displaymath}

\section{Special Function}

\noindent
The function

\begin{displaymath}
f_x(y) = \frac{1}{4x^3} \int_0^1 \biggl\{ 
(y-x^2)^2 \ln \biggl[ \frac{y - x^2 + 2x\xi}{y - x^2 -2x\xi} \biggr] +
(-y-x^2)^2 \ln \biggl[ \frac{- y - x^2 + 2x\xi}{-y - x^2 -2x\xi} \biggr]
\biggr\} \frac{d\xi}{\xi} + 2
\end{displaymath}

\noindent 
cannot be expressed in terms of elementary function but is nonetheless 
simple conceptually.  First of all, it is a (continuous) superposition 
of poles in the lower complex plane, since it is, up to contant factors, 
the difference of $_\omega$ and $<\rho_q 
\rho_{-q}>_\omega$, both of which have this property. As a check, we 
check the behavior at large y:

\begin{displaymath}
\lim_{y \rightarrow 0} f_x(y) = \frac{1}{2x^3} \int_0^1 \biggl\{
(y-x^2)^2 \ln \biggl[ (\frac{2x\xi}{y - x^2}) + \frac{1}{9}
(\frac{2x\xi}{y - x^2})^3\biggr]
\end{displaymath}

\begin{displaymath}
+ (-y-x^2)^2 \ln \biggl[ (\frac{2x\xi}{-y - x^2}) + \frac{1}{9}
(\frac{2x\xi}{-y - x^2})^3\biggr] \biggr\} \frac{d\xi}{\xi} + 2
= \frac{4}{9} \biggl[ \frac{1}{y - x^2} + \frac{1}{-y-x^2} \biggr]
\end{displaymath}

\noindent
Thus thus function is a spectrum of oscillators centered around 
frequency $x^2$ and having strength 4/9 independent of x. We can thus 
generate this function numerically by hilbert transforming
its much simpler imaginary part.  The latter is given by

\begin{displaymath}
{\rm Im} f_x(y) = a_x^+(y) + a_x^-(y)
\end{displaymath}

\begin{displaymath}
a_x^+(y) = \frac{1}{4x^3} \left[
\begin{array}{ccl}
\pi (y - x^2)^2 \ln | (y-x^2)/2x | 
 & \; \; \; \; \; \; &
|y - x^2 | < 2x \\
0 & & {\rm otherwise} \end{array} \right]
\end{displaymath}

\begin{displaymath}
a_x^-(y) = \frac{1}{4x^3} \left[
\begin{array}{ccl}
- \pi (- y - x^2)^2 \ln | (-y-x^2)/2x | 
 & \; \; \; \; \; \; &
|-y - x^2 | < 2x \\
0 & & {\rm otherwise} \end{array} \right]
\end{displaymath}

\noindent
The pole strength works out correctly, per

\begin{displaymath}
- \int_{-2x}^{2x} \frac{y^2}{4x^3} \ln \bigg| \frac{y}{2x} \bigg| dy
= -  \frac{1}{2x^3} \biggl\{ \frac{y^3}{3} \biggl[ \ln(\frac{y}{2x}
- \frac{1}{3} \biggr] \biggr\} \bigg|_0^{2x} = \frac{4}{9}
\end{displaymath} 

\begin{figure}
\input{fig}
\caption{Real and imaginary parts of longitudinal and transverse
susceptibilities $\chi_q(\omega)$ versus $y = \hbar \omega/E_f$ for the 
case x=1. The transverse susceptibility is the smoother of the two.}
\end{figure}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Electron Gas Gauge Invariance Relation}

\begin{displaymath}
\sum_{\mu \nu} q^\mu _\omega q^\nu =
\frac{4}{\hbar^2} \sum_{|k| < k_f} 
\frac{ [ (E_{k+q} - E_k)^2 - (\hbar \omega)^2 + (\hbar \omega^2)]
(E_{k+q} - E_k)}{(\hbar \omega)^2 - (E_{k+q} - E_k)^2}
+ \frac{N}{m} q^2
\end{displaymath}

\begin{displaymath}
= - \frac{N}{m} q^2
+ \omega^2 <\rho_q \rho_{-q}>_{\omega}
+ \frac{N}{m} q^2
= \omega^2 <\rho_q \rho_{-q}>_{\omega}
\end{displaymath}

\section{General Gauge Invariance Relation}

\noindent
This follows very generally from conservation of particles, distilled in 
the continuity equation

\begin{displaymath}
[ {\cal H}, \rho_q ] = \hbar \sum_\mu q^\mu j_q^\mu 
\; \; \; \; \; \; \;
\rho_q = \sum_j^N e^{i q \cdot r_j}
\; \; \; \; \; \; \;
j_q^\mu = \sum_j^N \frac{1}{2}
\{ \frac{\hbar}{m} \nabla^\mu, e^{i q \cdot r_j} \}
\end{displaymath}

\noindent
for the specific case of

\begin{displaymath}
{\cal H} = - \sum_j^N \frac{\hbar^2}{2m} \nabla_j^2
+ V(r_1 , ... , r_N)
\end{displaymath}

\noindent
We then have

\begin{displaymath}
\omega^2 <\rho_q \rho_{-q}>_{\omega}
= \omega^2 \sum_x <0|\rho_{-q}|x>  \biggl[
\frac{2 (E_x - E_0)}{(\hbar \omega)^2 - (E_x - E_0)^2} \biggr]
\end{displaymath}

\begin{displaymath}
= \frac{2}{\hbar^2} \sum_x 
<0|\rho_{-q}|x>  \biggl[
\frac{(\hbar \omega)^2 - (E_x - E_0)^2 + (E_x - E_0)^2}
{(\hbar \omega)^2 - (E_x - E_0)^2} \biggr]
\end{displaymath}

\begin{displaymath}
= \frac{2}{\hbar} \sum_\mu q^\mu \biggl\{ \sum_x 
<0|\rho_{-q}|x>  \biggl[ 1 +
\frac{(E_x - E_0)^2}{(\hbar \omega)^2 - (E_x - E_0)^2} \biggr] \biggr\}
\end{displaymath}

\begin{displaymath}
= \frac{2}{\hbar} \sum_\mu q^\mu \biggl\{
\frac{1}{2} <0|[\rho_{-q},j_q^\mu]|0>
- \sum_x 
<0|[{\cal H},\rho_{-q}]|x>  \biggl[
\frac{(E_x - E_0)}{(\hbar \omega)^2 - (E_x - E_0)^2} \biggr] \biggr\}
\end{displaymath}

\begin{displaymath}
= \sum_{\mu \nu} q^\mu \biggl\{ \frac{N}{m} \delta^{\mu \nu}
+ \sum_x 
<0|j_{-q}^\mu |x>  \biggl[
\frac{2(E_x - E_0)}{(\hbar \omega)^2 - (E_x - E_0)^2} \biggr]
\biggr\} q^\nu
\end{displaymath}

\begin{displaymath}
= \sum_{\mu \nu} q^\mu _\omega q^\nu
\end{displaymath}

\end{document}