Calculated covariance matrices for fission product yields using BeoH

,


Introduction
A number of applications rely on fission product yields (FPYs), such as reactor neutronics, spent fuel, dosimetry, and radiochemistry, and many of the major evaluated libraries [1][2][3] contain these data. However, only mean values and uncertainties are currently available-covariance matrices have not been developed for any of the libraries, and therefore, a format to store them has not been developed. Additionally, there has been no major update to at least the US evaluation, ENDF/B, since the 1994 evaluation of England and Rider [4], save the addition of a point at 2 MeV for neutron-induced fission of 239 Pu [5].
Recently, there has been a revitalization in the development of FPY modeling and experimental work. Modeling efforts have improved to the point where prompt and delayed fission observables-including FPYs-can be calculated consistently in a single framework [6,7], and large numbers of experimental efforts are underway to measure the energy-dependence of independent and cumulative FPYs [8][9][10] (and many more recent measurements that are still under analysis), including measurements of photon-induced fission [11,12]. Stochastic estimations of covariances for the independent and cumulative FPYs in the ENDF/B-VIII.0 and JEFF-3.3 evaluations have been developed [13]. Additionally, there has been a great deal of interest in using machine learning methods in the determination of FPYs [14][15][16][17][18], in particular where data are scarce. A recent International Workshop on Fission Product Yields was held in Santa Fe in 2019, and the Interna- * e-mail: lovell@lanl.gov tional Atomic Energy Agency (IAEA) currently has a Coordinated Research Project (CRP) focusing on modeling FPYs.
Here, we present first results of the cumulative fission product yield covariances calculated using BeoH and a Kalman filter optimization. In Sec. 2, we give a brief overview of BeoH, focusing particularly on the model parameters that are fitted, along with details of the Kalman filter optimization. Results from the covariances of 252 Cf spontaneous fission are presented in Sec. 3, and we briefly discuss covariances for fission product yields from neutron-induced reactions. Finally, we conclude in Sec. 4.

Theory
BeoH is a deterministic Hauser-Feshbach fission fragment decay model developed at Los Alamos National Laboratory [6,7]. It consistently calculates the emission of prompt neutrons and γ rays from excited fission fragments using the Hauser-Feshbach statistical theory, then performs a time-independent calculation to determine delayed observables. More details about these calculations can be found in [6,7], but here, we specifically want to discuss the parametrization of fission fragment initial conditions as these parameters go into the FPY optimization and impact the resulting covariance matrices.
The fission fragment initial conditions are distributions in mass, charge, total kinetic energy, spin, and parity, Y(A, Z, TKE, J, π) and functions of the incident neutron energy, E inc . The mass distribution, Y(A), is a sum of three and G 1,2 (A|E inc ) are the asymmetric Gaussians, In Eqs. (1) and (2), A c is the mass of the compound nucleus and W i , µ i , and σ i are the energy-dependent weights, means, and standard deviations of the Gaussians. The means and standard deviations are linear in their energy dependence, and W 1,2 are defined as for i = 1, 2, while W 0 is defined from normalization as is computed from the Wahl systematics [19], and we allow for scaling factors, f Z and f N to adjust the odd-even staggering.
The total kinetic energy, TKE, is a linear function of the incident neutron energy, and allows for a slope change at an incident energy of E 0 as indicated by some experimental data, e.g. [20], where a, b, d, and E 0 are fit to experimental data and c is determined from the continuity at E 0 . The TKE for each mass is a Gaussian with width and mean fitted to experimental data but scaled to reproduce the average TKE for the given incident neutron energy. The TKE is shared between the two fission fragments based on conservation of energy and momentum. The total excitation energy, TXE, of the fragment split is calculated from the Q value, TXE = Q − TKE, and is shared between the two fragments based on a ratio of temperatures, R 2 T = (U l a h (U h ))/(U h a l (U l )). The level density parameter, a(U), is energy dependent and calculated by the Gilbert-Cameron level density model [21]. This ratio can be taken to be dependent on the fission fragment mass, R T (A), to better reproduce the average neutron multiplicity as a function of mass, ν p (A).
The spin distribution is proportional to the available states in the level density formula, with an adjustable scaling factor, f , on the spin cutoff parameter, σ 2 l,h , where U is the excitation energy corrected by the pairing energy. Positive and negative parities are assumed to be equally probable. The initial population distribution is constructed for each fragment, and then the Hauser-Feshbach statistical decay is performed, de-exciting the fission fragment through neutron and γ-ray emission. The resulting population of each residual nucleus is collected.
The prompt observables are calculated as the weighted sums of the results of the Hauser-Feshbach calculations for each fission fragment. At this point in the calculation, we keep track of the isomeric states, M, and the independent yields have been calculated. The cumulative yields are calculated as where b jl are the branching ratios with L j total decay modes, N is the total number of nuclei that produce the i th nucleus, and δ jl, Evaluated decay data is used to determine the branching ratios, keeping consistency within the evaluated libraries. Delayed neutron emission is calculated by summing over all nuclei where the branching ratio includes a neutron emission.
For the calculations of neutron-induced fission, further information about the compound nucleus is needed that will be weighted by the multi-chance fission probabilities. These calculations include information about any pre-fission neutrons emitted, the excitation energy causing fission, and the total kinetic energy at the equivalent incident neutron energy. Again, the fission fragment initial conditions are determined for all fission fragment pairs that can be formed from each compound nucleus that will undergo fission. Each nucleus is then de-excited in the manner described above. Prompt and delayed observables are calculated consistently across incident neutron energies, providing additional constraints on the model, beyond what is available from a single observable alone.

Kalman filter
To perform the parameter optimization and calculate the fission product yield covariances, we use a Kalman filter routine [22], a type of Bayesian optimization where the likelihood and prior are taken to be multivariate Gaussian distributions. Often, a linear approximation is made, where changes in the model outputs, δf = f (x 1 ) − f (x 0 ), are assumed to be linearly related to changes in the model The parameter and output changes are connected through the sensitivity matrix, C, which is calculated for observables from BeoH by where the model difference is calculated for each cumulative fission product yield and the average prompt and delayed neutron multiplicities, ν p and ν d at each parameter x j . The parameter changes in Eq. (7) are determined through matrix algebra, where ϕ is the vector of experimental data to be fitted and the resulting parameter correlation matrix, P, is calculated Additionally, covariance matrices on the experimental data, V, and prior parameter values, X, are needed. In this work, we take the experimental covariance matrix to be diagonal, with the reported experimental uncertainties on the diagonal, however, a template of expected experimental uncertainties has been developed for fission product yields [23,24] which should be included as a future improvement. The parameter covariance matrix is also taken to be diagonal with a 1% uncertainty. Although this value is small compared to the uncertainty that we would expect on the model parameters, as discussed in the Sec. 2.2, the optimization is run in several subsequent steps, so the parameters can vary more than 1%. Once the change in parameter values and updated parameter covariance matrix are calculated, the updated models values can be calculated through Eq. (7) and compared to the exact model solution. A covariance matrix across the observables can be calculated, When the sensitivities are calculated for multiple observables, Eq. (11) provides a way to determine correlations between pairs of observables. In principle, the sensitivity matrices, C, in Eqs. (10) and (11) can be calculated for different observables but are done so in the same way, using Eq. (8).

Fitting details
We optimize BeoH using the cumulative FPYs, ν p , and ν d . The prompt neutron multiplicities constrain the magnitude of the total kinetic energy, spin cutoff factor, and Wahl scaling factors. The parameters in the pre-neutron emission mass distribution and R T (A) are additionally constrained. We fix the widths and means of the TKE distribution as a function of mass, and take all information about the compound nucleus (e.g. multi-chance fission probabilities, pre-fission neutron emission, etc.) as given by CoH 3 [25]. The multi-chance fission probabilities can have a large impact on the shape of energy-dependent quantities when crossing the openings of successive fission changes (for example see [7]); however, for the work discussed here, we take the default calculations from CoH 3 . Data for the cumulative FPYs, ν p , and ν d are taken from EXFOR.

Results
We first discuss result for the optimization of the fission product yields for the spontaneous fission of 252 Cf. Slightly fewer parameters are needed to be included in the fit than for the neutron-induced fission calculations because there is no incident energy dependence to take into account. Examples of the parameter sensitivities for the average prompt neutron multiplicity, ν p and one representative FPY, 147 Nd, are shown in Fig. 1 panels (a) and (b) respectively. We see strong sensitivity to the total kinetic energy, TKE, and spin cutoff parameter, f , in ν p , with some sensitivity to the parameters of the mass distribution.
Each FPY is also dependent on the parameters in the mass distribution, as expected; additionally, strong sensitivity to several of the R T parameters around the mass of the given fission fragment (here, A ∼ 147) is seen. These sensitivities are included in the Kalman filter optimization routine, along with the diagonal experimental covariance matrix, the initial parameter values, and a diagonal parameter covariance matrix. We perform two optimizations, one using experimental data taken from the EXFOR database [26] and one using the current ENDF/B-VIII.0 database [1]. For both fits, an iterative procedure is followed, first including ν p , ν d , and FPYs with yields of at least 5%, then including yields down to 1%, then 0.5%. For the fit to experimental data, a further step of including yields down to 0.2% was able to be performed due to the limited coverage of the experimental FPYs. Here, we will discuss both optimizations and the slight differences between them.
In Fig. 2, we show the parameter correlation matrices for (a) the fit to experimental data from EXFOR and (b) the fit to ENDF/B-VIII.0 data. Overall, the parameter correlations are very similar between the two optimizations: there are several strong correlations among the Y(A) parameters, along with TKE and f (through ν p ). There are less strong correlations along the diagonal of the correlation matrix, between the various R T values, that appear to be slightly different between the two optimizations. The differences between the two correlation matrices-P from the ENDF/B-VIII.0 fit subtracted from P determined by the fit to experimental data-are emphasized in Fig. 3. These differences are no more than 0.3 and mainly close to the diagonal in the correlations-between the R T values. This result is expected given the increased number of nuclei that are included in the fit to the ENDF/B-VIII.0 evaluation.
We then use Eq. (11) to calculate the observable covariance matrix between the various fission product yields for 252 Cf(sf). Figure 4 shows the FPY correlation matrix from the fit to experimental data. We note first that because the parameter correlation matrices for the two optimizations are very similar and they are sandwiched between the same parameter sensitivity matrix, the FPY correlation matrix shown here is very similar to the correlation matrix for the fit to ENDF/B-VIII.0 (some slight differences will be discussed later in this section). The FPYs along the x and y axes are ordered approximately by increasing ZAID (where ZAID = 1000 × Z + A), which correspond to the indices on the axes (ZAID is not given on the w a 1 w b   axes directly). We first notice that the FPYs are largely uncorrelated or have a very weak correlation between pairs, outside of some noticeable structure that will be discuss shortly. Most of the lack of strong correlations comes from the inclusion of R T (A) in the optimization procedure; when a fixed R T value is used, the FPY correlations are much stronger, although we do not show that difference in this contribution.
We can investigate the pieces of the matrix in Fig. 4 that do show significant correlations by taking a look at correlations across a mass chain and an isotopic chain. In Fig. 5(a), we show the correlations across the A = 147 chain of cumulative fission product yields. All of these fission products are highly correlated, a trend seen in almost all of the other mass chains, although these are not shown in this contribution. Figure 5 Figure 4. Correlation matrix between the cumulative fission product yields for 252 Cf(sf) using the fit to experimental data from EXFOR. Indices are ordered by increasing ZAID = 1000 × Z + A but do not correspond one-to-one to the numbers on the axes.
itive correlation. The heavier-mass isotopes are strongly correlated, much more than the lighter-mass isotopes, but there are almost no correlations between the heavier-and lighter-mass isotopes. Between the two optimization procedures, the A correlations vary less than 10% and the Z correlation vary less than 20%. The correlations between the cumulative fission product yields are therefore likely coming mainly from the model and not from the amount of data included in the optimization.
We are also developing covariances for fission product yields resulting from neutron-induced fission reactions, namely on 239 Pu, 235 U, and 238 U. For a given discrete incident neutron energy, the correlations are similar to those shown in Fig. 4, with some structure seen between isotopes but large regions of weak to no correlation. In addition to correlations between the cumulative fission product yields, we can also construct the covariance matrices between incident energies. The fission product yields are highly correlated across incident neutron energies, at least below the opening of second-chance fission, the covariances of which are still under development. However, we expect strong correlations in these higher energy regions, exactly the opposite of what is currently assumed with the lack of correlations in the current evaluated libraries. from ENDF/B-VIII.0 to constrain the model. The correlations between the FPY of various nuclei are overall weak, but we see more significant structure when considering the correlations between isotopes of the same mass or same charge.
Additionally, energy-dependent covariances are being developed for neutron-induced fission of 239 Pu, 235 U, and 238 U. Along with providing covariances between FPYs for each incident energy considered in the evaluation, covariances can also be constructed between incident energy points, providing cross-energy correlations. Because of the current optimization procedure of fitting parameters separately for each multi-chance fission channel, work must be done to insure that the covariances are reasonable as the multi-chance fission thresholds are crossed. Likely, sensitivities to all parameters will have to be calculated for the full incident energy range (thermal to 20 MeV), even if the fits are only performed within limited energy ranges. The construction of the energy-dependent covariances is still being investigated, along with a reasonable format for inclusion in a future ENDF release.