Verification of a Hill–Wheeler Statistical Fission Model Using ENDF/B-VIII.0 Data and a Reproducible Python Implementation for Charge-Yield Calculations

Verification of a Hill–Wheeler Statistical Fission Model Using ENDF/B-VIII.0 Data and a Reproducible Python Implementation for Charge-Yield Calculations H. Maruyama1 1Independent Researcher December 16, 2025 Abstract Fission-product charge-yield distributions are key inputs for reactor physics, fuelcycle analysis, and nuclear-data evaluation. However, accurate theoretical prediction remains challenging. In this work, we re-implement a Selective Channel Scission (SCS) statistical fission model augmented by the Hill–Wheeler penetrability and interpret the Hill–Wheeler logistic form as a quantum-statistical distribution analogous to the Fermi–Dirac function. We provide a fully reproducible Python workflow and validate it against evaluated charge-yield data from ENDF/B-VIII.0. We analyze nine actinides (232Th, 233U, 235U, 238U, 237Np, 239Pu, 240Pu, 242Pu, and 241Am) under three incident-neutron conditions (thermal, 500 keV, and 14 MeV, where applicable). The calculated charge yields Y(Z)reproduce the evaluated distributions with good agreement across the considered cases. We further extract an effective fission distance dw(Z)and a channel-dependent “Fermi energy” correction Ex(Z)from two-stage fitting, and we summarize the maximum Ex(Z)values at low energies for comparison with the fission-barrier scale. 1 Keywords: nuclear fission; Hill–Wheeler formula; Selective Channel Scission (SCS) model; charge yields; ENDF/B-VIII.0; Python; statistical model; Fermi level 1 Introduction Fission-product yield distributions (mass and charge yields) directly affect key quantities in reactor-core calculations (decay heat, activation, transmutation) and in spent-fuel composition evaluations. Although microscopic approaches continue to progress, firstprinciples prediction of detailed yield distributions remains difficult in routine applications. Consequently, evaluated nuclear-data libraries and phenomenological/statistical models play an essential role. The Selective Channel Scission (SCS) model [1, 2] treats fission as a set of discrete channels specified by fragment proton and neutron numbers and assigns statistical weights based on barrier-penetration arguments. In our recent work [3], the Hill–Wheeler penetrability [4] was reinterpreted as a quantum-statistical distribution of the same logistic form as the Fermi–Dirac function, enabling an analogy with semiconductor physics. The present work aims to (i) provide a clean, reproducible Python implementation, (ii) validate it on an independent evaluated dataset (ENDF/B-VIII.0), and (iii) document the full workflow (inputs/outputs/figures) to facilitate third-party reproduction. 2 Theoretical framework 2.1 Hill–Wheeler formula and its isomorphism to the Fermi–Dirac distribution For a parabolic barrier, the Hill–Wheeler penetrability is given by [4] P(E) = 1 1 + exp h2π(Ec−E) ℏωi.(1) 2 Here Eis the system energy, Ecis the barrier height, and ℏωcharacterizes the curvature near the barrier top. Equation (1) is algebraically identical to the Fermi–Dirac occupation function f(E) = 1 1 + exp hE−EF kBTi,(2) under the correspondence Ec↔EFand ℏω↔kBT[3]. In this paper, we use this logistic form to construct channel weights within an SCS-type statistical model. 2.2 Selective Channel Scission (SCS) model In the SCS model [1, 2], a fission channel is specified by fragment (Z1, N1)and (Z2, N2), constrained by conservation laws Z1+Z2=Ztotal,(3) N1+N2+νp=Ncompound,(4) where νpis the prompt-neutron multiplicity used as a representative value in the channel bookkeeping. For each channel, we define a logistic (Hill–Wheeler-like) weight p(Z1, N1) = 1 1 + exp hkC dw(Z1)−Q−Ex(Z1)i,(5) where C= 1.44 Z1Z2(MeV·fm) is the Coulomb term, dw(Z)(fm) is an effective scission (fission) distance, Q(MeV) is the Q-value computed from atomic masses, and Ex(Z) (MeV) is a channel-dependent “Fermi energy” correction. The scaling factor kis defined consistently with Ref. [3]. The charge yield is computed as Y(Z1) = PN1p(Z1, N1) PZ′ 1PN′ 1p(Z′ 1, N′ 1)×200%,(6) where the factor 200% corresponds to the total yield of two fragments. 3 2.3 Two-stage parameter estimation and low-yield filtering We estimate model parameters using a two-stage procedure similar to Ref. [3]: 1. Step 1: set Ex(Z)=0and fit dw(Z)independently for each Z(least squares in log space). 2. Step 2: approximate dw(Z)by a quadratic function (for stabilization), fix this dw(Z), and fit Ex(Z). In the low-yield tails of Y(Z), the estimation of Ex(Z)can become numerically unstable, producing non-physical spikes near the edges. To mitigate this, the Python implementation excludes points with Yexp(Z)< Ymin from the residual evaluation when fitting Ex(Z) (default Ymin = 10−3) and optionally applies a nuclide-dependent lower-Zcutoff Zmin so that the symmetric-side cutoff is Ztotal −Zmin. 3 Computational procedure 3.1 Data sources and preprocessing •Charge yields: evaluated independent yields from ENDF/B-VIII.0. •Atomic masses: AME2020 [5]. •Prompt neutron multiplicity: representative values based on ENDF (MT=456). The yields are normalized to match Eq. (6) for comparison and fitting. 3.2 Python implementation and optimization The workflow is implemented in Python 3 (tested with standard scientific packages: numpy, scipy,pandas,matplotlib,pyyaml). Parameter estimation uses scipy.optimize.least_squares (trust-region reflective) with bounds •dw(Z)∈[5,40] fm, 4 •Ex(Z)∈[−20,20] MeV. To avoid overflow in logistic computations, a numerically stable sigmoid (e.g., scipy.special.expit) is recommended. 3.3 Analysis cases and notation Table 1 summarizes the analyzed cases. Here Enis the incident-neutron kinetic energy (MeV), νpis the representative prompt-neutron multiplicity used in the channel bookkeeping, and Snis the neutron separation energy (MeV) of the compound nucleus. As a rough estimate, the initial excitation energy may be referenced as E∗≈En+Sn. Table 1: List of analyzed cases (ENDF/B-VIII.0). The Log RMSE values are shown rounded to 6 decimal places. Nuclide Incident energy En(MeV) νpSn(MeV) Log RMSE Converged 232Th 500 keV 0.52 4.7865 0.000000 Yes 232Th 14 MeV 14 4 4.7865 0.000000 Yes 233U thermal 2.53 ×10−82 6.8455 0.000000 Yes 233U 500 keV 0.52 6.8455 0.000000 Yes 233U 14 MeV 14 4 6.8455 0.000000 Yes 235U thermal 2.53 ×10−82 6.5456 0.000000 Yes 235U 500 keV 0.52 6.5456 0.000000 Yes 235U 14 MeV 14 4 6.5456 0.000000 Yes 238U 500 keV 0.52 4.8065 0.000000 Yes 238U 14 MeV 14 4 4.8065 0.000000 Yes 237Np thermal 2.53 ×10−82 5.4884 0.000000 Yes 237Np 500 keV 0.52 5.4884 0.000000 Yes 237Np 14 MeV 14 4 5.4884 0.000000 Yes 239Pu thermal 2.53 ×10−82 6.5343 0.000000 Yes 239Pu 500 keV 0.52 6.5343 0.000000 Yes 239Pu 14 MeV 14 4 6.5343 0.000000 Yes 240Pu thermal 2.53 ×10−82 5.2416 0.000000 Yes 240Pu 500 keV 0.52 5.2416 0.000000 Yes 240Pu 14 MeV 14 4 5.2416 0.000000 Yes 242Pu thermal 2.53 ×10−82 5.0337 0.000000 Yes 242Pu 500 keV 0.52 5.0337 0.000000 Yes 242Pu 14 MeV 14 4 5.0337 0.000000 Yes 241Am thermal 2.53 ×10−82 5.5377 0.000000 Yes 241Am 500 keV 0.52 5.5377 0.000000 Yes 241Am 14 MeV 14 4 5.5377 0.000000 Yes 5 4 Results 4.1 Analyzed nuclides Table 2 lists the nine nuclides and the representative νpvalues used for the thermal, 500 keV, and 14 MeV conditions. Thermal-neutron fission is not defined for 232Th and 238U, and thus thermal cases are not included for these nuclides. Table 2: Analyzed nuclides and representative prompt-neutron multiplicities used in this work. Nuclide Z Atarget thermal 500 keV 14 MeV (νp) (νp) (νp) 232Th 90 232 – 2 4 233U 92 233 2 2 4 235U 92 235 2 2 4 238U 92 238 – 2 4 237Np 93 237 2 2 4 239Pu 94 239 2 2 4 240Pu 94 240 2 2 4 242Pu 94 242 2 2 4 241Am 95 241 2 2 4 4.2 Comparison of charge-yield distributions We compare the calculated charge yields Y(Z)with evaluated ENDF/B-VIII.0 values for the three incident-energy conditions. Figures 1–3 show multi-panel comparisons (symbols: ENDF/B-VIII.0; lines: fitted SCS–HW model). 6 233U235U237Np 239Pu 240Pu 242Pu 241Am Figure 1: Comparison of charge yields Y(Z)at thermal neutron energy (symbols: ENDF/B-VIII.0, lines: SCS–HW fit). 232Th and 238U are excluded because thermalneutron fission is not defined for these nuclides. 7 232Th 233U235U 238U237Np 239Pu 240Pu 242Pu 241Am Figure 2: Comparison of charge yields Y(Z)at incident neutron energy 500 keV (symbols: ENDF/B-VIII.0, lines: SCS–HW fit) for nine nuclides. 8 232Th 233U235U 238U237Np 239Pu 240Pu 242Pu 241Am Figure 3: Comparison of charge yields Y(Z)at incident neutron energy 14 MeV (symbols: ENDF/B-VIII.0, lines: SCS–HW fit) for nine nuclides. 4.3 Estimated fission distance dw(Z) Figure 4 shows the inferred fission-distance parameter as the Step-1 fitted effective distance Reff(Z)(points). Following Eq. (5), the Coulomb term C= 1.44 Z1Z2(MeV·fm) enters as C/dw(Z1), which can be interpreted as an effective Coulomb energy at scission. If this effective Coulomb energy does not vary strongly across charge channels, then one expects the scaling dw(Z)∝Z1Z2=Z(Ztotal −Z),(7) i.e., a single-peaked (approximately parabolic) dependence on Zwith a maximum near symmetric charge split. The inferred Reff(Z)points in Fig. 4 exhibit a clear single9 [9] F. Poortmans, H. Weigmann, and H. Fischer, “Neutron resonance parameters of 242Pu,” Nucl. Phys. A 207, 342–352 (1973). [10] K. Godbey et al., “Isospin composition of fission barriers,” arXiv:2410.07259 (2024). 16