これは何？
Electronic Structure: Basic Theory and Practical Methods 2nd edition (Richard M. Martin著、
Cambridge University Press、2020)の行間を私なりに埋めたものです。内容を網羅している訳ではありません。
5.2 Noninteracting and Hartree-Fock Approximations
非相互作用近似では、シュレーディンガー方程式は
$$
-\frac{\hbar^2}{2m_e}\nabla^2\psi_i^\sigma(\mathbf{r})=\varepsilon_i^\sigma\psi_i^\sigma(\mathbf{r})
$$
となり、その解は、
$$
\psi_{k} = \left( \frac{1}{\Omega^{1/2}} \right) e^{i \mathbf{k} \cdot \mathbf{r}}, \epsilon_{k} = \frac{\hbar^{2}}{2m_{e}} k^{2}.
$$
である。
Density Matrix
密度行列
$$
\rho(\mathbf{r}, \sigma; \mathbf{r}', \sigma') = \delta_{\sigma, \sigma'} \sum_{i} \psi_{i}^{\sigma *}(\mathbf{r}) f_{i} \psi_{i}^{\sigma}(\mathbf{r}')
$$
にこれを代入すると、
\begin{align*}
\rho(\mathbf{r}, \mathbf{r}')
&=\rho(|\mathbf{r}-\mathbf{r}'|) \\
&=\frac{1}{\Omega}\frac{\Omega}{(2\pi)^3}\int\mathrm{d}\mathbf{k}
f\left(\beta(\varepsilon(\mathbf{k})-\mu)\right)
e^{i\mathbf{k}\mathbf{r}}\cdot e^{-i\mathbf{k}\mathbf{r}'}
\quad\left(\because \sum_i=\frac{\Omega}{(2\pi)^3}\int\mathrm{d}\mathbf{k}\right)
\\
&=\frac{1}{(2\pi)^3}\int\mathrm{d}\mathbf{k} f\left(\beta(\varepsilon(\mathbf{k})-\mu)\right) e^{i\mathbf{k} \cdot (\mathbf{r} - \mathbf{r}')} \\
\end{align*}
ここで、$|\mathbf{r} - \mathbf{r}'| = r$ とすると、$\rho(r)$は球対称で、
$$
\rho(r) = \frac{1}{(2\pi)^3} \int d^3k f\left( \beta\left(\frac{k^2}{2}-\mu\right) \right) e^{i\mathbf{k} \cdot (\mathbf{r} - \mathbf{r}')}
$$
球座標系で $d^3k = k^2 \sin\theta dk d\theta d\phi$、かつ$\mathbf{k} \cdot (\mathbf{r} - \mathbf{r}') = kr \cos\theta$ とすれば：
$$
\rho(R) = \frac{1}{(2\pi)^3} \int_0^\infty dk k^2 f\left( \beta\left(\frac{k^2}{2}-\mu\right) \right) \int_0^\pi d\theta \sin\theta e^{ikr \cos\theta} \int_0^{2\pi} d\phi
$$
ここで、
\begin{align*}
\int_0^{2\pi} d\phi &= 2\pi \\
\int_0^\pi \sin\theta\ e^{ikr \cos\theta} d\theta
&=-\int_1^{-1} e^{ikrx} dx \qquad(x=\cos\theta) \\
&= \left[\frac{e^{ikrx}}{ikr}
\right]_{-1}^1 \\
&=\frac{2\sin(kr)}{kr}
\end{align*}
よって、
$$
\rho(r) = \frac{1}{2\pi^2 r} \int_0^\infty dk k f\left( \beta\left(\frac{k^2}{2}-\mu\right) \right) \sin(kr)
$$
この形式は扱いにくい。具体的には、
- $f(\varepsilon)$は1から0へ緩やかに変化する関数であるのに対して、$f'(\varepsilon)$は低温でフェルミエネルギー近傍でのみ非ゼロになり、積分の有効範囲が狭まること
- $\sin(kr)$は奇関数なのに対し、$\cos(kr)$は偶関数なので積分範囲をフーリエ変換と同様の$-\infty\sim\infty$に拡張できること
から、部分積分を行う。
\begin{align*}
\int_0^\infty dk kf\left(\beta\left(\frac{k^2}{2}-\mu\right)\right) \sin(kr)
&=\left[k f\left(\cdots\right) \frac{-\cos(kr)}{r}\right]_0^\infty \\
&\qquad -\int_0^\infty dk \left\{f\left(\cdots\right)+\beta k^2f'\left(\cdots\right)\right\} \frac{-\cos(kr)}{r} \\
&=\frac{1}{r}\int_0^\infty dk \left\{f\left(\cdots\right)
+\beta k^2f'\left(\cdots\right)\right\}\cos(kr)
\end{align*}
この第1項についても部分積分を行うと、
\begin{align*}
\int_0^\infty dk f\left(\beta\left(\frac{k^2}{2}-\mu\right)\right)\cos(kr)
&=\left[f\left(\cdots\right) \frac{\sin(kr)}{r}\right]_0^\infty \\
&\qquad -\frac{1}{r}\int_0^\infty dk\beta k f' \left(\cdots\right)\sin(kr)\\
&=-\frac{1}{r}\int_0^\infty dk\beta k f' \left(\cdots\right)\sin(kr)
\end{align*}
ここで、
$$
f\left(\beta\left(\frac{k^2}{2}-\mu\right)\right)=\frac{1}{e^{\beta\left(\frac{k^2}{2}-\mu\right)}+1}
$$
より、
\begin{align*}
\frac{\partial}{\partial k}
f\left(\beta\left(\frac{k^2}{2}-\mu\right)\right)
&=-\frac{1}{\left(e^{\beta\left(\frac{k^2}{2}-\mu\right)}+1\right)^2}
\cdot e^{\beta\left(\frac{k^2}{2}-\mu\right)}\cdot\beta k \\
&=f'\left(\beta\left(\frac{k^2}{2}-\mu\right)\right)\cdot\beta k
\end{align*}
を用いた。
よって、
\begin{align*}
\rho(r)
&=\frac{1}{2\pi^2 r}
\left(
\frac{1}{r}\int_0^\infty dk \beta k^2f'\left(\cdots\right)\cos(kr)
-\frac{1}{r^2}\int_0^\infty dk\beta k f' \left(\cdots\right)\sin(kr)
\right) \\
&=\frac{\beta}{2\pi^2}\frac{1}{r}\int_0^\infty dk f'\left(\cdots\right)
\left(\frac{k^2\cos(kr)}{r}-\frac{k\sin(kr)}{r^2}\right) \\
&=\frac{\beta}{2\pi^2}\frac{1}{r}\frac{\partial}{\partial r}\int_0^\infty dk f'\left(\cdots\right) \left(\frac{k\sin(kr)}{r}\right) \\
&=-\frac{\beta}{2\pi^2}\frac{1}{r}\frac{\partial}{\partial r}\frac{1}{r}\frac{\partial}{\partial r}\int_0^\infty dk f'\left(\cdots\right)\cos(kr) \\
&=-\frac{\beta}{(2\pi)^2}\frac{1}{r}\frac{\partial}{\partial r}\frac{1}{r}\frac{\partial}{\partial r}\int_{-\infty}^\infty dk f'\left(\cdots\right)\cos(kr) \\
\end{align*}
（原著論文(https://arxiv.org/pdf/cond-mat/9804013) では先頭に"-"がついていおりテキスト誤植？）
ここで、$T=0$の場合を考える。$f'\left(\beta\left(\frac{k^2}{2}-\mu\right)\right)
=-\delta\left(\beta\left(\frac{k^2}{2}-\mu\right)\right)$で、任意の関数$f(k)$に対して、
\begin{align*}
&\beta\int_{-\infty}^\infty dk\delta\left(\beta\left(\frac{k^2}{2}-\frac{k_F ^2}{2}\right)\right)f(k) \\
&=\int_{-\frac{\beta}{2}k_F^2}^\infty\frac{dx}{\sqrt{\frac{2x}{\beta}+k_F^2}}\delta(x)f\left(\sqrt{\frac{2x}{\beta}+k_F^2}\right) \\
&\quad +\int_\infty^{-\frac{\beta}{2}k_F^2}\frac{dx}{-\sqrt{\frac{2x}{\beta}+k_F^2}}\delta(x)f\left(-\sqrt{\frac{2x}{\beta}+k_F^2}\right) \\
&\qquad\left(\text{ただし}\; x=\beta\frac{k^2-k_F^2}{2}\text{より}\beta dk =\frac{dx}{k}=\frac{dx}{\pm\sqrt{\frac{2x}{\beta}+k_F^2}}\right)\\
&=\frac{f(k_F)}{k_F}+\frac{f(-k_F)}{k_F}
\end{align*}
となるため、
\begin{align*}
\rho(r)
&=-\frac{\beta}{(2\pi)^2}\frac{1}{r}\frac{\partial}{\partial r}\frac{1}{r}\frac{\partial}{\partial r}\int_{-\infty}^\infty dk f'\left(\cdots\right)\cos(kr) \\
&=\frac{1}{(2\pi)^2}\frac{1}{r}\frac{\partial}{\partial r}\frac{1}{r}\frac{\partial}{\partial r}\int_{-\infty}^\infty dk \beta\delta\left(\cdots\right)\cos(kr) \\
&=\frac{1}{2\pi^2}\frac{1}{r}\frac{\partial}{\partial r}\frac{1}{r}\frac{\partial}{\partial r}
\frac{\cos(k_Fr)}{k_F} \\
&=\frac{1}{2\pi^2}\frac{1}{r}\frac{\partial}{\partial r}\frac{-\sin(k_Fr)}{r} \\
&=\frac{1}{2\pi^2}\frac{1}{r}\left\{
\frac{-1}{r^2}\left(-\sin(k_Fr)\right)
+\frac{k_F}{r}\left(-\cos(k_Fr)\right)
\right\} \\
&=\frac{k_F^3}{6\pi^2}
\left[
3
\frac{\sin(y)-y\cos(y)}{y^3}
\right]
\end{align*}
ただし、$y=k_Fr$とおいた。
Hartree-Fock Approximation
Hartree-Fock近似において、1電子軌道は次の方程式の解。（ただし、 $\psi_i:=\psi_i^{\sigma_i}(\mathbf{r}), \quad \hat{h}:=-\frac{1}{2} \nabla^2 + V_{\text{ext}}(\mathbf{r}) , \quad \hat{g}:=\frac{1}{|\mathbf{r}-\mathbf{r}'|}$ とする）
\left[
\hat{h}+\sum_{j,\sigma_j}\int d\mathbf{r}'\psi_j^*\psi_j\hat{g}\right]\psi_i
- \sum_j\int d\mathbf{r}'\psi_j^*\psi_i\hat{g}\psi_j = \varepsilon\psi_i
ここで、平面波がこの方程式の解になることを示せる。実際、
$$
\psi_{k}(\mathbf{r}) = \left( \frac{1}{\Omega^{1/2}} \right) e^{i \mathbf{k} \cdot \mathbf{r}}
$$
とすると、均一ガスの仮定から$\rho(\mathbf{r})=n$、つまり$V_{\text{ext}}(\mathbf{r})=-\int d\mathbf{r'}\frac{n}
{\left| \mathbf{r}-\mathbf{r}'\right|}$となるので、
\begin{align*}
\hat{h}\psi_k(\mathbf{r})
&=\frac{k^2}{2}\psi_k(\mathbf{r})
-\int d\mathbf{r'}\frac{n}
{\left| \mathbf{r}-\mathbf{r}'\right|}\psi_k(\mathbf{r})\\
\left[\sum_{k',\sigma_k'}
d\mathbf{r}'\psi_{k'}^*\psi_{k'}\hat{g}\right]\psi_k
&=\int d\mathbf{r'}\frac{n}
{\left| \mathbf{r}-\mathbf{r}'\right|}\psi_k(\mathbf{r})\\
\sum_{k'}\int d\mathbf{r}'\psi_{k'}^*\psi_k\hat{g}\psi_{k'}
&=\frac{1}{\Omega}\sum_{k'}\int d\mathbf{r}'
\frac{e^{i(\mathbf{k}-\mathbf{k}')\cdot\mathbf{r}'}\delta_{\sigma\sigma'}}{\left| \mathbf{r}-\mathbf{r}'\right|}
\frac{e^{i\mathbf{k}'\cdot\mathbf{r}}}{\sqrt{\Omega}} \\
&=\frac{1}{\Omega}\sum_{k'}\frac{4\pi}{|\mathbf{k}-\mathbf{k}'|^2}
\cdot e^{i(\mathbf{k}-\mathbf{k}')\cdot\mathbf{r}}
\frac{e^{i\mathbf{k}'\cdot\mathbf{r}}}{\sqrt{\Omega}} \\
&=\frac{1}{\Omega}\sum_{k'}\frac{4\pi}{|\mathbf{k}-\mathbf{k}'|^2}
\psi_k(\mathbf{r})
\\
\end{align*}
ここで、クーロンポテンシャルのフーリエ変換の式
\begin{align*}
\int d\mathbf{r}\frac{e^{-i\mathbf{k}\cdot\mathbf{r}}}{|\mathbf{r}|}=\frac{4\pi}{|\mathbf{k}|^2}
\end{align*}
より、
\begin{align*}
\int d\mathbf{r}'
\frac{e^{i(\mathbf{k}-\mathbf{k}')\cdot(\mathbf{r}'-\mathbf{r})}}
{|\mathbf{r}-\mathbf{r}'|}
=\frac{4\pi}{|\mathbf{k}-\mathbf{k}'|^2}
\end{align*}
が成り立つことを用いた。
よって、以下のように平面波はHFeq.の解になる。
\left[
\frac{k^2}{2}
-\frac{1}{\Omega}\sum_{k'}\frac{4\pi}{|\mathbf{k}-\mathbf{k}'|^2}
\right]\psi_k
=\varepsilon_k\psi_k
$k'$に関する和を積分で表記すると、
\left[
\frac{k^2}{2}
-\frac{1}{(2\pi)^3}\int d\mathbf{k}'\frac{4\pi}{|\mathbf{k}-\mathbf{k}'|^2}
\right]\psi_k
=\varepsilon_k\psi_k
\quad\left(\because \sum_k=\frac{\Omega}{(2\pi)^3}\int\mathrm{d}\mathbf{k}\right)
ここで、$x=k/k_F$として積分を実行すると以下が成り立つ。
4\pi\int_{k\le k_F} \frac{d\mathbf{k}'}{(2\pi)^3}\frac{1}{|\mathbf{k}-\mathbf{k}'|^2}
=\frac{k_F}{\pi} \left[ 1 + \frac{1 - x^2}{2x} \ln \left| \frac{1 + x}{1 - x}\right|\right]
よって、
\varepsilon_k=\frac{k^2}{2}-\frac{k_F}{\pi} \left[ 1 + \frac{1 - x^2}{2x} \ln \left| \frac{1 + x}{1 - x}\right|\right]
となる。
上記の積分の証明
最初の式では、次の積分を評価する必要があります：
$$\int_{k < k_F} \frac{d^3 k'}{(2\pi)^3} \frac{1}{|k - k'|^2}$$
球座標系で計算しましょう。まず、$\vec{k}$を$z$軸の方向に選ぶと、$|\vec{k} - \vec{k}'|^2 = k^2 + k'^2 - 2kk'\cos\theta$となります。ここで$\theta$は$\vec{k}$と$\vec{k}'$の間の角度です。
この積分は次のように書けます：
$$\int_{k < k_F} \frac{d^3 k'}{(2\pi)^3} \frac{1}{|k - k'|^2} = \int_0^{k_F} \frac{k'^2 dk'}{(2\pi)^3} \int_0^{\pi} \frac{2\pi \sin\theta d\theta}{k^2 + k'^2 - 2kk'\cos\theta}$$
角度積分から始めましょう：
$$\int_0^{\pi} \frac{\sin\theta d\theta}{k^2 + k'^2 - 2kk'\cos\theta}$$
変数置換 $u = \cos\theta$ とすると、$du = -\sin\theta d\theta$ なので：
$$\int_0^{\pi} \frac{\sin\theta d\theta}{k^2 + k'^2 - 2kk'\cos\theta} = -\int_1^{-1} \frac{du}{k^2 + k'^2 - 2kk'u} = \int_{-1}^{1} \frac{du}{k^2 + k'^2 - 2kk'u}$$
この積分を計算します：
$$\int_{-1}^{1} \frac{du}{k^2 + k'^2 - 2kk'u} = \frac{1}{2kk'}\ln\left|\frac{k^2 + k'^2 + 2kk'}{k^2 + k'^2 - 2kk'}\right| = \frac{1}{2kk'}\ln\left|\frac{(k+k')^2}{(k-k')^2}\right| = \frac{1}{kk'}\ln\left|\frac{k+k'}{k-k'}\right|$$
したがって、元の積分は次のようになります：
$$\int_{k < k_F} \frac{d^3 k'}{(2\pi)^3} \frac{1}{|k - k'|^2} = \frac{1}{4\pi^2} \int_0^{k_F} \frac{k'^2}{kk'}\ln\left|\frac{k+k'}{k-k'}\right| dk' = \frac{1}{4\pi^2 k} \int_0^{k_F} k' \ln\left|\frac{k+k'}{k-k'}\right| dk'$$
ここで変数置換 $y = \frac{k'}{k}$ とすると、$dk' = k dy$ および積分範囲は $0 \leq y \leq \frac{k_F}{k} = \frac{1}{x}$ となります（ここで $x = \frac{k}{k_F}$）：
$$\frac{1}{4\pi^2 k} \int_0^{k_F} k' \ln\left|\frac{k+k'}{k-k'}\right| dk' = \frac{1}{4\pi^2 k} \int_0^{1/x} k^2 y \ln\left|\frac{1+y}{1-y}\right| dy = \frac{k}{4\pi^2} \int_0^{1/x} y \ln\left|\frac{1+y}{1-y}\right| dy$$
この積分を部分積分法で解きます：
$$\int y \ln\left|\frac{1+y}{1-y}\right| dy = \frac{y^2}{2} \ln\left|\frac{1+y}{1-y}\right| - \int \frac{y^2}{2} \cdot \frac{2}{1-y^2} dy$$
$$= \frac{y^2}{2} \ln\left|\frac{1+y}{1-y}\right| - \int \frac{y^2}{1-y^2} dy$$
次に、
$$\int \frac{y^2}{1-y^2} dy = \int \frac{1+y^2-1}{1-y^2} dy = \int \left(-1 + \frac{1}{1-y^2}\right) dy = -y + \int \frac{1}{1-y^2} dy$$
$$= -y + \frac{1}{2}\ln\left|\frac{1+y}{1-y}\right|$$
元の積分に戻ると：
$$\int y \ln\left|\frac{1+y}{1-y}\right| dy = \frac{y^2}{2} \ln\left|\frac{1+y}{1-y}\right| - \left(-y + \frac{1}{2}\ln\left|\frac{1+y}{1-y}\right|\right)$$
$$= \frac{y^2}{2} \ln\left|\frac{1+y}{1-y}\right| + y - \frac{1}{2}\ln\left|\frac{1+y}{1-y}\right|$$
$$= y + \frac{y^2-1}{2}\ln\left|\frac{1+y}{1-y}\right|$$
これを上限 $\frac{1}{x}$ と下限 $0$ で評価します：
$$\left[y + \frac{y^2-1}{2}\ln\left|\frac{1+y}{1-y}\right|\right]_0^{1/x}$$
下限 $y=0$ での値：
- $y = 0$
- $\ln\left|\frac{1+y}{1-y}\right| = \ln(1) = 0$
よって $0$
上限 $y=\frac{1}{x}$ での値：
$$\frac{1}{x} + \frac{(\frac{1}{x})^2-1}{2}\ln\left|\frac{1+\frac{1}{x}}{1-\frac{1}{x}}\right| = \frac{1}{x} + \frac{\frac{1}{x^2}-1}{2}\ln\left|\frac{x+1}{x-1}\right|$$
したがって、元の積分は：
$$\frac{k}{4\pi^2} \left[\frac{1}{x} + \frac{\frac{1}{x^2}-1}{2}\ln\left|\frac{x+1}{x-1}\right|\right]$$
係数の整理：
$$\frac{k}{4\pi^2} \cdot \frac{1}{x} = \frac{k_F}{4\pi^2}$$
$$\frac{k}{4\pi^2} \cdot \frac{\frac{1}{x^2}-1}{2}\ln\left|\frac{x+1}{x-1}\right| = \frac{k_F}{4\pi^2} \cdot \frac{1-x^2}{2x}\ln\left|\frac{1+x}{1-x}\right|$$
合わせると：
$$\frac{k_F}{4\pi^2}\left[1 + \frac{1-x^2}{2x}\ln\left|\frac{1+x}{1-x}\right|\right]$$
よって、以下が成り立つことが証明できました。
4\pi\int_{k\le k_F} \frac{d\mathbf{k}'}{(2\pi)^3}\frac{1}{|\mathbf{k}-\mathbf{k}'|^2}
=\frac{k_F}{\pi} \left[ 1 + \frac{1 - x^2}{2x} \ln \left| \frac{1 + x}{1 - x}\right|\right]
ここで、
f(x) = -\left(1 + \frac{1 - x^2}{2x} \ln \left| \frac{1 + x}{1 - x}\right|\right)
と置くと、
\varepsilon_k=\frac{k^2}{2}+\frac{k_F}{\pi}f(x)
となるが、
- $f(x)$は任意の$x$に対して負である。つまり、HF近似によって固有値は減少する。
- バンドの底では$f(0)=-2$となる。
- フェルミ面近傍$(x\to 1)$で$f(x)$の勾配は対数発散する。一方で$f(x)\to -1$となる。
- よって、バンド幅は$W=\varepsilon_{k_F}-\varepsilon_0=\left(\frac{k_F^2}{2}-\frac{k_F}{\pi}\right)-\left(-\frac{2k_F}{\pi}\right)=\frac{k_F^2}{2}+\frac{k_F}{\pi}$つまり、$\Delta W=\frac{k_F}{\pi}$だけ交換相互作用がない系に対して増加する。
HF固有値をフェルミエネルギー$\left(\frac{k_F^2}{2}-\frac{k_F}{\pi}\right)$基準にすると、
\varepsilon_k=\frac{k_F^2}{2}\left\{(x^2-1)+\frac{2}{\pi k_F}[f(x)+1]\right\}
フェルミ面において、速度$d\varepsilon/dk$は発散するが、これは実験と明らかに矛盾することが示されている。よって、これはHFに付随する失敗であり、任意の金属に当てはまる。
The Exchange Energy and Exchange Hole
均一ガスモデルの、交換エネルギー及び交換ホールを求める。
電子あたりの全エネルギーは、交換エネルギーの和をとり、ダブルカウントを防ぐために1/2をかけることに注意すると、$f(x)$の平均が$-\frac{3}{2}$であることから、
\begin{align*}
\epsilon_x^\sigma
&=E_x^\sigma/N^\sigma \\
&=\frac{1}{2}\frac{k_F^\sigma}{\pi}\sum_x \frac{f(x)}{N^\sigma} \\
&=-\frac{3}{4\pi}k_F^\sigma \\
&=-\frac{3}{4}\left(\frac{6}{\pi}n^\sigma\right)^\frac{1}{3}
\end{align*}
平均の意味
\begin{align*}
\sum_x \frac{f(x)}{N^\sigma}
&=\frac{\int_0^{k_F} dk 4\pi k^2 f(k/k_F)}{\frac{4\pi k_F^3}{3}} \\
&=3\int_0^1 dx x^2 f(x) \\
&=-\frac{3}{2}
\end{align*}
ここで、$\left(\frac{6}{\pi}n^\sigma\right)^\frac{1}{3}=\frac{1}{\pi}\left(\frac{9\pi}{4}\right)^\frac{1}{3}/r_s$より、偏極していない場合は、$\epsilon_x \equiv \epsilon_x^\uparrow = \epsilon_x^\downarrow =-\frac{3}{4\pi}\left(\frac{9\pi}{4}\right)^\frac{1}{3}/r_s$
部分的に偏極している場合は、交換エネルギーは、全密度$n=n^\uparrow+n^\downarrow$と分極率$\zeta=\frac{n^\uparrow-n^\downarrow}{n}$を用いて、
\epsilon_x(n,\zeta)=\epsilon_x(n,0)
+\left[\epsilon_x(n,1)-\epsilon_x(n,0)\right] f_x(\zeta)
ただし、
f_x(\zeta)=\frac{1}{2}\frac{(1+\zeta)^{4/3}+(1-\zeta)^{4/3}-2}{2^{1/3}-1}
となる。
ここで、交換ホール$g_x$は
g_x(\mathbf{r}_1,\sigma_1;\mathbf{r}_2,\sigma_2)
=\delta_{\sigma_1,\sigma_2} g_x^{\sigma_1,\sigma_1}(|\mathbf{r}_1-\mathbf{r}_2|)
となり、$y=k_F^\sigma r$とすると、
\begin{align*}
g_x^{\sigma,\sigma}(|\mathbf{r}_1-\mathbf{r}_2|)
&=1-\frac{|\rho_\sigma(|\mathbf{r}_1-\mathbf{r}_2|)|^2}{n^\sigma(\mathbf{r}_1)n^\sigma(\mathbf{r}_2)} \\
&=1
-\frac{\left|
\frac{(k_F^\sigma)^3}{6\pi^2}\left[3\frac{\sin(y)-y\cos(y)}{y^3}\right]
\right|^2}
{\frac{(k_F^\sigma)^3}{6\pi^2}\cdot\frac{(k_F^\sigma)^3}{6\pi^2}} \\
&=1-\left[3\frac{\sin(y)-y\cos(y)}{y^3}\right]^2
\end{align*}
となる。均一ガスにおいて、フェルミオンの交換ホールは常に負になることが分かる。
5.3 Correlation Hole and Energy
多体系において、粒子には集団的な相関が生じ、粒子間の相互作用を低減する"遮蔽"が生じる。これにより系全体のエネルギーは低下する。
Thomas-Fermi Screening
フーリエ成分$k$の静的な外部電荷密度に対する電子ガスの応答を考え遮蔽を求める。波数$k$の応答は電子のエネルギー変化で決定されるが、Thomas-Fermi近似では、これが電子密度のみの関数であるとする。
まず、電荷密度$\rho^\text{ext}(\mathbf{r})$の正電荷のみにより生じる電位を$\phi^\text{ext}(\mathbf{r})$とする。$E=-\nabla \phi$と、Maxwell方程式より$\nabla E=\frac{\rho^\text{ext}}{\varepsilon}$から、次のポアソン方程式が成り立つ。（参照： https://ja.wikipedia.org/wiki/電磁気量の単位系 において、$\lambda=4\pi, \varepsilon_0=1$）
-\nabla^2 \phi^\text{ext}(\mathbf{r})=4\pi\rho^\text{ext}(\mathbf{r})
次に、正電荷と遮蔽電子が合わさって生じる電荷密度を$\rho(\mathbf{r})$、電位を$\phi(\mathbf{r})$とする。同様に以下が成り立つ。
-\nabla^2 \phi(\mathbf{r})=4\pi\rho(\mathbf{r})
ここで、
\rho^\text{ind}(\mathbf{r})=\rho(\mathbf{r})-\rho^\text{ext}(\mathbf{r})
は正電荷によってもたらされた電子ガスの密度変化である。
誘電体の理論との類似性、及び対称性から、$\phi$と$\phi^\text{ext}$に次の線形な関係式を仮定する。
\phi^\text{ext}(\mathbf{r})
=\int d\mathbf{r}' \epsilon(\mathbf{r}-\mathbf{r}')\phi(\mathbf{r}')
これをフーリエ変換して、
\phi^\text{ext}(\mathbf{k})=\epsilon(\mathbf{k})\phi(\mathbf{k})
が成り立つ。ここで、ポアソン方程式のフーリエ変換から、
\begin{align*}
k^2 \phi(\mathbf{k})&=4\pi\rho(\mathbf{k}) \\
k^2 \phi^\text{ext}(\mathbf{k})&=4\pi\rho^\text{ext}(\mathbf{k})
\end{align*}
となるが、$\rho^\text{ind}$と$\phi$に線形な関係、つまりフーリエ変換後に
\rho^\text{ind}(\mathbf{k})=\chi(\mathbf{k})\phi(\mathbf{k})
が成り立つと仮定すると、
\frac{k^2}{4\pi}(\phi(\mathbf{k})-\phi^\text{ext}(\mathbf{k}))
=\rho^\text{ind}(\mathbf{k})=\chi(\mathbf{k})\phi(\mathbf{k})
よって、
\phi(\mathbf{k})
=\phi^\text{ext}(\mathbf{k})/\left(1-\frac{4\pi}{k^2}\chi(\mathbf{k})\right)
つまり、
\epsilon(\mathbf{k})
=1-\frac{4\pi}{k^2}\chi(\mathbf{k})
=1-\frac{4\pi}{k^2}\frac{\rho^\text{ind}(\mathbf{k})}{\phi(\mathbf{k})}
ここで、$\chi(\mathbf{k})=-e^2\frac{\partial n_0}{\partial\mu}$より、
証明
n_0(\mu)
=\int \frac{d\mathbf{k}}{4\pi^3}\frac{1}{\exp[\beta(\hbar^2 k^2/2m-\mu)]+1}
とすると、
\rho^\text{ind}(\mathbf{r})
=-e[n_0(\mu+e\phi(\mathbf{r})-n_0(\mu)]
=-e^2\frac{\partial n_0}{\partial\mu}\phi(\mathbf{r})
より、$\chi(\mathbf{k})=-e^2\frac{\partial n_0}{\partial\mu}$が成立。
$k_\text{TF}^2=4\pi e^2\frac{\partial n_0}{\partial\mu}$とすると、
\epsilon(\mathbf{k})=1+\frac{k_\text{TF}^2}{k^2}
$k_\text{TF}^2$の影響を見るために、点電荷の場合を考える。つまり、
\phi^\text{ext}(\mathbf{r})=\frac{Q}{r},\quad
\phi^\text{ext}(\mathbf{k})=\frac{4\pi Q}{k^2}
遮蔽後の電位は、
\phi(\mathbf{k})
=\frac{1}{\epsilon(\mathbf{k})}\phi^\text{ext}(\mathbf{k})
=\frac{4\pi Q}{k^2+k_\text{TF}^2}
よって以下の変化を与える。
\frac{1}{k^2}\rightarrow \frac{1}{k^2+k_\text{TF}^2}
また、
\phi(\mathbf{r})
=\int \frac{d\mathbf{k}}{(2\pi)^3}e^{i\mathbf{k}\cdot\mathbf{r}}
\frac{4\pi Q}{k^2+k_\text{TF}^2}=\frac{Q}{r}e^{-\mathbf{k}_\text{TF}\cdot\mathbf{r}}
となり、$1/k_\text{TF}$のオーダー以上で無視できるポテンシャルであることが分かる。これは遮蔽ポテンシャル或いはYukawaポテンシャルと呼ばれる。
ここで、十分低温で$\frac{\partial n_0}{\partial\mu}=D(\varepsilon_F)=\frac{mk_F}{\hbar^2 \pi^2}$より、
k_\text{TF}
=\sqrt{4\pi e^2\frac{mk_F}{\hbar^2 \pi^2}}
=\sqrt{\frac{4k_F}{\pi a_0}}
=\left(\frac{12}{\pi}\right)^{1/3}r_s^{-1/2}\quad(\because a_0=1)
Correlation Energy
相関エネルギー及び相関ホールは解析的に求まらない。Wignarによって当初提案されたものは、相関エネルギーは
- 低密度極限では、電子は"Wigner結晶"を構築して、BCC格子上の点電荷による静電エネルギーになる
- 高密度極限では、定数に漸近する
との仮説に基づくもので、
\epsilon_c=-\frac{0.44}{r_s+7.8}
であった。これらは後に修正され、GellmannとBreucknerにより、高密度極限では、偏極がない場合に
\epsilon_c(r_s)=0.311 \ln(r_s) - 0.048 + r_s(A\ln(r_s)+C)+\cdots,
となり、低密度極限ではWigner結晶を構築してゼロ点振動により、
\epsilon_c(r_s)=\frac{a_1}{r_s}+\frac{a_2}{r_s^{3/2}}+\frac{a_3}{r_s^2}+\cdots ,
となるとされている。基底状態に関する結果として最も正確なのはquantum Monte Carlo (QMC)を用いたものである。
スピン偏極がある場合、PZによる単純な形式では相関エネルギーも交換エネルギーと同様に変化すると考える。
\epsilon_c(n,\zeta)=\epsilon_c(n,0)
+\left[\epsilon_c(n,1)-\epsilon_c(n,0)\right] f_x(\zeta)
少し複雑なVWNによる形式の方が、少しQMCの結果に近づくことが知られている。
ここで、交換相関"ホール"からどのように交換相関エネルギーが生まれるか見ておこう。それぞれの電子とホールとの相互作用によるポテンシャルエネルギーは次のように書ける。
\epsilon_\text{xc}^\text{pot}(r_s)
=E_\text{xc}^\text{pot}/N
=\frac{1}{N}\left[\left<\hat{V}_\text{int}\right>-E_\text{Hartree}(n)\right]
=\frac{1}{2n}e^2\int d^3 r\frac{n_\text{xc}(|\mathbf{r}|)}{|\mathbf{r}|}
ここで、$n_\text{xc}(|\mathbf{r}|)$は球対称で、電子密度の関数である。交換相互作用による反発の効果、相関相互作用による遮蔽の効果はいずれもエネルギーを低下させるので$\epsilon_\text{xc}^\text{pot}$は負である。しかしながら、電子相関によってエネルギーが低下すると運動エネルギーが増加するので、これが全て交換相関エネルギーではない。
5.4 Binding in sp- Bonded Metals
電子を均一ガスとして扱い、均一な背景としてのイオンのエネルギーを加えて電子あたりの全エネルギーを考えると、
\begin{align*}
\frac{E_\text{total}}{N}
&=t_0+\epsilon_x+\epsilon_c-\frac{1}{2}\frac{\alpha}{r_s}+\epsilon_R \\
&=\frac{1.105}{r_s^2}-\frac{0.458}{r_s}+\epsilon_c-\frac{1}{2}\frac{\alpha}{r_s}+\epsilon_R
\end{align*}
ただし、$\alpha$はMardelung定数、$\epsilon_R\left(=\frac{3}{2}\frac{R_c^2}{r_s^3}\right)$はコア同士の反発によるエネルギー。$\epsilon_c$を無視し、最密充填の代表的な値$\alpha=1.80$を代入すると、
\frac{r_s}{a_0}=0.814+\sqrt{0.899+3.31\left(\frac{R_c}{a_0}\right)^2}
となり、反発項を無視すると$r_s=1.76$と小さすぎるが、$R_c\approx 2a_0$とすると、$r_s\approx 4a_0$となり妥当な値を得る。この時、体積弾性率も数100GPaと妥当な値を得る。
5.5 Excitations and the Lindhard Dielectric Function
均一ガスの励起としては、電子数が変化（追加や除去）するものと、変化しないものがある。前者はバンドで評価できるが、Hartree-Fock法で求めたバンドは実験や高精度な計算とどの程度一致するだろうか？残念ながら、HF法によってバンド構造は$\frac{k^2}{2}$より曲率が大きかったが、実験（Na結晶に関するもの）では逆の傾向になることが分かっている。ただし良い出発点にはなる。
後者はLindhard Dielectric Functionとして計算することができる。低振動数での$\text{Im}\epsilon$のピークはDrude吸収と呼ばれ、実際に金属に存在することが知られている。一方、計算では$\omega>q v_F+\varepsilon_k$の時$\text{Im}\epsilon=0$となるのに対し、実際にはバンド間吸収があることが知られている。
参考文献
- Solid State Physics (Neil Ashcroft, N. Mermin 著 2021)
- https://arxiv.org/abs/cond-mat/9804013
- https://phy.ntnu.edu.tw/~changmc/Teach/SS/SSG_note/mchap09.pdf
- https://arxiv.org/pdf/1601.03544