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

Validation and Python Implementation of the Hill–Wheeler Statistical Nuclear Fission Model Using ENDF/B-VIII.0 Data: Reproducible Charge Yield Distribution Calculations H. Maruyama1 1Independent Researcher December 18, 2025 Abstract Fission-product charge-yield distributions are a key foundational dataset and an important input for reactor physics, the nuclear fuel cycle, and nuclear-data evaluation (and data assimilation), yet high-accuracy theoretical prediction remains challenging. In this study, we incorporate the Hill–Wheeler penetrability into the Selective Channel Scission (SCS) framework and adopt a statistical–mechanical interpretation in which the Hill–Wheeler formula is regarded as a “quantum-statistical distribution” isomorphic to the Fermi–Dirac function. The full computational workflow is re-implemented in Python in a reproducible manner. For validation, we use the evaluated nuclear-data library ENDF/B-VIII.0 and compute charge yields Y(Z)for nine nuclides (232Th, 233U, 235U, 238U, 237Np, 239Pu, 240Pu, 242Pu, and 241Am) under incident neutron-energy conditions (thermal, 500 keV, and 14 MeV). The calculated charge-yield distributions show good agreement with the charge-yield data contained in ENDF/B-VIII.0 for all considered conditions. The main contributions of this work are: (i) a transparent Python implementation, (ii) independent validation using ENDF/B-VIII.0, and (iii) a fully reproducible description including figures, tables, and workflow details. Keywords: nuclear fission; Hill–Wheeler formula; Selective Channel Scission (SCS) model; charge yield; ENDF/B-VIII.0; Python; statistical model; Fermi level 1 Introduction Fission-product (FP) yield distributions (mass and charge distributions) directly affect reactor-core calculations, including decay heat, activation, and transmutation; they also play a central role in evaluating nuclide inventories in spent fuel and in nuclear-data assimilation and evaluation. At the same time, treating the many-body quantum dynamics of nuclear fission from first principles remains difficult, and practical applications often rely on evaluated libraries and empirical or semi-empirical models. A classical framework for fission barriers and fission mechanisms was established by Bohr and Wheeler [1]. In recent 1 years, evaluation efforts have advanced not only for yields themselves but also for uncertainties (covariances) [2,3], and further developments in yield modeling and evaluation have been reported [4]. Machine-learning approaches for yield prediction and completion have also been explored [5,6]. Within this context, the SCS model [7,8] provides a statistical approach in which the fission process is decomposed into channels and yields are obtained by weighting each channel by a barrier-penetration probability. Furthermore, Ref. [9] proposed a framework that reinterprets the Hill–Wheeler formula [10] not merely as a penetrability for a parabolic barrier but as a quantum-statistical (logistic) distribution isomorphic to the Fermi–Dirac function, and demonstrated its effectiveness using JENDL-5.0 data. The SCS framework is, in essence, a channel-enumeration-based statistical model; in earlier studies it has been used mainly to discuss yields aggregated over mass number. To the author  s knowledge, an explicit and reproducible treatment and validation of evaluated chargeyield distributions Y(Z)in this statistical–mechanical interpretation has been presented in Ref. [9] and is extended and independently validated here using ENDF/B-VIII.0. The objectives of this work are to (a) reorganize the above framework as a Python implementation, (b) perform an independent validation using ENDF/B-VIII.0 [11], and (c) provide a reproducible description of the workflow, including inputs/outputs, figures, and tables. 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 [10] P(E) = 1 1 + exp [2π(Ec−E) ℏω],(1) where Eis the system energy, Ecis the barrier height, and ℏωcharacterizes the local curvature of the barrier (an effective quantum vibrational scale). The Fermi–Dirac distribution is f(E) = 1 1 + exp [E−EF kBT].(2) The functional forms are identical under the correspondence Ec↔EFand ℏω↔kBT, motivating an interpretation of Eq. (1) as an “occupation probability” (distribution function) [9]. In this work, we exploit this isomorphism to compute statistical weights for fission channels. 2.2 Overview of the Selective Channel Scission (SCS) model In the SCS model [7,8], nuclear fission is described as a set of channels characterized by fragment charge (Z1, Z2)and neutron number (N1, N2), subject to the conservation relations [7,8] Z1+Z2=Ztotal,(3) N1+N2+νp=Ncompound.(4) 2 Here νpis the prompt-neutron multiplicity. In this study we use representative values derived from ENDF/B-VIII.0 (MT=456) for each incident-energy condition. The weight of each channel is defined as a Hill–Wheeler-type logistic function [7,8,9] p(Z1, N1) = 1 1 + exp [k(C dw(Z1)−Q−Ex(Z1))],(5) where C= 1.44 Z1Z2(MeV·fm) is the Coulomb term, dw(Z)(fm) is the scission distance, Q(MeV) is the Q-value obtained from atomic-mass differences, and Ex(Z)(MeV) is an effective Fermi-level correction term. The scaling factor kfollows the definition in Ref. [9]. The charge yield is computed as [7,8] Y(Z1) = ∑N1p(Z1, N1) ∑Z′ 1∑N′ 1p(Z′ 1, N′ 1)×200%,(6) where 200% corresponds to the normalization for the sum of yields of the two fragments. 2.3 Parameter estimation (two-step fitting) Model parameters are estimated in two steps, following Ref. [9]: 1. Step 1: Set Ex(Z) = 0 and fit dw(Z)independently for each Z(least squares on logarithmic residuals). 2. Step 2: Approximate the Step-1 dw(Z)values by a quadratic polynomial to obtain a smoothed dw(Z), then fix this dw(Z)and estimate Ex(Z). This decomposition (i) stabilizes the estimation of the geometric scale dw(Z)and (ii) then lets Ex(Z)absorb channel-dependent corrections such as shell effects. In the low-yield tails of Y(Z), the estimation of Ex(Z)can become numerically unstable and may produce nonphysical excursions near endpoints. Accordingly, in our Python implementation the Step-2 residual evaluation excludes points with Yexp(Z)< Ymin (default Ymin = 10−3). In addition, we impose a nuclide-dependent lower bound Zmin (default 28, adjusted if needed) to cut both tails and improve robustness. 3 Computational Methods 3.1 Data sources and preprocessing •Charge-yield data: Independent fission-product yields in ENDF/B-VIII.0 [11] (processed into Y(Z)tables, e.g., using JANIS). •Atomic masses: AME2020 [12]. •Prompt neutron multiplicity: Representative values based on ENDF/B-VIII.0 [11] (MT=456). Yield data are normalized to match Eq. (6) and compared with model calculations. In general, fission-product yield evaluations in major libraries are established through a combination of numerous yield measurements, systematic evaluations, and model-based completion. 3 3.2 Python implementation and optimization The implementation assumes Python 3 (e.g., 3.10–3.12) and uses primarily numpy and scipy.optimize.least_squares for parameter estimation. We adopt the trust-region reflective algorithm with bounds •dw(Z)∈[5,40] fm, •Ex(Z)∈[−20,20] MeV. For numerical stability, a sigmoid implementation such as scipy.special.expit is useful to avoid overflow in logistic computations. 3.3 Analysis cases Table 1summarizes the analysis cases based on ENDF/B-VIII.0 (nuclide, energy condition, νp,Sn, etc.). Here Enis the incident neutron kinetic energy (MeV), νpis the representative prompt-neutron multiplicity used in the neutron-number conservation, and Sn is the neutron separation energy (MeV) of the compound nucleus formed after neutron capture. As a reference scale for initial excitation energy, we use E∗≈En+Sn. Table 1: List of analysis cases (ENDF/B-VIII.0). Log RMSE is shown rounded to 6 decimal places. Nuclide Incident energy regime 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 4 4 Results 4.1 Target nuclides Table 2lists the nine target nuclides and the incident neutron-energy conditions used in this work. Thermal-induced fission is not applicable (below threshold) for some nuclides in the present dataset; such thermal cases are omitted where they are not defined. Table 2: Target nuclides and representative calculation parameters 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 charge yields Y(Z)with the fitted SCS–HW model for incident neutron energies (thermal, 500 keV, and 14 MeV). Figure 1shows representative examples for thermal conditions (evaluated data: points; fitted model: lines). For 500 keV and 14 MeV, similarly good agreement is obtained; to avoid redundancy, full sets of figures and numerical outputs are provided in the reproducibility dataset (Zenodo). 5 233U235U237Np 239Pu 240Pu 242Pu 241Am Figure 1: Charge-yield comparison for thermal incident neutrons (points: ENDF/BVIII.0; lines: SCS–HW fit). Only nuclides with thermal fission data in the provided Y(Z)tables are shown. 4.3 Estimation of fission distance dw(Z) Figure 2shows the nuclide dependence of the fitted fission (effective scission) distance dw(Z). Points indicate the Step-1 estimates Reff (Z)obtained independently for each Z. In Step 2, we use a quadratic smoothing R(poly) eff (Z)as dw(Z)(the smoothing curve itself is not shown in the figure). For cases where a given energy condition is not defined in the provided Y(Z)tables, the corresponding series is absent. 6 232Th 233U235U 238U237Np 239Pu 240Pu 242Pu 241Am Figure 2: Estimated fission-distance points Reff (Z)from Step 1 for each nuclide, overlaid by incident-energy condition (where available). In Eq. (5), the Coulomb term C= 1.44 Z1Z2(MeV·fm) appears as C/dw(Z1), representing an effective Coulomb energy at scission. If the effective Coulomb energy at scission does not vary strongly among channels, one expects the scaling [9] dw(Z1)∝Z1Z2=Z1(Ztotal −Z1).(7) The fitted Reff(Z)in Fig. 2exhibits a single-hump behavior consistent with Z1(Ztotal − Z1), supporting the interpretation of the extracted parameter as a physically meaningful effective scission distance. For reproducibility, we express the Step-2 smoothing function as an equivalent quadratic form R(poly) eff (Z) = A+B(Z0−Z)Z[fm],(8) and summarize the parameters (A, B, Z0)in Table 3. Here Z0controls the symmetry of the parabola (maximum at Z=Z0/2); if the ideal scaling Z(Ztotal −Z)holds, one expects Z0≈Ztotal. In practice, the fitted Z0values are very close to Ztotal across all cases, indicating that the Step-1 estimates preserve the charge-conservation symmetry. 7 Table 3: Parameters of the effective fission distance R(poly) eff (Z)in Eq. (8). The Zrange indicates the charge-number range used in Step 1. Nuclide Incident energy A B Z0Zrange 232Th 500 keV -0.520783 0.000823809 89.9977 28–62 232Th 14 MeV -0.558394 0.000857709 89.9885 28–62 233Uthermal -0.550679 0.000837288 91.9948 26–66 233U500 keV -0.549342 0.000833038 92.0025 26–66 233U14 MeV -0.592277 0.000869987 91.9985 26–66 235Uthermal -0.544124 0.000820399 91.9997 28–64 235U500 keV -0.543429 0.000819832 91.9917 28–64 235U14 MeV -0.572065 0.000834452 92.0042 28–64 238U500 keV -0.525128 0.000795388 92.0020 28–64 238U14 MeV -0.571637 0.000834700 92.0059 28–64 237Np thermal -0.529254 0.000782988 92.9955 28–65 237Np 500 keV -0.536350 0.000793864 92.9990 28–65 237Np 14 MeV -0.564588 0.000815628 92.9976 28–65 239Pu thermal -0.535999 0.000775433 94.0010 28–66 239Pu 500 keV -0.547266 0.000793174 94.0041 28–66 239Pu 14 MeV -0.581857 0.000811046 94.0172 28–66 240Pu thermal -0.514895 0.000748634 93.9999 28–66 240Pu 500 keV -0.529123 0.000769088 94.0031 28–66 240Pu 14 MeV -0.571673 0.000805849 94.0053 28–66 242Pu thermal -0.524308 0.000755374 94.0032 30–64 242Pu 500 keV -0.531215 0.000765148 93.9997 30–64 242Pu 14 MeV -0.561641 0.000780989 94.0117 30–64 241Am thermal -0.513904 0.000730516 95.0032 28–67 241Am 500 keV -0.532349 0.000757247 95.0011 28–67 241Am 14 MeV -0.572137 0.000791832 94.9990 28–67 4.4 Estimation of Fermi energy Ex(Z) Figure 3shows the fitted Fermi-energy profiles Ex(Z). Table 4summarizes the maximum values Emax xextracted for thermal and 500 keV conditions and compares their low-energy mean ⟨Emax x⟩LE with literature ranges for the fission barrier energy Bf. The focus here is the low-excitation regime (thermal and sub-MeV), where penetrability near the barrier is expected to be dominant. The 14 MeV condition corresponds to a high-excitation regime; because of shell damping and possible mixing of multi-chance fission (with pre-fission neutron emission), it may be inappropriate to interpret Emax xin the same way as a barrier-related indicator. Therefore, 14 MeV is excluded from the low-energy mean ⟨Emax x⟩LE, which is computed using only thermal and 500 keV values (or only 500 keV if thermal is absent). 8 232Th 233U235U 238U237Np 239Pu 240Pu 242Pu 241Am Figure 3: Estimated Fermi-energy profiles Ex(Z)for each nuclide, overlaid by incidentenergy condition (where available). 9