A super-resolution technique to analyze single-crystal inelastic neutron scattering measurements using direct-geometry chopper spectrometers

Inelastic neutron scattering (INS) is a powerful probe of fundamental excitations in solids, including those of vibrational or magnetic origins. Neutron scattering from these excitations is characterized by the four-dimensional (4D) scattering function, S(Q, E), where Q is the momentum transfer and E is the energy transfer. In recent years, the use of highly pixelated detector arrays in conjunction with direct-geometry chopper spectrometer (DGS) instruments^1–7 has allowed for efficient measurements of single crystal 4D S(Q, E) functions over large ranges of Q and E. One single measurement of a single-crystal sample at a DGS instrument captures the sample scattering function S(Q, E) on a three-dimensional (3D) manifold in the 4D Q, E space. Typically, a sample is rotated around a single axis, while maintaining a single wavelength (i.e., monochromatic energy) of incident neutrons, to scan the 4D Q, E volume. A series of 3D manifolds measured during the scan are then combined into the volumetric 4D dataset, allowing for extracting 2D slices at high-symmetry Q directions by using software packages, such as Mantid,⁸ Horace,⁹ DAVE,¹⁰ Mslice,¹¹ and Utsusemi.¹² This is done by integrating the measured scattering intensity along two of the four dimensions, yielding a scattering intensity as a function of the remaining two dimensions. Most often, this is done to extract dispersion relations of fundamental excitations within crystalline solids by illustrating measured scattering intensity as a function of energy transfer along the vertical axis and a single direction in wave-vector transfer along the horizontal axis. These extracted dispersion data can be used to directly compare to model calculations of the materials dynamics in order to constrain or refine parameters in a model (i.e., the Hamiltonian). For magnetic systems, SpinW¹³ is often used to fit spin-wave models to the experimentally obtained dispersion along multiple high-symmetry directions simultaneously. This allows one to obtain exchange parameters or other energy dependent terms in a model Hamiltonian, which can represent the spin dynamics of the system.

Quantifying a crystalline material’s vibrational or magnetic dispersion is key to understanding the system’s dynamics. The accuracy of the measured dispersion (and that of the inferred quantities, which define the dynamics, such as exchange couplings or force constants) obtained from DGS experiments is bound by the instrument resolution. However, a significant number of measurements do not take into account instrumental resolution effects beyond the use of an analytical approximation of a single energy and wave-vector resolution of the instrument. DGS instrument resolution for the S(Q, E) scattering function is four-dimensional, and the effects of the resolution function are compounded by the slope and the curvature of the dispersion surface, making it cumbersome to accurately model. Given a set of instrument and experimental parameters, such as chopper settings and the sample shape, the resolution function for DGS neutron scattering instruments still varies at different (Q, E) points in the measured dynamical range. This resolution ellipsoid varies considerably across the detector array of the instrument and energy transfer. Furthermore, the resolution ellipsoid can have significant tilt that can lead to focusing and defocusing effects in the measured dispersion depending on the slope of the dispersion relative to the tilt. This can make it difficult to accurately extract a model based on the measured dispersions.

Recently, it was demonstrated¹⁴ that some super-resolution imagery techniques can be adapted to improve energy resolution in a phonon density of states measurement, g(E), based on neutron scattering techniques. This analysis was performed upon a measurement of the phonon density of states as a function of a single independent variable, energy transfer, which is a routine measurement for direct geometry chopper spectrometers measuring powder samples. In this work, we introduce a super-resolution technique to improve the ability of extracting dispersion from inelastic neutron scattering measurements of single crystal samples. We use techniques inspired by image correlation methodologies to improve extraction of information from measurements as a function of two independent variables, energy and wave-vector transfer.

Dispersion relations are obtained from single crystal measurements at DGS instruments by extracting slices of scattering intensity as a function of energy transfer, E, and wave-vector transfer, Q, along high-symmetry directions in reciprocal space. These slices are obtained by integrating a portion of the measured reciprocal space along the two orthogonal reciprocal space directions to the direction of interest in the chosen slice. Constant Q cuts through these 2D slices can then be extracted by integrating a portion of the slices along wave-vector transfer and plotting the resulting scattering intensity as a function of energy transfer. For each constant Q cut through a single 2D slice, centers of peaks in the energy spectrum are recorded as the excitation energies for this particular Q value. The energy of these peak centers plotted as a function of Q is the measured dispersions. An analytic or numerical model of the dispersion can then be directly compared to the experimentally determined dispersion values in order to determine model parameters. This is often done using non-linear curve fitting algorithms. This methodology of extracting dispersion parameters based on peak location does not account for instrumental energy or wave-vector resolution. Such quantities may shift or skew peak locations and, therefore, would serve to skew any model parameters determined from a comparison of peak positions to model dispersion. The method we present here calculates resolution-based corrections to the measured dispersion functions.

Figure 1 shows the steps of this dispersion correction workflow. An experimental slice is first obtained from data reduction [Fig. 1(a)]. A dispersion dataset is then obtained directly from the experimental slice, E_exp;0(q), following the normal steps outlined earlier. This quantity is shown as open symbols in Fig. 1(b). Those excitations are then fit to a dispersion model (in this case, the spinwave model generated in SpinW), and relevant parameters (in this case, exchange parameters) are obtained [Figs. 1(a)–1(d)]. The scattering cross section as a function of momentum transfer and energy transfer, including both the dispersion relation and the scattering amplitude, is now determined by the excitation model and relevant parameters. The next step is to obtain the instrument resolution function [Fig. 1(e)]. This can be achieved by using analytical calculations (this is often done for triple axis spectrometers^15–18) or Monte Carlo neutron ray tracing simulations. More details can be found in Sec. II B. By convolving the calculated scattering function with the instrument resolution function, we can obtain the modeled slice [Fig. 1(f)]. Then, we can compare the modeled slice and the experimental slice and obtain the energy shifts, ΔE(q) [Fig. 1(g)], required to match the modeled slice to the experimental slice. This step is key to the super-resolution dispersion determination and is further explained in Sec. II C. The energy shifts obtained can then be applied to the model dispersion curve (from first fit to the experimental data) to obtain the corrected dispersion curve [Fig. 1(h)],

The focus of this work is to demonstrate that measured dispersion relations can be corrected taking into account resolution effects, but a natural further step is to fit the dispersion model to the corrected dispersion, E_exp;1(q), to obtain the corrected parameters for the Hamiltonian.

The 4D resolution function for a DGS instrument in measurements of single-crystal samples has been modeled analytically^19,20 using the covariance matrix to simplify the treatment. It has also been calculated using Monte Carlo ray-tracing simulations^21–25 and was sometimes approximated using Gaussian functions when it was used in resolution convolution.²⁶ However, these corrections are not often taken into account in 4D DGS data analysis due to the complexity of the resolution convolution and the computing resources required for fitting the dispersion model obtained from a highly pixelated detector array.

Incorporating resolution into the single crystal dispersion fitting workflow is difficult because of the following:

• The resolution is 4D in nature, thus requiring a four-dimensional integration to convolve it with a model S(Q, E).

• The resolution function varies across the dynamical range depending on the momentum and energy transfer. For example, as the energy transfer increases, the energy resolution broadening decreases for direct geometry chopper spectrometers.

• The resolution ellipsoid can have a significant tilt that can lead to focusing and defocusing effects in the measured dispersion depending on the slope of the dispersion relative to the tilt of the resolution ellipsoid.

• For some DGS instruments, the energy resolution function is asymmetric as a result of the moderation process peculiar to neutron production in spallation neutron sources.²⁷ 

The last point here has very rarely been taken into account previously in single crystal dispersion fitting. Furthermore, to fit the dispersion model, incorporating the instrument resolution requires multiple evaluations of the dispersion model with varying parameters and resolution convolution with the model for each set of parameters. The number of iterations depends on the optimization algorithm used and how close the initial guess was to the optimal model parameters. Such an optimization procedure can be demanding in both programming and computing resources.

The super-resolution procedure outlined in this work is agnostic to the technique of resolution calculation and convolution. In this work, we have chosen to use a technique based on the Monte Carlo ray tracing simulation to illustrate the super-resolution methodology.

The MCViNE package^22,28 has been used to compute the resolution function for single crystal measurements at DGS instruments.^22,24,25 The procedure is reused here to simulate the resolution function (or point-spread function). The simulation starts with a beam simulation that matches experimental conditions, such as incident energy and chopper settings. In order to calculate the energy and wave-vector resolution function with MCViNE, we use a virtual sample that has the same geometric shape and lattice parameters of the real sample to scatter neutrons in the vicinity of a particular set of momentum and energy transfer h, k, l, E. In the simulation, we also take advantage of the measured UB matrix to orient the virtual sample just like what has been measured during the experiment. The virtual sample is also rotated around the vertical axis to a particular ω angle, similar to what happens in real measurements. The SEQUOIA detector system is simulated according to its specification, such as the positions and orientations of all detector packs, the ³He tube radius, length, and its spacing in the detector pack, the pressure of ³He gas in the detector tubes, and the detector pixel height. Only events that arrive at the particular detector pixel and time-of-flight bin corresponding to the nominal h, k, l, E are collected. These detector events are then reduced to $h̃,k̃,l̃$ and $Ẽ$⁠. Those $h̃,k̃,l̃,Ẽ$ values center around the nominal h, k, l, E, as expected. The differences between those $h̃,k̃,l̃,Ẽ$ from the nominal h, k, l, E are kept in a list of dh, dk, dl, dE. The Monte Carlo ray tracing approach captures details of the 4D resolution function, including, for example, the asymmetrical energy dependent line-shape mentioned previously.

To decrease the complexity of the resolution modeling and convolution, we perform them in two dimensions along the two axes of a slice. The details of the convolution will be presented later in this article. We calculate a list of dq from the saved list of dh, dk, dl, where q is along the high-symmetry Q direction of the slice of interest. The dq, dE events are then histogrammed into a profile of the point spread function (PSF) for the 2D slice. This procedure is repeated for the points on a grid on the (q, E) plane. Each PSF at one grid point is fit to a sheared 2D function, with one axis that is “energy-like” and has an asymmetric shape and another axis that is “momentum-like” and modeled as a Gaussian. An example of this parameterization is presented in Fig. 2. Then, an interpolation of the fitted parameters allows us to calculate the PSF function at any point in the q, E space.

Before convolution, the scattering intensities of the dispersion model are first integrated along the two Q directions perpendicular to the q direction for the slice of interest. Then, the integrated data are convolved with the 2D resolution function modeled earlier to obtain the convoluted slice. This convolution method is a good approximation of the full 4D convolution, and an example is shown in Sec. III B to demonstrate that.

The basic idea here is to compare the two slices, namely, the experimental slice [Fig. 1(a)] and the resolution-convolved modeled slice [Fig. 1(f)], to find the displacements between the dispersions in the two images, which can then be used to correct the model. Finding displacements (or disparity) in two images has been a long-standing challenge in image processing. Stereo imaging techniques that uncover depth information of a scene captured in two images by finding the disparity field between the images have been reviewed multiple times during the last few decades.^29–31 The techniques in traditional stereo vision find the disparity field d(x, y) that minimizes the difference between the left image (image taken by using a camera on the left) and the right image (image taken by using a camera on the right) warped by the disparity,

Used in this formula is the sum of the absolute differences, but often, the sum of square differences or normalized cross correlation is also used. I_l(x, y) and I_r(x, y) are left and right image intensities, and d(x, y) is the displacement along y. Local methods for dense disparity calculation minimize cost functions [Eq. (2)] for local patches, while global methods^32,33 and semi-global methods³⁴ take the smoothness of disparity into account by constraining the disparity d(x, y) using regularization. Advancements in computing techniques for disparity calculation have been proven useful in many fields of quantitative research, including remote-sensing and geophysics.³⁵ We aim to reuse the disparity calculation technique to estimate the displacement field between the experimental slice and the resolution-convolved model slice. In this first demonstrative work for application of image disparity calculation in the data analysis of neutron scattering measurements, we simplify the problem by limiting the slice to only one visible dispersion; therefore, the displacement field can be simplified to be independent of E, and the minimization problem becomes

The simplification of limiting the displacement field to be independent of E makes this optimization problem straightforward to program by using a SciPy³⁶ optimizer, for example.

In this work, we illustrate the super-resolution dispersion technique using a synthetic dataset that resembles a real experimental dataset, for which the full analysis will be reported elsewhere.³⁷ In the experiment, a Mn₃Si₂Te₆ single-crystal sample was measured at the SEQUOIA instrument³ with incident energy E_i = 60 meV in the high flux mode, and the sample was rotated at least 180° about the vertical axis, 1° per step, to cover a large volume in the reciprocal space. Slices along multiple high-symmetry directions show clear dispersions. The dispersions along those high-symmetry directions were obtained by finding centers of peaks in constant q cuts. They were fit to a Hamiltonian without on-site anisotropy, which allows for anisotropic interactions for each of the exchange interactions to account for the spin–orbit coupling,³⁷ 

Here, {J_i} are exchange coefficients, as shown in Fig. 3, while {Δ_i} introduce anisotropy.

First, we build a synthetic dataset for which we know the exact model and parameters for the dispersion surface. The synthetic dataset is obtained from a virtual neutron experiment performed by using the MCViNE software. MCViNE contains a scattering kernel that scatters neutrons according to a dispersion surface that a user can define by using arbitrary analytical functions.^22,23,28 Therefore, an analytical dispersion function similar to the dispersion surface of the spinwave model defined by Eq. (4) was employed. The spin coupling constants used in the spinwave model in Eq. (4) are J₁ = 1.663 meV, J₂ = 0.477 meV, J₃ = 0.835 meV, Δ₁ = 0.390, Δ₂ = −0.554, andΔ₃ = −0.431. Analytical functions were used to approximate the dispersion and scattering intensity in the vicinity of the (002) wave-vector, and they were parameterized as

where E_a = 11.6 meV, E_b = 9.05 meV, S0 = 14.8, Γ_h = Γ_l = 0.38 r.l.u., and Γ_k = 0.35 r.l.u.

Then, we performed a virtual experiment of a sample with this analytical dispersion function and scattering intensity using MCViNE. The simulation and data reduction consist of the following steps:

• Beam simulation. We can reuse the beam simulation performed earlier for the resolution calculation as it contains all the information regarding the instrument.

• Sample scattering simulation. The sample has the same geometrical shape as the real sample, with the dispersion defined in Eq. (6). Simulations with a series of ω rotation angles are performed, matching the real experiment.

• Detector simulation. The scattered neutrons are intercepted by the virtual SEQUOIA detector system, and the neutron events detected are saved in NeXus files, one for each ω angle.

• Reduction. The Mantid software⁸ is then used to reduce the NeXus data files in the same way as the real experiment, and corresponding slices are made.

Once we have the virtual experimental data, we perform the super-resolution workflow outlined in Fig. 1 upon it. The results are shown in Fig. 4. The corrected dispersion shows clear improvement in agreement with the model dispersion curve, compared to the original dispersion obtained directly from the experimental slice. The root-mean-square (rms) difference between the experimental dispersion data and the model dispersion data was reduced from 1.59 to 0.32 meV, showing a nearly fivefold improvement. It is worth noting that the nominal resolution of the SEQUOIA instrument at E_i = 60 meV for the high-flux mode is ∼2 meV. The fact that one single iteration of our super-resolution workflow can improve the dispersion data accuracy by five fold means our method is an efficient way to optimize the dispersion model while taking into account the instrument resolution effect.

Another way to check the improvement of the dispersion data is to observe the improvement of the fitting parameters in Eq. (6). Table I presents the model parameters for the exact model, the fitted model without resolution correction, and the fitted model with resolution correction. Without the correction, we obtained E_a = 12.93 meV and E_b = 9.53 meV while fitting the dispersion data to Eq. (6). Compared to the exact values of E_a = 11.6 meV and E_b = 9.05 meV, the fitted values are off by 1.33 and 0.48 meV, respectively. After correction, the fitting results are E_a = 11.80 and E_b = 8.78, and the errors are reduced to 0.20 and 0.27 meV.

An intuitive illustration of the improvement of the quality of the dispersion model is also presented in Fig. 5. Here, Fig. 5(a) shows the 00L slice obtained from the virtual experiment data. Figure 5(b) shows the resolution-convolved slice obtained from the dispersion model fitted to the original dispersion data from the virtual experiment without correction. It is clear that the dispersion is shifted upward in comparison to the virtual experimental data in Fig. 5(a). Figure 5(c) shows the resolution-convolved slice obtained from the dispersion model fitted to the corrected dispersion data, and this slice agrees much better with the experimental slice in Fig. 5(a). This agreement is also a validation of our resolution convolution procedure. The better agreement of the corrected dispersion model with the virtual experimental data is also evident in the residual plots shown in Figs. 5(d) and 5(e).

We also applied the super-resolution dispersion workflow on the experimental SEQUOIA dataset described in Sec. III A. In Sec. III B, for the virtual experimental data, we apply the super-resolution procedure to one single slice. In comparison, for the real experimental data, we treat multiple slices along different Q directions. The first step of the workflow remains the same; the dispersion data for each slice is first obtained without considering the resolution effect. Then, these dispersion data from multiple slices are fitted to the spin-wave model simultaneously to obtain the original set of the model parameters. The fitted model provides dispersion curves and scattering amplitudes for every slice, and they are convolved with resolution functions to obtained modeled slices. Each modeled slice is then compared to the corresponding experimental slice to obtain disparity curves, which are used to correct the dispersions. The corrected dispersions of all interested slices are then fit to the spin-wave model simultaneously again, yielding a new set of model parameters. One example of the dispersion correction is presented in Fig. 6 for the 00L dispersion. Energy corrections in the order of 1 meV are found near 002. The full dispersion correction data are reported elsewhere.³⁷ 

A new technique to obtain super-resolution dispersions along high-symmetry Q directions for single crystal measurements employing direct geometry neutron spectrometers is developed. This is done by computing the disparity curve between the resolution-convolved-model slice and the experimental slice and then applying the disparity to correct the dispersion. Here, the resolution-convolved slice was obtained by convolving the resolution with the scattering intensity of a dispersion model that was fit to the experimental dispersions obtained without any consideration of instrument resolution. The disparity of the slices was obtained by minimizing the difference between the experimental slice and the modeled slice warped by disparity, subjecting to the total variation regularization. The technique clearly shows improvements in the determination of dispersions, and it is computationally faster since classical methods would have required multiple iterations of model evaluation and resolution convolution with many more sets of model parameters. We show that this method can achieve fivefold super-resolution with respect to nominal resolution of the SEQUOIA³ instrument. The demonstration is facilitated by a MCViNE-based virtual experiment, which provides the virtual experimental data and the known target model to check the effectiveness of the super-resolution technique.

This super-resolution dispersion technique is limited by the signal-to-noise ratio as other imaging techniques. The 2D resolution convolution method used in this work can be updated to use 4D resolution convolution to improve the accuracy and the universality of this approach. More sophisticated disparity computation techniques can be adapted to remove the limit of single dispersion per slice. Finally, many image processing techniques may find various applications in neutron data analysis. Another potential application of the disparity calculation technique is to detect super-resolution variations in dispersions with respect to temperature/pressure under which the sample is measured.

This work focuses on the correction of the dispersion relation E(q) between the excitation energy E and its momentum q. It can be envisioned that the super-resolution techniques developed in the previous work for powder DGS data¹⁴ can be extended and combined with techniques developed in this work to reconstruct super-resolution I(q, E) slices, reducing the influence of the instrument broadening and providing information on the intrinsic linewidths of the excitations, Γ_E(q).

The authors thank Hillary Smith, Brent Fultz, Garrett Granroth, and Doug Abernathy for fruitful discussions. This work was partially supported by the Department of Energy, Laboratory Directed Research and Development SEED funding, under Contract No. DE-AC05-00OR22725. Work at the Spallation Neutron Source at Oak Ridge National Laboratory (ORNL) was supported by the Scientific User Facilities Division, Office of Basic Energy Sciences, U.S. Department of Energy (DOE). This research also used resources of the Spallation Neutron Source Second Target Station Project at ORNL. ORNL is managed by UT-Battelle LLC for DOE’s Office of Science, the single largest supporter of basic research in the physical sciences in the United States.

All authors declare that they have no conflicts of interest.

The data that support the findings of this study are available from the corresponding author upon reasonable request.