NCrystal : a library for thermal neutron transport

An open source software package for modelling thermal neutron transport is presented. The code facilitates Monte Carlo-based transport simulations and focuses in the initial release on interactions in both mosaic single crystals as well as polycrystalline materials and powders. Both coherent elastic (Bragg diffraction) and incoherent or inelastic (phonon) scattering are modelled, using basic parameters of the crystal unit cell as input. Included is a data library of validated crystal definitions, standalone tools and interfaces for C++, C and Python programming languages. Interfaces for two popular simulation packages, Geant4 and McStas, are provided, enabling highly realistic simulations of typical components at neutron scattering instruments, including beam filters, monochromators, analysers, samples and detectors. All interfaces are presented in detail, along with the end-user configuration procedure which is deliberately kept user-friendly and consistent across all interfaces. An overview of the relevant neutron scattering theory is provided, and the physics modelling capabilities of the software are discussed. Particular attention is given here to the ability to load crystal structures and form factors from various sources of input, and the results are benchmarked and validated against experimental data and existing crystallographic software. Good agreements are observed.


Introduction
The modelling of neutron transport through matter dates back to the efforts aimed at tackling neutron diffusion problems in the middle of the twentieth century, closely tied to the introduction of general purpose computers and the inception of the method of Monte Carlo simulations [1]. Since then, the needs for transport simulations involving neutrons have expanded, with a wide range of applications in radiation protection, reactor physics, radiotherapy, and scat-tering instruments at spallation or reactor sources. To facilitate such simulations a plethora of Monte Carlo simulation applications exists today, which for the purposes here can be divided into two categories: general purpose applications [2,3,4,5,6,7,8,9,10,11], capable of modelling a variety of particle types in flexible geometrical layouts, and specialised applications [12,13,14,15,16,17] dealing with the design of neutron scattering instruments. Applications in the latter group is focused on aspects of individual scatterings of thermal neutrons, but generally lack capabilities for dealing with arbitrarily complex geometries and the inclusion of physics of particles other than thermalised neutrons.
On the other hand, the general purpose Monte Carlo applications are typically oriented towards use-cases in fields such as reactor physics, high energy physics, and radiation protection, and as a rule neglect all material structure and effects of inter-atomic bindings -resorting instead to the approximation of treating all materials as a free gas of unbound atoms with no mutual interactions. This approximation is suitable at higher incident energies, but for neutrons it breaks down at the thermal scale where their wavelengths and kinetic energies are comparable to the typical distances and excitational energies resulting from inter-atomic bindings. The only presently available option for improved modelling of thermal neutron interactions involves the utilisation of specially prepared data files of differential cross section data. These "scattering kernels" are only readily available in neutron evaluation libraries, e.g. [18], for the dozen or so materials which have traditionally played an important role for shielding or moderation purposes at nuclear facilities. For other materials, scattering kernels must be carefully crafted using non-trivial procedures (like the application of delicate and computationally expensive molecular dynamics models [19]) and in practice this is rarely done. Furthermore, some aspects of thermal neutron scatterings are in practice ignored when modelling is exclusively based on scattering kernels: since nuclear reactor physics has historically mostly been focused on inelastic and multiple-scattering phenomena, it is usually not possible to include a precise model of the a priory highly significant process of Bragg diffraction into the setup for crystalline materials. At thermal energies, Bragg diffraction is the dominant process for many relevant materials, and although it is an elastic scattering process which does not change the neutron energy, it sends neutrons to characteristic solid angles and can therefore critically influence the geometrical reach of neutrons through a given geometry.
It is thus not currently straight-forward to carry out simulations which at the same time incorporate consistent and realistic physics models for neutrons at both high and low energy scales in general materials, while simultaneously supporting flexible geometrical layouts and the treatment of particles other than neutrons. Such simulations would nonetheless be remarkably useful, in particular when considering aspects of neutron scattering instruments [20,21]. Here, precise modelling of individual neutron scatterings is crucial in all beam-line components, while detailed understanding of background levels, the effects of shielding, or the workings of neutron detectors, require incorporation of detailed geometrical layouts and additional physics such as modelling of gammas, fast neutrons, and energy depositions resulting from secondary particles released upon neutron capture. Other examples include the design of advanced moderators or reflectors for novel reactors or neutron spallation sources, and proton or neutron-based radiotherapy, in which adverse irradiation from neutrons are abundant and an increasing concern [22,23].
The NCrystal toolkit presented here is aimed at remedying the situation.
Rather than introducing an entirely new application, it is designed as a flexible backend, able to plug the gaps where existing Monte Carlo applications are lacking in capabilities for treating thermal neutrons in structured materials. It was originally developed to facilitate the design and optimisation of neutron detectors for the European Spallation Source, presently under construction in Lund, Sweden [24,25,26,27,20], but is intended to be as widely useful as possible -not only as a service to the community, but also to ultimately ensure a higher level of quality for the toolkit itself due to feedback and validation from a larger potential userbase. Although dedicated utilities for calculating various properties of thermal neutrons exists [28,29,30,31], the scope of NCrystal is different, delivering at the same time an ambitious set of physics models, a flex-  [32]. 1 Included in the release is additionally tools, examples, data files, and a configuration file for CMake [33] with which everything can be built and installed. A website [34] has been set up for the NCrystal project, on which users will be able to interact with developers, locate future updates, and access relevant documentation.
Concerning NCrystal's capabilities for physics simulations, the work has so far focused on facilitating realistic simulations of both inelastic and elastic physics of thermal neutrons as they scatter in certain common bulk crystalline materials. Namely those which can be described in terms of a statistical arrangements of microscopic perfect crystals, and for which the neutron's interaction probabilities are not so strong as to break down the Born approximation. Additionally, material configuration can be carried out based simply on a brief description of the crystal structures. All together, the framework thus provides realistic thermal neutron physics in a large range of powdered, polycrystalline or mosaic single crystals readily used at neutron scattering facilities in beam filters, monochromators, analysers, samples, detectors, and shielding. What is currently not supported, given the mentioned constraints, is surface effects, materials that are either non-crystalline or textured, and effects related to a break-down of the Born approximation including the strong reflections found in neutron guide systems, or the strong diffraction effects encountered in certain macroscopic crystal systems such as synthetically grown Silicon. Future developments might address some of those deficiencies, depending on community 1 Optional components adding support for file formats discussed in Sections 4. interest and resources.
The present paper will present the NCrystal framework itself, including interfaces, configuration, and data input options. Although a brief overview of physics capabilities will be provided, detailed discussions of the implementation of actual thermal neutron scattering models are beyond the scope of the present paper, and will be the subject of dedicated publications at a later date. After a review of the relevant theoretical concepts in Section 2, the design and implementation of the core NCrystal code is presented in Section 3. Section 4 presents the capabilities for initialising relevant parameters of modelled crystals from various data sources, introduces the library of data files included with NCrystal and discusses how the loaded results have been validated. Section 5 presents the uniform method for material configuration intended for most endusers, and an overview of the existing language bindings and application interfaces is provided in Section 6. Finally, a discussion of future directions and planned improvements is carried out in Section 7.

Theoretical background
Rather than intending to be an exhaustive treatment, this section will provide a brief overview of crystals and neutron scattering, introducing relevant concepts, models and terminology needed for the purposes of the present discussions of NCrystal. Naturally, relevant references to more extensive background literature are provided for readers seeking more information.

Crystals
Although support for other types of materials is likely to be added eventually, the present scope of NCrystal is restricted to crystals. Consequently, a few basic concepts will be introduced in the following. Further details can be found in [35,36,37,38].
Informally, a crystal can be described as an arrangement of bound atoms which is built up by the periodic repetition of a basic element in all three spatial directions. This repeated element is known as the unit cell and is typically chosen such that it is the smallest such cell that reflects the symmetry of the structure. In an idealised setting, the crystal structure would be infinite, but in practice it is sufficient that the unit cell is much smaller than the entire structure. The choice of unit cell for a given crystal is not always unique, but it is always possible to select one which is a parallelepiped spanned by three linearly independent basis vectors a, b and c. 2 By convention they are defined by their respective lengths, a, b, and c, and the angles between them: α ≡ ∠( b, c), β ≡ ∠( a, c), and γ ≡ ∠( a, b). The crystal structure is completed by the additional specification of the contents of the unit cell, often provided in coordinates (x, y, z) relative to the crystal axes. The position of the ith constituent in the unit cell is thus given by: Crystals exhibit a discrete symmetry under spatial translation with the vectors: Apart from a trivial global offset, coordinates of all constituents in the entire crystal are given by R mno + p i for all combinations of mno and i. The translational symmetry means that any function representing features of the crystal system (such as particle densities) will obey: Additional symmetries under rotations, reflections or inversions are possible, depending on the detailed shape and contents of the unit cell. The symmetry properties of a crystal are described by the concept of space groups, and it is possible to describe any spatial crystal symmetry by one of only 230 such 2 Vectors are here and throughout the text donated with arrows ( a), while the absence of arrows indicate the corresponding scalar magnitudes (a ≡ | a|). Additionally, unit vectors are donated with hats (â ≡ a/a) and quantum mechanical operators are shown in bold (V, R).
groups. These groups can be divided into 7 general classes of crystal systems.
Loosely ordered from lowest to highest symmetry they are: triclinic, monoclinic, orthorhombic, tetragonal, trigonal, hexagonal, and cubic. The space groups are defined and described exhaustively in [38].
The space group of a given crystal structure naturally depends on the unit cell positions of constituents, p 1 , . . . , p N . Conversely, it is possible to construct the full list of these positions by applying the symmetry operators of a given space group to a smaller number of so-called Wyckoff positions [39]. It is indeed very common to find crystal structures defined solely in terms of their space group, unit cell shape and dimensions, and a list of elements and associated Wyckoff positions.
The set of points at positions R mno for all integers m, n, and o, constitutes the crystal lattice. In each such lattice, there exist an infinite number of families of evenly spaced parallel planes passing through the lattice points. Each such family of lattice planes is described by a common plane normaln and an interplanar distance d (also known as the d-spacing), with planes passing through the points jdn for all j ∈ Z. Particle diffraction in a crystal occurs in the form of reflections from such families of lattice planes. The task of finding and classifying these families is often done in the so-called reciprocal lattice in momentum-space, which is the Fourier transform of the direct lattice. The basis vectors of the reciprocal lattice are given by: Where V uc = a · ( b × c) is the unit cell volume. The points in the reciprocal lattice are thus given by: As will be motivated in Section 2.3, it can be shown that each point in the reciprocal lattice corresponds to a family of equidistant lattice planes in the direct lattice. The family of planes corresponding to τ hkl will be indexed by the hkl value (known as its Miller index), and has interplanar spacing d hkl = 2π/τ hkl and normaln =τ hkl = τ hkl /τ hkl .
In most real systems, regions in which the crystal lattice is continuous and the defining translational invariance of Eq. 3 is unbroken, exists only at the microscopic scale. 3 Macroscopic crystals are then built up from these independently oriented grains of perfect crystals (also known as "crystallites"), and the actual distribution of orientations will not only be related to various macroscopic material properties, but will influence interaction probabilities when the system is probed with impinging particles. As long as the grain sizes are small compared to the regions of material being probed, these effects can be accounted for in a statistical manner -for instance by integrating microscopic interaction cross sections over the distribution of crystallite orientations in order to arrive at effective macroscopic cross sections. Such integrations are particularly trivial to perform in the extreme case where the orientation of individual crystallites are completely independent and uniformly distributed over all solid angles. This model, known as the powder approximation, is not only suitable for modelling actual powdered crystalline samples, but can also be used to approximate interactions in polycrystalline materials like metals or ceramics -especially when the level of correlation in crystallite orientation ("texture") is low or when the setup involves sufficiently large geometries or spread in incoming particles that effects due to local correlations are washed out. An example of this would be the case where polycrystalline metal support structures are placed in various places throughout a neutron instrument. On the other hand, interactions in a highly textured polycrystalline sample placed in a tightly focused beam of probe particles might not be very well described by the powder approximation.
In another extreme, all crystallites are almost co-aligned in so-called mosaic crystals, with the exact distribution around the common axis of alignment given by the mosaicity distribution. In the simple isotropic case, this distribution is a Gaussian function of the angular displacement, with the corresponding width 3 Exceptions to this exist, like in some synthetic silicon crystals. However, the scattering theory discussed in Section 2.2 is in any case not strictly applicable to such systems.
(referred to as the mosaicity of the crystal) having typical values anywhere from a few arc seconds to a few degrees. In NCrystal, such Gaussian mosaic crystals are referred to as "single crystals", due to the large degree of crystallite alignment throughout the material. Some materials are better described by other mosaicity distributions, including anisotropic ones where the distribution is not symmetric with respect to all axes. A particular anisotropic distribution supported by NCrystal is one in which crystallites are aligned around one particular axis, but appearing with random rotations around the same axis. Such distributions occur in stratified or layered crystals, like pyrolytic graphite, in which certain planes of atoms tend to have strongly aligned normals in all crystallites, but no strong alignment for the rotation of the same planes around their normals.

Scattering theory
The theory of thermal neutron scatterings in condensed matter is thoroughly treated in the literature, see for instance [40,41,42]. The discussion in the present and following sections is to a large degree inspired by [40,41].
Scattering in materials by thermal neutrons, taken here to loosely mean energies below or at the eV scale or wavelengths above or at the nm scale, can be described via point-like neutron-nuclei interactions or magnetic dipole interactions between the neutron and any unpaired electrons in the target material.
For unpolarised neutrons and samples, the latter can be largely ignored [43] and the current discussion will not concern itself with such magnetic interactions, nor are they currently modelled in NCrystal -although there is nothing fundamental preventing their inclusion in the future.
Thermal neutrons carry such low energy that in practice no emission of secondary particles occurs during pure scattering interactions. As a consequence, scattering of unpolarised neutrons can be completely described in terms of the differential cross section for scattering the incident neutron with wavevector k i into the wavevector k f . Based on the Born approximation or Fermi's golden rule, this differential cross section is related to quantum mechanical transition amplitudes by the so-called master equation for thermal neutron scattering: Here m is the neutron mass, indices i and f denotes initial and final states respectively, η j represents a complete set of eigenstates of the target system, p(η i ) represents the initial occupation levels of target states, and V is the interaction total energy change of the neutron and target system respectively. At this point it might be useful to recall that in addition to wavenumber, k = | k|, the energy state of a non-relativistic neutron might equally be characterised in terms of energy E, momentum p, velocity v, or wavelength λ, related by p = k, 2E = mv 2 = p 2 /m, and λ = 2π/k.
Before continuing with the evaluation of Eq. 6, it is important to understand its validity and applicability. Based on the Born approximation, the perturbation experienced by the initial neutron wave function during the scattering is assumed to be small, so it is possible to neglect secondary scattering of waves created as a result of the primary scattering. The sum in Eq. 6 implies that the resulting cross sections will grow with the size of the target system, and therefore the validity of the equation will eventually break down as the size of the target system grows and the probability for secondary scattering becomes nonnegligible. In reality this issue is circumvented in the context of Monte Carlo transport simulations in which neutron transport is modelled in a series of steps between independent scatterings, with each simulated step length depending on the actual mean free path length predicted by the cross section in Eq. 6. 4 High cross sections will automatically lead to small step lengths, and therefore only a small part of the sample will be traversed for each application of Eq. 6.
Simulations employing such stepping will therefore to some degree account for multiple scattering phenomena like extinction effects in crystal diffraction. This effective subdivision of a target system into smaller decoherent systems is, however, only possible if the linear scale over which any symmetries in the material structure exists is small compared to the mean free path lengths involved, so wavefunctions originating from scattering in separate subsystems will be mutually incoherent. As described in Section 2.1 most real macroscopic crystals consists of independently aligned microscopic crystal grains, and in practice this ensures the necessary decoherence between different parts of the macroscopic system. The exception is the case of perfect macroscopic crystals (or mosaic crystals with vanishing mosaic spread) like some synthetic silicon crystals. For an excellent and detailed discussion about these issues, refer to [41,Sec. 11.5].
In order to proceed with the evaluation of Eq. 6, one must provide suitable representations of both interaction potential V and distributions of target states, p(η i ). When dealing with thermal neutrons and their relatively long wavelength, the short-range neutron-nuclei interactions can effectively be described as being point-like. An effective potential which describes them as such was introduced by Fermi [44]. For a collection of N nuclei located at positions R j it looks like: Here b j is the scattering length of the jth nucleus, an effective parameter capturing the details of the neutron-nucleus interaction which must be determined through measurements. The term scattering length stems from the fact that the corresponding cross section for scattering from a single fixed nucleus is 4πb 2 , equal to the classical cross section of a hard sphere of radius 2|b|. In the presence of nuclear resonances, the scattering length becomes strongly dependent on the neutron energy. However, at the energy scale of thermal neutrons, resonances only exist for a few rare isotopes -and even then mostly at the high end of the thermal scale [45]. For the present purposes, the scattering length thus attain constant values, specific to each type of nuclei. It is additionally possible to use non-zero imaginary components of scattering lengths to represent absorption physics, however in many contexts -including the present -the two types of interactions are dealt with separately and the scattering lengths will therefore be treated as real and constant numbers. Such separation between scattering and absorption processes is valid at the O(10 −4 ) level for thermal neutrons [46], which is considered acceptable. In fact, as the current scope of NCrystal does not cover physics of nuclear resonances, absorption processes will be dealt with using the simple -but in most cases accurate [45] -model of absorption cross sections being inversely proportional to the neutron velocity. This 1/v scaling can be intuitively interpreted as having the probability of absorption by a given nucleus proportional to the time spent by the neutron near the nucleus, but also follows from more careful reasoning [46,47,48].
It is possible to bring Eq. 6 into a form in which the sum over final states of the target system is removed and the potential is replaced with its Fourier transform. This is particularly convenient for potentials involving δ-functions like Eq. 7, as these transform into constants. The resulting form is: where Q is the momentum transfer, Q ≡ k f − k i , and S( Q, ω) is the scattering function defined by: The . . . notation is here used to donate operator expectation values in the and will here be abbreviated as: Leading to the shorter expression for the scattering function: This definition of the scattering function contains a sum over the target system under consideration (which as discussed must be of a linear size compatible with the applied Born approximation). Most such systems of interests can be divided into a number of statistically equivalent subsystems. That this is possible for crystals is obvious, since each unit cell forms such a subsystem, but subdivisions are generally possible in almost all systems, be they liquid, polymeric or gaseous in nature. The scattering function for a subsystem can be expressed as: (12) where now N refers to the subsystem size and might for instance represent the number of nuclei in a crystal unit cell, and b j b j represents an average performed over an ensemble of equivalent subsystems. Now, neutron-nuclei scattering has the peculiarity that the scattering lengths are isotope and spin-state dependent, whereas the positions of scatterers, the nuclei, are determined by chemical properties, depending on the nuclear charge but otherwise independent of isotope or spin state. 5 This independence means that the ensemble average of products of scattering lengths found at positions j and j obeys: 5 Of course, this statement becomes invalid for material temperatures near absolute zero where the energy difference between different nuclear spin states becomes comparable to thermal energies.
With this, Eq. 12 can be written as the sum of two distinct contributions: The terms in the coherent scattering function S coh ( Q, ω) directly depend on material structure as they involve pair correlations between nuclei at all po- the sum of the incoherent and coherent scattering cross sections is to a good approximation equal to the unbound scattering cross section, σ free , for which the total scattering cross section will converge at shorter wavelengths [46]. 6 The main difficulty in the evaluation of the scattering functions in Eqs. 15 and 16 is the integral over states in j , j implicit in Eq. 10, whose evaluation in principle relates to the potentially complicated time-dependent distribution 6 At very short wavelengths this conversion is ultimately disrupted by nuclear resonance physics or P-wave contributions.
of nuclei in the target system. In a crystal, the nuclear positions can be decomposed in terms of displacements u from equilibrium positions d: And thus: Under the assumption that displacements are small compared to inter-atomic distances, motion of nuclei can be described with potentials which are quadratic functions of the displacements. In this so-called harmonic approximation, nuclei are essentially described as being pulled towards their equilibrium positions by linear spring-like forces, with resulting harmonic vibrations. Working in this approximation, it is possible to show that: where the Debye-Waller function gives a measure of the time-independent average squared displacement of nuclei j along Q: The time dependence in Eq. 19 enters exclusively in the last exponential factor, which thus involves nuclear motion and can be expanded in a Taylor series as: The expectation value ( Q · u j (0))( Q · u j (t)) correlates linear displacements along Q of two nuclei at different times. It is possible to show that this is directly related to an interaction in which the neutron exchanges energy and momentum with a phonon state, and the expansion in Eq. 21 thus lends itself to physics interpretation in the phonon picture, with the nth term corresponding to n-phonon interactions. At lower displacements or Q values the expansion converges more rapidly, and accordingly multi-phonon physics will be less important at lower neutron energies or material temperatures.

Elastic scattering
The first term in the expansion of Eq. 21 gives rise to elastic ( k i = k f ) scattering when inserted into Eq. 15 or Eq. 16, since: With this, the partial differential cross section for incoherent elastic scattering can be immediately found: The only directional dependency enters here in the Debye-Waller factor, e −2Wj ( Q) , which approaches unity when the neutron wavelength is much larger than the atomic displacements, implying isotropic scattering. At higher neutron energies or temperatures this approximation breaks down, and the scattering will be increasingly focused in the forward direction at low values of the scatter angle θ, since for elastic scattering Q = 2k i sin(θ/2).
The coherent scattering terms involve correlations between pairs of nuclei and their evaluation will therefore be considerably more complex. In a crystal the sum j over all nuclei can be rewritten as a sum over unit cells first and unit cell contents second: l i . Here the index l assumes all R mno values from Eq. 2 and the index i runs over all nuclei in the unit cell. Thus, the equilibrium positions formerly denoted d j now becomes l + p i , with p i the positions defined in Eq. 1. Reordering sums appropriately and proceeding as for the incoherent case, the expression for coherent elastic scattering in crystals becomes: Where the form factor of the unit cell have been introduced: The parenthesised factor in Eq. 24 accounts for interference between different unit cells. Given that all l = R mno for suitable choice of integers m, n and o, the translational invariance of Eq. 3 implies: Where N uc is the number of unit cells in the crystal, which can be removed by adopting the convention of providing cross sections normalised per unit cell.
The remaining factor can be investigated further by expressing Q in terms of the basis vectors of the reciprocal lattice from Eq. 4, Q = q a τ a + q b τ b + q c τ c : Using: Eq. 27 becomes: Thus, coherent elastic scattering occurs only when the momentum transfer, Q, is exactly identical to one of the points in the reciprocal lattice, τ hkl , corresponding to interaction with the associated family of lattice planes. The fact that Q = τ hkl further supports the interpretation of reflection by lattice planes with normal along τ hkl , since elastic specular reflection by a mirror will always have the normal proportional to the transferred momentum. As ( τ a , τ b , τ c ) constitutes an orthogonal but not an orthonormal base: Where the last equality was found by inserting the definitions from Eq. 4 and evaluating the resulting cross product of two cross products with the BAC-CAB rule. With this result, the coherent elastic scattering function (normalised to the unit cell) can be written more compactly: The Debye-Waller factors will tend to suppress the form factors at high momenta, equivalent to large absolute values of h, k, and l. The summation range of these indices can therefore in practice be limited to a region around 0. This will be explored further in Section 4.1 where the impact on total coherent elastic cross sections due to a lower bound, d cut , on d hkl (corresponding to upper bounds on | τ hkl |) is investigated systematically for a large number of crystals.
For the hkl values selected for consideration it is then possible to pre-calculate the form factors, at the corresponding hkl indices. As will be discussed in sec- The scattering described by Eq. 31 is usually referred to as Bragg diffraction. Reflections from the family of lattice planes indexed by hkl requires and Q = 2k i sin(θ/2) where θ is the scattering angle. The maximal value of Q is found when θ = π and is 2k i . Thus, reflection is impossible unless (the Bragg condition): And the scattering angle θ will satisfy the Bragg equation: Often the above equation will be stated using the Bragg angle, defined as θ B ≡ θ/2. It also often contains on the left hand side an integer factor, n ≥ 1, denoting the so-called scattering order. However, this is merely a convenient manner in which to include multiple co-oriented lattice plane families in a single equation, utilising the fact that τ nh,nk,nl ∝ τ hkl and d nh,nk,nl = d hkl /n.
For a crystal powder, in which crystal grains appear with uniformly randomised orientations, the total coherent elastic cross section for scattering on a given family of lattice planes satisfying Eq. 32, can be found as an isotropic average over the orientation of τ hkl with respect to k i . Considering just the δ-functions, denoting the cosine of the angle between τ hkl and k i with µ, and dropping the hkl indices for brevity, this isotropic average gives: Where it was used that dΩδ( r − a) = 2a −1 δ(r 2 − a 2 ). Inserting the omitted factors from Eq. 31 in Eq. 34, it is evident that the complete coherent elastic cross section of a crystal powder can be written as the following sum over all lattice plane families satisfying Eq. 32: In case of a scattering event, the relative probability for it to happen on a particular hkl plane will depend on its contribution to the sum in Eq. 35, and the scattering angle θ will subsequently be determined by Eq. 33. The azimuthal scattering angle around the direction of k i is, however, not constrained and all such angles contribute equally to the cross section. Thus, k f will be distributed uniformly in a cone around k i with opening angle θ. Such cones are known as Debye-Scherrer cones, and are a prominent feature at any powder diffractometer.

Inelastic scattering
The n ≥ 1 terms in the expansion of Eq. 21 can be interpreted as inelastic scatterings in which the incoming neutron exchanges energy and momentum with n phonon states. The evaluation of these terms can in principle be very complicated depending on the material in question and the required precision.
For the purposes of the present publication, inelastic processes play only a minor role, and the details of the implementations of such processes in NCrystal is reserved for a future dedicated publication, with just a brief overview provided here.
For non-oriented materials like liquids and crystal powders, rotational symmetry implies that the scatter functions S( Q, ω) become independent of the direction of Q around k i , and thus become two-dimensional: S( Q, ω) = S(Q, ω).
As inelastic scattering terms do not include the factors of δ( ω) found in elastic terms, S(Q, ω) will generally be sufficiently smooth that it is possible to interpolate it from its values on a discrete grid in (Q, ω)-space. 7 Such grids of tabulated values, usually referred to as scattering kernels in the context of Monte Carlo applications like MCNP, can in principle provide a high degree of detail, but their construction and subsequent validation require significant efforts.
One appealing option is to directly measure S(Q, ω) in a neutron scattering experiment, which has obvious advantages but also disadvantages depending on experimental precision and coverage of (Q, ω) space, as well as access to an appropriate neutron instrument in the first place. The other option is theoretical calculations, which usually must employ some approximations due to the complexity involved in an ab initio approach. One such powerful if computationally expensive numerical approach is molecular dynamics simulations in which semiclassical modelling of atomic trajectories are used to provide expectation values for atomic displacements and correlations [19].
On the other end of the accuracy spectrum is the usage of various empirical closed form expressions [49,50,51] describing the total inelastic cross section as a function of neutron energy or wavelength. Not only do such formulas 7 While phonon terms contain δ-functions ensuring momentum and energy conservation, the convolution with phonon state densities and integration over grain orientations in the powder approximation eliminates them from the final scatter functions.
not provide details on scatter angles or energy transfers, they usually require additional tuned parameters and tend to be valid either at very long or very short neutron wavelengths. Although it is possible to patch two such formulas together [31,52] in order to cover both ends of the wavelength spectrum, the resulting behaviour at intermediate wavelengths tends in general to be highly inaccurate.
Returning to Eq. 21, it is possible to employ various approximations in order to extract results. This approach is in particular useful at longer neutron wavelengths, where the single-phonon (n = 1) term dominates, and indeed experimental techniques like neutron spectroscopy tend to focus on using the single-phonon term to access information about the dynamics of investigated samples, with contributions from multi-phonon terms (n ≥ 2) usually seen as unwanted background to be estimated and subtracted.
It is planned for NCrystal to support the modelling of inelastic scattering with enhanced data like scattering kernels or phonon density histograms when available, and prototype code with novel capabilities for precision sampling already exist [53]. It is nonetheless desirable that reasonably accurate modelling of inelastic scattering components should be provided by NCrystal even for materials where no more data than that required for Bragg diffraction is provided.
Thus, an iterative procedure suggested in [54] has been adopted for NCrystal, in which the contribution to the scatter function involving n phonons are estimated from the contribution involving n − 1 phonons, starting from singlephonon contributions determined by the Debye model (cf. Section 2.5). Additionally working under the so-called incoherent approximation [55], in which off-diagonal elements of j, j are assumed to cancel each other out in coherent inelastic scattering, allows the prescribed method to estimate both coherent and incoherent components of inelastic scattering.

The Debye model
It is clear from the discussion so far that any actual evaluation of scattering functions or cross sections involves evaluation of the Debye-Waller function defined in Eq. 20, which is not surprising since this function captures the dynamics of the system due to thermal fluctuations. In principle the evaluation requires highly non-trivial material-specific information, either based on theory, numerical work, measurements, or a combination of those. In order to be able to provide meaningful results for user defined materials which might lack such With α denoting a Cartesian coordinate and averaging over phonon polarisations, [57, in the notation of the present paper implies: Where the index i runs over all phonon states, N is the number of phonon states, M is the atomic mass and T is the material temperature. The sum over phonon states can be replaced by an integral with the state density, ρ( ω) = Now, as phonons in the Debye model propagate at constant velocity, the frequency of a phonon will be proportional to its momentum. Assuming that phonon states carry momenta distributed uniformly in momentum space, this implies ρ( ω) ∝ ω 2 . Fixing the normalisation Inserting into Eq. 37 and changing the integration variable from ω → u to u ≡ ω/k B T , this can be written as: Where: A result which is identical to [57,Eq. 11]. It is interesting to note that the displacements predicted by Eq. 38 do not vanish at 0 K, which is a quantum mechanical effect related to the wave functions in the ground states. For a given atomic mass, temperature and Debye temperature, it is straight-forward to evaluate f (T /T D ) numerically and thus approximate the Debye-Waller function as: where δ 2 j is the (isotropic) mean-squared displacement provided by Eq. 38. In principle the Debye temperature could be allowed to vary for each site j in the unit cell, but in practice it is usual to allow it to be specified for each type of element found in the crystal and NCrystal accordingly supports the specification of Debye temperatures either globally or per type of element. It is often the case that Debye temperatures for a given crystal can be found in the literature (e.g. [58,59]   of temperature for a number of crystals. Indeed, as required by the harmonic approximation, it is generally true that the predicted displacements are much smaller than typical inter-atomic distances, although the validity is strongest at lower temperatures.

Core framework overview
At its core, all capabilities of the NCrystal toolkit are implemented in an object oriented manner in a C++ library, providing both clearly defined interfaces for clients and internal separation between code implementing physics models, code loading data, and code providing infrastructure needed for integration into final applications. As will be discussed in further detail in Sections 5 and 6, it should be noted that while it is certainly possible to use the C++ classes discussed in the present section directly, typical users are expected to use other more suitable interfaces for their work -employing also a simpler and generic approach to material configuration.
In Figure 2 is shown the most important classes available in the release of NCrystal presented here, along with their internal inheritance relationships. At the root of the tree sits a few infrastructure classes not directly related to actual physics modelling, starting with the universal base class RCBase, which provides reference counted memory management for all derived classes. 8 One level below this sits the CalcBase class, from which all classes actually implementing physics modelling code must inherit. The main feature of CalcBase is currently that it provides derived classes with access to pseudo-random numbers, in a manner which can be configured either globally or separately for each CalcBase instance.
In order to do so, each CalcBase instance keeps a reference to an instance of a class deriving from RandomBase -a class which provides a generic interface for pseudo-random number generation. The primary purpose of this setup is to let NCrystal code use external sources of random numbers when embedded into existing frameworks managing their own random number generators, like those discussed in Sections 6.3 and 6.4. If no particular random generator is otherwise enabled, the system will fall back to using a xoroshiro128+ [60] generator implemented in the RandXRSR class.
Currently, the only class deriving directly from CalcBase is Process. This is an abstract class representing a physics process, specifying a general interface for calculation of cross sections as a function of incident neutron state. The returned cross section values are normalised to the number of atoms in the sample, and the neutron states are specified in terms of kinetic energy (E i ) and direction (k i ), which is equivalent to providing k i but using parameters typically available in 8 Although C++11 provides alternatives for such reference counting in the form of modern smart pointers, it is for the time being the aim of NCrystal to support also C++98, in which support for such is incomplete. Additional methods available via the Scatter interface:  processes are divided into those that are oriented in the sense that they actually depend onk i and those (isotropic) ones that do not. For the latter, one can access the cross section without specifying a direction. For reference, the precise methods available via the Process interface can be seen in Table 1. The Scatter class does on the other hand extend the Process interface, adding methods for random sampling of final states, as can also be seen in Table 1. Specifically, this sampling generates both energy transfers, E f − E i , and final direction of the neutron,k f . In the case of isotropic processes, only the energy transfer and the polar scattering angle, θ, fromk i tok f will depend on the actual physics implemented, and will itself be independent ofk i .
The azimuthal scattering angle will be independent and uniformly distributed Finally, the Info class is a data structure containing information about crystal structure at the microscopic scale of crystallites. As illustrated in Figure 3, the Info class serves the role of separating sources of crystal definitions from the physics models using the data as input. The data fields available on the Info class are provided in Table 2. Not all input sources will be able to provide all the data fields shown, nor will a particular physics modelling algorithm require all fields to be available. Thus, to ensure maximal flexibility for both providers and consumers of Info objects, all data fields are considered optional. If a given physics algorithm does not find the information it needs, it should simply indicate a configuration error, by throwing an appropriate C++ exception. When configuring materials via the recommended method for end-users (cf. Section 5), the factory code will generally try to avoid triggering undesired errors when the user intent seems obvious: if for instance the input data does not provide any information out of which it would be possible to provide cross sections for inelastic scattering processes, it is most likely that the user is simply only interested in elastic physics and no instantiation of inelastic processes will be attempted.
This scheme of decoupling data loading code from physics models allows NCrystal to easily support data read from a variety of input files, and makes it simple to add support for additional input sources in the future. If desired,

Contents Units
StructureInfo Basic info about unit cell: Lattice parameters (a, b, c, α,    In Listing 1 is shown an example of C++ code loading a crystal definition and constructing objects able to model absorption and scattering in a polycrystalline or powdered material. The example illustrates in practice how some of the classes discussed in this section are used, but also indicates the complexity arising out of the need to provide model-specific parameters. In Section 5, a more convenient method for initialisation and configuration will be introduced.

Physics models
At the heart of NCrystal is the actual modelling of specific physics processes, provided via the classes shown with blue outlines in Figure 2. The only mod-elling of absorption cross sections presently available is provided in the AbsOOV class, which implements the 1/v scaling model discussed in Section 2.2, by scaling the value of σ abs given at the reference velocity of 2200 m/s (cf. Table 2).
Despite being implemented with a simple model, absorption cross sections are nonetheless provided in NCrystal for completeness, in order to carry out validation against data from measurements of total cross sections like those presented in Section 4.4 or to facilitate the creation of plugins for Monte Carlo simulations like the one presented in Section 6.4. Additionally, it is important to reiterate that for most nuclear isotopes, the 1/v modelling is actually remarkably accurate at thermal neutron energies (cf. Section 2.2).
The remaining physics models all concern scattering processes. Three processes implement the coherent elastic physics of Bragg diffraction discussed in For reasons of scope, the three models of Bragg diffraction will be presented in detail in a future dedicated publication, covering details of both theory, implementation and validation.
Two models currently provide incoherent and inelastic physics: BkgdExtCurve and BkgdPhonDebye. 9 The former is a simple wrapper of externally provided σ bkgd (λ) cross section curves (cf. Table 2), and exists mostly for reference, allowing nxslib-provided cross section curves to be exposed via NCrystal when using .nxs files (cf. Section 4.2). Unless selected explicitly, the factories dis- 9 The term "Bkgd" is here used as a short-hand for inelastic and incoherent processes. It is merely used to lighten the notation, and of course is not meant to imply that e.g. inelastic neutron scattering is always to be considered "background" to some other signal (which would be incorrect). cussed in Section 5 will always prefer the more accurate model provided by BkgdPhonDebye. This improved model relies on the incoherent approximation and the iterative procedure discussed in Section 2.4, and will be presented in detail in a future dedicated publication. In the current release of NCrystal it does not provide any detailed sampling of scatter angles or energy transfers, and hence like BkgdExtCurve it derives from the ScatterXSCurve class. When asked to generate a scattering, ScatterXSCurve will scatter isotropically and either model energy transfers as absent (i.e. elastic), or as "fully thermalising".
In the latter case, the energy of the outgoing neutron will be sampled from a thermal (Maxwell) energy distribution specified solely by the temperature of the material simulated. Both of these choices use well defined distributions, but are clearly lacking in realism. It is the plan to improve this situation in the future, in order to provide a complete and consistent treatment of all components involved in thermal neutron scattering. For most materials, this will be achieved by adopting proper sampling models in BkgdPhonDebye, still based on the currently used iterative procedure for cross section estimation. As already discussed in Section 2.4, it is additionally planned to allow even higher accuracy and reliability when modelling inelastic scattering in select materials, by introducing data-driven models which are able to utilise pre-tabulated scatter-kernels or phonon state densities, if available.
databases [62,63], but ultimately abandoned due to the flexible nature of such files, many of which were found in practice to not contain suitable information for the purposes of NCrystal. Manually selected CIF files were, however, read with [64] and used to assemble a library of .ncmat files which is shipped with NCrystal. These files, providing the crystal structures listed in Table 3 Table 3 will be expanded in the future, in response to requests from the user community.

NCrystal material files (.ncmat)
Intended to be straight-forward to parse algorithmically, while at the same time being directly legible to scientists with crystallographic knowledge, the layout of NCrystal material files is deliberately kept both simple and intuitive, as illustrated with the sample file in Listing 2. These simple text files always begin with the keyword NCMAT and the version of the format. Next comes optionally a number of #-prefixed lines with free-form comments, intended primarily as a place to document the origin, purpose or suitability of the file. The data itself follows after this introduction, and is placed in four clearly denoted sections as shown in the listing. The @CELL and @ATOMPOSITIONS sections define the unit cell layout directly, and the definition must be compatible with the space group number which can optionally be provided in the @SPACEGROUP section. Finally,   The atomic positions are specified using relative lattice coordinates (cf. Eq. 1), and the name of an element is also required at each such position.
It is not necessary, nor currently possible, to specify element-specific data such as masses, cross sections or scattering lengths. Instead, NCrystal currently includes an internal database of such numbers for all natural elements with Z ≤ 92, based on [74,75,76]. Although somewhat inflexible, this scheme provides a high level of convenience, consistency and robustness by significantly lowering the amount of parameters required in each .ncmat file. It is envisioned that a future version of NCrystal would support more specialised use-cases by supporting not only the optional specification of element-and isotope-specific data in .ncmat files, but also to allow for specification of features like chemical disorder, impurities, doping or enrichment.
With the information parsed directly from the .ncmat file and the specification of a material temperature, the remaining information shown in Table 2 will be derived numerically at initialisation time (with the exception of σ bkgd (λ) haviour is much preferable to the O(1/d 3 cut ) behaviour of an implementation without such an early-abort.
In addition to the calculation of squared form-factors at all considered hkl points, points must also be sorted into families of points sharing d-spacing and squared form-factor values. This is done with a map-based O(log(d cut )) algorithm. All together, the final result is an hkl list initialisation algorithm which is fast enough that the trade-off between accuracy and initialisation time inherent in the choice of d cut is not very severe.
The value of d cut can always be set directly by the users, but in order to determine appropriate default values for the majority of users who are not expected to do so, the impact on Bragg diffraction cross sections in the powder approximation (cf. Eq. 35) was investigated. Although most powder diffraction experiments would concentrate on wavelengths longer than 2d cut , and therefore not be directly affected by it, the total cross section for diffraction in a powder is still a meaningful benchmark. After all, if a too high cut-off value is chosen, too many planes would be left out and the total cross section would be underestimated at shorter wavelengths -with corresponding degradation of realism in a simulation of for instance beam filters or shielding based on powders or polycrystalline materials. Thus, Figure 5 shows the relative impact of d cut on the cross section in the limiting case λ → 0. The impact levels gauged in this limit are clearly very conservative, since at this wavelength all omitted planes satisfy Eq. 32 resulting in a large relative impact, while the factor of λ 2 in Eq. or less, a number which is unlikely to cause concern for casual users, the default value is raised to d cut = 0.25 Å for materials with more than 40 atoms in the unit cell.
The resulting impact on the cross section curve of both thresholds applied by default, f cut = 10 −5 b and d cut = 0.1 Å or 0.25 Å, is shown in Figure 6. For most materials the effect is completely negligible, and even in the worst case of yttrium-oxide, the effect exists only at very short wavelengths and is hardly noticeable.  Relative impact on σ powder coh,el at λ

.nxs file loader
Routines in the nxs library [30,77,78], used in McStas [12,13], VitESS [14,15] and NXSG4 [31] to provide cross sections in crystal powders, comes with an associated text-based file format for crystal structure definition, recognised by the extension .nxs. The file format, described in [31], contains information similar to that in the .ncmat file of Listing 2 with a few notable differences.
The first one is that element-specific cross sections and masses must be specified file support would be exposed to these extra license requirements, the NCrystal code implementing the support for .nxs file loading is kept clearly separated from the rest, and it is straight-forward to build NCrystal without it. 11 Upon loading .nxs files, most information shown in Table 2  and subsequently ignoring all entries inevitably generated with a d-spacing beyond this value. As the nxs library code for grouping hkl entries into families is implemented with a relatively slow linear algorithm (resulting in an overall algorithmic complexity of O(N 4 max )), the d cut value selected by default for .nxs files is somewhat higher than for .ncmat files. It is chosen so as to correspond to N max = 20 but with d cut at most 0.5 Å and at least 0.1 Å. To prevent perceived programme lockups, it will result in an error if a user requests a d cut value which requires N max > 50.
Another important point is that the nxs code was not written to support single crystal modelling, and therefore the loaded HKLInfo objects will contain no lists of normals or hkl indices. However, when hkl family composition corresponds to symmetry equivalence groups (as is indeed the case for .nxs files), the NCrystal single crystal code is able to reconstruct such lists of normals on demand if absent, and it will therefore still be possible to use .nxs files for single crystal simulations in NCrystal. Finally, the nxs library provides estimates of inelastic/incoherent scattering cross sections based on various empirical formulas, and these will be provided in the σ bkgd (λ) field of Info objects when an .nxs file is loaded. By default, the cross section curve provided is similar to the one discussed in [31], representing an ad hoc combination of empirical formulas due to Freund [49] and Cassels [50]. As discussed in Section 3.1, these curves are provided for reference: the native NCrystal algorithms provide more robust predictions for inelastic/incoherent cross sections and is used by default also when working with .nxs files.

.laz and .lau file loader
Due to the widespread usage in various McStas components for modelling of Bragg diffraction, NCrystal also supports the loading of .laz file from [81] and .lau files from [82]. Both of these very similar text-based file formats can be generated from a CIF file by the cif2hkl application, part of iFit [83], and their most notable feature is that they directly contain hkl lists with d-spacings and form factors. Typically, .laz files are used in the context of powder diffraction, while .lau files can be used to deal with single crystal diffraction as well. This distinction implies that the latter files are larger, since they break down each hkl family into multiple lines of data, in order to provide enough information that all plane normals associated with a given family can be directly inferred.
The NCrystal code loading such files is rather simple, since no special calculations are needed to fill the HKLInfo section. By default, no thresholds are applied upon loading the hkl information from the file, but if desired it is of course possible to specify a custom d-spacing threshold in order to ignore some hkl families in the file. Other information in Table 2 loaded from information at the beginning of the files are StructureInfo, σ abs , and density. Finally, it is also possible to specify the temperature when loading the file, but it is important to note that since form factors are hard-coded in the file itself, their values are unaffected by the actual value provided. All in all, the loaded information is sufficient for modelling of Bragg diffraction, but absent is information which could be used to model inelastic or incoherent components. Thus, the factories described in Section 5 will create processes without such components for .laz or .lau files.

Validation
In order to simultaneously validate not only the code responsible for initialising crystal structure information from .ncmat and .nxs files, but also the individual data files provided with NCrystal, multiple approaches were pursued. One particular concern is the code responsible for creating hkl lists with  Table 3. Files marked with N are poly-atomic crystals with per-element Debye temperatures, which is not supported in .nxs files and which accordingly had to be compared using files in which a global Debye temperature was substituted.
Next, the hkl lists created from .ncmat files are validated again by using them to generate artificial powder diffraction spectrums and testing them with software which is normally used to decode crystal structures from such spectrums at real powder diffraction experiments. First, a simple neutron diffraction instrument was simulated with McStas (cf. Section 6.4), using NCrystal and the relevant .ncmat file to model a crystal powder sample. 12 As a result, a powder spectrum was produced for each of the four tested materials -marked with G in the last column of Table 3. These spectra were then used as input to GSAS-II [84], which in all cases managed to recover the crystal structure which was present in the original CIF file from which the relevant .ncmat file was produced. As an example, Figure 7 shows a refinement for a powder spectrum with 9.1 × 10 7 simulated neutrons collected in the detector array after interaction with a sapphire sample, and Table 4 shows the corresponding parameters extracted by GSAS-II, which in addition to space group and atomic (Wyckoff) positions also include atomic mean-squared displacements. The goodness-of-fit is provided by GSAS-II in terms of R-factor, which in all cases was less than 1.9%, indicating a very high degree of compatibility. Due to the usage of a full-scale simulation with actual NCrystal components enabled, these comparisons incidentally validate features of NCrystal beyond just the initialisation of crystal structure information -but dedicated future publications will document additional validations performed for these more thoroughly.
In addition to performing a full-blown instrument simulation with McStas, an idealised diffraction spectrum for a given neutron wavelength can be constructed directly from the d-spacings, multiplicities and squared form factors in a given 12 For reference, the instrument setup used was similar to the one shown in Listing 8, but with a few modifications carried out in order to increase the quality of the produced diffraction patterns. Thus, the sample size was reduced slightly, the number of detector bins increased, linear collimators were added before and after the monochromator, and a radial collimator was placed in front of the detector. Finally, the monochromator was changed to Germanium-511, in order to select wavelengths around 1.   are sufficiently realistic that it is possible to extract the crystal structure from them using FullProf [85]. For the 17 materials tested in this manner, marked with F in the last column of Table 3, FullProf was able to precisely recover their crystal structures -returning R-factors which were in all cases better than 2.7%. Figure 8 shows an example of a powder spectrum fitted with FullProf and Table 5 the corresponding extracted parameters.
As another independent validation of squared form factor predictions, mea-   2012009 in [62]), and directly compared against the ones predicted by NCrystal.
The result for d-spacings larger than 1.5 Å is shown in Figure 9, indicating a good agreement. Accordingly, this material is marked with R in the last column of Table 3.
As a final validation, experimental measurements of energy-dependent total cross sections were obtained from EXFOR [86] where available and tested against predictions from NCrystal. This tests not only .ncmat files and the associated crystal initialisation code, but also partly the implemented physics processes for powders. Most materials were validated in this manner, as indicated with T or T in the last column of Table 3  uranium oxide. Finally, the plot for tin in Figure 14 is interesting in that it confirms a good agreement of Bragg edges and single phonon scattering (despite some unclear result from one data set from an unpublished measurement which could be due to texture). But around 1.3 eV, effects of a nuclear resonance can be seen, which is not currently included in the modelling provided by NCrystal.

Factories and unified configuration
The object oriented implementation of NCrystal in a C++ class hierarchy described in Section 3, provide a high degree of flexibility. This flexibility is an advantage for experts wishing to extend or customise NCrystal, benefiting from low-level access to the various classes and utilities exposed to developers as part of the NCrystal API. However, as demonstrated by the example in Listing 1, this flexibility comes at a cost of complexity which it is desirable to contain. Not only in order to make the usage of NCrystal simpler and less error prone for     Due to this, the configuration of a single crystal material from that file will result in the anisotropic mosaicity distribution appropriate for pyrolytic graphite to be used rather than the standard isotropic one. Table 6 documents the configuration parameters most likely to be of interest to typical NCrystal users, but for clarity a few example strings and their implications will be presented in the following. 13 To make it easier to embed such configuration strings into other text-based data sources, it is mandated that they are only valid if they consist of simple ASCII characters (excluding control-codes) and without any of the following characters: Alternatively, the crystal direction can be specified in reciprocal lattice space rather than direct space as "@crys_hkl:ch,ck,cl@lab:lx,ly,lz". dir2 Secondary orientation of crystal, specified using the same syntax as for the dir1 parameter. Only components of dir2 orthogonal to dir1 will actually influence the orientation (but see dirtol). lcaxis Symmetry axis of anisotropic layered crystals similar to pyrolytic graphite (as vector, e.g. "0,0,1"). • "Al_sg225.ncmat", "Al_sg225.nxs", or "Al.laz" Polycrystalline or powdered aluminium at room temperature, using appropriate code for loading .ncmat, .nxs or .laz files respectively. Although in principle representing the same material, the realism will depend on the type of input file, as NCrystal automatically selects the most realistic modelling possible while also considering issues like load times. In this particular case, the scattering models based on the .laz file will not contain any incoherent or inelastic components, and the Bragg diffraction based on the .nxs file will contain fewer reflection planes at shorter d-spacings than the one based on the .ncmat file (cf. Section 4.2).
• "Be_sg194.ncmat;bkgd=0;temp=100K" The same beryllium again, but this time the only scattering being modelled is that of Bragg diffraction. While clearly decreasing the modelled realism, avoiding incoherent or inelastic scatterings is occasionally useful for reasons of clarity or validation.
in the filesystem.
Furthermore, the modelling was made faster at the cost of realism by increasing the d-spacing cutoff and disabling incoherent and inelastic components with bkgd=0.
Note that as discussed above, the lcaxis=0,0,1 parameter is embedded into the file itself, enabling the proper anisotropic mosaicity distribution.
This approximation could be useful in order to efficiently model a single crystal sapphire filter in a beam, assumed to be positioned such that the Bragg condition is unsatisfied for all planes.

Interfaces and bindings
The method of material configuration presented in Section 5 facilitates the usage of NCrystal in a variety of contexts, consistently employing the same configuration strings and data files everywhere. This allows the sharing of material configurations between applications and users, and implies a freedom of choice when tuning or validating such configurations -irrespective of which Monte Carlo application is ultimately used to study a given problem. Accordingly, Section 6.1 will present a tool for quick inspection and tuning of material configurations, while Section 6.2 will present language bindings available to advanced users needing direct access to NCrystal functionality -including  First of all, when providing just one configuration string, the corresponding material will be created as a powder, and two plots will be produced: the components of the resulting neutron cross sections and a distribution of randomly sampled scatter angles. For instance, the following command can be used to inspect a sapphire powder with default settings for temperature, packing-factor, d-spacing cut-off, etc.: The resulting plots are shown in Figure 15. Figure 15.a shows the components of the total interaction cross section, including absorption since the -a flag was specified. Figure 15.b shows a two-dimensional scatter-plot: at each neutron wavelength, scatterings are sampled with statistics proportional to the scattering cross section at that wavelength and the resulting scatter angles shown. By supplying the --dump flag, the graphical plotting will be replaced with a printout of loaded crystal information. For instance, the command: n c r y s t a l _ i n s p e c t f i l e --dump " C u 2 O _ s g 2 2 4 _ C u p r i t e . ncmat ; dcutoff =1 Aa " results in the printout shown in Listing 3. Figure 16 shows how specifying more than one configuration string at a time leads results in a plot comparing the resulting cross sections. Figure 16.a shows the total interaction cross section of thermal neutrons in a range of polycrystalline metals -which might for instance be considered for a support structure in some parts of a neutron instrumentand is the result of the command (all on one line):  n c r y s t a l _ i n s p e c t f i l e -a " Be_sg194 . ncmat ; temp =100 K " " Be_sg194 . ncmat ; temp =200 K " " Be_sg194 . ncmat ; temp =300 K "

Command-line inspection
Finally, the command ncrystal_inspectfile --test can be used to validate a given installation of NCrystal. Although ncrystal_inspectfile already provides a useful set of functionality, it is expected that additional standard tools will be added in the future. In particular it would be useful with more options for plot creation and with utilities for assisting with the creation of configuration strings involving the alignment of single crystals, for instance to assist in the configuration of simulations involving neutron monochromators.

C++, C and Python bindings
Using the configuration strings presented in Section 5, the somewhat complicated C++ code in Listing 1 can be significantly simplified as shown in Listing 4: not only is it fewer lines, but the only material-specific code is the one defining the configuration string variable named cfg. It is thus straight-forward to create C++ applications or plugins which use NCrystal as a backend, but lets users provide the configuration strings in whichever frontend is relevant to the task at hand. The lingua franca of software is, however, the C programming language, Listing 3: Printout produced by invoking ncrystal_inspectfile --dump on a single configuration string (here "Cu2O_sg224_Cuprite.ncmat;dcutoff=1Aa"). Refer to Table 2 for an explanation of the listed parameters.   Finally, it is possible to use NCrystal directly from Python [88], which mostly exists as a feature to support the usage of NCrystal for advanced scripting, analysis and plotting work. Python and C++ share concepts like classes and exceptions, and Python code using NCrystal thus ends up being similar and possibly even simpler than the corresponding C++ code, as can be seen in the Python equivalent of Listing 4, shown in Listing 6. Additionally, for efficiency and convenience, the Python interface supports vectorised access through Numpy [89] arrays. For a simple example using this, Listing 7 illustrates how to create a Matplotlib [90] plot of the scattering cross section as a function of neutron wavelengths between 0 Å and 10 Å.

Geant4 interface
Assuming Geant4 is configured to use a physics list which include the socalled HP (high precision) models for neutron physics, almost all comparisons between cross sections in Geant4 and NCrystal are qualitatively equivalent to the ones shown in Figure 17 for a magnesium powder. Starting around the keV scale, Geant4 provides detailed modelling of higher energy effects such as those related to nuclear resonances. At lower energies, the cross section curves are, however, completely smooth and featureless due to the free-gas approximation used. NCrystal, on the other hand, provides detailed structure-dependent scattering physics at the sub-eV scale -but has no capacity for modelling physics at the keV scale. At intermediate energy scales, neither nuclear resonance physics or material structure-dependent physics introduce significant features, and there is an overlap in predictions between NCrystal and Geant4. For absorption processes, material structure is unimportant, and the predictions from NCrystal and Geant4 are in perfect agreement over the range covered by NCrystal, due to the validity of the simple 1/v scaling of such cross sections (cf. Section 2.2).
However, unlike NCrystal, Geant4 provides detailed modelling of the secondary particles produced in absorption reactions. Consequently, the NCrystal plugin for Geant4 does not touch the absorption physics at all, instead focusing on replacing just the scattering physics for low energy neutrons. Presently this is done with a global hard-coded cross over point of 5 eV, but the exact transition model might be revisited in the future, since a few rare isotopes have nuclear resonances lower than this. As an example, the experimental data in Figure 14 indicates a resonance around 1.3 eV.
Geant4 user code requires very few changes in order to enable NCrystal modelling of scattering for low energy neutrons. First of all, the G4NCrystal.hh header file must be included. Next, materials for which it is desired to use NCrystal to provide scattering physics for thermal neutrons should be identified.
Instances of G4Material for these materials must then be created by providing appropriate NCrystal configuration strings: sures that those neutrons in the incoming white beam possessing a compatible wavelength, will be reflected at exactly 135°in the xy-plane, which is the defining feature of a single crystal neutron monochromator. Multiple scattering and geometrical boundaries are naturally accounted for, and the reflected neutrons exhibit a characteristic "zig-zag" walk, 14  Another example is shown in Figure 19, where it is illustrated how scattering of a monochromatic pencil beam of neutrons in an (untextured) polycrystalline sample changes qualitatively when NCrystal is enabled: instead of diffuse scattering due to the free-gas approximation, proper scattering into Debye-Scherrer cones by crystal planes is observed.

McStas interface
In McStas simulations, thermal neutrons are passed through an ordered list of components configured by users in a so-called instrument file. Components typically represent actual in-beam elements found at the modelled neutron instrument such as: source, optical guides, choppers, filters, monochromators, analysers, samples, or detectors. Each component is responsible for modelling both geometrical and physics effects, and in addition to scattering, absorption physics can be implemented either by trajectory termination or intensity reduction.
The NCrystal plugin for McStas is provided as a component with the name NCrystal_sample, but it can be used to model a variety of elements in addition to samples, including filters, monochromators, and analysers. For now, it have been extracted programmatically but for simplicity it was in this case determined by the user via the interactive tool described in Section 6.1. Figure 20 shows a 3D visualisation of the resulting simulation using McStas 2.4.1, while Figure 21 shows the corresponding diffraction pattern observed in the modelled detector array.
In addition to benefiting from future improvements to the NCrystal library, it is foreseen that the NCrystal_sample McStas component itself might also be further enhanced at a later state. In particular, it would be desirable to implement variance reduction techniques, at least in terms of making it possible to focus outgoing neutrons towards the next down-stream component. It is also likely that use-cases for more advanced geometrical layouts will arise, and the code has consequently been structured in a way which makes it straight-forward to add such features.

Outlook
The presented toolkit for thermal neutron transport is arguably unique in its attention to interfaces and capability for integration into various technical contexts and is already in version 1.0.0 very capable in terms of modelling of interactions in single crystals and crystal powders, and has already been used to enable a range of interesting studies (e.g. [20,53,91,92,93,94,95,96,97]).
Nonetheless, work has already begun on several improvements to both physics models and the framework itself. Firstly, as mentioned in Section 2.4, the capabilities for modelling of inelastic and incoherent scattering should see significant enhancements -with the possibility of supporting liquids or polymers to some extent. These models, and those implementing Bragg diffraction, will be described in detail in future dedicated publications.
Next, it is the plan to carry out framework extensions and refactorisations which will allow NCrystal to support enhanced material realism, by making the exact composition of materials and crystals more customisable. Once implemented, it should on one hand become possible to support multi-phase materials -needed for realistic multi-phase metal alloys or crystal powders suspended in liquids -and on the other hand the composition of each phase should become more flexible as well, allowing for enriched materials or chemical disorder. That would enable modelling of crystals in which some sites are not fully occupied in all cells -or occasionally occupied with elements playing the role of contaminants or dopants.
Several more technical developments are envisioned as well: planned interface extensions will enable better support of multi-threaded applications (such as ANTS2 or multi-threaded builds of Geant4), and several use-cases have been identified where it would be advantageous to be able to initialise crystal data directly from process memory rather than needing on-disk files. In the longer run, once all relevant platforms and applications support it, it is also planned to drop the support for the C++98 standard, in order to better benefit from modern C++ features and better cross platform support introduced in C++11 and beyond.
Beyond that, the future directions will depend on resources and community interest. At the very least the library of data files and the list of Monte Carlo applications with NCrystal support are both expected to expand. But given sufficient interest and contributions new ambitious physics models could be added, ranging from treatment of texture, bent crystals and new anisotropic mosaicity models to better facilities for dealing with nuclear resonances or branching into new areas like magnetic spin-dependent interactions or support of X-ray physics. Input, feedback, ideas or contributions are gratefully received via the NCrystal website [34].