Treatments of Thermal Neutron Scattering Data and Their Effect on Neutronics Calculations

In the conventional method to generate thermal scattering cross section of moderator materials, only one of the coherent elastic scattering and incoherent elastic scattering is considered in neutronics calculations. For the inelastic scattering, fixed incident energy grid is used in the nuclear data processing codes. The multipoint linearization method is used to refine the incident energy grid for inelastic scattering. We select ZrHx (zirconium hydride) as an example to analyze the effects of the above described treatments on the reactivity of several critical benchmarks. The numerical results show that the incident energy grid has an obvious effect on the effective multiplication factor (k eff) of the analyzed reactors; simultaneously considering the coherent and incoherent elastic scattering also affects k eff by tens of pcm.


INTRODUCTION
In the neutronics analysis of nuclear reactors, accurate prediction of the thermal neutron distribution has an important effect on behaviors of the reactors. It is necessary to provide accurate thermal neutron scattering cross sections for neutronics codes to simulate the neutron thermalization. In the thermal energy region, the neutron scattering is sensitive to the atomic structure and motion in a moderator. In the Evaluated Nuclear Data Files (ENDF), the thermal scattering law (TSL) data are provided for some moderator materials to describe the thermal scattering of the bound atoms. Currently, individual TSL data files are contained in the major ENDFs, such as ENDF/B-VIII.0 (Brown et al., 2018), JEFF-3.3 (Plompen et al., 2020), and JENDL-4.0 (Shibata et al., 2011). The modern ENDFs adopt a common format, namely, ENDF-6 (Trkov et al., 2012), to store the evaluated nuclear data. In ENDF-6 format libraries, a specific file (MF) is used to store a certain data type, and among the files, File 7 (MF 7) contains the TSL data for moderator materials.
Although several efforts have been made in Monte Carlo codes to directly use the TSL data (Čerba et al., 2013;Liu et al., 2016;Hart and Maldonado, 2017), for most cases, the TSL data should be first processed to calculate total scattering cross section and double-differential cross section and then converted to a specific format required by neutronics codes. For Monte Carlo codes, the obtained cross sections are converted to tabular data representing the energy and angle distributions of the secondary neutrons and stored in the ACE (A Compact ENDF) (Conlin and Romano, 2019) format library, whereas for deterministic-based codes adopting the multigroup approximation, the tabular data are further converted into multigroup cross sections and scattering matrices. The nuclear data processing codes, such as NJOY (Macfarlane et al., 2017) and NECP-Atlas (Zu et al., 2019), are designated for the above processing.
In the ENDF-6 format, the TSL data for elastic and inelastic scattering are provided in File 7 with different reaction number (MT), MT 2 for elastic scattering and MT 4 for inelastic scattering. According to the theory to generate the TSL data, there are two components in the elastic scattering, namely, coherent elastic scattering and incoherent elastic scattering (Squires, 2012). But, only one elastic scattering mode is given in the ENDF for a certain material, and the other one is ignored. In other words, the coherent and incoherent elastic scatterings are not simultaneously provided for a material. For example, the coherent elastic scattering is given for metal beryllium, and incoherent elastic scattering is given for the hydrogen bound in zirconium hydride (ZrH x ). The nuclear data processing codes just calculate corresponding elastic scattering cross section based on the data given in the ENDF. No works have been reported to show the effect of this treatment on neutronics calculations.
As for inelastic scattering, the ENDF provides the so-called S(α, β) data. The S(α, β) data are converted to a discrete tabular data for Monte Carlo calculations. The work by Conlin et al. (2012) showed that this representation can introduce noticeable deficiencies for differential calculations and recommended adopting a continuous S(α, β) table to represent the secondary energy and angle distributions. In the work by Hartling et al. (2018), it was found that in NJOY the inelastic scattering cross section is calculated on a fixed incident energy grid, and it has obvious effect on the Monte Carlo simulations. To solve this problem, an adaptive incident energy grid was implemented in the nuclear processing code NDEX (Wormald et al., 2020).
Recently, systematic researches have been done to calculate TSL data and thermal scattering cross sections in the nuclear data processing code NECP-Atlas. An advanced TSL data calculation module, called sab_calc , was developed. Using this module, accurate TSL data has been obtained for some materials, for example Be, ZrH x . The TSL data can be directly used by the therm_calc module (Zu et al., 2019) to calculate total thermal scattering cross section and double-differential cross section. In the present work, we will investigate the effects of the aforementioned treatments for thermal neutron scattering data on the neutronics calculations. The analysis is carried out with the Monte Caro code NECP-MCX .

Coherent Elastic Scattering Cross Section
According to the ENDF-6 format, the double-differential cross section of coherent elastic scattering is calculated as follows: where E is incident neutron energy; E′ is secondary neutron energy; T is temperature; μ LAB is the scattering cosine in the laboratory reference system; S(E, T) and μ i are obtained as follows: In Eqs 2, 3, E i is the Bragg edges. The variable S i (T) is not given in ENDFs, but S(E, T) is actually provided for each E i . The total scattering cross section can be calculated as follows: In the ACE format library, the tabular data for coherent elastic scattering consists of the Bragg edges, and corresponding S(E, T) are stored in the library.

Incoherent Elastic Scattering Cross Section
The double-differential cross section of coherent elastic scattering is calculated as follows: where σ b is the characteristic-bound cross section; W(T) is the Debye-Waller coefficient for temperature T. These two variables are given in ENDF.
In the ACE format library, the tabular data for incoherent scattering contain the energy grid, total scattering cross section, and outgoing angular distribution. The total incoherent scattering cross section is obtained as follows: For the angular distribution, the equally probable discrete cosine is stored in the library, which is calculated as follows: where N is the number of cosine bins; i is the index of cosine bins; μ i is calculated as follows: with μ 0 −1.

Inelastic Scattering Cross Section
For inelastic scattering, the double-differential cross section can be calculated using the TSL data provided by ENDF as follows: where k B is Boltzmann's constant; σ b is the characteristic-bound cross section for the target nuclide; S(α, β, T) is given in ENDF for Frontiers in Energy Research | www.frontiersin.org December 2021 | Volume 9 | Article 779261 temperature T; α and β are momentum transfer and energy transfer, respectively, which are calculated as follows: In the nuclear data processing code, the S(α, β, T) data are first used to calculate the double-differential cross sections at a certain incident energy grid according to Eq. 9, and then the doubledifferential cross sections are transferred to tabular data. The tabular data contain the scattering kernel P(E → E′) representing the probability that the neutron with incident energy E exits with energy E′ after scatter and the equally probable discrete cosines for each incident energy. P(E → E′) is calculated as follows: where σ(E → E′) is differential cross section obtained by integrating Eq. 9 with respect to the outgoing cosine over [−1, 1]; σ(E) is the total scattering cross section obtained by integrating Eq. 9 with respect to the outgoing energy and cosine.
In NJOY, an incident energy grid with 118 points is fixed in the source code. As mentioned previously, some works have found that the fixed incident energy grid shows obvious effect on the Monte Carlo simulations. Figure 1 shows the secondary distribution of the inelastic scattering from ZrH x . At low incident energy region, the secondary distribution slowly varies with incident energy E, and it seems reasonable to use fixed incident energy grid. However, in the higher energy region, the shapes of the secondary energy distribution show obvious differences for different incident energies. Using the cross section calculated at fixed incident energy grid to interpolate cross sections at other energies will introduce a large error.
In this work, we adopted the multipoint linearization method to refine the incident energy grid as follows. First, an initial incident energy grid is set, which can be selected from the fixed energies in NJOY. For each incident energy, the secondary energy of scattered neutrons is divided into two parts: the downscattering part and up-scattering part, as shown in Figure 2. Several secondary energy points are respectively set in the downscattering and up-scattering parts to check whether scattering probability can be linearly interpolated between two adjacent incident energies. In each secondary energy part, the set energy points divide the part into different intervals with identical logarithmic width. Therefore, for the down-scattering part, the energy points are determined as follows:  where E min ′ is the minimum secondary energy; N is the total number of energy points added in the down-scattering part and is decided by the users; M is the index of the energy points. Similarly, in the up-scattering part, the energy points are determined as follows: where E max ′ is the maximum secondary energy. In addition to the energies determined as Eqs 13, 14, the secondary energy that is equal to the incident energy is also used to check the convergence during the linearization procedure. The linearized procedure is described as follows.
For two adjacent incident energies and their midpoint, the scattering kernel, from incident energy E to each secondary energy E′ defined as Eqs 13, 14, is calculated according to Eqs 13, 14. The scattering kernel for the midpoint is calculated again by linearly interpolating as follows: where E 1 and E 2 are adjacent incident energies; E mid is the midpoint between E 1 and E 2 ; E 1 ′ , E 2 ′ , and E′ are the secondary energy points in the panel for the incident energies E 1 , E 2 and E mid , respectively. If the tolerance between the exact P(E mid → E′) calculated using Eq. 12 and the one obtained by linear interpolation is less than a preset criterion, it is considered that the scattering kernel can be linearly interpolated. In this case, the E 1 and E 2 will be included in the final grid, and the midpoint will be removed. Otherwise, the midpoint is added to the final incident energy grid, and the interval-halving technique (Cacuci, 2010) is used to subdivide the interval between E 1 and E 2 , until P(E mid → E′) can be linearly interpolated between two adjacent incident energies.

RESULTS AND DISCUSSION
We select ZrH x as an example material to analyze the treatments for thermal scattering data on neutronics calculations. Because the TSL data for coherent elastic scattering are not provided for ZrH x in the ENDFs, in this work, we calculated these data for ZrH x by the sab_calc module, which is developed based on the phonon expansion method (Squires, 2012;Wormald and Hawari, 2017). In the calculation of TSL data, the phonon density of states (PDOS) is the fundamental data. In this work, the PDOS of ZrH x is obtained as described in our previous work . Besides, the incoherent scattering cross sections can be exactly calculated according to the method described previously, and it was also calculated using the same PDOS and will not be discussed in the following parts.

Coherent Elastic Scattering Cross Section
We calculated the coherent elastic scattering cross sections of hydrogen in ZrH x at 300 K, based on the PDOS of δ-ZrH 1.5 . The coherent elastics scattering cross sections are shown in Figure 3. It can be seen that the coherent scattering cross section is larger than 1 b above the second Bragg energy. Compared with the incoherent scattering cross section, the value is negligible, so the TSL data for coherent scattering are not provided in ENDFs. The reason for the little coherent elastic scattering cross section is explained as follows. According the theory to calculate TSL data, the coherent and incoherent elastic scattering cross sections can be calculated as follows : where σ coh and σ inc are the characteristic-bound coherent cross section and incoherent cross section, which can be searched in the work by Sears (1992); S 0 (α, β) and S 0 s (α, β) are the scattering law, which can be calculated by the LEAPR module in NJOY and sab_calc module in NECP-Atlas.
For hydrogen, the bound coherent cross section is much less than the bound incoherent cross section as listed in Table 1,  which makes the final coherent scattering cross section much less than the incoherent scattering cross section. According to the work by Sears (1992), for most nuclides, the difference between the two bound scattering cross sections is very large. It seems reasonable to ignore the scattering mode with less bound cross sections. However, for some nuclides, the two bound cross sections are close, which will make the final coherent and incoherent scattering cross section comparable, for example, H-2, Li-6, and Li-7 listed in Table 1; Figure 4 shows the coherent and incoherent cross section of LiH. It can be seen that the coherent and incoherent elastic scattering cross sections of Li in LiH are comparable.

Inelastic Scattering Cross Sections
The previously described multipoint linearization method was implemented in the therm_calc module in this work. The incident energy points for ZrH x were generated with the multipoint linearization method, and the number of incident energy points between 1.0E-05 and 10 eV is 303. The inelastic scattering cross section of hydrogen in δ-ZrH 1.5 was calculated with the refined incident energy grid and compared with the results obtained using the fixed energy grid, as shown in Figure 5. The detailed cross section distribution for energies greater than 0.1 eV is given in Figure 6. The inelastic scattering cross section   shows oscillating along the energy. In the work by Whittemore (1964), the experimental and theoretical results revealed that the scattering cross section of hydrogen in ZrH x is oscillating, and the valley of each oscillation occurs at integer values of harmonic frequency hv. The value of hv for δ-ZrH 1.5 is 0.143 eV. In Figure 5, the inelastic scattering cross section is represented as a function of energy, whereas in Figure 6, it is represented as a function of the numbers of hv Figure 7 The refined grid obtained by multipoint linearization method captures more detailed variation of the cross section than the fixed grid. And the valley of the oscillation of the cross section appears at integers, which agrees with the experimental results by Whittemore (1964).
To show the effect of incident energy grid on the cross sections, we also analyzed the inelastic scattering cross sections of YH 2 , graphite, and H 2 O. Figure 8 shows the inelastic scattering cross section of YH 2 represented as a function of the numbers of hv. In the calculations, the TSL data of YH 2 are obtained from ENDF/B-VIII.0. For YH 2 , hv 0.119 eV. The inelastic scattering cross section distribution obtained using the refined incident energy grid captures more details of the oscillation. For graphite and H 2 O, the inelastic scattering cross sections vary smoothly with energy as shown in Figures 9, 10. The cross sections calculated based on fixed grid and refined grid are close to each other.

Results of Critical Benchmarks
To investigate the effects of the above treatments on the reactivity of reactors, a couple of critical benchmarks containing ZrH x as moderator was calculated using the Monte Carlo code NECP-MCX, which is a new Monte Carlo code developed by the Nuclear Engineering Computational Physics (NECP) Laboratory of Xi'an Jiaotong University. NECP-MCX is developed based on a hybrid Monte Carlo deterministic method, where the deterministic method is utilized to generate consistent mesh-based weightwindow and source-biasing parameters for the Monte Carlo method to reduce variance. The code has been verified against with various benchmarks Li et al., 2021).
The critical benchmarks are selected from ICSBEP benchmark (OECD-NEA, 2016) and listed in Table 2, including ICT003 benchmarks and HCM003 benchmarks. The ICT003 benchmark experiments were performed in a TRIGA Mark II reactor, which is a light-water reactor with an annular graphite reflector. The fuel in the reactor is a homogeneous mixture of uranium and ZrH x , with 12 wt% uranium of 20% enrichment. The temperature for all the materials in the benchmarks is 300 K. HCM003 benchmarks were performed on a reactor loaded with highly enriched uranium dioxide fuel (approximately 96% 235 U). The ZrH x is   December 2021 | Volume 9 | Article 779261 6 used as moderator. In the benchmarks, the temperature for all the materials is 300 K. The models for ICT003 and HCM003 benchmarks used by NECP-MCX were established according to the typical MCNP input given in the handbook of ICSBEP, without any simplification.
In our previous work, it was shown that the TSL data given in ENDF/B-VIII.0 and JEFF-3.3 could introduce larger errors into the reactivity of the TRIGA reactors, because the TSL is not obtained from a realistic crystal structure of ZrH x . Therefore, in this work, the TSL data calculated in the work by Zu et al. (2021) were adopted in the calculations: the thermal scattering cross sections obtained based on δ-ZrH 1.5 were adopted in ICT003 benchmarks; the thermal scattering cross sections obtained based on ε-ZrH 2 was adopted in HCM003 benchmarks. Except the TSL data, all the other nuclear data were extracted from the newly released CENDL-3.2. Ge et al. (2020) evaluated the nuclear data library.
In the thermal scattering library of ACE format, the flag IDPNC in the NXS array is used to indicate the elastic scattering mode for a material, IDPNC 4 for coherent elastic scattering and IDPNC 3 for incoherent elastic scattering. In order to use the two elastic scattering modes in the Monte Carol calculations, we extended the ACE library to include both the coherent and incoherent elastic scattering data (mixed elastic scattering) by setting IDPNC 5. The indices for these data were added to the JXS array in the ACE library. Besides, in the conventional Monte Carlo codes, only one elastic scattering model is sampled in the simulations. In this work, NECP-MCX was modified to simultaneously sample the coherent and incoherent elastic scattering.
In the calculations of the above benchmarks using NECP-MCX, the statistical uncertainties of the effective multiplication factor k eff were controlled within ±10 pcm. For ICT003, the calculations were run with 2,200 generations of 80,000 histories each, and the first 100 generations were excluded from statistics. For HCM003, 2,050 generations of 50,000 histories each were used, and the first 50 generations were excluded from statistics. Besides, in the calculations, only the scattering cross sections of ZrH x are generated with techniques mentioned previously, and for other materials, the scattering cross sections are generated using the conventional methods.
The k eff values calculated using different scattering cross sections are given in Table 3. The effect of the incident energy grid is given in the sixth column, which are values in the fourth column minus those in the third column. The effect of considering two elastic scattering modes is given in the last column, which is the value in the fifth column minus those in the fourth column. For HCM003 benchmarks, using the refined incident energy grid can reduce the k eff by a range from 111 to 141 pcm, and when the coherent elastic scattering is considered in the calculations, the k eff is further reduced by 38-66 pcm. Meanwhile, both the two factors make the k eff closer to the experimental results, whereas for the two ICT003 benchmarks, the refined incident energy grid gives a larger k eff of approximately 200 pcm than the fixed grid, and considering the coherent elastic scattering can predict a lesser k eff of approximately 30 pcm. Although it seems that the refined incident energy grid makes the k eff worse compared with the experimental, the uncertainties of experiment results of ICT003 are 560 pcm. The results of refined incident energy grid are still within the uncertainty range.
We also tested several assembly problems, including fuel pebble in HTR (She et al., 2021) and pressurized water reactor assembly benchmark VERA_2B (Godfrey, 2013). The results show that the incident energy grid has negligible effect on the k eff of these assemblies, because the inelastic scattering cross sections of graphite and H 2 O are smooth.

CONCLUSION
The treatments of thermal scattering cross sections are introduced in this article. The effects of ignoring one  elastic scattering mode in the evaluated nuclear data are analyzed using several critical benchmarks loaded with ZrH x . It is found that considering the coherent and incoherent elastic scattering simultaneously in the neutronics calculations can affect the effective multiplication factor by tens of pcm. The multipoint linearization method is adopted to refine the incident energy grid for inelastic scattering. The numerical results show that the incident energy grid has obvious effect on the effective multiplication factor.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.