Generation of thermal neutron scattering libraries for liquid para-hydrogen and ortho-deuterium using ring-polymer molecular dynamics

. In this paper we present results of combining ring-polymer molecular dynamics with the LEAPR and THEMR modules of NJOY to generate thermal neutron scattering libraries for liquid para-hydrogen and ortho-deuterium. We present the methodology and show that it produces results that are in good agreement with data from both recent available measurements and previous theoretical studies. We also present some simple benchmark Monte-Carlo simulations compared with other available libraries.


Introduction
The design of moderator and reflector systems at neutron scattering facilities are typically carried out using radiation transport software, which uses thermal neutron scattering libraries that are derived from the atomic and molecular properties of the materials. The generation of such libraries is usually done with the LEAPR/THERMR modules of the NJOY code system [1]. Of particular importance for the production of cold neutrons are libraries for para-hydrogen and ortho-deuterium. The current moderator system to be designed at the European Spallation Source (ESS) [2], for example, is based on a high-brightness liquid hydrogen moderator, while a highintensity liquid deuterium moderator is being designed [3].
The currently packaged models with NJOY for liquid para-hydrogen and ortho-deuterium are based on the work by Keinert and Sax [4]. Their models use analytical methods to treat the translational degrees of freedom, combined with the model from Young and Koppel [5] for quantum rotation, and the Vineyard approximation [6] for inter-molecular coherence. This model is the basis for the ENDF-VIII.0 evaluations [7] and was improved and updated by Granada et al. [8,9] and more recently was further refined through the implementation of the Sköld approximation [10] for the inter-molecular coherence. These recent libraries are included in the JEFF 3.3 [11] evaluations, and described in more detail in Ref [12]. The software for this implementation is available online in a Github repository under the NJOY-H2D2 branch [13]. In the above methodologies, the data used as input to the models are derived from theoretical and/or experimental information.
In a series of additional works, Guarini et al. [14,15] used ring-polymer molecular dynamics (RPMD) [16] to compute the neutron scattering cross-section for liquid para-hydrogen and ortho-deuterium, however neutron * e-mail: douglas.dijulio@ess.eu scattering libraries suitable for radiation transport calculations were not generated. Their approach offers the flexibility to derive temperature dependent input to models used for library generation and also to explore materials where limited previous experimental or theoretical treatment exits. In this work, we present results of combining the above mentioned methodologies, thus providing a means for generating thermal neutron scattering libraries based on quantum molecular dynamics techniques. The work presented here represents a culmination and complete update of the earlier studies carried out on these materials [17][18][19].

Methodology
For a homonuclear diatomic molecule, where the intramolecular degrees of freedom are considered to be completely decoupled from the translational degrees of freedom, the starting point is the double differential neutron scattering cross-section, which can be written as [20,21], where S d (Q, ω) and S s (Q, ω) are the center-of-mass distinct and self dynamic structure factors, respectively, k 1 and k 0 are the scattered and initial neutron wave vectors, Q and ω are the wave vector and energy transferred in the scattering process, and J and υ are the rotational and vibrational quantum numbers. The factor u(Q) is dependent on the coherent cross-sections while the function F(Q, J 0 , J 1 , υ 1 ) depends on the nuclear spin statistics and both the incoherent and coherent cross-sections. Both of these factors can be computed using the formalisms described in [5,20,22], which reduces the calculation of the double differential cross-section to determination of the self-and distinct dynamical structure factors. The self part can be directly linked to the output of the RPMD simulations through the intermediate scattering function, which is the time Fourier transform of S s (Q, ω), i.e. F s (Q, t) = e −Q 2 γ(t) , where we have used the Gaussian approximation [6,23] to write the intermediate scattering function in terms of the width function γ(t), which can further be expressed as [23] The frequency spectrum f (ω) provides the link to the simulations, which is given by where M is the mass of the particles, i.e. H 2 or D 2 in this case, k b being the Boltzmann constant, and T the temperature of the system [14,[24][25][26]. The function ⟨v(0) · v(t)⟩ (k) is the Kubo transformed velocity-autocorrelation function (KVACF) [27] and is the output of a RPMD simulation.
For the distinct part, we make use of the the Sköld approximation [10], which leads to where S (Q) is the center-of-mass static structure factor.

Dynamics
In NJOY-H2D2, the Gaussian approximation is calculated using the phonon expansion method [1] and is limited to solid-like frequency distributions, i.e. distributions that tend to zero as ω goes to zero. This is not true for liquids, which exhibit non-zero frequencies distributions at ω =0. Instead, in NJOY-H2D2 the frequency spectrum, related to f (ω), is described by the combination of a solid-like component and a diffusive component, where the Egelstaff-Schofield model [28] is used to describe the latter part and is given by [29] which is the diffusive weight within the system, and M diff is the effective diffusion mass of the particles. The usage of this model is convenient as it has an associated frequency spectrum, given by [29] In order to facilitate input into NJOY-H2D2, the frequency spectrum calculated from eq. (3), must be separated into diffusive and solid-like components. Following the prescription outlined in Ref. [29], this can be carried out through calculation of an appropriate f diff (ω) and the subtraction of it from f (ω), which yields a solid-like frequency spectrum. The two relevant parameters for this procedure are the diffusion constant, D, and the weight of the diffusive motion within the system, w t . The first parameter can be calculated directly from the output of RPMD simulations, i.e.
The parameter w t is however not directly available from the output of the simulations and instead we have made use of the two-phase model described in [30,31]. The diffusive weight is obtained from the the fluidicity factor, f , as defined in [30,31]. For input into NJOY-H2D2, the total weight of the solid-like and diffusive components should sum to 0.5 [1], which results in w t = f /2. For these calculations, we used the self-diffusion constants taken from the RPMD simulations.
In addition, in our calculations we used υ 1 =0 in Eq. (1), and we have followed the method given in [9] to include the intra-molecular vibrations, where we used a discrete oscillator with an energy of 0.546 eV, for H 2 , and 0.386 eV, for D 2 , with a weight of 1/6 to describe the vibration. The self dynamic structure factor is calculated as a convolution of the vibrational component, the diffusive component and the solid-like component, as described in [1,9].

Structure
The center-of-mass static structure factor is the remaining input to be provided to NJOY-H2D2. It can be computed directly from the output of RPMD simulations using [32,33] where n is the number density, h(r) = g(r) − 1 and g(r) is the center-of-mass radial distribution function, computed from the trajectories of the RPMD simulations. In order to avoid truncation effects during the integration of Eq. (7), we have extrapolated g(r) to higher values of r [33]. For this purpose, we used the function h(r) ≈ (A 0 /(r − r 0 ))e −(r−r 0 )/r 1 sin((r − r 0 )/r 2 ), with the free parameters r 0,1,2 and A 0 , fit to the simulation results starting around the position of the third zero of h(r) [34].

RPMD Simulations
For the RPMD simulations we used i-PI v2.0 [35] with the Silvera-Goldman potential [36] that is packaged with the software. When calculating the KVACF we set up a simulation box consisting of N = 512 molecules and equilibrated the system in the NVT ensemble for 50 ps and set the Trotter number to P = 64. The KVACF was calculated in the NVE ensemble and up to a maximum time lag of 2.5 ps and calculated from the trajectory files using the software included in the i-Pi distribution. The spectrum, f (ω), was computed by zero-padding and multiplying the KVACF by a Welch window using the pwtools toolkit [37]. For calculations of S (Q), we have followed the above equilibrium procedure but set the number of molecules to N = 4096. Calculation of the g(r) function from the trajectory files of the simulations was performed using the TRAVIS software package [38,39] and calculated from the bead positions of the simulations [40,41].
We have carried out the RPMD simulations for the temperatures of 15.7 K and 20 K for para-hydrogen and 19 K and 20 K for ortho-deuterium. The first two temperatures for each material were used as benchmarks against previously measured experimental data for the total scattering cross-sections, while the second temperatures for each material were used for comparisons of Monte-Carlo simulations and represent the temperatures of the moderator system at the ESS.  3 Results Fig. 1 shows the calculated KVACF functions from our RPMD simulations for the two different temperatures for each material. From these functions, Eq. (3) was used to calculate the frequency spectra and the results for H 2 at 15.7 K are shown in Fig. 2. The frequency spectra for diffusion, calculated using Eq. (6), and the resulting solidlike spectra is also shown. Table 1 shows the diffusion weights, w t , used for the NJOY-H2D2 calculations and the diffusion constants calculated from the RPMD simulations. The values D RPMD are those which are derived directly from the RPMD simulations, where it is known that they underestimate the experimental values [43]. Thus we used experimental values where possible and the self-diffusion constant for 19.0 K D 2 was scaled from the experimental value at 20.0 K based on our RPMD results. Figure 3 shows the calculated static structure factors compared to experimental data from [32,33,46]. For para-hydrogen, it can be seen that there is good agreement to the data of [33] at low Q and there is better agreement with the data from [46] at high Q. These differences are further discussed in [33]. On the other hand, it can be seen that there is good agreement for the ortho-deuterium results.

Total scattering cross-sections
Using the calculated data from the previous section, we prepared neutron scattering libraries with NJOY-H2D2. Each library was processed using the same number of outgoing angles, in continuous energy format, and to the indicated temperature. The calculated scattering crosssections for para-hydrogen at 15.7 K are shown in Figure 4 for the important region where the experimental data shows deviations. It can be seen that there is good agreement with both the calculations from Guarini et al. and the experimental data from Grammer et al. The cross-sections for 20 K are shown in Fig. 5 and compared against the two other libraries. The neutron absorption cross-section is also plotted as calculated from NCrystal [50]. Below 10 meV, the JEFF and RPMD data starts to deviate from the ENDF data and below 1 meV all three calculations deviate from each other. At this energy neutron absorption becomes dominant.
The total scattering cross-sections for ortho-deuterium . Calculated static structure factors compared to data digitized from [33,46] and [32].  . Calculated scattering cross-section for para-hydrogen compared to digitized data from previous calculations [14], experiments [47,48], and libraries from JEFF 3.3 [11]. Experimental data is taken from EXFOR [49] and corrected with the calculated absorption from NCrystal [50]. at 19 K are shown in Fig. 6, where it is seen that the agreement with the JEFF data and experimental data is in general good overall, with larger deviations from the ENDF data.

Monte-Carlo Simulations
To compare the libraries even further, we prepared simple test problems in MCNP6.2 [53], where a mono-energetic 1 MeV neutron source was placed in a sphere, comprising either para-hydrogen or ortho-deuterium, and we tallied the neutron flux leaving the sphere surface. We use a 10 cm radius sphere for the para-hydrogen case and a 30 cm radius sphere for the ortho-deuterium case. Figure 7 shows the results for ortho-deuterium, where   Calculated scattering cross-section for orthodeuterium compared to libraries from JEFF 3.3 [11], ENDF-VIII.0 [7], and data points [52] at 19 K. The data from [48,51] is taken from EXFOR and corrected for the absorption calculated by NCrystal.
good agreement between the current calculations and the JEFF libraries can be seen. At lower energies, the calculations based on the ENDF libraries deviate. For the para-hydrogen case, we a saw smaller maximum deviation of 10% across the same range compared to the ENDF libraries, while maximum deviations from the JEFF library were on the few percent level.

Conclusions
We have developed a methodology using ring-polymer molecular dynamics, combined with recent developments in neutron scattering library software, to produce thermal neutron scattering libraries for para-hydrogen and orthodeuterium. We compared the calculated total scattering cross-sections with experimental data and also carried out Monte-Carlo simulations in order to compare against other available libraries. Generally we find that there is good agreement with calculations based on the JEFF libraries.