Monte Carlo simulation of the resolution volume for the SEQUOIA spectrometer

Monte Carlo ray tracing simulations, of direct geometry spectrometers, have been particularly useful in instrument design and characterization. However, these tools can also be useful for experiment planning and analysis. To this end, the McStas Monte Carlo ray tracing model of SEQUOIA, the fine resolution fermi chopper spectrometer at the Spallation Neutron Source (SNS) of Oak Ridge National Laboratory (ORNL), has been modified to include the time of flight resolution sample and detector components. With these components, the resolution ellipsoid can be calculated for any detector pixel and energy bin of the instrument. The simulation is split in two pieces. First, the incident beamline up to the sample is simulated for 1 × 1011 neutron packets (4 days on 30 cores). This provides a virtual source for the backend that includes the resolution sample and monitor components. Next, a series of detector and energy pixels are computed in parallel. It takes on the order of 30 s to calculate a single resolution ellipsoid on a single core. Python scripts have been written to transform the ellipsoid into the space of an oriented single crystal, and to characterize the ellipsoid in various ways. Though this tool is under development as a planning tool, we have successfully used it to provide the resolution function for convolution with theoretical models. Specifically, theoretical calculations of the spin waves in YFeO3 were compared to measurements taken on SEQUOIA. Though the overall features of the spectra can be explained while neglecting resolution effects, the variation in intensity of the modes is well described once the resolution is included. As this was a single sharp mode, the simulated half intensity value of the resolution ellipsoid was used to provide the resolution width. A description of the simulation, its use, and paths forward for this technique will be discussed.


Introduction
Understanding the instrumental contribution to neutron experiments is critical for both interpreting the measured data and for planning new measurements. Traditionally this task has been carried out using analytical methods [1][2][3][4][5] but as instruments become more complex this becomes a significant effort. Another approach is to use Monte Carlo Ray tracing to simulate the instrument and the sample together. This is one of the standard methods of analysis for triple axis spectrometers with the Restrax package [6,7], but becomes more difficult on pulsed spallation source direct geometry spectrometers due to the asymmetric pulse shape and the large numbers of pixels. The advent of generalized Monte Carlo simulation packages such as McStas [8,9] McVine [10], Vitess [11][12][13] and NISP [14] make this case easier to realize and they have been particularly fruitful in several cases [15][16][17][18]. However sample kernels are hard to develop and the simulations can be very time intensive; especially if one wants to map out volumes of momentum transfer (Q) and energy transfer (ω) space.
In this contribution we describe a third route. This method builds on work for the European Spallation Source (ESS), where the resolution volume R(Q, ω) for direct geometry spectrometers was simulated for instrument design purposes [19], to an operating instrument at the a e-mail: granrothge@ornl.gov SNS. In both the ESS and this work, Monte Carlo Ray tracing was used to calculate R(Q, ω) for a specific detector pixel and ω bin. Specific to this work the matrix transformation was applied to put the ellipsoid in the same frame as the crystal. Once this is done two paths forward are discussed. One is to use this for guidance in experiment planning. The second is to calculate resolution ellipsoids along a dispersion surface for a sample and then convolute these ellipsoids with a model function for comparison with data. A few ways of characterizing the calculated resolution ellipsoid will also be discussed.

Monte Carlo Ray tracing
Detailed McStas [8,9] Monte Carlo simulations were developed during the design process of the SEQUOIA spectrometer [20,21]. These simulations were verified against initial measurements of the spectrometer showing excellent agreement [22]. More recent simulations were performed using McVine [10] for a sample of uranium nitride (UN) measured on SEQUOIA. The input to McVine was the same McStas simulation used in the instrument design phase. The results of these simulations show remarkable agreement with the data [18]. Therefore we are confident that this model is a good starting place for simulations to check the resolution volume of the spectrometer. For efficiency, the instrument was broken into two parts, where the dividing line is just before the EPJ Web of Conferences sample position. This is the position where the virtual monitor and the virtual source are placed for the front end and back end simulation, respectively. However the front end filters most of the neutrons and in the case below only about 1% of the neutrons generated by the moderator make it to the virtual source. The number of events needed as input to the back end simulations is 1 × 10 9 . Therefore simulations with 1 × 10 11 neutrons were run using a 30 core machine. The virtual source generated by this component is used for all subsequent simulations and thus makes the overall simulation relatively efficient. The instrument front end is the exact same code that was used in the earlier work up until the virtual monitor.
The back end of the instrument is a simple three component simulation consisting of a virtual source, a TOFRes sample component and a Res monitor component. These later two components work together to produce a list of all k i , k f , and probability that the neutron scatters into that detection pixel. This second simulation takes an angle about, and a position along, the vertical axis, and an energy transfer bin as an input. It then positions the Res monitor component at that location and tracks the corresponding time bin to calculate the resolution for a specific detector pixel. The focusing feature of McStas [23] is used so neutron events are calculated only for the desired detector pixel. The run time for a given resolution volume on a single core is less than 30 s. For instrument planning, this procedure is followed for a few individual pixels. For the measured excitation case points are calculated along the specific dispersion curve. The resultant neutron packets for each calculated resolution volume were then loaded into a customized NumPy [24] routine that applies UB [25] and goniometer matrices to convert from lab coordinates to crystal coordinates. The events are histogrammed into a 100 bins per non projected dimension during the projection step. When visualizing resolution ellipsoids, typically 2 dimensions are projected out to examine the volume. Tests for this case are detailed in the next section. However this class also provides methods to only project out 0 or 1 axis and thus the volume can be viewed in the Mantid extensions to paraview [26]. The one axis case is shown in Fig. 1. The zero axis case would involve a movie and is out of the scope of this contribution.

Characterizing the resolution shape
Once the points for the resolution volume are generated, the next step is to characterize its shape. Three methods will be discussed here: first taking the FWHM of the distribution, second using the convex hull [27] of the generated points at a certain probability level, and third taking the minimum volume ellipsoid around the simulated points. All these methods assume the fundamental data is characterized by a resolution volume normalized to unity. It is straightforward to look at the bounds of the resolution volume at a given probability level. The first approach is to look at anything contained within the FWHM of the probability distribution. This approach, though rather coarse, is quite effective for sharp modes as is demonstrated by the YFeO 3 example [28] given in Sect. 4. This is a reasonable approach if one has generated sufficient statistics to characterize the resolution volume width and if one is studying sharp modes with a relatively flat and weaker background. For the study of weak modes in the presence of nearby sharp modes, the intensity at the 1% level may be more important to consider. To that end, we take the approach of characterizing the resolution volume at various probability levels. One characterization is the convex hull [27] of the points. This is the npolytope determined by the extremal points of the n dimensional volume. It works well if sufficient statistics are generated at the probability level. However there is no simple parameterization of an arbitrary polychoron.
A simpler shape that gives a reasonable approximation to the resolution volume, minimizes the amount of simulation time and is easily parameterized is desirable. To this end, we investigated the minimum volume ellipsoid about the measured points at a specific probability level. The reasoning behind this choice is that generally resolution volumes are assumed to be ellipsoidal [5], one needs fewer simulated points, and the shape is easily described for use with other software. A full description of computing the minimum volume ellipsoid about a number of points is beyond the scope of this document and is well described in the literature [29][30][31]. Basically the idea is to expand the dimension of the problem by 1 to map it onto a hypersphere, next a variable transformation is used to put this into a problem that can be solved by linear least squares. Since only the outer most points in the volume matter in the process of calculating the ellipsoid, the convex hull is calculated first to reduce the size of the matrix in the least squares problem. This procedure was implemented in python by porting a matlab routine [31] available in the public domain. Figure 2 show contours for a resolution volume integrated over two of the Q dimensions. The 1% contours are represented in blue and the 50% contours are represented in red. The green and magenta curves are the minimum volume 03006-p.2 Figure 2. The minima volume ellipses around the 1% contour and the 50% contour of the resolution ellipsoid. The blue curve is the 1% contour and the green ellipse is the minimum volume ellipse for that contour. The red curve is the 50% contour and the magenta ellipse is it's minimum volume ellipse. ellipses calculated by the above procedure for the 1% and 50% contours, respectively. This procedure could be performed without integrating over the dimensions, however the two dimensional version is used for ease of explanation. Note the minima volume ellipse works well for the 50% contour, but does not work well for the 1% contour. Clearly the sharp rise and long tail of the Ikeda-Carpenter function [32], means the volume is non-ellipsoidal in this range. Therefore for sharp modes in a relatively flat background the 50% minimum volume ellipse is a reasonable approximation to the resolution volume. However the convex hull is a better representation for planning purposes where the mode is weak and strong sharp features are nearby.

Comparison to experimental data
The main reason for this work is to provide users with a method of estimating resolution volume for experiment planning. Nevertheless, the calculated resolutions volumes can be used to provide a R(Q, ω) to convolute with a specific model. As an example YFeO 3 was measured on SEQUOIA [28]. The result for the excitations in the 2 k3 direction are shown in Fig. 3a [28]. The spin rotation technique [33,34] was used to calculate the expected spin wave spectra and the model result is shown in Fig. 3b. The color scale in (b) was chosen to clearly show the dispersion that is visible in the measurement (in blueish tints) from that which is not (yellow to reddish tints). When comparing Figs. 3a and b, there are two major differences. The first difference is there is little intensity variation over the mode going from k = 1 to k = −1 in the data. However the model S(Q, ω) has more than an order of magnitude drop when going from |k| = 1 towards 0. The second difference is that there is a mode that is observed near |k| ∼ 1 and ω = 68 meV. To sort out if these observations are resolution effects, we calculate the resolution volumes at several points along the measured Figure 3. a) The measured data of the spin wave modes in the 2 k3 direction. b) The model for YFeO 3 as described in the text. Only the structure factor values represented by the blueish features are strong enough to be observed above the instrumental background. The color scale is reversed from the other plots to visually emphasize where the modes are observed. c) A plot of the resolution ellipsoids calculated at points along the measured dispersion in (a). The blue, green, and red colors are the 1%, 50% and 80% probability, respectively. The model in (b) convoluted with the resolution ellipsoids determined from the simulated resolution volumes in (c).

EPJ Web of Conferences
dispersion. An overview of these calculation results is provided in Fig. 3c. The first thing to note from this figure is that the 1% tails from the 68 meV do not overlap with the very weak mode that is at |k| ∼ 1 and ω ∼ 60 meV in the model. Therefore this signal comes from an extra mode that is described in detail in reference [28].
The other major result is best seen when the FWHM pulled from the 50% contours of Fig. 3c is convoluted with 3b to result in 3d. Now the intensity variation in the model matches what is observed in the data as |k| is varied from 1 to 0. Note that only points along the dispersion needed to be calculated. Therefore the simulation is rather efficient. In fact examination of Fig. 3b shows that the resolution ellipsoids change slowly with k and ω. Therefore fewer ellipsoids could be calculated with out affecting the result. A study of how sparsely the ellipsoids can be calculated is under way.

Conclusions
The use of Monte Carlo ray tracing to calculate the resolution volume for the SEQUOIA spectrometer has been demonstrated. Three methods for characterizing an approximate resolution volume has been presented. The FWHM method was used to accurately model the resolution for the YFeO 3 experiment, the minimum volume ellipsoid method has been described and will be used going forward for planning purposes and for sharp modes. The convex hull method was mentioned for completeness. For diffuse and weaker modes a more complete characterization of the resolution ellipsoid is under study. Work will continue to make these tools generally accessible to the broader user community.