IMPROVED THERMAL SCATTERING DATA PREPARATION

Convenient access to accurate nuclear data, particularly data describing low-energy neutrons, is crucial for trustworthy simulations of thermal nuclear systems. Obtaining the scattering kernel for thermal neutrons (i.e., neutrons with energy ∼1 eV or less) can be a difficult problem, since the neutron energy is not sufficient to break molecular bonds, and thus the neutrons must often interact with a much larger structure. The “scattering law” S(α, β), which is a function of unitless momentum α and energy β transfer, is used to relate the material’s phonon frequency distribution to the scattering kernel. LEAPR (a module of NJOY) and GASKET are two nuclear data processing codes that can be used to prepare the scattering law and use different approaches to approximate the same equations. LEAPR uses the “phonon expansion method” which involves iterative convolution. Iteratively solving convolution integrals is an expensive calculation to perform (to ease this calculation, LEAPR uses trapezoidal integration for the convolution). GASKET uses a more direct approach that, while avoiding the iterative convolutions, can become numerically unstable for some α, β combinations. When both methods are properly converged, they tend to agree quite well. The agreement and departure from agreement is presented here.


INTRODUCTION
The accuracy of a nuclear system simulation is highly dependent on the nuclear data available. In particular, the quality of thermal neutron scattering data greatly impacts the reliability of thermal reactor analysis and safety margin calculations, due to the important role that thermal neutrons play in many reactor designs. Thermal neutrons have energy on the order of ∼1 eV, which is not sufficient to break molecular bonds. Thus, a neutron that hits an atom effectively scatters from a much larger target. Additionally, thermal neutron wavelengths are on the order of interatomic spacing, which means neutrons can exist atop multiple nuclei at once, further complicating these interactions [1].
Thermal scattering behavior is represented by use of the scattering law * S(α, β), which is a function of unitless momentum exchange α and unitless energy exchange β. In the incoherent approximation, it can be used to calculate the double-differential scattering cross section, where σ b is the bound scattering cross section, µ is the scattering cosine, k b T is the temperature in eV, and the E, E are the initial/final neutron energies. The incoherent approximation is a major simplification that misrepresents scattering behavior for many materials. Note that small errors in S(α, β) due to phonon distributions resolution, processing methods, etc., will likely be significantly smaller than those induced by the incoherent approximation. The scattering law can be defined in terms of the phonon distribution ρ(β) (i.e., phonon density of states, vibrational frequency distribution): Phonon distributions are typically functions of energy change, but are presented as a function of unitless β (which is unitless energy exchange). This conversion simply requires scaling the domain by k b T . Historically, empirical models have been used to approximate and simplify the phonon spectra. These simplified spectra often look like the grey distribution in Fig. 1, where the continuous distribution is cutoff and higher energy peaks are represented as Dirac δ functions. Additionally, low-energy behavior is typically smoothed out and instead represented using a diffusive/translational term. Recently, however, advanced materials modeling tools are providing fully resolved phonon distributions (such as the red distribution in Fig. 1) that do not separate distributions into translational/continuous/Dirac δ components. To take advantage of these high quality phonon distributions, thermal scattering data preparation should accept non-uniform phonon distribution grids as to allow for adequate peak resolution without becoming prohibitively costly.

METHODS OF NUMERICALLY SOLVING FOR SCATTERING LAW
If a phonon distribution ρ(β) is provided to a processing code (such as the LEAPR module of NJOY [3] or GASKET [4]), the scattering law S(α, β) can be calculated for requested α, β values. There are two main approaches to solving Eqs. 2-3: the phonon expansion (LEAPR) approach, and the numerical integration (GASKET) approach.

Phonon Expansion Approach (LEAPR)
In the phonon expansion approach, the exponential of αγ(t) in Eq. 2 is expanded as a Taylor series, Figure 1: A full phonon distribution for water [2] is shown above in red, and a corresponding simplified distribution is shown in grey.
where the quantity in the square brackets of Eq. 5 is defined as e −β/2 W n (β), bringing that where subsequent W n (β) terms are obtained by a recursive convolution [3], In Eq. 6, the n th term in the sum corresponds to the creation/destruction of n vibrational modes (phonons), which is why this method is often called the "phonon expansion method" [5]. LEAPR approximates this sum by truncating to N terms, where N is set by default to 100 [3].

Discussion
The phonon expansion method of LEAPR was used to prepare the thermal scattering data files for the ENDF/B-VIII.0 release, indicating that it is a highly trusted approach for processing thermal scattering files [6]. The phonon expansion has been employed for decades and provides adequate results, but some of its shortcomings are certainly worth consideration, namely the cumbersome iterative convolution integral calculation.
In performing the phonon expansion, LEAPR iteratively performs trapezoidal integration, which causes error accumulation and can become especially problematic if the phonon spectrum is not adequately resolved (i.e., if the β spacing in Eq. 7 is not small enough).
Additionally, the phonon expansion method heavily favors input phonon frequency distributions that are provided on uniformly spaced grids, which limits flexibility and phonon peak resolution.
If a non-uniform phonon grid is used, proper selection of the β grid in Eq. 7 is difficult, as the domain of pertinent β values expands with increasing n. Additionally, a β grid that adequately resolves W n might not necessarily resolve W n+1 . While uniformity of the phonon grid is not a hard restriction, catering to non-uniformity is cumbersome and not supported in modern phonon expansion codes.

Numerical Integral Approach (GASKET)
When solving for the scattering law, the complex exponential in Eq. 3 can be expanded as trigonometric functions and separated into a real and an imaginary piece, where This simplification brings Eq. 2 to where again the complex exponential can be expanded as trigonometric functions, leaving Note that the corresponding i sin(βt + αF (t)) term does not appear in Eq. 10 since it is an odd function and cancels out.

Discussion
While the phonon expansion method favors phonon distributions with uniformly spaced grids, there exists no such preference with the GASKET approach (due to the lack of iterative convolution). GASKET does require, however, selection of a t grid upon which F (t) and H(t) will be computed. As will be seen in Sec. 4, inadequate choice of the t-grid may result in numerical instability. But once a time grid is selected, H(t) and F (t) can be precomputed, as they only depend on the phonon distribution and material temperature. Note that as t increases, both H(t) and F (t) will oscillate rapidly and the β integral will approach zero, meaning that F (t) and H(t) can be approximated by integrating to some t max .

AGREEMENT BETWEEN LEAPR AND GASKET
Since both LEAPR and GASKET approximate the same function (Eq. 2), we expect that they return similar results. To illustrate this agreement, the scattering law for H in H 2 O at 296 K is calculated using both methods and the results are presented in Fig. 2. The phonon distribution used for this calculation was based off of the CAB water model from [2]. The scattering laws shown in Fig. 2 illustrate significant agreement between LEAPR (solid line) and GASKET (dashed line). Note that the −β side of the scattering law is plotted -this is because S(α, −β) has more reasonably sized values than S(α, +β) and is thus easier to plot in detail. Similarly, the scattering law of Be in BeO at 296 K is calculated using both methods and shown in Fig. 3. The phonon distribution is taken from the ENDF-B/VIII release [6]. For most values, the methods produce very similar results. However, for low α values and large β values, the GASKET method (dashed lines) experiences numerical instability. To maintain clarity of the figure, this instability is illustrated by setting the scattering law value to an arbitrarily small number.

TIME CONVERGENCE OF THE GASKET METHOD
As seen in Fig. 3, the GASKET method is susceptible to numerical instability for some α, β combinations. In this section, we discuss potential solutions to avoiding this instability.

Behavior of F (t) and H(t)
The GASKET approximation requires selection of a t-grid. For simplicity, consider a uniformly spaced mesh that extends from t = 0 to t = t max with spacing ∆t. F (t) and H(t) are highly dependent on the phonon distribution, so here we consider two materials with dissimilar phonon spectra: H in H 2 O and Be in BeO (the phonon distributions of which are provided in Fig. 4). The F (t) and H(t) functions corresponding to each material (at 296 K) are shown in Fig. 5. The phonon distributions differ greatly, which causes F (t) and H(t) to differ as well. Clearly, proper selection of t max and ∆t is crucial and varies for a given material.

Treating Numerical Instability of GASKET Method
In the original GASKET documentation, there exists vague empirical guidelines which recommend what values of t max and ∆t might be appropriate, depending on where "significant structure" appears in the phonon distribution [4]. Instead of requiring special consideration for each distribution, we consider now alternate generalizable approaches to avoid numerical instability.

Adaptively Refining Time Grid
An adaptively refined time grid is a reasonable attractive solution, but comes with its own complications. If the time grid is selected to adaptively refine the integrand of Eq. 10, then a different t grid would exist for each α, β combination, and selection of so many adaptive grids would be prohibitively costly. Alternatively, the time grid could be selected to simply resolve F (t) and H(t) individually -since they do not have α or β dependence, this would be a single grid that would be precomputed once and used for all α, β values. This was considered but did not result in significant improvements, since the time oscillations of F (t) and H(t) will not necessarily match the oscillations of Eq. 10.

Careful choice of t max
Numerical instabilities occur when the time grid is unable to properly represent F (t) and H(t) vanishing. This can be resolved by carefully choosing a t max where both F (t) and H(t) are very close to zero. To illustrate the importance of t max , we consider the scattering laws produced using two similar cutoff values (t max = 96.4 and t max = 105.0). The first corresponds to a region where both F (t) and H(t) are very near zero, whereas the latter corresponds to a region where F (t) and H(t) are small but still mid-oscillation. These two sets of scattering laws are shown in Fig. 6, where the number of timesteps is kept the same. Note that regions affected by numerical instability and negative S(α, −β) values are represented by setting the scattering law to 10 −8 , to keep the figure clear.
Note that when t max increases slightly, the regions affected by numerical instability greatly increase (i.e., solution becomes worse), because F (t) and H(t) are not as vanished as they are when t max = 96.4. Choosing t max to minimize |F (t max )| and |H(t max )| greatly decreases instabilities.

CONCLUSIONS
Thermal neutron scattering is an important phenomena that must be properly represented in modern nuclear system simulations. Preparing inelastic thermal scattering data typically involves calculating the scattering law S(α, β). Solving for the incoherent scattering law (Eqs. 2-3) remains a difficult problem and two approaches are well-known: the phonon expansion (LEAPR) and the numerical integration approach (GASKET). The phonon expansion, which is the more popular of the two, requires iterative convolution which heavily favors phonon distributions with uniformly spaced energy grids. GASKET is more flexible in this regard and can easily cater to phonon distributions with non-uniform grids, but is susceptible to numerical instability especially for low α values. This instability can be combated with a fine time grid, or perhaps an adaptively refined mesh. By resolving the numerical instabilities of the GASKET method, non-uniform phonon grids Figure 6: S(α, −β) was calculated using 2 timegrids: one cutoff at t max = 96.4 and the other at t max = 105. When t max = 96.5, less instability is faced because F (t) and H(t) both go to zero (as seen on the right), which contrasts with t max = 105.
could be easily used, which would facilitate use of high-quality, non-empirical phonon distributions. Lastly, it is important to note that both of these methods are numerically equivalent and operate under the incoherent approximation, which for many materials will be a comparatively outstanding source of error.