pathSQE: an automated workflow for single-crystal inelastic neutron scattering data processing and analysis

Inelastic neutron scattering (INS) experiments utilizing modern time-of-flight spectrometers enable the comprehensive mapping of the energy (E)- and momentum (Q)-resolved dynamical structure factor of single crystals, probing both the lattice and magnetic excitations. Yet, the large size and complexity of four-dimensional INS data are challenging current analysis workflows, often resulting in an underutilization of the measured information. To help address this issue, this paper introduces new software interfaced with the Mantid framework, pathSQE, designed to streamline the processing, analysis and interpretation of 4D single-crystal INS data. By automating key tasks such as 1D/2D slicing, symmetrization, Brillouin zone folding, data visualization, prioritization and filtering, and comparisons with simulations, pathSQE facilitates and accelerates INS data analysis workflows. This paper outlines the features and implementation and provides several illustrations of the use of pathSQE on data collected on single crystals using direct-geometry time-of-flight spectrometers at the Spallation Neutron Source, including Ge, FeSi, MnO and SnS single-crystal measurements on the ARCS, HYSPEC and CNCS neutron spectrometers. Beyond streamlining post-experiment data processing, pathSQE establishes an automated and modular processing pipeline that could support future real-time experiment steering.

1. Introduction

With wavelengths and energies similar to atomic spacings and fundamental excitations in condensed matter, cold and thermal neutron beams constitute powerful direct probes of both the static structure and dynamics in materials. Inelastic neutron scattering (INS) has long been the workhorse technique to investigate the momentum and energy dependence of both lattice vibrations (phonons) and magnetic fluctuations (magnons) across reciprocal space (Placzek & Van Hove, 1954; Marshall & Lowde, 1968; Chatterji, 2005). The dispersions and linewidths – or more complex spectral functions – of such excitations can be directly and quantitatively inferred from the measured intensities, offering deep insights into the elementary excitations governing material behavior (Brockhouse & Iyengar, 1958; Pepy, 1974; Squires, 1978). Traditionally, INS experiments have been performed using triple-axis spectrometers at reactor-based neutron sources (Brockhouse, 1995), wherein measurements are sequentially acquired in predefined limited regions of interest for the wavevector transfer, Q, and energy transfer, E. Although triple-axis spectrometers remain widely utilized and are typically retrofitted with multiplexed analyzers to accelerate reciprocal-space mapping, the advent of high-power pulsed spallation neutron sources combined with highly multiplexed detector arrays, such as arrays of pixelated ³He tubes covering large solid angles subtending several steradians, has revolutionized INS experiments. A measurement of a single crystal with a fixed orientation on time-of-flight direct-geometry chopper spectrometer (DGS) instruments (Abernathy et al., 2012; Bew­ley et al., 2006; Bewley et al., 2011; Ehlers et al., 2011; Granroth et al., 2010; Itoh et al., 2011; Kajimoto et al., 2011) probes the dynamical structure factor S(Q, E) on a curved 3D surface in the 4D Q, E reciprocal space. By rotating the crystal (typically azimuthally) and combining measurements acquired over a series of crystal orientations, one maps the full 4D S(Q, E) over large Q, E volumes. A 3D subset (for a plane in Q space) of such a 4D dataset from an experiment performed on the ARCS spectrometer (Abernathy et al., 2012) at the Spallation Neutron Source (SNS), Oak Ridge National Laboratory (ORNL), is shown in Fig. 1 as an example. Depending on the incident neutron energy, the instrumental geometry and the crystal structure of the sample, the measured data volumes can span many Brillouin zones (BZs) – hundreds or even thousands of BZs in the case of thermal neutron beams. Furthermore, in the event mode of data collection employed at newer spallation sources, the details of each individual neutron detection are saved, and typical time-of-flight INS datasets can comprise billions of unique scattering events (Granroth et al., 2018).

Figure 1 Experimental time-of-flight INS S(Q, E) volume for a Ge single-crystal dataset measured with the ARCS spectrometer at the SNS. The out-of-plane [0, 0, L] direction is integrated between −0.1 and 0.1 reciprocal lattice units (r.l.u.)

Despite the wealth of information contained in time-of-flight INS datasets, extracting insights remains a challenge because of a combination of impeding factors in current analysis workflows, as illustrated in Fig. 2. First and foremost, the inherent complexity of properly handling event-based INS data requires software environments designed to handle these intricate yet important details. Specifically, to obtain a physically meaningful quantity like S(Q, E), the raw event data must undergo data reduction. This includes multiple processing steps such as transformation from instrument coordinates to the (Q, E) space of the sample and statistically rigorous combination, histogramming and normalization while carefully propagating uncertainties (Savici et al., 2022), as implemented in the Mantid package (Arnold et al., 2014). In the current practice, extracting information from reduced INS data involves painstakingly locating, processing and analyzing specific 1D cuts or 2D slices from the full 4D dataset in a labor-intensive manual process. Several software packages are commonly used for this purpose, such as Mantid (Arnold et al., 2014), Horace (Ewings et al., 2016), DAVE (Azuah et al., 2009) or Mslice (Coldea, 2000). Across all of this software, however, each time a new slice is desired, users must manually determine and specify descriptors defining the processing method and which subset of the full 4D dataset to process, such as the necessary Q projection axes and Q and E end points, bin sizes, and integration limits among other information. This generally creates a bottleneck in the data processing and analysis process, highlighted in red in Fig. 2. This manual and effort-intensive workflow often also requires prior knowledge of the expected scattering features and their Q and E location, at the risk of underutilizing the rich experimental datasets. Furthermore, its reliance on human intuition and trial-and-error exploration limits both efficiency and reproducibility (Noack et al., 2021). To rationalize the features probed in each experimental slice and visualization, it is often useful to also perform analogous theoretical simulations for comparison, such as computations of S(Q, E) based on lattice dynamics or spin-wave theory (Delaire et al., 2011a; Delaire et al., 2011b; Headings et al., 2010; Oh et al., 2013). However, accurately capturing experimental conditions requires significant time and expertise, limiting both the accessibility and throughput of this analysis (Han et al., 2024). Finally, combined with these technical limitations, the sheer scale of the dataset, in terms of both its high dimensionality and memory requirements, further hampers the ability of researchers to comprehensively analyze the INS data volume (Ratner et al., 2019; Doucet et al., 2021). These challenges have also sparked growing interest in applying machine learning (ML) methods to INS datasets – for tasks such as theoretical model optimization, data-driven exploration, and enhanced feature extraction and representation (Doucet et al., 2021; Han et al., 2025).

Figure 2 In the standard INS analysis workflow, experimental details and slice descriptors are used to process the event data and generate S(Q, E). Manually defining each slice is time and effort intensive, as shown in red, limiting both post-experimental data extraction and on-the-fly analysis for experimental steering.

Several codes have been developed to address some of these challenges. The recently developed program SHIVER (Savici, 2025) provides an intuitive graphical user interface – specifically designed to facilitate single-crystal time-of-flight INS analysis within Mantid – wherein users can take and visualize slices by populating fields rather than scripting. While extremely useful for fine-tuned analysis of individual slices, it lacks the capabilities to enable more automated, comprehensive analysis. The Phonon Explorer code (Reznik & Ahmadova, 2020) addresses some of these challenges by automating the Mantid processing and plotting of 1D S(E) cuts at symmetrically equivalent Q points throughout an INS data volume. While a valuable step toward automating INS data analysis, it remains limited with regards to 2D slicing capabilities or algorithmic filtering of data subvolumes, and simulation integration. On the simulation side, programs such as Phonopy (Togo, 2023; Togo et al., 2023), Simphonies (Bao et al., 2016), Euphonic (Fair et al., 2022), OCLIMAX (Cheng et al., 2019), SpinW (Toth & Lake, 2015) and SUNNY (Dahlbom et al., 2025) provide efficient methods for calculating INS spectra. Additionally, Euphonic and SpinW integrate with Horace, the former via the Horace–Euphonic interface. However, despite their ability to closely align simulations with experimental data, these programs still typically rely on manual processing of individual slices. To fully capitalize on the high information content of INS datasets from experiments using modern time-of-flight spectrometers with large detector arrays, there is a critical need for more efficient, systematic and automated analysis workflows that not only integrate effectively with established software like Mantid (Arnold et al., 2014) but also go beyond basic automation by incorporating data-driven selection strategies and facilitating direct experiment–theory comparisons (Han et al., 2025).

To address current bottlenecks in the exploration and analysis of single-crystal INS datasets from modern DGS instruments, we introduce pathSQE, a software tool designed to streamline the processing and interpretation of such INS datasets while integrating with the well-established INS data protocols and algorithms in Mantid. At its core, pathSQE automates data slicing and visualization using the auto-slice maker algorithm, requiring only the desired Q point coordinates to reduce, histogram and normalize the corresponding data. Additionally, using these same simple inputs, it systematically generates individual slices across all symmetrically equivalent Q points or segments, reducing the reliance on manual selection for 1D cuts and 2D slices. Beyond producing individual slices, the code implements a BZ folding algorithm that enhances measurement statistics by combining symmetrically equivalent regions (as illustrated in Fig. 3), harnessing the normalization and uncertainty propagation methods in Mantid to ensure statistical consistency. Finally, pathSQE integrates experimental and theoretical workflows by incorporating automated parallel phonon simulation capabilities that employ analogous Q sampling, histogramming, experimental masking and resolution blurring, enabling direct experiment–theory comparisons. Overall, pathSQE reduces the required manual effort across experimental and simulation workflows while leveraging the full 4D S(Q, E) volume provided by time-of-flight INS measurements. Beyond improving individual dataset analysis, pathSQE is designed to support broader advancements in INS by providing a scalable processing workflow for extracting insights from increasingly complex INS datasets.

Figure 3 Illustration of Brillouin zone (BZ) folding using a simple two-dimensional cubic lattice as an example. The unfolded dataset extends over multiple BZs, indicated by black boxes. Reciprocal lattice points () are shown as filled red circles, with the origin (O) marked by an open red circle. In the translational BZ stacking step, data points are shifted back into the first BZ (gray region) by applying the relation , where Q is the total momentum transfer and is the nearest reciprocal lattice vector. In the subsequent point-group folding step, the symmetries of the crystal are applied to further map the data into the irreducible wedge of the first BZ (orange region).

2. Features and implementation

2.1. Workflow overview

A high-level schematic of the pathSQE workflow is shown in Fig. 4. The code requires two primary inputs: a pathSQE input file (see Appendix A1) and a dataset configuration file. Once the relevant Q regions are identified, the code uses the auto-slice maker algorithm to process each in turn. From this core operation, several optional branches can be invoked, as illustrated in Fig. 4. These include slice filtering and evaluation, harmonic phonon S(Q, E) simulations, and BZ folding. All generated data products, both experimental and simulated, are saved using an intuitive directory structure with systematically named, searchable files. Additionally, summary Reports are automatically generated and saved to visualize the full set of outputs. In the sections that follow, each major feature from the workflow diagram is discussed in detail and demonstrated using a representative dataset collected on a Ge single crystal using ARCS at the SNS (a view of these data is shown in Fig. 1) and several other datasets from the CNCS and HYSPEC spectrometers, also at the SNS.

Figure 4 Summary of the pathSQE workflow, showing the main capabilities and their relationships. Bold bounding boxes denote inputs and outputs, and dashed gray arrows indicate conditional paths.

2.2. Auto-slice maker algorithm

2.2.1. Slice descriptor generation

The primary objective of the auto-slice maker algorithm is to minimize the need for human input by requiring only the Q point(s) of interest to automatically determine accurate slice descriptors. Using this information, it generates a corresponding slice description Python dictionary and calls the appropriate Mantid functions to process the data. The auto-slice maker is capable of producing a 1D S(E) spectrum at a specific point Q_in, or a 2D S(Q, E) spectrum along a 1D path spanning from Q_1,in to Q_2,in. Additionally, multiple Q segments or points of interest can be specified in the input file, all of which are sequentially processed, saved and visualized. This functionality serves as a fundamental building block upon which other features of the pathSQE workflow are built.

The slice descriptor determination process can be broken down into a few main steps. In the 2D S(Q, E) case, on the basis of the two points Q_1,in and Q_2,in that define the end points of the desired 1D segment, the difference vector ΔQ is calculated, normalized and used to define the first Q projection axis, Q_dim1. The second projection axis, Q_dim2, is then chosen to lie in a reciprocal plane that is parallel to the horizontal scattering plane. Finally, the third, Q_dim3, is simply the cross product of the first two to ensure that the voxels for histogramming and normalization are consistently rectangular. The input Q_1,in and Q_2,in are multiplied by the transformation matrix to express the segment end points in the new slice Q basis as Q_1,slice and Q_2,slice. The next step in fully describing the slice is to specify the Q end points, bin sizes and integration ranges along each dimension. For Q_dim1, the bins are , whereas for the integrated dimensions, the bins are and  − . For the simpler 1D S(E) case, all Q integration ranges are set as for i = 1, 2, 3. In both the 1D and 2D cases, the energy binning parameters are extracted from the input file. Finally, when the information needed to specify the desired slices or cuts has been determined, the MDNorm function (Savici et al., 2022) is called using the Mantid Python API to slice and normalize the data.

2.2.2. Demonstration on Ge data measured with ARCS

To illustrate the functionality of the auto-slice maker algorithm, we applied it to a Ge single-crystal dataset collected using ARCS, as shown in Fig. 5. A common approach to visualizing INS intensity involves constructing a 1D piecewise Q path composed of high-symmetry directions in reciprocal space. Germanium crystallizes in the diamond face-centered cubic (f.c.c.) structure (space group , No. 227). Its unit cell and BZ are shown in Figs. 5(a) and 5(b), respectively. Using the seeKpath package (Hinuma et al., 2017; Togo et al., 2024), pathSQE identifies a standardized high-symmetry path in the first BZ, defined by the coordinates of high-symmetry k-points in the primitive reciprocal-lattice basis (Setyawan & Curtarolo, 2010). For Ge, the resulting path and symmetry points are annotated in Fig. 5(b). This set of q-space coordinates was subsequently translated to various BZs and passed as input to the auto-slice maker, generating the processed data and corresponding aggregated visualizations in Figs. 5(c)–5(f). Following this, the code was re-executed using only the W points as input, yielding a set of representative 1D S(E) cuts [Fig. 5(g)], which are often used to support more quantitative analysis once features of interest have been identified.

Figure 5 Auto-slice maker workflow applied to a single-crystal dataset for Ge measured with the ARCS spectrometer at the SNS. (a) Conventional unit cell and (b) BZ of Ge. Automated 2D inelastic slices and visualizations of S(Q, E) along the high-symmetry q path translated to BZs (c) [1, 5, 1], (d) [3, 3, 1], (e) [1, 3, 1], and (f) [0, 2, 0]. (g) Overplotted S(Q = W, E) spectra, where colors correspond to the vertical colored lines in (c)–(f). Note how the INS data quality and information content vary significantly and somewhat unpredictably throughout the full 4D volume. All spectra were processed with Qdim1,step = 0.025 r.l.u, E = [0, 0.5, 40] meV and Qdim2,int = Qdim3,int = 0.05 r.l.u.

In addition to demonstrating the core functionality of the auto-slice maker, these examples highlight important challenges in the processing and interpretation of 4D S(Q, E) datasets. Namely, the measured 4D S(Q, E) often suffers from varying noise levels from different sources (Poisson counting statistics, electronics noise, background neutrons in instrument etc.), instrumental artifacts, blurring from Q–E-dependent 4D instrument resolution, and complex Q, E coverages from the scattering kinematic restrictions and instrument geometry. An illustrative example of such variation is evident in Fig. 5(c) versus Fig. 5(f): the former provides substantial coverage in the low- to mid-energy range with high data quality, while the latter shows irregularly shaped and sparse coverage in the high-energy region with a significantly degraded signal-to-noise ratio. Additionally, for phonon measurements, the observed intensity is modulated by the phonon polarization vectors. For example, the lowest transverse acoustic mode along K–Γ is absent in Fig. 5(d) due to polarization effects, but it is clearly visible in Fig. 5(c). Together, these factors can obscure valuable information contained in the dataset, make it more difficult to pinpoint features of interest, and limit the utility and statistical power of any individual region of data.

2.3. Automated data prioritization and quality filtering

The two comprehensive analysis methods discussed in the following sections, 1D systematic cuts and 2D systematic slices with BZ folding, are designed to automatically explore and utilize more of the full 4D dataset. However, given the vast and complex nature of the 4D data volume, along with variations in signal intensity and quality, we developed ways to automatically locate, prioritize, evaluate and filter data. Two primary methods are implemented to accomplish this: (1) BZ coverage determination and prioritization, and (2) slice evaluation and filtering.

The BZ coverage prioritization process identifies regions within the measured data volume that have greater Q–E coverage and ranks them for further processing. Furthermore, if a specific energy range is known to contain relevant features, this method can restrict analysis to BZs meeting a coverage threshold within that particular range. A 2D projection of this process is shown in Fig. 6 for the Ge dataset measured on ARCS. As shown in Fig. 6(a), the data initially reside in the (often conventional) Q basis used in Mantid. To perform the prioritization, a coarse 4D slice is generated by reducing, slicing and normalizing the dataset using Q projection axes derived from the inverse of the primitive-to-Mantid transformation matrix, ensuring each integer H, K, L point corresponds to a BZ center. By default, the Q extents span [−10, 10] in each primitive reciprocal lattice direction, encompassing up to 8000 potential BZs, though this range can be adjusted according to the crystal structure of the sample or the dynamical range of the instrument. A step size of 0.5 r.l.u. with a 0.25 r.l.u. offset ensures that each 1 r.l.u. cube is centered on a reciprocal-lattice node and subdivided into eight Q octants, as illustrated in Fig. 6(b). Energy binning is derived from the input file, using five bins across the full range (20% step size), yielding 40 Q–E bins per BZ. After the slice has been generated via the Mantid Python API, the fractional coverage for each BZ is computed as the ratio of valid bins (non-zero, non-NaN) to total bins. Fig. 6(b) shows the primitive-basis slice with BZ subdivisions, where green indicates signal and red indicates no data. BZs are then prioritized on the basis of these fractional coverage values [Fig. 6(c)], and those below a specified threshold are excluded to reduce processing time and avoid low-quality regions. This automated prioritization focuses computational resources on the most data-rich and experimentally relevant areas. Alternatively, specific BZs can be manually listed in the input file for targeted analysis.

Figure 6 2D projection of the BZ coverage determination and prioritization process applied to a Ge single-crystal dataset measured with the ARCS spectrometer. (a) Elastic Q map in the original conventional cubic Q basis. (b) Elastic Q map in the primitive Q basis such that each integer HKL is a BZ center. Each BZ (black lines) is partitioned into Q octants (gray lines), each with five energy bins, and the bin-wise Q, E coverage is determined as indicated by the red or green color. (c) Fractional Q, E coverages for each BZ within the dataset, ranked for further processing. The gray dashed line indicates the coverage threshold below which BZs are excluded (0.89 in this example).

To enhance the quality and utility of processed output data, pathSQE also includes automated, customizable slice evaluation and filtering capabilities. For example, when BZ folding, it is potentially useful to evaluate whether each slice should be included in the folding rather than incorporating all individual slices, as some may degrade the overall quality of the output. Similarly, for the 1D systematic cuts, when many potential Q points are explored, emphasizing those that meet certain criteria while deprioritizing others can improve both accuracy and efficiency. Filters may be designed as exclusive (e.g. removing noisy slices) or inclusive (e.g. retaining slices containing specific excitations). After each 1D or 2D slice is generated, it is evaluated against filter conditions specified in the input file, and only slices meeting all criteria are retained. The simplest predefined filter ensures adequate Q–E coverage by calculating the slice's fractional coverage and comparing it with a threshold. Beyond built-in filters, users can easily integrate custom ones; any Python function that takes a slice array as input and returns a Boolean can be used. This flexibility allows for tailored filtering to exclude distortions or artifacts or highlight features of interest. For example, in the FeSi ARCS dataset described later, a simple custom filter was used to automatically remove slices affected by a systematic experimental artifact.

2.4. 2D systematic slices and Brillouin zone folding

2.4.1. Procedural overview

pathSQE automates the slicing, normalization and visualization of all symmetrically equivalent 2D inelastic S(Q, E) slices within a time-of-flight INS dataset and additionally supports BZ folding as a complementary step that further integrates crystal symmetry. The BZ folding process consists of two key steps, a schematic example of which is shown in Fig. 3 using a simple 2D cubic crystal. BZ stacking exploits the lattice translational symmetries by shifting measurements from each BZ into the first BZ, consolidating the data that would otherwise span hundreds or thousands of BZs. Point-group folding, on the other hand, applies all point-group symmetry operations such that the data effectively occupy the irreducible wedge of the BZ to further symmetrize the data and improve statistics. In practice, pathSQE reverses the order of operations: point-group folding is performed first within each individual BZ, followed by BZ stacking. This ordering groups slices by BZ, yielding a more intuitive organizational structure for the processing and for visualization outputs (i.e. BZ Reports). Additionally, folding is applied to individual 2D slices, not on the full 4D volume, ensuring that all Q directions and symmetry operations align naturally with the slicing basis used for histogramming and normalization. Despite these modifications, the core symmetry principles underlying BZ folding remain unchanged. Because S(Q, E) can vary significantly across reciprocal space – e.g. with phonon scattering dominating at large |Q| and magnetic scattering at low |Q| – selectively folding regions with similar signal characteristics can help highlight specific features of interest and reduce the risk of signal dilution.

Upon execution, pathSQE identifies relevant BZs using the BZ coverage and prioritization algorithm (described previously). It retrieves the sample's point-group symmetry operations via the Mantid space-group factory (Arnold et al., 2014) and transforms the start and end Q points of each 2D path segment under all symmetry operations. The unique resulting Q point pairs are identified as symmetrically equivalent Q segments and are passed to the auto-slice maker, which generates the corresponding data slices. Optionally, the simulated harmonic phonon S(Q, E) counterparts are also produced (see Appendix A2). After slice generation, slice evaluation (if enabled) is applied to only retain slices satisfying all filter conditions. Crucially, within each BZ, all selected slices for a given path segment are stored as separate unnormalized data and normalization arrays, which are summed independently (Michels-Clark et al., 2016). For visualization purposes, the generated BZ Report file includes both the individual slices and a preliminary folded dataset, where normalization is performed at the BZ level. This visualization gives an in-depth view of the constituent data before folding and allows individual slices to be compared. Additionally, all slices are individually saved in a NeXus format using intuitive searchable filenames. The processing then shifts to the next BZ in the priority ranking. To ensure appropriate statistical weighting across the full dataset, normalization is only finalized after aggregating all unnormalized data and normalization arrays across all relevant BZs (Michels-Clark et al., 2016). Finally, the fully folded data are saved and visualized, and the corresponding dynamic susceptibility χ′′(Q, E) is also computed and saved.

2.4.2. Demonstration on Ge data measured with ARCS

The 2D systematic slicing and BZ folding workflow was applied to the Ge single-crystal dataset measured on ARCS, as illustrated in Figs. 7(a)–7(d). The selection of BZs for processing was informed by the BZ coverage determination and prioritization algorithm output, shown in Fig. 6(c). During intra-BZ folding (i.e. at the BZ Report level), symmetrically equivalent slices are processed; an example is shown in Fig. 7(a), which displays eight equivalent K–Γ slices in the [3, 1, 1] BZ. The folded experimental data in the same BZ along the full q path are shown in Fig. 7(b). Note that this particular BZ ranked ninth in the prioritization. The fully folded results, incorporating both point-group symmetry and BZ stacking, are presented in Fig.7(d), alongside the corresponding simulation in Fig. 7(c). Owing to the significant enhancement in data quality, the complete underlying dispersion relations become clearly resolved in the folded dataset, effectively highlighting the key excitations captured in the INS measurements.

Figure 7 Outputs from applying pathSQE 2D systematic slicing and BZ folding to a Ge single-crystal dataset measured on the ARCS spectrometer. (a) Excerpt from the [3, 1, 1] BZ Report, showing eight symmetrically equivalent slices along K–Γ segments. (b) Point-group-folded data in BZ [3, 1, 1]. Fully folded (c) simulation and (d) INS data along the high-symmetry q path. (e) Thin S(Q, E) shell from a one-minute measurement at fixed crystal orientation, with the out-of-plane [0, 0, L] direction integrated between −0.1 and 0.1 r.l.u. (f) Fully folded data from (e), projected along the high-symmetry q path. (g) The BZ is shown in the top left with the irreducible wedge outlined in black, a single facet of the wedge highlighted in gray and several high-symmetry q points marked in red. The schematic diagram illustrates 15 user-defined lower-symmetry q segments rastering the gray facet, indicated by the horizontal lines. (h) Fully folded data processed along each of the 15 q segments. The pink and blue boxes and corresponding q segments in (g) help indicate where on the rastered facet each segment resides. All spectra processed by pathSQE used Qdim1,step = 0.025 r.l.u. and E = [0, 0.5, 40] meV, while out-of-plane integration widths were set as Qdim2,int = Qdim3,int with values of (a–d) 0.05 r.l.u., (f) 0.1 r.l.u. and (h) 0.025 r.l.u.

We further demonstrate the utility of the workflow by applying it to two less conventional cases that were previously inaccessible using standard INS tools, as shown in Figs. 7(e)–7(h). As a first example, we start with extremely limited data from a single one-minute measurement at a fixed crystal orientation, as shown in Fig. 7(e). Using traditional manual INS data processing tools, extracting meaningful information from such a thin and complex S(Q, E) shell would be highly impractical. However, by applying our workflow, we obtain the fully folded result shown in Fig. 7(f), where a remarkable amount of information is recovered and aggregated from the sparse input. Importantly, pathSQE is also not restricted to high-symmetry q directions – arbitrary q segments can also be defined and processed equivalently. For instance, Fig. 7(g) shows 15 user-defined segments composed of low-symmetry directions that raster a facet of the irreducible wedge of the BZ. Folding this data using pathSQE yields the results in Fig. 7(h), demonstrating how using these often-underutilized lower-symmetry directions can enable effective extraction of the full dispersion surface. Together, the combination of systematic 2D slicing, BZ folding and detailed BZ Reports enables rapid and comprehensive exploration of complex 4D INS datasets, facilitating full characterization with unprecedented detail.

2.5. 1D systematic cuts

Similarly to the symmetry-aware 2D slicing capabilities described above, pathSQE also enables automated high-throughput extraction of 1D S(E) cuts across many symmetry-equivalent Q points. These 1D cuts are particularly useful for quantitative analysis of excitation energies, linewidths and spectral features, and are well suited for comparing measurements across experimental conditions (e.g. temperature or magnetic field) and simulations. Given an input Q_in, the workflow applies all point-group operations, maps the resulting vectors into relevant BZs using the prioritization output and identifies all unique Q points for processing. For each of these, a 1D S(E) cut is generated using the auto-slice maker algorithm. If multiple datasets under varied conditions are specified, analogous cuts are created for each in sequence. When enabled, harmonic phonon S(Q, E) simulations are also generated at each Q point using the automated simulation pipeline (see Appendix A2). All resulting data and simulations are saved and plotted to enable rapid visual comparison, and a summary scatter plot is produced to show the full set of probed Q points. This workflow was applied to the Ge dataset collected on ARCS at 10 K, using the K point (q = [0.75, 0.75, 0] in conventional coordinates) as the input q. The resulting scatter plot of all 196 probed Q_K points is shown in Fig. 8(a), and an excerpt of the generated Report containing the saved S(E) spectra is shown in Fig. 8(b). In addition to capturing variations in data quality and Q–E coverage, this approach facilitates the rapid identification of regions with prominent peaks and strong statistics for more detailed analysis.

Figure 8 1D systematic cuts applied to a Ge single-crystal dataset measured with the ARCS spectrometer. (a) 3D view of the locations of the 196 QK points identified within the INS data volume and processed by the workflow. (b) Excerpt from the K-point Report showing a subset of eight automated S(E) cuts at the listed QK points. All 1D spectra were processed with E = [0, 0.5, 40] meV, Qdim1,int = 0.025 r.l.u. and Qdim2,int = Qdim3,int = 0.05 r.l.u.

3. Additional examples

In the following subsections, we apply pathSQE to additional materials measured on various instruments at the SNS. Example input files for all of the cases shown herein – as well as an additional publicly available Si single-crystal dataset measured on ARCS – are provided in the pathSQE GitHub repository.

3.1. Folding, filtering and simulation comparison for FeSi measured with ARCS

The transition metal monosilicide FeSi exhibits rich physics due to the interplay of its lattice, electronic and magnetic degrees of freedom (Khan et al., 2022). Over the past several decades, it has been extensively studied for its thermoelectric potential, anomalous thermal properties and phonon softening, strong electron–phonon coupling, and topological phonon modes, among other intriguing behaviors (Jaccarino et al., 1967; Delaire et al., 2013; Miao et al., 2018; Khan et al., 2022). FeSi crystallizes in a noncentrosymmetric P2₁3 structure (space group No. 198), as shown in Fig. 9(a), with a simple cubic BZ and standardized high-symmetry q points, depicted in Fig. 9(b). Because of its eight-atom unit cell, FeSi possesses a complex phonon dispersion with 24 branches featuring closely spaced modes, degeneracies and crossings. In a previous study, Delaire et al. (2011b) collected an INS dataset on a large high-quality FeSi single crystal using ARCS at the SNS, measuring at 10 K with incident energy E_i = 40 meV and with 1° rotational steps spanning 35°. Complete experimental details are provided by Delaire et al. (2011b). In-plane and out-of-plane elastic scattering maps are shown in Figs. 9(c) and 9(d), respectively. Given its intricate phonon behavior, FeSi serves as an ideal test case for demonstrating the use of pathSQE to efficiently extract complex phonon dispersions from INS data.

Figure 9 pathSQE application to a FeSi single-crystal dataset measured with the ARCS spectrometer. (a) Unit cell and (b) BZ of FeSi. (c) In-plane and (d) out-of-plane elastic maps measured on ARCS at 10 K with Ei = 40 meV. pathSQE BZ-folded (e) data and (f) simulation across all 96 BZs above the 0.1 fractional coverage threshold. An experimental artifact is circled in red. (g) Running pathSQE with a custom user-defined slice evaluation and filter function automatically flags affected slices during processing and (h) systematically excludes them from the BZ folding. Inelastic spectra were processed with Qdim1,step = 0.025 r.l.u, E = [0, 0.25, 40] meV and Qdim2,int = Qdim3,int = 0.05 r.l.u.

The pathSQE BZ folding algorithm was applied to the dataset along the high-symmetry q path shown in Fig. 9(b). The fully folded experimental data and their equivalent simulation are displayed in Figs. 9(e) and 9(f), respectively, showing remarkable overall agreement. Notably, BZ folding significantly enhances measurement statistics, even clearly resolving fine features such as the topological acoustic and optical phonons at the M and R points, respectively (Miao et al., 2018). However, as highlighted in the red ellipses, an unexpected sharp signal appears along Γ–M at low energies (E = [0, 5] meV near Γ). Further investigation reveals that this is a systematic experimental artifact rather than a physical feature. Specifically, a saturation effect in the detector electronics, caused by strong Bragg peaks near the ends of the detectors, produces artificial inelastic tails. This artifact is especially pronounced in the vertical direction because of the geometry of the detector tubes. This effect is seen across all measured (T, E_i), examples of which are shown in Fig. 13 (see Appendix A3). In the inelastic data, this artifact manifests as spurious intensity extending along Γ–M due to the HHL alignment of the sample.

To circumvent this spurious scattering intensity, we implemented a custom automated filter to systematically identify and exclude affected slices from the pathSQE BZ folding. Specifically, the simple Python filter function flags a slice as contaminated if more than two pixels within  r.l.u. (i.e. the of the path segment closer to the M point) and E = [3, 6] meV exceed the 98th percentile intensity of the slice. This user-defined criterion is applied during the BZ folding procedure to automatically evaluate each slice and determine whether it should be included. The impact of this filter is illustrated in the BZ Reports, where affected slices are successfully identified and flagged, as shown in Fig. 9(g) alongside their corresponding simulations. The final fully folded dataset, shown in Fig. 9(h), demonstrates that the artifact has been effectively removed, yielding a more accurate experimental result that aligns well with the simulation in Fig. 9(f). This example demonstrates how pathSQE can improve both data quality and visualization, while also supporting systematic artifact detection and correction to enable more reliable interpretation of complex features in INS experiments.

3.2. Folding of mixed magnetic and nuclear scattering from MnO measured with HYSPEC

As the first system in which antiferromagnetism was experimentally observed (Shull et al., 1951), manganese(II) oxide (MnO) has long been investigated with neutron scattering. More recently, it has been studied for its intriguing structural, magnetic and transport properties (Pask et al., 2001; Ghosh, 2020). Below its Néel temperature, MnO adopts an f.c.c. rock salt structure (space group , No. 225), as illustrated in Fig. 10(a) (Roth, 1958). The corresponding BZ with annotated high-symmetry q points is shown in Fig. 10(b). An INS dataset was acquired on a MnO single crystal at 5 K using HYSPEC at the SNS (Winn et al., 2015), with an incident energy E_i = 25 meV. The instrument was operated in unpolarized mode in order to capture both nuclear and magnetic scattering contributions. Data were collected over a 280° angular range in 2° increments. In-plane and out-of-plane elastic scattering maps are shown in Figs. 10(c) and 10(d), respectively. Owing to the very limited vertical angular coverage of HYSPEC, the out-of-plane data span only 0.5 r.l.u. in total, as evidenced in Fig. 10(d). A representative radial 2D slice of S(Q, E) is shown in Fig. 10(e). Thus, this MnO dataset serves as an ideal test case for evaluating the performance of pathSQE under conditions of sparse Q, E coverage and mixed nuclear and magnetic scattering. Using pathSQE, the dataset was BZ folded along the high-symmetry q path defined in Fig. 10(b), yielding the fully folded spectrum shown in Fig. 10(f). In the folded data, both the acoustic phonon modes and magnon dispersions are clearly resolved. Beyond full BZ folding, pathSQE also enables the selective folding of specific BZs (e.g. only high- or low-Q regions) and supports custom slice evaluation and filtering, allowing nuclear and magnetic signals to be further isolated, enhanced and analyzed independently.

Figure 10 pathSQE application to a MnO single-crystal dataset measured with the HYSPEC spectrometer. (a) Unit cell and (b) BZ of MnO. (c) In-plane and (d) out-of-plane elastic maps measured on HYSPEC at 5 K with Ei = 25 meV. (e) Radial inelastic slice of S(Q, E) in the scattering plane. (f) pathSQE BZ-folded data across all 16 BZs above the 0.2 fractional coverage threshold. Inelastic spectra were processed using Qdim1,step = 0.025 r.l.u, E = [0, 0.2, 25] meV and Qdim2,int = Qdim3,int = 0.05 r.l.u.

3.3. Folding for lower-symmetry SnS measured with CNCS

In recent years, tin sulfide (SnS) has garnered significant attention for its performance in thermoelectric and optoelectronic applications (Norton et al., 2021). As a layered material featuring strong phonon anharmonicity, it is an intriguing system for studying lattice dynamics and phonon transport (Qin et al., 2016; Skelton et al., 2017; Lanigan-Atkins et al., 2020). Below the displacive phase transition temperature at ∼880 K, SnS crystallizes in the orthorhombic Pnma phase (space group No. 62), as shown in Fig. 11(a) (Chattopadhyay et al., 1986). Its corresponding BZ with labeled high-symmetry points is depicted in Fig. 11(b). In the low-temperature distorted phase, the SnS orthorhombic cell contains eight atoms and exhibits rather complex phonon dispersions with 24 branches featuring strong anisotropy (Lanigan-Atkins et al., 2020). Lanigan-Atkins et al. (2020) collected an INS dataset on a high-quality SnS single crystal using CNCS at the SNS (Ehlers et al., 2011), measured at 295 K with E_i = 17 meV and with 1° rotation steps about the [010] reciprocal direction spanning 90°. Lanigan-Atkins et al. (2020) give a more detailed discussion of the experimental details. In-plane and out-of-plane elastic scattering maps are provided in Figs. 11(c) and 11(d), respectively. Because of its more complex phonon dispersion, SnS presents an excellent test case for evaluating the ability of pathSQE to enhance information extraction from INS data. The SnS data were BZ folded using pathSQE along the high-symmetry q path shown in Fig. 11(b). The folded data from the highest-ranked individual BZ [−6, 0, 1] in the BZ prioritization step is displayed in Fig. 11(e). Finally, the fully folded data are shown in Fig. 11(f). Despite the phonon dispersion being quite complex and densely packed in the low-energy region, the underlying experimental dispersion is readily visible, enabling clearer and more comprehensive analysis of this challenging system.

Figure 11 pathSQE application to a SnS single-crystal dataset measured with the CNCS spectrometer. (a) Unit cell and (b) BZ of SnS. (c) In-plane and (d) out-of-plane elastic maps measured on CNCS at 295 K with Ei = 17 meV. pathSQE BZ-folded data (e) in top-ranked BZ [−6, 0, 1] and (f) across all 65 BZs above the 0.2 fractional coverage threshold. Inelastic spectra were processed using Qdim1,step = 0.025 r.l.u, E = [0, 0.2, 17] meV and Qdim2,int = Qdim3,int = 0.05 r.l.u.

4. Conclusions

Herein, we have presented the motivation for and implementation of pathSQE, an automated workflow for the processing and analysis of single-crystal INS datasets. The workflow addresses several bottlenecks and limitations of current INS tools, significantly accelerating data handling and unlocking the full potential of modern 4D time-of-flight INS datasets. In addition to streamlined 1D and 2D slicing core functionality, pathSQE seamlessly incorporates advanced features – such as intelligent data prioritization, automated slice evaluation and filtering, and BZ folding – that extend beyond the current state of the art. The workflow was successfully demonstrated on a range of materials spanning diverse symmetries and physical phenomena measured across multiple DGS instruments at the SNS. In all cases, pathSQE dramatically reduced the required manual effort, significantly enhanced the comprehensiveness of the analysis and enabled unprecedented exploration of the full INS datasets.

Beyond improving individual dataset analysis, pathSQE is designed to support broader advancements in INS. While pathSQE primarily focuses on data processing, its framework also facilitates direct comparison between experimental data and theoretical models. The inelastic neutron scattering cross section and, therefore, S(Q, E) are not translationally invariant in Q space. However, their support is (for crystalline systems, as assumed here); the dispersions (though not the intensities) periodically follow the tiling of reciprocal space by BZs and conform to point-group symmetries within each BZ. By applying the same folding operations to both the measured data and the simulated or fitted models, one obtains a consistent representation and comparison, enabling accurate model validation and refinement within the folded data space where experimental statistics are improved. The increasing use of machine learning in post-experimental analysis necessitates faster, more automated processing pipelines to generate large ML training datasets (Chen et al., 2021). Furthermore, real-time experimental steering requires immediate feedback (Doucet et al., 2021), a demand that pathSQE helps meet by reducing the time needed to process relevant data. Similarly, in order to perform on-the-fly experimental optimization using ML-driven or other algorithmic agents, an analytical tool, such as the auto-slice maker algorithm, is needed to rapidly process the data with minimal human interaction. Finally, with the advent of more advanced modeling capabilities, there has been growing interest in performing full virtual INS experiments using instrument digital twins (Islam et al., 2019; Lin et al., 2022; Han et al., 2025), which would greatly benefit from rapid parallel experimental processing and simulations to synchronously validate model accuracy. These considerations are especially relevant as new instruments with even greater data acquisition capabilities are being designed and constructed at facilities around the world, such as the Second Target Station at ORNL (Crawford et al., 2010) and the European Spallation Source (Andersen et al., 2020). Taken together, these developments underscore the growing need for robust scalable analysis workflows like pathSQE, capable of autonomously processing and extracting meaningful insights from increasingly information-rich experimental INS datasets.

APPENDIX A

A1. The input file

The input file defines how users interact with and control various aspects of the pathSQE program. The input file follows an intuitive, flexible and consistent format, structured into five main blocks. An example input file, which was used to process the Ge ARCS data shown in many of the figures in this work, can be seen in Fig. 12. In addition to this input file, a define_data.py file, which contains important details such as the INS raw event data file locations and the experimental UB matrix of the sample, must exist in the working directory. Finally, for analysis including automated phonon S(Q, E) simulations, a POSCAR crystal structure file and FORCE_CONSTANTS or FORCE_SETS file must be present.

Figure 12 pathSQE is controlled using an input file with five main blocks: (a) experimental and Q point information, (b) per-slice Q, E binning and integration ranges, (c) data processing and analysis options, (d) simulation settings, and (e) saving and output preferences. The file shown was used to process the data displayed in Figs. 6, 7(a)–7(d) and 8.

A2. Automated simulation pipeline

To interpret the excitations probed in an experiment, it is often useful to compute analogous S(Q, E) spectra for comparison. For meaningful results, simulations must account for key experimental parameters such as instrument resolution, energy range and Q-space sampling – requirements that typically involve significant setup time and computational expertise. To address this, pathSQE provides an integrated simulation pipeline that automates the generation of experimentally aligned harmonic phonon spectra using Phonopy. These simulations are designed to mirror the auto-slice maker output and support rapid reproducible comparisons without disrupting the broader workflow. Although currently limited to phonons via Phonopy, the pipeline is modular and can be extended to other tools such as SpinW.

The simulation pipeline consists of two main stages: (1) generating the base S(Q, E) spectrum, and (2) postprocessing to more closely match experimental conditions. As with slicing, the simulation workflow loops through all specified Q points or path segments automatically. Prior to this, the user must supply structural and force constant inputs (POSCAR and FORCE_CONSTANTS or FORCE_SETS in VASP + Phonopy format), which can be generated via standard density functional theory methods or obtained from online materials databases (Han et al., 2024; Choudhary et al., 2020; Jain et al., 2013). After the structural and force constant files have been loaded, phonon eigenvectors and mode energies are computed on a reciprocal mesh. Input Q points or path end points are transformed into the primitive reciprocal basis. For 2D simulations, a fine Q path is generated between Q_1,in and Q_2,in using steps four times smaller than Q_dim1,step. Temperature and sample composition are read from the input file, and isotope scattering lengths are retrieved from NIST data (Sears, 1992). Using these parameters – Q, temperature, eigenvectors and scattering lengths – the dynamical structure factor is computed, yielding phonon mode energies and their associated intensities.

Postprocessing then adapts the raw simulation output to resemble actual measurements. For 1D spectra, intensities are binned on a fine energy grid. In the 2D case, a high-resolution Q–E histogram is constructed, using bin sizes four times finer than the experimental resolution. Energy resolution is mimicked using Gaussian broadening along E, with FWHM estimated from elastic line width using the Windsor approximation. A fixed-resolution option is also available for computational efficiency. In two dimensions, Gaussian broadening is also applied in Q, followed by coarse rebinning to match experimental bins. Optionally, a mask reflecting the experimental Q–E coverage is applied to the simulation, limiting the final result to experimentally accessible regions. For 1D spectra, similar masking can restrict the energy range. After these steps, the processed experiment-matched S(E) or S(Q, E) spectrum is returned.

A3. FeSi Bragg peak tail artifacts

Fig. 13 illustrates systematic Bragg peak tail artifacts in several single-crystal FeSi datasets measured with the ARCS spectrometer.

Figure 13 (a) Inelastic slice of FeSi data with T = 300 K and Ei = 70 meV along [H, H, 6] including strong, slightly inelastic features at half-integer (i.e. M) points. (b) An in-plane Q map integrated over E = [1, 5] meV, illustrating that these signals appear systematically at M points throughout the full Q volume. (c) An out-of-plane Q map integrated over E = [0, 5], revealing the origin of these spurious features; namely, Bragg peaks near the detector ends are dramatically smeared and, in this case, reach the parallel Q plane below. (d–f) Further investigation indicates that similar artifacts appear across all measured (Ei, T). For example, out-of-plane Q maps integrated over E = [0, 5] are shown for T = 10 K and Ei = 40 meV with (d) L = [−1.1, −0.9], (e) L = [−0.1, 0.1], (f) L = [0.9, 1.1]. In inelastic Γ–M slices along these Bragg peak tails, the artifacts manifest as extended, slightly inelastic signals, as shown in Figs. 9(e)–9(g).