Recent developments in the point wise neutron treatment for FLUKA v4

. The present paper describes the new, original implementation of a point-wise neutron treatment that is made available with FLUKA v4.3-0 (https://fluka.cern). The development is based on the pre-processed neutron data libraries for Geant4 code (https://cern.ch/geant4). The implementation was written from scratch in C++ and is able to load information from all available databases (ENDF, JEFF, JENDL, CENDL, BROND), including the materials with S(α,β,T) nuclear binding effects. It can perform online Doppler broadening at any temperature defined by the user. It carries out all interactions in a fully correlated way, conserving quantities like energy and momentum with the exception of S(  ,T). Special attention has been given to minimize the CPU penalty, which is typically associated with neutron point-wise interactions.


Introduction
FLUKA [1,2,3] is a multipurpose Monte Carlo code for the simulation of particle transport and interaction with matter. FLUKA relies on dedicated models for hadronic interactions, with the notable exception of low-energy neutrons (below 20 MeV), for which a group-wise treatment based on evaluated nuclear data is employed for a sizable list of relevant isotopes. A limited pointwise treatment was made available only for a few isotopes and channels of practical relevance. The advantages of the group-wise neutron treatment are that it requires a very modest memory footprint and it is considerably faster than other methods, rendering this approach optimal for shielding calculations and deeppenetration problems. However, it suffers from several drawbacks, including insufficient precision in resonance regions, self-shielding issues, and a too coarse discretization of the energy and angular distribution of emitted secondaries. These drawbacks unfortunately limit the use of group-wise neutron transport for several applications.
In this work, a revised point-wise treatment for lowenergy neutron transport in FLUKA v4 (distributed by CERN) is described, relying on the processed evaluated neutron cross section libraries that are used by Geant4 [4,5,6]. Special attention is devoted to the implementation decisions adopted in order to ensure energy and momentum conservation at every interaction, to respect the supplied distributions of emitted particles, preserving correlations and avoiding discretizations to the extent possible, while trying to keep the implementation as performant as possible. * Corresponding author: Vasilis.Vlachoudis@cern.ch For performance reasons our implementation follows a hierarchical cross section grouping of the lowenergy neutron interactions. All cross sections are Doppler broadened [7] during initialization at the arbitrary temperatures requested by the user. For thermal neutron energies (En < 4 eV) the implementation accounts for chemical-binding effects. Molecular and crystal-lattice effects are accounted for via an S(α,β,T) treatment, including coherent and incoherent effects for a few selected materials. Finally, in view of the ongoing evolution of the FLUKA code, the implementation presented here is written as a C++ module that can be easily linked with other codes. Benchmarks against experimental data and other transport codes are performed.

Cross Sections
All-inclusive cross section data originally pre-processed with PREPRO for the Geant4 are freely available from IAEA [4,5,6]. The distributed cross sections are produced at 0 K starting from the ENDF-6 formatted databases using a linear interpolation with a maximum error of 0.1%.
Individual channel cross sections are read from the data files, while total cross sections (per natural element and per isotope) are calculated on initialization. The four main channels are i) Elastic Scattering, ii) Radiative Capture, iii) Inelastic (Reaction), iv) Fission. Below 4 eV, upon user demand, molecular/crystal binding effects can be accounted for by way of an S(α,β,T) treatment for selected materials. This treatment is composed of three channels: • inelastic scattering • coherent elastic scattering • incoherent elastic scattering The term inelastic refers here to the collective excitation of the scattering medium, not of an individual target nucleus. The nuclear inelastic reaction channel is further subdivided into 36 sub-channels. For some isotopes, the fission channel is further sub-divided into first, second, third, and fourth chance fission.
All cross sections are Doppler broadened during the initialization phase using the exact Doppler broadening formula, Eq.(1), with an adaptive numerical integration preserving the initial accuracy of the data Numerical integration is sufficiently fast: it is able to process all ~550 isotopes available in the databases in less than 20s on a modern 3-GHz single core computer. Fig.1 displays the 113 Cd(n,) cross section, Doppler-broadened to 296 K, yielding a curve (black) that is indistinguishable from that of JEFF 3.3 (green). The innovative cross section indexing algorithm (described below), together with the aforementioned hierarchical cross-section structure, greatly enhances the searching and interpolation of the neutron energy and cross section tuple, giving searching times equivalent to those for unionized-energy grid cross sections, with minimum memory requirements.

Indexing algorithm
When processing neutron cross-sections using a linear interpolation (requesting a fixed maximum inaccuracy of the result), an almost logarithmic energy grid is obtained. The resonance region requires higher density of energy points with respect to the lower (thermal) and higher energies. For the requested accuracy, some channels exhibit hundreds of thousands of pairs (E,σ). One of the worst cases is 238 U, with about 600'000 energy pairs). Binary searching will follow a log2 convergence, compromised by the delay in the almost random pre-fetching of the data from memory. To address this problem, we use an indexing table per channel that restricts the search in a small number of data pairs. It not only reduces the requested binary search steps, but also permits to pre-fetch the entire data subsection in the memory cache reducing the latency time from possible cache-misses. Profiling the algorithm revealed a significant gain in search time for all isotopes, almost up to a factor of 10 for some large data size cases like 238 U.
The indexing algorithm is based on an integer approximation of log2,, exploiting the IEEE754 format of floating point numbers. A single-precision (FP32) floating-point positive number is expressed as: = 2 −127 (1 + ) (4) where E is the exponent in base 2, and 0 ≤ < 1 is the mantissa. Taking the log2 at either side one obtains log 2 ( ) = − 127 + log 2 (1 + ) (5) which, can be approximated to first order in as log 2 ( ) ≅ + − 127 (6) Thus, performing a few of bit-shift operations, integer subtraction, and logical OR operations, we can create a unique integer index for each energy based on the top-most significant bits from the exponent representation and the remaining top-most bits from the mantissa, forcing it to start close to zero for the minimum available energy per channel (Fig.2). This integer index follows log2(E) representation of the energy which can be used to retrieve the tabular index. When a cross section is required for a specific energy E, the integer index is calculated, which is used to bracket the binary search in a reduced number of points from the lower point energy point (Ei,σi) based on the index up to the next indexed pair (Ei+1,σi+1). For optimum performance versus data size of the indexing table, a 12 significant bits = 4096 entries index was used, requiring an indexing table up to a maximum size of 16kb per channel.

Interactions
All interactions with the exception of S(α,β,T) in the novel FLUKA pointwise treatment are performed in a fully correlated way, respecting as much as possible the energy and angular distributions (if any), of all ejectiles, while conserving physical quantities like energy and momentum at every interaction. The interaction is performed by treating the N-body final state as subsequent 2-body emissions, where the kinematics gradually constrains the database distributions. Consequently, the first ejectile follows the database distributions, while subsequent ejectiles are slightly biased towards lower energies. The order in which ejectiles are emitted respects the order prescribed in the datasets, where the neutron is preferentially emitted first. All interactions are performed through a compound nucleus, which in the end might be left in an excited state (either a known level or in the continuum), promptly de-excited via FLUKA's nuclear evaporation model.

Elastic Scattering
Elastic scattering is the most relevant process for neutron thermalization. The attribute elastic implies that the neutron is scattered from the target nucleus without leaving the latter in an excited state. For high neutron energies above the limit of several kT (typically a few eV at room temperature), the target nucleus motion is negligible with respect to the neutron speed and, consequently, the code assumes the target nucleus is at rest. Instead, at lower energies, thermal scattering is applied (see below).
For this process, the polar scattering angle is sampled in the center of mass frame, provided by the databases as a function of the incoming neutron energy.

Thermal Scattering
Below the limit of a few eV, the target nucleus is no longer assumed to be at rest: its thermal motion has to be taken into account. There are two approaches implemented: Free gas approximation and molecular binding via a thermal scattering law S(α,β,T), especially relevant for light target nuclei,, for which the chemical or lattice binding is taken into account for neutrons below 4 eV.

Free gas approximation
In the free gas approximation, the target nucleus motion is sampled from the Maxwell-Boltzmann distribution. Sampling a target nucleus velocity is effectively altering the neutron energy in the rest frame of the target nucleus: this change in energy translates into a different cross section value, and hence an interaction rate. Moreover, it can arrive to kinematically impossible situations when the target nucleus is moving faster than the incoming neutron. There is a set of corrections in the code for this purpose, assuming a constant cross section model, which in general is a good approximation, except when strong resonances are present in the thermal region.

Molecular Binding S(α,β,T)
At very low neutron energies, when the neutron energies are comparable to the chemical binding energies or when the neutron wavelength becomes comparable to interatomic distances, the free gas approximation is no longer valid. Results can be highly inaccurate if these effects are not considered (see e.g. Fig.3). The effect of the local binding environment (molecular or crystalline) can be accounted for by means of a thermal scattering function S(α,β,T), where α and β are the respective momentum and energy transfer during the collision, divided by kT. The S(α,β,T) thermal scattering function is often divided into its inelastic and elastic parts, differing on the state in which the scattering system is left after the interaction.
In inelastic thermal scattering, the neutron energy can be either increased or reduced. Scattering data are encoded in table of (Ein, Eout, cosθ), where the outgoing energy Eout is represented with equiprobable energy bins. To avoid possible aliasing effects due to the presence of a sharp peak at Ein=Eout, the code is applying an energy rescaling / correction Instead, in elastic thermal scattering the neutron energy remains the same while the direction is changing. This channel is further divided into coherent scattering (e.g. on crystalline structures with the explicit presence of a Bragg edge) and incoherent scattering (typically in Hydrogenous materials).

Radiative Capture
In radiative capture (n,γ), the neutron is absorbed by the target, forming the compound nucleus which deexcites via the standard FLUKA evaporation model.

Inelastic Reaction
The inelastic channel in the ENDF databases is subdivided into 36 sub-channels, grouped in the following categories: • the composite channels (n,n'), (n,p), (n,d), (n,t), (n, 3 He), (n,α), where there is only one product (with the exception of photons). These channels additionally contain the the probability to populate a discrete level in residual nucleus (first 49 levels) or to populate the continuum (above the 50 th level). • all other known channels with multiple products (with the exception of photons) The various energy-and-angular distributions of emitted products are described with a plethora of different laws (e.g. Legendre expansion, tabulated distributions, etc.) in the ENDF format (files MF=4,5,6), often as a function of the incoming neutron energy. FLUKA is able to read and sample from all of them, using the information as stored in the ENDF databases without the need for further pre-processing. Most of the energy and angular sampling for emitted ejectiles is performed by rejection, using tight envelopes to improve the sampling efficiency.

Fission
The total number of fission neutrons is sampled from a Poisson distribution with a mean equal to the neutron multiplicity provided in the databases. These neutrons are further divided into pre-and post-scission neutrons as follows. Pre-scission neutrons are emitted from the compound nucleus, either using the information provided by the fission sub-channels if available (first chance (n,f), second (n,nf), third (n,2nf) and fourth chance (n,3f) fission), otherwise about ~10% [8] of the total fission neutrons are treated as pre-scission neutrons (if their emission is kinematically allowed). The remaining nucleus scissions into two fragments. The fragment mass and charge are sampled from the database distributions if available, otherwise the Wahl [9] systematics is employed. Fission fragments further emit a number of post-scission neutrons proportional to their neutron excess The final nuclear fragments are forced to de-excite using the photon distribution provided by the databases.

Nuclear recoils
All neutron interactions, with the exception of thermal scattering via the S(α,β,T) treatment described above, close the kinematics and push the recoiling residual nucleus to the main FLUKA particle stack. In the aftermath of an elastic scattering event, if the recoil has less than a few eV, its energy is deposited on the spot, in order to avoid flooding the main stack.

Benchmarking
An extensive benchmarking campaign took place within the FLUKAVAL framework [2] in order to validate the performances of the novel implementation reported here, both in terms of physics and of CPU time. The validation notably included single interaction checks on all isotopes in the available databases, as well as more complex setups and, if possible, comparisons with the previous group-wise implementation and the predictions of other simulation codes (for whom the license permits it). Overall, a remarkable agreement of the single interaction observables was obtained for all the databases. The present point-wise interaction scheme exhibits a modest slowdown with respect to the previous group-wise treatment, of the order of 5-30% for typical high-energy shielding applications (Fig.4), and up to a factor of 2-3 for applications where only neutron interactions are performed.

Summary
The newly implemented point-wise treatment in FLUKA offers a reliable and fast way of exploiting the information from all available neutron databases. The algorithm was heavily tested, yielding remarkable performances, both in terms of CPU time and of physics results. Special emphasis is put on the correlation among interaction products, preserving energy and momentum at each interaction (note, however, that this implies a small distortion of the final distributions). The fully correlated treatment allows one to perform event-byevent studies needed for detector simulations, thus overcoming an intrinsic limitation of the previous group-wise interaction treatment that FLUKA relied on.