Part II — Chapter 3

Two-Fluid Model

The Tisza-Landau two-fluid model describes superfluid helium as an intimate interpenetration of two fluid components — the superfluid and the normal fluid — each carrying its own velocity field. This framework accounts for the most striking transport phenomena of He II, including second sound, the fountain effect, and thermal counterflow.

The Tisza-Landau Two-Fluid Model

Historical Context

Tisza (1938) proposed that below $T_\lambda$, liquid $^4$He behaves as if composed of two interpenetrating fluids. Landau (1941) placed this on a firmer theoretical footing by identifying the normal component with the gas of elementary excitations (phonons and rotons) and the superfluid component with the condensate and its coherent ground-state motion.

The key postulate is that the total density $\rho$ decomposes as:

$$\rho = \rho_s + \rho_n$$

where $\rho_s$ is the superfluid density and $\rho_n$ is the normal density. Each component carries its own velocity field, $\mathbf{v}_s$and $\mathbf{v}_n$, with the total momentum density given by:

$$\mathbf{j} = \rho_s \mathbf{v}_s + \rho_n \mathbf{v}_n$$

Properties of the Two Components

The superfluid component carries zero entropy ($s_s = 0$), flows irrotationally ($\nabla \times \mathbf{v}_s = 0$), has zero viscosity, and dominates at $T \to 0$. Its velocity is set by the condensate phase: $\mathbf{v}_s = (\hbar/m)\nabla\phi$.

The normal component carries all the entropy, has finite viscosity $\eta_n$, consists of thermal excitations (phonons and rotons), and dominates at $T \to T_\lambda$.

Temperature Dependence of Density Fractions

Near the lambda transition, the superfluid fraction vanishes with a characteristic power law:

$$\frac{\rho_s}{\rho} \propto \left(1 - \frac{T}{T_\lambda}\right)^{\zeta}, \qquad \zeta \approx 0.6717$$

where the exponent $\zeta$ is related to the correlation-length critical exponent$\nu$ of the XY universality class. At $T = 0$, the system is entirely superfluid ($\rho_s/\rho = 1$), while at $T = T_\lambda \approx 2.177$ K the superfluid fraction drops to zero continuously.

Two-Fluid Equations of Motion

Conservation Laws

The linearized two-fluid equations (neglecting viscosity and second-order terms) consist of four coupled equations:

Mass conservation:

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho_s \mathbf{v}_s + \rho_n \mathbf{v}_n) = 0$$

Entropy conservation (reversible flow):

$$\frac{\partial (\rho s)}{\partial t} + \rho s \, \nabla \cdot \mathbf{v}_n = 0$$

Superfluid acceleration (London equation):

$$\frac{\partial \mathbf{v}_s}{\partial t} = -\nabla\!\left(\mu + \tfrac{1}{2}v_s^2\right) = -\frac{1}{\rho}\nabla P + s\nabla T$$

Momentum conservation:

$$\frac{\partial \mathbf{j}}{\partial t} = -\nabla P$$

Linearized Wave Equations

For small-amplitude oscillations, the secular equation for wave speed $u$ is:

$$u^4 - \left(c_1^2 + \frac{T s^2 \rho_s}{C_V \rho_n}\right)u^2 + c_1^2 \frac{Ts^2\rho_s}{C_P \rho_n} = 0$$

The two roots yield first sound ($u \approx c_1$) and second sound ($u \approx c_2$) as nearly independent modes.

First Sound: Density Waves

First sound is the ordinary pressure-density wave in which normal and superfluid components oscillate in phase. Both $\mathbf{v}_s$and $\mathbf{v}_n$ are parallel, and the wave propagates with velocity:

$$c_1^2 = \left(\frac{\partial P}{\partial \rho}\right)_s$$

In liquid $^4$He, $c_1 \approx 238$ m/s near $T = 0$, varying only weakly with temperature. This is essentially the adiabatic compressional sound speed, analogous to ordinary sound in a classical fluid.

During a first-sound oscillation, local temperature and entropy fluctuations are small because both components move together — there is no relative motion to transport entropy preferentially.

Second Sound: Temperature Waves

Second sound is a uniquely superfluid phenomenon: a propagating wave of temperature(or equivalently entropy) with nearly zero pressure variation. The normal and superfluid components oscillate in antiphase such that the total mass current$\mathbf{j} \approx 0$, while entropy is transported back and forth by the normal component:

$$c_2^2 = \frac{T s^2 \rho_s}{C_V \rho_n}$$

The second sound velocity vanishes at both ends of the temperature range:

As $T \to T_\lambda$: $\rho_s \to 0$, so $c_2 \to 0$
As $T \to 0$: $\rho_n \to 0$ but the entropy $s \to 0$ faster, so $c_2 \to 0$ as well

The maximum value of $c_2 \approx 20$ m/s occurs near $T \approx 1.65$ K. At very low temperatures where only phonons exist, $c_2 \to c_1/\sqrt{3}$, consistent with a phonon gas having $C_V = 3P/T$.

Experimental Detection

Peshkov (1944) first measured second sound by using a heater as a “loudspeaker” and a sensitive thermometer as a “microphone.” A periodic heat input launches temperature waves that propagate at $c_2$ through the superfluid, detected as temperature oscillations at a distant point. This experiment provided dramatic confirmation of the two-fluid model.

Fourth Sound in Narrow Channels

When superfluid helium is confined in channels narrow enough that the normal component is clamped by viscous forces ($\mathbf{v}_n = 0$), only the superfluid component can oscillate. The resulting wave mode is called fourth sound, with velocity:

$$c_4^2 = \frac{\rho_s}{\rho}\, c_1^2 + \frac{T s^2 \rho_s}{C_V \rho}$$

In practice, the second term is small, so approximately:

$$c_4 \approx \sqrt{\frac{\rho_s}{\rho}}\; c_1$$

Fourth sound has been measured in superleaks (fine powders, porous media, and narrow capillaries). At $T = 0$, $c_4 = c_1$ since $\rho_s = \rho$. Near $T_\lambda$, $c_4 \to 0$ as the superfluid fraction vanishes. There is also a third sound — surface waves on thin superfluid films — governed by a similar expression modified by van der Waals restoring forces.

Mutual Friction Between Components

In the idealized two-fluid model, the two components flow independently. In reality, interaction between quantized vortex lines and the normal fluid produces a mutual friction force $\mathbf{F}_{ns}$that couples the two velocity fields.

The Hall-Vinen mutual friction force per unit volume takes the form:

$$\mathbf{F}_{ns} = -\frac{B \rho_s \rho_n}{2\rho}\,\hat{\boldsymbol{\omega}} \times \bigl[\boldsymbol{\omega}_s \times (\mathbf{v}_n - \mathbf{v}_s)\bigr] - \frac{B' \rho_s \rho_n}{2\rho}\,\boldsymbol{\omega}_s \times (\mathbf{v}_n - \mathbf{v}_s)$$

where $\boldsymbol{\omega}_s = \nabla \times \mathbf{v}_s$ (nonzero only at vortex cores), and $B$, $B'$ are dimensionless mutual friction coefficients measured by second-sound attenuation in rotating helium. The $B$ term is dissipative and dominates the attenuation; the $B'$ term is reactive.

In rotating He II, this friction attenuates second sound proportionally to the vortex line density $L = 2\Omega/\kappa$, providing a standard probe of the vortex state in both rotating and turbulent superfluids.

Thermal Counterflow and Its Breakdown

The Counterflow Mechanism

When heat $\dot{Q}$ is applied to superfluid helium in a channel closed at one end, the steady state requires zero net mass flow: $\rho_s \mathbf{v}_s + \rho_n \mathbf{v}_n = 0$. The normal component flows away from the heater carrying entropy, while the superfluid flows toward it to maintain mass balance:

$$\dot{Q} = \rho s T \, v_n A, \qquad v_{ns} = v_n - v_s = \frac{\dot{Q}}{\rho_s s T A}\,\frac{\rho}{\rho_n}$$

where $A$ is the channel cross-section. This internal convection transports heat without net mass flow, giving He II its enormous effective thermal conductivity.

Critical Velocity and Breakdown

Above a critical counterflow velocity $v_{ns}^{\rm crit}$, quantized vortex lines nucleate and form a self-sustaining tangle. Three regimes are observed: TI (laminar, below $v_{ns}^{c1}$), TII (vortex tangle with $L \propto v_{ns}^2$), and TIII (coupled superfluid-classical turbulence at higher velocities).

Vinen Equation for Tangle Density

The evolution of the vortex line density $L$ in counterflow turbulence is described by the Vinen equation:

$$\frac{dL}{dt} = \alpha_V v_{ns} L^{3/2} - \beta_V \frac{\kappa}{2\pi} L^2$$

The first term represents vortex growth driven by the counterflow velocity, and the second represents vortex annihilation through reconnections. In steady state,$L_{\rm ss} = (\alpha_V / \beta_V)^2 (2\pi/\kappa)^2 v_{ns}^2 = \gamma^2 v_{ns}^2$, confirming the observed $L \propto v_{ns}^2$ scaling.

Detailed Derivation: Second Sound Velocity

Step 1: Linearized Two-Fluid Equations

We start from the four linearized two-fluid equations for small perturbations$\delta\rho$, $\delta s$, $\delta\mathbf{v}_s$, and $\delta\mathbf{v}_n$about equilibrium. Assuming plane-wave solutions $\propto e^{i(\mathbf{q}\cdot\mathbf{r} - \omega t)}$and restricting to the case where the total mass current vanishes ($\rho_s \delta\mathbf{v}_s + \rho_n\delta\mathbf{v}_n = 0$):

$$\rho_s\,\delta\mathbf{v}_s = -\rho_n\,\delta\mathbf{v}_n$$

This condition means no net mass flow — only entropy is transported.

Step 2: Entropy Oscillation Equation

From entropy conservation, $\partial(\rho s)/\partial t + \rho s\,\nabla\cdot\mathbf{v}_n = 0$, the entropy fluctuation satisfies:

$$-i\omega\,\rho\,\delta s + \rho s\,(i\mathbf{q}\cdot\delta\mathbf{v}_n) = 0 \quad \Longrightarrow \quad \delta s = \frac{s\,\mathbf{q}\cdot\delta\mathbf{v}_n}{\omega}$$

From the superfluid acceleration equation, $\partial\mathbf{v}_s/\partial t = -(1/\rho)\nabla P + s\nabla T$, and using the thermodynamic identity $\delta P = (\partial P/\partial T)_\rho\,\delta T + (\partial P/\partial\rho)_T\,\delta\rho$, with $\delta\rho = 0$ for a second-sound mode (no density change):

$$-i\omega\,\delta\mathbf{v}_s = -\frac{1}{\rho}\left(\frac{\partial P}{\partial T}\right)_\rho i\mathbf{q}\,\delta T + s\,i\mathbf{q}\,\delta T$$

Step 3: Combining to Obtain the Dispersion Relation

Using the thermodynamic relation $(\partial P/\partial T)_\rho = \rho s$ and combining the entropy equation with the superfluid acceleration equation, the temperature fluctuation satisfies a wave equation with dispersion:

$$\omega^2 = q^2 \cdot \frac{T s^2 \rho_s}{C_V \rho_n}$$

Therefore the second sound velocity is:

$$\boxed{c_2 = \sqrt{\frac{T s^2 \rho_s}{C_V \rho_n}}}$$

At very low temperatures (phonon-dominated regime), $C_V = \rho s T/3$ for a phonon gas, and$\rho_s \approx \rho$, so $c_2^2 = Ts^2\rho/(C_V\rho_n) = 3s^2\rho^2/(s\rho_n) \to c_1^2/3$, giving $c_2 \to c_1/\sqrt{3}$, the well-known low-temperature limit.

Derivation: The Coupled Sound Mode Secular Equation

Step 1: General Plane-Wave Ansatz

For the general case (allowing both density and temperature fluctuations), assume all perturbations vary as $e^{i(qx - \omega t)}$ with phase velocity $u = \omega/q$. The mass conservation equation gives:

$$-\omega\,\delta\rho + q(\rho_s\,\delta v_s + \rho_n\,\delta v_n) = 0$$

The momentum equation gives: $-\omega(\rho_s\,\delta v_s + \rho_n\,\delta v_n) + q\,\delta P = 0$.

Step 2: Eliminating Variables

Expressing $\delta P$ in terms of $\delta\rho$ and $\delta T$ via thermodynamics, and combining all four equations, one obtains a $2\times 2$ eigenvalue problem for$(\delta\rho, \delta T)$. The secular equation is:

$$u^4 - u^2\left[c_1^2 + \frac{Ts^2\rho_s}{C_V\rho_n}\right] + c_1^2\,\frac{Ts^2\rho_s}{C_P\rho_n} = 0$$

This is a quadratic in $u^2$ with two roots. Since $C_P \approx C_V$ in a nearly incompressible liquid (the difference $C_P - C_V = TV\alpha_P^2/\kappa_T$ is very small for He-II), the coupling between the two modes is weak. The two approximate roots are:

$$u_1^2 \approx c_1^2, \qquad u_2^2 \approx \frac{Ts^2\rho_s}{C_V\rho_n}$$

First sound is primarily a density-pressure wave with $\delta T \approx 0$, and second sound is primarily a temperature-entropy wave with $\delta P \approx 0$.

Derivation: Thermomechanical (Fountain) Pressure

Step 1: Chemical Potential Equilibrium

In equilibrium, the chemical potential must be uniform: $\nabla\mu = 0$. Using the Gibbs-Duhem relation for the superfluid component, $d\mu = -s\,dT + (1/\rho)\,dP$, we obtain:

$$\nabla\mu = -s\,\nabla T + \frac{1}{\rho}\,\nabla P = 0$$

Step 2: The London-Tisza Relation

Rearranging immediately gives the London-Tisza relation:

$$\nabla P = \rho\,s\,\nabla T$$

For a finite temperature difference $\Delta T$ across a superleak:

$$\Delta P = \rho\,s\,\Delta T$$

Near $T = 1.5$ K, $s \approx 750$ J/(kg$\cdot$K) and $\rho = 145$ kg/m$^3$, so $\Delta P/\Delta T \approx 1.09 \times 10^5$ Pa/K. A temperature difference of just 10 mK generates a pressure of about 1000 Pa (enough to drive a fountain several centimeters high).

Historical Context

The two-fluid model has two independent origins. Laszlo Tiszaproposed it in 1938, inspired by Fritz London's suggestion that the lambda transition was related to Bose-Einstein condensation. Tisza imagined the superfluid component as a “condensate” carrying no entropy, and predicted the existence of second sound — a bold prediction later confirmed experimentally.

Lev Landau (1941) independently developed the model from a different perspective, identifying the normal fluid with a gas of elementary excitations (phonons and rotons) rather than with “uncondensed atoms.” Landau's formulation was more rigorous and led to the correct hydrodynamic equations, the critical velocity criterion, and predictions for the excitation spectrum. Landau initially disputed Tisza's priority and rejected the BEC connection, leading to a prolonged controversy. Modern understanding recognizes both contributions as essential.

Vasily Peshkov first detected second sound in 1944 in Kapitsa's laboratory in Moscow, using a pulsed heater and bolometric detection. The measured velocity of about 19 m/s at 1.4 K confirmed the two-fluid prediction. The development of the Hall-Vinen mutual friction formalism in the 1950s extended the model to account for the coupling between the two fluids through quantized vortex lines, completing the hydrodynamic description of rotating and turbulent superfluids.

Applications of Two-Fluid Hydrodynamics

Cryogenic heat transfer: The two-fluid counterflow mechanism gives He-II an effective thermal conductivity $\sim 10^5$ times that of copper. This “thermal superconductivity” is exploited in cooling the superconducting magnets of the LHC and in cryogenic engineering for uniform temperature baths.
Second sound thermometry and imaging: Second sound can be used to detect and image quantized vortex lines, since vortices attenuate second sound via mutual friction. This technique enables non-invasive mapping of vortex tangles in quantum turbulence experiments.
Superfluid fountain pumps: The thermomechanical effect provides a pressure pump with no moving parts. Fountain-effect pumps are used in space-based helium dewars (e.g., the Gravity Probe B mission) where mechanical pumps would introduce vibration.
Analogue to relativistic field theories: The two-fluid equations are mathematically analogous to certain aspects of relativistic two-component field theories. The phonon-roton spectrum of He-II provides a condensed-matter analogue of the Higgs mechanism, and acoustic black holes can be realized using converging superfluid flows.
Dilution refrigeration: The two-fluid concept extends to$^3$He-$^4$He mixtures in dilution refrigerators, where the “superfluid” is $^4$He-rich and the “normal fluid” includes$^3$He quasiparticles. This enables cooling to millikelvin temperatures for quantum computing and low-temperature experiments.

Python Simulation: Sound Velocities & Density Fractions

First and second sound velocities alongside superfluid and normal density fractions as functions of temperature in He-4, computed from the two-fluid model.

Two-Fluid Model: Sound Velocities in He-4

Python

Temperature dependence of first sound, second sound, and density fractions

script.py110 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Two-Fluid Model: Sound Velocities & Density Fractions in He-4

T_lambda = 2.1768  # K, lambda transition temperature

# Temperature array (avoid T=0 and T=T_lambda exactly)
T = np.linspace(0.05, T_lambda - 0.01, 500)
t = T / T_lambda  # reduced temperature

# Density fractions (empirical power law)
rho_n_frac = t**5.6
rho_s_frac = 1.0 - rho_n_frac

# Specific heat model (phonon T^3 at low T, rising near T_lambda)
C_V = np.where(T < 0.6,
    0.0205 * T**3,
    0.0205 * 0.6**3 + 5.6 * (T - 0.6)**1.5)

# Entropy via cumulative integration of C_V/T
s = np.zeros_like(T)
dT = T[1] - T[0]
for i in range(1, len(T)):
    s[i] = s[i-1] + 0.5 * (C_V[i]/T[i] + C_V[i-1]/T[i-1]) * dT

# First sound velocity (empirical temperature dependence)
c1 = 238.3 * (1.0 - 0.035 * t**2 - 0.06 * t**6)

# Second sound: c2^2 = T s^2 rho_s / (C_V rho_n)
rho_n_safe = np.maximum(rho_n_frac, 1e-8)
C_V_safe = np.maximum(C_V, 1e-10)

c2_sq = T * s**2 * rho_s_frac / (C_V_safe * rho_n_safe)
c2_sq = np.maximum(c2_sq, 0)
c2 = np.sqrt(c2_sq)

# Scale to match experimental peak ~20 m/s
c2 = c2 * 20.0 / max(np.max(c2), 1e-10)

# ── Plotting ──
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.patch.set_facecolor('#0a0a1a')
fig.suptitle('Two-Fluid Model of Superfluid He-4',
             color='white', fontsize=16, fontweight='bold', y=0.98)

for ax in axes.flat:
    ax.set_facecolor('#0a0a1a')
    for spine in ax.spines.values():
        spine.set_edgecolor('#1e3a5f')

# Panel 1: Density fractions
ax = axes[0, 0]
ax.plot(T, rho_s_frac, color='#00e5ff', linewidth=2.5, label=r'$\rho_s/\rho$')
ax.plot(T, rho_n_frac, color='#ff6b6b', linewidth=2.5, label=r'$\rho_n/\rho$')
ax.axvline(T_lambda, color='#ffd700', linestyle='--', alpha=0.6, label=r'$T_\lambda$')
ax.set_xlabel('Temperature [K]', color='#94a3b8', fontsize=11)
ax.set_ylabel('Density Fraction', color='#94a3b8', fontsize=11)
ax.set_title('Superfluid & Normal Density Fractions', color='#00e5ff', fontsize=12)
ax.legend(fontsize=10, facecolor='#0d1117', edgecolor='#1e3a5f', labelcolor='white')
ax.tick_params(colors='#94a3b8')
ax.set_xlim(0, T_lambda)
ax.set_ylim(0, 1.05)

# Panel 2: First sound velocity
ax = axes[0, 1]
ax.plot(T, c1, color='#7c3aed', linewidth=2.5, label='First sound $c_1$')
ax.set_xlabel('Temperature [K]', color='#94a3b8', fontsize=11)
ax.set_ylabel('Velocity [m/s]', color='#94a3b8', fontsize=11)
ax.set_title('First Sound (Pressure/Density Waves)', color='#a78bfa', fontsize=12)
ax.legend(fontsize=10, facecolor='#0d1117', edgecolor='#1e3a5f', labelcolor='white')
ax.tick_params(colors='#94a3b8')
ax.set_xlim(0, T_lambda)
ax.set_ylim(220, 245)

# Panel 3: Second sound velocity
ax = axes[1, 0]
ax.plot(T, c2, color='#f59e0b', linewidth=2.5, label='Second sound $c_2$')
ax.set_xlabel('Temperature [K]', color='#94a3b8', fontsize=11)
ax.set_ylabel('Velocity [m/s]', color='#94a3b8', fontsize=11)
ax.set_title('Second Sound (Temperature/Entropy Waves)', color='#fbbf24', fontsize=12)
ax.legend(fontsize=10, facecolor='#0d1117', edgecolor='#1e3a5f', labelcolor='white')
ax.tick_params(colors='#94a3b8')
ax.set_xlim(0, T_lambda)
ax.set_ylim(0, 25)

# Panel 4: Fourth sound velocity
c4 = np.sqrt(rho_s_frac) * c1
ax = axes[1, 1]
ax.plot(T, c4, color='#22c55e', linewidth=2.5, label='Fourth sound $c_4$')
ax.plot(T, c1, color='#7c3aed', linewidth=1.5, alpha=0.5, linestyle='--', label='$c_1$ (reference)')
ax.set_xlabel('Temperature [K]', color='#94a3b8', fontsize=11)
ax.set_ylabel('Velocity [m/s]', color='#94a3b8', fontsize=11)
ax.set_title('Fourth Sound (Narrow Channels)', color='#4ade80', fontsize=12)
ax.legend(fontsize=10, facecolor='#0d1117', edgecolor='#1e3a5f', labelcolor='white')
ax.tick_params(colors='#94a3b8')
ax.set_xlim(0, T_lambda)

plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.savefig('two_fluid_model.png', dpi=130, bbox_inches='tight',
            facecolor=fig.get_facecolor())

print(f"Lambda temperature: T_lambda = {T_lambda} K")
print(f"First sound velocity at T=0: c1(0) ~ {c1[0]:.1f} m/s")
print(f"Peak second sound velocity: c2_max ~ {np.max(c2):.1f} m/s")
print(f"At T=1.0 K: rho_s/rho={np.interp(1.0,T,rho_s_frac):.4f}, rho_n/rho={np.interp(1.0,T,rho_n_frac):.4f}")
print(f"At T=2.0 K: rho_s/rho={np.interp(2.0,T,rho_s_frac):.4f}, rho_n/rho={np.interp(2.0,T,rho_n_frac):.4f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran Implementation: Wave Mode Solver

Numerical solution of the two-fluid wave equations to extract first sound, second sound, and fourth sound velocities across the superfluid temperature range.

Two-Fluid Wave Equation Eigenvalue Solver

Fortran

Tabular computation of sound mode velocities from linearized two-fluid hydrodynamics

program.f9083 lines

program two_fluid_wave_modes
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  integer, parameter :: N_T = 40
  integer :: i
  real(dp) :: T_lambda, T, t_red, dT
  real(dp) :: rho_s, rho_n, rho_total
  real(dp) :: c1_sq, c2_sq, c1, c2
  real(dp) :: s_ent, C_v, P_deriv

! Two-Fluid Wave Equations: Sound Mode Solver
  T_lambda = 2.1768_dp
  rho_total = 145.0_dp   ! kg/m^3

write(*,'(A)') '  Two-Fluid Model: Sound Velocity Computation'
  write(*,'(A,F8.4,A)') '  T_lambda = ', T_lambda, ' K'
  write(*,'(A,F8.1,A)') '  rho      = ', rho_total, ' kg/m^3'
  write(*,'(A)') ''
  write(*,'(A)') '---------------------------------------------------'
  write(*,'(A6,A10,A10,A10,A10,A10)') &
       'T[K]', 'rho_s/rho', 'rho_n/rho', 'c1[m/s]', 'c2[m/s]', 'c2/c1'
  write(*,'(A)') '---------------------------------------------------'

dT = (T_lambda - 0.1_dp) / real(N_T - 1, dp)

do i = 1, N_T
     T = 0.1_dp + real(i - 1, dp) * dT
     t_red = T / T_lambda

! Superfluid fraction (empirical power law)
     rho_n = rho_total * t_red**5.6_dp
     rho_s = rho_total - rho_n

! Specific heat and entropy (simplified models)
     if (T < 0.6_dp) then
        C_v = 20.5_dp * T**3
     else
        C_v = 20.5_dp * 0.6_dp**3 + 5600.0_dp * (T - 0.6_dp)**1.5_dp
     end if
     C_v = max(C_v, 1.0e-6_dp)

if (T < 0.6_dp) then
        s_ent = 6.833_dp * T**3
     else
        s_ent = 6.833_dp * 0.6_dp**3 + 3733.0_dp * (T - 0.6_dp)**1.5_dp / 1.5_dp
     end if
     s_ent = max(s_ent, 1.0e-8_dp)

! First sound: c1^2 = (dP/drho)_s
     P_deriv = 238.3_dp**2 * (1.0_dp - 0.035_dp*t_red**2 - 0.06_dp*t_red**6)
     c1_sq = P_deriv
     c1 = sqrt(max(c1_sq, 0.0_dp))

! Second sound: c2^2 = T s^2 rho_s / (C_v rho_n)
     if (rho_n > 1.0e-10_dp) then
        c2_sq = T * s_ent**2 * rho_s / (C_v * rho_n)
     else
        c2_sq = 0.0_dp
     end if

c2 = sqrt(max(c2_sq, 0.0_dp))
     c2 = min(c2, 25.0_dp)

write(*,'(F6.3, F10.4, F10.4, F10.2, F10.2, F10.4)') &
          T, rho_s/rho_total, rho_n/rho_total, c1, c2, c2/c1
  end do

write(*,'(A)') '---------------------------------------------------'
  write(*,'(A)') ''
  write(*,'(A)') '  First sound: in-phase oscillation, c1 ~ 238 m/s'
  write(*,'(A)') '  Second sound: counter-phase, c2 ~ 20 m/s peak'
  write(*,'(A)') ''

! ── Fourth sound (narrow channels: v_n = 0) ──
  write(*,'(A)') ''
  write(*,'(A)') '  Fourth sound: c4 = sqrt(rho_s/rho) * c1'
  write(*,'(A)') '  At T=0.5 K:'
  t_red = 0.5_dp / T_lambda
  rho_s = rho_total * (1.0_dp - t_red**5.6_dp)
  c1 = 238.3_dp * sqrt(1.0_dp - 0.035_dp*t_red**2 - 0.06_dp*t_red**6)
  write(*,'(A,F8.2,A)') '    c4 = ', sqrt(rho_s/rho_total)*c1, ' m/s'

end program two_fluid_wave_modes

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Summary of Sound Modes

Mode	Character	Speed
First sound	Pressure wave, $c_1^2 = (\partial P/\partial\rho)_s$	$\sim 238$ m/s
Second sound	Entropy wave, $c_2^2 = Ts^2\rho_s/(C_V\rho_n)$	$\sim 20$ m/s
Fourth sound	Superleak, $c_4^2 \approx (\rho_s/\rho)\, c_1^2$	$\sim 200$ m/s

← Chapter 2: Superfluid Helium Chapter 4: Gross-Pitaevskii →