Instrumental Resolution of the Chopper Spectrometer 4SEASONS Evaluated by Monte Carlo Simulation

TL;DR

This paper uses Monte Carlo simulations to evaluate and visualize the instrumental resolution of the 4SEASONS spectrometer, enhancing data analysis capabilities.

Contribution

The study introduces new simulation functionalities in Utsusemi for better visualization and prediction of the spectrometer's resolution function.

Findings

Simulations agree well with analytical calculations

Enhanced visualization of resolution functions implemented

Predictive capabilities for experimental configurations added

Abstract

We performed simulations of the resolution function of the 4SEASONS spectrometer at J-PARC by using the Monte Carlo simulation package McStas. The simulations showed reasonably good agreement with analytical calculations of energy and momentum resolutions by using a simplified description. We implemented new functionalities in Utsusemi, the standard data analysis tool used in 4SEASONS, to enable visualization of the simulated resolution function and predict its shape for specific experimental configurations.

Equations

$$\Delta E=2E_i{\left({\left\{\frac{\Delta t_c}{t_c}\left[1+\frac{L_1}{L_2}{\left(1-\frac{E}{E_i}\right)}^{\frac{3}{2}}\right]\right\}}^2+{\left\{\frac{\Delta t_m}{t_c}\left[1+\frac{L_3}{L_2}{\left(1-\frac{E}{E_i}\right)}^{\frac{3}{2}}\right]\right\}}^2+{\left[\frac{\Delta L_2}{L_2}\left(1-\frac{E}{E_i}\right)\right]}^2\right)}^{1/2}.$$

$$\begin{array}[]{rcl}Q_{x}&=&k_{i}\sin(\Phi_{i})-k_{f}\sin(\Phi_{f}),\\[3.0pt] Q_{z}&=&k_{i}\cos(\Phi_{i})-k_{f}\cos(\Phi_{f}),\end{array}$$

$$\Delta Q_x$$

$$\Delta Q_z={\left[(k_f\Delta {\Phi }_f\sin \Phi {)}^2+(\Delta k_i{)}^2+(\cos \Phi \Delta k_f{)}^2\right]}^{1/2}.$$

$$\Delta k_i=k_i\frac{{\left[(\Delta t_c{)}^2+(\Delta t_m{)}^2\right]}^{1/2}}{t_c},$$

$$\Delta k_f=k_f{\left[{\left(\frac{L_1}{L_1-L_3}\frac{\Delta t_c}{t_d-t_s}\right)}^2+{\left(\frac{L_3}{L_1-L_3}\frac{\Delta t_m}{t_d-t_s}\right)}^2+{\left(\frac{\Delta L_2}{L_2}\right)}^2\right]}^{1/2}.$$

$$\Delta Q_y={\left[(k_i\Delta {\Phi }_{i,v}{)}^2+(k_f\Delta {\Phi }_{f,v}{)}^2\right]}^{1/2}$$

Full text

1 INTRODUCTION

Estimating the instrumental resolution function of neutron scattering instruments is useful to determine physical quantities from experimental data. For inelastic neutron scattering instruments, the resolution function is expressed as an ellipsoid in the four-dimensional (4D) momentum and energy space [1]. For time-of-flight chopper spectrometers, the analytical method for calculating instrumental resolution has been established for a simple spectrometer [2]. This method has been extended to modern spectrometers [3]. However, precise calculations using analytical methods are complicated in the case of spectrometers with multiple optical components, advanced neutron guide shapes and multiple choppers. In this case, Monte-Carlo based numerical methods based can be used to evaluate the resolution function [4, 5, 6, 7]. Stimulated by preceding works, to estimate the resolution function of the 4SEASONS spectrometer, which is operating at J-PARC, we performed an McStas [4] simulation, based on the model described in Ref. [8].

4SEASONS is a time-of-flight direct geometry chopper spectrometer in the Materials and Life Science Experimental Facility (MLF) at J-PARC [8]. It is installed at the BL01 beam port for viewing the coupled moderator which is 18 m upstream from the sample position. Neutrons are transported by an elliptically converging guide tube coated with supermirrors, whose cross-section at the exit is 43.4 mm. The incident neutrons are monochromatized by a fast-rotating Fermi chopper positioned 1.7 m upstream from the sample position. In addition, 4SEASONS has a T0 chopper for suppressing fast neutrons and two disk choppers for band definition. Neutrons scattered by a sample are detected by the 19 mm-diameter and 2.5 m-long 3He position sensitive detectors placed cylindrically at 2.5 m from the sample position. 4SEASONS was designed for measuring spin and lattice dynamics in the 100–102 meV energy range [9]. It originally commissioned for a research project involving high-$T_{c}$ oxide superconductors. Now, its use has expanded to studies of other superconductors, magnetic materials, dielectrics, catalysts, and thermoelectric materials.

2 MONTE CARLO SIMULATION USING MCSTAS AND UTSUSEMI

2.1 Outline and Conditions of Simulation

The model of 4SEASONS for the McStas simulation is basically same as that used to simulate instrument performance [8], except for the sample and detector, which are replaced by the TOFRes_sample and the Res_monitor components, respectively. The original simulation model well reproduces the experimental neutron flux and energy resolution, as shown in Ref. [8]. The T0 chopper and the disk choppers were excluded from the simulation model. Considering that 2D magnetic excitations observed in materials such as copper oxide superconductors and iron-based superconductors account for most of the experiments performed using 4SEASONS [10, 11, 12, 13], in the present simulation, we assumed that the sample is a 2D system, and its 2D plane is perpendicular to $\mathbf{k}_{i}$. The sample is rod-shaped, and its diameter and height are denoted by $w_{s}$ and $h_{s}$, respectively. We studied samples of two sizes: one is $w_{s}\!=\!h_{s}\!=\!20\,\mathrm{mm}$ and the other is $w_{s}\!=\!h_{s}\!=\!40\,\mathrm{mm}$. The former was the size assumed in designing the instrument [8], and the latter is the typical size of the largest samples measured using the instrument. The assumed crystal structure is orthorhombic with lattice parameters $a\!\approx\!b=\!5.34\,\mathrm{\AA}$ and $c\!=\!13.24\,\mathrm{\AA}$. The 2D plane is defined by the $a$ and the $b$ axes. These values are typical for transition metal oxides with so-called K2NiF4 structures. The $[001]$ axis is parallel to the incident beam, and the $[100]$ axis is on the horizontal plane and perpendicular to the incident beam. This configuration is similar to those often employed in measurements of 2D systems [10, 11, 12, 13]. The incident energy ($E_{i}$) is 71 meV, and the rotating frequency of the Fermi chopper ($f$) is 250 Hz.

Simulations of the resolution function were performed using the TOFRes_sample and the Res_monitor components in McStas. The TOFRes_sample is a sample component for computation of the resolution function. It scatters neutrons isotropically within a specified angular range and selects the neutron energy uniformly so that neutron arrival times at the detector lie within a specified time bin. The TOFRes_sample is used together with the Res_monitor. It is a monitor component that records all scattering events. The output file of the Res_monitor is a list of “$\mathbf{k}_{i}$ (wave vector of incident neutrons), $\mathbf{k}_{f}$ (wave vector of scattered neutrons), position, intensity” of the neutrons that reach the detector at a time within the specified time bin [14]. Therefore, we need to specify the detector pixel position and the time corresponding to $\mathbf{Q}$ (momentum transfer) and $E$ (energy transfer) at which the resolution is estimated. The position of $\mathbf{Q}$ is specified by the 2D momentum transfer $\mathbf{Q}_{\mathrm{2D}}\!=\!(H,K)$ in reciprocal lattice units (r. l. u.). The component of $\mathbf{Q}$ parallel to $\mathbf{k}_{i}$ ($L$ in r. l. u.) and the scattering angle ($\Phi$) vary as functions of $E$ following momentum and energy conservations. We convert specified $\mathbf{Q}_{\mathrm{2D}}$ and $E$ to the corresponding detector position in real space and the neutron arrival time at the detector. The detector is a cylinder with a diameter of 19 mm and a height of 25 mm , which are the same dimensions as that of one detector pixel in 4SEASONS. The time bin at the detector is 1 $\mu$s.

To visualize the simulated resolution function, we converted the generated list of “$\mathbf{k}_{i}$, $\mathbf{k}_{f}$, position, and intensity” to a list of “$\mathbf{Q}$ (Å*-1*), $E$ (meV), intensity, and error”, and visualized it by using the $Q$-$E$ visualizer in Utsusemi for single-crystal experiments, VisContM [15]. Here $\mathbf{Q}$ and $E$ are given by $\mathbf{Q}=\mathbf{k}_{i}-\mathbf{k}_{f}$ and $E=E_{i}-E_{f}=\frac{\hbar^{2}}{2m_{n}}(k_{i}^{2}-k_{f}^{2})$, respectively, where $E_{i}$, $E_{f}$, and $m_{n}$ are incident and scattered neutron energies, and neutron mass. The $z$ axis is parallel to the incident beam, and the $x$ and $y$ axes are along the horizontal and vertical directions, respectively, in a right-handed coordinate system. By specifying crystallographic information in VisContM, momentum transfer can be visualized in r. l. u. in a manner similar to that in real experiments [Fig. 1(a)]. The original VisContM handles binary intensity data histogrammed as functions of detector pixel positions and $E$, and averages counts in a bin upon slicing data. However, it is not compatible with the list of resolution data which is a text file storing randomly generated (not histogrammed) $\mathbf{Q}$ and $E$ values of neutrons. Then, we added two new functions to VisContM: one is a function to read a text file in the form of [$\mathbf{Q}$ (Å*-1*), $E$ (meV), intensity, error] and the other is a new slice mode that simply sums all counts in a bin.

2.2 Results

Figures 1(b)–(d) show a simulated resolution function at $\mathbf{Q}_{\mathrm{2D}}\!=\!(1,0)$ and $E\!=\!0\,\mathrm{meV}$ for $E_{i}\!=\!71\,\mathrm{meV}$, $f\!=\!250\,\mathrm{Hz}$, and $w_{s}\!=\!h_{s}\!=\!20\,\mathrm{mm}$ when $\mathbf{c}\parallel\mathbf{k}_{i}$ and $\mathbf{a}$ is on the horizontal plane. To visualize the resolution function on 2D maps, we integrated it over two axes in the four momentum and energy axes [axes labeled “Thickness” in Fig. 1(a)] and plotted it as functions of the other two axes [axes labeled as “X-axis” and “Y-axis” in Fig. 1(a)]. The resolution function is generally expressed as an ellipsoid elongated and inclined with respect to the $H$, $K$, $L$, and $E$ axes. Figure 1(d) shows it is especially elongated along the $L$ or $E$ directions. This reflects the fact that the $L$ direction is parallel to the incident beam in the present simulation, and therefore $L$ and $E$ are correlated strongly. In other words, the resolution ellipsoid is elongated along the scan trajectory in the $\mathbf{Q}$-$E$ space [see the dotted line in Fig. 1(d)]. By contrast, when we see the resolution ellipsoid on the $H$-$K$ plane, it is symmetric with respect to the $H$ and $K$ axes. This means that in a scattering configuration with $\mathbf{k}_{i}$ and $\mathbf{k}_{f}$ in the horizontal plane, the horizontal and vertical components of the $Q$ resolution are not correlated.

To evaluate the simulated resolutions quantitatively, we projected the 4D resolution data onto one of the energy or momentum axes by integrating them over the other three axes. Then, we fit the resulting 1D resolution profiles to Gaussians to obtain their full widths at half maximum (FWHMs) as the resolution values. The closed circles in Fig. 2(a) show the energy resolution ($\Delta E$) for $w_{s}\!=\!h_{s}\!=\!20\,\mathrm{mm}$ as functions of $E$. The closed circles, open circles, and open diamonds in Fig. 2(b) show the momentum resolutions along the $H$, $K$, and $L$ axes, respectively, for $w_{s}\!=\!h_{s}\!=\!20\,\mathrm{mm}$ as functions of $E$. Given that in the present configuration the $H$, $K$, and $L$ axes are parallel to the $x$, $y$, $z$ axes, respectively, we represent their resolution components as $\Delta Q_{x}$, $\Delta Q_{y}$, and $\Delta Q_{z}$. The symbols in Figs. 2(c) and 2(d) show the energy and the momentum resolutions for the case of $w_{s}\!=\!h_{s}\!=\!40\,\mathrm{mm}$. As $E$ increases, both the energy and the momentum resolutions decrease (become better). For the momentum resolutions, $\Delta Q_{z}$ shows a faster decrease compared to $\Delta Q_{x}$ and $\Delta Q_{y}$. The increase in sample size from $w_{s}\!=\!h_{s}\!=\!20\,\mathrm{mm}$ to $w_{s}\!=\!h_{s}\!=\!40\,\mathrm{mm}$ increases $\Delta Q_{x}$ and $\Delta Q_{y}$, while it has little effect on the energy resolution and $\Delta Q_{z}$.

In addition to the widths of the resolution functions, Monte Carlo simulation is also useful for evaluating asymmetric tail in energy spectrum often observed at chopper spectrometers at a spallation source. The tail comes from the pulse tail at the moderator, and it is observed at the neutron energy gain side owing to the pin-hole camera effect of the chopper in the time-length space [16]. It is more visible on a spectrometer for viewing a coupled moderator such as 4SEASONS. The closed circles in Fig. 3(a) show the resolution function for elastic scattering integrated over all momentum axes, which clearly shows the tail on the neutron energy loss (positive $E$) side in contrast to the experiments [Tails can be seen similarly in Figs. 1(c) and 1(d).]. Here, we should recall that the observed intensity at a particular momentum transfer $\mathbf{Q}_{0}$ and $E_{0}$, $I(\mathbf{Q}_{0},E_{0})$ is a convolution of the scattering function $S(\mathbf{Q},E)$ and the resolution function $R(\mathbf{Q}-\mathbf{Q}_{0},E-E_{0})$, which is expressed as $I(\mathbf{Q}_{0},E_{0})=\int\!R(\mathbf{Q}-\mathbf{Q}_{0},E-E_{0})S(\mathbf{Q},E)\,d\mathbf{Q}dE$. If $S(\mathbf{Q},E)=\delta(E)$, as in the case of elastic incoherent scattering, $I(E_{0})=R(-E_{0})$. Then, the tail appears on the opposite side, that is, the negative $E$ side, in experiments as shown schematically in Fig. 3(b). The solid line in Fig. 3 shows the same resolution function inverted with respect to $E=0$, which coincides with the simulated scattering profile of vanadium (open circles).

3 COMPARISON WITH ANALYTICAL CALCULATION

To validate the resolutions evaluated by the Monte Carlo simulation, we performed analytical calculations of the energy and momentum resolution. In general, analytical calculation of the resolution ellipsoid in the 4D $\mathbf{Q}$-$E$ space is complicated. Here, instead, we consider a less general case of scattering in the horizontal plane, for which a streamlined calculation is used to obtain the energy and momentum resolution (FWHM). FWHM is proportional to its standard deviation if the distribution is a Gaussian.

3.1 Energy Resolution

We note here that the simplified expression used elsewhere [8, 17, 18, 19, 20] to calculate the energy resolution (FWHM) of the 4SEASONS spectrometer:

[TABLE]

$E_{i}$ and $E$ are the incident energy and energy transfer, respectively. $L_{1}$, $L_{2}$, and $L_{3}$ are the moderator-to-sample, sample-to-detector, and chopper-to-sample distances, respectively. For 4SEASONS, $L_{1}\!=\!18\,\mathrm{m}$, $L_{2}\!=\!2.5\,\mathrm{m}$, and $L_{3}\!=\!1.7\,\mathrm{m}$. $\Delta t_{c}$ and $\Delta t_{m}$ are the opening time of the Fermi chopper and the pulse width at the moderator, respectively. $\Delta t_{c}$ is effectively larger than its intrinsic value defined by chopper geometry, $\Delta t_{\mathrm{c0}}=d/(2{\pi}Df)$, because of angular divergence of the incident beam as $\Delta t_{c}=\Delta t_{\mathrm{c0}}\cdot p(\Delta\Phi_{i}^{\mathrm{max}})$. $D$ and $d$ are the diameter of the chopper rotor and the width of each of the slits, respectively, and $d/D\!=\!0.02$ for 4SEASONS. $\Delta\Phi_{i}^{\mathrm{max}}$ is the maximum angular divergence of the incident beam, and the function $p$ is described in Ref. [18]. $\Delta\Phi_{i}^{\mathrm{max}}$ has neutron energy dependence originating from neutron reflections by the supermirrors of the guide tube. It was estimated using the relationship between neutron wavelength and supermirror critical angle as described in Ref. [22]. $\Delta t_{m}$ varies with $E_{i}$ [23], and we used the numerical values of $\Delta t_{m}$ from Ref. [24]. $\Delta L_{2}$ is the uncertainty of $L_{2}$ resulting from the sample and the detector sizes. In general, $\Delta L_{2}$ depends on the sample shape, and the shape dependence should be considered seriously for high-resolution time-of-flight spectrometers [21]. Here, however, to consider the resolution of the middle-resolution spectrometer, we assumed $\Delta L_{2}=[(\pi w_{s}/4)^{2}+(\pi w_{d}/4)^{2}]^{1/2}$, where $w_{s}$ and $w_{d}$ are diameters of the sample and detector, respectively. $w_{d}$ is 19 mm for 4SEASONS. $\pi w_{s}/4$ ($\pi w_{d}/4$) is the average transmission length of a cylinder for a beam perpendicular to the cylinder. $t_{c}$ is the time at which neutrons with energy $E_{i}$ reach the Fermi chopper.

The solid lines in Figs. 2(a) and 2(c) show the $\Delta E$ calculated using Eq. (1) for $E_{i}\!=\!71\,\mathrm{meV}$ and $f\!=\!250\,\mathrm{Hz}$ in the cases of $w_{s}=20\,\mathrm{mm}$ and $w_{s}=40\,\mathrm{mm}$, respectively. They show excellent agreement with the results of the Monte Carlo simulations (closed circles). $\Delta E$ is dominated by the contribution of the first term in Eq. (1), and the other two terms have smaller contributions [17]. Therefore, the increase in sample size, which increases $\Delta L_{2}$, has little effect on the energy resolution, as shown in Figs. 2(a) and 2(c).

3.2 Momentum Resolution

Here, we concentrate on deriving the resolution components along the principal axes to compare them with the simulation results shown in Fig. 2(b). First, the horizontal components of the momentum resolution are calculated. The $x$ and $z$ components of $\mathbf{Q}$ are defined as

[TABLE]

where $\Phi_{i}$ and $\Phi_{f}$ are the angles of $\mathbf{k}_{i}$ and $\mathbf{k}_{f}$ relative to the primary beam centerline, respectively [Fig. 4(a)], and nominally $\Phi_{i}=0$. By differentiating Eq. (2) with respect to each parameter, we obtain the deviation of $Q_{x}$ from its nominal value as follows:

[TABLE]

where $\Delta\Phi_{i}$ ($\Delta\Phi_{f}$) is the angular divergence of the incident (scattered) beam on the horizontal plane, and $\Delta k_{i}$ ($\Delta k_{f}$) is the uncertainty in $k_{i}$ ($k_{f}$) [Fig. 4(a)]. Similarly, for $Q_{z}$, we obtain

[TABLE]

The contributions of each of the terms in Eqs. (3) and (4) can be understood geometrically from the scattering diagram in Fig. 4(a).

$\Delta k_{i}$ and $\Delta k_{f}$ are calculated by differentiating $k_{i}=\frac{m_{n}}{\hbar}(L_{1}-L_{3})/t_{c}$ and $k_{f}=\frac{m_{n}}{\hbar}L_{2}/(t_{d}-t_{s})$ with respect to $t_{c}$, $t_{s}$, and $L_{2}$, where $t_{s}$ ($t_{d}$) is the arrival time of neutrons at the sample (detector). Given that the uncertainties in $t_{c}$ and $t_{s}$ are $[(\Delta t_{c})^{2}+(\Delta t_{m})^{2}]^{1/2}$ and $\left\{[\Delta t_{c}L_{1}/(L_{1}-L_{3})]^{2}+[\Delta t_{m}L_{3}/(L_{1}-L_{3})]^{2}\right\}^{1/2}$, respectively,

[TABLE]

[TABLE]

By substituting Eqs. (5) and (6) into Eqs. (3) and (4), we obtain the expressions of $\Delta Q_{x}$ and $\Delta Q_{z}$ (they fall into the same expressions as written in Ref. [20] except for the $\Delta\Phi_{i}$ term in Eq. (3).). As for the angular divergences, $\Delta\Phi_{f}$ is determined by the sample size and the detector pixel size. By neglecting the shapes of the sample and detector, $\Delta\Phi_{f}=(w_{s}+w_{d})/(2L_{2})$. By contrast, it is difficult to describe $\Delta\Phi_{i}$ by using a simple function. We confirmed by simulations that although the maximum divergence of the incident beam can be well reproduced by the supermirror critical angle, which was used to calculate the $E$ resolution, the profile shape and the FWHM of the angular divergence deviate from the Gaussian and depends strongly on incident energy and sample size. Then, we fitted the simulated angular divergence of incident neutrons hitting the sample to a Gaussian, and used its FWHM as $\Delta\Phi_{i}$.

The vertical component of the momentum resolution $\Delta Q_{y}$ was obtained using the divergences of the incident and the scattered beams:

[TABLE]

where $\Delta\Phi_{i,v}$ and $\Delta\Phi_{f,v}$ are angular divergences of the incident and the scattered beams along the vertical direction, respectively. Given that the cross-section of the neutron guide tube of 4SEASONS has a square shape, $\Delta\Phi_{i,v}$ is equal to $\Delta\Phi_{i}$ in Eq. (3). As for $\Delta\Phi_{f,v}$, by using the sample height $h_{s}$ and detector pixel height $h_{d}$, $\Delta\Phi_{f,v}=(h_{s}+h_{d})/(2L_{2})$. $h_{d}\!=\!25\,\mathrm{mm}$ for 4SEASONS.

The solid, dotted, and broken lines in Figs. 2(b) and 2(d) show the calculated values of $\Delta Q_{x}$, $\Delta Q_{y}$, and $\Delta Q_{z}$, respectively, at $\mathbf{Q}_{\mathrm{2D}}\!=\!(1,0)$ ($a\!=\!5.34\,\mathrm{\AA}$) for $E_{i}\!=\!71\,\mathrm{meV}$ and $f\!=\!250\,\mathrm{Hz}$. They show reasonably good agreement with the simulations. The increase in sample size leads to an increase in $\Delta\Phi_{f}$ ($\Delta\Phi_{f,v}$) and $\Delta L_{2}$, resulting in an increase in the momentum resolution values. In the present calculations, the three terms in $\Delta Q_{x}$ have comparable contributions at $E\!=\!0$. As $E$ increases, the $\Delta\Phi_{f}$ and $\Delta k_{f}$ terms decrease because $k_{f}$ decreases and $\Phi$ increases, while the $\Delta\Phi_{i}$ term is independent of $E$. Therefore, the incident beam divergence is critical for the momentum resolution in such experiments, where $Q_{x}$ is of particular interest, which is often the case in experiments involving 1D or 2D systems. By contrast, in $\Delta Q_{z}$, the $\Delta k_{f}$ term is dominant, while the $\Delta\Phi_{f}$ term and the $\Delta k_{i}$ term are small and independent of $E$ (Note that $k_{f}\sin\Phi=Q_{x}$ is fixed in the present calculation.). The rapid decrease in $\Delta Q_{z}$ as a function of $E$ follows the behavior of $\Delta E$ in Fig. 2(a) and 2(c). This comes from the fact that $\Delta k_{f}$ is strongly correlated with the energy resolution. As for $\Delta Q_{y}$, it is determined only by the angular beam divergences, which are independent of $E$. $\Delta Q_{y}$ depends weakly on $E$ through the dependence of $k_{f}$ on $E$.

The difference between the simulations and analytical calculations probably comes from the rather simple modeling in the analytical case. As shown in Fig. 1, $\Delta Q_{x}$ and $\Delta Q_{z}$ are correlated strongly with the $E$ component of the resolution function. Neglecting these correlations may contribute to the discrepancy. Moreover, because $\Delta Q_{z}$ is dominated by the $\Delta k_{f}$ term, it can be affected by the asymmetric time structure of the neutron beam due to the moderator pulse tail. By contrast, in Fig. 2(b), $\Delta Q_{y}$ shows good agreement with the simulation results, possibly because $\Delta Q_{y}$ is simply expressed in terms of the angular divergence and the correlations with the other resolution components are small. However, as the sample size increases to $w_{s}\!=\!h_{s}\!=\!40\,\mathrm{mm}$, even the simulation and analytical results of $\Delta Q_{y}$ show a significant discrepancy. Eq. (7) for $\Delta Q_{y}$ is a square-root of the sum of the squares of the distributions of $k_{iy}$ ($\Delta k_{iy}$) and $k_{fy}$ ($\Delta k_{fy}$). Therefore, we assumes that $\Delta k_{iy}$ and $\Delta k_{fy}$ are not correlated. To check whether this assumption holds in the simulations, in Fig. 5, we plotted the simulated neutron distributions as functions of $k_{iy}$ and $k_{fy}$ for the two sample sizes. We found that a considerable correlation exists between $\Delta k_{iy}$ and $\Delta k_{fy}$, and there are fewer neutrons with $k_{iy}\sim k_{fy}$ in the $k_{iy}$-$k_{fy}$ map. Then, in the distribution of $Q_{y}=k_{iy}-k_{fy}$, the contribution of $Q_{y}\sim 0$ should be suppressed, while that of $|Q_{y}|>0$ should be enhanced, resulting in an effective increase in FWHM. This correlation comes from the fact that few neutrons that hit the sample upward (downward) and are scattered by the sample upward (downward) again can reach the detector pixel. As the sample size increases, the number of such unreached neutrons increases and this geometrical effect becomes significant in the momentum resolution. A similar geometrical correlation also exists between $k_{ix}$ and $k_{fx}$, and it should become one of the causes of the discrepancy in $\Delta Q_{x}$ between the simulated and the analytical results.

4 SUMMARY AND REMAINING ISSUES

We performed Monte Carlo simulations of the resolution function for the 4SEASONS spectrometer by using McStas. Visualization of the obtained 4D resolution functions by using Utsusemi, the data analysis suite for MLF, facilitates intuitive understanding of the resolution functions. Furthermore, we performed analytical calculations of the energy and momentum resolutions. Although they were performed only for a simple model, reasonably good agreements with the simulation results not only validate the simulations, but also prove that such a simple calculation is quite useful for obtaining a rough estimate of the resolution. We plan to develop the simulation toward more general cases of 3D systems and enhance compatibility with the data analysis software.

Finally, we mention aspects that are important for resolution analysis of experimental data but were not considered in the present study. In this study, we performed the simulation with the smallest time bin at the detector as possible to obtain the intrinsic energy resolution. In real experiments, however, data is often cut with a finite width of $E$ (roughly, 1–10 meV). To consider the width of the energy cut, we must use a finite time bin at the detector that corresponds to the width of the energy transfer. Another important aspect is sample mosaicity. The resolution obtained by McStas is the intrinsic resolution determined by the geometrical arrangements of the instrument components. In this sense, the effect of sample mosaicity is effectively included in the dynamic structure factor $S(\mathbf{Q},E)$ of the sample. However, it is often useful to combine the effects of sample mosaicity with the resolution to deduce the true $S(\mathbf{Q},E)$ representing physical properties [25]. However, it seems that the Res_sample and the TOFRes_monitor components in McStas cannot treat sample mosaicity directly, and therefore, some modifications to the simulation codes are required. It would be useful to consider the simple case that sample mosaicity is represented as a rotation of $\mathbf{Q}$ by $\Delta\theta$ [Fig. 4(b)]. Then, the resolution perpendicular to $\mathbf{Q}$ is most affected by the sample mosaicity. For example, when we are interested only in $Q_{x}$, as is the case in measurements of 1D or 2D systems, the term originating from the sample mosaicity $\Delta Q_{\mathrm{mos}}$ is $Q\Delta\theta\cos[\arcsin(Q_{x}/Q)]$. It has little effect on $\Delta Q_{x}$ when $E$ is sufficiently small, but it should be considered as $E$ increases.

5 ACKNOWLEDGMENTS

We thank K. Nakajima, M. Nakamura, G. E. Granroth, P. Willendrup and E. Farhi for valuable discussions and useful advices. R. K. and Y. I. were supported by JSPS KAKENHI Grant Number 15K04742.