Hydrogen Deuteride for Cold Neutron Production: A Model for the Double Differential Cross Section

Abstract

The present work deals with the modeling of the response to neutrons of heteronuclear diatomic liquids, with special interest in the case of hydrogen deuteride (HD), as a possible candidate for the moderation process required in the production of cold neutrons. Preliminary evaluations of the model giving the neutron double differential cross section of a heteronuclear vibrating rotor were performed in the recent past by using, as a first approximation, the ideal gas law for the center-of-mass translational dynamics. Here, the state-of-the-art methodology (based on the use of quantum simulations of the velocity autocorrelation function) for predicting the neutron response of moderately quantum fluids (like molecular hydrogen and deuterium at low temperatures) is applied to the heteronuclear form of this molecular liquid. The unavailability of the double differential cross section experimental data on liquid HD still compels us to test the calculations only at an integral level, i.e., against the only available measurements of the total neutron cross section of HD. Despite the well-tested and parameter-free computational approach, which includes proper consideration of the quantum effects, the present findings on HD indicate the evident need for more accurate measurements of its total cross section in extended ranges of incident energy, as well as of an experimental determination of the double differential cross section of this mild quantum liquid. For further applicative purposes, a very useful by-product of this study is the determination of the self diffusion coefficient D of the HD in the liquid phase.

1. Introduction

Liquid hydrogen and its isotopic forms are among the most renowned cryogenic fluids and, in the specific application to neutron techniques, are the most important low-temperature moderating materials used to realize cold (0.1–10 meV) [1] neutron sources. In recent years, quite an effort has been devoted to deeply refine the ability to predict the neutron scattering properties of these diatomic fluids and build up reliable databases collecting their double differential cross section (DDCS) in as wide kinematic ranges as possible. As far as hydrogen (${\mathrm{H}}_2$) and deuterium (${\mathrm{D}}_2$) are concerned, significant progress was possible some years ago when crucial experimental work was dedicated to accurate measurements of the total scattering cross section ${\sigma }_{\mathrm{s}}$ of para-${\mathrm{H}}_2$ [2] and, at the same time, quantum calculation methods have been employed, in place of classical computations of the translational dynamics, in the DDCS algorithms for these low-mass molecules undergoing nonnegligible quantum delocalization effects [3,4]. More recently, the latter approach, which gives the unique opportunity to get rid of adjustable parameters in DDCS evaluations, has been successfully adopted in moderator design and nuclear data processing codes [5,6] implemented at the European Spallation Source (ESS, Sweden). Recent experiments aimed at a better assessment of the coherent scattering properties of ${\mathrm{D}}_2$, and providing absolute scale determinations of its DDCS by inelastic neutron scattering [7], will also help further refinement of the scattering libraries for ${\mathrm{D}}_2$. Therefore, concerning the homonuclear representatives of the hydrogen family, we can be rather satisfied with the present computational and predicting capabilities, now permitting well-grounded and safe applications to forthcoming and existing neutron sources.

Nonetheless, another member among the hydrogen liquids, namely hydrogen deuteride (HD), has not received the same attention as its homonuclear partners in recent years, either from a scientific point of view or in applications, although some of its midway properties between ${\mathrm{H}}_2$ and ${\mathrm{D}}_2$ might also be conveniently exploited in neutron moderation processes, as explored in previous moderator experimental studies [8,9]. For instance, at the same (manageable) low temperatures of liquid ${\mathrm{H}}_2$ and ${\mathrm{D}}_2$ (∼20 K), HD provides stronger scattering than ${\mathrm{D}}_2$ and lower absorption than ${\mathrm{H}}_2$, while its molecular mass does not critically decrease the moderation efficiency, as shown by the successful use of an even heavier liquid, like ${\mathrm{D}}_2$, for cold neutron production. Also, mixtures of ${\mathrm{H}}_2$ and ${\mathrm{D}}_2$ (partly combining themselves in the form of HD molecules) could be envisaged. Indeed, it has been suggested that adding a small percentage of ${\mathrm{H}}_2$ to liquid ${\mathrm{D}}_2$ might improve the moderation efficiency in the case of small moderator volumes [10,11].

Nevertheless, when exploring the literature about experimental, theoretical and simulation work on liquid HD and, in particular, its neutron scattering properties, one is compelled to face an unexpected lack of published results, with the only exception being Ref. [9] which, however, requires a rather detailed inquiry about various unreported instrumental parameters needed for an accurate modeling of the experimental conditions. It is also quite surprising that the neutron response of HD, which belongs to an important class of systems anyway, has, for decades, escaped a detailed treatment like those elegantly devised, in the mid 1960s, on homonuclear diatomic molecules by Young and Koppel (YK) [12] and by Sears [13], or on other linear molecules by Lurie [14]. Indeed, only much later, the neutron cross section of HD has been briefly taken into consideration in the analysis of solid state data [15], and further years passed before a formal description of the heteronuclear diatomic case was tackled in studying the behavior of HD in the cages of clathrate hydrates [16,17]. However, like the papers by Sears [13] and by Lurie [14], these last works also deal with the scattering of cold and thermal neutrons from a low-temperature sample, therefore vibrations are not excited and the developed formalism is rightly limited to the specific case under consideration, where rotations alone are excited and only zero-point vibrational effects must be retained.

In order to extend the applicability of the model to the higher incident energies involved in moderation at neutron sources as well, we therefore found it important to provide a more general treatment of heteronuclear diatomic fluids, including harmonic vibrations, as formalized for the homonuclear case in Refs. [12,18], and as we performed in Refs. [3,4,19,20], and overcoming the use of vibrational Debye–Waller factors [15]. Here, we summarize the formalism for the DDCS of a heteronuclear (harmonically) vibrating rotor, and use it in connection with quantum calculations of the translational dynamics, following the same methodology described in Refs. [3,4]. The agreement of the present computational results for HD with the only available total cross section experimental data, collected by Seiffert in the 1970s [21,22], is only found to be very good by including systematic errors and by considering impurities. Future efforts regarding liquid HD should therefore be devoted to new neutron scattering experiments, taking advantage of the much higher performance of the instrumentation available nowadays.

2. Materials and Methods

Before summarizing the basic formalism, it is useful to recall under which hypotheses the adopted modeling of the DDCS holds [12,18,19]: (i) molecules are considered to be free vibrating rotors, i.e., the translational center of mass (CM) dynamics is assumed to be completely decoupled from the intramolecular motions; (ii) rotations are treated as independent of vibrations; and (iii) anharmonicity effects are neglected. The first hypothesis corresponds to consider an isotropic interaction among the molecules, with negligible orientational correlations. Assumptions (i) to (iii) are known to be very well satisfied for hydrogen, even in the (low pressure) solid phase [23]. However, since the centers of mass and charge do not coincide in the heteronuclear case, the full validity of hypothesis (i) for HD also can only be verified a posteriori. Moreover, for the case of low temperature liquids, an additional assumption is that all the molecules in the system initially (that is, before interaction with neutrons) lie in their ground vibrational state. Finally, only the case of unpolarized neutrons is considered here.

The above general assumptions for the treatment of the neutron cross sections of diatomic low temperature fluids are further accompanied by the important simplifications introduced by the heteronuclear nature of the molecule, which allows us to consider its nuclei as distinguishable particles. The absence of symmetry requirements for the total molecular wave function implies that there is no coupling between the total molecular spin and the rotational state of the molecule. Therefore, one can refer to the so-called “uncorrelated spin” case [19], where quantum-statistical averages involving nuclear and neutron spin variables can be carried out separately from those related to position variables.

In the mentioned hypotheses and conditions, the basic starting formulas leading to the double differential (per unit solid angle and exchanged frequency interval) cross section are those gathered in Equation (7) of Ref. [19] which, by omitting the superscript uncorr, we rewrite as:

$$\begin{array}{c}\hfill \frac{d^2\sigma }{d\Omega d\omega }=\frac{k_1}{k_0}{S}_{\mathrm{n}}(Q,\omega )=\frac{k_1}{k_0}\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }dt{e}^{-i\omega t}{F}_{\mathrm{n}}(Q,t)=\\ \hfill =\frac{k_1}{k_0}\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }dt{e}^{-i\omega t}[u\left(Q\right)F_{\mathrm{d}}(Q,t)+v(Q,t)F_{\mathrm{s}}(Q,t)],\end{array}$$

(1)

where $k_0$ and $k_1$ are the incident and scattered neutron wavevectors, and $S_{\mathrm{n}}(Q,\omega )$ is the total dynamic structure factor provided by neutron scattering, i.e., the time Fourier transform of $F_{\mathrm{n}}(Q,t)$, which is the neutron weighted combination of the distinct and self intermediate scattering functions $F_{\mathrm{d}}(Q,t)$ and $F_{\mathrm{s}}(Q,t)$, respectively. For molecular fluids, the latter are to be identified with the CM functions defined by

$$\begin{array}{c}\hfill F_{\mathrm{d}}(Q,t)=\frac{1}{N}\sum _{i=1}^N\sum _{j=1(j\ne i)}^N〈e^{-i\mathbf{Q}·{\mathbf{R}}_i\left(0\right)}{e}^{i\mathbf{Q}·{\mathbf{R}}_j\left(t\right)}〉\end{array}$$

(2)

$$\begin{array}{c}\hfill F_{\mathrm{s}}(Q,t)=\frac{1}{N}\sum _{i=1}^N〈e^{-i\mathbf{Q}·{\mathbf{R}}_i\left(0\right)}{e}^{i\mathbf{Q}·{\mathbf{R}}_i\left(t\right)}〉\end{array}$$

(3)

which provide the total CM intermediate scattering function $F(Q,t)=F_{\mathrm{d}}(Q,t)+F_{\mathrm{s}}(Q,t)$. In the above equations, N is the total number of molecules, ${\mathbf{R}}_i\left(0\right)$ is the CM position of the ith molecule in the (arbitrarily chosen) time origin, and ${\mathbf{R}}_j\left(t\right)$ is the CM position of a different molecule at a subsequent time t. The angle brackets denote the quantum canonical ensemble average. The isotropy of the fluid actually makes these functions depend only on the modulus Q of the exchanged wavevector $\mathbf{Q}={\mathbf{k}}_0-{\mathbf{k}}_1$. Furthermore, in the last member of Equation (1) the functions $u\left(Q\right)$ and $v(Q,t)$ play the role of inter- and intramolecular form factors weighting, respectively, the distinct and self CM dynamics. In analogy with the monatomic case, purely coherent scattering characterizes $u\left(Q\right)$, that is, only the so-called coherent part of the neutron cross section of the scattering unit probes the interparticle translational dynamics. As a consequence, $u\left(Q\right)$ exclusively contains the coherent scattering lengths of the various nuclei present in the molecule, and is independent of time. Differently, the intramolecular form factor is a function of time, and generally depends on both the coherent and incoherent nuclear scattering lengths. Explicit expressions for the form factors of a heteronuclear diatomic molecule at a low temperature are derived in Appendix A for direct implementation in computing code. There, the various parameters relevant to HD are also given.

The neutron DDCS is then obtained by inserting Equations (A14) and (A17) into Equation (1). Before discussing reasonable modelings of the CM translational dynamics, we wish to point out a few basic facts that are significant for the comparison of calculations with the (possibly available) neutron experimental spectra.

The first point regards what is actually accessed by experiments. It is well known that conventional neutron spectroscopy provides the Fourier transforms of space and time correlation functions. In Equations (2) and (3), we introduced the time autocorrelation of the space Fourier transform of the microscopic density, i.e., the total intermediate scattering function $F(Q,t)$, separated into its distinct and self parts. The latter were used in particular to describe the neutron version of the molecular scattering function as $F_{\mathrm{n}}(Q,t)=u\left(Q\right)F_{\mathrm{d}}(Q,t)+v(Q,t)F_{\mathrm{s}}(Q,t)$. Such a separation of the total $F_{\mathrm{n}}(Q,t)$, inherited from the formalism used to describe neutron scattering from monatomic fluids, i.e., $b_{\mathrm{coh}}^2{F}_{\mathrm{d}}(Q,t)+(b_{\mathrm{coh}}^2+b_{\mathrm{inc}}^2)F_{\mathrm{s}}(Q,t)$, has the merit of highlighting the relationship between the coherent scattering and the distinct dynamics. This also holds true for molecular liquids, since we saw that $u\left(Q\right)$ only contains the coherent scattering lengths of the constituent atoms. However, despite its conceptual significance, the mentioned separation has no feedback from reality, since neutrons provide a different combination and can only be the probe of “true” correlation functions, such as $F(Q,t)$ and $F_{\mathrm{s}}(Q,t)$, differently from $F_{\mathrm{d}}(Q,t)$. This means that the self and distinct contributions to the total dynamics cannot be disentangled in a neutron measurement on a totally coherent sample. Conversely, it is incoherent scattering that provides an (as remarkable as exclusive) pathway to the self dynamics. In other words, the output of a neutron experiment on a monatomic sample with nonzero coherent and incoherent scattering lengths is actually $b_{\mathrm{coh}}^{2}F(Q,t)+b_{\mathrm{inc}}^2{F}_{\mathrm{s}}(Q,t)$, the molecular version of which is

$$F_{\mathrm{n}}(Q,t)=u\left(Q\right)F(Q,t)+[v(Q,t)-u\left(Q\right)]F_{\mathrm{s}}(Q,t).$$

(4)

The second thing worth recalling concerns the spectral properties and the general features of the resulting DDCS. By completing the switching to Fourier $(Q,\omega )$ space, Equation (4) becomes

$$\begin{array}{c}\hfill S_{\mathrm{n}}(Q,\omega )=\frac{1}{2\pi }\{u\left(Q\right){\int }_{-\infty }^{+\infty }dt{e}^{-i\omega t}F(Q,t)+\\ \hfill +{\int }_{-\infty }^{+\infty }dt{e}^{-i\omega t}[v(Q,t)-u\left(Q\right)]F_{\mathrm{s}}(Q,t)\}=\\ \hfill =u\left(Q\right)[S(Q,\omega )-S_{\mathrm{s}}(Q,\omega )]+\left[\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }dt{e}^{-i\omega t}v(Q,t)F_{\mathrm{s}}(Q,t)\right].\end{array}$$

(5)

The first term is, therefore, related to distinct CM dynamic structure factor we are interested in when studying the translational collective dynamics of the system. Concerning the self properties, it is seen instead that the dependence of $v(Q,t)$ on time prevents one from expressing the second term as directly proportional to the self dynamic structure factor $S_{\mathrm{s}}(Q,\omega )$. However, Equation (A14) shows that time enters $v(Q,t)$ only in exponential form; therefore, it is possible to write

$$\begin{array}{c}\hfill {\int }_{-\infty }^{+\infty }dt{e}^{-i\omega t}v(Q,t)F_{\mathrm{s}}(Q,t)=\\ \hfill =\sum _{J_0{J}_1{v}_0{v}_1}{\int }_{-\infty }^{+\infty }dt{e}^{-i(\omega -{\omega }_{J_0{J}_1}-{\omega }_{v_0{v}_1})t}\mathcal{F}(Q;J_0{J}_1{v}_0{v}_1)F_{\mathrm{s}}(Q,t)\end{array}$$

(6)

where the explicit expression for $\mathcal{F}(Q;J_0{J}_1{v}_0{v}_1)$ in the low temperature case is given in Equation (A13) of Appendix A. In Equation (6), the rotational (${\omega }_{J_0{J}_1}$) and vibrational (${\omega }_{v_0{v}_1}$) transition frequencies, involving the initial (subscript 0) and final (subscript 1) quantum numbers, are seen to determine the time dependence of $v(Q,t)$. In the special case of a low temperature sample, all molecules are assumed to lie initially in the vibrational ground state ($v_0=0$); therefore, it is possible to write ${\omega }_{v_0{v}_1}=v_1{\omega }_{\mathrm{v}}$, with ${\omega }_{\mathrm{v}}$ indicating the frequency of the harmonic oscillator of reduced mass $\mu =\frac{m_1{m}_2}{M_{12}}$, with $m_1$ and $m_2$ the mass of the two (different) atoms in the molecule and $M_{12}=m_1+m_2$.

Consequently, the DDCS finally reads

$$\begin{array}{c}\hfill \frac{d^2\sigma }{d\Omega d\omega }=\frac{k_1}{k_0}{S}_{\mathrm{n}}(Q,\omega )=\frac{k_1}{k_0}\{u\left(Q\right)[S(Q,\omega )-S_{\mathrm{s}}(Q,\omega )]+\\ \hfill +\sum _{J_0{J}_1{v}_1}\mathcal{F}(Q;J_0{J}_{1}0v_1)S_{\mathrm{s}}(Q,\omega -{\omega }_{J_0{J}_1}-v_1{\omega }_{\mathrm{v}})\}.\end{array}$$

(7)

Equation (7) shows that the single-molecule contribution to the DDCS corresponds to a comb of lines centered at the frequencies of the possible rotovibrational transitions. These spectral components are therefore either central or shifted replicas of the lineshape describing the CM $S_{\mathrm{s}}(Q,\omega )$, with amplitudes ruled by the involved quantum numbers, the initial state probabilities and the nuclear scattering lengths.

Calculations of the DDCS, of course, require a modeling of both $S(Q,\omega )$ and $S_{\mathrm{s}}(Q,\omega )$. It is worth recalling that a comparison with the experiment is only possible if the model lineshapes obey the detailed balance principle. Therefore, if classical (i.e., symmetric) models are used for the dynamic structure factors, these must be duly asymmetrised via multiplication by the factor $\frac{\hslash \omega }{k_{\mathrm{B}}T}[n\left(\omega \right)+1]=\frac{\hslash \omega }{k_{\mathrm{B}}T}{[1-exp(-\frac{\hslash \omega }{k_{\mathrm{B}}T})]}^{-1}$ prior to their inclusion in Equation (7), where $n\left(\omega \right)$ is the Bose factor [24], with $k_{\mathrm{B}}$ and ℏ indicating the Boltzmann and Planck constants, respectively. Moreover, the finite energy resolution of the spectroscopic data needs to be taken into account by performing comparisons with calculations only after these have been properly broadened by the experimental resolution function.

As mentioned, the present work focuses on HD. It is well known that the low mass (between 2 and 6 a.m.u.) and relatively low temperatures (e.g., around 20 K) of molecular hydrogen and its isotopes in the liquid phase make the de Broglie thermal wavelength [25] $\Lambda =h/\sqrt{2\pi M{k}_{\mathrm{B}}T}$ reach values of the order of the molecular size, while remaining inferior to the average intermolecular distance [26]. Therefore, the quantum delocalization of individual particles affects the static and dynamic CM properties of these systems with respect to classical behavior, while indistinguishability can still be assumed to play a negligible role, thus justifying the use of Boltzmann statistics.

These overall assumptions are commonly summarized by saying that hydrogens are moderate quantum fluids, if compared to the paradigmatic case of helium. From a practical point of view, such a mild quantum nature has been the rationale behind the undiscouraged development of simulation algorithms still based on the possibility of defining trajectories in the phase space, but aimed at capturing the nonclassical effects of particle delocalization, at least on the simplest time correlation functions relevant to fluid dynamics. In this respect, several positive results were gathered in recent decades about the effectiveness of Centroid Molecular Dynamics (CMD) [27,28,29] and Ring Polymer Molecular Dynamics (RPMD) [30,31,32] simulation methods for the prediction of the CM velocity autocorrelation function (VAF) of the hydrogen homonuclear liquids. Among these confirmations, in the present context, the good performance of RPMD or CMD VAF calculations in estimating the total scattering neutron cross section of both ${\mathrm{H}}_2$ [3] and ${\mathrm{D}}_2$ [4] is of special relevance. In particular, the above algorithms were used to obtain quantum compliant evaluations of the VAF, which, combined with the Gaussian Approximation (GA) [33,34], are able to provide, at present, the most reliable determination of the CM $S_{\mathrm{s}}(Q,\omega )$ of the mentioned homonuclear liquids. A summary of the GA for $S_{\mathrm{s}}(Q,\omega )$ is provided in Appendix B.

Concerning the total dynamics, unfortunately no analogous way exists to derive accurate quantum-compliant CM $S(Q,\omega )$ from RPMD simulations. Indeed, such semiclassical techniques do not properly capture the quantum character of the system when dealing with nonlinear-operator time correlation functions like $F(Q,t)$, as shown in recent years by experimental determinations of the total dynamics of liquid ${\mathrm{D}}_2$ [7]. However, as we have shown in the case of ${\mathrm{D}}_2$ DDCS calculations [4], a simple, yet effective, solution is to use the Sköld approximation [35], which provides the total $S(Q,\omega )$ as a suitable modification of its self part alone, namely:

$$S(Q,\omega )\approx S\left(Q\right)S_{\mathrm{s}}\left(\frac{Q}{\sqrt{S\left(Q\right)}},\omega \right),$$

(8)

where $S\left(Q\right)$ is the CM static structure factor. In the absence of experimental values of $S\left(Q\right)$, a valid alternative is to resort to simulations of this quantity, which is also a by-product of the Path Integral (PI) RPMD calculations.

For the reasons explained above, we therefore performed a PI RPMD simulation for the case of HD at a temperature $T=$ 17 K (rather close to the HD triple point temperature $T_{\mathrm{tr}}$ = 16.60 K) and number density n = 24.37 ${\mathrm{nm}}^{-3}$. Note that $pVT$ data are not available for HD, with the only tabulated values being those on the ($n,T$) liquid–vapor coexistence curve reported in Ref. [36]. The chosen temperature is the one where the experimental total cross section data for HD are available [22]. The open source molecular dynamics code, i-PI v2.0 [37], was used for the RPMD simulations.

The simulations were carried out using a Trotter number (number of beads on each polymer) $P=64$, timestep of 1 fs, and with a system of $N=512$ polymers interacting through the Silvera–Goldman potential [38], which is the state-of-the-art interaction model for the hydrogen liquids. The system of polymers was equilibrated in the NVT ensemble for 50 ps using the white noise Langevin thermostat option “pile_g” of the software and with the coefficient $\tau =10$ fs. The VAF was calculated in the NVE ensemble up to a maximum time lag of 2.5 ps. For the calculations of $S\left(Q\right)$, the same equilibrium procedure was followed; however, the number of polymers was set to N = 4096. The pair distribution function $g\left(r\right)$ was obtained through the TRAVIS code [39,40] using the trajectory files of the simulations, and then used to compute $S\left(Q\right)$. Together with the canonical VAF, it was then possible to implement both the GA and the Sköld approximation in the neutron scattering law summarized in Equation (7).

3. Results

The results of the PI RPMD simulation of HD are reported in Figure 1, where both $S\left(Q\right)$ and the canonical VAF (see Appendix B) are displayed. The latter shows the expected behavior for a dense liquid, with a negative minimum witnessing the rattling motion of the tagged molecule within the temporary cage formed by its neighbors. In Appendix B, we show how an accurate estimate of the self diffusion coefficient D can be obtained from the simulated VAF, finding the value D = 3.2 $\times {10}^{-3}$ ${\mathrm{nm}}^2$ ${\mathrm{ps}}^{-1}$ for HD at the mentioned thermodynamic state. This is a significant outcome of this work, given the lack of literature values for D.

Specific calculations of the form factors of HD and of Equation (7) within the GA and Sköld approximations were performed by using the molecular parameters and neutron scattering lengths listed in

Table A1. Basic quantities used in the present calculations for HD.

Parameter	Description
B = 5.538 meV [44]	Rotational constant in the vibrational ground state
D = 0.003263 meV [44]	Centrifugal distortion coefficient in the vibrational ground state
$\hslash {\omega }_{\mathrm{v}}$ = 450.38 meV [44]	Quantum of vibrational energy
$r_{\mathrm{eq}}$ = 0.74142 Å [44]	Equilibrium internuclear distance
$m_1$ = 1.00794 a.m.u.	Mass of the proton
$m_2$ = 2.01410 a.m.u	Mass of the deuteron
${\gamma }_1=\frac{2}{3}$	Fraction of the internuclear distance pertaining to the H nucleus
${\gamma }_2=\frac{1}{3}$	Fraction of the internuclear distance pertaining to the D nucleus
$b_{\mathrm{coh},1}$ = −3.7406 fm [41]	Coherent scattering length of the H nucleus
$b_{\mathrm{inc},1}$ = 25.274 fm [41]	Incoherent scattering length of the H nucleus
$b_{\mathrm{coh},2}$ = 6.674 fm [41]	Coherent scattering length of the D nucleus
$b_{\mathrm{inc},2}$ = 4.033 fm [41]	Incoherent scattering length of the D nucleus

of Appendix A.

In Figure 2, we show the general features of the inter- and intramolecular form factors, $u\left(Q\right)$ and $v(Q,0)$, as a function of the exchanged wavevector Q. The various terms contributing, respectively, to Equation (A17) and, at $t=0$, to Equation (A12), are also displayed separately. Note that while $u\left(Q\right)$ does not depend on the incident energy, $v(Q,0)$ may vary with $E_0$ because, in practical implementations of the summation of Equation (A12), the number of rotovibrational transitions grows with incident energy. The example of Figure 2 refers, in particular, to an incident energy $E_0$ = 80 meV. As evident from Figure 2, the intramolecular form factor of HD overwhelms the intermolecular one, as expected.

Examples of the DDCS spectra according to Equation (7) for $E_0$ = 80 meV are shown in Figure 3 at three selected Q values. The typical asymmetry of the spectra of quantum systems, as well as the $J_0=0\to J_1=1$ rotational line at about $E=\hslash \omega =$ 11 meV, can be better appreciated with a growing Q. Instead, due to the rather low temperature, the detailed balance strongly depletes the anti-Stokes line $J_0=1\to J_1=0$, which is present despite the low population (the $J_0=1$ level of HD at 17 K corresponds to a 0.16% occupation probability), but not appreciable on the scale of the figure. The example value chosen for $E_0$ allows to span a reasonable exchanged-energy range, and is a plausible choice if one were interested in performing a DDCS measurement on liquid HD, similarly to what was performed in Ref. [7] for liquid ${\mathrm{D}}_2$. Of course, as mentioned in Section 2, before comparison with any (forthcoming) experimental DDCS result, the spectra of Figure 3 should be convoluted with the energy resolution function of the used instrument.

It is worth stating, finally, that the total DDCS of Figure 3, i.e., including the distinct dynamics, differs negligibly, at all Q values and at all incident energies, from its self component alone (not shown in the figure), as predictable from the very small values of $u\left(Q\right)$ of Figure 2. This simply means that HD is a “true” predominantly incoherent scatterer at all incident energies, even more than para-${\mathrm{H}}_2$. Indeed, the absence of spin correlations in HD does not introduce different scattering cross sections according to the parity of the levels involved in the transitions and the spin statistics of the nuclei in the molecule. Conversely, in para-${\mathrm{H}}_2$, if $E_0$ is insufficient to excite the first rotational line $J_0=0\to J_1=1$, then only the even–even $J_0=0\to J_1=0$ transition contributes, and spin correlations can be found [19] to weight such kinds of transitions only by the coherent cross section. Therefore, at some incident energies ($E_0<14.7$ meV), para-${\mathrm{H}}_2$ behaves coherently and the distinct dynamics plays a role. In HD, this is never the case, and the huge incoherent cross section of H dominates the scattering behavior independently of $E_0$.

In order to make a comparison with the experimental total cross section possible, the DDCS of HD was calculated at $T=17$ K as a function of scattering angle $\theta$ and exchanged energy $E=\hslash \omega$ at various values of the incident energy $E_0$ ranging between 1 and 85 meV. At such energies, only rotations are excited and the main contribution to the spectra comes from the elastic $J_0=0\to J_1=0$ line and the mentioned $J_0=0\to J_1=1$ Stokes transition (when excited).

Double integration of the DDCS over the exchanged energy and solid angle, according to:

$${\sigma }_{\mathrm{s}}={\int }_{\Omega }d\Omega {\int }_{-\infty }^{{\omega }_0}d\omega \frac{d^2\sigma }{d\Omega d\omega },$$

(9)

provides the total scattering cross section ${\sigma }_{\mathrm{s}}$ shown in Figure 4, in comparison with Seiffert’s data [21], from which the known absorption cross section ${\sigma }_{\mathrm{a}}\left(E_0\right)$ [41] has been subtracted: in formulas, ${\sigma }_{\mathrm{s}}^{expt.}={\sigma }^{expt.}-{\sigma }_{\mathrm{a}}$. The effect of the neglect of the distinct part in the calculation of ${\sigma }_{\mathrm{s}}$ is observed to be very small (compare the blue circles with the red stars in Figure 4), and slightly appreciable only at very low $E_0$ values.

Unexpectedly, agreement between the data and calculations is not as good as the one found for normal ${\mathrm{H}}_2$ and ${\mathrm{D}}_2$ [4,20]. In Ref. [21], the author reports an average systematic error of $\pm 8.1\%$ and an average statistical error of $2.8\%$ for HD. Figure 5 shows the possible effect of shifting Seiffert’s data upwards by the systematic error (empty circles), while statistical uncertainties are taken into account by the corresponding error bars. By doing so, a good agreement is recovered at all energies (given that the experimental data points are rather scattered), except in the low-energy range $1<E_0<20$ meV, a disagreement which is worth investigating and discussing.

4. Discussion

Trying to explore the reasons for the mismatch in the subthermal region, a possible effect of sample contamination might be considered. Information on this can be found in another paper by Seiffert [22], which reports the following composition determined by mass spectrometry: 94% HD, 5.5% ${\mathrm{H}}_2$ and 0.5% ${\mathrm{D}}_2$. However, it is not specified whether the systematic error estimated by the author already included the role played by the impurities, or if these have to be considered as an additional source of error. Moreover, normal ${\mathrm{H}}_2$ cannot be assumed to be part of the mixture, since even a small amount would provide total cross section values higher than those calculated for pure HD. This can be understood by considering the cross section monotonic behavior reported for n-${\mathrm{H}}_2$ in Refs. [20,21]. Conversely, the much lower cross section of para-${\mathrm{H}}_2$ in the range $1<E_0<20$ meV might well explain the discrepancy.

Preliminary calculations for the pure HD and the mentioned mixture were performed in past years using the Gaussian lineshape of the ideal gas for $S_s(Q,\omega )$, finding an overall disagreement even by accounting for the presence of impurities [42], and which could not be recovered, even by adding the declared systematic errors. Here, however, we also adopted the more refined RPMD-based GA approach for HD, so it is worth analyzing the possible improvements due to the model lineshape. Although the available ${\sigma }_{\mathrm{s}}$ data for ${\mathrm{D}}_2$ and para-${\mathrm{H}}_2$ are not exactly at 17 K, we combined the various calculations (all based on the GA and on quantum evaluations of the VAF) according to the declared composition, and report the results in Figure 6. The agreement is now very good at almost all energies. These results confirm both the sample contamination (likely not included in the estimate of systematic errors) but, more importantly, the improvements introduced by the use of quantum-compliant lineshapes, with respect to the earlier ideal gas modeling [42].

Nonetheless, the impressive result shown in Figure 6 has been obtained by assuming a well-defined direction of the systematic effects: an assumption the correctness of which we are, of course, not able to prove. A separate preliminary study [43], based on combining the earlier work of Ref. [14] with the approaches described in Ref. [6], indicated, on average, a roughly 10% lower cross section between 10 and 80 meV compared to the results presented here. Although the source of this discrepancy is under detailed investigation, it anyway emerges very clearly that new accurate ${\sigma }_{\mathrm{s}}$ measurements, as those performed in recent years for para-${\mathrm{H}}_2$, are also highly auspicable for HD.

5. Conclusions

This work completes the picture of the calculation of the neutron double differential cross section of diatomic molecules, providing formulas for a heteronuclear vibrating rotor at low temperature. The interest for a treatment of the heteronuclear case, and particularly for hydrogen deuteride, was triggered by possible applications of this fluid in neutron moderation. Unfortunately, the lack of experimental DDCS data for HD prevents one from a stringent test of the calculations. It is, therefore, clear that our capability to predict the neutron DDCS of this cryogenic liquid and the setting up of reliable scattering libraries would profit from accurate inelastic scattering experiments and improved quantum simulation methods. Similar importance would have experimental work aimed at the determination of the thermodynamic properties of HD, at least in the low temperature dense liquid region useful in applications.