Super-resolution energy spectra from neutron direct-geometry spectrometers

Inelastic neutron spectroscopy (INS) is a powerful probe of fundamental excitations in solids, including those of vibrational (phonon) or magnetic origins. Phonon measurements play a key role in understanding the physical properties of a solid, such as heat transportation and electrical conductivity,^1–4 superconductivity,^5–8 as well as thermodynamic quantities such as entropy, which influences the phase stability of materials.^9–13 Measurement of changes in phonon frequencies following changes in thermodynamic parameters such as temperature,^11,14–16 pressure,^17,18 and chemical composition⁸ provides important data for understanding the origin of anomalies in phonon behavior, such as anharmonicity. Measuring the phonon density of states (DOS), a reduced representation of the vibrational property of condensed matter, is usually the first step in determining the phonon properties of a material experimentally (see, e.g., Refs. 11 and 16).

Neutron direct geometry spectrometers (DGSs) are important and convenient tools for measuring phonon DOS of a powder sample. However, resolution functions in DGS measurements of phonon DOS are well known to be cumbersome to model.¹⁹ The resolution function of the ARCS instrument,²⁰ a typical DGS instrument, is not a simple Gaussian function, as usually assumed (justifiably) in most optical imaging techniques. For it and all single chopper DGS instruments at a spallation source, the resolution function is asymmetric. Such asymmetry arises from the moderation process peculiar to neutron production in spallation neutron sources.²¹ This complexity is neglected by most studies so far.

Another complexity associated with the resolution function of DGS instruments is it varies as a function of neutron energy transfer (E). As E increases, the energy resolution broadening decreases [Fig. 1(c)]. When the phonon energy range is not large, resolution variation is small enough so that studies can be done by comparing and/or fitting models (for example, Born–von Kármán models) and density functional theory (DFT) calculations to data, involving a convolution with simplified resolution functions,^6,11,14,22 ignoring finer features in data blurred by resolution broadening and statistical noise. However, for a material with a wide range of phonon energies, the resolution of a DGS instrument at lower energy transfers is unsatisfactory for higher incident energy [Fig. 1(c)]. One way to alleviate this problem is to measure the phonon spectrum using multiple incident energies. The highest E_i dataset has the largest dynamic range, r · E_i < E < E_i (r is a number that is approximately the relative resolution of the instrument at the elastic line, where there is no energy transfer; for example, for ARCS at SNS,²³ r ≈ 3%–5%), and covers almost the full phonon DOS spectrum, but the resolution of the lower energy transfers is not optimal. By using a lower E_i, the phonon DOS is measured over a smaller dynamic range but with finer energy resolution. All DOSs are stitched together to obtain one DOS curve. This is done by replacing the low E portion of the DOS from the higher E_i measurement with the partial DOS obtained from the lower E_i measurement. The resolution of the final stitched phonon DOS, however, is still uneven across the energy range.

Deconvolution techniques have been demonstrated for a few cases in neutron spectra^24,25 but have not been widely applied. The usual objection to deconvolution is that it “generates information from void.” As we will show below, however, the information content of the DGS data is actually not fully utilized. This is particularly true for the modern DGS instruments on spallation sources with finely pixelated cylindrical detector banks such as Merlin, ARCS, SEQUOIA, CNCS, LET, 4-SEASONS, and HRC.^20,26–31 Furthermore, advances in modeling the instrumental resolution function and in reconstruction have only become available recently. This contribution will show how these advances can use the full information content.

Specifically, we demonstrate a super-resolution technique for spectroscopy that builds upon techniques used in super-resolution imagery. These techniques work when (1) the sampling frequency in the measurement is higher than that of the nominal resolution; (2) sub-bin shifts (similar to subpixel displacements in multiframe imaging super-resolution) exist for multiple measurements (for a DGS instrument, one measurement per detector pixel) of the same data; and (3) the point spread function (PSF) has a sharp feature (contains high frequency components). When these conditions are met, the spectra from multiple measurements can be fused together, and an iterative reconstruction can be used to obtain data in finer resolution.

Super-resolution (SR) techniques and statistical methods have seen wide application in optical imaging techniques such as surveillance, forensic science, remote-sensing, and microscopy.^32–38 A typical multiframe imaging SR procedure³⁸ consists of an image registration step³⁹ (or “fusion”), followed by deconvolution of the point spread function. The fusion step essentially improves the sampling frequency of the “fused” image by combining multiple images of lower resolution shifted by subpixel displacements, as often seen in, for example, video footage of moving objects. This important step was not examined in the previous deconvolution work.^24,25 In direct-geometry spectrometers [see the schematic in Fig. 1(a)], the spectrum is measured in terms of energy transfer E, obtained by conversion from time-of-flight data: $E=Ei−Ef=Ei−12m(Lspt−ti)2$⁠, where E_i, E_f are initial and final energies of a neutron, L_sp is the sample to detector pixel distance, t is the total time of flight (TOF) measured, and t_i is the time-of-flight from the neutron moderator to the sample. For direct geometry spectrometers, E_i and t_i are fixed for an experiment. The time-of-flight at all detector pixels is clocked synchronously to a 0.1 μs precision. Therefore, as events are accumulated into TOF bins during reduction, the minimum bin size is 0.1 μs at SNS. As detector pixels are at varying distances from the sample [Fig. 1(a)], the constant time binning corresponds to an energy binning that varies across detector pixels [Fig. 1(b)]. The typical “event-mode” DGS reduction procedure accumulates events in one array of energy bins. Therefore, the data are fused together similar to multiframe SR imagery [Fig. 1(b)]. One obvious difference from multiframe SR imagery is that for the current work, the data are 1D (energy axis only). More nontrivial differences exist. First, the conversion from TOF to energy transfer E is not a linear transformation (e.g., affine transformation) as often used in optical SR image registration. Furthermore, in contrast to SR imagery, the initial bins (TOF) are usually very fine. When they are transformed to energy bins, they are usually much smaller than nominal energy resolution. The conventional wisdom in the DGS data reduction procedure is to use an E bin size that is a fraction of the nominal resolution of the instrument. For example, for ARCS, an E bin size (ΔE) at 1% of the incident energy is often used, while the nominal resolution is 3%–5%. For E_i = 300 meV, ΔE = 3 meV. Therefore, an E bin would typically encompass many TOF bins. At E = 150 meV, ΔE = 3 meV corresponds to ∼57 TOF bins. However, if we ever want to improve resolution, a higher sampling rate for E is needed. ΔE = 0.25 meV would correspond to ∼4 TOF bins. As the transformation is nonlinear, one energy bin almost always corresponds to a noninteger number of TOF bins. For example, as shown in Fig. 1(b), for E near 110 meV, one 0.25 meV energy bin corresponds to 3 or 4 TOF bins, depending on pixel; near 150 meV, one energy bin corresponds to 4 or 5 TOF bins; and near 70 meV, one energy bin corresponds to 2 or 3 TOF bins. The current event-mode reduction assumes that the events happen at exact time-of-flights. As events can happen near the boundaries of energy bins, there is always an error accumulating counts into energy bins for a single detector pixel, even if we perform a measurement for infinite time. However, because the variations in sample-pixel distances result in variations of energy bins that have “sub-E-bin” components, this kind of error is minimized as data from all pixels are combined, and we can safely use an E bin size similar to those corresponding to the 0.1 µs TOF bin. This benefit is enhanced by the cylindrical arrangement of the detector tubes in the detector array²⁷ and would have been severely limited if the traditional spherical configuration was used. With this realization, our first modification of the data reduction was to increase the sampling frequency in the energy axis to at least 3 times the usual binning frequency or ∼15 times finer than the nominal resolution. The high sampling rate in E provides a solid foundation for the next processing steps in achieving super-resolution.

Our next step is to find an accurate model of the “point spread function” (PSF) for a DGS energy spectrum. Previously, C₄H₂I₂S samples were used to measure the inelastic resolution functions experimentally at several energy transfers.²⁰ However, as such calibrants rely on the molecular structure, they are not tunable and thus modeling must be used to obtain the energy-dependent resolution function at an arbitrary energy transfer. Figure 1(a) shows a simple schematic of a DGS instrument. One major cause of resolution broadening in energy is the line shape of neutron spectra emitted from the moderator, which is asymmetric. As the line shape is a function of time and energy (velocity), filtering by the Fermi chopper selects a narrow band of time and produces a symmetric distribution. As the neutrons of different speed disperse while traveling to the detector, the shape becomes asymmetric again but with the tail now at short times. So, the chopper acts like a pin-hole in time of flight.¹⁹ 

More factors influence the instrument broadening, including the divergence of the neutron beam formed through neutron guides, the sample shape and size, and the detector geometry. Monte Carlo neutron ray tracing simulations with the MCViNE package^40,41 can capture these details and are used to calculate the PSF functions by taking relevant factors into account: modeling the neutron beam, the sample with a δ-function scattering kernel with a geometric shape and a scattering cross section same as the real sample, and the detector system in high accuracy.^20,41–43 In Fig. 1(c), examples of MCViNE-simulated, energy-dependent energy resolution functions are presented.

In order to obtain resolution functions for energy transfers across the dynamical range of the measurement, it is possible to use MCViNE to simulate at each energy bin a point spread function, but such computations are demanding in computing resources. We chose to first compute the resolution function using MCViNE at a few energy transfers and then fit each of these resolution functions to a revised⁴⁴ Ikeda-Carpenter function²¹ with parameters a, b, σ, and t₀, depending on energy transfer [Fig. 1(d)]. These parameters were found to vary slowly across the energy range of interest and can be interpolated to obtain a PSF at any point along the energy axis. One thing to remember is since the phonon DOS may be obtained from experiments measured using multiple incident energies, this procedure of simulation, modeling, and interpolation is required for all incident energies.

An example of the model fitting can be seen in Fig. 1(c). There, the crosses are MCViNE simulation results of PSF functions, and the solid lines are from fitted models, which agree very well with the simulation data points. An example of the fitted parameters as functions of energy transfer is presented in Fig. 1(d). Parameters a and b follow nearly linear relation with respect to E, while t₀ and σ increases quickly when E gets closer to the incident energy, as expected.²¹ 

In the following, we discuss the methods used to reconstruct a one-dimensional (1D) “super-resolution” signal from a DGS measurement with instrument broadening and noise. There are a few reports of the application of deconvolution methods on neutron spectra. Sivia, Silver, and Pynn showed that the maximum-entropy method was applicable to deconvolve the resolution function from neutron spectra,²⁴ and the asymmetric resolution function in pulsed source TOF spectrometers is beneficial in the deconvolution process due to higher frequency components of the sharp edge than symmetric, Gaussian-like resolution functions of the same full-width-at-half-maximum (FWHM). Weese et al. later reported that it was possible to deconvolve quasielastic neutron spectra using a Tikhonov regularization method.²⁵ In both efforts, the resolution function was assumed to be independent of energy transfer.

A spectrum measured at a DGS instrument can be in general written in terms of the resolution function and the true, resolution-free spectrum (we call it “ground-truth”) as

where f(E) is the measured spectrum, R(E′ − E; E) is the energy-dependent resolution function at energy transfer E, u(E) is the ground-truth spectrum, and n(E) is the noise term. A measurement is always discretized. Therefore, a discretized version of Eq. (1) is needed,

where f is the 1D data array for the measurement, R is the 2D resolution matrix, n is the noise array, and u is the ground truth array. This problem of obtaining u from f becomes ill-posed because of the noise and uncertainty in the PSF. Regularization techniques are usually employed to constrain the inverse problems like this to enforce some a priori knowledge about the ground truth u. The purpose of the reconstruction is, given the experimental data f and resolution function R, to find an estimate of the ground truth u that has improved resolution, by using signal-processing methods. Many image deblurring methods exist that can be leveraged for this purpose. We have adapted Lucy-Richardson (LR), one of the earliest deconvolution methods,^45,46 and more recent Linearized Bregman (LB),^47–49 and Split Bregman (SB)^50,51 reconstruction methods to the 1D datasets and taken into account the fact that the resolution function is energy-dependent. More mathematical and algorithmic details of the reconstruction techniques based on these three methods can be found in  Appendix E.

Figure 2(A) shows three common features in a phonon DOS curve: overlapping peaks, multiple separated peaks,⁵² and step functions. The reconstruction methods will help recover these features from noisy and smeared (convolved with resolution) DOS curves. The performance of the reconstruction methods is first tested against a dataset consisting of three synthetic DOS curves with these distinct, representative phonon DOS features, and the results are presented in Figs. 2(C) and 2(D). Shown in Fig. 2(C) are test results on synthetic datasets created by convolving common features with symmetric, Gaussian resolution functions and then adding Poisson noise. The first column is for the overlapping peaks. Since both peaks are relatively broad, the reconstructed data in panel (g) resemble well the original data in panel (a), regardless of the reconstruction algorithm. The second column is for multiple peaks with different widths. All algorithms perform well when the peak width is large. For sharper peaks, both the linearized Bregman algorithm and the split Bregman algorithm show oscillations at low intensity, while the Lucy-Richardson method performs well throughout the range. The third column is for step functions. Reconstruction results from all algorithms show some oscillations on top of the plateaus.

In direct geometry spectrometers at spallation sources, an important property of the instrument is that usually its energy resolution is asymmetric. The following test shows how this property can be taken advantage of. Shown in Fig. 2(D) are test results on synthetic datasets created by convolving common features with asymmetric resolution functions and then adding Poisson noise. The results are similar for broader peaks in columns one and two. Reconstructions in columns two and three for all algorithms show less oscillation than in Fig. 2(C). This can be attributed to the sharper edge of the asymmetric resolution function used here. The effects of the sharp edge of the resolution function in spallation neutron source data were originally discussed in Ref. 24.

In summary, the reconstruction algorithms work better for the synthetic overlapping-peak dataset in the first column than the synthetic multiple-peaks and step-function datasets in the second and third columns. The underlying reason, however, is that the multiple-peak dataset and the step-function datasets consist of higher-frequency components such as sharper peaks and edges, which are harder to resolve.

We then test the reconstruction algorithms with a synthetic DOS curve for graphite. This dataset is created by convolving the original DOS calculated from a DFT simulation with the ARCS instrument resolution function for E_i = 300 meV and adding Poisson noise. The results are shown in Fig. 3. The reconstructed DOS curves show clear improvements in recovering peaks at 200, 180, and 80 meV, and cliffs (step-functions) at 205 and 165 meV, compared to the resolution-smeared DOS curve. The split Bregman method does not work as well for peaks at 200 and 180 meV and is not the recommended method. A more quantitative analysis of the reconstruction efficacy of the algorithms can be found in  Appendix F.

We apply the SRDOS workflow to a phonon density of states dataset for a graphite sample measured at the ARCS instrument at SNS. Three different E_is, 30 meV, 130 meV, and 300 meV, were used for the measurement, and the resultant data were reduced to the phonon DOS using the standard procedure. A reconstruction algorithm is then applied to the three phonon DOS spectra measured using different E_is, taking into account the energy-dependent, asymmetric resolution functions to obtain three super-resolution DOS curves, which are then stitched together to obtain a full DOS curve. The results from three reconstruction algorithms are plotted in Fig. 4, along with the theoretical phonon DOS obtained from DFT simulation and the experimental phonon DOS. Compared to the experimental DOS at the bottom, the reconstructed DOS curves enhance peaks near 200, 155, and 105 meV, sharpen peaks near 180 and 80 meV, and recover more details in the energy range 120–160 meV. The DFT and measurement based phonon DOS curves agree to the level expected by the calculation technique,^53,54 demonstrating the power of the SR reconstruction. Among the three reconstruction algorithms, the Lucy-Richardson method and the linearized Bregman method show very similar results, while the split Bregman reconstruction shows an extra bump near 85 meV, which does not exist in the DFT DOS. The split Bregman reconstruction also shows slightly stronger fluctuation beyond 210 meV, the high energy cutoff for graphite phonons, than the other two reconstruction results.

Super-resolution phonon DOS spectra were obtained from S(Q, E) measured at ARCS using different incident energies. This is done by binning neutron event data in a sampling rate much higher than nominal instrument resolution, taking advantage of “sub-bin” shifts among detector pixels; modeling the resolution function with high accuracy; taking advantage of high frequency components in the sharp edge of the asymmetric DGS resolution function; and adapting and using reconstruction algorithms from SR imagery. For example, as long as enough counts are collected, the sampling rate for energy transfer can be safely increased to one bin per 0.1 meV, 1/100 of nominal energy resolution (10 meV) for an E_i = 300 meV measurement at the ARCS instrument, corresponding to ∼2 TOF bins of 0.1 μs at E = 150 meV. The sampling rate in the E axis can be further improved slightly when TOF data are converted to energy data by splitting events to neighboring bins or by interpolation. Another important factor is determined to be the sharper edge of the resolution function. For example, the sharper edge at E = 200 meV has a half width of ∼2 meV, and the effective resolution is smaller and estimated⁵⁵ to be ∼1.1 meV. It means that the resolution of the reconstructed data at E = 200 meV is ∼14 times better than the elastic resolution (∼15 meV) of the measured E_i = 300 dataset. With measurements using multiple incident energies, ∼5× resolution improvement across the energy range is readily accessible. This technique is limited by the signal-to-noise ratio as other reconstruction techniques, although for our synthetic datasets, a noise level less than ∼3% of maximum intensity is found to have little effect in reconstruction results. The difference of the resolution functions in datasets from different incident energies could be used to further refine the reconstruction results that are measured at multiple incident energies. This technique reduces the influence of instrument-specific resolution in the measured spectra and allows researchers to examine data closer to the “ground truth,” taking advantage of the existing research on image deconvolution and denoising. It makes use of information already in the current data, without going beyond the Nyquist frequency of the measurement. Similar techniques can be extended to higher dimensional DGS data, such as 2D S(Q, E) in powder measurements and 2D slice and 3D/4D volume data in single crystal S(Q, E) measurements, and possibly to other neutron scattering techniques, provided the three conditions are met: (1) a measurement sampling frequency higher than that of the nominal resolution; (2) sub-bin shifts for multiple measurements of the same data; and (3) the point spread function having a sharp feature.

The authors thank J. M. Borreguero Calvo, A. T. Savici, T. E. Proffen, M. K. Stoyanov, and T. R. Charlton for fruitful discussions. The graphite sample of grade G347A was funded by Strategic Planning Partnership between ORNL and Tokai Carbon Co., Ltd. (No. NFE-09-377 02345) with the U.S. Department of Energy. The DFT simulations were performed at the High Performance Facility in the University of Sharjah. This work was supported by the Department of Energy, Laboratory Directed Research and Development SEED funding, under Contract No. DE-AC05-00OR22725. This research used resources at the Spallation Neutron Source, a DOE Office of Science User Facility operated by the Oak Ridge National Laboratory.

The wide Angular-Range Chopper Spectrometer (ARCS) at the Spallation Neutron Source (SNS) at Oak Ridge National Laboratory²⁰ was used to measure the phonon density of states of the nuclear graphite grade G347A with a density of 1.85 g/cm³ at room temperature. The incident neutron energies E_i = 30, 130, and 300 meV were used to cover the whole energy transfer range (∼200 meV) in graphite. This measurement is a part of an extended study of the phonon DOS of nuclear graphite samples after irradiation by fast neutrons.⁵⁶ The full details of these measurements will be published elsewhere.⁵⁷ The neutron scattering data at 3 incident energies were collected as events and saved as NeXus⁵⁸ files, which are reduced by using Mantid⁵⁹ to S(Q, E) spectra. At each temperature, the multiphonon⁶⁰ package is used to iteratively improve a trial phonon DOS spectrum that fits the S(Q, E) spectrum best in the incoherent scattering approximation, taking multiphonon contribution into account. The three phonon DOS spectra can be “stitched” together by using the multiphonon⁶⁰ package to form one DOS spectrum. The stitching starts with the phonon DOS measured with the largest incident energy, $gEi=300(E)$⁠, covering the whole DOS spectrum. Moving on, the DOS curve is updated by replacing the low E portion of the DOS $gEi=300(E)$ with the partial DOS obtained from the E_i = 130 meV measurement, to form $gEi=130,300(E)$⁠. Similarly, we can perform the update for E_i = 30 meV measurement.

The vibrational properties of graphite are calculated using first-principles calculations as implemented by the VASP code.^61,62 A 6 × 6 × 1 supercell with 144 atoms was used to calculate the Hellman-Feynman forces. For the exchange-correlation functional, the generalized gradient approximation (GGA) of Perdew, Burke, and Ernzerhof formalism was used.⁶³ The integration over the Brillouin zone was performed using a 3 × 3 × 4 k-mesh⁶⁴ and an energy cutoff of 900 eV. The phonopy code⁵³ was used to calculate the phonon density of states.

In this work, we used the “dgsres” package⁶⁵ to calculate the resolution function at selected energy transfers by performing MCViNE simulations.⁴¹ One such simulation consists of an incident beam simulation, a simulation of neutron scattering from a powder-resolution sample that scatters neutrons only to a particular combination of Q, E (an example specification of the scattering kernel is presented in Fig. 5, and a simulation of detection of the scattered neutrons by the ARCS detector system that generates an event-mode NeXus file. The simulated data file is then reduced by Mantid⁵⁹ to obtain the I(E) spectrum, which would be the point-spread-function (PSF) at the Q, E point of choice. Since the energy resolution is nearly independent of Q, we only calculate the energy-dependent PSF for one Q value. Due to the shape of the dynamical range measured by a DGS instrument,⁶⁶ we choose a Q value to allow for calculation of larger E values.

The model⁴⁴ we choose to fit the resolution function is based on the Ikeda Carpenter function.²¹ In this model, the time distribution of neutron counts at the detector is written as

where

and

Here, l₁, l₂, and l₃ are moderator to Fermi-chopper, Fermi-chopper to sample, and sample to detector distances, respectively. The parameters in this model are a, b, R, which are related to the parameters in the original Ikeda-Carpenter model²¹ for the moderator, and σ for broadening caused by factors other than the moderator, including the sample and the detector, and t₀, a time offset parameter. It is then straightforward to transform the time distribution to energy distribution for comparison with MCViNE simulated PSF. We have chosen to keep the parameter R to be constant at 0.3.

The Lucy-Richardson (LR) method^45,46 is one of the classical methods for deconvolution. When it converges, it converges to the maximum-likelihood (ML) estimate of the solution of Eq. (2) of the main manuscript. It has found wide applications in deblurring images in, for example, astronomy.^67–69 

The adapted LR method uses the following iteration:

where R^* is the matrix for the flipped resolution functions. The initial condition for the iteration is

One of the state-of-the-art methods for restoring noisy and blurry images is the linearized Bregman method.⁴⁷ It is an iterative regularization procedure for solving the basis pursuit problem,

where ∥u∥₁ is the L₁ norm of u. It was proved^48,49 that the linearized Bregman iteration converges to the solution of

where ∥.∥ is the L₂ norm and μ is the regularizing parameter. The linearized Bregman iteration that we will adapt, analyze, and use here is generated by

where v^k is an auxiliary variable and δ is the step size of the iteration The shrink function

is the soft thresholding algorithm. The initial condition is u⁰ = v⁰ = 0.

It was proved⁴⁸ that the linearized Bregman iteration converges when $0<δ<δmax=2‖AAT‖$⁠. In general, a larger δ value converges faster than a lower δ value. Therefore, a δ value close to δ_max is usually used.

To compute the value of the parameter μ, we chose to balance the L1-norm and L2-norm regularization terms,

The stopping criterion for our denoising algorithm is when the residual is smaller than the error bar,

Here, σ is the error bar of the measurement f.

The Split Bregman method uses the Bregman iteration for minimizing the L₁ and L₂ terms in E4 separately by decoupling them. The decoupling is achieved by incorporating an alternative minimization scheme.^50,51 In this scheme, the L₁ term in E4 is replaced by an auxiliary variable “d” and the equation becomes

where ϕ is the linear operator of u. Equation (E10) is first minimized with respect to u by keeping d fixed, and in the next step, it is minimized with respect to d, while u is kept fixed. It has been shown that the alternate minimizing approach reaches the convergence with less number of iterations than the conventional approach.⁷⁰ The split Bregman method to solve the above problem is presented as follows:

where b is the Bregman vector and updated using E12 and λ is the regularizing parameter. Here, we use the gradient descend method and shrink function to solve u and d, respectively. Therefore, each iteration of the split Bregman algorithm consists of three steps:

where k is the iteration index and α is the step size. For the low values of α, the convergence is slow. μ is calculated using (E8) and the same principle is used to calculate λ from the following equation:

The iteration begins with the initial assumption of u⁰ = 0, d⁰ = 0, b⁰ = 0 and ends when Eq. (E9) is satisfied.

In order to further quantify the efficacy of the reconstruction methods, the root mean square error (RMSE) between the reconstructed dataset and the original dataset (ground truth) is calculated using the following equation:

where u_r, u_gt denote the reconstructed data and the ground truth, respectively. L is the size of the data. To make it easier to compare different datasets, each ground truth was scaled so that the maximum intensity is 1.

Figure 6 shows RMSE values for different datasets using different reconstruction algorithms. The RMSE values are used to measure how close the reconstructed results resemble the original data, and with Fig. 6, we can make the same conclusions as what we have observed by eye earlier in Figs. 2 and 3. It shows that RMSE is a good measure of reconstruction fidelity. For example, from Fig. 6, it can be seen that all three algorithms perform the best for the synthetic overlapping-peak dataset with small RMSE, and the disagreement (RMSE) increases for the synthetic multiple-peak and step-function datasets, regardless of the choice of resolution function (symmetric vs asymmetric). The Lucy-Richardson algorithm apparently performs better than the other two algorithms on the synthetic multiple-peak dataset. From Fig. 6, we also see that all algorithms perform similarly well for the synthetic graphite DOS without noise; however, the RMSE is in between the RMSE for the synthetic overlapping-peak dataset and that for the synthetic step-function dataset, demonstrating a more realistic DOS curve consisting of signals of varying features—both slowly varying valleys and peaks, as well as sharp peaks and cliffs. For the synthetic graphite DOS with noise, the split Bregman method yields a larger RSME than the other two methods.

Figure 7 and Table I present the computational costs. This figure shows the number of iterations and the computational times required to meet the convergence criteria for different datasets. The computational time is the “wall time,” which can be defined as the actual time for the program to finish from its start time, and it is dependent on the machine. This study was performed on a regular workstation with an Intel Xenon E5-2630 CPU. Roughly, the number of iterations and computation time follow a power law, and computation time increases with the number of iterations. It is clear from Fig. 7 that in general the split Bregman method needs fewer iterations to converge than the other two algorithms, meaning the split Bregman method is more efficient in each iteration getting closer to its solution. However, the computation time is in general still larger for the split Bregman method, meaning the efficiency in each iteration does not convert to overall efficiency. The linearized Bregman method and the Lucy-Richardson method show a similar number of iterations and computation time. All three algorithms are comparatively fast, requiring on the order of milliseconds.

In this section, we present a summary of the SRDOS technique. This section also serves as a quick guide of the SRDOS workflow. Figure 8 compares the SRDOS workflow to the multiframe SR imagery workflow. The first step of multiframe SR imagery is fusion. For SRDOS, this step can be achieved by choosing a E bin size smaller than the desired resolution in Mantid reduction. However, in order to gather enough statistics for a faithful reconstruction, the E bin size cannot be too small. For example, to achieve 3% error bar, it is necessary to have 1000 counts in one E bin.

Here, we give an example for estimating how long a measurement is needed for a SRDOS reconstruction. At E_i = 100 meV for the ARCS instrument, the flux on the sample is roughly 10⁵ counts/s cm⁻² MW⁻¹. If we assume the sample area is 10 cm² and the power is at 1 MW, there will be ∼10⁶ counts per second on the sample. Assuming a 10% scatterer, taking into account the solid angle coverage of the ARCS instrument, and the estimate that the inelastic scattering is often only 10% of the total scattering intensity, there will be ∼2000 counts/s for the DOS spectrum. If the desired resolution is 1% E_i (3–5× super-resolution), a E bin size of ∼0.25%E_i may be selected. This means on average there will be 5 counts/s per E bin. Collecting 200 s of data will be enough to bring an energy bin of average intensity to 1000 counts, and it means half of the bins will meet the statistics requirement. However, in order for the majority of the bins to reach the statistics requirement, and also in consideration of other sources of noise (from sample environments, for example), it is safer to use a longer (10×) time for data-collection. Furthermore, please note that the noise from the empty can measurement should be considered too if it is needed in data reduction. The time of data collection is in the order of an hour for this example, which is very much manageable, due to the fact that the ARCS instrument is a state-of-the-art high-flux DGS instrument.

Before reconstruction can be performed, resolution functions need to be calculated. For multiframe SR imagery, the resolution function is usually assumed to be a 2D Gaussian function with fixed widths across the image. The resolution is more complex in direct-geometry spectrometer data analysis. The measured intensity is often written as a convolution of a resolution function with the dynamical structure factor,

However, this definition hides the fact that R is a function that varies in form depending on the value of E (R is more asymmetric at larger E closer to E_i). A formula that explicitly shows this dependency can be written as

Another confusion about the resolution function in DGS is the direction of its long tail. In the usual convention [Eq. (H1)] for DGS, the long tail is on the positive side. Alternatively, the convention employed by the point spread function in imagery, namely,

could also be used. Such a change defines the long tail to be on the negative side. The reconstruction algorithm is particularly sensitive to the tail. Inconsistent usage of the resolution definition causes the reconstruction to produce results that are obviously incorrect.

After the DGS resolution function is calculated and fit as described in  Appendixes C and  D, the reconstruction can be performed using the methods described in  Appendix E.