Bose-Einstein Condensation Simulation

Bose-Einstein Condensation

Cool a trapped gas of N bosons in a 3D harmonic trap and watch a macroscopic fraction pile into the ground state below T_c — computed with an exact discrete sum over energy levels, not the continuum approximation.

Temperature (T / T_c)

Particle Number (N)

T / T_c

2.00

μ / ħω

N₀ / N

N_excited / N

Physics

$$E_n = n\hbar\omega, \qquad g_n = \frac{(n+1)(n+2)}{2}$$

$$\langle N_n \rangle = \frac{g_n}{e^{(E_n-\mu)/k_BT}-1}, \qquad N = \sum_{n=0}^{n_{\max}} \langle N_n \rangle$$

n: principal quantum number, n = n_x+n_y+n_z ≥ 0
ħω: trap energy quantum (ħ = reduced Planck's constant, ω = trap angular frequency)
E_n: energy of level n, measured from the ground state (E₀=0)
g_n: degeneracy of level n — how many (n_x,n_y,n_z) triples sum to n
μ: chemical potential — not set directly, but solved for (see "How to Calculate μ" below)
k_BT: thermal energy scale (k_B = Boltzmann constant, T = temperature)
⟨N_n⟩: expected number of particles occupying level n
N, n_max: total particle number (set by the slider); the level where the sum is cut off

How to Calculate μ?

μ is never set directly. For the current N and T, it's whatever value makes the discrete sum add up to exactly N:

$$N(\mu) = \sum_{n=0}^{n_{\max}} \frac{g_n}{e^{(n-\mu)/T}-1} \;\stackrel{!}{=}\; N_{\text{target}}$$

This has no closed-form inverse, so it's solved numerically. The key fact that makes it easy: every term increases as μ increases, so N(μ) is strictly increasing on (−∞, 0). It vanishes as μ → −∞, and diverges as μ → 0⁻ (the n=0 term alone blows up). So for any target N, exactly one μ solves the equation — and bisection finds it:

1. Start with lo = a safely negative guess, hi = −10⁻⁹ (just below 0).
2. If N(lo) is still > target, lo isn't negative enough — double it and repeat.
3. Repeat ~100×: check the midpoint m = (lo+hi)/2.
If N(m) < target, the answer is above m → set lo = m.
Otherwise the answer is below m → set hi = m.
4. After 100 halvings the bracket [lo, hi] has shrunk to machine precision — return μ = (lo+hi)/2.

Once μ is known, plugging it back into the n=0 term alone gives the condensate population N₀, and N−N₀ gives everything in the excited states.

Why Discrete Sum?

Textbook treatments often replace the sum over energy levels with a continuous integral over a density of states g(ε). That's a fine approximation for the excited states, but it silently discards the single quantum state at ε=0 — which has zero measure in a continuous integral, yet is exactly the state that becomes macroscopically occupied below T_c.

This is the exact sum the sliders on this page evaluate, term by term

$$N = \sum_{n=0}^{n_{\max}} \frac{g_n}{e^{(n-\mu)/T}-1}$$ $$= \underbrace{\frac{1}{e^{-\mu/T}-1}}_{n=0 \text{ term} \,=\, N_0} \;+\; \sum_{n=1}^{n_{\max}} \frac{g_n}{e^{(n-\mu)/T}-1}$$

This is what the continuum approximation uses instead

$$N \approx \int_0^{\infty} g(\epsilon)\,\frac{d\epsilon}{e^{(\epsilon-\mu)/T}-1}$$ $$g(\epsilon) = \frac{\epsilon^2}{2}$$

Notice the integral has no separate n=0 term at all — a single point contributes exactly zero to a ∫dε, so g(0)=0 there too. That's why the continuum approach needs a manually-inserted patch: "if μ would need to go positive, clamp it to 0 and dump the leftover particles into N₀ by hand."

The discrete sum never needs that patch. Its n=0 term, isolated above, diverges on its own as μ → 0⁻, so the ground state naturally absorbs however many particles the excited-state terms (n≥1) can't hold — condensation "just happens" as a consequence of the same formula, at every temperature.

Condensate Fraction vs. Temperature

Momentum-Space Density

Energy Level Populations (Harmonic Oscillator Ladder)