Nuclear level density, by counting the states

The exact integers the test compares against come from Euler's pentagonal-number recurrence, a different computation from this program's state count; posts/refs/nuclear-level-density-by-counting.py also enumerates the states of 12 fermions by brute force and finds the same numbers. The asymptotic ratios are evaluated there in 30-digit arithmetic.

What a level density is for

A heavy nucleus hit by a neutron is left with a few MeV of excitation energy, and at that energy it has an enormous number of states. How many states lie near a given energy, the level density $\rho (E)$, decides how likely it is to capture the neutron, and what it emits as it cools down. It enters reaction rates in reactors and in the formation of the heavy elements in stars.

The simplest model that gives the right trend is a Fermi gas: nucleons that do not interact with each other, sitting in single-particle levels. Here the levels are equally spaced, at energies $0,1,2,\dots$ in units of the spacing, and there is one kind of fermion, so each level holds at most one particle (the Pauli principle). With $N$ particles the lowest state fills levels $0$ to $N-1$. An excited state is any other choice of $N$ occupied levels, and its excitation energy $E$ is its total energy minus the lowest one. The program counts the choices at each $E$.

Counting without listing

Listing every choice of 100 levels out of 200 is out of the question, so the program counts with a table: after deciding about levels $0$ to $k$, entry $A[n,x]$ is the number of ways to have placed $n$ particles in them with energy $x$ above the least possible for $n$ particles, $n(n-1)/2$. Level $k$ is either empty, which changes nothing, or holds the $n$-th particle, which adds $k$ to the energy and $k-(n-1)$ to $x$. Going through $n$ from the top down uses each level at most once. Energies above $E_{\text{max}}$ are never needed, and no level above $N-1+E_{\text{max}}$ can be occupied at those energies, so the table is small.

The first counts can be checked by hand. $E=0$: one state, the ground state. $E=1$: one, the top particle up one level. $E=2$: two, the top particle up two, or the top two particles up one each. $E=3$: three. The pattern $1,1,2,3,5,\dots$ is the sequence of partition numbers $p(E)$, the number of ways to write $E$ as a sum of positive integers without regard to order ($4=3+1=2+2=2+1+1=1+1+1+1$).

The reason: number the particles from the top, and let ${\lambda }_i$ be how far the $i$-th one moved up. The Pauli principle keeps them in order, which forces ${\lambda }_1\ge {\lambda }_2\ge \cdots \ge 0$, and the ${\lambda }_i$ add up to $E$. So an excited state is a partition of $E$ into at most $N$ parts, and for $E\le N$ the limit on the number of parts never binds.

Bethe's formula

Hans Bethe worked out in 1936 how a Fermi gas's level density grows: $\rho (E)\propto e^{2\sqrt{aE}}$, where the level density parameter $a$ is ${\pi }^2/6$ times the number of single-particle levels per unit energy, which is 1 here. Hardy and Ramanujan had found the growth of the partition numbers in 1918, with the prefactor:

and $\pi \sqrt{2E/3}=2\sqrt{{\pi }^{2}E/6}$ is Bethe's exponent with $a={\pi }^2/6$. The two results are the same count.

Before looking: the formula is only the leading term of a series. Is it closer at $E=100$ or at $E=25$, and by how much?

At $E=25$ the count is 8.7% below the leading term and at $E=100$ it is 4.4% below (ratio_25, ratio_100): four times the energy, half the gap. The relative gap $1-p(E)/\text{HR}(E)$ times $\sqrt{E}$ is $0.437$ at $E=100$ (gap_times_root_e), and posts/refs/nuclear-level-density-by-counting.py follows it with exact integers to $E=6400$, where it is $0.4425$. The limit is the coefficient of the next term of the asymptotic series, $c=\frac{\sqrt{6}}{2\pi }+\frac{\pi }{24\sqrt{6}}=0.4433$ (C_NEXT), and $1-c/\sqrt{E}$ is the dashed curve in the second figure. At $E=100$ it misses the count by $0.0006$, and the miss falls like $1/E$.

Reading off the level density parameter

An experiment measures $\rho$, not $a$, and extracts $a$ from how fast $\ln \rho$ grows with $\sqrt{E}$. The exponent alone says the slope of $\ln \rho$ against $\sqrt{E}$ is $2\sqrt{a}=\pi \sqrt{2/3}=2.565$. The prefactor $1/E$ subtracts $2/\sqrt{E}$ from it. Before running: between $E=81$ and $E=100$, where $\sqrt{E}$ goes from $\sqrt{81}=9$ to $\sqrt{100}=10$, which of the two will the counts show?

The counts give a slope of $2.359$ (slope), not $2.565$: taking $2\sqrt{a}$ from them without the prefactor would put $a$ at $1.39$ instead of ${\pi }^2/6=1.645$, 15% low. With the $-2/\sqrt{E}$ term the prediction is $2.355$, within $0.005$. The prefactor matters at energies this low.