Monte Carlo Neutron Transport

Key Facts

Read this before modifying the MC solver.

• Power iteration MC: inactive cycles (source convergence) + active cycles (tallying)

• Woodcock delta-tracking: sample distance with \(\Sigma_{\max}\), reject with \(1 - \Sigma_t(r)/\Sigma_{\max}\)

• Analog absorption: on absorption, weight × \(\nu\Sigma_f/\Sigma_a\), reborn from \(\chi\)

• Russian roulette (w < threshold) and splitting (w > threshold) for population control

• keff estimator: \(k = \sum w_{\text{end}} / \sum w_{\text{start}}\) per cycle

• Statistical: \(\sigma \sim 1/\sqrt{N_{\text{active}}}\) (CLT). Use z-score for verification.

• Geometry: ConcentricPinCell (cylindrical) and SlabPinCell with periodic BC

• Gotcha: 1G homogeneous has σ = 0 (all neutrons see identical XS) — not a useful test

• Gotcha: pitch formula is 2 * r_outer, NOT 4 * area (ERR-017: 24% error + NaN)

• Direction sampling must be truly isotropic: cos(θ) = 1 - 2ξ, φ = 2πξ (ERR-018)

• Weight ratio keff = sum(w_end) / sum(w_start) must equal cycle keff (consistency check)

• Verification uses synthetic cross sections, not real nuclear data

Conventions

• Scattering matrix: Scattering Matrix Convention — SigS[g_from, g_to]

• Multi-group balance: (5) in Homogeneous Infinite-Medium Reactor

• Cross sections: Cross-Section Data Pipeline

• Verification: Cross-Section Library — regions A/B/C/D

• Eigenvalue: stochastic power iteration (distinct from deterministic)

Overview

The Monte Carlo (MC) method solves the neutron transport equation by simulating individual neutron histories stochastically. Rather than discretising the phase space (angle, energy, space) as deterministic methods do, MC tracks particles through a geometry, sampling collision events from probability distributions derived from the underlying cross sections.

ORPHEUS implements a power iteration Monte Carlo solver for the \(k\)-eigenvalue problem in a 2-D unit cell with periodic boundary conditions. Unlike the deterministic The Power Iteration Algorithm, MC power iteration uses stochastic sampling of neutron histories rather than matrix operations on discretised flux vectors. The key algorithmic features are:

• Woodcock delta-tracking — virtual-collision rejection that avoids distance-to-surface calculations, enabling simple geometry interfaces.

• Analog absorption with fission weight adjustment — on absorption, the neutron weight is multiplied by \(\nu\Sigma_f / \Sigma_a\) and reborn from the fission spectrum.

• Russian roulette and splitting — population control to maintain a stable neutron count across power iteration cycles.

• Batch statistics — the \(k\)-effective is estimated from the cycle-to-cycle weight ratio with Central Limit Theorem uncertainty.

Derivation sources. The analytical eigenvalues used for verification are computed independently by the derivation scripts. These are the source of truth for all equations in this chapter:

• derivations/mc.py — homogeneous eigenvalues from random walk probability (kinf_homogeneous()), heterogeneous eigenvalues from CP cylinder reference (kinf_from_cp()).

• derivations/_xs_library.py — synthetic cross-section library (regions A/B/C/D in {1G, 2G, 4G}).

Every numerical value cited in this chapter was produced by these scripts.

Solver Architecture

The MC solver is structured as a modular pipeline with five layers (MT-20260406-008), following the same separation-of-concerns pattern as the CP (CPMesh → CPSolver → solve_cp()) and SN (SNMesh → SNSolver → solve_sn()) solvers.

Mesh1D (edges, mat_ids, coord)          ← base geometry
    |
    v
MCMesh (point-wise material_id_at)      ← augmented geometry
    |
    v
Particle → Neutron                      ← per-particle state
    |
    v
NeutronBank                             ← population management
    |
    v
solve_monte_carlo()                     ← orchestrator
├── _precompute_xs()                       XS cache
├── for cycle:
│   ├── _random_walk()                     transport kernel
│   ├── _russian_roulette()                population control
│   └── _split_heavy()                     population control
└── MCResult

Layer 1: Geometry

The MC solver follows the same three-layer geometry pattern as CP and SN:

1. Base geometry — Mesh1D stores cell edges, material IDs, and the coordinate system.

2. Augmented geometry — MCMesh wraps a Mesh1D and adds the MC-specific point-wise material lookup needed by delta-tracking.

   • Cartesian — x position mapped to cell index via np.searchsorted(edges, x, side="right") - 1.

   • Cylindrical — radial distance \(r = \sqrt{(x - x_c)^2 + (y - y_c)^2}\) from the cell centre, then mapped to cell index.

3. Solver — solve_monte_carlo() receives an MCGeometry protocol object (satisfied by MCMesh, ConcentricPinCell, or SlabPinCell).

Design rationale (MT-20260406-001). Delta-tracking only needs the material at the collision point — no distance-to-surface calculation. This makes the geometry interface minimal:

@runtime_checkable
class MCGeometry(Protocol):
    pitch: float
    def material_id_at(self, x: float, y: float) -> int: ...

The MCGeometry protocol is runtime-checkable, enabling both MCMesh (wrapping Mesh1D) and direct geometry classes (ConcentricPinCell, SlabPinCell) to be used interchangeably. Verified by test_mc_properties.py::test_mcmesh_satisfies_protocol.

MCMesh vs standalone geometry classes. The standalone classes predate the unified geometry module. MCMesh integrates with the same Mesh1D factories used by CP and SN (e.g., pwr_pin_equivalent()), enabling cross-method comparisons on identical base geometry. Verified: test_mc_properties.py::test_mcmesh_cylindrical_matches_concentric shows zero mismatches over 10,000 random sample points.

Layer 2: Particle Hierarchy

Each transported particle carries its own phase-space state. The two-level hierarchy separates particle-type-independent state (position, weight) from neutron-specific state (energy group):

@dataclass
class Particle:
    """Base class for transported particles."""
    x: float
    y: float
    weight: float
    alive: bool = True

@dataclass
class Neutron(Particle):
    """A neutron with discrete energy group."""
    group: int = 0

Extensibility. Future gamma-ray transport would subclass Particle without group (photons may use continuous energy or a different discretization). The NeutronBank and population control functions operate on the particle arrays, so generalizing to a ParticleBank would require minimal changes.

Layer 3: NeutronBank

NeutronBank encapsulates the neutron population as parallel arrays (x, y, weight, group) with an integer n tracking the number of live particles. This design balances two goals:

1. Clean API — methods like normalize_weights(), compact(), grow() replace inline array manipulation scattered throughout the solver.

2. Performance — the inner transport loop accesses bank.x[i_n], bank.weight[i_n] directly (same speed as raw array indexing). Per-particle object creation is avoided in the hot path.

class NeutronBank:
    x: np.ndarray           # positions
    y: np.ndarray
    weight: np.ndarray      # statistical weights
    group: np.ndarray       # energy group indices
    n: int                  # number of live neutrons

    @classmethod
    def initialize(cls, n_neutrons, pitch, chi_cum, ng, rng): ...
    def normalize_weights(self, target): ...
    def save_start_weights(self) -> np.ndarray: ...
    def compact(self): ...           # remove dead (w=0) neutrons
    def grow(self, n_extra): ...     # extend arrays for splitting

Layer 4: Cross-Section Cache

_PrecomputedXS groups data that is constant across cycles:

@dataclass
class _PrecomputedXS:
    sig_t_max: np.ndarray       # (ng,) majorant per group
    sig_s_dense: dict           # mat_id -> (ng, ng) dense P0 scattering
    chi_cum: np.ndarray         # cumulative fission spectrum
    ng: int
    eg: np.ndarray              # energy group boundaries
    eg_mid: np.ndarray          # mid-group energies
    du: np.ndarray              # lethargy widths

Built by _precompute_xs(), which:

• Computes the majorant \(\Sigma_{\text{maj},g} = \max_m \Sigma_{t,m,g}\) (verified by test_majorant_computation).

• Converts sparse SigS[0] to dense arrays (P0 isotropic scattering only).

• Builds the cumulative fission spectrum for group sampling.

Layer 5: Orchestrator

solve_monte_carlo() is a ~60-line orchestrator that wires the layers together:

def solve_monte_carlo(materials, params) -> MCResult:
    xs = _precompute_xs(materials)
    bank = NeutronBank.initialize(...)

    for i_cycle in range(total_cycles):
        bank.normalize_weights(params.n_neutrons)
        w0 = bank.save_start_weights()

        for i_n in range(bank.n):
            _random_walk(bank, i_n, geom, materials, xs, rng, tally)

        _russian_roulette(bank, w0, rng)
        bank.compact()
        _split_heavy(bank, rng)
        keff_cycle = bank.weight[:bank.n].sum() / w0.sum()
        # ... accumulate statistics ...

    return MCResult(...)

Comparison with deterministic solvers. The CP and SN solvers use a shared power_iteration() loop that calls the EigenvalueSolver protocol (five methods: initial_flux_distribution, compute_fission_source, solve_fixed_source, compute_keff, converged). The MC solver cannot satisfy this protocol because \(k_{\text{eff}}\) emerges stochastically from cycle-to-cycle weight ratios, not from a deterministic flux/source iteration. Instead, the MC orchestrator implements its own cycle loop with population control (roulette, splitting) that has no deterministic analogue.

Geometry Mismatch: Square Cell vs Wigner-Seitz Cylinder (ERR-017)

When comparing MC (square cell, periodic BCs) against CP (Wigner-Seitz cylinder, white BCs), the cell areas must match:

(1)\[p^2 = \pi R_{\text{cell}}^2 \quad \Longrightarrow \quad p = R_{\text{cell}} \sqrt{\pi}\]

This is the convention used by pwr_pin_equivalent() (r_cell = pitch / sqrt(pi)).

ERR-017 investigation history. The pre-existing heterogeneous tests used the formula pitch = r_cell * sqrt(pi) * 2, which gives four times the correct cell area. The extra area was all moderator, which:

• For 1G 2-region: \(k_{\text{MC}} = 0.757\) vs \(k_{\text{ref}} = 0.990\) (24% systematic error).

• For 2G 2-region: the neutron population collapsed to zero (NaN) because the subcritical system with 4× moderator could not sustain 200 neutrons/cycle.

How it hid. All homogeneous tests passed (single material — pitch is irrelevant). The heterogeneous tests were marked @pytest.mark.slow and may not have been run regularly. The NaN z-scores failed the < 5.0 assertion without indicating the error direction.

Fix. Corrected to pitch = r_cell * sqrt(pi). Even with matched areas, the square corners contain extra moderator not present in the circle, introducing a ~3–6% systematic bias. The heterogeneous test tolerance accounts for this:

(2)\[|k_{\text{MC}} - k_{\text{ref}}| < 5\sigma_{\text{MC}} + 0.06 \, k_{\text{ref}}\]

Lesson (meta-lesson 8): When constructing geometry for cross-method comparison, verify the cell area/volume matches between the two methods. A factor-of-2 in a linear dimension is a factor-of-4 in area — large enough to change supercritical to subcritical.

The Monte Carlo Random Walk

The random walk is the heart of the MC solver. Each neutron undergoes a sequence of free flights and collisions until it is absorbed. The implementation in solve_monte_carlo() (lines 269–328 of monte_carlo.py) is annotated below.

Free-Flight Distance (Delta-Tracking)

Standard Monte Carlo tracks a neutron to the nearest material boundary, then samples a collision within the current region. Woodcock delta-tracking [Woodcock1965] eliminates the distance-to-surface calculation by introducing a fictitious majorant cross section:

(3)\[\Sigma_{\text{maj},g} = \max_m \Sigma_{t,m,g}\]

where the maximum is over all materials \(m\) for energy group \(g\). In the code:

sig_t_max = np.zeros(ng)
for mix in materials.values():
    sig_t_max = np.maximum(sig_t_max, mix.SigT)

Verified by test_mc_properties.py::test_majorant_computation: for 2G materials A (\(\Sigma_t = [0.50, 1.00]\)) and B (\(\Sigma_t = [0.60, 2.00]\)), the majorant is \([0.60, 2.00]\).

The free-flight distance is sampled from an exponential distribution with rate \(\Sigma_{\text{maj},g}\):

(4)\[s = -\frac{\ln \xi}{\Sigma_{\text{maj},g}}, \qquad \xi \sim U(0,1)\]

This gives \(E[s] = 1/\Sigma_{\text{maj},g}\) and \(\text{Var}[s] = 1/\Sigma_{\text{maj},g}^2\). Both moments are verified by test_mc_gaps.py::test_free_path_exponential.

In the code:

free_path = -np.log(rng.random()) / sig_t_max[ig]

Delta-Tracking Equivalence Proof

Claim: Sampling free flights from the majorant and rejecting virtual collisions recovers the correct transport kernel.

Proof. Decompose the total cross section at each point as real + virtual:

(5)\[\Sigma_{\text{maj},g} = \Sigma_{t,g}(\mathbf{r}) + \Sigma_{\text{virtual},g}(\mathbf{r})\]

The probability of a neutron at \(\mathbf{r}_0\) reaching distance \(s\) without any (real or virtual) collision is:

\[P(\text{no collision up to } s) = e^{-\Sigma_{\text{maj},g} \, s}\]

At each collision site \(\mathbf{r}_0 + s\hat{\Omega}\), the collision is real with probability:

\[P(\text{real} | \text{collision at } s) = \frac{\Sigma_{t,g}(\mathbf{r}_0 + s\hat{\Omega})}{\Sigma_{\text{maj},g}}\]

The joint probability of first real collision at distance \(s\) is:

\[p(s) \, ds = \Sigma_{\text{maj},g} \, e^{-\Sigma_{\text{maj},g} \, s} \cdot \frac{\Sigma_{t,g}(\mathbf{r}(s))}{\Sigma_{\text{maj},g}} \, ds = \Sigma_{t,g}(\mathbf{r}(s)) \, e^{-\Sigma_{\text{maj},g} \, s} \, ds\]

For a homogeneous medium (\(\Sigma_{t,g}\) constant along the path), this simplifies to the familiar exponential attenuation \(\Sigma_t e^{-\Sigma_t s}\) since the virtual collisions (which preserve direction and continue the flight) effectively reduce the rate from \(\Sigma_{\text{maj}}\) to \(\Sigma_t\).

Key property: Virtual collisions preserve the flight direction. This is essential — if direction were resampled on virtual collisions, the path would become a random walk instead of a straight line, and the exponential attenuation would be violated. Verified by test_mc_gaps.py::test_virtual_collision_preserves_direction.

In the code:

# Virtual or real collision?
if sig_v / sig_t_max[ig] >= rng.random():
    virtual_collision = True    # direction unchanged
else:
    virtual_collision = False   # process real collision

The virtual collision probability at the collision site is:

\[P_{\text{virtual}} = \frac{\Sigma_{\text{maj},g} - \Sigma_{t,g}} {\Sigma_{\text{maj},g}}\]

Verified by test_mc_properties.py::test_delta_tracking_virtual_probability and test_mc_properties.py::test_delta_tracking_homogeneous_no_virtual (homogeneous medium → zero virtual collisions).

Efficiency. The virtual collision fraction in material \(m\) is \(1 - \Sigma_{t,m,g}/\Sigma_{\text{maj},g}\). For the 2G library with materials A and B: in the fast group, B has \(\Sigma_t = 0.60 = \Sigma_{\text{maj}}\) (0% virtual), while A has \(\Sigma_t = 0.50\) (17% virtual). In the thermal group, A has \(\Sigma_t = 1.00\) (50% virtual). High virtual rates waste computation; see MT-20260403-006.

Direction Sampling

At each real collision that results in scattering, a new flight direction is sampled. The solver uses a 2-D projection:

(6)\[\begin{split}\theta &= \pi \, \xi_1, \qquad \phi = 2\pi \, \xi_2 \\ \Omega_x &= \sin\theta \, \cos\phi, \qquad \Omega_y = \sin\theta \, \sin\phi\end{split}\]

where \(\xi_1, \xi_2 \sim U(0,1)\). In the code:

theta = np.pi * rng.random()
phi = 2.0 * np.pi * rng.random()
dir_x = np.sin(theta) * np.cos(phi)
dir_y = np.sin(theta) * np.sin(phi)

Direction Sampling Limitation (ERR-018)

This is NOT isotropic sampling on the unit sphere. True isotropic sampling requires \(\theta = \arccos(1 - 2\xi_1)\) to account for the solid-angle Jacobian \(\sin\theta \, d\theta \, d\phi\).

Mathematical analysis. For uniform \(\theta\):

\[E[\sin^2\theta] = \frac{1}{\pi} \int_0^{\pi} \sin^2\theta \, d\theta = \frac{1}{2}\]

For true isotropic sampling:

\[E[\sin^2\theta]_{\text{iso}} = \frac{1}{2} \int_0^{\pi} \sin^3\theta \, d\theta = \frac{2}{3}\]

Therefore:

\[E[\Omega_x^2] = E[\sin^2\theta] \cdot E[\cos^2\phi] = \frac{1}{2} \cdot \frac{1}{2} = \frac{1}{4} \qquad\text{(uniform } \theta\text{)}\]

vs \(1/3\) for isotropic. The uniform-\(\theta\) formula shortens the average 2-D step length by \(\sqrt{1/4} / \sqrt{1/3} = \sqrt{3/4} \approx 0.87\), i.e., ~13% shorter.

Why it’s retained. The formula matches the original MATLAB monteCarloPWR.m. Since the solver tracks only the 2-D projection \((\Omega_x, \Omega_y)\) and uses periodic BCs on a square cell, the non-isotropic sampling changes the effective mean free path but does not invalidate the eigenvalue — the bias is absorbed into the geometry scaling. Fixing would break MATLAB compatibility and require re-establishing all reference values (MT-20260406-006).

Verified by test_mc_properties.py::test_direction_sampling: \(E[\Omega_x^2] = 0.25 \pm z < 5\).

Periodic Boundary Conditions

The unit cell is a square of side length \(p\) (pitch). When a neutron exits the cell, its position wraps via the modulo operation:

(7)\[x \leftarrow x \bmod p, \qquad y \leftarrow y \bmod p\]

In the code:

nx_ = nx_ % pitch
ny_ = ny_ % pitch

Python’s modulo returns non-negative for positive pitch. Verified by test_mc_properties.py::test_periodic_bc_wrapping with edge cases including negative positions and exact-boundary positions.

Collision Physics

At each real collision site, two outcomes are possible: scattering or absorption. The decision is based on the branching ratio:

(8)\[P(\text{scatter}) = \frac{\Sigma_{s,g}}{\Sigma_{t,g}}, \qquad P(\text{absorb}) = \frac{\Sigma_{a,g}}{\Sigma_{t,g}}\]

where \(\Sigma_{a,g} = \Sigma_{f,g} + \Sigma_{c,g} + \Sigma_{L,g}\). In the code:

sig_a = mat.SigF[ig] + mat.SigC[ig] + mat.SigL[ig]
sig_s_row = sig_s_dense[mat_id][ig, :]
sig_s_sum = sig_s_row.sum()
sig_t = sig_a + sig_s_sum

Verified by test_mc_properties.py::test_scattering_branching_ratio: for region A (1G), \(P(\text{scatter}) = 0.5/1.0 = 0.5\), confirmed by 100k samples. XS consistency (\(\Sigma_t = \Sigma_a + \Sigma_s\)) verified for all 12 materials by test_mc_gaps.py::test_xs_consistency_in_solver.

Cross-section preprocessing. The solver precomputes dense scattering rows for all materials:

sig_s_dense = {}
for mat_id, mix in materials.items():
    rows = np.array(mix.SigS[0].todense())  # (ng, ng)
    sig_s_dense[mat_id] = rows

The SigS[0] index selects the isotropic (P0) scattering matrix. Higher-order anisotropy (P1, etc.) is not used by the MC solver. The todense() conversion is necessary because the source storage is sparse (scipy.sparse.csr_matrix).

Scattering: Group Transfer

On scattering, the outgoing energy group \(g'\) is sampled from the cumulative distribution of the scattering row:

(9)\[P(g' \le G \mid \text{scatter from } g) = \frac{\sum_{g'=0}^{G} \Sigma_{s,g \to g'}} {\sum_{g'} \Sigma_{s,g \to g'}}\]

In the code:

cum_s = np.cumsum(sig_s_row)
ig = np.searchsorted(cum_s, rng.random() * sig_s_sum)
ig = min(ig, ng - 1)

Note

Convention (anti-ERR-002): The code uses sig_s_dense[mat_id][ig, :] which is row ig of the scattering matrix. The codebase convention is \(\Sigma_s[\text{from}, \text{to}]\), so row ig gives transfer from group ig to all groups.

For an asymmetric 2G matrix with zero upscatter (\(\Sigma_s[1,0] = 0\)), a neutron in group 1 can only scatter to group 1. If the code used the column (transpose), upscatter would appear. Verified by test_mc_properties.py::test_scattering_convention_no_upscatter: 10,000 samples from group 1, zero upscatter events.

Scattering fractions verified by test_mc_properties.py::test_scattering_cdf_sampling: for region A (2G), the expected fast-to-fast fraction is \(0.38/0.48 = 79.2\%\), confirmed by 100k samples within \(z < 5\).

Absorption: Fission Weight Adjustment

On absorption, rather than killing the neutron and sampling a new fission site, the solver uses analog absorption with weight adjustment. The neutron’s weight is multiplied by the expected number of fission neutrons per absorption:

(10)\[w \leftarrow w \cdot \frac{\nu\Sigma_{f,g}}{\Sigma_{a,g}}\]

In the code:

sig_p = mat.SigP[ig]           # = nu * Sigma_f
sig_a = mat.SigF[ig] + mat.SigC[ig] + mat.SigL[ig]
w *= sig_p / sig_a

1-group derivation from random walk probability. Consider a neutron undergoing a random walk in a homogeneous 1G infinite medium. At each collision:

1. The neutron scatters with probability \(c = \Sigma_s / \Sigma_t\).

2. After \(n\) scattering events, it is absorbed with probability \((1-c)\).

3. On absorption, it produces \(\nu\Sigma_f / \Sigma_a\) fission neutrons.

The probability of reaching the \(n\)-th collision is \(c^n\). The expected multiplication per generation is:

\[k = (1-c) \cdot \frac{\nu\Sigma_f}{\Sigma_a} \cdot \sum_{n=0}^{\infty} c^n = (1-c) \cdot \frac{\nu\Sigma_f}{\Sigma_a} \cdot \frac{1}{1-c} = \frac{\nu\Sigma_f}{\Sigma_a}\]

This is exact for 1G and equals the weight factor applied at each absorption. For region A (1G): \(\Sigma_a = 0.2 + 0.3 = 0.5\), \(\nu\Sigma_f = 2.5 \times 0.3 = 0.75\), so \(k = 1.5\). The numerical value is produced by derivations/mc.py via kinf_homogeneous(); the derivation above is presented for pedagogical context.

Verified by test_mc_properties.py::test_fission_weight_adjustment (hand-calculated weight factor = 1.5) and test_mc_properties.py::test_1g_homogeneous_deterministic (\(\sigma = 0\) because every neutron sees identical XS).

Non-fissile materials. When \(\nu\Sigma_f = 0\) (materials B, C, D), the weight goes to zero: \(w \leftarrow 0\). Russian roulette then terminates the neutron. Verified by test_mc_properties.py::test_fission_weight_non_fissile and test_mc_gaps.py::test_absorption_nonfissile_zeroes_weight.

Fission spectrum resampling. After weight adjustment, the neutron is reborn in a new energy group from the fission spectrum CDF:

(11)\[g_{\text{new}} = \text{searchsorted}\!\bigl( \text{cumsum}(\chi), \, \xi \bigr), \qquad \xi \sim U(0,1)\]

For 4G region A: \(\chi = [0.60, 0.35, 0.05, 0.00]\). Verified by test_mc_properties.py::test_chi_spectrum_sampling: 100k samples match expected group fractions within \(z < 5\) per group. Group 3 gets exactly zero samples.

Limitation: The fission spectrum \(\chi\) is taken from a single material (_any_mat), not the material at the absorption site. For heterogeneous problems with different \(\chi\) vectors, this is a simplification. All test materials use the same \(\chi\).

Population Control

Russian Roulette

After a complete random walk (free flights until absorption), neutrons with reduced weight are subjected to Russian roulette.

(12)\[P_{\text{kill}} = 1 - \frac{w}{w_0}\]

where \(w\) is the post-walk weight and \(w_0\) is the weight at cycle start. The outcome is:

(13)\[\begin{split}w_{\text{after}} = \begin{cases} 0 & \text{with probability } P_{\text{kill}} \\ w_0 & \text{with probability } 1 - P_{\text{kill}} \end{cases}\end{split}\]

In the code:

terminate_p = 1.0 - weight[i_n] / weight0[i_n]
if terminate_p >= rng.random():
    weight[i_n] = 0.0
elif terminate_p > 0:
    weight[i_n] = weight0[i_n]

Weight Conservation Proof

The expected weight after roulette equals the weight before:

(14)\[E[w_{\text{after}}] = (1 - P_{\text{kill}}) \cdot w_0 + P_{\text{kill}} \cdot 0 = \frac{w}{w_0} \cdot w_0 = w = w_{\text{before}}\]

Verified statistically by test_mc_properties.py::test_roulette_weight_conservation (100k neutrons with \(w/w_0 = 0.3\), mean weight after roulette matches \(0.3\) within \(z < 5\)) and by test_mc_properties.py::test_roulette_restore_weight (surviving neutrons have weight exactly \(w_0\)).

Supercritical edge case. When \(w > w_0\) (fission weight adjustment in supercritical system, e.g., \(\nu\Sigma_f/\Sigma_a = 1.5\)), then \(P_{\text{kill}} = 1 - w/w_0 < 0\). The terminate_p >= rng.random() condition is never true (rng returns \([0,1)\)), and terminate_p > 0 is also false, so the weight is unchanged at \(w\). This is correct — supercritical neutrons should not be rouletted. Verified by test_mc_gaps.py::test_roulette_supercritical_preserves_weight.

Splitting

Neutrons with weight \(w > 1\) are split into \(N\) copies, each with weight \(w/N\). The number of copies is:

(15)\[\begin{split}N = \begin{cases} \lfloor w \rfloor + 1 & \text{with probability } w - \lfloor w \rfloor \\ \lfloor w \rfloor & \text{with probability } \lfloor w \rfloor + 1 - w \end{cases}\end{split}\]

Splitting Weight Conservation Proof

Exact conservation: Each copy gets weight \(w/N\), so total weight = \(N \cdot (w/N) = w\). No stochastic noise.

Expected copies: The stochastic rounding gives:

\[E[N] = \lfloor w \rfloor \cdot (1 - (w - \lfloor w \rfloor)) + (\lfloor w \rfloor + 1) \cdot (w - \lfloor w \rfloor)\]

Let \(f = w - \lfloor w \rfloor\) and \(n = \lfloor w \rfloor\):

\[E[N] = n(1-f) + (n+1)f = n + f = w\]

So \(E[N] = w\), and total weight per copy is \(w/N\), meaning the expected number of particles times weight per particle equals \(w\).

Verified by test_mc_properties.py::test_splitting_weight_conservation (exact total weight) and test_mc_gaps.py::test_splitting_copy_count (statistical: \(P(N=4) = 0.7\) for \(w = 3.7\)).

Eigenvalue Estimation

keff from Weight Ratio

The \(k\)-effective for each cycle is estimated as:

(16)\[k_{\text{cycle}} = \frac{\sum_{n} w_n^{\text{end}}}{\sum_{n} w_n^{0}}\]

where \(w_n^{0}\) are the normalised starting weights and \(w_n^{\text{end}}\) are the weights after the random walk, roulette, and splitting. In the code:

keff_cycle = weight[:n_neutrons].sum() / weight0.sum()

Verified by test_mc_gaps.py::test_keff_cycle_estimator: \([1.5, 0.0, 1.2, 0.8, 2.0] / [1,1,1,1,1] = 5.5/5.0 = 1.1\).

Weight normalisation. At the start of each cycle, weights are rescaled so that \(\sum w_n = N_{\text{source}}\):

total_weight = weight[:n_neutrons].sum()
weight[:n_neutrons] *= params.n_neutrons / total_weight

Verified by test_mc_gaps.py::test_weight_normalization_consistency.

Batch Statistics

The final \(k\)-effective is the cumulative mean of \(M\) active cycle values:

(17)\[\bar{k}_M = \frac{1}{M} \sum_{m=1}^{M} k_m\]

The standard deviation of the mean:

(18)\[\sigma_M = \sqrt{\frac{1}{M(M-1)} \sum_{m=1}^{M} (k_m - \bar{k}_M)^2}\]

By the Central Limit Theorem, \(\sigma_M \sim 1/\sqrt{M}\). Verified by test_mc_convergence.py::test_sigma_scales_with_sqrt_n: the ratio \(\sigma(400)/\sigma(100)\) is within [0.25, 0.75] of the theoretical 0.5.

Verified term-by-term by test_mc_properties.py::test_batch_statistics_formula: for \(k_{\text{active}} = [1.0, 1.1, 0.9, 1.05]\), the hand-calculated cumulative means and sigmas match the solver’s formula to machine precision.

Source convergence. The first \(N_{\text{inactive}}\) cycles are discarded to allow the fission source to converge from the initial (uniform) guess. Only \(N_{\text{active}}\) cycles contribute. Verified by test_mc_convergence.py::test_inactive_cycles_reduce_bias.

Flux Tally (Known Limitation)

The solver accumulates a scattering detector during the random walk:

detect_s[ig] += w / sig_s_sum

at each scattering event in group \(g\) (line 307).

Warning

This is not a proper flux estimator (MT-20260406-007). A correct collision estimator would accumulate \(w / \Sigma_{t,g}\) at every real collision (both scattering and absorption). The current tally:

• Misses absorption events entirely

• Weights by \(1/\Sigma_s\) instead of \(1/\Sigma_t\)

• Produces negative flux_per_lethargy values because the energy grid is high-to-low (\(\Delta u = \ln(E_{g+1}/E_g) < 0\))

The eigenvalue \(k_{\text{eff}}\) is computed from the weight ratio (16) and is not affected by this tally.

Solver Parameters and Results

MCParams fields

Field	Default	Description
n_neutrons	100	Source neutrons per cycle
n_inactive	100	Inactive cycles (source convergence, not tallied)
n_active	2000	Active cycles (tallied for \(k_{\text{eff}}\))
pitch	3.6	Unit cell side length (cm)
seed	None	RNG seed (None = random)
geometry	None	MCGeometry (None → default slab)

MCResult fields

Field	Shape	Description
keff	scalar	Final cumulative mean \(\bar{k}_M\)
sigma	scalar	Standard deviation of the mean \(\sigma_M\)
keff_history	(n_active,)	Cumulative mean at each active cycle
sigma_history	(n_active,)	Cumulative sigma at each active cycle
flux_per_lethargy	(ng,)	Scattering detector / \(\Delta u\) (see limitation above)
eg_mid	(ng,)	Mid-group energies (eV)
elapsed_seconds	scalar	Wall-clock time

Analytical Verification

Homogeneous Infinite Medium

For a single material filling the entire cell, the eigenvalue is derived from the collision kernel.

1-group: From the random walk analysis above ((10)):

(19)\[k_\infty = \frac{\nu\Sigma_f}{\Sigma_a}\]

For region A (1G): \(k = 0.75/0.5 = 1.5\).

Multi-group: The eigenvalue is the dominant eigenvalue of

(20)\[k_\infty = \lambda_{\max}\!\left(\mathbf{A}^{-1} \mathbf{F}\right)\]

where \(\mathbf{A} = \text{diag}(\Sigma_t) - \Sigma_s^T\) (the transpose converts from the [from, to] storage convention to the in-scatter form needed by the balance equation; see scattering convention note in Monte Carlo Neutron Transport) and \(\mathbf{F} = \chi \otimes (\nu\Sigma_f)\). Computed by kinf_homogeneous().

Homogeneous eigenvalues

Case	\(k_{\text{ref}}\)	Tolerance
1G (mc_cyl1D_1eg_1rg)	1.500000	\(\sigma = 0\) (exact)
2G (mc_cyl1D_2eg_1rg)	1.875000	\(z < 5\)
4G (mc_cyl1D_4eg_1rg)	1.487762	\(z < 5\)

Heterogeneous Pin Cell

For multi-region problems, the reference eigenvalue comes from the CP cylinder derivation (kinf_from_cp()).

MC heterogeneous verification cases

Case	G	Reg	\(k_{\text{ref}}\)	Tolerance
mc_cyl1D_1eg_2rg	1	2	0.989750	\(5\sigma + 0.06 k\)
mc_cyl1D_2eg_2rg	2	2	0.739887	\(5\sigma + 0.06 k\)
mc_cyl1D_1eg_4rg	1	4	0.806751	\(5\sigma + 0.06 k\)
mc_cyl1D_2eg_4rg	2	4	0.502725	\(5\sigma + 0.06 k\)
mc_cyl1D_4eg_2rg	4	2	0.648637	\(5\sigma + 0.06 k\)
mc_cyl1D_4eg_4rg	4	4	0.446472	\(5\sigma + 0.06 k\)

Verification Suite

The MC solver is verified by 55 tests across four levels:

Verification test summary

Level	Count	Description
L0	31	Term-level isolation of each algorithmic component
L1	18	Eigenvalue: {1,2,4}G × {1,2,4}-region
L2	3	Convergence: \(\sigma \sim 1/\sqrt{N}\), bias, inactive cycles
XV	2	Cross-verification: MC vs CP (cylinder, slab)

Determinism. Verified by test_mc_gaps.py::test_seed_reproducibility (same seed → identical results) and test_different_seeds_differ (different seeds → different histories).

Failure Mode Coverage

Failure mode	Tests
Sign flip	L0-MC-004 (delta-tracking), L0-MC-016 (collision fraction)
Variable swap	L0-MC-005, L0-MC-014 (SigS^T)
Missing factor	L0-MC-006 (fission weight), L0-MC-015 (free-path)
Factor error	L0-MC-012 (direction moments)
Index error	L0-MC-005, L0-MC-007 (searchsorted)
Convention drift	L0-MC-013 (batch stats), XS consistency
1G degeneracy	L1: 2G/4G high-stats, 4G fast guard
Homogeneous degeneracy	L1: heterogeneous {2,4}-region
Geometry mismatch	ERR-017 fixed, pitch tested

Design Decisions

Delta-tracking vs surface-tracking (MT-20260403-001). Delta-tracking eliminates distance-to-surface calculations, simplifying the geometry interface to material_id_at(x, y). The cost is virtual collisions in low-density regions.

Analog absorption vs implicit capture (MT-20260403-002). Implicit capture reduces weight at every collision (\(w \leftarrow w \cdot \Sigma_s/\Sigma_t\)). Analog absorption terminates the walk at each absorption event, producing exactly one surviving particle with weight \(w \cdot \nu\Sigma_f/\Sigma_a\). Combined with roulette and splitting, this maintains a stable population.

Uniform theta sampling (ERR-018, MT-20260406-006). Retained from MATLAB original. See Direction Sampling Limitation (ERR-018).

1-group is degenerate. \(k = \nu\Sigma_f/\Sigma_a\) regardless of code correctness. The suite tests at 1, 2, AND 4 groups.

Particle/Neutron/Bank decomposition (MT-20260406-008). The solver was originally a single monolithic function (~220 lines). The restructuring separates concerns into five layers (see Monte Carlo Neutron Transport architecture section) following the same pattern as the CP and SN solvers. The Particle → Neutron hierarchy enables future gamma-ray transport without modifying the bank or population control code. The NeutronBank uses parallel arrays (not per-particle objects) in the hot path for performance parity with the original implementation. Verified: test_seed_reproducibility proves identical RNG sequences before and after restructuring.

Array-backed bank vs object-per-neutron. Each neutron’s random walk is inherently serial (next collision depends on current outcome), so per-particle Python overhead matters. The bank stores parallel arrays and the inner loop indexes them directly (bank.x[i_n]), giving the same performance as raw arrays. The Particle/Neutron dataclasses serve as the conceptual model and API documentation, not as the runtime representation in the hot path.

Limitations and Future Work

Tracking ID	Description
MT-20260403-004	Python neutron loop performance (inner loop not vectorised)
MT-20260403-005	Heterogeneous independent reference (high-statistics MC)
MT-20260403-006	Per-region majorant for efficiency
MT-20260406-005	Solver ignores Sig2 (n,2n) reactions
MT-20260406-006	Direction sampling not isotropic (ERR-018)
MT-20260406-007	Flux tally is not a proper estimator
MT-20260406-009	Event-based vectorised transport (process all neutrons per collision step using numpy batch operations; expected 10–50× speedup)
MT-20260406-010	Parallel neutron transport (multiprocessing with per-worker RNG via SeedSequence.spawn)

References

[Woodcock1965]

E.R. Woodcock, T. Murphy, P.J. Hemmings, and T.C. Longworth, “Techniques used in the GEM code for Monte Carlo neutronics calculations in reactors and other systems of complex geometry,” Proc. Conf. Applications of Computing Methods to Reactor Problems, ANL-7050, 1965.

[Lux1991]

I. Lux and L. Koblinger, Monte Carlo Particle Transport Methods: Neutron and Photon Calculations, CRC Press, 1991.

[Brown2005]

F.B. Brown, “Fundamentals of Monte Carlo Particle Transport,” LA-UR-05-4983, Los Alamos National Laboratory, 2005.