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 17, 2025 Abstract The charge yield distribution of fission products is an important input for reactor physics, nuclear fuel cycle analysis, and nuclear data evaluation, yet accurate theoretical prediction remains challenging. In this study, we reimplemented the computational procedure in Python based on a framework that incorporates the Hill–Wheeler transmission probability into the Selective Channel Scission (SCS) model and interprets it as a “quantum statistical distribution” isomorphic to the Fermi–Dirac distribution. For validation, we calculated the charge yield Y(Z)using the evaluated nuclear data library ENDF/B-VIII.0 for nine nuclides (232Th, 233U, 235U, 238U, 237Np, 239Pu, 240Pu, 242Pu, 241Am) under incident neutron energy conditions (thermal, 500 keV, 14 MeV). The results confirmed good reproducibility against the charge yield data recorded in ENDF/B-VIII.0 for all conditions. The main contributions of this paper are: (i) presentation of the Python implementation, (ii) independent validation using ENDF/B-VIII.0, and (iii) reproducible documentation including figures, tables, and procedures. Keywords: nuclear fission; Hill–Wheeler formula; Selective Channel Scission (SCS) model; charge yield; ENDF/B-VIII.0; Python; statistical model; Fermi level 1 Introduction The yield distribution of fission products (mass distribution and charge distribution) is directly related to decay heat calculations, activation, and transmutation in reactor core analysis, nuclide composition evaluation of spent fuel, and nuclear data assimilation and evaluation. However, treating the many-body quantum dynamics of nuclear fission from first principles is difficult, and in practice, one often relies on experimental data and empirical models recorded in evaluated libraries. The classical understanding of fission energy barriers and fission mechanisms was systematized by Bohr and Wheeler [1]. Furthermore, 1 recent progress has been made in the evaluation of fission yields and uncertainty (covariance) development [2,3], and efforts toward advanced modeling and evaluation of fission yields have also been reported [4]. Attempts at yield estimation and interpolation using machine learning have also been reported [5,6]. Against this background, the SCS model [7,8] is positioned as a statistical approach that decomposes the fission process into channels and provides yields through barrier transmission weighting for each channel. Furthermore, Ref. [9] proposed a framework that reinterprets the Hill–Wheeler formula [10] not merely as a transmission probability calculation formula but as a quantum statistical distribution isomorphic to the Fermi– Dirac distribution, and its effectiveness was demonstrated using JENDL-5.0 data. The objectives of this study are to (a) reorganize the above framework as a Python implementation, (b) perform independent validation using ENDF/B-VIII.0 [11], and (c) document the computational procedures and input/output in a form that allows others to reproduce the results. 2 Theoretical Framework 2.1 Isomorphism between the Hill–Wheeler Formula and the Fermi– Dirac Distribution The Hill–Wheeler transmission probability for a parabolic barrier is given by P(E) = 1 1 + exp [2π(Ec−E) ℏω](1) as described in Ref. [10]. Here, Eis the energy of the system, Ecis the barrier height, and ℏωcharacterizes the curvature near the barrier (effective quantum oscillation scale). On the other hand, the Fermi–Dirac distribution f(E) = 1 1 + exp [E−EF kBT](2) has the same functional form, and through the parameter correspondence (Ec↔EF, ℏω↔kBT), an interpretation of the Hill–Wheeler formula as an “occupation probability (distribution function)” is derived [9]. In this paper, we utilize this isomorphism for calculating the weighting (production probability) of fission channels. 2.2 Overview of the Selective Channel Scission (SCS) Model In the SCS model [7,8], fission is represented as a set of channels characterized by fragment charges (Z1, Z2)and neutron numbers (N1, N2), and the allowed channels are enumerated under the conservation laws [7,8] Z1+Z2=Ztotal,(3) N1+N2+νp=Ncompound (4) where νpis the prompt neutron emission number, and representative values from ENDF/BVIII.0 are used for each condition. 2 The weight of each channel is defined as a Hill–Wheeler type logistic function p(Z1, N1) = 1 1 + exp [k(C dw(Z1)−Q−Ex(Z1))] (5) following Refs. [7,8,9]. Here, C= 1.44Z1Z2(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 the effective Fermi level correction. kis a scaling factor following the definition in Ref. [9]. The charge yield is calculated as Y(Z1) = ∑N1p(Z1, N1) ∑Z′ 1∑N′ 1p(Z′ 1, N′ 1)×200% (6) following Refs. [7,8]. The 200% corresponds to normalization for the total yield of two fragments. 2.3 Parameter Estimation (Two-Step Fitting) Model parameters are estimated in the following two steps (following Ref. [9]): 1. Step 1: Setting Ex(Z) = 0, fit dw(Z)independently for each Z(least squares of logarithmic residuals). 2. Step 2: Approximate dw(Z)obtained in Step 1 with a quadratic polynomial, fix this dw(Z), and then estimate Ex(Z). This two-step approach assigns roles: (i) stable estimation of dw(Z)as a geometric parameter, and (ii) absorption of shell effects and other corrections into Ex(Z). Note that in the tail regions (low-yield regions) of the charge yield Y(Z), the estimation of Ex(Z)becomes numerically unstable, and non-physical spikes at the endpoints can occur. Therefore, in the Python implementation of this paper, points where Yexp(Z)is smaller than a threshold Ymin were excluded from the residual evaluation in Step 2 (the default value Ymin = 10−3was used in this study). Furthermore, a lower limit Zmin for Z used in fitting was set for each nuclide (default 28, adjusted as needed), and by cutting the left and right tail regions, instability caused by low-yield points was suppressed. 3 Computational Methods 3.1 Data Sources and Preprocessing •Charge yield data: Independent yield data (independent fission product yields) from ENDF/B-VIII.0 [11] obtained using JANIS, etc. •Atomic masses: AME2020 [12]. •Prompt neutron numbers: Representative values based on ENDF/B-VIII.0 [11] (MT=456). Yield data are normalized according to Y(Z)(Eq. (6)) for comparison with model calculations. Note that input data in the resonance region included in evaluated libraries are based on many experimental results. For example, neutron resonance parameter measurements for 242Pu [13] are an example of fundamental data for evaluations. 3 3.2 Python Implementation and Optimization The implementation assumes Python 3.x (e.g., 3.10–3.12) and uses mainly numpy and scipy.optimize.least_squares for parameter estimation. Optimization is performed using the trust-region reflective method with boundary conditions •dw(Z)∈[5,40] fm •Ex(Z)∈[−20,20] MeV From the viewpoint of numerical stability, the use of the sigmoid function (scipy.special.expit) for logistic calculations is effective for overflow avoidance. 3.3 Analysis Cases Table 1summarizes the analysis cases based on ENDF/B-VIII.0 (nuclide, energy, νp,Sn, etc.) treated in this paper. Here, Enis the kinetic energy of the incident neutron (MeV), νpis the prompt neutron emission number (representative value used for neutron number conservation in the model), and Snis the neutron separation energy (MeV) of the compound nucleus formed when the target nucleus captures a neutron. Hereafter, we refer to the initial excitation energy as approximately E∗≈En+Sn. Table 1: List of analysis cases (ENDF/B-VIII.0). Log RMSE is displayed 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 4 4 Results 4.1 Target Nuclides Table 2shows the main analysis targets (9 nuclides) and the evaluated incident neutron energy conditions. Table 2: Overview of target nuclides and calculation parameters (representative values) 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 224 241Am 95 241 2 2 4 Note: 232Th and 238U do not undergo thermal neutron fission (below threshold), so thermal neutron conditions are excluded. 4.2 Comparison of Charge Yield Distributions For the incident neutron energy conditions (thermal, 500 keV, 14 MeV), the charge yield Y(Z)was compared with the SCS–HW model fit results. Figure 1shows representative examples of Y(Z)for each nuclide under thermal neutron conditions (points: ENDF/BVIII.0 [11], lines: fit). Note that 232Th and 238U are excluded from the thermal neutron comparison figure because thermal neutron fission is not defined for these nuclides. Good agreement was also obtained for the 500 keV and 14 MeV conditions, but comparison figures showing all nuclides in parallel would be redundant, so they are omitted from the main text and included in the reproduction dataset (Zenodo) with figures and numerical data. 5 233U235U237Np 239Pu 240Pu 242Pu 241Am Figure 1: Comparison of charge yield Y(Z)at thermal incident neutron energy (points: ENDF/B-VIII.0, lines: SCS–HW fit). 232Th and 238U are excluded from this figure because thermal neutron fission is not defined for these nuclides. 4.3 Estimation of Fission Distance dw(Z) Figure 2shows the nuclide dependence of the fission distance (effective distance) dw(Z) obtained from fitting. The points in the figure are the effective fission distances Reff(Z)estimated independently for each Zin Step 1. In the subsequent Step 2, the Z-dependence of these discrete points Reff(Z)is smoothed with a quadratic polynomial function R(poly) eff (Z), which is used as dw(Z)(the smoothing curve is not shown in this figure). The points in each panel correspond to thermal, 500 keV, and 14 MeV conditions (232Th and 238U exclude thermal). 6 232Th 233U235U 238U237Np 239Pu 240Pu 242Pu 241Am Figure 2: Estimated fission distances (Reff(Z)estimates from Step 1) shown as points. Points in each panel correspond to thermal, 500 keV, and 14 MeV conditions (232Th and 238U exclude thermal). In Eq. (5), the Coulomb term C= 1.44 Z1Z2(MeV·fm) appears in the form C/dw(Z1), which represents the effective Coulomb energy at the scission point. If we assume that the effective Coulomb energy at scission does not vary significantly across channels, the scaling dw(Z1)∝Z1Z2=Z1(Ztotal −Z1)(7) is expected (Ref. [9]). The Reff (Z)shown in Fig. 2exhibits a unimodal (hill-shaped) behavior corresponding to Z1(Ztotal −Z1), supporting that the estimated fission distance parameters are obtained as quantities consistent with Coulomb geometry. Furthermore, from the viewpoint of reproducibility, the smoothing function used in Step 2 of this study is written as an equivalent quadratic expression R(poly) eff (Z) = A+B(Z0−Z)Z[fm] (8) and the coefficients (A, B, Z0)are summarized in Table 3. Here, Z0is a parameter that determines the symmetry of the parabola (maximum position at Z=Z0/2), and if the ideal scaling Z(Ztotal −Z)holds, Z0=Ztotal is expected. The actually obtained Z0values 7 were very close to Ztotal in all 25 cases, with |Z0−Ztotal|<0.02 at most (e.g., for 232Th, Z0= 89.9977 versus Ztotal = 90). This agreement indicates that the estimation results from Step 1 preserve the symmetry arising from charge conservation, and confirms that the Reff (and hence dw(Z)) extracted by this procedure functions as a physically reasonable effective fission distance. Table 3: Parameters of the effective fission distance R(poly) eff (Z)according to Eq. (8). The Zrange indicates the charge number range used for estimation 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 Fermi energy Ex(Z)obtained from fitting. Table 4summarizes the maximum values Emax xof Ex(Z)extracted under thermal and 500 keV conditions, the low-energy average ⟨Emax x⟩LE defined from these, and the literature value range of the fission barrier energy Bf. The main focus of this study is to compare Emax xin the low-excitation regime (thermal and sub-MeV region), where barrier-proximate transmission dominates, with the barrier scale. The 14 MeV condition corresponds to a high-excitation regime, where weakening of shell effects and (in some cases) multi-chance fission involving pre-fission neutron emission may occur, making it difficult to treat Emax xas a barrier indicator with the same interpretation as in the low-excitation regime. Therefore, the 14 MeV data were excluded from the average values in Table 4, and ⟨Emax x⟩LE was calculated from thermal and 500 keV 8 only. 232Th 233U235U 238U237Np 239Pu 240Pu 242Pu 241Am Figure 3: Nuclide dependence of the estimated Fermi energy Ex(Z). The three curves in each panel correspond to thermal, 500 keV, and 14 MeV conditions (232Th and 238U exclude thermal). 9