---

Categories: Perturbation, Physics, Quantum mechanics.

Jellium

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 is a version of the Fermi gas model, which we extend by including electron-electron interactions using time-independent perturbation theory.

0th order

Let us start with the 0th order of the perturbation expansion. Without interactions or potentials, this is simply a Fermi gas, so the single-electron wavefunctions are just plane waves. For mathematical convenience, we consider these waves in a cube of volume $V$ with periodic boundaries, leading to a discrete spectrum of allowed wavevectors $\vb{k}$, which becomes continuous for $V \to \infty$.

The unperturbed many-particle Hamiltonian $\hat{H}_0$ is given below in the second quantization, 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}$$

Which, at absolute zero $T = 0$, has an $N$-electron ground state known as the Fermi sphere $\Ket{\mathrm{FS}}$ that 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}$$

From our analysis of the Fermi gas, we have an important result for the wavenumber $k_F = |\vb{k}_{N/2}|$ of the highest filled orbital, as a function of the particle density $n = N / V$:

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

Now, let us calculate the total $N$-particle ground state energy $E^{(0)}$ of the unperturbed system:

$$\begin{aligned} E^{(0)} = \matrixel{\mathrm{FS}}{\hat{H}_0}{\mathrm{FS}} = 2 \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}$$

Where we have recognized the spin’s irrelevance, by replacing the sum over $s$ with a factor $2$. Next, we turn the sum over the allowed $\vb{k}$-values into an integral, which is a common trick enabled by our periodic boundary conditions, yielding:

$$\begin{aligned} E^{(0)} &= \frac{2 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}} \end{aligned}$$

The matrix element $\matrixel{\mathrm{FS}}{\hat{c}_{\vb{k}}^\dagger \hat{c}_{\vb{k}}}{\mathrm{FS}} \dd{\vb{k}}$ is either $0$ or $1$, depending on whether $\vb{k}$ is outside the Fermi sphere, or, equivalently, whether $|\vb{k}|$ is above or below $k_F$. We can write this fact by introducing a Heaviside step function $\Theta(k)$:

$$\begin{aligned} E^{(0)} &= \frac{\hbar^2 V}{8 \pi^3 m} \int_{-\infty}^\infty |\vb{k}|^2 \: \Theta(k_F - |\vb{k}|) \dd{\vb{k}} \end{aligned}$$

We evaluate this in spherical coordinates 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}{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 / V$:

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

Traditionally, this is rewritten using the Wigner-Seitz radius $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_0^2 m)$:

$$\begin{aligned} \frac{4 \pi}{3} (a_0 r_s)^3 \equiv \frac{1}{n} \qquad \implies \qquad r_s = \bigg( \frac{3}{4 \pi a_0^3 n} \bigg)^{1/3} \end{aligned}$$

Note that this is dimensionless due to our choice of $a_0$ as a unit. In the Fermi gas, we have:

$$\begin{aligned} r_s = \bigg( \frac{9 \pi}{4} \bigg)^{1/3} \frac{1}{a_0 k_F} \qquad \implies \qquad k_F = \bigg( \frac{9 \pi}{4} \bigg)^{1/3} \frac{1}{a_0 r_s} \end{aligned}$$

By inserting this into the ground state energy and using the definition of $a_0$, we can write:

$$\begin{aligned} \frac{E^{(0)}}{N} &= \frac{3 \hbar^2}{10 m} \bigg( \frac{9 \pi}{4} \bigg)^{2/3} \frac{1}{a_0^2 r_s^2} \\ &= \frac{3}{5} \bigg( \frac{9 \pi}{4} \bigg)^{2/3} \bigg( \frac{e_0^2}{8 \pi \varepsilon_0 a_0} \bigg) \frac{1}{r_s^2} \end{aligned}$$

Where the last parenthesized expression is the Rydberg unit of energy $\mathrm{Ry} \approx 13.6 \:\mathrm{eV}$, so:

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

This result is found in a lot of literature. The choice of Rydberg units is simply a tradition.

1st order

In the next term of the perturbation expansion, we start to include Coulomb interactions. Clearly, this will give better results when the interaction is relatively weak, but is that ever the case?

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}$ scales as $n^{1/3}$. We also know that the kinetic energy $E_\mathrm{kin} = E^{(0)}$ is proportional to $n^{2/3}$, meaning that it is reasonable to use perturbation theory as long as $1 \gg E_\mathrm{int} / E_\mathrm{kin} \propto n^{-1/3}$, i.e. 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_0^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_0^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_0^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_0^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_0^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_0^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_0^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_0^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_0^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_0^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_0^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_0^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_0^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_0^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_0^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_0^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_0^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_0^2 V}{16 \pi^4 \varepsilon_0} k_F^4 = \frac{- e_0^2 N}{16 \pi^4 \varepsilon_0 n} k_F^4 = -\frac{3 e_0^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_0^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_0^2}{16 \pi^2 \varepsilon_0} \frac{a_0 k_F}{a_0} = -\frac{3 e_0^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.

---