Step 1 — The valence space

Start by throwing almost everything away. ¹⁸O has 18 particles, but 16 of them form the closed ¹⁶O core — a filled-shell configuration so stable we declare it frozen and never let it move. What remains are 2 valence neutrons, and they are only allowed to sit in the three orbitals just above the core: quantum states a single particle can occupy, named d₅/₂, s₁/₂, d₃/₂ in spectroscopic notation (the letter is the orbital angular momentum $\ell$ — s means 0, d means 2 — and the subscript is the total angular momentum $j$, orbit and spin combined). Together these three are called the sd shell.

Why exactly these three? One fact closes off everything below, and one number makes everything above a small correction. Below: the core's 8 neutrons already fill the 1s and 1p orbitals completely, and the Pauli principle — no two identical fermions in the same state — leaves a valence neutron nothing to occupy down there. Above: the next orbitals (the pf shell) start about 5 MeV higher; that is the shell gap at $N=20$, the same kind of gap the spin-orbit demo opens before your eyes. Here is why a gap that size matters. A true ¹⁸O state is a blend of configurations — specific assignments of the two neutrons to orbitals — and how much a configuration contributes to a low-lying state is set by a ratio: how strongly the interaction couples it in, divided by how far uphill in energy it starts. The coupling is an MeV or so; a pf configuration starts 5+ MeV uphill; so its contribution enters at the few-percent level, squared in any probability. Small, but not zero — and here is the trick that makes the small space exact in practice: the interaction's numbers are chosen so that the sd-space calculation reproduces measured energies. The data already contain every cross-gap effect, so matrix elements tuned to the data carry those effects inside them. That is what effective interaction means, and it is why step 3's file is fit to experiment rather than derived from the bare force between two nucleons.

An orbital with angular momentum $j$ holds $2j+1$ m-states — its magnetic substates, one for each allowed projection $m$ of the angular momentum on the z-axis. Listing every (orbital, m) combination gives 12 "seats" for one species of particle, and that list is the entire universe of the calculation. A many-body state is a Slater determinant: the antisymmetric combination you get by filling seats, one particle per seat — which is the Pauli exclusion principle (no two identical fermions in the same state) enforced by construction. Since the force conserves the total z-projection $J_z$, we only need determinants whose m-values add up to the value we ask for — and counting them for ¹⁸O gives the first real number of the tutorial.

The program

"""Build your own shell model -- step 1: the valence space.

A shell-model calculation starts by throwing almost everything away. For a
nucleus like 18O we declare the 16O core FROZEN -- its 8 protons and 8 neutrons
never move -- and keep only the 2 leftover neutrons. Those two live in the
handful of orbitals just above the core: the sd shell, made of the orbitals
d5/2, s1/2, d3/2 (spectroscopic names: letter = orbital angular momentum
l = 2, 0, 2; subscript = total angular momentum j).

Each orbital with total angular momentum j holds 2j+1 magnetic substates,
m = -j ... +j. Listing every (orbital, m) combination for ONE species gives the
single-particle m-states -- the "seats" a valence neutron can occupy. That
list is the whole universe of this calculation.

The only other ingredient this step needs: a Slater determinant (the
antisymmetric many-body state built by filling seats, one particle per seat --
the Pauli principle as bookkeeping) is labeled by WHICH seats are filled.
Because the interaction conserves total Jz, we only need determinants whose
m-values add up to the Jz we ask for. Counting them for 18O is the first
real number the tutorial produces.
"""
from itertools import combinations

# label -> (name, l, 2j).  The numbering matches the interaction file coming
# in step 3 (usdb.int uses 1 = d3/2, 2 = d5/2, 3 = s1/2).
ORBITS = {1: ('d3/2', 2, 3), 2: ('d5/2', 2, 5), 3: ('s1/2', 0, 1)}

# every (orbital, m) seat for one species, in a fixed order
SEATS = []
for label, (name, l, j2) in ORBITS.items():
    for m2 in range(-j2, j2 + 1, 2):          # 2m runs -2j ... +2j in steps of 2
        SEATS.append((label, name, l, j2, m2))

print("the sd-shell single-particle m-states (one species):")
print(f"{'idx':>3} {'orbital':>8} {'l':>3} {'j':>5} {'m':>6}")
for i, (label, name, l, j2, m2) in enumerate(SEATS):
    print(f"{i:>3} {name:>8} {l:>3} {j2}/2 {m2:>4}/2")
print(f"total seats: {len(SEATS)}  (= (3+1) + (1+1) + (5+1) = 2j+1 summed)")

# 18O: two valence neutrons, total Jz = 0 -> pick 2 distinct seats whose
# m-values cancel. combinations() enforces the Pauli principle for free.
count = 0
for (s1, s2) in combinations(range(len(SEATS)), 2):
    if SEATS[s1][4] + SEATS[s2][4] == 0:
        count += 1
print(f"\n18O basis: 2 neutrons, Jz = 0  ->  {count} Slater determinants")
print("fourteen states. The whole 18O problem is a 14 x 14 matrix.")

browse the source