A simulational study of the indirect-geometry neutron spectrometer BIFROST at the European Spallation Source, from neutron source position to detector position

The incident detector rates that are anticipated for the indirect-geometry cold-neutron spectrometer BIFROST at the European Spallation Source are estimated, and the use of powerful simulation tools for the correct interpretation of neutron transport in crystalline materials is demonstrated.

The European Spallation Source (ESS) is intended to become the most powerful spallation neutron source in the world and the flagship of neutron science in upcoming decades. The exceptionally high neutron flux will provide unique opportunities for scientific experiments but also set high requirements for the detectors. One of the most challenging aspects is the rate capability and in particular the peak instantaneous rate capability, i.e. the number of neutrons hitting the detector per channel or cm 2 at the peak of the neutron pulse. The primary purpose of this paper is to estimate the incident rates that are anticipated for the BIFROST instrument planned for ESS, and also to demonstrate the use of powerful simulation tools for the correct interpretation of neutron transport in crystalline materials. A full simulation model of the instrument from source to detector position, implemented with the use of multiple simulation software packages, is presented. For a single detector tube, instantaneous incident rates with a maximum of 1.7 GHz for a Bragg peak from a single crystal and 0.3 MHz for a vanadium sample are found. This paper also includes the first application of a new pyrolytic graphite model and a comparison of different simulation tools to highlight their strengths and weaknesses.

Introduction
The European Spallation Source (ESS) ERIC (Peggs et al., 2013;Garoby et al., 2017;Andersen et al., 2020) is designed to operate using the most powerful spallation neutron source in the world and to provide unprecedentedly high neutron fluxes for instruments of various neutron techniques. One of these instruments is BIFROST (Ronnow et al., 2014;Freeman et al., 2015), a high-flux indirect-geometry cold spectrometer, optimized for small samples and extreme environments. BIFROST is primarily intended for single-crystal inelastic scattering studies, providing exceptional flux as it can operate in a white-beam mode. This flux allows for entirely new options for detailed investigations of complex multimode dynamics, hybrid modes, electro-magnons, spin wave continua and gap studies, under extreme conditions with controlled temperature, pressure and magnetic fields.
Harnessing the full ESS pulse by employing a polychromatic beam carries a high potential risk for detector rates that can saturate the detectors and therefore degrade the performance of the instrument. The chosen detector technology of ISSN 1600-5767 BIFROST, position-sensitive 3 He tubes are the 'gold standard' for neutron detection (Knoll, 2010). They are, however, quite rate limited. Non-position-sensitive tubes saturate at 100 kHz, but for position-sensitive 3 He tubes, operation at instantaneous rates above 30 kHz can be problematic. The exact rate capability of a detector is dependent on the details of the readout electronics. It is therefore essential to evaluate the rates anticipated for high-flux instruments (Kanaki, Klausz et al., 2018;Stefanescu et al., 2017;Piscitelli et al., 2017), in order to extract the respective detector requirements.
Monte Carlo simulation plays a key role in the development and characterization of instruments as a reliable, cheap and versatile tool . Feedback from simulations taken into account in the development of the instrument design can reduce the number of physical prototypes needed and also enables the quantification of otherwise unmeasurable properties. This is particularly the case for complicated instruments such as BIFROST. Development of complete and detailed instrument simulation models enables simulations from source to detectors, offering the opportunity to discover and decouple otherwise undetectable cumulative effects. These models can provide valuable input for developing calibration and correction routines for data reduction and analysis, and could later be used for experiment planning by users, to predict experimental conditions from specific proposed samples, i.e. sample size and composition.
To make detailed full-instrument simulations possible, advanced simulation tools have been developed, such as NCrystal , which enables Monte Carlo simulations of thermal neutrons in crystals, and Monte Carlo Particle List (MCPL)  https://mctools.github.io/mcpl/), which enables communication between different software packages. These tools can greatly enhance the capabilities of existing and widely used simulation software such as McStas (Lefmann & Nielsen, 1999;Willendrup et al., 2004) and Geant4 (Agostinelli et al., 2003;Allison et al., 2006Allison et al., , 2016, and enable implementation of full simulation models by connecting them to perform a chain of simulations, using each of them where they are the most capable. Such simulation of full instruments is a novel method that has been applied only in a handful of cases (Kanaki, Klausz et al., 2018;Stefanescu et al., 2019).
In this study, multiple Monte Carlo simulation tools are used together to implement a full simulation model of the BIFROST instrument from the neutron source to the detector position. This is the first application of the new NCrystal pyrolytic graphite material model. This full instrument model is used to estimate the incident detector rates that are anticipated in the case of the highest possible incident neutron intensity -which will be referred to as the 'highest case' -and normal-use scenarios, intended to serve as a basis for the determination of detector requirements for rate capability. For a part of the instrument both a McStas and a Geant4 model are implemented, facilitating the cross-validation of results and the comparison of the two simulation software packages.
In the following sections, the instrument and simulation model geometries and tools are introduced first, followed by the presentation of incident rates for elastic peaks for various instrument parameters, sample types, sizes and mosaicities, along with a demonstration of the differences between McStas and Geant4 simulation results. The study concludes with a demonstration of the elastic signal of a standard calibration sample.

BIFROST instrument and simulation model
BIFROST is a 162 m-long cold-neutron spectrometer intended to be built as a first-tranche instrument for ESS. It combines an indirect-geometry time-of-flight (ToF) front end and an angular-and energy-multiplexed crystal-analyserbased back end: a back end similar in principle to that installed recently at the CAMEA spectrometer at the Paul Scherrer Institut (Markó et al., 2018;. BIFROST is designed  to maximize the use of the ESS long pulse to enable measurements on small samples and study dynamic properties, transporting a 1.7 Å wavelength band to the sample and investigating an energy transfer range from À3 to +55 meV. The envisioned application fields include materials science, magnetism, life sciences and planetary sciences (Freeman et al., 2015).
The instrument consists of three main technical subsystems: the beam transport and conditioning system, the sample exposure system, and the scattering characterization system. A schematic model of the instrument is depicted in Fig. 1.
The beam transport and conditioning system is relatively simple. It has a curved guide section inside the bunker to lose the line of sight and four choppers as the only moving parts. Three of these choppers are placed inside the bunker and the fourth is placed in the middle of the instrument. The first one is a pulse-shaping chopper; this is the only chopper determining the energy resolution. The other three -two frame overlap and one bandwidth chopper -serve to sort out unwanted frames from the fast pulse-shaping chopper and to avoid pulse overlap at the sample position, respectively. The pulse-shaping chopper can reduce the ESS pulse width by a factor of up to 30 to match the best analyser resolution or allow the full pulse to reach the sample, which will result in a relaxed resolution but an order of 10 10 n s À1 cm À2 flux on sample. It is this mode which poses the greatest rate challenge for the detectors. The sample exposure system allows measurements with a strong magnetic field, high pressure and cryogenic temperatures. One of the main limitations today in single-crystal neutron spectroscopy is that measurements are only possible with large samples, which are not available for many sample types, but BIFROST will enable the study of sub-cubic millimetre samples thanks to its exceptional flux on sample and the efficient scattering characterization system. The scattering characterization system in Fig. 2 consists of the filtering system and the secondary spectrometer tank, covering a 90 scattering angle in the horizontal plane in two tank settings, in the 7-135 2 range. The filtering system, which is essential for background reduction on BIFROST, includes a cooled beryllium filter with roughly 90% transmission of neutrons with an energy below 5 meV (4.05 Å ) but very low transmission of neutrons with energies above, and coated lamellas as a radial collimator (Groitl, Rantsiou et al., 2016).
The secondary spectrometer tank houses multiple sets of analysers and detectors for different neutron energies, arranged in nine 'Q-channels'. Depending on the scattering angle of a neutron on the sample, it enters one of the Q-channels, which are separated by cross-talk shielding between them. In each Q-channel, several crystal analyser arcs are placed one behind the other to select different final energies by scattering neutrons vertically (down) towards the corresponding set of position-sensitive detectors, employing Rowland focusing [see Fig. 2(b)]. Further cross-talk shielding is applied to ensure that neutrons can reach the detectors only by scattering from the corresponding analyser arcs.
With this arrangement, BIFROST utilizes a variant of a novel analyser setup called CAMEA Markó et al., 2018), an acronym for continuous angle multiple energy analysis. Enabling multi-energy analysis in a single Q-channel by placing the analysers for higher neutron energies behind the ones for lower energies is possible owing to the high transparency of the 1 mm-thin highly oriented pyrolytic graphite blades (Mildner et al., 2001). The blades to be used have high mosaicity (60 arcmin), to allow the application of the prismatic analyser concept   The BIFROST scattering characterization system. 3D model of the secondary spectrometer tank (a) and side-view sketch of the sample and a single Qchannel with five energy channels in it (b).

Figure 1
Schematic model of BIFROST from source to detection position. From https://europeanspallationsource.se/instruments/bifrost. detector triplets for all five neutron energies chosen for BIFROST (2.7, 3.2, 3.8, 4.4 and 5.0 meV) in each Q-channel. According to the prismatic analyser concept, each of the three detectors of a triplet records a slightly different region of energy, as neutrons with different energies are scattered in slightly different directions.
In order to provide enough space for the detector tubes, the analysers and corresponding detectors in adjacent Q-channels are slightly shifted radially. The sample-analyser distances are increased or decreased by 4.6-7.5% in two out of three Q-channels. However, the analyser-detector distances are kept unchanged to keep the detectors of the same energy on the same vertical planes and hence keep the spectrometer tank geometry simple. As a result, the sample-analyser distance is shorter or longer than the analyser-detector distance in two out of three Q-channels, resulting in a slight asymmetry to the Rowland geometry. The three different types of Q-channels are repeated three times, giving the nine Q-channels a 'triple stagger' geometry.
The simulation of the BIFROST instrument is divided into two parts: the simulation of the long-beam transport and conditioning system from the moderator until the end of the last guide section before the sample, and the simulation of the sample and scattering characterization system together (see Fig. 3). The first part is done using McStas only; the second part is implemented and simulated in both McStas and Geant4, in order to compare the results of these simulation tools and to demonstrate why it is advantageous to use Geant4 for the back end of the instrument.
The transition between the two parts is facilitated by the MCPL tool. MCPL is a binary file format dedicated for storage and interchange of particles between various Monte Carlo simulation applications, like McStas, Geant4, McXtrace (Bergbä ck Knudsen et al., 2013) and MCNP (Werner, 2017). For the simulation of neutron transport in crystalline materials, NCrystal is used in both McStas and Geant4. In the next subsections all simulation tools and models are introduced.

McStas model
McStas (Lefmann & Nielsen, 1999;Willendrup et al., 2004) is a Monte Carlo simulation tool dedicated for simulation of neutron scattering instruments and experiments. It is user friendly, cross platform and open source, and uses a raytracing algorithm that enables fast neutron transport simula-tions over long distances and through many components, which is necessary for long instruments like BIFROST.
For the simulation of neutrons from the source to the end of the beam transport and conditioning system, a previously developed McStas model of the instrument is used (https:// bitbucket.org/europeanspallationsource/nosg-baselines/src/dev/ BIFROST/Optics/McStas/). This model, depicted in Fig. 4, contains the butterfly moderator source ('ESS_butterfly') (Andersen et al., 2018;Zanini et al., 2018), the four choppers, all guide sections and several McStas monitor components to characterize the beam at multiple locations along the guide. The source is used with the highest intensity, so the deduced rate numbers in this paper correspond to the maximum accelerator power of 5 MW. Expected rates will scale linearly with source power for constant proton energy. At the end of the last guide section all neutron data are saved in an MCPL file using the 'MCPL_output' McStas component. This file serves as input for both the McStas and the Geant4 simulation models of the second part of the instrument.
The McStas model of the sample and scattering characterization system (Fig. 4) contains a crystalline sample, one Q-channel including all five analyser arcs, and McStas monitor components at several places probing ToF, energy and position distribution of neutrons, in order to examine the change of the neutron beam. The analyser arcs consist of 7-9 blades using the 'NCrystal_sample' component with pyrolytic graphite material, described in more detail in Section 2.4. NCrystal is also used for all crystalline samples throughout this study. The simulation model does not contain the sample  Outline of simulation scheme.

Figure 4
McStas model of the beam transport and conditioning system (a) and the sample and scattering characterization system (b). The figures are at different scales. The neutron source-to-sample distance is 162 m and the sample-to-analyser distances are 1189-1622 mm. environment, eight out of the nine Q-channels, the filtering system, cross-talk shielding or the detectors.
Using a reduced geometry and excluding any model of the sample environment is intentional, aiming to get a conservative estimate in terms of the highest detector rates, but implementing only one Q-channel is the result of a limitation coming from the linearity of the McStas simulation process. In a McStas instrument definition file, the geometrical components like the source, guide sections, choppers, slits and sample are placed one after the other. McStas by default propagates neutrons from component to component in the exact order that they appear in this file. All neutrons that miss or do not interact with the component downstream are removed from further simulation. This process makes the simulation of long instruments fast, but on the other hand restricts the neutrons to follow one exact path, which does not allow the simulation of multiple Q-channels simultaneously. It is possible to change this behaviour by grouping components together, as described in the user and programmers' guide (Willendrup et al., 2020). This way it is possible to some extent to split the beam, by having a group of components as the potential next target of neutrons, but after the interaction with one of the group members, McStas tries to propagate all neutrons to the component that appears next in the instrument file. This means that if a user wants to split the beam and allow propagation in multiple directions through different components, consecutive groups must be implemented, all of which include the subsequent component in each direction. Extensive use of such grouping makes the instrument file immensely complex and still prohibits multiple interactions within one group, or back and forth propagation between groups of components.
For these reasons, only one Q-channel is implemented, within which this technique is used to handle the five sets (arcs) of analyser blades, all of which divide the beam into the partition that is scattered towards the detectors and the partition that is propagated towards the next set of analysers (or the beam stop behind the last set). In order to allow neutrons to proceed without interaction with a set of analyser blades, an extra virtual component is added to each group, which mimics an interaction without changing the neutron state and thereby prevents neutrons from being removed from the simulation. As mentioned, neutrons still have to follow the order of the groups, so backscattering or multiple scattering among blades of the same arc is still not possible in the simulations.
Although the cross-talk shielding between energy and Q-channels is not explicitly included in the model, as a consequence of the above-described process, a neutron can reach a particular detector tube only by scattering in one of the corresponding analyser blades. This is practically equivalent to an ideal cross-talk shielding absorbing all stray neutrons. The case is similar for the filtering system, which is replaced by a monitor component that transmits all neutrons below 7 meV energy and none above. As mentioned earlier, the transmission of beryllium drops sharply around 5 meV, which is in fact the highest of the five final energies selected by the analysers. Simulation of effects of this transition in the transmission is out of the scope of this article, and using ideal transmission in the 0-7 meV energy region keeps the rate estimates conservative.
The intent is to determine the incident detector rates, and therefore the simulation of the detection process is also out of scope. Detectors are modelled with McStas monitor components, and a neutron is counted as incident for a detector tube if it crosses the plane at the centre of the detectors within the outline of that particular tube. The sample-analyser distance is equal to the analyser-detector distance, meaning that symmetrical Q-channels are modelled.
As this subsection demonstrates, using McStas to model such a complex system as the analyser-detector system of BIFROST in detail is not straightforward and is subject to certain limitations.
McStas version 2.5 is used for the simulations.

Geant4 model
Geant4 ( Geant4 model of the sample and scattering characterization system including one Q-channel (a) or all Q-channels (b). The blue lines indicate simulated neutron paths. developed at CERN with applications in many fields, e.g. highenergy physics, nuclear physics, accelerators and medical physics. Its usability for simulation of neutron detectors has been greatly improved by the ESS Detector Group by building a framework Kittelmann et al., 2014) around it which adds several functionalities and integrates NCrystal and MCPL.
The Geant4 simulation model of the sample and scattering characterization system, depicted in Fig. 5, contains the same parts as the McStas model (crystalline sample, all five analyser arcs in a Q-channel) but with the option to simulate with all Q-channels included. To make the results comparable with the McStas model, the McStas monitor components are mimicked with empty volumes with exactly the same location and surface, in order to create histograms with the same predefined spatial, energy and ToF resolution.
Although the cross-talk shielding could be easily implemented in the model, for the same comparison purposes it is replaced by certain conditions applied to data at the analysis level. This means that neutrons cannot skip parts of the model. They can, however, scatter back and forth between the geometrical components many times, unlike in McStas. This gives the possibility to analyse effects of cross-talk on signals and to evaluate shielding strategies.
The Geant4 physics list used is QGSP_BIC_HP_EMZ, which uses high-precision models and cross sections for neutron energies lower than 20 MeV, and allows the correct treatment of thermal and cold neutrons when combined with NCrystal.

NCrystal
NCrystal  is a novel open-source software package for modelling thermal neutron transport in crystalline materials. It consists of a data library and associated tools which enable calculations for Monte Carlo simulations. It can be used together with McStas and Geant4 to enhance their capabilities for the correct treatment of neutron transport in typical components of neutron instruments, including beam filters, monochromators, analysers, samples and detectors. Physics modelled by NCrystal includes both coherent elastic (Bragg diffraction) and incoherent or inelastic (phonon) scattering. It treats all valid Bragg diffraction on each reflection plane explicitly and is able to use various models for inelastic scattering on phonons. Its data library (https:// github.com/mctools/ncrystal/wiki/Data-library) already contains the most popular crystals and the results are validated against the EXFOR database (https://www-nds.iaea.org/exfor/) and existing crystallographic software.
NCrystal focuses initially on scattering in single crystals or polycrystalline materials and powders. Most single-crystalline materials are appropriately modelled with crystallites orientated around some reference orientation with a Gaussian distribution that has a standard deviation of the mosaicity of the crystal. There are, however, single-crystalline materials with crystallite distributions so different from Gaussian that this approximation does not hold. One of these materials is pyrolytic graphite, which is widely used as a monochromator and analyser in neutron instruments. This is the case for BIFROST, where 369 highly oriented pyrolytic graphite analyser blades are used altogether in the nine Q-channels. Graphite has a layered structure, made up of graphene sheets in which carbon atoms are arranged in a hexagonal structure. In highly oriented pyrolytic graphite the crystallite axes orthogonal to the graphene sheets are distributed along a preferred direction, suitable for description with a Gaussian mosaicity distribution, but the orientation around this axis is completely random, resulting in powder-like features in neutron scattering. Recent developments have enabled NCrystal to handle materials with this kind of structure by using a specialized model for layered crystals.
NCrystal is used for the crystalline materials in both the McStas and the Geant4 models of BIFROST. The samples are modelled as single crystals, but for the pyrolytic graphite analysers the new layered crystal model is used. Demonstrating some of the NCrystal pyrolytic graphite properties, the components of total cross sections and the distribution of randomly sampled scattering angles are depicted in Fig. 6 for a powder sample, not the layered crystal distribution. As presented in the NCrystal data library, the cross sections are validated against experimental data (Walton et al., 1960). Until the latest version (1.0.0), NCrystal has treated absorption with the simple model of absorption cross sections being inversely proportional to the neutron speed. The absorption cross section for a particular neutron speed is calculated by scaling the value given at the reference speed of 2200 m s À1 . This applies for NCrystal materials used in McStas simulations, but not for those in Geant4 Components of the total interaction cross section (a) and two-dimensional scatter plot (b) for NCrystal pyrolytic graphite powder.
simulations. Geant4 models secondary particles produced in absorption, and therefore the NCrystal plugin does not interfere with the Geant4 absorption physics at all. As an example, the minor differences in absorption cross section for pyrolytic graphite are shown in Fig. 7.
All results of simulations have systematic and statistical uncertainties. Unlike systematic uncertainties, which are mainly due to imperfect modelling of the system, statistical uncertainties can be reduced by increasing the number of simulated particles. Throughout this work where the uncertainty is not significant, it is not indicated. Sources of systematic uncertainties are considered and are not expected to change the conclusions.
With the simulation models and tools described in this section, a detailed analysis is carried out, in order to give an estimate of the incident rates that are anticipated for detectors at BIFROST. The incident detector rates for elastic peaks using various samples and instrument parameters are presented in the next section, as is a comparison of results with McStas and Geant4 simulations.

Incident rates for coherent elastic peaks
Determination of anticipated detector rates for an instrument is a key part of defining requirements for the detectors to be used. It can prevent the detector rate capability from becoming the bottleneck of experiments or a source of performance degradation. For this reason, the simulation tools and models described in the previous section are used to determine the highest time-averaged and the highest instantaneous (peak) incident rates for the detector tubes. These rates are determined for the highest case and for more realistic operational conditions.
The highest-case incident rate for a single detector tube occurs when a strong Bragg peak from a single-crystal sample gets reflected to the detector. To get neutrons scattered on the sample onto the detectors, their energy has to match one of the energies selected by the analyser arcs. For the ESS source spectrum, the guide transmission and the energy resolution of the analysers, the highest incident rates are expected for the 5.0 meV (4.045 Å ) neutrons.
Regarding highly reflective materials that would result in the highest detector rates, the truly highest-case sample would be a pyrolytic graphite single crystal (d 002 = 3.3555 Å ), but to obtain results from a less unrealistic sample with large enough lattice parameter and strong Bragg peak, simulations are also performed using an yttrium oxide (Y 2 O 3 ) single crystal (d 2 2 2 2 2 = 3.0724 Å ).
Further parameters that influence the rates on detectors are the pulse-shaping chopper opening time, sample size and sample mosaicity. To get the highest possible rates, maximum flux mode is used, where the pulse-shaping chopper is fully open, resulting in a 10 10 n s À1 cm À2 flux on sample. The instrument is designed to facilitate measurements on small samples, but dimensions of up to 1.5 cm are possible. Therefore, cylindrical samples with a diameter and height of 1.5 cm are used in the simulations. The sample mosaicity resulting in the highest rate can depend on the divergence of the incident beam on the sample. However, as a rule of thumb the highest rates are expected when the mosaicity of the sample matches that of the analysers, so the sample mosaicity is set to 60 arcmin. It is shown later in Section 3.2 that this is a good assumption, and within the 60 AE 20 arcmin sample mosaicity range it has a less than 10% effect on the incident detector rates.
In order to realize the simulations, the samples are oriented to fulfil the Bragg condition for the incoming 5 meV neutrons on the selected scattering planes, and the single Q-channel modelled is rotated according to the resulting scattering angle. In BIFROST the whole scattering characterization system can be rotated around the sample, so having a Bragg peak from any sample in the exact direction of a single Q-channel is perfectly realistic.
The simulation with the parameters described above is done in two steps. First, the neutron transport from the source to the end of the beam transport and conditioning system is considered with the McStas model, saving neutron data at the end in an MCPL file. This file is then used as the source term for the simulation of the sample and scattering characterization system, which is done with both McStas and Geant4. The results are presented in Section 3.1. In the subsequent subsection the impact of different parameters like sample and analyser mosaicity, sample size and pulse-shaping chopper opening time is investigated; because of the good agreement of the McStas and Geant4 simulation results (demonstrated later in Section 3.1), this is only done using the Geant4 model. Comparison of neutron interaction cross sections of pyrolytic graphite powder in Geant4 and NCrystal. The lines for Geant4 neutron capture and NCrystal absorption cross sections are barely distinguishable, the former being higher in the whole energy range depicted.
Neutrons of a broad energy range -centred around 5.0 meV -are scattered on the sample toward the Q-channel, reaching the 2.7 meV analysers and therefore missing from the sample transmission spectrum. The change of the spectrum between the 2.7 and 5.0 meV analysers is caused by the spread of neutrons and absorption in the analysers. The neutrons selected by the 5.0 meV analysers are scattered towards the corresponding detector triplet and are therefore absent from the analyser transmission spectrum.
The time-averaged neutron intensities acquired by the integration of the energy spectra are shown in Table 1. The results of the McStas and Geant4 simulations agree with only minor differences.
The structure of the energy spectrum of the detector triplet is the result of summing the spectra of all three detectors, as depicted in Fig. 9. As expected, the analysers scatter neutrons with slightly different energies in slightly different directions -in accordance with the Bragg law -and as a result of this vertical spread, the three detectors of the triplet record slightly different regions of energy, in accordance with the prismatic analyser concept. The sample spreads the neutron beam similarly but in the horizontal plane. The combined effect of these processes is visible in Fig. 10, showing a diagonal shape in the time-averaged neutron intensities in the plane over the detector tubes.
Integrating the incident neutron intensities in Fig. 10 over the areas of the tubes gives the time-averaged incident rates for the tubes, presented in Table 2. The time-averaged incident rate for a single detector tube can be almost as high as 70 MHz.
Given that BIFROST is a ToF instrument at a spallation neutron  Table 1 Time-averaged neutron intensities at the sample and different parts of the scattering characterization system with a pyrolytic graphite single-crystal sample simulated with McStas and Geant4.

Figure 9
Time-averaged incident neutron energy spectra of the 5 meV detector tubes in McStas (a) and Geant4 (b) simulations with the pyrolytic graphite sample. The lines are only joining the points. source, the incident detector rate has a pulsed time structure. The ToF distribution of a single pulse on the 5 meV detectors is depicted in Fig. 11. By taking into account only those neutrons which arrive at the detectors at the peak of their ToF distribution in a short time range of 0.1 ms, the peak incident rates presented in Table 3 are acquired. Owing to the distinct energy range and therefore different ToF spectra of the tubes, the highest peak incident rate occurs at different times for each tube of a triplet. The results demonstrate that the peak incident rate on a single detector tube can be as high as 1.7 GHz. 3.1.2. Yttrium oxide sample. To give an impression of how the rates change with a different single crystal that is not the highest-case sample but also has a strong Bragg peak, the same simulation and analysis process is repeated using a Y 2 O 3 sample. Fig. 12 demonstrates the change of the timeaveraged neutron energy spectrum along the neutron's path at the sample in the scattering characterization system. The time-averaged neutron intensities acquired by the integration of the energy spectra are shown in Table 4.
The results of the McStas and Geant4 simulations agree with only minor differences, the only exception being the transmission spectrum of the sample. The transmission of the sample is 10% higher in McStas, with the same lowerintensity bands apparent in the spectrum caused by several crystal planes where the Bragg criterion is fulfilled for different neutron energies -including the hkl = 2 2 2 2 2 2 plane for 5.0 meV. The source of the discrepancy is the absorption process in the sample because different absorption cross sections are used in the two simulation tools, as described in Section 2.4.
The spectrum of the 2.7 meV analyser arc shows that, despite the presence of multiple strong Bragg peaks, it is only the 5.0 meV neutrons that are scattered toward the Q-channel. There are, however, three scattering planes (hkl = 2 3 3 3 3, hkl = 1 2 2 3 3 and hkl = 1  Table 2 Time-averaged incident neutron rates of the 5 meV detector tubes in McStas and Geant4 simulations, with a pyrolytic graphite sample.

Figure 11
ToF spectra of neutrons at the 5 meV detector triplet in McStas (a) and Geant4 (b) simulations with the pyrolytic graphite sample. The lines are only joining the points. Table 3 Peak incident rate of the 5 meV detector tubes in McStas and Geant4 simulations, with a pyrolytic graphite sample.

Detector tube
McStas (Hz) Geant4 (Hz) Top 1.67 Â 10 9 1.69 Â 10 9 Central 1.74 Â 10 9 1.75 Â 10 9 Bottom 1.56 Â 10 9 1.56 Â 10 9  which the 5.0 meV neutrons are Bragg-scattered not toward the Q-channel, causing the slight dip of the peak at 5 meV. The spread of neutrons after the sample is less significant than it was with pyrolytic graphite; fewer neutrons are lost on the way toward the 5.0 meV analysers. The time-averaged and peak incident neutron rates for the detector tubes in Tables 5 and 6 show that the neutrons are distributed more evenly between the tubes as a result of the flattened top of the energy spectrum. The maximum of the time-averaged incident rate for a single tube is found to be 41 MHz, with a peak incident rate of 1 GHz.
Comparing the time-averaged and peak rates with those acquired for pyrolytic graphite, both are lower by a factor of 1.7 but still of the order of 10 MHz for time-averaged and GHz for peak rates. This means that, even with a non-highest-case sample, the rates can be well above the capabilities of the standard 3 He detector tubes.
As demonstrated in this subsection, the McStas and Geant4 simulation results are in excellent agreement regarding the detector rates. For this reason, further simulations are only performed using the Geant4 model of the sample and the scattering characterization system. In the subsequent subsections multiple parameters are scanned in order to determine their effect on the incident detector rates, and to prove that the results presented above can be regarded as the highestcase incident rates. The changes in the incident detector rates due to modifying the studied parameters are expected to have the same trend for all single crystals, so all simulations are performed using the Y 2 O 3 sample.

Sample mosaicity
The mosaicity of the sample has multiple effects on the neutron beam Bragg-scattered on a selected scattered plane toward the Q-channel. A sample with higher mosaicity scatters neutrons of a wider energy range, as the higher spread of crystal plane orientations enables them to fulfil the Bragg criteria. This is also true for a wider incident angle range, meaning that neutrons of a divergent beam with higher incident angle have the possibility to be Braggscattered on the selected scattering plane. This higher spread of crystal plane orientations, on the other hand, lowers the probability of neutrons with energy and incident angle close to the ideal values being scattered. The cumulative effect is depicted in Fig. 13, showing the energy spectra of the scattered beam at the 2.7 meV analysers and the 5.0 meV detector triplet for a Y 2 O 3 sample with different mosaicities.  Table 4 Time-averaged neutron intensities at the sample and different parts of the scattering characterization system with a Y 2 O 3 single-crystal sample in McStas and Geant4 simulations.

Figure 12
Time-averaged neutron energy spectra at the sample and the scattering characterization system with a Y  The energy spectra of the beam at the 2.7 meV analyser arc show that the samples with higher mosaicities scatter neutrons of a wider energy range toward the analysers, as expected. It also shows that the intensities at this point of the instrument are getting higher for mosaicities up until 80 arcmin and therefore intensities for 80 arcmin are higher than those for 60 arcmin. The energy spectrum of neutrons hitting the detector triplet at 5.0 meV shows that these additional neutrons do not reach the detectors, as the highest intensities are found at 60 arcmin sample mosaicity. The resulting timeaveraged and peak incident rates of the central detector tube are presented in Table 7.
The results are in compliance with the expectation that the highest rates occur when the mosaicity of the sample matches that of the analysers, but also show that within the AE20 arcmin range it is a less than 10% effect.

Analyser mosaicity
The mosaicity of the analysers is a fixed value of 60 arcmin for BIFROST, but it is worth briefly investigating how it would affect the rate of the detector tubes. Fig. 14 depicts the neutron energy spectra of the 5.0 meV detector triplet for different analyser mosaicities with a Y 2 O 3 sample with a mosaicity of 60 arcmin. The spectra of the three tubes separately show that for mosaicities below 40 arcmin the two tubes on the sides are under-illuminated compared with the tube in the centre. In order to apply the prismatic analyser concept, the analyser mosaicity has to be large enough to sufficiently cover all used detectors. On increasing the mosaicity above 40 arcmin, the intensity in all three tubes slightly decreases. The resulting time-averaged and peak rates of the central detector tube are presented in Table 8.
The results show that the incident rate in a single detector tube could be 13-14% higher in the central tube with 20 arcmin mosaicity compared with the result for 60 arcmin, but the mosaicity has to be higher to apply the prismatic analyser concept, and in the range of 40-80 arcmin the change is less than 10%.

Sample size
Sample size is the limiting factor in many scientific cases, as it is not easy to grow large samples of some types. The beam delivery system of BIFROST is optimized for sample cross sections up to 15 Â 15 mm, but the realistic sample  Table 7 Time-averaged and peak incident neutron rates of the central 5 meV detector tube for different sample mosaicities with a Y 2 O 3 sample, where the mosaicity of the analysers is 60 arcmin.

Figure 14
Energy spectra of the neutrons hitting the detector triplet for the three tubes together (left) and separately (right) for 5.0 meV for different analyser mosaicities and a Y 2 O 3 sample with 60 arcmin mosaicity. The lines are only joining the points. Table 8 Time-averaged and peak incident neutron rates of the central 5 meV detector tubes for different analyser mosaicities and a Y 2 O 3 sample with 60 arcmin mosaicity.

Figure 13
Energy spectra of the neutron beam on the set of analysers for 2.7 meV neutrons (left) and of the neutrons hitting the detector triplet for 5.0 meV (right) for different sample mosaicities with a Y 2 O 3 sample. The mosaicity of the analysers is 60 arcmin. The lines are only joining the points.
sizes for the intended applications are much smaller than that, with an expected minimum sample size going down to 1 mm 3 . The height, width and thickness of the sample can have different effects on the incident detector rates. However, in this parameter scan their cumulative effects are investigated using cylindrical samples with equal diameter and height. The energy spectra of the scattered beam at the 2.7 meV analysers and the 5.0 meV detector triplet for Y 2 O 3 samples of different sizes are depicted in Fig. 15. The resulting time-averaged and peak incident rates of the central detector tube are presented in Table 9. As expected, the larger the sample, the higher the intensities. By reducing the sample size parameter (height and diameter) from 15 mm to 5 and 1 mm, the time-averaged incident rate of the centre tube drops by a factor of 10.5 and 900, respectively. Owing to the better resolutions in the case of smaller samples, the drop in the peak incident rate is lower, a factor of 10.3 for 5 mm and a factor of 650 for 1 mm.
Another effect of the better resolution is visible in the energy spectra of the detector triplets, where the three-peak structure is more apparent for smaller samples.

Pulse-shaping chopper opening time
The energy resolution of the instrument can be increased at the cost of neutron intensity by modifying the opening time of the pulse-shaping chopper. The fluxes on the sample and the detectors are both expected to drop significantly in the high-resolution setting, when the opening time is merely 0.1 ms, compared with the highflux mode achieved by an opening time of 5 ms.
The energy spectra of the neutron beam at the sample and the 5.0 meV detector triplet for a Y 2 O 3 sample for different pulse-shaping chopper opening times are depicted in Fig. 16, with the time-averaged intensities presented in Table 10. The timeaveraged and peak incident rates of the central 5.0 meV detector tube are presented in Table 11.
As expected, the time-averaged rates on the sample and on the detectors decrease with shorter pulse-shaping chopper opening times. The difference in time-averaged incident rates between the high-flux mode (5 ms) and the highresolution mode (0.1 ms) is approximately a factor of 20 for both the detector triplet and the central tube.
Regarding the peak rates in the central tube, however, this drop is less Energy spectra of the neutron beam on the set of analysers for 2.7 meV neutrons (left) and of the neutrons hitting the detector triplet for 5.0 meV (right) for a Y 2 O 3 sample of different sizes. The diameter and height of the cylindrical samples are equal, with the magnitude indicated in the legend. The mosaicity of both the sample and the analyser is 60 arcmin. The lines are only joining the points. Table 9 Time-averaged and peak incident neutron rates of the central 5 meV detector tube for Y 2 O 3 samples of different sizes, where the mosaicity of both the sample and the analyser is 60 arcmin.

Figure 16
Energy spectra of the neutron beam on the sample (left) and of the neutrons hitting the detector triplet for 5.0 meV (right) for a Y 2 O 3 sample for different pulse-shaping chopper opening times. The lines are only joining the points.

Table 10
Time-averaged neutron intensities at the sample and the 5.0 meV detector tubes with a Y 2 O 3 sample for different pulse-shaping chopper (PSC) opening times.
PSC opening time (ms) Sample (Hz) 5.0 meV detectors (Hz) 5.0 1.44 Â 10 10 1.21 Â 10 8 3.0 1.22 Â 10 10 1.14 Â 10 8 1.0 5.05 Â 10 9 4.38 Â 10 7 0.5 2.67 Â 10 9 2.27 Â 10 7 0.2 1.20 Â 10 9 1.01 Â 10 7 0.1 7.03 Â 10 8 5.93 Â 10 6 apparent. The highest rate for a 3 ms opening time is the same (within statistical uncertainty) as the rate for 5 ms, and the difference compared with 0.1 ms opening time is only a factor of 6.8. The reason for this difference is the better ToF resolution with shorter pulse-shaping chopper opening times. The higher time-averaged rates are distributed over a longer period of time on the detectors, as demonstrated in Fig. 17, showing the effect of the pulse-shaping chopper opening time on the ToF spectrum of neutrons hitting the detectors. The longer the opening time, the broader the ToF peaks. This directly affects the energy resolution and increases the dead time in the case of saturation. If one tube of a triplet is saturated, then none of the three can read out data, as they are connected in series. This means that for the detector triplet in the presented case for 5.0 ms opening time no data are recorded for more than 6 ms.

Elastic peak rates in representative operational conditions
The parameters chosen in Section 3.1 correspond to possible highest-case scenarios, and the rates acquired are far above the capabilities of 3 He tubes. However, the combination of a strongly scattering large sample and the highest flux mode is rather artificial, so it is worth evaluating a more representative operational scenario.
BIFROST is designed for small samples, as sample size is the limiting factor in many science cases. Hence, centimetre-sized crystals are not to be expected very often, only large samples with small magnetic moments and therefore small magnetic Bragg peak intensity. There is another parameter directly affecting the intensities but not discussed yet, the accelerator power of the ESS source. As mentioned in Section 2.2, the source power of 5 MW is used for the simulations. That is the eventual operational power of ESS, but it will initially operate at 2 MW. The intensities are expected to scale linearly with the source power.
For these reasons, the following parameters are selected to define the rates in a more representative operational case: 2 MW source power, 1 ms pulse-shaping chopper opening time, a Y 2 O 3 single-crystal sample with a height and diameter of 3 mm, and mosaicity of 60 arcmin. The time-averaged energy spectra of the neutron beam at the sample and different positions of the scattering characterization system acquired with these parameters are demonstrated in Fig. 18, and the integral values are presented in Table 12.
The combined effect of the lower source power, shorter pulse-shaping chopper opening time and smaller sample (cross section) decreased the time-averaged neutron intensity on the sample significantly, by a factor of 136 compared with the highest-case ToF spectrum of neutrons at the 5 meV detector triplet for the three tubes together (left) and separately (right) with a Y 2 O 3 sample for different pulse-shaping chopper opening times. The separated spectrum is not displayed in all cases in the right-hand figure to avoid the figure being overcrowded. The lines are only joining the points.

Figure 18
Time-averaged neutron energy spectra at the sample and the scattering characterization system with a Y 2 O 3 sample in the Geant4 simulation using a 2 MW source power, 3 mm height and diameter sample size, and 60 arcmin sample mosaicity. Incident beam on sample (in blue), beam transmitted through the sample (in orange), beam on the set of analysers for 2.7 meV neutrons (in green), beam on the set of analysers for 5.0 meV neutrons (in red), neutrons hitting the detector triplet for 5.0 meV (in purple), beam transmitted through all sets of analysers (in brown). The lines are only joining the points.

Table 11
Time-averaged and peak incident neutron rates of the central 5 meV detector tube with a Y 2 O 3 sample for different pulse-shaping chopper opening times.

PSC opening time (ms)
Time-averaged rate (Hz) Peak rate (Hz) 5.0 4.13 Â 10 7 1.04 Â 10 9 3.0 3.91 Â 10 7 1.04 Â 10 9 1.0 1.50 Â 10 7 8.76 Â 10 8 0.5 7.79 Â 10 6 5.31 Â 10 8 0.2 3.46 Â 10 6 2.52 Â 10 8 0.1 2.03 Â 10 6 1.50 Â 10 8 scenario with a Y 2 O 3 sample. Owing to the reduced sample thickness, the transmission through the sample is increased to 94 from 77%, as a result of lower absorption and weaker Bragg peaks. The lower incident intensity on the sample and weaker Bragg peak lead to a drop by a factor of 322-326 in the time-averaged neutron intensity on both the 2.7 and 5.0 meV analysers. The drop in the time-averaged neutron intensity on the 5.0 meV detector triplet is sightly lower, a factor of 295, because of the smaller divergence and better energy resolution of the neutron beam compared with the highest-case scenario.
The resulting time-averaged and peak incident neutron rates of the 5 meV detector tubes are presented in Table 13. The highest time-averaged incident neutron rate on a single tube is 0.15 MHz, which means a drop by a factor of 275, but the peak incident rate is 9.9 MHz, which is lower only by a factor of 105 compared with the highest-case results.
The numbers and reduction factors are in accordance with previous simulations in Sections 3.4 and 3.5, where the effect of sample size and pulse-shaping chopper opening time were investigated separately.
There are factors not taken into account in the current study that may further reduce the intensities slightly on the detectors, like the non-ideal transmission of the filtering system and the effect of the divergence jaws applicable for reducing the angular spread of neutrons.

Simulation with full scattering characterization system
The previous sections aimed to define the highest incident rates a detector tube can experience using different instrument and sample parameters in the case of a coherent elastic (Bragg) peak. This section presents the incident detector rates in the case of incoherent elastic peaks with a standard cali-bration sample, and demonstrates the use of the full simulation model of the scattering characterization system with all nine Q-channels.
The sample selected for this simulation is vanadium, which is assumed to be an incoherent elastic scatterer that scatters isotropically, and therefore it is used to calibrate the incident neutron intensity and the detector efficiencies in neutron spectrometers (Mayers, 1984). As mentioned earlier, in previous simulations the single Q-channel present in the model was rotated according to the Bragg angle of the sample for 5 meV neutrons ( = 37.067 for pyrolytic graphite and = 41.169 for Y 2 O 3 ). For vanadium the rotation of the nine Q-channels is arbitrarily selected so as to have 2 = 90 scattering angle for the central Q-channel, as depicted in Fig. 19.
The instrument and sample parameters are the same as for the highest-case scenario: 5 MW source power, 5 ms pulseshaping chopper opening time, 15 mm sample height and diameter. The time-averaged energy spectra of the neutron beam at the sample and different positions of the central Q-channel are demonstrated in Fig. 20, with the integral values presented in Table 14.
The transmission through the sample is only 23% -much lower than it is for pyrolytic graphite (80%) or Y 2 O 3 (77%)with no peaks missing from the spectrum, as expected from a sample scattering mainly incoherently. The wide spectrum of the 2.7 meV analysers also shows that neutrons are not coming from an elastic peak, but despite the wider energy range, a logarithmic scale is needed as the integrated intensity is more than two orders of magnitude lower than experienced with previous samples. The spectrum of the beam transmitted through all five sets of analysers clearly shows neutrons missing because they are selected by the analysers. These neutrons appear in the spectra of the detectors, which show that the peaks are narrower for lower energies as a result of the better energy resolution of the analysers for lower energies. The incident intensity on the 2.7 meV detectors is much Top-view schematic figure of the scattering characterization system model with all nine Q-channels. The red lines and corresponding angles indicate the scattering angle for the centre of each Q-channel.

Table 12
Time-averaged neutron intensities at the sample and different parts of the scattering characterization system with a Y 2 O 3 sample, using 2 MW source power, 1 ms PSC opening time, 3 mm height and diameter sample size, and 60 arcmin sample mosaicity.
Detector tube Time-averaged rate (Hz) Peak rate (Hz) Top 1.3 Â 10 5 9.3 Â 10 6 Central 1.5 Â 10 5 9.9 Â 10 6 Bottom 1.3 Â 10 5 9.4 Â 10 6 lower than in other detectors because of the energy range selected by the bandwidth chopper. The resulting timeaveraged incident rates are higher for higher energies, with the maximum of 26 AE 2 kHz for the 5.0 meV detector triplet. For a single detector tube the highest rates are found for the central 5.0 meV detector with a time-averaged intensity of 9 AE 1 kHz and peak intensity of 0.3 AE 0.1 MHz. The time-averaged incident neutron rates of all detector triplets in all Q-channels are presented in Table 15. The trends in the results demonstrate the combination of three effects. In each Q-channel the detector triplets for higher energies experience higher incident rates due to the wider energy ranges selected by the analysers, as shown for the central Q-channel earlier in this section. The second effect has roots in the 'triple stagger' geometry and the asymmetry of the Q-channels described in Section 2.1. The sample-analyser distances in Q-channels 1, 4 and 7 are shorter and in channels 3, 6 and 9 they are longer than the distances in channels 2, 5 and 8. The shorter distances increase the rates visibly because neutrons are not collimated by Bragg scattering on the sample and therefore their spread at longer distances becomes important. This effect on the rates is somewhat blurred by the third effect, caused by the anisotropy of the scattering cross section in vanadium. For the observed energies the scattering cross section of vanadium is slightly higher for higher scattering angles (Mayers, 1989), and more importantly, the neutron path length through the solid cylindrical sample is generally higher for neutrons scattered at lower angles. Therefore the absorption is higher for these neutrons. These two effects result in generally higher rates for Q-channels positioned for higher scattering angles, but owing to the asymmetry of the adjacent Q-channels, it is most apparent when comparing Q-channels of the same symmetry, like 2, 5 and 8.

Conclusions
BIFROST is an indirect ToF spectrometer at ESS and one of the first eight instruments to be constructed. One of the most challenging aspects of its operation is the rate capability and in particular the peak instantaneous rate capability, i.e. the number of neutrons hitting the detector per channel at the peak of the neutron pulse. There is no intent to measure the intensity of elastic peaks as they are considered background for this instrument. However, it is vital that the detectors are not degraded by such intensity and remain capable of measuring weak inelastic signals, as soon as possible after saturation. This implies that the detector aspects of recovery time and high rate tolerance have to be carefully evaluated by measurements to prove that scientific performance will be intact.
A detailed methodology for acquiring the results is presented. The full simulation of the instrument from source to detector position has been carried out using multiple simulation software packages. Flexible models of the sample and the scattering characterization system of BIFROST were implemented in both McStas and Geant4 and a comparison of their strengths and weaknesses is presented. The capability of both simulation tools is enhanced by the NCrystal library and associated tools. The first application of the special NCrystal pyrolytic graphite is presented, demonstrating its capabilities for modelling analysers for neutron scattering applications.  Table 14 Time-averaged neutron intensities at the sample and different parts of the scattering characterization system with a vanadium sample in a Geant4 simulation.

Figure 20
Time-averaged neutron energy spectra at the sample and in the central Q-channel with a vanadium sample in the Geant4 simulation. Incident beam on sample (in blue), beam transmitted through the sample (in orange), beam on the set of analysers for 2.7 meV neutrons (in green), neutrons hitting the detector triplets for energies 2.7-5.0 meV (in red-grey), beam transmitted through all sets of analysers (in mustard). The lines are only joining the points.

Table 15
Time-averaged neutron intensities of the five detector triplets in all nine Q-channels with a vanadium sample. McStas is capable of simulating instruments as long as 160 m, and even handling beam splitting to some extent, to treat simulations with multiple sets of analysers. However, the latter comes with great complexity and some limitations, as it is not within the natural usage of this simulation software. Geant4, on the other hand, is not suited to simulating the beam transport system of an instrument, but with the use of NCrystal, it is an entirely appropriate tool for a scattering characterization system with any level of geometrical complexity, and even offers the possibility to include parts like a filtering system and cross-talk shielding and to take into account back-scattering. The results of the McStas and Geant4 models of the scattering characterization system were compared using various single-crystal samples. The results show perfect agreement, the only exception being the transmission through the sample where a difference of less than 10% is found in one case, due to the more detailed modelling of absorption in Geant4.
With this knowledge at hand a choice was made to combine the McStas model of the beam transport system and the Geant4 model of the sample and the scattering characterization system using the MCPL tools. Using this model, the incident detector rates anticipated at the BIFROST instrument for different configurations are presented. The impact of sample type, sample and analyser mosaicity, sample size, and pulse-shaping chopper opening time was studied on the incident detector rates. For instrument configurations and sample parameters representing highest-case conditions, the peak rate can reach the value of 1-1.7 GHz for a single detector tube with time-averaged rates of 40-70 MHz. These tubes are expected to reach saturation well below that, at 50-100 kHz. These tubes will also be saturated for a minimum of 5 ms, but the saturation deadtime for detecting signals is more like 6 ms because the counting detector tubes are coupled in triplets.
To overcome challenges caused by these rates, an operational evaluation of a measurement strategy will be the key to the successful operation of this instrument. More 'everyday' realistic samples give a lower rate challenge. However, these samples will still saturate detectors.
A simulation with the full analyser system is presented using a common calibration sample. This model can now be used to predict experimental conditions from specific proposed samples, i.e. sample size and composition for experiment planning purposes for users.
The results here show the potential power of source-todetector simulation for neutron scattering. These simulations are feasible as a result of tools recently developed. It is now possible to realistically simulate very complex systems.