Homogeneous Infinite-Medium Reactor

Key Facts 

Read this before modifying the homogeneous solver.

Balance: \(\mathbf{A}\phi = \frac{1}{k}\mathbf{F}\phi\) where \(\mathbf{A} = \text{diag}(\Sigma_t) - \Sigma_s^T\), \(\mathbf{F} = \chi \otimes (\nu\Sigma_f)\)
1-group: \(k = \nu\Sigma_f / \Sigma_a\) (exact, no iteration)
Multi-group: \(k = \lambda_{\max}(\mathbf{A}^{-1}\mathbf{F})\) via numpy.eigvals
This is the reference eigenvalue for ALL solvers on homogeneous problems
Tolerance: < 1e-12 (limited only by FP arithmetic on small dense matrices)
Gotcha: this eigenvalue is flux-shape independent — it tests nothing about spatial or angular discretization

The infinite homogeneous medium is the simplest model in reactor physics. All spatial dependence vanishes (infinite geometry), all angular dependence integrates out (isotropic medium), and the neutron transport equation reduces to a pure energy balance. The only unknowns are the neutron energy spectrum \(\phi(E)\) and the infinite multiplication factor \(\kinf\).

Despite its simplicity, the homogeneous model is the foundation on which all other solvers build:

It is the first module students encounter in the ORPHEUS curriculum, introducing the multi-group eigenvalue problem and the power iteration algorithm.
The cross-section preparation pipeline — isotope loading, sigma-zero self-shielding, interpolation, macroscopic summation — is exercised here and reused unchanged by every subsequent solver (SN, MoC, CP, Monte Carlo, diffusion).
Analytical eigenvalues for 1-, 2-, and 4-group homogeneous media serve as verification benchmarks for all deterministic solvers.

This chapter derives the infinite-medium eigenvalue problem from first principles, describes the cross-section preparation pipeline, and presents the power iteration algorithm used to compute \(\kinf\) and \(\phi(E)\).

The solver is implemented in HomogeneousSolver, which satisfies the EigenvalueSolver protocol. The convenience wrapper solve_homogeneous_infinite() runs the full calculation and returns a HomogeneousResult.

From the Boltzmann Equation to the Infinite Medium 

The Boltzmann Transport Equation 

The starting point is the steady-state neutron transport equation in its integro-differential form [Duderstadt1976]:

(1)\[\hat{\Omega} \cdot \nabla \psi(\mathbf{r}, \hat{\Omega}, E) + \Sigma_\mathrm{t}(\mathbf{r}, E) \, \psi(\mathbf{r}, \hat{\Omega}, E) = \int_0^\infty \!\!\int_{4\pi} \Sigma_\mathrm{s}(\mathbf{r}, E' \!\to\! E, \hat{\Omega}' \!\to\! \hat{\Omega}) \, \psi(\mathbf{r}, \hat{\Omega}', E') \, d\Omega' \, dE' + \frac{\chi(E)}{4\pi \, k} \int_0^\infty \nu\Sigma_\mathrm{f}(\mathbf{r}, E') \, \phi(\mathbf{r}, E') \, dE'\]

Here \(\psi(\mathbf{r}, \hat{\Omega}, E)\) is the angular flux, \(\phi(\mathbf{r}, E) = \int_{4\pi} \psi \, d\Omega\) is the scalar flux, \(\chi(E)\) is the fission spectrum, and \(k\) is the multiplication factor eigenvalue.

Simplification for the Infinite Homogeneous Medium 

Three physical conditions dramatically simplify Eq. (1):

Infinite geometry — no boundaries, so the flux is spatially uniform: \(\nabla \psi = 0\). The streaming term vanishes entirely, and with it all leakage.
Homogeneous medium — all cross sections are independent of position: \(\Sigma_x(\mathbf{r}, E) = \Sigma_x(E)\).
Isotropy — in an infinite homogeneous medium with isotropic sources, the angular flux is isotropic: \(\psi(\hat{\Omega}, E) = \phi(E) / 4\pi\). The scattering kernel reduces to its \(P_0\) (isotropic) component.

After integrating over all directions, the transport equation collapses to a one-dimensional energy balance:

(2)\[\Sigt{} \phi(E) = \int_0^\infty \Sigma_{\mathrm{s},0}(E' \!\to\! E) \, \phi(E') \, dE' + \frac{\chi(E)}{k} \int_0^\infty \nu\Sigma_\mathrm{f}(E') \, \phi(E') \, dE'\]

where \(\Sigma_{\mathrm{s},0}\) is the isotropic scattering kernel.

Note

In the infinite homogeneous medium, scattering merely redistributes neutrons in energy. It does not change the total production-to-loss ratio, so \(\kinf\) depends only on the fission and absorption cross sections. Scattering does, however, determine the shape of the neutron spectrum \(\phi(E)\) — specifically the 1/E slowing-down region and the thermal Maxwellian peak.

Multi-Group Energy Discretisation 

Group-Averaged Cross Sections 

The continuous energy variable is discretised into \(G\) groups. Group 1 carries the highest energies (fast neutrons), group \(G\) the lowest (thermal neutrons). The energy boundaries \(E_0 > E_1 > \cdots > E_G\) define the grid stored in Mixture.eg.

The group flux is the integral over the group’s energy interval:

(3)\[\phi_g = \int_{E_g}^{E_{g-1}} \phi(E) \, dE\]

Group-averaged cross sections are flux-weighted averages:

(4)\[\Sigt{g} = \frac{1}{\phi_g} \int_{E_g}^{E_{g-1}} \Sigt{}(E) \, \phi(E) \, dE\]

In practice, these averages are pre-computed and stored in the 421-group HELIOS library that ships with ORPHEUS. The library provides cross sections tabulated at several background cross section (\(\sigma_0\)) values; the sigma-zero iteration (see Sigma-Zero Self-Shielding) selects the appropriate value for each isotope and group.

The Multi-Group Neutron Balance 

Substituting group-averaged quantities into Eq. (2) gives the multi-group neutron balance for group \(g\):

(5)\[\Sigt{g} \, \phi_g = \sum_{g'=1}^{G} \Sigs{g' \to g} \, \phi_{g'} + \frac{\chi_g}{k} \sum_{g'=1}^{G} \nSigf{g'} \, \phi_{g'}\]

The first term on the right is in-scattering from all groups (including self-scattering \(g' = g\)), and the second is the fission source weighted by the fission spectrum \(\chi_g\).

Matrix Form 

Collecting all \(G\) group equations into vectors and matrices:

(6)\[\mathbf{A} \, \boldsymbol{\phi} = \frac{1}{k} \, \mathbf{F} \, \boldsymbol{\phi}\]

where the removal matrix and fission matrix are:

(7)\[\mathbf{A} = \mathrm{diag}(\Sigt{g}) - \boldsymbol{\Sigma}_{\mathrm{s}}^T - 2 \, \boldsymbol{\Sigma}_2^T\]

(8)\[\mathbf{F} = \boldsymbol{\chi} \otimes \bigl(\nu\boldsymbol{\Sigma}_\mathrm{f} + 2 \, \text{colsum}(\boldsymbol{\Sigma}_2)\bigr)\]

Here \(\boldsymbol{\Sigma}_{\mathrm{s}}\) is the \(G \times G\) scattering transfer matrix (\(P_0\) component) and \(\boldsymbol{\Sigma}_2\) is the \((n,2n)\) transfer matrix.

Note

The \((n,2n)\) reaction appears in both matrices. In the removal matrix, \(-2\boldsymbol{\Sigma}_2^T\) removes the incident neutron and accounts for the two emitted neutrons entering the scattering system. In the fission matrix, the column-sum of \(\boldsymbol{\Sigma}_2\) adds the net production of one extra neutron per \((n,2n)\) event.

See HomogeneousSolver.__init__ for the construction of \(\mathbf{A}\) and HomogeneousSolver.compute_fission_source() for the production term.

The eigenvalue \(k = \kinf\) is the largest eigenvalue of the generalised problem (6). By the Perron–Frobenius theorem [Hebert2009], the dominant eigenvector \(\boldsymbol{\phi}\) is the unique non-negative solution — the fundamental mode — which is the physically meaningful neutron spectrum.

Scattering Matrix Convention 

The scattering transfer matrix \(\boldsymbol{\Sigma}_{\mathrm{s}}\) is stored in the from-row, to-column convention:

(9)\[(\boldsymbol{\Sigma}_{\mathrm{s}})_{g',g} = \Sigs{g' \to g}\]

That is, row \(g'\) gives the source group and column \(g\) gives the destination group. A downscatter-only matrix is therefore lower-triangular: non-zero entries only below (or on) the diagonal, because \(\Sigs{g' \to g} = 0\) when \(g < g'\) (no neutrons scatter from thermal to fast).

The neutron balance (5) requires the in-scattering into group \(g\) from all groups \(g'\):

\[\sum_{g'} \Sigs{g' \to g} \phi_{g'} = \bigl(\boldsymbol{\Sigma}_{\mathrm{s}}^T \cdot \boldsymbol{\phi}\bigr)_g\]

This is why the removal matrix (7) uses the transpose \(\boldsymbol{\Sigma}_{\mathrm{s}}^T\): the transpose converts column \(g\) (destination) to row \(g\), making the matrix-vector product give the in-scattering rate per group.

Warning

Getting the transpose wrong is a common source of bugs (see ERR-002 in the error catalog). For symmetric scattering matrices (e.g., 1-group self-scatter), the transpose is invisible, and the bug only manifests in multi-group problems with asymmetric down-scatter. This is why verification must always include \(\geq 2\) groups.

The Mixture stores SigS as a list of \(G \times G\) sparse matrices, one per Legendre order. The \(P_0\) component SigS[0] is used by the homogeneous solver; the higher orders are used by transport solvers with anisotropic scattering (SN, MoC).

Analytical Solutions 

One-Group Theory 

For a single energy group, the matrices reduce to scalars. The scattering terms cancel (a neutron scattered in group 1 remains in group 1), and the eigenvalue problem gives immediately:

(10)\[\kinf = \frac{\nu \Sigf{}}{\Siga{}}\]

This is the most fundamental result in reactor physics. It states that \(\kinf\) is the ratio of neutron production to neutron absorption, which is the definition of the infinite multiplication factor [Stacey2007].

For the connection to the four-factor formula: in a single-material homogeneous medium the thermal utilisation \(f = 1\), the resonance escape probability \(p = 1\) (no spatial heterogeneity), and the fast fission factor \(\varepsilon = 1\), so \(\kinf = \eta \cdot f \cdot p \cdot \varepsilon = \eta\).

Numerical example (from derivations.homogeneous.derive_1g()): \(\Sigt{} = 1.0\), \(\Sigma_\mathrm{c} = 0.2\), \(\Sigf{} = 0.3\), \(\nu = 2.5\), \(\Sigs{} = 0.5\) cm^-1:

\[\kinf = \frac{2.5 \times 0.3}{0.2 + 0.3} = 1.500000\]

Two-Group Theory 

For two energy groups (fast and thermal) with downscatter only (\(\chi = [1, 0]\), no upscatter from thermal to fast), the matrices are:

(11)\[\begin{split}\mathbf{A} = \begin{pmatrix} \Sigt{1} - \Sigs{1 \to 1} & 0 \\ -\Sigs{1 \to 2} & \Sigt{2} - \Sigs{2 \to 2} \end{pmatrix}\end{split}\]

(12)\[\begin{split}\mathbf{F} = \begin{pmatrix} \nu_1 \Sigf{1} & \nu_2 \Sigf{2} \\ 0 & 0 \end{pmatrix}\end{split}\]

Note that \(\mathbf{A}\) is lower-triangular because there is no upscatter (\(\Sigs{2 \to 1} = 0\)). This makes the inverse analytical:

(13)\[\begin{split}\mathbf{A}^{-1} = \begin{pmatrix} \dfrac{1}{\Sigma_{\mathrm{r},1}} & 0 \\[8pt] \dfrac{\Sigs{1 \to 2}}{\Sigma_{\mathrm{r},1} \, \Sigma_{\mathrm{r},2}} & \dfrac{1}{\Sigma_{\mathrm{r},2}} \end{pmatrix}\end{split}\]

where \(\Sigma_{\mathrm{r},g} = \Sigt{g} - \Sigs{g \to g}\) is the removal cross section for group \(g\) (total minus in-group scattering = absorption + out-scattering).

The eigenvalue matrix \(\mathbf{M} = \mathbf{A}^{-1}\mathbf{F}\) is:

(14)\[\begin{split}\mathbf{M} = \begin{pmatrix} \dfrac{\nu_1 \Sigf{1}}{\Sigma_{\mathrm{r},1}} & \dfrac{\nu_2 \Sigf{2}}{\Sigma_{\mathrm{r},1}} \\[8pt] \dfrac{\Sigs{1 \to 2}\, \nu_1 \Sigf{1}} {\Sigma_{\mathrm{r},1}\,\Sigma_{\mathrm{r},2}} + 0 & \dfrac{\Sigs{1 \to 2}\, \nu_2 \Sigf{2}} {\Sigma_{\mathrm{r},1}\,\Sigma_{\mathrm{r},2}} + \dfrac{\nu_2 \Sigf{2}}{\Sigma_{\mathrm{r},2}} \end{pmatrix}\end{split}\]

The characteristic equation \(\det(\mathbf{M} - \lambda\mathbf{I}) = 0\) gives a quadratic in \(\lambda\):

(15)\[\lambda^2 - \bigl(M_{11} + M_{22}\bigr)\lambda + \bigl(M_{11}M_{22} - M_{12}M_{21}\bigr) = 0\]

whose roots are:

(16)\[\lambda_{\pm} = \frac{(M_{11} + M_{22}) \pm \sqrt{(M_{11} - M_{22})^2 + 4 M_{12} M_{21}}}{2}\]

The dominant root \(\lambda_+\) is \(\kinf\).

Worked numerical example (from derivations.homogeneous.derive_2g()):

\(g\)	\(\Sigt{}\)	\(\Sigma_\mathrm{c}\)	\(\Sigf{}\)	\(\nu\)	\(\Sigs{g \to g}\)	\(\Sigs{1 \to 2}\)
1	0.50	0.01	0.01	2.50	0.38	0.10
2	1.00	0.02	0.08	2.50	0.90	—

The removal cross sections are \(\Sigma_{\mathrm{r},1} = 0.50 - 0.38 = 0.12\) and \(\Sigma_{\mathrm{r},2} = 1.00 - 0.90 = 0.10\).

The eigenvalue matrix entries are:

\[\begin{split}M_{11} &= \frac{2.50 \times 0.01}{0.12} = 0.208\overline{3} \\[4pt] M_{12} &= \frac{2.50 \times 0.08}{0.12} = 1.6\overline{6} \\[4pt] M_{21} &= \frac{0.10 \times 2.50 \times 0.01}{0.12 \times 0.10} = 2.083\overline{3} \\[4pt] M_{22} &= \frac{0.10 \times 2.50 \times 0.08}{0.12 \times 0.10} + \frac{2.50 \times 0.08}{0.10} = 1.6\overline{6} + 2.0 = 3.6\overline{6}\end{split}\]

Substituting into the quadratic formula (16):

\[\kinf = \lambda_+ = 1.8750000000\]

The second eigenvalue \(\lambda_- = 2.0833\overline{3}\) is also positive but smaller. The dominance ratio is \(|\lambda_-/\lambda_+| = 1.11\), which for this particular set of cross sections exceeds unity — indicating that the spectrum is initially far from the fundamental mode but converges rapidly because the power iteration eigenvalue still converges monotonically.

Note

The large \(\kinf\) in these analytical benchmarks reflects the synthetic cross sections chosen for verification, not a physical reactor. The cross sections are deliberately simple to enable exact symbolic solutions.

Four-Group Theory 

For four groups (fast, epithermal, thermal-1, thermal-2) with a full downscatter cascade and fission in all groups, the characteristic polynomial is degree 4 and has no convenient closed form. The analytical eigenvalue is computed numerically by SymPy’s symbolic eigenvalue solver applied to the \(4 \times 4\) matrix \(\mathbf{A}^{-1} \mathbf{F}\).

Result (from derivations.homogeneous.derive_4g()):

\[\kinf = 1.4877619048\]

Warning

The analytical eigenvalues are computed from the same matrix structure as the numerical solver. This is a code-verification test (does the code correctly implement the matrix algebra?), not a physics-validation test. Independent validation requires comparison to a different code or to experimental data (see the MATLAB reference values in the demo scripts).

Cross-Section Preparation 

Before any solver can run, the macroscopic cross sections must be assembled from isotopic data. This section describes the pipeline implemented in data.macro_xs, which is exercised for the first time in the homogeneous module and reused by all subsequent solvers.

Pipeline Overview 

The cross-section preparation follows five steps:

Load isotopes — read the 421-group microscopic cross-section library for each nuclide at the desired temperature.
Compute number densities — convert mass densities and compositions to number densities in the library’s unit system.
Sigma-zero iteration — find the self-consistent background cross section for each isotope and group (self-shielding).
Interpolate — evaluate microscopic cross sections at the converged sigma-zero values.
Sum to macroscopic — weight by number densities and sum over isotopes to obtain the Mixture.

This pipeline is encapsulated in compute_macro_xs().

Number Densities 

The atomic number density of species \(i\) (in \(1/(\text{barn} \cdot \text{cm})\)) is:

(17)\[N_i = \frac{\rho_i}{m_u \, A_i}\]

where \(\rho_i\) is the partial mass density in \(\text{g}/\text{cm}^3\), \(m_u = 1.660538 \times 10^{-24}\) g is the atomic mass unit, and \(A_i\) is the atomic weight. The factor \(10^{-24}\) converts the natural units (\(\text{cm}^{-3}\)) to the library units (\(1/(\text{barn} \cdot \text{cm})\)).

For aqueous solutions, the water density is obtained from the IAPWS-IF97 steam tables via pyXSteam. See aqueous_uranium() and pwr_like_mix().

Sigma-Zero Self-Shielding 

Cross sections in the resonance region depend strongly on the background cross section \(\sigma_{0,i,g}\) — a measure of how “dilute” isotope \(i\) is relative to its neighbours. The background cross section is defined as [Bondarenko1964]:

(18)\[\sigma_{0,i,g} = \frac{\Sigma_\mathrm{escape} + \displaystyle\sum_{j \ne i} N_j \, \sigma_{\mathrm{t},j,g}}{N_i}\]

where \(\Sigma_\mathrm{escape}\) is the escape cross section (zero for an infinite homogeneous medium) and the sum runs over all other isotopes in the mixture.

Physical meaning: when \(\sigma_0\) is large (dilute limit or strong moderator), the resonance peaks are fully resolved and the effective cross section is close to the infinite-dilution value. When \(\sigma_0\) is small (concentrated heavy absorber), the neutron flux is depressed at resonance energies — self-shielding — and the effective cross section is reduced.

The definition (18) is implicit: \(\sigma_{\mathrm{t},j,g}\) itself depends on \(\sigma_{0,j,g}\) through the library interpolation tables. The solution is obtained by fixed-point iteration:

Initialise \(\sigma_0\) to a large value (\(10^{10}\) barns, the infinite-dilution limit).
Interpolate \(\sigma_{\mathrm{t},j,g}\) from the library at the current \(\sigma_0\).
Recompute \(\sigma_0\) from Eq. (18).
Repeat until \(\|\sigma_0^{(n)} - \sigma_0^{(n-1)}\| < 10^{-6}\).

Convergence is fast (typically 3–5 iterations) because the dependence of \(\sigma_\mathrm{t}\) on \(\sigma_0\) is weak and monotonic. This is implemented in solve_sigma_zeros().

Note

For an infinite homogeneous medium, \(\Sigma_\mathrm{escape} = 0\). The sigma-zero depends only on the other isotopes in the mixture. For heterogeneous cells (fuel pins), the escape cross section \(\Sigma_e = \Sigma_\mathrm{pot} / \bar{\ell} \approx S/(4V)\) accounts for spatial self-shielding via the equivalence theory of Bondarenko [Bondarenko1964].

Cross-Section Interpolation 

The 421-group library tabulates microscopic cross sections at discrete \(\sigma_0\) base points (e.g., \(10^0, 10^1, \ldots, 10^{10}\) barns). Once the sigma-zero iteration converges, the cross section at the converged \(\sigma_0\) is obtained by log-linear interpolation in \(\log_{10}(\sigma_0)\) space:

(19)\[\sigma_{x,g}(\sigma_0) \approx \sigma_{x,g}(\sigma_0^{(a)}) + \frac{\log_{10} \sigma_0 - \log_{10} \sigma_0^{(a)}} {\log_{10} \sigma_0^{(b)} - \log_{10} \sigma_0^{(a)}} \bigl[\sigma_{x,g}(\sigma_0^{(b)}) - \sigma_{x,g}(\sigma_0^{(a)})\bigr]\]

where \(\sigma_0^{(a)}\) and \(\sigma_0^{(b)}\) are the bracketing base points. This is performed by interp_xs_field() for scalar cross sections and interp_sig_s() for scattering matrices.

Macroscopic Summation 

The macroscopic cross section for reaction \(x\) in group \(g\) is the density-weighted sum over all isotopes:

(20)\[\Sigma_{x,g} = \sum_{i=1}^{I} N_i \, \sigma_{x,i,g}\]

The following reaction types are assembled:

Attribute	Reaction	Notes
`SigC`	\((n,\gamma)\) capture	Radiative capture
`SigL`	\((n,\alpha)\) loss	Charged-particle emission
`SigF`	\((n,f)\) fission	Fission cross section
`SigP`	Production	\(\nu\Sigf{}\), summed over fissile isotopes only
`SigS`	Scattering matrices	One \(G \times G\) sparse matrix per Legendre order
`Sig2`	\((n,2n)\) matrix	\(G \times G\) sparse transfer matrix
`SigT`	Total	\(\Sigma_\mathrm{c} + \Sigma_\mathrm{L} + \Sigma_\mathrm{f} + \text{rowsum}(\Sigma_\mathrm{s}^{P_0}) + \text{rowsum}(\Sigma_2)\)
`chi`	Fission spectrum	Taken from first fissile isotope (simplification)

The absorption cross section used in the eigenvalue update (26) is not stored directly but computed as a derived property (absorption_xs):

(21)\[\Siga{g} = \Sigf{g} + \Sigma_{\mathrm{c},g} + \Sigma_{\mathrm{L},g} + \text{rowsum}(\boldsymbol{\Sigma}_{2,g})\]

This includes fission (neutron is absorbed to produce fission fragments), radiative capture \((n,\gamma)\), charged-particle emission \((n,\alpha)\), and the \((n,2n)\) reaction (where one neutron is “absorbed” and two are emitted, for a net gain of one).

The result is stored in a Mixture dataclass, which is the universal input to all ORPHEUS solvers.

Neutron Spectrum Physics 

The shape of the neutron energy spectrum \(\phi(E)\) in a homogeneous medium is controlled by the competition between moderation (slowing-down) and absorption. Three distinct energy regions are visible in the spectrum plots:

Fast Region (\(E > 0.1\) MeV)

Neutrons are born in fission with a spectrum peaked around 2 MeV. At these energies, scattering is nearly isotropic in the centre-of-mass frame and the mean logarithmic energy loss per collision with hydrogen is \(\xi = 1\). The fission source produces the characteristic fast peak.

For heavy nuclei like ²³⁸U, the elastic energy loss per collision is very small (\(\xi \approx 2/A\)), so the spectrum in the fast range is close to the fission spectrum \(\chi(E)\).

Slowing-Down Region (\(1\;\text{eV} < E < 0.1\;\text{MeV}\))

In this intermediate range, neutrons are slowed by elastic scattering (primarily with hydrogen). In the absence of absorption, the slowing-down equation yields the well-known 1/E flux law:

(22)\[\phi(E) = \frac{S}{\xi \Sigt{}} \cdot \frac{1}{E}\]

where \(S\) is the slowing-down source (neutrons entering from above) and \(\xi\) is the mean logarithmic energy decrement. On a flux-per-lethargy plot, the 1/E region appears as a horizontal plateau:

\[\frac{\phi}{du} = \frac{\phi(E) \cdot E}{\Delta u} \propto \frac{1}{E} \cdot E = \text{const}\]

This is why flux-per-lethargy is the standard representation: it makes the slowing-down region flat, and deviations (resonance dips from ²³⁸U, thermal peak) are immediately visible.

Resonance absorption (²³⁸U capture resonances) creates dips in the spectrum throughout this range. The sigma-zero self-shielding (see Sigma-Zero Self-Shielding) accounts for the flux depression in the resonance peaks.

Thermal Region (\(E < 1\) eV)

Below about 1 eV, neutrons reach thermal equilibrium with the moderator atoms. The thermal flux approaches a Maxwell–Boltzmann distribution at the moderator temperature \(T\):

(23)\[\phi_\mathrm{th}(E) \propto E \, \exp\!\left(-\frac{E}{k_B T}\right)\]

which peaks at \(E_\mathrm{peak} = k_B T\). At room temperature (294 K), \(k_B T = 0.0253\) eV, producing the characteristic thermal peak near 0.025 eV.

At higher moderator temperatures (e.g., 600 K in a PWR), the peak shifts to higher energies and broadens — this is Doppler broadening of the moderator distribution, which affects the thermal spectrum shape and hence the thermal cross sections.

Absorber poisons (e.g., boron) selectively remove thermal neutrons, depressing the thermal peak. This is clearly visible comparing the aqueous spectrum (no boron, strong thermal peak) with the PWR-like spectrum (4000 ppm B, suppressed thermal peak) in the Example Problems section below.

The Power Iteration Algorithm 

Algorithm 

The eigenvalue problem (6) is solved by power iteration [Hebert2009], converging to the dominant eigenvalue \(\kinf\) and the fundamental mode \(\boldsymbol{\phi}\). The algorithm implements the EigenvalueSolver protocol:

Initialise: \(\boldsymbol{\phi}^{(0)} = [1, 1, \ldots, 1]^T\), \(k^{(0)} = 1.0\).
Fission source (HomogeneousSolver.compute_fission_source()):

(24)\[\mathbf{Q}_f^{(n)} = \frac{\boldsymbol{\chi}}{k^{(n-1)}} \bigl(\boldsymbol{\Sigma}_\mathrm{p} + 2 \, \text{colsum}(\boldsymbol{\Sigma}_2)\bigr) \cdot \boldsymbol{\phi}^{(n-1)}\]
Fixed-source solve (HomogeneousSolver.solve_fixed_source()):

(25)\[\mathbf{A} \, \boldsymbol{\phi}^{(n)} = \mathbf{Q}_f^{(n)}\]

This is a single sparse direct solve via scipy.sparse.linalg.spsolve(). The removal matrix \(\mathbf{A}\) is constant across iterations (pre-built in HomogeneousSolver.__init__()).
Eigenvalue update (HomogeneousSolver.compute_keff()):

(26)\[k^{(n)} = \frac{\text{production}}{\text{absorption}} = \frac{\bigl(\boldsymbol{\Sigma}_\mathrm{p} + 2 \, \text{colsum}(\boldsymbol{\Sigma}_2)\bigr) \cdot \boldsymbol{\phi}^{(n)}} {\boldsymbol{\Sigma}_\mathrm{a} \cdot \boldsymbol{\phi}^{(n)}}\]

There is no leakage term in an infinite medium.
Convergence (HomogeneousSolver.converged()): stop when \(|k^{(n)} - k^{(n-1)}| < 10^{-10}\) after at least 3 iterations.

The generic loop is implemented in power_iteration().

Note

Unlike spatially-dependent solvers (SN, CP, diffusion), the homogeneous solver has no inner iteration: the removal matrix \(\mathbf{A}\) is constant and inverted directly. The only iteration is the outer power iteration on \(k\). This makes the homogeneous solver extremely fast — convergence in 3–5 iterations is typical.

Convergence Properties 

Power iteration converges to the dominant eigenvalue at a rate governed by the dominance ratio \(\rho = |k_1 / k_0|\), where \(k_0\) and \(k_1\) are the two largest eigenvalues of \(\mathbf{A}^{-1}\mathbf{F}\). After \(n\) iterations, the eigenvalue error decays as:

(27)\[|k^{(n)} - k_0| \sim \rho^n\]

For the 421-group industrial problems, the dominance ratio is very small (the spectrum is dominated by a single fundamental mode), so convergence is rapid. The following table shows the iteration history for the aqueous uranium case:

Iteration	\(k^{(n)}\)	\(\|k^{(n)} - k^{(n-1)}\|\)
1	1.03596	—
2	1.03596	\(< 10^{-5}\)
3	1.03596	\(< 10^{-10}\)

The Perron–Frobenius theorem guarantees that the converged eigenvector is the unique non-negative solution. This is critical for physical interpretability: the neutron flux must be non-negative everywhere in energy space, and the fundamental mode is the only eigenvector with this property.

Why homogeneous converges faster than spatially-dependent solvers:

In spatially-dependent problems (SN, diffusion), the dominance ratio approaches unity as the mesh is refined and the optical thickness increases, requiring hundreds of outer iterations. The homogeneous problem has no spatial mesh — the only degrees of freedom are the \(G\) energy groups — so the eigenvalue spectrum of \(\mathbf{A}^{-1}\mathbf{F}\) is typically well-separated, giving \(\rho \ll 1\).

Flux Normalisation 

The eigenvector \(\boldsymbol{\phi}\) is determined only up to a scalar multiple. After convergence, solve_homogeneous_infinite() normalises the flux so that the total neutron production rate is 100 n/cm³/s:

(28)\[\boldsymbol{\phi} \leftarrow \boldsymbol{\phi} \times \frac{100}{\bigl(\boldsymbol{\Sigma}_\mathrm{p} + 2\,\text{colsum}(\boldsymbol{\Sigma}_2)\bigr) \cdot \boldsymbol{\phi}}\]

Post-processing computes two spectral representations stored in HomogeneousResult:

Flux per unit energy: \(\phi_g / \Delta E_g\) (HomogeneousResult.flux_per_energy)
Flux per unit lethargy: \(\phi_g / \Delta u_g\) (HomogeneousResult.flux_per_lethargy)

where \(\Delta E_g = E_{g-1} - E_g\) and \(\Delta u_g = \ln(E_{g-1} / E_g)\).

Example Problems 

Aqueous Uranium Solution Reactor 

The simplest physical problem: water with dissolved uranium-235 (1000 ppm) at room temperature (294 K) and atmospheric pressure. This models a bare, infinite, aqueous homogeneous reactor — a configuration historically important for early criticality experiments.

The mixture contains only three isotopes: ¹H, ¹⁶O, and ²³⁵U. Water provides the moderation (hydrogen down-scatter) and ²³⁵U the fission source. The water density is obtained from the IAPWS-IF97 steam tables.

See aqueous_uranium().

import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

from data.macro_xs.recipes import aqueous_uranium
from homogeneous import solve_homogeneous_infinite

mix = aqueous_uranium(temp_K=294, pressure_MPa=0.1, u_conc_ppm=1000.0)
result = solve_homogeneous_infinite(mix)

fig, ax = plt.subplots()
ax.semilogx(result.eg_mid, result.flux_per_lethargy, 'b-', linewidth=1.2)
ax.set_xlabel('Energy (eV)')
ax.set_ylabel(r'Flux per unit lethargy $\phi / \Delta u$')
ax.set_title(
    rf'Aqueous U Solution — $k_\infty$ = {result.k_inf:.5f}'
)
ax.set_xlim(1e-3, 1e7)
ax.grid(True, alpha=0.3)
plt.tight_layout()

PWR-Like Homogenised Cell 

A more realistic problem: a PWR unit cell (UO₂ fuel, Zircaloy cladding, borated water) volume-homogenised into a single mixture. This is not a physically realisable configuration, but it exercises the full cross-section pipeline with 12 isotopes, self-shielding of ²³⁸U resonances, and boron absorption.

The geometric homogenisation uses volume fractions from the pin-cell geometry:

\[f_\mathrm{fuel} = \frac{r_\mathrm{fuel}^2}{r_\mathrm{cell}^2}, \quad f_\mathrm{clad} = \frac{r_\mathrm{clad,out}^2 - r_\mathrm{clad,in}^2} {r_\mathrm{cell}^2}, \quad f_\mathrm{cool} = \frac{r_\mathrm{cell}^2 - r_\mathrm{clad,out}^2} {r_\mathrm{cell}^2}\]

where \(r_\mathrm{cell} = p / \sqrt{\pi}\) is the Wigner–Seitz equivalent radius for a square lattice of pitch \(p\).

The mixture includes: ²³⁵U, ²³⁸U, ¹⁶O (fuel), five Zr isotopes (^{90,91,92,94,96}Zr), ¹H, ¹⁶O (coolant), ¹⁰B, ¹¹B.

See pwr_like_mix().

import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

from data.macro_xs.recipes import pwr_like_mix
from homogeneous import solve_homogeneous_infinite

mix = pwr_like_mix()
result = solve_homogeneous_infinite(mix)

fig, ax = plt.subplots()
ax.semilogx(result.eg_mid, result.flux_per_lethargy, 'r-', linewidth=1.2)
ax.set_xlabel('Energy (eV)')
ax.set_ylabel(r'Flux per unit lethargy $\phi / \Delta u$')
ax.set_title(
    rf'PWR-Like Homogenised Cell — $k_\infty$ = {result.k_inf:.5f}'
)
ax.set_xlim(1e-3, 1e7)
ax.grid(True, alpha=0.3)
plt.tight_layout()

Comparison 

Property	Aqueous U Solution	PWR-Like Mixture
\(\kinf\)	1.03596	1.01357
Fuel	Dissolved ²³⁵U (1000 ppm)	UO₂ (3% enrichment)
Moderator	Light water (294 K)	Borated water (600 K, 4000 ppm B)
Isotopes	3	12
Self-shielding	Negligible (²³⁵U dilute)	Significant (²³⁸U resonances)
MATLAB reference	1.03596	1.01357

Verification 

The homogeneous solver is verified against analytical eigenvalues derived symbolically with SymPy (see derivations.homogeneous). The same VerificationCase objects serve both the documentation (LaTeX equations in Verification Suite) and the test suite.

Benchmark	Groups	Regions	\(\kinf\) (analytical)	Error
`homo_1eg`	1	1	1.5000000000	\(< 10^{-12}\)
`homo_2eg`	2	1	1.8750000000	\(< 10^{-12}\)
`homo_4eg`	4	1	1.4877619048	\(< 10^{-12}\)

All three benchmarks achieve machine-precision agreement (\(< 10^{-12}\)), confirming that the solver correctly implements the matrix eigenvalue algebra.

Additionally, the two 421-group industrial problems (aqueous uranium and PWR-like mixture) are verified against the original MATLAB implementation to 5 significant digits.

Run the verification suite:

pytest tests/test_homogeneous.py -v

Comparison with Spatially-Dependent Solvers 

The homogeneous infinite-medium solver sits at the simplest end of the solver hierarchy. The following table compares it with the spatially-dependent solvers available in ORPHEUS:

Aspect	Homogeneous	Collision Probability	Discrete Ordinates	Diffusion
Spatial dependence	None	Region-averaged	Mesh-resolved	Mesh-resolved
Angular dependence	None (isotropic)	Integrated out	Discrete ordinates	Fick’s law
Transport operator	\(\mathbf{A}^{-1}\) (direct)	\(P_\infty\) matrix	Diamond-difference sweep	Implicit solve
Inner iterations	None	None	Scattering source	None
Typical convergence	3–5 outer	10–20 outer	20–50 outer	100+ outer
Eigenvalue computed	\(\kinf\)	\(\kinf\) (lattice)	\(\kinf\) (lattice)	\(\keff\) (core)
Implementation	`HomogeneousSolver`	`CPSolver`	`SNSolver`	`DiffusionSolver`

References 

[Duderstadt1976]

J.J. Duderstadt and L.J. Hamilton, Nuclear Reactor Analysis, Wiley, 1976.

[Bondarenko1964] (1,2)

I.I. Bondarenko et al., Group Constants for Nuclear Reactor Calculations, Consultants Bureau, 1964.

[Hebert2009] (1,2)

A. Hebert, Applied Reactor Physics, Presses internationales Polytechnique, 2009.

[Stacey2007]

W.M. Stacey, Nuclear Reactor Physics, 2nd ed., Wiley-VCH, 2007.