Part II · Chapter 2

Superfluid Helium

1. The Lambda Transition

Liquid $^4$He undergoes a remarkable phase transition at the lambda point, $T_\lambda = 2.17\;\text{K}$, at saturated vapour pressure. The name arises from the characteristic shape of the specific heat curve $C(T)$, which resembles the Greek letter $\lambda$. Above $T_\lambda$, the liquid is called He-I and behaves as a classical viscous fluid. Below $T_\lambda$, it enters the He-II phase and exhibits superfluidity.

The transition is continuous (second-order) and belongs to the three-dimensional XY universality class. The order parameter is a complex scalar field$\Psi(\mathbf{r}) = |\Psi|e^{i\phi}$, where $|\Psi|^2$ is proportional to the superfluid density. The critical exponent for the specific heat divergence is $\alpha \approx -0.0127$, measured with extraordinary precision in the space-shuttle experiment (Lipa et al., 2003).

The specific heat near the lambda point follows:
$$C(T) \sim A_\pm |t|^{-\alpha}, \quad t = \frac{T - T_\lambda}{T_\lambda}$$
where $A_+/A_- \approx 1.054$ is the amplitude ratio and $\alpha \approx -0.0127$. The negative exponent means $C$ develops a cusp rather than a true divergence.

At atmospheric pressure, $^4$He remains liquid down to absolute zero — it never solidifies unless pressurized above $\sim 25\;\text{atm}$. This is a direct consequence of the large zero-point energy relative to the weak van der Waals binding, a quantum mechanical effect that makes helium the only element to exhibit this behaviour.

2. Properties of He-II

Zero Viscosity Below Critical Velocity

He-II flows through narrow capillaries and porous media without any measurable viscous drag, provided the flow velocity remains below a critical velocity $v_c$. In Kapitza’s classic experiment (1938), liquid helium flowed through a thin gap between glass discs with no detectable pressure drop. Yet simultaneously, a torsional oscillator immersed in He-II measures a finite viscosity from the normal component. This paradox is resolved by the two-fluid model(Chapter 3).

Thermal Superconductivity

The thermal conductivity of He-II is effectively infinite — orders of magnitude larger than that of copper. Heat transport occurs not by diffusion but by convective counterflow: the superfluid component flows toward the heat source while the normal component (carrying entropy) flows away. This mechanism makes it impossible to maintain a temperature gradient in bulk He-II, so the liquid cannot boil in the ordinary sense. Instead, evaporation occurs only at the free surface.

Film Flow and the Rollin Film

He-II spontaneously creeps up the walls of its container as a thin film (thickness$\sim 30\;\text{nm}$), driven by van der Waals attraction. If the container is partially immersed in a bath, the film continuously flows over the rim until the levels equalise. This Rollin film dramatically demonstrates the absence of viscosity for the superfluid component.

Key Numbers for He-II

Lambda temperature: $T_\lambda = 2.1768\;\text{K}$
Number density: $n = 2.18 \times 10^{28}\;\text{m}^{-3}$
Atomic mass: $m = 6.646 \times 10^{-27}\;\text{kg}$
Speed of first sound: $c_1 \approx 238\;\text{m/s}$
Speed of second sound (at $1.5\;\text{K}$): $c_2 \approx 20\;\text{m/s}$
Quantum of circulation: $\kappa = h/m = 9.97 \times 10^{-8}\;\text{m}^2\text{/s}$

3. The Landau Criterion for Superfluidity

Landau’s 1941 argument for superfluidity is rooted in energy-momentum conservation. Consider a superfluid flowing through a capillary at velocity $\mathbf{v}$. In the rest frame of the fluid, the capillary walls move at $-\mathbf{v}$. For dissipation to occur, the walls must create an excitation of energy $\epsilon(p)$and momentum $\mathbf{p}$.

Energy conservation in the laboratory frame requires:

$$E' = E + \epsilon(\mathbf{p}) + \mathbf{p}\cdot\mathbf{v}$$
For this to be energetically favourable we need $\epsilon(\mathbf{p}) + \mathbf{p}\cdot\mathbf{v} < 0$, which is minimised when $\mathbf{p}$ is antiparallel to $\mathbf{v}$. Thus dissipation is only possible when:
$$v > v_c = \min_p \frac{\epsilon(p)}{p}$$

This is the celebrated Landau criterion. If the excitation spectrum has no states with arbitrarily small $\epsilon/p$, the fluid is superfluid below $v_c$. For a free-particle spectrum$\epsilon = p^2/2m$, one obtains $v_c = 0$ — no superfluidity. This is why the form of the excitation spectrum is crucial.

For He-II, the minimum of $\epsilon(p)/p$ occurs at the roton minimum, giving a critical velocity$v_c \approx 58\;\text{m/s}$. Experimentally, observed critical velocities are typically much lower ($\sim 1\;\text{cm/s}$ to $\sim 1\;\text{m/s}$) because dissipation is dominated by vortex nucleation rather than roton creation.

4. The Excitation Spectrum: Phonons and Rotons

The elementary excitations of He-II form a single continuous branch known as the phonon-roton curve. At small wavevectors, the spectrum is linear (phonon-like):

$$\epsilon(k) \approx \hbar c_s k \quad \text{(phonon region, } k \lesssim 0.8\;\text{Å}^{-1}\text{)}$$
where $c_s \approx 238\;\text{m/s}$ is the speed of first sound. At higher wavevectors, the spectrum passes through a maximum and then descends to a local minimum near $k_0 \approx 1.92\;\text{Å}^{-1}$ — the roton minimum:
$$\epsilon(k) \approx \Delta + \frac{\hbar^2(k - k_0)^2}{2\mu}$$
with roton gap $\Delta/k_B \approx 8.6\;\text{K}$ and effective mass$\mu \approx 0.16\,m_{\text{He}}$.

The physical nature of the roton has been debated since Landau introduced it. Feynman showed that rotons correspond to highly localised, quantized backflow patterns — essentially microscopic vortex rings. Neutron scattering experiments by Henshaw and Woods (1961) provided the first direct measurement of the full dispersion curve, confirming Landau’s prediction.

The Landau-Feynman theory predicts the excitation spectrum from the static structure factor $S(k)$:

$$\epsilon(k) = \frac{\hbar^2 k^2}{2m\,S(k)}$$
Since $S(k)$ has a peak at $k \approx 2\;\text{Å}^{-1}$(reflecting the short-range order of the liquid), the energy has a corresponding dip — the roton minimum. This elegant result connects the static structure of the liquid directly to its dynamical excitations.

5. Fountain Effect and Thermomechanical Effect

In the two-fluid picture, a temperature difference across a superleak (a porous plug that passes only the superfluid component) drives a pressure difference. This is the thermomechanical effect, governed by the London-Tisza relation:

$$\nabla P = \rho\,s\,\nabla T$$
where $\rho$ is the total density and $s$ is the specific entropy. A temperature gradient drives superfluid flow toward the warmer region.

The fountain effect is its most spectacular manifestation. When a vessel packed with fine powder (the superleak) is heated, superfluid rushes in through the powder to equalise the chemical potential. The resulting pressure build-up can drive a jet of helium several centimetres into the air — a fountain powered solely by a temperature difference of a fraction of a kelvin.

The inverse process also occurs: forcing superfluid through a superleak cools the source reservoir and heats the receiving one (the mechanocaloric effect). These effects provided early evidence that He-II carries zero entropy in its superfluid component.

6. Superfluid Fraction vs Temperature

The superfluid density $\rho_s(T)$ vanishes at $T_\lambda$ and equals the total density at $T = 0$. Defining the superfluid fraction$f_s = \rho_s/\rho$, the two-fluid constraint is:

$$\rho = \rho_s(T) + \rho_n(T), \qquad f_s + f_n = 1$$
Near $T_\lambda$, $\rho_s$ vanishes as a power law:
$$\rho_s(T) \propto (T_\lambda - T)^{\zeta}, \quad \zeta \approx 0.6717$$
where $\zeta = \nu \approx 2/3$ is the correlation length exponent of the 3D XY model.

At low temperatures, the normal-fluid density arises from thermally excited phonons and rotons. The phonon contribution dominates below $\sim 0.6\;\text{K}$:

$$\rho_n^{\text{phonon}} = \frac{2\pi^2}{45} \frac{(k_B T)^4}{\hbar^3 c_s^5}$$
while the roton contribution dominates for $1\;\text{K} \lesssim T \lesssim T_\lambda$:
$$\rho_n^{\text{roton}} = \frac{2k_0^4}{3(2\pi)^{3/2}} \left(\frac{\hbar^2}{\mu k_B T}\right)^{1/2} e^{-\Delta/(k_B T)}$$

Experimental measurements of $\rho_s/\rho$ via the Andronikashvili torsional oscillator method show excellent agreement with these predictions. At $T = 0$, the condensate fraction (fraction of atoms in the zero-momentum state) is only about 7% due to strong interactions, yet the superfluid fraction is 100% — a striking distinction between condensation and superfluidity.

7. Helium-3: A Fermionic Superfluid

$^3$He atoms are fermions (nuclear spin $I = 1/2$) and cannot undergo Bose-Einstein condensation directly. Instead, they form Cooper pairs at millikelvin temperatures, analogous to electrons in BCS superconductors. The critical temperature is$T_c \approx 2.5\;\text{mK}$ at 34 bar — three orders of magnitude below the $^4$He lambda point.

The pairing is in the p-wave channel ($L = 1$) with spin-triplet ($S = 1$), giving a total angular momentum that can couple in multiple ways. This leads to multiple superfluid phases:

A-phase (ABM state): an equal-spin-pairing state with $\hat{d} \cdot \hat{l}$ order. The gap has point nodes along $\hat{l}$, and the phase spontaneously breaks time-reversal symmetry.
B-phase (BW state): an isotropic p-wave state with full $J = 0$ pairing ($L = 1, S = 1$ coupled to$J = 0$). The gap is isotropic in magnitude, and this phase fills most of the low-temperature, low-pressure phase diagram.
A$_1$-phase: appears in a magnetic field; only one spin species pairs.

The order parameter of superfluid $^3$He is a $3\times 3$ complex matrix $A_{\mu i}$ relating spin ($\mu$) and orbital ($i$) degrees of freedom:

$$\Delta_{\mu\nu}(\mathbf{k}) = i(\sigma_\mu \sigma_y)_{\mu\nu}\,\sum_i A_{\mu i}\,\hat{k}_i$$
The B-phase has $A_{\mu i} = \Delta_B\,R_{\mu i}(\hat{n},\theta)$ with a rotation matrix, while the A-phase has $A_{\mu i} = \Delta_A\,\hat{d}_\mu(\hat{m}_i + i\hat{n}_i)$.

Superfluid $^3$He is one of the richest systems in condensed matter physics. The A-phase is a topological superfluid hosting Majorana fermions at its surface, and the B-phase features a fully gapped topological order classified by an integer invariant ($\mathbb{Z}$). The 1996 Nobel Prize was awarded to Lee, Osheroff, and Richardson for the discovery of superfluidity in $^3$He.

Detailed Derivation: Landau Critical Velocity

Step 1: Galilean Transformation of Energy

Consider a superfluid of mass $M$ flowing through a tube at velocity $\mathbf{v}$. In the fluid rest frame, the fluid is at rest and the tube walls move at $-\mathbf{v}$. Suppose the walls create an excitation with momentum $\mathbf{p}$ and energy $\epsilon(p)$ in the fluid frame. In the lab frame, the fluid momentum changes from $M\mathbf{v}$ to $M\mathbf{v} + \mathbf{p}$, so the new kinetic energy is:

$$E_{\text{lab}} = \frac{(M\mathbf{v} + \mathbf{p})^2}{2M} + \epsilon(p) = \frac{Mv^2}{2} + \mathbf{p}\cdot\mathbf{v} + \frac{p^2}{2M} + \epsilon(p)$$

Step 2: Energy Change and Dissipation Condition

The energy change due to creating the excitation is (dropping the negligible $p^2/(2M)$ recoil term for macroscopic $M$):

$$\Delta E = \epsilon(p) + \mathbf{p}\cdot\mathbf{v}$$

Dissipation is energetically favourable only if $\Delta E < 0$. This is minimized when$\mathbf{p}$ is antiparallel to $\mathbf{v}$, giving:

$$\Delta E_{\min} = \epsilon(p) - pv$$

For dissipation to occur, we need $v > \epsilon(p)/p$ for at least some $p$. Hence superfluidity persists below the Landau critical velocity:

$$v_c = \min_p \frac{\epsilon(p)}{p}$$

Step 3: Application to He-4 Excitation Spectrum

For the phonon-roton spectrum, $\epsilon(p)/p$ must be evaluated at the roton minimum. Near the roton minimum at $p_0 = \hbar k_0$, the spectrum is$\epsilon(p) = \Delta + (p - p_0)^2/(2\mu)$. Setting $d(\epsilon/p)/dp = 0$:

$$\frac{d}{dp}\left(\frac{\epsilon}{p}\right) = \frac{p\,\epsilon'(p) - \epsilon(p)}{p^2} = 0 \quad \Longrightarrow \quad \epsilon'(p_{\min}) = \frac{\epsilon(p_{\min})}{p_{\min}}$$

At the roton minimum itself ($p = p_0$, where $\epsilon' = 0$), the ratio$\epsilon/p = \Delta/p_0$. Numerically, with $\Delta/k_B = 8.62$ K and$p_0/\hbar = 1.92$ $\text{\AA}^{-1}$:

$$v_c = \frac{\Delta}{p_0} = \frac{8.62 \times 1.381 \times 10^{-23}}{1.92 \times 10^{10} \times 1.055 \times 10^{-34}} \approx 58\;\text{m/s}$$

Detailed Derivation: Phonon Contribution to Normal Fluid Density

Step 1: Landau's Formula for Normal Fluid Density

Landau derived the normal fluid density by considering the momentum carried by thermal excitations. In a frame where the normal fluid moves at velocity $\mathbf{v}_n$, the momentum density of the excitation gas is $\rho_n \mathbf{v}_n$. This gives:

$$\rho_n = -\frac{1}{3}\int \frac{p^2}{\epsilon(p)}\,\frac{\partial n_B}{\partial \epsilon}\,\frac{4\pi p^2\,dp}{(2\pi\hbar)^3}$$

where $n_B(\epsilon) = (e^{\epsilon/k_BT} - 1)^{-1}$ is the Bose distribution.

Step 2: Evaluating for Phonons ($\epsilon = c_s p$)

For phonons, $\epsilon = c_s p$, so $p^2/\epsilon = p/c_s$. Using$-\partial n_B/\partial\epsilon = (1/k_BT)\,e^{\epsilon/k_BT}/(e^{\epsilon/k_BT}-1)^2$and substituting $x = c_s p/(k_BT)$:

$$\rho_n^{\text{ph}} = \frac{1}{3} \cdot \frac{4\pi}{(2\pi\hbar)^3} \cdot \frac{1}{c_s k_BT} \int_0^{\infty} p^3 \,\frac{e^{c_sp/(k_BT)}}{(e^{c_sp/(k_BT)}-1)^2}\,dp$$

Substituting $p = k_BT\,x/c_s$ and $dp = (k_BT/c_s)\,dx$:

$$\rho_n^{\text{ph}} = \frac{4\pi}{3(2\pi\hbar)^3}\,\frac{(k_BT)^4}{c_s^5}\int_0^{\infty}\frac{x^3 e^x}{(e^x - 1)^2}\,dx$$

The integral evaluates to $4!\,\zeta(4)\cdot(1 - 2^{-3}) = 4\pi^4/15$. After simplification:

$$\rho_n^{\text{ph}} = \frac{2\pi^2}{45}\,\frac{(k_BT)^4}{\hbar^3 c_s^5}$$

This $T^4$ scaling dominates below about 0.6 K, where roton contributions are exponentially suppressed.

Derivation: Feynman's Variational Excitation Spectrum

Step 1: Trial Wave Function

Feynman (1954) proposed a trial excited state of the form:

$$|\mathbf{k}\rangle = \frac{1}{\sqrt{N S(k)}}\,\rho_{\mathbf{k}}\,|0\rangle, \qquad \rho_{\mathbf{k}} = \sum_{j=1}^{N} e^{i\mathbf{k}\cdot\mathbf{r}_j}$$

where $|0\rangle$ is the ground state and $\rho_{\mathbf{k}}$ is the density fluctuation operator. The normalization involves the static structure factor$S(k) = \langle 0|\rho_{-\mathbf{k}}\rho_{\mathbf{k}}|0\rangle/N$.

Step 2: Variational Energy

Computing $\epsilon(k) = \langle\mathbf{k}|H - E_0|\mathbf{k}\rangle$ using the f-sum rule$\sum_n (\epsilon_n - E_0)|\langle n|\rho_{\mathbf{k}}|0\rangle|^2 = N\hbar^2 k^2/(2m)$:

$$\epsilon(k) = \frac{\hbar^2 k^2}{2m\,S(k)}$$

This is an upper bound on the true excitation energy. Since $S(k)$ has a prominent peak near $k \approx 2\;\text{\AA}^{-1}$ (reflecting the short-range order of the liquid, analogous to the first diffraction peak in the pair correlation function), the energy $\epsilon(k)$has a corresponding dip — the roton minimum. At small $k$, $S(k) \to \hbar k/(2mc_s)$(compressibility sum rule), recovering the phonon branch $\epsilon(k) = \hbar c_s k$.

Historical Context

Superfluidity in helium was discovered independently by Pyotr Kapitsain Moscow and by John F. Allen and Don Misener in Cambridge, both publishing in January 1938. Kapitsa coined the term “superfluidity” by analogy with superconductivity. He received the 1978 Nobel Prize in Physics for this work, though the omission of Allen and Misener remains controversial.

Fritz London (1938) made the conceptual leap connecting superfluidity to Bose-Einstein condensation, proposing that the lambda transition was a BEC modified by interactions. Laszlo Tisza (1938) introduced the two-fluid model as a phenomenological description, proposing that the superfluid component carried zero entropy. Lev Landau (1941) placed the theory on a rigorous foundation with his excitation spectrum approach and the critical velocity criterion, receiving the 1962 Nobel Prize. Landau initially rejected London's BEC connection, though modern theory vindicates London's insight.

Richard Feynman (1953–1955) provided microscopic understanding by deriving the excitation spectrum from the structure factor and explaining quantized vortices. The phonon-roton spectrum was directly confirmed by Henshaw and Woods (1961) using inelastic neutron scattering at Chalk River, Canada. The discovery of superfluidity in $^3$He by Lee, Osheroff, and Richardson (1972; Nobel 1996) opened an entirely new chapter in the field.

Applications of Superfluid Helium

Particle accelerator cooling: Superfluid $^4$He at 1.9 K is used to cool the superconducting magnets of the Large Hadron Collider (LHC) at CERN. Its extraordinary thermal conductivity prevents local hot spots that would quench the magnets. The LHC contains about 96 tonnes of superfluid helium.
Space-based infrared telescopes: The Infrared Astronomical Satellite (IRAS), Spitzer Space Telescope, and Herschel Space Observatory all used superfluid helium cryostats to cool their detectors to below 2 K, suppressing thermal noise for sensitive infrared measurements.
Precision tests of fundamental physics: The superfluid helium gyroscope (using persistent superfluid flow) achieves rotation sensitivity comparable to laser ring gyroscopes. The GP-B experiment used superfluid helium to cool its gyroscopes for testing frame-dragging predicted by general relativity.
Quantum turbulence research: Superfluid helium is the premier system for studying quantum turbulence — turbulence composed entirely of quantized vortex lines. This provides insight into classical turbulence phenomena such as the Kolmogorov energy cascade, but with the added constraint of quantized circulation.
Neutron moderation: Superfluid helium is used as an ultracold neutron source. Neutrons scatter from phonons in the superfluid, losing energy until they reach ultra-low energies ($< 300$ neV), enabling precision measurements of the neutron electric dipole moment and neutron lifetime.

Simulation 1: He-4 Excitation Spectrum (Python)

This simulation computes and plots the phonon-roton dispersion curve for superfluid$^4$He using the empirical parametrisation, then evaluates the Landau critical velocity from $v_c = \min(\epsilon(p)/p)$.

He-4 Phonon-Roton Excitation Spectrum & Critical Velocity

Python

script.py92 lines

import numpy as np

# ── Physical constants ──
hbar = 1.0546e-34       # J·s
kB   = 1.3806e-23       # J/K
m_He = 6.6464e-27       # kg (He-4 mass)
angstrom = 1e-10        # m

# ── Roton parameters (Donnelly & Barenghi) ──
c_s    = 238.3           # speed of first sound, m/s
k0     = 1.92 / angstrom # roton minimum wavevector, 1/m
Delta  = 8.62 * kB       # roton gap in Joules
mu_rot = 0.16 * m_He     # roton effective mass

# ── Build the dispersion curve ──
# Use empirical Landau-type parametrisation
k_arr = np.linspace(0.01, 3.5, 800) / angstrom   # wavevector in 1/m
k_disp = k_arr * angstrom                          # for display in 1/Angstrom

# Phonon branch at low k
eps_phonon = hbar * c_s * k_arr

# Roton parabola near minimum
eps_roton = Delta + hbar**2 * (k_arr - k0)**2 / (2.0 * mu_rot)

# Smooth interpolation: use the lower of the two branches
# with a smooth blending function
k_cross = 0.9 / angstrom  # approximate crossover region
sigma = 0.3 / angstrom
weight = 0.5 * (1.0 + np.tanh((k_arr - k_cross) / sigma))
eps = (1.0 - weight) * eps_phonon + weight * eps_roton

# ── Enforce correct behaviour ──
# At very low k, the phonon branch must dominate
# At high k, beyond the roton minimum, eps rises again
# Cap by the roton parabola to avoid artefacts
eps = np.minimum(eps, eps_roton)
eps = np.maximum(eps, 0.0)

# Ensure phonon region is purely linear for k < 0.6 Å^-1
phonon_mask = k_disp < 0.6
eps[phonon_mask] = hbar * c_s * k_arr[phonon_mask]

# ── Compute critical velocity ──
ratio = eps / (hbar * k_arr)
ratio[ratio < 0] = np.inf
i_min = np.argmin(ratio)
v_c = ratio[i_min]
k_c = k_disp[i_min]
print("=" * 60)
print("  He-4 EXCITATION SPECTRUM: PHONON-ROTON CURVE")
print("=" * 60)
print()

# Print dispersion data for plotting
print("  k (1/Å)     ε/kB (K)      ε/p (m/s)")
print("  " + "-"*45)
# Sample every 40th point
for i in range(0, len(k_arr), 40):
    e_K = eps[i] / kB
    vv  = eps[i] / (hbar * k_arr[i])
    print(f"  {k_disp[i]:8.3f}     {e_K:8.3f}      {vv:10.2f}")

print()
print(f"  Roton gap:     Δ/kB = {Delta/kB:.2f} K")
print(f"  Roton k0:      k0   = {k0*angstrom:.2f} Å⁻¹")
print(f"  Roton mass:    μ    = {mu_rot/m_He:.2f} m_He")
print()
print(f"  ── Landau Critical Velocity ──")
print(f"  v_c = min(ε/p) = {v_c:.2f} m/s")
print(f"  occurs at k    = {k_c:.2f} Å⁻¹")
print(f"  ε(k_c)/kB      = {eps[i_min]/kB:.2f} K")
print()

# Compare with speed of sound
print(f"  Speed of sound:    c_s   = {c_s:.1f} m/s")
print(f"  Landau v_c:        v_c   = {v_c:.1f} m/s")
print(f"  Ratio v_c/c_s:             {v_c/c_s:.3f}")
print()
print("  Note: Experimental critical velocities are typically")
print("  much lower (~1 cm/s to ~1 m/s) due to vortex nucleation.")
print()

# Print the full dispersion for chart rendering
print("=" * 60)
print("  DISPERSION CURVE DATA (for plotting)")
print("=" * 60)
print(f"  {'k(1/A)':>10}  {'E/kB(K)':>10}")
for i in range(0, len(k_arr), 8):
    print(f"  {k_disp[i]:10.4f}  {eps[i]/kB:10.4f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Simulation 2: Superfluid Fraction & Specific Heat (Fortran)

This program computes the superfluid fraction $\rho_s/\rho$ as a function of temperature using the phonon and roton contributions to the normal-fluid density, and evaluates the specific heat near the lambda transition.

Superfluid Fraction ρ_s/ρ and Specific Heat Near Tλ

Fortran

program.f90147 lines

program superfluid_he4
  implicit none

! ── Physical constants ──
  double precision, parameter :: kB     = 1.3806d-23     ! J/K
  double precision, parameter :: hbar   = 1.0546d-34     ! J·s
  double precision, parameter :: pi     = 3.14159265358979d0
  double precision, parameter :: m_He   = 6.6464d-27     ! kg
  double precision, parameter :: angstrom = 1.0d-10      ! m

! ── He-4 parameters ──
  double precision, parameter :: T_lambda = 2.1768d0     ! K
  double precision, parameter :: c_s      = 238.3d0      ! m/s (speed of sound)
  double precision, parameter :: rho_tot  = 145.0d0      ! kg/m^3 (liquid density)
  double precision, parameter :: k0_roton = 1.92d10      ! 1/m (roton wavevector)
  double precision, parameter :: Delta_rot = 8.62d0      ! K (roton gap / kB)
  double precision, parameter :: mu_rot   = 0.16d0 * m_He ! roton effective mass
  double precision, parameter :: nu_exp   = 0.6717d0     ! 3D XY exponent

integer :: i, N
  double precision :: T, dT, rho_n_ph, rho_n_rot, rho_n, f_s
  double precision :: t_red, C_plus, C_minus, alpha_exp, C_val
  double precision :: A_plus, A_minus

! ── Critical exponent for specific heat ──
  alpha_exp = -0.0127d0
  A_plus  = 1.0d0    ! amplitude above T_lambda (normalised)
  A_minus = A_plus / 1.054d0  ! amplitude below T_lambda

! ============================================================
  !  PART 1: Superfluid fraction vs temperature
  ! ============================================================
  write(*,'(A)') '============================================================'
  write(*,'(A)') '  SUPERFLUID FRACTION rho_s/rho vs TEMPERATURE'
  write(*,'(A)') '============================================================'
  write(*,'(A)') ''
  write(*,'(A,A)') '  T(K)        rho_n_ph/rho   rho_n_rot/rho  ', &
                   'rho_s/rho'
  write(*,'(A)') '  ' // repeat('-', 56)

N = 100
  dT = (T_lambda - 0.05d0) / dble(N)

do i = 0, N
    T = 0.05d0 + dble(i) * dT
    if (T >= T_lambda) exit

! Phonon contribution to normal fluid density
    ! rho_n^ph = (2 pi^2 / 45) (kB T)^4 / (hbar^3 c_s^5)
    rho_n_ph = (2.0d0 * pi**2 / 45.0d0) &
               * (kB * T)**4 / (hbar**3 * c_s**5)

! Roton contribution to normal fluid density
    ! rho_n^rot = (2 k0^4 / (3 (2pi)^(3/2))) * sqrt(hbar^2/(mu kB T)) * exp(-Delta/(kB T))
    if (T > 0.3d0) then
      rho_n_rot = (2.0d0 * k0_roton**4) &
                  / (3.0d0 * (2.0d0*pi)**1.5d0) &
                  * sqrt(hbar**2 / (mu_rot * kB * T)) &
                  * exp(-Delta_rot / T)
    else
      rho_n_rot = 0.0d0
    end if

rho_n = rho_n_ph + rho_n_rot

! Near T_lambda, use the power-law form instead
    t_red = (T_lambda - T) / T_lambda
    if (t_red < 0.15d0) then
      ! rho_s/rho ~ (1 - T/T_lambda)^nu
      f_s = t_red**nu_exp
    else
      f_s = max(0.0d0, 1.0d0 - rho_n / rho_tot)
      f_s = min(1.0d0, f_s)
    end if

! Print every 5th point
    if (mod(i, 5) == 0) then
      write(*,'(2X,F6.3,6X,ES12.4,3X,ES12.4,3X,F8.5)') &
        T, rho_n_ph/rho_tot, rho_n_rot/rho_tot, f_s
    end if
  end do

write(*,'(A)') ''
  write(*,'(A)') '  At T = 0 K: rho_s/rho = 1.000 (entire fluid is superfluid)'
  write(*,'(A,F6.3,A)') '  At T = ', T_lambda, ' K: rho_s/rho = 0 (lambda transition)'
  write(*,'(A)') ''

! ============================================================
  !  PART 2: Specific heat near the lambda transition
  ! ============================================================
  write(*,'(A)') '============================================================'
  write(*,'(A)') '  SPECIFIC HEAT NEAR THE LAMBDA TRANSITION'
  write(*,'(A)') '============================================================'
  write(*,'(A)') ''
  write(*,'(A,F10.6)') '  Critical exponent alpha = ', alpha_exp
  write(*,'(A,F10.4)') '  T_lambda = ', T_lambda
  write(*,'(A,F10.4)') '  Amplitude ratio A+/A-  = ', A_plus/A_minus
  write(*,'(A)') ''
  write(*,'(A)') '  T(K)        |t|           C/A       (side)'
  write(*,'(A)') '  ' // repeat('-', 50)

! Below T_lambda
  do i = 1, 12
    t_red = 10.0d0**(-0.5d0 * dble(i))
    T = T_lambda * (1.0d0 - t_red)
    C_val = A_minus * abs(t_red)**(-alpha_exp)
    write(*,'(2X,F8.5,4X,ES10.3,4X,F10.4,4X,A)') &
      T, t_red, C_val, 'below'
  end do

write(*,'(2X,A)') '  --- lambda point ---'

! Above T_lambda
  do i = 12, 1, -1
    t_red = 10.0d0**(-0.5d0 * dble(i))
    T = T_lambda * (1.0d0 + t_red)
    C_val = A_plus * abs(t_red)**(-alpha_exp)
    write(*,'(2X,F8.5,4X,ES10.3,4X,F10.4,4X,A)') &
      T, t_red, C_val, 'above'
  end do

write(*,'(A)') ''
  write(*,'(A)') '  Note: C(T) shows a characteristic lambda-shaped cusp.'
  write(*,'(A)') '  Since alpha < 0, C does not truly diverge but has a'
  write(*,'(A)') '  finite cusp at T_lambda.'
  write(*,'(A)') ''

! ============================================================
  !  PART 3: Summary of key results
  ! ============================================================
  write(*,'(A)') '============================================================'
  write(*,'(A)') '  SUMMARY'
  write(*,'(A)') '============================================================'
  write(*,'(A)') ''
  write(*,'(A)') '  Superfluid He-4 Key Properties:'
  write(*,'(A,F8.4,A)') '    Lambda temperature:   ', T_lambda, ' K'
  write(*,'(A,F8.1,A)') '    Speed of sound:       ', c_s, ' m/s'
  write(*,'(A,F8.2,A)') '    Roton gap Delta/kB:   ', Delta_rot, ' K'
  write(*,'(A,F8.2,A)') '    Roton wavevector k0:  ', k0_roton*angstrom, ' 1/Å'
  write(*,'(A,F8.4)')   '    XY exponent nu:       ', nu_exp
  write(*,'(A,F8.4)')   '    Specific heat alpha:  ', alpha_exp
  write(*,'(A)') ''
  write(*,'(A)') '============================================================'

end program superfluid_he4

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← Chapter 1: BEC Chapter 3: Two-Fluid Model →