Improvement of Geant4 Neutron-HP package: Unresolved Resonance Region description with Probability Tables
M. Zmeškal L. Thulliez P. Tamagno E. Dumonteil






1,2]M. Zmeškalcor1
marek.zmeskal@cvrez.cz
3]L. Thulliez
loic.thulliez@cea.fr
4]P. Tamagno
3]E. Dumonteil


[1]Research Centre Rez, Hlavni 130, Husinec-Rez, 250 68, Czech Republic
[2]Dept. of Nuclear Reactors, Faculty of Nuclear Sciences and Physical Engineering, Czech Technical University in Prague, V Holesovickach 2, Prague, 180 00, Czech Republic
[3]Université Paris-Saclay, CEA, Institut de Recherche sur les Lois Fondamentales de l'Univers, 91191, Gif-sur-Yvette, France
[4]CEA, DES, IRESNE, DER, SPRC, Cadarache, F-13108 Saint-Paul-Lez-Durance, France

[cor1]Corresponding author





Whether for shielding applications or for criticality safety studies, solving the neutron transport equation with good accuracy requires to take into account the resonant structure of cross sections in part of the Unresolved Resonance Region (URR). In this energy range even if the resonances can no longer be resolved experimentally, neglecting them can lead to significant numerical biases, namely in flux-based quantities. In Geant4, low energy neutrons are transported using evaluated nuclear data libraries handled by the Neutron High-Precision (Neutron-HP) package. In the version 11.01.p02 of the code, the URR can only be described by average smooth cross sections that do not take into account the statistical resonant structure of the cross sections. To overcome this shortcoming, the treatment of the URR with the use of the probability table method has been implemented in Geant4 and successfully validated with the reference Monte Carlo neutron transport codes (version 6.2) and (version 12). These developments will be taken into account in the next release of . All the validations of have been performed with probability tables generated from both the NJOY and CALENDF pre-processing tools. Therefore now has this unique feature to study the relative impact of the strategies involved during the production of probability table by the two pre-processing codes. This has been used to show that self-shielding is important also for inelastic cross sections in the example of ^238U. The tool to generate probability tables usable by either from  NJOY or from CALENDF is made available on a dedicated GitLab repository and will be included in Geant4. 







Neutron-HP Neutron Unresolved Resonance Region Probability tables NJOY CALENDF 










§ INTRODUCTION







Monte Carlo codes find widespread use in particle transport applications, enabling the individual tracking of particles and the preservation of correlations between them. To this purpose, has been developed by an international collaboration with the aim of supporting high-energy physics applications as well as other domains such as space science, medical physics, engineering, and nuclear physics <cit.>. In , the individual paths and reactions can be described by nuclear models or by evaluated nuclear data depending on the energy range of the particles. This second approach is used in its Neutron High-Precision (Neutron-HP) package describing low-energy neutrons, i.e. below 20 MeV, and in reference neutron transport codes such as <cit.> or <cit.>. With the increasing needs for neutron transport capabilities within the broader context of multi-particle transport codes, required by applications such as accelerator, medical and fundamental physics in addition to nuclear industry, the Neutron-HP package has recently been under huge improvements <cit.> bringing low-energy neutron transport physics almost as precise as reference neutron transport codes such as or . The last known drawback of the Neutron-HP package is its treatment of cross-sections in the Unresolved Resonance Region (URR), improved in this work by the introduction of probability tables.
For low energy neutrons (under 20 MeV), the isotope cross-sections are usually divided into 2 to 4 energy ranges as sketched in Figure <ref>, with energy boundaries varying from one isotope to another – in this case it is described for elastic cross section of ^238U. In the lowest energy part, say below a few eV called the thermal energy range the cross-section does not have any resonances, therefore it is a smooth function of energy. The cross section here is described either with the special thermal scattering law or together with the next region named Resolved Resonance Region (RRR). Then from a few eV to a few keV, the cross-section has huge amplitude variations because of the nuclear states of the compound nucleus formed by the neutron with the target. This is the Resolved Resonance Region (RRR), where the resonances are experimentally well observed, which allows to assign a set of resonance parameters (energy position, width, spin and parity) by fitting the resonance parameters with the help of the R-matrix formalism <cit.>. Above this energy, the resonances can no more be measured (resolved) because the distance between two resonances becomes small compared to the neutron beam experimental energy resolution and because the resonance widths are such that the resonances begin to overlap each other, this is the Unresolved Resonance Region (URR). In the following the URR is further defined as the region where self-shielding is important and needs special treatment in Monte-Carlo particle transport. There, the cross section is obtained from adjusted average resonance parameters. The ENDF files compile the  distribution and average of the partial widths and spacings of the resonances obtained from this adjustment. This region is also a transition between the R-Matrix theory describing the RRR and Hauser-Feshbach theory which describes the smooth cross-section above this region in continuum, in which the individual resonances fully overlap and are no longer distinguishable <cit.>. At an energy higher than the URR, the cross section is often described by an optical model potential.









The cross sections are often represented as point-wise tabulated data with prescribed interpolation laws satisfying a given precision reconstruction by a pre-processing code such as PREPRO <cit.> or NJOY <cit.>. In the URR, a smooth point-wise average cross section can be used as in – the only option available until now – or the so-called probability table (PT) method can be used allowing a better accounting of self-shielding effects <cit.>. Probability tables are usually obtained before the particle transport calculations, also using Monte-Carlo techniques, by aggregation of sampled resonant cross sections. In this latter case, a set of resonances characterizing the cross-section has to be defined. First the resonance energies, i.e. their spacing in a given energy range, are sampled from a Wigner distribution knowing the average resonance spacing provided in the ENDF file. Then for each sampled resonance energy, a resonance width is sampled in a χ^2 distribution with N degrees of freedom whose parameters are also defined in the ENDF file. This step is repeated for all the spins. Finally the PT is computed by combining all the sampled and reconstructed cross sections contributions. 
Then during a Monte Carlo simulation, when the neutron has an energy inside the URR, the corresponding PT is fetched and a cross-section is randomly sampled from it. 
Taking into account resonances in the URR by means of a PT is important for example for criticality studies since the use of PTs can increase the reactivity by up to 500 pcm for fast and medium systems with low enrichment due to self-shielding effect from ^238U <cit.> and also for shielding applications since the dose rate behind a ^238U shielding can increase by up to a factor 2.42 <cit.>. Thus omitting the statistical resonance structure in the URR could lead to non-conservative results. 
Until now, has only used smooth average cross sections to describe the URR. Therefore this paper reports the implementation of the URR within the Neutron-HP package. The benchmark methodology and the associated tools used in this work are detailed in Section <ref>. Then the different pre-processing codes, along with their specificities, used to build probability tables for are presented in Section <ref> along with a description of the technical implementation of the PT in . In Section <ref> the validation of the URR description by the PT method in is performed with the reference neutron transport codes and . Finally the influence of the choice of the NJOY or CALENDF pre-processing tool, used to generate the PT, on Monte Carlo simulation results is studied. 




§ BENCHMARK METHODOLOGY AND TOOLS



This study has been performed with version 11.01.p02. The  developments regarding the use of PT have been validated against the reference neutron transport codes  (version 12) and  (version 6.2). is developed since the mid-1990s at CEA-Saclay (France) and is used as a reference code by the main French nuclear companies <cit.>. is developed at Los Alamos National Laboratory (United States) and is used worldwide as a reference neutron transport code for many applications involving neutrons <cit.>. Both codes are used for criticality-safety studies, reactor physics and shielding applications. They benefit from a very large qualification database gathering more than 1000 experimental benchmarks and are regularly validated using inter-code comparisons <cit.>.
The evaluated nuclear data library ENDF/B-VII.1 <cit.> is used in this work.
In the following, the developments performed in are validated with a macroscopic benchmark that will be referred to as the sphere benchmark. It consists of a homogeneous sphere with a 50 cm radius from a given material – ^238U or tungsten isotopes – in which the neutron flux is recorded using a track length estimator. A 1 cm thick outer shell made of vacuum is used to record the outgoing flux, still with a track length estimator. An isotropic point-like neutron source is placed at its center. Two types of initial neutron energy spectra are used. The first one is a 1/E energy spectrum from 10 keV to 300 keV used for the ^238U isotope allowing to simulate a typical reactor spectrum. The second one is a uniformly distributed energy spectrum from 1 keV to 200 keV for the tungsten isotopes, mainly used in shielding applications. The energy ranges of the input energy spectra have been chosen to cover the whole URR range of the studied isotopes as shown in Table <ref>. This benchmark is an adaptation of the already existing benchmarks that can be found on the following GitLab repository <cit.>. Each simulation has been performed with 10^9 initial neutrons and normalized to one source neutron per second. The results of two codes are compared relative to three times the statistical uncertainties (3σ_stat).
Hereafter the validations of the use of PT generated by either NJOY <cit.> or CALENDF <cit.> are respectively performed with and . It has to be pointed out that for these benchmarks the exact same probability tables are always used. and are compared with PT from NJOY and and with PT from CALENDF, following the usual calculation scheme of both reference codes. It is worth noticing that is here the only code able to deal with both NJOY and CALENDF probability tables. 
In the Doppler broadening of the 0 K cross sections can be activated or deactivated in using the command line “/process/had/particle_hp/neglect_Doppler_broadening true/false”. If it is deactivated, uses the 0 K cross sections. When it is activated the Doppler broadening is performed stochastically on-the-fly, at any temperature, allowing to accommodate the different physics fields requirements. However this comes at a huge computational cost. For the physics cases dealt with in this work, this prevents to perform the simulations in a decent amount of time, i.e. in less than 30 days. Therefore in the following, a workaround has been used, which is now explained. Pre-Doppler broadened cross sections have been prepared at 294 K with the BROADR module of NJOY with a reconstruction tolerance of 0.1 %. After that, the 0 K cross sections have been replaced by the newly generated 294 K ones in the cross section database. Then in the simulation the on-the-fly Doppler broadening is deactivated resulting in the use of the pre-Doppler broadened cross sections at 294 K. This workaround leads to decrease of the computational time by at least a factor 100 for ^238U. It has to be mentioned that both and usually use pre-Doppler broadened cross sections. The Doppler Broadening Rejection Correction (DBRC) option is turned off in <cit.> and <cit.>, it is not available in .







§ PROBABILITY TABLES IN 






 §.§ NJOY and CALENDF processing tools


The two processing codes used in this work to generate probability tables are NJOY and CALENDF.   
uses probability tables generated by the PURR module of NJOY written in the ACE format, while uses the one from CALENDF. In evaluated nuclear data files, reactions are labeled by an MT number, for instance MT=1 for total, MT=2 for elastic scattering, etc. The related partial cross sections in CALENDF are however gathered for some usages as macro-reaction cross sections (MTREP). In this work, a special attention has been paid to the partial cross section grouping in CALENDF as presented in Table <ref> to agree with their Geant4 handling.  Although both, NJOY and CALENDF, codes create probability tables with the same approach of the sampling from ladders, there are some differences that could give rise to discrepancies, which are now discussed. 










Firstly, both codes generate PT for total (MT=1), elastic (MT=2), fission (MT=18) and radiative capture (MT=102) cross sections. While CALENDF also generates PT for the inelastic channels, if any, NJOY uses the smooth inelastic average cross sections.
Secondly, in CALENDF the user can define the energy bin limits and has to make sure that the energy grid is thin enough. Then PT are computed for each energy bin of the grid. Since one PT cross section corresponds to one realization, it will vary for different random generator seeds. This leads to randomly scattered cross section values around the average smooth cross section as can be seen in Figure <ref>. The average smooth cross section is recovered when applying a moving average algorithm as shown in Figure <ref>. On the other hand, NJOY computes PT in a much smaller number of energy points, at which the averaged cross section reconstructed from the PT and averaged smooth cross section are always equal. During the neutron transport simulation, for a given neutron energy, the corresponding PT cross sections are obtained by linear interpolation between the closest energy points at which PT have been computed. 
To get an accurate URR description a trade-off has to be done between a dense energy bin structure with less samples per bin (4 for ipreci=4 to 16 with ipreci=6) as in CALENDF and less energy points and more samples per point (minimum 16) as in NJOY. The difference between cross sections for different precision criteria in CALENDF can be seen in Figure <ref>. For ^184W in , Figure <ref> shows that the differences can be up to 20 % and that the neutron spectrum is smoother with the precision criteria ipreci=6 than with ipreci=4. The cross section calculated with higher precision is less scattered around the smooth one and so is closer to it. 
The influence of the chosen number of generated ladders in NJOY on the simulation results can be seen in Figure <ref>. Overall increasing the minimum value of 16 ladders to 32 and then to 64 ladders in NJOY shows that the differences are not greater than 5 % with the studied materials.
A tool to convert ENDF files to probability table libraries directly usable by has been written and made available to the community on GitLab <cit.>. The parameter settings used in this paper are now detailed for the two processing codes. In NJOY the cross section reconstruction tolerance is set to 0.1 % to Doppler broaden the cross section with the RECONR and BROADR modules and σ_0 to 10^10 barns for an infinite dilution. Moreover the number of probability bins is set to 20 and the number of resonance ladders to 16.
In CALENDF, MTREP have been defined as presented in Table <ref>. LCORSCT is set to false, meaning that the total cross-sections is the sum of the partial cross sections and LFORMRF is also set to false, meaning that the formalism to describe RRR is the one recommended by the evaluator, yet this should not impact the PT computed in the URR. 
The energy grid is the one from the example named “inu238e” from <cit.> having 11232 energy groups. 
The infinite dilution value is 10^10 barns as well and the precision parameter is set to ipreci=4 (see above discussion). The PTs for are then prepared from the .tpc file, i.e. after regrouping them as recommended in <cit.>, and all negative cross-sections are set to 0. This tool has been successfully validated in comparing its generated PT with the one used by as can be seen in Figure <ref>. The differences are larger than the statistical uncertainties but they are still less than 3 %, which can be explained by the different CALENDF settings to prepare the data. This tool has not been benchmarked against previous PT from NJOY since no PT data set is distributed with . Therefore this tool has been used in this work to prepare PT from NJOY, in ACE format, for both and . 














 §.§ Technical implementation in Neutron-HP package




For each reaction cross section _XXX_ ∈ {Elastic, Fission, Capture, Inelastic} two new classes named G4ParticleHP_XXX_DataPT and G4ParticleHP_XXX_URR have been written. These classes leverage the new class G4ParticleHPIsoProbabilityTable, which is created for each isotope and temperature. 
Depending on the PT data either from NJOY or CALENDF, the derived class from G4ParticleHPIsoProbabilityTable, respectively G4ParticleHPIsoProbabilityTable_NJOY or G4ParticleHPIsoProbabilityTable_CALENDF, are instantiated. In these classes the probability table file is read and the proper selection of the cross-section is performed. The most important feature to correctly handle probability table data during the simulation is to keep the same random number to sample the cross section in different spatial regions for a given isotope and for a given neutron trajectory having a given energy. An overview of the implementation of the new PT classes in the Neutron-HP package is presented in Figure <ref>.














§ VALIDATIONS OF PROBABILITY TABLES IN 



The validation procedure is made in three steps. Firstly, results for a single isotope sphere benchmark from , and are compared allowing to validate the good treatment of probability tables in . Secondly, a benchmark with a material having multiple isotopes are performed allowing to validate that the URR limits of all isotopes are correctly handled with respect to each other. Thirdly, a case where the sphere is made of one single isotope with different temperatures is studied to validate the correct handling of multiple temperatures. Theses different validation steps are now presented.





 §.§ Single isotope material


This first validation step consists in studying a sphere made of one single isotope allowing to validate the correct treatment of the PT. For that purpose a benchmark with ^238U and ^184W, without and with PT, are performed whose results are respectively presented in Figures <ref> and <ref>. 
 First of all the comparison of neutron fluxes between and or shows differences less than the statistical uncertainties. Therefore the probability table method implementation in is validated. It has to be noticed that the ^184W URR upper energy limit at 100 keV is clearly visible in most cases and gives rather non-physical results caused by the cross section definition in the ENDF file. 
Secondly, for the inside and escaping sphere fluxes, there are large differences between the use of smooth cross section and PT in the URR. 
For example, behind a uranium shielding the outgoing flux increases up to a factor 2 when using PT description of the URR compared to the use of smooth cross section as can be seen in Figure <ref>. These results are in agreement with the one from ref. <cit.>. The differences are even larger for the ^184W isotope case (Figure <ref>) since the outgoing neutron flux is around 15× higher with NJOY PT and 9× higher with CALENDF PT compared to the use of smooth cross sections.











 §.§ Multiple isotope material



This second validation step consists of studying a sphere made of a natural tungsten element, i.e. with the different tungsten isotopes with their natural abundances as presented in Table <ref>. This allows to validate that isotopes with different URR energy limits are correctly taken into account with respect to each other in the simulation. The results presented in Figure <ref> again show that the difference between simulations without and with PT are significant. and agree well with each other within the statistical uncertainties. This validates the correct handling of the different isotopes PT. However when comparing and , differences larger than the statistical uncertainties arise. This is clearly visible in the energy ranges 2.2-8.5 keV (corresponding to the lower limit of the single isotope URR range) and 45-100 keV (corresponding to the upper limit of the single isotope URR range) mainly for the flux inside the sphere, i.e. in the energy range where the cross section from point-wise and PT cross sections are combined together for the different isotopes.

Since the and multiple-isotope calculations agree with each other and that all the single tungsten isotope calculations agree with each other for both /and /as can be seen in the appendix respectively in Figures <ref> and <ref>, this indicates that the source of discrepancies could come from the isotope mixing in . This trend is also observed with another MC code named LAST <cit.> that uses the same libraries and input files as as can be seen in Figure <ref>. This is still under investigation.









 §.§ Same isotope with different temperatures


This third validation step consists of a hollow sphere made of a single isotope with 50 cm radius at 100 K containing an inner sphere of 25 cm radius at 10 K. Mixing the two temperatures allows to validate the correct handling of multiple temperatures for the same isotope. These two temperatures have been chosen to be low because in this case only the on-the-fly Doppler broadening of the 0 K cross-section in can be used with a limited computational cost. Indeed the workaround made in to use pre-Doppler broadened cross sections cannot be used here (cf. Section <ref>) since two temperatures are now involved. 
The neutron flux inside the whole sphere is presented in Figure <ref> and exhibits a difference larger than the statistical uncertainty below 20 keV for both cases with and without PT. This hints that it is due to the stochastic on-the-fly Doppler broadening in . In fact, as can be seen in Table <ref>, there are discrepancies up to 2.5 % between the value of the cross section obtained with the SIGMA1 algorithm of NJOY and the average value from the stochastic on-the-fly Doppler broadening algorithm of . Despite these effects, and results agree with each other in the URR energy range validating the correct handling of different temperatures of PT for the same isotope. No comparison has been made here with because in the official release there are no cross sections at 10 K and 100 K.








§ COMPARISONS BETWEEN IMPACT OF NJOY AND CALENDF PROBABILITY TABLES



can now use PT from both NJOY and CALENDF. Therefore it is now a powerful tool to study the impact of the different strategies used during the production of probability tables from URR parameters in each pre-processing code. To our knowledge, this is the first time that this is possible for a neutron transport code <cit.>.  
The comparison between fluxes using NJOY and CALENDF PT can be seen for the neutron flux inside the sphere of ^238U in Figure <ref> and in Figure <ref> for ^184W and tungsten element ^natW. 
For the three compositions, the choice of the pre-processing code used to generate the PT is shown to have a significant impact, up to 40% for the ^184W case, well above the statistical uncertainty. The differences for ^natW and ^238U are below 6 % for most energy bins. Additionally the more stochastic nature of the cross section in CALENDF's PT can be seen on the ^184W flux. A deeper investigation of the differences has been performed for the ^238U case.
Calculations have been performed using (blue curve) or not using (green curve) CALENDF's PT for the inelastic reaction. Using smooth inelastic cross section results are in a closer agreement – withing 2 % – with the simulation using NJOY's PT.
The remaining difference on the flux above 300 keV is due to fission neutrons (the initial neutron source spans from 10 keV to 300 keV (see Tab. <ref>).
This may indicate an impact of the PT related to fission which is to first order the main contributor to the fission rate.
Adding PT for the inelastic reaction leads to an overestimation above 45 keV and an underestimation below this energy which corresponds to the (n,n') energy threshold.
This tilt is expected, since in scattering reactions a neutron is moved from an energy to another, and so, even if the estimation of the scattering reaction rate changes, the slowing down process tends to preserve the amount of neutrons.
This shows that the self-shielding of inelastic cross sections is relevant.









§ CONCLUSION



The treatment of the Unresolved Resonance Region with the probability table method has been successfully implemented into . The validation performed with the help of a macroscopic benchmark with the reference neutron transport codes and shows very good agreement within the statistical uncertainties in all cases, except a small discrepancy which is still under investigation and could come from when mixing multiple isotopes. With this work can now be used to study the influence of the different strategies used in the pre-processing codes NJOY and CALENDF to generate probability tables which makes, to our knowledge, a unique tool.
This new treatment of the Unresolved Resonance Region with the probability table method will be incorporated in the next release. This makes the Neutron-HP package as precise as other reference neutron transport codes such as and in terms of physics description. The last drawback of regarding the neutron transport compared to others dedicated codes is its computational time which is still prohibitive, as underlined in this work. This mainly comes from its stochastic on-the-fly Doppler broadening method. In this work, a workaround has been found to use pre-Doppler broadened cross sections from the BROADR module of NJOY, allowing to overcome this last difficulty. It has also been shown that there are small differences between the cross section values obtained with the SIGMA1 and the on-the-fly Doppler broadening, which call for more investigations. Soon the use of pre-Doppler broadened cross sections at a given temperature in will be addressed to speed-up the simulations.



§ ACKNOWLEDGEMENT


This work was supported by the French government under the France 2030 program (P2I - Graduate School Physics) under reference ANR-11-IDEX-0003. The presented results were obtained using the CICRR infrastructure, which is financially supported by the Ministry of Education, Youth and Sports - project LM2023041.
The authors wish to sincerely thank Alberto Ribon head of the Hadronic Working Group for fruitful discussions regarding . 
 is a registered trademark of CEA. 














§ APPENDIX





unsrtnat





