Thermal neutron scattering evaluation framework

A neutron scattering kernel data evaluation framework for computation of model-dependent predictions and their uncertainties is outlined. In this framework, model parameters are fitted to doubledifferential cross section measurements and their uncertainties. For convenience, the initial implementation of this framework uses the molecular dynamics model implemented in the GROMACS code. It is applied to light water using the TIP4P/2005f interaction model. These trajectories computed by GROMACS are then processed using nMOLDYN to compute the density of states, which is then used to calculate the scattering kernel using the Gaussian approximation. Double differential cross sections computed from the scattering kernel are then fitted to double-differential scattering data measured at the Spallation Neutron Source detector at Oak Ridge National Laboratory. The fitting procedure is designed to yield optimized model-parameters and their uncertainties in the form of a covariance matrix, from which new evaluations of thermal neutron scattering kernel will be generated. The Unified Monte Carlo method will be used to fit the simulation data to the experimental data.


SUMMARY
In this work, a thermal neutron scattering data evaluation framework is presented that combines measured scattering data and computer simulations to evaluate the dynamic structure factor (DSF), double differential cross section (DDCS), and their uncertainties.
The original parameter set of a given interaction model is randomly sampled according to interaction parameters' prior probability distribution function. For each set of perturbed parameters, a corresponding DSF and DDCS are computed, and a weight associated with this set of perturbed parameters is obtained using a Unified Monte Carlo (UMC) method from the differences between simulated and measured data. Using these weights, the best estimate of the DSF and its uncertainty is computed as a weighted average of DSF values of all perturbed parameters sets. This is the first time thermal neutron scattering kernel There's an ideal in computer science that the accuracy of a computer simulation can only be as good as input data it is given. This theorem holds true for nuclear engineering as well in the form of cross section data. With the rise in fidelity of neutron transport codes, the primary source of uncertainty is moving away from uncertainties in the solution method of transport codes and shifting towards uncertainties in the nuclear data.
This is especially concerning for thermal scattering cross section data specifically, as there are no available uncertainties in this energy region.

Thermal Scattering
Neutron cross sections can be categorized into three energy regions: thermal, epithermal, and fast. A schematic showing these energy regions is shown in Figure 1. The exact bounds between thermal, epithermal, and fast are debatable, but their general locations are shown in the plot. In the epithermal and fast energy regions, the neutron is energetic enough to render the vibrational energy of the target nucleus as well as the binding energies of a target molecule or crystalline structure as negligible. In the thermal region, however, the neutron energy is comparable to these vibrational and binding energies, meaning they must be considered when considering what the neutron cross section is at these energies.

Figure 1. Energy groupings demonstrated on 235 U
As with epithermal and fast systems, the desire for accurate nuclear data for thermal systems is crucial [1]. With the rise in interest of GEN-IV reactor systems, specifically very high temperature and molten salt reactors, there has been a need for newer, more accurate thermal scattering data. In addition to GEN-IV reactors, current light water reactors that are applying for license extensions need high fidelity cross sections and uncertainties to better quantify whether they can operate safely for another 20 years. In addition, thermal moderator data plays a key role in nuclear criticality safety analyses.
Currently, there are very limited thermal moderator data for materials that are of interest to nuclear criticality safety (e.g., Lucite, paraffin, hydrofluoric acid, etc.). The lack of uncertainties or covariance data for thermal scattering materials means that there is no way of quantifying the effects of thermal scattering uncertainties in quantities of interest in reactor systems, though there have been recent efforts to try and quantify these covariances, [2]. Additionally, there currently does not exist a method for storing the uncertainties or covariances in the ENDF file format.

Water
Water is notoriously difficult to model computationally [3]. There currently exist more than 20 unique models that describe the location of the atoms in the molecule, the distribution of charge in the molecule, and how the molecule interacts with other molecules [4]. The earliest grouping of models is based on empirical data, meaning the models are non-polarizable, use point charges to represent electrostatic forces, and a Lennard-Jones potential for dispersion and repulsion. As computers became more powerful, water models were created to better characterize the polarization effects seen in water-water interactions. Finally, with the rise of ab initio method, highly detailed models of water can be generated for use in other ab initio code systems.
Attempts to determine the thermal scattering cross section of water analytically date back to the 1960's, when the first analytical model for the double differential cross section of water was developed by Nelkin [5]. This model made several assumptions, including: approximating the normal modes of motions in terms of torsional oscillations and translational motions of a rigid water molecule plus the internal vibrations of the molecule, replacing the hindered molecular rotation with a single torsional oscillation, and small collision times. The model was the basis of the ENDF/B-III evaluation of the thermal scattering of hydrogen in light water, where no attempt was made to estimate uncertainties [6]. This evaluation was generated using the code GASKET, which improved on the Nelkin model by replacing the single torsional oscillator by a broad band of distributed modes.
No changes were made in the evaluation of the thermal scattering kernel of light water until 1994 with the release of ENDF/B-VI Release 2 [7], which kept much of the physical model from ENDF/B-III, but extended the α and β grids, which correspond to momentum and energy transfer, respectively. This new release was also evaluated using the LEAPR code, now found in NJOY [8]. This model was kept until 2006, when ENDF/B-VII was released. This evaluation was generated at Institute for Nuclear Technology and Energy Systems (IKE) [9] using NJOY, where the α and β grids were again extended and physical constants were updated to match more recent hydrogen and oxygen evaluations.

Motivation and Goals
This dissertation is motivated by the lack of thermal scattering data for various materials previously mentioned and the need to estimate their uncertainties. The recent improvements to the thermal scattering data for water made by the CAB model also shows that there is still room for development for ways to generate thermal scattering data [10]. The fact that many of the available thermal scattering data is generated using a variety of methods instead of a one-size-fits-all method is also a reason for investigation.
The goal of this dissertation is to provide a generalized framework for generating thermal scattering data and to validate this framework against available experimental data as well as experimental benchmarks to prove their improved results in real-world applications. The expression for double differential thermal scattering cross sections and Unified Monte Carlo are derived in Chapter 2. The generalized framework, as well as results comparing these new cross sections against experimental data, is discussed in Chapter 3. Validation of these new cross sections using benchmark problems are presented in Chapter 4. Finally, conclusions and future work are given in Chapter 5.

CHAPTER 2. THEORY
For completeness, the full derivation of the double differential cross section from the beginnings of a simple scattering experiment will be outlined here. This derivation is detailed in several references [11] [12], and will be summarized below. Although this derivation is based on quantum mechanics, some approximations are made in order to use classical molecular dynamic atomic trajectories.

Static Target
A thermal neutron scattering experiment can be described in the following way.
Suppose a neutron of energy , spin , and momentum ℏ is traveling towards a target as depicted below in Figure 2. The target is assumed to be static in this example, but it is not a required limitation.

Figure 2: Diagram of Neutron Scattering Experiment
If the target is assumed to remain static during the experiment, then the number of neutrons (or the count ) that hit the detector can be defined as where is the detector resolution, Φ is the neutrons crossing per unit area per unit time (also called the incident flux of neutrons), is the number of atoms in the target, and / Ω is the differential scattering cross section with respect to angle. This assume a single scattering event occurs in the target. The energy dependence of the cross section will be taken into consideration later.
The next step is to write out what the incident and scattering wavefunctions are. Since the energies of interest are on the order of several meV, the corresponding neutron wavelength is on the order of 10 -10 m, as shown below which is significantly larger than the 10 -15 -10 -14 m range associated with the nuclear forces that cause scattering. Because of this, the incident wave function can be assumed to be composed entirely of s waves, meaning that the scattering is spherically symmetric.
Assuming the neutron is scattering along the z-axis, the incident wavefunction can be written as where is the magnitude of the scattering direction . Since scattering is assumed to be spherically symmetric, the scattered wavefunction is (2.6) The differential cross section from Eq. 2.1 can be defined as the scattered flux over the incident flux: The variable is known as the scattering length, and is complex. A large imaginary part indicates that the scattering interaction results in the formation of a compound nuclear resonance in the thermal or epithermal energy range, which results in large radiative capture relative to scattering. Since only few nuclides have large imaginary components near the energies of interest for thermal scattering, imaginary scattering lengths can be ignored for this analysis.

Non-Static Target
In practical situations, however, the target does not remain static, but rather changes when interacting with a neutron. Taking this into account, the differential cross section can be rewritten as (2.14) The scattering experiment does not measure the target state or the neutron spin.
Because of this, the next step is to sum over all final target states and neutron spin states , then average over all initial over all initial target states and neutron spin states using probabilities of targets state and neutron spin states . This is done below ( ) Expanding the bracketed part of Eq. 2.15 is done using the definition of bra-ket notation (which can be found in [13]): Since there are many nuclei in the target, the potential term for the scattering system is rewritten as This can also be thought of as the momentum transfer, since ℏ is the momentum.
Additionally, the following terms are defined This is all done so that Eq. 2.18 can be rewritten as This relation cannot be further simplified until an expression for the potential is found.
The Fermi Pseudopotential that reproduces the measured bound scattering length described in the next section is what will be used as this potential.

Fermi Pseudopotential
To further evaluate Eq. 2.22, suppose there is only 1 fixed nucleus as the target. This means that = 1, and (after assuming the target nucleus is fixed at the origin so that 1 = 0 and = ) that Eq. 2.18 can be rewritten as since the target wavefunction is normalized per unit volume, where the total volume equals 3 . Then, since = , the above can be substituted into 2.12 to give (2.24) From here, an approximation of ( ) is required. Since the potential range of the nuclear forces that cause neutron scattering are on the order of 10 -15 m relative to interatomic distances on the order of 10 -10 m, the potential can be approximated by a Dirac delta function, where is a real constant. Substituting this into Eq. 2.24 gives where the last equality follows from Eq. 2.7. Therefore, Using this, Eq. 2.25 becomes This is known as the Fermi Pseudopotential. It is worth noting here that this potential is not the actual potential of the nucleus. It was derived using Fermi's Golden Rule, which is similar to the Born approximation in that they are both first-order perturbation approximations. It is the potential that, when used with Fermi's Golden Rule (Eq. 2.9), gives the required result of isotropic scattering for a single fixed nucleus.
Reconsidering the target with many nuclei in the target, Eq. 2.28 can be generalized to multiple nuclei in the target, Inserting this result into Eq. 2.20 gives which means that, when combined with Eqs. 2.15 and 2.22 From here, the Dirac delta function needs to be evaluated. To do this, it will be expanded in the time domain.

Time Domain
Before continuing, it is worthwhile to take a step back and review a bit of quantum mechanics that will be important later. Using the bra-ket notation introduced in Eq. 2.16, where ̂ is the Hamiltonian of the scattering system. This Hamiltonian can be thought of as the sum of the kinetic and potential energy operators in quantum mechanics. This can be extended to show by expanding the exponential function in a Taylor series. The closure relation [13] is also introduced here With both of these relations, the derivation can continue. In order to take the energy component into consideration in Eq. 2.31, the Dirac Delta function is expanded using (2.35) Additionally, the term in the brackets can be expanded The latter is because the scattering length is assumed to be real. Inserting this relation (along with Eq. 2.35) into Eq. 2.31 gives (2.37) From here, the bra-ket notation relation from Eq. 2.33 is applied to give (2.38) Before the final equations are derived, a few more definitions and conventions are required. First, the initial probability of states in the target, , is assumed to be given by a Boltzmann distribution: Then the Heisenberg operator is defined for convenience: It's worth noting that (0) = . Finally, the thermal operator is defined as Combining these conventions into Eq. 2.38 gives The last step involves dealing with the scattering lengths, ′ & . It is impossible to know the exact scattering length of each individual nuclide in the target due to the random distribution of spins for each isotope and of isotopes in a target. Because of this, it is easiest to simply take an average over all spin states and isotopes in target (2.43)

Coherent and Incoherent Scattering
While Eq. 2.43 sufficiently describes the double differential scattering cross section, it is worthwhile to break it down into terms that are easier to grasp. The first step is therefore to define the dynamic structure factor (DSF) of nuclide and ′ as This DSF has often times been calculated as ′ ( , ), where is the frequency given by = /ℎ, where ℎ is Planck's Constant. This work will deal with the energy to avoid confusion between the two. Next, the spin states are assumed to be uncorrelated from each other, meaning that which can be rewritten as (2.47) The first term of above equation is known as the coherent cross section, and the second term is known as the incoherent cross section. The coherent cross section can be thought of as the correlation between positions of pairs of nuclei at different times. This gives rise to interference effects, and is strongly dependent on the relative arrangement of atoms in a structure. It is this reason why coherent scattering is much more important in crystalline solids, such as graphite or polyethylene. On the contrary, the incoherent cross section depends on the correlations at different times; it does not give interference effects.
To further condense everything down, the following terms are defined In the field of nuclear engineering, the so-called "scattering law", ( , ), is defined by the unit-less variables , which correlates to momentum transfer, and , which correlates to energy transfer. The scattering law, as well as the unit-less variables and are defined as

Van Hove Theory
The crux of the Van Hove Theory is that scattering can be determined by scattering functions [14]. The space-time correlation function is defined as as well as a loss of the relation of odd moments

Gaussian Approximation
An alternative to the Van Hove theory is to calculate the DSF directly using the Gaussian approximation. This approximation only applies to the incoherent component of the double differential scattering cross section, and is sometimes referred to as the "Incoherent approximation". In this scenario, the space-time correlation function is defined as where Γ( ) is the width function, which is interpreted as the mean square departure of the particle from the origin after time t. The associated intermediate structure factor is therefore The width function can be thought of as the mean square departure of a particle at time t from its location at t=0. There are multiple ways to calculate this width function [19]. In this framework, the width function will be calculated using the frequency distribution (also known as the density of states) as shown below where ( ) is the frequency distribution. This particular width function is valid for an arbitrary target, regardless of whether it's a solid or liquid.
This process of obtaining the incoherent scattering law is what is used by NJOY [8].
It can be shown that the subsequent DSF satisfies the universal detailed balance relation mentioned in Eq. 2.62, as well as the first odd moment in Eq. 2.63.

Unified Monte Carlo
Nuclear data at energy ranges above the thermal neutron energy group has traditionally been evaluated using the generalized least squares (GLS) method. This method is based on the principle of maximum entropy and is applicable to many situations, but does have the drawback of requiring sensitivities that would require modifications to the molecular dynamics code [20]. Because of this limitation, a newer, less restrictive model was deemed necessary for evaluating thermal scattering data. The Unified Monte Carlo (UMC) method was first described by Smith [21] and follows from Bayes Theorem and the Principle of Maximum Entropy [22]. In the following description, an experiment containing experimental data points is represented by , it's associated covariance matrix is , simulation results (calculated using a nuclear model, or some computer code) containing data points is represented by , and it's associated covariance matrix is . Bayes Theorem gives the posterior probability density function (PDF) ( ) in the following form where is a normalization constant, ( , | ) is a likelihood PDF (dependent on experimental data), and ( | ) is the prior PDF (dependent on the simulation data), and is a collection of random variables. The normalization constant is chosen so that the posterior PDF integrates to unity when integrated over the entire domain space. Using this notation, the mean value of each random variable in the collection of random variables and the elements of its covariance matrix are defined as There are currently two different UMC methodologies that are used: UMC-G and UMC-B.

UMC-G
The Principle of Maximum Entropy states that, if a collection of random variables is summarized by only their mean values and covariance matrix, the optimal choice for the PDF is a multivariate Gaussian function. This leads to the following likelihood function shown, however, that Metropolis algorithm will converge on an answer several orders of magnitude quicker than a brute force method would [21].
The UMC-G method has the benefit of creating an analytic approximation for the posterior that can be readily sampled. One of its main shortcomings is the need for the simulation covariance matrix. It also has issues with higher-order distribution moments, leading to biases in cases where non-linear effects and distribution skewness and kurtosis are present.

UMC-B
An alternate to the UMC-G method is the UMC-B method [23]. This method came about from the realization that certain analyses excessively rely on nuclear modelling and A benefit to the UMC-B method, aside from not needing to calculate the simulation covariance matrix, is that all information in the prior function is preserved, including the non-linear terms neglected by the UMC-G case. A drawback, however, is that the sampling range for the nuclear model parameters must be sufficiently large to ensure that there are no biases. This means that there will be model parameters sampled that can lead to un-physical results and must be rejected, which wastes computational resources.

Water Models
As mentioned in the introduction, water is very difficult to model. Here the empirically-based models of water models will be discussed. These models, rather than the more accurate polarization or ab initio methods, are used because they will be used in a classical molecular dynamics (MD) code system that cannot handle ab initio models. In  , where is the permittivity of free space, is the dielectric constant, and are the charges of particles and , and is the distance between particles and . Because there is no limit to the range of this potential, most MD codes will cutoff interactions longer than a given length. Long range electrostatics are handled using particle-mesh Ewald method [24].
The Lennard-Jones potential describes how neutrally-charged particles and molecules interact with each other. This term has the function form of ( ) = 4 (( ) 12 − ( ) 6 ), where is the depth of the potential well, is the finite distance between particles & at which the potential changes from repulsive to attractive, as seen in Figure 4. The potential has a repulsive effect (corresponding to a positive value) for small distances < , and then an attractive effect (corresponding to a negative value) for distances > while asymptotically approaching zero. Like the Coulomb potential, there is no maximum range for the Lennard-Jones potential, so computer codes use approximations for when to stop calculating its contributions to the potential of a particle.

Figure 4. Lennard-Jones Potential
Bonded terms, on the other hand, deal with the interactions within an atom. In water, there are 2 dominant bonded interactions that are used: bond stretching and bond angle.
The most common bond stretching form that is commonly used in water is the harmonic potential, which has the form where is the bond strength between particles and , and , is the equilibrium length between particles and . There are a few potentials, however, where anharmonic bond stretching is required. For these cases, the Morse potential is used: where is the depth of the potential well and defines the steepness of the well. A plot comparing the two bond potentials is shown below in Figure 5. It can be seen that, for increasing radius, the Morse potential does exhibit slightly anharmonic behavior.
where is the bond strength, is the angle between the three atoms, and 0 is the equilibrium angle between the three atoms.

CHAPTER 3. FRAMEWORK OVERVIEW & RESULTS
Here, the specific details of the framework are presented and discussed. The light water experimental dataset used to validate the framework is first described. The process for generating thermal scattering kernels is detailed, and the double differential cross sections resulting from the framework are then validated against other experimental data not used in the original UMC fitting procedure.

Experimental Data
The experimental data was gathered from the Fine-Resolution Fermi Chopper

Spectrometer (SEQUOIA) detector at the ORNL Spallation Neutron Source (SNS) in
2005 by a research group from Rensselaer Polytechnic Institute (RPI). The data was collected at incident energies of 55, 160, 250, 600, 1000, 3000, and 5000 meV between scattering angles of 3°-58° with 1° increments. The energy resolution for the double differential cross sections are 0.5 meV for the 55 meV case, 1 meV for the 160 meV case, and 2 meV for the remaining 5 cases. Each experiment was carried out at a temperature of 300 K.
The SEQUOIA detector is a time-of-flight spectrometer, which works by limiting the incident neutron energy to one specific energy using metal cylinders called choppers.
These choppers have a wedge holed out of them that, when rotating at a specific speed, allows only energy of neutron through. These neutrons then hit the sample, which is a 0.1mm thick sample of light water in an aluminum can. The sample was chosen to be thin to reduce the effects of multiple scattering. The scattered neutrons then travel to one of the detectors, which tally the time when the neutrons arrived, meaning that the scattered energy can be calculated. Two separate experimental runs were performed and averaged to improve the accuracy of the results. A third run where the water is removed is then done to get the effects of the aluminum can, so that it can be subtracted off the previous two runs.

Framework
To generate the mean DSF, the UMC method will be used to compare against the experimental DDCS described above. The decision to use DDCS was based on the fact that it is an experimentally measurable quantity, while the DSF cannot be experimentally measured. There is no reason why total cross section could not be used instead, but the DDCS measurements should give a better understanding into the totality of the DSF, while the total cross section would give an understanding to its integral properties, which may hide some underlying that would be found by comparing to the DDCS.
The first step of the framework is to obtain the trajectory data for water. For this, the code GROMACS [25] was used due to its highly customizable input parameters. In addition, the ability to create new parameter files means that any feasible model for water can be used. For this case, the TIP4P/2005f potential [26] will be used. This was chosen because the potential was originally fitted to yield a more accurate frequency distribution  From this equation, a total of eight parameters will be perturbed based on a Gaussian distribution: , , , , , , and . Not shown in the above equation is the parameter that dictates the distance the dummy particle is from the oxygen atom, , which will also be modified. The published values of these variables are shown below in Table 1. The molecular dynamics simulation is broken up into 4 steps, as outlined in [27]. Once the MD trajectories are saved, the next step is to use these trajectories to calculate the DSF. For this application of the framework, the Gaussian approximation (described in Section 2.1.7) is used. This approximation is valid for light water because the incoherent scattering cross section for light water (160.54 b) is much greater than the coherent cross section (7.7486 b). To calculate the DSF using the Gaussian approximation, the density of state is required. This is calculated using the velocity autocorrelation as shown in Eq To better compare these simulated results against the experimental data, a simplified model of the SNS detector was modeled in MCNP, where a monoenergetic beam of neutrons was fired at a cube of water, and the scattering results were tallied at rings meant to represent how the SEQUOIA detector tallied the scattering events. Since MCNP is being used, the DSF is converted to ACE format using NJOY [8]. The SEQUOIA detector resolution was applied afterwards. Finally, the UMC procedure is done. In this framework, the UMC-B method will be used, as it has the benefit of not requiring the calculation of the covariance matrix of the simulation data.
Before the framework was fully implemented, a small set of ensembles were run where the TIP4P/2005f parameters were not perturbed. Instead, only the initial position of the atoms in the simulation was changed. This was done to see how the weighting function might change between ensembles were the only difference should be statistical noise and the random movements of the molecules. This run yielded a surprising conclusion; the overall variation of the functions is not very large, but the magnitude of the functions is quite small, which points to a potential issue. Due to the energy spacing in the data (0.5 meV for 55 meV data, 1 meV for 160 meV data, and 2 meV for the other incident energies), the sheer volume of data meant that the UMC procedure could not be carried out directly. Because of this, a replacement for the weighting function in Eq. 2.74 was required. Previously, this value has been too large, there have been a couple of recommended options [22]. The options include: verifying that the simulation model can sufficiently describe the data, labeling certain experimental data points as 'questionable' and enhancing their specific uncertainties, or spreading the discrepancy across all data points equally. This work chose the latter of the options, as there were many data points that were found to be discrepant, and changing the model would be infeasible.
The framework was run twice; once varying the parameters by 5% each, and then by varying the parameters by the amounts given from applying UMC to those results. The where the parameters are varied by 5%, are listed below in Table 2. These values were then used as a prior for the second iteration of the framework, where the parameters were perturbed based on the standard deviations in Table 2. With this implementation, a further restriction was applied; of the 3464 simulations that ran to completion, only 60 reasonable results for property quantities of interest of light water (density, relative static dielectric constant, isothermal compressibility, dipole moment, and diffusion coefficient). Of these 5, the most important for the purpose of thermal scattering is the diffusion coefficient. This is because the diffusion coefficient is calculated directly from the frequency distribution, which is the primary input for calculating the intermediate structure factor from Eq. 2.66.

TIP4P/2005f Parameters and Properties
First, the UMC method was used to determine what the updated potential values for the TIP4P/2005f model should be. These new results and uncertainties are shown in Table 3, along with the original parameters of the TIP4P/2005f potential (from Table 1) and its percent difference from the original values. A plot of the distribution of these parameters is shown in Figure 6. The correlation matrix of the parameters is also  Almost every percent standard deviation of the values in Table 3 decrease from those   in Table 2, with the exception of the potential well in the Lennard Jones potential ε. The large uncertainty on the equilibrium scattering angle θo may indicate how insensitive it is in thermal scattering. This is backed up by the fact that, between various other potential models for water (TIP3P, SPC, etc.) the scattering angle can vary from 104.52°-109.47° [28]. The values shown in Table 3 Table 2.
The distribution of the parameter is mostly Gaussian, which is to be expected, as the parameters were sampled using a Gaussian distribution. One of the interesting outliers is how the distance from the hydrogen and oxygen molecule, bo, appears to be decently larger than the published value. There are other parameters (namely Kθ) that differ by a wider margin than bo, but the fact that only 2 simulations produced values smaller than the experimental value of the hydrogen-oxygen distance is interesting.  The correlation matrix reveals some interesting properties of the TIP4P potential.
There is a strong inverse correlation between the equilibrium angle θo and the depth of the potential well in the Lennard Jones potential ε, which is not inherently obvious by inspecting the formula. Additionally, there is almost no correlation between the potential well ε and the distance between particles σ. This is especially odd, as they are both in the same equation for the Lennard Jones potential in Eq. 2.78. When combined with the fact that the potential well ε had a jump in percent standard deviation as previously discussed, this may point towards an issue with the value assigned to this specific parameter.
In addition to the TIP4P parameters, the properties mentioned before were calculated using the UMC-B method. These results are shown in Table 5 with a plot of the distribution of these properties in Figure 7. The biggest differences come from the dipole moment and relative static dielectric constant. While these disagreements are quite large compared to the experimental values, this is expected when using the TIP4P/2005f potential, as shown in Table 6 [26]. Comparing the two tables, the UMC results give better agreement for the dipole moment, relative static dielectric constant, and diffusion coefficient by noticeable margins. The UMC analysis does do worse than the original in calculating density and isothermal compressibility, but not by excessively large margins. Regarding the distribution of properties shown in Figure 7, the dipole moment, isothermal compressibility, and temperature appear to follow a Gaussian distribution, while the relative static dielectric constant, density, and diffusion coefficient do not. The dielectric constant and density appear to follow more of a beta distribution, though the reason for why they do is unclear. In addition, there doesn't appear to be any sort of distribution that the diffusion coefficient follows, which does not make intuitive sense.
This warrants future investigation into why these properties appear to follow their respective distributions.

Double Differential Cross Section
As a first step of validating the new thermal scattering data, they are plotted against the experimental data from the SNS, which can be found in APPENDIX A. In these plots, the green band represents the simulation data +/-1 standard deviation away from the mean value. These were created by perturbing the mean DSF by the uncertainty of the DSF multiplied by either +/-1. Unfortunately, this is synonymous to assuming that the correlation matrix is a full matrix of ones, which is not true. The alternative would be to randomly vary the individual values of the DSF by their uncertainties, which implies a correlation matrix of ones along the diagonal and zeros elsewhere. Doing this, however, yields unphysical results that cause NJOY to break. Because of this, the full correlation matrix is assumed for now until a better covariance matrix can be calculated. The data were generated by running a simplified MCNP model meant to recreate the SEQUOIA detector, with the detector resolution function applied afterwards.
Overall, the simulation results seem to agree favorably with the ENDF/B-VIII.β3, which makes sense since the same molecular dynamic code and light water model was used for both. The differences come about from the modification of the parameters made in this work. There are a couple of weird quirks in the data that should be mentioned. In the 55 meV data plots, the experimental results at 35° disagree significantly with the simulation results, which contrasts the general agreement they've shown in the other scattering angles. This is most likely due to an error in the data, as the results at 34° and 36° agree much more favorably. Another instance occurs with the 160 meV plot at 35°, where there is an odd indentation at the peak of the experimental data which is also not exhibited at 34° and 36°.
Additionally, the thermal scattering data is plotted against independently gathered experimental data [29]. The data from this source did not include any uncertainties, so no attempt is made to assume what they may be here. The plots of these cross sections are shown in Figure 8 and Figure 9. Both of these sets have simulation data that are convoluted with a Gaussian resolution function to best approximate the detector resolution function. The 151 meV data assumes a Gaussian function with σ=7 meV, and the 304 meV data has σ=9.6 meV. As with the SNS data set, the simulation data appears to agree with both the ENDF/B-VII.1 and ENDF/B-VIII.β3 data. All sets of these data, however, appear to diverge from the experimental results at scattering angles greater than 90°, which corresponds to a back-scattering event. This may point to a slight deficiency in how the underlying theoretical models handle back-scattering.

Total Cross Section
Finally, the total cross section was plotted against several other experimentally gathered cross sections [30] [31] and is plotted in Figure 10. As with the DDCS in the previous section, the green band represents the simulation data +/-1 standard deviation away from the mean value. Here, however, the uncertainties were calculated using two different methods. The first method (Method 1) was calculated by perturbing the mean DSF by the uncertainty of the DSF multiplied by either +/-1 and calculating the total scattering cross section. The second method (Method 2) uncertainties were generated by applying the UMC-B method of directly on the total scattering cross section (using the UMC-B weights that were previously calculated with the DDCS), finding its associated mean and uncertainty, then perturbing the mean total cross section by the uncertainty of the total cross section multiplied by either +/-1. In both methods, the absorption of hydrogen and the total oxygen cross section were added.
Overall, the simulation results for the first method are slightly less accurate at low energies (less than 1 meV) compared against the ENDF/B-VIII.β3 results. The simulation results are also noticeably larger than both ENDF results in the 1-10 meV range. This is particularly strange, as the double differential cross sections seem to agree favorably with the ENDF results. The results from Method 2, however, are noticeably better than the Method 1 results. Specifically, the low-energy region agrees more favorably with the experimental data, and the uncertainties in Method 1 are significantly larger than the uncertainties in Method 2. To understand why there are such significant differences between the simulation results and ENDF results, a plot of the normalized DSF is shown in Figure 11, where the top plot has a linear y-axis and the bottom plot has a logarithmic y-axis. It is clear that, the further from the peak the structure factors get, the more they diverge. Specifically, the simulation results appear to be larger than the ENDF/B-VIII.β3 results. This is most likely due to the method which is used to calculate the DSF. The ENDF (both ENDF/B-

PST-033-003
The  Figure 13, and a XY view of the model is shown in Figure   14.
Material compositions for the model are shown given in Table 7.   Table 7. The experimental criticality was found to be 1.000. There are no biases for the benchmark problems, but there are several sources of experimental uncertainties.
Uncertainties in the temperature of the solution account for 16 Figure 15, and a XY view of the model is shown in Figure 16.
Material compositions for the model are shown given in Table 8.   Table 8.

HCT-006-003
The HCT-006-003 experiment is a water-moderated hexagonally pitched lattice with  Figure 17, and a one-quarter view from the XY plane of the model is shown in Figure 18. Material compositions for the model are shown given in Table 9.  Table 9.  Table 9. The experimental criticality was found to be 1.000. There are some uncertainties in the material composition and exact experimental setups that contribute to some uncertainties. Specifically, uncertainties in the pitch, fuel cross section area, and length account for 28 pcm uncertainty, while uncertainties in the fuel mass, enrichment, and clad mass contribute 44 pcm uncertainty to the criticality. For this framework, the simplified model is used. The simplifications in this model mean that, instead of the cross-shaped fuel rods, cylindrical fuel rods are used. This introduced a bias in the criticality that reduced the criticality by 0.0231. The spiraling in the fuel rods was also removed for the model, but it was determined that this introduced a negligible difference in the criticality, and as such is ignored. The uncertainties mentioned before were also slightly increased in the simplified model to account for statistical uncertainties in the estimated change of the criticality. With these factors, the benchmark-model criticality was found to be 0.9769 ± 0.0049.

Benchmark Results
To validate the new cross sections, the benchmarks were run using the ENDF/B-VII.1, and ENDF/B-VIII.β3 libraries as well as the new simulated cross sections. Since the uncertainties cannot be propagated through MCNP, the simulated data was perturbed by the simulated uncertainties and run. Specifically, the DSF was created by perturbing the mean DSF by the uncertainty of the DSF multiplied by either a normally distributed random number. This should give a rudimentary first approximation of the sensitivity of the system to thermal scattering cross sections. The simulations were all run using MCNP6.1, with each input being run such that the stochastic uncertainty would be 4 pcm.
A general observation about all three benchmarks is that the simulation results show a greater keff than the ENDF/B-VII.1 or ENDF/B-VIII.β3 results. This can be attributed to the fact that, as pointed out in the total cross section plot in Figure 10, the simulation cross section is greater than the either ENDF libraries below 10 meV. The benchmarks were originally chosen since they each exhibited a negative sensitivity to perturbations in the hydrogen cross section in the 1-100 meV range, and it was thought that The PST benchmark shows that the simulation results get closer to the benchmark results, and they are not overly sensitive to thermal scattering, based on the change of eigenvalue between the maximum and minimum values (26 pcm). The difference between the ENDF/B-VIII.β3 and simulation results is 61 pcm, which is a greater difference than the previously mentioned 26 pcm range of the 6 variations to the cross sections. The LCT benchmark results also show that the simulation results get closer to the benchmark results, but the ENDF/B-VIII.β3 and simulation results both fall within the standard deviation of the benchmark uncertainty, so it is not as meaningful as with the PST benchmark. Even with an eigenvalue difference of 78 pcm between the maximum and minimum variations, the varied simulations fall within the uncertainty of the benchmark. The HCT results in show the simulation seems to do worse than the ENDF/B-VIII.β3 library at calculating the eigenvalue, and it is the most sensitive, exhibiting a difference of 103 pcm between the maximum and minimum variations. This specific benchmark has a very large uncertainty (490 pcm), and the ENDF/B-VIII.β3 and simulation results are both just outside this uncertainty threshold.
As a further check, the benchmark problems were also run with the ACE files from the 60 accepted ensembles. These plots are shown in Figure 19 for the PST benchmark, Figure 20 for the LCT benchmark, and Figure 21 for the HCT benchmark. In each of these plots, the first data point represents the benchmark, and the next 60 are the simulation runs. Of the 60 simulations, 8 runs did not finish due to MCNP losing track of the particles. These runs are shown as a red x equal to the arithmetic average of the other 52 runs. and 109 pcm, respectively. These are much more closely packed together than the results from Table 10 to Table 12, which were 26, 78, and 103 pcm, respectively. This may indicate that the previous method of perturbing the DSF, adding the uncertainty times a Gaussian random number, is not a good way to evaluate the sensitivities.

CHAPTER 5. CONCLUSION AND FUTURE WORK
A new methodology for generating thermal neutron scattering kernels has been developed that combines information from experiments and computer simulations which evaluates not only thermal scattering kernels, but their associated uncertainties and covariance matrices. The application of the UMC method shows that, while previously only used in fast-spectrum data, it can be used for thermal scattering data. This is the first While incoherent approximation used in this work is justified for water, other methods of evaluating the DSF, such as the Van Hove theory or various ab initio models can be used due to the generality of the evaluation framework. Since the Van Hove theory is currently limited to being calculated after applying classical approximations, generating the atomic trajectories using more accurate molecular dynamics methods (such as ab initio molecular dynamics) are expected to yield more accurate results. There has been a recent attempt to combine the Gaussian approximation with the Van Hove theory [35], which could be expanded to include coherent contributions.
This work presents a Monte-Carlo estimate of sensitivities of several integral benchmark experiments to uncertainties in thermal neutron scattering kernel. This is the first step in addressing the absence of methods of propagating thermal neutron scattering covariance data through neutron transport codes and the related calculation of sensitivities of integral benchmark experiments to thermal neutron scattering data. It is hoped that this work will encourage development of other innovative methods for generating thermal scattering data and its covariance. While the framework allows for the construction of a covariance matrix, one was not calculated in this current implementation, and will be included in future publications. The presented framework can be applied to any material, including other materials of interest for thermal reactor applications (such as graphite, Silicon Carbide, or FLiBe), as well as materials important for criticality safety applications (such as lucite, Teflon, and polyethylene).