Chapter 4

Nuclear Sizes & Densities

Measuring Nuclear Sizes

The nucleus, with a radius on the order of femtometers ($1\;\text{fm} = 10^{-15}\;\text{m}$), is far too small to probe with visible light. Instead, we use high-energy electron scattering, which provides a clean electromagnetic probe of the nuclear charge distribution. The key insight of Robert Hofstadter's Nobel Prize-winning experiments (1961) was that the diffraction pattern of scattered electrons directly encodes the nuclear form factor.

The fundamental empirical result is that nuclear radii scale as:

$$R = r_0 A^{1/3}, \quad r_0 \approx 1.2\text{--}1.3 \;\text{fm}$$

This $A^{1/3}$ scaling implies that nuclear matter has roughly constant density, a property known as nuclear saturation. The volume scales as $V \propto A$, meaning each nucleon occupies approximately the same volume regardless of the nucleus size.

Electron Scattering and Form Factors

In Born approximation, the differential cross section for elastic electron scattering from a nucleus deviates from the point-nucleus (Mott) cross section by the square of the form factor:

$$\frac{d\sigma}{d\Omega} = \left(\frac{d\sigma}{d\Omega}\right)_{\text{Mott}} |F(q)|^2$$

The Mott cross section for a spin-1/2 electron scattering from a spinless point charge $Ze$ is:

$$\left(\frac{d\sigma}{d\Omega}\right)_{\text{Mott}} = \frac{Z^2 \alpha^2 (\hbar c)^2}{4 E^2 \sin^4(\theta/2)} \cos^2\!\left(\frac{\theta}{2}\right)$$

The form factor $F(\mathbf{q})$ is the Fourier transform of the normalized charge distribution. For a spherically symmetric distribution $\rho(r)$:

$$F(q) = \frac{4\pi}{Ze} \int_0^\infty \rho(r)\,\frac{\sin(qr)}{qr}\,r^2\,dr$$

where $q = |\mathbf{q}| = 2k\sin(\theta/2)$ is the magnitude of the three-momentum transfer and $k = E/(\hbar c)$ is the electron wave number. The minima in $|F(q)|^2$ correspond to diffraction minima and directly reveal the nuclear size.

Low-q Expansion

Expanding the form factor for small momentum transfer yields the model-independent determination of the RMS charge radius:

$$F(q) \approx 1 - \frac{q^2}{6}\langle r^2 \rangle + \cdots, \qquad \langle r^2 \rangle = \frac{\int \rho(r)\,r^4\,dr}{\int \rho(r)\,r^2\,dr}$$

Fermi (Woods-Saxon) Charge Distribution

Experimental data from electron scattering are well described by the two-parameter Fermi distribution (also called the Woods-Saxon distribution when used for the nuclear potential):

$$\rho(r) = \frac{\rho_0}{1 + \exp\!\left(\dfrac{r - c}{a}\right)}$$

The parameters have clear physical meaning:

Central density $\rho_0$

The density at $r = 0$. For nuclear matter (total nucleon density),$\rho_0 \approx 0.17\;\text{nucleons/fm}^3$. For the charge distribution,$\rho_{0,\text{ch}} \approx 0.07\text{--}0.08\;e/\text{fm}^3$. Remarkably, this is approximately constant for all nuclei with $A \gtrsim 20$.

Half-density radius $c$

The radius where $\rho(c) = \rho_0/2$. This scales as$c \approx 1.07\,A^{1/3}$ fm. It is closely related to (but not identical with) the RMS radius. The relationship between $c$ and $\langle r^2 \rangle^{1/2}$ depends on the diffuseness parameter.

Diffuseness $a$

Controls the surface thickness. The distance over which the density falls from 90% to 10% of $\rho_0$ is the surface thickness $t = 4a\ln 3 \approx 4.4a \approx 2.3$ fm. This is nearly constant across all nuclei with $a \approx 0.52\text{--}0.59$ fm.

The RMS radius for the Fermi distribution can be computed analytically in the limit $a \ll c$:

$$\langle r^2 \rangle = \frac{3}{5}c^2 + \frac{7}{5}\pi^2 a^2$$

Neutron Skin and Neutron Distributions

Electron scattering probes the charge (proton) distribution. The neutron distribution is more difficult to measure, requiring hadronic probes (proton scattering, pion scattering) or parity-violating electron scattering (PREX experiment at Jefferson Lab).

In neutron-rich nuclei, the neutron distribution extends beyond the proton distribution, forming a neutron skin. The neutron skin thickness is defined as:

$$\Delta r_{np} = \langle r^2 \rangle_n^{1/2} - \langle r^2 \rangle_p^{1/2}$$

The PREX-II experiment on $^{208}$Pb measured $\Delta r_{np} = 0.283 \pm 0.071$ fm, which has profound implications for the equation of state of neutron-rich matter and neutron star structure. The neutron skin thickness is closely correlated with the slope of the symmetry energy $L$ at saturation density:

$$S(\rho) \approx S_0 + \frac{L}{3}\left(\frac{\rho - \rho_0}{\rho_0}\right) + \cdots$$

A larger $L$ (stiffer symmetry energy) pushes neutrons outward, producing a thicker neutron skin and predicting larger neutron star radii.

Isotope Shifts and Muonic Atoms

Complementary methods for measuring nuclear charge radii include optical isotope shifts and muonic X-ray spectroscopy.

Optical Isotope Shifts

Atomic energy levels are sensitive to the nuclear charge radius through the finite nuclear size correction. For s-electrons in a hydrogen-like atom:

$$\Delta E_{\text{fs}} = \frac{2}{3}\frac{Z\alpha^2}{n^3}\left(\frac{Zm_e c}{2\hbar}\right)^2 \langle r^2 \rangle \cdot m_e c^2$$

The change in this correction between two isotopes gives the differential mean-square charge radius $\delta\langle r^2 \rangle^{A,A'}$. Modern laser spectroscopy achieves precisions of $\sim 10^{-3}$ fm$^2$ on these differences.

Muonic Atoms

The muon ($m_\mu \approx 207\,m_e$) orbits much closer to the nucleus, with Bohr radius$a_\mu = a_0\,m_e/m_\mu \approx 256$ fm for hydrogen. For heavy atoms, the muonic 1s orbit penetrates deep into the nucleus, making muonic X-ray transitions extremely sensitive to the charge distribution. The 2p → 1s transition energy shifts are:

$$\Delta E \propto Z^4 \alpha^4 m_\mu c^2 \left(\frac{R}{a_\mu}\right)^2$$

The 2010 proton radius puzzle arose when muonic hydrogen measurements yielded$r_p = 0.84184(67)$ fm, significantly smaller than the electron scattering value of $r_p = 0.8768(69)$ fm. This discrepancy has since been largely resolved, with the muonic value now accepted as more precise.

Nuclear Matter Density and Saturation

The near-constancy of nuclear central density is one of the most fundamental properties of nuclear matter. The saturation density is:

$$\rho_0 \approx 0.17\;\text{nucleons/fm}^3 \approx 2.8 \times 10^{14}\;\text{g/cm}^3$$

This enormous density (comparable to the density inside neutron stars) results from the balance between attractive nuclear forces and short-range repulsion. The saturation property means the binding energy per nucleon and the density are approximately independent of $A$ for $A \gtrsim 20$.

The nuclear equation of state near saturation can be expanded as:

$$\frac{E}{A}(\rho) \approx \frac{E}{A}(\rho_0) + \frac{K_\infty}{18}\left(\frac{\rho - \rho_0}{\rho_0}\right)^2 + \cdots$$

where $K_\infty \approx 230 \pm 20$ MeV is the nuclear incompressibility, determined from the energies of isoscalar giant monopole resonances (the nuclear breathing mode). This parameter constrains models of supernovae and neutron star mergers.

Python Simulation: Charge Distributions and Form Factors

Computes Fermi charge distributions, form factors via spherical Bessel transform, and RMS radii for several nuclei. Demonstrates the $R \propto A^{1/3}$ scaling.

Nuclear Charge Distributions

Python

Fermi distributions, form factors |F(q)|^2, and RMS charge radii

script.py130 lines

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

# ──────────────────────────────────────────────
# Nuclear Charge Distributions and Form Factors
# Electron scattering analysis
# ──────────────────────────────────────────────

r = np.linspace(0.001, 12.0, 2000)  # radial distance in fm
q = np.linspace(0.01, 4.0, 500)     # momentum transfer in fm^-1

# Fermi (Woods-Saxon) charge distribution parameters for various nuclei
nuclei = {
    'He-4':   {'Z': 2,  'A': 4,   'c': 0.964, 'a': 0.322, 'color': '#22c55e'},
    'C-12':   {'Z': 6,  'A': 12,  'c': 2.355, 'a': 0.522, 'color': '#3b82f6'},
    'Ca-40':  {'Z': 20, 'A': 40,  'c': 3.766, 'a': 0.586, 'color': '#f59e0b'},
    'Zr-90':  {'Z': 40, 'A': 90,  'c': 4.272, 'a': 0.523, 'color': '#ef4444'},
    'Pb-208': {'Z': 82, 'A': 208, 'c': 6.624, 'a': 0.549, 'color': '#a855f7'},
}

def fermi_distribution(r_arr, rho0, c, a):
    """Fermi (Woods-Saxon) charge distribution"""
    return rho0 / (1.0 + np.exp((r_arr - c) / a))

def form_factor_numerical(q_arr, r_arr, rho_arr):
    """Compute form factor F(q) from charge distribution via spherical Bessel transform"""
    dr = r_arr[1] - r_arr[0]
    F = np.zeros(len(q_arr))
    # Normalization: integral of rho(r) * 4*pi*r^2 dr
    norm = 4.0 * np.pi * np.trapezoid(rho_arr * r_arr**2, r_arr)
    for i, qi in enumerate(q_arr):
        if qi < 1e-10:
            F[i] = 1.0
        else:
            integrand = rho_arr * r_arr**2 * np.sin(qi * r_arr) / (qi * r_arr)
            F[i] = 4.0 * np.pi * np.trapezoid(integrand, r_arr) / norm
    return F

def rms_radius(r_arr, rho_arr):
    """Compute RMS charge radius from distribution"""
    num = np.trapezoid(rho_arr * r_arr**4, r_arr)
    den = np.trapezoid(rho_arr * r_arr**2, r_arr)
    return np.sqrt(num / den)

fig, axes = plt.subplots(1, 3, figsize=(18, 5.5))
fig.patch.set_facecolor('#0a0a1a')

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

print("Nuclear Charge Distribution Analysis")
print("=" * 55)
print(f"{'Nucleus':<10} {'c (fm)':<10} {'a (fm)':<10} {'R_rms (fm)':<12} {'R_rms/A^(1/3)'}")
print("-" * 55)

for name, params in nuclei.items():
    c = params['c']
    a = params['a']
    A = params['A']
    col = params['color']

# Normalize rho0 so that integral of rho * 4*pi*r^2 dr = Z
    rho_raw = fermi_distribution(r, 1.0, c, a)
    norm_integral = 4.0 * np.pi * np.trapezoid(rho_raw * r**2, r)
    rho0 = params['Z'] / norm_integral
    rho = fermi_distribution(r, rho0, c, a)

# RMS radius
    r_rms = rms_radius(r, rho)

# Form factor
    F = form_factor_numerical(q, r, rho)

# Plot charge distribution (normalized to central density)
    rho_central = rho0
    axes[0].plot(r, rho / rho_central, color=col, linewidth=2, label=name)

# Plot form factor
    F_abs = np.abs(F)
    F_abs[F_abs < 1e-6] = 1e-6
    axes[1].semilogy(q, F_abs**2, color=col, linewidth=2, label=name)

# RMS radius vs A^(1/3)
    axes[2].plot(A**(1.0/3.0), r_rms, 'o', color=col, markersize=10, zorder=5)
    axes[2].annotate(name, (A**(1.0/3.0), r_rms), textcoords="offset points",
                     xytext=(5, 5), color=col, fontsize=8)

print(f"{name:<10} {c:<10.3f} {a:<10.3f} {r_rms:<12.3f} {r_rms/A**(1.0/3.0):.3f}")

# Linear fit for R = r0 * A^(1/3)
A_vals = np.array([p['A'] for p in nuclei.values()])
r_rms_vals = np.array([rms_radius(r, fermi_distribution(r, 1.0, p['c'], p['a']))
                        for p in nuclei.values()])
A13 = A_vals**(1.0/3.0)
r0_fit = np.polyfit(A13, r_rms_vals, 1)[0]
A13_line = np.linspace(1.0, 6.5, 100)
axes[2].plot(A13_line, r0_fit * A13_line, '--', color='#94a3b8', linewidth=1.5,
             label=f'$R = {r0_fit:.2f} A^{{1/3}}$ fm')

axes[0].set_xlabel('r [fm]', color='#94a3b8', fontsize=12)
axes[0].set_ylabel(r'$\rho(r) / \rho_0$', color='#94a3b8', fontsize=12)
axes[0].set_title('Charge Distributions', color='white', fontsize=13)
axes[0].set_xlim(0, 10)
axes[0].legend(fontsize=8, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[0].tick_params(colors='#94a3b8')

axes[1].set_xlabel('q [fm$^{-1}$]', color='#94a3b8', fontsize=12)
axes[1].set_ylabel('$|F(q)|^2$', color='#94a3b8', fontsize=12)
axes[1].set_title('Form Factors', color='white', fontsize=13)
axes[1].set_ylim(1e-6, 1.5)
axes[1].legend(fontsize=8, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[1].tick_params(colors='#94a3b8')

axes[2].set_xlabel('$A^{1/3}$', color='#94a3b8', fontsize=12)
axes[2].set_ylabel('$R_{rms}$ [fm]', color='#94a3b8', fontsize=12)
axes[2].set_title('RMS Radii vs $A^{1/3}$', color='white', fontsize=13)
axes[2].legend(fontsize=9, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[2].tick_params(colors='#94a3b8')

plt.tight_layout()
plt.savefig('charge_distributions.png', dpi=130, bbox_inches='tight',
            facecolor=fig.get_facecolor())

print(f"\nBest-fit r0 = {r0_fit:.3f} fm")
print(f"Expected r0 ~ 1.2 fm for charge radii")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran Implementation: Nuclear Size Systematics

Computes RMS charge radii and surface thickness for a range of nuclei using Fermi distribution parameters from electron scattering data.

Nuclear Size Calculator

Fortran

Tabulates nuclear charge radii, central densities, and surface parameters

nuclear_sizes.f9086 lines

program nuclear_sizes
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  real(dp), parameter :: pi = 3.14159265358979323846_dp
  integer, parameter :: nr = 2000
  real(dp) :: r(nr), rho(nr), dr
  real(dp) :: c, a, rho0, A_mass, Z_charge
  real(dp) :: integral_r2, integral_r4, r_rms, r_half
  real(dp) :: surface_thickness
  integer :: i, inuc

! Nuclear parameters: c (half-density radius), a (diffuseness)
  integer, parameter :: n_nuclei = 6
  real(dp) :: c_arr(n_nuclei), a_arr(n_nuclei), A_arr(n_nuclei), Z_arr(n_nuclei)
  character(len=8) :: names(n_nuclei)

names = (/'He-4    ', 'O-16    ', 'Ca-40   ', 'Ca-48   ', 'Sn-120  ', 'Pb-208  '/)
  c_arr = (/0.964_dp, 2.608_dp, 3.766_dp, 3.742_dp, 4.652_dp, 6.624_dp/)
  a_arr = (/0.322_dp, 0.513_dp, 0.586_dp, 0.588_dp, 0.563_dp, 0.549_dp/)
  A_arr = (/4.0_dp, 16.0_dp, 40.0_dp, 48.0_dp, 120.0_dp, 208.0_dp/)
  Z_arr = (/2.0_dp, 8.0_dp, 20.0_dp, 20.0_dp, 50.0_dp, 82.0_dp/)

dr = 15.0_dp / real(nr, dp)
  do i = 1, nr
    r(i) = (real(i, dp) - 0.5_dp) * dr
  end do

write(*,'(A)') '=== Nuclear Sizes and Densities ==='
  write(*,'(A)') ''
  write(*,'(A)') 'Fermi distribution: rho(r) = rho0 / (1 + exp((r-c)/a))'
  write(*,'(A)') ''
  write(*,'(A,A8,A10,A10,A12,A12,A14)') ' ', 'Nucleus', 'c [fm]', 'a [fm]', &
       'R_rms [fm]', 'R/A^(1/3)', 't_90-10 [fm]'
  write(*,'(A)') ' ' // repeat('-', 70)

do inuc = 1, n_nuclei
    c = c_arr(inuc)
    a = a_arr(inuc)
    A_mass = A_arr(inuc)
    Z_charge = Z_arr(inuc)

! Compute Fermi distribution
    do i = 1, nr
      rho(i) = 1.0_dp / (1.0_dp + exp((r(i) - c) / a))
    end do

! Normalize: integral of rho * 4*pi*r^2 dr = Z
    integral_r2 = 0.0_dp
    do i = 1, nr
      integral_r2 = integral_r2 + rho(i) * r(i)**2 * dr
    end do
    rho0 = Z_charge / (4.0_dp * pi * integral_r2)

do i = 1, nr
      rho(i) = rho0 / (1.0_dp + exp((r(i) - c) / a))
    end do

! RMS radius
    integral_r2 = 0.0_dp
    integral_r4 = 0.0_dp
    do i = 1, nr
      integral_r2 = integral_r2 + rho(i) * r(i)**2 * dr
      integral_r4 = integral_r4 + rho(i) * r(i)**4 * dr
    end do
    r_rms = sqrt(integral_r4 / integral_r2)

! Surface thickness: 90% to 10% density points
    ! t = 4 * ln(3) * a
    surface_thickness = 4.0_dp * log(3.0_dp) * a

write(*,'(A,A8,F10.3,F10.3,F12.3,F12.3,F14.3)') ' ', names(inuc), c, a, &
         r_rms, r_rms / A_mass**(1.0_dp/3.0_dp), surface_thickness
  end do

write(*,'(A)') ''
  write(*,'(A)') 'Key results:'
  write(*,'(A)') '  - R = r0 * A^(1/3)  with r0 ~ 1.2 fm for charge radii'
  write(*,'(A)') '  - Central density rho0 ~ 0.08 e/fm^3 (charge)'
  write(*,'(A)') '  - Surface thickness t ~ 2.3 fm (nearly constant)'
  write(*,'(A)') '  - Neutron skin in Ca-48: ~0.12-0.15 fm'
  write(*,'(A)') ''
  write(*,'(A)') 'Neutron skin comparison (Ca-40 vs Ca-48):'
  write(*,'(A,F8.3,A)') '  Ca-40 R_ch = ', c_arr(3), ' fm (N=Z, minimal skin)'
  write(*,'(A,F8.3,A)') '  Ca-48 R_ch = ', c_arr(4), ' fm (8 extra neutrons)'
  write(*,'(A)') '  Neutron skin arises from excess neutrons extending beyond protons'
end program nuclear_sizes

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← Nuclear Models The Deuteron →