Thermal Relics and Abundance Calculations

A thermal relic remembers a competition between number-changing collisions and cosmic expansion. The useful conserved variable is not the physical number density but the number per unit bath entropy. This page derives its evolution, solves a representative cold-relic problem, and contrasts it with relativistic decoupling and later entropy dilution. The benchmark is a controlled single-species calculation, not a claim about the identity of dark matter.

Required background. States, Processes, and the Laws of Thermodynamics supplies Thermodynamics and Statistical Mechanics for equilibrium state variables, entropy, temperature, and the first law. Statistical Ensembles and Partition Functions and Quantum Statistical Mechanics and Density Operators develop the equilibrium distributions used below, while Thermodynamics in an Expanding Universe fixes the reference bath and effective-degree-of-freedom conventions.

Helpful background. The Einstein Equation supplies the gravitational convention, while The Friedmann Equations supplies the radiation-era $H(T)$. Relativistic Matter and Compact Objects supplies Relativistic Matter and Compact Objects thermodynamic and stress-energy context. Classical Field Modes and Canonical Variables supplies the Flat-Spacetime Quantum Field Theory mode convention, and Free Fields, Creation Operators, and Fock Space supplies occupation-state counting.

From number density to yield

Use $\hbar=c=k_{\rm B}=1$, cosmic time $t$, and $H=\dot a/a$. Consider a stable, self-conjugate species $\chi$ of mass $m$ and internal multiplicity $g_\chi$. Let $n$ sum over all those internal states, assume a homogeneous and isotropic population with negligible equilibrium chemical potential, and suppose kinetic equilibrium holds with a bath at temperature $T$. The annihilation products are assumed to remain equilibrated. For identical incoming particles, the event rate per physical volume is $R_{\rm ann}=\tfrac12\langle\sigma v_{\rm rel}\rangle n^2$; each event removes two $\chi$ particles, so the two factors cancel in the convention developed on Reaction Rates, Decoupling, and Freeze-Out:

$$\dot n+3Hn =-\langle\sigma v_{\rm rel}\rangle \left(n^2-n_{\rm eq}^2\right).$$

The plus sign in the dilution term follows from $a^{-3}d(a^3n)/dt$; the minus sign printed in Baudis and Profumo 2026, §27.3, printed p. 3, Eq. (27.2) (Official PDF) is a typographical error. Define

$$Y\equiv\frac{n}{s}, \qquad x\equiv\frac{m}{T}, \qquad s=\frac{2\pi^2}{45}g_{*s}(T)T^3.$$

Here $s$ is the entropy density of the dominant coupled reference bath; the tracked $\chi$ contribution is neglected in $g_{*s}$. Adiabatic bath evolution gives $\dot s+3Hs=0$, and hence

$$\dot Y=-s\langle\sigma v_{\rm rel}\rangle \left(Y^2-Y_{\rm eq}^2\right).$$

Changes in entropy degrees of freedom alter the clock. From $g_{*s}T^3a^3=\text{constant}$,

$$\frac{dx}{dt}=\frac{Hx}{\mathcal G(T)}, \qquad \mathcal G(T) =1+\frac13\frac{d\ln g_{*s}}{d\ln T}.$$

During radiation domination,

$$H=\sqrt{\frac{4\pi^3g_*}{45}}\frac{T^2}{M_{\rm Pl}} \simeq1.66\,g_*^{1/2}\frac{T^2}{M_{\rm Pl}}, \qquad M_{\rm Pl}=G^{-1/2}=1.22089\times10^{19}\ {\rm GeV},$$

where $g_*$ counts energy-density degrees of freedom and $M_{\rm Pl}$ is the nonreduced Planck mass Olive and Peacock 2026, §22.3.3, printed p. 13, Eq. (22.47) (Official PDF). Combining the equations gives

$$\boxed{ \frac{dY}{dx} =-\sqrt{\frac{\pi}{45}}\,\mathcal G(T) \frac{g_{*s}}{g_*^{1/2}} \frac{M_{\rm Pl}m}{x^2} \langle\sigma v_{\rm rel}\rangle \left(Y^2-Y_{\rm eq}^2\right) }.$$

The exact coefficient is $\sqrt{\pi/45}=0.2642218198\ldots$. The often-quoted $0.264$ is a useful rounded estimate; numerical values below retain the exact coefficient.

For a classical Maxwell–Boltzmann species with zero chemical potential, the equilibrium density is exact in $x$:

$$\begin{aligned} Y_{\rm eq}(x) &=\frac{45}{4\pi^4}\frac{g_\chi}{g_{*s}}x^2K_2(x),\\ Y_{\rm eq}(x) &\sim0.145\frac{g_\chi}{g_{*s}}x^{3/2}e^{-x} \left[1+\frac{15}{8x}+O(x^{-2})\right] \quad (x\gg1). \end{aligned}$$

This derivation shows exactly where adiabaticity, radiation domination, kinetic equilibrium, and the tabulated $g_*,g_{*s}$ functions enter. The phase-space reduction and identical-particle convention are derived in Gondolo and Gelmini 1991, §2, printed pp. 147–152, Eqs. (2.4)–(2.24) (author-uploaded full text).

The cold-relic asymptote

First take constant $g_*$ and $g_{*s}$ and an $s$-wave cross section $\langle\sigma v_{\rm rel}\rangle=\sigma_0$. Here $\sigma_0$ is the constant thermal-average coefficient, distinct from the scale factor $a(t)$. Then $\mathcal G=1$ and

$$\frac{dY}{dx} =-\frac{\lambda}{x^2}(Y^2-Y_{\rm eq}^2), \qquad \lambda=\sqrt{\frac{\pi}{45}}\frac{g_{*s}}{g_*^{1/2}}M_{\rm Pl}m\sigma_0.$$

After $Y_{\rm eq}$ becomes negligible, the equation is exactly separable:

$$\boxed{ \frac{1}{Y_\infty} =\frac{1}{Y_f}+\frac{\lambda}{x_f} }.$$

If residual annihilation dominates over the matching term, $Y_\infty\simeq x_f/\lambda$. The present density is

$$\Omega_\chi h^2 =\frac{m s_0Y_\infty}{\rho_c/h^2} \propto\frac{x_f}{\langle\sigma v_{\rm rel}\rangle}.$$

The mass cancels at leading order for constant $s$-wave annihilation, while $x_f$ retains a logarithmic mass and cross-section dependence. This is the origin—and the limit—of the inverse-cross-section rule. The hot- and cold-relic estimates, plus the dependence on entropy history and background expansion, are summarized in Baudis and Profumo 2026, §27.3, printed pp. 2–4, eqs. (27.3)–(27.6) (Official PDF).

A solved single-species benchmark

Take

$$m=100\ {\rm GeV},\quad g_\chi=2,\quad g_*=g_{*s}=90,\quad \sigma_0=2.00\times10^{-9}\ {\rm GeV}^{-2},\quad M_{\rm Pl}=1.22089\times10^{19}\ {\rm GeV}.$$

The thermal-average coefficient corresponds to $2.3347\times10^{-26}\ {\rm cm^3\,s^{-1}}$. Starting from $Y=Y_{\rm eq}$ at $x=1$, integrate the displayed ODE with the exact $K_2$ equilibrium density to $x_{\max}=10^4$. Use an implicit Radau solver with relative tolerance $10^{-10}$, absolute tolerance $10^{-22}$, and maximum step $5$, and locate $Y/Y_{\rm eq}=2$ with a dense-output root. Because $Y_{\rm eq}$ is negligible beyond the endpoint, complete the remaining $s$-wave tail analytically:

$$\frac{1}{Y_\infty} =\frac{1}{Y(x_{\max})}+\frac{\lambda}{x_{\max}}.$$

For the density conversion use $s_0=2891.2\ {\rm cm^{-3}}$ and $\rho_c/h^2=1.05375\times10^{-5}\ {\rm GeV\,cm^{-3}}$ for the final conversion, where $h\equiv H_0/(100\,\mathrm{km\,s^{-1}\,Mpc^{-1}})$. Here $s_0$ is the present entropy density of the standard bath; its photon and neutrino contributions follow from Olive and Peacock 2026, §22.3.1, printed p. 10, Eq. (22.40), and §22.3.3, printed p. 13, Eqs. (22.49)–(22.50), with the present photon temperature in §22.4.1, printed p. 17, Eq. (22.62) (Official PDF). The quoted critical density is a fixed conversion input defined by $\rho_c=3H_0^2/(8\pi G)$.

$\sigma_0$ in ${\rm GeV}^{-2}$	$x$ where $Y/Y_{\rm eq}=2$	Tail-corrected $Y_\infty$	$\Omega_\chi h^2$
$1.00\times10^{-9}$	$22.614$	$6.9999\times10^{-12}$	$0.19206$
$2.00\times10^{-9}$	$23.286$	$3.6104\times10^{-12}$	$0.099059$
$4.00\times10^{-9}$	$23.960$	$1.8604\times10^{-12}$	$0.051046$

Starting instead at $x=2$, $5$, or $10$ changes the central $Y_\infty$ by less than the declared $10^{-10}$ relative target, confirming that the solver tracks the equilibrium attractor before departure. Tightening to relative tolerance $10^{-11}$, absolute tolerance $10^{-24}$, and maximum step $2$, and applying the analytic tail after endpoints from $10^3$ through $3\times10^4$, leaves every tabulated digit unchanged. Halving and doubling $\sigma_0$ nearly invert the abundance; the small failure of exact inverse scaling is the predicted logarithmic shift of the departure point. Matching the exact post-freeze-out solution at $x_f=23.3$ gives $Y_\infty=3.58\times10^{-12}$, within one percent of the reduced-ODE integration.

This benchmark freezes $g_*$ and $g_{*s}$ and uses classical Maxwell–Boltzmann statistics. It numerically tests the reduced one-equation model and its scaling, not a particle model, a full precision freeze-out calculation, or a measured dark-matter abundance.

Relativistic decoupling and entropy dilution

If a species decouples while relativistic, annihilation never produces the nonrelativistic Boltzmann chase above. Its equilibrium yield simply freezes:

$$Y_\infty= \frac{45\zeta(3)}{2\pi^4} \frac{g_\chi}{g_{*s,d}} \begin{cases} 1, & \text{boson},\\ 3/4, & \text{fermion}, \end{cases}$$

Here $g_{*s,d}$ counts the still-coupled reference bath used in the denominator $s$, not the tracked relic. Thus the coefficients are $0.278\,g_\chi/g_{*s,d}$ for a boson and $0.208\,g_\chi/g_{*s,d}$ for a fermion. A sector already at temperature $T_\chi=\xi_dT$ acquires an additional factor $\xi_d^3$. Subsequent entropy redistribution within the coupled bath lowers $T_\chi/T$ but leaves the yield unchanged because the bath’s comoving entropy is conserved. In this hot-relic limit $\Omega h^2\propto m$, rather than $1/\langle\sigma v_{\rm rel}\rangle$.

If a later process increases the bath comoving entropy by $\Delta_S=S_{\rm after}/S_{\rm before}$ after $\chi$ has chemically decoupled, while $n_\chi a^3$ remains fixed, then

$$Y_{\rm after}=\frac{Y_{\rm before}}{\Delta_S}, \qquad \Omega_{\chi,\rm diluted}h^2 =\frac{\Omega_{\chi,\rm adiabatic}h^2}{\Delta_S}.$$

The density relation compares the final relic abundance with an adiabatic reference history at the same specified late-time bath entropy density and critical-density normalization. It does not say that the surviving particles lose mass or disappear at the instant entropy is injected: at a fixed scale factor their comoving number is unchanged, while their number per bath entropy is reduced.

That dilution is new physics relative to the adiabatic benchmark and must be stated rather than hidden in an “effective” annihilation cross section.

Where the one-equation treatment fails

Coannihilating species require a coupled system or a derived effective cross section with equilibrium weights. Resonances, thresholds, Sommerfeld enhancement, forbidden channels, an asymmetry, cannibal reactions, separate sector temperatures, kinetic nonequilibrium, nonstandard $H(T)$, or entropy injection can all change the scaling. A numerical calculation must supply the actual thermal average and $g_*,g_{*s}$ tables and demonstrate stability under the integration range and tolerances.

Exercises

Suppose the leading annihilation is $p$-wave, $\langle\sigma v_{\rm rel}\rangle=b/x$, with constant degrees of freedom. Derive the post-freeze-out abundance. Then find the effect of a later entropy increase $\Delta_S$, assuming it neither regenerates $\chi$ nor restarts annihilation, so $n_\chi a^3$ remains fixed.

Solution

Define

$$\lambda_b=\sqrt{\frac{\pi}{45}}\frac{g_{*s}}{g_*^{1/2}}M_{\rm Pl}mb.$$

After $Y_{\rm eq}$ is negligible,

$$\frac{dY}{dx}=-\frac{\lambda_b}{x^3}Y^2.$$

Integrating from $x_f$ to infinity gives

$$\boxed{ \frac{1}{Y_\infty} =\frac{1}{Y_f}+\frac{\lambda_b}{2x_f^2} }.$$

When the integral term dominates, $Y_\infty\simeq2x_f^2/\lambda_b$, so the leading density scales as $x_f^2/b$. Under the stated no-regeneration and no-reannihilation assumptions, a later entropy increase does not alter the already conserved particle number; it changes the yield to $Y_\infty/\Delta_S$ and dilutes $\Omega h^2$ by the same factor.

References

• Baudis, Laura, and Stefano Profumo. “Dark Matter.” In F. Takahashi et al. (Particle Data Group), Review of Particle Physics, International Journal of Modern Physics A 41, 2630011 (2026). Review revised August 2025; issue dated 1 June 2026. Official PDF.
• Gondolo, Paolo, and Graciela Gelmini. “Cosmic Abundances of Stable Particles: Improved Analysis.” Nuclear Physics B 360, 145–179 (1991). doi:10.1016/0550-3213(91)90438-4; author-uploaded full text.
• Olive, K. A., and J. A. Peacock. “Big-Bang Cosmology.” In F. Takahashi et al. (Particle Data Group), Review of Particle Physics, International Journal of Modern Physics A 41, 2630011 (2026). Review revised August 2025; issue dated 1 June 2026. Official PDF.