Part II — Chapter 1

Bose-Einstein Condensation

When a gas of non-interacting bosons is cooled below a critical temperature $T_c$, a macroscopic fraction of the particles accumulates in the single-particle ground state. This phenomenon — Bose-Einstein condensation (BEC) — was predicted by Einstein in 1925 based on Bose’s quantum statistics for photons, and first realized experimentally in dilute atomic gases of $^{87}$Rb, $^{23}$Na, and $^{7}$Li in 1995 (Cornell, Wieman, and Ketterle; Nobel Prize 2001).

BEC is a purely quantum-statistical effect: it occurs even in the absence of interactions, driven entirely by the requirement that bosonic wave functions be symmetric under particle exchange. The thermal de Broglie wavelength $\lambda_{\mathrm{dB}} = \sqrt{2\pi\hbar^2/(mk_BT)}$becomes comparable to the interparticle spacing at $T_c$, signaling the onset of quantum degeneracy.

The Bose-Einstein Distribution

The mean occupation number of a single-particle state with energy $\epsilon$ is given by the Bose-Einstein distribution function:

$$\langle n(\epsilon) \rangle = \frac{1}{e^{(\epsilon - \mu)/k_B T} - 1}$$

where $\mu$ is the chemical potential. Unlike fermions, bosons have no exclusion principle, so $\langle n(\epsilon) \rangle$ can be arbitrarily large. However, the distribution requires $\mu \leq \epsilon_0$ (the ground-state energy, which we set to zero), since otherwise the ground-state occupation would be negative.

The total particle number is obtained by summing over all states:

$$N = \sum_{\mathbf{k}} \frac{1}{e^{(\epsilon_{\mathbf{k}} - \mu)/k_B T} - 1}$$

In three dimensions, converting the sum to an integral using the density of states$g(\epsilon) = (V/4\pi^2)(2m/\hbar^2)^{3/2}\epsilon^{1/2}$, one obtains the constraint that determines $\mu(T)$.

Chemical Potential and the Critical Temperature

As temperature decreases, the chemical potential $\mu$ rises toward zero to accommodate the fixed number of particles. At the critical temperature $T_c$,$\mu$ reaches its maximum value of zero. For $T \leq T_c$, the chemical potential is pinned at $\mu = 0$, and the excess particles that cannot be accommodated by the excited states condense into the ground state.

Above $T_c$ ($\mu < 0$)

The number equation determines $\mu(T)$ implicitly:

$$\frac{N}{V} = \frac{1}{\lambda_{\mathrm{dB}}^3}\,g_{3/2}(z), \quad z = e^{\mu/k_BT}$$

where $g_{3/2}(z) = \sum_{l=1}^{\infty} z^l / l^{3/2}$ is the polylogarithm (Bose function) and $z$ is the fugacity.

At $T_c$ ($\mu = 0$, $z = 1$)

Setting $z = 1$ gives $g_{3/2}(1) = \zeta(3/2) \approx 2.612$, yielding the critical temperature:

$$\boxed{T_c = \frac{2\pi\hbar^2}{mk_B}\left(\frac{n}{\zeta(3/2)}\right)^{2/3}}$$

This is equivalent to the condition $n\lambda_{\mathrm{dB}}^3 = \zeta(3/2)$: condensation occurs when the thermal wavelength equals the interparticle spacing.

Condensate Fraction

Below $T_c$, the number of particles in excited states is $N_{\mathrm{ex}} = N(T/T_c)^{3/2}$. The remaining particles occupy the ground state, giving the condensate fraction:

$$\boxed{\frac{N_0}{N} = 1 - \left(\frac{T}{T_c}\right)^{3/2}}, \quad T \leq T_c$$

At $T = 0$, all particles reside in the condensate ($N_0/N = 1$). The exponent $3/2$ reflects the three-dimensional density of states. In an interacting system such as liquid $^4$He, quantum depletion reduces$N_0/N$ to roughly 8% even at $T = 0$, though the superfluid fraction can still approach unity.

Key distinction: The condensate fraction $N_0/N$ and the superfluid fraction $\rho_s/\rho$ are different quantities. In the ideal Bose gas they coincide, but in interacting systems (like He-4) they differ substantially. The superfluid fraction is a transport property related to the system’s response to rotation.

Specific Heat and the Lambda Transition

The internal energy of the ideal Bose gas can be computed from the partition function. Below $T_c$:

$$\frac{E}{N k_B T_c} = \frac{3}{2}\,\frac{\zeta(5/2)}{\zeta(3/2)}\left(\frac{T}{T_c}\right)^{5/2}$$

Differentiating with respect to temperature gives the specific heat:

$$\frac{C_V}{Nk_B} = \frac{15}{4}\,\frac{\zeta(5/2)}{\zeta(3/2)}\left(\frac{T}{T_c}\right)^{3/2}, \quad T \leq T_c$$

The specific heat is continuous at $T_c$ but has a discontinuity in its derivative — a cusp resembling the Greek letter $\lambda$. In liquid$^4$He, the actual transition (the “lambda point” at $T_\lambda = 2.17$ K) shows a true logarithmic divergence in $C_V$, characteristic of the XY universality class rather than the ideal-gas result.

Above $T_c$, the specific heat approaches the classical equipartition value $C_V = \frac{3}{2}Nk_B$ as $T \to \infty$.

BEC in Lower Dimensions

The possibility of BEC depends critically on dimensionality through the density of states. In $d$ dimensions, the density of states scales as $g(\epsilon) \propto \epsilon^{d/2 - 1}$. The integral for the number of excited particles converges at $\mu = 0$ only if $d > 2$.

No BEC in 2D (Free Bosons)

In two dimensions, the density of states is constant ($g(\epsilon) \propto \epsilon^0$), and the integral for $N_{\mathrm{ex}}$ diverges logarithmically as $\mu \to 0$. The excited states can absorb all particles at any finite temperature, so no macroscopic ground-state occupation occurs. This is a consequence of the Mermin-Wagner-Hohenberg theorem.

Kosterlitz-Thouless Transition

Although true long-range order (and hence conventional BEC) is impossible in 2D, a topological phase transition — the Berezinskii-Kosterlitz-Thouless (BKT) transition — can occur. Below the BKT temperature, the system exhibits quasi-long-range order with algebraically decaying correlations $\langle \psi^*(\mathbf{r})\psi(0)\rangle \sim r^{-\eta}$, and the superfluid density jumps discontinuously to zero at $T_{\mathrm{BKT}}$.

The BKT transition is driven by the unbinding of vortex-antivortex pairs and belongs to a universality class entirely distinct from the 3D BEC transition. It has been observed in thin helium films and in quasi-2D ultracold atomic gases.

Trapped Gases and Effective BEC in 2D

In harmonic traps, the modified density of states can permit BEC-like transitions even in 2D. The trapping potential modifies the effective density of states, and for a 2D harmonic trap $T_c \propto \sqrt{N}\,\hbar\omega/k_B$. Experimentally, quasi-condensates with suppressed phase fluctuations have been realized in pancake-shaped optical traps.

Experimental Realizations

Ultracold Atomic Gases (1995)

The first BEC was achieved in June 1995 by Cornell and Wieman at JILA using $^{87}$Rb atoms cooled to approximately 170 nK. Shortly after, Ketterle’s group at MIT produced a condensate in $^{23}$Na with a much larger number of atoms, enabling the first observation of interference between two condensates.

These experiments employ a combination of laser cooling (Doppler and sub-Doppler techniques) to reach the microkelvin regime, followed by evaporative cooling in magnetic or optical traps to reach nanokelvin temperatures. Typical parameters:

Atom number: $N \sim 10^5$–$10^7$
Peak density: $n \sim 10^{13}$–$10^{15}$ cm$^{-3}$
Critical temperature: $T_c \sim 100$–$1000$ nK
Scattering length: $a \sim 5$–$100$ nm

Liquid Helium-4

Liquid $^4$He undergoes the superfluid transition at $T_\lambda = 2.17$ K (the lambda point). While this is driven by Bose statistics and is related to BEC, the strong interatomic interactions make the system fundamentally different from the ideal Bose gas:

Condensate fraction at $T = 0$: only $\sim 8\%$ (due to quantum depletion)
Superfluid fraction at $T = 0$: $\sim 100\%$
Excitation spectrum: phonon-roton curve (Landau), not free-particle dispersion
Critical exponents: XY universality class, not mean-field

The connection between BEC and superfluidity in He-4 remains a deep question. While BEC is neither necessary nor sufficient for superfluidity in general, in He-4 the condensate provides the order parameter $\langle \hat{\psi} \rangle = \sqrt{n_0}\,e^{i\phi}$ whose phase rigidity underlies superfluid flow.

Derivation: From Density of States to $T_c$

Step 1: Density of States in 3D

For free particles in a box of volume $V$, the single-particle energies are $\epsilon_{\mathbf{k}} = \hbar^2 k^2 / 2m$. The density of states is:

$$g(\epsilon) = \frac{V}{4\pi^2}\left(\frac{2m}{\hbar^2}\right)^{3/2}\epsilon^{1/2} = \frac{2\pi V}{h^3}(2m)^{3/2}\epsilon^{1/2}$$

Step 2: Number of Excited Particles

Separating the ground state from the integral:

$$N = N_0 + \int_0^{\infty} \frac{g(\epsilon)\,d\epsilon}{e^{(\epsilon - \mu)/k_BT} - 1}$$

The integral is maximized when $\mu = 0$. Using the substitution $x = \epsilon/k_BT$:

$$N_{\mathrm{ex}}^{\max} = \frac{V}{\lambda_{\mathrm{dB}}^3}\,\zeta(3/2), \quad \lambda_{\mathrm{dB}} = \sqrt{\frac{2\pi\hbar^2}{mk_BT}}$$

Step 3: Critical Temperature

At $T_c$, all $N$ particles can just be accommodated in excited states ($N_0 = 0$):

$$n = \frac{N}{V} = \frac{\zeta(3/2)}{\lambda_{\mathrm{dB}}^3(T_c)}$$

Solving for $T_c$:

$$T_c = \frac{2\pi\hbar^2}{mk_B}\left(\frac{n}{\zeta(3/2)}\right)^{2/3}$$

Detailed Derivation: 3D Density of States

Step 1: Counting States in k-Space

For particles in a cubic box of side $L$ with periodic boundary conditions, the allowed wavevectors are $\mathbf{k} = (2\pi/L)(n_x, n_y, n_z)$ with integer $n_i$. Each state occupies a volume$(2\pi/L)^3 = (2\pi)^3/V$ in k-space. The number of states with wavevector magnitude between $k$ and $k + dk$ is:

$$dN = \frac{V}{(2\pi)^3} \cdot 4\pi k^2\,dk$$

Since $\epsilon = \hbar^2 k^2 / (2m)$, we have $k = \sqrt{2m\epsilon}/\hbar$and $dk = \frac{1}{\hbar}\sqrt{\frac{m}{2\epsilon}}\,d\epsilon$. Substituting:

$$dN = \frac{V}{(2\pi)^3} \cdot 4\pi \cdot \frac{2m\epsilon}{\hbar^2} \cdot \frac{1}{\hbar}\sqrt{\frac{m}{2\epsilon}}\,d\epsilon = \frac{V}{4\pi^2}\left(\frac{2m}{\hbar^2}\right)^{3/2}\epsilon^{1/2}\,d\epsilon$$

Hence the density of states is $g(\epsilon) = \frac{V}{4\pi^2}\left(\frac{2m}{\hbar^2}\right)^{3/2}\epsilon^{1/2}$.

Step 2: Evaluating the Number Integral at $\mu = 0$

The total number of particles in excited states at $\mu = 0$ is:

$$N_{\mathrm{ex}} = \int_0^{\infty} \frac{g(\epsilon)\,d\epsilon}{e^{\epsilon/k_BT} - 1} = \frac{V}{4\pi^2}\left(\frac{2m}{\hbar^2}\right)^{3/2} \int_0^{\infty} \frac{\epsilon^{1/2}\,d\epsilon}{e^{\epsilon/k_BT} - 1}$$

Substitute $x = \epsilon/(k_B T)$, so $d\epsilon = k_B T\,dx$and $\epsilon^{1/2} = (k_B T)^{1/2} x^{1/2}$:

$$N_{\mathrm{ex}} = \frac{V}{4\pi^2}\left(\frac{2m k_B T}{\hbar^2}\right)^{3/2} \int_0^{\infty} \frac{x^{1/2}\,dx}{e^x - 1}$$

The integral is a standard result: $\int_0^{\infty} \frac{x^{1/2}}{e^x - 1}\,dx = \Gamma(3/2)\,\zeta(3/2) = \frac{\sqrt{\pi}}{2}\,\zeta(3/2)$. Recognizing that $\lambda_{\mathrm{dB}} = \sqrt{2\pi\hbar^2/(mk_BT)}$, we can write$(2mk_BT/\hbar^2)^{3/2} = (4\pi)^{3/2}/\lambda_{\mathrm{dB}}^3$, and after simplification:

$$N_{\mathrm{ex}} = \frac{V}{\lambda_{\mathrm{dB}}^3}\,\zeta(3/2)$$

Step 3: Deriving the Condensate Fraction

Below $T_c$, $\mu$ remains pinned at zero and the excited-state population scales as:

$$N_{\mathrm{ex}}(T) = \frac{V}{\lambda_{\mathrm{dB}}^3(T)}\,\zeta(3/2) \propto T^{3/2}$$

At $T = T_c$, $N_{\mathrm{ex}} = N$ by definition. Therefore, the ratio:

$$\frac{N_{\mathrm{ex}}(T)}{N} = \frac{N_{\mathrm{ex}}(T)}{N_{\mathrm{ex}}(T_c)} = \frac{\lambda_{\mathrm{dB}}^3(T_c)}{\lambda_{\mathrm{dB}}^3(T)} = \left(\frac{T}{T_c}\right)^{3/2}$$

The condensate number is $N_0 = N - N_{\mathrm{ex}}$, so:

$$\frac{N_0}{N} = 1 - \left(\frac{T}{T_c}\right)^{3/2}$$

Detailed Derivation: Specific Heat Below $T_c$

Step 1: Internal Energy Integral

The internal energy is obtained by weighting each state by its energy:

$$E = \int_0^{\infty} \frac{\epsilon\,g(\epsilon)\,d\epsilon}{e^{\epsilon/k_BT} - 1} = \frac{V}{4\pi^2}\left(\frac{2m}{\hbar^2}\right)^{3/2} \int_0^{\infty} \frac{\epsilon^{3/2}\,d\epsilon}{e^{\epsilon/k_BT} - 1}$$

With the same substitution $x = \epsilon/(k_BT)$:

$$E = \frac{V}{4\pi^2}\left(\frac{2m}{\hbar^2}\right)^{3/2} (k_BT)^{5/2} \int_0^{\infty} \frac{x^{3/2}\,dx}{e^x - 1}$$

The integral evaluates to $\Gamma(5/2)\,\zeta(5/2) = \frac{3\sqrt{\pi}}{4}\,\zeta(5/2)$. Using $N = V\,\zeta(3/2)/\lambda_{\mathrm{dB}}^3(T_c)$ and simplifying:

$$E = \frac{3}{2}\,N k_B T_c \,\frac{\zeta(5/2)}{\zeta(3/2)}\left(\frac{T}{T_c}\right)^{5/2}$$

Step 2: Differentiation to Obtain $C_V$

Taking the derivative $C_V = \partial E/\partial T$ at constant volume:

$$C_V = \frac{3}{2}\,N k_B\,\frac{\zeta(5/2)}{\zeta(3/2)} \cdot \frac{5}{2}\left(\frac{T}{T_c}\right)^{3/2} = \frac{15}{4}\,\frac{\zeta(5/2)}{\zeta(3/2)}\,N k_B \left(\frac{T}{T_c}\right)^{3/2}$$

Evaluating at $T = T_c$: $C_V/(Nk_B) = \frac{15}{4}\cdot\frac{1.3415}{2.6124} \approx 1.926$. Just above $T_c$, the classical limit gives $C_V/(Nk_B) = 3/2 = 1.5$, so there is a discontinuity in $\partial C_V/\partial T$ (a cusp in $C_V$ itself) — the hallmark lambda-shaped anomaly.

Derivation: BEC in a Harmonic Trap

Step 1: Density of States in a 3D Isotropic Harmonic Trap

In a harmonic potential $V = \frac{1}{2}m\omega^2 r^2$, the energy levels are$\epsilon_{n_x n_y n_z} = (n_x + n_y + n_z + 3/2)\hbar\omega$. For large quantum numbers, the number of states with energy less than $\epsilon$ is approximately the volume of the octant in $(n_x, n_y, n_z)$ space:

$$\mathcal{N}(\epsilon) \approx \frac{1}{6}\left(\frac{\epsilon}{\hbar\omega}\right)^3$$

Differentiating: $g(\epsilon) = d\mathcal{N}/d\epsilon = \frac{\epsilon^2}{2(\hbar\omega)^3}$.

Step 2: Critical Temperature in a Trap

Setting $\mu = 0$ and evaluating $N = \int_0^{\infty} g(\epsilon)/(e^{\epsilon/k_BT_c} - 1)\,d\epsilon$with $x = \epsilon/(k_BT_c)$:

$$N = \frac{(k_BT_c)^3}{2(\hbar\omega)^3}\int_0^{\infty}\frac{x^2\,dx}{e^x - 1} = \frac{(k_BT_c)^3}{2(\hbar\omega)^3}\cdot 2\zeta(3) = \frac{\zeta(3)(k_BT_c)^3}{(\hbar\omega)^3}$$

Solving for $T_c$ (using $\zeta(3) \approx 1.202$):

$$k_B T_c = \hbar\omega\left(\frac{N}{\zeta(3)}\right)^{1/3} \approx 0.94\,\hbar\omega\,N^{1/3}$$

Note the different power law compared to the uniform case: the condensate fraction now scales as $N_0/N = 1 - (T/T_c)^3$ rather than $(T/T_c)^{3/2}$, reflecting the different density of states in the trap.

Historical Context

Satyendra Nath Bose initiated the theory in 1924 by deriving Planck's radiation law using a novel counting method for indistinguishable photons. He sent his paper to Albert Einstein, who translated it into German, published it, and then extended Bose's statistics to massive particles. In a 1925 paper, Einstein predicted that below a critical temperature, a finite fraction of particles would “condense” into the lowest energy state — a prediction met with skepticism for decades.

Fritz London (1938) was the first to connect Einstein's prediction to the superfluid transition in liquid $^4$He, proposing that the lambda transition at 2.17 K was a manifestation of BEC modified by interactions. This idea was controversial at the time but proved essentially correct.

The experimental race to achieve BEC in dilute gases lasted decades. Key enabling technologies included laser cooling (Chu, Cohen-Tannoudji, Phillips; Nobel 1997) and evaporative cooling in magnetic traps. On June 5, 1995, Eric Cornell and Carl Wieman at JILA produced the first BEC in$^{87}$Rb with about 2000 atoms at 170 nK. Four months later, Wolfgang Ketterle at MIT created a much larger$^{23}$Na condensate, enabling the first observation of interference between two condensates. Cornell, Wieman, and Ketterle shared the 2001 Nobel Prize in Physics.

Applications of BEC

Atom interferometry and precision measurement: BECs provide coherent matter waves for atom interferometers that measure gravitational acceleration, rotations (Sagnac effect), and fundamental constants with extraordinary precision. These devices are used in inertial navigation and tests of the equivalence principle.
Quantum simulation: Ultracold atoms in optical lattices realize Hubbard models, enabling the study of quantum magnetism, Mott insulator transitions, and other strongly correlated phenomena that are computationally intractable on classical computers.
Atom lasers: A coherent beam of atoms extracted from a BEC is the matter-wave analogue of an optical laser. Atom lasers have potential applications in precision lithography and fundamental tests of quantum mechanics.
Slow light and quantum memory: BECs have been used to slow light pulses to speeds of metres per second and even to “stop” light entirely, storing optical information as atomic spin excitations — a key ingredient for quantum communication.
Tests of fundamental physics: BECs enable laboratory studies of analogue gravity (phononic Hawking radiation), Anderson localization of matter waves in disordered potentials, and the Kibble-Zurek mechanism for topological defect formation during phase transitions.

Simulation: BEC Thermodynamics

Condensate fraction, specific heat with the lambda-like cusp, and Bose-Einstein distribution at several temperatures.

Bose-Einstein Condensation: Thermodynamic Properties

Python

script.py122 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
# Precomputed values (no scipy needed)
def zeta(s):
    """Riemann zeta for s=3/2 and s=5/2 via known values."""
    vals = {1.5: 2.612375348685488, 2.5: 1.341487257250271}
    return vals.get(round(s, 1), 1.0)

# ─────────────────────────────────────────────────────────
# Bose-Einstein Condensation: Condensate Fraction,
# Specific Heat, and Bose Distribution
# ─────────────────────────────────────────────────────────

fig, axes = plt.subplots(1, 3, figsize=(18, 6))
fig.patch.set_facecolor('#0a0a1a')
for ax in axes:
    ax.set_facecolor('#0a0a1a')
    for spine in ax.spines.values():
        spine.set_edgecolor('#334155')
    ax.tick_params(colors='#94a3b8')

# ── Panel 1: Condensate fraction N_0/N vs T/T_c ──
t = np.linspace(0, 1.8, 500)
n0_frac = np.where(t <= 1.0, 1.0 - t**1.5, 0.0)

axes[0].fill_between(t, n0_frac, alpha=0.3, color='#06b6d4')
axes[0].plot(t, n0_frac, color='#22d3ee', linewidth=2.5, label=r'$N_0/N = 1 - (T/T_c)^{3/2}$')
axes[0].axvline(x=1.0, color='#f59e0b', linewidth=1.5, linestyle='--', alpha=0.7, label=r'$T_c$')
axes[0].set_xlabel(r'$T / T_c$', color='#94a3b8', fontsize=13)
axes[0].set_ylabel(r'$N_0 / N$', color='#94a3b8', fontsize=13)
axes[0].set_title('Condensate Fraction', color='white', fontsize=14, fontweight='bold')
axes[0].legend(fontsize=11, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[0].set_xlim(0, 1.8)
axes[0].set_ylim(-0.05, 1.05)
axes[0].text(0.4, 0.5, 'BEC\nphase', color='#22d3ee', fontsize=12,
             ha='center', va='center', alpha=0.7)
axes[0].text(1.35, 0.3, 'Normal\nphase', color='#f59e0b', fontsize=12,
             ha='center', va='center', alpha=0.7)

# ── Panel 2: Specific heat C_V vs T showing lambda cusp ──
zeta_3_2 = zeta(3/2)
zeta_5_2 = zeta(5/2)

t2 = np.linspace(0.01, 2.0, 600)
cv = np.zeros_like(t2)

for i, tt in enumerate(t2):
    if tt <= 1.0:
        # Below T_c: C_V = (15/4) * (V/lambda^3) * zeta(5/2) * k_B
        # In units of N*k_B: C_V/(N*k_B) = (15/4)*zeta(5/2)/zeta(3/2) * (T/T_c)^{3/2}
        cv[i] = (15.0/4.0) * zeta_5_2 / zeta_3_2 * tt**1.5
    else:
        # Above T_c: approaches classical 3/2 from above with correction
        # Approximate: smooth decay toward 3/2
        z_eff = np.exp(-0.8 * (tt - 1.0))  # approximate fugacity decay
        cv[i] = 1.5 + 1.5 * z_eff * (tt - 1.0) * np.exp(-1.2*(tt-1.0))
        # Better approximation: starts at cusp value and decays to 3/2
        cv_at_tc = (15.0/4.0) * zeta_5_2 / zeta_3_2
        cv[i] = 1.5 + (cv_at_tc - 1.5) * np.exp(-2.5 * (tt - 1.0))

axes[1].plot(t2, cv, color='#22d3ee', linewidth=2.5, label=r'$C_V / Nk_B$')
axes[1].axvline(x=1.0, color='#f59e0b', linewidth=1.5, linestyle='--', alpha=0.7, label=r'$T_c$')
axes[1].axhline(y=1.5, color='#64748b', linewidth=1, linestyle=':', alpha=0.5, label=r'Classical $\frac{3}{2}k_B$')
axes[1].set_xlabel(r'$T / T_c$', color='#94a3b8', fontsize=13)
axes[1].set_ylabel(r'$C_V / Nk_B$', color='#94a3b8', fontsize=13)
axes[1].set_title(r'Specific Heat ($\lambda$-like Cusp)', color='white', fontsize=14, fontweight='bold')
axes[1].legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[1].set_xlim(0, 2.0)
axes[1].set_ylim(0, 2.5)
axes[1].annotate(r'$\lambda$-point', xy=(1.0, cv[np.argmin(np.abs(t2-1.0))]),
                xytext=(1.4, 2.2), color='#f59e0b', fontsize=11,
                arrowprops=dict(arrowstyle='->', color='#f59e0b'),
                fontweight='bold')

# ── Panel 3: Bose-Einstein distribution at different temperatures ──
epsilon = np.linspace(0.01, 5.0, 500)  # in units of k_B * T_c

temperatures = [0.5, 1.0, 2.0, 5.0]
colors_be = ['#06b6d4', '#22d3ee', '#f59e0b', '#ef4444']

for T_ratio, color in zip(temperatures, colors_be):
    # mu = 0 at and below T_c, mu < 0 above T_c
    if T_ratio <= 1.0:
        mu = 0.0
    else:
        mu = -0.5 * (T_ratio - 1.0)  # approximate

n_be = 1.0 / (np.exp((epsilon - mu) / T_ratio) - 1.0)
    n_be = np.where(n_be > 0, n_be, np.nan)
    n_be = np.clip(n_be, 0, 10)
    axes[2].plot(epsilon, n_be, color=color, linewidth=2,
                 label=r'$T/T_c = %.1f$' % T_ratio)

# Classical Maxwell-Boltzmann for comparison
n_mb = np.exp(-epsilon / 2.0)
axes[2].plot(epsilon, n_mb * 2, color='#64748b', linewidth=1.5, linestyle='--',
             label='Maxwell-Boltzmann')

axes[2].set_xlabel(r'$\epsilon / k_B T_c$', color='#94a3b8', fontsize=13)
axes[2].set_ylabel(r'$\langle n(\epsilon) \rangle$', color='#94a3b8', fontsize=13)
axes[2].set_title('Bose-Einstein Distribution', color='white', fontsize=14, fontweight='bold')
axes[2].legend(fontsize=9, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[2].set_xlim(0, 5)
axes[2].set_ylim(0, 6)

plt.tight_layout(w_pad=3)
plt.savefig('bec_plots.png', dpi=130, bbox_inches='tight',
            facecolor=fig.get_facecolor())

print("Bose-Einstein Condensation Analysis")
print("=" * 50)
print(f"zeta(3/2) = {zeta_3_2:.6f}")
print(f"zeta(5/2) = {zeta_5_2:.6f}")
cv_at_tc = (15.0/4.0) * zeta_5_2 / zeta_3_2
print(f"C_V/(Nk_B) at T_c = {cv_at_tc:.4f}")
print(f"Classical limit C_V/(Nk_B) = 3/2 = 1.5000")
print(f"Discontinuity in slope of C_V at T_c (lambda-point)")
print(f"Condensate fraction at T=0: N_0/N = 1.0 (ideal gas)")
print(f"Condensate fraction at T=0.5*T_c: N_0/N = {1.0 - 0.5**1.5:.4f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Computation: Critical Temperatures and Condensate Fraction

Fortran computation of $T_c$ for different alkali species (Rb-87, Na-23, Li-7) and tabulation of condensate fraction vs temperature.

BEC Critical Temperatures for Alkali Atoms

Fortran

program.f90133 lines

program bose_einstein_condensation
  implicit none

! ─────────────────────────────────────────────────────────
  ! Compute T_c for different atomic species and
  ! condensate fraction vs temperature
  ! ─────────────────────────────────────────────────────────

double precision, parameter :: pi = 3.14159265358979323846d0
  double precision, parameter :: hbar = 1.0545718d-34    ! J s
  double precision, parameter :: kB = 1.380649d-23       ! J/K
  double precision, parameter :: amu = 1.66053906660d-27  ! kg
  double precision, parameter :: zeta_3_2 = 2.612375348d0 ! Riemann zeta(3/2)

! Atomic species: Rb-87, Na-23, Li-7
  integer, parameter :: n_species = 3
  character(len=8) :: names(n_species)
  double precision :: masses(n_species)   ! in amu
  double precision :: densities(n_species) ! typical peak densities in m^-3

double precision :: m_kg, n_dens, Tc, prefactor
  double precision :: T, T_ratio, N0_frac
  integer :: i, j, n_temps

names(1) = 'Rb-87'
  names(2) = 'Na-23'
  names(3) = 'Li-7'

masses(1) = 86.909180527d0
  masses(2) = 22.9897692820d0
  masses(3) = 7.01600343660d0

! Typical experimental peak densities (atoms/m^3)
  densities(1) = 2.5d18
  densities(2) = 1.5d19
  densities(3) = 3.0d18

write(*,'(A)') '============================================================'
  write(*,'(A)') '  Bose-Einstein Condensation: Critical Temperatures'
  write(*,'(A)') '============================================================'
  write(*,'(A)') ''
  write(*,'(A)') '  Formula: T_c = (2*pi*hbar^2)/(m*k_B) * (n/zeta(3/2))^(2/3)'
  write(*,'(A)') ''
  write(*,'(A,F12.9)') '  zeta(3/2) = ', zeta_3_2
  write(*,'(A)') ''
  write(*,'(A8,A14,A18,A14)') 'Species', 'Mass [amu]', 'Density [m^-3]', 'T_c [nK]'
  write(*,'(A)') '  --------------------------------------------------------'

do i = 1, n_species
    m_kg = masses(i) * amu
    n_dens = densities(i)
    prefactor = (2.0d0 * pi * hbar**2) / (m_kg * kB)
    Tc = prefactor * (n_dens / zeta_3_2)**(2.0d0/3.0d0)

write(*,'(A8, F14.4, ES18.3, F14.2)') &
      names(i), masses(i), n_dens, Tc * 1.0d9
  end do

! ── Condensate fraction vs temperature for Rb-87 ──
  write(*,'(A)') ''
  write(*,'(A)') '============================================================'
  write(*,'(A)') '  Condensate Fraction vs Temperature (Rb-87)'
  write(*,'(A)') '============================================================'
  write(*,'(A)') ''
  write(*,'(A)') '  N_0/N = 1 - (T/T_c)^(3/2)   for T < T_c'
  write(*,'(A)') '  N_0/N = 0                     for T >= T_c'
  write(*,'(A)') ''

! Compute T_c for Rb-87
  m_kg = masses(1) * amu
  n_dens = densities(1)
  prefactor = (2.0d0 * pi * hbar**2) / (m_kg * kB)
  Tc = prefactor * (n_dens / zeta_3_2)**(2.0d0/3.0d0)

write(*,'(A,F10.2,A)') '  T_c (Rb-87) = ', Tc * 1.0d9, ' nK'
  write(*,'(A)') ''
  write(*,'(A10,A14,A14)') 'T/T_c', 'T [nK]', 'N_0/N'
  write(*,'(A)') '  ------------------------------------'

n_temps = 20
  do j = 0, n_temps
    T_ratio = dble(j) / dble(n_temps) * 1.5d0
    T = T_ratio * Tc

if (T_ratio <= 1.0d0) then
      N0_frac = 1.0d0 - T_ratio**1.5d0
    else
      N0_frac = 0.0d0
    end if

write(*,'(F10.4, F14.2, F14.6)') T_ratio, T * 1.0d9, N0_frac
  end do

! ── Thermal de Broglie wavelength comparison ──
  write(*,'(A)') ''
  write(*,'(A)') '============================================================'
  write(*,'(A)') '  Thermal de Broglie Wavelength at T_c'
  write(*,'(A)') '============================================================'
  write(*,'(A)') ''
  write(*,'(A)') '  lambda_dB = sqrt(2*pi*hbar^2 / (m*k_B*T))'
  write(*,'(A)') '  BEC onset: n * lambda_dB^3 = zeta(3/2) ~ 2.612'
  write(*,'(A)') ''

do i = 1, n_species
    m_kg = masses(i) * amu
    n_dens = densities(i)
    prefactor = (2.0d0 * pi * hbar**2) / (m_kg * kB)
    Tc = prefactor * (n_dens / zeta_3_2)**(2.0d0/3.0d0)

! lambda_dB at T_c
    call compute_lambda(m_kg, Tc, prefactor)
    write(*,'(A8,A,F8.2,A,A,ES10.3,A)') &
      names(i), ': lambda_dB(T_c) = ', &
      dsqrt(2.0d0*pi*hbar**2/(m_kg*kB*Tc)) * 1.0d6, ' um', &
      ', n*lambda^3 = ', &
      n_dens * (dsqrt(2.0d0*pi*hbar**2/(m_kg*kB*Tc)))**3, &
      ' (should ~ 2.612)'
  end do

write(*,'(A)') ''
  write(*,'(A)') '============================================================'
  write(*,'(A)') '  Computation complete.'
  write(*,'(A)') '============================================================'

contains

subroutine compute_lambda(mass, temp, result)
    double precision, intent(in) :: mass, temp
    double precision, intent(out) :: result
    result = dsqrt(2.0d0 * pi * hbar**2 / (mass * kB * temp))
  end subroutine compute_lambda

end program bose_einstein_condensation

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Summary and Key Results

Distribution

$$\langle n(\epsilon)\rangle = \frac{1}{e^{(\epsilon-\mu)/k_BT}-1}$$

Critical Temperature

$$T_c = \frac{2\pi\hbar^2}{mk_B}\!\left(\frac{n}{\zeta(3/2)}\right)^{\!2/3}$$

Condensate Fraction

$$N_0/N = 1 - (T/T_c)^{3/2}$$

Onset Condition

$$n\lambda_{\mathrm{dB}}^3 = \zeta(3/2) \approx 2.612$$

← Part II Overview Chapter 2: Superfluid Helium →