Generation of thermal scattering files with the CINEL code

. The CINEL code dedicated to generate the thermal neutron scattering files in ENDF-6 format for solid crystalline, free gas materials and liquid water is presented. Compared to the LEAPR module of the NJOY code, CINEL is able to calculate the coherent and incoherent elastic scattering cross sections for any solid crystalline materials. Specific material properties such as anharmonicity and texture can be taken into account in CINEL. The calculation of the thermal scattering laws can be accelerated by using graphics processing unit (GPU), which enables to remove the short collision time approximation for large values of momentum transfer. CINEL is able to generate automatically the grids of dimensionless momentum and energy transfers. The Sampling the Velocity of the Target nucleus (SVT) algorithm capable of determining the scattered neutron distributions is implemented in CINEL. The obtained distributions for free target nuclei such as hydrogen and oxygen are in good agreement with analytical results and Monte-Carlo simulations when incident neutron energies are above a few eV. The introduction of the effective temperature and the rejection step to the SVT algorithm shows improvements to the neutron up-scattering treatment of hydrogen bound in liquid water.


Introduction
Neutron transport through materials is crucial for nuclear reactor core design, neutron detection and shielding. If the incident neutron energy is much larger than the atomic binding energy, the target materials can be approximated by a collection of unbounded atoms. In this case, the treatment of the neutron scattering with materials can be well approximated by the two-body collision kinetics. However, when neutron energies lie below a few eV, the atomic vibration behaviors and/or crystalline structure play a predominant role in the neutron scattering reactions. Low-energy neutrons have been used to probe the crystalline structure of oxides of uranium [1][2][3][4] and other materials [5,6] since the sixties.
The low energy (below a few eV) neutron scattering theory is well detailed in the literature [7][8][9]. The key quantity is the double differential neutron scattering cross section (DDXS) which represents the probability of finding a scattered neutron with a given energy at a specific solid angle. The DDXS has been expressed as a function of the dynamic structure factor S( Q, ω) by Van Hove [10] in 1950s, where Q and ω are respectively the momentum transfer and energy transfer of neutrons. Under the incoherent and Gaussian approximations, the * e-mail: gilles.noguere@cea.fr dynamic structure factor S( Q, ω) can be calculated from the phonon density of states (PDOS) of materials. The PDOS representing the structure dynamics of materials can be obtained from ab initio density functional theory (DFT) or molecular dynamics (MD) calculations. S( Q, ω) are later converted to dimensionless thermal scattering laws (TSL) or S(α, β). S(α, β) are stored in ASCII files by following the ENDF-6 format requirements [11], to prepare the neutron scattering data libraries.
The TSL in ENDF format have been generated by using the LEAPR module [12] of the nuclear data processing code NJOY [13] which is an industry standard of the same class codes. Stemming from the incoherent and Gaussian approximations, LEAPR is able to calculate the TSL for various kinds of materials. For solid crystalline materials, S(α, β) are further divided in an elastic and inelastic parts, with coherent and incoherent terms. However, the calculations of the coherent elastic scattering cross sections in LEAPR are limited to crystalline materials with specific symmetries. In addition, to reduce the computational time and/or for coding convenience, LEAPR adopts the short collision time (SCT) approximation and the atom site approximation in which the mean-squared displacement (MSD) are supposed the same for different types of atoms in the same material.
To overcome the limitations of the LEAPR module, new codes have been developed in the scientific community: NJOY+NCrystal [14], FLASSH [15] and OCLIMAX [16,17]. They are able to calculate the coherent elastic scattering cross sections of any kind of solid crystalline materials. NJOY+NCrystal allows to generate the coherent and incoherent elastic cross sections in the same ENDF file (mixed elastic format [18]). FLASSH and OCLIMAX improve the calculations of TSL by introducing the coherent one-phonon correction [8] to the incoherent approximation and eliminating the cubic approximation.
In this work we present the CINEL code which is dedicated to generate the thermal neutron scattering files in ENDF-6 format for solid crystalline, free gas materials and liquid water. Compared to the LEAPR module of the NJOY code, CINEL is able to compute the coherent elastic scattering cross sections of any crystals. The generated ENDF files support the mixed elastic format. CINEL is able to generate automatically the α and β grids used for the calculations of S(α, β). Specific properties such as anharmonicity and texture can be taken into account in CINEL. The SCT approximation is removed thanks to the use of graphics processing unit (GPU) acceleration. CINEL is able to use the PDOS obtained by MD or DFT calculations.
The main expressions involved in the low-energy neutron scattering formalism are reported in Section 2. The implementation of CINEL is briefly presented in Section 3. CINEL is composed of three modules named Cubic, INELastic and SVT, respectively. Validations of Cubic and INELastic by comparing with other codes and databases are presented in Section 4. The automatic generation of α and β grids, mixed elastic cross section, taking into account the anharmonicity and texture are illustrated with examples. Section 5 presents the validations of the SVT module by comparing with Monte-Carlo calculations and analytical results. Conclusions and perspectives are given in the last section.

General formalisms
As mentioned above, the key quantity in low-energy neutron scattering is the double differential cross section (DDXS). For unpolarised neutrons and samples, based on Born approximation or Fermi's Golden rule [7], the DDXS can be obtained by: where k i and k f are respectively the wave vectors of incident and scattered neutrons, k i and k f represent the wavenumbers with k = | k|, Q is the momentum transfer of neutron defined by Q ≡ ( k i − k f ) in which is the Planck constant divided by 2π, ω is the energy transfer with ω ≡ E i − E f 1 , and S( Q, ω) is the dynamic structure 1 E i represents the incident neutron energy, the subscript i is omitted in this work for the sake of clarity. factor or scattering function defined by: In equation (2), N is the number of particles in the scattering system under consideration, b j represents the scattering length of the particle j, b j is the average value, A represents the operator expectation value in the scattering system, and j , j (notation taken from Ref. [9], which is referred as intermediate function in Ref. [8]), represents the correlation between the position of the particle j at time t and the position of the particle j at time 0: where R(t) is the time-dependent Heisenberg operator. It is assumed that there is no correlation between the scattering lengths and the positions of the particles in the scattering system [7], then Equation (4) enables to decompose the dynamic structure factor in equation (2) into two distinct parts: where and (7) S coh ( Q, ω) is the coherent scattering function, representing the correlation between the positions of the same particle at different time (when j = j ) and the correlation between the positions of different particles at different time (when j = j ). Therefore, S coh ( Q, ω) gives interference effects [8]. S inc ( Q, ω) is the incoherent scattering function, which represents the correlation between the positions of the same particle at different time. The incoherent scattering does not give interference effects (there is no terms j , j for j = j in Eq. (7)).
Based on equation (2), the value of the scattering function S( Q, ω) increases with the number of particles in the scattering system. In practice, the normalization is performed with respect to the number of particles given by the chemical formula or the unit cell [7]. In addition, the principle of detailed balance must be fulfilled for the dynamic structure factor S( Q, ω) [7]: where k B is the Boltzmann constant, T is the temperature of the material under consideration.

Harmonic, incoherent and cubic approximations
For solid crystalline materials, the operator representing the position of atom R(t) in the intermediate function (Eq. (3)) can be represented as the displacement of the atom u(t) from its equilibrium position d at time t: In addition, the equilibrium position d is a simple vector, which commutes with all the operators. Hence, the intermediate function can be obtained by: For practical purposes, it is desired to permute the exponential and the operator expectation in exp(−i Q · u j (0)) exp(i Q · u j (t)) in equation (10). To this end, the harmonic approximation is supposed, in which atoms are assumed to have harmonic vibrations, i.e., the interatomic forces are linear with respect to the displacements of atoms and all the higher terms related to the anharmonic vibrations are neglected. The intermediate function j , j in equation (10) becomes: where W j ( Q) is the Debye-Waller function, which measures the mean-squared displacement (MSD) of the atom j along the direction Q: Note that corrections related to the anharmonic effects can be lately introduced to the intermediate function in equation (11), in the case of small atomic displacements compared to interatomic distances.
In equation (11), ( Q · u j (0))( Q · u j (t)) represents the correlation of linear displacement along Q of two atoms at different time. Though the correlation for distinct atoms j = j is well formulated, the calculation is computationally cumbersome. Thus, in this work, a widely used approximation is utilized, called incoherent approximation, in which the correlation for distinct atoms is neglected for inelastic scattering, i.e., For j = j , the last term in equation (11) does not vanish. Its computation is quite difficult ( ( Q · u j (0))( Q · u j (t)) ). To circumvent this difficulty, the cubic approximation is used in this work. It assumes that for atoms of a solid crystalline material under consideration, the interatomic forces along all directions are isotropic. Under the cubic approximation, we have: where M j is the mass of the atom j, ρ j (ω) represents the phonon density of states (PDOS) of the atom j satisfying The cubic approximation has the considerable merit that it enables to connect the PDOS of atoms to their dynamics. By utilizing the symmetry of the PDOS versus ω, i.e., ρ j (−ω) = ρ j (ω), equation (14) can be rewritten as: An important condition emerges from equations (16) and (17): ρ j (ω) must vary as ω 2 as ω goes to zero, to ensure the convergence of the temperature-dependent function P j (ω). The integral of P j (ω) provokes the Debye-Waller coefficient Λ j (T ) [13]: Based on equations (12) and (16), the Debye-Waller function W j ( Q) represents the dynamics of the atom j at t = 0. In addition to equation (18), W j ( Q) can be obtained from the corresponding Debye-Waller coefficient Λ j (T ) as: Thanks to the incoherent and cubic approximations, the dynamic term ( Q · u j (0))( Q · u j (t)) and the Debye-Waller function W j ( Q) in the intermediate function (Eq. (11)) can be obtained with equations (16) and (19), respectively. The following step consists of expanding exp( ( Q · u j (0))( Q · u j (t)) ) in equation (11) in a Taylor series: Equation (20) is referred to phonon expansion, originally introduced by Sjölander [19]. The first term n = 0 corresponds to elastic scattering (no energy exchange between the incident neutron and the system under consideration). The terms n ≥ 1 represent inelastic scattering contributions.
The cross sections of the coherent elastic scattering, incoherent elastic scattering and inelastic scattering are presented in the following sections.

Coherent elastic scattering
Important crystallographic quantities used for the calculations of the coherent elastic scattering cross section for solid crystalline materials are presented. Any crystal can be characterized by a periodic unit cell containing N uc atoms. In the direct lattice, it is defined by a set of unit vectors a, b, c of lengths a, b, c and volume V uc = a · ( b × c). The angles between them are conventionally denoted α, β and γ. The position of the jth atom located at the point (x j , y j , z j ) is given by: In the reciprocal lattice, the unit vectors become: and indices h, k, l denote planes (hkl) orthogonal to the reciprocal lattice vector: Given that the Fourier transform of a constant function is a Dirac delta function: and thanks to equations (6), (11) and (20), the coherent elastic dynamic structure factor S el coh ( Q, ω) is given by: Considering a crystal with the previous notations, S el coh ( Q, ω) can be rewritten as: The neutron diffraction will occur from the planes (hkl) that are oriented at the correct angle to fulfill the Bragg condition. The coherent elastic cross section per atom emerges from the sum of all the neutron scattering contributions over the N uc atoms of the unit cell and plans (hkl). The analytical expression taking into account the anharmonicity and the texture properties for crystal powder is given by: in which m is the neutron mass, E hkl = 2 τ 2 hkl /(8m) represents the Bragg edges, d hkl = 2π/τ hkl stands for the distance between adjacent planes (hkl) and b j is the bound coherent scattering length.
Compared to the reported equation in the literature [9,14], in our work, the texture or preferred orientation correction is taken into account by introducing the cylindrically symmetric Pole-Density Distribution Function (PDDF) P hkl (Θ hkl ) which depends on the orientation angle Θ hkl , i.e., angle between the preferred orientation vector PO and the plan vectors τ hkl : Various types of PDDFs reported in the literature [20] are investigated, such as the models proposed by March-Dollase [21,22], Altomare et al. [23],Černý, Valvoda, and Chládek [24]. The impact of introducing the PDDFs to the calculations of σ el coh (E) will be illustrated in Section 4.4. The third-cumulant coefficient c j 123 is introduced in the form factor to account for the anharmonicity of atoms [25]. The influence of the anharmonicity of oxygen atoms in uranium dioxide on σ el coh (E) will be investigated in Section 4.3.

Incoherent elastic scattering
S el inc, j (Q, ω) representing the incoherent elastic dynamic structure factor of the atom j is obtained thanks to equations (7), (11), (19), (20) and (24), which is given by: The incoherent elastic scattering cross section σ el inc, j (E) is obtained by integrating the dynamic structure factor S el inc, j (Q, ω) in energy and scattering angle. The analytical result is given by: where A j = M j /m is ratio of the mass of the atom j to the neutron mass and σ inc, j = 4π b 2 j − (b j ) 2 is the bound incoherent scattering cross section of the atom j. The incoherent elastic cross section can thus be determined by the Debye-Waller factor Λ j (T ).

Inelastic scattering
Under the incoherent and cubic approximations, the inelastic dynamic structure factor S inel (Q, ω) can be obtained by summing the coherent and incoherent inelastic components. Thanks to equations (6), (7), (11) and (20), S inel j (Q, ω) representing the inelastic dynamic structure factor of the atom j is given by: where T j, n (ω) has the generic form (when n ≥ 2): The one-phonon term is given by:

Scattering models for free gas materials
The momentum transfer Q and the energy transfer ω are large when neutrons scatter with free gas particles. Hence, the interaction time is short and the time dependent Heisenberg operator R(t) which represents the position of the particle in the intermediate function (Eq. (3)) can be approximated by its first order term [26]: where v is the speed of the particle.
Since there is no inter-particle interaction for a free gas material, the incoherent approximation is automatically valid. Thus, the intermediate function j, j is given by: The free gas material is assumed to be monatomic. From the statistical point of view, the speeds of the free gas atoms v follow the Maxwell- where v 2 p = 2k B T /M is the square of the most probable speed, T represents the temperature of the gas and M is the mass of the gas atom.
By replacing the M-B distribution (Eq. (37)) in equation (36), thanks to the Gaussian integral property and the Fourier transform of the Gaussian function, the intermediate function for monatomic free gas is given by: where w(t) is the width function of the free gas: (39) should be corrected to account for quantum effects [27,28]: Based on equations (2), (38), (39) and (40), the dynamic structure factor for monatomic free gas can be obtained by:

Scattering models for liquid water
The intermediate function j , j (Eq. (3)) of liquid materials obtained by a Fourier transform of the generalized pair distribution function G( r, t) is introduced by Van Hove [10].
By distinguishing the contribution from the same particle to different particles, Van Hove split G( r, t) into a self part G s ( r, t) and a distinct part G d ( r, t): For a given particle at position r and time t, G s ( r, t)d rdt represents the probability of finding the same particle between r and r + d r in a time interval t and t + dt. Similarly, G d ( r, t)d rdt is the probability of finding a distinct particle between r and r + d r in a time interval t and t + dt.
Since it is difficult to predict theoretically the contribution from the distinct part, we use the incoherent approximation, as for solid crystalline materials: To obtain an analytical form of G s ( r, t), Vineyard [29] introduces the Gaussian approximation that assumes a weak time-dependent atomic position coupling: where 3w 2 (t) represents the mean square displacement of the particle after time t, in which w(t) is called the width function.
Under the incoherent and Gaussian approximations, the dynamic structure factor S( Q, ω) for liquid materials is given by: The Q 2 dependence of equation (45) indicates that S( Q, ω) is isotropic with respect to the momentum transfer, leading to: The width function w(t) plays a central role in the calculation of the dynamic structure factor. In this work, we focus on liquid water which is widely used as moderator in pressurized water reactors. An analytical expression of w(t) was proposed in the sixties [30] for liquid water. However, introducing such an expression in equation (45) leads to numerical issues at zero energy transfer. These mathematical problems can be solved by decoupling the mean square displacement into a diffusion part w diff (t) and a solid-like vibrational part w vib (t): The corresponding diffusion S diff (Q, ω) and solid-like S vib (Q, ω) parts of the dynamic structure factor are calculated with equation (45) by introducing different models for w diff (t) and w vib (t). According to the convolution theorem for Fourier transforms, the composite dynamic structure factor is calculated as follows: where /b 2 is a scale coefficient, the values of S diff (Q, ω) and S vib (Q, ω) depend on the density of states ρ(ω) of liquid water, which is partitioned in two parts: The diffusion w d and solid-like vibrational w v weights satisfy the condition: The diffusion weight w d can be calculated with the model involving the fluidicity [31,32], proposed by Marquez Damian (private communications, 2021). The roto-translational diffusion model with random jump diffusion correction for S diff (Q, ω) has been investigated in reference [33]. The formalisms are recalled here.
In view of using a descriptive analytical model, the diffusion term in equation (48) related to liquid water is approximated by a translational diffusion motion corrected for rotational diffusion contributions: The translational diffusion behavior of the water molecule was interpreted with the Egelstaff-Schofield diffusion model [34], in which the width function is given by: where M is the mass of the scattering atom and c is a dimensionless diffusion parameter. Its meaning has been revisited in reference [35] by introducing the Singwi-Sjolander residence time τ 0 [36] to account for the non-continuous motions of the water molecules: with in which D represents the self-diffusion coefficient of the liquid water. The dynamic structure factor S trans (Q, ω) can be derived analytically by introducing equation (52) in equation (45) with correction to fulfill the principle of detailed balance: in which K 1 (x) is the modified Bessel function of the second kind, with: The rotational diffusion correction in equation (51) is related to the dynamics of hydrogen bonds. Their contributions is approximated with the Sears expansion [37] expressed in terms of spherical Bessel functions j l (QR cm ) of order l: This model assumes a continuous and isotropic rotation of the hydrogen atoms arround the center of mass of the molecule with a rotational relaxation time τ R . The rotation of the water molecule is confined to a spherical surface of radius R cm , whose value is often taken equal to the O-H intramolecular distance.
The computation of the solid-like S vib (Q, ω) part of the dynamic structure factor is equivalent to the inelastic part of dynamic structure factor for the solid crystalline materials.

Converting S(Q, ω) into thermal scattering laws
The dynamic structure factor for nuclear materials (e.g., beryllium, iron, liquid water) is routinely used to simulate the neutron behaviors in light water reactors. The neutron transport formalism implemented in dedicated Monte-Carlo or deterministic code systems relies on the symmetric form of the thermal scattering laws S(α, β) [12,38], which is defined as a function of the dimensionless parameters α: and β: where µ = cos(θ) is the cosine of the scattering angle θ in the laboratory system and A is the ratio of the mass M of the scattering atom to the neutron mass. The change of sign in equation (59) is introduced to keep neutron energy gains positive in neutron transport calculations. The relationship between S(α, β) and S(Q, ω) is given by the following expression: Note that the thermal scattering laws S(α, β) are symmetric in energy.

Implementation of CINEL
CINEL is developed in Python by using JupyterLab [39]. JupyterLab is a new generation of user interface [40] which is an open-source tool enabling to mix live code, texts, mathematical equations and interactive graphics. An example of JupyterLab interface is given in Figure 1 for illustration. This interactive development environment facilitates the verification and visualization of calculations.
CINEL consists of three modules: Cubic, INELastic and SVT. The data flow presenting the generations of the thermal scattering laws (TSL) in ENDF-6 format by using the Cubic and INELastic modules, and the scattered neutron energy distributions via the SVT module are shown in Figure 2. These modules will be briefly presented in the following.

Cubic and INELastic modules
Various physical quantities are necessary for the generation of the TSL such as scattering length, positions of atom in the unit cell, general settings and temperaturedependent parameters including the lattice parameters, PDOS and the phonon expansion order. Optional parameters are preferred orientation of crystal PO, cylindrically symmetric PDDF P hkl (Θ hkl ) and the third-cumulant coefficient c O 123 . The Cubic module calculates the Debye-Waller coefficients Λ(T ) and the coherent elastic scattering cross sections σ el coh (E). Λ(T ) is computed from the PDOS (Eq. (18)) since the cubic approximation is assumed. σ el coh (E) is calculated by directly looping through the  indexes hkl until the interplanar distance d hkl of the corresponding plan is smaller than 0.1Å. The default threshold value is chosen to be 0.1Å because the contribution from plans with d hkl < 0.1Å to σ el coh (E) can be neglected, as discussed in reference [9]. Nonetheless, the threshold value can always be modified by the users. During the loop over hkl, plans with identical interplanar distance d hkl and form factor |F (τ hkl )| are grouped together since they contribute equally to σ el coh (E) (Eq. (27)). The number of identical plans (hkl) is referred to multiplicity M hkl . These physical quantities are calculated by CINEL and stored in ASCII file, as illustrated in Figure 3.
If the incident neutron energy E is provided, the diffraction angle θ hkl corresponding to the Bragg edge E hkl is given by: The orientation angle Θ hkl and P hkl (Θ hkl ) can be obtained if the preferred orientation of crystal and the PDDF model are given in input. To improve the performance of Cubic and INELastic, a Just-in-Time (JIT) compiler named Numba [41] is used. Numba allows to reduce the computational time by directly adding simple Python syntaxes to the initial functions without rewriting them in low-level languages. In addition, Numba supports Compute Unified Device Architecture (CUDA) graphics processing unit (GPU) programming, which enables these two modules to benefit from the powerful computing capability of GPU.
GPU programming, introduced in the early 2000s, has been developed tremendously during the last two decades. In our work, the use of GPU enables to accelerate the phonon expansion calculation. To test the speedup with GPU, the TSL of hydrogen bound in water molecule (H in H 2 O) at room temperature were calculated with a phonon expansion order N phonon = 2000. Two kinds of GPU are used: NVIDIA ® QUADRO ® K620 with 2 gigabyte memory and 384 cores; NVIDIA ® QUADRO ® K6000 with 12 gigabyte memory and 2880 cores. The comparisons of the computational time are presented in Table 1. High performing GPU allows to significantly reduce the computational time down to 2 minutes. Approximations adopted in codes are summarized in Table 2. For strongly coherent crystalline materials such as pyrolytic graphite, the introduction of the coherent one-phonon correction [8] to the incoherent approximation enables to better reproduce the experimental inelastic scattering cross sections as shown in reference [42]. The calculations of one-phonon coherent inelastic scattering 2 For coherent one-phonon correction to incoherent approximation. function necessitates the dispersion relations of materials. This correction will be implemented in the future version of CINEL. The cubic approximation will also be eliminated since it only necessitates the partial PDOS which can be calculated by DFT codes.

SVT module
The SVT module is implemented with the Sampling the Velocity of the Target nucleus (SVT) algorithm [43,44], in which the neutron-nucleus scattering is approximated by the two-body kinetics as illustrated in Figure 4. The neutron mass and the mass of target nucleus are respectively m and M with M = Am. The speeds of neutron and target nucleus are respectively v and V . Angles between two speeds are represented by θ. The subscripts of v or θ indicate their belonging frame and the superscript represents the physical quantities after scattering. For example, v TR and v TR are respectively the speeds of neutron in the target-at-rest (TR) frame before and after scattering. For the sake of clarity, the subscripts of quantities in the laboratory (LAB) frame are omitted.
As presented in the data flow of CINEL (Fig. 2), the incident neutron energy E and the velocities of the target nuclei V are necessary for the calculation of the scattered neutron energy distributions. V are either sampled from trajectory files obtained by molecular dynamics (MD) calculations, or from a Maxwell-Boltzmann (M-B) distribution M T (V ) which is determined by the temperature or effective temperature of the nuclei. In the latter case, the cosine of the angle between the incident neutron and the target nucleus µ = cos(θ) before collision is sampled uniformly on [−1; 1].
The equations implemented in the SVT module to calculate the scattered neutron energy E and the cosine of the scattering angle µ = cos(θ ) are summarized in the following. The first step consists of calculating the relative velocity between the incident neutron and the target where v = 2E/m represents the incident neutron velocity. This step enables to transform the neutron-target system from the LAB frame to the TR frame, as shown in Figure 4, from the top plot to the bottom left plot. Next step is to transform the studied system from the TR frame to the center-of-mass (CM) frame in which the neutrontarget collision is assumed to be isotropic. The speed of the CM frame in the TR frame c TR is computed by: Then the speeds of neutron and target in the CM frame v CM and V CM can be obtained respectively by subtracting the speed of the CM frame: and Since the neutron-target collision is assumed to be isotropic in the CM frame, the cosine of the scattering angle µ CM = cos(θ CM ) is sampled uniformly on the interval [−1; 1]. The total momentum and the kinetic energy of the neutron-target system are conserved before and after the collision. It can be easily shown that the scalar velocities of the neutron and the target remain invariant in the CM frame, i.e., v CM = v CM and V CM = V CM , as illustrated in the bottom right plot of Figure 4.
The following step is to transform back the neutrontarget system from the CM frame to the TR frame. Based on geometric relations, speeds in the TR and CM frames lie on the same plan (P ). The scalar velocity of the postcollision neutron v TR and the cosine of the scattering angle µ TR = cos(θ TR ) in the TR frame can be respectively calculated by: and The final step consist of transforming the neutrontarget system back to the LAB frame to calculate v and µ = cos(θ ). The target speed V lying on the plan (P ) which does not coincide with the plan (P ) except for µ = ±1. Thus, an azimuthal angle φ CM must be sampled. φ CM is sampled uniformly on the interval ] − π; π] based on the assumption that the neutron scattering is isotropic in the CM frame. In summary, there are four variables sampled in the SVT algorithm [45]: V and µ are related to the motion of target nuclei, and µ CM and φ CM are used in the collision kinematics. µ r = cos(θ r ) is the cosine of the angle between the relative speed v − V and the speed of the target nucleus V : The speed of the scattered neutron in the TR frame v TR is decomposed into two components with one parallel to the relative speed v − V and the other perpendicular to it, i.e., v TR = v TR, + v TR, ⊥ . Equally, the speed of the target nucleus is decomposed into V = V + V ⊥ . Hence, the speed of the scattered neutron in the LAB frame v can be obtained by: Based on geometric relations, the velocity of the scattered neutron in the LAB frame v can be calculated by: where µ V = cos(θ V ) represents the cosine of the angle between v TR and V , which is obtained thanks to the decompositions of v TR and V in equation (69) and µ V can be computed thanks to the calculations of µ TR and µ r from equations (67) and (68). Since v TR can be obtained from equation (66), the velocity of scattered neutron in the LAB frame v in equation (70) can be determined. Finally, the cosine of the scattering angle µ = cos(θ ) can be calculated by:

Numerical validations of the cubic and INELastic modules
As presented in the data flow of CINEL in Figure 2, the Debye-Waller factor Λ(T ) and the coherent elastic scattering cross section σ el coh (E) computed from the Cubic module, and the symmetric TSL (S(α, β)) calculated from the INELastic module are merged into a file in ENDF-6 format. The ENDF files are later processed by using the THERMR module of NJOY to obtain the inelastic scattering cross sections. The numerical validations of the Cubic and INELastic modules of CINEL were performed by comparing the calculated scattering cross sections of solid crystalline materials with the ENDF/B-VIII.0 database and the NJOY+NCrystal library.

Automatic generation of α and β grids: UO 2
Uranium dioxide (UO 2 ) is the major component of the UOX fuel. In the ENDF/B-VIII.0 database, UO 2 as ideal fluorite structure with Fm3m symmetry is illustrated in Figure 5a. The atom positions in the unit cell are presented in Table 3.
The PDOS of 238 U in UO 2 and 16 O in UO 2 for the latest ENDF/B-VIII.0 database are generated by using the ab initio lattice dynamics methods [48]. The two peaks of the PDOS of 238 U in UO 2 in Figure 6 correspond to the acoustic modes of the uranium atoms. The peak at nearly 10 meV corresponds to the transverse acoustic (TA) mode while the peak at nearly 20 meV corresponds to the longitudinal acoustic (LA) mode. The peaks at nearly 30 meV, 55 meV and 70 meV are dominated by the optical modes of the oxygen atoms. The peak at roughly 30 meV stems from the transverse and longitudinal optical (TO1, LO1) modes. The peak at nearly 55 meV corresponds to the transverse optical mode (TO2) and the last peak at nearly 70 meV corresponds to the longitudinal optical mode (LO2).
CINEL is able to generate automatically the α and β grids for the calculations of the S(α, β) tables. Important steps are summarized in the following. The first step consists of determining the minimum and the maximum values of the α and β grids. According to equation (58), the minimum momentum transfer α min is obtained when µ = 1: The scattered neutron energy is regarded as having an energy gain δE compared to the incident energy, i.e., E = E + δE, then equation (73) can be rewritten as: The momentum transfer achieves its maximum value α max when µ = −1: Since the incident neutron energy has a threshold E thr , α max in equation (75) is obtained: In the case of the energy transfer β, based on equation (59), the minimum value β min and the maximum value β max are respectively:  and The second step is to fix the number of grids N α and N β and the distribution modes: logarithmic, linear or mixed.
In the case of UO 2 , the value of δE is set to be 2.8 meV. The energy threshold E thr = 5 eV, N α = 150 and N β = 200. Note that these parameters can be modified by the users.
Results are presented in Figure 7. The total scattering cross sections of UO 2 for a wide temperature range (from room temperature to 1200 K) obtained by using the Cubic and INELastic modules are in very good agreement with the ENDF/B-VIII.0 database.

Mixed elastic format: MgH 2 and MgD 2
Magnetism hydride (MgH 2 ) is investigated as a potential candidate for cold neutron reflectors [51]. Magnetism deuteride (MgD 2 ) shows improvements on neutron slowing down since the neutron capture cross section of deuterium is lower than that of hydrogen. MgH 2 and MgD 2 share the same crystal structure which is illustrated in Figure 5b. The PDOS of components in MgH 2 and MgD 2 generated by Campi and Bernasconi from University of Milano-Bicocca within the density functional perturbation theory [49,50], are presented in Figure 8.
To overcome the shortage of the ENDF-6 format in which the coherent and incoherent elastic components cannot be encompassed in the same file, the scientific community has proposed recently a mixed elastic format.   24 Mg in MgD2 and D in MgD2, generated by using the density functional perturbation theory [49,50].
This new option is supported in CINEL when generating the ENDF files. The mixed elastic cross sections are illustrated with MgD 2 in Figure 9. The comparisons of the CINEL calculated total cross sections of MgH 2 and MgD 2 with the NJOY+NCrystal library are presented in Figure 10. An excellent agreement with the NJOY+NCrystal library is obtained.
The total cross sections (elastic and inelastic) are in excellent agreement with the ENDF/B-VIII.0 database and the NJOY+NCrystal library, which validates the Cubic and the INELastic modules. The symmetric S(α, β) tables generated by using the INELastic module are in excellent agreement with the ENDF/B-VIII.0 database, as illustrated in Xu's thesis [52].

Anharmonicity of oxygen
The anharmonic motion of the oxygen atoms in UO 2 with fluorite symmetry (illustrated in Fig. 5a) is taken into account by using the third-cumulant coefficients c j 123 [25]. c U 123 are zero due to the symmetry of UO 2 . The impact of introducing the non-zero coefficient c O 123 on the coherent elastic cross sections is illustrated in Figure 11. The positions of the Bragg diffraction peaks remain unchanged. However, the impact on the coherent elastic cross section σ el coh (E) increases with the incident neutron energy E until σ el coh (E) reaches the asymptotic part (when E is larger than about 0.1 eV). The impact of c O 123 on the Monte-Carlo simulations of uranium dioxide is discussed in reference [53].

Texture: zirconium
Texture or preferred orientation presented in solid crystalline materials impacts the neutron diffraction compared to materials with random orientation. In this work, the texture is taken into account by the cylindrically symmetric Pole-Density Distribution Function (PDDF) P hkl (Θ hkl ). The PDDFs proposed by March-Dollase [21,22] and Altomare et al. [23] are respectively given by: and P Al hkl (Θ hkl ) = exp(P 1 cos(2Θ)), where P 1 is given by the users. The comparison of the PDDFs P MD hkl (Θ hkl ) and P Al hkl (Θ hkl ) is illustrated in Figure 12. These PDDFs are applied to the zirconium crystalline with preferred orientation PO = [0, 0, 1] at room temperature. The obtained coherent elastic scattering cross sections of zirconium are compared to the calculations with random orientation in Figure 13. The positions of Bragg diffraction peaks are invariant. The modifications of the peak intensities correspond to the introduced PDDFs presented in Figure 12. where The scattered neutron energy distributions obtained by using the SVT module of CINEL and Tripoli-4 ® are compared in Figure 14. Overall good agreement is obtained for all the angles, which validates the two-body kinematic equations implemented in the SVT module of CINEL.

Comparison with analytical results of the free gas model
The DDXS for free gas materials is obtained by converting the dynamic structure factor S( Q, ω) (Eq. (41)) to dimensionless S(α, β): The analytical results for free 1 H with incident neutron energies E = 1 and 5 eV are compared with the scattered neutron energy distributions obtained by using the SVT module in Figure 15. Excellent agreement is obtained for E = 5 eV, because the SVT algorithm satisfies the assumptions made in the analytical free gas model: the momentum and energy transfers are large, and the velocities of target nuclei follow a M-B distribution, as presented in Section 2.3. Discrepancies are observed for large scattering angles when E = 1 eV. Since in this case, the first assumption is not satisfied. Similar results are found for free 16 O, as presented in Figure 16. These results enable to validate the sampling of the velocities from a M-B distribution in the SVT module.
The comparison with Tripoli-4 ® calculations enables to validate the two-body kinematic equations implemented in the SVT module (as presented in Sect. 3.2). We remark in Section 2.3 that a simplified S(α, β) formula for free gas model is obtained thanks to the assumption of M-B distribution. Thus, the comparison with the analytical free gas models allows to validate the sampling of the velocities from a M-B distribution in the SVT module. Therefore, we think that these two comparisons are necessary.

Comparison with velocities sampled from MD trajectory files
As presented in the data flow of CINEL (Fig. 2), the SVT module is able to sample the velocities from the MD calculated trajectory files. In this work, the trajectory files of liquid water generated by using the MD code GROMACS [55] are used to validate the velocity sampling.
The generation of the trajectory files for liquid water by using GROMACS is described in Scotta's thesis [56]. The physical model used in the MD simulation consists of in total 512 water molecules within a cubic box with side length 2.48 nm. We calculated the number volume density of water at 294 K by ρ V = N A M H2O ρ H2O = 33.5 nm −3 where N A is the Avogadro constant, M H2O and ρ H2O are respectively the molar mass and mass density of water. Then the size of the box is determined by (512/33.5) 1/3 = 2.48 nm. TIP4P/2005f [57], a flexible water potential model taking into account the intermolecular and intramolecular interactions, is adopted to quantify the forces between atoms. Thanks to this potential, the positions and velocities of the hydrogen and oxygen atoms in the water molecules are obtained for a time step of 0.6 fs and they are stored in the trajectory files with a time duration of 100 ps. Each water molecule is composed of two hydrogen atoms and one oxygen atom. For each time step, the partial velocities along the x, y and z directions are respectively V j,x , V j,y and V j,z , with j representing the hydrogen atom H or oxygen atom O. The velocities V j is thus computed via the partial components as follows: The distributions of velocities of H and O are presented in the left hand plot of Figure 17. The velocity µ j representing the cosine of the angle between the z components V j,z and the velocities V j for oxygen and hydrogen atoms, is calculated by: µ j is shown in the right hand plot of Figure 17, which is nearly an uniform distribution in [−1;1]. The velocities V j and the corresponding cosine of angle µ j are sampled from the MD trajectory files by using the SVT module. The scattered neutron energy distributions for H and O are compared with the results obtained with the velocities sampled from a M-B distribution at 294 K, as presented in Figure 18. Excellent agreement is obtained, which validates the sampling of the velocities from the MD trajectory files implemented in the SVT module. Compared to a M-B distribution, velocities in the MD trajectories are strongly correlated between nearby time steps. These correlations can be illustrated with the velocity auto-correlation function (VACF), as shown in reference [56]. While the correlations between the velocities in the trajectories cannot be directly taken into account by the SVT algorithm. We have tried to incorporate these correlations into the SVT algorithm but our attempts did not succeed.

Discussions on rejection option to improve the neutron up-scattering
The free gas model or SVT treatment for neutron-nucleus scattering at thermodynamic temperature T are not able to take into account the vibration dynamics or atom binding effects of materials, but the effective temperature T eff must be involved, as discussed in the literature [60][61][62][63]. In this work, special attentions are given to the neutron upscattering: when neutron energy lies below a few eV, the neutron-nucleus scattering may lead to an energy gain for the incident neutron. Based on the short collision approximation (SCT), after introducing the effective temperature T eff to the DDXS of free gas materials (Eq. (83)), the upscattering part of the DDXS must be further corrected to fulfill the principle of detailed balance [62]. The obtained DDXS involving T and T eff is given by: where T eff is calculated from the PDOS [13,61]: Compared to the DDXS obtained by using the free gas model at T eff , the DDXS based on the SCT can be regarded as multiplying a factor for the up-scattering part (when E < E ): It is confirmed in Section 5.2 that the SVT algorithm is equivalent to the analytical free gas model for large incident neutron energies. Based on the relation in equation (88), the SCT can be introduced in the SVT module by implementing an additional rejection step in case of neutron up-scattering: when E < E , a variable R is sampled uniformly in [0;1]; by comparing R and exp is larger, the scattered neutron will be rejected.
To validate the rejection step implemented in the SVT module, neutron scattering with 1 H in 1 H 2 O and 16 O in U 16 O 2 are investigated in this work. The first step consists of calculating their effective temperatures from the PDOS of the CAB model [58] and the PDOS reported in reference [59], respectively. The second step is to employ these effective temperatures in the SVT module and use the rejection option for the up-scattering part. The final step  [58]. T eff of 16 O in UO 2 is obtained via the PDOS of [59].
Target nucleus E (eV) T (K) T eff (K) Scattering treatments 1  Tripoli-4 ® with S(α, β) is to perform Monte-Carlo simulations with the Tripoli-4 ® code by using the S(α, β) tables. This step is able to provide accurate scattered neutron energy distributions for comparison. Tested models are summarized in Table 4. The comparison of results obtained by using the SVT module and the Tripoli-4 ® code for 1 H in 1 H 2 O is illustrated in the top plot of Figure 19. At small scattering angle (θ = 10 • ), the SVT at T eff with/without the rejection both fail to reproduce the quasi-elastic peak and the associated structures. However, the SVT treatment with rejection (green curve) reproduces better the scattered neutron distribution related to the up-scattering part. At large θ , the neutron up-scattering can be neglected. Therefore, these two SVT treatments produce the same distributions. Same comparisons are performed in the case of 16 O in U 16 O 2 . The obtained results are presented in the bottom plot of Figure 19. At small θ , both SVT treatments are incapable of reproducing the structures related to the atomic motions of the oxygen atoms bound in UO 2 . This result is in consistency with the case of 1 H in 1 H 2 O. The effective temperature of 16 O in U 16 O 2 (T eff = 381 K) is close to T = 298 K. Therefore, the up-scattering part produced by using the SVT algorithm with rejection does not show significant improvements compared to the SVT treatment at T eff .
These results show interest to add a new rejection option in the Tripoli-4 ® code that allows to correct the up-scattering by taking into account T and T eff in the simulation. This new option will be of great interest to test the impact of different scattering models on integral benchmarks.

Conclusions and prospectives
In this work we present a new code CINEL dedicated to generate the thermal scattering files in ENDF-6 format for solid crystalline, free gas materials and liquid water. The low-energy neutron scattering theory is reviewed by presenting the involved approximations and the equations to calculate the elastic (coherent and incoherent) scattering cross section and S(α, β). The implementation of the three modules of CINEL (Cubic, INELastic and SVT) is presented. The Cubic and INELastic modules are validated by comparing the CINEL results with the ENDF/B-VIII.0 database and the NJOY+NCrystal library. Examples are given to illustrate the automatic generation of α and β grids, mixed elastic format, anharmonicity and texture in CINEL. The SVT module is validated by comparing the scattered neutron distributions with analytical and Monte-Carlo simulations for free target nuclei such as hydrogen and oxygen. Though CINEL presents improvements compared to the LEAPR module of the NJOY code, the incoherent approximation adopted in the CINEL code is reported to be inadequate for strongly coherent crystalline materials such as pyrolytic graphite. Therefore, the coherent Fig. 19. Comparison of scattered neutron energy distributions of 1 H in 1 H2O and 16 O in U 16 O2, obtained by using the SVT module of CINEL at T eff with rejection, and Tripoli-4 ® with S(α, β) for incident neutron energy E = 5 eV and 6.67 eV, respectively. one-phonon correction will be implemented in the future version of CINEL to improve the incoherent approximation. The cubic approximation used in the current version of CINEL will also be eliminated. CINEL will also be extended to calculate the thermal scattering files for liquid materials other than water.