---
title: "Jellium"
date: 2021-11-23
categories:
- Physics
- Quantum mechanics
- Perturbation
layout: "concept"
---

**Jellium**, also called the **uniform** or **homogeneous electron gas**,
is a theoretical material where all electrons are free,
and the ions' positive charge is smeared into a uniform background "jelly".
This simple model lets us study electron interactions easily.


## Without interactions

Let us start by neglecting electron-electron interactions.
This is clearly a dubious assumption, but we will stick with it for now.
For an infinitely large sample of jellium,
the single-electron states are simply plane waves.
We consider an arbitrary cube of volume $V$,
and impose periodic boundary conditions on it,
such that the single-particle orbitals are (suppressing spin):

$$\begin{aligned}
    \Inprod{\vb{r}}{\psi_{\vb{k}}}
    = \psi_{\vb{k}}(\vb{r})
    = \frac{1}{\sqrt{V}} \exp(i \vb{k} \cdot \vb{r})
    \qquad \quad
    \vb{k} = \frac{2 \pi}{V^{1/3}} (n_x, n_y, n_z)
\end{aligned}$$

Where $n_x, n_y, n_z \in \mathbb{Z}$.
This is a discrete (but infinite) set of independent orbitals,
so it is natural to use the
[second quantization](/know/concept/second-quantization/)
to write the non-interacting Hamiltonian $\hat{H}_0$,
where $\hbar^2 |\vb{k}|^2 / (2 m)$ is the kinetic energy
of the orbital with wavevector $\vb{k}$, and $s$ is the spin:

$$\begin{aligned}
    \hat{H}_0
    = \sum_{s} \sum_{\vb{k}} \frac{\hbar^2 |\vb{k}|^2}{2 m} \hat{c}_{s,\vb{k}}^\dagger \hat{c}_{s,\vb{k}}
\end{aligned}$$

Assuming that the temperature $T = 0$,
the $N$-electron ground state of this Hamiltonian
is known as the **Fermi sea** or **Fermi sphere** $\Ket{\mathrm{FS}}$,
and is constructed by filling up the single-electron states
starting from the lowest energy:

$$\begin{aligned}
    \Ket{\mathrm{FS}}
    = \prod_{s} \prod_{j = 1}^{N/2} \hat{c}_{s,\vb{k}_j}^\dagger \Ket{0}
\end{aligned}$$

Because $T = 0$, all the electrons stay in their assigned state.
The energy and wavenumber $|\vb{k}|$ of the highest filled orbital
are called the **Fermi energy** $\epsilon_F$ and **Fermi wavenumber** $k_F$,
and obey the expected kinetic energy relation:

$$\begin{aligned}
    \boxed{
        \epsilon_F
        = \frac{\hbar^2}{2 m} k_F^2
    }
\end{aligned}$$

The Fermi sea can be visualized in $\vb{k}$-space as a sphere with radius $k_F$.
Because $\vb{k}$ is discrete, the sphere's surface is not smooth,
but in the limit $V \to \infty$ it becomes perfect.

Now, we would like a relation between the system's parameters,
e.g. $N$ and $V$, and the resulting values of $\epsilon_F$ or $k_F$.
The total population $N$ must be given by:

$$\begin{aligned}
    N
    = \sum_{s} \sum_{\vb{k}} \matrixel{\mathrm{FS}}{\hat{c}_{s,\vb{k}}^\dagger \hat{c}_{s,\vb{k}}}{\mathrm{FS}}
    = \sum_{s} \frac{V}{(2 \pi)^3} \int_{-\infty}^\infty \matrixel{\mathrm{FS}}{\hat{c}_{s,\vb{k}}^\dagger \hat{c}_{s,\vb{k}}}{\mathrm{FS}} \dd{\vb{k}}
\end{aligned}$$

Where we have turned the sum over $\vb{k}$ into an integral with a constant factor,
by using that each orbital exclusively occupies a volume $(2 \pi)^3 / V$ in $\vb{k}$-space.

At zero temperature, this inner product can only be $0$ or $1$,
depending on whether $\vb{k}$ is outside or inside the Fermi sphere.
We can therefore rewrite using a
[Heaviside step function](/know/concept/heaviside-step-function/):

$$\begin{aligned}
    N
    = \sum_{s} \frac{V}{(2 \pi)^3} \int_{-\infty}^\infty \Theta(k_F - |\vb{k}|) \dd{\vb{k}}
    = 2 \frac{V}{(2 \pi)^3} \int_{-\infty}^\infty \Theta(k_F - |\vb{k}|) \dd{\vb{k}}
\end{aligned}$$

Where we realized that spin does not matter,
and replaced the sum over $s$ by a factor $2$.
In order to evaluate this 3D integral,
we go to [spherical coordinates](/know/concept/spherical-coordinates/)
$(|\vb{k}|, \theta, \varphi)$:

$$\begin{aligned}
    N
    &= \frac{V}{4 \pi^3} \int_0^{2 \pi} \int_0^\pi \int_0^\infty \Theta(k_F - |\vb{k}|) |\vb{k}|^2 \sin(\theta) \dd{|\vb{k}|} \dd{\theta} \dd{\varphi}
    \\
    &= \frac{V}{4 \pi^3} 4 \pi \int_0^{k_F} |\vb{k}|^2 \dd{|\vb{k}|}
    = \frac{V}{\pi^2} \bigg[ \frac{|\vb{k}|^3}{3} \bigg]_0^{k_F}
    = \frac{V}{3 \pi^2} k_F^3
\end{aligned}$$

Using that the electron density $n = N/V$,
we thus arrive at the following relation:

$$\begin{aligned}
    \boxed{
        k_F^3
        = 3 \pi^2 n
    }
\end{aligned}$$

This result also justifies our assumption that $T = 0$:
we can accurately calculate the density $n$ for many conducting materials,
and this relation then gives $k_F$ and $\epsilon_F$.
It turns out that $\epsilon_F$ is usually very large
compared to the thermal energy $k_B T$ at reasonable temperatures,
so we can conclude that thermal fluctuations are negligible.

Now, $\epsilon_F$ is the highest single-electron energy,
but about the total $N$-particle energy $E^{(0)}$?

$$\begin{aligned}
    E^{(0)}
    = \matrixel{\mathrm{FS}}{\hat{H}_0}{\mathrm{FS}}
    = \sum_{s} \sum_{\vb{k}} \frac{\hbar^2 |\vb{k}|^2}{2 m} \matrixel{\mathrm{FS}}{\hat{c}_{s,\vb{k}}^\dagger \hat{c}_{s,\vb{k}}}{\mathrm{FS}}
\end{aligned}$$

Once again, we turn the sum over $\vb{k}$ into an integral,
and recognize the spin's irrelevance:

$$\begin{aligned}
    E^{(0)}
    &= \sum_{s} \frac{V}{(2 \pi)^3} \int_{-\infty}^\infty \frac{\hbar^2 |\vb{k}|^2}{2 m}
    \matrixel{\mathrm{FS}}{\hat{c}_{\vb{k}}^\dagger \hat{c}_{\vb{k}}}{\mathrm{FS}} \dd{\vb{k}}
    \\
    &= \frac{\hbar^2 V}{8 \pi^3 m} \int_{-\infty}^\infty |\vb{k}|^2 \: \Theta(k_F - |\vb{k}|) \dd{\vb{k}}
\end{aligned}$$

In spherical coordinates,
we evaluate the integral and find that $E^{(0)}$ is proportional to $k_F^5$:

$$\begin{aligned}
    E^{(0)}
    &= \frac{\hbar^2 V}{8 \pi^3 m} \int_0^{2 \pi}
    \int_0^\pi \int_0^\infty \Big( |\vb{k}|^2 \: \Theta(k_F - |\vb{k}|) \Big) |\vb{k}|^2 \sin(\theta) \dd{|\vb{k}|} \dd{\theta} \dd{\varphi}
    \\
    &= \frac{\hbar^2 V}{8 \pi^3 m} 4 \pi \int_0^{k_F} |\vb{k}|^4 \dd{|\vb{k}|}
    = \frac{\hbar^2 V}{2 \pi^2 m} \bigg[ \frac{|\vb{k}|^5}{5} \bigg]_0^{k_F}
    = \frac{\hbar^2 V}{10 \pi^2 m} k_F^5
\end{aligned}$$

In general, it is more useful to consider
the average kinetic energy per electron $E^{(0)} / N$,
which we find to be as follows, using that $k_F^3 = 3 \pi^2 n$:

$$\begin{aligned}
    \boxed{
        \frac{E^{(0)}}{N}
        = \frac{3 \hbar^2}{10 m} k_F^2
        = \frac{3}{5} \epsilon_F
    }
    \:\sim\: n^{2/3}
\end{aligned}$$

Traditionally, this is expressed using a dimensionless parameter $r_s$,
defined as the radius of a sphere containing a single electron,
measured in Bohr radii $a_0 \equiv 4 \pi \varepsilon_0 \hbar^2 / (e^2 m)$:

$$\begin{aligned}
        \frac{4 \pi}{3} (a_0 r_s)^3
        = \frac{1}{n}
        = \frac{3 \pi^2}{k_F^3}
        \quad \implies \quad
        r_s
        = \Big( \frac{3}{4 \pi a_0^3 n} \Big)^{1/3}
        = \Big( \frac{9 \pi}{4} \Big)^{1/3} \frac{1}{a_0 k_F}
\end{aligned}$$

Such that the ground state energy can be rewritten in Rydberg units of energy like so:

$$\begin{aligned}
    \frac{E^{(0)}}{N}
    = \frac{3 \hbar^2}{10 m} \frac{4 \pi \varepsilon_0 e^2}{4 \pi \varepsilon_0 e^2} \frac{a_0^2 k_F^2}{a_0^2}
    = \frac{3 e^2}{40 \pi \varepsilon_0} \Big( \frac{9 \pi}{4} \Big)^{2/3} \frac{1}{a_0 r_s^2}
    \approx \frac{2.21}{r_s^2} \; \mathrm{Ry}
\end{aligned}$$


## With interactions

To include Coulomb interactions, let us try
[time-independent pertubation theory](/know/concept/time-independent-perturbation-theory/).
Clearly, this will give better results when the interaction is relatively weak, if ever.

The Coulomb potential is proportional to the inverse distance,
and the average electron spacing is roughly $n^{-1/3}$,
so the interaction energy $E_\mathrm{int}$ should scale as $n^{1/3}$.
We already know that the kinetic energy $E_\mathrm{kin} = E^{(0)}$ scales as $n^{2/3}$,
meaning perturbation theory should be reasonable
if $1 \gg E_\mathrm{int} / E_\mathrm{kin} \sim n^{-1/3}$,
so in the limit of high density $n \to \infty$.

The two-body Coulomb interaction operator $\hat{W}$
is as follows in second-quantized form:

$$\begin{aligned}
    \hat{W}
    = \frac{1}{2 V} \sum_{s_1 s_2} \sum_{\vb{k}_1 \vb{k}_2} \sum_{\vb{q} \neq 0} \frac{e^2}{\varepsilon_0 |\vb{q}|^2}
    \hat{c}_{s_1, \vb{k}_1 + \vb{q}}^\dagger \hat{c}_{s_2, \vb{k}_2 - \vb{q}}^\dagger \hat{c}_{s_2, \vb{k}_2} \hat{c}_{s_1, \vb{k}_1}
\end{aligned}$$

The first-order correction $E^{(1)}$ to the ground state (i.e. Fermi sea) energy
is then given by:

$$\begin{aligned}
    E^{(1)}
    = \matrixel{\mathrm{FS}}{\hat{W}}{\mathrm{FS}}
    = \frac{e^2}{2 \varepsilon_0 V} \sum_{s_1 s_2} \sum_{\vb{k}_1 \vb{k}_2} \sum_{\vb{q} \neq 0} \frac{1}{|\vb{q}|^2}
    \matrixel{\mathrm{FS}}{
        \hat{c}_{s_1, \vb{k}_1 + \vb{q}}^\dagger \hat{c}_{s_2, \vb{k}_2 - \vb{q}}^\dagger \hat{c}_{s_2, \vb{k}_2} \hat{c}_{s_1, \vb{k}_1}
    }{\mathrm{FS}}
\end{aligned}$$

This inner product can only be nonzero
if the two creation operators $\hat{c}^\dagger$
are for the same orbitals as the two annihilation operators $\hat{c}$.
Since $\vb{q} \neq 0$, this means that $s_1 = s_2$,
and that momentum is conserved: $\vb{k}_2 = \vb{k}_1 \!+\! \vb{q}$.
And of course both $\vb{k}_1$ and $\vb{k}_1 \!+\! \vb{q}$
must be inside the Fermi sphere,
to avoid annihilating an empty orbital.
Let $s = s_1$ and $\vb{k} = \vb{k}_1$:

$$\begin{aligned}
    E^{(1)}
    &= \frac{e^2}{2 \varepsilon_0 V} \sum_{s} \sum_{\vb{k}} \sum_{\vb{q} \neq 0} \frac{1}{|\vb{q}|^2}
    \matrixel{\mathrm{FS}}{
        \hat{c}_{s, \vb{k} + \vb{q}}^\dagger \hat{c}_{s, \vb{k}}^\dagger \hat{c}_{s, \vb{k} + \vb{q}} \hat{c}_{s, \vb{k}}
    }{\mathrm{FS}}
    \\
    &= \frac{- e^2}{2 \varepsilon_0 V} \sum_{s} \sum_{\vb{k}} \sum_{\vb{q} \neq 0} \frac{1}{|\vb{q}|^2}
    \matrixel{\mathrm{FS}}{
        \big( \hat{c}_{s, \vb{k} + \vb{q}}^\dagger \hat{c}_{s, \vb{k} + \vb{q}}\big) \big(\hat{c}_{s, \vb{k}}^\dagger \hat{c}_{s, \vb{k}}\big)
    }{\mathrm{FS}}
    \\
    &= \frac{- e^2}{2 \varepsilon_0 V} \sum_{s} \sum_{\vb{k}} \sum_{\vb{q} \neq 0} \frac{1}{|\vb{q}|^2}
    \Theta(k_F - |\vb{k}|) \:\Theta(k_F - |\vb{k} \!+\! \vb{q}|)
\end{aligned}$$

Next, we convert the sum over $\vb{q}$ into an integral in spherical coordinates.
Clearly, $\vb{q}$ is the "jump" made by an electron from one orbital to another,
so the largest possible jump
goes from a point on the Fermi surface to the opposite point,
and thus has length $2 k_F$.
This yields the integration limit, and therefore leads to:

$$\begin{aligned}
    E^{(1)}
    &= \frac{- e^2}{(2 \pi)^3 \varepsilon_0} \sum_{\vb{k}}
    \int_0^{2 \pi} \!\!\int_0^\pi \!\!\int_0^\infty \Theta(k_F \!-\! |\vb{k}|) \: \Theta(k_F \!-\! |\vb{k} \!+\! \vb{q}|) \frac{|\vb{q}|^2}{|\vb{q}|^2}
    \sin(\theta_q) \dd{|\vb{q}|} \dd{\theta_q} \dd{\varphi_q}
    \\
    &= \frac{- e^2}{2 \pi^2 \varepsilon_0} \sum_{\vb{k}}
    \int_0^{2 k_F} \Theta(k_F \!-\! |\vb{k}|) \: \Theta(k_F \!-\! |\vb{k} \!+\! \vb{q}|) \dd{|\vb{q}|}
\end{aligned}$$

Where we have used that the direction of $\vb{q}$,
i.e. $(\theta_q,\varphi_q)$, is irrelevant,
as long as we define $\theta_k$ as
the angle between $\vb{q}$ and $\vb{k} \!+\! \vb{q}$
when we go to spherical coordinates $(|\vb{k}|, \theta_k, \varphi_k)$ for $\vb{k}$:

$$\begin{aligned}
    E^{(1)}
    &= \frac{- e^2 V}{16 \pi^5 \varepsilon_0} \int_0^{2 k_F} \!\!\!\!\int_0^{2 \pi} \!\!\!\int_0^\pi \!\!\!\int_0^\infty
    \!\Theta(k_F \!-\! |\vb{k}|) \: \Theta(k_F \!-\! |\vb{k} \!+\! \vb{q}|)
    \: |\vb{k}|^2 \sin(\theta_k) \dd{|\vb{k}|} \dd{\theta_k} \dd{\varphi_k} \dd{|\vb{q}|}
    \\
    &= \frac{- e^2 V}{16 \pi^5 \varepsilon_0} \int_0^{2 k_F} \!\!\!\!\int_0^{2 \pi} \!\!\!\int_0^\pi \!\!\!\int_0^{k_F}
    \!\Theta(k_F \!-\! |\vb{k} \!+\! \vb{q}|)
    \: |\vb{k}|^2 \sin(\theta_k) \dd{|\vb{k}|} \dd{\theta_k} \dd{\varphi_k} \dd{|\vb{q}|}
\end{aligned}$$

Unfortunately, this last step function is less easy to translate into integration limits.
In effect, we are trying to calculate the intersection volume of two spheres,
both with radius $k_F$, one centered on the origin (for $\vb{k}$),
and the other centered on $\vb{q}$ (for $\vb{k} \!+\! \vb{q}$).
Imagine a triangle with side lengths $|\vb{k}|$, $|\vb{q}|$ and $|\vb{k} \!+\! \vb{q}|^2$,
where $\theta_k$ is the angle between $|\vb{k}|$ and $|\vb{k} \!+\! \vb{q}|$.
The *law of cosines* then gives the following relation:

$$\begin{aligned}
    |\vb{k}|^2
    = |\vb{q}|^2 + |\vb{k} \!+\! \vb{q}|^2 - 2 |\vb{q}| |\vb{k} \!+\! \vb{q}| \cos(\theta_k)
\end{aligned}$$

We already know that $|\vb{k}| < k_F$ and $0 < |\vb{q}| < 2 k_F$,
so by isolating for $\cos(\theta_k)$,
we can obtain bounds on $\theta_k$ and $|\vb{k}|$.
Let $|\vb{k}| \to k_F$ in both cases, then:

$$\begin{aligned}
    \cos(\theta_k)
    = \frac{|\vb{k} \!+\! \vb{q}|^2 + |\vb{q}|^2 - |\vb{k}|^2}{2 |\vb{k} \!+\! \vb{q}| |\vb{q}|}
    &\:\:\underset{|\vb{q}| \to 0}{>}\:\:\: \frac{k_F^2 + |\vb{q}|^2 - k_F^2}{2 k_F |\vb{q}|}
    = \frac{|\vb{q}|}{2 k_F}
    \\
    &\underset{|\vb{q}| \to 2 k_F}{<}\:\: \frac{k_F^2 + 4 k_F ^2 - k_F^2}{2 k_F 2 k_F}
    = 1
\end{aligned}$$

Meaning that $0 < \theta_k < \arccos{|\vb{q}| / (2 k_F)}$.
To get a lower limit for $|\vb{k}|$, we "cheat" by artificially demanding
that $\vb{k}$ does not cross the halfway point between the spheres,
with the result that $|\vb{k}| \cos(\theta_k) > |\vb{q}|/2$.
Then, thanks to symmetry (both spheres have the same radius),
we just multiply the integral by $2$,
for $\vb{k}$ on the other side of the halfway point.

Armed with these integration limits, we return to calculating $E^{(1)}$,
substituting $\xi \equiv \cos(\theta_k)$:

$$\begin{aligned}
    E^{(1)}
    &= \frac{- e^2 V}{16 \pi^5 \varepsilon_0} 2 \int_0^{2 k_F} \!\!\!\int_0^{2 \pi} \!\!\int_0^{\arccos{|\vb{q}| / (2 k_F)}}
    \!\!\int_{|\vb{q}|/(2 \cos{\theta_k})}^{k_F} |\vb{k}|^2 \sin(\theta_k) \dd{|\vb{k}|} \dd{\theta_k} \dd{\varphi_k} \dd{|\vb{q}|}
    \\
    &= \frac{e^2 V}{8 \pi^5 \varepsilon_0} 2 \pi \int_0^{2 k_F} \!\!\!\int_1^{|\vb{q}| / (2 k_F)}
    \!\!\int_{|\vb{q}|/(2 \xi)}^{k_F} |\vb{k}|^2 \frac{\sin(\theta_k)}{\sin(\theta_k)} \dd{|\vb{k}|} \dd{\xi} \dd{|\vb{q}|}
    \\
    &= \frac{- e^2 V}{4 \pi^4 \varepsilon_0} \int_0^{2 k_F} \!\!\!\int_{|\vb{q}| / (2 k_F)}^1
    \!\!\int_{|\vb{q}|/(2 \xi)}^{k_F} |\vb{k}|^2 \dd{|\vb{k}|} \dd{\xi} \dd{|\vb{q}|}
\end{aligned}$$

Where we have used that $\varphi_k$ does not appear in the integrand.
Evaluating these integrals:

$$\begin{aligned}
    E^{(1)}
    &= \frac{- e^2 V}{4 \pi^4 \varepsilon_0} \int_0^{2 k_F} \!\!\!\int_{|\vb{q}| / (2 k_F)}^1
    \bigg[ \frac{|\vb{k}|^3}{3} \bigg]_{|\vb{q}|/(2 \xi)}^{k_F} \dd{\xi} \dd{|\vb{q}|}
    \\
    &= \frac{- e^2 V}{4 \pi^4 \varepsilon_0} \int_0^{2 k_F} \!\!\!\int_{|\vb{q}| / (2 k_F)}^1
    \bigg( \frac{k_F^3}{3} - \frac{|\vb{q}|^3}{24 \xi^3} \bigg) \dd{\xi} \dd{|\vb{q}|}
    \\
    &= \frac{- e^2 V}{4 \pi^4 \varepsilon_0} \int_0^{2 k_F}
    \bigg[ \frac{k_F^3}{3} x + \frac{|\vb{q}|^3}{48 \xi^2} \bigg]_{|\vb{q}| / (2 k_F)}^1 \dd{|\vb{q}|}
    \\
    &= \frac{- e^2 V}{4 \pi^4 \varepsilon_0} \int_0^{2 k_F}
    \bigg( \frac{k_F^3}{3} + \frac{|\vb{q}|^3}{48} - \frac{k_F^2 |\vb{q}|}{4} \bigg) \dd{|\vb{q}|}
    \\
    &= \frac{- e^2 V}{4 \pi^4 \varepsilon_0} \bigg[ \frac{k_F^3 |\vb{q}|}{3} + \frac{|\vb{q}|^4}{192} - \frac{k_F^2 |\vb{q}|^2}{8} \bigg]_0^{2 k_F}
    \\
    &= \frac{- e^2 V}{16 \pi^4 \varepsilon_0} k_F^4
    = \frac{- e^2 N}{16 \pi^4 \varepsilon_0 n} k_F^4
    = -\frac{3 e^2 N}{16 \pi^2 \varepsilon_0} k_F
\end{aligned}$$

Per particle, the first-order energy correction $E^{(1)}$
is therefore found to be as follows:

$$\begin{aligned}
    \boxed{
        \frac{E^{(1)}}{N}
        = -\frac{3 e^2}{16 \pi^2 \varepsilon_0} k_F
    }
\end{aligned}$$

This can also be written using the parameter $r_s$ introduced above, leading to:

$$\begin{aligned}
    \frac{E^{(1)}}{N}
    = -\frac{3 e^2}{16 \pi^2 \varepsilon_0} \frac{a_0 k_F}{a_0}
    = -\frac{3 e^2}{16 \pi^2 \varepsilon_0} \Big( \frac{9 \pi}{4} \Big)^{1/3} \frac{1}{a_0 r_s}
\end{aligned}$$

Consequently, for sufficiently high densities $n$,
the total energy $E$ per particle is given by:

$$\begin{aligned}
    \boxed{
        \frac{E}{N}
        \approx \bigg( \frac{2.21}{r_s^2} - \frac{0.92}{r_s} \bigg) \; \mathrm{Ry}
    }
\end{aligned}$$

Unfortunately, this is as far as we can go.
In theory, the second-order energy correction $E^{(2)}$ is as shown below,
but it turns out that it (and all higher orders) diverge:

$$\begin{aligned}
    E^{(2)}
    = \sum_{\Psi_n \neq \mathrm{FS}} \frac{\big| \matrixel{\mathrm{FS}}{\hat{W}}{\Psi_n} \big|^2}{E^{(0)} - E_n}
\end{aligned}$$

The only cure for this is to go to infinite order,
where all the infinities add up to a finite result.



## References
1.  H. Bruus, K. Flensberg,
    *Many-body quantum theory in condensed matter physics*,
    2016, Oxford.