Module 4 · Week 7 · Graduate

Out-of-Equilibrium Physics: Boltzmann, Decoupling, Phase Transitions

The third Sakharov condition — departure from thermal equilibrium — is the hardest to satisfy and the most model-dependent. This module covers the machinery: equilibrium distributions, the Boltzmann equation, decoupling criteria, and first-order phase transitions with bubble nucleation.

1. Thermal Equilibrium & Detailed Balance

Equilibrium phase-space distributions for bosons (+) and fermions (−) are

\[f_{eq}(E,T,\mu) = \frac{1}{e^{(E-\mu)/T} \mp 1}.\]

Number and energy densities for a single-species gas:

\[n_{eq} = \frac{g}{2\pi^2}\int_m^\infty \frac{\sqrt{E^2-m^2}\,E}{e^{(E-\mu)/T}\mp 1}\,dE.\]

Two important limits:

Relativistic (T >> m): \(n_{eq} = \frac{\zeta(3)}{\pi^2} g T^3\) (fermions with factor 3/4).
Non-relativistic (T << m): \(n_{eq} = g \left(\frac{mT}{2\pi}\right)^{3/2} e^{-m/T}\) (Boltzmann-suppressed).

In equilibrium, detailed balance forces the net baryon charge to zero (by the CPT argument of Module 1). So baryogenesis requires a dynamical trigger that pushes some species out of equilibrium.

2. The Boltzmann Equation in an Expanding Universe

The covariant Liouville operator in FLRW spacetime acting on a homogeneous, isotropic distribution gives

\[\frac{df}{dt} = \frac{\partial f}{\partial t} - H\, p\,\frac{\partial f}{\partial p} = C[f]\]

where \(C[f]\) is the collision integral. Integrating over momentum gives the number-density Boltzmann equation (Kolb-Turner form):

\[\frac{dn}{dt} + 3H\,n = -\langle \sigma v\rangle\bigl(n^2 - n_{eq}^2\bigr).\]

The 3Hn term is dilution from expansion; the right-hand side is the interaction rate that drives n toward n_eq. Switching to the yield \(Y = n/s\) (which removes expansion for entropy-conserving processes):

\[\frac{dY}{dx} = -\frac{\langle\sigma v\rangle s}{xH}\,\bigl(Y^2 - Y_{eq}^2\bigr), \quad x\equiv m/T.\]

3. The Decoupling Criterion \(\Gamma < H\)

A process with rate \(\Gamma\) remains in equilibrium so long as it fires faster than the universe expands:

\[\text{equilibrium} \iff \Gamma \gtrsim H.\]

When \(\Gamma \lesssim H\), the species “decouples” or “freezes out” — its comoving number density stops evolving. For heavy particles of mass M whose decay rate is \(\Gamma_X\), the ratio at \(T = M\) is (with \(H \sim T^2/M_{Pl}\))

\[K \equiv \frac{\Gamma_X}{H(T=M)} = \frac{\Gamma_X\, M_{Pl}}{1.66\sqrt{g_*}\, M^2}.\]

\(K\) is the “decay parameter”. \(K\ll 1\): strong departure from equilibrium (ideal for baryogenesis); \(K\gg 1\): X stays in equilibrium and any asymmetry is washed out. Efficient baryogenesis requires \(K\sim 1\)– \(10\).

4. First-Order Phase Transitions & Bubble Nucleation

An alternative way to go out of equilibrium is a first-order cosmological phase transition. The free-energy potential develops two minima separated by a barrier; at the critical temperature they are degenerate, but below T_c the true vacuum is lower. Bubbles of true vacuum nucleate and expand through the false-vacuum plasma.

Bubble nucleation rate

By the thin-wall approximation, the nucleation rate per unit volume is

\[\Gamma/V \sim T^4\, e^{-S_3/T}, \qquad S_3 = \int d^3x\left[\frac{1}{2}(\nabla\phi)^2 + V_T(\phi)\right].\]

The critical bubble has O(3)-symmetric profile \(\phi(r)\) satisfying \(\phi'' + (2/r)\phi' = dV_T/d\phi\). When the integrated nucleation probability per Hubble volume reaches O(1), bubbles percolate and the transition completes.

Why bubbles matter for EW baryogenesis

Bubble walls are out of equilibrium: they sweep through the plasma at speed \(v_w\), creating sharp gradients in the Higgs field. CP-violating reflection at the wall produces chiral asymmetries in the symmetric (outside) phase, which sphalerons — still active in the symmetric phase — convert to baryon number. The baryon asymmetry is then protected inside the bubble, where sphalerons are exponentially suppressed.

5. Diagram: Decoupling and Bubble Nucleation

Two mechanisms for the Sakharov out-of-equilibrium condition: (A) decoupling of a heavy species as \(\Gamma\) falls below \(H\), and (B) bubble nucleation during a first-order phase transition.

6. Application: Right-Handed Neutrino Decays

The archetypal baryogenesis Boltzmann system tracks the heavy right-handed neutrino \(N\) and the generated lepton asymmetry \(N_{B-L}\). Using the dimensionless variable \(z = M_N/T\):

\[\frac{dN_N}{dz} = -D\,(N_N - N_N^{eq}),\qquad \frac{dN_{B-L}}{dz} = \varepsilon\, D\,(N_N - N_N^{eq}) - W\, N_{B-L}.\]

Here D is the decay term, W is the washout term, and \(\varepsilon\) is the CP asymmetry. This is the equation we solve numerically in the simulation below and again in Module 7 with more detail.

Simulation: Boltzmann Equation & Bubble Expansion

The simulation below (a) numerically integrates the standard freeze-out Boltzmann equation for various decay parameters K, (b) solves the leptogenesis-style coupled system showing the B−L asymmetry build-up and washout, and (c) draws a bubble profile \(\phi(r)\)from a numerical Coleman-de Luccia calculation in a simple scalar potential.

Python

script.py209 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from scipy.integrate import odeint

plt.rcParams['font.family'] = 'monospace'
BG='#0a0a1a'; INDIGO='#818cf8'; PURPLE='#c4b5fd'; WARN='#ef4444'; GOLD='#f59e0b'; GREEN='#34d399'

fig = plt.figure(figsize=(14, 12), facecolor=BG)
fig.suptitle('Out-of-Equilibrium: Boltzmann evolution and phase-transition bubbles',
             color='white', fontsize=16, fontweight='bold', y=0.98)
gs = gridspec.GridSpec(3, 2, figure=fig, hspace=0.45, wspace=0.35)

# ============================================================
# Plot 1: Freeze-out Boltzmann for a heavy decaying particle
# dY/dx = -K/x^2 * (Y - Y_eq),  x = M/T
# Equilibrium yield for a non-rel species: Y_eq ~ x^(3/2) exp(-x)
# ============================================================
def Y_eq(x):
    return 0.145 * x**1.5 * np.exp(-x)  # normalized roughly

def boltzmann_rhs(Y, x, K):
    yeq = Y_eq(x)
    return -K/x**2 * (Y - yeq)

x_arr = np.linspace(0.1, 30, 2000)
Y_eq_arr = Y_eq(x_arr)

ax1 = fig.add_subplot(gs[0, :], facecolor=BG)
ax1.semilogy(x_arr, Y_eq_arr, color='white', lw=1.5, ls='--', label='Y_eq (equilibrium)')

K_values = [0.01, 1.0, 10, 100]
colors = [GREEN, PURPLE, GOLD, WARN]
for K, col in zip(K_values, colors):
    Y0 = Y_eq(x_arr[0])
    Y_sol = odeint(boltzmann_rhs, Y0, x_arr, args=(K,))
    ax1.semilogy(x_arr, Y_sol[:,0], lw=2.2, color=col, label=f'K = {K}')

ax1.set_xlabel('x = M / T', color=PURPLE, fontsize=10)
ax1.set_ylabel('Y = n / s', color=PURPLE, fontsize=10)
ax1.set_title('Freeze-out of a heavy particle for various decay parameters K',
              color='white', fontsize=12)
ax1.legend(fontsize=10, facecolor='#0f172a', labelcolor='white', edgecolor='#312e81')
ax1.set_ylim(1e-8, 1)
ax1.tick_params(colors=PURPLE, labelsize=9)
ax1.grid(True, alpha=0.15, color=PURPLE)
for sp in ax1.spines.values(): sp.set_color('#312e81')

# ============================================================
# Plot 2: Coupled leptogenesis Boltzmann
# dN_N/dz = -D * (N_N - N_N^eq)
# dN_{B-L}/dz = eps * D * (N_N - N_N^eq) - W * N_{B-L}
# D ~ K z K_1(z)/K_2(z),   W ~ K z^3 K_1(z)/4
# ============================================================
from scipy.special import kv

def D_W(z, K):
    Bz = kv(1, z) / kv(2, z)
    D_val = K * z * Bz
    W_val = 0.25 * K * z**3 * kv(1, z)
    return D_val, W_val

def N_eq_fn(z):
    # non-rel equilibrium for N
    return 0.75 * z**2 * kv(2, z) / 2

eps = 1e-6   # CP asymmetry
z_arr = np.linspace(0.1, 25, 1500)

ax2 = fig.add_subplot(gs[1, 0], facecolor=BG)

for K, col in zip([1, 10, 100], [GREEN, GOLD, WARN]):
    def rhs(y, z, K):
        NN, NBL = y
        D_val, W_val = D_W(z, K)
        Neq = N_eq_fn(z)
        dNN = -D_val * (NN - Neq)
        dNBL = eps * D_val * (NN - Neq) - W_val * NBL
        return [dNN, dNBL]
    y0 = [N_eq_fn(z_arr[0]), 0]
    sol = odeint(rhs, y0, z_arr, args=(K,))
    ax2.plot(z_arr, sol[:,0], lw=2, color=col, label=f'K={K}  N_N')
ax2.plot(z_arr, [N_eq_fn(z) for z in z_arr], color='white', lw=1.3, ls='--', label='N_N^eq')
ax2.set_xlabel('z = M_N / T', color=PURPLE, fontsize=10)
ax2.set_ylabel('N_N', color=PURPLE, fontsize=10)
ax2.set_yscale('log'); ax2.set_ylim(1e-8, 1.5)
ax2.set_title('N_1 population evolution', color='white', fontsize=11)
ax2.legend(fontsize=8, facecolor='#0f172a', labelcolor='white', edgecolor='#312e81')
ax2.tick_params(colors=PURPLE, labelsize=9)
ax2.grid(True, alpha=0.15, color=PURPLE)
for sp in ax2.spines.values(): sp.set_color('#312e81')

ax3 = fig.add_subplot(gs[1, 1], facecolor=BG)
for K, col in zip([1, 10, 100], [GREEN, GOLD, WARN]):
    def rhs(y, z, K):
        NN, NBL = y
        D_val, W_val = D_W(z, K)
        Neq = N_eq_fn(z)
        dNN = -D_val * (NN - Neq)
        dNBL = eps * D_val * (NN - Neq) - W_val * NBL
        return [dNN, dNBL]
    y0 = [N_eq_fn(z_arr[0]), 0]
    sol = odeint(rhs, y0, z_arr, args=(K,))
    ax3.plot(z_arr, np.abs(sol[:,1]/eps), lw=2, color=col, label=f'K={K}')
ax3.set_xlabel('z = M_N / T', color=PURPLE, fontsize=10)
ax3.set_ylabel('|N_(B-L)| / \u03B5', color=PURPLE, fontsize=10)
ax3.set_yscale('log')
ax3.set_title('B-L asymmetry vs K (efficiency)', color='white', fontsize=11)
ax3.legend(fontsize=8, facecolor='#0f172a', labelcolor='white', edgecolor='#312e81')
ax3.tick_params(colors=PURPLE, labelsize=9)
ax3.grid(True, alpha=0.15, color=PURPLE)
for sp in ax3.spines.values(): sp.set_color('#312e81')

# ============================================================
# Plot 3: Coleman-de Luccia bubble profile (thin-wall toy)
# phi'' + 2/r phi' = dV/dphi
# V(phi) = lam (phi^2 - v^2)^2 / 4 - eps phi (asymmetry)
# ============================================================
lam_s = 1.0
v = 1.0
eps_v = 0.2

def dVdphi(phi):
    return lam_s * phi * (phi**2 - v**2) - eps_v

def bounce_ode(y, r):
    phi, dphi = y
    if r < 1e-4:
        ddphi = dVdphi(phi)/3  # regularize origin
    else:
        ddphi = dVdphi(phi) - 2*dphi/r
    return [dphi, ddphi]

# shoot for initial phi_0 such that phi(r -> inf) -> false vacuum (phi = -v-ish)
# Simple shooting loop:
def shoot(phi0, r_max=20, steps=3000):
    r = np.linspace(1e-4, r_max, steps)
    sol = odeint(bounce_ode, [phi0, 0], r)
    return r, sol

# scan
best_phi0 = None
best_gap = 1e9
for phi0_try in np.linspace(0.8*v, 1.5*v, 200):
    r, sol = shoot(phi0_try, 15, 2000)
    # If overshoots (goes negative and oscillates), phi0 too large;
    # if undershoots (stays positive), too small.
    # Measure minimum absolute value far out
    tail = sol[-100:, 0]
    gap = abs(tail[-1] - (-v))
    if gap < best_gap:
        best_gap = gap; best_phi0 = phi0_try

r_best, sol_best = shoot(best_phi0, 15, 3000)

ax4 = fig.add_subplot(gs[2, 0], facecolor=BG)
ax4.plot(r_best, sol_best[:,0], color=INDIGO, lw=2.5, label='\u03C6(r)')
ax4.axhline(v,   color=GREEN, lw=1, ls=':', alpha=0.7)
ax4.axhline(-v,  color=WARN,  lw=1, ls=':', alpha=0.7)
ax4.text(0.5, 1.05, 'true vacuum',  color=GREEN, fontsize=9)
ax4.text(10,  -1.1, 'false vacuum', color=WARN,  fontsize=9)
ax4.set_xlabel('radius r [arb. units]', color=PURPLE, fontsize=10)
ax4.set_ylabel('\u03C6(r)', color=PURPLE, fontsize=10)
ax4.set_title('Critical bubble profile (shooting solution)', color='white', fontsize=11)
ax4.tick_params(colors=PURPLE, labelsize=9)
ax4.grid(True, alpha=0.15, color=PURPLE)
for sp in ax4.spines.values(): sp.set_color('#312e81')

# Bubble expansion  r(t) = r_c + v_w * t
ax5 = fig.add_subplot(gs[2, 1], facecolor=BG)
t = np.linspace(0, 20, 300)
for vw, col, lab in [(0.3, INDIGO, 'v_w = 0.3c'), (0.5, PURPLE, '0.5c'),
                     (0.7, GOLD, '0.7c'), (0.9, WARN, '0.9c')]:
    r_c = 1.0
    r = r_c + vw * t
    ax5.plot(t, r, lw=2, color=col, label=lab)
ax5.set_xlabel('time [arb.]', color=PURPLE, fontsize=10)
ax5.set_ylabel('bubble radius', color=PURPLE, fontsize=10)
ax5.set_title('Bubble expansion for various wall velocities', color='white', fontsize=11)
ax5.legend(fontsize=9, facecolor='#0f172a', labelcolor='white', edgecolor='#312e81')
ax5.tick_params(colors=PURPLE, labelsize=9)
ax5.grid(True, alpha=0.15, color=PURPLE)
for sp in ax5.spines.values(): sp.set_color('#312e81')

# ============ console output ============
print('=== Boltzmann freeze-out: final yields ===')
for K in K_values:
    Y0 = Y_eq(x_arr[0])
    Y_sol = odeint(boltzmann_rhs, Y0, x_arr, args=(K,))
    print(f'  K = {K:>6.2f}   ->  Y(inf) = {Y_sol[-1,0]:.4e}')

print('\n=== Leptogenesis: final |N_(B-L)|/eps ===')
for K in [1, 10, 100]:
    def rhs(y, z, K):
        NN, NBL = y
        D_val, W_val = D_W(z, K)
        Neq = N_eq_fn(z)
        return [-D_val*(NN - Neq), eps*D_val*(NN - Neq) - W_val*NBL]
    y0 = [N_eq_fn(z_arr[0]), 0]
    sol = odeint(rhs, y0, z_arr, args=(K,))
    print(f'  K = {K:>4d}   ->  |N_(B-L)|/eps = {abs(sol[-1,1])/eps:.4e}')

print(f'\n=== Critical bubble: initial phi_0 = {best_phi0:.4f}  (v = {v:.2f}) ===')

plt.savefig('output.png', dpi=110, bbox_inches='tight', facecolor=BG)
print('\nSaved output.png')

Click Run to execute the Python code

Code will be executed with Python 3 on the server

7. Deriving the Boltzmann Equation

Starting from the relativistic Boltzmann equation in FLRW spacetime,

\[p^\mu \frac{\partial f}{\partial x^\mu} - \Gamma^\mu_{\rho\sigma}\,p^\rho p^\sigma\frac{\partial f}{\partial p^\mu} = C[f],\]

integration over the physical momentum gives (after some work)

\[\dot n + 3 H n = -\int C[f]\,\frac{d^3 p}{(2\pi)^3}.\]

For a 2\u21922 process \(1+2 \leftrightarrow 3+4\) the collision integral becomes

\[\int C\,\frac{d^3 p_1}{(2\pi)^3} = \int d\Pi_1 d\Pi_2 d\Pi_3 d\Pi_4 (2\pi)^4\delta^4(p)\,|\mathcal M|^2(f_1 f_2 - f_3 f_4).\]

In the Maxwell-Boltzmann limit this reduces to the familiar \(-\langle\sigma v\rangle(n^2 - n_{eq}^2)\) form used throughout the course.

6b. Freeze-Out vs Freeze-In

Two qualitatively different out-of-equilibrium regimes for a species X:

Freeze-out (strong coupling)

X starts in equilibrium at high T and later decouples when \(\Gamma < H\). Final abundance is set by freeze-out temperature; final yield decreases with increasing coupling. Canonical for WIMP dark matter and heavy decaying X in GUT baryogenesis.

Freeze-in (weak coupling)

X never reaches equilibrium; its abundance slowly accumulates from rare interactions with the thermal bath. Final yield increases with coupling. Relevant for light sterile neutrinos, FIMP dark matter, and some baryogenesis scenarios with very weakly coupled portal fields.

Appendix: Entropy Production and Dilution

A late-decaying heavy particle can inject significant entropy into the plasma, diluting any pre-existing asymmetry. If a species of mass M and abundance Y_X decays at \(T_{dec}\), the entropy increase is

\[\frac{s_{\mathrm{after}}}{s_{\mathrm{before}}} \sim 1 + \frac{M\, Y_X}{T_{\mathrm{dec}}}.\]

For the generated baryon asymmetry, \(Y_B^{\mathrm{final}} = Y_B^{\mathrm{initial}}/\Delta\) where \(\Delta = s_{\mathrm{after}}/s_{\mathrm{before}}\). This is why moduli problems in SUSY cosmology and gravitino decays can be fatal for high-scale baryogenesis: they dilute \(\eta\) below the observed value.

8. Preview: Affleck-Dine Baryogenesis

One elegant alternative evades the usual out-of-equilibrium requirements by using coherent oscillations of a scalar field carrying baryon number. In supersymmetric extensions of the SM, flat directions of the scalar potential can hold large VEVs during inflation; after inflation, CP-violating terms in the potential cause the field to acquire an angular velocity in the complex plane, which corresponds to a baryon number density:

\[n_B = \mathrm{Im}(\phi^*\dot\phi) = |\phi|^2\,\dot\theta.\]

This mechanism (Affleck-Dine 1985) is an example of “coherent” baryogenesis — no decay, no phase transition, just classical scalar field dynamics. It can produce very large asymmetries naturally, often needing dilution rather than enhancement.

References

Kolb, E. W. & Turner, M. S., The Early Universe, Chapter 5 (Boltzmann equation).
Lee, B. W. & Weinberg, S., “Cosmological lower bound on heavy-neutrino masses”, Phys. Rev. Lett. 39, 165 (1977).
Gondolo, P. & Gelmini, G., “Cosmic abundances of stable particles”, Nucl. Phys. B 360, 145 (1991).
Coleman, S., “Fate of the false vacuum”, Phys. Rev. D 15, 2929 (1977).
Coleman, S. & de Luccia, F., “Gravitational effects on and of vacuum decay”, Phys. Rev. D 21, 3305 (1980).
Quiros, M., “Finite temperature field theory and phase transitions”, hep-ph/9901312 (1999).
Buchmuller, W., Di Bari, P. & Plumacher, M., “Leptogenesis for pedestrians”, Annals Phys. 315, 305 (2005).

← Module 3 Module 5: GUT Baryogenesis →