Coherent Neutron Scattering in Polycrystalline Deuterium and its Implications for Ultracold Neutron Production

# Coherent Neutron Scattering in Polycrystalline Deuterium and its Implications for Ultracold Neutron Production

## Abstract

This paper presents a calculation of the neutron cross-sections in solid materials (used in practical neutron sources) with a large coherent scattering contribution. In particular, the dynamic structure function S(Q, ) of polycrystalline ortho-D is evaluated using a Monte-Carlo calculation that performs an average over scattering angles relative to crystal axes in random orientations. This method uses an analytical dispersion function with force constants derived from neutron scattering data of single crystal D in the framework of an axially symmetric force tensor. The resulting two dimensional map of S(Q, ) captures details of the phonon branches as well as the molecular rotations, that can be compared directly to data from inelastic neutron scattering on polycrystalline D. This high resolution information is used to calculate the absolute cross-sections of production and upscattering loss of ultracold neutron (UCN). The resulting scattering cross-sections are significantly different, especially for UCN upscattering, from the previous predictions using the approach centered on the incoherent approximation.

Ultracold Neutron; UCN; Solid Deuterium; Coherent Scattering; Incoherent Approximation

## I Introduction

Neutron moderation is a well established technique, developed for applications in fission reactorsGlasstone and Edlund (1952) and nuclear weapons. Interest in this topic has been revived with the development of new neutron scattering facilities for investigation of materials, new sources of ultracold neutrons for fundamental physics, as well as new deep-underground laboratories which require extensive shielding. In the latter case, neutrons generated in the surrounding rock via the spallation process initiated by high energy cosmic rays are a major source of background. On the pure academic front, neutron moderation can be put in a broader context of phase space compression, developed for several different types of particle beams in various areas of physics, from stochastic electron cooling in high energy particle accelerators on Tera-eV scales to laser cooling of neutral atoms to form Bose-Einstein condensates on pico-eV scales. The phase space compression of neutrons is especially challenging. Because the electromagnetic interaction of the neutron is relatively weak due to its zero charge and small magnetic moment, the usual mechanisms of phase space compression for charged particles do not apply. Instead, the nuclear force couples neutrons to individual nuclei in the interacting medium. Upon scattering, the nuclear recoil leads to rapid moderation for fast neutrons, but stops to be effective for neutrons with energy below a few tens of milli-eV. For neutrons of particular interests to fundamental research, such as cold, very-cold, and ultra-cold neutrons (with energies ranging from milli-eV to nano-eV), the only efficient way to date to increase neutron phase space density is through inelastic collisions in a cold medium, so-called neutron moderator.

To study the process of neutron cooling and estimate the neutron flux for the broad range of applications mentioned above, standard simulation codes are available including the widely-used Monte Carlo N-Particle Transport Code (MCNP)Labratory (2008). The code was developed in the 1960s at Los Alamos National Laboratory, and has since been extensively bench-marked with experimental data; the physics models and scattering kernels used have been refined with every new release. In MCNP, moderations of neutrons in a medium are modeled by the Boltzmann equation that details the evolution of the phase space distribution function of neutrons. The Boltzmann equation includes neutron diffusion, neutron absorption loss, and inelastic processes in which neutrons downscatter with energy loss or upscatter with energy boost. In this formalism, the scattering kernel determines the probability for inelastic processes. For fast neutrons (with energies from eV to MeV), nuclear resonances strongly influence the scattering and absorption processes. Away from the resonances, the neutron energy loss through nuclear recoil is calculated with a simple free gas model, in which the time scale of collisions is so short that nuclei are treated as free. For neutrons with wavelength larger than the separation between scattering nuclei, however, thermal neutron scattering kernels that include coherent diffraction and collective excitations in bulk materials have only been incorporated for limited materials MacFarlane (1994). Whereas nuclear cross-sections for fast neutron reactions are detailed for almost every stable isotope, there exist only a handful of materials with extensive treatments that include condensed matter effects beyond the simple nuclear recoil in the free gas model.

Hydrogen and deuterium in liquid and solid forms are of particular interest, as they are the materials of choice to construct cold neutron moderators used in research reactors and spallation neutron sources. The large nuclear cross-section, together with the maximum energy loss in recoil kinematics of n-H scattering, leads to rapid cool-down of incident neutrons. However, the H(n,)H reaction results in significant loss of the population of the already moderated, low-energy neutrons. Neutron absorption loss in H (i.e., D) is much less of an issue, so even though the cross-section is an order of magnitude smaller, D compares favorably to H as a good neutron moderator. Since the ground state of hydrogen forms a diatomic molecule, the scattered wavefunctions of the incident neutron off the two identical nuclei within a single molecule interfere coherently. Together with the spin-dependent nuclear force, the inter-molecular interference leads to different scattering amplitudes between states with distinct nuclear spin. Furthermore, spin statistics modify the molecular form factor, giving rise to selection rules governing the coupling to discrete sets of rotational states Young and Koppel (1964). In fact, applying Fermi-Dirac statistics in the diatomic H molecule was the key to solve the puzzle of discrepancies between measured values of nuclear cross-section of neutron-proton scattering and theoretical predictions using the nuclear bound state of D Schwinger (1940); Hamermesh and Schwinger (1946), during the early construction of nuclear structure theory. To this end, most cross-section evaluations have properly included this coherent spin effect MacFarlane (1994).

On the other hand, all reported works (including NJOY njo (1999), Liu et al. (2000), Atchison et al. (2007), Granada (2009)) that evaluate inelastic cross-sections in the thermal and cold regimes for use in MCNP employ the incoherent approximation (IA) that assumes negligible interference between wavefunctions scattered from different lattice sites in the intra-molecular scale. In the IA, the amplitude of coherent scattering is evaluated using the same algorithm used in evaluating the incoherent scattering. The use of IA is popular because the inelastic incoherent cross-sections can be readily calculated via standard phonon expansions, only requiring limited information on density of states (or the frequency distribution, as referred in some literature) of the collective excitations as an input. While IA works quite well for moderators consisting of hydrogenous materials (=1.7583 b, =80.27 b), in the case of D, the contribution of the coherent scattering is not small (=5.592 b, =2.050 b) and thus the adoption of IA should be examined carefully. Especially for low energy neutrons with wavelengths comparable to and larger than the lattice constants, the use of IA seems particularly deficient. To address the shortcomings of this assumption and all previous calculations that invoked this assumption, we have developed the first full treatment of neutron scattering for neutrons with wavelength larger than angstroms. This calculation includes both the coherent and the incoherent scattering which excite single and multiple phonons, as well as the processes which excite rotational transitions. The elastic Bragg scattering is also included. In this paper, we will show that this full model, when applied to solid D, significantly modifies the cross-sections of ultracold neutrons (UCN), which are free neutrons with E/ 3 mK, where is the Boltzmann constant.

The scattering process responsible for UCN production differs from the thermalization process in the thermal and cold neutron moderators, in which multiple scatterings take place while the energy spectrum evolves continuously until the neutron gas establishes thermal equilibrium. In contrast, superthermal UCN production occurs as a single scattering event during which the incident neutron comes to a near-full stop, giving up its energy and momentum to a matching quasi-particle created in the interacting media. Even though the phase space of this inelastic scattering process is limited, the cross-section is non-negligible, provided that the interacting medium has excitations that coincide with the energy and momentum of incident cold neutrons. Subsequent scattering of the down-converted UCN (mostly upscatterings) is suppressed inside a superthermal UCN source by simply reducing the thermal population of the quasi-particles. Cooling the moderator effectively pumps out these quasi-particles. For the typical application of neutron thermalization down to a few milli-eV, scattering kernels evaluated based on IA give reasonable predictions of the neutron yields. However, many criteria of IA break down when applied to the extreme low-energy regime (E 350 neV) of UCN physics. Re-evaluating the UCN cross-sections in solid ortho-D using the full model shows that UCN cross-sections differ significantly from previous estimates using the IA approach Liu et al. (2000).

## Ii Dynamics in solid D2

Calculation of the coherent inelastic scattering in solid D requires detailed knowledge of the dynamics and the dispersion of the collective energy excitations. Fortunately, the dynamics of D in a molecular lattice were studied extensively in the early days of neutron scattering experiments Squires and Stewart (1955); Egelstaff et al. (1967); Elliott and Hartmann (1967); Schott (1970); Diehl and Biem (1975); Danchuk et al. (2004); Nielsen and MÃ¸ller (1971); Nielsen (1973); Schmidt et al. (1984). Due to the large lattice constants resulting from the large zero-point motion, molecular rotation remains free in the solid matrix and the translational degrees of freedom can be decoupled and the coupling strength estimated using an isotropic Lennard-Jones potential. The crystal is simple enough that the dispersion energy and the polarization vector of distinct phonon branches can be calculated with a tensor-force model, without the use of molecular dynamic simulations that require intensive numerical computations. Dewames et al.  Dewames et al. (1965); Lehman et al. (1962); Collins (1962) developed a Born-von Karman model, using an axially symmetric (A-S) two-body potential including up to the third nearest-neighbors, for systems with the hexagonal close packed (HCP) structure. This dynamical matrix model was applied successfully to HCP terbium Houmann (1970), and recently extended to HCP Tb, Sc, Ti and Co Vaks and Khromov (2008). Using the force constants measured by Nielsen Nielsen (1973) in the A-S model framework, we re-constructed the full energy dispersion that is a function of the three dimensional momentum vector of phonons propagating in solid ortho-D. This angular-dependent dispersion function is then used in a Monte-Carlo code to construct the inelastic coherent scattering cross-section for polycrystalline solid D. For each scattering event sampled by the Monte-Carlo, the scattering kinematics is determined and the coherent scattering amplitude is calculated following the standard textbook formalism Lovesey (1984). The major difference in evaluating the coherent and incoherent scattering cross-sections lies in the kinematic conditions of momentum conservation. For the incoherent process, the scattered wavefunctions are not summed coherently as the phase coherence is destroyed by fluctuations of the scattering sites. This leads to a relaxation of the momentum conservation law typically applied to two-particle systems. Instead the approach of IA uses the density of states of the phonon modes to estimate the relative contribution of the scattering amplitude for phonons of different energy. The energy of phonon equates to the energy transfer of the scattered neutron a result of energy conservation. As a consequence, cross-sections evaluated using IA has a diffuse smooth dependence without any localized enhanced intensity, as would have been expected in coherent scatterings.

Solid hydrogen and deuterium are considered quantum solids due to the light mass of individual molecules. The quantum nature is determined by its thermal wavelength, defined as  Lovesey (1984). For D, it is 4.6 Å at 5 K and 3 Å at 18 K. For a light molecule (and atom), its large thermal wavelength exdending beyond the inter-molecular spacing causes direct overlap of the wavefunctions of adjacent molecules. This could lead to quantum inteference. Provided that the adjacent molecules are identical particles, the Bose-Einstein statistics could significantly modify the collective behavior of the molecules, as in the case of superfluid helium. Here, the ortho-state of D molecule has possible nuclear spin of 0 and 2 and distinct sub-states of =(0,0), (2,2), (2,1), (2,0), (2,-1), (2,-2) that are equally populated. Given that the seperation between molecules on the HCP lattice is , the mean seperation of particles with identical states is , which is larger than the thermal wavelength by 40%. This suggests that at temperatures higher than 5 K, solid D can be treated as a classical system to the first order approximation. Many ongoing research Danchuk et al. (2004); Colognesi et al. (2009) are looking for the evidence of quantum interference in H, HD and possibly D.

In a classical system, the double differential cross-section of neutron scattering is conventionally expressed as:

 d2σdΩdε=kk0S(→Q,ω), (1)

where arises from the ratio of phase space before and after scattering, and the dynamic structure function, , contains all the detailed physics of the neutron interactions. The dynamic structure function is related to the scattering kernel used in the neutron-moderation Boltzmann equation, in which the rate of neutron moderation is directly proportional to the dynamic structure function of the interacting medium. For the remaining of the paper, we examine neutron scattering in a classical solid lattice, where the scattering can be described as

 S(→Q,ω) = b2e−2W{∑→Gδ(→Q−→G)δ(ω)+∑j,→qℏ2(→Q⋅→ej)22Mωj(→q) (2) [∑→G(nj(→q)+1)δ(→Q−→q−→G)δ(ω−ωj(→q)) +∑→Gnj(→q)δ(→Q+→q−→G)δ(ω+ωj(→q))] +O(Q4)+...}.

The scattering length is determined by the strength of interaction between the neutron and the scattering site. The Debye-Waller factor is a result of the de-localization of the atoms on the lattice sites due to thermal motion (as well as zero-point motion), smearing out the sharpness of constructive interferences of the scattered wavefunctions. On a lattice with simple cubic symmetry, at low temperatures, reduces to

 2W→34ℏ2Q22MωDas T→0, (3)

where is the mass of the lattice site. For solid D, is the total mass of the D molecule. The Debye temperature for solid D is around 110 K.

The first term in Eq.(2) is the coherent elastic scattering, during which the energy transfer of the neutron is zero but the momentum transfer can be absorbed by the whole lattice if the Bragg condition () is satisfied. The next term describes the inelastic scattering where the neutron energy is released to create a single phonon in the mode with energy and polarization vector . The occupation number of the created phonon is , which satisfies the Bose-Einstein statistics. The delta functions enforce energy and momentum conservation. Any change in momentum of the incident neutron during collision is picked up by a phonon of momentum . The total momentum conservation is adjusted by the appropriate inverse lattice vector if the change of momentum is beyond the first Brillouin zone. The third term leads to energy and momentum gain by the neutron through absorption of a single phonon. Higher-order terms describe multi-phonon processes, with rapidly decreasing probabilities for small momentum transfer . Finally, with D molecules arranged in a solid lattice with lattice vector , the inverse lattice vectors along the principle axes are

 →g1=2π→a2×→a3V,% →g2=2π→a3×→a1V, →g3=2π→a1×→a2V.

where is the volume of the unit cell.

## Iii Elastic Scattering and Form Factors

Even though the elastic scattering is well understood, comparing the experimental data of cross-sections with the calculation allows for consistency checks on the occasionally non-trivial molecular form factors. This is particularly useful when calculating cross-sections in solid oxygen, where the spin form factor is less well-known than the form factors of atomic nuclei, with additional complications due to preferential alignments of the molecular axis relative to the crystal axes. For D, the situation is somewhat simpler because the rotational motion remains free in the solid lattice, and the molecular wavefunction can be described by spherical harmonics without preferred orientations Young and Koppel (1964). The total cross-section for elastic scattering (integrating the first term in Eq.(2)) can be separated into coherent Bragg part and incoherent diffuse part Egelstaff and Pease (1954); Lovesey (1984):

 σcoh0 = 4πλ28πV∑hkldhkl|Fcohhkl|2e−2W, and (4) σinc0 = 4π|Finc|2e−2W, (5)

where is the inter-planar spacing between the [hkl] planes in the lattice space, i.e.,

 dhkl=2π|→G|, where →G=h→g1+k→g2+l→g3. (6)

The form factor within the unit cell differs depending on whether the coherence of the scattered wave is preserved after scattering. For coherent processes,

 Fcohhkl = Nbasis∑i=1bcohiei→G⋅→Ri if λ⩾2dhkl (7) = 0 otherwise.

For incoherent processes,

 Finc =  ⎷Nbasis∑i=1(binci)2 . (8)

In a non-Bravais crystal, the number of atoms in the unit cell, , is larger than 1.

On each lattice site in the D crystal resides a diatomic molecule, the additional molecular form factor of which significantly modifies the scattering length. For D, the bosonic nature of the two identical deuteron nuclei leads to unique couplings to a set of molecular rotational states of either even or odd symmetry, so that the symmetry of the whole system is preserved under permutation of identical particles. It follows that there exist two different types of D molecules, i.e., ortho-D with even nuclear spin (=0, 2) coupling to symmetric molecular rotational wavefunctions (with orbital angular momentum =0,2,4,…), and para-D with odd nuclear spin (=1) coupling to anti-symmetric molecular wavefunctions (of =1,3,5,..). By convention the ortho-state designates the group with the highest spin multiplicity. On account of the spin-dependence of the nuclear force, neutron scattering lengths of ortho and para-D are different. In addition, when interacting with one species of diatomic molecule, say ortho-D molecules (with =0,2 and =0 in the ground state), the neutron scattering length also depends on the final state of the target molecule Young and Koppel (1964), i.e.,

 b2D2∣∣J=0→0 = (a2coh+58a2inc)C2(000;00)|A00|2, (9) b2D2∣∣J=0→1 = 38a2inc(2⋅1+1)C2(011;00)|A01|2% , (10) b2D2∣∣J=0→2 = 38a2inc(2⋅2+1)C2(022;00)|A02|2% , (11)

where  fm and  fm, and the molecular form factors are:

 A00=2j0(Qa2),A01=2j1(Qa2),A02=2j2(Qa2), (12)

with as the spherical Bessel function of order . The separation of the two deuteron nuclei is . All Clebsch-Gordan coefficients used in the above equations are 1. The process resulting in the ortho- to para- transition requires that the neutron to transfer energy equal to the rotational transition energy, and thus it is only present in the inelastic scattering.

The energy dependence of the total elastic scattering cross-section of polycrystalline samples is determined by the static structure function of single-crystal, which is HCP for D Bostanjo.O and Kleinsch.R (1967). However, some Raman scattering data Collins et al. (1996); Stein et al. (1972) indicate an FCC component depending on the temperature and the pressure of the solid. Fig. 2 shows the calculations together with experimental data. Note the molecular form factor and the Debye-Waller factor together significantly reduce the elastic scattering amplitude, as both suppress the neutron scattering amplitude associated with large momentum transfers (see Fig. 5 and more discussion in Sec.IV). Note also that the experimental data do not show the large diffraction peak due to scattering off the [011] plane. There is evidence that the solid D grown in a cold cell tends to self-anneal, becoming a few single crystals of large size. In these experiments, neutrons scatter through diffractions from several single crystals of large size oriented at different angles. As a result, the total cross-section deviates from the prediction in the powder limit described by Eq.(4).

The presented total cross-section is extracted from neutron transmission measurements. Excluding the [011] scattering alone brings the prediction on top of the experimental data (see the light blue dashed curve in Fig. 2). In these measurements, the finite area of the neutron detector would lead to reduced scattering amplitudes, because neutrons scattered at low angles are still detected. However, the surface coverage of the detector has to be unphysically large in order to account for the missing amplitude. This rules out detector geometry as the cause of the missing scattering amplitude. On the other hand, missing the major [011] peak suggests that the crystal is grown preferably with the c-axis perpendicular to the direction of the neutron beam. This would lead to reduced Bragg scattering for small momentum transfer centered around the direction of the c-axis. The reduced scattering cross-section could also be explained by multiple scattering processes, in which some neutrons that are already scattered once scatter again back into the beam and thus enhances the neutron transmission. The presence of multiple scattering leads to a reduction in the cross-section if the analysis algorithm does not include a correction for multiple scattering. Due to this complication, the data of neutron transmission, which was collected along with the data of UCN production, cannot be used to distinguish between the HCP and FCC lattice structures. Finally, the discrepancy between the experimental data and the theoretical predictions at neutron energies higher than 10 meV is accounted by including contributions from inelastic scattering, as explained in the following section.

## Iv Inelastic Scattering

The inelastic scattering involving creations and annihilations of single phonon in Eq.(2) is calculated numerically for coherent and incoherent processes with corresponding scattering lengths given by Eqs.(9) and (11). To evaluate coherent scattering in a polycrystalline sample beyond the IA prescription, we developed a Monte-Carlo algorithm to sample the dispersion curve of each single crystallite oriented at random angles, strictly following the kinematics of energy and momentum conservation in each coherent scattering event Liu and Young (2004). Taking the angular average of the momentum transfer vector simplifies the orientation-dependent scattering amplitude into a two dimensional map of S(Q, ). The result of one such calculation (a high resolution map with 500200 grids in the (Q, ) space) is shown in Fig. 3. The calculation includes 1-phonon and 2-phonon contributions from both the coherent and incoherent process, calculated independently for the ortho ortho () and the ortho para () transition. Note that the neutron scattering resulting in ortho para transition involves spin flip, and thus is a purely incoherent process. The relative weight of the scattering amplitude of the ortho-ortho and ortho-para processes is dictated by Eqs.(9) and (11). In addition, the calculation also includes inelastic scattering leading to the rotational excitation , where incident neutrons lose energy in ortho-D through inter-molecular excitations without activating translational phonons via a nuclear recoil. In upscattering, this process is absent in ortho-D for neutrons with energy less than the transition energy of 7.1 meV.

The average algorithm used in polycrystalline samples prescribes that, for every given magnitude of momentum transfer, a random angle is assigned to determine the momentum transfer vector. Setting this momentum transfer to that of the phonon, the eigenfrequencies and eigenvectors are calculated uniquely using the force tensor matrix. For the HCP structure, the result consists of six eigenvalues, each corresponding to a translational degree of freedom originated from the two molecules in a unit cell. The energy transfer bin that contains the calculated eigen-frequency is then augmented by a numerical value of the scattering amplitude estimated using Eq.(2). For each magnitude of momentum transfer, up to 10 different angles are sampled isotropically over . In the end, a two-dimensional map of the dynamic structure function is evaluated for positive energy transfer, which corresponds to neutron downscattering. A second map is evaluated independently for negative energy transfer (neutron upscattering). The main difference between upscattering and downscattering arises from the dependence on the phonon occupation number . Upscattering occurs when the neutron encounters a phonon, and thus the scattering probability scales with the thermal population of phonons. In contrast, downscattering of neutrons through phonon creation takes place only when the kinematics of energy and momentum transfer is favorable, as determined by the existing modes of phonon excitation. The downscattering amplitude scales as (1 + ), where weakly enhances the transition probability due to the bosonic nature of phonons.

Next, to evaluate the incoherent contributions, we construct the density of states by integrating the coherent dynamic structure function derived in the previous step, after correcting all the dynamical factors, over the whole accessible range of momentum transfer:

 Z(ω)=∫∞Q=12Qfreeωℏ2Q2/2me2W(n(ω)+1)Scoh(Q,ω)dQ. (13)

However, simply integrating over all possible values leads to a non-vanishing density of states as approaches zero. This disagrees with the Debye model, which prescribes the density of states for phonons as , where is the speed of sound in the solid. To fix this discrepancy, we set the lower bound of the integration to half of the momentum carried by a free neutron, for the following reasons: As illustrated in Fig. 3, the scattering amplitude peaks beyond the first Brillouin zone. The accessible kinematic range of free neutrons in the cold regime (with incident energy less than 10 meV) simply cannot produce excitations within the first Brillouin zone, except for those close to the zone boundaries. In other words, the sound speed of phonons exceeds the group velocity of the neutron wave-packet, leading to no coupling between neutrons and phonons of low momentums. Excluding the rather large amplitude for coherent creation of acoustic phonons with momentum smaller than half of the free neutron momentum (in the first Brillouin zone) restores the expected dependence at low momentum transfers (as shown in Fig.4). Only the coherent process is included in evaluating .

We have also calculated using the classical approach by integrating the derivative of the dispersion curve over the hypersurface with constant , following a method described in Raubenheimer and Gilat (1967). The result is consistent with the method discussed above. With , the incoherent scattering law can be calculated using the standard textbook formula, without the need for the rather computation-intensive algorithm developed for the coherent process.

For a final consistency check, using the derived density of states, we can calculate the Debye-Waller factor using the polarization vectors in the HCP structure Lovesey (1984):

 W(Q)=3∫ωm0dωZ(ω)ℏ24Mω(2n(ω)+1)|→Q⋅→e|2av, (14)

and examine how well the approximation of the cubic symmetry of Eq.(3) applies to the HCP structure. We found that the Debye-Waller factor of solid D at 18 K agrees remarkablly well with the simple cubic approximation, and that at 5 K is slightly above the value of the cubic structure (see Fig.5). For the range relevent to UCN production (Q ), the agreement is within 5%. For completeness, the molecular form factors are plotted to compare with the Debye-Waller factor. For the transition, the form factor suppresses the high momentum transfer, similiar to the effect of the Debye-Waller factor. For the and transitions, the form factors start at zero and grow with increasing . Processes involving the transition (E15 meV) are highly suppressed in magnitudes and thus we don’t expect it to contribute significaly to the cross-section.

After combining independent results for coherent and incoherent cross-section, the major features of the inelastic scattering are clearly visible in the dynamic structure function evaluated with our full model (see Fig. 3). The acoustic phonon branches that extend from elastic Bragg peaks into high energy transfers can be clearly identified. Only phonons created along the symmetry axes have energy dispersion that is periodic in momentum; phonons created along other directions have less symmetry, but their energy can be calculated analytically using the force tensor. We find that the clustering of the scattering amplitude around 5 meV is due to mostly the phonon modes in the basal plane. Including scatterings with directions out of the basal plane, the phase space increases resulting in large scattering amplitude associated with these in-plane phonon modes. The narrow band of large scattering amplitude around 9 meV is due to the longitudinal optical branch of phonons associated with scatterings along the c-axis. Multi-phonon contribution increases the scattering amplitude for energies larger than 10 meV and momentum transfer larger than 3/.

Overall, the scattering amplitude is significantly suppressed for momentum transfers larger than 5/, owing to the combined effect of the molecular form factor and the Debye-Waller factor. The same suppression due to these form factors was observed in the elastic scattering as discussed in the previous section. In contrast to the interpretations given in Gutsmiedl et al. (2009), our study shows that the large scattering amplitude cannot be simply associated with phonons on the zone-boundaries along the symmetry axes, because the phase space for interaction with these phonons is very small compared with that available for inelastic scattering. Finally, the pure transition (without phonon creation) is visible at a fixed energy transfer of 7.1 meV. The molecular form factor dictates the dependence of the scattering amplitude in the transition, as it increases from zero with increasing . By comparison, the dynamic structure function calculated in the IA approach (as shown in Fig. 1) is grossly over-simplified.

In order to account for the scattering amplitude beyond the energy of the most energetic phonon available (i.e., the transverse optical phonon moving in the [001] direction), we have also calculated the 2-phonon contributions. The coherent dynamic structure function involving two phonons can be expressed as:

 12∑j,→q,j′,→q′ℏ2(→Q⋅→ej)22Mωj(→q)ℏ2(→Q⋅→ej′)22Mωj′(→q′)∑→G[(nj(→q)+1) ×(nj′(→q′)+1)δ(→Q−→q−→q′−→G)δ(ω−ω(→q)−ω(→q′)) +(nj(→q)+1)nj′(→q′)δ(→Q−→q+→q′−→G)δ(ω−ω(→q)+ω(→q′)) +nj(→q)(nj′(→q′)+1)δ(→Q+→q−→q′−→G)δ(ω+ω(→q)−ω(→q′)) +nj(→q)nj′(→q′)δ(→Q+→q+→q′−→G)δ(ω+ω(→q)+ω(→q′))], (15)

where the first term describes the creation of two phonons through a single scattering process, the second and the third term describe the simultaneous creation of one phonon and annihilation of another phonon of different energy, and the final term describes the annihilation of two phonons in a single scattering event. The combined amplitude is smaller than that of the single phonon process. Nevertheless, the two-phonon process enlarges the scattering phase space to include those at higher energy transfers, and thus neutrons with energies higher than the Debye temperature participate in neutron downscattering through multi-phonon processes.

In constructing the for the two phonon process, we implement the following steps: For a given magnitude of momentum transfer , first throw the dice (i.e., the random number generator in the Monte-Carlo code) to determine one angle. This random angle together with determines the vector of momentum transfer. Throw the dice again to determine the magnitude of the momentum of the first phonon, and again to determine the angle of this momentum vector. Once the two vectors are specified, the momentum of the second phonon is determined by the momentum conservation law. Now with the momentum vectors specified, we can calculate the energy of each of the two phonons and sum them up as the total energy transfer. The energy bin corresponding to the total energy transfer for this scattering event is incremented by the amplitude described in Eq.(15). Note that the momentum conservation and energy conservation are strictly applied in every scattering event as dictated by the delta functions in Eq.(15).

To address the effect of preferred directions of crystal orientation as indicated by the total cross-section data (see Sec.III), we have also calculated the scattering laws by restricting the scattering angle relative to the inverse lattice vectors. A few examples of the resulting dynamic structure functions are shown in Fig. 6. To accentuate the differences, only the coherent scattering contributions are plotted. Comparison with recent measurements Gutsmiedl et al. (2009) suggests that the experimental data are best described by scatterings confined to within 45 of the basal plane. The detailed information contained in these inelastic scattering maps strongly suggests that the missing peak in the total cross-section (Fig.2) could be explained by a crystal orientation where the c-axis is aligned preferentially perpendicular to direction of neutron propagation.

With the detailed dynamic structure function, calculation of the total cross-section becomes a straightforward integration of Eq.(1) over the appropriate phase space:

 σtot(E0) = 12∫dΩ∫dωkk0S(Q,ω) (16) = 12k20∫dϕ∫dω∫QupperQlowerdQQS(Q,ω).

For neutron downscattering, (upper graph in Fig.3) is integrated over the energy transfer from 0 to , where is the energy of incident neutrons. For neutron upscattering, (lower graph in Fig.3) is integrated over energy transfer from 0 to . For every , the range of integration on is bounded by the momentum transfer in forward scattering and backward scattering:

 |Qlower| = √2mℏ2∣∣√E0−ω−√E0∣∣ |Qupper| = √2mℏ2(√E0−ω+√E0). (17)

The range of integration is illustrated by the white lines in Fig.3 for several incident neutron energies. Since the dynamic structure function derived for polycrystalline target is already angle-averaged, the integration in Eq.(16) over the azimuthal angle simply yields . In the end, to report the total inelastic scattering cross-section per molecule, the result is divided by the number of molecules contained in each unit cell, which is two in HCP D (shown as the pre-factor of 1/2 in Eq.(16)). The same normalization is applied to the incoherent component. According to Eq.(8), the lattice form factor is . The form factor together with the normalization brings the incoherent amplitude back to the scattering cross-section of one particle. The total cross-section is plotted in Fig. 7. For solid ortho-D at 5 K, the upscattering is negligible for cold neutrons over the plotted energy range from 0.1 to 20 meV, whereas the downscattering amplitude becomes appreciable for incident neutron energy larger than 5 meV. Excluding the missing [011] Bragg peak ranging from 3 to 6 meV, the inclusion of contributions from inelastic scattering produces reasonable agreements to the total cross-section data derived from simple transmission measurements.

## V UCN cross-sections

The cross-sections for UCN production and UCN upscattering can be evaluated simply by finding the value of where , which is the dispersion curve of free neutron (the parabolic curve shown in Fig. 3). Here, the integration of Eq.(16) simplifies because the energy scale of UCN ( 350 neV) is several orders of magnitude smaller than the energy scale of cold neutrons (1 meV). For inelastic scattering events leading to UCN productions, the momentum transfer equals the initial momentum of the neutron due to the negligible . Similarly, the energy transfer . In estimating the UCN production cross-section, the integration of the double differential cross-section Eq.(1) becomes:

 σ(Ei→Eucn)=12∫dEucn∫dΩfkucnkiSdown(Q,+ω).

A change of integration variable , where , further simplifies the above integral to

 2πℏ22mnEi∫dω∫dQQSdown(Q,+ω).

For scattering with the final state of interest, i.e., UCN with energy , the range of integration is a narrow band around the parabola of the free neutron, , with , and thus the above integral can be reduced to:

 2πℏ22mnEiQ∗Sdown(Q∗,+ω∗)∫EiEi−Eucndω∫Q∗+kucnQ∗−kucndQ
 =2πℏ22mnEiQ∗Sdown(Q∗,+ω∗)Eucn2Δkucn, (18)

given that S(Q, ) is a smooth function. Next, to estimate the production of UCN that can be confined in the experimental apparatus, in which the Fermi-potential of the container material sets the maximum energy of UCN (), an additional step is required to integrate the UCN energy from 0 to :

 ΔkucnEucn = ∫kucn0dkucnℏ2k2ucn2mn = 13ℏ2k3ucn2mn=13kucnEucn.

As a result, the cross-section of UCN production for neutrons with incident energy is

 σ(Ei→(0−Eucn))=124πℏ22mnEikiSdown(ki,Ei)13kucnEucn, (19)

to produced storable UCN with from 0 to .

Two sets of experimental data on UCN production are compared to the predictions of the full model. First, note that the production cross-section reported in Atchison et al. (2007) is normalized to each atom and their model does not include either the spin statistics or the molecular form factor of molecular D; the cross-section independently reported in Muller (2008); Gutsmiedl et al. (2009) should be corrected by a factor of two to properly account for the range of integration on the allowed momentum transfer through to . To allow for direct comparisons, we multiply both data sets by a factor of two without any further corrections. Overall, the full model gives fair agreement with both the results of indirect UCN production extracted from the data of cold neutron scattering Muller (2008) and the results of direct UCN production Atchison et al. (2007)(as shown in Fig. 8). Around the peak production at 6 meV, the full model gives a striking agreement with the experimental data, whereas the IA approach using a somewhat realistic density of states misses the peak by 1 meV. The simple Debye model fails to capture any details of energy dependence. The transition at 7.1 meV is a delta function in the model, but the finite energy resolution of the scattering instrument integrating around the transition energy reduces the amplitude. In spite of it, the data suggests that the transition energy is slightly higher than 7.1 meV, as used in the calculation. The second peak around 9 meV predicted by the full model is absent in the Münich data Muller (2008), however, is present in the PSI data. The disappearance of the 9 meV peak could be a result of the preferred crystal orientation as discussed in Sec.IV. By restricting the scattering angle around the basal plane, the scattering intensity of the 9 meV peak can be adjusted to a smaller value (as shown in Fig. 6). For low scattering intensities at energy range below 4 meV and higher than 11 meV, the full model predicts production cross-sections smaller than the reported experimental data. This might simply due to the artifact of insufficient background subtraction in the experimental data. In particular, the subtraction of non-vanishing tails of the large elastic peaks is difficult to carry out. As for the high energy end, multiple scatterings couple to neutrons of energies higher than the Debye temperature, and thus enhance the scattering amplitudes. Neither of the data points is corrected for effects of multiple scattering. The very different dimension of the target used in these two experiments also hints that the origin of the high energy excess might be multiple scattering in origin.

UCN production can be optimized by varying the energy spectrum of the incident cold neutrons. Integrating the product of the energy-dependent cross-section for UCN production and the incident neutron flux over the energy spectrum yields a UCN production cross-section,

 σ(Tcn)=12∫dσ(Ei→Eucn)dEiEi(kBTcn)2e−Ei/kBTcndEi, (20)

that depends on the temperature of incident cold neutrons. Here we assume that the cold neutron flux has a thermalized energy spectrum with temperature . As shown in Fig. 9, the UCN production is the greatest when coupled to a flux with a 40 K Maxwell-Boltzmann energy spectrum, whereas the optimized cold neutron spectrum predicted by the IA is 33 K with a cross-section 32% larger than that predicted by the full model. This difference comes from the fact that the IA approach over-estimates the contributions below 5 meV, where the free neutron parabola does not intersect with most phonon branches at low (see Sec. IV). Even though the differential cross-section for UCN production peaks at 6 meV, the overlap with a flux of cold neutrons with a 70 K spectrum does not result in more UCN production, because many neutrons in the Boltzmann distribution are spread out over the higher energy range resulting in a reduction of peaked flux. Using a colder flux of neutrons, a higher percentage of neutrons are directly under the production peak.

The achievable UCN density does not only depend on the production cross-section, but also scales with the lifetime of UCN, which sets the limit on the time duration during which the UCN density can accumulate without loss inside the source. Upscattering is a source of loss that needs to be controlled. Following the same approach for UCN production, the cross-section for UCN upscattering is evaluated by finding the along the free neutron parabola as illustrated in the lower graph in Fig. 3. Note that the upscattering amplitude is directly proportional to the occupation number of phonons, which obey the Bose-Einstein statistics. The occupation number peaks at low and falls off in magnitude as increases. The amplitude of UCN upscattering depends even more strongly on the physical presence of these phonons than that of the downscattering. For small momentum transfers within the first Brillouin zone, the free neutrons do not intersect with any branches of acoustic phonon (see the lower graph in Fig.3). Without direct coupling to the acoustic phonons, there is no coherent process that upscatters UCN though phonon annihilation. For solids at low temperatures, the non-zero upscattering amplitude comes from the remaining incoherent 1-phonon and multi-phonon processes (see Fig. 10). As the temperature increases, the contribution from coherent scattering increases with increasing thermal population of phonons beyond the first Brillouin zone (shown in the inset plot in Fig. 10), and the coherent phonon annihilation take place. Fig. 10 presents the evolution of the differential cross-section of UCN upscattering with increasing temperature from 5 K to 18 K.

The total cross-section for UCN upscattering is calculated by integrating (lower graph in Fig.3) following:

 σ(Eucn→Ef)=12∫dEf∫dΩkfkucnSup(Q,−ω). (21)

Upon upscattering, UCN scatters into a narrow band of phase space with the final energy and momentum around the free neutron parabola (as illustrated in the lower graph in Fig. 3). The integration simplifies into:

 σ(Eucn→all possible Ef)= 12∫dϕ∫ω∗+EucnEucndω∫Q∗+kucnQ∗−kucndQQ∗1k2iSup(Q∗,−ω∗) =2π2∫∞0dEfQ∗1k2ucnSup(Q∗,−ω∗)2kucn =4π21kucn∫∞0dEfkfSup(kf,−Ef). (22)

In Fig. 11 plots upscattering cross-section for UCN with a wavelength of 500 . For UCN with a different wavelength, the upscattering is rescaled linearly with its initial wavelength. The full model predicts a UCN upscattering cross-section 24 times smaller than that estimated by IA (see Fig. 11) depending on the temperature. The excess in the IA prediction is due to the use of the density of states that over-estimates the number of modes in small energy transfers within the first Brillouin zone (see details shown in Fig. 10). This correction is by no means small. With a reduced upscattering cross-section, UCN can live longer inside the source at higher temperatures.

The ultimate limit on the UCN lifetime in solid D comes from the small, yet non-vanishing neutron absorption of each deuteron nuclide. In a source where the escape time of UCN is comparable to (or larger than) the absorption time, the UCN yield saturates when the upscattering loss is reduced to the same level as the nuclear absorption loss. Reduction of the upscattering by further cooling would not increase the UCN yield of UCN by more than a factor of two. In practice, the gain is smaller than two because of the additional sources of loss, such as the hydrogen contamination and the para-D contamination Liu et al. (2000); Morris et al. (2002). The condition for saturation in UCN output sets the practical operational temperature of the source, and determines the requirements on cryogenics engineering. The direct consequence of reduced upscattering cross-section in solid ortho-D is that the predicted saturation temperature is 7.5 K, which is 2 degrees higher than that predicted previously using the IA. The higher operating temperature makes the implementation of cryogenic source somewhat less challenging as most refrigerator systems have higher cooling power at elevated temperatures.

For a practical sample of converted D gas containing a residual 3% of para-D, the total loss cross-section is increased by 1 b, because the temperature-independent upscattering cross-section of pure para-D is 31 b Liu et al. (2000). In this case, the full model predicts the saturation temperature for UCN production to be 11 K, and the IA model predicts it to be 8 K. Using un-converted normal D gas with 33% para-D, the loss cross-section is about 10 b. According to the upscattering cross-sections presented in Fig. 11, the UCN yield would reach saturation at temperatures as high as 18 K, and no superthermal gain would be observed by cooling the UCN converter.

Finally, the updated upscattering cross-section predicted by the full model of for o-D is applied to understand the UCN production data from our recent measurements using solid D Lavelle et al. (2010). The Monte-Carlo simulation using the upscattering cross-sections calculated using the IA approach fails to reproduce the temperature dependence measured in the experiment. Details of the Monte-Carlo simulation can be found in Lavelle et al. (2010). We also investigated other effects that lead to increased elastic scattering inside the source, however, we find that the higher saturation temperature can only be explained by reduced loss. The coherent scattering leading to smaller upscattering cross-sections predicted using our full model is a very likely candidate to provide the required modification.

We then extract upscattering cross-sections using the same set of experimental data reported in Lavelle et al. (2010). The upscattering cross-sections of solid D are derived from another Monte-Carlo study, in which the upscattering cross-section is varied until the simulated result on UCN yield best fits the experimental data for each different temperature. Results of this study are shown in Fig. 12, compared to theoretical predictions based on the full model and the IA. The Monte-Carlo code includes the full range of UCN energy spectrum from 0 up to 1 eV. To facilitate direct comparison with previous calculations, the extracted upscattering cross-section are for UCN with a wavelength of 500 Å. Below the triple point, the upscattring cross-sections agree quite well with the full model prediction. Above the triple point, the mechanism for UCN upscattering (in the liquid phase) is significantly different from that in the solid, and none of the models discussed in this paper apply. Data points of small upscattering cross-sections fluctuates at low temperatures due to low statistics of UCN signal in general. At temperatures higher than 10 K, where the UCN upscattering cross-section is larger than 1 b, the incoherent model is excluded at the 2 level and higher.

We also surveyed experimental results on UCN production using solid deuterium reported independently by various research groups. Data sets are plotted in Fig. 13 together with the results of simulations. While the results of simulation using the updated upscattering cross-section agree quite well with our data and the data set measured at PSI, the data sets from LANL and Mainz groups show a much steeper temperature dependence for temperatures higher than 10 K. The difference comes from the different source configurations. In these two experiments, solid D was condensed from vapor at the end of the UCN guide, which was cooled below the solidification temperature of D. The source was designed to reduce the transmission loss by eliminating the vacuum window that is typically installed to contain the volatile D. The windowless source worked quite well at low temperatures, however, for temperatures higher than 10 K, the entire UCN guide was filled with D gas at the saturated vapor pressure. The large upscattering cross-section from the D vapor, together with a large volume where the vapor permeated, resulted in a steeper temperature dependence than the simple prediction where UCN upscatters through absorption of phonon in the confined region of the source. Both PSI and our experiment used D contained in a target cell, in which the vapor region is considerably smaller. None of the data sets show a low saturation temperature, as predicted by the IA. Instead, the temperature dependence of experimental data points agrees better with results of the simulation using upscattering cross-sections predicted by the full model.

## Vi Conclusion

In this paper, we have presented calculations of the neutron cross-sections in polycrystalline ortho-D. All known physics on spin statistics, rotational excitations, elastic Bragg scattering, coherent inelastic, and incoherent inelastic scattering have been implemented in our full model for the first time. The presented work addresses the shortcomings of the widely used IA in estimating neutron cross-sections at energies at thermal and cold regime. For applications on modern neutron sources for UCN production using solid D and other materials with large coherent scattering lengths, such as helium, oxygen and nitrogen, the interference between scattered neutron wavefunctions across all lattice sites is important and should not be ignored. For the case of solid ortho-D, we have shown that UCN cross-sections are significantly altered under the new treatment. In particular, the UCN upscattering cross-section is a factor of 24 smaller than previous predictions throughout the temperature range of interest to UCN production. The same modifications will impact designs of sources (beyond UCN sources) that aim to produce long wavelength neutrons using materials other than hydrogen. The Monte-Carlo algorithm developed in this work can be applied to all materials with significant coherent scattering components, provided that the energy dispersion of the quasi-particles is known.

###### Acknowledgements.
The work was supported by NSF grants 0457219, 0758018.

### References

1. S. Glasstone and M. Edlund, The Elements of NUCLEAR REACTOR THEORY (D. Van Nostrand Company, Inc., New York, 1952).
2. L. A. N. Labratory, Mcnp5 monte carlo code, ccc-730 (2008).
3. R. MacFarlane, New thermal neutron scattering files for endf/b-vi release 2 (1994).
4. J. A. Young and J. U. Koppel, Physical Review a-General Physics 135, A603 (1964).
5. J. Schwinger, Physical Review 58, 1004 (1940).
6. M. Hamermesh and J. Schwinger, Physical Review 69, 145 (1946).
7. Njoy – nuclear data processing system (1999).
8. C. Y. Liu, A. R. Young, and S. K. Lamoreaux, Physical Review B 62, R3581 (2000).
9. F. Atchison, B. Blau, K. Bodek, B. van den Brandt, T. Brys, M. Daum, P. Fierlinger, A. Frei, P. Geltenbort, P. Hautle, et al., Physical Review Letters 99, 262502 (2007).
10. J. R. Granada, EPL 86, 66007 (2009).
11. P. A. Egelstaff, B. C. Haywood, and F. J. Webb, Proceedings of the Physical Society of London 90, 681 (1967).
12. R. J. Elliott and W. M. Hartmann, Proceedings of the Physical Society of London 90, 671 (1967).
13. H. W. Diehl and W. Biem, Zeitschrift Fur Physik B-Condensed Matter 20, 137 (1975).
14. V. V. Danchuk, N. N. Galtsov, M. A. Strzhemechny, and A. I. Prokhvatilov, Low Temperature Physics 30, 118 (2004).
15. G. L. Squires and A. T. Stewart, Proceedings of the Royal Society of London Series a-Mathematical and Physical Sciences 230, 19 (1955).
16. W. Schott, Zeitschrift Fur Physik 231, 243 (1970).
17. M. Nielsen and H. B. MÃ¸ller, Physical Review B 3, 4383 (1971).
18. M. Nielsen, Physical Review B 7, 1626 (1973).
19. J. W. Schmidt, M. Nielsen, and W. B. Daniels, Physical Review B 30, 6308 (1984).
20. R. E. Dewames, T. Wolfram, and G. W. Lehman, Physical Review 138, A717 (1965).
21. G. W. Lehman, R. E. Dewames, and T. Wolfram, Physical Review 128, 1593 (1962).
22. M. F. Collins, Proceedings of the Physical Society of London 80, 362 (1962).
23. J. C. G. Houmann, Physical Review B 1, 3943 (1970).
24. V. G. Vaks and K. Y. Khromov, Journal of Experimental and Theoretical Physics 106, 495 (2008).
25. S. Lovesey, Theory of Neutron Scattering from Condensed Matter, vol. 1 (Oxford University Press, Oxford, 1984).
26. D. Colognesi, F. Formisano, A. J. Ramirez-Cuesta, and L. Ulivi, Physical Review B 79, (2009).
27. P. A. Egelstaff and R. S. Pease, Journal of Scientific Instruments 31, 207 (1954).
28. Bostanjo.O and Kleinsch.R, Journal of Chemical Physics 46, 2004 (1967).
29. G. W. Collins, W. G. Unites, E. R. Mapoles, and T. P. Bernat, Physical Review B 53, 102 (1996).
30. H. Stein, H. Stiller, and Stockmey.R, Journal of Chemical Physics 57, 1726 (1972).
31. W.-D. Seiffert, Tech. Rep., Euroaäische atomgemeinschaft euratom (1970).
32. C. Lavelle, W. Fox, G. Manus, P. McChesney, D. Salvat, Y. Shin, M. Makela, C. Morris, A. Saunders, A. Couture, et al., Ultracold neutron production in a pulsed neutron beam line (2010).
33. C. Y. Liu and A. Young, arxiv:nucl-th/0406004 (2004), submitted to Phys. Rev. B.
34. Z. C. Yu, S. S. Malik, and R. Golub, Zeitschrift Fur Physik B-Condensed Matter 62, 137 (1986).
35. L. J. Raubenheimer and G. Gilat, Physical Review 157, 586 (1967).
36. E. Gutsmiedl, A. Frei, A. R. MÃ¼ller, S. Paul, M. Urban, H. Schober, C. Morkel, and T. Unruh, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 611, 256 (2009).
37. A. Muller, Ph.D. thesis (2008).
38. L. Thorsten, Ph.D. thesis (2010).
39. C. L. Morris, J. M. Anaya, T. J. Bowles, B. W. Filippone, P. Geltenbort, R. E. Hill, M. Hino, S. Hoedl, G. E. Hogan, T. M. Ito, et al., Physical Review Letters 89, 272501 (2002).
40. A. Serebrov, V. Mityukhlyaev, A. Zakharov, A. Kharitonov, V. Shustov, V. Kuz’minov, M. Lasakov, R. Tal’daev, A. Aldushchenkov, V. Varlamov, et al., Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 440, 658 (2000).
41. F. Atchison, B. Blau, B. van den Brandt, T. BryÃâº, M. Daum, P. Fierlinger, P. Hautle, R. Henneck, S. Heule, K. Kirch, et al., Physical Review Letters 95, 182502 (2005).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters