Neutron Diffusion Theory & Fermi Age Slowing Down

1. Phase Space and Angular Flux

A complete statistical description of a neutron gas requires a six-dimensional phase space: position $\vec{r} = (x, y, z)$ and velocity $\vec{v} = v \vec{\Omega}$ (or energy $E = \frac{1}{2} m v^2$ and unit direction vector $\vec{\Omega} \in S^2$). The angular neutron density $n(\vec{r}, \vec{\Omega}, E, t)$ is defined such that: $$n(\vec{r}, \vec{\Omega}, E, t) \, d^3 r \, d\Omega \, dE$$ is the expected number of neutrons in spatial volume $d^3 r$ around $\vec{r}$, with kinetic energies in $(E, E+dE)$, traveling in direction cone $d\Omega$ around $\vec{\Omega}$ at time $t$. The fundamental quantity of transport theory is the angular neutron flux $\psi$: $$\psi(\vec{r}, \vec{\Omega}, E, t) \equiv v \, n(\vec{r}, \vec{\Omega}, E, t)$$ Integrating over all $4\pi$ steradians yields the scalar neutron flux $\phi$: $$\phi(\vec{r}, E, t) = \int_{4\pi} \psi(\vec{r}, \vec{\Omega}, E, t) \, d\Omega = v \, n(\vec{r}, E, t)$$ The net neutron current density vector $\vec{J}$ is the first directional moment: $$\vec{J}(\vec{r}, E, t) = \int_{4\pi} \vec{\Omega} \, \psi(\vec{r}, \vec{\Omega}, E, t) \, d\Omega$$

2. The Steady-State Boltzmann Neutron Transport Equation

Balancing neutron gains and losses in an infinitesimal phase-space cell $d^3 r \, d\Omega \, dE$: $$\text{Streaming Loss} + \text{Collision Removal} = \text{In-Scattering Gain} + \text{Fission / External Source}$$ $$\mathbf{\vec{\Omega} \cdot \nabla \psi(\vec{r}, \vec{\Omega}, E) + \Sigma_t(\vec{r}, E) \psi(\vec{r}, \vec{\Omega}, E) = \int_{0}^{\infty} dE' \int_{4\pi} d\Omega' \, \Sigma_s(E' \to E, \vec{\Omega}' \to \vec{\Omega}) \psi(\vec{r}, \vec{\Omega}', E') + S(\vec{r}, \vec{\Omega}, E)}$$ This is the linear Boltzmann Integro-Differential Transport Equation. Because exact analytical solutions exist only for idealized infinite half-spaces, reactor physics employs systematic angular approximations—principally the Diffusion ($P_1$) Approximation.

1. The Spherical Harmonics $P_1$ Expansion

In media where scattering dominates absorption ($\Sigma_s \gg \Sigma_a$) and locations are several mean free paths away from boundaries and localized sources, the angular flux is nearly isotropic with a small directional bias: $$\psi(\vec{r}, \vec{\Omega}) \approx \frac{1}{4\pi} \phi(\vec{r}) + \frac{3}{4\pi} \vec{\Omega} \cdot \vec{J}(\vec{r})$$ This two-term expansion in Legendre polynomials ($P_0$ and $P_1$) is the $P_1$ approximation.

2. Derivation of Fick's Law of Diffusion

Substitute the $P_1$ expansion into the monoenergetic transport equation: $$\vec{\Omega} \cdot \nabla \left[ \frac{1}{4\pi} \phi + \frac{3}{4\pi} \vec{\Omega} \cdot \vec{J} \right] + \Sigma_t \left[ \frac{1}{4\pi} \phi + \frac{3}{4\pi} \vec{\Omega} \cdot \vec{J} \right] = \frac{1}{4\pi} \Sigma_s \phi + \frac{3}{4\pi} \Sigma_s \bar{\mu}_0 \vec{\Omega} \cdot \vec{J} + \frac{S}{4\pi}$$ To extract the vector current equation, multiply the entire transport equation by $\vec{\Omega}$ and integrate over all solid angles $\int_{4\pi} d\Omega$: Using the solid angle identities: $$\int_{4\pi} \vec{\Omega} \, d\Omega = 0, \qquad \int_{4\pi} \Omega_i \Omega_j \, d\Omega = \frac{4\pi}{3} \delta_{ij}$$ The streaming term becomes: $$\int_{4\pi} \vec{\Omega} (\vec{\Omega} \cdot \nabla \phi) \, d\Omega = \frac{4\pi}{3} \nabla \phi$$ The collision and scattering terms become: $$\frac{3}{4\pi} (\Sigma_t - \bar{\mu}_0 \Sigma_s) \int_{4\pi} \vec{\Omega} (\vec{\Omega} \cdot \vec{J}) \, d\Omega = (\Sigma_t - \bar{\mu}_0 \Sigma_s) \vec{J} \equiv \Sigma_{\text{tr}} \vec{J}$$ Equating both sides yields: $$\frac{1}{3} \nabla \phi + \Sigma_{\text{tr}} \vec{J} = 0 \implies \mathbf{\vec{J}(\vec{r}) = - \frac{1}{3 \Sigma_{\text{tr}}} \nabla \phi(\vec{r}) = - D \nabla \phi(\vec{r})}$$ This is Fick's Law of Neutron Diffusion! The diffusion coefficient $D$ is: $$\mathbf{D = \frac{1}{3 \Sigma_{\text{tr}}} = \frac{1}{3 \Sigma_s (1 - \bar{\mu}_0)} = \frac{\lambda_{\text{tr}}}{3}}$$

3. Validity Conditions for Fick's Law

Fick's law is valid under four strict physical conditions:

1. The medium is weakly absorbing: $\Sigma_a \ll \Sigma_s$.
2. Scattering is linearly anisotropic in the LAB frame ($\bar{\mu}_0 = 2/(3A)$).
3. The point of observation is at least $2\text{ to }3$ mean free paths away from strong localized sources.
4. The point of observation is at least $2\text{ to }3$ mean free paths away from vacuum boundaries or material interfaces.

1. The Continuity Equation

Consider an arbitrary volume $V$ bounded by surface $S$. At steady state, the rate of neutron loss must balance the rate of neutron production: $$\int_S \vec{J} \cdot d\vec{A} + \int_V \Sigma_a \phi \, dV = \int_V S \, dV$$ Applying Gauss's divergence theorem to the surface leakage integral: $$\int_V \left( \nabla \cdot \vec{J} + \Sigma_a \phi - S \right) dV = 0$$ Since this holds for any arbitrary volume, the differential neutron continuity equation is: $$\mathbf{\nabla \cdot \vec{J}(\vec{r}) + \Sigma_a \phi(\vec{r}) = S(\vec{r})}$$

2. The Steady-State Diffusion Equation

Substituting Fick's Law $\vec{J} = -D \nabla\phi$ into the continuity equation (assuming uniform diffusion coefficient $D$): $$\nabla \cdot (-D \nabla \phi) + \Sigma_a \phi = S(\vec{r})$$ $$\mathbf{-D \nabla^2 \phi(\vec{r}) + \Sigma_a \phi(\vec{r}) = S(\vec{r})}$$ Dividing through by $D$: $$\mathbf{-\nabla^2 \phi(\vec{r}) + \frac{1}{L^2} \phi(\vec{r}) = \frac{S(\vec{r})}{D}}$$ where $L$ is the fundamental thermal neutron diffusion length: $$\mathbf{L \equiv \sqrt{\frac{D}{\Sigma_a}} = \sqrt{\frac{1}{3 \Sigma_{\text{tr}} \Sigma_a}}}$$ In source-free regions ($S = 0$), the homogeneous diffusion equation reduces to the screened Poisson / Helmholtz form: $$\mathbf{\nabla^2 \phi(\vec{r}) - \frac{1}{L^2} \phi(\vec{r}) = 0}$$

1. Interface Continuity Conditions

At an interface between two distinct multiplying or moderating media (Medium 1 and Medium 2):

1. Continuity of Neutron Flux: $$\phi_1(\vec{r}_{\text{int}}) = \phi_2(\vec{r}_{\text{int}})$$ (Neutron density cannot jump discontinuously across a mathematical plane).
2. Continuity of Normal Current Density: $$-D_1 \left( \nabla \phi_1 \cdot \hat{n} \right) = -D_2 \left( \nabla \phi_2 \cdot \hat{n} \right)$$ (Neutrons are conserved; no neutrons can accumulate on a boundary of zero thickness).

2. Vacuum Boundary and Extrapolated Distance

At a vacuum boundary (outer surface of a reactor core adjacent to vacuum or air), no neutrons can enter the reactor from the outside: $$J_-(\vec{r}_s) = \int_{\vec{\Omega} \cdot \hat{n} < 0} |\vec{\Omega} \cdot \hat{n}| \psi(\vec{r}_s, \vec{\Omega}) \, d\Omega = 0$$ Using the $P_1$ partial currents formula: $$J_{\pm} = \frac{\phi}{4} \mp \frac{D}{2} \frac{d\phi}{dn} = 0 \implies \frac{\phi}{4} + \frac{D}{2} \frac{d\phi}{dn} = 0 \implies \frac{\phi}{|d\phi/dn|} = 2D = \frac{2}{3} \lambda_{\text{tr}}$$ Linear extrapolation predicts that the asymptotic flux vanishes at a distance $d$ outside the physical surface: $$d = \frac{2}{3} \lambda_{\text{tr}} \approx 0.667 \, \lambda_{\text{tr}} \quad (P_1\text{ approximation})$$

3. Exact Transport Result: Milne's Integral Problem

A rigorous transport-theoretic solution of the Milne problem (Wiener-Hopf method) demonstrates that the true asymptotic neutron flux vanishes outside the physical boundary at: $$\mathbf{d = 0.710446 \, \lambda_{\text{tr}} \approx 0.71 \, \lambda_{\text{tr}} = \frac{0.71}{\Sigma_{\text{tr}}}}$$ The extrapolated boundary $\tilde{R} = R + d$ allows reactor physicists to apply the simple Dirichlet boundary condition: $$\mathbf{\phi(\tilde{R}) = \phi(R + d) = 0}$$

1. Point Isotropic Source in an Infinite Medium

Consider a point source emitting $S$ monoenergetic neutrons per second at the origin $\vec{r} = 0$ in an infinite homogeneous moderator. In spherical coordinates with radial symmetry, the source-free diffusion equation for $r > 0$ is: $$\frac{1}{r^2} \frac{d}{dr}\left( r^2 \frac{d\phi}{dr} \right) - \frac{1}{L^2} \phi(r) = 0$$ Making the substitution $w(r) = r \phi(r)$: $$\frac{d^2 w}{dr^2} - \frac{1}{L^2} w = 0 \implies w(r) = A e^{-r/L} + C e^{+r/L}$$ Since physical flux must remain bounded as $r \to \infty$, $C = 0$: $$\phi(r) = \frac{A}{r} e^{-r/L}$$ To evaluate $A$, enforce source conservation as $r \to 0$: $$\lim_{r \to 0} 4\pi r^2 J(r) = \lim_{r \to 0} \left[ -4\pi r^2 D \frac{d\phi}{dr} \right] = 4\pi D A = S \implies A = \frac{S}{4\pi D}$$ The exact Green's function for a point source in 3D diffusion theory is: $$\mathbf{\phi(r) = \frac{S}{4\pi D r} e^{-r/L}}$$

2. Planar Infinite Sheet Source

For an infinite planar sheet source at $x = 0$ emitting $S_{\text{plane}}$ neutrons/$\text{cm}^2\cdot\text{s}$ into an infinite medium: $$\frac{d^2\phi}{dx^2} - \frac{1}{L^2} \phi = 0 \implies \mathbf{\phi(x) = \frac{S_{\text{plane}} L}{2 D} e^{-|x|/L}}$$

3. Physical Interpretation of Diffusion Length $L$

The mean square distance $\langle r^2 \rangle$ that a thermal neutron travels from its point of birth (thermalization) to its point of ultimate absorption is: $$\langle r^2 \rangle = \frac{\int_0^\infty r^2 \, \Sigma_a \phi(r) \cdot 4\pi r^2 dr}{\int_0^\infty \Sigma_a \phi(r) \cdot 4\pi r^2 dr} = \frac{\int_0^\infty r^3 e^{-r/L} dr}{\int_0^\infty r e^{-r/L} dr} = \frac{3! L^4}{1! L^2} = 6 L^2$$ Therefore: $$\mathbf{L^2 = \frac{1}{6} \langle r^2 \rangle \implies L = \sqrt{\frac{\langle r^2 \rangle}{6}}}$$ $L$ is directly proportional to the root-mean-square net displacement of a thermal neutron before absorption!

1. The Fermi Continuous Slowing Down Model

In intermediate and heavy moderators, neutrons undergo many collisions with small fractional energy loss, behaving as a continuous fluid slowing down from birth energy $E_0$ toward thermal energy $E_{\text{th}}$. Combining the slowing down density $q(\vec{r}, E) = \xi \Sigma_s E \phi(\vec{r}, E)$ with the spatial leakage $-D(E) \nabla^2 \phi$: $$\nabla \cdot \vec{J}(\vec{r}, E) + \frac{\partial q(\vec{r}, E)}{\partial u} = 0 \implies -D(E) \nabla^2 \phi + \frac{\partial q}{\partial u} = 0$$ Since $\phi = \frac{q}{\xi \Sigma_s E}$, we obtain: $$\frac{D(E)}{\xi \Sigma_s(E)} \nabla^2 q(\vec{r}, E) = -E \frac{\partial q}{\partial E}$$

2. Definition of Fermi Age $\tau$

Enrico Fermi defined the age parameter $\tau(E)$ (having dimensions of $\text{length}^2 = \text{cm}^2$): $$\mathbf{\tau(E) \equiv \int_{E}^{E_0} \frac{D(E')}{\xi \Sigma_s(E')} \frac{dE'}{E'} = \int_{0}^{u} \frac{D(u')}{\xi \Sigma_s(u')} du'}$$ By construction, $d\tau = - \frac{D(E)}{\xi \Sigma_s E} dE$. Substituting into the slowing down equation yields the celebrated Fermi Age Equation: $$\mathbf{\nabla^2 q(\vec{r}, \tau) = \frac{\partial q(\vec{r}, \tau)}{\partial \tau}}$$ This equation is mathematically isomorphic to Fourier's classical heat conduction equation $\nabla^2 T = \frac{1}{\alpha} \frac{\partial T}{\partial t}$, where Fermi age $\tau$ plays the role of time!

3. Gaussian Spatial Slowing Down Kernel

For a point burst of $S$ fission neutrons born at the origin at age $\tau = 0$ ($q(\vec{r}, 0) = S \delta(\vec{r})$), the solution of the Fermi age equation is the 3D Gaussian distribution: $$\mathbf{q(r, \tau) = \frac{S}{(4\pi \tau)^{3/2}} \exp\left( - \frac{r^2}{4\tau} \right)}$$ The mean square slowing-down distance from fission to thermal energy $\tau_{\text{th}}$ is: $$\langle r_s^2 \rangle = \frac{\int_0^\infty r^2 q(r, \tau) 4\pi r^2 dr}{\int_0^\infty q(r, \tau) 4\pi r^2 dr} = 6 \tau_{\text{th}}$$ Hence: $$\mathbf{\tau_{\text{th}} = \frac{1}{6} \langle r_s^2 \rangle \equiv L_s^2}$$ where $L_s$ is the slowing-down length of the moderator!

1. The Migration Area $M^2$

A neutron born in fission travels a distance while slowing down to thermal energy (characterized by Fermi age $\tau$), and then travels an additional distance as a thermal neutron before being absorbed (characterized by diffusion area $L^2$). The total mean square displacement from fission birth to ultimate absorption is the sum of the two independent random walks: $$\langle r_{\text{total}}^2 \rangle = \langle r_{\text{slowing}}^2 \rangle + \langle r_{\text{diffusion}}^2 \rangle = 6 \tau + 6 L^2$$ We define the Migration Area $M^2$: $$\mathbf{M^2 \equiv L^2 + \tau = \frac{1}{6} \langle r_{\text{total}}^2 \rangle}$$ The Migration Length $M$ is the square root: $$\mathbf{M \equiv \sqrt{M^2} = \sqrt{L^2 + \tau}}$$

2. Moderating Material Comparison

Key Physical Distinction: In Light Water, slowing down dominates migration ($M^2 \approx \tau = 27\text{ cm}^2$ vs $L^2 = 8\text{ cm}^2$). In contrast, in Heavy Water, thermal diffusion overwhelmingly dominates migration ($L^2 = 29{,}240\text{ cm}^2 \gg \tau = 131\text{ cm}^2$), because thermal neutrons wander for nearly two meters before being absorbed!